Solving the Stefan Problem with WhenEvent

Question

The formulation of the problem:

I tried to solve it with MOL and the method of V.R. Voller:

The script:

Clear["Global`*"]
n = 100; (*Number of Elements*)
H = 0.05; (*Length*)
dH = H/n; (*Length of an Element*)

\[Rho]s = 917;
\[Rho]l = 1/(1.003*^-3);
cs = 0.185*^3;
cl = 4.179*^3;
ks = 1.16*1.91;
kl = 613*^-3;

\[Alpha]l = kl/(\[Rho]l*cl);
\[Alpha]s = ks/(\[Rho]s*cs);

L = 333;

q = 40000;
\[Theta]i = 271;
\[Theta]m = 273;
Z = Table[Subscript[z, i], {i, 0, n}];
For[i = 0, i <= n, i++, Subscript[z, i] = i*dH];

tmax = 10;

(*Defining Tables *)
\[CapitalTheta][t_] = Table[Subscript[\[Theta], i][t], {i, 0, n}];
\[CapitalAlpha][t_] = Table[Subscript[\[Alpha], i][t], {i, 0, n}];
F[t_] = Table[HeavisideTheta[Subscript[\[Theta], i][t]], {i, 0, n}];

dd\[CapitalTheta] = 
  NDSolve`FiniteDifferenceDerivative[2, 
   dH*Range[0, n], \[CapitalTheta][t]];
d\[CapitalTheta] = 
  NDSolve`FiniteDifferenceDerivative[1, 
   dH*Range[0, n], \[CapitalTheta][t]];

dvbls = Head[#] & /@ \[CapitalAlpha][t];

(*eqn*)
\[Theta]eqn = 
  Thread[D[\[CapitalTheta][t], t] - \[CapitalAlpha][t]*
      dd\[CapitalTheta] + \[Rho]l*L*D[F[t], t] == 0];

(*BC*)
bc0 = d\[CapitalTheta][[1]] == -\[Rho]l*cl*q/kl;
bcH = Subscript [\[Theta], n][t] == \[Rho]s*cs*(\[Theta]i - \[Theta]m);

\[Theta]eqn[[1]] = bc0;
\[Theta]eqn[[-1]] = bcH;

(*IC*)
\[Theta]ic = 
  Thread[\[CapitalTheta][0] == cs*\[Rho]s*(\[Theta]i - \[Theta]m)];
\[Alpha]ic = Thread[\[CapitalAlpha][0] == \[Alpha]s];

(*Events*)
event1 = WhenEvent[Subscript[\[Theta], #][t] <= 0, 
     Subscript[\[Alpha], #][t] -> \[Alpha]s] & /@ Range[0, n];
event2 = WhenEvent[Subscript[\[Theta], #][t] > 0, 
     Subscript[\[Alpha], #][t] -> \[Alpha]l] & /@ Range[0, n];

system = Join[\[Theta]eqn, \[Theta]ic, \[Alpha]ic, event1, event2];

lines = NDSolveValue[system, 
  Join[\[CapitalTheta][t], \[CapitalAlpha][t]], {t, 0, tmax}, 
  DiscreteVariables -> dvbls, 
  StepMonitor :> (laststep = thisstep; thisstep = t; 
    stepsize = thisstep - laststep;)]

If I try to solve it without WhenEvent - assuming constant diffusing coefficient, I get excellent results. But if I include the events, I get a lot of errors:

I didn't find any good documentation to those errors. It would be very helpful if someone had a reference to that kind of error, or knows where my code went wrong.

Thanks!

Answer 1

With a little modification this code works also with the message

NDSolveValue::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended.

Nevertheless we have output and it looks reasonable for alpha

Clear["Global`*"]
n = 50;(*Number of Elements*)H = 0.05;(*Length*)dH = 
 1/n;(*Length of an Element*)\[Rho]s = 917;
\[Rho]l = 1/(1.003*^-3);
cs = 0.185*^3;
cl = 4.179*^3;
ks = 1.16*1.91;
kl = 613*^-3;

\[Alpha]l = kl/(\[Rho]l*cl);
\[Alpha]s = ks/(\[Rho]s*cs);

L = 333;

q = 40000;
\[Theta]i = 271;
\[Theta]m = 273;
Z = Table[Subscript[z, i], {i, 0, n}];
For[i = 0, i <= n, i++, Subscript[z, i] = i*dH];

tmax = 10;

(*Defining Tables*)
\[CapitalTheta][t_] = Table[Subscript[\[Theta], i][t], {i, 0, n}];
\[CapitalAlpha][t_] = Table[Subscript[\[Alpha], i][t], {i, 0, n}];
F[t_] = Table[
   ArcTan[10^2 Subscript[\[Theta], i][t]]/Pi + 1/2, {i, 0, n}];

dd\[CapitalTheta] = 
  NDSolve`FiniteDifferenceDerivative[2, 
   dH*Range[0, n], \[CapitalTheta][t]];
d\[CapitalTheta] = 
  NDSolve`FiniteDifferenceDerivative[1, 
   dH*Range[0, n], \[CapitalTheta][t]];

dvbls = Head[#] & /@ \[CapitalAlpha][t];

(*eqn*)
\[Theta]eqn = 
  Thread[D[\[CapitalTheta][t], t] - \[CapitalAlpha][t]*
      dd\[CapitalTheta]/H^2 + \[Rho]l*L*D[F[t], t] == 0];

(*BC*)
bc0 = d\[CapitalTheta][[1]] == -H \[Rho]l*cl*q/kl;
bcH = Subscript[\[Theta], n][t] == \[Rho]s*cs*(\[Theta]i - \[Theta]m);

\[Theta]eqn[[1]] = bc0;
\[Theta]eqn[[-1]] = bcH;

(*IC*)
\[Theta]ic = 
  Thread[\[CapitalTheta][0] == cs*\[Rho]s*(\[Theta]i - \[Theta]m)];
\[Alpha]ic = Thread[\[CapitalAlpha][0] == \[Alpha]s];

(*Events*)

event1 = WhenEvent[Subscript[\[Theta], #][t] <= 0, 
     Subscript[\[Alpha], #][t] -> \[Alpha]s] & /@ Range[0, n];
event2 = WhenEvent[Subscript[\[Theta], #][t] > 0, 
     Subscript[\[Alpha], #][t] -> \[Alpha]l] & /@ Range[0, n];


system = Join[\[Theta]eqn, \[Theta]ic, \[Alpha]ic, event1, event2];


lines = NDSolveValue[system, 
   Join[\[CapitalTheta][t], \[CapitalAlpha][t]], {t, 0, tmax}, 
   DiscreteVariables -> dvbls];

Note that we map solution on the unit interval and use continues function for F. Visualization

{ListLinePlot[
  Table[Drop[lines, -n - 1][[1 ;; 20]]/(\[Rho]s*cs), {t, 1, 10, 1}], 
  PlotRange -> All, PlotLegends -> Automatic], 
 ListLinePlot[Table[Drop[lines, -n - 1]/(\[Rho]s*cs), {t, 1, 10, 1}], 
  PlotRange -> {-10, 5}, Frame -> True, PlotLabel -> T], 
 ListLinePlot[Table[Drop[lines, n + 1], {t, 1, 10, 1}], 
  PlotRange -> All, PlotLabel -> "\[Alpha]"]}