A finite-element toolbox for the simulation of solid–liquid phase-change systems with natural convection

Abstract

We present and distribute a new numerical system using classical finite elements with mesh adaptivity for computing two-dimensional liquid–solid phase-change systems involving natural convection. The programs are written as a toolbox for FreeFem++ (www3.freefem.org), a free finite-element software available for all existing operating systems. The code implements a single domain approach. The same set of equations is solved in both liquid and solid phases: the incompressible Navier–Stokes equations with Boussinesq approximation for thermal effects. This model describes naturally the evolution of the liquid flow which is dominated by convection effects. To make it valid also in the solid phase, a Carman–Kozeny-type penalty term is added to the momentum equations. The penalty term brings progressively (through an artificial mushy region) the velocity to zero into the solid. The energy equation is also modified to be valid in both phases using an enthalpy (temperature-transform) model introducing a regularized latent-heat term. Model equations are discretized using Galerkin triangular finite elements. Piecewise quadratic (P2) finite-elements are used for the velocity and piecewise linear (P1) for the pressure. For the temperature both P2 and P1 discretizations are possible. The coupled system of equations is integrated in time using a second-order Gear scheme. Non-linearities are treated implicitly and the resulting discrete equations are solved using a Newton algorithm. An efficient mesh adaptivity algorithm using metrics control is used to adapt the mesh every time step. This allows us to accurately capture multiple solid–liquid interfaces present in the domain, the boundary-layer structure at the walls and the unsteady convection cells in the liquid. We present several validations of the toolbox, by simulating benchmark cases of increasing difficulty: natural convection of air, natural convection of water, melting of a phase-change material, a melting-solidification cycle, and, finally, a water freezing case. Other similar cases could be easily simulated with this toolbox, since the code structure is extremely versatile and the syntax very close to the mathematical formulation of the model.

Programm summary

Program Title: PCM-Toolbox-2D

Program Files doi: http://dx.doi.org/10.17632/phby62rhgv.1

Licensing provisions: Apache License, 2.0

Programming language: FreeFem$++$(free software, www3.freefem.org)

Nature of problem: The software computes 2D configurations of liquid–solid phase-change problems with convection in the liquid phase. Natural convection, melting and solidification processes are illustrated in the paper. The software can be easily modified to take into account different related physical models.

Solution method: We use a single domain approach, solving the incompressible Navier–Stokes equations with Boussinesq approximation in both liquid and solid phases. A Carman–Kozeny-type penalty term is added to the momentum equations to bring the velocity to zero into the solid phase. An enthalpy model is used in the energy equation to take into account the phase change. Discontinuous variables (latent heat, material properties) are regularized through an intermediate (mushy) region. Space discretization is based on Galerkin triangular finite elements. Piecewise quadratic (P2) finite-elements are used for the velocity and piecewise linear (P1) for the pressure. For the temperature both P2 and P1 discretizations are possible. A second order Gear scheme is used for the time integration of the coupled system of equations. Non-linear terms are treated implicitly and the resulting discrete equations are solved using a Newton algorithm. A mesh adaptivity algorithm is implemented to reduce the computational time and increase the local space accuracy when (multiple) interfaces are present.

Introduction

Solid–liquid phase-change problems are encountered in numerous practical applications, ranging from metal casting and thermal energy storage (phase-change materials) to food freezing and cryosurgery. Melting and solidification are also fundamental processes in geophysical problems, such as Earth’s mantle formation, lava lakes or magma chambers. In many of these problems, temperature gradients induce buoyancy forces in the liquid (melted) phase generating a significant convective flow. Consequently, the appropriate mathematical description of the liquid phase is the usual model for the natural convection flow: the incompressible Navier–Stokes system of equations with Boussinesq approximation for thermal (buoyancy) effects. In this model, the energy conservation equation is written as a convection–diffusion equation for the temperature. In the solid phase, conduction is the main phenomenon and the appropriate model is the classical heat equation. The main modelling difficulty is to link these two models by taking into account the separation of the two phases by a sharp interface, across which thermodynamic properties are discontinuous. We offer below a short description of the two main approaches suggested in the literature to deal with this problem. For a comprehensive review of models for phase-change problems with convection, see [1]. Note that a different category of models was recently suggested in the literature, based on the Lattice Boltzmann Method [2], [3]. Such methods based on non-deterministic models are not discussed in this introduction.

A first modelling approach, usually referred to as the multi-domain (or deforming-grid) method, is based on the classical Stefan two-phase model. Solid and liquid domains are separated and the corresponding conservation equations are solved in each domain. Boundary conditions at the interface between domains are obtained by imposing the Stefan condition (balance of heat fluxes at the interface). The position of the solid–liquid interface is tracked and moved explicitly using either front tracking or front fixing methods. The former method uses deforming grids to reconstruct the interface, while the latter is based on a time-depending coordinate transform, mapping the physical domain into a fixed computational domain. For a detailed description of these methods, see e.g. [4], [5], [6], [7]. The main drawback of deforming-grid methods is their algorithmic complexity, which makes difficult to accurately capture solid–liquid interfaces of complicated shape or structure (e.g. with mushy regions between solid and liquid phases). Configurations with multiple interacting interfaces (see the solidification of a phase-change material presented in this paper) are also difficult to simulate with these methods (see also [8]).

The second modelling approach avoids to impose explicitly the Stefan condition at the solid–liquid interface and therefore uses a single-domain (or fixed-grid) model. The same system of equations is solved in both liquid and solid phases. The energy balance at the interface is implicitly taken into account by the model. Consequently, the position of the interface is obtained a posteriori by post-processing the computed temperature field. Phase-field methods [9] and enthalpy methods [10], [11] are the most commonly used single-domain models. In phase-field methods, a supplementary partial differential equation for the evolution of the order parameter (a continuous variable taking the value 0 in the solid and 1 in the liquid) has to be solved, coupled with the conservation laws [12]. This new equation is model dependent and its numerical solution could lead to diffuse solid–liquid interfaces. For recent contributions in this area, see [13], [14], [15]. We focus below on enthalpy methods, which are the most widely used single-domain models due to their algorithmic simplicity.

The main idea behind enthalpy models is to formulate the energy conservation law in terms of enthalpy and temperature and include latent heat effects in the definition of the enthalpy. The obtained equation applies to both liquid and solid phases and implicitly takes into account the separation of the phases. Another advantage of enthalpy methods, when compared to previously described models, is to remove the limitation of the phase-change occurring at a fixed temperature. The presence of mushy regions can be easily modelled with these methods. Two types of formulations of enthalpy methods exist in the literature, depending on the main variable used to solve the energy equation: enthalpy or temperature-based methods. In enthalpy-based formulations the main variable is the enthalpy [16], [17], [18]. Temperature is computed from the temperature–enthalpy coupling model. An iterative loop is necessary to solve the energy equation, formulated with both enthalpy and temperature variables. For a review of different iterative techniques to solve the energy equation, see [19]. A second variety of enthalpy-based formulations consists in rewriting the energy equation with enthalpy terms only [20], [21]. In temperature-based formulations, the energy equation is formulated in terms of temperature only. The latent heat is treated either by deriving an apparent heat capacity coefficient to define the total enthalpy [22], [23], [24] or by introducing a source term in the energy equation [19], [25]. Advantages and drawbacks of each approach are discussed in detail in [26].

Single-domain methods are very appealing for numerical implementations. The same system of equations is solved in the entire computational domain, making possible algorithmic or computer-architecture optimizations. If enthalpy models offer an elegant solution to deal with the same energy conservation equation in both phases, a last modelling problem has to be solved. It concerns the extension of the Navier–Stokes–Boussinesq equations from the liquid to the solid phase. Different techniques to bring the velocity to zero in the solid region were suggested. The most straightforward is the switch-off technique, which decouples solid and liquid computational points and overwrites the value of the velocity by setting it to zero in the solid region. Different implementations of this technique with finite-volume methods are presented in [27], [28]. In variable viscosity techniques [29], [30], [31], the fluid viscosity depends on the temperature and is artificially increased to huge values in the solid regions through a regularization or mushy zone. To avoid blow-up or numerical inconsistencies, the large gradients of viscosity must be correctly resolved in the mushy region. This is naturally achieved in finite-element methods with dynamical mesh adaptivity [32], while in finite-volume methods with fixed grids, the time step has to be adapted to the space resolution [27]. Versions of the variable viscosity approach were theoretically studied by Aldbaissy et al. [33]; Woodfield et al. [34]. A new formulation of the phase-change problem using as variables the pseudo-stress, strain rate and velocity for the Navier–Stokes–Brinkman equations was recently suggested by Alvarez et al. [35].

A third technique used to ensure a zero velocity field in the solid phase is the so-called enthalpy-porosity model [36]. A penalization source term is introduced in the momentum equation to allow the switch from the full Navier–Stokes equations in the liquid phase to a Darcy equation for porous media. The mushy region is thus regarded as a very dense porous medium that sharply brings the velocity to zero in the solid region. The expression of the penalization source term generally follows the Carman–Kozeny model for the permeability of a porous medium [21], [37], [38], but other mathematically equivalent expressions were suggested [15], [39]. Different formulations and implementations of the enthalpy-porosity model are presented in [8], [40], [41].

Concerning the space discretization of these models, finite difference (FD) or finite volume (FV) methods are generally used in the literature. When single-mesh models are used, the general strategy to capture the solid–liquid interface is to dramatically increase the mesh resolution in the whole domain. This results in a considerable increase of the computational time, even for two dimensional cases. Finite element (FE) methods offer the possibility to dynamically refine the mesh only in the regions of the domain where sharp phenomena take place (e.g. solid–liquid interface, recirculation zones). Different FE approaches were suggested, from enthalpy-type methods (e.g. [42]) to front-tracking methods (e.g. [43]). Recently, adaptive FE methods were suggested for classical two-phase Stefan problem and phase-change systems with convection. In [44] a hierarchical error estimator was used for 2D Stefan problems (without convection). A different mesh adaptivity method, based on the definition of edge length from a solution dependent metric, was used in [38] to deal with the same Stefan problems in 3D and also with 2D phase-change systems with convection. The time-dependent mesh adaptivity strategy based on metrics control suggested by Danaila et al. [32] proved very effective in refining the mesh near walls and interfaces, when simulating melting and solidification phase-change systems (see also [45]). Zimmerman and Kowalski [46] implemented the variable viscosity model suggested by Danaila et al. [32] in a different finite-element framework using FEniCS and an AMR (Adaptive Mesh Refinement) technique based on a dual-weighted residual method. In a very recent contribution, Belhamadia et al. [47] derived a time-dependent adaptive remeshing method for phase-change problem with convection based on an error estimator applicable to second or higher order variables.

However, to the best of our knowledge, no adaptive finite-element programs exist in the CPC Program Library for phase-change problems. The purpose of this paper is therefore to distribute a finite-element toolbox to solve two-dimensional solid–liquid phase-change problems.

The present toolbox is based on a single-domain enthalpy-porosity model for solid–liquid phase change problems with convection. For the energy conservation equation, a temperature-based formulation takes into account the latent heat by introducing a discontinuous source term. For the mass and momentum conservation equations, we solve in the entire domain the incompressible Navier–Stokes equations with Boussinesq approximation for buoyancy effects. To bring the velocity to zero in the solid phase, we introduce in the momentum equation a penalty term following the Carman–Kozeny model. The coupled system of momentum and energy equations is integrated in time using a second-order Gear scheme. All the terms are treated implicitly and the resulting discretized equations are solved using a Newton method [32]. The advantage of this formulation is to permit a straightforward implementation of different types of non-linearities. For the space discretization we use Taylor–Hood triangular finite elements, i.e. P${}_2$ for the velocity and temperature and P${}_1$ for the pressure. Temperature is discretized using P${}_2$ or P${}_1$ finite elements. The choice of using Taylor–Hood finite elements for the fluid flow and P${}_2$ for the temperature was also made in [34], [47]. It offers second order accuracy and (inf–sup) stability of the space discretization of the fluid flow. Other possible choices for the finite-element discretization could be tested for the implementation of the algorithm presented in this paper. For instance, the mini-element (see [48], [49]), offering a global linear convergence, could be an interesting alternative to replace the Taylor–Hood element and thus optimize the computational time. This possibility could be explored with the present toolbox, since the mini-element exists in FreeFem$++$.

Discontinuous variables (latent heat, thermal diffusivity, etc.) at the solid–liquid interface are regularized through an intermediate artificial mushy region. Single domain methods require a refined mesh near the interface, where large enthalpy gradients have to be correctly resolved. An optimized dynamical mesh adaptivity algorithm allows us to adapt the mesh every time step and thus accurately capture the evolution of the interface. Mesh adaptivity feature of the toolbox offers the possibility to deal with complicated phase-change cases, involving multiple solid–liquid interfaces. There are two main novelties in the present numerical approach, when compared to Danaila et al. [32]: (i) we use the Carman–Kozeny model to bring the velocity to zero inside the solid phase, instead of a viscosity penalty method (imposing a large value of the viscosity in the solid); (ii) we increase the time accuracy of the algorithm by replacing the first-order Euler scheme with the second-order Gear (BDF2) scheme (see also [38]).

The programs were built and organized as a toolbox forFreeFem$++$ [50], [51], which is a free software (under LGPL license). FreeFem$++$¹ offers a large variety of triangular finite elements (linear and quadratic Lagrangian elements, discontinuous P${}_1$, Raviart–Thomas elements, etc.) to solve partial differential equations. It is an integrated product with its own high level programming language and a syntax close to mathematical formulations, making the implementation of numerical algorithms very easy. Among the features making FreeFem$++$ an easy-to-use and highly adaptive software we recall the advanced automatic mesh generator, mesh adaptation, problem description by its variational formulation, automatic interpolation of data, colour display on line, postscript printouts, etc. The FreeFem$++$ programming framework offers the advantage to hide all technical issues related to the implementation of the finite element method. It becomes then easy to use the present toolbox to code new numerical algorithms for similar problems with phase-change.

The paper is organized as follows. Section 2 introduces the governing equations and the regularization technique used in the enthalpy-porosity model. Section 3 presents the Newton algorithm for the Navier–Stokes–Boussinesq system of equations. The mesh adaptivity technique is also discussed in this section. A description of the programs is given in Section 4. The accuracy of the numerical method is tested in Section 5 using manufactured solutions. Finally, Section 6 is devoted to extensive numerical validations of the method. The robustness of the algorithm is demonstrated by comparing our results with well defined classical benchmarks. The capabilities of the toolbox to deal with complex geometries are also illustrated. The main features of the software and possible extensions are summarized in Section 7.