3
$\begingroup$

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:

sf=500-default method sf=500-StiffnessSwitching method

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):

interfacePosition

$\endgroup$
  • $\begingroup$ Are you sure the underlying model itself is correct? Since you've mentioned "moving interface", I guess this problem is somewhat related to this one, right? Then I guess $s(\tau)$ is probably between $0$ and $1$? Nevertheless, with nbrPoints = 25; scalingFactor = 1; and MaxSteps -> 10 10^5 option I find $s(\tau)$ hits zero at about t = 0.005882695744315565. BTW the Simplify in NDSolveValue can be taken away, it only slows down the code. $\endgroup$ – xzczd Feb 27 at 11:53
  • $\begingroup$ @xzczd I'm quite sure of the model. And yes the problem is related to the one you mentioned. That's where I have found the useful tools ptoo and ptoode. The interface moves between 0 and 1. But I just tried with your modifications and I still have an error: "NDSolveValue::ndsz: At [Tau] == 1.1510886780232115`*^-9, step size is effectively zero; singularity or stiff system suspected." Have you changed something else? $\endgroup$ – Mariel Feb 27 at 13:31
  • $\begingroup$ No, I don't make other modification. This seems to be another backslide of v11. I can reproduce the issue in v11.2, but my previous test is done in v9.0.1. So far I haven't found a way to adjust v11.2 to produce the result of v9.0.1. $\endgroup$ – xzczd Feb 27 at 14:28
  • 1
    $\begingroup$ The heat equation, and the other equations, are in spherical coordinates, explaining why it is $r^2$. I've done the variable change manually. I agree it's not the safest way but I did it so many times that I should have eliminate all the errors now. $\endgroup$ – Mariel Mar 1 at 11:43
  • $\begingroup$ How do you transform the integro-differential equation to pure PDE? $\endgroup$ – xzczd Mar 1 at 13:35
0
$\begingroup$

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"}]}

fig1

$\endgroup$
  • $\begingroup$ I'm sorry, but the solution actually hasn't converged. If we choose tn = tm/8000, we'll see the calculation stops at about τ=0.0012, and with tn = tm/16000, the calculation stops at about τ=0.0008, with tn = tm/32000, stops at about τ=0.0005 $\endgroup$ – xzczd Feb 28 at 5:45
  • $\begingroup$ This is only the first step of the algorithm. In the scheme we need to add a variable step and increase the accuracy of the calculation Ti[x]. $\endgroup$ – Alex Trounev Feb 28 at 11:59
  • $\begingroup$ Well, the ODE solver of NDSolve is already adaptive (and long-tested), but as shown in OP's question, it just fails to go to τ=0.22. I really doubt if a self-made adaptive ODE solver will behave better. Though OP claims the underlying model is correct, I still suspect something is wrong with the equation system itself. $\endgroup$ – xzczd Feb 28 at 12:18
  • $\begingroup$ I agree that in the Stefan problem, the equation relates $ds/dt$ and heat flow from both sides proportionally$\lambda \nabla T$. In this case, we see the connection $ds/dt$ and $\partial T/\partial t$ at the mobile border. $\endgroup$ – Alex Trounev Feb 28 at 13:27
  • $\begingroup$ Thanks for your answer. I think there's an error in the discretisation of ds/dt: it should be 1/(2 y (1 - 2 s[i - 1])) and not 1/(2y). With this correction the result is less beautiful and there is still no convergence. $\endgroup$ – Mariel Feb 28 at 15:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.