Adrien Drouillet
email: adrien.drouillet@math.u-bordeaux.fr

Romain Le Tellier

Raphaël Loubère

Mathieu Peybernes

Louis Viot

Multi-dimensional simulation of phase change by a 0D-2D model coupling via Stefan condition

Keywords: Simulation of phase change, Fusion, Solidification, 0D multi-layer model, 2D heat conduction model, Model coupling

Considering phase changes associated with a high-temperature molten material cooled down from the outside, this work presents an improvement of the modelling and the numerical simulation of such processes for an application pertaining to the safety of light water nuclear reactors. Postulating a core meltdown accident, the behaviour of the core melt (aka corium) into a steel vessel is of tremendous importance when evaluating the vessel integrity. Evaluating correctly the heat fluxes requires the numerical simulation of the interaction between the liquid material and its solid counterpart which forms during the solidification process, but also may melt back. To simulate this configuration, encountered in various industrial applications, one considers a bi-phase model constituted by a liquid phase in contact and interaction with its solid phase. The liquid phase may solidify in presence of low energetic source, while the solid phase may melt due to an intense heat flux from the high-temperature liquid. In the frame of the in-house legacy code, several simplifying assumptions (0D multi-layer discretization, instantaneous heat transfer via a quadratic temperature profile in solids) are made for the modelling of such phase changes. In the present work, these shortcomings are illustrated and further overcome by solving a 2D heat conduction model in the solid by a mixed Raviart-Thomas finite element method coupled to the liquid phase due to heat and mass exchanges through Stefan condition. The liquid phase is modeled with a 0D multi-layer approach. The 0D-liquid and 2D-solid models are coupled by a Stefan like phase change interface model. Several sanity checks are performed to assess the validity of the approach on 1D and 2D academical configurations for which exact or reference solutions are available. Then more advanced situations (genuine multidimensional phase changes and an "industrial-like scenario") are simulated to verify the appropriate behavior of the obtained coupled simulation scheme.

Context and introduction

Numerical resolution of solidification and melting in coupled multi-scale and multi-physics problems is a substantial issue arising in several engineering fields. Because of the lack of full physical, phenomenological knowledge and because of computational resources limitation, heterogeneous models coupling is usually required. Some examples are multidimensional models (coupling of mixed 0D, 1D and 2D models for instance) coupled through an interface which purpose is to filter the information passed from one model or scale to another one. As such the modelling and numerical resolution of the phase change at the interface between materials is a major issue, even if recent advanced unified models treating all continuum dealing with phase changes are emerging [GLCTPRD21]. In particular, in this paper, we are interested in the simulation of solidification and melting phenomena involving a hot liquid (resp. cold solid container) which solidifies (resp. melts) on contact with a cold solid container (resp. hot liquid). This configuration is studied in various industrial applications such as continuous casting [AJM15], latent heat thermal energy storage [ZL13] and nuclear Severe Accidents (SAs) [Jac13,Seh12]. The latter is our target application even if our approach may adapt to any other one. During such an accident, the core of the nuclear reactor heats up leading to the formation of so-called 'corium' (i.e a mixed of molten core oxide and metallic liquid materials). This liquid releases heat volumetrically because of the fission products nuclear decay. Then, this decay heat is transferred by convection and conduction across the reactor pressure vessel and can lead to the ablation, melting and subsequent failure of the vessel, see figure 1 for a sketch. In nuclear power safety, various strategies have been proposed in order to mitigate a reactor SA and prevent the release of corium into the environment. One of them, named In-Vessel-Retention (IVR), consists in sustaining the integrity of the reactor pressure vessel by external reactor vessel cooling Associated simulations lead to coupled problems exhibiting tremendously different scales in time and space. Among all the complex and coupled transient phenomena, the behavior of the interfaces between the liquid and solid phases and the container has been identified to be a key phenomenon. Requiring fast-running models in particular when Monte-Carlo based sensitivity analyses are performed, most SA scenario codes that are commonly used for IVR strategy evaluation [LTSB15, SCS + 95, FI94, DSC + 09] use 0D models, also called "lumped parameter" models. In particular the phase change and the associated liquid/solid moving interface are mainly solved by 0D stationary models after many simplified and questionable modelling assumptions are adopted [CFS18, CFB + 19]. In this paper we propose to point out the main drawbacks of the modelling available in classical lumped parameter models where so-called "correlation-based closures" are needed.

The aim of this work is to propose an improvement based on a multidimensional coupling of the models by considering a transient 2D heat conduction equation model in the solid domain which solves the energy balance and the phase change. The 2D heat conduction model for the solid is coupled to the liquid via the classical Stefan condition allowing mass and heat exchanges via solidification or fusion. Several numerical methods have been proposed for solving Stefan-type problems classically considered for phase change modelling. Originally, fixed-grid methods are the most popular ones due to their computational efficiency and a relative simple implementation in multi-dimensions. They are mostly based on a diffusive approximation of the interface using an enthalpy method from a monolithic energy conservation formulation of the latent heat on the whole domain, see [Vol97] for more details. However these methods do not rely on an explicit representation of the interface which could be not accurate enough when different interfaces are considered and small scale geometrical entities appear. This is for instance the case in SA scenario when a "thin" metal phase on top of the liquid corium is formed, and leads to an intense heat flux concentration at that point (the so-called focusing effect [TLA + 97, RKCA98]). Such a situation among others requires an accurate geometrical description and capture of the phase front. In this context, a pure moving-grid method (with an explicit Lagrangian description of the phase front) may then seem more appropriate. But it would be inoperative or particularly complex for complex interface geometries. Level set methods dedicated to the tracking of sharp interface [CMOS97] do not seem to be suitable for multiple interfaces as they require several regularization steps which may lead to a complex algorithm. Moreover in our case the interfaces between the liquid and solids may appear ex-nihilo or disappear, see figure 1 (ablation and failure), which is a difficult situation to be handled by Level-Set techniques.

The goal of this work is then to track explicitly a genuine 2D phase front expecting to capture more accurately the heat flux between the phases. In this work, the liquid phase is assumed to be solved classically with a 0D model through correlation-based closures to evaluate the heat fluxes exerted on the solid. Contrarily the 2D heat equation will be solved in the solid by a mixed Raviart-Thomas Finite Element (FE) discretization on unstructured grid [RT77]. This method allows to access explicitly the heat fluxes. Those heat fluxes are needed to solve the moving phase front problem by solving the Stefan condition. As such the 2D heat conduction is solved in a moving solid domain which requires a remeshing-projection algorithm but does not require anymore any closure hypothesis. Such modelling improvement is evaluated by measuring the impact of 0D closures in the lumped parameter models on relevant test cases. Attention is paid on the effective coupling between those truly different models: (multi-layer) stationary closure-based 0D for the liquid, 2D non-stationary heat conduction for the solid and the container. This paper is structured as follows. In section 2, we give a brief presentation of the modelling of liquid/solid phase changes and the associated coupling strategy. Then, in section 3, we propose an improvement of the liquid/solid coupling modelling and an associated numerical resolution strategy based on Finite Element solver, remeshing and projection. In section 4, numerical results are presented in order to verify and validate the new modelling, coupling and implementation. Numerical simulations are used to assess the validity of the overall modelling approach on 1D, 2D academical test cases and more advanced coupled numerical test cases. Finally, concluding remarks and perspectives are drawn in section 5.

Modelling and model coupling

In this section we describe our modelling approach for liquid/solid phase changes, recalling that the corium fusion/solidification is only one possible application and other extensions to different contexts are feasible under minor adaptations. This model is only one part of the overall corium modelling; the way such a phase change model is combined with other models in this context is rapidly described.

Radiative heat transfer

Volumetrically heated

Heat conduction natural convection

Light liquid metal layer

Heavy metal

Vessel

Crust ablation failure Molten oxides pool Fig. 1 Sketch of a nuclear Severe Accident -After melting, the core of the nuclear reactor forms the so-called 'corium', i.e a mixed of molten core oxides and metallic liquids (red/pink color) possibly solidified into a crust (green). Because of the fission products nuclear decay, the liquid corium releases heat volumetrically, further convected and conducted across the reactor pressure vessel (blue), and may lead to the ablation, melting and failure of the vessel.

Ideal multi-dimensional model of phase change

The system of PDEs modelling the physical phenomena is composed of two phases called "liquid" and "solid", which are ideally characterised by two moving 3D domains denoted by Ω l (t) and Ω s (t) and separated by an interface Γ (t) = ∂Ω s (t) ∩ ∂Ω l (t) also referred to as the phase front. The outward pointing unit normal to the boundary ∂Ω is denoted by n. Moreover n Γ represents the unit normal to Γ (t) arbitrarily oriented from the liquid domain to the solid one. The general modelling for the phases is ruled by a set of conservation laws namely for the mass, momentum, and energy as function of the position x and time t. This dependency is implicitly always assumed and omitted for now on. One peculiarity of this work is that we consider that the liquid and solid phases can be made of different materials with different phase change properties. For an homogeneous liquid phase, we make the hypothesis that the conservation laws are written under the Boussinesq approximation in terms of the liquid velocity v and temperature T as the incompressible Navier-Stokes (NS) equations for all

x ∈ Ω l (t) ∈ R 3 and time t ∈ R + ∇. v = 0, (1) ∂ v ∂t + v.∇ v = -1 ρ ∇p + ν∆ v -B(T -T ref ) g, (2)
ρC p ( ∂T ∂t + v • ∇T ) = ∇ • (λ∇T ) + f, (3)
where ρ is the constant fluid mass density, T the temperature, p the pressure, C p is the specific heat capacity, ν the kinematic viscosity, λ the thermal conductivity, B the volumetric thermal expansion coefficient, f is the source term (including residual mass power), g is the gravity and T ref a reference temperature of the liquid.

For the solid phase, we assume that the temperature is governed by the heat conduction equation such that for all position

x ∈ Ω s (t) ∈ R 3 ρC p ∂T ∂t = ∇ • (λ∇T ) + f, (4)
where C p , λ, ρ, ν are material dependent parameters which may be functions of space and time. For instance in the case of a corium like material, the phases are composed of many different species (related to suboxidized metal elements) and their associated parameters usually dependent on the temperature. In this paper, the calculation of these parameters is not detailed and we consider them as given (see [LTHB19] for instance in case of corium). The liquid and solid domains Ω l (t), Ω s (t) are evolving in time and space, separated by the phase change front Γ (t). The boundary conditions (BCs) associated to the liquid and solid domains away from Γ (t) are given for all time t by where T imp , ϕ imp s and ϕ imp l are given and the heat fluxes ϕ s , ϕ l will be specified in the following.

∀ x ∈ ∂Ω l \ Γ (t), v • n = 0, T = T imp or ϕ l • n = ϕ imp l , (5)
∀ x ∈ ∂Ω s \ Γ (t), T = T imp or ϕ s • n = ϕ imp s , (6)
The interface Γ simulates the moving phase front. The model describing its displacement is based on the conservation of mass and energy across this phase front. In this work we consider the so-called Stefan condition [Ste90, Kam61,Gup18] for the front velocity u for all x ∈ Γ (t)

u = u n Γ = 1 ρ Γ ∆H Γ (ϕ l -ϕ s ) n Γ (7)
where ∆H Γ is the latent heat (fusion or solidification) and ρ Γ is the mass density at the interface and u the magnitude of the velocity in the normal direction. Under the local thermal equilibrium assumption at the interface between the solid and liquid phases, the temperature is continuous across the interface in such a way that a unique interface temperature T Γ = T s = T l can be defined for any x ∈ Γ (t) and all t > 0. Heat fluxes ϕ l and ϕ s are given by the Fourier's law as

ϕ k = ϕ k • n Γ , ϕ k = -λ k ∇T k k = l, s , (8)
where λ l , T l (resp. λ s , T s ) are the conductivity and temperature of the liquid (resp. solid) phase in the vicinity of the front where n Γ is considered. Remark that 1 ρ Γ ∆H Γ is always positive during the phase change and defined for all x ∈ Γ (t) by

(ρ Γ , ∆H Γ ) = (ρ s , ∆H fus ) in fusion, (ρ l , ∆H sol ) in solidification, (9)
where ρ s (resp. ρ l ) and ∆H fus (resp. ∆H sol ) are determined according to the composition of solid (resp. liquid) phase at the interface.

According to (7) if the liquid and solid heat fluxes are equal, i.e ϕ s = ϕ l , then the front must remain static and we observe a pure heat conduction situation between the phases. Contrarily in a solidification/fusion situation then ϕ lϕ s is strictly negative/positive. We assume that two extreme temperatures at the interface can be defined depending on the materials' composition: T fus the temperature of fusion of the solid phase, and, T sol the temperature of solidification of the liquid phase. The general case where the solid and liquid can have different compositions is considered. Thus, the solid fusion temperature and the liquid solidification temperature can be different and the condition T fus ≥ T sol is assumed. Then we consider three possible boundary conditions at the front.

The fusion situation where we set T Γ = T fus and ϕ l > ϕ s is assumed, the solidification for which T Γ = T sol is set and ϕ l < ϕ s is assumed, and the pure conduction where ϕ s = ϕ l . In this last situation, encountered when T sol = T fus , the front velocity u is zero and the interface temperature T Γ verifies T sol < T Γ < T fus , see figure 2 for a logical diagram.

Geometrical simplification of the model and use in multiphase corium pool modelling

The model presented in the previous section has no restriction for the dimension nor the geometry of the domain. However in many engineering applications the vessel/container can be assumed to have an axisymmetrical geometry. Moreover, we focus on the phase interactions at a relative lower scale than the full characteristic length. Consequently one often adopts a 2D Cartesian x, z planar domain which partly covers the liquid and the vessel, see figure 3 left panel for a sketch. Doing so the previously described model is maintained alike in a 2D x, z coordinate system with symmetry Boundary Conditions (BCs) at the liquid (West boundary) and Dirichlet BCs at the solids (East boundary) which remain to be defined. The liquid domain has finally an output BC on the top (North boundary). Vessel. The vessel is a metallic solid characterised by its composition and temperature. We assume that only conductive thermal effects play a role within the metal. It could melt if its temperature exceeds the so-called fusion temperature T fus characterising this solid. As a consequence the vessel can exchange mass and heat at the contact with the liquid or solid facing it. Moreover the external cooling of the vessel is considered at given fixed temperature T ves which is always cooler than the fusion temperature of the vessel. Liquid phase. The description of a multi-component liquid is generally complex due to its mixed and multi-scale nature. The conservation equations (1), ( 2) and (3) can only describe the thermal-hydraulics behaviour of a single phase liquid. This is the case when the corium pool is only made of molten oxide materials from the reactor core. For such an homogeneous pool, this model can capture complex convective flow motions (due to the volumetric heat generation associated with decay heat associated with such a low-viscosity fluid) [SLA] and determines the unsteady heat flux distribution on the vessel, impacting its failure or its resistance. When considering that the molten pool can also contains metallic materials, the situation is more complex because the liquid is now composed of two immiscible phases that tends to stratify with gravity. In addition, because of thermochemically induced mass transfers at the interface between these two phases, local modifications of the density can trigger Rayleigh-Taylor instabilities in such a way that the liquid phases spatial distribution is time-dependent. A discussion of modelling approaches for such a multiphase configuration is out of the scope of this paper and the reader is referred to, for instance, [ZL21] for more information and temptative simulations of stratification transients.

Here we drastically simplify the modelling by following a stratification description given in figure 3. Such an idealized configuration is typically considered in lumped parameter models in severe accident codes where the time evolution of the stratification into different layers is described by a simplified inter-layer mass transfer model (e.g. [LTSB15]). This model calculates mass source terms for each layer mass conservation equation that describe in a very simplistic way the Rayleigh-Taylor induced mass transfers. Then, for each layer, an energy conservation equation, obtained by integration of (3) over the layer volume, describes the evolution of the average layer temperature as a function of the heat exchanges through its boundary. Nusselt correlations are typically used to close the system of equations. These different equations are coupled because of interlayer heat exchanges. The reader is referred to [LTSB15, VTP20] for more details. Under such a simplifying modelling, the interface conditions (7) discussed in section 2.1 can still be used for each layer in order to describe fusion/solidification processes; in this case, ϕ l , the liquid heat flux at the liquid/solid interface differs from one layer to the other and is given by the solution of the layer-wise integral energy conservation equations.

Solid phase. When the liquid temperature reaches the solidification temperature, T sol , then a so-called solid crust is generated. The crust composition therefore matches the liquid one at the time of its creation and we assume that it is a solid material. The main characteristics are its temperature and composition, and, only heat conduction is considered. Usually the solid vessel has a lower temperature than the solid crust and has a tendency to cool it down, even if melting could occur. Phase front. The phase front Γ is a 1D curve evolving in the 2D computational domain, see figure 3 left panel. This front motion relies on Stefan condition. Because the solid crust may play the role of protection shield for the vessel, the precise location of this front in space and time is of paramount important. Therefore its multidimensional description and evolution is one key design principle of our modelling.

Although several important simplifications have been made to model the liquid phase for instance, the whole model is based on a multi-dimensional fluid-like model for the liquid phase and a multi-dimensional heat conduction in the solid phases, both coupled and interacting through the interface Stefan condition. Fundamentally it remains somewhat close with respect to the one in section 2.1.

Coupling methodology

In this section we briefly recall the coupling strategy developed into the PROCOR simulation code ("PROpagation de CORium") [VTP20]. PROCOR is an integrated code, composed by multiple multi-dimensional models. Coupling different models and their associated numerical methods is rarely a well-posed problem. Our context is no exception to the rule but we ought to develop a flexible and evolutionary framework [VSDV18]. Although all physical phenomena occur simultaneously, we adopt in this study an explicit time coupling scheme with a coupling through the domain interfaces via boundary conditions (BCs), i.e mass, energy or flux exchanges, see Fig. 4-left for a sketch. The simulation framework is based on models (eventually sub-models), called 'partitions', which are communicating through their own boundary conditions. This framework offers two types of time

T m T m+1 T m+2 t k t k-1 t n t n+1 T m-1 t ∆T ∆t ∆t ∆ t model A model B model C
Fig. 4 Left: Sketch of the coupling methodology. A set of interconnected models is constructed forming the simulation coupling scheme. The models exchange information, i.e mass, energy, temperatures, fluxes, etc. (straight lines). Some links may or may not exist (dashed lines). -Right: Sketch of a possible time dependency between three models (A red, B black and C green). The macro-time-step ∆T determines the meeting points T m to exchange data (vertical plain arrows). During ∆T each model may evolve at its own local time-step ∆t < ∆T employing the most updated available data computed by the neighbor models (colored arrows).

discretisation: explicit and implicit schemes (see [Vio18, VSDV18] for reference). Let us assume that the time interval [0, T ] is discretised into time-steps ∆T = T m+1 -T m > 0. The explicit scheme requires only one call of each solver one after the other during each time step. This type of scheme has a restrictive condition on the time step.

On the contrary the implicit scheme is based on an iterative fixed point method used to solve the problem derived from the model coupling via boundary conditions, and is free from time-step restriction. Nonetheless the multiplication of exchanged variables between the models renders its solve particularly difficult.

The coupling time step ∆T , also called "macro-time" step, defines the different meeting points T m when models exchange information with each others. At the beginning of the macro time step, from the most updated available information from the linked models, each model is solved with its own internal integration scheme and updates its own variables. The restricted time-step of the internal integration scheme of a given model, say ∆t = t kt k-1 , may be smaller than the macro-time-step ∆T , so multiple sub-time-steps may be demanded to a model before reaching the meeting point and the re-coupling with the neighbor models, see Fig. 4-right panel for an illustration.

3 Improved solid/liquid modelling and simulation

In this section we present our modelling and simulation strategy focusing particularly on the solid parts. Indeed the liquid part of the simulation still follows the in-house 0D multi-layer modelling and the coupling and time discretisation uses the explicit coupling described in the previous section. In this work we focus more specifically on the solid phase and its interaction with the liquid one.

2D heat conduction modelling and 2D interface motion

One considers the solid phase Ω s as a whole 2D-(x, z) solid material which form evolves in reaction to the heat flux coming from the liquid phase. Moreover the heat exchange in Ω s core is modelled by the 2D time-dependent heat equation ( 4) subject to initial condition T 0 , temperature BCs on top Γ top and bottom Γ bot , and flux or temperature conditions on the left (liquid side) and right (solid vessel side) boundaries:

ρC p ∂T ∂t -∇ • (λ∇T ) = f, on Ω s (t), T (x, z, t) = T 0 , at t = 0, T | γ ext or ϕ| γ ext , T | Γ or ϕ| Γ , T | Γ top , T | Γ bot given on ∂Ω s (t), ( 10
)
where f is a source term corresponding to a decay heat source in the case of corium. The motion of the phase front Γ is subject to Stefan condition ( 7) recalled here

u = u n Γ = 1 ρ Γ ∆H Γ ϕ liq -ϕ in n Γ . ( 11
)
where ϕ in = λ∇T • n Γ in the vicinity of the front, i.e where n Γ is considered, and, ϕ liq is the input heat flux from the liquid phase. This approach relies on system (10) which is solved on the moving domain Ω s (t). Its motion is defined by the motion of the interface Γ between the phases and Stefan condition (7).

The liquid configuration usually evolves in time, possibly leading numerically to a solid cell facing different liquid layer composition, generating a retro-action between the phases. We numerically ensure in the following that any solid cell in contact with Γ faces one and only one liquid layer. In that way, any smearing effect of the heat flux from the pool is avoided, in particular when a high flux jump occurs between two layers. The evolution of the liquid configuration may generate situations for which a solid cell is not facing any liquid layer anymore then an adiabatic boundary condition is imposed. We denote by P k , 1 ≤ k ≤ N p a liquid layer, see Fig. 3-right. Let us introduce some notation: The solid domain Ω s (t) is a polygonal domain defined on the left (resp. on the right) by a broken line Γ (t) (resp. γ(t) ≡ γ ext (t)) and straight line and fixed boundaries Γ top and Γ bot . It is paved with an unstructured triangular mesh without gap nor hole and subject to BCs at the boundaries and some initial condition at t = 0. We denote by (Ω s ) h (t) its unstructured triangular 2D mesh made of N j > 0 cells. Any given solid triangular cell of index j is referred to as ω j (t) ∈ (Ω s ) h (t). In each cell ω j (t) the temperature is assumed to be constant in space and denoted T j (t). The boundary conditions at the liquid/solid phase front are given by an heat flux ϕ liq = ϕ liq • n Γ , while at the solid container on γ ext by a Dirichlet temperature T ext or prescribed heat flux ϕ ext = ϕ ext • n γ . Moreover the front is defined by a chained list of N i vertical or horizontal segments

I i+1/2 (t) = [ x i (t), x i+1 (t)] with unit normal direction n Γ i+1/2 (t): Γ (t) = I i+1/2 (t) = [ x i (t), x i+1 (t)], x i (t) = (x i (t), z i (t)), 1 ≤ i ≤ N i (t) + 1 ( 12
)
Let us denote the temperature, heat flux and velocity in the normal direction n Γ i+1/2 at each interface segment I i+1/2 by T in i+1/2 , ϕ in i+1/2 and u in i+1/2 respectively. The front points move in x-direction and z-direction with velocity u i (t) = (u i (t), w i (t)) following the trajectory equation

d x i dt = u i , x i (t = 0) = x 0 , (13)
where the velocity u i is deduced from Stefan condition. As already mentioned in figure 2, three possible situations can be encountered for any segment I i+1/2 : fusion, conduction or solidification. In the conduction case we impose the continuity of heat flux without interface motion, while in fusion/solidification cases the equality of temperatures with interface motion are considered. According to the phase change, each segment I i+1/2 of Γ is displaced. If this interface is in the fusion or solidification state but with (ϕ liqϕ in i+1/2 ) cancelling, then the phase front is stopped and the interface state transitions towards a conduction state (see Table 1). On the contrary, if the interface is in the conduction state and its temperature reaches the fusion or solidification temperature (T in i+1/2 = T sol or T fus ), then the interface state adapts accordingly. The algorithm at t n , the beginning of a time-step, considers the following data:

States at phase front Γ State: Fusion ←→ Conduction ←→ Solidification Interface velocity: u i+1/2 . n Γ i+1/2 > 0 u i+1/2 . n Γ i+1/2 = 0 u i+1/2 . n Γ i+1/2 < 0 Interface temperature: T in i+1/2 = T fus T sol < T in i+1/2 < T fus T in i+1/2 = T sol Heat flux: ϕ in i+1/2 < ϕ liq ϕ in i+1/2 = ϕ liq ϕ in i+1/2 > ϕ liq
Table 1 Summarise of the phase changes at the liquid/solid interface. a discrete interface Γ n , i.e interface points x n i and segments I n i+1/2 ; -a solid domain Ω n s , its associated triangular mesh (Ω s ) n h and associated BCs; -an associated piece-wise constant temperature field (T n j ) 1≤j≤Nj . The algorithm proceeds as: solve the heat equation ( 10) for ∆t to get the temperature and heat fluxes (Solve step), then the front is advanced with (11) to get Γ n+1 (Lagrangian interface step) which implicitly updates the solid domain Ω n+1 s , next a new mesh (Ω s ) n+1 h is constructed (Remesh step), and, finally one projects the temperature field T n+1 j from (Ω s ) n h onto (Ω s ) n+1 h (Projection step). The time-step is then complete. Let us describe these steps in details in the following sections.

0 1 1 • T • T × ϕ × ϕ × ϕ × ϕ × ϕ × ϕ

Solve step -Mixed Raviart-Thomas finite element discretization

To solve system (10) we rely on the mixed Raviart-Thomas Finite Element (FE) discretization [RT77]. More precisely we consider Raviart-Thomas element of zero order (RT0) over the fixed unstructured triangular 2D mesh (Ω s ) n h paving the solid domain Ω n s at time t n subject to BCs at boundaries Γ s and initial condition at t = 0. The mixed FE method is a formulation such that the problem discretization results in a linear algebraic system of the form

A B t B 0 T ϕ = f g , (14)
with A, B matrices and T , ϕ, f , g vectors. It is an alternative method which combines some advantages of finite volume and finite element methods, for instance the method is locally conservative, the fluxes are continuous between neighbour cells, unstructured grids can be handled, an-isotropic discontinuous conductivity are consistently dealt with. Here unknown vectors T and ϕ refer respectively to the degrees of freedom of the discrete temperature and heat flux on the mesh. While discrete temperature is located at cell center, the normal heat flux is at face center. A such the mixed FE choice explicitly computes the heat fluxes at each face, in particular those ones constituting the interface front Γ . Thus, the displacement of each interface segment is computed by solving its associated Stefan problem (see section 3.3) without any artificially heat flux interpolation. We do not detail this genuinely classical mixed FE method and its properties but only refer to one text book [AdVBL17] (although a vast literature is available on the topic). Dirichlet type of boundary conditions are considered in the fusion/solidification states, and Neumann type boundary conditions in the conduction state. As such the mixed FE method provides the temperature (T n+1 j ) j=1,...,N n j and heat flux vectors at the cell faces, (j, f ):

ϕ n+1 j = ϕ n+1 j,f f =1,2,3 = (λ∇T ) n+1 j,f • n j,f f =1,2,3
, where n j,f is the outward pointing unit normal to face (j, f ), at discrete time t n+1 in any of the N n j cell ω n j ∈ (Ω s ) n h ( n j,f =n Γ j on faces laying on Γ ). This finite element problem is solved explicitly under the CFL condition of the form ∆t solver ≤ CFL ∆x 2 2D , where CFL < 1 is a safety constant, and D a measure of the maximal diffusion coefficient. In this work we have used the open-source computing platform FEniCS for solving the partial differential equations [ABH + 15, LMW + 12].

2D Lagrangian interface motion step -Solidification, conduction or fusion

The motion of phase front Γ is driven by equation ( 11) which determines the front velocity in the normal direction n Γ i+1/2 of each segment I n i+1/2 (coinciding to a face f of a solid cell ω j ) during the time-step ∆t:

u n+1 i+1/2 = u n+1 i+1/2 n Γ i+1/2 = 1 ρ n+1 i+1/2 ∆H n+1 i+1/2 ϕ liq,n+1 k -ϕ in,n+1 i+1/2 n Γ i+1/2 , ( 15
)
where ϕ in,n+1 i+1/2 is the heat flux computed in section 3.2 for the face f of the solid cell ω j coinciding with

I n i+1/2 , that is ϕ in,n+1
i+1/2 ≡ ϕ n+1 j,f and n Γ i+1/2 =n j,f the normal at face f on Γ . Moreover ϕ liq,n+1 k is the heat flux computed by the liquid model for the liquid layer k. Next ρ n+1 i+1/2 and ∆H n+1 i+1/2 are the interface parameters computed depending on its phase state:

ρ n+1 i+1/2 = ρ n+1 j , ∆H n+1 i+1/2 = ∆H fus j if ϕ liq,n+1 k -ϕ in,n+1 i+1/2 > 0, ( 16
)
ρ n+1 i+1/2 = ρ sol k , ∆H n+1 i+1/2 = ∆H sol k if ϕ liq,n+1 k -ϕ in,n+1 i+1/2 < 0. ( 17
)
Where ρ sol k corresponds to the facing liquid layer density. ∆H fus j (resp. ∆H sol k ) is determined according to the composition of the corresponding cell solid (resp. layer liquid). The velocity of the front from (15) allows to compute the averaged mass flux across this interface segment during ∆t. Once the segment velocity u n+1 i+1/2 is determined, then, the horizontal (resp. vertical) segment I i+1/2 is translated vertically (resp. horizontally) according to the motion of its center:

x n+1 i+1/2 = x n i+1/2 + ∆t u n+1 i+1/2 . ( 18
)
Notice that the segment velocity can be zero if ϕ in,n+1 i+1/2 = ϕ liq,n+1 k , meaning that this interface segment is indeed locally static for the current time-step, equivalently the interface segment is in pure conduction state. The displacement of the N i segments from t n to t n+1 leads to a new geometric discretization of the interface Γ n+1 based on a new chain list of segments, see figure 6. Because the size of any interface segment I i+1/2 may become arbitrarily small during the phase change process, a smallest allowed segment size is assumed with parameter δx = δz > 0, depending on a characteristics mesh length. It implies that Γ is captured up to δx, and smaller segments are discarded. This small parameter also avoids that the CFL condition becomes artificially restricted. This interface displacement method needs a restrictive time step condition which reads as

∆t I ≤ CFL ∆x max i u n i+1/2 .

Remesh step

The remeshing step consists in determining a mesh of the solid phase given the updated interface segments I n+1 i+1/2 ∈ Γ and the boundaries Γ top , Γ bot and γ. The updated interface segments I n+1 i+1/2 provide an implicit coarse meshing of Ω n+1 s by extending them in x or z directions, see figure 6-left. Their intersections become new mesh points, which are further gathered into large quadrangular mesh cells. If those cells are larger than user-given characteristics lengths ∆x > 0 and ∆z > 0, then they are further split into four sub-quadrangles. The obtained mesh is therefore an unstructured mesh of Ω n+1 s made of quadrangles. In order to apply the mixed Raviart-Thomas spatial scheme, at last each of these quadrangles is split into two triangles by their top-left/bottom-right diagonal and the obtained mesh made of triangles ω n+1 j is referred to as (Ω s ) n+1

Γ = I i+1/2 = [ x i , x i+1 ], s i = (x i , z i ), 1 ≤ i ≤ N i + 1 .
δx and δz are the smallest allowed measures for a segment/face/cell -Right: Illustration of the assumption of instantaneous quadratic temperature profile T j within the crust layer C j knowing boundary and mean temperatures.

Projection step

The projection operator is made as simple as possible and consists of a first-order interpolation of the temperature field from the triangular mesh of the solid phase at t n onto the new one at t n+1 . If, the center of a new cell ω n+1 k is located in the old domain Ω n s , then an interpolation using RT0 finite element shape functions is made to compute the temperature T n+1 k at the new cell center. Contrarily, if the cell center of ω n+1 k is outside the old domain then its temperature is uniquely associated to the temperature of the closest old edge of Γ n of (Ω s ) n h . This latter configuration can occur only in solidification regime and, therefore the solidification temperature T sol is set at this cell. This relative crude algorithm allows to remap the old temperature onto the new mesh consistently and without any ambiguity nor further approximation. Notice that more advanced algorithms could be considered: secondorder accurate ones with profile reconstruction, conservative methods for specific enthalpy, etc. [LS05, KSW03].

3.6 Summary of the algorithm -Solve, Interface, Remesh and Projection As a summary, the solid domain Ω n+1 s is modified following the phase front Γ n+1 such that in general Ω n+1 s = Ω n s . Therefore a remesh step must provide (Ω s ) n+1 h , followed by a projection step which transfers the temperature profile (T n+1 j ) j=1,...,N n j defined on (Ω s ) n h onto the new mesh (Ω s ) n+1 h . Index j is reserved for the mesh cell at time t n , while k for t n+1 cells determined by the new interface Γ n+1 . Therefore starting at t n from cell centered temperatures T n j defined on the cell mesh ω n j ∈ (Ω s ) n h implicitly defined by front Γ n , our algorithm relies on the succession of four operators: a FE Solve, an Interface motion, a Remesh and a Projection. They can be recast into:

S : (Ω s ) n h , (T n j ) j=1,...,N n j , BCs, -→ (T n+1 j , ( ϕ n+1 j ) j=1,...,N n j ) (19)
I : Γ n , ( ϕ j ) j=1,...,N n j , BCs, -→ Γ n+1 (20) R : Ω n s , (Ω s ) n h , Γ n+1 -→ Ω n+1 s , (Ω s ) n+1 h , (21)
P : (Ω s ) n+1 h , (Ω s ) n h , (T n+1 j ) j=1,...,N n j -→ (T n+1 k ) k=1,...,N n+1 k , (22)
where T n+1 k are defined on the mesh (Ω s ) n+1 h at t n+1 implicitly defined by front Γ n+1 . With the notation 'BCs' we include the data emanating from the interacting models (liquid, solid container, true boundary conditions, etc.), see Fig. 4-left.

As already mentioned, we re-employ the explicit time discretization provided by the simulation platform PROCOR, see section 2.3 and [Vio18, VTP20] for more details. So, during ∆t we solve sequentially the models: S → I → R → P → BCs. In order to ensure mass conservation during the interface displacement, a directional splitting is implemented. Only vertical displacements are taken into account during a first resolution, then the horizontal ones are dealt with. A switch between these two directional steps is done from one time step to the other to reduce artificial mesh imprint. The time step restriction is built upon the restriction from the FE solver and the interface model supplemented with a coupling constrain imposed by the explicit scheme, which reads

∆t S ≤ CFL ∆x ∆z 2λ , ∆t I ≤ CFL ∆x max i u n i+1/2 , ( 23
)
where CFL = 1/2 is our safety parameter. At last the time-step is determined by taking the minimal value ∆t = min(∆t I , ∆t S ). This complete the description of our solver.

3.7 An alternative approach based on 0D modelling of phase change

In this section we briefly present the numerical strategy developed into PROCOR simulation code. In the 0D modelling philosophy the liquid and solid are assumed to be stratified along the vertical direction (z axis). Each layer of the solid is a rectangle aligned in the x-direction, and we denote by V j its surface. Likewise for the 2D description, any solid layer is assumed to face a unique liquid layer. A solid layer C j is entirely determined by its state vector composed by its mass m j = ρ j V j where ρ j is the density, and its averaged temperature T j . Some other given constitutive quantities characterise the type of material, these are the thermal conductivity λ j , T sol j (resp. T fus j ) and ∆H sol j (resp. ∆H fus j ) the solidification (resp. fusion) temperature and latent heat of solidification (resp. fusion) of the liquid facing C j (resp. solid in C j ). In this 0D modelling, only the vertical segments of C j in contact with phase front Γ , denoted Γ j , are moving horizontally following the Stefan condition. As such vertical displacements are neglected. Moreover we define the averaged conductive heat fluxes associated to the boundaries of C j with Γ and γ ext respectively as

ϕ in j (t) = 1 ∆z j Γj (t) -λ j ∂ T j ∂x (z, t)dz, ϕ ext j (t) = 1 ∆z j γ ext j (t) -λ j ∂ T j ∂x (z, t)dz, (24)
where ∆z j is the vertical space discretization and T j is a continuous approximation of the temperature in the solid layer C j which will be determined below. The conservation laws applied to any solid layer C j govern the evolution of its mass m j and temperature T j as

dm j dt = -ṁΓ j , (25)
C p m j dT j dt + ṁΓ j (T fus -T j ) = ∆z j ϕ in j -ϕ ext j + f, (26)
where ṁΓ j is the mass flow rate subsequent to the 1D horizontal displacement of the phase front. Written as such the solid model neglects the possible axial conduction between two adjacent solid layers, even if some more advanced models do so [PLTV + 18]. In the same way, axial fluxes imposed by the liquid phase on the horizontal segments of Γ are neglected. As already mentioned three possible situations can be encountered for any solid layer C j : fusion, conduction or solidification. The mass flow rate ṁΓ j on Γ in j is determined by those transitions which can be summarized here by the following conditions:

Fusion: T in j = T fus j =⇒ ∆H fus j ṁΓ j = ∆z j ϕ liq j -ϕ in j , T in j = T fus j , (27)
Conduction:

T sol j < T in j < T fus j =⇒ ṁΓ j = 0, ϕ in j = ϕ liq j , (28)
Solidification:

T in j = T sol j =⇒ ∆H sol j ṁΓ j = ∆z j ϕ liq j -ϕ in j , T in j = T sol j . (29)
Each vertical segment Γ j is displaced in the x-direction for a fusion or solidification situation. On the contrary if ṁΓ j cancels, then the phase front is stopped and the cell transitions towards a pure conduction state. Notice that the value of ϕ liq j in (27-29) is set by the liquid 0D multi-layer model. The remaining entities ϕ in j and ϕ ext j (24) are the heat fluxes at the liquid and external/solid boundaries respectively and they depend on the form of T j . In respect with a 0D modelling, a quadratic temperature profile in x direction (independent of z) within any solid layer is assumed, i.e. the form of T j is T j (x) = ax 2 + bx + c, [TSS17]. The parameters a, b, c are computed in order to match the averaged solid temperature T j and the two given flux/temperature boundary conditions at Γ and γ ext . From this "ad-hoc" definition of T j , we compute the value of the heat fluxes (24) which then close the 0D solid model (25)(26). Although well fitted to PROCOR we can point several drawbacks of the 0D modelling approach, namely: Multi-layer solid description. The horizontal size of the solid layer evolves solely according to the facing liquid layer at Γ j , and, the contact with the vessel at γ ext j . No displacement is considered for horizontal inter-layer segments. The vertical size ∆z j of a layer is fixed when created. Therefore C j has no mass or heat exchanges with the top and bottom solid layers (nor through boundary conditions). Although this assumption is valid for mass transfer between solid layers, it is obviously a crude approximation that completely neglects axial heat conduction; 0D modelling of solid layer via a quadratic temperature profile. The integral solid layer model instantaneously averages the boundary effects from the liquid/container along the entire solid layer. As such there can be no transient heat conduction inside a solid layer. Moreover the assumption that the temperature T j adopts a quadratic profile is not adequate in transient regimes, although it is certainly valid for steady state solution.

With such a 0D model closed with quadratic temperature profile, see figure 6-right, the multi-dimensional heat conduction is not genuinely solved in the solid phase.

As a consequence the quality of the phase front Γ may suffer from those before-mentioned approximations although its precise description and evolution is generally of paramount importance.

In the numerical section we will use PROCOR as a possible simulation code to compare with in the specific context of SA scenarios.

Numerical results

This numerical section presents validation tests for the new interface modelling and solid/liquid coupling as described in the previous sections. The methodology of testing first focuses on the sanity validations in 1D and 2D to check the implementation and modelling assumptions. Secondly, simplified test cases are simulated and the results are compared to the legacy PROCOR simulations; the goal being to numerically prove that our approach brings some benefits in certain configurations. At last we present a more representative physical situation for which the coupling with an evolving liquid phase is done. Notice that, apart from the sanity tests and a genuine 2D test case, we discard the axial displacement of the phase front to allow a fair comparison with PROCOR simulation tool. A full 2D test case is used to prove that our approach can handle such axial phase front displacements and, as a consequence, can be associated to a genuine multi-dimensional modelling of the liquid phase.

Sanity check: 1D Stefan problem

In order to validate our approach we use the so-called 1D Stefan problem. This one dimensional problem possesses an analytical solution [Ste90, Kam61] to which we can observe the convergence of the the spacetime evolving numerical solution as the mesh is refined. This test case simulates the 1D solidification (resp. fusion) of an homogeneous material, such that the solid and liquid phases have the same material properties. The computational domain corresponds to the solid Ω s = [0, s(t)]. A liquid domain (not computed) stands on the right of this domain: Ω l (t) = [s(t), +∞] at the constant solidification temperature T sol , where s(t) is the location of the moving front in-between at any time t > 0. An external temperature T ext < T sol is enforced at the left boundary, and T sol at the right one. The aim is to compute the position of the interface s(t) and the temperature profile T (x, t) in the solid state.

The governing equations of the 1D Stefan problem are:

∂T ∂t (x, t) = α ∂ 2 T ∂x 2 (x, t), β ds dt (s(t), t) = -∇T (s(t), t) , (30)
where α = λ ρC p , β = ρ∆H sol λ , and x ∈ Ω s = [0, s(t)], t ∈ [0, t final ]. The boundary and initial condition are given by: T (0, t) = T ext , T (s(t), t) = T sol , and s(t 0 ) = 0. The unique solution of this problem is, see [Ste90]:

T (x, t) = T ext + (T sol -T ext ) erf x 2 √ αt erf(k) , (31)
s(t) = 2 k √ αt, ke k 2 erf(k) = 1 βαπ T sol -T ext , ( 32
)
where erf is the error function. k is approximated by Newton's method and further substituted into (31) and ( 32). For our simulation we choose t 0 = 0.16256, α = 1.0, β = 1.0, T sol = 0.0 and that k 0.62 leading to an initial front position at s(t 0 ) = 0.5. The simulation runs with our 2D code for any (x, z) ∈ Ω = [0, s(t 0 )] × [0, 1], z-direction being the ignored direction. The mesh is constituted of N x × N z quadrangular cells each split into two triangles, and N t timessteps are employed to reach the final time leading to ∆x = 1 Nx , ∆z = 1 Nz and t final = t 0 +N t ∆t. Here we employ ∆z = 1. Initially the solid cells ω j are such that x j + ∆x 2 < s(t 0 ), all others are liquid ones. The solidification front moves in the right direction and is located at s(t) for all time t > t 0 , see figure 7-top. As the front moves, some cells of size ∆x transition from the liquid to the solid state during the re-meshing step. The final time is set to t final = 100. In figure 7 we plot the position of the front s(t) as a function of time t for several ∆t and ∆x (coloured symbols) against the exact solution in black straight line. The numerical solution almost perfectly reproduces the exact one when time and space step decrease. The temperature profiles are quasi-linear ones. As expected the differences between the approximate and exact solutions are minor and we omit these figures. Next we present in table 2 the errors and orders of convergence in L 1 , L 2 and L ∞ norms for the interface position s(t) to assess that the scheme is converging to the exact solution. This 1D sanity check shows that our approach is capable of capturing the solution of Stefan 1D problem. A liquid material at temperature T liq is filling Ω = R 2 /B at t = 0 where B = (x, y) ∈ R 2 , x ≤ 1/20, or y ≤ 1/20 at t = 0. The liquid solidifies in presence of a cold solid in B initially at T sol cooled down from the semi-infinite axes which are maintained below the solidification temperature, i.e T ext < T sol . As time increases, 0 < t 1 < t 2 < t 3 < t 4 , the phase front adopts a 2D curved shape, while, away from the origin, the 1D Stefan theory applies.

T ext = -1.0, which implies ∆t L 1 error L 1 order L 2 error L 2 order L ∞ error L ∞ order 1.0 2.87 × 10 -1 - 3.11 × 10 -1 - 8.33 × 10 -1 - 0.

Sanity check: Solidification generated by an infinite wedge shape

In this section we ought to validate the 2D implementation by considering a classical 2D test case, for which there exist no exact but only approximate solutions [RA71,BUD73]. This problem features a liquid material at temperature T sol filling Ω = R 2,+ /B at t = 0 where B = (x, y) ∈ R 2,+ x ≤ 1/20, or y ≤ 1/20 . Domain B is a solid initially at solidification temperature T sol . An external temperature T ext < T sol is enforced at the wedge x = 0 and y = 0, as such it gradually cools down and freezes the liquid phase, see Fig. 8-left for an illustration. The evolving solidification front has a 2D shape close to the corner, while, along the axes away from the origin, the classical 1D theory applies, i.e the solidification fron10.1098/rsta.2020.0130t s(t) has a well determined 1D motion parallel to the axe, see the previous test case. In this test case we are interested in the simulation of the 2D phase front shape as a function of time. We consider the specific boundary conditions T ext = -1 for the semi-infinite axes, and, the initial condition T sol = 0 for the solid at time t = 0 and the liquid phase is also considered at the solidification temperature T sol . The simulation is run up to time t = 2.4 × 10 -3 by increments of size δt = 0.2 × 10 -3 . All other physical parameters, like the conduction, are fixed as the 1D Stefan test case. In Fig. 9 we present the numerical results obtained by our 2D modelling approach. On the left panel we present the phase front in-between the solid/liquid materials for different times t = 0.4k × 10 -3 for k = 1, 2, 3, 4, 5, 6 and a zoom on the corner location. Notice that the solid domain increases in size (number of cells and surface) as the solidification front freezes the liquid. The phase front presents a perfect symmetry with respect to the y = x axe as expected. Moreover the location of the phase front and the temperature profile away from the corner match the exact 1D Stefan solution. On the right panel of Fig. 9 is plotted the temperature profile in the solid domain at time t = 2.4 × 10 -3 . We verify that T ext ≤ T ≤ T sol and observe true multidimensional effects close to the wedge, while the expected 1D temperature profile is retrieved away from it. At last in this panel the unstructured triangular mesh is plotted, and we observe a clear adaptation to the 2D phase front shape. Notice that the 2D phase front shape in the vicinity of the wedge is in accordance with other published results [WG20, YG10, RA71].

Academic test: an idealised example of solid fusion

In this section we consider yet another idealised test case for the fusion of a solid vessel, see figure 10 for a sketch of the configuration and [PLTV + 18]. This test is built from a typical in-vessel configuration of interest for light water reactor severe accident analysis. In such a case, the maximum effect of axial conduction is expected to be at the top of the liquid (corium) where the flux onto the vessel undergoes its largest discontinuity. We build a simplified configuration reproducing such a flux discontinuity here. A rectangular domain Ω = [0, 1] × [0, 10] is assumed to represent the solid initially at constant temperature T 0 . We define a part of the boundary as Γ * = {(x, z), s.t. x = 0, 0 ≤ z ≤ z * }, with z * = 2, and, γ ext = {(x, z), s.t. x = 1, 0 ≤ z ≤ 10}. Adiabatic BCs are considered except on Γ * where we impose a constant positive heat flux ϕ * , and, on Γ ext where a Dirichlet cooling boundary condition T 0 = T ext is imposed. Moreover we place five gauges, A, B, C, D, E, on γ ext at position z A < z B < z C < z D < z E (resp. at 0.1 < 1.1 < 1.9 < 2.1 < 4.1). Because of the positive heat flux at Γ * we expect the diffusion of heat in the solid, possibly coupled with its fusion if the local cell temperature exceeds the fusion temperature, i.e T j > T fus . The computational mesh is constituted of N x × N z × 2 cells. A quadrangular mesh made of N x ×N z cells is first built with cell size ∆x×∆z and each quadrangle is further split by one of its diagonal into two triangles. A "coarse" mesh employs ∆x = 0.02 and ∆z = 0.37. An "intermediate" (resp. "fine) mesh is then built considering ∆x/M and ∆/M with M = 2 (resp.M = 4). N t times-steps are employed to reach the final time t final > 0. The non-dimensional equations can be parametrised by their Biot number Bi spanning different physical regimes, see Appendix A for details. Here the Biot number adopts the values Bi = 0.8 and 1.2, which are representative of two physical behaviours: for Bi = 0.8 a pure conduction regime is sustained, contrarily for Bi > 1 fusion must occur. For the sake of simplicity, we consider T ext = T 0 = 0 and T fus = 1. In the following we compare our 2D modelling results against the classical 0D results from the in-house code PROCOR.

Conduction situation: Bi = 0.8 The comparison between the models only highlights the 2D effect on the conduction, see figure 11 where we present the temperature profiles obtained using the intermediate mesh size ∆x = 0.01 up to final time t final = 4 (three times are presented). We observe genuine 2D effects because the cells above position z * are clearly heated, while it is not observed for the 0D model (right panels). Also at time t = 0.03 the 0D model produces nonphysical negative temperatures (purple color), a clear disadvantage of the quadratic profile temperature hypothesis. Contrarily the 2D model maintains the physical admissibility of the temperature. At late time, beyond t = 1 a steady-state solution is attained by both models. While different, in particular over heated on the Γ * boundary and unchanged above z * for the 0D model, they are in qualitative agreement. Notice that the results with coarser and finer meshes present the same general behaviours therefore we do not present the figures here. Next figure 12-top depicts the outgoing heat flux on γ ext at the five gauges Fusion: Bi = 1.2 In figure 13 we present the temperature profile and the fusion front for the 2D and 0D models for three times: t = 0.025, 1, and 4. At the beginning of the transient (left panel) only the heat conduction takes place and no fusion has yet occurred. On the contrary at intermediate and final times, the fusion is clearly visible as well as the heat conduction. At later time, the steady-state solution is attained. At early time (left panel) the 0D model again produces non-physical negative temperatures visible (purple colour), while our approach is free from such an artefact. This is again a drawback of the 0D model linked to the quadratic profile hypothesis. For later times, we can see on the 2D model results that, consistently with the Bi = 0.8 situation, the thermal energy diffuses towards the top layers beyond z * . This effect can not be simulated by the 0D model as observed on the right panels. As a consequence the temperature of the top layer remains constant equal to 0 due to the multi-layer 0D modelling assumption which is in violation with the physical continuity of the temperature. Obviously the shape of the fusion front is drastically different for the two models, in particular for the 2D model the 2D diffusion of heat leads to a slower phase front velocity. For a more quantitative comparison we propose in figure 14 Even for z ≥ 2.8 the temperature is certainly non zero. This justifies the use of genuine multidimensional heat conduction for the solid phase modelling. Next we present in figure 16 the time evolution of the ratio between the initial and final mass of a solid x-aligned layer at gauge A, B and C. The 0D model results are in straight red line while in blue symbols for the 2D model. The mass ratio decreases when the fusion starts and the two models already disagree on when it occurs. Moreover while the 0D model predicts the same evolution for the three mass layers reaching about 26% of melted mass, the 2D models predicts a different evolution and final steady states for each layer; 24%, 16% and 0% for gauges A, B and C respectively. In term of total melted mass, the 0D model ablated 236% more than the 2D one while the minimal residual solid thickness at steady-state is 0.86 (resp. 0.83) for the 2D (resp. 0D model). Notice that this idealised fusion test case has been simulated with smaller and larger Biot numbers leading to the same qualitative conclusions.

Academical test 2: an idealised example of vessel fusion with 2D ablation

On the previous tests the axial displacements of the phase front was not allowed on purpose to allow a fair comparison with the 0D modelling in PROCOR code. In other words our 2D model was restricted in its capabilities, the ablation/melting was occurring only in the horizontal direction. Here we allow such 2D ablation and re-run the previous solid fusion test. The goal is to observe if this no-axial motion simplification assumption has a macroscopic impact on the phase front evolution and on the steady state solution. In order to provide a fair comparison with the previous test, this test case is set up by considering an input thermal power on the whole liquid/solid interface Γ * equal to the one of the previous test case. Therefore, a constant input thermal power P * is assumed during this test case deduced from the previous one. Indeed, in the last case, the input heat flux ϕ * was applied on a constant front length equal to l * = (z *z min ) = z * (the flux was injected only through the vertical segments), and so P * was constant. Here, horizontal and vertical displacements being allowed, the phase front length l * under azimuth z * does vary in time. To ensure that the same heat power is always transferred, the heat power P * is considered constant and uniformly distributed on all phase font segments under z * . As a consequence the heat flux ϕ * might vary during time in function of interface shape. In figure 17 we present the temperature profile and phase front location for the mesh size ∆x = 0.05 and for two different times t = 1 and t = 2. Although qualitatively the phase fronts are in agreement, they nonetheless have some differences, justifying that axial displacement of the phase front may be yet another important modelling feature to take into account. However this feature is directly linked to liquid modelling because this later must provide the heat flux along all segments of Γ * . Currently in the PROCOR environment this capability is lacking, hence our hypothesis concerning the heat power to be injected.

More advanced test case in relation with an industrial like application

The last section has illustrated some improvements brought by the 2D models on simplified situations. In this section we present a more advanced test case in relation with an industrial like application to simulate an IVR strategy. In this case the liquid, solid and vessel models are computing and coupled all-together, retro-acting with each others leading to a genuine fully integrated simulation.

Physical situation -Test case presentation

The physical situation and simulation setup is thoroughly described in [VTP20], however we recall the main steps. The purpose of this test is to simulate the crust formation between the vessel and the liquid corium during an IVR strategy and observe its impact on the thermal load onto the vessel. The crust is a solid obtained from the solidification of the liquid corium. A special attention is paid to its composition during the solidification process, mainly made of oxides (ZrO 2 , U O 2 , F eO) and metal (F e), in order to differentiate the oxide crust with a lower conductivity from the metallic one. Indeed, due to the transient stratification of the liquid corium, an oxide crust could face a metallic corium layer and protect the vessel from the possible heat flux focusing effect taking place in the corium metallic layer. Vice-versa a metallic crust may be in contact with liquid oxide corium. The test is designed with the typical geometry of a light water reactor. The vessel is a spherical cap with a radius and a height of 2 m and made of steel. It contains the liquid corium heated in volume with a total power assumed to be constant at 14 MW. The external vessel wall is maintained at 1600 K by an external cooling and we distinguish two specific solidification/fusion temperatures: 3000 K for the oxide material and, 1600 K for the metallic one. This configuration is designed to challenge the new modelling against complex and various situations encountered in severe accident scenarios. Both short-and long-term behaviours occur in a single calculation; 1. Firstly, the initialisation is made with a pure oxide liquid corium of 82.5 tons at temperature of 2900 K. At this stage the formation of an oxide crust between the vessel and the liquid corium (up to its stabilisation) is expected. 2. In a second part, starting at t = 6000 s, molten steel at 1600 K is added on top of the liquid corium at constant mass flow rate of 10 kg/s during 2000 s. Because of thermo-chemical effects coupled with gravityinduced liquid separation, this addition of molten steel leads to a three-layer liquid stratification: light/heavy metal and oxide layers. 3. In the last stage, starting at t = 10000 s, the total power is reduced to 0.1 MW and the external temperature of the vessel drops to 1200 K, to emulate the long term conditions when the system cools down and solidifies.

Numerical set up and results

The coupling time step in the explicit time discretization, ∆t coupl. , is equal to 10 s and further reduced to 5 s when highly discontinuous phenomena are expected. We consider the multi-layer stratified liquid 0D model for the corium pool (presented in section 3). Both 0D and 2D models presented in section 3 for the crust are considered and compared here. The 2D crust mesh matches the vertical discretization of the 0D model, i.e 50 vertical cells, see section 3.7. A 2D mesh paving the 0D crust domain is built.

The main observable diagnostics are the thickness of the crust, the thermal load imposed by the crust onto the vessel, and the temperature profile in the crust. Some of these variables are compared with the results of PROCOR which uses the 0D thermal model for the crust presented in corium liquid computed by our 2D simulation coupling scheme. The colours correspond to the mass fraction of the main oxide species, red colour indicating an oxide composition. In the first stage the liquid and solid corium are solely composed of oxide materials. In the second stage, in the middle panel, the liquid corium is stratified in three layers occurring after the metallic liquid input. The top and bottom liquid layers do not contain any oxide species, they are composed essentially with metallic materials (blue colour). The top crust is entirely melted and one oxide crust is facing a metallic liquid layer. In the last stage the crust has grown due to the solidification induced by the lowering of the internal power. On the top we can observe the creation of a solidified metallic crust. These results are in qualitatively in agreement with the in-house code PROCOR results, see [VTP20] Figure 19 presents the total mass of each layer of the liquid corium and the solidified corium crust mass. It presents both, the results for the 0D model and for the 2D model. The 0D and 2D models are in rather good agreement with some slight differences however. The right panel of figure 19 represents the difference of the total solidified crust mass between the 0D and 2D models. The 2D crust is heavier than the 0D one (up to 1400 kg, i.e about 5 to 10 percents). The oscillations observed during the second stage are explained by a model reaction delay to the occurrence of the stratification. Many more diagnostics can be used to point out the differences between the 2D and 0D modelling approaches[VTP20], but we omit them here. This industrial-like test case shows that the new 2D modelling of the solid material (corium crust) can be coupled and retro-act with a 0D multi-layer stratified liquid model (corium pool) and a metallic vessel. Our approach can reproduce the general behaviours obtained by the legacy 0D in-house model except for some localised but important differences in terms of crust thickness (and heat fluxes, not shown here). Those differences are mainly due to the conduction in the vertical direction between two horizontal layers which is taken into account with our 2D model. This is particularly the case when two crust cells, in solidification state, are facing two different liquid layers having different solidification temperatures, typically metal and oxide liquid layers. In this configuration, the two crust cells have two drastically different boundary temperatures and the 2D diffusion approach becomes essential to dissipate the high temperature gradient in-between these cells.

Conclusion and perspectives

In this paper we have presented several improvements for the modelling and simulation of a phase change front in-between a "hot" liquid and a "cold" solid. Such situation occurs in many diverse physical and engineering applications. While our work is not restricted to the In-Vessel-Retention application scenario for nuclear severe accident situations, our target application is the simulation of the phase changes of a molten core oxide and metallic liquid material (corium). The improvements are brought into the simulation code PROCOR dedicated to compute the interaction between a hot liquid corium pool and its container, an externally cooled reactor metallic vessel, relying on a multi-layer 0D solid model supplemented with several simplification assumptions (independence of the layers, neglecting the transient heat conduction, ad-hoc quadratic temperature profile). In this work we consider a 2D heat conduction model for the solid phase coupled to a liquid phase via 2D interface motion driven by Stefan condition. As such the solid domain is evolving in time and 2D space as a function of the phase changes, this requires a remeshing procedure and a subsequent projection step. The 2D heat equation is solved with a numerically suited mixed finite element method and the solid domain is determined via the resolution of Stefan condition on the edges at the solid/liquid contact. The 2D solid model is semi explicitly coupled with the in-house liquid and solid vessel/container models.

We have performed several test cases to validate the approach and to point out several drawbacks of the original 0D modelling. (Notice that such modelling is available in several classical severe accident scenario codes.) First sanity checks involving only fusion/solidification fronts are simulated showing that the 2D heat conduction approach can capture the exact or analytical phase front for the 1D and 2D Stefan problem. Secondly, we have run several idealised 2D fusion situations to emphasise the new simulation capabilities brought by our approach: 2D heat conduction, improved front shape and appropriate physical behaviours (positive temperature, reasonable transient heat flux, etc.). We have observed that the 2D approach improves the prediction capability of the whole simulation coupling scheme when supplemented with a 0D liquid model and external boundary conditions. In a last test we have proposed a more industrial-like situation involving a liquid corium which composition and stratification evolve in time. The corium solidified crust then reacts to the liquid corium evolution and viceversa. We have compared the multi-layer 0D and 2D models and observed that the new approach produces close but different results. More work is now required to assess that the proposed improvements are well suited within the whole simulation framework by a genuine verification campaign. Finally, a coupling with an incompressible Navier-Stokes simulation solver is in progress to improve the modelling of the liquid phase and the subsequent interaction with the solids.

A Dimensionless equations

In . As a consequence we have θ fus = 1 and θ ext = 0. We observe that if θ ∈ [θ ext , θ fus ] reaches the value θ fus = 1, the crust ablation. We denote by φ the normal component of Φ that we will also call as 'heat-flux'. The last remaining free parameter is the Biot number Bi.

Fig. 2

 2 Fig.2Physical states T Γ , φ l,s are used to determine the configuration of the front Γ : solidification, fusion or conduction when T fus = T sol . K sol , K fus , T fus , T sol are constants dependent on the material and geometrical situations. Notice that the change from solidification to fusion must transition by the conduction situation.

Fig. 3

 3 Fig.3Left: sketch of the stratified liquid, solid crust and the surrounding vessel in a 2D x, z geometry. Γ is the phase front -Right: Sketch of the proposed 2D model where the 2D heat equation is solved in the solid phase. The solid domain Ωs(t) evolves in time due to solidification or fusion at the liquid/solid phase front Γ while the liquid state evolves in composition and location with a stratified 0D multi-layer model. The solid is meshed with an unstructured mesh and the heat equation is solved subject to flux/temperature boundary conditions from the liquid and solid neighbour domains (and top/bottom external world)

Fig. 5

 5 Fig.5Sketch of the finite element RT 0 -P 0 -The temperature is located on the element center (•), while the heat flux is located on the face center and normal to the face. (×).

Fig. 6

 6 Fig. 6 Left: Discretization of the phase front Γ= I i+1/2 = [ x i , x i+1 ], s i = (x i , z i ), 1 ≤ i ≤ N i + 1 .δx and δz are the smallest allowed measures for a segment/face/cell -Right: Illustration of the assumption of instantaneous quadratic temperature profile T j within the crust layer C j knowing boundary and mean temperatures.

Fig. 7

 7 Fig.7Stefan problem -Top: Sketch of the configuration with the exact initial position s(t 0 ) (left) and s(t) at time t > t 0 (right) -Left/Right: Position of the front s(t) as a function of time t for the current simulations (coloured symbols) for different time-step (left) and step sizes (right) vs the exact solution (black straight line), and zoom to observe the differences.

Fig. 8

 8 Fig.8Sketch of the solidification generated by an infinite corner shape. A liquid material at temperature T liq is filling Ω = R 2 /B at t = 0 where B = (x, y) ∈ R 2 , x ≤ 1/20, or y ≤ 1/20 at t = 0. The liquid solidifies in presence of a cold solid in B initially at T sol cooled down from the semi-infinite axes which are maintained below the solidification temperature, i.e T ext < T sol . As time increases, 0 < t 1 < t 2 < t 3 < t 4 , the phase front adopts a 2D curved shape, while, away from the origin, the 1D Stefan theory applies.

Fig. 9 Fig. 10

 910 Fig. 9 Solidification generated by an infinite corner shape -Left: numerical results of the phase front shape for different intermediate times -Right: temperature profile associated to the time t = 2.4 × 10 -3 .

Fig. 11

 11 Fig. 11 Vessel fusion test case -Conduction: Bi = 0.8 -Temperature profile (colours) for the 2D (left column) and 0D models (right column) -Intermediate mesh size ∆x = 0.01 -Left panel: short time t = 0.03. Middle panel: time t = 1. Right panel: final time t = 4.

  -top the output heat flux at the five gauges as a function of time for the 2D model. Coherent physically admissible heat-flux profiles are observed. Next on bottom panel of figure 14 we compare the 0D and 2D models on gauges A, C and D. For the 0D model, the output heat flux profiles are the same at gauges A and C because this model neglects the multi-dimensional conduction. For the same reason, beyond z * , where the input heat flux is 0, the facing output heat flux is constant equal to 0 (in D and E). Another artefact is again the spurious positive value of the heat flux observed at short time, before time 0.15. Also the change in profile continuity marks the beginning of the fusion observed around time t 0.5. The 2D model captures the 2D conduction and its effects on the heat flux along the z-direction still. It maintains the physically admissible flux value between 0 and -1.It is important to notice than even on gauge A, the two models already differ by 30% on the heat flux value. In figure15we display the temperature and fusion front results at final time obtained for the coarse, intermediate and fine meshes. They are better captured with finer meshes and no spurious phenomenon is observed when the mesh gets refined. A visual investigation convinced that the temperature beyond z * reaches values of the order 0.3 -0.6 for z * ≤ z ≤ 2.4, and, 0.1 -0.3 for 2.4 ≤ z ≤ 2.8.

Fig. 12

 12 Fig. 12 Vessel fusion test case -Conduction: Bi = 0.8 -Intermediate mesh size ∆x = 0.01 -Top: Heat fluxes ϕout at gauges A, B, C, D, E (colours and symbols) as a function of time -Bottom: Heat fluxes comparison at gauges A, C and D between 0D (red curve) and 2D models (blue symbols).

Fig. 13

 13 Fig. 13 Vessel fusion test case -Fusion: Bi = 1.2 -Temperature profile (colours) and the fusion front for the 2D model (left column) and 0D model (right column) -Intermediate mesh size ∆x = 0.01 -From left to right: short time t = 0.03, t = 1 and final time t = 4.

Fig. 14

 14 Fig. 14 Vessel fusion test case -Fusion: Bi = 1.2 -Intermediate mesh size ∆x = 0.01 -Top: Heat fluxes at gauges A, B, C, D, E in colours as a function of time -Bottom: Heat fluxes comparison at gauges A, C and D between 0D (red curve) and 2D models (blue symbols).

Fig. 15

 15 Fig. 15 Vessel fusion test case -Fusion: Bi = 1.2 -Temperature profile (colours) and the fusion front for the 2D model -Final time t final = 4 -Left to right panels: coarse, intermediate and fine mesh sizes.

Fig. 16

 16 Fig. 16 Vessel fusion test case -Fusion: Bi = 1.2 -Intermediate mesh size ∆x = 0.01 -Mass ratio as a function of time for the different x-layers at gauges A, B and C in straight red line for the 0D model, in blue symbols for the current 2D one.

Fig. 17

 17 Fig. 17Full 2D test case -Fusion: Bi = 1.2 at t = 1 and t = 2 -Mesh size ∆x = 0.05 -2D conduction and 2D ablation on the left, 2D conduction and 1D horizontal ablation on the right -Temperature in colour.

  Fig. 17Full 2D test case -Fusion: Bi = 1.2 at t = 1 and t = 2 -Mesh size ∆x = 0.05 -2D conduction and 2D ablation on the left, 2D conduction and 1D horizontal ablation on the right -Temperature in colour.

Fig. 18

 18 Fig. 18 Industrial-like test case -Representation of the liquid corium on the left and zoom on the crust on the right for the different stages (1) (t = 5950 s), (2) (t = 8400 s), (3) (t = 20500 s) for the 2D model -In colour: the mass fraction of main oxide species : ZrO 2 , U O 2 , F eO.

Fig. 19

 19 Fig.19Industrial-like test -Comparison 0D and 2D modelling approaches -Left: mass of layers in the liquid corium and crust as a function of time. Black: oxide liquid layer. Red: light metal liquid layer. heavy metal liquid layer. Purple: solidified crust. Curves with crosses (×) correspond to the results of 2D model and those with circles (•) to the 0D model -Right: differences of the total solidified crust mass ∆m between 0D and 2D models as a function of time.

  the following we consider a dimensionless version of the equations with the following scaling quantities for time, space, temperature and heat flux: t * , x * , z * , T * , Φ * . The associated dimensionless variables are then τ = the temperature, and, Φ = ϕ Φ * , Φ in = ϕ in Φ * for the heat flux. From those variables we can deduce the following dimensionless numbers: Fourier number Fo = λt * ρCpx 2 * , Biot number Bi = Φ * x * λT * , and, the aspect ratio of the problem Γ = x * z * . As a consequence, the heat equation with dimensionless numbers BiΦ in for χ = 0, ξ ∈ [0, ξ 0 ], 0 if χ = 0, ξ ∈ [ξ 0 , ξmax] or ξ = 0 or ξ = ξmax(34)θ = θ ext for χ = χmax,(35)We arbitrarily choose Fo = 1 which implies t * = x 2 * ρCp λ . Moreover we take x * = xmax, z * = z 0 = 2, Φ * = ϕ in , T * = T fus . At last, for sake of simplicity we choose T ext = T 0 = 0. For this test case we maintain a constant aspect ratio Γ = x

Table 2

 2 Stefan problem -Errors in L 1 , L 2 and L∞ norms at final time for a refinement in time with ∆x = 0.1 (top) and in space with ∆t = 0.3 (bottom).

= T = T sol T T

Acknowledgements

This work has been carried out within the framework of the PROCOR platform development funded by CEA, EDF and Framatome. RL acknowledges the financial and scientific support of CEA Cadarache.