Influence of Natural Convection and Volume Change on Numerical Simulation of Phase Change Materials for Latent Heat Storage

For the numerical simulation of a heat storage based on phase change materials (PCMs) an enhanced model is presented, considering the physical effects of convection flow in the liquid phase as well as the volume change during phase change. A modified heat capacity is used to realize the phase change. The phase change material is initially defined as a liquid with temperaturedependent material properties. A volume force is added to the Navier-Stokes equations to allow a circulating flow field in the liquid phase and prevent flow motion in the solid phase. The volume change is implemented with the Arbitrary Lagrangian-Eulerian method. A laboratory phase change experiment was performed using the PCM RT42 with a melting temperature of 42 ◦C. The laboratory experiment was calculated numerically using the enhanced model to evaluate the numerical model and to investigate the influence of the simulation parameters on the thermal behavior of the PCM. The thermal conductivity is determined as the main influencing parameter. A good agreement of the simulated melting front throughout a major part of the laboratory experiment has been shown. COMSOL Multiphysics provides a default model for phase change, which neglects convection flow and volume change. Compared to the default model, the enhanced model achieves more accurate results but requires more computational cost for complex latent heat storage systems. Using the default model without convection can be reasonable, considering that the heat storage design is either over-dimensioned or a suitable correction factor must be applied.


Introduction
Phase change materials are well-suited for heat storage and heat management systems since latent heat is absorbed or released at a constant temperature during the phase change. Therefore, a significantly higher heat quantity is available compared to mere sensible heating. Technical applications of PCM as heat storage are, for example, air conditioning or peak load reduction. PCMs can also be used as buffer storage for solar systems and waste heat utilization, as well as for passive thermal management of battery systems [1,2].
Numerical simulation is an important tool for designing systems with PCM in order to perform optimization without time-consuming and cost-intensive tests. Numerical models of a PCM during phase change require high computational cost due to the occurring physical effects. To reduce the computational effort, simplifications can be made by neglecting certain physical effects. This study investigates the influence of the most used simplifications in numerical phase change models on the accuracy of the numerical results by comparing the results of a laboratory experiment with two different numerical models.
The physical effects during the phase change of a PCM are shown in Figure 1 as a schematic illustration. It shows the section through a cube filled with solid PCM, which is heated from one side with a heat source. The melting occurs starting from this heat source as soon as the temperature of the PCM exceeds the melting temperature. The illustrated physical effects during phase change apply to PCMs regardless of their specific melting temperature. Heat conduction occurs through the entire system towards the lower temperature. Besides, free convection and heat radiation to the environment takes place. Small differences in density in the liquid phase of the PCM occur due to different temperatures. Hence, a convection flow in the liquid fraction of the PCM develops. In this example, a clockwise convection flow is created, resulting in a faster melting of the top areas. A significant change in volume of up to 15 % occurs during the phase change process of a PCM due to the differences in density in the solid and liquid phases. In [3], the modeling of the volume change during phase change is presented with the software ANSYS Fluent with the Volume-of-Fluid method, which is used with a finite-volume discretization. This approach cannot be implemented with the software COMSOL Multiphysics, since COMSOL is based on the finite element method. The volume change is neglected in all present studies regarding numerical modeling of phase change with the software COMSOL Multiphysics. The Marangoni effect describes a temperature-dependent surface tension, resulting in a fluid flow along the surface. This effect can be neglected [4]. Accounting for the effect of the natural convection in the liquid phase requires a coupled system of equations for heat transfer and fluid flow. This results in large computational time [3]. One of the first numerical investigations of phase change with consideration of natural convection was performed by Voller et al. [5]. The numerical model must allow a fluid flow in the liquid phase of the PCM and prevent the flow in the solid part. An easy approach would be to increase the dynamic viscosity in areas with temperatures smaller than the phase change temperature to values which correspond to a solid material. The work of [6] shows that no convergence can be achieved with this approach using COMSOL Multiphysics. The so-called enthalpy-porosity method describes a solution for implementing phase change with natural convection and was first published in [7]. The method is applied in various recent work, for example in [8][9][10][11][12][13] as well as in this paper.
The software COMSOL Multiphysics includes a PCM module which neglects convection flow and volume change. This model will be called default model. However, the physical effects of convection and volume change have a decisive influence on the melting process, as shown above.
In this study, an enhanced numerical PCM model for COMSOL Multiphysics 5.5 is presented, which takes the convection flow in the liquid phase of the PCM and the volume change into account. This model will be called enhanced model in the following. An evaluation of the enhanced model is performed by a numerical recalculation of a performed laboratory experiment. The results are compared to the default model and the application of the numerical models for designing heat storage with PCM is discussed.

PCM Selection
The PCM RT42 with a melting temperature of 42 • C is used. The manufacturer's specifications of the Paraffin RT42 [14] and its thermal volume expansion coefficient [15] are listed in Table 1.  Figure 2a shows the setup of the laboratory experiment consisting of an open upward cylindrical container filled with the PCM RT42, an outer wall made of plexiglass and an aluminum bottom plate. A steel pipe runs vertically centered within the housing and is connected to an upward water flow. The water is tempered to approximately 60 • C at the inlet. The PCM has a height of 18 cm, the air layer has a height of 8.5 cm, and the outer radius of the enclosure is 2.9 cm. In preparation for the laboratory experiment, the PCM is completely melted and cooled back down to room temperature to minimize air inclusions, which are a result of filling the solid PCM in the cylinder. The ambient temperature is constant at 20 • C. In the beginning, the PCM is entirely in the solid phase. During the melting process a measuring template is used to determine the position of the melting front by measuring distances from the top surface of the PCM at specific times. Only the melting front, which is accessible from the top opening of the cylinder, is measurable with this measuring method. The time for a complete phase change is 255 min and a final volume expansion of 3.2 % is measured afterwards. The manufacturer specifies a volume change of 13.6 % for a temperature range between 15 • C and 80 • C, whereas a significantly smaller temperature range is investigated in the laboratory experiment between 20 • C and 60 • C. Nevertheless, the discrepancy to the material data is large and the volume change is smaller than expected. The results of the experiment such as the melting front are described later. Figure 2b shows the setup of the numerical model. The geometry represents the laboratory experiment, which is modeled two-dimensional and rotationally symmetric. The water flow is represented by a time-dependent temperature boundary condition, which corresponds to the measured and averaged water temperature in the laboratory experiment. To recalculate the laboratory experiment, a density difference is used which corresponds to the measured volume expansion of 3.2 %. Therefore, the assumption of a linear characteristic of density is made between the given values of 15 • C and 80 • C. The influence of different densities while maintaining a constant density difference for the solid and liquid phase is investigated. No significant differences in the results were found. Hence, the densities were set to 825 kg/m 3 for the liquid phase and 852 kg/m 3 for the solid phase corresponding to temperatures of 45 • C for the liquid phase and 30 • C for the solid phase.

Numerical Model
COMSOL Multiphysics provides modules for modeling heat transfer and fluid flow which were used in this study. The boundary conditions correspond to the laboratory experiment. The heat transfer coefficient α describes the heat flux from surface areas to the environment. Free convection and radiation are considered according to DIN EN ISO 6946, and in this case α is defined as 7.69 W/(m 2 ·K).
A Navier-Slip wall condition was applied to the horizontal outer walls of the PCM as required in [16] for walls next to moving meshes. For the computation, a surface tension σ needs to be specified and the values for Paraffin are between 19 mN/m and 35 mN/m [17]. Preliminary studies have shown that the impact of the magnitude of the surface tension in the given range is neglectable. Thus, the surface tension σ is set to 30 mN/m. Slip wall condition was applied to all remaining walls.
The mesh is generated with triangular elements. The phase transition region needs to be resolved by a sufficient amount of small elements which leads to high computational cost. A grid independency study was performed for three different mesh sizes regarding the overall melting time. The results are shown in Figure 3. The grid with 34,697 elements is chosen due to its good results with a reasonable computation time.
Melting is a time-dependent process, so therefore a transient numerical modeling is necessary. Adaptive time stepping and a fully coupled approach was used, which solves the nonlinear system of equations as a single large system of equations within one iteration [18].
The fundamental theory described in the following section can be found in [18]. The description of the temperature field follows from conservation of energy to Equation (1). The numerical model for solving the flow field is based on the basic equations of fluid mechanics according to [19]. The Navier-Stokes equations for an incompressible flow with a Newtonian fluid are shown in Equations (2) and (3). An additional term −S(T) · u, which is described later, is necessary to model the transition area between the solid and liquid phases.
The buoyancy force F B causes the natural convection in the liquid phase and is calculated using the Boussinesq approximation according to Equation (4), with the magnitude of the acceleration of gravity g of 9.81 m/s 2 .
The phase change is characterized by the melting fraction ϕ according to Equation (5). The parameter phase change temperature range ∆T pc defines the range of the transition area.
To describe the convection flow in the liquid phase, the so-called enthalpy porosity method is used to allow convection flow in the liquid phase while preventing flow in the solid PCM. Therefore, an additional source term S(T) is defined according to Equation (6), which contains the Carman-Kozeny relation. The so-called mushy zone factor Am is important for numerical modelling and its magnitude has a decisive influence on the results and is described below. The constant = 0.001 prevents division by zero for the case ϕ(T) = 0.
The source term S(T) approaches zero for liquid areas with a melting fraction of ϕ(T) = 1. This results in a flow corresponding to the standard Navier-Stokes equations. Solid areas and the transition area have a melting fraction ϕ(T) < 1, which leads to a great value of the source term S(T) depending on the magnitude of the mushy zone factor Am. Therefore, the fluid flow is prevented in the solid area and reduced in the transition area between solid and liquid phase.
The best results are achieved by using great magnitudes for the mushy zone factor. Preliminary studies have shown that there is no significant change in the results for A m ≥ 1·10 4 kg/(m 3 ·s) for incompressible flows. In this paper, the mushy zone factor A m = 1·10 4 kg/(m 3 ·s) is used, thus convergence is achieved for compressible flow.
The volume change is solved using the Arbitrary Lagrangian-Eulerian method. It describes the motion of a surface which divides the air and the PCM [18]. In COMSOL Multiphysics, this can be added by the "Moving Mesh" module with a deformable geometry for the two domains PCM and air.
The physics of phase transformation is implemented according to [11] with userdefined functions and material properties. The PCM is treated numerically as a liquid. User-defined functions are used to describe the material properties of the solid and liquid PCM as well as the flow in the transition zone. The melting fraction ϕ(T) is used to describe temperature-dependent material properties. The density ρ(T) is obtained according to Equation (7) and the thermal conductivity λ(T) similarly, according to Equation (8).
To define the phase-dependent viscosity µ(T), the source term S(T) is used as well as an auxiliary constant ξ = 1 m 2 , as shown in Equation (9).
The heat capacity method defines the phase change using a modified specific heat capacity. The latent heat is the characteristic quantity for the phase change. A Gaussian function according to Equation (10) is used, which is centered around the phase change temperature T pc .
The modified specific heat capacity c p (T) is calculated according to Equation (11), maintaining the energy balance during phase change.
The validation of the enhanced model was done by recalculating the experiment presented in [20]. Figure 4 shows a good agreement of the overall melting time.

Results
The numerical results of the velocity v, the temperature ϑ, and the melting fraction ϕ are shown in Figure 5 for the enhanced model. Only the PCM domain is visualized. The influence of the convection flow in the upper region of the PCM can be seen in the curved phase boundary between the solid and liquid phase, resulting in a faster melting of the upper regions of the PCM. In addition, due to the high thermal conductivity of the aluminum bottom plate, melting starts at the bottom. The volume expansion of the PCM into the air domain results from the density difference between the solid and liquid PCM. Additionally, the measured position of the melting front of the laboratory experiment for three different time points are marked with white triangles in Figure 6. The number of measured points differ between different times due to the measurement method described in Section 2.2. Melting starting from the bottom plate is represented by one measuring point at the bottom. The position of the melting front in between the top and bottom measuring points is not measurable. Until 120 min approximately 75 % of the PCM is melted and there is a good agreement in the position and shape of the melting front between the laboratory experiment and the simulation. The melting of the last 25 % of the PCM required 135 min in the laboratory experiment but only 65 min in the simulation.
A parameter study was performed to investigate the influence of the following simulation parameters on the overall melting time t pc : the specific heat capacity c p,s and c p,l , the phase change temperature range ∆T pc , the thermal expansion coefficient β, the latent heat L, the heat transfer coefficient α, and the thermal conductivity λ. The investigated parameters are individually increased or decreased by 20 %. The change of the melting time ∆t pc compared to the enhanced model with a melting time t pc of 184 min is studied. The latent heat L, the heat transfer coefficient α, and the thermal conductivity λ were found to have the greatest influence, which is shown in Table 2. The values used for the latent heat and the heat transfer coefficient can be considered as reliable data. The thermal conductivity λ has the greatest influence on the melting time t pc . In general, there is a temperature and phase dependency for the thermal conductivity [21], which is unknown for the used PCM RT42. Moreover, small air inclusions occur in the PCM during the solidification process in the laboratory experiment, which leads to a decreased thermal conductivity. Therefore, a deviation of the manufacturer's data from the real thermal conductivity is possible. In this case, fast charging of the heat storage is achieved by melting 75 % of the PCM and an insulation around the heat storage is recommended additionally.
Another numerical calculation of the laboratory experiment is made using the default model, which is available in COMSOL and neglects convection flow and volume change. Figure 7 shows the numerical results of the default model (a) and of the enhanced model (c) after 120 min. In the default model the melting front is vertical to the heat source, and no complete melting of the PCM occurs. This is a discrepancy to the results of the laboratory experiment and to the enhanced model. Therefore, a latent heat storage designed with the default model leads to an over-dimensioned thermal energy storage system. The default model and the enhanced model show the same results during the solidification process.

Discussion
The numerical calculation of a performed laboratory experiment shows that the physics of the phase change is well represented with the enhanced model until 75 % of the PCM has been melted. The melting duration of the last 25 % of the PCM is not represented correctly in the enhanced model. Therefore, in this setup it is recommended to avoid melting the last 25 % of the PCM and to add an insulation to achieve faster melting. The limitation of the enhanced model is large computational costs. The implementation of the enhanced model leads to an increase of the computational time by a factor of 60 compared to the default model due to complex flow calculation.
The default model is suited for the design of heat storages with large dimensions because it is easy to implement for complex geometries and has low computational costs. For melting, the melting fraction is represented inaccurately due to the neglection of the convection flow in the liquid PCM and the volume change. This leads to the design of over-dimensioned latent heat storage with the default model. To compensate for the neglected physical effects of the default model, the use of a correction factor is recommended to avoid an over-dimensioned storage. The aim of using a factor is a good agreement of the average melting fraction of the default model and the laboratory experiment. The thermal conductivity is suitable because it has the greatest influence on the melting process. The value of the factor is not universal but depends on the operating conditions like the PCM and the geometry. A further parameter study found that a multiplication of the thermal conductivity by the factor 2 results in a good representation of the laboratory experiment with the PCM RT42. Figure 7b shows the melting fraction of the default model with factor 2 after 120 min. The increased thermal conductivity leads to faster melting. This will not result in a correct shape of the melted PCM, but it compensates the neglected influence of the convection flow on the heat transfer rate. As a result, the melting fraction of the system and the required melting time can be calculated. Further investigations of the correction factor for PCMs with different melting temperatures are suggested. Acknowledgments: This work has been supported by the Fraunhofer Institute for Manufacturing Technology and Advanced Materials IFAM in Dresden (Fraunhofer IFAM Dresden) within the project »PCM90Plus« as part of a diploma thesis. Thanks go to the colleagues at Fraunhofer IFAM Dresden for their valuable support.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Vector quantity