Melting over a wavy surface in a rectangular cavity heated from below

Energy 64 (2014) 212e219

Contents lists available at ScienceDirect

Energy journal homepage: www.elsevier.com/locate/energy

Melting over a wavy surface in a rectangular cavity heated from below T. Kousksou a, *, M. Mahdaoui a, A. Ahmed a, A. Ait Msaad b a 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 b Ecole Nationale Supérieure d’Arts et Métiers, ENSAM Marjane II, BP-4024 Meknès Ismailia, Morocco

a r t i c l e i n f o

a b s t r a c t

Article history: Received 2 March 2013 Received in revised form 7 November 2013 Accepted 11 November 2013 Available online 15 December 2013

The current numerical study is conducted to analyze melting in a rectangular closed enclosure by subjecting the bottom wavy surface to a uniform temperature. The cavity vertical walls and the top wall are insulated while the bottom wall is maintained at temperature TB ¼ 38.3  C. The enclosure was filled by solid Gallium initially at temperature Ti ¼ 28.3  C. A numerical code is developed using an unstructured finite-volume method and an enthalpy porosity technique to solve for natural convection coupled to solideliquid phase change. The validity of the numerical code used is ascertained by comparing our results with previously published results. The effect of the amplitude of the wavy surface on the flow structure and heat transfer characteristics is investigated in detail. It is found that the rate of the melting increases with the elevation in the magnitude of the amplitude value of the wavy surface. Ó 2013 Elsevier Ltd. All rights reserved.

Keywords: Phase change material Wavy surface Natural convection Heat transfer

1. Introduction The problem of solideliquid phase change with natural convection in the liquid phase has continued to be an important research area for several decades, due to its application in many technologically important processes such as metal casting, food freeze drying, electronics cooling, and in particular, thermal energy storage [1e7]. The solideliquid transition provides also a suitable technique for controlling temperature in systems, which are subjected to periodic heating. This process allows, for periodic heating, the conversion of temperature oscillations into oscillations of the melting interface, with a significant damping of the perturbation. Furthermore, the energy stored during melting can be recovered during freezing, with significant energetic opportunities. Although many experimental and numerical studies have been dedicated to convection-dominated melting of PCMs (Phase Change Materials) for various geometry configurations, e.g., along a vertical wall, inside as well as around a horizontal cylinder, etc. [8e 11], little effort has been reported on the melting of a PCM heated from below [8,10,12e17]. It is well known that when a horizontal layer of fluid is heated from below, a cellular form of natural convection may occur inside the liquid, similar to that observed during classical RayleigheBénard convection. As a result of the interaction between the cellular convection and the melting process, the phase-change interface, which is assumed to be at the equilibrium * Corresponding author. E-mail address: [email protected] (T. Kousksou). 0360-5442/$ e see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.energy.2013.11.033

temperature, tends to get wavy. For the melting from below, few works exist in the literature. Yen and Galea [17,18] and Seki et al. [19] experimentally studied the melting of a horizontal ice slab heated from below. Hale and Viskanta [20] performed experiments for melting from below and solidification from top of n-octadecane in a rectangular cavity. They did not present flow patterns and phase change interface shapes, however, Gau et al. [21] presented flow visualization for melting from below of an n-octadecane slab in a rectangular cavity. Diaz and Viskanta [22] extended the experiments of Gau et al. [21] to morphology observation of the liquidesolid interface. Prud’Homme et al. [14] studied numerically the melting of a pure substance within a vertical cylindrical enclosure heated isothermally from below and assuming adiabatic conditions at the vertical side. The governing equations for the convective flow in the melt are solved using computer generated body-fitted curvilinear coordinates. In their study, the authors found that the critical Rayleigh number for the onset of convection based on the melt layer thickness is 2197. The flow patterns show that the initial multicellular regime is replaced ultimately by a single Bénard cell. Lacroix [12] performed a series of numerical simulations using a variable grid method for the melting of a low Prandtl substance. The author analyzed the effect of the temperature difference between the bottom wall and the initial temperature of the solid. The effect of the Rayleigh number, on the flow patterns during the melting of a pure phase change material in square cavity isothermally heated from below and which is adiabatic at the rest of the walls, is studied numerically by Gong and Mujumdar [13]. The governing equations for the convective flow in

T. Kousksou et al. / Energy 64 (2014) 212e219

Nomenclature a AP b c f ! g LF L Nu p Pr Ra S Ste T t ! u V

amplitude of the wavy surface (m) linearized coefficient number of undulation specific heat capacity (J kg1 K1) liquid fraction acceleration of gravity vector (m s2) melting heat (J kg1) length (m) Nusselt number pressure (Pa) Prandtl number Rayleigh number heat source (W m3) Stefan number temperature ( C) time (s) Vitesse vector (m s1) control volume (m3)

the melt are solved using the Streamline Upwind/Petrov Galerkin finite element method in combination with a fixed grid primitive variable method. More recently, Fteiti and Nasrallah [10] investigated the melting process in rectangular enclosure heated from below and cooled from above. The authors used the enthalpye porosity method to track the motion of the solideliquid interface and the governing equations are numerically solved using the control volume based finite element method. They concluded that the melting process is more rapid as the aspect ratio of the enclosure is small, but the volume fraction melted at the steady state is less important. It is interesting to note that, flow and heat transfer from irregular surfaces are often encountered in many engineering applications to enhance heat transfer such as micro-electronic devices, flat plate solar collectors and flat-plate condensers in refrigerators, etc. Roughened surfaces could be used in latent storage systems where the wall heat flux is known. One of the reasons why a roughened surface is more efficient in heat transfer is its capability to promote fluid motion near the surface; in this way a complex wavy surface is expected to promote a larger heat transfer rate than a flat plate. This complex geometry will promote a correspondingly complicated motion in the fluid near the surface; this motion is described by the nonlinear boundary-layer equations. A vast amount of literature about convection along a sinusoidal wavy surface is available for different heating conditions and various kinds of fluids [23e25]. The major conclusion is that the total heat-transfer rate for a wavy surface of any kind is, in general, greater than that of the corresponding flat plate, and may be a function of the ratios of amplitudes to wavelengths of the surface. On the other hand, the result that the average heat-transfer rate per unit of wetted surface for a wavy surface is less than that of a flat plate is confirmed. The total heat-transfer rate is the more important factor in designing a heat-transfer surface. To the best of the authors’ knowledge, no attention has been paid to the problem of the heat transfer during melting in a rectangular cavity that is heated from a horizontal wavy surface. The objective of the present study is to examine the melting process in a rectangular cavity by subjecting the bottom wavy surface to a relatively higher temperature than the melting temperature of the PCM. We aim to investigate the impact of the wavy surface amplitude on the kinetics melting of the PCM inside the enclosure. The results are shown in terms of parametric presentations of streamlines and isotherms.

W x

213

width (m) coordinate (m)

Greek symbols coefficient of volumetric thermal expansion (K1) thermal conductivity (W m1 K1) transport variable diffusion coefficient dynamic viscosity (kg m1 s1) density (kg m3) stream function (m2 s1) under-relaxation factor

b l f G m r j u

Subscripts C cold H hot i initial m melting nb neighboring

2. Mathematical formulation 2.1. Physical model and basic equations Unsteady two-dimensional melting of PCM is governed by the basic laws represented by the 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. - The Boussinesq approximation is valid, i.e., liquid density variations arise only in the buoyancy source term, but are otherwise neglected. - The liquid is Newtonian. - Viscous dissipation is neglected. - Fluid motion in the melt is laminar and two-dimensional. Since the present formulation deals with solutions on unstructured grids, it is essential to represent the conservation laws in their respective integral forms.

Z

!! u : n dS ¼ 0

(1)

S

d dt

Z

r! u dV þ

Z

V

!! r! u u : n dS ¼ 

Z

/

V

S

Z

V p dV þ

Z þ

s :! n dS

S

! B u dV

V

(2) d dt

Z

r cp T dV þ

V

Z

! r cp T ! u : n dS ¼

S

Z S

/

l V T:! n dS þ

Z V

rLF

vf vt (3)

! where u is the velocity vector, p the pressure and T the temperature. s is the viscous stress tensor for a Newtonian fluid :





s ¼ m Vu þ Vu

T 

(4)

214

T. Kousksou et al. / Energy 64 (2014) 212e219

Z

X Z 

 Z vf dSi ¼ Sf dV vxi

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:

d dt

! ! ! B u ¼ rbðT  Tm Þ g þ B u

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 [8,13]. The coupling of the dependent variables was obtained on the basis of the iterative SIMPLE algorithm [10,14]. 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

(5)

! where b 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 [26]:

Cð1  f Þ2 B ¼  3 f þε

(6)

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  1015 kg/m3 s and ε ¼ 0:001 are used. The shape of the horizontal wavy surface profile is assumed to mimic the following pattern

Y ¼ a½1  cosð2bpXÞ

(7)

where a is the amplitude of the wavy surface and b is the number of undulation (see Fig. 1).

2.2. Numerical procedure The conservation Eqs. (1)e(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 Refs. [1,2]. The present code has a two dimensional unstructured finitevolume 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

rfdV þ

rui f  Gf

i

V

AP fP ¼

X

Sj

Anb fnb þ bf

(8)

V

(9)

nb

The coefficients Anb contain contributions of the neighboring (nb) CVs, arising out of convection and diffusion fluxes as defined by Eqs. (1)e(3). The central coefficient AP on the other hand, includes the contributions from all the neighbors and the transient term. In some of the cases, where sources term linearization was applied, it also contained part of the source terms. b4 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 [27] where the new value of f at each outer iteration n þ 1 in cell P is calculated as follows:

fPnþ1 ¼ fPn þ

AnP Dt ðT  Tm Þ rLF V P

(10)

This update is followed by an overshoot/undershoot correction:

( fPnþ1

¼

0 1

if fPnþ1 < 0 if fPnþ1 > 1

(11)

Inner iterations are stopped when the 1-norm of residual is ! reduced by a factor set to 108 for u and energy equations, and to 106 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 4nþ1  4n  and the difference between the first two iterates. The time steps used for the calculations vary from Dt ¼ 1 s down to 0.01 s. Under-relaxation factors are used for outer iterations: uu ¼ 0:65,uT ¼ 0:8, uP ¼ 0:9 where subscripts refer to the corresponding equation’s main variable.

3. Test of the numerical method

Fig. 1. Layout of the physical model.

In view to test and validate the numerical accuracy of the present modeling of a PCM melting process heated from below the

T. Kousksou et al. / Energy 64 (2014) 212e219

215

Table 1 Physical properties of pure gallium. Density (liquid) Reference density Reference temperature thermal expansion coefficient of liquid Thermal conductivity Melting temperature Latent heat of fusion Specific heat capacity Dynamic viscosity Prandtl number

6093 kg m3 6095 kg m3 29.78  C 1.2 104 K1 32 W m1 K1 29.78  C 80,160 J kg1 381.5 J kg1 1.81 103 kg m1 s1 2.16 102

Fig. 3. Streamlines ((a), (b), (A), (B)) and corresponding isotherms at t ¼ 342.5 s. Small letters: present model; capital letters: Gong model.

Fig. 2. Numerical and experimental melting front position.

numerical procedure is submitted to the two tests described in the following analysis.

3.1. Melting of a pure metal in Gau and Viskanta study [28] Before conducting such numerical investigations, the physical model was first used for validation against experimental data obtained by Gau and Viskanta’s [28], for the two-dimensional problem of melting of PCM with the presence of natural convection in a rectangular enclosure. The dimension of the cavity used in the experiments is L ¼ 6.35 cm and W ¼ 8.89 cm. The left hot wall and the right cold wall are maintained at temperatures, TH ¼ 38.35  C and TC ¼ 28.3  C, respectively, and it was filled by solid PCM initially at temperature TC. The horizontal walls are insulated. The physical properties for pure gallium used in this investigation are given in Table 1. The physical property values were chosen at 32  C, a temperature that is representative of the temperature range of the experiment. The Prandtl number for liquid gallium at this temperature is 0.0216 and the situation shown in the rectangular cavity corresponds to a Stefan number ðSte ¼ cðTH  Tm Þ=LF Þ of 0.039 and

a Rayleigh number ðRa ¼ r2 cgbL3 ðTH  Tm Þ=mlÞ of 6.105. Once the calculation is started, the solid gallium melts. Fig. 2 shows the shape and location of the solideliquid interface at several times during the melting process. The black and red lines indicate the experimental [28] and calculated data respectively. Before a time of 2 min, the shape of the interface is nearly flat because convection is still weak and melting is driven by conduction. After 2 min, the interface becomes wavy due to the circular flow inside the molten region. The position of the melt front near the top surface in the calculation before a time of 17 min is over-estimated compared to experiment, and after 10 min, it is underestimated. However, the overall trend shows good agreement with experiment and we can

Table 2 Physical properties and non-dimensional parameters. Aspect ratio (A ¼ L/W) Rayleigh number Prandtl number Stefan number Bottom temperature Melting temperature Initial temperature Ratio of solid/liquid specific heat capacity Ratio of solid/liquid conductivity

1 2.844  106 46.1 0.138 29.96  C 17  C 16.7  C 0.964 2.419

Fig. 4. Streamlines ((a), (b), (A), (B)) and corresponding isotherms at t ¼ 420 s. Small letters: present model; capital letters: Gong model.

216

T. Kousksou et al. / Energy 64 (2014) 212e219

safely said that our model and code are fairly validated to track the melting boundary satisfactorily. The discordance observed may be explained by the difficulty encountered by Gau and Viskanta’s [28] to ensure the stability in vertical walls temperatures during the experiment. The same remarks were made by other authors when validating their numerical code. 3.2. Melting of a paraffin in Gong and Mujumdar study [13]

Fig. 5. Local dimensionless heat flux distribution along the heated surface at different time steps.

In the second problem used to test the accuracy of our numerical procedures, we compare numerical results obtained by Gong and Mujumdar [13]. It concerns the melting of a PCM (n-octadecane) in a rectangular cell heated from below. The top and the two vertical walls are assumed to be adiabatic. The physical properties and nondimensional parameters used in this investigation are given in Table 2. Numerical investigations were conducted using 23,000 cells and the time step of 103 s was found to be sufficient to give accurate results. Figs. 3 and 4 show the computed streamlines and isotherms at different time steps. In the conductive zone (t < 200 s), our results are similar to those obtained by Gong and Mujumdar [13]. However, in the convective regime, the present model provides the same results as Gong model, but at time

Fig. 6. -a: Stream function and temperatures contours at time 100 s for various amplitudes. 6-b: Stream function and temperatures contours at time 200 s for various amplitudes. 6c: Stream function and temperatures contours at time 400 s for various amplitudes. 6-d: Stream function and temperatures contours at time 600 s for various amplitudes. 6-e: Stream function and temperatures contours at time 1000 s for various amplitudes.

T. Kousksou et al. / Energy 64 (2014) 212e219

sections slightly shifted. For example, streamlines and isotherms obtained by Gong for t ¼ 342.5 s are identical to those obtained by the present model at time t ¼ 420 s. The temporal evolution of the dimensionless heat flux distribution along the heated wall obtained by the two models at different time steps is almost identical, as shown in Fig. 5. 4. Results and discussion The characteristics of the flow and temperature fields in the rectangular cavity are examined by exploring the effects of the amplitude a of the wavy vertical surface on the melting process. In the current numerical investigation, the dimension of the cavity used is L ¼ 6 cm and W ¼ 6 cm and it was filled by solid PCM (Gallium) initially at temperature Ti ¼ 28.3  C. The top and the two vertical walls are assumed to be adiabatic. The bottom wall is maintained at temperature TB ¼ 38.3  C. The following parametric domains of the wavy surface are considered: 0 mm  a  14 mm. It is interesting to note that for melting from below and over a wavy surface, the Rayleigh number calculated using the cavity height is

217

not the true Rayleigh number. The true Rayleigh number should be calculated using the actual melt height instead of the height of the cavity. The effect of the amplitude on the streamlines and temperature contours for several times is gauged through the result illustrated in Fig. 6. The number of undulation is maintained at b ¼ 2. It appears from these results that varying the amplitude between 0 mm and 14 mm disturbs the heat transfer in the vicinity of the hot wall. At early times, heat transfer in the melt zone predominantly by conduction and the liquidesolid interface mimics the wall’s profile. However, as melting progresses, the cavity expands and convection takes over conduction. Note that, in the present study three-dimensional convection is neglected since a two dimensional model is used. However, the duration of the three dimensional convection is very short compared to the whole melting process. After the conduction dominating stage, a complex structure of the fluid dynamic is characterized by a multi-cellular flow patterns. The number and the size of the Bénard cells are dependent on the height of the molten phase.

Fig. 6. (continued).

218

T. Kousksou et al. / Energy 64 (2014) 212e219

Fig. 9. Liquid fraction versus time for various amplitudes. Fig. 7. Nusselt number at time 200 s for various amplitudes.

At t ¼ 400 s, four counter-rotating cells are captured for a ¼ 0 mm and three counter-rotating cells are captured for a ¼ 4, 8, 14 mm as shown in Fig. 6-c. As the liquid layer increases, the rolls at the center disappear to let those at the bottom corners to increase. So at t ¼ 600 s, when the liquid fills nearly the half of the domain, two counter-rotating rolls with the same size and intensity are obtained for 0 mm  a  14 mm. What can be noted here as a new finding, is the fact that, for various amplitudes, the warm fluid rises through the vertical walls and falls from the center of the cavity. For this reason, the interface shows a trough in the center. In some other numerical works, as in Ref. [13], the same number of rolls is obtained, but the fluid circulation is opposite to the present and the liquidesolid interface has a crest in the center. As the melting progresses, the liquid domain is more extended, the symmetry of the flow is broken and one of the two rolls becomes arbitrarily dominating. The dominating cell increases in size in contrast to the other as shown at t ¼ 1000 s in Fig. 6-e. The liquidesolid interface is more advanced in the proximity of the dominating cell. The local Nusselt number at the hot wall is a good indicator of how convection affects overall conduction through the cavity. Dimensionless local Nusselt number at the hot side can be defined as below:

NuX ¼

vT vY X TH  Tm

roll pairings. The drops are more pronounced for the wavy surface. We note that the amplitude of the wavy surface impacts the distribution of the local heat flux along the wavy. We can also note that the magnitude of these peaks increases as melting progresses. From an engineering point of view, the rate of melting is one of the most important parameters of the problem. The time evolution of the total liquid fraction in the cavity (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. 9 displays the time evolution of the total fraction of the liquid in the cavity. It appears from this figure that the increase of the amplitude leads to an increase in the rate melting inside the cavity. At early times, heat transfer in the melt zone predominantly by conduction and the evolution of liquid fraction versus time does not seem to vary by varying the magnitude of the amplitude value of the wavy surface. However, as melting progresses, the convection takes over conduction and the rate of the melting increases with the elevation in the magnitude of the amplitude value of the wavy surface. This can be explained by the fact that the increase in the amplitude value of the wavy surface ameliorates the heat transfer at the bottom part of the cavity and leads to an increase in the melting rate in this region.

(12)

The local Nusselt number was calculated from the converged temperature field at each time step. The impact of the amplitude on the Nusselt number for several times is shown in Figs. 7 and 8. Sudden drops of the Nusselt plot appear at times corresponding to

5. Conclusion Melting process along a vertical wavy surface with uniform surface temperature has been studied numerically. It is found that the rate of the melting increases with the elevation in the magnitude of the amplitude value of the wavy surface. Since the use of the wavy surface allows the increase in the exchange surface between the wavy surface and the PCM, it seems interesting to study the melting process over a wavy surface under a constant heat flux condition. Along these lines, present investigations are now focusing on the optimization of the ratio of amplitudes to wavelengths of the surface wavy surface as well as determining how to integrate efficiently such a device into a real system. References

Fig. 8. Nusselt number at time 600 s for various amplitudes.

[1] Dincer I, Rosen MA. Thermal energy storage, systems and applications. New York: Wiley; 2002. [2] El Omari K, Kousksou T, Le Guer Y. Impact of shape of container on natural convection and melting inside enclosures used for passive cooling of electronic devices. Appl Therm Eng 2011;31:3022e35. [3] Wang YH, Yang YT. Three-dimensional transient cooling simulations of a portable electronic device using PCM (phase change materials) in multi-fin heat sink. Energy 2011;36:5214e24.

T. Kousksou et al. / Energy 64 (2014) 212e219 [4] Cheng WL, Mei BJ, Liu YN, Huang YH, Yuan XD. A novel household refrigerator with shape-stabilized PCM (phase change material) heat storage condensers: an experimental investigation. Energy 2011;36:5797e804. [5] Kousksou T, Bruel P, Cherreau G, Leoussoff V, El Rhafiki T. PCM storage for solar DHW: from an unfulfilled promise to a real benefit. Sol Energy 2011;85: 2033e40. [6] Kandasamy R, Wang X, Mujumdar AS. Transient cooling of electronics using phase change material (PCM)-based heat sinks. Appl Therm Eng 2008;28: 1047e57. [7] Nayak KC, Saha SK, Srinivasan K, Dutta P. A numerical model for heat sinks with phase change materials and thermal conductivity enhancers. Int J Heat Mass Transf 2006;49:1833e44. [8] Gong ZX, Devahastin S, Mujumdar AS. Enhanced heat transfer in free convection-dominated melting in a rectangular cavity with an isothermal vertical wall. Appl Therm Eng 1999;19:1237e51. [9] Khodadadi JM, Hosseinizadeh SF. Nanoparticle-enhanced phase change materials (NEPCM) with great potential for improved thermal energy storage. Int Commun Heat Mass Transf 2007;34:534e43. [10] Fteiti M, Ben Nasrallah S. Numerical study of interaction between the fluid structure and the moving interface during the melting from below in a rectangular closed enclosure. Comput Mech 2005;35:161e9. [11] Hernandez-Guerrero A, Aceves SM, Cabrera-Ruiz E, Romero-Mendez R. Effect of cell geometry on the freezing and melting processes inside a thermal energy storage cell. J Energy Resour Technol 2005;127:9. [12] Lacroix M. Numerical simulation of melting and resolidification of a phase change material around two cylindrical heat exchangers. Numer Heat Transf Part A 1993;24:143e60. [13] Gong ZX, Mujumdar AS. Flow and heat transfer convection-dominated melting in a rectangular cavity heat from below. Int J Heat Mass Transf 1998;41:2573e80. [14] Prud’Homme M, Hung Nguyen T, Wu YK. Simulation numérique de la fusion à l’intérieur d’un cylindre adiabatique chauffé par le bas. Int J Heat Mass Transf 1991;34:2275e86.

219

[15] Lacroix M. Etude numérique de la fusion dans une enceinte rectangulaire chauffée par le bas. Trans Can Soc Mech Eng 1992;16:165e83. [16] Kousksou T, Arid A, Majid J, Zeraouli Y. Numerical modeling of doublediffusive convection in ice slurry storage tank. Int J Refrig 2010;33:1550e8. [17] Yen YC, Galea F. Onset of convection in a layer of water formed by melting of ice from below. Phys Fluid 1968;11:1263e9. [18] Yen YC. Free convection heat transfer characteristics in a melt water layer. J Heat Transf 1980;103:550e6. [19] Seki N, Fukosako S, Sugawara M. A criteria on onset of free convection in a horizontal melted water layer with free surface. Trans ASME Ser C J Heat Transf 1977;99:92e8. [20] Hale NW, Viskanta R. Solid-liquid phase change heat transfer from above or below. Int J Heat Mass Transf 1980;23:283e92. [21] Gau C, Viskanta R, Ho CJ. Flow visualization during solid-liquid phase change heat transfer II. Melting in rectangular cavity. Int Commun Heat Mass Transf 1983;10:183e90. [22] Diaz LA, Viskanta R. Visualization of the solid-liquid interface morphology formed by natural convection during melting of a solid from below. Int Commun Heat Mass Transf 1984;11:35e43. [23] Tashtoush B, Al-Odat M. Magnetic field effect on heat and fluid flow over a wavy surface with a variable heat flux. J Magn Magn Mater 2004;268:357e63. [24] Wang C-C, Chen C-K. Forced convection in a wavy wall channel. Int J Heat Mass Transf 2002;45:2587e95. [25] Cheng C-Y. Natural convection heat and mass transfer near a vertical wavy surface with constant wall temperature and concentration in a porous medium. Int Commun Heat Mass Transf 2000;27:1143e54. [26] Kim S, Kim MC, Lee B. Numerical analysis of convection driven melting and solidification in a rectangular enclosure. J Ind Eng Chem 2002;8:185e90. [27] Voller V. Fast implicit finite-difference method for the analysis of phase change problems. Numer Heat Transf Part B Fundam 1990;17:155e69. [28] Gau C, Viskanta R. Melting and solidification of a pure metal on a vertical wall. J Heat Transf 1986;108:174e81.