# Solving the Stefan Problem with WhenEvent

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] =
NDSolveFiniteDifferenceDerivative[2,
dH*Range[0, n], \[CapitalTheta][t]];
d\[CapitalTheta] =
NDSolveFiniteDifferenceDerivative[1,
dH*Range[0, n], \[CapitalTheta][t]];

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

(*eqn*)
\[Theta]eqn =
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 =

(*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!

• Did you seen this post mathematica.stackexchange.com/questions/184920/… ? Mar 15, 2021 at 20:44
• @AlexTrounev I had - it's a really good one so thanks. But I also try to understand why my script does not work, as it applies a different method.
– Ymir
Mar 16, 2021 at 6:25
• You code is fine, but it not for numerical computation. You can compare it with Tim Laska approach to compute solidification problem with enthalpy method. Did you seen this paper emerald.com/insight/publication/issn/0961-5539 Mar 16, 2021 at 15:20

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] =
NDSolveFiniteDifferenceDerivative[2,
dH*Range[0, n], \[CapitalTheta][t]];
d\[CapitalTheta] =
NDSolveFiniteDifferenceDerivative[1,
dH*Range[0, n], \[CapitalTheta][t]];

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

(*eqn*)
\[Theta]eqn =
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 =

(*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]"]}