Contents lists available at ScienceDirect

Case Studies in Thermal Engineering

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

Melting inside a horizontal cylindrical capsule

M. Hlimia, S. Hamdaouib, M. Mahdaouic'*, T. Kousksoud, A. Ait Msaadb, A. Jamila, A. El Bouardib

a École Supérieure de Technologie, Université Sidi Mohamed Ibn Abdelah, Fès, Morocco b Laboratoire Énergétique, Université Abdelmalek Essaâdi, Tétouan, Morocco c LIMSI-CNRS, Université Paris Saclay, B.P. 133, 91403 Orsay, France

d Laboratoire des Sciences de l'Ingénieur Appliquées à la Mécanique et au Génie Electrique (SIAME), Université de Pau et des Pays de l'Adour - IFR - A. Jules Ferry, 64000 Pau, France

CrossMarlc

ARTICLE INFO

ABSTRACT

Keywords:

Melting

Natural convection Horizontal cylinder

The role of natural convection on solid-liquid interface motion during constrained melting within a horizontal cylindrical capsule was investigated. A numerical code is developed using an unstructured finite-volume method and an enthalpy porosity technique to solve for natural convection coupled to solid-liquid phase change. Flow patterns for different Rayleigh numbers are presented. The resulting melt shapes and the temperature in the PCM provide conclusive evidence of the importance of natural convection on heat transfer in the melt region.

1. Introduction

Heat transfer during melting and/or freezing of a phase change material (PCM) has attracted considerable attention over the past several decades due to its relevance to many technological applications such as latent-heat energy storage systems, solar energy systems, casting and crystal growth processes [1,2].

Divers papers have examined various aspects of melting and freezing phenomena both from the fundamental as well as applications point of view [3,4]. The non-linearity of the governing energy equation and a wide variety of geometric and thermal boundary conditions provide a fertile ground for challenging basic research problems [5,6]. Also, numerous industrial applications in diverse industries provide the necessary incentive for engineering research and development [7,8]. Time-dependent boundary conditions, under some conditions can lead to interesting and unique multiple moving boundaries as well.

Although a number of experimental and numerical studies have been devoted to convection dominated melting of a PCM for various geometric configurations, particular attention is given to melting in a horizontal cylinder as a model for thermal energy storage system [9-16]. A number of numerical and analytical studies have been performed in an attempt to model the melting phenomenon in a horizontal cylinder based on the Boussinesq approximation [17-23]. Saitoh and Hirose [17] presented experimental results which show convexed solid-liquid interface at the bottom of the unmelted solid during the melting process. Rieger et al. [18], Ho and Viskanta [19] and Yoo and Ro [20] investigated experimentally the evolution of the solid-liquid interface during melting of a PCM contained in a horizontal cylinder. They showed concaved interface numerically and experimentally.

Prasad and Sengupta [21]studied numerically the unconstrained melting of a PCM inside an isothermal horizontal tube. The model evaluated the irregular, temporal shape of the solid-liquid interface, the downward traveling motion of solid PCM due to density difference, and natural convection in the liquid phase.

A numerical investigation of melting of a phase change material inside a horizontal cylindrical capsule that is heated isothermally

* Corresponding author. E-mail address: mustapha_mahdaoui@yahoo.fr (M. Mahdaoui).

http://dx.doi.org/10.1016/j.csite.2016.10.001 Received 9 September 2016; Accepted 7 October 2016 Available online 08 October 2016

2214-157X/ © 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by/4.0/).

Nomenclature Greek symbols

c specific heat capacity (J kg-1 K-1) T the viscous stress tensor

f liquid fraction À thermal conductivity (W m-1 K-1)

~g the acceleration of gravity vector ( m s-2) r diffusion coefficient

Lf melting heat (J kg-1) ß dynamic viscosity (kg m-1 s-1)

Nu Nusselt number P density(kg m-3)

p pressure (Pa) ß the coefficient of volumetric thermal expansion

S surface (m2) (K-1)

T temperature (°C)

t time (s) Subscripts

~u velocity vector ( m s-1)

V control volume (m3) i initial

x coordinate (m) m melting

was presented by Park et al. [22]. An initial disturbance was introduced for the sake of examining its effect on the solution for two values of the Rayleigh number. The results exposed that at Ra=106, the flow was unicellular, but when the Rayleigh number was increased to 8.106, two branching solutions (bi-cellular and tri-cellular flow) were obtained according to the type of initial disturbances. Ro and Kim's [23] showed, by using the numerical analysis, that both interface shapes are possible.

Chung et al. [24] analyzed numerically by the enthalpy method the multi-cell structure and the thermal instability at the early stage of the melting process in a horizontal cylinder heated isothermally for a relatively wide range of the Rayleigh numbers. The results exposed that at the low Rayleigh number, the flow in the liquid gap is in the stable state as the viscous force dominates. The intermediate Rayleigh number brought about the onset of the Benard convection that depends on the strength of the base flow under development and interaction between the two flows determines the flow pattern afterwards. At the high Rayleigh number, the Benard convection shows an orderly behavior without being affected by the base flow.

It is interesting to note that few report [23] for the melting process inside a horizontal cylinder for the high Rayleigh number has been found even though the thermal instability at the initial stage depends strongly on the Rayleigh number.

This study aims to investigate numerically the natural convection dominated melting of a PCM within a horizontal cylindrical capsule for Rayleigh numbers in the range 104<Ra<108. The enthalpy-porosity method is employed where the interface solid-liquid is implicitly determined in the calculation domain. A numerical code is developed using an unstructured finite-volume method and to solve for natural convection coupled to solid-liquid phase change.

2. Physical model and basic equations

As depicted in Fig. 1, the problem considered is the melting process of a PCM at initial temperature of Ti in a horizontal cylinder whose wall temperature is Tw above the melting temperature Tm. The gallium is chosen as the test PCM. The thermophysical properties are listed in Table 1. Unsteady two-dimensional melting of PCM is governed by the basic laws represented by the

Fig. 1. Schematic diagram of the inward melting process in a horizontal cylinder.

Table 1

Physical properties of pure gallium.

Density (liquid) Reference temperature

Volumetric thermal expansion coefficient of liquid

Thermal conductivity

Melting temperature (Tm)

Latent heat of fusion

Specific heat capacity

Dynamic viscosity

Prandtl number

6093 kg m-3 29.78 °C 1.2 10-4K-1 32 W m-1 K-1 29.78 °C 80,160 J kg-1 381.5 J kg-1

1.81 10-3 kg m-1 s-1 2.16 10-2

Fig. 2. Melting front positions at several instants.

400 600

Time (s)

Fig. 3. Liquid fraction versus time. continuity, momentum and energy equations and by the following assumptions:

- The thermophysical properties of the PCM are constant but may be different for the liquid and solid phases.

Fo = 2.20

Fo= 11.013

Fo= 16.519

Fo= 27.533

Fig. 4. Stream lines, temperatures contours and interfaces liquid-solid at various times (Ra=3-104).

- The Boussinesq approximation is valid, i.e., liquid density variations arise only in the buoyancy source term, but are otherwise neglected.

- The melting of PCM will be symmetric.

- The liquid is Newtonian and incompressible.

Fo= 0.91

Fo= 2.73

Fo= 5.461

Fo= 7.281

Fig. 5. Stream lines, temperatures contours and interfaces liquid-solid at various times (Ra=5-106). The flow is laminar and two-dimensional.

Then the conservation equations in their respective integral forms are:

Fo= 0.072

Fo= 0.182

Fo= 0.364

Fo= 0.782

Fig. 6. Stream lines, temperatures contours and interfaces liquid-solid at various times (Ra=4-107).

/ ~u . J S

' dS = 0

Nusselt

Fig. 7. Local Nusselt number at various times (Ra=3-104).

10 15 20 25

Nusselt

Fig. 8. Local Nusselt number at various times (Ra=5-106).

— f P U dV + f P ~u~u . ndS = - f V p dV + f T .ndS + f~Au dV

dt J v J s Jv J S Jv

~ f P cpT dV + f p cpT U.~ndS = f ÀV T.~ndS + f pLf —

dt Jv p Js p Js Jv F dt

where ~u is the velocity vector, p the pressure and T the temperature. t is the viscous stress tensor for a Newtonien fluid:

T = fi (Vu + (Vu)T)

The integration occurs over a control volume V surrounded by a surface S, which is oriented by an outward unit normal vector ~n . The source term in Eq. (2) contains two parts:

Case Studies in Thermal Engineering 8 (2016) 359-369 40

Nusselt

Fig. 9. Local Nusselt number at various times (Ra=4-107).

4000 6000

Time (s)

Fig. 10. Liquid fraction for two Rayleigh numbers.

Au = pP (T - Tm)J + A S

where p is the coefficient of volumetric thermal expansion and g the acceleration of gravity vector. The first part of the term source represents the buoyancy forces due to the thermal dilatation. Tm is the melting temperature of the PCM. The last term is added to account for the velocity switch-off required in the solid region. During the solution process of the momentum field, the velocity at the computational cell located in the solid phase should be suppressed while the velocities in the liquid phase remain unaffected. Also, as the solid region melts the mass in the corresponding computational cell should begin to move. The present study adopts a Darcy-like momentum source term to simulate the velocity switch-off [24,25]

C (1 - f )2

f3 + £

The constant C has a large value to suppress the velocity at the cell becomes solid and £ is small number used to prevent the division by zero when a cell is fully located in the solid region, namely f = 0. In this work, C = 1 x 1015 and £ = 0.001 are used.

3. Numerical procedure

The conservation Eqs. (1)-(3) are solved by implementing them in an in house code. This code has been successfully validated in several situations involving flow and heat transfer as in [26-29]. The present code has a two dimensional unstructured finite-volume framework that is applied to hybrid meshes. The variables values are stored in cell centers in a collocated arrangement. All the conservation equations have the same general form. By taking into account the shape of control volumes, the representative

conservation equation to be discretized may be written as

- f py dV + ^ f p uy — rdSi = f SvdV

- Jy ^JSi dxi Jv

dt ■

which has four distinctive parts: rate of change, convection, diffusion and source (the mass conservation equation, however, does not have a diffusion term). The rate-of-change and source terms were integrated over the cell volume, whereas the convection and diffusion terms formed the sum of fluxes through the control volume faces.

Solidification and melting are generally transient phenomena, where the explicit schemes are too restrictive owing to stability limitations. Hence implicit schemes are often preferred and the simplest choice is the first order Euler scheme. The cell face values, appearing in the convective fluxes, were obtained by blending the UDS (upwind differencing scheme) and the CDS (central difference scheme) using the differed correction approach [26]. The coupling of the dependent variables was obtained on the basis of the iterative SIMPLE algorithm [30,31].

Summation of the fluxes through all the faces of a given CV (Control Volume) results in an algebraic equation which links the value of the dependent variable at the CV center with the neighboring values. The equation may also be written in a conventional manner as

apVp = ^ A"bVnb + bv

nb (8)

The coefficients A^ contain contributions of the neighboring CVs, nb, arising out of convection and diffusion fluxes as defined by Eqs. (1)-(3). The central coefficient AP on the other hand, includes the contributions from all the neighbours and the transient term. In some of the cases, where sources term linearization was applied, it also contained part of the source terms. b^contains all the terms those are treated as known (source terms, differed corrections and part of the unsteady term).

The evaluation of the source term in the energy equation has been made using the new source algorithm proposed by Voller [32] where the new value of f at each outer iteration n+1 in cell P is calculated as follows:

fP + PALV (TP — Tm) (9)

This update is followed by an overshoot/undershoot correction:

,n+1_ J° f fP+l < 0

fn+1 =

P " \l if fn+l > 1 (10)

Inner iterations are stopped when the 1-norm of residual is reduced by a factor set to 10-8 for u and energy equations, and to 10-6 for the pressure correction equation. The convergence criterion for outer iterations is based on a total number of outer iterations rather than a prescribed value for the residual. The choice of total number of outer iterations is based on the monitoring of the outer iteration errors (ratio of the norms of the difference between two consecutives iterates — ifin\ and the difference between the first two iterates. The time steps used for the calculations vary from At = 1 s down to 0.001 s.

4. The code validation

The present code has been successfully validated in several situations involving flow and heat transfer as in [26-29]. We focus here on the validation of the code in the case of the melting of a pure PCM coupled to the natural convection in the melt. The chosen test case is the two-dimensional numerical benchmark presented by Hannoun et al. [33] that concerns the melting of tin in a square cavity subject at one side to a temperature higher than the melting temperature. The authors have presented extensive results obtained with second order accurate finite volume method, in structured meshes as fine as 600x600. We carry out a numerical simulation in the same conditions as in [33] and with the same material properties, using a uniform grid 350x350 of size. In Fig. 2 we plot the melting front at different times.

We present also in Fig. 3 the total liquid fraction in the whole square cavity. All of these results are compared to those of [33] and exhibit a satisfactory correspondence.

5. Results and discussion

Using the above-described model, simulations were carried out for melting of a PCM (Gallium) within a horizontal cylinder. Numerical investigations were conducted using 85,000 cells and the time step of 10-3 s was found to be sufficient to give accurate results.

The effect of the melting process on the streamlines, temperature contours and the solid-liquid interface for different Fourier number is gauged through the result illustrated in Fig. 4. During the early phase the liquid is confined between the rigid heated cylinder and a concentric moving solid-liquid interface. Inspection of the figures reveals that at early times the melt regions are similar in shape and that when heat transfer to the gallium is predominantly by conduction the melt region is symmetrical about the axis of the cylinder. After some time natural convection develops and intensifies, influencing the melt shape in general and the melt region above the cylinder in particular. Natural convection conveys the hot liquid to the upper part of the melt zone and in this

manner continues to support the upward movement of the solid-liquid interface.

At low Rayleigh number (Ra=3-104) the solid-liquid interface is approximately concentric inside the cylinder (see Fig. 4). We can also note that during the melting process the flow changes into a totally different pattern as in Fig. 4 which shows the flow transition from the base flow (single-cell) to the three-cell flow at the intermediate stage and finally to the two cell flow as time elapses.

Figs. 5 and 6 display the flow behavior during the melting process for Rayleigh numbers of 5-106 and 4-107 respectively. After the conduction dominating stage, a complex structure of the fluid dynamics characterized by a multi-cellular flow patterns. The number and the size of the cells are dependent on the size of the molten phase. At Fo=5.461, four counter-rotating cells are captured for Ra=5-106 (see Fig. 5). As the liquid layer increases, the number of rolls decreases. So at Fo=7.281, when the liquid fills nearly the half of the cylinder, two counter-rotating rolls are obtained. This phenomenon is very similar to the Benard convection due to the thermal instability because the solid-liquid interface and the cylinder bottom wall can be approximated as two flat plates as far as the liquid gap between these two surfaces is narrow enough. This phenomenon leads then in the sequel of calculations to the development of the abovementioned roll cells. The establishment of these vortices influences the melting kinetics of the whole system. Especially the heat transfer in the lower part of the melting gap is greatly improved which results in a "moonshaped" melting contour. We can also observe that the flow in the liquid zone remains in the stable state at low Rayleigh number. At the high Rayleigh number, the Benard convection presents an orderly behavior without been affected by the base flow.

The local Nusselt number at the hot walls is an ideal indicator showing how convection influences overall conduction through the cylindrical capsule:

h. Lref

NU =— (11) where Lref is the reference length, X is the thermal conductivity and h is the heat transfer coefficient:

Tw - Tref (12)

q is the convective heat flux, Tw is the wall temperature, and Tref is the reference temperature.

The effect of the Rayleigh number on the local Nusselt number at various angular positions along the cylinder surface for several Fourier numbers is shown in Figs. 7, 8 and 9. For low Rayleigh number (Ra=3-104), the Nusselt number is approximately constant at the early stage of the melting process. This result confirms that at early times, heat transfer in the melt zone predominantly by conduction. Sudden drops of the Nusselt number appear at times corresponding to roll pairings. The drops are more pronounced for Ra=5-106and Ra=4-107.These results can be explained by the fact that the higher Rayleigh number is associated with strong natural convection velocities resulting in higher heat transfer rates along the cylinder surface.

It is interesting to note that the rate of melting is one of the most important parameters of the problem. The time evolution of the total liquid fraction (ratio of volume of melt to volume cavity) is a factor that has been widely used as a monitoring parameter in earlier publications.

From the liquid fraction versus time plot, one can get both the rate melting (slope of the tangent line at a given time) and the average melting rate (ratio of current liquid fraction and time). Fig. 10 shows the time evolution of the total fraction of the liquid in the cylindrical capsule for two Rayleigh numbers. We note that as Rayleigh number increases the PCM inside the cylinder melts faster. We can also note that, once the convection is established, the variation of the liquid fraction with time seems to be linear.

6. Conclusion

This paper investigates numerically the effect of the Rayleigh number on the natural convection during the melting process within a horizontal cylinder. It is found that the melting process of PCM is dominated by conduction at early stages indicated by concentric solid-liquid interfaces parallel to the circular heated wall. At later times, natural convection is augmented and this affects the shape of the melting front and consequently results in unsymmetrical melting about the horizontal axis of the capsule. The numerical results indicate also that for low Rayleigh number, the flow in the melting zone remains in the stable state as the viscous force is dominant. At the high Rayleigh number, the Rayleigh-Benard convection indicates an orderly behavior without been affected by the base flow.

References

[1] T. Kousksou, P. Bruel, A. Jamil, T. El Rhafiki, Y. Zeraouli, Energy storage: applications and challenges, Sol. Energy Mater. Sol. Cells 120 (2014) 59-80.

[2] A. Sharma, V.V. Tyagi, C.R. Chen, D. Buddhi, Review on thermal energy storage with phase change materials and applications, Renew. Sustain. Energy Rev. 13 (2009) 318-345.

[3] J.M. Khodadadi, Y. Zhang, Effects of buoyancy-driven convection on melting within spherical containers, Int. J. Heat Mass Transf. 44 (2001) 1605-1618.

[4] A. Felix Regin, S.C. Solanki, J.S. Saini, Heat transfer characteristics of thermal energy storage system using PCM capsules: a review, Renew. Sustain. Energy Rev. 12 (2008) 2438-2458.

[5] A. Arid, T. Kousksou, S. Jegadheeswaran, A. Jamil, Y. Zeraouli, Numerical simulation of ice melting near the density inversion point under periodic thermal boundary conditions, Fluid Dyn. Mater. Process. 305 (2012) 1-19.

[6] T. Kousksou, A. Arid, J. Majid, Y. Zeraouli, Numerical modeling of double-diffusive convection in ice slurry storage tank, Int. J. Refrig. 33 (2010) 1550-1558.

[7] Y.H. Wan, Y.T. Yang, Three-dimensional transient cooling simulations of a portable electronic device using PCM (phase change materials) in multi-fin heat sink, Energy 36 (2011) 5214-5224.

[8] W.L. Cheng WL, B.J. Mei, Y.N. Liu, Y.H. Huang, X.D. Yuan, A novel household refrigerator with shape-stabilized PCM (phase change material) heat storage condensers: an experimental investigation, Energy 36 (2011) 5797-5804.

[9] S. Kalaiselvam, M. Veerappan, A. Arul Aaron, S. Iniyan, Experimental and analytical investigation of solidification and melting characteristics of PCMs inside cylindrical encapsulation, Int. J. Therm. Sci. 47 (2008) 858-874.

[10] N.S. Dhaidan, J.M. Khodadadi, T.A. Al-Hattab, S.M. Al-Mashat, Experimental and numerical study of constrained melting ofn-octadecane with CuO nanoparticle dispersions in a horizontalcylindrical capsule subjected to a constant heat flux, Int. J. Heat Mass Transf. 67 (2013) 523-534.

[11] E. Assis, L. Katsman, G. Ziskind, R. Letan, Numerical and experimental study of melting in a spherical shell, Int. J. Heat Mass Transf. 50 (2007) 1790-1804.

[12] F.L. Tan, S.F. Hosseinizadeh, J.M. Khodadadi, Liwu Fan, Experimental and computational study of constrained melting of phase change materials (PCM) inside a spherical capsule, Int. J. Heat Mass Transf. 52 (2009) 3464-3472.

[13] D.B. Khillarkar, Melting of Phase Change Material in Horizontal Annuli (Thesis), McGill University Montreal, Quebec, Canada, 1998.

[14] A. Felix Regin, S.C. Solanki, J.S. Saini, Latent heat thermal energy storage using cylindrical capsule: numerical and experimental investigations, Renew. Energy 31 (2006) 2025-2041.

[15] T. Kousksou, M. Mahdaoui, A. Ahmed, J. Batina, Numerical simulation of PCM melting over a wavy surface, Int. J. Numer. Methods Heat Fluid Flow 24 (2014) 1660-1669.

[16] K. El Omari, T. Kousksou, Y. Le Guer, Impact of shape of container on natural convection and melting inside enclosures used for passive cooling of electronic devices, Appl. Therm. Eng. 31 (2011) 3022-3035.

[17] T. Saitoh, K. Hirose, High Rayleigh number solutions to problems of latent heat thermal energy storage in a horizontal cylinder capsule, ASME J. Heat Transf. 104 (1982) 545-553.

[18] H. Rieger, U. Projahn, M. Bareiss, H. Beer, Heat transfer during melting inside a horizontal tube, ASME J. Heat Transf. 105 (1983) 226-234.

[19] C.J. Ho, R. Viskanta, Heat transfer during inward melting in a horizontal tube, Int. J. Heat Mass Transf. 27 (1984) 705-716.

[20] H.S. Yoo, S.T. Ro, Numerical analysis of the phase change processes by coordinate transformations, Trans. KSME 10 (1986) 585-592.

[21] A. Prasad, S. Sengupta, Nusselt number and melt time correlations for melting inside a horizontal cylinder subjected to an isothermal wall temperature condition, ASME J. Heat Transf. 110 (1988) 340-345.

[22] C.E. Park, S.S. Kim, K.S. Chang, Branching solutions to inward melting problem in a horizontal tube, Int. Commun. Heat Mass Transf. 18 (1991) 343-350.

[23] S.T. Ro, C.J. Kim, Bifurcation phenomenon during the fixed solid mode melting inside a horizontal cylinder, Int. J. Heat Mass Transf. 37 (1994) 1101-1109.

[24] J.D. Chung, J.S. Lee, Thermal instability during the melting process in an isothermally heated horizontal cylinder, Int. J. Heat Mass Transf. 40 (1997) 3899-3907.

[25] B.J. Jones, D. Sun, S. Krishnan, S.V. Garimella, Experimental and numerical study of melting in a cylinder, Int. J. Heat Mass Transf. 49 (2006) 2724-2738.

[26] T. Kousksou, M. Mahdaoui, A. Ahmed, A. Ait Msaad, Melting over a wavy surface in a rectangular cavity heated from below, Energy 64 (2014) 212-219.

[27] M. Mahdaoui, T. Kousksou, S. Blancher, A. Ait Msaad, T. El Rhafiki, M. Mouqallid, A numerical analysis of solid-liquid phase change heat transfer around a horizontal cylinder, Appl. Math. Model. 38 (2014) 1101-1110.

[28] K. Tarik, M. Mustapha, H. Mohamed, J. Abdelmajid, T. El Rhafiki, Heat transfer from horizontal cylinder with fins embedded in PCM, in: Proceedings of the 2014 International Renewable and Sustainable Energy Conference (IRSEC), Ouarzazate, 2014, pp. 442-446

[29] K. Tarik, M. Mustapha, H. Mohamed, J. Abdelmajid, T. El Rhafiki, Heat transfer from multiple cylinders during melting process: solar applications, in: Proceedings of the 2014 International Renewable and Sustainable Energy Conference (IRSEC), Ouarzazate, 2014, pp. 141-146

[30] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington DC, 1980.

[31] S.V. Patankar, D.B. Spalding, A calculation of heat, mass and momentum transfer in three dimensional parabolic flows, Int. J. Heat Mass Transf. 15 (1972) 1787-1806.

[32] V. Voller, Fast implicit finite-difference method for the analysis of phase change problems, Numer. Heat Transf. Part B - Fundam. 17 (1990) 155-169.

[33] N. Hannoun, V. Alexiades, T.Z. Mai, A reference solution for phase change with convection, Int. J. Numer. Methods Fluids 48 (2005) 1283-1308.