A relaxation scheme for a hyperbolic multiphase flow model. Part I: barotropic eos

This article is the first of two in which we develop a relaxation finite volume scheme for the convective part of the multiphase flow models introduced in the series of papers [10, 9, 4]. In the present article we focus on barotropic flows where in each phase the pressure is a given function of the density. The case of general equations of state will be the purpose of the second article. We show how it is possible to extend the relaxation scheme designed in [8] for the barotropic Baer-Nunziato two-phase flow model to the multiphase flow model with N-where N is arbitrarily large-phases. The obtained scheme inherits the main properties of the relaxation scheme designed for the Baer-Nunziato two phase flow model. The approximated phase fractions and phase densities are proven to remain positive and a discrete energy inequality is also proven under a classical CFL condition. For the same level of refinement, the relaxation scheme is shown to be much more accurate than Rusanov's scheme, and for a given level of approximation error, the relaxation scheme is shown to perform much better in terms of computational cost than Rusanov's scheme. Moreover, contrary to Rusanov's scheme which develops strong oscillations when approximating vanishing phase solutions, the numerical results show that the relaxation scheme remains stable in such regimes.

cope with arbitrarily small values of the statistical phase fractions, which are proven to remain positive as well as the phase densities. Finally, a fully discrete energy inequality is also proven under a classical CFL condition.
In [11], the relaxation scheme has shown to compare well with two of the most popular existing schemes available for the full Baer-Nunziato model (with energy equations), namely Schwendeman-Wahle-Kapila's Godunov-type scheme [24] and Tokareva-Toro's HLLC scheme [27]. Still for the Baer-Nunziato model, the relaxation scheme also has shown a higher precision and a lower computational cost (for comparable accuracy) than Rusanov's scheme. Regarding the multiphase flow model considered in the present paper, the relaxation scheme is compared to Rusanov's scheme, the only scheme presently available. As expected, for the same level of refinement, the relaxation scheme is shown to be much more accurate than Rusanov's scheme, and for a given level of approximation error, the relaxation scheme is shown to perform much better in terms of computational cost than Rusanov's scheme. Moreover, the numerical results show that the relaxation scheme is much more stable than Rusanov's scheme which develops strong oscillations in vanishing phase regimes.
The relaxation scheme described here is restricted to the simulations of flows with subsonic relative speeds, i.e. flows for which the difference between the material velocities of the phases is less than the monophasic speeds of sound. For multiphase flow simulations in the nuclear industry context, this is not a restriction, but it would be interesting though to extend the present scheme to sonic and supersonic flows.
For the sake of concision and simplicity, this work is only concerned with barotropic equations of states. However, as it was done for the two phase Baer-Nunziato model in [11], an extension of the relaxation scheme to the full multiphase flow model with energy equations is within easy reach. This will be the purpose of a companion paper.
The paper is organized as follows. Section 2 is devoted to the presentation of the multiphase flow model. In Section 3 we explain how to extend the relaxation Riemann solver designed in [12] to the multiphase flow model and Section 4 is devoted to the numerical applications on the three phase flow model. In addition to a convergence and CPU cost study in Test-case 1, we simulate in Test-cases 2 and 3 vanishing phase configurations where two of the three phases have nearly disappeared in some space region. In particular, Test-case 3 is dedicated to the interaction of a gas shock wave with a lid of rigid particles.

The multiphase flow model
We consider the following system of partial differential equations (PDEs) introduced in [17] for the modeling of the evolution of N distinct compressible phases in a one dimensional space: for k = 1, . . . , N , x ∈ R and t > 0: l =k P kl (U)∂ x α l = 0. (2.1c) The model consists in N coupled Euler-type systems. The quantities α k , ρ k and u k represent the mean statistical fraction, the mean density and the mean velocity in phase k (for k = 1, . . . , N ). The quantity p k is the pressure in phase k. We assume barotropic pressure laws for each phase so that the pressure is a given function of the density ρ k → p k (ρ k ) with the classical assumption that p k (ρ k ) > 0. The mean statistical fractions and the mean densities are positive and the following saturation constraint holds everywhere at every time: N k=1 α k = 1. Thus, among the N equations (2.1a), N − 1 are independent and the main unknown U is expected to belong to the physical space: Ω U = U = (α 1 , . . . , α N −1 , α 1 ρ 1 , . . . , α N ρ N , α 1 ρ 1 u 1 , . . . , α N ρ N u N ) T ∈ R 3N −1 , such that 0 < α 1 , . . . , α N −1 < 1 and α k ρ k > 0 for all k = 1, . . . , N .
Following [17], several closure laws can be given for the so-called interface velocity V I (U) and interface pressures P kl (U). Throughout the whole paper, we make the following choice: V I (U) = u 1 , and for k = 1, P 1l (U) = p l (ρ l ), for l = 2, . . . , N for k = 1, P kl (U) = p k (ρ k ), for l = 1, . . . , N, l = k. (2.2) With this particular choice, observing that the saturation constraint gives N l=1,l =k ∂ x α l = −∂ x α k for all k = 1, . . . , N the momentum equations (2.1c) can be simplified as follows: Remark 2.1. When N = 2, system (2.1) is the convective part of the Baer-Nunziato two phase flow model [1]. This model is thus an extension of the Baer-Nunziato two phase flow model to N (possibly ≥ 3) phases. As for the Baer-Nunziato model, in the areas where all the statistical fractions α k are constant in space, system (2.1) consists in N independent Euler systems weighted by the statistical fraction of the corresponding phase. These Euler systems are coupled through non-conservative terms which are active only in the areas where the statistical fractions gradients are non zero.
Remark 2.2. The choice V I (U) = u 1 is classical for the two phase flow model when phase 1 is dispersed and phase 2 prevails in the fluid. It is then natural to take an interfacial velocity which is equal to the material velocity of the dispersed phase. For three phase flows, the choice V I (U) = u 1 has been made in [6]. When simulating the steam explosion phenomenon [3,21], it corresponds to taking an interfacial velocity equal to the material velocity of the corium particles.
Proof. The proof can be found in [18].
Remark 2.4. The system is not hyperbolic in the usual sense because when (2.5) is not satisfied, the right eigenvectors do not span the whole space R 3N −1 . Two possible phenomena may cause a loss of the strict hyperbolicity: an interaction between the linearly degenerate field of velocity u 1 with one of the acoustic fields of the phase k for k = 2, . . . , N , and vanishing values of one of the phase fractions α k , k = 1, . . . , N . In the physical configurations of interest in the present work (such as three phase flows in nuclear reactors), the flows have strongly subsonic relative velocities, i.e. a relative Mach number much smaller than one: so that resonant configurations corresponding to wave interaction between acoustic fields and the u 1 -contact discontinuity are unlikely to occur. In addition, following the definition of the admissible physical space Ω U , one never has α k = 0. However, α k = 0 is to be understood in the sense α k → 0 since one aim of this work is to construct a robust enough numerical scheme that could handle all the possible values of α k ∈ (0, 1), k = 1, . . . , N , especially, arbitrarily small values.
An important consequence of the closure law V I (U) = u 1 is the linear degeneracy of the field associated with the eigenvalue σ 1 (U) = . . . = σ N −1 (U) = u 1 . This allows to define solutions with discontinuous phase fractions through the Riemann invariants of this linear field. Indeed, as proven in [18], there are 2N independent Riemann invariants associated with this field which is enough to parametrize the integral curves of the field since the multiplicity of the eigenvalue is N − 1. This can be done as long as the system is hyperbolic i.e. as long as (2.5) is satisfied, which prevents the interaction between shock waves and the non conservative products in the model.
An important consequence of the closure law (2.2) for the interface pressures P kl (U) is the existence of an additional conservation law for the smooth solutions of (2.1). Defining the specific internal energy of phase k, e k by e k (ρ k ) = p k (ρ k )/ρ 2 k and the specific total energy of phase k by E k = u 2 k /2 + e k (ρ k ), the smooth solutions of (2.1) satisfy the following identities: In [18], the mappings U → (α k ρ k E k )(U) are proven to be (non strictly) convex for all k = 1, . . . , N . Since, as long as the system is hyperbolic, the gradients of α k are supported away from shock waves, it is natural to use theses mappings as mathematical entropies of system (2.1) and select the physical non-smooth weak solutions of (2.1) as those which satisfy the following entropy inequalities in the weak sense: If a shock appears in phase 1, inequality (2.9) is strict and if a shock appears in phase k for some k ∈ {2, . . . , N } inequality (2.10) is strict. Summing for k = 1, . . . , N , the entropy weak solutions of (2.1) are seen to satisfy the following total energy inequality: Obviously, for smooth solutions, (2.11) is an equality.

A relaxation approximate Riemann solver
As for the Baer-Nunziato two phase flow model, system (2.1) has genuinely non-linear fields associated with the phasic acoustic waves, which makes the construction of an exact Riemann solver very difficult. Following similar steps as in [12], we introduce a relaxation approximation of the multiphase flow model (2.1) which is an enlarged system involving N additional unknowns T k , associated with linearizations of the phasic pressure laws. These linearizations are designed to get a quasilinear enlarged system, shifting the initial non-linearity from the convective part to a stiff relaxation source term. The relaxation approximation is based on the idea that the solutions of the original system are formally recovered as the limit of the solutions of the proposed enlarged system, in the regime of a vanishing relaxation coefficient ε > 0. For a general framework on relaxation schemes we refer to [5,9,10].
We propose to approximate the solutions of (2.1) by the solutions of the following Suliciu relaxation type model (see [5,25,26]) in the limit ε → 0: where the relaxation unknown is now expected to belong to the following enlarged phase space: such that 0 < α 1 , . . . , α N −1 < 1, α k ρ k > 0 and α k ρ k T k > 0 for k = 1, . . . , N , and where: The saturation constraint is still valid: For each phase k = 1, . . . , N , τ k = ρ −1 k is the specific volume of phase k and the pressure π k is a (partially) linearized pressure π k (τ k , T k ), the equation of state of which is defined by: where τ → P k (τ ) := p k (τ −1 ) is the pressure of phase k seen as a function of the specific volume τ = ρ −1 . Accordingly with (2.2) the relaxation interface pressure Π kl (W) is defined by: When N = 2, system (3.1) is exactly the same relaxation approximation introduced in [12] for the Baer-Nunziato two phase flow model. In the formal limit ε → 0, the additional variable T k tends towards the specific volume τ k , and the linearized pressure law π k (τ k , T k ) tends towards the original non-linear pressure law p k (ρ k ), thus recovering system (2.1) in the first 3N − 1 equations of (3.1). The constants (a k ) k=1,...,N in (3.3) are positive parameters that must be taken large enough so as to satisfy the following sub-characteristic condition (also called Whitham's condition): for all the values T k encountered in the solution of (3.1). Performing a Chapman-Enskog expansion, one can see that Whitham's condition expresses that system (3.1) is a viscous perturbation of system (2.1) in the regime of small ε. At the numerical level, a fractional step method is commonly used in the implementation of relaxation methods: the first step is a time-advancing step using the solution of the Riemann problem for the convective part of (3.1): (3.6) while the second step consists in an instantaneous relaxation towards the equilibrium system by imposing T k = τ k in the solution obtained by the first step. This second step is equivalent to sending ε to 0 instantaneously. As a consequence, we now focus on constructing an exact Riemann solver for the homogeneous convective system (3.6). Let us first state the main mathematical properties of system (3.6). The linearization (3.3) is designed so that system (3.6) has only linearly degenerate fields, thus making the resolution of the Riemann problem for (3.6) easier than for the original system (2.1). We have the following proposition: is weakly hyperbolic on Ω W in the following sense. It admits the following 4N − 1 real eigenvalues: All the characteristic fields are linearly degenerate and the corresponding right eigenvectors are linearly independent if, and only if, Proof. The proof is given in Appendix A.
Remark 3.2. Here again, one never has α k = 0 for W ∈ Ω W . However, α k = 0 is to be understood in the sense α k → 0.
We also have the following properties: The smooth solutions as well as the entropy weak solutions of (3.6) satisfy the following phasic energy equations: Summing for k = 1, . . . , N , the smooth solutions and the entropy weak solutions of (3.6) are seen to satisfy an additional conservation law: Under Whitham's condition (3.5), to be met for all the T k under consideration, the following Gibbs principles are satisfied for k = 1, . . . , N : (3.11) where E k (u k , τ k ) = u 2 k /2 + e k (ρ k ).
Proof. The proof of (3.8) and (3.9) follows from classical manipulations. From the phasic mass and momentum equations we first derive the evolution equation satisfied by u 2 k /2. We then derive an equation for π k (τ k , T k ), using the mass equation of phase k and the advection equation of T k . Combining these two equations and the fact that T k is advected, we obtain (3.8) and (3.9). The proof of Gibbs principle follows from an easy study of the function T k → E k (u k , τ k , T k ).
Remark 3.4. Since all the characteristic fields of system (3.6) are linearly degenerate, the mixture energy equation (3.10) is expected to be satisfied for not only smooth but also weak solutions. However, as we will see later when constructing the solutions of the Riemann problem, in the stiff cases of vanishing phases where one of the left or right phase fractions α k,L or α k,R is close to zero for some k ∈ {1, . . . , N }, ensuring positive values of the densities of phase k requires an extra dissipation of the mixture energy by the computed solution.

The relaxation Riemann problem
Let (W L , W R ) be two constant states in Ω W and consider the Riemann problem for system (3.6) with the following initial condition: 3.1.1. Definition of the solutions to the Riemann problem Following Proposition 3.1, a solution to the Riemann problem (3.6)-(3.12) is expected to be a self-similar function composed of constant intermediate states separated by waves which are contact discontinuities associated with the system's eigenvalues. Since the phase fractions are transported by the material velocity of phase 1, the non-conservative products involving the phase fraction gradients are only active across this wave and the phases are independent away from this wave. In particular, for a fixed k in {2, . . . , N }, the phase k quantities may change across the contact discontinuities associated with the eigenvalues u k − a k τ k , u k , u k − a k τ k and u 1 and are constant across the other waves. For the applications aimed at by this work, we are only interested in solutions which have a subsonic wave ordering: Consequently, in the self-similar Riemann solution, the propagation velocity u 1 lies in-between the acoustic waves of all the other phases. Moreover, ensuring the positivity of the phase 1 densities also requires the material wave u 1 to lie in between the acoustic waves of phase 1.
We now introduce the definition of solutions to the Riemann problem (3.6)-(3.12), which is a straightforward extension of the definition of solutions in the case of two phase flows (see [12], Def. 4.1).
Definition 3.5. Let (W L , W R ) be two states in Ω W . A solution to the Riemann problem (3.6)-(3.12) with subsonic wave ordering is a self-similar mapping W(x, t) = W Riem (x/t; W L , W R ) where the function ξ → W Riem (ξ; W L , W R ) satisfies the following properties: (i) W Riem (ξ; W L , W R ) is a piecewise constant function, composed of 3N waves (associated with the eigenvalues u k − a k τ k , u k and u k + a k τ k for k = 1, . . . , N ) separating 3N + 1 constant intermediate states belonging to Ω W , and such that: satisfies the following system of PDEs in the distributional sense: also satisfies the following energy equations in the distributional sense: The solution has a subsonic wave ordering in the following sense: Remark 3.6. Following (3.4), there are N − 1 interface pressures corresponding to the phase pressures (π 2 , . . . , π N ). Moreover, the saturation constraint (3.2) gives N l=1,l =k ∂ x α l = −∂ x α k for all k = 1, . . . , N , which justifies the simplified form of the non-conservative product D * (W L , W R )δ 0 (x − u * 1 t). Remark 3.7. Equation (3.16) is a relaxed version of (3.9) in which the energy of phase k is allowed to be dissipated across the u 1 -wave despite the linear degeneracy of the associated field. As explained in [12] for the relaxation approximation of the Baer-Nunziato two phase flow model, allowing such a dissipation may be necessary when an initial phase fraction α k,L or α k,R is close to zero, in order to ensure the positivity of all the intermediate densities.

The resolution strategy: a fixed-point research
Following the ideas developed in [12], the resolution of the Riemann problem (3.6)-(3.12) is based on a fixed-point research which formally amounts to iterating on a two step procedure involving the pair (u * 1 , Π * ) ∈ R × R N −1 . We first remark that system (3.14) can be written in the following form: and for k = 2, . . . , N : First step: The family of interface pressures Π * = (π * 2 , . . . , π * N ) ∈ R N −1 defining the non-conservative products . . , N are assumed to be known. Hence, system (S 1 ), which gathers the governing equations for phase 1, is completely independent of the other phases since the non-conservative terms can now be seen as a known source term and the Riemann problem for (S 1 ) can therefore be solved independently of the other phases. Observe that system (S 1 ) is very similar to System (4.16) of [12] in the two phase flow framework. This Riemann problem is easily solved since (S 1 ) is a hyperbolic system with a source term which is a Dirac mass supported by the kinematic wave of velocity u * 1 . Hence, there is no additional wave due to the source term.
Consequently, knowing a prediction of the interface pressures Π * = (π * 2 , . . . , π * N ) ∈ R N −1 , one can explicitly compute the value of the kinematic speed u * 1 by solving the Riemann problem associated with phase 1. This first step enables to define a function: (3.18) Second step: The advection velocity u * 1 of the phase fractions α k is assumed to be a known constant. Thus, for all k = 2, . . . , N , the governing equations for phase k, gathered in system (S k ), form a system which is independent of all the other systems (S l ) with l = 1, . . . , N , l = k. In addition to the kinematic velocity u k and the acoustic speeds u k ± a k τ k , the Riemann problem for (S k ) involves an additional wave whose known constant velocity is u * 1 . This wave is weighted with an unknown weight π * k ∆α k (only for the momentum equation) which is calculated by solving the Riemann problem for (S k ) and then applying Rankine-Hugoniot's jump relation to the momentum equation for the traveling wave of speed u * 1 . Here again, we observe that system (S k ) is the exact same system already solved in the two phase flow framework (see [12], Syst. (4.20)). This is what justifies the straightforward extension of the relaxation scheme to the multiphase flow model. Solving all these Riemann problems for (S k ) with k = 2, . . . , N , this second step allows to define a function: Fixed-point research: Performing an iterative procedure on these two steps actually boils down to the following fixed-point research.
The interval in which u * 1 must be sought corresponds to the subsonic condition (3.17) on the one hand, and to the positivity of the intermediate states of phase 1 on the other hand.
Let us now introduce some notations which depend explicitly and solely on the initial states (W L , W R ) and on the relaxation parameters (a k ) k=1,...,N . For k = 1, . . . , N : Remark 3.8. Observe that with these definitions, u k,L − a k τ k,L = u k − a k τ k,L and u k,R + a k τ k,R = u k + a k τ k,R .

First step of the fixed-point procedure: solving phase 1
In this first step, the interface pressures Π * = (π * 2 , . . . , π * N ) ∈ R N −1 defining the non-conservative products . . , N are first assumed to be known and one solves the Riemann problem for (S 1 ) with the initial condition where W 1 = (α 1 , α 1 ρ 1 , α 1 ρ 1 u 1 , α 1 ρ 1 T 1 ) 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. Observe that the non-conservative products π * k ∂ x α k are not ambiguous here since for all k = 2, . . . , N , π * k is a known constant. System (S 1 ) is very similar to system (4.16) of [12] encountered in the two phase flow framework, and the resolution of the corresponding Riemann problem follows from the exact same steps. Therefore, we only state the main results and the reader is referred to Section 4.4 of [12] for the proofs.
The following proposition characterizes the convective behavior of system (S 1 ).
Proposition 3.9. System (S 1 ) is a hyperbolic system with linearly degenerate fields associated with the eigenvalues u 1 − a 1 τ 1 , u 1 and u 1 + a 1 τ 1 . The eigenvalue u 1 has multiplicity 2.
Proof. The proof is similar to that of Proposition 3.1 which is given in Appendix A.
We have the following well-posedness result for the governing equations of phase 1. The Riemann problem (S 1 )-(3.22) differs from the Riemann problem in the two phase case in Proposition 4.4 of [12] only by the value of the source term N l=2 π * l ∂ x α l . Hence, the proof follows from similar steps as the proof of Proposition 4.4 from [12].
The intermediate densities ρ − 1 and ρ + 1 are positive if and only if Moreover, this unique solution satisfies the following energy equation in the usual weak sense: Proof. The proof consists in solving a system with the unknowns W − 1 and W + 1 . All the fields are linearly degenerated. The components of these state vectors are linked to each other through the Riemann invariants of the u * 1 -wave. The quantity N l=2 π * l ∂ x α l is a given source term in the right hand side of the Rankine-Hugoniot relation associated with the momentum equation for this wave. W − 1 and W 1,L are connected through the Riemann invariants of the {u 1 − a 1 τ 1 }-wave. W + 1 and W 1,R are connected through the Riemann invariants of the {u 1 + a 1 τ 1 }-wave. We refer to Proposition 4.4 of [12] for more details.
The expression of u * 1 given in equation (3.24) A convenient reformulation of (3.24) is the following: 3.1.4. Second step of the fixed-point procedure: solving phase k, for k = 2, . . . , N In this second step, the transport velocity u * 1 of the phase fractions α k is assumed to be known, while the vector of interface pressures Π is an unknown that must be calculated by solving the N − 1 independent Riemann problems for (S k ), for k = 2, . . . , N with the initial condition where W k = (α k , α k ρ k , α k ρ k u k , α k ρ k T k ) denotes the state vector for phase k, and (W k,L , W k,R ) are the restriction of the complete initial data (W L , W R ) to the phase k variables. System (S k ) is the exact same system as system (4.20) of [12] encountered in the two phase flow framework, and the resolution of the corresponding Riemann problem follows from the exact same steps. Here again, we only state the main results and the reader is referred to Section 4.5 of [12] for the detailed proofs. Once the resolution is done, applying Rankine-Hugoniot's jump relation to the momentum equation of (S k ) gives the expression of π * k ∆α k .
Proposition 3.11. For every k = 2, . . . , N , system (S k ) admits four real eigenvalues that are u k − a k τ k , u k , u k +a k τ k and u * 1 . All the fields are linearly degenerate and the system is hyperbolic if, and only if |u k −u * 1 | = a k τ k .
Proof. The proof is similar to that of Proposition 3.1 which is given in Appendix A.
We search for Riemann solutions which comply with the subsonic relative speed constraint (3.17). Such solutions are of three types depending on the relative wave ordering between the eigenvalues u k and u * 1 : We denote ξ → W k (ξ; W k,L , W k,R ) the self-similar mapping defining this solution and we may now recall the following result stated in Proposition 4.8 of [12] in the framework of the Baer-Nunziato two phase flow model. Proposition 3.12. Assume that a k is such that τ k,L > 0 and τ k,R > 0 and define for (ν, ω) ∈ R * + × R * + the two-variable function: which can be extended by continuity to ω = 1 by setting M 0 (ν, 1) = 0.
These solutions are parametrized by a real number M and the intermediate states are given by: These solutions satisfy the following energy identity: -The Riemann problem (S k )-(3.29) admits solutions with the subsonic wave ordering u k − a k τ k < u k < u * 1 < u k + a k τ k , if and only if u k − u * 1 ≤ 0 and u k − a k τ k,L < u * 1 . (3.34) By the Galilean invariance of system (S k ) written in the mowing frame of speed u * 1 (see [12], Syst. (4.33)), such a solution is given by ξ → VW k (2u * 1 − ξ; VW k,R , VW k,L ) where the operator V changes the relative velocities u k − u * 1 into their opposite values The solution also satisfies an energy identity similar to (3.33).
The intermediate states are obtained by passing to the limit as M * k → 0 in the expressions given in the case of the wave ordering u * 1 < u k . The solution also satisfies an energy identity similar to (3.33).
Proof. The proof consists in solving a system with the unknowns W − k and W + k and W k,R * (or W k,L * ). All the fields are linearly degenerated. The components of W − k and W + k are linked to each other through the Riemann invariants of the u * 1 -wave. W − k and W k,L are connected through the Riemann invariants of the {u k − a k τ k }wave. W + k and W k,R * are connected through the Riemann invariants of the u k -wave. W k,R * and W k,R are connected through the Riemann invariants of the {u k + a k τ k }-wave. We refer to Proposition 4.8 of [12] for more details. Remark 3.13. As mentioned in Remark 3.7, taking Q k (u * 1 , W L , W R ) > 0 may be necessary when an initial phase fraction α k,L or α k,R is close to zero, in order to ensure the positivity of all the intermediate states densities. In the case of the wave ordering u k − a k τ k < u * 1 < u k < u k + a k τ k for instance, ensuring the positivity of τ k,R * in the regime ν k 1 (i.e. α k,R → 0) may require taking 0 < M < M 0 (ν k , ω k ) which implies Q k (u * 1 , W L , W R ) > 0. The precise choice of M in the interval (0, M 0 (ν k , ω k )) made in Section 4.5.2 of [12] to ensure the positivity of τ k,R * is still valid and will not be detailed here.
Now that the solution of the Riemann problem (S k )-(3.29) has been computed, we may apply Rankine-Hugoniot's jump relation to the momentum equation of (S k ) in order to determine the expression of the non-conservative product π * k ∆α k with respect to the given parameter u * 1 , thus defining the function . We obtain the following expression which is directly taken form equation (4.51) of [12]: where the function θ k is defined by: (3.37) Remark 3.14. This function θ k corresponds to an energy preserving solution (with Q k (u * 1 , W L , W R ) = 0) assuming all the intermediate densities are positive. If one has to dissipate energy in order to ensure the positivity of the densities, function θ k must be slightly modified (see [12], Sect. 4.5.2 for more details).

Solution to the fixed-point problem
Solving the fixed-point (3.20) amounts to re-coupling phase 1 with the other phases which have been decoupled for a separate resolution. This is done by equalizing the two expressions obtained for N k=2 π * k ∆α k in the first step (3.27) on the one hand and in the second step (3.36) (after summation over k = 2, . . . , N ) on the other hand. We obtain that u * 1 must solve the following scalar fixed-point problem : where the function Θ is defined by Θ(u) = θ 1 (u) + . . . + θ N (u). We have the following theorem which states a necessary and sufficient condition for the existence of solutions to the Riemann problem (3.6)-(3.12). One important fact is that this condition can be explicitly tested against the initial data (W L , W R ).
Theorem 3.15. Let be given a pair of admissible initial states (W L , W R ) ∈ Ω W × Ω W and assume that the parameter a k is such that τ k,L > 0 and τ k,R > 0 for all k = 1, . . . , N . The Riemann problem (3.6)-(3.12) admits solutions in the sense of Definition 3.5 if, and only if, the following condition holds: Remark 3.16. When simulating real industrial applications, the relaxation Riemann solver used for the convective effects will be associated with another step for the treatment of zero-th order terms enforcing the return to pressure, velocity (and possibly temperature) equilibrium between the phases. Hence, the pressure disequilibrium between the phases in the initial states is usually expected to be small, which yields small values of the quantities π 1 − π k . Hence, in most applications, condition (3.39) is expected to be satisfied. However, even away from pressure equilibrium, it is easy to observe that assumption (3.39) is always satisfied if the parameters (a k ) k=1,...,N are taken large enough. Indeed, denoting a = (a 1 , . . . , a N ), one can prove that: where C L and C R are two positive constants depending on (W L , W R ). For u * 1 ∈ (c L , u k ), we have, denoting ν k = . Hence, we obtain: It is not difficult to prove that the right hand side of this equality is positive which implies that θ k is strictly increasing on the interval (c L , u k ). Actually this exact computation has already been done in the framework of the Baer-Nunziato model (see [12], Eq. (4.65) for the details). A similar computation proves that θ k is also positive on the interval (u k , c R ). We obtain that all the functions θ k are continuous and strictly increasing on the open interval (c L , c R ) and so is Θ = θ 1 + . . . + θ N . The result of Theorem 3.15 follows from the intermediate value theorem.

The relaxation finite volume scheme and its properties
We now derive a finite volume scheme for the approximation of the entropy weak solutions of a Cauchy problem associated with system (2.1). For simplicity in the notations, we assume constant positive time and space steps ∆t and ∆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. We also introduce the discrete intermediate times t n = n∆t, n ∈ N. The approximate solution at time 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 : The Finite Volume relaxation scheme reads (see [12] for more details): where the numerical fluxes are computed thanks to the exact Riemann solver W Riem (ξ; W L , W R ) constructed for the relaxation system: The non-conservative part of the flux D * (W L , W R ) is defined in Definition 3.5 and the mappings M and P are given by: At each interface x j+ 1 2 , the relaxation Riemann solver W Riem (ξ; M(U n j ), M(U n j+1 )) depends on the family of relaxation parameters (a k ) k=1,...,N which must be chosen so as to ensure the conditions stated in the existence Theorem 3.15, and to satisfy some stability properties. Observe that one might take different relaxation parameters (a k ) k=1,...,N for each interface, which amounts to approximating system (2.1) by a different relaxation approximation at each interface, which is more or less diffusive depending on how large are the local parameters. Further discussion on the practical computation of these parameters is postponed to the appendix B, as well as the detailed description of the computation of the numerical fluxes F ± (U L , U R ).
Remark 3.17 (The method is valid for all barotropic e.o.s.). The Riemann solution W Riem (ξ; W L , W R ) only depends on the quantities u k , π k , τ k,L and τ k,R defined in (3.21) and on the left and right phase fractions α k,L and α k,R for k = 1, . . . , N . Indeed, the solution of the fixed-point problem (3.38) only depends on these quantities and so do the intermediate states (see (3.23) and (3.32)). Therefore, the dependence of the Riemann solution W Riem (ξ; W L , W R ) on the barotropic equation of state occurs only through the computation of π k (τ k,L , T k,L ) and π k (τ k,R , T k,R ). For (W L , W R ) = (M (U n j ), M(U n j+1 )), we have T k,L = (τ k ) n j , and T k,R = (τ k ) n j+1 and thus π k (τ k,L , T k,L ) = p k ((ρ k ) n j ) and π k (τ k,R , T k,R ) = p k ((ρ k ) n j+1 ) for all k = 1, . . . , N . These quantities can be computed for any barotropic e.o.s. at the beginning of each time step.
We may now state the following theorem, which gathers the main properties of this scheme, and which constitutes the main result of the paper. ∈ Ω U for all j ∈ Z).
-Conservativity. The discretizations of the partial masses α k ρ k , k = 1, . . . , N , and the total mixture momentum N k=1 α k ρ k u k are conservative. -Discrete energy inequalities. Assume that the relaxation parameters (a k ) n j+ 1 2 , k = 1, . . . , N satisfy Whitham's condition at each time step and each interface, i.e. that for all k = 1, . . . , N , n ∈ N, j ∈ Z, (a k ) n j+ 1 2 is large enough so that for all T k in the solution ξ → W Riem (ξ; M(U n j ), M(U n j+1 )). Then, the values U n j , j ∈ Z, n ∈ N, computed by the scheme satisfy the following discrete energy inequalities, which are discrete counterparts of the energy inequalities (2.9) and (2.10) satisfied by the exact entropy weak solutions of the model:

43)
and for k = 2, . . . , N: where for j ∈ Z, (α k ρ k E k u k + α k π k u k ) n ) is the right hand side trace of the phasic energy flux evaluated at x j+ 1 2 . Proof. A classical reformulation of the approximate Riemann solver allows to see that U n+1 j is the cell-average over C j of the following function at t = t n+1 : Hence, the positivity property on the phase fractions and phase densities is a direct consequence of Theorem 3.15 and Definition 3.5 which states the positivity of the densities in the relaxation Riemann solution. For this purpose, energy dissipation (3.16) across the u 1 -contact discontinuity may be necessary for enforcing this property when the ratio (or its inverse) is large for some j ∈ Z. The conservativity of the discretization of the mass equation is straightforward. That of the discretization of the total momentum is a consequence of the fact that u * 1 is the solution of the fixed-point problem (3.38). Let us now prove the discrete energy inequalities (3.44) for phases k = 2, . . . , N satisfied by the scheme under Whitham's condition (3.42). Assuming the CFL condition (3.40), the solution of (3.6) over [ ,xj ] (x)  1 ∆x ≥ 0. Since the initial data is at equilibrium: W(x, t n ) = M(U n j ) for all x ∈ C j (i.e. (T k ) n j is set to be equal to (τ k ) n j ) one has (α k ρ k E k )(M (U n j )) = (α k ρ k E k )(U n j ) according to Proposition 3.3. Applying the Rankine-Hugoniot jump relation to (3.47) across the line {(x, t), x = x j+ 1 2 , t > 0}, yields: Hence, since (Q k ) n j+ 1 2 ≥ 0, for the interface x j+ 1 2 , taking the trace of (α k ρ k E k u k + α k π k u k ) at 0 + instead of 0 − in (3.48) only improves the inequality. Furthermore, assuming that the parameter a k satisfies Whitham's condition (3.42), the Gibbs principle stated in (3.11) holds true so that: 1 ∆x Invoking the convexity of the mapping U → (α k ρ k E k )(U) (see [18]), Jensen's inequality implies that which yields the desired discrete energy inequality for phase k. The proof of the discrete energy inequality for phase 1 follows similar steps.

Numerical results
In this section, we present three test cases on which the performances of the relaxation scheme are illustrated. We only consider the three phase flow model (i.e. N = 3). In the first two test cases, the thermodynamics of the three phases are given by ideal gas pressure laws for k = 1, 2, 3: and we consider the approximation of the solutions of two different Riemann problems. In the third test case, we consider the simulation of a shock tube apparatus, where a gas shock wave interacts with a lid of rigid particles. This third case is also simulated with the three phase flow model although it is a two phase flow. The thermodynamics of the particle phase is given by a barotropic stiffened gas e.o.s. We recall that the scheme relies on a relaxation Riemann solver which requires solving a fixed-point problem in order to compute, for every cell interface x j+ 1 2 , the zero of a scalar function (see Eq. (3.38)). Newton's method is used in order to compute this solution. Usually, convergence is achieved within three iterations.

Test-case 1: A Riemann problem with all the waves
In this test-case, the thermodynamics of all three phases are given by barotropic e.o.s. (4.1) with the parameters given in Table 1. The wave pattern for phase 1 consists of a left-traveling rarefaction wave, a phase fraction discontinuity of velocity u 1 and a right-traveling shock. For phase 2 the wave pattern is composed of a left-traveling shock, the phase fraction discontinuity, and a right-traveling rarefaction wave. Finally, the wave pattern for phase 3 is composed of a left-traveling shock, the phase fraction discontinuity, and a right-traveling shock. The u 1 -contact discontinuity separates two regions denoted − and + respectively on the left and right sides of the discontinuity (see Fig. 1).  The relaxation scheme is compared with Rusanov's scheme, which is the only numerical scheme presently available for the three-phase flow model (see [6]). In Figure 2, 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). 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. In addition, the velocity of the contact discontinuity is not well estimated for the phase 2 variables with such a coarse mesh. We can also see that for the phase 2 and phase 3 variables, there are no oscillations as one can see for Rusanov's scheme: the curves are monotone between the intermediate states. The intermediate states for the phases 2 and 3 are not captured by Rusanov's scheme whereas the relaxation scheme gives a rather good estimation, even for the narrow state in phase 3 between the contact discontinuity and the right-traveling shock. These observations assess that, for the same level of refinement, the relaxation method is much more accurate than Rusanov's scheme.
A mesh refinement process has also been implemented in order to check numerically the convergence of the method, as well as its performances in terms of CPU-time cost. For this purpose, we compute the discrete L 1 -error between the approximate solution and the exact one at the final time T max = N ∆t = 0.05, normalized by the discrete L 1 -norm of the exact solution: where φ is any of the non-conservative variables (α 1 , α 2 , ρ 1 , u 1 , ρ 2 , u 2 , ρ 3 , u 3 ). The calculations have been implemented on several meshes composed of 100 × 2 n cells with n = 0, 1, . . . , 10 (knowing that the domain size is L = 1). In Figure 3, the error E(∆x) at the final time T max = 0.05, is plotted against ∆x in a log − log scale for both schemes. We can see that all the errors converge towards zero with at least the expected order of ∆x 1/2 . Note that ∆x 1/2 is only an asymptotic order of convergence, and in this particular case, one would have to implement the calculation on more refined meshes in order to reach the theoretically expected order of ∆x 1/2 . Figure 4 displays the error on the non-conservative variables with respect to the CPU-time of the calculation expressed in seconds for both the relaxation scheme and Rusanov's scheme. Each point of the plot corresponds to one single calculation for a given mesh size. One can see that, for all the variables except ρ 1 and u 1 , if one prescribes a given level of the error, the computational time needed to reach this error with Rusanov's scheme is higher than that needed by the relaxation scheme. On some variables, the gain of time can be spectacular. For instance, for the same error on the phase 1 fraction α 1 , the gain in computational cost is forty times when using the relaxation method rather than Rusanov's scheme which is a quite striking result. Indeed, even if Rusanov's scheme is known for its poor performances in terms of accuracy, it is also an attractive scheme for its reduced complexity. This means that the better accuracy of the relaxation scheme (for a fixed mesh) widely compensates for its (relative) complexity.

Test-case 2: A Riemann problem with a coupling between a single phase region and a mixture region
In this test-case, the thermodynamics of all three phases are still given by barotropic e.o.s. (4.1) with the parameters given in Table 3.
Here, we consider a Riemann problem in which two phases vanish in one of the initial states, which means that the corresponding phase fractions are equal to zero. For this kind of Riemann problem, the u 1 -wave separates a mixture region where the three phases coexist from a single phase region with the remaining phase.  Figure 3. Test-case 1: L 1 -Error with respect to ∆x for the relaxation scheme and Rusanov's scheme.  The solution is composed of a {u 3 − c 3 }-shock wave in the left-hand side (LHS) region where only phase 3 is present. This region is separated by a u 1 -contact discontinuity from the right-hand side (RHS) region where the three phases are mixed. In this RHS region, the solution is composed of a {u 1 + c 1 }-shock, a {u 2 + c 2 }-shock and a {u 3 + c 3 }-rarefaction wave.
In practice, the numerical method requires values of α 1,L and α 2,L that lie strictly in the interval (0, 1). Therefore, in the numerical implementation, we take α 1,L = α 2,L = 10 −10 . 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 1 and 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 ρ 1,L , u 1,L , ρ 2,L , u 2,L to the values on the right of the u 1 -contact discontinuity, which is coherent with the preservation of the Riemann invariants of this wave, and avoids the formation of fictitious acoustic waves for phases 1 and 2 in the LHS region. For the relaxation scheme, this choice enables to avoid oscillations of the phases 1 and 2  Figure 4. Test-case 1: L 1 -Error with respect to computational cost (in seconds) for the relaxation scheme (bullets, red line) and Rusanov's scheme (squares, blue line). density and velocity in the region where these phases are not present, as seen in Figure 5. However, some tests have been conducted that assess that taking other values of (ρ 1,L , u 1,L , ρ 2,L , u 2,L ) has little impact on the phase 3 quantities as well as on the phases 1 and 2 quantities where these phases are present.
As expected, we can see that for the same level of refinement, the relaxation method is more accurate than Rusanov's scheme. As regards the region where phases 1 and 2 do not exist, we can see that the relaxation scheme 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 α 1 or α 2 , since the solution approximated by Rusanov's scheme develops quite large oscillations. Moreover, the amplitude of these oscillations increases with the mesh refinement.

Test-case 3: Interaction of a gas shock wave with a lid of rigid particles
In this test-case, we consider a configuration where an incoming gas shock wave hits a cloud of spherical rigid particles. A sketch of the shock tube apparatus is given in Figure 6. The tube ranges from x = 0 m to x = 3.75 m and is closed at both ends. At the initial time t = 0 s, the cloud of particles lies between x = 2.97 m and x = 3.37 m. This test-case is adapted from the experimental setup presented in [7,8].
A gas pressure disequilibrium is initiated at x = 0.75 m. This initial pressure disequilibrium produces a left going gas rarefaction wave and a right-going gas shock wave. The rarefaction wave is soon reflected by the left wall while the shock first hits the lid of particles (producing a small left-going wave) and is then reflected by the right wall after crossing the lid of particles. In order to retrieve this behavior, four pressure transducers are located at stations S n for n = 1, . . . , 4 (see Fig. 6).
We decide to simulate this experiment with Rusanov's scheme and the relaxation scheme for the barotropic three-phase flow model. The particle phase has label 1 while the gas phase has label 2. Phase 3 is an "absent" phase, the statistical fraction of which is set to α 3 = 10 −10 everywhere. At time t = 0, the particle phase is residual outside the lid of particles and we set α 1 = 10 −10 for x / ∈ (2.97, 3.37). In the lid of particles, i.e. for x ∈ (2.97, 3.37) we set α 1 = 0.0104. At the initial time, all phases are at rest and in relative pressure equilibrium, with a space pressure disequilibrium at x = 0.75 m, thus for k = 1, . . . , 3:  Figure 7 displays the numerical approximations on the time interval (0, 0.01) of the total mean pressure P = 3 k=1 α k p k at stations S n for n = 1, . . . , 4. The obtained computations are run with a 5000 cell mesh with both the relaxation scheme and Rusanov's scheme. For both schemes, the observed curves reflect the expected behavior of the total mean pressure. At station S 1 , the pressure first jumps to a value P * when the right-going gas shock wave reaches x = 2.7 m around time t = 0.0037 s. The pressure at this location then remains steady until the reflection of the left-going rarefaction wave meets this position around time t = 0.0058 s which causes a decrease of the pressure until a second jump occurs (around time t = 0.0088 s) due to the reflection of the right-going gas shock wave which has hit the right wall boundary after crossing the lid of particles. Similar features are observed at stations S 2 and S 3 with the (expected) difference that the more rightward the station is located, the later the first pressure jump occurs and the sooner the second pressure jump occurs. If we now turn to station S 4 , a first pressure jump to the value P * is recorded around time t = 0.0056 s due to the right-going gas shock wave, soon followed by a second pressure jump (around t = 0.0059 s) to a value P * * due to the reflection of this wave on the right wall boundary. The pressure then remains steady until the location of station S 4 is reached by the reflection of the left-going rarefaction wave on the left wall boundary, causing a decrease in the pressure.
In Table 5, we compare the values of P * and P * * obtained respectively in the experiment [7,8], to those obtained with Rusanov's scheme and with the Relaxation scheme. We can see that the pressure values obtained with the barotropic three-phase flow model are slightly over estimated, a behavior that has already been observed in [6] on the same test-case.

Conclusion
We have extended to the N -phase compressible model developed in [17] the relaxation finite volume scheme designed in [12] for the barotropic Baer-Nunziato two phase flow model. The obtained scheme inherits the main properties of the scheme designed for the two phase framework. It applies to general barotropic equations of state (see Rem. 3.17). It is able to cope with arbitrarily small values of the statistical phase fractions. The approximated phase fractions and phase densities are proven to remain positive and a discrete energy inequality is proven under a classical CFL condition.
For N = 3, three test cases have been implemented which assess the good behaviour of the relaxation scheme. For the same level of refinement, the relaxation scheme is shown to be much more accurate than Rusanov's scheme, and for a given level of approximation error, the relaxation scheme is shown to perform much better in terms of computational cost than Rusanov's scheme. Moreover, contrary to Rusanov's scheme which develops strong oscillations when approximating vanishing phase solutions, the numerical results show that the relaxation scheme remains stable in such regimes. Given that Rusanov's scheme is the only numerical scheme presently available for the considered three phase flow model, the present paper therefore constitutes an improvement in this area.
Several natural extensions to this work can be considered. First of all, the scheme can be easily extended to the multidimensional framework. Indeed, the multidimensional version of the N -phase model (see [6] for the multi-D three phase model) is invariant under Galilean transformation and under frame rotation. Thus, the one-dimensional relaxation Riemann solver can still be used to obtain a finite volume scheme on two and three dimensional unstructured grids. An update of the cell unknown is obtained through a convex combination of 1D updates associated with the cell faces. These 1D updates are computed with the relaxation 1D scheme by considering local 1D Riemann problems in the orthogonal direction to the grid faces. Thanks to the convex combination, the positivity and entropy properties of the scheme are still valid for the multidimensional scheme under a natural CFL condition. We refer to [23] where this extension is detailed for the Baer-Nunziato two phase flow model, and where 2D-test cases have already been conducted. Another natural generalization is the extension of the scheme to higher order. A formally order two scheme can be obtained by considering a classical minmod reconstruction on the symmetrizing variable and a second order Runge-Kutta time scheme (see [13] for the two phase model). Such a procedure however does not ensure the preservation of the discrete energy inequality. Designing entropy satisfying second order numerical schemes is an open topic which is still under investigation. Finally, the extension of the relaxation numerical scheme to the multiphase flow model with non barotropic equations of state will be considered in a forthcoming paper. The key property of the relaxation scheme which allows this relatively simple extension is the existence of the discrete energy inequalities (3.43), (3.44). Thanks to these and to the second principle of thermodynamics which connects the phasic energies and the transported phasic entropies, one is able to extend the present Riemann solver to the full model with general e.o.s. through minor adaptations. This work has already been done for the Baer-Nunziato two phase model in [11]. The obtained scheme was shown to compare well with two of the most popular existing available schemes, namely Schwendeman-Wahle-Kapila's Godunov-type scheme [24] and Tokareva-Toro's HLLC scheme [27].
Appendix A. Eigenstructure of the relaxation system Proposition A.1. System (3.6) is weakly hyperbolic on Ω W in the following sense. It admits the following 4N − 1 real eigenvalues: We write π k instead of π k (τ k , T k ) in order to ease the notations. Defining M k = (u k − u 1 )/(a k τ k ), for k = 2, . . . , N , the matrices A, B 1 , . . . , B N and C 1 , . . . , C N are given as follows. A is a (N − 1) × (N − 1) diagonal matrix with u 1 on the diagonal. B 1 , . . . , B N are 3 × (N − 1) matrices and C 1 , . . . , C N are 3 × 3 matrices defined by: where δ p,q is the Kronecker symbol: for p, q ∈ N, δ p,q = 1 if p = q and δ p,q = 0 otherwise. Since A is diagonal and C k is R-diagonalizable with eigenvalues u k − a k τ k and u k + a k τ k and u k , the matrix A (W) admits the eigenvalues u 1 (with multiplicity N ), u 1 − a 1 τ 1 , u 1 + a 1 τ 1 and u k − a k τ k , u k + a k τ k and u k for k = 2, . . . , N . In addition, A (W) is R-diagonalizable provided that the corresponding right eigenvectors span R 4N −1 . The right eigenvectors are the columns of the following block matrix: where A is a (N − 1) × (N − 1) diagonal matrix defined by A = diag(1 − M 2 2 , . . . , 1 − M 2 N ). B 1 , . . . , B N are 3 × (N − 1) matrices and C 1 , . . . , C N are 3 × 3 matrices defined by: The first N − 1 columns are the eigenvectors associated with the eigenvalue u 1 . For k = 1, . . . , N , the N + 2(k − 1) -th, N + (2k − 1) -th and N + 2k -th columns are the eigenvectors associated with u k − a k τ k , u k + a k τ k and u k respectively. We can see that R(W) is invertible if and only if M k = 1 for all k = 2, . . . , N i.e. if and only if |u k − u 1 | = a k τ k for all k = 2, . . . , N . In particular, if for some k = 2, . . . , N , one has u k = u 1 , R(W) is still invertible. Denote (R j (W)) 1≤j≤4N −1 the columns of R(W). If 1 ≤ j ≤ N − 1, we can see that the (N + 1)-th component of R j (W) is zero. This implies that for all 1 ≤ j ≤ N − 1 R j (W) · ∇ W (u 1 ) = 0. Hence, the field associated with the eigenvalue u 1 is linearly degenerated. Now we observe that for all k = 1, . . . , N : R N +2(k−1) (W) · ∇ W (u k − a k τ k ) = 0, R N +(2k−1) (W) · ∇ W (u k + a k τ k ) = 0, This means that all the other fields are also linearly degenerated. ). The parameters (a k ) k=1,...,N , must be chosen large enough so as to satisfy several requirements: -In order to ensure the stability of the relaxation approximation, a k must satisfy Whitham's condition (3.42).