Solving Stefan's solidification problem - for the case of 3 regions
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.
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__] := NDSolve`FiniteDifferenceDerivative@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,
MapThread[
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,
arguments, heads, newVarsSolved},
pos = Flatten[
Outer[Position, functions, oldVars], {{1}, {2}, {3, 4}}];
heads = functions[[All, 0]];
arguments = List @@@ functions;
newVarsSolved = newVars /. Solve[transformations, newVars][[1]];
functionsReplacements =
Map[Function[i,
heads[[i]] -> (Function[#, #2] &[arguments[[i]],
ReplacePart[functions[[i]],
Thread[pos[[i]] -> newVarsSolved]]])],
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[]];
transformation = Thread[newVars == mapping];
{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[
Developer`ToPackedArray@#["ValuesOnGrid"] & /@ sollst3 //
Transpose, {Flatten@sollst3[[1]]["Grid"], grid}];
sollst1 =
NDSolveValue[{eqs, eq0}, P1 /@ grid // Flatten, {t, 0, tend}];
sol1 = ListInterpolation[
Developer`ToPackedArray@#["ValuesOnGrid"] & /@ sollst1 //
Transpose, {Flatten@sollst1[[1]]["Grid"], grid}];
sollst2 =
NDSolveValue[{eqs, eq0}, P2 /@ grid // Flatten, {t, 0, tend}];
sol2 = ListInterpolation[
Developer`ToPackedArray@#["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", ""}]}
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["NDSolve`FEM`"]
(*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 = NDSolve`VariableData[{"DependentVariables" -> {T},
"Space" -> {x}, "Time" -> t}];
nr = ToNumericalRegion[FullRegion[1], {{0, .01}}];
sd = NDSolve`SolutionData[{"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]]*
IdentityMatrix[1]}}, "LoadCoefficients" -> {{ll}}];
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},
NDSolve`SetSolutionDataComponent[sd, "Time", t];
NDSolve`SetSolutionDataComponent[sd, "DependentVariables", u];
l = sdpde["LoadVector"];
s = sdpde["StiffnessMatrix"];
d = sdpde["DampingMatrix"];
tdpde = DiscretizePDE[initCoeffs, methodData, sd, "Transient"];
tbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd,
"Transient"];
l += tdpde["LoadVector"];
s += tdpde["StiffnessMatrix"];
d += tdpde["DampingMatrix"];
nldpde = DiscretizePDE[initCoeffs, methodData, sd, "Nonlinear"];
nlbcs =
DiscretizeBoundaryConditions[initBCs, methodData, sd,
"Nonlinear"];
l += nldpde["LoadVector"];
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},
NDSolve`SetSolutionDataComponent[sd, "Time", t];
NDSolve`SetSolutionDataComponent[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!