1 Introduction
Melting and solidification processes in phasechange materials are relevant to many areas of engineering and scientific research. Examples include the casting of metals, storage of solar and thermal energy as latent heat, and modeling phasechange processes in the cryosphere, e.g. as relevant for climate projections. In some regimes, the phasechange process is substantially affected by convection in the liquid phase, which has been demonstrated in a number of experiments Sparrow et al. (1978, 1979); Okada (1984); Kowalewski and Rebow (1999); Schüller et al. (2017).
Any spatiotemporally resolved simulation of this multiphysics process is a mathematical and computational challenge. Common modeling strategies couple systems of multiparameter nonlinear partial differential equations resulting from the balance laws of mass, momentum, and energy. For realistic scenarios, solving the system requires efficient and robust numerical methods. The design and validation of models for different regimes and the formulation of accurate solution methods are areas of ongoing research.
The current work focuses on pure materials at macroscopic scale, for which melting and freezing occur at one specific temperature, namely the melting/freezing point . The corresponding phasechange process is commonly referred to as isothermal phasechange Voller et al. (1987). Isothermal phasechange produces well defined, distinct phase interfaces at small spatial scales. There are many other materials that do not fall into this class, e.g. alloys or sea ice, and that produce substantial mushy regions. For such other materials, melting and freezing processes are bounded by solidus and liquidus temperatures which can both vary depending on the local material composition. Isothermal phasechange can hence be interpreted as the limiting case, in which the solidus and liquidus temperatures coincide and are constant.
Many computational fluid dynamics (CFD) applications are concerned with a physical setting that consists of entirely gas or entirely liquid materials, sometimes extending into multiphase gasliquid applications. In these cases, solids exist only as boundaries and are usually fixed in space and time. Simulating convectioncoupled melting and freezing requires the extension of CFD into multiphase liquidsolid applications. This involves modeling the spatiotemporally evolving interfaces between liquids and solids. For this, one can use either an interface tracking or an interface capturing method Voller et al. (1987); Alexiades and Solomon (1992). Interface tracking methods explicitly track the phase interfaces, solve separate systems of governing equations in each phase domain, and enforce coupling constraints at the interfaces. Interface capturing methods instead solve one system of governing equations on a single domain that is occupied by both phases. Interface capturing handles phase interfaces implicitly. The phase interfaces can be postprocessed as an explicit function of the solution.
The current work uses an interface capturing method that relies on an enthalpy formulation of the phasechange process. This is often also referred to as an enthalpy method Voller et al. (1987). Enthalpy methods write the energy balance in terms of both the temperature and the enthalpy, and hence need to be closed with an equation relating these two quantities. For incompressible materials, this relation is wellunderstood and can be easily phrased in terms of the liquid and solid phase volume fractions. Enthalpy methods can hence be closed by assuming an equation of state relating the phase volume fractions to to the material’s effective temperature. Considering isothermal phasechange, the equation of state reduces to an indicator type functional relation. This choice results in an energy balance taking the form of a standard heat equation extended by a source term which accounts for the gain or loss of latent heat at the phase interface. Enthalpy methods have been applied extensively to convectioncoupled phasechange Voller et al. (1987); Voller and Prakash (1987); Brent et al. (1988); Alexiades and Solomon (1992); Giangi et al. (2000); Evans and Knoll (2007); Belhamadia et al. (2012); Danaila et al. (2014); Zimmerman and Kowalski (2017); Rakotondrandisa et al. (2019), but there remain many opportunities for improvement.
Enthalpy methods are often applied to regimes where the velocity vanishes in solid regions. Since interfaces between phases are not tracked explicitly, typical velocity boundary conditions, such as the commonly used noslip condition, cannot be applied directly. Rather, different modeling techniques are used to guarantee a vanishing velocity in the solid phase. These are sometimes referred to as a solid velocity correction Wang et al. (2010). Different approaches are compared in Voller et al. (1987, 1990). Three primary approaches to solid velocity correction have prevailed:

modifying the algorithm to directly set the velocity to zero in the solid regions Wang et al. (2010),
The first approach provides a pragmatic brute force solution and allows no velocity dampening close to the phase interface. The second and third approaches both rely on additional numerical parameters that can affect the evolution of phase interfaces. Critically, these new parameters, which initially do not have any physical significance, must be calibrated based on experimental data, as seen in Rakotondrandisa et al. (2019) or Danaila et al. (2014). This need for calibration may undermine the predictiveness of the models, because it is not clear to what extent the calibrated parameters also apply to other physical regimes or geometric settings. Furthermore, recent applications of the second approach to isothermal convectioncoupled phasechange, e.g. Belhamadia et al. (2012); Schüller et al. (2017); Rakotondrandisa et al. (2019), have used the KozenyCarman relation in the momentum source term. This contrasts with the explaination in Voller et al. (1987); Brent et al. (1988) that the KozenyCarman relation is most physically relevant for nonisothermal mushy region phasechange. A simpler formulation of the velocity relaxation can be assumed for isothermal phasechange of pure materials Brent et al. (1988); Evans and Knoll (2007) and will be leveraged in the current article.
For any chosen approach to solid velocity correction, the step change in velocity at phase interfaces disrupts nonlinear solver convergence Brent et al. (1988). A common regularization approach is to introduce a temperature range in which the phasechange occurs rather than restricting it to exactly the melting/freezing point. There seems to be a strong analogy between this type of numerical method and the physical situation of nonisothermal phasechange for nonpure materials. In both situations, the community tends to refer to the lower and upper temperature bounds respectively as solidus and liquidus temperatures. It is, however, important to distinguish cases in which solidus and liquidus temperatures are physically motivated versus cases in which they are introduced only to stabilize an otherwise unstable computational model. Simulations in the current work consider sensitivity to the phase interface regularization, and reduce the regularization until it does not create unphysical mushy layers or substantially affect the shape of phase interfaces.
Enthalpy methods define a system of partial differential equations which must be discretized in space and time. Many approaches have been applied for spatial discretization, including finite differences Voller et al. (1987); Voller and Prakash (1987); Brent et al. (1988), finite volumes Giangi et al. (2000); Evans and Knoll (2007); Wang et al. (2010); Schüller et al. (2017), the finite element method (FEM) with operator splitting Belhamadia et al. (2012), and monolithic mixed FEM Danaila et al. (2014); Zimmerman and Kowalski (2017); Rakotondrandisa et al. (2019); Woodfield et al. (2019). The formulations in Danaila et al. (2014); Zimmerman and Kowalski (2017); Rakotondrandisa et al. (2019); Woodfield et al. (2019)
all used TaylorHood elements simultaneously with the pressure penalty method to stabilize the mixed finite element system. TaylorHood elements increase the degree of the velocity basis, adding a degree of freedom per cell in each dimension. The current work demonstrates successful stabilization using the pressure penalty method alone. In the context of enthalpy formulated convectioncoupled phasechange, the stability of mixed FEM schemes was analyzed for the first time in
Woodfield et al. (2019). Also in this context, per the current authors’ knowledge, the only formal accuracy verifications of implemented solvers are in Woodfield et al. (2019) and in the current work. Both use the method of manufactured solutions to verify the orders of accuracy of the spatial and temporal discretizations. Typically, temporal discretization is performed with finite differences. Secondorder accurate fully implicit finite difference temporal discretization was advocated in Evans and Knoll (2007) and subsequently used in Belhamadia et al. (2012); Rakotondrandisa et al. (2019).Based on the approach in Danaila et al. (2014), an open source implementation was published in Zimmerman and Kowalski (2017). Both used mixed finite elements in space and the backward Euler method in time. While the results in both Danaila et al. (2014) and Zimmerman and Kowalski (2017) were promising, the nonlinear solvers were not robust. Convergence required carefully chosen time step sizes, dynamic adaptive mesh refinement (AMR), and artificial initialization of the new phase in oneway melting or solidification processes. AMR requires an a posteriori
error estimate, which in turn requires first obtaining a solution on a coarse or manually refined mesh
Zimmerman and Kowalski (2017). Furthermore, it has not been shown that the benefits of AMR outweigh the costs. The benefit of reduced degrees of freedom competes with the cost of dynamic remeshing. The phase interfaces are moving in time, and the locations of these interfaces are where the mesh must be locally refined. The possible benefit of perelement accuracy from locally refined meshes competes with interpolation errors between nonmatching meshes of different time steps. At least for twodimensional geometries, more accurate solutions might be obtained with less computational cost on static uniform meshes.
This article is structured as follows: Section 2 presents an enthalpy formulation for convectioncoupled phasechange processes with a generalized approach to phaseinterface regularization. Section 3 presents a monolithic solver approach with mixed finite elements, implicit time discretization, and Newton’s method. Motivated by Brent et al. (1988); Zimmerman and Kowalski (2017), a robust continuation procedure for solving the nonlinear problem is introduced. Section 4 introduces an opensource implementation and describes the opensource libraries on which it depends. The solver’s second order accuracy is verified both in space and time. Section 5 compares simulation results both for the melting of octadecane and for the freezing of water to experimental data from the literature.
2 Physical and mathematical model
2.1 Physical regime
The current work focuses on simulating melting and solidification processes that are strongly influenced by convection of the liquid phase. Particularly, the focus is on regimes where there is a strong twoway coupling between the fluid dynamics and the evolving shapes of phase interfaces. The liquid is assumed as an incompressible Newtonian fluid undergoing natural convection. Following a Boussinesq assumption, the fluid density is assumed to only depend on the temperature. In this case, the flow is driven by buoyancy resulting from temperature gradients. Furthermore, the focus is on pure materials such as octadecane and distilled water. Such materials undergo isothermal phasechange Voller et al. (1987); Brent et al. (1988). This means that melting and solidification occur at a single temperature.
2.2 Phase interface regularization
Consider a phasechange material which can either be solid, liquid, or a mixture of the two phases. Any small reference volume of this material has a liquid volume fraction and a solid volume fraction . Assuming there are no other phases, e.g. no gas inclusions, . When the material is fully liquid, . When it is fully solid, . In this case, the value of can be used as a substitute for the phase state.
Isothermal phasechange processes are characterized by having distinct interfaces between phases. At these interfaces, the temperature of the material is at the melting/freezing point . In liquid regions, . In solid regions, . This means that the liquid volume fraction can be reconstituted from the temperature alone, i.e.
(1) 
With this model, phase interfaces are theoretically infinitesimally thin lines in twodimensional space or surfaces in threedimensional space. The phase interfaces can be regularized by convoluting with a Gaussian kernel, e.g.
(2) 
yielding
(3) 
This introduces a parameter, , which is always positive. As approaches zero, approaches . Alternative regularizations of phase interfaces that rely on similar parameters have been proposed in other literature Danaila et al. (2014); Zimmerman and Kowalski (2017); Rakotondrandisa et al. (2019). With too large of , substantial regions may develop where . Such regions are numerical artifacts. The parameter exists solely to make solutions computable. Therefore, for results in Section 5, was reduced until regions where were small and until the shapes of the contours were not affected.
2.3 Governing equations
Convectioncoupled phasechange can be simulated in enthalpy form Voller et al. (1987) using the following governing equations, which are balances of mass, momentum, and energy
(4)  
(5)  
(6) 
The independent unknowns are pressure , velocity , and temperature . The mass equation (4) is the continuity equation for an incompressible material. The momentum equation (5) begins with the incompressible NavierStokes momentum equation with two additional terms. The first additional term, , is a Boussinesq buoyancy term with temperaturedependent . Under the Boussinesq approximation, density variations are only considered with respect to the temperature. The incompressible flow considers a reference density
. The symmetric part of the rateofstrain tensor is
.The second additional term to the momentum equation, , with , forces the velocity toward zero in the solid region. This is sometimes referred to as the solid velocity correction Wang et al. (2010). The regularized solid volume fraction is , with . This solid velocity correction term is similar to that proposed in Brent et al. (1988) for isothermal phasechange. In other recent work applied to isothermal phasechange Belhamadia et al. (2012); Danaila et al. (2014); Schüller et al. (2017); Rakotondrandisa et al. (2019), the KozenyCarman relation was instead used. In addition to being a nonlinear function of , the KozenyCarman relation adds a new parameter which is needed to avoid division by zero. The KozenyCarman relation was designed to account for the effects of permeability in porous media. For this reason, it is sometimes used in enthalpy methods applied to nonisothermal phasechange processes that involve substantial mushy layers. Given no physical justification for applying the KozenyCarman relation to isothermal phasechange, the simpler term is instead used in the current work. The same choice was made in Evans and Knoll (2007), where a pure material was also simulated.
The energy balance is in enthalpy form Voller et al. (1987). This invokes a reference temperature and produces a source term accounting for the phasechange with latent heat . The volumetric heat capacity and thermal conductivity can both depend on the phase. The volumeaveraged values are and .
2.4 Nondimensionalization
Based on the approach in Danaila et al. (2014), the independent variables were nondimensionalized with respect to characteristic scales. These include spatial scale , speed scale , and temperature scale . The liquid material properties were used as reference values, e.g. the liquid kinematic viscosity . Additionally, the reference temperature was chosen as the isothermal melting and freezing temperature, i.e. . The nondimensional variables were chosen to be
(7) 
In this case, the nondimensional melting and freezing temperature is always . The phasedependent material properties were normalized with respect to the liquid values, i.e.
(8)  
(9) 
In nondimensional form, the governing equations (4), (5), and (6) then read
(10)  
(11)  
(12) 
where the Grashof, Stefan, and Prandtl numbers are defined as
(13) 
In the solid velocity correction term, is the nondimensional velocity relaxation factor. In the Boussinesq buoyancy term, represents a general temperaturedependent density model which must be specified for a given liquid. In cases where this density model is nonlinear, e.g. when accounting for the density anomaly observed in water, a constant reference value must be chosen for the thermal expansion coefficient. The liquid thermal diffusivity is .
When referring to the governing equations in the remainder of this article, only the nondimensional system (10, 11, 12) is considered. The scaling highlights the regimes of natural convection and phasechange. For example, Section 5 presents simulations of laboratory scale experiments for the melting of octadecane and freezing of water, and describes the cases in terms of their similarity parameters. The code presented in Section 4 solves the nondimensional equations.
3 Numerical methods
The timedependent problem was solved as a sequence of initial boundary values problems. Mixed finite elements were used for the spatial discretization. Backward difference formulas were used for the temporal discretization. For the current work, nonhomogeneous Dirichlet and homogeneous Neumann boundary conditions were applied to the temperature field, homogeneous Dirichlet boundary conditions were applied to the velocity field. The nonlinear system was solved with Newton’s method. Robustly and accurately solving the nonlinear problem required special attention to its regularization, for which a continuation procedure was developed.
3.1 Spatial discretization
Following a similar approach to Danaila et al. (2014); Zimmerman and Kowalski (2017), mixed finite elements were used to approximate the system (10, 11, 12) with basis functions , yielding the weak form residual
(14)  
where and .
Piecewise linear polynomials were used for all basis functions. Since these do not satisfy the LBB compatibility condition for the incompressible NavierStokes system Donea and Huerta (2003), the system was stabilized via the pressure penalty method with . This also ensures the uniqueness of the pressure solution Woodfield et al. (2019). For similar problems, a combination of both the pressure penalty method and TaylorHood elements had been used in Danaila et al. (2014); Zimmerman and Kowalski (2017); Woodfield et al. (2019). This means that, where a piecewise linear basis approximates the velocity field in the current work, Danaila et al. (2014); Zimmerman and Kowalski (2017) used piecewise quadratic approximations without improving the order of accuracy. With cells of a mesh in dimensional space, the equalorder basis has less degrees of freedom than the TaylorHood basis. In Section 4, this spatial discretization’s second order accuracy is verified.
3.2 Temporal discretization
The time derivatives were discretized with the family of order constant time step size () backward difference formulas (BDF) Ascher and Petzold (1998). For some unknown at time step , these can be written as
(15) 
where the coefficients , shown in Table 1 depend only on . This requires storing solutions from the previous time steps. For example, refers to the already computed solution from the previous time step.
1  1  1  

2  3/2  2  1/2  
3  11/6  3  3/2  1/3 
3.3 Nonlinear solution with continuation
Many terms in the residual (14
) are highly nonlinear. The nonlinear system was solved with Newton’s method as follows. Consider the vectorvalued system solution
for any discrete time step . The nonlinear problem’s solution is approximated for the next time step by solving a sequence of linear problems(16) 
whose iterates converge to . Here, stands for the Gateaux derivative. It can be derived analytically, which was done for a similar problem in Zimmerman and Kowalski (2017). Or as done in the current work, it can be automatically computed with the symbolic capabilities of the software in Section 4. Each iteration of (16) is a linear system which can be solved robustly with a direct solver.
When trying to solve the nonlinear problem for the solution at time step , it was first attempted to use the known solution from the previous time step as an initial guess. With small regularization parameter in (3), it turned out that was not a suitable initial guess. In Danaila et al. (2014); Zimmerman and Kowalski (2017), this issue was resolved via adaptive mesh refinement, careful selection of the time step size, and suitably tailored initialization of a melt layer in melting problems or a frozen layer in freezing problems. In the current work, the nonlinear problem was solved robustly with a continuation procedure. Initial guesses were improved by solving intermediate problems with an increased value of the regularization parameter . Whenever a solution was found, it was used as an initial guess for solving the same system again, albeit with a smaller value of . The challenge hence is to find a sequence of values resulting in converging intermediate steps that reach the originally specified value. Figure 1 explains the procedure in more detail.
The resulting sequence of intermediate regularizations depends on the exact nonlinear solver method used, and on its options such as convergence criteria and maximum iterations. It was found that using the relative residual as a convergence criteria would allow for significant accumulation of error in the absolute residual between time steps. For this reason, the relative tolerance convergence criteria was disabled. The nonlinear solver was limited to twentyfour maximum iterations to reach the convergence criteria. Section 5 reports the total number of Newton iterations for each of the selected simulations. Despite having to solve a sequence of nonlinear problems for each time step, the total numbers of Newton iterations were not excessively high, because the successfully solved intermediate problems converge optimally.
4 Implementation and verification with opensource libraries
In order to reduce the required coding effort, leverage stateoftheart methods and solvers, and increase the usability of the software and reproducibility of results, the code was implemented using the open source finite element library Firedrake Rathgeber et al. (2016). On top of Firedrake, for this work, the open source Sapphire Zimmerman (2019) framework was developed. Its purpose is to construct time dependent simulations using the methods from Section 3. The convectioncoupled phasechange simulations used to obtain the results in Section 5 were implemented within the Sapphire framework.
4.1 Firedrake
In recent years, a technology has been developed to automatically implement finite element discretizations of variational problems, the Unified Form Language (UFL) Alnæs et al. (2012). UFL provides a domain specific language for users to abstractly define their problem, and automatically writes and compiles efficient C++ code. From the existing libraries based on UFL, this work used Firedrake Rathgeber et al. (2016). In addition to finite element discretization, Firedrake was used to generate meshes of simple geometries, to interface with efficient linear and nonlinear solvers from PETSc Balay et al. (2018), and to output solutions to a standard format. Firedrake is designed for performance portability Balay et al. (1997); Dalcin et al. (2011). For the results presented in Section 5, all linear systems were solved with MUMPS Amestoy et al. (2001, 2006), and all nonlinear systems were solved with PETSc’s SNES (Scalable Nonlinear Equations Solvers) line search method, both via the Firedrake interfaces.
By default, Firedrake automatically determines the quadrature degree for a given problem. Alternatively, the user can specify a degree. In the current work, when integrating the regularized liquid volume fraction (3) with small regularization parameter value , Firedrake automatically determined excessively large quadrature degrees. Computations were performed more efficiently by reducing the quadrature degree, which reduced the number of function evaluations needed for each degree of freedom when performing finite element assembly. This essentially introduced a new numerical parameter to the model, the quadrature degree, . Therefore, is included in the sensitivity studies of Section 5.
Firedrake depends on a stack of open source scientific software libraries. A benefit is that ongoing technological developments from a variety of open source scientific software projects are leveraged. A downside is that the software can be difficult to configure and compile. To manage this complexity, the Firedrake project maintains an installation script. Furthermore, to facilitate the reproduction of scientific results, Firedrake provides a command line tool which installs a specified set of versions and configurations throughout the software stack. The exact versions of Firedrake and its software stack used for the results in the current work were documented with a DOI on Zenodo zenodo/Firedrake20190620.0 (2019), ensuring reproducibility.
4.2 Sapphire
Sapphire was developed for the purposes of the current work. Sapphire provides a Python class for time dependent simulations that are governed by PDE’s. The PDE’s must be discretized in space with finite elements and in time with finite differences. Spatial discretization is handled by Firedrake, with the PDE being defined in variational form using UFL. The time discretization module includes backward different formulas up to sixth order. Five arguments are required to instantiate a Sapphire simulation:

a mesh, e.g. from Firedrake’s builtin mesh functions or converted from an external meshing library,

a finite element or mixed finite element defined by UFL,

the residual of the governing equations in variational form, e.g. (14), defined by UFL,

a (possibly empty) list of Firedrake Dirichlet boundary conditions,

initial values as a Firedrake finite element function.
With these defined, a simulation is ready to run.
In Section 4.3, the convectioncoupled phasechange solver’s accuracy is verified using the method of manufactured solutions (MMS) according to the general guidelines in Roache (2002). Manually deriving MMS source terms creates many opportunities for errors, especially for systems of PDE’s involving vector calculus. Often symbolic math libraries are used to derive the source terms. Python has such a symbolic math library which does not depend on UFL. With Sapphire already depending on UFL, it was more convenient to use UFL itself for symbolically deriving the MMS source terms. An additional module in Sapphire hence uses UFL to automatically derive the source terms. The module also automatically augments the variational problem which needs to be solved for verification. As part of its test suite, Sapphire uses the MMS module to verify the accuracy of solutions to many different PDE’s.
Sapphire is a flexible tool, is publicly available, and is freely licensed. Simulations have been implemented in Sapphire which solve the heat equation, convectiondiffusion, steady and unsteady incompressible NavierStokes, natural convection, and phasechange in enthalpy form without convection. This provided a strong base for implementing convectioncoupled phasechange in enthalpy form. For the convectioncoupled phasechange simulation, the most significant additional work was to implement the continuation procedure from Figure 1. The version of Sapphire used to produce the results in the current work is documented with a DOI on Zenodo Zimmerman (2019).
4.3 Solver verification
Numerical solutions to the weak residual (14) approximate solutions to the governing equations (10, 11, 12). As there is no sufficiently complex analytical reference solution available, spatial and temporal convergence were verified via the method of manufactured solutions (MMS). MMS is a general method for verifying PDE solvers. There is an additional benefit to MMS when approximating solutions to strong form governing equations via weak formulations. In addition to verifying the code, the weak formulation is also verified, because the weak form is tested against the strong form.
A manufactured solution can be any function which is sufficiently differentiable in space and time to exercise the spatial and temporal derivates in the governing equations. Substituting a manufactured solution, which is not an exact solution, into the governing equations yields a source term in each equation. Adding the source terms to the governing equations creates an augmented problem. To verify a solver, it is applied to the augmented problem, and the resulting solution must approximate the manufactured solution with the expected orders of accuracy. Manufactured solutions do not necessarily need to have physical interpretations Roache (2002). The purpose is to reveal coding errors and to verify the solver.
To apply the MMS procedure for verifying a solver’s spatial or temporal orders of accuracy, Sapphire requires the user to

implement a manufacture solution with UFL,

and set a single time step size and a list of meshes with decreasing cell sizes for spatial verification, or a single mesh and a list of for temporal verification.
Separate manufactured solutions can be used to independently verify the spatial and temporal discretizations.
The strong form governing equations and manufactured solutions are defined in a Python module, and the verification process is run with a Python script. Then, the module is passed to one of Sapphire’s MMS verification functions. Sapphire automatically

derives the source term for each component of the manufactured solution, multiplies them with appropriate test functions, and adds them to the weak form,

sets initial values and boundary conditions as given by the solution,

runs simulations for specified lists of or ,

tabulates errors between the approximate finite element solutions and the manufactured solutions,

reports convergence rates with respect to or ,

and asserts that the reported convergence rate matches the expected order of accuracy within a specified tolerance.
Verification of the current convectioncoupled phasechange solver was performed on a unit square domain with Dirichlet boundary conditions. The manufactured solution was
where and are the components of the position .
For the buoyancy term in (11), a classical linear Boussinesq model was used by setting . The similarity parameters were set to , , and . The material property parameters were set to , , and . The regularization parameter was set to and the solid velocity relaxation parameter was set to . The quadrature degree was automatically set by Firedrake.
To verify spatial discretization accuracy independently of any time discretization error, the manufactured solution was evaluated at time , and no further time dependency as assumed. The resulting problem was solved in a single pseudotime step. This was repeated on a series of uniformly refined meshes, corresponding to the mesh shown in Figure 3 with a different cell edge length, , for each mesh. To independently verify the temporal discretization accuracy, a sufficiently refined mesh was used such that the spatial discretization error would not dominate the total error. The unsteady problem was solved repeatedly with a sequence of time step sizes . Figure 2 shows the observed convergence orders with respect to and . Having applied second order accurate discretization methods, the effective orders of accuracy appear to be asymptotically approaching the theoretical values. The entire Sapphirebased verification procedure takes approximately a minute to run on a single CPU. It is included into a continuous integration process that runs automatically whenever the code’s master version on GitHub Zimmerman (2019) is modified.
5 Validation with benchmark experiments
Two experimental data sets from the literature were considered in order to validate the implemented model. First, octadecane melting was simulated for comparison to an experiment from Okada (1984). Second, water freezing was simulated for comparison to experiments from Kowalewski and Rebow (1999). Both simulations were performed on unit square geometric domains, and both used the same configuration of boundary conditions shown in Figure 3, on the left. The top and bottom boundaries were assumed to be adiabatic, while the left and right walls were kept respectively at constant hot temperature and cold temperature . No slip velocity boundary conditions were applied on each wall. For both benchmarks, gravity points vertically downward, i.e. . Also in both cases, a uniform triangular mesh with cell edge length was used. An exemplary mesh with is shown in Figure 3, on the right.
Each of the two following subsections is devoted to validation against one of the experiments. In both cases, the subsection includes

a description of the experimental setup.

a description of the simulation setup. This includes the scaling and similarity parameters, the phasedependent material properties, the buoyancy model, the procedure for obtaining initial values, and the boundary condition values,

a sensitivity study involving five numerical parameters. These are the mesh cell size , the time step size , the phase interface regularization parameter , the solid velocity relaxation parameter , and the quadrature degree . With respect to a nominal set of parameters values, each parameter is varied individually and the resulting phase interfaces are evaluated at a fixed time.

a selected simulation result. The selection is based on the sensitivity study. The numerical parameters for the selected case are reported. The temperature field and velocity streamlines are visualized at multiple time steps. The simulated phase interface is compared to the experiment at multiple times.
5.1 Melting octadecane
The melting of octadecane was simulated for comparison to an experiment from Okada (1984). In the experiment, a slab of octadecane was insulated on the top and bottom. The material was initially at its melting temperature. One side was heated to melt the material.
5.1.1 Simulation setup
The height of the experimental cavity as taken from Okada (1984) was used as the length scale, i.e. . A realistic value for the liquid kinematic viscosity, as taken from Wang et al. (2010), furthermore yields the characteristic time scale . Based on physical parameters provided in Danaila et al. (2014) and according to the previously introduced definitions (13), the similarity parameters are Gr = 5820, Pr = 56.2, and Ste = 0.0450. Octadecane has approximately constant volumetric heat capacity and thermal conductivity between phases. Therefore the normalized values for the energy balance (12) are . Finally, a linear Boussinesq model with constant thermal expansion coefficient is assumed. Given the defined Grashof number (13), this yields for the momentum balance (11).
For the nondimensional temperature corresponding to the melting point, i.e. , the regularized liquid volume fraction as defined in (3) is evaluated to be . This is true for any chosen regularization parameter . In order to capture the experimental conditions, in which the material at time zero had been solid with a temperature close to the melting temperature, the initial temperature for the simulation was offset by a small amount, i.e. . This approach follows earlier work, e.g. Danaila et al. (2014); Zimmerman and Kowalski (2017); Rakotondrandisa et al. (2019). For this particular choice of the initial temperature , setting yields an initial liquid volume fraction , meaning that more than 99% of the material is initially solid.
During the experiment, the left boundary was kept at a constant temperature well above the melting point, while the right boundary was kept at a temperature close to the melting point. For the simulation boundary conditions shown in Figure 3, this translates into and . The high temperature at the left boundary causes melting from left to right.
5.1.2 Sensitivity to numerical parameters
For this sensitivity study, the nominal parameter values were , , , , and . From their nominal values, each parameter was varied independently. The resulting phase interfaces were compared at time . The results are shown in Figure 4. The time is also the final simulated time presented in Figure 5.
The mesh cell size, , was varied from a maximum of to a minimum of . The upper part of the phase interface was most sensitive to . Between the smallest values, and , the position of the phase interface in the upper portion of the domain still visibly changed; but overall there appeared to be an asymptotic convergence behavior. The time step size, , was varied from a maximum of to a minimum of . For the largest size, , there was a significant decrease in the melting near the top wall. For , the result was no longer significantly changed by reducing .
The regularization parameter, , was varied from a maximum of to a minimum of . For the largest value, , the entire phase interface shifted rightward. This was because, at the initial temperature , the initial liquid volume fraction was . This meant that about ten percent of the material was already melted at the initial time, shifting the entire phase interface in the direction of melting. For , less than one percent of the material was melted at the initial time, and this effect no longer appeared to dominate. Between the three smaller values of ranging from to , there appeared to be little sensitivity.
The solid velocity relaxation parameter, , was varied from a maximum of to a minimum of . For the largest value, , significantly more melting occured near the top wall. Between the three smallest values of ranging from to , there appeared to be little sensitivity. The quadrature degree, , was varied from a minimum of to a maximum of . For the smallest value, , significantly more melting occured near the top wall. Between the three largest values of ranging from to , there appeared to be little sensitivity.
5.1.3 Results and discussion
Figure 5 shows the temperature field, streamlines of the velocity field, and the postprocessed position of the phase interface at times and . Convection in the liquid phase caused a sharper temperature gradient near the top of the phase interface than towards its bottom. This in turn caused a faster melting rate in the upper region. This effect was more pronounced at the later time.
At the final time, after a total of seventynine time steps, 6971 Newton iterations were used in total. This averages to about ninety Newton iterations per time step. This includes all iterations from intermediate problems during the continuation procedure sketched in Figure 1. For all time steps, the largest intermediate regularization parameter value was . For the first time step, eight additional values ( = 0.257, 0.1295, 0.06575, 0.033875, 0.0179375, 0.013953125, 0.0119609375, 0.00996875) were needed between and the targeted value, . At the final time, an additional four values of ( = 0.097625, 0.0816875, 0.07371875, 0.005984375) were needed.
About halfway through the run time, at , the simulated and experimental interface positions showed good agreement both near the top and near the bottom walls. There were minor deviations in the interior. At the final time, there was good agreement for a large portion of the phase interface, but it proceeded more quickly in the experiment than in the simulation near the top and bottom walls. The previously conducted sensitivity study, shown in Figure 4, suggests that the phase interface near the top wall could agree better with the experimental data by reducing the mesh cell size , but no such improvement should be expected near the bottom wall. It has already been reported in Rakotondrandisa et al. (2019) that there is some uncertainty about the insulation of the top and bottom walls in the experiment. The phase interface as observed in the experiment is clearly not orthogonal to the top wall. A further analysis of the impact of various formulations for the boundary conditions, however, is outside of the scope of the current work. On the interior of the domain, the results compare well to the experiment. One possible way forward could be to apply a heat flux boundary condition on the top wall instead of the adiabatic boundary condition which was used. For further validation, the experimental conditions should be examined in more detail, especially with regards to the boundary conditions.
5.2 Freezing water
The freezing of distilled water was simulated for comparison to benchmarking experiments in Kowalewski and Rebow (1999), where a cube of liquid water was frozen from one side. For waterice, the density, heat capacity, and thermal conductivity vary significantly between phases. Furthermore, liquid water’s thermal expansion coefficient is nonlinear, with the sign inverting at the temperature of water’s density anomaly.
5.2.1 Simulation setup
The edge length of the experimental cube’s test section, taken from Kowalewski and Rebow (1999), was used as the length scale, i.e. . A realistic value for the liquid kinematic viscosity, as taken from Michałek and Kowalewski (2003), furthermore yields the characteristic time scale . Based on physical parameters provided in Danaila et al. (2014) and according to the previously introduced definitions (13), the similarity parameters are Gr = , Pr = 6.99, and Ste = 0.125. Phasedependent thermal conductivity and volumetric heat capacity were included with and in (8) and (9). The nonlinear water density model from Danaila et al. (2014) was used. First proposed in Gebhart and Mollendorf (1977), the key feature is the density anomaly of water. For example, the density anomaly caused the double convection pattern shown in Figure 6. The density is expanded around its maximum which occurs at .
The simulation initial values correspond to the “warm start” initial conditions from Kowalewski and Rebow (1999). For this, the left and right walls were respectively kept constant at hot and cold temperatures and . A steady state convection was reached. For the simulation, this initial temperature range was used for the temperature scale, i.e. . With this scale, the dimensionless boundary temperatures were and . The steady state problem was solved directly. Therefore, the liquid volume fraction was set to a constant and time derivatives were set to zero. In this case, the momentum (11) and energy (12) equations reduce to
(17)  
(18) 
These equations were solved with the same mixed finite elements and Newton method as from Section 3. The high Grashof number, Gr, caused the nonlinear solver to diverge. To handle this, the relaxation and continuation algorithm from Figure 1 was applied to Gr instead of . For this problem, Newton’s method converged optimally when solving in sequence Gr/16, Gr/8, Gr/4, Gr/2, and finally Gr. In the Sapphire code Zimmerman (2019), both Gr and continuation are handled by the same routine which allows the regularization parameter to be specified. To begin the time dependent simulation, the cold wall was dropped to while the hot wall remained at . This caused freezing from right to left. Note that the temperature scaling was unchanged, so the physical temperature of the cold wall was , as it was in the experiment.
5.2.2 Sensitivity to numerical parameters
For this sensitivity study, the nominal parameter values were , , , , and . From their nominal values, each parameter was varied independently. The resulting phase interfaces were compared at time . The results are shown in Figure 7. The time is also the final simulated time presented in Figure 8.
The mesh cell size, , was varied from a maximum of to a minimum of . The largest size, , significantly degraded the solution. The solution appeared to converge asymptotically as was reduced. Between the smallest two values, and , there was almost no change near either the top or bottom wall, and only a minor difference in the center. This result is better than was seen for the octadecane melting simulation in Section 5.1, where there was a much larger sensitivity to when using the same values. The time step size, , was varied from a maximum of to a minimum of . Nearly no sensitivity was observed between these settings. The regularization parameter, , was varied from a maximum of to a minimum of . The solution appeared to be approaching an asymptotic limit. The difference between the two smallest values, and , was barely visible. The solid velocity relaxation parameter was varied from a maximum of to a minimum of . For the largest value , significantly more freezing occured near the center and less freezing occured near the bottom wall. The solution appeared to be approaching an asymptotic limit. Between the two smallest values of , and , the change is barely visible. The quadrature degree, , was varied between , , and . This had no significant effect on the solution.
5.2.3 Results and discussion
The sensitivity study revealed that large time step sizes can be used without significantly changing the resulting phase interface. To most easily compare to experimental results in Kowalewski and Rebow (1999), solutions were obtained at every onehundred physical seconds, corresponding to a simulated time step size of about . Figure 8 visualizes results at times , , and . As freezing proceeded from right to left, the two circulating regions of natural convection were maintained, translating to the left along with the phase interface. The bottom circulating region showed a lower temperature and velocity. The freezing front was nearly planar and proceeded more rapidly in this region. The top circulating region had higher temperature and velocity. The freezing front was curved and proceeded less rapidly near the top wall. Also, near the phase interface, note the sharper temperature gradient near the top wall. This slowed the freezing process in that region.
At , after a total of twentythree time steps, 1812 Newton iterations were used in total. This averages to about eighty Newton iterations per time step. This includes all iterations from intermediate problems during the continuation procedure sketched in Figure 1. For all time steps, the largest intermediate regularization parameter value was . For the first time step, four intermediate values of ( = 2.05, 1.027, 0.5155, 0.25975) were needed between and the targeted value, . At , an additional four values of ( = 0.131875, 0.09990625, 0.0679375, 0.03596875) were needed.
Figure 9 compares the simulation results to a group of experimental runs from Kowalewski and Rebow (1999). The simulation compares remarkably well to experimental run #1; but that run reportedly used a “cold start”. The simulation used a “warm start”, which corresponds to the experimental runs #4 and #5. It was also shown in Kowalewski and Rebow (1999) that the freezing process was sensitive to the thermal control of the test section walls. Experimental runs #1 and #4 had a constant air flow around the walls. For run #5, the test section was submerged in a water bath. The importance of simulating threedimensional heat transfer in the side walls was further demonstrated in Giangi et al. (2000). There, it was also noted that the viscosity of water substantially increases near its freezing temperature, and that this in return substantially affects the freezing process. Further investigating the effects of a temperature dependent viscosity or heat transfer through the test section walls was outside the scope of this work.
6 Conclusions
An enthalpy method was used to simulate convectioncoupled isothermal phasechange on a single geometric domain. The governing equations were discretized in space with mixed finite elements and in time with backward difference formulas. The resulting system of nonlinear equations was solved with Newton’s method. As new contributions, the current work

simplified the stabilization of the mixed finite element formulation. In Danaila et al. (2014); Zimmerman and Kowalski (2017), both TaylorHood elements and the pressure penalty method were used in combination. In the current work, equal order elements were stabilized with the pressure penalty method alone.

designed an efficient continuation procedure capable of robustly solving the nonlinear problem for a variety of physical and numerical parameters. Without this procedure, previous literature Danaila et al. (2014); Zimmerman and Kowalski (2017) relied on adaptive mesh refinement, carefully selected time step sizes, specifying a oneway melting or solidification process, and initializing the new phase in part of the domain.

verified the solver’s accuracy to second order in space and time via a convergence study, using the method of manufactured solutions.

presented simulations for the melting of octadecane and one for the freezing of water. For each, a sensitivity study was performed with respect to five numerical parameters. These were the mesh cell size , the time step size , the phase interface regularization parameter , the solid velocity relaxation parameter , and the quadrature degree . In both cases, it was demonstrated that the results could be obtained which were not sensitive to , , , or . For the water freezing case, sensitivty to was also minor between the two smallest values. For the octadecane melting case, results were still visibly sensitive to , though an asymptotic convergence behavior was observed. Further reducing was impractical, because it would increase the number of degrees of freedom in the linear system quadratically. This would be especially computationally intensive for the employed direct linear system solver.

compared results from the two simulation cases to experimental data sets from the literature, namely Okada (1984) for melting octadecane and Kowalewski and Rebow (1999) for freezing water. For the case of melting octadecane, the simulated phase interface shows good agreement to the experiment in the interior, while there were large discrepancies near the adiabatic walls. The water freezing simulation was compared to multiple experimental runs. The runs use varied procedures for controlling initial conditions and thermal regulation of the test section walls. While the simulated phase interface compared well to one of the runs, that run reportedly used a different procedure for initialization which is not considered by the simulation. More information is needed to accurately define the initial conditions for the water freezing case, and adiabatic wall assumptions are too ideal for both cases.

shared the open source code in a new Python packaged called Sapphire Zimmerman (2019), using the finite element library Firedrake Rathgeber et al. (2016). Sapphire’s test suite covers many of the results of this paper. Reproducing exact results is only a matter of setting the described parameter values. The versions of Sapphire and Firedrake were documented with DOI’s on Zenodo, respectively at Zimmerman (2019) and zenodo/Firedrake20190620.0 (2019).
The implementation in Sapphire Zimmerman (2019) is flexible, and can readily be applied to problems with other geometries, parameters, initial values, and boundary conditions. By using UFL Alnæs et al. (2012) for defining the mathematical model symbolically, the model is easily modifiable and extensible. As a next step, we want to apply the approach to nonisothermal phasechange involving mushy layers, e.g. melting and refreezing cycles of salt water. The continuation procedure has been robust for convectioncoupled isothermal phasechange, and there is reason to believe that it will also work well, if not better, for the case of nonisothermal phasechange. Still, there may exist a more general approach to regularizing such nonlinear problems, and this is also worth exploring. There are also promising routes for reducing computational costs. One route would be the application of iterative linear system solvers, which would require the development of a preconditioner. Such an iterative solver would make threedimensional geometries tractable, and would also open opportunities for addressing inverse problems.
Acknowledgments
The authors were funded in part by the Excellence Initiative of the German Federal and State Governments through grant GSC 111. The work was furthermore supported by the Federal Ministry of Economic Affairs and Energy, on the basis of a decision by the German Bundestag (50 NA 1502).
References
References
 Sparrow et al. (1978) E. Sparrow, R. R. Schmidt, J. W. Ramsey, Experiments on the role of natural convection in the melting of solids, Journal of Heat Transfer 100 (1978).
 Sparrow et al. (1979) E. Sparrow, J. W. Ramsey, R. G. Kemink, Freezing controlled by natural convection, Journal of Heat Transfer 101 (1979).
 Okada (1984) M. Okada, Analysis of heat transfer during melting from a vertical wall, International Journal of Heat and Mass Transfer 27 (1984) 2057–2066.
 Kowalewski and Rebow (1999) T. A. Kowalewski, M. Rebow, Freezing of water in a differentially heated cubic cavity, International Journal of Computational Fluid Dynamics 11 (1999) 193–210.
 Schüller et al. (2017) K. Schüller, B. Berkels, J. Kowalski, Integrated modeling and validation for phase change with natural convection, in: M. Schäfer, M. Behr, M. Mehl, B. Wohlmuth (Eds.), Recent Advances in Computational Engineering, volume 124 of Lectures Notes in Computational Science and Engineering (LNCSE), Springer, 2017, pp. 127–144.
 Voller et al. (1987) V. Voller, M. Cross, N. Markatos, An enthalpy method for convection/diffusion phase change, International Journal for Numerical Methods in Engineering 24 (1987) 271–284.
 Alexiades and Solomon (1992) V. Alexiades, A. Solomon, Mathematical modeling of melting and freezing processes, Bristol, PA (United States); Hemisphere Publishing, 1992.
 Voller and Prakash (1987) V. R. Voller, C. Prakash, A fixed grid numerical modelling methodology for convectiondiffusion mushy region phasechange problems, International Journal of Heat and Mass Transfer 30 (1987) 1709–1719.
 Brent et al. (1988) A. Brent, V. Voller, K. Reid, Enthalpyporosity technique for modeling convectiondiffusion phase change: application to the melting of a pure metal, Numerical Heat Transfer, Part A: Applications 13 (1988) 297–318.
 Giangi et al. (2000) M. Giangi, T. A. Kowalewski, F. Stella, E. Leonardi, Natural convection during ice formation: numerical simulation vs. experimental results, Computer Assisted Mechanics and Engineering Sciences 7 (2000) 321–342.
 Evans and Knoll (2007) K. J. Evans, D. A. Knoll, Temporal accuracy analysis of phase change convection simulations using the JFNKSIMPLE algorithm, International Journal for Numerical Methods in Fluids 55 (2007) 637–653.
 Belhamadia et al. (2012) Y. Belhamadia, A. S. Kane, A. Fortin, An enhanced mathematical model for phase change problems with natural convection, International Journal of Numerical Analysis and Modeling 3 (2012) 192–206.
 Danaila et al. (2014) I. Danaila, R. Moglan, F. Hecht, S. Le Masson, A Newton method with adaptive finite elements for solving phasechange problems with natural convection, Journal of Computational Physics 274 (2014) 826–840.
 Zimmerman and Kowalski (2017) A. G. Zimmerman, J. Kowalski, Monolithic simulation of convectioncoupled phasechange: Verification and reproducibility, in: M. Schäfer, M. Behr, M. Mehl, B. Wohlmuth (Eds.), Recent Advances in Computational Engineering, volume 124 of Lectures Notes in Computational Science and Engineering (LNCSE), Springer, 2017, pp. 177–197.
 Rakotondrandisa et al. (2019) A. Rakotondrandisa, I. Danaila, L. Danaila, Numerical modelling of a meltingsolidification cycle of a phasechange material with complete or partial melting, International Journal of Heat and Fluid Flow 76 (2019) 57–71.
 Wang et al. (2010) S. Wang, A. Faghri, T. L. Bergman, A comprehensive numerical model for melting with natural convection, International Journal of Heat and Mass Transfer 53 (2010) 1986–2000.
 Voller et al. (1990) V. Voller, C. Swaminathan, B. G. Thomas, Fixed grid techniques for phase change problems: a review, International Journal for Numerical Methods in Engineering 30 (1990) 875–898.
 Woodfield et al. (2019) J. Woodfield, M. Alvarez, B. GómezVargas, R. RuizBaier, Stability and finite element approximation of phase change models for natural convection in porous media, Journal of Computational and Applied Mathematics 360 (2019) 117–137.
 Donea and Huerta (2003) J. Donea, A. Huerta, Finite element methods for flow problems, John Wiley & Sons, 2003.

Ascher and Petzold (1998)
U. M. Ascher, L. R. Petzold, Computer methods for ordinary differential equations and differentialalgebraic equations, volume 61, SIAM, 1998.
 Rathgeber et al. (2016) F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. T. McRae, G.T. Bercea, G. R. Markall, P. H. J. Kelly, Firedrake: automating the finite element method by composing abstractions, ACM Transactions on Mathematical Software 43 (2016) 24:1–24:27.
 Zimmerman (2019) A. G. Zimmerman, geofluiddynamics/sapphire: Convectioncoupled phasechange verification and preliminary validation, https://doi.org/10.5281/zenodo.3250317, 2019.
 Alnæs et al. (2012) M. S. Alnæs, A. Logg, K. B. Ølgaard, M. E. Rognes, G. N. Wells, Unified form language: A domainspecific language for weak formulations of partial differential equations, ACM Transactions on Mathematical Software 40 (2012) 9.
 Balay et al. (2018) S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, PETSc Users Manual, Technical Report ANL95/11  Revision 3.9, Argonne National Laboratory, 2018.
 Balay et al. (1997) S. Balay, W. D. Gropp, L. C. McInnes, B. F. Smith, Efficient management of parallelism in object oriented numerical software libraries, in: E. Arge, A. M. Bruaset, H. P. Langtangen (Eds.), Modern Software Tools in Scientific Computing, Birkhäuser Press, 1997, pp. 163–202.
 Dalcin et al. (2011) L. D. Dalcin, R. R. Paz, P. A. Kler, A. Cosimo, Parallel distributed computing using Python, Advances in Water Resources 34 (2011) 1124–1139. New Computational Methods and Software Tools.
 Amestoy et al. (2001) P. Amestoy, I. Duff, J. L’Excellent, J. Koster, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM Journal on Matrix Analysis and Applications 23 (2001) 15–41.
 Amestoy et al. (2006) P. Amestoy, A. Guermouche, J. L’Excellent, S. Pralet, Hybrid scheduling for the parallel solution of linear systems, Parallel Computing 32 (2006) 136–156.
 zenodo/Firedrake20190620.0 (2019) zenodo/Firedrake20190620.0, Software used in ’Simulating convectioncoupled phasechange in enthalpy form with mixed finite elements’, https://doi.org/10.5281/zenodo.3250636, 2019.
 Roache (2002) P. J. Roache, Code verification by the method of manufactured solutions, Journal of Fluids Engineering 124 (2002) 4–10.
 Michałek and Kowalewski (2003) T. Michałek, T. A. Kowalewski, Simulations of the water freezing process–numerical benchmarks, Task Quarterly 7 (2003) 389–408.
 Gebhart and Mollendorf (1977) B. Gebhart, J. C. Mollendorf, A new density relation for pure and saline water, Deep Sea Research 24 (1977) 831–848.
Comments
There are no comments yet.