Frédéric Coquel

Jean-Marc Hérard

Khaled Saleh

Nicolas Seguin

A Robust Entropy-Satisfying Finite Volume Scheme for the Isentropic Baer-Nunziato Model

Keywords: Two-phase flows, entropy-satisfying methods, relaxation techniques, Riemann problem. AMS subject classifications : 76T05, 35L60, 35F55

We construct an approximate Riemann solver for the isentropic Baer-Nunziato two-phase flow model, that is able to cope with arbitrarily small values of the statistical phase fractions. The solver relies on a relaxation approximation of the model for which the Riemann problem is exactly solved for subsonic relative speeds. In an original manner, the Riemann solutions to the linearly degenerate relaxation system are allowed to dissipate the total energy in the vanishing phase regimes, thereby enforcing the robustness and stability of the method in the limits of small phase fractions. The scheme is proved to satisfy a discrete entropy inequality and to preserve positive values of the statistical fractions and densities. The numerical simulations show a much higher precision and a more reduced computational cost (for comparable accuracy) than standard numerical schemes used in the nuclear industry. Finally, two test-cases 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] or the re-flooding phase following a LOCA. 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. The derived schemes are expected to ensure important stability properties such as the positivity of the densities and discrete entropy inequalities. In addition, as explicit schemes are needed for the simulation of these potentially highly unsteady phenomena, one major challenge is the control of the time step. In this context, the aim of this work is to design a robust and entropy-satisfying scheme for the numerical approximation of two-phase flows with vapor or liquid fractions arbitrarily close to zero.

The model concerned by this paper is the isentropic version of the two-fluid model introduced by Baer and Nunziato 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] 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] (see also [START_REF] Saurel | A multiphase godunov method for compressible multifluid and multiphase flows[END_REF] for a related framework). The model consists in 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 to assume 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 may be added in order to account for this pressure relaxation phenomenon as well as friction terms for the relaxation of the phasic velocities towards an equilibrium. However, this work is mainly concerned with the convective effects and these relaxation source terms are not considered here (see [START_REF] Coquel | Closure laws for a two-fluid two pressure model[END_REF] for some modeling choices of these terms and [START_REF] Hérard | A fractional step method to compute a class of compressible gas-luiquid flows[END_REF] for their numerical treatment).

Contrary to the single pressure models, the Baer-Nunziato model provides a pleasant property which is the weak hyperbolicity of its convective part. Indeed, unlike single pressure models, where the characteristic eigenvalues may be complex, the Baer-Nunziato model admits seven real eigenvalues and the associated right eigenvectors form a basis unless the relative velocity between the phases equals the speed of sound in the liquid (see [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF]). A first consequence is that the initial value problem is not a priori ill-posed for this model as long as the system remains strictly hyperbolic, that is to say as long as no interaction occurs between the acoustic waves and the material wave transporting the phase fraction. Another important consequence is that the definition of the non-conservative term in the model is not ambiguous, still as long as the system remains strictly hyperbolic. Indeed, the non-conservative term in the Baer-Nunziato model is transported by the linearly degenerate wave fraction wave. Hence, as long as the system is hyperbolic, this non-conservative product is naturally defined by the Riemann invariants associated with this wave (see [START_REF] Lefloch | Shock waves for nonlinear hyperbolic systems in nonconservative form[END_REF][START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF]). One of these Riemann invariants expresses the conservation of the total mixture energy across the phase fraction wave in the domain of hyperbolicity of the model. When the strict hyperbolicity is not satisfied, one speaks of a resonance phenomenon and the definition of the non-conservative product is no longer straightforward. In these cases, one must define regularization models in order to prescribe the behavior of such resonant solutions. One usual model used in a related framework [START_REF] Isaacson | Convergence of the 2 × 2 Godunov method for a general resonant nonlinear balance law[END_REF][START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF] that encompasses the Baer-Nunziato model, assumes a monotonic evolution of the phase fraction within the resonant wave. Other regularizations processes have been introduced in the context of non-conservative interface coupling between hyperbolic systems [START_REF] Ambroso | The coupling of homogeneous models for two-phase flows[END_REF][START_REF] Boutin | Coupling nonlinear hyperbolic equations (iii). A regularization method based on thick interfaces[END_REF]. In [START_REF] Sainsaulieu | Contribution à la modélisation mathématique et numérique des écoulements diphasiques constitués d'un nuage de particules dans un écoulement de gaz[END_REF], viscous regularization methods are used for related twophase flow models and in the works of Yong [START_REF] Yong | Entropy and global existence for hyperbolic balance laws[END_REF], Kawashima-Yong [START_REF] Kawashima | Dissipative structure and entropy for hyperbolic systems of balance laws[END_REF], Chen et al. [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF] and Hanouzet et al. [START_REF] Hanouzet | Global existence of smooth solutions for partially dissipative hyperbolic systems with a convex entropy[END_REF], the authors investigate the stabilization effects of zero-th order source terms in the context of strictly hyperbolic balance equations. One common consequence of all these stabilization procedures is that they imply a dissipation of the system's energy either globally (for viscous and zero-th order stabilizations) or locally around the resonant wave. In this latter case, the regularization is associated with a so-called kinetic relation which is a generalized Rankine-Hugoniot type relation on the system's energy expressing its dissipation through the resonant linearly degenerate wave. For a general review on kinetic relations, we refer to [START_REF] Berthon | Why many theories of shock waves are necessary: kinetic relations for non-conservative systems[END_REF].

In our context of interest, namely nuclear liquid-vapor flows, the resonance due to wave interaction between acoustic fields and the phase transport equation is unlikely to arise since it would imply sonic or supersonic values of the relative velocity between the phases. However, the aim of this work is to capitalize on the above considerations concerning the stabilization of resonant solutions in order to propose an original framework for the stable computation of vanishing phase solutions. The method is the following. Among all the energy dissipative mechanisms, we seek those that might ensure stable approximate solutions for arbitrarily small values of the phase fractions. For instance, an interesting question would be to know if in the Riemann problem, weakening the energy conservation across the phase fraction wave through a suitable kinetic relation, may ensure stable solutions in the regimes of small phase fractions. Unfortunately, despite some interesting works in that sense ( [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF][START_REF] Deledicque | A conservative approximation to compressible two-phase flow models in the stiff mechanical relaxation limit[END_REF][START_REF] Andrianov | The Riemann problem for the Baer-Nunziato two-phase flow model[END_REF]), calculating the exact solution of the Riemann problem for the Baer-Nunziato model with any initial data is so far out of reach. One main obstruction is that the characteristic eigenvalues of the system are not naturally ordered, and no method has been found yet that could determine a priori their ordering, with respect to the initial data. In addition, the strong non-linearities of the pressure laws make even more difficult the derivation of an exact Riemann solver. Following the pioneering work of Harten, Lax and van Leer [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF], other approaches consider approximate Riemann solvers [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], and the method described in the present paper enters this category. Let us mention some other schemes grounded on 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] Liu | [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] Thanh | Numerical approximation for a Baer-Nunziato model of twophase flows[END_REF].

Actually, the method considered here relies on a relaxation approximation of the model, similar to that in Ambroso, Chalons, Coquel and Galié [START_REF] Ambroso | Relaxation and numerical approximation of a two-fluid two-pressure diphasic model[END_REF]. The idea consists in introducing a larger system, in which the pressure laws have been linearized, and which relaxes towards the actual system of Baer-Nunziato in the regime of a small relaxation parameter (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]). While in [START_REF] Ambroso | Relaxation and numerical approximation of a two-fluid two-pressure diphasic model[END_REF], the authors calculate an approximate solution of the relaxation system thanks to a prediction of the non-conservative product, one major contribution of our work is the effective resolution of the relaxation Riemann problem in the framework of subsonic wave ordering, by calculating exact solutions. In particular, we provide explicit conditions on the initial data that enable the a priori determination of the relative ordering of the waves. The proof relies on a fixed-point procedure that consists in iteratively considering the evolution equations of each phase, taken separately. An outstanding property is that this iterative procedure boils down to a fixed-point research on a monotonic scalar function. This property is crucial for the numerical computations using the relaxation Riemann solver since it avoids heavy computational costs as assessed by the numerical tests. In order to enforce the stability of the solutions when the initial phase fractions have arbitrarily small values, the solutions are allowed in these regimes to dissipate the total energy of the system across the phase fraction wave, through the definition of a suitable kinetic relation. Hence, despite the linear degeneracy of all the waves implied by the relaxation approximation procedure, some extra dissipation is introduced and one may speak of dissipative relaxation. This original idea has been introduced in order to control the time step in the vanishing phase regimes. Indeed, another way of stabilizing the numerical solutions in these regimes is to take larger values of the relaxation parameters, thus introducing more diffusion to the relaxation approximation. However, too large values of the relaxation parameters may involve too small time steps because of the CFL restrictions, since these parameters directly control the size of the wave fan in the approximate Riemann solver.

The resulting scheme is proved to preserve positive densties and to satisfy a discrete entropyinequality under a sub-characteristic condition (Whitham's condition). To our knowledge, there exists no other scheme that is proved to satisfy these two properties. In addition, for the same level of refinement, the scheme is shown to be much more accurate than the Rusanov 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. Actually, comparing with Lax-Friedrichs type schemes is quite significant since for such stiff configurations as vanishing phase cases, these schemes are commonly used in the industrial context because of their known robustness [START_REF] Hérard | A fractional step method to compute a class of compressible gas-luiquid 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] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF] for example of high order schemes). Nevertheless this work is focused on the design of stable computations for vanishing phase configurations. In these regimes indeed, even high order methods may develop instabilities.

The paper is organized as follows. Sections 2 and 3 are devoted to the presentation of the Baer-Nunziato model and its relaxation approximation. Some issues concerning the resonance, kinetic relations and vanishing phase regimes are also discussed. Section 4 is the core of the paper. It displays an existence theorem for the relaxation Riemann problem, for a generalized class of dissipative solutions. For the sake of understanding, the constructive proof is explained in detail. In section 5, we use the relaxation exact Riemann solver to derive an approximate solver in the sense of Harten, Lax and van Leer [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] and thus, a Finite Volume scheme. Finally, section 6 is devoted to the numerical tests. In addition to a convergence and CPU cost study, two test-cases assess that the scheme provides a robust numerical treatment of vanishing phase solutions.

The isentropic model of Baer-Nunziato

The isentropic Baer-Nunziato model is a two-phase flow model formulated in Eulerian coordinates where balance equations account for the evolution of mass and momentum of each phase. The velocities of each phase are denoted u i , i ∈ {1, 2}, while the densities are denoted ρ i , i ∈ {1, 2}. Each phase has a statistical phase fraction α i , i ∈ {1, 2}, with the saturation constraint α 1 + α 2 = 1. The model reads:

∂ t U + ∂ x f(U) + c(U)∂ x U = 0, x ∈ R, t > 0, (2.1)
with

U =       α 1 α 1 ρ 1 α 1 ρ 1 u 1 α 2 ρ 2 α 2 ρ 2 u 2       , f(U) =       0 α 1 ρ 1 u 1 α 1 ρ 1 u 2 1 + α 1 p 1 (ρ 1 ) α 2 ρ 2 u 2 α 2 ρ 2 u 2 2 + α 2 p 2 (ρ 2 )       , c(U)∂ x U =       u 2 0 -p 1 (ρ 1 ) 0 +p 1 (ρ 1 )       ∂ x α 1 .
(2.2)

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

Ω = U ∈ R 5 , 0 < α 1 < 1 and α i ρ i > 0 for i ∈ {1, 2} . (2.3)
We assume barotropic pressure laws for each phase ρ i → p i (ρ i ), i ∈ {1, 2} with smooth dependence on the density, and which satisfy the following natural assumptions for all ρ i > 0:

p i (ρ i ) > 0, p i (ρ i ) > 0, lim ρi→0 p i (ρ i ) = 0, lim ρi→+∞ p i (ρ i ) = +∞. (2.4)
We define the mapping τ → P i (τ ) := p i (τ -1 ) which is the phasic pressure seen as a function of the specific volume τ = ρ -1 . In the whole paper, this smooth function is assumed to be strictly convex:

P i (τ i ) > 0, for all τ i > 0, i ∈ {1, 2}. (2.5)

Main mathematical properties

The following proposition characterizes the fields of this system: Proposition 2.1. System (2.1) is weakly hyperbolic since it admits the following real eigenvalues

σ 1 (U) = u 2 , σ 2 (U) = u 1 -c 1 (ρ 1 ), σ 3 (U) = u 1 + c 1 (ρ 1 ), σ 4 (U) = u 2 -c 2 (ρ 2 ), σ 5 (U) = u 2 + c 2 (ρ 2 ), (2.6
) where c i (ρ i ) = p i (ρ i ) is the speed of sound for phase i. The corresponding right eigenvectors are linearly independent if, and only if,

α 1 = 0, α 2 = 0, |u 1 -u 2 | = c 1 (ρ 1 ).
(2.7)

When (2.7) is not satisfied, the system is said to be resonant. The characteristic fields associated with σ 2 , σ 3 , σ 4 and σ 5 are genuinely non-linear, while the characteristic field associated with σ 1 is linearly degenerate.

Remark 2.1. Actually, following the definition of the admissible physical space (2.3), one never has α 1 = 0 or α 2 = 0. However, α i = 0 is to be understood in the sense α i → 0 since the aim of this work is to construct a robust enough Riemann solver that could handle all the possible values of α i , i ∈ {1, 2}, especially, arbitrarily small values.

Proof . Denoting U = (α 1 , ρ 1 , u 1 , ρ 2 , u 2 ) T , the smooth solutions of system (2.1) satisfy the equivalent following system

∂ t U + A(U)∂ x U = 0,
where

A(U) =        u 2 0 0 0 0 ρ1 α1 (u 1 -u 2 ) u 1 ρ 1 0 0 0 p 1 (ρ1) ρ1 u 1 0 0 0 0 0 u 2 ρ 2 p1(ρ1)-p2(ρ2) (1-α1)ρ2 0 0 p 2 (ρ2) ρ2 u 2        .
This matrix admits five real eigenvalues that are provided in (2.6). Denoting M = u1-u2 c1 the relative Mach number, the corresponding right eigenvectors can be chosen as

r 1 (U) = α 1 (1 -α 1 )(M 2 -1), -(1 -α 1 )ρ 1 M 2 , (1 -α 1 )c 1 M, α 1 p 2 (ρ 2 ) -p 1 (ρ 1 ) c 2 2 (M 2 -1), 0 T , r 2 (U) = -2 p 1 (ρ 1 ) p 1 (ρ 1 )ρ 1 + 2p 1 (ρ 1 ) 0, ρ 1 , -p 1 (ρ 1 ), 0, 0 T , r 3 (U) = 2 p 1 (ρ 1 ) p 1 (ρ 1 )ρ 1 + 2p 1 (ρ 1 ) 0, ρ 1 , p 1 (ρ 1 ), 0, 0 T , r 4 (U) = -2 p 2 (ρ 2 ) p 2 (ρ 2 )ρ 2 + 2p 2 (ρ 2 ) 0, 0, 0, ρ 2 , -p 2 (ρ 2 ) T , r 5 (U) = 2 p 2 (ρ 2 ) p 2 (ρ 2 )ρ 2 + 2p 2 (ρ 2 ) 0, 0, 0, ρ 2 , p 2 (ρ 2 ) T .
Hence, the system is hyperbolic, i.e. the five eigenvectors r k (U), k in {1, .., 5} span R 5 if and only if

α 1 = 0, α 2 = 0 and (u 1 -u 2 ) 2 = c 2 1 (ρ 1
). Moreover, one can easily verify that

∇ U σ 1 (U).r 1 (U) = 0, ∇ U σ k (U).r k (U) = 1, k = 2, 3, 4, 5
which proves that the first characteristic field is linearly degenerate while the four others are genuinely non-linear.

The following proposition states the existence of two phasic energy equations satisfied by the smooth solutions of system (2.1).

Proposition 2.2. Defining E

i := E i (u i , τ i ) = u 2 i
2 + e i (τ i ), i ∈ {1, 2} with τ → e i (τ ) an antiderivative of τ → -P i (τ ), the smooth solutions of system (2.1) satisfy the following phasic energy equations:

∂ t (α i ρ i E i ) + ∂ x (α i ρ i E i + α i p i (ρ i )) u i -u 2 p 1 (ρ 1 )∂ x α i = 0, i ∈ {1, 2}.
(2.8)

Summing over i = 1, 2 yields the following additional conservation law, which expresses the total mixture energy conservation by the smooth solutions of system (2.1):

∂ t (α 1 ρ 1 E 1 + α 2 ρ 2 E 2 ) + ∂ x ((α 1 ρ 1 E 1 + α 1 p 1 (ρ 1 )) u 1 + (α 2 ρ 2 E 2 + α 2 p 2 (ρ 2 )) u 2 ) = 0. (2.9)
Proof . The proof is classical. It consists in multiplying the phasic momentum equations (third and fifth components of system (2.1)) respectively by u i and re-arranging the equations using the definition of e i (τ i ) and the mass conservation equations (second and fourth components of (2.1)).

As regards the non-smooth weak solutions of (2.1), one has to add a so-called entropy criterion in order to select the relevant physical solutions. Thus, an entropy weak solution of (2.1) is a function U(x, t) that satisfies (2.1) in the sense of distributions as well as the following entropy inequality:

∂ t (α 1 ρ 1 E 1 + α 2 ρ 2 E 2 ) + ∂ x ((α 1 ρ 1 E 1 + α 1 p 1 (ρ 1 )) u 1 + (α 2 ρ 2 E 2 + α 2 p 2 (ρ 2 )) u 2 ) ≤ 0.
(2.10)

When the solution contains shock waves, inequality (2.10) is strict in order to account for the physical loss of energy due to viscous phenomena that are not modeled in system (2.1).

Riemann weak solutions in the case of strict hyperbolicity

The definition of weak solutions to system (2.1) is not straightforward because of the non conservative product p 1 ∂ x α 1 . Indeed, considering a discontinuity of the phase fraction α 1 across which the pressure p 1 may be discontinuous, it is clear that the product p 1 ∂ x α 1 is not a priori defined in the sense of distributions. However, since α 1 discontinuities coincide with a linearly degenerate field, there is no ambiguity in the definition of the product as long as the system is hyperbolic, i.e. as long as the three conditions (2.7) are met. In this case indeed, there are two equivalent ways for defining the product p 1 ∂ x α 1 . Let us describe briefly these too equivalent alternatives in the context of a Riemann problem. System (2.1) is supplied with an initial condition of the type

U(x, t = 0) = U L if x < 0, U R if x > 0. (2.11)
In most cases (hyperbolic cases, see thereafter), if the initial data is such that α 1,L = α 1,R , the solution of the Riemann problem is a self-similar function composed of intermediate states separated by waves associated with the acoustic fields (rarefaction or shock waves) or by a contact discontinuity associated with the eigenvalue u 2 . Actually, the phase fractions α i only jump through this contact discontinuity from α i,L to α i,R . Away from this wave, the phase fractions are constant and the system behaves as two independent (isentropic) Euler systems associated with both phases since the term ∂ x α i vanishes.

u 1 + c 1 u 2 + c 2 u * 2 u 1 -c 1 u 2 -c 2 u 1 -c 1 α 1,R α 1,L U - U +
For the actual computation of the Riemann solution, one needs jump relations across each wave in order to link the different intermediate states. For the acoustic waves, the jump relations are the same as for the Euler isentropic equations associated with each phase. For the u 2 -wave however, the two phases are coupled and there are two equivalent ways of defining the relation between two intermediate states U -and U + separated by such a discontinuity.

Definition through Riemann invariants

When the system is hyperbolic, there are four Riemann invariants associated with the linearly degenerate field σ 1 (U) = u 2 , that is to say four independent functions Φ k (U), k = 1..4 satisfying ∇ U Φ k (U).r 1 (U) = 0, for all k = 1..4. In the hyperbolic case, two states U -and U + are separated by a u 2 -contact discontinuity if and only if Φ k (U -) = Φ k (U + ) for all k = 1..4. For a proof of this result in the context of the full Baer-Nunziato model with energies, we refer to [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF]. The derivative ∂ x α 1 corresponds to a Dirac measure of the form

(α + 1 -α - 1 )δ x-u * 2 t
, where u * 2 is the propagation speed of the α 1 -discontinuity, and the non-conservative product p 1 ∂ x α 1 identifies with a measure

p 1 ∂ x α 1 := p * 1 (α + 1 -α - 1 )δ x-u * 2 t
(2.12)

whose weight p * 1 (α + 1 -α - 1
) is obtained by applying Rankine-Hugoniot's jump relation equivalently to one of the momentum equations:

p * 1 (α + 1 -α - 1 ) = -u * 2 [α 1 ρ 1 u 1 ] x t =u * 2 + [α 1 ρ 1 u 2 1 + α 1 p 1 ] x t =u * 2 = --u * 2 [α 2 ρ 2 u 2 ] x t =u * 2 + [α 2 ρ 2 u 2 2 + α 2 p 2 ] x t =u * 2 .
(2.13)

Here, [X] x t =u * 2 = X + -X -denotes the difference between the values taken by the quantity X on the right and on the left of the u 2 -contact discontinuity. As the eigenvalue u 2 is a Riemann invariant, the second line of (2.13) gives

p 1 ∂ x α 1 = p * 1 (α + 1 -α - 1 )δ x-u * 2 t = -[α 2 p 2 ] x t =u * 2 δ x-u * 2 t .
(2.14)

Definition through Rankine-Hugoniot jump relations

In the hyperbolic case, an equivalent way of defining the non-conservative product p 1 ∂ x α 1 across a u 2 -contact discontinuity is through Rankine-Hugoniot jump relations. A first relation is given by the continuity of the eigenvalue u 2 across this linearly degenerate wave, and we get two more independent jump relations by applying Rankine-Hugoniot's formula to the mass conservation equation of phase 1 and to the total momentum conservation:

[u 2 ] x t =u * 2 = 0, (2.15)
-u * 2 [α 1 ρ 1 ] x t =u * 2 + [α 1 ρ 1 u 1 ] x t =u * 2 = 0, (2.16
)

-u * 2 [α 1 ρ 1 u 1 + α 2 ρ 2 u 2 ] x t =u * 2 + α 1 ρ 1 u 2 1 + α 1 p 1 + α 2 ρ 2 u 2 2 + α 2 p 2 x t =u * 2 = 0.
(2.17)

Observe that the mass conservation equation for phase 2 gives no additional information since the corresponding jump relation is naturally satisfied because u 2 is constant through the discontinuity.

Hence, there is one missing information in order to fully define the discontinuity and therefore the non-conservative product. If we assume the hyperbolicity of the system, this last jump relation is obtained by applying Rankine-Hugoniot's formula to the total energy conservation (2.9) which yields

-u * 2 [α 1 ρ 1 E 1 + α 2 ρ 2 E 2 ] x t =u * 2 + [(α 1 ρ 1 E 1 + α 1 p 1 ) u 1 + (α 2 ρ 2 E 2 + α 2 p 2 ) u 2 ] x t =u * 2 = 0. (2.18)
Indeed, a theoretical result (see [START_REF] Lefloch | Shock waves for nonlinear hyperbolic systems in nonconservative form[END_REF]) states that for a hyperbolic system, any additional conservation law satisfied by the smooth solutions of the system is also satisfied in the weak sense through contact discontinuities. Actually, a sufficient condition for the proof of this result is the existence of a right eigenvector basis. The non-conservative product is then obtained by applying Rankine-Hugoniot's relation to the momentum equations just as previously, which yields (2.14).

Local energy equations around the α 1 -discontinuity

In the case of hyperbolicity, the phasic energy equations (2.8) are satisfied in the weak sense in the neighborhood of a u 2 -contact discontinuity. For phase 2, this is a direct consequence of the fact that u 2 is constant across this wave according to (2.15). Indeed, one may write

-u * 2 [α 2 ρ 2 E 2 ] x t =u * 2 + [(α 2 ρ 2 E 2 + α 2 p 2 ) u 2 ] x t =u * 2 -u * 2 [α 2 p 2 ] = 0, (2.19)
which by the definition (2.14) of the non-conservative product

p 1 ∂ x α 1 , is the Rankine-Hugoniot relation of ∂ t (α 2 ρ 2 E 2 ) + ∂ x (α 2 ρ 2 E 2 + α 2 p 2 (ρ 2 )) u 2 -u 2 p 1 (ρ 1 )∂ x α 2 = 0. (2.20)
As for phase 1, the local phasic energy equation is a consequence of the constancy (2.15) of u 2 and of the conservation of the total mixture energy (2.18). Thus in the hyperbolicity areas of the system, the weak solutions satisfy

-u * 2 [α 1 ρ 1 E 1 ] x t =u * 2 + [(α 1 ρ 1 E 1 + α 1 p 1 ) u 1 ] x t =u * 2 -u * 2 [α 2 p 2 ] = 0, (2.21)
which is the weak formulation of the following PDE satisfied locally near the α 1 -discontinuity:

∂ t (α 1 ρ 1 E 1 ) + ∂ x (α 1 ρ 1 E 1 + α 1 p 1 (ρ 1 )) u 1 -u 2 p 1 (ρ 1 )∂ x α 1 = 0. (2.22)
Let us emphasize again, that this last property (2.21)-(2.22) is a consequence of the assumed hyperbolicity of the system since, this is the very property that allows the total energy conservation to be used in the weak sense across a u 2 -contact discontinuity. If the system is not strictly hyperbolic, this property is not guaranteed as explained in the next section.

Riemann weak solutions in the resonant case

When (2.7) is not satisfied, for instance if at some point of the flow, the physical quantities satisfy |u 1u 2 | = c 1 (ρ 1 ), the hyperbolicity is lost and the system is said to be resonant. In that case, the existence of four independent Riemann invariants is no longer guaranteed and the jump relation associated with the total mixture energy conservation is potentially lost. Indeed, the existing theoretical results stating that for a hyperbolic system, any additional conservation law satisfied by the smooth solutions is also satisfied in the weak sense through contact discontinuities, assume the existence of a right eigenvector basis. Therefore, in these resonant cases, if the total mixture energy can no more be conserved across the u 2 -contact discontinuity, enforcing its dissipation appears to be necessary for stability reasons. To this end, one has to model regularization processes along the u 2 -contact discontinuity by adding some information on the expected behavior of such resonant solutions. With any regularization process, is associated a so-called kinetic relation, which is an additional Rankine-Hugoniot type relation, that enables the actual computation of jumps across a resonant contact discontinuity, for instance when solving the Riemann problem. In our context, a kinetic relation takes the form

-u * 2 [α 1 ρ 1 E 1 + α 2 ρ 2 E 2 ] x t =u * 2 + [(α 1 ρ 1 E 1 + α 1 p 1 ) u 1 + (α 2 ρ 2 E 2 + α 2 p 2 ) u 2 ] x t =u * 2 = -Q(u * 2 , U L , U R ), (2.23) where Q(u * 2 , U L , U R
) is a positive number measuring the amount of dissipated energy across the u 2wave. In practice, the definition of such a kinetic relation has to obey several natural requirements. First of all, in the areas of the phase space Ω where the system is hyperbolic, Q(u * 2 , U L , U R ) is expected to vanish in order to restore the total mixture energy conservation. On the other hand, in the resonance areas, Q(u * 2 , U L , U R ) has to be designed so as to ensure the existence of stable solutions of the Riemann problem. Of course, several choices could be made for the regularization process, leading to several possible definitions of Q(u * 2 , U L , U R ). Thus, the uniqueness of solutions can only be recovered for a given choice of the kinetic relation.

To exemplify the energy dissipation across a u 2 -contact discontinuity in the resonant case, one may consider the case of a {u 1 -c 1 }-shock superimposing with the u 2 -contact discontinuity. If one considers the jump relation of the total energy through such a discontinuity, it reads

-u * 2 [α 1 ρ 1 E 1 + α 2 ρ 2 E 2 ] x t =u * 2 + [(α 1 ρ 1 E 1 + α 1 p 1 ) u 1 + (α 2 ρ 2 E 2 + α 2 p 2 ) u 2 ] x t =u * 2 < 0. (2.24)
Actually, on the one hand, the energy equation for phase 2 (2. [START_REF] Deledicque | A conservative approximation to compressible two-phase flow models in the stiff mechanical relaxation limit[END_REF])-(2.20) is still satisfied. On the other hand, the energy of phase 1 is strictly dissipated in this case due to the shock lying inside the discontinuity. The jump relation (2.24) is equivalent to

-u * 2 [α 1 ρ 1 E 1 ] x t =u * 2 + [(α 1 ρ 1 E 1 + α 1 p 1 ) u 1 ] x t =u * 2 -u * 2 [α 2 p 2 ] < 0, (2.25)
which is the weak formulation of

∂ t (α 1 ρ 1 E 1 ) + ∂ x (α 1 ρ 1 E 1 + α 1 p 1 (ρ 1 )) u 1 -u 2 p 1 (ρ 1 )∂ x α 1 < 0. (2.26)
Considering now the physically relevant case of nearly vanishing phases, where the phase fractions may be arbitrarily close (but not equal) to zero, the strict hyperbolicity of the system is not actually lost in the sense that the right eigenvector basis still exists as long as α i = 0, i ∈ {1, 2}. Nevertheless, the matrix composed of the right eigenvectors becomes singular as one of the phase fractions tends to zero. Obviously, this may cause instabilities of the model in these vanishing phase regimes. Yet, in many application contexts where one of the phases is naturally expected to disappear, one is interested in stable (approximate) computations of such stiff regimes. The purpose of this paper is to construct an approximate Riemann solver, for which the solutions are controlled in the regimes α i → 0. This Riemann solver relies on a relaxation approximation which is presented in next section.

Remark 2.2. In the physical configurations aimed at in the present work (such as two-phase flows in nuclear reactors), resonant configurations corresponding to wave interaction between acoustic fields and the u 2 -contact discontinuity are unlikely to occur. Indeed, for these applications, 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. (2.

27)

Having set aside the sonic (M = 1) and supersonic (M > 1) flows, the main purpose of this work is to tackle the issue of vanishing phase instabilities.

A relaxation approximation

We now introduce a relaxation approximation of system (2.1). The new system possesses the same structure as the initial model (2.1) but the differential part is linearly degenerate while the non-linear part is postponed to the source term. The goal is to construct and analyze the linearly degenerate differential part. In general, the Riemann problem for such systems is simple to solve but here, we inherit the non-conservative and resonant structure from (2.1). As a consequence, we have to introduce a decrease of the associated mixture energy in some extreme cases by the use of a so-called kinetic relation in order to obtain admissible solutions.

Denoting W = (α 1 , α 1 ρ 1 , α 1 ρ 1 u 1 , α 2 ρ 2 , α 2 ρ 2 u 2 , α 1 ρ 1 T 1 , α 2 ρ 2 T 2 )
T the relaxation state vector, we propose the following relaxation approximation for system (2.1):

∂ t W ε + ∂ x g(W ε ) + d(W ε )∂ x W ε = 1 ε R(W ε ), (3.1)
where

g(W) =           0 α 1 ρ 1 u 1 α 1 ρ 1 u 2 1 + α 1 π 1 (τ 1 , T 1 ) α 2 ρ 2 u 2 α 2 ρ 2 u 2 2 + α 2 π 2 (τ 2 , T 2 ) α 1 ρ 1 T 1 u 1 α 2 ρ 2 T 2 u 2           , d(W)∂ x W =           u 2 0 -π 1 (τ 1 , T 1 ) 0 +π 1 (τ 1 , T 1 ) 0 0           ∂ x α 1 , R(W) =           0 0 0 0 0 α 1 ρ 1 (τ 1 -T 1 ) α 2 ρ 2 (τ 2 -T 2 )           . (3.2)
For each phase i in {1, 2} the linearized pressure π i (τ i , T i ) is a function defined as

π i (τ i , T i ) = P i (T i ) + a 2 i (T i -τ i ). (3.3)
We can see that in the formal limit ε → 0, the additional variable T i tends towards the specific volume τ i , and the linearized pressure π i tends towards the original non-linear pressure p i , thus recovering system (2.1) in the first five equations of (3.1). From this point and to ease the notation, we will omit the superscript ε . In the sequel, the original system (2.1) will be referred to as the equilibrium system as opposed to the relaxation system. The constants a i in (3.3) are two constant positive parameters that must be taken large enough sa as to satisfy the so-called Whitham condition:

a 2 i > max Ti∈Ki {-P i (T i )} (3.4)
where K i is a connected compact set of R that contains all the values of τ i and T i of the solution.

Main mathematical properties

Let us now focus on the convective part of system (3.1):

∂ t W + ∂ x g(W) + d(W)∂ x W = 0. (3.5)
The solutions are sought in the domain of positive densities ρ i and positive T i :

Ω r = W ∈ R 7 , 0 < α 1 < 1, α i ρ i > 0, α i ρ i T i > 0, i ∈ {1, 2} . (3.6)
Proposition 3.1. System (3.5) is weakly hyperbolic since it admits the following real eigenvalues

σ 1 (W) = u 1 -a 1 τ 1 , σ 2 (W) = u 1 , σ 3 (W) = u 1 + a 1 τ 1 , σ 4 (W) = u 2 -a 2 τ 2 , σ 5 (W) = σ 6 (W) = u 2 , σ 7 (W) = u 2 + a 2 τ 2 . (3.7)
The corresponding right eigenvectors are linearly independent if, and only if

α 1 = 0, α 2 = 0, |u 1 -u 2 | = a 1 τ 1 . (3.8)
All the characteristic fields associated with these eigenvalues are linearly degenerate.

Proof . The proof is similar to that of Proposition 2.1. It is left to the reader.

Remark 3.1. Here again, one never has α 1 = 0 or α 2 = 0. However, α i = 0 is to be understood in the sense α i → 0.

Unlike system (2.1), one remarkable property of the relaxation system (3.5) is the linear degeneracy of all the characteristic fields. This has the helpful consequence that jump relations can be easily derived through each wave. The relaxation approximation is therefore a pleasant way to get around the difficulties due to non-linearity (discrimination between shocks and rarefaction waves, jump relations for shocks...) which arise when solving the Riemann problem for (2.1).

In a similar way to that for the equilibrium system, we have balance equations on the phasic energies as well as a total mixture energy conservation equation satisfied by the smooth solutions of system (3.5): Proposition 3.2. The smooth solutions of system (3.5) satisfy the following phasic energy equations:

∂ t (α i ρ i E i ) + ∂ x (α i ρ i E i + α i π i ) u i -u 2 π 1 (τ 1 , T 1 )∂ x α i = 0, i ∈ {1, 2}. (3.9)
where the phasic energies are defined by

E i := E i (u i , τ i , T i ) = u 2 i 2 + e i (T i ) + π 2 i (τ i , T i ) -P 2 i (T i ) 2a 2 i , i ∈ {1, 2}. (3.10)
Summing over i = 1, 2 yields the following additional conservation law, also satisfied by the smooth solutions of system (3.5):

∂ t (α 1 ρ 1 E 1 + α 2 ρ 2 E 2 ) + ∂ x ((α 1 ρ 1 E 1 + α 1 π 1 ) u 1 + (α 2 ρ 2 E 2 + α 2 π 2 ) u 2 ) = 0. (3.11)
Proof . The proof is left to the reader.

Jump relations and kinetic relation for vanishing phases

Being given a pair of initial states (W L , W R ), we supply the homogeneous system (3.5) with the following initial condition

W(x, t = 0) = W L if x < 0, W R if x > 0. (3.12)
As for the equilibrium system, the solution of the Riemann problem is a self-similar function composed of intermediate states separated by contact discontinuities associated with the eigenvalues of the system. The phase fractions α i only jump through the u 2 -contact discontinuity from α i,L to α i,R . Away from this wave, the phase fractions are constant and the system behaves as two independent relaxation systems for the isentropic Euler equations of each phase

u 1 + a 1 τ 1 u 2 + a 2 τ 2 u *
For all the discontinuities except the one associated with the eigenvalue u 2 , the system is locally conservative (the product π 1 ∂ x α i locally vanishes) and the jump conditions are simply obtained by the Rankine-Hugoniot relations applied to the conservative equations of the system.

On the contrary, for the u 2 -wave, we have ∂ x α i = 0. In fact ∂ x α i identifies with a Dirac measure and the pressure π 1 may be discontinuous across this wave. The product π 1 ∂ x α i is determined by the Riemann invariants of this linearly degenerate field, or equivalently by the Rankine-Hugoniot jump relations. A first relation is given by the continuity of the eigenvalue u 2 , and we get three more independent jump relations by applying Rankine-Hugoniot's formula to the conservative equations of phase 1 and to the total momentum conservation:

[u 2 ] ξ=u * 2 = 0, (3.13) -u * 2 [α 1 ρ 1 ] ξ=u * 2 + [α 1 ρ 1 u 1 ] ξ=u * 2 = 0, (3.14)
-u * 2 [α 1 ρ 1 T 1 ] ξ=u * 2 + [α 1 ρ 1 u 1 T 1 ] ξ=u * 2 = 0, (3.15)
-u * 2 [α 1 ρ 1 u 1 + α 2 ρ 2 u 2 ] ξ=u * 2 + α 1 ρ 1 u 2 1 + α 1 π 1 + α 2 ρ 2 u 2 2 + α 2 π 2 ξ=u * 2 = 0. (3.16)
Finally, if the system is hyperbolic, a last jump relation (recall that the eigenvalue u 2 has multiplicity 2) is obtained by applying Rankine-Hugoniot's formula to the total energy preservation (3.11) which yields

-u * 2 [α 1 ρ 1 E 1 + α 2 ρ 2 E 2 ] ξ=u * 2 + [(α 1 ρ 1 E 1 + α 1 π 1 ) u 1 + (α 2 ρ 2 E 2 + α 2 π 2 ) u 2 ] ξ=u * 2 = 0. (3.17)
In the sequel, when solving the Riemann problem for given initial data (W L , W R ) and a given pair of parameters (a 1 , a 2 ), it appears that the exact solution may have non-positive densities if the ratio of initial phase fractions

α 1,L
α 1,R (or its inverse) is too large. One way to avoid this problem is to take larger values of the parameters (a 1 , a 2 ) and solve the problem again. At the numerical level, taking larger values of these parameters introduces some more diffusion to the relaxation approximation, thus stabilizing the simulation. However, too large values of these parameters may involve too small time steps because of the CFL restrictions, since (a 1 , a 2 ) directly control the size of the wave fan in the approximate Riemann solver. In this work, in order to control the time step, we propose to add the needed extra diffusion directly at the PDE level in the relaxation approximation, in the form of an energy dissipation. This is actually an original idea in the context of relaxation approximations, since usually, the first step of the method (treatment of the convective system (3.5)) is energy preserving because of the linear degeneracy of all the fields. Instead, the total energy conservation is replaced in the regimes of vanishing phases by a kinetic relation of the form

-u * 2 [α 1 ρ 1 E 1 + α 2 ρ 2 E 2 ] x t =u * 2 + [(α 1 ρ 1 E 1 + α 1 π 1 ) u 1 + (α 2 ρ 2 E 2 + α 2 π 2 ) u 2 ] x t =u * 2 = -Q(u * 2 , W L , W R ), (3.18) where Q(u * 2 , W L , W R
) is a positive scalar function measuring the amount of dissipated energy across the u 2 -wave. This kinetic relation allows to completely define the jump between two states W -and W + separated by a u 2 -contact discontinuity. The non-conservative product π 1 ∂ x α 1 then identifies with a Dirac measure

π * 1 ∆α 1 δ x-u * 2 t , ∆α 1 := α 1,R -α 1,L , (3.19)
whose weight is obtained by applying Rankine-Hugoniot jump relation to the momentum equations:

π * 1 ∆α 1 = -u * 2 [α 1 ρ 1 u 1 ] ξ=u * 2 + α 1 ρ 1 u 2 1 + α 1 π 1 ξ=u * 2 = --u * 2 [α 2 ρ 2 u 2 ] ξ=u * 2 + α 2 ρ 2 u 2 2 + α 2 π 2 ξ=u * 2 = -[α 2 π 2 ] ξ=u * 2 .
(3.20)

Solving the Riemann problem for the relaxation system

The aim of this section is to solve the Riemann problem associated with the homogeneous part of the relaxation system. Being given a pair of initial states (W L , W R ), we seek solutions of the following Cauchy problem:

∂ t W + ∂ x g(W) + d(W)∂ x W = 0, (4.1)
with the initial condition

W(x, t = 0) = W L if x < 0, W R if x > 0. (4.2)

Definition of the solutions to the Riemann problem

The solution is sought in the form of a self-similar function only depending on the variable ξ = x t , that is to say W(x, t) = W r (x/t). As all the fields are linearly degenerate, if the solution remains in the domain of hyperbolicity, the function W r (ξ) is a piecewise constant function, composed of (at most) seven constant states separated by (at most) six contact discontinuities, associated with the six eigenvalues u 1 ± a 1 τ 1 , u 2 ± a 2 τ 2 , u 1 and u 2 (see Lax's theory for Riemann problems [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]). More precisely, we give the following definition for the solutions of the Riemann problem (4.1)-(4.2). 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 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 i∈{1,2} {u i,L -a i τ i,L } =⇒ W r (ξ; W L , W R ) = W L , ξ > max i∈{1,2} {u i,R + a i τ i,R } =⇒ W r (ξ; W L , W R ) = W R . (4.3)
(ii) There exists a real number u * 2 (depending on (W L , W R )) with the following properties. Denoting

Ξ -= {ξ ∈ R such that ξ < u * 2 } and Ξ + = {ξ ∈ R such that ξ > u * 2 }, α 1 (ξ) equals α 1,L on Ξ - and α 1,R on Ξ + and the function W r (ξ; W L , W R ) is a weak solution of -ξW r + g(W r ) = 0, ξ ∈ Ξ -∪ Ξ + . ( 4

.4)

(iii) There exists a real number π * 1 (depending on (W L , W R )) such that the following jump relations across the u * 2 -discontinuity hold

[u 2 ] ξ=u * 2 = 0, (4.5)
-u * 2 [W r ] ξ=u * 2 + [g(W r )] ξ=u * 2 + D * (W L , W R ) = 0, (4.6
)

with D * (W L , W R ) = ∆α 1 (u * 2 , 0, -π * 1 , 0, π * 1 , 0, 0) T and ∆α 1 = α 1,R -α 1,L . (iv) The energy jump of phase 1 across the u * 2 -discontinuity satisfies -u * 2 [α 1 ρ 1 E 1 ] ξ=u * 2 + [α 1 ρ 1 E 1 u 1 + α 1 π 1 u 1 ] ξ=u * 2 -u * 2 π * 1 ∆α 1 ≤ 0. (4.7)
If (4.7) is a strict inequality, the solution is said to be energy-dissipative. Otherwise the solution is energy-preserving.

(v) 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 . (4.8)
Some comments on Definition 4.1:

1. In the above definition, u * 2 is the the propagation velocity of the phase fraction wave. 2. The non-conservative product π 1 ∂ x α 1 identifies with the Dirac measure π * 1 ∆α 1 δ x-u * 2 t . 3. Item (ii) states that, away from the u 2 -wave, the system behaves as two independent relaxation systems for the isentropic Euler equations.

4. Items (iii) and (iv) define the re-coupling relations between the two phases at the void fraction wave u 2 . In particular, summing the third and fifth components of (4.6) yields the conservation of the total momentum. Note that since u 2 remains constant across this wave (4.5), the fourth (mass conservation of phase 2), and seventh (convection of α 2 ρ 2 T 2 ) components of (4.6) are automatically satisfied.

5. Inequality (4.7) expresses that the phase 1 energy is either conserved (when (4.7) is an equality) or dissipated (when (4.7) is a strict inequality) through the u 2 -contact discontinuity. As for the phase 2 energy, since u 2 is constant through this wave (4.5), and by the fifth component of (4.6), it is exactly conserved in the following sense:

-u * 2 [α 2 ρ 2 E 2 ] ξ=u * 2 + [α 2 ρ 2 E 2 u 2 + α 2 π 2 u 2 ] ξ=u * 2 + u * 2 π * 1 ∆α 1 = 0. (4.9)
Of course, speaking of the "conservation" of the phasic energies in an abuse of the usual terminology since neither of the two phasic energies has a conservative equation. Nevertheless, in the whole paper, we use this term when the phasic energy E i satisfies its natural equation (3.9) in the weak sense across the u 2 -contact discontinuity. Summing (4.7) and (4.9), we see that the total mixture energy is either conserved as in (3.17) or dissipated across this discontinuity. In that, the conservation (resp. dissipation) of the phase 1 energy is equivalent to the conservation (resp. dissipation) of the total mixture energy across the u 2 -field. As this field is linearly degenerate, the total energy is actually expected to be preserved when the system is hyperbolic. In the sequel indeed, we will see that in most cases, preserving the total energy (3.17) through the u 2 -contact discontinuity is possible for constructing admissible solutions of the Riemann problem. In this case, one may speak of energy-preserving solutions. However, it appears that when the ratio

α 1,L α 1,R
(or its inverse) is large, constructing solutions with positive densities while maintaining the exact energy conservation (3.17) across the u 2 -contact discontinuity is impossible. It will be shown that in these vanishing phase regimes, solutions with positive densities cannot be obtained unless one authorizes some dissipation of the total energy through the u 2 -contact, and one speaks of energy-dissipative solutions in that particular case.

6. For the applications aimed at by 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 void fraction α 1 lies in-between the acoustic waves of phase 1 which is what is required in item (v). Observe that this requirement prevents the loss of hyperbolicity due to wave interactions. In the sequel, these solutions are classified in three categories depending on the ordering between the u 1 -contact discontinuity, and the u 2 -contact discontinuity.

x t u1 + a1τ1 u2 + a2τ2 u * 1 u * 2 u2 -a2τ2 u1 -a1τ1 A solution with u * 2 < u * 1 . x t u1 + a1τ1 u2 + a2τ2 u * 2 = u * 1 u2 -a2τ2 u1 -a1τ1 A solution with u * 2 = u * 1 . x t u1 + a1τ1 u2 + a2τ2 u * 2 u * 1 u1 -a1τ1 u2 -a2τ2 A solution with u * 2 > u * 1 .
7. The considered solutions are allowed to have phasic supersonic speeds |u i | > a i τ i . 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 .

4.2

The resolution strategy: a fixed-point procedure

The key challenge to solving the Riemann problem for the relaxation system consists in determining this non-conservative product π * 1 ∆α 1 δ x-u * 2 t in the case α 1,L = α 1,R , i.e. in determining the values of u * 2 and π * 1 satisfying the constraints of Definition 4.1. Indeed, if α 1,L = α 1,R , the two phases are decoupled and the resolution is straightforward as stated in the following lemma.

1,L = α 1,R . Then, a self-similar solution ξ → W r (ξ; W L , W R ) is such that α 1 (ξ) = cst = α 1,L = α 1,R
, so that the non-conservative product π 1 (τ 1 , T 1 )∂ x α i vanishes. As a consequence, the evolutions of the two phases are completely decoupled and the intermediate states for each phase are given for i ∈ {1, 2} by

u i := 1 2 (u i,L + u i,R ) - 1 2a i (π i,R -π i,L ), (4.10)
π i := 1 2 (π i,R + π i,L ) - a i 2 (u i,R -u i,L ), (4.11)
τ i,L := τ i,L + 1 a i (u i -u i,L ) = τ i,L + 1 2a i (u i,R -u i,L ) - 1 2a 2 i (π i,R -π i,L ), (4.12
)

τ i,R := τ i,R - 1 a i (u i -u i,R ) = τ i,R + 1 2a i (u i,R -u i,L ) + 1 2a 2 i (π i,R -π i,L ). (4.13) x t u1,R + a1τ1,R u 1 u1,L -a1τ1,L τ1,L, u1,L, π1,L τ 1,L u 1 π 1 τ 1,R u 1 π 1 τ1,R, u1,R, π1,R x t u2,R + a2τ2,R u 2 u2,L -a2τ2,L τ2,L, u2,L, π2,L τ 2,L u 2 π 2 τ 2,R u 2 π 2 τ2,R, u2,R, π2,R
Phase solutions in the case of constant initial phase fractions αi,L = αi,R. In each case, T i,L = Ti,L and T i,R = Ti,R.

The value of u * 2 is given by u 2 and the solution has a subsonic wave ordering if and only if, these quantities satisfy the following constraint:

u 1 -a 1 τ 1,L < u 2 < u 1 + a 1 τ 1,R ⇐⇒ -a 1 τ 1,R < u 1 -u 2 < a 1 τ 1,L . (4.14)
Observe that the quantities defined in (4.10) to (4.13) are independent of the phase fractions α 1,L = α 1,R . On the contrary, if α i,L = α i,R , the evolutions of both phases are affected by the u 2wave (which has multiplicity 2) and the physical quantities of the two phases are coupled through this wave.

x t u1,R + a1τ1,R u * 1 u * 2 u1,L -a1τ1,L x t u2,R + a2τ2,R u * 2 u2,L -a2τ2,L
Phase solutions in the case of non-constant initial phase fractions αi,L = αi,R.

Starting from the known solution in the decoupled case |α 1,L -α 1,R | = 0, we seek to construct a branch of solutions with subsonic wave ordering, in the non-conservative cases |α 1,L -α 1,R | = 0. Actually the aim is to expose a subsonic type condition, similar to (4.14) which accounts for the subsonic ordering requirement

u 1 -a 1 τ 1,L < u * 2 < u 1 + a 1 τ 1,R , (4.15)
and ensures the existence of such a solution. Of course, the main difficulty here is that the value of u * 2 is not a priori known with respect to the initial data unlike in the case α 1,L = α 1,R . However, the analysis carried out in this paper will expose a very simple generalization of condition (4.14) valid for the case α 1,L = α 1,R and that can be explicitly tested with respect to the initial data just as (4.14).

For this purpose, we recall the following key remark, made for instance by Ambroso et al. [START_REF] Ambroso | Relaxation and numerical approximation of a two-fluid two-pressure diphasic model[END_REF] and which is the cornerstone of the whole resolution strategy.

Key remark: Consider the case α 1,L = α 1,R . If one is able to make a prediction of the pressure π * 1 that defines the non-conservative product π 1 ∂ x α i and therefore shift it to a known right hand side of the system, then one can see that the governing equations for phase 2 are completely independent of the phase 1 quantities, namely ρ 1 , u 1 and T 1 .

In [START_REF] Ambroso | Relaxation and numerical approximation of a two-fluid two-pressure diphasic model[END_REF], the authors calculate an approximate solution of the Riemann problem (4.1)-(4.2) by actually formulating a prediction π * 1 (W L , W R ) of the non-conservative product. Our approach is different. It consists in performing a fixed-point procedure in order to determine the "right" value of π * 1 and thus calculate an exact solution to the Riemann problem (4.1)-(4.2). Formally, the fixed-point procedure consists in iterating on the pair (u * 2 , π * 1 ) by alternately considering each one of the two phases, as described hereunder.

First step: The pressure π * 1 defining the non-conservative product

π 1 ∂ x α 1 = π * 1 ∂ x α 1 = -π * 1 ∂
x α 2 is first assumed to be known. Hence, the governing equations for phase 2 are completely independent of the phase 1 quantities, and one determines u * 2 by solving the following Riemann problem for the phase 2 system:

       ∂ t α 2 + u 2 ∂ x α 2 = 0, ∂ t (α 2 ρ 2 ) + ∂ x (α 2 ρ 2 u 2 ) = 0, ∂ t (α 2 ρ 2 u 2 ) + ∂ x (α 2 ρ 2 u 2 2 + α 2 π 2 (τ 2 , T 2 )) -π * 1 ∂ x α 2 = 0, ∂ t (α 2 ρ 2 T 2 ) + ∂ x (α 2 ρ 2 T 2 u 2 ) = 0, (4.16)
with the initial condition

W 2 (x, t = 0) = W 2,L if x < 0, W 2,R if x > 0, (4.17)
where W 2 = (α 2 , α 2 ρ 2 , α 2 ρ 2 u 2 , α 2 ρ 2 T 2 ) T denotes the state vector for phase 2 and (W 2,L , W 2,R ) are the restriction of the complete initial data (W L , W R ) to the phase 2 variables. The sought solutions of (4.16)-(4.17) are asked to obey an additional equation on the energy which reads:

∂ t (α 2 ρ 2 E 2 ) + ∂ x (α 2 ρ 2 E 2 + α 2 π 2 ) u 2 -u 2 π * 1 ∂ x α 2 = 0. (4.18)
Solving this Riemann problem enables to calculate the value u * 2 of the phase fraction propagation speed with respect to the assumed value of π * 1 , and thus to define a function

F [W L , W R ; a 2 ] : R -→ R π * 1 -→ u * 2 . (4.19)
Second step: The advection velocity u * 2 of the phase fraction α 1 is then assumed to be known. As a consequence, the evolution of the phase 1 quantities is independent of phase 2, and computing them amounts to solving the following phase 1 Riemann problem:

       ∂ 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 π 1 (τ 1 , T 1 )) -π 1 ∂ x α 1 = 0, ∂ t (α 1 ρ 1 T 1 ) + ∂ x (α 1 ρ 1 T 1 u 1 ) = 0, (4.20)
with the initial condition

W 1 (x, t = 0) = W 1,L if x < 0, W 1,R if x > 0, (4.21)
where

W 1 = (α 1 , α 1 ρ 1 , α 1 ρ 1 u 1 , α 1 ρ 1 T 1 )
T denotes the state vector for phase 1, and (W 1,L , W 1,R ) are the restriction of the complete initial data (W L , W R ) to the phase 1 variables. Here, in accordance with (4.7), the sought solutions are asked to obey the following inequality on the phase 1 energy:

∂ t (α 1 ρ 1 E 1 ) + ∂ x (α 1 ρ 1 E 1 + α 1 π 1 ) u 1 -u * 2 π 1 ∂ x α 1 ≤ 0. (4.22)
Once the Riemann problem for (4.20)-(4.21) is solved, applying Rankine-Hugoniot's jump relation to the momentum equation for the phase fraction wave of speed u * 2 , allows to compute the weight π * 1 of the non-conservative product π 1 ∂ x α 1 . Hence, this second step enables to define a function

G [W L , W R ; a 1 ] : R -→ R u * 2 -→ π * 1 .
(4.23)

Fixed-point: Performing an iterative procedure on these two steps actually boils down to the following fixed-point research.

Find u * 2 in (u 1 -a 1 τ 1,L , u 1 + a 1 τ 1,R ) (u 2 -a 1 τ 2,L , u 2 + a 2 τ 2,R ) such that u * 2 = F [W L , W R ; a 2 ] • G [W L , W R ; a 1 ] (u * 2 ). (4.24)
The interval where u * 2 must be sought corresponds to the subsonic wave ordering condition (4.8) on the one hand, and to the positivity of the intermediate states of phase 2 on the other hand (see Proposition (4.4)).

Section 4.3 is devoted to presenting and commenting the main result of the paper, which is an existence theorem for the Riemann problem (4.1)-(4.2). In section 4.4, the first step of the iterative process is performed and we give the explicit formula of function F [W L , W R ; a 2 ] defined in (4.19). In section 4.4, the second step of the iterative procedure is performed. We restrict the presentation to the research of solutions with the wave configurations u * 2 < u * 1 . The results for the other wave configurations u * 2 = u * 1 and u * 2 > u * 1 can be obtained through the same process, or can be inferred from the Galilean invariance of the equations. This second step allows us to define an explicit formula for function G [W L , W R ; a 1 ] introduced in (4.23). Finally, in section 4.6, we prove that, under some explicitly formulated assumptions on the initial data (W L , W R ) (see Theorem 4.2), there exists a unique energy-preserving solution to the fixed-point problem (4.24). This unique solution corresponds to the exact conservation of the total energy across the u 2 -contact discontinuity. It is shown however that in some cases where the ratio α 1,L α 1,R is large, this solution may have non-positive densities. By defining an appropriate kinetic relation which amounts to relaxing the conservation of the total energy, we recover the existence of positive solutions.

An existence theorem for solutions with subsonic wave ordering

We may now state the existence theorem for the Riemann problem (4.1)-(4.2). We refer to equations (4.10) to (4.13) for the definition of the quantities used in the theorem, an we define the following number which solely depends on the initial phase fractions:

Λ α := α 2,R -α 2,L α 2,R + α 2,L . (4.25)
Theorem 4.2. Let be given a pair of admissible initial states (W L , W R ) ∈ Ω r × Ω r and assume that the parameter a i is such that τ i,L > 0 and τ i,R > 0 for i in {1, 2}. There exists solutions with subsonic wave ordering to the Riemann problem (4.1)-(4.2) in the sense of Definition 4.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 .
Moreover, condition (A) can be decomposed into the three following conditions defining the wave configuration:

(i) Either

(A1) 0 < u 1 -u 2 -1 a2 Λ α (π 1 -π 2 ) 1 + a1 a2 |Λ α | < a 1 τ 1,L , -
and the solutions have the wave configuration

u * 2 < u * 1 . (ii) Or (A2) -a 1 τ 1,R < u 1 -u 2 -1 a2 Λ α (π 1 -π 2 ) 1 + a1 a2 |Λ α | < 0,
and the solutions have the wave configuration

u * 2 > u * 1 . (iii) Or (A3) u 1 -u 2 - 1 a 2 Λ α (π 1 -π 2 ) = 0,
and the solutions have the wave configuration

u * 2 = u * 1 .
The proof of this theorem follows from the steps described in the three following sections 4.4, 4.5 and 4.6. Before giving the details of these steps, let us first make some comments on this result. Some comments on Theorem 4.2:

1. Assumption (A) (actually (A1), (A2) or (A3)) can be explicitly tested in terms of the initial data and the parameters a i , i ∈ {1, 2}. Of course, there is no similar result concerning the Riemann problem for the isentropic equilibrium Baer-Nunziato system (2.1).

2. Assumption (A) reduces to (4.14) when α 1,L = α 1,R since in this case Λ α = 0. In this sense, assumption (A) is a generalization of (4.14) for the non-conservative case α 1,L = α 1,R .

3. 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 a 2 Λ α (π 1 -π 2 ) 1+ a 1 a 2 |Λ α |
, 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 .

4. One may formulate a more geometrical interpretation of Theorem 4.2. Assuming that there exists a solution with subsonic relative speeds when |α 1,R -α 1,L | = 0 (i.e. assuming (4.14)), the theorem shows that if |α 1,R -α 1,L | = 0 is sufficiently small, then the Riemann problem still admits energy-preserving subsonic solutions. Provided that one allows some energy-dissipation across the u 2 -wave, this branch of solutions can be followed for Riemann problems in which |α 1,R -α 1,L | increases (holding the other quantities in the initial left and right data fixed) until assumption

(A) is violated, or until |α 1,R -α 1,L | = 1.

Positivity of phase 1 densities. If the ratio

α 1,L α 1,R is in a neighborhood of 1, condition (
A) is a necessary and sufficient condition for the existence of a unique energy-preserving solution. If

α 1,L
α 1,R is too large (under (A1)), or too small (under (A2)), depending on the wave ordering, ensuring positive densities for phase 1 may require 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 < 0. (4.26)
This, combined with the energy conservation in phase 2 (4.18) (see Proposition 4.4), implies that the total mixture energy is dissipated. In sections 4.5.2 and 4.6.2, we propose a kinetic relation for the determination of one solution, among all the admissible dissipative solutions given by the theorem.

6. Positivity of phase 2 densities. Assumption (A) allows to compute the value of the wave propagation velocity u * 2 (see section 4.6). 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 . (4.27)
In the numerical applications using this Riemann solver (see sections 5 and 6), it will always be possible to ensure property (B) by taking a large enough value of the relaxation parameter a 2 .

First step of the fixed-point procedure

In this first step, the pressure π * 1 defining the non-conservative product

π 1 ∂ x α 1 = π * 1 ∂ x α 1 = -π * 1 ∂
x α 2 is assumed to be known, and one solves the Riemann problem (4.16)-(4.17) in order to compute u * 2 , the propagation speed of the phase fraction α 2 . Observe that the non-conservative product π * 1 ∂ x α 2 is not ambiguous here since π * 1 is a known constant. The following proposition characterizes the convective behavior of system (4.16).

Proposition 4.3. System (4.16) is a hyperbolic system with linearly degenerate fields associated with the eigenvalues u 2 -a 2 τ 2 , u 2 and u 2 + a 2 τ 2 . The eigenvalue u 2 has multiplicity 2.

Proof . The proof is left to the reader.

We have the following well-posedness result for the governing equations of phase 2: Proposition 4.4. Assume that the parameter a 2 is such that τ 2,L > 0 and τ 2,R > 0. Then the Riemann problem (4.16)-(4.17) admits a unique solution whose intermediate states are defined by:

x t u2 + a2τ2 u * 2 u2 -a2τ2 W2,L τ2,L * ,T2,L * u2,L * τ2,R * ,T2,R * u2,R * W2,R τ 2,L * = τ 2,L + ∆α 1 a 2 2 π 2 -π * 1 α 2,L + α 2,R , u 2,L * = u * 2 = u 2 + ∆α 1 a 2 π 2 -π * 1 α 2,L + α 2,R , T 2,L * = T 2,L , (4.28)
τ 2,R * = τ 2,R - ∆α 1 a 2 2 π 2 -π * 1 α 2,L + α 2,R , u 2,R * = u * 2 , T 2,R * = T 2,R . (4.29)
The intermediate densities ρ 2,L * and ρ 2,R * are positive if and only if

u 2 -a 2 τ 2,L < u * 2 < u 2 + a 2 τ 2,R . (4.30)
Moreover, this unique solution satisfies the following energy equation in the usual weak sense:

∂ t (α 2 ρ 2 E 2 ) + ∂ x (α 2 ρ 2 E 2 + α 2 π 2 ) u 2 -u 2 π * 1 ∂ x α 2 = 0. (4.31)
Proof . We only sketch the proof. The expressions of the intermediate states directly follow from classical manipulations of Rankine-Hugoniot's jump relations. The only non classical relation is the jump relation across the u 2 -wave for the momentum equation, where the non-conservative product is taken into account:

-u * 2 [α 2 ρ 2 u 2 ] + α 2 ρ 2 u 2 2 =0 + [α 2 π 2 ] = π * 1 ∆α 2 = -π * 1 ∆α 1 . ⇐⇒ [α 2 π 2 ] = π * 1 ∆α 2 = -π * 1 ∆α 1 .
(4.32)

The densities ρ 2,L * and ρ 2,R * are positive if and only if u * 2 -a 2 τ 2,L * < u * 2 < u * 2 + a 2 τ 2,R * . As the fields are linearly degenerate, the corresponding eigenvalues are Riemann invariants and we have u * 2 -a 2 τ 2,L * = u 2,L -a 2 τ 2,L = u 2 -a 2 τ 2,L . In the same way, u * 2 +a 2 τ 2,R * = u 2,R +a 2 τ 2,R = u 2 +a 2 τ 2,R . For the energy equation (4.31), the proof consists in verifying that the associated Rankine-Hugoniot jump relation is satisfied for all of the three waves u 2 -a 2 τ 2 , u 2 and u 2 + a 2 τ 2 . Across the two extreme waves where ∂ x α 2 vanishes, the result directly follows from the decoupling of the two phases into two (relaxed) Euler systems for which the result is well-known (see for instance [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF]). As for the u 2 -wave, thanks to (4.32), the energy jump relation reads

-u * 2 [α 2 ρ 2 E 2 ] ξ=u * 2 + [α 2 ρ 2 E 2 u 2 ] ξ=u * 2 =0 + [α 2 π 2 u 2 ] ξ=u * 2 = u * 2 [α 2 π 2 ] = u * 2 π * 1 ∆α 2 = -u * 2 π * 1 ∆α 1 .
Remark

Second step of the fixed-point procedure

In this step, the velocity u * 2 of the wave supporting the α 1 discontinuity is assumed to be known, while the pressure π * 1 defining the non-conservative product

π 1 ∂ x α 1 = π * 1 ∆α 1 δ x-u * 2 t
is an unknown that must be calculated by solving the Riemann problem (4.20)-(4.21) and applying Rankine-Hugoniot's jump relation to the momentum equation.

A convenient change of variables

In order to solve the Riemann problem (4.20)-(4.21), it is judicious to rewrite the equations in the frame moving at the known constant speed u * 2 . For this purpose, we perform the following change of variables: (x, t) → (y, t) = (x -u * 2 t, t). Any function W of the variables (x, t), is associated with a function W of the variables (y, t) such that W(y, t) = W(x, t) i.e. W(x, t) = W(x -u * 2 t, t). Denoting w 1 = u 1 -u * 2 the fluid velocity of phase 1 in the frame of the u * 2 -wave, system (4.20) rewrites

       ∂ t α 1 = 0, ∂ t (α 1 ρ 1 ) + ∂ y (α 1 ρ 1 w 1 ) = 0, ∂ t (α 1 ρ 1 w 1 ) + ∂ y (α 1 ρ 1 w 2 1 + α 1 π 1 (τ 1 , T 1 )) -π 1 ∂ y α 1 = 0, ∂ t (α 1 ρ 1 T 1 ) + ∂ y (α 1 ρ 1 T 1 w 1 ) = 0, (4.33)
while the initial conditions (4.21) become, in the non-conservative variables:

W 1 (x, t = 0) = W 1,L = (α 1,L , ρ 1,L , w 1,L = u 1,L -u * 2 , T 1,L ), y < 0, W 1,R = (α 1,R , ρ 1,R , w 1,R = u 1,R -u * 2 , T 1,R ), y > 0. (4.34)
The solutions of (4.33)-(4.34) are asked to obey the following conservative energy inequality obtained by applying the same change of variables to the energy inequality (4.22):

∂ t α 1 ρ 1 E 1 + ∂ y α 1 ρ 1 E 1 w 1 + α 1 π 1 w 1 ≤ 0, (4.35)
where

E 1 = w 2 1 2 + e 1 (T 1 ) + 1 2a 2 1 (π 2 1 (τ 1 , T 1 ) -P 2 1 (T 1 )
). A solution which satisfies (4.35) with a strict inequality is an energy-dissipative solution. In the equality case, the solution is said to be energypreserving. System (4.33) is nothing but the relaxation system introduced in [START_REF] Coquel | A Robust and Entropy-Satisfying Numerical Scheme for Fluid Flows in Discontinuous Nozzles[END_REF] for the approximation of nozzle flows, and for which the associated Riemann problem has been fully resolved. Hence, we actually calculate a solution W(y, t) of the Riemann problem associated with system (4.33), and the solution for the original Riemann problem (4.20)-(4.21) is obtained by W(x, t) = W(x -u * 2 t, t), and by adding u * 2 to the velocities w 1 . The following proposition characterizes the convective behavior of system (4.33). Proposition 4.5. System (4.33) admits four real eigenvalues that are w 1 -a 1 τ 1 , w 1 , w 1 + a 1 τ 1 and 0. All the fields are linearly degenerate and the system is hyperbolic if, and only if |w 1 | = a 1 τ 1 .

Proof . The proof can be easily recovered following closely related steps developed in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] for the Suliciu relaxation system. The details are left to the reader. The u * 2 -wave in (4.20) corresponds to the standing wave in (4.33) and the subsonic condition on the relative speed (4.8) now reads

w 1,L -a 1,L τ 1,L < 0 < w 1,R + a 1,R τ 1,R . (4.36)
For the resolution, we only consider solutions with the subsonic wave ordering u * 2 < u * 1 since the other possible wave orderings can be obtained by the Galilean invariance of the equations. This ordering for (4.20) corresponds to the wave ordering w 1 -a 1 τ 1 < 0 < w 1 < w 1 + a 1 τ 1 for (4.33), which in [START_REF] Coquel | A Robust and Entropy-Satisfying Numerical Scheme for Fluid Flows in Discontinuous Nozzles[END_REF] is referred to as the < 1, 2 > wave configuration:

x t u1 -a1τ1 u * 2 u1 u1 + a1τ1 W1,L W - 1 W + 1 W1,R * W1,R y t w1 -a1τ1 0 w1 w1 + a1τ1 W1,L W - 1 W + 1 W1,R * W1,R ⇐⇒
Following [START_REF] Coquel | A Robust and Entropy-Satisfying Numerical Scheme for Fluid Flows in Discontinuous Nozzles[END_REF], we introduce the following notations, which are used for the definition of the solutions to the Riemann problem (4.33)-(4.34).

ν = α 1,L α 1,R , w = u 1 -u * 2 , M * L = w a 1 τ 1,L , ω = 1 -M * L 1 + M * L . (4.37)
where the quantites are defined in (4.10)-(4.11)-(4.12)-(4.13). We may now recall the main result (see Proposition (3.4) of [START_REF] Coquel | A Robust and Entropy-Satisfying Numerical Scheme for Fluid Flows in Discontinuous Nozzles[END_REF]) which states that one can build a one-parameter family of solutions with the subsonic wave ordering w 1 -a 1 τ 1 < 0 < w 1 < w 1 + a 1 τ 1 for the Riemann problem (4.33)-(4.34) (i.e. the wave ordering u 1 -a 1 τ 1 < u * 2 < u * 1 < u 1 + a 1 τ 1 for the Riemann problem (4.20)-(4.21)), and the dissipation of energy across the standing wave is directly driven by the underlying parameter. Proposition 4.6. Assume that a 1 is such that τ 1,L > 0 and τ 1,R > 0. Then the Riemann problem (4.33)-(4.34) admits solutions with the subsonic wave ordering w

1 -a 1 τ 1 < 0 < w 1 < w 1 + a 1 τ 1 , if and only if 0 < M * L < 1. (4.38)
These solutions can be parametrized by

M := M -= w - 1 a1τ - 1
, the Mach number of the state on the left of the standing wave, and the intermediate states are given by:

τ - 1 = τ 1,L 1 -M * L 1 -M , w - 1 = a 1 Mτ - 1 , T - 1 = T 1,L , (4.39)
τ + 1 = τ 1,L 1 + M * L 1 + νM , w + 1 = νa 1 Mτ + 1 , T + 1 = T 1,L , (4.40
)

τ 1,R * = τ 1,R + τ 1,L M * L -νM 1 + νM , w 1,R * = νa 1 Mτ + 1 , T 1,R * = T 1,R . (4.41)
Besides, there exists a critical value ν in (1, +∞] independent of (α 1,L , α 1,R ) and possibly infinite such that the following alternative holds.

• If ν < ν , M belongs to the interval (0, M 0 (ν, ω)] (0, min(1, 1/ν)) with

M 0 (ν, ω) = 1 2   1 + ω 2 1 -ω 2 1 + 1 ν - 1 + ω 2 1 -ω 2 2 1 + 1 ν 2 - 4 ν   . (4.42)
The value M = M 0 (ν, ω) gives the unique energy-preserving solution and for 0 < M < M 0 (ν, ω), the solution is energy-dissipative.

• If ν ≥ ν , no energy-preserving solution has positive densities. The initial data is such that 0 < M * L ν < M 0 (ν, ω) < min(1, 1/ν) where M 0 (ν, ω) is given by (4.42). M must be strictly less than M 0 (ν, ω), and by taking M close enough to M * L ν it is always possible to ensure that all the densities remain positive.

In both cases, the choice of the value of M determines the mass, momentum and energy jumps across the standing wave through

[α 1 ρ 1 w 1 ] y t =0 = 0, (4.43
)

α 1 ρ 1 w 2 1 + α 1 π 1 y t =0 = π 1 ∆α 1 -a 2 1 ((α 1,R + α 1,L )M * L -2α 1,L M) τ 1,L , (4.44
)

α 1 ρ 1 E 1 w 1 + α 1 π 1 w 1 y t =0 = 1 2 (w + a 1 τ 1,L ) 2 Q 0 (M) (1 + νM)(1 -M) ϕ (M; ν, ω) , (4.45)
where Proof . The proof is mainly given in [START_REF] Coquel | A Robust and Entropy-Satisfying Numerical Scheme for Fluid Flows in Discontinuous Nozzles[END_REF] with slightly different notations. Indeed, it is easy to see that condition (3.20) in [START_REF] Coquel | A Robust and Entropy-Satisfying Numerical Scheme for Fluid Flows in Discontinuous Nozzles[END_REF] is equivalent to "0 < M L < 1" which here corresponds to 0 < M * L < 1. The proof consists in proving that 0 < M * L < 1 is a necessary and sufficient condition for the existence of solutions M in the interval (0, min(1, 1/ν)) to the second order polynomial in-equation ϕ (M; ν, ω) ≤ 0, thus enforcing a non-positive energy jump across the standing wave according to (4.45). The value M 0 (ν, ω) given in (4.42) is the unique solution of ϕ (M; ν, ω) = 0 that lies in the interval (0, min(1, 1/ν)). As explained in [START_REF] Coquel | A Robust and Entropy-Satisfying Numerical Scheme for Fluid Flows in Discontinuous Nozzles[END_REF], the existence of ν is related to the expression of τ 1,R * in (4.41) which is the only intermediate specific volume that may be non-positive. It is possible to show that the function

Q 0 (M) = α 1,L ρ - 1 w - 1 = α 1,R ρ + 1 w + 1 >
ν → τ 1,R * (ν) = τ 1,R -τ 1,L νM 0 (ν, ω) -M * L 1 + νM 0 (ν, ω) , (4.47)
is a non-increasing function that may become negative for large values of ν. Observe that for ν = 1, we have νM 0 (ν, ω) = M * L which implies that the pathological values of ν are larger than one (i.e. ν > 1). In such pathological cases, in order to impose the positivity of τ 1,R * we must no longer exactly conserve the energy at the standing wave (by taking M = M 0 (ν, ω)) but dissipate it by taking M smaller than M 0 (ν, ω). Indeed, ϕ (M; ν, ω) ≤ 0 for all M ∈ (0, M 0 (ν, ω)]. The expression of τ 1,R * clearly shows that if M is taken close enough to M * L ν (remember that ν > 1), we have τ 1,R * close to τ 1,R which is positive by hypothesis.

Let us however prove (4.44) which is not proved in [START_REF] Coquel | A Robust and Entropy-Satisfying Numerical Scheme for Fluid Flows in Discontinuous Nozzles[END_REF]. As T 1 is constant across the standing wave and is equal to T 1,L for both states W - 1 and W + 1 , one may write

α 1 ρ 1 w 2 1 + α 1 π 1 y t =0 = (α + 1 ρ + 1 w + 1 2 -α - 1 ρ - 1 w - 1 2 ) -a 2 1 (α + 1 τ + 1 -α - 1 τ - 1 ) + (P 1 (T 1,L ) + a 2 1 T 1,L )∆α 1 , Moreover, α - 1 = α 1,L and α + 1 = α 1,R . Thus α 1 ρ 1 w 2 1 + α 1 π 1 y t =0 = a 2 1 α 1,R w + 1 2 a 2 1 τ + 1 2 τ + 1 -α 1,L w - 1 2 a 2 1 τ - 1 2 τ - 1 -a 2 1 (α 1,R τ + 1 -α 1,L τ - 1 ) +(P 1 (T 1,L ) + a 2 1 T 1,L )∆α 1 = a 2 1 α 1,R ((νM) 2 -1)τ + 1 -α 1,L (M 2 -1)τ - 1 + (P 1 (T 1,L ) + a 2 1 T 1,L )∆α 1 .
With the expressions of ν (4.37), τ - 1 (4.39), and τ + 1 (4.40) we get

α 1 ρ 1 w 2 1 + α 1 π 1 y t =0 = -a 2 1 (1 -(νM) 2 ) τ + 1 τ 1,L -ν(1 -M 2 ) τ - 1 τ 1,L α 1,R τ 1,L + (P 1 (T 1,L ) + a 2 1 T 1,L )∆α 1 = -a 2 1 ((1 -νM)(1 + M * L ) -ν(1 + M)(1 -M * L )) α 1,R τ 1,L + (P 1 (T 1,L ) + a 2 1 T 1,L )∆α 1 =(P 1 (T 1,L ) + a 2 1 T 1,L )∆α 1 -a 2 1 (1 -ν + (1 + ν)M * L -2νM) α 1,R τ 1,L =(P 1 (T 1,L ) + a 2 1 (T 1,L -τ 1,L ))∆α 1 -a 2 1 ((1 + ν)M * L -2νM) α 1,R τ 1,L , with P 1 (T 1,L ) + a 2 1 (T 1,L -τ 1,L ) = π 1 which yields (4.44).

Kinetic relation

As seen in Proposition 4.6, the admissible solutions are parametrized by a real number M lying in the interval (0, M 0 (ν, ω)]. The aim of this section is to propose an "instruction manual" for the practical choice of the parameter M in this interval, with respect to the data of the problem (W 1,L , W 1,R , u * 2 ). The rule consists in enforcing a lower bound on τ 1,R * , the only specific volume that may be nonpositive:

τ 1,R * ≥ µτ 1,R , (4.48)
where µ is a fixed real number in (0, 1). If the energy-preserving value M = M 0 (ν, ω) is such that (4.48) is satisfied, then one chooses this value. Otherwise, one chooses the value of M which verifies

τ 1,R * = µτ 1,R . As the function M → τ 1,R * (M) = τ 1,R -τ 1,L νM-M * L 1+νM
, is a non-increasing function, this amounts to take M = M(ν, M * L ) where

M(ν, M * L ) := min (M 0 (ν, ω) , M µ (ν, M * L )) , ω = 1 -M * L 1 + M * L , (4.49)
with

M µ (ν, M * L ) := 1 ν M * L + (1 -µ) τ 1,R τ 1,L 1 -(1 -µ) τ 1,R τ 1,L . (4.50) Note that if ν ≤ 1, then M 0 (ν, ω) < M µ (ν, M * L )
, which means that no dissipative correction is added since the lower-bound on τ 1,R * is already satisfied by the energy-preserving choice M 0 (ν, ω). In this section, we prove that condition (A1) is a necessary and sufficient condition for the existence of solutions to the following fixed-point problem :

Jump relations in the original frame

Defining

π * 1 ∆α 1 := π 1 ∆α 1 -a 2 1 ((α 1,R + α 1,L )M * L -2α 1,L M) τ 1,L , (4.51) Q(u * 2 , W L , W R ) := - 1 2 (w + a 1 τ 1,L ) 2 Q 0 (M) (1 + νM)(1 -M) ϕ (M; ν, ω) , ( 4
-u * 2 [α 1 ρ 1 ] x t =u * 2 + [α 1 ρ 1 ] x t =u * 2 = 0, (4.53) -u * 2 [α 1 ρ 1 u 1 ] x t =u * 2 + α 1 ρ 1 u 2 1 + α 1 π 1 x t =u * 2 -π * 1 ∆α 1 = 0, (4.54) -u * 2 [α 1 ρ 1 E 1 ] x t =u * 2 + [α 1 ρ 1 E 1 u 1 + α 1 π 1 u 1 ] x t =u * 2 -u * 2 π * 1 ∆α 1 = -Q(u * 2 , W L , W R ) ≤ 0. ( 4
Find u * 2 in (u 1 -a 1 τ 1,L , u 1 + a 1 τ 1,R ) (u 2 -a 1 τ 2,L , u 2 + a 2 τ 2,R ) such that u * 2 = F [W L , W R ; a 2 ] • G [W L , W R ; a 1 ] (u * 2 ), with u * 2 < u * 1 , (4.56)
and therefore, for the existence of solutions to the Riemann problem (4.1)-(4.2) with the subsonic wave ordering

u 1,L -a 1 τ 1,L < u * 2 < u * 1 < u 1,R + a 1 τ 1,R .
Let us first introduce some non-dimensional numbers built on the quantities defined in (4.10)-(4.11)-(4.12)-(4.13) :

M L := u 1 -u 2 a 1 τ 1,L , P L := π 1 -π 2 a 2 1 τ 1,L . ( 4

.57)

Solving the fixed-point (4.56) amounts to re-coupling the two phases that have been decoupled for a separate resolution. We start by rewriting the expression of π * 1 ∆α 1 obtained for phase 2 in (4.28):

π * 1 ∆α 1 = ∆α 1 π 2 + a 2 (α 2,L + α 2,R ) u 2 -u * 2 = ∆α 1 π 2 + a 2 (α 2,L + α 2,R ) u 1 -u * 2 + u 2 -u 1 . (4.58)
Hence, solving the fixed-point problem (4.56) amounts to seeking u * 2 such that the two expressions of π * 1 ∆α 1 given in (4.51) and (4.58) are equal, i.e. such that

π 1 ∆α 1 -a 2 1 ((α 1,R + α 1,L )M * L -2α 1,L M) τ 1,L = π 2 ∆α 1 + a 2 (α 2,L + α 2,R )a 1 τ 1,L M * L -a 2 (α 2,L + α 2,R )(u 1 -u 2 ).
(4.59)

The energy-preserving case

We first look for solutions that exactly preserve the energy equality across the u * 2 -wave. Therefore, we take

M := M 0 α 1,L α 1,R , 1-M * L 1+M * L
, where M 0 (ν, ω) is defined in (4.42). Introducing the non-dimensional quantities M L , P L and Λ α , equation (4.59) re-writes as

M L - a 1 a 2 Λ α P L = M * L + a 1 a 2 1 α 2,L + α 2,R (α 1,R + α 1,L )M * L -2α 1,L M 0 α 1,L α 1,R , 1 -M * L 1 + M * L . (4.60) Now, considering the change of variables u * 2 -→ M * L = u 1 -u * 2 a1τ 1,L
, solving the fixed-point problem (4.56) is equivalent to finding M * L such that equation (4.60) holds. Observe that by Proposition 4.6, the solution has the subsonic wave ordering u 1 -a 1 τ 1 < u 2 < u 1 if and only if M * L belongs to (0, 1). Defining the function

Ψ 0 :    (0, 1) -→ R m -→ m + a 1 a 2 1 α 2,L + α 2,R (α 1,R + α 1,L )m -2α 1,L M 0 α 1,L α 1,R , 1-m 1+m , (4.61)
the following proposition proves that condition (A1) is equivalent to the existence of a unique solution M * L in (0, 1) to our fixed-point problem: Proposition 4.7. Function m → Ψ 0 (m) is a differentiable increasing function on the interval (0, 1), whose limit values are

lim m→0 Ψ 0 (m) = 0, lim m→1 Ψ 0 (m) = 1 + a 1 a 2 |Λ α |. (4.62)
Hence, if the following condition, which is equivalent to (A1) holds,

0 < M L - a 1 a 2 Λ α P L < 1 + a 1 a 2 |Λ α |, (4.63)
then there exists a unique M * L in (0, 1) such that

Ψ 0 (M * L ) = M L - a 1 a 2 Λ α P L . (4.64)
Proof . The function Ψ 0 is clearly differentiable on the interval (0, 1). Differentiating w.r.t m, one gets

Ψ 0 (m) = 1 + a 1 a 2 α 1,R + α 1,L α 2,L + α 2,R - a 1 a 2 2α 1,L α 2,L + α 2,R d dm M 0 α 1,L α 1,R , 1 -m 1 + m
= 1 + a 1 a 2 α 1,R + α 1,L α 2,L + α 2,R - a 1 a 2 2α 1,L α 2,L + α 2,R ∂M 0 ∂ω α 1,L α 1,R , 1 -m 1 + m • dω dm ,
where ω = 1-m 1+m . We have dω dm = -2 (1+m) 2 , hence

Ψ 0 (m) = 1 + a 1 a 2 1 α 2,L + α 2,R α 1,R + α 1,L + 2α 1,L ∂M 0 ∂ω α 1,L α 1,R , 1 -m 1 + m • 2 (1 + m) 2 .

Therefore, again denoting

ν = α 1,L α 1,R and ω = 1-m
1+m , it is sufficient for the derivative Ψ 0 (m) to be positive, that

1 + ν + ν(1 + ω) 2 ∂M 0 ∂ω (ν, ω) ≥ 0,
for all ω in (0, 1),

⇐⇒ 1 + (1 + ω) 2 ∂ ∂ω ν 1 + ν M 0 (ν, ω) ≥ 0, for all ω in (0, 1). (4.65)
With the expression of M 0 (ν, ω) in (4.42), one gets

ν 1 + ν M 0 (ν, ω) = 1 2   1 + ω 2 1 -ω 2 - 1 + ω 2 1 -ω 2 2 - 4ν (1 + ν) 2   .
Differentiating this with respect to ω yields

∂ ∂ω ν 1 + ν M 0 (ν, ω) = 1 2   4ω (1 -ω 2 ) 2 - 4ω (1 -ω 2 ) 2 1 + ω 2 1 -ω 2 1 + ω 2 1 -ω 2 2 - 4ν (1 + ν) 2 -1/2   = 2ω (1 -ω 2 ) 2   1 -1 - 4ν (1 + ν) 2 1 -ω 2 1 + ω 2 2 -1/2   . (4.66)
Casting this in (4.65), the sufficient condition for the function Ψ 0 to be increasing becomes

1 + 2ω (1 + ω) 2 (1 -ω 2 ) 2   1 -1 - 4ν (1 + ν) 2 1 -ω 2 1 + ω 2 2 -1/2   ≥ 0, for all ω in (0, 1). (4.67)
Now, isolating the terms in ν and those in ω, (4.67) is equivalent to

4ν (1 + ν) 2 ≤ 1 + ω 2 1 -ω 2 2   1 - 1 1 + (1-ω 2 ) 2 2ω(1+ω) 2 2    , for all ω in (0, 1). (4.68)
A calculation shows that the right-hand side term of (4.68) is independent of ω and equals 1. Hence, a sufficient condition for the function Ψ 0 to be strictly increasing is

4ν (1 + ν) 2 ≤ 1, (4.69)
which is true for any ν in R + . As for the limit values of Ψ 0 , observe that the function M 0 (ν, ω) is such that lim ω→0 M 0 (ν, ω) = min 1, 1 ν and lim ω→1 M 0 (ν, ω) = 0, hence the limits (4.62) as m tends to 0 and 1. Finally, Proposition 4.7 follows from the intermediate value theorem.

Thus, provided positive values of the densities, Proposition 4.7 proves that (A1) is a necessary and sufficient condition for the existence and uniqueness of an energy-preserving solution. If the phase 1 densities are not positive, one must authorize some energy dissipation as detailed hereunder.

The energy-dissipating case

According to the proof of Proposition 4.6, it may occur when the ratio ν =

M = M 0 α 1,L α 1,R , 1-M * L 1-M * L
) is non-positive. In such pathological case, the unique energy-preserving solution of the previous section is not admissible and one has to authorize energy dissipation introducing some kinetic relation. As explained in section 4.5.2, this kinetic relation is obtained by prescribing a lower-bound on τ 1,R * , namely τ 1,R * ≥ µτ 1,R where µ ∈ (0, 1), through the definition of a new value of the Mach number M = M(ν, M * L ) given by (4.49). The fixed-point research must now be performed with the new function

Ψ :    (0, 1) -→ R m -→ m + a 1 a 2 1 α 2,L + α 2,R (α 1,R + α 1,L )m -2α 1,L M α 1,L α 1,R , m . (4.70)
Observe that if m is such that

M 0 ν, 1-m 1+m ≤ M µ (ν, m) then Ψ(m) = Ψ 0 (m). In particular, if α 1,L ≤ α 1,R
, then Ψ identifies with Ψ 0 on the whole interval (0, 1). We have the following proposition which shows that, provided an appropriate choice of the parameter µ ∈ (0, 1), there still exists a unique solution M * L ∈ (0, 1) under condition (A1). Proposition 4.8. If the parameter µ ∈ (0, 1) is close enough to one, the function m → Ψ(m) is a Lipschitz-continuous increasing function on the interval (0, 1), whose limit values are

lim m→0 Ψ(m) = 0, lim m→1 Ψ(m) = 1 + a 1 a 2 |Λ α |. (4.71)
Hence, if condition (A1) holds, then there exists a unique M * L in (0, 1) such that

Ψ(M * L ) = M L - a 1 a 2 Λ α P L . (4.72)
Proof . If α 1,L ≤ α 1,R then Ψ ≡ Ψ 0 and the result follows from the energy-preserving case. Let us turn to the case α 1,L > α 1,R . As the minimum of two differentiable functions, M(ν, m) is Lipschitzcontinuous and so is Ψ. Actually, Ψ is almost everywhere differentiable on (0, 1). For the limit values of Ψ, we know from the energy-preserving case that lim

m→0 M 0 ν, 1-m 1+m = 0, lim m→1 M 0 ν, 1-m 1+m = min(1, 1/ν)
, and if µ is close enough to one, we have

M µ (ν, 0) = 1 ν (1 -µ) τ 1,R τ 1,L 1 -(1 -µ) τ 1,R τ 1,L > 0, M µ (ν, 1) = 1 ν 1 + (1 -µ) τ 1,R τ 1,L 1 -(1 -µ) τ 1,R τ 1,L > min 1, 1 ν . (4.73)
Hence, Ψ and Ψ 0 share the same limit values at 0 and 1. As for the monotony of function Ψ, we may write that for almost every m in (0, 1):

Ψ (m) ≥ min   Ψ 0 (m), 1 + a 1 a 2 1 α 2,L + α 2,R   α1,R + α 1,L -2 α 1,R 1 -(1 -µ) τ 1,R τ 1,L       . (4.74) For µ = 1, this expression gives Ψ (m) ≥ min Ψ 0 (m), 1 + a 1 a 2 |α 1,R -α 1,L | α 2,L + α 2,R since α 1,L > α 1,R . As
Ψ 0 (m) > 0 by the study of the energy-preserving case, this proves that Ψ is strictly increasing if µ is close enough to one, which concludes the proof.

Proof of Theorem 4.2 :

We may now complete the proof of Theorem 4.2. If the ratio

α 1,L α 1,
R is close to one, Proposition 4.6 concerning phase 1, asserts that no energy dissipation is needed for ensuring the positivity of the densities. Hence, for α 1,L α 1,R in a neighborhood of 1, by Proposition 4.7, condition (A1) is a necessary and sufficient condition for the existence of a unique solution to the fixed-point problem (4.60), i.e. for the existence and uniqueness of an energy-preserving solution to the Riemann problem (4.1)-(4.2) with the subsonic wave ordering u * 2 < u * 1 . For large values of the ratio

α 1,L
α 1,R , Proposition 4.6 shows that ensuring positive densities for phase 1 may require strict energy dissipation across the u 2 -wave. In this case, still assuming condition (A1), Proposition 4.8 proves that using the kinetic relation (4.49) defining M with respect to the pair

α 1,L α 1,R , M *
L , it is always possible to ensure the existence of a unique solution with positive densities for phase 1 by dissipating the total energy. Note that uniqueness is achieved once the kinetic relation is prescribed. If another choice had been made for the kinetic relation (4.49), another dissipative solution could be found.

Finally, thanks to the Galilean invariance of system (4.1), one can prove that the symmetric waveconfiguration u * 2 > u * 1 is implied by (A2) by exchanging the subscripts L and R and changing the velocities to their opposite values. As for condition (A3), it can be obtained by passing to the limit in (A1). The corresponding M * L is equal to zero, and we obtain the u * 2 = u * 1 configuration. This observation concludes the proof of Theorem 4.2.

Expression of the Riemann solution

In this section, we construct the solution W(x, t; W L , W R ) = W r (x/t; W L , W R ) for a given pair of initial conditions (W L , W R ) in Ω r and two parameters a 1 and a 2 meeting the conditions of Theorem 4.2. We distinguish the three different cases corresponding respectively to (A1), (A2) and (A3).

• If (A1) holds, the phasic solutions have the following form:

x t u1,L -a1τ1,L u * 2 u * 1 u1,R + a1τ1,R W1,L W - 1 W + 1 W1,R * W1,R x t u2,R + a2τ2,R u * 2 u2,L -a2τ2,L W2,L W2,L * W2,R * W2,R
The values u * 2 and π * 1 are calculated as follows. First use an iterative method (Newton's method or a dichotomy algorithm for instance) to compute M * L such that

Ψ(M * L ) = M L - a 1 a 2 Λ α P L . (4.75)
According to section 4.6, M * L always exists under (A1) and is unique if µ is close enough to one. We then obtain u * 2 by u *

2 = u 1 -a 1 τ 1,L M * L ,
: u - 1 = w - 1 + u * 2 , u + 1 = w + 1 + u * 2 , u 1,R * = w 1,R * + u * 2 .
• If (A2) holds, we exploit the Galilean invariance of the equations. The solution is obtained by the transformation

W r (ξ; W L , W R ) := VW r (-ξ; VW R , VW L ), (4.76)
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 1 , x 2 , -x 3 , x 4 , -x 5 , x 6 , x 7 ). (4.77)
Of course, the function W r (-ξ; VW R , VW L ) is computed through the first case, since for these new initial data (VW R , VW L ), it is condition (A1) that holds.

• If (A3) holds, u * 2 is equal to u 1 (i.e. M * L = 0). The intermediate states for phase 2 are obtained through equations (4.28) and (4.29) in Proposition 4.4, and the intermediate states for phase 1 are computed by passing to the limit as M * L goes to zero (i.e. ω → 1). We obtain the intermediate states in equations (4.39) to (4.41), for M = 0.

A positive and entropy-satisfying Finite Volume scheme

In this section, the exact Riemann solver W r (ξ; W L , W R ) build in section 4 for the relaxation system 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 system (2.1). This approximate Riemann solver is proved to satisfy important stability properties such as the preservation of the densities positivity, and a discrete entropy inequality which is a discrete counterpart of the energy inequality (2.10) satisfied by the exact weak solutions of the model.

The HLL approximate Riemann solver

In order to define the approximate Riemann solver, let us introduce the mappings M :

R 5 -→ R 7 (x k ) k=1,..,5 -→ (x 1 , x 2 , x 3 , x 4 , x 5 , x 1 , 1 -x 1 ).

(5.1)

P : R 7 -→ R 5 (x k ) k=1,..,7 -→ (x 1 , x 2 , x 3 , x 4 , x 5 ).
(5.2)

An approximate solution of the Riemann problem (2.1)-(2.11) is obtained by

U app (ξ; U L , U R ) := PW r (ξ; M (U L ), M (U R )). (5.3)
Or course, W r (ξ; M (U L ), M (U R )) depends on the parameters (a 1 , a 2 ) which must be chosen so as to satisfy some stability properties (see section 5.5 below). With this approximate Riemann solver, are associated the following two interface numerical fluxes:

F -(U L , U R ) := f(U L ) - 0 -∞ {U app (ξ; U L , U R ) -U L } dξ, (5.4)
F + (U L , U R ) := f(U R ) + +∞ 0 {U app (ξ; U L , U R ) -U R } dξ. (5.5)
The integrals in (5.4) and (5.5) are well-defined thanks to (4.3) and to the property

P • M = Id R 5 . Observe that F -(U L , U R ) = F + (U L , U R ) unless α 1,L = α 1,R
because of the non-conservative part of system (4.1). In the sequel F ± k denotes the k th component of vector F ± .

The Finite Volume relaxation scheme

We now derive a Finite Volume scheme for the approximation of the entropy weak solutions of a Cauchy problem associated with system (2.1):

∂ t U + ∂ x f(U) + c(U)∂ x U = 0, x ∈ R, t > 0, U(x, 0) = U 0 (x),
x ∈ R.

(5.6)

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 . We assume that ∆t and ∆x satisfy the CFL condition

∆t ∆x max i∈{1,2},j∈Z max |(u i -a i τ i ) n j |, |(u i + a i τ i ) n j+1 | < 1 2 .
(5.7)

The Finite Volume relaxation scheme reads

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 ) .
(5.8)

We do not precise at this point the computation of the constant a i . This discussion is postponed to section 5.5.

Remark 5.1. The numerical fluxes formula that are used in practice in the implementation are

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} , 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}
,

where D * (W L , W R ) := (α 1,R -α 1,L )(u * 2 (W L , W R ), 0, -π * 1 (W L , W R ), 0, π * 1 (W L , W R ), 0, 0) T .

Basic properties of the scheme

The following theorem states some important properties of the relaxation scheme.

Theorem 5.1. The approximate Riemann solver (5.3) has the following properties:

• Consistency. The numerical fluxes satisfy F -(U, U) = F + (U, U) = f (U).
• Positivity. Under the CFL condition (5.7), if the initial condition x → U 0 (x) is in Ω, then the values (U n j ) j∈Z,n∈N computed by the scheme are such that, 0 < α n i,j < 1, and (α i ρ i ) n j > 0, for all i in {1, 2}, j in Z and n in N.

(5.9)

• Mass conservation. The fluxes (5.4) and (5.5) satisfy

F - k (U L , U R ) = F + k for k in {2, 4}.
• Total momentum conservation. The fluxes (5.4) and (5.5) satisfy

F - 3 (U L , U R ) + F - 5 (U L , U R ) = F + 3 (U L , U R ) + F + 5 (U L , U R ).
Proof . Following [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] or [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF], the approximate Riemann solver is a Godunov type scheme where U n+1 j is the cell-average over C j of the function 1 x<xj U app (x/∆t; U n j-1 , U n j ) + 1 x≥xj U app (x/∆t; U n j , U n j+1 ). Hence, the positivity property is a direct consequence of Theorem 4.2 which ensures the positivity of the phase fractions (which are simply transported) and of the densities in the relaxation Riemann solver. For this purpose, entropy dissipation across the u 2 -contact discontinuity may be necessary for enforcing this property when the ratio α 1,L α 1,R (or its inverse) is large. The proof of the other properties involves no particular difficulties. It consists in simple verifications using the expressions of the intermediate states.

Non-linear stability

Let us now turn to the discrete entropy inequality property. In order to ease the notations, we denote by η and F η the energy-energy flux pair of the Baer-Nunziato model: η(U) :=

2 i=1 α i ρ i E i and F η (U) := 2 i=1 α i (ρ i E i + p i (ρ i ))u i .
The entropy weak solutions of (2.1) satisfy the following energy inequality in the weak sense:

∂ t η(U) + ∂ x F η (U) ≤ 0.
(5.10)

We also recall the following property which is proved in [START_REF] Coquel | Two properties of two-velocity two-pressure models for two-phase flows[END_REF][START_REF] Saleh | Analyse et Simulation Numérique par Relaxation d'Ecoulements Diphasiques Compressibles. Contribution au Traitement des Phases Evanescentes[END_REF]:

Proposition 5.2. The mapping η : Ω -→ R U -→ η(U) is non-strictly convex.
In the same way, we introduce the energy-energy flux pair for the relaxation system: η r (W) :=

2 i=1 α i ρ i E i and F r η (W) := 2 i=1 α i (ρ i E i + π i (τ i , T i ))u i .
The weak solutions of (4.1)-(4.2) in the sense of Definition 4.1 satisfy the inequality

∂ t η r (W) + ∂ x F r η (W) = -Q(u * 2 , W L , W R )δ x-u * 2 t , (5.11)
where Q(u * 2 , W L , W R ) is a non-negative number defined in (4.52). In the sequel, we prove that if a so-called sub-characteristic condition (also known as Whitham's condition, see [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]) is verified by the parameters a 1 and a 2 , then the relaxation Riemann solver satisfies a discrete entropy inequality, as defined in Theorem 5.4.

Definition 5.1. Consider (U L , U R ) ∈ Ω × Ω and let (W L , W R ) = (M (U L ), M (U R )) ∈ Ω r × Ω r
be the corresponding relaxation initial data. Let ∆x and ∆t be two space and time steps satisfying the CFL condition (5.7). Denoting τ i (ξ) = ρ -1 i (ξ) the specific volumes in the solution W r (ξ; W L , W R ) of the Riemann problem (4.1)-(4.2), the parameters (a 1 , a 2 ) are said to satisfy Whitham's condition for

(U L , U R ) if for i in {1, 2}, a 2 i > -P i (τ i (ξ)), for almost every ξ in - ∆x 2∆t , ∆x 2∆t .
(5.12)

Lemma 5.3. Take the same notations as in Definition 5.1 and denote U app (ξ) and W r (ξ) so as to ease the notation. If (a 1 , a 2 ) satisfy Whitham's condition for (U L , U R ), then the relaxation approximate Riemann solver satisfies an interface entropy inequality (see [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]) in the sense that

η( U L ) -η(U L ) + 2∆t ∆x F r η W r (0 + ) -F η (U L ) ≤ 0, (5.13) η( U R ) -η(U R ) + 2∆t ∆x F η (U R ) -F r η W r (0 + ≤ 0, (5.14)
where Proof . We only prove inequality (5.13) (the proof of (5.14) is similar). By Jensen's inequality, the convexity of the map U → η(U) (see Proposition 5.2) implies that it is sufficient to prove

U L = 2 ∆x 0 -∆x 2 U app (x/∆t)dx = 2∆t ∆x 0 -∆x 2∆t U app (ξ)dξ, (5.15)
2∆t ∆x 0 -∆x 2∆t η (U app (ξ)) dξ -η(U L ) + 2∆t ∆x F r η W r (0 + ) -F η (U L ) ≤ 0, (5.17
2∆t ∆x 0 -∆x 2∆t η r (W r (ξ))dξ -η r (W L ) + 2∆t ∆x F r η (W r (0 -)) -F r η (W L ) ≤ 0. (5.18) Now, as (W L , W R ) = (M (U L ), M (U R )) are at equilibrium, we have η r (W L ) = η(U L ) and F r η (W L ) = F η (U L ). Moreover, the Riemann solution is constructed such that F r η (W r (0 + )) -F r η (W r (0 -)) ≤ 0 (indeed, we have F r η (W r (0 + )) -F r η (W r (0 -)) = 0 unless u * 2 = 0 in which case F r η (W r (0 + )) - F r η (W r (0 -)) = -Q(u * 2 , W L , W R ) ≤ 0). Replacing in (5.18) this yields -η(U L ) + 2∆t ∆x F r η W r (0 + ) -F η (U L ) ≤ - 2∆t ∆x 0 -∆x 2∆t η r (W r (ξ)) dξ. (5.19)
Hence, a sufficient condition for (5.17) (and thus for (5.13)) to hold true is

2∆t ∆x 0 -∆x 2∆t {η (U app (ξ)) -η r (W r (ξ))} dξ ≤ 0.
(5.20)

Now, for almost every ξ in -∆x 2∆t , 0 , we have

η (U app (ξ)) -η r (W r (ξ)) = 2 i=1 (α i ρ i )(ξ) e(τ i (ξ)) -e i (T i (ξ)) - 1 2a 2 i π 2 i (τ i (ξ), T i (ξ)) -P 2 i (T i (ξ)) .
(5.21) Omitting the dependence on ξ, we have for i = 1, 2:

π 2 i (τ i , T i ) -P 2 i (T i ) = (π i (τ i , T i ) -P i (T i )) (π i (τ i , T i ) + P i (T i )) = a 2 i (T i -τ i ) 2P i (T i ) + a 2 i (T i -τ i ) = -2a 2 i e i (T i )(T i -τ i ) + a 4 i (T i -τ i ) 2 ,
since e i = -P i . Hence,

e i (τ i ) -e i (T i ) - 1 2a 2 i π 2 i (τ i , T i ) -P 2 i (T i ) = e i (τ i ) -e i (T i ) -e i (T i )(τ i -T i ) - a 2 i 2 (T i -τ i ) 2 . (5.22)
A Taylor expansion with integral remainder gives

e i (τ i ) -e i (T i ) -e i (T i )(τ i -T i ) = (T i -τ i ) 2 1 0 e i (sτ i + (1 -s)T i )(1 -s)ds. (5.23)
Then, replacing in (5.22) and observing that e i = -P i we get a sufficient condition for (5.17) (and thus for (5.13)):

2 1 0 -P i (sτ i (ξ) + (1 -s)T i (ξ))(1 -s)ds -a 2 i ≤ 0 for a.e. ξ in - ∆x 2∆t , 0 . (5.24)
Noticing that in the solution T i (ξ) = τ i,L or τ i,R and using the strict convexity of τ → P i (τ ), we get for a.e. ξ in -∆x 2∆t , 0 :

2 1 0 -P i (sτ i (ξ) + (1 -s)T i (ξ))(1 -s)ds ≤ max s∈[0,1] {-P i (sτ i (ξ) + (1 -s)T i (ξ))} 2 1 0 (1 -s)ds ≤ ess sup ξ∈[-∆x 2∆t , ∆x 2∆t ] { -P i (τ i (ξ))} < a 2 i (5.25)
by Whitham's condition. This concludes the proof of inequality (5.13) under Whitham's condition.

We may now prove the following theorem which states that under the CFL condition (5.7), Whitham's condition (5.12) guarantees a discrete entropy inequality for the relaxation finite volume scheme. This is a classical corollary of Lemma 5.3. We refer to [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] for the proof.

Theorem 5.4. Assume the CFL condition (5.7) and suppose that for all n ∈ N and j ∈ Z, the pair (a 1 , a 2 ) j+ 1 2 satisfies Whitham's condition for (U n j , U n j+1 ). Then the relaxation scheme satisfies the following discrete entropy inequality:

η(U n+1 j ) -η(U n j ) + ∆t ∆x (H(U n j , U n j+1 ) -H(U n j-1 , U n j )) ≤ 0, (5.26)
where the numerical entropy flux is given by H

(U L , U R ) = F r η (W r (0 + ; M (U L ), M (U R ))
). This can be seen as a stability condition because if one considers the discrete L 1 -norm of the total mixture energy at time t n : j∈Z η(U n j )∆x , then summing inequality (5.26) over the cells yields

j∈Z η(U n+1 j )∆x ≤ j∈Z η(U n j )∆x, for all n in N, (5.27)
which means that the total mixture energy is decreasing in time.

5.5 Practical choice of the pair (a 1 , a 2 )

The pair of parameters (a 1 , a 2 ), which is computed at each interface x j+ 1 2 must be chosen large enough so as to satisfy several requirements:

• In order to ensure the stability of the relaxation approximation, a i must satisfy Whitham's condition (5.12). For simplicity however, we do not impose Whitham's condition everywhere in the solution of the Riemann problem (2.1)-(2.11), but only for the left and right initial data at each interface: for i in {1, 2}, a 2 i > max(-P i (τ i,L ), -P i (τ i,R )).

(5.28)

In practice, no instabilities were observed during the numerical simulations due to this simpler Whitham-like condition.

• The specific volumes τ i,L and τ i,R must be positive. As the initial conditions of the local Riemann problems are such that for which the exact solution is composed of a {u 1 -c 1 }-shock wave, followed by a {u 2 -c 2 }-rarefaction wave, followed by a u 2 -contact discontinuity, followed by a {u 2 + c 2 }-shock and finally followed by a {u 1 + c 1 }-rarefaction wave (see Figure 6.1). The intermediate states are given by: At each interface x j+ 1 2 , Newton's method is used in order to compute the solution M * L of (4.75). The iterative procedure is stopped when the error is less than 10 -12 .

(W L , W R ) = (M (U L ), M (U R )), we have T i,L = τ i,L and T i,R = τ i,R . Thus τ i,L = τ i,L + 1 2a i (u i,R -u i,L ) - 1 2a 2 i (P i (τ i,R ) -P i (τ i,L )) , (5.29)
τ i,R = τ i,R + 1 2a i (u i,R -u i,L ) + 1 2a 2 i (P i (τ i,R ) -P i (τ i,L )) . ( 5
U 1 = (0.
In Figure 6.1, the approximate solution computed with the relaxation scheme is compared with both the exact solution and the approximate solution obtained with 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 results show that unlike Rusanov's scheme, the relaxation method correctly captures 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 method than by Rusanov's scheme for which the numerical diffusion is larger. We can also see that for the phase 1 variables, there are no oscillations as one can see for Rusanov's scheme: the curves are monotone between the intermediate states. As for phase 2, the intermediate states are captured by the relaxation method 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 is 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 , normalized by the discrete L 1 -norm of the exact solution:

E(∆x) = cellsj |ϕ n j -ϕ ex (x j , T )|∆x cellsj |ϕ ex (x j , T )|∆x , (6.2)
where ϕ is any of the conservative variables (α 1 , α 1 ρ 1 , α 1 ρ 1 u 1 , α 2 ρ 2 , α 2 ρ 2 u 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 6.1, the error E(∆x) at the final time T = 0.14, is plotted against ∆x in a log -log scale. Only the error on the phase fraction α 1 converges towards zero with the expected order of ∆x 1/2 , while the other variables seem 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 much more refined meshes in order to reach the expected order of ∆x 1/2 . Figure 6.1 also displays the error on the conservative variables with respect to the CPU-time of the calculation expressed in seconds. Each point of the plot corresponds to one single calculation for a given mesh size (going from 400 to 102400 cells for the relaxation scheme and from 800 to 102400 cells for Rusanov's scheme). One can see that, for all the variables except α 1 ρ 1 u 1 , 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 instance, for the same error on the phase fraction α 1 , the gain in computational cost is more than 13 when using the relaxation method rather than Rusanov's scheme. For the variable α 1 ρ 1 u 1 , even if the two methods seem to provide similar results, the relaxation method seems to give slightly better results for mesh sizes beyond 10000 cells.

x t u 2 u 2 + c 2 u 1 -c 1 U L U R U 1 U 2 U 3 U 4 u 1 + c 1 u 2 -c 2
Wave structure of the exact Riemann solution

Wave structure of the exact Riemann solution

Test-case 2: a vanishing phase case

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. Assuming for instance that α 1,L = 1 and 0 < α 1,R < 1, the right state is a mixture of both phases while the left initial state is composed solely of phase 1. This type of vanishing-phase Riemann solution is introduced in [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] in the more general context of the complete Baer-Nunziato system with energy equations. The construction of an exact solution with vanishing phase fraction is not classical and we choose to follow the natural approach given in [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF].

α 1 α 1 ρ 1 α 1 ρ 1 u 1 α 2 ρ 2 α 2 ρ 2 u 2 sqrt(∆x) 10
For more details, see appendix A. The solution is composed of a {u 1 -c 1 }-shock wave from U L to U 1 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 righthand side (RHS) region where the two phases are mixed. In this RHS region, the solution is composed of a {u 2 + c 2 }-rarefaction wave connecting U 2 to U 3 followed by a {u 1 + c 1 }-rarefaction wave from U 3 to U R (see Figure 6.1). In practice, the numerical method requires values of α 1,L and α 1,R that lie strictly in the interval (0, 1). Therefore, in the numerical implementation, we take α 1,L = 1 -10 -9 . 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 and u 2,L to the values on the right of the u 2 -contact discontinuity. For the relaxation scheme, this choice enables to avoid oscillations of phase 2 quantities in the region where phase 2 in not present. However, some tests have been conducted that assess that taking other values of (ρ 2,L , u 2,L ) has little impact on the phase 1 quantities as well as on the phase 2 quantities where this phase is present.

At each interface x j+ 1 2 , a dichotomy algorithm is used in order to compute the solution M * L of (4.75). Indeed for such a vanishing phase test-case, Newton's method fails to converge. The dichotomy algorithm is stopped when the error is less than 10 -12 . As for the first test-case, we can see that for the same level of refinement, the relaxation method is 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 relaxation method is much more stable than Rusanov's scheme. Indeed, the relaxation scheme behaves better than Rusanov's scheme when it comes to divisions by small values of α 2 , since the solution approximated by Rusanov's scheme develops important oscillations.

Test-case 3: Coupling between two pure phases

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). The intermediate states are given by U L = (1.0, 0.861773876012754, 3.552800564555003, 4.641588833612778, 1.0), U 1 = (1.0, 2.154434690031884, 1., 4.641588833612778, 1.), U 2 = (0., 2.154434690031884, 1., 4.641588833612778, 1.), U R = (0., 2.154434690031884, 1.0, 6.962383250419167, 1.767119653712349).

(6.3)

The exact solution is composed of a {u 1 -c 1 }-shock wave from U L to U 1 in the LHS region where only phase 1 is present. U 1 is connected to U 2 by a u 2 -contact discontinuity separating the two pure phase regions. In the RHS region, where only phase 2 exists, U 2 is connected to U R by a {u 2 + c 2 }-rarefaction wave. For more details on how the exact solution is constructed see appendix A.2.

In the numerical implementation, we set α 1,L = 1 -10 -9 and α 1,R = 10 -9 . A dichotomy procedure is used in order to compute the solution M * L of (4.75) at each interface x j+ 1 2 . The dichotomy algorithm is stopped when the error is less than 10 -12 .

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, with Rusanov's scheme, strong oscillations pop up in the regions where a phase vanishes. Observe also that unlike Rusanov's scheme, the relaxation scheme does not fail to correctly approximate the evolution of the phase fraction α 1 .

Conclusion and further works

The work performed in this paper provides some interesting teachings. First of all, the relaxation approximation introduced in (3.1) seems to be a reasonable choice of continuous approximation for the isentropic Baer-Nunziato model. Indeed, thanks to the linear degeneracy of all the fields, one has been able to exactly solve the associated relaxation Riemann problem in the framework of solutions with subsonic wave ordering. More surprising is the fact that we have been able to provide explicit necessary and sufficient conditions on the initial left and right data for the existence of such solutions with subsonic wave ordering.

One natural question is the possibility of extending this result to all the wave-orderings that may be considered for the solution, namely, supersonic and sonic wave orderings. The eventual aim would be to construct a partition of the whole space of initial conditions Ω r × Ω r , where each element of the partition would be characterized by explicit conditions similar to (A) that indicate the waveordering of the solution for a given pair (W L , W R ). Unfortunately, performing the study of all the wave-orderings is not a trivial task because it may involve δ-shock type solutions due to the linear resonance phenomena (see [START_REF] Bouchut | Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness[END_REF]) when a (linearized) acoustic field interacts with the transport wave u 2 . As a matter of fact, Dirac solutions have already been observed for the relaxation system (4.33) in the context of Euler's equations in nozzle (see [START_REF] Coquel | A Robust and Entropy-Satisfying Numerical Scheme for Fluid Flows in Discontinuous Nozzles[END_REF]). For this simpler system, which is a major building block in the resolution of the relaxation Riemann problem (4.1)-(4.2) for the Baer-Nunziato model, the Riemann problem has been completely solved for any initial data, and δ-shock solutions have been introduced for sonic flows.

As for Euler's equations in a nozzle, the relaxation approximation introduced here provides a simple Riemann solver of Harten, Lax and van Leer, which enables a stable treatment of vanishing phase cases. Indeed, thanks to the introduction of a generalized class of solutions for which energy dissipation through the linearly degenerate void fraction wave has been allowed, Riemann problems with arbitrarily small initial phase fractions may be computed by the solver. Introducing energy dissipation in the first step of the relaxation approximation is an original idea since all the fields of the relaxation system are linearly degenerate. Another more natural choice would be to increase the relaxation parameters (a 1 , a 2 ) thus introducing some dissipation in the method but this, contrary to our approach, would restrict the time-step to potentially small values in the regimes of vanishing phases. In addition, the relaxation Riemann solver is proved to preserve positive densities and to satisfy a discrete entropy inequality which guarantees the stability of the derived scheme. To our knowledge, there exists no other scheme that is proved to satisfy these two properties. The numerical tests prove the good behavior of the scheme in vanishing phase cases. Moreover, for a given error on the solution, it is shown that the relaxation scheme compares much favorably with Rusanov's scheme in terms of computational cost.

An important sequel of this work is the extension of the method to the full Baer-Nunziato model with energies thanks to the energy/entropy duality. Indeed, thanks to the second principle of thermodynamics which connects the phasic energies and the transported phasic entropies, one is able to extend this Riemann solver to the full model with energies through minor adaptations (see [START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF] for a related framework). This is explained in [START_REF] Saleh | Analyse et Simulation Numérique par Relaxation d'Ecoulements Diphasiques Compressibles. Contribution au Traitement des Phases Evanescentes[END_REF] and is the purpose of the forthcoming article [START_REF] Coquel | An entropy-eatisfying relaxation scheme for the full 3D Baer-Nunziato model with energies[END_REF], together with the multi-dimension extension of the method.

Definition 4 . 1 .

 41 Let (W L , W R ) be two states in Ω r . A solution to the Riemann problem (4.1)-(4.2)

Lemma 4 . 1 .

 41 Consider the Riemann problem (4.1)-(4.2) with α

Remark 4 . 2 .

 42 0 is the constant mass flux across the standing wave andϕ(M; ν, ω) = ω 2 (νM + 1)(M + 1) -(νM -1)(M -1). (4.46) As M lies in the interval (0, min(1, 1/ν)), one has (1 + νM)(1 -M) > 0.Hence, the sign of the energy-jump across the standing wave is determined by the sign of ϕ (M; ν, ω). If M is such that ϕ (M; ν, ω) = 0, then the solution is energy-preserving. If M is such that ϕ (M; ν, ω) < 0, the solution is energy-dissipative. The value M 0 (ν, ω) given in (4.42) is the solution of ϕ (M; ν, ω) = 0.

  If ν > 1, for µ close enough to one, one has M µ (ν, M * L ) close to M * L ν and then M(ν, M * L ) ∈ (0, M 0 (ν, ω)] which implies that the energy is now dissipated since ϕ (M; ν, ω) ≤ 0 for all M ∈ (0, M 0 (ν, ω)]. Formula (4.49) is referred to as the kinetic relation since the choice of M prescribes an energy-dissipation rate through (4.45) and determines a unique solution to the Riemann-problem (4.33)-(4.34) (or equivalently to the Riemann problem (4.20)-(4.21)).

  .52) easy manipulations using the relation w 1 = u 1 -u * 2 show that the jump relations (4.43)-(4.44)-(4.45) are equivalent to the following jump relations in the original frame:

4. 6

 6 Solution to the fixed-point problem and proof of Theorem 4.2

α 1,L α 1 ,

 1 R is large, that the solution M * L of the fixed-point problem (4.64) is such that the specific volume τ 1,R * (see (4.41) with

Figure 1 :

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

Figure 2 :

 2 Figure 2: Test-case 1: L 1 -Error with respect to ∆x and L 1 -Error with respect to computational cost (in seconds), for the conservative variables.

3 u 2 + c 2 u 1 + c 1

 32211 Wave structure of the exact Riemann solutionWave structure of the exact Riemann solution

Figure 3 :

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

Figure 4 :

 4 Figure 4: Test-case 3: Structure of the solution and space variations of the physical variables at the final time T = 0.07. Mesh size: 1000 cells.

  4.1. The expression of u * 2 given in equation (4.28) defines the function F [W L , W R ; a 2 ] introduced in (4.19), since u * 2 is expressed as a function of π * 1 . It clearly appears that if α 1,L = α 1,R , the non-conservative product vanishes and the resolution of the Riemann problem yields u * 2 = u 2 as seen in Lemma 4.1.

  .55) Finally, observe that expression (4.51) defines the function G [W L , W R ; a 1 ] introduced in (4.23) since it expresses π * 1 as a function of u * 2 (through the relative Mach number M * L ).

  while π * 1 is obtained through (4.58). Then the intermediate states for phase 2 are given by equations (4.28) and (4.29) in Proposition 4.4. Once prescribed the value M := M(ν, M * L ) according to (4.49), the intermediate states for phase 1 are given in equations (4.39) to (4.41) of Proposition 4.6, except for the velocities to which one must add u * 2

  1. Consider the intermediate states

	U L = (1.0, 1.8, 0.747051068928543, 3.979765198025580, 0.6),
	U 1 = (1.0, 2.0, 0.4, 3.979765198025580, 0.6),
	U 2 = (0.4, 1.982040094756841, 0.095469338564172, 3.979765198025580, 0.6),
	U 3 = (0.4, 1.9820400948, 0.0954693386, 5.1736947574, 1.0690676047),
	U R = (0.4, 2.081142099494683, 0.267119045902047, 5.173694757433254, 1.069067604724276).

u * 1 u 1 -a 1 τ 1 u 2 -a 2 τ 2 u 1 -a 1 τ 1 α 1,R α 1,L W -W +

Acknowledgements. The authors would like to thank Jean-Paul Daniel for his helpful remarks. This work has been partially funded by ANRT and EDF through an EDF-CIFRE contract 529/2009. The forth author is partially supported by the LRC Manon (Modélisation et Approximation Numérique Orientées pour l'énergie Nucléaire -CEA DM2S/LJLL).

Equations (5.29) and (5.30) are two second order polynomials in a -1 i , and by taking a i large enough, one can guarantee that τ i,L > 0 and τ i,R > 0, since the initial specific volumes τ i,L and τ i,R are positive.

• Finally, (a 1 , a 2 ) must be chosen so as to meet the condition (A) of Theorem 4.2 as well as the positivity condition of the phase 2 densities (B).

Thereafter, we propose an iterative algorithm for the computation of the parameters (a 1 , a 2 ) at each interface. Fixedpoint(a 1 , a 2 ) is a subroutine that computes a numerical approximation of the solution u * 2 of the fixed-point problem (4.24), using some numerical method such as Newton's method or a dichotomy algorithm. The notation not(P) is the negation of the logical statement P.

Algorithm for the choice of (a 1 , a 2 ):

Choose κ a (small) parameter in the interval (0, 1). For all interface x j+ 1 2 , j ∈ Z, calculate a 1,j+ 1 2 and a 2,j+ 1 2 as follows • For i in {1, 2} initialize a i,j+ 1 2 :

It is possible to prove that this algorithm always converges in the sense that there is non infinite looping due to the while-conditions. 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 (5.7) 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 6).

Numerical tests

In this section, we present Riemann-type test-cases on which the performances of the relaxation scheme are tested. The phasic equations of state are given by the following ideal gas pressure laws:

2 , with κ 2 = 1 and γ 2 = 1.5.

(6.1)

In the sequel, U = (α 1 , ρ 1 , u 1 , ρ 2 , u 2 ) T denotes the state vector in non-conservative variables.

A Appendices

A.1 Vanishing-phase solution of test-case 2

The vanishing phase solution of test-case 2 is constructed as follows. First of all, two states U -and U + are constructed so as to be connected by a u 2 -contact discontinuity with α - 1 = α 1,L = 1 and α + 1 = α 1,R = 0.4 (U -resp. U + is denoted U 1 resp. U 2 in (6.2)). For this purpose, the jump relations (2.15)-(2.16)-(2.17)-(2.18) associated with the contact discontinuity reduce to

since α - 1 = 1 and α - 2 = 0. Thanks to these jump relations, given the values ρ - 1 = 2.0, u - 1 = 0.4 and u - 2 = 0.6 (no value of ρ - 2 is needed), we compute the values of ρ + 1 , u + 1 , ρ + 2 and u + 2 which are given in the intermediate state U 2 in (6.2). The value of ρ - 2 is then imposed to be equal to ρ + 2 . Then the state U -= U 1 is connected on its left with the state U L through a {u 1 -c 1 }-shock. The state U 2 is connected to U 3 through a {u 2 + c 2 }-rarefaction wave, and finally U 3 is connected to U R through a {u 1 + c 1 }-rarefaction wave.

A.2 Coupling solution of test-case 3

The same procedure is implemented for the construction of the exact solution of test-case 3, which corresponds to a coupling between a pure phase 1 (α - 1 = α 1,L = 1) on the right and a pure phase 2 (α + 1 = α 1,R = 0) on the left. This time, the jump relations of the u 2 -contact discontinuity reduce to

A solution is given by u

p) where p = p - 1 = p + 2 = 10 is the common pressure. The values of ρ + 1 and u + 1 , which are of no importance since phase 1 in not present on the right of the contact discontinuity, are taken equal to ρ - 1 and u - 1 . Similarly, we set ρ - 2 := ρ + 2 . This concludes the construction of the intermediate states U -= U 1 and U + = U 2 . Then the state U 1 is connected on its left with the state U L through a {u 1 -c 1 }-shock and the state U 2 is connected to U R through a {u 2 + c 2 }-rarefaction wave.