c Indian Academy of Sciences

Sadhana Vol. 38, Part 6, December 2013, pp. 12591285. 

Modified enthalpy method for the simulation of melting

and solidification
NIRANJAN N GUDIBANDE and KANNAN N IYER
Department of Mechanical Engineering, Indian Institute of Technology Bombay,
Mumbai 400 076, India
e-mail: [email protected]; [email protected]
MS received 18 February 2012; revised 19 February 2013; accepted 6 March 2013
Abstract. Enthalpy method is commonly used in the simulation of melting and
solidification owing to its ease of implementation. It however has a few shortcomings.
When it is used to simulate melting/solidification on a coarse grid, the temperature
time history of a point close to the interface shows waviness. While simulating melting
with natural convection, in order to impose no-slip and impermeability boundary conditions, momentum sink terms are used with some arbitrary constants called mushy
zone constants. The values of these are very large and have no physical basis. Further,
the chosen values affect the predictions and hence have to be tuned for satisfactory
comparison with experimental data. To overcome these deficiencies, a new cell splitting method under the framework of the enthalpy method has been proposed. This
method does not produce waviness nor requires mushy zone constants for simulating
melting with natural convection. The method is then demonstrated for a simple onedimensional melting problem and the results are compared with analytical solutions.
The method is then demonstrated to work in two-dimensions and comparisons are
shown with analytical solutions for problems with planar and curvilinear interfaces.
To further benchmark the present method, simulations are performed for melting in a
rectangular cavity with natural convection in the liquid melt. The solidliquid interface obtained is compared satisfactorily with the experimental results available in
literature.
Keywords.

Melting; enthalpy method; wavy interface; mushy zone constant

1. Introduction
The study of melting and solidification offers insights in the design of casting, welding, latent
thermal energy storage systems, etc., and in the analysis of melting of glaciers, solidification of
lava, thawing of food, etc. The motivation for the present work arose from a need to analyse

For correspondence

1259

1260

Niranjan N Gudibande and Kannan N Iyer

melting of lead in a metallic cask due to fire. Numerical tools are valuable in such a study.
Numerical techniques for simulating melting and solidification must be capable of handling the
presence of a moving solidliquid interface (an internal boundary) across which the properties
may be discontinuous. They must also be capable of handling natural convection in the molten
liquid. Several numerical methods exist for simulating melting/solidification and phase change
phenomena in general. These methods can be classified into deforming grid, front tracking and
enthalpy methods.
In a deforming grid method, meshes are created for the solid and the liquid domains separately
with the solidliquid interface as a common boundary. Governing equations are solved for temperature and velocity fields with appropriate boundary conditions at the solidliquid interface.
Stefan condition (an energy balance at the solidliquid interface) is used to obtain the velocity
with which the interface moves as the material undergoes melting/solidification. This velocity
is integrated in time to obtain the new location of the interface. Deforming grid method requires
creation of grids at each time step. This method has been used for the simulation of melting
of pure metals with natural convection in the molten liquid by Lacroix (1989), Cerimele et al
(2002) and Kumar et al (2006).
In a front tracking method, the interface is represented by a set of marker points connected
by curves (in 2d) or triangular elements (in 3d). These marker particles move on a fixed grid
as the material undergoes melting or solidification. The governing equations are solved on the
fixed grid and for nodal points which are close to the interface finite difference stencils which
involve nodal points of the same phase and one or more points on the interface are used. The
velocity of the marker points are obtained by applying the Stefan condition. As the material
undergoes melting or solidification, these marker points may form clusters or move far apart
and thus redistribution of the marker particles becomes necessary. Udaykumar et al (1996) and
Li et al (2003) have demonstrated this method for melting and solidification. Both the above
mentioned methods require special treatment when the interface undergoes a topological change
due to merging or splitting of interface.
In the enthalpy method, the energy equations of the solid and the liquid domain and the
boundary conditions on the interface are combined into a single equation called the enthalpy
equation. By solving the enthalpy equation the temperature and the liquid fraction fields ( f ) can
be obtained. Liquid fraction is the ratio of volume of liquid in a control volume (CV) to the
total volume of the CV. Thus the solid and the liquid control volumes have f = 1 and f = 0,
respectively and control volumes which have f (0, 1) contain the solidliquid interface. The
enthalpy method is both simple and robust. It is simple in the sense that no complicated algorithms are required for recreation of grids (as in deforming grid method) and in finding suitable
stencils for finite difference approximation (as in front tracking method). It is also robust as it
does not require any special treatment when the interface undergoes topological changes. This
method has been used by several researchers such as Gartling (1980), Voller et al (1987), Brent
et al (1988), Cao & Faghri (1990).
The enthalpy method however has some shortcomings. When enthalpy method is used to simulate melting/solidification on a coarse grid and the temperature of a point close to the interface
is plotted as a function of time, it shows a wavy character. This has been attributed to the fact
that the temperature of the control volume containing the interface is imposed to be at the melting
temperature rather than at the actual location of the solidliquid interface. Since the discovery
of this phenomena by Price & Slack (1954) several methods have been proposed in the past
to eliminate the waviness. These include the implicit time stepping method of Voller & Cross
(1981), explicit enthalpy method of Tacke (1985), centroidal temperature correction method
of Date (1994) and the effective conduction length method of Kim & Anghaie (2001). The

Modified enthalpy method

1261

implicit time stepping and the effective conduction length model are limited for one-dimensional
problems. In the centroidal temperature correction method the interface becomes diffuse.
In simulating melting with natural convection, the impermeability and the no-slip condition
has to be imposed at the solidliquid interface. Since the exact location of the solidliquid interface is unknown in the enthalpy method, these boundary conditions are enforced indirectly.
Authors such as Voller et al (1987), Brent et al (1988) have used momentum sink terms in the
momentum equations, while Gartling (1980) and Cao & Faghri (1990) have used a variable viscosity method. In variable viscosity method, viscosity is written as a function of liquid fraction.
This function ensures that the viscosity of the solid is very large and approaches that of the
molten liquid when f = 1.
The sink terms and variable viscosity models use large numbers varying from 105 107 as
constants (termed mushy constants). Very few studies are available on the sensitivity of the
solution to these constants. In a recent study by Shmueli et al (2010), the effect of these mushy
constants on the solution was examined. The authors concluded that there exists an optimum
value for the mushy constant which would match the experimental results well. Thus, there exists
some arbitrariness in specifying these constants. To overcome these deficiencies, Niranjan &
Iyer (2011) had recently proposed a cell splitting method. This paper provides complete details
of the cell splitting method.

2. The governing equation

The continuity, momentum and the energy equations can be used to obtain the temperature and
velocity fields in the liquid and the energy equation can be used to obtain the temperature field in
the solid. The molten liquid is assumed to be incompressible and Boussinesq approximation for
natural convection is assumed to be valid. The integral forms of the governing equations along
with the appropriate boundary conditions at the solidliquid interface are listed below.
Continuity equation


V.nd
A = 0.

(1)

Momentum equation








A = P.nd
A + (V).nd
A+
g T Tr e f d V .
(V)d V + VV.nd
t
V

(2)
Energy equations
In liquid domain

V

In solid domain

(C P T )d V +
t

V


C P T V.nd
A=

(C P T )d V =
t

kT.nd
A.

(3)


kT.nd
A.
A

(4)

1262

Niranjan N Gudibande and Kannan N Iyer

At the solidliquid interface, the velocity will be zero and the temperature will be set at the
melting point. An additional equation can be derived by applying the energy balance at the
interface, (known as the Stefan condition)




(5)
kT.n S kT.n L = |Vn | H M .
Here, Vn denotes the velocity of the solidliquid interface. The direction of this velocity is
normal to the interface. The enthalpy method combines the energy equations (3) and (4) along
with the boundary conditions into a single equation for enthalpy



d
hd V + hV.nd
A = kT.nd
A.
(6)
dt
V

Here, h denotes the enthalpy. The enthalpy h can be related to the temperature through the
following relation,

h<0
CPS
0 < h < H M
T TM = 0
(7)

hHM
h > H M .
CPL

Here enthalpy of the solid at its melting point is chosen as the reference enthalpy and is assigned
a value of zero.

3. Overall methodology
In the present method, the control volumes (CV) whose liquid fraction ( f ) lie in between 0 and 1
are identified. These CVs will henceforth be called as mushy CVs. Using a suitable interface
construction method and the f distribution, a piece-wise linear interface is constructed. The constructed interface is used to split the mushy CV into two single phase CVs as shown in figure 1.
Thus, if there are m, n, r, solid, liquid and mushy CVs, respectively, then the splitting procedure
will divide the domain  into liquid sub-domain  L containing m +r CVs and solid sub-domain
 S containing n + r CVs. By solving the mass, momentum and the energy equations in  L ,

Figure 1. The demonstration of the cell splitting method.

Modified enthalpy method

1263

u, v, and TL are obtained, while TS is obtained by solving the energy equation in  S . The constructed interface will form a boundary for these equations. For the energy equation, a Dirichlet
type of boundary condition is imposed by specifying the interface to be at the melting temperature. The no-slip and the impermeability boundary conditions are imposed for the momentum
and the mass equations on the interface. An additional equation is necessary to obtain the evolution of liquid fraction, f , with time. This is obtained by solving the enthalpy form of the energy
equation for the non-split control volumes.

4. Interface reconstruction
Interface reconstruction is an important step in the present method in which an approximate
solidliquid interface is constructed from the liquid fraction data. In volume of fluid method
(VOF) too, interface reconstruction plays a vital role. There are several interface construction
methods which have been proposed by researchers to simulate two phase flows using VOF methods. A detailed description of several of these methods can be found in Rider & Kothe (1998)
and Pilliod & Puckett (2004). A typical interface constructed by VOF method, however, produces an interface which is discontinuous at the faces of the control volumes. In VOF method
the constructed interface is used to geometrically advect the liquid fractions and the discontinuity does not pose any problem. The present cell splitting method however relies on the existence
of a continuous solidliquid interface and hence a VOF interface reconstruction method cannot
be directly used. To produce a continuous piece-wise linear interface the following procedure is
employed.
The control volumes sharing a common face and having liquid fractions between 0 and 1 are
identified. The solidliquid interface will pass through the common face of these control volume such as the control volume A and B shown in figure 2. The objective is to find the point
of intersection of the interface with this common face. It is assumed that the interface can be
approximated by a straight line segment in the neighbourhood of the common face. To obtain a
unique line representing the solidliquid interface two constraints have to be imposed. By imposing the conditions that the ratios of volume of liquid to that of the solid in the control volumes A

Figure 2. The interface reconstruction method.

1264

Niranjan N Gudibande and Kannan N Iyer

and B in figure 3 are in the ratios of f A : 1 f A and f B : 1 f B , respectively two constraints can
be imposed. Thus a unique straight line can be obtained. This procedure is used for obtaining the
intersection points in both horizontal and the vertical directions. The interface is then obtained
by joining the intersection points by piece-wise straight lines. Boundary conditions are utilized
for finding the boundary intersection point such as point B1 in figure 2. For e.g., if the boundary
condition is temperature specified, the temperature gradient in x direction can be obtained from
the known temperatures at L and M. Similarly, the gradient in y direction can be obtained from
the temperatures at L and N. Using the two components, the direction of gradient vector can be
found. As the interface is an isotherm, its direction will be normal to the gradient vector. Thus,
the slope of the interface can be estimated using which a linear interface can be constructed such
that it divides the control volume in the ratio of f : 1 f . The point of intersection of this constructed interface with the boundary will determine the point B1. Thus, all the points that will
constitute an interface are obtained. The line joining these points defines the interface. The procedure followed for the reconstruction is second order accurate and this is derived in Appendix A.
The constructed interface is then used to obtain geometric properties such as the centroids, the
distances of centroids from the interface, the area of the interface, etc.
The solidliquid (SL) interface moves as the material undergoes melting or solidification. The
SL interface is one of the faces of the split control volume such as the face XY of the control
volume ABXYD shown in figure 3. The movement of the face XY results in the expansion/
contraction of the control volume ABXYD. Thus the governing equation has to be modified to
account for this movement. Following Donea et al (1982), a generic transport equation for a
transport variable  in case of moving grid (interface) can be written as




d
A = q D + S,
d V + V Vg .nd
(8)
dt
V

where q D and S denote the diffusive transport and volumetric source of , respectively. The
time derivative term outside the integral accounts for the expansion/contraction of the control
volume and the Vg .n term (Vg is the velocity of the grid) accounts for the convective contribution
due to the movement of the grids. The surface area A can be split into the original Cartesian
part AC and the interface part A I , for e.g., AC for the control volume ABXYD in figure 3 is
AB BX YD DA and the interfacial area A I is XY. At the solidliquid interface due to the

Figure 3. Expansion of control volume due to motion of the solidliquid interface.

Modified enthalpy method

1265

impermeable boundary condition V.n is zero. Further for the split control volumes the Vg .n term
is non-zero only at the solid/liquid interface.
Incorporating the above mentioned simplifications, equation (8) can be written as,
d
dt

d V
V

Vg .nd
A+
AI

V.nd
A = q D + S.

(9)

AC

The grid movement should also satisfy the space conservation equation. For a split control
volume this can be expressed as,
d
dt


dV
V

Vg .nd
A = 0.

(10)

AI

at the solidliquid interface which represents interface temperature or the velocity components
has a constant value determined by the boundary condition. Let I represent this constant value.
Multiplying the equation (10) by I and subtracting it from the equation (9) results in
d
dt


( I )d V +
V

V.nd
A = q D + S.

(11)

AC

The temperatures and the velocity fields in the liquid domain and the temperature in the solid
domain can be obtained by substituting suitably for in the equation (11). An additional equation is required to obtain the evolution of the interface with time. This is obtained by solving
the enthalpy form of the energy equation. The enthalpy equation is solved on the whole control
volume such as ABCD in figure 3.
d
dt


hd V +
V


hV.nd
A=

AC

kT.nd
A.

(12)

AC

By solving the above equation h

t can be obtained. The enthalpy can be written as a sum of
sensible and the latent heat parts, thus
h = f [C P (T TM ) + H M ] + (1 f ) [C P (T TM )] .

(13)

From the enthalpy field the liquid fraction f can be obtained.

5. Numerical discretisation
The governing equations are discretised on a Cartesian mesh. The grids are non-orthogonal (the
line joining the centroids of the neighbouring control volumes are non-orthogonal to their common face) in the mushy CVs and in the CVs which are close to the solidliquid interface. The calculation of the diffusive flux in these CVs requires treatment which is similar to that in a complex
domain. The details of discretisation of diffusive fluxes in a complex domain are discussed in
Mathur & Murthy (1997). Consider the calculation of the diffusive flux for the face e of the CV
in figure 4. This involves the calculation of the derivative x on the face of the control volume.

1266

Niranjan N Gudibande and Kannan N Iyer

Figure 4. Sample control volume for discretisation.

Let denote a coordinate along the line joining the centroids of CV, viz., P and E. x can
be evaluated as
y y
.
(14)
x =
x
The finite difference for various terms in the above equation can be written as
x =

xE xP
,


y =

yE y P
,


E P
,


y =

A B
,
yA yB

(15)


where  = (x E x P )2 + (y E y P )2 .
It can be seen from the equation (14) that calculation of y involves the values of at points
A and B. The value of B can be calculated from the boundary condition at the solidliquid
interface. The value of A is obtained by inverse distance weighted interpolation from the neighbouring control volumes. For the mushy CV, the existence of a solidliquid interface demands
an additional diffusive flux calculation at the interface. This diffusive flux requires the evaluation of the term n where n is a coordinate normal to the interface. This is obtained by dropping
a perpendicular from the centroid P to the interface and evaluating n as,
P I P
n = 
.
(x P x P I )2 + (y P y P I )2

(16)

The total diffusive transport for a given control volume can thus be written in a standard form as
D=

dnb nb + d I I + SS d P P ,

(17)

E,W,N ,S

where SS denotes the source term due to secondary flux and d I I denotes the contribution from
the solidliquid interface. It is possible that some of the dnb terms may be zero, for e.g., for the

Modified enthalpy method

1267

control volume in figure 4, d S is zero. The convective terms are evaluated by first order upwind
method. The final discrete form of the equation (11) can be written as
( f ( P I )) ( f ( P I )) O
V =
t

anb nb + a I I + SS + S a P P .

(18)

E,W,N ,S

In the above equation, SS and S denote the source term due to secondary flux and other body
source terms, respectively.
5.1 Solution of the NavierStokes equation
The NavierStokes equations are solved on a collocated Cartesian grid. The solution of these
equations requires the calculation of pressure force in addition to the solution of the convection diffusion equation. The final form of the momentum equation after discretisation of the
convective and diffusive terms read as

anb u nb + S + SS + FP X
(19)
aP u P =
E,W,N ,S

aP vP =

anb vnb + S + SS + FPY ,

(20)

E,W,N ,S

where FP X and FPY denotes the pressure force terms in the x and the y directions, respectively.
To calculate the pressure force the pressures at the face of the control volumes are used. The
pressure on a face such as face e in figure 4 is obtained by linear interpolation of the pressures
from centroids P and E. The pressure at the solidliquid interface is obtained by applying
the NavierStokes equation at the interface which results in
P
2un
= 2 ,
n
n

(21)

where n is a coordinate normal to the interface. By finite difference,

PP I = PP + n

2un
.
n 2

(22)

The pressure at the interface can be obtained. The pressure force in the x and y directions are
FP X = Pw Aw Pe Ae (Ae Aw ) PP I .

(23)

FPY = Ps As Pn An (As An ) PP I .

(24)

The areas (Ae Aw ) and (As An ) are the projected areas of the interface in the positive x and
y directions, respectively.
5.2 Momentum interpolation
The governing equations are solved on a collocated grid. Interpolating the velocities at the centroids to obtain the velocities at the faces result in pressure checkerboarding. To overcome this
difficulty, the momentum interpolation of Rhie & Chow (1983) has been used.

1268

Niranjan N Gudibande and Kannan N Iyer

The momentum equations can be written as

u = u +

1
FP X .
ap

(25)

v = v +

1
FPY .
ap

(26)

In order to obtain the velocity at the face e of the control volume P in figure 4 momentum
interpolation is followed. The face velocity can be written as
u f e = u f e +
where


FP X =

1
apf e

FP X ,

(27)

P
P
dV
V f e .
x
x

(28)

P V f e
.
x apf e

(29)

Thus,
u f e = u f e

The value of u f e , a p f e , V f e are obtained by linear interpolation of the corresponding values

from P and E control volumes (figure 4). The derivative term Px is evaluated in a manner
similar to that involved in the calculation of the diffusive term described in section 5. The discrete
form of the above equation is
u f e = u f e

PE PP V f e
+ SS P ,

x a p f e

(30)

where SSP contains the secondary pressure gradient term arising due to the non-orthogonality
of the grid. The pressure at the centroids of the control volumes are obtained by the continuity
equation by using SIMPLE method of Patankar (1980). The face velocity correction term can be
written as,
P  PP V f e
u f e = E
,
(31)

x a p f e
where P  is the pressure correction term. By substituting the face velocities in the continuity equation and incorporating the correction terms, an equation for the pressure correction is
obtained. This pressure correction equation is solved and the pressures and velocities are corrected. The enthalpy equation is solved on the whole control volume. To calculate the diffusive
and convective fluxes, temperatures at the centroids are used. This prevents any interpolation of
temperature across the solidliquid interface.
5.3 Summary of the numerical procedure
Assuming the velocity, temperature and the liquid fraction are available at a previous time, the
following steps are followed
(a) Solve the momentum equations (equations (19) and (20)) and obtain u and v.

(b) Obtain the face velocities by momentum interpolation using equation (30).

Modified enthalpy method

1269

(c) Solve the pressure correction equation and hence correct the face and the centroidal
velocities (equation (31)).
(d) Solve the energy equations to obtain the temperatures in the liquid and the solid domain
using equation (18).
(e) Solve the enthalpy equation to obtain the enthalpy field and hence obtain liquid fraction f
from equation (13).
(f) Using the newly obtained liquid fraction reconstruct the interface.
(g) Solve the above steps until convergence.
The form taken by the system of governing equations for the special case of one-dimension
is mentioned in appendix B. Appendix B also contrasts the significant difference between the
present method and the explicit enthalpy method of Tacke (1985) a popular method which also
eliminates waviness.

6. Results and discussions

To benchmark the cell splitting method and to demonstrate that it can eliminate waviness produced by the enthalpy method, a simple 1-D problem is chosen. The problem consists of a
semi-infinite slab (0 x < ), which is initially solid at its melting temperature. The left end
is raised to a temperature TW which is greater than the melting point of the material TM . The
solidliquid interface moves monotonously to the right.
The governing equations for this situation is given by
T
2T
= 2
t
x

for 0 < x < s(t),

(32)

where s(t) denotes the position of the solidliquid interface at time t. At the solidliquid
interface the boundary conditions are
T = TM

and,
k

(33)

ds(t)
T
= H M
.
x
dt

(34)

By introducing the following non-dimensional variables,

X =

x
,
LS

S=

s
,
LS

t
,
L 2S

(T TM )
.
(TW TM )

(35)

The governing equation becomes

2

.
=

X 2

(36)

=0

(37)

Along with

1 d S(t)
=
.
X
St d

(38)

1270

Niranjan N Gudibande and Kannan N Iyer

At the solidliquid interface here Stefan number St is defined by

St =

C P (TW TM )
.
H M

(39)

The analytical solution for the above problem was given by Neumann and has been reproduced
in Crank (1984). The solution in the non-dimensional form is


X
er f 2

=1
(40)
0 < X < S( ).
er f ()
The position of the interface is given by

S( ) = 2 ,

(41)

where is obtained by solving the equation

St
exp(2 )er f ()

= .

(42)

A one-dimensional finite Cartesian computational domain of length equal to the length scale
L S was chosen and was discretised into 10 equal spaced control volumes. The right end was
imposed to be at the melting temperature. The temperatures of the points at non-dimensional
distances of X * = 0.25, X * = 0.5 and X * = 0.75 are plotted as a function of time in figures 57.
These figures show a comparison of the results obtained from the classical enthalpy method and
the present cell-splitting method with the analytical solutions. It can be seen that the cell splitting
method has eliminated the waviness.
Few comments can be made on the waviness of the enthalpy method. In the enthalpy method,
the centroid of the control volume which is undergoing a phase change is anchored at the melting point. In this duration, the solid and the liquid domains experience constant temperature
boundary conditions and tend to move towards a steady state. This explains the flattening of the
temperature and thus the waviness. The waviness will be more pronounced if the interface takes
more time to cross a control volume. For a given Stefan number the interface velocity decreases
with time. Thus, the waviness for a given Stefan number is more at the location X * = 0.75 as
compared to the location at X * = 0.25. This can be seen from figures 57. Similarly, at a fixed
location, the velocity with which a solidliquid interface moves through this location increases
with the Stefan number. This explains the decrease in the waviness with increase in the Stefan
number.
It can be seen from figures 57 that the classical enthalpy method seems to constantly under
predict the temperature for higher Stefan numbers. This can be explained as follows. Let a ,
b denote the temperature of a control volume k predicted by the classical enthalpy method at
the start ( a ) and the end ( b ) of a wavy cycle, respectively (figure 5a). This wavy cycle will
correspond to the melting of a control volume i + 1 (figure 5b). a and b will be roughly equal
to the steady state temperatures reached by the control volume k when centroids of the control
volumes i and i + 1 are set at the melting temperature. For a duration a < < *, where * is
the time at which the interface is at the centroid of the control volume at i + 1, the analytical
results at control volume k will be lower than b and vice versa for * < < b . Further, if **

Modified enthalpy method

1271

0.8
Analytical
Classical enthalpy
Modified enthalpy

0.6

X*=0.25

0.4

b
X*=0.75

0.2

X*=0.5

(a)

(b)
Figure 5. (a) Comparison of the temperature at the mid point of the domain as a function of time for
Stefan number of 0.1. (b) The control volumes with interface for the explanation of waviness.

denotes the time taken for the control volume k to reach the pseudo steady state corresponding to
melting temperature being imposed at the centroid of the control volume i + 1, for ** < < b
the temperature of the control volume will be b . If ( ** a ) is less than ( * a ), for some
duration the temperature obtained from the enthalpy method will be over predicted. If however

1272

Niranjan N Gudibande and Kannan N Iyer

0.8
Analytical
Classical enthalpy
Modified enthalpy

0.6
X*=0.25

0.4
X*=0.50

0.2
X*=0.75

0.2

0.4

0.6

0.8

Figure 6. Comparison of the temperature at the mid point of the domain as a function of time for Stefan
number of 0.5.

0.8

Analytical
Classical enthalpy
Modified enthalpy

0.6

X*=0.25

0.4
X*=0.5

0.2
X*=0.75

0.1

0.2

0.3

0.4

0.5

Figure 7. Comparison of the temperature at the mid point of the domain as a function of time for Stefan
number of 1.0.

Modified enthalpy method

1273

( ** a ) is sufficiently greater than ( * a ) the temperature obtained from the enthalpy

method will under predict for the entire duration a << b .
For very low Stefan numbers (such as St = 0.1), the velocity for the solidliquid interface
is lower. This implies that the time taken to reach steady state in a control volume k is significantly lower than the time taken for the control volume i + 1 to undergo phase change, i.e.,
( ** a )<<( b a ). Thus, for some duration, the enthalpy method over-predicts the temperature (figure 5). However, for larger Stefan numbers (such as St = 1) the time taken to reach
steady state is larger than the time taken for the control volume to undergo phase change. This
may be inferred from figure 7 where the curvilinear profile does not reach an asymptotic straight
line as seen in figure 5. This explains the reason behind the under prediction of the temperature
by the enthalpy method in comparison to the analytical solution for the position X * = 0.25 in
figure 7.
It may be noted that in 2 or 3 dimensions when the interface is not oriented along a coordinate axis the waviness in the enthalpy method will be considerably reduced. In such cases,
the control volumes which are undergoing phase change melt at different rates thereby avoiding
synchronisation of the imposition of the artificial boundary condition on the centroids of all the
control volumes undergoing phase change. This will reduce waviness. To demonstrate the reduction in the waviness of the classical enthalpy method and to further benchmark the cell splitting
method, 2-d conduction melting problems are simulated. Since no true 2-d analytical solutions
are available, two pseudo 2-d problems were created from 1-d exact solutions. The first problem
consists of a 2-dimensional domain embedded in the semi-infinite slab melting problem (considered previously) as shown in figure 8. The Cartesian domain has a non-dimensional width
and length of 1. The temperatures on the boundary of the Cartesian domain are calculated from
the analytical solution and are imposed as the boundary conditions. The solidliquid interface in
such a situation will be oriented at an angle of with respect to the x axis of the 2-d domain.
The domain is discretised into 10 10 control volumes. The temperature of the centre of the 2-d
domain is plotted as a function of time in figure 9 for 2 different orientations of the solidliquid
interface.
To benchmark the present method when the solidliquid interface is curved, melting of an
infinite cylinder due to a line heat source is considered (figure 10). The problem consists of an

Figure 8. The 2-d domain chosen for melting with straight interface.

Niranjan N Gudibande and Kannan N Iyer

Analytical
Modified enthalpy
Classical enthalpy

0.4

0.2

An
gle
1

0.3

An
gle
6

1274

0.1

0.5

1.5

Figure 9. Temperature time plot at the centre of a two-dimensional domain for various angle of
orientation of the interface. Stefan number =1.0.

infinite cylinder initially solid at its melting point. A line heat source aligned along the z direction
is switched on at time t > 0. The solidliquid interface which will be an infinite cylinder with axis
along the z direction will expand as time progresses. The governing equation for the situation
will be


T
1
T
=
r
0 < r < R(t),
(43)
t
r r
r
where R is the location of the solidliquid interface. At the origin


T
= Q,
Lim 2r k
r 0
r

Figure 10. The 2-d domain chosen for melting with circular interface.

(44)

Modified enthalpy method

1275

where Q is the strength of the line heat source. At the solidliquid interface
T = TM ,

(45)

d R(t)
T
= H M
.
r
dt

(46)

By introducing the following variables,

r =

r
,
LS

R =

R
,
LS

t
,
L 2S

(T TM )
,
Tr e f

The governing equation in non-dimensional form is




0 < r < R ( ).
= r

r r
r
At the origin

Tr e f =

Q
k

(47)

(48)




Lim 2r
= 1.
r 0
r

(49)

= 0,

(50)

At the solidliquid interface

1
=
r
St

d R ( )
d

(51)

where

C P Tref
.
(52)
H M
The analytical solution for this situation was given by Paterson (1952) which in non-dimensional
form is

 2 
1
r
2
Ei( ) Ei
0 < r < R( ).
=
(53)
4
4
St =

The position of the interface is given by

1

R ( ) = 2 2 ,

(54)

where is obtained by solving the equation

1 2
2
= .
e
4
St

(55)

In the above equations the Ei is the exponential integral defined as

x
Ei (x) =

et
dt.
t

(56)

A 2-d Cartesian domain which has a unit non-dimensional width and length is carved out as
shown in figure 10. For this domain, the temperature at the boundaries of the domain is calculated
from the analytical solutions and is imposed as boundary conditions.

1276

Niranjan N Gudibande and Kannan N Iyer

The domain is discretised into 10 10 control volumes. The temperature of the centre of
the 2-d domain is plotted as a function of time in figure 11. It can be seen that the temperature
predicted by the cell-splitting method shows no waviness and matches well with the analytical
solutions. It can also be seen that the waviness in the enthalpy method has been reduced.
To demonstrate the applicability of the method to the simulation of melting with natural convection, two bench mark problems were chosen, the first involved melting of gallium a low
Prandtl substance (Pr = 0.021) and second involved melting of n-octadecane a high Prandtl
substance (Pr = 57). Both these simulation involved melting in a rectangular container, where
the substance was initially solid at its melting point. At time t = 0 the left wall was raised to
a temperature TW > TM . The top and the bottom walls are insulated. The density difference
between the hot and the cold fluids causes convection currents in the fluid. These convective
currents alter the temperature distribution by creating higher temperatures near the top of the
cavity. This increases the melting rate at the top and causes a distortion in the shape of the solid
liquid interface. In order to model the effect of density difference, Boussinesq approximation
was used. Figure 12a shows the solidliquid interface at various times obtained from the present
method and the enthalpy porosity method of Brent et al (1988) with a mushy zone constant of
105 and 1010 for the melting of a gallium slab. It can be seen that although the deviation between
the solidliquid interface obtained from enthalpy-porosity method and the present cell-splitting
method are small, the solidliquid interface predicted by the present method seems to be sandwiched between the enthalpy-porosity method for mushy zone constants of 105 and 1010 (in later
times). Both the simulations were performed on a 80 80 grid. Thus, by tuning these mushy
zone constants better agreements between the present method and the enthalpy-porosity method
can be obtained. It may be emphasized that although the predictions of the enthalpy-porosity
method seems less sensitive to the mushy zone constants in this situation, significant sensitivity to mushy zone constants has been reported in the melting of a vertical cylinder by Shmueli

0.08
Analytical
Modified enthalpy
Classical enthalpy

0.06

0.04

0.02
1

Figure 11. Temperature time plot at the centre of the two-dimensional domain for Stefan number =1.0.

Modified enthalpy method

0.6

1277

2 min

6 min

17 min

10 min

0.5

Y/L

0.4
0.3
0.2
Present method
Enthalpy porosity1010
Enthalpy porosity105

0.1
0.2

0.4

0.6

0.8

X/L
(a)
+
+
+
+
+
+
+
+
+
+

0.6
0.5

Y/L

0.4
0.3

0.2

+
+

0.1

Enthalpy-porosity
Present method
2 mins
6 mins
10 mins
17 mins

0.2

0.4

0.6

0.8

X/L
(b)
Figure 12. Comparison of the numerical prediction with the (a) enthalpy-porosity model of Brent et al
(1988) for mushy zone constants of 105 and 1010. (b) Enthalpy-porosity model of Brent et al (1988)
for mushy zone constants of 105 and experimental results of Gua & Viskanta (1986) for an aspect ratio
of 0.714.

1278

Niranjan N Gudibande and Kannan N Iyer

1

0.8

y/h

0.6

0.4
Present method
Fo=0.033
Fo=0.065
Fo=0.097
Fo=0.13

0.2

0.5

1.5

x/h
Figure 13. Comparison of the numerical prediction with the experimental results of Okada (1984).
Ra = 1.39 106 , St = 0.192.

et al (2010). The present method eliminates the need for these mushy zone constant. Figure 12b
shows the predictions of the solidliquid interface present method with the experimental results
of Gua & Viskanta (1986). It can be seen that the results agree reasonably well. Figure 13 shows
the solidliquid interface for the melting of n-octadecane. The numerical results carried out on
a 80 80 grid compare well with experimental results of Okada (1984).

7. Conclusions
In this study, a cell splitting method based on enthalpy method for simulation of melting and
solidification is presented. The cell splitting method is demonstrated to eliminate waviness
present in the classical enthalpy method. The method also eliminates the need for any mushy
zone constants for simulating melting with natural convection. Thus for pure materials, the
method has eliminated the need for empirical constants. The method has been successfully
benchmarked against both one and two-dimensional melting problems with and without natural
convection.

Appendix A
In this section, we prove that the interface constructed is second order accurate. Let (x, y) = 0
denote the equation of the actual interface. Consider the situation shown in figure 14. We assume
that in the neighbourhood of x = 0 the curve (x, y) can be written as, y = f (x). Let y =
g(x) = mx + c denote the straight line approximation to this curve.

Modified enthalpy method

1279

Figure 14. The error analysis for interface construction method.

The error e will be defined as the distance between the points A and B in figure 14. The error
is thus
e = |g(0) f (0)| = |c f (0)| .
(57)
The approximate interface g(x) will be chosen by the present interface construction method
such that ratio in which it divides the control volumes P and E are same as that by f (x).
The above mentioned conditions will hold true only if f (x) and g(x) intersect each other at
least once in P and E.
Let xa and xb denote these points of intersection such that 0 < xb < h and h < xa < 0,
where h denotes the dimension of the control volume in the x direction. By Taylor series
f (x) = f (0) +
At the points xa and xb

f  (0)
f  (0) 2
x+
x + ......
1!
2!

f (xa ) = g(xa )
f (xb ) = g(xb ).

Thus,
mxa + c = f (0) +

(59)

f  (0)
f  (0) 2
xa +
x + ......
1!
2! a

f  (0)
f  (0) 2
xb +
x + ......
1!
2! b
Multiplying equation (60) by xb and equation (61) by xa and subtracting,
mxb + c = f (0) +

(xb xa ) (c f (0)) =

(58)



f  (0)
f n (0)
xa xb (xb xa ) + ...... +
xa xb xbn xan .
2!
n!

(60)
(61)

(62)

The above equation can be simplified into,

n2


f n (0)
f  (0)
xar xbn2r .
xa xb + ...... +
xa xb
c f (0) =
2!
n!

(63)

Taking the modulus of the function


n2


 
 f n (0)



 f (0)


n2r
r
|c f (0)| 
xa xb
xa xb  + ...... + 
xa xb
 + ......
 n!

2!
0

(64)

1280

Niranjan N Gudibande and Kannan N Iyer

Since
|xa | h

and

|xb | h,

(65)


 n
  
 f (0)  n
 f (0)  2
 h + ...... .



|c f (0)| 
h + ...... + 
2! 
(n 1)! 

(66)

Thus from equations (57) and (66)

e = O(h 2 ).

(67)

Appendix B
In this section, the form taken by the present cell splitting method for one-dimension situation
is presented. The significant differences between the present scheme and the explicit enthalpy
scheme of Tacke (1985) which eliminate waviness are also highlighted. Consider a stencil of
control volume as show in figure 15. The solidliquid interface is assumed to be present in the
control volume P while control volume E is liquid and W is solid. The centroids of the
liquid and the solid control volumes are denoted by X L and X S , respectively. The liquid fractions
at the east and the west face of the control volume are denoted by f e and f w , respectively. For
the sake of clarity, the temperatures in the solid and the liquid domain are distinguished by the
superscripts S and L, respectively.
The differential form of the energy equation for the liquid domain obtained by suitably
substituting for in equation (11) is





d
.
C P T L TM d V = kT L .nd
A + kT L .nd
A.
(68)
dt
V

AI

Equation (68) can be discretised for the control volume shown in figure 15 as
L L
TW + a IL TM + b,
a PL TPL = a EL TEL + aW

(69)

where
a EL = f e

k
,
x EL x PL


aI =

k
x I x PL

L
aW
= fw

k
,
L
x PL x W

0 < fp < 1
f P = 0, f P = 1

Figure 15. The stencil for the control volume discretisation.

(70)

(71)

Modified enthalpy method

1281

L
a PL = a EL + aW
+ a IL +

b=

C P f P x P
t


C P f PO x P O L
C P f P x P
TM .
T
TM +
t
t

(72)

(73)

In the above equations the superscript O denotes the values in the previous time steps and
x P denotes the width of the control volume P in figure 15. In the situation shown in figure 15
where the liquid is not present on the east face a EL = 0. A similar equation can be derived to
obtain the temperature in the solid domain.
The enthalpy form of the energy equation (12) is



d
hd V + hV.nd
A = kT.nd
A.
(74)
dt
V

AC

AC

which can be written in a discrete form of as








hP hO
L
S
P
TWL TPL + aW
TWS TPS , (75)
= a EL TEL TPL + a ES TES TPS + aW
t

S the coefficients obtained by solving the energy equation in the solid domain.
where a ES and aW
From the enthalpy equation the liquid fraction can be recovered from,




(76)
h P = C P f P TPL TM + C P (1 f P ) TPS TP + f P H M ,

which asserts that the enthalpy of a control volume is composed of its sensible and latent heat
parts.
Tacke (1985) had proposed an explicit enthalpy scheme which eliminates waviness in the
enthalpy method. The significant difference between the schemes of the present cell splitting
method with that of Tacke (1985) are
(a) In the method by Tacke (1985) the diffusive flux on the west face of the control volume for
e.g., was discretised as
Q W = k

TM TW
,
X I XW

(77)

while in the present method the diffusive flux is discretised as

Q W = k

TPL TW
X PL X W

(78)

(b) In the method by Tacke (1985) the temperature of the solid and liquid domains which are
necessary to calculate the sensible heat in the equation (76) are obtained by linear interpolation of the temperatures TM , TE and TM , TW . While in the present method these temperature
are obtained by solving the energy equation (69).

1282

Niranjan N Gudibande and Kannan N Iyer

Nomenclature
A
AI
AC
a
b
CP
C
D
d
e
f
f (x)
FP X FPY
g
g(x)
H
h
k
LS
m
N
n
P
P
Q
qD
r
r
R
R
S
S( )
SS
SS P
s(t)
T
Tr e f
TM
TW
t

Area of a face of a control volume, m2

Area of the interface part of a control volume, m2
Area of the Cartesian part of a control volume, m2
Coefficients appearing in a transport equation
Source term in discretised equation
Specific heat, J/kg-K
y intercept of the constructed interface, m
Net diffusion into a control volume
Coefficients appearing in a diffusion equation
Error between the actual and the predicted
location of the interface
Liquid fraction
Equation of the interface
Force due to pressure in x and y direction, N
Acceleration due to gravity, m/s2
Equation of the constructed interface
Enthalpy, J/kg
Width of the control volume in appendix A, m
Thermal conductivity, W/m-K
A characteristic length scale, m
Slope of the constructed interface
Direction normal to the solidliquid interface, m
Unit vector normal to a face
Pressure, N/m2
Pressure correction, N/m2
Strength of a line heat source, W/m
Net diffusion into a control volume
Radial coordinate, m
Non-dimensional radial coordinate
Radial location of the solidliquid interface, m
Radial location of the solidliquid interface
in non-dimensional form
Volumetric source terms in a transport equation
Non-dimensional position of the solidliquid
interface along the x direction
Source term due to secondary gradient in
transport equation
Source term due to secondary gradient in
pressure, m/s
Position of the solidliquid interface along
the x direction, m
Temperature, K
A reference temperature, K
Melting point temperature, K
Temperature at the wall, K
Time, s

Modified enthalpy method

(u,v)
u,
v
V
V
Vg
X
(x,y)
(xa ,xb )
H M
V

Velocity coordinates, m/s

Tilde velocity used in momentum interpolation, m/s
Volume of a control volume, m3
Velocity vector, m/s
Velocity of the grid, m/s
Non-dimensional coordinate along x axis
Coordinate axis, m
x location of the points at which the
constructed interface intersects
the actual interface
Latent heat of melting, J/kg
Volume of a control volume, m3

Non-dimensional numbers
Pr
Ra
St

Prandtl number
Rayleigh number
Stefan number

Greek Symbols

(x, y)

a , b

a , b
*
**



Thermal diffusivity, m2 /s
Volumetric expansion coefficient, 1/K
Equation of the interface
Non-dimensional temperature
Non-dimensional temperatures at the start
and end of a wavy cycle respectively
A constant appearing in equation (55)
Dynamic viscosity, Pa-s
Kinematic viscosity, m2/s
A coordinate along the line joining the
centroids of the control volumes, m
Density, kg/m3
Non-dimensional time
Non-dimensional times at the start and end
of a wavy cycle respectively
Non-dimensional time at which the interface
is at the centroid of the control volume
Non-dimensional time when the one dimension
melting system reaches pseudo steady
state as predicted by classical enthalpy method
A generic transport variable
The domain

Subscripts
E,W ,N ,S
e,w,n,s

Tags used to address the neighbours of chosen control volume

Tags used to address the faces of chosen control volume

1283

1284
f
I
nb
P
n
x,y,
S
L

Niranjan N Gudibande and Kannan N Iyer

Face
Interface
Neighbours
A tag used to address a chosen control volume
Normal
Partial derivatives along the corresponding directions
Solid
Liquid

Subscripts
S
L

Solid
Liquid

Abbreviations
CV
SL

Control Volume
Solidliquid

References
Brent A D, Voller V R and Reid K J 1988 Enthalpy-porosity technique for modelling convection-diffusion
phase change: Application to the melting of a pure metal. Numer. Heat Transfer 13(3): 297318
Cao Y and Faghri A 1990 Numerical analysis of phase-change problems including natural convection.
J. Heat Trans. 112(3): 812816
Cerimele M M, Mansutti D and Pistella F 2002 Numerical modelling of liquid/solid phase transitions
analysis of a gallium melting test. Comput. Fluids 31(47): 437451
Crank J 1984 Free and moving boundary problems (Oxford: Oxford University Press)
Date A W 1994 A novel enthalpy formulation for multidimensional solidification and melting of a pure
substance. Sadhana 19(5): 833850
Donea J, Giuliani S and Halleux J P 1982 An arbitrary Lagrangian-Eulerian finite element method for
transient dynamic fluid-structure interactions. Comput. Method Appl. M. 33(13): 689723
Gartling D K 1980 Finite element analysis of convective heat transfer problems with change of phase.
Computer Methods in Fluids, Morgan K, Taylor C, Brebbia C A (eds), Pentech, London, 257284
Gua C and Viskanta R 1986 Melting and solidification of a pure metal on a vertical wall. J. Heat Trans.
108(1): 174181
Kim S and Anghaie S 2001 An effective conduction length model in the enthalpy formulation for the Stefan
problem. Int. Commun. Heat Mass 28(6): 733741
Kumar V, Durst F and Ray S 2006 Modeling moving-boundary problems of solidification and melting
adopting an arbitrary Lagrangian-Eulerian approach. Numer. Heat Tr. B- Fund. 49(4): 299331
Lacroix M 1989 Computation of heat transfer during melting of a pure substance from an isothermal wall.
Numer. Heat Tr. B- Fund. 15(2): 191210
Li C Y, Garimella S V and Simpson J E 2003 Fixed-grid front-tracking algorithm for solidification
problems, part I: Method and validation. Numer. Heat Tr. B- Fund. 43(2): 117141
Mathur S R and Murthy J Y 1997 A pressure-based method for unstructured meshes. Numer. Heat Tr. BFund. 31(2): 195215

Modified enthalpy method

1285

Niranjan G N and Iyer K 2011 Modified enthalpy method for the simulation of melting and solidification.
Proc. 21st National and 10th ISHMT-ASME Heat and Mass Transfer Conf., Chennai, India, December
2730. Paper number ISHMT_16_017
Okada M 1984 Analysis of heat transfer during melting from a vertical wall. Int. J. Heat Mass Tran. 27(11):
20572066
Patankar S V 1980 Numerical Heat Transfer and Fluid Flow (Washington, DC: Taylor and Francis)
Paterson S 1952 Propagation of a boundary of fusion. Glasgow Math. J. 1(1): 4247
Pilliod Jr J E and Puckett E G 2004 Second-order accurate volume-of-fluid algorithms for tracking material
interfaces. J. Comput. Phys. 199(2): 465502
Price P H and Slack M R 1954 The effect of latent heat on numerical solutions of the heat flow equation.
Br. J. Appl. Phys. 5(11): 285287
Rhie C M and Chow W L 1983 Numerical study of the turbulent flow past an airfoil with trailing edge
separation. AIAA J. 21(11): 15251532
Rider W J and Kothe D B 1998 Reconstructing volume tracking. J. Comput. Phys. 141(2): 112152
Shmueli H, Ziskind G and Letan R 2010 Melting in a vertical cylindrical tube: Numerical investigation and
comparison with experiments. Int. J. Heat Mass Trans. 53(1920): 40824091
Tacke K 1985 Discretization of the explicit enthalpy method for planar phase change. Int. J. Numer. Meth.
Eng. 21(3): 543554
Udaykumar H S, Shyy W and Rao M M 1996 ELAFINT: A mixed Eulerian-Lagrangian method for fluid
flows with complex and moving boundaries. Int. J. Numer. Meth. Fl. 22(8): 691712
Voller V and Cross M 1981 Accurate solutions of moving boundary problems using the enthalpy method.
Int. J. Heat Mass Tran. 24(3): 545556
Voller V R, Cross M and Markatos N C 1987 Enthalpy method for convection/diffusion phase change.
Int. J. Numer. Meth. Eng. 24(1): 271284