8
$\begingroup$

I have a system of PDEs in the following form,

$$\frac{\partial T_1}{\partial t}=\frac{\partial^2 T_1}{\partial x^2},\,\,0<x<S(t)$$ with $T_1(x,0)=-10,\,\,T_1(0,t)=-10,\,\, T_1(S(t),t)=10$, and $$\frac{\partial T_2}{\partial t}+\frac{dS(t)}{dt}\frac{\partial T_2}{\partial x}=\frac{\partial^2 T_2}{\partial x^2},\,\,x>S(t)$$ with $T_2(x,0)=10,\,\,T_2(S(t),t)=10,\,\, T_2(20,t)=20$.

The above two equations are zip together by $$\frac{dS(t)}{dt}=\frac{\partial T_1}{\partial x}\Biggl|_{x=S(t)}-\frac{\partial T_2}{\partial x}\Biggl|_{x=S(t)}$$ with $S(0)=0$.

Eq1 = D[T1[x, t], t] == D[T1[x, t], x, x]
Eq2 = D[T2[x, t], t]+D[S[t], t] D[T2[x, t], x]== D[T2[x, t], x, x]
T1[x, 0] == -10,T1[0,t] == -10, T1[S(t), t] == 10,T2[x, 0] == 10,T2[S(t), t] == 10,T2[20, t] == 20

Is there a way to solve the above system in Mathematica?

$\endgroup$
13
$\begingroup$

This problem is quite similar to PDE with Stefan Conditions, a.k.a variable boundary, and can be solved more easily because now we have DChange and pdetoode.

We first use DChange to transform $T_1(x_1,t)$ and $T_2(x_2,t)$ to $T_1(\xi,t)$ and $T_2(\xi,t)$, where

$$\frac{x_1}{S(t)}=\xi$$ $$\frac{x_2-S(t)}{20-S(t)}=\xi$$

xR = 20;
eqL = D[T1[x, t], t] == D[T1[x, t], x, x];
eqR = D[T2[x, t], t] + D[S[t], t] D[T2[x, t], x] == D[T2[x, t], x, x];
{icL, icR} = {T1[x, 0] == -10, T2[x, 0] == 10};
{bcL, bcR} = {{T1[0, t] == -10, T1[S@(t), t] == 10}, 
              {T2[S@(t), t] == 10, T2[xR, t] == 20}};
slopeL = (D[T1[x, t], x] /. x -> S[t]);
slopeR = (D[T2[x, t], x] /. x -> S[t]);
bcmidfunc = S'[t] == # - #2 &;
With[{eps = 10^-10}, icmid = S[0] == eps];
(* Definition of DChange isn't included in this code piece,
   please find it in the link above *)
changeL = DChange[#, x/S@t == ξ, x, ξ, T1[x, t]] &;
changeR = DChange[#, (x - S[t])/(xR - S[t]) == ξ, x, ξ, T2[x, t]] &;

Clear[toode]
toode[expr_Equal] := With[{sf = 100}, sf # + D[#, t] & /@ expr]
toode[expr_List] := toode /@ expr

{neweqL, newicL, newbcL, newslopeL} = 
 changeL@{eqL, icL, toode@bcL, slopeL} /. S[0] -> S[t]
{neweqR, newicR, newbcR, newslopeR} = 
 changeR@{eqR, icR, toode@bcR, slopeR} /. S[0] -> S[t] // Simplify

Remark

  1. The /. S[0] -> S[t] is added because currently DChange can't handle coefficients involving the independent variable very well.

  2. Check this post if you don't understand why we need toode.

After the change of variable, $T_1$ and $T_2$ are both defined in $[0,1]\times[0,t_{end}]$, while $S$ is defined in $[0,t_{end}]$, so NDSolve still can't solve the equation set directly. We need to discretize the PDEs to a set of ODEs. I'll use pdetoode for the task:

points = 25;
xdifforder = 4;
{ξL, ξR} = domain = {0, 1};
grid = Array[# &, points, domain];

(* Definition of pdetoode isn't included in this code piece,
   please find it in the link above *)
ptoo = pdetoode[{T1, T2}[ξ, t], t, grid, xdifforder];
del = #[[2 ;; -2]] &;

{odeL, odeicL, odebcL, odeslopeL, odeR, odeicR, odebcR, odeslopeR} = 
    MapAt[del, ptoo /@ {neweqL, newicL, newbcL, newslopeL, neweqR, newicR, newbcR, 
          newslopeR}, {{1}, {5}}];

odebcmid = bcmidfunc[odeslopeL, odeslopeR]
odeicmid = icmid;
tend = 40;
sollst = NDSolveValue[{odeL, odeR, odeicL, odeicR, odebcL, odebcR, 
        odebcmid, odeicmid}, {T1 /@ grid, T2 /@ grid, S}, {t, 0, tend}];

solL = rebuild[#, grid, -1]&@sollst[[1]];
solR = rebuild[#, grid, -1]&@sollst[[2]];
solS = sollst[[-1]];
sol = {x, t} \[Function] 
      Piecewise[{{solL[x/solS@t, t], ξL <= x/solS@t <= ξR}}, 
        solR[(x - solS[t])/(xR - solS[t]), t]];

Manipulate[Plot[sol[x, t], {x, 0, xR}, PlotRange -> {-10, 20}], {t, 0, tend}]
(* Plot3D[sol[x, t], {x, 0, xR}, {t, 0, tend}, PlotRange -> {-10, 20}, 
   PlotPoints -> 50] *)

enter image description here

Mathematica graphics

$\endgroup$
  • $\begingroup$ Thanks for your efforts. I am unable to run your code without error. $\endgroup$ – zhk Nov 12 '16 at 15:54
  • 1
    $\begingroup$ @mmm Have you already execute the definition of DChange and pdetoode? $\endgroup$ – xzczd Nov 12 '16 at 16:01
  • $\begingroup$ Now it is working. Thanks $\endgroup$ – zhk Nov 12 '16 at 16:23
  • $\begingroup$ Is it possible to plot S vs t? $\endgroup$ – zhk Nov 12 '16 at 16:25
  • 1
    $\begingroup$ @mmm Sure, just run Plot[solS[t], {t, 0, tend}] $\endgroup$ – xzczd Nov 13 '16 at 2:24

Your Answer

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

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