Impact of shape of container on natural convection and melting inside enclosures used for passive cooling of electronic devices

Impact of shape of container on natural convection and melting inside enclosures used for passive cooling of electronic devices


The present paper numerically analyzes a passive cooling system using enclosures with different geometries filled with thermal conductivity-enhanced phase change material (PCM). A numerical code is developed using an unstructured finite-volume method and an enthalpy-porosity technique to solve for natural convection coupled to a solid-liquid phase change. Five geometries containing the same volume of PCM are compared while cooling the same surface. The unsteady evolution of the melting front and the velocity and temperature fields is detailed. Other indicators of cooling efficiency are monitored, including the maximum temperature reached at the cooled surface. The computational results show the high impact of varying geometry: a maximum temperature difference as high as 40 C is observed between two of the enclosures. The best efficiency is obtained for an enclosure shifted vertically relative to the cooled surface. Other findings and recommendations are made for the design of PCM-filled enclosures.

Laboratoire des Sciences de l’Ingénieur Appliquées à la Mécanique and et au Génie Electrique (SIAME) and Université de Pau et des Pays de l’Adour (UPPA) and Bat. D’Alembert, rue Jules Ferry, BP. 7511, 64075 Pau cedex - France 1

Keywords: Passive cooling, phase change material, melting, natural convection, enthalpy-porosity method, unstructured finite-volume method.

1 Introduction

The melting of phase change materials (PCM) coupled to natural convection in enclosures has been studied extensively. This situation is encountered in many technical applications, such as latent heat storage systems [1, 2] and thermal insulation for buildings [3, 4]. The melting of PCM is also used to control the temperature of the surface of electronic components that release instantaneous or periodic high density heat fluxes to moderate the need for classical cooling devices. There has been increasing interest in this type of passive cooling for electronic circuits, such as chipsets, laptop processors or graphics cards [5, 6]. These elements are continuously becoming smaller, and their released heat densities are increasing; therefore, they require efficient, economic, and silent cooling systems. A particular application for such systems is the cooling of electrical devices inside shelters located in difficult-to-access sites and subject to periodic temperature changes. In these situations, cooling solutions that do not use moving mechanical elements, such as fans, are more suitable to avoid the need for frequent maintenance.

Several authors have studied phase changes and natural-convection-dominated melting in enclosures. Beyond the development and application of numerical models and experimental techniques, these authors have proposed several methods and concepts to enhance the heat transfer and melting rate. Gong et al. [7] demonstrated the positive effect of inverting a square PCM container heated from the side during the melting process when thermal stratification in the melt slows the melting rate. For the same geometry, Khodadadi and Hosseinizadeh [8] numerically investigated the effect of adding nanoparticles to the PCM. They showed that when the enclosure is cooled laterally, the rate of freezing increases because of the presence of nanoparticles. Starting from a study of PCM melting in a square enclosure heated from below, Fteïti and Ben Nasrallah [9] examined the impact of the aspect ratio of the PCM-filled enclosure and found that flat enclosures exhibit faster melting but have a lower asymptotic limit of the total melt ratio (they imposed a temperature boundary condition equal to the melting temperature at the top side of the enclosure and a lower temperature at the bottom). Hernandez-Guerrero et al. [10] also studied the impact of the aspect ratio with a different numerical model and found similar qualitative results for the case of tall enclosures. It is worth noting that the two latter works compared differently shaped enclosures that contained varying quantities of PCM.

Several researchers have considered the cooling of a vertical surface that dissipates heat flux using a tall enclosure filled with PCM. Binet and Lacroix [11] performed numerical and experimental studies to analyze the impact of the positions of three heat sources, their sizes and the aspect ratio of the enclosure itself. A similar numerical study was conducted by Krishnan and Garimella [12] on the conjugate heat transfer through the enclosure walls. Pal and Joshi [13] studied a uniformly dissipating heat source using experiments and numerical modeling, and they established correlations for the temporal evolution of the parietal heat transfer and the melt fraction. Huang et al. [14] used a three dimensional (3D) numerical model to investigate the cooling of photovoltaic cells by a PCM enclosure equipped with metallic internal fins.

A large amount of research has been conducted to investigate the applicability of cooling by PCM for electronic and electric devices. For example, the work in [15] focused on the cooling of a mobile device using embedded PCM enclosures, [16, 17, 18, 6] focused on PCM-based heat sinks, while in [19] PCM was used to control the temperature of Li-ion batteries.

Unlike for thermal insulation, the efficiency of passive cooling with PCM is closely related to the rapidity of its melting: the temperature of the heat-dissipating surface does not rise rapidly because of the absorption of the latent heat of fusion.

The speed with which the PCM melts in a laterally heated enclosure is controlled by the conductivity of the solid PCM during the early stages and by the natural convection currents in the melt during the later stages. This natural convection flow depends on the thermo-physical properties of the PCM and the geometry of the enclosure. Thus, the shape and extent of the enclosure containing the PCM have an important effect on the kinetics of the melting front and, therefore, on the temperature of the cooled wall.

In this work, we are interested in the cooling of a vertical surface releasing a fixed heat flux by means of a PCM-filled enclosure that entirely covers the surface. The aim of this study is to propose and examine different geometric shapes and relative positions of the enclosure. Numerical simulation is a well-suited tool for this purpose, especially in the early stages of design. It allows different test situations to be compared before an experimental prototype is constructed. A numerical model based on a fixed grid method and applied to an unstructured finite-volume formulation has been developed for this purpose.

In the following sections, the studied system is described with the different considered geometries and the flow conditions are then detailed. Next, the physical model and the numerical method are presented, and the implemented code is validated. Afterwards, detailed results are presented concerning the evolution of the velocity and the temperature fields, the melting front and the rate of fusion. Finally, the geometric impact on the coupled effects of natural convection and phase change is analyzed. Conclusions are then drawn based on these results.

2 The studied system

In this study, we are interested in the case of a vertical wall releasing a uniform heat flux at a high rate but for a limited duration. We aim to maintain this surface (designated ) at a sufficiently low temperature to avoid damaging the electrical component behind it. To achieve this objective, we provide this surface with an enclosure filled with a PCM. The surface to be cooled constitutes one face, or a part of a face, of this enclosure while the other faces are exposed to the external air and are subject to natural convection. When the heat release from the surface begins, its temperature rises until it reaches the fusion temperature of the PCM, which causes the material to melt. Natural convection currents take place progressively in the melt, and the generated flow ensures that heat is transferred from the surface to the solid PCM and maintains the melting process. The strength of the fluid flow involved in this process is closely related to the shape of the formed melted area, the extent of which is limited by the boundaries of the enclosure. Therefore, we expect that the efficiency of the heat transfer and the rate at which PCM melts, strongly depend on the shape of the enclosure. This simple passive cooling device is appropriate for intermittent heat release and is studied here for a uniform heat flux rate of . For continuous heat release and higher heat flux densities, more sophisticated solutions could be considered such as heat pipes [20].

The “natural” or “intuitive” choice of the shape of the PCM container can be regarded as the rectangular geometry which has been extensively studied in the literature for its simplicity. It can also be considered that the circular shape has the advantage of the absence of corners that may retain non melted PCM. Furthermore, the “intuitive” choice of the relative position of the container to the heating surface would be a centered position. Thus, in this study we aim to verify the relevance of these different choices. Consequently, we proceeded to investigate the impact of different enclosure shapes on the natural convection flow coupled to the solid-liquid phase change. The five different geometries studied contained the same volume of PCM. Fig. 1 shows the computational meshes used to model the enclosures (the chosen meshes densities are described in section 4.4). In this figure, the location of the surface is indicated in gray. The area of is the same for all of the enclosures. These geometries (detailed in the next section) were not all considered at the beginning of this work but were designed progressively given the results obtained from the first tested shapes. However, the results obtained for the five enclosures will be presented together. Some of the enclosures are rectangular, whereas others have rounded corners. Barlett et al. [21] have studied natural convection in both rectangular and rounded-corner enclosures and have shown that the latter shape substantially improves heat transfer. Thus, we considered rounded shapes in the present study and investigated their potential role in situations involving a phase change.

A known drawback of the considered cooling system is the relatively low thermal conductivity of common PCMs (e.g., wax paraffin), which may cause a high temperature peak at the surface in the case of high heat flux early-on in the heat release, when the heat transfer occurs only by conduction. To overcome this limitation, several solutions have been investigated in the literature, such as utilizing fins [6], including a high thermal conductivity metallic foam [6, 22] or adding nanoparticles. Extensive reviews of all these solutions can be found in [23, 24]. In our study, we consider a solution that includes graphite nanoparticles with high thermal conductivity in the PCM. These types of materials have experienced extensive development and exhibit high potential for use in several applications. Many researchers [25, 26] have demonstrated that the thermal conductivity of these products can reach unity (SI units). As indicated above, in the present study, we focus on the impact of the enclosure shapes on the melting kinematics, natural convection strength and temperature extremes at the cooled surface. Therefore, we choose not to consider the details of the mixture of PCM and nanoparticles and the relationship between the percentage of nanoparticles and the thermophysical properties of the obtained PCM, such as viscosity, effective conductivity and specific heat capacity. Instead we consider a model PCM with constant thermophysical properties but enhanced thermal conductivity (see next section for details).

3 Flow conditions and PCM properties

As stated above, in this work, we study the melting of PCM inside five enclosures of different shapes that surround the same volume. We limit our investigation to two-dimensional configurations that can be modeled in a vertical plane. One of the sides of the first two enclosures on the left in Fig 1 (a and b) exactly matches the surface . The first enclosure is rounded (a half disc), and the second is rectangular. The three other enclosures (c, d and e) have sides that vertically exceeded the surface , and thus, their widths are smaller than those of enclosures (a) and (b). Enclosure (c) is oblong with rounded corners, whereas enclosure (d) is rectangular. The last enclosure (e) has the same geometry as enclosure (d), but it is translated upward such that its left and right sides exceed the surface only at the top. The surface has a fixed height of , whereas the dimensions of the five enclosures are adjusted to contain the same volume. These dimensions are given in Tab. 1 in terms of maximum width and height.

(a) 2.5000 5.0000
(b) 1.9635 5.0000
(c) 1.5742 6.5742
(d) 1.4933 6.5742
(e) 1.4933 6.5742
Table 1: Maximum widths and heights of the five enclosures (in cm).

The rate of heat flux at the surface is fixed to a constant value, , for all cases. The heat transfer to the external air at the other boundaries is considered to occur by natural convection and is modeled by a heat transfer coefficient , whereas the external air temperature is assumed to be equal to the initial temperature of the PCM, C.

We consider a model PCM in this study composed of a pure paraffin wax and graphite nanoparticles. To simplify the analysis in the study, the properties of this PCM are considered to be identical in the solid and liquid phases and non-temperature dependent (except in the body forces term). The values of the adopted properties are listed in Tab. 2. A maximum temperature difference characterizing the flow can be obtained by considering the final or maximum admissible temperature that reaches the surface , which corresponds to the temperature at which an electronic device is damaged (i.e., C): C. The Rayleigh number based on this temperature difference and on the height of the surface is , which corresponds to a laminar flow regime. The Stefan number corresponding to this limit situation is .

Dynamic viscosity ()
Density () 800
Thermal conductivity () 1
Specific heat () 2500
Latent heat ()
Thermal dilatation coefficient ()
Prandtl’s number ()
Melting temperature () C
Table 2: Thermophysical properties of the model enhanced-conductivity PCM .
Figure 1: The geometries of the five studied enclosures: a,b, c, d, and e (from left to right) and corresponding meshes (in gray: the position of the surface ).

4 Computational modeling

4.1 The governing equations

The unsteady equations governing the flow of incompressible non-isothermal fluids are solved over the entire computational domain. The three-dimensional equations of conservation of momentum, mass and energy (in terms of temperature) are considered valid for both solid and liquid phases, which are distinguished by a liquid volume fraction that takes values 0 and 1, for solid and liquid phases, respectively. These conservation equations are considered in their integral forms:


where is the velocity vector, the pressure and the temperature. , is the viscous stress tensor for a Newtonian fluid:


The integration occurs over a volume surrounded by a surface , which is oriented by an outward unit normal vector . The source term in the momentum conservation equation (Eq. (1)) contains two parts:


where is the coefficient of volumetric thermal expansion and the acceleration of gravity vector. The first part of this source term represents the buoyancy forces due the thermal dilatation. For sake of simplicity, the Boussinesq’s approximation is used in this comparative study that focuses on the impact of different geometries. The temperature was chosen as the mean temperature of the PCM liquid phase and was recalculated at each time step. Therefore, is more representative of the temperatures in the liquid phase throughout the whole process, especially when all of the PCM has melted. The second part of the source term (Eq. 5) is a penalization term that ensures zero velocity in the computational cells where the PCM is solid, i.e., where . The penalization coefficient , based on the Carman-Koseny relation for a porous medium, is written as a function of as follows:


where , and . This coefficient ensures a smooth transition between solid and liquid media. In the case of pure PCM, the phase change front is sharp, and the transition takes place over only one computational cell in the direction of the front displacement.

In Eq. (3), the last term of the right hand side (RHS) introduces the effect of the latent heat of phase change into the energy equation. This effect is accounted for by considering the variation in time of the liquid fraction .

This general three-dimensional (3D) model can be restricted to deal with 2D computations (e.g., in the () plan) without any change by considering a single layer of computational cells (in the direction) and by neutralizing the top and bottom faces (with respect to ). It is this approach that was used for all the computations presented in this study, thus, the considered meshes are composed of one layer of hexahedra with or without prisms.

4.2 The numerical method

The conservation equations (1, 2 and 3) are solved by implementing them in an in-house code called Tamaris. This code has a three-dimensional unstructured finite-volume framework that is applied to hybrid meshes. The variable values (, , and ) are stored in cell centers in a collocated arrangement. The cell shapes can vary (e.g., tetrahedral, hexahedral, prismatic or pyramidal).

Figure 2: A computational cell and one of its neighbors .

To describe the discretization method used in the code, we can write Eqs. (1) and (3) in the generic convection-diffusion form with respect to a conserved variable :


where is a diffusion coefficient and a source term. The spatial schemes used to approximate the diffusive and convective fluxes are both second-order accurate. The diffusion term is discretized by approximating the surface integrals with a sum over all cell faces (Fig. 2):


where is the area of face . For unstructured meshes, orthogonality is an exception, and it needs to be handled correctly. Thus, the normal gradient is decomposed into an implicit contribution that uses the values of at the centers of the two cells sharing the face (the first term on the RHS of Eq. (9)) and a non-orthogonality correction term treated explicitly by a deferred approach to preserve the second-order accuracy of the centered differencing. We use the over-relaxed decomposition suggested by [27] to enhance the convergence properties of the discretization of the diffusive term:


is the vector joining the centers of the two adjacent cells (see Fig. 2). The average gradient is interpolated from the gradients of these neighboring cells.

The gradients of the variables at the cell centers are computed by Gauss’ theorem:


where is the mean value of the variable interpolated using the values at the centers of two cells sharing face :


Once the gradient is calculated for all computational cells, the values are used to determine a new estimate of as follows:


These new values of are used to re-compute the gradients more accurately using Eq. (10) [28].

Convection terms are also transformed into a sum over faces by decomposing the surface :


where the face values require appropriate interpolation to be accurate and bounded. Thus, we use the non-linear high-resolution (HR) bounded scheme CUBISTA by Alves et al. [29] in the formulation of Ng et al. [30], where they expressed is a function of the upwind (UP) value of and its centered differencing (CD) value:


The coefficient is determined for each face based on the local shape of the flow solution using the normalized variable diagram (NVD) framework and observing the convection boundedness criterion (CBC) [31]. The first term of the RHS of Eq. (14) is accounted for implicitly, whereas the second term is treated explicitly with the deferred-correction practice.

The pressure-velocity coupling is ensured by the SIMPLE algorithm [32], whereas the mass fluxes at the cell faces are evaluated using the Rhie-Chow interpolation [33] to avoid pressure checkerboarding. The implicit three-time-step Gear’s scheme [34] of second-order accuracy is used to discretize the unsteady terms:


The superscript stands for the current time step and is the time step. The RHS of Eq. (7) is taken at time . By applying the former discretizations, the generic conservation Eq. (7) transforms into the algebraic form:


Within each iteration of the SIMPLE algorithm, after the resolution of the momentum equation and of the Poisson equation for the pressure correction [28, 32], the energy equation is solved, and the fluid fraction is updated. These last two steps are repeated until the variation of is sufficiently small; then, the next SIMPLE iteration starts unless the convergence for , and is achieved, in which case a new time step is considered. A diagram of this algorithm is given in Fig. 3. The resolution of the energy equation is integrated in the SIMPLE iteration to take into account the high level of its coupling with the momentum equation through the body forces term. Additionally, the liquid fraction correction is iterated with the energy equation resolving to tightly couple the temperature field with the liquid fraction field.

Figure 3: The liquid fraction updating procedure as included in the SIMPLE algorithm.

In this study, the liquid fraction is updated by the “new source” algorithm proposed by Voller [35], where the new value of at iteration and in cell is calculated as follows:


where is the melting temperature of the PCM, and is the coefficient of in the discretized Eq. (16) that corresponds to temperature. This update is followed by an overshoot/undershoot correction:


Following the “new source” algorithm [35], the energy equation is penalized in the computational cells belonging to the phase change front () to ensure that the temperature at these cells is equal to the melting temperature . This procedure is performed by adding a penalization source term, equal to , to the energy equations corresponding to these cells. This practice accelerates the convergence of the updating algorithm, and the number of iterations needed to reach convergence () may be lowered to 1 or 2 depending on the size of the mesh and time-steps. The value of is fixed in this work to . At each iteration, the discretization technique presented above leads to a linear system of algebraic equations in the form of Eq. (16) with a non-symmetric sparse matrix for each variable. These linear systems are solved using an ILU-preconditioned GMRES procedure implemented in the IML++ library [36]. In the scope of this work, all of the computational meshes were generated using the open-source software Gmsh [37].

4.3 The code validation

The present code has been successfully validated in some situations involving flow and heat transfer, as in [38, 39]. An additional code validation concerned with pure natural convection is given in A. We focus here on validating the code when applied to the case of melting of a pure PCM coupled to natural convection in the melt. The chosen test case is the 2D numerical benchmark presented by Hannoun et al. [40], which involves melting tin in a square enclosure subject on one side to a temperature higher than the melting temperature. The authors have presented extensive results obtained using a second-order accurate finite-volume method in structured meshes as fine as 600 600.

Figure 4: Comparison of the obtained melting front positions at several instants (left) and of the evolution of the total liquid fraction in the enclosure (right) with the benchmark results of Hannoun et al. [40].

We carry out a numerical simulation in the same conditions as those described in [40] with the same material properties using a uniform grid 200 200 in size. In Fig. 4, we plot the melting front at different times. These lines correspond to the isolines of . The melting fronts at times and present a wave-like shape that is due to the presence of three and two fluid recirculations (or rolls) respectively. The top two rolls have merged into a single one at about . We also present the total liquid fraction in the entire square enclosure, which is calculated as follows:


where is the volume of a computational cell , and the summation is over all cells in the computational domain.

All these results are compared to those in [40] and exhibit satisfactory correspondence.

4.4 Mesh size-dependence study

To choose the mesh size with the best compromise between accuracy and computational cost, we conduct a mesh size-dependence study to determine the number of computational cells of the mesh needed to achieve satisfactory accuracy. The rectangular enclosure (b) undergoing an unsteady melting process in the same conditions as those described in section 3 is chosen as an example to explain the mesh selection process. Three meshes (, and ) of increasing size that contain 6000, 11660 and 24000 cells, respectively, are considered. The results obtained with the three meshes are compared based on the time evolution of the mean values and on the instantaneous spatial fields (i.e., melting front position, temperature and liquid fraction).

Figure 5: Comparison of results obtained by three meshes at : (a) melting front positions and (b) temperature along the cutline in the liquid phase.

Fig. 5(a) gives the form of the melting front obtained by the three meshes at time . These fronts are obtained by plotting the isovalue line of the liquid fraction field. The fronts obtained by the three meshes are almost superimposed on each other, which indicates that the mesh size has little influence on the front position. More significant differences are observed in the temperature field, as seen in Fig. 5(b), which shows the variation of the temperature along a horizontal axis passing through the center of the enclosure. For clarity, the plot is limited to the liquid region, where the temperature is not constant. However, the results of the three meshes are comparable. To achieve a quantitative comparison, we consider 50 points along this cutline and calculate the local error in the temperature at each point as follows: , where stands for meshes and ; the results are compared to those of the finest mesh . Thus, the mean value of the error is , and the maximum error is . They are both reported in Tab. 3.

2.2% 8.1%
1.2% 6.1%
Table 3: Relative error of the local values of along the central horizontal cutline at using the results given by the finest mesh as the baseline.

In addition to the spatial results, we compare the time evolution of the total liquid fraction in the enclosure (Eq. (19)) and the mean temperature of the surface , which is calculated as follows:


where is the area of a mesh face located on the surface . Relative errors of meshes and compared to are calculated at each time step for and over . The time average and maximum values of these errors are reported in Tab. 4

0.11% 0.75% 0.68% 4.34%
0.06% 0.40% 0.09% 1.43%
Table 4: Mean and maximum values of the relative error of and for a duration of . .

It is evident from Tabs. 3 and 4 that the results given by mesh are sufficiently close to those given by mesh , which is two times finer. Thus, the size of mesh was adopted to perform the computations presented in this work, and all the meshes generated for the five enclosures have approximately 11,000 cells.

5 Results

The computations for the five enclosures were conducted assuming the same conditions and using identical numerical parameters. The time step size was fixed to a relatively small value , and at each time step, the convergence criteria of the SIMPLE algorithm were fixed at for the velocity residuals and for the temperature residual. The latter criterion of convergence was fixed to a lower value in order to achieve a higher level of energy conservation in its both sensible and latent forms as they are of different orders of magnitude, furthermore, the released energy is tightly related to the amount of melted PCM. The residuals are calculated from Eq. (16) as follows:


where is the number of cells in the computational domain. Within each SIMPLE iteration, the liquid fraction update algorithm is stopped when the difference between two successive values is less than . At the first SIMPLE iteration of a time step, a few iterations are needed for to converge (less than 3), and in the remaining SIMPLE iterations, only one iteration is needed. For numerical stability reasons, the imposed heat flux rate at the surface is progressively raised during the first 5 seconds to reach its constant value . This practice enables easier convergence for the early time steps.

5.1 Velocity field and melting front evolution

We show the flow patterns represented by velocity vectors in the five studied enclosures in Fig. 6. This figure also gives the position of the melting front (or the solid phase zones) at five different instants during the process. At , a thin vertical layer of PCM has melted, and it has roughly the same size and the same form in all the enclosures. The flow field has a regular form with ascendant and descendant parallel fluid currents. Up to this time, the five enclosures have the same behavior, and the geometry does not play a role. At , more PCM has melted, especially near the upper part of the surface . The behavior of the enclosures then differs visibly: the enclosures (a) and (b) that do not exceed the surface have a melting front that extends more horizontally due to the presence of the upper boundary. Thus, the shapes of melted zones are comparable between enclosures (a) and (b) and among (c), (d) and (e).

At , the melted PCM zone in the enclosure (e) extends upward with a round shape, and as a consequence, more solid PCM is exposed to the flow current of the liquid for this enclosure than the others (the liquid-solid interface is the most extended for enclosure (e)). From this instant on, the melting fronts in the five enclosures present an inclined curve because the melting is more advanced in the upper regions due to the thermal stratification (see Fig. 8). The flow patterns in the enclosure (e) show the formation of recirculation above the surface . This recirculation, visible at (Fig. 6), is in fact an unsteady oscillatory phenomenon with a period less than 2 s, which is why it is not visible in the figure when . These oscillations are induced by the formation of a large liquid space above the surface . At the upper extremity of , the hot flow stream rising along the surface oscillates between two orientations: a vertical orientation, as observed at , and a roughly horizontal orientation, as observed at or , which results in a secondary recirculation in the upper-right corner. This unsteady behavior starts as early as for enclosure (e) but much later () for enclosure (d) because the upper liquid space is narrower. The flows in the three other enclosures do not show any oscillatory characteristics. At , the melting in all the enclosures is almost complete. We can visually observe that enclosures (a) and (e) contain the least amount of solid PCM. This observation is confirmed quantitatively in Fig. 11.

To give more quantitative details about the velocity field in the enclosures, we plot the vertical component of the velocity vector along the horizontal axis passing through the center of the surface in Fig. 7. Initially (up to ), the velocities are similar, so they are not shown. At , we observe substantial differences in the maximum value of the velocity near the boundary , where enclosure (e) exhibits the highest value (), and enclosure (a) exhibits the lowest value (). The positions of the positive maximum are nearly the same, whereas in the fluid zone near the melting front, the downward negative velocities have different profiles due to the differences in the position of the front. As the process continues ( and ), the velocity level decreases progressively, and a large zone with almost zero vertical velocity appears in the center of the enclosures. The maximum value in enclosure (e) at is . This velocity decrease is due to the reduction of the quantity of solid PCM, which is at the low temperature , in the enclosures compared to the initial amounts. As a result, a weaker temperature gradient exists in the liquid which results in weaker buoyancy forces.

Figure 6: Velocity vectors and melting fronts for the five enclosures after 50 and 100 s.
Figure 6: (Cont.) Velocity vectors and melting fronts for the five enclosures after 200 and 300 s.
Figure 6: (Cont.) Velocity vectors and melting fronts for the five enclosures after 500 s.
Figure 7: Vertical velocity components along the horizontal axis passing through the center of the surface for the five enclosures after 200 s (a), 300 s (b) and 500 s (c).

5.2 Temperature field in the enclosures

Figure 8: Temperature fields in the five enclosures after 50, 100, 200, 300, 400 and 500 s.

Fig. 8 shows the temperature field at six different times during the melting process. As early as , we observe the establishment of thermal stratification in the melt, which is slightly perturbed in the upper region by the upward fluid currents and near the lateral boundaries by the convective heat exchange with the external media. Higher temperatures were observed in the liquid of enclosure (a), and lower temperatures were observed in enclosure (e). Quantitative changes in temperature are given in the next section. As for the melting front position, we observe similarities in the temperature fields between enclosures (a) and (b) and between enclosures (c) and (d), whereas enclosure (e) exhibits a distinct temperature pattern.

5.3 Evolution of the parietal temperature

The evolution of , the mean temperature of the surface (Eq. (20)), gives an overview of the impact of the flow and melting kinetics on the global temperature level (Fig. 9). During the first minute, the mean temperature evolves identically and linearly for all the enclosures. This first phase is governed essentially by heat conduction; thus, it gives the same results for all five enclosures. After this phase of steep augmentation, continues to increase in enclosure (a) but with a weaker slope, but it decreases in all the other enclosures. The strongest decrease is observed for enclosure (e), and it lasts for approximately 250 s before increasing again. For enclosures (b), (c) and (d), this period is approximately 170 s.

Another important temperature to monitor is the maximum temperature of the surface because this temperature should be lower than the damage temperature of the device to be cooled. Fig. 10 plots the maximum temperature for the five enclosures. This maximum is always located at the top of , as in Fig. 8. Its evolution shows differences between the enclosures as early as . For enclosure (a), exhibits two linear phases with positive slopes with inflexions at about and . The evolution of for enclosure (b) has a horizontal plateau at , whereas for the other three enclosures, this plateau is lower and is slightly decreasing. The plateau for enclosure (e) ( long) is the longest.

The plateau in the and curves emerges because of the establishment of natural convection currents and the presence of the front of solid PCM. The solid PCM absorbs the thermal energy transported by the flow from the hot surface in the form of latent heat, and its melting releases cold liquid at , which increases the temperature gradient. The circular shape of enclosure (a) induces a larger distance between the solid PCM and the surface . This distance is less important for enclosure (b), and its rectangular shape induces a plateau in the temperature profile and a lower value of . When the enclosures are flatter, as for (c) and (d), the results are better. The shape of the enclosure has another effect in terms of the extent of the melting front, which is more extended for the enclosures with the lowest aspect ratio. The more geometrically extended this melting front is in the enclosure, the lower and longer the plateau. By comparing the curves of enclosures (c) and (d) with (e), we can conclude that an efficient way to increase the extent of the front is to locate more PCM at the top of the enclosure where it is in contact with oscillatory natural convection streams. The higher temperature liquid is located in this zone because of thermal stratification, which increases the melting rate.

Figure 9: Evolution of the mean temperature of the surface .
Figure 10: Evolution of the maximum temperature of the surface .

5.4 Evolution of the global PCM melting

The total liquid fraction defined by Eq. (19) can be monitored to follow the evolution of the melting in the five enclosures, as shown in Fig. 11. We observe that the curve of is linear and identical in all the enclosures until . After this time, the rate of melting decreases in all the enclosures except for enclosure (e), for which the curve inflexion occurs as late as . This enclosure has the fastest rate of melting, and the solid PCM completely disappears at . Enclosure (a) is also completely melted at this time, and even though it shows the highest maximum surface temperature (Fig. 10), its mean surface temperature is lower than those of enclosures (b), (c) and (d) after (Fig. 9). Fig. 6 can explain these differences. The bottom left corner of the rectangular enclosure (b) retains the PCM solid and limits its melting. In enclosures (c) and (d), more solid PCM is trapped in the lower part of the enclosure located below the surface . This zone is less influenced by the hot liquid streams recirculating inside the enclosure, and thus, it is more difficult to melt. As a consequence, enclosures (c) and (d) are the last ones to completely melt, approximately after enclosures (a) and (e) are completely melted.

Figure 11: Evolution of the total liquid fraction in the five enclosures.

5.5 Parietal heat transfer

To quantify the heat transfer at the cooled surface, we define a global Nusselt number as follows:


Fig. 12 gives the evolution of the Nusselt number for the five enclosures during the melting process. After a fast transition period, where the values are important (), and before they drop to a value of , we observe an increase toward a peak value followed by a progressive decrease toward an asymptotic value of . The characteristic “bumps” in the curves obtained for the different enclosures have different sizes and correspond to the coupled impact of the natural convection and latent heat of the melting of the PCM on the temperature of the cooled surface . In agreement with the results of velocity (Fig. 7) and parietal temperature (Fig. 9), the most important value is obtained for enclosure (e), and the lowest is obtained for enclosure (a) up to . After , enclosure (a) performs better than enclosures (b), (c) and (d) when the rate of fusion in these enclosures slows down. However, the most interesting behaviors expected from this mode of passive cooling are the fast reaction and limited maximum temperature. With regard to these two features, enclosure (a) is the worst choice and enclosure (e) is the best choice. We can also notice that for all the efficiency indicators (, , and ), the oblong enclosure (c) performs slightly better than the rectangular enclosure (d).

Figure 12: Evolution of the parietal Nusselt number for the five enclosures.

6 Conclusions

In this study, we numerically modeled the natural convection-dominated melting of a PCM inside an enclosure used to control the temperature of a surface releasing a heat flux. We investigated the impact of the shape of this enclosure and the relative position of the cooled surface on the flow and heat transfer. We developed a numerical model using a fixed grid enthalpy-porosity technique coupled with a three-dimensional flow solver based on an unstructured finite-volume method involving second-order accurate spatial and temporal numerical schemes.

Five geometries containing the same quantity of PCM were compared. The geometries were rounded or rectangular, thick or thin, centered relative to the cooled surface or shifted upward vertically. The unsteady behaviors of the five enclosures were analyzed by examining the evolution of the liquid-solid interface and considering the forms and strengths of the natural convection flow. The fluid velocity increases for thinner enclosures that have a space above the cooled surface; thus, in the enclosure shifted upward (enclosure (e)) the maximum velocity is achieved. In spite of the presence of this flow, strong thermal stratification exists in the melt, which causes a hot spot located at the upper extremity of the surface to appear.

Several efficiency indicators were monitored, including the mean and maximum temperatures, the total liquid fraction and the parietal Nusselt number. All of these indicators demonstrate the importance of the effect of the geometry on cooling efficiency. If we focus on the maximum temperature indicator, we see that at as early as , the difference between the best and worst geometry choices ((e) and (a), respectively) is equal to C, and this difference reaches C when , which is of great importance for the thermal protection of the heat-dissipating devices. We can summarize the most important findings of this work as follows:

  • In thin enclosures, liquid PCM zones expand upward, while for wider enclosures, this expansion tends to be horizontal. Consequently, for the former enclosures, the solid cold PCM is in contact with ascending hot liquid streams, which contributes to the lowering of the global temperature.

  • The last zones to melt in the enclosures are located at the bottom and on the side opposite to the cooled surface. Thus, PCM enclosures should be designed without corners and without zones located below this surface.

  • A portion of the PCM should be placed above the cooled surface.

  • The use of rounded corners has a slight positive effect on the efficiency (e.g., and ).

Based on these findings, more efficient geometries can be designed; however, in the future, the impact of the geometry choice on the inverse process, i.e., solidification, should be taken into account and analyzed. Furthermore, the relations between the fraction of nanoparticles introduced and the thermophysical properties of the PCM such as the viscosity of the liquid phase, should be modeled and their impacts on hydrodynamics and heat transfer analyzed. These aspects constitute some of the goals of a future work.

Appendix A Test case of natural convection between a cylinder and a square duct

Here we use the benchmark test case proposed by Demirdzic et al. [41] to validate the implementation of natural convection in the code. A cylinder is placed inside a square enclosure of dimensions with . The cylinder has a diameter of and its center is displaced vertically by . The cylinder is heated at , while the vertical walls of the duct are cooled at and the horizontal walls are assumed to be adiabatic. The Rayleigh number based on is , and the Prandtl number is . The benchmark solution of Demirdzic et al. [41] was obtained using a structured mesh of . Because of the symmetry of the problem, only one-half of the domain needs to be modeled. For our computation we used a mesh composed of cells, with quadrilateral cells around the walls and triangular cells elsewhere as showed in Fig. 14. In this figure are presented the calculated streamlines and temperature contours that compare well with those of [41]. The local Nusselt number along the vertical wall of the duct is presented in Fig. 14. A very good agreement is found between our results and those of the benchmark. Furthermore, the calculated mean Nusselt number at this wall is while the benchmark value is which gives a relative difference of less than .

Figure 13: Calculated results of the natural convection test case of Demirdzic et al. [41]. Right half: the used grid and temperature contours (values from to with a step of ). Left half: streamline pattern.
Figure 14: Distribution of the local Nusselt number along the vertical cold wall compared to the benchmark data of Demirdzic et al. [41].


  1. footnotetext: Corresponding author:
    Preprint published in Applied Thermal Engineering, Volume 31, Issue 14-15, October 2011, Pages 3022-3035.


  1. K. El Omari, J.-P. Dumas, Crystallization of supercooled spherical nodules in a flow, International Journal of Thermal Sciences 43 (12) (2004) 1171–1180.
  2. T. Kousksou, J.-P. Bédécarrats, J.-P. Dumas, A. Mimet, Dynamic modelling of the storage of an encapsulated ice tank, Applied Thermal Engineering 25 (10) (2005) 1534–1548.
  3. E. Alawadhi, Thermal analysis of a building brick containing phase change material, Energy and Buildings 40 (3) (2008) 351–357.
  4. M. Ahmad, A. Bontemps, H. Sallée, D. Quenard, Thermal testing and numerical simulation of a prototype cell using light wallboards coupling vacuum isolation panels and phase change material, Energy and Buildings 38 (6) (2006) 673–681.
  5. R. Kandasamy, X. Wang, A. S. Mujumdar, Transient cooling of electronics using phase change material (PCM)-based heat sinks, Applied Thermal Engineering 28 (8-9) (2008) 1047 – 1057.
  6. K. C. Nayak, S. K. Saha, K. Srinivasan, P. Dutta, A numerical model for heat sinks with phase change materials and thermal conductivity enhancers, International Journal of Heat and Mass Transfer 49 (11-12) (2006) 1833–1844.
  7. Z. Gong, S. Devahastin, A. Mujumdar, Enhanced heat transfer in free convection-dominated melting in a rectangular cavity with an isothermal vertical wall, Applied Thermal Engineering 19 (12) (1999) 1237–1251.
  8. J. Khodadadi, S. Hosseinizadeh, Nanoparticle-enhanced phase change materials (NEPCM) with great potential for improved thermal energy storage, International Communications in Heat and Mass Transfer 34 (5) (2007) 534–543.
  9. M. Fteïti, S. Ben Nasrallah, Numerical study of interaction between the fluid structure and the moving interface during the melting from below in a rectangular closed enclosure, Computational Mechanics 35 (3) (2005) 161–169.
  10. A. Hernández-Guerrero, S. Aceves, E. Cabrera-Ruiz, R. Romero-Mendez, Effect of cell geometry on the freezing and melting processes inside a thermal energy storage cell, Journal of Energy Resources Technology 127 (2) (2005) 95–102.
  11. B. Binet, M. Lacroix, Melting from heat sources flush mounted on a conducting vertical wall, International Journal of Numerical Methods for Heat and Fluid Flow 10 (3) (2000) 286–306.
  12. S. Krishnan, S. Garimella, Analysis of a phase change energy storage system for pulsed power dissipation, IEEE transactions on components and packaging technologies 27 (1) (2004) 191–199.
  13. D. Pal, Y. Joshi, Melting in a side heated tall enclosure by a uniformly dissipating heat source, International Journal of Heat and Mass Transfer 44 (2) (2001) 375–387.
  14. M. Huang, P. Eames, B. Norton, Comparison of a small-scale 3D PCM thermal control model with a validated 2D PCM thermal control model, Solar Energy Materials and Solar Cells 90 (13) (2006) 1961–1972.
  15. F. Tan, C. Tso, Cooling of mobile electronic devices using phase change materials, Applied thermal engineering 24 (2-3) (2004) 159–169.
  16. V. Shatikian, G. Ziskind, R. Letan, Numerical investigation of a PCM-based heat sink with internal fins: Constant heat flux, International Journal of Heat and Mass Transfer 51 (5-6) (2008) 1488–1493.
  17. R. Akhilesh, A. Narasimhan, C. Balaji, Method to improve geometry for heat transfer enhancement in PCM composite heat sinks, International Journal of Heat and Mass Transfer 48 (13) (2005) 2759–2770.
  18. X. Wang, C. Yap, A. Mujumdar, A parametric study of phase change material (PCM)-based heat sinks, International Journal of Thermal Sciences 47 (8) (2008) 1055–1068.
  19. R. Kizilel, R. Sabbah, J. Selman, S. Al-Hallaj, An alternative cooling system to enhance the safety of Li-ion battery packs, Journal of Power Sources 194 (2) (2009) 1105–1112.
  20. D. Reay, P. Kew, Heat Pipes – Theory, Design and Applications, 5th Edition, Butterworth-Heinemann, Oxford, 2007, Ch. 8: Cooling of Electronic Components, pp. 319–341.
  21. A. Barletta, E. Nobile, F. Pinto, E. Rossi di Schio, E. Zanchini, Natural convection in a 2D-cavity with vertical isothermal walls: Cross-validation of two numerical solutions, International Journal of Thermal Sciences 45 (9) (2006) 917–922.
  22. O. Mesalhy, K. Lafdi, A. Elgafy, K. Bowman, Numerical study for enhancing the thermal conductivity of phase change material (PCM) storage using high thermal conductivity porous matrix, Energy Conversion and Management 46 (6) (2005) 847–867.
  23. S. Jegadheeswaran, S. Pohekar, Performance enhancement in latent heat thermal storage system: A review, Renewable and Sustainable Energy Reviews 13 (9) (2009) 2225–2244.
  24. 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 (LHTESS), Renewable and Sustainable Energy Reviews 14 (2) (2010) 615–628.
  25. S. Kim, L. Drzal, High latent heat storage and high thermal conductive phase change materials using exfoliated graphite nanoplatelets, Solar Energy Materials and Solar Cells 93 (1) (2009) 136–142.
  26. A. Sari, A. Karaipekli, Thermal conductivity and latent heat thermal energy storage characteristics of paraffin/expanded graphite composite as phase change material, Applied Thermal Engineering 27 (8-9) (2007) 1271–1277.
  27. H. Jasak, Error analysis and estimation for the finite volume method with applications to fluid flows, Ph.D. thesis, University of London (1996).
  28. J. H. Ferziger, M. Perić, Computational methods for fluid dynamics, Springer New York, 2002.
  29. M. A. Alves, P. J. Oliveira, F. T. Pinho, A convergent and universally bounded interpolation scheme for the treatment of advection, International Journal for Numerical Methods in Fluids 41 (2003) 47–75.
  30. K. C. Ng, M. Z. Yusoff, E. Y. K. Ng, Higher-order bounded differencing schemes for compressible and incompressible flows, International Journal for Numerical Methods in Fluids 53 (1) (2007) 57–80.
  31. P. H. Gaskell, A. K. C. Lau, Curvature-compensated convective transport: SMART, a new boundedness-preserving transport algorithm, International Journal for Numerical Methods in Fluids 8 (6) (1988) 617–641.
  32. S. V. Patankar, Numerical heat transfer and fluid flow, Taylor & Francis, 1980.
  33. C. Rhie, W. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA Journal 21 (11) (1983) 1525–1532.
  34. C. W. Gear, Numerical initial value problems in ordinary differential equations, Prentice-Hall Series in Automatic Control, Prentice-Hall, Englewood Cliffs, N.J., 1971.
  35. V. Voller, Fast implicit finite-difference method for the analysis of phase change problems, Numerical Heat Transfer Part B-Fundamentals 17 (1990) 155–169.
  36. J. Dongarra, A. Lumsdaine, R. Pozo, K. Remington, A sparse matrix library in C++ for high performance architectures, in: Proceedings of the Second Object Oriented Numerics Conference, 1994, pp. 214–218.
  37. C. Geuzaine, J. F. Remacle, Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (11) (2009) 1309–1331.
  38. K. El Omari, Y. Le Guer, A numerical study of thermal chaotic mixing in a two rod rotating mixer, Computational Thermal Science 1 (2009) 55–73.
  39. K. El Omari, Y. Le Guer, Thermal chaotic mixing of power law fluids in a mixer with alternately-rotating walls., Journal of Non-Newtonian Fluid Mechanics 165 (2010) 641–651.
  40. N. Hannoun, V. Alexiades, T. Mai, A reference solution for phase change with convection, International Journal for Numerical Methods in Fluids 48 (11) (2005) 1283–1308.
  41. I. Demirdžić, Ž. Lilek, M. Perić, Fluid flow and heat transfer test problems for non-orthogonal grids: Bench-mark solutions, International Journal for Numerical Methods in Fluids 15 (3) (1992) 329–354.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description