# How can I use Mathematica to solve a Stefan Problem using an explicit scheme?

I would like to plot the temperature distribution for a particular case. the problem is stated in the paper, behind a paywall here, summarized as

Consider a one-dimensional container of length l, full of liquid with a freezing temperature $T_m$. Suppose the initial temperature of the liquid $T_L$ is higher than $T_m$, and one end of the liquid x = 0 is maintained at temperature $T_S$ (< $T_m$) for t > 0, whereas the other end $x = l$ is insulated. The solidification process consequently starts from $x = 0$, and extends over increasing intervals as the time $t$ increases (a well-known Stefan problem). We assume that the material density $\rho$ is constant; and the thermophysical properties are the latent heat $L$, the respective specific heats of the liquid and solid $c_L$ and $c_S$, and the respective thermal conductivities of the liquid and solid $k_L$ and $k_S$.

The initial and boundary conditions are given by,

The temperature distribution follows

where $E(x,t)$ is the enthalpy at time $t$ and position $x$. The enthalpy at time zero is $$E(x,0) = \rho c_L (T_L-T_m) + \rho L$$

and for later times is given by

Taking $T_L=37^\circ \mathrm{C}$, $T_S=-200^\circ \mathrm{C}$, and $l=0.1\mathrm{m}$, and the other constants defined below, how can I calculate the temperature distribution using Mathematica?

Where is the problem? If I'd like to have a temperature (T) distribution at t=1000s, then i need data at t=999s. Which makes it quite awkward to write code in Mathematica. Can anyone help me on this?

This is the work I have done so far:

Clear["Global*"];
T[i_, n_] :=
If[En[i, n] <= 0,
Tm + En[i, n]/(ρ cs),
If[0 < En[i, n] < ρ L, Tm,
If[ρ L <= E[i, n],
Tm + (En[i, n] - ρ L)/(ρ cl)]]]
En[i_, n_ + 1] :=
En[i, n + 1] =
En[i, n] + Δt/Δx*(-ks (
T[i, n] - T[i - 1, n])/Δx -
kl (T[i, n] - T[i + 1, n])/Δx)
T[0, n_] := T1[t] /. t -> n Δt;
M = 5;
T[M, n_] := T2[t] /. t -> n Δt;
T[i_, 0] := T0[x] /. x -> i Δx
En[i_, 0] := ρ cl (T1[0] - Tm) + ρ L;

L = 0.1; ks = 0.00266; kl = 0.006; \
cs = 1.7; cl = 4.1868; ρ = 1000;
Δt = 0.1; Δx = L/M; Tm = 0;

T0[x_] = -200;
T1[t_] = 37;
T2[t_] = -200;

Table[ListPlot[Table[{Δx i, T[i, n]}, {i, 0, M}],
PlotJoined -> True, PlotRange -> {0, 0.4}, AxesLabel -> {"x", ""},
PlotLabel -> "T[x,t], t=" <> ToString[Δt n]], {n, 0,
80, 20}]


Would it be possible to model both types of phase transitions - not just freezing but also melting?

• I need to know a little bit more. So the material is length L - do you have it divided into subsections? If so, how many? What are E and rho? Both E and T have two different forms - a continuous form, $E(x,t)$ and $T(x,t)$, and a form with sub and superscripts - $E^n_i$ and $T^n_i$. How are they related? You also have $T_m$. Be more explicit, or provide a link to where the problem is laid out in more detail. Commented Sep 24, 2015 at 9:05
• @JasonB I edited the question. Thanks in advanced for any feedback! Commented Sep 24, 2015 at 9:20
• do you have some good reason not to use NDSolve ? Commented Sep 24, 2015 at 12:09
• just march it out, Nest[ <code to update E,T> , <initial E> , num time steps ]. (This is the point of an explicit scheme you don't need to solve simultaneously for all time at once ) Commented Sep 24, 2015 at 14:07
• Commented Oct 30, 2015 at 2:27

Here is a function that can return the temperature dynamics for a one-dimensional fluid for a given amount of time. The entire fluid is initially held at temperature T=Tinitial, one end is insulated and stays at this temperature for all time. The other end is in contact with an infinite reservoir at temperature T=Tsource, which can be lower or higher than Tinitial. The accuracy depends on the timestep, Δt and the number of gridpoints.

temperaturelist[Tinitial_, Tsource_, Δt_, tfinal_, ngridpoints_] := Module[{
L = 333.73,
ks = 0.00266,
kl = 0.0006,
cs = 1.7,
cl = 4.1868,
ρ = 1000,
Tm = 0,
tmax, ntpoints, ntdatapoints, Δx,
temperature, tempfunction, templist,
enthalpy, Δenthfunc
},

Δx = 0.1/ngridpoints;
tmax = tfinal*60;
ntpoints = Round[tmax/Δt];
ntdatapoints = 300;

(*Constants defined above, next define the temperature as an instantaneous local function of the enthalpy*)
tempfunction[enth_] := Tm + Which[enth <= 0,
enth/(ρ cs),
0 < enth < ρ L,
0,
enth >= ρ L,
(enth - ρ L)/(ρ cl)
];

(*Next define the change in enthalpy, which depends on the instantaneous temperature for the grid point in question and
also on the nearest neighbor grid points*)

Δenthfunc[temp_] := Module[{klist, qminus, qplus},
klist = If[# < Tm, ks, kl] & /@ temp;
qminus = -((temp - RotateRight[temp])/(.5 Δx (1/klist +1/RotateRight[klist])))[[2 ;; -2]];
qplus = -((RotateLeft[temp] - temp)/(.5 Δx (1/klist +1/RotateLeft[klist])))[[2 ;; -2]];
Δt/Δx (qminus - qplus)
];
(* Define the initial temperature and enthalpy *)
temperature = ConstantArray[Tinitial, ngridpoints];
temperature[[1]] = Tsource;
enthalpy =
ConstantArray[
If[Tinitial > Tm, ρ cl (Tinitial - Tm) + ρ L, ρ cs (Tinitial -Tm)], ngridpoints - 2];

(* Next run the temperature dynamics simulation.  The value of Δt is critical here, or else the results are nonsense *)
PrintTemporary[
"Running dynamics, number of timepoints = " <>
IntegerString[ntpoints]];
Monitor[
templist = Reap[
Do[
temperature[[2 ;; -2]] = tempfunction /@ enthalpy;
enthalpy += Δenthfunc[temperature];
If[Mod[n - 1, Round[ntpoints/ntdatapoints]] == 0,
Sow[temperature]];
, {n, ntpoints + 1}]][[2, 1]];
, n];
(* Return the result *)
templist
];


tlist1 = temperaturelist[37., -200., 0.1, 50, 128];
tlist2 = temperaturelist[-37., 200., 0.1, 50, 128];


and here are plots of the results, with the red line marking the boundary between liquid and solid phases.

Grid[{(Show[
ListDensityPlot[#, PlotLegends -> Automatic,
DataRange -> {{0, .1}, {0, 50}},
FrameLabel -> {"x (m)", "t (min)"}, ImageSize -> 400],
ListContourPlot[#, ContourShading -> None, Contours -> {0},
DataRange -> {{0, .1}, {0, 50}}, ContourStyle -> {{Thick, Red}}]
] & /@ {tlist1, tlist2})}]


• Great work! Many thanks! This is very useful and other approach to the problem. The ListDensityPlot is not working properly, due to one ']' too many, but I'll figure it out eventually. Thank you! Commented Sep 29, 2015 at 17:25
• Corrected . Show[ was added Commented Sep 29, 2015 at 17:31
• Do you perhaps know why it doesn't work for the opposite case when defining initial conditions, e.g. Tl = 10+ 273.15; Ts = -10 + 273.15; and initial temperature is -10+273` - Where is that defined in your code? Commented Sep 29, 2015 at 20:42
• Now you are changing the problem after the fact - you want to not only be able to model freezing of a liquid but also melting of a solid. This physical model can do both, but the code I wrote only does freezing. I've just changed it to be more versatile, so I will post that here shortly. Commented Sep 30, 2015 at 10:16
• But I think that you need more experience with programming, before tackling such a problem. It doesn't need to be Mathematica, this would be a fun exercise in C++ or python. Become more familiar with lists, with loops, etc. Also, the equations you posted originally were wrong - specifically the enthalpy equations. When you combined them from the original paper, you made some errors that would have doomed you even if your programming were perfect - so I had to go to the source paper and figure it out myself. Commented Sep 30, 2015 at 10:19