An Optimum Enthalpy Approach for Melting and Solidification with Volume Change

Classical numerical methods for solving solid–liquid phase change assume a constant density upon melting or solidification and are not efficient when applied to phase change with volume expansion or shrinkage. However, solid–liquid phase change is accompanied by a volume change and an appropriate numerical method must take this into account. Therefore, an efficient algorithm for solid–liquid phase change with a density change is presented. Its performance for a one-dimensional solidification problem and for the quasi two-dimensional melting of octadecane in a cubic cavity was tested. The new algorithm requires less than 1/9 of the iterations compared to the source based method in one dimension and less than 1/7 in two dimensions. Moreover, the new method is validated against PIV measurements from the literature. A conjugate heat transfer simulation, which includes parts of the experimental setup, shows that parasitic heat fluxes can significantly alter the shape of the phase front near the bottom wall.


Introduction
Melting and solidification processes receive scientific interest as they are fundamental to many phenomena in nature and technical applications. Generally, solid-liquid phase change extends over multiple length and time scales, but for most engineering purposes a macroscopic continuum approach is adequate. Examples are latent energy thermal heat storage [1,2] and metal processing [3]. Multiple approaches exist to solve solid-liquid phase change in the framework of continuum theory numerically. Mostly, these approaches are classified by the way they treat the moving phase boundary into deforming and fixed grid schemes [4,5]. It has turned out that fixed grid schemes offer more flexibility and are easier to implement than deforming grid schemes. Nevertheless, numerical methods for melting and solidification on fixed grids are still challenging [6]. The energy equation is highly nonlinear and most problems are intertwined with fluid flow and some even with solid movement [7]. Furthermore, an additional variable, called phase fraction, is required in fixed grid schemes. At the time when the first fixed grid numerical methods emerged, only diffusive phase change was considered [8,9]. Later, multiple approaches were developed to describe convective phase change on fixed grids [10][11][12]. Great effort was put into the development of efficient algorithms, especially by Voller and his coauthors [13,14]. This resulted in a numerical scheme, based on Newton linearization, with a minimized number of iterations. The numerical scheme was named the optimum approach and the idea of linearizing phase change problems has triggered further research [15]. However, back then computers were not nearly as powerful as today. To keep the computational effort under control, the volume change upon melting or solidification was disregarded and the variable density in the liquid was described by the Boussinesq approximation. This is still the standard today, despite a drastic increase in computational power and the fact that solid-liquid phase change is always accompanied by a volume change. Moreover, the validity range of the Boussinseq approximation is unknown for solid-liquid phase change.
Nowadays, researchers start to include density change in their simulations [16][17][18], e.g., Hassab et al. [19] studied the effect of volume expansion for different heating boundary conditions and Dallaire and Gosselin [20] the solidification of water in an elastic capsule. However, they all lack a customized algorithm for solid-liquid phase change accompanied by density change and have to use slowly converging methods. To tackle this problem, we developed a new algorithm for melting and freezing of phase change material (PCM) with density change, which minimizes the number of iterations required and is based on the optimum approach.
This paper describes the new algorithm and compares its performance to the source based method. The source based method converges slowly but is quite stable and often used. It is implemented in most commercial codes and is used as a benchmark. Thereby, the original source based method was adjusted to be able to solve for a density change upon melting and solidification. Test cases were a one-dimensional solidification problem, where an analytical solution exists, and a two-dimensional melting problem with natural convection. Furthermore, we compared the numerical solution of the two-dimensional melting problem to experimental velocity fields obtained by Faden et al. [21].

Mathematical Model
Liquid and Solid phase are treated on the same grid and are tracked by the phase fraction α. With the aid of the phase fraction, the velocity inside the solid phase is set to zero using a Darcy source term: where D is the Darcy constant and a small numerical constant to avoid dividing by zero. The flow in the liquid phase is assumed to be laminar and Newtonian. For a Newtonian fluid, the viscous stress tensor reads: where η is the dynamic viscosity. Furthermore, the PCM is assumed to be incompressible, hence density is solely depending on the temperature, i.e., ρ = ρ(T). Pressure and viscous terms in the energy equation are not taken into account. Considering these assumptions the resulting mass, momentum and energy equations are: where the mixture heat conductivity is defined as a linear blend between its solid and liquid value: This system of equations requires an equation of state to be closed. Mostly, this relation between the enthalpy and the temperature is not given directly but split up in a sensible and a latent part. To be able to use material properties from the literature, we follow this approach and define a piecewise enthalpy-temperature relation. Moreover, instead of a melting temperature, we define a melting range: where T S is the solidus, T L the liquidus temperature and L the latent heat of fusion. The solidus and the liquidus enthalpy are defined as: Additionally, the melting range is thin, so we can assume a constant heat capacity in the mushy region: This simplifies the inversion of the enthalpy-temperature relation considerably, which is necessary for the numerical solution later on: The relation between the enthalpy and the phase fraction follows from the mixture approach [22]: with Combining these two equations, the phase fraction can be expressed as a function of the enthalpy: Alternatively, α can be expressed as a function of the temperature: Apparently, the phase fraction is a non-linear function of the temperature. However, with this definition of the phase fraction, the enthalpy (Equation (7)) is still a linear function of the temperature in each interval.
It should be noted that the methods described in the next section work for arbitrary density and enthalpy temperature relations and not only for the specific cases presented here.

Numerical Solution
The mathematical model was implemented into the finite volume, CFD toolbox OpenFOAM (version 2.2.2) [23], whereby the solver was based on the standard OpenFOAM solver buoyantPimpleFoam. The system of partial differential equations was solved in a sequential manner, whereas strong couplings were treated through inner iterations and weak couplings through outer iterations ( Figure 1). The focus here was on the T-h coupling, hence, it was assumed that the mass flux is known. A detailed description on the pressure velocity coupling is given in the Appendix A.

Optimum Approach
The original optimum approach was developed for phase change with sensible enthalpy transport and constant thermophysical properties. Its main idea is to update the enthalpy directly on the T-h curve. The starting point is the semi-discretized form of the energy equation: where the index o refers to the old time step, m to the known mass flux, k to the old iteration step inside the T-h iteration and k + 1 to the new iteration step. The star * indicates an intermediate value.
The new iteration step of the enthalpy h k+1 is approximated via a Taylor expansion: and inserted into the semi-discretized energy equation: Equation (17) is linear in T k+1, * and can be solved by a linear matrix solver. It is important to note that the second term on the left hand side is zero when the solution converges, since in this case T k = T k+1, * . The new iteration step of the enthalpy is calculated by inserting T k+1, * into Equation (16). However, it is possible that cells skip the sharp enthalpy change with this procedure. To overcome this problem, the temperature is forced back onto the T-h curve by Equation (10) once the new enthalpy value is calculated:

Optimum Density Approach
If the density is constant, the optimum approach leads to a minimized number of iterations. However, for most PCM, the density is discontinuous at the solid-liquid phase change. Moreover, for the outflow boundary to work, the total enthalpy has to be convected. The optimum density approach addresses these two issues. Here, the starting point is the semi-discretized energy equation with variable physical properties and total enthalpy transport: A Taylor expansion of the product consisting of density and enthalpy accounts for the discontinuous density: The term involving both derivatives approaches zero quadratically and is therefore not taken into account. With the enthalpy in the phase change region defined as the derivative of h with respect to T is while the derivative of the density with respect to the temperature follows from Equations (11) and (14): Second, we Taylor expand the convection term: Here, only the enthalpy is expanded to be consistent with the PISO algorithm, which gives the conservative mass fluxes, i.e., the product of the density and the velocity at the cell faces ρ m u m , at the iteration level m. Inserting Equations (20) and (24) into the semi-discretized energy Equation (19) yields: The second term and the summation of the fifth and sixth terms are zero at convergence. After Equation (25) is solved by a linear matrix solver, the temperature is forced back onto the T-h curve, the new phase fraction value is calculated by Equation (13) and the material properties are updated. The iteration process is ended if the residuum of α is reached, i.e., max(|α k+1 − α k |) < 10 −6 .

Source Based Method
The source based method is based on the idea of splitting the enthalpy into a sensible and a latent part and treating the latent part as a known source term. This yields: This equation is linear in T k+1 and can be solved readily. Next, the phase fraction α is updated under the assumption that the enthalpy inside a cell is the same before and after the update step where ω is an under relaxation factor and To assure boundedness of the phase fraction, the calculated α k+1, * has to be corrected: In a final step, the material properties are updated, the residuum is calculated and a new iteration step is started if the residuum is not fulfilled.

One-Dimensional Solidification
The first test case was a one-dimensional solidification Stefan problem with a time dependent convective boundary condition and a density change. The particular nature of this problem is its analytical solution found by Tarzia [24]. For this test case, we did without explicit units, but assumed that they are consistent.
The one-dimensional simulation domain is depicted in Figure 2. The infinite analytical domain was cut off at a length of l = 10, where a fixed pressure outlet boundary condition is applied. Initially, the PCM was liquid with a temperature of T ini = 1. The material properties did not differ between the two phases, except for the density: λ = c p = ρ l = 1, ρ s = [0.9, 0.8, 0.7]. Both the Darcy constant and the viscosity were set to zero: D = η = 0. The freezing temperature T f of the material was 0, which was approximated by solidus and liquidus temperatures of T S = −0.01 and T L = 0.01, respectively. At t > 0, a convective cooling condition with a time-dependent heat transfer coefficient was applied at the left boundary: ∂T ∂x with Subsequently, the material started to solidify and the increase in volume pushed the liquid through the outlet at the right. The domain was divided into 800 cells (∆x = 0.0125) and the time step was 0.001. The simulation ended at a time of 8.
The analytical solution of this problem is rather long and is not repeated here. It can be found in Tarzia [24].

Two-Dimensional Melting with Convection
The geometry of the second test case was a cubic cavity with an edge length of l = 40 mm and a velocity outlet at the top left corner (Figure 3). The outlet velocity was obtained by the incompressibility assumption, i.e., ρ = ρ(T). After each solution of the nonlinear T-h-coupling, the mass inside the capsule was calculated via a sum over all cells m = ∑ ρ m V and subtracted from the mass at the beginning of the time step. This mass difference had to leave the capsule in this time step, as the density was not affected by the pressure: Rearranged, this yielded a formula for the outlet velocity: The correct outlet velocity was calculated iteratively with the outer iteration loop. The pressure at the outlet was allowed to find its own value: This outflow boundary was necessary, because a zero gradient velocity boundary condition assumes a fully developed flow. This assumption does not hold if the outflow boundary is simply prescribed at the edge of the cavity. In that case, a recirculation area develops. Another solution would be a geometry extension at the outlet, e.g., a short pipe.
At the beginning, the PCM inside the cavity was completely solid at a temperature of T ini = 298.15 K. With the start of the simulation, the temperature of the left cavity wall was raised to T h = 308.15 K. The temperature of the right cavity wall was kept at the initial temperature T c = T ini . At the two remaining cavity walls and at the outlet boundary, the normal gradient of the temperature was set to zero:  Table 1; solid density and viscosity were measured in-house. The spatial discretization of the convection term has to be chosen carefully, due to the steep gradient of most variables at the phase boundary. Generally, first-order schemes are far too diffusive, while second-order schemes cannot predict the phase boundary accurately enough. To overcome this problem, the convection term was discretized by a high resolution scheme. For the two Taylor expanded convection terms appearing in Equation (24), a simple upwind discretization sufficed, as they canceled each other out at convergence. The spatial discretization resulted in both symmetric and asymmetric matrix systems. The former were solved by a preconditioned conjugate gradient solver with a diagonal incomplete Cholesky preconditioner method, whereas the latter by a preconditioned biconjugate gradient solver with a diagonal incomplete lower upper preconditioner. The residual for the phase fraction was set to 10 −6 in 1D and 10 −7 in 2D.

Validation Case
For the validation of the numerical model, we used PIV and phase front measurements by Faden et al. [21]. Their experimental setup consists of a quadratic, transparent polycarbonate pipe, with a length of 40 mm and a wall thickness of l poly = 1 mm, which is clamped between a copper and an aluminum plate. The PCM is inside the pipe and initially at T ini = 298.15 K. With the start of the experiment, the temperature of the copper plate is raised to T h = 308.15 K, whereas the aluminum plate is kept at the initial temperature. During the experiment, the setup is insulated by foamed polystyrene, whereby the insulation on top of the pipe stands out a bit. Excess material can leave the capsule through a passage way in the copper plate.
We mimicked the experimental setup and extended the simulation domain to include the top and bottom capsule walls and the insulation (Figure 4), but reduced the problem to two dimensions. The surrounding air and the cork plate on which the experimental setup was mounted were modeled as a thermal resistance. The extension of the simulation domain was necessary because the thermal conductivity of the polycarbonate capsule walls was higher than that of liquid octadecane. This led to a parasitic heat flux, which altered the shape of the phase front, especially at the bottom of the capsule, where the convective motion of the liquid PCM hampered melting. Moreover, the insulation was clamped between the copper and the aluminum plate, which led to an additional parasitic heat flux.  The solid capsule walls were treated computationally, similar to the solid PCM. A Darcy source term suppressed the velocity, and the thermal properties were implemented via two additional variables ζ an θ. The heat transfer coefficient between the experimental setup and the surrounding air was calculated by using the VDI Heat Atlas [27]. This yielded a heat transfer coefficient of k 1 ≈ 7 W/m 2 K at the top and k 2 ≈ 3 W/m 2 K at the bottom. The ambient temperature was set to T ∞ = 295.15 K. The time step was determined by a maximal Co-number of 1.0 with a maximum time step of 1.0 s. The solver settings and the other boundary conditions were the same as in the previous section. The thermophysical properties of polycarbonate and the polystyrene insulation are summarized in Table 2.

Mesh Influence for the Convection Test Case
Three meshes with 10,000, 40,000 and 90,000 cells were created and compared. The relative differences in the global liquid fraction with respect to the fine mesh were highest at the beginning of the simulation with 60% for the coarse and 20% for the medium mesh. This huge difference can be explained easily: With an increasing number of cells, the cell size decreased and therefore the point in time where the first cell was liquid was earlier. Throughout the rest of the simulation, the maximum relative difference between the coarse and the fine mesh was approximately 6% and 1% for the medium mesh ( Figure 5). A relative difference of 1% was acceptable, thus we chose the medium mesh for all further computations, including the validation case. Here, the total number of cells was higher (58,000), because of the extended simulation domain, which included the insulation and the polycarbonate pipe walls.

One-Dimensional Solidification
For the one-dimensional solidification test case, both approaches yielded almost the exact same results, i.e., they did not differ up to the sixth decimal place in the global liquid fraction at the end of the simulation. For that reason, the figures displaying the liquid fraction or the temperature only show one solution method. The temperature distribution obtained by the numerical solution closely matches the analytical solution and it is hard to distinguish the two lines by eye ( Figure 6). Hence, both methods give correct results and can be employed for solving solid-liquid phase change accompanied by a density change. However, both methods differed by the required number of nonlinear iterations (Figure 7). The source base method needed approximately nine times as many iterations compared to the optimum density approach, slightly depending on the density ratio.
For both approaches, the number of nonlinear iterations increased with decreasing density ratio, but the increase was moderate.  The density ratio not only has a computational influence but also a physical one. The liquid fraction at the end of the simulation rose with the density ratio, since the volume expansion upon solidification pushed the liquid through the outlet and therefore reduced the amount of liquid inside the simulation domain. However, the volume change also resulted in a higher thermal resistance. Therefore, the functional relation between the density ratio and the liquid fraction was somewhat less than linear (Figure 8). A phenomenon inherent to all enthalpy methods during diffusive dominated phase change is a stepwise progression of the phase boundary. This results in an oscillating velocity for density varying phase change (Figure 9). These oscillations have been observed before, but it seems that they were stronger here than those found by Alexiades and Drake [22]. This can be explained by the effect that segregated pressure velocity algorithms generally tend to produce oscillations at discontinuous surfaces. However, the influence of these oscillations on the temperature and phase fraction is small, as in the finite volume method they are convected by the mass flux on the cell faces and these mass fluxes do not oscillate. If two-or three-dimensional problems are studied and natural convection is taken into account, the progression of the phase front is no longer stepwise and the velocity oscillations vanish.

Two-Dimensional Melting of Octadecane with Natural Convection
As previously, the optimum density approach required less computational effort than the source based method, but the source based method did not perform as poorly as in the one-dimensional case. This was also found by Voller and Swaminathan for a comparison without density change in the original optimum approach paper [14]. Depending on the Courant number, the source based method required 4-7 times as many iterations compared to the optimum density approach (Figure 10).
It was obvious that the maximum Courant number and thus the time step greatly influenced the iterations required. A small time step meant more time steps and therefore the T-h coupling had to be solved more often, but fewer iterations per time step were necessary. However, two difficulties arise, if the time step is increased: First, the solution methods for the T-h coupling can become unstable. This happened here for a Courant number greater than 2 for the source based method. This may be resolved by applying underrelaxation. Second, the number of outer iterations per time step has to be increased, because the coupling between the momentum and the energy equations is not solved accurately enough. This happened here for the optimum density approach at max Co > 4. Here, the solution would be an increased number of outer iterations. However, in both cases, the computational effort was increased despite a higher time step. That is why we did not increase the time step any further.
The difference in the iteration number affected the required computational time strongly. For a maximum Courant number of 2, the source based method needed four times as long as the optimum density approach with the same maximum Courant number. When the stability of the solver was also taken into account, the source based method with a Courant number of 2 needed approximately 5.5 times as long as the optimum density approach with a Courant number of 4. This comparison seemed reasonable, since the influence of the time step on the liquid fraction after 4 h was small, e.g., for the optimum density approach, the difference in the global liquid fraction between the largest and the smallest time step/maximum Co-number combination was below 0.02 %.  At the end of the simulation, the source based method required over five times as many nonlinear T-h iterations compared to the optimum density approach and six times as many calls to the linear matrix solver. The convergence history inside a time step emphasizes this behavior ( Figure 12). Clearly, the residual of the optimum density approach decreased more quickly. The position of the phase front is shown in Figure 13 (left) for the optimum density approach. Showing only the solution obtained by the optimum density approach suffices, because the liquid fraction and the position of the phase front did not differ much compared to the source based method. The optimum density approach gave a global liquid fraction at the end of the simulation of α = 0.643, and α = 0.641 with the source based approach. The reason for this small difference was the different discretization of the convection term. This is an issue the authors plan to investigate further in a following paper.

Validation
All simulations presented in this section were performed with the optimum density solver on the extended simulation domain. The measured [21] and the numerically obtained phase front positions agree well ( Figure 13).
After 60 min, the influence of natural convection was clearly visible as the phase front was tilted to the right and melting at the top of the capsule was clearly faster than at the bottom. With increasing time, the influence of convection became even more dominant. After 4 h, only a small streak of solid PCM separated the right capsule wall and the liquid PCM at the top of the capsule, while at the bottom of the capsule the liquid PCM progressed to about a quarter of the capsule length for the experimental results and a bit less than a quarter for the numerical ones. The kink in the phase front shape at the bottom of the capsule, resulting from the parasitic heat flux through the capsule walls, was also present in the numerical solution on the extended simulation domain. The greatest deviations arose at the top of the capsule, where the numerical solution predicted a larger already melted area. Possible reasons for this deviation are: increased natural convection due to the two-dimensionality assumption, uncertainties in thermo-physical properties, heat loss while the top insulation was removed for taking pictures required for PIV and the presence of an air bubble in the experiment, which hampered melting at the top of the capsule. This air bubble was visible in the experimental velocity field in the top left corner, as an area with zero velocity (Figure 14). Similar to the phase front positions, the velocity fields matched well. Of course, the measured velocity field was not as smooth as the numerically obtained one. Basically, the flow was a large vortex driven by buoyancy forces, whereas the highest velocities were at the left cavity wall and at the phase boundary. The fluid inside this vortex was at rest. In both cases, the boundary layer thickness towards the solid boundaries was considerably thinner than towards the fluid at rest. The temperature distribution of the extended simulation domain (Figure 15), including the PCM, the walls and the insulation, confirmed the above-mentioned parasitic heat fluxes. At the bottom of the capsule, the wall and the complete insulation were clamped between the hot copper and the cold aluminum plate. For this reason, the temperature inside the insulation and the wall was higher than that of the PCM at the same horizontal position and a parasitic heat flux flowed in the PCM. Thereby, the influence of the wall was higher, as the thermal conductivity was higher. Nevertheless, the kink in the phase boundary could not be explained by including only the wall and applying an adiabatic boundary condition. Although the insulation at the top stood out slightly, its general influence on the wall and the insulation was the same.
Transparent materials with a distinctly lower thermal conductivity than polycarbonate are hard to find. Thus, the only way of decreasing the parasitic heat flux through the wall is to make them thinner. The influence of the insulation, however, could be removed if the experiments were conducted in an environment with controlled temperature. This would also simplify the necessary boundary condition for a validation simulation. Figure 15. Temperature profile in the PCM, the pipe walls and the insulation after: 1 h (left); and 3 h (right). The black lines were inserted to guide to eye and mark the pipe walls. The melting temperature is white.

Conclusions
Taylor expanding the density and the enthalpy around the melting temperature can significantly reduce the computational effort needed to solve the non-linear energy equation in solid-liquid phase change with volume expansion or shrinkage. The algorithm obtained by following this approach needs a minimized number of iterations and gives the same results as the slowly converging source based method. Moreover, the number of iterations increases only slightly with an increasing density ratio.
The new algorithm was validated against experiments and our numerical results show that parasitic heat fluxes can alter the shape of the phase front near the bottom wall significantly. It is therefore necessary to include certain parts of the experimental setup, e.g. the capsule walls in the simulation domain.
We anticipate that volume expansion upon melting or solidification will be included in more and more simulations regarding latent thermal heat storage, since density changes of 10% or more, which are typical for PCM, cannot simply be ignored in many applications, e.g. macro-capsules in which the volume expansion leads to stress in the capsule walls. Even more importantly the energy conservation is violated.
Interestingly, the influence of this violation and of the validity of the Boussinesq approximation for melting process is not known, although it is widespread. This is an issue we will address with this new algorithm in a forthcoming paper.
By solving this equation, the pressure field is adjusted in a way that the mass fluxes at the cell faces obey continuity. This is not equivalent to a divergence free velocity field at the cell centers.

Appendix A.2. PISO-Algorithm
The pressure equation derived above is fundamental for the PISO-Algorithm [29]. Hereafter, the index m corresponds to either the old time step or the previous outer iteration and the index l to the PISO-loop itself. The PISO-algorithm contains the following steps (at first iteration l = m): 1. Momentum predictor: C m u m, * = r m − ∇p m − ( g · x∇ρ) m (A11) 2. Update the matrix coefficients: 3. Solve the pressure equation: H( u l )) = r m − (H ) m u l (A15) 4. Calculate the conservative mass fluxes: 5. Velocity corrector: 6. Go back to Step 3, if the residuum or the maximum number of iterations has not yet been reached.
The velocity corrector does not involve the solution of a system of linear equations, because the left and right side are evaluated at different iteration steps (operator splitting). This means that the velocity is only corrected by the pressure term −(A −1 ) m ∇p (l+1) . To take the influence of the velocity term (A −1 ) m H( u l ) into account, the PISO-loop has to be repeated twice [30]. The buoyancy term (A −1 g · x∇ρ) m is updated through outer iterations. Moreover, the coefficients of H and A −1 are not updated inside the correction loop, although they depend on the mass flux of the last time step φ. This is done, because it is believed that the pressure-velocity coupling is more important then the nonlinearity of the convection term [31]. However, the coefficients are of course updated within the outer iterations.