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.

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

fig1

Your Answer

 
discard

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.