Applied Thermal Engineering 181 (2020) 115990

Contents lists available at ScienceDirect

Applied Thermal Engineering

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

Buoyancy-driven melting and heat transfer around a horizontal cylinder in

square enclosure filled with phase change material
Azim Memon a, Garima Mishra a, b, Anoop K. Gupta a, *
a
Department of Chemical and Biochemical Engineering, Indian Institute of Technology Patna, Patna 801106, Bihar, India
b
Department of Chemical Engineering, Indian Institute of Technology Kanpur, Kanpur 208016, U.P., India

A R T I C L E I N F O A B S T R A C T

Keywords: This work primarily focuses on designing a cost effective method yielding high rate of heat transfer from a heated
PCM object to ensure the faster melting of the phase change material (PCM). A numerical investigation has been
Natural convection performed on the laminar natural convection from a two-dimensional heated circular cylinder confined in a
Melting
square enclosure filled with a phase change material, namely, the lauric acid. In particular, the coupled mo­
Circular cylinder
Square enclosure
mentum and energy equations are solved to delineate the influence of the geometric position of the cylinder
within the square enclosure on the melt and heat transfer characteristics of the phase change material. The solid-
liquid phase change is formulated using the enthalpy–porosity technique. The detailed velocity and temperature
fields are visualized in terms of velocity vectors, isotherm contours, melt fraction contours, melting rate, energy
storage and surface averaged Nusselt number. The results reported herein are restricted to the six distinct
geometrical positions (i.e., four symmetric and two asymmetric) of the cylinder within the enclosure. The rate of
convective heat transport and hence the melt fraction is found to decrease with the increasing distance of the
cylinder from the base of the enclosure. In addition, the influence of the mushy zone parameter (Amush) on the
melting behaviour has also been explored. The reported numerical results for the melt fraction have been
correlated as a function of the melting time and the geometric positions of the heated cylinder. Finally, one can
conclude that it is possible to realize an enhanced rate of melting and energy storage simply by adjusting the
location of the cylinder under confinement and thus the performance of a thermal energy storage (TES) system
can remarkably be improved.

resources into the energy system [2,3].

Amongst the other thermal energy storage methods, the latent heat
1. Introduction
thermal energy storage (LHTES) systems are one of the most promising
ways of storing thermal energy due to their high energy storage density
The unrelenting decrease in the availability of the conventional en­
(5–14 times more energy per unit volume compared to sensible thermal
ergy resources such as fossil fuels, and the continuous increase in the
energy devices), and quasi-isothermal melting/solidification process
demand and supply chain gap have been a driving force for the utili­
[4]. Phase change materials (PCMs) are latent heat storage substances
zation of renewable energy resources. Thus, for many years, the capa­
that absorb and release large amount of heat at nearly constant tem­
bilities and capacity of these non-conventional form of energy resources
perature corresponding to the phase transition temperature of the ma­
in saving the world from energy crisis have been explored at large. As of
terial. The development and use of PCMs for efficient energy utilization
2016, they contribute to the 24% of electricity generation globally [1].
is currently the most active research field. Several comprehensive re­
However, the intermittent and fluctuating supply of energy (mainly
views exist in the literature critically examining the potential of PCMs in
solar and wind energy) has been a great challenge for the conventional
hybrid/integrated energy storage devices [2–10]. The large number of
application of energy derived from these sources. Thermal energy stor­
PCMs are now available over wide ranges of operating temperature
age (TES) systems have been a key technology to overcome the vari­
which makes them a commercially viable option for a number of varied
ability and uncertainty associated with the non-conventional energy
applications such as heat and cold storage devices [6,7,9], thermal
resources [2]. They offer a stable, efficient and environment friendly
processing of food items and temperature sensitive materials [11],
means of waste heat recovery for integrating the renewable energy

* Corresponding author.
E-mail addresses: azim_2021cb03@iitp.ac.in (A. Memon), garima.mishra.12@gmail.com (G. Mishra), anoopg@iitp.ac.in (A.K. Gupta).

https://doi.org/10.1016/j.applthermaleng.2020.115990
Received 16 May 2020; Received in revised form 29 July 2020; Accepted 29 August 2020
Available online 7 September 2020
1359-4311/© 2020 Elsevier Ltd. All rights reserved.
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

1
Nomenclature Vx, Vy x-, y-components of velocity, ms−
x, y Cartesian coordinates, m
Amush mushy zone parameter, kg.m− 3s− 1
Cp specific heat capacity, J.kg− 1K− 1 Greek symbols
d diameter of the cylinder, m β coefficient of volume expansion, K− 1
g acceleration due to gravity, ms− 2 λ melt fraction, dimensionless
H specific enthalpy, J.kg− 1 δ distance between centre of the enclosure and centre of the
k thermal conductivity, W.m− 1K− 1 cylinder, m
L side length of the square enclosure, m μ viscosity of PCM, Pa.s
Lf latent heat of fusion of the PCM, J.kg− 1 ρ density of fluid, kg.m− 3
Nu surface-average Nusselt number, dimensionless Subscripts
p pressure, Pa x x-axis
Pr Prandtl number, dimensionless y y-axis
q′′ heat flux from the surface of the cylinder, W.m− 2 ref reference conditions
Ra Rayleigh number, dimensionless s solid phase
Ste* modified Stefan number, dimensionless l liquid phase
t time, s
TC initial temperature of the PCM in the enclosure, K Abbreviation
TH temperature at the surface of the cylinder, K PCM phase change material
Tm melting temperature of the PCM, K TES thermal energy storage
V velocity vector, ms− 1

cooling of electronic components [12,13], solar heating systems [33] and Ismail and Stuginsky [34] were among the early researchers to
[14,15], thermal storage in buildings [16,17], cooling helmets [18], determine the influence of natural convection on the melting rate and
textiles [19], HVAC units [20], biomedical applications [21], electric also to delineate the comparison among the various solution strategies
vehicles [22], greenhouses [23], etc. Some commonly used and employed. The solution methods for the numerical modelling of the
commercially available PCMs are paraffin wax (Tm ≈ 57 ◦ C–65 ◦ C), melting/solidification process have been divided in two categories,
lauric acid (Tm ≈ 46 ◦ C), liquid gallium (Tm ≈ 29.8 ◦ C), RT-25 (Tm ≈ namely, the fixed grid methods and the adaptive meshing methods. The
26 ◦ C), capric acid (Tm ≈ 32 ◦ C), n-octadecane (Tm ≈ 28 ◦ C), dodecanoic adaptive meshing methods or the moving mesh methods offer an
acid (Tm ≈ 43 ◦ C), and many salt based PCMs such as CaCl2⋅6H2O (Tm ≈ advantage of good precision of the numerical results. However, the
27 ◦ C–29 ◦ C), etc. complex formulation limits their application to simple geometries
However, for the effective utilization of the PCMs for a specific [35,36]. The fixed grid methods employ the enthalpy formulation for
application, the physical and themo-chemical properties, availability the energy equation. These methods can efficiently model the melting
and economic viability of the PCM are some of the key selection pa­ process by considering either an isothermal melting process or by
rameters. A PCM should have a suitable phase change temperature as simplifying the phase change process by creating a mushy zone (semi-
per the application, high latent heat, large thermal conductivity both in solid, two-phase zone) in between the two-phases. The isothermal
the liquid and solid phases, negligible sub-cooling, low vapour pressure modelling strategy gives accurate solutions but has a prerequisite for the
for small change in volume during phase transition, chemical and interface boundary to solve the energy equation. On the other hand, the
physical stability and non-toxicity, etc. [2,7]. mushy zone approach reduces the computational efforts with similar
Besides, detailed research efforts need to be directed to examine the accuracy levels and can feasibly be used to model complex multi-
melt characteristics for efficiently designing a latent heat thermal energy dimensional problems with ease [36–40]. The flow through the
storage unit. The analysis of heat and mass transfer in phase transition mushy-zone is modelled as flow through a porous medium using a
problems is defined as a separate class of moving boundary problems. modified form of Carman-Kozeny equation. This avoids numerical in­
During the process of melting or solidification, the solid-liquid interface stabilities in reaching sharp interface boundaries and considers a smooth
changes as the heat is absorbed or released into the system. The differ­ transition from a solid to a liquid state or vice versa. However, the grid
ence in the thermo-physical properties of the two phases and the diffi­ resolution in resolving the gradients in the mushy zone and the selection
culty in tracking the moving interface complicate the study of melting of mushy zone parameter (Amush) is crucial to the efficient modelling of
and solidification process [8]. This class of problems were first examined the melting phenomenon [41,42]. The influence of mushy zone
by Stefan in 1889 for ice formation [24]. Henceforth, several solutions parameter has been numerically examined in detail by Fadl and Eames
(analytical and numerical) were made available for the Stefan’s problem [42] for melting inside the vertically and horizontally oriented thermal
using different techniques. The various analytical and numerical energy system.
methods used to solve this class of problems have been extensively The melting of PCM within the enclosures of various shapes such as
reviewed in the literature [4,8,9]. Very few analytical solutions are ducts, pipes, annulus, spherical shells, etc., has been extensively exam­
available in the closed form, and that too only for idealized one- ined over the years. Excellent reviews on the numerical modelling of the
dimensional geometries with simplified boundary conditions [25,26]. PCM melting inside in geometries of various shapes exist in literature
Numerical techniques based on finite element and finite volume [9,43]. Amidst all, examining the melting of PCM in rectangular en­
methods have been reported in solving such moving boundary problems closures forms a fundamental configuration of great relevance and in­
with improved accuracy [27–30]. Initially, the problem was approached dustrial importance [44–46]. Kamkari et al. [46] experimentally studied
by considering conduction to be the sole mode of heat transfer. Such a the effect of inclination angle on the melt characteristics for the PCM
simplification does not correctly capture the intricacies in a complex (lauric acid) in a rectangular enclosure. They found that the melting rate
process such as phase change. Hence, later research in the field was is high for a horizontally placed enclosure compared to that for a vertical
focused to incorporate and examine the influence of natural convection enclosure. Moreover, recently Elsayed [47] has investigated the effect of
on the melting process. Gobin [31], Hahn and Ozisik [32], Eckert et al. hot air stream on the PCM (paraffin wax) melting inside the triangular

2
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

Adiabatic wall

Top asymmetric PCM at TC

L = 4 cm
(δx = –L/4, δy = +L/4)
Top
(δy = +L/4)

Adiabatic wall
y
Adiabatic wall
Centre
TH (δ = 0)
Left x
(δx = –L/4) d
Bottom
(δy = –L/4)
Bottom asymmetric
(δx , δy = –L/4) L = 4 cm

Adiabatic wall
Fig. 1a. Schematics of the computational domain and grid structure.

containers of constant volume but different apex angles (24◦ , 60◦ and pure Gallium as the phase change material using finite-volume method.
100◦ ) using ANSYS-Fluent. They observed the higher heat storage ca­ They examined the configuration for the two specific boundary condi­
pacity in the acute triangular cylinders in comparison to the obtuse and tions of constant wall temperature and constant heat flux on the surface
equilateral cylinders. of the cylinder. They reported that the large part of energy is transmitted
Next, free convection heat transfer from a cylinder in a solid enclo­ by means of natural convection from the upper part rather than the
sure forms the prototype case for a broad range of melting/solidification lower part of the heated cylinder. Besides, several strategies have been
applications utilizing PCMs such as in heat exchangers, mixing vessels, employed to improve the heat transfer and melting rate in a phase
solar cells, recyclable batteries, power plants, etc. Here, the heat transfer change process by increasing the thermal conductivity of PCMs [58,59],
rate and so the melting is limited by the interaction between the cylinder by using extended surfaces such as finned cylinders/walls [60] or by
surface and the enclosure walls, and the thermal nature of confining micro-encapsulation of the PCM [61]. Another strategy can be due to
walls (isothermal, adiabatic, etc.) in addition to the properties of phase simple modifications to the geometrical configuration which leads to an
change materials. White et al. [48] and Yao and Chen [49] were amongst alternative but cost efficient method to achieve the higher charging/
the first researchers to examine the problem of melting of the PCM melting rate in the PCM. Notwithstanding the existence of literature
around a circular cylinder using perturbation analysis. However, these dealing the melting/solidification of the PCM around a horizontal cyl­
studies have limitations of applicability to only initial melting phase inder, to the best of authors’ knowledge, there is no comparative study
where conduction is the dominant mode of heat transfer and buoyancy- available focusing on delineating the influence of various symmetric and
induced currents are weak. Bathelt and Viskanta [50] initiated the asymmetric locations of the cylinder within the square enclosure filled
experimental investigation for melting of paraffins around a heated with a PCM. Thus, the present work aims to fill this gap and examines
circular cylinder. They recorded the melt fraction picture at different the effects of various possible geometrical adjustments in the location of
time instants of melting and reported the heat transfer coefficient at the circular cylinder placed inside a square enclosure on the melting
solid-liquid interface. Cheng et al. [51] conducted experiments to study process of the lauric acid as the PCM. The work has been carried out for
the effect of natural convection on ice formation around a cold circular the six different geometric placements/locations of the cylinder within
cylinder. They reported the values of heat transfer coefficients and the the enclosure, namely, centre, left, top, bottom, top asymmetric and
average Nusselt number behaviour at the ice-water interface. bottom asymmetric. The influence of the mushy zone parameter (Amush)
Later, Rieger et al. [52] used the body fitted coordinates to numer­ on the melting rate, melting interface and on the absorbed energy has
ically compute this problem and emphasized on the importance of the also been explored. Apart from this, the work has also been extended for
convective currents in the melting process. Prusa and Yao [53,54] seven other shapes of the cylinder to develop an insight into the melting
studied the enhancement in the melting rate for a heated horizontal behaviour around variously shaped 2D heated objects and thus to
cylinder by means of the numerical and perturbation methods for the enhance the applicability of the work in a broader spectrum.
conditions of constant heat flux and constant wall temperature pre­
scribed on the cylinder surface. The solutions were valid for short times 2. Problem statement and mathematical formulation
after the initiation of melting process. Sasaguchi et al. [55] also used the
body-fitted coordinates to examine the melting around a single and two Consider an isothermal (at temperature, TH) circular cylinder (long
horizontally and vertically placed circular cylinders in an enclosure. in z-direction) of diameter ‘d’ placed horizontally in a square enclosure
They examined the influence of the initial temperature on the melting of side length L = 4 cm. The confinement ratio (d/L) is taken to be 0.1.
rate and compared the rate of melting for the two cases of a single and The enclosure is filled with a solid PCM, namely, the lauric acid kept at a
double cylinders inside the enclosure. temperature TC where TC < TH. The effect of geometric placement of the
Additionally, Sugawara et al. [56] numerically examined the impact cylinder within the square enclosure has been analysed for both sym­
of the vertical position of a circular cylinder on freezing/melting of metric (centre, left, top, bottom) and asymmetric (top asymmetric, bottom
water/ice around the cylinder in a square cavity. They concluded that symmetric) positioning of the cylinder with respect to both the horizontal
melting significantly decreases with increasing cylinder distance from and vertical centreline of the enclosure, as shown in Fig. 1a. In this work,
the base of the enclosure. Mahdaoui et al. [57] numerically studied the δ = 0 refers to the centre position, δx = − L/4 to left position, δy = +L/4
melting process around a horizontal cylinder in an enclosure filled with to top position, δy = − L/4 to the bottom position. Two asymmetric

3
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

positions of the cylinder considered in this study, namely, top asym­ ∫ T

metric and bottom asymmetric, can be referred as (δx = − L/4, δy = +L/ h = href + Cp dT (5)
4) and (δx = − L/4, δy = − L/4), respectively. Owing to the prevailing Tref

temperature difference, the solid lauric acid acquires the heat from the The latent enthalpy used in Eq. (4) is related to the liquid or melt
hot circular cylinder and melts as the temperature increases beyond the fraction (λ) during the phase change as follows:
average melting temperature, Tm. Consequently, the rise in temperature
results in a decrease in the density of the fluid melt thereby giving rise to ΔH = λ(T)Lf (6)
buoyancy-induced flow currents in the upward direction; thus facili­
tating heat transfer by natural convection. During the melting process, where, href is the enthalpy at reference temperature Tref , Lf is the latent
the phase change material transits from solid to liquid phase via a heat and λ(T) is the melt fraction defined as follows:
transition zone known as the “mushy zone”. Mushy zone can be ⎧
0 if T < (Tm − ΔT)
⎪
considered as a semi-solid zone between the purely solid and liquid ⎪
⎪
⎨
T − Tm + ΔT
phases. The transition temperature range (or melting range) considered λ(T) = if (Tm − ΔT)⩽T⩽(Tm + ΔT) (7)
⎪ 2ΔT
in the present work is nearly 4.7 K (=2ΔT) where (Tm–ΔT) is considered ⎪
⎪
⎩
to be the “solidus”, and (Tm + ΔT) to be the “liquidus” temperature. For 1 if T > (Tm + ΔT)
the solid-liquid transition region inside the PCM, the Enthalpy-Porosity The melt fraction λ(T) varies between 0 and 1. It is zero when the
approach [39,62,63] has been employed. This technique does not track PCM is solid phase whereas it is 1 for the liquid phase. The energy source
the solid-liquid interface explicitly. Instead, it treats the mushy region as term, Sh, used in Eq. (3) is given as:
porous with the porosity equal to the melt fraction (or volume fraction)
associated with each cell in the computational domain. The melt fraction Sh =
∂(ρ ΔH)
+ ∇⋅(ρ VΔH) (8)
indicates the fraction of the cell volume that is in liquid form. In order to ∂t
keep the level of the complexity to a tractable level, the following The momentum source term, Sm, in Eq. (2) describes the flow in a
pertinent assumptions have been made: porous medium and can be defined as a function of melt fraction as
follows,
• The melted PCM behaves as a Newtonian fluid.
• The flow is two-dimensional, laminar and incompressible (except for Sm = − A(T) V (9)
the body force term). The function A(T) is defined using a modified form of Carman-
• The properties of PCM are different for solid and liquid phases. Kozeny equation [39,63] as follows:
• Both the solid and liquid phases are homogeneous, isotropic and are
in thermal equilibrium with the interface. Amush (1 − λ(T))2
A(T) = (10)
• Thermal radiation and viscous dissipation effects are negligible. λ3 (T) + ε
• The density stratification within the fluid elements are well governed
In the above expression, the function A(T) accounts for the velocity
by the Boussinesq approximation given as, ρl − ρ = ρl β(T − Tm ),
switch-off in solid PCM. When the entire PCM melts (for T > Tm + ΔT),
where β is the coefficient of thermal expansion of PCM.
the melt fraction λ(T)⟶ 1 and the source term (Sm) in the momentum
equation becomes zero approaching the general momentum conserva­
The thermo-physical properties of the PCM (density, ρpcm; thermal
tion equation applicable for a Newtonian fluid flow. However, when the
conductivity, kpcm; and specific heat capacity, Cp,pcm) except the viscosity
temperature of PCM is lower than the solidus temperature (i.e., T < Tm –
are the functions of temperature due to their dependence on the melt
ΔT), the entire PCM is solid and hence λ(T)⟶ 0. In this case, the flow
fraction.
velocity is nearly zero thereby leading to minimal contributions from the
source term. Thus, the source term shows its prominent contributions
during the melting process in the mushy zone only (i.e., within melting
2.1. Governing equations regime). In this regime, the mushy zone parameter Amush plays a major
role in controlling the heat transfer and the melting characteristics of
Within the framework of above-mentioned assumptions, the coupled phase change material. The values of this parameter depends on the
velocity and temperature fields are governed by the time-dependent morphology of the PCM and significantly influence the location of the
continuity, momentum and energy equations, which are written as solid-liquid melt interface and its progression during the melting pro­
follows: cess. Sufficiently large values of Amush suppress the flow in the regions
Continuity equation: where the PCM is solid and thus correctly predict the solid-liquid
∂ρ interface. Low values of Amush lead to unrealistic predictions in delin­
+ ∇.(ρV) = 0 (1) eating the moving solid-liquid interface whereas excessively values of
∂t
Amush cause the results to oscillate. Alternatively, very large values of
Momentum equation:
Amush assume more volume of the PCM in the mushy region which is
∂(ρV) mostly static in nature. This under-estimates the actual behaviour of the
+ ∇.(ρV) = − ∇p + μ∇2 V + ρg + Sm (2)
∂t PCM under natural convection, however, very small values of Amush
over-predict the extent of melting. Hence, a careful selection of Amush
Energy equation:
parameter is crucial in envisaging the numerical modelling of the
∂(ρH) melting process. Values between 104 and 107 kg/m3s are the recom­
+ ∇.(ρVH) = k∇2 T + Sh (3)
∂t mended values for most calculations [42,62]. These aspects have been
In the aforementioned Eqs. (1)–(3), V is the velocity vector, p is critically reviewed and discussed in ref. [42]. One may also find a
pressure, T is temperature, Sm is the momentum source term and H is the detailed analysis and comparison between the choice of various values
specific enthalpy. The specific enthalpy (H) is defined as the sum of of Amush and experimental data in the study of Kheirabadi and Groulx
sensible enthalpy (h) and latent enthalpy (ΔH) as follows: [64] who recommended Amush = 106 kg/m3s to ensure the best overall
match between their numerical predictions and the experimental re­
H = h + ΔH (4) sults. In the present work, Amush = 105 kg/m3s has been employed.
Moreover, the influence of Amush on the melting rate and the progression
where the sensible enthalpy can be expressed as:

4
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

Table 1
• Specific heat capacity
Thermo-physical properties of the lauric acid as the PCM [42].
Property Solid Liquid The specific heat capacity for a PCM can be written as:
Liquidus temperature, Tl = Tm + ΔT (K) – 321.35
Solidus temperature, Ts = Tm–ΔT (K) 316.65 – Cp,pcm (T) = Cp,s + (Cp,l − Cp,s )λ(T) (13)
Melting temperature, Tm (K) 319
The thermo-physical properties of the PCM are calculated at the
Density, ρ (kg/m3) 940 885
Specific heat, Cp (J/kg.K) 2180 2390 liquidus temperature (Tl) for the evaluation of the Rayleigh number,
Thermal conductivity, k (W/m.K) 0.16 0.14 Rad = βg(TH − Tm ) d3 /υl αl , Prandtl number, Pr = Cp,l μ/kl , and modified
Viscosity, µ (kg/m.s) – 0.008 Stefan number, Ste* = {Cp,s (Tm − TC ) + Cp,l (TH − Tm )}/Lf , which are
Thermal expansion coefficient, β (1/K) 0.0008
valued as ~2 × 104, ~140 and ~0.55, respectively. In the end, the
–
Latent heat of fusion, Lf (J/kg) 187,210 –
problem definition is completed by specifying initial condition and
physically realistic boundary conditions for this configuration. For the
of solid-liquid interface have also been examined in the present work for velocity field, the usual no-slip boundary condition (V = 0) is prescribed
a broad range of 104 ≤ Amush ≤ 106 kg/m3s. Therefore, the range of on the enclosure walls and the cylinder surface. For temperature, the
values of Amush considered herein are consistent with that reported enclosure walls are thermally insulated, i.e., ∂T/∂ns = 0 whereas the
elsewhere [42,62,64]. A factor ‘ε’ is introduced into the Eq. (10) to ac­ cylinder surface acts as a heat source and is maintained at a fixed tem­
count for the numerical applicability of the equation in the solid phase, i. perature of TH = 343.15 K (70 ◦ C) for all time t. Initially (t = 0), the
e., at λ = 0. A constant value of ε = 0.001 has been used here. temperature of the PCM filled in the enclosure is maintained at 298.15
K.
2.2. Thermo-physical properties Finally, the governing equations subjected to the above boundary
conditions are solved in terms of the primitive variables (V, p, T) and the
The thermo-physical properties of PCM (lauric acid) are listed in results are post-processed to evaluate the local and the derived quanti­
Table 1. As mentioned earlier, the thermo-physical properties of the ties to discuss the momentum and heat transfer aspects. The spatial and
PCM (except viscosity) are considered to be temperature-dependent and temporal evolution of the velocity and temperature field is studied using
thus are functions of the melt fraction (λ) of the PCM. This functional velocity vectors, temperature profiles, melt fraction contours, melting
dependence has been approximated in the present work as follows: and energy storage rate, and surface average Nusselt number.
For temperatures beyond the melting regime, i.e., T < (Tm–ΔT) and T
• Density > (Tm + ΔT), the energy is gained by the PCM through sensible heating
only whereas within the melting regime, the latent heat also contributes
The density variation of the PCM can be written as follows: to the energy transfer. The total energy stored within the PCM is
calculated by integrating the enthalpy over the control volume.
ρpcm (T) = ρs + (ρl − ρs )λ(T) (11)
The surface average Nusselt number is calculated by integrating the
temperature gradient over the surface of cylinder as follows:
• Thermal conductivity ∫ ( ⃒ )
hd d ∂T ⃒⃒
Nu = = − dAs (14)
kl As (TH − Tm ) As ∂ns ⃒s
Similar to the above expression for density, the thermal conductivity
can be defined by the following expression:
where, h is the surface averaged heat transfer coefficient defined as
kpcm (T) = ks + (kl − ks )λ(T) (12) follows:

Fig. 1b. Grid sensitivity test for centre (δ = 0) and bottom (δy = − L/4) position of the cylinder: (i) melt fraction; (ii) average temperature of the PCM.

5
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

Fig. 2a. Comparison of time variation of melt fraction for the PCM filled inside the enclosure: (i) horizontal enclosure; (ii) vertical enclosure. Solid lines: present
results, broken lines: Fadl and Eames [42], symbols: Kamkari et al. [46].

q′′ (maximum for top and top asymmetric positions).

h= (15)
(TH − Tm ) The computational domain in the present case is fixed by the ge­
ometry of the problem itself. However, it is crucial to conduct a sys­
3. Numerical solution methodology tematic grid analysis for our results to be free from numerical artefacts.
Three non-uniform grid structures, namely, G1 (~20,500 CVs), G2
In present work, numerical simulations have been carried out using (~40,000 CVs) and G3 (~70,200 CVs), varying in terms of the number
the finite volume method (FVM) based commercial CFD code ANSYS of grid points on the cylinder surface, and the smallest grid size in the
Fluent (version 19.2). The geometry and mesh creation was done using computational domain were tested. The relative performance of the
ANSYS Workbench-19.2. The unstructured quadrilateral cells with non- grids was assessed by comparing the time evolution of melt fraction and
uniform spacing were generated in the computational domain. Due to average temperature of the PCM profiles (Fig. 1b). Negligible differences
the expected steep gradients of velocity and temperature, a very fine in the results between the predictions of grid G2 and G3 were seen.
mesh was formed in the close proximity of the cylinder. The mesh was Therefore, grid G2 consisting of 200 grid points on the surface of the
progressively made coarser away from the cylinder using a specific cylinder, ~ 40,000 CVs in the entire computational domain, the smallest
growth rate parameter. A two-dimensional double precision (2ddp), grid size of ~6 × 10− 5 m and face areas in the range (~5 × 10− 5 and 4.8
transient, laminar and pressure based segregated solver was used with × 10− 4 m2) was considered to be sufficient to capture the prevailing
second order upwind scheme for discretizing the convective terms in the steep velocity and temperature gradients in the boundary layers and for
momentum and energy equations. To incorporate the natural convection the present numerical results to be free from grid effects. Subsequently, a
effects, the gravity vector was set as − 9.81 m/s2 in the y-direction. The systematic time-step test was conducted at the four different time-step
central differencing scheme was selected for the diffusion terms in the sizes, namely, Δt = 0.05 s, 0.1 s, 0.2 s and 0.5 s. The time-step of 0.1 s
governing equations. The semi-implicit method for the pressure linked was found to be suitable to efficiently capture the melt interface and
equations (SIMPLE) was used for pressure-velocity coupling and the temporal evolutions in the melting regions. Hence, the results reported
PRESTO (Pressure Staggering Option) scheme was selected for the herein are based on grid G2 and a time-step size of 0.1 s. The adequacy
spatial discretization and interpolation of pressure terms [65]. A second- for the selection of these computational parameters is further demon­
order accurate, implicit scheme was used for time-discretization. The strated in the ensuing section by making some benchmark comparisons.
under-relaxation factors for density, momentum, pressure correction,
thermal energy and liquid fraction were taken as 1, 0.7, 0.3, 1 and 0.9 4. Results and discussion
respectively. The maximum number iterations per time step was set as
400. The simulations were deemed to stop when the residuals of prim­ In this study, extensive numerical results have been reported for
itive variables satisfy the relative convergence (tolerance) criterion of melting of the lauric acid as PCM around a heated horizontal circular
⃒( n+1 ⃒
⃒ Γ − Γn )/Γn+1 ⃒⩽10− 5 for velocity and pressure and ⩽10− 8 for tem­ cylinder confined in a square enclosure. The geometrical position of the
perature at each time step between the two consecutive iterations, cylinder has been systematically varied as δ = 0 (centre), δx = − L/4
where n is the iteration number and Γ stands for the independent vari­ (left), δy = ±L/4 (top/bottom) for the four symmetric locations; δx = − L/
ables (V, p and T). Within this framework, the corresponding velocity 4, δy = +L/4 (top asymmetric) and δx = − L/4, δy = − L/4 (bottom
and temperature field were stabilized at least up to the four significant asymmetric) for the two asymmetric locations with respect to the sur­
digits. No notable changes in the results were observed on further rounding walls of the enclosure. These distinct locations inside the
improving (lowering) this criterion and warranted excessive computa­ enclosure delineate the effects of position of the cylinder on the melt
tion time. Computations were carried out on a workstation with Intel characteristics and the rate of convective heat transfer within the liquid
Xeon W-2155 processor (3.3 GHz frequency) having 64 GB of RAM. The melt. The value of the non-dimensional numbers, namely, the Rayleigh
simulations were run until the whole PCM inside the enclosure is melt number and Prandtl number are within the range for the steady and
and the total computational time to achieve the complete melting (λ = 1) laminar flow conditions to prevail. It is noteworthy to recall here that
varied from ~190 h (minimum for bottom position) to ~520 h due to the imposition of the identical boundary conditions on the walls

6
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

Mahdaoui et al. [57] Present results cylinder (at TH = 311 K) immersed inside a rectangular enclosure filled
with the Gallium as the PCM (Tm = 302.93 K). Our results are in close
match with those reported by Mahdaoui et al. [57], see Fig. 2b. Fig. 2a
t = 60 sec t = 60 sec
also suggests the existence of an optimum value of the mushy zone
parameter (Amush ≈ 2 × 105 and 5 × 105) which closely matches with the
experimental results. This value not only depends on the properties of
the PCM and physics of the flow but also on the orientation/configu­
ration of the computational setup. Therefore, a prudent selection of the
value of Amush has been the subject of numerous discussions to treat the
fluid flow governed by the natural convection within the mushy region
[64]. The impact of varying the mushy-zone constant (Amush) on the
predicted numerical results has been presented in ensuing sections.
Based on the preceding comparisons, the numerical results presented in
the ensuing sections are believed to be reliable within ±3–4% deviation.
t = 200 sec t = 200 sec
4.2. Temperature contours and velocity vectors

It is customary to visualize the velocity and temperature variations

during the melting process to delineate the recirculation regions,
strength of flow (strong/weak) and hot and cold spots within the fluid
melt. Fig. 3 and Fig. 4 show the representative results of temperature
contours and velocity vectors, respectively, for various symmetric lo­
cations of cylinder within the square enclosure at different time instants
of the melting process. In the early stages of heating, the PCM is in solid
state and therefore conduction is the predominant mode of heat transfer
(t = 1 min, Fig. 3). As time progresses, temperature of the PCM close to
t = 800 sec t = 800 sec the cylinder surface approaches the melting range and the melting starts.
The temperature of the liquid PCM (melt) is higher in the vicinity of the
cylinder as compared to the liquid PCM slightly away from the cylinder.
This temperature difference induces the density stratification within the
liquid PCM and thus the hot PCM melt moves upwards and forms a rising
plume above the cylinder as can be seen in Fig. 3. The ascending hot
PCM melt is then replaced by the PCM melt at relatively low tempera­
ture from the sides of the cylinder. At this stage, the buoyant forces
overcome the viscous forces and the natural convection initiates within
the PCM melt. Consequently, weak advective currents in form of two
counter rotating cells symmetrical about the y-axis are visible above the
Fig. 2b. Comparison of the structure of solid-liquid interface around a circular surface of the heated circular cylinder (see Fig. 4 at t = 5 min). As the
cylinder embedded in a rectangular enclosure filled with the lauric acid as PCM melting progresses, the buoyancy is intensified and large volume of the
with the results reported by Mahdaoui et al. [57]. liquid PCM is collected above the cylinder which is evident from the
distorted shape of the temperature contours. Once the thermal plume
of the enclosure, the corresponding positions of the cylinder on the right touches the upper wall of the enclosure, the accumulated liquid PCM at
half part of the computational domain, such as (δx = +L/4), (δx = +L/4, the top of the enclosure starts descending. Thereafter, the descending
δy = +L/4) and (δx = +L/4, δy = − L/4), are analogous to their left PCM melt comes in contact with the surrounding solid PCM which is at a
counterparts, therefore, these have not been included here for the sake lower temperature. The surrounding PCM in its solid state acts as a cold
of brevity. Prior to presenting new results, it is useful to ensure the ac­ wall for the ascending fluid melt. Owing to the feeble convection cur­
curacy and reliability of the numerical scheme and results by presenting rents in the initial stage (t = 5 min), velocity vector and temperature
some benchmark comparisons. contour plots show similar results irrespective of the geometric position
of the cylinder within the enclosure, i.e., for all values of δ. The weak
advection at this stage is also manifested in form of symmetric flow and
4.1. Validation of results temperature contours even for the asymmetric placement of the cylinder
within the enclosure, i.e. for δx = − L/4 (left). Similar velocity vector
As noted earlier, while there are limited results for the configuration plots were observed for asymmetric cases (top asymmetric and bottom
studied herein, reliable results are available for the melting of PCM in­ asymmetric), at least qualitatively, thus these are not shown here for the
side a vertically and horizontally oriented rectangular enclosure. Fig. 2a sake of brevity.
shows the comparison of melting rate of lauric acid in two different As time increases, the melting rate increases which is accompanied
orientations (horizontal and vertical) of the rectangular enclosure with by the progressive expansion in the liquid PCM region. Here, the
the experimental results of Kamkari et al. [46] and the numerical study convective currents become stronger leading to an improved heat
of Fadl and Eames [42]. Two different values of Amush, i.e., 105 and 106 transfer rate. The asymmetric effects become visible from the velocity
have been considered. Evidently, the present results are in fair agree­ vector profiles and temperature contours as the expansion of the PCM
ment with the reported experimental (~5–6% error) and numerical melt increases and melt interface touches the left wall of the enclosure.
studies (~1–2% error). The small differences in the results can be safely At the stage, the symmetry in the velocity and temperature profiles
attributed to different mesh refinement and grid structures used in Fadl about the y-axis is lost and more melting and energy transfer occurs left
and Eames [42]. Similarly, the results of solid-liquid interfaces (melting to the cylinder as compared to the right. These trends are depicted in
front) have been validated with the numerical predictions of Mahdaoui Fig. 3 and Fig. 4 for the asymmetric placement of the cylinder close to
et al. [57] who investigated melting and heat transfer around a circular the left wall of the enclosure (δx = − L/4) at t = 20 min. For other cases, a

7
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

Fig. 3. Temporal evolution of isotherm contours for different symmetric positions of the cylinder within the enclosure.

symmetric spread of the fluid region and the thermal plume can be that the rate of hear transfer and thus the fraction of the PCM melted in
observed. the enclosure is anticipated to be the highest for this case (δy = − L/4),
As expected, the buoyancy induced flow currents are seen to be however, the melt characteristics have been discussed in more detail in
stronger in the region above the heated cylinder, whereas the PCM filled the ensuing section. Furthermore, once the temperature of the PCM in
in the enclosure below the cylinder remains in the solid state at a tem­ the upper half of enclosure becomes higher than the liquidus tempera­
perature below the melting range where conduction dominates over ture (T + ΔT), the whole PCM above the cylinder is melted and there­
convection. Therefore, the more is the phase change material above the after the energy is transferred from the heated cylinder to the PCM by
cylinder, the greater is the extent and strength of the thermal plume and sensible heating only as illustrated in Fig. 3 for t = 80 min.
thus the higher is the rate of melting. Clearly, the extent of the fluid
circulation in the PCM melt is seen to be maximum for the case δy = − L/
4.3. Melt fraction
4, i.e., when the cylinder is close to the bottom surface of the enclosure.
Intuitively, the trends in the temperature and velocity contours suggest
Fig. 5 and Fig. 6 illustrate the temporal evolutions in the melt

8
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

Fig. 4. Representative velocity vectors for different symmetric positions of the cylinder inside the enclosure at different time instants (a) t = 5 min; (b) t = 20 min; (c)
40 min.

fraction contours (constant liquid fraction lines) for different geometric interface remains uniform and parallel to the cylinder surface. Also,
placements of the cylinder inside the square enclosure, i.e., for the four there exists a fore and aft symmetry about the x-axis for all values of δ.
symmetric configurations of the cylinder, namely, δ = 0 (centre), δy = As the temperature of the PCM in the close proximity of the cylinder
+L/4 (close to the top wall), δy = − L/4 (close to the bottom wall) and δx increases and comes in the melting range (Tm − ΔT) ≤ T ≤ (Tm + ΔT),
= − L/4 (close to the left wall) as shown in Fig. 5, and the two asym­ the melting commences and the liquid fraction in a cell takes a value
metric configurations, namely top asymmetric (δx = − L/4, δy = +L/4) between 0 ≤ λ ≤ 1 depending upon the value of the temperature. It is
and bottom asymmetric (δx = − L/4, δy = − L/4) as shown in Fig. 6. These interesting to mention here that at this stage, the liquid PCM is in
results can be visualized similar to that of temperature contours dis­ transition region and still stagnant due to high viscous forces. Once the
cussed in the aforementioned section. This is simply due to the fact that temperature of liquid PCM increases beyond the liquidus temperature
in the melting process predicted by the enthalpy-porosity technique, the where λ = 1, the natural convection starts and upward moving
temperature is a prime factor in determining the position of the solid- convective currents becomes stronger. This causes a concomitant
liquid interface and hence the melt fraction. Due to this, only a key expansion in the solid-liquid interface in the direction of buoyancy, i.e.,
features have been recapitulated here. the extend of liquid zone increases. For δx = − L/4 (close to left wall), the
At early phase (t ≈ 1 min, Fig. 5), heat transfer from the heated growing liquid PCM capsule intersect with the left wall and the sym­
cylinder to the surrounding solid PCM is mainly due to conduction and metry about the y-axis is lost at this instant (t = 20 min). The fluid melt
melt fraction contours are concentric circles and closely follow the shape then slowly drifts towards the top and then to the right thereby
of the cylinder. Until this mode of heat transfer prevails, the solid-liquid encompassing the entire region above the cylinder (t = 80 min). When

9
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

Fig. 4. (continued).

the whole PCM above the cylinder is fully melted (λ = 1), the energy is δy = − L/4, i.e., for the cylinder placed close to the bottom wall of the
transferred to the PCM beneath the cylinder by conduction were the enclosure. It achieves the completely charged condition (λ = 1) at
feeble convective currents lead to slower progression of the solid-liquid about t = 130 min. As the cylinder moves upwards, the time required to
interface and finally the whole PCM inside the enclosure comes in liquid reach the fully-melt state increases thereby showing the reduced rate
state. For instance, at t ≈ 130 min, the entire PCM confined in the of melting. For example, the slowest melting rate is observed for the
enclosure melts for the bottom configuration of the cylinder (δy = − L/4), case of δy = +L/4 (top) where the PCM takes ~480 min to reach the
as can be seen in Fig. 5. It is worthwhile to notice here that these results fully liquid state. In other words, at any instant, the melt fraction of the
are in line with those reported by Mahdaoui et al. [57] for a circular PCM is higher for the bottom positioning of the cylinder inside the
cylinder positioned inside the rectangular enclosure, at least qualita­ enclosure, such as at t ≈ 130 min, the melt fraction of the PCM for the
tively. In continuation to this, the trends displayed in Fig. 6 for the two configuration δy = − L/4 (bottom) is equal to 1 (complete liquid state)
asymmetric positions of the cylinder inside the enclosure considered which is ~80% and ~25% higher than the corresponding values of the
herein are the blend of the observations seen for the top, left and bottom melt fraction for δy = +L/4 (top); and δ = 0 (centre) or δx = − L/4 (left),
configurations of the cylinder. At the early stages (t ~ < 5 min), both respectively. Therefore, the vertical distance of the cylinder from the
asymmetric cases predict the same melting and thermal behaviour as base of the enclosure plays a key role in deciding the progression of the
also observed in case of all the four symmetric positions of the cylinder solid-liquid interface and determining the time needed for the entire
such as centre, left, top and bottom (see Fig. 5). The major changes PCM to melt. Interestingly, the results shown by top asymmetric
embark once the time progresses beyond 5 min. For instance, in the top configuration is the same as shown by the top configuration. Similarly,
asymmetric case, the thermal plume and thus the melting interface first the bottom asymmetric and bottom configurations of the cylinder also
touches the top wall of the enclosure and then spreads along the wall and resulted the same results, thereby one may conclude that different
moves downward. Underneath the cylinder, the PCM melts primarily by horizontal positions on the cylinder produce nearly indistinguishable
conduction and thus melting is very sluggish. On the other hand, in case results and thus one may place the cylinder at any horizontal location
of bottom asymmetric position, the rising thermal plume (and thus the (left, centre or right). It is the vertical position of the cylinder which
melting interface) first touches the left wall and continues to rise up until mainly governs the overall melting and heat transfer. These observa­
it touches the top wall and then spreads in the right part of the domain. tions can be justified simply due to the presence of same amount of the
Since the PCM has more volume above the cylinder here, melting is PCM above the cylinder. It is also worthwhile to mention that this
accelerated due to natural convection until the liquid PCM occupies the study only deals with the melting of the PCM around an isolated cyl­
entire space above the cylinder. After this, conduction becomes domi­ inder placed inside the enclosure. When more than one cylinders are
nant and therefore, a slower rate of melting is observed. encountered, the inter-spacing between the two cylinders also plays a
Fig. 7a shows the variation in the values of melt fraction (averaged vital role in describing the physics of the flow due to the varying levels
over the entire domain under consideration) of the PCM with time. of interaction among the melting interfaces [66].
Such plots are necessary to depict the rate of melting in the phase From an application standpoint, the present numerical results for
transition process. The fastest melting rate is observed for the case of melt fraction (λ) are correlated as a function of the charging (melting)

10
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

Fig. 5. Temporal evolution of melt fraction contours for different symmetric positions of the cylinder within the enclosure.

time t as follows: 4.4. Energy storage

b
at
λ = if 0⩽t⩽tcrit ; (16a) The amount of energy storage within the PCM has been calculated at
c + tb
a particular time instant by subtracting the initial energy content within
λ = 1 if t > tcrit (16b) the PCM from the energy content at any referred time. These results are
important in elucidating the storage capacity and performance evalua­
In the above equations, tcrit is the critical time (in min) corresponding tion of a PCM. Fig. 7b shows the time variation in the amount of stored
to λ = 1 (or total time required for complete melting of the PCM). Table 2 energy for different positions of the circular cylinder. During early stages
summarizes the values of the critical time and the empirically fitted of melting, the amount of energy stored in the PCM increases with the
constants in Eq. (16a) for all six distinct configurations of the cylinder same rate irrespective of the location of the cylinder inside the enclo­
inside the enclosure. sure. Evidently, the rate of energy transfer to the PCM is the fastest for δy
= − L/4 (bottom) configuration and the slowest for δy = +L/4 (top)

11
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

Fig. 6. Asymmetric positions of the cylinder within the enclosure: (a) melt fraction contour; (b) temperature contour.

configuration. As the cylinder moves upwards in the enclosure, the rate 4.5. Average Nusselt number
of energy storage within the PCM decays. Therefore, the time taken by
the PCM to reach the same energy levels is ~4 times and ~2 times higher Notwithstanding with the fact that the cylinder is isothermally
for the cylinder at position δy = +L/4 (top) and δx = − L/4 (left) as heated, the temperature gradient and thus the rate of heat transfer varies
compared to that taken at the position δy = − L/4. Similarly, at any fixed locally along the surface of the cylinder. The surface-averaged Nusselt
time instant, the bottom configuration stores more energy in comparison number is expected to be dependent on both the melting time and
to all other locations of the cylinder. For instance, at time t ≈ 130 min, geometrical positioning of the cylinder within the enclosure. Fig. 8
the amount of energy stored in the liquid PCM for the cylinder position shows the time variation of the surface average Nusselt number for the
at δy = − L/4 is ~65% and ~20% greater than that stored for the posi­ six different positions of the cylinder within the enclosure considered
tions δy = +L/4; and δ = 0 or δx = ±L/4, respectively. The plots for δx = herein. During melting, the convective heat transfer rate is limited by
− L/4 (left) and δ = 0 (centre) are very close to each other due to the the rate of melting of the PCM. As postulated, during early phase of the
same distance of the cylinder from the bottom surface of the enclosure melting process, i.e., for time t < 5 min., the conduction dominates over
for these cases. Moreover, the trends on energy storage rate for asym­ convection. As the melting progresses, the buoyancy-induced flow
metric configurations (top asymmetric and bottom symmetric) resemble strengthens and the convection contributes increasingly to the overall
to that already discussed in previous section on the rate of melting in heat transfer rate. The Nusselt number values increase as the convective
Fig. 7a. currents becomes stronger in the region above the heated cylinder and
attain their maxima at a time corresponding to the maximum rate of
melting above the cylinder. For instance, the maxima occur at time t ≈

12
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

attain an asymptotic constant value of zero corresponding to the negli­

gible temperature-difference available in the liquid PCM. This steady
state condition is quickly achieved in case of bottom configuration of the
cylinder (δy = − L/4) because of the highest melting rate. A thorough
inspection of Fig. 8 suggests that the Nusselt number curve for δy = − L/4
(bottom) at t ≈ 130 min crosses all other curves and takes the minimum
value. This switch over occurs when Nu ≈ 1 showing the comparable
strengths of both convective and conductive heat transfer. All in all,
suffice it to conclude here that the lower is the distance between the hot
cylinder and the bottom of the enclosure, the greater will be the heat
transfer rate during the melting process. The Nusselt number results for
the two asymmetric configurations, namely, top asymmetric and bottom
asymmetric are seen to approach top and bottom configuration results,
respectively, thereby confirming our previous observations that hori­
zontal locations give rise to similar results due to the same height from
the bottom wall of the enclosure.

4.6. Effect of mushy zone parameter (Amush)

Fig. 7a. Time variation of melt fraction for different geometric positions of the
cylinder inside the enclosure.
Suffice it to reiterate here that the results presented in this work are

Table 2
Values of constants used in Eq. (16a) along with the corresponding errors.
Configuration tcrit (min) a b c % Error

Average Maximum

Top 480 2.433 0.640 74.80 1.5 16.5

Left 270 1.154 1.279 222.6 4.0 16.7
Centre 255 1.138 1.339 261.4 3.2 19.7
Bottom 130 1.133 2.290 7790.7 3.6 21.3
Top Asymmetric 480 2.399 0.645 75.01 1.8 13.5
Bottom Asymmetric 130 1.255 1.823 1636.7 3.1 19.6

Fig. 7b. Temporal variation of the total energy stored in the enclosure at
different positions of cylinder inside the enclosure. Fig. 8. Time history of the average Nusselt number for different geometric
placements of the cylinder within the square enclosure.

18 min for the top (δy = +L/4) and left (δx = − L/4) positions; at t ≈ 38
min and t ≈ 45 min for the centre (δ = 0) and bottom (δy = − L/4) po­ based on the sufficiently large value of the mushy zone constant Amush =
sitions of the cylinder, respectively. Due to the upward buoyant cur­ 105 kg/(m3s), therefore, it is wise to study the influence of Amush on the
rents, the melting in the PCM above the cylinder is mainly governed by development of solid-liquid interface and the rate of melting. A few
natural convection (fast heat transfer) in contrast to the PCM below the simulations have also been run to briefly depict this dependence. Figs. 9,
cylinder where conduction (slow heat transfer) solely causes the melting 10a and 10b show the representative results for the central location of
to occur. Therefore, once most of the PCM above the cylinder is melted, the cylinder (δ = 0) at the three different values of the mushy zone
the Nusselt number values monotonically decrease due to the continu­ parameter, namely, Amush = 104, 105 and 106 kg/(m3s). Clearly, the
ously decreasing temperature difference in the PCM which deteriorates melting rate and hence the convective heat transfer rate slows down (i.
the strength of buoyant currents (the velocity Vy ~ ΔT1/2) and thus the e., delayed melting) as the value of Amush gradually increases from 104 to
rate of heat transfer. Finally, the Nusselt number values are seen to 106 kg/(m3s). In other words, the smaller values of the Amush over predict

13
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

Fig. 9. Effect of mushy zone parameter (Amush) on the solid-liquid interface for δ = 0 (centre) at the two different instants (a) t = 20 min; (b) t = 40 min.

Fig. 10a. Influence of mushy zone parameter (Amush) on the time variation of Fig. 10b. Influence of the mushy zone parameters (Amush) on the time variation
melt fraction for δ = 0. of the total energy stored in the PCM for δ = 0.

Fig. 11. Cylinders of various cross-sections (of same perimeter) considered in this study.

14
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

Fig. 12. Temporal evolution of melt fraction contours for different shapes of the cylinder within the square enclosure for δ = 0 (cylinder located at the centre).

the extent of melting in the PCM under otherwise identical conditions, as The aforementioned trends are consistent with the findings of Fadl
evident in Fig. 9. However, the effects of the different values of the and Eames [42] and can be reasonably justified with the fact that the
mushy zone parameter are more prominent for the time range where momentum source term, Sm, in Eq. (2) acts as a damping term in the
convection heat transfer subsides the conduction. Significant effects of mushy region. The increasing values of Amush add more volume into the
Amush are visible in the range of t ≈ 5 to 130 min. For time instants t < 5 mushy region where the melted PCM is mostly static (small fluid ve­
min (Fig. 10a), the conduction is the dominant mode of heat transfer and locities). This adversely affects the natural convection and thus reduces
hence no effects of the Amush are visible. Also, for t > ~130 min. the the overall heat transfer and melting extent in the domain. Therefore, a
effects of Amush subside due to the exhausted melting rate in the PCM proper selection of Amush is crucial for the accurate prediction of melting
above the cylinder. In the PCM below the cylinder, conduction governs and heat transfer in the latent heat thermal energy storage (LHTES)
the rate of heat transfer and accordingly the Amush shows virtually no systems utilizing PCMs. Essential is to state here that the numerical
influence on the solid-liquid interface in this region. The results for the predictions of results by a sufficiently large value of Amush should be
total energy storage within the liquid PCM (plotted in Fig. 10b) are compared with the available experimental data to ascertain the viability
analogous to those observed in Fig. 10a. of the solution.

15
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

energy within the PCM are observed at any time instants. Numerical
predictions show that among all the cases, the equilateral triangular
cylinder witnesses the maximum melting and the energy storage rate in
the PCM whereas the horizontally oriented rectangular cylinder results
the minimum values. This is due to the abrupt changes in the body
contour present on the surface of the triangular cylinder which causes
thermal boundary layers to bend thereby increasing the rate of hear
transfer.
It is readily acknowledged that at sufficiently large Rayleigh
numbers, the effects of cylinder body contours are expected to be sig­
nificant due to the relatively strong buoyant currents and distortion of
the thermal plumes. Therefore, this part is a matter of detailed discus­
sion and a comprehensive analysis over these effects will be presented in
the future work. In addition to this, unlike the heat transfer fluids where
the wall effects (size ratio of the cylinder to the enclosure) also play a
vital role in free-convection heat transfer problems, in case of PCM not
only the strength of buoyancy (Grashof number) may change with the
blockage ratio but also the amount (or mass) of PCM. The smaller is the
Fig. 13a. Time variation of melt fraction for variously shaped cylinders within area available in the annulus, the lesser is the amount of the PCM filled
the square enclosure for δ = 0 (centre). inside it and thus the faster is the melting. At a constant value of the
Grashof number (based on cylinder diameter) and for a range of
confinement (blockage) levels, the melting behaviour is straightforward.
However, the melting behaviour of the PCM under various levels of
confinement generates a different class of problems where both the
Grashof number and blockage ratio may change simultaneously. A
thorough and rigorous analysis in this regard is thus required. This
important section will be presented in detail in our future studies.

5. Conclusions

The current study presented extensive new results and a through

discussion on the melt characteristics and thermal development around
an isothermal 2D circular cylinder placed in a square enclosure filled
with the PCM (lauric acid) under the laminar natural convection heat
transfer regime. Time-dependent governing equations of momentum
and energy based on enthalpy-porosity technique were solved numeri­
cally using ANSYS Fluent (version 19.2). The onset of melting was
modelled by incorporating a momentum source term in the Navier-
Stokes equation. The influence of the location of the cylinder within
Fig. 13b. Time variation of energy stored in the PCM for variously shaped the enclosure was examined for the four symmetric configurations,
cylinders within the square enclosure for δ = 0 (centre). namely, at centre (δ = 0), close to the left wall (δx = − L/4), close to the
bottom wall (δy = − L/4), close to the top wall (δy = +L/4); and two
4.7. Effect of shape of the cylinder asymmetric configurations, namely, top asymmetric (δx = − L/4, δy =
+L/4) and bottom asymmetric (δx = − L/4, δy = − L/4). The represen­
Further attempts have been made to elucidate the effects of cylinder tative results of velocity and temperature fields were visualized in terms
shape on the rate of melting and energy storage. A close comparison has of velocity vectors, temperature contours, solid-liquid interface, rate of
been addressed here for seven shapes of the cylinder other than the melting and energy storage and the average Nusselt number. Broadly,
circular, as listed in Fig. 11. It is worthwhile to mention here that the the melting in the PCM above the cylinder is governed by convection
shape of the cylinder is varied keeping the perimeter as constant so to whereas melting in the PCM beneath the cylinder is primarily because of
achieve the same heat transfer area of the cylinder. The position of conduction only. A horizontal movement in the cylinder position along
cylinder is also kept fixed at the centre of the enclosure (δ = 0). Only the x-axis exhibited virtually the same (within ±2%) rate of melting and
limited simulations have been carried out maintaining the same number energy storage. At any instant, the highest melting and energy storage
of grid points on the cylinder and nearly the same number of cells in the within the PCM was noticed for δy = − L/4, i.e., the configuration where
computational domain as obtained for the case of a circular cylinder the hot cylinder is closest to the bottom of the enclosure. Alternatively,
(refer section 3). Figs. 12, 13a and 13b show the representative plots for the larger is the PCM available above the heated object, the higher is the
the development of melting front, melt fraction and total energy storage rate of heat transfer fostering the rate of melting. Thus, it is possible to
in the PCM for a range of shapes of the cylinder, respectively. Owing to achieve a significant improvement in the rate of melting and thermal
the weak buoyancy at moderately low value of the Rayleigh number energy storage within the PCM using simple geometric modifications in
considered here, the shape of the heated cylinder only shows a little the latent heat thermal energy storage (LHTES) systems. Moreover, the
influence on the results. Initially, due to the slower rate of heat transfer larger values of Amush were seen to reduce the melting time (delayed
(primarily by conduction only), the plots are seen to overlap (Figs. 13a melting), however, these effects were merely visible at those instants
and 13b) and morphology of the solid-liquid interface is nearly indis­ where heat transfer is dominated by convection. Finally, the shape of the
tinguishable from each other (Fig. 12). As the melting gradually in­ cylinder also affects the melting behaviour, though a little variation in
creases, these effects are visible and ~5–6% differences between the the results were observed due to small Rayleigh number considered
maximum and the minimum value of the melt fraction and the absorbed herein. A close agreement was seen to exist between the present results
and those reported in the literature.

16
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

Declaration of Competing Interest [28] G. Comini, S.D. Guidice, R.W. Lewis, O.C. Zienkiewicz, Finite element solution of
non-linear heat conduction problems with special reference to phase change, Int. J.
Numer. Methods Eng. 8 (1974) 613–624.
The authors declare that they have no known competing financial [29] W. Rolph III, K.J. Bathe, An efficient algorithm for analysis of nonlinear heat
interests or personal relationships that could have appeared to influence transfer with phase change, Int. J. Numer. Methods Eng. 18 (1982) 119–134.
the work reported in this paper. [30] J.S. Hsiao, B.T.F. Chung, An efficient algorithm for finite solution to two
dimensional heat transfer with melting and freezing, ASME paper 84-HT-2,
presented at the 22nd national heat transfer conference, Niagara Falls, August,
Acknowledgements 1984.
[31] D. Gobin, Role de la convection thermique dans les processus de fusion-
solidification, Ecole d_ete, GUT-CET, Mod_elisation num_erique en thermique,
The corresponding author gratefully acknowledges the Department Institut d_etudes scientifiques de Carg_ese (1992).
of Science & Technology (DST), Govt. of INDIA for the INSPIRE Faculty [32] D.W. Hahn, M.N. Ozisik, Phase Change Problems, Chapter 12. Heat Conduction,
3rd ed., John Wiley & Sons Inc., New Jersey, 2012.
Research Grant for the period 2018-2023 to carry out this work.
[33] E.R.G. Eckert, R.J. Goldstein, W.E. Ibele, S.V. Patankar, T.W. Simon, P.
J. Strykowski, K.K. Tamma, T.H. Kuehn, A. Bar-Cohen, J.V.R. Heberlein, D.
References L. Hofeldt, J.H. Davidson, J. Bischof, F. Kulacki, Heat transfer––a review of 1994
literature, Int. J. Heat Mass Transfer 40 (1997) (1994) 3729–3804.
[34] K.A.R. Ismail, R. Stuginsky Jr., A parametric study on possible fixed bed models for
[1] https://www.c2es.org/content/renewable-energy/.
PCM and sensible heat storage, Appl. Therm. Eng. 19 (1999) 757–788.
[2] H. Nazir, M. Batool, F.J.B. Osorio, M. Isaza-Ruiz, X. Xu, K. Vignarooban, P. Phelan,
[35] J. Yoo, B. Rubinsky, Numerical computation using finite elements for the moving
A.M.K. Inamuddin, Recent developments in phase change materials for energy
interface in heat transfer problems with phase change transformation, Numer. Heat
storage applications: a review, Int. J. Heat Mass Transfer 129 (2019) 491–523.
Transfer 6 (1983) 209–222.
[3] A. Hussain, S.M. Arif, M. Aslam, Emerging renewable and sustainable energy
[36] R.M. Furzerland, A comparative study of numerical methods for moving boundary
technologies: state of the art, Renew. Sustain. Energy Rev. 71 (2017) 12–28.
problems, J. Inst. Math. Appl. 26 (1980) 411–429.
[4] A. Sharma, V.V. Tyagi, C.R. Chen, D. Buddhi, Review on thermal energy storage
[37] N. Shamsunder, E. Sparrow, Analysis of multidimensional phase change via the
with phase change materials and applications, Renew. Sustain. Energy Rev. 13
enthalpy model, J. Heat Transfer Trans. ASME 19 (1975) 333–340.
(2009) 318–345.
[38] S.E. Hibbert, N.C. Markatos, V.R. Voller, Computer simulation of moving interface,
[5] Y.B. Tao, Y.-L. He, A review of phase change material and performance
convective, phase change process, Int. J. Heat Mass Transfer 31 (1988) 1785–1795.
enhancement method for latent heat storage system, Renew. Sustain. Energy Rev.
[39] A.D. Brent, V.R. Voller, K.J. Reid, Enthalpy-porosity technique for modeling
93 (2018) 245–259.
convection- diffusion phase change: application to the melting of a pure metal,
[6] J.P. da Cunha, P. Eames, Thermal energy storage for low and medium temperature
Numer. Heat Transfer 13 (1988) 297–318.
applications using phase change materials– a review, Appl. Energy 177 (2016)
[40] L.W. Hunter, J.R. Kuttler, The enthalpy method for heat conduction problems with
227–238.
moving boundaries, J. Heat Transfer Trans. ASME 111 (1989) 239–242.
[7] M.M. Farid, A.M. Khudhair, S.A.K. Razack, S. Al-Hallaj, A review on phase change
[41] S. Arena, E. Casti, J. Gasia, L.F. Cabeza, G. Cau, Numerical simulation of a finned
energy storage: materials and applications, Energy Convers. Manage. 45 (2004)
tube LHTES system: influence of the mushy zone constant on the phase change
1597–1615.
behaviour, Energy Proc. 126 (2017) 517–524.
[8] B. Zalba, J.M. Marı ́n, L.F. Cabeza, H. Mehling, Review on thermal energy storage
[42] M. Fadl, P.C. Eames, Numerical investigation of the influence of mushy zone
with phase change: materials, heat transfer analysis and applications, Appl. Therm.
parameter Amush on heat transfer characteristics in vertically and horizontally
Eng. 23 (2003) 251–283.
oriented thermal energy storage system, Appl. Therm. Eng. 151 (2019) 90–99.
[9] Y. Dutil, D.R. Rousse, N.B. Salah, S. Lassue, L. Zalewski, A review on phase-change
[43] A.A. Al-abidi, S.B. Mat, K. Sopian, M.Y. Sulaiman, A.Th. Mohammed, CFD
materials: mathematical modeling and simulations, Renew. Sustain. Energy Rev.
applications for latent heat thermal energy storage: a review, Renew. Sustain.
15 (2011) 112–130.
Energy Rev. 20 (2013) 353–363.
[10] H. Akeiber, P. Nejat, M.Z.A. Majid, M.A. Wahid, F. Jomehzadeh, I.Z. Famileh, J.
[44] K. Kant, A. Shukla, A. Sharma, Biwole, Melting and solidification behaviour of
K. Calautit, B.R. Hughes, S.A. Zaki, A review on phase change material (PCM) for
phase change materials with cyclic heating and cooling, J. Energy Storage 15
sustainable passive cooling in building envelopes, Renew. Sustain. Energy Rev. 60
(2018) 274–282.
(2016) 1470–1497.
[45] L. Zeng, J. Lu, Y. Li, W. Li, S. Liu, J. Zhu, Numerical study of the influences of
[11] E. Alehosseini, S.M. Jafari, Micro/nano-encapsulated phase change materials
geometry orientation on phase change material’s melting process, Adv. Mech. Eng.
(PCMs) as emerging materials for the food industry, Trends Food Sci. Technol. 91
9 (2017) 1–11.
(2019) 116–128.
[46] B. Kamkari, H. Shokouhmand, F. Bruno, Experimental investigation of the effect of
[12] F.L. Tan, C.P. Tso, Cooling of mobile electronic devices using phase change
inclination angle on convection-driven melting of phase change material in a
materials, Appl. Therm. Eng. 24 (2004) 159–169.
rectangular enclosure, Int. J. Heat Mass Transfer 72 (2014) 186–200.
[13] X.Q. Wang, A.S. Mujumdar, C. Yap, Effect of orientation for phase change material
[47] A.O. Elsayed, Numerical investigation on PCM melting in triangular cylinders,
(PCM)-based heat sinks for transient thermal management of electric components,
Alexandria Eng. J. 57 (2018) 2819–2828.
Int. Commun. Heat Mass Transfer 34 (2007) 801–808.
[48] R.D. White, A.G. Bathelt, W. Leidenfrost, R. Viskanta, Study of heat transfer and
[14] A. Sari, Thermal reliability test of some fatty acids as PCMs used for solar thermal
melting from a cylinder imbedded in a phase change material, ASME Paper No. 77-
latent heat storage applications, Energy Convers. Manage. 44 (2003) 2277–2287.
HT-42, 1977.
[15] M. Kenisarin, K. Mahkamov, Solar energy storage using phase change materials,
[49] L.S. Yao, F.F. Chen, Effects of natural convection in the melted region around a
Renew. Sustain. Energy Rev. 11 (2011) 1913–1965.
heated horizontal cylinder, J. Heat Transfer 102 (1980) 667–672.
[16] D. Zhou, C.Y. Zhao, Y. Tian, Review on thermal energy storage with phase change
[50] A.G. Bathelt, R. Viskanta, Heat transfer at the solid–liquid interface during melting
materials (PCMs) in building applications, Appl. Energy 92 (2012) 593–605.
from a horizontal cylinder, Int. J. Heat Mass Transfer 23 (1980) 1493–1503.
[17] A. Sharma, A. Shukla, C.R.R. Chen, S. Dwivedi, Development of phase change
[51] K.C. Cheng, H. Inaba, R.G. Gilpin, Effect of natural convection on ice formation
materials for building applications, Energy Build. 64 (2013) 403–407.
around an isothermally cooled horizontal cylinder, J. Heat Transfer 110 (1988)
[18] S. Ghani, E.M.A.A. Elbialy, F. Bakochristou, S.M.A. Gamaledin, M.M. Rashwan, The
931–937.
effect of forced convection and PCM on helmets’ thermal performance in hot and
[52] H. Rieger, V. Projahan, H. Beer, Analysis of the heat transport mechanisms during
arid environments, Appl. Therm. Eng. 111 (2017) 624–637.
melting around a horizontal circular cylinder, Int. J. Heat Mass Transfer 25 (1982)
[19] S. Mondal, Phase change materials for smart textiles – an overview, Appl. Therm.
137–147.
Eng. 28 (2008) 1536–1550.
[53] J. Prusa, L.S. Yao, Melting around a horizontal heated cylinder: Part I–
[20] A.A. Al-Abidi, S.B. Mat, K. Sopian, M. Sulaiman, C.H. Lim, A. Th, Review of thermal
perturbation and numerical solutions for constant heat flux boundary condition,
energy storage for air conditioning systems, Renew. Sustain. Energy Rev. 16 (2012)
J. Heat Transfer 106 (1984) 376–384.
5802–5819.
[54] J. Prusa, L.S. Yao Yao, Melting around a horizontal heated cylinder: Part
[21] A. Shukla, A. Sharma, M. Shukla, C.R.R. Chen, Development of thermal energy
II–numerical solution for isothermal boundary condition, J. Heat Transfer 106
storage materials for biomedical applications, J. Med. Eng. Technol. 39 (2015)
(1984) (1984) 469–474.
363–368.
[55] K. Sasaguchi, K. Kusano, R. Viskanta, A numerical analysis of solid-liquid phase
[22] N. Javani, I. Dincer, G.F. Naterer, G.L. Rohrauer, Modeling of passive thermal
change heat transfer around a single and two horizontal, vertically spaced
management for electric vehicle battery packs with PCM between cells, Appl.
cylinders in a rectangular cavity, Int. J. Heat Mass Transfer 40 (1997) 1343–1354.
Therm. Eng. 73 (2014) 307–316.
[56] M. Sugawara, Y. Komatsu, H. Beer, Melting and freezing around a horizontal
[23] F. Berroug, E.K. Lakhal, M. El Omari, M. Faraji, H. El Qarnia, Thermal performance
cylinder placed in a square cavity, Heat Mass Transfer 45 (2008) 83–92.
of a greenhouse with a phase change material north wall, Energy buildings 43
[57] M. Mahdaoui, T. Kousksou, S. Blancher, A. Ait Msaad, T. El Rhafiki, M. Mouqallid,
(2011) 3027–3035.
A numerical analysis of solid–liquid phase change heat transfer around a horizontal
[24] J. Stefan, Uber einige probleme der theorie der warmeleitung, Sber Akad Wiss
cylinder, Appl. Math. Model. 38 (2014) 1101–1110.
Wien 98 (1889) 473–484.
[58] K. Kant, A. Shukla, A. Sharma, P.H. Biwole, Heat transfer study of phase change
[25] J. Crank, Free and Moving Boundary Problems, Clarendon Press, Oxford, 1984.
materials with graphene nanoparticle for thermal energy storage, Sol. Energy 146
[26] J.M. Hill, One-dimensional Stefan problems: An Introduction, Longman Scientific
(2017) 453–463.
& Technical, Harlow, England, 1987.
[27] C. Bonacina, G. Comini, A. Fasano, M. Primicerio, Numerical solution of phase
change problems, Int. J. Heat Mass transfer 16 (1973) 1825–1832.

17
A. Memon et al. Applied Thermal Engineering 181 (2020) 115990

[59] J.M. Khodadadi, S.F. Hosseinizadeh, Nanoparticle-enhanced phase change [64] A.C. Kheirabadi, D. Groulx, The effect of the mushy zone constant on simulated
materials (NEPCM) with great potential for improved thermal energy storage, Int. phase change heat transfer, In: Proceedings of CHT-15, ICHMT International
Commun. Heat Mass Transfer 34 (2007) 534–543. Symposium on Advances in Computational Heat Transfer, May 25-29, 2015, NJ,
[60] P. Wang, H. Yao, Z. Lan, Z. Peng, Y. Huang, Y. Ding, Numerical investigation of USA. doi: 10.1615/ICHMT.2015.IntSympAdvComputHeatTransf.460.
PCM melting process in sleeve tube with internal fins, Energy Convers. Manage. [65] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington
110 (2016) 428–435. DC, 1980.
[61] S. Jegadheeswaran, S.D. Pohekar, Performance enhancement in latent heat thermal [66] T. Kousksou, M. Mahdaoui, M. Hlimi, R. El Alaiji, T. El Rhafiki, Latent energy
storage system: a review, Renew. Sustain. Energy Rev. 13 (2009) 2225–2244. storage: Melting process around heating cylinders, Case Stud. Therm. Eng. 8 (2016)
[62] ANSYS Inc., ANSYS Fluent user’s guide 19.2, 2019. 128–140.
[63] V.R. Voller, C. Prakash, A fixed grid numerical modelling methodology for
convection-diffusion mushy region phase-change problems, Int. J. Heat Mass
Transfer 30 (1987) 1709–1719.

18