Olivier Hurisse
email: olivier.hurisse@edf.fr

Numerical simulations of steady and unsteady two-phase flows using a homogeneous model

Keywords: Two-phase flow, homogeneous model, mass transfer, mixture sound speed, depressurization, heat transfer

The capability of a homogeneous model to simulate steady and unsteady two-phase flows is investigated. The latter is based on the Euler set of equations supplemented by a complex equation of state describing the thermodynamical behavior of the mixture. No equilibrium assumption is made except for the kinematic equilibrium. The return to the thermodynamical equilibrium is ensured by three source terms that comply with the second law of thermodynamics. The numerical code built on the basis of this model has been verified and some validation results are discussed here. The speed of propagation of a pressure signal is first studied and compared with experimental measurements. Then a more complex situation is investigated: SUPERCANON experiment which corresponds to a sudden depressurization of heated water (associated to a Loss Of Coolant Accident, or LOCA). At last, the results of a numerical experiment of heating of flowing water in a pipe are compared to those obtained with an industrial code.

Introduction

When two different phases of the same component coexist in a stable manner, the thermodynamical theory states that their temperatures, pressures and chemical potentials must be equal [START_REF] Kubo | Thermodynamics: an advanced course with problems and solutions[END_REF]. Conversely, when the equilibrium between the two phases is not yet reached, these temperatures, pressures and chemical potentials may differ. Most of the models that are used to perform two-phase flow simulations involving mass transfer are based on one or more equilibrium assumptions. This can be restrictive when dealing with highly unsteady flows for which these assumptions may not be relevant. The homogeneous model used in the present work [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF][START_REF] Faccanoni | Modelling and simulation of liquid-vapor phase transition in compressible flows based on thermodynamical equilibrium[END_REF][START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF][START_REF] Jung | Numerical simulations of two-fluid flow on multicores accelerator[END_REF] only makes the assumption that the velocities of the two phases are equal and the full thermodynamical disequilibrium is accounted for. This hypothesis of equal velocities can show some limitations -for instance when gravity plays an important role, or when simulating jets, etc... -but it remains acceptable for a lot of industrial applications with forced convection in pipes, as those encountered for nuclear power plants. One can note that some two-fluid models are developed without equilibrium assumptions; for instance those in [START_REF] Bdzil | Two-phase modeling of deflagrationto-detonation transition in granular materials: A critical examination of modeling issues[END_REF][START_REF] Berry | Relap-7 theory manual[END_REF][START_REF] Coquel | Modelling compressible multiphase flows[END_REF][START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF][START_REF] Jin | Compressible two-pressure two-phase flow models[END_REF][START_REF] Kapila | Two-phase modeling of DDT: Structure of the velocity-relaxation zone[END_REF][START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flashing flows[END_REF] which are based on 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], or in [START_REF] Gavrilyuk | Mathematical and numerical modeling of two-phase compressible flows with micro-inertia[END_REF]. It should also be mentioned that homogeneous models have been studied with different velocities for the two phases [START_REF] Bilicki | The effect of three closures on critical conditions in two-phase flow with unequal phase velocities[END_REF].

The homogeneous model in the present work was first proposed in [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF]. It is built on the Euler set of equations supplemented by three fractions which allow to account for the temperature, pressure and chemical potential disequilibria between the two phases. The return to equilibrium is then ensured by three source terms which are defined in order to fulfill the second law of thermodynamics. The resulting model possesses a complex Equation Of State (EOS) for the mixture and is hyperbolic provided that each phasic entropy is concave with respect to the phasic specific volume and to the phasic specific internal energy, and that the mixture temperature is positive [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF]. Moreover, the transition between single-phase flow and two-phase flow is managed quite naturally. This is a crucial feature since a lot of industrial situations involve such transitions.

Among other transient phenomena involving two-phase flows in confined or pressurized devices, one may for instance cite: water hammers [START_REF] Gale | Water hammer in elastic pipes[END_REF][START_REF] Streeter | Hydraulic transients[END_REF], Boiling Liquid Expanding Vapor Explosion (BLEVE) [START_REF] Baker | Guidelines for vapor cloud explosion, pressure vessel burst, bleve, and flash fire hazards[END_REF][START_REF] Xie | Thermodynamic and Gasdynamic Aspects of a Boiling Liquid Expanding Vapour Explosion[END_REF], breaches in pressurized pipes, erosion due to the collapse of cavitation bubbles [START_REF] Lauterborn | Experimental investigations of cavitation-bubble collapse in the neighbourhood of a solid boundary[END_REF], etc... In such situations, pressure waves are produced that may damage the installation. Hence, the prediction of the propagation of such waves is important for safety studies. It is well-known that the propagation of pressure waves in two-phase flows strongly depends on the gas fraction [START_REF] Karplus | The velocity of sound in a liquid containing gas bubbles[END_REF]. Roughly speaking, the speed of sound in a medium depends on the ratio of its specific volume to its compressibility. In a homogeneous mixture of liquid and gas, the compressibility is close to the compressibility of the gas, whereas the specific volume is close to the specific volume of the liquid. As a consequence, the sound of speed will be lower than that in the pure phases [START_REF] Feldman | The speed of sound in single component two phase fluids: Theoretical and experimental[END_REF]. This behavior has been investigated with experimental measurements for liquid-gas mixtures (without mass transfer). The aim of the present paper is to evaluate the capability of a code to reproduce some schematic situation encountered in nuclear applications. We first compare the speed of propagation of a pressure pulse for different air-water mixture fractions to the measurements of [START_REF] Karplus | The velocity of sound in a liquid containing gas bubbles[END_REF]. A more complex situation of depressurization of heated water is also investigated. In the latter, waves are generated by the sudden depressurization and travel in both the pure liquid and in the steam-water mixture. We also investigate the case of the heating of water flowing in a pipe.

The whole code that is used in this study possesses very few physical parameters, and the time-scale associated to the return to the thermodynamical equilibrium is the most tricky to determine. This timescale should be chosen on the basis of physical phenomena but this is a difficult task and very few work are available on that subject when considering a complete two-phase flow model [START_REF] Bilicki | Evaluation of the relaxation time of heat and mass exchange in the liquid-vapour bubble flow[END_REF][START_REF] Gavrilyuk | The structure of pressure relaxation terms: one-velocity case[END_REF].

The schemes that provide the unsteady numerical approximations are described in [START_REF] Helluy | Verification of a two-phase flow code based on an homogeneous model[END_REF][START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF]. An important verification process has been done and it has been reported in [START_REF] Helluy | Verification of a two-phase flow code based on an homogeneous model[END_REF][START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF]. In these references, the verification test cases have been chosen to mimic the main situations encountered in the nuclear domain. Two kinds of problems have thus been studied: Riemann problems which are classical unsteady test problems, and a steady-state problem involving the heating of a two-phase mixture. The former (resp. latter) are representative problems of the situation induced by the experiment proposed in section 6 (resp. 5). The case of section 5 corresponds to the computation of a steady-state. In order to obtain accurate approximations with little CPU-time, we propose in section 3 a steady-state algorithm based on the idea of [START_REF] Hurisse | A new verification test case for the compressible module of Code Saturne based on a non trivial steady-state solution for the Euler system with energy[END_REF].

An overview of the model and the numerical schemes is first proposed in section 2 and 3. In section 4, we assess the computation of the propagation of a pressure-pulse by our code. This is done on the basis of the measurements of the speed of propagation of a pressure-pulse in air-water reported in [START_REF] Karplus | The velocity of sound in a liquid containing gas bubbles[END_REF]. Section 5 is then devoted to the study of a case of heating of water in a pipe involving small pressure variations. This case is inspired from the work [START_REF] Bartolomei | Heat and Mass Transfer IV (in russian), Minsk[END_REF][START_REF] Bartolomei | An experimental investigation of true volumetric vapour content with subcooled boiling in tubes[END_REF]. It is schematic of the steam production in a steam generator of a Pressurized Water Reactor (PWR) and we compare results obtained by the model to the results obtained with an industrial code. We then focus in section 6 on a more complex situation : the SUPERCANON [START_REF] Riegel | Contribution à l'étude de la décompression d'une capacité en régime diphasique[END_REF] experiment, which could be assimilated to a two-phase shock-tube with mass transfer. These experiments where conducted in the late 70s on this experimental facility to measure the depressurization of heated water. They were designed to be representative a of a Loss Of Coolant Accident (LOCA) in the primary circuit of a PWR nuclear power plant, as the Edwards'pipe blowdown experiments [START_REF] Edwards | Studies of phenomena connected with the depressurization of water reactors[END_REF].

A homogeneous model for two-phase flows

The aim of this section is to provide a quick overview of the model and of the numerical schemes. For a more detailed presentation, the reader could refer to [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF][START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF][START_REF] Jung | Numerical simulations of two-fluid flow on multicores accelerator[END_REF] for the model, to [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF] for the numerical schemes and to [START_REF] Helluy | Verification of a two-phase flow code based on an homogeneous model[END_REF][START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF] for the verification process. In the following, a subscript v will denote a vapor or gas quantity, and a subscript l a liquid one.

The quantities describing the mixture of the two phases are: the specific volume τ = 1/ρ where ρ is the density, U is the velocity and e is the specific internal energy. The three fractions define the way the two phases are mixed in terms of: the volume through the vapor volumic fraction α v , the mass through the vapor mass fraction y v and the energy through the vapor energy-fraction z v . In fact these fractions allow to express the phasic quantities in terms of the mixture quantities:

τ v = α v τ y v , τ l = α l τ y l , e v = z v e y v , e l = z l e y l , U v = U, U l = U, (1)
where

α l = 1 -α v , y l = 1 -y v and z l = 1 -z v .
The set of equations for the homogeneous model is then:

                     ∂ ∂t (ρY ) + ∂ ∂x (ρU Y ) = ρΓ Y ∂ ∂t (ρ) + ∂ ∂x (ρU ) = 0 ∂ ∂t (ρU ) + ∂ ∂x ρU 2 + P = 0 ∂ ∂t (ρE) + ∂ ∂x (ρU E + U P ) = 0, (2)
where Y = (α v , y v , z v ), and where E stands for the specific total energy:

E = e + U 2 /2.
The thermodynamical closure of the system is ensured by the definition of the mixture entropy:

s = y v s v (τ v , e v ) + y l s l (τ l , e l ), (3)
where the phasic entropies s k (τ k , e k ) must be specified by the user. The phasic pressures P k and the phasic temperatures T k are deduced from the phasic entropies through the phasic Gibbs relations

T k ds k = de k + P k dτ k , (4)
which lead to the definitions:

T -1 k = ∂ ∂e k (s k (τ k , e k )) |τ k , P k = T k ∂ ∂τ k (s k (τ k , e k )) |e k . (5)
The pressure law P and the temperature law T for the mixture are obtained by writing the Gibbs relation [START_REF] Bdzil | Two-phase modeling of deflagrationto-detonation transition in granular materials: A critical examination of modeling issues[END_REF] for the mixture, which involves the mixture entropy (3):

T ds = de + P dτ + ∂s ∂α v |e,τ,yv,zv dα v + ∂s ∂y v |e,τ,αv,zv dy v + ∂s ∂z v |e,τ,αv,yv dz v . (6)
We then have the definitions:

1 T = ∂ ∂e (s) |τ,Y and P T = ∂ ∂τ (s) |e,Y , (7)
for the pressure and the temperature of the mixture. Using the phasic pressure P k and temperature T k , they can be written:

P (Y, τ, e) = α l T l P l + αv Tv P v z l T l + zv Tv and 1 T = z l T l + z v T v . (8)
The three source terms Γ Y (one for each fraction) rule the thermodynamical exchange between the phases and allow the system to return to the thermodynamical equilibrium. Hence they must be chosen to comply with the second principle of thermodynamics. They are written here as in [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF][START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF]:

Γ Y = Y eq -Y λ , (9)
where λ is a characteristic time-scale, and Y eq = (α v,eq , y v,eq , z v,eq ), [START_REF] Jin | Compressible two-pressure two-phase flow models[END_REF] defines the equilibrium fractions. This equilibrium state Y eq maximizes the specific entropy of the mixture s for a given specific internal energy e and a given specific volume τ . When the two phases coexist Y eq belongs to ]0, 1[ 3 and the derivative of the entropy with respect to the mixture fractions Y is then null. A simple calculus of this derivative shows that this means that the temperatures, pressures and chemical potentials are equal.

The sound speed of system ( 2) is defined as:

c 2 = -τ 2 ∂ ∂τ (P ) |s = -τ 2 ∂ ∂τ (P ) |e + τ 2 P ∂ ∂e (P ) |τ , (11)
and, using the formula (8) for the mixture pressure P , it can be written:

- c 2 T τ 2 = 1 y v (-α v , P z v ) .(d 2 s v ). -α v P z v + 1 y l (-α l , P z l ) .(d 2 s l ). -α l P z l (12)
where d 2 s k stands for the Hessian matrix of the phasic entropies s k (τ k , e k ):

d 2 s k =       ∂ 2 (s k ) ∂τ k ∂τ k ∂ 2 (s k ) ∂τ k ∂e k ∂ 2 (s k ) ∂τ k ∂e k ∂ 2 (s k ) ∂e k ∂e k .       (13)
We recall that the phasic sound speeds are defined as:

c 2 k T k τ 2 k = -(-1, P k ) .s k . -1 P k (14)
It must be emphasized that this mixture celerity c is thus not a barycenter of the phasic celerities c k . This point will be discussed in section 4 on the basis of numerical results. However, if (i) each phasic entropy is concave with respect to the phasic specific volume and to the phasic specific internal energy, (ii) and the mixture temperature is positive, then c 2 is non-negative, and hence c is real. Thus the conditions (i) and (ii) provide a sufficient condition for system (2) to be hyperbolic.

We restrict hereafter to Stiffened Gas EOS, that is:

s k (τ k , e k ) = C v,k ln (e k -Q k -Π k τ k )τ γ k -1 k + s 0 k ( 15
)
where C v,k is the heat capacity, Q k is a reference value of the enthalpy h k = e k + P k τ k , the parameter -Π k corresponds to the minimal pressure1 , γ k is the adiabatic coefficient and s 0 k is the reference entropy. This phasic entropy allows to define the pressure and temperature laws through the phasic Gibbs relations (4):

P k (τ k , e k ) = e k -Q k τ k (γ k -1) -Π k γ k , (16)
and

T k (τ k , e k ) = e k -Q k -Π k τ k C v,k . ( 17
)
3 Numerical schemes

We propose in this section a quick overview of the numerical schemes which are detailed in D and C.

Numerical schemes for the unsteady simulations

The overall numerical scheme is a fractional step method, as proposed in [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF]: the convection part of the system is first solved using a Godunov type explicit scheme; and the source terms are then integrated. As proposed in [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF] a VFRoe-ncv scheme is used to compute the numerical fluxes at the interface between two cells of the mesh. The linearised problem at the interface is solved by considering the variables (Y, τ, U, P ) [START_REF] Buffard | A sequel to a rough godunov scheme: application to real gases[END_REF] and, in order to avoid non-entropic waves to be computed, this algorithm is supplemented with the entropic correction proposed in [START_REF] Helluy | A simple parameter-free entropy correction for approximate riemann solvers[END_REF]. Moreover, a partial WFRoe scheme [START_REF] Gallouët | A new approximate godunov scheme with application to dense gassolid flows[END_REF] is performed to improve the prediction of the speed of the contact wave in the linearised problem at the interface between two cells. In order to improve the accuracy of this scheme (particularly for regular solutions), a modification has been proposed in [START_REF] Helluy | Verification of a two-phase flow code based on an homogeneous model[END_REF] on the basis of the idea of [START_REF] Dellacherie | Analysis of godunov type schemes applied to the compressible euler system at low mach number[END_REF]. The interfacial pressure issued from the linearized Riemann problem is thus mixed with a centered pressure chosen as the mean of the pressure of the left and right cells. The details of the scheme are given in D, and it is used in section 4 and 6.

Numerical scheme for the steady-state computations

We propose here a numerical scheme devoted to the computation of the approximations of the onedimensionnal steady-state problem associated with [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF], that is ∀x ∈ [0, L] (L > 0):

                       d dx (ρU Y ) = ρΓ Y d dx (ρU ) = 0 d dx ρU 2 + P = ρF d dx (ρU E + U P ) = Φ. ( 18
)
The source terms F and Φ respectively denote the volumic forces acting on the system and the power received by the system. In particular, the term Φ gathers the power associated with the forces F and the external heating sources Ψ: Φ = ρU F + Ψ. The details of the scheme are given in C and it is used in section 5.

4 The speed of propagation of a pressure-pulse in a mixture of air and water

When considering a two-phase flow model, an important question that arises is to investigate if the speed of propagation of a pressure-pulse in the mixture is physically relevant. This is indeed an important point for two-phase flow simulation because this speed determines the speed of propagation of the pressure overshoots that could damage industrial devices. It is well-known that the speed of propagation of a pressure-signal in a mixture is not the barycenter of the phasic sound speeds. In particular, this speed of propagation can be far less than the phasic sound speed in the gas phase when considering intermediate volume fractions [START_REF] Karplus | The velocity of sound in a liquid containing gas bubbles[END_REF]. We thus focus here on the air-water configuration for which experimental results are available in [START_REF] Karplus | The velocity of sound in a liquid containing gas bubbles[END_REF]. In order to make such comparisons, the test case proposed in [START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF][START_REF] Liu | Contribution to the verification and the validation of an unsteady two-phase flow model[END_REF] is reproduced to calculate the speed of propagation of a small pressure-perturbation.

The test case consists in calculating the speed of propagation of a pressure perturbation in a mixture of air and water at different volume-fractions α v ∈ [0, 1]. The computational domain is 1 m long and the mixture is at 1 bar and 293 K. For this pressure and temperature, the sound speed in pure water (α v = 0) is 1480 m/s and 343 m/s in pure air (α v = 1). The initial conditions are thus:

α v ∈ [0, 1], T = 293.0 K, P = 10 5 + 10 3 exp -(x-0.5) 2 σ 2 P a; (19)
with σ = 0.01 m. The other variables are deduced from the previous one [START_REF] Xie | Thermodynamic and Gasdynamic Aspects of a Boiling Liquid Expanding Vapour Explosion[END_REF] by imposing the temperature and pressure equilibrium (see appendix A for details about the modifications of the model of section 2 when the mass transfer is not accounted for). With these initial conditions [START_REF] Xie | Thermodynamic and Gasdynamic Aspects of a Boiling Liquid Expanding Vapour Explosion[END_REF], two pressure peaks travel towards the left and towards the right outlets. Their speed C peak is computed by measuring at time t the distance d(t) from the peak to the middle of the domain (x = 0.5, which is the initial location of the peak) and with the obvious formula:

C peak = d(t)/t.
The results are plotted on figure 1 when considering different configurations for the relaxation timescale. Simulations have been carried for different relaxation time-scales λ = {10 -20 s, 10 -7 s, 10 -6 s, 10 -5 s, 10 -4 s}.

The "frozen" sound speed corresponds to the value of the celerity [START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flashing flows[END_REF]. This sound speed corresponds to the celerity computed for the convective part of system 2, thus omitting the effect of the thermodynamical relaxation due to the source terms. With the 2000-cell mesh used for the simulation, the time-step always remains in the range [1.5 10 -7 s, 4.0 10 -7 s]. The smaller time-scale thus corresponds here to an instantaneous relaxation to the thermodynamical equilibrium. When the time-scale is increased above 10 -5 s, three peaks can be observed for intermediate volume fractions. Figure 2 shows the pressure at time 10 -4 s for different values of the time-scale λ. The central peak that can be observed for λ ≤ 10 -5 s can afterwards separates into two other peaks while the two small peaks tend to decrease. The transition between the two kinds of solutions is fast with respect to λ. Roughly speaking, λ can be gathered into two classes : those that lead to a speed of propagation close to that of λ = 10 -20 , and those that lead to a speed of propagation close to the frozen celerity.

The heating of water flowing in a pipe

We are interested here in the computation of steady states by using the method proposed in C, which allows to obtain approximate solutions of the system [START_REF] Baker | Guidelines for vapor cloud explosion, pressure vessel burst, bleve, and flash fire hazards[END_REF]. For this purpose, we choose to mimic the BARTOLOMEI test case [START_REF] Bartolomei | Heat and Mass Transfer IV (in russian), Minsk[END_REF][START_REF] Bartolomei | An experimental investigation of true volumetric vapour content with subcooled boiling in tubes[END_REF] by the 1D test case described on figure [START_REF] Faccanoni | Modelling and simulation of liquid-vapor phase transition in compressible flows based on thermodynamical equilibrium[END_REF]. Liquid water at T 0 = 495 K is injected in the pipe with a flow rate of D 0 = 1500 m 3 /s with a pressure of P 0 = 68.9 bars. The fluid is heated until x = 1 m and, depending on the value of the external power heating Ψ 0 , steam may be created. The external forces F are set to zero here. The thermodynamical parameters of the phasic EOS used in this section are given in B.3. In section 5.1 we propose to study the influence of the mesh refinement and of the relaxation time scale on the results computed using the method of C. Then in section 5.2, we compare these results with those obtained with the industrial code TH1D [START_REF] Baurin | Module TH1D V1.0 : Note de principe[END_REF]. The 1D test case proposed here can not reproduce the complexity of the real BARTOLOMEI experiment, thus comparisons with the experimental data would be meaningless. We thus choose a code-to-code comparison.

Mesh refinement and relaxation to the thermodynamical equilibrium

We first study the evolution of the approximated solutions when the mesh is refined. We set Ψ 0 = 5.5 10 8 W/m 3 and λ = 10 -8 s and we compare the approximated solutions computed for 40, 80, 160, 320, 640 and 1280 cells. Results are plotted on Figure 4. It can be seen that the approximations barely vary for meshes with more than 160 cells, at least at the visualization scale. In the following, we will thus use a mesh with 300 cells which gives a very satisfactory accuracy. Moreover, in this case, U < 5 m/s and the mesh size associated with λ is then:

d λ = U λ = 5 10 -8 m.
If a cell is greater than d λ , then the fluid remains in this cell during a time that is larger than λ. Hence the thermodynamical equilibrium is achieved in this cell. The dimension d λ = 5 10 -8 m corresponds to a uniform mesh of the domain [0, 1.4 m] that contains more than 28 millions of cells. As a consequence, the thermodynamical equilibrium is always achieved for the meshes proposed on Figure 4.

We refer in the following lines to Figure 4. Pure liquid enters the domain at x = 0 and a two-phase mixture flows to the right outlet at x = 1.4. We can observe three zones. This would be different with a greater time-scale as shown on Figure 5. It can be observed that with λ = 10 -4 s the steam production is slightly delayed and that it continues after x = 1 m. It should also be noticed that small time-scales λ lead to less smooth approximate solutions.

Comparisons with the results obtained with TH1D

We compare here the results obtained by the method proposed in C and the results obtained by using the TH1D (TermoHydraulic in 1D) module [START_REF] Baurin | Module TH1D V1.0 : Note de principe[END_REF]. The latter is a one-dimensionnal steady-state version of the code THYC [START_REF] Aubry | THYC Code de Thermohydraulique des Coeurs de Racteurs Version 1 .0. Modélisation et Méthodes Numériques[END_REF][START_REF] Aubry | The THYC three-dimensional thermal-hydraulic code for rod bundles: recent developments and validation tests[END_REF][START_REF] Guelfi | THYC (ThermoHYdraulique des Composants) Version 4.1 -Note de Principe[END_REF][START_REF] Le Coq | The THYC computer code. A finite volume approach for 3 dimensional two-phase flows in tube bundles[END_REF] (TermoHydraulic of the Core) and it is based on the THETIS [START_REF] Rascle | Spécifications fonctionnelles de THETIS. Calcul par interpolation des fonctions thermodynamiques des fluides diphasiques[END_REF] tabulated thermodynamical laws. The THYC code is based on a two-phase flow model (of the so-called Homogeneous Relaxed Model type) for which the mixture law is obtained by assuming the pressure equilibrium and by assuming that the vapor phase is always at saturation, as in [START_REF] Bilicki | Relaxation models for wave phenomena in liquid-vapor bubble flow in channels[END_REF][START_REF] Bilicki | The effect of three closures on critical conditions in two-phase flow with unequal phase velocities[END_REF]. The TH1D computations presented here are part of the work reported in [START_REF] Baurin | Consolidation et validation d'un module de thermo-hydraulique pour les calculs des coeurs de centrales REP[END_REF]. We then set the relaxation time-scales to zero for both codes, and all the sub-models of TH1D are disabled, in particular :

• the relative velocity between the phase is maintained equal to zero;

• the specific source term used to produce steam in the mass fraction equation in the case of wall heating is disabled (the sole volumic heating source term is kept);

The results are then compared for different heating powers: 4.0 10 8 W/m 3 , 4. The global behavior of the two codes is the same, but some quantitative differences can be noted.

It must be emphasized that in TH1D the reference pressure 6.89 10 5 P a is imposed at x = 1.4 m, whereas this pressure is imposed at x = 0 for the code of C. This difference leads to a shift in the two pressure profiles. The influence of this shift on the thermodynamical results (density, temperature) is not very important. Indeed, the relative variations of the pressure are small and the associated variation of the saturation temperature (the dashed lines on the figures below) is negligible with respect to the variation of the temperature due to the heating.

For the low heating powers 4.0 10 8 W/m 3 and 4.5 10 8 W/m 3 (see resp. Figures 6 and7), no vapor is created by the two codes. The results are almost identical, which tends to prove that the liquid states are correctly represented by the Stiffened Gas EOS. In particular the two temperatures are very close which implies that the saturation is reached at almost the same coordinate x (when the heating power is sufficient). Nevertheless, some differences can be noted in the amount of vapor created and therefore in the mixture density (see . It can also be noted that the differences tend to decrease for high heating powers (see Figures 12).

The quantitative difference observed with these numerical simulations is a consequence of the two different thermodynamical laws used to described the two-phase mixture. The comparison proposed here between the method of C and TH1D allows to estimate the difference between the two thermodynamical representations of the two-phase flows. A difference can be noted for the vapor volume fraction and the mean density, at least for the intermediate heating powers. This may indicate that the Stiffened Gas EOS used here are satisfactory, at least for this test case which involves very small variations of the pressure.

Sudden depressurization of heated water

The present test case is associated to 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 13 is a sketch of the facility. A tube (100 mm of inner diameter and 4.389 m long) is filled with water and 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 and is supposed 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 13). 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 on figure 14 for two different experimental runs.

In 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 λ. Unfortunately, we are actually not able to derive a law for λ on the basis of physical considerations. For this reason, we propose here two different laws for the relaxation time-scale λ. These laws are built following the observations reported in [START_REF] Bilicki | Evaluation of the relaxation time of heat and mass exchange in the liquid-vapour bubble flow[END_REF], on the basis of an experiment of a flashing (steady) flow in a divergent duct. These laws combined two terms :

• one term that diminishes the relaxation time-scale when the flow is too far from equilibrium,

• and one term that decreases the relaxation time-scale when the steam fraction increases.

The laws that are proposed here are also built on such terms, but their exact form is modified to account for strong disequilibria. Hence the laws below are not lying on mechanical assumptions and they must be seen as "toy laws" whose parameters have been chosen to fit the experiments. The three laws are the following.

(i) λ is uniform and equal to 10 -30 s. With this choice the two phases always remain at the thermodynamical equilibrium and the whole model behaves like a Homogeneous Equilibrium Model (HEM).

It will be denoted in the following by "Equilibrium".

(ii) λ is defined as:

λ = λ 0 e - |αv -αv,eq | dα 2 e - αv αc 3 , (20)
with the parameters: λ 0 = 4 10 -4 s, dα = 3.16 10 -2 et α c = 0.35, and where α v,eq stands for the vapor volume fraction at equilibrium. This law for λ is denoted by "disequilibrium 1" in the following.

(iii) λ is defined as:

λ = λ 0 e - |αv -αv,eq | dα e - αv αc 4 , (21)
with the parameters: λ 0 = 9 10 -4 s, dα = 5 10 -4 et α c = 0.42. This law for λ is denoted by "disequilibrium 2" in the following.

In laws (ii) and (iii), the first exponential determines the "maximum distance to the saturation", whereas the second exponential involves a critical volume fraction of vapor α c above which the flow remains near the equilibrium. These two terms have been plotted on figure 15 for the two disequilibrium laws.

The gravity effect is classicaly neglected for this test case because the whole depressurisation is rather rapid. Moreover, the computational domain is x ∈ [0, 15 m]. The domain x < 4.389 m is filled with high pressurized water. Since we are dealing with a two-phase flow model, the domain x ∈ [4.389, 15 m] is filled with steam at a pressure of 1 bar and a temperature of 300 • C. In the experiment, this domain was filled with air at 1 bar and ambient temperature. The right boundary condition is a supersonic exit, which means that the fluxes at the boundary are computed considering the flow inside the pipe. This choice for the right boundary condition is legitimated by the high speed flows in the vapor domain (it reaches 600 m/s).

The numerical results obtained with the three laws (i), (ii) and (iii) are compared to some experiments on figures 16 et 17. It can be noted that for the pressure, the two different experiments show very different Figure 15: Plots of the two exponential terms involved in the laws for "disequilibrium 1" (black) and "disequilibrium 2" (red). Top : the first exponential term involving |α v -α v,eq | is plotted with respect to (α v -α v,eq ) , bottom : the second exponential term involving α v /α c is plotted with respect to α v . patterns for the time less than 0.075 s, see figure 16. This corresponds to the fact that, at point P 1 , the flow is out of equilibrium until that time. Yet the sudden initial pressure drop in the pipe generates vapor by nucleation, and depending on the water purity, the vapor generation includes homogeneous and heterogeneous nucleation. Each impurity in the water can promote nucleation and then vapor generation. As a consequence, the flow remains closer to the equilibrium, whereas with few impurities, the vapor generation occurs through homogeneous nucleation. This could explain the great variability in the measurements which can also be observed at point P t when considering the volume fraction of vapor, see figure 17.

The equilibrium law for λ gives results that are not satisfactory, which was expected since the flow is far from 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]. On the contrary, the results are in good agreement with the experiments when considering the laws (ii) and (iii) for λ, at least until 0.1 s. The predicted final pressure drop occurs too early with respect to the measurements (∼ 0.105 s versus ∼ 0.15 s). It could be due to the Stiffened Gas EOS which is built around the reference point P 0 = 85 bar, and thus it could be not accurate enough for such thermodynamical variations (from 150 bar to 1 bar). Moreover, the code does not account for head losses or pipe deformation [START_REF] Allievi | Teoria del colpo dariete (the theory of waterhammer)[END_REF][START_REF] Gale | Water hammer in elastic pipes[END_REF][START_REF] Liu | Contribution to the verification and the validation of an unsteady two-phase flow model[END_REF] which is known to slow down the propagation of the waves. This could delay the final pressure drop in our simulations. Nevertheless, it can be observed that the two "toy laws" can be fitted to retrieve the two different behaviors of the two experiments for times smaller than 1.0 s. This might indicate that the model is well-suited for the computation of sudden depressurization, provided the use of: a more physical law for λ, and probably a more complex set of phasic EOS.

Conclusion

The two-phase flow model presented in section 2 takes into account the thermodynamical behavior of each phase through a simple stiffened gas EOS. It was shown in section 4 that the predicted speed of propagation of a pressure-perturbation in a mixture is in good agreement with some experimental data. The comparisons of section 5 between the model of section 2 (using Stiffened Gas EOS) and the industrial code TH1D allow to estimate the accuracy of the Stiffened Gas EOS used for each phase. The agreement is quite good between the two codes, which tends to show that Stiffened Gas EOS can be accurate, at least when small variations of the pressure are involved. When large variations of the pressure and the temperature are observed, as in the LOCA case of section 6, the agreement between simulations and experimental results remains satisfactory even if some discrepancies are observed. The algorithms proposed to account for the thermodynamical disequilibrium (see D.2) are not restricted to Stiffened Gas EOS and it would be very interesting to extend the present simulations with more realsitic EOS for the LOCA test case. Furthermore, the sound speed can be very different along the computational domain: high sound speed in the pure liquid and low sound speed in two-phase regions. In such situations, the explicit numerical scheme used in the present work can lead to small time-steps (chosen with respect to the maximum sound-speed in the whole domain) which may lead to poor approximations in the two-phase domain. Improving the convection scheme by using semi-implicit schemes or low-mach schemes could improve the simulations for these situations. At last, modeling work has to be done to define relaxation time-scales on the basis on mechanical considerations. A Modifications of the model for two-phase flows whithout mass transfer.

When the mass transfer does not occur, the mass fraction of the mixture must remain constant. In order to disable the mass transfer in the model of section 2, the source term in the equation on y is then set to zero. Moreover, the thermodynamical equilibrium is only defined by the temperature and pressure equilibrium, since the chemical potential equilibrium does not make sense. If the mass transfer is enabled, the equilibrium fractions α v,eq , y v,eq and z v,eq are defined as the state that reaches the maximum of the mixture entropy s(Y, τ, e) for frozen e and τ . When mass transfer is disabled, the equilibrium fractions α v,eq and z v,eq are defined as the state that reaches the maximum of the mixture entropy s(Y, τ, e) for frozen e, τ and y.

Since the entropy is strictly concave [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF][START_REF] Jung | Numerical simulations of two-fluid flow on multicores accelerator[END_REF][START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF] with respect to Y , it is also concave with respect to (α v , z v ) for a given y ∈]0, 1[. Thus the equilibrium fractions are also defined in a unique manner when no mass transfer is accounted for. In this case, and for Stiffened Gas EOS, the computation of the temperature-pressure equilibrium leads to the resolution of a polynomial of degree two (for Perfect Gas, i.e. Π k = 0 for k = l, v, the computation is straightforward). It should be noted that when considering more complex phasic EOS, the computation of the equilibrium fractions may be more complex.

B Values of the parameters of the phasic EOS.

The parameters of the phasic EOS [START_REF] Bilicki | The effect of three closures on critical conditions in two-phase flow with unequal phase velocities[END_REF] are determined on the basis of a reference point (P 0 , T 0 ) and with the help of an industrial tool based on the IAPWS formulations. For each test case, the parameters are chosen with a different reference point (P 0 , T 0 ). Since we have a restricted number of parameters, we cannot fit all the thermodynamical properties. The choice of the relevant properties depends on each test case. The reference entropies s 0 k are only useful when the mass transfer between the phases is active. In that case, they play a role in the definition of the chemical potentials µ k = h k -T k s k , where h k = e k + P k τ k is the specific enthalpy. The reference entropies are then chosen to ensure that the saturation curve computed from the pressure-temperature-potential equilibrium fits the reference (IAPWS) saturation curve at P 0 .

B.1 Parameters for the air-water experiment

At (P 0 , T 0 ), the density of the liquid ρ 0 l , the enthalpies h 0 l and h 0 v and the celerities c 0 l and c 0 v are computed. They are then used to choose C v,k , Q k , Π k and γ k with the following formulae:

γ k = 1 + (c 0 k ) 2 h 0 k , C v,k = h 0 k T 0 γ k , Π l = ρ 0 l C v,l T 0 (γ l -1) -P 0 , Π v = 0. (22)
The reference entropies s 0 k and the energy of activation Q k are only useful when the mass transfer between the phases is active, hence they are set here to zero. The values used in section 4 are gathered in table 18.

P 0 (P a)

B.2 Parameters for the LOCA experiment

We use here the method proposed in [START_REF] Daude | Numerical experiments using a hllc-type scheme with ale formulation for compressible two-phase flows five-equation models with phase transition[END_REF]

B.3 Parameters for the volumic heating experiment

The aim of the test case of section 5 is to compute a steady state associated to the heating of the flow, so that no waves are computed. Unlike in the two previous sections, B.1 and B.2, the EOS parameters are chosen here to render a good thermal behavior of the liquid phase by fitting the heat capacities. The phasic celerities are not considered here for the liquid phase. For the EOS parameters of the vapor phase we apply the method proposed in B.1.

One of the main features of the test case of section 5 is that the pressure remains almost constant (i.e. the pressure variations remain small with respect to the reference pressure). The method for choosing the liquid EOS-parameters is the following. For Stiffened Gas EOS, the liquid enthalpy reads:

h l = Q l + C p,l T l , (23)
where C p,l = γ l C v,l . The parameters Q l and C p,l are chosen using a linear regression on tabulated data of the liquid enthalpy for T ∈ [495 K, 554 K] and at P = 68.9 bars. Then we choose a representative tabulated value of C v,l in the range T ∈ [495 K, 554 K]. Hence we have Q l , C v,l and γ l = C p,l /C v,l . It remains to choose Π l ; the latter is chosen so that the specific volume for the Stiffened Gas law,

τ l = (γ l -1)C v,l T P + Π l , (24)
fits the tabulated specific volume on the range T ∈ [495 K, 554 K] and at P = 68.9 bars. Finally, the reference entropies s 0 l and s 0 v are chosen to recover the saturation T sat (68.9 bars) = 554 K through the chemical potential equilibrium. This procedure leads to the parameters gathered in table 20.

C Numerical scheme for the steady-state computations

We propose here a numerical scheme devoted to the computation of the approximations of the onedimensionnal steady-state problem associated with [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF], that is ∀x ∈ [0, L] (L > 0):

                       d dx (ρU Y ) = ρΓ Y d dx (ρU ) = 0 d dx ρU 2 + P = ρF d dx (ρU E + U P ) = Φ. (25)
The source terms F and Φ respectively denote the volumic forces acting on the system and the power received by the system. In particular, the term Φ gathers the power associated with the forces F and the external heating sources Ψ: Φ = ρU F + Ψ. Thanks to the second equation of [START_REF] Helluy | Verification of a two-phase flow code based on an homogeneous model[END_REF], the mass flow rate D 0 = (ρU )(x) is uniform. We assume that D 0 > 0, and we then have to provide an inlet boundary condition at x = 0, which is denoted by: Y 0 , ρ 0 , U 0 , e 0 . Since the mass flow rate is uniform, we have:

∀x ∈ [0, L], D 0 = (ρU )(x) = ρ 0 U 0 . (26)
By replacing U = τ D 0 in [START_REF] Helluy | Verification of a two-phase flow code based on an homogeneous model[END_REF] we finally have to solve on [0, L] the problem:

               D 0 d dx (Y ) = ρΓ Y d dx D 2 0 τ + P = ρF D 0 d dx (E + τ P ) = Φ (27)
with Y (x = 0) = Y 0 , ρ(x = 0) = ρ 0 , U (x = 0) = U 0 , e(x = 0) = e 0 . The problem [START_REF] Bartolomei | Heat and Mass Transfer IV (in russian), Minsk[END_REF] only depends on the thermodynamical variables since the velocity U is (afterwards) deduced from τ by the relation

U = τ D 0 .
When the fractions Y are uniform, the system (27) can be solved analytically for F = 0 and provided that Ψ can be analytically integrated [START_REF] Helluy | Verification of a two-phase flow code based on an homogeneous model[END_REF]. For non-uniform fractions, we propose the following discrete integration.

We first assume that an approximation of the solution of ( 27) is known at x = x a ∈ [0, L[. System [START_REF] Bartolomei | Heat and Mass Transfer IV (in russian), Minsk[END_REF] is then integrated on the interval x ∈ [x a , x b = x a + δ x ] by using an implicit first-order integration for the source terms:

       D 0 (Y b -Y a ) = δ x ρ b Y eq,b -Y b λ b D 2 0 (τ b -τ a ) + (P b -P a ) = δ x ρ b F b D 0 (E b -E a ) + D 0 (τ b P b -τ a P a ) = δ x Φ b . (28)
The subscripts a and b respectively denote the values at x = x a and x = x b . System of equations ( 28 We now assume that the source term F only depends on τ and U , which thanks to the uniform mass flow rate is equivalent to assume that it only depends on τ . The second equation of ( 28) then gives explicitly the pressure P b as a function of τ b :

P b = P a + δ x ρ b F b -D 2 0 (τ b -τ a ). ( 29
)
The equilibrium fractions Y eq are functions of τ and e. Since e is a function of Y , τ and P through [START_REF] Coquel | Modelling compressible multiphase flows[END_REF], the fractions Y eq can also be seen as functions of Y , τ and P . Hence, the first equation of ( 28) can be written:

D 0 (Y b -Y a ) = δ x ρ b Y eq (Y b , τ b , P b ) -Y b λ(Y b , τ b , P b ) . (30)
Thus, when τ b is known, equation [START_REF] Riegel | Contribution à l'étude de la décompression d'une capacité en régime diphasique[END_REF] gives P b and the fraction Y b is implicitly defined by [START_REF] Edwards | Studies of phenomena connected with the depressurization of water reactors[END_REF]. At last, the third equation of ( 28) provides an estimation of e b :

D 0 (e b -e a ) + D 0 (τ b -τ a ) P a - D 2 0 (τ b -τa) 2 = δ x (Φ(Y b , τ b , P b ) -D 0 F b (τ b )) = δ x Ψ(Y b , τ b , P b ). (31)
The solution of system (28) can then be determined by finding the value of τ b such that the internal energy e b computed through the chain rules ( 29), ( 30) and (31) equals the internal energy e(Y b , τ b , P b ), where the law for e is given by (8).

To summarize, we have to find a zero of:

G(τ b ) = e b -e(Y b , τ b , P b ), (32)
where:

• P b is computed through (29) from τ b (explicit computation);
• Y b is computed through [START_REF] Edwards | Studies of phenomena connected with the depressurization of water reactors[END_REF] from τ b and P b (implicit computation);

• e b is computed through [START_REF] Buffard | A sequel to a rough godunov scheme: application to real gases[END_REF] from Y b , τ b and P b (explicit computation);

• (Y, τ, P ) → e(Y, τ, P ) corresponds to (8).

In the following, a zero of τ → G(τ ) is computed using a Newton algorithm (which is initialized with τ = τ a ), and at each iteration the fraction Y b is obtained using a fixed-point algorithm for equation [START_REF] Edwards | Studies of phenomena connected with the depressurization of water reactors[END_REF].

As an example, the function G, α b and P b are plotted with respect to τ b on figure 21 and 22 considering the initial conditions of section 5.2. The heating power is equal to 7.0 10 8 W and δx = 0.8 m (see figure 12). Two time-scales have been used λ = 10 -4 s (red dashed line) and λ = 10 -8 s (black plain line). It can be noticed on figure 21 that the function G can have a discontinuous derivative which corresponds to the phase transition from pure liquid to two-phase flow. We insist on the fact that, a priori, we are not ensured that there exists a solution or that this solution is unique when it exists. Nevertheless, since the approximate solution at x b is computed with an initial guess equal to the solution at x a (τ b,ini = τ a ), the whole procedure is rather efficient. The main drawback of this method of construction of the whole approximation on [0, L] is that the error in the approximation is transmitted from one point to the next one. Hence a great attention has to be paid to the accuracy of the non-linear solvers (Newton and fixed point methods).

In the algorithm described above, no linearization of the thermodynamic laws is used and the only approximation is carried by the source terms. These terms are discretized using a first order formula, but it could be possible to use higher order integration formulas.

At last, the computation of the equilibrium fractions Y eq is explicitely described in D.2.

D The numerical scheme used in the unsteady computations

The overall numerical scheme is a fractional step method, as proposed in [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF]: the convection part of the system is first solved using a Godunov type explicit scheme described in D.1; and the source terms are then integrated using an semi-implicit Euler scheme as depicted in D.2.

D.1 The numerical scheme for the convective part

The numerical scheme used to compute approximate solutions of the convection part of the model relies on the WFRoe-ncv scheme [START_REF] Gallouët | A new approximate godunov scheme with application to dense gassolid flows[END_REF] and VFRoe-ncv scheme using the variable (Y, τ, U, P ) [START_REF] Buffard | A sequel to a rough godunov scheme: application to real gases[END_REF]. Let us first recall the latter algorithm.

For the sake of simplicity, we restrict to regular meshes of size ∆x such that: ∆x = x i+ 1 2 -x i-1 2 , i ∈ Z. We denote ∆t the time step, where ∆t = t n+1 -t n , n ∈ N. In order to approximate solutions of the exact solution W ∈ R p of the conservative hyperbolic system:

∂ ∂t (W ) + ∂ ∂x (F (W )) = 0 W (x, 0) = W 0 (x) (33)
with F (W ) in R p . Let W n i be the approximate value of

W n+1 i = W n i - ∆t ∆x φ n i+ 1 2 -φ n i-1 2 (34)
The numerical flux φ n depends on W n i and W n i+1 when restricting to first-order schemes. The approximate Godunov flux φ(W L , W R ) is obtained by solving exactly the following linear 1D Riemann problem:

     ∂ ∂t (Z) + B( Ẑ) ∂ ∂x (Z) = 0 Z(x, 0) = Z L if x < 0 Z R otherwise (35)
with the initial condition: Z L = Z(W i ) and Z R = Z(W i+1 ). The matrix:

B(Z) = (W ,Z (Z)) -1 A(W (Z)) W ,Z (Z) (36)
(A(W ) is the Jacobian matrix of flux F (W )). Once the exact solution Z * ( x t ; Z L , Z R ) of this approximate problem ( 35) is obtained, the numerical flux is defined as: φ(W L , W R ) = F (W (Z * (0; Z L , Z R ))) [START_REF] Baurin | Module TH1D V1.0 : Note de principe[END_REF] Let us set l k , λ k and r k , k = 1, ..., p, left eigenvectors, eigenvalues and right eigenvectors of matrix B( Ẑ) respectively. The solution Z * ( x t ; Z L , Z R ) of the linear Riemann problem is defined everywhere (except along x t = λ k ):

Z * x t ; Z L , Z R = Z L + x t > λ k ( t l k .(Z R -Z L )) r k (38) = Z R - x t < λ k ( t l k .(Z R -Z L )) r k (39)
or in a slightly different form:

Z R -Z L = k=1,p ( t l k .(Z R -Z L )) r k = k=1,p αk r k (40)
setting:

αk = t l k .(Z R -Z L ) (41)
The only remaining unknown is the mean Ẑ which must comply with the condition:

Ẑ(Z l = Z 0 , Z r = Z 0 ) = Z 0 (42)
The standard average which is used is:

Ẑ(Z L , Z R ) = (Z L + Z R )/2 (43)
The explicit form of the approximate Godunov scheme will be under conservative form:

W n+1 i -W n i + ∆t ∆x F (W (Z * (0; Z n i , Z n i+1 ))) -F (W (Z * (0; Z n i-1 , Z n i ))) = 0 (44)
This concludes the definition of the VFRoe-ncv scheme [START_REF] Buffard | A sequel to a rough godunov scheme: application to real gases[END_REF].

A different prediction may be obtained using instead:

Ẑ(Z L , Z R ) = Z * (0; Z L , Z R ) ( 45
)
where the approximate value at the interface Z * (0; Z L , Z R ) is obtained solving [START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF] with:

Ẑ(Z L , Z R ) = (Z L + Z R )/2 (46)
This corresponds to the WFRoe-ncv scheme [START_REF] Gallouët | A new approximate godunov scheme with application to dense gassolid flows[END_REF].

The sheme used in this study is based on the variable change Z = (Y, τ, U, P ). Moreover, the velocity value of the average state Ẑ(Z L , Z R ) = (Z L + Z R )/2 is replaced by the value of the velocity U * (0; Z L , Z R ) which is computed through [START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF] as in the WFroe-ncv scheme.

At last, the pressure P * (0; Z L , Z R ) obtained by solving [START_REF] Crouzet | Validation of a two-fluid model on unsteady liquid-vapor water flows[END_REF] is blended with the centred pressure (P L + P R )/2 as proposed in [START_REF] Dellacherie | Analysis of godunov type schemes applied to the compressible euler system at low mach number[END_REF]

where C L and C R are respectively the left and right sound speeds. We have chosen the function:

f b (M ) = 1-cos(ΠM/M lim ) 2
, if M ≤ M lim , 1, otherwise.

(48)

The modified interfacial pressure P * LR used to compute the numerical fluxes is then:

P * LR = f b (M )P * (0; Z L , Z R ) + (1 -f b (M )) P L + P R 2 ; ( 49
)
where P L and P R are the left and right pressures. The parameter M lim is the limit mach number above which the pressure correction is not activated.

D.2 The numerical scheme for the source terms

The source terms of the system (2) are accounted for by discretizing the system of ODEs: We first remark that it can be written in an equivalent manner: It can thus be noticed that if the parameter λ is constant, the system (51) can be integrated exactly. We then approach system (51) by replacing λ(t) by its value λ(0). It yields for the fractions:

                                 ∂ ∂t (
                                  
               ∂ ∂t (α(t)) =
α eq (τ (0), e(0)) -α(t) λ(0) ∂ ∂t (y(t)) = y eq (τ (0), e(0)) -y(t) λ(0) ∂ ∂t (z(t)) = z eq (τ (0), e(0)) -z(t) λ(0)

The numerical approximation W n+1 is the exact solution of the system (52) at time t = ∆t and with the initial condition W n+1, * , the values obtained after the convection step. The final approximation W n+1

Figure 1 :Figure 2 :

 12 Figure1: Speed of propagation of the pressure-perturbation C peak (in m/s) with respect to the air volumefraction α v , comparison with the measurements of[START_REF] Karplus | The velocity of sound in a liquid containing gas bubbles[END_REF]. The figure on the bottom is a zoom around the experimental measurement points from[START_REF] Karplus | The velocity of sound in a liquid containing gas bubbles[END_REF].

Figure 3 :

 3 Figure 3: Description of the test case .

  (i) On [0, 0.74 m], the temperature of the liquid increases and the pressure slightly decreases. The liquid reaches the saturation temperature at 0.74 m.(ii) On [0.74 m, 1 m], steam is instantaneously created at saturation temperature and the pressure decreases more rapidly (since the time scale λ is small).(iii) Beyond the heating zone [1 m, 1.4 m], nothing happens and all the variables remain close to the value they have reached at x = 1 m.

Figure 4 :

 4 Figure 4: Approximated solutions along the domain x ∈ [0, 1.4 m] for different meshes: 40, 80, 160, 320, 640 and 1280 cells.

Figure 5 :

 5 Figure 5: Approximated solutions along the domain x ∈ [0, 1.4 m] for 300 cells and different values of the time-scale : λ = 10 -8 s (black) and λ = 10 -4 s (red). On the temperature plot, the dashed lines represent the saturation temperature.

  5 10 8 W/m 3 , 5.0 10 8 W/m 3 , 5.5 10 8 W/m 3 , 6.0 10 8 W/m 3 , 6.5 10 8 W/m 3 and 7.0 10 8 W/m 3 . Figures 6 to 12 show the differences between the two codes.

Figure 6 :

 6 Figure6: Comparison of the approximate solutions obtained with the method of C (black), and those obtained with TH1D (red). The heating power is here equal to 4.0 10 8 W/m 3 . The volumic steam fraction is equal to zero on the whole domain for both codes.

Figure 7 :

 7 Figure7: Comparison of the approximate solutions obtained with the method of C (black), and those obtained with TH1D (red). The heating power is here equal to 4.5 10 8 W/m 3 .

Figure 8 :

 8 Figure 8: Comparison of the approximate solutions obtained with the method of C (black), and those obtained with TH1D (red). The heating power is here equal to 5.0 10 8 W/m 3 .

Figure 9 :

 9 Figure9: Comparison of the approximate solutions obtained with the method of C (black), and those obtained with TH1D (red). The heating power is here equal to 5.5 10 8 W/m 3 .

Figure 10 :

 10 Figure10: Comparison of the approximate solutions obtained with the method of C (black), and those obtained with TH1D (red). The heating power is here equal to 6.0 10 8 W/m 3 .

Figure 11 :

 11 Figure11: Comparison of the approximate solutions obtained with the method of C (black), and those obtained with TH1D (red). The heating power is here equal to 6.5 10 8 W/m 3 .

Figure 12 :

 12 Figure12: Comparison of the approximate solutions obtained with the method of C (black), and those obtained with TH1D (red). The heating power is here equal to 7.0 10 8 W/m 3 .

Figure 13 :Figure 14 :

 1314 Figure 13: 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 .

Figure 16 :

 16 Figure 16: Pressure at point P1 with respect to the time: two different runs are plotted (black lines with circles and squares), and the numerical results for the code are shown with different laws for the relaxation time-scale λ.

Figure 17 :

 17 Figure 17: vapor volume fraction at point Pt with respect to the time: the results of different runs are gathered (black circles), and the numerical results for the code are shown with different laws for the relaxation time-scale λ.

Figure 18 :

 18 Figure 18: EOS parameters used for the simulation of 4.

Figure 20 :

 20 Figure 20: EOS parameters used for the simulation of 5.

  ) is a non-linear system whose unknowns are Y b , τ b and P b (or e b ). The velocity U b is deduced from the mass flow rate and the specific volume: U b = D 0 τ b .

Figure 21 :

 21 Figure 21: G (top), α b (middle) and P b (bottom) with respect to τ b for the initial conditions of section 5.2 and a heating power equal to 7.0 10 8 W . The red dashed lines correspond to the time-scale λ = 10 -4 s and the black plain line to λ = 10 8 s. A zoom around the zero of G is proposed on figure 22

Figure 22 :

 22 Figure 22: G (top), α b (middle) and P b (bottom) with respect to τ b for the initial conditions of section 5.2 and a heating power equal to 7.0 10 8 W . The red dashed lines correspond to the time-scale λ = 10 -4 s and the black plain line to λ = 10 8 s. Zoom around the zero of G.

i+ 1 2 through the interface {x i+ 1 2 }

 22 × [t n ; t n+1 ] is defined below. The time step must agree with a CFL conditions as detailed below. The flux φ n i+ 1 2

  . The blending function f b (M ) for the pressure is based on the local Mach numberM = |U * (0; Z L , Z R )|/max(C L , C R ),

  α) = α eq (τ, e) -α λ ∂ ∂t (y) = y eq (τ, e) -y λ ∂ ∂t (z) = z eq (τ, e) -

  t)) = α eq (τ (0), e(0)) -α(t) λ(t) ∂ ∂t (y(t)) = y eq (τ (0), e(0)) -y(t) λ(t) ∂ ∂t (z(t)) = z eq (τ (0), e(0)) -z(t) λ

  to evaluate the EOS parameters at a reference pressure equal to 85 bars. The values used in section 6 are gathered in table19.

	P 0 (P a)	85.0 10 5
	T 0 (K)	573.15
	Cv,v	6.0337 10 2
	Qv	2.28830 10 6
	γv	1.49
	Πv	0.0
	s 0 v C v,l	3.951133539778530 10 4 3.61217 10 3
	Q l	-1.5307 10 6
	γ l	1.38
	Π l	5.707983952 10 8
	s 0 l	0
	Figure 19: EOS parameters used for the simulation of 6.

The phasic entropy and the phasic sound speeds are defined for P k > -Π k ; and the phasic temperature is positive forP k > -Π k

Acknowledgements. The present work is supported by the EDF R&D project "Local Thermohydraulic Platform". The numerical results have been obtained with Code Saturne which is developed by EDF R&D and freely distributed under a GNU GPL version 2 licence (http://code-saturne.org). The author thanks G. Baurin, J. Coulet and O. Touazi for providing their results with the TH1D code presented in section 5.

then reads:

α n+1 = e (-∆t/λ n+1, * ) α n+1, * -α n+1, * eq (e (-∆t/λ n+1, * ) -1); y n+1 = e (-∆t/λ n+1, * ) y n+1, * -y n+1, * eq (e (-∆t/λ n+1, * ) -1); z n+1 = e (-∆t/λ n+1, * ) z n+1, * -z n+1, * eq (e (-∆t/λ n+1, * ) -1);

It now remains to compute the equilibrium fractions α eq , y eq and z eq for a given couple (τ, e). The equilibrium fractions Y eq = (α eq , y eq , z eq ) are defined as the fraction Y such that the mixture entropy (3) reaches its maximum value at a given specific energy e and specific volume τ . This maximum is indeed uniquely defined thanks to the concavity of the mixture entropy, see [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF][START_REF] Jung | Numerical simulations of two-fluid flow on multicores accelerator[END_REF]. When the equilibrium fraction lies in ]0, 1[ 3 , the derivative of the mixture entropy with respect to y at Y = Y eq vanishes. This condition can be written:

In order to compute a numerical approximation of Y eq , we choose to solve (54). For Stiffened Gas laws, a method has been proposed in [START_REF] Hurisse | Application of an homogeneous model to simulate the heating of two-phase flows[END_REF]. In this method, two among the fractions are eliminated and the system (54) is solved using one non-linear equation for one fraction. We use here the same algorithm as the algorithm proposed in [START_REF] Faccanoni | Modelling and simulation of liquid-vapor phase transition in compressible flows based on thermodynamical equilibrium[END_REF] which can be applied to all the phasic equations of states. The only requirement is to have the following functions defined on the pressure/temperature thermodynamical plane :

In the pressure/temperature thermodynamical plane, system (54) can be written :

e = y eq e l (P l , T l ) + (1 -y eq )e v (P v , T v ), τ = y eq τ l (P l , T l ) + (1 -y eq )τ v (P v , T v ),

where e and τ are given. We set T l = T v = T and P l = P v = P which leads to:

   e = y eq e l (P, T ) + (1 -y eq )e v (P, T ), τ = y eq τ l (P, T ) + (1 -y eq )τ v (P, T ), µ l (P, T ) = µ v (P, T ).

(57)

For a given pressure P , the last equation of (57) defines the saturation temperature P → T sat (P ) independently of the first and second equations. The specific energy e and the specific volume τ are given, and two liquid mass fractions y l,e and y l,τ can be defined: e = y l,e e l (P, T sat (P )) + (1 -y l,e )e v (P, T sat (P )),

and τ = y l,τ τ l (P, T sat (P )) + (1 -y l,τ )τ v (P, T sat (P )).

These fractions y l,e et y l,τ are functions of the pressure P , and are equal if the pressure P corresponds to the thermodynamical equilibrium. Hence, finding the thermodynamical equilibrium is equivalent to find the zero of the function: P → y l,e (P ) -y l,τ (P ).

(60)

Once the equilibrium pressure P eq found, the equilibrium mass fraction is known and we compute:

α eq = y eq τ l (P eq , T sat (P eq )) τ and z eq = y eq e l (P eq , T sat (P eq )) e .

Once the equilibrium pressure P eq is computed,

• if the fraction y eq is negative, the equilibrium state is the pure phase v, that is y eq = 0;

• if the fraction y eq is greater than 1, the equilibrium state is the pure phase l, that is y eq = 1;

• otherwise the equilibrium state is the mixture of phase l and v defined by y eq .

For the results presented in the paper, the pressure P eq is obtained using a dichotomy algorithm. This algorithm is robust and accurate, and in practice we have already applied it to compute equilibrium fractions using the IAPWS97 functions.