Numerical simulation of a phase change material melting process

. Storage processes are usually integrated in solar energy systems applications due to daily variation of this energy source availability. Among different thermal storage solutions, phase change materials (PCM) lately became more extensively used covering a wide range of operating temperatures. In this regard, a numerical simulation of a PCM melting process is performed under ANSYS CFD environment. A particular configuration is considered consisting in a 2m length annular tube having a 5.48 cm external diameter. The tube is filled with paraffin chosen as PCM. A concentric interior tube of 2.54 cm diameter is used for transporting the heat transfer fluid (HTF) from the solar collector. Heat is transferred through the 1mm thick pipe wall to the PCM placed all around the HTF tube. The numerical results reveal the melting process of the PCM at different instances and tube sections. The time variation of the PCM liquid fraction is emphasized. The results describe the dynamic behavior of a PCM melting process and might be further integrated in any solar power plant storage charging process simulation.


Introduction
Energy storage represents an essential component of solar energy systems, that lately became more and more used to transform thermal energy of Sun into mechanical, electric and/or thermal one to the user. The solar energy source has the important advantage of being free and renewable. Nevertheless, it has a certain variability, not only from day to day, but even along the same day. More than that, cloudy days could affect the system predicted operation due to lack of input energy over some periods. Thus, a backup in terms of thermal energy source is required in order to maintain the system continuous operation at desired parameters. The main role of a thermal energy storage component (TES) is to ensure continuously the right levels of temperature and heat flux input to the converting system. Secondly, this component improves the overall system operation, maintaining constant parameters in operation, extending the operation interval and improving the overall system efficiency and reliability.
Thermal energy may be stored either by storing the sensible heat, as water tanks or rock bed, or the latent one, as phase change materials (PCM) use. The latest is a viable and efficient solution since it offers a high energy density and a constant temperature level corresponding to the phase change process.
The current literature offers an important number of publications reporting PCM research and results. An extensive review of 100 papers on TES is reported in [1]. More than one third is dedicated to PCM researches on modelling and experimental results. Sharma et.al. [2] presents a review containing 155 papers only on PCM systems and applications. Stages in PCM systems development, properties and criteria for PCM selection and their classification from different points of view are presented. As grouped in organic, inorganic and eutectic, PCM are employed over a wide range of temperatures. Among the most important criteria for PCM choice, one may consider the melting point, latent heat and volumetric latent heat storage capacity, as these are constraining the system operation parameters and size. Organic PCM are paraffin wax and non-paraffin. They are not corrosive, safe in use and have a very long freeze-melt cycle, so they are a reliable choice for the system. Their melting point varies from 5.5°C to 75.9°C for paraffin and reaches 127°C for non-paraffin [2] so that they cover a quite large temperature range of applications. The volumetric latent heat storage capacity for the organic compounds lies on 128-200 kg/dm 3 [2]. A double value is characteristic to inorganic PCM. These are salt hydrates and metallics. Their melting point varies between 14°C to 117°C for salt hydrates and even higher, 1300°C for carbonate salts [3], and 30°C to 125°C for metallics. The biggest problems encountered with inorganic PCM could be a more reduced chargedischarge cycle due to salt solubility in water decrease in time and necessity to discharge at a lower temperature than the fusion one, respectively the weight problem for metallics. Eutectic PCM represents a mixture of two or more compounds melting and freezing congruently. Their melting point varies between 13°C to 767°C for LiF/CaF2 (80.5/19.5 mol%) [3]. Alva et.al. [3] presents storage materials used in solar energy applications, for different temperature levels, and also taking into account the unitary costs, which lies between 0.1$/kg and 100$/kg. Different PCM encapsulations can be encountered, namely bulk storage in tank heat exchangers, requiring extensive heat exchange, macro-encapsulation and microencapsulation [1]. The macro-encapsulation consists in including the PCM inside a tube, sphere, cylinder, while the micro-encapsulation consists in including the PCM in a microsphere incorporated at its turn in a compatible material. The macro-encapsulation is the most widespread method [1]. More data on commercially produced PCM and geometries are reported by [4].
Regarding mathematical modeling of PCM systems, also known as latent heat thermal energy storage system (LHTES), a review is reported by Verma et.al. [5]. The authors presented a synthesis of proposed methods by grouping them in two categories, namely methods based on First Law of Thermodynamics, and respectively methods based on Second Law. It was emphasized that the methods based on energy conservation have been validated, most of them, by experimental results, while the other category based on exergy analysis takes into account also the effect of time duration through which heat is supplied and the temperature of heat in storage, but they were not validated with experimental data.
Other numerical modelling of PCM behavior are based on software simulations, such as TRNSYS or ANSYS CFD. Al-Abidi et.al. [6] presents a review regarding CFD applications for LHTES in electronic cooling technology, building thermal storage, heating, ventilation and air conditioning, and other similar areas. Spherical, rectangular and cylindrical capsule geometries are discussed. One dimensional (1D) and two dimensional (2D) models are overseen, and only one 3D model is presented. Among these, for paraffin PCM, there are two 2D models that used power law for discretization and SIMPLE algorithm adopted for pressure-velocity coupling for a cylindrical PCM, and one 3D model that used the backward Euler scheme for the temperature discretization in time for a PCM bags in duct geometry. Results as melt fraction over time are presented.
Nayak et.al. [7] presents a comparison between three PCM materials (paraffin wax, sodium acetate tri-hydrate and phenolphthalein) used for an automobile application, for cold start of the engine. In their study, the PCM material is filled inside the PCM container placed in the middle of the engine cooling water heat exchanger. The simulation is done with Gambit and Fluent software tools and validated experimentally. The charge and discharge processes are presented. For example, paraffin charging period is presented for 1300s, starting from 27°C and reaching 48.1°C, while coolant water temperature varies from 70°C to 59°C.
Another application of paraffin wax is reported by Dakhil et.al. [8] for a solar heat pump. The simulation was performed in Fluent, as a 2D model using SIMPLE algorithm for discretization and k-ε model of turbulence. The PCM properties were supposed independent on temperature (melting temperature of paraffin 51°C and latent heat of 140 kJ/kg), initially being in solid state close to the melting temperature. A comparison in terms of temperature plot versus time is presented for CFD and experimental results and authors reported good agreement between them.
Li et.al. [9] studied the melting process of a paraffin blend RT27 inside a horizontal annulus considering water inside the inner tube by applying a CFD simulation. Sensitivity studies were developed to determine the optimal geometric configuration.
Among PCM, paraffins are widely used in applications involving solar energy systems. The choice of the TES component is very important. In our previous studies, we have modelled an absorption cooling system (ACS) fueled by solar energy using water tanks for thermal energy storage. As the water tank temperature was very sensitive to ACS load, instability in system daily operation was emphasized. Thus, the importance of ensuring a constant temperature for the ACS desorber was found as being crucial not only for the system performances, but also for the system operation [10][11][12]. PCM storage could be used instead to maintain a constant input temperature to the ACS desorber.
Consequently, the target of the paper is to present the results for a numerical simulation of a PCM melting process performed under ANSYS CFD environment. A 3D model is presented for a tube filled in with paraffin and surrounded by an exterior one filled with a HTF (water) from solar collector. An inside understanding of the melting process is also performed, analyzing the advancement of the melting process in the PCM volume all along the PCM tube. Solid, mushy and liquid zones are emphasized at different sections along the PCM tube as time passes. Delays in melting time between two consecutive radial sections are put into evidence. The correct computation of PCM melting time in each section is essential for system evaluation accuracy, as it involves heat transfer phenomena.

System description
One of the most used geometries for a LHTES is shown in fig.1a and consists of many tubes, arranged in a regular manner under a pitch p, surrounded by a cylindrical shell. Inside de tubes the HTF is flowing. The space among the tubes and the cylindrical shell is filled in by the PCM. The number of tubes and length of the storage tank depend on the PCM quantity, which have to be used for energy storage. In the case of massive storage tank, the geometry of the system is too complicated for simulations. Frequently, one chooses a simplified geometry, consisting of only one tube and an annular enclosing piece of PCM, whose external diameter is D pcm = p.
The simplified geometry used in this simulation is presented in fig.1b. The inner diameter of the tube through which the HTF is flowing is Din = 0.0254 m, while the external diameter of the annular enclosure filled in with PCM has a value of D pcm = 0.0508 m. The length of the simplified storage system is L=2 m. Additionally, another five radial test sections were considered in order to analyze the melting process advancement. The first and the last ones are placed at 0.01 m for the inflow and outflow sections, respectively (z=0.01m and z=1.99m). The other three sections are positioned at z = 0.5, 1.0 and 1.5 m measured from the inflow section. Water has used as HTF, and paraffin-wax as PCM. The water inlet temperature and water mass flow rate are set to T htf,in = 60°C and ṁ htf =0.005 kg/s, respectively. The water flow is laminar and experiences the hydrodynamic and thermal developing phenomena on a certain length from the inlet section.
The melting latent heat and melting temperature of paraffin-wax are Δh sl =164 kJ/kg and T melt = 50°C respectively. The thermophysical properties of solid and liquid phases are shown in table 1. However, from numerical reasons, one has considered that the density of both solid and liquid phases have the same value,  = 818 kg/m 3 and the solid phase starts to melt at T ms = 49.5°C and finishes melting at T ml = 50.5°C. Between these temperature values, the liquid and solid phases coexist. The heat transfer process between the HTF and PCM is obviously unsteady. Since the dynamics of the melting process is investigated, the initial temperature values used for starting the computations are set at 47°C for both the HTF and PCM domains.

Mathematical model and numerical procedures
The 3D mathematical model relies on the continuity, momentum and energy equations for the HTF domain, written for laminar flow and on energy equation for the PCM domain. These equations were solved in the ANSYS FLUENT v14.0 environment where the solidification & melting model was selected. This model neglects the natural convection occurring in the PCM liquid phase but considers a mushy zone as transition between the liquid and solid PCM phases. The enthalpy-porosity technique treats the mushy region (partially solidified region) as a porous medium. The first order implicit time discretization scheme and second order upwind spatial discretization algorithm were employed in calculations. For pressure-velocity coupling, the SIMPLE algorithm was set. Fig. 2 presents some details of the final grid used in computations.

Fig. 2. Details of final grid used in numerical simulations.
This grid was obtained after establishing the independency of the numerical solution with respect to spatial discretization and contains about 627000 faces and 230000 nodes. The time step used for advancing in time the numerical solution was Δ = 20 seconds.

Results and discussions
Figs.3-6 present the temperature profiles and the distributions of liquid mass fractions in the considered sections of the system at certain time moments. Let analyze first the temperature distributions in HTF domain, and after that the dynamic of both, temperature change and melting process evolution occurring in PCM region.
In the entrance region of the duct, the HTF flow experiences hydrodynamic and thermal developing phenomena. Due to the higher Pr value of water, the thermal developing length is more important than the hydrodynamic one, which may represent an advantage concerning the PCM melting occurring in this first TES length. Along the time, tube walls temperatures slightly increase due to rise of the PCM temperature near the outer wall of the tube. As a result, the heat flux transferred from the HTF to PCM decreases as the time increases. Whatever the time is, the HTF temperature profiles at z = 0.01 and z = 0.5 correspond to thermal developing flow regions, which ends somewhere above this last section. On the upper radial sections, the temperature profiles reveal that the convection heat transfer takes place in the thermal developed regime, meaning that the convection heat transfer coefficient is constant.
In what follows, the time variation of temperature change inside the PCM domain and the time advancement of melting process is explained. After 600 seconds (see fig.3), the PCM temperature already increased in all testing sections over the value of the initial temperature (i.e. 320 K), but differently. As expected, the highest temperature gap is happening at z = 0.01 m, where the temperature lies between 332.5 K on the outer wall of the tube and 322 K on the outside diameter of PCM domain, approximately. As a result, in this section, the melting process is already developing. The PCM liquid phase extends until a radius of about r = 0.016m, where the temperature is 323.5 K. At lower radius values, the liquid phase of PCM is obviously superheated. At higher radius, the mushy zone extends until the PCM temperature reaches a value of 322.5 K, and above this radius, the solid region still exists. However, in all solid PCM region, the temperature increased above the initial temperature value.
One may observe that at z = 1.0 m, the melting process is in its early stage and the extent of mushy zone is lower than in the previous analyzed section, while at z = 1.99 m only the preheat of solid PCM is present.
As the time increases, heat spreads on the larger PCM area and the melting process advances, but its velocity decreases along the axis. After 1200 seconds, the temperature was further increased in all PCM volume and the liquid PCM phase lies until a radius of about 0.02 m at z = 0.01 m, and 0.0170 at z = 1.99m, respectively.
As shown in fig. 4, at t=8100 seconds, the melting is quite finished on z = 0.01 surface and has substantially advanced in all other considered sections. One also has to note that at z = 1.0 m only the liquid phase and the mushy zone exist, meaning that the melting process has already taken place in the most faraway area of PCM, whereas on z = 1.99 m a thin layer of solid PCM still exists.
As fig. 5 present, after a time of t = 10800 seconds, the melting was ended on about the first half (till z= 1 m) of the PCM volume. On the other half, excepting a short volume near the end of the tube, the mushy zone already reaches the PCM external surface. However, starting from this moment, one needs another 2100 seconds to fully melt the remaining part of solid and mushy PCM. As revealed by fig. 6, at this time the liquid PCM is entirely superheated, but the superheated degree decreases along the HTF flow direction.  The time dynamics of PCM melting process is shown in fig. 7. There, the area-weighted averaged of mass liquid fraction in considered test sections is plotted vs. time. The radial variation shapes are sharper on sections placed at lower axial coordinate than those existing at higher ones. The most important shape difference appears between the radial sections of z = 0.01 m and z = 0.5 m. It is due to the fact that, on this length, both the convection heat transfer coefficient and the bulk temperature of HTF decrease with the increasing length. Once the fully thermal developed regime of HTF flow is reached, only the decrease of HTF bulk temperature along the flow direction is responsible for heat transfer decrease, so that the time shapes of area-weighted averaged of mass liquid fraction in the following sections become closer to each other. the following two successive sections is reduced to 1056 seconds and further it decreases to 724 seconds and 650 seconds, respectively, as presented in Table 2. This phenomenon occurs due to the fact that the HTF convection heat transfer and HTF bulk temperature decrease along the flow direction.

Conclusions
The study performed in the present paper emphasizes the melting process of a paraffin at different radial sections and time moments. The simulation was performed in CFD environment, considering a 3D model consisting in two concentric tubes, the inner one containing the water as HTF, while the outer one contains paraffin as PCM. The results obtained put into evidence the advancement of the melting process over time, analyzing the dimensions of the solid, mushy and liquid regions at different radial sections along the tube and at different time moments from the start process. As expected, the melting process is not uniform all along the tube, as the HTF convection heat transfer and HTF bulk temperature decrease along the flow direction. The delay melting time between the inlet and outlet sections is considerably high, reaching a value of 4300s, which is more than one hour. Over this period, the melt PCM in the inlet section is superheated. The work can be further integrated in a solar energy conversion system, considering time dependent HTF temperature existing a solar collector and consequently simulating the stored energy and temperature level in the PCM component over time.