Received: 25 September 2018 Revised: 22 November 2018 Accepted: 25 November 2018

DOI: 10.1002/er.4363

SPECIAL ISSUE PAPER

Analysis of solidification in a finite PCM storage with

internal fins by employing heat balance integral method

Rohit Kothari1 | Sreetam Das2 | Santosh K. Sahu1 | Shailesh I. Kundalwal1

1
Discipline of Mechanical Engineering,
Summary
Indian Institute of Technology Indore,
Indore, India Here, a simplified analytical model has been proposed to predict solid fraction,
2
Department of Mechanical Engineering, solid–liquid interface, solidification time, and temperature distribution during
College of Engineering and Technology, solidification of phase change material (PCM) in a two‐dimensional latent heat
Bhubaneswar, Bhubaneswar, India
thermal energy storage system (LHTES) with horizontal internal plate fins.
Correspondence Host of boundary conditions such as imposed constant heat flux, end‐wall tem-
Rohit Kothari, Discipline of Mechanical
perature, and convective air environment on the vertical walls are considered
Engineering, Indian Institute of
Technology Indore, Khandwa Road, for the analysis. Heat balance integral method was used to obtain the solution.
Simrol, Indore 453552, India. Present model yields closed‐form solution for temperature variation and solid
Email: [email protected]
fraction as a function of various modeling parameters. Also, solidification time
Funding information of PCM, which is useful in optimum design of PCM‐based thermal energy stor-
Department of Science and Technology, ages, has been evaluated during the analysis. The solidification time was found
Government of India, Grant/Award
Number: DST/TMD/MES/2k17/65;
to be reduced by 93% by reducing the aspect ratio from 8 to 0.125 for constant
DST‐INSPIRE Fellowship program, heat flux boundary condition. While, for constant wall temperature boundary
Government of India, Grant/Award
condition, the solidification time reduces by 99% by changing the aspect ratio
Number: IF170534
from 5 to 0.05. In case of convective air boundary surrounding, the solidifica-
tion time is found to reduce by 88% by reducing the aspect ratio from 8 to
0.125. Based on the analytical solution, correlations have been proposed to pre-
dict solidification time in terms of aspect ratio and end‐wall boundary
condition.

KEYWORDS
constant temperature phase transition, interface of solid and liquid, latent heat thermal energy
storage (LHTES), phase change material (PCM), solidification time

Nomenclature:
English Symbols: Ts, Temperature of solid PCM, °C; T f , Temperature of fin, °C; Lp, PCM's heat of fusion, J/kg; Tm, Solidification temperature of
PCM, °C; cs, Solid PCM's specific heat, J/kg/K; lc, PCM cell's half height, m; q″w , Wall heat flux, W/m2; t, Time, s; X(t), Interface of solid and liquid
(direction X), m; Y(t), Interface of solid and liquid (direction Y), m; St , St , cs ΔT ; cl ΔT (Stefan's number); k, Thermal conductivity of PCM, W/m°C;
s l Lp Lp
h h
k f , Thermal conductivity of fin material, W/m C; h, Coefficient of heat transfer during convection, W/m2/K; l f , Length of fin, m; Hs,
°
; H f , ; x,
ks kf
y, r, Coordinates
Greek Symbols: αs, αl, Solid and liquid PCM's thermal diffusivity, m2/s; α f , Fin material's thermal diffusivity, m2/s; ρ, PCM's density, kg/m3; δ, Fin's
lf
half thickness, m; σ, (Aspect ratio)
lc
Subfixes: l, Liquid; p, PCM; s, Solid; f , Fin

Int J Energy Res. 2019;1–23. wileyonlinelibrary.com/journal/er © 2019 John Wiley & Sons, Ltd. 1
2 KOTHARI ET AL.

1 | INTRODUCTION simulated the transient behavior of encapsulated PCM in

a container. Good agreement was obtained between the
The demand of energy in various sectors such as indus- results and test data.
trial, utility, and commercial sectors vary with time.1 Low thermal conductivity is the most significant
Thermal energy storage (TES) systems have the ability drawback of PCM. Numerous investigations were
and are promising candidates for smoothing the tempo- reported in literature to enhance the thermal conductiv-
ral variation of energy demands. In addition, TES can ity of PCMs with high‐conductivity metal particles,
be used to utilize the renewable energy sources and metal foam, and fins.27-30 Fins are widely used to
improve the energy efficiency. TES can be achieved by enhance the thermal conductivity of PCMs and heat
using sensible or latent heat storage. The latter is more transfer in LHTES system. Several experimental and
preferable to the former because of constant temperature numerical studies were reported on heat transfer
phase transition and high storage density.2,3 Phase enhancement of LHTES systems with internal fins.31,32
change materials (PCMs) based Latent heat thermal During the experimental study performed by Velraj
energy storage (LHTES) for storing thermal energy has et al30 with paraffin wax as PCM inside a cylindrical‐
various applications such as photovoltaic (PV) system, vertical shell type heat exchanger, the solidification time
solar water heating, electronic cooling and thermal com- was found to be inversely proportional to the number of
fort. Important advancements in LHTES systems, mate- internal longitudinal fins. Costa et al31 numerically ana-
rials and applications have recently been reviewed by lyzed the thermal performance of LHTES system with
various researchers.4-9 Also, novel PCMs such as pork and without internal fins using enthalpy‐based method.
fat have recently been used by Nizetic et al10 for PV It was reported that the fins significantly enhance the
systems. thermal conductivity of PCMs. Also, the authors
Heat transfer in LHTES materials, namely, PCMs, is reported the melt fraction and energy stored of the
complex in nature because of nonlinear problem and PCM. The melting and solidification process of PCM
moving solid–liquid interface. The moving solid–liquid storage was investigated through both experimentally
problem is commonly termed as moving boundary prob- and numerically by Lamberg et al.32 The investigations
lem (MBP).1-3 The analytical solution for the one‐ were performed with and without fins, and the heat
dimensional (1D) MBP termed as Stefan's problem.11,12 capacity method was found to be more accurate com-
Analytical solutions of 1D MBP for different boundary pared with the enthalpy method.
conditions were reported by Alexiades and Solomon13 Several analytical studies were reported for analyz-
and Ozisik.12 These include the perturbation method, ing the PCM melting and solidification time with
the Megerlin method, the quasi‐stationary approxima- fins.1-3,33-38 The analytical and numerical solutions for
tion, Kantorovich method, and the heat balance integral the location of interface of solid and liquid and temper-
method (HBIM). ature variation during solidification3 and melting33 of
Several theoretical14-22 and numerical23-26 models were PCM with finned two‐dimensional storage were
proposed to obtain the PCM temperature variation during reported for constant end‐wall temperature condition.
solidification and melting. Solomon et al14 reported a solu- It was reported that geometry of PCM and material
tion for rectangular TES with convective air boundary con- properties affect the PCM solidification time. Talati
dition by using quasi‐stationary approximation. Lu15 et al2 reported the solution for PCM solidification in
presented 1D conduction‐based heat transfer model to two‐dimensional PCM storage having internal fins by
investigate melting of PCM in Cartesian coordinates. employing both analytical and numerical techniques
Imposed heat flux was taken on the PCM top surface, while for constant heat flux wall condition. Position of inter-
the PCM bottom surface was taken as convective air envi- face of solid and liquid and temperature variation of
ronment. Authors reported PCM temperature distribution fin in a rectangular‐finned PCM storage for convective
in terms of Biot number and time. Kothari et al17 developed air environment at the boundaries were predicted by
a mathematical model for obtaining the PCM temperature Mosaffa et al.1 Bauer37 reported analytical solution for
distribution in an annulus by using Variational formula- PCM's solidification in fin geometries using effective
tion or quasi‐steady approximation. In addition to this, medium approach with constant end‐wall temperature
several numerical methods were reported to obtain the boundary condition. The model is useful to optimize
solution for moving boundary problem. In general, the finned TES. Some of the studies that report the analyti-
enthalpy method is commonly used because of the elimi- cal solution for PCM storage having internal fins are
nation of explicit treatment of boundary conditions at summarized in Table 1. This includes the solution for
solid–liquid interface.23 Zivkovich and Fuji24 introduced location of interface of solid and liquid and solution
enthalpy formulation‐based mathematical model and for various regions.
TABLE 1 Solid–liquid interface locations and solutions for regions 1 and 2 for rectangular phase change material storage with internal fins reported by various researchers
KOTHARI

Boundary
ET AL.

S. No Reference Problem Adapted Conditions Solution for Region 1 (along X Direction) Solution for Region 2 (along Y Direction) Validation

1 "
  
1. Mosaffa et al Solidification in a Constant imposed ∞ cm exp½− λ2m þ η2 τ ψ Yes (numerical)
hð T m − T ∞ Þ θðζ ; tÞ ¼ ∑m¼1 cos ðλ m ζ Þ þ sin ð λ m ζ Þ þ
finite rectangular convective air T s ðx; tÞ ¼ N ð λm Þ λ2m þ η2 λm
ks þ hX ðt Þ
PCM storage with environment "
∞ β X ðt Þ cosðβ X ðt ÞÞ − sinðβ X ðt ÞÞ ψ ψþη
m A1 expðηζ Þ þ A2 expð−ηζ Þ where A1 ¼ − A2 ,
internal fins ½x þ 2 ∑ m
m  ψ−η ψ−η
n¼1 H s þ β2m þ H 2s X ðt Þ
2
# ψη½1 þ expðηÞ − ψ2 ½1 − expðηÞ 2ψλm

 eð−αs βm tÞ sinðβm ðX ðt Þ − x Þ A2 ¼
2  , tanðλm Þ ¼ − 2 ,
× β2m þ H 2s × 2 ψ þ η2 sinhðηÞ þ 4ψη coshðηÞ ψ − λ2m
β2m sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
 
ks T m þ hX ðtÞT ∞ 2ks ðT m − T ∞ Þ
þ Y ðt Þ ¼ t where Y(t) is solid–liquid interface,
ks hX ðt Þ ρs L
where βm cot (βmX(t)) = − Hs k s l2f h
η¼ ,ψ ¼ H f lf ¼ lf .
Y ðt ÞDk f kf
"
 #
2. Talati et al2 Solidification in a Constant imposed −8q″w ∞ exp −αs β2n t    Yes (numerical)
T s ðx; t Þ ¼ ∑ cos ð β x Þ −2ψ 1 2

n θ¼ expð−ξτ Þ þ ψ þ þ ψ η2 − η
finite rectangular end‐wall heat ks π n¼1 ð2n−1Þ2 ξ 6 ξ
 where (
PCM storage with flux ″
 ″

q x q X ðt Þ ∞ 2ψ ð−1Þnþ1 − 1


− w þ w þ Tm ; þ ∑n¼1 exp − ξ þ n2 π 2 τ
internal fins ks ks ξþn π 2 2

ð2n − 1Þπ

βn ¼ 2ψξ ð−1Þnþ1 − 1
2X ðt Þ þ 2 2 2 2
cosðnπηÞ;
n π ðξ þ n π Þ
hpffiffiffiffiffiffiffiffiffiffi i12 pffiffiffiffiffiffi
−1
ks l2f x
Y ðt Þ ¼ 2 1þ2ste 2 αs t where Y(t) is solid–liquid interface, ξ ¼ ,η ¼ .
Y ðt ÞDkf lf
2
3. Lamberg et al3 Solidification in a Constant imposed Sx¼0 ¼ 2λpffiffiffiffiffiffi
αs t λeλ where Sx = 0 = solid–liquid T = Tm + θ(Tw − Tm), pffiffiffi Yes (numerical)
coshððη − 0:5Þ νÞ 4
finite rectangular end‐wall St l cp ðT w − T m Þ θ¼ pffiffiffi − υτ
erf ðλÞ ¼ pffiffiffi ¼ pffiffiffi coshð0:5 νÞ πe
π −L π 2 2
PCM storage with temperature
interface location, λ = root of transcendental ∞ ð−1Þn eð−ð2nþ1Þ π τ Þ
internal fins ∑n¼0 

ð2n þ 1Þ½1 þ ν= ð2n þ 1Þ2 π 2
equation
cosðð2n þ 1Þπ ðη − 0:5ÞÞ;
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
γ ¼ 2ξkλ2 θτ where γ = solid–liquid interface location, ξ = modified Stefan
number
33 ffi 2
pffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffi
4. Lamberg et al Melting in a semi‐ Constant imposed Sx ðt Þ ¼ 2λpffiffiffiffiffi
αl t λeλ where Sx(t) = solid–liquid B 2 1 B 2 Yes (numerical)
ðT w − T m Þ eBt− X =A − eBt− X =A
infinite rectangular end‐wall St l cp ð Tw − TmÞ 2
erf ðλÞ ¼ pffiffiffi ¼ pffiffiffi " pffiffiffiffiffiffiffiffiffiffi#
π −L π pffiffiffiffiffiffiffiffiffiffiffi
PCM storage with temperature B 2 x BX 2 t
interface location, λ = root of transcendental 1 − e2 X =A þ erf pffiffiffiffiffi −
internal fins 2 At x
equation
pffiffiffiffiffiffiffiffiffiffi
ffi pffiffiffiffiffiffiffiffiffiffi!!)
2 BX2
=A erf x BX 2 t
þe pffiffiffiffiffi −
2 At x h
T f ðx; t Þ ¼ T m þ where B ¼

eBt ρcp f D

(Continues)
3
 4

TABLE 1 (Continued)

Boundary
S. No Reference Problem Adapted Conditions Solution for Region 1 (along X Direction) Solution for Region 2 (along Y Direction) Validation

2



5. Mosaffa et al34 Solidification in a Constant imposed Ts − Tm r 

∞ exp −αs γ m t ∞ 1 Yes (numerical)
¼ ln − ∑m¼1 θðr; t Þ ¼ ∑m¼1 exp − γ 2m − ν2 τ ψo ðγ m ; ηÞ
shell and tube convective air T∞ − Tm a N ðγ m Þ N ðγ m Þ
ηout where, ψo(γm, η) = Yo(γmη)
PCM storage with environment Ω2 ½J o ðγ m r ÞY o ðγ m ΩÞ−Y o ðγ m r ÞJ o ðγ m ΩÞ × × ∫ − ηφðηÞψo ðγ m ; ηÞdη
1 ηin
internal fins  
½Y o ðγ m ΩÞI 1 − J o ðγ m ΩÞI 2 g
Ω J1(γmηout) − Jo(γmη)Y1(γmηout)
1=ðH s Rin Þ þ ln
Rin
where Jo and Yo is Bessel function of zero order.
Where U0Y0(γmΩ) − W0J0(γmΩ) = 0

Abbreviation: PCM, phase change material.

KOTHARI
ET AL.
KOTHARI ET AL. 5

In case of TES systems, the melting and solidification

time of PCMs affect the thermal performance. In general,
the solidification time is nearly three times greater than
the melting time.3 The solidification behavior of PCM is
found to depend on various thermal and geometrical
parameters including the aspect ratio, fin thickness, and
thermophysical properties of fin and PCM.37 It is evident
from the literature that numerous studies have been
made that analyze the melting/solidification behaviors
of PCMs with internal fins inside storage system involv-
ing rectangular and cylindrical geometries.1-3,33-38 Most
of the two‐dimensional analytical studies employ separa-
tion of variable (SV) method for the analysis. These tech-
niques compute the solution by truncating the series of
FIGURE 1 A phase change material (PCM) storage with PCM
complex function. The truncation of infinite series of a
and internal horizontal fins
function may lead to errors in the solution. Therefore,
there is a need to develop theoretical models during solid-
ification of PCM with other methods involving important not affect the solidification of PCM as reported by
parameters. HBIM is one of many analytical methods Lamberg and Siren.3,33 In Region 2, transfer of heat to
used to solve conduction problems. HBIM has been used the PCM takes place by the wall as well as the fin. After
to solve Stefan problems because of its simplicity, short short duration, heat is mainly transferred by the fin. In
computation time while providing the solutions with rea- such case, the interface of solid and liquid is assumed
sonable accuracy.39 to move one dimensionally in direction Y. In solidifica-
Here, an effort has been made to obtain the behavior tion process, conduction is assumed to be primary mode
of PCM during solidification process in two‐dimensional of transfer of heat. The natural convection, at the inter-
rectangular PCM storage system with horizontal internal face of solid and liquid because of the difference in tem-
fins subjected to various end‐wall boundary conditions, perature of liquid PCM, has minimal effect on location
namely, constant heat flux, end‐wall temperature, and of interface of solid and liquid compared with heat con-
convective air environment. HBIM is employed to obtain duction in solid PCM.1 Therefore, the natural convection
the solution. The objective of the present work is to is not considered in the analysis. The heat transfer prob-
evaluate the temperature variation in PCM, predict loca- lem is complex because of the nonlinear and transient
tion of interface of solid and liquid, solidification time nature. The following assumptions are made for the
for various geometrical and thermal boundary condi- analysis.1-3
tions. In addition to this, correlations have been pro-
posed to predict solidification time as a function of (i) Initially, the temperature of PCM and fin is consid-
aspect ratio and end‐wall boundary condition for all ered to be the same as that of solidification temper-
the cases. ature (Tm).
(ii) The temperature variation of fin is 1D because of
small thickness, geometry and high thermal
2 | THE PHYSICAL PROBLEM conductivity.
(iii) Various properties of fin, PCM, and heat transfer
Figure 1 represents the schematic of PCM storage unit fluid (HTF) are considered to be constant irrespec-
with internal horizontal fins. Here, left and right walls tive of temperature.
are subjected to constant temperature, constant heat (iv) PCM is assumed to be homogeneous, and solidifica-
flux, and constant convective air environment condi- tion is considered to take place isothermally.
tions. A simplified 1D model is developed to estimate
the location of interface of solid and liquid and temper-
ature variation of fin. Because of symmetry of PCM stor-
age unit, a single symmetry cell was considered for the 3 | MATHEMATICAL
analysis and is shown in Figure 2. Two regions, namely, FORMULATION
Region 1 and Region 2, is considered for the analysis. In
Region 1, the transfer of heat is 1D, toward direction X. The present analytical model was developed in two parts.
Here, heat is transferred only by walls, and the fin does The problem can be treated as Stefan's problem13 in
6 KOTHARI ET AL.

3.1 | Constant imposed heat flux at the

walls

3.1.1 | Region 1, {0 ≤ x ≤ X(t)}

The transient heat conduction equation in one dimension

for the PCM (Figure 2A) can be expressed as follows:

∂2 T s 1 ∂T s
¼ ; 0 ≤ x ≤ X ðt Þ; t ¼ 0: (1)
∂x 2 αs ∂t

The initial and boundary conditions are given by:

T s ðx; 0Þ ¼ T m ; 0 ≤ x ≤ lf ; (2a)

dX ðt Þ ∂T s ðX ðtÞ; t Þ
−ρLp ¼ −ks ; (2b)
dt ∂x
∂T s
−k s ¼ q″w q″w < 0: (2c)
∂x
Here, X (t) is the distance of interface of solid and liquid
in direction X.
On integrating the governing Equation 1 from end‐
wall to solid–liquid interface with respect to x and utiliz-
ing the initial and boundary conditions, Equations 2a to
2c and 1 can be rewritten as follows:

d X ðtÞ
∫ ρcs ½T s ðx; t Þ − T m dx ¼ −ρLp X ′ ðt Þ − q″w ðt Þ; 0 < x < X ðt Þ:
dt 0
(3)

In order to solve the above equation and boundary condi-

tions (Equations 2a‐2c), a guess temperature profile is
assumed and is given as follows:

T s ðx; t Þ ¼ T m þ Aðt Þ½x − X ðt Þ þ Bðt Þ½x−X ðt Þ2 ; (4)

where, A (t) and B (t) are the constants; given by:

" sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi#
Lp 4cs q″w X ðt Þ
A ðt Þ ¼ −1 þ 1 − ; (5)
2cs X ðt Þ k s Lp

FIGURE 2 Schematic of different regions in a two‐dimensional " sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi#2

phase change material (PCM) storage with different boundary Lp 4cs q″w X ðtÞ
B ðt Þ ¼ 2 −1 þ 1− ; (6)
wall conditions (A) constant heat flux, (B) constant temperature, 8cs fX ðt Þg k s Lp
and (C) constant convective air environment [Colour figure can be
viewed at wileyonlinelibrary.com] Utilizing Equations 4 to 6, polynomial form solution for
various parameters such as cs, q″w , ks, ρ, Lp, X(t), and t.
This can be written in the following form.
Region 1, which is specifically solvable MBP. In Region
2, transfer of heat is assumed to follow the direction Y, A1 fX ðtÞg4 þ B1 fX ðtÞg3 þ C 1 fX ðtÞg2 þ D1 fX ðt Þg þ E 1 ¼ 0
and energy balance is required for the analysis. The (7)
problem is further divided based on different boundary where,
conditions. This has been elaborated in subsequent
2
sections. A1 ¼ cs q″w ; (8)
KOTHARI ET AL. 7

in which δ is the half thickness of fin and Y (t) is the dis-
B1 ¼ 16k s cs Lp q″w ρ2 − 1 ; (9)
tance of interface of solid and liquid in direction Y.
h i The initial and boundary conditions of the fin are

2
C1 ¼ 4k s 16k s L2p − 3cs q″w t − k s ρ2 L2p ; (10) given by:

T f ðx; 0Þ ¼ T m 0 ≤ x ≤ lf ; (17a)
D1 ¼ 96k 2s Lp q″w t; (11)



2 ∂T f ð0; t Þ ∂T f lf ; t
E 1 ¼ 6ks q″w t ; (12) −k f ¼ kf ¼ q″w : (17b)
∂x ∂x
In order to obtain X (t), we need values of various proper-
On integrating the governing equation along the half‐
ties of PCM such as cs, q″w , ks, ρ,Lp, and t. The values of
length of fin with respect to x, Equation 15 can be rewrit-
various properties of PCM are summarized in Table 3.
ten as follows:

3.1.2 | Region 2, {0 ≤ x ≤ Y(t)} l =2  lf =2

1 d
f

 ∂T f ðx; t Þ
∫ T f − T m dx ¼ (18)
Here, the interface of solid and liquid moves only in one αf dt 0 ∂x 0
direction that is in direction Y. In such case, one can con- l =2
f


sider the heat transfer takes place because of conjugate − M ∫ T f − T m dx:
effect of heat transfer through fin. Figure 3 shows the 0

arbitrary differential element dx, which is detached from

The solution of Equation 15 for different boundary condi-
the PCM storage. An energy balance for the element dx
tions has been reported by several researchers.1-3 While,
can be written as follows:
for HBIM method, one needs to assume the guess temper-
−q″x þ q″xþdx þ q″in ¼ E_ st ; (14) ature profile, and this can be written as follows:
h pffiffiffiffi pffiffiffiffi i
where, E_ st is the rate of change of energy storage per
T f ðx; tÞ ¼ T m þ t Aðt Þe M x þ Bðt Þe− M x þ Cðt Þ : (19)
unit area of the fin. Substituting rate equations, the
energy balance equation (Equation 14) can be rewritten
Applying the initial and boundary conditions given by
as follows:
Equations 17a and 17b of the fin in Equation 18, and
∂2 T f
 1 ∂T f substituting in Equation 19, one can obtain the tempera-
− M Tf − Tm ¼ ; (15)
∂x 2 αf ∂t ture variation of fin and is given by:
where,
q″w
T f ðx; t Þ ¼ T m þ tP1 þ pffiffiffiffiffi
ks kf M
M¼ ; (16) " pffiffiffiffi pffiffiffiffi pffiffiffiffi pffiffiffiffi #
k f δY ðt Þ
e M x þ e− M x þ e M ðx−lf Þ þ e− M ðx−lf Þ
pffiffiffiffi pffiffiffiffi ; (20)
e M lf − e− M lf

where,

2α q″ Q  pffiffiffiffi 
P1 ¼
f w  − p1ffiffiffiffiffi e M lf − 1
lf kf 1 þ Mαf t lf M
 p ffi 
ffiffiffi (21a)
R1
− pffiffiffiffiffi 1 − e− M lf ;
lf M
" pffiffiffiffi #
q″w 1 þ e− M lf
Q1 ¼ pffiffiffiffiffi pffiffiffiffi pffiffiffiffi ; (21b)
k f t M e M lf − e− M lf

" pffiffiffiffi #
FIGURE 3 Energy balance in a differential element of the phase q″w 1 þ e M lf
change material storage [Colour figure can be viewed at R1 ¼ pffiffiffiffiffi pffiffiffiffi pffiffiffiffi (21c)
k f t M e M l f − e− M l f
wileyonlinelibrary.com]
8 KOTHARI ET AL.

d X ðt Þ

ks ∫ ρcs T j ðx; t Þ − T m dx
M¼ : (21d) dt 0
k f δY ðt Þ ¼ −ρLp X ′ ðt Þ − k s T s ð0; t Þ for 0 < x < X ðt Þ: (25)

The location of interface of solid and liquid in Region 2, A(t) and B(t) are the constants in assumed temperature
that is in direction Y, Y (t) is defined following the profile, and their values are obtained as follows:
Megerlin method13 as below.
Lp pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1=2 A ðt Þ ¼ −1 þ 1 − 2st s ; (26)

1 þ 2Ste−1 pffiffiffiffiffiffi c s X ðt Þ
Y ðt Þ ¼ 2 αs t Megerlin method; (22)
2
Lp pffiffiffiffiffiffiffiffiffiffiffiffiffi
2
where, Ste denotes the Stefan number and is provided by: B ðt Þ ¼ 2 −1 þ 1−2st s : (27)
2cs fX ðt Þg


cs T m − T f
Ste ¼ : (23) Utilizing Equations 4 and 25 to 27, the position of inter-
Lp
face of solid and liquid is given as follows:
It may be noted by using Equations 20 to 23, one can
obtain temperature variation of the fin and the location "sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi #
3ð−1 þ 2st s − 1 − 2st s Þ pffiffiffiffiffiffi
of interface of solid and liquid Y (t). X ðt Þ ¼ 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi αs t : (28)
−7 − 2st s þ 1 − 2st s

3.2 | Constant imposed temperature at

the walls
3.2.2 | Region 2, {0 ≤ x ≤ Y(t)}
Figure 2B shows the location of interface of solid and liq-
uid during solidification in two different regions for con- It may be noted that solution of Region 2 for this model is
stant imposed temperature at the end‐walls. obtained by considering similar governing equation
(Equation 15), initial condition (Equation 17a) and
assumed temperature profile (Equation 19).
3.2.1 | Region 1, {0 ≤ x ≤ X(t)} Apart from above equations, following boundary
conditions for constant imposed end‐wall heat flux
Here, the solution of the model is obtained by considering (Equation 17b) is replaced with constant end‐wall tem-
similar governing equation (Equation 1), initial condition perature boundary condition and is given as follows:
(Equations 2a and 2b), and assumed temperature profile
(Equation 4).

T f ð0; t Þ ¼ T f lf ; t ¼ T w : (29)
In addition to this, boundary condition for constant
imposed end‐wall heat flux (Equation 2c) is replaced with Applying the initial and boundary conditions of the fin
constant end‐wall temperature boundary condition and is (Equations 17a and 29) in Equation 18 and substituting
given as follows: in Equation 19, one can obtain the temperature variation

 of fin and is given by:
T s ð0; t Þ ¼ T s lf ; t ¼ T w < T m : (24)
T f ðx; t Þ ¼ T w
On integrating the governing equation (Equation 1) from h pffiffiffiffi pffiffiffiffi  pffiffiffiffi i
end wall to interface of solid and liquid with respect to x þtP2 e M ðx−lf Þ þ e− M x − 1 þ e− M lf ; (30)
and utilizing the initial and boundary conditions given by
Equations 2a, 2b, 24, and 1 can be rewritten as follows: where,

pffiffiffiffiffi
ðT m − T w Þlf M
P2 ¼ 2  pffiffiffiffi 3 (31a)
pffiffiffiffi n pffiffiffiffi o n pffiffiffiffi o l pffiffiffiffi ffin
p ffiffiffi
ffi o Mα t e − M lf
− 1
f M f
2t 4e− M lf eð M lf =2Þ − 1 þ 1 − eð− M lf =2Þ − 1 þ e− M lf þ 5
2 1 þ Mαf t
KOTHARI ET AL. 9

The value of constants A(t) and B(t) assumed in tempera-

ks
M¼ (31b) ture profile (Equation 4) are obtained as follows:
k f δY ðt Þ
Lp pffiffiffiffiffiffiffiffiffiffiffiffiffi

Here, the interface of solid and liquid Y (t) is defined sim- A ðt Þ ¼ 1 þ st s − 1 ; (37)
cs X ðt Þ
ilar to the Equation 22. Also, variation in temperature of
the fin and the location of interface of solid and liquid Y
(t) can be obtained using Equations 22 23 and 30‐(31a).
−Lp pffiffiffiffiffiffiffiffiffiffiffiffiffi
2
B ðt Þ ¼ 2 1 þ st s −1 ; (38)
3.3 | Convective air environment at the 2cs fX ðt Þg
walls where, sts is Stefan number for solid PCM and is given by
cs ðT m − T ∞ Þ
The location of interface of solid and liquid during solid- st s ¼ .
Lp
ification in two different regions for convective air envi-
Utilizing Equations 4, 35, 37, and 38, the solid–liquid
ronment boundary condition is shown in Figure 2C.
interface position is given as follows:
 pffiffiffiffiffiffiffiffiffiffiffiffiffi 
3.3.1 | Region 1, {0 ≤ x ≤ X(t)} 3ht st s − 2f1 − 1 þ st s g
pffiffiffiffiffiffiffiffiffiffiffiffiffi :
X ðt Þ ¼ (39)
ρcs 7 − st s − 1 þ st s
In addition to governing Equation 1, initial and boundary
conditions (Equations 2a‐2b) and assumed temperature
profile (Equation 4), constant air environment boundary 3.3.2 | Region 2, {0 ≤ x ≤ Y(t)}
condition is used instead of constant imposed end‐wall
heat flux (Equation 2c) to obtain the solution. And the In addition to Governing Equation 15, initial condition
solution is given as follows: (Equation 17a) and assumed temperature profile (Equa-
tion 19), constant air environment boundary condition
−∂T s ð0; t Þ
þ H s fT s ð0; t Þ − T ∞ g ¼ 0: (34) is used instead of constant imposed end‐wall heat flux
∂x
(Equation 17b) for obtaining the solution. The solution
The governing Equation 1 can be rewritten by integrating is given as follows:
it with respect to x from end wall to interface of solid and
liquid, and substituting the initial and boundary condi- −∂T f ð0; t Þ  
þ H f T f ð0; t Þ − T ∞
tions (Equations 2a, 2b, and 34), we can obtain the fol- ∂x

∂T f lf ; t 
 
lowing relation: ¼ þ H f T f lf ; t − T ∞ ¼ 0: (40)
∂x
d X ðt Þ
∫ ρcs ½T s ðx; t Þ − T m dx ¼ −ρLp X ′ ðt Þ The temperature variation of fin is obtained by substitut-
dt 0 (35)
ing Equations 17a, 18, and 40 in Equation 19 given as fol-
− h½T ∞ − T ð0; t Þ for 0 < x < X ðt Þ:
lows:
Assuming, h pffiffiffiffi pffiffiffiffi
T f ðx; t Þ ¼ T ∞ þ tP3 Q3 e M x þ e− M x
k s þ hX ðt Þ ( !
≈ 1; (36a) pffiffiffiffiffi pffiffiffiffiffii#
2k s þ hX ðt Þ M M
þ Q3 × −1 −1− ; (41)
Hf Hf
2hX ðt Þ
≈ 1: (36b) where,
k s þ hX ðt Þ

ðT m − T ∞ Þlf
P 3 ¼ "( pffiffiffiffiffi ! pffiffiffiffiffi) ( pffiffiffiffi ! )#; (42a)
M M   n −
pffiffiffiffi o e M lf − 1  pffiffiffiffi 
t R3 −1 −1− lf þ 2S3 þ R3 þ e M lf pffiffiffiffiffi þ S3 1 þ e M l f
Hf Hf M
10 KOTHARI ET AL.

pffiffiffiffiffi pffiffiffiffiffi
M −pffiffiffi ffi M pffiffiffiffi solid and liquid and solid fraction of PCM obtained from
e M lf
þ þ 1 − e− M lf the present prediction with the two‐dimensional numeri-
Hf Hf
Q3 ¼ pffiffiffiffiffi pffiffiffiffiffi ; (42b) cal solution of Talati et al,2 Lamberg et al,3 and Mosaffa
M M lfp ffiffiffi
ffi M pffiffiffiffi
e þ − 1 þ e M lf et al1 for constant heat flux, end‐wall temperature, and
Hf Hf
convective air boundary condition, respectively, and the
αf H f t results of the validation are shown in Figure 4A‐C and
R3 ¼ ; (42c) Table 2A‐C. For the sake of brevity, the location of solid
1 þ Mαf t
and liquid interface in directions X and Y are compared
ks at σ<1. The average error, between the present predic-
M¼ ; (42d) tions and two‐dimensional numerical results, in
k f δY ðt Þ
obtaining the location of interface of solid and liquid in
Here, the interface of solid and liquid Y(t) is defined sim-
direction X is found to be 0.69, 1.73, and 1.06 mm for con-
ilar to the Equation 22. Also, temperature variation of fin
stant heat flux, end‐wall temperature, and convective air
and location of interface of solid and liquid Y(t) can be
boundary condition, respectively. While, the average
obtained using the Equations 22, 23, and 41 to 42d.
error in obtaining the location of interface of solid and
liquid in direction Y is found to be 0.19, 0.39, and
4 | C O M P A R I S O N O F P R E S EN T 0.57 mm in direction Y for constant heat flux, end‐wall
PREDICTION WITH NUMERICAL temperature, and convective air boundary condition,
RESULTS respectively.
The solid fraction of PCM during solidification for
In order to verify the present analytical model, efforts varied range of aspect ratio and different boundary condi-
have been made to compare the location of interface of tion are compared with the two‐dimensional numerical

FIGURE 4 Comparison of solid–liquid interface locations obtained by present model with the existing numerical solutions for (A) σ = 2,
t = 800 seconds, constant heat flux boundary condition; (B) σ = 5, t = 1085 seconds, constant end‐wall temperature boundary condition; and
(C) σ = 2, t = 1000 seconds, constant convective air environment boundary condition [Colour figure can be viewed at wileyonlinelibrary.com]
KOTHARI ET AL. 11

TABLE 2A Comparison of present analytical model and two‐ TABLE 3 Thermal properties of phase change material and fin
dimensional numerical model2 for constant heat flux boundary material2
condition in different test cases
Property Paraffin Aluminum Fin
Solid
Density, ρ (kg/m3) 1450 2770
Fraction Solid Fraction
Predicted Obtained from Specific heat, c (J/kg‐K) 3600 875
from Two‐dimensional Thermal conductivity, k (W/m‐K) 0.6 177
Aspect Present Numerical Error
Latent heat of fusion, Lp (J/kg) 162 000 …
Ratio Time Study (εanal) Model,2 (εnum) (I‐II)
(σ) (s) (I) (II) × 100% Melting/solidification 28 …
temperature, Tm (°C)
0.5 400 0. 43 0.35 6.7
1 400 0. 30 0.25 3.3
2 400 0. 26 0.21 1.4 5 | R ESULTS A ND DISCUSSION

Here, an analytical model was developed to analyze the

temperature distribution in both PCM and fin during
TABLE 2B Comparison of present analytical model and two‐ the solidification of PCM inside the rectangular container
dimensional numerical model3 for constant temperature boundary with horizontal internal fins. The fraction of the fin
condition in different test cases length to the half‐height of symmetry section (l f /lc) is
considered as the cell aspect ratio (σ) (see Figure 1). Var-
Solid Fraction
ious cases with different values of σ are considered for the
Solid Obtained from
Fraction Two‐ analysis. For all the cases, half of the thickness of fin is
Predicted dimensional considered to be 0.5 mm. Initially, the liquid PCM is at
Aspect from Present Numerical Error the solidification temperature (Tm) inside the PCM stor-
Ratio Time Study (εanal) Model,3 (εnum) (I‐II) age unit. HBIM was employed to solve the conduction
(σ) (s) (I) (II) × 100% equation associated with different boundary conditions.
Paraffin, 0.2 724 0.90 0.79 10.6 Different guess profiles for temperature distribution were
Paraffin, 1 1085 0.29 0.23 6.1 used for the analysis. Furthermore, the present predic-
Paraffin, 5 1085 0.61 0.53 7.1 tions were compared with the results obtained by existing
analytical and numerical techniques.

TABLE 2C Comparison of present analytical model and two‐

dimensional numerical model1 for convective air boundary condi-
5.1 | Constant imposed heat flux at the
tion in different test cases
wall

Solid Solid Fraction The PCM and fin material considered for the present
Fraction Obtained from study are salt hydrate Climsel C28 and aluminum, respec-
Predicted Two‐dimensional tively. Their properties are summarized in Table 3. The
Aspect from Present Numerical Error different values of σ are taken as 0.5, 1, and 2 for the anal-
Ratio Time Study (εanal) Model,1 (εnum) (I‐II)
ysis. Dimensions of PCM storage unit are shown in
(σ) (s) (I) (II) × 100%
Table 4. Constant heat flux, q″w = −1000 W/m2, was taken
0.5 350 0.27 0.40 −13.0
for solidification of PCM in all the cases.
1 500 0.27 0.36 −9.4 Figures 5A‐F depict the comparison of location of
2 500 0.20 0.23 −3.0 interface of solid and liquid obtained from the present
mathematical model with the 1D analytical and two‐

TABLE 4 Dimensions of phase change material storage unit

solutions and is summarized in Table 2A. Present predic- Dimensions Case 1 Case 2 Case 3
tion exhibits good agreement with the numerical solu-
Length of fin, l f (mm) 10 15 20
tion.1-3 The maximum error between the present
Half height of cell, lc (mm) 20 15 10
prediction and the two‐dimensional numerical solution
is found to be ±13%. Therefore, the model can be useful lf 0.5 1 2
Cell aspect ratio, σ ¼
in designing the PCM‐based TES systems. lc
12 KOTHARI ET AL.

FIGURE 5 Comparison of solid–liquid interface locations obtained by present model with the one‐dimensional analytical and two‐
dimensional numerical solutions by Talati et al2 withq″w = −1000 W/m2 for (A) σ = 0.5, t = 400 seconds; (B) σ = 0.5, t = 800 seconds; (C)
σ = 1, t = 400 seconds; (D) σ = 1, t = 800 seconds; (E) σ = 2, t = 400 seconds; and (F) σ = 2, t = 800 seconds [Colour figure can be viewed at
wileyonlinelibrary.com]

dimensional numerical solutions of Talati et al2 for differ- directions X and Y. Figure 5 shows the comparison of
ent values of time steps (t) and σ. Present predictions results obtained by using integral method in the present
exhibit good agreement with the 1D results obtained by study with various available mathematical models. Here,
employing SV method and two‐dimensional results the results of 1D analytical and two‐dimensional numer-
obtained by finite difference method.2 ical models of Talati et al2 are considered for comparison.
It is evident that geometry of PCM storage changes In all the cases, sharp corners are formed on the interface
the solidification rate. When σ is very less than unity, of solid and liquid at the point near both the wall and the
heat flows only from wall to interface of solid and liquid, fin surfaces. This is because of the high temperature at
in direction X. When σ is much larger than unity, heat the corners. The maximum error in solid–liquid interface
flows only from fin to interface of solid and liquid, in between the two sets of results was found to be 0.96 mm
direction Y. In case of σ equal to unity, heat flows from in direction X and 0.40 mm in direction Y at σ = 2 and
both wall and fin to interface of solid and liquid, in the t = 400 seconds. However, present model exhibits
KOTHARI ET AL. 13

excellent agreement with the two‐dimensional numerical

model when the value of σ = 0.5.
The temperature distributions of fin obtained from
the present analytical model are compared with the
results obtained by employing SV and finite difference
methods,2 and are shown in Figures 6A‐C. It is observed
from Figure 6A that the present predictions exhibit good
agreement with the results obtained by SV method. The
maximum error in temperature difference between both
sets of results was found to be very negligible (~0.095°C)
at σ = 2 and t = 800 seconds. It has been reported that
1D heat transfer occurs in case of short fins when the
value of σ is much smaller than unity. The temperature
variation along the fin length is found to be low because
of high thermal conductivity and constant heat flux.
However, the temperature variation increases as l f
increases. The current results exhibit excellent agreement
with the numerical solutions obtained by Talati et al.2
Also, the present predictions agree well with the results
obtained by employing SV method.
Solid fraction of PCM usually denotes the amount of
solidified PCM in the symmetric storage cell for a given
duration. The solid fraction of PCM changes from entirely
solid (ε = unity) to completely liquid (ε = 0) and is
expressed as follows1-3:

2X ðt Þ½lc − δ−Y ðt Þ þ lf Y ðt Þ
ε¼ : (43)
½lc − δlf

Time of PCM solidification in a rectangular storage hav-

ing fins for each symmetry cell was obtained using the
analytical formulation for different cases and is presented
in Figure 7. It is noticed from Figure 7 that the PCM
solidification is fast for smaller values of σ. Initially, the
PCM solidification is faster because of small thickness of
solid PCM layer. However, as the time progresses, the
thickness of solid PCM layer increases, which increases
the thermal resistance and solidification rate decreases.
Figure 8 shows variation in PCM's total solidification
time for different aspect ratio at constant volume of
PCM and TCE. It may be noted from the figure that solid-
ification time is less for lower aspect ratio of TES system
for constant volume of PCM and TCE. This may be due
to the fact that with the decrease in aspect ratio, distance
FIGURE 6 Comparison of fin temperatures obtained by present
between fins increases and a huge amount of heat is lib- model with the one‐dimensional analytical and two‐dimensional
erated from the end walls. The PCM's total solidification numerical solutions by Talati et al2 with q″w = −1000 W/m2 for (A)
time is found to be less for lower aspect ratio σ = 0.5, t = 400 seconds and t = 800 seconds; (B) σ = 1,
(σ = 0.125) compared with σ = 8 (Figure 8). The reduc- t = 400 seconds and t = 800 seconds; and (C) σ = 2, t = 400 seconds
tion in solidification time was found to be approximately and t = 800 seconds [Colour figure can be viewed at
93% by reducing the aspect ratio from 8 to 0.125 for a wileyonlinelibrary.com]
given imposed heat flux at the walls. In addition to this,
with the increase in the imposed heat flux from 500 to
14 KOTHARI ET AL.

FIGURE 8 Variation in total solidification time of phase change

material for different imposed heat flux at the end walls [Colour
figure can be viewed at wileyonlinelibrary.com]

FIGURE 7 (A) Variation of solid fraction of phase change

material (PCM) with time, (B) solidification time of PCM for FIGURE 9 Comparison of computed solidification temperature
different cell aspect ratios [Colour figure can be viewed at with those predicted by the proposed correlation [Colour figure
wileyonlinelibrary.com] can be viewed at wileyonlinelibrary.com]

2000 W/m2, the solidification time reduces by approxi- 5.2 | Constant imposed temperature at
mately 60% for the same aspect ratio. the walls
In the present study, a correlation has been proposed
for solidification time (ts) in terms of aspect ratio and In this study, pure n‐octadecane paraffin is used as the
end‐wall boundary condition and expressed in PCM with aluminum as fins and commercial grade salt
Equation 44. hydrate (Climsel 23) used as the PCM with steel as fins
for the analysis. Table 5 summarizes the properties of fin

−0:5976 materials and PCM. Different values of σ such as 0.2, 1,
t s ¼ 125077:867ðσ Þ0:6768 q″w : (44)
and 5 are considered for the analysis. Table 6 shows the
dimensions of PCM storage unit for which the constant
The proposed correlation for ts is valid for the range of wall temperature was considered (Figure 2).
0.125 ≤ σ ≤ 8 and the error band to anticipate 99% test It is observed from Figure 10 that the aspect ratio (σ)
data of the proposed correlation for solidification time of PCM storage unit affects the rate of solidification. The
(ts) is found to be ±5% as shown in Figure 9. heat flows only in direction X, ie, from wall to interface of
KOTHARI ET AL. 15

TABLE 5 Thermal properties of phase change materials and fin materials3

Paraffin Salt Hydrate Aluminum Steel

Property n‐Octadecane Climsel23 Fin Fin

Density, ρ (kg/m3) 777 1480 2713 7854

Specific heat, c (J/kg‐K) 2660 2660 960 434
Thermal conductivity, k (W/m‐K) 0.149 0.6 180 60.5
Latent heat of fusion, Lp (J/kg) 241 360 148 000 … …
Melting/solidification temperature, Tm (°C) 28 23 … …
Viscosity, μ (kg/m‐s) 0.00385 … … …
−1
Thermal expansion coefficient, α (K ) 0.001 … … …

TABLE 6 Dimensions of phase change material storage, end‐wall temperature, and materials used in different test cases

Cases Case 1 Case 2 Case 3 Case 4 Case 5 Case 6

PCM material Paraffin Paraffin Paraffin Salt hydrate Salt hydrate Salt hydrate
Fin material Aluminum Aluminum Aluminum Steel Steel Steel
End‐wall temperature, Tw, (°C) 13 13 13 8 8 8
Length of fin, l f (mm) 10 50 50 10 50 50
Half height of cell, lc(mm) 50 50 10 50 50 10
lf 0.2 1 5 0.2 1 5
Cell aspect ratio, σ ¼
lc

solid and liquid when the value of σ is less than unity. be −1.389°C at σ = 1 and t = 4226 seconds. The present
However, the heat flows only in direction Y, ie, from wall results agree well with the existing numerical results.3
to interface of solid and liquid when the reverse is true. In The variation of solid fraction with time for two differ-
case of σ equals to unity, heat flows from both wall to ent PCMs for σ = 0.2 and σ = 5 was evaluated using Equa-
solid–liquid interface and fin to solid–liquid interface, tion 43 and is demonstrated in Figure 12. It may be noted
ie, both in X and Y directions. Figures 10A‐E depict the that for the same aspect ratio the solidification rate is
comparison of locations of interface of solid and liquid higher for salt hydrate as compared with paraffin because
obtained from the present 1D analytical predictions and of better heat conductivity of salt hydrate compared with
two‐dimensional numerical solutions of Lamberg et al3 paraffin. This may be due to the difference in density and
for different values of t and σ. The maximum error in latent heat of fusion of two PCMs. The reduction in respec-
predicting interface of solid and liquid by present model tive solidification rate is found to be 73.4% and 50% for
was found to be 1.82 and 0.37 mm in X and Y directions, σ = 0.2 and σ = 5 for paraffin and salt hydrate.
respectively, at σ = 5 and t = 1085 seconds for Total solidification time of PCM for solidification time
n‐octadecane. While, the maximum error for salt hydrate of paraffin for different aspect ratio at constant volume of
was found to be 1.23 and 1.43 mm in X and Y directions, PCM and TCE is presented in Figure 13. The solidification
respectively, at σ = 1 and t = 4226 seconds. However, time is found to be less for lower aspect ratio of TES sys-
present model exhibits excellent agreement with the tem. This may be because the distance between fins
two‐dimensional numerical model for σ = 0.2. increases with the decrease in aspect ratio and a large
Here, efforts were also made to estimate the variation of amount of heat is discharged from the end walls. The
temperature in fin, and the same are presented in Figure 11 PCM's total solidification time is remarkably less for lower
. The temperature variation of fin obtained from the pres- value of aspect ratio (σ = 0.05) compared with σ = 5
ent analytical model was compared with the results (Figure 13). The reduction in solidification time was found
obtained by employing Fourier's ring theory and finite dif- to be approximately 99% by reducing the aspect ratio from
ference method.3 The maximum error in temperature by 5 to 0.05 for same imposed heat flux at the walls. In addi-
the present model was found to be very negligible tion to this, with the decrease in the imposed end wall
(0.25°C) at σ = 5 and t = 1085 seconds for n‐octadecane. from 18 to 8°C, the solidification time reduces by approx-
While, the maximum error for salt hydrate was found to imately 54% for the same aspect ratio.
16 KOTHARI ET AL.

FIGURE 10 Comparison of solid–liquid interface locations obtained by present model with the one‐dimensional analytical, one‐
dimensional numerical and two‐dimensional numerical solutions by Lamberg et al3 for (A) σ = 0.2, t = 724 seconds; (B) σ = 5,
t = 1085 seconds; (C) σ =0.2, t = 169 seconds; (D) σ = 1, t = 4226 seconds; and (E) σ = 5, t = 1127 seconds [Colour figure can be viewed at
wileyonlinelibrary.com]

In the present investigation, a correlation has been 5.3 | Convective air environment
proposed for solidification time (ts) in terms of aspect boundary conditions
ratio and end‐wall boundary condition and expressed in
Equation 46. In this section, paraffin and aluminum are taken as PCM
and fin material, respectively. Their properties are sum-
marized in Table 7. Here, different values of σ such as
t s ¼ 329:3757ðσ Þ1:0021 ðT w Þ1:0512 : (46)
0.5, 1, and 2 are considered. The dimensions of PCM stor-
age unit are listed in Table 8. Constant air boundary con-
The proposed correlation for ts is valid for the range of dition was considered for PCM storage unit (Figure 2).
0.05 ≤ σ ≤ 8 and the error band to anticipate 99% test HTF temperature and coefficient of heat transfer during
data of the proposed correlation for solidification time convection for PCM solidification were assumed constant
(ts) is found to be ±10% as shown in Figure 14. in all the cases (Figure 2C). Air was chosen as HTF
KOTHARI ET AL. 17

FIGURE 11 Comparison of fin temperatures obtained by present model with the one‐dimensional analytical, one‐dimensional numerical,
and two‐dimensional numerical solutions by Lamberg et al3 for (A) σ = 0.2, t = 724 seconds; (B) σ = 5, t = 1085 seconds; (C) σ = 0.2,
t = 169 seconds; (D) σ = 1, t = 4226 seconds; and (E) σ = 5, t = 1127 seconds [Colour figure can be viewed at wileyonlinelibrary.com]

having the temperature of 10°C and the value of coeffi- in the directions X and Y) at σ = 1. The comparison of
cient of heat transfer during convection at the wall was location of interface of solid and liquid obtained from
taken as 65 W/m2/K1.1 the present analytical investigation with the 1D mathe-
Rate of solidification is influenced by the aspect ratio matical and two‐dimensional numerical solutions of
of PCM storage unit. It can be seen from Figure 15 that Mosaffa et al1 for different values of t and σ is presented
when σ is less than unity, heat flows only in direction in Figure 15A‐F. The maximum error in interface of
X, that is, from wall to solid–liquid interface. However, solid–liquid by the present model was found to be
heat flows only in direction Y, that is, from fin to inter- 1.1477 mm in direction X at σ = 1 and t = 1000 seconds
face of solid and liquid, when the value of σ is greater and 0.973 mm in direction Y at σ = 2 and t = 1000 seconds.
than unity. While, heat flows from both wall to solid– The present model agrees well with the two‐dimensional
liquid interface and fin to interface of solid–liquid (both numerical model for σ = 0.5.
18 KOTHARI ET AL.

FIGURE 13 Variation in total solidification time of phase change

material for different imposed temperature at the end walls [Colour
figure can be viewed at wileyonlinelibrary.com]

FIGURE 12 (A) Variation of solid fraction of phase change

material (PCM) with time, (B) solidification time of PCM for
different cell aspect ratios and PCMs [Colour figure can be viewed
at wileyonlinelibrary.com] FIGURE 14 Comparison of computed solidification temperature
with those predicted by the proposed correlation [Colour figure can
be viewed at wileyonlinelibrary.com]

Here, the fin's temperature variation was determined

and compared with the results obtained by employing TABLE 7 Thermal properties of phase change material and fin
SV and enthalpy method,1 and is presented in Figure 15 material1
A‐C. It is noticed from Figure 16 that the present predic- Aluminum
tions exhibit good agreement with the result obtained by Property Paraffin Fin
SV method. The maximum error in temperature differ- Density, ρ (kg/m3) 830 2770
ence by the present model was found to be negligible at
Specific heat, c (J/kg‐K) 1920 (solid) 875
σ = 2 and t = 1000 seconds. It may also be observed that
3260 (liquid)
the length of fin significantly changes the temperature
Thermal conductivity, 0.514 (solid) 177
variation. For the small values of l f , the temperature dif-
k (W/m‐K) 0.224 (liquid)
ference is small, while it increases with the values of l f .
This is attributed to the 1D heat transfer in the small Latent heat of fusion, 251 000 …
Lp (J/kg)
length of a fin and its higher thermal conductivity. The
present predictions agree well with the existing numerical Melting/solidification 32 …
temperature, Tm (°C)
solutions.1
KOTHARI ET AL. 19

The variation of solid fraction with time for σ = 0.5 Figure 18 shows variation in PCM total solidification time
and σ = 2 was evaluated using Equation 43 and is demon- for different aspect ratio at constant volume of PCM and
strated in Figure 17. It is noticeable from Figure 17 that TCE. Different values of temperature of HTF and coeffi-
the PCM solidification is fast for smaller values of σ. cient of heat transfer during convection has been taken
for the analysis. Here, the reduction in solidification time
TABLE 8 Dimensions of phase change material storage unit
was found to be approximately 88% by reducing the
Dimensions Case 1 Case 2 Case 3 aspect ratio from 8 to 0.125 for a given imposed convec-
tive environment at the walls. Also, with the increase in
Length of fin, l f (mm) 10 15 20
the coefficient of heat transfer from 55 to 75 W/m2/K1,
Half height of cell, lc (mm) 20 15 10
the solidification time reduces by approximately 26% at
lf 0.5 1 2 T∞ = 10°C for the same aspect ratio. However, the reduc-
Cell aspect ratio, σ ¼
lc
tion in solidification time is obtained as approximately

FIGURE 15 Comparison of solid–liquid interface locations obtained by present model with the one‐dimensional analytical and two‐
dimensional numerical solutions by Maosaffa et al1 with h = 65 W/m2/K and T∞ = 10°C for (A) σ = 0.5, t = 350 seconds; (B) σ = 0.5,
t = 700 seconds; (C) σ = 1, t = 500 seconds; (D) σ = 1, t = 1000 seconds; (E) σ = 2, t = 500 seconds; and (F) σ = 2, t = 1000 seconds [Colour
figure can be viewed at wileyonlinelibrary.com]
20 KOTHARI ET AL.

FIGURE 17 (A) Variation of solid fraction of phase change

material (PCM) with time, (B) solidification time of PCM for
different cell aspect ratios at h = 65 W/m2/K and T∞ = 10°C
[Colour figure can be viewed at wileyonlinelibrary.com]

55% by reducing T∞ from 20to 10°C at h = 65 W/m2/K1

for the same aspect ratio. Here, it may be noted the
decrease of HTF temperature (T∞) has more effect on
solidification time compared with increase in coefficient
of heat transfer during convection (h).
The proposed correlation for solidification time (ts) in
terms of aspect ratio and end‐wall boundary condition is
FIGURE 16 Comparison of fin temperatures obtained by present expressed in Equation 47.
model with the one‐dimensional analytical and two‐dimensional
numerical solutions by Mosaffa et al1 with h = 65 W/m2/K and
T∞ = 10°C for (A) σ = 0.5, t = 350 seconds and t = 700 seconds; (B)
t s ¼ 125077:867ðσ Þ0:5014 ðhÞ−0:5814 ðT m −T ∞ Þ−1:0002 : (47)
σ = 1, t = 500 seconds and t = 1000 seconds; and (C) σ = 2, The proposed correlation for ts is valid for the range of
t = 500 seconds and t = 1000 seconds [Colour figure can be viewed 0.125 ≤ σ ≤ 8 and Tm > T∞. The error band to anticipate
at wileyonlinelibrary.com] 99% test data of the proposed correlation for solidifica-
tion time (ts) is found to be ±10% as shown in
Figure 19.
KOTHARI ET AL. 21

• Closed form expressions were obtained for various

parameters such as temperature distribution, solid–
liquid interface, and solid fraction.
• Present predictions agree well with the results
obtained by employing SV and finite difference
methods.
• The solidification time is found to reduce by 93% and
99% by reducing the aspect ratio from 8 to 0.125 for
constant heat flux boundary condition and reducing
the aspect ratio from 5 to 0.2 for constant wall temper-
ature boundary condition, respectively.
• For the same aspect ratio, the solidification time is
found to be more for salt hydrate compared with paraf-
fin. By reducing the aspect ratio from σ = 0.2 to σ = 5,
the solidification time was found to reduce by 50% and
FIGURE 18 Variation in total solidification time of phase change 73% for paraffin and salt hydrate, respectively. The
material for different convective air conditions at the end walls solidification time was found to decrease by 88% for
[Colour figure can be viewed at wileyonlinelibrary.com]
paraffin by reducing the aspect ratio from 8 to 0.125
for constant convective air boundary condition.
• Correlations for the solidification time is given by

−0:5976
t s ¼ 125077:867ðσ Þ0:6768 q″w , ts = 329.3757
(σ)1.0021(Tw)1.0512, and ts = 125077.867(σ)0.5014(h)
−0.5814
(Tm − T∞)−1.0002 for constant heat flux, end‐
wall temperature, and convective boundary condi-
tions, respectively. The error band to anticipate 99%
test data of the proposed correlation for solidification
time (ts) is found to be ±10%.

The present analytical model can contribute to better

understanding of solidification of PCM in a finite storage
with internal fins. The present study reports the effect of
various geometrical and thermal boundary conditions on
solidification time. The present analytical model can be
useful in optimum design of PCM‐based TES system with
internal fins.
FIGURE 19 Comparison of computed solidification temperature
with those predicted by the proposed correlation [Colour figure can
be viewed at wileyonlinelibrary.com] ACKNOWLEDGEMENT
The financial support received from DST, Government of
India (DST Sanction No.: DST/TMD/MES/2k17/65) is
6 | CONCLUSION gratefully acknowledged by the authors. Also, the finan-
cial support received from DST‐INSPIRE Fellowship pro-
A theoretical model has been proposed to estimate vari- gram (IF170534) is admired by the first author.
ous parameters such as solid fraction, interface of solid
and liquid, solidification time, and temperature variation
ORCID
during the solidification of PCM in a two‐dimensional
LHTES with horizontal internal fins. Various boundary Rohit Kothari https://orcid.org/0000-0002-7678-6690
conditions, namely, constant heat flux, end‐wall tempera-
ture, and convective air boundary conditions, were con-
RE FER EN CES
sidered at the walls of the storage system. HBIM was
used to obtain the solution. The conclusions obtained 1. Mosaffa AH, Talati F, Rosen MA, Tabirizi HB. Approximate
from the present study are elaborated below. analytical model for PCM solidification in a rectangular finned
22 KOTHARI ET AL.

container with convective cooling boundaries. Int Comm Heat 20. Kalaiselvam S, Veerappan M, Aaron AA, Iniyan S. Experimental
Mass Tran. 2012;39(2):318‐324. and analytical investigation of solidification and melting charac-
2. Talati F, Mosaffa AH, Rosen MA. Analytical approximation for teristics of PCMs inside cylindrical encapsulation. Int J Therm
solidification process in PCM storage with internal fins: imposed Sci. 2008;47(7):858‐874.
heat flux. Heat Mass Tran. 2011;47(4):369‐376. 21. Saha SK, Dutta P. Performance analysis of heat sinks with phase
3. Lamberg P, Siren K. Approximate analytical model for solidifi- change materials subjected to transient and cyclic heating. Heat
cation in a finite PCM storage with internal fins. App Math Tran Eng. 2015;36(16):1349‐1359.
Model. 2003;27(7):491‐513. 22. Kalaiselvam S, Parameshwaran R, Harikrishnan S. Analytical
4. Dincer I, Rosen MA. Thermal Energy Storage: System and Appli- and experimental investigations of nano particles embedded
cations. 2nd ed. Chichester: John Wiley & Sons; 2003. phase change materials for cooling applications in modern
buildings. Renew Energy. 2012;39(1):375‐387.
5. Zalba B, Marin JM, Cabeza LF, Mehling LF. Review on thermal
energy storage with phase change materials, heat transfer anal- 23. Voller VR. Fast implicit finite‐difference method for analysis of
ysis and applications. Appl Therm Eng. 2003;23(3):251‐283. phase change problems. Numer Heat Tr B Fund.
1990;17(2):155‐169.
6. Farid MM, Khudhair AM, Cabeza LF, Mehling H. A review on
phase change energy storage: materials and applications. Energ 24. Zivkovich B, Fujii I. An analysis of isothermal phase change of
Conver Manage. 2004;45(9‐10):1597‐1615. phase change material within rectangular and cylindrical con-
tainers. Sol Energ. 2001;70(1):51‐61.
7. Sharma SD, Sagara K. Latent heat storage materials and sys-
tems: a review. Int J Green Energ. 2005;2(1):1‐56. 25. Esen M. Thermal performance of solar aided latent heat store
8. Baljit SSS, Chan H, Sopian K. Review of building integrated used for space heating by heat pump. Sol Energ.
applications of photovoltaic and solar thermal systems. J Clean 2000;69(1):15‐25.
Prod. 2016;137:677‐689. 26. Arıcı M, Bilgin F, Nižetić S, Papadopoulos AM. Phase change
9. ISarbu I, Dorca A. Review on heat transfer analysis in thermal material based cooling of photovoltaic panel: a simplified
energy storage using latent heat storage systems and phase numerical 2 model for the optimization of the phase change
change materials. Int J Energ Res. 2018;1‐36. https://doi.org/ material layer and general economic 3 evaluation. J Clean Prod.
10.1002/er.4196 2018;189:738‐745. https://doi.org/10.1016/j.jclepro.2018.04.057

10. Nižetić S, Arıcı M, Bilgin F, Grubišić‐Čabo F. Investigation of 27. Mettawee EBS, Assassa GMR. Thermal conductivity enhance-
pork fat as potential novel phase change material for passive ment in a latent heat storage system. Sol Energ.
cooling applications in photovoltaics. J Clean Prod. 2007;81(7):839‐845.
2018;170:1006‐1016. 28. Siahpush A, O'brien J, Crepeau J. Phase change heat transfer
11. Carslaw H, Jaeger J. Conduction of Heat in Solids. 2nd ed. New enhancement using copper porous foam. ASME Trans Heat
York: Oxford University Press; 1959. Trans. 2008;130:1614‐1616.
12. Ozisik MN. Heat Conduction. 2nd ed. New York: John Wiley 29. Lacroix M. Study of heat transfer behavior of a latent heat ther-
and Sons; 1993. mal energy storage unit with a finned tube. Heat Mass Tran.
1993;36(8):2083‐2092.
13. Alexiades V, Solomon AD. Mathematical Modelling of Melting
and Freezing Processes. Washington DC: Hemisphere Publish- 30. Velraj R, Seeniraj RV, Hafner B, Faber C, Schwarzer K. Heat
ing Corporation; 1993. transfer enhancement in latent heat storage system. Sol Energ.
14. Solomon AD, Wilson DG, Alexiades V. The quasistationary 1999;62:19‐28.
approximation for the Stefan problem with the convective 31. Costa M, Buddhi D, Oliva A. Numerical simulation of latent
boundary condition. Int J Math Math Sci. 1993;7:549‐563. heat thermal menergy storage system with enhanced heat con-
15. Lu TJ. Thermal management of high power electronics with phase duction. Energ Conver Manage. 1998;39(3‐4):319‐330.
change cooling. Int J Heat Mass Tran. 2000;43(13):2245‐2256. 32. Lamberg P, Lehtiniemi R, Henell AM. Numerical and experi-
16. Chakraborty S. Dutta, P. (2003) analytical solution for heat mental investigation of melting and freezing processes in
transfer during cyclic melting and freezing of phase change phase change material storage. Int J Therm Sci.
material used in electronic of electrical packaging. J Electron 2004;43(3):277‐287.
Packag. 2003;125(1):126‐133. 33. Lamberg P, Siren K. Analytical model for melting in a semi‐
17. Kothari R, Sahu SK, Kundalwal SI. Comprehensive analysis of infinite PCM storage with an internal fin. Heat Mass Tran.
melting and solidification of a phase change material in an 2003;39(2):167‐176.
annulus. Heat Mass Tran. 2018. https://doi.org/10.1007/s00231‐ 34. Mosaffa AH, Talati F, Tabrizi HB, Rosen MA. Analytical model-
018‐2453‐9.(In press) ling of PCM solidification in a shell and tube finned thermal
18. Xia Q, Chen Y, Yang C, Zhang T, Zang Y. A new model of phase storage for air conditioning systems. Energ Buildings.
change process for thermal energy storage. Int J Energ Res. 2012;49:356‐361.
2018;42(12):3877‐3887. 35. Lamberg P. Approximate analytical model for two phase solidi-
19. Jiji LM, Gaye S. Analysis of solidification and melting of PCM fication problem in a finned phase‐change material storage.
with energy generation. Appl Therm Eng. 2006;26(5‐6):568‐575. Appl Energy. 2004;77(2):131‐152.
KOTHARI ET AL. 23

36. Mosaffa AH, Talati F, Rosen MA, Tabrizi HB. Phase change 39. Sahu SK, Das PK, Bhattacharyya S. A three‐region conduction‐
material solidification in a finned cylindrical shell thermal controlled rewetting analysis by the heat balance integral
energy storage: an approximate analytical approach. Therm method. Int J Therm Sci. 2009;48(11):2100‐2107.
Sci. 2013;17(2):407‐418.
37. Bauer T. Approximate analytical solutions for the solidification
of PCMs in fin geometries using effective thermophysical prop- How to cite this article: Kothari R, Das S, Sahu
erties. Int J Heat Mass Tran. 2011;54(23‐24):4923‐4930. SK, Kundalwal SI. Analysis of solidification in a
38. Taghilou M, Talati F. Analytical and numerical analysis of finite PCM storage with internal fins by employing
PCM solidification inside a rectangular finned container with heat balance integral method. Int J Energy Res.
time‐dependent boundary condition. Int J Therm Sci. 2019;1–23. https://doi.org/10.1002/er.4363
2018;133:69‐81.