Frédéric Coquel
email: frederic.coquel@cmap.polytechnique.fr

Khaled Saleh
email: saleh@ann.jussieu.fr

Nicolas Seguin
email: nicolas.seguin@upmc.fr

A Robust and Entropy-Satisfying Numerical Scheme for Fluid Flows in Discontinuous Nozzles

Keywords: Discontinuous nozzle flows, relaxation techniques, Riemann problem. AMS subject classifications : 76S05, 35L60, 35F55

We propose in this work an original finite volume scheme for the system of gas dynamics in a nozzle. Our numerical method is based on a piecewise constant discretization of the crosssection and on a approximate Riemann solver in the sense of Harten, Lax and van Leer. The solver is obtained by the use of a relaxation approximation that leads to a positive and entropy satisfying numerical scheme for all variation of section, even discontinuous with arbitrary large jumps. To do so, we introduce in the first step of the relaxation solver a singular dissipation measure superposed on the standing wave which enables us to control the approximate speeds of sound and thus, the time step, even for extreme initial data.

Introduction

The design of stable and accurate numerical schemes for hyperbolic systems is still a difficult problem and the challenge becomes much more difficult in presence of stiff source terms. Such an issue may occur in the frame of flows which are influenced by external effects, due for instance to the surrounding domain, another fluid, external forces... We are interested here in the numerical approximation of the solutions of a model describing one-dimensional barotropic flows in a nozzle. In this model, ρ and w are respectively the density and the velocity of the fluid while α stands for the cross-section of the nozzle, which is assumed to be constant in time. Under the classical assumption that α (and its variations) is small with respect to a characteristic length in the mainstream direction, the flow can be supposed to be one-dimensional and described by the following set of 1 partial differential equations:

∂ t (αρ) + ∂ x (αρw) = 0, ∂ t (αρw) + ∂ x (αρw 2 + αp(τ )) -p(τ )∂ x α = 0, (1.1)
where τ = ρ -1 is the specific volume and τ → p(τ ) is a barotropic pressure law. The first equation is the classical conservation of mass and the second equation governs the dynamics of the horizontal mean momentum.

The set of equations (1.1) can also model the dynamics of a fluid flow in a porous medium. In such a case, α represents the local porosity in which the fluid evolves. While classical derivations in the context of a flow in a nozzle assume the cross-section α to be smooth, the porosity can achieve large discontinuities. Moreover, this model can be seen as a first step in the design of a robust numerical scheme for two-phase flows, see [START_REF] Ambroso | Relaxation and numerical approximation of a two-fluid two-pressure diphasic model[END_REF] and [START_REF] Coquel | A Robust Entropy-Satisfying Finite Volume Scheme for the Isentropic Baer-Nunziato Model[END_REF], where α stands for the void fraction. As a consequence, we are interested in a numerical scheme which complies with the following requirements, even for discontinuous α with arbitrary large jumps:

1. robustness: it has to preserve the positivity of the density ρ, 2. consistency and nonlinear stability: it has to satisfy discrete entropy inequalities, 3. efficiency: the two previous requirements have to be achieved in such a way that the time step remains acceptable.

It is clear that the two first points are not so easy to obtain, but, in our opinion, the third point is probably the most challenging. Moreover, we want to obtain a numerical scheme which is easy to implement. Let us emphasize that we do not pay too much attention to accuracy issues in this paper, the scheme we propose is clearly first-order, we only focus on these stability properties.

Nevertheless, an increase of the accuracy could be proposed using high order reconstruction methods for instance (note that such techniques cannot help us to answer to the difficulties mentioned above).

We here propose a new numerical scheme for gas dynamics in a nozzle which meets all these requirements. Up to our knowledge, this is the only scheme whith such properties, except the Godunov scheme based on the exact Riemann problem, which cannot be considered for industrial applications [START_REF] Lefloch | The Riemann problem for fluid flows in a nozzle with discontinuous cross-section[END_REF][START_REF] Han | Exact Riemann solutions to compressible Euler equations in ducts with discontinuous cross-section[END_REF]. The three main ingredients which enable us to obtain them are: a piecewise constant discretisation of the cross-section α, the construction of a simple Riemann problem to compute the numerical fluxes, a dissipative modification of the underlying system to allow to deal with severe cases.

The idea of taking a constant-by-cell cross-section goes back to the works of LeRoux and co-workers [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF][START_REF] Gosse | Un schéma-équilibre adapté aux lois de conservation scalaires non-homogènes[END_REF] and to the paper of Isaacson and Temple [START_REF] Isaacson | Convergence of the 2 × 2 Godunov method for a general resonant nonlinear balance law[END_REF], which have been extended by Gosse and co-workers [START_REF] Gosse | A well-balanced scheme using non-conservative products designed for hyperbolic systems of conservation laws with source terms[END_REF][START_REF] Gosse | Localization effects and measure source terms in numerical schemes for balance laws[END_REF][START_REF] Amadori | Godunov-type approximation for a general resonant balance law with large data[END_REF]. The consequence of such a discretisation is to concentrate the source term at the interfaces of the mesh and to ease the construction of well-balanced schemes. In all these pioneering works, the numerical fluxes are obtained by solving each interfacial Riemann problem exactly, which is not an easy task because of the presence of a singular source term. Several attempts to simplify this Riemann solver have been proposed (see for instance [START_REF] Gallouët | Some approximate Godunov schemes to compute shallow-water equations with topography[END_REF]), but the overall resulting numerical scheme may lack for stability properties for severe test cases.

The numerical scheme we construct is based on an approximate Riemann solver, in the spirit of Harten, Lax and van Leer [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]. It is only composed by constant states separated by discontinuities, which makes its practical implementation easy. In order to ensure the positivity of the density and the decrease of the total energy, we interprete this simple solver as the exact solver of a relaxation approximation of system (1.1), following [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF][START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF]. Let us mention that several numerical schemes have been developed on the basis of exact Riemann solvers but without the constant-by-cell discretisation of the cross-section. Generalized Riemann problems and/or non constant-by-cell discretisations of the unknown have to be considered [START_REF] Glaz | The asymptotic analysis of wave interactions and numerical calculations of transonic nozzle flow[END_REF][START_REF] Glimm | A generalized Riemann problem for quasi-onedimensional gas flows[END_REF][START_REF] Ben-Artzi | An upwind second-order scheme for compressible duct flows[END_REF][START_REF] Huang | A conservative, piecewise-steady difference scheme for transonic nozzle flow[END_REF][START_REF] Chen | Global solutions to the compressible Euler equations with geometrical structure[END_REF].

The cornerstone of this scheme is the resolution of the Riemann problem associated with the homogeneous relaxation model for arbitrary data. Even if the relaxation approximation provides a linearly degenerate system, the resonance phenomenon persists since the source term is singular (the cross-section is discontinuous at the interface). In few words, resonance in hyperbolic systems consists in the superimposition of an acoustic wave on the discontinuity of the cross-section, also called the standing wave, leading the associated eigenvectors to be colinear (as a consequence, the system is no longer hyperbolic). In the frame of the original model of gas dynamics in a nozzle, the resonance causes nonuniqueness of the solution of the Riemann problem, as proved in [START_REF] Isaacson | Convergence of the 2 × 2 Godunov method for a general resonant nonlinear balance law[END_REF] and [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF] (see also [START_REF] Hong | The generic solution of the Riemann problem in a neighborhood of a point of resonance for systems of nonlinear balance laws[END_REF][START_REF] Hong | An extension of Glimm's method to inhomogeneous strictly hyperbolic systems of conservation laws by "weaker than weak" of the Riemann problem[END_REF][START_REF] Lefloch | The Riemann problem for fluid flows in a nozzle with discontinuous cross-section[END_REF][START_REF] Han | Exact Riemann solutions to compressible Euler equations in ducts with discontinuous cross-section[END_REF]). In such cases, the superposed shock wave introduces a singular dissipation concentrated in the standing wave. Here, the troubles are different. Global existence still remains true but, for particular initial data, measure solutions have to be considered. They naturally appear when resonance occurs, as parts of the solution in some limit regimes for given patterns of solutions, but we are able to circumvent these solutions by slightly perturbing the relaxation coefficient a which governs the acoustic part of the relaxation model. As a result, the solutions we consider for the final numerical scheme belong the the classical setting of piecewise constant solutions seperated by linearly degenerate waves. Nevertheless, we authorize ourselves to consider a singular dissipation term concentrated in the standing wave, as it may occurs in the original setting. It is woth noting that it only intervenes for severe initial data: resonant cases, high jumps of cross-section, low densities. The main consequence of the use of this possible singular dissipation is a strong control of the time step for all data. This last feature is the most original part of our work.

The outline of this paper is the following. The next section is devoted to the presentation of the main features of the model for gas dynamics in a nozzle. Section 3 is the core of this work: the relaxation approximation is presented and the associated Riemann problem is solved. With the help of this analysis, the numerical approximation is studied in Section 4. Properties of the scheme are described with a special care about the positivity of the density and the non-linear stability. A preliminary numerical test is also presented at the end. Two appendices about some technical developments complete this work.

Let us stress that the purpose of this work is to propose an entropy satisfying numerical scheme with positivity properties. The main consequence is that we do not focus on the uniqueness of the solution of the relaxed Riemann problem (and actually, it fails in the presence of resonance as expected [START_REF] Isaacson | Convergence of the 2 × 2 Godunov method for a general resonant nonlinear balance law[END_REF]). Moreover, we make an extensive description of the Riemann solver and its properties, therefore only one numerical test is presented. A careful numerical study of our scheme is the subject of a forthcoming work.

2 The Euler equations in a nozzle with variable cross-section

Presentation and main properties

The model describing one-dimensional barotropic flows in a nozzle can be described by the following set of partial differential equations:

   ∂ t α = 0, ∂ t (αρ) + ∂ x (αρw) = 0, ∂ t (αρw) + ∂ x (αρw 2 + αp(τ )) -p(τ )∂ x α = 0, (2.1)
where τ = ρ -1 is the specific volume and τ → p(τ ) is a barotropic pressure law. The first equation expresses the constancy of the section α, while the second and the third equations are respectively the mass and the momentum equations. All along this paper, we assume that the pressure p is a smooth function of τ satisfying the following classical properties. For all τ > 0, p(τ ) > 0, p (τ ) < 0, with lim τ →0 p(τ ) = +∞ and lim τ →+∞ p(τ ) = 0. An example of such a pressure law is an ideal gas barotropic pressure law p(τ ) = Sτ -γ with S > 0 and γ > 1. System (2.1) takes the following condensed form:

∂ t U + ∂ x f(U) + c(U)∂ x U = 0, (2.2)
where U = (α, αρ, αρw) T is the vector of unknowns and the functions f and c are given by

f(U) =   0 αρw αρw 2 + αp(τ )   , c(U)∂ x U =   0 0 -p(τ )∂ x α   . (2.3)
In practice, the constant section α is determined once and for all by the initial condition, and thus it is not properly speaking an unknown function. However, the section α appears in the mass and momentum equations, especially in the pressure terms. Therefore, in the numerical simulations, where the solutions of system (2.1) are approximated by a Finite Volume Method, it is more appropriate to consider α as an unknown function, since it allows us to use the convenient machinery of hyperbolic systems theory. In particular, we will be able to construct (approximate) self-similar solutions to system (2.1) (i.e. solution depending only on x/t). The following proposition holds, that characterizes the fields of this system. Proposition 2.1. For any U in the phase space Ω defined by

Ω = U = (α, αρ, αρw) T ∈ R 3 , α > 0, αρ > 0 , (2.4)
system (2.1) admits the three following eigenvalues

σ 0 (U) = 0, σ 1 (U) = w -c(τ ), σ 2 (U) = w + c(τ ), (2.5)
where c(τ ) = τ -p (τ ) is the speed of sound. The system is hyperbolic on Ω (i.e. the corresponding right eigenvectors span R 3 ) if, and only if (w -c(τ ))(w + c(τ )) = 0. If the latter condition is not fulfilled, the system is said to be resonant. Moreover, the characteristic field associated with σ 0 is linearly degenerate, while the characteristic fields associated with σ 1 and σ 2 are genuinely non-linear.

Proof . The proof is classical (see for instance [START_REF] Lefloch | The Riemann problem for fluid flows in a nozzle with discontinuous cross-section[END_REF]).

The phase space Ω introduced in (2.4) is the physically relevant domain where the solutions of (2.1) have to lie. Indeed, the section α has to be positive (which is trivially imposed by the initial condition) as well as the fluid density ρ. In the sequel Ω will be referred to as the phase space of positive solutions. As regards the smooth solutions of system (2.1), we have the following property: Proposition 2.2. The smooth solutions of (2.1) obey the following additional conservation law

∂ t (αρE) + ∂ x (αρEw + αp(τ )w) = 0, (2.6
)

where E = w 2 2 + e(τ )
is the total energy and where the function τ → e(τ ) is given by e (τ ) = -p(τ ).

Moreover, the function (α, αρ, αρw) → αρE is convex.

Proof . Equation (2.6) follows from classical manipulations of system (2.1). The convexity of αρE as a function of U is a consequence of the convexity of the internal energy e with respect to τ .

When one considers non-smooth weak solutions of system (2.1), it is well known that there is no uniqueness of such solutions and one has to add a so-called entropy selection criterion in order to select the relevant physical solutions of (2.1). Definition 2.1. A weak solution of system (2.1) is said to be an entropy solution if it satisfies the following inequality in the weak sense of the distributions

∂ t (αρE) + ∂ x (αρEw + αp(τ )w) ≤ 0. (2.7)
As the function (α, αρ, αρw) → αρE is convex, this selection criterion can be formally justified by the vanishing viscosity method (see for example [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]). When the solution contains strong shocks, inequality (2.7) is strict, and this accounts for the loss of energy due to viscosity.

Standing wave and resonance

For the sake of numerical applications, one has to consider the case of discontinuous cross-sections α, the simpler example of which is given by a Riemann-type initial condition α(x) = α L if x < 0 and α(x) = α R if x > 0. Since α is constant throughout time, this gives rise to a standing discontinuity across which one has to define jump relations. The main difficulty lies in the treatment of the non-conservative product p(τ )∂ x α since this product cannot be represented in the usual sense of distributions. Nevertheless, in the region of hyperbolicity of system (2.1), i.e. when w = ±c(τ ), this non-conservative product is supported by the standing wave associated with the linearly degenerate field σ 0 = 0, and the natural definition of p(τ )∂ x α is drawn from the conservation of the Riemann invariants associated with σ 0 . A well-known result on hyperbolic systems, see for instance [START_REF] Lefloch | Shock waves for nonlinear hyperbolic systems in non-conservative form[END_REF][START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF], states that for a linearly degenerate field, the Rankine-Hugoniot jump relations are valid for any additional conservation law which is satisfied by smooth solutions. Therefore, the two Riemann invariants are obtained by applying the Rankine-Hugoniot jump conditions to the mass conservation equation and to the conservative equation on the velocity

∂ t w + ∂ w 2 2 + e(τ ) + τ p(τ ) = 0, (2.8)
which leads to

[αρw] 0 = w 2 2 + e(τ ) + τ p(τ ) 0 = 0, (2.9)
where [X] 0 denotes the jump of any quantity X across the standing wave. An easy consequence of these relations is the preservation of the energy flux across the standing wave: [αρEw + αp(τ )w] 0 = 0. In the case of steady states at rest, the jump relations reduce to w = 0 and [ρ] 0 = 0.

(2.10)

When the resonance phenomenon appears, i.e. when there exists in the solution a state (α, αρ, αρw) such that w = ±c(τ ), the hyperbolicity of the system is lost which means that the right eigenvectors are no longer linearly independent, and the standing wave superimposes with a non-linear field associated with one of the extreme eigenvalues σ 1 or σ 2 . In this particular case, defining the non-conservative product is difficult and the uniqueness of solutions is lost in general (even with the entropy criterion given by Definition 2.1), see [START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF][START_REF] Isaacson | Convergence of the 2 × 2 Godunov method for a general resonant nonlinear balance law[END_REF]. Besides, if the standing wave superimposes with a stationary shock, while the Rankine-Hugoniot jump relation

[αρw]

0 = 0 remains valid, the second one of (2.9) is no longer true. As a consequence, the energy flux is no longer preserved across the wave and we rather have

[αρEw + αp(τ )w] 0 < 0 (2.11)
since the energy strictly decreases through the shock.

Numerical approximation and Riemann solvers

One of the most classical approaches for the numerical approximation of the solutions of (2.1) is the so-called well-balanced approach (see [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF][START_REF] Gosse | Un schéma-équilibre adapté aux lois de conservation scalaires non-homogènes[END_REF][START_REF] Chinnayya | A well-balanced numerical scheme for the approximation of the shallow-water equations with topography: the resonance phenomenon[END_REF]) which relies on the construction of the exact solution of system (2.1) for the particular case where the initial condition is given by a constant state U L for x < 0 and a constant state U R for x > 0 (one speaks of a Riemann problem):

U 0 (x) = U L if x < 0, U R if x > 0.
(2.12)

Unfortunately, the exact solution of this Riemann problem is quite uneasy to obtain (see [START_REF] Lefloch | The Riemann problem for fluid flows in a nozzle with discontinuous cross-section[END_REF][START_REF] Andrianov | On the solution to the Riemann problem for the compressible duct flow[END_REF][START_REF] Han | Exact Riemann solutions to compressible Euler equations in ducts with discontinuous cross-section[END_REF]) due to the non-linearities of the pressure law and to the difficulties linked with the resonance phenomenon (definition of the non-conservative product, non-uniqueness...). Therefore, an other approach is preferred, where solving the Riemann problem for system (2.1) is replaced by solving an easier Riemann problem for an enlarged system obtained by a relaxation approximation method.

3 Relaxation approximation

The relaxation system and its main properties

In this section, we propose a suitable relaxation approximation of the entropy weak solutions of system (2.1). For this purpose, we first recall that the genuine non-linearity of the two extreme fields (also referred to as the {σ 1 , σ 2 }-fields in the sequel) is due to the non-linearities of the pressure law τ → p(τ ). In the spirit of [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF], we consider an enlarged system involving an additional unknown T associated with a linearization π of the pressure law. This linearization is designed to get a quasilinear enlarged system, shifting the initial non-linearity from the convective part to a stiff relaxation source term. The relaxation approximation is based on the idea that the solutions of the original system are formally recovered as the limit of the solutions of the proposed enlarged system, in the regime of a vanishing relaxation coefficient ε > 0. As a relaxation approximation of (2.1), we propose the following system:

         ∂ t α ε = 0, ∂ t (αρ) ε + ∂ x (αρw) ε = 0, ∂ t (αρw) ε + ∂ x (αρw 2 + απ(τ, T )) ε -π(τ, T ) ε ∂ x α ε = 0, ∂ t (αρT ) ε + ∂ x (αρT w) ε = 1 ε (αρ) ε (τ -T ) ε , (3.1)
where the linearization of the pressure law is given by

π(τ, T ) = p(T ) + a 2 (T -τ ). (3.2)
System (3.1) takes the following condensed form:

∂ t W ε + ∂ x g(W ε ) + d(W ε )∂ x W ε = 1 ε R(W ε ), (3.3)
where W = (α, αρ, αρw, αρT ) T is the vector of unknowns and the functions g, d and R are given by

g(W) =     0 αρw αρw 2 + απ αρT w     , d(W)∂ x W =     0 0 -π∂ x α 0     , R(W) =     0 0 0 αρ(τ -T )     . (3.4)
To ease the notation hereafter, we will omit the superscript ε . From this point, we will refer to the original system (2.1) as the equilibrium system, while system (3.1) will be referred to as the relaxation system. We can see that in the formal limit ε → 0, the additionned variable T tends towards the specific volume τ , and the linearized pressure π tends towards the original non-linear pressure p, thus recovering the equilibrium system (2.1) in the first three equations of (3.1). The constant a in (3.2) is a constant positive parameter that must be taken large enough to prevent system (3.1) from instabilities in the regime of small values of ε. This will be clarified in section 4.4.

It is relevant to focus on the convective part of system (3.1) since a fractional step method is commonly used in the implementation of relaxation methods: the first step is a time-advancing step using the solution of the Riemann problem for the convective part of (3.1):

       ∂ t α = 0, ∂ t (αρ) + ∂ x (αρw) = 0, ∂ t (αρw) + ∂ x (αρw 2 + απ(τ, T )) -π(τ, T )∂ x α = 0, ∂ t (αρT ) + ∂ x (αρT w) = 0, (3.5)
while the second step consists in an instantaneous relaxation towards the equilibrium system by imposing T = τ in the outcome of the first step. This second step is equivalent to sending ε to 0 instantaneously (see section 4 for details).

We now state the main property that motivates the introduction of the proposed relaxation system: Proposition 3.1. For any W in the phase space Ω r defined by

Ω r = W = (α, αρ, αρw, αρT ) T ∈ R 4 , α > 0, αρ > 0, αρT > 0 , (3.6)
system (3.5) admits the four following eigenvalues

σ r 0 (W) = 0, σ r 1 (W) = w -aτ, σ r 2 (W) = w, σ r 3 (W) = w + aτ, (3.7)
and is hyperbolic on Ω r (i.e. the corresponding right eigenvectors span R 4 ) if, and only if (waτ )(w+aτ ) = 0. If the latter condition is not fulfilled, the system is said to be resonant. Moreover, all the characteristic fields associated with {σ r i } i=0..3 are linearly degenerate.

Proof . The proof can be easily recovered following closely related steps developped in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]. Details are left to the reader.

Note that the crucial property here is the linear degeneracy of the two extreme waves. This enables us to easily define jump relations across these two waves. More precisely, the first equation of (3.5) shows that for any solution of the Riemann problem, the jump of α only occurs through the σ 0 standing wave, therefore α is a Riemann invariant for both acoustic fields. Similarly, equation four in (3.5) shows that T is also a Riemann invariant for these fields and the last Riemann invariant is determined by observing that for any linearly degenerate wave, the eigenvalue is also constant through this field (any other invariant of the field can be expressed as a continuous function of these three Riemann invariants).

Remark 3.1. System (3.5) could be studied for itself without relaxation consideration, i.e. without considering that it is precisely designed to approximate the natural physical system (2.1). In that case, there is no reason to impose the positivity of the density in the solutions and the phase space for (3.5) turns to be larger than Ω r defined in (3.6). For our relaxation approximation purposes though, we ask the solutions of (3.5) to stay within the phase space Ω r . The positivity of the additioned variable T is necessary in order for the relaxed pressure π(τ, T ) to be well-defined. Subsequently, any vector W is said to be positive if it satisfies W ∈ Ω r , and any solution (x, t) → W(x, t) is said to be a positive solution if for all (x, t) in R x × R + t , W(x, t) belongs to Ω r .

Jump relations across the standing wave

We now focus on the definition of jump relations across the standing wave in the PDE model (3.5). Applying the Rankine-Hugoniot jump relation to the mass conservation equation as well as to the transport equation of T yields two Riemann invariants for the standing wave provided that system (3.5) is hyperbolic (see hereafter). But as the non-conservative product π(τ, T )∂ x α is not well defined across the standing wave (π(τ, T ) may not be continuous across this wave), we cannot apply the Rankine-Hugoniot relation to the momentum conservation equation. Moreover, the variable w no longer satisfies a conservative equation as in the equilibrium case (2.8). Instead, we seek an additional conservation law satisfied by the smooth solutions of (3.5) eventually leading to a full set of jump relations. We have the following statement:

Proposition 3.2. The smooth solutions of (3.5) obey the following additional conservation law

∂ t (αρE) + ∂ x (αρEw + απ(τ, T )w) = 0, (3.8)
where the total energy is

E = w 2 2 + e(T ) + π 2 (τ, T ) -p 2 (T ) 2a 2 . (3.9)
Proof . This follows from classical manipulations. The details are left to the reader.

For a hyperbolic conservative system, the conservation of energy (3.8) holds true in the weak sense for any solution presenting only contact discontinuities, and the Riemann invariant obtained by applying the Rankine-Hugoniot jump relation to equation (3.8) can be expressed as a continuous function of the other Riemann invariants. Nevertheless, system (3.5) is not conservative in the neighborhood of the standing wave and this is the reason why applying the Rankine-Hugoniot relation to (3.8) yields a new jump relation. Note that there are no theoretical results that impose relation (3.8) to be exactly maintained across the standing wave when the resonance occurs (i.e. when w = ±aτ ), and we will see that, if equation (3.8) is exactly satisfied in the weak sense, we will not be able to impose the invariance of the domain Ω r . Indeed, it will be proved that keeping domain Ω r invariant requires the decrease of the energy in general. This is related to the fact that Ω r is not the natural space for the solutions of system (3.5) (see Remark 3.1). These considerations motivate the construction of solutions to the Riemann problem where the energy decreases (in the weak sense) across the standing wave:

∂ t (αρE) + ∂ x (αρEw + απ(τ, T )w) ≤ 0, in D , (3.10)
as it may happen for the equilibrium system (see (2.11)).

Solving the Riemann problem for the relaxation system

Definition of the solutions and existence theorem

Let be given W L and W R , two positive states in Ω r . We are now interested in solving the Riemann problem for system (3.5), i.e. we seek solutions satisfying the initial condition

W 0 (x) = W L if x < 0, W R if x > 0. (3.11) Definition 3.1. A solution of the Riemann problem (3.5)-(3.11) is a function W : (x, t) ∈ R x × R + t → W(x, t)
∈ Ω r satisfying the following properties:

1. W is a self similar mapping W(x, t) = W r (x/t; W L , W R ).
2. W is made of constant intermediate states separated by waves whose constant speeds are eigenvalues {σ r i } i=0..3 of the system, and each eigenvalue σ r i appears at most once in W.

3. W is a weak solution on R * x × R + t of    ∂ t (αρ) + ∂ x (αρw) = 0, ∂ t (αρw) + ∂ x (αρw 2 + απ(τ, T )) = 0, ∂ t (αρT ) + ∂ x (αρT w) = 0.
(3.12)

Defining [X]

0 = X(0 + ) -X(0 -) for any function X of x/t, the dissipation of energy across the standing wave is non-positive in the sense that

[αρwE + απw] 0 ≤ 0. (3.13) 5. W satisfies, for all T > 0 and R > 0 such that for all x ≥ R, W(-x, T ) = W L and W(x, T ) = W R , T [αρw] 0 = R -R (αρ)(T, x)dx -R((αρ) L + (αρ) R ) + T ((αρw) R -(αρw) L ). (3.14)
Moreover, we impose the following alternative

• if [αρw] 0 = 0 then [αρwT ] 0 = 0, • if [αρw] 0 = 0 then [T ] 0 = 0.
Let us provide some comments on this definition. The three first points are classical and they also imply that (3.8) is satisfied in R *

x × R + t . The remaining points specify the behavior of the solution through the standing wave.

Denoting f = -[αρwE + απw] 0 , a solution satisfies ∂ t (αρE) + ∂ x (αρEw + απ(τ, T )w) = -f δ 0 .
As mentioned before, we allow a non-zero energy dissipation through the standing wave. In some cases, an appropriate choice of the dissipation rate f will enable us to guarantee the positivity of the intermediate states (c.f. remark 3.1).

Point 5 expresses that a Dirac measure corresponding to a mass concentration is allowed to appear at x = 0. But thanks to (3.14), the global conservation of mass is ensured. Solutions which involve Dirac measures have already been introduced in many works, see for instance [START_REF] Bouchut | Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness[END_REF].

Here, this phenomenon is due to the linear denegeracy of the extreme fields together with the resonance pathology. Indeed, such singular solutions do not appear in the equilibrium problem (see for example [START_REF] Isaacson | Convergence of the 2 × 2 Godunov method for a general resonant nonlinear balance law[END_REF][START_REF] Goatin | The Riemann problem for a class of resonant hyperbolic systems of balance laws[END_REF][START_REF] Hong | The generic solution of the Riemann problem in a neighborhood of a point of resonance for systems of nonlinear balance laws[END_REF][START_REF] Hong | An extension of Glimm's method to inhomogeneous strictly hyperbolic systems of conservation laws by "weaker than weak" of the Riemann problem[END_REF] and previously quoted references).

When [αρw]

0 = 0, the classical conservative jump relation on the PDE for T is no longer available and the alternative which ends point 5 is a natural way define T .

Let us now define some notations depending only on the physical data V L := (ρ L , w L , T L ) and V R := (ρ R , w R , T R ) and that will be useful afterwards:

w := 1 2 (w L + w R ) - 1 2a (π R -π L ), (3.15)
π := 1 2 (π R + π L ) - a 2 (w R -w L ), (3.16
)

τ L := τ L + 1 a (w -w L ) = τ L + 1 2a (w R -w L ) - 1 2a 2 (π R -π L ), (3.17)
τ R := τ R - 1 a (w -w R ) = τ R + 1 2a (w R -w L ) + 1 2a 2 (π R -π L ). (3.18)
In fact, these quantities are respectively the speed, the linearized pressure, and the specific volumes of the solution obtained with a constant initial section α L = α R , provided that the specific volumes τ L and τ R are positive. Let us stress from now on that a will be chosen large for stability matters (see section 4.4) and in particular large enough to enforce the positivity of τ L and τ R . In the sequel, we always assume that the constant parameter a is such that τ L and τ R are positive. This requirement is equivalent to imposing the natural ordering of the waves w L -aτ L < w < w R + aτ R . Thereafter, the self-similar function depicted above will be referred to as the constant section solution. We also introduce the Mach numbers of the intermediate states for the constant section solution:

M L := w L aτ L , M L := w aτ L , M R := w aτ R , M R := w R aτ R . (3.19)
The main result of this section is the following existence theorem for the Riemann problem.

= α L α R is. Remark 3.2.
1. By choosing a large enough, it is always possible to have τ L > 0 and τ R > 0. By the same way, one can also ensure that [αρw] 0 = [αρEw + απw] 0 = 0, which means that there is neither mass concentration nor entropy dissipation at the standing wave.

2. This procedure for the choice of a should be attractive since the solution of the Riemann problem becomes easier to describe (and to compute). However, in the numerical simulations, too much large values of a could imply very small time steps ∆t through the CFL condition and unnecessary numerical diffusion.

3. It is well-known that the Riemann problem for the equilibrium system (2.1) may admit up to three self-similar solutions in the context of resonance. Concerning the relaxation system, several solutions exist near the resonance, in particular when the energy flux is not constant across the standing wave. Actually, since we only aim at constructing a numerical scheme, we did not focus on the uniqueness issue.

4. The continuity of the solution with respect to the data is a very difficult question. On the one hand, the lack of uniqueness is an obvious obstacle to prove any continuity property. On the other hand, in the resonant cases, the solution admits singular parts at the standing wave, which have to be accounted for in the topology used for defining the continuity, which is still unclear for the authors. However, we construct resonant solutions as limits of non-resonant solutions through the computation of the respective intermediate states, when the speed of any wave tends to zero.

The proof of this existence theorem follows from an actual construction of the solution for every given initial data W L and W R . For strictly hyperbolic systems of conservation laws, the characteristic eigenvalues are naturally ordered (see for example Euler's equations). Consequently, if all the characteristic fields are linearly degenerate, the solution is sought in the form of constant states separated by contact discontinuities whose speeds are equal to the corresponding eigenvalues. For system (3.5), the eigenvalues are not naturally ordered because of the existence of a standing wave, and a resonance phenomenon does appear for sonic flows (i.e. flows with vanishing (waτ )(w +aτ )). Therefore, the classical proof must be slightly modified. We first focus on a particular non-resonant ordering of the eigenvalues (for instance w -aτ < 0 < w < w + aτ ) and we determine sufficient conditions (that sometimes appear to be necessary) on the initial states W L and W R for the solution to have this particular ordering. We do the same for the other possible non-resonant orderings (that may be supersonic). Resonant solutions are then studied as limits of non resonant solutions as the speed of any wave tends to zero. Eventually, we check a posteriori that the determined conditions totally cover the entire domain of initial conditions Ω r × Ω r . We show that the conditions that give the ordering of the wave speeds can be expressed in terms of the physical data (V L ,V R ) and of the ratio of left and right sections : ν = α L α R . In addition, for certain values of ν (large or small values depending on the flow direction) the solution may have to dissipate energy in the standing wave in order to preserve the positivity of the densities (again, see Remark 3.1).

The following figure provides a schematic representation of the solution given by Theorem 3.3. It represents the map of the admissible solutions with respect to the initial states W L and W R . The right part of the chart corresponds to the solutions with positive material speed, while the left part depicts the symmetric configurations with negative material speed. The blue lines represent the solutions whose structure needs to refer to a measure concentrated at x = 0.

ν = αL αR ML MR ν = 1 ν = 0 ML = 0 ML = 1 MR = 0 MR = -1 w < 0 w > 0 M R = -ν νM L = 1
In the sequel, a solution of the Riemann problem is said to have signature < i, j > with i and j in {0, 1, 2, 3} if it is composed with i left-going waves and j right-going waves. For example, the solution with the ordering of the eigenvalues w -aτ < 0 < w < w + aτ is said to have signature < 1, 2 >.

Non-resonant solutions

Solutions with signature < 1, 2 >:

Let us now turn on to solutions with the subsonic non-resonant ordering of the eigenvalues w -aτ < 0 < w < w + aτ i.e. solutions with the wave signature < 1, 2 >.

x t w -aτ 0 w w + aτ WL W - W + W3 WR Signature < 1, 2 >
The following result shows that, under a sufficient (and necessary) condition on the initial states, one can build a one-parameter family of solutions (in the sense of Definition 3.1) with signature < 1, 2 >, and the dissipation of energy across the standing wave is directly driven by the underlying parameter. where w and M L are defined respectively in (3.15) and (3.19). These solutions can be parametrized by M := M -= w - aτ -, the Mach number of the state on the left of the standing wave, and the intermediate states are given by:

τ -= τ L 1 -M L 1 -M , w -= aMτ -, T -= T L , (3.21)
τ + = τ L 1 + M L 1 + νM , w + = νaMτ + , T + = T L , (3.22
)

τ 3 = τ R + τ L M L -νM 1 + νM , w 3 = νaMτ + , T 3 = T R . (3.23)
Besides, there exists a critical value ν in (1, +∞] depending only on the physical data (V L , V R ) and possibly infinite such that the following alternative holds:

• Either ν < ν , and in this case, M belongs to the interval (0, M 0 (ω, ν)] ⊂ (0, min(1, 1/ν)) with

M 0 (ω, ν) = 1 2   1 + ω 2 1 -ω 2 1 + 1 ν - 1 + ω 2 1 -ω 2 2 1 + 1 ν 2 - 4 ν   , (3.24)
where

ω = 1 -M L 1 + M L ∈ (0, 1). ( 3

.25)

The value M = M 0 (ω, ν) gives the unique solution that exactly preserves the energy equality (3.8) across the standing wave, and for 0 < M < M 0 (ω, ν), the energy is dissipated.

• Or ν ≥ ν , and in that case, no positive solution can preserve the energy equality (3.8). The initial data is such that 0 < M L ν < M 0 (ω, ν) where M 0 (ω, ν) is given by (3.24). M must be strictly less than M 0 (ω, ν), and by taking M close enough to

M L ν
we ensure that all the densities remain positive.

In both cases, the mass is conserved across the standing wave [αρw] 0 = 0 and the choice of the value of M determines the dissipation of energy across the standing wave through

[αρwE + απw] 0 = 1 2 (w + aτ L ) 2 Q 0 (M)Ψ(M; ν, ω) ≤ 0, (3.26)
where Q 0 (M) = α L ρ -w -= α R ρ + w + > 0 is the constant mass flux across the standing wave and Ψ is a non-positive function defined by

Ψ(M; ν, ω) = νM -1 νM + 1 -ω 2 M + 1 M -1 , with ω = 1 -M L 1 + M L . (3.27)
Proof . The proof relies on lengthly but easy calculations and therefore, we only sketch it. We look for a solution that satisfies the mass conservation across the standing wave [αρw] 0 = 0. We first focus on energy preserving solutions and we express the jump relation corresponding to the energy conservation equation (3.8) across the standing wave in terms of the left and right states W -and W + , which reads

[αρwE + απw] 0 = -f with f = 0.
This, combined with the mass conservation implies that

τ + 2 (νM) 2 -1 -τ -2 M 2 -1 = 0. (3.28)
In addition, the solution must satisfy equations (3.12) outside a neighborhood of x = 0 which results in a full set of classical Rankine-Hugoniot jump relations. Using these jump relations through the other waves, we can wind up the information to the initial left and right states, showing that (3.28) is equivalent to

Ψ(M; ν, ω) = νM -1 νM + 1 -ω 2 M + 1 M -1 = 0 ⇐⇒ (νM -1)(M -1) -ω 2 (νM + 1)(M + 1) = 0, (3.29)
where the expression of ω is given in (3.27). Then we observe that for the solution to be of signature < 1, 2 >, w + has to be positive and so has to be w -(by the mass conservation). Thus, M = M -must be positive. Moreover, the {w -aτ }-wave must be negative which means that w L -aτ L = w --aτ -< 0, i.e. M L < 1 and M < 1. By (3.28), this implies that M < 1/ν. Consequently, M must be sought in the interval (0, min(1, 1/ν)). Let us now check that (3.29) has a (unique) root in (0, min(1, 1/ν)) if and only if w > 0, and that this root is given by (3. The existence of ν is related to the expression of τ 3 in (3.23) which is the only intermediate specific volume that may be non-positive. It is possible to show that for fixed V L and V R , the function

ν → τ 3 (ν, M 0 (ω, ν)) = τ R + τ L M L -νM 0 (ω, ν) 1 + νM 0 (ω, ν) (3.30)
is a non-increasing function that may become negative for large values of ν. Then, in order to impose the positivity of τ 3 we must no longer exactly conserve the energy at the standing wave (by taking M = M 0 (ω, ν)) but dissipate it by taking M smaller than M 0 (ω, ν). The expression of τ 3 clearly shows that if M is taken close enough to M L ν , we have τ 3 close to τ R which is positive. Finally, the global mass conservation (3.14) is naturally true since the Rankine-Hugoniot relation associated to the mass equation is satisfied for all the waves. Remark 3.3. We can compute explicitly the expression of ν : it is the value of ν which cancels τ 3 in equation (3.30) If we introduce

τ ∞ 3 = lim ν→+∞ τ 3 (ν, M 0 (ω, ν)) = τ R -τ L M L 1 -M L 1 + M L , (3.31)
we can prove that

ν =          +∞ if τ ∞ 3 ≥ 0, M L + τ R τ L 1 - τ R τ L (1 - τ R τ L )(1 + ω 2 ) -(1 -ω 2 )(M L + τ R τ L ) (1 - τ R τ L )(1 -ω 2 ) -(1 + ω 2 )(M L + τ R τ L ) > 1 if τ ∞ 3 < 0. (3.32)
Moreover, for ν ≥ ν , appendix A gives a procedure to choose the value of M and determine the corresponding energy dissipation.

Solutions of signature < 0, 3 >:

We now seek solutions with the supersonic non-resonant ordering of the eigenvalues 0 < w-aτ < w < w + aτ i.e. solutions of signature < 0, 3 >. Again, a one-parameter family of solutions is built, the involved parameter being directly related to the energy dissipation across the standing wave. (3.33)

These solutions can be parametrized by a real parameter θ ∈ (0, 1] measuring the dissipation of energy across the standing wave, and the intermediate states are given by:

τ + = θτ L M 2 L -1 ν 2 M 2 L -1 , w + = νaM L τ + , T + = T L , (3.34)
τ 2 = τ L + τ L 2 (M L -1) 1 -θ (M L + 1)(νM L -1) (M L -1)(νM L + 1) , w 2 = w + + a(τ 2 -τ + ), T 2 = T L , (3.35)
τ 3 = τ R + τ L 2 (M L -1) 1 -θ (M L + 1)(νM L -1) (M L -1)(νM L + 1) , w 3 = w 2 , T 3 = T R .
(3.36)

Besides, there exists a critical value ν ∈ (1, +∞] depending only on the physical data (V L , V R ) and possibly infinite such that the following alternative holds

• Either ν < ν , and in this case the value θ = 1 gives the unique solution that exactly preserves the energy equality (3.8) across the standing wave, and for 0 < θ < 1, the energy is dissipated.

• Or ν ≥ ν , and in that case, no positive solution can preserve the energy equality (3.8). The parameter θ must be strictly less than 1, and by taking θ close enough to 0, we ensure that all the densities remain positive.

In both cases, the mass is conserved across the standing wave [αρw] 0 = 0 and the choice of the value of θ determines the dissipation of energy across the standing wave through

[αρwE + απw] 0 = 1 2 a 2 τ L 2 (M L 2 -1)(θ 2 -1)α L ρ L w L ≤ 0. (3.37)
Proof . Here again, we only sketch the proof. We look for a solution satisfying the mass conservation across the stationary wave [αρw] 0 = 0. The jump relation corresponding to the energy inequality (3.10) leads to nearly the same equation as (3.28):

τ + 2 (νM L ) 2 -1 -τ L 2 M L 2 -1 ≤ 0.
Hence there exists θ in (0, 1) such that

τ + = θτ L M 2 L -1 ν 2 M 2 L -1
, and the value θ = 1 corresponds to the exact preservation of energy. The intermediate states are then computed thanks to the Rankine-Hugoniot jump relations. Eventually, we observe that the functions

ν → τ 2 (ν, θ) = τ L + τ L 2 (M L -1) 1 -θ (M L + 1)(νM L -1) (M L -1)(νM L + 1) , ν → τ 3 (ν, θ) = τ R + τ L 2 (M L -1) 1 -θ (M L + 1)(νM L -1) (M L -1)(νM L + 1) ,
with θ identically equal to 1, are non-increasing functions that may become negative for large values of ν. Then, in order to impose the positivity of τ 2 and τ 3 , we must no longer exactly preserve the energy at the standing wave but dissipate it by taking θ in the interval (0, 1) close enough to 0. Finally, the global mass conservation (3.14) is naturally true since the Rankine-Hugoniot relation associated to the mass equation is satisfied for all the waves.

Remark 3.4. We can compute explicitly the expression of ν . If we introduce

τ ∞ = min lim ν→+∞ τ 2 (ν, θ = 1), lim ν→+∞ τ 3 (ν, θ = 1) = min(τ L , τ R )- τ L 2 (M L -1) M L + 1 M L -1 -1 , (3.38) we can prove that ν =            +∞ if τ ∞ ≥ 0, 1 M L M 2 L -1 + 2 min(τ L ,τ R ) τ L + M L -1 2 M 2 L -1 - 2 min(τ L ,τ R ) τ L + M L -1 2 > 1 if τ ∞ < 0. (3.39)
Moreover, for ν ≥ ν , appendix B describes a procedure to choose the value of θ and determine the corresponding energy dissipation.

Resonant solutions

We now study resonant solutions that are obtained by formally passing to the limit in the nonresonant configuration < 1, 2 > when letting the material wave speed w tend to zero.

Solutions with signature < 1, 1 >:

We first seek solutions with the subsonic non-resonant ordering of the eigenvalues w -aτ < 0 = w < w + aτ i.e. solutions with signature < 1, 1 >. The intermediate states of this solution are given by:

τ -= τ L , w -= 0, T -= T L , (3.41)
τ + = τ R , w + = 0, T + = T R . (3.42)
The mass and energy fluxes are exactly preserved across the standing wave:

[αρw] 0 = [αρwE + απw] 0 = 0. (3.43)
Proof . We assume that w = 0. Let us prove that the intermediate states given by equations (3.41)-(3.42) determine a positive solution of signature < 1, 1 >. For the left-going acoustic wave, we have α = cst = α L and T = cst = T L . Besides, we have w L -aτ L = w -aτ L = -aτ L = -aτ -= w --aτ -. Thus the Rankine-Hugoniot jump relations for (3.12) are clearly satisfied. Similarly, we prove that the jump relations corresponding to the right-going acoustic wave are also satisfied. As for the 0-w wave, the fact that w -= w + = 0, clearly yields

[αρw] 0 = 0, [αρwT ] 0 = 0 and [αρwE + απw] 0 = 0. (3.44)
Here again, the global mass conservation (3.14) is naturally true since the Rankine-Hugoniot relation associated to the mass equation is satisfied for all the waves. Thus Definition 3.1 is satisfied.

Remark 3.5. As mentioned in Section 3.2, the valid jump relations are those of (3.44). In particular, the relaxed pressure π is not constant across the standing wave.

Let us now study resonant solutions that are obtained by formally passing to the limit in nonresonant configurations < 1, 2 > or < 0, 3 > when letting the acoustic wave speed w -aτ tend to zero. We distinguish the case of a divergent section α R > α L i.e. ν < 1 and the case of a convergent section α R < α L i.e. ν > 1.

Acoustic resonance for ν < 1:

We consider initial left and right states W L and W R such that w > 0 and M L < 1 which means that the corresponding solution is of signature < 1, 2 > according to Proposition 3.4. Then, we study the formal limit of the solution as M L goes to 1 -which is equivalent to sending w L -aτ L to 0 -. We expect the {w -aτ }-wave to cross the standing wave and reappear on its right, thus letting the solution to shift from signature < 1, 2 > to signature < 0, 3 >. However, Proposition 3.5 shows that signature < 0, 3 > is possible only if νM L > 1. This implies that in the case of a divergent section ν < 1, there exists a range of values of M L , namely [1, 1 ν ], for which the acoustic wave w -aτ does not appear in the solution. Therefore, we are brought to study the resonant signature < 0, 2 > represented in the figure below. The intermediate states are given by

τ 2 = 2τ L + τ L (M L -1) 1 + νM L , w 2 = νaM L τ 2 , T 2 = T L , (3.46)
τ 3 = τ R + τ L 1 -νM L 1 + νM L + τ L M L -1 1 + νM L , w 3 = w 2 = νaM L τ 2 , T 3 = T R . (3.47)
The mass is conserved across the standing wave [αρw] 0 = 0, while the energy dissipation across this wave is given by

[αρwE + απw] 0 = 1 2 a 2 (2τ L + τ L (M L -1)) 2 νM L -1 νM L + 1 -a 2 τ 2 L (M 2 L -1) α L ρ L w L ≤ 0. (3.

48)

Proof . There are eight unknowns since we have to determine only two intermediate states. Thus we need eight independent jump relations in order to calculate these intermediate states. For the discontinuities located at x t = w 2 and x t = w 3 + aτ 3 , we use the classical Rankine-Hugoniot jump relations associated with system (3.5) which provides us with six independent equations (three for each wave). And for the stationary discontinuity, we use the mass conservation equation (i.e.

[αρw] 0 = 0) as well as the conservation equation of T (i.e. [αρwT ] 0 = 0) which provides us with two more equations. It is then possible to verify that for given W L and W R , there exists a unique solution of signature < 0, 2 > given by equations (3.46)- (3.47). Eventually, we calculate the flux related to the energy equation (3.10) on the standing wave:

[αρwE + απw] 0 = a 2 2 τ 2 2 (νM L ) 2 -1 -τ L 2 M L 2 -1 α L ρ L w L .
This dissipation is wholly determined by W L and W R , and is non-positive since M L ≥ 1 and νM L ≤ 1. Finally, the global mass conservation (3.14) is naturally true since the Rankine-Hugoniot relation associated to the mass equation is satisfied for all the waves.

Acoustic resonance for ν > 1:

Following similar steps as previously, we consider initial left and right states W L and W R such that w > 0 and M L < 1.The corresponding solution is of signature < 1, 2 >, and we study the formal limit of the solution as M L goes to 1 -which is equivalent to sending w L -aτ L to 0 -. It is easy to verify that for M L close to 1, ν = +∞ which means that the Mach number of the state on the left of the standing wave is given by M = M 0 (ω, ν) (see equation (3.24)). Simple calculations show that lim

M L →1 -M 0 (ω, ν) = lim M L →1 - M 0 (ω, ν) = lim ω→0 + M 0 (ω, ν) = 1 ν
as soon as ν > 1. This implies that the specific volume on the left of the standing wave tends to zero:

lim M L →1 -τ -= lim M L →1 - τ L 1 -M L 1 -M 0 (ω, ν) = 0,
which means that the partial mass tends to infinity:

lim M L →1 -α -ρ -= +∞.
However, the Lebesgue measure of the cone supporting this intermediate state tends to zero as M L goes to 1 -: µ (x, t), w L -aτ L < x t < 0 -→ M L →1 -0, and we expect a Dirac measure to appear whose weight is given by lim

M L →1 - 0 w L -aτ L (α -ρ -)(ξ)dξ = lim M L →1 - 0 w --aτ - (α -ρ -)(ξ)dξ = lim M L →1 --(w --aτ -)α L ρ - = lim M L →1 --a(M 0 -1)α L = -a 1 ν -1 α L = -a(α R -α L ) > 0.
Therefore, we are brought to study the resonant signature < 0, 2 > + δ 0 , with little abuse in the notation, depicted in the figure below,

x t 0 w w + aτ WL W2 W3 WR Signature < 0, 2 > + δ0
This non classical solution may be represented by a function that however does not satisfy the mass conservation across the standing wave. The missing mass between the states at the left and right of the standing wave is precisely supported by a Dirac measure on the half-line {x = 0, t > 0} represented in blue in the above figure. Thus Proposition 3.8 shows that there exists a solution in the sense of Definition 3.1 which is a piecewise constant function with non-zero mass flux across the standing wave The intermediate states are given by

[αρw] 0 = a(α R -α L ) = 0.
τ 2 = τ L , w 2 = w = aτ L , T 2 = T L , (3.50)
τ 3 = τ R , w 3 = w = aτ L , T 3 = T R . (3.51)
The mass flux across the standing wave is non-zero:

[αρw] 0 = a(α R -α L ). (3.52)
Moreover, this solution dissipates energy across the standing wave, and the dissipation is completely determined by the initial condition:

[αρwE + απw] 0 = a(α R -α L ) a 2 T 2 L 2 + e(T L ) + p L T L , (3.53) which is negative since α R < α L .
Proof . The intermediate states are obtained by passing to the limit as

M L → 1 -in the expressions (3.21)-(3.23) of the intermediate states of signature < 1, 2 >. Note that when M L → 1 -, we have M L → 1 -and M = M 0 → 1 ν .
Easy manipulations show that the jump relations corresponding to the w and w + aτ waves are satisfied. Indeed As for the jump relations across the stationary wave, one easily checks that

π 2 = π(τ L , T L ) = π = π(τ R , T R ) = π 3 ,
[αρw] 0 = α R w τ L -α L w L τ L = a(α R -α L ) = 0, and T 2 = T L .
It remains to verify that the global mass conservation (3.14) is satisfied in this context where [αρw] 0 = 0. For any T > 0 and R > 0 such that for all x ≥ R, W(-x, T ) = W L and W(x, T ) = W R , we have

R -R (αρ)(T, x)dx = R(α L ρ L ) + T w α R τ L + T (w + aτ R -w ) α R τ R + (R -T (w R + aτ R ))(α R ρ R ), with w = aτ L . Noticing that (αρw) L = aα L , we get R -R (αρ)(T, x)dx -R((αρ) L + (αρ) R ) + T ((αρw) R -(αρw) L ) = T a(α R -α L ) = [αρw] 0 .

Proof of the main result Theorem 3.3

In sections 3.3.2 and 3.3.3, we constructed solutions for w ≥ 0. These solutions correspond to a material wave with positive speed w ≥ 0. We can also construct the symmetric solutions for w < 0 which are denoted by < 2, 1 >, < 3, 0 >, < 2, 0 > and < 2, 0 > + δ 0 with clear notations. Thanks to the Gallilean invariance of system (3.5) (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]), the intermediate states of these symmetric solutions are obtained by exchanging the subscripts L and R and by applying the mapping (α, αρ, αρw, αρT ) → (α, αρ, -αρw, αρT ) to the solutions constructed above. The details are left to the reader. Finally the proof of theorem 3.3 is straightforward. If the constant a is such that τ L > 0 and τ R > 0, then Propositions 3.6, 3.4, 3.5, 3.7 and 3.8 as well as their symmetric counterparts show that, for all positive initial sates W L and W R , there exists a solution in the sense of Definition 3.1. Indeed, the conditions stated in the propositions cover the whole domain of initial conditions Ω r × Ω r .

Numerical approximation

In this section, we use the relaxation approximation defined in section 3 in order to derive a numerical scheme for approximating the entropy weak solutions of the equilibrium system (2.1). We consider a Cauchy problem

∂ t U + ∂ x f(U) + c(U)∂ x U = 0, x ∈ R, t > 0, U(x, 0) = U 0 (x).
For simplicity in the notations, we assume a constant positive time step ∆t and a constant space step ∆x > 0 and we define λ = ∆t ∆x . We introduce a partition of the space R = j∈Z [x j-1 2 , x j+ 1 2 [ where x j+ 1 2 -x j-1 2 = ∆x for all j in Z. We also introduce the discrete intermediate times

t n = n∆t, n ∈ N. The approximate solution at time t n , x ∈ R → U λ (x, t n ) ∈ Ω is a piecewise constant function whose value on each cell C j = [x j-1 2 , x j+ 1 2 [ is a constant value denoted by U n j : U λ (x, t n ) = U n j , for all x in C j , j in Z, in n ∈ N.
Endwise, we denote by x j = 1 2 (x j-1 2 + x j+ 1 2 ) the center of each cell C j . At time t = 0, we use the initial condition U 0 to define the sequence (U 0 j ) j∈Z by

U 0 j = 1 ∆x x j+ 1 2 x j-1 2 U 0 (x)dx, j in Z.

The relaxation method

We now describe the two-step splitting method associated with the relaxation system (3.1) in order to calculate U λ (•, t n+1 ) from U λ (•, t n ). The first step consists in a time-advancing step for the convective part of the relaxation system (3.1), and the second step takes into account the relaxation source term. We first introduce the piecewise constant approximate solution at time t n of system

(3.5) x → W λ (x, t n ) = W n j in C j with W n j =     α n j (αρ) n j (αρw) n j (αρT ) n j     .
At time t = 0, W 0 j is set at equilibrium which means that (αρT ) 0 j = α 0 j . The two steps are defined as follows.

Step 1: Evolution in time (t n → t n+1-)

In the first step, the following Cauchy problem is exactly solved for t ∈ [0, ∆t] with ∆t small enough (see condition (4.2) below)

∂ t W + ∂ x g( W) + d( W)∂ x W = 0, W(x, 0) = W λ (x, t n ). (4.1)
Since x → W λ (x, t n ) is piecewise constant, the exact solution of (4.1) is obtained by gluing together the solutions of the Riemann problems set at each cell interface x j+ 1 2 , provided that these solutions do not interact during the period ∆t, i.e. provided the following classical CFL condition

∆t ∆x max W |σ r i (W)| < 1 2 , i ∈ {0, ..., 3} , (4.2)
for all the W under consideration. More precisely,

If (x, t) ∈ [x j , x j+1 ] × [0, ∆t], then W λ (x, t) = W r x -x j+1/2 t ; W n j , W n j+1 , where (x, t) → W r x t ; W L , W R is the solution of the Riemann problem        ∂ t W + ∂ x g(W) + d(W)∂ x W = 0, W(x, 0) = W L if x < 0, W R if x > 0, (4.3)
constructed in section 3. In order to define a piecewise constant approximate solution at time t n+1-, the solution W λ (x, t) is averaged on each cell C j at time ∆t:

W λ (x, t n+1-) = W n+1- j =     α n+1- j (αρ) n+1- j (αρw) n+1- j (αρT ) n+1- j     = 1 ∆x x j+ 1 2 x j-1 2 W λ (x, ∆t)dx, ∀x ∈ C j , ∀j ∈ Z. (4.4)
Step 2: Instantaneous relaxation (t n+1-→ t n+1 )

In the second step, we solve at time t n + ∆t the ordinary differential equation

∂ t W ε = 1 ε R(W ε ), (4.5)
in the asymptotic regime ε → 0. As an initial condition, we take the function W λ (x, t n+1-) obtained at the end of the first step. Using the definition (3.4) of the relaxation term R, we see that this amounts to imposing T n+1 j := τ n+1 j , thus we have

W n+1 j =     α n+1- j (αρ) n+1- j (αρw) n+1- j α n+1- j     ,
and the new cell value at time t n+1 of the approximate solution U λ (•, t n+1 ) is given by

U n+1 j =   α n+1- j (αρ) n+1- j (αρw) n+1- j   .
This completes the description of the two-step relaxation method.

Remark 4.1 (Choice of the parameter a). In the first step, the solution of the Riemann problem at each interface x j+ 1 2 always exists if the constant a is chosen large enough. As a matter of fact, at each interface, since W L and W R are set to equilibrium, we have T L = τ L and T R = τ R . Thus

τ L = τ L + 1 2a (w R -w L ) - 1 2a 2 (p(τ R ) -p(τ L )) , (4.6)
τ R = τ R + 1 2a (w R -w L ) + 1 2a 2 (p(τ R ) -p(τ L )) , (4.7)
and if a is taken large enough, we ensure that τ L > 0 and τ R > 0 since τ L and τ R are strictly positive. Besides, a can be chosen locally at each interface x j+ 1 2 since the Riemann problems do not interact under the CFL condition (4.2), and it can be chosen so as to avoid any non classical solution with a mass concentration at one of the interfaces, this in order to guarantee the conservativity of the partial mass αρ and the non-linear stability of the method.

Finite volume formulation

It this section, we show that the two-step relaxation method described in the previous section can be written in the form of a non conservative finite volume scheme

U n+1 j = U n j - ∆t ∆x F - j+ 1 2 -F + j-1 2 ,
where F -

j+ 1 2 = F -(U n j , U n j+1 ) and F + j-1 2 = F + (U n j-1 , U n j )
are the left and right numerical fluxes at the cell interfaces x j-1 2 and x j+ 1 2 . Here, the left and right fluxes F -and F + are two distinct functions in order to take into account the non conservative product (see also [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]).

The first step of the relaxation method shows that W λ (x, t) is the exact solution of

∂ t W + ∂ x g(W) + d(W)∂ x W = 0, on R × [t n , t n+1 ] with the initial data W λ (x, t n ) = W n j for all x in C j , with j in Z. Integrating on the rectangle C j × [t n , t n+1 ], we get W n+1- j = W n j - ∆t ∆x g(W r (0 -; W n j , W n j+1 )) -g(W r (0 + ; W n j-1 , W n j )) , (4.8)
since α = α n j is constant on C j × [t n , t n+1 ] so that the product d(W)∂ x W identically vanishes within C j . We then recall that the initial values W n j are set to equilibrium which means that

W n j = (U n j , α n j ), i.e. W n j = M (U n j )
where the mapping M is defined as

M : R 3 -→ R 4
(x, y, z) -→ (x, y, z, x). (4.9)

This mapping, which happens here to be linear, maps U to its so-called maxwellian equilibrium M (U) according to the terminology used in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]. Moreover, the relaxation step shows that U n+1 j = PW n+1- j where P is the linear operator

P : R 4 -→ R 3 (x, y, z, t) -→ (x, y, z). (4.10)
Eventually, when applying operator P to equation (4.8) (note that P • M = Id R 3 ) we obtain the finite volume formulation of our scheme

U n+1 j = U n j - ∆t ∆x F -(U n j , U n j+1 ) -F + (U n j-1 , U n j ) , (4.11)
with

F ± (U L , U R ) = P g W r 0 ± ; M (U L ), M (U R ) . (4.12)
In the sequel, F ± α , F ± αρ and F ± αρw are respectively the first, the second and the third coordinates of the fluxes vectors F -and F + . In practice, it is the finite volume formulation that is used to implement the numerical simulation. Subsequently, we denote by (RS) the relaxation scheme described in sections 4.1-4.2, and whose finite volume form is given by equations (4.11)-(4.12). In the following two sections, we state the main properties of the relaxation scheme.

Basic properties of the scheme

The relaxation approximation method provides a very convenient framework for the L 1 -stability of finite volume methods since the preservation of the phase space Ω by the scheme is almost straightforward. Indeed, the following property states the positivity of the approximated values of the section α n j as well as the the positivity of the partial masses (αρ) n j .

Property 4.1 (L 1 -stability). Under the CFL condition (4.2), the relaxation Finite Volume scheme (RS) preserves positive values for the section and for the density. Indeed, if the initial condition x → U 0 (x) is in Ω, then the values (U n j ) j∈Z,n∈N computed by the scheme are such that, α n j = α 0 j > 0, (αρ) n j > 0, for all j in Z and all n in N, (

that is to say, the section α is preserved throughout time at the level, and the piecewise constant approximate solution U λ (x, t) is also in Ω.

Proof . The first line of equation (4.11) reads α n+1 j = α n j for all j in Z and all n in N. Thus, if α 0 j > 0, this gives the result on α n j . For the positivity of the partial masses (αρ) n j , it is more convenient to consider the two-step splitting formulation of the scheme. The second line of equation (4.4) shows that (αρ) n+1 j is the P 0 projection of the partial mass in the solution W λ (x, ∆t) of the relaxation system. Under the CFL condition (4.2), this solution is obtained by gluing together the Riemann solutions arising from each interface x j+1/2 . Since these solutions are positive according to Theorem 3.3, this concludes the proof.

We also have the following classical consistency property for the relaxation scheme (RS) which guarantees that the constant solutions of system (2.1) are exactly computed. Property 4.2 (Consistency). The relaxation Finite Volume scheme (RS) is consistent in the sense that, for all U in the phase space Ω, the numerical fluxes F -and F + satisfy

F -(U, U) = F + (U, U) = f (U), (4.14)
where f (U), which is defined in (2.3), is the conservative part of the exact flux of the equilibrium system (2.1).

Proof . The proof is almost straightforward, denoting W = M (U), we immediately see that W r (0 ± ; W, W) = W (see equations (3.41) and (3.42)). And P g (W) = f (U) since W = M (U) is at equilibrium.

In addition, under some condition on the choice of the numerical parameter a, the relaxation method is conservative for the mass equation: . If for each interface x j+ 1 2 , the local value of the parameter a = a j+ 1 2 is chosen large enough so that

           ν j+ 1 2 > 1 =⇒ M L,j+ 1 2 = w n j a j+ 1 2 τ n j = 1, ν j+ 1 2 < 1 =⇒ M R,j+ 1 2 = w n j+1 a j+ 1 2 τ n j+1 = -1, (4.15)
then the relaxation scheme (RS) is conservative for the partial mass αρ, in the sense that

F - αρ (U n j , U n j+1 ) = F + αρ (U n j , U n j+1
) for all j in Z. (4.16)

Proof . We have

F + αρ (U n j , U n j+1 ) -F - αρ (U n j , U n j+1 ) = [αρw] 0 W r ( . ; M (U n j ), M (U n j+1 )) . (4.17)
If a j+ 1 2 is chosen as in (4.15), then the solution W r ( . ; M (U n j ), M (U n j+1 )) of the relaxation Riemann problem is a classical solution without mass concentration at x = x j+ 1 2 . Hence [αρw] 0 = 0. Property 4.4 (Well-balanced property). The relaxation scheme (RS) exactly preserves the steady states at rest : w = 0 and ρ = cst. Indeed, if there exists w 0 ∈ R and ρ 0 > 0 such that (αρw) 0 j = 0 and

Non-linear stability

The aim of this section is to exhibit a sufficient condition on the parameter a, which so far is still not determined, that ensures a discrete entropy inequality of the form

(αρE)(U n+1 j ) -(αρE)(U n j ) + ∆t ∆x (G(U n j , U n j+1 ) -G(U n j-1 , U n j )) ≤ 0, (4.19)
which is a discrete counterpart of the energy inequality (2.7) verified by the exact solutions, thus assessing the stability of the method. Here, the numerical entropy flux G(U L , U R ) is to be determined. This can be seen as a stability condition because if one considers the discrete L 1 -norm of the total energy at time t n : j∈Z (αρE)(U n j )∆x, then summing inequality (4.19) over the cells yields

j∈Z (αρE)(U n+1 j )∆x ≤ j∈Z (αρE)(U n j )∆x, for all n in N, (4.20)
which means that the total energy is decreasing in time.

In the first step of the relaxation method, the solution of the Riemann problem (3.5)- (3.11) is computed at each interface x j+ 1 2 , and therefore, one must determine the values of the numerical parameter a. Observe that, under the CFL condition (4.2), the Riemann problems do not interact and the parameter a can be chosen locally at each cell interface x j+ 1 2 . In this section, we prove that if a so-called sub-characteristic condition (also known as Whitham's condition, see [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]) is verified by the parameter a at each cell interface, then the discrete values computed by the relaxation scheme satisfy a discrete entropy inequality.

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

(U L , U R ) if a 2 > -p (τ (ξ)), for almost every ξ in - ∆x 2∆t , ∆x 2∆t . ( 4

.21)

In order to ease the presentation, let us first define some notations for the entropy inequalities of both the equilibrium and the relaxed system. We denote η(U) = αρE the energy of the equilibrium system and F η (U) = αρwE +αwp the associated entropy flux. Hence, the energy inequality satisfied by the entropy weak solutions of (2.1) reads

∂ t η(U) + ∂ x F η (U) ≤ 0. (4.22)
We recall that the mapping η : U → η(U) is convexe on Ω. In the same way we denote η r (W) = αρE the energy of the relaxation system and the associated entropy flux reads F r η (W) = αρwE + αwπ. With these notations, we have the following theorem which proves the non-linear stability of the scheme under Whitham's condition. Theorem 4.5. Assume the CFL condition (4.2) and suppose that for all n ∈ N and j ∈ Z, the parameter a j+ 1 2 satisfies Whitham's condition for (U n j , U n j+1 ). Then the relaxation scheme satisfies the following discrete entropy inequality:

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

(U L , U R ) = F r η (W r (0 + ; M (U L ), M (U R ))).
Remark 4.2. In [START_REF] Saleh | Analyse et approximation numérique par relaxation d'écoulements diphasiques compressibles. Contribution au traitement des phases évanescentes[END_REF], we exhibit a less restrictive condition on the parameter a that still ensures the discrete entropy inequality (4.23) while also handling nearly resonnant solutions for which the specific volumes may be very small. This so-called weak Whitham's condition is based on averaged entropy inequalities instead of a pointwise study of the self-similar Riemann solution.

In order to prove this theorem, let us first prove the following lemma.

Lemma 4.6. With the same notations as in Definition 4.1 and denoting W r (ξ) so as to ease the notation, if a satisfies Whitham's condition for (U L , U R ), then the relaxation approximate Riemann solver satisfies a discrete entropy inequality by interface (see [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]) in the sense that

η( U L ) -η(U L ) + 2∆t ∆x F r η W r (0 + ) -F η (U L ) ≤ 0, (4.24) η( U R ) -η(U R ) + 2∆t ∆x F η (U R ) -F r η W r (0 + ≤ 0, (4.25)
where

U L = 2 ∆x 0 -∆x 2 PW r (x/∆t)dx = 2∆t ∆x 0 -∆x 2∆t PW r (ξ)dξ, (4.26)
U R = 2 ∆x ∆x 2 0 PW r (x/∆t)dx = 2∆t ∆x ∆x 2∆t 0 PW r (ξ)dξ. (4.27)
Proof . We only prove inequality (4.24) (the proof of (4.25) is similar). By Jensen's inequality, the convexity of the map U → η(U) implies that it is sufficient to prove

2∆t ∆x 0 -∆x 2∆t η (PW r (ξ)) dξ -η(U L ) + 2∆t ∆x F r η W r (0 + ) -F η (U L ) ≤ 0, (4.28
∂ t η r (W r ) + ∂ x F r η (W r ) = -f δ 0 , (4.29)
in the weak sense, where f δ 0 = -[αρEw + απw] 0 δ 0 is a positive measure. Integrating this equation over ] -∆x 2 , 0[×]0, ∆t[, and dividing by ∆x 2 , we get

2∆t ∆x 0 -∆x 2∆t η r (W r (ξ))dξ -η r (W L ) + 2∆t ∆x F r η (W r (0 -)) -F r η (W L ) ≤ 0. (4.30) Now, as (W L , W R ) = (M (U L ), M (U R )) are at equilibrium, we have η r (W L ) = η(U L ) and F r η (W L ) = F η (U L ). Moreover, the Riemann solution is constructed such that F r η (W r (0 + )) -F r η (W r (0 -)) ≤ 0. Replacing in (4.30) this yields -η(U L ) + 2∆t ∆x F r η W r (0 + ) -F η (U L ) ≤ - 2∆t ∆x 0 -∆x 2∆t η r (W r (ξ)) dξ. (4.31)
Hence, a sufficient condition for (4.28) (and thus for (4.24)) to hold true is

2∆t ∆x 0 -∆x 2∆t {η (PW r (ξ)) -η r (W r (ξ))} dξ ≤ 0. (4.32)
Now, for almost every ξ in -∆x 2∆t , 0 , we have

η (PW r (ξ))-η r (W r (ξ)) = αρ(ξ) e(τ (ξ)) -e(T (ξ)) - 1 2a 2 π 2 (τ (ξ), T (ξ)) -p 2 (T (ξ)) . (4.33)
Omitting the dependence on ξ, we have: Then, replacing in (4.34) and observing that e = -p we get a sufficient condition for (4.28) (and thus for (4.24)):

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

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

We may now prove theorem 4.5:

Proof . The proof is given in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF], but for the sake of completeness, we reproduce it here. Defining the averages for each half-cell [x j-1 2 , x j ] and [x j , x j+ 1 2 ]:

U R j-1 2 = 2 ∆x xj x j-1 2 PW r x/∆t; M (U n j-1 ), M (U n j ) dx, (4.38)
U L j+ 1 2 = 2 ∆x x j+ 1 2 xj PW r x/∆t; M (U n j ), M (U n j+1 ) dx, (4.39)
we have, under the CFL condition (4.2):

U n+1 j = 1 2 U R j-1 2 + 1 2 U L j+ 1 2
, and as η is convex

η U n+1 j ≤ 1 2 η U R j-1 2 + 1 2 η U L j+ 1 2 . ( 4

.40)

As a j-1 2 satisfies Whitham's condition for (U n j-1 , U n j ), we can apply inequality (4.25) of Lemma 4.6 with U L = U n j-1 and U R = U n j , which yields where the numerical entropy flux is given by G(U L , U R ) = F r η (W r (0 + ; M (U L ), M (U R ))). Finally, observe that G is consistent with the exact entropy flux F η since G(U, U) = F η (U) for all U in Ω. Indeed, for any W in Ω r , we have W r (0 + ; W, W) = W. And if W = M (U) is at equilibrium, we get F r η (W) = F η (U). This concludes the proof of Theorem 4.5.

η U R j-1 2 -η(U n j ) + 2∆t ∆x F η (U n j ) -F r η W r (0 + ; M (U n j-1 ), M (U n j ) ≤ 0. ( 4

Numerical results

In this section, we present two preliminary numerical illustrations in order to assess the good behaviour of the method. An extensive comparison of our method with several schemes of the literature and with exact solutions will be provided in a forthcoming paper [START_REF] Saleh | Thorough numerical comparisons of approximate Riemann solvers for nozzle flows[END_REF]. In particular, important issues such as accuracy, robustness and CPU efficiency will be dealt with.

In both test-cases, the chosen pressure law is an ideal gas pressure law p(τ ) = τ -γ , γ = 3. (4.44)

In the sequel, we denote U = (α, ρ, w) the vector of non-conservative variables.

A subsonic Riemann problem

We consider the following Riemann initial data: the solution is composed of a standing wave separating the states U 1 and U 2 , a left-going σ 1rarefaction wave separating U L and U 1 , and a right-going σ 2 -shock separating U 2 and U R . This test is said to be subsonic in the sence that |w| < c(τ ) everywhere in the solution. Figure 1 displays the relaxation approximation of some classical quantities at the final time T = 1.0 and for two different mesh sizes of 10 2 and 10 3 cells, as well as the exact solution. We observe that the relaxation scheme has a satisfactory behaviour even for the coarse 10 2 -cell mesh.

A Riemann problem with a large α jump

We now consider a Riemann problem with a large jump in the initial cross-section:

U(x, t = 0) = U L = (1., 0.988056834959612, 0.125759712385390), x < 0, U R = (100., 1.01, 0.018403108075689), x > 0. the solution is composed of a standing wave separating the states U 1 and U 2 , a left-going σ 1 -shock separating U L and U 1 , and a right-going σ 2 -rarefaction wave separating U 2 and U R . Such a Riemann solution is actually difficult to approximate in practice because of the very large jump in the initial section. Figure 2 displays the relaxation approximation of some classical quantities at the final time T = 0.15 and for two different mesh sizes of 10 2 and 10 3 cells, as well as the exact solution. We can see that even for the coarse mesh composed with 10 2 cells, the relaxation scheme provides a rather good numerical approximation of the exact solution. For the sake of comparison, the same Riemann solution is computed with Rusanov's scheme, which is known for its robusteness due to its very diffusive feature. The results, provided in Figure 3, show a very unstable behaviour even for a quite refined mesh of 10 3 cells. For the coarse mesh of 10 2 cells, the density values may reach 10 -5 in some cells. Thus, despite its very diffusive nature, Rusanov's scheme is not robust enough for the computation of nozzle flows with strongly discontinuous cross-section. Moreover, to our knowledge, there exists no similar test-cases in the related literature which exemplifies a satisfactory behaviour with such large initial cross-section jumps.

  τL, wL, πL τ L , w , π τ R , w , π τR, wR, πRSelf-similar solution in the case of an initial data with αL = αR.

Theorem 3 . 3 .

 33 Let W L and W R be two positive states in Ω r . Assume that a is such that τ L > 0 and τ R > 0. Then the Riemann problem (3.5)-(3.11) admits a positive solution in the sense of Definition 3.1, whatever the ratio ν

Proposition 3 . 4 .

 34 Let W L and W R be two positive states in Ω r . The Riemann problem (3.5)-(3.11) admits positive solutions in the sense of Definition 3.1 with signature < 1, 2 >, if and only if w > 0 and M L < 1, (3.20)

  [START_REF] Hong | The generic solution of the Riemann problem in a neighborhood of a point of resonance for systems of nonlinear balance laws[END_REF]. Defining ϕ(M) = (νM -1)(M -1) -ω 2 (νM + 1)(M + 1), its first derivative reads ϕ (M) = ν(M -1) + (νM -1) -ω 2 (ν(M + 1) + νM + 1), which is negative on the interval (0, min(1, 1/ν)). In addition, we have ϕ(min(1, 1/ν)) = -2ω 2 (1 + min(ν, 1/ν)) < 0. Hence, by the intermediate value theorem, ϕ has a unique root in (0, min(1, 1/ν)) if and only if ϕ(0) = 1 -ω 2 > 0. From the definition (3.25) of ω, we have ω 2 < 1 ⇔ M L > 0 ⇔ w > 0. The expressions of the intermediate states follow from the Rankine-Hugoniot jump relations. Conversely, if w > 0, then (3.29) has a unique root M 0 in (0, min(1, 1/ν)), and formulas (3.21)-(3.23) give a positive solution of signature < 1, 2 > provided that w L -aτ L < 0 i.e. M L < 1.

Proposition 3 . 5 .

 35 Let W L and W R be two positive states in Ω r . The Riemann problem (3.5)-(3.11) admits positive solutions in the sense of Definition 3.1 with signature < 0, 3 > if and only if w > 0, M L > 1 and νM L > 1.

Signature < 1 , 1 >Proposition 3 . 6 .

 1136 We have the following result: Let W L and W R be two positive states in Ω r . The Riemann problem (3.5)-(3.11) admits a positive solution in the sense of Definition 3.1 with signature < 1, 1 >, if w = 0.(3.40)

2 >Proposition 3 . 7 .

 237 Proposition 3.7 shows that in this acoustic resonance case where ν < 1, one can build a dissipative solution in the sense of Definition 3.1 where the {w -aτ }-wave does not appear. Let W L and W R be two positive states in Ω r . The Riemann problem (3.5)-(3.11) admits a unique positive solution in the sense of Definition 3.1 with signature < 0, 2 >, if w > 0, M L ≥ 1 and νM L ≤ 1.(3.45)

Proposition 3 . 8 .

 38 Let W L and W R be two positive states in Ω r . The Riemann problem (3.5)-(3.11) admits a solution in the sense of Definition 3.1 with signature < 0, 2 > + δ 0 , if ν > 1, w > 0 and M L = 1.(3.49)

  and therefore we have [α, w, π] w = 0 and [α, w + aτ, T ] w+aτ = 0.

Property 4 . 3 (

 43 Conservativity). Denote ν j+1

= ρ 0 ,= ρ 0 = ρ 0

 000 for all j in Z and all n in N.(4.18) Proof . Let us assume that at time t n , (αρw) n j for all j in Z, i.e. w n j = 0 and ρ n j = ρ 0 for all j in Z. At each interface, one has w L = w R = 0 and τ L = τ R = 1/ρ 0 . Hence w = 0 and the solution has the signature < 1, 1 > with w -= w + = 0 and τ L = τ R = 1/(ρ 0 ) i.e. all the intermediate states are at equilibrium (w = 0 and ρ = ρ 0 ). After averaging the solution we obtain for all j in Z. The proof follows from an induction argument.

  ) under Whitham's condition (4.21). The solution W r (ξ) of the Riemann problem (3.5)-(3.11) satisfies

1 0e

 1 .34) A Taylor expansion with integral remainder gives e(τ ) -e(T ) -e (T )(τ -T ) = (T -τ ) 2 (sτ + (1 -s)T )(1 -s)ds. (4.35)

U

  (x, t = 0) = U L = (0.3, 0.206052848877390, -0.003218270138816),x < 0,U R = (0.4, 0.099, -0.015876669673295), x > 0. (4.45) Setting U 1 = (0.3, 0.144092901312860, 0.104099507077253), U 2 = (0.4, 0.15, 0.075), (4.46)

U 1 =

 1 (1., 0.998037207029911, 0.108472909864928), U 2 = (100., 1., 0.0010826),(4.48)

Figure 1 :

 1 Figure 1: Relaxation approximation of the subsonic Riemann problem (4.45) at time T = 1.0, for two mesh sizes of 10 2 and 10 3 cells.

Figure 2 :

 2 Figure 2: Relaxation approximation of the subsonic Riemann problem (4.47) at time T = 0.15, for two mesh sizes of 10 2 and 10 3 cells.

Figure 3 :

 3 Figure 3: Rusanov's scheme approximation of the subsonic Riemann problem (4.47) at time T = 0.15, for two mesh sizes of 10 2 and 10 3 cells.

Acknowledgments. The second author receives a financial support by ANRT through an EDF-CIFRE contract 529/2009. The third author is partially supported by the LRC Manon (Modélisation et Approximation Numérique Orientées pour l'énergie Nucléaire -CEA DM2S/LJLL).

Appendix

A Choice of M for signature < 1, 2 > and corresponding dissipation

When ν < +∞ and ν ≥ ν , M must be chosen in the open interval (0, M 0 (ω, ν)), small enough so as to guarantee the positivity of τ 3 . Being given a fixed real number µ in (0, 1), we may choose M by prescribing τ 3 to a fixed strictly positive value

for every ν ≥ ν c where ν c is the only value of ν that satisfies

and whose expression is

Hence, for ν < ν c , we take M = M 0 (ω, ν) and for ν ≥ ν c , the chosen value of M is obtained by evaluating the inverse function of M → τ 3 at τ 3 = µτ R , which gives

and the corresponding dissipation reads

B Choice of θ for signature < 0, 3 > and corresponding dissipation

When ν < +∞ and ν ≥ ν , the parameter θ must be chosen in the open interval (0, 1) small enough so as to guarantee the positivity of τ 2 and τ 3 . For the sake of clarity, let us assume that τ L ≤ τ R . The case τ L ≥ τ R is straightforward. Being given a fixed real number µ in (0, 1), we may choose θ by prescribing τ 2 to a fixed strictly positive value

for every ν ≥ ν c where ν c is the only value of ν that satisfies

and whose expression is

Hence, for ν < ν c , we take θ = 1 and for ν ≥ ν c , the chosen value of θ is obtained by evaluating the inverse function of θ → τ 2 at τ 2 = µτ L , which gives