Mixed finite elements for convection-coupled phase-change in enthalpy form: Open software verified and applied to 2D benchmarks

TL;DR

This paper presents a verified open-source mixed finite element software for accurately simulating convection-coupled phase-change processes in 2D, with a novel continuation approach and validation against experimental data.

Contribution

It introduces a separation of mesh refinement from nonlinear regularization and a continuation procedure for robust solutions, along with a verified open-source implementation.

Findings

Code accurately solves 2D phase-change problems with verified convergence.

Benchmark simulations match experimental data and literature results.

Numerical errors are effectively reduced through parameter sensitivities.

Abstract

Melting and solidification processes are often affected by natural convection of the liquid, posing a multi-physics problem involving fluid flow, convective and diffusive heat transfer, and phase-change reactions. Enthalpy methods formulate this convection-coupled phase-change problem on a single computational domain. The governing equations can be solved accurately with a monolithic approach using mixed finite elements and Newton's method. Previously, the monolithic approach has relied on adaptive mesh refinement to regularize local nonlinearities at phase interfaces. This contribution instead separates mesh refinement from nonlinear problem regularization and provides a continuation procedure which robustly obtains accurate solutions on the tested 2D uniform meshes. A flexible and extensible open source implementation is provided. The code is formally verified to accurately solve the governing equations in time and in 2D space, and convergence rates are shown. Two benchmark simulations are presented in detail with comparison to experimental data sets and corresponding results from the literature, one for the melting of octadecane and another for the freezing of water. Sensitivities to key numerical parameters are presented. For the case of freezing water, effective reduction of numerical errors from these key parameters is successfully demonstrated. Two more simulations are briefly presented, one for melting at a higher Rayleigh number and one for melting gallium.

Tables

Table 1. Table 1: Coefficients for backward difference formulas up to third order, taken from [ 22 ] . For example, the backward Euler method is given by k = 1 𝑘 1 k=1 , and BDF2 is given by k = 2 𝑘 2 k=2 .

$k$	$a_0$	$a_1$	$a_2$	$a_3$
1	1	-1
2	3/2	-2	1/2
3	11/6	-3	3/2	-1/3

Table 2. Table 2: Empirical convergence of second-order spatial and temporal discretizations. The errors for discrete pressure p h subscript 𝑝 ℎ p_{h} , velocity 𝐮 h subscript 𝐮 ℎ \mathbf{u}_{h} , and temperature T h subscript 𝑇 ℎ T_{h} are computed in their natural norms with respect to the manufactured solutions p M subscript 𝑝 𝑀 p_{M} , 𝐮 M subscript 𝐮 𝑀 \mathbf{u}_{M} , and T M subscript 𝑇 𝑀 T_{M} to compute convergence rates r h , p subscript 𝑟 ℎ 𝑝 r_{h,p} , r h , 𝐮 subscript 𝑟 ℎ 𝐮 r_{h,\mathbf{u}} , and r h , T subscript 𝑟 ℎ 𝑇 r_{h,T} . The time discretization errors are computed in the L 2 superscript 𝐿 2 L^{2} -norm to determine convergence rates r Δ ​ t , 𝐮 subscript 𝑟 Δ 𝑡 𝐮 r_{\Delta t,\mathbf{u}} and r Δ ​ t , T subscript 𝑟 Δ 𝑡 𝑇 r_{\Delta t,T} .

$h$	${\parallel p_h-p_M\parallel }_{L^2}$	$r_{h,p}$	${\parallel {𝐮}_h-{𝐮}_M\parallel }_{H^1}$	$r_{h,𝐮}$	${\parallel T_h-T_M\parallel }_{H^1}$	$r_{h,T}$
1/32	4.748e-01		1.910e-02		1.534e-03
1/64	3.041e-02	3.965	4.537e-03	2.074	3.791e-04	2.017
1/128	1.912e-03	3.991	1.120e-03	2.018	9.447e-05	2.005
1/256	1.205e-04	3.988	2.796e-04	2.002	2.360e-05	2.001

Equations

$$f_l=\{\begin{array}{ll}0,& T<T_m\\ 1,& T>T_m\end{array}$$

$$\zeta (T)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(\frac{-T^2}{2{\sigma }^2}\right)$$

$${\varphi }_l=(f_l\ast \zeta )(T)=\frac{1}{2}\left(1+\mathrm{erf}\left(\frac{T-T_m}{\sigma \sqrt{2}}\right)\right)$$

$$\nabla \cdot \text{mathpzc}\mathbf{u}$$

$${\rho }_0\left({\partial }_{\text{mathpzc}t}\text{mathpzc}\mathbf{u}+\nabla \text{mathpzc}\mathbf{u}\cdot \text{mathpzc}\mathbf{u}+\frac{{\varphi }_s}{d}\text{mathpzc}\mathbf{u}\right)+\nabla \text{mathpzc}p-2\nabla \cdot (\mu \mathrm{sym}\nabla \text{mathpzc}\mathbf{u})-\rho \mathbf{g}$$

$${\partial }_{\text{mathpzc}t}\left(\text{mathpzc}C(T-T_0)\right)+{\rho }_{l}L{\partial }_{\text{mathpzc}t}{\varphi }_l+\text{mathpzc}\mathbf{u}\cdot \nabla (\text{mathpzc}C(T-T_0))-\nabla \cdot (\text{mathpzc}k\nabla \text{mathpzc}T)$$

$$\mathbf{x}=\frac{\text{mathpzc}\mathbf{x}}{\mathrm{X}},t=\frac{\text{mathpzc}t\mathrm{U}}{\mathrm{X}},p=\frac{\text{mathpzc}p}{{\rho }_l{\mathrm{U}}^2},\mathbf{u}=\frac{\text{mathpzc}\mathbf{u}}{\mathrm{U}},T=\frac{\text{mathpzc}T-\text{mathpzc}{T}_m}{\delta \mathrm{T}},$$

$$\kappa$$

$$C$$

$$\nabla \cdot \mathbf{u}=0,$$

$${\partial }_t\mathbf{u}+\nabla \mathbf{u}\cdot \mathbf{u}+\nabla p-\frac{2}{\mathrm{Re}}\nabla \cdot \mathrm{sym}\nabla \mathbf{u}+\frac{\mathrm{Ra}}{\mathrm{Pr}}b\hat{\mathbf{g}}+\frac{{\varphi }_s}{\tau }\mathbf{u}=0,$$

$${\partial }_t(CT)+\frac{1}{\mathrm{Ste}}{\partial }_t{\varphi }_l+\mathbf{u}\cdot \nabla (CT)-\frac{1}{\mathrm{RePr}}\nabla \cdot \left(\kappa \nabla T\right)=0$$

$$\mathrm{Re}=\frac{\mathrm{UX}}{{\nu }_l},\mathrm{Ra}=\frac{g{\beta }_0{\mathrm{X}}^3\delta \mathrm{T}}{{\nu }_l{\alpha }_l},\mathrm{Pr}=\frac{{\nu }_l}{{\alpha }_l},\mathrm{Ste}=\frac{c_l\delta \mathrm{T}}{L},$$

$$\mathcal{F}=\left({\psi }_p,\nabla \cdot \mathbf{u}\right)+\left({\mathbf{\psi }}_u,{\partial }_t\mathbf{u}+\nabla \mathbf{u}\cdot \mathbf{u}+\frac{\mathrm{Ra}}{\mathrm{Pr}}b\hat{\mathbf{g}}+\frac{1}{\tau }{\varphi }_s\mathbf{u}\right)$$

$$-\left(\nabla \cdot {\mathbf{\psi }}_u,p\right)+\frac{2}{\mathrm{Re}}\left(\mathrm{sym}\nabla {\mathbf{\psi }}_u,\mathrm{sym}\nabla \mathbf{u}\right)$$

$$+\left({\psi }_T,{\partial }_t(CT)+\frac{1}{\mathrm{Ste}}{\partial }_t{\varphi }_l+\mathbf{u}\cdot \nabla (CT)\right)+\frac{1}{\mathrm{RePr}}\left(\nabla {\psi }_T,\kappa \nabla T\right)$$

$${\partial }_{t}w=\frac{1}{\Delta t}\sum _{i=0}^k{a}_i{w}^{n-i}$$

$$D_{\mathbf{w}}\mathcal{F}({\mathbf{w}}_j)({\mathbf{w}}_{j+1}-{\mathbf{w}}_j)=\mathcal{F}({\mathbf{w}}_j),$$

$${\mathbf{u}}_M$$

$$p_M$$

$$T_M$$

$$\nabla \mathbf{u}\cdot \mathbf{u}+\nabla p-\frac{2}{\mathrm{Re}}\nabla \cdot \mathrm{sym}\nabla \mathbf{u}+\frac{\mathrm{Ra}}{\mathrm{Pr}}b\hat{\mathbf{g}}=0,$$

$$\mathbf{u}\cdot \nabla T-\frac{1}{\mathrm{RePr}}\nabla \cdot \left(\kappa \nabla T\right)=0$$

Full text

††journal: arXiv

1 Introduction

Melting and solidification processes in phase-change 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 phase-change processes in the cryosphere, e.g. as relevant for climate projections. In some regimes, the phase-change process is substantially affected by convection in the liquid phase, which has been demonstrated in a number of experiments [1, 2, 3, 4, 5].

Any spatiotemporally resolved simulation of this multi-physics process is a mathematical and computational challenge. Common modeling strategies couple systems of multi-parameter 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 $T_{m}$. The corresponding phase-change process is commonly referred to as isothermal phase-change [6]. Isothermal phase-change 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, i.e. regions containing liquid and solid phases. 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 phase-change 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 multi-phase gas-liquid applications. In these cases, solids exist only as boundaries and are usually fixed in space and time. Simulating convection-coupled melting and freezing requires the extension of CFD into multi-phase liquid-solid 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 [6, 7]. 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 post-processed as an explicit function of the solution.

The current work uses an interface capturing method that relies on an enthalpy formulation of the phase-change process. This is often referred to as an enthalpy method [6]. Generally, enthalpy methods write the energy balance in terms of both the temperature and the enthalpy, and hence must be closed with an equation relating these two quantities. For incompressible materials, this relation is well-understood and can be easily phrased in terms of the liquid and solid phase volume fractions. For isothermal phase-change processes with liquid and solid regions respectively above and below the melting/freezing temperature, the phase can be written as a function of only the temperature. In these cases, the enthalpy can be eliminated from the system as an unknown. This choice results in an energy balance taking the form of a standard convection-diffusion 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 convection-coupled phase-change [6, 8, 9, 7, 10, 11, 12, 13, 14, 15, 16], but there remain many opportunities for improvement, e.g. regarding their dependency on numerical parameters and computational feasibility in three spatial dimensions.

In the context of enthalpy methods, a variety of techniques have been applied to modeling the velocity in the solid phase. These are sometimes referred to as a solid velocity correction [17]. Different approaches are compared in [6, 18]. Three primary approaches to solid velocity correction have prevailed:

modifying the algorithm to set the velocity to zero in the solid [17], 2. 2.

adding a phase-dependent source term to the momentum equation, forcing the velocity to relax to zero in the solid [12, 5, 15], and 3. 3.

prescribing a phase-dependent viscosity with orders of magnitude larger values in the solid than in the liquid [13, 14].

Of these, only the third approach can physically model solid regions which are not attached to stationary boundaries. The first approach provides a 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 additional parameters must be calibrated based on experimental data, as seen in [13, 15, 19]. 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.

For any chosen approach to solid velocity correction, the abrupt change in velocity at phase interfaces disrupts nonlinear solver convergence [9]. A common regularization approach is to introduce a temperature range in which the phase-change 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 non-isothermal phase-change for non-pure 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 as numerical parameters to stabilize an otherwise unstable computational model. The current work includes sensitivity studies where the mesh size is held constant and the phase interface regularization is reduced until it does not substantially affect the size of the artificial mushy region. The minimum extent of the artificial mushy region will always be constrained by the size of the mesh cells.

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 [6, 8, 9], finite volumes [10, 11, 17, 5], the finite element method (FEM) with operator splitting [12], and monolithic mixed FEM [13, 14, 15, 20, 16, 21, 19]. Mixed FEM methods come both in purely primal forms [13, 14, 15, 20, 16, 19] and also in forms with auxiliary variables [16, 21]. The cited primal formulations all used Taylor-Hood elements for the pressure-velocity system. Of these, [20, 16, 19] used quadratic elements for the temperature, making the spatial discretization method theoretically second order accurate for all of the solution fields. Only recently has formal analysis of the stability, and numerical accuracy verification, of mixed FEM schemes for enthalpy formulated convection-coupled phase-change begun [16, 21]. Both articles use the method of manufactured solutions to verify the spatial order of accuracy. Temporal accuracy is also verified using manufactured solutions in [16]. Typically, temporal discretization is performed with finite differences. Second-order accurate fully implicit finite difference temporal discretization was advocated in [11] and subsequently used in [12, 15, 16, 19].

Based on the primal mixed finite element approach in [13], an open source implementation was published in [14]. Both [13] and [14] used adaptive mesh refinement (AMR), albeit with alternative formulations. While results in [14] were promising, it was later difficult to apply its presented implementation to solidification problems. Furthermore, while that implementation could adaptively refine meshes, adaptive coarsening was not implemented. The lack of mesh coarsening capability was computationally demanding for simulations run over a long simulated time with the phase interface moving large distances through the domain.

There are some drawbacks to any AMR approach. For example, the benefit of reduced degrees of freedom competes with the cost of dynamic re-meshing, especially given that the phase interface moves throughout the spatial domain in time. Furthermore, the benefit of per-cell accuracy from locally refined meshes competes with accumulating interpolation errors which come from transferring solutions between non-matching meshes of different time steps. Still, application of AMR methods to convection-coupled phase-change is a promising, successful, and active (e.g. [15, 20, 19]) field of research,.

Another key limitation, noted in [14], is that AMR requires an a posteriori error estimate, which in turn requires first obtaining a solution on a coarse or manually refined mesh. Convergence of the nonlinear solution methods presented in [13] and [14] both required careful preparation of the problem setting, in particular relying on manual local refinement of the initial mesh, carefully chosen time step sizes, and artificial initialization of the new phase in one-way melting or solidification processes. The current paper sets aside AMR to focus on more reliably obtaining solutions on uniform grids, including relatively coarse grids which could be used, for example, to automatically initialize AMR methods.

This article is structured as follows: Section 2 presents an enthalpy formulation for convection-coupled phase-change processes with a generalized approach to phase-interface regularization. Section 3 presents a monolithic solver approach with mixed finite elements, implicit time discretization, and Newton’s method. Motivated by [9, 14], a continuation procedure is introduced with the goal of robustly solving the nonlinear problem without local mesh refinement, which is made difficult by the irregularity at the phase interface. Section 4 introduces a new open-source implementation and describes the open-source libraries on which it depends. Furthermore, empirical convergence rates are shown for the spatial and temporal discretizations using the method of manufactured solutions. Section 5 presents benchmark simulation results for the melting of octadecane and freezing of water, both on unit square geometries. For both, sensitivities to numerical parameters are reported and the resulting time evolutions of the phase interface are compared to experimental data from the literature. Furthermore, brief comparisons are made to the simulation results in [19]. Section 6 briefly presents two additional example simulations on rectangular geometries, one for melting octadecane at a higher Rayleigh number and one for melting gallium.

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 two-way coupling between the fluid dynamics and the evolving shapes of phase interfaces. It is assumed that all solid regions are attached to stationary boundaries. 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 phase-change [6, 9]. This means that melting and solidification occur at a single temperature.

2.2 Phase interface regularization

Consider a bulk phase-change material which can either be fully liquid, fully solid, or have both liquid and solid regions. Any small reference volume of this material has a liquid volume fraction $f_{l}$ and a solid volume fraction $f_{s}$. Assuming there are no other phases, e.g. no gas inclusions, $f_{l}+f_{s}=1$. When the material is fully liquid, $f_{l}=1$. When it is fully solid, $f_{l}=0$. In this case, the value of $f_{l}$ can be used as a substitute for the phase state.

Isothermal phase-change processes are characterized by having distinct interfaces between phases. At these interfaces, the temperature of the material is at the melting/freezing point $T=T_{m}$. In liquid regions, $T>T_{m}$. In solid regions, $T<T_{m}$. This means that the liquid volume fraction $f_{l}$ can be reconstituted from the temperature alone, i.e.

[TABLE]

With this model, phase interfaces are theoretically infinitesimally thin lines in two-dimensional space or surfaces in three-dimensional space. The phase interfaces can be regularized by convoluting $f_{l}$ with a Gaussian kernel, e.g.

[TABLE]

yielding

[TABLE]

This introduces a parameter, $\sigma$, which is the standard deviation of the Gaussian and is always positive. As $\sigma$ approaches zero, $\phi_{l}$ approaches $f_{l}$. With too large of $\sigma$, substantial regions may develop where $0<\phi_{l}<1$. In the context of isothermal phase-change with distinct phase interfaces, such regions are numerical artifacts and will be referred to as artificial mushy regions. Alternative regularizations of phase interfaces that rely on similar parameters have been proposed in other literature [13, 14, 15]. These approaches often highlight that the artificial mushy region is constricted to a temperature interval around $T_{m}$. Using (3), that interval can be approximated as $\left[-2\sigma,2\sigma\right]$, i.e. within two standard deviations of $T_{m}$. It is important to note that $\sigma$ is a numerical parameter and does not carry physical meaning. Therefore, for results in Section 5, $\sigma$ was reduced, for a given mesh, until regions where $0<\phi_{l}<1$ were small and until the shapes of the $T=T_{m}$ contours were not affected.

2.3 Governing equations

Convection-coupled phase-change can be simulated in enthalpy form [6] using the following governing equations, which are balances of mass, momentum, and energy

[TABLE]

The independent unknowns are pressure $\mathpzc{p}$, velocity $\mathbf{\mathpzc{u}}$, and temperature $\mathpzc{T}$. The unknown liquid volume fraction $\phi_{l}$, solid volume fraction $\phi_{s}$, and Boussinesq density $\rho$ always depend on $\mathpzc{T}$. The thermal conductivity $\mathpzc{k}$ and volumetric heat capacity $\mathpzc{C}$ can both depend on the phase and therefore depend on $\mathpzc{T}$. The volume-averaged values are $\mathpzc{k}=\phi_{l}\mathpzc{k}_{l}+\phi_{s}\mathpzc{k}_{s}$ and $\mathpzc{C}=\phi_{l}\rho_{l}\mathpzc{c}_{l}+\phi_{s}\rho_{s}\mathpzc{c}_{s}$.

The mass equation (4) is the continuity equation for an incompressible material. The momentum equation (5) is the incompressible Navier-Stokes momentum equation with Boussinesq force $\rho\mathbf{g}$ and an additional solid velocity correction [17] term. This term, $\phi_{s}\mathbf{\mathpzc{u}}/d$, with $d\ll 1$, relaxes the velocity toward zero in the solid region, where $\phi_{s}$ is the regularized solid volume fraction and $\phi_{l}+\phi_{s}=1$. This particular choice of solid velocity correction was used for example in [9] and [11]. Under the Boussinesq approximation, density variations are only considered with respect to the temperature. The incompressible flow considers a reference density $\rho_{0}$. The symmetric part of the rate-of-strain tensor is denoted by $\mathrm{sym}\nabla\mathpzc{u}$. The energy balance is in enthalpy form [6]. This invokes a reference temperature $\mathpzc{T_{0}}$ and produces a source term accounting for the phase-change with latent heat $L$.

2.4 Nondimensionalization

The independent variables were nondimensionalized with respect to characteristic scales similar to [13]. These include spatial scale $\mathrm{X}$, speed scale $\mathrm{U}$, and temperature scale $\delta\mathrm{T}$. The liquid material properties were used as reference values, e.g. the liquid kinematic viscosity $\nu_{l}$. Additionally, the reference temperature was chosen as the isothermal melting and freezing temperature, i.e. $\mathpzc{T}_{0}=\mathpzc{T}_{m}$. The nondimensional variables were chosen to be

[TABLE]

In this case, the nondimensional melting and freezing temperature is always $T_{m}=0$. The phase-dependent material properties were normalized with respect to the liquid values, i.e.

[TABLE]

In nondimensional form, the governing equations (4), (5), and (6) then read

[TABLE]

where the Reynolds, Rayleigh, Prandtl, and Stefan numbers are defined as

[TABLE]

In the solid velocity correction term, $\tau=d\mathrm{U}/\mathrm{X}$ is the nondimensional velocity relaxation factor. The Boussinesq buoyancy force, with unit vector $\hat{\mathbf{g}}$ in the direction of gravity, contains a general term $b=b(T)$ which must be specified for a given liquid density model. In cases where this density model is nonlinear, e.g. when accounting for the density anomaly observed in water, a constant reference value $\beta_{0}$ must be chosen for the thermal expansion coefficient. The liquid thermal diffusivity is $\alpha_{l}=\kappa_{l}/(\rho_{l}c_{l})$.

The nondimensional system of equations (10, 11, 12) constitute the governing equations for the remainder of this article. The scaling highlights the regimes of natural convection and phase-change. 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 length and speed scales $\mathrm{X}$ and $\mathrm{U}$ are chosen separately for each simulation. As explained in [13], common choices for convection-coupled melting and solidification are $\mathrm{U}=\nu_{l}/\mathrm{X}$ and $\mathrm{U}=\alpha_{l}/\mathrm{X}$. The former yields $\mathrm{Re}=1$ while the latter yields $\mathrm{Re}=1/\mathrm{Pr}$.

3 Numerical methods

The time-dependent problem was solved as a sequence of initial boundary values problems. Mixed finite elements were used for the spatial discretization and backward difference formulas were used for the temporal discretization. Non-homogeneous Dirichlet and homogeneous Neumann boundary conditions were applied to the temperature field, while homogeneous Dirichlet boundary conditions were applied to the velocity field. The nonlinear system was solved with Newton’s method. Reliably 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 [13, 14, 16], mixed finite elements were used to approximate the system (10, 11, 12) with basis functions $\psi_{p},\boldsymbol{\psi}_{u},\psi_{T}$, yielding the weak form residual

[TABLE]

where $(u,v)=\int uvd\mathbf{x}$ and $(\mathbf{u},\mathbf{v})=\int\mathbf{u}\cdot\mathbf{v}d\mathbf{x}$.

Piece-wise linear polynomials were used for approximating the pressure $p$ while piece-wise quadratic polynomials were used for the velocity $\mathbf{u}$ and temperature $T$. This forms a mixed element which is analogous to Taylor-Hood elements for the pressure-velocity subspace. The same element was used in [16], while only $\mathbf{u}$ (and not $T$) was quadratic in [13] and [14].

Homogeneous Dirichlet boundary conditions were always applied to the velocity, which reflect a stationary solid whenever it connects to the boundary, and a no-slip boundary for the liquid phase. Both non-homogeneous Dirichlet and homogeneous Neumann boundary conditions were applied to the temperature. No boundary conditions were applied to the pressure. Therefore, the pressure solution was only defined up to a constant. Successful solution required informing the linear solver of the appropriate nullspace, which was automated by the software in Section 4.1. To solution was made unique by subtracting the mean pressure (i.e. $\int p\ d\mathbf{x}$), which also therefore guaranteed zero mean pressure in the post-processed solution.

3.2 Temporal discretization

The time derivatives were discretized with the BDF2 method, i.e. the second order method from the family of $k^{\mathrm{th}}$ order constant time step size ($\Delta t$) backward difference formulas (BDF-$k$) [22]. For some unknown $w^{n}$ at time step $n$, the BDF-$k$ discretization can be written as

[TABLE]

where the coefficients $a_{i}$, shown in Table 1 depend only on $k$. This requires storing solutions from the previous $k$ time steps.

In addition to the independent unknowns, (15) was applied directly to $\phi_{l}$ and to $CT$ in (12). For example, applying the second order formula (BDF2) to $\phi_{l}$ yields $\partial_{t}\phi_{l}=(3\phi_{l}(T^{n})-4\phi_{l}(T^{n-1})+\phi_{l}(T^{n-2}))/\left(2\Delta t\right)$.

3.3 Nonlinear solution

Many terms in the residual (14) are highly nonlinear. The nonlinear system was solved with Newton’s method as follows. Consider the vector-valued system solution $\mathbf{w}^{n}=\begin{pmatrix}p^{n}&\mathbf{u}^{n}&T^{n}\end{pmatrix}^{\mathsf{T}}$ for any discrete time step $n$. The nonlinear problem’s solution is approximated for the next time step by solving a sequence of linear problems

[TABLE]

whose iterates $\mathbf{w}_{j+1}$ converge to $\mathbf{w}^{n}$. Here, $D_{\mathbf{w}}\mathcal{F}$ stands for the Gateaux derivative of $\mathcal{F}$. It can be derived analytically, which was done for a similar problem in [14]. Or as done in the current work, it can be symbolically computed using the software in Section 4. Each iteration of (16) is a linear system which can be solved robustly with a direct solver. At each new time step $n$ with unknown solution $\mathbf{w}^{n}$, the initial guess for the nonlinear solver was first taken as the known solution $\mathbf{w}^{n-1}$ from the previous time step $n-1$.

It is not simple to generate suitable initial guesses to reliably converge the nonlinear solver. Convergence is particularly difficult when using a small enough regularization parameter $\sigma$ in (3) such that the solution is not highly sensitive to the regularization. A continuation procedure, detailed in the next section, was developed to automatically search for better initial guesses when the nonlinear solver fails. The key idea of the procedure is to recognize that the regularity of nonlinear problem’s solution space is dominated by $\sigma$. Only $\sigma$, no other aspects of the problem, was modified during continuation. No data from previously solved time steps was modified. Entirely automated as part of the implementation in Section 4, the procedure successfully solved multiple benchmark problems on coarse uniform meshes with large time step sizes using small $\sigma$.

3.4 Continuation

Initial guesses were improved by solving intermediate problems with an increased value of $\sigma$. Increasing $\sigma$ improves the regularity of the solution, making the solver more likely to converge even for a poor-quality initial guess. Whenever a solution was found, it was used as an initial guess for solving the same system again, albeit with a smaller value of $\sigma$. The challenge hence was to find a sequence of $\sigma$ values resulting in solvable problems, with the final solution using the originally specified $\sigma$. The search process had two modes: Either the search was bounded by two values of $\sigma$ or it was unbounded. In the unbounded case, $\sigma$ was always doubled until the problem could be solved. In the bounded case, new $\sigma$ values were always inserted into the sequence to bisect the nearest solved and unsolved values. Figure 1 shows a detailed example of this process.

When applying this approach, the resulting sequence of intermediate $\sigma$ values depends on the exact nonlinear solver method used, including on its options such as convergence criteria and maximum iterations. It was found that using the relative residual as a convergence criterion would allow for significant accumulation of error in the absolute residual between time steps. For this reason, the relative tolerance convergence criterion was disabled, forcing the nonlinear solver to reach a low absolute residual or otherwise fail. Any single call to the nonlinear solver was only allowed twenty-four maximum iterations to reach the convergence criteria, meaning again that the solver would otherwise report failure.

4 Implementation and verification with open-source libraries

In order to reduce the required coding effort, leverage state-of-the-art 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 [23]. On top of Firedrake, for this work, the open source Sapphire [24] framework was developed which implements the methods from Section 3. Sapphire also provides an automated procedure for verifying implementations with the method of manufactured solutions (MMS). In addition to a general overview of Firedrake and Sapphire, the results of applying Sapphire’s MMS procedure to verify the current implementation are also presented here.

4.1 Firedrake

Firedrake [23] is a Python package which automates many aspects of finite element method software development. It uses the Unified Form Language (UFL) [25] to provide users with an abstract interface for rigorously specifying variational problems and their finite element discretizations. With the problem and discretization specified by the user, Firedrake employs the Two-Stage Form Compiler (TSFC) [26] to automatically write efficient code for the finite element assembly. The linear and nonlinear solvers are provided by PETSc [27, 28] via PETSc for Python [29]. For the results in this work, all linear systems were solved with MUMPS [30, 31, 32], 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 used during finite element assembly 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 $\sigma$, 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, $q$. Therefore, $q$ 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 archived with a DOI on Zenodo [33], ensuring the reproducibility of these results.

4.2 Sapphire

Sapphire is a Python package that was developed for the purposes of the current work. Sapphire provides a Python class for time dependent simulations that are governed by partial differential equations (PDEs). The PDEs 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, as written in Section 3.2, up to sixth order. Five arguments are required to instantiate a Sapphire simulation:

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

• 2.

a finite element or mixed finite element defined by UFL,

• 3.

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

• 4.

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

• 5.

initial values as a Firedrake finite element function.

With these defined, a simulation is ready to run.

The solver’s accuracy can be verified using the method of manufactured solutions (MMS), e.g. as described in [34]. Manually deriving MMS source terms creates many opportunities for errors, especially for systems of PDEs involving vector calculus. Sapphire provides a MMS module which uses UFL to automatically derive the source terms and augment 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 its solutions to many different PDEs.

Sapphire is a flexible tool, is publicly available, and is freely licensed. Simulations have been implemented in Sapphire which solve the heat equation, convection-diffusion, steady and unsteady incompressible Navier-Stokes, natural convection, phase-change in enthalpy form without convection, and convection-coupled phase-change in enthalpy form, the last of which is the focus of this article. The version of Sapphire used to produce the results in the current work was archived with a DOI on Zenodo at [24].

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).

4.3.1 MMS procedure

MMS is a general method for verifying PDE solvers. Manufactured solutions do not necessarily need to have physical interpretations [34]. Their purpose purpose rather is to reveal coding errors and to verify the solver. Therefore, a manufactured solution can be any function which is sufficiently differentiable in space and time to exercise the spatial and temporal derivatives in the governing equations. 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.

Substituting a manufactured solution, which is typically 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 which the manufactured solution is an exact solution. 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.

To apply the MMS procedure for verifying a solver’s spatial or temporal orders of accuracy, Sapphire requires the user to

implement the strong form governing equations (10, 11, 12) with UFL,

• 2.

implement a manufactured solution with UFL,

• 3.

and set a single time step size $\Delta t$ and a list of meshes with decreasing cell sizes $h$ for spatial verification, or a single mesh and a list of $\Delta t$ 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,

• 2.

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

• 3.

runs simulations for specified lists of $h$ or $\Delta t$,

• 4.

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

• 5.

reports convergence rates with respect to $h$ or $\Delta t$,

• 6.

and asserts that the reported convergence rate matches the expected order of accuracy within a specified tolerance.

4.3.2 Verification of convection-coupled phase-change

Verification of the current convection-coupled phase-change solver was performed on a unit square domain. Dirichlet boundary conditions were applied to the velocity and temperature. The manufactured solution was

[TABLE]

where $x$ and $y$ are the components of the position $\mathbf{x}$. The definitions of $p^{*}$ and $p_{M}$ were separated so that $p^{*}$ could be chosen freely while $p_{M}$ would maintain zero mean. For the buoyancy term in (11), a classical linear Boussinesq model was used by setting $b=T$. The similarity parameters were set to Re = 20, $\mathrm{Ra}=2.5\times 10^{6}$, $\mathrm{Pr}=7.0$, and $\mathrm{Ste}=0.13$. The material property parameters were set to $\rho_{s}/\rho_{l}=0.92$, $c_{s}/c_{l}=0.50$, and $\kappa_{s}/\kappa_{l}=3.8$. The numerical parameters were set to $\sigma=0.1$, $\tau=10^{-6}$, and $q=4$.

To verify spatial discretization accuracy independently of any time discretization error, the manufactured solution was evaluated at time $t=1$, and no further time dependency was assumed. The resulting problem was solved in a single pseudo-time step. This was repeated on a series of uniformly refined meshes. 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 $\Delta t$. Table 2 shows the observed convergence orders with respect to $h$ and $\Delta t$. As expected, the temperature and velocity solutions converge quadratically. On the other hand, the pressure error shows super-convergence. Given that the pressure solution is not of primary interest in this work, its super-convergence was not further investigated.

5 Validation with benchmark experiments

Two experimental data sets from the literature were considered in order to validate the implemented model. Simulations were run for each, with the goal of accurately modeling the time evolution of the phase interface, i.e. the $T=T_{m}$ isotherm. First, octadecane melting was simulated for comparison to an experiment from [3]. Second, water freezing was simulated for comparison to experiments from [4]. Both cases used a unit square geometric domain with boundary conditions shown in Figure 2. The top and bottom boundaries were assumed to be adiabatic, while the left and right walls were kept respectively at constant hot temperature $T_{h}$ and cold temperature $T_{c}$. Zero velocity boundary conditions were applied on each wall. For both benchmarks, gravity points vertically downward, i.e. $\mathbf{\hat{g}}=(0,-1)^{\mathrm{T}}$. Also in both cases, a uniform triangular mesh with cell edge length $h$ was used. An example mesh with $h=1/6$ is shown in Figure 2, on the right.

Each of the two following subsections is devoted to one of the two cases. In both cases, the subsection includes

a description of the experimental setup. 2. 2.

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

a sensitivity study involving the mesh cell size $h$, time step size $\Delta t$, phase interface regularization parameter $\sigma$, solid velocity relaxation parameter $\tau$, and quadrature degree $q$. With respect to a nominal set of parameters values, each parameter is varied individually and the resulting phase interfaces are compared at a fixed time. 4. 4.

selected simulation results compared to experimental data and to corresponding simulation results in [19].

Where compute times are reported, the simulations ran on a single core of an Intel(R) Core(TM) i5-7500 CPU. The exact versions of Sapphire and Firedrake used to generate all results in this section are archived respectively at [24] and [33].

5.1 Melting octadecane

The melting of octadecane was simulated for comparison to an experiment from [3]. 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 set-up

The speed scale was chosen as $\mathrm{U}=\nu_{l}/X$. The height of the experimental cavity as taken from [3] was used as the length scale, i.e. $\mathrm{X}=0.015\ \mathrm{m}$. A realistic value for the liquid kinematic viscosity, $\nu_{l}=10^{-5}\ \mathrm{m^{2}/s}$ as taken from [17], furthermore yields the characteristic time scale $\mathrm{X}^{2}/\nu_{l}=22.5\ \mathrm{s}$. Based on physical parameters provided in [13] and according to the previously introduced definitions (13), the similarity parameters are Ra = 3.27$\times 10^{5}$, 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 $C=\kappa=1$. Finally, a linear Boussinesq model with constant thermal expansion coefficient is assumed. In this case, the buoyancy function in the nondimensional momentum balance (11) is simply $b=T$.

For the nondimensional temperature corresponding to the melting point, i.e. $T=T_{m}=0$, the regularized liquid volume fraction as defined in (3) is evaluated to be $\phi_{l}=0.5$. This is true for any chosen regularization parameter $\sigma$. 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. $T_{0}=-0.01$. This approach follows earlier work, e.g. [13, 14, 15]. For this particular choice of the initial temperature $T_{0}$, setting $\sigma=0.004$ yields an initial liquid volume fraction $\phi_{l}=0.006$, 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 2, this translates into $T_{h}=1$ and $T_{c}=T_{0}$. 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 $h=1/56$, $\Delta t=1$, $\sigma=0.004$, $\tau=10^{-12}$, and $q=4$. From their nominal values, each parameter was varied independently. The resulting phase interfaces were compared at time $t=79$. The results are shown in Figure 3. The time $t=79$ is also the final simulated time presented in Figure 4.

The mesh cell size, $h$, was varied from a maximum of $h=1/14$ to a minimum of $h=1/112$. The upper part of the phase interface was most sensitive to $h$. Between the smallest values, $h=1/56$ and $h=1/112$, 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, $\Delta t$, was varied from a maximum of $\Delta t=1.975$ to a minimum of $\Delta t=0.5$. For the largest size, $\Delta t=1.975$, there was a significant decrease in the melting near the top wall. For $\Delta t\leq 1$, the result was no longer significantly changed by reducing $\Delta t$. The regularization parameter, $\sigma$, was varied from a maximum of $\sigma=0.008$ to a minimum of $\sigma=0.002$. For the largest value, $\sigma=0.008$, the entire phase interface shifted rightward. This was because, at the initial temperature $T_{0}=-0.01$, the initial liquid volume fraction was $\phi_{l}=0.1056$. 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 $\sigma=0.004$, less than one percent of the material was melted at the initial time, and this effect no longer appeared to dominate. Between the two smallest values, $\sigma=0.004$ and $\sigma=0.002$, the change is barely visible. The solid velocity relaxation parameter, $\tau$, was varied from a maximum of $\tau=10^{-8}$ to a minimum of $\tau=10^{-14}$. For the largest value, $\tau=10^{-8}$, significantly more melting occurred near the top wall. An asymptotic behavior was observed while reducing $\tau$ by orders of magnitude. The change when reducing from $\tau=10^{-12}$ to $\tau=10^{-14}$ was negligible. The quadrature degree, $q$, was varied from a minimum of $2$ to a maximum of $8$. For the smallest value, $q=2$, significantly more melting occurred near the top wall. Between $q=4$ and $q=8$ the change was negligible.

5.1.3 Results and discussion

Figure 4 shows the temperature field, streamlines of the velocity field, and the post-processed position of the phase interface at times $t=40$ and $t=79$. The numerical parameters were $h=1/56$, $\Delta t=1$, $\sigma=0.004$, $\tau=10^{-12}$, and $q=4$. 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 ($t=79$), after a total of seventy-nine time steps ($\Delta t~{}=~{}1$), 2847 total Newton iterations were used which ran altogether in a total of one hour (wall-time). This averages to about thirty-six 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 $\sigma=0.256$. For the first time step, only a single intermediate value, $\sigma=0.016$, was required. By the final time, the intermediate values were ($\sigma$ = 0.256, 0.13, 0.067, 0.0355).

The simulation under-predicts melting with respect to the experiment. Qualitatively, the largest deviation is near the top wall at the end of the simulation. The previously conducted sensitivity study, shown in Figure 3, suggests that the phase interface near the top wall could agree better with the experimental data by reducing the mesh cell size $h$, but no such improvement should be expected near the bottom wall or interior.

It has already been reported in [15] 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 is outside of the scope of the current work. Qualitatively, especially away from the top wall, 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.

It is also possible to better match the experiment by adjusting the regularization parameter $\sigma$ and velocity relaxation factor $\tau$. This also allows a more direct comparison between the current approach and the latest results from [19]. Therefore, another Sapphire simulation was run for the current case, only changing $\sigma$ and $\tau$ from the result in Figure 4. Figure 5 compares the result to [19] and to the experiment. The simulation from [19] used 2900 mesh cells which were adaptively re-meshed at every time step, clustering cells near the phase interface. Sapphire’s simulation used 3136 uniformly distributed mesh cells (i.e. $h=1/56$). To approximate the regularization in [19], $\sigma$ was set to $\sigma=0.00875$. Similar to the approach [19], simulations were run with $10^{-6}\leq\tau\leq 10^{-8}$.

Figure 5 shows the result with $\tau=10^{-8}$. The result compares very well to [19] and to the experiment. A total of 949 Newton iterations were used and the simulation required twenty-five minutes of compute time. A run time of one hour and nine minutes was reported in [19]. Their smaller time step size of $\Delta t=0.1$ could account for their longer run time. Many more computer hardware and software variables would need to be controlled to make a detailed performance comparison. Still, it is remarkable that similar results were found with uniform and adapted meshes using similar number of mesh cells and the same finite elements.

5.2 Freezing water

The freezing of distilled water was simulated for comparison to benchmark experiments in [4], where a cube of liquid water was frozen from one side. For water-ice, 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 greatest density. This is sometimes referred to as the density anomaly of water.

5.2.1 Simulation set-up

The speed scale was chosen as $\mathrm{U}=\nu_{l}/X$. The edge length of the experimental cube’s test section, taken from [4], was used as the length scale, i.e. $\mathrm{X}=0.038\ \mathrm{m}$. A realistic value for the liquid kinematic viscosity, $\nu_{l}=1.0032\times 10^{-6}\ \mathrm{m^{2}/s}$ as taken from [35], furthermore yields the characteristic time scale $\mathrm{X}^{2}/\nu_{l}=1440\mathrm{s}$. Based on physical parameters provided in [13] and according to the previously introduced definitions (13), the similarity parameters are Ra = $2.52\times 10^{6}$, Pr = 6.99, and Ste = 0.125. Phase-dependent thermal conductivity and volumetric heat capacity were included with $\kappa_{s}/\kappa_{l}=3.767$ and $(\rho_{s}c_{s})/(\rho_{l}c_{l})=0.4867$ in (8) and (9). The nonlinear water density model from [13] was used. First proposed in [36], the model accounts for the density anomaly of water by expanding the density around its maximum $999.972\ \mathrm{kg/m^{-3}}$ which occurs at $4.0293\ \mathrm{{}^{\circ}C}$.

The simulation initial values correspond to the “warm start” initial conditions from [4]. For this, the left and right walls were respectively kept constant at hot and cold temperatures $\mathpzc{T}_{h}=10\ ^{\circ}\mathrm{C}$ and $\mathpzc{T}_{c}=0\ ^{\circ}\mathrm{C}$. A steady state convection was reached. For the simulation, this initial temperature range was used for the temperature scale, i.e. $\delta\mathrm{T}=\mathpzc{T}_{h}-\mathpzc{T}_{c}=10\ ^{\circ}\mathrm{C}$. With this scale, the dimensionless boundary temperatures were $T_{h}=1$ and $T_{c}=0$. The steady state problem was solved directly. Therefore, the liquid volume fraction was set to a constant $\phi_{l}=1$ and time derivatives were set to zero. In this case, the momentum (11) and energy (12) equations reduce to

[TABLE]

These equations were solved with the same mixed finite elements and Newton method as from Section 3.

The high Rayleigh number caused the nonlinear solver to diverge. Conveniently, when using the continuation module of Sapphire [24], which implements the procedure from Section 3.4, any scalar parameter can be chosen as the continuation parameter instead of $\sigma$. The Rayleigh number was chosen as the continuation parameter, and bounded search mode was used, bounded between values of zero and Ra. The resulting continuation sequence was 0, Ra/16, Ra/8, Ra/4, Ra/2, and finally Ra.

To begin the time dependent simulation, the cold wall temperature was dropped to $T_{c}=-1$ while the hot wall temperature remained at $T_{h}=1$. This caused freezing to proceed from right to left. Note that the temperature scaling was unchanged. Therefore, the physical temperature of the cold wall was $\mathpzc{T}_{c}=-10\ ^{\circ}\mathrm{C}$, as it was in the experiment.

5.2.2 Sensitivity to numerical parameters

For this sensitivity study, the nominal parameter values were $h=1/56$, $\Delta t=0.2$, $\sigma=0.004$, $\tau=10^{-10}$, and $q=4$. From their nominal values, each parameter was varied independently. The resulting phase interfaces, shown in Figure 7, were compared at time $t=1.6$. The mesh cell size, $h$, was varied from a maximum of $h=1/14$ to a minimum of $h=1/112$. Using the largest size significantly degraded the accuracy. The solution appeared to converge asymptotically as $h$ was reduced. Between the smallest two values, 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 $h$ when using the same values.

The time step size, $\Delta t$, was varied from a maximum of $\Delta t=0.8$ to a minimum of $\Delta t=0.1$. Only the largest time step size significantly affected the result. The regularization parameter, $\sigma$, was varied from a maximum of $\sigma=0.016$ to a minimum of $\sigma=0.001$. The difference between the two smallest values ($\sigma=0.002$ and $\sigma=0.001$) was negligible. The solid velocity relaxation parameter $\tau$ was varied from a maximum of $\tau=10^{-6}$ to a minimum of $\tau=10^{-12}$. For the largest value $\tau=10^{-6}$, significantly more freezing occurred near the center and less freezing occurred near the bottom wall. The solution appeared to be approaching an asymptotic limit. Between the two smallest values of $\tau$ ($\tau=10^{-10}$ and $\tau=10^{-12}$), the change was negligible. The quadrature degree, $q$, was varied between $q=2$, $q=4$, and $q=8$, and the effects were negligible.

Figure 8 further demonstrates the effects of varying $\sigma$ and $\tau$ on a constant coarse mesh. The size of the artificial mushy region (where $0<\phi_{l}<1$) varied with $\sigma$. For the case with largest $\sigma$ and smallest $\tau$, the velocity solution was disturbed so greatly that the convection-cell in the lower-right of the liquid region disappeared entirely. The one-dimensional profiles on the bottom of Figure 8 further demonstrate the effect of the velocity relaxation in the artificial mushy region. For both cases with $\sigma=0.005$, the artificial mushy region was contained by a single cell (though for the $\tau=10^{-6}$ case, a point with $0<\phi_{l}<1$ is visible, because there are two temperature degrees of freedom in each direction of each cell). With $\sigma=0.08$, the artificial mushy region extended over multiple cells. For both cases with $\tau=10^{-12}$, the velocity was approximately zero at the left edge of the artificial mushy region and remained zero throughout the solid. For both cases with $\tau=10^{-6}$, the velocity was substantially above zero at the left edge of the artificial mushy region (though for the $\sigma=0.005$ case, this is harder to see). Furthermore, for the $\tau=10^{-6}$ cases, while the velocity reached approximately zero at the left edge of the solid region, it steadily increased while moving further into the solid region, until the rightmost cell where it again dropped to zero to meet the boundary condition.

5.2.3 Results and discussion

The sensitivity study revealed that large time step sizes $\Delta t$ can be used without significantly changing the resulting phase interface. To most easily compare to experimental results in [4], solutions were obtained at every one-hundred physical seconds, corresponding to a simulated time step size of about $\Delta t=0.0695$. The other numerical parameters were $h=1/112$, $\sigma~{}=~{}0.004$, $\tau=10^{-10}$, and $q=4$. Figure 9 visualizes results at the initial steady state and two later times, $t=0.3475$ (500 seconds) and $t=1.5985$ (2300 seconds). The temperature field is colored and the velocity streamlines are black with thickness proportional to the velocity magnitude. 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.

By the final time, after a total of twenty-three time steps, 1910 Newton iterations were used in total, i.e. about 83 Newton iterations per time step. This includes all iterations from intermediate problems during the continuation procedure sketched in Figure 1. For the first time step, eight intermediate values of $\sigma$ = 0.512, 0.258, 0.131, 0.0675, 0.03575, 0.0278125, 0.019875, and 0.0119375 were used for continuation. At $t=0.417$ an additional intermediate value of $\sigma=0.09925$ was used, which sufficed for the rest of the simulation. As shown in Figure 1, each step in continuation corresponds to solving a nonlinear system with Newton’s method. For any given time step in this simulation, the average number of Newton iterations per continuation step was between eight and nine.

Figure 10 compares the simulation results to a group of experimental runs from [4]. The simulation compares remarkably well to experimental run #1; but that run reportedly used a “cold start”. The simulation instead used a “warm start”, which corresponds to the experimental runs #4 and #5. It was also shown in [4] 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 three-dimensional heat transfer in the side walls was further demonstrated in [10]. 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.

Figure 11 compares a simulation using Sapphire to results presented in [19]. As when comparing to their octadecane melting simulation, the regularization parameter was set to $\sigma=0.00875$ and simulations were ran with $10^{-6}<\tau<10^{-8}$. In the simulation from [19], the mesh was adaptively refined at each time step, resulting in meshes of less than 3000 vertices. That simulation required several days [19] of compute time. Sapphire’s simulation, using 2500 uniformly distributed mesh vertices (i.e. $h=0.02$), required fifteen minutes of compute time (having solved a total of 784 Newton iterations) and yielded a nearly identical result.

6 Additional examples

Two additional examples demonstrate the flexibility of the method and code: a higher Rayleigh number octadecane melting example based on [37] and a gallium (lower Prandtl number) melting example based on [20]. Both simulation setups were mostly similar to the octadecane melting simulation in Section 5.1. The same code was run on the same CPU as in Section 5.1.

6.1 Melting octadecane at a higher Rayleigh number

This example is characterized primarily by the higher Rayleigh number of Ra = $10^{7}$, where Ra = 3.27 $\times 10^{5}$ was used in Section 5.1. The other similarity parameters were Pr = 50 and Ste = 0.1. The length scale was X = 0.1 meters and the time scale was $\mathrm{X}^{2}/\nu_{l}=1,000$ seconds. A rectangular domain was used with $L_{x}=0.5$. The numerical parameters were $h=1/80$, $\Delta t=0.05$, $\sigma=0.005$, $\tau=10^{-10}$, and $q=4$.

By the final time, i.e. after one-hundred time steps, a total of 2900 Newton iterations were used including all continuation steps. Up to four intermediate regularizations were required for continuation at each time step. The simulation required one hour and three minutes of compute time. The result is shown in Figure 12 and nearly matches the accepted result from [37] that predicted the least melting.

6.2 Melting gallium

This example is characterized primarily by a smaller Prandtl number of Pr = 0.0216. The other similarity parameters were Ra = $7\times 10^{5}$ and Ste = 0.046. Unlike the other simulations presented in this work, the velocity scale was chosen as $\mathrm{U}=\alpha_{l}/\mathrm{X}$ and therefore the Reynolds number appearing in the dimensionless governing equations was Re = 1/Pr. The nondimensional temperature initial value and cold wall boundary condition was $T_{c}=-0.1546$. The length scale was X = 0.0635 meters and the time scale was $\mathrm{X}^{2}/\alpha_{l}=292.9$ seconds. A rectangular domain was used with $L_{x}=0.25$. Using $h=1/160$, and again a uniform mesh, there were a total of 12,800 mesh cells. The other numerical parameters were $\Delta t=0.001$, $\sigma=0.01$, $\tau=10^{-10}$, and $q=4$.

By the final time, i.e. after 409 time steps, a total of 5700 Newton iterations were used including all continuation steps. Up to five intermediate regularizations were required for continuation at each time step. The simulation required three hours and fifty minutes of compute time. Figure 13 shows the result. Compared to the results from [20], where adapted meshes with an average of 12,000 mesh cells were used on a rectangular domain with $L_{x}=0.5$, the largest discrepancy is with the upper convection cells at time $t=0.29$. An additional convection cell developed in the [20] result. Overall, the results compare remarkably well.

7 Conclusions

An enthalpy method was used to simulate convection-coupled isothermal phase-change 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

presented a new continuation procedure that reliably converged Newton’s method for the tested 2D isothermal phase-change benchmarks. Without this procedure, previous literature such as [13] and [14] relied on carefully selected time step sizes, specifying a one-way melting or solidification process, initializing the new phase in part of the domain, multiple approaches to adaptive mesh refinement, and manually refining the initial mesh, in order to achieve convergence of Newton’s method. The key idea of the new continuation procedure is to recognize that the regularity of the nonlinear problem’s solution space is dominated by a global scalar parameter, the regularization parameter $\sigma$. The procedure highlights a path forward for regularizing this class of nonlinear problems.

• 2.

shared the open source code in a new Python packaged called Sapphire [24], using the finite element library Firedrake [23]. Sapphire’s test suite covers many of this paper’s results. These tests can be used as scripts to reproduce the results. The tests are separated into verification tests and validation tests. The validation tests include both benchmark problems from this paper, using coarse meshes and large time steps so that they run quickly. The versions of Sapphire and Firedrake were documented with DOIs on Zenodo, respectively at [24] and [33].

• 3.

verified the implementation via a formal convergence study using the method of manufactured solutions. Second order spatial accuracy was verified in 2D for the velocity and temperature. The pressure’s spatial discretization error showed super-convergence, which was not further explored in this contribution. Second order temporal accuracy was verified for all solution components which were discretized in time, i.e. the velocity and the temperature.

• 4.

presented benchmark simulations for melting octadecane and freezing water in square cavities. For each, sensitivities were studied with respect to the mesh cell size $h$, time step size $\Delta t$, regularization parameter $\sigma$, solid velocity relaxation parameter $\tau$, and quadrature degree $q$. For the freezing water case, results were obtained which were not largely sensitive to the five numerical parameters $h$, $\Delta t$, $\sigma$, $\tau$, or $q$. For the octadecane melting case, the refined result was still largely sensitive to $h$. Further reducing $h$ was impractical, because it would increase the number of degrees of freedom in the linear system quadratically, which is particularly prohibitive for the direct linear system solver. Reducing $h$ could affect sensitivity in the other parameters, especially $\sigma$.

• 5.

compared the benchmark simulation results to experimental data sets from the literature, namely [3] for melting octadecane and [4] for freezing water. For the melting octadecane case, the simulation under-predicts melting with respect to the experiment. The largest discrepancy is near the top wall at the end of the simulation. Reducing the mesh cell size $h$ should reduce the discrepancy near the top wall, but no such improvement should be expected near the bottom wall or interior. For the water freezing case, multiple experimental runs were considered. The runs varied with respect to initial conditions and thermal regulation of the test section walls. While the simulation compared well to one of the runs, that run 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.

• 6.

compared the benchmark simulation results to those from [19] which used adaptive mesh refinement. The results compare favorably, even when using uniform meshes with a similar number of cells.

• 7.

presented two additional example simulations, one for melting octadecane with a higher Rayleigh number, and another for melting gallium. In both cases, results compare favorably with those in the literature.

The application scope of the Sapphire software [24] is potentially much larger than demonstrated in the current work. Having used the finite element method, application to more complicated geometries should be straightforward. By using UFL [25] for symbolically defining the mathematical model and Firedrake [23] for automating most of the implementation, the model and implementation are easily modifiable and extensible. The continuation procedure has been an effective tool for reliably solving the nonlinear problems in this work for a large range of parameters. Still, there should exist a more general, less ad hoc, approach to regularizing this class of nonlinear problems. Finding such an approach could further decrease computational costs and further improve robustness. There are other 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. Implementing the iterative method and preconditioner should be straightforward using Firedrake. Successfully applying an iterative solver would make it practical to solve problems with three-dimensional geometries, 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).