Hamza Boukili
email: boukili.hamza@gmail.com

Jean-Marc Hérard
email: jean-marc.herard@edf.fr

Simulation and preliminary validation of a three-phase flow model with energy

Keywords: Three-phase flows, vapour explosion, hyperbolic systems, finite volumes, relaxation time scales

 which focuses on vapour explosion simulation.

Introduction

In fluid dynamics, a multiphase flow is a simultaneous flow of two or more components in the same media, or the flow of one single component but with different chemical properties or thermodynamic states. It has a wide range of applications in the modern industry that cover different flow configurations, for instance: gas-liquid-solid mixtures, non miscible liquidliquid mixtures, change of state liquid-gas mixtures, etc.

In the nuclear industry, many configurations of multiphase flows arise, either in the normal operation (for instance in the steam generator : liquid water -water vapor flow) or the accidental operation (for instance: Loss of Coolant Accident [START_REF] Glossary | Loss of coolant accident (LOCA)[END_REF] or Reactivity Initiated Accident [START_REF]Nuclear Energy Agency Committee on the Safety of Nuclear Installations, Reactivity Initiated Accident (RIA) Fuel Codes Benchmark Phase-II[END_REF]). Vapor Explosion [START_REF] Berthoud | Vapor explosions[END_REF] arises as one of the multiphase flows that could take place in an accidental scenario. It could happen in the case of a meltdown of the reactor core, if the corium (i.e. mixture of molten nuclear fuel, fission products, control rods, reactor vessel structure materials, etc.) would get in touch with liquid water. The contact between these two fluids is generally characterized by an intense and rapid heat transfer that leads to a strong evaporation of the liquid water. During the expansion of water vapor, pressure shock waves could be formed, which might damage the surrounding structures. Thus, the modeling and the numerical simulation of Vapor Explosion is a completely relevant and highly challenging task to undertake.

Many efforts have been devoted to the modeling and simulation of Vapor Explosion. This requires at least a three-phase flow model to account for the corium phase, the liquid water phase, and the water vapor phase; it should be noted that the latter three are not miscible. Some models consider more than three fields, this could be useful for instance to take into account the different gaseous phases or fission products (see for instance [START_REF] Meignen | The challenge of modelling fuel-coolant interaction. Part II: steam explosion[END_REF]). However, existing multiphase models generally suffer from typical mathematical short-comings, one of which is the loss of hyperbolicity, and also the non-uniqueness of jump conditions.

For the sake of a correct representation of shock and rarefaction waves, compressible multiphase flow models with hyperbolicity criterion and unique jump conditions are required. The hyperbolicity provides an assurance of the well-posedness of the initial-value problem, while the uniqueness of jump conditions is mandatory, for non conservative first-order systems, otherwise computed shock patterns totally depend on the chosen scheme, which of course does not make sense (see [START_REF] Berthon | Contribution à l'analyse numérique des equations de Navier-Stokes compressibles à deux entropies spécifiques[END_REF][START_REF] Hou | Why nonconservative schemes converge to wrong solutions: error analysis[END_REF][START_REF] Guillemaud | Modélisation et simulation numérique des écoulements diphasiques par une approche bifluide à deux pressions[END_REF]).

In the two-phase flow framework, Baer and Nunziato [START_REF] Baer | A two phase mixture theory for the deflagration to detonation transition (DDT) in reactive granular materials[END_REF] proposed a rather appealing model that allows to describe the strong shock and rarefaction waves, and also the compressible effects, while preserving the hyperbolicity condition. This model has been studied in [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF] and extended by A.K. Kapila et al. [START_REF] Kapila | Two phase modeling of deflagration to detonation transition : srtucture of the velocity relaxation zone[END_REF] and S. Gavrilyuk and R. Saurel [START_REF] Gavrilyuk | Mathematical and numerical modelling of two-phase compressible flows with micro inertia[END_REF] among others. Another category of models exists in the two-phase flow literature [START_REF] Drew | Theory of Multicomponent Fluids[END_REF][START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flows[END_REF], it consists of considering one pressure field for the two fluids, or in other words with instantaneous pressure equilibrium between all phases, but this approach does not ensure the preservation of the mathematical properties we mentioned above, namely the hyperbolicity and the existence and uniqueness of jump conditions.

In [START_REF] Hérard | A three-phase flow model[END_REF] an extension of Baer-Nunziato model to the three-phase flow framework for non miscible components was suggested. The main difficulty that arises is the description of what happens at the interface between the three phases. In fact, in the two-phase mixture there is only one interface velocity and one interface pressure that need to be defined, whereas in the three-phase mixture more closure laws need to be introduced at the statistical interface. The paper [START_REF] Hérard | A three-phase flow model[END_REF] addresses the main modeling choices that allow to guarantee a relevant entropy inequality, and also the uniqueness of jump conditions. Other recent works have allowed to get more results considering the same class of two or three-phase flow models [START_REF] Müller | Closure conditions for non-equilibrium multi-component models[END_REF][START_REF] Hérard | Some mathematical properties of a hyperbolic multiphase flow model[END_REF][START_REF] Hantke | Closure conditions for a one temperature non-equilibrium multicomponent model of Baer-Nunziato type[END_REF][START_REF] Hérard | A three-phase flow model with two miscible phases[END_REF].

In the present work, we consider the same class of three-phase flow models as in [START_REF] Hérard | A three-phase flow model[END_REF], and we especially take a deeper look at the source terms that account for the relaxation processes between the three phases. The present paper is organized as follows. In Section 2, we present the set of equations of the three phase flow model, as well as the different closure laws and the mathematical properties of the closed model. In this section, we also discuss the properties of the relaxation processes: velocity, pressure, temperature and mass transfer. These four distinct relaxation effects indeed have a significant impact on the dynamics of the three-phase mixture. In Section 3, we focus on the discrete framework, the global resolution strategy actually consists of the fractional step method [START_REF] Hérard | A fractional step method to compute a class of compressible gas-liquid flows[END_REF]. We give the details of the numerical scheme that was considered to take the convective subset into account, since the latter includes non conservative terms. We also present the numerical schemes for the relaxation effects. In Section 4, we present the different numerical results that consist of:

• Verification of the convective subset: pure convection test cases, which involve a Riemann problem;

• Verification of each relaxation subset;

• Validation of the dynamic effects on a shock tube apparatus [START_REF] Chauvin | Étude expérimentale de l'atténuation d'une onde de choc par un nuage de gouttes et validation numérique[END_REF][START_REF] Chauvin | Experimental investigation of the propagation of a planar shock wave through a two-phase gas-liquid medium[END_REF];

• Preliminary validation of the whole model (with dynamic and thermal effects) on a vapour explosion test case, based on the KROTOS experiment [START_REF] Huhtiniemi | KROTOS 38 to KROTOS 44[END_REF].

Seven appendices devoted to technical details complete the whole.

2. Governing equations and main properties of the three-dimensional threephase flow model

Governing equations

We consider the following system of partial differential equations for the modeling of a three-phase flow with two non miscible components (water, liquid metal), assuming that water may be present under two distinct states: liquid or vapor phase. For k = 1, .., 3 and t > 0:

       ∂α k ∂t + V i (W).∇α k = S α k (W) ∂m k ∂t + ∇.(m k U k ) = S m k (W) ∂m k U k ∂t + ∇.(m k U k ⊗ U k + α k p k Id) + 3 l=1,l =k Π kl (W)∇α l = S U k (W) ∂α k E k ∂t + ∇.(α k E k U k + α k p k U k ) -3 l=1,l =k Π kl (W) ∂α l ∂t = S E k (W) (1)
The quantities

α k ∈ [0, 1], ρ k , m k = α k ρ k , U k , p k , e k (p k , ρ k
) and E k represent respectively the statistical fraction, the mean density, the partial mass, the mean velocity, the mean pressure, the mean internal energy and the mean total energy of phase k, k = 1, 2, 3, where:

E k = 1 2 ρ k U k .U k + ρ k e k (p k , ρ k ) (2)
Since the liquid water, water vapor and liquid metal are not miscible, we have the following constraint on statistical fractions:

α 1 + α 2 + α 3 = 1 (3)
The state variable W ∈ R 17 denotes the following vector:

W = (α 2 , α 3 , m 1 , m 2 , m 3 , m 1 U 1 , m 2 U 2 , m 3 U 3 , α 1 E 1 , α 2 E 2 , α 3 E 3 ) t (4)
Some additional thermodynamic variables need to be defined. We set:

c 2 k (p k , ρ k ) = p k ρ 2 k - ∂e k (p k , ρ k ) ∂ρ k ∂e k (p k , ρ k ) ∂p k -1 (5)
s k (p k , ρ k ), the specific entropy of phase k, is defined such that:

c 2 k ∂s k (p k , ρ k ) ∂p k + ∂s k (p k , ρ k ) ∂ρ k = 0 ( 6
)
The temperature is given by:

1 T k = ∂s k ∂p k ∂e k ∂p k -1 (7)
And µ k denotes the Gibbs potential:

µ k = e k + p k ρ k -T k s k (8)
The first and second equations of [START_REF] Ambroso | Relaxation and numerical approximation of a two-fluid two-pressure diphasic model[END_REF] give the evolution of statistical fractions and partial masses, while the third and fourth equation stand for the momentum balance and energy balance equations.

In this work, the interface velocity is defined by:

V i (W) = U 1 (9)
Thus, following [START_REF] Hérard | A three-phase flow model[END_REF], the interface pressures are given by:

Π 12 (W) = Π 21 (W) = Π 23 (W) = p 2 Π 13 (W) = Π 31 (W) = Π 32 (W) = p 3 (10)
Remark 1:

The choice made in (9) is not unique. In fact, an other possibility consists of considering V i (W) as a convex combination of velocities U k , k = 1, 2, 3:

V i (W) = 3 k=1 m k U k 3 k=1 m k
This leads to a different set of interface pressures Π kl (W) which is uniquely defined, and in agreement with the entropy inequality (see paragraph 2.2.1 and Appendix G of [START_REF] Hérard | A three-phase flow model[END_REF]).

Source terms should be such that:

3 k=1 S α k (W) = 3 k=1 S m k (W) = 3 k=1 S E k (W) = 0 ; 3 k=1 S U k (W) = 0
since they only take into account the internal transfers between phases.

Closure laws to account for the different relaxation effects involved in (1) are given by:

• Pressure relaxation:

S α k (W) = 3 l=1,l =k K P kl (W )(p k -p l ) (11)
where:

K P kl (W) = α k α l τ P kl (W)Π 0 , τ P kl (W) = τ P lk (W)
is symmetric positive, and represents the pressure-relaxation time scale between phases k and l (see [START_REF] Gavrilyuk | The structure of pressure relaxation terms: the one-velocity case[END_REF]), and Π 0 is a positive reference pressure.

• Mass transfer:

S m k (W) = 3 l=1,l =k Γ kl (W) = 3 l=1,l =k K m kl (W) µ l T l - µ k T k (12)
where:

K m kl (W) = 1 τ m kl (W)Γ 0 m k m l m k +m l .
The symmetric positive function τ m kl (W) represents the characteristic mass transfer time scale between phases k and l, and Γ 0 is a positive constant (with dimension µ/T ).

• Momentum interfacial transfer term:

S U k (W) = 3 l=1,l =k D kl (W) + 3 l=1,l =k Ũkl Γ kl (W) (13)
where D kl represents the drag effect between phases k and l:

D kl (W) = e kl (W)(U l -U k ) (14)
The terms e kl are chosen under the form:

e kl (W) = m k m l τ U kl (W)M 0 (15)
where τ U kl (W) is symmetric positive, and accounts for the velocity-relaxation time scale. M 0 = m 1 + m 2 + m 3 is the total mass. The velocity Ũkl is chosen under the symmetric form:

Ũkl (W) = U k + U l 2 (16)
• Energy balance source term:

S E k (W) = 3 l=1,l =k V kl (W).D kl (W) + 3 l=1,l =k ψ kl (W) + 3 l=1,l =k Hkl Γ kl (W) (17)
where:

V kl = 1 2 (U k + U l ) ( 18
)
and Hkl is given by: Hkl

(W) = U k .U l 2 (19)
The term ψ kl accounts for the heat transfer between phases k and l, it is defined by:

ψ kl (W) = K T kl (W)(T l -T k ) (20)
with:

K T kl (W) = 1 τ T kl (W) m k m l C V,k C V,l m k C V,k + m l C V,l (21)
Here again τ T kl (W) is symmetric positive and accounts for the heat transfer characteristic time between phases k and l, and C V,k denotes the volumetric heat capacity of phase k.

Main properties of the three-phase flow model 2.2.1. Entropy

We define η(W), the mixture entropy, and F η (W), the mixture entropy flux, such that:

η(W) = -3 k=1 m k Log(s k ) F η (W) = -3 k=1 m k Log(s k )U k (22)
We have the following property:

Property 1:
Considering closure laws (9-21), the following entropy inequality holds for smooth solutions of system (1):

∂η(W) ∂t + ∇.F η (W) ≤ 0 (23)
Proof:

For regular solutions of the system (1), the governing equation of η(W) reads:

∂η(W) ∂t + ∇.F η (W) = - 3 k=1 1 T k p k S α k (W) + l =k Π kl (W)S α l (W) - 3 k=1 1 T k 3 l=1,l =k ((V kl (W) -U k ).D kl (W) + ψ kl (W)) + 3 k=1 µ k T k 3 l=1,l =k Γ kl (W) (24)
Now, on the basis of the interface pressure definitions [START_REF] Coquel | A positive and entropy-satisfying finite volume scheme for the Baer Nunziato model[END_REF] and the closure laws presented in [START_REF] Crouzet | Validation of a two-fluid model on unsteady water-vapour flows[END_REF][START_REF] Drew | Theory of Multicomponent Fluids[END_REF][START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF][START_REF] Gavrilyuk | The structure of pressure relaxation terms: the one-velocity case[END_REF][START_REF] Gavrilyuk | Mathematical and numerical modelling of two-phase compressible flows with micro inertia[END_REF][START_REF] Gelfand | Droplet breakup phenomena in flows with velocity lag[END_REF][START_REF] Guillemaud | Modélisation et simulation numérique des écoulements diphasiques par une approche bifluide à deux pressions[END_REF][START_REF] Han | Efficient and robust relaxation procedures for multi-component mixtures including phase transition[END_REF][START_REF] Hantke | Closure conditions for a one temperature non-equilibrium multicomponent model of Baer-Nunziato type[END_REF][START_REF] Hérard | A three-phase flow model[END_REF][START_REF] Hérard | A fractional step method to compute a class of compressible gas-liquid flows[END_REF] a straightforward calculus provides:

3 k=1 1 T k p k S α k (W) + l =k Π kl (W)S α l (W) = 1 T 1 3 k=1 p k S α k (W) = 1 T 1 3 k,l=1,k =l α k α l τ P kl Π 0 (p k -p l ) 2 (25)
3 k=1 1 T k l =k (V kl (W) -U k ).D kl (W) = 1 2 3 k=1 1 T k l =k e kl (W)(U l -U k ) 2 (26) 3 k=1 1 T k l =k ψ kl (W) = 1≤l<k≤3 K T kl (W) T k T l (T k -T l ) 2 (27) 3 k=1 µ k T k 3 l=1,l =k Γ kl (W) = 3 k=1 µ k T k 3 l=1,l =k K m kl (W) µ l T l - µ k T k = - 1≤l<k≤3 K m kl (W) µ l T l - µ k T k 2 (28)
It obviously follows that the regular solutions of system (1) comply with the inequality [START_REF] Hérard | Some mathematical properties of a hyperbolic multiphase flow model[END_REF].

We also note that the right-hand side of (24) vanishes as soon as pressure, velocity, temperature and Gibbs potential equilibria are reached.

Hyperbolicity and structure of fields

We focus here on the convective subset -left hand side-of (1). Let n be a unit vector, and τ 1 ,τ 2 such that (n, τ 1 , τ 2 ) defines an orthonormal basis of the 3D space. Considering the invariance of equations (1) under frame rotation, and neglecting transverse derivatives of all components, the associated one-dimensional problem in the n direction writes:

               ∂ t α k + (U 1 .n)∂ xn α k = 0 ∂ t m k + ∂ xn (m k (U k .n)) = 0 ∂ t (m k (U k .n)) + ∂ xn (m k (U k .n) 2 + α k p k ) + 3 l=1,l =k Π kl (W)∂ xn α l = 0 ∂ t (α k E k ) + ∂ xn (α k (E k + p k )(U k .n)) -3 l=1,l =k Π kl (W)∂ t α l = 0 ∂ t (m k (U k .τ 1 )) + ∂ xn (m k (U k .n)(U k .τ 1 )) = 0 ∂ t (m k (U k .τ 2 )) + ∂ xn (m k (U k .n)(U k .τ 2 )) = 0 (29)
Property 2: 2.1 The system (29) admits the following real eigenvalues:

λ 1,2,3,4,5 (W) = U 1 .n ; λ 6,7,8 (W) = U 2 .n ; λ 9,10,11 (W) = U 3 .n λ 12,13 (W) = U 1 .n ± c 1 ; λ 14,15 (W) = U 2 .n ± c 2 ; λ 16,17 (W) = U 3 .n ± c 3 ( 30
)
Associated right eigenvectors span the whole space R 17 if:

(U 1 -U k ).n = ±c k (31)
2.2 Fields λ k for k ∈ 1, .., 11 are Linearly Degenerated ; other fields are Genuinely Non Linear.

Proof:

For smooth solutions the system (29) can be written as:

∂ t Z + A(Z)∂ xn Z = 0 (32)
The computation of eigenvalues of the matrix A is rather easy when choosing the variable:

Z = (α 2 , α 3 , s 1 , s 2 , s 3 ,U 1 .n, U 1 .τ 1 , U 1 .τ 2 , U 2 .n, U 2 .τ 1 , U 2 .τ 2 , U 3 .n, U 3 .τ 1 , U 3 .τ 2 , p 1 , p 2 , p 3 ) t (33)
(See more detailed calculations in the one-dimensional framework in [START_REF] Hérard | A three-phase flow model[END_REF] Appendix A).

We also check that: ∂λ k (Z) ∂Z .r k (Z) = 0 for k ∈ 1, .., 11, and ∂λ k (Z) ∂Z .r k (Z) = 0 otherwise.

Additional properties in the one-dimensional framework

We consider a pure one-dimensional problem, thus system (29) can be written as:

       ∂ t α k + u 1 ∂ x α k = 0 ∂ t m k + ∂ xn (m k u k ) = 0 ∂ t (m k u k ) + ∂ x (m k u 2 k + α k p k ) + 3 l=1,l =k Π kl (w)∂ x α l = 0 ∂ t (α k E k ) + ∂ xn (α k (E k + p k )u k ) -3 l=1,l =k Π kl (w)∂ t α l = 0 ( 34
)
where the 1D state variable is denoted w ∈ R 11 :

w = (α 2 , α 3 , m 1 , m 2 , m 3 , m 1 u 1 , m 2 u 2 , m 3 u 3 , α 1 E 1 , α 2 E 2 , α 3 E 3 ) t ( 35
)
Property 3:

• The convective system (34) admits eleven real eigenvalues which read:

λ 1,2,3 (w) = u 1 ; λ 4 (w) = u 2 ; λ 5 (w) = u 3 λ 6,7 (w) = u 1 ± c 1 ; λ 8,9 (w) = u 2 ± c 2 ; λ 10,11 (w) = u 3 ± c 3 ( 36
)
The field associated with λ k (w) for k ∈ {1, .., 5} is Linearly Degenerated, while fields associated with λ 6-11 (w) are Genuinely Non Linear.

• Regarding the 1D Riemann problem associated with [START_REF]Nuclear Energy Agency Committee on the Safety of Nuclear Installations, Reactivity Initiated Accident (RIA) Fuel Codes Benchmark Phase-II[END_REF], the LD field λ 1,2,3 (w) admits the following eight Riemann invariants:

I 1 1,2,3 (w) = m 2 (u 2 -u 1 ) ; I 2 1,2,3 (w) = m 3 (u 3 -u 1 ) I 3 1,2,3 (w) = s 2 ; I 4 1,2,3 (w) = s 3 ; I 5 1,2,3 (w) = u 1 I 6 1,2,3 (w) = 3 k=1 α k p k + m 2 (u 2 -u 1 ) 2 + m 3 (u 3 -u 1 ) 2 I 7 1,2,3 (w) = e 2 + p 2 ρ 2 + 1 2 (u 2 -u 1 ) 2 ; I 8 1,2,3 (w) = e 3 + p 3 ρ 3 + 1 2 (u 3 -u 1 ) 2 (37)
• We note ∆(φ) = φ r -φ l . For each isolated GNL wave, the following exact jump conditions hold for phase index k = 1, 2, 3, through any discontinuity separating states l, r and moving with speed σ:

       ∆(α k ) = 0 ∆(ρ k (u k -σ)) = 0 ∆(ρ k u k (u k -σ) + p k ) = 0 ∆(E k (u k -σ) + p k u k ) = 0 (38)
Proof:

The proof of these properties is classical and left to the reader.

The latter Riemann invariants and jump conditions are particularly important because they will enable us to build exact solutions of the 1D Riemann problem for the system [START_REF]Nuclear Energy Agency Committee on the Safety of Nuclear Installations, Reactivity Initiated Accident (RIA) Fuel Codes Benchmark Phase-II[END_REF], and verify the convergence of algorithms (see paragraph 4.1). We mention that these properties are the exact counterpart of the two-phase Baer-Nunziato model [START_REF] Baer | A two phase mixture theory for the deflagration to detonation transition (DDT) in reactive granular materials[END_REF].

Admissibility of thermodynamic quantities

In this subsection, we study the admissibility of thermodynamic variables in the onedimensional framework. We consider smooth solutions of the system (34), and we make the following assumptions: Assumption 1: For k = 1, 2, 3 we assume admissible initial and boundary conditions, i.e.: ∀x ∈ Ω :

   0 ≤ α k (x, 0) 0 ≤ m k (x, 0) 0 ≤ s k (x, 0)
and ∀t ∈ [0, T ] :

   0 ≤ α k (x Γ , t) 0 ≤ m k (x Γ , t) 0 ≤ s k (x Γ , t) (39)
where Γ denotes the boundary of the domain Ω, and T is the simulation final time.

Assumption 2:

We also assume that functions a k (defined in (??)), u k and ∂u k ∂x remain in L ∞ (Ω × [0, T ]).

We have then the following result:

Property 4:
The regular solutions of the system (34) are consistent with the physical requirements:

∀(x, t) ∈ (Ω × [0, T ]) :    0 ≤ α k (x, t) 0 ≤ m k (x, t) 0 ≤ s k (x, t) (40)
Proof:

It is classical, and based on building the evolution equations of the α k , m k and s k as smooth solutions of the system (34) (See [START_REF] Hérard | A three-phase flow model[END_REF] for some more details).

We can prove that a similar admissibility result holds for pressure variables p k , when we consider the specific case of Stiffened Gas equations of state (EOS). In that case, for each phase k = 1, 2, 3:

p k + γ k Π k = (γ k -1)ρ k e k (41)
where γ k > 1 and Π k > 0 are thermodynamic constants of phase k. We denote also:

P k = α k (p k + Π k )/(γ k -1) (42)
For regular solutions of the system (34), the evolution equations of P k read:

   ∂ t P 1 + ∂ x (u 1 P 1 ) + (γ 1 -1)P 1 ∂ x u 1 = 0 ∂ t P 2 + ∂ x (u 2 P 2 ) + (γ 2 -1)P 2 (∂ x u 2 + (u 2 -u 1 )∂ x Log(α 2 )) = 0 ∂ t P 3 + ∂ x (u 3 P 3 ) + (γ 3 -1)P 3 (∂ x u 3 + (u 3 -u 1 )∂ x Log(α 3 )) = 0 (43)
Based on the same type of arguments as those invoked in the proof of Property 4, we conclude that ∀(x, t) ∈ (Ω × [0, T ]) : 0 ≤ P k (x, t), provided that, in addition, (u ku 1 )∂ x Log(α k )) remains bounded.

Of course, the latter results can easily be extended to the 3D framework.

Velocity relaxation

In this subsection, we focus on the relaxation mechanisms embedded in our PDE system. More particularly, we study what happens in a continuous framework, considering an homogeneous flow such that the spatial derivatives are null. Thus, only time derivatives and source terms are taken into account. The studied system is:

       ∂ t α k = S α k (W) ∂ t m k = S m k (W) ∂ t (m k U k ) = S U k (W) ∂ t (α k E k ) -3 l=1,l =k Π kl (W) ∂α l ∂t = S E k (W) (44)
The system (44) is itself split into four subsystems, in order to study separately the different relaxation effects. In this paragraph, we study the velocity relaxation effects, the corresponding PDE system is obtained from (44) by considering only the velocity-related source terms. The velocity relaxation system then writes:

                       ∂ t α k = 0 ∂ t m k = 0 ∂ t (m k U k ) = 3 l=1,l =k D kl (W) ∂ t (α k E k ) = 3 l=1,l =k V kl (W).D kl (W) (45a) (45b) (45c) (45d)
We have the following noteworthy property:

Property 5:
The velocity relaxation step is consistent with the admissibility of the internal energy, i.e. this step keeps the internal energy e k in the admissible range.

Proof:

Starting with:

α k E k = m k e k + 1
2 m k U 2 k , using the closures ( 14)- [START_REF] Han | Efficient and robust relaxation procedures for multi-component mixtures including phase transition[END_REF] and combining (45c) -(45d) we obtain:

m k ∂ t e k = 1 2 3 l=1,l =k e kl (W)(U l -U k ) 2 ≥ 0 (46)
This means that the internal energy is non-decreasing, regardless of the chosen EOS.

Thus, the integration method we adopt consists, first, of obtaining the velocity variation over time, and then using the result to update the energy. The velocity variation over time is obtained by studying the velocity differences, as shown in what follows.

Let (n, τ 1 , τ 2 ) be an orthonormal basis of the 3D space. If we set X n 12 = U 1 .n -U 2 .n, X n 13 = U 1 .n-U 3 .n, and X n = (X n 12 , X n 13 ) t , the momentum equation of (45) could be written as:

∂ t X n = -A(W)X n (47)
where:

A(W) = a 11 a 12 a 21 a 22        a 11 = e 12 ( 1 m 1 + 1 m 2 ) + e 23
m 1 + 1 m 3 ) + e 23 m 3 (48a) (48b)
We obtain formally then:

X n (t) = X n (t 0 )exp - t t 0 A(s)ds (49)
A similar result stands for X τ 1 and X τ 2 , when considering the normal velocities in the directions τ 1 and τ 2 , with the same matrix A(W). This comes from the structure of (45c), which remains unchanged regardless of the projection direction. This velocity variation mechanism (49) is identical to what happens in the barotropic framework, which has been studied in a previous work [START_REF] Boukili | Relaxation and simulation of a barotropic three-phase flow model[END_REF]. Equation (49) gives then the evolution of velocity differences, the total kinetic energy conservation allows, formally also, to retrieve the variation over time of velocities U k . Details are given in the numerical scheme subsection 3.3.1.

Pressure relaxation

In this paragraph, we focus on the pressure relaxation effects. The concerned PDE system, extracted from (44), writes:

                 ∂ t α k = S α k (W) ∂ t m k = 0 ∂ t (m k U k ) = 0 ∂ t (α k E k ) - 3 l=1,l =k Π kl (W)∂ t α l = 0 (50a) (50b) (50c) (50d)
In order to understand the underlying relaxation mechanism, we study, in this sub-step, the evolution of the pressure differences. To do so, we define the following notations:

• The coefficients: A k = ρ k c 2 k α k and b 1 = 1 ρ 1 ∂e 1 ∂p 1 -1 ;
• The pressure differences: y 12 = p 1 -p 2 , y 13 = p 1 -p 3 , y 23 = -y 32 = y 13 -y 12 and Y = (y 12 , y 13 ) t .

We have the following result:

Property 6:
The evolution of the pressure differences writes:

∂ t Y = - 1 Π 0 B(W)Y (51)
where the matrix

B(W) = b 11 b 12 b 21 b 22
is given by:

               b 11 = A 1 + A 2 -b 1 α 1 y 12 α 1 α 2 τ 12 + A 2 -b 1 α 1 y 32 α 2 α 3 τ 23 b 12 = A 1 -b 1 α 1 y 13 α 1 α 3 τ 13 -A 2 -b 1 α 1 y 32 α 2 α 3 τ 23 b 21 = A 1 -b 1 α 1 y 12 α 1 α 2 τ 12 -A 3 -b 1 α 1 y 23 α 2 α 3 τ 23 b 22 = A 1 + A 3 -b 1 α 1 y 13 α 1 α 3 τ 13 + A 3 -b 1 α 1 y 23 α 2 α 3 τ 23 (52)
Proof:

On the one hand, equation (50b) yields:

∂ t ρ k = - ρ k α k ∂ t α k = - ρ k α k S α k (W)
On the other hand, combining (50c) and (50d) gives:

m k ∂ t e k = 3 l=1,l =k Π kl (W)∂ t α l = 3 l=1,l =k Π kl (W)S α l (W) Knowing that ∂ t e k = ∂e k ∂p k ∂ t p k + ∂e k ∂ρ k ∂ t ρ k , we can write: ∂e k ∂p k ∂ t p k + ∂e k ∂ρ k ∂ t ρ k = 1 m k 3 l=1,l =k Π kl (W)S α l (W)
By using the evolution equation of ρ k above we get:

∂ t p k = ∂e k ∂p k -1 1 m k 3 l=1,l =k Π kl (W)S α l (W) + ρ k α k ∂e k ∂ρ k S α k (W)
By expliciting the interface pressures Π kl (W), given in [START_REF] Coquel | A positive and entropy-satisfying finite volume scheme for the Baer Nunziato model[END_REF], we obtain:

   ∂ t p 1 = -A 1 S α 1 (W) + b 1 α 1 3 k=1 p k S α k (W) ∂ t p 2 = -A 2 S α 2 (W) ∂ t p 3 = -A 3 S α 3 (W)
Finally, a straightforward, though cumbersome, calculation gives (51) and (52), using the closure of S α k (W) given in [START_REF] Crouzet | Validation of a two-fluid model on unsteady water-vapour flows[END_REF].

Remark 2:

The matrix B(W) defined in (51)-(52) depends in particular of the pressure differences themselves. Actually, as shown in (52), each coefficient of B(W) contains a first degree polynomial wrt Y. Thus, a threshold effect may arise, as detailed in Appendix 2 while focusing on two-phase flows. We emphasize that this effect does not exist in the barotropic framework.

Turning to the three-phase flow model (1), a necessary condition to be verified for initial conditions is that the trace of B(W) should be strictly positive (whatever B(W) admits real or complex eigenvalues). Otherwise the return to the pressure equilibrium cannot be guaranteed.

Since the trance of B reads:

tr(B(W)) = b 11 + b 22 = 1 τ P 12 Π 0 α 2 ρ 1 c 2 1 + α 1 ρ 2 c 2 2 -α 2 b 1 y 12 + 1 τ P 13 Π 0 α 3 ρ 1 c 2 1 + α 1 ρ 3 c 2 3 -α 3 b 1 y 13 + 1 τ P 23 Π 0 α 2 ρ 3 c 2 3 + α 3 ρ 2 c 2 2
We may set:

|y| = max(|y 12 |, |y 13 |)
Hence, a sufficient condition on |y| to ensure a positive trace will be:

|y| < 1 b 1 α 2 τ P 12 + α 3 τ P 13 -1 k<l 1 τ P kl (α k ρ l c 2 l + α l ρ k c 2 k ) (t = 0) (53)
Eventually, when assuming vanishing phase k = 3, this allows to retrieve the result of Appendix 2 for two-phase flows.

Moreover, as emphasized in [START_REF] Boukili | Relaxation and simulation of a barotropic three-phase flow model[END_REF], stable oscillations wrt time may occur in three-phase flow models, owing to the structure of B(W).

From a practical point of view, Property 6 allows to retrieve the pressure differences, yet this is not sufficient to fully compute the pressure relaxation step. For this purpose, we need to define the evolution equations through the pressure relaxation step for more variables:

Property 7:

In the pressure relaxation step (50), we have:

• The specific entropies evolution is governed by:

   ∂ t s 1 = a 1 m 1 3 k=1 p k S α k (W) ∂ t s 2 = 0 ∂ t s 3 = 0 (54)
and specific entropies s k remain in the admissible range.

• We set Π = α 1 α 2 α 3 , δ kl = α k α l , we have:        ∂ t Π = (a(α 2 -α 1 ) + b(α 3 -α 2 ) + c(α 1 -α 3 )) Π ∂ t δ 12 = (a(α 2 -α 1 ) + (b -c)α 3 ) δ 12 ∂ t δ 13 = (c(α 1 -α 3 ) + (a -b)α 2 ) δ 13 ∂ t δ 23 = (b(α 3 -α 2 ) + (c -a)α 1 ) δ 23 (55) where: a = p 1 -p 2 τ P 12 Π 0 b = p 2 -p 3 τ P 23 Π 0 c = p 3 -p 1 τ P 13 Π 0 .

Sketch of proof:

We have:

s k = s k (p k , ρ k ), this gives: ∂ t s k = ∂s k ∂p k ∂ t p k + ∂s k ∂ρ k ∂ t ρ k .
We use then the evolution equations of p k and ρ k to obtain the result (54).

Assuming that the initial data is admissible, it is obvious that s 2 and s 3 are in the admissible range. For s 1 , we use [START_REF] Crouzet | Validation of a two-fluid model on unsteady water-vapour flows[END_REF] and get:

∂ t s 1 = a 1 m 1 1≤l<k≤3 α k α l τ P kl Π 0 (p k -p l ) 2 ≥ 0
s 1 is then, as well, admissible. Eventually, (55) is obtained from (50a), through direct computations.

Therefore, numerical schemes is built in order to preserve the previous properties, more details are given in the numerical scheme subsection 3.3.2.

Temperature relaxation

In this paragraph, we examine the heat transfer subsystem. Similarly to paragraphs 2.2.5 and 2.2.6, we consider only the time derivatives and the heat transfer source terms. The studied system writes:

       ∂ t α k = 0 ∂ t m k = 0 ∂ t (m k U k ) = 0 ∂ t (α k E k ) = 3 l=1,l =k K T kl (W)(T l -T k ) (56)
Taking into account the invariance of the partial masses and the partial kinetic energies, the energy balance of (56) can be rewritten as:

m k ∂ t e k = 3 l=1,l =k K T kl (W)(T l -T k ) (57)
For the Stiffened Gas EOS we know that:

T k = 1 C V,k e k - Π k ρ k (58)
where Π k is the constant such that:

p k + γ k Π k = (γ k -1)e k ρ k .
This means that:

∂ t e k = C V,k ∂ t T k (59) since ∂ t ρ k = 0.
Therefore, (57) becomes:

∂ t T k = 3 l=1,l =k K T kl (W) m k C V,k (T l -T k ) (60)
We apply the same process that we used in paragraph 2.2.5 (Velocity relaxation); by considering the temperature differences T 1 -T 2 and T 1 -T 3 we have:

∂ t T 1 -T 2 T 1 -T 3 = -A(W) T 1 -T 2 T 1 -T 3 (61)
where the matrix A(W) is defined by:

A(W) = K T 12 m 1 C V,1 + K T 12 m 2 C V,2 + K T 23 m 2 C V,2 K T 13 m 1 C V,1 - K T 23 m 2 C V,2 K T 12 m 1 C V,1 - K T 23 m 3 C V,3 K T 13 m 1 C V,1 + K T 13 m 3 C V,3 + K T 23 m 3 C V,3 (62)
We can then formally write:

T 1 -T 2 T 1 -T 3 (t) = exp - t t 0 A(s)ds T 1 -T 2 T 1 -T 3 (t 0 ) (63)
It should be noted that the determinant and the trace of the matrix A are positive:

det(A) > 0 tr(A) > 0
A direct calculation proves it. This means that matrix A eigenvalues are either real positive, or complex with positive real part, which ensures that temperature differences decrease over time, their decrease is monotonous when eigenvalues are real, and oscillating when eigenvalues are complex.

Mass transfer

In this paragraph, we focus on the mass transfer, or Gibbs potential relaxation. In practice, we study what happens in the continuous framework when the space derivatives are null, and only the mass transfer source terms are at stake. The concerned PDE subsystem writes as follows:

         ∂ t α k = 0 ∂ t m k = 3 l=1,l =k Γ kl (W) ∂ t (m k U k ) = 3 l=1,l =k Ũkl Γ kl (W) ∂ t (α k E k ) = 3 l=1,l =k Hkl Γ kl (W) (64)
In the sequel of this work, only the mass transfer between phases 2 and 3 will be considered. This is driven by the fact that our objective is the simulation of vapour explosion, where a liquid metal gets in contact with liquid water, a part of which is eventually evaporated. The phase 1 in our study represents the liquid metal, while phases 2 and 3 represent respectively the liquid water and water vapour. Therefore, the mass transfer can happen only between phases 2 and 3. In other words:

Γ 12 = Γ 13 = 0
In practice, this means that phase k = 1 state variables don't move within this step, the effective PDE system to solve is the following:

                         ∂ t α k = 0 (k = 1, 2, 3) ∂ t m 2 = Γ 23 (W) ∂ t (m 2 U 2 ) = Ũ23 Γ 23 (W) ∂ t (α 2 E 2 ) = H23 Γ 23 (W) ∂ t (m 2 + m 3 ) = 0 ∂ t (m 2 U 2 + m 3 U 3 ) = 0 ∂ t (α 2 E 2 + α 3 E 3 ) = 0 (65a) (65b) (65c) (65d) (65e) (65f) (65g)
It is important to notice here that the mass transfer term Γ 23 (W) could be expressed only in function of the densities and internal energies:

Γ 23 (W) = Γ 23 (ρ 2 , e 2 , ρ 3 , e 3 ) (66)
Moreover, considering the closures of Ũ23 and H23 given in ( 16) and ( 19), we can easily prove that:

∂ t (m 2 e 2 ) = ∂ t (m 3 e 3 ) = 0 (67)
The result (67) means that we can express Γ 23 as a function of only one variable m 2 , since:

         ρ 2 = m 2 α 2 (t=0) e 2 = (m 2 e 2 )(t=0) m 2 ρ 3 = (m 2 +m 3 )(t=0)-m 2 α 3 (t=0) e 3 = (m 3 e 3 )(t=0) (m 2 +m 3 )(t=0)-m 2 (68) Therefore: Γ 23 (ρ 2 , e 2 , ρ 3 , e 3 ) = m 2 (M -m 2 ) τ m 23 M Γ23 (m 2 ) (69)
where M = m 2 + m 3 is constant, and noting:

Γ23 (m 2 ) = 1 Γ 0 µ 2 T 2 - µ 3 T 3 (m 2 ) ( 70
)
This approach is the one introduced in [START_REF] Crouzet | Validation of a two-fluid model on unsteady water-vapour flows[END_REF][START_REF] Hurisse | Simulation des écoulements industriels diphasiques compressibles[END_REF].

As shown in subsection 2.2.1, an entropy inequality exists for the global PDE system (1). It is notable that a similar entropy result holds for each one of the studied subsystems, namely:

• the convective subsystem (29);

• the velocity relaxation subsystem (45);

• the pressure relaxation subsystem (50).

• the temperature relaxation subsystem (56).

• the mass transfer subsystem (65).

Moreover, each one of the subsystems ((45), ( 50), ( 56) and ( 65)) guarantees the conservation of the mass, the total momentum and the total energy. This conservation is a substantial point of the relaxation processes, and will be used to build the approximation algorithms detailed in the next section.

Numerical method

We consider a classical Finite Volumes formulation, where the computational domain is meshed using unstructured 3D cells denoted Ω i , the volume of which is denoted ω i . S ij stands for the surface of the interface between cells Ω i and Ω j , n ij is the normal vector pointing from Ω i towards Ω j . We define ∆t n the time step such that: t n+1 = t n + ∆t n .

Fractional step method

In the spirit of [START_REF] Hérard | A fractional step method to compute a class of compressible gas-liquid flows[END_REF], the time scheme is the following:

• Step 1: Evolution step
For a given initial condition W n i we compute an approximate solution of W at time t n+1 , namely W n+1,- i , by solving the homogeneous part of the system (1):

       ∂α k ∂t + V i (W).∇α k = 0 ∂m k ∂t + ∇.(m k U k ) = 0 ∂m k U k ∂t + ∇.(m k U k ⊗ U k + α k p k Id) + 3 l=1,l =k Π kl (W)∇α l = 0 ∂α k E k ∂t + ∇.(α k E k U k + α k p k U k ) + V i (W). 3 l=1,l =k Π kl (W)∇α l = 0 (71) • Step 2: Relaxation step Taking W n+1,- i
as initial data, we compute W n+1 i an approximate solution of the subsystem (44).

Computing the evolution step

The scheme we consider in order to compute the evolution step is the Rusanov scheme [START_REF] Rusanov | Calculation of interaction of non steady shock waves with obstacles[END_REF]. Thus, at each interface i, j separating cells Ω i and Ω j we define the numerical normal fluxes:

       F α k n (W, n ij ) = 0 F m k n (W, n ij ) = m k U k .n ij F U k n (W, n ij ) = U k .n ij m k U k + α k p k n ij F E k n (W, n ij ) = α k (E k + p k )U k .n ij (72)
and:

       2G α k n (W i , W j , n ij ) = -r ij ((α k ) j -(α k ) i ) 2G m k n (W i , W j , n ij ) = F m k n (W i , n ij ) + F m k n (W j , n ij ) -r ij ((m k ) j -(m k ) i ) 2G U k n (W i , W j , n ij ) = F U k n (W i , n ij ) + F U k n (W j , n ij ) -r ij ((m k U k ) j -(m k U k ) i ) 2G E k n (W i , W j , n ij ) = F E k n (W i , n ij ) + F E k n (W j , n ij ) -r ij ((α k E k ) j -(α k E k ) i ) (73)
where r ij is defined by:

r ij = max k=1,2,3 ((|U k .n ij | + c k ) i , (|U k .n ij | + c k ) j ) (74)
Hence, using the standard notation φ ij = (φ i + φ j )/2 whatever φ is, the solution of the evolution step is obtained in each cell Ω i by the following update:

                                                                       ω i (α k ) n+1,- i -(α k ) n i +∆t n   j∈V (i) G α k n (W i , W j , n ij )S ij   +∆t n (U 1 ) n i .   j∈V (i) (α k ) n ij n ij S ij   = 0 ω i (m k ) n+1,- i -(m k ) n i + ∆t n j∈V (i) G m k n (W i , W j , n ij )S ij = 0 ω i (m k U k ) n+1,- i -(m k U k ) n i +∆t n   j∈V (i) G U k n (W i , W j , n ij )S ij   +∆t n 3 l=1,l =k Π kl (W) n i   j∈V (i) (α l ) n ij n ij S ij   = 0 ω i (α k E k ) n+1,- i -(α k E k ) n i +∆t n   j∈V (i) G E k n (W i , W j , n ij )S ij   +∆t n (U 1 ) n i . 3 l=1,l =k Π kl (W) n i   j∈V (i) (α l ) n ij n ij S ij   = 0 (75)
and we have the expected result:

Property 7:
The evolution step guarantees positive values of partial masses and statistical fractions if the time step complies with the constraint:

∆t n   j∈V (i) r ij S ij   ≤ 2ω i (76)
Proof: It is classical and omitted.

Computing the relaxation step

In this step, we compute approximate solutions of the relaxation step (44), which is itself split into four sub-steps: the velocity relaxation step, the pressure relaxation step, the temperature relaxation step and the mass transfer step. It is important to note that these relaxation effects are interconnected, mainly via the energy balance and the mass transfer. However, the numerical method we adopt consists in treating them separately, in the order given above.

In what follows, we give the concrete details of the algorithms used for each relaxation sub-step.

Velocity relaxation: approximate solutions of (45)

Two algorithms were considered. The difference between them lays in the approximation of the equation (47). The first algorithm uses an exact approximation of a linearised form of (47), while the second considers an implicit Euler approximation. The other algorithm steps are identical.

In what follows we present the details of each algorithm:

Algorithm 3.3.1.1
In each cell Ω i , starting with W n+1,- i , we compute W * i by following the sequence: a) Initialize the velocity differences: X n with (X n 12 , X n 13 ) n+1,- i where:

X n 1k = (U 1 -U k ).n
Initialize also the velocity relaxation matrix defined in (48):

A(W n+1,- i );
b) Compute the exact solution (X n ) * i of the linear ODE:

∂ t X n = -A(W n+1,- i )X n (77)
at time t = t n + ∆t n , using the initial condition defined in a);

c) Compute (X τ 1 ) * i = (X τ 1 12 , X τ 1 13 ) * i and (X τ 2 ) * i = (X τ 2 12 , X τ 2 13
) * i by solving the same EDO (77), using the initial conditions (X τ 1 12 , X τ 1 13 ) n+1,- i and (X τ 2 12 , X τ 2 13 ) n+1,- i respectively;

d) Compute (U 1 ) * i by the total momentum conservation:

(U 1 ) * i = 3 k=1 (m k U k ) n+1,- i + (m 2 ) n+1,- i (X 12 ) * i + (m 3 ) n+1,- i (X 13 ) * i (m 1 + m 2 + m 3 ) n+1,- i (78)
where:

(X kl ) * i = (X n kl ) * i n + (X τ 1 kl ) * i τ 1 + (X τ 2 kl ) * i τ 2 e) Update (U 2 ) * i and (U 3 ) * i : (U 2 ) * i = (U 1 ) * i -(X 12 ) * i ; (U 2 ) * i = (U 1 ) * i -(X 13 ) * i (79)
f) Finally, update the energies by the discrete counterpart of (45d):

(α k E k ) * i = (α k E k ) n+1,- i + ∆t n 3 l=1,l =k e kl (W n+1,- i ) 2 ((U l ) * i ) 2 -((U k ) * i ) 2 (80) Algorithm 3.3.1.2
In each cell Ω i , starting with W n+1,- i , we compute W * i by following the sequence: a) Initialize the velocity differences: X n with (X n 12 , X n 13 ) n+1,- i , and the velocity relaxation matrix defined in (48):

A(W n+1,- i ); b) Compute (X n ) * i such that: (Id + ∆t n A(W n+1,- i ))(X n ) * i = (X n ) n+1,- i (81)

Remark 3:

This numerical scheme is consistent with the mass conservation, the total momentum conservation and the total energy conservation.

Remark 4:

Obviously, when considering instantaneous velocity relaxation, (77)/ (81)-( 78)-( 79)-(80) degenerates into:

(U k ) * i = 3 l=1 (m l U l ) n+1,- i 3 l=1 (m l ) n+1,- i (82)

Pressure relaxation: approximate solutions of (50)

Here we present two different algorithms that were implemented, in order to compute approximate solutions of (50). The basic unknown is (α 1 , α 2 , p 1 , p 2 , p 3 ).

∂ t Y = - 1 Π 0 B(W * i )Y (
(s 2 ) n+1 i = (s 2 ) * i (s 3 ) n+1 i = (s 3 ) * i (84)
and (s 1 ) n+1 i as approximate solution at time t = t n + ∆t n of the ODE:

∂ t s 1 = (a 1 ) * i (m 1 ) * i 1≤l<k≤3 α k α l τ P kl Π 0 * i (p k -p l ) n+1 i 2 (85)
using Euler scheme.

d) Compute the pressure (p 1 ) n+1 i solution of g(x) = 1, where:

g(x) = (m 1 ) * i ρ 1 (x, (s 1 ) n+1 i ) + (m 2 ) * i ρ 2 (x -(y 12 ) n+1 i , (s 2 ) n+1 i ) + (m 3 ) * i ρ 3 (x -(y 13 ) n+1 i , (s 3 ) n+1 i ) (86)
The equation g(x) = 1 is nothing but the saturation condition of the statistical fractions (3); e) Update the pressures of phases, for k = 2, 3 compute:

(p k ) n+1 i = (p 1 ) n+1 i -(y 1k ) n+1 i (87)
f) Update the statistical fractions, for k = 1, 2, 3 compute:

(α k ) n+1 i = (m k ) * i ρ k ((p k ) n+1 i , (s k ) n+1 i ) (88)
g) Finally, update the total energies by combining (50a) and (50d):

(α k E k ) n+1 i = (α k E k ) * i + ∆t n 3 l=1,l =k (Π kl (W)) n+1 i (S α l (W)) n+1 i (89)
where (Π kl (W)) n+1 i and (S α l (W)) n+1 i are expressed, owing to [START_REF] Coquel | A positive and entropy-satisfying finite volume scheme for the Baer Nunziato model[END_REF] and [START_REF] Crouzet | Validation of a two-fluid model on unsteady water-vapour flows[END_REF], in terms of (p k ) n+1 i and (α k ) n+1 i . We note that:

k (m k e k ) n+1 -(m k e k ) * = ∆t n (p 2 + p 3 ) n+1 ((S α 1 ) n+1 + (S α 2 ) n+1 + (S α 3 ) n+1 ) = 0
The well-posedness of this algorithm relies on the following property:

Property 8:

The equation g(x) = 1 is well-posed and admits a unique solution in the admissible range.

Proof:

Let D k be the definition domain of the density ρ k (p k , s k ) wrt its first parameter p k . Thus, if we denote D g the definition domain of the function g we have:

D g = {x / x ∈ D 1 and x -(y 12 ) n+1 i ∈ D 2 and x -(y 13 ) n+1 i ∈ D 3 }
For any type of EOS, and for any x ∈ D g we have:

g (x) = - 3 k=1 m k ρ 2 k ∂ρ k ∂p k s k
We assume that the EOS complies with the following constraints:

• Positive derivative wrt the pressure:

∂ρ k ∂p k s k > 0 • Asymptotic limit: lim p k →+∞ ρ k (p k , s k ) = +∞ • Null density: ∃p 0 k ∈ D k /ρ k (p 0 k , s k ) = 0
Thus we have on the one hand:

∀x ∈ D g : g (x) < 0
On the other hand: lim

∂ t Y = - 1 Π 0 B(W * i )Y ( 90
)
at time t = t n + ∆t n , using the initial condition Y * i . c) Compute the volume fraction variables Π and δ kl with respect to (55), i.e. Π n+1 i and (δ kl ) n+1 i as the exact solutions at time t = t n + ∆t n of the ODEs:

       ∂ t Π = [a((α 2 ) * i -(α 1 ) * i ) + b((α 3 ) * i -(α 2 ) * i ) + c((α 1 ) * i -(α 3 ) * i )] Π ∂ t δ 12 = [a((α 2 ) * i -(α 1 ) * i ) + (b -c)(α 3 ) * i ] δ 12 ∂ t δ 13 = [c((α 1 ) * i -(α 3 ) * i ) + (a -b)(α 2 ) * i ] δ 13 ∂ t δ 23 = [b((α 3 ) * i -(α 2 ) * i ) + (c -a)(α 1 ) * i ] δ 23 (91)
where: a =

(y 12 ) n+1 i τ P 12 Π 0 b = (y 13 ) n+1 i -(y 12 ) n+1 i τ P 23 Π 0 c = -(y 13 ) n+1 i τ P 13 Π 0 d) Update the volume fractions (α k ) n+1 i and the densities (ρ k ) n+1 i : (α k ) n+1 i = Π n+1 i (δ jl ) n+1 i j, l = k (92) (ρ k ) n+1 i = (m k ) * i (α k ) n+1 i (93)
e) Compute the pressure (p 1 ) n+1 i solution of h(x) = 0 where:

h(x) =(m 1 ) * i e 1 (x, (ρ 1 ) n+1 i ) + (m 2 ) * i e 2 (x -(y 12 ) n+1 i , (ρ 2 ) n+1 i ) + (m 3 ) * i e 3 (x -(y 13 ) n+1 i , (ρ 3 ) n+1 i ) - 3 k=1 (m k ) * i (e k ) * i (94)
The equation h(x) = 0 is nothing but the conservation of the total internal energy during the pressure relaxation sub-step; f) Update the pressures of phases, for k = 2, 3 compute:

(p k ) n+1 i = (p 1 ) n+1 i -(y 1k ) n+1 i (95)
g) Finally, update the total energies by combining (50a) and (50d):

(α k E k ) n+1 i = (α k E k ) * i + ∆t n 3 l=1,l =k (Π kl (W)) n+1 i (S α l (W)) n+1 i ( 96
)
where (Π kl (W)) n+1 i and (S α l (W)) n+1 i are expressed, owing to [START_REF] Coquel | A positive and entropy-satisfying finite volume scheme for the Baer Nunziato model[END_REF] and [START_REF] Crouzet | Validation of a two-fluid model on unsteady water-vapour flows[END_REF], in terms of (p k ) n+1 i and (α k ) n+1 i . Remark 5: Algorithms 3.3.2.1 and 3.3.2.2. are consistent with the mass conservation, the total momentum conservation and the total energy conservation. Additionally, alternative pressure relaxation algorithms could be built when replacing ODEs (83) and (90) by their first order implicit Euler counterparts, similarly to velocity relaxation algorithms (where (81) is the implicit Euler counterpart of ODE (77)).

Remark 6:

In the asymptotic situation where the user would assume infinite drag effects between phases, thus neglecting relative velocities (U 1 = U 2 = U 3 = U where U is defined in (82)), we may consider the limit case where p 1 = p 2 = p 3 = P where P is given in appendix A1.

Temperature relaxation: approximate solutions of (56)

Taking into account the continuous framework study presented in paragraph 2.2.7, we consider the following algorithm to compute approximate solutions of (56). In each cell Ω i , starting from a state W * i , and for a given time step ∆t n , we follow the sequence:

• Compute (T 1 -T 2 ) n+1 i and (T 1 -T 3 ) n+1 i
, the exact solutions of the ODE:

∂ t T 1 -T 2 T 1 -T 3 = -A(W * ) T 1 -T 2 T 1 -T 3 (97)
at time t n + ∆t n ;

• Compute (T 1 ) n+1 i , by solving the total energy conservation equation:

(T 1 ) n+1 i = 1 3 k=1 (m k ) * i C V,k 3 k=1 (m k ) * i C V,k (T k ) * i + 3 k=2 (m k ) * i C V,k (T 1 -T k ) n+1 i ( 98
)
and deduct (T 2 ) n+1 i and (T 3 ) n+1 i by: (T k )

n+1 i = (T 1 ) n+1 i -(T 1 -T k ) n+1 i , k = 2, 3 ;
• Update the total energies by:

(α k E k ) n+1 i = 1 2 m k U k .U k * i + (m k e k ) n+1 i ( 99
)
where :

(m k e k ) n+1 i = (m k ) * i e k ((ρ k ) * i , (T k ) n+1 i
) as the masses and the densities are constant in the temperature relaxation step.

Mass transfer: approximate solutions of (65)

For computing approximate solutions of (65), the numerical scheme that we consider is the one which was presented in [START_REF] Crouzet | Validation of a two-fluid model on unsteady water-vapour flows[END_REF][START_REF] Hurisse | Simulation des écoulements industriels diphasiques compressibles[END_REF]. In each cell Ω i , starting from a state W * i , and for a given time step ∆t n , we compute W n+1 i by following the sequence:

• Compute (m 2 ) n+1 i , solution of: (m 2 ) n+1 i = (m 2 ) * i + ∆t n (M τ m 23 ) * (m 2 ) n+1 i (M -(m 2 ) n+1 i ) Γ23 ((m 2 ) n+1 i ) ( 100
)
where:

M = (m 2 ) * i + (m 3 ) * i • Compute (m 3 ) n+1 i : (m 3 ) n+1 i = M -(m 2 ) n+1 i (101)
• Compute (U 2 ) n+1 i and (U 3 ) n+1 i by solving the linear system wrt (U 2 , U 3

) n+1 i : (m 2 U 2 ) n+1 i = (m 2 U 2 ) * i + ( Ũ23 ) n+1 i ∆tn (M τ m 23 ) * (m 2 ) n+1 i (M -(m 2 ) n+1 i ) Γ23 ((m 2 ) n+1 i ) (m 3 U 3 ) n+1 i = (m 2 U 2 ) * i + (m 3 U 3 ) * i -(m 2 U 2 ) n+1 i (102)
• Update the total energies:

(α 2 E 2 ) n+1 i = (α 2 E 2 ) * i + ( H23 ) n+1 i ∆tn (M τ m 23 ) * (m 2 ) n+1 i (M -(m 2 ) n+1 i ) Γ23 ((m 2 ) n+1 i ) (α 3 E 3 ) n+1 i = (α 2 E 2 ) * i + (α 3 E 3 ) * i -(α 2 E 2 ) n+1 i ( 103
)
Remark 7:

• For more details about the resolution of equation (100) we refer the reader to the previous works [START_REF] Crouzet | Validation of a two-fluid model on unsteady water-vapour flows[END_REF][START_REF] Hurisse | Simulation des écoulements industriels diphasiques compressibles[END_REF] ;

• By construction, this scheme ensures the conservation of the total mass, the total momentum and the total energy.

Numerical results

This section is dedicated to the presentation of the different numerical test cases that were conducted using the three-phase flow model presented in Section 2 and the numerical schemes and algorithms given in Section 3.

In the first subsection, we present a verification test case where we examine the computation of the evolution step, this consists in studying the convergence rates in terms of the mesh size on a Riemann problem. In the second subsection, we focus on the relaxation step and we present a series of verification test cases which correspond to each relaxation subsystem. Then, in section 4.3 we consider a two-phase shock-tube test case, which studies the impact of a planar shock wave on a cloud of deformable droplets. Finally, in subsection 4.4 a more complete test case representative of a vapour explosion scenario is presented.

Verification of the evolution step

This subsection is devoted to the verification of the convective subset of (1) using the numerical scheme (75). We consider a Riemann problem where two waves are at stake: the contact discontinuity associated with the eigenvalue λ 1,2,3 (w) = u 1 and a right-going shock wave of phase 1.

In order to initialize the Riemann problem, we start with a given left state w L , then we use the Riemann invariants introduced in (37) to build the intermediate state w int such that: We emphasize that for this test case, all three phase EOS are chosen to be Ideal Gas EOS with the following thermodynamic parameters:

I i 1,2,3 (w L ) = I i 1,
Phase 1 Phase 2 Phase 3 γ k 1.4 1.7 1.6
Calculations are run with 6 different mesh sizes h: 1/5000, 1/10000, 1/25000, 1/50000, 1/100000 and 1/200000. The computational domain is the unit cell [0, 1] and the Riemann problem is initialized in x = 0.5.

Figure 1 gives the error in L 1 norm measured on the system's principal variables. We note that for coarse meshes, the u 1 curve shows a high convergence rate close to h 1 , which is the asymptotic rate for shock waves. As the mesh is refined, the convergence rate tends towards h 1/2 . That is the rate corresponding to the contact discontinuity associated with the coupling wave λ = u 1 . For the other variables, the h 1/2 rate is obtained more quickly.

In order to get a better understanding of the structure of fields involved in this test case, figure 2 gives the spacial representation of state variables at the final time of the simulation T = 0.8ms, while focusing on the 1/50000 mesh results.

We can observe the two jumps related to:

• The contact discontinuity located around x contact = 0.74. In fact, considering the initial data, the coupling wave speed is:

σ contact = u L 1 = 300.0 m/s (104)
Knowing that the Riemann problem was initialized at x Riemann = 0.5, the observed jump's position x contact is in agreement with the theoretically expected result: x Riemann + T σ contact ;

x contact = 0
• The phase 1 shock wave observed around x shock = 0.95. Actually, the initial data give a right-going shock the speed of which is:

σ shock = (ρ 1 u 1 ) R -(ρ 1 u 1 ) int ρ R 1 -ρ int 1 = 559.68442 m/s (105)
We check that all variables of phases 2 and 3 do not jump at this position, which is in agreement with shock relations. Of course, the statistical fraction α 1 does not jump, as expected.

Verification of the relaxation schemes 4.2.1. Velocity relaxation verification test cases

This paragraph provides the results of two verification test cases of the velocity relaxation sub-step. It consists in verifying the algorithms introduced in the subsection (3.3.1), to find approximate solutions of (45). Appendix A3 gives the details of the analytical solution in each case.

• Test case A3.1 (see figure 3)

The analytical solution in this case is in the discrete approximation space of Algorithm 3.3.1.1. In fact, by choosing the coefficients e kl (introduced in ( 14)) constants, the approximation (77) becomes an exact counterpart of (129), which means that Algorithm 3.3.1.1 in this case computes the exact analytical solution. The error is at its least when the mesh is coarse, here its value is around 10 -16 (round-off error). Then, the error starts to increase progressively as the mesh is refined. This is due to the cumulation of the rounding errors throughout the different algorithm steps.

Turning then to Algorithm 3.3.1.2, an expected convergence at the rate h 1 is observed. • Test case A3.2 (see figure 4)

In this case, the coefficients e kl are chosen to be time-dependent. This means that none of the two algorithms computes the exact solution, both compute approximations, thus the error is decreasing as the mesh is refined. Once more, the convergence rate tends towards the expected h 1 rate.

Pressure relaxation verification test case

We turn in this paragraph to the verification of the schemes introduced in (3.3.2), to obtain approximate solutions of the pressure relaxation sub-system (50). In appendix A4 we present the analytical solution that was considered, and in figure 5 we give the obtained numerical results with both algorithms (3.3.2.1) and (3.3.2.2).

We note that the Algorithm 3.3.2.2 shows a convergent behaviour, starting by a constant error level for coarse meshes, then it arrives at a point where the mesh gets sufficiently refined (mesh size dt ≈ 10 -4 ) and thus the convergence could be observed. This computational result was obtained with τ0 = 10 2 (τ 0 is defined in (140)), which corresponds to:

dt τ0 /(γK 1 ) ≈ 1 (106)
(See the property A4.1 of Appendix A4 for the definition of K 1 ). Other numerical simulations were run with different τ0 , γ and K 1 values, we noticed that (106) hold for all of them, this means that the convergence start to be observable when the mesh size dt becomes small enough compared to the pressure relaxation characteristic time included in τ0 γK 1 .

Temperature relaxation verification test case

In this paragraph, we implement a test case that aims at the verification of the scheme introduced in subsection 3.3.3 to compute approximate solutions of (56). Appendix A7 provides the considered analytical solution, as well as the initial dataset that was used in the computations. Figure 6 shows that the variables converge at the expected h 1 rate.

Mass transfer verification test case

The mass transfer scheme was studied previously in [START_REF] Crouzet | Validation of a two-fluid model on unsteady water-vapour flows[END_REF][START_REF] Hurisse | Simulation des écoulements industriels diphasiques compressibles[END_REF]. For more details about the analytical solution considered in this verification test case, we refer the reader to [START_REF] Crouzet | Validation of a two-fluid model on unsteady water-vapour flows[END_REF] and its Appendix A. This test's EOS framework is the Stiffened Gas EOS:

p k + γ k Π k = (γ k -1)ρ k e k s k = C V k log e k -Π k ρ k ρ 1-γ k k + q k (107)
In order to be able to build the analytical solution, one possibility is that phases 2 and 3, that are really concerned by the mass transfer, must comply with:

γ 2 C V 2 = γ 3 C V 3 (108)
Table 1 gives the numerical data that was used in the test case. We mention that the test's final time is T f inal = 1.0 s, and the mass transfer relaxation time scale is fixed to: We check also that EOS of phases 2 and 3 satisfy (108).

τ m 23 = 5.0 × 10 +3 Phase k = 1 Phase k = 2 Phase k = 3 u k 0.
Turning to the results, figure 7 gives the error in function of the mesh size. We check that the expected h 1 convergence rate is retrieved.

Mass transfer verification

First validation test case

In this subsection we focus on the shock-tube experiment presented in [START_REF] Chauvin | Étude expérimentale de l'atténuation d'une onde de choc par un nuage de gouttes et validation numérique[END_REF][START_REF] Chauvin | Experimental investigation of the propagation of a planar shock wave through a two-phase gas-liquid medium[END_REF]. It consists in the analysis of the effects of a planar air shock wave on a cloud of water droplets. A highpressure chamber filled with air generates the planar shock wave, which moves rightwards into the air at atmospheric pressure. The cloud of droplets is placed in the low-pressure chamber and undergoes the impact of the incident shock wave.

In order to get a correct understanding of the shock wave / droplets interaction, several pressure sensors have been placed in the low-pressure chamber. This allows to record the pressure histories in different tube positions, and therefore highlight the undergoing physical mechanisms. Figure 8 gives the geometrical details of the shock-tube apparatus, as well as the droplets and the pressure sensors' positions. In [START_REF] Chauvin | Étude expérimentale de l'atténuation d'une onde de choc par un nuage de gouttes et validation numérique[END_REF][START_REF] Chauvin | Experimental investigation of the propagation of a planar shock wave through a two-phase gas-liquid medium[END_REF], many experimental settings have been presented: a single-phase shock-tube without droplets, and three different two-phase gas-liquid mixtures where different cloud sizes were considered. It should be noted that the thermal effects don't play the primary role in this experiment, instead it is the dynamics (including velocity and pressure effects) that have the most important influence. Therefore, temperature relaxation and mass transfer will be neglected in the numerical simulations of this experimental setup.

In fact, when put together, the convective subsystem, the velocity and the pressure relaxation effects allow to get an adequate simulation of the experiment's dynamics. In addition to that, an interfacial area equation was included in order to take into account the atomization effect. The details of the interfactial area equation can be found in [START_REF] Boukili | Relaxation and simulation of a barotropic three-phase flow model[END_REF]. In the continuity of the latter work, we adopt the same initialization of the numerical test cases (namely for the pressures, the velocities and the statistical fractions). The difference meanwhile concerns the EOS: while [START_REF] Boukili | Relaxation and simulation of a barotropic three-phase flow model[END_REF] focuses on the barotropic framework, in the present work we adopt a non-barotropic EOS. Appendix A5 presents the different EOS as well as the test case initialization datasets.

In the sequel, we discuss the comparison between the experimental and the obtained numerical results:

• Single-phase flow configuration: This is a reference test case where we consider a single-phase shock tube configuration. No particle lid is included (α 1 is set to 10 -10 in all the computational domain) and the mean pressure is denoted: Figure 9 gives the pressure history recorded in station 4. We denote P * the pressure level just behind the right-going shock wave observed around t = 3.7 ms, and P * * the pressure level after the reflection of this shock wave on the right wall boundary of the tube.

P mix = 3 k=1 α k P k (109)
In Table 2 we give the comparison between the different experimental, analytical and numerical results, including the barotropic case that was presented in [START_REF] Boukili | Relaxation and simulation of a barotropic three-phase flow model[END_REF]. We note, on the one hand, that the numerical results on a sufficiently refined mesh (50000 cells per meter here) are in agreement with the analytical solutions, both for the barotropic and the full Euler model (full here means the model with energy equations), which is expected and mandatory. On the other hand, we note that the full Euler model enables to obtain a better approximation of the experimental results. This has been already mentioned in [START_REF] Boukili | Relaxation and simulation of a barotropic three-phase flow model[END_REF], in figure 9 we give the corresponding numerical result when using the model with energy. The shock wave celerity is also an interesting aspect to observe. In Table 2 we note that the experimental and the numerical (or exact) values display a 6% difference. Nonetheless, it is not sufficient to explain the gap shown in figure 9, where the numerical reflected wave returns back more quickly than the experimental wave, more than what could correspond to the latter 6% celerity discrepancy. This delay of the experimental reflected wave could be partially justified by the friction on the tube's sidewall and end-wall, which is not taken into account in our model. It could also be the consequence of a geometric difference, between the numerical computational domain (figure 18) and the real tube used in the experimental setup [START_REF] Chauvin | Étude expérimentale de l'atténuation d'une onde de choc par un nuage de gouttes et validation numérique[END_REF][START_REF] Chauvin | Experimental investigation of the propagation of a planar shock wave through a two-phase gas-liquid medium[END_REF].

Another notable conclusion of this single-phase comparison is the fact that the full Euler model provides results that are already slightly different from the experimental observations. Hence, we may expect some greater discrepancies when we run two-phase or three-phase simulations.

• Gas-liquid two-phase flow configuration:

This test corresponds to the interaction of a shock wave with a lid of liquid deformable droplets. The initialization is identical to the gas-solid case, but we take into account the interfacial area equation A as presented in [START_REF] Boukili | Relaxation and simulation of a barotropic three-phase flow model[END_REF]. Figure 10 gives the numerical results as well as the experimental data. We note that, at the beginning, the numerical simulations show a similar pattern to the gas-solid case, with an incident shock wave of 2.4 bar amplitude on station 2, and a smaller amplitude on station 3. But the major difference comes from the sudden decrease that follows the arrival of the incident shock, which is clearly visible on station 3. This is the characteristic signature of the droplets atomization [START_REF] Gelfand | Droplet breakup phenomena in flows with velocity lag[END_REF] under the effect of the in-coming shock wave: each droplet breaks into smaller droplets until reaching a certain diameter, where the particle cloud reaches an equilibrium with the incident shock wave. Afterwards, the pressure starts to build up and a compression is observed once more.

On station 2, the fragmentation signature can hardly be observed, and the compression following the incident shock occurs quickly: rapidly a pressure plateau is observed from t ≈ 2.8 ms to t ≈ 6 ms. The end of the plateau is marked by a pressure decrease due the arrival from the left wall boundary of the reflected rarefaction wave.

As expected from the single-phase comparison, the numerical pressure levels are greater than the measured experimental results. Nevertheless, it seems to be remarkable that the qualitative pressure behaviour is rather good, as it takes into account the major dynamical phenomena that occur in this fragmentation experiment. The pressure levels also are closer to the experimentally measured values than those arising from the barotropic model computation (as observed in [START_REF] Boukili | Relaxation and simulation of a barotropic three-phase flow model[END_REF]). In order to get a more complete view on this fragmentation test case, we provide here some additional numerical results, though the comparison with experimental data is not available. Figure 11 gives the evolution of the relative velocity wrt time in stations 2 and 3. We can notice that the relative velocity jumps and it reaches a maximum value right at the arrival of the incident shock wave, which in turn activates the fragmentation process. After that, it decreases rather quickly. In both stations, the relative velocity decrease could be put in match with the pressure build-up after the incident shock wave. Both the pressure and the relative velocity reflect the nature of the fragmentation.

It is worth to be noted that the phase 1 momentum behaviour is consistent with the evolution of α 1 and A, as shown in figure 12. This allows to understand that the particle cloud, in addition to the fragmentation, tends to move rightwards under the effect of the incident right-going shock wave. This movement is not uniform, the left side of the cloud (located initially around x = 3.0) undergoes a more significant shifting than the right side (initially at x = 3.4).

Finally, figure 13 gives the comparison between the numerical results obtained on different meshes. We notice that the gap between the results gets smaller and smaller as the mesh is refined. Therefore, the discussed 1/50000 mesh size results above are reasonably representative of the model solution.

Second validation test case: the KROTOS experimental setup

The KROTOS experimental program was conducted in order to study some aspects of fuel coolant interactions (FCI), by performing research tests on specific phenomena involved in a FCI scenario. Vapour explosion is one of the aspects that the KROTOS experiment examined, by considering several fuel-coolant mixtures where different types of nuclear reactor prototypic materials were involved (zirconium dioxide, aluminium oxide, uranium dioxide, tin, etc.), the liquid water served as a standard coolant. Various conditions have been tested, including fuel masses between 1 kg and 10 kg at temperatures around 3000 K, while the liquid water was taken either at the saturation temperature, or with low subcooling (≈ 10 K less than the saturation temperature). For more details about these tests, we refer the reader to [START_REF] Huhtiniemi | KROTOS 38 to KROTOS 44[END_REF] where detailed descriptions and results are given.

The particular configuration that we focus on is the KROTOS-44 test, where alumina melt was poured into water with low subcooling, melt droplets are formed and progress slowly inside the water volume. An external pressure pulse was used to trigger the explosion before the first droplet hits the bottom of the water pool, at the triggering instant we consider that all the melt droplets are located inside the water volume. This configuration was identified as a numerical benchmark to test and compare different FCI codes, including MC3D [START_REF] Meignen | The challenge of modelling fuel-coolant interaction. Part II: steam explosion[END_REF] which is commonly used in France for vapour explosion simulations. [33] gives the details of the imposed geometry and the melt-coolant distribution, at the explosion's initiation instant. Based on this KROTOS-44 experiment, we have adopted the setup given in figure 14 for our numerical simulation.

In the considered setup of figure 14, the phases distribution in the low pressure chamber is identical to what is given in [33]: the chamber is filled of liquid water at 363 K which is present as a pure phase between x = 2.00 and x = 2.15, then there is a 70 cm thick mixed zone called the interaction zone, where the three phases (melt droplets, liquid water and its vapour) are present, after which there is the plug composed of liquid and vapour water, and finally for x > 3.23 a pure gaseous volume covering the whole. The low pressure chamber was equipped with pressure stations at different positions, in order to record pressure histories over time and allow the examination of the pressure's evolution.

The dataset that was used to initialize the numerical simulation is:

= 0.026 α 2 = 1 -2 × 10 -6 α 2 = 1 -2 × 10 -6 α 2 = 0.884 α 2 = 0.835 -10 -6 α 3 = 1 -2 × 10 -6 α 1 = α 3 = 10 -6 α 1 = α 3 = 10 -6 α 3 = 0.09 α 3 = 0.165 α 1 = α 2 = 10 -6 T k = 1000.0 K T k = 363.0 T 1 = 2500.0 K α 1 = 10 -6 T 3 = 700.00 K p k = 150 × 10 5 P a T 2 = 363.00 K T k = 363.00 K T 1 = T 2 = 363.00 K u k = 0
T 3 = 1000.0 K u k = 0 and p k = 1.0 × 10 5 P a In terms of EOS, we consider the General Stiffened Gas framework:

       s k = C V,k log e k -Q k -Π k ρ k 1 ρ γ k -1 k + q k p k + γ k Π k = (γ k -1)ρ k (e k -Q k ) C V,k T k = e k -Q k -Π k ρ k (110)
The adopted numerical values are given in table 4, where phases 1, 2 and 3 represent respectively the melt droplets, the liquid water and the vapour water. The initial droplet diameter is set to 15 mm, as stated in [33]. We mention that in this test case the diameter remains constant, as we did not take the fragmentation into account; the interfactial area A equation (see [START_REF] Boukili | Relaxation and simulation of a barotropic three-phase flow model[END_REF] for more details) was not included in the realized numerical simulations. In order to treat the relaxations effects, we have chosen instantaneous relaxation algorithms for both the pressure and the velocity, that is:

Phase k γ k Π k Q k C V,k q k 1 (
τ P kl = 0 τ U kl = 0 (111)
while for mass transfer, we have considered the value given in [START_REF] Lochon | Comparison of two-fluid models on steamwater transients[END_REF], and an equivalent value for heat exchange:

τ m 23 = 5.0 × 10 -5 τ T kl = τ T 0 = 10 -5 (112)
We mention also that we set for all numerical tests CF L = 1/2.

Turning to numerical results, a first aspect to examine is the mesh refinement. In fact, for a final simulation time of less than 2 ms, several meshes were tested. Figure 15 provides the comparison between the different meshes results in term of velocity, at the same time t = 0.7 ms. From figure 15 we can tell that the discrepancies between the results are reduced as mesh is refined. Thus, in the sequel we will focus on the finest mesh (375000 cells) to examine the pressure and velocity results. In figure 16, pressure histories in stations S 1 , .., S 4 are given. It presents the mean pressure P mix (defined previously in (109)) wrt time at each station. For station S 1 , we notice that the behaviour of the pressure signal at the early time instants, namely t ≈ 0.05 ms, is similar to what happens in the one-phase flow framework: a sudden pressure jump, followed by a plateau.

Turning to stations S 2 , S 3 and S 4 located inside the interaction zone, we remark that as the computation starts, the pressure increases in a uniform way over the three stations, which is translated by the fact that the three pressure curves are superimposed for t ∈ [0, 0.25] (respectively for S 3 and S 4 for t ∈ [0, 0.5]). This is the signature of the occurring evaporation inside the interaction zone (see the flow in a box analysis given in Chapter 3 of [START_REF] Boukili | Schémas de simulation d'un modèle à trois phases immiscibles pour application à l'explosion vapeur[END_REF]), which is a direct consequence of the initial thermodynamic setup of this zone. In fact, giving the temperature and Gibbs potential differences at the initial instant, mass transfer acts and creates more and more vapour, which drives the pressure to increase. At t ≈ 0.25 ms, station S 2 observes the arrival of the shock wave, while S 3 and S 4 continue the evaporation process.

Focusing on the pressure jumps created by the shock wave, from figure 16 we can notice that the jump created in station S 1 is δP ≈ 9 M P a, while in stations S 3 and S 4 it reaches the values δP ≈ 14 M P a and δP ≈ 15 M P a respectively. This shows that the pressure jump induced by the shock wave increases as this wave travels inside the interaction zone.

In terms of total pressure P mix , at t ≈ 0.5 ms, when the shock wave reaches station S 3 , the total pressure reaches P mix ≈ 17 M P a, which is higher than the triggering initial pressure 15 M P a. This behaviour is retrieved in station S 4 as well, with even a higher total pressure value P mix ≈ 18 M P a. We recall that in subsection 4.3 we have seen pressure jumps resulting from the dynamic effects (velocity and pressure relaxations) but without exceeding the triggering high pressure. Here, what is involved -in addition to dynamicsis thermal effects, namely thermal exchange and the resulting evaporation, and the result is higher pressure jumps that are increased, rather than attenuated, when the shock wave travels inside the interaction zone.

In figure 17, we provide complementary results of this test case, by presenting the pressure and velocity at five distinct time instants. The progress of the right-going shock wave can be observed, at t = 1.0 ms the front of this wave is located around x = 2.90 m, which means that at this time the shock wave has completely crossed the interaction zone. Pressure's values can be also checked, we notice that for t = 0.55 ms and t = 0.70 ms the reached pressure values are higher than 15 M P a, which confirms the observed result in stations S 3 and S 4 . In this figure, we can also see the progress of the left-going rarefaction wave, which -shortly before t = 1.0 ms reaches the left tube's boundary where it gets reflected.

Finally, it is relevant to draw a comparison between the numerical and experimental results concerning the KROTOS-44 setup. As provided in [START_REF] Huhtiniemi | KROTOS 38 to KROTOS 44[END_REF], the measured pressure upper bound around stations S 3 and S 4 locations is of the order of 50 M P a to 60 M P a, while our model's result is smaller and of the order of 18 M P a, this was expected as the breakup phenomenon was not taken into account. Nevertheless, we mention that the numerical benchmark given in [33] shows an important dispersion between the different vapour explosion codes, pressure values vary from 10 M P a to more than 100 M P a. In turn, the code

Conclusion

This paper is the continuation of research efforts concerning the modeling and the simulation of vapour explosion (VE ). A previous study of a barotropic three-phase flow model was presented in [START_REF] Boukili | Relaxation and simulation of a barotropic three-phase flow model[END_REF]. In the present paper, we have discussed a three-phase flow model including the energy balance. The model properties were highlighted, with special focus on the hyperbolicity, the entropy inequality and the jump conditions. A fractional step method was considered, we have shown that on the continuous level it allows to remain consistent with the entropy equation, and in the discrete framework it enables to get satisfactory results on fine enough meshes. The comparison with the results of the experimental setup presented in [START_REF] Chauvin | Étude expérimentale de l'atténuation d'une onde de choc par un nuage de gouttes et validation numérique[END_REF][START_REF] Chauvin | Experimental investigation of the propagation of a planar shock wave through a two-phase gas-liquid medium[END_REF] is satisfactory, which means that the model is fairly able to account for the mechanical effects that take place in such gas-droplet flow configurations. The study of a VE test case, on the basis of the KROTOS-44 experiment [START_REF] Huhtiniemi | KROTOS 38 to KROTOS 44[END_REF], is also given. The combination of dynamic effects and thermal exchanges is highlighted, and qualitatively fair enough results are retrieved.

Nevertheless, there are still some points that require more work in order to progress towards a better three-phase flow VE oriented model:

• The interfacial area equation needs to be taken into account, to translate the droplets break-up and eventually allow a better estimation of the pressure's evolution;

• It would be worth pursuing the investigation of accuracy of the non-instantaneous relaxation algorithms (for pressure and velocity). This would allow to consider a wider range of time scales;

• There is a need for more accurate and robust Riemann solvers (see for instance [START_REF] Ambroso | Relaxation and numerical approximation of a two-fluid two-pressure diphasic model[END_REF][START_REF] Ambroso | A Godunov-type method for the seven-equation model of compressible two-phase flow[END_REF][START_REF] Coquel | A positive and entropy-satisfying finite volume scheme for the Baer Nunziato model[END_REF][START_REF] Saleh | A relaxation scheme for a hyperbolic multiphase flow model[END_REF] in that direction).

Acknowledgements

The first author receives financial support by ANRT and EDF through an EDF/CIFRE grant number 2016/0611. Computational facilities were provided by EDF. Both authors would like to thank Olivier Hurisse for his help and for various useful discussions on the subject.

Appendix A1 -Instantaneous pressure relaxation

Here in this appendix, we consider the case of a pressure relaxation system where the pressures differences p k -p l are neglected. In this case, we get at once:

         p 1 = p 2 = p 3 = P ∂ t m k = 0 ∂ t (m k U k ) = 0 ∂ t (α k E k ) + P∂ t α k = 0 (113a) (113b) (113c) (113d)
The energy balance (113d) can be rewritten as:

∂ t (m k e k ) + P∂ t α k = 0
In order to compute a discrete approximate solution, we apply an implicit Euler scheme:

(m k e k ) -(m k e k ) * ∆t + P α k -α * k ∆t = 1 ∆t [(m k e k ) -(m k e k ) * + P(α k -α * k )] = 0 (114)
where the * exponent indicates the given initial data.

For further development, we consider a Stiffened Gas EOS:

(γ k -1)m k e k = α k p k + α k γ k Π k (115)
This gives:

(α k P + α k γ k Π k ) -(α k p k + α k γ k Π k ) * + (γ k -1)P(α k -α * k ) = 0 (116)
By denoting P * k = (α k p k + α k γ k Π k ) * , we express the statistical fractions in function of the pressure:

α k = P * k + (γ k -1)α * k P γ k (P + Π k ) (117)
Now, to obtain the pressure P, we solve the saturation condition:

3 k=1 α k = 1 (118)
Thus, we consider the function:

h(P) = 3 k=1 P * k + (γ k -1)α * k P γ k (P + Π k ) (119) which definition domain is D =] -x 0 , +∞[, where x 0 = Π k 0 = min k=1,2,3 (Π k ).
We have the following key property:

Property:

The equation h(P) = 1 admits a unique solution in the domain D.

Proof:

For any x ∈ D we have:

h (x) = - 3 k=1 α * k (p * k + Π k ) γ k (x + Π k ) 2 < 0 Moreover, we have: lim x→+∞ h(x) = 3 k=1 α * k (γ k -1) γ k = 1 - 3 k=1 α * k γ k < 1
and:

lim x→-x 0 + h(x) = lim x→-x 0 + 3 k=1,k =k 0 P * k + (γ k -1)α * k (-Π k 0 ) γ k (-Π k 0 + Π k ) + P * k 0 + (γ k 0 -1)α * k 0 x γ k 0 (x + Π k 0 )
A straightforward computation gives:

P * k 0 + (γ k 0 -1)α * k 0 (-Π k 0 ) = α * k 0 (p * k 0 + Π k 0 ) > 0
The three momentum equations are equivalent to:

   m 1 ∂ t u 1 = e 12 (w)(u 2 -u 1 ) + e 13 (w)(u 3 -u 1 ) m 2 ∂ t u 2 = e 12 (w)(u 1 -u 2 ) + e 23 (w)(u 3 -u 2 ) ∂ t (m 1 u 1 + m 2 u 2 + m 3 u 3 ) = 0 (123)
where the phase 3 equation has been replaced by the total momentum conservation equation.

If we denote Q 0 = m 1 u 1 (t = 0) + m 2 u 2 (t = 0) + m 3 u 3 (t = 0) the total momentum in the velocity relaxation sub-step, we can integrate the phase 3 equation and write:

u 3 (t) = 1 m 3 (Q 0 -m 1 u 1 (t) -m 2 u 2 (t)) (124)
We use (124) to replace u 3 (t) in (123), this gives: where

  Phase k = 1 Phase k = 2 Phase k = 3 u k 10.
Q 1 A 1 = Q 2 A 2 .
The velocity u 3 is given by (124), and the energies are obtained by direct integration of the energy balance of (123).

The initial dataset that was considered for this case is given in Table 5. The time relaxation scales are chosen such that: This gives:

∂ t u 1 + α(t + 1)u 1 = β(t + 1) ∂ t u 2 + α(t + 1)u 2 = β(t + 1) (132)
where α and β are time-independent constants:

α = 1 m 1 m 2 + 1 m 1 m 3 + 1 m 2 m 3 β = Q 0 m 1 m 2 m 3 (133)
This gives:

u 1 (t) = β α + u 1 (t = 0) -β α exp(-α(t 2 /2 + t)) u 2 (t) = β α + u 2 (t = 0) -β α exp(-α(t 2 /2 + t)) (134)
The velocity u 3 and the energies are obtained similarly to the first case.

Practically, for this case, the initial dataset is identical to the test case A3.1 (see Table 5).

We have also:

                         Π 0 ∂ t α 2 = - K 1 τ P 12 α 2 + γ 1 -γ 2 γ 2 -1 K 2 τ P 12 α 2-γ 2 2 + K 2 τ P 12 α 1-γ 2 2 + 1 τ P 23 - 1 τ P 12 K 2 α 1-γ 2 2 α 3 + γ 1 -1 γ 3 -1 1 τ P 12 - 1 τ P 23 K 3 α 2 α 1-γ 3 3 Π 0 ∂ t α 3 = - K 1 τ P 13 α 3 + γ 1 -γ 3 γ 3 -1 K 3 τ P 13 α 2-γ 3 3 + K 3 τ P 13 α 1-γ 3 3 + 1 τ P 23 - 1 τ P 13 K 3 α 1-γ 3 3 α 2 + γ 1 -1 γ 2 -1 1 τ P 13 - 1 τ P 23 K 2 α 3 α 1-γ 2 2 (138)
Proof:

(137a) and (137b) are nothing but the conservation of the specific entropies of phases 2 and 3. In fact, we have seen in (54) that:

∂ t s 2 = ∂ t s 3 = 0 This means that: ∂ t (p 2 ρ -γ 2 2 ) = ∂ t (p 3 ρ -γ 3 3 ) = 0 Knowing that m k = α k ρ k and ∂ t m k =
0 in the pressure relaxation step, we conclude that (137a) and (137b) hold.

(137c) is simply obtained by writing the conservation of the total energy:

∂ t (m 1 e 1 + m 2 e 2 + m 3 e 3 ) = 0
Therefore, we use the latter results (137) to develop (135a) and (135b), this allows to obtain (138) through direct calculations.

In order to explicitly develop the analytical solution, we add more assumptions:

• The three phases have an identical EOS:

γ 1 = γ 2 = γ 3 := γ > 1 (139)
• The pressure relaxation time scales are all equal to the same constant value:

∀k, l : τ kl (W) := τ P 0 and we denote: τ0 = τ P 0 Π 0 (140) Under these assumptions, we have :

given by the mathematical modelling, notably the relaxation characteristic times. Thus, we focus on the following three configurations:

• Configuration 1 :

We consider Algorithm 3.3.1.1 for the velocity relaxation, with:

   e 12 (W) = 1 8 Aρ 2 C 12 |U 1 -U 2 | e 13 (W) = m 3 m 2 e 12 (W) e 23 (W) = m 3 m 1 e 12 (W) (142)
The Closure of e 12 (see Chapter 1 subsection 4.2) gives an implicit definition of the velocity relaxation time scale τ U 12 , the chosen e 13 and e 23 allow to retrieve the same relaxation time scale τ U 12 = τ U 13 = τ U 23 = τ U while being consistent with the general form:

e kl (W) = 1 τ U kl m k m l m 1 + m 2 + m 3 (143)
For the pressure relaxation, we consider Algorithm 3.3.2.1 with the same relaxation time scale for all the phases:

τ P kl (W)Π 0 = τ P 0 = 10 -5 (144)
• Configuration 2 :

Here we consider the same relaxations algorithms of the Configuration 1, but with the characteristic times:

τ U Conf ig.2 = τ U kl /10 N τ P Conf ig.2 = τ P kl /10 N (145)
Several values of N were tested in order to understand the impact of these relaxation time scales on the relaxation procedure, the value we present in this Configuration 2 is N = 10.

• Configuration 3 :

For this configuration, we consider the instantaneous relaxation algorithms, both for the pressure and the velocity.

In order to comment the computational results, we restrict the scope of the analysis to the following three aspects: the mesh convergence, the comparison of Configuration 1 and Configuration 2, and the comparison of Configuration 2 and Configuration 3.

(i) Mesh convergence

Figure 19 gives the comparison in terms of velocity of the results on different meshes for the three configurations. Four meshes were tested, from 1000 elements per meter up to 100000. The conclusion that stands for the three configurations is, while the 1000 element mesh allows to get the correct order of magnitude of the velocity, it is not fine enough to catch all the intermediate states, at it can be seen on Configuration 1 velocity U 3 figure. When the mesh is refined, the intermediate states are well retrieved, the shocks are straight, and the gap between the the meshes gets reduced with each refinement. A noteworthy observation lays in the fact that we have shown two distinct velocities for Configuration 1, we can see that the profiles of U 2 and U 3 are quite different, and this is not surprising since the PDE model is built on the basis of distinct velocity fields, and the numerical scheme behind respects this setup. For Configuration 2 though, we did not mention in figure 19 which velocity was represented, because the three of them are superimposed. This will be analysed with more details in the next point (comparison between Configuration 1 and Configuration 2 ). For Configuration 3, only one velocity profile is shown, since we used the instantaneous relaxation algorithm, which ensures that at the end of each time step the velocities (and the pressures) are equal.

In the following comparisons, we consider the 50000 element per meter mesh. Owing to the results presented in figure 19, this mesh gives fair approximation of the model's solutions. For the velocities, we can notice that Configuration 1 's U 1 is not visible in the figure, this is due to the fact that the values taken by U 1 are significantly smaller than those of U 2 and U 3 , and therefore they could not be put in the same scale. Meanwhile, Configuration 2 's three velocity fields are overlapped at the scale of the figure. This means that dividing the relaxation time scales τ U by 10 N with N = 10 makes the velocities relax almost instantaneously, in a way such that the differences ∆U kl = U k -U l are almost null at the end of each time step. This is a significant point, its importance relies in the fact that it shows how these relaxation time scales drive the relaxation process's speed. The different intermediate values of N that were tested (and are not shown here) confirm this relaxation behaviour: the larger N is, the smaller ∆U kl relative velocities are at the end of each time step.

In figure 20's pressure profile, no difference between the pressures of Configuration 1 is visible. This suggests that this case's setup do not allow to bring to light the pressures differences, even if there are visible on the velocities level.

Finally, we turn the reader's attention to the fact that the change in the relaxation time scales did not affect only the velocity differences, it has also impacted the position of the shock waves. Indeed, at t = 1.2 ms while the right-going shock wave is observed around x = 2.75 for Configuration 1, it is located around x = 2.15 for Configuration 2. In figure 21 we give the comparison of Configuration 2 and Configuration 3 results at t = 1.6 ms. As previously indicated, for Configuration 2 there is no need to specify which velocity or pressure field is represented, because all the three are overlapped. For Configuration 3 this is an obvious result, considering the involved instantaneous relaxations. The notable point in figure 21 though is the fact that there is a persistent lag between the results of the two configurations. In other words, this means that the general relax-ation procedures (employed in Configuration 1 ) do not converge to the instantaneous relaxation algorithms used in Configuration 3, even if the relaxation time scales τ P and τ U are significantly small, such is the case of Configuration 2.

Appendix A7 -An analytical solution for the heat transfer subsystem

We recall the heat transfer sub-system (56):

       ∂ t α k = 0 ∂ t m k = 0 ∂ t (m k U k ) = 0 ∂ t (α k E k ) = 3 l=1,l =k K T kl (W)(T l -T k )
We focus particularly on the energy balances that have been written under the form (60):

∂ t T k = 3 l=1,l =k K T kl (W) m k C V,k (T l -T k )
We expand the three equations, and replace the last one by the total conservation of the internal energies, this gives:

     ∂ t T 1 = K T 12 (W) m 1 C V,1 (T 2 -T 1 ) + K T 13 (W) m 1 C V,1 (T 3 -T 1 ) ∂ t T 2 = K T 12 (W) m 2 C V,2 (T 1 -T 2 ) + K T 23 (W) m 2 C V,2 (T 3 -T 2 ) ∂ t (m 1 C V,1 T 1 + m 2 C V,2 T 2 + m 3 C V,3 T 3 ) = 0 (146)
We define:

E 0 = m 1 C V,1 T 1 (t = 0) + m 2 C V,2 T 2 (t = 0) + m 3 C V,3 T 3 (t = 0) (147)
We can therefore write:

T 3 (t) = 1 m 3 C V,3 [E 0 -m 1 C V,1 T 1 (t) -m 2 C V,2 T 2 (t)] (148)
We then replace in the first two equations of (146) to get: This allows to fully decouple the equations in (149), and their resolution becomes easier. Practically, we choose a set of time-dependent coefficients K T kl such that:

  Phase k = 1
m 1 C V,1 K T 23 = m 2 C V,2 K T 13 = m 3 C V,3 K T 12 = t + 1 (151)
Therefore, if we denote:

α = 1 m 1 C V,1 m 2 C V,2 + 1 m 1 C V,1 m 3 C V,3 + 1 m 2 C V,2 m 3 C V,3 β = E 0 m 1 C V,1 m 2 C V,2 m 3 C V,3 (152)
The system (149) becomes:

∂ t T 1 + α(t + 1)T 1 = β(t + 1) ∂ t T 2 + α(t + 1)T 2 = β(t + 1) (153)
The solution is then given by: T 1 (t) = β α + T 0 1 -β α exp(-α(t 2 /2 + t)) T 2 (t) = β α + T 0 2 -β α exp(-α(t 2 /2 + t))

T 3 (t) is given by (148). This determines completely the analytical solution. In practice, we considered the following initial dataset:

m 2 a 12 = e 13 m 1 -e 23 m 2 a 21 = e 12 m 1 -e 23 m 3 a 22 =

 21212211322 e 13 ( 1

Algorithm 3 . 3 . 2 . 1

 3321 In each cell Ω i , starting with W * i we compute W n+1 i by following the sequence: a) Initialize the pressure differences Y * i = Y(W * i ) and the pressure relaxation matrix defined in (52): B(W * i ); b) Compute Y n+1 i , the exact solution of the ODE:

  This allows to conclude that the equation g(x) = 1 admits a unique solution in D g . Algorithm 3.3.2.2 Again, in each cell Ω i , starting with W * i we compute W n+1 i by following the sequence: a) Initialize the pressure differences Y * i = Y(W * i ) and the pressure relaxation matrix defined in (52): B(W * i ); b) Compute Y n+1 i , the exact solution of the ODE:

Figure 1 :

 1 Figure 1: L 1 of the error for ρ k , p k , α k and u k wrt the mesh size h

Figure 2 :

 2 Figure 2: Spatial representation of variables at the final simulation time T = 0.8 ms

Figure 3 :

 3 Figure 3: Test case A3.1 convergence rates for the algorithms 3.3.1.1 (left) and 3.3.1.2 (right)

Figure 4 : 2 Figure 5 :

 425 Figure 4: Test case A3.2 convergence rates for the algorithms 3.3.1.1 (left) and 3.3.1.2 (right)

Figure 6 :

 6 Figure 6: Heat transfer verification test case: convergence rates on the temperature variables T k

Figure 7 :

 7 Figure 7: Mass transfer verification test case: convergence rates of the variables m 2 and m 3

Figure 8 :

 8 Figure 8: Shock-tube apparatus (bed hight 400mm)

Figure 9 :

 9 Figure 9: Shock wave in a single-phase medium: experimental data (blue) and numerical results with the energy equations (red)

Figure 10 :

 10 Figure 10: Shock wave through gas-liquid two-phase medium -Mean pressure in stations 2 & 3 -Comparison of experimental and numerical results

Figure 11 :

 11 Figure 11: The relative velocity u 2 -u 1 in the stations 2 and 3

Figure 12 :

 12 Figure 12: Phase 1 statistical fraction and interfactial area at three distinct instants

Figure 13 :

 13 Figure 13: Stations 2 and 3 pressure histories for different meshes

Figure 15 :

 15 Figure 15: Velocity at t = 0.7 ms on various meshes

Figure 16 :

 16 Figure 16: Pressure histories in stations S k , k = 1, .., 4

Figure 17 :

 17 Figure 17: Pressure (left) and velocity (right) at 5 time instants

em 2 e 12 ;•

 12 12 = 1.00 ; e 13 = m 3 Test case A3.2: time-dependent e kl In this case we choose the coefficients e kl such that: m 1 e 23 (w) = m 2 e 13 (w) = m 3 e 12 (w) = t + 1 (131)

Figure 19 :

 19 Figure 19: Velocities on different meshes for the three configurations at different time instants

  (ii) Comparison of Configuration 1 and Configuration 2

Figure 20

 20 Figure 20 presents the comparison of the pressure and the velocity profiles of Configuration 1 and Configuration 2 at the time instant t = 1.2 ms.

Figure 20 :

 20 Figure 20: Comparison of Configurations 1 and 2 -Velocity (left) and pressure (right) profiles at t = 1.2 ms

Figure 21 :

 21 Figure 21: Comparison of Configurations 2 and 3 -Velocity (left) and pressure (right) profiles at t = 1.6 ms

  2,3 (w int ) for i = 1, .., 8. After that, we consider the phase 1 right-going shock wave such that ρ R 1 /ρ int 1 = 1/2, and use the exact connections through shock waves to determine the state w R . The following table gives the adopted numerical values:Left state w L Intermediate state w int Right state w R

	α 1	3/10	3/5	3/5
	α 2	7/20	1/5	1/5
	α 3	7/20	1/5	1/5
	ρ 1	1/4	1/10	5/100
	ρ 2	0.14675324	0.14368748	0.14368748
	ρ 3	0.16666666	0.14535470	0.14535470
	p 1	1.00 × 10 4	1.0597086 × 10 4	3.8534858 × 10 3
	p 2	1.13 × 10 4	1.0901632 × 10 4	1.0901632 × 10 4
	p 3	1.20 × 10 4	9.6407273 × 10 3	9.6407273 × 10 3
	u 1	300.0	300.0	40.315571
	u 2	350.0	389.36692	389.36692
	u 3	400.0	560.65856	560.65856

Table 1 :

 1 Mass transfer verification test case: initialisation and EOS dataset

		00	0.00	0.00
	p k	7.0 × 10 5	7.0 × 10 5	4.0 × 10 5
	ρ k	998.1	5.0	20.0
	α k	0.3	0.5	0.2
	γ k	27.07	10.0	10.0
	Π k	8.06 10 7	0.00	0.00
	C V,k	10.58	700.0	700.0
	q k	0.00	-1836.1	-1836.1

Table 2 :

 2 Single-phase flow experimental, analytical and numerical results

	One-phase	Experimental Euler model (exact solution) Euler model (numerical simulation)
	flow	results	Barotropic	With energy	Barotropic	With energy
	P * (bar)	≈ 2.4	2.78	2.46	2.78	2.46
	P * * (bar)	≈ 5.0	6.85	5.42	6.85	5.42
	Celerity σ (m/s)	≈ 466		497		497

Table 3 :

 3 Shock tube initialization dataset

Table 5 :

 5 Velocity relaxation initialisation dataset

		00	150.0	200.0
	p k	7.10 5	7.10 5	4.10 5
	ρ k	10.00	1.00	1.50
	α k	0.30	0.50	0.20
	γ k	27.07	3.00	2.00
	Π k	8.06 10 7	0.00	0.00

Table 8 :

 8 Phase k = 2 Phase k = 3 Heat transfer initialisation dataset

	u k	2.00	3.00	7.00
	p k	7.10 3	1.10 4	4.10 3
	ρ k	5.20	2.70	0.30
	α k	0.30	0.35	0.35
	γ k	27.07	1.40	2.40
	Π k	8.06 10 7	0.00	0.00

This allows to obtain the limit: lim x→-x 0 + h(x) = +∞ and therefore the desired result.

Remark:

• The statistical fractions defined by (117) are in the admissible range [0, 1]. In fact, (117) ensures that α k ≥ 0, moreover, P satisfies (118), hence 0 ≤ α k ≤ 1.

• This instantaneous pressure relaxation algorithm can be defined for a N -phase flow where N > 3. The fractions α k and the pressure P are defined in the same way in the Stiffened Gas EOS case.

Appendix A2 -Threshold effect

In the two-phase framework, we consider the counterpart of (50), which represents the pressure relaxation sub-step:

We denote:

Using the same method of the paragraph 2.2.6, we prove that the pressure difference y 12 = p 1 -p 2 satisfies the equation:

In order to simplify the analysis of the equation, we denote:

We denote:

We notice that the differential equation verified by x 12 is:

By solving this equation, we have the general solution of (121):

The behaviour of y 12 (t) therefore depends on the initial condition y 0 12 = y 12 (0):

• If y 0 12 < 0: The solution y 12 (t) is defined on R + , and we can check that it is decreasing, which ensures the relaxation behaviour.

• If y 0 12 > 0 and y 0 12 < y * 12 where y * 12 is such that: 1 = y * 12 ∞ 0 b(s)e -µ(s) ds: y 12 (t) is defined on R + , and we can prove in the case a(t) = a 0 and b(t) = b 0 that the relaxation behaviour is ensured: lim t→∞ y 12 (t) = 0.

• If y 0 12 > y * 12 : In this case there exists a certain t 0 ∈ R * + such that:

The function y 12 is then not well defined on R + and the relaxation can't be ensured.

Appendix A3 -Two particular analytical solutions of the velocity relaxation system

In the velocity relaxation sub-step, we recall that the concerned PDE system is (45). In the pure 1D framework it can be rewritten as:

Appendix A4 -A particular analytical solution of the pressure relaxation system

We recall that the considered PDE system in the pressure relaxation sub-step is (50):

where the partial masses and momentums are constants:

Thus, the variables are: (α 2 , α 3 , α 1 p 1 , α 2 p 2 , α 3 p 3 ).

In order to build analytical solution for the system (135), we consider the case of the Ideal Gas EOS:

Property A4.1: Under the assumption (136), we have the following result:

where K 1 , K 2 and K 3 are positive integration constants:

The solution of (135) with assumptions (139) and (140) reads:

where A k = K k K 1 , α 0 2 and α 0 3 are the initial values.

Proof:

Under the assumptions (139) and (140), the system (138) becomes:

To integrate these ODEs we denote for k = 2, 3:

The equation verified by Ψ k is:

By solving this equation we obtain:

This allows to retrieve the result (141).

In Table 6 we give the initial dataset that was used for the pressure relaxation verification computations:

Appendix A5 -Initial data of the shock tube test cases of paragraph 4.3

In this appendix, we present the initial datasets that were considered in the shock tube numerical simulations.

Single-phase flow configuration:

Gas-solid and gas-liquid two-phase flow configurations:

97 < x < 3.37] and 10 -10 if not

The EOS were chosen such that:

• Phase 1: Stiffened Gas EOS (liquid water)

with:

• Phase 2: Perfect Gas EOS (air)

with γ 2 = 1.40 ;

• Phase 3: identical to phase 2.

Appendix A6 -Study of the pressure-velocity relaxation effects

In this appendix, we focus on the pressure-velocity relaxation effects embedded in the global PDE system [START_REF] Ambroso | Relaxation and numerical approximation of a two-fluid two-pressure diphasic model[END_REF]. In Section 3, different pressure-velocity relaxation algorithms were presented, and a two-phase flow test case was given in Section 4. Here we consider a threephase shock-tube configuration (see figure 18), where a shock-wave carried by the phase 2 travels along the low-pressure side of the tube. It crosses two distinct domains where phases 1 and 3 are present at different proportions.

With: α 1 = 0.01 for x ∈ [2.0, 2.2] and α 1 = 10 -10 elsewhere. For k = 3, α 3 = 0.01 for both x ∈ [2.0, 2.2] and x ∈ [3.08, 3.12], and α 3 = 10 -10 elsewhere. The chosen EOS are a mix of Stiffened Gas and Perfect Gas EOS, as shown below in Table 7. As presented in subsection 3.3, we have many simulation possibilities considering the different pressure and velocity relaxation algorithms, as well as the parameters that are not

Phase k γ