Transient behavior analysis of the melting of nanoparticle-enhanced phase change material inside a rectangular latent heat storage unit

Accepted Manuscript Research Paper Transient behavior analysis of the melting of nanoparticle-enhanced phase change material inside a rectangular latent heat storage unit Radouane Elbahjaoui, Hamid El Qarnia PII: DOI: Reference:

S1359-4311(16)32566-2 http://dx.doi.org/10.1016/j.applthermaleng.2016.10.115 ATE 9317

To appear in:

Applied Thermal Engineering

Received Date: Revised Date: Accepted Date:

13 July 2016 14 October 2016 19 October 2016

Please cite this article as: R. Elbahjaoui, H. El Qarnia, Transient behavior analysis of the melting of nanoparticleenhanced phase change material inside a rectangular latent heat storage unit, Applied Thermal Engineering (2016), doi: http://dx.doi.org/10.1016/j.applthermaleng.2016.10.115

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Transient behavior analysis of the melting of nanoparticle-enhanced phase change material inside a rectangular latent heat storage unit Radouane Elbahjaoui*, Hamid El Qarnia Cadi Ayyad University, Faculty of Sciences Semlalia, Department of Physics, P.O 2390, Fluid Mechanics and Energetic Laboratory (affiliated to CNRST, URAC 27), Marrakesh, Morocco

*

Corresponding author

Email: [email protected] Phone: +212605706053 Fax: +212524436769

1

Transient behavior analysis of the melting of nanoparticle-enhanced phase change material inside a rectangular latent heat storage unit Abstract Melting of Paraffin Wax (P116) dispersed with Al2O3 nanoparticles in a rectangular latent heat storage unit (LHSU) is numerically investigated. The storage unit consists of a number of vertical and identical plates of nanoparticle-enhanced phase change material (NEPCM) separated by rectangular channels through which a heat transfer fluid flows (HTF: Water). A two-dimensional mathematical model is considered to numerically investigate the heat and flow characteristics of the LHSU. The heat transfer and fluid flow during the melting process were formulated using the enthalpy-porosity method. The finite difference forms of the governing equations are obtained using the finite volume approach. The numerical model has been validated by experimental, theoretical and numerical results available in the literature. The effects of the aspect ratio of NEPCM slabs, volumetric fraction of nanoparticles, Reynolds number and Rayleigh number on the flow characteristics and thermal performance of the storage unit were investigated. A correlation including all these investigated control parameters has been developed to predict the time required for the complete melting of NEPCM. Keywords: Nanoparticle-enhanced phase change material (NEPCM), Phase change material (PCM), Nanoparticles, Latent heat storage unit (LHSU), Thermal energy storage, Melting, CFDs.

2

Nomenclature A aspect ratio, H/ (d/2) cp specific heat at constant pressure (J/kg. K) d thickness of the PCM slabs (m) dn diameter of nanoparticles (m) dimensionless diameter of nanoparticles, dn/lo dn Dh hydraulic diameter (m), Dh=2w e / (w + e)) dimensionless hydraulic diameter, Dh/lo Dh e depth of the storage unit (e = 1m) f liquid fraction g acceleration of gravity (m/ s2) h specific sensible enthalpy (J/ kg) H height of the PCM slabs (m) k thermal conductivity (W/ m. K) dimensionless thermal conductivity, k/km,l k l0 characteristic length, l0  H d / 2 p pressure (Pa) q dimensionless heat transfer rate Q dimensionless heat Ra Rayleigh number Re Reynolds number Ste Stefan number Su source term for u momentum equation

Sv Sh T t u, v U, V x,y w

source term for v momentum equation source term for energy equation in terms of enthalpy temperature (K) time (s) velocity components (m/s) dimensionless velocity components Cartesian coordinates (m) thickness of the HTF channels (m)

Greek symbols thermal diffusivity (m2/ s),   k/ρcp  dimensionless thermal diffusivity, α/αm.l  ϕ volumetric fraction of nanoparticles β coefficient of thermal expansion of liquid PCM (K-1) θ dimensionless temperature τ dimensionless time μ dynamic viscosity (N. s/ m2) ν kinematic viscosity (m2/ s) Δh latent heat (J/ kg) ρ density (kg/ m3) ε storage efficiency Ψ dimensionless stream function Subscripts e f

inlet heat transfer fluid (HTF) 3

l lat melt n nm m o P s sen

liquid latent melting nanoparticle Nanoparticle-enhanced phase change material (NEPCM) PCM outlet node point solid sensible

Superscripts k iteration k 0

old value at previous time step 1

Introduction

Due to their high energy storage capacity, the latent heat storage units using phase change materials (PCMs) have gained great attention over the last three decades. However, the phase change materials (PCMs) are characterized by a low thermal conductivity which limits the heat exchange rate and slowed their melting/solidification rate. To overcome this drawback and improve the heat transfer rate during the solid-liquid phase change process, several methods have been suggested in the literature. These methods include the dispersion of nanoparticles with a very high conductivity in PCMs [1], integration of metal matrix in PCMs [2], use of multiple PCMs [3, 4] and porous matrix [5]. Among these methods, the dispersion of high conductivity nanoparticles in PCM is the main topic of this study. In recent years, a number of experimental, numerical and theoretical works have been performed to investigate the thermal performance of nanoparticles enhanced phase change materials. Dhaidan [6] presented a review on the experimental, numerical and theoretical investigations of the melting of NEPCM inside various shape containers including rectangular, spherical and cylindrical enclosures. The author reported the effect of the geometrical parameters and operational conditions on the heat transfer and melting characteristics. The results show that the increase of the amount of nanoparticles dispersion in PCMs enhances heat transfer rate and reduces the melting time. Khodadadi and Hosseinizadeh [7] studied numerically the improving of the functionality of phase change material (PCM: water) by dispersion of nanoparticles (copper) during the solidification process with natural convection in a square cavity. The results of this study showed that the use of NEPCM improves the heat release compared to a basic PCM. Ho and Gao [8] experimentally studied the melting in a square enclosure filled with PCM (n-octadecane) dispersed with alumina nanoparticles (Al2O3). They reported that the influence of the solid subcooling through the un-melted PCM can be further enhanced by dispersion of high conductivity nanoparticles in the enclosure. They also indicated that the natural convection heat transfer in the liquid region degrades significantly with the increase of the mass fraction of nanoparticles when compared with natural convection corresponding to a base PCM. Sebti 4

et al. [9] conducted a 2D numerical study to investigate the heat transfer enhancement during melting process in a square enclosure through dispersion of nanoparticles. These authors reported that the dispersion of nanoparticles decreases the melting time. They also reported that the increase of the difference between the melting point and the hot wall temperature leads to a significant decrease of the melting time. Hossain et al. [10] developed a two dimensional thermal model to study the melting process of NEPCM inside a rectangular porous medium. The results show that the NEPCM melts rapidly within the lower porosity medium. Shuja et al. [11] performed an experimental and numerical investigations on the thermal characteristics of PCM (n-octadecane) integrated with metallic meshes. They examined the effect of the mesh geometry on the melting time for two different metallic mesh arrangements. The results show that the melting time of PCM is the shortest for the triangular mesh geometry then followed by the rectangular and hexagonal meshes. Sciacovelli et al. [12] conducted a numerical analysis of the melting of NEPCM in a shell-and-tube latent heat storage unit (LHSU). They demonstrate that the melting time is reduced by 15% when the volume fraction of dispersed nanoparticles is 4%. Kashani et al. [13] numerically studied the effects of the nanoparticles volume fraction and cold wall temperature during the solidification of NEPCM filled in a rectangular enclosure. They found that a large solid fraction is obtained for a smaller wall temperature and a high nanoparticles volume fraction. Hosseinizadeh et al. [14] conducted a numerical study of the melting of NEPCM inside a spherical container. The simulation results of this study show a pronounced potential of using phase change material with dispersed nanoparticles. Jourabian et al. [15] numerically examined the melting of Cu/water NEPCM inside a cylindrical horizontal annulus by using the enthalpy-based Lattice Boltzmann method. Their results show that owing to the enhancement of the thermal conductivity, the effect of the volume fraction of nanoparticles on the heat transfer rate is more significant. Their results also indicate that the temperature distribution in the NEPCM is faster than that of a base PCM. Ranjbar et al. [16] numerically studied the heat transfer enhancement in the LHSU through dispersion of nanoparticles during the solidification process. They reported that the increase in the volume fraction of the dispersed nanoparticles substantially increases the heat transfer rate. They also found that the dispersion of high conductivity nanoparticles promotes the heat transfer by conduction in NEPCM.Elbahjaoui et al. [17]investigated numerically the melting of PCM dispersed with alumina nanoparticles in a thermal storage unit composed of a number of vertical and identical slabs separated by HTF fluid flow. They evaluated the effect of the volumetric fraction of nanoparticles, HTF inlet temperature and mass flow rate on the heat transfer and performance of the storage unit. The results show that the operational conditions should be designed to achieve a significant enhancement of the thermal performance of the storage unit.Valan et al. [18] numerically investigated the performance enhancement of the paraffin wax dispersed with nano-alumina particles in a concentric double pipe heat exchanger. Their results indicate that the nano-alumina particles enhance the low thermal conductivity of the paraffin wax and thus greatly improve the charge and discharge rates of thermal energy. Mahdi and Nsofor [19] investigated numerically the solidification of PCM dispersed with alumina nanoparticles inside a triplex-tube storage unit. They evaluated the effect of the nanoparticles’ dispersion on the total solidification time, flow structure and heat transfer characteristics. The results show that the dispersion of nanoparticles with a volumetric concentration of 3% and 8% reduces the total solidification time by 8% and 20%, 5

respectively.Fan and Khodadadi [20] conducted an experimental investigation of enhanced thermal conductivity of cyclohexane-based NEPCM with copper oxide nanoparticles during the freezing process. Their results show that the measured thermal conductivity of the liquid phase is enhanced by the increase of the nanoparticles concentration, whereas the thermal conductivity data of the solid phase exhibit a non-monotonic enhancement when the concentration exceeds 2%. Das et al. [21] investigated numerically the melting behavior of PCM dispersed with graphene nano-sheets in a vertical shell-and-tube latent heat storage unit. Water acting as a HTF flows in the inner tube and transfers heat to the PCM filled in the shell space. The results show that the inclusion of graphene nano-sheets decreases the melting time due to the enhancement of the thermal conductivity of PCM. Lokesh et al. [22] experimentally investigated the melting and solidification processes of NEPCM in a stainless steel container. Their results showed higher enhancement of NEPCM thermal conductivity both in the liquid and solid phases. They also indicate that a reduction in the melting and solidification times of 42% and 26% can be achieved by using NEPCM with a concentration of 0.9 % and 0.3 %, respectively. Bechiri and Mansouri [23] developed an analytical solution for the melting and solidification of nano-enhanced PCM filled in a cylindrical latent heat storage unit. They used the variables separation and the exponential integration function to solve the energy equation based on pure conduction in PCM. They also conducted a parametric study to examine the effect of the volumetric heat generation on the thermal performance of the storage unit during both melting and solidification processes. To our knowledge, there are no comprehensive studies published in the literature for the melting of the nanoparticle-enhanced phase change material in rectangular enclosures heated by heat transfer fluid under laminar forced convection. In the present work, the modeling and numerical investigation of a latent heat storage unit (LHSU) composed of several vertical slabs filled with Al2O3/Paraffin Wax is presented. A laminar heat transfer fluid flow (water) circulates between the nanoparticle-enhanced PCM slabs and transfers heat to the storage unit. This study is motivated by the need to enhance the thermal performance of such thermal energy storage units due to their pronounced potential use in solar energy applications. The first objective of the present study is to investigate the effects of the nanoparticles’ volumetric fraction, aspect ratio of NEPCM slabs, Reynolds and Rayleigh numbers on the thermal behaviour and performance of the thermal energy storage unit. The second objective is to develop correlations including all these investigated control parameters to predict the time required for the complete melting of NEPCM. Such correlations are useful for engineering applications as they present guidelines to design latent heat storage units and optimize their thermal performances. 2 Problem formulation 2.1 Description of the studied physical model As shown in Fig 1a, the latent heat storage unit (LHSU) consists of a number of vertical and identical plates filled with NEPCM (Al2O3/P116). A heat transfer fluid (HTF: water), with 6

mass flow rate m and inlet temperature Tf,e>Tmelt, flows between the plates and transfers heat by laminar forced convection to NEPCM. For reasons of symmetry, the study of the overall system can be reduced to the study of a half NEPCM plate and a half HTF channel as shown in Fig. 1b. The height and thickness of the rectangular NEPCM plates are H and d, respectively. The thickness of the HTF channels is w. Initially, the plates contain solid NEPCM at a temperature equal to the melting point, Tmelt. At τ >0, the HTF starts to flow between the NEPCM plates, which triggers the melting process. The water inlet temperature, Tfe, is fixed at a constant value higher than the melting point, Tmelt. 2.2 Governing equations The modeling of the thermal and flow characteristics of the LHSU is made by using the governing conservation equations of energy for HTF, continuity, momentum and energy equations for NEPCM. The flows of HTF and liquid NEPCM are laminar and twodimensional. The HTF and the liquid phase of NEPCM are incompressible and Newtonian. For the HTF (water) circulating in channels, the flow is assumed hydro-dynamically fully developed. The thermo-physical properties are assumed to be constant except for the density difference according to the Boussinesq approximation. The enthalpy–porosity method was used for the formulation of the phase change problem [24].In this method, the solid-liquid interface and the two phases of NEPCM (solid and liquid) are modeled as a porous medium. Consequently, the transient and two-dimensional set of equations can be expressed as follows: For HTF

Tf (vf Tf )  T T    ( f f )  ( f f ) t y x x y y where vf 

(1)

3 m x 2 (1  ( ) ). 2 N c f w e w/2

For NEPCM

(nm u) (nm v)  0 x y

(2)

(nm u) (nm uu) (nm vu) p  u  u      ( nm )  ( nm )  Su (3a) t x y x x x y y (nm v) (nm uv)  (nm vv) p  v  v      (nm )  ( nm )  Sv t x y y x x y y (nm h) (nm uh) (nm vh)  k nm h  k h    ( )  ( nm )  Sh t x y x cp,nm x y cp,nm y with 7

(3b)

(4)

(1  f ) 2 Su  C 3 u (f  b)

Sv  C

(5a)

(1  f )2 v  nm gnm (Tnm  Tmelt ) (f 3  b)

(5b)

h   c p ,nm dTnm

Sh  nm hnm

(5c)

f t

(5d)

In these relations u and v are the velocity components, f stands for the liquid fraction, b is a constant introduced to avoid the division by zero (b = 0.001), C is a mushy zone constant that extinguishes velocity in the solid zone (C = 1.6×106 kg m-3 s-1) and h stands for the sensible enthalpy. As explained in paper [24] dealing with the solid-liquid phase change process, the momentum equation is similar to that of the porous medium with a porosity of zero value in the solid NEPCM and porosity equal to value 1 in the liquid NEPCM. The terms Su and Sv in the momentum equations are the source terms used to cancel the velocity in solid NEPCM. They gradually reduce the velocity components from a finite value in the liquid (the local liquid fraction, f = 1) to zero in the full solid (f = 0). The source term, S h, in the energy equation is the driving element resulting from the phase change process and the local liquid fraction reflects its evolution. It is equal to zero in the full liquid and solid phases and takes a finite value in the cells undergoing a phase change process. The density of the NEPCM is given by:

nm  n  (1  )m

(6)

The specific heat and Boussinesq term of the NEPCM are defined as follows: (ρc p )nm  ( ρc p )n  (1  )( ρc p )m

(7)

()nm  (1  )()m  ()n

(8)

The latent heat of the NEPCM is evaluated according to the work of Khodadadi and Zhang [25] as follows :

(ρh)nm  (1  )( ρh)m

(9)

where ϕ is the nanoparticles’ volumetric fraction and subscripts n, m and nm stand for nanoparticles, PCM and NEPCM, respectively. The viscosity of the NEPCM containing a small spherical suspension of solid particles is given as follows[25]:

nm 

m (1   )2.5

(10)

8

Whereas the thermal conductivity of the NEPCM is defined as follows:

knm  keff  kd

(11)

where keff is the thermal conductivity of the stagnant NEPCM and kd corresponds to the conductivity enhancement term which represents the Brownian motion. They are expressed as follows:

keff 

kn  2km  2(km  kn ) km kn  2km  (km  kn )

(12)

kd  C  (  c p )nm u 2  v 2 dn

(13)

where dnisthe diameter of the nanoparticles (dn= 45 nm) and C  is an empirical constant evaluated according to the work of Wakao and Kaguei [26] ( C  =0.1). The governing equations are converted into dimensionless form using the following dimensionless variables and parameters: X

 t T  Tmelt x y u v , Y  ,  m,l2 , U  ,V  ,  , l0 l0 l0  m,l / l0  m,l / l0 Tf ,e  Tmelt

Ra 

g m l03 (Tf ,e  Tmelt )

 m m,l

, Prm 

m  p , Prf  f , P  ,  m,l f  nm ( m,l / l0 ) 2

 nm

c p,m (Tf ,e  Tmelt )   2m w  nm ,  f  f ,Ste  , Re  ,w  ,  m,l  m,l h m f N c (w  e) l0

d

k k   d H , H  , k f  f , k nm  nm , nm  nm and  nm  nm . l0 l0 k m,l k m,l m m

(14)

where Tf e is the HTF inlet temperature, Tmelt is the melting temperature, Ra is the Rayleigh number, Ste is the Stefan number and l0 is the characteristic length ( its square is proportional to the mass of NEPCM. It should be noted that, in the present study, the value of l0 is kept constant at 0.07062 m). The expression of the characteristic length is given by:

l0  H

d 2

(15)

The dimensionless height and thickness of the NEPCM slabs can be expressed as a function of the aspect ratio, A, as follows: H A

(16)

2 A

(17)

d

9

By using Eqs. (6)-(11), the expressions of nm ,  nm and  nm can be written as follows:

nm 

1   ( n n  1) 1   ( n  1)

(18)

nm 

(1  )2.5 1  (n  1)

(19)

 nm 

kn  2km  2(km  kn ) km  C * U 2  V 2 d n 1   ( n c p ,n  1)   kn  2km  (km  kn ) 

(20)

knm 

kn  2km  2(km  kn ) km  C * 1  (n c p ,n  1)  U 2  V 2 d n k n  2 k m  ( k m  k n )

(21)

where n 

cp,n n  k k d , cp,n  , n  n , k n  n , k m  m  f  (1  f )k m,s and d n  n m cp,m m k m,l k m,l l0

The set of the governing equations are written as follows: For HTF

f (Vf f )       (f f )  (f f )  Y X X Y Y where Vf 

(22)

3 Re X 2 Prf f (1  ( ) ). 2 Dh w/2

For NEPCM Continuity: U V  0 X Y

(23)

Momentum: U UU VU P  U  U     (Prm nm ) (Prm nm )  Su  X Y X X X Y Y

(24)

V UV VV P  V  V     (Prm nm ) (Prm nm )  Sv  X Y Y X X Y Y

(25)

Energy:

10

(nm ) (Unm ) (Vnm )        ( nm nm )  ( nm nm )  S  X Y X X Y Y

(26)

where S u  C

(1  f )3 (1  f )3 1 f 1 U , S   C V  Ra Prm  nm nm , S   ( ) v 3 3 ( f  b) ( f  b) Ste  1   (  n c p ,n  1)

and C  C

l0 2

 nm  m,l

2.3 Boundary and initial conditions The initial and boundary conditions of the present computational domain are given by the following relations:

  0: nm   f  0, U  V  0 X  0:

 f X

0

X  w / 2: k f

 f

Y  0:  f  1,

X

 f Y



(28)

 knm

nm , U V  0 X

(29)

nm  0, U  V  0 Y

X  w / 2  d / 2:

Y  H:

(27)

(30)

nm V  0, U  0, 0 Y X

(31)

nm  0, U  V  0 Y

(32)

As it can be seen in afore mentioned equations, the control parameters of the studied system are given as follows:

A, Re, Ra,  , w, Prm , Ste,  f , Pr f , n , c p,n , kn , km,s , k f , n and dn 3

Computational methodology

All of the above established equations have been solved using the control volume method. The power law scheme has been used to treat the convective term, while an explicit scheme has been used for the time integration. The SIMPLE algorithm was used to treat the pressurevelocity coupling in NEPCM momentum equations. The resulting algebraic equations have been solved using the Tri-diagonal Matrix iterative method (TDMA). The numerical scheme has been implemented into an in-house simulation code for numerical calculation. The

11

iterative calculation process is continued for every time step, until convergence is achieved when the following criteria are met.

Max

( f ) p k  ( f ) p k 1 ( f ) p

k 1

 106 (33)

Max in (i, j )  out (i, j )  108

Q f  Qnm

(34)

 5 103

Qf

(35)

where

in (i, j )  U (i, j )Y  V (i, j )X   out (i, j )  U (i  1, j )Y  V (i, j  1)X 

(36)

The dimensionless heat transmitted by HTF to NEPCM, Qf, and the dimensionless heat stored in NEPCM, Qnm, are defined as follows: Qf 

M1 N Re  f Prf R f (f ,e  f ,o )   R f (f (i, j)  f 0 (i, j))XY 4 i 1 j1

Qnm 

M



N

 (nm (i, j)  nm0 (i, j))X Y 

i  M 11 j 1

M

N

  ( f (i, j)  f

i  M 11 j 1

0

(i, j ))

(37)

X Y 1 Ste 1   ( n c p ,n  1) (38)

The different terms constituting the expressions of Q f and Qnm are defined as follows: 



Re  f Pr f

R f ( f ,e   f ,o ) : dimensionless energy lost (exchanged) by HTF during 4 its passage between the inlet and outlet sections. M1 N

 R i 1 j 1

f

( f (i, j )   f 0 (i, j ))X Y : dimensionless energy stored in HTF by sensible

heat. 

M

N

  (

i  M 11 j 1

nm

(i, j )   nm 0 (i, j ))X Y : dimensionless energy stored in NEPCM by

sensible heat. 

M

N

  ( f (i, j)  f

0

(i, j ))

i  M 11 j 1

X Y 1 : dimensionless energy stored in Ste 1   ( n c p ,n  1)

NEPCM by latent heat. with

12

M1

 f ,o 

V  i 1

f

f

( X , Y  H , )X M1

V X i 1

, Rf 

 f c p, f  nm c p ,nm

(39)

f

In order to check the grid size dependence results, the numerical simulations were conducted for three grid sizes which are 20×30, 40×90 and 45×110 with ϕ = 2 %, A = 6, Re = 40 and Ra = 7.95x107 (Ste = 0.1444). The values of the other control parameters are displayed in Table 1.It should be noted that the Rayleigh and Stefan numbers depend on the HTF inlet temperature. So, the variation of the HTF inlet temperature causes the simultaneous variation of the Rayleigh and Stefan numbers. Numerical simulations were also conducted for four time steps which are   1.195x10-4 (5 s), 2.39x10-4 (10 s), 4.78x10-4 (20 s) and 9.56x10-4 (40 s) to test the time step dependency of the results. The effects of the various grid sizes and time steps on the melting fraction at   0.04972 (t = 34.66 min) are shown in Tables 2 and 3, respectively. These results show that the relative deviation of the melting fraction between the grid sizes 40x90 and 45x110 is only 0.23 %. The relative deviation of the melting fraction between the dimensionless time steps 2.39x10-4 (10 s) and 1.195x10-4 (5 s) is only 0.033%. Therefore, the grid size 40x90 and time step 2.39x10-4 (10 s) are adopted in our numerical calculations. It should be noticed that this grid size is only valid for an aspect ratio fixed at A = 6.For others aspect ratios, other grid size dependence studies have been performed. In order to test the reliability of the computational code, the numerical predictions of the current work are compared with the numerical, analytical and experimental results published in literature. The first validation concerns the melting of Gallium in a rectangular enclosure differentially heated on vertical surfaces. The horizontal walls were adiabatic. The melting front position obtained by the current modelis compared with that of the numerical results of Brent et al. [27] and the experimental data of Gau and Viskanta [28] at t = 2 min, 6 min, 10 min and 17 min. The heated left wall and cold right wall are kept at temperatures of 311.15 K and 301.45 K, respectively. The comparative results are illustrated in Fig.2. The analysis of these results shows that the average relative error of melting front position between the current model and the experiment data is 13.2%, 10.9%, 7.23% and 4.5% at t = 2 min, 6 min, 10 min and 17 min, respectively . Consequently, a reasonably good agreement can be considered between the numerical and experimental results. A second validation was carried out regarding the solidification of NEPCM (water-copper) inside a differentially heated square enclosure with insulated horizontal walls. The temperature of the left wall is maintained at the freezing temperature (0°C) while the cold right wall temperature was fixed at 10°C lower than the freezing temperature of the NEPCM. The freezing process of NEPCM was considered for three volumetric fractions of nanoparticles: 0%, 10% and 20%. The time variation of the liquid fraction of the NEPCM was illustrated in Fig. 3 to compare the results obtained by the current model with those produced by Khodadadi and Hosseinizadeh [7]. The average relative error of liquid fraction between the present numerical model and that of Khodadadi and Hosseinizadeh [7] is 5.7%, 8.3% and 3.1% for a nanoparticles’volumetric fraction of 0%, 10% and 20%, respectively. Because the rectangular channel through which HTF flows is a part of the LHSU, another validation was carried out regarding forced convection heat transfer in a plane channel. Such a 13

validation allows testing the robustness of the computational code to solve the heat transfer by forced convection in this part of the LHSU. Two boundary conditions have been considered: the isothermal wall and the constant heat flux imposed on the wall. The thermo-physical properties of the HTF (water) circulating in the plane channel are presented in Table 4. The values of the width, height and depth of the plane channel are 0.0254 m, 20 m and 1 m, respectively. The mass flow rate of HTF is 0.29 kg/s (Re = 1000) and Tf,e = 38 °C. Fig.4 shows the evolution of the local Nusselt number, as predicted numerically by the current computational code, for both boundary conditions. This figure shows that when the flow becomes thermally developed, the local Nusseltnumber achieves the values 7.54 and 8.23 for an isothermal wall and a constant heat flux wall, respectively. These results are in good agreement with those published in literature [29]. 4

Results and discussions

Numerical simulations were carried out to investigate the thermal characteristics and performance of the NEPCM storage unit heated by laminar forced HTF flow. The HTF, PCM and nanoparticles used in the calculations are water, paraffin Wax P116 and alumina (Al2O3), respectively. Their corresponding thermal properties are summarized in Table 4. In this study, the effects of the nanoparticles’ volumetric fraction, ϕ, aspect ratioof NEPCM slabs, A, Reynolds number, Re, and Rayleigh number, Ra, on the flow and thermal characteristics of the LHSU were examined. The volumetric fraction of nanoparticles is varied between 0% and 8%, the aspect ratioof NEPCM slabs is varied between 0.5 and 16, the Reynolds number is changed between 20 and 2000, the Rayleigh number is varied between 7.95x107 and 2.63x108, and the Stefan number is comprised between 0.1444 and 0.47756. The other governing control parameters are fixed at the same values listed in Table 1. 4.1 Effect of the volumetric fraction of nanoparticles on streamlines and isotherms The streamlines describing the flow field in the liquid NEPCM at different times:   0.01, 0.0502, 0.1075, 0.1291, 0.1578, 0.1893 and 0.2668 for   0 % and   8 % are shown in Figs. 5a-5g and Figs. 6a-6g, respectively, with A = 6, Re = 40 and Ra = 7.95x107 (Ste = 0.1444). With the start of the HTF circulation between the NEPCM slabs, a part of the heat lost by the HTF is transmitted to NEPCM thereby causing the beginning of the melting process. These figures show that natural convection begins to settle near the heat exchange wall separating HTF and NEPCM. Clockwise cells of ascending warm and descending cold liquid are created in the NEPCM liquid. The heated NEPCM liquid moves up toward the top part of slabs before impinging on the solid-liquid interface. At the beginning of the melting process (   0.01), the effect of the natural convection is generally weak and the main heat transfer mode is the thermal conduction. At   0.0502, the NEPCM liquid flow is monocellular and the flow intensity is substantially the same for both cases of volumetric fraction of nanoparticles (  min = 120). As time progresses, the effect of the natural convection becomes more important which affects the morphology of the solid-liquid interface. For both cases of volumetric fractions of nanoparticles, the displacement of the hot NEPCM liquid towards the top portion of the slabs causes an accelerated melting in this region. At   0.1075, the NEPCM liquid flow becomes bi-cellular and two clockwise cells are formed. With the progress of time (   0.1291, 0.1578 and 0.1893), the two clockwise cells formed 14

widen in volume and the convective flow appears more intensified for   0% (the maximum absolute value of the dimensionless streamline function is obtained for   0% (  min = 110, 70 and 50 at   0.1291, 0.1578 and 0.1893, respectively)). At   0.2668, the whole NEPCM has become liquid where the streamlines covered the entire slabs. It should be noted that for both cases of nanoparticles’ volumetric fractions, the two central cells previously formed combine into a single cell occupying the whole volume of the slab. The temperature distributions in the liquid phase of the NEPCM at the same times and conditions as previously, for   0% and   8% are illustrated in Figs. 7a-7g and 8a-8g, respectively. At the early stages of the melting process (   0.01), the heat transfer within NEPCM is controlled by thermal conduction and therefore the isotherms are almost in alignment with the y-axis. As time progresses, natural convection starts to affect the melting process and the NEPCM liquid begins to occupy a substantial portion of the slab. For both cases of nanoparticles’ volumetric fractions, the heated liquid phase in contact with the exchange wall moves up toward the upper part of the slab which causes deformation of the temperature contours. The maximum temperature is reached at the upper part of the NEPCM slab. These figures also show that the greatest temperature gradients are occurred close the heat exchange wall. It should be noted that the comparison between Figs. 7 and 8 shows that the maximum temperature is achieved for the case of   8%. It is also shown that the melting front ( nm  0 ) progresses more rapidly toward the right side for the case of   8% as shown in Figs. 8a-8g owing to the increase of the effective thermal conductivity of the NEPCM. 4.2 Parametric study In order to quantify the thermal performance of the NEPCM storage unit heated by laminar HTF flow, two kinds of storage efficiencies have been defined. The first concerns the thermal energy storage by sensible heat. It is expressed as the ratio between the sensible energy stored in the NEPCM storage unit and the maximum energy (latent and sensible energy) that could be stored. It is given as follows:

 sen 

(1   ( n c p ,n  1)) Ste 

w /2  d /2

X  w /2



H

Y 0

 nm dX dY

1    Ste(1   ( n c p ,n  1))

(40)

The second kind of efficiency concerns the energy storage by latent heat. It is defined as the ratio between the latent heat stored and the maximum energy that could be charged in the NEPCM. It is expressed as follows:

 lat 

f (1   ) 1    Ste(1   ( n c p ,n  1))

(41)

It should be noted that the latent storage efficiency reaches its maximum value when the entire NEPCM is melted (f=1). This maximum value of the latent storage efficiency, as can be concluded from Eq. (40), is lower than the value 1. It should be noted that the sensible storage efficiency achieves its maximum value when the liquid NEPCM is heated to the HTF inlet temperature. At this moment, the thermal equilibrium between the NEPCM and HTF is 15

established and the sum of the both storage efficiencies becomes equal to the value 1 (  sen   lat  1). 4.2.1 Effect of the aspect ratio of NEPCM slabs The effect of the NEPCM’s aspect ratio, A, on the latent storage efficiency for different volumetric fractions of nanoparticlesis illustrated in Fig. 9 for Re = 40 and Ra = 7.95x107. As can be seen, the latent storage efficiency quickly increases to achieve its maximum value and remains then constant. The maximum value of the latent storage efficiency is reached when all the NEPCM is completely melted. This maximum value of the latent storage efficiency is reduced with the dispersion of nanoparticles. Its value is about 0.874, 0.871 and 0.861 for 0%, 2% and 8%, respectively, for the all aspect ratios. The decrease of this maximum can be explained by the fact that the replacement of a volume fraction of PCM by nanoparticles leads to a reduction in the amount of PCM and hence to a decrease of the amount of energy stored in latent form. The effect of the NEPCM’s slab aspect ratioon the melting timefor different volumetric fractions of nanoparticlesis illustrated in Fig.10 for Re = 40 and Ra = 7.95x107. As can be seen, the melting time decreases with the increase of the aspect ratio. This is dueto the increase of the heat exchange surface between the NEPCM and HTF and thus the improvement of the heat transfer as the aspect ratio increases. This figure clearly shows that the melting time is reduced with the increase in the volumetric fraction of nanoparticles. The results also show that the effect of the suspension of nanoparticles in PCM slightly decreases with the increase in the aspect ratio of the NEPCM slabs. Indeed, when the volumetric fraction of nanoparticles ranges from 0 to 8%, the melting time is reduced by 14.25 %, 14 % and 13.91 %, for A= 0.5, 2 and 16, respectively. The effect of the aspect ratio of NEPCM slabs on the time wise variation of the sensible storage efficiency for various volumetric fractions of nanoparticles is shown in Fig.11 for Re = 40 and Ra = 7.95x107. This figure shows two distinct behavior stages for the sensible storage efficiency variation. In the first stage, the temperature of the molten NEPCM in contact with the heat exchange wall increases with the time progress, which linearly increases the amount of energy stored by sensible heat. Therefore, the sensible storage efficiency undergoes a linear increase over time. It should be noted that during this stage, the amount of energy stored in NEPCM by latent heat is predominant. In the second stage, the sensible storage efficiency undergoes a sudden increase to reach its maximum value and then remains constant. This sudden increase is due to the fact that during this stage, the NEPCM is almost completely melted, and heat transfer by natural convection in liquid NEPCM becomes more important which causes the quickly increase in the molten NEPCM temperature. The maximum value of the sensible storage efficiency is achieved when the NEPCM temperature becomes equal to the HTF inlet temperature, i.e. when the steady state is reached. For A = 0.5, the time required for the sensible storage efficiency to reach its maximum value is   0.84, 0.815 and 0.76 for   0%, 2% and 8%, respectively. For A = 2, the sensible storage efficiency reaches its maximum value at   0.53, 0.51 and 0.47 for   0%, 2% and 8%, respectively. For A = 16, the time corresponding to the thermal equilibrium between the NEPCM and HTF is about   0.283, 0.274 and 0.256 for   0%, 2% and 8%, respectively. It 16

is interesting to note that, for the same Rayleigh number, the maximum value of the sensible and latent efficiencies is a function of specific property of the NEPCM that only varies with the volumetric fraction of nanoparticles. Therefore, the sensible efficiency is independent on the aspect ratio and augments as the nanoparticles’ volumetric fraction increases. This result can be explained by the fact that the nanoparticles dispersed in the PCM stores energy only in sensible form which contributes to the total sensible energy storage. These results also show that the time corresponding to the thermal equilibrium between the NEPCM and HTF (charging time) decreases with the increase of the aspect ratio. This result is due to the fact that the increase of the aspect ratio allows to the increase in the heat exchange surface which enhances in turn the heat transfer between the HTF and NEPCM. 4.2.2 Effect of the Reynolds number The effect of the Reynolds number on the time wise variation of the latent storage efficiency for several volumetric fractions of nanoparticles is shown in Fig.12. The aspect ratio of the NEPCM slabs is A = 6 and the Rayleigh number is Ra = 7.95x107. The results reveal that the latent storage efficiency reaches its maximum value in a shorter time when the Reynolds number increases which causes the decrease in the time required for the complete melting. This maximum value of the latent storage efficiency slightly decreases with the increase of the volumetric fraction of NEPCM and its value is independent of the Reynolds number. The effect of the Reynolds number on the melting time is shown in Figs. 13a and 13b for A = 6 and 16, respectively. As concluded in the previous section, the melting time decreases with the increase in the Reynolds number. The decrease in the melting time is more pronounced for larger aspect ratio and volumetric fraction of nanoparticles. It is also interesting to note that a very lower reduction in the melting time is obtained for a Reynolds number exceeding the value 1000 (Re>1000). Consequently, it is recommended to avoid the use of very high HTF flow rate (Re>1000) during the storage process, so as not wasting too much pumping energy. The effect of the Reynolds numberon the time wise variation of the sensible storage efficiency, for different volumetric fractions of nanoparticles, is illustrated in Fig.14 for A= 6 and Ra = 7.95x107. As can be seen in this figure the sensible storage efficiency increases nonlinearly to reach its maximum value when the thermal equilibrium is established between the HTF and NEPCM. As it is expected, such values are independent of the Reynolds number. This figure also shows that the charging time of the storage unit decreases with the increase of the volumetric fraction of nanoparticles. Therefore, the sensible heat is early stored in NEPCM for a high nanoparticles’ volumetric fraction. It can be also concluded, from this figure, that the increase in the Reynolds number leads to the decrease in the charging time. Indeed, the increase in the Reynolds number intensifies the heat transfer rate between the HTF and NEPCM through the exchange surface. The effect of the Reynolds number on the time evolution of the dimensionless heat transfer rate under different volumetric fractions of nanoparticles is displayed in Fig.15 for A = 6 and Ra = 7.95x107. The heat transfer rate behavior shows three different stages. In the first stage, the heat transfer rate increases until achieve its maximum value. This maximum value increases significantly with the increase in the Reynolds number and slightly augments with the increase in the nanoparticles’ volumetric fraction. The heat transfer is mainly dominated 17

by conduction mode during this stage. In the second stage, the heat transfer rate undergoes a sudden decrease in the early moments, and then increases to reach its second maximum value (for Re = 40 and 100) or takes a virtually constant value (for Re = 20). The decrease in the heat transfer rate at the early moments is due to the molten thermal resistance formed near the heat exchange wall. However, the second increase in heat transfer rate is explained by the natural convection which starts to develop in liquid NEPCM and therefore intensifies NEPCM melting process. It is interesting to note that during this stage, the increase in the Reynolds number and nanoparticles’ volumetric fraction causes higher heat transfer rate between HTF and NEPCM. The third stage of heat transfer rate begins when the NEPCM melting process approaches to its end. In this stage, the heat transfer rate decreases gradually and achieves the zero value when the thermal equilibrium is achieved between NEPCM and HTF. 4.2.3 Effect of the Rayleigh number The effect of the Rayleigh number on the time wise variation of the latent storage efficiency, for different nanoparticles’ volumetric fractions, with A = 6 and Re = 40 is displayed in Fig.16. The latent storage efficiency undergoes a quickly increase to reach its maximum value when the NEPCM is completely melted. This maximum value of the latent efficiency depends on both the Rayleigh number (Stefan number) and the volumetric fraction of the nanoparticles as given in Eq. (41). It decreases with the increase in the Rayleigh number. In fact, when this latest is augmented, the Stefan number also increases via the HTF inlet temperature which leads to a reduction in the maximum value of the latent efficiency according to its definition. Physically this means that the maximum amount of energy that could be stored by sensible form increases with the increase of the Rayleigh number, and because the amount of latent heat stored is the same for different Rayleigh numbers (constant nanoparticles’ volumetric fraction and hence a constant mass of PCM), the maximum value of the latent storage efficiency decreases as the Rayleigh number is augmented. It can be also shown from this figure that the increase in the volumetric fraction of nanoparticles dispersed in pure PCM leads to the decrease in the maximum value of the latent storage efficiency, owing to the decrease in the PCM mass. The results also show that for Ra = 7.95x107, the maximum value of the latent storage efficiency is achieved at   0.2651, 0.256 and 0.227 for   0%, 2% and 8%, respectively. For Ra = 1.4x108, the maximum value of the latent storage efficiency is achieved at   0.14, 0.136 and 0.124 for   0%, 2% and 8%, respectively. For Ra = 2.63x108, the maximum amount of the latent heat is stored in the NEPCM at   0.0722, 0.0705 and 0.065 for   0%, 2% and 8%, respectively. From these results, it can be concluded that the melting time decreases with the increase in both the Rayleigh number and volumetric fraction of nanoparticles. Indeed, this is due to the fact that the increase of the Rayleigh number leads to the increase in the temperature difference between the HTF inlet and the melting point, while the increase in the volumetric fraction of nanoparticles enhances the PCM thermal conductivity and reduces the amount of base PCM. The effect of the Rayleigh number on the time wise variation of the sensible storage efficiency, for several volumetric fractions of nanoparticles, with A = 6 and Re = 40 is illustrated in Fig.17. As can be seen, the sensible storage efficiency increases to reach its 18

maximum value and then remains constant when the thermal equilibrium is achieved. In contrast to the results of the effects of the aspect ratio and Reynolds number presented above, the maximum value of the sensible storage efficiency depends on the Rayleigh number. This is obvious because the Rayleigh number represents the HTF inlet temperature on which the sensible energy storage also depends. Because the Rayleigh number increases with increasing the inlet HTF temperature, a large amount of sensible heat is stored in the NEPCM for a high Rayleigh number which explains the augmentation in the maximum value of the sensible storage efficiency with the Rayleigh number. It can be also concluded from these results that for the various volumetric fractions of nanoparticles, the charging time of the storage unit (time required by the sensible efficiency to reach its maximum value) decreases significantly with the increase in the Rayleigh number. Such a result is due to the fact that the increase in the Rayleigh number (which also represents an increasing in the HTF inlet temperature) leads to the increase of the temperature difference between the HTF inlet temperature and the melting point. The results also show that the charging time also decreases with the suspension of nanoparticles in PCM because of the enhancement of the heat transfer and the reduction in the amount of base PCM as the volumetric fraction of nanoparticles is augmented. 4.2.4 Correlation for the melting time Due to the significant amount of energy stored in the LHSU by latent heat, the melting time is considered to be the main outcome investigated by researchers. This parameter is of great interest for the design optimization of the LHSU. For this reason, the development of correlations that can estimate the melting time is highly requested. Therefore, based on the simulation results, a combined correlation for the melting time can be developed. The obtained correlation is given as a function of the volumetric fraction of nanoparticles, ϕ, the aspect ratio of the NEPCM slabs, A, the Reynolds number, Re, and the Rayleigh number, Ra. It is expressed as follows:

 melt 

 0  1 A  1    2 2   1 A 1   2 A  3    3 A2   4  2   2 A

(42)

The variables α1, α2, α3, β1, β2,β3,β4, γ1and γ2are formulated as a function of the Reynolds and Rayleigh numbers. The value of  0 is provided using the following expression: Ra Re     7  6.91265 10 3282 

 0  0.14524  1.5821exp 

(43a)

The values of the variables α1, α2, α3, β1, β2,β3,β4, γ1and γ2are given by using the following formula: (αi, βi, γi) 

a0  a1 Ra  b1 Re b2 Re2  b3 Re3  c1 Ra Re 1  a2 Ra  b4 Re a3 Ra 2  b5 Re2  a4 Ra 3  c2 Ra Re

(43b)

wherein the values of a0 , a1 , a2 , a3 , a4 , b1 , b2 , b3 , b4 , b5 , c1 , c2 are summarized in Tables 5, 6 and 7 for the calculation of variables group (α1, α2, α3), (β1, β2, β3) and (γ1, γ2), respectively.

19

The developed correlation in expression (42) is valid for the range of control parameters studied in the previous sections which are 0% ≤ ϕ ≤ 8%, 0.5 ≤ A ≤ 16, 20 ≤ Re ≤ 1000 and 7.95x107 ≤ Ra ≤ 2.63x108 (0.1444≤ Ste ≤ 0.47756). The melting times obtained by the numerical simulations are compared with those calculated by using the developed correlation. Figs. 18(a)-(d) show, for Re = 40, the effect of the aspect ratio, A, and Rayleigh number, Ra, on the melting time of the NEPCM for ϕ = 0%, 2%, 4% and 8%, respectively. For Re = 100, the effect of the aspect ratio, A, and Rayleigh number, Ra, on the melting time for ϕ = 0%, 2%, 4% and 8% is displayed in Figs. 19(a)-(d), respectively. Lastly, the effect of the aspect ratio and Rayleigh number on the melting time is shown in Figs. 20(a)-(d) for Re = 1000. From these figures, it can be clearly seen that the developed correlation estimates, with high precision,the time required for the complete melting of NEPCM inside the latent heat storage unit. In addition, the average relative error between the melting time obtained by the numerical simulations and that calculated from the developed correlation is 3.55 %. 5

Conclusion

In the present work, the melting of nanoparticle-enhanced phase change material in a rectangular latent heat storage unit (LHSU) heated by laminar heat transfer fluid (HTF) flow has been studied numerically. A mathematical model based on the conservation equations of mass, momentum and energy has been elaborated and validated to investigate the effect of the volumetric fraction of nanoparticles dispersed in pure phase change material, the aspect ratio of the NEPCM enclosures, the Reynolds number and the Rayleigh number on the flow characteristics and thermal performance of the LHSU. Correlation for the melting time of the NEPCM encompassing the volumetric fraction of nanoparticles, aspect ratio of the NEPCM slabs, Reynolds and Rayleigh numbers has been developed. Based on the findings of the present computational study, the following conclusions are drawn: 

 

  

The dispersion of high conductivity nanoparticles in a pure phase change material causes the decrease in the melting time and this effect rises with the increase in the volumetric fraction of nanoparticles (ϕ ≤ 8%); The increase in the Rayleigh number strongly decreases both the melting and charging times, and enhances the energy stored by sensible heat; The increase of the aspect ratio of the NEPCM slabs increases the heat exchange surface and therefore causes reduction in the time required for both sensible and latent heat storage; The increase in the Reynolds number (Re ≤ 1000) significantly decreases the melting time and the charging time ( the time required for complete sensible heat storage); The increase in the volumetric fraction of nanoparticles improves the sensible heat storage and slightly affects the latent heat storage; The comparison between the melting time calculated from the proposed correlation and that obtained by the numerical investigations showed a good agreement (with an average relative error of 3.55 %) and thus, the developed correlation can be used to estimate the time required for the complete melting with high accuracy. 20

References [1] T. Li, Y.V. Rogovchenko, Asymptotic Behavior of Higher-Order Quasilinear Neutral Differential Equations, Abstract and Applied Analysis, 2014 (2014) 11. [2] J.A.K. Xinglin Tong , M. RuhulAmin, Enhancement of heat transfer by inserting a metal matrix into a phase change material. [3] H. Ait Adine, H. El Qarnia, Numerical analysis of the thermal behaviour of a shell-andtube heat storage unit using phase change materials, Applied Mathematical Modelling, 33 (2009) 2132-2144. [4] Y. Tian, C.Y. Zhao, Thermal and exergetic analysis of Metal Foam-enhanced Cascaded Thermal Energy Storage (MF-CTES), International Journal of Heat and Mass Transfer, 58 (2013) 86-96. [5] J.M. Marín, B. Zalba, L.F. Cabeza, H. Mehling, Improvement of a thermal energy storage using plates with paraffin–graphite composite, International Journal of Heat and Mass Transfer, 48 (2005) 2561-2570. [6] N.S. Dhaidan, Nanostructures assisted melting of phase change materials in various cavities, Applied Thermal Engineering, 111 (2017) 193-212. [7] J.M. Khodadadi, S.F. Hosseinizadeh, Nanoparticle-enhanced phase change materials (NEPCM) with great potential for improved thermal energy storage, International Communications in Heat and Mass Transfer, 34 (2007) 534-543. [8] C.J. Ho, J.Y. Gao, An experimental study on melting heat transfer of paraffin dispersed with Al2O3 nanoparticles in a vertical enclosure, International Journal of Heat and Mass Transfer, 62 (2013) 2-8. [9] S. Sebti, M. Mastiani, H. Mirzaei, A. Dadvand, S. Kashani, S. Hosseini, Numerical study of the melting of nano-enhanced phase change material in a square cavity, Journal of Zhejiang University SCIENCE A, 14 (2013) 307-316. [10] R. Hossain, S. Mahmud, A. Dutta, I. Pop, Energy storage system based on nanoparticleenhanced phase change material inside porous medium, International Journal of Thermal Sciences, 91 (2015) 49-58. [11] S.Z. Shuja, B.S. Yilbas, M.M. Shaukat, Melting enhancement of a phase change material with presence of a metallic mesh, Applied Thermal Engineering, 79 (2015) 163-173. [12] A. Sciacovelli, F. Colella, V. Verda, Melting of PCM in a thermal energy storage unit: Numerical investigation and effect of nanoparticle enhancement, International Journal of Energy Research, 37 (2013) 1610-1623. [13] S. Kashani, A.A. Ranjbar, M.M. Madani, M. Mastiani, H. Jalaly, Numerical study of solidification of a nano-enhanced phase change material (NEPCM) in a thermal storage system, Journal of Applied Mechanics and Technical Physics, 54 (2013) 702-712. [14] S.F. Hosseinizadeh, A.A.R. Darzi, F.L. Tan, Numerical investigations of unconstrained melting of nano-enhanced phase change material (NEPCM) inside a spherical container, International Journal of Thermal Sciences, 51 (2012) 77-83. [15] M. Jourabian, M. Farhadi, A.A. Rabienataj Darzi, Outward melting of ice enhanced by Cu nanoparticles inside cylindrical horizontal annulus: Lattice Boltzmann approach, Applied Mathematical Modelling, 37 (2013) 8813-8825. [16] A.A. Ranjbar, S. Kashani, S.F. Hosseinizadeh, M. Ghanbarpour, Numerical heat transfer studies of a latent heat storage system containing nano-enhanced phase change material, Thermal Science, 15 (2011) 169-181. [17] R. Elbahjaoui, H. El Qarnia, M. El Ganaoui, Melting of nanoparticle-enhanced phase change material inside an enclosure heated by laminar heat transfer fluid flow, The European Physical Journal Applied Physics, 74 (2016) 24616. 21

[18] A.A. Valan, A.P. Sasmito, A.S. Mujumdar, Numerical performance study of paraffin wax dispersed with alumina in a concentric pipe latent heat storage system, Thermal Science, 17 (2013) 419-430. [19] J.M. Mahdi, E.C. Nsofor, Solidification of a PCM with nanoparticles in triplex-tube thermal energy storage system, Applied Thermal Engineering, 108 (2016) 596-604. [20] L. Fan, J.M. Khodadadi, An experimental investigation of enhanced thermal conductivity and expedited unidirectional freezing of cyclohexane-based nanoparticle suspensions utilized as nano-enhanced phase change materials (NePCM), International Journal of Thermal Sciences, 62 (2012) 120-126. [21] N. Das, Y. Takata, M. Kohno, S. Harish, Melting of graphene based phase change nanocomposites in vertical latent heat thermal energy storage unit, Applied Thermal Engineering, 107 (2016) 101-113. [22] S. Lokesh, P. Murugan, A. Sathishkumar, V. Kumaresan, R. Velraj, Melting/solidification characteristics of paraffin based nanocomposite for thermal energy storage applications, Thermal Science, (2015) 170-170. [23] M. Bechiri, K. Mansouri, Analytical study of heat generation effects on melting and solidification of nano-enhanced PCM inside a horizontal cylindrical enclosure, Applied Thermal Engineering, 104 (2016) 779-790. [24] V.R. Voller, C. Prakash, A fixed grid numerical modelling methodology for convectiondiffusion mushy region phase-change problems, International Journal of Heat and Mass Transfer, 30 (1987) 1709-1719. [25] J. Khodadadi, Y. Zhang, Effects of buoyancy-driven convection on melting within spherical containers, International Journal of Heat and Mass Transfer, 44 (2001) 1605-1618. [26] N. Wakao, S. Kaguei, Heat and Mass Transfer in Packed BedsGordon and Breach Science Publications, New York, (1982). [27] A.D. Brent, V.R. Voller, K.J. Reid, Enthalpy-porosity technique for modeling convection-diffusion phase change: Application to the melting of a pure metal, Numerical Heat Transfer, 13 (1988) 297-318. [28] C. Gau, R. Viskanta, Melting and solidification of a metal system in a rectangular cavity, International Journal of Heat and Mass Transfer, 27 (1984) 113-123. [29] W.M. Kays, M.E. Crawford, Convective Heat and Mass Transfer, McGraw-Hill, New York, 1970.

22

Table 1: Values of the fixed control parameters Control parameters Value

w

Prf

Prm

f

n

cp,n

kn

k m,s

kf

n

dn

0.085

3.77

19.87

1.3

4.49

0.305

150

1.5

2.666

0.02

6.37×10-7

Table 2: Effect of the grid size on the melting fraction at   0.04972 (34.66 min) with ϕ = 2 %, A = 6, Re = 40 and Ra = 7.95x107(Ste = 0.1444). Grid size, M×N 20 x 30 40 x 90 45 x 110

f 0. 2802 0. 3003 0. 2996

Deviation (%) 7.17 0.23

Table 3: Effect of the time step on the melting fraction at   0.04972 (34.66 min) with ϕ = 2%, A = 6, Re = 40 and Ra = 7.95x107(Ste = 0.1444). Time step, ∆τ 9.56x10-4 (40s) 4.78x10-4 (20s) 2.39x10-4 (10s) 1.195x10-4 (5s)

f 0.3056 0.3018 0. 3003 0.3004

Deviation (%) 1.24 0.497 0.033

Table 4: Thermal properties of HTF, PCM and alumina nanoparticles Property 3 Density (kg/ m ) Specific heat (J/ kg K) Thermal conductivity (W/ m K) Dynamic viscosity (kg/ m s) Latent heat (kJ/ kg) Melting temperature (° C) Thermal expansion (1/K)

HTF (water) 989 4180 0.64 577x10-6 -

23

paraffin Wax 802 2510 0.36(s), 0.24(l) 1.3x10-3 226 47 -4 5×10

Alumina (Al2O3) 3600 765 36 -5 1×10

Table 5: values of a0,a1,a2,a3,a4,b1,b2,b2,b3,b4,b5,c1 and c2 for the calculation of the variables α1, α2 and α3.

a0

α1 -175.6

α2 -2202.23

α3 -42.272

a1

9.2828 x10-7

1.31 x 10-5

6.248 x 10-7

a2

-1.049 x10-5

-3.802 x 10-5

-1.52 x 10-5

a3

4.76 x10-14

1.053 x 10-13

3.856 x 10-14

a4

0

0

0

b1

-5.4553

70.624

2.3713

b2

8.955 x10-5

0.08826

0.00203

b3

0

0

0

b4

-16.0012

2.45818

13.514

b5

6.826 x 10

-4

c1

2.2654 x 10

c2

7.02 x 10-8

-8

0.04734 -1.9911 x 10

0.02335 -6

-6.458 x 10-7

24

-4.465 x 10-8 -3.1756 x 10-7

Table 6: values of a0,a1,a2,a3,a4,b1,b2,b2,b3,b4,b5,c1 and c2 for the calculation of the variables β1, β2, β3and β4.

a0

β1 2116.24

β2 -23.396

a1

4.768 ×10-6

2.9744 ×10-7 -5

β4 -218439.28

0.00119

0.00147 -5.46 ×10-5

-2.6515 ×10

a3

1.783 ×10-13

-2.8725 ×10-16

6.64 ×10-13

3.866 ×10-13

a4

0

2.0208 ×10-24

0

0

b1

-184.155

0.01028

4912.235

-2646.083

0.54895

0.00339

0

0

568.685

-160.763

0.15781

-1.0163 ×10

b3

0

2.5895 ×10

b4

195.331

0.00141

b5

0.04929

-6.154 ×10 -8

c1

3.411 ×10

c2

-1.7781 ×10

-6

-5

-9

-7

-6.445 ×10

-5

a2

b2

-2.682 ×10

-9

β3 -69207.85

0.2816

0.00157

0

-6.5224 ×10

0

-7.73 ×10

25

-6

-5

3.01265 ×10-5 1.2031 ×10

-6

Table 7: values of a0,a1,a2,a3,a4,b1,b2,b2,b3,b4,b5,c1 and c2 for the calculation of the variables γ1 and γ2. γ1

γ2

a0

-347.748

6944.872

a1

2.422 ×10-6

-9.3331 ×10-5

a2

-1.1941 ×10-5

2.0562 ×10-5

a3

2.2613 ×10-13

-2.0651 ×10-13

a4

0

0

b1

46.5323

33.14678

b2

-0.00242

-0.0942

b3

0

0

b4

-73.1915

-39.7328

b5

-0.0013

-0.00471

c1

-3.2376 ×10-7

-9.7479 ×10-8

c2

4.83983 ×10-7

2.937 ×10-7

26

Fig.1a Latent heat storage unit (LHSU)

Fig. 1b Computationaldomain

27

2 min

6 min

10 min

17 min





y (m)







Current model Expr. Gau and Viskanta Brent et al.















x (m) Fig.2 Melting front at different times: (blue line) current model, (black dashed line) experimental results of Gau and Viskanta [28] and (red dashed line) numerical results of Brent et al. [27].

28

 

Current model    Current model    Current model    Khodadadi and Hosseinizadeh [7]    Khodadadi and Hosseinizadeh [7]    Khodadadi and Hosseinizadeh [7]   



Liquid fraction

       















Time (s)

Fig. 3 Instantaneous variation of the liquid fraction of NEPCM : comparison between the current results and those obtained by Khodadadi and Hosseinizadeh [7].

29

     

Nu

 

Isothermal wall Constant heat flux

   

8.23



7.54

 





 y (m)



Fig.4 Local Nusselt number (H = 20m)

30



  

  

c

e



 0

0

-8 -40

-5



-5

-120

0

-8

-5

0

-8

-8





  

d

-8

0

-4 0

-1 2

b

-15

  

  

a

0 -8 0

-12 -40











-125

-80 -8

Y

Y

Y

Y

-5

Y

-95

-15 -40













 -70



-1 1



0





   

    X

X

min = - 5 ;

  

    X

min = -120 ; max = 0

max = 0



min = -125 ; max = 0





    X

min = -110 ; max = 0



    X

min = -70 ; max = 0

  

f

g



-8

-15



-30

Y

Y

0

-40

-5



-45





-20



-15 -5

-50



0





    X

min = -50 ; max = 0



    X

min = -45 ; max = 0

Fig. 5 Streamlines for   0% at (a)   0.01, (b)0.0502, (c) 0.1075, (d) 0.1291, (e) 0.1578, (f) 0.1893 and (g) 0.2668 with A = 6, Re= 40 and Ra =7.95×107(Ste = 0.1444).

31

   -40 -80

b

  

  

c

-8

e

-5



-5

0

-8



-15



-5 -8





  

d

-120

0

-12

  

a

0

0

-40 -80









 -12

0

-80

Y

Y

Y

Y

Y

-40

-125









 -95

-15

-40 -5





   



X

    X

max = 0

min = - 5 ;



    X

min = -120 ; max = 0

  





min = -125 ; max = 0





    X

min = -95 ; max = 0



-60



    X

min = -60 ; max = 0

  

f

g



-35





-40

-35



-20

0

-30



Y

Y

-5

-8

-15





 -8 -5

-15

0

-35



    X

min = -35 ; max = 0



    X

min = -40 ; max = 0

Fig. 6 Streamlines for   8% at (a)   0.01, (b)0.0502, (c) 0.1075, (d) 0.1291, (e) 0.1578, (f) 0.1893 and (g) 0.2668 with A = 6, Re = 40 and Ra = 7.95x107 (Ste = 0.1444).

32

  

  

a

  

  

  





















Y

Y

e

Y

d

Y

c

Y

b























   



   

X   

  









Y

g

Y

f











    X



X



    X



    X



    X

           

    X

Fig.7 Isotherms for   0% at (a)   0.01, (b)0.0502, (c) 0.1075, (d) 0.1291, (e) 0.1578, (f) 0.1893 and (g) 0.2668 with A = 6, Re = 40 and Ra =7.95x107 (Ste = 0.1444)

33

  

  

  

b

  

  





















Y

Y

e

Y

d

Y

c

Y

a























   



   

X   

  









Y

g

Y

f











    X



X



    X



    X



    X

           

    X

Fig. 8 Isotherms for   8% at (a)   0.01, (b)0.0502, (c) 0.1075, (d) 0.1291, (e) 0.1578, (f) 0.1893 and (g) 0.2668 with A = 6, Re = 40 and Ra = 7.95×107(Ste = 0.1444)

34







A = 16 



 lat

A=2 



A = 0.5



        





















Fig. 9 Effect of the aspect ratio, A, on the time wise variation of the latent storage efficiency for various volumetric fractions of nanoparticles with Re = 40 and Ra = 7.95×107

35

  

        



 melt

      



















A Fig. 10 Effect of the aspect ratio, A, of the NEPCM’s slabs on the melting time for various volumetric fractions of nanoparticles with Re = 40 and Ra =7.95×107

36







A = 16

 sen



A=2 



A = 0.5

        

























Fig. 11 Effect of the aspect ratio, A, on the time wise variation of the sensible storage efficiency for various volumetric fractions of nanoparticles with Re = 40 and Ra = 7.95×107

37





Re = 1000 



Re = 40

 lat





Re = 20





        























Fig. 12 Effect of the Reynolds number on the time wise variation of the latent storage efficiency for various volumetric fractions of nanoparticles with A = 6 and Ra=7.95x107

38

(a)

(b) 





        



        



melt

melt



























Re











Re

Fig. 13 Effect of the Reynolds number on the melting time for various volumetric fractions of nanoparticles for (a) A = 6 and (b) A = 16 with Ra = 7.95x107

39









Re = 1000

 sen





Re = 40



Re = 20 



         













Fig. 14 Effect of the Reynolds number on the time wise variation of the sensible storage efficiency for various volumetric fractions of nanoparticles with A= 6 and Ra=7.95×107

40





  

q



Re = 1000 

Re = 40



Re = 20 













Fig. 15 Effect of the Reynolds number on the time wise variation of the dimensionless heat transfer rate for various volumetric fractions of nanoparticles with A= 6 and Ra=7.95x107

41





        





Ra = 2.63 108

 lat



Ra = 1.4 10



8



Ra = 7.95 10

7























Fig.16 Effect of the Rayleigh number on the time wise variation of the latent storage efficiency for various volumetric fractions of nanoparticles with A = 6 and Re = 40

42





Ra = 2.63 10



        

8

 sen





Ra = 1.4 10

8







Ra = 7.95 10











7









Fig.17 Effect of the Rayleigh number on the time wise variation of the sensible storage efficiency for various volumetric fractions of nanoparticles with A = 6 andRe = 40

43

a



b

Numerical simulation Developed Correlation

Re = 40



Numerical simulation Developed Correlation

Re = 40



  

Ra = 7.95 10



  

7





m

m

Ra =7.95 10



Ra = 1.4 10



8

Ra =1.4 10



8



Ra = 2.63 10

8





7

Ra =2.63 10



























8





A



d

Numerical simulation Developed Correlation

Re = 40









Numerical simulation Developed Correlation

Re = 40



  



  



Ra = 7.94 107

m

Ra = 7.95 107

m





c



Ra = 1.4 10











8

Ra = 2.63 10







A



Ra = 1.4 10

8

8















A



Ra = 2.63 10







8











A

Fig.18 Comparison between the melting time calculated using the developed correlation and numerical predictions for (a) ϕ=0% (b) 2% (c) 4% and (d) 8% with Re = 40

44





a



b

Numerical simulation Developed Correlation

Numerical simulation Developed Correlation



Re = 100

Re = 100

  





Ra = 7.95 10

Ra = 7.95 10

7

7



m



m

  





Ra = 1.4 10



8

Ra = 1.4 10



8



Ra = 2.63 10



8



Ra = 2.63 10

8



































A

c











d

Numerical simulation Developed Correlation



Numerical simulation Developed Correlation

Re = 100

Re = 100

  



Ra = 7.95 10



  



7

Ra = 7.95 10

7



m



m



A





Ra = 1.4 10



8

Ra = 1.4 10

8







Ra = 2.63 108

Ra = 2.63 10

8























A





















A

Fig.19 Comparison between the melting time calculated using the developed correlation and numerical predictions for (a) ϕ=0% (b) 2% (c) 4% and (d) 8% with Re = 100

45





a

b

Numerical simulation Developed Correlation 

Numerical simulation Developed Correlation



Re = 1000

Re = 1000

  

  





Ra = 7.95 140

7

Ra = 7.95 10



Ra = 1.4 140



8



Ra = 1.4 10





8







Ra = 2.63 140 

7



m

m







8



Ra = 2.63 10 



















8



A











A





c

d

Numerical simulation Developed Correlation 

Numerical simulation Developed Correlation 

Re = 1000

Re = 1000

   

Ra = 7.95 10

   

7

Ra = 7.95 107 

m

m





Ra = 1.4 10 

Ra = 1.4 10





8



Ra = 2.63 10

8

Ra = 2.63 10







8

8





















A





















A

Fig. 20 Comparison between the melting time calculated using the developed correlation and numerical predictions for (a) ϕ = 0% (b) 2% (c) 4% and (d) 8% with Re = 1000

46

Highlights 

Melting of nanoparticle-enhanced phase change in a rectangular latent heat storage unit is investigated numerically



Melting time decreases with increasing volumetric fraction of nanoparticles



Melting time decreases with increasing Rayleigh number



Melting time decreases with increasing aspect ratio of NEPCM slabs



Melting time decreases with increasing Reynolds number

47