10 1016-j Rser 2014 02 032

Mathematical solutions and numerical models employed
for the investigations of PCMs' phase transformations

Shuli Liu
a,n
, Yongcai Li
a,b
, Yaqin Zhang
a
a
Department of the Civil Engineering, Architecture and Buildings, Faculty of Engineering and Computing, Coventry University, CV1 5FB, UK
b
Faculty of the Urban Construction & Environment Engineering, Chongqing University, China
a r t i c l e i n f o
Article history:
Received 5 March 2013
Received in revised form
9 January 2014
Accepted 22 February 2014
Available online 15 March 2014
Keywords:
Stephan problem
Neumann problem
Heat transfer mechanisms
Mathematical solutions
Numerical models
a b s t r a c t
Latent heat thermal energy storage (LHTES) has recently attracted increasing interest related to the thermal
applications, and has been studied by researchers using theoretical and numerical approaches. The heat
transfer mechanisms during the charging and discharging periods of the phase change materials (PCMs) and
two basic problems for phase transformations have been discussed in this paper. 1D, 2D and 3D popular
mathematical solutions based on the heat transfer mechanisms of conduction and/or convection have been
analyzed. Then, various numerical models for encapsulated PCMs in terms of macro-encapsulated PCM and
micro-encapsulated PCM were investigated. The numerical simulation of heat transfer problems in phase
change processes is complicated and the achieved results are approximate according to a number of conditions.
The accuracy of numerical solution depends on the assumptions that are made by the authors. Therefore this
review will enable the researchers to have an overall view of the mathematical and numerical methods used
for PCM's phase transformations. It offers the researcher a guideline of selecting the appropriate theoretical
solutions according to their researches' purposes.
& 2014 Elsevier Ltd. All rights reserved.
Contents
1. Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 660
2. Heat transfer mechanisms during the phase transformations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 661
2.1. Conduction and convection heat transfer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 661
2.2. Stefan problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 661
2.3. Neumann problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 662
2.4. Other possibility mechanisms in phase change process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 662
3. Mathematical solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 663
3.1. Conduction acting as the only heat transfer mechanism. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 663
3.1.1. Fixed grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 663
3.1.2. Adaptive grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 665
3.2. Conduction and convection heat transfer mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 665
4. Numerical models for various PCMs capsules and packages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 667
4.1. Macro-encapsulated PCMs numerical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 667
4.1.1. Rectangular encapsulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 667
4.1.2. Cylindrical containment models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 669
4.1.3. Spherical containment models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 670
4.2. Microencapsulated PCM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 670
4.2.1. Incorporation with building components . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 670
4.2.2. Microencapsulated phase change slurry (MPCS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 671
4.3. Summary of numerical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 672
5. Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 672
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 672
Contents lists available at ScienceDirect
journal homepage: www.elsevier.com/locate/rser
Renewable and Sustainable Energy Reviews
http://dx.doi.org/10.1016/j.rser.2014.02.032
1364-0321 & 2014 Elsevier Ltd. All rights reserved.
n
Corresponding author. Tel.: 44 24 7765 7822.
E-mail addresses: [email protected], [email protected] (S. Liu).
Renewable and Sustainable Energy Reviews 33 (2014) 659674

1. Introduction
Due to the attractive features of latent heat storage, phase
change materials (PCMs) are mainly used to store energy at a xed
temperature (melting/solidication temperature) with high
energy density [1]. PCMs have been applied in various systems
and aspects such as energy storage systems, free ventilation, free
heating and cooling for buildings, spacecraft, food, medicine
conservation, smoothing the peaks of exothermic temperature in
chemical reactions etc. PCMs are a better choice comparing to the
sensible heat storage in applications, because of its nearly iso-
thermal storing mechanism and high storage density. PCMs absorb
abundant energy through the phase transition and release the
stored energy. Table 1 compares the typical properties of different
thermal storage materials tested at the room temperature.
It indicates that PCMs can save 92.8% of mass and up to 90% space
to store the same amount of thermal energy comparing to the
sensible thermal materials such like concrete and water.
Therefore, PCMs attract the researchers' interesting in practically
individual and incorporation applications in different engineering
elds. The possibility, feasibility, thermal performance and economic
analysis of using PCMs call a series of theoretical and experiential
investigations, the mainly two methods used to research the thermo-
physical parameters of PCMs and interrelated systems.
The experimental approaches offer a better indication of the
actual PCM behavior and performance in comparison to theore-
tical analysis, as the experimental tests can present the PCMs'
behaviors more directly, visibly and credibly without any pre-set
assumptions however, the experiments are unachievable in some
cases, such as the large scale or unsteady around environment, so
there are still a few points need to be considered:
Time and cost consuming will be higher than a theoretical

approach, hence, the budget and processes needs to be estab-
lished and scheduled.
The scales of the experiment test rig have to be considered, and

then suitable laboratory location and space needs to be
determined to house the rig.
Test rig needs to be constructed and operated properly.
Proper environment parameters have to be simulated and

adjusted to imitate the practical environment.
Relevant parameters need to be monitored, measuring appara-

tuses need to be calibrated, and the failed data need to be
eliminated.
Except these agonizing experimental matters, there are still
some unavoidable testing errors. However the theoretical methods
can avoid all these weakness and predicate the PCMs performance
conveniently. The major advantage of the theoretical/numerical
approaches is that various conditions can be carried out by
changing the variables in a numerical model. Therefore, more
and more investigators prefer to study the phase change problems
by mathematical solutions and numerical simulations. There are
only eight review papers on PCMs that involves the aspect of
mathematical solutions and/or numerical modeling of latent heat
thermal energy storage (LHTES) have been published during
20022012. Table 2 summaries these eight relative review papers
and comments their mathematical and/or numerical aspects of
LHTES. However, few logistic and comprehensive reviews on the
mathematical solutions and numerical modeling for melting and
solidication processes of PCMs were found in published litera-
tures even though there are many individual publications in phase
change problems, particularly in heat transfer mechanisms and
simulations. Due to the complexity of phase transformations,
Nomenclature
f melt fraction
c specic heat at constant pressure
H total volumetric enthalpy (J/m
3
)
h average heat transfer coefcient (J/kg)
k thermal conductivity (W/mK)
L latent heat (J/kg)
St Stefan number
t time (s)
T temperature (K)
T
m
melting temperature (K)
T
s
initial temperature of solid PCM (K)
T
0
constant temperature imposed on x0 (K)
P present nodal point
t position of liquidsolid interface
x; y Cartesian coordinates
Greek symbols
difference
p thermal diffusivity of PCM (m
2
/s)
dimensionless number used in derivations as a tem-
porary substitution
dimensionless number in solution to Neumann pro-
blem, liquid fraction
density (kg/m
3
)
Subscripts
l liquid
s solid
i ith spatial step
j jth time step
ef f effective
Table 1
Properties comparison of different thermal storage materials.
Property Unit Materials
Concrete (sand
and gravel)
Water Organic
PCM
Inorganic
PCM
Density [kg=m
3
] 2240 1000 800 1600
Specic heat
capacity
[kJ=kg K] 1.1 4.2 2.0 2.0
Latent heat [kJ=kg K] 190 230
Storage mass for
10
6
kJ, avg
[kg] 60,000 16,000 5300 4350
Storage volume for
10
6
kJ, avg
[m
3
] 26.8 16 6.6 2.7
Relative storage
mass
13.79 4 1.25 1.0
Relative storage
volume
10 6 2.5 1.0
1. Storage mass and volume are calculated for storing 106 kJ energy with a
temperature rise of 15 K for concrete (sand and gravel) and water.
2. Relative mass and volume are based on the inorganic phase change materials.
S. Liu et al. / Renewable and Sustainable Energy Reviews 33 (2014) 659674 660

a good understanding of the heat transfer mechanisms, phase
change characteristics and the differences among various mathe-
matical and numerical simulation methods is required before
researchers starting their theoretical study. The aim of this paper
is to comprehensively study the mathematical and simulation
modeling applied for LHTES. Firstly it discusses the issues of the
heat transfer mechanisms during the charging and discharging
processes and the Stephan problem and Neumann problem.
Secondly, it sums the strong and weak aspects of various mathe-
matical solutions when conduction and convection heat transfer
occurring inside the PCMs by the considering xed grid method
and the adaptive grid method. Finally, numerical modeling of
various PCM packages' geometries in terms of macro-encapsulated
PCMs and micro-encapsulated PCMs are studied.
2. Heat transfer mechanisms during the phase
transformations
2.1. Conduction and convection heat transfer
During the charging and discharging processes, the possible
heat transfer mechanisms are conduction and convection. How-
ever, the issue of which heat transfers mechanism takes the main
role in different stages of phase transformation has been argued
for decades. When PCMs are used to store or release thermal
energy, conduction is usually believed playing the most important
role on the heat transfer during the phase transformation pro-
cesses [24]. As early as 1831, Lam and Clapeyron have conducted
the rst study on phase change problem by only considering the
pure conduction. However, some researchers persist that natural
convection is a more important mechanism in the phase change
process especially in the melt region. Lazaridis [100] proposed a
study of the relative importance of conduction and convection in
1970. A pioneer study performed by Sparrow et al. in 1977 [41],
they concluded in their study that natural convection could not
be ignored in the analysis of phase change problems. In 1994,
Hasan [5] concluded that the convection heat transfer active an
important role in the melting process, and a simplied model by
only considering the conduction heat transfer does not describe
the melting process properly. Later, Lacroix et al. [6] reported
the similar ndings in their research that natural convection is
the main heat transfer mechanism during the melting process.
In 1999, Velraj [7] obtained the same conclusion in their research
and reported that during the melting process natural convection
occurs in the melt layer which results in the heat transfer rate
increase comparing to the solidication process. Buddhi et al. [8]
proposed an explanation for this phenomenon that the density
differences between the solid and liquid PCM induce the buoy-
ancy, which causes the convective motions in the liquid phase.
However, Zhang and Yi [9] believed that with the solid PCM
melting into liquid, the PCM volume keeps increasing, which
results in the convection becoming the predominant heat transfer
mechanism rather than conduction. Sari et al. [10] found that the
heat transferred from a heat exchanger to a PCM of stearic acid is
largely inuenced by natural convection at the melting layer, in
addition to conduction and forced convection heat transfer.
For the melting process, the PCM changes its phase from solid
to a mushy state, and then liquid, which is reversible during the
solidication process. Hence there are six stages to nish a
charging and discharging cycle. Therefore it is possible in certain
stages of the phase transformation process there are more than
one kind of heat transfer mechanisms acting at work, but how to
weight the percentages of conduction and convection heat transfer
in each stage has been the main challenge for the researchers
currently. This paper will review the most common research
methods used for the conduction and convection heat transfer
inside the PCM packages.
2.2. Stefan problem
Moving boundary problem named Stefan problem is another
issue to develop a numerical modeling of PCMs [11,12]. The
simplest of the Stefan problems is the one-phase Stefan problem
since only one-phase involved. The term of one-phase designates
only the liquid phases active in the transformation and the solid
phase stay at its melting temperature. Stefan's solution with
constant thermophysical properties shows that the rate of melting
or solidication in a semi-innite region is governed by a dimen-
sionless number, known as the Stefan number (St),
St
C
l
T
l
T
m
L
1
where C
l
is the heat capacity of the liquid PCMs, L is the latent
heat of fusion, and T
l
and T
m
are the surrounding and melting
temperatures, respectively.
Table 2
Some review papers related to numerical study of phase change problems.
Reference Published
year
Journal Comment
Zalba et al. [1] 2003 Applied thermal
engineering
The review was divided into three parts: materials, heat transfer and applications. Heat transfer
was considered mainly from a theoretical point of view, considering different simulation techniques.
Verma et al. [2] 2008 Renewable and sustainable
energy reviews
Mathematical models of a LHTES were reviewed to optimize the material selection and to assist
in the optimal designing of the systems. Some important characteristics of different models and their
assumptions used were presented and discussed as well.
Jegadheeswaran
and Pohekar [3]
2009 Renewable and sustainable
energy reviews
Various thermal conductivity enhancement techniques reviewed in this paper, and the issues related
to mathematical modeling of LHTES with enhancement techniques are also discussed.
Zhu et al. [4] 2009 Energy conversion and
management
This paper presented an overview on dynamic characteristics and energy performance of buildings
employing PCMs by three research methods used, i.e., simulation, experiment, combined
simulation and experiment
Agyenim et al. [5] 2010 Renewable and sustainable
energy reviews
This paper provided the formulation of the phase change problem. In terms of problem formulation,
it was concluded that the common approach has been the use of enthalpy formulation.
Kuznik [6] 2011 Renewable and sustainable
energy reviews
The review was the rst comprehensive review of the integration of phase change materials in building
walls. Various numerical studies concerning the integration were summarized.
Dutil et al. [7] 2011 Renewable and sustainable
energy reviews
The review presented models based on the rst lawand on the second lawof thermodynamics. This paper
tried to enable one to start his/her research with an exhaustive overview of the subject.
Zhou [8] 2012 Applied energy This paper reviewed previous works on latent thermal energy storage in building applications, including
PCMs, current building applications and their thermal performance analyses, as well as numerical
simulation of buildings with PCMs.

The most published solutions are to solve the one-dimensional
cases subjected to an innite or semi-innite region with simple
initial and boundary conditions. Whereas with the heat transfer
continuing, the interface boundary is constantly moving as the
liquid and solid phases shrinking and growing, which disable the
prediction of the boundary location [13]. Because of that the solid
liquid interface is not xed, but moving with time, the heat
transfer mechanisms during a PCM phase transformation process
are complex. Therefore the phase change transition is difcult to
analyze owing to the three reasons: the solidliquid interface is
moving; the interface location is non linear; it consists of thermal
conduction and natural convection heat transfer mechanisms.
Since these three factors, the non-linearity of the governing
equations is introduced to the moving boundary, and the precise
analytical solutions are only possible for a limited number of
scenarios [14]. This view is shared by Krkl et al. [15] who
pointed out that the mathematical models have been proposed in
each research only applied to very specic boundary conditions,
hence are not feasible to more complex practical applications. This
Stefan problem is further complicated by the fact that most of the
methods previously proposed by researchers only involve one
moving boundary, whereas it actually consist of more than one
interface location [16].
2.3. Neumann problem
The Stefan problem was extended to the two-phase problem,
the so-called Neumann problem which is more realistic [17].
In Neumann problem, the initial state of the PCM is assumed to
be solid, during the melting process, its initial temperature does
not equal to the phase change temperature, and the melting
temperature does not maintain at a constant value. If the melting
happens in a semi-innite slab (0oxo1), the solid PCM is
initially at a uniform temperature T
s
oT
m
, and a constant
temperature T
0
is imposed on the slab surface x 0, with the
assumptions of constant thermo-physical properties of the PCM,
the problem can be mathematically expressed as follows:
Heat conduction in solid or liquid region
C
T
t
k
2
T
x
2
2
where is the density, C is the specic heat, k is the thermal
conductivity, and t and x are the time and space coordinates
respectively.
The heat uxes transferring from the liquid phase to the solid
liquid interface, as well as the latent heat absorbing rate by the
melting PCM, the movement of the solid/liquid interface can be
determined through the following energy balance:
L
x
t
t
k
l
T
l
x
k
s
T
s
x
3
where L is the latent heat of fusion of the PCM and is the density
of liquid PCM.
Initial condition
Tx40; t 0 T
s
oT
m
4
Solidliquid interface temperature
Tx t; t 40 T
m
5
With the following boundary conditions:
T0; t T
0
4T
m
for t 40 6
Tx; t T
s
for x-1; t 40 7
where t is the position of the solidliquid interface (melting
front). Fig. 1 clearly illustrates this two-phase Stefan problem.
Analytical solution to such a problem was obtained by Neu-
mann in term of a similarity variable

t
2

p
l
p 8
The nal Neumann's solution can be written as
Interface position
t 2

p
l
t
p
9
Temperature of the liquid phase
Tx; t T
l
T
l
T
m
erf x=2

p
l
t
p

erf
10
Temperature of the solid phase
Tx; t T
s
T
m
T
s
erf cx=2

p
s
t
p

erf c

p
l
=p
s
p

11
The in Eq. (9)(11) is the solution to the transcendental
equation
St
l
exp
2
erf
St
s
p
s
p
p
l
p
expp
l
2
=p
s
erf c

p
l
=p
s
p

p
12
where
St
l

C
l
T
l
T
m
L
13
St
s

C
s
T
m
T
s
L
14
However, the Neumann's solution is applicable only for moving
boundary problems in the rectangular coordinate system.
2.4. Other possibility mechanisms in phase change process
Furthermore, there are several other issues with the use of a
theoretical approach in the study of PCMs. Alexiades [18] pointed
out that there were many mechanisms involved during a PCM
phase transition, such like a change in volume, density, thermal
conductivity, specic heat capacity, super-cooling, etc. Conse-
quently, accurate reection of the proposed mathematical solution
and numerical model need to consider the dynamic properties of
the phase change process. Another major issue with PCMs is that
they act as self-insulating materials. When PCM solidication
occurs from the top of the heat surface, solid insulating layer will
be developed which moves inward during the whole solidication
process. With the increase in the size and thickness of the solid
layer, the heat transfer rate from the liquid PCM to the heat
exchanger surface decreases until it becomes so small that will not
Fig. 1. Schematic illustration of the two-phase Stefan problem.

be possible to maintain at an acceptable heat transfer rate. The
self-insulating characteristic of the PCM becomes a strong problem
when bulk containerization is used. For example the large cylind-
rical used as containers lled with PCM, the self-insulating issue
can become very serious. Large heat exchange area will accelerate
the formation of the crust and heavily decrease the heat transfer
rate. Radhakrishnan et al. [14] analyzed that the method to lighten
the crust issue is minimizing the depth of the PCM pack. So that
even a large majority of the PCM performs solidication, the
thermal resistance is reduced to allow the removal of the latent
heat released from the solidied PCM. The self-insulating char-
acteristic of the PCM during discharging process needs to be
considered in any mathematical model as this can affect thermal
energy discharging time and rate. Furthermore, the self-insulating
has a direct effect on the phase transition boundary, where the
solid and liquid region meets and co-exists, which is linked back to
the Stephan Problem and Neumann's method again.
Logical heat transfer mechanisms during the charging and
discharging periods of the PCM are the essential elements to
develop any accurate mathematical model. However, as described
there are quite a few issues to properly describe the heat transfer
mechanisms: conduction and natural convention during the phase
change process. This paper reviews most of the researches
employing the mathematical calculation and numerical simulation
published in the articles and reports, but does not nd any
consensuses on which one is the dominant heat transfer mechan-
ism, which one can be ignored and how to combine these two
transfer mechanisms in each modeling. Different researchers have
developed their specic methods to evaluate the thermal energy
transfer during the phase transition process for various situations
[1923].
3. Mathematical solutions
The analyzes of the heat transfer mechanisms in phase change
processes are complex owing to the movable solidliquid bound-
ary, which depends on the latent heat absorbing or releasing rate
at the interface. In order to simplify the simulation of the phase
change problems, assumptions have been made by various
researchers in their studies. In this section, different numerical
simulations in terms of conduction heat transfer and conduction
and/or convection heat transfer are discussed and analyzed.
3.1. Conduction acting as the only heat transfer mechanism
To develop a mathematical model, the dominant heat transfer
mechanism whether conduction or convection, during the char-
ging and discharging processes, need to be analyzed and deter-
mined. Although there are quite a few issues in the heat transfer
mechanisms during the phase transformation, some researchers
simplify this process by only considering conduction in pure sub-
stances. In this case, the mathematical solutions can be classied
as follows: xed grid and adaptive grid.
3.1.1. Fixed grid
Fixed grid method does not require explicit treatment of
conditions on the phase change boundary, and just requires the
use of xed-grid schemes able to cope with strong nonlinearities.
Therefore this method is able to utilize standard solution proce-
dures for the energy equations. They are amenable to physical
interpretation and easy to implement, especially for multidimen-
sional problems and for mufti-interface problems.
In the xed grid method, the heat ow equation is approximated
by nite difference replacements for the derivatives in order to
calculate the values of temperature T
i;j
, at any time t
n
nt, the
moving boundary will be located between two adjacent grid points,
for instance, between mx and m1x at a location of x;
where 0oo1, as illustrated in Fig. 2.
The mathematical solution of the one-phase ice-melting pre-
sents a simple illustration of this method. The instantaneous
temperature at point i; j 1 is calculated from the temperatures
of the previous step on the basis of the following formulation:
T
i;j 1
T
i;j
t
x
2

T
i 1;j
2T
i;j
T
i 1;j
i 0; 1; m1: 15
In terms of three-point Lagrange interpolation instead of
Eq. (14), the temperature at x mx is computed as,
T
m;n1
T
m;n
2t
x
2

1
n1
T
m1;n
n
T
m;n

16
The variation of the location of the moving boundary is
n1

n
t
x
2

n
n1
T
m1;n
n
T
m;n

17
As shown, the numerical solution of this xed grid method is
carried out on a space grid remaining at a xed value throughout
the mathematical calculation. To approximate the Stefan condi-
tions both for the moving boundary and the partial differential
equation at the adjacent grid points, various mathematical solu-
tions have been proposed. Murray and Landis [24] proposed two
ctitious temperatures: one is achieved by quadratic extrapolation
from the temperatures in the solid PCM and the other one is from
the temperatures in the liquid PCM. The melting temperature and
the position of the moving interface are incorporated in the
ctitious temperatures, which are then substituted into a stan-
dard approximate to calculate the temperature near the inter-
face instead of using special formulae. Ciment and Guenther [25]
introduced a method of spatial mesh renement on both sides of
the moving boundary. With this method, Lazaridis [26] applied the
explicit nite difference approximations to solve two-phase soli-
dication problems in both two and three space dimensions.
The major advantage of the xed grid method is that the
mathematical calculation of the phase transition can be achieved
through simple modications of existing heat transfer numerical
methods and/or software. As such, they can be used for modeling of
a variety of complex moving boundary problems. Basu and Date [27]
and Voller et al. [28] reviewed the applications of the xed grid
methods for phase change problem. The two mainly xed-grid
methods that have been proposed are enthalpy-based methods
and the effective heat capacity method. However, the xed grid
method sometimes cannot work efciently when the boundary
Fig. 2. Position of the moving boundary in a xed grid.

moves along a distance larger than a space increment in a time step,
depending on the velocity of the moving boundary. Therefore, the
adaptive grid is developed by researchers and will be discussed in
Section 3.1.2.
3.1.1.1. Enthalpy based method. The existence of nonlinearity of the
liquidsolid interface and the unknown location of the moving
boundary are the two main challenges in simulating the phase
change problems [29]. To deal with these, a new model formu-
lation called the enthalpy method was introduced [30]. In the
enthalpy models, the enthalpy is used as a dependent variable
value along with temperature. By introducing an enthalpy method,
the phase-change problem becomes much easier since the
governing equations are the same for the two phases. The
interface conditions are automatically achieved and a mushy
zone between the two phases is created. This zone avoids sharp
discontinuities that may create some numerical instability. Hence,
the enthalpy mathematical methods are most attractive and
commonly used to handle phase change problems where a
solidliquid mush zone is present between two phases. Hunter
[31] conrmed that the enthalpy method is most suitable for
typical applications, the reason for that is this method does not
require explicit treatment of the conditions on the phase change
boundary [32].
Enthalpy-based method can be illustrated by considering a
one-dimensional heat conduction controlled phase transforma-
tion. For a phase change process, energy conservation can be
expressed in terms of total volumetric enthalpy and temperature.
This relationship between the enthalpy and temperature is usually
assumed to be a step function for isothermal phase change problems.
Fig. 3 shows the enthalpytemperature curves for isothermal phase
change.
For isothermal phase change, the conservation of energy can be
expressed in terms of temperature and the total volumetric
enthalpy [33]
H
t
kT 18
where the total volumetric enthalpy H is the sum of sensible and
latent heat, and can be expressed in terms of temperature of the
PCM as follows:
HT HT
l
f T 19
The rst term on the right side of Eq. (19) accounts for the sensible
heat of the PCM while the second term accounts for the latent heat
of the PCM.
And where
H
Z
T
Tm
cdT 20
To be able to calculate the total enthalpy, the liquid fraction is
given as follows:

0 T oT
m
solid
0; 1 T T
m
mushy
1 T 4T
m
liquid
8
>
<
>
:
21
Integrating the Eqs. (18) and (20), an alternate form for one-
dimensional heat transfer in the PCM was obtained:
H
t

x

H
x

l

f
t
22
where is the thermal diffusivity.
3.1.1.2. Effective heat capacity method (energy based method). The
effective heat capacity method is also a common method used in
modeling phase change problems except for the enthalpy method.
The effective heat capacity method was rst proposed in [34,35]. In
the effective heat capacity method, the latent heat is approximated
by a large heat insensible form over the phase change temperature
interval, (T
2
T
1
) [36]. The effective heat capacity of the PCM (c
ef f
)
is directly proportional to the energy stored and released during the
phase change but inversely proportional to the interval of the
melting or solidication temperature range. During the phase
change the heat capacity of the PCM is
c
ef f

L
T
2
T
1
c
s
23
where T
1
is the temperature at which melting or solidication
begins and T
2
is the temperature at which the PCM is totally melted
or solidied. The following is a denition of the effective heat
capacity for each phase change period.
c
ef f

c
s
T oT
1
solid
L
T
2
T
1
c
s
T
1
rT rT
2
mushy
c
l
T 4T
2
liquid
8
>
<
>
:
24
In terms of the denition of the effective heat capacity, the
energy equation for one dimension can be written as follows:
c
ef f
T
t
k
2
T
x
2
25
Although the enthalpy based method and the effective heat
capacity method achieve much simpler mathematical operations
with reasonable accuracy, it is worth noting that each of the
methods has its inherent disadvantages. Enthalpy method is
difcult to handle super-cooling and temperature oscillation
problems. Whilst the effective heat capacity method is difcult
in resolving the phase change problems when the phase change
temperature range is small, and is not applicable for the cases
where phase change occurs at xed temperature.
Two approaches, nite-difference and nite-element techniques,
are rstly used to solve the phase change problems numerically.
The rst nite-difference studies were carried out by Lam and
Clapeyron in 1831 and Stefan in 1891, concerning ice formation.
In 1912 the results of Neuman were published, establishing the precise
solutions to more general phase change problems. In 1943, London
and Seban [37] analyzed the process of ice formation for different
geometries (cylinder, sphere, and at plate). However, Shamsundar
Fig. 3. Relationships between enthalpy and temperature for isothermal phase
change.

and Srinivasan [38] asserted that London's one-dimensional formula-
tion led to errors by increasing the progress of the solidication
process and proposed a two-dimensional formulation for cylinders. A
one-dimensional formulation was presented by Bonacina et al. to
assess the accuracy of the effective heat capacity method [39]. Rabin
and Korin presented a multidimensional numerical solution based on
the effective heat capacity method for solving transient phase change
problems by using the nite-difference method [40]. The proposed
mathematical method showed a good agreement with two exact
solutions and two numerical solutions based on nite-difference and
nite-element methods.
Although the nite-difference method has been traditionally
employed for solving phase change problems in the past, oscilla-
tions in the front position and/or temperature occur due to a poor
heat balance when the phase change front lies between nodes in a
nite difference solution [41,42]. Nowadays the nite-element
method has been more attractive, because of the ability of the
nite-element method to handle complex coupled thermomechanical
mechanisms with various and complex boundary conditions [43].
Consequently, nite-element methods can either accurately repre-
sent the temperature history or track the phase front, or both. A
few earlier studies based on nite-element formulations were
carried out by Rubinsky and Cravahlo in 1981 [44] and Ettouney
and Brown in 1982 [45], concerning on the location of the solid
liquid interface. To deal with oscillations, Comini et al. [46] and
McAdie et al. [47] developed nite-element based enthalpy
formulations, and applied the three level PredictorCorrector
and Regula-Falsi methods to account for the sudden jumps or
sharp changes in the enthalpytemperature relationships. Bhatta-
charya et al. developed a new enthalpy based method to study
phase change in a multicomponent material using the Galerkin
nite-element method. This method can efciently capture multi-
ple thawing fronts which may originate at any spatial location
within the sample [48]. Hence, nite-element methods are pre-
ferred to be used in enthalpy approaches and the effective heat
capacity method.
Beside the nite-difference and nite-element methods, the
nite-volume and control-volume- nite-element methods were
also employed to handle the phase change problems [39].
3.1.2. Adaptive grid
The problem associated with the xed grid method can be deal
with by utilizing the adaptive grid method to improve the
computational efciency, using the so-called adaptive grid (or
front tracking) numerical schemes to capture the moving bound-
ary [49]. In the adaptive grid method, the exact location of the
moving boundary is evaluated on a grid at each step. There are two
main approaches used to achieve this: the interface-tting grids
and the variable-space grids. The interface-tting grids (also
referred to as the variable time step methods), where a uniform
spatial grid but a non-uniform time step are used, has been
repeatedly employed to solve two-phase and one-dimensional
problems. The another approach is variable-space grids, also
known as the dynamic grids, where the number of spatial intervals
is kept constant and the spatial intervals are adjusted in such a
manner so that the moving boundary lies on a particular grid
point. Thus, in these methods the spatial intervals are a function of
time. The applications of the variable grid method to study the
phase change problem can be found in [5054].
Lacroix and Voller [55] performed a study on simulation of
diffusion/convection controlled solidication processes in a rec-
tangular cavity by using the nite-difference techniques. They
concluded that the xed grid must be ner for solving the melting
problem of a material with a unique melting temperature. Com-
parison between the xed and adaptive grid has also been done by
Bertrand et al. [56]. They found that the adaptive grid method is
better adapted to the melting problem than the xed grid method.
However the adaptive grid method may fail to simulate the
situations where the transition from the liquid to the solid phase
is not a macroscopic surface.
Table 3 presents some references of the mathematical and
numerical studies undertaken with the only consideration of the
conduction heat transfer.
3.2. Conduction and convection heat transfer mechanisms
During the phase transformation especially melting process,
the temperature and concentration gradients in the liquid phase of
PCM keep varying, which results in the movement of the liquid
PCM, named convection occurring under the action of buoyancy
forces due to the density gradients. In the previous numerical
investigations, the convection heat ow in the liquid phase
received less attention than conduction owing to the limited
computational capabilities of computer and the mathematical
complexities to formulate the convection heat transfer during
the phase transformation. The inuence of the natural convection
on the phase change problems were initially considered in 1977.
The study of the possible role of the natural convection in the melt
region has been conducted by Sparrow et al. in 1978 [67]. The
ndings of this study indicated that the natural convection is rst-
order important and has to be accounted in the phase transforma-
tion analysis. A number of researchers [6871] have also reported
that the convection affects not only the rate of melting or
solidication but also the resulting of the distribution of the liquid
phase of the PCM in a multi-component system. One of the
pioneer numerical studies was carried out by Yao and Chen in
1980 [72], by using an approximate solution. They studied the
effect of the natural convection on the phase change process and
concluded that it strongly depends on the Rayleigh number. The
equation of k
e
=k
l
cRa
n
includes an effective of thermal conduc-
tivity was used to describe the inuence of natural convection on
the melting process [73,74].
The inuence of the convection heat transfer on the melting
and freezing processes has been studied by Lamberg et al. [64],
they found that when the effect of natural convection is omitted in
the simulation. The numerical result was poor compared with the
experimental result during the melting process, whilst it showed a
good estimation during the freezing process. Hence it is essential
to model natural convection in the liquid PCM during the melting
process. Besides the numerical methods mentioned in Section 3.1,
the integral method [75], boundary xing method [76], unstruc-
tured nite-element method [77], enthalpy-porosity approach
[7880], coordinate transformation method [81,82] and equivalent
thermal capacity method [83] have been developed for natural
convection simulation.
Yao and Cherney [84] studied the effect of the natural convec-
tion on the melting of a solid PCM around a hot horizontal cylinder
by using the integral method. The results demonstrated that the
integral solution had surprising accuracy when it was compared
with the quasi-steady solution when Stefan number was small.
Rieger et al. [85] investigated numerically the melting process of a
PCM inside a heated horizontal circular cylinder. Both heat con-
duction and convection have been taken into account to treat this
moving boundary problem, and the complex structure of the time-
wise changing physical domain (melt region) have been success-
fully overcome by applying body tted coordinate technique.
Ismail and Silva [86] developed a 2D mathematical model to study
the melting of a PCM around a horizontal circular cylinder
considering the presence of the natural convection in the melt
region. A coordinate transformation technique was used to x the
moving front. The numerical predictions were compared with
available results to establish the validity of the model, and a

Table 3
Numerical studies by only considering the conduction heat transfer.
Reference(s) Dimensional Numerical
simulation method
Assumptions Comment Validation
Silva et al.
[57]
1D Enthalpy based
method
1-D conduction controlled the phase
change process.
The uid ow was fully developed with
temperature independent properties.
Thermal radiation between the wall
and uid was neglected.
The model can be used to predict the dynamic
performance of this kind of vertical rectangular
heat exchange unit with reasonable accuracy
Yes
Krkl et al.
[15]
2D Energy based
equation
The phases were homogeneous.
Conduction heat transfer was
only considered.
Heat loss or gain from the store was
neglected.
Calculating the temperature of air and PCM
at the formerly dened volume nodes.
Yes
Costa et al.
[3]
2D Energy based
equation.
The thermophysical properties of the PCM
and n were independent of temperature.
The PCM was initially in the solid phase.
The PCM was homogenous and isotropic.
Calculations have been made for the melt fraction
and energy is stored by conduction only
Numerical
[33,58]
Fully implicit
nite-difference
method.
Gong et al.
[59]
2D Enthalpy based
method.
The heat-transfer uid was incompressible
and viscous dissipation was negligible.
The uid ow was radially uniform and
the axial velocity was an independent
parameter.
Thermal losses through the outer wall
of the PCM were negligible.
Conduction heat transfer was the
mainly transfer method.
For mode1, simulation has been done by
introducing
hot and cold uids from the same end of the tube.
While for mode 2, simulation has been done by
introducing hot and cold uids from the different
ends of the tube.
No
Standard Galerkin
nite-element
Sharma et al.
[60]
2D Enthalpy based
method
The thermo-physical properties of PCM's and
n material were independent of
temperature but different for solid and liquid
phases.
The PCM was initially in solid phase. The
PCM was homogeneous and isotropic.
The mode of heat transfer was
conduction only.
Calculations were made to study the effect of
thermal conductivity and thermal capacity of fatty
acids, and conductivity of heat exchanger materials
and their effect on melt fraction.
No
Fully implicit
nite-difference
method
Dwarka and
Kim [61]
3D Implicit enthalpy
based method.
Air temperature was constant.
Initial temperatures were the same through
the board.
There was no energy loss to surroundings.
Conduction heat transfer was only
considered.
All thermo-physical properties were
constant except for the heat capacity.
Comparison of the thermal performance of
randomly mixed and laminated PCM drywall
system
No
Dolado et al.
[62]
1D Finite-difference
method
Only conduction heat transfer was
considered inside the PCM plate.
The temperature of the airow was analyzed
in a one dimensional way.
A model was developed to simulate the
performance
of a LHTES, and analyze the heat transfer between
the air and a commercial macroencapsulated PCM.
Yes
Fukai et al.
[63]
2D Enthalpy based
method
The surfaces of the heat exchanger were
adiabatic.
The cross-sections of the tubes was square
The convective heat transfer term was
neglected.
A numerical model predicted well the experimental
outlet uid temperatures and the local
temperatures
in the composite.
Yes
Control-volume
method
Lamberg et al.
[64]
2D Enthalpy based
method and Effective
heat capacity
method
The heat conductivity and density of the PCM
and container were constant.
The physical properties chosen for the PCM
were the average values of solid and
liquid PCM.
The convection heat transfer was negligible
during solidication process.
Both numerical methods gave good estimations for
the temperature distribution. Particularly, when
the temperature range was 2 1C, the effective heat
capacity method was the most precise numerical
method.
Yes
Zivkovic and
Fujii [65]
1D Enthalpy based
method
The effects of natural convection within the
melt were negligible and can be ignored.
The PCM behaved ideally.
The PCM was assumed to have a denite
melting temperature.
The cylindrical container required nearly twice of
the melting time as for the rectangular container of
the same volume and heat transfer area
Yes

satisfactory agreement was found. They concluded that the
numerical model was adequate to represent the physical situation
of the proposed system.
Furzerland [87] investigated the enthalpy method and the
coordinate transformation method through the comparison of
the solutions of specic problems of one dimensional heat transfer
by considering pure convection. One of his conclusions was that
the enthalpy method is very attractive owing to these: it is easy to
program and there are no computational overheads associated
with tracking the moving interface within a specic range of
fusion temperatures.
A summary for various numerical solutions on phase change
problem considering both conduction and convection is presented
in Table 4.
4. Numerical models for various PCMs capsules and packages
It can be easily found that solidliquid PCMs have been widely
used for studies on phase change problems compared to solid
solid, liquidgas and solidgas PCMs. This is due to their large heat
storage capacity and little volume change during the phase change
process. Since solidliquid PCMs experience a phase transition
during the charging process and/or charging process, encapsula-
tion is required for holding the liquid and/or solid phases of the
PCM and keeping it isolate from the surrounding. This ensures the
exibility in frequent phase change processes, and an increase in
heat transfer rate. It is worth noting that the heat transfer
characteristics in the encapsulated PCM would change signi-
cantly depending on the parameters of various encapsulations.
Hence, vast of numerical studies on phase change problems have
extended to engineering applications. In this section, a compre-
hensive review on numerical models for encapsulated PCMs in
terms of macro-encapsulated PCM and microencapsulated PCM is
presented.
4.1. Macro-encapsulated PCMs numerical models
Macro-encapsulation is a common form of encapsulation used
for thermal storage applications. The shape of macro-encapsulate
varies from rectangular panels, spheres to cylinders.
4.1.1. Rectangular encapsulation
The most intensely investigated LHS containment is the rec-
tangular encapsulation because to its simple boundary conditions
and easy work principle. The rst numerical study on rectangular
geometry has been conducted by Shamsundar and Sparrow
[107,108] in 1975 and 1976 by applying the nite-difference
method to investigate the solidication of a at plate. Shamsundar
also used the same method to the case of a square geometry [109]
in 1978.
Hamdan and Elwerr [110] presented a simple 2-D numerical
model to study the melting process of a solid PCM within a
rectangular enclosure. In this model, the rate of melting depends
essentially on the properties of the PCM, such as the thermal
diffusivity, viscosity, conductivity, latent heat of fusion and specic
heat. Natural convection was treated as the dominant heat transfer
mode within the melt region. Conduction was only assumed to
take place within the layer very close to the solid boundary. The
numerical results showed a high agreement with previously
published experimental results [111,112]. Therefore, it can be
effectively used to predict the energy storage performance of the
PCMs contained within the rectangular encapsulation.
Further research was carried out by Lacroix [113] in 2001, who
presented a mathematical model based on the energy conserva-
tion equation to study the contact melting of a sub-cooled PCM
inside a rectangular cavity including the natural convection effect.
Numerical results demonstrated that the melted fraction at the
bottom of the cavity is larger, by an order of magnitude than that
of the conduction dominated melting at the top. The melting
process is essentially governed by the magnitude of the Stefan
number and strongly inuenced by the lateral dimensions of the
cavity. Jiji and Gaye [114] also analytically examined one-
dimensional solidication and melting of a slab with uniform
volumetric energy generation. They found that the low thermal
conductivity of the PCMs present a signicant challenge in the
design of PCM-based electronic cooling systems.
In order to address the inherent drawback of the low thermal
conductivity of PCM, various numerical solutions have been
applied for PCM with different thermal conductivity enhancers
(TCEs). Lacroix and Benmadda [115] numerically analyzed the
melting of PCMs in a rectangular enclosure with horizontal ns.
A 2-D enthalpy model was used to solve the phase change
problem using a xed-gird method. The effects of the number
and length of ns on the melting rate were considered in this
study. It was concluded that a few longer ns were more effective
in the melting rate than a large number of shorter ns. The effect
of vertical ns on melting rate was studied by the same team [116]
by using a xed-grid enthalpy model in 1998. It was found that the
onset of natural convection in the melted zone was delayed when
the distance between the ns was decreased. According to their
results, the optimum n spacing decreases as the Rayleigh number
increases. Akhilesh et al. [117] numerically studied a rectangular
module with vertical ns, which was heated from above. The
numerical model was solved by using nite-difference method.
Only conduction heat transfer was included in this study. The
analysis showed that more ns increase the rate of energy stored
in the melting PCM until a critical value. There is no obvious
performance enhancement by only increasing the number of ns
beyond the value. Similar results were obtained in another
research of Gharebagi and Sezai [118] in 2005 considering an
enclosure with vertical ns added to a horizontal heated wall. The
results indicated that the heat transfer rate to the melting PCM can
be improved by adding ns. In addition, vertically heated walls
with horizontal ns exhibit better performance than horizontally
heated walls with vertical ns.
Table 3 (continued )
Reference(s) Dimensional Numerical
simulation method
Assumptions Comment Validation
Saman et al.
[66]
2D Enthalpy based
equation
The thermo physical properties of PCM
can be assumed constant.
The PCM was initially solid and its
temperature was assumed at a certain value
below the melting point.
The natural convection was taken
into account during charging period.
The numerical model was able to predict quite
accurately the melting time and the heat transfer
rate during the melting with the geometry being
considered
Yes

Table 4
Numerical studies by considering both the conduction and convection.
Reference(s) Dimensional Numerical solution method Assumptions Comment Validation
Costa et al.
[3]
1D Energy based equation The thermophysical properties of the
PCM and n material were independent
of temperature.
The PCM was initially in the solid phase.
The PCM was homogenous and isotropic.
Calculations have been made for the melt
fraction and energy stored under conduction
plus convection condition.
Numerical
[30,58] Fully implicit nite-
difference method
Halawa et al.
[88]
2D Enthalpy based method For the PCM melting process, the PCM
was solid and its temperature was
assumed at a certain value below the
melting point.
For freezing process, the PCM was
initially liquid and its temperature was
assumed at certain value above the
melting point.
Typical characteristics of the melting/freezing
of PCM slabs were discussed. Considerations
in the design of the LHTES were also given.
Numerical [89]
Wang et al.
[90]
2D Energy based method The PCM was pure and homogeneous.
The melting front was an isothermal
interface in reality.
The thermophysical properties of the
PCM were constant, with the exception
of the linear densitytemperature
relation.
The natural convection ow of the liquid
phase was incompressible and laminar
with no viscous dissipation.
The density of PCM was constant in
this study.
Calculations were performed to study the
effect of thermal conductivity and thermal
capacity of fatty acids and conductivity of
heat exchanger materials and their effect
on melt fraction.
No
Finite-volume method
Jana [91] 2D Enthalpy-porosity based
method
The density was considered to be
constant in the unsteady and convective
terms and was allowed to vary only in
the body-force term,
The phase change occurred at a single
temperature,
The temperatures of liquid and solid at
the interface were equal.
Prediction of solidication and melting
processes of pure materials
using moving grids.
Numerical and
Experimental
[92] Finite-volume method
Hamdan and
Al-Hinti
[93]
3D The solidliquid interface was
smooth, plane.
The solid phase was initially at the
melting temperature.
The effect of inertia inside the thermal
boundary layers was negligible relative
to that of buoyancy and viscosity.
The heat loss from the walls of the
enclosure was negligible.
The PCM material was pure, homogenous
and isotropic.
Boussinesq approximation was used in
this analysis.
The melting process of a solid PCM contained
in a rectangular enclosure heated from
a vertical side at a constant heat ux
was investigated analytically.
Experimental
[94]
Zhao et al.
[95]
1D Energy based equation The heat transfer process inside the
molten layer was one-dimensional, and
the tangential temperature gradient was
neglected.
The thermo-physical properties of PCM
were constant.
The process was quasi-steady, and the
acting forces on PCM were balanced at
every moment of the melting.
The expression of melting velocity was
obtained by using the lubrication theory.
Numerical and
Experimental
[96]
Liu et al. [97] 2D Enthalpy porosity method The two-phase mixed region can be
described with the porosity.
The effective thermal conductivity of the
mixture can be calculated as the volume
average of the conductivities of porous
matrix material and PCM
The porous media and uid were not
assumed to be in thermal equilibrium.
The effects of the structural parameters of the
porous media and inlet conditions of HTF on
the thermal performances of LHTES were
numerically analyzed.
Experimental
[98] 3D
Shmueli et al.
[99]
2D Enthalpy porosity method Both solid and liquid phases were
homogeneous and isotropic.
The melting process was axisymmetric.
The ndings of the present study made it
possible to dene the heat transfer
mechanisms from the tube wall to the liquid
Experimental
[100]

Apart from applying ns to PCMs, the effect of metal foams on
heat transfer enhancement in rectangular encapsulated PCMs
was investigated by Tian and Zhao [98] in 2011. The numerical
investigation was based on the two-equation non-equilibrium
heat transfer model, in which the coupled heat conduction and
natural convection were considered at the phase transition and
liquid zones. The numerical results were validated by experimen-
tal data. The main ndings of the investigation were that the heat
conduction rate increases signicantly by using metal foams, due
to their high thermal conductivities, and that natural convection is
suppressed owing to the large ow resistance in metal foams.
In spite of this suppression caused by the metal foams, the overall
heat transfer performance is improved when metal foams are
embedded into PCM. This implies that the enhancement of heat
conduction offsets or exceeds the natural convection loss. The
results also indicated that for different metal foam samples, heat
transfer rate can be further increased by using metal foams with
smaller porosities and bigger pore densities.
4.1.2. Cylindrical containment models
Three modes of cylindrical PCM container congurations are
classied. The rst model is the PCM lling the shell and the heat
transfer uid (HTF) owing through a single tube (pipe model).
In the second mode the PCMlls the tube and the HTF ows parallel
to the tube (tube-tube model). The third cylinder model is the PCM
encapsulated in a shell and HTF ows through the tube outside the
shell to improve the heat transfer of PCMs (tube-shell model).
Ismail and Alves [119] numerically analyzed a pipe model with
the PCM encapsulated in a shell while water as the HTF owing
through the tube. Radial conduction was assumed to be dominated
during the process of solidication and the energy equation for the
PCM was solved in conjunction with the equation governing
the uid bulk temperature as a function of time by employing
the control-volume method. The effect of Biot number, ratio of
diameters and inlet uid temperature on the thermal performance
of the storage unit were presented. Lacroix [105] developed an
enthalpy based model to predict the transient thermal perfor-
mance characteristics of a pipe LHTES with the PCM on the shell
side and the HTF circulating inside the tube. Including the axial
conduction effect in the PCM and employing an enthalpy based
method, the coupled conservation equations governing the heat
transfer in the PCM and HTF were solved iteratively using a nite-
volume based nite-difference technique. Numerical results were
presented for various thermal and geometric parameters such as
shell radius, mass ow rate and inlet temperature of the HTF.
It was concluded that for a given type of PCM, these parameters
must be selected carefully in order to optimize the performance of
the storage unit. Then the tube-tube model was theoretically
studied by Esen and Durmus [120] in 1998. A theoretical model
based on an enthalpy method was developed to predict the effects
of various thermal and geometric parameters, congurations of
the cylindrical container on the whole melting time for different
PCMs. The theoretical results showed that the whole melting time
of PCM depends on not only thermal and geometric parameters of
the container, but also on the thermophysical properties of the
PCM. It was found that a thicker PCM pack would result in longer
melting time.
During 2006, Regin et al. [121] experimentally and numerically
studied the third model that the melting behavior of parafn wax
encapsulated in a cylindrical capsule. The numerical analysis
carried out by using enthalpy method and the results were veried
with the experimental data. The experiments have been done by
the visualization technique without disturbing the actual process
Table 4 (continued )
Reference(s) Dimensional Numerical solution method Assumptions Comment Validation
The molten PCM and the air were
incompressible Newtonian uids,
and laminar ow was assumed in both.
PCM and then to the solid phase at various
locations and instants.
Ye et al.
[101]
2D Enthalpy porosity
approach
Density and dynamic viscosity of the
PCM depended on temperature.
Variation of property values in liquid and
solid phase PCM was imposed using
piecewise linear interpolation.
Constant thermophysical properties
were specied for aluminum.
The numerical study provided a detailed
knowledge regarding interface heat transfer
rate provides a deeper understanding the
heat
transfer mechanisms.
Experimental
[102] and
Numerical
[103]
Tao and He
[104]
2D Enthalpy method The axial heat conduction and viscous
dissipation in the HTF was negligible.
The ow of HTF was treated as one
dimensional uid ow.
The PCM region was treated as an
axisymmetric model.
The effect of natural convection of PCM
during melting was taken into account
with an effective thermal conductivity of
the liquid phase of PCM.
The effects of the non-steady inlet conditions
of HTF on melting time, melting fraction, TES
capacity, solidliquid interface, heat ux on
tube surface and HTF outlet temperature
were analyzed.
Experimental
[105]
Adine and
Qarnia
[106]
2D Conservation energy
equations
The ow of HTF was dynamically
developed.
The axial conduction and viscous
dissipation in the uid were negligible.
The thermophysical properties of the
HTF and PCMs were independent of the
temperature.
The thermophysical properties of solid
and liquid phases of PCM were equal,
except their thermal conductivities.
This parametric study could provide
guidelines
for system thermal performance and design
optimization.
Experimental
[105]
Control volume approach

of melting. Numerical results indicated that the melting process is
chiey governed by the magnitude of the Stefan number, phase
change temperature range and the capsule radius. The analysis
showed that the agreement between the analytical and experimental
results is signicantly improved when the results are obtained by
considering the phase change temperature range and the natural
convection in the liquid phase instead of only considering the
conduction heat transfer.
4.1.3. Spherical containment models
Generally, the PCM in the sphere is subjected to constrained
(the solid and liquid phases have the same density) and uncon-
strained melting (the solid portion drops out of the melting layer
due to its higher density).
The constrained melting of the PCM within a spherical container
was investigated by Khodadadi and Zhang in 2001 [4]. The
computations were based on an iterative, nite-volume numerical
procedure by using the primitive-dependent variables, whereby the
time-dependent continuity, momentum and energy equations in
the spherical coordinate system were solved. With the buoyancy-
driven convection gradually dominated the heat transfer mechan-
ism due to the growth of the melt zone, the solid PCM melted faster
in comparison to the bottom zone. When the buoyancy effects
became remarkable, as many as three time-dependent re-circulat-
ing vortices were observed from the numerical simulation. The
computational ndings were veried through qualitative con-
strained melting experiments using a high-Prandtl number wax.
An experimental and computational investigation directed at
understanding the role of buoyancy-driven convection during con-
strained melting of the PCMinside a spherical capsule was reported by
Tan et al. in 2009 again [122]. The computations were based on an
iterative, nite-volume numerical procedure that incorporated a
single-domain enthalpy formulation for simulation of the phase
change phenomenon. A Darcy's Law-type porous media treatment
linked the effect of the phase change on convection. The computa-
tional predictions pointed to the strong thermal stratication in the
upper half of the sphere that resulted from rising of the molten PCM
along the inner surface of the sphere thus displaced the colder uid.
The measured temperature data and computational results near the
bottom indicate the establishment of an unstable uid layer that
promotes chaotic uctuations and is responsible for waviness of the
bottom of the PCM. Regin et al. [123] also discussed an experimental
and numerical study of a spherical capsule lled with a parafn wax
(melting temperature range of 52.961.6 1C) as the constrained PCM
and placed in a convective environment during the melting process.
Time-resolved temperature readings were obtained by using a rake of
nine thermocouples that were placed on a horizontal plane passing
through the center of the sphere. Keeping the initial temperature of
the PCM constant, the temperature of the convective heat transfer
uid varied from 70 to 82 1C corresponding to Stefan numbers of
0.11430.2501. A 1D model for phase change by employing the
enthalpy method and nite-difference formulation was developed.
The model was then used to investigate the effect of capsule size and
the Stefan number on the instantaneous heat ux, liquid fraction and
melting time.
Considering the spherical geometries, apart from several ana-
lytical investigations of diffusion controlled phase change pro-
cesses reported before 1980s, the rst study on the unconstrained
melting of the PCM within spherical container was carried out by
Moore and Bayazitoglu [124] in 1982. A mathematical model was
developed by assuming that the liquid remains at the fusion
temperature during the whole period and solved by employing
the perturbation method. The solidliquid interface positions and
the temperature proles were determined for various Stefan and
Fourier numbers. The experimental data agreed well with the
simulation results. During 1987, Roy and Sengupta [125] per-
formed a similar study and found an analytical solution for the
melting rate at the lower surface of the solid PCM in contact with
the heated surface. Later, Roy and Sengupta [126] performed
another study on the unconstrained melting process inside a
sphere. Two zones of the melting were identied: a thin melting
layer at the bottom and a thicker layer at the top of the sphere.
They found that a signicant amount of melting took place at the
upper surface due to the signicance of natural convection,
whereas in previous research [124,125], the effect of natural
convection on the melting process was neglected.
4.2. Microencapsulated PCM
Micro-capsulation is a kind of tiny capsules (around 0.1
1000 m in diameter) [127]. Micro-encapsulation has a higher
heat transfer rates as compared to that of macro-encapsulation
[128,129] due to a substantially higher surface area to volume
ratio. Meanwhile, microencapsulation can tolerate the change in
volume during phase change process [130]. Microencapsulation
hence is a prominent technology for use of PCM in building
materials. With the advent of PCM implemented in gypsum board,
plaster, concrete and/or other wall covering materials, thermal
storage can be part of the building structure even for light weight
buildings.
4.2.1. Incorporation with building components
PCMs can be incorporated with gypsum board, plaster, concrete
and/or other wall covering materials for the temperature control
application purpose.
In general, storage density and phase change temperatures are
very important parameters, since they determine the storage
system size and human thermal comfort. Further, accurate estima-
tion of the PCM enthalpy variation in the working temperature
range is essential for correct mathematical modeling of the storage
system [131].
The rst research on the application of PCMs in buildings has
been carried out by Telkes in 1975 [132]. Na
2
S0
4
U10 H
2
0 with a
melting temperature of 32 1C was incorporated with walls, parti-
tions, ceilings and oors to serve as temperature regulators.
Determination criteria for the optimal phase change temperature
and storage density were provided by Zhou [133] as following:
T
m;opt
T
r
Q=ht
stor
26
T
r
t
d
T
d
t
n
T
n
=t
d
t
n
27
D
opt
t
n
hT
m;opt
T
n
=L 28
where T
m;opt
is optimal phase change temperature and D
opt
is
optimal thickness of the wall; Q represents the heat absorbed by
unit area of the room surface; T
r
is the average room temperature;
t
stor
is the heat storage time; the subscripts n and d represent night
and daytime respectively.
Kuznik et al. [134] developed a one-dimensional model to
optimize the thickness of phase change material only considering
pure conduction. The simulation was based on a lightweight wall
and interior/exterior temperature evolutions within a period of
24 h. The results showed that an optimal PCM thickness exists.
With the optimal PCM thickness of 1 cm, the thermal inertia of the
building was doubled and the PCM can reduce the temperature
uctuations inside rooms. To provide guidelines useful in selecting
an optimal PCM for building wallboard, a numerical model was
developed to determine the optimal melting temperature of PCM
by Neeper [135]. The results indicated that the optimal melting
temperature depends on the average room temperature. When the

melting temperature of PCM closed to the average room tempera-
ture, the maximum diurnal energy storage was achieved. Xiao
et al. [136] established a simplied theoretical model to optimize
an interior PCM for energy storage in a lightweight room. In this
paper energy conservation equations were presented to calculate
the optimal phase change temperature and the total amount of
latent heat capacity. The analytical results agreed with Neeper's
nding that the optimal melting temperature depends on the
average indoor air temperature. Further, and it also depends on
the radiation absorbed by the PCM panels.
A 1-D enthalpy model for analyzing the thermal performance
of a shape-stabilized PCM oor was developed by Xu et al. [137].
By using the model, the inuence of various parameters (melting
temperature of PCM, thickness of PCM layer, latent heat of fusion,
and thermal conductivity of PCM) on the thermal performance
was analyzed by using a fully implicit the nite-difference method.
It was found that for a given position or weather condition, the
suitable melting temperature of PCM is roughly equal to the
average indoor air temperature of sunny winter days, the heat of
fusion and the thermal conductivity of PCM should be larger than
120 kJ/kg and 0.5 W/(m K), respectively, and the thickness of shape
stabilized PCM plate used under the oor should not be larger
than 20 mm. Zhou et al. [138] extended the investigation of the
inuence of various parameters from oor to walls and ceiling by
using a veried enthalpy model. The model was validated by
experimental work, and simulated results agreed well with the
measured data. The numerical results showed that for the present
conditions, the optimal melting temperature is about 20 1C and
the heat of fusion should not be less than 90 kJ/kg, thin PCM plates
with large areas are advantageous and the effect of PCM plates
located at the inner surface of interior wall was superior to that of
exterior wall (the south wall). A similar result to Xu et al. [137] was
also obtained for the condition of the PCM wall and ceiling that the
thermal conductivity should be larger than 0.5 W/(mK).
Athienitis et al. [139] experimentally and numerically studied
the application of PCM in building envelope components for
thermal storage. A 1-D non-linear enthalpy model was developed
to simulate the transient heat transfer process in the walls using
explicit nite difference scheme. The numerical result showed that
the utilization of the PCM gypsum board may reduce the max-
imum room temperature by about 4 1C during the day time and
can reduce the heating load at night signicantly.
The numerical simulation of passive heating with PCMs was
carried out by Chen et al. [140] by using a 1-D nonlinear
mathematical model. The effective heat capacity method was used
to simulate the PCM heat transfer problem by using the software
MATLAB. The result indicated that applying proper PCM to the
inner surface of the north wall in the ordinary room can not only
enhance the indoor thermal comfort dramatically, but also
increase the utilization rate of the solar radiation. Consequently
the heating energy consuming is decreased and the goal of saving
energy has been achieved. The energy-saving rate of heating
season can get to 17% or higher.
4.2.2. Microencapsulated phase change slurry (MPCS)
Microencapsulated phase change slurry (MPCS) as a new
technique has been getting more and more attention since this
technique serves as both the heat transfer uids and energy
storage media. Consequently it can increase the thermal storage
capacity and improve the energy efciency of the thermal system.
Due to the MPCS' increasingly important role, some review papers
on the properties, heat transfer characteristics and application of
MPCS have been published in recent years [127,141143]. As a
kind of suspension, the heat transfer and uid ow characteristics
of MPCS are different from those of the marco-encapsulated PCMs.
In accordance to the review of Delgado et al. [142], the heat
transfer characteristics of the MPCS can be classied into forced
convection and natural convection.
4.2.2.1. Laminar forced convection heat transfer. The MPCS used to
enhance convective heat transfer is generally laminar due to the high
viscosity. Charunyakorn et al. [144] rstly conducted a numerical
simulation of MPCS ow in circular tubes under different boundary
conditions. The energy equation was formulated by taking into
consideration both the heat absorption (and release) during phase
transition and solved using an implicit nite difference method. Their
parametric study showed that the bulk Stefan number and volumetric
particle concentration are the two dominant parameters. Their model
also showed that the heat transfer enhancement ratio of the uid
could be 24 times of that of water. Zhang and Faghri [145] modied
the model of Charunyakorn [144] to include the effects of the
microcapsule walls, initial sub-cooling and the melting temperature
range. A temperature transforming model was used to solve the
melting in the microcapsule by instead of a quasi-steady model.
The numerical results showed a very good agreement with the
experimental results of Goel et al. [146]. Their numerical results
suggested that the differences between the experimental results of
Goel et al. and numerical predictions of Charunyakom could be
decreased greatly if these additional factors are considered. The
authors however did not give any correlation or other similar
criteria that could be used for future design.
Complicated source terms are applied to the theoretical models
above to account for the phase change effects. In order to simplify
the simulation process, Alisetti and Roy [147] developed a simpler
effective specic heat model to study the heat transfer in MPCS. The
results showed that the exact form of the specic heat function is not
critical as long as the latent heat is incorporated correctly within a
nite melting temperature range. Based on the concept of effective
specic heat, Roy and Avanic [148] developed a model to analyze the
laminar forced convection heat transfer to a PCM suspension in a
circular duct under constant wall heat ux condition. The model has
been veried with experimental data. The results demonstrated that
Stefan number is the only parameter that has a signicant impact on
the thermal performance. Sub-cooling effect only can affect the heat
transfer performance at very low heat uxes or when the inlet
temperature is much lower than the melting point. A simple
correlation to predict the wall temperature rise as a function of the
tube length was also developed for future design.
Hu and Zhang [149] introduced a model to provide a novel insight
for the forced convective heat transfer enhancement of MPCS owing
through a circular tube with constant heat ux. The model was
developed based on effective specic heat capacity method and used
to analyze the inuence of various factors on heat transfer enhance-
ment. A modied Nusselt number was introduced since the conven-
tional Nusselt number correlations for internal ow of single phase
uids are not suitable for accurately describing the heat transfer
enhancement with MPCS. The modied local Nusselt number is
inversely proportional to the dimensionless wall temperature. The
numerical results indicated that the exact nature of the phase change
process strongly affects the degree of convective heat transfer
enhancement. The Stefan number and the PCM concentration
resulted to be the most important parameters on the improvement
of heat transfer in MPCS. Zeng et al. [150] analyzed experimentally
and numerically the enhanced convective heat transfer mechanism
of MPCS in the thermal fully developed range. An enthalpy model
was proposed and developed for simulating the forced convective
heat transfer characteristics of laminar ow with MPCS in a circular
tube under constant wall heat ux. The results showed that in
the phase change heat transfer region the Stefan number and the
dimensionless phase change temperature range number are the

most important parameters inuencing the Nusselt number uctua-
tion prole and the dimensionless wall temperature. The bulk uid
Reynolds number, particle diameter and the microcapsule volumetric
concentration also inuence the Nusselt number prole and the
dimensionless wall temperature but are independent of the phase
change process.
4.2.2.2. Turbulent forced convection heat transfer. A majority of past
studies have considered laminar forced convection heat transfer to
MPCM and only a few numerical investigations on turbulent
forced convection heat transfer have been reported. However,
the turbulent ows occur in majority of engineering applications
of MPCS. In order to investigate the turbulent heat transfer to PCM
suspensions, Roy and Avanic [151] presented an effective specic
heat capacity model for PCM suspensions in a circular tube with
constant heat ux under turbulent ow condition. The numerical
results were found to agree with previously published
experimental data and showed that this effective specic heat
model can effectively model turbulent heat transfer with PCM
suspensions. The most primary parameter in heat transfer was the
Stefan number. The melt temperature range and degree of sub-
cooling are other two important parameters. Royon and Guiffant
[152] developed a numerical model to describe the thermal
behavior of a slurry of PCM slurries ow in a circular duct under
different ow parameters (Reynolds number). The simulation
results demonstrated that the inuence of Reynolds number on
the minimum length of the heat exchanger is relatively weak.
4.2.2.3. Natural convection heat transfer. Inaba et al. [153]
developed a 2-D model to study the uid ow and heat-transfer
characteristics for RayleighBnard natural convection of non-
Newtonian PCM slurry. They found that Rayleigh number,
Prandtl number and aspect ratio could be the primary factors for
most of Newtonian and non-Newtonian uids to evaluate the
natural convection. A modied Stefan number was dened in the
paper to evaluate the natural convection in a PCM slurry. Later,
Inaba et al. [154] presented another 2-D model to study the
natural convection heat transfer of a MPCS. It was found that the
natural convection effect and heat transfer enhancement were due
to the contribution of the latent heat transfer.
4.3. Summary of numerical models
It can be seen that the numerical models collected above for
encapsulated PCMs all have their limitations. As the practical
phase transformations in different applications are complicated,
and the thermal conditions are not ideal. Various assumptions
were set up for each numerical solution according to different
numerical simulation purposes. There is none of numerical models
can describe the phase change problem properly. Hence, the
numerical modeling of LHTES only approximately approaches
but not accurately presents the practical thermal performance of
LHTES. Therefore, more investigation should be carried out to
study the phase change problems so as to provide a more detailed
insight of the actual heat transfer behaviors inside PCMs during
the phase transformations.
5. Conclusions
This paper presented a comprehensive review of mathematical
and numerical methods applied to the solutions of phase change
problems. Firstly heat transfer mechanisms during the charging and
discharging processes were discussed in this review. The heat
transfer mechanisms play the most important role on the numerical
results. Consequently, it is important to weight the percentages of
conduction and convection heat transfer taken in each stage. Then,
it presented the fundamental mathematical descriptions of the
phase change phenomenon, the Stephan problem and Neumann
problem. At last, a description of the numerical solutions for phase
change problems based on considering only conduction and con-
sidering also convection. The PCMs have been applied in thermal
energy storage applications widely due to their attractive features.
Various numerical models for encapsulated PCMs in terms of
macro-encapsulated and micro-encapsulated PCM were also ana-
lyzed in this paper. Nevertheless, the incorporation of PCMs in a
particular application calls for an analysis that will enable the
researcher to understand the heat transfer principle inside PCM
correctly. From this survey, some further works are recommended:
Further investigation should be carried out to obtain adequate

knowledge of the phase change problems so as to minimize
the assumptions made in the numerical models.
Development of new numerical model to broaden its applica-

tion in phase change problems and improve the accuracy of its
results.
Assessment of the stability and convergence on the numerical

results is always required, and the numerical predictions should
always be validated by using appropriate experimental data.
Unied dimensionless parameters needs to be developed to

make the gained knowledge applicable to other cases.
Life cycle analysis, both economical and energetic feasibility of

PCM application should be performed as LHTES is a more
expensive form of thermal storage, but few scientic publica-
tions covered this matter.
References
[1] Trelles JP, Duy JJ. Numerical simulation of a porous latent heat thermal
energy storage for thermoelectric cooling. Appl Therm Eng 2003;23:
164764.
[2] Ho CJ, Chu CH. Numerical simulation of heat penetration through a vertical
rectangular phase change material/air composite cell. Int J Heat Mass Transf
1996;39(9):178595.
[3] Costa M, Buddhi D, Oliva A. Numerical simulation of a latent heat thermal
energy storage system with enhanced heat conduction. Energy Convers
Manag 1998;39:31930.
[4] Khodadadi JM, Zhang Y. Effects of buoyancy-driven convection on melting
within spherical containers. Int J Heat Mass Transf 2001;44(8):160518.
[5] Hasan A. Phase change material energy storage system employing palmitic
acid. Sol Energy 1994;52:14354.
[6] Lacroix M, Duong T. Experimental improvements of heat transfer in a latent
heat thermal energy storage unit with embedded heat sources. Energy
Convers Manag 1998;39(8):70316.
[7] Velraj R, Seeniraj RV, Hafner B, Faber C, Schwarzer K. Heat transfer
enhancement in a latent heat storage system. Sol Energy 1999;65(3):17180.
[8] Buddhi D, Bansal NK, Sawney RL, Sodha MS. Solar thermal storage systems
using phase change materials. Int J Energy Res 1988;12:54755.
[9] Zhang Y, Yi J. A simple method, the T-history method of determining the heat
of fusion specic heat and thermal conductivity of phase change materials.
Meas Sci Technol 1999;10:2015.
[10] Sari A, Kaygusuz K. Thermal energy storage system using stearic acid as a
phase change material. Sol Energy 2001;71(6):36576.
[11] Zhang Z, Bejan A. The problem of time-dependant natural convection
melting with conduction in the solid. Int J Heat Mass Transf 1989;23:
244757.
[12] Hasnain SM. Review on sustainable thermal energy storage technologies,
part 1: heat storage materials and techniques. Energy Convers Manag
1998;39(11):112738.
[13] Sasaguchi K, Viskanta R. Phase change heat transfer during melting and
re-solidication of the melt area around cylindrical heat source (s)/sink(s).
J Energy Resour Technol 1989;111:439.
[14] Radhakrishnan KB, Balakrishnan AR. Heat transfer analysis of thermal energy
storage using phase change materials. Heat Recovery Syst CHP 1992;12(5):
42735.
[15] Krkl A, Wheldon A, Hadley P. Mathematical modeling of the thermal
performance of a phase-change material store: cooling cycle. Appl Therm
Eng 1996;16(7):61323.
[16] Cheng TF. Numerical analysis of non linear multiphase Stefan problems.
Comput Str 2000;75:22533.

[17] zisik MN. Heat conduction. New York: Wiley; 1980.
[18] Alexiades V. Mathematical modeling of melting and freezing processes.
Washington: Hemisphere Publishing Corporation; 1993.
[19] Lynch DR. Unied approach to simulation on deforming elements with
application to phase change problems. J Comput Phys 1982;47(3):387410.
[20] Lynch DR, Sullivan JM. Heat conservation in deforming element phase
change simulation. J Comput Phys 1985;57:30317.
[21] Zabaras N, Ruan Y. A deforming nite element method analysis of inverse
Stefan problems. Int J Numer Methods Eng 1989;28:295313.
[22] Zabaras N, Ruan Y. Moving and deforming nite element simulation of
two-dimensional Stefan problems. Commun Appl Numer Methods 1990;6:
495506.
[23] Ruan Y, Zabaras N. An inverse nite-element technique to determine the
change of phase interface location in 2-dimensional melting problems.
Commun Appl Numer Methods 1991;7:32545.
[24] Murray WD, Landis F. Numerical and machine solutions of transient heat-
conduction problems involving melting or freezing. Heat Transf 1959;81:
10612.
[25] Ciment M, Guenther RB. Numerical solution of a free boundary value
problem for parabolic equations. Appl Anal 1974;4:3962.
[26] Lazaridis A. A numerical solution of the multidimensional solidication (or
melting) problem. Int J Heat Mass Transf 1970;13:145977.
[27] Basu B, Date AW. Numerical modeling of melting and solidication pro-
blems: a review. Sadhana 1988;13:169213.
[28] Voller VR, Swaminathan CR, Thomas BG. Fixed grid techniques for phase
change problems: a review. Int J Numer Methods Eng 1990;30:87598.
[29] Kim C, Kaviani M. A numerical method for phase change problems. Int J Heat
Mass Transf 1990;33(12):272134.
[30] Voller VR. Numerical heat transfer. Part B 1990;17:155.
[31] Hunter LW, Kuttler JR. The enthalpy method for heat conduction problems
with moving boundaries. J Heat Transf Trans ASME 1989;111:23942.
[32] Voller VR. Fast implicit nite-difference method for the analysis of phase
change problems. Numer Heat Transf Part B 1990;17:15569.
[33] Sharma A, Tyagi VV, Chen CR, Buddhi D. Review on thermal energy storage
with phase change materials and applications. Renew Sustain Energy Rev
2009;13(2):31845.
[34] Tihonov AN, Samarskii AA. The equations of mathematical physics [2nd ed.].
Moscow: Gittl; 1953.
[35] Albasiny EL. The Solution of Nonlinear Heat Conduction Problems on the
Pilot ACE. Proc. IEE Part B: Radio Electron Eng 1956;103(1):15862.
[36] Peippo K, Kauranen P, Lund P. A multicomponent PCM wall optimized for
passive solar heating. Energy Build 1991;17(4):25970.
[37] London AL, Seban RA. Rate of ice formation. Trans. ASME 1943;65:7718.
[38] Shamsundar N, Srinivasan R. Analysis of energy storage by phase change
with an array of cylindrical tubes. New York: American Society of Mechanical
Engineers; 1978; 3540.
[39] Bonacina C, Comini G, Fasano A, Primicerio M. Numerical solution of phase-
change problems. Int J Heat Mass Transf 1973;16:182532.
[40] Rabin Y, Korin E. An efcient numerical solution for the multidimensional
solidication (or melting) problem using a microcomputer. Int J Heat Mass
Transf 1993;36(3):67383.
[41] Bell GE. On the performance of the enthalpy method. Int J Heat Mass Transf
1982;25:5879.
[42] Voller VR, Cross M. Accurate solutions of moving boundary problems using
the enthalpy method. Int J Heat Mass Transf 1981;24:54556.
[43] Nedjar B. An enthalpy-based nite element method for nonlinear heat
problems involving phase change. Comput Str 2002;80:921.
[44] Rubinsky B, Cravahlo EG. A nite element method for the solution of one-
dimensional phase change problems. Int J Heat Mass Transf 1981;24:19879.
[45] Ettouney HM, Brown RA. Finite-element methods for steady solidication
problems. J Comput Phys 1982;49:11850.
[46] Comini G, Giudice SD, Saro O. A conservative algorithm for multidimensional
conduction phase change. Int J Numer Methods Eng 1990;30:697709.
[47] McAdie RL, Cross JT, Lewis RW, Gethin DT. A nite element enthalpy
technique for solving coupled nonlinear heat conduction/mass diffusion
problems with phase change. Int J Numer Methods Heat Fluid Flow
1995;5:90721.
[48] Bhattacharya M, Basak T, Ayappa KG. A xed-grid nite element based
enthalpy formulation for generalized phase change problems: role of super-
cial mushy region. Int J Heat Mass Transf 2002;45:488198.
[49] arler B, Kuhn G. Dual reciprocity boundary element method for convective-
diffusive solidliquid phase change problems, part 1. Formul Eng Anal Bound
Elem 1998;21:5363.
[50] Heitz WL, Westwater JW. Extension of the numerical method for melting
and freezing problems. Int J Heat Mass Transf 1970;13:13715.
[51] Rathjen KA, Jiji LM. Heat conduction with melting or freezing in a corner.
J Heat Transf 1971;93:1014.
[52] Crank J, Gupta RS. A method for solving moving boundary problems in heat
ow using cubic splines or polynomials. J Inst Math Appl 1972;10:296304.
[53] Gupta RS. Moving grid method without interpolations. Comput Methods
Appl Mech Eng 1974;4:14352.
[54] Zerroukat M, Chatwin CR. Computational moving boundary problems. New
York: Wiley; 1994.
[55] Lacroix M, Voller VR. Finite difference solutions of solidication phase
change problems: transformed versus xed grids. Numer Heat Transf Part
B 1990;17:2541.
[56] Bertrand O, Binet B, Combeau H, Couturier S, Delannoy Y, Gobin D, et al.
Melting driven by natural convection. A comparison exercise: rst results. Int
J Ther Sci 1999;38:526.
[57] Silva PD, Goncalves LC, Pires L. Transient behavior of a latent-heat thermal
energy store: numerical and experimental studies. Appl Energy 2002;73:
8398.
[58] Goodrich LE. Efcient numerical technique for one-dimensional thermal
problems with phase change. Int J Heat Mass Transf 1978;21:61521.
[59] Gong ZX, Mujumdar AS. Finite-element analysis of cyclic heat transfer in a
shell-and-tube latent heat energy storage exchanger. Appl Therm Eng
1997;17(6):58391.
[60] Sharma A, Wona LD, Buddhi D, Parka JU. Technical note: numerical heat
transfer studies of the fatty acids for different heat exchanger materials on
the performance of a latent heat storage system. Renew Energy 2005;30:
217987.
[61] Dwarka K, Kim JS. Thermal analysis of composite phase change drywall
systems. Trans ASME J Sol Energy Eng 2005;127:3526.
[62] Dolado P, Lazaro A, Marin JM, Zalba B. Characterization of melting and
solidication in a real scale PCM-air heat exchanger: Numerical model and
experimental validation. Energy Convers Manag 2011;52:1890907.
[63] Fukai J, Hamada Y, Morozumi Y, Miyatake O. Improvement of thermal
characteristics of latent heat thermal energy storage units using carbon-
ber brushes: experiments and modeling. Int J Heat Mass Transf 2003;46:
451325.
[64] Lamberg P, Lehtiniemi R, Henell AM. Numerical and experimental investiga-
tion of melting and freezing processes in phase change material storage. Int J
Therm Sci 2004;43:27787.
[65] Zivkovic B, Fujii I. An analysis of isothermal phase change of phase change
material within rectangular and cylindrical containers. Sol Energy
2001;70:5161.
[66] Saman W, Bruno F, Halawa E. Thermal performance of PCM thermal storage
unit for a roof integrated solar heating system. Sol Energy 2005;78(2):3419.
[67] Sparrow EM, Patankar SV, Ramadhyani S. Analysis of melting in the presence
of natural convection in the melt region. ASME J Heat Transf 1977;99:5206.
[68] Ramsey JW, Sparrow EM. Melting and natural convection due to a vertical
embedded heater. J Heat Transf 1978;100:36870.
[69] Hale NW, Viskanta R. Photographic observation of the solidliquid interface
motion during melting of a solid heat from an isothermal vertical wall. Lett
Heat Mass Transf 1978;5:32937.
[70] Bathelt AG, Viskanta R, Leidenfrost W. An experimental investigation of
natural convection in the melted region around a heated horizontal cylinder.
J Fluid Mech 1979;90:22739.
[71] Harrison C, Weinberg F. The inuence of convection on heat transfer in liquid
tin. Metall Trans B 1985;16:3557.
[72] Yao LS, Chen FF. Effects of natural convection in the melted region around a
heated horizontal cylinder. ASME J Heat Transf 1980;102:66772.
[73] Farid MM, Husian RM. An electrical storage heater using the phase change
method of heat storage. Energy Convers Manag 1990;30(3):21930.
[74] Farid MM, Hamad FA, Abu-Arabi M. Melting and solidication in multi
dimensional geometry and presence of more than one interface. Energy
Convers Manag 1998;39(8):80918.
[75] zisik MN. Phase-change problems. Heat conduction. Wiley-Interscience;
1993.
[76] Fomin SA, Saitoh TS. Melting of unxed material in spherical capsule with
non-isothermal wall. Int J Heat Mass Transf 1999;42:4197205.
[77] Omari KEl, Kousksou T, Guer YL. Impact of shape of container on natural
convection and melting inside enclosures used for passive cooling of
electronic devices. Appl Therm Eng 2011;31:302235.
[78] Voller VR, Cross M, Markatos NC. An enthalpy method for convection/
diffusion phase change. Int J Numer Methods Eng 1987;24:27184.
[79] Brent AD, Voller VR, Reid KJ. Enthalpy-porosity technique for modeling
convectiondiffusion phase change: application to the melting of a pure
metal. Numer Heat Transf 1988;13:297318.
[80] Shatikian V, Ziskind G, Letan R. Numerical investigation of a PCM-based heat
sink with internal ns. Int J Heat Mass Transf 2005;48:3689706.
[81] Prusa J, Yao LS. Melting around a horizontal heated cylinder: part I
perturbation and numerical solutions for constant heat ux boundary
condition. ASME J Heat Transf 1984;106:37684.
[82] Prusa J, Yao LS. Melting around a horizontal heated cylinder: part II
numerical solution for isothermal boundary condition. ASME J Heat Transf
1984;106:46972.
[83] Morgan K. A numerical analysis of freezing and melting from a vertical wall.
Comput Methods Appl Mech Eng 1981;28:27584.
[84] Yao LS, Cherney W. Transient phase-change around a horizontal cylinder. Int
J Heat Mass Transf 1981;24:197181.
[85] Rieger H, Projahn U, Beer H. Analysis of the heat transport mechanisms
during melting around a horizontal circular cylinder. Int J Heat Mass Transf
1982;25:13747.
[86] Ismail KAR, da Silva MGE. Melting of PCM around a horizontal cylinder with
constant surface temperature. Int J Therm Sci 2003;42:114552.
[87] Furzerland RM. A comparative study of numerical methods for moving
boundary problems. J Inst Math Appl 1980;26:41129.
[88] Halawa E, Bruno F, Saman W. Numerical analysis of a PCM thermal storage
system with varying wall temperature. Energy Convers Manag 2005;46:
2592604.

[89] Vakilaltojjar SM. Phase change thermal storage system for space heating and
cooling [Ph.D. thesis]. School AME: University of South Australia; 2000.
[90] Wang S, Faghri A, Bergman TL. A comprehensive numerical model for
melting with natural convection. Int J Heat Mass Transf 2010;53:19862000.
[91] Jana S. A numerical method to compute solidication and melting processes.
Appl Math Model 2007;31:93119.
[92] Wolff F, Viskanta R. Solidication of a pure metal at a vertical wall in the
presence of liquid superheat. Int J Heat Mass Transf 1988;31:173544.
[93] Hamdan MA, Al-Hinti I. Analysis of heat transfer during the melting of a
phase-change material. Appl Therm Eng 2004;24:193544.
[94] Bnard C, Gobin D, Martinez F. Melting in rectangular enclosures: experi-
ments and numerical simulations. J Heat Transf 1985;107:794802.
[95] Zhao Y, Chen W, Sun F. Effects of convection and density difference on
contact melting around a cylinder. J Energy Convers Manag 2010;51:2537.
[96] Moallemi MK, Viskanta R. Melting around a migrating heat source. J Heat
Transf 1985;107:4519.
[97] Liu Z, Yao Y, Wu H. Numerical modeling for solidliquid phase change
phenomena in porous media: Shell-and-tube type latent heat thermal
energy storage. Appl Energy 2013;112:122232.
[98] Tian Y, Zhao C. A numerical investigation of heat transfer in phase change
materials (PCMs) embedded in porous metals. Energy 2011;36:553946.
[99] Shmueli H, Ziskind G, Letan R. Melting in a vertical cylindrical tube:
numerical investigation and comparison with experiments. Int J Heat Mass
Transf 2010;53:408291.
[100] Katsman L. Investigation of phase change in cylindrical geometry with internal
ns [M.Sc. thesis]. Heat Transfer Laboratory, Department of Mechanical Engi-
neering. Ben-Gurion University of the Negev: Beer-Sheva, Israel; 2006.
[101] Ye W, Zhu D, Wang N. Fluid ow and heat transfer in a latent thermal energy
unit with different phase change material (PCM) cavity volume fractions.
Appl Therm Eng 2012;42:4957.
[102] Gau C, Viskanta R. Melting and solidication of a pure metal on a vertical
wall. J Heat Transf 1986;108:17481.
[103] Lacroix M. Predictions of natural convection-dominated phase change
problems by the vorticitye-velocity formulation of the Naviere-Stokes
equations. Numer Heat Transf Part B: Fundam: Int J Comput Method
1992;22:7993.
[104] Tao YB, He YL. Numerical study on thermal energy storage performance of
phase change material under non-steady-state inlet boundary. Appl Energy
2011;88:41729.
[105] Lacroix M. Numerical simulation of a shell-and-tube latent heat thermal
energy storage unit. Sol Energy 1993;50(4):35767.
[106] Adine HA, Qarnia HE. Numerical analysis of the thermal behavior of a shell-
and-tube heat storage unit using phase change materials. Appl Math Model
2009;33:213244.
[107] Shamsundar N, Sparrow EM. Analysis of multidimensional conduction phase
change via the enthalpy model. Int J Heat Mass Transf 1975;97:33340.
[108] Shamsundar N, Sparrow EM. Effect of density change on multidimensional
conduction phase change. Int J Heat Mass Transf ASME 1976;98:5507.
[109] Shamsundar N. Comparison of numerical methods for diffusion problems
with moving boundaries. In: Wilson DG, Solomon AD, Boggs PT, editors.
Moving boundary problems. New York: Academic Press; 1978.
[110] Hamdan MA, Elwerr FA. Thermal energy storage using a phase change
material. Sol Energy 1996;56(2):1839.
[111] Webb B, Viskanta R. Analysis of heat transfer during melting of a pure metal
from an isothermal vertical wall. Numer Heat Transf 1986;9:53958.
[112] Yueng W. Engineering analysis of heat transfer during melting in vertical
rectangular enclosures. Int J Heat Mass Transf 1989;32:68996.
[113] Lacroix M. Contact melting of a phase change material inside a heated
parallelepedic capsule. Energy Convers Manag 2001;42:3544.
[114] Jiji LM, Gaye S. Analysis of solidication and melting of PCM with energy
generation. Appl Therm Eng 2006;26:56875.
[115] Lacroix M, Benmadda M. Numerical simulation of natural convection domi-
nated melting and solidication from a nned vertical wall. Numer Heat
Transf Part A Appl 1997;31(1):7186.
[116] Lacroix M, Benmadda M. Analysis of natural convection melting from a
heated wall with vertically oriented ns. Int J Numer Methods Heat Fluid
Flow 1998;8(4):46578.
[117] Akhilesh R, Narasimhan A, Balaji C. Method to improve geometry for heat
transfer enhancement in PCM composite heat sinks. Int J Heat Mass Transf
2005;48(13):275970.
[118] Gharebaghi M, Sezai I. Enhancement of heat transfer in latent heat storage
modules with internal ns. Numer Heat Transf Part A Appl 2008;53
(7):74965.
[119] Ismail KAR, Alves CLF. Analysis of shell-tube PCM storage system. In:
Proceedings of the eighth international heat transfer conference, San Fran-
cisco; 1986. p. 178186.
[120] Esen M, Durmus A. Geometric design of solar-aided latent heat store
depending on various parameters and phase change materials. Sol Energy
1998;62(1):1928.
[121] Regin AF, Solanki SC, Saini JS. Latent heat thermal energy storage using
cylindrical capsule: numerical and experimental investigations. Renew
Energy 2006;31:202541.
[122] Tan FL, Hosseinizadeh SF, Khodadadi JM, Fan L. Experimental and computa-
tional study of constrained melting of phase change materials (PCM) inside a
spherical capsule. Int J Heat Mass Transf 2009;52:346472.
[123] Regin AF, Solanki SC, Saini JS. Experimental and numerical analysis of melting
of PCM inside a spherical capsule. In: Proceedings of the AIAA-2006-3618,
9th AIAA/ASME joint thermophysics and heat transfer conference. California:
San Francisco; June 58, 2006.
[124] Moore FE, Bayazitoglu Y. Melting within a spherical enclosure. J. Heat Transf
ASME 1982;104:1923.
[125] Roy SK, Sengupta S. The melting process within spherical enclosures. J. Heat
Transf ASME 1987;109:4602.
[126] Roy SK, Sengupta S. Gravity-assisted melting in a spherical enclosure: effects
of natural convection. Int J Heat Mass Transf 1990;33:113547.
[127] Zhang P, Ma ZW, Wang RZ. An overview of phase change material slurries:
MPCS and CHS. Renew Sustain Energy Rev 2010;14:598614.
[128] Kasza KE, Chen MM. Improvement of the performance of solar energy or
waste heat utilization systems by using phase-change slurry as an enhance-
ment heat transfer storage uid. J Sol Energy Eng 1985;107:22936.
[129] Hawlader MNA, Uddin MS, Khin MM. Microencapsulated PCM thermal-
energy storage system. Appl Energy 2003;74:195202.
[130] Tyagi VV, Kaushik SC, Tyagi SK, Akiyama T. Development of phase change
materials based microencapsulated technology for buildings: a review.
Renew Sustain Energy Rev 2011;15:137391.
[131] Kousksou T, Arid A, Jamil A, Zeraouli Y. Thermal behavior of building material
containing microencapsulated PCM. Thermochim Acta 2012;550:427.
[132] Telkes, M. Thermal storage for solar heating and cooling. In: Proceeding of
the workshop on solar energy storage subsystems for the heating and cooling
of buildings, Charlottesville: Virginia, USA; 1975.
[133] Zhou D, Zhao CY, Tian Y. Review on thermal energy storage with phase
change materials (PCMs) in building applications. Appl Energy 2012;92:
593605.
[134] Kuznik F, Virgone J, Noel J. Optimization of a phase change material
wallboard for building use. Appl Therm Eng 2008;28:12918.
[135] Neeper DA. Thermal dynamics of wallboard with latent heat storage. Sol
Energy 2000;68(5):393403.
[136] Xiao W, Wang X, Zhang Y. Analytical optimization of interior PCM for energy
storage in a lightweight passive solar room. Appl Energy 2009;86:20138.
[137] Xu X, Zhang Y, Lin K, Di H, Yang R. Modeling and simulation on the thermal
performance of shape-stabilized phase change material oor used in passive
solar buildings. Energy Build 2005;37:108491.
[138] Zhou G, Zhang Y, Lin K, Xiao W. Thermal analysis of a direct-gain room with
shape-stabilized PCM plates. Renew Energy 2008;33:122836.
[139] Athienitis AK, Liu C, Hawes D, Banu D, Feldman D. Investigation of the
thermal performance of a passive solar test-room with wall latent heat
storage. Build Environ 1997;32:40510.
[140] Chen C, Guo H, Liu Y, Yue H, Wang C. A new kind of phase change material
(PCM) for energy-storing wallboard. Energy Build 2008;40:88290.
[141] Chen Z, Fang G. Preparation and heat transfer characteristics of microencap-
sulated phase change material slurry: a review. Renew Sustain Energy Rev
2011;15:462432.
[142] Delgado M, Lzaro A, Mazo J, Zalba B. Reviewon phase change material emulsions
and microencapsulated phase change material slurries: Materials, heat transfer
studies and applications. Renew Sustain Energy Rev 2012;16:25373.
[143] Youssef Z, Delahaye A, Huang L, Trinquet F, Fournaison L, Pollerberg C, et al.
State of the art on phase change material slurries. Energy Convers Manag
2013;65:12032.
[144] Charunyakorn P, Sengupta S, Roy SK. Forced convection heat transfer in
microencapsulated phase change material slurries: ow in circular ducts. Int
J Heat Mass Transf 1991;34:81933.
[145] Zhang Y, Faghri A. Analysis of forced convection heat transfer in micro-
encapsulated phase change material suspensions. J Thermophys Heat Transf
1995;9:72732.
[146] Goel M, Roy SK, Sengupta S. Laminar forced convection heat transfer in
microencapsulated phase change material suspensions. Int J Heat Mass
Transf 1994;37(4):593604.
[147] Alisetti EL, Roy SK. Forced convection heat transfer to phase change material
slurries in circular ducts. J Thermophys Heat Transf 2000;14:1158.
[148] Roy SK, Avanic BL. Laminar forced convection heat transfer with phase
change material suspensions. Int Commun Heat Mass Transf 2001;28(7):
895904.
[149] Hu X, Zhang Y. Novel insight and numerical analysis of convective heat
transfer enhancement with microencapsulated phase change material slur-
ries: laminar ow in a circular tube with constant heat ux. Int J Heat Mass
Transf 2002;45:316372.
[150] Zeng R, Wang X, Chen B, Zhang Y, Niu J, Wang X, et al. Heat transfer
characteristics of microencapsulated phase change material slurry in laminar
ow under constant heat ux. Appl Energy 2009;86:266170.
[151] Roy SK, Avanic BL. Turbulent heat transfer with phase change material
suspensions. Int J Heat Mass Transf 2001;44:227785.
[152] Royon L, Guiffant G. Forced convection heat transfer with slurry of phase
change material in circular ducts: a phenomenological approach. Energy
Convers Manag 2008;49:92832.
[153] Inaba H, Dai C, Horibe A. Numerical simulation of RayleighBnard convec-
tion in non-Newtonian phase change-material slurry. Int J Therm Sci
2003;42:47180.
[154] Inaba H, Zhang Y, Horibe A, Haruki N. Numerical simulation of natural
convection of latent heat phase-change-material microcapsulate slurry
packed in a horizontal rectangular enclosure heated from below and cooled
from above. Heat Mass Transf 2007;43:45970.

10 1016-j Rser 2014 02 032

Uploaded by

Copyright:

Available Formats

10 1016-j Rser 2014 02 032

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

10 1016-j Rser 2014 02 032

Uploaded by

Copyright:

Available Formats

Mathematical solutions and numerical models employed

for the investigations of PCMs' phase transformations

Time and cost consuming will be higher than a theoretical

The scales of the experiment test rig have to be considered, and

Test rig needs to be constructed and operated properly.

Proper environment parameters have to be simulated and

Relevant parameters need to be monitored, measuring appara-

Further investigation should be carried out to obtain adequate

Development of new numerical model to broaden its applica-

Assessment of the stability and convergence on the numerical

Unied dimensionless parameters needs to be developed to

Life cycle analysis, both economical and energetic feasibility of

You might also like