O Hurisse

L Quibel

Simulations of liquid-vapor water flows with non-condensable gases on the basis of a two-fluid model

Keywords: Two-phase flows, blowdown, out-of-equilibrium, mass transfer, two-fluid model, non-condensable gas

Two phase flows including liquid water, vapor and non-condensable gases are often encountered in industrial applications. In this paper, an extension of the classical Baer-Nunziato two-fluid model is presented in order to account for the non-condensable gases. In this model, the vapor and non-condensable gases form a miscible mixture of gases at temperature equilibrium and kinematic equilibrium. The complete model then allows to account for the full kinematic and thermodynamic disequilibrium between the liquid phase and the gas phase. Moreover, a new form for the masstransfer source-term has been retained. It leads to a more straightforward numerical integration. When focusing on the numerical scheme, a new algorithm is proposed here for the pressure relaxation effect. At last, some numerical experiments are performed in order to assess the behavior of the model with respect to the relaxation time-scales.

Introduction

Several industrial applications require to handle simulations where water phase change occurs in presence of non-condensable gases. This is the case when cryogenic flows, fuel combustion, cavitation and condensation phenomena are at stake. More specifically in the nuclear domain, numerous safety studies require to deal with pressurized hot water in presence of non-condensable gases. One may for instance quote two of them. First the steam-explosion phenomenon may occur during Reactivity Insertion Accident (RIA). In the latter, some very hot particles and hot noncondensable gases are expelled into hot pressurized water, transferring a great amount of energy to the fluid. The consequence is a rapid vaporization of the liquid into vapor which is accompanied by the generation of strong pressure waves that travel into the devices. The second situation of interest in this work is the Loss Of Coolant Accident (LOCA). In LOCA scenarii for Pressurized Water Reactor (PWR) power plants, high pressurized liquid water at 155 bars/300 • C vaporizes through a breach into a room that contains air at 1 bar/20 • C. As for steam-explosion, the estimation of the pressure waves that travel through the domain are of prime interest. Different numerical strategies have been proposed in order to simulate water phase change in presence of non-condensable gases, as for instance the Smooth Particle Hydrodynamics method [START_REF] Yang | Smoothed particle hydrodynamics method for evaporating multiphase flows[END_REF], the Lattice-Boltzamnn method [START_REF] Zheng | Single droplet condensation in presence of non-condensable gas by a multi-component multi-phase thermal lattice boltzmann model[END_REF], level-set methods [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF][START_REF] Ménard | Coupling level set/vof/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet[END_REF], Volume of Fluid methods [START_REF] Ménard | Coupling level set/vof/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet[END_REF][START_REF] Nikolopoulos | A numerical investigation of the evaporation process of a liquid droplet impinging onto a hot substrate[END_REF][START_REF] Strotos | Numerical investigation on the evaporation of droplets depositing on heated surfaces at low weber numbers[END_REF]. Nevertheless, we focus here on strategies based on finite volume techniques and on diffuse interface models. They seem appropriate to simulate complex scenarii involving strong shock waves traveling in large domain, as those described above. Many diffuse-interface two-phase flow models have been proposed in order to perform numerical simulations, see among many others [START_REF] Downar-Zapolski | The non-equilibrium relaxation model for one-dimensional flashing liquid flow[END_REF][START_REF] Bilicki | Relaxation models for wave phenomena in liquid-vapor bubble flow in channels[END_REF][START_REF] Faucher | Computation of flashing flows in variable cross-secglimm1998two,tion ducts[END_REF][START_REF] Bilicki | A thermodynamically consistent relaxation model for turbulent binary mixture undergoing phase transition[END_REF][START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF][START_REF] Lochon | Computation of fast depressurization of water using a two-fluid model: revisiting Bilicki modelling of mass transfer[END_REF][START_REF] Helluy | Assessment of numerical schemes for complex two-phase flows with real equations of state[END_REF][START_REF] Hurisse | Simulations of a simplified LOCA scenario with a non-equilibrium homogeneous model[END_REF][START_REF] Saurel | Modelling phase transition in metastable liquids: Application to cavitating and flashing flows[END_REF]. LOCA scenarii have in particular been studied for instance in [START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF][START_REF] Lochon | Computation of fast depressurization of water using a two-fluid model: revisiting Bilicki modelling of mass transfer[END_REF][START_REF] Helluy | Assessment of numerical schemes for complex two-phase flows with real equations of state[END_REF][START_REF] Hurisse | Simulations of a simplified LOCA scenario with a non-equilibrium homogeneous model[END_REF]. In all these references, the models do not account for the incondensable gases and the air has thus been replaced by hot vapor. Even if good agreement has been found with respect to experimental measurements, it was assumed that the liquid of the primary circuit did not contain air or incondensable gases. Yet, a mixture of air and water can have thermodynamical properties that are very different from those of the pure water or pure air, see [START_REF] Karplus | The velocity of sound in a liquid containing gas bubbles, armour research foundation[END_REF][START_REF] Feldman | The speed of sound in single-component two-phase fluids: theoretical and experimental[END_REF][START_REF] Kieffer | Sound speed in liquid-gas mixtures: Water-air and water-steam[END_REF] for instance when considering the speed of propagation of a pressure wave. At last, as shown for instance in [START_REF] Bartak | A study of the rapid depressurization of hot water and the dynamics of vapour bubble generation in superheated water[END_REF][START_REF] Helluy | Assessment of numerical schemes for complex two-phase flows with real equations of state[END_REF][START_REF] Hurisse | Simulations of a simplified LOCA scenario with a non-equilibrium homogeneous model[END_REF], accounting for out-of-equilibrium states seems mandatory for the LOCA simulations with high pressure and high temperature. In the present paper, a model dedicated to the simulation of flows with these features is used to simulate a LOCA scenario.

Since these scenarii involve strong shock waves and high speed flows, first-order conservationlaw models are promising. Unfortunately, there exists few of these mathematical models in the literature that allow to account for water/air flows including heat, mass and energy transfer and that are able to cope with out-of-equilibrium situations. This is particularly true if one focuses on those that are built on a rigorous thermodynamical background and that possess a minimal set of mathematical properties. Let us be more explicit on what the terms "thermodynamical background" and "minimal set of mathematical properties" contain. Obviously a model involving thermodynamical aspects should be closed using thermodynamical closures that are in agreement both with the Gibbs relation and with the second law of thermodynamics. Moreover, mixing miscibility (for vapor and incondensable gases) and immiscibility (between liquid and the "gas" phase) constraints has to be introduced rigorously in the model. It is an interesting point to note that the second law of thermodynamics is also of importance for the "minimal set of mathematical properties". Indeed, it may help to define an entropy inequality that permits to uniquely define the Rankine-Hugoniot jump relations, and thus to get a unique definition of the shocks. A second important mathematical property that is classically required for first-order conservation-law models is the hyperbolicity. The latter is directly linked to the notion of stability in time of the analytical solutions of the model.

Among the models available in the literature, one can for example quote the following references [START_REF] Mathis | A thermodynamically consistent model of a liquid-vapor fluid with a gas[END_REF][START_REF] Hérard | A three-phase flow model with two miscible phases[END_REF][START_REF] Hérard | A four-field three-phase flow model with both miscible and immiscible components[END_REF] or [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows. application to evaporation phenomena[END_REF][START_REF] Furfaro | Modeling droplet phase change in the presence of a multi-component gas mixture[END_REF][START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF][START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF][START_REF] Deng | Diffuse interface modelling of reactive multi-phase flows applied to a sub-critical cryogenic jet[END_REF]. The articles of the former set are based on the same thermodynamical closure assumptions and they meet the thermodynamical and mathematical properties mentioned in the previous paragraph even if they are based on slightly different sets of Partial Differential Equations (PDE). The miscible phases share the same volume which allows to recover the Dalton's law in a very natural manner, whatever the EOS (Equation Of States) for each phases are. In the following, this assumption is chosen which means that miscible phases share the same volume. The details of the thermodynamical closures are explained in Appendix A. We thus focus in this work on the references [START_REF] Mathis | A thermodynamically consistent model of a liquid-vapor fluid with a gas[END_REF][START_REF] Hérard | A three-phase flow model with two miscible phases[END_REF][START_REF] Hérard | A four-field three-phase flow model with both miscible and immiscible components[END_REF]. It should be noted that a different set of thermodynamical closures has been proposed in [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows. application to evaporation phenomena[END_REF][START_REF] Furfaro | Modeling droplet phase change in the presence of a multi-component gas mixture[END_REF][START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF][START_REF] Boivin | A thermodynamic closure for the simulation of multiphase reactive flows[END_REF][START_REF] Deng | Diffuse interface modelling of reactive multi-phase flows applied to a sub-critical cryogenic jet[END_REF] and it has been applied to several models, including two-fluid models and homogeneous models. In these works, it is assumed that each phase occupies its own volume, which is a classical choice for instance when dealing with immiscible phases. An explanation of the difference between these two sets of closures can be found in [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF].

In [START_REF] Hérard | A three-phase flow model with two miscible phases[END_REF][START_REF] Hérard | A four-field three-phase flow model with both miscible and immiscible components[END_REF], the multifluid approach has been chosen. These models are built as extensions of the classical two-phase flow model proposed in [START_REF] Baer | A two-phase mixture theory for the Deflagration-To-Detonation transition (DDT) in reactive granular materials[END_REF] or in many other works, see [START_REF] Powers | Theory of two-phase detonation-part I: modeling[END_REF][START_REF] Glimm | Two-pressure two-phase flow[END_REF][START_REF] Gavrilyuk | Mathematical and numerical modeling of two-phase compressible flows with micro-inertia[END_REF][START_REF] Coquel | Closure laws for a two-fluid two-pressure model[END_REF][START_REF] Drew | Theory of multicomponent fluids[END_REF] for instance. For multifluid models [START_REF] Hérard | A three-phase flow model with two miscible phases[END_REF][START_REF] Hérard | A four-field three-phase flow model with both miscible and immiscible components[END_REF], as in the Baer-Nunziato model [START_REF] Baer | A two-phase mixture theory for the Deflagration-To-Detonation transition (DDT) in reactive granular materials[END_REF], each phase possesses its own velocity field, so that even the vapor and the air may have different velocities. The return to the thermodynamical equilibrium is ensured by a "classical" form of source terms for which all the effects are independent, in the sense that each one is associated with its own time-scale. In the homogeneous model proposed in [START_REF] Mathis | A thermodynamically consistent model of a liquid-vapor fluid with a gas[END_REF], which extends the model proposed in [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF], it is assumed that all the phases share the same velocity field. On the contrary to the multifluid models cited above, the return to the thermodynamical equilibrium is ensured through source terms that are coupled and that are associated with a unique time-scale. This model has been used to perform numerical simulations of a simplified LOCA scenario in [START_REF] Helluy | Assessment of numerical schemes for complex two-phase flows with real equations of state[END_REF][START_REF] Hurisse | A homogeneous two-phase flow model with non-condensable gas using Noble-Able-Chemkin EOS[END_REF][START_REF] Quibel | Simulation of water-vapor two-phase flows with non-condensable gas[END_REF][START_REF] Helluy | Simulation of a liquid-vapour compressible flow by a Lattice Boltzmann Method[END_REF]. It should be noted that homogeneous models can also be used with independent source terms as in [START_REF] Ghazi | Modélisation d'écoulements compressibles avec transition de phase et prise en compte des états métastables[END_REF][START_REF] Ghazi | Vapour-liquid phase transition and metastability[END_REF] for instance. Moreover, many homogeneous approaches used for two-phase flow modeling are based upon the model proposed in [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF]. The latter is built by reducing the Baer-Nunziato model [START_REF] Baer | A two-phase mixture theory for the Deflagration-To-Detonation transition (DDT) in reactive granular materials[END_REF] to a five-equation model. Some applications require to compute a relative velocity between the liquid phase and the gas phase, for instance when dealing with the fragmentation of droplets for steam-explosion as in [START_REF] Boukili | Simulation and preliminary validation of a three-phase flow model with energy[END_REF][START_REF] Boukili | Simulation schemes of a non-miscible three-phase flow model for vapor explosion applications[END_REF]. For homogeneous models, the relative velocity is not computed and it could only be reconstructed using correlations. Besides, the multifluid approach seems more appealing for these situations since the relative velocity can be computed using the velocity field of each phase. Despite many accurate schemes are available in the literature for two-fluid two-phase flows [START_REF] Coquel | A positive and entropy-satisfying finite volume scheme for the Baer-Nunziato model[END_REF][START_REF] Ambroso | Relaxation and numerical approximation of a two-fluid two-pressure diphasic model[END_REF][START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF][START_REF] Lochon | HLLC-type riemann solver with approximated two-phase contact for the computation of the Baer-Nunziato two-fluid model[END_REF], there exists less possibilities when turning to a higher number of components [START_REF] Saleh | A relaxation scheme for a hyperbolic multiphase flow model[END_REF]. Moreover, it has been shown in [START_REF] Boukili | Simulation and preliminary validation of a three-phase flow model with energy[END_REF][START_REF] Boukili | Simulation schemes of a non-miscible three-phase flow model for vapor explosion applications[END_REF] that the numerical simulation of the relaxation source terms of multifluid models requires complex numerical schemes whereas some robust schemes are available for twofluid models [START_REF] Hérard | A fractional step method to compute a class of compressible gasliquid flows[END_REF][START_REF] Hérard | Computing two-fluid models of compressible water-vapour flows with mass transfer[END_REF][START_REF] Liu | Contribution to the verification and the validation of an unsteady two-phase flow model[END_REF][START_REF] Lochon | Modélisation et simulation d'écoulements transitoires eau-vapeur en approche bifluide[END_REF][START_REF] Gallouët | Numerical modeling of two-phase flows using the twofluid two-pressure approach[END_REF]. Considering these points, we choose as a modeling basis a two-fluid model with pressure relaxation as proposed in [START_REF] Baer | A two-phase mixture theory for the Deflagration-To-Detonation transition (DDT) in reactive granular materials[END_REF]. Hence, we may expect to inherit from all the numerical schemes that have been developed during the past 15 years for two-fluid models. The two-fluid model proposed in this paper is an extension of the Baer-Nunziato model [START_REF] Baer | A two-phase mixture theory for the Deflagration-To-Detonation transition (DDT) in reactive granular materials[END_REF] where the two phases represent: a liquid water phase, and a gas phase that corresponds to an homogeneous miscible mixture of water vapor and incondensable gas. In the gas phase, we assume that the temperature equilibrium always holds between the vapor and the incondensable gas. Hence, the classical set of PDE of the Baer-Nunziato model just needs to be supplemented by a sole balance equation for the mass of incondensable gas. Finally, the system of equations associated with the model proposed here corresponds to the system of equations of the Baer-Nunziato model supplemented by an advection equation for the mass of incondensable gas and with a mixture pressure law for the gas phase that depends on the mass of incondensable gas. The modeling of the thermodynamical behavior of the gas phase is then obtained following the classical thermodynamical principles used to obtain homogeneous models as in [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF][START_REF] Mathis | Theoretical and numerical study of phase transition flows[END_REF][START_REF] Jung | Numerical simulations of two-fluid flow on multicores accelerator[END_REF][START_REF] Quibel | Simulation of water-vapor two-phase flows with non-condensable gas[END_REF].

The paper is organized as follows. In section 2, the model is presented. The convective part of the model directly inherits its properties from the underlying Baer-Nunziato model and some properties of the relaxation source terms are exhibited. It should be mentioned that following [START_REF] Hurisse | Various choices of source terms for a class of two-fluid two-velocity models[END_REF], the source term for the mass transfer is not in a classical form. The overall numerical scheme proposed in section 3 is based on a splitting of the set of PDE into the convective part and the source terms. The simulation of the convective part is based on a VFRoe-ncv scheme [START_REF] Buffard | A sequel to a rough Godunov scheme: application to real gases[END_REF], whereas new numerical schemes are proposed in order to deal with the thermodynamical relaxation effects. At last, some numerical results of the simulation of a simplified LOCA scenario are proposed in section 4.

Extension of the Baer-Nunziato model to account for incondensable gases

System of equations and closure laws

Let us consider a mixture of liquid water, vapor and incondensable gas. Each of these phases are respectively denoted by a subscript l, v and a. The mixture of phases v and a will be referred in the following as the gas phase, and it will be denoted by a subscript g. The steam and the incondensable gas are assumed to be miscible, and to share the same velocity U g and the same temperature T g . We then consider a model for the liquid phase l and the gas phase g which is directly derived from the classical two-fluid Baer-Nunziato model [START_REF] Baer | A two-phase mixture theory for the Deflagration-To-Detonation transition (DDT) in reactive granular materials[END_REF] and which is inspired from [START_REF] Hérard | A three-phase flow model with two miscible phases[END_REF]:

                       ∂ t (α g ) + U g ∂ x (α g ) = Φ g , ∂ t (α a ρ a ) + ∂ x (α a ρ a U g ) = 0, ∂ t (α g ρ g ) + ∂ x (α g ρ g U g ) = -Γ l , ∂ t (α g ρ g U g ) + ∂ x α g (ρ g U 2 g + P g ) -P l ∂ x (α g ) = -S u l , ∂ t (α g ρ g E g ) + ∂ x (α g U g (ρ g E g + P g )) + P l ∂ t (α g ) = -S E l , ∂ t (α l ρ l ) + ∂ x (α l ρ l U l ) = Γ l , ∂ t (α l ρ l U l ) + ∂ x (α l (ρ l U 2 l + P l )) -P l ∂ x (α l ) = S u l , ∂ t (α l ρ l E l ) + ∂ x (α l U l (ρ l E l + P l )) + P l ∂ t (α l ) = S E l .
(1)

The system above has been written in a form that mimics the Baer-Nunziato model. Theses eight equations involve several quantities which need to be defined, and some closures have to be specified. We consider here the set of the variable:

W = (α g , ρ a , ρ g , U g , e g , ρ l , U l , e l ),
where α g is the gas fraction and it is related to α l through a statistical conservation property of the whole mixture α l + α g = 1 [START_REF] Hérard | Une approche bifluide statistique de modélisation des écoulements diphasiques à phases compressibles[END_REF][START_REF] Liu | Contribution to the verification and the validation of an unsteady two-phase flow model[END_REF][START_REF] Lochon | Modélisation et simulation d'écoulements transitoires eau-vapeur en approche bifluide[END_REF], U k are the velocities, ρ k the densities, e k = E k -U 2 k /2 are the specific internal energies and E k the total specific energies. Then, the pressure P g , P l and the sources terms have to be defined with respect to W . The pressures and the temperatures of the miscible mixture of the vapor and the incondensable gas, i.e. phase g, are obtained through a classical approach described in details in Appendix A. The main closure relations for the gas phase derived in Appendix A are:

α g = α a = α v , ρ g = ρ a + ρ v ,
ρ g e g = ρ a e a + ρ v e v , and P g = P a (ρ a , e a ) + P v (ρ v , e v ), [START_REF] Zheng | Single droplet condensation in presence of non-condensable gas by a multi-component multi-phase thermal lattice boltzmann model[END_REF] where the internal energies e a and e v are linked through the temperature equilibrium assumption:

T a (ρ a , e a ) = T v (ρ v , e v ).
Closures (2) are a direct consequence of the assumption that the vapor and the incondensable gas share the same volume. The gas pressure-law, i.e. last relation in [START_REF] Zheng | Single droplet condensation in presence of non-condensable gas by a multi-component multi-phase thermal lattice boltzmann model[END_REF], then obeys the classical Dalton's law (or Gibbs-Dalton's law). Moreover, it should be noted that the entropy of phase g depends on ρ a , ρ g and e g ; whereas the liquid entropy only depends on ρ l and e l . The source terms are defined below on the basis of the second law of thermodynamics.

Let us focus here on the definition of the closure laws for the source terms. For that purpose, we follow a classical way. Since the source terms Φ g , Γ l , S E l and S u l account for interfacial transfer between the phases, they should agree with the second law of thermodynamics. In order to proceed, we define the mixture entropy s as the sum of the phasic entropies:

ρs = m g s g (ρ a , ρ g , e g ) + m l s l (ρ l , e l ), (3)
where ρ = m g + m l and m k = α k ρ k . The equation for the entropy of the gas phase is given in Appendix A by equation (A.10). Thanks to the closure relation for P g (A.11), it reads:

∂ t (α g ρ g s g ) + ∂ x (α g ρ g U g s g ) = Φg Tg (P g -P l ) + Γ l µv Tv -1 Tg S E l -U g S u l + U 2 g 2 Γ l , (4)
and the equation for the evolution of the liquid entropy can be written from system (1):

∂ t (α l ρ l s l ) + ∂ x (α l ρ l U l s l ) = -Γ l µ l T l + 1 T l S E l -U l S u l + U 2 l 2 Γ l . (5)
Hence, the equation of evolution of the mixture specific entropy s is obtained by summing equations ( 4) and ( 5):

∂ t (ρs) + ∂ x (α l ρ l U l s l + α g ρ g U g s g ) = Φg Tg (P g -P l ) + Γ l µv-U 2 l /2 Tv -µ l -U 2 v /2 T l + S E l 1 T l -1 Tg + S u l Ug Tg -U l T l . (6)
The source term S u l corresponds to the exchange of momentum between the phase g and the phase l and it is composed of a part which is directly due to the drag force and a part due to the mass transfer. In the same way, the term S E l acts on the total energies which combine internal energies and kinetic energies. Hence, S E l gathers the energy exchange due to the pure heat exchange and the amount of heat exchange due to the momentum exchange. We choose here to write S u l and S E l in a classical manner [START_REF] Coquel | Closure laws for a two-fluid two-pressure model[END_REF][START_REF] Coquel | Modelling compressible multiphase flows[END_REF][START_REF] Hérard | Une approche bifluide statistique de modélisation des écoulements diphasiques à phases compressibles[END_REF] which allows to separate the different contributions:

S u l = D l + Γ l U l + U g 2 , (7)
S E l = Γ l U l U g 2 + D l U l + U g 2 + Ψ l , (8)
where D l stands for the drag force of the phase g on the phase l, and Ψ l stands for the pure heat transfer. The equation for the mixture entropy given by (6) then reads:

∂ t (ρs) + ∂ x (α l ρ l U l s l + α g ρ g U g s g ) = Φg Tg (P g -P l ) + Γ l µv Tv -µ l T l + Ψ l 1 T l -1 Tg + D l 1 Tg + 1 T l (U l -U g ). (9)
Following for instance [START_REF] Hurisse | Various choices of source terms for a class of two-fluid two-velocity models[END_REF], it can be proved the following result when excluding single phase situations: the mixture entropy s is concave with respect to the conservative thermodynamical variables α k , m k and m k e k . Therefore, to ensure the increase of the entropy s through the source terms in equation ( 9), one possible choice is to close the terms Φ g , Γ l , Ψ l and D l such that all the terms on the right hand side are positive. As proposed in [START_REF] Hurisse | Various choices of source terms for a class of two-fluid two-velocity models[END_REF], we thus choose:

Γ l = m l -m l λ m , (10)
Φ g = α g (1 -α g ) (P g -P l ) λ P , (11)
Ψ l = m g T g m l T l m g T 2 g + m l T 2 l (T g -T l ) λ T , (12)
D l = m g m l m g + m l (U g -U l ) λ U , (13)
where λ P , λ m , λ T and λ U are positive but may depend on W . The choices [START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF], [START_REF] Lochon | Computation of fast depressurization of water using a two-fluid model: revisiting Bilicki modelling of mass transfer[END_REF] and ( 13) are classical, whereas the mass transfer source term [START_REF] Bilicki | A thermodynamically consistent relaxation model for turbulent binary mixture undergoing phase transition[END_REF] is not classical [START_REF] Hurisse | Various choices of source terms for a class of two-fluid two-velocity models[END_REF]. The mass transfer source term [START_REF] Bilicki | A thermodynamically consistent relaxation model for turbulent binary mixture undergoing phase transition[END_REF] involves an equilibrium liquid mass, denoted by m l , which is defined below. A more classical form for the mass transfer could be:

Γ l = K m µ v T v - µ l T l ,
with K m > 0. The non-classical form [START_REF] Bilicki | A thermodynamically consistent relaxation model for turbulent binary mixture undergoing phase transition[END_REF] retained here simply arises from the remark that:

µ v T v - µ l T l = ∂s ∂m l |αg,ma,mgeg,m l e l ,m l +mg , (14)
and that s m : m l → s(α g , m a , e a , m 0 -m a -m l , m g e g , m l , m l e l ) is strictly concave for m l ∈]0, m 0m a [ with m 0 = m l + m g + m a . Then there exists a unique maximum for s m on ]0, m 0 -m a [, this maximum is reached at m l = m l :

∀m l ∈]0, m 0 -m a [, s m (m l ) ≤ s m (m l ). ( 15
)
Moreover, due to the concavity of s m we also have the classical result for m l ∈]0, m 0 -m a [:

∀m l ∈]0, m 0 -m a [, s m (m l ) ≤ s m (m l ) + ∂s m (m l ) ∂m l (m l -m l ), (16)
which, by setting m l = m l in (16) and using [START_REF] Saurel | Modelling phase transition in metastable liquids: Application to cavitating and flashing flows[END_REF], leads to:

∀m l ∈]0, m 0 -m a [, (m l -m l )
∂s ∂m l |αg,ma,mgeg,m l e l ,m l +mg ≥ 0.

Hence, the second term on the right hand side of equation ( 9), which can be written:

Γ l µ v T v - µ l T l = Γ l ∂s ∂m l |αg,ma,mgeg,m l e l ,m l +mg
is positive for the form of the mass transfer proposed here by equation [START_REF] Bilicki | A thermodynamically consistent relaxation model for turbulent binary mixture undergoing phase transition[END_REF]. This proves that choice [START_REF] Bilicki | A thermodynamically consistent relaxation model for turbulent binary mixture undergoing phase transition[END_REF] for the mass transfer term implies an increase of the entropy of the mixture s, and that it agrees with the second law of thermodynamics. The reader is referred to [START_REF] Hurisse | Various choices of source terms for a class of two-fluid two-velocity models[END_REF] for a more detailed discussion. It will be shown in section 3.2 that this non-classical form for the mass transfer term helps in defining robust numerical schemes for the mass transfer.

Some properties of the model

The model described in section 2.1 is very close to the Baer-Nunziato model and it thus shares a lot of its properties. We propose here a quick overview of some of these properties.

In order to avoid complex computations, we introduce in this section some specific EOS for the phases k ∈ {a, l, v}. We consider the classical Stiffened Gas EOS whose entropy reads:

s k (τ k , e k ) = C V,k ln (e k -Q k -Π k τ k )τ (γ k -1) k + s 0 k , (17)
where

C V,k , Q k , Π k , γ k , s 0 k are constant parameters.
We introduce here the mass fraction y a of phase a as:

y a = ρ a /ρ g , (18)
which represents the mass fraction of incondensable into the sole gas phase. Before going further in the study of the properties of the model, let us remark that in this case the thermodynamic behavior of phase g follows a Stiffened Gas EOS for a fixed y a . Indeed, the specific entropy of the gas phase is:

s g (y a , τ g , e g ) = C V,g (y a ) ln (e g -Q g (y a ) -Π g (y a )τ g )τ (γg(ya)-1) g + s 0 g (y a ), (19)
where the parameters C V,g , Q g , Π g , γ g and s 0 g only depend on y a and on the parameters of each phase a and v. This feature has been developed in detail in Appendix A and it has already been used for instance in [START_REF] Bachmann | A simple model for cavitation with noncondensable gases[END_REF].

Some properties of the convective part of the model

The eigenstructure is studied with the non-conservative set of variables:

Y = (α g , y a , U g , P g , s g , U l , P l , s l ).
When neglecting the source terms, system (1) can be written for regular solutions:

∂ t (Y ) + A(Y )∂ x (Y ) = 0,
with the 8 × 8 convection matrix:

A(Y ) =             U g 0 0 0 0 0 0 0 0 U g 0 0 0 0 0 0 (P g -P l )/m g 0 U g τ g 0 0 0 0 0 0 c 2 g /τ g U g 0 0 0 0 0 0 0 0 U g 0 0 0 0 0 0 0 0 U l τ l 0 ρ l c 2 l (U g -U l )/α l 0 0 0 0 c 2 l /τ l U l 0 0 0 0 0 0 0 0 U l             . ( 20
)
The sound speed c k are defined through the relations:

c 2 l (τ l , e l ) = τ 2 l P l ∂P l ∂e l |τ l - ∂P l ∂τ l |e l , (21)
and

c 2 g (y a , τ g , e g ) = τ 2 g P g ∂P g ∂e g |τg,ya - ∂P g ∂τ g |eg,ya . (22)
The strict concavity of the phasic entropies s k with respect to (τ k , e k ) and the positivity of the temperatures T k guarantee that c 2 k > 0 for k ∈ {a, l, v}. Consequently, thanks to relation (A. [START_REF] Bartak | A study of the rapid depressurization of hot water and the dynamics of vapour bubble generation in superheated water[END_REF]), we also have: c 2 g > 0. This implies that the eigenvalues λ k of the matrix A(Y ) belong to R and that the set of the associated right eigenvectors r k spans R 8 . They read:

λ 1 = U g , r 1
= (1, 0, 0, 0, 0, 0, 0, 0, 0) ; λ 2 = U g , r 2 = (0, 1, 0, 0, 0, 0, 0, 0, 0) ; λ 3 = U g -c g , r 3 = (0, 0, 1, -c g /τ g , 0, 0, 0, 0, 0) ; λ 4 = U g , r 4 = (0, 0, 0, 0, 1, 0, 0, 0, 0) ; λ 5 = U g + c g , r 5 = (0, 0, 1, c g /τ g , 0, 0, 0, 0, 0) ; λ 6 = U l -c l , r 6 = (0, 0, 0, 0, 0, 0, 1, -c l /τ l , 0) ; λ 7 = U l , r 7 = (0, 0, 0, 0, 0, 0, 0, 0, 1) ;

λ 8 = U l + c l , r 8 = (0, 0, 0, 0, 0, 0, 1, c l /τ l , 0) .
System ( 1) is thus hyperbolic provided that resonance does not occur: it occurs when for |U l -U g | = c l . From a practical point of view, when introducing the drag force, the situations we are interesting in always lead to small relative velocities with respect to the sound speeds. The resonance phenomena is therefore very unlikely to appear. The fields associated with the eigenvalues λ k for k ∈ {1, 2, 4, 7} are linearly degenerate since ∇ Y λ k .r k = 0 for k ∈ {1, 2, 4, 7} ; and the fields associated with the remaining eigenvalues λ k are genuinely non-linear. It should be mentioned here that the field λ 1 , associated to α g , is linearly degenerate thanks to the choice of U g for transporting α g [START_REF] Coquel | Closure laws for a two-fluid two-pressure model[END_REF]. According to [START_REF] Coquel | Closure laws for a two-fluid two-pressure model[END_REF], two other velocities transport for α g allow to recover such a property:

U l and (m l U l + m v U v )/(m l + m v )
. The choice between these three closures does not make consensus, but it has been noticed in [START_REF] Lochon | Modélisation et simulation d'écoulements transitoires eau-vapeur en approche bifluide[END_REF] that they give very close numerical results when the drag force is accounted for.

When examining the set of eigenvectors, it can be noted that the fraction α g remains constant when crossing the waves associated with λ k , k ∈ {3, 5, 6, 7, 8}; and that conversely the velocity U g and the pressure P l respectively do not change in the waves associated with λ k , k ∈ {3, 5} and λ k , k ∈ {6, 8}. Then, when crossing a wave α g and P l do not jump simultaneously, neither do α g and the product P l U g . As a consequence, the non-conservative products that appear in the momentum and energy equations of system (1) allow to define unique jump relations for the shocks in the genuinely non-linear waves. The uniqueness of the definition of shocks is due to the linearly degenerate nature of the field associated to λ 1 , see [START_REF] Coquel | Closure laws for a two-fluid two-pressure model[END_REF].

As shown in [START_REF] Coquel | Two properties of two-velocity two-pressure models for two-phase flows[END_REF] for the Baer-Nunziato model, system (1) is symmetrizable in the variable Y following the same argument. Indeed, the difference between the convection matrix A(Y ), given by [START_REF] Mathis | A thermodynamically consistent model of a liquid-vapor fluid with a gas[END_REF], and the convection matrix of the corresponding Baer-Nunziato model only consists in the components associated with y a . In matrix A(Y ) this corresponds to the second line and second column, which only contain a non zero term on the diagonal of A(Y ).

Some properties of the source terms

From now the properties that have been highlighted referred to the convective part of the model. We now describe some properties that concern the source terms (10)-( 12). We will focus in particular on positivity properties for the fraction α g , for the partial masses m l and m g , and for the temperatures T l and T g . Since some EOS, as the Stiffened Gas EOS, allow the internal energy to become negative, we retain here the temperature and the density as the thermodynamical quantities to be positive. We then consider uniform solutions of system (1) which are regular in time and bounded. It is also assumed in the following that initial conditions are such that:

α g (t = 0) ∈]0, 1[, m g (t = 0) > 0, m l (t = 0) > 0, T g (t = 0) > 0 and T l (t = 0) > 0.
The source terms Φ g , Γ l and Ψ l are studied individually.

Pressure relaxation Φ g .

The Ordinary Differential Equation (ODE) system associated with the term Φ g , which is associated with the pressure relaxation effect, reads:

       ∂ t (α g ) = Φ g , ∂ t (m g E g ) + P l ∂ t (α g ) = 0, ∂ t (m l E l ) + P l ∂ t (α l ) = 0, ∂ t (m a ) = ∂ t (m g ) = ∂ t (m g U g ) = ∂ t (m l ) = ∂ t (m l U l ) = 0. ( 23
)
From first equation of [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows. application to evaporation phenomena[END_REF] we can obtain that:

∂ t (α g (1 -α g )) = (1 -2α g )Φ g = α g (1 -α g ) (1 -2α g )(P g -P l ) λ P .
As a consequence, provided that (1 -2α g )(P g -P l )/λ P remains bounded, α g (t) lies in ]0, 1[. We straightforwardly get from the mass equations that: m g (t) > 0 and m l (t) > 0. In order to maintain positive temperatures, the terms λ P has to be modified. Indeed, by setting:

1 λ P = T g T l λ P (T g + T l ) ,
the energy equations, i.e. second and third equations of ( 23), ensure that the temperatures are positive. Since we are dealing with regular solutions, these equations can be written:

m k ∂e k ∂T k |ρ k ∂ t (T k ) + ∂ρ k ∂T k |e k ∂ t (ρ k ) = -P l ∂ t (α k ) ,
and by using the mass equations, it yields:

m k ∂e k ∂T k |ρ k ∂ t (T k ) = ρ 2 k ∂ρ k ∂T k |e k -P l ∂ t (α k ) .
By combining the equations for l and g, we get:

∂ t (T l T g ) = T l T g (A g T l + A l T g ), (24)
where:

A k = α g (1 -α g ) ρ 2 k ∂ρ k ∂T k |e k -P l (P g -P l ) λ P m k ∂e k ∂T k |ρ k (T g + T l )
.

Hence, provided that (A g T l + A l T g ) remains bounded, the temperatures are positive: T l (t) > 0 and T g (t) > 0.

Mass transfer Γ l .

The ODE system associated with the mass transfer, i.e. the source term Γ l , is:

                   ∂ t (α g ) = ∂ t (m a ) = 0, ∂ t (m l ) = Γ l , ∂ t (m l U l ) = Γ l U l +Ug 2 , ∂ t (m l E l ) = Γ l U l Ug 2 , ∂ t (m l ) + ∂ t (m g ) = 0, ∂ t (m l U l ) + ∂ t (m g U g ) = 0, ∂ t (m l E l ) + ∂ t (m g E g ) = 0, (25)
where we have:

Γ l = m l -m l λ m .
The liquid mass at equilibrium m l belongs to ]0, m l (t = 0) + m g (t = 0)[ (see section 2), hence, second equation of system [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF] implies that m l (t) also belongs to ]0, m l (t = 0) + m g (t = 0)[. Thanks to the conservation of the total mass, i.e. fifth equation of ( 25), we deduce that m g (t) ∈ ]0, m l (t = 0) + m g (t = 0)[. Moreover, first equation of system [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF] straightforwardly gives that the fraction α g belongs to ]0, 1[. Let us now turn to the temperatures. It should first be noted that we get from (25) that:

∂ t (m k e k ) = 0, k ∈ {l, g},
so that when combined with the mass equations, we obtain:

m k ∂ t (T k ) = ρ k ∂T k ∂ρ k |e k -e k ∂T k ∂e k |ρ k ∂ t (m k ) .
Therefore, setting:

1 λ m = T g T l λ m (T g + T l ) 2 ,
and combining the temperature equation for each phase leads to:

∂ t (T l T g ) = T l T g (m l B l -m g B g ) . (26)
with:

B k = T k (m l -m l ) m g m l λ m (T g + T l ) 2 ρ k ∂T k ∂ρ k |e k -e k ∂T k ∂e k |ρ k .
As a consequence, provided that the term (m l B l -m g B g ) remains bounded, the temperatures are positive: T l (t) > 0 and T g (t) > 0.

Pure heat transfer Ψ l .

The ODE system associated with the term Ψ l is:

       ∂ t (α g ) = ∂ t (m a ) = ∂ t (m l ) = ∂ t (m g ) = 0, ∂ t (U l ) = ∂ t (U g ) = 0, ∂ t (m l e l ) = Ψ l , ∂ t (m g e g ) = -Ψ l . (27)
It can be straightforwardly obtained that m k (t) > 0 and α g (t) ∈]0, 1[. As for the source terms Ψ g and Γ l , we write the temperature equations, which read:

m l ∂e l ∂T l |ρ l ∂ t (T l ) = T g T l m g m l m g T 2 g + m l T 2 l (T g -T l ) λ T = -m g ∂e g ∂T g |ρg ∂ t (T g ) ,
hence the evolution of the product of the pressure is:

∂ t (T l T l ) = T g T l (T g -T l ) λ T (m g T 2 g + m l T 2 l ) m g T g ∂T l ∂e l |ρ l -m l T l ∂T g ∂e g |ρg . ( 28
)
As a consequence, provided that the terms in the parenthesis on the right hand side of ( 28) is bounded, the temperatures remain positive: T l (t) > 0 and T g (t) > 0.

Numerical schemes

The overall numerical scheme is built on a fractional step approach using a Lie-Trotter splitting. First, the convective part of the model is accounted for using a classical first-order finite-volume approach. Numerical fluxes between the different cells are computed here thanks to the VFRoencv scheme [START_REF] Buffard | A sequel to a rough Godunov scheme: application to real gases[END_REF] using the variables Y already defined in section 2.2.1. This approximate Godunov solver is depicted in detail in [START_REF] Buffard | A sequel to a rough Godunov scheme: application to real gases[END_REF]. In order to assess the implementation of this scheme, a verification test case is proposed on the basis of a Riemann problem in Appendix C.

In the second step the source terms ( 11)-( 13) are solved sequentially. In the following, the numerical scheme associated with the drag force [START_REF] Helluy | Assessment of numerical schemes for complex two-phase flows with real equations of state[END_REF] and for the heat exchange [START_REF] Lochon | Computation of fast depressurization of water using a two-fluid model: revisiting Bilicki modelling of mass transfer[END_REF] are classical and they have been presented in several references [START_REF] Hérard | A fractional step method to compute a class of compressible gasliquid flows[END_REF][START_REF] Hérard | Computing two-fluid models of compressible water-vapour flows with mass transfer[END_REF][START_REF] Liu | Contribution to the verification and the validation of an unsteady two-phase flow model[END_REF][START_REF] Lochon | Modélisation et simulation d'écoulements transitoires eau-vapeur en approche bifluide[END_REF], we thus focus in this section on the thermodynamical source terms [START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF] and [START_REF] Bilicki | A thermodynamically consistent relaxation model for turbulent binary mixture undergoing phase transition[END_REF]. From now, we assume that the time-step ∆t is given by the scheme used to solve the convective part, and that the approximated solution at time t n is known. Since the time step is imposed by the convective scheme, we need to built implicit schemes for the source terms in order to maintain a good accuracy and stability of numerical approximations. We detail below how to compute the approximated solution at time t n+1 = t n + ∆t through the fractional steps associated with the source terms (11) (section 3.1) and (10) (section 3.2).

We insist on the fact that the numerical schemes are not restricted to the EOS used in section 2.2. Moreover, it should be stressed out that the splitting of the system associated with the fractional step algorithm complies with the entropy inequality without any constraint. Furthermore, in the schemes dedicated to each source term, conservation properties have been imposed to the numerical approximated solutions when they are required for the associated continuous solutions (i.e. total mass conservation, total momentum conservation and total energy conservation).

At last, we retain here the source terms [START_REF] Bilicki | A thermodynamically consistent relaxation model for turbulent binary mixture undergoing phase transition[END_REF], [START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF], [START_REF] Lochon | Computation of fast depressurization of water using a two-fluid model: revisiting Bilicki modelling of mass transfer[END_REF], [START_REF] Helluy | Assessment of numerical schemes for complex two-phase flows with real equations of state[END_REF] without applying the normalization by the temperature as proposed in 2.2.2. In order to get similar magnitudes in the choice of the time-scale, the following normalizations are applied : λ P = P 0 λ P and λ T = λ T /C V,0 , where P 0 is a reference pressure and C V,0 is a reference heat capacity. With these choices, λ P and λ T are time-scales, as λ m and λ U .

An algorithm for the pressure relaxation term

For this step we need to compute approximated solutions of system of ODE [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows. application to evaporation phenomena[END_REF]. Algorithms have been proposed for instance in [START_REF] Gallouët | Numerical modeling of two-phase flows using the twofluid two-pressure approach[END_REF][START_REF] Hérard | A fractional step method to compute a class of compressible gasliquid flows[END_REF][START_REF] Hérard | Computing two-fluid models of compressible water-vapour flows with mass transfer[END_REF][START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF] for the pressure relaxation. The scheme proposed here relies on the same kind of techniques but it takes advantage of the liquid entropy equation to limit the number of terms that have to be taken explicit. In fact we propose an implicit Euler scheme for all the quantities. The sole relaxation parameter λ P is explicit. System [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows. application to evaporation phenomena[END_REF] straightforwardly provides the solutions:

m a (t n+1 ) = m a (t n ), m g (t n+1 ) = m g (t n ), m l (t n+1 ) = m l (t n ), ( 29
) U g (t n+1 ) = U g (t n ), U l (t n+1 ) = U l (t n ). ( 30
)
It should be noted that thanks to (5), we get with Γ l = D l = Ψ l = 0 that the liquid entropy is constant for system [START_REF] Métayer | Dynamic relaxation processes in compressible multiphase flows. application to evaporation phenomena[END_REF]. Hence we have:

s l (t n+1 ) = s l (t n ).
The first three equations of system ( 23) can be written in an equivalent manner:

∂ t (α g ) = α g (1 -α g ) (P g -P l ) λ P , (31)
∂ t (s l ) = 0, ( 32
) ∂ t (m g e g ) + ∂ t (m l e l ) = 0. (33)
Then, by integrating equations ( 32) and ( 33) between t n and t > t n , we get:

s l (t) = s l (t n ), ( 34
) m g (t n )(e g (t) -e g (t n )) + m l (t n )(e l (t) -e l (t n )) = 0. ( 35
)
When considering the thermodynamical plane (ρ l , s l ) for the liquid phase, we get from [START_REF] Gavrilyuk | Mathematical and numerical modeling of two-phase compressible flows with micro-inertia[END_REF] and [START_REF] Hurisse | A homogeneous two-phase flow model with non-condensable gas using Noble-Able-Chemkin EOS[END_REF]:

∂ t (α g ) = α g (1 -α g ) λ P P g α g m g (t n ) , e g (t) -P l m l (t n ) (1 -α g ) , s l (t n ) , (36)
e g (t) = e g (t n ) - m l (t n ) m g (t n ) e l m l (t n ) (1 -α g ) , s l (t n ) -e l (t n ) . (37)
By substituting e g (t) given by equation [START_REF] Helluy | Simulation of a liquid-vapour compressible flow by a Lattice Boltzmann Method[END_REF] into equation [START_REF] Quibel | Simulation of water-vapor two-phase flows with non-condensable gas[END_REF], we obtain an equation with a sole unknown: α g (t). Note that λ P also only depends on α g thanks to relations ( 29), ( 30), [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF] and [START_REF] Helluy | Simulation of a liquid-vapour compressible flow by a Lattice Boltzmann Method[END_REF]. We have chosen here to approximate λ P in an explicit manner: λ P ∼ λ P (α g (t = 0)).

Even with this simplification, solving explicitly equation [START_REF] Quibel | Simulation of water-vapor two-phase flows with non-condensable gas[END_REF] in α g can be difficult. Therefore, we compute an approximated solution α * g for α g at time t n+1 with the implicit Euler scheme:

   α * g -αg(t n ) ∆t = α * g (1-α * g ) λ P (t n ) P g α * g mg(t n ) , e * g -P l m l (t n ) (1-α * g ) , s l (t n ) e * g = e g (t n ) -m l (t n ) mg(t n ) e l m l (t n ) (1-α * g ) , s l (t n ) -e l (t n ) . (38)
Obviously, ( 38) is non-linear even for very simple equations of states. From a practical point of view, dichotomy algorithm or quasi-Newton methods are used to find α * g . Once the new fraction is computed, α n+1 g = α * g , all the other quantities are updated explicitly using relations [START_REF] Powers | Theory of two-phase detonation-part I: modeling[END_REF] for the masses, [START_REF] Glimm | Two-pressure two-phase flow[END_REF] for the velocities, [START_REF] Ghazi | Modélisation d'écoulements compressibles avec transition de phase et prise en compte des états métastables[END_REF] for the energies and the gas fraction:

α n+1 g = α * g , m n+1 a = m a (t n ), m n+1 g = m g (t n ), m n+1 l = m l (t n ), U n+1 g = U g (t n ), U n+1 l = U l (t n ), e n+1 l = e l m l (t n ) (1 -α * g ) , s l (t n ) , e n+1 g = e g (t n ) - m l (t n ) m g (t n ) e n+1 l -e l (t n ) .
Properties of the algorithm.

It can be proved that (38) possesses a unique solution for α * g on ]0, 1[ provided that the pressure gap in equation ( 38) remains bounded. The details of the proof may be found in Appendix B. By construction, the algorithm ensures the conservation of the total mass, the total momentum, and the total energy. The partial masses are constant for this algorithm and they are thus positive.

Accounting for the mass transfer

The system of equations that allows to account for the mass transfer arises from the closures ( 7) and ( 8) with D u = Ψ l = 0. The equations for the total energies of system (25) can be written in terms of the internal energies:

∂ t (m l e l ) = 0, ( 39
) ∂ t (m g e g ) + ∂ t (m l e l ) = 0, ( 40
)
which leads to the integration over [t n , t]:

m l e l (t) = m l e l (t n ), (41)
m g e g (t) = m g e g (t n ). ( 42
)
Moreover, from the conservation of the total mass, we obtain:

m l (t) + m v (t) = m l (t n ) + m v (t n ) = m lv (t n ). ( 43
)
The equilibrium liquid mass m l is defined as the liquid mass that realizes the maximum of the entropy s m defined in section 2. If we exclude single-phase flow situations, the maximum of s m is reached when its derivatives vanish. Hence, according to [START_REF] Hurisse | Simulations of a simplified LOCA scenario with a non-equilibrium homogeneous model[END_REF] and using ( 41)-( 43), m l is the unique solution of:

µ l T l α l (t n ) m l , m l e l (t n ) m l - µ v T v α v (t n ) m lv (t n ) -m l , m g e g (t n ) m lv (t n ) -m l = 0. ( 44
)
From a practical point of view, the equation above can be numerically solved by the mean of a dichotomy algorithm or quasi-Newton methods. An important point to be quoted is that m l only depends on the quantities at time t n and it is thus constant when considering system [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF]. The liquid mass equation issued from system (25) is then:

∂ t (m l (t)) = m l (t n ) -m l (t) λ m (t) . ( 45
)
If the time-scale λ m does not depend of the partial masses m v and m l , it is also constant and equation ( 45) can be straightforwardly integrated on [t n , t]. In the general case, when the form of λ m does not permit an exact integration, an approximated solution for system [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF] is obtained by computing the exact solution of system:

∂ t (m l (t)) = m l (t n ) -m l (t) λ m (t n ) . ( 46
)
which corresponds to system [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF] where the time-scale λ m has been fixed to its value at time t n . We thus get:

m l (t) = m l (t n ) + (m l (t n ) -m l (t n )) 1 -exp -(t-t n )/λm(t n ) ,
and, hence we set:

m n+1 l = m l (t n + ∆t) = m l (t n ) + (m l (t n ) -m l (t n )) 1 -exp -∆t/λm(t n ) . (47)
Once m n+1 l is computed, m v , m g e g and m l e l are updated according to relations ( 41)-( 43):

α n+1 g = α g (t n ), m n+1 a = m a (t n ), m n+1 v = m l (t n ) + m v (t n ) -m n+1 l , (m l e l ) n+1 = m l e l (t n ), (m g e g ) n+1 = m g e g (t n ).
It remains to update the velocities. Since the mass equations have been solved, the term Γ l in the momentum equation of system (25) is known. Let us set Γ * l ∆t = m n+1 l -m l (t n ). The subsystem for updating the velocities is then:

m l ∂ t (U l ) = Γ * l U g -U l 2 , (48)
m g ∂ t (U g ) = -Γ * l U g -U l 2 . ( 49
)
We then approximate the factor m k on the left hand side by a constant value mk which can be: the value at time t n , or the value at time t (which arises from the integration of the mass equation), or a barycentric combination of these two values. Since mk is constant, we now have to integrate the equations:

∂ t (U l ) = Γ * l ml U g -U l 2 , ( 50
) ∂ t (U g ) = - Γ * l mg U g -U l 2 . (51)
Equations ( 50)-( 51) can be exactly integrated by considering the relative velocity U l -U g and the momentum conservation. Indeed the equation for the relative velocity reads:

∂ t (U l -U g ) = -Γ * l 1 ml + 1 ml U l -U g 2 , (52)
and it gives the solution:

U l (t) -U g (t) = (U l (t n ) -U g (t n )) exp -(t -t n )Γ * l 1 ml + 1 ml .
Then, when combined with the equation of conservation of the total momentum:

m l (t)U l (t) + m g (t)U g (t) = m l (t n )U l (t n ) + m g (t n )U g (t n ),
one obtains a 2 × 2 linear system in (U l (t), U g (t)) allowing to update the velocities. It should be noted that the relative velocity may increase or decrease depending on the sign of Γ * l .

Properties of the algorithm.

Thanks to its definition [START_REF] Ambroso | Relaxation and numerical approximation of a two-fluid two-pressure diphasic model[END_REF], m l (t n ) belongs to ]0, m v (t n ) + m l (t n )[. Then, since the update for the liquid mass ( 47) is a barycentric formula, it implies that m n+1 l ∈]0, m v (t n ) + m l (t n )[. Thanks to the total mass conservation, we get that m n+1 v also belongs to ]0, m v (t n ) + m l (t n )[. The fraction α g and the energies m k e k are constant. Moreover, once the term Γ * l has been computed, the energies m k e k and the velocities are updated using exact integrations. Hence these updates are in agreement with the total energy equation of ( 25), and thus they agree with the total energy conservation for the mixture. At last, momentum conservation for the mixture is enforced by construction.

Simulation of the blowdown of a pressurized pipe containing hot water

We are interested here in the blowdown experiments reported in [START_REF] Riegel | Contribution à l'étude de la décompression d'une capacité en régime diphasique[END_REF]. These experiments are representative of the Loss Of Coolant Accident (LOCA). LOCA scenarii are classical for the safety studies in the nuclear domain. They represent the situation where a breach occurs on the primary circuit of a Pressurized Water Reactor (PWR). In these scenarii, the coolant in the primary circuit contains water at 300 • C which is pressurized at 150 bars. Hence, when a breach occurs in this circuit a vaporization front is propagated into the circuit, and a jet expels through the breach the steam created into the circuit. Initial condition for such a situation corresponds to single phases inside (water) and outside (air) the circuit. When using a two-velocity model, as the Bear-Nunziato model for instance, one has to initialize the computation with a mixture of water and steam inside and outside the circuit. This is due to the set of equations and to the variables that are not defined for pure single phase situations. Then two drawbacks can be noted. First, the outside of the circuit is filled with hot steam (with a small amount of liquid water) which may not behaves exactly as air. Secondly, in the pipes one has to set initial conditions such that the amount of vapor is sufficient to avoid its complete condensation before the arrival of the depressurization waves, see [START_REF] Lochon | Computation of fast depressurization of water using a two-fluid model: revisiting Bilicki modelling of mass transfer[END_REF] for instance. Unfortunately, it is known that a small amount of gas in the liquid can modify the speed of propagation of the wave in the mixture, see [START_REF] Liu | Contribution to the verification and the validation of an unsteady two-phase flow model[END_REF][START_REF] Hurisse | Numerical simulations of steady and unsteady two-phase flows using a homogeneous model[END_REF] or the results of section 4.1. By using the present model with a mixture of air, liquid and vapor in the pipes, equilibrium states with very small amount of gas can be used as initial conditions. Hence mass transfer is avoided before the arrival of the depressurization waves without modifying their speed of propagation.

The test case investigated here is associated with the experimental facility SUPERCANON [START_REF] Riegel | Contribution à l'étude de la décompression d'une capacité en régime diphasique[END_REF]. It was set up to measure the sudden depressurization of heated water from 150 bars to 1 bar, which is representative of a LOCA in the primary circuit of a PWR. Figure 1 closed with a cap. The water is heated and when the operating conditions are reached (in pressure and temperature), the cap is released by a system based on an exploding cordon. This release is assumed to be almost instantaneous with respect to the fluid phenomena. The pressure in the tube is measured at six different locations P 1,..,6 , and the vapor fraction is measured at the point P t (see figure 1). Three different initial temperatures have been tested for the liquid water in the tube: 280 • C, 300 • C and 320 • C, which respectively correspond to: the temperature at the inlet of the core, the mean temperature in the core and the temperature at the outlet of the core. Moreover, different breach diameters have been used at the outlet of the tube, but we only focus here on the case with a fully opened tube. We also restrict ourselves to the initial temperature of 300 • C.

The scenario of the experiment is the following [START_REF] Edwards | Studies of phenomena connected with the depressurization of water reactors[END_REF]. When the cap is released, a "saturation" wave travels from the cap location to the end of the tube. Due to this wave, the temperature in the pipe remains almost constant and the pressure drops to the saturation pressure at the initial temperature: P = P sat (573.15 K) = 86 10 5 P a. Then, a vaporisation front travels into the pipe and the vapor fraction starts to increase. The vaporisation front is a two-phase phenomenon and it travels much slower than the "saturation" wave which occurs in pure liquid water. Through the vaporisation front, the vapor fraction increases and both the pressure and the temperature drop again. The sketch of the time evolution of the pressure at points P 1 is plotted in figure 2 for two different experimental runs based on the "same conditions". The blue curve represents the typical pressure profile at P 1 that is obtained when enforcing the thermodynamical equilibrium in the simulations, which means that all the relaxation time scales in the source terms defined in section 2.1 are set to zero. For this experimental setting, the numerical results obtained by considering the thermodynamical equilibrium seems not in agreement with the experimental data. Indeed, for this experiment, the two phases are likely to be out of the thermodynamical equilibrium [START_REF] Bartak | A study of the rapid depressurization of hot water and the dynamics of vapour bubble generation in superheated water[END_REF]. As a consequence, the response of the model should be sensitive to the choice of the time-scale arising in the relaxation source terms, as shown for instance in [START_REF] Helluy | Assessment of numerical schemes for complex two-phase flows with real equations of state[END_REF] (when considering another model) or in [START_REF] Liu | Contribution to the verification and the validation of an unsteady two-phase flow model[END_REF][START_REF] Lochon | Modélisation et simulation d'écoulements transitoires eau-vapeur en approche bifluide[END_REF]. liquid fraction close to 1 and a small amount of vapor and air ; whereas the room-domain contains a small quantity of liquid. The left boundary condition is a wall, while the right boundary condition is an outlet. We are interested in the pressure values at x = 0.5 m, which is near the wall at the end of the pipe. The thermodynamical behavior of the liquid is modeled by a Stiffened Gas EOS, while air and vapor are modeled by a perfect gas law. Their parameters are gathered in table 1. Several simulations are performed using different times scales. Several simulations are performed in order to investigate the influence on the results of the different parameters of the model. We first study the influence of the gas/liquid fraction used as initial condition in the pipe. Then several relaxation times are considered for the different source terms. Since the model contains four relaxation source terms, performing a complete matrix of test cases leads to a huge amount of computations. We then restrict ourselves by changing one time-scale while keeping the others fixed. Obviously, this choice is arbitrary and it cannot be representative all the possible configurations, but it allows to highlight some influences of the time-scales on the results. In these simulations, a particular focus is set on the pressure undershoot that is observed for the experimental data with the squares in figure 2. This undershoot is associated with a non-equilibrium situation and it strongly depends on the relaxation time-scales as already shown by [START_REF] Helluy | Assessment of numerical schemes for complex two-phase flows with real equations of state[END_REF][START_REF] Hurisse | Simulations of a simplified LOCA scenario with a non-equilibrium homogeneous model[END_REF][START_REF] Quibel | Simulation of water-vapor two-phase flows with non-condensable gas[END_REF][START_REF] Liu | Contribution to the verification and the validation of an unsteady two-phase flow model[END_REF][START_REF] Lochon | Computation of fast depressurization of water using a two-fluid model: revisiting Bilicki modelling of mass transfer[END_REF][START_REF] Hurisse | Numerical simulations of steady and unsteady two-phase flows using a homogeneous model[END_REF]. Moreover, it is known that with stiffened gas EOS, the vaporization front travels too fast and that more accurate EOS are required in order to catch this front accurately, see [START_REF] Helluy | Assessment of numerical schemes for complex two-phase flows with real equations of state[END_REF][START_REF] Hurisse | Simulations of a simplified LOCA scenario with a non-equilibrium homogeneous model[END_REF][START_REF] Quibel | Simulation of water-vapor two-phase flows with non-condensable gas[END_REF].

k = l k = v k = a C V,k (J/K/
Remark. From a practical point of view, instantaneous relaxation (i.e. when a time scale is set to zero) is achieved by setting the associated numerical time-scale to 10 -35 s. Considering that for the test cases the time steps are much larger than 10 -35 s, this allows to recover a numerical relaxation up to the round-off error.

Influence of the gas fraction in the pipe

We first study the results of simulation for different values of the initial gas fraction in the pipe. We thus choose the following liquid fractions: α l = {0.99995, 0.9995, 0.995, 0.95}. The gas phase is then composed of the mixture of air and vapor ensuring the thermodynamical equilibrium with a mixture pressure of 150 bars and a temperature of 573.15 K. The thermodynamical equilibrium is then enforced during all the simulations by choosing λ P = λ m = λ T = λ U = 0.

The pressure at point P 1 , which is at a distance of 0.5 m of the end of the left wall, is plotted with respect to the time in figure 3. It can first be noticed that the gas fraction has a great influence on the results at a given mesh size. Indeed, the less gas there is in the liquid, the more the waves present stiff profiles. The smoothing of the stiff waves clearly arises from the higher compressibility of the liquid/gas mixture when the gas fraction increases. Moreover, it clearly appears on figure 3 that the first depressurization wave travels faster when the gas fraction is low, which is in accordance with [START_REF] Liu | Contribution to the verification and the validation of an unsteady two-phase flow model[END_REF][START_REF] Hurisse | Numerical simulations of steady and unsteady two-phase flows using a homogeneous model[END_REF]. When compared with the experimental results, it seems that the SUPERCANON experiments were performed with a mixture that contained a very low fraction of gas. So, in the following, we focus on the case with α l = 0.99995.

The different results for the three different meshes (see the bottom right plot in figure 3) clearly show that the convergence of the approximated solution (at the scale of the plots) require a smaller mesh size when the liquid fraction is important. This is obviously due to the stiffer wave profiles arising for low gas fractions. In the following section, we choose the mesh with 2000 cells for the different computations.

Influence of the mass transfer

In this section the liquid fraction α l in the pipe is chosen equal to 0.99995 and the influence of the mass transfer time-scale λ m is investigated. The temperature, pressure and velocity equilibrium are imposed by choosing λ P = λ T = λ U = 0. The computations are then performed with a timescale λ m = {10 -4 s, 10 -5 s, 10 -6 s, 10 -7 s} and the results are reported in figure 4. The pressure undershoot is very sensitive to the time-scale for the mass transfer. For λ m < 10 -5 s the pressure profile is very close to the results for λ m = 0 and for λ m > 10 -4 s the pressure undershoot is too important with respect to the experimental data. The range of values for λ m seems thus very restricted. From now, λ m is chosen equal to 10 -5 s which appears in figure 4 to give a value of the pressure undershoot close to the one observed in the experiments.

Influence of the temperature relaxation

We set here α l = 0.99995 and λ m = 10 -5 s. The pressure and velocity equilibrium are imposed by choosing λ P = λ U = 0. Several values for the temperature relaxation time-scale are used: λ T = {10 -4 s, 10 -5 s, 10 -6 s, 10 -7 s}. The pressure at point P 1 is plotted in figure 5 for the different values of λ T . As for the mass transfer relaxation time-scale, the time-scale for pressure relaxation has a great influence on the pressure undershoot. For λ T < 10 -5 s the pressure profile is very close to the results for λ T = 0 and for λ T > 10 -4 s the pressure undershoot is too important with respect to the experimental data. The time-scale retained for the next simulations is then equal to 10 -5 s.

Influence of the pressure relaxation

We set here α l = 0.99995 and λ m = λ T = 10 -5 s. The velocity equilibrium are imposed by choosing λ U = 0. Several values for the temperature relaxation time-scale are used: λ P = {10 -5 s, 10 -6 s, 10 -7 s}. The results are very sensitive with respect to λ P . The time-scale λ P = 10 -4 s leads to a divergence of the simulations because of a loss of hyperbolicty. The latter is due to a negative temperature because of the parameters Q k and Π k . Since the numerical scheme is based on the variable e k and ρ k (in order to compute easily the correct shocks), it may happen that an approximated value (ρ k , e k ) leads to a negative temperature, see also Appendix A. This drawback is associated with the Stiffened Gas EOS and our choice of a convection scheme based on the variablese k and ρ k .

Realistic model for the relaxation time scales

In the previous sections we have only considered constant and uniform time-scales. Some previous results clearly show that for this experiment, more complex laws for the time-scales seem to allow to provide more realistic results, see for instance [START_REF] Lochon | Computation of fast depressurization of water using a two-fluid model: revisiting Bilicki modelling of mass transfer[END_REF][START_REF] Helluy | Assessment of numerical schemes for complex two-phase flows with real equations of state[END_REF][START_REF] Hurisse | Simulations of a simplified LOCA scenario with a non-equilibrium homogeneous model[END_REF]. In [START_REF] Lochon | Computation of fast depressurization of water using a two-fluid model: revisiting Bilicki modelling of mass transfer[END_REF], a relaxation time-scale has been proposed for simulating the SUPERCANON experiment with a liquid-vapor Baer-Nunziato model. It has been adapted from the time-scale proposed in [START_REF] Downar-Zapolski | The non-equilibrium relaxation model for one-dimensional flashing liquid flow[END_REF] for the simulation of the flashing of superheated liquid flow in a nozzle. It should indeed be noted that in [START_REF] Downar-Zapolski | The non-equilibrium relaxation model for one-dimensional flashing liquid flow[END_REF] the model was different. The latter was a classical homogeneous model enforcing the pressure equilibrium and assuming that the vapor is always at saturation. We propose here to modify the parameters of the time-scale law in order to recover at least qualitatively the behavior observed on the experimental data. We thus set λ U = λ P = λ T = 0 and the law for the mass transfer time-scale is: , otherwise. [START_REF] Mathis | Theoretical and numerical study of phase transition flows[END_REF] In figure 7 the simulation results are compared with the experiment data : pressure at point P 1 and volume fraction of vapor at point P t are compared while using either λ m = 0 (equilibrium simulation) or modified correlation [START_REF] Mathis | Theoretical and numerical study of phase transition flows[END_REF]. It clearly appears that, in comparison with the equilibrium simulation, the vaporization is delayed when using the law given by [START_REF] Mathis | Theoretical and numerical study of phase transition flows[END_REF]. With the latter, the pressure undershoot after the first arrival of the first depressurization can be observed on the numerical results, even if it qualitatively does not perfectly match with the pressure undershoot of the experimental data. In the same way, the arrival of the vaporization front arises too early with respect to the experimental data. That point can be at least partially explained by the use of Stiffened Gas laws. Indeed, it has been shown in [START_REF] Hurisse | Simulations of a simplified LOCA scenario with a non-equilibrium homogeneous model[END_REF][START_REF] Quibel | Simulation of water-vapor two-phase flows with non-condensable gas[END_REF] that the use of more accurate phasic equations of states highly improves qualitatively the results. In particular the Stiffened Gas EOS seems to lead to a vaporization front that travels too fast towards the end of the tube.

λ m =    6.

Conclusion

The model studied in the present paper allows to deal with incondensable/steam/water flows. The main additional assumptions with respect to the model proposed in [START_REF] Hérard | A three-phase flow model with two miscible phases[END_REF] is that the two miscible phases, i.e. the steam and the incondensable gas, are at thermal equilibrium and at kinematic equilibrium. The resulting set of equations is then very close to the classical Baer-Nunziato model. This enables to inherit from all the work that has been done in the past decade to build numerical schemes for the convection part of this class of model.

Considering the interfacial exchanges, a mass transfer source term as proposed in [START_REF] Hurisse | Various choices of source terms for a class of two-fluid two-velocity models[END_REF] has been tested. Moreover, a new robust and efficient numerical scheme for the pressure relaxation effect has been introduced. This scheme can be applied with any EOS (provided that the latter is in agreement with the result of uniqueness of Appendix B) and it is strongly related to the form of the exact solution associated with the pressure relaxation subsystem.

In [START_REF] Lochon | Computation of fast depressurization of water using a two-fluid model: revisiting Bilicki modelling of mass transfer[END_REF], the time-scale for the mass transfer and the temperature relaxation are assumed to be linked. In our numerical tests, we did not consider such an assumption. It could be interesting to test such a link by using other source term models proposed in [START_REF] Hurisse | Various choices of source terms for a class of two-fluid two-velocity models[END_REF]. Indeed, in the latter source terms involving a unique time-scale have been proposed for the thermodynamical relaxation effects (mass transfer, pressure relaxation and temperature relaxation). This proposition of source term is clearly inspired from the homogeneous model as studied for instance in [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF][START_REF] Mathis | Theoretical and numerical study of phase transition flows[END_REF][START_REF] Jung | Numerical simulations of two-fluid flow on multicores accelerator[END_REF][START_REF] Hurisse | Numerical simulations of steady and unsteady two-phase flows using a homogeneous model[END_REF][START_REF] Quibel | Simulation of water-vapor two-phase flows with non-condensable gas[END_REF] among many others. Moreover, in [START_REF] Hurisse | A homogeneous two-phase flow model with non-condensable gas using Noble-Able-Chemkin EOS[END_REF][START_REF] Quibel | Simulation of water-vapor two-phase flows with non-condensable gas[END_REF] such an homogeneous model for steam-water flows including incondensable gases has been tested together with complex EOS. The results seem satisfactory and these EOS could be applied to the present model.

The specific entropies are supposed to be in C(R 2 + ) and strictly concave with respect to (τ k , e k ) on (R + ) 2 . Moreover, we assume that:

∂s k ∂e k |τ k > 0.
The phasic temperature T k and the phasic pressure P k of each phase is then defined thanks to the phasic entropy s k and to the Gibbs relation:

T k ds k = de k + P k dτ k , which yields: 1 T k (τ k , e k ) = ∂s k ∂e k |τ k and P k (τ k , e k ) T k (τ k , e k ) = ∂s k ∂τ k |e k .
Equation (A.3) straightforwardly gives ρ v with respect to W ; whereas equations (A.4) and (A.5) form a 2×2 non-linear subsytem that allows to get e a and e v from W . Once solved, all the densities ρ k and the specific internal energies e k , k ∈ {a, l, v}, are known with respect to W , and thus are known: the phasic specific entropies s k , the phasic temperatures T k and pressures P k , k ∈ {a, l, v}. Subsystem (A.4) and (A.

F(e a ) = T a (ρ a , e a ) -T v ρ g -ρ a , ρ g e g -ρ a e a ρ g -ρ a ,
where e a,m (ρ a ) is the minimal energy of phase a above which T a ≤ 0 and e a,M (ρ a ) is the minimal energy of phase a above which T v ≤ 0. Solving equation (A.6) is then equivalent to find e a such that F(e a ) = 0. The latter has a unique solution when the following sufficient conditions are fulfilled:

∂T k ∂e k |τ k > 0, k = {a, v}, (A.7) lim ea→ea,m(ρa) T a (ρ a , e a ) -T v ρ g -ρ a , ρ g e g -ρ a e a ρ g -ρ a < 0, (A.8) lim ea→e a,M (ρa)
T a (ρ a , e a ) -T v ρ g -ρ a , ρ g e g -ρ a e a ρ g -ρ a > 0. (A.9) Condition (A.7) ensures that e a → F(e a ) is a strictly increasing function. Since F is continuous, conditions (A.8) and (A.9) then guarantees that F(e a ) = 0 has a unique solution on ]e a,m (ρ a ), e a,M (ρ a )[. It remains to define the pressure P g and the temperature T g for the gas phase.

For that purpose, we focus on the entropy of phase g.

The entropy of the gas phase s g is defined as the sum of the entropies of phases a and v by (3):

ρ g s g = ρ a s a (ρ a , e a ) + ρ v s v (ρ v , e v ),
where ρ g and e g have been defined above by equations (A.1) and (A.2). When using system (1), one can get the equation of time and space evolution for entropy s g :

∂ t (α g ρ g s g ) + ∂ x (α g ρ g U g s g ) + αg Tg (P g -P a -P v )∂ x (U g ) = Φg Tg (P a + P v -P l ) + Γ l µv Tv -1 Tg S E l -U g S u l + U 2 g 2 Γ l , (A.10)
where the temperature equilibrium assumption T g = T a = T v holds. Then, for the entropy to be conserved along the streamlines of phase g (when no source terms are accounted for), we have to set P g -P a -P v = 0, which corresponds to the classical Dalton's law:

P g = P a (τ a , e a ) + P v (τ v , e v ). (A.11)
The temperature T g is generally not explicitly known with respect to W . Indeed, it arises from the solution of non-linear equation (A.6), but it only depends on W .

with γ g (y a )C V,g (y a ) = (y a γ a C V,a + (1 -y a )γ v C V,v ). Then, when substituting in equation (A. [START_REF] Karplus | The velocity of sound in a liquid containing gas bubbles, armour research foundation[END_REF]) the temperature T g by (A.15), we obtain:

P g = (γ g (y a ) -1) e g -Q g (y a ) τ g -γ g (y a )Π g (y a ), (A.17)
We finally turn to the entropy of phase g, and after some calculus, entropy s g reads:

s g (y a , τ g , e g ) = C V,g (y a ) ln (e g -Q g (y a ) -Π g (y a )τ g )τ (γg(ya)-1) g + s 0 g (y a ). (A.18)

For a fixed mass fraction y a , the mixture g behaves as a Stiffened Gas. We have indeed with (A. The sound speed in the phase g is defined through relation [START_REF] Hérard | A four-field three-phase flow model with both miscible and immiscible components[END_REF] and in each pure phase the sound speed is given by definition [START_REF] Hérard | A three-phase flow model with two miscible phases[END_REF]. With the Stiffened Gas EOS, the sound speed in each pure phase k ∈ {a, l, v} is given by equation (A.13). Then, thanks to definition [START_REF] Hérard | A four-field three-phase flow model with both miscible and immiscible components[END_REF] and to the temperature equilibrium assumption, one gets for phase g:

c 2 g (y a , τ g , e g ) γ g (y a ) = y a c 2 a γ a + (1 -y a ) c 2 v γ v . (A.19)
As a consequence, when the phasic sound speed c a and c v are in R + , then c g is also in R + .

Sufficient conditions (A.7)-(A.9) have been introduced in order to ensure that the energies e a and e v are defined in a unique manner through the closures (A.3)-(A.5). In the specific case of the Stiffened Gas EOS, these conditions are superfluous. Indeed temperature equilibrium can be solved explicitly by using z a (A.14). Depending on the choice of the different EOS parameters, z a may not belong to [0, 1]. This is a discrepancy which can be associated with the Stiffened Gas EOS. In particular, this EOS authorizes the energy to be non-positive if the parameter Q k is negative.

Appendix B. Uniqueness of the solution for the pressure relaxation algorithm

The aim of this section is to prove that equation ( 38) has a unique solution on ]0, 1[ provided that the EOS fulfill some assumptions. We recall that [START_REF] Ghazi | Modélisation d'écoulements compressibles avec transition de phase et prise en compte des états métastables[END_REF] is the key equation for the algorithm of pressure relaxation described in section 3.1. Once it is solved, all the other quantities are explicitly computed. Hence, if equation ( 38) has a unique solution on ]0, 1[, the whole algorithm possesses a unique solution for the whole conservative variable.

The first equation of (38) can be written:

α * g -α g (t n ) = α * g (1 -α * g )A n (α * g ), (B.1)
where the liquid EOS is written in the (τ l , e l )-plane:

A n (α , s l (t n ) -e l (t n ) .

In section (3.1), the liquid EOS is written in the (ρ l , e l )-plane, the use of the (τ l , e l )-plane here is equivalent.

First, when A n, * = 0, the solution of (B.1) is obviously unique. So, let us assume that A n, * = 0. Let us define on α * g ∈]0, 1[ the function α * g → F n (α * g ) as:

F n (α * g ) = α * g -α g (t n ) α * g (1 -α * g ) -A n (α * g ). (B.
2)

The derivative of F with respect to α * g is then: we get that α * g → F n (α * g ) is strictly increasing on ] ∈]0, 1[. It should be noted that for Stiffened Gas EOS, these inequalities are equivalent to: P g + γ g (y a )Π(y a ) ≥ 0; P l ≥ 0. Moreover, if the pressure laws for the phases g and l are such that: ∀e g , lim τg→0 + P g (τ g , e g ) = +∞, and ∀s l , lim τ l →0 + P l (τ l , s l ) = +∞, (B.5)

dF n (α * g ) dα * g = (α * g ) 2 -
we get that: lim through the five following Riemann invariants:

I Ug = {U g , s l , α l ρ l (U l -U g ), α l P l + α g P g + α l ρ l (U l -U g ) 2 , e l + P l /rho l + (U l -U g ) 2 /2}. (C.2)

Computing the different states is classical and it is not recalled here. The error between the exact solution Z(t, x) at time t and the approximate solution obtained with the scheme of section 3 is measured according to the relative L 1 -norm:

err = N i=1 |Z(t, x i ) -Z n i | N i=1 |Z(t, x i )|
, where x i denotes the center of the cell i = {1, .., N }, and Z n i the approximate value in cell i at iteration n. In the following, the CFL number is equal to 1/2. The computational domain is [0, 1], and it is discretized by N uniform cells. The logarithm of the error err at time t = 2 10 -4 s is plotted in figure C.9 with respect to the logarithm of the mesh size 1/N for N from 100 to 100000. Figure C.10 shows the exact solution at time t = 2 10 -4 s and the approximate solutions for meshes with 400 and 2000 cells. Due to the stiffness of the liquid EOS, strong spurious bumps can be observed at the location of the U l -c l shock and at the location of the U l + c l "ghost wave". The amplitude of these patterns tends to decrease when the mesh is refined. The left column corresponds to variables for the phase g, whereas the right column corresponds to the variables for the phase l. On the last line, the plot on the left corresponds to ya and the plot on the right to α g . In the U l -C l shock wave, velocity an pressure undergo a jump whose amplitude is small with respect to the spurious bumps (see table C.3). At the scale of the plot they can not be seen.

Figure 2 :

 2 Figure 2: Sketch of the time evolution of the pressure at the point P 1 . The black curves represent experimental measurements for two different runs (for the same conditions). The blue curve represents the pressure profile simulated when considering the thermodynamical equilibrium.

Figure 3 :

 3 Figure 3: Comparisons for the pressure at point P 1 for different liquid fractions and different mesh sizes.

Figure 4 :

 4 Figure 4: Comparisons for the pressure at point P 1 for α l = 0.99995 and for different time-scales for the mass transfer λ m . The figure on the right represents a zoom on the pressure undershoot.

Figure 5 :

 5 Figure 5: Comparisons for the pressure at point P 1 for α l = 0.99995 and for different time-scales for the temperature relaxation λ T . The figure on the right represents a zoom on the pressure undershoot.

Figure 6 :

 6 Figure 6: Comparisons for the pressure at point P 1 for α l = 0.99995 and for different time-scales for the pressure relaxation λ P . The figure on the right represents a zoom on the pressure undershoot.

  51 10 -4 (α v ) -0.257 Psat(T )-P Pc-Psat(T ) -1.76 , if P < 10 6 P a;3.84 10 -10 (α v ) -0.84 Psat(T )-P Psat(T ) -5

Figure 7 :

 7 Figure7: Comparisons for the pressure at point P 1 and for the void fraction at point P t , considering the experiment (black crosses +), the equilibrium simulation for λ m = 0 (red line) and the simulation using the correlation (53) (blue line). Initial condition in the pipe corresponds to a volume fraction of liquid α l = 0.99995.

5 )

 5 can be reduced to a sole equation by introducing ρ v from equation (A.3) and e v from (A.4) into equation (A.5). It yields a non-linear equation for finding e a : T a (ρ a , e a ) = T v ρ g -ρ a , ρ g e g -ρ a e a ρ g -ρ a . (A.6) Let us denote e a → F(e a ) the function defined on ]e a,m (ρ a ), e a,M (ρ a )[ by:

g

  →0 + A n (α * g ) = +∞, and limα * g →1 -A n (α * g ) = -∞.Conditions (B.5) thus imply that limα * g →0 + F n (α * g ) = -∞,and lim α * g →1 -F n (α * g ) = +∞. Assuming conditions (B.4), (B.5) and that α n g ∈]0, 1[, the theorem of the intermediate values applied to function α * g → F n (α * g ) on ]0, 1[ then gives the results that (B.1) admits a unique solution on ]0, 1[. It should be noted that conditions (B.5) are fulfilled by the Stiffened Gas EOS.

Figure C. 8 :

 8 Figure C.8: Sketch of the Riemann problem.

Figure C. 9 :Figure C. 10 :

 910 Figure C.9: Convergence curves for the different variables.

  is a sketch of the facility. A tube (100 mm of inner diameter and 4.389 m long) is filled with water and Figure 1: Sketch of the SUPERCANON configuration. The pressure is measured at the points P 1 , ..., P 6 and the vapor fraction is measured at the point P t .

		tube							cap
				P t
	00000000000000000000000000000 00000000000000000000000000000 00000000000000000000000000000 00000000000000000000000000000 00000000000000000000000000000 00000000000000000000000000000 11111111111111111111111111111 11111111111111111111111111111 11111111111111111111111111111 11111111111111111111111111111 11111111111111111111111111111 11111111111111111111111111111
		P 1	P 2		P 3	P 4	P 5 P 6	X (m)
	0	0,5	1,356	1,502	1,648	2,6925	3,885	4,259	4,389

Table 1 :

 1 Parameters for the EOS for the three components.

	kg) 1.452904592629688 10 3 4.441148752333071 10 3	7.18 10 2
	γ k (-)	1.614924811807376	1.085507894797296	1.4000231
	Π k (P a)	3.563521398523755 10 8	0	0
	Q k (J)	0	0	0
	s 0 k (J/K/kg)	0	-4.769786773517021 10 4	0

  g = e g (t n ) -m l (t n ) m g (t n ) e l (1 -α * g ) m l (t n )

* g ) = ∆t λ P (t n ) P g α * g m g (t n ) , e * g -P l (1 -α * g ) m l (t n )

, s l (t n ) , and e *

  2α g (t n )α * g + α g (t n ) 2 (α * g (1 -α * g )) 2For α g (t n ) ∈]0, 1[, the first term on the left hand side of equation B.3 is strictly positive. Then, if one assumes that the EOS of phases l and g fulfill the inequalities:

	-	∆t λ P (t n )	1 m n g	∂P g ∂τ g	n, * |eg	+	1 m n g	∂P g ∂e g	n, * |τg	∂e l ∂τ l	n, * |s l	+	1 m n l	∂τ l ∂P l	|s l n, *	.	(B.3)
	∂P g ∂τ g |eg	≤ 0;	∂P g ∂e g |τg	≤ 0;		∂e l ∂τ l |s l	≥ 0;	∂P l ∂τ l |s l	≤ 0;		(B.4)

  Table C.3 gives the different states that compose the Riemann problem retained here. They have been obtained with the EOS parameters of table C.2. The U l -c l and U g + c g shock waves and the contact wave travel respectively with a speed: σ U l -c l = -1347.01835941669 m/s, σ Ug+cg = 404.737069698175 m/s, σ Ug = 5 m/s. Table C.2: Parameters for the EOS for the three components. Table C.3: Values for the different states of the Riemann problem.

			k = l			k = v	k = a
	C V,k (J/K/kg) 1.659732071941701 10 3 6.626564746983661 10 3	7.18 10 2
		γ k (-)	6.636214111922141	1.083834328358209	1.4000231
		Π k (P a)	3.348508243030720 10 8		0	0
		Q k (J)	0			0	0
	s 0 k (J/K/kg)	10 4	-8.646479253448585 10 4	0
		Z L		Z 1		Z 2	Z R
	y a	0.205	0.205		0.2	0.2
	α g	0.495	0.495		0.5	0.5
	ρ g (kg/m 3 ) 0.619780775226864 0.619780775226864	0.650769813988207	0.620914399831763
	U g (m/s)	5		5		5	-14.2205491976928
	P g (P a)	1.0 10 5	1.0 10 5		9.99994503590090 10 4 9.49994778410586 10 4
	ρ l (kg/m 3 ) 1221.42181799682 1221.42184547160		1221.42178476433	1221.42178476433
	U l (m/s)	2.00003034495258	2		1.96999984940252	1.96999984940252
	P l (P a)	9.995 10 4	1.0 10 5		9.98895218378901 10 4 9.98895218378901 10 4

Acknowledgements

The last author received financial support by ANRT through an EDF/CIFRE grant number 2017/0476. Computational facilities were provided by EDF. The authors thank Jean-Marc Hérard for fruitful discussions.

Appendix A. Mixture EOS for the phase g

In this appendix, we deduce the EOS of the gas phase from the EOS of phases a and v and from temperature equilibrium assumption, as for instance in [START_REF] Gallouët | Some recent finite volume schemes to compute euler equations using real gas EOS[END_REF][START_REF] Helluy | Interpolated pressure laws in two-fluid simulations and hyperbolicity[END_REF][START_REF] Bachmann | A simple model for cavitation with noncondensable gases[END_REF]. We first consider the general case and the specific case of Stiffened Gas EOS is then examined. For this particular EOS, it can be proved that the mixture EOS for the gas phase corresponds to a Stiffened Gas EOS for a fixed incondensable fraction y a = ρ a /ρ g .

Appendix A.1. The mixture EOS for the gas phase: general case Before all, we assume that the thermodynamic behavior of each phase k = {l, v, a} is given by the mean of its specific entropy s k as a function of the specific internal energy e k and the specific volume

We examine now the properties of the mixture of the phases.

The gas fraction α g is related to the vapor fraction α v , to the incondensable fraction α a and to the liquid fraction α l through the relations: α g = α v = α a and α l + α g = 1. The former arises from the miscibility of the phases v and a, see for instance [START_REF] Mathis | A thermodynamically consistent model of a liquid-vapor fluid with a gas[END_REF][START_REF] Hérard | A three-phase flow model with two miscible phases[END_REF]; whereas the latter represents a statistical conservation property of the whole mixture [START_REF] Hérard | Une approche bifluide statistique de modélisation des écoulements diphasiques à phases compressibles[END_REF][START_REF] Liu | Contribution to the verification and the validation of an unsteady two-phase flow model[END_REF][START_REF] Lochon | Modélisation et simulation d'écoulements transitoires eau-vapeur en approche bifluide[END_REF].

If we denote respectively by M k (in kg), V k (in m 3 ) and E k (in J) the mass, volume and energy of phase k ∈ {a, v}, we have for the gas phase:

The latter equalities directly traduce that the phases a and v are immiscible in a perfect manner, so that they occupy the same volume [START_REF] Mathis | A thermodynamically consistent model of a liquid-vapor fluid with a gas[END_REF][START_REF] Hérard | A three-phase flow model with two miscible phases[END_REF]. Thus the density of phase g is:

and, in the same way, we get that the volumetric internal energy of phase g, ρ g e g = E g /V g with e g = E g /M g the specific internal energy, reads:

To these relations, one must add the temperature equilibrium assumption between phases a and v which reads:

Let us be more explicit on all these variables and the closures of the model. If one assumes that the set of the conservative variables:

is the unknown of system (1), all the other quantities have to be defined with respect to W . We focus here on the thermodynamical part of the model. Following the assumptions proposed above, we have:

(A.5)

Appendix A.2. The mixture EOS for the gas phase: Stiffened Gas phasic EOS We consider from now that the two phases a and v follow a Stiffened Gas EOS, and thus that their specific entropies are defined by [START_REF] Feldman | The speed of sound in single-component two-phase fluids: theoretical and experimental[END_REF]. For k ∈ {a, l, v}, the phasic pressure and temperature is then:

The pressure can also be written in terms of the temperature and specific volume:

Moreover, for k ∈ {a, l, v} the sound speed c k is given by:

Let us define the mass fraction y a of phase a and the energy fraction z a of phase a: y a = ρ a /ρ g , z a = ρ a e a /(ρ g e g ) = y a e a /e g .

We recall that the quantities of the phase g are given through the relations:

First of all, it should be noted that the energy fraction z a is not an unknown: it can be expressed as a function of y a , τ g and e g . Indeed, the temperature equilibrium relation (A.6) can be written using the relations above:

and, hence, when τ a , τ g and e g are given, z a is obtained through the temperature equilibrium. When considering the Stiffened Gas EOS, z a can be explicitly written. It reads:

When summing the temperature formulae and accounting for the temperature equilibrium

The pressure P g can be obtained by summing the formula (A.12) for the pressures P k with T g = T a = T v , we then get:

Appendix C. Verification test case on the basis of a Riemann problem

The aim of this section is to present a verification test case for the convective part of the model proposed in section 2. The associated system is hyperbolic and shock waves are defined in a unique manner. We thus propose here to solve a Riemann problem. The solution of this Riemann problem is then used for comparison with the approximate solutions computed by the VFRoe-ncv scheme [START_REF] Buffard | A sequel to a rough Godunov scheme: application to real gases[END_REF] using the variable Y (see section 2.2.1).

Building Riemann problems for two-fluid models based on the Baer-Nunziato model can be tricky when considering a general problem defined by the right and left states. We thus follow a classical way to build a Riemann problem: one of the states is given and the other states are computed by choosing the complete waves structure. The computation of a new state from a known state (across a wave) is performed on the basis of jump relations or Riemann invariants depending on the nature of the wave that separates these two states. We restrict here to shock waves and contact waves. Indeed, rarefaction waves require more complex computations which may lead to larger computation errors and inaccurate comparisons. The Riemann problem chosen here is composed of two intermediate states Z The properties of the convective part of the model are recalled in section 2.2.1. In the shock waves, the fraction α g is constant so that the classical single-phase Rankine-Hugoniot relations are recovered for each phase l and g:

for k = {l, g}. For the contact wave, the connection between the states Z 1 and Z 2 is performed