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,

enter image description here

The temperature distribution follows

enter image description here

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

enter image description here 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?

enter image description here

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:

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?

  • $\begingroup$ 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. $\endgroup$
    – Jason B.
    Sep 24, 2015 at 9:05
  • $\begingroup$ @JasonB I edited the question. Thanks in advanced for any feedback! $\endgroup$
    – energyMax
    Sep 24, 2015 at 9:20
  • 2
    $\begingroup$ do you have some good reason not to use NDSolve ? $\endgroup$
    – george2079
    Sep 24, 2015 at 12:09
  • 1
    $\begingroup$ 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 ) $\endgroup$
    – george2079
    Sep 24, 2015 at 14:07
  • 1
    $\begingroup$ Related:mathematica.stackexchange.com/a/58622/1871 $\endgroup$
    – xzczd
    Oct 30, 2015 at 2:27

1 Answer 1


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,
      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 = 
     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 *)
    "Running dynamics, number of timepoints = " <> 
    templist = Reap[
         temperature[[2 ;; -2]] = tempfunction /@ enthalpy;
         enthalpy += Δenthfunc[temperature];
         If[Mod[n - 1, Round[ntpoints/ntdatapoints]] == 0, 
         , {n, ntpoints + 1}]][[2, 1]];
    , n];
   (* Return the result *)

Here are two examples, one for freezing, the other for melting:

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.

      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})}]

enter image description here

  • $\begingroup$ 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! $\endgroup$
    – energyMax
    Sep 29, 2015 at 17:25
  • $\begingroup$ Corrected . Show[ was added $\endgroup$
    – energyMax
    Sep 29, 2015 at 17:31
  • $\begingroup$ 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? $\endgroup$
    – energyMax
    Sep 29, 2015 at 20:42
  • 1
    $\begingroup$ 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. $\endgroup$
    – Jason B.
    Sep 30, 2015 at 10:16
  • 1
    $\begingroup$ 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. $\endgroup$
    – Jason B.
    Sep 30, 2015 at 10:19

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.