NDSolve gives unexpected results when using the method of lines

Question

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:

Answer 1

Yes! The problem was that some information was missing in the problem I stated above. Specifically, due to the discontinuity of the function $\alpha$ at $r_{out}$, the solution fails to be differentiable at that point, and multiple solutions exist. NDSolve gives one of these following some criterium related to fluxes (the heat flux at one side of the discontinuity has to be equal to the heat flux at the other side). Notice that the flux depends on the thermal conductivity $k$, which is not explicitly stated in the problem above. In general the thermal diffusivity $\alpha$ is equal to the conductivity $k$ over the heat capacity $c$ (i.e. $\alpha = k / c$). In order to use properly the NDSolve function in this case (in the sense that the continuity of the flux is preserved at $r_{out}$) one should avoid working with diffusion coefficient and use $k$ and $c$ instead, that is one should write

With[{u = u[t, r]}, eqn = k[r] ((1/r) D[r D[u, r], r]) - c[r] D[u, t];
 robinbc = NeumannValue[p*k[rin], r == rin];
 bc = DirichletCondition[u == 0, r == rmax];
 ic = u == 0 /. {t -> 0};]