I've been trying to solve the following PDE using the NDSolve function but it seems something is not working properly. The PDE is a the heat equation on polar coordinates and assuming angular simmetry:
$u_t (t, r) = \alpha(r) \frac{1}{r} \partial_r (r u_r (t, r))$
with boundary conditions $u_t(t, r_{in}) = p$ and $u(t, r_{max}) = 0$, and initial condition $u(0,r)=0$ for $r\in(r_{in}, r_{max})$. The function $\alpha$ represents the diffusion coefficient of the temperature in the media.
The idea is that $r_{max}$ is large enough so that $u(t,r)$ is alsmost flat and equal to zero close to $r_{max}$ through all the time window integrated.
I have tried different $\alpha(r)$ functions, which are
alpStep[r_] :=If[r < rout, aA, aB];
alpLin[r_] :=If[r < rout, aA + (aB - aA)(r-rin)/(rout - rin), aB];
alpA[r_] := aA;
alpB[r_] := aB;
where $r_{out}\in(r_{ín},r_{max})$ is a critical radius beyond which the diffusion remains constant.
The unexpected result is the following. I set aA > aB, and I integrate the different PDEs (with the different alp functions) by means of NDSolve. It turns out that the solution associated to the alpStep function takes values much more higher than the other functions, and I would expect it to be somewhere in between the solutions associated to the functions alpA and alpB. May be the problem is due to the discontinuous diffusion, but I can't see why.
The code is the following:
rin = 0.05;
rout = 0.15;
rmax = 5;
p = 0.01;
tend = 24*10;
aA = 0.01;
aB = 0.001;
alpStep[r_] := If[r < rout, aA, aB];
alpLin[r_] := If[r < rout, aA + (aB - aA) (r - rin)/(rout - rin), aB];
alpA[r_] := aA;
alpB[r_] := aB;
opts = Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}};
With[{u = u[t, r]}, eqn = alpStep[r] ((1/r) D[r D[u, r], r]) - D[u, t];
robinbc = NeumannValue[p*alpStep[rin], r == rin];
bc = DirichletCondition[u == 0, r == rmax];
ic = u == 0 /. {t -> 0};]
solStepAB =
NDSolveValue[{eqn == robinbc, bc, ic},
u, {t, 0, tend}, {r, rin, rmax}, opts];
With[{u = u[t, r]}, eqn = alpLin[r] ((1/r) D[r D[u, r], r]) - D[u, t];
robinbc = NeumannValue[p*alpLin[rin], r == rin];
bc = DirichletCondition[u == 0, r == rmax];
ic = u == 0 /. {t -> 0};]
solLinAB =
NDSolveValue[{eqn == robinbc, bc, ic},
u, {t, 0, tend}, {r, rin, rmax}, opts];
With[{u = u[t, r]}, eqn = alpA[r] ((1/r) D[r D[u, r], r]) - D[u, t];
robinbc = NeumannValue[p*alpA[rin], r == rin];
bc = DirichletCondition[u == 0, r == rmax];
ic = u == 0 /. {t -> 0};]
solConA =
NDSolveValue[{eqn == robinbc, bc, ic},
u, {t, 0, tend}, {r, rin, rmax}, opts];
With[{u = u[t, r]}, eqn = alpB[r] ((1/r) D[r D[u, r], r]) - D[u, t];
robinbc = NeumannValue[p*alpB[rin], r == rin];
bc = DirichletCondition[u == 0, r == rmax];
ic = u == 0 /. {t -> 0};]
solConB =
NDSolveValue[{eqn == robinbc, bc, ic},
u, {t, 0, tend}, {r, rin, rmax}, opts];
Grid[{{
Plot[{solStepAB[t, rin], solLinAB[t, rin], solConA[t, rin],
solConB[t, rin]}, {t, tend/1000, tend}],
Plot[{solStepAB[t, rout], solLinAB[t, rout], solConA[t, rout],
solConB[t, rout]}, {t, tend/1000, tend}],
Plot[{solStepAB[t, rout + (rout - rin)],
solLinAB[t, rout + (rout - rin)],
solConA[t, rout + (rout - rin)],
solConB[t, rout + (rout - rin)]}, {t, tend/1000, tend}]
}}]
The picture with the temperature profiles for the different solutions alpStep (blue), alpLin (Orange), alpA (Green) and alpB (Red) at $r = r_{in}, r_{out}, 2 r_{out} - r_{in}$ from left to right is the following, where it can be seen that the profile associated to alpStep is much higher that the others:
FEMDocumentation/tutorial/FiniteElementBestPractice#588198981
in the help system or online here $\endgroup$