International Communications in Heat and Mass Transfer 134 (2022) 105993

Contents lists available at ScienceDirect

International Communications in Heat and Mass Transfer

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

Effect of mushy zone constant on the melting of a solid-liquid PCM under

hyper-gravity conditions
Vivek Kumar Singh a, *, Akshat Patel b
a
Thermal Engineering Division, Space Applications Centre, ISRO, Ahmedabad, Gujarat 380015, India
b
Department of Mechanical Engineering, Institute of Technology, Nirma University, Ahmedabad, Gujarat 382481, India

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

Keywords: Mushy zone constant is an important parameter of the enthalpy-porosity method which defines how the flow
Enthalpy-porosity method evolves within the cavity. The present numerical study investigates the effect of mushy zone constant (C = 1 ×
Mushy zone constant 104, 1 × 106 & 1 × 108 kg m− 3 s− 1) on the melting phenomena of a phase change material (PCM) under hyper-
Hyper-gravity
gravity conditions (2g, 5g, 10g and 20g). Value of ‘C’ influences the fluid velocity within the liquid region of PCM
Phase-change material
and its mobility during phase change process. The melting rate is enhanced by 32% as value of ‘C’ is reduced
from 1 × 108 to 1 × 104 at normal gravitational conditions while heat transfer rate at hot wall improves as value
of ‘C’ is decreased from 1 × 108 to 1 × 104 with maximum discrepancy being 83%. Under the combined effect of
value of ‘C’ and hyper-gravity, the movement of liquid region within the cavity becomes distinct in comparison
to the normal gravity condition. The melting rate is enhanced by 74% at ‘C’ = 1 × 104 & by 54% at ‘C’ = 1 × 104
as the convection within the liquid region also enhances.

been extensively studied for aerospace and launch vehicle applications

1. Introduction since the 1970s by NASA and have also been implemented in thermal
control modules for several spaceflight missions such as Apollo 15, 16 &
Avionic systems of aerospace and launch vehicles consist of several 17 [6].
electronic subsystems which perform various important functions of the Since then, PCMs have been subjected to several experimental and
vehicle. With recent technological advancements in technology, more numerical studies to evaluate their performance under different loading
powerful electronic devices can be accommodated in smaller sizes. This conditions and configurations and also to integrate them into existing
increases the power densities of these devices and hence increases the thermal control setups of spacecraft [7–12]. However, PCMs have a poor
overall heat load on the cooling system, subsequently complicating the thermal conductivity which leads to poor heat penetration into the bulk
thermal management of such systems [1,2]. If the devices aren’t cooled of PCM. To compensate for this, PCMs are usually employed within
appropriately, their reliability and accuracy reduce subsequently metallic sinks with internal fins or through encapsulation within
reducing the overall efficiency of the system. Traditional methods for metallic sphere [13] which provides uniform heat penetration into the
thermal management of aerospace and launch vehicles include air & entirety of PCM [14–19]. Another approach to improve heat penetration
liquid cooling, pump-driven multiphase cooling, vapour cooling sys­ is to enhance PCM’s thermal conductivity by using nanoparticles. The
tems, etc. [3–5]. However, these systems have low efficiency, are heavy, addition of nanoparticles improves the thermal conductivity of PCM
and occupy large spaces within the vehicles. significantly which subsequently enhances the overall heat transfer into
Phase change materials (PCM) possess a very large latent heat stor­ the PCM [20–23]. Such strategies have been extensively used for avi­
age capacity and have drawn attention in the field of thermal control as onics thermal management of several aircraft and guided missiles
a passive method of cooling. Such storage capacity allows PCM to store [24,25].
energy about their phase transition temperature. This also damps the Aerospace and launch vehicles are exposed to a set of complicated
temperature rise of the heat-generating device. PCMs are also light­ external conditions when they take off or are launched from the ground
weight and occupy less volume within the vehicle. Thus, PCMs have level. These external conditions comprise of the vibrations induced due

Abbreviations: PCM, Phase Change Material; PISO, Pressure-Implicit with Splitting of Operators; SIMPLE, Semi-Implicit Method for Pressure Linked Equations;
PRESTO!, Pressure Staggering Option; QUICK, Quadratic Upwind Interpolation for Convective Kinetics; NITA, Non-Iterative Time Advancement.
* Corresponding author.
E-mail address: singhvivek84@gmail.com (V.K. Singh).

https://doi.org/10.1016/j.icheatmasstransfer.2022.105993

Available online 20 April 2022

0735-1933/© 2022 Elsevier Ltd. All rights reserved.
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

Nomenclature dT Differential Temperature Change (K)

Q̇ Heat transfer rate (W)
A Porosity function (kg m− 3 s− 1)
g Gravitational acceleration (m s− 2) Greek symbols
C Mushy Zone Constant (kg m− 3 s− 1) ρ Density (kg/m3)
k Thermal conductivity (W/mK) Φ Liquid Fraction
W Width (m) μ Dynamic Viscosity (kg m− 1 s− 1)
H Height (m) ε Small constant for porosity function (=0.001)
L Length (m)
Subscripts
T Temperature (K)
h Hot wall
ΔH Latent heat of fusion (J kg− 1)
c Cold wall
Cp Specific Heat (J kg− 1 K− 1)
pc Phase Change
u, v Velocity components in x & y direction (m s− 1)
i Initial
S Momentum sink term
l Lower temperature
P Fluid Pressure (Pa)
u Upper temperature
x, y Cartesian Co-ordinates
ref Reference
h Fluid Enthalpy (J kg− 1)
m mid
Pr Prandtl No
so Solid
Ra Rayleigh No
lqd Liquid
q̇ Heat flux (W m− 2)
t Time (s) Superscripts
tm Melt Time (s) Average value
n Outward normal at wall and interfaces

to drag force, external aerodynamic heating and hyper-gravity along external magnetic field. They used the Arbitrary Lagrangian-Eulerian
with other factors such as internal heat generation due to the operation method with deformed mesh to track the melt front. Such an approach
of various subsystems. The hyper-gravity condition corresponds to a required only energy balance across the melt front and didn’t need any
situation wherein an object experiences acceleration greater than the mathematical function for velocity transition. They also extended the
acceleration due to gravity. Hyper-gravity conditions exist when aero­ deformed mesh approach to investigate fluid-structure interactions for
space and launch vehicles accelerate and when aircraft and missiles flexible partitions [36,37]. Gartling [38] attempted to account for the
undertake maneuvering actions. Such varying gravity conditions may velocity transition for stationary meshes by varying the dynamic vis­
influence the performance of the thermal control systems wherein PCMs cosity of fluid according to the temperature such that its value became
have been employed and also subsequently may influence the operation zero as the temperature falls below the melting temperature. Morgan
of the overall system. Thus, the effect of such varying gravity conditions et al. [39] introduced a new method wherein the velocity was modu­
on the phase transition phenomena of PCM becomes an area of interest. lated according to the mean latent enthalpy of the cell. However, these
Filippeschi et al. [26] subjected a PCM-aluminium foam composite to methods were difficult to implement for PCMs which change their phase
hyper-gravity conditions using a large diameter centrifuge to experi­ over a range of temperature and thus weren’t completely accurate. Brent
mentally explore the influence of the same. They observed the evolution et al. [40] introduced a new parameter A(Φ) which would modify the
of the dynamic melt front through an IR camera and found that the source terms for momentum equations such that the Carman-Kozeny
melting rate was enhanced by 12% as gravity conditions were changed equations [30] for porous medium would be numerically enforced on
from 5g to 10g. Li et al. [27] employed the Lattice Boltzmann Method to the flow equations within the mushy zone. This parameter is a non-
explore the phase transition of a paraffin wax‑copper foam composite linear function of the cell porosity which in turn is computed using
under the gravity conditions 3g, 5g & 7g. They noted that the average the cell enthalpy. Thus, the parameter was such that it would be zero for
melt fraction was increased by 62% at 7g when compared to normal fully liquid cells and thus had no influence while for cells undergoing
ground conditions. Xu et al. [28] undertook a similar numerical study to phase change it would be of the same order of the other transport terms.
demonstrate that hyper-gravity condition enhances the melting rate of a For fully solid cells, the value of parameter would be sufficiently high
paraffin wax RT44 under constant heat flux conditions. such that any velocity prediction would be effectively forced to zero.
The enthalpy porosity method is the most widely used method to The parameter also contains a constant known as the mushy zone
numerically model the phase change phenomena in such studies constant ‘C’ which captures the effect of mushy zone’s morphology. The
[29,30]. The other methods to achieve the same rely on modulating the value of ‘C’ significantly influences the overall phase change process and
specific heat of the PCM about the phase transition temperature the extent of interaction of the different phases. However, value of ‘C’
[31–33]. However, the enthalpy-porosity formulation is a more general may be different for different materials. The value of ‘C’ used for
approach that can account for all types of melting phenomena while the Rubitherm RT82 range from 1 × 105 [41] to 1.6 × 106 [42] while for
latter methods are more suitable to conduction dominated phase gallium it ranges from 1 × 105 to 1 × 1015 [43,44]. Certain studies also
change. The computational region exhibiting phase transition is correlated the value of ‘C’ with the solid particle diameters within the
considered to be pseudo porous where cell porosity varies according to mushy zone constant [45–49]. Thus, the choice of mushy zone constant
the enthalpy value associated with it as enthalpy-porosity formulation is becomes very important as the behaviour predicted using the enthalpy
a stationary mesh method. However, as multiple phases are present, porosity method for the same material can be considerably different at
employing a stationary mesh technique complicates the treatment of different values of ‘C’.
velocity gradient across the stationary solid phase and mobile liquid Hong et al. [50] explored the charging and discharging of a nano­
phase. To avoid such treatment of velocity, Izadi et al. [34,35] employed particle based PCM under microgravity conditions for different values of
a deformed mesh approach based on the Galerkin finite element method ‘C’ and concluded that melting process is greatly influenced by the value
to numerically investigate the melting behaviour of PCM under an of gravity and mushy zone constant with maximum discrepancy for

2
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

liquid fraction being 88%. However, the same did not affect the dis­ thermally active walls. The overall physical model is adapted from Hong
charging process. The phase change process under hyper-gravity may et al. [50] & Arici et al. [51] so that the results of the present study can
also be affected by the type of PCM and the corresponding value of be compared directly with the results for microgravity as reported in
mushy zone constant ‘C’ which appropriately describes its behaviour. these studies.
Even though experimental studies are available in open literature for the
effect of hyper-gravity on phase change process of a PCM, the focus is on
2.2. Mathematical model
the combined performance of PCM and metal foam and effect of mushy
zone constant has been excluded.
As mentioned earlier, the enthalpy porosity method is widely used
To the best of author’s knowledge, no research paper is available in
for several numerical studies of PCM based energy systems and is also
open literature that focuses on melting characteristics of PCM under
used for the present study. The method can be implemented on a fixed
hyper-gravity conditions and at different values of mushy zone constant.
grid which allows for a reduction in computational efforts without
The present study aims to bridge this gap by numerically investigating
foregoing the level of accuracy. The mathematical model is subjected to
the phase change phenomena of a PCM at different values of mushy zone
the below assumptions to simplify the problem:
constant across varying hyper-gravity conditions and in the absence of
any metal foam or filler material which may enhance PCM’s thermal
1. PCM is isotropic and homogenous.
conductivity and consequently interfere with its heat transfer charac­
2. The liquid phase is a Newtonian liquid and its flow is laminar.
teristics. This study will highlight the effects of hypergravity on PCMs
3. Volumetric expansion of PCM during the phase change is neglected.
with different values of ‘C’ that describes their phase change behaviour.
4. The remaining dimension of the cavity is of a much larger order and
The comparison is carried out by comparing the melting and heat
its effect can be neglected.
transfer characteristics using plots of average liquid fraction and heat
transfer rate. The discussion is further supported by contour plots of
The x-y coordinate system is defined for the computational domain
average liquid fraction and velocity which highlight how the flow field
as shown in Fig. 1 and accordingly, the governing equations describing
and melt front evolves within the cavity. Comparisons are also made
the overall performance of heat storage by the PCM within the cavity can
between behaviour under hypergravity and ground gravity conditions to
be written as:
highlight and explore how different PCMs may behave when exposed to
such hypergravity conditions. This will also give a qualitative under­
Continuity Equation:
standing how different PCMs would behave when employed within
practical heat sinks for avionics of aerospace and launch vehicles. ∂ρ ∂(ρu) ∂(ρv)
+ + =0 (1)
∂t ∂x ∂y
2. Model description and solution procedure

2.1. Physical model Momentum Equations:

( ) ( )
∂(ρu) ∂(ρuu) ∂(ρuv) ∂ ∂u ∂ ∂u ∂P
To explore the melting phenomena of a PCM, a rectangular cavity is + + = μ + μ − + Sx (2a)
∂t ∂x ∂y ∂x ∂x ∂y ∂y ∂x
selected whose entire interior is assumed to be filled with PCM. The side
( ) ( )
walls of the cavity are insulated while the bottom and top wall are only ∂(ρv) ∂(ρuv) ∂(ρvv) ∂ ∂v ∂ ∂v ∂P
+ + = μ + μ − − ρg + Sy (2b)
insulated along 50% of their length. The remaining portion of the walls ∂t ∂x ∂y ∂x ∂x ∂y ∂y ∂y
is present symmetrically about the centreline and is maintained at a
constant temperature. The bottom wall (hot wall) is kept at a constant
temperature of Th = 345 K while the top wall (cold wall) is maintained at
a constant temperature of Tc = 300 K. The dimensions of the cavity are
20 mm (H) × 20 mm (W).
The PCM considered for the present study is paraffin wax whose
properties are enhanced using nanoparticles of Al2O3 and CuO. This
PCM is chosen for the present study as its properties have been properly
documented [50–53]. The PCM has a lower melting point of 319.15 K
and an upper melting point of 321.15 K. Thermophysical properties of
this PCM are listed in Table 1. The entire domain is initialized from a
temperature of Ti = 300 K and thus initial liquid fraction is zero. Thus,
for t > 0, the PCM begins to store energy as it begins interacting with the

Table 1
Thermophysical Properties of paraffin wax [50].
Property Value

Phase Change Temperature 319.15 K (Tpc,l)–321.15 K (Tpc,u)

Latent Heat of fusion (ΔH) 173,400 J kg− 1
Density (ρ) 750 if T ≤ Tpc,l
750
if T > Tpc,l
0.001(T − 319.15) + 1
Thermal Conductivity (k) 0.21 W m− 1 K− 1 if T ≤ Tpc,l
( )
0.09 T − Tpc,l
0.21 − W m− 1 K− 1 if Tpc,l < T ≤ Tpc,l
Tpc,u − Tpc,l
− 1 − 1
0.12 W m K if T > Tpc,l
Specific Heat (CP) 2890 J kg− 1 K− 1
( )
Dynamic Viscosity (μ) 1790
− 4.25 +
0.001e T kg m− 1 s− 1
Fig. 1. Schematic diagram of the PCM cavity.

3
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

Energy Equation: (2) At t > 0 for xm − Wh

≤ x ≤ xm + W2h ∩ y = 0:-
2
( ) ( )
∂(ρh) ∂(ρuh) ∂(ρvh) ∂ ∂T ∂ ∂T
+ + = k + k (3) T = Th
∂t ∂x ∂y ∂x ∂x ∂y ∂y
The source terms for momentum equations are further represented in (3) At t > 0 for xm − Wc
≤ x ≤ xm + W2c ∩ y = H:-
2
the term of the porosity function A(Φ) as below [40]:
T = Tc
Sx = − A(Φ)u (4a)
Sy = − A(Φ)v (4b)
(4) At t > 0 for the remaining walls:-
Here, Φ is the cell porosity and highlights the cell liquid fraction.
Further, the porosity function is a non-linear function of the cell porosity ∂T/∂n = 0
which is defined as below [40]:

C(1 − Φ)2 (5) For solid-liquid interface:-

A(Φ) = (5)
Φ3 + ε
∂Tso ∂Tlqd
Here, ‘C’ is the mushy zone constant that controls the velocity Tso = Tlqd ; kso = klqd
∂n ∂n.
transition across the phase interface. To avoid division by zero when cell
porosity reduces to zero (in the solid domain), a small constant ‘ε’
(=0.001) is also introduced [51,53].
Further, the liquid fraction is defined as 2.4. Solution procedure
⎧
⎪ 0 T ≤ Tpc,l
⎪
⎪
⎪ The initial and boundary conditions mentioned in Section 2.3 are
⎨ T− T
Φ=
pc,l
Tpc,l < T ≤ Tpc,u (6) now imposed on the governing equations described in Section 2.2 and
⎪
⎪
⎪
Tpc,u − Tpc,l are solved using ANSYS Fluent. Pressure based solver is used along with
⎪
⎩
1 T > Tpc,u double precision. PISO algorithm is used for pressure velocity coupling.
PISO scheme belongs to the SIMPLE family of algorithms and also uti­
Hence, the liquid fraction is zero for a cell whose temperature is lizes the neighbour and skewness correction along with velocity
below the lower melting temperature of PCM and is unity when the correction to improve the accuracy. PRESTO is used to interpolate the
temperature is above the upper melting temperature. Between these 2 pressure at faces while momentum and energy equations are discretized
temperatures i.e., for cells undergoing the phase change, the liquid using QUICK scheme. Temporal discretization is carried out using the
fraction is between zero and unity. Also, the specific enthalpy of the cell Non-Iterative Time Advancement (NITA) scheme which utilizes only a
depends on the temperature and cell liquid fraction and is defined as single iteration per time step and doesn’t require any iterations within a
follows: time step. The Non-iterative solver relaxation factors for pressure, mo­
( )
∫ T mentum and energy equations are 0.6, 0.8, and 0.8 respectively. The
h= href + CP dT + Φ(ΔH) (7) solution is said to be converged when the absolute residuals for the time
step fall below 1 × 10− 3, 1 × 10− 3 & 1 × 10− 8 for continuity, momentum
Tref

The term within the bracket corresponds to the sensible heat and energy equations respectively.
enthalpy whereas the latent heat enthalpy owing to phase change is
represented by the second term. Further, to interpret the natural con­ 2.5. Mesh and time step independence study
vection and phase change phenomena, several other data quantities are
also extracted and derived using the extracted data. The average To ensure that the results obtained are independent of the mesh size
instantaneous heat flux at the hot wall can be obtained as: and time step used, the liquid fraction variation w.r.t. flow time is
∫ Wh compared at various mesh sizes and time steps for the case of g and ‘C’ =
104. Fig. 2 shows the above-mentioned liquid fraction variation w.r.t.
xm + 2
q(t) = q̇(x, t)dx (8)
flow time at different combinations of mesh sizes and time steps.
Wh
xm − 2

Quadrilateral elements are used to generate the mesh for the computa­
Here q̇(x, t) is the instantaneous local heat flux at the hot wall, q(t) is
tional domain. The liquid fraction profile varies by a maximum of 12%
the instantaneous hot wall heat flux. Further, q(t) can be used to eval­
as mesh size is decreased from 0.5 mm to 0.1 mm and within 1% as the
uate the heat transfer rate at the hot wall using the below expression:
mesh size is decreased from 0.1 mm to 0.05 mm. Thus, a mesh size of 0.1
Q̇(t) = q(t) × Wh (9) mm is found to be considerably accurate and is selected for the present
study. Further, computational results at various time steps of 0.025 s,
Apart from the aforementioned quantities, several other quantities
0.01 s and 0.005 s are also examined and compared at the mesh size of
and variables are also exported to interpret the heat transfer and flow
0.1 mm. The variation in liquid fraction profile by a maximum of 2.15%
during the phase transition.
as the time step is reduced from 0.25 s to 0.01 s and by 0.13% as time
step is further reduced from 0.01 s and 0.005 s. Thus, a time step of 0.01
2.3. Initial and boundary conditions s is selected for the present study.

The initial and boundary conditions for the computational domain is 2.6. Validation of numerical model
as described below. Apart from the below conditions, no-slip condition is
also imposed at all walls. Phase transition phenomena of PCM without filler material like fins
and foam within hyper-gravity environment is not thoroughly investi­
(1) At t = 0 for entire domain:- gated and there is a lack of experimental data in literature for validation.
u = v = Φ = 0, T = TC Thus, the present numerical model is validated with the experimental
results of Gau & Viskanta [54] at earth gravity conditions. The experi­
mental domain consists of rectangular enclosures of different aspect

4
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

Fig. 2. Plots of Liquid fraction w.r.t flow time at various mesh size and time steps.

ratios filled with gallium metal which has a phase change temperature of change phenomena and flow within the cavity. As the mushy zone
302.78 K. The top and bottom wall of the enclosure are kept adiabatic constant, ‘C’, is an artificial constant, using different values of ‘C’ for the
while left and right wall of the enclosure are kept at constant tempera­ same PCM under the same values of Prandtl Number, Pr, and Rayleigh
tures of 311 K and 301.3 K respectively. The entire domain was Number, Ra, can give different numerical predictions. The flow and
initialized from temperature equal to that of cold wall i.e., 301.3 K. The temperature fields within the cavity and heat transfer rates at the hot
hot wall begins raising the temperature of gallium within the cavity wall, etc. can be drastically different at different values of ‘C’. Hence,
which subsequently melts and undergoes phase change. The melt before describing the combined effect of mushy zone constant and
interface formed then propagates towards cold wall. hyper-gravity conditions, it becomes vital to explore the flow patterns
The results of evolution of melt interface at different times as pre­ and the effect of ‘C’ under normal gravity conditions.
dicted by the present model is compared with the experimental results The plot of liquid fraction and heat transfer rate w.r.t. time for
for aspect ratios of 0.714, 0.5 & 0.286 as shown in Fig. 3. The mushy different values of ‘C’ is shown are Fig. 4(a) and (b) respectively.
zone constant value used for validation was C = 1.6 × 106 [40]. As Consider the average liquid fraction and heat transfer rate profile for ‘C’
evident from Fig. 3, it is evident that the present numerical model can = 1 × 104 (highlighted in blue) in Fig. 4. As the constant temperature
predict the melt front evolution and the phase change process of gallium boundary conditions at hot and cold wall are imposed at t = 0 s, the PCM
metal with considerable accuracy. The maximum error in the position of begins storing heat in the vicinity of hot wall. When the PCM layers
melt interface is approximately 6.4% (for aspect ratio of 0.5 at 17 mins). adjacent to hot wall reaches the lower phase change temperature (Tpc,l),
Hence the present numerical model is found to be considerably accurate they begin to exhibit phase transition and slowly form a thin liquid layer
to model phase change process and is further used for the present study. in the vicinity of the same. Point 1 in Fig. 4 corresponds to this state.
However, this liquid layer remains stagnant till sufficient amount of
3. Results and discussion PCM has undergone phase change so that sufficient space is available for
the convection current to develop.
The present study aims to understand the effect of mushy zone This occurs at point 2 and this point also consequently marks the
constant on the phase change process under hyper-gravity conditions. onset of natural convection within the cavity. Upto point 2, predominant
Hence to carry out the same, the numerical model is solved at mushy mode of heat transport is diffusion and thus similar amount of PCM has
zone constant values of 104, 106 & 108 and gravitational values of 2g, 5g, undergone phase transition for all values of ‘C’. This is highlighted in the
10g and 20g. Further, the mushy zone constant values selected represent plot of liquid fraction where in the plots coincide with each other for all
the recommended range of values of ‘C’ [55]. values of ‘C’. Correspondingly, heat transfer rate decreases as the system
Section 3.1 first describes and establishes the effect of the mushy progresses to point 2 from point 1. This is because the hot wall now
zone constant under normal gravitational value and how it affects the dissipates heat energy via diffusion into the liquid layer which has
phase change process within the cavity. Section 3.2 then extends the smaller value of thermal conductivity. After point 2, natural convection
discussion to hyper-gravity conditions and describes the combined effect begins to develop within the liquid region and plots of liquid fraction
of mushy zone constant and hyper-gravity conditions. become distinct implying that the different values of ‘C’ have different
effects on the convection intensity within the fluid layer. The heat
transfer rate also begins to increase after experiencing a local-minima at
3.1. Effect of mushy zone constant under normal gravity condition (g)
point 2. The plots liquid fraction begins increasing at an increasing rate
as the natural convection within the fluid layer becomes progressively
This section describes how the mushy zone constant affects the phase

5
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

Fig. 3. Comparison of melt front of Solid-Liquid interface at different times for present numerical model and experimental results of Gau & Viskanta [54] at aspect
ratios of 0.714 (top), 0.5 (middle) and 0.286 (bottom).

more intense. Correspondingly, the liquid layer begins to grow in both 3. Point 3 denotes the state when liquid region has completely filled up
directions (predominantly in +y direction due to buoyancy effect the bottom of the cavity. The heat transfer rate also reaches its maximum
induced by the gravitational force). value at this point owing to a mixed effect of heat storage in latent form
The liquid fraction keeps on increasing at an increasing rate till point and natural convection in the bottom region. Beyond point 3, liquid

6
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

Fig. 4. (a) Liquid fraction variation w.r.t time for ‘C’ = 1 × 104, 1 × 106 & 1 × 108 (kg m− 3 s− 1) (b) Heat transfer rate as a function of time at hot wall Melt fraction
contours at (1) t = 6 s (2) t = 60 s (3) t = 730 s (4) t = 1200 s (5) t = 1590 s (Liquid and solid regions are indicated by red and blue respectively). (For interpretation of
the references to colour in this figure legend, the reader is referred to the web version of this article.)

fraction profile increases at a decreasing rate while the heat transfer rate Convection current within the liquid region now begins penetrating
starts to decrease. This is because further phase change occurs away the upper stagnant solid PCM. The liquid region is able to penetrate the
from the hot wall and the predominant mode of heat transfer in the solid PCM only upto a certain height in the vicinity of the cold wall. This
bottom region is now only natural convection. state corresponds to point 4 on the liquid fraction and heat transfer rate

7
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

profile. Beyond point 4, the liquid region begins penetrating the top Similar behaviour is also observed at different values of ‘C’. The plots
corners of the cavity and only expands predominantly in the x direction. liquid fraction w.r.t time for ‘C’ = 1 × 104, 1 × 106 & 1 × 108 are shown
Also, the heat transfer rate increases momentarily and reaches another in Fig. 4(a). The melting rates are maximum for ‘C’ = 1 × 104 and are the
local-maxima at point 4 as the liquid region also begins dissipating heat least for ‘C’ = 1 × 108. This implies that the smaller values of ‘C’ result in
into the remaining solid PCM layer in the vicinity of cold wall which is greater melting rates under a similar convection setup. The melting rate
thin enough to conduct away this heat into the cold wall. reduces by a maximum of 23.6% as the ‘C’ value changes from 1 × 104 to
The liquid fraction profile further increases at a decreasing rate and 1 × 106 and reduces by a maximum of 10.53% as the value of ‘C’
ultimately reaches its steady-state value (~0.95) at point 5 as the PCM in changes from 1 × 106 to 1 × 108. The corresponding heat transfer rate is
the top corners also begins changing phase ultimately resulting in a also maximum for lower values of ‘C’ which is depicted in Fig. 4(b).
liquid fraction contour corresponding to point 5. A small amount of solid The velocity (top) and stream function (bottom) contours at flow
PCM always exists near the cold wall and as the temperature of the cold time when the average liquid fraction is 0.25,0.5 and 0.75 for ‘C’ = 1 ×
wall is well below the lower phase change temperature, this solid PCM 104, 1 × 106 & 1 × 108 is shown in Fig. 5. The velocity contours indicate
never exhibits phase transition. that several counter-rotating convection cells exist within the liquid

Fig. 5. Velocity (top) and stream function (bottom) contours at tm/4, tm/2 & 3tm/4 at ‘C’ = 1 × 104, 1 × 106 & 1 × 108.

8
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

region that tends to further penetrate upwards into solid PCM owing to gravity conditions for different values of ‘C’. The melt time reduces as
density changes at ‘C’ = 1 × 104 whereas for ‘C’ = 1 × 108, only a single the value of gravity is increased because the convection intensity within
central pair of such counter-rotating convection cell exists. Further, the the cavity also increases. The decrease in melt time from g to 20g is
contours of stream function highlight how these convection cells within 73.5% for ‘C’ = 1 × 104, 62.72% for ‘C’ = 1 × 106 and 54.46% for ‘C’ =
the liquid region tend to rotate and further propagate into the solid PCM 1 × 108. Further, the change in melt time also decreases by different
layer with time. As several convection cells exist for ‘C’ = 1 × 104, magnitudes at different hyper-gravity values. The change is greater at a
several such convection currents also exist within the cavity as smaller value of gravity, while it is smaller at a greater value of gravity.
compared to central convection current like in the case of ‘C’ = 1 × 108. This is also observed in the plots of liquid fraction w.r.t. time under
This results in greater convection intensity and mobility of the fluid different values of g for different values of ‘C’ as shown in Fig. 7. The
region near the hot wall and thus the heat transfer rates for ‘C’ = 1 × 104 are vertical dimension of the cavity is finite and hence even though buoy­
also superior as compared to that of ‘C’ = 1 × 108. Such multiple convection ancy force increases with the gravity value, the flow can only rise upto a
cells also exist for ‘C’ = 1 × 108 at very early stages but they tend to merge certain height within the cavity after which it starts to expand in the x-
into the central convection cell at a faster rate. The presence of multiple direction. Thus, the flow reaches the cold wall at similar times across
convection cells and convection currents allows for greater mobility of the different gravity values and the corresponding reduction in melt time is
overall melt front and thus improves the natural convection within the also smaller. Thus, the change in melt time decreases as gravitational
cavity allowing for better heat transfer rates. Hence with an increase in acceleration increases.
value of ‘C’, the time period for which multiple convection currents and The variation of the liquid fraction at higher values of gravitational
counter-rotating convection cells exist decreases. This in turn weakens the acceleration (5g, 10g & 20g) is different than the variation at normal
natural convection and subsequently reduces the heat transfer. gravity conditions as described in Section 3.1. The liquid fraction in­
creases at a linear rate at higher gravity values for ‘C’ = 1 × 104 while it
first increases at an increasing rate and then at a decreasing rate at lower
3.2. Effect of hyper-gravity gravity values for ‘C’ = 1 × 106 & 1 × 108. As the hyper-gravity value is
increased for smaller values of ‘C’, the plot of liquid fraction tends to
The gravitational force is responsible for inducing the natural con­ become linear indicated as indicated by the plots in Fig. 7. Such linear
vection within the cavity. As the fluid in the vicinity of the hot wall rate of increase is also observed at ‘C’ = 1 × 106 & 1 × 108 but only at
becomes hotter, it also becomes less dense and begins to expand and greater values of g.
experiences greater buoyancy force. Thus, the hotter fluid rises upwards The liquid fraction and velocity contours across different gravity
against the gravity force and is replaced by the colder fluid which now values for ‘C’ = 1 × 104 is shown in Fig. 8 at the instances when the
interacts with the hot wall. This colder fluid now takes away the heat average liquid fraction is 0.25, 0.5, and 0.75. As the value of gravity
from the hot wall and subsequently becomes hotter and rise upwards increases, a greater buoyancy force is exerted over the fluid and hence
and the entire cycle continues. Thus, under varying gravity conditions, the tendency of the fluid to rise upwards is also greater. Thus, the flow
the tendency of the hotter fluid to rise upwards will also differ leading to now expands in the vertical direction by a greater magnitude and the
different convection patterns and heat transfer characteristics. Further, tendency of flow to expand in the x-direction is smaller. Hence, for the
as PCM also undergoes phase change, the convection within the cavity same level of liquid fraction, the liquid region is now elongated along
will also be affected by the value of mushy zone constant along with the +y direction and becomes narrower for 5g,10g, and 20g. The liquid
hyper-gravity value. Hence, the same is investigated under hyper- fraction contours for 10g & 20g also highlight how the liquid region
gravity values of 2g, 5g, 10g, and 20g. penetrates the solid PCM in a manner different than that of the normal
Fig. 6 shows the comparison of time taken for the average liquid gravity conditions. The liquid region first reaches the cold wall & top
fraction to reach its steady-state (melt time) value under varying hyper-

Fig. 6. Melt time comparison under varying hyper-gravity conditions for ‘C’ = 1 × 104, 1 × 106 & 1 × 108.

9
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

Fig. 7. Plots of liquid fraction and heat transfer rate w.r.t time for ‘C’ = 1 × 104, 1 × 106 & 1 × 108.

corners and then expands in the x-direction at 10g & 20g owing to a Fig. 9 shows the liquid fraction and velocity contours across different
greater magnitude of buoyancy force along the vertical direction. The gravity values for ‘C’ = 1 × 104, 1 × 106 & 1 × 108 at the instances when
liquid layer at normal gravity conditions first penetrates the bottom the average liquid fraction is 0.5. As the gravity value is increased, the
corners of the cavity and then reaches the top corners and cold wall as liquid region also rises upto greater heights under the effect of increased
described in Section 3.1. Thus, the plot of liquid fraction is steeper at 5g, buoyancy force. However, the liquid region for ‘C’ = 1 × 104 rise upto
10g, and 20g as the phase change occurs in the bottom region for greater greater heights as compared to liquid region for ‘C’ = 1 × 106 & 1 × 108
durations. The abovementioned behaviour is observed at ‘C’ = 1 × 104 even for the similar values of liquid fraction. The liquid fraction contour
kg m− 3 s− 1. However, the behaviour at ‘C’ = 1 × 106 & 1 × 108 is of ‘C’ = 1 × 104 at greater values of g also has uniform width throughout
different than the behaviour at ‘C’ = 1 × 104. its height whereas the liquid fraction contours of ‘C’ = 1 × 106 & 1 × 108

10
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

Fig. 8. Liquid Fraction (left half) and Velocity (right half) contours at tm/4, tm/2 & 3tm/4 & g, 2g, 5g, 10g and 20g for ‘C’ = 1 × 104.

have a diamond-like shape. This is because, at smaller values of ‘C’, the ‘C’ and g is shown in the right column of Fig. 7. With an increase in the
magnitude of velocities within the mushy zone obtained through the gravitational acceleration, the intensity of natural convection increases,
numerical solution is greater as compared to the magnitude of velocities and heat transfer also improves.
obtained for larger values of ‘C’. Thus, the melt front has greater Heat transfer rate initially increases as the convection within the
mobility at smaller values of ‘C’ and the melt front can expand in the +y liquid region becomes progressively stronger and then decreases as the
direction for smaller values of ‘C’ and thus the liquid region rises to a melt front rises upto its maximum height and begins to expand in the x
greater height. But as the value of ‘C’ increases, the mobility of the melt direction. Once the liquid fraction reaches its steady state value, the heat
front also reduces, and thus the tendency of the liquid region to rise also transfer rate tend to fluctuate about a mean value. This is because the
decreases. The melt front rather also expands along the x direction thus flow instability due to the presence of thermally unstable layers (hotter
forming diamond like shape. Thus, melt front at greater values of ‘C’ and and sparser hot layer below the colder denser layer). Hence, the heat
g eventually reaches the bottom corners of the cavity first and further transfer rate is shown by the linear trendline after the average liquid
penetrates the solid PCM in the same way as that of normal gravity fraction reaches its steady state.
conditions. The variation of heat transfer with time for various values of The heat transfer rate is maximum for ‘C’ = 1 × 104 and 20g, as this

11
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

Fig. 9. Liquid Fraction (left half) and Velocity (right half) contours at tm/2 & g, 2g, 5g, 10g and 20g for ‘C’ = 1 × 104, 1 × 106 & 1 × 108.

combination results in maximum mobility of melt front and fastest function of cell liquid fraction. This porosity function in turn modu­
melting rate due to stronger convection within the liquid region. Also, lates source term of flow equations to account for velocity gradient
the due to high mobility of melt front, the liquid region has greater across the phase interface. To carry out the numerical investigation, a
tendency to rise upwards and hence the phase change process occurs in physical model consisting of paraffin wax filled in a square cavity with
the bottom of the cavity for longer durations. Thus, the heat dissipation partially thermally active walls is considered. Different values of ‘C’
within the cavity near the hot wall is a combination of heat storage in under a similar numerical setup can yield different fluid flow and heat
latent form and convection resulting in very high melting rates and transfer characteristics. Furthermore, values of ‘C’ for various PCMs can
better heat transfer rates. be different. Hence, the appropriate choice of ‘C’, for any future nu­
merical study becomes very important.
4. Conclusion The major conclusions are as follows:

This study explores the influence of mushy zone constant, ‘C’, on the • The value of ‘C’ influences the melt front mobility and its shape.
phase change phenomena under hyper-gravity conditions. The value of Lower values of ‘C’, result in a greater magnitude of velocities within
‘C’ modulates the value of porosity function which is a non-linear the mushy zone leading to better mobility of the melt front. Thus, the

12
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

liquid region can expand and penetrate the solid PCM at a quicker Declaration of Competing Interest
rate. This results in a faster onset of convection and greater con­
vection intensity within the liquid region. Subsequently, the melting None.
rate is maximum for ‘C’ = 1 × 104 and reduces by ~32% as value of
‘C’ is increased to 1 × 108. The corresponding heat transfer rates are References
maximum for ‘C’ = 1 × 104 (with the maximum heat transfer rate
being ~110 W) and reduces by 83% as the value of ‘C’ is increased to [1] M. Dooley, C. Lui, R. Newman, Efficient propulsion, power, and thermal
management integration, in: 49th AIAA/ASME/SAE/ASEE Jt. Propuls. Conf. 1
‘C’ = 1 × 108. PartF, 2013, pp. 1–8, https://doi.org/10.2514/6.2013-3681.
• The shape of the melt front has greater curvature at lower values of [2] A. Akin, H.S. Kahveci, An optimization study for rotorcraft avionics bay cooling,
‘C’ as multiple convection cells apart from the central cell are pre­ Aerosp. Sci. Technol. 90 (2019) 1–11, https://doi.org/10.1016/j.ast.2019.04.029.
[3] T. Dixit, I. Ghosh, Review of micro- and mini-channel heat sinks and heat
sent. These convection cells eventually merge into the central con­ exchangers for single phase fluids, Renew. Sust. Energ. Rev. 41 (2015) 1298–1311,
vection cell but at much later stages (~900 s for ‘C’ = 1 × 104 and g) https://doi.org/10.1016/J.RSER.2014.09.024.
and tend to exist for longer durations at lower values of ‘C’. Melt [4] Z. He, Y. Yan, Z. Zhang, Thermal management and temperature uniformity
enhancement of electronic devices by micro heat sinks: a review, Energy 216
front mobility is greater at lower values of ‘C’ and the liquid region (2021), 119223, https://doi.org/10.1016/J.ENERGY.2020.119223.
penetrates the solid PCM predominantly in the vertical direction. But [5] I. Mudawar, Assessment of high-heat-flux thermal management schemes, IEEE
the mobility reduces with an increase in value of ‘C’ and the liquid Trans. Components Packag. Technol. 24 (2001) 122–141, https://doi.org/
10.1109/6144.926375.
region penetrates the solid PCM in the vertical direction by a lesser
[6] J.P. Collette, P. Rochus, R. Peyrou-Lauga, O. Pin, N. Nutal, M. Larnicol, J. Crahay,
degree. Subsequently, the melt front is more elongated the vertical Phase change material device for spacecraft thermal control, in: 62nd Int.
direction and has smaller width at a smaller value of ‘C’ while its Astronaut. Congr. 2011, IAC 2011 7, 2011, pp. 6020–6031.
shape it diamond-like at higher values. [7] D.G. Gilmore, Spacecraft Thermal Control Handbook Volume I, 2002. http://app.
knovel.com/hotlink/toc/id:kpSTCHVFT2/spacecraft-thermal-control.
• An increase in the value of gravity increases the buoyancy force [8] A. Abhat, M. Groll, Investigation of Phase Change Material /PCM/ Devices for
being exerted on the fluid. This results in a greater tendency of the Thermal Control Purposes in Satellites, 1974, https://doi.org/10.2514/6.1974-
hotter fluid to rise upwards which also improves the convection 728.
[9] T.Y. Kim, B.S. Hyun, J.J. Lee, J. Rhee, Numerical study of the spacecraft thermal
within the liquid region. With an increase in the gravitational ac­ control hardware combining solid-liquid phase change material and a heat pipe,
celeration, PCM’s melting rate enhances across the different values Aerosp. Sci. Technol. 27 (2013) 10–16, https://doi.org/10.1016/j.
of ‘C’. However, the nature of increase is different for different values ast.2012.05.007.
[10] G. Quinn, E. Hodgson, R. Stephan, Phase Change Material Trade Study : A
of ‘C’. The melting rate is approximately linear for ‘C’ = 1 × 104 Comparison Between Wax and Water for Manned Spacecraft, 2019, pp. 1–13.
under hyper-gravity conditions of 5g, 10g and 20g while for larger [11] K. Yamada, H. Nagano, Y. Kobayashi, T. Totani, H.S. Panel, P.C. Material, C. Fiber,
values of ‘C’, the rate first increases and then decreases. This is R. Polymer, D.S. Calorimetry, Encapsulated in a High - thermal conductivity CFRP
for, 2014.
because owing to greater mobility at ‘C’ = 1 × 104, the liquid region [12] H. Dong, C. Yang, B. Gong, K. Zi, Thermal Control of CCD Focal Plane Assembly for
penetrates the upper portion of the cavity first while at ‘C’ = 1 × 106 Attitude-varied Space Camera, 2010 Int. Conf. Comput. Mechatronics, Control
& ‘C’ = 1 × 108, it fills up the bottom half first. Thus, in the former, Electron. Eng. C 2, 2010, pp. 597–599, https://doi.org/10.1109/
CMCE.2010.5610400.
phase change occurs in the vicinity of hot wall for longer duration,
[13] Q. Xiong, H.M. Alshehri, R. Monfaredi, T. Tayebi, F. Majdoub, A. Hajjar,
while in the latter, phase change occurs away from the hot wall after M. Delpisheh, M. Izadi, Application of phase change material in improving trombe
a certain time. wall efficiency: an up-to-date and comprehensive overview, Energy Build. 258
• Melting rate increases with increase in the gravitational acceleration (2022), 111824, https://doi.org/10.1016/j.enbuild.2021.111824.
[14] M. Izadi, M. Sheremet, H.M. Alshehri, T. Ambreen, M.H. Doranehgard, Numerical
by a greater magnitude at lower values of g (g, 2g & 5g) than the study on charging process inside a grid-structure thermal storage, J. Energy
increase at higher values (10g & 20g). This is because even though Storage 45 (2022), 103522, https://doi.org/10.1016/j.est.2021.103522.
the buoyancy force increase with the value of gravity, the liquid [15] H.M. Ali, A. Arshad, Experimental investigation of n-eicosane based circular pin-fin
heat sinks for passive cooling of electronic devices, Int. J. Heat Mass Transf. 112
region can rise upto a certain height only owing to finite vertical (2017) 649–661, https://doi.org/10.1016/j.ijheatmasstransfer.2017.05.004.
dimensions of the cavity. [16] S.A. Shehzad, B. Alshuraiaan, M.S. Kamel, M. Izadi, T. Ambreen, Influence of fin
• Heat transfer rates improve with an increase in the value of g as the orientation on the natural convection of aqueous-based nano-encapsulated PCMs in
a heat exchanger equipped with wing-like fins, Chem. Eng. Process. 160 (2021),
melt front expands at a quicker rate allowing for better convection 108287, https://doi.org/10.1016/j.cep.2020.108287.
within the liquid region thus resulting in better heat transfer rates. [17] C.R. Raj, S. Suresh, R.R. Bhavsar, V.K. Singh, K.A. Govind, Influence of fin
Additionally, heat transfer rates for ‘C’ = 1 × 104 are maximum as configurations in the heat transfer effectiveness of solid solid PCM based thermal
control module for satellite avionics: numerical simulations, J. Energy Storage 29
the heat transfer in the vicinity of the hot wall due to a combined (2020), https://doi.org/10.1016/j.est.2020.101332.
effect of latent heat storage and natural convection for longer [18] A.N. Desai, H. Shah, V.K. Singh, Novel inverted fin configurations for enhancing
durations. the thermal performance of PCM based thermal control unit: a numerical study,
Appl. Therm. Eng. 195 (2021), 117155, https://doi.org/10.1016/J.
• An increase in the melting and heat transfer rate also imply that PCM
APPLTHERMALENG.2021.117155.
based heat sinks would perform better under hyper-gravity condi­ [19] A.N. Desai, A. Gunjal, V.K. Singh, Numerical investigations of fin efficacy for phase
tions given that sufficient amount of PCM is available to compensate change material (PCM) based thermal control module, Int. J. Heat Mass Transf. 147
for the higher melting rate. Subsequently, the performance of PCMs (2020), 118855, https://doi.org/10.1016/J.
IJHEATMASSTRANSFER.2019.118855.
with a smaller value of ‘C’ such as Rubitherm RT42 (‘C’ = 1 × 105 [20] C.R. Raj, S. Suresh, V.K. Singh, R.R. Bhavsar, M. Chandrasekar, V. Archita, Life
[50]) would deviate by a greater magnitude under hyper-gravity cycle assessment of nanoalloy enhanced layered perovskite solid-solid phase
conditions than PCMs with a greater value of ‘C’ such as Paraffin change material till 10000 thermal cycles for energy storage applications,
J. Energy Storage 35 (2021), 102220, https://doi.org/10.1016/j.est.2020.102220.
wax (‘C’ = 1 × 108 [50]) when compared to their respective earth [21] C.R. Raj, S. Suresh, V.K. Singh, R.R. Bhavsar, M.V. C, S. Vasudevan, V. Archita,
gravity performances. Experimental investigation on nanoalloy enhanced layered perovskite PCM tamped
in a tapered triangular heat sink for satellite avionics thermal management, Int. J.
Therm. Sci. 167 (2021), 107007, https://doi.org/10.1016/J.
CRediT authorship contribution statement IJTHERMALSCI.2021.107007.
[22] Q. Xiong, A. Hajjar, B. Alshuraiaan, M. Izadi, S. Altnji, S.A. Shehzad, State-of-the-
Vivek Kumar Singh: Conceptualization, Supervision, Methodology, art review of nanofluids in solar collectors: a review based on the type of the
dispersed nanoparticles, J. Clean. Prod. 310 (2021), 127528, https://doi.org/
Investigation, Writing – review & editing. Akshat Patel: Methodology, 10.1016/j.jclepro.2021.127528.
Software, Validation, Formal analysis, Investigation, Visualization, [23] Q. Xiong, S. Altnji, T. Tayebi, M. Izadi, A. Hajjar, B. Sundén, L.K.B. Li,
Writing – original draft. A comprehensive review on the application of hybrid nanofluids in solar energy
collectors, Sustain. Energy Technol. Assess. 47 (2021), https://doi.org/10.1016/j.
seta.2021.101341.

13
V.K. Singh and A. Patel International Communications in Heat and Mass Transfer 134 (2022) 105993

[24] D. Pal, Y.K. Joshi, Thermal management of an avionics module using solid-liquid [40] A.D. Brent, V.R. Voller, K.J. Reid, Enthalpy-porosity technique for modeling
phase-change materials, J. Thermophys. Heat Transf. 12 (1998) 256–262, https:// convection-diffusion phase change: application to the melting of a pure metal,
doi.org/10.2514/2.6329. Numer. Heat Transf. 13 (1988) 297–318, https://doi.org/10.1080/
[25] A.J. Fossett, M.T. Maguire, A.A. Kudirka, F.E. Mills, D.A. Brown, Avionics passive 10407788808913615.
cooling with microencapsulated phase change materials, J. Electron. Packag. [41] A.M. Abdulateef, J. Abdulateef, S. Mat, K. Sopian, B. Elhub, M.A. Mussa,
Trans. ASME 120 (1998) 238–242, https://doi.org/10.1115/1.2792628. Experimental and numerical study of solidifying phase-change material in a
[26] S. Filippeschi, M. Mameli, P. Di Marco, Experimental Analysis of the Melting triplex-tube heat exchanger with longitudinal/triangular fins, Int. Commun. Heat
Process in a PCM / Aluminum Foam Composite Material in Hypergravity Mass Transf. 90 (2018) 73–84, https://doi.org/10.1016/j.
Conditions 6, 2019, pp. 451–467. icheatmasstransfer.2017.10.003.
[27] X. Li, Z. Zhu, Z. Xu, T. Ma, H. Zhang, J. Liu, Q. Wang, Effect of supergravity on heat [42] K. Nithyanandam, R. Pitchumani, Computational modeling of dynamic response of
transfer characteristics of PCM with the pore-scale lattice Boltzmann method, a latent thermal energy storage system with embedded heat pipes, J. Sol. Energy
Energy Proc. 158 (2019) 4641–4647, https://doi.org/10.1016/j. Eng. Trans. ASME 136 (2014) 1–9, https://doi.org/10.1115/1.4024745.
egypro.2019.01.742. [43] X.H. Yang, S.C. Tan, J. Liu, Numerical investigation of the phase change process of
[28] Y. Xu, Z. Zhu, S. Li, J. Wang, Numerical investigation on melting process of a phase low melting point metal, Int. J. Heat Mass Transf. 100 (2016) 899–907, https://
change material under supergravity, J. Therm. Sci. Eng. Appl. 13 (2021). doi.org/10.1016/j.ijheatmasstransfer.2016.04.109.
[29] V.R. Voller, M. Cross, N.C. Markatos, An enthalpy method for convection/diffusion [44] T. Kousksou, M. Mahdaoui, A. Ahmed, A.A. Msaad, Melting over a wavy surface in
phase change, Int. J. Numer. Methods Eng. 24 (1987) 271–284, https://doi.org/ a rectangular cavity heated from below, Energy. 64 (2014) 212–219, https://doi.
10.1002/nme.1620240119. org/10.1016/j.energy.2013.11.033.
[30] V.R. Voller, C. Prakash, A fixed grid numerical modelling methodology for [45] M. Mbaye, E. Bilgen, Phase change process by natural convection-diffusion in
convection diffusion mushy region phase change problems, Int. J. Heat Mass rectangular enclosures, Heat Mass Transf. Und Stoffuebertragung 37 (2001) 35–42,
Transf. 30 (1978) 1709–1719. https://doi.org/10.1007/s002310000095.
[31] H. Hu, S.A. Argyropoulos, Mathematical modelling of solidification and melting: a [46] O. Ben-David, A. Levy, B. Mikhailovich, A. Azulay, 3D numerical and experimental
review, Model. Simul. Mater. Sci. Eng. 4 (1996) 371–396, https://doi.org/ study of gallium melting in a rectangular container, Int. J. Heat Mass Transf. 67
10.1088/0965-0393/4/4/004. (2013) 260–271, https://doi.org/10.1016/j.ijheatmasstransfer.2013.07.058.
[32] J.H.N. Ehms, R. De Césaro Oliveski, L.A.O. Rocha, C. Biserni, M. Garai, Fixed grid [47] T. Tabassum, M. Hasan, L. Begum, 2-D numerical investigation of melting of an
numerical models for solidification and melting of phase change materials (PCMs), impure PCM in the arbitrary-shaped annuli, Int. J. Therm. Sci. 114 (2017)
Appl. Sci. 9 (2019), https://doi.org/10.3390/app9204334. 296–319, https://doi.org/10.1016/J.IJTHERMALSCI.2017.01.006.
[33] V.R. Voller, C.R. Swaminathan, B.G. Thomas, Fixed grid techniques for phase [48] C. Beckermann, R. Viskanta, Natural convection solid/liquid phase change in
change problems: a review, Int. J. Numer. Methods Eng. 30 (1990) 875–898, porous media, Int. J. Heat Mass Transf. 31 (1988) 35–46, https://doi.org/10.1016/
https://doi.org/10.1002/nme.1620300419. 0017-9310(88)90220-7.
[34] S.A.M. Mehryan, A. Tahmasebi, M. Izadi, M. Ghalambaz, Melting behavior of phase [49] W. Libeer, F. Ramos, C. Newton, M. Alipanahrostami, C. Depcik, X. Li, Two-phase
change materials in the presence of a non-uniform magnetic-field due to two heat and mass transfer of phase change materials in thermal management systems,
variable magnetic sources, Int. J. Heat Mass Transf. 149 (2020), 119184, https:// Int. J. Heat Mass Transf. 100 (2016) 215–223, https://doi.org/10.1016/J.
doi.org/10.1016/j.ijheatmasstransfer.2019.119184. IJHEATMASSTRANSFER.2016.04.076.
[35] M. Izadi, M. Ghalambaz, S.A.M. Mehryan, Location impact of a pair of magnetic [50] Y. Hong, W.B. Ye, J. Du, S.M. Huang, Solid-liquid phase-change thermal storage
sources on melting of a magneto-Ferro phase change substance, Chin. J. Phys. 65 and release behaviors in a rectangular cavity under the impacts of mushy region
(2020) 377–388, https://doi.org/10.1016/j.cjph.2020.03.002. and low gravity, Int. J. Heat Mass Transf. 130 (2019) 1120–1132, https://doi.org/
[36] A.B. Shahrestani, B. Alshuraiaan, M. Izadi, Combined natural convection-FSI inside 10.1016/j.ijheatmasstransfer.2018.11.024.
a circular enclosure divided by a movable barrier, Int. Commun. Heat Mass Transf. [51] M. Arıcı, E. Tütüncü, M. Kan, H. Karabay, Melting of nanoparticle-enhanced
126 (2021), 105426, https://doi.org/10.1016/j.icheatmasstransfer.2021.105426. paraffin wax in a rectangular enclosure with partially active walls, Int. J. Heat Mass
[37] S.A.M. Mehryan, M. Ghalambaz, R. Kalantar Feeoj, A. Hajjar, M. Izadi, Free Transf. 104 (2017) 7–17, https://doi.org/10.1016/j.
convection in a trapezoidal enclosure divided by a flexible partition, Int. J. Heat ijheatmasstransfer.2016.08.017.
Mass Transf. 149 (2020), https://doi.org/10.1016/j. [52] W.R. Humphries, E.I. Griggs, A Design Handbook for Phase Change Thermal
ijheatmasstransfer.2019.119186. Control and Energy Storage Devices, Nasa, 1977, p. 255.
[38] D.K. Gartling, Finite element analysis of convective heat transfer problems with [53] H. Shmueli, G. Ziskind, R. Letan, Melting in a vertical cylindrical tube: numerical
change of phase, in: Conf. Numer. Methods Laminar Turbul. Flow, Swansea, Uk, investigation and comparison with experiments, Int. J. Heat Mass Transf. 53 (2010)
1978. 4082–4091, https://doi.org/10.1016/J.IJHEATMASSTRANSFER.2010.05.028.
[39] K. Morgan, A numerical analysis of freezing and melting with convection, Comput. [54] C. Gau, R. Viskanta, Melting and solidification of a pure metal on a vertical wall,
Methods Appl. Mech. Eng. 28 (1981) 275–284, https://doi.org/10.1016/0045- J. Heat Transf. 108 (1986) 174–181, https://doi.org/10.1115/1.3246884.
7825(81)90002-5. [55] Fluent 15.0 User Guide, Fluent Inc, USA, 2013.

14