Stéphane Clain
email: clain@mip.ups-tlse.fr

David Rochette
email: rochette@moniut.univ-bpclermont.fr

First-and Second-order finite volume methods for the one-dimensional nonconservative Euler system

Keywords: Euler system, nonconservative, hyperbolic, finite volumes, Rusanov, MUSCL, resonant, Riemann problem

Gas flow in porous media with a nonconstant porosity function provides a nonconservative Euler system. We propose a new class of schemes for such a system for the one-dimensional situations based on nonconservative fluxes preserving the steady-state solutions. We derive a second-order scheme using a splitting of the porosity function into a discontinuous and a regular part where the regular part is treated as a source term while the discontinuous part is treated with the nonconservative fluxes. We then present a classification of all the configurations for the Riemann problem solutions. In particularly, we carefully study the resonant situations when two eigenvalues are superposed. Based on the classification, we deal with the inverse Riemann problem and present algorithms to compute the exact solutions. We finally propose new Sod problems to test the schemes for the resonant situations where numerical simulations are performed to compare with the theoretical solutions. The schemes accuracy (first-and second-order) is evaluated comparing with a nontrivial steady-state solution with the numerical approximation and convergence curves are established.

Introduction

Compressible fluid flows in porous media arise in many natural phenomena, for example the dust gas in volcanic eruption [START_REF] Pelanti | High-resolution finite volume methods for dusty gas jets and plumes[END_REF]. In industrial processes, porous media are also widely used for filtration or protection purposes. For example, grids (assimilated as a porous media) are employed to reduce gas velocity and absorb energy during an explosion in electrical switchgears consecutive to an electrical accident [START_REF] Rochette | Numerical investigations on the pressure wave absorption and the gas cooling interacting in porous filter, during an internal arc fault in a medium voltage cell[END_REF]. Based on the space-average methods, the representative elementary volume concept [START_REF] Bear | Introduction to modeling of transport phenomena in porous media[END_REF], the mathematical modelling derives from the Euler equations homogenization where a new force P ∂ x φ appears to take the porosity variation into account leading to a so-called nonconservative problem. The homogenizated Euler system is completed with the trivial equation ∂ t φ = 0 to provide an augmented nonconservative hyperbolic problem of the form

∂ t U + ∂ x F (U ) = G(U )∂ x φ.
(1)

Other problems, where a nonconservative term appears, cast in this general system as multiphase and multifluid flow problems [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF] where the phase fraction α replaces the porosity, compressible gas flow in duct [START_REF] Andrianov | On the solution to the Riemann problem for the compressible duct flow[END_REF] where the variational cross section identifies to the porosity and shallow-water [START_REF] Vázquez-Cendón | Improved treatment of source terms in upwind schemes for the shallow water equations in channels with irregular geometry[END_REF] where the variational bed elevation and breadth is similar to the porosity. It is well-known that the nonconservative term induced by the porosity variation leads to mathematical and numerical difficulties. On one hand, the nonconservative term gives rise to the product of a distribution with a discontinuous function in presence of shocks and the classical Rankine-Hugoniot condition does not hold any longer. On the other hand, the nonconservative term is responsible of specifics steady-state situations which have to be preserved by the numerical scheme leading to the so-called well-balanced property. An other difficulty is linked to the relative positions of the eigenvalues since the problem is not strictly hyperbolic. In particular, the eigenvalue deriving from the nonconservative term can be superposed with the eigenvalues of the Euler system leading to the resonant situation. It results a complex composition of shocks inside the interface defined by the characteristic field curves corresponding to the two superposed eigenvalues. Dal [START_REF] Maso | Definition and weak stability of a non-conservative product[END_REF] introduce a general framework to define the nonconservative product A(U )∂ x U for any variation bounded function U = U (x) where A(U ) is a n × n matrix, continuous with respect to U . Such a product is obtained using a path function ψ = ψ(s, U -, U + ) defined in the phase space linking any admissible state U -to any admissible state U + and the Rankine-Hugoniot condition at a shock of velocity σ writes (see also [START_REF] Le Floch | shock wave for nonlinear hyperbolic system in nonconservative form[END_REF])

1 0 A(φ(s, U -, U + )∂ s ψ(s, U -, U + ) ds = σ(U + -U -).
If A(U ) corresponds to the jacobian matrix of a flux function F then we recover the classical Rankine-Hugoniot condition and we have

1 0 A(φ(s, U -, U + )∂ s ψ(s, U -, U + ) ds = F (U + ) -F (U -).
independently of the choice of the path ψ. In the particular situation of system (1), the nonconservative part G(U )∂ x φ is relevant if the function φ is not constant.

In a precursor work, [START_REF] Toumi | A weak formulation of Roe's approximate Riemann solver[END_REF] proposes a Roe-like scheme for conservative system introducing the Roe average matrix between two states U L and U R

A(U L , U R )(U L -U R ) = 1 0 A(U L + s(U R -U L ))(U R -U L ) ds
where A(U ) = ∂ U F (U ) is the jacobian matrix of the flux F . An extension of the Toumi formula for nonconservative system has been proposed by [START_REF] Gosse | A well-balanced flux-vector splitting scheme designed for hyperbolic system of conservation laws with source terms[END_REF] and [START_REF] Gosse | A well-balanced scheme using non-conservative product designed for hyperbolic systems of conservative laws with source terms[END_REF] for the Euler system and [START_REF] Parés | On the well-balance property of Reo's method for nonconservative hyperbolic systems[END_REF] for the shallow-water problem where A(U ) is constituted of the jacobian matrix of F (U ) and the nonconservative contribution G(U ). Such methods where the jacobian matrix and the nonconservative flux are combined (the Q-methods [START_REF] Castro | A Q-scheme for a class of systems of coupled conservation laws with source term. Applications to a two-layer 1-D shallow water system[END_REF]) are employed in numerous applications like shallow-water [START_REF] Vázquez-Cendón | Improved treatment of source terms in upwind schemes for the shallow water equations in channels with irregular geometry[END_REF] or porous media [START_REF] Rochette | Numerical scheme to compute a compressible gas flow in variable porosity media[END_REF][START_REF] Rochette | Two-dimensional computation of gas flow in a porous bed characterized by a porosity jump[END_REF]. Extensions to higher order scheme have been proposed using a ENO-WENO reconstruction or the MUSCL technique [START_REF] Gallardo | High order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. Application to shallow-water systems[END_REF][START_REF] Gascón | Construction of second-order TVD schemes for nonhomogeneous hyperbolic conservative laws[END_REF][START_REF] Vukovic | ENO and WENO schemes with exact conservation property for one-dimensional shallow water equations[END_REF].

In the problems mentioned above, it appears that the nonconservative term corresponds to a linearly degenerate simple wave and the Riemann invariant is employed in place to the Rankine-Hugoniot condition. In particular, the Riemann problem resolution can be performed using in one hand the classical simple waves associated to the eigenvalues of the conservative part and in the other hand the Riemann invariants associated to the nonconservative term. Based on this strategy, the Riemann problem resolution for the scalar hyperbolic equation with source term has been done by [START_REF] Noussair | Analysis of nonlinear resonance in conservation law with point sources and well-balanced scheme[END_REF]. The shallow-water problem with variable bed has been done by [START_REF] Alcrudo | Exact solutions to the Riemann problem of the shallow water eqautions with a bottom step[END_REF], [START_REF] Chinnayya | A new general Riemann solver for the shallow water equations, with frinction and topography[END_REF], [START_REF] Chinnayya | A well-balanced numerical scheme for shallow-water equations with topography: resonance phenomena[END_REF] and [START_REF] Noussair | Riemann problem with nonlinear resonance effects and wellbalanced godunov scheme for shallow fluid past an obstacle[END_REF]. The Riemann problem for the compressible duct flow has been considered by [START_REF] Andrianov | Analytical and numerical investigation of two-phase flow[END_REF] and [START_REF] Andrianov | On the solution to the Riemann problem for the compressible duct flow[END_REF] (see also [START_REF] Le Floch | The riemann problem for fluid flows in a nozzle with discontinuous cross-section[END_REF] for the isentropic case). The Riemann problem for the two-phase flow has been studied by [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF] and [START_REF] Deledicque | An exact Riemann solver for compressible two-phase flow models containing non-conservative products[END_REF]. We also mention the work of [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF] where the authors study the resonant situation in a general framework.

The delicate point is the construction of the simple wave associated to the nonconservative term. For a given left state, there exists zero, one or two right states such that the Riemann invariants are preserved (see [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF] for the general case). In the situation when two solutions are possible, we usually have a subsonic (subcritical) flow and a supersonic (supercritical) flow thus a new criterion (the continuation criterion [START_REF] Chinnayya | A well-balanced numerical scheme for shallow-water equations with topography: resonance phenomena[END_REF] or the admissible criterion [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF]; [START_REF] Isaacson | Nonlinear resonance in systems of conservation laws[END_REF]) is required to select a solution. Moreover, the particular situation where two eigenvalues merge gives rise to a complex situation where three or more states coexist in the interface at the same moment. A complete description is proposed by [START_REF] Chinnayya | A new general Riemann solver for the shallow water equations, with frinction and topography[END_REF], [START_REF] Chinnayya | A well-balanced numerical scheme for shallow-water equations with topography: resonance phenomena[END_REF] for the shallow-water case while Le [START_REF] Le Floch | The riemann problem for fluid flows in a nozzle with discontinuous cross-section[END_REF] give the solution of the Euler system with isentropic flow problem for the resonant situation. Up to our knowledge, there is not a complete description of the admissible configurations for the Riemann problem for the Euler system.

From a numerical point of view, an important step has been done with the introduction by [START_REF] Bermúdez | Upwind methods for hyperbolic conservation laws with source term[END_REF] of the conservative property: the C-property while [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF] propose a similar notion : the well-balanced condition, for the scalar hyperbolic problem with a source term. We find also the same idea in [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF] for the two-phase flow problem where the authors propose a numerical nonconservative contribution which preserve the pressure and the velocity of the steady-state solutions. All these conditions tend to the following criterion : the Riemann invariants associated to the nonconservative term have to be numerically preserved. A particular case is the fluid at rest where the steady-state configuration has to be maintained numerically. [START_REF] Parés | On the well-balance property of Reo's method for nonconservative hyperbolic systems[END_REF] propose a general definition of well-balanced scheme with order k based on the scheme capacity to preserve the Riemann invariant up to and order ∆x k where ∆x is the characteristic length of the mesh.

In this paper, we propose a new class of scheme based on the criterion of well-balanced scheme introduced by [START_REF] Parés | On the well-balance property of Reo's method for nonconservative hyperbolic systems[END_REF]. The basic idea is to employ a generic numerical flux to solve the conservative part and we construct the nonconservative flux such that the well-balanced criterion is satisfied. To this end, we propose the notion of intermediate state we use to construct the numerical methods and we apply the technique to the classical Rusanov flux. A second-order scheme based on a MUSCL resconstruction is proposed where the porosity is decomposed into regular and discontinuous functions. The regular contribution is treated as a classical source term while the discontinuous contribution is taken into account with the nonconservative flux.

To test the method, we have considered the Riemann problem and we present an new approach based on the configuration identification where we propose a classification of the Riemann problem solutions. Indeed, the main difficulty with the nonconservative problems is the nonstrictly hyperbolicity character of the system and eigenvalues can cross or merge. We develop a technique where we describe all the admissible configurations which can appear and we use the classification to solve the inverse Riemann problem introduced by [START_REF] Andrianov | On the solution to the Riemann problem for the compressible duct flow[END_REF]. In particularly we are able to compute complex situations including the resonant cases and prove that rarefaction can only reach the sonic point from the lower porosity side. New Sod tube tests corresponding to particular difficult situations are proposed to check the solvers. Numerical tests have been performed to compare the approximated solution with the exact solution computed with the inverse Riemann problem algorithm. Other simulations of a nontrivial steady-state solution are also proposed to measured the method accuracy.

The paper is organized as follows. In section 2, we present the nonconservative Euler system and the properties of the steady-state solutions. Section 3 and 4 are devoted to the nonconservative schemes based on the Rusanov flux and its extension to the second-order using the MUSCL method and the porosity splitting principle. In section 5 and 6, we establish a classification of all the configurations of the Riemann problem, including the resonant particular situations, and a generic method to compute the exact solutions: the Riemann inverse problem. Finally in section 7, we present in a first part a selection of Sod problems, scanning all the configurations, to test the numerical schemes in comparison with the exact solutions and in a second part, we evaluate the schemes accuracy in the framework of a steady-state solution with regular porosity.

Nonconservative Euler system

Notations

We consider the one-dimensional nonconservative Euler system:

∂ t   φρ φρu φE   + ∂ x   φρu φρu 2 + φP φu(E + P )   =   0 P ∂ x φ 0   , (2)
where φ stands for the porosity, ρ the gas density, u the velocity, P the pressure and E the total energy composed of the internal energy e and the kinetic energy: E = ρ( 1 2 u 2 + e). In addition, to close the system, we use perfect gas law P = (γ -1)ρe with γ ∈]1, 3[. Note that the numerical scheme we shall present also deals with more complex pressure law given by function P = P (ρ, e) which depends on the density and the internal energy.

The conservative quantities (or conservative state) are represented by vector U = (φρ, φρu, φE) which belongs to the conservative variable phase space Ω c ⊂ R + × R × R + while the physical (primitive) variables vector V = (φ, ρ, u, P ) belongs to the physical variable phase space

Ω p ⊂ R + × R + × R × R + . We have a one to one mapping (φ, U ) → V (φ, U ) such that V (x, t) = V (φ(x), U (x, t)) with inverse function V → (φ, U (V )) such that U (V (x, t)) = U (x, t).
In the sequel, we shall drop the hat symbol for the sake of simplicity.

The flux vector is given by

F (U ) =   F ρ (U ) F u (U ) F e (U )   =   φρu φρu 2 + φP φu(E + P )   , (3)
while the nonconservative term writes

G(U ) =   G ρ (U ) G u (U ) G e (U )   =   0 P 0   . (4) Since U = U (V ) is a function of V , we adopt the notation F α (V ) = F α (U (V )) and G α (V ) = G α (U (V )) with α = ρ, u, e.

Steady-state solutions

The main challenge for the numerical schemes to solve nonconservative system is the steady-state solutions preservation. Let us consider a regular stationary solution U (x) of system (2), we then have

d dx F (U (x)) = G(U (x)) dφ dx . (5)
Assume that φ = φ(x) is a strictly monotone function on the interval [x -, x + ] with respect to x and ranges between φ -= φ(x -) and φ + = φ(x + ) and we denote U -= U (x -) and U + = U (x + ). We change the variable x by the variable φ where we state

U (x) = U (φ(x)) = U (φ). Function U (φ) is then solution of the system d dφ F ( U (φ)) = G( U (φ)). (6)
We drop the tilde symbol for the sake of simplicity and we deduce that U belongs to an integral curve parameterized with φ. For a given U -and φ -, we have (at least locally) a unique curve W(φ; φ -, U -) = W(φ; V -) solution of (6) with W(φ -; V -) = U -. Since U + belongs to the integral curve, we also have

W(φ + ; V -) = U + .
The main advantage to use φ as a variable is that relation (6) still holds even if φ(x) is a discontinuous function of x. The ability to handle the φ discontinuity is of crucial importance for the Riemann problem. Indeed, function φ is constant in the open sets x < 0 and x > 0 and the nonconservative problem turns to a conservative one as noticed by [START_REF] Andrianov | On the solution to the Riemann problem for the compressible duct flow[END_REF]. The nonconservative part only acts at the interface x = 0 where φ jumps from φ -to φ + . We smooth the φ discontinuity using the following linear approximation (see [START_REF] Chinnayya | A well-balanced numerical scheme for shallow-water equations with topography: resonance phenomena[END_REF] or [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF])

φ ε (x) =    φ ε (x) = φ -, x < -ε, φ ε (x) = ((x -ε)φ -+ (x + ε)φ + )/(2ε), x ∈ [-ε, +ε], φ ε (x) = φ + , x > +ε.
Substituting φ by φ ε , we get a smooth steady-state solution U ε (x) for equation ( 5) and an equivalent smooth solution U ε (φ) for equation ( 6) for any ε > 0. The point is that we can take the limit in relation ( 6) when ε converges to 0.

In the case of the nonconservative Euler problem, system (6) is equivalent to the three relations

(a) φρu = D, (b) d dφ (φρu 2 + φP ) = P, (c) u 2 + 2γ γ -1 P ρ = H, (7)
where D and H are constants. If u = 0, we get that P and ρ are constant.

Assuming now that u = 0 then we eliminate u in the two last equations and equation (7-b) provides

- D 2 (φρ) 3 d dφ (φρ) + 1 ρ dP dφ = 0, while differentiation of equation (7-c) yields - D 2 (φρ) 3 d dφ (φρ) + γ γ -1 ρ dP dφ -P dρ dφ ρ 2 = 0.
We combine the two expressions and we get

dP dφ = γ P ρ dρ dφ
and we obtain the relation P ρ γ = S by integration where S > 0 is a constant. For a given initial state V -, the curve W(φ, V -) is implicitly given by the three relations

φρu = D, u 2 + 2γ γ -1 P ρ = H, P ρ γ = S, (8)
where constants D, H and S are determined by the initial condition V -. We deduce from relation (8) an implicit relation P = P (φ) given by

D 2 φ 2 + 2γ γ -1 P P S 1 γ = H P S 2 γ . ( 9
)

Nonconservative scheme

Domain R is uniformly divided with cells

K i = [x i-1/2 , x i+1/2 ], i ∈ Z of
length ∆x where we set x i-1/2 = i∆x and

x i = x i-1/2 + ∆x 2 is the cell center.
We denote by

φ i = 1 |K i | K i φ(x) dx
, by α n i an approximation of the mean value of α on cell K i at time t n for the other physical variables α = ρ, u, P and U n i (resp. V n i ) for an approximation of U (resp. V ) at time t n on the cell K i . We consider generic schemes of the form

U n+1 i = U n i - ∆t ∆x F i+1/2 + G - i+1/2 -F i-1/2 -G + i-1/2 . ( 10
) Flux F i+1/2 = F(V i , V i+1 ) is the conservative numerical flux across interface x i+1/2 while G - i+1/2 = G -(V i , V i+1 ) and G + i-1/2 = G + (V i-1 , V i )
represent the nonconservative contribution across the interface consecutive to the φ space variation. Real functions F α , G -,α and G +,α α = ρ, u, e represent the three components of the flux vectors.

Remark 1. Since G ± (V L , V R ) represent the nonconservative contributions due to the φ variation, the numerical nonconservative fluxes vanish when the porosity is constant hence we state that for any physical state V 0

G -(V 0 , V 0 ) = G + (V 0 , V 0 ) = 0.
For any physical state V j ∈ Ω p , corresponding to the cells K j , j = i-1, i, i+1, we introduce the three points operator H

H(V i-1 , V i , V i+1 ) = F(V i , V i+1 ) + G -(V i , V i+1 ) -F(V i-1 , V i ) -G + (V i-1 , V i )
(11) hence the three points scheme reads

U n+1 i = U n i - ∆t ∆x H(V n i-1 , V n i , V n i+1 ) with V n j = V (φ j , U n j ), j = i -1, i, i + 1.

Well-balanced flux criterion

Following the idea of [START_REF] Parés | On the well-balance property of Reo's method for nonconservative hyperbolic systems[END_REF] and [START_REF] Parés | Numerical methods for nonconservative hyperbolic systems: a theoritical framework[END_REF], we introduce a definition of well-balanced schemes. We emphase that the definitions we propose are based on the φ parameterization of the intergral curve instead of the x parameterization as it is proposed by [START_REF] Parés | On the well-balance property of Reo's method for nonconservative hyperbolic systems[END_REF]. The motivation is that we want to define well-balanced schemes even if φ is discontinuous in space.

Definition 1 (well-balanced scheme). Let V i ∈ Ω p and W(φ; V i ) the parameterization of the integral curve with φ. For any φ j , j = i-1 and j = i+1, we set U j = W(φ j ; V i ) and V j = V (φ j , U j ).

(a) The scheme is exactly well-balanced for the state V i if for any φ j , j = i -1 and j = i + 1

H(V i-1 , V i , V i+1 ) = 0. (b)
The scheme is well-balanced with order p for the state V i if for any φ j , j = i -1 and j = i + 1

H(V i-1 , V i , V i+1 ) = O(max(|φ i -φ i-1 |, |φ i -φ i+1 |) p+1 .
(c) The scheme is exactly well-balanced (resp. well-balanced with order p) if it is exactly well-balanced (resp. well-balanced with order p) for any

V i ∈ Ω p .
Note that we recover the definition proposed by [START_REF] Parés | On the well-balance property of Reo's method for nonconservative hyperbolic systems[END_REF] in the case where φ is a local regular function in space i.e. φ j = φ(x j ), j = i-1, i, i+ 1. Indeed, in this case, we have

|φ i -φ i-1 | = O(∆x) and |φ i -φ i+1 | = O(∆x), hence H(V (x i-1 ), V (x i ), V (x i+1 )) = O(∆x) p+1 .
The well-balanced definition induces a nonconservative flux consistency with the nonconservative term for steady-state solution. To this end, we assume that the numerical flux F(V 1 , V 2 ) and the physical flux F (U ) are regular functions (at least C 1 ) and we denote by

∇ 1 F(V 1 , V 2 ) and ∇ 2 F(V 1 , V 2 )
the Jacobian matrices with respect to V 1 and V 2 respectively.

Proposition 1. Let φ(x) be a regular function on the neighbouring of x i , i.e. on the intervals [x i-1 , x i+1 ] for ∆x small enough, and U i ∈ Ω c a given conservative state. We denote by U (x) = W(φ(x); V i ) the local regular solution of equation ( 5) and we set φ j = φ(x j ),

U j = U (x j ) = W(φ j ; V i ), V j = V (U j , φ j ) for j = i -1, i + 1.
Assume that H is well-balanced with order p ≥ 1 then we have

1 ∆x G -(V i , V i+1 ) -G + (V i-1 , V i ) = -G(x i )∂ x φ(x i ) + O(∆x). (12)
Proof. To prove the proposition, we recall the notation where we distinguish the vector U (x) with the function U (V ). Thanks to the consistency condition F(V, V ) = F ( U (V )), we deduce

∇ 1 F(V, V ) + ∇ 2 F(V, V ) = ∇ U F ( U (V )).∇ V U (V ).
On the one hand, since we have a p th order scheme with p ≥ 1, we can write

1 ∆x F(V i , V i+1 ) -F(V i-1 , V i ) + 1 ∆x G -(V i , V i+1 ) -G + (V i-1 , V i ) = O(∆x) p .
On the other hand, we introduce vector

F (V i , V i ) such that 1 ∆x F(V i , V i+1 ) -F(V i-1 , V i ) = 1 ∆x F(V i , V i+1 ) -F(V i , V i ) + F(V i , V i ) -F(V i-1 , V i ) .
Writing the Taylor expansion for

F(V i , V i+1 ) = F(V (x i ), V (x i + ∆x)) yields 1 ∆x F(V i , V i+1 ) -F(V i , V i ) = ∇ 2 F(V i , V i )∂ x V (x) |x=x i + O(∆x).
In the same way, we have

1 ∆x F(V i , V i ) -F(V i-1 , V i ) = ∇ 1 F(V i , V i )∂ x V (x) |x=x i + O(∆x).
Hence we get with the consistency property

1 ∆x F(V i , V i+1 ) -F(V i-1 , V i ) = ∇ U F ( U (V )) |V =V i ∇ U (V ) |V =V i ∂ x V (x) |x=x i = ∂ x (F (U (x)) |x=x i + O(∆x).
We then deduce

1 ∆x G -(V i , V i+1 ) -G + (V i-1 , V i ) + ∂ x (F (U (x)) |x=x i = O(∆x).
Since U (x) belongs to the integral curve, relation ( 5) is satisfied and we get

∂ x F (U (x)) |x=x i = G(U ) |U =U i ∂ x φ(x) |x=x i ,
thus relation ( 12) holds.

We now turn to the discontinuous case for the Riemann problem with initial conditions U (x, 0) = U L , φ(x) = φ L for x < 0 and U (x, 0) = U R , φ(x) = φ R for x > 0. To deal with a steady-state problem, we assume that U R = W(φ R ; V L ) i.e. U L and U R belong to the same integral curve. Since we have a stationary solution, we reduce the study on the cells K 0 = [-∆x, 0] and K 1 = [0, ∆x] which share the interface x = 0 where U

(x) = U L on K 0 while U (x) = U R on K 1 .
Proposition 2. The following consistancy relations hold: (a) If H is a p th order well-balance scheme, then

G + (V L , V R ) -G -(V L , V R ) = φ R φ L G(W(φ; V L )) dφ + O(|φ R -φ L |) p+1 . (13) (b) If H is an exact well-balanced scheme, then G + (V R , V L ) -G -(V L , V R ) = φ R φ L G(W(φ; V L )) dφ. ( 14
)
Proof. We only deal with the situation where H is a p th order well-balanced scheme, the case of an exact well-balanced scheme is similar and provide relation ( 14). On the cell K 0 , we write

H(V L , V L , V R ) = F(V L , V R ) + G -(V L , V R ) -F(V L , V L ) -G + (V L , V L ), while we have on cell K 1 H(V L , V R , V R ) = F(V R , V R ) + G -(V R , V R ) -F(V L , V R ) -G + (V L , V R ).
Using remark 1 and adding the two relations, the definition of a p th order well-balanced scheme yields

F(V L , V R ) + G -(V L , V R ) -F(V L , V L ) + F(V R , V R ) -F(V L , V R ) -G + (V L , V R ) = O(|φ R -φ L |) p+1 .
Since U L and U R belong to the same integral curve, we have by integration of equation ( 6)

F(V R , V R ) -F(V L , V L ) = F (U R ) -F (U L ) = φ R φ L G(W(φ; V L )) dφ.
Hence relation (13) holds.

The intermediate state

To construct explicitly the nonconservative schemes, we propose a new technique (up to our knowledge) based on intermediate states which provide expressions for the nonconservative flux G ± . More precisely, let us consider three physical states

V i-1 = V (φ i-1 , U i-1 ), V i = V (φ i , U i ) and V i+1 = V (φ i+1 , U i+1 ) given on cell K i-1 , K i and K i+1 respectively such that U i-1 = W(φ i-1 ; V i ), U i+1 = W(φ i+1 ; V i ).
For the couple (φ i-1 , φ i ), we introduce an intermediate value φ * i-1/2 and its corresponding conservative state

U * i-1/2 = W(φ * i-1/2 ; V i ) and physical state V * i-1/2 = V (φ * i-1/2 , U * i-1/2 ) at the interface x i-1/2 .
In the same way, we define φ * i+1/2 , U * i+1/2 and V * i+1/2 at the interface x i+1/2 . The goal is to choose

φ * i-1/2 and φ * i+1/2 such that F(V i , V i±1 ) -F (V * i±1/2 ) = O(|φ i -φ i±1 | p+1
H(V i-1 , V i , V i+1 ) = (15) F(V i , V i+1 ) -F (V * i+1/2 ) + G -(V i , V i+1 ) + φ * i+1/2 φ i G(W(φ; V i )) dφ - F(V i-1 , V i ) -F (V * i-1/2 ) + G + (V i-1 , V i ) + φ * i-1/2 φ i G(W(φ; V i )) dφ .
Proof. Since the conservative states U * i-1/2 and U * i+1/2 belong to the same integral curve W(φ; V i ), we have

F (V * i-1/2 ) -F (V i ) = F (U * i-1/2 ) -F (U i ) = φ * i-1/2 φ i G(W(φ; V i )) dφ, F (V * i+1/2 ) -F (V i ) = F (U * i+1/2 ) -F (U i ) = φ * i+1/2 φ i G(W(φ; V i )) dφ.
On the other hand, definition of the numerical scheme H writes

H(V i-1 , V i , V i+1 ) = F(V i , V i+1 ) + G -(V i , V i+1 ) -F(V i-1 , V i ) -G + (V i-1 , V i ) = F(V i , V i+1 ) -F (V i ) + G -(V i , V i+1 ) - F(V i-1 , V i ) -F (V i ) + G + (V i-1 , V i ) .
Substitution of F (V i ) in the last relation provides relation (15). Note that no regularity assumption has been made to obtain relation (15). Assume that physical states V L and V R correspond to a steady-state discontinuous solution such that U R = W(φ R ; V L ). Relation (15) then writes

H(V L , V L , V R ) = F(V L , V R ) -F (V * ) + φ * φ L G(W(φ; V i )) dφ + G -(V L , V R ) - F(V L , V L ) -F (V L ) + φ L φ L G(W(φ; V i )) dφ + G + (V L , V L )
where φ * is an intermediate value between φ L and φ R and

U * = W(φ * ; V L ), V * = V (φ * , U * ).
After simplification, we have on cells K 0 and K 1

H(V L , V L , V R ) = F(V L , V R ) -F (V * ) + G -(V L , V R ) + φ * φ L G(W(φ; V i )) dφ, (16)
H(V L , V R , V R ) = F(V L , V R ) -F (V * ) + G + (V L , V R ) + φ * φ R G(W(φ; V i )) dφ.
(17)

The Rusanov nonconservative scheme

We elaborate a well-balanced scheme based on a very simple flux initially proposed by [START_REF] Rusanov | Calculation of interaction of nonsteady-shock waves with obstacles[END_REF]. We then construct the nonconservative flux such that the well-balanced property is achieved using an adapted intermediate state. [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF] use the same flux for the multiphase flows while [START_REF] Andrianov | Analytical and numerical investigation of two-phase flow[END_REF] introduces the Rusanov flux for the Baer-Nunziato problem. The main advantages are the robustness and the simplicity of such a flux which can be used with complex pressure law. The main drawbacks are the important diffusion generated by the diffusion term to stabilize the scheme and contact discontinuities are not well-determined since the solver does not take into account the intermediate waves [START_REF] Toro | restoration of the contact surface in the HLL-Riemann solver[END_REF]. It is noticeable that second-order schemes like MUSCL technique strongly reduces the diffusion effect and the Rusanov flux usage is relevant in a second-order scheme perspective.

The Rusanov conservative flux F

We first begin with the expression of the conservative part of the flux which is an adaptation of the original Rusanov flux for the Euler problem with porosity. The original Rusanov flux F(V L , V R ) for the classical Euler system (i.e. with φ = 1) associated to the physical flux F (V ) writes

F(V L , V R ) = F (V L ) + F (V R ) 2 -λ(U R -U L ) where λ = sup(c L + |u L |, c R + |u R |),
with c L (resp. c R ) the sound velocity associated to state V L (resp. V R ). We shall modify the classical Rusanov scheme for the following reason. Consider a Riemann problem with V L = (φ L , ρ, 0, P ) and V R = (φ R , ρ, 0, P ) with ρ, P ∈]0, +∞[. Since U L and U R belong to the same integral curve, we have a steady-state solution.

The first and third components of the nonconservative contribution G are null therefore the first and third components of the nonconservative numerical flux G ± vanish and, for the mass equation, the scheme writes on cell K 0

ρ 1 0 = ρ - ∆t ∆x c(φ R ρ -φ L ρ)
with c the sound velocity associated to the configuration. If the porosity is not constant i.e. φ L = φ R , the density changes due to the viscosity contribution and the steady-state solution is not preserved. A similar situation also happens with the energy equation. To preserve the steady-state solution, we introduce the following Rusanov flux adapted to the porosity variation situation:

F(V L , V R ) = F (V L ) + F (V R ) 2 -V(V L , V R ), (18)
with

V(V L , V R ) = λφ LR   ρ R -ρ L ρ R u R -ρ L u L E R -E L   where φ LR = max(φ L , φ R ). Since φ ∈]0, 1],
we can also choose φ LR = 1.

The Rusanov nonconservative flux G ±

We now construct the nonconservative flux G ± in order to preserve the steady-state solutions. Since the nonconservative contribution only acts on the impulsion equation, we have G ±,α = 0 for α = ρ and e. We now turn to the construction of the nonconservative flux G ±,u . Relation ( 15) suggests that we have to define the intermediate states such that the conservative contribution writes

F(V i , V i+1 ) -F (V * i+1/2 ) - F(V i-1 , V i ) -F (V * i-1/2 ) = O(∆φ) p+1 , while the nonconservative contribution writes φ * i+1/2 φ i G(W(φ; V i )) dφ + G -(V i , V i+1 ) - φ * i-1/2 φ i G(W(φ; V i )) dφ + G + (V i-1 , V i ) = O(∆φ) p+1 .
To explicite the scheme, we have to fix the intermediate values φ * i-1/2 and φ * i+1/2 . We propose here to set

φ * i-1/2 = φ i-1 + φ i 2 , φ * i+1/2 = φ i+1 + φ i 2 . ( 19
)
and we define the Rusanov nonconservative fluxes by:

G -(V i , V i+1 ) = - φ i+1 -φ i 2   0 P i 0   , G + (V i-1 , V i ) = + φ i -φ i-1 2   0 P i 0   .
(20)

Nonconservative Rusanov scheme order : the regular case

To study the nonconservative scheme order, we first deal with the regular situation where we assume that (φ(x), U (x)) is a regular steady-state solution.

Lemma 1. We have the following error order

F(V i , V i+1 ) -F (V * i+1/2 ) -F(V i-1 , V i ) -F (V * i-1/2 ) = O(∆φ 2 ) ( 21
)
with ∆φ = max(|φ i+1 -φ i |, |φ i -φ i-1 |).
Proof. The conservative part of the flux is composed of a difference of flux

A i+1/2 = F (V i+1 ) + F (V i ) 2 -F (V * i+1/2 )
and the viscous term

B i+1/2 = -λ i+1/2   ρ i+1 -ρ i ρ i+1 u i+1 -ρ i u i E i+1 -E i   with U * i+1/2 = W(φ * i+1/2 ; V i ) and V * i+1/2 = V (φ * i+1/2 , U * i+1/2 ). Since U i+1 = W(φ i+1 ; V i ), we have F α (V i ) = F α (V i+1 ) = F α (V * i+1/2
) for α = ρ, e and we obtain the simpler expression

A i+1/2 =    0 F u (V i ) + F u (V i+1 ) 2 -F u (V * i+1/2 ) 0    . We denote ∆φ i+1/2 = φ i+1 -φ i . The Taylor expansion of function W(φ; V i ) gives U * i+1/2 -U i = W φ i + φ i+1 2 ; V i -W(φ i ; V i ) = - dW dφ φ i + φ i+1 2 ; V i ∆φ i+1/2 2 + O(∆φ 2 i+1/2 ), U * i+1/2 -U i+1 = W φ i + φ i+1 2 ; V i -W(φ i+1 ; V i ) = + dW dφ φ i + φ i+1 2 ; V i ∆φ i+1/2 2 + O(∆φ 2 i+1/2 ). We deduce that U * i+1/2 -U i = O(∆φ i+1/2 ), U * i+1/2 -U i+1 = O(∆φ i+1/2
) and the impulsion flux difference writes

F u (V i+1 ) + F u (V i ) 2 -F u (V * i+1/2 ) = 1 2 [F u (U i+1 ) -F u (U * i+1/2 )] + 1 2 [F u (U i ) -F u (U * i+1/2 )] = 1 2 ∇ U F u (U * i+1/2 )(U i+1 -U * i+1/2 ) + 1 2 ∇ U F u (U * i+1/2 )(U i -U * i+1/2 ) + O(∆φ 2 i+1/2 ).
We finally obtain

F u (V i+1 ) + F u (V i ) 2 -F u (V * i+1/2 ) = O(∆φ 2 i+1/2 ), (22)
hence we get A i+1/2 = O(|∆φ| 2 ). In the same way, we have

A i-1/2 = O(|∆φ| 2 ).
We now focus on the viscous contribution

B i+1/2 -B i-1/2 . A simple calculation gives B i+1/2 -B i-1/2 = -λ i+1/2   ρ i+1 -ρ i ρ i+1 u i+1 -ρ i u i E i+1 -E i   + λ i-1/2   ρ i -ρ i-1 ρ i u i -ρ i-1 u i-1 E i -E i-1   = (λ i-1/2 -λ i+1/2 )   ρ i+1 -ρ i ρ i+1 u i+1 -ρ i u i E i+1 -E i   - λ i-1/2   ρ i+1 -2ρ i + ρ i-1 ρ i+1 u i+1 -2ρ i u i + ρ i-1 u i-1 E i+1 -2E i + E i-1   ,
where

λ i+1/2 = max(|u i | + c i , |u i+1 | + c i+1 ) and λ i-1/2 = max(|u i | + c i , |u i-1 | + c i-1 ). Since function max(α, β) is a 1-lipschitzian function, we get (λ i-1/2 - λ i+1/2 ) = O(∆φ).
On the other hand, we have

α i -α i-1 = O(∆φ) and α i+1 -2α i + α i-1 = O(∆φ 2 ) for α = ρ, ρu, E.
We finally obtain

B i+1/2 -B i-1/2 = O(∆φ 2 ) ( 23
)
and estimate (21) holds from relations ( 22)-( 23).

Remark 2. Note that relation ( 22) holds even if the steady-state is not regular while estimate ( 23) is obtained under the regularity assumption. It is of numerical interest to see that the order discrepancy comes from the viscosity stability terms B i+1/2 and B i-1/2 when we deal with irregular solutions .

We now treat the nonconservative numerical flux.

Lemma 2. We have the following error order

φ * i+1/2 φ i G(W(φ; V i )) dφ + G -(V i , V i+1 ) - φ * i-1/2 φ i G(W(φ; V i )) dφ + G + (V i-1 , V i ) = O(∆φ) 2 . ( 24
)
Proof. Since G ±,α = 0 and G α = 0 for components α = ρ and α = e, we only deal with the impulsion equation. Using the midpoint approximation formula for the integral, we have

C i = φ * i+1/2 φ i G u (W(φ; V i )) dφ + G -,u (V i , V i+1 ) - φ * i-1/2 φ i G u (W(φ; V i )) dφ + G +,u (V i-1 , V i ) = φ * i+1/2 φ i P (φ) dφ -P i φ i+1 -φ i 2 - φ * i-1/2 φ i P (φ) dφ + P i φ i -φ i-1 2 = φ * i+1/2 φ * i-1/2 P (φ) dφ -P i φ i+1 -φ i-1 2 = P (φ i )(φ * i+1/2 -φ * i-1/2 ) -P i φ i+1 -φ i-1 2 + O(∆φ 2 ) = O(∆φ 2 ).
Hence we deduce estimate (24) . Combining the two lemmas, we have the following theorem for the regular steady-state solutions.

Theorem 1. The scheme H using the modified Rusanov conservative flux (18) and the associated Rusanov nonconservative flux ( 20) is a first-order scheme following definition 1.

Nonconservative Rusanov scheme order : the discontinuous case

We now turn to the case of a discontinuous steady-state solution where φ(x) jumps from φ L to φ R at x = 0. We define the steady-state solution using the φ parameterization U (φ) = W(φ; V L ) of the integral curve of equation ( 6). We denote by K 0 = [-∆x, 0] and K 1 = [0, ∆x] the cells on both side of the discontinuity and we study the order of

H(V L , V L , V R ) on K 0 and H(V L , V R , V R ) on K 1 . In the sequel, we set φ * = φ L +φ R 2 , U * = W(φ * ; V L ) and V * = V (φ * , U * ).
Proposition 4. Assume that the velocity is null then the scheme is exactly well-balanced:

H(V L , V L , V R ) = H(V L , V R , V R ) = 0. ( 25
)
Proof. If u L = 0, we deduce that u R = 0 thanks to the parameterization (8) and conclude that ρ L = ρ R and E L = E R . The viscosity contribution vanishes and we have

H(V L , V L , V R ) = F(V L , V R ) + G -(V L , V R ) -F(V L , V L ) = F (V L ) + F (V R ) 2 -F (V L ) + G -(V L , V R ) = 1 2   0 P R φ R -P L φ L 0   -    0 φ R -φ L 2 P L 0    = 0.
Then the scheme is exactly well-balanced.

Remark 3. If the velocity does not vanish, we loose the first-order due to the viscosity contribution. On the other hand, the nonconservative flux writes

φ * φ L P (φ) dφ -P L φ R -φ L 2 = O(∆φ)
which is not enough to provide a first-order scheme in the sense of definition 1. A solution consists in using a more precise numerical integration formula like the trapezoidal formula but the intermediate value P * has to be determined.

Second-order nonconservative scheme

Second-order schemes provide more precise approximation of the solution and less viscosity in the vicinity of shocks. We propose a second-order method based on the MUSCL reconstruction to provide better accuracy of the numerical approximation. Our method derives from the following remark, if φ is regular (at least differentiable) then ∂ x φ is a volumic source term. On the contrary, if φ presents a discontinuity at point x 0 , then ∂ x φ is a dirac distribution supported by the point x 0 , hence a punctual source term. The second-order method we proposed is based on the following porosity decomposition

φ(x) = φ d (x) + φ r (x)
where φ d and φ r are respectively the discontinuous part and the regular part (at least C 1 ). From a numerical point of view, the regular contribution is approximated by a volumic term while the discontinuous part is taken into account with the nonconservative flux. If we consider the first-order scheme, φ is a constant piecewise function and φ r = 0. The nonconservative contribution is exclusively evaluated with the nonconservative numerical flux. We now turn to the second-order situation where we use the classical MUSCL technique.

The MUSCL reconstruction

For each interface at point x i+1/2 , we aim to compute new values V - i+1/2 and V + i+1/2 on both side of the interface. We then substitute the first-order numerical flux F(V i , V i+1 ) with the second-order numerical flux F(V - i+1/2 , V + i+1/2 ) to compute the conservative flux across the interface x i+1/2 . We proceed in the same way with the nonconservative flux.

For any approximation (α i ) i∈Z with α = φ, ρ, u, P , we define the slopes

σ α i+1/2 = α i+1 -α i ∆x
which represent an approximation of ∂ x α on the interfaces x i+1/2 , i ∈ Z. To evaluate an approximation σ α i of the derivative on cell K i , we set

σ α i = ψ(σ α i-1/2 , σ α i+1/2 ) (26)
where the function ψ is a limiter function. For example, we use the minmod limiter

ψ(p, q) = minmod(p, q) = 0 if pq ≤ 0, p |p| min(|p|, |q|) if pq > 0. ( 27
)
We then propose the decomposition of φ = φ d (x) + φ r (x) such that φ d (x) is constant piecewise, φ r (x) is continuous linear piecewise where ∂ x φ r (x) = σ φ i+1/2 on each cell K i . We define the new values on the interface x i+1/2 setting for α = φ, ρ, u, P

α - i+1/2 = α i + σ i-1 ∆x 2 , α + i+1/2 = α i+1 -σ i ∆x 2 .

Second-order nonconservative scheme

Using the generic finite volume scheme proposed in (10), we consider the second-order scheme

U n+1 i = U n i - ∆t ∆x F(V n,- i+1/2 , V n,+ i+1/2 ) + G -(V n,- i+1/2 , V n,+ i+1/2 ) - F(V n,- i-1/2 , V n,+ i-1/2 ) + G + (V n,- i-1/2 , V n,+ i-1/2 ) +S(V n,- i+1/2 , V n,+ i-1/2 ) ( 28
)
with

S(V n,- i+1/2 , V n,+ i-1/2 ) =   0 P n i σ φ i ∆x 0   =     0 P + i-1/2 + P - i+1/2 2 (φ - i+1/2 -φ + i-1/2 ) 0     . (29)
The term S corresponds to the contribution of the regular part of P ∂ x φ.

Remark 4. If φ is a constant piecewise function, we does not recover exactly the first-order scheme since the MUSCL reconstruction provides new values for ρ, u and P at the interfaces. On the other hand, if φ is regular, the φ jump at the interface φ + i+1/2 -φ - i+1/2 is very small, of order ∆x 2 , hence the principal nonconservative contribution is furnished by the source term S.

The Riemann problem configurations

In this section, we present a complete description of the admissible configurations, i.e. the succession of simple waves separated by constant states, for the nonconservative Euler system. Two major situations arise whether all the waves are distinct or two waves are superposed: the resonant cases. The resonant cases have been analysed by [START_REF] Noussair | Analysis of nonlinear resonance in conservation law with point sources and well-balanced scheme[END_REF] in the framework of a scalar nonconservative problem with a source term. The works of [START_REF] Noussair | Riemann problem with nonlinear resonance effects and wellbalanced godunov scheme for shallow fluid past an obstacle[END_REF] and [START_REF] Chinnayya | A well-balanced numerical scheme for shallow-water equations with topography: resonance phenomena[END_REF] deal with the resonant configurations for the shallow-water problem while Le Floch & Thanh (2003) treat the resonant case for the isentropic nonconservative Euler problem. The configurations in the resonant case for a general nonconservative hyperbolic system are presented in [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF] where a sharp study of the interactions between a genuinely nonlinear characteristic field with the linearly degenerated field associated to the nonconservative term is performed (see also the pioneer works of [START_REF] Isaacson | Nonlinear resonance in systems of conservation laws[END_REF]).

We do not intend to solve the Riemann problem leading to huge complex nonlinear problems as [START_REF] Andrianov | On the solution to the Riemann problem for the compressible duct flow[END_REF]; [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF]; Le [START_REF] Le Floch | The riemann problem for fluid flows in a nozzle with discontinuous cross-section[END_REF]; [START_REF] Noussair | Riemann problem with nonlinear resonance effects and wellbalanced godunov scheme for shallow fluid past an obstacle[END_REF]; [START_REF] Chinnayya | A well-balanced numerical scheme for shallow-water equations with topography: resonance phenomena[END_REF] but we propose here a construction of all the configurations based on three criterions: the configuration stability criterion, the sign criterion and the Mach criterion we shall define in the sequel. A study of the eigenvalue sign and the criterions allow to disqualify the noncompatible configurations. Moreover, a specific study of the linearly degenerated field associated to the nonconservative term brings new informations to construct the configurations for the resonant case.

Eigenvalues and eigenvectors

To study the nonconservative Euler system which describes gas flow in porous media, we add equation ∂ t φ = 0 to system (2), named the augmented system (see [START_REF] Andrianov | On the solution to the Riemann problem for the compressible duct flow[END_REF]; [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF] for justifications) leading to a problem with four unknowns. For regular solution V , the augmented system equipped with the perfect gas law P = (γ -1)ρe with γ ∈]1, 3[ can be rewritten in the form

∂ t V + A(V )∂ x V = 0, with A(V ) =      0 0 0 0 ρu φ u ρ 0 0 0 u 1 ρ γP u φ 0 γP u      .
We get a non-strictly hyperbolic system with the eigenvalues

λ 0 = 0, λ 1 = u -c, λ 2 = u, λ 3 = u + c,
and the associated eigenvectors

R 0 =     φ(u 2 -c 2 ) -ρu 2 uc 2 -ρu 2 c 2     , R 1 =     0 -ρ c -ρc 2     , R 2 =     0 1 0 0     , R 3 =     0 ρ c ρc 2     .
The characteristic fields (or simple waves) associated to eigenvalues λ 1 and λ 3 are genuinely nonlinear while the characteristic fields associated to λ 0 and λ 2 are linearly degenerated. We say that V is a subsonic, sonic or supersonic state if we have u 2 < c 2 , u 2 = c 2 or u 2 > c 2 respectively.

Eigenvalue λ 0 = 0 characterizes the change of porosity and a new difficulty arises since we are no longer able to order the four eigenvalues. In particular, for sonic state V , we have only three distinct eigenvalues (λ 0 = λ 1 or λ 0 = λ 3 ).

Eigenvectors R 1 , R 2 and R 3 correspond to the classical Euler system eigenvectors with constant porosity. It results that the simple waves parameterization (shock, rarefaction, contact discontinuity) associated to each eigenvalue is identical to the classical Euler system ones.

The stability configuration criterion

In the sequel, we denote by k -w the simple wave associated to the eigenvalue λ k . For k = 1 or k = 3, k -r represents a rarefaction whereas k -s represents a shock. For k = 0 and k = 2 the simple waves are contact discontinuities we still denote k -w. For given left and right states V L and V R , we seek an autosimilar solution for the Riemann problem and we shall describe using the following definition.

Definition 2 (configuration). A configuration is the succession of simple waves from left to right which separate constant states.

As an example, the configuration C = {1 -r; 0 -w; 2 -w; 3 -s} means that we have four simple waves separating five constant states : V L is linked to the state V l by the 1-r rarefaction, V l is linked to V r by the contact discontinuity 0 -w while V r is linked to V a by the contact discontinuity 2 -w. At last V a is linked to V R by the 3 -s shock (see Fig. 1).

We now introduce a criterion to select the configurations which are of interest from a numerical point of view. Each simple wave k -w corresponds to a one-parameter characteristic curve in the phase space. Existence of a solution for the Riemann problem means that there exists a path in the phase space from V L to V R following the characteristic curves given by the configuration. Since V is a vector of R 4 , we need four independent characteristic curves to link V L to V R .

Definition 3 (configuration stability criterion). Assume that for given left and right states V L , V R we have a solution for the Riemann problem with a given configuration X. We say that the configuration is stable if there exists a small enough ε > 0 such that : for any physical states V R and

1-r 0-w 2-w 3-s t x V L V l V r V a V R
V L with | V L -V L | < ε and | V R -V R | < ε
, there exists a solution for the Riemann problem associated to the initial states V L and V R with the same configuration X.

For example, let us consider the Riemann problem with V L = (φ L = 1, ρ L = 1, u L = 10, P L = 10 5 ) and V R = (φ R = 1, ρ R = 5, u R = 10, P R = 10 5 ), the associated configuration is constituted of a unique contact discontinuity X = {2 -w}. If one slightly modifies the pressure, the velocity or the porosity, new simple waves appear and configuration X no longer holds. The goal of the stability criterion is to only focus on the configurations which are important from a numerical point of view.

The sign criterion

Let us denote by P -= {(x, t); x < 0, t > 0} and P + = {(x, t); x > 0, t > 0} the two half-planes. As noticed [START_REF] Andrianov | On the solution to the Riemann problem for the compressible duct flow[END_REF], the augmented system reduces to the classical Euler system with constant porosity in each half-plane. In particular, there exists at most three ordered simple waves in each half-plane leading to an important number of potential configurations. Nevertheless, the admissible configuration number will be drastically cut down thanks to the following sign criterion.

Definition 4 (sign criterion). The sign velocity does not change across the interface x = 0.

We apply the Rankine-Hugoniot condition to the conservative mass conservation equation of system (2) for the contact discontinuity associated to λ 0 .

Since the shock velocity is null, we get [φρu] = 0 at interface x = 0, hence the sign of u is preserved for positive density and porosity.

Proposition 5. The following successions of waves and constant states

(a) {3 -w, 0 -w, 1 -w}, (b) {3 -w, 0 -w, 2 -w}, (c) {2 -w, 0 -w, 1 -w}, (d) {2 -w, 0 -w, 2 -w}, are not admissible.
Proof. We present the proof for the case (a). Let denote by V l and V r the constant states situated on the left and right of the interface x = 0. Since we deal with rarefactions or entropic shocks, the Lax entropy condition yields

λ 3 (V l ) < 0, λ 1 (V r ) > 0.
It results that u l < 0 while u r > 0 which is a contradiction with the sign criterion, thus the succession (a) is not admissible. The other successions are treated in the same way.

We list in Fig. 2 the configurations which satisfy both the sign and the stability configuration criterions. Configurations A, B, C and D involve the four simple waves only one time (see [START_REF] Andrianov | On the solution to the Riemann problem for the compressible duct flow[END_REF]). The two last configurations F and G correspond to the situation when the 1 -w or the 3 -w waves are splitted by the interface. We have not presented here the resonant situations when the 1 -w or the 3 -w wave is superposed with the 0 -w wave. To characterize more precisely the two last situations, we study in the next subsection the transition across the interface when the porosity changes.

The MACH criterion

Across the porous interface x = 0, the state V l changes to the state V r in function of the porosity variation. Since mass and energy equations in system (2) are written in a conservative form, we deduce the Rankine-Hugoniot condition across the contact discontinuity λ 0 :

[φρu] = 0, [u(φE + φP )] = 0.

(30)

We can not use the classical Rankine-Hugoniot condition for the impulsion equation since there is no conservative form but since we deal with a contact discontinuity, we use the Riemann invariants associated to the simple wave 0 -w (see [START_REF] Andrianov | On the solution to the Riemann problem for the compressible duct flow[END_REF], p. 881):

1-w 0-w 2-w 3-w 1-w 0-w 2-w 3-w 1-w 0-w 2-w 3-w 1-w 0-w 1-w 0-w 2-w 3-w 3-w 2-w 3-w 0-w 1-w 2-w V L V l V r V a V R V L V L V R V R V R V R V R V L V L V L V a V a V r V r V a V b V l V a V a V b V l V r V b V b V l V r V l configuration A configuration B configuration C configuration D configuration G configuration F
w 1 0 (V ) = φρu, w 2 0 (V ) = P ρ γ , w 3 0 (V ) = u 2 + 2γP (γ -1)ρ . ( 31
)
We note that w 1 0 and w 3 0 correspond to relations (30) while the entropy invariant w 2 0 comes from the impulsion equation. Moreover, the three invariants provide the same equations which characterize the steady-state solutions (8) (see [START_REF] Chinnayya | A well-balanced numerical scheme for shallow-water equations with topography: resonance phenomena[END_REF] p.12 for the shallow-water problem).

For a given left state V l of porosity φ L , we seek a right state of porosity φ R which preserves the Riemann invariants (31). To this end, we denote by

D = φ L ρ l u l , S = P l ρ γ l , H = u 2 l + 2γP l (γ -1)ρ l , (32)
and we seek ρ r , u r and P r such that w 1 0 (V r ) = D, w 2 0 (V r ) = S, w 3 0 (V r ) = H. After some algebric manipulations, the problem consists in finding the density as a solution of the scalar equation

g(ρ; S, H) = - D φ R 2 ( 33
)
where we have defined

g(ρ; S, H) := 2γS γ -1 ρ γ+1 -ρ 2 H. (34)
If we manage to compute ρ, we obtain the pressure P and the velocity u which preserve the entropy S and the flux D. We now deal with the calculation of ρ. Differentiation of function g gives

g ′ (ρ) = 2ρ γ(γ + 1) γ -1 Sρ γ-1 -H .
The derivative only vanishes at point ρ son = ρ son (H, S) given by γρ γ-1 son = (γ -1)H (γ + 1)S . Since γ ∈]1, 3[, function g ′ is positive for ρ > ρ son and negative for ρ < ρ son . To sum up, we have the proposition.

Proposition 6. Function g admits a minimum at point ρ son which only depends on S and H characterized by

γρ γ-1 son = (γ -1)H (γ + 1)S and g(ρ son ) = 1 -γ 1 + γ Hρ 2 son < 0.
g

• if ρ > ρ son , then u 2 < c 2 (subsonic branch),
• if ρ < ρ son , then u 2 > c 2 (supersonic branch). For H and S given, function g reaches the minimum at the sonic point ρ son . Consequently, for a given D, we deduce the minimum value φ min = φ min (D, S, H) such that equation (34) has a solution given by

γ -1 γ + 1 Hρ 2 son = D φ min 2 .
Hence we get

φ 2 min = (γ + 1)D 2 (γ -1)Hρ 2 son . ( 35
)
The solutions of equation ( 34) are detailed in the following corollary.

29

Corollary 1. Let (φ, D, S, H) be fixed and ρ son be the density associated to S and H. We denote by

χ(φ, D, S, H) = D φ 2 + g(ρ son ) = D φ 2 + 1 -γ 1 + γ Hρ 2 son .
We have the three situations:

1. If χ(φ, D, S, H) < 0 (i.e. φ > φ min ), there are two solutions ρ sup < ρ son and ρ sub > ρ son for equation ( 33). 2. If χ(φ, D, S, H) > 0 (i.e. φ < φ min ), there is no solution for equation ( 33). 3. If χ(φ, D, S, H) = 0 (i.e. φ = φ min ), there is the unique solution ρ son for equation ( 33).

In case (1) two solutions are available and the question of the choice is crucial.

For example, assume that we have a prescribed left subsonic state V l = (φ L , ρ l , u l , P l ). For φ R in the vicinity of φ L , we have two possible solutions but only one belongs to the same subsonic branch ρ > ρ son . For continuity reason, it is judicious that the solution ρ of equation g

(ρ) = -D φ R 2
belongs to the subsonic branch. Chinnayya et al. (2004) (p. 13) give a clever justification for the shallow-water problem using a regularisation of the topography while Goatin & Le Floch (2004) (p. 891) state a similar criterion.

Definition 5 (MACH criterion). The subsonic or supersonic regimes are preserved across the interface. If V l and V r are the states on the left and right sides of the contact discontinuity 0 -w then they are both subsonic or supersonic.

Note that the MACH criterion does not apply in case of a sonic state. For example, if V l is sonic, V r could be subsonic or supersonic. A direct consequence of the MACH criterion concerns the splitting of the 1 -w or 3 -w waves.

Proposition 7. The simple waves 1 -w or 3 -w can not be splitted with two constant states on both sides of the interface x = 0.

Proof. We give the proof when the 1 -w wave is splitted into two parts leading to the configuration F = {1-w, 0-w, 1-w, 2-w, 3-w} represented

1-w 0-w 2-w 3-w t x V L V a V b V l V R V r 1-w constant states
Figure 4: The 1w is splitted with two constant states on both sides of the interface. Such a situation is not available since the MACH and sign criterions are not respected.

in Fig. 4. The 1 -w wave is splitted by the 0 -w wave with two constant states on both sides of the interface with V l on the left and V r on the right. Due to the presence of the 1 -w in the half-plane P + , state V r is supersonic with u r > c r whether the simple wave is a rarefaction or an entropic shock thanks to the Lax condition. In consequence, the sign and MACH criterions yield that u l > c l . We now deal with the 1 -w which separates the two constant states V L and V l in P -. Since V l satisfies u l > c l , the 1 -w can not be a rarefaction because V l belongs to the P -half-plane. On the other hand, the Lax condition says u L -c L > σ 1 > u l -c l > 0 hence the 1 -s shock velocity σ 1 is positive which is a contradiction with the assumption that the shock belongs to P -.

The splitting wave configurations LR and RR

Let us consider the situation where the 1 -w simple wave is splitted into two waves by the interface x = 0. Proposition 7 says that there is at most a constant state on one side of the interface. Consequently, the 1 -w wave situated on the other side is a rarefaction which touches the interface at x = 0. We present in Fig. 5 the two situations for the 1 -w wave whether the rarefaction is on the left or on the right leading to the configuration LR 1 (Left Rarefaction for the 1 -w wave) and RR 1 (Right Rarefaction for the 1 -w wave). Similar configurations denoted by LR3 and RR 3 hold for the 31 3 -w wave splitting. Such configurations have been proposed by [START_REF] Noussair | Analysis of nonlinear resonance in conservation law with point sources and well-balanced scheme[END_REF] (p. 324, case 6 and 8) for the scalar nonconservative problem, by [START_REF] Noussair | Riemann problem with nonlinear resonance effects and wellbalanced godunov scheme for shallow fluid past an obstacle[END_REF] [START_REF] Lowe | Two-phase shock-tube porblems and numerical methods of solution[END_REF] for the two-phase flow problem, by [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF] (p. 897, case 2a(A) and case 2b(A')) for a general nonconservative hyperbolic system.

1-w 0-w 2-w 3-w t x V L V a V b V l V R 1-r 1-r 0-w 2-w 3-w t x V L V a V b V R V r 1-w V r V l Configuration RR 1 Configuration LR 1 Figure 5:
The 1w is splitted into a rarefaction on one side and a constant state on the other side which provides the following configurations whether the rarefaction is on the left (LR 1 configuration) or on the right (RR 1 configuration) of the interface.

We now analyse such a situation where we shall prove the important result: the rarefaction takes place only on the lower porosity side.

Proposition 8. The 0 -w wave parts the 1 -w or 3 -w simple waves in the following way (see Fig. 6 for the configuration designations).

Rarefaction from P -. We assume there is a rarefaction in the P - half-plane which reaches the interface x = 0 to a limit state V l such that u l = ±c l and jumps to a constant state V r on the right:

• If φ L > φ R , there is no solution; • If φ L < φ R ,
there is two possibilities: * a 1 -r rarefaction with u l , u r positive and V r supersonic (case LR 1 ), * a 3-r rarefaction with u l , u r negative and V r subsonic (case LR 3 ).

Rarefaction from P + . We assume there is a rarefaction in the P + half-plane which reaches the interface x = 0 to a limit state V r such that u r = ±c r and jumps to a constant state V l on the left:

• If φ L < φ R ,
there is no solution;

• If φ L > φ R , there is two possibilities: * a 1-r rarefaction with u l , u r positive and V l subsonic (case RR 1 ), * a 3 -r rarefaction with u l , u r negative and V l supersonic (case RR 3 ).

1-r 0-w 2-w 3-w t x 2-w 3-w t x 1-w 0-w t 3-r 2-w 3-w t x 1-w t x 1-w 0-w 2-w 3-w x x x φ L <φ R 3-r t x 1-r 0-w 2-w 3-w x 1-w LR 1 RR 1 LR 3 RR 3 V L V L V L V L V a V r V a V l V a V b V a V b V b V r V b V l V R V R V R V R V r V l V l V r φ L >φ R
Figure 6: Admissible configurations where the 1-w or the 3-w simple waves are splitted by the 0-w contact discontinuity. Configurations LR 1 and RR 1 correspond to the 1-w wave with φ L < φ R and φ L > φ R respectively while configurations LR 3 and RR 3 correspond to the 3w wave with φ L < φ R and φ L > φ R respectively.

Proof. We only deal with the left rarefaction. We reach the state V l with a 1 -r or a 3 -r rarefaction from the left such that u 2 l = c 2 l . For the limit state, we associate the values of D, S and H given by relations (32) and we seek a solution ρ of the equation

g(ρ) + D φ R 2 = 0. ( 36
) • If φ R < φ L then D φ R 2 > D φ L 2
. Since the left state is sonic, function g attains its global minimum at point ρ l and we have g(ρ)

+ D φ L 2 ≥ g(ρ l ) + D φ L 2 = 0 for all ρ > 0. It results that g(ρ) + D φ R 2
> 0 and equation ( 36) has no solution.

• Assume now that φ R > φ L , we then have

D φ R 2 < D φ L 2 . Since g(ρ l ) + D φ L 2 = 0, we deduce that g(ρ l ) + D φ R 2 < 0.
From proposition 6 there exists two solutions ρ sub > ρ son = ρ l and ρ sup < ρ son = ρ l for equation (36). We now distinguish the situations whether we have a 1 -r or a 3 -r rarefaction.

* For a 1 -r rarefaction (configuration LR 1 ). We have, on one hand, a 1 -r rarefaction in the half-plane P -from the V L state to a V l state located at x = 0 where V l is prescribed with the condition u l = c l and we have, on the other hand, a 1 -w simple wave in the half-plane P + . We cross the interface x = 0 using the parameter φ R to provide state V r on the right side of the interface. Since ρ l is the minimum of function g and φ R > φ L , there is two solutions: ρ sub for the subsonic branch and ρ sup for the supersonic branch. Since there exists a 1 -w wave in the half-plane P + , then V r must be supersonic whether the wave is a rarefaction or an entropic shock. Hence, we have to choose ρ = ρ sup . * For a 3 -r rarefaction (configuration LR 3 ). We have, on one hand, a 3 -r rarefaction in the half-plane P -from the V b state to a V l state located at x = 0 and we have, on the other hand, a 3 -w simple wave in the half-plane P + . As in the previous case, we have two possible solutions for the state V r = V R situated on the right side of the interface. Since there exists a 3 -w wave in the half-plane P + , then V R must be subsonic and we have to choose ρ = ρ sub . To sum up, we propose the following criterion when a rarefaction touches the porous interface.

Definition 6 (rarefaction in porous media criterion). A rarefaction can only reach the interface from the lower porosity side.

Remark 5. To our knowledge, such a criterion has not been clearly brought to the fore before. For the shallow-water problem, Chinnayya et al. (2004) (p. 19) show that the transition with a rarefaction reaching the interface only occurs on the lower vertical height side. Moreover, in the general case of nonconservative hyperbolic systems studied by [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF], all the rarefactions reaching the interface satisfy the criterion in definition 6. Remark 6. When a k -w wave splitting occurs, we have shown that the first part of the wave situated on the lower porous side must be a rarefaction. The second part of the wave (in the higher porosity side) is not a priori a rarefaction and could also be a shock. We shall present examples of various configurations in the numerical test section where the second part of the k -w wave is a rarefaction or a shock.

The resonant configurations

Since the eigenvalues are not strictly ordered due to the presence of the eigenvalue λ 0 = 0, the question arises when a k -s shock crosses the interface leading to the so-called resonant situation. We have to treat two situations whether the k -w wave is splitted (with a rarefaction on one side) or not.

The simpler situation concerns a simple k -s shock which crosses the interface. The k -s shock splits the 0 -w contact discontinuity into two parts where an intermediate porosity φ s is introduced (see [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF] p.892 case 1a (C) for example). We denote by R 1 the resonant configuration with the 1 -s shock while R 3 represents the resonant configuration with the 3 -s shock. A more complex situation arises when the k -w is splitted into two parts: a rarefaction on the lower porosity side and a stationary shock of null velocity (see [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF] p. 897 case 2a (C) and [START_REF] Chinnayya | A well-balanced numerical scheme for shallow-water equations with topography: resonance phenomena[END_REF] p. 22 section 3.3.4 for examples). The shock itself splits the contact discontinuity and we shall one more time introduce an intermediate porosity. We denote by LRR 1 and LRR 3 the Left Rarefaction and Resonant configurations associated to the 1 -w and 3 -w respectively when φ L < φ R while RRR 1 and RRR 3 are the Right Rarefaction and Resonant configurations when φ L > φ R .

Configurations R 1 and R 3

The resonant configuration R 1 appears when the 1-s shock and the 0-w are superposed. To illustrate the phenomena, let us consider the situation where we have a 1 -s shock on the left of the interface (configuration C in Fig. 7). We increase the left velocity u L until the shock reaches the interface x = 0. At that very moment, the shock touches the left side of the contact discontinuity. If we increase a little bit the velocity, the shock shares the interface in two parts introducing an intermediate porosity φ s ∈ [φ L , φ R ] leading to the following configuration {0 -w, 1 -s, 0 -w} at the same point x = 0 (R 1 configuration). We increase one more time the velocity until the 1-s shock reaches the right side of the interface, i.e. φ s = φ R . Configuration R 1 respects the stability configuration criterion since the 1 -s shock stays on the interface for small perturbations of V L and V R but the intermediate porosity φ s changes. The same situation arises with a 3 -w shock leading to the R 3 resonant configuration. Such configurations have been proposed by [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF] (p. 896, case 2b(C')) for a general nonconservative hyperbolic system. We now detail the construction of configuration R 1 composed of two 0-w contact discontinuities shared by a stationary 1 -s shock. We assume that the physical states are V l and V r on the left and right of the interface. To construct the transition, we introduce an intermediate porosity φ s and we consider the following waves succession linking V l to V r (see Fig. (8)).

1-s 0-w 2-w 3-w 0-w 2-w 3-w 2-w 0-w 1-s V L V L V R V R V R V L V a V a V a V r V l V l V r V l V r configuration C configuration D Configuration R 1 0-w 3-w 1-s V b
-From V l to V s,l : we have a 0 -w contact discontinuity from φ L to an intermediate porosity

φ s ∈ [φ L , φ R ].
-From V s,l to V s,r : we have a stationary 1 -s shock of velocity σ 1 = 0.

-From V s,r to V r : we have a 0 -w contact discontinuity from φ s to φ R . The porosity φ s has to be chosen such that the intermediate states V s,l and V s,r define a stationary shock. Since we require an 1 -s entropy shock, the Lax condition implies u s,l -c s,l > 0 > u s,r -c s,r hence V s,l is supersonic while state V s,r is subsonic. Since V s,l is a supersonic state, the MACH criterion says that the left state V l must be a supersonic state while V s,r and V r have to be subsonic states.

φ R V r V s,l V s,r V l φ s φ L x=0 { 0-w 0-w 1-s
To compute the two states and the intermediate porosity, we proceed in the following way. For a given V l , we define the unique ρ s,l = ρ s,l (φ s ) on the supersonic branch of g (see Fig. 3) such that we satisfy the relation

g(ρ s,l ; S l , H l ) = 2γS l γ -1 ρ γ+1 s,l -ρ 2 s,l H l = - D l φ s 2 , ( 37
)
where H l , S l and D l are computed with the left state V l . In the same way for a given V r , we define a unique ρ s,r = ρ s,r (φ s ) on the subsonic branch of g such that we satisfy the relation

g(ρ s,r ; S r , H r ) = 2γS r γ -1 ρ γ+1 s,r -ρ 2 s,r H r = - D r φ s 2 , ( 38
)
where H r , S r and D r are computed with the right state V r . With ρ s,l and ρ s,r in hand, we compute the velocities u s,l , u s,r and the pressures P s,l , P s,r and we fix the intermediate porosity φ s using the stationary shock condition ρ s,l u s,l = ρ s,r u s,r .

Of course φ s is implicitely given by relation ( 39) and an iterative algorithm should be employed to compute an approximation of φ s such that we satisfy the three relations ( 37), ( 38) and ( 39).

Configurations LRR 1 and LRR 3

We now deal with a more complex situation when a k -w wave is splitted into a k -r rarefaction on the left side and a k -s shock superposed with the 0 -w contact discontinuity. For such a situation, we must have φ L < φ R such that the rarefaction occurs on the left side. We then obtain a Left Rarefaction and Resonant configuration LRR 1 or LRR 3 for the 1 -w or the 3 -w respectively. Configuration LRR has been studied by Chinnayya et al. (2004) (p. 19, case c) for the shallow-water problem, by [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF] (p. 897, case 2a(C)) for a general nonconservative hyperbolic system. We detail the LRR 1 configuration as an example (see Fig. 9 for the notation).

φ R V r V s,l V s,r φ s φ L V l V L x=0 { 0-w 0-w 1-s 1-r
-A 1 -r rarefaction takes place on the left of the interface from V L to the sonic state V l -From V l to V s,l : we have a 0 -w contact discontinuity from φ L to an intermediate stage

φ s ∈ [φ L , φ R ] such that V s,l is supersonic.
-From V s,l to V s,r : we have a stationary 1 -s shock of velocity σ 1 = 0 and V s,r is subsonic.

-From V s,r to V r : we have a 0 -w contact discontinuity from φ s to φ R .

Like the R 1 configuration, one has to choose the intermediate porosity such that V s,l and V s,r are linked with a stationary shock of null velocity and respect the entropy Lax condition u 1,l -c 1,l > 0 > u 1,r -c 1,r . To compute the two states and the intermediate porosity, we have to fix φ s such that relations (37), ( 38) and ( 39) are satisfied. The only difference with the configuration R 1 is that V l is the sonic state which links V L with a 1 -r rarefaction.

Configurations RRR 1 and RRR 3

We now analyse a similar situation when the k -w wave is splitted into a k -r rarefaction on the right side and a k -s shock superposed with the 0 -w contact discontinuity. For such a situation, we must have φ L > φ R such that the rarefaction occurs on the right side. We then obtain a Right Rarefaction and Resonant configuration RRR 1 or RRR 3 for the 1 -w or the 3 -w respectively.

We detail the RRR 1 configuration as an example (see Fig. 10 for the notation).

-From V l = V L to V s,l : we have a 0 -w contact discontinuity from φ L to an intermediate stage

φ s ∈ [φ L , φ R ].
-From V s,l to V s,r : we have a stationary 1 -s shock of velocity σ 1 = 0.

-From V s,r to V r : we have a 0 -w contact discontinuity from φ s to φ R such that V r is a sonic state.

-A 1 -r rarefaction takes place on the right of the interface from V r to V a . Like the LRR 1 configuration, one has to choose the intermediate porosity such that V s,l and V s,r are linked with a 1 -s stationary shock of null velocity and respect the entropy Lax condition u 1,l -c 1,l > 0 > u 1,r -c 1,r . The only difference with the configuration LRR 1 is that the rarefaction is on the right and φ s has to be chosen such that V r is a sonic state. This last point leads to a more complex problem from the numerical point of view when we solve the inverse Riemann problem.

V s,l V s,r V L V r V a φ s V l φ L φ R x=0 { 0-w 0-w 1-s 1-r

Parameterization of the 1 -w wave

We consider two states V 1,l and V 1,r separated by the simple wave 1 -w which is a shock (P 1,r > P 1,l ) or a rarefaction ( P 1,r < P 1,l ). For a given left state V 1,l and a prescribed pressure of the right state P 1,r = P 1 , we deduce the other components of V 1,r with the following relation.

-If P 1,r > P 1,l , we have a shock and we compute

ρ 1,r = ρ 1,l P 1,l (γ -1) + P 1,r (γ + 1) P 1,l (γ + 1) + P 1,r (γ -1) , u 1,r = u 1,l - (P 1,r -P 1,l )(ρ 1,r -ρ 1,l ) ρ 1,r ρ 1,l .
The discontinuity moves with the velocity

σ 1 = ρ 1,l u 1,l -ρ 1,r u 1,r ρ 1,l -ρ 1,r .
-If P 1,r < P 1,l , we have a rarefaction and we compute

ρ 1,r = ρ 1,l P 1,r P 1,l 1 γ , c 1,r = γP 1,r ρ 1,r , u 1,r = u 1,l + 2 γ -1 (c 1,l -c 1,r )
with the left state sound velocity c 1,l = γP 1,l ρ 1,l . The rarefaction area is characterized by

u 1,l -c 1,l ≤ x t ≤ u 1,r -c 1,r .
Note that the porosity does not change hence φ 1,l = φ 1,r .

Parameterization of the 3 -w wave

We consider two states V 3,l and V 3,r separated by the simple wave 3 -w which is a shock (P 3,r < P 3,l ) or a rarefaction ( P 3,r > P 3,l ). For the given left state V 3,l and a prescribed pressure of the right state P 3,r = P 3 , we deduce the other components with the following algorithm.

-If P 3,r < P 3,l , we have a shock and we compute ρ 3,r = ρ 3,l P 3,l (γ -1) + P 3,r (γ + 1)

P 3,l (γ + 1) + P 3,r (γ -1) , u 3,r = u 3,l - (P 3,r -P 3,l )(ρ 3,r -ρ 3,l ) ρ 3,r ρ 3,l .
The discontinuity moves with the velocity σ = ρ 3,r u 3,r -ρ 3,l u 3,l ρ 3,r -ρ 3,l .

-If P 3,r > P 3,l , we have a rarefaction and we compute

ρ 3,r = ρ 3,l P 3,r P 3,l 1 γ , c 3,r = γP 3,r ρ 3,r , u 3,r = u 3,l - 2 γ -1 (c 3,l -c 3,r )
with the left state sound velocity c 3,l = γP 3,l ρ 3,l . The rarefaction area is given by

u 3,l + c 3,l ≤ x t ≤ u 3,r + c 3,r .
Note that the porosity does not change hence φ 3,l = φ 3,r .

Parameterization of the 2 -w wave

For a given left state V 2,l , we have u 2,r = u 2,l , P 2,r = P 2,l since pressure and velocity are unchanged across the contact discontinuity moving with the velocity σ = u 2,l = u 2,r . We prescribe the right density ρ 2,r = ρ 2 and the sound velocity for the right state is given by c 2,r = γP 2,r ρ 2,r .

6.1.4. Parameterization of the 0 -w wave For a given left state V 0,l , we compute

D = φ L ρ 0,l u 0,l , S = P 0,l ρ γ 0,l , H = u 2 0,l + 2γP 0,l (γ -1)ρ 0,l .
We then seek the solutions of the nonlinear equation

2γS γ -1 ρ γ+1 -ρ 2 H = - D φ R 2
using for example an iterative Newton algorithm. Let γρ γ-1 son = (γ -1)H (γ + 1)S be the density corresponding to the sonic state, following corollary 1, we compute χ(φ R , D, S, H) and we have the three situations:

• If χ(φ R , D, S, H) > 0, there is no solution.

• If χ(φ R , D, S, H) = 0, there is the unique solution ρ 0,r = ρ son .

• If χ(φ R , D, S, H) < 0, there are two solutions and we choose ρ 0,r = ρ sup if V 0,l is a supersonic state and ρ 0,r = ρ sub if V 0,l is a subsonic state.

We then obtain the two other components with

u 0,r = D φ R ρ 0,r
, P 0,r = P 0,l ρ 0,r ρ 0,l γ .

Resolution of the inverse Riemann problem

Let V L be the left state and (P 1 , ρ 2 , P 3 , φ R ) be the set of prescribed intermediate values. We aim to construct all the intermediate states and V R which respect a given configuration C. We proceed with a "Try and Check" technique in the following manner. We start from the left side with V L , then we compute the first intermediate constant state V a situated just after the 1 -w. Then we check if the wave we have produced corresponds to the pattern of the configuration we are supposed to respect. If the wave is wrong (for example V a is supersonic whereas configuration C requires a subsonic state), we stop the procedure and no solution is available for such a configuration. If the 1 -w wave is admissible, we continue with the second wave and so on.

Configurations A, B, C, D

We only detail the algorithm for the configuration C = {1 -w, 2 -w, 0w, 3 -w}, the other cases are similar. We shall employ the notations introduced in Fig. 2. Let V L be the left state and (P 1 , ρ 2 , P 3 , φ R ) the set of intermediate values.

-If P L < P 1 we have a shock of velocity σ 1 otherwise we have a rarefaction. We then compute the new state V a . In case of a shock, we have to check that σ 1 < 0 whereas in case of a rarefaction we have to check that u a -c a < 0. If the condition is not satisfied, the C configuration is not available and we stop the algorithm.

-Furthermore the velocity u a has to be positive. Indeed, the velocity sign does not change across the interface x = 0 and the configuration C stipulates that λ 2 > 0. If u a < 0 then the C configuration is not available.

-Using φ R , we compute the intermediate state V b where we choose the subsonic solution using a Newton algorithm to solve equation (34). If there is no solution, the configuration C is not available.

-We compute the state V c after the contact discontinuity prescribing the density ρ 2 .

-If P 3 < P c we have a shock of velocity σ 3 whereas we have a rarefaction if P 3 > P c . We compute the state V d using the parameterization. In case of a shock, one has to check that σ 3 > u c .

The state V d corresponds to the right state V R we shall use in the simulations.

Configurations LR 1 and LR 3

The configurations LR are a little bit more complex. We only deal with the LR 1 case since the LR 3 configuration is computed in the similar manner. We have to construct a rarefaction which joins the left state V L to the sonic point. We shall employ the notations introduced in Fig. 6. The algorithm is the following.

-We first check that φ L < φ R and that u L -c L < 0. If one of these conditions is not satisfied, configuration LR 1 is not available.

-We compute the sonic state V l using the 1 -r rarefaction such that u l = c l .

-We compute the state V r solving relation (34). We take the solution ρ sup corresponding to the supersonic state.

-If P 1 < P r we have a rarefaction while we have a shock of velocity σ 1 if P 1 > P r and we compute an intermediate supersonic state V a with pressure P 1 . In case of a shock, we have to check that σ 1 > 0. If σ 1 < 0 the LR 1 configuration is not available.

-We compute the state V b using the contact discontinuity where we prescribe the density ρ 2 .

-We compute the state V c using the prescribed pressure P 3 . In case of a shock, we have to check that σ 3 > u b .

The state V c corresponds to the right state V R we shall use in the simulations.

Configurations LRR 1 and LRR 3

We now deal with the LRR configuration (φ L < φ R ) where a k -w wave is splitted into a rarefaction on the left side of the interface and a stationary shock. We shall employ the notations introduced in Fig. 9.

First of all, the V L state has to be subsonic since we have a 1-r rarefaction in the half-plane P -. From V L , we compute the sonic state V l linked to V L by the 1 -r rarefaction such that u l = c l . We proceed computing the stationary shock.

Computation of φ s and the V r state. Note that state V l is sonic and we have to compute a supersonic V s,l state and a subsonic V s,r state. To this end, for a given φ s , we compute the supersonic state V l,s linking V l when the porosity changes from φ L to φ s solving equation ( 37). With V l,s in hand, we compute V r,s using relation (40). We then deduce ρ s,r and u s,r , hence the state V s,r . Finally we compute V r solving equation ( 38) when the porosity changes from φ s to φ R . The intermediate porosity has to be fixed such that P r = P 1 where P 1 is the prescribed pressure. To this end, we introduce the function φ → f (φ) = P r -P 1 and we seek φ s such that f (φ s ) = 0. We use a Lagrange method where we initalise the algorithm with φ = φ L and φ = φ R .

Computation of the other states. We compute the state V a after the 2 -w contact discontinuity using the ρ 2 density. We finally determine state V R using the pressure P 3 . If P 3 > P a , we have a rarefaction while we have a shock if P 3 < P a . In the last case, we have to check that σ 3 > u a to obtain an admissible LRR 1 configuration.

Configurations RRR 1 and RRR 3

We treat the RRR configuration (φ L > φ R ) where a k -w wave is splitted into a rarefaction on the right side of the interface and a stationary shock. We shall employ the notations introduced in Fig. 10. The configuration is more complex to solve than the LRR configuration since we have to compute φ s such that V r is a sonic state.

Computation of φ s and the V r state. To obtain such a φ s we proceed as follows. Let φ s ∈ [φ L , φ R ], we link the state V l = V L to the state V s,l solving equation (37) when the porosity changes from φ L to φ s where we choose the supersonic solution. We then compute the V s,r subsonic state using relation (40) of the 1 -s stationary shock. Finally we determine with equation ( 35) the porosity φ min such that we obtain a sonic state V r . The point is that V r should be sonic with φ min = φ R so we consider the function

φ s → f (φ s ) = φ min -φ R ,
and we seek φ s such that f (φ s ) = 0. We use a Lagrange algorithm where we initialise the algorithm with φ = φ L and φ = φ R . Assume that we have determined φ s and V r , then we check that P r < P 1 such that we can construct a rarefaction.

Computation of the other states. We compute the state V a after the 2 -w contact discontinuity using the ρ 2 density. We finally determine state V R using the pressure P 3 . If P 3 > P a , we have a rarefaction while we have a shock if P 3 < P a . In this last case, we have to check that σ 3 > u a to obtain an admissible RRR 1 configuration.

Numerical Results

Numerical investigations have been carried out to test the numerical scheme based on the Rusanov flux and the nonconservative flux given in section 3. Two sets of tests are proposed: the first set of tests aims to check the ability of the numerical method to solve the Riemann problem for several characteristics situations (rarefaction and resonant configurations) while a second set of tests is dedicated to the comparison between the first-and second-order schemes with a regular porosity function to check the performance of the decomposition into nonconservative flux and source term. Computations are performed using the OFELI library of Touzani (1998Touzani ( -2003) ) to handle the mesh.

New Sod tests

We consider a Sod tube on domain [0,2] where the initial condition discontinuity is located at point x = 0.8. All the simulations have been performed using a uniform subdivision of 800 elements and we adapt the time step to guarantee the scheme stability using the CFL condition deriving from the conservative part of the flux:

∆t ≤ ∆x 2 max i (|u i | + c i ) . (41)
First-order schemes are stable if the CFL condition is respected but we have to cut by two (and sometimes by three) the time step to carry out simulations with the second-order scheme to preserve stability. We face to an important number of configurations and we have selected a representative set of situations we shall compare to the exact solution obtained with the inverse Riemann problem. Configurations A, B, C and D are the simpler cases where the four waves are clearly isolated. Such a situation have been yet studied [START_REF] Andrianov | On the solution to the Riemann problem for the compressible duct flow[END_REF] and we only present the comparison between the numerical solution and the exact solution with the configuration C. The main point is to check that the contact discontinuity consecutive to porosity jump is well calculated. We then proceed with the LR and RR configurations when the 1 -w is splitted into a rarefaction reaching the sonic point on one side and a second contribution of the 1 -w wave on the other side. We consider the two configurations whether φ L < φ R or φ L > φ R and for each case, two simulations have been carried out where the second part of the 1 -w wave is a rarefaction or a shock. We finally deal with the resonant configuration where a 1 -s stationary shock splits the 0 -w contact discontinuity. We have first considered the R 1 situation where the 1 -w is reduced to the 1 -s stationary shock. We then proceed with more complex configurations like LRR 1 and RRR 1 . Numerical simulations have been carried out with φ L < φ R to obtain a left rarefaction reaching the sonic point and a 1 -s stationary shock sharing the 0 -w discontinuity at an intermediate porosity φ s ∈ [φ L , φ R ]. A similar configuration has been studied when φ L > φ R . 1 gives the states obtained with the inverse Riemann problem while Fig. 11 shows the comparison between the exact solution and the numerical approximations for the C configuration. Viscosity effects of the Rusanov flux are strongly reduced by the MUSCL procedure and we obtain an accurate approximation of the solution. In particularly, we note that the 0w contact discontinuity situated at x = 0.8 is correctly resolved.

Configuration C

Configuration LR 1 with a rarefaction for the second part of the wave

Comments. Table 2 gives the states obtained with the inverse Riemann problem while Fig. 12 shows the comparison between the exact solution and the numerical approximations for the LR1 configuration. We note that the constant state situated just after the contact discontinuity is not well-approximated whereas all the other states are well-evaluated. Our point of view is that the rarefaction on the left side of the interface x = 0.8 is not well-evaluated and does not reach the sonic point inducing a poor approximation of the constant state on the right side.

All the numerical experiences we have carried out for such a situation indicate a high sensitivity of the right state V r with respect to the approximation of the left rarefaction accuracy. Comments. Table 3 gives the states obtained with the inverse Riemann problem while Fig. 13 shows the comparison between the exact solution and the numerical approximations for the LR 1 configuration. The goal is to numerically check that the second part of the 1w wave can be a rarefaction or a shock. As in the previous situation, the poor approximation of the rarefaction to the sonic point induces a poor approximation of the constant state situated on the right of the interface.

Configuration R 1

Comments. Table 5 gives the states obtained with the inverse Riemann problem while Fig. 15 shows the comparison between the exact solution and the numerical approximations. We observe that the numerical approximation suits very well with the exact solution. We also note that the scheme does not detect the intermediate states V sl and V sr we have introduced to compute the 1s stationary shock but reproduces the state V l and V r . Comments. Table 6 gives the states obtained with the inverse Riemann problem while Fig. 16 shows the comparison between the exact solution and the numerical approximations for the LRR1 configuration. Like in the LR1 case, the major difficulty is the computation of the rarefaction situated on the left of the interface leading to a poor approximation of the constant state located at the right of the interface, the other states being well-approached. Such a configuration is difficult to numerically solve and shows the limit of the Rusanov scheme we have proposed. New schemes have to be considered to correctly solve such a configuration.

Configuration RRR 1

Comments. Table 7 gives the states obtained with the inverse Riemann problem while Fig. 17 shows the comparison between the exact solution and the numerical approximations for the RRR 1 configuration. Like in the LRR 1 the main difficulty is the rarefaction approximation. We observe a different behaviour of the rarefaction between the LRR and the RRR configuration.

General comments on the Sod tests

In the major case, the contact discontinuity 0 -w is well-computed and the second-order scheme manages to strongly reduce the diffusion effect of the Rusanov flux. Since the porosity is constant on both side of the interface, the porosity splitting φ = φ r + φ d is reduced to φ = φ d and we keep a firstorder scheme in the vicinity of the interface whereas we have a full secondorder scheme on both side of the interface. We would like to highlight the difficulty to compute an accurate solution when the configuration involves a rarefaction which reaches the sonic point at the interface (configuration LR, RR, LRR and RRR). A small error on the rarefaction approximation induces an important error on the constant state situated on the other side of the interface. We also mention that the scheme fails to compute an accurate solution when the porosity jump is to high (more precisely when the rate of porosity is far from one). Indeed, we have proved that the scheme consistency error is of order |∆φ| 1 which leads to a poor accuracy with too large porosity discontinuities even if we use a finer mesh.

= {0 -w, 1 -s, 0 -w, 1 -r, 2 -w, 3 -s}. φ ρ (kg.m -3 ) u (m.s -1 ) P (Pa) Mach V L 1.0 1.3 2 1 1.92725 V s,l 0.

Numerical tests with regular porosity

The second set of numerical tests consists in computing approximations of the stationary solution using a regular porosity. We check the ability of the scheme to correctly reproduce stationary solutions satisfying the wellbalanced criterion, and we measure the scheme accuracy.

A stationary solution V (x) has to respect the steady-state equation (6) which is equivalent to the three invariants:

For a given state V L = (φ(0), ρ L , u L , P L ), the steady-state solution at any point x is given by V (x) = (φ(x), ρ(x), u(x), P (x)) where ρ(x) is a solution of g(ρ; S, H) = -D φ(x) 2 , coefficients D, H and S being computed with V L . If V L is a subsonic state, we use the subsonic branch of g while we seek a supersonic solution if V L is a supersonic state.

To perform the simulation, we use a regular mesh T h of spatial step h = 1 I+1 where I is the cell number. For a given left state V L at point x = 0, we compute the right state V R at point x = 1 solving equation g. In the test we propose, we have immediately V R = V L since φ(0) = φ(1). We then compute the solution using the finite volume scheme till we converge to a stationary solution using the left and right condition and we denote by V i the approximation of the steady-state solution on cell K i .

To measure the scheme accuracy and the convergence to the steady-state solution, we compute D i , H i and S i for each cells and we set

D err = max i |D L -D i |, H err = max i |H L -H i |, S err = max i |S L -S i |.
Of course, if one has D err = H err = S err , we have exactly obtain the steadystate solution.

In a first test, we perform the simulation for a subsonic and a supersonic configurations. For each configuration, we have used the first-order scheme and the second-order scheme using the minmod and the vanleer limiter. We have also performed a second-order simulation cancelling the limiting procedure. We would like to mention that we enforce the Dirichlet condition on cells 1, 2, I -1 and I to obtain the second-order accuracy. Indeed, if one only enforces the states at cells 1 and I, the second-order accuracy can not be attained.

Subsonic stationary case

We first deal with the subsonic configuration where we use the left state ρ L = 1, u L = 100 and P L = 100000 and use four different meshes of size h = 0.01, 0.005, 0.0025 and 0.00125 respectively. Table 8 presents the D and H errors with respect to the mesh. We obtain a slow convergence (lower than one) with the first-order scheme but a very good accuracy with the second-order schemes. We note important differences with the choice of the limiter and super convergence effect when the limiter procedure is cancelled (remember that we use a uniform mesh). The great difference of accuracy between the first-and second-order is noticeable and proves that (in this case at least) the nonconservative term splitting is efficient. We draw the convergence curves using the logarithm scale in Fig. 18 where the errors uniformly decrease with respect to h. We mention that the steady-state solution is rapidly obtained with a first-order scheme whereas we have to compute the unstationary solution for a very larger time interval to reach the solution with a second-order scheme. For example, the time computation of the steady-state solution without limiter procedure is more than one hundred longer that the time to compute the steady-state solution with the first-order scheme.

Supersonic stationary case

We now deal with the supersonic configuration where we use the left state ρ L = 1, u L = 2000, P L = 100000. The supersonic case is easier to compute from a numerical point of view and steady-state are reached more rapidly than in the subsonic case. Like in the subsonic situation, we obtain in Table 9 a convergence order of one with the first-order scheme while the secondorder schemes converge with an order greater than two. We find again a super convergence effect with the scheme where the limiter procedure is cancelled. Like in the subsonic case, we observe in Fig. 19 a regular convergence of the D and H errors.

Conclusions

In this paper, we have presented a generic method to construct wellbalanced scheme for the nonconservative Euler system based on a conservative flux and a nonconservative flux such that the steady-state solutions are preserved. We have also introduced the porosity splitting principle to construct second-order scheme. The technique has been applied to the Rusanov conservative flux where we deduce the nonconservative flux to satisfy the well-balanced criterion. The second part of the article is dedicated to a complete description of the Riemann problem where we focus on the resonant situations. We have listed 14 configurations and proposed a technique to solve the inverse Riemann problem. At last, we have presented two sets

Figure 1 :

 1 Figure 1: An example of configuration : C = {1r; 0w; 2w; 3 -s}.

Figure 2 :

 2 Figure 2: Admissible configurations respecting the sign and the stable configuration criterions.

Figure 3 :

 3 Figure 3: Function g : the minimum is attained at point ρ son .

Figure 7 :

 7 Figure 7: Configuration R 1 is a transition between configurations C and D with a 1s shock.

Figure 8 :

 8 Figure 8: Resonant configuration R 1 . The stationary 1-s shock occurs at an intermediate porosity φ s .

Figure 9 :

 9 Figure 9: Left Rarefaction and Resonant configuration LRR 1 with φ L < φ R . The 1w is constituted of a rarefaction on the left and a 1s stationary shock sharing the 0w contact discontinuity with an intermediate porosity φ s .

Figure 10 :

 10 Figure 10: Right Rarefaction and Resonant configuration RRR 1 with φ L > φ R . The 1w is constituted of a rarefaction on the right and a 1s stationary shock sharing the 0w contact discontinuity with an intermediate porosity φ s .

Figure 11 :

 11 Figure 11: Configuration of type C = {1r, 0w, 2w, 3 -s}.

Figure 13 :

 13 Figure 13: An example of configuration : LR1 = {1r, 0w, 1s, 2w, 3 -s}.

Figure 14 :

 14 Figure 14: An example of configuration : RR1 = {1r, 0w, 1r, 2w, 3 -s}.

Figure 15 :

 15 Figure 15: An example of configuration : R1 = {0w, 1r, 0w, 2w, 3 -s}.

  (a) D = φρu, (b) H = u 2

Figure 16 :

 16 Figure 16: An example of configuration : LRR1 = {1r, 0w, 1s, 0w, 2w, 3 -s}.

Figure 17 :

 17 Figure 17: An example of configuration : RRR1 = {0w, 1s, 0w, 1r, 2w, 3 -s}.

Figure 18 :

 18 Figure 18: Subsonic stationary case : D and H errors.

Figure 19 :

 19 Figure 19: Supersonic stationary case : D and H errors.

  is a decreasing function on interval [0, ρ son ] and an increasing function on interval [ρ son , +∞[. Let denote by u son , P son and c son the respective velocity, pressure and sound velocity associated to ρ son , then u 2 son = c 2 son . Moreover, we have

Table 1 :

 1 Configuration of type C = {1r, 0w, 2w, 3 -s}.

	φ ρ (kg.m -3 ) u (m.s -1 ) P (Pa)	Mach
	V L 0.9	3.6	100	300000 0.29277
	V l 0.9	2.69478	196.113 200000 0.608399
	V r 1.0	2.82888	168.135 214071 0.51656
	V a 1.0	3.4	168.135 214071 0.566308
	V R 1.0	3.23885	153.785 200000 0.523034
	Comments. Table

Table 2 :

 2 Configuration of type LR1 = {1r, 0w, 1r, 2w, 3 -s}. Configuration of type LR 1 with shock for the second part of the wave

	φ ρ (kg.m -3 ) u (m.s -1 ) P (Pa)	Mach
	V L 0.8	5	250	400000 0.747018
	V l 0.8	4.03113	320.553 295869	1
	V r 0.8	2.26274	456.86	131823 1.59971
	V a 0.8	1.58382	555.19	80000 2.08778
	V b 1.0	1.68	555.19	80000 2.15024
	V R 1.0	2.3764	647.909 130000 2.3412
	7.1.3.

Table 3 :

 3

	Mach

Configuration of type LR1 = {1r, 0w, 1s, 2w, 3 -s}.

φ ρ (kg.m -3 ) u (m.s -1 ) P (Pa)

Table 4 :

 4 Configuration of type RR1 = {1r, 0w, 1r, 2w, 3 -s}.

	φ ρ (kg.m -3 ) u (m.s -1 ) P (Pa)	Mach
	V L 1.0	3.6	0	400000	0
	V a 1.0	2.50105	138.545 240219 0.377818
	V s 1.0	1.70111	339.491 140043 0.999999
	V b 1.0	1.33739	419.224 100000 1.29572
	V c 0.6	0.7	419.224 100000 0.937414
	V R 0.6	0.57037	329.131	75000	0.7671

Table 5 :

 5 Configuration of type R1 = {0w, 1s, 0w, 2w, 3 -s}.

		φ	ρ (kg.m -3 ) u (m.s -1 ) P (Pa)	Mach
	V L	1.0	1.0	500	10000 1.33631
	V s,l 0.979139	1.05179	485.508 107325 1.28454
	V s,r 0.979139	1.56585	326.119 188717 0.793928
	V r	0.95	1.45327	362.159 170000 0.89492
	V a	0.95	1.2	362.159 170000 0.813207
	V R	0.95	1.34771	414.462 200000 0.909292
	7.1.6. Configuration LRR 1

Table 6 :

 6 Configuration of type LRR1 = {1r, 0w, 1s, 0w, 2w, 3 -s}.

		φ	ρ (kg.m -3 ) u (m.s -1 ) P (Pa)	Mach
	V L	1.3	1.862	0.826	2.4583 0.607559
	V l	1.3	1.32768	1.27062 1.53106	1
	V s,l 1.31102	1.21448	1.37737 1.35148 1.10351
	V s,r 1.31102	1.42713	1.17214 1.69479 0.909049
	V r	1.6	1.78931	0.766029 2.32609 0.567819
	V a	1.6	2	0.766029 2.32609 0.60032
	V R	1.6	1.79564	0.629806	2	0.504356

Table 7 :

 7 Configuration of type RRR1

Table 8 :

 8 Subsonic stationary case : D and H errors.

		h	10. 10 -3 5. 10 -3 2.5 10 -3 1.25 10 -3 order
		first-order	19.2	12.31	7.48	4.31	0.72
	(D error)	minmod	3.97	0.996	0.231	0.0561	2.04
		vanleer	1.03	0.149	0.027	0.0080	2.33
		no limiter	0.45	0.061	0.0084	0.0013	2.81
		h	10. 10 -3 5. 10 -3 2.5 10 -3 1.25 10 -3 order
		first-order 6454.6 3465.05 1844.8	969.4	0.91
	(H error)	minmod	126.45	32.606	13.77	5.29	1.52
		vanleer	64.92	12.48	2.25	0.303	2.58
		no limiter	68.10	10.68	1.24	0.157	2.92

Table 9 :

 9 Supersonic stationary case : D and H errors.

		h	10. 10 -3 5. 10 -3 2.5 10 -3 1.25 10 -3 order
		first-order 162.16	80.92	40.16	19.98	1.00
	(D error)	minmod	19.31	4.89	1.227	0.306	1.99
		vanleer	8.62	2.094	0.516	0.128	2.02
		no limiter	1.34	0.242	0.0515	0.012	2.26
		h	10. 10 -3 5. 10 -3 2.5 10 -3 1.25 10 -3 order
		first-order 24381.5 10845.4 5047.2	2424.2	1.10
	(H error)	minmod	2975.1	665.2	158.3	38.612	2.09
		vanleer	244.2	27.8	3.331	0.407	3.07
		no limiter	180.0	20.72	2.459	0.299	3.07

Acknowledgments

The authors would like to thank the Schneider Electric company for its financial support.

The inverse Riemann problem

Let us consider the Riemann problem with initial left and right states V L and V R and assume that there exists an autosimilar solution. We characterize the solution by its configuration C and the set (P 1 , ρ 2 , P 3 , φ R ) where P 1 is the pressure of the constant state situated on the right of the 1 -w wave, ρ 2 is the density of the constant state situated on the right of the 2 -w wave and P 3 is the pressure of the constant state situated on the right of the 3 -w wave. We then define an application

For the inverse Riemann problem, we proceed in a different way (see [START_REF] Andrianov | Analytical and numerical investigation of two-phase flow[END_REF], [START_REF] Andrianov | On the solution to the Riemann problem for the compressible duct flow[END_REF]). We assume that the left state V L is known and we want to determine the right state V R (if it is possible) with the prescribe configuration C and intermediate values

Note that for a given set (P 1 , ρ 2 , P 3 , φ R ), we can obtain different solutions for different configurations and some configurations have no solution.

Parameterization of the simple waves

To compute each intermediate state for the inverse Riemann problem, we shall use a parameterization of the k -w waves in the simpler cases but configurations like R k , LR k , RR k , LRR k and RRR k require some iterative algorithm we present in the sequel. We introduce some new notations to provide a parameterization of the simple waves. For a given k -w wave, we denote by V k,l and V k,r the left and the right states such that V k,l is linked to V k,r by the one-parameter characteristic curve in the phase space. We construct the inverse Riemann problem solution in the following way:

• For the 1 -w and 3 -w waves, the right state is fixed by the pressure setting P 1,r = P 1 or P 3,r = P 3 .

• For the 2 -w wave, the right state is fixed by the density setting ρ 2,r = ρ 2 .

• For the 0 -w wave, the right state is fixed by the porosity setting φ 0,r = φ R .

Configurations RR 1 and RR 3

The configurations RR are more complex than configurations LR because we have to guess the intermediate subsonic state V l on the left of the interface such that V r is exactly a sonic state.

Computation of the V l and V r states. A Lagrange-like iterative procedure has to be employed to compute the intermediate subsonic state V l . To this end we introduce the intermediate pressure P l and denote by V l the corresponding state obtained with a 1 -r rarefaction if P l < P L or a 1 -s shock if P l > P L . For such a V l , we compute the associated V r on the other side of the interface corresponding to the right porosity φ R where we solve equation ( 34) taking the subsonic solution ρ r = ρ sub . We then consider the function

If one has f (P l ) = 0, the V r state is sonic and we obtain an admissible state V l . To compute P l we adapt the Lagrange algorithm to function f . We take P L and P 1 as initial guess to perform the computation. If the algorithm converges we obtain an approximation of V l and V r .

Computation of the other states. Once we have P r , we check that P 1 < P r since we must have a rarefaction and we compute the intermediate state V a . If P 1 > P r configuration RR 1 is not available. The next operations are the following.

-We compute the state V b after the contact discontinuity prescribing the density ρ 2 .

-If P 3 < P b we have a shock of velocity σ 3 whereas we have a rarefaction if P 3 > P b . We compute the state V c using the parametrization. In case of a shock, one has to check that σ 3 > u b .

The state V c corresponds to the right state V R we shall use in the simulations.

Configurations R 1 and R 3

We now deal with the resonant configurations. We first consider the simpler cases when a k -s shock splits the 0 -w wave. We shall employ the notations introduced in Fig. 8. The main difficulty is to determine the intermediate porosity such that we have a stationary shock between V l,s and V r,s . We only treat the case R 1 but the R 3 configuration is similar. Computation of φ s and the V r state. First of all, if state V l is not supersonic then configuration R 1 is not avalaible. Now, for a given φ s , we compute the supersonic state V l,s linking V l when the porosity changes from φ L to φ s solving equation ( 37). With V l,s in hand, we compute V r,s using the stationary 1 -s shock parameterization :

After some algebric manipulation we get

On the other hand, we have the relation ρ s,r = ρ s,l P s,l (γ -1) + P s,r (γ + 1) P s,l (γ + 1) + P s,r (γ -1) ,

and combining the two equations, we obtain

We then deduce ρ s,r and u s,r , hence the state V s,r . We compute the density ρ r taking the subsonic solution of equation ( 38). We then deduce the state V r . The intermediate porosity has to be fixed such that P r = P 1 where P 1 is the prescribed pressure. To this end, we introduce the function

and we seek φ s such that f (φ s ) = 0. We use a Lagrange method where we initialise the algorithm with φ = φ L and φ = φ R .

Computation of the other states. Assume that state V r is well-calculated, we compute the state V a after the 2 -w contact discontinuity using the ρ 2 density. We finally determine state V R using the pressure P 3 . If P 3 > P a we have a rarefaction while we have a shock if P 3 < P a . In this last case, we have to check that σ 3 > u a to obtain an admissible R 1 configuration.

Configuration of type RR 1

Comments. Table 4 gives the states obtained with the inverse Riemann problem while Fig. 14 shows the comparison between the exact solution and the numerical approximations for the RR 1 configuration. The numerical simulations we have carried out seem to indicate that the solution is less sensitive to the rarefaction approximation situated on the right of the interface. The case φ L > φ R provides a better numerical solution accuracy than the case φ L < φ R . As in the LR 1 case, we can numerically check that we can also have a rarefaction or a shock for the second part of the 1w when φ L > φ R .

of numerical tests to show the scheme ability to compute approximations for the Riemann problem and we have measured the convergence order in case of a steady-state solution with a smooth porosity function.