Frédéric Coquel

Jean-Marc Hérard

Khaled Saleh

A Positive and Entropy-Satisfying Finite Volume Scheme for the Baer-Nunziato Model

Keywords: Compressible multi-phase flows, Hyperbolic PDEs, Energy-entropy duality, Entropysatisfying methods, Relaxation techniques, Riemann problem, Riemann solvers, Finite volumes AMS subject classifications : 76T10, 65M08, 35L60, 35F55

. The relaxation scheme also shows a higher precision and a lower computational cost (for comparable accuracy) than a standard numerical scheme used in the nuclear industry, namely Rusanov's scheme. Finally, we assess the good behavior of the scheme when approximating vanishing phase solutions.

Introduction

The modeling and numerical simulation of two-phase flows is a relevant approach for a detailed investigation of some patterns occurring in water-vapor flows such as those encountered in nuclear power plants. The targeted applications are the normal operating mode of pressurized water reactors as well as incidental configurations such as the Departure from Nucleate Boiling (DNB) [START_REF] Nrc | Glossary. Departure from Nucleate Boiling (DNB)[END_REF], the Loss of Coolant Accident (LOCA) [START_REF] Nrc | Glossary. Loss of Coolant Accident (LOCA)[END_REF], the re-flooding phase following a LOCA or the Reactivity Initiated Accident (RIA) [START_REF]Reactivity Initiated Accident (RIA)[END_REF]. In the normal operating mode, the flow in the primary circuit is quasi monophasic as there is a priori no vapor in the fluid. In the incidental configurations however, the vapor statistical fraction may take values ranging from zero to nearly one if some areas of the fluid have reached the boiling point. The modeling as well as the numerical simulation of such phenomena remains challenging since both models that can handle phase transitions and robust numerical schemes are needed. While in the normal operating mode of pressurized water reactors, homogeneous models assuming thermodynamical equilibrium between the phases are used (in practice, only the liquid phase is present), the simulation of incidental configurations requires more detailed two-phase flow models accounting for thermodynamical disequilibrium. Naturally, as opposed to the numerical approximation of homogeneous models, explicit schemes are needed for the simulation of these potentially highly unsteady phenomena, and one major challenge therefore is the control of the numerical time step. In addition, the derived schemes are expected to ensure important stability properties such as the positivity of the densities and internal energies, as well as discrete entropy inequalities. In this context, the aim of this work is to design a robust positivity-preserving and entropy-satisfying scheme for the numerical approximation of two-phase flows with vapor or liquid fractions arbitrarily close to zero.

This paper is concerned with the Baer-Nunziato two-phase flow model introduced in [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF], and studied in various papers [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF][START_REF] Andrianov | The Riemann problem for the Baer-Nunziato two-phase flow model[END_REF][START_REF] Coquel | Closure laws for a two-fluid two pressure model[END_REF][START_REF] Gavrilyuk | Mathematical and numerical modeling of two-phase compressible flows with micro-inertia[END_REF][START_REF] Kapila | Two-phase modeling of DDT: Structure of the velocity-relaxation zone[END_REF]. The model consists of two sets of partial differential equations accounting for the evolution of mass, momentum and total energy for each phase, in addition to a transport equation for the phase fraction. The evolution equations of the two phases are coupled through first order non-conservative terms depending on the phase fraction gradient. A major feature of the Baer-Nunziato model is the assumption of two different velocities and two different pressures for the two phases. This approach is not genuinely usual in the nuclear industry where the commonly implemented methods assume the same pressure for the two phases at every time and everywhere in the flow. This latter assumption is justified by the very short time-scale associated with the relaxation of the phasic pressures towards an equilibrium. In the two-fluid two-pressure models (such as Baer-Nunziato's), zero-th order source terms are added in order to account for this pressure relaxation phenomenon as well as a drag force for the relaxation of the phasic velocities towards an equilibrium. Other source terms can also be included in order to account for the relaxation of phasic temperatures and chemical potentials. However, this work is mainly concerned with the convective effects and these zero-th order relaxation terms are not considered in the present paper. We refer to [START_REF] Coquel | Closure laws for a two-fluid two pressure model[END_REF] for some modeling choices of these zero-th order relaxation terms and to [START_REF] Hérard | A fractional step method to compute a class of compressible gas-liquid flows[END_REF][START_REF] Lochon | Comparison of two-fluid models on steam-water transients[END_REF] for their numerical treatment. We also refer to the Conclusion section 6 for some explanation on how the treatment of these terms will affect the numerical method presented in this paper. Various models exist that are related to the Baer-Nunziato model. One may mention various closure laws for the interfacial velocity and pressure [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF][START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF][START_REF] Glimm | Renormalization group solution of two-phase flow equations for Rayleigh-Taylor mixing[END_REF] or extensions to multi-phase flows [START_REF] Hérard | A three-phase flow model[END_REF][START_REF] Han | Modeling of multi-component flows with phase transition and application to collapsing bubbles[END_REF][START_REF] Müller | Closure conditions for non-equilibrium multi-component models[END_REF].

Various approaches were considered to approximate the admissible weak solutions of the first order Baer-Nunziato model. One may mention exact Riemann solvers [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] or approximate Riemann solvers [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] Ambroso | A Godunov-type method for the seven-equation model of compressible two-phase flow[END_REF]. Let us mention some other schemes grounded on flux of operator splitting techniques [START_REF] Chalons | Large time-step numerical scheme for the seven-equation model of compressible two-phase flows[END_REF][START_REF] Coquel | A splitting method for the isentropic Baer-Nunziato two-phase flow model[END_REF][START_REF] Daude | On the Computation of the Baer-Nunziato Model Using ALE Formulation with HLL-and HLLC-type Solvers Towards Fluid-structure Interactions[END_REF][START_REF] Thanh | Numerical approximation for a Baer-Nunziato model of twophase flows[END_REF][START_REF] Liu | Contribution à la vérification et à la validation d'un modèle diphasique bifluide instationnaire[END_REF][START_REF] Saleh | Analyse et Simulation Numérique par Relaxation d'Ecoulements Diphasiques Compressibles. Contribution au Traitement des Phases Evanescentes[END_REF][START_REF] Thanh | A robust numerical method for approximating solutions of a model of two-phase flows and its properties[END_REF][START_REF] Tokareva | A flux splitting method for the BaerNunziato equations of compressible two-phase flow[END_REF]. Let us also mention the original work of [START_REF] Abgrall | Large time-step numerical scheme for the seven-equation model of compressible two-phase flows[END_REF][START_REF] Dallet | A comparative study of numerical schemes for the Baer-Nunziato model[END_REF] where two staggered grids are used (one for the scalar unknowns and the other for the velocities) and where the internal energies are discretized instead of the total energies.

The finite volume scheme we describe in the present paper for the convective part of the Baer-Nunziato model relies on two main building blocks. The first block is a relaxation finite volume scheme previously designed in [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF] for the isentropic version of the Baer-Nunziato model (the phasic entropies remain constant in both time and space along the process), a scheme which was proved to ensure positive densities and to satisfy discrete energy dissipation inequalities. The second building block is a duality principle between energy and entropy which, according to the second principle of thermodynamics states that, keeping all the other thermodynamic variables constant, the mathematical entropy is a decreasing function of the total energy. This duality principle was already used in previous works to extend schemes designed for the isentropic Euler equations to the full Euler equations (see [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] and [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]), and in this work, we apply these techniques to the Baer-Nunziato two-phase flow model. In [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF], a relaxation Riemann solver was designed for the isentropic Baer-Nunziato model. The main properties of this scheme are firstly, to compute positive densities thanks to an energy dissipation process, secondly to satisfy discrete energy inequalities for each phase, and finally to compute robust approximations of vanishing phase cases where one (or both) of the phase fractions are arbitrarily close to zero in some areas of the flow. The fact that the phasic entropies are simply advected for smooth solutions of the Baer-Nunziato model, combined with the energy-entropy duality principle, actually allows us to use the very same Riemann solver designed in [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF], provided that one supplements it with a correction step which consists in recovering the energy conservation and entropy dissipation for each phase. Concerning the neglected zero-th order source terms, there exist methods that allow their numerical treatment in accordance with the total entropy (the sum of both phasic entropies) dissipation (see [START_REF] Hérard | A fractional step method to compute a class of compressible gas-liquid flows[END_REF][START_REF] Lochon | Comparison of two-fluid models on steam-water transients[END_REF] and the Conclusion section 6).

Nevertheless, we draw the reader's attention on the fact that the relaxation scheme for the isentropic model, and its extension to the full model described here, are restricted to the simulations of flows with subsonic relative speeds, i.e. flows for which the difference between the material velocities of the phases is less than the speed of sound in the dominating phase, which would be the liquid phase in the usual operating of a nuclear power plant. For the simulation of nuclear liquid-vapor flows, this is not a restriction, but it would be interesting though to extend the present scheme to sonic and supersonic flows. An interesting work on this subject is done in [8].

The resulting scheme is proven to preserve positive phase fractions, densities and internal energies, as well as real values of the phasic speeds of sound for stiffened gas and ideal gas e.o.s.. In addition, it is proven to satisfy a discrete entropy inequality for each phase, under a sub-characteristic condition (Whitham's condition). To our knowledge, there exists no other scheme that is proved to satisfy these properties altogether. The relaxation scheme compares well with two of the most popular existing schemes available for the Baer-Nunziato model, namely Schwendeman-Wahle-Kapila's Godunov-type scheme [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] and Toro-Tokareva's HLLC scheme [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]. In addition, for the same level of refinement, the scheme is shown to be much more accurate than Rusanov's scheme, and for a given level of approximation error, the relaxation scheme is shown to perform much better in terms of computational cost than this classical scheme. This is an important result because the approximate Riemann solver designed in [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF] and re-used here relies on a fixed-point research for an increasing scalar function defined on the interval (0, 1). Hence, the numerical tests assess that no heavy computational costs are due to this fixed-point research. Actually, comparing with Rusanov's scheme is quite significant since for such stiff configurations as vanishing phase cases, this scheme is commonly used in the industrial context because of its known robustness and simplicity [START_REF] Hérard | A fractional step method to compute a class of compressible gas-liquid flows[END_REF]. Our relaxation scheme is first-order accurate and an interesting further work is the extension to higher orders (see [START_REF] Dumbser | FORCE schemes on unstructured meshes II: Non-conservative hyperbolic systems[END_REF][START_REF] Franquet | Rungekutta discontinuous galerkin method for the approximation of Baer and Nunziato type multiphase models[END_REF][START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF][START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF] for examples of high order schemes).

The paper is organized as follows. Section 2 is devoted to the presentation of the first order Baer-Nunziato model. In Section 3, an auxiliary two-phase flow model is introduced, where the phasic entropies are conserved and the phasic total energies are dissipated. We explain how to extend the relaxation scheme designed in [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF] to this auxiliary model. For the sake of completeness, the fully detailed Riemann solution is given in Section 7.1 of the appendix. In Section 4, we give the correction step which relies on the energy-entropy duality principle, and the resulting finite volume scheme for the Baer-Nunziato model is fully described. Finally, Section 5 is devoted to the numerical tests. The relaxation finite volume scheme is compared with Schwendeman-Wahle-Kapila's Godunov-type scheme [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF], Toro-Tokareva's HLLC scheme [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF] and Rusanov's scheme. In addition to a convergence and CPU cost study, one test case simulates a near-vacuum configuration, and two test-cases assess that the scheme provides a robust numerical treatment of vanishing phase solutions. For the reader who is eager to rapidly implement the numerical scheme, we refer to Section 7.2 of the Appendix, where the procedure for computing the finite volume numerical fluxes is fully described.

The first order Baer-Nunziato model

The Baer-Nunziato model is a non-viscous two-phase flow model formulated in Eulerian coordinates and describing the evolution of the mass, momentum and total energy of each phase. Each phase is indexed by an integer k ∈ {1, 2}, the density of phase k is denoted ρ k , its velocity u k , and its specific total energy E k . At each point x of the space and at each time t, the probability of finding phase k is denoted α k (x, t). We assume the saturation constraint α 1 + α 2 = 1. In one-space dimension, the convective part of the model introduce in [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] reads:

∂ t U + ∂ x F (U) + C(U)∂ x U = 0, x ∈ R, t > 0, (1)
where

U =           α 1 α 1 ρ 1 α 2 ρ 2 α 1 ρ 1 u 1 α 2 ρ 2 u 2 α 1 ρ 1 E 1 α 2 ρ 2 E 2           , F (U) =           0 α 1 ρ 1 u 1 α 2 ρ 2 u 2 α 1 ρ 1 u 2 1 + α 1 p 1 α 2 ρ 2 u 2 2 + α 2 p 2 α 1 ρ 1 E 1 u 1 + α 1 p 1 u 1 α 2 ρ 2 E 2 u 2 + α 2 p 2 u 2           , C(U)∂ x U =           u 2 0 0 -p 1 +p 1 -p 1 u 2 +p 1 u 2           ∂ x α 1 . (2)
In the complete model, zero-th order source terms are added in order to account for the pressure relaxation phenomenon as well as a drag force for the relaxation of the phasic velocities towards an equilibrium. Other source terms can also be included in order to account for the relaxation of phasic temperatures and chemical potentials. However, this work is mainly concerned with the convective effects and these zero-th order relaxation terms are not considered in the present paper. We refer to the Conclusion section 6 for some explanation on how to treat these relaxation terms without deteriorating the properties of the numerical method presented in this paper.

The state vector U is expected to belong to the natural physical space:

Ω U = U ∈ R 7 , α 1 ∈ (0, 1), α k ρ k > 0, and α k ρ k (E k -u 2 k /2) > 0 for k ∈ {1, 2} . (3)
For each k ∈ {1, 2}, p k denotes the pressure of phase k. Defining e k := E k -u 2 k /2 the specific internal energy of phase k, the pressure p k = p k (ρ k , e k ) is given by an equation of state (e.o.s.) as a function defined for all positive ρ k and all positive e k .

We assume that, taken separately, the two phases follow the second principle of thermodynamics so that for each phase k ∈ {1, 2}, there exists a positive integrating factor T k (ρ k , e k ) such that the following differential form 1

T k p k ρ 2 k dρ k -de k , (4)
is the exact differential of some strictly convex function s k (ρ k , e k ), called the (mathematical) entropy of phase k.

The following proposition characterizes the wave structure of this system: Proposition 2.1. For all U ∈ Ω U , the Jacobian matrix F ′ (U) + C(U) admits the following seven eigenvalues: 1) is weakly hyperbolic on Ω U in the following sense: all the eigenvalues are real and the corresponding right eigenvectors are linearly independent if, and only if,

σ 1 (U) = σ 2 (U) = u 2 , σ 3 (U) = u 1 σ 4 (U) = u 1 -c 1 (ρ 1 , e 1 ), σ 5 (U) = u 1 + c 1 (ρ 1 , e 1 ) σ 6 (U) = u 2 -c 2 (ρ 2 , e 2 ), σ 7 (U) = u 2 + c 2 (ρ 2 , e 2 ), (5)
where c k (ρ k , e k ) 2 = ∂ ρ k p k (ρ k , e k ) + p k (ρ k , e k )/ρ 2 k ∂ e k p k (ρ k , e k ). If c k (ρ k , e k ) 2 > 0, then system (
α 1 = 0, α 2 = 0, |u 1 -u 2 | = c 1 (ρ 1 , e 1 ). (6)
When ( 6) is not satisfied, the system is said to be resonant. The characteristic fields associated with σ 4 , σ 5 , σ 6 and σ 7 are genuinely non-linear, while the characteristic fields associated with σ 1,2 and σ 3 are linearly degenerate.

Remark 2.1. The condition c k (ρ k , e k ) 2 > 0 is a classical condition that ensures the hyperbolicity for monophasic flows. In general, assuming U ∈ Ω U is not sufficient to guarantee that c k (ρ k , e k ) 2 > 0.

For the stiffened gas e.o.s. for instance, where the pressure is given by

p k (ρ k , e k ) = (γ k -1)ρ k e k -γ k p ∞,k , (7)
where γ k > 1 and p ∞,k ≥ 0 are two constants, a classical calculation yields

ρ k c k (ρ k , e k ) 2 = γ k (γ k - 1)(ρ k e k -p ∞,k ).
Hence, the hyperbolicity of the system requires a more restrictive condition than simply the positivity of the internal energy which reads : ρ k e k > p ∞,k . For the stiffened gas e.o.s., the relaxation scheme proposed in this article will be proven to preserve this condition at the discrete level.

Remark 2.2. The system is not hyperbolic in the usual sense because when ( 6) is not satisfied, the right eigenvectors do not span the whole space R 7 . Two possible phenomena may cause a loss of the strict hyperbolicity: an interaction between the advective field of velocity u 2 with one of the acoustic fields of phase 1, and vanishing values of one of the phase fractions α k . In the physical configurations of interest in the present work (such as two-phase flows in nuclear reactors), the flows have strongly subsonic relative velocities, i.e. a relative Mach number much smaller than one:

M = |u 1 -u 2 | c 1 (ρ 1 , e 1 ) << 1, (8)
so that resonant configurations corresponding to wave interaction between acoustic fields and the u 2contact discontinuity are unlikely to occur. In addition, following the definition of the admissible physical space Ω U , one never has α 1 = 0 or α 2 = 0. However, α k = 0 is to be understood in the sense α k → 0 since one aim of this work is to construct a robust enough numerical scheme that could handle all the possible values of α k , k ∈ {1, 2}, especially, arbitrarily small values.

A simple computation shows that the smooth solutions of (1) also obey the following additional conservation laws on the phasic entropies:

∂ t (α k ρ k s k ) + ∂ x (α k ρ k s k u k ) = 0, k ∈ {1, 2}. (9)
As regards the non-smooth weak solutions of (1), one has to add a so-called entropy criterion in order to select the relevant physical solutions. In view of the convexity of the entropy s k (ρ k , e k ), an entropy weak solution is a weak solution of (1) which satisfies the following entropy inequalities in the usual weak sense:

∂ t (α k ρ k s k ) + ∂ x (α k ρ k s k u k ) ≤ 0, k ∈ {1, 2}. (10)
When the solution contains shock waves, inequalities [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] are strict in order to account for the physical loss of entropy due to viscous phenomena that are not modeled in system [START_REF]Atelier de vérification de schémas pour la simulation des modèles diphasiques[END_REF].

The existence of the phasic entropy conservation laws ( 9) and ( 10) will play a central role in the numerical approximation of the solutions of the Baer-Nunziato model. They permit an energyentropy duality principle which allows a natural extension to the non-isentropic model (1) of the energy-dissipative relaxation scheme designed for the isentropic model in [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF].

For the sake of completeness, let us recall the system of PDEs corresponding to the first order isentropic model: for x ∈ R, t > 0:

∂ t α 1 + u 2 ∂ x α 1 = 0, ∂ t (α 1 ρ 1 ) + ∂ x (α 1 ρ 1 u 1 ) = 0, ∂ t (α 1 ρ 1 u 1 ) + ∂ x (α 1 ρ 1 u 2 1 + α 1 p 1 (τ 1 )) -p 1 (τ 1 )∂ x α 1 = 0, ∂ t (α 2 ρ 2 ) + ∂ x (α 2 ρ 2 u 2 ) = 0, ∂ t (α 2 ρ 2 u 2 ) + ∂ x (α 2 ρ 2 u 2 2 + α 2 p 2 (τ 2 )) -p 1 (τ 1 )∂ x α 2 = 0. (11)
In this case, the phasic pressures are functions solely of the phasic specific volumes p k (τ k ), where τ k = ρ -1 k , and the admissible weak solutions are seen to dissipate the phasic energies according to:

∂ t (α k ρ k E k ) + ∂ x (α k ρ k E k u k + α k p k (τ k )u k ) -u 2 p 1 (τ 1 )∂ x α k ≤ 0, k ∈ {1, 2}, (12)
with

E k = u 2 k /2 + e k (τ k
) where e k is an anti-derivative of -p k . In a previous work [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF], a relaxation scheme was designed for this isentropic Baer-Nunziato model. This scheme was proved to satisfy desirable properties such as maintaining positive phase fractions and densities, ensuring discrete counterparts of the energy inequalities [START_REF] Coquel | Closure laws for a two-fluid two pressure model[END_REF], and finally computing with robustness solutions where some phase fractions are arbitrarily close to zero.

Approximating the weak solutions of an auxiliary model

As an intermediate step towards the purpose of approximating the entropy weak solutions of (1), let us introduce the following auxiliary system

∂ t U + ∂ x F(U) + C(U)∂ x U = 0, x ∈ R, t > 0, (13)
where

U =           α 1 α 1 ρ 1 α 2 ρ 2 α 1 ρ 1 u 1 α 2 ρ 2 u 2 α 1 ρ 1 s 1 α 2 ρ 2 s 2           , F(U) =           0 α 1 ρ 1 u 1 α 2 ρ 2 u 2 α 1 ρ 1 u 2 1 + α 1 P 1 α 2 ρ 2 u 2 2 + α 2 P 2 α 1 ρ 1 s 1 u 1 α 2 ρ 2 s 2 u 2           , C(U)∂ x U =           u 2 0 0 -P 1 +P 1 0 0           ∂ x α 1 . (14)
Compared to the classical Baer-Nunziato model (1), the phasic energy equations have been replaced by the two conservation laws for the phasic entropies. Hence, α k ρ k s k now play the role of independent conservative variables whose evolution is governed according to their own conservative equations. The phasic pressures P k are now seen as functions of the phasic specific volumes τ k = ρ -1 k and the phasic entropies s k so that P k = P k (τ k , s k ). These pressure functions are computed as follows: by the second law of thermodynamics, one has:

∂s k ∂e k (ρ k , e k ) = - 1 T k (ρ k , e k ) , with T k (ρ k , e k ) > 0.
Hence, the mapping e → s k (ρ k , e) is monotone and thus invertible for all ρ k > 0. We denote by s → e k (τ k , s) the inverse of this mapping, which is a positive function. The dependency on the density ρ k has been replaced here by a dependency on the specific volume τ k . The pressure function P k (τ k , s k ) is then defined as follows :

P k (τ k , s k ) = p k (τ -1 k , e k (τ k , s k )).
The phasic total energy is recovered by computing E k (u k , τ k , s k ) = u 2 k /2 + e k (τ k , s k ). The auxiliary state vector U is now expected to belong to the physical space:

Ω U = U ∈ R 7 , α 1 ∈ (0, 1), α k ρ k > 0, and α k ρ k e k (τ k , s k ) > 0 for k ∈ {1, 2} . (15)
We have the following property:

Proposition 3.1. The two following equivalent assertions are satisfied :

(i) The mapping

(α k ρ k s k ) : Ω U -→ R U -→ (α k ρ k s k )(U) satisfies ∂ α k ρ k E k (α k ρ k s k )(U) = -1/T k and is convex.
(ii) The mapping

(α k ρ k E k ) : Ω U -→ R + U -→ (α k ρ k E k )(U) satisfies ∂ α k ρ k s k (α k ρ k E k )(U) = -T k and is convex.
Proof. In order to compute the partial derivative of (α k ρ k s k )(U) with respect to α k ρ k E k , let us calculate the differential of α k ρ k s k . Invoking the second law of thermodynamics T k ds k = -de k + p k ρ -2 k dρ k and the definition e k = E k -u 2 k /2 of the internal energy, we obtain:

T k d(α k ρ k s k ) = (α k ρ k )T k ds k + T k s k d(α k ρ k ) = -(α k ρ k )de k + p k ρ -2 k (α k ρ k )dρ k + T k s k d(α k ρ k ) = -(α k ρ k )de k -p k dα k + (p k ρ -1 k + T k s k )d(α k ρ k ) = -(α k ρ k )dE k + (α k ρ k u k )du k -p k dα k + (p k ρ -1 k + T k s k )d(α k ρ k ) = -(α k ρ k )dE k + u k d(α k ρ k u k ) -u 2 k d(α k ρ k ) -p k dα k + (p k ρ -1 k + T k s k )d(α k ρ k ) = -d(α k ρ k E k ) + u k d(α k ρ k u k ) -p k dα k + (E k -u 2 k + p k ρ -1 k + T k s k )d(α k ρ k ). For k ∈ {1, 2}, ∂ α k ρ k E k (α k ρ k s k )(U) is the derivative of α k ρ k s k with respect to α k ρ k E k when keep- ing constant the variables (α i , α i ρ i , α i ρ i u i ) for i ∈ {1, 2}, and the variable α 3-k ρ 3-k E 3-k . Hence, ∂ α k ρ k E k (α k ρ k s k )(U) = -1/T k and the same computation proves that ∂ α k ρ k s k (α k ρ k E k )(U) = -T k .
The proof of the convexity of these two mappings relies on the convexity of the function s k (ρ k , e k ). It follows lengthy calculations (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]). We admit this result.

Of course, smooth solutions of (13) also solve (1) in the classical sense, which implies that they share the same hyperbolic structure, but entropy weak solutions of (1) and ( 13) do differ. Indeed, following Proposition 3.1, since U → (α k ρ k E k )(U) is convex, while the entropy weak solutions of (1) are defined so as to dissipate the phasic entropies, it is natural to select weak solutions of the hyperbolic model [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] according to the differential inequalities:

∂ t (α k ρ k E k ) + ∂ x (α k ρ k E k u k + α k P k (τ k , s k )u k ) -P 1 (τ 1 , s 1 )u 2 ∂ x α k ≤ 0, k ∈ {1, 2}. ( 16
)
Observe that for constant initial entropies s k (x, 0) = s 0 k , the auxiliary model [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] reduces to the isentropic model [START_REF] Chalons | Large time-step numerical scheme for the seven-equation model of compressible two-phase flows[END_REF], with the pressure laws τ k → P(τ k , s 0 k ). Therefore, in the case of constant entropies, extending the relaxation scheme designed in [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF] to the auxiliary model ( 13) is straightforward. Furthermore, even for non constant initial entropies, the derivation of the self-similar solutions for ( 13) is very close to the isentropic setting because the specific entropies are now just advected by the corresponding phase velocity:

∂ t s k + u k ∂ x s k = 0, k ∈ {1, 2}. ( 17
)
For this reason, the Riemann solutions of the auxiliary model [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] are simpler to approximate than those of (1). But again, if smooth solutions of (1) and ( 13) are the same, their shock solutions are distinct. Hence, a numerical scheme for advancing in time discrete solutions of the original PDEs (1)- [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] based on solving a sequence of Riemann solutions for the auxiliary model ( 13)-( 16) must be given a correction which enforces an energy discretization which is consistent with the original model (1), while ensuring discrete entropy inequalities consistently with [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]. The required correction step is in fact immediate because of the general thermodynamic assumptions made on the complete equation of state. It relies on a duality principle in between energy and entropy, which, according to Proposition 3.1, states that

α k ρ k s k is a decreasing function of α k ρ k E k .
In the present section, we provide a relaxation scheme for approximating the energy dissipating weak solutions of the auxiliary system [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF]. This relaxation scheme is straightforwardly obtained as an extension of the relaxation scheme designed in [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF] for the isentropic Baer-Nunziato model and consequently inherits its main properties (positivity of the phase fractions and densities, numerical energy dissipation, robustness for vanishing phase fractions). Again this extension is made possible thanks to the advective equations [START_REF] Dallet | A comparative study of numerical schemes for the Baer-Nunziato model[END_REF] on the entropies. In Sections 3.1 and 3.2, we define the relaxation approximation for system [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] and state the existence theorem for the corresponding Riemann solver. This existence result, as it directly follows from the isentropic case, is not proven here. We refer the reader to [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF] for the complete proof. In Section 3.3, we derive, thanks to this approximate Riemann solver, the numerical scheme for the auxiliary model [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF]. Finally, in Section 4, we explain how to obtain a positive and entropy-satisfying scheme for the original model (1), thanks to the above mentioned duality principle between energy and entropy.

3.1 Relaxation approximation for the auxiliary model [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] System [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] shares the same hyperbolic structure as system [START_REF]Atelier de vérification de schémas pour la simulation des modèles diphasiques[END_REF]. Therefore, it has four genuinely non-linear fields associated with the phasic acoustic waves, which make the construction of an exact Riemann solver very difficult. In the spirit of [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF], the relaxation approximation consists in considering an enlarged system involving two additional unknowns T k , associated with linearizations π k of the phasic pressure laws. This linearization is designed to get a quasilinear enlarged system, shifting the initial non-linearity from the convective part to a stiff relaxation source term. The relaxation approximation is based on the idea that the solutions of the original system are formally recovered as the limit of the solutions of the proposed enlarged system, in the regime of a vanishing relaxation coefficient ε > 0. For a general framework on relaxation schemes we refer to [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF][START_REF] Coquel | Relaxation of fluid systems[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF].

We propose to approximate the Riemann problem for (13) by the self-similar solution of the following Suliciu relaxation model:

∂ t W ε + ∂ x g(W ε ) + d(W ε )∂ x W ε = 1 ε R(W ε ), x ∈ R, t > 0, (18)
with state vector

W = (α 1 , α 1 ρ 1 , α 2 ρ 2 , α 1 ρ 1 u 1 , α 2 ρ 2 u 2 , α 1 ρ 1 s 1 , α 2 ρ 2 s 2 , α 1 ρ 1 T 1 , α 2 ρ 2 T 2 ) T and g(W) =               0 α 1 ρ 1 u 1 α 2 ρ 2 u 2 α 1 ρ 1 u 2 1 + α 1 π 1 α 2 ρ 2 u 2 2 + α 2 π 2 α 1 ρ 1 s 1 u 1 α 2 ρ 2 s 2 u 2 α 1 ρ 1 T 1 u 1 α 2 ρ 2 T 2 u 2               , d(W)∂ x W =               u 2 0 0 -π 1 π 1 0 0 0 0               ∂ x α 1 , R(W) =               0 0 0 0 0 0 0 α 1 ρ 1 (τ 1 -T 1 ) α 2 ρ 2 (τ 2 -T 2 )               . ( 19
)
For each phase k in {1, 2} the pressure π k is a (partially) linearized pressure π k (τ k , T k , s k ), the e.o.s. of which is defined by:

π k (τ k , T k , s k ) = P k (T k , s k ) + a 2 k (T k -τ k ). ( 20
)
In the formal limit ε → 0, the additional variable T k tends towards the specific volume τ k , and the linearized pressure law π k (τ k , T k , s k ) tends towards the original non-linear pressure law P k (τ k , s k ), thus recovering system [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] in the first seven equations of [START_REF] Daude | On the Computation of the Baer-Nunziato Model Using ALE Formulation with HLL-and HLLC-type Solvers Towards Fluid-structure Interactions[END_REF]. The solution of ( 18) should be parametrized by ε. However, in order to ease the notation, we omit the superscript ε in W ε . The constants a k in [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF] are two positive parameters that must be taken large enough so as to satisfy the following sub-characteristic condition (also called Whitham's condition):

a 2 k > -∂ τ k P k (T k , s k ), k in {1, 2}, (21)
for all T k and s k encountered in the solution of [START_REF] Daude | On the Computation of the Baer-Nunziato Model Using ALE Formulation with HLL-and HLLC-type Solvers Towards Fluid-structure Interactions[END_REF]. Performing a Chapman-Enskog expansion, we can see that Whitham's condition expresses that system ( 18) is a viscous perturbation of system [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] in the regime of small ε. In addition, there exists two energy functionals E k (u k , τ k , T k , s k ), which under Whitham's condition, provide an H-theorem like result as stated in Proposition 3.2. The smooth solutions of (18) satisfy the following energy equations

∂ t (α k ρ k E k ) + ∂ x (α k ρ k E k u k + α k π k u k ) -u 2 π 1 ∂ x α k = 1 ε α k ρ k a 2 k + ∂ τ k P k (T k , s k ) (τ k -T k ) 2 , ( 22
)
where

E k := E k (u k , τ k , T k , s k ) = u 2 k 2 + e k (T k , s k ) + π 2 k (τ k , T k , s k ) -P 2 k (T k , s k ) 2a 2 k , k ∈ {1, 2}. ( 23
)
Under Whitham's condition [START_REF] Flåtten | Relaxation two-phase flow models and the subcharacteristic condition[END_REF], to be met for all the (T k , s k ) under consideration, the following Gibbs principles are satisfied for k ∈ {1, 2}:

τ k = arg min T k {E k (u k , τ k , T k , s k )}, and E k (u k , τ k , τ k , s k ) = E k (u k , τ k , s k ), (24)
where

E k (u k , τ k , s k ) = u 2 k /2 + e k (τ k , s k ).
At the numerical level, a fractional step method is commonly used in the implementation of relaxation methods: the first step is a time-advancing step using the solution of the Riemann problem for the convective part of ( 18):

∂ t W + ∂ x g(W) + d(W)∂ x W = 0, (25)
while the second step consists in an instantaneous relaxation towards the equilibrium system by imposing T k = τ k in the solution obtained by the first step. This second step is equivalent to sending ε to 0 instantaneously. As a consequence, we now focus on constructing an exact Riemann solver for the homogeneous convective system [START_REF] Glimm | Renormalization group solution of two-phase flow equations for Rayleigh-Taylor mixing[END_REF]. Let us first state the main mathematical properties of the convective system [START_REF] Glimm | Renormalization group solution of two-phase flow equations for Rayleigh-Taylor mixing[END_REF], the solutions of which are sought in the domain of positive densities ρ k and positive T k :

Ω W = W ∈ R 7 , 0 < α 1 < 1, α k ρ k > 0, α k ρ k T k > 0, for k ∈ {1, 2} . (26)
Proposition 3.3. System (25) is weakly hyperbolic on Ω W in the following sense. For all W ∈ Ω W , the Jacobian matrix g ′ (W) + d(W) admits the following real eigenvalues

σ 1 (W) = σ 2 (W) = σ 3 (W) = u 2 , σ 4 (W) = σ 5 (W) = u 1 , σ 6 (W) = u 1 -a 1 τ 1 , σ 7 (W) = u 1 + a 1 τ 1 , σ 8 (W) = u 2 -a 2 τ 2 , σ 9 (W) = u 2 + a 2 τ 2 . (27)
All the characteristic fields associated with these eigenvalues are linearly degenerate and the corresponding right eigenvectors are linearly independent if, and only if

α 1 = 0, α 2 = 0, |u 1 -u 2 | = a 1 τ 1 . ( 28
)
The smooth solutions of system (25) satisfy the following phasic energy equations:

∂ t (α k ρ k E k ) + ∂ x (α k ρ k E k u k + α k π k u k ) -u 2 π 1 ∂ x α k = 0. ( 29
)
Summing over k ∈ {1, 2}, the smooth solutions are seen to conserve the total mixture energy:

∂ t 2 k=1 α k ρ k E k + ∂ x 2 k=1 (α k ρ k E k u k + α k π k u k ) = 0. ( 30
)
Remark 3.1. In the definition of Ω W , the space of admissible states for the solutions of system (25), no positivity requirement has been given for the phasic specific entropies s k . However, since all the waves are linearly degenerate, the weak solutions are expected to obey a maximum principle on the specific entropies since these two quantities are simply advected:

∂ t s k + u k ∂ x s k = 0, for k = {1, 2}. ( 31
)
Remark 3.2. We look for subsonic solutions which are solutions that remain in the domain of

Ω W where |u 1 -u 2 | < a 1 τ 1 .
Here again, one never has α 1 = 0 or α 2 = 0. However, α k = 0 is to be understood in the sense α k → 0.

Remark 3.3. Since all the characteristic fields of system (25) are linearly degenerate, the mixture energy equation ( 30) is expected to be satisfied for not only smooth but also weak solutions. However, in the stiff cases of vanishing phases where one of the left or right phase fractions α k,L or α k,R is close to zero, ensuring positive values of the densities requires an extra dissipation of the mixture energy by the computed solution (see the comments on Definition 3.1 below).

The relaxation Riemann problem

Let (W L , W R ) be two elements of Ω W . We now consider the Cauchy problem for [START_REF] Glimm | Renormalization group solution of two-phase flow equations for Rayleigh-Taylor mixing[END_REF] with the following Riemann type initial data:

W(x, 0) = W L if x < 0, W R if x > 0. ( 32
)
Extending the relaxation Riemann solution computed in [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF]Section 3] for the isentropic setting to the present Riemann problem ( 25)- [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] follows from the crucial observation that both the relaxation specific volume T k and the specific entropy s k are advected in the same way by the phasic flow velocity u k :

∂ t T k + u k ∂ x T k = 0, ∂ t s k + u k ∂ x s k = 0.
Therefore, for self-similar initial data, the Riemann solution, as soon as it exists, necessarily obeys

T k (ξ) = T k,L , ξ < u * k T k,R , u * k < ξ, s k (ξ) = s k,L , ξ < u * k s k,R , u * k < ξ, (33)
where ξ = x/t is the self-similar variable, and u * k is the effective propagation speed associated with the eigenvalue u k in the Riemann solution. Furthermore, any given combination of these variables, say φ(T k , s k ), is also advected by u k . Hence, we obtain from [START_REF] Kapila | Two-phase modeling of DDT: Structure of the velocity-relaxation zone[END_REF], that the non-linear laws arising from the equation of state evolve in the Riemann solution, virtually the same way as within the isentropic setting. Indeed, the entropies s k in the relaxation model [START_REF] Glimm | Renormalization group solution of two-phase flow equations for Rayleigh-Taylor mixing[END_REF] and in the associated energies [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF], are systematically involved in non-linear functions already depending on the variable T k : namely P k (T k , s k ) and e k (T k , s k ). Such functions are solely evaluated on the left and right states in the selfsimilar initial data and hence always contribute to any given jump conditions in terms of

P k (T k,L , s k,L ), e k (T k,L , s k,L ), P k (T k,R , s k,R ) or e k (T k,R , s k,R
). For instance, computing the value of the linearized pressure π k (ξ) at some point ξ of the Riemann fan goes as follows:

π k (ξ) = π k (τ k (ξ), T k (ξ), s k (ξ)) = P k (T k (ξ), s k (ξ)) + a 2 k (T k (ξ) -τ k (ξ)),
where

P k (T k (ξ), s k (ξ)) = P k (T k,L , s k,L ) if ξ < u * k and P k (T k (ξ), s k (ξ)) = P k (T k,R , s k,R
) otherwise, whereas in the isentropic setting, one would have P k (T k,L ) or P k (T k,R ). The same observations can be made for the internal energy e k (T k , s k ) when computing the total energy

E k (u k (ξ), τ k (ξ), T k (ξ), s k (ξ)).
Hence, compared to the isentropic case, it is just as if the relaxation unknown T k is replaced by a two-dimensional vector (T k , s k ). We formalize these observations in 32) admits a solution if, and only if, the isentropic Riemann problem obtained when taking constant initial entropies s k,L = s k,R , k ∈ {1, 2} while keeping the other initial data unchanged, admits a solution. When such a solution exists, the mathematical formulae for defining the phasic quantities τ k , u k , π k , e k and the void fraction α k within the Riemann fan read exactly the same as in the isentropic framework [16, Section 3], provided the following replacements:

Proposition 3.4. Let (W L , W R ) ∈ Ω W × Ω W . The Riemann problem (25)-(
P k (T k,L ) -→ P k (T k,L , s k,L ), P k (T k,R ) -→ P k (T k,R , s k,R ), e k (T k,L ) -→ e k (T k,L , s k,L ), e k (T k,R ) -→ e k (T k,R , s k,R ). ( 34
)
In the following definition, we recall the main features of a solution to the Riemann problem ( 25)- [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF]. Definition 3.1. Let (W L , W R ) be two states in Ω W . A solution to the Riemann problem ( 25)- [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] with subsonic wave ordering is a self-similar mapping W(x, t) = W r (x/t; W L , W R ) where the function ξ → W r (ξ; W L , W R ) satisfies the following properties:

(i) W r (ξ; W L , W R
) is a piecewise constant function, composed of (at most) seven intermediate states belonging to Ω W , separated by (at most) six contact discontinuities associated with the eigenvalues

u 1 ± a 1 τ 1 , u 2 ± a 2 τ 2 , u 1 , u 2 and such that ξ < min k∈{1,2} {u k,L -a k τ k,L } =⇒ W r (ξ; W L , W R ) = W L , ξ > max k∈{1,2} {u k,R + a k τ k,R } =⇒ W r (ξ; W L , W R ) = W R . (35)
(ii) There exists two real numbers u * 2 and π * 1 (depending on (W L , W R )) such that the function W(x, t) = W r (x/t; W L , W R ) satisfies the following PDEs in the distributional sense: for k ∈ {1, 2},

∂ t α k + u * 2 ∂ x α k = 0, ( 36
) ∂ t (α k ρ k ) + ∂ x (α k ρ k u k ) = 0, ( 37
)
∂ t (α k ρ k u k ) + ∂ x (α k ρ k u 2 k + α k ) -π * 1 ∂ x α k = 0, ( 38
)
∂ t (α k ρ k s k ) + ∂ x (α k ρ k s k u k ) = 0, ( 39
) ∂ t (α k ρ k T k ) + ∂ x (α k ρ k T k u k ) = 0, ( 40
)
where ∂ x α k identifies with the Dirac measure

∆α k δ x-u * 2 t , with ∆α k = α k,R -α k,L . (iii) Furthermore, the function W(x, t) = W r (x/t; W L , W R
) also satisfies the following energy equations in the distributional sense:

∂ t (α 2 ρ 2 E 2 ) + ∂ x (α 2 ρ 2 E 2 u 2 + α 2 π 2 u 2 ) -u * 2 π * 1 ∂ x α 2 = 0, ( 41
) ∂ t (α 1 ρ 1 E 1 ) + ∂ x (α 1 ρ 1 E 1 u 1 + α 1 π 1 u 1 ) -u * 2 π * 1 ∂ x α 1 = -Q(u * 2 , W L , W R )δ x-u * 2 t , (42)
where

Q(u * 2 , W L , W R ) is a non-negative number. (iv)
The solution has a subsonic wave ordering in the following sense:

u 1,L -a 1 τ 1,L < u * 2 < u 1,R + a 1 τ 1,R . (43)
Before stating the existence theorem for subsonic solutions proved in [16, Section 3], let us introduce some notations built on the initial states (W L , W R ) and on the relaxation parameters (a 1 , a 2 ). For k in {1, 2},

u ♯ k := 1 2 (u k,L + u k,R ) - 1 2a k (π k (τ k,R , T k,R , s k,R ) -π k (τ k,L , T k,L , s k,L )) , (44)
π ♯ k := 1 2 (π k (τ k,R , T k,R , s k,R ) + π k (τ k,L , T k,L , s k,L )) - a k 2 (u k,R -u k,L ) , (45)
τ ♯ k,L := τ k,L + 1 a k (u ♯ k -u k,L ), ( 46
) τ ♯ k,R := τ k,R - 1 a k (u ♯ k -u k,R ). ( 47
)
We also introduce the following dimensionless number that only depends on the initial phase fractions:

Λ α := α 2,R -α 2,L α 2,R + α 2,L . (48)
We may now state the existence result for the Riemann problem ( 25)- [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF], which is directly inferred from the existence Theorem for the solutions to the relaxation Riemann problem for the isentropic case designed in [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF]Section 3]. The construction of the self-similar solution is fully provided in Appendix 7.1. Theorem 3.5. Given a pair of admissible initial states

(W L , W R ) ∈ Ω W × Ω W , assume that the parameter a k is such that τ ♯ k,L > 0 and τ ♯ k,R > 0 for k in {1, 2}.
Then there exists solutions with subsonic wave ordering to the Riemann problem (25)- [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF], in the sense of Definition 3.1, if the following condition holds:

(A) -a 1 τ ♯ 1,R < u ♯ 1 -u ♯ 2 -1 a2 Λ α (π ♯ 1 -π ♯ 2 ) 1 + a1 a2 |Λ α | < a 1 τ ♯ 1,L .
Proof. Following Proposition 3.4, see [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF]Section 3] for a constructive proof and the remarks below. See Appendix 7.1 for the expressions of the intermediate states of the solution.

Some comments on Definition 3.1 and Theorem 3.5:

1. Assumption (A) can be explicitly tested in terms of the initial data and the parameters a k , k ∈ {1, 2}. The quantities a 1 τ ♯ 1,L and a 1 τ ♯ 1,R can be seen as two sound propagation speeds, while the quantity

(u ♯ 1 -u ♯ 2 -1 a2 Λ α (π ♯ 1 -π ♯ 2 ))/(1 + a1 a2 |Λ α |)
, which has the dimension of a velocity, measures the difference between the pressures and kinematic velocities of the two phases, in the initial data. Observe that if the initial data is close to the pressure and velocity equilibrium between the two phases, this quantity is expected to be small compared to a 1 τ ♯ 1,L and a 1 τ ♯ 1,R . This is actually the case when, in addition to the convective system (1), zero-th order source terms are added to the model in order to account for relaxation phenomena that tend to bring the two phases towards thermodynamical (T 1 = T 2 ), mechanical (u 1 = u 2 and p 1 = p 2 ) and chemical equilibria (see [START_REF] Coquel | Closure laws for a two-fluid two pressure model[END_REF][START_REF] Flåtten | Relaxation two-phase flow models and the subcharacteristic condition[END_REF] for the models and [START_REF] Hérard | A fractional step method to compute a class of compressible gas-liquid flows[END_REF][START_REF] Lochon | Comparison of two-fluid models on steam-water transients[END_REF] for adapted numerical methods).

The quantity u *

2 is the propagation velocity of the phase fraction wave. It is computed as the zero of a monotone real function z → Ψ (WL,WR) (z) on a bounded interval. Assumption (A) is a sufficient and necessary condition for this function Ψ (WL,WR) to have a unique zero (i.e. a unique number u * 2 satisfying Ψ (WL,WR) (u * 2 ) = 0). Hence, solving this fixed-point problem enables to locate the phase fraction wave by coupling two monophasic systems. Let us stress again on the fact this fixed-point problem is very easy to solve numerically, since it boils down to searching the zero of a strictly monotone function on a bounded interval. We refer the reader to equation (80) in Appendix 7.1 and to the paper [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF] for more details. The phase fraction derivative ∂ x α 1 identifies with the Dirac measure ∆α 1 δ 0 (x-u * 2 t). This means that for all open subset ω ⊂ R and for any self-similar function g(x, t) = g r (x/t), one has :

(ξ,t)∈ω×R + ∂ x α 1 (ξt, t)g(ξt, t)dξ dt = g r (u * 2 ), if u * 2 ∈ ω, 0, otherwise.
Item (ii) implies that, away from the u 2 -wave, the system behaves as two independent relaxation systems, one for each phase.

3. Positivity of phase 1 densities. If the ratio α1,L α1,R is in a neighborhood of 1, the solution computed thanks to condition (A) has positive densities and satisfies the phasic energy equations [START_REF] Hérard | A three-phase flow model[END_REF] in the weak sense. In this case, the solution is said to be energy-preserving and the total mixture energy is also conserved according to the conservative equation [START_REF] Hérard | A fractional step method to compute a class of compressible gas-liquid flows[END_REF]. If α1,L α1,R is too large, or too small, depending on the wave ordering between u * 2 and u * 1 , the solution computed thanks to condition (A) may have non-positive densities in phase 1. In such stiff cases, ensuring positive densities for phase 1 is recovered by allowing a strict dissipation of the phase 1 energy:

∂ t (α 1 ρ 1 E 1 ) + ∂ x (α 1 ρ 1 E 1 + α 1 π 1 ) u 1 -u * 2 π * 1 ∂ x α 1 = -Q(u * 2 , W L , W R )δ x-u * 2 t , (49)
where

Q(u * 2 , W L , W R ) < 0. The function Q(u * 2 , W L , W R
) is a kinetic relation which is chosen large enough so as to impose the positivity of all the phase 1 densities. The value of Q(u * 2 , W L , W R ) parametrizes the whole solution and the choice of Q(u * 2 , W L , W R ) prescribes a unique solution.

4. Positivity of phase 2 densities. Assumption (A) allows to compute the value of the wave propagation velocity u * 2 (see comment 2). With this value, one has to verify that the following property, which is equivalent to the positivity of the phase 2 densities, is satisfied:

(B) u ♯ 2 -a 2 τ ♯ 2,L < u * 2 < u ♯ 2 + a 2 τ ♯ 2,R . (50)
In the numerical applications using this Riemann solver (see Section 5), it will always be possible to ensure property (B) by taking a large enough value of the relaxation parameter a 2 (see Appendix 7.2). Note that this condition is a monophasic condition which is not related to the two-fluid modeling. Indeed, the same condition is required when approximating Euler's equations with a similar relaxation scheme.

5. Maximum principle for the entropies. The phasic entropies s k , k ∈ {1, 2} satisfy a maximum principle in the solution since they are simply advected by the phasic velocities according to [START_REF]Reactivity Initiated Accident (RIA)[END_REF].

6. For the applications envisioned for this work, such as nuclear flows, we are only interested in solutions which have a subsonic wave ordering, i.e. solutions for which the propagation velocity u * 2 of the phase fraction α 1 lies in-between the acoustic waves of phase 1, which is what is required in item (iv). However, the considered solutions are allowed to have phasic supersonic speeds |u k | > a k τ k . Indeed, the subsonic property considered here is related to the relative velocity u 1 -u 2 with respect to the phase 1 speed of sound a 1 τ 1 .

The relaxation scheme for the auxiliary model

In this section, the exact Riemann solver W r (ξ; W L , W R ) for the relaxation system ( 25) is used to derive an approximate Riemann solver of Harten, Lax and van Leer [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] for the simulation of the auxiliary system [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF]. The aim is to approximate the admissible weak solution of a Cauchy problem associated with system (13):

∂ t U + ∂ x F(U) + C(U)∂ x U = 0, x ∈ R, t > 0, U(x, 0) = U 0 (x), x ∈ R, (51)
with a discretization which provides discrete counterparts of the energy inequalities [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF] satisfied by the exact solutions of the auxiliary model. As expected, the numerical scheme is identical to the relaxation scheme designed in [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF] for the isentropic model.

We define a time and space discretization as follows: for simplicity in the notations, we assume constant positive time and space steps ∆t and ∆x, and we define λ = ∆t ∆x . The space is partitioned into cells R = j∈Z C j where

C j = [x j-1 2 , x j+ 1 2 [ with x j+ 1 2 = (j + 1 2
)∆x for all j in Z. The centers of the cells are denoted x j = j∆x for all j in Z. We also introduce the discrete intermediate times t n = n∆t, n ∈ N. The approximate solution at time t n , x ∈ R → U λ (x, t n ) ∈ Ω is a piecewise constant function whose value on each cell C j is a constant value denoted by U n j . Since U λ (x, t n ) is piecewise constant, the exact solution of the following Cauchy problem at time t n

∂ t U + ∂ x F(U) + C(U)∂ x U = 0, x ∈ R, t > 0, U(x, 0) = U λ (x, t n ), x ∈ R, (52)
is obtained by juxtaposing the solutions of the Riemann problems set at each cell interface x j+ 1 2 , provided that these Riemann problems do not interact. The relaxation approximation is an approximate Riemann solver which consists in defining:

U n+1 j := 1 ∆x x j+ 1 2 x j-1 2 U app (x, ∆t) dx, j ∈ Z,
where U app (x, t) is the following approximate solution of (52):

U app (x, t) := j∈Z PW r x -x j+ 1 2 t ; M (U n j ), M (U n j+1 ) 1 [xj,xj+1] (x), (53)
where 1 [xj ,xj+1] is the characteristic function of the interval [x j , x j+1 ] and the mappings P and M are defined by: M : R 7 -→ R 9 (x k ) k=1,..,7 -→ (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 , x 1 , 1 -x 1 ).

(54)

P : R 9 -→ R 7 (x k ) k=1,..,9 -→ (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 ). ( 55
)
For a given vector U, W = M (U) is the relaxation vector obtained by keeping α k , α k ρ k , α k ρ k u k and α k ρ k s k unchanged, while setting T k to be equal to τ k . One says that W ∈ Ω W is at equilibrium if there exists U ∈ Ω U such that W = M (U). For a given relaxation vector W, U = PW is the projection of W which consists in dropping the relaxation unknowns T k .

In order for the interface Riemann problems not to interact and thus for U app (x, t) to be a correct approximate solution of (52) at time ∆t, the time step ∆t is chosen small enough so as to satisfy the CFL condition ∆t ∆x max

k∈{1,2},j∈Z max |(u k -a k τ k ) n j |, |(u k + a k τ k ) n j+1 | < 1 2 . ( 56
)
Of course, at each interface x j+ 1 2 , the relaxation Riemann solver W r ξ; M (U n j ), M (U n j+1 ) depends on two parameters (a k ) n j+ 1 2 , k ∈ {1, 2} which must be chosen so as to ensure the conditions stated in the existence Theorem 3.5, and to satisfy some stability properties. Observe that one might take different relaxation parameters a k , k ∈ {1, 2} for each interface, which amounts to approximating the equilibrium system (13) by a different relaxation approximation at each interface, which is more or less diffusive depending on how large are the local parameters (a k ) n j+ 1 2 , k ∈ {1, 2}. Further discussion on the practical computation of these parameters is postponed to Section 7.2 of the Appendices.

Since W r (ξ; W L , W R ) is the exact solution of the relaxation Riemann problem ( 25)-( 32), the updated unknown U n+1 j may be computed by a non-conservative finite volume formula as stated in Proposition 3.6. Provided the CFL condition (56) is satisfied, the updated unknown U n+1 j is given by:

U n+1 j = U n j - ∆t ∆x F -(U n j , U n j+1 ) -F + (U n j-1 , U n j ) . ( 57
)
where the numerical fluxes read

F -(U L , U R ) = Pg W r 0 -; M (U L ), M (U R ) + PD * (M (U L ), M (U R )) 1 {u * 2 <0} , (58)
F + (U L , U R ) = Pg W r 0 + ; M (U L ), M (U R ) -PD * (M (U L ), M (U R )) 1 {u * 2 >0} , (59)
with D * (W L , W R ) := (α 1,R -α 1,L )(u * 2 (W L , W R ), 0, 0, -π * 1 (W L , W R ), π * 1 (W L , W R ), 0, 0, 0, 0) T . The quantity 1 {u * 2 <0} (resp. 1 {u * 2 >0}
) equals one when u * 2 < 0 (resp. u * 2 > 0) and zero otherwise.

Proof. Under the CFL condition (56), the exact solution of [START_REF] Glimm | Renormalization group solution of two-phase flow equations for Rayleigh-Taylor mixing[END_REF] with the piecewise constant initial data W(x, 0) :

= j∈Z M (U n j )1 [xj ,xj+1] (x) is the function: W(x, t) := j∈Z W r x -x j+ 1 2 t ; M (U n j ), M (U n j+1 ) 1 [xj,xj+1] (x),
since the interface Riemann problems do not interact. In addition, under (56), ( 25) may be written:

∂ t W + ∂ x g(W) + j∈Z D * (M (U n j ), M (U n j+1 ))δ 0 x -x j+ 1 2 -(u * 2 ) n j+ 1 2 t = 0,
where D * (W L , W R ) is defined in the proposition. Integrating this PDE over (x j-1 2 , x j+ 1 2 ) × [0, ∆t] and dividing by ∆x, one obtains:

1 ∆x x j+ 1 2 x j-1 2 W(x, ∆t) dx = M (U n j ) - ∆t ∆x g W r 0 -; M (U n j ), M (U n j+1 ) -g W r 0 + ; M (U n j-1 ), M (U n j ) - ∆t ∆x D * (M (U n j ), M (U n j+1 ))1 (u * 2 ) n j+ 1 2 <0 - ∆t ∆x D * (M (U n j-1 ), M (U n j ))1 (u * 2 ) n j-1 2

>0

.

Applying operator P to this equation yields (57).

This approximate Riemann solver is proved to ensure a conservative discretization of the partial masses, partial entropies and total mixture momentum and to satisfy important stability properties such as the preservation of the densities positivity, a maximum principle for the entropies, and hence the positivity of the phasic internal energies, and discrete energy inequalities which are discrete counterparts of the energy inequalities (29) satisfied by the exact weak solutions of the model. Indeed, we have the following result: Proposition 3.7. The numerical scheme (57) for the auxiliary model has the following properties:

• Positivity: Under the CFL condition (56), the scheme preserves positive values of the phase fractions, densities and internal energies:

for all n ∈ N, if U n j ∈ Ω U for all j ∈ Z, then 0 < (α k ) n+1 j < 1, (α k ρ k ) n+1 j > 0, and (α k ρ k e k (τ k , s k )) n+1 j
> 0 for k = 1, 2 and all j ∈ Z, i.e. U n+1 j ∈ Ω U for all j ∈ Z. Moreover, if the thermodynamics of phase k follows an ideal gas or a stiffened gas e.o.s.(see [START_REF] Archambeau | Code Saturne: A Finite Volume Code for the computation of turbulent incompressible flows -Industrial Applications[END_REF]), then the finite volume scheme (57) preserves positive values of the quantity ρ k e k (τ k , s k ) -p ∞,k :

(ρ k e k (τ k , s k )) n j -p ∞,k > 0, ∀j ∈ Z =⇒ (ρ k e k (τ k , s k )) n+1 j -p ∞,k > 0, ∀j ∈ Z.
• Phasic mass conservation: Denoting F ± i the i th component of vector F ± , the fluxes for the phasic partial masses α k ρ k are conservative:

F - i (U L , U R ) = F + i (U L , U R ) for i in {2, 3}. Hence, denoting (α k ρ k u k ) n j+ 1 2 = F ± 1+k (U n j , U n j+1
) for k = 1, 2, one has:

(α k ρ k ) n+1 j = (α k ρ k ) n j - ∆t ∆x (α k ρ k u k ) n j+ 1 2 -(α k ρ k u k ) n j-1 2 . ( 60
)
• Phasic entropy conservation. The fluxes for the phasic entropies α k ρ k s k are conservative:

F - i (U L , U R ) = F + i (U L , U R ) for i in {6, 7}. Hence, denoting (α k ρ k s k u k ) n j+ 1 2 = F ± 5+k (U n j , U n j+1 ) for k = 1, 2, one has: (α k ρ k s k ) n+1 j = (α k ρ k s k ) n j - ∆t ∆x (α k ρ k s k u k ) n j+ 1 2 -(α k ρ k s k u k ) n j-1 2 . ( 61
)
• Total momentum conservation. The fluxes for the mixture momentum

k=1,2 α k ρ k u k are conservative: k=1,2 F - 3+k (U L , U R ) = k=1,2 F + 3+k (U L , U R ). Hence, denoting ( k=1,2 α k ρ k u 2 k + α k π k ) n j+ 1 2 = k=1,2 F ± 3+k (U n j , U n j+1 ) for k = 1, 2, one has: 2 k=1 (α k ρ k u k ) n+1 j = 2 k=1 (α k ρ k u k ) n j - ∆t ∆x k=1,2 α k ρ k u 2 k + α k π k n j+ 1 2 + ∆t ∆x k=1,2 α k ρ k u 2 k + α k π k n j-1 2 . ( 62
)
• Discrete energy inequalities. Assume that the relaxation parameters (a k ) n j+ 1

2

, k = 1, 2 satisfy Whitham's condition at each time step and each interface, i.e. that for all n ∈ N, j ∈ Z, (a k ) n j+ 1 2 , k = 1, 2 are large enough so that

((a k ) n j+ 1 2 ) 2 > -∂ τ k P k (T k , s k ), (63)
for all T k and s k in the solution ξ → W r ξ; M (U n j ), M (U n j+1 ) . Then, the values U n j , j ∈ Z, n ∈ N, computed by the scheme satisfy the following discrete energy inequalities:

(α k ρ k E k )(U n+1 j ) ≤ (α k ρ k E k )(U n j ) - ∆t ∆x (α k ρ k E k u k + α k π k u k ) n j+ 1 2 -(α k ρ k E k u k + α k π k u k ) n j-1 2 + ∆t ∆x 1 (u * 2 ) n j-1 2 ≥0 (u * 2 π * 1 ) n j-1 2 (α k ) n j -(α k ) n j-1 + ∆t ∆x 1 (u * 2 ) n j+ 1 2 ≤0 (u * 2 π * 1 ) n j+ 1 2 (α k ) n j+1 -(α k ) n j , ( 64
) where for j ∈ Z, (α k ρ k E k u k + α k π k u k ) n j+ 1 2 = (α k ρ k E k u k + α k π k u k ) W r 0 + ; M (U n j ), M (U n j+1
) is the right hand side trace of the phasic energy flux evaluated at x j+ 1 2 . Note that (60) and (61) are updating formulae for the next time step unknown U n+1 j whereas (62) and the energy inequalities (64) are properties satisfied by the values U n j , j ∈ Z, n ∈ N, computed by the numerical scheme.

Proof of Prop. 3.7. The approximate Riemann solver is a Godunov type scheme where U n+1 j is the cell-average over C j of the function U app (x, t). Hence, the positivity property on the phase fractions and phase densities is a direct consequence of Theorem 3.5. For this purpose, energy dissipation [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF] across the u 2 -contact discontinuity may be necessary for enforcing this property when the ratio α1,j α1,j+1

(or its inverse) is large for some j ∈ Z.

The positivity of the phasic internal energies is more intricate. Under the CFL condition (56), (s k ) n+1 j is a convex combination of (s k ) n j-1 , (s k ) n j and (s k ) n j+1 since the phasic entropies are advected in the solutions of the local Riemann problems for [START_REF] Glimm | Renormalization group solution of two-phase flow equations for Rayleigh-Taylor mixing[END_REF]. Let us define jmax ∈ {j -1, j, j + 1} such that (s k ) n jmax = max i=j-1,j,j+1

(s k ) n i . Since s → e k ((τ k ) n+1 j , s
) is a positive decreasing function by the second law of thermodynamics, one has:

e k (τ k ) n+1 j , (s k ) n+1 j ≥ e k (τ k ) n+1 j , (s k ) n jmax > 0.
Note that the quantities e k (τ k ) n+1 j , s are well defined (and positive) since (τ k ) n+1 j > 0. In a similar way, we prove that, if the thermodynamics of phase k follows a stiffened gas e.o.s. according to [START_REF] Archambeau | Code Saturne: A Finite Volume Code for the computation of turbulent incompressible flows -Industrial Applications[END_REF], then the numerical scheme (57) preserves positive values of the quantity ρ k e k (τ k , s k )p ∞,k . It follows from the fact that for a stiffened gas e.o.s., one has:

ρ k e k (τ k , s k ) -p ∞,k = ρ γ k k exp s k,0 -s k C v k ,
where s k,0 is a constant reference entropy and C v k is the (constant) heat capacity at constant volume.

Hence, the function s → ρ k e k (τ k , s) -p ∞,k is also a positive and decreasing function whenever ρ k > 0. This yields:

(ρ k e k (τ k , s k )) n+1 j -p ∞,k ≥ (ρ k ) n+1 j e k (τ k ) n+1 j , (s k ) n jmax -p ∞,k .
The right hand side of this inequality reads:

(ρ k ) n+1 j e k (τ k ) n+1 j , (s k ) n jmax -p ∞,k = ((ρ k ) n+1 j ) γ k exp s k,0 -(s k ) n jmax C v k = (ρ k ) n+1 j (ρ k ) n jmax γ k ((ρ k ) n jmax ) γ k exp s k,0 -(s k ) n jmax C v k = (ρ k ) n+1 j (ρ k ) n jmax γ k (ρ k e k (τ k , s k )) n jmax -p ∞,k
and is therefore positive since at time t n , we have (ρ k e k (τ k , s k )) n j -p ∞,k > 0 for all j ∈ Z. The proof of (60), ( 61) and (62) involves no particular difficulties. It is a direct consequence of equations ( 37), [START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF] and [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] satisfied by the relaxation Riemann solutions at each interface.

Let us prove the discrete energy inequalities (64) satisfied by the scheme under Whitham's condition (63). Assuming the CFL condition (56), the solution of ( 25) over [x j-1 2 , x j+ 1 2 ]×[t n , t n+1 ] is the function

W(x, t) := W r x -x j-1 2 t -t n ; M (U n j-1 ), M (U n j ) 1 [x j-1 2 ,xj] (x) + W r x -x j+ 1 2 t -t n ; M (U n j ), M (U n j+1 ) 1 [xj ,x j+ 1 2 ] (x). ( 65
)
According to Theorem 3.5, this function satisfies the phase 1 energy equation:

∂ t (α 1 ρ 1 E 1 ) + ∂ x (α 1 ρ 1 E 1 u 1 + α 1 π 1 u 1 ) -u * 2 π * 1 ∂ x α 1 = -Q n j-1 2 δ 0 x -x j-1 2 -(u * 2 ) n j-1 2 (t -t n ) -Q n j+ 1 2 δ 0 x -x j+ 1 2 -(u * 2 ) n j+ 1 2 (t -t n ) , (66)
where for i ∈ Z, we have denoted

Q n i-1 2 = Q (u * 2 ) n i-1 2 , M (U n i-1 ), M (U n i ) . Integrating this equation over ]x j-1 2 , x j+ 1 2 [×[t n , t n+1
] and dividing by ∆x yields:

1 ∆x

x j+ 1 2 x j-1 2 (α 1 ρ 1 E 1 )(W(x, t n+1 )) dx ≤ (α 1 ρ 1 E 1 )(M (U n j )) - ∆t ∆x (α 1 ρ 1 E 1 u 1 + α 1 π 1 u 1 ) W r 0 -; M (U n j ), M (U n j+1 ) + ∆t ∆x (α 1 ρ 1 E 1 u 1 + α 1 π 1 u 1 ) W r 0 + ; M (U n j-1 ), M (U n j ) + ∆t ∆x 1 (u * 2 ) n j-1 2 ≥0 (u * 2 π * 1 ) n j-1 2 (α 1 ) n j -(α 1 ) n j-1 + ∆t ∆x 1 (u * 2 ) n j+ 1 2 ≤0 (u * 2 π * 1 ) n j+ 1 2 (α 1 ) n j+1 -(α 1 ) n j , (67) because Q n j-1 2 ≥ 0 and Q n j+ 1 2 ≥ 0. Since the initial data is at equilibrium: W(x, t n ) = M (U n j
) for all x ∈ C j ( i.e. (T 1 ) n j is set to be equal to (τ 1 ) n j ) one has (α 1 ρ 1 E 1 )(M (U n j )) = (α 1 ρ 1 E 1 )(U n j ) according to Proposition 3.2. Applying the Rankine-Hugoniot jump relation to (66) across the line {(x, t), x = x j+ 1 2 , t > 0}, yields:

(α 1 ρ 1 E 1 u 1 + α 1 π 1 u 1 ) W r 0 -; M (U n j ), M (U n j+1 ) = (α 1 ρ 1 E 1 u 1 + α 1 π 1 u 1 ) W r 0 + ; M (U n j ), M (U n j+1 ) + Q n j+ 1 2 1 (u * 2 ) n j+ 1 2 =0 .
Hence, since Q n j+ 1 2 ≥ 0, for the interface x j+ 1 2 , taking the trace of (α 1 ρ 1 E 1 u 1 + α 1 π 1 u 1 ) at 0 + instead of 0 -in (67) only improves the inequality. Furthermore, assuming that the parameter a 1 satisfies Whitham's condition (63), the Gibbs principle stated in [START_REF] Gavrilyuk | Mathematical and numerical modeling of two-phase compressible flows with micro-inertia[END_REF] holds true so that:

1 ∆x x j+ 1 2 x j-1 2 (α 1 ρ 1 E 1 )(U app (x, t n+1 )) dx ≤ 1 ∆x x j+ 1 2 x j-1 2 (α 1 ρ 1 E 1 )(W(x, t n+1 )) dx.
Invoking the convexity of the mapping U → (α 1 ρ 1 E 1 )(U) (see Prop. 3.1), Jensen's inequality implies that

(α 1 ρ 1 E 1 )(U n+1 j ) ≤ 1 ∆x x j+ 1 2 x j-1 2 (α 1 ρ 1 E 1 )(U app (x, t n+1 )) dx,
which yields the desired discrete energy inequality for phase 1. The proof of the discrete energy inequality for phase 2 follows similar steps.

4 A positive and entropy-satisfying scheme for the first order Baer-Nunziato model

In the previous section we have designed a numerical scheme for an auxiliary two-phase flow model where the exact solutions conserve the phasic entropies while the phasic energies are dissipated by shock solutions. The scheme has been proven to satisfy discrete counterparts of these features while ensuring the positivity of the relevant quantities.

In the present section, we describe the correction to be given to the auxiliary scheme in order to conserve the phasic energies while dissipating the phasic entropies. We end up with a numerical scheme which is consistent with the entropy weak solutions of any Cauchy problem associated with the original Baer-Nunziato model (1):

∂ t U + ∂ x F (U) + C(U)∂ x U = 0, x ∈ R, t > 0, U(x, 0) = U 0 (x), x ∈ R. (68)
We keep the same time and space discretization as described in Section 3.3. The approximate solution at time t n , x ∈ R → U λ (x, t n ) ∈ Ω is a piecewise constant function whose value on each cell C j is a constant value denoted by U n j . The updated value U n+1 j is computed through a two-step algorithm described hereunder:

A fractional step algorithm

• Step 1: updating the auxiliary unknown. Given

U n j = (α 1 ) n j , (α 1 ρ 1 ) n j , (α 2 ρ 2 ) n j , (α 1 ρ 1 u 1 ) n j , (α 2 ρ 2 u 2 ) n j , (α 1 ρ 1 E 1 ) n j , (α 2 ρ 2 E 2 ) n j T ,
we begin with setting the auxiliary unknown U n j as follows:

U n j = (α 1 ) n j , (α 1 ρ 1 ) n j , (α 2 ρ 2 ) n j , (α 1 ρ 1 u 1 ) n j , (α 2 ρ 2 u 2 ) n j , (α 1 ρ 1 s 1 )(U n j ), (α 2 ρ 2 s 2 )(U n j ) T ,
where (α k ρ k s k )(U n j ) is the partial entropy of phase k, computed from U n j , knowing the density ρ k , the total energy E k and the kinetic energy u 2 k /2:

(α k ρ k s k )(U n j ) := α k ρ k s k (ρ k ) n j , (E k -u 2 k /2) n j .
Observe that, with this definition of U n j , one has (α

k ρ k E k )(U n j ) = (α k ρ k E k ) n j .
We then compute U n+1,- j by applying the relaxation scheme designed for the auxiliary model:

U n+1,- j = U n j - ∆t ∆x F -(U n j , U n j+1 ) -F + (U n j-1 , U n j ) . (69)
According to Proposition 3.7 the phasic energies are dissipated at the discrete level following:

(α k ρ k E k )(U n+1,- j ) ≤ (α k ρ k E k ) n j - ∆t ∆x (α k ρ k E k u k + α k π k u k ) n j+ 1 2 -(α k ρ k E k u k + α k π k u k ) n j-1 2 + ∆t ∆x 1 (u * 2 ) n j-1 2 ≥0 (u * 2 π * 1 ) n j-1 2 (α k ) n j -(α k ) n j-1 + ∆t ∆x 1 (u * 2 ) n j+ 1 2 ≤0 (u * 2 π * 1 ) n j+ 1 2 (α k ) n j+1 -(α k ) n j .
(70)

• Step 2: Exchanging energy and entropy. This final step is a correction step which aims at enforcing conservative updates for the energies of the original unknown U n+1 j . It simply consists in keeping unchanged the updates of the phase fractions, partial masses and momentum:

(α 1 ) n+1 j := (α 1 ) n+1,- j , (α k ρ k ) n+1 j := (α k ρ k ) n+1,- j , (α k ρ k u k ) n+1 j := (α k ρ k u k ) n+1,- j , k = 1, 2,
(71) while enforcing energy conservation by defining the energies updates as:

(α k ρ k E k ) n+1 j := (α k ρ k E k ) n j - ∆t ∆x (α k ρ k E k u k + α k π k u k ) n j+ 1 2 -(α k ρ k E k u k + α k π k u k ) n j-1 2 + ∆t ∆x 1 (u * 2 ) n j-1 2 ≥0 (u * 2 π * 1 ) n j-1 2 (α k ) n j -(α k ) n j-1 + ∆t ∆x 1 (u * 2 ) n j+ 1 2 ≤0 (u * 2 π * 1 ) n j+ 1 2 (α k ) n j+1 -(α k ) n j . (72)

Finite volume formulation of the scheme

In practice, in the implementation, when performing the first step of the method, i.e. when applying the relaxation scheme to the auxiliary variable U, one does not update the last two variables which are the phasic entropies. Indeed, computing the phasic entropies (α k ρ k s k ) n+1,- j is not needed for the update of the phasic energies which is performed in the second step. Therefore, the two step algorithm described in the previous section can be reformulated as a classical non-conservative finite volume scheme. Indeed, we have the following result: Proposition 4.1. The two step algorithm described in equations (69)-( 71)-( 72) is equivalent to the following non-conservative finite volume scheme:

U n+1 j = U n j - ∆t ∆x F -(U n j , U n j+1 ) -F + (U n j-1 , U n j ) . ( 73
)
where the first five components of F ± (U L , U R ) coincide with the first five components of

F ± (U L , U R ) with U computed from U by imposing α k ρ k s k = (α k ρ k s k )(U).
The last two components of F ± (U L , U R ) are given for k ∈ {1, 2} by:

F - 5+k (U L , U R ) = (α k ρ k E k u k + α k π k u k ) W r 0 + ; M (U L ), M (U R ) + 1 {u * 2 ≤0} (u * 2 π * 1 ) ((α k ) R -(α k ) L ) , F + 5+k (U L , U R ) = (α k ρ k E k u k + α k π k u k ) W r 0 + ; M (U L ), M (U R ) -1 {u * 2 ≥0} (u * 2 π * 1 ) ((α k ) R -(α k ) L ) .
Proof. The proposition follows from elementary verifications using equations ( 69)-( 71)-( 72) and the expressions of the energy numerical fluxes (α

k ρ k E k u k + α k π k u k ) j+ 1 2 given in Proposition 3.7.
For the reader who is eager to rapidly implement the numerical scheme, we refer to Appendix 7.2 where the expressions of the numerical fluxes F ± (U L , U R ) are given in detail.

Recasting the scheme in a finite volume formulation with two interface fluxes F ± (U L , U R ) is very interesting since it allows a nearly straightforward extension of the scheme to the 2D and 3D versions of the Baer-Nunziato model on unstructured meshes. Indeed, the multi-dimensional Baer-Nunziato model is invariant by Galilean transformations. Therefore, by assuming that, in the neighborhood of a multi-D cell interface, one has a local 1D Riemann problem in the orthogonal direction to the interface, it is possible to use the very same fluxes F ± (U L , U R ).

Main properties of the scheme

We may now state the following theorem, which gathers the main properties of this scheme, and which constitutes the main result of the paper.

Theorem 4.2. The finite volume scheme (73) for the Baer-Nunziato model has the following properties:

• Positivity: Under the CFL condition (56), the scheme preserves positive values of the phase fractions, densities and internal energies:

for all n ∈ N, if (U n j ∈ Ω U for all j ∈ Z), then 0 < (α k ) n+1 j < 1, (α k ρ k ) n+1 j > 0, and (E k -u 2 k /2) n+1 j > 0 for k = 1
, 2 and all j ∈ Z, i.e. (U n+1 j ∈ Ω U for all j ∈ Z). Moreover, if the thermodynamics of phase k follows an ideal gas or a stiffened gas e.o.s., then the finite volume scheme (73) preserves real values for the speed of sound of phase k:

for all n ∈ N, if c k (ρ k ) n j , (e k ) n j 2 > 0, for all j ∈ Z, then c k (ρ k ) n+1 j , (e k ) n+1 j 2 > 0,
for all j ∈ Z (see Proposition 2.1 and Remark 2.1 for the definition of c k (ρ k , e k ) 2 ).

• Conservativity: The discretizations of the partial masses α k ρ k , k ∈ {1, 2}, the total mixture momentum α 1 ρ 1 u 1 + α 2 ρ 2 u 2 and the total mixture energy α

1 ρ 1 E 1 + α 2 ρ 2 E 2 , are conservative.
• Discrete entropy inequalities. Assume that the relaxation parameters (a k ) n j+ 1

2

, k = 1, 2 satisfy Whitham's condition at each time step and each interface, i.e. that for all n ∈ N, j ∈ Z, (a k ) n j+ 1 2 , k = 1, 2 are large enough so that

((a k ) n j+ 1 2 ) 2 > -∂ τ k P k (T k , s k ), ( 74
)
for all T k and s k in the solution ξ → W r ξ; M (U n j ), M (U n j+1 ) . Then, the values U n j computed by the scheme satisfy the following discrete entropy inequalities: for k = 1, 2:

(α k ρ k s k )(U n+1 j ) ≤ (α k ρ k s k )(U n j ) - ∆t ∆x (α k ρ k s k u k ) n j+ 1 2 -(α k ρ k s k u k ) n j-1 2 , ( 75
)
where the entropy fluxes

(α k ρ k s k u k ) n j+ 1 2
are defined in Proposition 3.7. These inequalities are discrete counterparts of the entropy inequalities [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] satisfied by the admissible weak solutions of the Baer-Nunziato model [START_REF]Atelier de vérification de schémas pour la simulation des modèles diphasiques[END_REF].

Note that preserving positive values of the phase fractions, the phasic densities and also the phasic internal energies altogether is an unprecedented result. Furthermore, to our knowledge, this scheme is the first scheme approximating the solutions of the Baer-Nunziato model, for which discrete entropy inequalities [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] are proven.

Proof. The positivity of the phase fractions (α k ) n+1 j and partial masses (α k ρ k ) n+1 j follows directly from Proposition 3.7. To check that the proposed algorithm preserves the positivity of the internal energies, namely (e k ) n+1 j = (E k -u 2 k /2) n+1 j > 0, it is sufficient to notice that before the update of the energy in the second step, the internal energy (e k ) n+1,- j := e k ((τ k ) n+1,- j , (s k ) n+1,- j ) is positive according to Proposition 3.7. The second step results in increasing the total energies and one

has (α k ρ k E k ) n+1 j ≥ (α k ρ k E k )(U n+1,- j ) by (70)-(72), while (α k ρ k ) n+1 j = (α k ρ k ) n+1,- j , which yields (E k ) n+1 j ≥ E k (U n+1,- j
). Now, as the kinetic energy ((u k ) n+1,- j ) 2 /2 is unchanged by the second step and E k (U n+1,- j ) = (e k ) n+1,- j -((u k ) n+1,- j ) 2 /2, we infer that (e k ) n+1 j ≥ (e k ) n+1,- j > 0 and hence the required positivity property for the internal energies.

In the same way, for an ideal gas or a stiffened gas e.o.s., we prove that the scheme preserves real values of the speed of sound. One has We now prove the discrete entropy inequalities (75). The first step provides a conservative update of the phasic entropy equations. Indeed, the last two components of the vector equation ( 69) yield

ρ k c k (ρ k , e k ) 2 = γ k (γ k -1)(ρ k e k -p ∞,k ) (with p ∞,k = 0 for an ideal gas). If c k (ρ k ) n j , (e k ) n
(α k ρ k s k ) n+1,- j = (α k ρ k s k )(U n j ) - ∆t ∆x (α k ρ k s k u k ) n j+ 1 2 -(α k ρ k s k u k ) n j-1 2 . ( 76
)
Then, thanks to the thermodynamics assumptions, we know from Proposition 3.1 that

∂ α k ρ k E k (α k ρ k s k )(U) = -1/T k .
Consequently, we infer from (α

k ρ k E k ) n+1 j ≥ (α k ρ k E k )(U n+1,- j ) that (α k ρ k s k )(U n+1 j ) ≤ (α k ρ k s k ) n+1,- j .

Injecting in (76) yields the discrete entropy inequalities (75).

For most of equations of state, that are given as a function p k (ρ k , e k ), the quantity α k ρ k s k cannot be expressed as an explicit function of U, which makes it even impossible to compute the time step initial values of the entropies (α k ρ k s k )(U n j ). Still, this does not prevent the discrete inequalities (75) from holding true.

The impossibility, for many equations of state, to express the entropies s k explicitly does not prevent the computation of the numerical fluxes F ± i (U n j , U n j+1 ), i = 1, .., 7, for the updates of (α 1 ) n+1 j , (α k ρ k ) n+1 j , (α k ρ k u k ) n+1 j in the first step, and the updates of (α k ρ k E k ) n+1 j in the second step. Indeed, even though these numerical fluxes involve terms of the form e k (T k , s k ) and π k (τ k , T k , s k ) to be evaluated on the relaxation Riemann solution of the first step, the discussion of Section 3.2 shows that these functions are solely evaluated on the piecewise constant initial data in terms of e k ((T k ) n j , (s k ) n j ) and P k ((T k ) n j , (s k ) n j ). But observe that by the thermodynamics, e k ((

T k ) n j , (s k ) n j ) is nothing else but (e k ) n j = (E k -u 2 k /2) n j and P k ((T k ) n j , (s k ) n j ) is equal to p k ((ρ k ) n j
, (e k ) n j ). Hence, even though the entropy is used for the analysis of the numerical method, one may implement this scheme even for general and possibly incomplete equations of state that are given as a function p k (ρ k , e k ), since the numerical fluxes may still be computed at each time step in terms of the initial unknowns U n j , j ∈ Z. In particular, this allows the use of tabulated equations of state.

Numerical tests

In this section, we present Riemann-type test-cases on which the performance of the relaxation scheme is tested and compared with that of three other schemes: Schwendeman-Wahle-Kapila's first order Godunov-type scheme [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF], Toro-Tokareva's finite volume HLLC scheme [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF] and Rusanov's scheme (a Lax-Friedrichs type scheme, see [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF]). The thermodynamics follow either an ideal gas or a stiffened gas law:

p k (ρ k , e k ) = (γ k -1)ρ k e k -γ k p ∞,k ,
where γ k > 1 and p ∞,k ≥ 0 are two constants. The e.o.s. parameters of each test-case are given in Table 1 as well as the initial discontinuity position, the final time of the simulation and the CFL number. The initial and intermediate states of each solution are given in Tables 2 to 6. The u 2 -contact discontinuity separates two regions denoted -and + respectively on the left and right sides of the discontinuity. If the u 1 -contact discontinuity has non-zero strength, an additional region L * or R * also exists according to the sign of u 2 -u 1 as described in Figure 1.

We recall that the scheme relies on a relaxation Riemann solver which requires solving a fixed point in order to compute, for every cell interface x j+ 1 2 , the zero of a scalar function (see eq. ( 80) in Appendix 7.2). A dichotomy (bisection) method is used in order to compute this solution. The iterative procedure is stopped when the error is less than 10 In this first test-case, both phases follow an ideal gas e.o.s. (see Table 1). The wave pattern for phase 1 consists of a left-traveling shock, a material contact discontinuity u 1 , a phase fraction discontinuity of velocity u 2 and a right-traveling rarefaction wave. For phase 2 the wave pattern is composed of a left-traveling rarefaction wave, the phase fraction discontinuity, and a right-traveling shock.

x t u 2 u 1 u 1 -c 1 u 2 -c 2 u 1 + c 1 u 2 + c 2 R L + - L * x t u 1 u 2 u 1 -c 1 u 2 -c 2 u 1 + c 1 u 2 + c 2 R L R * + -

Results for Test-case 1

In Figure 2, the approximate solution computed with the relaxation scheme is compared with the exact solution, and with the approximate solutions obtained with the Godunov-type scheme, the HLLC scheme and Rusanov's scheme. The results show that unlike Rusanov's scheme, the three other methods, which give very similar results, correctly capture the intermediate states even for this rather coarse mesh of 100 cells. This coarse mesh is a typical example of an industrial mesh, reduced to one direction, since 100 cells in 1D correspond to a 10 6 -cell mesh in 3D. It appears that the contact discontinuity is captured more sharply by the relaxation scheme, the Godunov-type scheme and the HLLC scheme than by Rusanov's scheme for which the numerical diffusion is larger. However, we observe that the Godunov-type and HLLC schemes seem to be slightly more accurate that the relaxation scheme when capturing the u 2 -contact discontinuity for the phase 1 variables. Indeed, for the relaxation scheme, there is one more point within the contact discontinuity for these variables. We can also see that for the phase 2 variables, there are no oscillations as one can see for Rusanov's scheme: the curves are monotone between the intermediate states. For phase 1, the intermediate states are captured by the relaxation, the Godunov-type and the HLLC methods, while with Rusanov's scheme, this weak level of refinement is clearly not enough to capture any intermediate state. These observations assess that, for the same level of refinement, the relaxation method (as well as the Godunov-type scheme and the HLLC scheme) is much more accurate than Rusanov's scheme.

A mesh refinement process has also been implemented in order to check numerically the convergence of the method, as well as it's performances in terms of CPU-time cost. For this purpose, we compute the discrete L 1 -error between the approximate solution and the exact one at the final time T max = N ∆t, normalized by the discrete L 1 -norm of the exact solution:

E(∆x) = j |φ N j -φ ex (x j , T max )|∆x j |φ ex (x j , T )|∆x , (77)
where φ is any of the non conservative variables (α 1 , ρ 1 , u 1 , p 1 , ρ 2 , u 2 , p 2 ). The calculations have been implemented on several meshes composed of 100 × 2 n cells with n = 0, 1, .., 10 (knowing that the domain size is L = 1). In Figure 3, the error E(∆x) at the final time T max = 0.15, is plotted against ∆x in a log -log scale. We can see that all the errors converge towards zero with the expected order of ∆x 1/2 , except the error for u 2 which seems to converge with a higher rate. However, ∆x 1/2 is only an asymptotic order of convergence, and in this particular case, one would have to implement the calculation on more refined meshes in order to reach the theoretically expected order of ∆x 1/2 . Figure 3 also shows the error on the non conservative variables with respect to the CPU-time of the calculation expressed in seconds for both the relaxation scheme and Rusanov's scheme. Each point of the plot corresponds to one single calculation for a given mesh size. One can see that, if one prescribes a given level of the error, the computational cost of Rusanov's scheme is significantly higher than that of the relaxation method for all the variables. For instance, for the same error on the phase 1 density ρ 1 , the gain in computational cost is more than a hundred times when using the relaxation method rather than Rusanov's scheme which is a quite striking result. Indeed, even if Rusanov's scheme is known for its poor perfrmances in terms of accuracy, it is also an attractive scheme for its reduced complexity. This means that the better accuracy of the relaxation scheme (for a fixed mesh) widely compensates for its (relative) complexity.

Results for Test-case 2

The second test-case was taken from [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]. Phase 1 follows an ideal gas e.o.s. while phase 2 follows a stiffened gas e.o.s. (see Table 1). From left to right, the solution for phase 1 consists of a lefttraveling rarefaction wave, the phase fraction discontinuity, a material contact discontinuity u 1 , and a right-traveling shock. For phase 2 the wave pattern is composed of a left-traveling rarefaction wave, the phase fraction discontinuity, and a right-traveling shock.

As the jump of initial pressures is very large, strong shocks are generated in each phase. The distance between the right shock and contact waves is small in phase 1, which makes it difficult for all the schemes to capture the intermediate states at this weak level of refinement (100 cells). We observe however that the Godunov-type scheme, the HLLC scheme and the relaxation scheme remain more accurate than Rusanov's scheme. We also observe that the narrow intermediate state for ρ 1 between the u 1 -contact discontinuity and the {u 1 + c 1 }-shock is better captured with the Godunov-type and the HLLC schemes than by the relaxation scheme. For phase 2, Rusanov's scheme fails to correctly capture the speed of the right-going shock due to the large difference between the pressures before and after the shock. On the contrary, the other schemes capture the shock with the correct speed.

A convergence study has also been performed for this test-case. The observed convergence rate is slightly larger than ∆x 1/2 , and the error v.s. CPU plots show a smaller computational cost for the relaxation scheme than for Rusanov's scheme. However, on this test-case, the observed gain in the computational time is less than in the first test-case and is not the same for all the variables. This test was also taken from [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]. Both phases consist of two symmetric rarefaction waves and a stationary u 2 -contact discontinuity. As the region between the rarefaction waves is close to vacuum, this test-case is useful to assess the pressure positivity property. Note that the positivity of the pressures is expected here since both phases follow an ideal gas e.o.s.. The results are given in Figure 6. We can see that the computed pressures are positive for all the schemes. In addition, all the schemes have similar results, except for the resolution of the phase fraction discontinuity which appears to be very diffused by Rusanov's scheme while it is exactly captured by the other three schemes. For the relaxation scheme, this is a property satisfied by the relaxation Riemann solver for this type of discontinuities. It is also naturally satisfied by the Godunov-type scheme and by the HLLC scheme (see [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]). We now consider a Riemann problem in which one of the two phases vanishes in one of the initial states, which means that the corresponding phase fraction α 1 or α 2 is equal to zero. For this kind of Riemann problem, the u 2 -contact separates a mixture region where the two phases coexist from a single phase region with the remaining phase. Various examples of such problems were introduced in [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF], [START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF][START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF] or [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF].

Results for Test-case 3

Results for Test-case 4

The solution is composed of a {u 1 -c 1 }-shock wave in the left-hand side (LHS) region where only phase 1 is present. This region is separated by a u 2 -contact discontinuity from the right-hand side (RHS) region where the two phases are mixed. In this RHS region, the solution is composed of a u 1 -contact discontinuity, followed by a {u 2 + c 2 }-rarefaction wave and a {u 1 + c 1 }-shock (see Figure 7). In practice, the numerical method requires values of α 1,L and α 1,R that lie strictly in the interval (0, Therefore, in the numerical implementation, we take α 1,L = 1 -10 -4 . The aim here is to give a qualitative comparison between the numerical approximation and the exact solution. Moreover, there is theoretically no need to specify left initial values for the phase 2 quantities since this phase is not present in the LHS region. For the sake of the numerical simulations however, one must provide such values. We choose to set ρ 2,L , u 2,L and p 2,L to the values on the right of the u 2 -contact discontinuity, which is coherent with the preservation of the Riemann invariants of this wave, and avoids the formation of fictitious acoustic waves for phase 2 in the LHS region. For the relaxation scheme, this choice enables to avoid oscillations of phase 2 quantities in the region where phase 2 is not present. However, some tests have been conducted that assess that taking other values of (ρ 2,L , u 2,L , p 2,L ) has little impact on the phase 1 quantities as well as on the phase 2 quantities where this phase is present.

We can see that for the same level of refinement, the relaxation method, the Godunov-type method and the HLLC method are more accurate than Rusanov's scheme, which can be seen especially for phase 1. As regards the region where phase 2 does not exist, we can see that the three other methods are much more stable than Rusanov's scheme. Indeed, theses schemes behave better than Rusanov's scheme when it comes to divisions by small values of α 2 , since the solution approximated by Rusanov's scheme develops quite large values.

Results for Test-case 5

The last test-case considers the coupling between two pure phases. A left region, where only phase 1 exists (α 1,L = 1), is separated by a u 2 -contact discontinuity from a right region, where only phase 2 is present (α 1,R = 0). In the existence region of phase, the solution is composed of a shock (phase 1) or a rarefaction wave (phase 2).

In the numerical implementation, we set α 1,L = 1 -10 -9 and α 1,R = 10 -9 . In addition, in the LHS region, where phase 2 is absent, we choose to set ρ 2,L , u 1 in the RHS region: we set

ρ 1,R = ρ - 1 , u 1,R = u - 1 and p 1,R = p - 1 .
Another choice could have been made for the initialization of the absent phase by imposing an instantaneous local thermodynamical equilibrium between the phases at time t = 0. This would be coherent with the relaxation zero-th order source terms that are usually added to the model when simulating practical industrial configurations.

The results are given in Figure 8. One can see that, in the LHS region, the quantities of the only present phase 1 are correctly approximated while the quantities of the vanishing phase 2 remain stable despite the division by small values of α 2 . The same observation can be made for the RHS region. On the contrary, Rusanov's scheme fails to approximate such a vanishing phase solution.

Conclusion

The work performed in [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF] and in the present paper provides an accurate and robust finite volume scheme for approximating the entropy dissipating weak solutions of the Baer-Nunziato two-phase flow model. The scheme relies on an exact Riemann solver for a relaxation approximation à la Suliciu of the convective part of the Baer-Nunziato model. To our knowledge, this is the only existing scheme for which the approximated phase fractions, phase densities and phase internal energies are proven to remain positive without any smallness assumption on the data or on the phase fraction gradient. In addition, it is the only scheme for which discrete counterparts of the entropy inequalities satisfied by the exact solutions of the model are proven for all thermodynamically admissible equations of state, under a fully computable CFL condition.

The scheme is well-adapted for subsonic flows (in terms of the relative velocity between the phases) and flows for which the phases are close to the thermodynamical and mechanical equilibrium, a state which is characterized by the equality of pressures, velocities and temperatures of both phases. This numerical method is therefore a natural candidate for simulating the convective part of the complete two-phase flow model, where zero-th order source terms are added to account for the relaxation phenomena that tend to bring the two phases towards thermodynamical, mechanical and chemical equilibria. When these source terms are added, the relaxation scheme can be implemented within a fractional step procedure in order to account for all the physical effects. In such a procedure, the first step is the treatment of the first order part of the Baer-Nunziato model thanks to the relaxation scheme, while the following steps consist in successive ODE solvers where the various relaxation effects are treated separately. To fix ideas, let us recall the general form of the full Baer-Nunziato model with relaxation source terms:

∂ t α k + u 2 ∂ x α k = Φ k , ∂ t (α k ρ k ) + ∂ x (α k ρ k u k ) = Γ k , ∂ t (α k ρ k u k ) + ∂ x (α k ρ k u 2 k + α k p k ) -p 1 ∂ x α k = D k + U Γ k , ∂ t (α k ρ k E k ) + ∂ x (α k ρ k E k u k + α k p k u k ) -u 2 p 1 ∂ x α k = Q k + U D k + H Γ k -p 1 Φ k , (78)
where U = 1 2 (u 1 + u 2 ) and H = 1 2 u 1 u 2 . The quantities Φ k , Γ k , D k and Q k account respectively for the relaxation of pressures, chemical potentials, velocities and temperatures according to:

Φ k = Θ p (p k -p 3-k ), Γ k = Θ µ (µ 3-k -µ k ), D k = Θ u (u 3-k -u k ), Q k = Θ T (T 3-k -T k ).
We refer to [START_REF] Lochon | Comparison of two-fluid models on steam-water transients[END_REF] for the precise definition of the chemical potentials µ k and that of the positive quantities Θ p , Θ µ , Θ u and Θ T . In [START_REF] Lochon | Comparison of two-fluid models on steam-water transients[END_REF], a fractional step method is described for the treatment of these source terms. It is proven that, provided stiffened gas e.o.s. for both phases, every ODE-type step of this method is well posed in the sense that existence and uniqueness of the considered quantities are guaranteed in the relevant intervals. Moreover, at the semi-discrete level in time, each one of these steps is compatible with the total entropy inequality satisfied by (78):

∂ t k=1,2 α k ρ k s k + ∂ x k=1,2 α k ρ k s k u k ≤ - Θ p T 2 (p 1 -p 2 ) 2 -Θ µ (µ 1 -µ 2 ) 2 -Θ u T 1 + T 2 2T 1 T 2 (u 1 -u 2 ) 2 - Θ T T 1 T 2 (T 1 -T 2 ) 2 .
The relaxation scheme was specially designed for the simulation of vanishing phase solutions, where in some areas of the flow, the fluid is quasi monophasic i.e. one of the phases has nearly disappeared. In particular, the scheme has been proven to robustly handle sharp interfaces between two quasimonophasic regions as assessed by the results of Test-case 5 (see Section 5.5). Simulating vanishing phase solutions is a crucial issue for a detailed investigation of incidental configurations in the nuclear industry such as the Departure from Nucleate Boiling (DNB) [START_REF] Nrc | Glossary. Departure from Nucleate Boiling (DNB)[END_REF], the Loss of Coolant Accident (LOCA) [START_REF] Nrc | Glossary. Loss of Coolant Accident (LOCA)[END_REF] or the Reactivity Initiated Accident (RIA) [START_REF]Reactivity Initiated Accident (RIA)[END_REF]. Some numerical methods had already been proposed for the approximation of vanishing phase solutions ( [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF][START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]). The work performed in [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF] and in the present paper provides a detailed theoretical and numerical answer to the robustness issues rising up when attempting to simulate vanishing phase solutions.

Despite a relatively complex theory aiming at constructing the underlying approximate Riemann solver, and at analyzing the main properties of the numerical method (positivity, discrete entropy inequalities), the proposed scheme is a rather simple scheme as regards its practical implementation as explained in the appendices 7. The scheme applies for all equations of state for which the pressure is a given function of the density and of the specific internal energy. In particular, this allows the use of incomplete or tabulated equations of state.

It appears that the relaxation scheme has similar performances as two of the most popular existing schemes available for the Baer-Nunziato model, namely Schwendeman-Wahle-Kapila's first order accurate Godunov-type scheme [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] and Toro-Tokareva's finite volume HLLC scheme [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF]. In addition, the scheme compares very favorably with Lax-Friedrichs type schemes that are commonly used in the nuclear industry for their known robustness and simplicity. As a matter of fact, the relaxation finite volume scheme was proved to be much more accurate than Rusanov's scheme for the same level of refinement. In addition, for a prescribed level of accuracy (in terms of the L 1 -error for instance), the computational cost of the relaxation scheme is much lower than that of Rusanov's scheme. Indeed, for some test-cases, reaching the same level of accuracy on some variables may require more than a hundred times more CPU-time to Rusanov's scheme than to the relaxation scheme! In a recent benchmark on numerical methods for two-phase flows [START_REF]Atelier de vérification de schémas pour la simulation des modèles diphasiques[END_REF], the relaxation scheme was proven to compare very well with various other schemes in terms of CPU-time performances as well as robustness [START_REF] Dallet | A comparative study of numerical schemes for the Baer-Nunziato model[END_REF].

Thanks to the invariance of the Baer-Nunziato model under Galilean transformations, the finite volume formulation of the relaxation scheme allows a straightforward extension to 2D and 3D unstructured meshes. As a matter of fact, the scheme has already been implemented in a proprietary module for 3D two-phase flows developed by the French national electricity company EDF within the framework of the industrial CFD code Code Saturne [START_REF] Archambeau | Code Saturne: A Finite Volume Code for the computation of turbulent incompressible flows -Industrial Applications[END_REF]. The scheme has been successfully applied within nuclear safety studies, for numerical simulations of the primary circuit of pressurized water reactors. A forthcoming paper is in preparation, where the relaxation scheme is used for the simulation of 3D industrial cases.

Wave structure Figure 3: Test-case 1: L 1 -Error with respect to ∆x for the relaxation scheme and L 1 -Error with respect to computational cost (in seconds) for the relaxation scheme (straight line) and Rusanov's scheme (dashed line).

α1 x t u 2 u 2 + c 2 u 1 -c 1 u 1 u 1 + c 1 u 2 -c 2 0.
Wave structure Figure 5: Test-case 2: L 1 -Error with respect to ∆x for the relaxation scheme and L 1 -Error with respect to computational cost (in seconds) for the relaxation scheme (straight line) and Rusanov's scheme (dashed line).

α1 x t u 1 + c 1 u 2 + c 2 u 2 u 1 u 2 -c 2 u 1 -
Wave structure 3. Use an iterative method (e.g. Newton's method or a dichotomy (bisection) method) to compute M * L ∈ (0, 1) such that

α 1 x t u 2 u 1 -c 1 u 2 -c 2 u 1 + c 1 u 2 + c 2 0.
Ψ(M * L ) = M ♯ L - a 1 a 2 Λ α P ♯ L . (80)
According to [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF], M * L always exists under (A1) and is unique if µ is close enough to one. In practice, the iterative method is initialized at m 0 = max(0, min(M ♯ L , 1)). The non-conservative product d(W)∂ x W:

When α 1,L = α 1,R , the non-conservative product d(W)∂ x W identifies with a Dirac measure propagating at the constant velocity u * 2 . This Dirac measure is given by

D * (W L , W R )δ x-u * 2 t ,
where D * (W L , W R ) := (α 1,R -α 1,L ) (u * 2 , 0, 0, -π * 1 , +π * 1 , 0, 0, 0, 0) T . The pressure π * 1 is defined for α 1,R = α 1,L by

π * 1 := π ♯ 2 -a 2 α 2,R + α 2,L α 1,R -α 1,L (u * 2 -u ♯ 2 ).

Practical implementation of the relaxation finite volume scheme

In this appendix, we describe in detail the practical implementation of the scheme. We recall the space and time discretization: we assume a positive space step ∆x and the time step ∆t is dynamically updated through the CFL condition. The space is partitioned into cells R = j∈Z [x j-1 2 , x j+ 1 2 [ with x j+ 1 2 = (j + 1 2 )∆x for all j in Z. The centers of the cells are denoted x j = j∆x for all j in Z. We also introduce the discrete intermediate times t n = n∆t, n ∈ N.

The solution of the Cauchy problem:

∂ t U + ∂ x F (U) + C(U)∂ x U = 0, x ∈ R, t > 0, U(x, 0) = U 0 (x),
x ∈ R, is approximated at time t n by U n j on the cell [x j-1 2 , x j+ 1 2 [. The values of the approximate solution are inductively computed as follows:

Initialization:

U 0 j = 1 ∆x
x j+ 1 2

x j-1 2 U 0 (x) dx.

Time evolution:

U n+1 j = U n j - ∆t ∆x F -(U n j , U n j+1 ) -F + (U n j-1 , U n j ) . (83)
At each cell interface x j+ 1 2 , the numerical fluxes F ± (U n j , U n j+1 ) are computed thanks to the relaxation approximate Riemann solver. They depend on the states (U n j , U n j+1 ) but also on the local values of the relaxation parameters a n k,j+ 1 2 , k = 1, 2. Denoting U L = U n j and U R = U n j+1 and a k , k = 1, 2 for simplicity, the fluxes F ± (U L , U R ) are computed through the following steps.

1. Local choice of the pair (a 1 , a 2 ). The pair of parameters (a 1 , a 2 ), must be chosen large enough so as to satisfy several requirements:

• In order to ensure the stability of the relaxation approximation, a k must satisfy Whitham's condition (63). For simplicity however, we do not impose Whitham's condition everywhere in the solution of the Riemann problem ( 25)- [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] (which is possible however), but only for the left and right initial data at each interface:

for k in {1, 2}, a k > max (ρ k,L c k (ρ k,L , e k,L ), ρ k,R c k (ρ k,L , e k,L )) ,

where c k (ρ k , e k ) is the speed of sound in phase k. In practice, no instabilities were observed during the numerical simulations due to this simpler Whitham-like condition. • In order to compute the solution of the relaxation Riemann problem, the specific volumes τ ♯ k,L (U L , U R ; a k ) and τ ♯ k,R (U L , U R ; a k ) defined in (79) must be positive. The expressions of τ ♯ k,L (U L , U R ; a k ) and τ ♯ k,R (U L , U R ; a k ) are two second order polynomials in a -1 k whose constant terms are respectively τ k,L and τ k,R . Hence, by taking a k large enough, one can guarantee that τ ♯ k,L (U L , U R ; a k ) > 0 and τ ♯ k,R (U L , U R ; a k ) > 0, since the initial specific volumes τ k,L and τ k,R are positive.

• Finally, in order for the relaxation Riemann problem ( 25)- [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] to have a positive solution, (a 1 , a 2 ) must be chosen so as to meet condition (A) of Theorem 3.5 as well as the positivity condition of the phase 2 densities (B) (see the comments after Theorem 3.5).

Thereafter, we propose an iterative algorithm for the computation of the parameters (a 1 , a 2 ) at each interface. The notation not(P) is the negation of the logical statement P.

• Choose η a (small) parameter in the interval (0, 1).

• For k in {1, 2} initialize a k : In this algorithm, the computation of u * 2 requires the computation of the solution of the fixedpoint problem (80), using some numerical method such as Newton's method or a dichotomy (bisection) algorithm. It is possible to prove that this algorithm always converges in the sense that there is no infinite looping due to the while-conditions. Indeed, it is easy to observe that assumptions (A) and (B) are always satisfied if the parameters (a 1 , a 2 ) are taken large enough. Moreover, this algorithm provides reasonable values of a 1 and a 2 , since in all the numerical simulations, the time step obtained through the CFL condition (56) remains reasonably large and does not go to zero. In fact, the obtained values of a 1 and a 2 are quite satisfying since the relaxation scheme compares very favorably with Rusanov's scheme, in terms of CPU-time performances (see Section 5).

2. Calculation of the numerical fluxes. Once the relaxation parameters are known, one may give the expressions of the numerical fluxes F ± (U L , U R ). Observe that, as a by-product of the above algorithm for the computation of (a 1 , a 2 ), the propagation velocity u * 2 is already known, and one does not need to redo the fixed-point procedure. Given the solution ξ → W(ξ; U L , U R ; a 1 , a 2 ) of the relaxation Riemann problem ( 25)- [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] (see appendix 7.1 for the expression of the intermediate states), which we denote W(ξ) for the sake of simplicity, the numerical fluxes are computed as follows:

F ± (U L , U R ) =                     0 (α 1 ρ 1 u 1 ) W(0 ± ) (α 2 ρ 2 u 2 ) W(0 ± ) (α 1 ρ 1 u 2 1 + α 1 π 1 ) W(0 ± ) (α 2 ρ 2 u 2 2 + α 2 π 2 ) W(0 ± ) (α 1 ρ 1 E 1 u 1 + α 1 π 1 u 1 ) W(0 + ) (α 2 ρ 2 E 2 u 2 + α 2 π 2 u 2 ) W(0 + )                     +                     (u * 2 ) ± 0 0 - (u * 2 ) ± u * 2 π * 1 (u * 2 ) ± u * 2 π * 1 -(u * 2 ) ± π * 1 (u * 2 ) ± π * 1                     (α 1,R -α 1,L ),
where u * 2 is already known as a result of the first step (choice of the pair (a 1 , a 2 )) and the expression of π *

1 is given at the end of appendix 7.1. In the above expression of the numerical fluxes, we have denoted (u *

2 ) + = max(u * 2 , 0), (u * 2 ) -= min(u * 2 , 0) and the functions x → (x) ± x are extended by 0 at x = 0. Finally, the time step is computed so as to satisfy the CFL condition:

∆t ∆x max k∈{1,2},j∈Z max |(u k -a k τ k ) n j |, |(u k + a k τ k ) n j+1 | < 1 2 ,
and the scheme (83) can be now applied to update the values of the unknown U n+1 j for j ∈ Z.

j 2 >

 2 0, then, before the update of the energy in the second step, the quantity (ρ k e k -p ∞,k ) n+1,- j is positive by Proposition 3.7. The update of the energy in the second step amounts to increasing the internal energy (e k ) n+1 j > (e k ) n+1,- j while keeping the density unchanged (ρ k ) n+1 j = (ρ k ) n+1,- j , which yields (ρ k e k -p ∞,k ) n+1 j > 0, hence the positivity of c k (ρ k , e k ) 2 .

Figure 1 :

 1 Figure 1: Structure of the Riemann solutions, notations for the intermediate states.

Figure 6 :

 6 Figure 6: Test-case 3: Structure of the solution and space variations of the physical variables at the final time T max = 0.15. Mesh size: 100 cells.

Figure 8 :

 8 Figure 8: Test-case 5: Structure of the solution and space variations of the physical variables at the final time T max = 0.05. Mesh size: 1000 cells. Straight line: exact solution, circles: relaxation scheme.

Figure 9 :a 1 a 2 α 1

 921 Figure 9: Intermediate states of the exact solution of the Riemann problem (25)-[START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF] with the wave ordering u * 2 < u * 1 .

4 . 1 -a 1 τ ♯ 1 ,L M * L . 5 . 2 + νa 1 τ ♯ 1 6 .

 4115216 The velocity u * 2 is obtained by u * 2 = u ♯ The velocity u * 1 is obtained by u * 1 = u * The intermediate states for phase 1 are given by• Phase fractions: α - 1 = α 1,L , α + 1 = α 1,R * = α 1,R . • Specific volumes: u - 1 = u * 2 + a 1 τ ♯ 1,L M(M * L ) 1 -M * L 1 -M(M * L ) , u + 1 = u 1,R * = u * 1 .

Figure 10 :

 10 Figure 10: Intermediate states of the exact solution of the Riemann problem (25)-(32) with the wave ordering u * 2 = u * 1 .

a•

  k := (1 + η) max (ρ k,L c k (ρ k,L , e k,L ), ρ k,R c k (ρ k,R , e k,R )). For k in {1, 2}: do {a k := (1 + η)a k } while τ ♯ k,L (U L , U R ; a k ) ≤ 0 or τ ♯ k,R (U L , U R ; a k ) ≤ 0 . • do { a 2 := (1 + η)a 2 ,do {a 1 := (1 + η)a 1 } while not(A) , compute the value of u * 2 in the solution W(U L , U R ; a 1 , a 2 ), do } while not(B) .

Table 1 :

 1 -12 . E.O.S. parameters, initial discontinuity position, final time, Courant-Friedrichs-Lewy number.

		Test 1	Test 2	Test 3	Test 4	Test 5
	γ 1	1.4	1.4	1.4	3	3
	p ∞,1	0	0	0	0	0
	γ 2	1.4	3	1.4	1.4	1.4
	p ∞,2	0	100	0	0	0
	x 0	0	0.8	0.5	0	0
	T max	0.15	0.007	0.15	0.15	0.05
	CFL	0.45	0.45	0.45	0.45	0.45

Table 2 :

 2 Test-case 1: Left, right and intermediate states of the exact solution.

		Region L	Region L *	Region -	Region +	Region R
	α 1	0.2	0.2	0.2	0.7	0.7
	ρ 1	0.21430	0.35	0.698	0.90583	0.96964
	u 1	-0.02609	-0.7683	-0.7683	-0.11581	-0.03629
	p 1	0.3	0.6045	0.6045	0.87069	0.95776
	ρ 2	1.00003	1.00003	0.9436	1.0591	0.99993
	u 2	0.00007	0.00007	0.0684	0.0684	-0.00004
	p 2	1.0	0.9219	0.9219	1.08383	1.0

Table 3 :

 3 Test-case 2: Left, right and intermediate states of the exact solution.

		Region L	Region -	Region +	Region R *	Region R
	α 1	0.3	0.3	0.8	0.8	0.8
	ρ 1	1.0	0.4684	0.50297	5.9991	1.0
	u 1	-19.59741	6.7332	-1.75405	-1.75405	-19.59741
	p 1	1000.0	345.8279	382.08567	382.08567	0.01
	ρ 2	1.0	0.7687	1.6087	1.6087	1.0
	u 2	-19.59716	-6.3085	-6.3085	-6.3085	-19.59741
	p 2	1000.0	399.5878	466.72591	466.72591	0.01

Table 4 :

 4 Test-case 3: Left, right and intermediate states of the exact solution.

		Region L	Region -	Region +	Region R
	α 1	0.2	0.2	0.5	0.5
	ρ 1	0.99988	0.0219	0.0219	0.99988
	u 1	-1.99931	0.0	0.0	1.99931
	p 1	0.4	0.0019	0.0019	0.4
	ρ 2	0.99988	0.0219	0.0219	0.99988
	u 2	-1.99931	0.0	0.0	1.99931
	p 2	0.4	0.0019	0.0019	0.4

Table 5 :

 5 Test-case 4: Left, right and intermediate states of the exact solution.

		Region L	Region -	Region +	Region R *	Region R
	α 1	1.0	1.0	0.4	0.4	0.4
	ρ 1	1.6	2.0	1.84850	2.03335	1.62668
	u 1	0.80311	0.4	0.91147	0.91147	0.55623
	p 1	1.3	2.6	2.05277	2.05277	1.02638
	ρ 2	-	-	4.0	4.0	7.69667
	u 2	-	-	0.1	0.1	0.74797
	p 2	-	-	2.45335	2.45335	6.13338

Table 6 :

 6 2,L and p 2,L to the values on the right of the u 2 -contact discontinuity i.e. to the values ρ + 2 , u + 2 and p + 2 . The symmetric choice is made for phase Test-case 5: Left, right and intermediate states of the exact solution.

		Region L	Region -	Region +	Region R
	α 1	1.0	1.0	0.0	0.0
	ρ 1	1.6	2.0	-	-
	u 1	1.79057	1.0	-	-
	p 1	5.0	10.0	-	-
	ρ 2	-	-	2.0	2.67183
	u 2	-	-	1.0	1.78888
	p 2	-	-	10.0	15.0

  c 1

												α 1
		10 -1
									0.9 10 -2		Rusanov HLLC
	E(∆x)	10 -3 10 -2					∆x 1/2 α 1 ρ 1 u 1 p 1 ρ 2 u 2 p 2		0.5 0.6 0.7 0.8 10 -3		Godunov-type Relaxation Exact
		10 -4	10 -5		10 -4	∆x	10 -3	10 -2	0.3 0.4 10 0		10 1	10 2	10 3	10 4	10 5
									0.2
					u 1			0		0.2	0.4 u 2	0.6	0.8	1
	10 -1			u1							u2
		10	Rusanov							Rusanov
		5	Relaxation HLLC Godunov-type					-5 10 -2	Relaxation HLLC Godunov-type
				Exact							Exact
	0 10 -2							-10
		-5							10 -3
	-10							-15
	-15 10 -3
	-20 10 0		10 1	10 2	10 3	10 4	10 5	-20 10 -4 10 0		10 1	10 2	10 3	10 4	10 5
		0		0.2	0.4	0.6	0.8	1	0		0.2	0.4	0.6	0.8	1
					ρ1 ρ 1						ρ2 ρ 2
	6 10 -1	Rusanov HLLC Godunov-type					1.6		Rusanov HLLC Godunov-type
		5	Relaxation Exact					1.4		Relaxation Exact
		4							10 -2
									1.2
	3 10 -2
									1
		2
		1							0.8 10 -3
	0 10 -3 0 10 0		0.2 10 1	0.4 10 2	0.6 10 3	0.8 10 4	1 10 5	0.6 10 0 0		0.2 10 1	0.4 10 2	0.6 10 3	0.8 10 4	1 10 5
		1200			p1 p 1			1200			p2 p 2
		1000							1000
	800 10 -2							800 10 -2
		600							600
		400							400
	200 10 -3		Rusanov HLLC					200 10 -3		HLLC Rusanov
			Godunov-type							Godunov-type
		0	Relaxation Exact					0		Relaxation Exact
	10 -4 10 0	0	0.2 10 1	0.4 10 2	0.6 10 3	0.8 10 4	1 10 5	10 -4 10 0	0	0.2 10 1	0.4 10 2	0.6 10 3	0.8 10 4	1 10 5
		Figure 4: Test-case 2: Structure of the solution and space variations of the physical
		variables at the final time T max = 0.007. Mesh size: 100 cells.

  Structure of the solution and space variations of the physical variables at the final time T max = 0.15. Mesh size: 100 cells.

	u 1 ρ 1 p 1 u 2 0 0 0 u 1 0 ρ 1 0 p 1 0 Figure 7: Test-case 4: Wave structure -0.4 -0.2 No phase 2 Godunov-type 0.2 Rusanov HLLC Relaxation Exact -0.4 -0.2 0.2 No phase 2 Rusanov HLLC Godunov-type Relaxation Exact -0.4 -0.2 0.2 No phase 2 Rusanov HLLC Godunov-type Relaxation Exact t u 1 -c 1 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 only phase 1 only phase 2 0.4 0.4 0.4 u 2 + c 2 x 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 -0.4 -0.2 0.2 0.4 No phase 2 No phase 1 1.6 1.7 1.8 1.9 2 2.1 -0.4 -0.2 0.2 0.4 No phase 2 No phase 1 4 5 6 7 8 9 11 -0.4 -0.2 0.2 0.4 No phase 2 No phase 1 10	1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 0 0.2 0.4 0.6 0.8 1 1 1.2 1.4 1.6 1.8 2 1.8 2 2.2 2.4 2.6 2.8 3 9 10 11 12 13 14 15 17 16	-0.4 -0.4 -0.4 No phase 2 -0.2 No phase 2 -0.2 No phase 2 -0.2 -0.4 -0.2 No phase 2 -0.4 -0.2 No phase 2 -0.4 -0.2 No phase 2 -0.4 -0.2 No phase 2 -0.4 -0.2 No phase 2	0 0 0 0 0 0 u 2 ρ 2 0 p 2 0 α 1 u 2 ρ 2 p 2	0.2 Rusanov HLLC Godunov-type Relaxation Exact 0.2 Rusanov HLLC Godunov-type Relaxation Exact 0.2 Rusanov HLLC Godunov-type Relaxation Exact 0.2 Rusanov HLLC Godunov-type Relaxation Exact 0.2 No phase 1 0.4 0.4 0.4 0.4 0.4 0.2 0.4 No phase 1 0.2 0.4 No phase 1 0.2 0.4 No phase 1

Acknowledgements. The authors would like to thank Nicolas Seguin who is a co-author of [16] upon which the present paper is based. The authors would like to thank him for his thorough reading of this paper which has constantly led to improving the text. The authors are grateful to Eleuterio Toro and Svetlana Tokareva who kindly agreed to pass the test-cases of the present paper with their HLLC scheme and to share the results. The authors would also like to warmly thank Donald Schwendeman and Michael Hennessey who also agreed to pass the test-cases with the Schwendeman-Wahle-Kapila's Godunov-type scheme and to share the results. In both cases, their remarkable reactivity has contributed to improving this paper. Finally, the authors would like to thank the reviewers for their constructive remarks. This work has been partially funded by ANRT and EDF through an EDF-CIFRE contract 529/2009.

Appendices

7.1 Construction of the solution to the Riemann problem (25)- [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF].

Given (W L , W R , a 1 , a 2 ) (satisfying T k,L = τ k,L and T k,R = τ k,R for k ∈ {1, 2}) such that the conditions of Theorem 3.5 are met, we give the expression of the piecewise constant solution of the Riemann problem ( 25)- [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF]. For the sake of simplicity, the solution will be expressed in non conservative variables W = (α 1 , τ 1 , τ 2 , u 1 , u 2 , π 1 , π 2 , E 1 , E 2 ).

In practice, when implementing the numerical scheme, the relaxation Riemann solution of ( 25)-( 32) is used to compute the numerical fluxes at each interface between two states (U L , U R ) and the relaxation states (W L , W R ) are actually computed from these two states (U L , U R ). For this reason, the solution will be denoted

We recall the following notations built on the initial states (W L , W R ) (and therefore depending on (U L , U R )) and on the relaxation parameters (a 1 , a 2 ), which are useful for the computation of the solution.

For k in {1, 2}:

We also recall the dimensionless number of equation ( 48):

and define as U ♯ (U L , U R ; a 1 , a 2 ) the central expression of assumption (A) of Theorem 3.5:

Later in this section, we will omit the dependency of these quantities on (U L , U R ; a 1 , a 2 ). Following Theorem 3.5, if a 1 and a 2 are such that τ

L holds true, then there exists a self-similar solution to the Riemann problem ( 25)- [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF]. Following [START_REF] Coquel | A robust entropy-satisfying finite volume scheme for the isentropic Baer-Nunziato model[END_REF], we distinguish three different cases corresponding to different orderings of the kinematic waves,

With each one of these wave configurations is associated a different expression of assumption (A) depending on the sign of U ♯ .

Solution with the wave ordering u

-The intermediate states, which are represented in Figure 9, and the velocities u * 1 (U L , U R ; a 1 , a 2 ) and u * 2 (U L , U R ; a 1 , a 2 ) (simply denoted u * 1 and u * 2 hereafter) are computed through the following steps performed in the very same order.

• Relaxation pressures π 1 (τ 1 , T 1 , s 1 ):

). . 7. The intermediate states for phase 2 are then given by

2 ). .

Solution with the wave ordering u

For the determination of the wave velocities and the intermediate states, the simplest thing to do is to exploit the Galilean invariance of the equations. In this case indeed, the solution is obtained by the transformation

where the operator V changes the velocities into their opposite values:

V : (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 , x 8 , x 9 ) → (x 1 , x 2 , x 3 , -x 4 , -x 5 , x 6 , x 7 , x 8 , x 9 ).

(82)

Of course, the function ξ → W(ξ; VU R , VU L ; a 1 , a 2 ) is computed through the first case, since for these new initial data (VU R , VU L ), it is condition (A1) that holds.

Solution with the wave ordering u

The kinematic velocities are given by u

The intermediate states for phase 2 are obtained by the same formulae as in the case u * 2 < u * 1 , while the intermediate states for phase 1 (see Figure 10) read