Eur. Phys. J.

Special Topics 226, 1169–1176

© EDP Sciences, Springer-Verlag 2017 THE EUROPEAN
DOI: 10.1140/epjst/e2016-60207-7 PHYSICAL JOURNAL
SPECIAL TOPICS
Regular Article

Enhancement of heat transfer rate on phase

change materials with thermocapillary flows
Santiago Madruga1,a and Carolina Mendoza2
1
Department of Applied Mathematics to the Aerospace Engineering, School of Aerospace
Engineering, Universidad Politécnica de Madrid (UPM), Plaza Cardenal Cisneros 3,
28040 Madrid, Spain
2
Instituto de Ciencias Matemáticas, CSIC-UAM-UC3M-UCM, C/ Nicolás Cabrera 15,
Campus Cantoblanco UAM, 28049, Madrid, Spain

Received 30 June 2016 / Received in final form 29 July 2016

Published online 2 May 2017

Abstract. We carry out simulations of the melting process on the phase

change material n-octadecane in squared geometries in the presence of
natural convection and including thermocapillary effects. We show how
the introduction of thermocapillary effects enhances the heat transfer
rate, being the effect especially relevant for small Bond numbers. Thus
induction of Marangoni flows results in a useful mechanism to enhance
the typical slow heat transfer rate of paraffin waxes in applications of
energy storage or passive control management.

1 Introduction

The important environment problems of traditional sources of energy have promoted

an intense and increasing interest in energy efficiency and development of renewable
and clean sources of energy. It has realized that to achieve this goal producing energy
must be accompanied with systems to store it, like thermal energy storage systems.
The usage of Phase Change Materials (PCM) is a low cost, environmentally
friendly and safe solution for thermal energy storage. PCM use the large latent heat
of the solid/liquid phase transition to store a large amount of energy during melt-
ing or release it during solidification, changing very little the temperature during
the transition. Currently, there are hundredths of PCM known -as organic, inorganic
and eutectics- with a wide range of melting points (from temperatures below zero to
hundredths of Celsius degrees).
Among PCM the paraffin waxes are very well suited for thermal storage at room
temperatures, and they are widely used for improving energy efficiency in buildings [1].
Interesting, they are used as well within space crafts in microgravity for the conser-
vation of samples, food and electronic thermal energy management [2, 3]. However, a
major issue in thermal regulation with most PCM, including paraffins, is their low
conductivity, usually < 1 W m−1 K−1 . This leads to very long times during the heat
storage and discharge phases, reducing their applicability.
a
e-mail: santiago.madruga@upm.es
1170 The European Physical Journal Special Topics

On ground applications, the main choice to reduce the problem of low conductiv-
ity is to promote convective motions within the liquid phase of the PCM. Convective
motions driven by gradients of density induced by differences of temperature enhance
the heat transfer rate about an order of magnitude with respect to conductive heat
transport. Another approach to accelerate the heat transfer is to place a large area
of PCM in contact with high conductivity materials, like metallic fins or metallic
foams [4]. Whereas this solution is applicable under microgravity conditions, it in-
creases the mass and size of the devices, and the absence of convective driving may
not be fully compensated.
Prompted by the above considerations a mechanism to enhance the heat transfer
on PCMs without increasing mass and volume, is to maximize the Marangoni flow
induced by thermal gradients of surface tension. We aim at studying the melting
dynamics of the paraffin n-octadecane, which exhibits a solid/liquid phase change at
26.1 ◦ C and is widely used due to its stability and suffers no undercooling.
In Section 2 we present the governing equations of the PCM model and the
geometries selected to study the melting process under natural convection and ther-
mocapillarity. The results of our simulations for natural convection and superposed
thermocapillary driving are presented in Section 3, where we compare results for
two Bond numbers and imposed temperatures. Finally, conclusions are provided in
Section 4.

2 Governing equations and geometry

We have used two squared domain with sizes h = 1 and 2 cm. The left side of the
squares is conductive held at a constant temperature Th , greater than the melting
temperature of the n-octadecane (Tl = 26.1 ◦ C), which is held initially at a solid phase
(Ti = 25 ◦ C). The rest of the sides of the squared geometry are adiabatic. The top
side is considered to be a free surface where the thermocapillary effects are acting.
This configuration corresponds to a classical lateral heating.

2.1 Momentum equation

We consider the flow laminar, two-dimensional and incompressible. The governing

equation expressing the balance of momentum has the vectorial form
 
∂u C(1 − fl )2
ρ + (u · ∇) u = −∇p + μ∇2 u − ρ g [1−α(T −Tref )] ey + u, (1)
∂t δ + fl3

where fl is the local volume liquid fraction of the liquid phase of the PCM in a repre-
sentative control volume. ∇ = (∂x , ∂y ) and ∂t are the spatial and temporal operators.
T the temperature of a control volume that can contain pure solid (fl = 0), liquid
(fl = 1) or a mixture of both phases (0 < fl < 1); u = (u, v) is the fluid velocity; ρ is
the density; μ is the dynamic viscosity, p is the pressure, g the magnitude of gravity
acceleration, ey an unit vector pointing in the vertical direction upwards, Tref is a
reference temperature where physical properties are given and α the thermal expan-
sion coefficient. The bulk physical properties are supposed to be constant within the
range of temperatures studied, with the exception of the density in the buoyancy
term. We have neglected the viscous dissipation in the momentum equation since the
dissipation number of natural convection satisfies (gα/cl )h  1 [9] for the properties
of n-octadecane (c.f. Tab. 1) and scale h used in this work.
 IMA8 – Interfacial Fluid Dynamics and Processes 1171

Table 1. Thermophysical properties of n-octadecane. Values for the solid and liquid states
are listed when distinguished in the equations of the PCM model of this work.
ρ (Kg m−3 ) [5] 776
μ (N s m−2 ) [5] 3.6 × 10−3
cs |cl (J Kg−1 K−1 ) [6] 1934|2196
λs |λl (W m−1 K−1 ) [5] 0.358|0.13
γ (Nm−1 K−1 ) [7] 8.4 × 10−5
Ts |Tl (K) [5] 298.25|299.65
L (J Kg−1 ) [5] 243.5 × 103
α (K−1 ) [8] 9.1 × 10−4

The last term in the momentum equation uses the Carman-Kozeny equation to
model the solid/liquid interface as a porous mushy layer, whose porosity is given by
the liquid fraction fl . In this term δ  1 is a tiny constant to avoid division by zero
without physical meaning, and C is a constant for the mushy region that depends on
the PCM. We set the Darcy coefficient C = 1.6 × 106 kg m−3 s, in compliance with
previous works ([10]). When a control volume is completely liquid (fl = 1) this term is
null, like in a single phase fluid, when it is completely solid (fl = 0) the term diverges
and the velocity of the liquid becomes null, like in a solid. For intermediate values of
fl the PCM is within the mushy zone. This formulation allows to use the momentum
equation in the whole domain without the complication of tracking the solid/liquid
interface ([11]).

2.2 Energy equation

The thermal energy of the system comes from the contribution of the usual sensible
heat, due to changes in temperature in the solid and liquid phases of the PCM, and
from the latent heat content. Assuming the same density in each phase it can be
expressed as a function of temperature as follows
 
∂
+ u · ∇ ρ ((1 − fl )Cpcm,s + fl Cpcm,l ) T
∂t
∂fl
= ∇ ((1 − fl )λpcm,s + fl λpcm,l ) ∇T −ρ L (2)
∂t

where Cpcm,s (Cpcm,l ) and λpcm,s (λpcm,l ) are the specific heats and conductivities of
the PCM in the solid (liquid) phase, averaged with the liquid fraction, and L is the
latent heat of the solid/liquid phase change of the PCM.
The latent heat released by a control volume during the solid to liquid phase
change depends on the melted PCM given by liquid fraction as fl ρ L. As a conse-
quence the coupling between the energy and momentum equation is given through
the liquid fraction field fl , which in turn depends on the temperature, the master
variable of the phase change process. We model the liquid fraction in the mushy zone
using a linear relationship between the solidus and liquidus temperatures
⎧
⎨ 0, T ≤ Ts
⎪
fl = 1, T ≥ Tl (3)
⎪
⎩
(T − Ts )/(Tl − Ts ), Ts < T < Tl .
1172 The European Physical Journal Special Topics

2.3 Boundary conditions

The dependence of interfacial tension with the temperature at the free surface on top
is approximated by the linear state equation

σ = σ0 − γ(T − Tref ) (4)

where σ0 is the interfacial tension at the reference temperature Tref , and γ =

− ∂σ/∂T |Tref accounts for the dependence of interfacial tension with temperature.
The lateral heating along the free surface generates at Marangoni flow from the
balance between shear force and surface tension

μ ∂z u = −γ∂x T at free surface. (5)

There is no penetration of the liquid on the free surface and the free surface is
supposed to be adiabatic.

w = 0 at free surface (6)

∂T
= 0 at free surface. (7)
∂z
We use the open source software OpenFOAM, based on finite volumes, to sim-
ulate the equations of the PCM model discussed above together the boundary con-
ditions. We have validated the implementation of the Marangoni effect with results
by Bergmann and Keller [12], who simulated the combined effect of buoyancy and
thermocapillarity. We find a grid of 400 × 400 cells is enough to guarantee convergence
of solutions.

3 Results
We have carried out simulations for both geometries with and without thermocap-
illary effects to evaluate their impact on the melting dynamics. We can estimate a
Bond number Bo = Ra/M a = αρ g h2 /γ, where Ra and M a are the Rayleigh and
Marangoni numbers respectively, assuming the domain filled with liquid PCM: ∼ 8.3
for the small square and ∼ 33 for the large one. Thus a major contribution of ther-
mocapillarity is expected on the smaller domain.
Figure 1 shows the evolution of the volume fraction of melted PCM (global liq-
uid fraction) for both geometries. In presence of natural convection, melting time
at Th = 50 ◦ C for 1 cm square is 1378.1 s, and it is reduced to 1130 s (18%) when
Marangoni driving is superposed. At higher temperature Th = 100 ◦ C melting time
for natural convection is 426 s, and is reduced to 332 s (28.3%) when Marangoni
driving is superposed. The melting time for Th = 100 ◦ C exhibits an improved rela-
tive difference of 10.3% compared to Th = 50 ◦ C when thermocapillarity is included.
This shows that is not only the ratio between natural convection and thermocapillary
convection, expressed by the Bond number, that dictates the melting rate: it is as well
the average velocity of the liquid at the free surface that contributes to the magnitude
of the advance of the solid/liquid interface. Thus, the speed of the melting front is
affected by the rate of latent heat loss at the interface, the strength of convective
motions driven by density gradients and driven by surface tension gradients.
The lower conductivity of the liquid phase with respect to the solid phase (c.f.
Tab. 1) leads to a quick advance of the overall liquid fraction in the first stages of
melting due to enhanced conductive transport within the solid phase of PCM. How-
ever, once the liquid phase of the PCM becomes more abundant the melting slows
 IMA8 – Interfacial Fluid Dynamics and Processes 1173

Fig. 1. Evolution of the global liquid fraction of n-octadecane for a square of size 1 cm (left)
and 2 cm (right). Dashed lines correspond to Th = 323 K, solid to Th = 373 K. NC at the
legends denotes Natural Convection.

down. Thus, for instance, it takes 59% of total melting time to liquefy the remain-
ing 50% of solid PCM of 1 cm square at Th = 50 ◦ C when only natural convection
acts. The melting of the first 50% is even faster when thermocapillary effects are
included, taking 69% of the total melting time to liquefy the remaining 50% of PCM
volume. This difference on heat transfer rate performance is enhanced by the fact
that convective motions due to thermocapillarity appear quicker than convective mo-
tions due to natural convection, which requires about 1.6 mm of melted n-octadecane
to destabilize the conductive state by buoyancy. Indeed for 2 cm square the relative
strength of Marangoni effects is weaker, and the time difference between buoyancy
and superposed to thermocapillarity cases to melt the first 50% is reduced to 4%. At
Th = 100 ◦ C similar relative differences are found.
Figure 2 exhibits snapshots of streamlines and temperature fields across represen-
tative times of the melting process of n-octadecane for 1 cm square at Th = 50 ◦ C and
Th = 100 ◦ C. Times have been chosen to roughly exhibit 0.1, 0.5 and 0.8 of the global
liquid fraction. For natural convection driving, the liquefied PCM flows from the hot
left side and later cool down along the solid/liquid interface, returning and creating
a single large convective cell. This global structure is conserved until complete melt-
ing of the PCM. Ascending hot melted PCM creates a faster advance of the interface
near the hot adiabatic free surface, combined with a slow motion at the cooler bottom
leads to an inclined interface which as time advances confines progressively the solid
PCM on the right bottom corner. The detachment of the solid/liquid interface from
the free surface is produced at both temperatures ∼ 58% of the complete melting
time. This melting dynamic is similar for Th = 50 ◦ C and Th = 100 ◦ C.
As thermocapillary driving is superposed to natural convection (c.f. third and
fourth rows of left block of Fig. 2) not only enhanced heat transfer rates result, as
explained before, but also a qualitative change in the melting dynamics. First, there
is a quick advance of the solid/liquid interface at the free surface that detaches from
it at 130 s (12% of complete melting time) for Th = 50 ◦ C, and 25 s (7.5% of complete
melting time) for Th = 100 ◦ C, liquefying the PCM at the upper half of the square.
The reduction of time to detachment, from 12% to 7.5% is a consequence of the
stronger thermocapillary flow with increasing Th . Second, the convective patterns
become more involved, with the creation of a thin convective cell below the free
surface near of the left hot side at the early stages of melting, and the creation of two
large counter-rotating cells at intermediate times. The smaller convective cell at the
lower part carries cooler liquid from the top upwards along the solid/liquid interface.
Figure 3 exhibits snapshots for the temperature field and streamlines at represen-
tative times of 2 cm square, with Bo = 33. Convective motions are more involved with
respect to 1 cm square due to greater buoyancy effects. Thus streamlines become less
1174 The European Physical Journal Special Topics

Fig. 2. Square of side 1 cm. Snapshots at representative times of the temperature field
(second and fourth rows) and streamlines (first and third rows) for Th = 50 ◦ C (left block)
and Th = 100 ◦ C (right block). Melting evolution is shown when only natural convection
is acting (first and second rows), and when natural convection plus thermocapillarity are
acting together (third and fourth columns).

Fig. 3. Square of side 2 cm. Snapshots at representative times of the temperature field
(second and fourth rows) and streamlines (first and third rows) for Th = 50 ◦ C (left block)
and Th = 100 ◦ C (right block). Melting evolution is shown when only natural convection
is acting (first and second rows), and when natural convection plus thermocapillarity are
acting together (third and fourth columns).
 IMA8 – Interfacial Fluid Dynamics and Processes 1175

smooth at all times, and at the latter stages of evolution new convective cells appear
with respect to 1 cm square. However, overall dynamics of the solid/liquid front is
similar as the discussed above for smaller Bond number. The weakening of the im-
pact of thermocapillary flows at this large Bo on performance of heat transfer rate is
appreciated comparing the effect of increasing Th with respect to the previous case.
Now the introduction of thermocapillary flows decreases the melting time to just 4.2%
for Th = 50 ◦ C, and 6.3% for Th = 100 ◦ C, and the relative difference between both
temperatures is 2.1%, compared with 10.3% for Bo = 8.3. In addition, the time to
detachment of the solid/liquid interface from the free surface when thermocapillarity
is included is 11% of the total melting time for Th = 50 ◦ C, and 8% for Th = 100 ◦ C.
Thus the difference is reduced to 3%, compared to 4.5% for Bo = 8.3.

4 Conclusions
We have simulated the melting dynamics of the phase change material n-octadecane
using an enthalpy-porosity formulation for the Navier-Stokes coupled to energy equa-
tions. Two squared geometries of 1 cm and 2 cm subjected to lateral heating, with
respective Bond numbers 8.3 and 33 have been used to evaluate the influence of
thermocapillary effects on heat transfer rate performance. We have calculated the
global liquid fraction curves when convective motions are driven only by buoyancy
and when thermocapillary flows are included. We find for Bo = 8.3 that adding ther-
mocapillary effects decreases the total melting time 18% with Th = 50 ◦ C and 28%
with Th = 100 ◦ C. At Bo = 33 the influence of thermocapillary effects weakens with a
decrease of the total melting time of 4.2% with Th = 50 ◦ C and 6.3% with Th = 100 ◦ C.
Snapshots of stream functions and temperature fields have allowed following the
evolution of melting dynamics. We find that overall evolution of solid/liquid inter-
face is similar for both geometries and imposed external temperatures for buoyancy
driving convection. The same happens when thermocapillary driving is superposed.
However there are relevant qualitative and quantitative differences when buoyancy
only, and added thermocapillarity are compared. The last case produces a very quick
detachment of the solid/liquid front from the upper free surface.
From our results, we conclude that thermocapillary effects at small Bo are a
very relevant mechanism to enhance the heat transfer rate on n-octadecane, espe-
cially when subjected to high external temperatures. Which happens, for instance,
in applications of thermal control management on electronic devices. Importantly,
this improvement of the heat transfer rate performance can be accomplished without
adding complexity like embedding high conductivity metallic structures within the
PCM, how is usually done.

Santiago Madruga acknowledges support by Erasmus Mundus EASED programme (Grant

2012-5538/004-001) coordinated by Centrale Supélec, and the Spanish Ministerio de
Economı́a y Competitividad under Projects No. TRA2013-45808-R, No. ESP2013-45432-P
and No. ESP2015-70458-P. Carolina Mendoza acknowledges support by the Ministerio de
Economı́a y Competitividad under Project No. MTM2014-56392-R. We thank Gonzalo
Sáez-Mischlich for code development and Pedro Mongelos for support with postprocessing.

References
1. S.A. Memon, Phase change materials integrated in building walls: A state of the art
review, Renew. Sustain. Energy Rev. 31, 870 (2014)
2. X. Gui, D. Tang, S. Liang, B. Lin, X. Yuan, Influence of void ratio on thermal perfor-
mance of heat pipe receiver, Int. J. Heat Fluid Flow 33, 109 (2012)
1176 The European Physical Journal Special Topics

3. T.Y. Kim, B.S. Hyun, J.J. Lee, J. Rhee, Numerical study of the spacecraft thermal
control hardware combining solid-liquid phase change material and a heat pipe, Aerosp.
Sci. Technol. 27, 10 (2013)
4. D. Fernandes, F. Pitié, G. Cáceres, J. Baeyens, Thermal energy storage: How previous
findings determine current research priorities, Energy 39, 246 (2012)
5. C.J. Ho, J.Y. Gao, Preparation and thermophysical properties of nanoparticle-in-paraffin
emulsion as phase change material, Int. Commun. Heat Mass Transf. 36, 467 (2009)
6. E.M. Alawadhi, Thermal analysis of a building brick containing phase change material,
Energy Build. 40, 351 (2008)
7. D.R. Lide, ed., CRC Handbook of Chemistry and Physics, 84th edn. (CRC Press, 2003)
8. N.S. Dhaidan, J.M. Khodadadi, T. Al-Hattab, S.M. Al-Mashat, Experimental and
numerical investigation of melting of phase change material/nanoparticle suspensions
in a square container subjected to a constant heat flux, Int. J. Heat Mass Transf. 66,
672 (2013)
9. B. Gebhart Effects of viscous dissipation in natural convection, J. Fluid Mech. 14, 225
(1962)
10. S.S. Sebti, M. Mastiani, H. Mirzaei, A. Dadvand, S. Kashani, S.A. Hosseini, Numerical
study of the melting of nano-enhanced phase change material in a square cavity,
J. Zhejiang Univ. Sci. A 14, 307 (2013)
11. C. Beckermann, R. Viskanta, Natural convection solid/liquid phase change in porous
media, Int. J. Heat Mass Transf. 31, 35 (1988)
12. T.L. Bergman, J.R. Keller. Combined buoyancy, surface tension flow in liquid metals.
Numerical Heat Transfer 13, 49 (1988)