JOURNAL OF

COMPUTATIONAL AND
APPLIED MATHEMATICS

ELSEVIER Journal of Computational and Applied Mathematics 81 (1997) 135-144

The numerical solution of one-phase classical Stefan problem

S. K u t l u a y * , A.R. Bahadir, A. Ozde~
Department of Mathematics, Faculty of Arts and Science, University of in6nii, Malatya 44069, Turkey
Received 22 December 1995

Abstract

In this paper, variable space grid and boundary Immobilisation Techniques based on the explicit finite difference are
applied to the one-phase classical Stefan problem. It is shown that all the results obtained by the two methods are in good
agreement with the exact solution, and exhibit the expected convergence as the mesh size is refined.

Keywords: Stefan problem; Variable space grid; Boundary immobilisation; Finite differences

AMS classification: 65N06

1. Introduction

A large number of problems in various areas of applied science appear as movin9 boundary or
phase change problems. These can arise in heat conduction situations in conjunction with a change
of phase and initial and moving boundary conditions, and need to be solved in a time-dependent
space domain with a moving boundary condition. Since the moving boundary is a function of time
and its location has to be determined as a part of the solution, such problems are inherently
nonlinear. In general, the nonlinearity associated with the moving boundary significantly compli-
cates the analysis of this class of problems. A common example is the problem of meltin9 of ice that
was first treated by Stefan [9] and after whom such problems are widely referred to as Stefan
problems.
In recent years, these problems have motivated considerable research into numerical solution
methods. A wide range of numerical methods applied to Stefan problems has been reported with
selected references in Crank's comprehensive book [2].
In this study, as a simple example of a moving boundary problem, the one-dimensional melting ice
problem is taken into consideration with one phase. We present a brief comparison of a few well-
known standard numerical solution techniques based on explicit finite-difference approximations.

* Corresponding author.

0377-0427/97/$17.00 O 1997 Elsevier Science B.V. All rights reserved

PI1 S 0 3 7 7 - 0 4 2 7 ( 9 7 ) 0 0 0 3 4 - 4
136 S. Kutluay et al./Journal of Computational and Applied Mathematics 81 (1997) 135-144

2. One-phase classical Stefan problem

Here we are interested in the temperature distribution U(x, t) in the liquid region, 0 < x < s(t)
and in the location of the liquid/solid (moving) interface. Sometimes we omit the arguments x and t.
A typical dimensionless model of the one-phase Stefan problem is in the following form (see, e.g.,
[4]):
The temperature U is governed by the heat conduction equation
(~U •2U
O<x<s(t), t>O, (1)
& cgx2 '
subject to the boundary conditions

~U
--= -exp(t), x=O, t>0,
~x
(2)
U=0, x=s(t), t>0,
The location of the liquid/solid interface is given by the heat balance equation known as Stefan
condition
ds _ 0U ]
dt ¢?x I' x = s(t), t > 0 , (3)
U=0

Initially (at t -- 0) there is no liquid region which implies the condition

s(0)=0, t--0. (4)

Referring to [3, 4], this problem has the following exact solution for the liquid temperature
distribution U(x, t) and the interface location s(t), respectively:

U(x,t)=exp(t-x)-l, O<.x<.s(t), 0 < t < l , (5)

s(t) = t, t >- O.
The exact solution (5) can be used to compare the numerical solutions and can also be used to
initialise the numerical schemes.
Here we deal with the finite-difference solution of the dimensionless model problem given by
Eqs. (1)-(4). Before proceeding further with this numerical method, it is worth tracing or fixing the
time-dependent moving boundary since its location is unknown a priori. To do this, several
numerical techniques based on finite-difference approximations have been used extensively for the
numerical solution of the Stefan problems. A survey of these can be found in [2, 7].
We are only concerned with variable space grid and boundary immobilisation techniques.
These are reviewed in detail in the subsequent sections and applied to the problem given by Eqs.
(1)-(4).
 S. Kutluay et al./Journal of Computational and Applied Mathematics 81 (1997) 135-144 137

3. Solution of the model problem

3.1. A variable space grid (VSG) method f o r the problem

The number of space intervals between a fixed boundary x = 0 and a moving boundary x = s(t)
was kept constant and equal to N so that the moving boundary always lay on the N t h grid [6]. By
tracking particular grid lines, as opposed to constant x, and differentiating with respect to time
t the following expression was obtained for ith grid point,

0U O_~xUt d~t ~ OU
?7 = +-h-'
and was also assumed that the node xi is moved by the expression

dxi xi ds
dt s(t) d t '
in which the suffices t, i, and x are to be kept constant during the differentiation processes and
omitted for clarity below. Thus, in the dimensionless model problem, the heat conduction equation
(1) takes the form

c~U xi ds cgU__ --c32U O < x < s(t), t > 0 , (6)

&- s dt c?x + ~x 2 '
subject to the conditions (2). Eq. (3), subject to condition (4), remains unchanged.
The key to the method is to note that the grid size Ax = s(t)/N varies with time t in each time step
At, since there is a fixed number N of grid points.
An explicit numerical solution based on finite differences to the above problem is easily obtained
by replacing the time and temperature derivatives at nodes (x'f, tin) by the usual forward and
central differences, respectively. The first derivatives with respect to the time variable, c3U / & and
ds/dt, are represented by the usual forward differences (see, e.g., [8, Section 2]). A suitable
replacement for the temperature gradient at the moving interface (x = s(t) = N A x ) , examined in
[3], is given by the three-term backward difference

OU x 3UN - 4UN-1 + UN- 2

0x =~ - 2Ax + O(Ax2). (7)

An explicit finite difference approximation to Eq. (6) is

kx m gm
u'f +1 = u'r +--2-~Sm (U'P+I -- Up-1) + r(u'f+l -- 2U m + Up-1) (8)

with a truncation error of O ( k ) + O(h2). In the above equations uT' '-" U(xT', t,,), xT' = ihm,
tm = to + ink. hm( =- Ax") is the grid size at the ruth time step, k( -- At) is the time step, to is a starting
time, and r = k/h 2.
The boundary condition at x = 0, in terms of central differences, can be written as
um
- 1 = uT + 2hexp(tm). (9)
138 S. Kutluay et al./Journal of Computational and Applied Mathematics 81 (1997) 135-144

Eq. (8) also includes the fictitious temperature u m

_ 1 at i = 0, namely,
kx"d ~m ,
u '~+1 = u '~ + 2-~S-S~Smt Ul -- Urn-a) + r (u "; -- 2 U'~ + Urn-,). (10)

By eliminating u m
_ a between Eqs. (9) and (10) the model problem is expressed by the following
explicit forms:
Heatflow in 0 < x < s(t):
. ..
u7,+1 = ( 1 - - 2 r ) u ~ ' + 2 r u m + l + 2rh -~ )exp(t,,), i=O, m=O, 1,2,...,

kxT' gm ,,
uT'+l=uT+ 2--~s (Ui+x--UT_,)+r(uT'+I--2um+u']'_I), i=l,...,N--1, m=0,1,2,...,

u}"=0, i=N, m=0,1,2,....

Heat balance at x = s(t):
k
s,,+l = s ~ - ~ - ~ ( 3 u ~ - 4 u ~ _ l +u~_2), m=0,1,2,... (11)

S0~0.

The updated interface location Sin+1 and grid size hm+ ~ are calculated at each time from Eq. (11)
and hm + ~ = Sm+,/N, respectively.
F o r stability it is convenient to use Von N e u m a n n ' s m e t h o d of analysis (see, e.g., [8, Section 2])
which gives the following b o u n d on the size of the time step k:
2h 2
k~< (12)
4 + (hg) 2

in which g = ds/dt.

3.2. A bounda~ immobilisation (BIM) technique for the problem

It is also possible to fix the moving interface by using a fixed co-ordinate system in space for the
moving b o u n d a r y problem. This is accomplished by using the spatial co-ordinate transformation
X
= (13)
s(t)

which was proposed in [5] and first applied to a finite-difference scheme in [1]. U n d e r this
transformation the moving interface x = s(t) is fixed at ¢ = 1 in the new co-ordinate system (~, t).
We now re-formulate the dimensionless model problem defined by Eqs. (1)-(4) in the co-ordinate
system (4, t). F r o m the usual partial differential relations we have
~U 10U OZU 1 ~2U
~x-s~' ~x 2 - s 2 ~ 2 ,
 S. Kutluay et aL /Journal of Computational and Applied Mathematics 81 (1997) 135-144 139

and

OU x OU ~ 3__~t ~ x as a u ~-~t~
-& - - ~ -& + - s 2 dt O ~ + '

where the suffices x and t are to be kept constant during the differentiation processes and omitted
for simplicity below.
Thus, we obtain the following description of the dimensionless model problem under the
transformation (13):
H e a t f l o w in 0 < ~ < 1:
t?U ~ ds OU 1 (~2U
- -

0t s d t t 3 ~~+ s 2 ~ 2 , ~ 0<~<1, t>O, (14)

0U
-- - sexp(t), ~ = O,
(15)
U(1, t)=O, 4=1.
Heat balance at ~ = 1:
ds 10U
4=1, t>0, (16)
dt- s ~3~'
s(0)=0, t=0. (17)
An explicit numerical solution to the above 'immobilised' problem is easily obtained using usual
finite differences. As done before the finite-difference replacements of 02U/~3~ 2 and aU/O¢ at the
nodes (¢~, tm) is respectively replaced by the usual central differences. The first derivatives with
respect to the time variable, OU/& and ds/dt, are replaced by the usual forward difference. An
appropriate replacement for the temperature gradient at the moving interface s(t) is given by
Eq. (7).
One explicit finite difference representation of Eq. (14) is
rh~i Sm r
u7'+1 = u m + ~ (uT'+x - uT'- 1) + ZYs,(u?+l
, - 2u7' + uP-l) (18)

with a truncation error of O ( k ) + O(h2). In the above equations uT' ~ U(~, tr,), ~ = ih,
t,, = to + ink, to is the time at which the numerical process is initialised, s,. ..~ s(t,.), g = ds/dt,
r = k/h 2, h(-- A~) and k ( - At) are the mesh size and time step, respectively. N is the number of
space intervals (in the direction of 4).
Using central difference for the left boundary condition we obtain
vim-1 = l,ln~ "31- 2hs exp(tm). (19)
The finite difference Eq. (18) at ( = 0 becomes
rh~o gin, m r
/'/3 +1 = US "[- ~ t/'/1 -- /Am-1) + ~mm(U~ -- 2u~' + u"- x). (20)
140 S. Kutluay et al. /Journal of Computational and Applied Mathematics 81 (1997) 135-144

Elimination of the fictitious value u ~- 1 between Eqs. (19) and (20) yields

( 2r'~ m 2r (2hr k~oSm)eXp(tm).

+ u'i' +,sin --
Introducing these replacements into the model problem given by Eqs. (14)-(17) and using the
change of variable z(t) = sE(t) for the sake of algebraic simplicity leads to the following finite-
difference form:

u7' + 1 = ( 1 - 2 rZm/lu~~ " + 2Zm

ru~'+l+ 4hr-k~i£m
2,/~z~
exp(tm), i 0,
= m=
0,1,2,
... ,

rh~izm r
U m+l = U m + ""T'~ (Urn+ 1 -- UP-l) "1- - - (Urn+ I - - 2 U 7 ' + uT'- 1),
,~z,, Zm
m=0,1,2,..., i=1,..., N-l,
u~'=0, i=N; m=0,1,2,...,

Zm+I =Zm--~(3U~--4U~-x +U~v-2), m = 0 , 1 , 2 , . . . ,

ZO=0.

in which Zm ~" z(t,,), ~ = dz/dt.

For stability analysis we can again use Von-Neumann's approach to obtain the bound on the
size of the time step k. It can be obtained as
k <<.½hZz. (21)

4. Numerical results

All calculations were performed in double precision arithmetic on a PC-486 Intel processor
running under DOS 5.0 and using Prospero F O R T R A N compiler.
In both methods considered in this study, the numerical process was initialised at to = 0.1 using
the exact solution given by Eq. (5).

4.1. Results for the VSG algorithm

Tables 1 and 2 display, respectively, the numerical results for the temperature distribution and
the interface movement at a final time of tf = 0.5. It is observed that all the results are in good
agreement with the exact solution, and exhibit the expected convergence as the mesh size is refined.
The interface predictions shown in Table 3 are also in good agreement with the exact interface
location at some internal points between the initial time, to = 0.1 and the final time, tf -- 0.5. For
example, at t f = 0 . 5 the percentage error decreases from 0 . 5 2 9 2 2 x 1 0 - 1 % ( N = 1 0 ) to
0.85 × 1 0 - 3 % (N = 80).
There is obviously good agreement between numerical and exact results. Table 4 presents the
values obtained by applying Richardson's extrapolation to the value of the weighted 1-norm
S. Kutluay et al./Journal of Computational and Applied Mathematics 81 (1997) 135-144 141

Table 1
Variable space grid: Values of the temperature distribution as predicted by the
numerical (explicit finite-difference) and exact solutions at a final time of tf = 0.5

x/s Numerical solution Exact

solution
N = 10 N = 20 N = 40 N = 80

0.0 0.650614 0.649198 0.648841 0.648751 0.648721

0.1 0.569972 0.568728 0.568416 0.568338 0.568312
0.2 0.493253 0.492182 0.491914 0.491847 0.491825
0.3 0.420281 0.419372 0.419144 0.419087 0.419068
0.4 0.350873 0.350113 0.349922 0.349875 0.349859
0.5 0.284853 0.284232 0.284077 0.284038 0.284025
0.6 0.222052 0.221565 0.221443 0.221413 0.221403
0.7 0.162314 0.161954 0.161864 0.161842 0.161834
0.8 0.105487 0.105250 0.105191 0.105176 0.105171
0.9 0.051428 0.051310 0.051281 0.051274 0.051271
1.0 0.0 0.0 0.0 0.0 0.0

Table 2
Variable space grid: Values of the location and speed of the moving interface as predicted
by the numerical (explicit finite-difference) and exact solutions at a final time of te = 0.5

N Numerical solution

Sm % Error S,, % Error

10 0.50026461 0.52922 × 10 -1 0.99986499 0.13501 x 10 -1

20 0.50006726 0.13452 × 10 -1 0.99997042 0.29580 x 10 -2
40 0.50001693 0.33860 x 10 -2 0.99999304 0.69600 x 10 -3
80 0.50000425 0.85000 × 10 -3 0.99999831 0.16900 x 10 -3
Exact solution 0.5 Exact solution 1.0

Table 3
Location of the moving boundary

t,, Numerical value, S,, Exact value

s(t,.)
N = 10 N = 20 N = 40 N = 80

0.1 . . . . 0.1
0.2 0.20005774 0.20001448 0.20000362 0.20000091 0.2
0.3 0.30022646 0.30005682 0.30001422 0.30000356 0.3
0.4 0.40039488 0.40009928 0.40002487 0.40000622 0.4
0.5 0.50026461 0.50006726 0.50001693 0.50000425 0.5
142 s. Kutluay et al./Journal of Computational and Applied Mathematics 81 (1997) 135-144

defined b y

Ile[ll = ~ 1 Ni~0
-11 g (xi, tm) ] e = [eo ... e N - 1 ] T

which gives an a p p r o x i m a t e rate o f c o n v e r g e n c e of 1.9947. This agrees with the theoretical

e x p e c t a t i o n o f O(h2).

4.2. Results for the BIM algorithm

Results o b t a i n e d b y the B I M m e t h o d are d i s p l a y e d in T a b l e s 5 a n d 6 with a c o m p a r i s o n of the

exact solution at a final time o f tf = 0.5. It is clearly seen t h a t all the results are r e a s o n a b l y in g o o d
a g r e e m e n t with the exact solution. It is o b s e r v e d t h a t the n u m e r i c a l solution c o n v e r g e s to the exact
o n e as the n u m b e r of intervals is increased.
At s o m e internal points b e t w e e n the initial time to = 0.1 a n d the final time, tf = 0.5 the
e s t i m a t i o n o f the interface l o c a t i o n s h o w n in T a b l e 7 is also in very g o o d a g r e e m e n t with the exact
one, e.g., at t f = 0 . 5 the p e r c e n t a g e e r r o r decreasing f r o m 0.16624 x 1 0 - 1 % ( N = 10) to
0.25960 x 1 0 - 3 % (N = 80).

Table 4
Values of Ilelll for the numerical
predictions shown in Table 1

h Ilelll

0.1 0.002942
0.05 0.000738
0.025 0.000185

Table 5
A boundary immobilisation: Values of the temperature distribution as predicted by the
numerical (explicit finite difference) and exact solutions at tf = 0.5

Numerical solution Exact

solution
N = 10 N = 20 N = 40 N = 80

0.0 0.649102 0.648820 0.648746 0.648728 0.648721

0.1 0.568657 0.568399 0.568334 0.568318 0.568312
0.2 0.492116 0.491898 0.491843 0.491829 0.491825
0.3 0.419311 0.419129 0.419083 0.419071 0.419068
0.4 0.350059 0.349909 0.349871 0.349862 0.349859
0.5 0.284186 0.284066 0.284035 0.284028 0.284025
0.6 0.221527 0.221434 0.221411 0.221405 0.221403
0.7 0.161925 0.161857 0.161840 0.161836 0.161834
0.8 0.105230 0.105186 0.105175 0.105172 0.105171
0.9 0.051300 0.051278 0.051273 0.051272 0.051271
1.0 0.0 0.0 0.0 0.0 0.0
 S, Kutluay et al./Journal of Computational and Applied Mathematics 81 (1997) 135-144 143

Table 6
A boundary immobilisation: Values of the location of the moving interface as predicted by
the numerical (explicit finite-difference) and exact solutions at a final time of location of
tf = 0.5

N Numerical solution

Sm % Error S,. % Error

10 0.49916879 0.16624 0.99976347 0.23653 x 10 -1

20 0.49979212 0.41576 x 10 -1 0.99994376 0.56240 x 10 -2
40 0.49994807 0.10386 x 10 -1 0.99998630 0.13700 x 10 -2
80 0.49998702 0.25960 x 10 -2 0.99999662 0.33800 x 10 -3
E x a c t s o l u t i o n 0.5 Exact solution 1.0

Table 7
Location of the moving boundary

tm Numerical value, S,, Exact value

s(t,,)
N= 10 N=20 N=40 N=80

0.1 . . . . 0.1
0.2 0.19985926 0.19996478 0.19999119 0.19999780 0.2
0.3 0.29979494 0.29994872 0.29998718 0.29999680 0.3
0.4 0.39965955 0.39991493 0.39997875 0.39999469 0.4
0.5 0.49916879 0.49979212 0.49994807 0.49998702 0.5

Table 8
Values of Ilelll for the numerical
predictions shown in Table 5

h Ilelll
0.1 0.000575
0.05 0.000144
0.025 0.000036

It is clearly observed that all the numerical predictions are in good agreement with the analytic
solution and there is certainly evidence of a convergent scheme. Table 8 shows the values of II e II 1
for the predicted temperature distributions displayed in Table 5. In fact, the value of Ilelll
converges to the exact value at a rate of about 1.9966. This is again consistent with the theoretical
expectation of O(h2).
Summarising all the numerical results obtained by using the schemes described above, based on
the explicit finite difference, show reasonably good agreement with the exact solution, and exhibit
the expected convergence as the number of space intervals is increased. The BIM, with respect to
144 S. Kutluay et al./Journal of Computational and Applied Mathematics 81 (1997) 135-144

the temperature distribution, generates slightly more accurate solutions than the VSG. However,
with respect to the interface movement (location and speed), the VSG produces very accurate
solutions than the BIM.
In conclusion, when the two methods described above are applied to the Stefan type problems,
the predicted solutions are in reasonably good agreement with the exact ones for many practical
purposes.

References

[1] J. Crank, Two methods for the numerical solution of moving-boundary problems in diffusion and heat flow, Quart.
J. Mech. Appl. Math. X (2) (1957) 220-231.
[2] J. Crank, Free and Moving Boundary Problems, Clarendon Press, Oxford, 1984.
[3] R.M. Furzeland, A comparative study of numerical methods for moving boundary problems, J. Inst. Maths. Appl. 26
(1980) 411-429.
[4] K.H. Hoffman, Fachbereich Mathematik, Vol. I-III, Berlin Freie Universitat, 1977.
I-5] H.G. Landau, Heat conduction in a melting solid, Quart. J. Appl. Math. 8 (1950) 81-94.
[6] W.D. Murray, F. Landis, Numerical and machine solutions of the transient heat conduction problems involving
melting or freezing, J. Heat Transfer 81(C) (1959) 106-112.
I-7] J.R. Ockendon, W.R. Hogkins, Moving Boundary problems in Heat and Diffusion, 2nd ed., Clarendon Press,
Oxford, 1975.
I-8] G.D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, 3rd ed., Clarendon
Press, Oxford, 1987.
[91 J. Stefan, Uber die theorie der eisbildung inbesondee uber die eisbindung im polarmeere, Ann. Phys. U. Chem. 42
(1891) 269-286.