Solving Stefan's solidification problem - for the case of 3 regions

This question heavily related to this question, where the case of two PDE's are solved along with a zipping condition that is a function of time.

Using the link in the code I have solved this set of equations $$$$\frac{\partial T_1}{\partial t}=\alpha_1\frac{\partial^2 T_1}{\partial x^2},\quad x\in\left[S\left(t\right),L\right]$$$$

$$$$\frac{\partial T_1}{\partial t}=\alpha_2\frac{\partial^2 T_1}{\partial x^2},\quad x\in\left[0,S\left(t\right)\right]$$$$

$$$$\frac{\partial T_1}{\partial x}-\frac{\partial T_2}{\partial x}=a\frac{\partial S}{\partial t}$$$$

Along with the BC's and IC's: $$$$T_1\left(0,t>0\right)=1, \quad T_1\left(S\left(t\right),t\right)=T_s, \quad T_2\left(S\left(t\right),t\right)=T_s, \quad$$$$ $$$$T_2\left(L,t\right)=0, \quad T_1\left(x,0\right)=0, \quad T_2\left(x,0\right)=0$$$$

where $$\alpha_1,\alpha_2,L,a$$ ,and $$T_s$$ are known constants.

This equation is known as Stefan's solidification problem and the heat coefficients are different in both solid ($$\alpha_2$$) and liquid ($$\alpha_1$$) stats.

I wish to solve the case where after some time, say $$t=t^*$$,the boundary condition at $$x=0$$ is changed to $$T_1(0,t>t^*)=0$$. The solution assumes $$T_1>T_s$$ and $$T_2 at any time, which is now wrong.

I think there are two possible ways to construct that (none of them worked for me):

1) expand the solution, presented by @xzczd in the link above to consist 3 different PDEs, on for each region $$(\alpha_1,\alpha_2,\alpha_1)$$.

2) use the same solution @xzczd presented and adding a piece-wise function to alpha as $$\alpha=\alpha_1$$ when $$T>Ts$$ and $$\alpha=\alpha_2$$ when $$T.

If anyone has experience with a related problem or any idea for how to make it work I will be grateful!

Following the comments to this post I add the post to contain the code I have used for this case:

xR = 10; tend = 200; Pxs = 0.5; \[Kappa] = -10;
hyd1 = 2; h1 = 2;
w = 1;
hyd2 = 1; h2 = 1;
k1 = 1; k2 = hyd2^3/hyd1^3; k3 = k1;
\[Alpha]1 = 1; \[Alpha]2 = hyd1^3/hyd2^3;
eqL = D[P1[x, t], t] == \[Alpha]2 D[P1[x, t], x, x];
eqM = D[P2[x, t], t] == \[Alpha]2 D[P2[x, t], x, x];
eqR = D[P3[x, t], t] == \[Alpha]2 D[P3[x, t], x, x];
{icL, icM, icR} = {P1[x, 0] == 1, P2[x, 0] == 0, P3[x, 0] == 0};
{bcL, bcM,bcR} = {{P1[0, t] == 1, P1[S1@(t), t] == Pxs}, {P2[S1@(t), t] == Pxs, P2[S2@(t), t] == Pxs}, {P3[S2@(t), t] == Pxs, P3[xR, t] == 0}};
slopeL = ( k1 D[P1[x, t], x] /. x -> S1[t]);
slopeM1 = ( k2 D[P2[x, t], x] /. x -> S1[t]);
bcmidfunc1 = S1'[t] == (# - #2)/\[Kappa] &;
slopeM2 = ( k2 D[P2[x, t], x] /. x -> S2[t]);
slopeR = ( k3 D[P3[x, t], x] /. x -> S2[t]);
bcmidfunc2 = S2'[t] == - ((# - #2)/\[Kappa]) &;
With[{eps = 10^-10}, icmid1 = S1[0] == eps];
With[{eps = 10^-10}, icmid2 = S2[0] == 1 + eps];
changeL = DChange[#, x/S1@t == \[Xi], x, \[Xi], P1[x, t]] &;
changeM = DChange[#, (x - S1[t])/(S2[t] - S1[t]) == \[Xi], x, \[Xi], P2[x,t]] &;
changeR = DChange[#, (x - S2[t])/(xR - S2[t]) == \[Xi], x, \[Xi], 3[x,t]] &;
toode = With[{sf = 100}, Map[sf # + D[#, t] &, #, {2}] &];
{neweqL, newicL, newbcL, newslopeL} = changeL@{eqL, icL, toode@bcL,slopeL} /. S1[0] -> S1[t];
{neweqM1, newicM1, newbcM1, newslopeM1} = changeM@{eqM, icM, toode@bcM, slopeM1} /. {S1[0] -> S1[t]} // Simplify;
{neweqM2, newicM2, newbcM2, newslopeM2} = changeM@{eqM, icM, toode@bcM, slopeM2} /. {S2[0] -> S2[t]} // Simplify;
{neweqR, newicR, newbcR, newslopeR} = changeR@{eqR, icR, toode@bcR, slopeR} /. S2[0] -> S2[t] // Simplify;
points = 25;
xdifforder = 4;
{\[Xi]L, \[Xi]R} = domain = {0, 1};
grid = Array[# &, points, domain];
ptoo = pdetoode[{P1, P2, P3}[\[Xi], t], t, grid, xdifforder];
del = #[[2 ;; -2]] &;
{odeL, odeicL, odebcL, odeslopeL, odeM1, odeicM1, odebcM1, odeslopeM1, deM2, odeicM2, odebcM2, odeslopeM2, odeR, odeicR, odebcR, odeslopeR} = MapAt[del, ptoo /@ {neweqL, newicL, newbcL, newslopeL, neweqM1, newicM1, newbcM1, newslopeM1, neweqM2, newicM2, newbcM2, newslopeM2, neweqR, newicR, newbcR, newslopeR}, {{1}, {5}}];
odebcmid1 = bcmidfunc1[odeslopeL, odeslopeM1];
odebcmid2 = bcmidfunc2[odeslopeM2, odeslopeR];
odeicmid1 = icmid1;
odeicmid2 = icmid2;
sollst = NDSolveValue[{odeL, odeM1, odeM2, odeR, odeicL, odeicM1, odeicM2, odeicR, odebcL, odebcM1, odebcM2, odebcR, odebcmid1, odebcmid2, odeicmid1, odeicmid2}, {P1 /@ grid, P2 /@ grid, P3 /@ grid, S1, S2} // Flatten, {t, 0, tend}];


The Error i get is

    NDSolveValue::overdet: There are fewer dependent variables,...., than equations, so the system is overdetermined.


Thanks, Ofek.

• Well, I have to say the question is rather unclear, for example, where's the equation for $T_3$? Also, please learn to format the question properly by using e.g. LaTeX. Oct 30 '18 at 12:09
• Hi @xzczd, thanks for the comments, I've now changed the text into latex format, and fixed some typos that might have caused it to be unclear. Oct 30 '18 at 12:47
• For the other thing you have mentioned: I'm not sure that needs to be an equation for region 3. The issue is that in the link I have mentioned, the tempreture in region 1 is $T>T_s$ and therefore $\alpha=\alpha_1$ and in region 2 $T<T_s$ and therefore $\alpha=\alpha_2$. In the case I have presented, a change in the boundary condition will lower the temperature in region 1 and therefore we get 3 regions: $(\alpha_1,\alpha_2,\alpha_1)$ hope that's a bit clearer now, and thanks once again for your help! Oct 30 '18 at 12:53
• There must exist multiple typos in your post, please double check it. For example, where's $\alpha_2$? Oct 30 '18 at 13:06
• @OfekPeretz Thanks for the explanation. Apparently, it is necessary to introduce a second moving boundary. Oct 30 '18 at 13:36

There is a solution to the problem when the value of the function P1[x,t] on the border does not fall to zero, but to a critical value Pxs. Then the second moving boundary is not needed, and the solution can be obtained using NDSolve (my algorithm)

Nmax = 500; xR = 10; t0 = 1/3; tm = 100; Pxs = 1/2; s0 = 1/10; u0 =
2/100;

S1[0] = s0; S1[-1] = S1[0];
f[t_] := If[t < tm, 1, Pxs]
f1[t_] := If[t < 10^-6, 0, f[t]]
f2[t_] := If[t < 10^-6, 0, Pxs]
p1[0] = NDSolveValue[{D[P1[x, t], t] == D[P1[x, t], x, x],
P1[0, t] == f1[t], P1[S1[0], t] == f2[t], P1[x, 0] == 0},
P1, {x, 0, S1[0]}, {t, 0, t0}];
p2[0] = NDSolveValue[{D[P2[x, t], t] == D[P2[x, t], x, x],
P2[xR, t] == 0, P2[S1[0], t] == f2[t], P2[x, 0] == 0},
P2, {x, S1[0], xR}, {t, 0, t0}];
p3[0][x_, t_] := 0

Table[T[i] = t0 + i*t0, {i, 0, Nmax}];

Do[p1[i] =
NDSolveValue[{D[P1[x, t], t] == D[P1[x, t], x, x],
P1[0, t] == f1[t], P1[S1[i - 1], t] == f2[t],
P1[x, T[i - 1]] ==
If[x < S1[i - 2], p1[i - 1][x, T[i - 1]], f2[T[i - 1]]]},
P1, {x, 0, S1[i - 1]}, {t, T[i - 1], T[i]}] ;
p2[i] =
NDSolveValue[{D[P2[x, t], t] == D[P2[x, t], x, x], P2[xR, t] == 0,
P2[S1[i - 1], t] == f2[T[i - 1]],
P2[x, T[i - 1]] == p2[i - 1][x, T[i - 1]]},
P2, {x, S1[i - 1], xR}, {t, T[i - 1], T[i]}];
S1[i] = S1[i - 1] -
u0*( Derivative[1, 0][p1[i]][S1[i - 1], T[i]] -
Derivative[1, 0][p2[i]][S1[i - 1], T[i]]), {i, 1,
Nmax}] // Quiet
P12 = Interpolation[
Flatten[Table[{{x, T[i]},
If[x < S1[i - 1], p1[i][x, T[i]], p2[i][x, T[i]]]}, {x, .01,
xR, .01}, {i, 0, Nmax}], 1], InterpolationOrder -> 3]; // Quiet

pp1[i_] :=
Table[{x, If[x < S1[i - 1], p1[i][x, T[i]], p2[i][x, T[i]]]}, {x, 0,
xR, .01}];

{ListLinePlot[Table[pp1[i], {i, 0, 290, 10}], PlotRange -> All,
AxesLabel -> {"x", "P"}],
ListLinePlot[Table[pp1[i], {i, 300, 320, 1}], PlotRange -> All,
AxesLabel -> {"x", "P"}],
ListLinePlot[Table[pp1[i], {i, 320, Nmax, 10}], PlotRange -> All,
AxesLabel -> {"x", "P"}]}
{ListPlot[Table[{T[i], S1[i]}, {i, 0, Nmax}],
AxesLabel -> {"t", "S1"}],
Plot3D[P12[x, t], {x, 0, xR}, {t, 0, T[Nmax]}, PlotRange -> {0, 1},
Mesh -> None, ColorFunction -> "TemperatureMap", PlotPoints -> 40,
AxesLabel -> {"x", "t", "P"}] // Quiet}


Unexpectedly, I managed to debug the code that @Ofek Peretz wrote. Corrections had to be done a lot. The result is the curve that I received in part using my code, and Tim Laska received in full. The solution uses two wellknown functions that I will present for the convenience of readers. The first function is called pdetoode introduced by @xzczd and DChange introduced by @Kuba

Clear[fdd, pdetoode, tooderule]
fdd[{}, grid_, value_, order_, periodic_] := value;
fdd[a__] := NDSolveFiniteDifferenceDerivative@a;

pdetoode[funcvalue_List, rest__] :=
pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]],
rest];
pdetoode[{func__}[var__], rest__] :=
pdetoode[Alternatives[func][var], rest];
pdetoode[front__, grid_?VectorQ, o_Integer, periodic_: False] :=
pdetoode[front, {grid}, o, periodic];

pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer,
periodic : True | False | {(True | False) ..} : False] :=
With[{pos = Position[{var}, time][[1, 1]]},
With[{bound = #[[{1, -1}]] & /@ {grid},
pat = Repeated[_, {pos - 1}],
spacevar = Alternatives @@ Delete[{var}, pos]},
With[{coordtoindex =
Function[coord,
Piecewise[{{1, # === #2[[1]]}, {-1, # === #2[[-1]]}},
All] &, {coord, bound}]]},
tooderule@
Flatten@{((u : func) |
Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat,
t_, x2___] :> (Sow@coordtoindex@{x1, x2};

fdd[{dx1, dx2}, {grid},
Outer[Derivative[dt][u@##]@t &, grid],
"DifferenceOrder" -> o,
PeriodicInterpolation -> periodic]),
inde : spacevar :>
With[{i = Position[spacevar, inde][[1, 1]]},
Outer[Slot@i &, grid]]}]]];

tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] :=
Equal[tooderule[rule][a - b], 0] //.
eqn : HoldPattern@Equal[_, _] :> Thread@eqn;
tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@
Reap[expr /. rule]

Clear@pdetoae;
pdetoae[funcvalue_List, rest__] :=
pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoae[{func__}[var__], rest__] :=
pdetoae[Alternatives[func][var], rest];

pdetoae[func_[var__], rest__] :=
Module[{t},
Function[pde, #[
pde /. {Derivative[d__][u : func][inde__] :>
Derivative[d, 0][u][inde, t], (u : func)[inde__] :>
u[inde, t]}] /. (u : func)[i__][t] :> u[i]] &@
pdetoode[func[var, t], t, rest]]
ClearAll[DChange];

DChange[expr_, transformations_List, oldVars_List, newVars_List,
functions_List] :=
Module[{pos, functionsReplacements, variablesReplacements,
pos = Flatten[
Outer[Position, functions, oldVars], {{1}, {2}, {3, 4}}];
arguments = List @@@ functions;
newVarsSolved = newVars /. Solve[transformations, newVars][[1]];
functionsReplacements =
Map[Function[i,
ReplacePart[functions[[i]],
Range@Length@functions];
variablesReplacements = Solve[transformations, oldVars][[1]];
expr /. functionsReplacements /. variablesReplacements //
Simplify // Normal];

DChange[expr_, functions : {(_[___] == _) ..}] :=
expr /. Replace[
functions, (f_[vars__] == body_) :> (f ->
Function[{vars}, body]), {1}]

DChange[expr_, x___] :=
DChange[expr, ##] & @@
Replace[{x}, var : Except[_List] :> {var}, {1}];

DChange[expr_, coordinates : Verbatim[Rule][__String], oldVars_List,
newVars_List, functions_] :=
Module[{mapping, transformation},
mapping =
Check[CoordinateTransformData[coordinates, "Mapping", oldVars],
Abort[]];
{DChange[expr, transformation, oldVars, newVars, functions],
transformation}];


Below is the debugged code from which @Ofek Peretz began the discussion

xR = 10; tend = 286.288; t0 = 200; Pxs = 0.5; \[Kappa] = 10;
hyd1 = 1; h1 = 2;
w = 1;
hyd2 = 1; h2 = 1;
k1 = 1; k2 = hyd2^3/hyd1^3; k3 = k1;
\[Alpha]1 = 1; \[Alpha]2 = hyd1^3/hyd2^3;
fp10[t_] := (1 - Tanh[20*(t - t0)])/2
fp11[t_] := (1 + Pxs - (1 - Pxs)*Tanh[20*(t - t0)])/2
fs1[t_] := (1 + Tanh[20*(t - t0)])/2
eqL = D[P1[x, t], t] == \[Alpha]2 D[P1[x, t], x, x];
eqM = D[P2[x, t], t] == \[Alpha]2 D[P2[x, t], x, x];
eqR = D[P3[x, t], t] == \[Alpha]2 D[P3[x, t], x, x];
{icL, icM, icR} = {P1[x, 0] == 1, P2[x, 0] == 1, P3[x, 0] == 0};
{bcL, bcM,
bcR} = {{P1[0, t] == 1, P1[1, t] == 1}, {P2[0, t] == 1,
P2[1, t] == Pxs}, {P3[0, t] == Pxs, P3[1, t] == 0}};

With[{eps = 10^-3}, icmid1 = S1[0] == eps];
With[{eps = 2*10^-3}, icmid2 = S2[0] == eps];
changeL = DChange[#, x/S1@t == \[Xi], x, \[Xi], P1[x, t]] &;
changeM =
DChange[#, (2*x - S1[t] - S2[t])/(S2[t] - S1[t])/2 + 1/2 == \[Xi],
x, \[Xi], P2[x, t]] &;
changeR =
DChange[#, (x - S2[t])/(xR - S2[t]) == \[Xi], x, \[Xi], P3[x, t]] &;
toode = With[{sf = 100}, Map[sf # + D[#, t] &, #, {2}] &];
neweqL = changeL@eqL // Simplify;
neweqM1 = changeM@eqM // Simplify;

neweqR = changeR@eqR // Simplify;
bc = {newbcL, newbcM1, newbcR} = {toode@bcL, toode@bcM, toode@bcR};
points = 25; h = 1/(points - 1);
xdifforder = 2;
{\[Xi]L, \[Xi]R} = domain = {0, 1};
grid = Array[# &, points, domain];
ptoo = pdetoode[{P1[\[Xi], t], P2[\[Xi], t], P3[\[Xi], t]}, t, grid,
xdifforder];
del = #[[2 ;; -2]] &;
{newicL, newicM1,
newicR} = {{P1[\[Xi], 0] == 1}, {P2[\[Xi], 0] ==
1}, {P3[\[Xi], 0] == 0}};
{odeL, odeicL, odebcL, odeM1, odeicM1, odebcM1, odeR, odeicR,
odebcR} = {neweqL // ptoo // del, del /@ ptoo@newicL,
newbcL // ptoo, neweqM1 // ptoo // del, del /@ ptoo@newicM1,
newbcM1 // ptoo, neweqR // ptoo // del, del /@ ptoo@newicR,
newbcR // ptoo};
eqS1 = S1'[t] ==
fs1[t]*(k1*Derivative[1, 0][P1][1, t]/S1[t] -
k2*Derivative[1, 0][P2][0, t]/(S2[t] - S1[t]))/\[Kappa];
eqS2 = S2'[
t] == (k3*Derivative[1, 0][P3][0, t]/(xR - S2[t]) -
k2*Derivative[1, 0][P2][1, t]/(S2[t] - S1[t]))/\[Kappa];
odeicmid1 = icmid1;
odeicmid2 = icmid2;
op = {
\!$$\*SuperscriptBox[\(P1$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[1, t] -> (P1[1][t] - P1[1 - h][t])/h,
\!$$\*SuperscriptBox[\(P2$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[0, t] -> (P2[h][t] - P2[0][t])/h,
\!$$\*SuperscriptBox[\(P2$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[1, t] -> (P2[1][t] - P2[1 - h][t])/h,
\!$$\*SuperscriptBox[\(P3$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[0, t] -> (P3[h][t] - P3[0][t])/h};
eq0 = {P1[0][0] == 1, P2[0][0] == 1, P3[0][0] == Pxs, P1[1][0] == 1,
P2[1][0] == Pxs, P3[1][0] == 0};

eqsop = {odeL, odeM1, odeR, odeicL, odeicM1, odeicR, odebcL, odebcM1,
odebcR, odeicmid1, odeicmid2, eqS1, eqS2} /. op;
op1 = {eqsop[[
7]] -> {-100*fp10[t] + 100 P1[0][t] + Derivative[1][P1[0]][t] ==
0, -100*fp11[t] + 100 P1[1][t] + Derivative[1][P1[1]][t] == 0},
eqsop[[8]] -> {-100*fp11[t] + 100 P2[0][t] +
Derivative[1][P2[0]][t] ==
0, -50. + 100 P2[1][t] + Derivative[1][P2[1]][t] == 0}};

eqs = eqsop /. op1;

sollst3 =
NDSolveValue[{eqs, eq0}, P3 /@ grid // Flatten, {t, 0, tend}];

sol3 = ListInterpolation[
DeveloperToPackedArray@#["ValuesOnGrid"] & /@ sollst3 //
Transpose, {Flatten@sollst3[[1]]["Grid"], grid}];
sollst1 =
NDSolveValue[{eqs, eq0}, P1 /@ grid // Flatten, {t, 0, tend}];
sol1 = ListInterpolation[
DeveloperToPackedArray@#["ValuesOnGrid"] & /@ sollst1 //
Transpose, {Flatten@sollst1[[1]]["Grid"], grid}];
sollst2 =
NDSolveValue[{eqs, eq0}, P2 /@ grid // Flatten, {t, 0, tend}];
sol2 = ListInterpolation[
DeveloperToPackedArray@#["ValuesOnGrid"] & /@ sollst2 //
Transpose, {Flatten@sollst2[[1]]["Grid"], grid}];

{ContourPlot[sol1[t, x], {t, 0, tend}, {x, 0, 1}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> {"t", "\[Xi]"},
PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "P1"],
ContourPlot[sol2[t, x], {t, 0, tend}, {x, 0, 1}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> {"t", "\[Xi]"},
PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "P2"],
ContourPlot[sol3[t, x], {t, 0, tend}, {x, 0, 1}, Contours -> 20,
ColorFunction -> Hue, FrameLabel -> {"t", "\[Xi]"},
PlotLegends -> Automatic, PlotRange -> All, PlotLabel -> "P3"]}
S11 = NDSolveValue[{eqs, eq0}, S1, {t, 0, tend}]; S21 =
NDSolveValue[{eqs, eq0}, S2, {t, 0, tend}];
Plot[{S11[t], S21[t]}, {t, 0, tend}, AxesLabel -> {"t", "S1,S2"}]


Add a couple of pictures for clarity of what P1,P2,P3 looks like in the coordinates x, t.

 F[t_, x_] :=
If[x >= S21[t], sol3[t, (x - S21[t])/(xR - S21[t])],
If[x <= S11[t], sol1[t, x/S11[t]],
sol2[t, (2*x - S11[t] - S21[t])/(S21[t] - S11[t])/2 + 1/2]]]
l1 = Table[{t, x, F[t, x]}, {t, 0, 1.9, .05}, {x, 0, 5, .1}];
l2 = Table[{t, x, F[t, x]}, {t, 2, tend, .5}, {x, 0, 5, .1}];
ff = Interpolation[Join[Flatten[l1, 1], Flatten[l2, 1]],
InterpolationOrder -> 4];

{ContourPlot[ff[t, x], {t, 0, tend}, {x, 0, 5}, Contours -> 38,
PlotPoints -> 50, FrameLabel -> {"t", "x"},
PlotLegends -> Automatic, ColorFunction -> "TemperatureMap"],
ListPlot3D[Join[Flatten[l1, 1], Flatten[l2, 1]], Mesh -> None,
InterpolationOrder -> 4, ColorFunction -> Hue,
AxesLabel -> {"t", "x", ""}]}


Important Update

Through a copy-paste mistake, I dropped the temperature dependence on the effective heat capacity, which was the point of the exercise. I have modified the post to correct that mistake.

Update #2

The difference between liquid and solid thermal conductivity can be quite large (4x for water), so I added support for the temperature dependent thermal conductivity.

Modified Original Post

You might find it easier to introduce the concept of a "mushy" zone versus tracking the phase boundaries. The typical enthalpy versus temperature curve for the melting of ice is shown in the following.

We can approximate this curve with continuous functions and not lose much in the simulation. Typically, this is done by defining an apparent heat capacity like so

$${\hat C_{p}(T)}={\hat C_{p,s}}\left ( 1-\epsilon(T) \right )+\epsilon(T){\hat C_{p,l}} + \Delta {\hat H_f}\frac{d \epsilon}{dT}$$

Where $$\epsilon$$ is the phase fraction of liquid. Any nice shaped S-curve will do so let's try a normalized Tanh function that ranges from 0 to 1. $$T_m$$ is the meltpoint and $$\delta$$ is the width of the phase transition.

$$\epsilon(T)=\frac{{1 + \operatorname{Tanh} \left( {8\frac{{\left( {T - {T_m}} \right)}}{\delta }} \right)}}{2}$$

And the temperature derivative of $$\epsilon$$ is

$$\frac{d\epsilon}{dT}=\frac{{4\operatorname{Sech} {{\left( {8\frac{{\left( {T - {T_m}} \right)}}{\delta }} \right)}^2}}}{\delta }$$

Now, we can integrate $${\hat C_{p}(T)}$$ to see that we recover the enthalpy phase transition curve.

us[m_, t_] := (1 + Tanh[m (t)])/2(* Unit Step *)
eps[tm_, deltaT_, t_] := (1 + Tanh[8 (t - tm)/deltaT])/2
deps[tm_, deltaT_, t_] := (4 Sech[(8 (t - tm))/deltaT]^2)/deltaT
cpeff[deltaHf_, cpl_, cps_, tm_, deltaT_, t_] :=
cps (1 - eps[tm, deltaT, t]) + cpl eps[tm, deltaT, t] +
deltaHf deps[tm, deltaT, t]
keff[kl_, ks_, tm_, deltaT_, t_] :=
ks (1 - eps[tm, deltaT, t]) + kl  eps[tm, deltaT, t]
deltah[deltaHf_, cpl_, cps_, tm_, deltaT_, t_] :=
1/16 (8 (cpl + cps) t + (cpl - cps) deltaT Log[
Cosh[(8 (t - tm))/deltaT]] + 8 deltaHf Tanh[(8 (t - tm))/deltaT])
pulseSeq[a_, fac_, m_, t0_, t_] := (
a (1 + (1 + fac) us[m, t0 - t] (-1 + us[m, -t]) + (-1 + fac) us[
m, -t]))/fac


Using this phase fraction approach, we only have one temperature equation to solve as shown.

$$\frac{\partial T}{\partial t}=\frac{k}{\rho \hat C_{p}(T)} \frac{\partial^2 T}{\partial x^2}$$

In the following, I start with water at $$5^\circ C$$ and "instantaneously" imposes a $$-5^\circ C$$ wall. After 1800 seconds, the wall returns to $$5^\circ C$$. The Mathematica code is shown below (note that earlier version of this post contained an error in cpeff).

eqn = D[u[t, x], t] -
0.555/(1000 cpeff[334000., 4180, 2060, 0, 0.2, u[t, x]])
D[u[t, x], {x, 2}] == 0;
opts = Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 500}};
ics = {u[0, x] == 5};
bcs = {u[t, 0] == pulseSeq[5, 1, 10, 1800, t], u[t, 0.05] == 5};
eqns = {eqn} ~ Join ~ ics ~ Join ~ bcs;
ifun = NDSolveValue[eqns, {u}, {t, 0, 3600}, {x, 0, 0.05}, opts];
plts = Plot3D[ifun[[1]][t, x], {t, 0, 3600}, {x, 0, 0.0075},
PlotRange -> All, PlotPoints -> {200, 200}, MaxRecursion -> 4,
ColorFunction -> (ColorData["MintColors"][#3] &),
MeshFunctions -> {If[-0.01 < #3 < 0.01, #3, 0] &,
If[#3 < 0, #3, 0] &, If[#3 > 0, #3, 0] &}, Mesh -> 21,
AxesLabel -> Automatic,
MeshStyle -> {{Black, Thick}, {Blue, Thick, Dashed}, {Red, Thick,
DotDashed}}, Boxed -> False, ImageSize -> Large];
sz = 300;
Grid[{{Show[{plts}, ViewProjection -> "Orthographic",
ViewPoint -> Front, ImageSize -> sz,
Background -> RGBColor[0.84, 0.92, 1.], Boxed -> False],
Show[{plts}, ViewProjection -> "Orthographic", ViewPoint -> Right,
ImageSize -> sz, Background -> RGBColor[0.84, 0.92, 1.],
Boxed -> False]}, {Show[{plts}, ViewProjection -> "Orthographic",
ViewPoint -> Top, ImageSize -> sz,
Background -> RGBColor[0.84, 0.92, 1.], Boxed -> False],
Show[{plts}, ViewProjection -> "Perspective",
ViewPoint -> {Above, Left, Front}, ImageSize -> sz,
Background -> RGBColor[0.84, 0.92, 1.], Boxed -> False]}},
Dividers -> Center]


The model seems to behave reasonably well. There is nice parabolic growth in the ice layer initially. Once the heating is turned back on, the ice first melts at the wall as expected. When all the ice is gone, the system starts to heat more rapidly as on would expect. An elegant way to produce the phase boundary was suggested by @Alex Trounev.

ContourPlot[ifun[[1]][t, x] == 0, {t, 0, 3600}, {x, 0, 0.01}, MaxRecursion -> 4]


Enhancement for Temperature Dependent Thermal Conductivity

Ice is about 4x more thermally conductive than liquid water. In my initial analysis, the thermal conductivity, $$k$$, was presumed to be constant. The original energy balance was given by

$$\frac{{\partial T}}{{\partial t}} - \alpha \frac{{{\partial ^2}T}}{{\partial {x^2}}} = 0$$

To use the "mushy" zone, $$\alpha$$ is no longer constant, but a function of temperature so we really need to recast it in coefficient form as shown below.

$$m\frac{{{\partial ^2}}}{{\partial {t^2}}}u + d\frac{\partial }{{\partial t}}u + \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = 0$$

All terms except $$d$$ and $$c$$ are zero so the new energy balance becomes:

$$\rho {{\hat C}_p}(T)\frac{\partial u }{{\partial t}} + \nabla \cdot \left( { - k(T)\nabla u} \right) = \rho {{\hat C}_p}(T)\frac{\partial u }{{\partial t}} + \frac{\partial }{{\partial x}}\left( { - k(T)\frac{\partial u }{{\partial x}}} \right) = 0$$

Where $$k(T)$$ will be given by the same s-shaped curve $$\epsilon(T)$$ or

$$k(T) = k_s\left ( 1-\epsilon(T) \right )+k_l\epsilon(T)$$

The Mathematica code for a $$-4^\circ C$$ to $$4^\circ C$$ pulse sequence (sans plotting that is the same as above) is shown below. I needed to supply a small MaxStepFraction to capture the change in wall temperature at 1800s. This can lead to large file sizes. Perhaps, someone can comment on better/more robust settings.

rho = 1000;
dhf = 334000.;
cpl = 4180.0;
cps = 2060;
kl = 0.55575;
ks = 2.22;
tm = 0;
delta = 0.25;
eqn = rho cpeff[dhf, cpl, cps, tm, delta, u[t, x]] D[u[t, x], t] -
D[keff[kl, ks, tm, delta, u[t, x]] D[u[t, x], x], x] == 0;
opts = Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 200}};
ics = {u[0, x] == 4};
bcs = {u[t, 0] == pulseSeq[4, 1, 0.0625, 1800, t], u[t, 0.01] == 4};
eqns = {eqn} ~ Join ~ ics ~ Join ~ bcs;
ifun = NDSolveValue[eqns, {u}, {t, 0, 3600}, {x, 0, 0.01}, opts,
MaxStepFraction -> 1/2000];


Variable Conductivity

The "mushy" zone is fuzzy so I created a contour plot of the $$\pm 30\%$$ bands using the following Mathematica code.

ContourPlot[{ifun[[1]][t, x],
ifun[[1]][t, x] == -delta/11.54156032711171, ifun[[1]][t, x] == 0,
ifun[[1]][t, x] == delta/11.54156032711171}, {t, 0, 3600}, {x, 0,
0.01}, MaxRecursion -> 4]


We expect the spacing between solid and median temperature band to be about 4x that of the liquid and median band because there is 4x the resistance through the liquid layer. The plot suggests that this is the case.

Constant Conductivity

I changed the pulse sequence, boundary conditions, and the domain dimensions from the previous example. For completeness, I re-ran the constant liquid conductivity case with the geometry and conditions.

Because the liquid conductivity is lower than the solid, the film does not grow as fast as the variable conductivity case.

• Tim, you created excellent code. Can be used to highlight the phase transition boundary ContourPlot[ifun[[1]][t, x] == 0, {t, 0, 3600}, {x, 0, 0.05}] Dec 11 '18 at 16:45
• @AlexTrounev Thank you very much Alex! I appreciate your work. Your method to display the phase transition boundary is much much simpler. I presume that you won't mind if I update the post. Dec 11 '18 at 16:59
• Tim, thank you! Dec 11 '18 at 17:55
• @AlexTrounev I screwed up the original code and dropped the temperature dependence of the heat capacity (a significant error). I have corrected the code in the post. Dec 12 '18 at 12:41
• Ok! Now it looks like my solution. But I have to add a second moving border. Dec 12 '18 at 12:50

Inspired by Tim's excellent answer I was wondering if it is not possible to use the FEM for this. In Version 12.0 there is no way to do this from top level. However, we still can solve this phase change problem with the low level FEM functions, which I'd like to show here. Why bother? In this specific case we look at a 1D spatial plus time PDE and there is not much to be gained by making use of the FEM. In two spatial dimensions, things are different as the FEM will then allow us to model phase change in arbitrary shaped regions.

Since Tim covered the theory so well in his post I am going to focus on the FEM part. The basic idea is to write a function discretizeResidual the returns the residual of the PDE and then have NDSolve do the time integration. Variants of this approach are documented in the FEMProgramming tutorial.

We start by setting up the parameters.

Needs["NDSolveFEM"]
(*latent heat*)
Lm = 3.335*10^5;
(* melt point *)
Tm = 0;
(*temperature interval for phase change*)
dT = 0.25;
Tc = -5;
Th = 5;
r1 = 918;
r2 = 997;
Cp1 = 2052;
Cp2 = 4179;


Next we set up some derived parameters and some helper functions:

SmoothedStepFunction[fmin_, fmax_, ts_, m_] :=
Function[t, (fmin*Exp[ts*m] + fmax*Exp[m*t])/(Exp[ts*m] + Exp[m*t])]

theta = SmoothedStepFunction[1, 0, 0, 50];

rRatio = r1/r2;
(*convertion from space coord. to material coord*)
r = r2*rRatio;

k1 = 2.31;
(*convertion from space coord. to material coord*)
k2 = 0.613*rRatio;
k[T_] := theta[T]*k1 + (1 - theta[T])*k2


Setting up CL with a mass fraction a

a[T_] := (1 - 2 theta[T])/2;
CL[T_] := Evaluate[Lm*D[a[Temp], Temp] /. {Temp -> T}]


An alternative method is to use a regularized delta function

(* RegularizedDelta[g_,T_,Tm_]:=Piecewise[{{1/(4*g)(1+Cos[Pi/(2*g)(T-Tm)]),RealAbs[T-Tm]\[LessEqual]2*g},{0,True}}];
CL[T_]:=Lm*RegularizedDelta[dT/4,T,Tm]
*)


We perform a basic check that the integral over CL matches the latent heat.

(* check *)
NIntegrate[CL[T], {T, Tc, Th}] == Lm


Set up Cp with latent heat.

Cp[T_] := theta[T]*Cp1 + (1 - theta[T])*Cp2 + CL[T]


For completeness we have Cp without latent heat

(* Cp[T_]:=theta[T]*Cp1+(1-theta[T])*Cp2 *)


Set up the temperature profile on the left boundary:

TL[t_] := Piecewise[{{SmoothedStepFunction[Th, Tc, 0.5, 20][t], t <= 1}, {Tc,
1 < t <= 1800}, {SmoothedStepFunction[Tc, Th, 1800.5, 20][t],
1800 < t <= 1801}, {Th, True}}];

(*Plot[TL[t],{t,0,3600},PlotTheme\[Rule]"Detailed",PlotLabel\[Rule]\
"TLeft",FrameLabel\[Rule]{"time [s]","temperature \
[\[Degree]C]"},PlotLegends\[Rule]None,LabelStyle\[Rule]Directive[14],\
ImageSize\[Rule]Medium,Exclusions\[Rule]None]*)


Now we start to set up the low level FEM functions.

dampCoef[t_, T_] := r*Cp[T]
diffuCoef[t_, T_] := k[T]
ll = 0;
ic = T[0, x] == Th;
bc = {DirichletCondition[T[t, x] == TL[t], x == 0.],
DirichletCondition[T[t, x] == Th, x == 0.01]};
tEnd = 3600;


Initialize the BCs and method data:

vd = NDSolveVariableData[{"DependentVariables" -> {T},
"Space" -> {x}, "Time" -> t}];
nr = ToNumericalRegion[FullRegion[1], {{0, .01}}];
sd = NDSolveSolutionData[{"Space" -> nr, "Time" -> 0.}];
initBCs = InitializeBoundaryConditions[vd, sd, {bc}];
methodData =
InitializePDEMethodData[vd, sd,
Method -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.01/50}}];


We initialize the coefficients:

initCoeffs = InitializePDECoefficients[vd, sd,
"DampingCoefficients" -> {{dampCoef[t, T[t, x]]}},
"DiffusionCoefficients" -> {{-diffuCoef[t, T[t, x]]*


We discretize the stationary components once before everything else and make use of these during the time integration such that we do not have to recompute them at every step.

sdpde = DiscretizePDE[initCoeffs, methodData, sd];
sbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd];


What follows is the heart of the idea. The function discretizePDEResidual gets the time t, dependent variable u and the time derivative of u as an argument. From there we update the solution data with the current time and dependent variable. We extract the stationay components computed previously and add the transient and nonlinear components. After that we deploy the boundary conditons and finally compute the residual.

discretizePDEResidual[t_?NumericQ, u_?VectorQ, dudt_?VectorQ] :=

Module[{l, s, d, tdpde, tbcs, nldpde, nlbcs, sdTemp},

NDSolveSetSolutionDataComponent[sd, "Time", t];
NDSolveSetSolutionDataComponent[sd, "DependentVariables", u];

s = sdpde["StiffnessMatrix"];
d = sdpde["DampingMatrix"];

tdpde = DiscretizePDE[initCoeffs, methodData, sd, "Transient"];
tbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd,
"Transient"];
s += tdpde["StiffnessMatrix"];
d += tdpde["DampingMatrix"];

nldpde = DiscretizePDE[initCoeffs, methodData, sd, "Nonlinear"];
nlbcs =
DiscretizeBoundaryConditions[initBCs, methodData, sd,
"Nonlinear"];
s += nldpde["StiffnessMatrix"];
d += nldpde["DampingMatrix"];

DeployBoundaryConditions[{l, s, d}, nlbcs];
DeployBoundaryConditions[{l, s, d}, tbcs];
DeployBoundaryConditions[{l, s, d}, sbcs];

d.dudt + s.u - l
]


To allow for a fast time integration we write a small help function to construct a damping matrix at a specific point in time. From this damping matrix we extract the sparsity pattern. Details can be found in the FEMProgramming tutorial.

initializeDampingMatrix[t_, u_?VectorQ] := Module[
{d, tdpde, nldpde},

NDSolveSetSolutionDataComponent[sd, "Time", t];
NDSolveSetSolutionDataComponent[sd, "DependentVariables", u];

d = sdpde["DampingMatrix"];

tdpde = DiscretizePDE[initCoeffs, methodData, sd, "Transient"];
d += tdpde["DampingMatrix"];

nldpde = DiscretizePDE[initCoeffs, methodData, sd, "Nonlinear"];
d += nldpde["DampingMatrix"];

d
]


One thing to keep in mind with this function is that we should choose a time that allows transient components to be picked up. The exact time does not matter and it could be t=0 just make sure that possible transient components are picked up.

sparsity = initializeDampingMatrix[10^-6, init]["PatternArray"];


Create the initial conditions:

init = #[[1]] & /@ methodData["ElementMesh"]["Coordinates"];


Do the time integration:

measure =
AbsoluteTiming[
MaxMemoryUsed[
Monitor[tTfun =
NDSolveValue[{discretizePDEResidual[t, T[t], T'[ t]] == 0,
T[0] == init}, T, {t, 0, tEnd}
, Method -> {"EquationSimplification" -> "Residual"}
, Jacobian -> {Automatic, Sparse -> sparsity}
, EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])],
monitor]]/(1024.^2)];
Print["Time -> ", measure[[1]], "\nMemory -> ", measure[[2]]]


On my computer this takes about 130 seconds and 55MB of memory.

Create final the interpolating function:

Tfun = ElementMeshInterpolation[{tTfun["Coordinates"][[1]],
methodData["ElementMesh"]}, Partition[tTfun["ValuesOnGrid"], 1]];


Visualize the result:

TRange = {Tc, Th};
legendBar =
BarLegend[{"TemperatureMap", TRange}, 50,
LegendLabel -> Text[Style["[C]", Opacity[0.6]]]];
options = {PlotRange -> TRange,
ColorFunction -> ColorData[{"TemperatureMap", TRange}],
ColorFunctionScaling -> False, ContourStyle -> Opacity[0.1],
Contours -> 30, PlotPoints -> 50,
FrameLabel -> {"time [s]", "x [m]"}, LabelStyle -> Directive[14],
PlotLabel -> Text[Style["Temperature: T(t,x)", 18]],
ImageSize -> Medium};
tempField =
Legended[ContourPlot[Tfun[t, x], {t, 0, 3600}, {x, 0, 0.01},
Evaluate[options]], legendBar]; phaseBoundary =
ContourPlot[{Tfun[t, x] == 0}, {t, 0, 3600}, {x, 0, 0.01},
MaxRecursion -> 5]; Show[tempField, phaseBoundary]


As many of you know, posts like these, were I get an opportunity to translate the stuff into the low level FEM functions, form the basis for what will be doable in NDSolve`. Once more I'd like to take the opportunity to thank Tim for his great post (and this is certainly not the only one!) Well done and thanks!

• This is a good job and working code (+1). Mar 11 '20 at 13:18