Applied Thermal Engineering 124 (2017) 1123–1133

Contents lists available at ScienceDirect

Applied Thermal Engineering

journal homepage: www.elsevier.com/locate/apthermeng

Research Paper

Melting dynamics of a phase change material (PCM) with dispersed

metallic nanoparticles using transport coefficients from empirical and
mean field models
Santiago Madruga ⇑, Gonzalo S. Mischlich
Department of Applied Mathematics to Aerospace Engineering, School of Aeronautics and Aerospace Engineering, Universidad Politécnica de Madrid (UPM), Plaza Cardenal Cisneros
3, 28040 Madrid, Spain

a r t i c l e i n f o a b s t r a c t

Article history: We study the melting process of n-octadecane with dispersed Al2 O3 nanoparticles in a semicircle. The
Received 28 September 2016 effective transport coefficients of the resulting nanofluid are modeled with (i) mean field models due
Revised 26 May 2017 to Maxwell-Garnett for the conductivity and Brinkmann for viscosity, and (ii) an empirical model based
Accepted 17 June 2017
on a least square fit to experimental data due to Corcione (2011). In both cases, we consider a uniform
Available online 20 June 2017
nanoparticle distribution in the liquid and solid phases and incorporate as well the change of conductiv-
ity in the latter phase. We carry out simulations with the transport coefficients predicted by both models
Keywords:
and find that Maxwell & Brinkmann overestimates heat transfer rates compared to the empirical fit for
Nanofluid
PCM
most of the ranges of nanoparticle concentration, size, and temperature. However, the proper selection
Convection of nanoparticles attending to their size and temperature can lead to enhanced heat transfer, even beyond
Nanoparticles of mean field model predictions. We show how the effective Prandtl number is the single most important
parameter that determines the dynamics and duration of the melting process, and how predictions of our
simulations agree with recent experiments (Ho and Gao, 2009).
Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction nanometers at small volume fractions, forming the named Nano-

enhanced PCM (NePCM). In contrast to micrometer size particles,
Phase Change Materials (PCM) exploit the latent heat of the they do not sediment, clog or produce considerable pressure drops
solid/liquid phase change to store a large amount of heat energy [18].
during melting or release it to the environment during solidifica- Heat transfer performance and dynamics of the resulting nano-
tion, hardly changing the temperature during the process. This sta- fluid (base fluid plus nanoparticles) depends critically on its
bility with respect to temperature changes and thermal storage effective thermal conductivity and viscosity. Because of that, a
capacity is key in many industrial applications of these materials large number of experiments have been carried out to measure
such as electronic cooling, air conditioning in buildings, waste heat these values in aqueous and non-aqueous fluids, specially on
recovery or to compensate the time offset between energy produc- Cu;Al2 O3 ; CuO and TiO2 nanoparticles [10,40].
tion and consumption in solar power plants. In addition, sophisticated theoretical models have been devel-
A major issue in thermal regulation with PCMs is the low con- oped to predict the conductivity for a given combination liquid/
ductivity typical of these materials, which leads to long times dur- nanoparticle at rest (Maxwell-Garnett, Hamilton-Crosser, Brugge-
ing the heat storage and release phases. Many solutions have been man, Jeffrey, Davis, etc.) and viscosity (Einstein, Brinkmann,
proposed to remedy this shortcoming by introducing high conduc- Batchelor, etc). These models tend to underestimate the value of
tivity metals within the PCM, often as metallic fins, or embedding these transport coefficients as a function of nanoparticle concen-
the PCM in a metallic porous foam [3,22]. Another promising trations, and they are not completely satisfactory even when
solution which holds high potential to increase the heat transfer micro-mixing effects due to Brownian motion are included [27].
rate, without adding complexity to the device, is the dispersion A different approach to determining the effective transport
in the PCM of metallic particles of sizes in the range of tens of coefficients is to fit available experimental data with a least square
regression. Corcione [12] uses experimental data of suspensions of
⇑ Corresponding author. a large set of metallic nanoparticles in water and ethylene glycol at
E-mail address: santiago.madruga@upm.es (S. Madruga). volume fractions ranging from 0.002 to 0.09 and temperatures

http://dx.doi.org/10.1016/j.applthermaleng.2017.06.097
1359-4311/Ó 2017 Elsevier Ltd. All rights reserved.
1124 S. Madruga, G.S. Mischlich / Applied Thermal Engineering 124 (2017) 1123–1133

from 294 K to 324 K to predict the thermal conductivity. Suspen- change materials. For the case of spherical and cylindrical geome-
sions in water, ethanol, ethylene glycol and propylene at volume tries, the voids break the mid-plane symmetry and importantly
fractions from 0.0001 to 0.071 and temperatures from 293 K to affect the heat transfer complicating the thermal boundary condi-
333 K are used to predict the viscosity. Although given the com- tions, since the upper free surface is open to air and the curved part
plexity of base fluid and nanoparticles interactions no current the- is in contact with the enclosure. In this work, a bidimensional
ories or single correlations can predict the broad spectrum of geometry in the form of a semicircle will be used, with a conduc-
experimental data, predictions from this fit are better than pro- tive curved part and adiabatic upper surface like happens in voids
posed theoretical models at a wide range of volume fractions. within spherical or cylindrical containers.
The potential of metallic nanoparticles in heat transfer applica- Most of the simulations on nanoparticle enhanced PCMs are
tions with PCMs has prompted the study of these systems during based on commercial CFD software and a lesser amount of in-
most of the last decade. Experiments with Cu nanoparticles of size house codes. However, Open Source code is an attractive long term
25 nm at 1 wt% on paraffin on a conductive state have found an option to compare and validate results on involved phase change
enhanced melting/solidification rates about 33% [43]. Ho and Gao simulations. From these considerations, we will use the Open
[19] measured the viscosity and conductivity of Al2 O3 nanoparti- Source code OpenFoam to carry out the simulations presented in
cles in n-octadecane and found the relative increase of effective this work.
viscosity was larger than that of thermal conductivity, casting From the whole previous discussion, we aim at studying the
doubt on the heat transfer enhancement of this nanofluid. Motahar melting of n-octadecane in half a disk with dispersed Al2 O3
et al. [30] studied as well the thermal conductivity and viscosity on nanoparticles of sizes 25 nm, 50 nm and 100 nm . We follow two
n-octadecane but with mesoporous silica particles of size about models to calculate the transport coefficients: (i) the mean field
300 nm at 1 wt.% and 5 wt.% showing a boost of non-Newtonian models due to Maxwell-Garnett for conductivity and Brinkmann
viscosity of 60% and conductivity only 5%. Jesumathy et al. [23] for viscosity, (ii) an empirical fit to experimental data due to
studied paraffin wax with 40 nm size CuO particles at 2, 5 and Corcione. To carry out the numerical simulations, we use the Open
10 wt.% in a complex annular heat storage system with a heat Source code OpenFoam based on the finite volume method.
transfer fluid flowing through the inner cylinder and found a In Section 2 we present the equations of the model, with a dis-
reduced melting time of 13, 21 and 35% for these volume fractions cussion of the models selected for the effective conductivity and
accompanied by a substantial decrease of latent heat. Reduction of viscosity. The next Section 3 explains the algorithms we follow
subcooling in water with Al2 O3 nanoparticles has been as well to solve the equations and validate our code with experimental
reported [42]. A recent experimental work [31] on the thermal and numerical results available in the literature of PCMs. The
behavior of a paraffin wax with dispersed alumina nanoparticles results of our simulations are presented in Section 4, where we
shows that not only volume fractions and conductivity of nanopar- observe how the most frequently used Maxwell-Garnett and
ticles determine the heat transfer capability of this nano-enhanced Brinkmann models underestimate the duration of the melting
PCM, but as well temperature and diameters of nanoparticles. The process compared to use transport coefficients derived from an
overall interplay between conduction, where metallic nanoparti- empirical fit in most of the parameter space. We show how the
cles always speed up the heat transfer, and convection, where vis- effective Prandtl number can be used to estimate the heat transfer
cosity and conductivity are competitive, determine the final heat performance of nanoparticle/liquified PCM combinations. This sec-
transfer performance of nano-enhanced PCMs. tion includes as well a comparison with a simpler conductive
Simulations of nanoparticle enhanced PCMs have been carried phase change process and the effect of volume fraction and size
out using a modified Maxwell-Garnett model for conductivity of nanoparticles in the melting process. Finally, conclusions of
and Brinkmann for viscosity in water with copper nanoparticles our work are provided in Section 5.
in a rectangle subjected to lateral heating [25] or within an enclo-
sure with lateral wavy surfaces [1]. Other materials like RT27 with 2. Basic equations and modeling
copper nanoparticles have been used inside a spherical container
filled up to 85% of capacity to simulate a void [20], paraffin wax The geometry consists of a semicircle of radius 2 cm (c.f. Fig. 1).
with copper nanoparticles in a square geometry [32], or n- The internal area (6.3 cm2) is within the order of other simulations
octadecane with CuO nanoparticles within annular, square and
cylinder geometries with a finite element method by Dhaidan
et al. [16,14,15].
Variations from Maxwell-Garnett and Brinkmann models have
been carried out replacing Brinkmann’s model for viscosity. For
instance, Kashani et al. [24] simulate copper nanoparticles dis-
persed in n-hexadecane in a rectangular geometry subjected to lat-
eral heating using Corcione’s empirical fit for viscosity. As well R
SOLID
Arasu and Mujumdar [5] uses an exponential model for the viscos-
ity to simulate the melting of a paraffin wax with Al2 O3 nanoparti-
cles in a square geometry subjected to lateral heating, and find a
faster melting rate for small volume fractions of the nanoparticles. y
The spherical geometry is specially interesting in micro-
encapsulated phase change materials, where the PCMs are encap- x
sulated in a polymeric spherical shell. It has been thoroughly
studied analytically, numerically and experimentally, both in melt-
ing and solidification processes [11,26,34,21]. However, to account
for the thermal expansion of the PCM, due to the difference in den- Fig. 1. Sketch of the computational domain. The inner semicircle represents the
solid phase of the PCM, and the region between it and the domain boundary the
sities between the liquid and solid phases, it is usual to leave a
liquid phase. The curved boundary is subjected to a constant temperature T h , and
space in the container generating unfilled spherical geometries the upper flat part of the boundary is adiabatic. Idealized convective cells appear as
[6,33]. For this same reason, voids appear in all encapsulated phase long arrows.
 S. Madruga, G.S. Mischlich / Applied Thermal Engineering 124 (2017) 1123–1133 1125

Table 1 keff represents the effective thermal conductivity. The enthalpy H

Thermophysical properties of n-octadecane (subscript s for the solid state and l for the
of the system comes from the contribution of the usual sensible
liquid state) and Al2 O3 nanoparticles (subscript np). Values for the solid and liquid
states are listed for the cases distinguished in the equations.
heat, due to changes of temperature in the solid and liquid phases
of the pcm and nanoparticles, and from the latent heat content,
q (kg m3 ) Ho and Gao [19] 776
which is the distinctive feature of PCM materials, i.e.
l (N s m2 ) Ho and Gao [19] 3:6  103 Z Z
T T
1 0 0
cs jcl (J kg K1 ) Alawadhi [4] 1934j2196
H ¼ Hs þ Hl þ Hnp ¼ ð1  /Þð1  f l Þqs cs ðT 0 ÞdT þ ð1  /Þf l ql cl ðT 0 ÞdT
ks jkl (W m1 K1 ) Ho and Gao [19] 0:358j0:13 T ref T ref

T s jT l (K) Ho and Gao [19] 298:25j299:65 Z T

0
L (J Kg1 ) Ho and Gao [19] 243:5  103 þ ð1  /Þf l ql L þ / qnp cnp ðT0ÞdT ð4Þ
T ref
1 4
a (K ) Dhaidan et al. [15] 9:1  10

knp (W m1 K1 ) Lin and Violi [28] 46 where cs (qs ), cl (ql ), cnp (qnp ) are the specific heats (densities) of the
1
cnp (J kg K1 ) Lin and Violi [28] 850 PCM in the solid phase, liquid phase, and nanoparticles, respec-
qnp (kg m3 ) Lin and Violi [28] 3900 tively. f l is the volume liquid fraction, and L is the latent heat of
anp (K1 ) Lin and Violi [28] 1:67  105 the solid/liquid phase change of the PCM. Notice that Eq. (4) pro-
vides a consistent thermodynamic formulation to add the effect of
nanoparticles on the enthalpy of the system. It can be readily gen-
eralized to nonuniform particle distributions and incorporates as
on PCMs [34,13,6]. The curved part of the domain is conductive and well the effect of the nanoparticles on the solid phase of PCM. On
held at a constant temperature T h ¼ 100  C, greater than the melt- supposing constant physical properties in each phase ql ¼ qS ¼ q,
ing temperature of the n-octadecane (T l ¼ 26:1  C), which is ini- the enthalpy Eq. (2) can be recast in an energy equation in terms
tially in the solid state (T i ¼ 25  C). The flat part of the boundary of the temperature
is adiabatic, behaving like a free surface without heat exchange
h
@ i
with the surroundings. The thermophysical properties of þ u  r qð1  /Þðð1  f l Þcs þ f l cl Þ þ qnp / cnp T
n-octadecane are listed in Table 1, @t
We consider the flow laminar, two-dimensional, and incom- @f l
pressible. For simplicity, viscous dissipation is neglected, nanofluid ¼ rðkeff rTÞ  ð1  /Þq L : ð5Þ
@t
considered a single phase, and nanoparticles and surrounding PCM
When T l  T s is a significant fraction of the range of tempera-
are supposed to be in local thermal equilibrium. The physical prop-
ture T h  T i , a relevant consideration in the form of the energy
erties are supposed to be constant within the range of tempera-
equation comes from the morphology of the solid/liquid interface,
tures studied, except the density in the buoyancy term, following
called the mushy region, as thoroughly discussed in Voller et al. [37]
the Boussinesq approximation. The origin of a Cartesian frame is
and shown numerically in Voller et al. [38]. At this scenario, the
placed in the center of the circular section of the domain. The gov-
limit case ul ¼ us (solid and liquid fully dispersed in the mushy
erning equations expressing the balance of momentum and energy,
region) corresponds to u ¼ ul ¼ us and leads to the appearance of
together with the continuity equation, read
an advective term for the latent heat in the energy equation; how-

 h i
@u ever, the opposed limit case us ¼ 0 (solidus fixed and liquid moving
qnf þ ðu  rÞu ¼ rp þ leff r2 u  g qnf  ðqaÞnf ðT  T ref Þ ey
@t around) corresponds to u ¼ f l ul and there is not such a term. How-
2 ever, for thin mushy regions, like in this work where
Cð1  f l Þ
 3
u; ð1Þ ðT l  T s Þ=ðT h  T i Þ  0:019 (c.f. Table 1) leads to mushy regions of
dþ fl just hundredths of micrometers, the advective term for the case

 of solid fully dispersed in the liquid phase is neglectable.
@H The last term in the momentum equation provides an empirical
þ r  ðu HÞ ¼ r  ðkeff rTÞ; ð2Þ
@t proportionality relationship, due to Darcy, between the pressure
gradient in a porous medium and the fluid velocity within it. The
and
functional form corresponds to the Carman-Kozeny equation and
r  u ¼ 0; ð3Þ is a convenient way to damp the velocity flow within the mushy
region and suppress it at the solid phase. We set in this term the
respectively. Where r ¼ ð@ x ; @ y Þ and @ t are the spatial and temporal Darcy coefficient C ¼ 1:6  106 kg m3 s, in compliance with previ-
operators, T is the volume-averaged temperature of a representative ous works [32,1], and d  1 a tiny constant to avoid division by
elementary volume that can contain solidified nanofluid, liquid zero without physical meaning. When the cell is completely liquid
nanofluid or a mixture of both phases in local thermal equilibrium; (f l ¼ 1) the Darcy term is null, like in a single phase fluid, when it is
u ¼ ðu; v Þ is the velocity, being u and v the horizontal and vertical completely solid (f l ¼ 0) the Darcy term becomes dominant, and
components; leff is the dynamic effective viscosity of the nanofluid, the velocity of the liquid becomes null. For intermediate values
p is the pressure, g is the magnitude of the gravity acceleration in of f l (0 < f l < 1), the PCM is in the mushy zone. Through this for-
the buoyancy term responsible of thermal convective motions, ey mulation, the Darcy term allows to use the momentum equation
a unit vector pointing in the vertical direction upwards, T ref is a ref- in the whole domain and model the phase change when fluid
erence temperature where physical properties are given. The den- motion is present without the complication of tracking the solid/
sity of the nanofluid qnf and the thermal expansion coefficient a liquid interface [35,37].
are averaged by the volume fraction of nanoparticles / according Finally, the latent heat contained in an elementary volume
to qnf ¼ ð1  /Þql þ / qnp and ðqaÞnf ¼ ð1  /Þðq aÞf þ /ðq aÞnp depends on the liquid fraction as DH ¼ f l q Lnf ¼ f l q ð1  /Þ L. Then
respectively, where the subscripts l and np refer to the liquid phase the coupling between the energy and momentum equation is given
of the PCM and the nanoparticles. The changes in density between through the liquid fraction f l , which in turn depends on the tem-
the liquid and solid phases are neglected. perature, the master variable of the phase change process. The liq-
The energy balance (2) has been written regarding the enthalpy uid fraction in the mushy zone is modeled by a linear relationship
H, with units of energy per volume, and the transport coefficient between the solidus and liquidus temperatures
1126 S. Madruga, G.S. Mischlich / Applied Thermal Engineering 124 (2017) 1123–1133

8
>
<0 if T 6 T s leff 1
¼ ð9Þ
f l ¼ DH=ðq Lnf Þ ¼ 1 if T P T l ð6Þ ll 1  34:87dnp =df 0:3 /1:03
>
:
ðT  T s Þ=ðT l  T s Þ if T s < T < T l
where dnp is the nanoparticle diameter, and df the equivalent diam-
DH ranges from 0 (PCM solid) to q Lnf (PCM liquid). The liquid and eter of a pure fluid molecule. This fit does not underestimate the
solid phases coexist for intermediate values. viscosity of nanofluids like Brinkmann’s mean field model does.
As for the thermal conductivity, the fit is more involved with depen-
2.1. Transport coefficients of the nanofluid dence on the Reynolds number Re, Prandtl number Pr and temper-
ature according to
2.1.1. Maxwell-Garnett and Brinkmann models  10  0:03
keff T knp
The models more frequently used to predict the thermal con- ¼ 1 þ 4:4 Re0:4 Pr0:66 /0:66 ð10Þ
kl Tl kf
ductivity and viscosity of nanoparticles dispersed in liquids are
the Maxwell-Garnett and Brinkmann models, respectively. Both More details on these expressions can be found in Corcione [12].
are mean field microscopic models, and a large number of works We notice that M&B is a very simplified model for the calcula-
on nanofluids have used these models, from now on M&B (Maxwell tion of transport coefficients on nanofluids. It neglects important
and Brinkmann), and are a reference to compare predictions with physical features like the size of nanoparticles, temperature, and
experimental results. properties of the base fluid. The empirical model LSC brings all
The Maxwell-Garnett model reads these physical parameters into consideration. Although other
physical effects are neglected, it narrows the gap between oversim-
keff knp þ 2kl  2/ðkl  knp Þ plified and overcomplex models incorporating available experi-
¼ ð7Þ
kl knp þ 2kl þ /ðkl  knp Þ mental quantities.

here keff ; kl and knp are the effective thermal conductivities of the
nanofluid at rest, the liquid phase of the PCM, and nanoparticles, 3. Numerical methods and validation of the model
respectively [29]. For particles embedded in solids, like the
nanoparticles in the solid phase of the PCM, this is a reliable model 3.1. Numerical methods
and will be used in this work. Notice how it does not depend on the
temperature or physical size of the nanoparticles. More complex To solve the model presented above, momentum equation mod-
formulations amend this model adding a term to include contribu- ified with a Darcy porous term coupled with an energy equation
tions to the thermal dispersion due to the Brownian motion. The including a source of latent heat, we use the open source software
value of the empirical factor of this coefficient is not stated explic- under GPL license OpenFoam [39]. It is written in C++ making pos-
itly in the PCM literature, and other factors like aggregation or sible a simple top-level description of partial differential equations
fluid/particle molecular interactions make this amended model still by operator overloading. Accompanying the basic suite a large col-
lacking to predict effective conductivities [27]. lection of tutorials and solvers applicable to a wide range of prob-
The Brinkmann model for the effective viscosity of the nano- lems in computational fluid dynamics in the laminar and turbulent
fluid treats the nanoparticles as rigid spheres and reads regimes are provided. It is designed with parallelization capability
for general 3D geometries and can be adapted straightly to 2D and
lf 1D geometries. OpenFoam discretizes partial derivative equations
leff ¼ ; ð8Þ
ð1  /Þ2:5 using the Finite Volume Method and has been thoroughly tested
and validated during more than a decade in solving mass, energy
where leff is the effective viscosity and lf the viscosity of the pure and momentum conservation equations. The conservation of fluxes
fluid [9]. This model is known to underestimate the viscosity of insured by the finite volume method makes this discretization spe-
nanofluids when compared with experiments. cially advantageous when nonlinear advective terms appear as in
convection dominated PCM dynamics.
2.1.2. Least-Square fit to experimental data The convective terms are discretized using a second order
We use an empirical correlation for the effective thermal con- upwind scheme. The time integration is carried out using a second
ductivity and viscosity developed recently by Corcione, from now order Crank-Nicolson scheme. The source term in the energy equa-
on LSC (Least-Square Corcione), which exhibits a good agreement tion requires particular attention since it couples the temperature
with Cu;Al2 O3 ; CuO and TiO2 nanoparticles dispersed in water, and liquid fraction fields strongly. Following Voller and Swami-
ethylene glycol and propylene [12]. nathan [36], it is linearized as a function depending on tempera-
There is not a single correlation able to predict the whole array ture, split in an explicit (zero order term) and implicit part (first
of experimental data, as mentioned in the introduction. However, order term), and the liquid fraction updated at every iteration.
it is enlightening to compare LSC predictions to experimental data A Rhie-Chow interpolation method is used to solve the enthalpy
from Ho and Gao [19] on Al2 O3 dispersed in n-octadecane at porosity equations to avoid checkerboard solutions. The momen-
weight fraction 5%. We obtain for the more sensitive thermal con- tum and continuity equation are solved using the PIMPLE algo-
ductivity keff ¼ 0:1326 for dnp ¼ 159:6 nm at T ¼ 30  C, and rithm, which ensures a right pressure-velocity coupling by
keff ¼ 0:137 at T ¼ 60  C. Both below 10% of discrepancy with combining SIMPLE and PISO algorithms [2]. The temperature equa-
respect to the experimental data reported by Ho and Gao. This dif- tion is solved for each PIMPLE outer-iteration, ensuring the conver-
ference is acceptable taking the perspective that Ho and Gao and gence of the velocity, pressure, and temperature fields. To improve
another recent set of experimental data by Motahar et al. [30] on convergence under-relaxation factors are used in the velocity,
pure n-octadecane show discrepancies about 10% in the thermal pressure, and temperature with values 0.7, 0.3 and 0.5 respectively.
conductivity. The good prediction of LSC is remarkable since it After the discretization and linearization of the equations, the
has been developed on metallic nanoparticles dispersed in low original PDE problem is transformed in a system of linear equa-
molecular weight liquids compared to the high value of paraffins. tions solved at each iteration using a multi-grid solver with toler-
The least-square fit to experimental data by Corcione leads to ance 108 for the pressure, velocity and temperature fields, and
the following prediction for the effective dynamic viscosity 106 for the liquid fraction.
 S. Madruga, G.S. Mischlich / Applied Thermal Engineering 124 (2017) 1123–1133 1127

3.2. Validation of the model a narrow parallelepiped enclosure of 63.5 mm height, 88.9 mm
length, and 38.1 mm depth is filled with gallium at a uniform inter-
We validate the stability and accuracy of the solver with a clas- nal temperature of 302.95 K. The left hot wall is kept at a constant
sical experimental setup [17], finite volume codes of commercial temperature T ¼ 311 K, the right cold at T ¼ 301:3 K, and horizon-
origin [13], in-house [41,32], and based on a finite difference con- tal top and bottom walls are adiabatic. The advance of the solid/
trol volume method [8]. liquid interface is shown in Fig. 2(b) at t ¼ 2; 6; 15, and 19 min.
A high quality structured hexahedral mesh created with the The quantitative agreement between our simulation and the
open-source mesher GMSH is used. The sharp corners between experiment is reasonable for all the times, capturing the regions
the curved and plane boundaries make the geometry very demand- of different slopes. Different between simulations and experimen-
ing to secure convergence and special attention has been put to tal data come from sources like sub-cooling, maintenance of
have a fine enough mesh to capture the high velocity and temper- constant temperature at the walls in experiments, or even 3D
ature gradients that appear near the boundaries and moving solid/ effects, which have been shown to play a role in the melting of gal-
liquid interface. lium even if the depth of the enclosure is small [7]. Fig. 2(b) shows
The dependence of the solution in the grid was tested from as well a comparison with the finite difference control volume code
coarse to fine grids using 10,000, 30,000, 50,000, 70,000, 90,000, of Brent et al. [8]. We obtain a better fit with experimental data at
110,000 and 130,000 cells up to leveling off the residuals of the the upper half part of the domain and later melting times, which is
1 probably due to the finer mesh allowed by current processing
norm L2 ð½0; 1Þ of the function f l -to ensure common support for
the curves- between consecutive cell sizes. Fig. 2(a) shows the evo- power. A comparison with Wittig and Nikrityuk [41] shows a more
lution of the overall liquid fraction in n-octadecane for the geome- stable front position in our numerical results.
try studied in this work at T h ¼ 100  C and initial temperature To validate the model and simulations with curved geometries
T i ¼ 25  C. Increasing the number of cells reduces the variation of with and without mid-plane symmetry we compare with Darzi
the melting time up to less than 0.5% from 90,000 cells progres- et al. [13]. They carry out simulations on melting of n-eicoseane
sively. A grid with this number of cells provides sufficient accuracy between two cylinders of diameter 20 mm and 40 mm in configu-
for the problems simulated in this work. ration concentric (circles in Fig. 3(a)) and eccentric (squares in 3
The code is as well validated with the benchmark experiment (a)) with 5 mm center to center distance. Darzi et al. carried out
on gallium melting by Gau and Viskanta [17]. In this experiment the simulations with the commercial code Fluent 6.0 using a grid

Fig. 2. (a) Evolution of the overall liquid fraction in n-octadecane for the geometry of Fig. 1 at T h ¼ 100  C for seven meshes with increasing cell count. (b) Evolution of the
solid/liquid interface as a function of time for the benchmark experiment on melting of gallium by Gau and Viskanta [17]. Up triangles correspond to the experimental data,
dashed line to simulations by Brent et al. [8], the dash-dotted line to the bidimensional case of Wittig and Nikrityuk [41], and the solid line to our simulations.

Fig. 3. (a) Evolution of the liquid fraction of n-eicosane as a function of time in an annulus. The solid and dashed lines correspond to our results. Circles to Array A of Darzi
et al. [13], and triangles to Array B of Darzi et al. [13]. (b) Evolution of the liquid fraction as a function of time in a square enclosure filled with paraffin wax and dispersed
copper nanoparticles. Lines correspond to our results and symbols to simulations of Sebti et al. [32]. / ¼ 0 corresponds to thick line (down triangles), / ¼ 0:025 to dashed
lines (squares), and / ¼ 0:05 to dash-dot line (up triangles).
1128 S. Madruga, G.S. Mischlich / Applied Thermal Engineering 124 (2017) 1123–1133

of 6470 cells. The agreement between our OpenFoam based model concentrations below Pr of the base fluid. In light of this only
(thick and dashed lines) is good for both geometries, particularly in enhanced heat transfer, at fixed temperature, is expected within
the linear regime at the beginning of melting and sublinear at the a narrow range of /. In short, from Fig. 4 one expects faster melting
final stages. time using M&B at most of the nanoparticle concentrations com-
Finally, we compare our results with the simulation of a PCM pared to LSC, leading to an overestimation of the heat transfer per-
enhanced with metallic nanoparticles carried out by Sebti et al. formance in PCMs with nanoparticles compared to that predicted
[32] with an in-house fortran finite volume solver following as well by an empirical fit of the transport coefficients.
an enthalpy porosity formulation. They used a Brinkman model for
viscosity and Maxwell-Garnett for the conductivity. The enclosure 4.1. Conductive transport
is a square cavity of length 10 mm filled with paraffin wax and dis-
persed copper nanoparticles at volume fractions 0.025 and 0.5. The A reference to the overall effect of natural convection in the
initial temperature is T i ¼ 301:3 K, the left hot wall is held at dynamics of PCM melting is provided taking into account first only
T c ¼ 306:3 K and the rest of walls are adiabatic. Our results showed transport by conduction, i.e. using only Eq. (5). Evolution of the
in Fig. 3(b) superpose to Sebti et al., indicating an excellent fit dur- overall liquid fraction as a function of time is shown in Fig. 5,
ing the whole melting process at all volume fractions. Notice that where the conductive hot curved wall is held at T h ¼ 100  C and
simulations have been validated from paraffin waxes, which exhi- an initial temperature of 25  C is used. Inclusion of nanoparticles
bit a slow conductive and convective transport, to the liquid metal modifies the conductivity of the solid and stagnant melted part,
gallium with very high heat transport rate. but there is not particle diffusion and a Maxwell-Garnett model
is suitable to calculate the effective thermal conductivity. The
4. Results and discussion melting time for the base fluid is 3442 s, nanoparticles at a volume
fraction / ¼ 0:01 reduce it about 3.9% to 3308 s, 19.4% for / ¼ 0:06,
The addition of metallic nanoparticles increases the effective and 32.2% for / ¼ 0:11. This improved heat storage rate is expected
thermal conductivity and viscosity of the base fluid in convection. due to the much higher conductivity of the metallic nanoparticles.
These two effects are competitive as for the final heat transfer The lower conductivity of the liquid phase with respect to the
capability of the nanofluid. If the first is dominant a faster melting solid phase leads to a quick advance of the overall liquid fraction
of the PCM will be promoted with nanoparticles, if the second is in the first states of melting, requiring about 85% of the melting
dominant then a slower melting is expected. The relative strength time to melt the last 50% of the PCM where the liquid phase is
of both effects can be evaluated using the effective Prandtl number dominant. The slowing down of the melting rate is observed as a
Preff ¼ meff =jeff , where meff is the effective kinematic viscosity and loss of linearity of the overall liquid fraction curve at f l ¼ 0:5,
jeff the effective thermal diffusivity. where a sub-linear regime appears at a log-log scale. This reduction
Features of the nanofluid, such as small nanoparticle concentra- of the melting time under conduction has been observed as well in
tion, small particle diameter or high temperature promote the experiments of Cu nanoparticles in paraffin, although physical
underestimation of conductivity by Maxwell-Garnett with respect properties of the used paraffin are not provided [43]. At larger par-
to LSC. Opposite features favor the inverse situation. On the other ticle concentrations / ¼ 0:06 and 0:11 we observe a similar
hand, Brinkman’s model underestimates the viscosity for all advance of the overall liquid fraction, scaling with f l ðtÞ  t 0:47 ,
nanoparticle sizes, temperatures, and nanoparticle concentrations and a faster decay in the final sub-linear regime.
with respect to LSC, being this effect more pronounced for small
particles. Thus LSC predicts a lower effective Prandtl number
4.2. Convective transport
(and increased heat transfer rate) than M&B for very small
nanoparticles at high temperatures and low concentrations. For
We present in this section the results including as well the con-
the rest of the cases, M&B predicts a lower effective Prandtl num-
vective motion driven by buoyancy. The impact on the heat trans-
ber. This behavior is illustrated in Fig. 4 for selected temperatures
fer process of the nanofluid is huge since overall the convective
and nanoparticle sizes of Al2 O3 dispersed in n-octadecane. M&B
motion decrease the melting time by a factor ten. Not surprisingly,
shows a monotonic decrease of Pr eff with the volume fraction, thus
predicting a continuous improvement of the heat transfer rate of
the nanofluid with increasing /. However, LSC shows a monotonic
increase of Pr eff from either / ¼ 0 or a local minimum at low

Fig. 5. Evolution of the liquid fraction as a function of time in log-log scale for
n-octadecane when only conductive heat transport is considered. Effective
conductivity of the phases solid and quiescent liquid are calculated using the
Fig. 4. Effective Prandtl number as a function of the nanoparticle volume fraction / Maxwell-Garnett model. Conductive curved boundary is held at T h ¼ 100  C. (i)
predicted by Corcione’s least-square empirical fit (symbols) and Maxwell&Brink- (thick line) / ¼ 0:0, (ii) (dashed line) / ¼ 0:01, (iii) (dot-dashed line) / ¼ 0:06, (iv)
mann (solid line) models. (dotted line) / ¼ 0:11
 S. Madruga, G.S. Mischlich / Applied Thermal Engineering 124 (2017) 1123–1133 1129

favoring convective motion is frequently used to reduce the heat

storage and release times in applications of PCMs.
We can calculate an effective Rayleigh number during the melt-
3
ing in the region below the solidus Raeff ¼ g a DTh =ðjeff meff Þ, where
DT ¼ T h  T i and h is the thickness between the boundary at the
lowest part of the domain and the solidus along the vertical sym-
metry axis. A maximum Raeff can be estimated when h reaches
the half of the depth along the vertical symmetry axis. Thus for
pure n-octadecane Ra extends up to 1:9  106 . Increasing /, M&B
shrinks slightly the range up to Raeff  1:8 for / ¼ 0:06, and LSC
more strongly to Raeff ¼ 3:3  105 for / ¼ 0:06.
The two main regimes observed in the overall liquid fraction
curve can be derived from the classic Rayleigh-Bénard theory of
convection in a horizontal liquid layer heated from below. When
h is small Raeff is small and the viscous and thermal dissipation Fig. 7. Evolution of the liquid fraction of n-octadecane for T h ¼ 100  C. (i) (solid
line) base fluid, (ii) (dashed line) M&B model with / ¼ 0:01, (iii) LSC model with
stronger than the destabilizing buoyancy and heat transport takes
/ ¼ 0:01 for dnp ¼ 25 nm (circles) and dnp ¼ 100 nm (squares).
place by conduction. This first conductive regime appears as an ini-
tial linear growth of the Rayleigh number as a function of time in
the log-log scale of Fig. 6 (up to 40 s at / ¼ 0). As the solid/liquid
interface advances Raeff increases and the buoyancy force becomes
dominant, leading to the destabilization of the quiescent liquid at
the bottom of the disk and the appearance of a layer of convective
cells. As a reference, for a horizontal layer of n-octadecane bounded
by rigid plates, the destabilization takes places when Ra reaches
the critical value 1708 at a thickness of 1.5 mm. This second con-
vection dominated regime is faster, starts at the end of the first lin-
ear growth stage of Fig. 6 (from 40 s at / ¼ 0), and leads to an
almost linear growth of the overall liquid fraction curve up to the
latest phase of PCM melting at 85% of the overall liquid fraction.
There the strong shrinkage of the more conductive solid phase
declines the slope of the overall liquid fraction curve (c.f. Figs. 7
and 8). Within the convection dominated regime, there are addi-
tional sub-regimes that depend on specifics physical parameters.
The other primary mechanism behind the dynamics of melting Fig. 8. Evolution of the liquid fraction of n-octadecane for T h ¼ 100  C. (i) (solid
is the lateral heating dominant nearby the upper free surface due line) base fluid, (ii) (dashed line) M&B model with / ¼ 0:06, (iii) (squares) LSC
to the symmetry of the disk heated from the curved boundary. model for / ¼ 0:06 and dnp ¼ 50 nm, (iv) (triangles down) LSC model for / ¼ 0:06
and dnp ¼ 25 nm.
Interestingly, from the zone between lateral and vertical dominant
cold plumes emerge from the colder solid side and may lead to
strong temperature gradients near the conductive wall. Lateral 4.2.1. Volume fraction / ¼ 0:01
heating leads quickly to the formation of a single large convective We first carry out simulations in the region where M&B predicts
cell with liquid moving upwards adjacent to the conductive wall, larger Preff than LSC. Fig. 7 exhibits the evolution of the overall liq-
roughly parallel to the free surface and falling near the colder uid fraction of n-octadecane during melting with T h ¼ 100  C. M&B
solid/liquid interface. Since fluid velocity is stronger near the upper provides a single curve, independent of the nanoparticle size. As
flat surface there is a more substantial advance of the solid/liquid expected, the melting time is the longest in the absence of metallic
interface at this region that breaks the circular symmetry of the nanoparticles (390 s). The inclusion of Al2 O3 nanoparticles at
solidus (c.f. Figs. 9 and 11). / ¼ 0:01 leads to a small decrease of the melting time using

Fig. 6. Dy3 (space dependence of the Rayleigh number) as a function of the time measured across the vertical axis of symmetry of the semicircle. The solid line corresponds to
LSC for dnp ¼ 50 nm, and the dashed-lines with symbols to M&B.
1130 S. Madruga, G.S. Mischlich / Applied Thermal Engineering 124 (2017) 1123–1133

Fig. 9. Snapshots at t ¼ 100 s; 200 s and 300 s of the temperature field (third and fourth rows) and streamlines (first and second rows) for volume fraction / ¼ 0:01. First and
third rows exhibit the results for the transport coefficients calculated according to LSC for dnp ¼ 50 nm, and second and fourth to M&B.

M&B (380 s). Instead, if the empirical LSC is used to calculate the At t ¼ 300 s different dynamics between M&B and LSC begins to
transport coefficients, we find a larger reduction of the melting appear. LSC shows a reconnection of the upper layer of cells merg-
time: 372 s, 368 s and 360 s for particles diameters of 100 nm, ing as well with the large lateral cell, leaving just a single pair of
50 nm, and 25 nm nanometers, respectively. For the smallest par- small cells at the bottom split from the general circulation. The
ticles, dnp ¼ 25 nm, LSC leads to a melting time reduction of 5.6% large laterally dominated flow of the upper part leads to a plume
with respect to M&B. of cold fluid descending through the symmetry axis of the domain
A sequence of snapshots for streamlines and temperature fields up to touching the hot plume rising from the bottom. Meanwhile,
at t ¼ 100; 200 and 300 s for M&B and LSC is shown at Fig. 9. At M&B shows a state where reconnection of cells has not happened.
t ¼ 100 s three pairs of counter-rotating convective cells are gener- This makes M&B more inefficient in the heat transport with respect
ated at the bottom of the disk after a Rayleigh-Bénard instability at to LSC where large convective cells carry the heat from the hot wall
the onset of convection. The snapshots of the temperature field to the colder solidus. The effect of this process is the reduced size of
show a hot plume wherever a pair of counter-rotating cells moves the solidus in the LSC model compared to M&B.
upwards. This hotter ascending liquid leads to a deeper advance of Interestingly, the overall dynamics of the melting process can
the solid/liquid interface at these regions generating three concave be tracked from the motion of cold and hot plumes which can be
parts at the bottom of the solidus boundary. located through the local Nusselt number, defined as the local
As melting continues at t ¼ 200 s the dynamics of the LSC and dimensionless heat flux:
M&B is qualitatively the same. The increased gap of the melted part
at the bottom gives room to the creation of two superposed layers R
NuðhÞ ¼ jeff ðhÞn  rT: ð11Þ
of convective cells. This happens through the shift of the outer jf DT
pairs of cells at 100 s and their superposition over the central pair.
All that is accompanied by an increase in the size of all cells and is where R is the semicircle radius, DT is the difference between the
observed as well in the temperature field, where the central plume melting temperature and T h ; jf is the n-octadecane thermal con-
becomes caught between the ones outside. The superposition of ductivity, jeff is the effective thermal conductivity at the wall, and
the convective cells is as well observed in Fig. 6 from 110 s to n is the outward unit vector normal to the boundary. Fig. 10 exhi-
 40 s with a slowdown of the speed of propagation of the solid/ bits the variation of the local Nusselt number at the conductive
liquid front. boundary as a function of the angle in polar coordinates. At
 S. Madruga, G.S. Mischlich / Applied Thermal Engineering 124 (2017) 1123–1133 1131

Fig. 10. Local Nusselt number at / ¼ 0:01 for t ¼ 100 s; 200 s and 300 s calculated at the curved boundary of the domain as a function of the angle in polar coordinates. Solid
line corresponds to M&B, line with symbols (circles) to LSC for dnp ¼ 50 nm.

t ¼ 100 s there is an absolute maximum at h  3 p=8 and a local at this volume fraction LSC predicts a surge to Preff ¼ 72:6 for
minimum closer to the vertical axis. The position of the maximum dnp ¼ 50 nm, and Preff ¼ 97 for dnp ¼ 25 nm. Thus a decrease in
corresponds to plumes of cold liquid approaching the hot wall, the heat transfer rate is expected for the latter model. This is con-
and minimum corresponds to hot plumes emerging from the wall firmed in Fig. 8, where the liquid fraction for M&B is above LSC
(c.f. third and fourth rows of Fig. 9). At t ¼ 200 s the maximum once convection sets in. Furthermore, the melting time for M&B
and minimum shift towards lower angles, but overall curve is the is 345:1 s, while for LSC is 371:4 s (dnp ¼ 50 nm) and 406 s
same. However, at latest times t ¼ 300 s local Nusselt curves for (dnp ¼ 25 nm). This scenario where effective Prandtl number is
M&B and LSC begin to separate as expected from the very different lower for M&B than LSC is the most frequent situation in simula-
temperature field of Fig. 9. tions and experiments, and that leads to an overestimation of the
melting rate of M&B, as shown in the liquid fraction curves.
4.2.2. Volume fraction / ¼ 0:06 A sequence of snapshots at t ¼ 100; 200 and 300 s for streamli-
Prandtl number for base n-octadecane is 60.8 according to nes and temperature fields for M&B and LSC is displayed at Fig. 11.
Fig. 4. M&B predicts a decline to Pr eff ¼ 50:8 for / ¼ 0:06. However, At t ¼ 100 s only a pair of well formed counter-rotating convective

Fig. 11. Snapshots at t ¼ 100 s; 200 s and 300 s of the temperature field (third and fourth rows) and streamlines (first and second rows) for volume fraction / ¼ 0:06. First and
third rows exhibit the results for the transport coefficients calculated according to LSC for dnp ¼ 50 nm, and second and fourth for M&B.
1132 S. Madruga, G.S. Mischlich / Applied Thermal Engineering 124 (2017) 1123–1133

Fig. 12. Local Nusselt number at / ¼ 0:06 for t ¼ 100 s; 200 s and 300 s calculated at the curved boundary of the domain as a function of the angle in polar coordinates. Solid
line corresponds to M&B, line with symbols (circles) to LSC for dnp ¼ 50 nm.

cells are generated at the bottom of the domain for LSC, mean- 5. Conclusions
while, for M&B we have three pairs of convective cells (and three
ascending hot plumes) and more advanced solid/liquid interface, We have carried out a numerical study on melting dynamics of
showing the enhanced heat transfer rate predicted by M&B for the PCM n-octadecane with Al2 O3 nanoparticles uniformly dis-
/ ¼ 0:06, and decreased for LSC with respect to / ¼ 0:01. persed at small concentrations. A semicircular geometry of radius
At further times 200 and 300 s, LSC overall circulation continues 2 cm has been selected, whose curved part is held at constant tem-
to be the same, with a single hot plume arising at the bottom. How- perature and free surface is adiabatic. The phase change problem,
ever, M&B simulations evolve as LSC at / ¼ 0:01, with cold and hot with the complication of a moving solid/liquid interface, has been
plumes touching in the symmetry axis. This reflects the close Preff modeled using an enthalpy-porosity formulation. The numerical
exhibited in both cases. Comparing solidus of both of LSC and M&B procedure is based on finite volumes using the Open Source soft-
models for / ¼ 0:06 at 300 s shows a minor area of solidus left for ware OpenFoam.
M&B. The effect of metallic nanoparticles on the solid phase of the
As in case / ¼ 0:01, the overall dynamics of the melting process PCM has been modeled via a Maxwell-Garnett relation for the con-
with / ¼ 0:06 can be tracked from the motion of cold and hot ductivity. The enhanced conductivity of the solid phase due to the
plumes using the local Nusselt number. In the latest case, the local metallic nanoparticles has been shown to be important in the
Nusselt curves (Fig. 12) for M&B and LSC differ for all shown times reduction of the overall melting time of the PCM. The modeling
because of the different configuration of convective cells resulting of the effective conductivity and viscosity of the liquid phase of
from M&B and LSC model. the PCM has been carried out through mean field models,
Fig. 4 shows that final melting time of the PCM depends on the Maxwell-Garnett for conductivity and Brinkmann for viscosity,
interplay between conductive and convective processes. The melt- and an empirical model based on a fit to experimental data pro-
ing time with LSC for / ¼ 0:06 and dnp ¼ 50 nm is slightly shorter posed by Corcione. Predictions of both models are compared in
than without nanoparticles in spite of the increased effective the article. Nanoparticles of small size, high temperature and small
Prandtl number for LSC. The first factor explaining this behavior volume fractions favor a lower effective Prandtl number -and
is the enhanced conductivity of the solid phase by the presence enhanced rate of melting- for Corcione’s, the opposite situation
of metallic nanoparticles leading to a faster heat transfer rate while promotes a lower effective Prandtl number for M&B. The most fre-
the conduction dominates the melting process, up to f l  0:2 in quent cases studied numerically in the literature are found within
Fig. 8. In this regime, the overall liquid fraction curve of the base the favorable M&B regime, which leads to an overestimation of the
fluid is below the cases with nanoparticles. The second factor deals reduction of melting times for the convective regime compared
with the configuration of convective cells and plumes at the latter with the empirical model due to Corcione. However, we have
stages of melting, which determines the evolution of the very small shown that is possible to choose the size of alumina nanoparticles,
fraction of solidus area at this stage. The LSC model exhibits two volume fractions, and temperatures where the NePCM exhibits an
convective cells at the bottom, one on each side of the symmetry improved heat transfer rate with respect to what is predicted by
axis generating a reduced viscous dissipation compared to more the simpler M&B model. In these situations, NePCM become even
involved configurations, and a hot plume in the symmetry axis more attractive for thermal management and energy storage
enhancing the heat transfer in the zones nearby. This improvement applications.
in the heat transfer along the symmetry axis is seen as well in Fig. 6 Our results agree with experiments carried out by Ho and Gao
(b), where the lack of the transition regime of superimposed cells in n-octadecane with dispersed Al2 O3 nanoparticles, where no
for LSC with / ¼ 0:06 accelerates the motion of the solid/liquid improvement of the heat transfer rate in the convective regime
front. was forecasted based upon measures of effective conductivity
This analysis illustrates that while the effective Prandtl number and viscosity.
has the strongest influence in the melting process (as observed in There exists no universal fit to experimental data of the effec-
LSC for dnp ¼ 25 nm during the whole melting), one must take into tive conductivity and viscosity and LSC is more likely not the best
account more factors to explain the melting dynamics; such as the possible option. However, our results show that predictions of heat
duration and extension of conductive dominated regime, or the transfer rates of nanofluid based PCMs with classic mean field
state at the latter phase of melting. models like M&B are generally overestimated, and a careful
 S. Madruga, G.S. Mischlich / Applied Thermal Engineering 124 (2017) 1123–1133 1133

combination of fluid and nanoparticles must be selected to [18] L. Godson, B. Raja, D. Mohan Lal, S. Wongwises, Enhancement of heat transfer
using nanofluids-an overview, Renew. Sustain. Energy Rev. 14 (2) (2010) 629–
improve the heat transfer rate.
641.
[19] C.J. Ho, J.Y. Gao, Preparation and thermophysical properties of nanoparticle-in-
Acknowledgements paraffin emulsion as phase change material, Int. Commun. Heat Mass Transfer
36 (5) (2009) 467–470.
[20] S.F. Hosseinizadeh, A.A.R. Darzi, F.L. Tan, Numerical investigations of
Santiago Madruga acknowledges support by Erasmus Mundus unconstrained melting of nano-enhanced phase change material (NEPCM)
EASED programme (Grant 2012-5538/004-001) coordinated by inside a spherical container, Int. J. Therm. Sci. 51 (1) (2012) 77–83.
Centrale Supélec, the Spanish Ministerio de Economía y Competitivi- [21] S.F. Hosseinizadeh, A.A. Rabienataj Darzi, F.L. Tan, J.M. Khodadadi,
Unconstrained melting inside a sphere, Int. J. Therm. Sci. 63 (2013) 55–64.
dad under Projects No. TRA2016-75075-R, No. ESP2013-45432-P [22] S. Jegadheeswaran, S.D. Pohekar, Performance enhancement in latent heat
and No. ESP2015-70458-P, and the computer resources and thermal storage system: a review, Renew. Sustain. Energy Rev. 13 (9) (2009)
technical assistance provided by the Centro de Supercomputación 2225–2244.
[23] S. Jesumathy, M. Udayakumar, S. Suresh, Experimental study of enhanced heat
y Visualización de Madrid (CeSViMa). transfer by addition of CuO nanoparticle, Heat Mass Transfer 48 (6) (2012)
965–978.
References [24] S. Kashani, A.A. Ranjbar, M.M. Madani, M. Mastiani, H. Jalaly, Numerical study
of solidification of a nano-enhanced phase change material (NEPCM) in a
thermal storage system, J. Appl. Mech. Tech. Phys. 54 (5) (2013) 702–712.
[1] M. Abdollahzadeh, M. Esmaeilpour, Enhancement of phase change material
[25] J.M. Khodadadi, S.F. Hosseinizadeh, Nanoparticle-enhanced phase change
(PCM) based latent heat storage system with nano fluid and wavy surface, Int.
materials (NEPCM) with great potential for improved thermal energy
J. Heat Mass Transfer 80 (2015) 376–385.
storage, Int. Commun. Heat Mass Transfer 34 (5) (2007) 534–543.
[2] H.J. Aguerre, S. Marquez Damian, J.M. Gimenez, N.M. Nigro, Modeling of
[26] J.M. Khodadadi, Y. Zhang, Effects of buoyancy-driven convection on melting
compressible fluid problems with openfoam using dynamic mesh technology,
within spherical containers, Int. J. Heat Mass Transfer 44 (8) (2001) 1605–
Mecánica Comput. XXXII (2013) 995–1011.
1618.
[3] F. Agyenim, P. Eames, M. Smyth, A comparison of heat transfer enhancement in
[27] C. Kleinstreuer, Y. Feng, Experimental and theoretical studies of nanofluid
a medium temperature thermal energy storage heat exchanger using fins, Sol.
thermal conductivity enhancement: a review, Nanoscale Res. Lett. 6 (1) (2011)
Energy 83 (9) (2009) 1509–1520.
439.
[4] E.M. Alawadhi, Thermal analysis of a building brick containing phase change
[28] K.C. Lin, A. Violi, Natural convection heat transfer of nanofluids in a vertical
material, Energy Build. 40 (3) (2008) 351–357.
cavity: Effects of non-uniform particle diameter and temperature on thermal
[5] A.V. Arasu, A.S. Mujumdar, Numerical study on melting of paraffin wax with
conductivity, Int. J. Heat Fluid Flow 31 (2) (2010) 236–245.
Al2O3 in a square enclosure, Int. Commun. Heat Mass Transfer 39 (1) (2012) 8–
[29] J. Maxwell, A Treatise on Electricity and Magnetism, third ed., Dover
16.
Publications, 1954.
[6] E. Assis, G. Ziskind, R. Letan, Numerical and experimental study of
[30] S. Motahar, N. Nikkam, A.A. Alemrajabi, R. Khodabandeh, M.S. Toprak, M.
solidification in a spherical shell, J. Heat Transfer 131 (2) (2009) 024502.
Muhammed, A novel phase change material containing mesoporous silica
[7] O. Ben-David, A. Levy, B. Mikhailovich, A. Azulay, 3D numerical and
nanoparticles for thermal storage: a study on thermal conductivity and
experimental study of gallium melting in a rectangular container, Int. J. Heat
viscosity, Int. Commun. Heat Mass Transfer 56 (2014) 114–120.
Mass Transfer 67 (2013) 260–271.
[31] M. Nourani, N. Hamdami, J. Keramat, A. Moheb, M. Shahedi, Thermal behavior
[8] A. Brent, V. Voller, K. Reid, Enthalpy-porosity technique for modelling
of paraffin-nano-Al2O3 stabilized by sodium stearoyl lactylate as a stable
convection-diffusion phase change: application to the melting of a pure
phase change material with high thermal conductivity, Renew. Energy 88
metal, Numer. Heat Transfer 13 (1988) 297–318.
(2016) 474–482.
[9] H.C. Brinkman, The viscosity of concentrated suspensions and solutions, J.
[32] S.S. Sebti, M. Mastiani, H. Mirzaei, A. Dadvand, S. Kashani, S.A. Hosseini,
Chem. Phys. 20 (1952), pp. 571–571.
Numerical study of the melting of nano-enhanced phase change material in a
[10] J. Buongiorno, D.C. Venerus, N. Prabhat, T. McKrell, J. Townsend, R.
square cavity, J. Zhejiang Univ. Sci. A 14 (5) (2013) 307–316.
Christianson, Y.V. Tolmachev, P. Keblinski, L.W. Hu, J.L. Alvarado, I.C. Bang, S.
[33] L. Solomon, A.F. Elmozughi, A. Oztekin, S. Neti, Effect of internal void
W. Bishnoi, M. Bonetti, F. Botz, A. Cecere, Y. Chang, G. Chen, H. Chen, S.J. Chung,
placement on the heat transfer performance Encapsulated phase change
M.K. Chyu, S.K. Das, R. Di Paola, Y. Ding, F. Dubois, G. Dzido, J. Eapen, W. Escher,
material for energy storage, Renew. Energy 78 (2015) 438–447.
D. Funfschilling, Q. Galand, J. Gao, P.E. Gharagozloo, K.E. Goodson, J.G.
[34] F.L. Tan, S.F. Hosseinizadeh, J.M. Khodadadi, L. Fan, Experimental and
Gutierrez, H. Hong, M. Horton, K.S. Hwang, C.S. Iorio, S.P. Jang, A.B. Jarzebski,
computational study of constrained melting of phase change materials
Y. Jiang, L. Jin, S. Kabelac, A. Kamath, M.A. Kedzierski, L.G. Kieng, C. Kim, J.H.
(PCM) inside a spherical capsule, Int. J. Heat Mass Transfer 52 (15–16)
Kim, S. Kim, S.H. Lee, K.C. Leong, I. Manna, B. Michel, R. Ni, H.E. Patel, J. Philip,
(2009) 3464–3472.
D. Poulikakos, C. Reynaud, R. Savino, P.K. Singh, P. Song, T. Sundararajan, E.
[35] V. Voller, C. Prakash, A fixed grid numerical modelling methodology for
Timofeeva, T. Tritcak, A.N. Turanov, S. Van Vaerenbergh, D. Wen, S. Witharana,
convection-diffusion mushy region phase-change problems, Int. J. Heat Mass
C. Yang, W.H. Yeh, X.Z. Zhao, S.Q. Zhou, A benchmark study on the thermal
Transfer 30 (8) (1987) 1709–1719.
conductivity of nanofluids, J. Appl. Phys. 106 (9) (2009).
[36] V. Voller, C. Swaminathan, General source–based method for solidification
[11] J. Caldwell, C.-c. Chan, Spherical solidi Ò cation by the enthalpy method and the
phase change, Numer. Heat Transfer Part B 19 (1991) 175–189.
heat balance integral method, Appl. Math. Model. 24 (2000) 45–53.
[37] V.R. Voller, A.D. Brent, C. Prakash, The modelling of heat, mass and solute
[12] M. Corcione, Empirical correlating equations for predicting the effective
transport in solidification systems, Int. J. Heat Mass Transfer 32 (9) (1989)
thermal conductivity and dynamic viscosity of nanofluids, Energy Convers.
1719–1731.
Manage. 52 (1) (2011) 789–793.
[38] V.R. Voller, A.D. Brent, C. Prakash, Modelling the mushy region in a binary
[13] A.R. Darzi, M. Farhadi, K. Sedighi, Numerical study of melting inside concentric
alloy, Appl. Math. Model. 14 (6) (1990) 320–326.
and eccentric horizontal annulus, Appl. Math. Model. 36 (9) (2012) 4080–
[39] H.G. Weller, G. Tabor, A tensorial approach to computational continuum
4086.
mechanics using object-oriented techniques, Comput. Phys. 12 (6) (1998)
[14] N.S. Dhaidan, J.M. Khodadadi, T.A. Al-Hattab, S.M. Al-Mashat, Experimental
620–631.
and numerical investigation of melting of NePCM inside an annular container
[40] D. Wen, G. Lin, S. Vafaei, K. Zhang, Review of nanofluids for heat transfer
under a constant heat flux including the effect of eccentricity, Int. J. Heat Mass
applications, Particuology 7 (2) (2009) 141–150.
Transfer 67 (2013) 455–468.
[41] K. Wittig, P.a. Nikrityuk, Three-dimensionality of fluid flow in the benchmark
[15] N.S. Dhaidan, J.M. Khodadadi, T.a. Al-Hattab, S.M. Al-Mashat, Experimental and
experiment for a pure metal melting on a vertical wall, IOP Conf. Ser. Mater.
numerical investigation of melting of phase change material/nanoparticle
Sci. Eng. 27 (2012) 012054.
suspensions in a square container subjected to a constant heat flux, Int. J. Heat
[42] S. Wu, D. Zhu, X. Li, H. Li, J. Lei, Thermal energy storage behavior of Al2O3-H2O
Mass Transfer 66 (2013) 672–683.
nanofluids, Thermochim. Acta 483 (1-2) (2009) 73–77.
[16] N.S. Dhaidan, J.M. Khodadadi, T.a. Al-Hattab, S.M. Al-Mashat, Experimental and
[43] S.Y. Wu, H. Wang, S. Xiao, D.S. Zhu, An investigation of melting/freezing
numerical study of constrained melting of n-octadecane with CuO
characteristics of nanoparticle-enhanced phase change materials, J. Therm.
nanoparticle dispersions in a horizontal cylindrical capsule subjected to a
Anal. Calorim. 110 (3) (2012) 1127–1131.
constant heat flux, Int. J. Heat Mass Transfer 67 (2013) 523–534.
[17] C. Gau, R. Viskanta, Melting and solidification of a pure metal on a verticall
wall, J. Heat Transfer. 108 (1) (1986) 174–181.