Finite element moving mesh analysis of phase change

problems with natural convection

R.T. Tenchev
a,
*
, J.A. Mackenzie
a
, T.J. Scanlon
b
, M.T. Stickland
b
a
Department of Mathematics, University of Strathclyde, Glasgow G1 1XH, UK
b
Department of Mechanical Engineering, University of Strathclyde, Glasgow G1 1XJ, UK
Accepted 13 March 2005
Available online 3 June 2005
Abstract
This paper discusses the application of an r-renement, moving mesh technique for the solution of heat transfer problems with
natural convection and phase change. The moving mesh technique keeps the number of elements and their connectivity xed and
clusters the nodes towards the phase change front at the expense of the solution of an extra dierential equation. The governing
dierential equations describing the physical problem are modied to account for the mesh movement between time steps. The
energy conservation equation uses the apparent heat capacity method to take into account the latent heat of phase change. The nite
element discretization of all equations is presented. Several test problems are solved and the moving mesh FEM results are in a very
good agreement with those in the published literature. The sensitivity of the results to variations of some user-denable computa-
tional parameters is found to be low, which means that the moving mesh method may be used without extensive previous experience.
Its basic advantage is that less elements may be used to achieve accurate results.
2005 Elsevier Inc. All rights reserved.
Keywords: Moving mesh; Phase change; Natural convection
1. Introduction
Phase change problems (solidication or melting) are
also known as Stefan problems or moving boundary
problems. Their main characteristic feature is that, in
addition to the xed boundaries of the domain, there
is a boundary across which the phase change takes place
(Crank, 1984; Ockendon and Hodgkins, 1975). This
phase change boundary is time dependent and is not
known a priori. Thermodynamic equilibrium conditions
are to be satised on it. The problems are highly non-lin-
ear and analytical solutions do not exist except in some
very simple cases. Due to their importance many numer-
ical approximations have been developed to solve these
problems (Dalhuijsen and Segal, 1981; Lewis and Rob-
erts, 1987; Pardo and Weckman, 1990; Salcuden and
Abdullah, 1988; Thomas et al., 1984), which can be clas-
sied in two main categories: front tracking methods
and xed grid methods.
With front tracking methods (Gupta, 2000; Pardo
and Weckman, 1990; Voller et al., 1990) the discrete
phase change front is tracked continuously and treated
as a moving boundary between the liquid and solid
phase. Dierent dierential equations may be used in
these phases. The latent heat involved in the phase
change is treated explicitly (hence accurately) as an
internal boundary condition. The change in size and
shape of the computational domain requires either grid
movement techniques or co-ordinate system transforma-
tions. The front tracking method is applicable for prob-
lems with isothermal phase change. It is generally not
suitable for problems where the phase change takes
0142-727X/$ - see front matter 2005 Elsevier Inc. All rights reserved.
doi:10.1016/j.ijheatuidow.2005.03.003
*
Corresponding author. Tel.: +44 141 548 3668; fax: +44 141 552
8657.
E-mail address: [email protected] (R.T. Tenchev).
www.elsevier.com/locate/ijh
International Journal of Heat and Fluid Flow 26 (2005) 597612
place within a temperature interval and involves the for-
mation of a so called mushy region.
With xed grid methods (Crivelli and Idelsohn, 1986;
Fachinotti et al., 1999; Morgan et al., 1978; Nedjar,
2002; Rolph and Bathe, 1982; Tamma and Namburu,
1990; Voller and Cross, 1981) the phase change front
is not tracked explicitly but is instead recovered a poste-
riori from the computed temperature eld. These meth-
ods are also known as single domain methods because
the same dierential equation can be used for the solid
and liquid region. Hence, their basic advantage is that
the numerical solution can be achieved through simple
modications of existing heat transfer numerical meth-
ods. The essential feature, however, is the way the latent
heat evolution is treated: either by the use of an appar-
ent heat capacity coecient or by the use a heat source/
sink term. These approaches are collectively known as
the enthalpy method, because they can be derived from
the energy conservation equation written in terms of the
enthalpy, which is the sum of the apparent and latent
heat. Enthalpy methods are the natural choice when
the phase change occurs over a temperature interval.
In the case of isothermal phase change there is a discon-
tinuity in the enthalpy across the phase change front,
which is normally smoothed by assuming that the phase
change occurs over a temperature interval.
When the apparent heat capacity method is used the
heat capacity of the mushy region is increased by a term
which is directly proportional to the latent heat and is
inversely proportional to the phase change interval.
There will be loss of accuracy if a region changes its state
from solid to liquid or vice versa without passing
through the mushy state. The latter may be avoided
by either decreasing the computational time step or
increasing the phase change interval. The argument is
also valid for the source/sink term formulation.
Enthalpy methods either require a phasewise exact
integration strategy (Nigro et al., 2000) or a ne mesh
near the phase change front in order to capture the large
enthalpy gradient in the mushy region. The smaller the
phase change interval the narrower the mushy region
is and the more rened the mesh should be. If it is
expected that the phase change front will transverse
the whole domain then a rened mesh should be used
everywhere, although it is required only close to the
phase change front, and near domain boundaries to pick
up boundary layer eects. This serious disadvantage can
be overcome if adaptive mesh renement is employed.
Either an h or r-adaptive strategy could be used for this
purpose. One advantage of the r-renement or moving
mesh approach is that the element connectivity usually
remains unchanged and hence the code is much easier
to implement than an h-renement code which requires
a more complicated data structure. Examples of such
moving mesh methods include: (Beckett et al., 2001,
2002; Cao et al., 1999; Huang, 1999; Lynch, 1985; Lynch
and ONeill, 1981; Miller and Miller, 1981; Mackenzie
and Robertson, 2000).
The idea behind the moving mesh method is that a
local mesh renement is applied only where and when
it is required. The total number of elements and their
connectivity stay constant throughout the solution
process. A new mesh is generated at every time step.
The computed solution from the previous time step
(e.g. the phase change isotherm) is used to determine
where the mesh should be rened. The method is based
on the single domain formulation although the grid is
not xed. The method is not a front tracking one
because the moving mesh technique does not ensure that
element nodes lie on the phase change isotherm thus
tracking it explicitly.
The present paper investigates the application of the
moving mesh technique in the nite element solution
of two dimensional phase change problems when the
apparent heat capacity method is used and the uid ow
is driven by natural convection. Both melting and freez-
ing problems are considered. The sensitivity of the
results to the variation of some computational para-
meters is studied.
2. Governing equations
The governing dierential equations are presented in
the order they are solved at each time step of the FEM
solution. First, the energy equation is solved (using com-
puted velocities from the previous time step), which
gives the temperature distribution and the position of
the phase change front. Second, the moving mesh equa-
tions are solved to adapt the mesh towards the phase
Nomenclature
G moving mesh monitor function
h latent heat of phase change
T
m
phase change temperature
b coecient of volume expansion
e
1
, e
2
lower and upper limits of the phase change
temperature interval e = e
1
+ e
2
l
1
, l
2
coecients in G, governing mesh density at
the phase change front and at the boundaries
n, g computational coordinates in 2D
f, s computational and physical coordinates in
1D
s dimensionless time
598 R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612
change front. Since the time step and the movement of
the front are small there is no need to iterate between
the two analyses. Finally, the NavierStokes equations
are solved to determine the velocity of the uid due to
the natural convection.
2.1. Energy equation
The energy equation for transient heat transfer by
both conduction and convection is
o(qcT)
ot
= \ (k\T) \ (qcTu) (1)
and it should be modied to account for the mesh
movement.
When the mesh (i.e. the dierential volume for which
Eq. (1) is derived) is moving with velocity dx/dt the total
time derivative is required:
D(qcT)
Dt
=
o(qcT)
ot
\(qcT)
dx
dt
(2)
The notation D()/Dt means dierentiation with respect
to time along the grid trajectory (x(t), y(t)). If the mesh
is xed D()/Dt = o()/ot. Substituting the partial time
derivative from Eq. (1) into Eq. (2) and rearranging,
results:
D(qcT)
Dt
= \ (k\T) \ (qcTu) \(qcT)
dx
dt
(3)
When the uid is incompressible ($ u = 0 and
$ (qcTu) = u $(qcT)), Eq. (3) can be written as:
D(qcT)
Dt
= \ (k\T) u
dx
dt
_ _
\(qcT) (4)
The above equation can be derived directly from Eq. (1)
if the dierential volume (i.e. the mesh), which is moving
with dx/dt, is considered xed in space and the uid
velocity instead of u is u dx/dt.
The phase change is assumed to take place within a
temperature interval. The latent heat h is accounted
for by using an apparent heat capacity:
qc =
qc
solid
if T
m
e
1
>T
qc
solid
qc
liquid
2

qh
e
1
e
2
if T
m
e
1
6T 6T
m
e
2
qc
liquid
if T > T
m
e
2
_

_
(5)
2.1.1. FEM approximation
Galerkins FEM discretization results in a system of
ordinary dierential equations:
C
T
DT
Dt
(K
T
U
T
)T = F
T
(6)
where C
T
=

C
e
T
, K
T
=

K
e
T
, etc. In future, for sim-
plicity, the superscript e will be omitted in all element
matrices. Matrices C
T
and K
T
are the same as in the case
of a stationary mesh:
C
T
=
_ _
X
N
T
qcNdX;
K
T
=
_ _
X
k
oN
T
ox
oN
ox

oN
T
oy
oN
oy
_ _
dX; (7)
and N = N
1
N
2
N
n
| is the row vector of the
shape functions for an element with n nodes and are also
used as weighting functions. In the present studies four
node bi-linear elements are used. For the examples con-
sidered later we have found that the convective term is
not dominant and there is no need to use upwinded
weighting functions. The subscript T is used to indicate
matrices and vectors used in the energy equation, i.e.
when temperature is computed.
The term in Eq. (3) involving the mesh velocity dx/dt
may be treated in two ways. Implicitly, which is the
more accurate formulation if the mesh velocity is known
at the current mesh conguration and time step. This
formulation contributes to U
T
:
U
T
=
_ _
X
N
T
qc u
dx
dt
_ _
oN
ox
v
dy
dt
_ _
oN
oy
_ _
dX;
F
T
= 0 (8a)
Explicitly, it will contribute to F
T
, using the temperature
computed at the previous time step and mesh.
U
T
=
_ _
X
N
T
qc u
oN
ox
v
oN
oy
_ _
dX;
F
T
=
_ _
X
N
T
qc\T
i1

dx
dt
dX (8b)
The mesh velocity can be computed only approximately,
(itime, Dttime step) as:
dx
dt
~
x
i
x
i1
Dt
; (9)
which may be considered as the averaged velocity at
time i
1
2
. So both implicit and explicit formulations
are approximate. In the present study the explicit formu-
lation is used for simplicity. Making use of nodal values
and shape functions, the force term in Eq. (8b) is
F
T
=
_ _
X
N
T
qc
oN
j
ox
i1
T
i1
j
_ _
(N
j
x
i
j
) (N
j
x
i1
j
)
Dt
_

oN
j
oy
i1
T
i1
j
_ _
(N
j
y
i
j
) (N
j
y
i1
j
)
Dt
_
dX (10)
where the superscript i denotes at which time the vari-
able is computed. Repeated index j within the brackets
means summation over the element nodes.
The ordinary dierential equations, Eq. (6), are
discretized by the backwards Euler method (i.e. the h
method, with h = 1). The additional advantage, compared
to other methods with h = 1/2 or h = 3/4, is that the load
R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612 599
vector from the previous time step is not required. This
vector is available, but it is computed for a dierent
mesh and its use will bring about complications.
The time derivative in Eq. (6), for node j, is discret-
ized as DT=Dt ~ (T
i
j
T
i1
j
)=Dt, where T
i
j
and T
i1
j
are
the computed nodal values of the temperature at times
i and i 1, respectively. The position of node j changes
from time i 1 to time i due to the mesh movement, i.e.
the temperature at the node is function of its spatial
position, x, as well as time, t. This is exactly what is im-
plied by the total derivative and its computation comes
naturally in the presented FEM discretization.
2.2. Moving mesh equations
A graphical interpretation of the adaptive mesh gen-
eration is shown for the one dimensional case by the
mapping in Fig. 1. A uniform computational mesh is
dened along f. The real physical mesh is to be gener-
ated along s. Let at s
m
be the location of the phase
change isotherm at a given time. A rened mesh is
required towards point s
m
. If a function f(s) is con-
structed such that it has a high gradient in the vicinity
of s
m
the mapping represented by the dotted lines will
give the required mesh renement.
The monitor function G that governs the mesh rene-
ment is dened as:
G(s) = C
df
ds
(11)
The construction of f(s) is achieved by the principal of
equidistribution (i.e. the areas below G(s) between con-
secutive mesh points to be equal) where G(s) is a known
function. Eq. (11) can be transformed to a boundary
value problem by dierentiation:
d
ds
G
1
df
ds
_ _
= 0 (12)
In 2D the mapping is dened by minimization a func-
tional of the form, (Cao et al., 1999):
U(n; g) =
_
X
(\n
T
G
1
\n \g
T
G
1
\g) dxdy (13)
The EulerLagrange equations characterizing the extre-
mum of Eq. (13) take the form:
\ (G
1
\n) = 0
\ (G
1
\g) = 0 (14)
where n, g are coordinates in the computational 2D do-
main and G is a matrix of monitor functions:
G =
G
x
(x; y) 0
0 G
y
(x; y)
_ _
(15)
To have an explicit denition of the physical coordi-
nates Eqs. (12) and (14) have to be written in terms of
os/of and ox/on and a pseudo transient formulation is
introduced to smooth out the evolution of the mesh
(Beckett et al., 2001; Cao et al., 1999; Huang, 1999). A
diagonal monitor function allows to decoupled the 2D
equations.
os
o
~
t
=
1
P
s
A
s
o
2
s
of
2
D
s
os
of
_ _
(16)
ox
o
~
t
=
1
P
A
o
2
x
on
2
B
o
2
x
onog
C
o
2
x
og
2
D
ox
on
E
ox
og
_ _
(17a)
oy
o
~
t
=
1
P
A
o
2
y
on
2
B
o
2
y
onog
C
o
2
y
og
2
D
oy
on
E
oy
og
_ _
(17b)
where A, B, . . . , E are non-linear functions of the physi-
cal coordinates, P and P
s
are normalizing, scaling coef-
cients, which are evaluated using the last known
physical coordinates (see Appendix A). The steady state
solution is required (i.e. os=o
~
t = 0, etc.) and the transient
formulation (
~
t is a ctitious time) is used only to indi-
cate a solution procedure that would be similar to the
one used to solve the transient dierential equations of
the physical problem. Eq. (16) is to be solved on each
continuous boundary edge whose two end points are
xed and do not move. The computed nodal coordinates
on the boundary edges are to be used as Dirichlet
boundary conditions for Eqs. (17).
The denition of the monitor functions in Eq. (15) is
based on the work of Beckett et al. (2001):
G
x
= 1

N
X
i=1
l
i

l
2
i
x x
m;i
[ [
2
1
_
_

_
_

_l
L
(T)
G
y
= 1

N
Y
i=1
l
i

l
2
i
[x x
m;i
[
2
1
_
_

_
_

_l
L
(T) (18)
s

G (s) =C
s

s
m
s
G
= (s)

m
Fig. 1. Mapping from a uniform computational mesh f to a physical
mesh s via a monitor function G(s). Mesh renement at point s
m
.
600 R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612
The term [x x
m,i
[ is the distance of a point x = (x, y) to
a curve i, which is dened discretely by the set of points
x
m,i
= (x
m,i
, y
m,i
), toward which mesh renement is
required. The total number of such curves is N
x
and
N
y
for mesh movement in the x and y directions, respec-
tively. Coecients l
i
governs the mesh renement to-
wards curve i which is either the phase change front or
the xed boundaries where thermal and velocity layers
are to be resolved. These boundaries are dierent for
the x and y mesh movement equations. For example,
in the melting problem described later in Fig. 2, the
mesh is rened only towards the left boundary in the x
equation (17a). There is no need to rene towards the
right boundary because the material near it is expected
to be in solid state throughout the whole solution.
Renement towards the top and bottom boundaries is
used in the y equation (17b). Alternatively, instead of
including the boundaries into the monitor function,
the computational mesh can be rened at the corre-
sponding boundaries. This, however, will require the
actual generation of a computational mesh. Coecient
l
L
(T), allows the mesh to be rened in the liquid region
at the expense of the solid one: l
L
= 1 if T < T
m
and
l
L
> 1 otherwise.
2.2.1. FEM approximation
Galerkins discretization of Eqs. (17), on an element
level, results in:
_ _
X
N
T
NdX
ox
ot

_ _
X
N
T
1
P
A
o
2
N
on
2
C
o
2
N
og
2
B
o
2
N
onog
_
D
oN
on
E
oN
og
_
dXx = 0 (19)
Note: The computational coordinates n, g should not
be confused with the local element coordinates. How-
ever, when the computational mesh is uniform, i.e. there
are no mesh renements towards the xed boundaries,
the computational domain may be chosen in such a
way, so that the computational element to be identical
to the isoparametric, non-dimensional, nite element.
In this case there is no need to construct a computa-
tional mesh.
Integrating by parts and ignoring the boundary line
integrals, because Dirichlet boundary conditions will
be used on all boundaries, results in:

_ _
X
N
T
A
P
o
2
N
on
2
dX =
_ _
X
o(N
T
A=P)
on
oN
on
dX
=
_ _
X
N
T
o(A=P)
on
oN
on
dX
_ _
X
oN
T
on
A
P
oN
on
dX
(20)

_ _
X
N
T
C
P
o
2
N
og
2
dX =
_ _
X
o(N
T
C=P)
og
oN
og
dX
=
_ _
X
N
T
o(C=P)
og
oN
og
dX
_ _
X
oN
T
og
C
P
oN
og
dX
(21)

_ _
X
N
T
B
P
o
2
N
onog
dX =
1
2
_ _
X
o(N
T
B=P)
on
oN
og
dX
_

_ _
X
o(N
T
B=P)
og
oN
on
dX
_
(22)
The averaging procedure is used in Eq. (22) to elimi-
nate any bias in the integration by parts. The remaining
terms need not be modied. After the nite element spa-
tial discretization, the resulting system of ordinary dif-
ferential equations takes the form:
M
dx
d
~
t
(A B C D E)x = 0; (23)
where the matrices M, A, . . . , E, are given in Appendix
A. To solve (23) the system is rst linearised by evaluat-
ing the matrices A, . . . , E using the mesh and the solution
at the current time step. The resulting linear set of equa-
tions is solved using an ILU-preconditioned BICGstab
routine. As the current mesh is usually a good initial
guess only a nite number (usually around 10) BICG-
stab iterations are performed.
The reason for using only Dirichlet boundary condi-
tions (which can be determined by solving the 1Dmoving
mesh equations on all boundaries) is that the Neumann
boundary conditions corresponding to the boundary
integrals are not known. In the nite element approxima-
tion it is possible to assume symmetry boundary condi-
tions along each boundary curve. In this case the 1D
moving mesh equations need not be solved. However,
in some cases, Eq. (23) may become ill-conditioned and
the moving mesh to be unstable. On the other hand the
1D mesh movement on a boundary edge is not identical
to the mesh movement that would be calculated for this
edge by the 2D formulation. The dierences die away
quickly, but it is possible to present some numerical
problems, if a proper temporal smoothing (i.e. suitable
ctitious time step D
~
t) is not used.
y
x
L
y
dT/dn=0
dT/dn=0
T
0
T
0
L
x
T
w
Fig. 2. Layout of test problems. Melting problems: L
x
= L
y
= 1;
T
w
= 1; T
0
= 0; T
m
= 0. Freezing problem: L
x
= 0.3 m; L
y
= 0.15 m;
T
w
= 10 C; T
0
= 10.2 C; T
m
= 0 C.
R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612 601
2.3. NavierStokes equations
Assuming density to be constant except in the buoy-
ancy term (Boussinesq approximation) and that the only
body force is due to the gravity acceleration in the y
direction, the NavierStokes equations for a stationary
mesh are:
ou
ot
u \u \
l
q
0
\u
_ _
=
1
q
0
op
ox
(24)
ov
ot
u \v \
l
q
0
\v
_ _
=
1
q
0
op
oy
qg
_ _
(25)
ou
ox

ov
oy
= 0 (26)
The incompressibility condition, Eq. (26), may be re-
laxed by using a penalty formulation for the pressure
(Heinrich and Pepper, 1996):
p = p
s
k
ou
ox

ov
oy
_ _
(27)
The derivatives of the hydrostatic pressure p
s
are:
op
s
ox
= 0;
op
s
oy
= q
0
g (28)
Substituting Eqs. (27) and (28) into Eqs. (24) and (25),
gives:
ou
ot
u \u \
T
l
q
0
\u
_ _

k
q
0
o
ox
ou
ox

ov
oy
_ _
= 0
ov
ot
u \v \
T
l
q
0
\v
_ _

k
q
0
o
oy
ou
ox

ov
oy
_ _
=
(q
0
q)g
q
0
(29)
Taking into account the mesh movement in a similar
way as in the energy equation, the nal form of the
equations to be solved by FEM is
Du
Dt
u \u
dx
dt
\u \
T
l
q
0
\u
_ _

k
q
0
o
ox
ou
ox

ov
oy
_ _
= 0
Dv
Dt
u \v
dx
dt
\v \
T
l
q
0
\v
_ _

k
q
0
o
oy
ou
ox

ov
oy
_ _
=
(q
0
q)g
q
0
(30)
The nite element approximation is given in Appendix B.
3. Numerical examples
A schematic presentation of the test problems is given
in Fig. 2.
3.1. Melting problems
The melting problems solved are those considered by
Gobin and Le Quere (2000). They compared the numer-
ical solutions presented by 13 contributors using dier-
ent methods and software. They make no comparison
with experimental results because, to the best of their
knowledge, reliable results are not available. Here we
solve three of the four cases presented by Gobin and
Le Quere (2000), keeping their original numbering.
The case that is not solved has a very low Ra number
and melting is dominated by conduction. It is very close
to the Stefan problem, in which convection is neglected.
The latter has been solved for all considered test cases as
a rst check for the moving mesh formulation.
The transition from the dimensionless numbers to the
material data used in the presented equations is as
Table 1
Melting problems: material data (in terms of dimensionless numbers)
Pr Ste Ra s
Case 2 0.02 0.01 2.5 10
5
0.04
Case 3 50 0.1 1.0 10
7
0.01
Case 4 50 0.1 1.0 10
8
0.01
Table 2
Melting without convection. Comparison of FEM moving mesh with analytical results
Mesh l
1
Time steps Error % at s Max. abs. error%
Case 2 Case 3 Case 4 Case 2 Case 3 Case 4
15 1 200 300 0.065 0.076 0.076 0.896 0.987 0.987
600 0.224 0.055 0.055 0.987 0.641 0.641
600 300 0.206 0.037 0.037 0.301 0.288 0.288
600 0.205 0.026 0.026 0.301 0.388 0.388
45 1 200 300 0.082 0.076 0.076 0.598 0.383 0.383
600 0.016 0.007 0.007 0.515 0.313 0.313
600 300 0.267 0.050 0.050 0.622 0.254 0.254
600 0.029 0.005 0.005 0.331 0.195 0.195
602 R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612
follow: q
0
=

Ra=Pr
_
[kg/m
3
]; c = Pr [J/kg C]; h = Pr/
Ste [J/kg]; t =

RaPr
_
Ste
s [s]; l = 1 [kg/m s]; b = 1 [1/C];
k = 1[J/(m s C)]; g = 1 [m/s
2
]; T
w
= 1 C; T
0
= 0 C;
L = 1 [m]; q = q
0
(1 b(T
w
T)) [kg/m
3
]. These num-
bers are given in Table 1. The time s is the time when
results are compared.
3.1.1. Melting without convection
The test problems are rst solved neglecting convec-
tion. The results are compared with the one-dimensional
analytical solution of the Stefan problem, used also by
Mackenzie and Robertson (2000). The nite element
mesh is a single row of rectangular elements in the x-
direction. The error in the position of the melting front
at time s and the maximum error from all time steps are
presented in Table 2.
Error% =
x
FEM
x
Analytical
L
x
100% (31)
To smooth out the initial temperature discontinuity
between the hot wall and the interior, smaller time steps
are used in the beginning of the analysis. The rst 100
time steps increase in a geometric progression with a fac-
tor 1.07. The remaining time steps are equal to the 100th
one. The phase change temperature interval (PCTI) is
e = e
1
+ e
2
with e
1
= e
2
= 0.01. The initial temperature
in the solid and the boundary condition on the right wall
are dropped down by 1.1e
1
in order to accommodate the
whole PCTI.
Fig. 3. Case 2: Mesh, temperature and velocity at time s = 0.04 (left) and s = 0.08 (right).
R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612 603
The FEM moving mesh results are in very good
agreement with the analytical solution and do not
depend on the computational parameters when they
are varied within reasonable bounds. The latter are
determined by two combinations of the parameters:
(i) The number of elements, the level of the mesh
renement governed by l
1
, and the magnitude of
e should ensure that there are at least 34 elements
across the width of the PCTI. Thus the numerical
integration of the large discontinuity in the appar-
ent heat capacity can be done suciently accu-
rately. The error occurs in the elements through
which the isotherms (T
m
e
1
) or (T
m
+ e
2
) pass.
Increasing the number of the Gauss integration
points for these elements does not help much
because the convergence rate is very low.
(ii) The magnitudes of the time step and e should
ensure that the area under the PCTI progresses
by overlapping its previous position and does not
leave any volumes that change phase without
being exposed to the eect of the latent, i.e. being
in a mushy state.
The PCTI is the most important parameter. If it has
small values it is more likely not to meet the above
requirements and the solution would be more inaccu-
rate. However, using large values for the PCTI modies
Fig. 4. Case 3: Mesh, temperature and velocity at time s = 0.01 (left) and s = 0.02 (right).
604 R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612
the actual physical problem, which can become unac-
ceptable, especially, when the phase change is isother-
mal. The need for a proper choice of the PCTI may be
avoided if the apparent heat capacity is dened by
means of the enthalpy H
qc =
dH
dT
~
dH=dx
dT=dx
(32)
where the approximation is written for the 1D formula-
tion of the problem. Having to approximate dH/dx and
dT/dx discretely means that the phase change cannot be
omitted. Furthermore, in a nite element context this
approach takes into account the inuence of all points
in the element and this has a smoothing eect which
again can be reduced if the elements close to the melting
isotherm are small.
However, a rened mesh at the phase change front
will require small time steps to avoid the large jumps
of the mushy phase, as discussed in (ii). It is possible
to have good engineering solutions (error < 5%) of 1D
test problems with a suitable combination of large time
steps and small element size in which a large volume
(more than 50%) undergoes phase transition without
passing through the mushy state. This accuracy can be
achieved because the energy error accumulated in this
volume can be partially cancelled out by the energy
error in the mushy volume, which has been in a mushy
state for longer than necessary.
Fig. 5. Case 4: Mesh, temperature and velocity at time s = 0.005 (left) and s = 0.01 (right).
R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612 605
In the case of general 2D phase change problems var-
ious spatial approximations of Eq. (32) have been pro-
posed (Dalhuijsen and Segal, 1981). Our numerical
studies have shown that in these cases it is more dicult
to get the above error canceling eect and that, on aver-
age, better accuracy is achieved when the apparent heat
capacity is given in terms of the phase change tempera-
ture interval, Eq. (5).
3.1.2. Melting with natural convection
The moving mesh, temperature and velocity at two
time steps for Cases 24 are shown in Figs. 35, respec-
tively. In Case 2 the velocity eld is dominated by the
presence of recirculation rolls. There are more rolls at
the early stages of the melting process when the aspect
ratio y/x of the liquid region is bigger. As the melting
front propagates towards the interior the rolls merge
and decrease in number. These rolls cause the wavy
shape of the melting front and they are numerically pre-
dicted by most analyses presented by Gobin and Le
Quere (2000). It has been established by a careful mesh
renement study (Hannoun et al., 2003) that this veloc-
ity pattern is the correct solution of the mathematical
problem, although the existence of such rolls has not still
been reported experimentally. In Case 3 recirculation
rolls are observed at the very early stages of the melting
and they soon collapse into one circulation. For Case 4,
with the meshes and time steps used, recirculation rolls
are not observed. These two cases are characterized by
high velocities (strong convective contribution in the
transport equation) and very thin thermal boundary lay-
ers. Computationally they are more demanding (ner
meshes and smaller time steps) and more dicult to
solve than Case 2.
In Fig. 6 the computed phase change front using the
moving mesh technique (solid lines) and the range of
variation of the numerical results (dotted lines) given
by Gobin and Le Quere (2000) are plotted. They dene
a range of best solutions, which is presented in Fig. 6 by
the dashed lines. For the three cases the moving mesh re-
sults lie within the range of the best solutions. Only two
results show such accuracyby Le Quere, results No 10,
and by Wintru, results No 13. The former uses a one-
domain, enthalpy formulation and the latter uses a
front tracking technique. Two more results using the
one-domain formulation are within the best solutions
for Case 2 and Case 3, but have not solved Case 4.
The meshes NX NY and the time steps used are pre-
sented in Table 3.
The additional user denable data used in the FEM
analyses are give in Table 4. The initial and the right wall
boundary temperature T
0
is lowered by 1.5e
1
in order to
accommodate the phase change interval. Numerical
studies have shown that better accuracy is achieved when
the phase change interval lies entirely within the temper-
ature of the phase that is going to change phase, i.e. into
0 0.1 0.2 0.3 0.4 0.5 0.6
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Case 2
= 0.04
Case 3
= 0.01
Case 4
= 0.01
Y
Mesh
40x60
Mesh
50x60
Mesh
40x60
X X X
Fig. 6. Position of the phase change front at time s. Solid lineFEM moving mesh solution. Area between dotted linesrange of all numerical
solutions by Gobin and Le Quere (2000); Area between dashed linesrange of best numerical solutions by Gobin and Le Quere (2000).
Table 3
Mesh and time steps used in FEM moving mesh analyses and in the best solutions given by Gobin and Le Quere (2000)
Moving mesh Le Quere (No. 10) Couturier-Sadat (No. 6) Medale (No. 11)
Mesh Ds Mesh Ds Mesh Ds Mesh Ds
Case 2 30 30 2.5 10
5
128 192 2 10
7
202 202 1 10
5
100 50 1 10
4
Case 3 40 40 6.6 10
6
No data 122 122 2 10
6
200 100 1 10
5
Case 4 50 60 2.5 10
6
192 192 1 10
7
Not solved Not solved
606 R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612
the solid for melting problems and into the liquid for
solidication problems. The values of l
1
may be inter-
changed without aecting the accuracy, which is later
shown in the sensitivity study. Case 4 has narrower ther-
mal boundary layers and thus needs more mesh rene-
ment towards the boundaries, i.e. higher values for l
2
.
In Fig. 7 the evolution of the averaged Nusselt num-
ber on the hot wall is presented. The numerical results
presented by Gobin and Le Quere (2000) show again a
large scatter. The moving mesh results are close to the
empirical correlation for Case 3 and Case 4. For Case
2 they are very close to Le Queres results, which form
the upper bound of the numerical solutions.
In Figs. 8 and 9 the sensitivity of the position of the
phase change front on the variation of some user den-
able parameters can be observed for Case 2 and Case 3,
respectively. Similar low sensitivity behaviour is com-
puted for Case 4 and that is why these results are not
presented. The time for the numerical simulation s is
given in Table 1 and the computational parameters that
are not varied explicitly are those given in Table 4. To
extend the sensitivity study dierent number of times
steps, N
s
, are used as followCase 2: N
s
= 1000 for
problems in Fig. 8A, N
s
= 300 for Fig. 8B and C; Case
3: N
s
= 1500 for Fig. 9A and C, N
s
= 700 for Fig. 9B.
It can be seen that the parameters can be varied in a
wide range without having a considerable eect on the
position of the phase change front. The interaction
between the parameters and the reasonable bounds of
their variation are the same as discussed previously in
0 0.002 0.004 0.006 0.008 0.01
10
20
30
40
50
60
0 0.002 0.004 0.006 0.008 0.01
10
16
22
28
34
40
0 0.01 0.02 0.03 0.04
3
5
7
9
11
13
Nu Nu Nu

Fig. 7. Evolution of the averaged Nusselt number on the hot wall. Solid lineFEM results; Dashed linesrange of variation of numerical results by
Gobin and Le Quere (2000). Dotted linetarget result (empirical correlation) by Gobin and Le Quere (2000).
Table 4
Computational parameters used in the FEM moving mesh analyses
e
1
e
2
l
1
l
2
l
L
Case 2 0.02 0.0 300 50 1.2
Case 3 0.02 0.0 100 50 2.0
Case 4 0.05 0.0 200 100 1.5
0.24 0.28 0.32 0.36 0.4
0
0.2
0.4
0.6
0.8
1
0.24 0.28 0.32 0.36 0.4 0.24 0.28 0.32 0.36 0.4 0.24 0.28 0.32 0.36 0.4
D A B C
Y
X X X X
Fig. 8. Case 2: Position of the phase change frontsensitivity study. (A) 1: mesh 40 60, 2: mesh 30 30, 3: mesh 20 25. (B) 1: l
1
= 500, 2:
l
1
= 100, 3: l
1
= 20. (C) 1: e
1
= 0.02, 2: e
1
= 0.10, 3: e
1
= 0.01. (D) 1: N
s
= 1200, 2: N
s
= 600, 3: N
s
= 200. Graphs: 1: solid lines, 2: dashed
lines, 3: dotted lines.
R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612 607
(i) and (ii). The presence of uid ow imposes two addi-
tional limitations. First, there should be enough ele-
ments in the uid phase to model correctly the uid
ow. For the same reason the mesh should not be over
rened at the phase change interval. Second, the time
step should be small enough to allow correct modeling
of the uid ow. When a small number of time steps
are used the velocity becomes more unstable between
time steps (e.g. it has a wavy pattern along the top
boundary). This velocity instability seems to be the pri-
mary reason for the variation in the position of the
phase change front in Fig. 8D. When the uid ow is
ignored (see results in Table 2, computed with a small
number of time steps) the errors are very small.
Increasing the PCTI in Fig. 8C from e
1
= 0.01 to
e
1
= 0.02 has almost no inuence on the predicted posi-
tion of the phase change front. A decrease of e
1
(not
shown on the graph, but observed during the numerical
sensitivity studies) causes the front to propagate faster
because the number of elements in the mushy region is
reduced and the latent heat is not accounted for prop-
erly. A large e
1
reduces the accuracy, too, because the
physical problem has been changed (e
1
= 0.1 is 10% of
the temperature range in the problem).
0 0.1 0.2 0.3 0.4 0.5
0
0.2
0.4
0.6
0.8
1
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
D A B C
Y
X X X X
Fig. 9. Case 3: Position of the phase change frontsensitivity study. (A) 1: mesh 70 110, 2: mesh 40 40, 3mesh 30 60. (B) 1: l
1
= 200, 2:
l
1
= 100, 3: l
1
= 50. (C) 1: e
1
= 0.02, 2: e
1
= 0.10, 3: e
1
= 0.01. (D) 1: N
s
= 1000, 2: N
s
= 500, 3: N
s
= 250. Graphs: 1: solid lines, 2: dashed
lines, 3: dotted lines.
Fig. 10. Water freezing: Mesh and temperature at time t = 1 h and 5 h and velocity at t = 1 h.
608 R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612
3.2. Freezing problem
Freezing of water on the vertical wall of a rectangular
cavity is studied. The outline of the problem is presented
in Fig. 2. It has been studied experimentally by Braga
and Viskanta (1992) and numerically by Scanlon and
Stickland (2004). The water density as a function of T
[C] is approximated as (Scanlon and Stickland, 2004):
q(T)
= 999.972(1 9.2793 10
6
[T 4.0293[
1.894816
) [kg=m
3
[
(33)
The moving mesh of 60 by 40 elements (along x and
y, respectively), the temperature isotherms at times
t = 1 h and t = 5 h, and the velocity at t = 1 h are shown
in Fig. 10. The density inversion at 4.0293 C is respon-
sible for the ow structure with a big anti-clockwise cir-
culation where T > 4.0293 C and a small clockwise one
at the bottom left where T < 4.0293 C. The same pat-
tern has been observed experimentally by Braga and
Viskanta (1992) and numerically by Scanlon and Stick-
land (2004).
The position of the phase change front at several
times is shown in Fig. 11. At time t = 1 h and t = 2 h
there is a very good agreement with the experimental re-
sults by Braga and Viskanta (1992). At t = 5 h the agree-
ment is generally good with the exception that the
computed ice thickness is larger. It is possible that the
error is due to computed lower velocities near the freez-
ing boundary (in Braga and Viskanta (1992) there is no
data of the magnitude of the experimental velocities).
The temperature histories at the location of two ther-
mocouples, TC1 (x = 0.05 m, y = 0.02 m) and TC2 (x =
0.15 m, y = 0.075 m) are shown in Fig. 12. There is
good agreement with the experimental results by Braga
and Viskanta (1992), and the numerical results by
Scanlon and Stickland (2004). It should be noted that
the latter are achieved by nite volume method (Fluent
software) using a mesh of 300 200 elements which was
obtained via a grid renement study.
In Fig. 13 the results of a sensitivity study are shown.
The error of the position of the phase change front is
computed by comparing the FEM and experimental x-
coordinates of the freezing front at time t = 60 min, at
50 pre-dened, uniformly distributed locations along y.
Error =

50
i=1
(x
FEM
x
exp
)
2
_ _1
2

50
i=1
x
exp
_ _
1
100% (34)
It can be seen that the accuracy decreases when
l < 200 and it is almost independent on the number time
steps and number of elements used. The investigation of
the phase change temperature interval e shows that it
has an optimal value of about 1.5 C and that better
accuracy is achieved when e
1
= 0, i.e. when e is situated
entirely into the phase that is going to change phase.
4. Conclusions
An r-adaptive moving mesh technique has been suc-
cessfully applied to the solution of phase change problems
with natural convection. The generation of the moving
mesh requires the solution of two extra partial dierential
equations at each time step. The main advantage of the
method is that less number of elements may be used.
The accuracy is preserved because the mesh is kept su-
ciently rened around the phase change at all times. The
numerical experiments in this paper have been computed
on rectangular domains. For more complicated non-
0 0.02 0.04 0.06
0
0.03
0.06
0.09
0.12
0.15
Freezing front
2h
5h
Y

[
m
]

1h
X [m]
Fig. 11. Water freezing: Position of the phase change front at times
t = 1 h, 2 h and 5 h. Solid line: FEM results (moving mesh 60 40).
Dotted line (d): Experimental data by Braga and Viskanta (1992).
0 10 20 30 40 50 60
0
2
4
6
8
10
12
T

[

C
]
TC2
TC1
Time [min]
Fig. 12. Water freezing. Comparison with published results. Temper-
ature at thermocouples TC1 and TC2. Solid lineFEM results
(moving mesh 60 40); Dots oExperimental data by Braga and
Viskanta (1992); Dashed lineNumerical results by Scanlon and
Stickland (2004), FLUENT analysis (nite volume method) with a
xed mesh 300 200 elements.
R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612 609
convex domains the moving mesh method can be ex-
tended as shown by Cao et al. (1999), and we leave work
in this direction for further investigation.
Acknowledgement
The authors would like to acknowledge the EPSRC
nancial support, grant No. GR/R2603/01.
Appendix A
Moving mesh equations on the boundary edges:
os
o
~
t
=
1
P
s
o
of
G
os
of
_ _
=
1
P
s
G
o
2
s
of
2

oG
on
os
of
_ _
=
1
P
s
A
s
o
2
s
of
2
D
s
os
of
_ _
(A:1)
A
s
= G; D
s
=
oG
on
; P
s
=

A
2
s
D
2
s
_
Moving mesh equations in the interior:
ox
o
~
t
=
1
P
A
o
2
x
on
2
B
o
2
x
onog
C
o
2
x
og
2
D
ox
on
E
ox
og
_ _
(A:2)
Coecients:
A =
1
J
2
G
ox
og
ox
og

oy
og
oy
og
_ _
;
C =
1
J
2
G
ox
on
ox
on

oy
on
oy
on
_ _
(A:3)
B =
2
J
2
G
ox
on
ox
og

oy
on
oy
og
_ _
(A:4)
D =
1
J
2
G
2
ox
og
ox
og

oy
og
oy
og
_ _
oG
on
_

ox
on
ox
og

oy
on
oy
og
_ _
oG
og
_
(A:5)
E =
1
J
2
G
2
ox
on
ox
on

oy
on
oy
on
_ _
oG
og
_

ox
on
ox
og

oy
on
oy
og
_ _
oG
on
_
(A:6)
P =

A
2
C
2
D
2
E
2
_
; J =
ox
on
oy
og

ox
og
oy
on
(A:7)
0 2000 4000
0
0.5
1
1.5
0 500 1000 1500 2000
0
0.5
1
1.5
0 200 400 600
0
0.5
1
1.5
Error % Error % Error %
2
=1
1
=0
2
=2
1
=0
=600
=300
=600
=300
of time steps of elements
0 1 2 3
0
0.5
1
1.5
0 1 2 3
0
0.5
1
1.5
0 1 2 3
0
0.5
1
1.5
Error % Error %
Error %
1
=0
1
=
2 2
=0
=600
=600
=600
=300
=300
=300
2
(
1
+
2
)
1

Fig. 13. Sensitivity study of the accuracy of the position of the freezing front at t = 60 min; e
1
and e
2
dene the variation of the phase change interval:
T
m
e
1
6 T 6 T
m
+ e
2
; l is a parameter in the monitor function G, Eq. (18), dening the mesh density towards the phase change front (the bigger
l the more rened is the mesh).
610 R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612
Finite element approximation (on an element level):
M
dx
dt
(A B C D E)x = 0 (A:8)
where
M =
_ _
X
N
T
NdX (A:9a)
A =
_ _
X
oN
T
on
A
P
oN
on
dX (A:9b)
B =
1
2
_ _
X
oN
T
on
B
P
oN
og
dX
_ _
X
oN
T
og
B
P
oN
on
dX
_ _
(A:9c)
C =
_ _
X
oN
T
og
C
P
oN
og
dX (A:9d)
D =
_ _
X
N
T
1
2
o(B=P)
og

o(A=P)
on

D
P
_ _
oN
on
dX (A:9e)
E =
_ _
X
N
T
1
2
o(B=P)
on

o(C=P)
og

E
P
_ _
oN
og
dX (A:9f)
Appendix B
NavierStokes equations:
Du
Dt
u \u
dx
dt
\u \
T
l
q
0
\u
_ _

k
q
0
o
ox
ou
ox

ov
oy
_ _
= 0
Dv
Dt
u \v
dx
dt
\v \
T
l
q
0
\v
_ _

k
q
0
o
oy
ou
ox

ov
oy
_ _
=
(q
0
q)g
q
0
(B:1)
Finite element approximation (on an element level):
C
V
ou
ot
ov
ot
_ _
(K
V
U
V
L
V
)
u
v
_ _
=
F
x
F
y
_ _
(B:2)
C
V
=
_ _
X
N
T
NdX 0
0
_ _
X
N
T
NdX
_ _
(B:3)
K
V
=
_ _
X
l
q
0
oN
T
ox
oN
ox

oN
T
oy
oN
oy
_ _
dX 0
0
_ _
X
l
q
0
oN
T
ox
oN
ox

oN
T
oy
oN
oy
_ _
dX
_

_
_

_
(B:4)
U
V
=
_ _
X
N
T
u
oN
ox
v
oN
oy
_ _
dX 0
0
_ _
X
N
T
u
oN
ox
v
oN
oy
_ _
dX
_

_
_

_
(B:5)
L
V
=
k
q
0
_ _
X
oN
T
ox
oN
ox
dX
_ _
X
oN
T
ox
oN
oy
dX
_ _
X
oN
T
oy
oN
ox
dX
_ _
X
oN
T
oy
oN
oy
dX
_
_
_
_
(B:6)
F
x
F
y
_ _
=
_ _
X
N
T dx
dt
\udX
_ _
X
N
T (q
0
q)g
q
0

dx
dt
\v
_ _
dX
_
_
_
_
_
_
(B:7)
Notes:
(i) First only element matrices K
V
and U
V
are com-
puted and assembled. Let z be the value of the
maximum by absolute value element in the assem-
bled matrix. Let p be the number of digits in the
computer presentation of a real number. In order
half of the signicant digits of

(K
V
+ U
V
) to be
preserved after assembling all L
V
, the penalty
number k is computed as k = exp(log(z) + p/2).
(ii) The global penalty matrix

L
V
needs to be singu-
lar so matrices L
V
are computed by reduced inte-
gration, one point Gauss quadrature.
(iii) For a fully solid element (one for which the maxi-
mum nodal temperature is less than the melting
temperature, i.e. maxT
N
< T
m
) the matrices in
Eq. (30) are not computed. For a partially solid
element (i.e. minT
N
< T
m
< maxT
N
) an articial
viscosity l
solid
= 1000l is used if the Gauss point
is in a solid state.
(iv) All nodes for which T < T
m
(i.e. they are in a solid
state) are prescribed zero velocities, in addition to
the boundary nodes with no-slip boundary condi-
tions.
References
Beckett, G., Mackenzie, J.A., Robertson, M.L., 2001. A moving mesh
nite element method for the solution of two-dimensional Stefan
problems. J. Comput. Phys. 168, 500518.
Beckett, G., Mackenzie, J.A., Ramage, A., Sloan, D.M., 2002.
Computational solution of two-dimensional unsteady PDEs using
moving mesh methods. J. Comput. Phys. 182, 478495.
Braga, S.L., Viskanta, S., 1992. Eect of density extremum on the
solidication of water on a vertical wall of a rectangular cavity.
Exp. Therm. Fluid Sci. 5, 703713.
Cao, W., Huang, W., Russell, R.D., 1999. An r-adaptive nite element
method based upon moving mesh PDEs. J. Comput. Phys. 149,
221244.
Crank, J., 1984. Free and Moving Boundary Problems. Clarendon
Press, Oxford.
Crivelli, L.A., Idelsohn, S.R., 1986. A temperature-based nite element
solution for phase-change problems. Int. J. Numer. Methods Eng.
23, 99119.
Dalhuijsen, A.J., Segal, A., 1981. Comparison of nite techniques for
solidication problems. Int. J. Numer. Methods Eng. 23, 8196.
Fachinotti, V.D., Cardona, A., Huespe, A.E., 1999. A fast convergent
and accurate temperature model for phase change heat conduction.
Int. J. Heat Mass Transfer 44, 18631884.
Gobin, D., Le Quere, P., 2000. Melting from an isothermal vertical
wall. Synthesis of a numerical comparison. Comput. Assist. Mech.
Eng. Sci. 7, 289306.
R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612 611
Gupta, S.C., 2000. A moving grid numerical scheme for multi-
dimensional solidication with transition temperature range.
Comput. Methods Appl. Mech. Eng. 189, 525544.
Hannoun, N., Alexiades, V., Mai, T.Z., 2003. Resolving the contro-
versy over tin and gallium melting in a rectangular cavity heated
from the side. Numer. Heat Transfer, Part B 44, 253276.
Heinrich, J.C., Pepper, D.W., 1996. The Intermediate Finite Element
Method. Taylor and Francis.
Huang, W., 1999. Practical aspects of formulation and solution of
moving mesh partial dierential equations. Mathematics Research
Report No. 99-11-02, University of Kansas.
Lewis, R.W., Roberts, P.M., 1987. Finite element simulation of
solidication problems. Appl. Sci. Res. 44, 6192.
Lynch, D.R., 1985. Unied approach to simulation on deforming
elements with application to phase change problems. J. Comput.
Phys. 57, 303317.
Lynch, D.R., ONeill, K., 1981. Continuously deforming nite
elements for the solution of parabolic problems with and without
phase change. Int. J. Numer. Methods Eng. 17, 8196.
Mackenzie, J.A., Robertson, M.L., 2000. The numerical solution of
one-dimensional phase change problems using an adaptive moving
mesh method. J. Comput. Phys. 161, 537557.
Miller, K., Miller, R.N., 1981. Moving nite elements. SIAM J.
Numer. Anal. 18, 10191057.
Morgan, K., Lewis, R.W., Zienkiewicz, O.C., 1978. An improved
algorithm for heat conduction problems with phase change. Int. J.
Numer. Methods Eng. 12, 11911195.
Nedjar, B., 2002. An enthalpy-based nite element method for non-
linear heat problems involving phase change. Comput. Struct. 80,
921.
Nigro, N., Huespe, A., Fachinotti, V., 2000. Phasewise numerical
integration of nite element method applied to solidication
processes. Int. J. Heat Mass Transfer 43, 10531066.
Ockendon, J.R., Hodgkins, W.R. (Eds.), 1975. Moving Boundary
Problems in Heat Flow and Diusion. Oxford University Press.
Pardo, E., Weckman, D.C., 1990. A xed grid nite element technique
for modeling phase change in steady state conductionadvection
problems. Int. J. Numer. Methods Eng. 29, 969984.
Rolph, W.D., Bathe, K.J., 1982. An ecient algorithm for analysis of
nonlinear heat transfer with phase change. Int. J. Numer. Methods
Eng. 18, 119134.
Salcuden, M., Abdullah, Z., 1988. On the numerical modeling of heat
transfer during solidication processes. Int. J. Numer. Methods
Eng. 28, 445473.
Scanlon, T.J., Stickland, M.T., 2004. Anumerical analysis of buoyancy-
driven melting and freezing. Int. J. Heat Mass Transfer 47, 429436.
Tamma, K.K., Namburu, R.R., 1990. Recent advances, trends and
new perspectives via enthalpy-based nite element formulations for
applications to solidication problems. Int. J. Numer. Methods
Eng. 30, 803820.
Thomas, B.G., Samarasekara, I.V., Brimacombe, J.K., 1984. Com-
parison of numerical modeling techniques for complex two-
dimensional, transient heat conduction problems. Metall. Trans.
B 15, 307318.
Voller, V., Cross, M., 1981. Accurate solutions of moving boundary
problems using the enthalpy method. Int. J. Heat Mass Transfer
24, 545556.
Voller, V.R., Swaminathan, C.R., Thomas, B.G., 1990. Fixed grid
techniques for phase change problems: a review. Int. J. Numer.
Methods Eng. 30, 875898.
612 R.T. Tenchev et al. / Int. J. Heat and Fluid Flow 26 (2005) 597612