One term in non-linear equations prevents NDSolve to work

Question

I am currently working on the stratification of the core of the planet Mercury, meaning the formation of a conductive layer at the top of the convective core, with a moving interface between both layers. After some variable changes to simplify my equations, here is the system of equations I want to solve:

$$\frac{\partial T}{\partial \tau}(\tau,x) = \frac{1}{(s(\tau)-1)^2}\frac{\partial^2 T}{\partial x^2}(\tau,x) +\left(\frac{2}{1+x(s(\tau)-1)}\frac{1}{s(\tau)-1}+\frac{x}{s(\tau)-1}\frac{\mathrm{d}s}{\mathrm{d}\tau}(\tau)\right) \frac{\partial T}{\partial x}(\tau,x)$$

$$\left(s(\tau) \left(\frac{T_{c0}}{T_{s0}-T_{c0}} + T_s(\tau)\right) s'(\tau) + \frac{1}{2 y}T_s'(\tau)\right)\left(2s(\tau) - \mathrm{e}^{y s^2(t)} \sqrt{\frac{\pi}{y}} \mathrm{erf}(\sqrt{y}s(\tau))\right) = 4 y s^3(\tau) \left(\frac{T_{c0}}{T_{s0}-T_{c0}}+T_s(\tau)\right)$$

$$\frac{1}{1-s(\tau)}\frac{\partial T}{\partial x}(\tau,0) = \frac{r_c}{k(T_{s0}-T_{c0})}q_c\left(\frac{\rho C_p r_c^2}{k}\tau+t_0\right)$$

$$T_s(\tau) = T(\tau,1)$$ $$\frac{1}{s(\tau)-1}\frac{\partial T}{\partial x}(\tau,1) = -2 y s(\tau)\left(\frac{T_{c0}}{T_{s0}-T_{c0}}+T_s(\tau)\right)$$

with $T$ the temperature profile in the conductive layer, $s$ the interface position, $T_s$ the interface temperature, $T_{c0}$, $T_{s0}$, $r_c$, $k$, $\rho$, $C_p$, $t_0$ and $y = \frac{g_c \alpha r_c}{2 C_p}$ constants.

I have discretised the spatial part of these equations in order to get a system of ODE's using the functions ptoo and ptoode (see here). Then I have used the function Solve in order to rewrite equations in the form '$\frac{\mathrm{d}...}{\mathrm{d}\tau}=...$' for all variables. And finally solve the equations using NDSolveValue. But I have got an error from NDSolve (see below).

If I delete the term $\frac{\mathrm{d}s}{\mathrm{d}\tau}$ in the right hand side of the first equation, everything goes fine and NDSolve solves my equations without complaining.

Is there something I can do in order to make NDSolve solving the system with the problematic term? I have tried to rearrange the equations in order to give simplified equations to NDSolve, changing the method (StiffnessSwitching, FixedStep, StartingStepSize or increasing the maximum number of steps) and I always have errors like 'max number of points reached' or 'stiff system'.

Here is my code:

(*constants*)
rc = 2050 10^3;
cp = 850;
rho = 7200;
alpha = 5 10^-5;
gc = 4;
k = 40;
y = (gc alpha rc)/(2 cp);

(*parameters*)
s0 = 2049 10^3;
Tc0 = 2100;
Ts0 = Exp[(-alpha gc)/(2 cp rc) (s0^2 - rc^2)] Tc0;
t0 = 0.099 10^9 365.25 24 3600;


qcmb[t_] = 
 With[{a = 0.004891658583550395, b = 0.34057028569554804, 
   c = 1.0021984665846737`*^-15}, a + b E^(- c t)]

(*equations*)
energyAdim = ( 
     s[τ] (Tc0/(Ts0 - Tc0) + Ts[τ]) s'[τ] + 
      1/(2 y) Ts'[τ]) (2 s[τ] - 
      E^(y s[τ]^2) Sqrt[π/y]
        Erf[Sqrt[y] s[τ]]) == 4 y s[τ]^3 (Tc0/(Ts0 - Tc0) + Ts[τ]);
tempContinuityAdim = Ts[τ] == T[τ, 1];
heAdim = D[T[τ, x], τ] == 
   1/(s[τ] - 1)^2 D[
      T[τ, x], {x, 
       2}] + (2/(1 + x (s[τ] - 1)) 1/(s[τ] - 1) + 
       x/(s[τ] - 1) D[s[τ], τ]) D[T[τ, x], x];
bc1Adim = 
  1/(1 - s[τ]) D[T[τ, x], x] == rc/(k (Ts0 - Tc0)) qcmb[(rho cp rc^2)/k τ + t0] /. x -> 0;
bc2Adim = 
  1/(s[τ] - 1) D[T[τ, x], x] == 
    -2 y s[τ] (Tc0/(Ts0 - Tc0) + Ts[τ]) /. x -> 1;

(*initial conditions*)
iniTAdim = 
  T[0, x] == (Exp[-y (2 x (s0/rc - 1) + x^2 (s0/rc - 1)^2)] - 1)/(
   Ts0/Tc0 - 1);
iniTsAdim = Ts[0] == 1;
inisAdim = s[0] == s0/rc;

(*parameters for transforming PDE's in ODE's*)
nbrPoints = 100;
scalingFactor = 1000;
xDiffOrder = 2;

{xL, xR} = domain = {0, 1};
grid = Array[# &, nbrPoints, domain];

ptoo = pdetoode[T[τ, x], τ, grid, xDiffOrder];
toode[expr_Equal] := 
  With[{sf = scalingFactor}, sf # + D[#, τ] & /@ expr];
toode[expr_List] := toode /@ expr;

energyODE = energyAdim;
tempContinuityODE = ptoo[toode[tempContinuityAdim]];
heODE = ptoo[heAdim]; del = #[[2 ;; -2]] &; heODE = del[heODE];
{bc1ODE, bc2ODE} = {ptoo[toode[bc1Adim]], ptoo[toode[bc2Adim]]};

iniTODE = ptoo[toode[iniTAdim]];
iniTsODE = iniTsAdim;
inisODE = inisAdim;

(*rewritting the equations like 'd.../dτ = ...'*)
dTdtau[τ] = 
  Flatten@{ptoo@Derivative[1, 0][T][τ, x], Ts'[τ], 
    s'[τ]};
solveDerivative = 
  Solve[Flatten@{Collect[heODE, ptoo@T[τ, x]], 
     Collect[bc1ODE, 
      Flatten[{ptoo@T[τ, x], 
        ptoo@Derivative[1, 0][T][τ, x]}]], 
     Collect[bc2ODE // Simplify, 
      Flatten[{ptoo@T[τ, x], 
        ptoo@Derivative[1, 0][T][τ, x]}]], tempContinuityODE, 
     energyAdim}, dTdtau[τ]];

Solving the equations using the default method gives:

result = 
  NDSolveValue[{Thread[
     dTdtau[τ] == (dTdtau[τ] /. solveDerivative[[1]])], 
    iniTODE // Simplify, iniTsODE, inisODE}, {T /@ grid, Ts, s} // 
    Flatten, {τ, 0, 
    k/(rho cp rc^2) (10^9 365.25 24 3600 4.5 - t0)}];
(*NDSolveValue::mxst : Maximum number of 10000 steps reached at the point τ == 3.640731908397398`*^-12.*)

Using the StiffnessSwitching method, I am going a bit further and the error is different but we are still far from the end value $\tau_{end} = 0.22086$:

result = 
  NDSolveValue[{Thread[
     dTdtau[τ] == (dTdtau[τ] /. solveDerivative[[1]])], 
    iniTODE // Simplify, iniTsODE, inisODE}, {T /@ grid, Ts, s} // 
    Flatten, {τ, 0, 
    k/(rho cp rc^2) (10^9 365.25 24 3600 4.5 - t0)}, 
   Method -> "StiffnessSwitching"];
(*NDSolveValue::ndsz: At τ == 6.36963610146291`*^-11, step size is effectively zero; singularity or stiff system suspected.*)

Changing the number of points in the space grid (nbrPoints) or the scaling factor (scalingFactor) does not help: results are not converging:

And for a number of points larger than 250, I got the message error:

Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{\"EquationSimplification\"->\"Residual\"}.

---

Edit: the model

• heat equation in the conductive layer: $\rho C_p \frac{\partial T}{\partial t}(t,r) = \frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2 k \frac{\partial T}{\partial r}(t,r)\right)$

• energy budget at the interface $s(t)$: cooling of the convective layer equals the heat conducted through the interface: $-\rho C_p \int_0^{s(t)} \frac{\partial T_a}{\partial t}(t,r) \mathrm{d}V = -4\pi s^2(t) k \frac{\partial T_a}{\partial r}(t,s(t))$

• given heat flux at the top of the core $r_c$: $-k \frac{\partial T}{\partial r}(t,r_c) = q_c(t)$

• temperature and heat flux continuity at the interface: $T_s(t) = T(t,s(t))$ and $-k \frac{\partial T}{\partial r}(t,s(t)) = -k \frac{\partial T_a}{\partial r}(t,s(t))$

with $\rho$, $C_p$, $\alpha_c$, $g_c$ and $k$ constants (density, specific heat, thermal expansivity, gravity and thermal conductivity respectively), $T(t,r)$ the temperature profile in the conductive layer, $T_a(t,r)$ the temperature profile in the convective layer, $T_s(t) = T_a(t,s(t))$ the temperature at the interface $s(t)$ and $\frac{\partial T_a}{\partial r}(t,r) = -\frac{\alpha_c g_c}{r_c C_p}r T_a(t,r)$ the adiabatic gradient.

In order to simplify the equations, I have considered the following variable changes: $$\tau = \frac{k}{\rho C_p r_c^2}(t-t_0)$$

$$x = \frac{r}{s(t)} \mathrm{ if} r \leq s(t) \mathrm{and} \frac{r-r_c}{s(t)-r_c} \mathrm{if} r\geq s(t)$$

$$T(\tau,x) = \frac{T(t,r) - T_{c0}}{T_{s0}-T_{c0}}$$

$$s(\tau) = \frac{s(t)}{r_c}$$ with $T_{c0}$, $T_{s0}$ and $t_0$ constants.

---

Edit 2: version 8

Using version 8 of mathematica, I can solve the equations but I do not have spatial convergence (example with scaling factor = 1, nx = nbrPoints):

Answer 1

This system of equations can be solved by explicit Euler in time.

(*constants*)rc = 2050 10^3;
cp = 850;
rho = 7200;
alpha = 5 10^-5;
gc = 4;
k = 40;
y = (gc alpha rc)/(2 cp);

(*parameters*)
s0 = 2049 10^3;
Tc0 = 2100;
Ts0 = Exp[(-alpha gc)/(2 cp rc) (s0^2 - rc^2)] Tc0;
t0 = 0.099 10^9 365.25 24 3600;
tm = k/(rho cp rc^2) (10^9 365.25 24 3600 4.5 - t0);

qcmb[t_] := 
  With[{a = 0.004891658583550395, b = 0.34057028569554804, 
    c = 1.0021984665846737`*^-15}, a + b E^(-c t)];


T[0][x_] := (Exp[-y (2 x (s0/rc - 1) + x^2 (s0/rc - 1)^2)] - 
     1)/(Ts0/Tc0 - 1);
T[-1][x_] := (Exp[-y (2 x (s0/rc - 1) + x^2 (s0/rc - 1)^2)] - 
     1)/(Ts0/Tc0 - 1);
s[0] = s0/rc;
n = 200; tn = tm/4000; 
Do[s[i] = 
   s[i - 1] + 
    tn*(4 y s[
           i - 1]^3 (Tc0/(Ts0 - Tc0) + T[i - 1][1])/(2 s[i - 1] - 
            E^(y s[i - 1]^2) Sqrt[\[Pi]/y] Erf[Sqrt[y] s[i - 1]]) - 
        1/(2 y) (T[i - 1][1] - T[i - 2][1])/tn)/(s[
         i - 1]*(Tc0/(Ts0 - Tc0) + T[i - 1][1])); np = i; 
  If[s[i] <= 0, Break[]]; 
  T[i] = NDSolveValue[{(Ti[x] - T[i - 1][x])/tn == 
      1/(s[i - 1] - 1)^2 *
        Ti''[x] + (2/(1 + x (s[i - 1] - 1)) 1/(s[i - 1] - 1) + 
          x/(s[i - 1] - 
              1) ((4 y s[
                  i - 1]^3 (Tc0/(Ts0 - Tc0) + 
                   T[i - 1][1])/(2 s[i - 1] - 
                   E^(y s[i - 1]^2) Sqrt[\[Pi]/y] Erf[
                    Sqrt[y] s[i - 1]]) - 
               1/(2 y) (T[i - 1][1] - T[i - 2][1])/tn)/(s[
                i - 1]*(Tc0/(Ts0 - Tc0) + T[i - 1][1])))) Ti'[x], 
     1/(1 - s[i - 1]) Ti'[0] == 
      rc/(k (Ts0 - Tc0)) qcmb[(rho cp rc^2)/k tn*i + t0], 
     1/(s[i - 1] - 1) Ti'[1] == -2 y s[
        i - 1] (Tc0/(Ts0 - Tc0) + T[i - 1][1])}, Ti, {x, 0, 1}];, {i, 
   1, n}] // Quiet

T3 = Table[{tn*i, x, T[i][x]}, {i, 0, np - 1}, {x, 0, 1, .02}];

T2 = Interpolation[Flatten[T3, 1]];

{Plot3D[T2[t, x], {t, 0, tn*(np - 1)}, {x, 0, 1}, Mesh -> None, 
  ColorFunction -> Hue, AxesLabel -> {"t", "x", ""}, PlotLabel -> "T",
   PlotRange -> All], 
 ListLinePlot[Table[{tn*i, s[i]}, {i, 0, np - 1}], PlotRange -> All, 
  AxesLabel -> {"t", "s"}]}