NUMERICAL SIMULATION OF NATURAL CONVECTION MELTING IN 2D AND 3D ENCLOSURES

Natural convection melting in 2D and 3D enclosures with a local heater is studied numerically. The present research is related to a development of effective cooling system for the electronic devices using the phase change material that is essentially important nowadays. The domain of interest includes vertical cold walls, adiabatic horizontal walls and a discrete heater of constant high temperature located on the bottom adiabatic wall. The cavity is filled with a phase change material (PCM) in solid state at the beginning of the process. During the heating from the heat source PCM is melting. Numerical solution of the present problem has been conducted using the dimensionless transformed variables such as stream function and vorticity for 2D cavity and vector potential functions and vorticity vector for 3D cavity with appropriate initial and boundary conditions. The developed numerical technique has been verified comprehensively. Obtained results have shown a potential of the used methods for 2D and 3D problems. It has been found that, melting process is more intensive in 3D case and the heat transfer rate at the heater is greater for 2D in comparison with 3D data.


INTRODUCTION
Nowadays, one of effective cooling methods for modern electronic devices is a usage of phase change materials (PCM). These materials have a moderate melting temperature, high latent melting energy and high heat capacity. Heat generation or absorption inside these materials due to the phase transformation with constant temperature is widely used not only for electronic cooling but also in solar energy storages, building temperature control systems and in other systems [1][2][3][4][5][6][7][8][9][10][11][12][13][14].
Thus, Fan et al. [2] experimentally have analyzed the effects of melting temperature and the presence of internal fins on behavior of PCM for thermal management of electronics. It has been found that application of a PCM with a higher melting temperature can extend a longer time of protection of the devices from over-heating. Gharbi et al. [3] have studied experimentally the heat transfer behavior of PCM in a cavity with a local heater having a radiator system. The authors have revealed that an addition of PCM extends twice the critical time of the heat sink and the radiator system consisting of solid fins essentially changes the thermal control. Numerical simulation of PCM melting within an enclosure with internal horizontal rectangular fins has been performed by Sharifi et al. [4]. It has been found that utilization of horizontal rectangular fins initially promotes rapid melting, followed by slower melting once PCM in the inter-fin regions has been liquefied. Gharbi et al. [5] have studied experimentally the thermal performance of PCM inside vertical cavity with three discrete heaters on vertical wall. The flow structures and heat transfer patterns have been analyzed for different heat flux cycles. Taghilou and Talati have examined numerically [6] a melting process within a finned rectangular cavity filled with PCM using the lattice Boltzmann method. The obtained results have shown that using of heat-conducting solid fins can lead to more essential melting of PCM. Prieto et al. [7] have investigated numerically the thermal behavior of PCM heat exchangers for heating systems. The RT60 paraffin has been used for analysis. The authors used a full CFD model in Ansys Fluent 14.5 and found that melting and solidification in the vertical arrangement are slower than in the horizontal one. Madruga and Mendoza [8] have studied numerically the melting of n-octadecane under the thermocapillary effect at a free surface inside the circular sections. It has been shown that thermocapillary effect on the free surface has a powerful influence on the heat transfer rate, decreasing a total melting time. Wang et al. [9] have experimentally analyzed the thermal behavior of the vertical cylindrical battery with composite paraffin and fin structure. Ren and Chan [10] have studied numerically the PCM melting process in a cavity with internal fins under the lateral heating using the lattice Boltzmann method. They have found that using internal fins could enhance the heat transfer process in the PCM cavity, namely, the PCM melting speed can be increased with the length of the internal fins. Moreover, the melting speed is enhanced when the PCM cavity is inclined towards the counter clockwise direction. Wang et al. [11] have investigated numerically an opportunity to optimize the latent thermal energy storage system using the PCM and internal fins. Obtained results have shown that fin-length, fin-ratio and angle between neighbor fins have affected the heat transfer enhancement inside the PCM storage system. Tabassum et al. [12] have examined the melting process of PCM in the annular gap using the finite volume method. Mathematical model has been formulated on the basis of the enthalpy-porosity approach using the primitive variables. Comprehensive analysis has revealed that the eccentric position of the inner tube in the inverse triangular annulus is an excellent configuration for the storage system among the studied heat exchangers. Yang et al. [15] have explored the thermal performance of the low melting point metal (LMPM) based internal finned heat sink using the experimental setup and commercial Fluent package. It has been found that LMPMs can be a very good candidate as a PCM for thermal energy storage and thermal management.
Abovementioned analysis of several works on PCM melting shows that in the case of numerical study the primitive variables are the main variables for investigation and 2D study is a wide-spread approach, while in the case of 3D simulation the commercial packages or open access codes are used.
The objective of the present study is to numerically analyze unsteady melting driven by natural convection in 2D and 3D cavities having local heaters using finite difference method and transformed variables such as stream function and vorticity for 2D cavity and vector potential functions and vorticity vector for 3D cavity. It should be noted that, it the first time for usage of vector potential functions and vorticity vector for numerical analysis of natural convection melting within cubical cavity with a local heater. 2D analysis of the natural convection melting in a square cavity has been performed earlier by Bondareva and Sheremet [16,17] using also the uniform magnetic field. 3D analysis has been conducted also in the case of uniform magnetic field by Bondareva and Sheremet [18,19]. The present work is a more detailed comparison of the obtained data for 2D and 3D cases with highlighting an opportunity to use non-primitive variables for simulation of PCM melting process in 3D domains.

GOVERNING EQUATIONS AND NUMERICAL TECHNIQUE
We consider 2D and 3D cavities ( Figure 1) with a heat source of constant temperature mounted on the bottom wall. At the beginning of the process the domain of interest is filled with PCM in solid state having fusion temperature. Two opposite vertical walls of the enclosure have a low constant temperature, while other walls are adiabatic. Melt is considered to be incompressible, radiatively non-participating and the Boussinesq approximation is valid. The melt motion and heat transfer in the cavity are assumed to be two-dimensional for the square cavity ( Figure 1a) and three-dimensional for the cubical enclosure ( Figure 1b). It should be noted that 2D cavity can be considered as a middle cross-section of 3D cavity. The melt flow and heat transfer are described using Navier-Stokes equations, including the energy equation for the melt formulated on the basis of enthalpy approach. The governing partial differential equations can be written as follows in tensor form for the liquid state of material (inside the melt zone) [16][17][18][19]: It should be noted that here we used the enthalpy formulation for the energy equation.
For the solid state of PCM we utilized the heat conduction equation in following form To combine the energy equations for melt and solid material and to delete the miscoordinations for enthalpy function at solid-liquid interface the following smoothing function was utilized: It is worth noting that solution of the fundamental problems is more useful to perform in dimensionless variables. Therefore, we introduce the following dimensionless variables in Cartesian coordinates: In the case of 3D cavity (Figure 1b) we define the dimensionless vector potential functions in the following form [18][19][20]: and vorticity vector in the form [18][19][20]:
In the case of 3D cavity we have the following initial and boundary conditions [18,19]

ZZ  
Governing equations with corresponding initial and boundary conditions for 2D cavity Eqs. (6)- (8) and (16) and for 3D cavity Eqs. (9)- (15) and (17) have been solved by finite difference method of the second order accuracy [16][17][18][19][20][21][22]. The diffusive terms have been approximated by central differences. The convective terms have been discretized by the second order Samarsky monotonic difference scheme. The parabolic equations have been solved on the basis of Samarsky locally one-dimensional scheme. The obtained systems of algebraic equations have been solved by Thomas algorithm. The partial differential equations for the stream function or vector potential functions have been discretized by means of the five/seven-points difference scheme on the basis of central differences for the partial derivatives. The obtained linear discretized equations have been solved by the successive over relaxation method. It should be noted that some difficulties have been met in definition of vorticity vector components at the phase line. This definition has been performed using the additional nodes and more complicated approximation of the Poisson equations for vector potential functions.
The comprehensive verification has been carried out including grid independence test and some natural convection problems. Detailed description of this validation part can be found in [16][17][18][19][20].
In the case of 3D grid independency test, we analyzed the position of solid-liquid interface at the middle cross-section Y = 0.5 for Ra = 3.210 6 with three time moments such as  = 520.2,  = 867,  = 1387.2 and using two different grids of 60×60×60 points (solid lines) and 100×100×100 points (dashed lines). Figure 2 shows a location of solid-liquid interface for the considered grids.
Taking into account the non-significant changes in the obtained locations, a uniform grid of 60×60×60 points has been used for analysis. Figure 3 shows a good comparison for the position of phase line with experimental data of Gau and Viskanta [23] in the case of natural convection melting of pure gallium in a rectangular cavity. The dimension of the cavity used in the experiments is 6.35 cm and 8.89 cm. The left hot wall and the right cold wall are maintained at temperatures 38.35C and 28.3C, respectively, and the cavity was initially filled by solid PCM at temperature 28.3C. The horizontal walls are insulated. The physical properties of pure gallium were chosen at 32C, a temperature that is representative of the temperature range of the experiment. One can find in Figure 3 a location of phase line at different time moments. At the beginning of the melting process (2 min) the main heat transfer mechanism is a heat conduction, therefore the phase line is parallel to the vertical isothermal wall. Further, convection develops with a formation of convective cell in the upper part of the cavity, where phase line bends. Comparing the location of the calculated phase line with experimental data it is possible to conclude that main differences are for initial time level (for some points at 6 min the relative discrepancy is about 20%). While for more intensive melting this discrepancy is low (at 17 min the relative discrepancy is less than 10%). It should be noted that the main reason for such differences is an anisotropic gallium used in the experiments.  Figure 4 demonstrates an evolution of streamlines and isotherms in 2D cavity and in the middle crosssection Y = 0.5 of 3D cavity for Ra = 510 4 . It should be noted that, small value of buoyancy force does not lead to differences between 2D and 3D formulations. One can find only some differences in a location of isotherm  = -0.2 near the upper wall. At the same time an increase in time leads to a diminution of this difference.
Regardless of the dimensionless time value, two convective cells are formed inside the melt zone, illustrating an ascending flow in central part and descending flows near the phase transition line. An increase in time leads to a growth of the melting zone and a development of convective thermal plume over the heater. More essential melting occurs in vertical direction in comparison with horizontal direction. Figure 5 demonstrates distributions of streamlines and isotherms for 2D and 3D cases at Ra = 510 7 ,  = 219. It is an initial time level for high Rayleigh number value. One can find a presence of some differences in melting flow structure and as a result in temperature profiles. Main differences are in convective cells over the heater due to different heat fluxes for 2D and 3D surfaces. More detailed comparison between 2D and 3D results can be performed on the basis of vertical velocity and temperature profiles presented in Figure 6 for Ra = 410 5 and different values of the dimensionless time. First of all it is possible to analyze an evolution of the melting process in the following form. Heating from the local heat source reflects a formation of ascending convective flow in central part and two descending flows near the phase line. Motion velocity in these flows increases with time and melting distributes also in horizontal direction taking into account the zero velocity values. At initial time, melt circulation in 2D case is more essential in comparison with 3D case, but this time these differences vanish. A displacement of the central ascending flow in 3D case from the symmetry line can be explained by a formation of the transient 3D regime with unstable thermal plume.
Temperature profiles (see Figure 6b) also reflect more essential heating in 2D case at initial time, while further on temperature in cubical cavity is greater then temperature in a square cavity.  Evolution of the local Nusselt number at Ra = 410 5 for 2D and 3D cases is shown in Figure 7. Behavior of Nu along the heater surface is similar for 2D and 3D cases. Along the vertical surfaces the local Nusselt number decreases from the bottom part till upper surface due to a formation of expanding thermal boundary layer near these surfaces. Along the upper heater wall one can find a maximum Nu values near the edges of the upper surface and a minimum value in the central part. The main reason for such behavior is a formation of a thermal plume in the central part where the temperature gradient decreases. The local Nusselt number increases with time due to more essential heating. 2D case illustrates great values of Nu in comparison with 3D case.

CONCLUDING REMARKS
Natural convection melting of paraffin in a cubical cavity and square cavity with a local heater along the bottom wall has been analyzed numerically. Governing equations formulated in dimensionless "vector potential functions -vorticity vector" for 3D case and "stream function -vorticity" for 2D case have been solved by finite difference method of the second order accuracy. Comparison between 2D and 3D data illustrates an efficiency of the developed mathematical model and numerical technique. The presented approach to numerical simulation of 2D and 3D natural convection melting problem can be used for analysis of effective cooling system for electronic devices and new versions of the thermal energy storages. Comparison between 2D and 3D cases has shown that at initial time level melting in 2D occurs intensively, while with time this behavior changes. Also local Nusselt number along the heater surface is greater for 2D case in comparison with 3D data.

ACKNOWLEDGMENTS
This work was conducted as a government task of the Ministry of Education