6
$\begingroup$

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 \begin{equation} \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] \end{equation}

\begin{equation} \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] \end{equation}

\begin{equation} \frac{\partial T_1}{\partial x}-\frac{\partial T_2}{\partial x}=a\frac{\partial S}{\partial t} \end{equation}

Along with the BC's and IC's: \begin{equation} 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 \end{equation} \begin{equation} T_2\left(L,t\right)=0, \quad T_1\left(x,0\right)=0, \quad T_2\left(x,0\right)=0 \end{equation}

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<T_s$ 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<Ts$.

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.

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – xzczd Oct 30 '18 at 12:09
  • $\begingroup$ 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. $\endgroup$ – Ofek Peretz Oct 30 '18 at 12:47
  • $\begingroup$ 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! $\endgroup$ – Ofek Peretz Oct 30 '18 at 12:53
  • 2
    $\begingroup$ There must exist multiple typos in your post, please double check it. For example, where's $\alpha_2$? $\endgroup$ – xzczd Oct 30 '18 at 13:06
  • 1
    $\begingroup$ @OfekPeretz Thanks for the explanation. Apparently, it is necessary to introduce a second moving boundary. $\endgroup$ – Alex Trounev Oct 30 '18 at 13:36
6
$\begingroup$

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}

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

fig2 fig3 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", ""}]}

fig4

$\endgroup$
12
$\begingroup$

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. IcePhaseCurve

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

Phase Transition Simulated

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]

Simulation Results Corrected

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]

Phase boundary plot

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

Variable Thermal Conductivity Plot

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]

Mushy boundaries

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.

Constant Thermal Conductivity Plot

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

$\endgroup$
  • $\begingroup$ 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}] $\endgroup$ – Alex Trounev Dec 11 '18 at 16:45
  • $\begingroup$ @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. $\endgroup$ – Tim Laska Dec 11 '18 at 16:59
  • $\begingroup$ Tim, thank you! $\endgroup$ – Alex Trounev Dec 11 '18 at 17:55
  • $\begingroup$ @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. $\endgroup$ – Tim Laska Dec 12 '18 at 12:41
  • $\begingroup$ Ok! Now it looks like my solution. But I have to add a second moving border. $\endgroup$ – Alex Trounev Dec 12 '18 at 12:50

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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