Mathematical Modeling and Associated Numerical Simulation of Fusion/Solidification Front Evolution in the Context of Severe Accident of Nuclear Power Engineering

Abstract

Considering transient processes where liquid/solid phase change occurs, this paper focuses on the associated modeling and numerical treatment in the frame of “Computational Fluid Dynamics” simulations. While being of importance in many industrial applications involving solidification and melting of mixed materials, including power and manufacturing engineering, the first application of this work pertains to the analysis of severe accidents in a nuclear reactor. Indeed, in this context, the molten core materials (a.k.a. corium) can form a high-temperature multiphase liquid pool at the boundary of which fusion and solidification phenomena are of prime importance. In this context, even if materials at play are treated as pure components, it is mandatory to distinguish two different phase change temperatures with a solid fusion temperature and a liquid solidification temperature. Accordingly, in the frame of a sharp interface representation, the paper introduces non-classical heterogeneous conditions at the liquid/solid boundary in such a way that both moving interface (through Stefan conditions associated with fusion or solidification) and static interface (imposing heat flux continuity) are supported at the same time on different spatial locations along this boundary. Within a monolithic resolution of Navier–Stokes and heat conduction equations, this interface is explicitly tracked with combined Front-Tracking and VOF methods. In order to ensure zero velocity in the solid phase, an Immersed Boundary Method and a direct forcing penalization are also introduced. The main relevant features of this combination of numerical methods are discussed along with their implementation in the TrioCFD code taking advantage of the pre-existing code capabilities. Numerical simulations including both verification tests and a case of interest for our industrial application are reported and demonstrate the applicability of the proposed triptych model+methods+code to treat such problems. The numerical tools and the simulation code developed in this work could be used not only in the several accident context but also to simulate melting, solidification and fusion processes occurring in aerodynamics, hypersonic reentry vehicles and laser applications to cite but a few.

1. Context and Introduction

A vast amount of physical processes involve solidification and melting of mixed materials, including power and manufacturing engineering. For instance, in the context of nuclear power engineering, melting processes are of prime interest during a severe accident progression as a consequence of the insufficient nuclear reactor core cooling. While not restricted to this context, we shall nonetheless employ it as a toy context for our subsequent developments.

Such ‘phase change’, from solid to liquid and vice-versa, implies drastic thermo-mechanical and chemical modifications of the materials under scrutiny. The transition area, also referred to as the ‘interface’, is of paramount importance for the appropriate description of such rapid modifications. In this work, we ought to improve the mathematical modeling and numerical simulations of such phase changes. The general physical model involves an incompressible hot liquid and a solid material exchanging mass and heat through their common interface. The motion of this interface is determined by these phase changes, i.e., the melting of the solid or the solidification of the liquid. Meanwhile, the phase changes are induced by the temperature gradients with respect to the material characteristics. The temperature gradients are modified by the liquid convection, the cooling/heating from the external world, and the transfer of mass and heat during the phase changes. The corresponding mathematical modeling must involve those coupled effects. As a consequence, the model must pay a special attention at capturing and tracking the interface region in-between the liquid and solid phases, and the different materials involved. Moreover, in the context of nuclear reactor core meltdown, the material under consideration is the so-called ‘corium’, a complex mixture of nuclear fuel, fission products, control rods, structural materials of the reactor, products of their chemical reaction with air, water and steam, etc. This liquid corium is generally contained within a metallic container or vessel, the integrity of which is to be ensured. The liquid corium pool in contact with the cold vessel can solidify, creating a protective ‘solid corium crust’ in-between the vessel and the liquid. On the contrary, large convection phenomena or new material addition may also change the thermal condition in contact with the crust which may melt back into the corium liquid, removing, as such, the thermal protection of the vessel. The determination of the transient conditions, in order to maintain the integrity of the container, must resort mainly to numerical simulation because such extreme situations are difficult to reproduce experimentally at large scale. As such, an appropriate physical model supplemented with robust and efficient numerical simulation tools are needed.

The simplified model of Partial Differential Equations (PDEs) under consideration involves a conducting domain (liquid and solid phases and the container) and the incompressible liquid phase coupled through a moving interface between the liquid and solid phases The conduction in the domain is solved by the heat equation subject to thermal boundary conditions on a fixed grid relying on a Ghost Fluid method [1]. The liquid phase is modeled by the incompressible Navier–Stokes equations subject to gravity within a closed solid pool with only one open boundary on the top, from which, new material can be supplemented. The Navier–Stokes equations are solved by a projection/correction method on a fixed staggered Cartesian mesh [2,3] on the whole domain. The interface is the location where the phase changes may occur, i.e., solidification of the liquid or fusion of the solid. In our context, this moving phase boundary problem reduces to the well known Stefan problem [4]. In the Stefan problem, the moving phases interface and the associated conditions are non-linear functions of the unknown location of the interface changing in time. The source of non-linearity lies in the phase change itself, and, the temperature dependent properties of the material. Moreover the solids are considered as motionless, and, they only can gain or loose mass and energy during the phase changes.

In this work, an explicit tracking of this sharp interface is considered. Indeed, in our case of application, this interface remains continuous as the phase change does not involve complicated phase configurations (e.g., with a dispersed phase embedded in a continuous one), in such a way that a front tracking type modeling [5] is readily applicable. In the frame of multiphase corium simulation, this choice is to be distinguished from:

• an enthalpy-based numerical method [6] that is the most common approach reported in the literature for simulating the solidification at the interface of an homogeneous liquid corium pool (see, for instance, [7]). The work reported in this paper aims at providing an alternative and a priori more accurate approach with an explicit tracking of such a liquid/solid interface;
• the diffuse interface approach followed in [8] where liquid phase stratification involving droplet detachment and coalescence processes in a two-phase corium pool is simulated. This work and the current one are complementary as they pursue the same goal: providing models and numerical methods for simulating the thermalhydraulic behaviour of a corium pool composed of different chemically reactive liquid phases with, possibly, liquid/solid phase change at their boundaries.

This paper focuses on several important characteristics of the numerical methods at play that are related to important features of corium pool simulation. First of all, the complexity of the corium mixture in terms of phase segregation for a multicomponent materials requires, in principle, evolution equations for the phase composition so that, phase thermophysical properties can be evaluated based on the local composition. At this stage of our R&D effort, in the restricted framework of a single liquid and a single solid phase, a simplified assumption is considered where both liquid and solid materials are treated as pure materials with different (but fixed in time) properties. As such, contrarily to the classical Stefan problem, two different phase change temperatures need to be defined: a solidification (resp. fusion) temperature $T_{\mathrm{sol}}$ (resp. $T_{\mathrm{fus}}$) for the liquid (resp. solid) with $T_{\mathrm{sol}}<T_{\mathrm{fus}}$. These temperatures are determined according to the composition of the solid and liquid materials. In between these temperatures, the liquid and solid materials are in pure ‘conduction’ phase: they are neither gaining nor loosing mass nor energy, they simply conduct heat. This conduction situation has to be taken into account in the modeling and the numerical method design, along, obviously with the solidification and fusion situations.

Secondly, our choice of a monolithic resolution of Navier–Stokes and heat conduction implies that some velocity will be inexorably but artificially be created in the solids. We will rely on an Immersed Boundary Method and a direct forcing penalization method to cure such a situation [9,10,11]. Third, the implementation within the TrioCFD platform [12] demands the adaptation of the previous modifications, their validation and verification against analytical or known solutions. These are the main goals of this work.

This paper is organized as follows. After this introduction, we present the context of modeling, the governing equations, focusing on the conditions at the phase front. Then, the numerical methods employed to solve the coupled system of non-linear PDEs are presented in Section 3. A validation and verification test suite is presented in the but-last section, where the triptych model+methods+code is tested on academical or more advanced problems. At last, conclusions and perspectives are drawn.

2. Modeling, Governing Equations

2.1. General Context of Modeling

The context of this work is related to the modeling of the coupling of a hot liquid material and a cold solid material in contact through an interface. The main purpose of our modeling is to simulate the phase change (solidification and fusion) of these materials considering a sharp moving interface in-between the phases. The phases then exchange heat and mass. Specific to our target materials is the existence of an intermediate situation which occurs when the interface temperature is in-between the material prescribed fusion and solidification temperatures. Then no phase change must occur and the materials only conduct heat. This difference between fusion and solidification temperatures occurs usually when the material compositions of the two phases are drastically different. As already mentioned, the liquid material is embedded into a cold container opened from its top and cooled down elsewhere. Via its contact with the rest of the world, some new material mass may be added and heat transfer (cooling) usually occurs.

2.2. Governing Equations

The system of PDEs modeling the physical phenomena is composed of a system defined in a spatial domain $\mathsf{\Omega }$ which may have two phases called “liquid” and “solid”, which are ideally characterized by two moving 3D domains denoted by $\mathsf{\Omega }\left(t\right)={\mathsf{\Omega }}_l\left(t\right)\bigcup {\mathsf{\Omega }}_s\left(t\right)$ and separated by an interface $\mathsf{\Gamma }\left(t\right)=\partial {\mathsf{\Omega }}_s\left(t\right)\cap \partial {\mathsf{\Omega }}_l\left(t\right)$ also referred to as the phase front. The outward pointing unit normal to the boundary $\partial \mathsf{\Omega }$ is denoted by $\mathit{\nu }$. Moreover ${\mathit{n}}^{\mathsf{\Gamma }}$ represents the unit normal to $\mathsf{\Gamma }\left(t\right)$ arbitrarily oriented from the liquid domain to the solid one. Each phase is further identified by the color function $\chi$ of value 1 in the solid while 0 in the liquid. The normal vector to $\mathsf{\Gamma }\left(t\right)$ is related to the color function by:

$$\begin{array}{cc}\hfill {\mathit{n}}^{\mathsf{\Gamma }}{\delta }^{\mathsf{\Gamma }}& =\frac{\nabla \chi }{\parallel \nabla \chi \parallel }.\hfill \end{array}$$

(1)

where ${\delta }^{\mathsf{\Gamma }}$ is the Dirac delta function at the interface. Physical properties (conductivity, density…) within each part of the phase are assumed constant, and each phase is incompressible. For any variable $\varphi$ defined in $\mathsf{\Omega }\left(t\right)={\mathsf{\Omega }}_l\left(t\right)\cup {\mathsf{\Omega }}_s\left(t\right)$ the following relation holds:

$$\begin{array}{cc}\hfill \varphi (\mathit{x},t)& =\chi (\mathit{x},t){\varphi }_s(\mathit{x},t)+(1-\chi (\mathit{x},t)){\varphi }_l(\mathit{x},t),\hfill \end{array}$$

(2)

where ${\varphi }_k$ is the value of function $\varphi$ in each phase ${\mathsf{\Omega }}_k$ for $k=l,s$.

The general modeling is ruled by a set of conservation laws namely for the mass, momentum, and energy as function of the position $\mathit{x}$ and time t under gravity $\mathit{g}$. This dependency is always implicitly assumed and omitted for now on. The materials are supposed to be conductive and possibly in motion for the fluid part. Hence we introduce the velocity of the material, $\mathit{u}$, and its temperature, T, defined in $\mathsf{\Omega }$. Several material dependent parameters such as the density $\rho$, latent heat, conductivity, etc. will be later introduced as needed. The material temperature is computed by the resolution of the heat equation subject to boundary temperatures. The material velocity is zero in the solid parts and must remain so. For the liquid phase, we make the hypothesis that the conservation laws are written under the Boussinesq approximation in terms of the liquid velocity $\mathit{u}$ and temperature variation $\delta T$ as the incompressible Navier–Stokes (NS) equations for all $\mathit{x}\in {\mathsf{\Omega }}_l\left(t\right)$ and time $t\in {\mathbb{R}}^{+}$. Nonetheless the possible phase changes do generate mass and heat exchanges between liquid and solid parts as their interface moves with an unknown phase front velocity $\mathit{v}$. We model the phase front velocity by the so-called Stefan condition which is ruled by the difference between heat fluxes given by Fourier’s law.

Let us briefly describe our monolithic formulation of the Navier–Stokes equations coupled with heat transfer. Considering such a monolithic formulation, the unknowns ($\mathit{u}$, P, $\delta T$, $\chi$) are defined in the whole domain $\mathsf{\Omega }$ and follow the PDEs:

$$\begin{array}{ccc}\hfill \nabla ·\mathit{u}& =-\dot{m}\left(\frac{1}{{\rho }_l}-\frac{1}{{\rho }_s}\right){\delta }^{\mathsf{\Gamma }},\hfill & \hfill \mathrm{in}\mathsf{\Omega }\left(t\right),\end{array}$$

(3)

$$\begin{array}{ccc}\hfill \rho \frac{\partial \mathit{u}}{\partial t}+\rho \mathit{u}·\nabla \mathit{u}& =-\nabla P+\rho \mathit{g}+\mathcal{D}-{\rho }_l\beta (T-T_0)(1-\chi )\mathit{g},\hfill & \hfill \mathrm{in}\mathsf{\Omega }\left(t\right),\end{array}$$

(4)

$$\begin{array}{ccc}\hfill \frac{\partial \rho C_p\delta T}{\partial t}+\nabla ·\left(\rho C_p\mathit{u}\delta T\right)& =\nabla ·(\lambda \nabla T)+\dot{m}\mathcal{H}{\delta }^{\mathsf{\Gamma }},\hfill & \hfill \mathrm{in}\mathsf{\Omega }\left(t\right),\end{array}$$

(5)

$$\begin{array}{ccc}\hfill \frac{\partial \chi }{\partial t}+\mathit{v}·\nabla \chi & =0,\hfill & \hfill \mathrm{in}\mathsf{\Omega }\left(t\right),\end{array}$$

(6)

where $\mathcal{D}$ is the diffusion term, $\mathit{u}$ is the fluid velocity, ${\rho }_l$ (resp. ${\rho }_s$) the constant liquid density (resp. solid density), ${\delta }^{\mathsf{\Gamma }}$ the Dirac delta function at the interface, P the pressure and $\mathit{g}={(0,-g)}^t$ the gravity vector, where $g=9.81$. $T_0$ corresponds to a reference temperature and $\beta$ refers to the thermal expansion coefficient. Recall that $\rho$ refers to the relation (2), i.e., $\rho =\chi {\rho }_s+(1-\chi ){\rho }_l$. As already mentioned $\delta T$ is the temperature difference between the material temperature T and the reference temperature $T_0$, that is $\delta T=T-T_0$. At the interface between the phases, represented by the delta Dirac function ${\delta }^{\mathsf{\Gamma }}$, we define $\dot{m}$ the mass transfer rate due to the phase change and $\mathit{v}$ the corresponding interface velocity. $\mathcal{H}$ corresponds to the latent heat (${\mathcal{H}}^{\mathrm{fus}}$ in case of fusion or ${\mathcal{H}}^{\mathrm{sol}}$ for solidification). The material dependent parameters are $C_p$ the thermal capacity, $\rho$ the density, $\mu$ the dynamic viscosity and $\lambda$ the conductivity. The classical diffusive term for Navier–Stokes equations recasts into:

$$\begin{array}{cc}\hfill \mathcal{D}& =\nabla ·\left(\mu (\nabla \mathit{u}+{\nabla }^T\mathit{u})\right).\hfill \end{array}$$

(7)

System (3)–(6) models several physical phenomena. First, Equations (3) and (4) without source term correspond to the incompressible Navier–Stokes equations, while the source term describes the change in volume due to the possible phase change and mass transfer. Equation (5) models the heat equation and the last term, $\dot{m}\mathcal{H}{\delta }^{\mathsf{\Gamma }}$, corresponds to an extra heat exchange term associated to the mass transfer. At last Equation (6) models the motion of the color function, representing the phase front motion at velocity $\mathit{v}$. This velocity field only exists at the phase front $\mathsf{\Gamma }$ and we extend it to zero elsewhere. In Equation (4), the zero velocity condition in the solid is not continuously set. This condition will be ensured numerically and detailed in Section 3.4.

The Boundary Conditions (BCs) in $\partial \mathsf{\Omega }$ can be of Dirichlet or Neumann types, that is

$$\begin{array}{cccc}\hfill \mathit{u}& ={\mathit{u}}^{\partial {\mathsf{\Omega }}^D},\hfill & \hfill T& =T^{\partial {\mathsf{\Omega }}^D},\hfill \end{array}$$

(8)

$$\begin{array}{cccc}\hfill \frac{\partial \mathit{u}}{\partial \mathit{n}}& =\nabla u^{\partial {\mathsf{\Omega }}^N}·\mathit{n},\hfill & \hfill \frac{\partial T}{\partial \mathit{n}}& =\partial T^{\partial {\mathsf{\Omega }}^N},\hfill \end{array}$$

(9)

where ${\mathit{u}}^{\partial \mathsf{\Omega },D},T^{\partial \mathsf{\Omega },D},\nabla u^{\partial \mathsf{\Omega },N},\partial T^{\partial \mathsf{\Omega },N}$ are user-prescribed values and $\partial {\mathsf{\Omega }}^D/\partial {\mathsf{\Omega }}^N$ correspond to the location on the boundary where Dirichlet/Neumann BCs are prescribed. Equation (6) models the motion of the color function, representing the phase front motion at velocity $\mathit{v}$. This velocity field only exists at the phase front $\mathsf{\Gamma }$ and we extend it to zero elsewhere.

To close the system, the thermal and motion conditions at the phase front $\mathsf{\Gamma }$ have to be modeled and are the object of the following section.

2.3. Conditions at the Phase Front

2.3.1. Stefan Condition at the Moving Phase Front

The interface $\mathsf{\Gamma }$ modeling the moving phase front, is represented by the color function $\chi$, and, is located at ${\delta }^{\mathsf{\Gamma }}$ with normal ${\mathit{n}}^{\mathsf{\Gamma }}$ and is moving at velocity $\mathit{v}$. The model describing its displacement is based on the conservation of mass and energy across it: the so-called Stefan condition [4,13,14]. Indeed the heat fluxes $\mathit{\phi }$ on each side of $\mathsf{\Gamma }$ and its magnitude in the normal direction are given by Fourier’s law as, $k=l,s$

$$\begin{array}{ccc}\hfill {\mathit{\phi }}_k& =-{\lambda }_k\nabla T_k,{\phi }_k& ={\mathit{\phi }}_k·{\mathit{n}}^{\mathsf{\Gamma }},\hfill \end{array}$$

(10)

where $\lambda$, T are the conductivity and temperature in the vicinity of the front. The heat fluxes on both sides (liquid and solid) of the phase front are a priori different in terms of direction and magnitude, and their difference is directly linked to the phase change rate through the so-called Stefan condition:

$$\begin{array}{cc}\hfill \dot{m}& {\displaystyle =\frac{{\phi }_l-{\phi }_s}{{\mathcal{H}}^{*}},}\end{array}$$

(11)

where $\dot{m}$ is the mass flux per unit area crossing the interface due to the phase change and ${\mathcal{H}}^{*}$ is the strictly positive latent heat. ${\mathcal{H}}^{*}$ takes either the value of the latent heat of fusion/solidification, ${\mathcal{H}}^{\mathrm{fus}}$/${\mathcal{H}}^{\mathrm{sol}}$ according to the sign of ${\phi }_l-{\phi }_s$. If ${\phi }_l-{\phi }_s<0$ ($\dot{m}$ is negative), the liquid is solidifying and ${\mathcal{H}}^{*}={\mathcal{H}}^{\mathrm{sol}}$. Otherwise the solid is melting and ${\mathcal{H}}^{*}={\mathcal{H}}^{\mathrm{fus}}$. Thus:

$$\begin{array}{cc}\hfill {\mathcal{H}}^{*}& =\left\{\begin{array}{c}{\mathcal{H}}^{\mathrm{fus}}\mathrm{if}{\phi }_l-{\phi }_s>0,\hfill \\ {\mathcal{H}}^{\mathrm{sol}}\mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$

(12)

Moreover written in terms of density and velocity, we can deduce the phase front velocity as a function of the mass flux $\dot{m}$ and the local density $\rho$ as

$$\begin{array}{c}\hfill {\displaystyle \mathit{v}=-\frac{\dot{m}}{\rho }{\mathit{n}}^{\mathsf{\Gamma }}.}\end{array}$$

(13)

The local density is determined according to the composition of the solid/liquid phase at the vicinity of the interface. Assuming an immobile solid (${\mathit{u}}_s=0$) leads to consider that the thermal flux jump between the phases which generate some motion has to be taken into account in the liquid.

The velocity of the interface is then deduced as being

$$\begin{array}{cc}\hfill \mathit{v}& ={\mathit{u}}_k-\frac{\dot{m}}{{\rho }_k}{\mathit{n}}^{\mathsf{\Gamma }},\hfill \end{array}$$

(14)

where ${\mathit{u}}_k$ is the velocity of the phase k at the vicinity of the interface. Here, notice that the phase change imposes a velocity at the interface only if the densities of the phases are different, and, if the mass flux is non zero. Equation (14) is verified in both phases at the vicinity of the phase front, therefore the following condition is verified:

$$\begin{array}{c}\hfill {\displaystyle \mathit{v}=-\frac{\dot{m}}{{\rho }_s}{\mathit{n}}^{\mathsf{\Gamma }}={\mathit{u}}_l-\frac{\dot{m}}{{\rho }_l}{\mathit{n}}^{\mathsf{\Gamma }}⟹{\mathit{u}}_l=\dot{m}\left(\frac{1}{{\rho }_l}-\frac{1}{{\rho }_s}\right){\mathit{n}}^{\mathsf{\Gamma }}.}\end{array}$$

(15)

Notice that, if the liquid and solid normal heat fluxes are equal, i.e., ${\phi }_s={\phi }_l$, then the front remains static and we observe a pure heat conduction situation. Contrarily in a solidification/fusion situation then ${\phi }_l-{\phi }_s$ is strictly non-zero and the phase front moves.

2.3.2. Interface States

Let us assume that two extreme temperatures at the interface can be reached depending on the materials’ composition: $T^{\mathrm{fus}}$ the temperature of fusion of the solid phase, and, $T^{\mathrm{sol}}$ the temperature of solidification of the liquid phase such that $T^{\mathrm{fus}}\ge T^{\mathrm{sol}}$. Defining the phase front temperature $T\equiv T^{\mathsf{\Gamma }}$ for any $\mathit{x}\in \mathsf{\Gamma }\left(t\right)$ and all $t>0$, we consider three possible so-called ‘states’ leading to three conditions at the front:

$$\begin{array}{cccc}\hfill \mathrm{Solidification}\mathrm{state}& \iff \dot{m}\le 0\hfill & \hfill \Rightarrow & T^{\mathsf{\Gamma }}=T^{\mathrm{sol}},\hfill \end{array}$$

(16)

$$\begin{array}{cccc}\hfill \mathrm{Fusion}\mathrm{state}& \iff \dot{m}\ge 0\hfill & \hfill \Rightarrow & T^{\mathsf{\Gamma }}=T^{\mathrm{fus}},\hfill \end{array}$$

(17)

$$\begin{array}{cccc}\hfill \mathrm{Conduction}\mathrm{state}& \iff T^{\mathrm{sol}}<T^{\mathsf{\Gamma }}<T^{\mathrm{fus}}\hfill & \hfill \Rightarrow & \dot{m}=0.\hfill \end{array}$$

(18)

Notice that the phase front is moving for the first two situations while it is static for the last one, see Figure 1 for an illustration.

From those considerations we set the following thermal closures:

$$\begin{array}{cc}\hfill T=T^{\mathrm{sol}},& \mathrm{if}\dot{m}\le 0,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\mathit{\phi }}_s·{\mathit{n}}^{\mathsf{\Gamma }}={\mathit{\phi }}_l·{\mathit{n}}^{\mathsf{\Gamma }},& \mathrm{if}{T}^{\mathrm{sol}}<T^{\mathsf{\Gamma }}<T^{\mathrm{fus}},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill T=T^{\mathrm{fus}},& \mathrm{if}\dot{m}>0.\hfill \end{array}$$

(21)

Notice that the conduction equation ${\mathit{\phi }}_s·{\mathit{n}}^{\mathsf{\Gamma }}={\mathit{\phi }}_l·{\mathit{n}}^{\mathsf{\Gamma }}$ is a consequence of the relation $\dot{m}=0$, see (11). In the conduction state, the temperature must be deduced from the conduction Equation (20). These conditions are recalled in the logical diagram on Figure 2.

As a final remark of this modeling section let us point out that in the case of homogeneous liquid and solid density, ${\rho }_l={\rho }_s$ then the right hand side of Equation (3) is identically null, and it amounts to solve a classical incompressible condition. Otherwise, when ${\rho }_l\ne {\rho }_s$, the right hand side implies a volume change along with mass transfer between the phases. Finally system (3)–(6) along with BCs (8)–(9) is closed with phase front conditions (19)–(21). The numerical method designed to solve this problem is presented in the next section.

3. Numerical Methods

3.1. Rationale

The numerical method employed to solved the previous system of PDEs is inspired, partly re-used from Bois et al. [15,16] and implemented into the open-source code TrioCFD [12]. The code TrioCFD has been originally designed to deal with liquid/vapor applications. Contrarily, our context involves solid/liquid phase changes, and, therefore, the possible creation and presence of solids with the specificity that solid and liquid can be different pure materials. Consequently the original numerical methods in TrioCFD demand evolution and adaptation to fit with this solid/liquid phase change context. Namely, the one-fluid Navier–Stokes (NS) equations are solved by a prediction-correction scheme. For our context this demands to be penalized in the solids using for instance some Immersed-Boundary-Conditions (IBC) technique [9]. The phase change front evolution is captured by a front-tracking (FT) method [5] where the interfaces are tracked with a $d-1$ moving surface mesh, where d is the space dimension. This moving mesh technique is further coupled with a Volume-Of-Fluid (VoF) method employing a Level-Set (LS) distance function [17] in order to ensure a semi-local mass conservation. At last a Ghost Fluid Method [1] is used to solve the heat transfer within each phase separately. In this section we present the main aspects of the numerical methods developed in TrioCFD and their evolution to fit our context.

3.2. Time Discretization and Model Coupling

For the sake of clarity, an Euler explicit time discretization scheme is considered in this work. The time segment $[0,T]$ is split into successive time-steps $[t^n,t^{n+1}]$ with $\mathsf{\Delta }t\equiv \mathsf{\Delta }{t}^n=t^{n+1}-t^n>0$ and $n\in \mathbb{N}$. In what follows we identify any variable or parameter $\varphi$ at time $t^n$ as ${\varphi }^n$. Velocity and position of the materials are known at $t^n$. Then the models describing the fluids, interface and solids are solved successively in a mildly coupled manner. Namely the velocity is first updated by the monolithic penalized NS equations. Secondly, the phase front is displaced due to the new temperature profile generating some phase changes so that the color function $\chi$ must be updated. Thirdly, the temperature field is computed within each phase by solving the heat equation over the time-step. Note that the extension to high order of accuracy in time via Runge-Kutta schemes, recast as successive explicit Euler steps, is feasible within TrioCFD but out of the scope of this work.

3.3. Mesh and Front Representation

Computational Mesh

The computational domain at time $t^n$ is an evolving polygon/polyhedron $\mathsf{\Omega }\left(t^n\right)$ in 2D/3D. For the sake of clarity we present the 2D discretization reminding that the 3D extension follows the same concepts. $\mathsf{\Omega }\left(t^n\right)$ is covered without gap nor hole by an unstructured mesh made of fixed quadrangles, also called cells, and denoted by ${\omega }_c$ for $c=1,\dots ,N_c$. This mesh is referred to as ${\mathcal{M}}_{\mathsf{\Omega }}$.

The phase front of dimension $d=1$, represented by the color function $\chi$, is meshed by a list of ordered points, called markers, each referred to by ${\mathit{\xi }}_k$, $k=1,\dots ,N_k$. This 1D mesh is called ${\mathcal{M}}_{\chi }$ and evolves within ${\mathcal{M}}_{\mathsf{\Omega }}$. The number of markers $N_k$ may increase or decrease to ensure an acceptable geometrical quality of the phase front. At each marker we compute the normal as the sum of the two normal vectors of the front segments impinging at the marker.

As such any quadrangular cell can be either ‘pure’ or ‘mixed’ if ${\mathcal{M}}_{\chi }$ crosses ${\omega }_c$ or not. If ‘pure’, the cell is entirely included in the liquid phase or in the solid one. If ‘mixed’ then at least one or two pieces interface segments cross the cell splitting it into pure liquid and solid states, see Figure 3.

3.4. Navier–Stokes Discretization via a Penalized Prediction-Correction Method

The Navier–Stokes equations are solved using a so-called prediction-projection method on a staggered grid with a Marker-And-Cell (MAC) discretization [3]. The velocity components are located at the cell faces, while the scalar transport variables (pressure, temperature, level set function) and the thermodynamic properties (density, viscosity, thermal conductivity, heat capacity) are computed at the cell center. This method computes the velocity field on the whole domain, and so, a penalized method is used in combination to ensure a vanishing velocity in the solid. For this purpose, the penalized method introduced initially in Belliard [11] without phase change has been adapted to our context.

The time-dependent incompressible Navier–Stokes equations are solved explicitly in time using the standard projection method [2]: from the quantity known at $t^n$, the pressure $P^{n+1}$ and velocity ${\mathit{u}}^{n+1}$ are deduced using the prediction-correction method [16], leading to the discrete momentum equation. The spatial discrete operators considered in this study correspond to those initially implemented in TrioCFD [18] and are identified with index h. In particular the convection operator refers to the QUICK scheme (Quadratic Upstream Interpolation for Convective Kinematics [19]) while the diffusion operator corresponds to a usual second order finite difference scheme. Let us describe these two prediction and correction steps in the following.

Penalized prediction: This step consists in using an Euler explicit time scheme to compute a predicted velocity ${\mathit{u}}^{*}$ neglecting the pressure term in (4). The zero-velocity constraint on the solid phase is dealt with by adding a penalization term which is controlled by an arbitrarily small parameter $\eta ={10}^{-12}\ll 1$.

Introducing the notation:

$$\begin{array}{cc}\hfill \mathcal{T}& ={\nabla }_h·\left({\mathit{u}}^n{\mathit{u}}^n\right)+\frac{1}{{\rho }^n}\nabla ·\left({\mu }^n({\nabla }_h{\mathit{u}}^n+{\nabla }_h^T{\mathit{u}}^n)\right)+\mathit{g}-\frac{{\rho }_l}{{\rho }^n}\beta (T^n-T_0)(1-{\chi }^n)\mathit{g},\hfill \end{array}$$

(22)

the predicted velocity, ${\mathit{u}}^{*}$, is computed as:

$$\begin{array}{cc}\hfill {\mathit{u}}^{*}& {\displaystyle ={\mathit{u}}^n+\mathsf{\Delta }t\mathcal{T}-\frac{\chi }{\eta }{\mathit{u}}^{*}.}\end{array}$$

(23)

Reformulating (23) as

$$\begin{array}{c}\hfill {\displaystyle {\mathit{u}}^{*}=\frac{1}{{\displaystyle 1+\frac{\chi }{\eta }}}\left({\mathit{u}}^n+\mathsf{\Delta }t\mathcal{T}\right),}\end{array}$$

(24)

we notice that, numerically, the predicted velocity ${\mathit{u}}^{*}$ vanishes in the solid phase, where $\chi$ equals 1, as $(1+\frac{\chi }{\eta })\underset{\eta \to 0}{⟶}\infty$, while it remains unchanged in the liquid phase where $\chi =0$.

Penalized correction: Following [11], the penalized correction consists in computing the corrected velocity which fulfills the constraint (3) as

$$\begin{array}{cc}\hfill {\mathit{u}}^{n+1}& {\displaystyle ={\mathit{u}}^{*}+\frac{\mathsf{\Delta }t}{{\rho }_f^n}{\nabla }_h{P}^{n+1}+\frac{\chi }{\eta }({\mathit{u}}^{*}-{\mathit{u}}^{n+1}),}\hfill \end{array}$$

(25)

where $P^{n+1}$ has to be preliminary determined and ${\rho }_f^n$ corresponds to an averaged density at the current face usually equal to ${\rho }_l$ in the liquid phase. The correction step presented in [11] does not take into account any phase which violates the non-zero divergence Equation (3). On the other hand, the numerical discretization of this term introduced initially in TrioCFD [15] uses a signed distance (defined in Section 3.5), and, has been derived without considering any penalized method in the context of water vapor/liquid applications, considering therefore a ratio of density of about 1000 between both phases. Our context, and others in engineering, witness density ratio of order unity and the presence of solid or quasi-solids.

Therefore, we propose here to integrate in any mixed cell the term, $\dot{m}\left(\frac{1}{{\rho }_l}-\frac{1}{{\rho }_s}\right)$ appearing in (3), along the interface to get a discrete volume term as:

$$\begin{array}{cc}\hfill {\zeta }^n& \equiv \zeta \left(t^n\right)=\frac{∥{\mathsf{\Gamma }}^n\cap {\omega }_c∥}{∥{\omega }_c∥}{\dot{m}}^n\left(\frac{1}{{\rho }_l}-\frac{1}{{\rho }_s}\right).\hfill \end{array}$$

(26)

Here ${\displaystyle \frac{∥{\mathsf{\Gamma }}^n\cap {\omega }_c∥}{∥{\omega }_c∥}}$ is the length (or surface in 3D) of the front in the given mixed cell, and, ${\dot{m}}^n$ corresponds to the mass flux across this front (detailed later in Section 3.5).

Next, the classical correction step [15] is reformulated regarding the contribution of the penalized term introduced in (23). The pressure term $P^{n+1}$ is then deduced from

$$\begin{array}{cc}\hfill \frac{\eta }{\eta +\chi }{\nabla }_h·\left(\frac{1}{{\rho }_f^n}{\nabla }_h{P}^{n+1}\right)& =-\frac{1}{\mathsf{\Delta }t}{\nabla }_h·{\mathit{u}}^{*}+\frac{1}{\mathsf{\Delta }t}{\zeta }^n.\hfill \end{array}$$

(27)

However, the term ${\zeta }^n$ is a priori defined in a mixed cell. Hence the parameter $\eta$, intended to go to zero in cells where the solid is present, generates that the numerical contribution of the phase change ${\zeta }^n$ in (27) is numerically canceled. To overcome this undesirable behavior, we made the choice to affect this discrete volume change term to the closest pure liquid cells. Recall that these liquid cells are not subject to any penalization term. This choice is somewhat justified by the fact that only the liquid has to endure the whole volume change constraint, while the solid has to maintain a zero-velocity. This numerical algorithm allowing to transfer the mixed cell contribution of ${\zeta }^n$ to the closest pure liquid cells is not detailed here but relies on crude interpolation techniques. As such the contribution of ${\zeta }^n$ for a given mixed cell is consistent with the corresponding global volume change computed in the corresponding pure liquid cells.

3.5. Interface Treatment: Motion and Geometrical Regularization

As already mentioned the phase front is represented by a sharp interface, that is a surface mesh made of straight segments in 2D (and co-planar triangles in 3D).

3.5.1. Displacement

The displacement of the markers drives the motion of the interface. A level-set method is adopted for capturing the phase front interface dynamics. $\mathsf{\Gamma }$ is represented using the zero level of a continuous surface Level Set (LS) function $\varphi$.

$$\begin{array}{c}\hfill \mathsf{\Gamma }\left(t\right)=\left\{\mathit{x}\in \mathsf{\Omega },\mathrm{s}.\mathrm{t}.\varphi (\mathit{x},t)=0\right\},\end{array}$$

(28)

where the level set function is represented by a signed distance function $\varphi (\mathit{x},t)=\pm \left|d\right|$, where d represents the shortest distance from any grid point to the interface with the signs $\varphi >0$ in the liquid phase and $\varphi <0$ in the solid one. The advection of the level set function is performed by the non-conservative equation

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \varphi }{\partial t}+\mathit{v}·\nabla \varphi =0,}\end{array}$$

(29)

where $\mathit{v}$ is physically valid only at the interface $\mathsf{\Gamma }$ but is needed in the small region close to it. To ensure that the signed distance function property of $\varphi$ is preserved during the computations, $\varphi$ can be periodically re-initialized, loosing as such the conservation of the phases. To counter-fight this effect, a Volume-of-Fluid (VoF) method is supplemented to compute the conservative evolution of the colour function $\chi$. The discrepancy between the conservative VoF volumes and the LS ones is used in the regularization/smoothing algorithms.

3.5.2. Regularization/Smoothing

In order to preserve sufficient geometrical quality of the interface mesh during the simulation, regularization and smoothing algorithms are employed within TrioCFD. These geometrical modifications are performed under the constraint that the mass of each phase is conserved as best as possible. First, the remeshing procedure maintains the mesh edge length within a prescribed range. It relies on deleting markers on small edges, or adding some at the middle of too long ones. As such the connectivity of the interface mesh does change. Second, the smoothing procedure consists in a regular redistribution the markers without connectivity changes along/on the interface mesh in 2D/3D. Somewhat it corresponds to a Laplacian smoothing on the marker position under the constraint that the associated changes of volume/mass remain controlled. The detailed description of these algorithms can be found in [15,16].

3.6. Ghost-Fluid Approach for Heat Transfer Equation

The ghost-fluid approach is employed to solve the energy equation in each phase governed by (5) in terms of temperature variations. The moving phase front renders this resolution difficult because the velocity and heat flux are discontinuous across it. To leverage this difficulty the Ghost-Fluid approach is employed as it replaces the discontinuity by a continuous extension of Eulerian fields across the interface. The temperature being defined at the cell-centers, so will be the temperature gradient normal to the interface $\nabla T·\mathit{n}$ and the phase-change rate $\dot{m}$.

3.6.1. Temperature Gradient and Rate of Change

In the vicinity of the interface, first, based on linear extrapolation of the liquid temperature from pure liquid cells, the ‘liquid ghost temperatures’ are extrapolated in the ghost cells. From the interface temperature, $T^{\mathsf{\Gamma }}$, assumed to be given (see later for its determination) and the distance function, d, the normal component of the temperature gradient at the interface, $\nabla T·\mathit{n}$, is readily obtained for both phases as

$$\begin{array}{cccc}\hfill (\stackrel{g}{\nabla }{T}_s·\mathit{n})& =\frac{T^{\mathsf{\Gamma }}-T_s}{d},\hfill & \hfill (\stackrel{g}{\nabla }{T}_l·\mathit{n})& =\frac{T^{\mathsf{\Gamma }}-T_l}{d}.\hfill \end{array}$$

(30)

As such these cell-centered fields represent the normal heat fluxes at the interface and can be computed for any liquid or solid cell. Once mixed and pure solid cells get a gradient, it is possible to assign a ‘liquid ghost temperature’ in the solid phase and, vice-versa, a ‘solid ghost temperature’ in the liquid phase which are further used for the next time-step.

Using the ghost gradients in each phase, we are able to calculate $\dot{m}$ corresponding to a discrete mass rate of change due to phase changes in mixed cells:

$$\begin{array}{cc}\hfill \dot{m}& =\frac{{\lambda }_l(\stackrel{g}{\nabla }{T}_l·\mathit{n})-{\lambda }_s(\stackrel{g}{\nabla }{T}_s·\mathit{n})}{{\mathcal{H}}^{*}},\hfill \end{array}$$

(31)

where ${\mathcal{H}}^{*}$ is the enthalpy of fusion or solidification which depends on the interface state (see next paragraph). $\dot{m}$ in (31) is a priori only defined in mixed cells, but a prolongation is carried on to extend this phase change rate $\dot{m}$ in the vicinity of the interface, that is in pure solid cells.

3.6.2. Interface Temperature/Interface State

In our model, and as shown in Figure 1, the interface temperature can only vary between $T^{\mathrm{fus}}$ and $T^{\mathrm{sol}}$. As such the ‘interface state’ becomes a variable of interest and we follow it with an Eulerian cell-centered field which takes values $S/C/F$ for solidification/conduction/fusion cases. At each time step, the state of the interface is updated in mixed cells and in their vicinity, following the logigram Figure 2. The interface temperature in mixed cell is deduced using the following relation

$$\begin{array}{cc}\hfill T^{\mathsf{\Gamma }}& =\left\{\begin{array}{cc}{T}^{\mathrm{fus}}\hfill & \mathrm{if}\mathrm{the}\mathrm{mixed}\mathrm{cell}\mathrm{is}\mathrm{in}\mathrm{fusion},\hfill \\ T^{\mathrm{cond}}\hfill & \mathrm{if}\mathrm{the}\mathrm{mixed}\mathrm{cell}\mathrm{is}\mathrm{in}\mathrm{conduction},\hfill \\ T^{\mathrm{sol}}\hfill & \mathrm{if}\mathrm{the}\mathrm{mixed}\mathrm{cell}\mathrm{is}\mathrm{in}\mathrm{solidification}.\hfill \end{array}\right.\hfill \end{array}$$

(32)

where $T^{\mathrm{cond}}$ is the last unknown when conduction occurs. To determine $T^{\mathrm{cond}}$ in this work we rely on the heat flux equality through the interface. A simple linear interpolation is performed between two local averaged temperatures/distances of the solid and liquid phase, $\overline{T_s}$, $\overline{T_l}$, $\overline{d_s}$, and $\overline{d_l}$. While in 1D the situation is simple and there is no need to compute any of those averaged values, see Figure 4, in 2D and 3D there is no unique temperature and distance to the interface to rely on. In practice one computes the average of pure cell values considering the $3\times 3$ stencil around the current cell. At last the interface temperature is obtained thanks to the formula

$$\begin{array}{cc}\hfill T^{\mathrm{cond}}& {\displaystyle =\frac{{\lambda }_s\overline{T_s}\overline{d_l}-{\lambda }_l\overline{T_l}\overline{d_s}}{{\lambda }_s\overline{d_l}-{\lambda }_l\overline{d_s}}.}\hfill \end{array}$$

(33)

This equation ends the description of the numerical method which has been partly reemployed and adapted within the simulation code TrioCFD. In the next section a validation and verification campaign is carried on.

4. Validation and Verification Test Suite

4.1. Stefan Problems

In order to validate our numerical approach we simulate the so-called 1D Stefan problem. This test case has an analytical solution [4,13] towards which we can observe the convergence of the space-time evolving numerical solution as the mesh/time steps are refined. The problem consists in simulating the 1D solidification (resp. fusion) of an homogeneous material, for which the solid and liquid phases have the same material properties. The semi infinite computational domain corresponds to $\mathsf{\Omega }={\mathsf{\Omega }}_s\left(t\right)\bigcup {\mathsf{\Omega }}_l\left(t\right)$ with the solid represented by ${\mathsf{\Omega }}_s=[0,s\left(t\right)]$ where $s\left(t\right)$ is the location of the moving front in-between for any time $t>0$. A semi infinite liquid domain should stand on the right of ${\mathsf{\Omega }}_s$ at constant solidification temperature $T^{\mathrm{sol}}$ but we restrict it to ${\mathsf{\Omega }}_l=[s\left(t\right),10]$ so that the right boundary is too far to interfere with the computation. An external temperature $T^{\mathrm{ext}}<T^{\mathrm{sol}}$ is enforced at the left boundary of ${\mathsf{\Omega }}_s$. The aim is to compute the position of the interface $s\left(t\right)$ and the temperature profile $T(x,t)$. The governing equations of the 1D Stefan problem are:

$$\begin{array}{c}\hfill \frac{\partial T}{\partial t}(x,t)=\alpha \frac{{\partial }^{2}T}{\partial x^2}(x,t),\beta \frac{ds}{dt}\left(s\left(t\right),t\right)=-\nabla T\left(s\left(t\right),t\right),\end{array}$$

(34)

where ${\displaystyle \alpha =\frac{\lambda }{\rho C_p}}$, ${\displaystyle \beta =\frac{\rho \mathsf{\Delta }{\mathcal{H}}^{\mathrm{sol}}}{\lambda }}$, and for all $x\in \mathsf{\Omega }$ and $t\in [0,t_{\mathrm{final}}]$. The boundary and initial condition are given by: $T(0,t)=T^{\mathrm{ext}}$, $T(s\left(t\right),t)=T^{\mathrm{sol}}$, and $s(t=0)=0$. The unique solution of this problem is given by [4]:

$$\begin{array}{ccc}\hfill T(x,t)& =& {\displaystyle T^{\mathrm{ext}}+(T^{\mathrm{sol}}-T^{\mathrm{ext}})\frac{\mathrm{erf}\left({\displaystyle \frac{x}{2\sqrt{\alpha t}}}\right)}{\mathrm{erf}\left(k\right)},}\hfill \end{array}$$

(35)

$$\begin{array}{ccc}\hfill s\left(t\right)& =& {\displaystyle 2k\sqrt{\alpha t},k{e}^{k^2}\mathrm{erf}\left(k\right)=\frac{1}{\beta \alpha \pi }\left(T^{\mathrm{sol}}-T^{\mathrm{ext}}\right),}\hfill \end{array}$$

(36)

where erf is the error function. k is approximated by Newton’s method on the last equation and further substituted into (35) and (36)-left. We consider a rectangular domain ${\mathsf{\Omega }}_s$ of 4 units in length by $0.5$ in height. This domain is paved by ${\mathcal{N}}_x$ cells in length and ${\mathcal{N}}_y=9$ cells in height, ${\mathcal{N}}_x$ is subject to take increasing values.

Since our modeling does not handle the creation of solid initially, our simulations start with a pre-existing solid. Hence we assume an initial front position at $s\left(t_{\mathrm{init}}\right)=0.5$ and the parameter values $\alpha =1.0$, $\beta =1.0$ and $T^{\mathrm{ext}}=-1.0$. The initial temperature of the solid is given by the exact solution: ${\displaystyle T_s(x,t_{\mathrm{init}})=-1+\frac{1}{\mathrm{erf}\left(k\right)}\mathrm{erf}\left(\frac{x}{2\sqrt{t_{\mathrm{init}}}}\right)}$ where $t_{\mathrm{init}}=0.2893$. The initial temperature and velocity of the fluid are set to 0. Then we set the following constant values for the physical parameters: $\lambda =1$, $C_p=1$, $\mu =1$ and ${\mathcal{H}}^{\mathrm{fus}}={\mathcal{H}}^{\mathrm{sol}}=1$. Moreover we consider two different densities for the liquid and solid: ${\rho }_l=2$ and ${\rho }_s=1$. Then two configurations are simulated: one where $T^{\mathrm{fus}}=T^{\mathrm{sol}}$, then one for which $T^{\mathrm{fus}}\ne T^{\mathrm{sol}}$.

4.1.1. Configuration $T^{\mathrm{fus}}=T^{\mathrm{sol}}=0$

The results in

Table 1. Stefan problem with configuration ${\rho }_l\ne {\rho }_s$ with fixed time step $\mathsf{\Delta }t={10}^{-4}$—Error and order of convergence in $L^1$, $L^2$ and $L^{\infty }$ norms—$\mathsf{\Delta }x$ space step varying between $0.08$ and $0.008$.

${\mathit{N}}_{\mathit{x}}$	$\mathsf{\Delta }\mathit{x}$	${\mathit{L}}^1$ Error	${\mathit{L}}^1$ Order	${\mathit{L}}^2$ Error	${\mathit{L}}^2$ Order	${\mathit{L}}^{\mathit{\infty }}$ Error	${\mathit{L}}^{\mathit{\infty }}$ Order
50	0.08	$3.728E-3$	–	$4.072E-3$	–	$5.844E-3$	–
75	0.053	$2.531E-3$	$0.94$	$2.752E-3$	$0.95$	$3.950E-3$	$0.95$
100	0.04	$1.903E-3$	$1.01$	$2.068E-3$	$1.02$	$2.992E-3$	$0.99$
150	0.027	$1.255E-3$	$1.05$	$1.368E-3$	$1.05$	$2.015E-3$	$1.01$
200	0.02	$9.235E-4$	$1.02$	$1.011E-3$	$1.01$	$1.477E-3$	$1.03$
400	0.01	$4.713E-4$	$0.97$	$5.137E-4$	$0.98$	$7.48E-4$	$0.98$
500	0.008	$3.842E-4$	$0.92$	$4.172E-4$	$0.93$	$6.05E-4$	$0.95$
–	–	$L^1$ average →	0.97	$L^2$ average →	0.99	$L^{\infty }$ average →	0.99

gather the errors in norm $L^1$, $L^2$ and $L^{\infty }$ obtained by our numerical modeling for decreasing values of $\mathsf{\Delta }x$. Each simulation is performed until the final time $t^{\max }=0.5$ with a fixed time step equal to $\mathsf{\Delta }t={10}^{-4}$. The results in Table 1 show a reduction of the errors, proving numerically that the numerical method converges at order 1 as expected.

Next, Figure 5-left shows the position of the interface over time for each mesh compared to the exact position. Each simulation correctly approximates the exact solution, and, the finer the mesh the better the approximation (see the zoom). Figure 5-right shows the position of the interface over time for different time steps with a fixed space step $\mathsf{\Delta }x=0.01$ corresponding to ${\mathcal{N}}_x=400$. We observe that our numerical results fit relatively well the exact solution. If we further increase the time step, the method diverges because the coupling of the equations requires a small time step for the numerical method to remain stable. This restricted time step leads to a situation for which the error produced by the discretization in time is much smaller than the spacial error.

4.1.2. Configuration $T^{\mathrm{fus}}\ne T^{\mathrm{sol}}$

Here, the solidification temperature is set to $T^{\mathrm{sol}}=-0.5$ and the melting temperature to $T^{\mathrm{fus}}=0.5$. In this situation, while the liquid is heated up on the right boundary, it cools down at the interface $s\left(t\right)$. Hence, one expects a scenario for which a cooling period at the beginning of the test case, corresponding to a pure conduction period, is followed by a solidification period. No analytical solution is available. It is possible to assess the starting time of the solidification which corresponds to the time when the temperature at the interface reaches $T^{\mathrm{sol}}$ by $t=0.2748$. In Figure 6 we plot the interface position $s\left(t\right)$ as a function of time. We clearly observe no motion of the interface (conduction state) at the beginning, then the displacement of the interface and its progressive stop (solidification state), and no motion anymore (back to conduction state). We can observe that the time of transition to solidification is well approximated by $t\approx 0.27$. The solidification stops around $t\approx 0.8$, but there is no analytical solution to compare this result with. Stefan problem was a sanity check used to prove that the model solved by our adapted numerical method in TrioCFD is able to capture the phenomenology generated by a moving interface driven by solidification/fusion and conduction periods. More sanity checks have been conducted but we omit them in this paper.

4.2. Two-Dimensional Corner Stefan-like Test Case

This double Stefan-like test case consists in the solidification of a semi-infinite 2D domain with a corner, see a analysis and numerical studies in [20,21,22]. The fixed computational domain is $\mathsf{\Omega }={[0;1]}^2$ which initially is split into a liquid phase at rest ${\mathsf{\Omega }}_l=[E_x(t=0);1]\times [E_y(t=0);1]$ at solidification temperature $T_{\mathrm{sol}}$ and a solid one ${\mathsf{\Omega }}_s=\mathsf{\Omega }{\mathsf{\Omega }}_l$. The initial corner shaped solid thickness is set to $E_x(t=0)=E_y(t=0)=0.1$. The mesh is set to $110\times 110$ cells and the time step is fixed at $\mathsf{\Delta }t={10}^{-4}$. Initial temperature is constant at $T=0$. The temperature at boundaries $x=0$ (bottom) and $y=0$ (left) is fixed at a temperature $T^{\mathrm{ext}}=-1$ lower than the solidification temperature, $T^{\mathrm{sol}}=-0.5$, so that the solidification phenomenon starts because of these cold boundaries. Notice that in this test $T^{\mathrm{sol}}\ne T^{\mathrm{fus}}$. The liquid domain is closed but for the top-right corner where an opened BC is set for the liquid to escape if appropriate, the rest of the boundaries are of slip-wall type. This initial configuration is represented in Figure 7-left along with its associated physical parameters. Far away from the corner, the computed solution should be close to the classical 1D Stefan solution in x and y directions. Contrarily, in the vicinity of the corner, genuine 2D effects are expected.

Figure 8 shows the results obtained for four successive times $t=0.005$, $0.01$, $0.015$, and $0.25$. At time $t=0.005$ (top-left), the interface has not yet reached the solidification temperature, and, consequently, the interface state is in pure conduction, there is no phase change, nor any movement of the front. On the other hand, we observe a cooling phenomenon occurring from the boundaries. More importantly a non-constant temperature profile is observed on the liquid side proving that heat conduction takes place across the interface. At time $t=0.01$ (top-right) in the vicinity of the interface, close to the corner, solidification starts, while, far from the corner, the pure conduction situation still occurs. We observe that this part of the interface close to the corner is set in motion, moves, deforms, and retroactively generates a velocity on the liquid because of the change of volume linked to the difference of densities. Also one observes that the liquid tends to escape from the top-right open boundary as expected. On the figure at time $t=0.015$ (bottom-left) the solidification state propagates along the phase change front symmetrically upward and rightward. Finally on the panel at time $t=0.025$ (bottom-right) the solidification state has reached the whole length of the front, since the latter is at a temperature $T^{\mathrm{sol}}=-0.5$. As such the whole solidification front is set in motion.

Thus, this test case allows to simulate different features of the implemented numerical method: pure conduction through the front, solidification of the liquid setting in motion of the front and changing the temperature and velocity fields, and the regular transition between the states. Moreover, the coupling with the penalization method allows to maintain and ensure a zero-velocity in the solid. At last we have verified that far away from the corner the exact 1D Stefan solution is retrieved until the effects of the 2D corner interfere.

At last, a mesh independency study is presented in Figure 9 by comparing results for a similar test case (only modified by using $T^{\mathrm{sol}}=T^{\mathrm{fus}}=0$ and thus leading to less CPU-time consuming) on three successively Cartesian refined meshes: a coarse ($50\times 50$ cells), intermediate ($100\times 100$ cells) and a finer mesh ($150\times 150$ cells). As expected the refined meshes capture more accurately the flow, the fileds and the interface location. The numerical results seem to converge appropriately as can be seen on the fourth panel where the location of the interface is plotted in purple/black/white lines for the coarse/intermediate/fine mesh respectively. Some numerical comparisons against other methods can be also found in [21,22] for instance.

4.3. Focusing Effect within a Thin Metallic Layer

In this last test case, we propose to use our numerical modeling approach in a more industrial-like application. In the context of nuclear severe accident, one important research interest is the simulation of a thin metallic liquid layer on the top of a corium pool. This configuration occurs when an In-Vessel Retention (IVR) strategy is considered which consists in sustaining the integrity of the reactor pressure vessel by External Reactor Vessel Cooling [23,24,25,26,27]. To estimate the success of IVR strategy, an accurate evaluation of the heat fluxes that the corium pool imposes to the vessel lower head is needed. Due to the ablation of the vessel wall and to thermo-chemical effects, a thin metallic layer, in direct contact with the vessel wall, may lay on top of the corium pool. Such a thin layer may drastically increase the heat flux on the vessel up to possible failure of its integrity [28]: this is the so-called heat flux ‘focusing effect’. In the present work, simulation of the interaction between the thin metal layer and the vessel wall are considered in order to evaluate the effects of the phase change phenomena compared to the work done in [29]. Indeed, in [29] the vessel was not simulated and only single-phase thermalhydraulics simulations were performed (with a temperature set equal to the steel fusion temperature at the boundary). Thus, the objective of this current test case is to genuinely simulate the vessel in order to observe the impact of the melting phenomenon. Two-dimensional and three-dimensional simulations are proposed.

4.3.1. Test Case Description

A sketch of the configuration is presented in Figure 10. The global IVR context is illustrated on the left side. The actual geometry considered in the present simulations is depicted on the right of the figure. This geometry focuses on the coupling between the thin metallic layer and the vessel. The fixed global length (x-axis) is set at $L=200$ cm, initially 15 cm for the solid vessel and 185 cm for the liquid metal layer. A fixed height (z-axis) of $H=5$ cm is assumed while the width is in the range [$0,13$ cm]. Regarding thermal boundary conditions, the bottom surface of the liquid ($z=0$) is subjected to a constant heat flux of ${10}^6$ W/m${}^2$. The front, back (perpendicular to the y axis) and one lateral (at $x=2$ m) walls are adiabatic as well as for the solid bottom boundary. Contrarily the left lateral wall is set at a fixed temperature 400 K representative of the external vessel wall cooling by water nucleate boiling. At the top surface, fixed temperature is imposed at 1960 K for the liquid and adiabatic condition is assumed for the solid. For mechanical boundary conditions, no slip conditions are assumed. The properties of steel are assumed to be constant as in [29] and gathered in

Table 2. Physical properties of steel—Values taken from [29].

Parameter	Value
Density $\rho$	6720 ${\mathrm{kg}/\mathrm{m}}^3$
Coefficient of thermal expansion, $\beta$	$3\times {10}^{-5}$
Thermal conductivity, $\lambda$	20 W/mK
Fusion temperature, $T^{\mathrm{fus}}$	1658 K
Enthalpy of fusion, $\mathsf{\Delta }l{H}^{\mathrm{fus}}$	$2.76\times {10}^5$ J/kg
Heat capacity, $C_p$	674 J/kg K
Kinematic viscosity, $\mu$	$4.5696\times {10}^{-3}$ ${\mathrm{m}}^2/\mathrm{s}$

. Note that this rectangular box setup differs from the actual axisymmetric vessel geometry but was considered in [29] because it corresponds to the geometry of the BALI-Metal experiment [30] used for validation purpose.

4.3.2. Two-Dimensional Simulations

First of all, results of 2D simulations ($z-x$ plane) are presented in Figure 11 and Figure 12 in terms of the temperature and velocity fields computed in the thin metal layer and the vessel at different intermediate times and using a computational mesh constituted of $1200\times 100$ cells. We observe the displacement of the interface in the transient regime and the establishment of a multi-dimensional flow related to thermal gradients in the liquid. It is interesting to note that at the beginning of the transient, because of the way the temperature field is initialized, the interface starts moving to the right: the liquid slightly solidifies (see Figure 11). It is only for $t>50$ s that the vessel fusion occurs (see Figure 12). Near the interface, the vertical thermal gradients are more pronounced, reflecting the cooling effects near the interface. They are associated with a re-circulation that can be observed at different times in the figure. We also observe the setting up of Rayleigh-Bénard type convection cells far from the interface (and not influenced by the physical phenomena at the interface). This structure of the flow corresponds to the one expected and observed in the BALI experiment [30] and the simulations without taking into account the vessel presented in [26,29].

4.3.3. Three-Dimensional Simulations

Then, we present in Figure 13 results obtained for the 3D case with the same 2D input parameters given in Table 2 and Figure 10. At the top panel, the temperature field is presented while the velocity vectors and fusion front are given in the middle panel. The position of the interface for different slices orthogonal to the y axis is also depicted and compared to the interface computed in the previous 2D test case. According to the adiabatic boundary conditions applied on the front and back boundaries (perpendicular to the y axis), the main differences compared to the 2D case appear during the transient regime. Therefore the numerical results are presented at the very beginning of the simulation for two different times: 10 s and 30 s. As expected, the difference between 2D and 3D calculations remains limited in terms of the interface position because of these adiabatic boundary conditions on the front and back boundaries. Furthermore, these numerical results are intended to prove the ability of the code to handle the 3D displacements of the interface with phase change in the context of a complex flow. This proof of concept now demands more advanced and quantitative analyses. In terms of computation cost, this 3D simulation uses a computational mesh constituted of $1200\times 50\times 50$ cells and has required ten days of computing time using 64 CPU cores. The parallelization of our numerical development relies on the MPI parallel version of TrioCFD.

5. Conclusions and Perspectives

In this paper, we focus on the modeling and numerical simulation of two-phase flow involving solid-liquid phase changes. In our context, the fusion and solidification temperatures may be different which implies that a transition zone may exist in between the phase changes. This situation may be encountered when considering corium as a pure (liquid or solid) material for the study of severe accidents in nuclear power plants for instance. The modeling, numerical methods, and simulation code, must be adapted. In this work, we have developed, adapted and extended from [11,16] a model and code based

• on a monolithic Navier–Stokes system of equations penalized in the solid phases and adapted to the specificity of such phase changes,
• on a Front-Tracking approach to follow the phase change front supplemented with a Volume-of-Fluid method to counter-fight possible lack of mass conservation,
• on monolithic heat equation supplemented with the Ghost-Fluid approach to deal with the discontinuity at the front.

The numerical methods used to solve such a model have been presented in this paper. This numerical model has been further implemented/adapted in the code TrioCFD [12,18].

A numerical validation and verification campaign has been carried out, and, some results have been gathered such as the Stefan 1D and 2D validation tests and a more industrial-type 2D/3D test case corresponding to the analysis of the so-called ‘focusing effect’ pertaining to corium behavior in the realm of nuclear reactor safety.

The simulation tool, as such produced for general applications coupling liquid/solid with phase change, seems to be adequate to be used at a more operational level and this is one of our goal in a near future. From the numerical method point of view we plan to analyze the time-step determination, possibly considering implicit schemes. So far we have only used a semi-explicit time scheme with an ad-hoc time-step restriction. Moreover, a more in-depth study should be carried on to evaluate the accuracy of the thermal fluxes through the interface related to the use of the Ghost-Fluid method. This work takes a step forward in our general ongoing R&D effort to provide models and numerical methods for simulating the thermal-hydraulic behavior of a corium pool composed of different chemically reactive liquid phases with, possibly, liquid/solid phase change at their boundaries. The presented CFD results and, more importantly, the future ones will be used for engineering applications following an up-scaling methodology (see [31] for instance) in order to enhance modeling in the industrial integral code PROCOR [25].