- Email: [email protected]
A numerical investigation of melting phase change process via the enthalpy-porosity approach: Application to hydrated salts
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Contents lists available at ScienceDirect
International Communications in Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ichmt
A numerical investigation of melting phase change process via the enthalpyporosity approach: Application to hydrated salts
MARK
Zohir Younsia,b,c, Hassane Najib,c,⁎ a b c
FUPL, Hautes Etudes d'Ingénieur (HEI), LGCgE (EA 4515), 13 rue de Toul, F-59000 Lille, France Univ. Artois, Laboratoire Génie Civil & géo-Environement (LGCgE EA 4515), Technoparc Futura, F-62400 Béthune, France Université Lille Nord de France, LGCgE (EA 4515), Lille F-59000, France
A R T I C L E I N F O
A B S T R A C T
Keywords: Phase change material Natural convection Numerical simulation Finite volume method Enthalpy-porosity formulation Melting Fusion
The convection dominated melting of hydrated salts (PCM 27) is investigated numerically using a pressure based finite volume method (FVM) with an enthalpy porosity technique, which allows a fixed-grid solution of the coupled momentum and energy equations, and trace the solid-liquid interface without resorting to other equations or transformations. This investigation addresses effects of Rayleigh number, the aspect ratio and Prandtl number on the effectiveness of melting process and solid-liquid interface. In addition, the average Nusselt number is set to complete the parametric study conducted. Through the results achieved, the enthalpyporosity has been found to converge rapidly while producing accurate results for both the position and morphology of the melt front at different times. The results are presented in terms of temperature profiles, streamlines, isotherms, moving interface position, liquid fraction and average Nusselt number. From this survey, it appears that our numerical findings support the experimental observations. These results indicated that, in terms of efficiency of the melting process, the best configuration of the rectangular container is the cavity having an aspect ratio of 0.19 (H = 21 cm ; W = 4 cm). Therewith, it is found that the time required to melt the material is proportional to the Rayleigh number and the convection heat transfer. Finally, this study demonstrates the feasibility for simulating macro-encapsulated hydrated salts via the enthalpy-porosity-based approach.
1. Introduction Storage of thermal energy is essential for reducing the continuous increase of the gap between energy demand and supply. There are three different ways one can store thermal energy in a storage material: form of sensible heat, latent heat and thermochemical energy. In the thermochemical storage system the corresponding storage materials used are still too expansive and their use requires a technical and economic feasibility study [1]. The sensible heat thermal storage (SHTS) system is very simple and the temperature of the storage material changes with the quantity of energy added or removed from the system. On the other hand, in a latent-heat thermal energy storage (LHTES) system, the density of energy storage is higher allowing the design of more compact systems with less heat losses and allowing energy storage at a nearly constant temperature [1,2]. Several studies have shown that using phase change materials (PCMs) as promising LHTES materials allows absorbing or releasing latent heat during a phase change from solidstate to liquid-state (melting process) or liquid-state to solid-state (freezing process). This system found several potential engineering ⁎
applications, such as thermal energy storage, thermal management, heating and cooling systems, electronic products, drying technology, waste heat recovery, solar air collectors, and solar cookers, and so on. One disadvantage of LHTES systems is the lower thermal conductivities of PCMs used, which slowed the melting/solidification rate. A faster fusion/solidification rate is highly sought after in the field of engineering [3]. Among the techniques for improving the thermal conductivity of PCMs and the heat transfer are microencapsulated PCM [4], porous metals [5], use of stationary highly conductive structures [6] and multiple PCMs [7]. In addition to these approaches, optimization of the conditioning system of PCMs along with the choice of type of PCM is the main interest of this research. These parameters modify the phase change process, as well as the heat transfer rate, which will be discussed in the present work. Since organic and inorganic PCMs are widely used in various engineering fields, they have been extensively studied in literature. For low-temperature residential applications, the hydrated salts (inorganic PCMs) are the most commonly used because of their low cost, low volumic expansion, stability of melting point and high volumetric
Corresponding author at: Univ. Artois, Laboratoire Génie Civil & géo-Environement (LGCgE EA 4515), Technoparc Futura, F-62400 Béthune, France. E-mail address: [email protected] (H. Naji).
http://dx.doi.org/10.1016/j.icheatmasstransfer.2017.05.012
0735-1933/ © 2017 Elsevier Ltd. All rights reserved.
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
ΔV x, y
Nomenclature A B C c f g H h k L W H Nu Pr p R Ra S Ste T T0 Tm TC TH t u v
matrix coefficients porosity function morphological constant (see Relationship (5)) specific heat (J/kg ∙K) liquid fraction gravity vector (m/s2) total enthalpy (J/m3) specific enthalpy (J/m3) thermal conductivity (W/m ∙K) latent heat of fusion (J/kg) length of rectangular enclosure in x-direction (m) length of rectangular enclosure in y-direction (m) Average Nusselt number Prandtl number pressure (Pa) aspect ratio (H/W) Rayleigh number source term Stefan number temperature (K) initial temperature (K) melting point of materials (K) isothermal right wall temperature (K) isothermal left wall temperature (K) time (s) velocity in the x direction (m/s) velocity in the y direction (m/s)
control volume Cartesian coordinates (m)
Greek symbols α β ε ∇ ν Δt Γ μ ρ φ ω
thermal diffusivity (m2/s) expansion coefficient (K− 1) arbitrary computational constant Nabla operator kinematic viscosity (m2/s) time step (s) diffusion coefficient dynamic viscosity (kg/m ∙s) density (kg/m3) generalized variable relaxation factor
Superscripts/subscripts c interface position k iteration level l liquid phase m melting P present node ref reference temperature value S solid phase W, E, P, N, S west, east, center, north and south nodes 0 initial value
of the cavity. Brent et al. [25] were among first authors to develop the enthalpy-porosity technique to trace the solid-liquid interface during the melting and solidification process. This approach has been validated by replicating the results of experimental work [24] using a 2D rectangular cavity. Shamsundar and Sparrow [26] proposed a multi-dimensional conduction phase change model for a solar energy storage system using the enthalpy-based method. Hamdan and Elwerr [27] have developed a mathematical approach to handle solid-liquid interface propagation during the two-dimensional melting process. Dhaidan et al. [28] have explored experimentally and numerically the melting of n-octadecane with CuO nanoparticle suspensions in a square enclosure. Their findings showed a significant effect of natural convection and an increase of the Rayleigh number on the shape of the solid-liquid interface. One of the most successful applications for PCM is thermal energy storage for air heating. Younsi et al. [29] have provided a new air heating system with a PCM energy store called Trombe wall improved, and allowing reducing energy consumption by 30% compared to conventional wall systems. As shown in Fig.1, the Trombe wall improved
energy storage density [8]. In this context, the calcium chloride hexahydrate CaCl2 ∙ 6H2O is one of potential materials for the low-temperature thermal storage, and thereby it has been explored extensively in the literature research [9–12]. Younsi et al. [12] have investigated experimentally the melting characteristics of salt hydrates in a rectangular cavity heated through its sides. Their results indicated that the melting efficiency has great potential for storing heat. Tyagi and Buddhi [13] have experimentally studied the thermal stability of CaCl2 ∙6H2O after any thermal cycling tests. They concluded that CaCl2 ∙ 6H2O is the material which has a promising potential for heating and cooling applications with small variations in the material latent heat of fusion during thermal cycling [8]. Macro-encapsulation of PCMs is a promising technic for storing the latent heat energy where the material is sealed in an enclosure which can be of any shape (rectangular, cylindrical, and spherical). It should be noted that the choice of the enclosure's geometry is a key parameter for improving the heat transfer rate. The rectangular and cylindrical containers are the most used due to their facility and low production cost [14]. To assess the potential use of macro-encapsulated PCMs in rectangular containers, many numerical, analytical and experimental studies have been conducted over the last decade [15–19]. The melting of PCMs aroused great interest because of its important role in the thermal energy storage. Several analytical, numerical and experimental studies were carried on phenomena associated with melting in different geometry configurations [20–23]. These works have concluded that the onset of melting is clearly dominated by diffusion, while buoyancy-driven convection has a paramount role in the melting process. However, little effort has been made on the optimization of the melting of a PCM packaged in a rectangular envelope and lay out vertically. Gau and Viskanta [24] have investigated experimentally solid-liquid interface motion and heat transfer during melting and solidification in a rectangular cavity filled with Gallium (Ga). They found that the heat transfer at the beginning of melting is dominated by conduction as time progresses, and that the convection develops on top
Fig. 1. Illustration of a latent heat thermal energy storage system [29].
13
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
system utilizes a thermal energy storage unit and an intelligently controlled air distribution unit. The principle of such a system can be outlined briefly as follows: the air between the glazing and the PCM is heated by greenhouse effect and allows charging the PCM. The PCM is discharged by passing air through the space between the wall and the insulation allowing air to be heated and providing thermal comfort for occupied spaces. Based on the review conducted, several approaches have been developed to handle the problem of the phase change, and many researchers have developed numerical models to assess the effectiveness of phase change materials packaged in rectangular containers. To go further in that sense, here, we aim to achieve the design of a latent heat thermal energy storage system by numerically addressing the problem of melting natural convection in a rectangular container heated from one of its sides. Such a problem is handled by using the enthalpy-porosity-based approach. To assess the accuracy of this formulation, main findings are compared with available results of the literature results [24,30–33] on ratios' aspect. In addition, simulations have been carried out to assess effects of relevant parameters (e.g. Rayleigh number, aspect ratio and Prandtl number) on the melting process. Likewise, the fusion time of the hydrated salts in a rectangular container is compared to these parameters. A complete physical overview of the dynamic and thermal fields in terms of streamlines, isotherms and the average Nusselt number along the hot wall are presented. In the light of the results achieved, it should be stated that the required time for good PCM melting inside the wall is highly controlled both by the container shape and heat flux imposed. The remainder of this paper is organized as follows. The problem description and mathematical formulation are introduced in Section 2. Subsequently, the computational methodology is briefly outlined in Section 3. In this section, the enthalpy FV-based approach with the discretization of the governing equation along with the Thomas algorithm is briefly reported. Afterwards, in Section 4, a validation of the adopted numerical approach via available results has been performed. Section 5 presents and comments our findings for a range of relevant parameters, such as Rayleigh number, aspect ratio and Prandtl number. Finally, the main conclusions are addressed in Section 6.
Table 1 Thermophysical properties of the PCM 27. Material
Solid phase Liquid phase
Thermophysical properties Density (kg/m3)
Specific heat (J/kg ∙K)
Thermal conductivity (W/ m ∙K)
Latent heat of fusion (kJ/ kg)
1710 1530
1751.5 2225
0.577 0.813
– 172.42
hydrated salt (PCM 27) having a melting temperature of 27 °C. Such a model has been studied experimentally by [34]. The horizontal walls are perfectly insulated (adiabatic) while the two left and right vertical walls are maintained isothermal at constant temperatures TC (cold) and TH (hot), respectively. The PCM in its fluid phase is assumed to be viscous, homogeneous; isotropic, incompressible and Newtonian. Its flow is two-dimensional, laminar and unsteady, without viscous dissipation or radiative heat transfer. In addition, the liquid motion induced by volumetric variation during melting is assumed to be neglected, thereby pointing that the solid density and that of the liquid are here equal. For the purpose of the numerical simulation, a no-slip condition holds at all the walls. Moreover, all thermo-physical properties of the considered PCM are assumed constant except the density variation in the buoyancy term which is modeled by the Boussinesq approximation for involving thermal buoyancy. These thermo-physical properties of these components are gathered in Table 1. Note that, in this analysis, it is assumed that the initial temperature of the PCM is equal to TC. Thereby, thermal energy will start penetrating through the left wall due to the temperature difference and will initiate the melting process. 2.2. Mathematical modeling Based on the aforementioned assumptions, the governing equations for the 2D transient laminar flow including heat transfer melting process and buoyancy-driven convection may thus be written as follows: Solid region:
2. Problem description and governing equations
∂T = αs ∇2 T ∂t
2.1. Geometry
(1)
Liquid region: The model enclosure is a 2-D square cavity with dimensions W × H on the x–y plane as depicted in Fig. 2. This enclosure is filled with a
→ ∇⋅ V = 0
(2) Fig. 2. Configuration under study (left) along with coordinates and boundary conditions (right).
14
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
→ → → → → ∂V + ( V ⋅∇) V = −∇ (p ρl ) + νl ∇2 V + β (T − Tm )→ g + S→ V ∂t
(3)
∂T + (V ⋅∇) T = αl ∇2 T + Sh ∂t
(4)
In Eqs. (3)–(4), source term S is, respectively given by:
→→ → L ∂f S V = −B (f ) V and Sh = − cl ∂t
h=
∫T
T
ρcdT
→ → → → → → g ∂V + ( V ⋅∇) V = −∇P + Pr∇2 V + Ra Pr θ → + BV ∂t͠ |g |
(11)
→ 1 ∂f ∂θ + ( V ⋅∇) θ = ∇2 θ − Ste ∂t͠ ∂t͠
(12)
→ 1 ∂f ∂θ + ( V ⋅∇) θ = α͠ ∇2 θ − Ste ∂t͠ ∂t͠
• left wall: • right wall: • bottom wall: • top wall: • solid walls:
(7)
Pr = ν α, Ra = gβH 3 (TH − Tm) Ste = cl (TH − Tm) Lf and R = W H αν , (14) The average number of Nusselt, which is a dimensionless measure of heat transfer through the hot wall can be computed as follows:
Nu = u = v = 0 and T = TH u = v = 0 and T = TC u = v = 0 and ∂ T/∂y = 0
1 R (Tm − TH )
H
∫ (∂T
∂x )x = 0 dy (15)
0
Prior to implementing the numerical algorithm, Eqs. (10)–(13) have been recast in the following general equation of a conservation equation for a general variable φ:
u = v = 0 and ∂ T/∂y = 0 u = v = 0 and T = Tm
→ ∂ (φ) + ∇⋅( V φ) = ∇⋅(Γφ ∇φ) + Sφ ∂t
(16)
In this equation, the expression of diffusion coefficients Γφ and the source term Sφ of each transported variable φ are gathered in Table 2.
→ → T(x, y, 0) = Tini = Tc and V = 0 at solid walls.
where TH > Tm > TC and s is the horizontal distance between the left isothermal wall and the melting front. Note that TC is kept constant in all simulations performed here. To cast the above equations in a dimensional form, we employ the following dimensionless variables:
3. Computational procedure Since the numerical method is adopted from the existing literature, we give here only a brief outline of the major features of the method used. The interested reader can refer for example to Ref. [41] for complete details on such a method. Integrating Eq. (16) over a control
t͠ = tαl H 2, (X , Y ) = (x H , y H ), (U , V ) = (uH αl , vH αl ), θ = (T − Tm) (TH − Tm), P = pH 2 ραl2 , α = αs αl, B∗ = BH 2 ραl
Table 2 Diffusion coefficient Γφ and source term Sφ for variable φ in the generic transport equation.
(8) Based on these dimensionless parameters, the aforementioned Eqs. (1)–(5) are then written as follows: Solid region
∂θ = α∇2 θ ∂t͠
θ|X = 0 = 1 θ|X = W/H = (TC − Tm)/(TH − Tm) ∂ θ/∂Y |Y = 0 = 0 ∂ θ/∂Y |Y = 1 = 0 → → V (U , V ) = 0
The dimensionless initial condition is given by: → → V = 0 , θ = 0 at t͠ = 0 , 0 < (X, Y) ≤ 1. It should be noted that such a problem involves the Prandtl number, the Rayleigh number, the Stefan number and the aspect ratio, defined by:
where ρfL is latent heat, which varies between zero for solid to ρL for liquid. Note that using this formulation, the absorption of latent heat during melting is isolated in the source term Sh. In the past, this method has been validated many times for convection-dominated melting problems in a rectangular container with different PCMs (Gallium, tin, n-Octadecane, caprillic acid) [25,42–44]. The required boundary and initial conditions for the considered problem are
• left wall (x = 0): • right wall (x = W): • bottom wall (y = 0): • top wall (y = H): • melting front (x = s): • initial conditions:
(13)
where α͠ = f + (1 − f ) α . In this way, in the liquid zone (θ > 0, f = 1), α͠ = 1, and in the solid zone (θ < 0, f = 0), α͠ = α . Note that the term (S͠ h = −Ste−1∂f ∂t͠ ) is the latent heat term (source) due to phase change. These equations are subjected to the following boundary conditions:
(6)
m
(10)
Note that Eqs. (9) and (12) can be written as follows:
(5)
where B(f) = C(1 − f)2/(f3 + ε) causes the gradual change in velocity from a finite value in the liquid to zero in the solid in the computational cells, f being the volume liquid fraction. Here, the multiplier C is the mushy zone constant (=109) reflecting the morphology of the melting/ solidification and ε is an arbitrary computational constant, set (= 3 × 10− 4 here) for avoiding a possible zero denominator (when → f = 0 corresponding to solid medium case); S → V is the momentum dissipation source term used to suppress the velocity in the mushy and solid zone, and Sh is the latent heat source term (caused by absorption or release). Gravity acts down (in the negative y-direction). As implied, an enthalpy-based method is adopted here for modeling the phase change problem [35–38]. Thereby, the energy equation follows Voller's formulation [39–41] for which the total enthalpy is split into sensible and latent enthalpy parts as
H (T ) = h (T ) + ρfL
→ ∇⋅ V = 0
(9)
Liquid region 15
φ
Γφ
Sφ
1 U V θ
0 Pr Pr f + (1 − f)αs/αl
0 −∇XP + BU Ra Pr θ − ∇YP + BV
−
∂f ∂t͠
Ste
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Ap φp =
Table 3 Thermophysical properties of Gallium.
∑
Ai φi + Sφ
i = W , E , S, N
cp (J/ kg ⋅ K)
k (W/ m ⋅ K)
Lf (kJ/kg)
ρ (kg/m )
υ (m /s)
β (K
302.93
381.5
32
80.16
6093
1.38 × 10− 5
1.2 × 10− 4
2
(17)
The iterative solution continues until convergence of the flow and energy fields at every time step is reached. Convergence is declared only when the three following conditions are met: 1) the liquid fraction field is stabilized, 2) the dimensionless residual for the enthalpy equation is less than 10− 5, and 3) the dimensionless residual for the mass conservation equation is less than 10− 4. These tolerances are chosen so that the final solution did not show perceptible changes. The pertinent feature of the present enthalpy approach for modeling the phase change problem is the source term Sh for the energy equation (Eq. (13)). It is useful to recall that this term keeps track of latent heat evolution, and its driving element is the local liquid fraction. This term is an indicator which is 1 in fully liquid regions, 0 in completely solid regions, while remaining within the interval [0, 1] close to the phase front. In addition, a time-varying distribution of the local liquid fraction is achieved iteratively from the solution of the energy equation via the following expressions:
−1
Tm (K)
3
)
f Pk + 1 = f Pk + ωΔtAPk hPk ρLΔV
(18)
with Fig. 3. Progress of the melting front at different times. Comparison with Refs. [24,30,32] at Pr = 0.0216, Ste = 0.0391, R = 2, and Ra = 2.1 × 105.
f Pk + 1 =
k+1 ⎧1 if f P > 1 ⎨ 0 if f k + 1 < 0 P ⎩
(19)
where ω is an under-relaxation factor. It is worth noting that, after the kth-iteration, the liquid fraction update is applied at every node to reflect the evolution of the phase change. 4. Validation of the enthalpy-based approach Prior to computations, checks have been conducted to validate our simulations compared to other results which have been reported in literature [24,30,31–33]. For this purpose, two benchmark cases have been carried out to demonstrate the approach herein considered. The first test case is one that had been handled by Gau et al. [24], Gong and Mujumdar [30], and Lacroix [32]. This is two-dimensional problem of melting of Gallium with the presence of natural convection in a rectangular enclosure at Pr = 0.0216, Ste = 0.0391, and Ra = 2.1 × 105. The dimensions of the cavity are H = 0.0445 m and W = 0.089 m (aspect ratio R = 2). The left hot wall and the right cold wall are maintained at temperatures, TH = 303.35 K and TC = 302.5 K, respectively. The horizontal walls are adiabatic. The relevant thermophysical properties of Gallium are listed in Table 3. Several numerical investigations were conducted to check the grid size and time step dependence results. A uniform grid size (50 × 50), time step (Δ t = 1 s), and under-relaxation factors for the momentum and energy equation (ωh = 0.6 and ωu = ωv = 0.8) were found sufficient to give accurate results. Several works were performed on the fusion of gallium and shows that most authors specify that the results are independent of the mesh. Viswanath [31] showed that the mesh refinement of 30 × 20–50 × 40 leads to a change of 3% of the final solution and concluded that this change is acceptable. The melting front is plotted at several times during the melting process, i.e., 5, 8 and 12 min (see Fig. 3), and compared favorably with the results of Gau and Viskanta [24], Gong and Mujumdar [30] and Lacroix [32]. In other words, it appears that experimental melt front positions are quite well reproduced with a slight discrepancy. Such a difference may be attributed to 1) the ignorance of the solid phase subcooling (2 °C under the melting temperature) by the current model, 2) the low wall temperature (8 °C above the melting) and the high metal conductivity (k = 32 W ∙ m− 1 ∙ K− 1), for which the metal sub-cooling seems to influence the experiment running (see [15])[, and 3]) the difficulty of maintaining vertical walls at a desired temperature during
Fig. 4. Progress of the melting front at different times. Comparison among results of Refs. [24,31,33] at Pr = 0.0216, Ste = 0.0391, R = 1.4 and Ra = 6.057 × 105. Table 4 Dimensionless parameters for the melting of PCM hydrated salts. R
Ra
Pr
Ste
cs/cl
ks/kl
0.12
1.34 × 108
274.7
0.2968
0.79
0.71
volume and estimation of each term of the integral equation in terms of discrete value on the nodal point lead to the production of a system of algebraic equations. To insure stability, a first order hybrid scheme [45] is, however, used for the discretization of the convection-diffusion terms in the momentum and energy equations. A fully implicit Euler method is used for the time-stepping procedure. Algebraic systems of equations originating from spatial and temporal discretizations are solved iteratively (throughout the physical domain considered) using the Tri-Diagonal Matrix Algorithm (by sweeping the cavity from left to right). On the other hand, the Semi-Implicit Algorithm for Pressure Linked Equations (SIMPLER) algorithm has been chosen to handle the pressure-velocity coupling. Note that the discretized equations for the dependent variable φ are written for each node P as:
16
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Fig. 5. Isotherms, streamlines and melting interface at different instants and various Rayleigh numbers with R = 0.12; a) Ra = 7.8 × 107; b) Ra = 1.34 × 108; and c) Ra = 1.98 × 108.
a) Ra =7.8×107 ; TH = 40 °C
b) Ra =1.34×108 ; TH = 50 °C
17 min
60 min
120 min 180 min 8 c) Ra = 1.98×10 ; TH = 60 °C
17
240 min
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
temperature-based method. At 8 min, the gap decreases between results, especially between our current results and the Gong and Mujumdar [30] experimental results. As time proceeds (at 12 min), our findings reproduce well those of Gau and Viskanta [24], except for the bottom region. It should be pointed out that, at the enclosure top end, all results tend to converge. The second benchmark test deals with the two-dimensional problem of the fusion of pure Gallium a rectangular cavity having an aspect ratio of 1.4. The hot and cold walls are isothermal (TH = 311.15 K and TC = 301.45 K), while horizontal walls remain adiabatic. In addition, the initial Gallium temperature was equal to TC. Fig. 4 gathers the comparisons concerning to the melting front progression at different times, viz. 5, 10 and 19 min. It can be seen that the comparison of the present results with those of the selected authors remains quite favorable. More precisely, at 5 min and 10 min, the agreement is excellent while at t = 19 min, some deviation is exhibited between all numerical simulations and experimental data. At this time, the present model simulations appear to be slightly offset as the melting front is behind that proposed by other solutions at the cavity top, while it is in front at the cavity bottom. To support this, it should be noted that Kim et al. [46] stated that the delayed heat-up affects the melting process and thereby can lead to a discrepancy between numerical and experimental results. Indeed, considering a delayed heat-up, they obtain a better agreement, especially at 19 min. To sum up, it seems that under these benchmark cases, the current model corroborates available results in literature. Thereby, it can be stated that our findings give confidence to the built code for the problem in this study.
Fig. 6. Temporal evolution of liquid fraction with different Ra (melting of PCM 27).
5. Results and discussion This section is dedicated to present and discuss the melting process in the enclosure considered in terms of influence of key parameters, such as the Rayleigh number, the aspect ratio and the Prandtl number on the dynamic and thermal fields, the melting fraction, and the average Nusselt number. In the following, some significant results showing effects of these parameters are comprehensively presented and discussed. To achieve this, a macro-encapsulated PCM 27 (hydrated salts) in a rectangular container is considered as a thermal energy storage system as shown in Fig. 2. The hexahydrate of calcium chloride (CaCl2 ∙6H2O) having a 27 °C melting temperature, is chosen as thermal latent energy storage material. Recall that the melting mode is investigated via the enthalpy-porosity approach. It should be pointed out that the main thermo-physical properties and dimensionless parameters of the sample are listed in Tables 1 and 4, respectively. In all cases, a 150 × 150 grid was used. Initially, the storage PCM is solid with a uniform temperature of 15 °C. Then, the temperature of the left vertical wall is suddenly raised to a prescribed temperature higher than the melting point, resulting in melting of hexahydrate of calcium. It should be noted that the efficiency of heat transfer with a PCM is thoroughly related to its rate melting. Therefore, the efficiency of the heat transfer and the rate at which the PCM melts greatly depends on the enclosure shape. Note that, various geometries containing the same PCM's amount were considered. These were either thick, thin, high or small. In this study, several efficiency indicators can be monitored, including mean and maximum cavity temperatures, total liquid fraction, rate at which PCM melts, and the Nusselt number. However, in the present study, we only retained the last 3 parameters, viz., the total liquid fraction, the PCM melting rate, and the Nusselt number.
Fig. 7. Temperature profiles at mid-height (y = H/2) of cavity for different Ra at 240 min.
Fig. 8. Temporal evolution of the average Nusselt number for various Ra numbers at R = 0.12.
the experiment. It should be noted that the sub-cooling is of approximately 2 °C in the solid phase. As for the discrepancy exhibited with Lacroix results, they can be imputed to the use of an enthalpy-porosity method with a fixed grid, instead of the temperature-based method involving a mobile grid to handle phase change effects. Despite this, it is observed that, at 5 min, the two approaches show similar results to the two experiments for the enclosure's bottom. Regarding the brick's top, the present model results seem to be in good agreement with those of Lacroix obtained with the
5.1. Influence of Raleigh number Recall that the Rayleigh number assesses the competition between convective and conduction heat transfer rates. To demonstrate the effect of this number, different values have been considered, viz., 7.8 × 107, 1.34 × 108 and 1.98 × 108, which correspond respectively 18
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Fig. 9. Isotherms and the melting interface at different times and several aspect ratios for Ra = 1.34 × 108.
17 min
60 min
120 min
180 min
R = 0.05
R = 0.12
240 min R = 0.19
R = 0.26
19
R = 0.33
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Note that the melt fraction is also targeted here. Fig. 6 depicts temporal variation of such a quantity at R = 0.12 and for different values of the Rayleigh number. It can be seen that the melt fraction exhibited similar profiles in the early times when the heat transfer is primarily dominated by conduction. Over time, the effect of Rayleigh number becomes more obvious from Ra = 1.34 × 108 since convection is the dominating heat transfer mechanism whose intensity increases the melting rate and the melt fraction of 18% for Ra = 1.34 × 108 and of 23% for Ra = 1.98 × 108. Based on our salient results, one can state that, at R = 0.12 and for Rayleigh numbers up to 7.8 × 107, the maximum value of melt fraction is of 48%. However, for a large Rayleigh number and at this ratio, the maximum melting fraction is 70%. We can already conclude that the latent heat thermal energy storage (LHTES) system, presented in Fig. 1, can only be effective if the external temperature is much higher than the melting temperature of the PCM. Further results are illustrated in Fig. 7. This shows the temperature profiles at the mid-height of the cavity. A careful review of this figure shows that the profile of temperature presents a change when the heat transfer mechanism passes from conduction to convection. For low Rayleigh numbers up to 7.8 × 107, the profile corresponds to a conduction regime. However, for a larger Rayleigh number, the curve becomes almost straight at the center, pointing a stratified flow. Fig. 8 shows the temporal evolution of the average Nusselt number parameterized by the Rayleigh number along the hot wall at R = 0.12. Initially, for all Rayleigh numbers, Nusselt numbers are very high at the initiation of the heat transfer process and decreases sharply during the conduction regime. It should be noted that high initial values of Nusselt number are due to the low thermal resistance of the very thin liquid layers at the fusion start up, which raises the heat transfer. At the end of the conduction regime, Nusselt numbers reach local minimum values. Furthermore, it can be seen that, for a lower Rayleigh number, the flow rate is more stable. At high Ra (= 1.34 × 108), a small oscillation occurs showing the formation of convection currents in the liquid PCM. This is due to the crucial role of convection on heat transfer. Then, after a transition period between conduction and convection, the Nusselt numbers attain quasi-stable values indicating that heat transfer is dominated by convection during the fusion process.
Fig. 10. Temporal evolution of the total liquid fraction.
Fig. 11. Temporal evolution of the average Nusselt number for various aspect ratios at Ra = 1.34 × 108.
5.2. Influence of geometry configuration to temperatures of 40, 50 and 60 °C of the hot wall, while keeping R = 0.12. The steady-state values of isotherms and streamlines at several times ranging from 17 to 240 min are shown in Fig. 5. By reviewing isotherms carefully, it appears that the intensity of the dominant heat transfer mechanism changes as Ra increases. For low Rayleigh numbers up to 7.8 × 107, the solid-liquid interface is almost flat since the convection is still weak and the heat is transferred primarily by conduction between hot and cold walls. With a large Rayleigh number (Ra = 1.34 × 108) and after 17 min, the convection process induced by buoyancy prevails. It can be noted that convective rollers move clockwise toward the melting interface. When these convective rolls reach the interface, thermal energy is transferred toward the solid and the liquid cools, which starts descending downwardly along the interface to the bottom of the cavity. Over time, the convection region extends downwardly. This leads to an irregular morphology of the melting front. Furthermore, increasing the Rayleigh number enhances the progress speed of melting interface and heat transfer rate of PCM 27 with a reduction in the melting time. In addition, one can notice that streamlines confirm, firstly, the presence of a one recirculation in the melting zone, and secondly that the heat transfer is stronger with a large Rayleigh number. This single structure is initially located at the top of the cavity and, over time (i.e. with the increase of the Fourier number), it goes to the bottom of the cavity. It is also observed that the melting rate occurs toward the bottom of the cavity, and that the liquid cools down as it descends along the interface. This expected result is mainly due to the presence of strong natural convection currents in the upper region of the enclosure.
Here, we seek to study the influence of the geometry of the cavity on the melting process and the heat transfer rate of PCM considered with reduction in the melting time. Thereby, we can determine the best configuration that can ensure good effectiveness of the melting process. This influence is investigated initially by varying the container's thickness, and then its height. 5.2.1. Influence of aspect ratio (with x-variation) Fig. 9 depicts the evolution of isotherms and melting interface at various instants and several aspect ratios R (0.05, 0.12, 0.19,0.26, 0.33), and for Ra = 1.34 × 108, the height being fixed at 21 cm. Through this figure, some important conclusions can be drawn: for low aspect ratios, heat transfer is primarily by conduction because of the small thickness of the cavity. One can state that the more the aspect ratio increases, the more heat transfer is dominated by convection and the greater is the rate of melt. Note that, initially, the heat transfer is governed by mixed conduction and convection. As time progresses, the initial conduction-convection dominated heat transfer process is gradually replaced by convection. The more the aspect ratio number increases, the more the whole enclosure becomes thermally stratified (isotherms are almost horizontal and parallel) at the steady state. Accordingly, such thermal stratification by convection decreases sharply and, as a result, the heat transfer rate also decreases. Furthermore, we find that the melting is faster at the top and slower at the bottom. In addition, the advancement rate of the solid-liquid interface decreases with the increase in aspect ratio, and the best configuration appears to 20
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Fig. 12. Isotherms and the melting interface at different times and several aspect ratios for Ra = 1.34 × 108.
17 min
60 min
120 min
180 min
R = 0.08
R = 0.1
240 min R = 0.12
R = 0.16
R = 0.23
the rectangular container thickness. In addition, since the conditioned material is intended to be integrated into the latent heat thermal energy storage system, and as it will undergo several melting-freezing cycles in the day, a compromise should be found between the melting rate and
match to R = 0.19. Temporal variations of the liquid fraction parameterized by the aspect ratio (R) are depicted in Fig. 10. It appears clearly that the liquid fraction and the melting rate are improved by a moderate increase in 21
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Fig. 13 presents the predicted temporal evolution of the liquid fraction with the aspect ratio as parameter. It can be seen that the liquid fraction profile remains the same for the conduction regime. Then, when the convection develops, the liquid fraction becomes a decreasing function of the aspect ratio. In addition, when the aspect ratio increases, the fusion velocity increases, and therefore the fraction of the melted material increases. Fig. 14 shows the temporal evolution of the heat transfer through the hot wall of the rectangular container. At the beginning of the melting process, the average Nusselt number at the left wall is very large for all aspect ratios due to the conduction mode. The more the process continues, the more the average Nusselt number decreases and becomes rather stable. Moreover, it is observed that changes in the average Nusselt number are similar, and this is for any aspect ratio. However, the heat transfer decreases when the aspect ratio number is high. In addition, for all values of the aspect ratio, the stable regime is reached very quickly, thus demonstrating that the melting is a quick process.
Fig. 13. Temporal evolution of liquid fraction with different aspect ratios.
5.3. Influence of Prandtl number To investigate the effect of Prandtl number on the melting process, we have performed simulations for different Pr ranging from Pr = 2.747 to Pr = 274.7, the Rayleigh number and the aspect ratio being fixed at 1.34 × 108, 0.12, respectively. It is worth recalling that the Prandtl value depends on the dynamic viscosity. The PCM's supplier provides a value of 10− 1 kg/(m ∙ s). So, it is this value that has been used in our previous publications [12,15]. Nevertheless, we aimed the effect of using dirty salts with a low viscosity value during heat transfer. Fig. 15 depicts the evolution of temperature contours and the melting interface according to the Prandtl number. A careful review of this figure shows that the effect of Pr on melting process is obvious for low values of Pr. At an early stage of the melting (17 min), regardless of the value of Prandtl, the heat is transferred mainly by conduction in the bottom of the cavity and by convection at the top of the cavity (especially for Pr = 274.7). As time progresses, the convection effect on the melting process is amplified and the amount of molten material is increased. It is clearly seen that the melting rate is more pronounced and that the melting interface is well advanced for Pr = 2.747 when compared to others (note that, no difference in enhancement has been observed at large Prandtl numbers). Based on these findings, we may state that the use of latent thermal energy storage material having a medium dynamic viscosity from 10− 3 may improve progression of the melting interface and temperature field. Fig. 16 depicts temporal variations of the liquid fraction at different Prandtl values. It is found that the fusion rate intensity increases when the Prandtl number decreases. This confirms our previous analysis on Fig. 15. It appears that the convection is better developed at low Prandtl numbers. Fig. 17 illustrates temporal evolutions of the average Nusselt numbers along the heated wall for various Prandtl numbers. As can be seen, shapes of the local Nusselt curve for different Prandtl numbers are overall similar. In addition, we highlight the two regimes of the melting process, namely conduction and convection. However, the presence of a trough is due to the high heat transfer, which causes a rise of the temperature gradient.
Fig. 14. Temporal evolution of the average Nusselt number for various aspect ratios and at Ra = 1.34 × 108.
the melting time. Thereby, the rectangular container with aspect ratio R = 0.19 offers best choice of the melting process. The time developments of heat transfer calculated along the heated wall for different aspect ratios when Ra = 1.34 × 108 are shown in Fig. 11. Initially, the average Nusselt number is very high for all aspect ratios because the heat transfer is governed by conduction. Then, the heat transfer decreases sharply followed by an oscillation in the transitional stage (with the exception for R = 0.05, where the heat transfer is mainly carried out via conduction). The intensity of this oscillation increases with increasing aspect ratio. This is due to the intensity of the convection. Then, the average Nusselt number decreases gradually until the state becomes completely stable. Furthermore, we notice that the heat transfer strongly depends on the aspect ratio, and that the average Nusselt number reduces as the aspect ratio increases. 5.2.2. Influence of aspect ratio (with y-variation) The temperature contours and melting interface at different times are depicted in Fig. 12 for five values of aspect ratio (the thickness being fixed at 2.5 cm), viz. R = 0.08 , 0.1 , 0.12 , 0.16 , 0.23. The aim of these simulations is to investigate by how much the heat transfer rate can be enhanced by increasing the height of the PCM container. As expected, a convective flow develops over time while being the main heat transfer mechanism. At the beginning of the melting, the heat transfer rate increases in the upper part of the cavity with the aspect ratio, and certainly, it has a great influence on the liquid-solid interface morphology. As time progressed, convection propagates down the cavity, especially for high values of aspect ratio. Likewise, an increasing aspect ratio number improves the melting heat transfer rate and increases the melting velocity.
6. Conclusions This study was motivated by the desire to numerically investigate the melting process of a macro-encapsulated hydrated salt phase change material in a rectangular container differentially-heated by the sides to simulate the solar irradiation received by the vertical wall. Intended results have been conducted via a finite volume/enthalpyporosity method termed hybrid model. An important implementation issue is to accurately solve the moving interface. Such an approach was 22
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Fig. 15. Isotherms and the melting interface at different times and various Prandtl numbers, a) Pr = 2.747; b) Pr = 27.47; c) Pr = 274.7; Ra = 1.34 × 108 and R = 0.12.
(a)
(b)
17 min
60 min
(c) 120 min
180 min
Fig. 16. Temporal evolution of the liquid fraction at different Prandtl numbers.
240 min
Fig. 17. Temporal evolution of the average Nusselt number for various Pr at Ra = 1.34 × 108 and R = 0.12.
validated first by comparing our findings with numerical and experimental data from the literature. Then, a parametric study was performed using the model developed to study effects of key parameters such as the Rayleigh number, the aspect ratio, and the Prandtl number on the effectiveness of the melting process. Through this study, the following conclusions and recommendations can be drawn:
•
• The heat transfer characteristics of PCM 27 are assessed by in-
specting the isotherms, the streamlines, location of melting 23
interface, the average Nusselt number, and the liquid fraction. From the obtained findings, it appears that the fusion process is promoted and the efficiency of heat transfer is improved by increasing the thickness of the material, and by decreasing its viscosity (low values of Pr). At the early stage of melting, the heat transfer is mainly governed by conduction indicated by a solid-liquid interface parallel to the vertical heated wall. At later time instants, convection heat transfer is
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
• • • •
(LHTESS), Renew. Sust. Energ. Rev. 14 (2010) 615–628. [15] A. Joulin, Z. Younsi, L. Zalewski, S. Lassue, D.R. Rousse, J.P. Cavrot, Experimental and numerical investigation of a phase change material: thermal-energy storage and release, Appl. Energy 88 (2011) 2454–2462. [16] M. Costa, D. Buddhi, A. Oliva, Numerical simulation of a latent heat thermal energy storage system with enhanced heat conduction, Energy Convers. Manag. 39 (3–4) (1998) 319–330. [17] E. Halawa, F. Bruno, W. Saman, Numerical analysis of a PCM thermal storage system with varying wall temperature, Energy Convers. Manag. 46 (2005) 2592–2604. [18] G. Hed, R. Bellander, Mathematical modelling of PCM air heat exchanger, Energ. Buildings 38 (2006) 82–89. [19] D. Buddhi, N.K. Bansal, R.L. Sawhney, M.S. Sodha, Solar thermal storage systems using phase change materials, Int. J. Energy Res. 12 (3) (1988) 547–555. [20] Z. Younsi, A. Joulin, L. Zalewski, S. Lassue, D.R. Rousse, Phase change materials: a numerical method for the behavior predictions, Proceedings of the Fourth International Conference on Thermal Engineering: Theory and Applications, Abu Dhabi, UAE, January 12–14, 2009. [21] F. Rosler, D. Bruggemann, Shell-and-tube type latent heat thermal energy storage: numerical analysis and comparison with experiments, Heat Mass Transf. 47 (8) (2011) 1027–1033. [22] M. Jourabian, M. Farhadi, A.A.R. Darzi, Outward melting of ice enhanced by Cu nanoparticles inside cylindrical horizontal annulus: lattice Boltzmann approach, Appl. Math. Model. 37 (20 − 21) (2013) 8813–8825. [23] Z.X. Gong, S. Devahastin, A.S. Mujumdar, Enhanced heat transfer in free convection-dominated melting in a rectangular cavity with an isothermal vertical wall, Appl. Therm. Eng. 19 (12) (1999) 1237–1251. [24] C. Gau, R. Viskanta, Melting and solidification of a pure metal on a vertical wall, J. Heat Transf. 108 (1986) 174–181. [25] 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, Numer. Heat Transf., Part A 13 (1988) 297–318. [26] N. Shamsundar, E.M. Sparrow, Analysis of multidimensional conduction phase change via the enthalpy model, J. Heat Transf. ASME 97 (1975) 333–340. [27] M.A. Hamdan, F.A. Elwerr, Thermal energy storage using a phase change material, Sol. Energy 56 (2) (1996) 183–189. [28] N.S. Dhaidan, J.M. Khodadadi, T.A. Al-Hattab, S.M. Al-Masha, Experimental and numerical investigation of melting of phase change material/nanoparticle suspensions in a square container subjected to a constant heat flux, Int. J. Heat Mass Transf. 66 (2013) 672–683. [29] Z. Younsi, Etude expérimentale et numérique du comportement thermique de matériaux à changement de phase. Intégration dans un composant solaire passif pour l'habitat, Ph.D. Thesis Université d'Artois, 2008. [30] Z.X. Gong, A.S. Mujumdar, Flow and heat transfer in convection-dominated melting in a rectangular cavity heated from below, Int. J. Heat Mass Transf. 41 (17) (1998) 2573–2580. [31] R. Viswanath, Y. Jaluria, A comparison of different solution methodologies for melting and solidification problems in enclosures, Numer. Heat Transfer, Part B 24 (1) (1993) 77–105. [32] M. Lacroix, Computation of heat transfer during melting of a pure substance from an isothermal wall, Numer. Heat Transf., Part B 15 (1989) 191–210. [33] C.P. Dessai, K. Vafai, A unified examination of the melting process within a twodimensional rectangular cavity, ASME Trans. J. Heat Transfer. 115 (4) (1993) 1072–1075. [34] Z. Younsi, S. Lassue, L. Zalewski, D.R. Rousse, A. Joulin, A novel technique for the experimental thermophysical characterization of phase change materials, Int. J. Thermophys. 32 (2010) 674–692. [35] R. Viskanta, Phase change heat transfer, Solar Heat Storage: Latent Heat Materials, Edited by G.A. Lane, CRC Press, Boca Raton, Florida, 1983, pp. 153–222. [36] R. Viskanta, Natural convection melting and solidification, Natural Convection: Fundamentals and Applications, Hemisphere Publishing Company, Washington DC, 1985. [37] Y. Zhang, Z. Chen, Q. Wang, Q. Wu, Melting in an enclosure with discrete heating at constant rate, Exp. Thermal Fluid Sci. 6 (1993) 196–201. [38] P. Brousseau, M. Lacroix, Study of the thermal performance of a multi-layer PCM storage unit, J. Energy Convers. Manag. 37 (5) (1996) 599–609. [39] R. Viskanta, Heat transfer during melting and solidification of metals, J. Heat Transfer. Trans. ASME 110 (1988) 1205–1219. [40] P. Prescott, F.P. Incropera, Convection heat and mass transfer in alloy solidification, Adv. Heat Tran. 28 (1996) 231–338. [41] M. Lacroix, V.R. Voller, Finite difference solutions of solidification phase change problems: transformed vs. fixed grids, Numer. Heat Transfer 17 (1990) 25–41. [42] S.P. Ketkar, M. Parang, R.V. Arimilli, Experimental and numerical study of solidification and melting of pure materials, J. Thermophys. 5 (1991) 40–45. [43] R. Viswanath, Y. Jaluria, Numerical study of conjugate transient solidification in an enclosed region, Numer. Heat Transfer, Part A 27 (1995) 519–536. [44] M. Costa, A. Oliva, C.D.P. Segarra, R. Alba, Numerical simulation of solid-liquid phase change phenomena, Comput. Methods Appl. Mech. Eng. 91 (1991) 1123–1134. [45] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, 1980. [46] S. Kim, M.C. Kim, S.B. Lee, Prediction of melting process driven by conductionconvection in a cavity heated from the side, Korean J. Chem. Eng. 18 (5) (2001) 593–598.
promoted and this leads to an irregular morphology of the melting front. The time required for material melting is proportional to the Rayleigh number and the heat transfer by convection. For a fixed number of Rayleigh, the more the aspect ratio increases, the more the heat transfer is intense. In addition, heat transfer is higher with a large Rayleigh number, the aspect ratio being fixed. The significant rise in the thickness of the container induces the establishment of a thermal stratification and the liquid remains within the top portion of the enclosure. Under dimensions adopted herein, the best choice of the melting process is achieved in a cavity whose aspect ratio is 0.19 (H = 21 cm; W = 4 cm).
To sum up, our findings show that the efficiency of the macro-encapsulated PCM in a rectangular container can be improved. The incorporation of such a composite material within the latent heat thermal energy storage system can improve comfort, and the obtained results should be useful in helping to reduce the energy consumption in buildings. Conflict of interest statement The authors declare no conflict of interests regarding authorship and/or publication of this paper. Acknowledgments The authors are grateful to anonymous referees for their insightful review and rewarding comments that helped in improving the original manuscript. Likewise, we further thank the editor-in-chief for his continuing incentives. References [1] A. Sharma, V.V. Tyagi, C.R. Chen, D. Buddhi, Review on thermal energy storage with phase change materials and applications, Renew. Sust. Energ. Rev. 13 (2) (2009) 318–345. [2] B. Zalba, J.M. Marin, L.F. Cabeza, H. Mehling, Review on thermal energy storage with phase change: materials, heat transfer analysis and applications, Appl. Therm. Eng. 23 (3) (2003) 251–283. [3] M.M. Farid, A.M. Khudhair, S.A.K. Razack, S. Al-Hallaj, A review on phase change energy storage: materials and applications, Energy Convers. Manag. 45 (2004) 1597–1615. [4] R. Jacob, F. Bruno, Review on shell materials used in the encapsulation of phase change materials for high temperature thermal energy storage, Renew. Sust. Energ. Rev. 48 (2015) 79–87. [5] Y. Tian, C.Y. Zhao, A numerical investigation of heat transfer in phase change materials (PCMs) embedded in porous metals, Energy 36 (2011) 5539–5546. [6] B. Kamkari, H. Shokouhmand, F. Bruno, Experimental investigation of the effect of inclination angle on convection-driven melting of phase change material in a rectangular enclosure, Int. J. Heat Mass Transf. 72 (2014) 186–200. [7] Y. Tian, C.Y. Zhao, Thermal and exergetic analysis of metal foam enhanced cascaded thermal energy storage (MF-CTES), Int. J. Heat Mass Transf. 58 (2013) 86–96. [8] M. Jradi, M. Gillott, S. Riffat, Simulation of the transient behaviour of encapsulated organic and inorganic phase change materials for low-temperature energy storage, Appl. Therm. Eng. 59 (2013) 211–222. [9] A. Koca, H. Oztop, T. Koyun, Y. Varol, Energy and exergy analysis of a latent heat storage system with phase change material for a solar collector, Renew. Energy 33 (4) (2008) 567–574. [10] A. Khalifa, E. Abbas, A comparative performance study of some thermal storage materials used for solar space heating, Energ. Buildings 41 (2009) 407–415. [11] J.N.W. Chiu, V. Martin, Submerged finned heat exchanger latent heat storage design and its experimental verification, Appl. Energy 93 (2012) 507–516. [12] Z. Younsi, S. Lassue, L. Zalewski, D.R. Rousse, A. Joulin, A novel technique for the experimental thermophysical characterization of phase change materials, Int. J. Thermophys. 32 (2010) 674–692. [13] V.V. Tyagi, D. Buddhi, Thermal cycle testing of calcium chloride hexahydrate as a possible PCM for latent heat storage, Sol. Energy Mater. Sol. Cells 92 (2008) 891–899. [14] F. Agyenim, N. Hewitt, P. Eames, M. Smyth, A review of materials, heat transfer and phase change problem formulation for latent heat thermal energy storage systems
24
Contents lists available at ScienceDirect
International Communications in Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ichmt
A numerical investigation of melting phase change process via the enthalpyporosity approach: Application to hydrated salts
MARK
Zohir Younsia,b,c, Hassane Najib,c,⁎ a b c
FUPL, Hautes Etudes d'Ingénieur (HEI), LGCgE (EA 4515), 13 rue de Toul, F-59000 Lille, France Univ. Artois, Laboratoire Génie Civil & géo-Environement (LGCgE EA 4515), Technoparc Futura, F-62400 Béthune, France Université Lille Nord de France, LGCgE (EA 4515), Lille F-59000, France
A R T I C L E I N F O
A B S T R A C T
Keywords: Phase change material Natural convection Numerical simulation Finite volume method Enthalpy-porosity formulation Melting Fusion
The convection dominated melting of hydrated salts (PCM 27) is investigated numerically using a pressure based finite volume method (FVM) with an enthalpy porosity technique, which allows a fixed-grid solution of the coupled momentum and energy equations, and trace the solid-liquid interface without resorting to other equations or transformations. This investigation addresses effects of Rayleigh number, the aspect ratio and Prandtl number on the effectiveness of melting process and solid-liquid interface. In addition, the average Nusselt number is set to complete the parametric study conducted. Through the results achieved, the enthalpyporosity has been found to converge rapidly while producing accurate results for both the position and morphology of the melt front at different times. The results are presented in terms of temperature profiles, streamlines, isotherms, moving interface position, liquid fraction and average Nusselt number. From this survey, it appears that our numerical findings support the experimental observations. These results indicated that, in terms of efficiency of the melting process, the best configuration of the rectangular container is the cavity having an aspect ratio of 0.19 (H = 21 cm ; W = 4 cm). Therewith, it is found that the time required to melt the material is proportional to the Rayleigh number and the convection heat transfer. Finally, this study demonstrates the feasibility for simulating macro-encapsulated hydrated salts via the enthalpy-porosity-based approach.
1. Introduction Storage of thermal energy is essential for reducing the continuous increase of the gap between energy demand and supply. There are three different ways one can store thermal energy in a storage material: form of sensible heat, latent heat and thermochemical energy. In the thermochemical storage system the corresponding storage materials used are still too expansive and their use requires a technical and economic feasibility study [1]. The sensible heat thermal storage (SHTS) system is very simple and the temperature of the storage material changes with the quantity of energy added or removed from the system. On the other hand, in a latent-heat thermal energy storage (LHTES) system, the density of energy storage is higher allowing the design of more compact systems with less heat losses and allowing energy storage at a nearly constant temperature [1,2]. Several studies have shown that using phase change materials (PCMs) as promising LHTES materials allows absorbing or releasing latent heat during a phase change from solidstate to liquid-state (melting process) or liquid-state to solid-state (freezing process). This system found several potential engineering ⁎
applications, such as thermal energy storage, thermal management, heating and cooling systems, electronic products, drying technology, waste heat recovery, solar air collectors, and solar cookers, and so on. One disadvantage of LHTES systems is the lower thermal conductivities of PCMs used, which slowed the melting/solidification rate. A faster fusion/solidification rate is highly sought after in the field of engineering [3]. Among the techniques for improving the thermal conductivity of PCMs and the heat transfer are microencapsulated PCM [4], porous metals [5], use of stationary highly conductive structures [6] and multiple PCMs [7]. In addition to these approaches, optimization of the conditioning system of PCMs along with the choice of type of PCM is the main interest of this research. These parameters modify the phase change process, as well as the heat transfer rate, which will be discussed in the present work. Since organic and inorganic PCMs are widely used in various engineering fields, they have been extensively studied in literature. For low-temperature residential applications, the hydrated salts (inorganic PCMs) are the most commonly used because of their low cost, low volumic expansion, stability of melting point and high volumetric
Corresponding author at: Univ. Artois, Laboratoire Génie Civil & géo-Environement (LGCgE EA 4515), Technoparc Futura, F-62400 Béthune, France. E-mail address: [email protected] (H. Naji).
http://dx.doi.org/10.1016/j.icheatmasstransfer.2017.05.012
0735-1933/ © 2017 Elsevier Ltd. All rights reserved.
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
ΔV x, y
Nomenclature A B C c f g H h k L W H Nu Pr p R Ra S Ste T T0 Tm TC TH t u v
matrix coefficients porosity function morphological constant (see Relationship (5)) specific heat (J/kg ∙K) liquid fraction gravity vector (m/s2) total enthalpy (J/m3) specific enthalpy (J/m3) thermal conductivity (W/m ∙K) latent heat of fusion (J/kg) length of rectangular enclosure in x-direction (m) length of rectangular enclosure in y-direction (m) Average Nusselt number Prandtl number pressure (Pa) aspect ratio (H/W) Rayleigh number source term Stefan number temperature (K) initial temperature (K) melting point of materials (K) isothermal right wall temperature (K) isothermal left wall temperature (K) time (s) velocity in the x direction (m/s) velocity in the y direction (m/s)
control volume Cartesian coordinates (m)
Greek symbols α β ε ∇ ν Δt Γ μ ρ φ ω
thermal diffusivity (m2/s) expansion coefficient (K− 1) arbitrary computational constant Nabla operator kinematic viscosity (m2/s) time step (s) diffusion coefficient dynamic viscosity (kg/m ∙s) density (kg/m3) generalized variable relaxation factor
Superscripts/subscripts c interface position k iteration level l liquid phase m melting P present node ref reference temperature value S solid phase W, E, P, N, S west, east, center, north and south nodes 0 initial value
of the cavity. Brent et al. [25] were among first authors to develop the enthalpy-porosity technique to trace the solid-liquid interface during the melting and solidification process. This approach has been validated by replicating the results of experimental work [24] using a 2D rectangular cavity. Shamsundar and Sparrow [26] proposed a multi-dimensional conduction phase change model for a solar energy storage system using the enthalpy-based method. Hamdan and Elwerr [27] have developed a mathematical approach to handle solid-liquid interface propagation during the two-dimensional melting process. Dhaidan et al. [28] have explored experimentally and numerically the melting of n-octadecane with CuO nanoparticle suspensions in a square enclosure. Their findings showed a significant effect of natural convection and an increase of the Rayleigh number on the shape of the solid-liquid interface. One of the most successful applications for PCM is thermal energy storage for air heating. Younsi et al. [29] have provided a new air heating system with a PCM energy store called Trombe wall improved, and allowing reducing energy consumption by 30% compared to conventional wall systems. As shown in Fig.1, the Trombe wall improved
energy storage density [8]. In this context, the calcium chloride hexahydrate CaCl2 ∙ 6H2O is one of potential materials for the low-temperature thermal storage, and thereby it has been explored extensively in the literature research [9–12]. Younsi et al. [12] have investigated experimentally the melting characteristics of salt hydrates in a rectangular cavity heated through its sides. Their results indicated that the melting efficiency has great potential for storing heat. Tyagi and Buddhi [13] have experimentally studied the thermal stability of CaCl2 ∙6H2O after any thermal cycling tests. They concluded that CaCl2 ∙ 6H2O is the material which has a promising potential for heating and cooling applications with small variations in the material latent heat of fusion during thermal cycling [8]. Macro-encapsulation of PCMs is a promising technic for storing the latent heat energy where the material is sealed in an enclosure which can be of any shape (rectangular, cylindrical, and spherical). It should be noted that the choice of the enclosure's geometry is a key parameter for improving the heat transfer rate. The rectangular and cylindrical containers are the most used due to their facility and low production cost [14]. To assess the potential use of macro-encapsulated PCMs in rectangular containers, many numerical, analytical and experimental studies have been conducted over the last decade [15–19]. The melting of PCMs aroused great interest because of its important role in the thermal energy storage. Several analytical, numerical and experimental studies were carried on phenomena associated with melting in different geometry configurations [20–23]. These works have concluded that the onset of melting is clearly dominated by diffusion, while buoyancy-driven convection has a paramount role in the melting process. However, little effort has been made on the optimization of the melting of a PCM packaged in a rectangular envelope and lay out vertically. Gau and Viskanta [24] have investigated experimentally solid-liquid interface motion and heat transfer during melting and solidification in a rectangular cavity filled with Gallium (Ga). They found that the heat transfer at the beginning of melting is dominated by conduction as time progresses, and that the convection develops on top
Fig. 1. Illustration of a latent heat thermal energy storage system [29].
13
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
system utilizes a thermal energy storage unit and an intelligently controlled air distribution unit. The principle of such a system can be outlined briefly as follows: the air between the glazing and the PCM is heated by greenhouse effect and allows charging the PCM. The PCM is discharged by passing air through the space between the wall and the insulation allowing air to be heated and providing thermal comfort for occupied spaces. Based on the review conducted, several approaches have been developed to handle the problem of the phase change, and many researchers have developed numerical models to assess the effectiveness of phase change materials packaged in rectangular containers. To go further in that sense, here, we aim to achieve the design of a latent heat thermal energy storage system by numerically addressing the problem of melting natural convection in a rectangular container heated from one of its sides. Such a problem is handled by using the enthalpy-porosity-based approach. To assess the accuracy of this formulation, main findings are compared with available results of the literature results [24,30–33] on ratios' aspect. In addition, simulations have been carried out to assess effects of relevant parameters (e.g. Rayleigh number, aspect ratio and Prandtl number) on the melting process. Likewise, the fusion time of the hydrated salts in a rectangular container is compared to these parameters. A complete physical overview of the dynamic and thermal fields in terms of streamlines, isotherms and the average Nusselt number along the hot wall are presented. In the light of the results achieved, it should be stated that the required time for good PCM melting inside the wall is highly controlled both by the container shape and heat flux imposed. The remainder of this paper is organized as follows. The problem description and mathematical formulation are introduced in Section 2. Subsequently, the computational methodology is briefly outlined in Section 3. In this section, the enthalpy FV-based approach with the discretization of the governing equation along with the Thomas algorithm is briefly reported. Afterwards, in Section 4, a validation of the adopted numerical approach via available results has been performed. Section 5 presents and comments our findings for a range of relevant parameters, such as Rayleigh number, aspect ratio and Prandtl number. Finally, the main conclusions are addressed in Section 6.
Table 1 Thermophysical properties of the PCM 27. Material
Solid phase Liquid phase
Thermophysical properties Density (kg/m3)
Specific heat (J/kg ∙K)
Thermal conductivity (W/ m ∙K)
Latent heat of fusion (kJ/ kg)
1710 1530
1751.5 2225
0.577 0.813
– 172.42
hydrated salt (PCM 27) having a melting temperature of 27 °C. Such a model has been studied experimentally by [34]. The horizontal walls are perfectly insulated (adiabatic) while the two left and right vertical walls are maintained isothermal at constant temperatures TC (cold) and TH (hot), respectively. The PCM in its fluid phase is assumed to be viscous, homogeneous; isotropic, incompressible and Newtonian. Its flow is two-dimensional, laminar and unsteady, without viscous dissipation or radiative heat transfer. In addition, the liquid motion induced by volumetric variation during melting is assumed to be neglected, thereby pointing that the solid density and that of the liquid are here equal. For the purpose of the numerical simulation, a no-slip condition holds at all the walls. Moreover, all thermo-physical properties of the considered PCM are assumed constant except the density variation in the buoyancy term which is modeled by the Boussinesq approximation for involving thermal buoyancy. These thermo-physical properties of these components are gathered in Table 1. Note that, in this analysis, it is assumed that the initial temperature of the PCM is equal to TC. Thereby, thermal energy will start penetrating through the left wall due to the temperature difference and will initiate the melting process. 2.2. Mathematical modeling Based on the aforementioned assumptions, the governing equations for the 2D transient laminar flow including heat transfer melting process and buoyancy-driven convection may thus be written as follows: Solid region:
2. Problem description and governing equations
∂T = αs ∇2 T ∂t
2.1. Geometry
(1)
Liquid region: The model enclosure is a 2-D square cavity with dimensions W × H on the x–y plane as depicted in Fig. 2. This enclosure is filled with a
→ ∇⋅ V = 0
(2) Fig. 2. Configuration under study (left) along with coordinates and boundary conditions (right).
14
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
→ → → → → ∂V + ( V ⋅∇) V = −∇ (p ρl ) + νl ∇2 V + β (T − Tm )→ g + S→ V ∂t
(3)
∂T + (V ⋅∇) T = αl ∇2 T + Sh ∂t
(4)
In Eqs. (3)–(4), source term S is, respectively given by:
→→ → L ∂f S V = −B (f ) V and Sh = − cl ∂t
h=
∫T
T
ρcdT
→ → → → → → g ∂V + ( V ⋅∇) V = −∇P + Pr∇2 V + Ra Pr θ → + BV ∂t͠ |g |
(11)
→ 1 ∂f ∂θ + ( V ⋅∇) θ = ∇2 θ − Ste ∂t͠ ∂t͠
(12)
→ 1 ∂f ∂θ + ( V ⋅∇) θ = α͠ ∇2 θ − Ste ∂t͠ ∂t͠
• left wall: • right wall: • bottom wall: • top wall: • solid walls:
(7)
Pr = ν α, Ra = gβH 3 (TH − Tm) Ste = cl (TH − Tm) Lf and R = W H αν , (14) The average number of Nusselt, which is a dimensionless measure of heat transfer through the hot wall can be computed as follows:
Nu = u = v = 0 and T = TH u = v = 0 and T = TC u = v = 0 and ∂ T/∂y = 0
1 R (Tm − TH )
H
∫ (∂T
∂x )x = 0 dy (15)
0
Prior to implementing the numerical algorithm, Eqs. (10)–(13) have been recast in the following general equation of a conservation equation for a general variable φ:
u = v = 0 and ∂ T/∂y = 0 u = v = 0 and T = Tm
→ ∂ (φ) + ∇⋅( V φ) = ∇⋅(Γφ ∇φ) + Sφ ∂t
(16)
In this equation, the expression of diffusion coefficients Γφ and the source term Sφ of each transported variable φ are gathered in Table 2.
→ → T(x, y, 0) = Tini = Tc and V = 0 at solid walls.
where TH > Tm > TC and s is the horizontal distance between the left isothermal wall and the melting front. Note that TC is kept constant in all simulations performed here. To cast the above equations in a dimensional form, we employ the following dimensionless variables:
3. Computational procedure Since the numerical method is adopted from the existing literature, we give here only a brief outline of the major features of the method used. The interested reader can refer for example to Ref. [41] for complete details on such a method. Integrating Eq. (16) over a control
t͠ = tαl H 2, (X , Y ) = (x H , y H ), (U , V ) = (uH αl , vH αl ), θ = (T − Tm) (TH − Tm), P = pH 2 ραl2 , α = αs αl, B∗ = BH 2 ραl
Table 2 Diffusion coefficient Γφ and source term Sφ for variable φ in the generic transport equation.
(8) Based on these dimensionless parameters, the aforementioned Eqs. (1)–(5) are then written as follows: Solid region
∂θ = α∇2 θ ∂t͠
θ|X = 0 = 1 θ|X = W/H = (TC − Tm)/(TH − Tm) ∂ θ/∂Y |Y = 0 = 0 ∂ θ/∂Y |Y = 1 = 0 → → V (U , V ) = 0
The dimensionless initial condition is given by: → → V = 0 , θ = 0 at t͠ = 0 , 0 < (X, Y) ≤ 1. It should be noted that such a problem involves the Prandtl number, the Rayleigh number, the Stefan number and the aspect ratio, defined by:
where ρfL is latent heat, which varies between zero for solid to ρL for liquid. Note that using this formulation, the absorption of latent heat during melting is isolated in the source term Sh. In the past, this method has been validated many times for convection-dominated melting problems in a rectangular container with different PCMs (Gallium, tin, n-Octadecane, caprillic acid) [25,42–44]. The required boundary and initial conditions for the considered problem are
• left wall (x = 0): • right wall (x = W): • bottom wall (y = 0): • top wall (y = H): • melting front (x = s): • initial conditions:
(13)
where α͠ = f + (1 − f ) α . In this way, in the liquid zone (θ > 0, f = 1), α͠ = 1, and in the solid zone (θ < 0, f = 0), α͠ = α . Note that the term (S͠ h = −Ste−1∂f ∂t͠ ) is the latent heat term (source) due to phase change. These equations are subjected to the following boundary conditions:
(6)
m
(10)
Note that Eqs. (9) and (12) can be written as follows:
(5)
where B(f) = C(1 − f)2/(f3 + ε) causes the gradual change in velocity from a finite value in the liquid to zero in the solid in the computational cells, f being the volume liquid fraction. Here, the multiplier C is the mushy zone constant (=109) reflecting the morphology of the melting/ solidification and ε is an arbitrary computational constant, set (= 3 × 10− 4 here) for avoiding a possible zero denominator (when → f = 0 corresponding to solid medium case); S → V is the momentum dissipation source term used to suppress the velocity in the mushy and solid zone, and Sh is the latent heat source term (caused by absorption or release). Gravity acts down (in the negative y-direction). As implied, an enthalpy-based method is adopted here for modeling the phase change problem [35–38]. Thereby, the energy equation follows Voller's formulation [39–41] for which the total enthalpy is split into sensible and latent enthalpy parts as
H (T ) = h (T ) + ρfL
→ ∇⋅ V = 0
(9)
Liquid region 15
φ
Γφ
Sφ
1 U V θ
0 Pr Pr f + (1 − f)αs/αl
0 −∇XP + BU Ra Pr θ − ∇YP + BV
−
∂f ∂t͠
Ste
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Ap φp =
Table 3 Thermophysical properties of Gallium.
∑
Ai φi + Sφ
i = W , E , S, N
cp (J/ kg ⋅ K)
k (W/ m ⋅ K)
Lf (kJ/kg)
ρ (kg/m )
υ (m /s)
β (K
302.93
381.5
32
80.16
6093
1.38 × 10− 5
1.2 × 10− 4
2
(17)
The iterative solution continues until convergence of the flow and energy fields at every time step is reached. Convergence is declared only when the three following conditions are met: 1) the liquid fraction field is stabilized, 2) the dimensionless residual for the enthalpy equation is less than 10− 5, and 3) the dimensionless residual for the mass conservation equation is less than 10− 4. These tolerances are chosen so that the final solution did not show perceptible changes. The pertinent feature of the present enthalpy approach for modeling the phase change problem is the source term Sh for the energy equation (Eq. (13)). It is useful to recall that this term keeps track of latent heat evolution, and its driving element is the local liquid fraction. This term is an indicator which is 1 in fully liquid regions, 0 in completely solid regions, while remaining within the interval [0, 1] close to the phase front. In addition, a time-varying distribution of the local liquid fraction is achieved iteratively from the solution of the energy equation via the following expressions:
−1
Tm (K)
3
)
f Pk + 1 = f Pk + ωΔtAPk hPk ρLΔV
(18)
with Fig. 3. Progress of the melting front at different times. Comparison with Refs. [24,30,32] at Pr = 0.0216, Ste = 0.0391, R = 2, and Ra = 2.1 × 105.
f Pk + 1 =
k+1 ⎧1 if f P > 1 ⎨ 0 if f k + 1 < 0 P ⎩
(19)
where ω is an under-relaxation factor. It is worth noting that, after the kth-iteration, the liquid fraction update is applied at every node to reflect the evolution of the phase change. 4. Validation of the enthalpy-based approach Prior to computations, checks have been conducted to validate our simulations compared to other results which have been reported in literature [24,30,31–33]. For this purpose, two benchmark cases have been carried out to demonstrate the approach herein considered. The first test case is one that had been handled by Gau et al. [24], Gong and Mujumdar [30], and Lacroix [32]. This is two-dimensional problem of melting of Gallium with the presence of natural convection in a rectangular enclosure at Pr = 0.0216, Ste = 0.0391, and Ra = 2.1 × 105. The dimensions of the cavity are H = 0.0445 m and W = 0.089 m (aspect ratio R = 2). The left hot wall and the right cold wall are maintained at temperatures, TH = 303.35 K and TC = 302.5 K, respectively. The horizontal walls are adiabatic. The relevant thermophysical properties of Gallium are listed in Table 3. Several numerical investigations were conducted to check the grid size and time step dependence results. A uniform grid size (50 × 50), time step (Δ t = 1 s), and under-relaxation factors for the momentum and energy equation (ωh = 0.6 and ωu = ωv = 0.8) were found sufficient to give accurate results. Several works were performed on the fusion of gallium and shows that most authors specify that the results are independent of the mesh. Viswanath [31] showed that the mesh refinement of 30 × 20–50 × 40 leads to a change of 3% of the final solution and concluded that this change is acceptable. The melting front is plotted at several times during the melting process, i.e., 5, 8 and 12 min (see Fig. 3), and compared favorably with the results of Gau and Viskanta [24], Gong and Mujumdar [30] and Lacroix [32]. In other words, it appears that experimental melt front positions are quite well reproduced with a slight discrepancy. Such a difference may be attributed to 1) the ignorance of the solid phase subcooling (2 °C under the melting temperature) by the current model, 2) the low wall temperature (8 °C above the melting) and the high metal conductivity (k = 32 W ∙ m− 1 ∙ K− 1), for which the metal sub-cooling seems to influence the experiment running (see [15])[, and 3]) the difficulty of maintaining vertical walls at a desired temperature during
Fig. 4. Progress of the melting front at different times. Comparison among results of Refs. [24,31,33] at Pr = 0.0216, Ste = 0.0391, R = 1.4 and Ra = 6.057 × 105. Table 4 Dimensionless parameters for the melting of PCM hydrated salts. R
Ra
Pr
Ste
cs/cl
ks/kl
0.12
1.34 × 108
274.7
0.2968
0.79
0.71
volume and estimation of each term of the integral equation in terms of discrete value on the nodal point lead to the production of a system of algebraic equations. To insure stability, a first order hybrid scheme [45] is, however, used for the discretization of the convection-diffusion terms in the momentum and energy equations. A fully implicit Euler method is used for the time-stepping procedure. Algebraic systems of equations originating from spatial and temporal discretizations are solved iteratively (throughout the physical domain considered) using the Tri-Diagonal Matrix Algorithm (by sweeping the cavity from left to right). On the other hand, the Semi-Implicit Algorithm for Pressure Linked Equations (SIMPLER) algorithm has been chosen to handle the pressure-velocity coupling. Note that the discretized equations for the dependent variable φ are written for each node P as:
16
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Fig. 5. Isotherms, streamlines and melting interface at different instants and various Rayleigh numbers with R = 0.12; a) Ra = 7.8 × 107; b) Ra = 1.34 × 108; and c) Ra = 1.98 × 108.
a) Ra =7.8×107 ; TH = 40 °C
b) Ra =1.34×108 ; TH = 50 °C
17 min
60 min
120 min 180 min 8 c) Ra = 1.98×10 ; TH = 60 °C
17
240 min
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
temperature-based method. At 8 min, the gap decreases between results, especially between our current results and the Gong and Mujumdar [30] experimental results. As time proceeds (at 12 min), our findings reproduce well those of Gau and Viskanta [24], except for the bottom region. It should be pointed out that, at the enclosure top end, all results tend to converge. The second benchmark test deals with the two-dimensional problem of the fusion of pure Gallium a rectangular cavity having an aspect ratio of 1.4. The hot and cold walls are isothermal (TH = 311.15 K and TC = 301.45 K), while horizontal walls remain adiabatic. In addition, the initial Gallium temperature was equal to TC. Fig. 4 gathers the comparisons concerning to the melting front progression at different times, viz. 5, 10 and 19 min. It can be seen that the comparison of the present results with those of the selected authors remains quite favorable. More precisely, at 5 min and 10 min, the agreement is excellent while at t = 19 min, some deviation is exhibited between all numerical simulations and experimental data. At this time, the present model simulations appear to be slightly offset as the melting front is behind that proposed by other solutions at the cavity top, while it is in front at the cavity bottom. To support this, it should be noted that Kim et al. [46] stated that the delayed heat-up affects the melting process and thereby can lead to a discrepancy between numerical and experimental results. Indeed, considering a delayed heat-up, they obtain a better agreement, especially at 19 min. To sum up, it seems that under these benchmark cases, the current model corroborates available results in literature. Thereby, it can be stated that our findings give confidence to the built code for the problem in this study.
Fig. 6. Temporal evolution of liquid fraction with different Ra (melting of PCM 27).
5. Results and discussion This section is dedicated to present and discuss the melting process in the enclosure considered in terms of influence of key parameters, such as the Rayleigh number, the aspect ratio and the Prandtl number on the dynamic and thermal fields, the melting fraction, and the average Nusselt number. In the following, some significant results showing effects of these parameters are comprehensively presented and discussed. To achieve this, a macro-encapsulated PCM 27 (hydrated salts) in a rectangular container is considered as a thermal energy storage system as shown in Fig. 2. The hexahydrate of calcium chloride (CaCl2 ∙6H2O) having a 27 °C melting temperature, is chosen as thermal latent energy storage material. Recall that the melting mode is investigated via the enthalpy-porosity approach. It should be pointed out that the main thermo-physical properties and dimensionless parameters of the sample are listed in Tables 1 and 4, respectively. In all cases, a 150 × 150 grid was used. Initially, the storage PCM is solid with a uniform temperature of 15 °C. Then, the temperature of the left vertical wall is suddenly raised to a prescribed temperature higher than the melting point, resulting in melting of hexahydrate of calcium. It should be noted that the efficiency of heat transfer with a PCM is thoroughly related to its rate melting. Therefore, the efficiency of the heat transfer and the rate at which the PCM melts greatly depends on the enclosure shape. Note that, various geometries containing the same PCM's amount were considered. These were either thick, thin, high or small. In this study, several efficiency indicators can be monitored, including mean and maximum cavity temperatures, total liquid fraction, rate at which PCM melts, and the Nusselt number. However, in the present study, we only retained the last 3 parameters, viz., the total liquid fraction, the PCM melting rate, and the Nusselt number.
Fig. 7. Temperature profiles at mid-height (y = H/2) of cavity for different Ra at 240 min.
Fig. 8. Temporal evolution of the average Nusselt number for various Ra numbers at R = 0.12.
the experiment. It should be noted that the sub-cooling is of approximately 2 °C in the solid phase. As for the discrepancy exhibited with Lacroix results, they can be imputed to the use of an enthalpy-porosity method with a fixed grid, instead of the temperature-based method involving a mobile grid to handle phase change effects. Despite this, it is observed that, at 5 min, the two approaches show similar results to the two experiments for the enclosure's bottom. Regarding the brick's top, the present model results seem to be in good agreement with those of Lacroix obtained with the
5.1. Influence of Raleigh number Recall that the Rayleigh number assesses the competition between convective and conduction heat transfer rates. To demonstrate the effect of this number, different values have been considered, viz., 7.8 × 107, 1.34 × 108 and 1.98 × 108, which correspond respectively 18
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Fig. 9. Isotherms and the melting interface at different times and several aspect ratios for Ra = 1.34 × 108.
17 min
60 min
120 min
180 min
R = 0.05
R = 0.12
240 min R = 0.19
R = 0.26
19
R = 0.33
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Note that the melt fraction is also targeted here. Fig. 6 depicts temporal variation of such a quantity at R = 0.12 and for different values of the Rayleigh number. It can be seen that the melt fraction exhibited similar profiles in the early times when the heat transfer is primarily dominated by conduction. Over time, the effect of Rayleigh number becomes more obvious from Ra = 1.34 × 108 since convection is the dominating heat transfer mechanism whose intensity increases the melting rate and the melt fraction of 18% for Ra = 1.34 × 108 and of 23% for Ra = 1.98 × 108. Based on our salient results, one can state that, at R = 0.12 and for Rayleigh numbers up to 7.8 × 107, the maximum value of melt fraction is of 48%. However, for a large Rayleigh number and at this ratio, the maximum melting fraction is 70%. We can already conclude that the latent heat thermal energy storage (LHTES) system, presented in Fig. 1, can only be effective if the external temperature is much higher than the melting temperature of the PCM. Further results are illustrated in Fig. 7. This shows the temperature profiles at the mid-height of the cavity. A careful review of this figure shows that the profile of temperature presents a change when the heat transfer mechanism passes from conduction to convection. For low Rayleigh numbers up to 7.8 × 107, the profile corresponds to a conduction regime. However, for a larger Rayleigh number, the curve becomes almost straight at the center, pointing a stratified flow. Fig. 8 shows the temporal evolution of the average Nusselt number parameterized by the Rayleigh number along the hot wall at R = 0.12. Initially, for all Rayleigh numbers, Nusselt numbers are very high at the initiation of the heat transfer process and decreases sharply during the conduction regime. It should be noted that high initial values of Nusselt number are due to the low thermal resistance of the very thin liquid layers at the fusion start up, which raises the heat transfer. At the end of the conduction regime, Nusselt numbers reach local minimum values. Furthermore, it can be seen that, for a lower Rayleigh number, the flow rate is more stable. At high Ra (= 1.34 × 108), a small oscillation occurs showing the formation of convection currents in the liquid PCM. This is due to the crucial role of convection on heat transfer. Then, after a transition period between conduction and convection, the Nusselt numbers attain quasi-stable values indicating that heat transfer is dominated by convection during the fusion process.
Fig. 10. Temporal evolution of the total liquid fraction.
Fig. 11. Temporal evolution of the average Nusselt number for various aspect ratios at Ra = 1.34 × 108.
5.2. Influence of geometry configuration to temperatures of 40, 50 and 60 °C of the hot wall, while keeping R = 0.12. The steady-state values of isotherms and streamlines at several times ranging from 17 to 240 min are shown in Fig. 5. By reviewing isotherms carefully, it appears that the intensity of the dominant heat transfer mechanism changes as Ra increases. For low Rayleigh numbers up to 7.8 × 107, the solid-liquid interface is almost flat since the convection is still weak and the heat is transferred primarily by conduction between hot and cold walls. With a large Rayleigh number (Ra = 1.34 × 108) and after 17 min, the convection process induced by buoyancy prevails. It can be noted that convective rollers move clockwise toward the melting interface. When these convective rolls reach the interface, thermal energy is transferred toward the solid and the liquid cools, which starts descending downwardly along the interface to the bottom of the cavity. Over time, the convection region extends downwardly. This leads to an irregular morphology of the melting front. Furthermore, increasing the Rayleigh number enhances the progress speed of melting interface and heat transfer rate of PCM 27 with a reduction in the melting time. In addition, one can notice that streamlines confirm, firstly, the presence of a one recirculation in the melting zone, and secondly that the heat transfer is stronger with a large Rayleigh number. This single structure is initially located at the top of the cavity and, over time (i.e. with the increase of the Fourier number), it goes to the bottom of the cavity. It is also observed that the melting rate occurs toward the bottom of the cavity, and that the liquid cools down as it descends along the interface. This expected result is mainly due to the presence of strong natural convection currents in the upper region of the enclosure.
Here, we seek to study the influence of the geometry of the cavity on the melting process and the heat transfer rate of PCM considered with reduction in the melting time. Thereby, we can determine the best configuration that can ensure good effectiveness of the melting process. This influence is investigated initially by varying the container's thickness, and then its height. 5.2.1. Influence of aspect ratio (with x-variation) Fig. 9 depicts the evolution of isotherms and melting interface at various instants and several aspect ratios R (0.05, 0.12, 0.19,0.26, 0.33), and for Ra = 1.34 × 108, the height being fixed at 21 cm. Through this figure, some important conclusions can be drawn: for low aspect ratios, heat transfer is primarily by conduction because of the small thickness of the cavity. One can state that the more the aspect ratio increases, the more heat transfer is dominated by convection and the greater is the rate of melt. Note that, initially, the heat transfer is governed by mixed conduction and convection. As time progresses, the initial conduction-convection dominated heat transfer process is gradually replaced by convection. The more the aspect ratio number increases, the more the whole enclosure becomes thermally stratified (isotherms are almost horizontal and parallel) at the steady state. Accordingly, such thermal stratification by convection decreases sharply and, as a result, the heat transfer rate also decreases. Furthermore, we find that the melting is faster at the top and slower at the bottom. In addition, the advancement rate of the solid-liquid interface decreases with the increase in aspect ratio, and the best configuration appears to 20
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Fig. 12. Isotherms and the melting interface at different times and several aspect ratios for Ra = 1.34 × 108.
17 min
60 min
120 min
180 min
R = 0.08
R = 0.1
240 min R = 0.12
R = 0.16
R = 0.23
the rectangular container thickness. In addition, since the conditioned material is intended to be integrated into the latent heat thermal energy storage system, and as it will undergo several melting-freezing cycles in the day, a compromise should be found between the melting rate and
match to R = 0.19. Temporal variations of the liquid fraction parameterized by the aspect ratio (R) are depicted in Fig. 10. It appears clearly that the liquid fraction and the melting rate are improved by a moderate increase in 21
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Fig. 13 presents the predicted temporal evolution of the liquid fraction with the aspect ratio as parameter. It can be seen that the liquid fraction profile remains the same for the conduction regime. Then, when the convection develops, the liquid fraction becomes a decreasing function of the aspect ratio. In addition, when the aspect ratio increases, the fusion velocity increases, and therefore the fraction of the melted material increases. Fig. 14 shows the temporal evolution of the heat transfer through the hot wall of the rectangular container. At the beginning of the melting process, the average Nusselt number at the left wall is very large for all aspect ratios due to the conduction mode. The more the process continues, the more the average Nusselt number decreases and becomes rather stable. Moreover, it is observed that changes in the average Nusselt number are similar, and this is for any aspect ratio. However, the heat transfer decreases when the aspect ratio number is high. In addition, for all values of the aspect ratio, the stable regime is reached very quickly, thus demonstrating that the melting is a quick process.
Fig. 13. Temporal evolution of liquid fraction with different aspect ratios.
5.3. Influence of Prandtl number To investigate the effect of Prandtl number on the melting process, we have performed simulations for different Pr ranging from Pr = 2.747 to Pr = 274.7, the Rayleigh number and the aspect ratio being fixed at 1.34 × 108, 0.12, respectively. It is worth recalling that the Prandtl value depends on the dynamic viscosity. The PCM's supplier provides a value of 10− 1 kg/(m ∙ s). So, it is this value that has been used in our previous publications [12,15]. Nevertheless, we aimed the effect of using dirty salts with a low viscosity value during heat transfer. Fig. 15 depicts the evolution of temperature contours and the melting interface according to the Prandtl number. A careful review of this figure shows that the effect of Pr on melting process is obvious for low values of Pr. At an early stage of the melting (17 min), regardless of the value of Prandtl, the heat is transferred mainly by conduction in the bottom of the cavity and by convection at the top of the cavity (especially for Pr = 274.7). As time progresses, the convection effect on the melting process is amplified and the amount of molten material is increased. It is clearly seen that the melting rate is more pronounced and that the melting interface is well advanced for Pr = 2.747 when compared to others (note that, no difference in enhancement has been observed at large Prandtl numbers). Based on these findings, we may state that the use of latent thermal energy storage material having a medium dynamic viscosity from 10− 3 may improve progression of the melting interface and temperature field. Fig. 16 depicts temporal variations of the liquid fraction at different Prandtl values. It is found that the fusion rate intensity increases when the Prandtl number decreases. This confirms our previous analysis on Fig. 15. It appears that the convection is better developed at low Prandtl numbers. Fig. 17 illustrates temporal evolutions of the average Nusselt numbers along the heated wall for various Prandtl numbers. As can be seen, shapes of the local Nusselt curve for different Prandtl numbers are overall similar. In addition, we highlight the two regimes of the melting process, namely conduction and convection. However, the presence of a trough is due to the high heat transfer, which causes a rise of the temperature gradient.
Fig. 14. Temporal evolution of the average Nusselt number for various aspect ratios and at Ra = 1.34 × 108.
the melting time. Thereby, the rectangular container with aspect ratio R = 0.19 offers best choice of the melting process. The time developments of heat transfer calculated along the heated wall for different aspect ratios when Ra = 1.34 × 108 are shown in Fig. 11. Initially, the average Nusselt number is very high for all aspect ratios because the heat transfer is governed by conduction. Then, the heat transfer decreases sharply followed by an oscillation in the transitional stage (with the exception for R = 0.05, where the heat transfer is mainly carried out via conduction). The intensity of this oscillation increases with increasing aspect ratio. This is due to the intensity of the convection. Then, the average Nusselt number decreases gradually until the state becomes completely stable. Furthermore, we notice that the heat transfer strongly depends on the aspect ratio, and that the average Nusselt number reduces as the aspect ratio increases. 5.2.2. Influence of aspect ratio (with y-variation) The temperature contours and melting interface at different times are depicted in Fig. 12 for five values of aspect ratio (the thickness being fixed at 2.5 cm), viz. R = 0.08 , 0.1 , 0.12 , 0.16 , 0.23. The aim of these simulations is to investigate by how much the heat transfer rate can be enhanced by increasing the height of the PCM container. As expected, a convective flow develops over time while being the main heat transfer mechanism. At the beginning of the melting, the heat transfer rate increases in the upper part of the cavity with the aspect ratio, and certainly, it has a great influence on the liquid-solid interface morphology. As time progressed, convection propagates down the cavity, especially for high values of aspect ratio. Likewise, an increasing aspect ratio number improves the melting heat transfer rate and increases the melting velocity.
6. Conclusions This study was motivated by the desire to numerically investigate the melting process of a macro-encapsulated hydrated salt phase change material in a rectangular container differentially-heated by the sides to simulate the solar irradiation received by the vertical wall. Intended results have been conducted via a finite volume/enthalpyporosity method termed hybrid model. An important implementation issue is to accurately solve the moving interface. Such an approach was 22
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
Fig. 15. Isotherms and the melting interface at different times and various Prandtl numbers, a) Pr = 2.747; b) Pr = 27.47; c) Pr = 274.7; Ra = 1.34 × 108 and R = 0.12.
(a)
(b)
17 min
60 min
(c) 120 min
180 min
Fig. 16. Temporal evolution of the liquid fraction at different Prandtl numbers.
240 min
Fig. 17. Temporal evolution of the average Nusselt number for various Pr at Ra = 1.34 × 108 and R = 0.12.
validated first by comparing our findings with numerical and experimental data from the literature. Then, a parametric study was performed using the model developed to study effects of key parameters such as the Rayleigh number, the aspect ratio, and the Prandtl number on the effectiveness of the melting process. Through this study, the following conclusions and recommendations can be drawn:
•
• The heat transfer characteristics of PCM 27 are assessed by in-
specting the isotherms, the streamlines, location of melting 23
interface, the average Nusselt number, and the liquid fraction. From the obtained findings, it appears that the fusion process is promoted and the efficiency of heat transfer is improved by increasing the thickness of the material, and by decreasing its viscosity (low values of Pr). At the early stage of melting, the heat transfer is mainly governed by conduction indicated by a solid-liquid interface parallel to the vertical heated wall. At later time instants, convection heat transfer is
International Communications in Heat and Mass Transfer 86 (2017) 12–24
Z. Younsi, H. Naji
• • • •
(LHTESS), Renew. Sust. Energ. Rev. 14 (2010) 615–628. [15] A. Joulin, Z. Younsi, L. Zalewski, S. Lassue, D.R. Rousse, J.P. Cavrot, Experimental and numerical investigation of a phase change material: thermal-energy storage and release, Appl. Energy 88 (2011) 2454–2462. [16] M. Costa, D. Buddhi, A. Oliva, Numerical simulation of a latent heat thermal energy storage system with enhanced heat conduction, Energy Convers. Manag. 39 (3–4) (1998) 319–330. [17] E. Halawa, F. Bruno, W. Saman, Numerical analysis of a PCM thermal storage system with varying wall temperature, Energy Convers. Manag. 46 (2005) 2592–2604. [18] G. Hed, R. Bellander, Mathematical modelling of PCM air heat exchanger, Energ. Buildings 38 (2006) 82–89. [19] D. Buddhi, N.K. Bansal, R.L. Sawhney, M.S. Sodha, Solar thermal storage systems using phase change materials, Int. J. Energy Res. 12 (3) (1988) 547–555. [20] Z. Younsi, A. Joulin, L. Zalewski, S. Lassue, D.R. Rousse, Phase change materials: a numerical method for the behavior predictions, Proceedings of the Fourth International Conference on Thermal Engineering: Theory and Applications, Abu Dhabi, UAE, January 12–14, 2009. [21] F. Rosler, D. Bruggemann, Shell-and-tube type latent heat thermal energy storage: numerical analysis and comparison with experiments, Heat Mass Transf. 47 (8) (2011) 1027–1033. [22] M. Jourabian, M. Farhadi, A.A.R. Darzi, Outward melting of ice enhanced by Cu nanoparticles inside cylindrical horizontal annulus: lattice Boltzmann approach, Appl. Math. Model. 37 (20 − 21) (2013) 8813–8825. [23] Z.X. Gong, S. Devahastin, A.S. Mujumdar, Enhanced heat transfer in free convection-dominated melting in a rectangular cavity with an isothermal vertical wall, Appl. Therm. Eng. 19 (12) (1999) 1237–1251. [24] C. Gau, R. Viskanta, Melting and solidification of a pure metal on a vertical wall, J. Heat Transf. 108 (1986) 174–181. [25] 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, Numer. Heat Transf., Part A 13 (1988) 297–318. [26] N. Shamsundar, E.M. Sparrow, Analysis of multidimensional conduction phase change via the enthalpy model, J. Heat Transf. ASME 97 (1975) 333–340. [27] M.A. Hamdan, F.A. Elwerr, Thermal energy storage using a phase change material, Sol. Energy 56 (2) (1996) 183–189. [28] N.S. Dhaidan, J.M. Khodadadi, T.A. Al-Hattab, S.M. Al-Masha, Experimental and numerical investigation of melting of phase change material/nanoparticle suspensions in a square container subjected to a constant heat flux, Int. J. Heat Mass Transf. 66 (2013) 672–683. [29] Z. Younsi, Etude expérimentale et numérique du comportement thermique de matériaux à changement de phase. Intégration dans un composant solaire passif pour l'habitat, Ph.D. Thesis Université d'Artois, 2008. [30] Z.X. Gong, A.S. Mujumdar, Flow and heat transfer in convection-dominated melting in a rectangular cavity heated from below, Int. J. Heat Mass Transf. 41 (17) (1998) 2573–2580. [31] R. Viswanath, Y. Jaluria, A comparison of different solution methodologies for melting and solidification problems in enclosures, Numer. Heat Transfer, Part B 24 (1) (1993) 77–105. [32] M. Lacroix, Computation of heat transfer during melting of a pure substance from an isothermal wall, Numer. Heat Transf., Part B 15 (1989) 191–210. [33] C.P. Dessai, K. Vafai, A unified examination of the melting process within a twodimensional rectangular cavity, ASME Trans. J. Heat Transfer. 115 (4) (1993) 1072–1075. [34] Z. Younsi, S. Lassue, L. Zalewski, D.R. Rousse, A. Joulin, A novel technique for the experimental thermophysical characterization of phase change materials, Int. J. Thermophys. 32 (2010) 674–692. [35] R. Viskanta, Phase change heat transfer, Solar Heat Storage: Latent Heat Materials, Edited by G.A. Lane, CRC Press, Boca Raton, Florida, 1983, pp. 153–222. [36] R. Viskanta, Natural convection melting and solidification, Natural Convection: Fundamentals and Applications, Hemisphere Publishing Company, Washington DC, 1985. [37] Y. Zhang, Z. Chen, Q. Wang, Q. Wu, Melting in an enclosure with discrete heating at constant rate, Exp. Thermal Fluid Sci. 6 (1993) 196–201. [38] P. Brousseau, M. Lacroix, Study of the thermal performance of a multi-layer PCM storage unit, J. Energy Convers. Manag. 37 (5) (1996) 599–609. [39] R. Viskanta, Heat transfer during melting and solidification of metals, J. Heat Transfer. Trans. ASME 110 (1988) 1205–1219. [40] P. Prescott, F.P. Incropera, Convection heat and mass transfer in alloy solidification, Adv. Heat Tran. 28 (1996) 231–338. [41] M. Lacroix, V.R. Voller, Finite difference solutions of solidification phase change problems: transformed vs. fixed grids, Numer. Heat Transfer 17 (1990) 25–41. [42] S.P. Ketkar, M. Parang, R.V. Arimilli, Experimental and numerical study of solidification and melting of pure materials, J. Thermophys. 5 (1991) 40–45. [43] R. Viswanath, Y. Jaluria, Numerical study of conjugate transient solidification in an enclosed region, Numer. Heat Transfer, Part A 27 (1995) 519–536. [44] M. Costa, A. Oliva, C.D.P. Segarra, R. Alba, Numerical simulation of solid-liquid phase change phenomena, Comput. Methods Appl. Mech. Eng. 91 (1991) 1123–1134. [45] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, 1980. [46] S. Kim, M.C. Kim, S.B. Lee, Prediction of melting process driven by conductionconvection in a cavity heated from the side, Korean J. Chem. Eng. 18 (5) (2001) 593–598.
promoted and this leads to an irregular morphology of the melting front. The time required for material melting is proportional to the Rayleigh number and the heat transfer by convection. For a fixed number of Rayleigh, the more the aspect ratio increases, the more the heat transfer is intense. In addition, heat transfer is higher with a large Rayleigh number, the aspect ratio being fixed. The significant rise in the thickness of the container induces the establishment of a thermal stratification and the liquid remains within the top portion of the enclosure. Under dimensions adopted herein, the best choice of the melting process is achieved in a cavity whose aspect ratio is 0.19 (H = 21 cm; W = 4 cm).
To sum up, our findings show that the efficiency of the macro-encapsulated PCM in a rectangular container can be improved. The incorporation of such a composite material within the latent heat thermal energy storage system can improve comfort, and the obtained results should be useful in helping to reduce the energy consumption in buildings. Conflict of interest statement The authors declare no conflict of interests regarding authorship and/or publication of this paper. Acknowledgments The authors are grateful to anonymous referees for their insightful review and rewarding comments that helped in improving the original manuscript. Likewise, we further thank the editor-in-chief for his continuing incentives. References [1] A. Sharma, V.V. Tyagi, C.R. Chen, D. Buddhi, Review on thermal energy storage with phase change materials and applications, Renew. Sust. Energ. Rev. 13 (2) (2009) 318–345. [2] B. Zalba, J.M. Marin, L.F. Cabeza, H. Mehling, Review on thermal energy storage with phase change: materials, heat transfer analysis and applications, Appl. Therm. Eng. 23 (3) (2003) 251–283. [3] M.M. Farid, A.M. Khudhair, S.A.K. Razack, S. Al-Hallaj, A review on phase change energy storage: materials and applications, Energy Convers. Manag. 45 (2004) 1597–1615. [4] R. Jacob, F. Bruno, Review on shell materials used in the encapsulation of phase change materials for high temperature thermal energy storage, Renew. Sust. Energ. Rev. 48 (2015) 79–87. [5] Y. Tian, C.Y. Zhao, A numerical investigation of heat transfer in phase change materials (PCMs) embedded in porous metals, Energy 36 (2011) 5539–5546. [6] B. Kamkari, H. Shokouhmand, F. Bruno, Experimental investigation of the effect of inclination angle on convection-driven melting of phase change material in a rectangular enclosure, Int. J. Heat Mass Transf. 72 (2014) 186–200. [7] Y. Tian, C.Y. Zhao, Thermal and exergetic analysis of metal foam enhanced cascaded thermal energy storage (MF-CTES), Int. J. Heat Mass Transf. 58 (2013) 86–96. [8] M. Jradi, M. Gillott, S. Riffat, Simulation of the transient behaviour of encapsulated organic and inorganic phase change materials for low-temperature energy storage, Appl. Therm. Eng. 59 (2013) 211–222. [9] A. Koca, H. Oztop, T. Koyun, Y. Varol, Energy and exergy analysis of a latent heat storage system with phase change material for a solar collector, Renew. Energy 33 (4) (2008) 567–574. [10] A. Khalifa, E. Abbas, A comparative performance study of some thermal storage materials used for solar space heating, Energ. Buildings 41 (2009) 407–415. [11] J.N.W. Chiu, V. Martin, Submerged finned heat exchanger latent heat storage design and its experimental verification, Appl. Energy 93 (2012) 507–516. [12] Z. Younsi, S. Lassue, L. Zalewski, D.R. Rousse, A. Joulin, A novel technique for the experimental thermophysical characterization of phase change materials, Int. J. Thermophys. 32 (2010) 674–692. [13] V.V. Tyagi, D. Buddhi, Thermal cycle testing of calcium chloride hexahydrate as a possible PCM for latent heat storage, Sol. Energy Mater. Sol. Cells 92 (2008) 891–899. [14] F. Agyenim, N. Hewitt, P. Eames, M. Smyth, A review of materials, heat transfer and phase change problem formulation for latent heat thermal energy storage systems
24