Solid-Liquid Phase Change in the Presence of Gas Phase: Numerical Modeling and Validation

In the present work a fundamental numerical model is implemented in the open source framework, finite volume method (FVM) based computational fluid dynamics (CFD) program OpenFOAM®, to tackle the solid-liquid phase change phenomena in the presence of gas phase. The volume of fluid (VoF) method is employed to distinguish between the phase change material (PCM) and the gas phase, and the enthalpyporosity method is adopted to capture the moving boundary phase change phenomena in the PCM. The discussed model is validated against the well documented cases in the literature i.e., 1D Stefan problem, melting of gallium and free surface thermocapillary flow with phase change. The simulation results are in very good agreement with those found in the literature. Finally, experiments are performed with partially filled paraffin wax (Rubitherm® RT31) in the presence of air in an enclosure. The left and right walls of the enclosure are made of aluminum plates and maintained at different temperatures while the remaining walls are made of Plexiglas. The temporal evolution of the melting front during melting obtained from the discussed numerical model is compared with experiments.


Introduction
The solid-liquid phase change in the presence of third phase (liquid or gas) is of the utmost importance in many industrial processes such as casting, welding, thermal storage systems, powder atomization, etc. Therefore, the understanding of solidification and melting, its initial and thermal conditions plays a crucial role to make the process more efficient. Over the past few decades numerical simulations of solid-liquid phase change is gaining popularity in the industry not only to improve the quality of their products, but also to make cost effective manufacturing decisions [1]. A common approach for modeling such multi physical phenomena is coupling of Volume of Fluid (VoF) method with a solidliquid phase change model in Eulerian frame.
Since the introduction of the enthalpy method by Voller et al. [2], several numerical methods, formulations and models have been found in the literature to tackle the solid-liquid phase change phenomena and especially in handling of the velocity field in the solid phase. Hu and Argyropoulos [3] reviewed various methods and stated the pros and cons of the discussed methods according to the nature of the problem and computational expenses. Rösler and Brüggemann [4] developed a numerical model to study the melting process of phase change material (PCM) in latent heat storage systems and incorporated an error function to describe the melt fraction. Miehe [5] implemented and compared the viscosity and Darcy model to investigate the twin-roll casting process of the magnesium alloy. Kasibhatla et al. [6] adopted the variable viscosity model to study the melting and settling of PCM in a capsule. More recently König-Haagen et al. [7] performed a comprehensive benchmarking of fixed-grid methods and presented a qualitative and quantitative analysis of most common energy formulations for the modeling of melting.
The VoF approach proposed by Hirt and Nichols [8] is the well-known interface-capturing technique, which uses an indicator function to mark different fluids. This approach has also been previously used for the solid-liquid phase change simulations in a multiphase system. Assis et al. [9] performed parametric CFD analysis of melting process of PCM (Rubitherm® RT27) in a partially filled spherical shell. The VoF method was employed to account for PCM-air system and the enthalpy-porosity method was adopted for the solid liquid phase change. A temperature dependent density function ϱ(T) was used for air phase, whereas a variable density based on Boussinesq approximation was implemented in the liquid PCM.
Shmueli et al. [10] investigated the influence of the large mushy zone constant C in vertical cylindrical tube, partially filled with PCM (RT27) with a flat surface. A temperature dependent viscosity function μ(T) and variable density ϱ were used in the liquid PCM. Two sets of thermophysical properties were examined, at first uniform values of specific heat capacity (c pl = c ps ) and thermal conductivity (κ l = κ s ) for solid and liquid phase was assigned, later a small variation of thermal conductivity between phases was considered (i.e. (κ l ≠ κ s )). Rösler [11] investigated close contact melting of a PCM with air above it in a capsule. The PCM was assumed to be an incompressible fluid, whereas the gas phase (air) was modeled as a compressible fluid. Similarly, Kasibhatla et al. [12] also investigated the close contact melting in a capsule using a multiphase approach and variable viscosity model to enforce the velocity in the solid region to zero.
Furthermore, Tan et al. [13] performed simulations in small scale geometries to investigate the effects of thermocapillary convection and the dynamics of the gas phase during the solidification of pure bismuth with an argon gas on the top using a multiphase approach. Saldi [14] studied the Marangoni driven free surface hydrodynamics in liquid steel weld pools during solid-liquid phase change with free surface deformation. Finally, Richter et al. [15] performed mold filling and solidification simulations; the gas phase was treated as an incompressible fluid with constant properties. A special treatment was employed to handle the Darcytype source term at the interface between PCM and air.
As evident from the above literature research, the coupling of VoF Method with the enthalpy porosity method is an effective way to tackle such multi-physical phenomena numerically. Therefore, in this study a fundamental numerical model of solid-liquid phase change in the presence of gas phase is presented. The PCM and gas phase are considered as two-phase immiscible incompressible fluid and are modeled using the VoF [8] approach based on interface capturing method. A fixed grid enthalpy-porosity method [16,17] is employed for the solid-liquid phase change. Thermophysical properties of both the fluids are constant and independent of temperature. The Boussinesq approximation is used for both the fluids (liquid PCM and air).

1. Governing Equations
The solid-liquid phase change with a free surface is handled by solving the governing equations discussed in this section, using a single domain formulation. The conservation equations accounts for both the differences in fluid properties as well as surface tension forces acting at the interface. The mass, momentum and energy conservation equations are given by [13,14,16,17]: where = ( , , ) is the velocity vector, t the time, p the pressure, the density, the dynamic viscosity, ℎ the specific enthalpy, T the temperature, the thermal conductivity, the buoyancy source term, the surface tension force term, the Darcy type source term.

Volume of Fluid (VoF) method
The VoF uses an indicator function α to mark different fluids: in the present work fluid A is the PCM and fluid B is the gas phase. The volume fraction α transport equation can be written as: where r = A − B , is the relative velocity between two fluids. On the contrary to the classical VoF method originally proposed by Hirt and Nichols [8], the VoF method in OpenFOAM® uses an additional convective term also referred to as the artificial compression term (the third term in Eq. (5)) to minimize the numerical diffusion at the interface [18]. This term is only active within the cells containing interface due to its dependence on (1 − ) where ∈ (0, 1). The relative velocity also known as the compression velocity r is computed by: where φ is the fluid flux, surface area and C α compression coefficient. In the present work C α is set to 1. The thermophysical properties of each fluid are expressed in terms of volume fraction α: where Ψ is either , , , . The surface tension force term S σ in Eq. (2) is the continuum surface force (CSF) model proposed by Brackbill et al. [19] and is expressed as: where σ is the surface tension and the surface curvature. The surface curvature is defined as:

3. Solid-liquid phase change
The specific enthalpy ℎ in Eq. (3) can be defined as the sum of sensible heat and latent heat, such that: where c p is the specific heat, T Ref the reference temperature, α restricts the term in fluid A (PCM), γ(T) the melt fraction as a function of temperature and L the latent heat of fusion.
By combining Eq. (3) and Eq. (10), we get the energy conservation equation in terms of temperature, as shown below: The melt fraction γ(T) ∈ [0, 1] is defined as a linear function of temperature: where sol is the solidus temperature and liq the liquidus temperature.
The Darcy type source term D in Eq. (2) is modeled using the Carman-Kozeny equation for flows in porous media [16] and is defined as: where is the large constant that describes the morphology of the melting front in the mushy zone and = 1 × 10 -3 is the small computational constant to avoid "divide by zero" error. In Eq. (13), when = 0 the source term dominates in the cells containing solid phase and enforces the velocity in that cell to zero, and when = 1 the source term vanishes in the cells containing liquid phase.
The thermophysical properties of each phase are expressed in terms of melt fraction 1 , as shown in Eq. (14), where is either , or . The buoyancy source term in Eq. (2) is described by the Boussinesq approximation, such that: where = ( , , ) is the gravity vector. Figure 1 illustrates an overview of the complete model and is implemented in the open source FVM based CFD software OpenFOAM®. For the pressure-velocity coupling, the PIMPLE (PISO and SIMPLE) algorithm is employed. The VoF method is solved explicitly using a modified FCT (Flux Corrected Transport) method called MULES (Multidimensional Universal Limiter with Explicit Solution) and is based on the work of Zalesak [20] with an iterative calculation of the weighting factors [21]. A first order implicit Euler scheme is employed for the temporal discretization. For the gradient terms a second order Gauss linear scheme is applied. A variety of HRS (High Resolution Schemes) discretization schemes are used to solve the divergence terms of the equations.

1. Stefan problem
The one-dimensional one-phase Stefan problem is a well-known transient heat transfer problem for directional solidification [1]. Figure 2 shows the schematic diagram of the problem.  Table 1. The analytical solution for the determination of interface location x m (t) at time t is given by: where D is the thermal diffusivity, and is the nonlinear eigenvalue and is determined by solving the below transcendental equation using iterative methods.  adiabatic; an empty 2 boundary condition is applied on the front and back boundaries of the domain. At time t = 0, solidification commences from the left boundary and continues to propagate towards the right boundary. Figure 3 shows a comparison of numerical simulation and analytical solution for the one-dimensional Stefan problem. The interface locations predicted by the simulation are in excellent agreement with the analytical solution.

2. Melting of gallium
The melting of gallium in a rectangular enclosure is a well-documented solid-liquid phase change experiment performed by Gau and Viskanta [22]. The computational domain is shown in Figure 4, the size of the domain is 88.9 mm × 63.5 mm and is divided into 42 × 32 hexahedral numerical cells. The left and right walls are maintained at constant temperatures h = 311.15 K, c = 301.45 K respectively. The top and bottom walls are assumed to be adiabatic. The initial temperature of the system is init = 301.45 K. The thermophysical properties of gallium used in the simulations are given in Table 2 Furthermore, the mushy zone constant C and the computational constant ε in Eq. (13) are set to 1.6 × 10 6 kg m −3 s −1 and 1 × 10 −3 respectively, see Brent et al [16]. The temporal evolution of the melting front obtained from simulation exhibits good agreement with experiments of Gau and Viskanta [22] as shown in Figure  5. The slight deviations between the simulation and experiments can be explained as follows: a. the melting fronts are obtained from twodimensional simulation, whereas the experiments performed by Gau and Viskanta [22] were three-dimensional (for further information refer [5,23]), b. temperature independent thermophysical properties are used in the simulation, c. the experiment was stopped multiple times and the melt was drained in order to track and record the melting front, which may have influenced the results.

Free surface thermocapillary flow with phase change
A free surface thermocapillary flow with phase change in a two-dimensional rectangular container under microgravity conditions investigated by Tan et al. [13] is used for validation. For the sake of validation, the source term S σ in Eq. (8) is resolved into the normal f σ,n and the tangential f σ,t components [13,14,24] and is expressed by: where is the surface normal vector and the thermal dependence of surface tension. Figure 6 shows the computational domain of size 15 mm × 5 mm and is discretized into 140 × 62 hexahedral cells [13,14]. The bismuth melt is filled up to height Bi = 4 mm and the rest is filled with argon gas. The left and right walls are kept at constant temperatures h = 552.55 K and c = 540.55 K respectively. The top and bottom walls are enforced with a linear temperature distribution from h to c , such that, the isoline of melting temperature m = 544.55 K is at a distance m from the left wall. The dimensionless numbers corresponding to bismuth and argon gas are listed in Table 3 [13]:  [13], and the simulations by Saldi [14] is shown in Figure 7. An excellent agreement is achieved in comparison with literature.

Experimental setup
Experiments are conducted to further investigate the introduced model. Figure 8 illustrates the schematic of the experimental setup adapted from [11,22,26]. The domain consists of a square enclosure with inside dimensions of 50 mm × 50 mm × 50 mm made of Plexiglas. The left and right walls are made of aluminum plates which serve as heat source and heat sink respectively. A couple of multi-purpose riser pipes are screwed on the top wall of the enclosure which allows free flow of air to compensate for volume expansion/contraction during phase change. The enclosure is initially filled with paraffin wax (RT31) up to a height of 40 mm in liquid form through one of the riser pipes and is allowed to solidify for a certain period of time. All the walls are insulated except for the front and back walls. A DSLR camera is placed at an appropriate distance from the front wall to capture the sharp images of the temporal evolution of the melting front. A LED diffuser is placed at a suitable distance from the back wall in order to differentiate the solid and liquid phase during the melting process of the wax. Finally, all the images are processed and automated using ImageJ [25] for the evaluation of melting fronts and melt volume fraction of the paraffin wax.

Results and discussions
The thermophysical properties of RT26 to RT35 provided by the manufacturer [27] are comparable, the specific heat capacity and the thermal conductivity of the above mentioned PCMs are same and have uniform values for each (solid/liquid) phase. Shmueli et al. [10] compared their numerical results with experiments using the thermophysical properties provided by the manufacturer with a slightly varied thermal conductivities for solid and liquid phase. The thermophysical properties used in this study are adapted from the above literature and are tabulated in Table 4. The grid convergence study is performed to investigate the influence of cell size on the melt fraction. Different grid resolutions of 2D computational domain have been investigated, from coarse (200 × 200 × 1) cells to fine (400 × 400 × 1) cells increased in steps of 25 cells in each x and y direction. A grid resolution of (275 × 275 × 1) is selected in this study. The relative error of the melt fraction between the selected grid to the finest grid resolution is less than 2%. An adaptive time stepping with Courant Number ( ) of 0.75 is adopted. Figure 9 shows the experimental and numerical, temporal evolution of the melting front at 15, 30, 45 and 60 min respectively, for the temperature at hot and cold walls h = 313.15 K and c = 295.15 K respectively. The liquid paraffin close to the hot wall rises upwards and flow horizontally towards the cold melting front, where it releases heat at the upper region of the solid paraffin and descends along the cold melting front, thereby advancing the melting front rapidly in the upper region than in the lower region. The streamlines in the liquid paraffin shows that there is no movement approximately at the center of the melt, this is because of the upwards and downwards flow of liquid paraffin close to hot wall and cold melting front (indicated by solid arrows in the bottom row of Figure 9). It is observed that the air flow above liquid paraffin also progresses along with the melting front in a circular motion 3 . A secondary air flow is also observed in the region above the solid paraffin. However, the strength of this secondary air flow is minuscule. In the experiment it is notable that the liquid paraffin-air interface is thicker than in simulations. This is because of the contact angle between liquid paraffinair interface and the Plexiglas wall, which leads to the formation of meniscus near the Plexiglas wall. Furthermore, the melting fronts in the simulation appear to be sharper than in the experiment near the paraffinair interface. This is because in experiments the melting fronts are obtained by the evaluation of photographs taken through front Plexiglas wall. Nevertheless, it is noteworthy to mention that the surface tension effects are neglected during simulations. Figure 10 shows the comparison of melting fronts between experiments and simulation. The height of the liquid paraffin in the experimental domain varies due to the volume expansion during phase change. Therefore, for simplicity the temporal evolution of melting fronts has been normalized to the width of the experimental domain in x direction and to the height of the liquid paraffin in y direction. The melting fronts obtained by the numerical model are in good agreement with experiments.
However, the melting fronts at the lower region advances faster in experiments than in simulations, this is because the thermal conductivity of the Plexiglas walls also has an effect on the melting of the paraffin, however the influence of the Plexiglas walls is not considered during simulations.

Conclusions
In this study, a fundamental numerical model to tackle the solid-liquid phase change phenomena in the presence of gas phase has been investigated by means of well documented case studies in the literature and experimental validation. The model is based on the use of VoF approach to account for PCM and gas phase and enthalpy-porosity method to account for solid-liquid phase change in the PCM. The model has been validated using multiple references cases. The simulated results are in good agreement with the results described in the literature. To further investigate the limitations of the model, experiments with paraffin wax (RT31) are conducted. Experiments are performed in a square enclosure partially filled with paraffin wax in the presence of air. The left and right walls of the enclosure are maintained at temperatures h = 313.15 K and c = 295.15 K respectively, whereas the remaining walls are assumed to be adiabatic. The temporal evolution of the melting front obtained from simulation show a good agreement with the experiments. The small discrepancies in the results can be accounted for, as follows:  In experiments it is observed that the side Plexiglas walls also influence the melting phenomena within the enclosure, whereas in simulations these walls were not considered and are modeled as perfectly insulated.  Temperature independent thermophysical properties are used in the simulations for PCM (solid/liquid) and gas phase.  Only 2 dimensional simulations are carried out and the influence of the third dimension is not yet investigated.  Finally, the errors in evaluation of photographs obtained from experiments due to volume expansion and surface tension effects.