domain: spi phys

text:
Claude Marmignon

Fabio Naddei

Florent Renac
email: florent.renac@onera.fr

Energy relaxation approximation for compressible multicomponent flows in thermal nonequilibrium

Keywords: Compressible multicomponent flows, thermal nonequilibrium, relaxation of hyperbolic systems, finite volume method, relaxation scheme Mathematics Subject Classification (2000) 65M12, 65M70, 76T10

teaching and research institutions in France or abroad, or from public or private research centers.

Energy relaxation approximation for compressible multicomponent flows in thermal nonequilibrium

Introduction

We are here interested in finite volume methods to simulate inviscid hypersonic high-temperature flows. Such simulations are of strong significance in many applications (e.g., hypersonic air vehicles [START_REF] Knight | Assessment of cfd capability for prediction of hypersonic shock interactions[END_REF], reentry vehicles [START_REF]Hypersonic and High Temperature Gas Dynamics[END_REF], meteoroid entry into atmosphere [START_REF] Henneton | Numerical simulation of sonic boom from hypersonic meteoroids[END_REF]) and scientific topics (e.g., weakly ionized gases, heat transfer [START_REF] Prakash | High-order shock-fitting methods for direct numerical simulation of hypersonic flow with chemical and thermal nonequilibrium[END_REF], boundary layer stability [START_REF] Knisely | Sound radiation by supersonic unstable modes in hypersonic blunt cone boundary layers. ii. direct numerical simulation[END_REF], shock propagation [START_REF] Honma | Weak spherical shock-wave transitions of n-waves in air with vibrational excitation[END_REF]) related to hypersonic flows. For such flows, effects of thermal and chemical nonequilibria are important and cannot be modeled by the monocomponent compressible Euler equations for a polytropic gas. Real gas models usually include multiple temperatures, chemical reaction rates and vibrational relaxation effects [START_REF] Giovangigli | Multicomponent Flow Modeling[END_REF][START_REF]Hypersonic and High Temperature Gas Dynamics[END_REF]. We here focus on issues related to the numerical treatment of the convective fluxes due to their hyperbolic nature and on the capture of associated features such as strong shocks. We therefore consider thermal nonequilibrium only and neglect chemical nonequilibrium and relaxation of vibration energies that are associated to numerical issues of different nature.

The numerical analysis of hypersonic flows is usually challenging because the characteristic time scales of the chemical reactions and molecular vibrations may be quite different from the characteristic time scale of the flow field. Taking into account the variations in the chemical composition and internal energy modes of a fluid requires to resolve the mass fractions and vibration energies. The thermodynamic properties then depend on these variables which complicates the design of numerical schemes with desirable properties such as robustness (i.e., that keeps positivity of partial densities and internal and vibration energies), stability from a discrete entropy inequality, maximum principle on the mass fractions, etc.

The design of numerical schemes for the approximation of the compressible multicomponent Euler equations has been an active field of research over the past decades. Park proposed an implicit time marching associated to central differencing of ionized flows [START_REF] Park | On convergence of computation of chemically reacting flows[END_REF], while finite volume discretizations have been widely developed with flux splitting techniques [START_REF] Colella | Efficient solution algorithms for the Riemann problem for real gases[END_REF][START_REF] Candler | Computation of weakly ionized hypersonic flows in thermochemical nonequilibrium[END_REF][START_REF] Liu | Nonequilibrium flow computations. i. an analysis of numerical formulations of conservation laws[END_REF][START_REF] Shuen | Inviscid flux-splitting algorithms for real gases with non-equilibrium chemistry[END_REF], Jacobian based methods such as the Roe method [START_REF] Coquel | A Roe-type linearization for the Euler equations for weakly ionized multi-component and multi-temperature gas[END_REF][START_REF] Glaister | An approximate linearised Riemann solver for the three-dimensional Euler equations for real gases using operator splitting[END_REF], the AUSM scheme [START_REF] Gaitonde | An Assessment of CFD for Prediction of 2-D and 3-D High-Speed Flows[END_REF], relaxation based approximate Riemann solvers (ARS) [START_REF] Rouzaud | Development of a Relaxation Scheme for Weakly Ionised Gases[END_REF], etc. High-order extensions have been proposed with the secondorder MUSCL method [START_REF] Druguet | Effects of numerics on Navier-Stokes computations of hypersonic double-cone flows[END_REF], ENO and WENO reconstructions [START_REF] Ton | Improved shock-capturing methods for multicomponent and reacting flows[END_REF][START_REF] Drikakis | Advances in turbulent flow computations using high-resolution methods[END_REF], interface capturing schemes [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach[END_REF][START_REF] Karni | Hybrid multifluid algorithms[END_REF]. Shock fitting techniques have also been addressed in [START_REF] Prakash | High-order shock-fitting methods for direct numerical simulation of hypersonic flow with chemical and thermal nonequilibrium[END_REF]. In this work we will consider the design of explicit first-order finite volume schemes based on ARS. Note that such schemes are also used as building blocks in many other discretizations including high-order finite volume schemes [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF], spectral discontinuous methods such as the discontinuous Galerkin method [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws ii: general framework[END_REF], conservative elementwise flux differencing schemes [START_REF] Fisher | High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains[END_REF], etc.

To ensure entropy stability and robustness when using ARS such as the HLL [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF], Roe [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF], Rusanov [START_REF] Rusanov | Calculation of interaction of non-steady shock waves with obstacles[END_REF], relaxation [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF][START_REF] Coquel | Relaxation of fluid systems[END_REF] schemes, etc., one needs an estimation from above of the maximum wave speeds in the Riemann problem. However, fast estimates such as the two-rarefaction approximation [51, Ch. 9], the iterative algorithm from [START_REF] Guermond | Fast estimation from above of the maximum wave speed in the Riemann problem for the Euler equations[END_REF], or the one based on eigenvalues of the Roe linearisation [START_REF] Einfeldt | On Godunov-type methods near low densities[END_REF] will require time-consuming Newton-Raphson iterations when the equation of state (EOS) differs in the left and right states due to different species compositions. In [START_REF] Dellacherie | Relaxation schemes for the multicomponent Euler system[END_REF] a relaxation technique is applied to the multicomponent Euler system which allows the use of monocomponent schemes for each component and associated EOS and the scheme inherits properties from the monocomponent scheme.

However, this technique requires to compute as many monocomponent schemes as there are species which can become time consuming. Moreover, the entropy of the relaxation system is proved to be convex for constant mass fractions only which is valid for isolated shocks, but fails for interactions of shocks with material interfaces. Here we consider the energy relaxation technique introduced in [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF] for the approximation of the monocomponent compressible Euler equations with a general EOS. In this method, one considers a decomposition of the internal energy including the energy for a polytropic gas thus relaxing the general EOS. The method then allows the design of numerical schemes by using classical numerical fluxes for polytropic gases coupled to instantaneous relaxation of the energy.

In this work, we extend this method to our model and show how to define a numerical scheme from a scheme for the polytropic gas dynamics through a simple formula (equation (4.23)) which corresponds to a splitting of hyperbolic and relaxation operators. In the limit of instantaneous relaxation we show that the solution of the energy relaxation approximation formally converges to a unique and stable equilibrium solution to the multicomponent Euler equations which justifies the splitting. By defining the adiabatic exponent of the polytropic gas as an upper bound of the possible values of adiabatic exponent of the mixture, the scheme for the multicomponent system inherits the properties of the scheme for the monocomponent system: discrete entropy inequality, positivity of the partial densities and internal energies, discrete maximum principle on the mass fractions, and discrete minimum principle on the entropy. An attempt to apply the energy relaxation approximation to the multicomponent Euler system for a fluid mixture in thermal equilibrium has been made in [START_REF] Renac | Entropy stable, robust and high-order DGSEM for the compressible multicomponent Euler equations[END_REF] for the design of a high-order discontinuous Galerkin scheme. However, this work did not provide a general framework to build numerical schemes. The closure laws for the fluid mixture indeed prevent the derivation of a strictly convex entropy for the relaxation system which in turn prevents to apply stability theorems [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Dellacherie | Relaxation schemes for the multicomponent Euler system[END_REF] to the relaxation process. As a consequence the well-posedness of the energy relaxation approximation has not been investigated either. On the other hand, the present work successfully addresses this property and may use any polytropic monocomponent numerical fluxes.

The paper is organized as follows. Section 2 presents the multicomponent compressible Euler system in thermal nonequilibrium. The unstructured finite volume scheme and associated three-point scheme are described in section 3. We introduce and analyze the relaxation in energy approximation in section 4 and use it in section 5 to derive three numerical fluxes for the finite volume scheme. These three schemes are then assessed by numerical experiments in one and two space dimensions in section 6 and concluding remarks about this work are given in section 7.

Model problem

Governing equations and thermodynamic model

Let Ω ⊂ R d be a bounded domain in d space dimensions, we consider the multispecies and multi-temperature model for flows in thermal nonequilibrium [START_REF] Park | Nonequilibrium Hypersonic Aerothermodynamics[END_REF]. Let the IBVP described by the multicomponent compressible Euler system for a mixture of n s species

∂ t u + ∇ • f (u) = 0, in Ω × (0, ∞), (2.1a) u(•, 0) = u 0 (•), in Ω, (2.1b)
with some boundary conditions to be prescribed on ∂Ω (see section 5.4). Here

u =     ρ ρv ρE ρe v     , f (u) =     ρv ρvv + pI (ρE + p)v ρe v v     , (2.2)
denote the conserved variables and the convective fluxes with ρ = (ρ 1 , . . . , ρ n s ) the vector of densities of the n s species, while ρ, v in R d , and E denote the density, velocity vector, and total specific energy of the mixture, respectively. The n d ≤ n s first species are diatomic molecules that are in thermal nonequilibrium and ρe v = (ρ 1 e v 1 , . . . , ρ n d e v n d ) denotes the vector of the associated partial vibration energies. Taking translation-vibration into account and neglecting rotationvibration coupling and anharmonic contributions, each partial vibration energy is linked to the associated vibration temperature T v β through

e v β (T v β ) = r β ϑ v β exp ϑ v β T v β -1 , r β = R M β , 1 ≤ β ≤ n d , (2.3)
where ϑ v β is the characteristic harmonic oscillator temperature, r β is the gas constant of the βth species, R is the universal gas constant, and M β is the molecular weight of the species.

The mixture density, pressure and vibration energy are defined from quantities of the individual species through

ρ = n s α=1 ρ α = ρ n s α=1 Y α , p = n s α=1 p α , ρe v = n d β=1 ρ β e v β , (2.4)
where Y α = ρ α ρ denotes the mass fraction of the αth species, so we have

n s α=1 Y α = 1. (2.5)
The specific total energy of the mixture reads

E = h 0 + e t + e v + e c , h 0 = n s α=1 Y α h 0 α , e t = n s α=1 Y α e t α , e c = 1 2 v • v, (2.6)
where h 0 α ≥ 0 is the enthalpy of formation of species α, e t α = C t v α T denotes the internal translation-rotation energy with C t v α = 3 2 r α for a monoatomic species and C t v α = 5 2 r α for diatomic molecules. The EOS for the mixture pressure in (2.4) is given by the Dalton's law and the partial pressures are assumed to obey polytropic ideal gas EOSs, so

p = n s α=1 ρ α r α T = ρr(Y)T, r(Y) = Z(Y)R, Z(Y) = n s α=1 Y α M α , (2.7)
where Y = (Y 1 , . . . , Y n s ) . Note that the pressure may also be written as

p(Y, ρ, e t ) = γ(Y) -1 ρe t , (2.8)
with

γ(Y) = r(Y) C v t (Y) + 1, r(Y) (2.7) = n s α=1 Y α r α , C v t (Y) = n s α=1 Y α C t v α .
(2.9)

Recalling that either

C t v α = 3 2 r α , or C t v α = 5 2 r α , we have 2 5 C v t (Y) ≤ r(Y) ≤ 2 3 C v t (Y)
by linearity, which, from (2.9), induces the following bounds on γ(Y):

7 5 ≤ γ(Y) ≤ 5 3 ∀0 ≤ Y 1≤α≤n s ≤ 1.
(2.10) System (2.1a) is hyperbolic in the direction n in R d over the set of states [START_REF] Liu | Nonequilibrium flow computations. i. an analysis of numerical formulations of conservation laws[END_REF]] Finally note that we are assuming in (2.11) that the partial densities are positive which would prevent vanishing phases: ρ α = 0 for some α. When such situation occurs the partial velocities, pressure and energies of the species also vanish and this is equivalent to removing the species in the model (2.1) so ρ α > 0 in (2.11) is justified and do not exclude vanishing phases.

Ω a = {u ∈ R n s +n d +d+1 : ρ 1≤α≤n s > 0, v ∈ R d , e t > 0, e v 1≤β≤n d > 0}, (2.11) with eigenvalues λ 1 = v • n -c ≤ λ 2 = • • • = λ n s +n d +d = v • n ≤ λ n s +n d +d+1 = v • n + c,

Entropy pair

Solutions to (2.1) should satisfy an entropy inequality

∂ t η(u) + ∇ • q(u) ≤ 0 (2.13)
for some entropy -entropy flux pair (η, q) with η(•) a strictly convex function and η (u) f i (u) = q i (u) for 1 ≤ i ≤ d, where η (u) = ∂ u η(u) is a column vector with entries η (u) i = ∂ u i η(u) (see (2.18)). In this section we recall the entropy pair for (2.1) derived in [START_REF] Flament | Entropy and entropy production in thermal and chemical non-equilibrium flows[END_REF] and then prove convexity of the entropy which extends the result from [START_REF] Gouasmi | A minimum entropy principle in the compressible multicomponent Euler equations[END_REF] that proved convexity of the entropy for the system in thermal equilibrium, i.e., neglecting the vibrational energy model in (2.1) and (2.2). Following [START_REF] Flament | Entropy and entropy production in thermal and chemical non-equilibrium flows[END_REF], the entropy for a mixture with internal degrees of freedom in nonequilibrium is the sum of associated entropies defined by their differential forms

Tds t α = de t α + p α dτ α ∀1 ≤ α ≤ n s , (2.14a)
T v β ds v β = de v β ∀1 ≤ β ≤ n d , (2.14b)
with τ α = 1 ρ α the covolume of the species α. The entropy pair in (2.13) reads

η(u) = -ρs(u), q(u) = -ρs(u)v, s ≡ n s α=1 Y α s t α + n d β=1 Y β s v β .
(2.15)

Neglecting rotation-vibration coupling and anharmonic contributions, the specific entropies read [START_REF] Flament | Entropy and entropy production in thermal and chemical non-equilibrium flows[END_REF] (up to some additive constants)

s t α (τ α , e t α ) = C t v α ln(e t α ) + r α ln(τ α ) ∀1 ≤ α ≤ n s , (2.16a)
s v β (e v β ) = r β ln(e v β ) + r β 1 + e v β r β ϑ v β ln 1 + r β ϑ v β e v β ∀1 ≤ β ≤ n d . (2.16b)
Note that for smooth solutions, manipulations of (2.1) together with (2.14) show that these entropies satisfy the following conservation laws

∂ t n s α=1 ρ α s t α + ∇ • n s α=1 ρ α s t α v = 0, ∂ t (ρ β s v β ) + ∇ • (ρ β s v β v) = 0 ∀1 ≤ β ≤ n d .
Proposition 2.1 The entropy in (2.15) is a strictly convex and twice differentiable function of u in Ω a .

Proof Twice differentiability is straightforward from (2.16). To prove the convexity we use the trick introduced in [START_REF] Harten | Convex entropies and hyperbolicity for general Euler equations[END_REF] and also used in [START_REF] Gouasmi | A minimum entropy principle in the compressible multicomponent Euler equations[END_REF] to prove that the Hessian of the entropy H η is congruent to the following strictly convex diagonal matrix:

∂u ∂Z H η ∂u ∂Z = ∂u ∂Z ∂η (u) ∂Z = diag (r α τ α ) 1≤α≤n s , (ρθ) 1≤i≤d , ρC v t (Y)θ 2 , -ρ β s v β (e v β ) 1≤β≤n d , (2.17)
where

θ = 1 T , s v β (e v β ) = - r β e v
β (e v β +r β ϑ v β ) < 0 from (2.16b), and Z(u) = (ρ , v , T, e v ) denotes a one-to-one change of variables. Indeed, with some slight abuse in the notation we have

∂u ∂Z =                      1
0 1 0 0 0 0 0 v • • • v • • • v ρI d 0 0 0 0 ρE 1 • • • ρE n d • • • ρE n s ρv ρC v t (Y) ρ 1 • • • ρ n d e v 1 0 • • • 0 0 0 ρ 1 0 . . . . . . . . . . . . . . . . . . 0 e v n d 0 0 0 0 ρ n d                     
, where I d is the identity matrix of size d and

ρE α = ∂ ρ α ρE = C t v α T+h 0 α +ψ α e v α +e c with ψ α = 1 if 1 ≤ α ≤ n d and ψ α = 0 if n d < α ≤ n s . So det(∂ Z u) = ρ d+1 C v t (Y)Π n d β=1 ρ β > 0. Let g t
α = e t α + p α τ α -Ts t α be the Gibbs free energy of the αth species and

g v β = e v β -T v β s v β ,
= d(ρE) -d(ρe v ) -v • d(ρv) - n s α=1 (g t α + h 0 α -e c )dρ α , n d β=1 d(ρ β s v β ) = n d β=1 s v β (e v β ) -e v β θ v β dρ β + θ v β d(ρ β e v β ), with θ v β = 1 T v β = s v β (e v β )
, so the entropy variables read

η (u) := ∂η(u) ∂u =               C t v 1 + r 1 -s t 1 + ψ 1 θ v 1 g v 1 + (h 0 1 -e c )θ . . . C t v ns + r n s -s t n s + ψ n s θ v n s g v n s + (h 0 n s -e c )θ θv -θ θ -θ v 1 . . . θ -θ v n d               , (2.18)
and we obtain for ∂η (u)

∂Z                       r 1 ρ 1 0 -θv -C t v 1 θ -(h 0 1 -e c )θ 2 e v 1 s v 1 (e v
ρ ns -θv -C t v ns θ -(h 0 n s -e c )θ 2 0 • • • 0 0 • • • • • • • • • 0 θI d -θ 2 v 0 • • • 0 0 • • • • • • • • • 0 0 θ 2 0 • • • 0 0 • • • 0 • • • 0 0 -θ 2 -s v 1 (e v 1 ) 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 • • • 0 • • • 0 0 -θ 2 0 -s v n d (e v n d )                      
, so it may be easily checked that (2.17) holds true.

Finally, let τ = 1 ρ be the covolume of the mixture. Using [START_REF] Druguet | Effects of numerics on Navier-Stokes computations of hypersonic double-cone flows[END_REF]) the entropy of the mixture in (2.15) may be written as

e t α C t vα = e t C v t (Y) = T and Y α τ α = τ , for all α, in (2.
s(Y, τ, e t , e v ) = n s α=1 Y α C t v α ln C t v α C v t (Y) e t + Y α r α ln τ Y α + s v (Y, e v ) = C v t (Y) ln e t + r(Y) ln τ + K(Y) + s v (Y, e v ), (2.19a)
K(Y) = n s α=1 Y α C t v α ln C t vα C v t (Y) -r α ln Y α , (2.19b)
s v (Y, e v ) = n d β=1 Y β s v β (e v β ).
(2.19c)

Finite volume method

We consider finite volume schemes for unstructured meshes Ω h ⊂ R d of the form

U n+1 κ -U n κ + ∆t (n) |κ| e∈∂κ |e|h(U n κ , U n κ + e , n e ) = 0 ∀κ ∈ Ω h , n ≥ 0, (3.1)
for the discretization of (2.1a). Here U n+1 κ approximates the averaged solution in the cell κ at time t (n+1) = t (n) + ∆t (n) , ∆t (n) > 0 is the time step, n e is the unit outward normal vector on the face e in ∂κ, and κ + e the neighboring cell sharing the interface e (see fig. 3.1). We assume that each element is shape-regular in the sense of [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]: the ratio of the radius of the largest inscribed ball to the diameter is bounded by below by a positive constant independent of the mesh. The initial condition for (3.1) reads Designing schemes of the form (3.1) on general unstructured grids may be a difficult task and one way to overcome this relies on the use of more simple schemes as building blocks from which we deduce sufficient conditions to inherit their properties. Motivated by the fact that the numerical flux in (3.1) usually solves one-dimensional problems in the direction n e at interfaces e, it is convenient to also consider three-point numerical schemes of the form

U 0 κ = 1 |κ| κ u 0 (x)dV ∀κ ∈ Ω h .
U n+1 j -U n j + ∆t (n) ∆x h(U n j , U n j+1 , n) -h(U n j-1 , U n j , n) = 0, (3.2)
where U n j approximates the averaged solution in the jth cell at time t (n) , ∆x is the space step, and n ∈ R d is a unit arbitrary vector. In particular we are looking for schemes (3.2) that have the following properties under a given condition on the time step

∆t (n) ∆x max j∈Z |λ(U n j )| ≤ 1 2 , (3.3)
where |λ(•)| corresponds to the maximum absolute value of the wave speeds (and will be defined in section 5): the scheme is (i) consistent with (2.1a) and conservative which requires the numerical flux to be consistent:

h(u, u, n) = f (u) • n ∀u ∈ Ω a , (3.4)
and conservative:

h(u -, u + , n) = -h(u + , u -, -n) ∀u ± ∈ Ω a ; (3.5)
(ii) uses a Lipschitz continuous numerical flux;

(iii) entropy stable (ES) for the pair (η, q) in (2.13): it satisfies the inequality

η(U n+1 j ) -η(U n j ) + ∆t (n) ∆x Q(U n j , U n j+1 , n) -Q(U n j-1 , U n j , n) ≤ 0, (3.6)
with some conservative and consistent entropy numerical flux

Q(u -, u + , n) = -Q(u + , u -, -n), Q(u, u, n) = q(u)
• n; (iv) robust: the solution remains in the set of states (2.11): U n j∈Z in Ω a implies U n+1 j∈Z in Ω a ; (v) and it satisfies a discrete maximum principle on the mass fractions:

min(Y n α j-1 , Y n α j , Y n α j+1 ) ≤ Y n+1 α j ≤ max(Y n α j-1 , Y n α j , Y n α j+1 ) ∀1 ≤ α ≤ n s , (3.7)
(vi) together with a minimum principle on the specific entropy in (2.15) [START_REF] Tadmor | A minimum entropy principle in the gas dynamics equations[END_REF][START_REF] Gouasmi | A minimum entropy principle in the compressible multicomponent Euler equations[END_REF]:

s(U n+1 j ) ≥ min s(U n j-1 ), s(U n j ), s(U n j+1 ) . (3.8)
Then it is a classical matter (see e.g. [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] and references therein) that the finite volume scheme (3.1) with the same numerical flux also enjoys properties (i) to (vi) since under the following condition on the time step

∆t (n) max κ∈Ω h |∂κ| |κ| max e∈∂κ |λ(U n κ ± )| ≤ 1 2 , |∂κ| := e∈∂κ |e|, (3.9)
the scheme is a convex combination of updates of three-point schemes (3.2):

U n+1 κ = e∈∂κ |e| |∂κ| U n κ - ∆t (n) |∂κ| |κ| h(U n κ , U n κ + e , n e ) -h(U n κ , U n κ , n e ) , (3.10)
with weights |e| |∂κ| . In particular, the scheme (3.1) satisfies the entropy inequality

η(U n+1 κ ) -η(U n κ ) + ∆t (n) |κ| e∈∂κ |e|Q(U n κ , U n κ + e , n e ) ≤ 0, (3.11)
consistent with (2.13).

In the following we first design two-point numerical fluxes for the one-dimensional scheme (3.2) to satisfy properties (i) to (vi) under the CFL condition (3.3) and then use these numerical fluxes in the finite volume schemes (3.1) with the CFL condition (3.9).

Energy relaxation approximation

In this section we derive a general framework that allows the use of standard numerical schemes for the classical gas dynamics with a polytropic ideal gas EOS. The main results are summarized in theorem 4.3 and show how to build a twopoint numerical flux for (3.1) to satisfy properties (i) to (vi) in section 3 from a numerical flux for the compressible Euler equations with a polytropic EOS.

We here extend the energy relaxation approximation [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF] for the multicomponent Euler system including the vibration energies (section 4.1) and introduce a convex entropy in section 4.2. Section 4.3 is devoted to the analysis of solutions to the relaxation system close to equilibrium. In the limit of instantaneous relaxation, we prove that:

solutions to the relaxation system formally converge to a unique and stable equilibrium solution to the multicomponent Euler equations (2.1a) (theorem 4.1); this equilibrium corresponds to a global minimum of the relaxation entropy which satisfies a minimization principle (lemma 4.2); small perturbations close to the equilibrium are associated to dissipative processes in (2.1a) (theorem 4.2).

These results are then used to infer a numerical scheme for (2.1a) from one for the relaxation system (section 4.4) based on a splitting of the hyperbolic and relaxation operators.

Energy relaxation system

Following the energy relaxation method introduced in [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF], we consider the system

∂ t w + ∇ • g(w ) = - 1 w -M(w ) , (4.1)
and we will denote by (4.1) →∞ the system in homogeneous form, i.e., with → ∞.

For the sake of clarity, we now remove the superscript in w unless needed. Here

w =       ρ ρv ρE p t ρe v ρe r t       , g(w) =       ρv ρvv + p γ (ρ, e p t )I ρE p t + p γ (ρ, e p t ) v ρe v v ρe r t v       , w-M(w) =       0 0 ρ F (Y, e p t ) -e r t 0 ρ e r t -F (Y, e p t )      
, with > 0 the characteristic relaxation time scale, and

p γ (ρ, e p t ) = (γ -1)ρe p t , e p t = E p t -e c , (4.2)
where exponents p and r stand for polytropic and relaxation, respectively. Solutions to (4.1) satisfy the additional conservation law

∂ t ρ + ∇ • (ρ v ) = 0, (4.3)
for the mixture density so the variables ρ, ρv and ρE p t are uncoupled from Y, since n s α=1 Y α = 1, and from e v by (4.2), but are coupled to e r t through the relaxation source terms. This is an important aspect of the model (4.1) and also allows to interpret E p t , e p t , and p γ as the energies and pressure of a polytropic EOS. From (2.10) we set γ as

γ > max 0≤Y 1≤α≤ns ≤1 γ(Y) = 5 3 , (4.4)
which constitutes the subcharacteristic condition for (4.1) to relax to an equilibrium as ↓ 0 [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF] as will be clarified in the results below. System (4.1) is hyperbolic over the set of states

Ω r = w ∈ R n s +n d +d+2 : ρ 1≤α≤n s > 0, v ∈ R d , e p t > 0, e v 1≤β≤n d > 0, e r t > 0 . (4.5) Let w = lim ↓0 w , in this limit, one formally recovers (2.1a) with u = Lw, w = M(w), f (u) = Lg(P(u)), (4.6)
with the operators L : Ω r → Ω a and P : Ω a → Ω r defined by

Lw = ρ ,
E(Y, e r t ) = γ(Y)-1 γ-γ(Y) e r t , T (Y, s γ , e r t ) = γ-γ(Y) γ(Y)-1 (-s γ ) γ e r t 1 γ-1 , (4.11b
∂ τ ζ = -r(Y) τ , ∂ e p t ζ = -r(Y) (γ-1)e p t , ∂ e r t ζ = γ(Y)-γ γ-1 C v t (Y) e r t . (4.13)
Then, direct manipulations (see proof of theorem 4.1 below) show that for smooth solutions of (4.1) the quantity ζ = ζ(w ) satisfies the relation

∂ t ρ ζ + ∇ • (ρ ζ v ) = - 1 ρ r(Y ) (γ -1)e p t e r t e r t -F (Y , e p t ) 2 ,
and defines a convex entropy as stated in the next lemma. Proof This proof has been moved to appendix A for the sake of readability.

Properties of the relaxation system close to equilibrium

We first prove the following minimization principle which states that the equilibrium (4. with 0 < Y 1≤α≤n s ≤ 1, τ > 0, e v 1≤β≤n d > 0, and the minimum is reached at a unique global equilibrium which is solution to (4.8).

Proof Note that (4.14) corresponds to the minimization of a strictly convex function (see lemma 4.1) in a convex set under the linear constraint e p t + e r t = e t , so we only need to find a local minimum for ζ which satisfies (4.14). Let prove that ς in (4.12) is non negative and vanishes at equilibrium (4.8). We rewrite

ς = C v t (Y) ln f (α, x) with f (α, x) = (1-α)(1+x) x αx 1-α α , x =
α = γ(Y)-1 γ-1
in (0, 1) from (4.4). We have

∂ x f (α, x) = 1-α x 2 (αx + α -1), thus ∂ x f (α, x) < 0 for 0 < x < x min := 1-α α , ∂ x f (α, x) > 0 for x > x min , and ∂ x f (α, x min ) = 0. Since f (α, x min ) = 1, ς vanishes at the global minimum αx min = 1 -α ⇔ γ(Y)-1 γ-1 e r t e p t = 1 -γ(Y)-1 γ-1
which indeed corresponds to the equilibrium (4.8): e r t = F (Y, e p t ). The next result concerns the spatially homogeneous system in (4.1):

∂ t w = - 1 w -M(w ) , (4.15)
and is analogue to the H-theorem for kinetic equations. The result below shows that in the limit of instantaneous relaxation ↓ 0 the solution to (4.1) will converge to the equilibrium (4.8).

Theorem 4.1 The vector of variables u = (ρ , ρv , ρE, ρe v ) , with E = E p t + e r t + h 0 + e v , is an invariant of (4.15) and the entropy ρζ decreases in time and reaches a unique minimum which corresponds to the equilibrium (4.8). This equilibrium is stable in the sense of Lyapunov.

Proof From (4.15), we directly obtain that ρ, ρv, and ρe v are constant so ∂ t Y = 0 and ∂ t e v = 0. Then summing the ρE p t and ρe r t equations gives ∂ t (ρE p t + ρe r t ) = 0 so u in theorem 4.1 is constant.

Then, for smooth solutions of (4.15) we get ) which corresponds to the equilibrium (4.8) which in turn corresponds to the global minimum of ζ from lemma 4.2. We therefore conclude that the system is stable by applying the Lyapunov stability criterion with the Lyapunov function w → ζ(w + w 0 ) -ζ(w 0 ) where w 0 = P(u) corresponds to the equilibrium (4.8) projected onto Ω r with P defined in (4.7). Finally note that the partial energies are given explicitly by e t α = C t v α e t /C v t (Y) which confirms that the equilibrium corresponds to a unique state.

∂ t ζ(Y ,
The last result describes the first-order asymptotic analysis of small perturbations in the relaxation process in the neighborhood of the equilibrium (4.8) by performing a formal Chapman-Enskog expansion [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF]. This result extends [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF]Prop. 2.4] to multicomponent flows and allows to understand the relaxation process close to equilibrium as a viscous perturbation to (2.1a) consistent with (2.1a) when ↓ 0. Theorem 4.2 In the limit ↓ 0, small perturbations in (4.1) to the equilibrium (4.8), of the form w = w 0 + w 1 + • • • , obey the following first order asymptotic expansion in :

∂ t u + ∇ • f (u ) = ∇ • f v (u , ∇u ), f v (u , ∇u ) = µ(Y , ρ , e t )     0 (∇ • v )I d (∇ • v )v 0     ,
with w 0 = P(u), u solution to (2.1a) in the limit ↓ 0, and µ(Y , ρ , e t ) =

(γ-γ(Y ))(γ(Y )-1)
∂ t Y + v • ∇Y = 0, (4.17a)
∂ t e p t + v • ∇e p t + (γ -1)e p t ∇ • v = e r t -F (Y , e p t ) , (4.17b)
∂ t e r t + v • ∇e r t = - e r t -F (Y , e p t ) , (4.17c)
from which we deduce

∂ t F (Y , e p t ) + v • ∇F (Y , e p t ) = γ -γ(Y ) γ(Y ) -1 e r t -F (Y , e p t ) -(γ -1)e p t ∇ • v . (4.18)
Plugging (4.16) into either (4.17b), or (4.17c), the order O( -1 ) imposes

e r 0 t = F (Y , e p 0 t ) = γ -γ(Y ) γ(Y ) -1 e p 0 t ,
while the constraint e t = e p t + e r t = e p 0 t + F (Y , e p 0 t ) in (4.14) gives

e p 0 t = γ(Y ) -1 γ -1 e t , e r 0 t = γ -γ(Y ) γ -1 e t , e p k t + e r k t = 0 ∀k ≥ 1.
Plugging again (4.16) into (4.17c) and (4.18), we obtain at leading order

∂ t e r 0 t + v • ∇e r 0 t = -e r 1 t -F (Y , e p 1 t ) , ∂ t F (Y , e p 0 t ) + v • ∇F (Y , e p 0 t ) = γ -γ(Y ) γ(Y ) -1 e r 1 t -F (Y , e p 1 t ) -(γ -1)e p 0 t ∇ • v ,
and since e r 0 t = F (Y , e p 0 t ), we get

-e r 1 t -F (Y , e p 1 t ) = γ -γ(Y ) γ(Y ) -1 e r 1 t -F (Y , e p 1 t ) -(γ -1)e p 0 t ∇ • v ,
and using the above expressions for e p 0 t and e r 0 t gives

e p 1 t = -e r 1 t = -γ -γ(Y ) γ(Y ) -1 γ -1 2 e t ∇ • v .
Finally, in (4.1) consider the momentum equation and add up the equations for ρE p t , ρe r t , ρe v together with an equation for ρh 0 = α ρ α h 0 α . We then obtain up to order O( )

∂ t ρ v + ∇ • ρ v v + p γ (ρ , e p 0 t + e p 1 t ) = 0, ∂ t ρ E + ∇ • ρ E + p γ (ρ , e p 0 t + e p 1 t ) v = 0,
and we conclude by observing that p γ (ρ , e p 0 t ) = p(Y , ρ , e t ) from (4.9) and by using the expression for e p 1 t .

General framework for the design of two-point numerical fluxes

We now clarify the form of the numerical flux for (2.1a) that we deduce from a numerical flux for (4.1) in homogeneous form. The former flux will satisfy the properties (i) to (vi) in section 3 providing that the latter satisfies similar properties. The three-point scheme for (4.1) →∞ reads

W n+1 j -W n j + ∆t (n) ∆x H(W n j , W n j+1 , n) -H(W n j-1 , W n j , n) = 0, (4.19)
with H(w, w, n) = g(w) • n. We assume that under some CFL condition on the time step (see section 5), (4.19) enjoys the properties (i) to (vi) in section 3. In particular we have

ρζ(W n+1 j ) ≤ ρζ(W n j ) -∆t (n) ∆x Z(W n j , W n j+1 , n) -Z(W n j-1 , W n j , n) , (4.20
min(Y n α j-1 , Y n α j , Y n α j+1 ) ≤ Y n+1 α j ≤ max(Y n α j-1 , Y n α j , Y n α j+1 ) ∀1 ≤ α ≤ n s , (4.21)
and the specific entropy

ζ(W n+1 j ) ≤ max ζ(W n j-1 ), ζ(W n j ), ζ(W n j+1 ) , (4.22)
and is robust, W n≥0 j∈Z ∈ Ω r . If (4.4) holds, the three-point numerical scheme (3.2) with the Lipschitz, consistent and conservative numerical flux

h(u -, u + , n) = LH P(u -), P(u + ), n , (4.23)
h X (u -, u + , n) = H X P(u -), P(u + ), n , X ∈ {ρ, ρv, ρe v },

h ρE (u -, u + , n) =
≤ ρζ(W n+1 j ) (4.20) ≤ ρζ(W n j ) -∆t (n) ∆x Z(W n j , W n j+1 , n) -Z(W n j-1 , W n j , n) = ρζ P(U n j ) -∆t (n) ∆x Z P(U n j ), P(U n j+1 ), n -Z P(U n j-1 ), P(U n j ), n = η(U n j ) -∆t (n) ∆x Q(U n j , U n j+1 , n) -Q(U n j-1 , U n j , n) .
Finally, plugging ζ(W n j ) = -s(U n j ) into (4.22) for all j and using (4.14) we obtain (3.8), while (3.7) holds because (4.21) and the components associated to ρ in (4.23) remain unaffected.

Since the pressure in (4.1) obeys a polytropic ideal gas EOS and the variables (Y, e v , e r t ) are purely transported, one may use many methods for (4.19) such as, e.g., the Godunov [START_REF] Godunov | A difference scheme for numerical computation of discontinuous solutions of equations of fluid dynamics[END_REF], Rusanov [START_REF] Rusanov | Calculation of interaction of non-steady shock waves with obstacles[END_REF], HLL [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF], or Roe [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF] schemes, though the latter method does not guaranty robustness [START_REF] Einfeldt | On Godunov-type methods near low densities[END_REF]. In the next section we will consider some of these schemes.

In the definition of the numerical flux (4.23), the L operator consists in adding up some components of H to build the numerical flux for the total energy, ρE, while the P operators consist in taking data at equilibrium. This last operation is equivalent to applying instantaneous relaxation, i.e., to consider (4.1) →∞ , through a splitting of hyperbolic and relaxation operators in (4.1) [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF]. Note that instantaneous relaxation is here justified by the analysis in section 4.3.

Remark 4.1 We note that theorem 4.3 may be directly applied to multidimensional schemes (3.1) instead of three-point schemes (3.2) for which one may obtain a less restrictive constraint on the time step since, from (3.10), (3.9) is only a sufficient condition for (3.1) to satisfy properties (i) to (vi). This is for instance the case of simple solvers, see [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]Sec. 3] and [3, Sec. 2.3]. In contrast, considering three-point schemes (3.2) with the CFL condition (3.9) may facilitate the analysis and may allow to include more schemes in this framework.

Examples of two-point numerical fluxes

In this section we consider examples of two-point numerical fluxes for the homogeneous energy relaxation system (4.1) →∞ to illustrate theorem 4.3. Such numerical fluxes define numerical fluxes for (2.1a) through (4.23). As already noticed, other numerical fluxes may be used since we use a simple polytropic EOS in (4.1).

We here consider Riemann type solvers with numerical fluxes in (4.19) of the form

H(w -, w + , n) = g W(0; w -, w + , n) • n, (5.1)
where W(•; w L , w R , n) is used to approximate the solution to the Riemann problem (4.1) →∞ with initial data w 0

(x) = w L if x := x • n < 0 and w 0 (x) = w R if x • n > 0.
We then build two-point numerical fluxes for (2.1a) by simply applying (4.23).

The Godunov method

As noticed in [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF] it is possible to apply the exact Riemann solver [START_REF] Godunov | A difference scheme for numerical computation of discontinuous solutions of equations of fluid dynamics[END_REF] for polytropic gas to (4.19) where W corresponds to the exact entropy weak solution to the Riemann problem. Consider the compressible Euler equations

∂ t w + ∇ • g( w) = 0, w =   ρ ρv ρE p t   , g( w) =   ρv ρvv + p γ I (ρE p t + p γ )v   , (5.2)
with p γ defined from (4.2) and γ satisfying (4.4). Any variable ψ in {Y, e v , e r t } is uncoupled from the w variables and is purely transported in (4.1) →∞ . Noting that the intermediate states are (ψ L , ψ L , ψ R , ψ R ), the entropy weak solution is made of the Riemann solution for the Euler equations with variables w and fluxes g( w) plus the states for ψ [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF]Lemma 4.6]. The Godunov method is thus ES and guarantees robustness of (4.19) as well as the minimum and maximum principles (4.21) and (4.22) under some standard CFL condition.

Let ρ, u = v • n and p γ be the solution to the Riemann problem for (5.2) with initial data w0 (

x) = ρ L , u L , p(Y L , ρ L , e t L ) if x := x • n < 0 and w0 (x) = ρ R , u R , p(Y R , ρ R , e t R ) if x
• n > 0 and let u be the velocity in the star region.

Note that p γ X = p X given by (2.8) for X = L, R since data are taken at equilibrium in (4.23). Then the numerical flux for (2.1a) reads

h God (u L , u R , n) =     ρ( L Y L + R Y R )u ρu(un + L v ⊥ L + R v ⊥ R ) + p γ n ρ L E L + R E R u + p γ u ρ( L e v L + R e v R )u     , (5.3)
where

E X = h 0 (Y X ) + e p t + γ-γ(Y X ) γ(Y X )-1 e p t + e v X + v ⊥ X •v ⊥ X +u 2 2 , e p t = p γ (γ-1)ρ , v ⊥ X = v X -(v X • n)n, X = L, R, L = 1 if u >
0 and 0 else, and R = 1 -L . Finally, the condition on the time step is (3.3) with λ(•) the maximum wave speed in the Riemann problem for (5.2).

The HLL numerical flux

The HLL ARS [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] for (4.19) reads

W hll ( x t ; w L , w R , n) =    w L , x t < S L , S R w R -S L w L +g(w L )-g(w R ) S R -S L , S L < x t < S R , w R , S R < x t , (5.4)
and is ES [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] and robust [START_REF] Einfeldt | On Godunov-type methods near low densities[END_REF] under the CFL condition

(3.3), with λ = max(|S L |, |S R |),
providing that S L (resp. S R ) is a lower (resp. upper) bound of the speed of the leftmost (resp. rightmost) wave in the exact Riemann solution. Applying (4.23), the numerical flux for (2.1a) reads

h hll (u L , u R , n) =    f (u L ) • n, x t < S L , S R f (u L )•n-S L f (u R )•n+S L S R (u R -u L ) S R -S L , S L < x t < S R , f (u R ) • n, S R < x t , (5.5)
and we evaluate the wave speeds from the two-rarefaction approximation [51, Ch. 9]:

S L = v L • n -c L , S R = v R • n + c R , and
c X = c γ (ρ X , p X ) 1 + γ + 1 2γ p tr p X -1 + , X = L, R, where (•) + = max(•, 0) denotes the positive part, c γ (ρ, p) = γp/ρ, p X = p(Y X , ρ X , e t X )
given by (2.8), γ satisfying (4.4), and

p tr = c γ (ρ L , p L ) + c γ (ρ R , p R ) + γ-1 2 (v L -v R ) • n c γ (ρ L , p L )p -γ-1 2γ L + c γ (ρ R , p R )p -γ-1 2γ R 2γ γ-1 .
Remark 5.1 The two-rarefaction approximation holds for the compressible Euler equations with a polytropic EOS for an adiabatic exponent 1 < γ ≤ 5 3 [START_REF] Guermond | Fast estimation from above of the maximum wave speed in the Riemann problem for the Euler equations[END_REF]. The strict inequality in (4.4) may thus prevent the bound estimates with this approach. However, the analysis in [START_REF] Guermond | Fast estimation from above of the maximum wave speed in the Riemann problem for the Euler equations[END_REF]Lemma 4.2] shows that this may occur only for moderate shock strengths so the scheme remains ES for strong shocks as expected in practice. For instance, we use γ = 1.01× 5 3 in the numerical experiments of section 6 for which the above estimates are valid when either p tr ≤ p X , or p tr 1.05p X .

Pressure relaxation-based numerical flux

We now consider the numerical flux based on relaxation of pressure [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF]Prop. 2.21]. The ARS for (4.19) reads

W r ( x t ; w L , w R , n) =        w L , x t < S L , w L , S L < x t < u , w R , u < x t < S R , w R , S R < x t , (5.6)
where w

X = (ρ X Y X , ρ X v X , ρ X E p, t,X , ρ X e v X , ρ X e r t,X ) , for X = L, R, and v X = v ⊥ X + u n, X = L, R, u = a L u L + a R u R + p γ,L -p γ,R a L + a R (5.7a) p = a R p γ,L + a L p γ,R + a L a R (u L -u R ) a L + a R , (5.7b
)

1 ρ L = 1 ρ L + u -u L a L , 1 ρ R = 1 ρ R + u R -u a R , (5.7c)
E p, t,L = E p t,L - p u -p γ,L u L a L , E p, t,L = E p t,R - p γ,R u R -p u a R , (5.7d) where v ⊥ X = v X -u X n, u X = v X • n, ρ X = n s α=1 ρ α X , Y X = 1
ρ X ρ X , and p γ,X = p γ (ρ X , e p t,X ) defined by (4.2). The wave speeds in (5.6) are evaluated from S L = u L -a L /ρ L and S R = u R + a R /ρ R where the approximate Lagrangian sound speeds [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF] are defined by

   a L ρ L = c γ (ρ L , p γ,L ) + γ+1 2 p γ,R -p γ,L ρ R c γ (ρ R ,p γ,R ) + u L -u R + a R ρ R = c γ (ρ R , p γ,R ) + γ+1 2 p γ,L -p γ,R a L + u L -u R + , if p γ,R ≥ p γ,L , (5.8a)
   a R ρ R = c γ (ρ R , p γ,R ) + γ+1 2 p γ,L -p γ,R ρ L c γ (ρ L ,p γ,L ) + u L -u R + a L ρ L = c γ (ρ L , p γ,L ) + γ+1 2 p γ,R -p γ,L a R + u L -u R + , else, (5.8b)
with γ defined from (4.4). This numerical scheme is based on a relaxation approximation using evolution equations for a relaxation pressure in place of p γ and for a in (5.7) in place of the Lagrangian sound speed ρc γ (ρ, p γ ). The Riemann solution contains only linearly degenerate fields and (5.6) follows from projection of the initial data onto an equilibrium manifold. We refer to [3, Sec. 2.4] or [START_REF] Coquel | Relaxation of fluid systems[END_REF] for complete introductions and in-depth analyses. In particular, the analysis in [3, Sec. 2.4] proves the ES, robustness and the minimum principle on entropy by reversing the roles of energy conservation and entropy inequality [START_REF] Coquel | Some New Godunov and Relaxation Methods for Two-Phase Flow Problems[END_REF]. This technique also applies to the entropy ρζ(w) and we may consider ρE p t = ρe p t (Y, τ, ζ, e r t , e v ) + 1 2 ρv • v as an entropy for the system defined by conservation laws for (ρ , ρv , ρζ, ρe v , ρe r t ) . Indeed, the convexity of ρE

h r (u L , u R , n) =        f (u L ) • n, x t < S L , f L , S L < x t < u , f R , u < x t < S R , f (u R ) • n, S R < x t , (5.9)
with

f X =     ρ X u Y X ρ X u v X + p n (ρ X E X + p )u ρ X u e v X     , X = L, R, and
u = a L u L + a R u R + p L -p R a L + a R , p = a R p L + a L p R + a L a R (u L -u R ) a L + a R , E L = E L - p u -p L u L a L , E R = E R - p R u R -p u a R , p X = p(Y X , ρ X , e t X ),
the other quantities being defined in (5.7) and the wave speed estimates are defined from (5.8) with p γ = p from (4.23).

Wall boundary conditions

Let consider the case of an impermeability condition, u := v • n = 0, at a wall Γ w ⊂ ∂Ω h which is commonly imposed through the use of a mirror state u + = (ρ, ρv + , ρE, ρe v ) with v + = v -2un. For elements κ adjacent to a wall, we modify (3.1) in the following way

U n+1 κ = U n κ - ∆t (n) |κ| e∈∂κ\Γ w |e|h(U n κ , U n κ + e , n e )+ e∈∂κ∩Γ w |e|h r (0, U n κ , U n + κ , n e ) ,
where h corresponds to one of the above numerical fluxes, h r is pressure relaxationbased flux (5.9), and the exponent + denotes the mirror state. The above scheme still can be written as a convex combination of updates of three-point schemes (3.2) as in (3.10), so the entropy inequality (3.11) holds.

Using the mirror state we have from (5.7) that a L = a R = a, p L = p R = p given by (2.8), Y L = Y R = Y so the left and right states have the same thermodynamics. We thus obtain h r (0, u, u + , n) = (0, 0, p n, 0, 0) with p = p + av • n, and a = γ(Y)ρp + (γ(Y) + 1)ρ(v • n) + . This boundary condition is consistent with the impermeability condition and enforces the pressure through the characteristic associated to the eigenvalue v • n + c(Y, e t ).

Note that from theorem 4.3 the entropy flux vanishes at wall boundary interfaces since by Q(u, u + , n) = Z P(u), P(u + ), n = η(u)v • n evaluated at W r 0; P(u), P(u + ), n for which v • n = u = 0. Assuming either compactly supported solutions, or using ES boundary conditions from [START_REF] Svärd | Entropy-stable schemes for the Euler equations with far-field and wall boundary conditions[END_REF] at far-field boundaries, we end with the following global estimate on the entropy:

κ∈Ω h |κ|η(U n+1 κ ) ≤ κ∈Ω h |κ|η(U n κ ),
and using the strict convexity of the entropy η(u), one may use Dafermos' argument to prove L 2 stability of the solution [START_REF] Dafermos | Hyperbolic Conservation Laws in Continuum Physics. Grundlehren der mathematischen Wissenschaften[END_REF] (see e.g. [48, Th. 2.6]).

Numerical experiments

In this section we present numerical experiments, obtained with the CFD code Aghora developed at ONERA [START_REF] Renac | Aghora: A High-Order DG Solver for Turbulent Flow Simulations[END_REF], on problems in one and two space dimensions in order to illustrate the performance of the schemes derived in this work. We restrict ourselves to first-order finite volume schemes of the form (3.1) and use γ = 1.01 × 5 3 in (5.3) to ensure the inequality in (4.4), while we set γ = 5 3 in (5.8) and increase the wave speed estimates S X=L,R by a factor 1.01 in (5.5) and (5.9). The time step is evaluated through (3.3) with λ defined in sections 5.1 to 5.3. For 2D simulations, we impose the freestream values at supersonic inlets and extrapolate variables at supersonic outlets, while we apply the impermeability boundary condition in section 5.4 at walls. Steady computations are obtained by using local time stepping until the 2 norm of the vector of residuals has decreased by a factor 10 10 . Additional results obtained for a monocomponent perfect gas with an equivalent adiabatic exponent are also reported for the sake of comparison: we use either the Roe solver [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF] with entropy fix [START_REF] Harten | Self adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF] (referred to as ROE-PG), or the HLL solver [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] with the two-rarefaction approximation [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF]Ch. 9] for computing the wave speeds (referred to as HLL-PG).

One-dimensional shock-tube problems

We first consider the convection of a material interface separating air (ρ

L = 3.607655, Y 1,L = 1-Y 2,L = 1, e v 1 ,L = 1.8070291, γ 1 = 1.4) from helium (ρ R = 0.5, Y 1,R = 1 -Y 2,R = 0, e v 1 ,R = 0, γ 2 = 5
3 ) in a flow with pressure p L = p R = 1 and velocity u L = u R = 1. Note that a uniform temperature profile is imposed across the interface to avoid spurious oscillations inherent to discretely conservative methods [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach[END_REF], since we are are here interested in the resolution capabilities of the present approach. Results are shown in fig. 6.1 and highlight convergence of the three schemes with some more smearing of the contact by the HLL scheme as expected. We now consider a shock tube problem adapted from [START_REF] Shuen | Inviscid flux-splitting algorithms for real gases with non-equilibrium chemistry[END_REF] initially separating regions with large pressure and temperature ratios: u L = u R = 0, p L = 100p R = 100bars, and T L = 30T R = 9000K. We consider air in thermal equilibrium with a 5 species model with a uniform composition Y N 2 = 0.7543, Y O 2 = 0.2283, Y N O = 0.01026, Y N = 6.5 × 10 -7 , and Y O = 0.00713. We neglect the enthalpies of formation so the gas is a perfect gas with an equivalent adiabatic exponent γ(Y) = 1.402 and we compare our results to the Roe solver for a perfect gas with an adiabatic exponent of 1.402 (ROE-PG). Results in fig. 6.2 show that all solvers provide similar results and converge to the entropy weak solution. We stress that in spite of the crude assumption γ > 5 3 in the numerical fluxes from section 5 they offer similar accuracy as the Roe solver.

Hypersonic flow over a sphere

We now consider the 2D hypersonic flow over a 1 4 inch diameter sphere with the freestream conditions of Lobb's experiments [START_REF] Lobb | Chapter 26 -experimental measurement of shock detachment distance on spheres fired in air at hypervelocities[END_REF]. The freestream Mach number is in fig. 6.5. As expected, while some differences can be identified for underresolved simulations, the results are almost perfectly overlapping for sufficiently fine resolutions.

M ∞ = u ∞ c ∞ = 15.
We are now interested in comparing results obtained with the different schemes and analysing their convergence under grid refinement. To this end we compare the convergence of shock distance from the sphere in fig. 6.6. We use grids with N = 20×20, 40×40, 80×80 and 160×160 elements for the simulation (see fig. 6.3), while the reference distance x ref is evaluated with the Godunov numerical flux on a fine mesh with N = 320 × 320. The results confirm convergence of the shock position and highlight close values obtained with the three different schemes.

Hypersonic flow over a double cone

We finally consider the 2D hypersonic flow over a double cone with angles 25 and 55 deg. adapted from [START_REF] Druguet | Effects of numerics on Navier-Stokes computations of hypersonic double-cone flows[END_REF][START_REF] Knight | Assessment of cfd capability for prediction of hypersonic shock interactions[END_REF] and made of molecular and atomic nitrogen with mass fractions Y N 2 = 0.99, Y N = 0.01. The freestream Mach number is M ∞ = 11.3 with ρ ∞ = 1.34 × 10 -3 kg/m 3 and T ∞ = 303K. The freestream vibration temperature of the molecular nitrogen is taken at T N 2 ∞ = 3085K. We use a series of five unstructured grids (see fig. 6.3). A symmetry condition is imposed at the bottom boundary. Contours of Mach number obtained with the three different schemes on the second finest mesh are displayed in fig. 6.7. Compared to references [START_REF] Druguet | Effects of numerics on Navier-Stokes computations of hypersonic double-cone flows[END_REF][START_REF] Knight | Assessment of cfd capability for prediction of hypersonic shock interactions[END_REF], we observe a strong overestimation of the distance of the bow shock to the wall due to the absence of chemical reactions. However, the results with the three schemes are in good agreement. As done for the previous configuration, we compare in fig. 6.8 the obtained solution to that corresponding to the use of an equivalent perfect gas with adiabatic exponent γ = 1.4032 (HLL-PG) on a fine grid. Once again a very good agreement is obtained. Finally in fig. 6.9, we display the pressure distribution at the wall obtained with the schemes on the five grids. The first pressure peak corresponds to the reflexion of the separated shock at the wall, while the second peak corresponds to rapid pressure variations due to the geometrical transition between the cones. We observe convergence of the solution as the mesh is refined and a close agreement between results from the three schemes on the finest grids.

Concluding remarks

We introduce a general framework to design finite volume schemes for the compressible multicomponent Euler equations in thermal nonequilibrium. The framework allows to define a numerical scheme for its discretization from a scheme for the discretization of the monocomponent polytropic gas dynamics through a simple linear formula. Moreover, the numerical scheme inherits the properties of the scheme for the polytropic gas dynamics under a subcharacteristic condition on the adiabatic exponent of the polytropic gas.

This framework relies on the extension of the relaxation of energy for the gas dynamics equations [START_REF] Coquel | Relaxation of energy and approximate Riemann solvers for general pressure laws in fluid dynamics[END_REF] to the model under consideration in this work. Three different numerical fluxes are constructed with this framework from the polytropic Godunov exact Riemann solver [START_REF] Godunov | A difference scheme for numerical computation of discontinuous solutions of equations of fluid dynamics[END_REF], the HLL numerical flux [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF], and the pressure-based relaxation solver [START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources[END_REF]. They are assessed through numerical simulations of flows in one and two space dimensions with discontinuous solutions and complex wave interactions. The results highlight robustness, nonlinear stability, convergence of the present method, as well as similar performances of the three schemes.

Other numerical fluxes may be deduced from this framework. We also stress that the numerical fluxes designed in this framework can be used as building blocks in second-order finite volume schemes with slope reconstruction [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF], or in the general framework of conservative elementwise flux differencing schemes [START_REF] Fisher | High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains[END_REF]. Such schemes can be shown to inherit properties of the first-order scheme [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF][START_REF] Ranocha | Comparison of some entropy conservative numerical fluxes for the Euler equations[END_REF][START_REF] Zhang | Stability analysis and a priori error estimates of the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws[END_REF][START_REF] Renac | Entropy stable, robust and high-order DGSEM for the compressible multicomponent Euler equations[END_REF]. Future work will consider the use of discontinuous Galerkin schemes for the discretization of the compressible multicomponent Euler equations in thermal nonequilibrium.

A Convexity of the entropy for the energy relaxation system

The object of this appendix is the proof of lemma 4.1.

A.1 Preliminaries

Without loss of generality we define ns as the one corresponding to one species that satisfies rn s = minα rα. Likewise, the mapping w → (Y, τ, e p t , e r t , ev) is surjective in Ω r , so we may rewrite ζ = ζ(w) as a function of these variables. This latter change of variables is also motivated by the result below. Rather than proving convexity of ρζ(w) in terms of the w variables, it is indeed advantageous to consider ζ(Y, τ, e p t , e r t , ev).

Lemma A.1 Given twice differentiable functions f and g such that f (w) = ρg(Y, τ, e p t , e r t , ev), f is strictly convex iff. g is strictly convex in Ω r .

Proof Convexity being invariant under linear maps, the convexity of f is equivalent to that of f (w) = f 1 (ρ 1 , . . . , ρ ns-1 , ρ, ρv , ρE p t , ρe r t , ρev). Then, it is a classical matter that the convexity of f 1 and f 2 with f 1 (ρ, y) = ρf 2 Using (A.2), we get ∂ k r = r k -rn s ≥ 0 since by assumption rn s = minα rα, so we rewrite

Y k + 1 + 1 γ -1 -1 ( x k ∂ k r) 2 r(Y) - n d k=1 s v k (e v k ) Y k (x v k -e v k x k ) 2 + ( x k (γ -1)∂ k Cv t -x k ∂ k r) 2 γ -γ(Y) (γ -
r k x 2 k Y k = rn s + ∂ k r Y k x 2 k = rn s Y k x 2 k + rn s + l Y l ∂ l r r(Y) ∂ k r Y k x 2 k ,
and hence obtain

Q(x) = rn s ( x k ) 2 Yn s + rn s Y k 1 + ∂ k r r(Y) x 2 k + ns-1 k,l=1 ∂ k r∂ l r r(Y) Y l Y k x 2 k -x k x l = rn s ( x k ) 2 Yn s + rn s Y k 1 + ∂ k r r(Y) x 2 k + 1 2 ns-1 k,l=1 ∂ k r∂ l rY k Y l r(Y) x k Y k - x l Y l 2 ,
which is positive, providing that the x 1≤k<ns are not all zero since ∂ k r ≥ 0, so x H ζ x > 0 and we conclude that ζ(Y, τ, e p t , e r t , ev) is strictly convex.

Fig. 3 . 1 :

 31 Fig. 3.1: Notations for the mesh for d = 2.

Lemma 4 . 1

 41 Under the assumption (4.4), the function ρζ(w) defined by (4.11) is a strictly convex entropy for (4.1) in Ω r .

  p t (Y, τ, ζ, e r t , e v ) is equivalent to the convexity of ζ(Y, τ, e p t , e r t , e v ) since from (4.13) ∂ e p t ζ < 0 [23, chap. 2]. The Bouchut scheme guarantees positivity of ρ and e p t under the CFL condition (3.3) with λ = max(|S L |, |S R |). Positivity of ρ = ρY, e r t and e v then follows by cell-averaging the Riemann solution (5.6). Likewise, the discrete minimum maximum principles (4.21) and (4.22) holds for the same reason. Applying (4.23), the numerical flux for (2.1a) reads

  Fig. 6.1: Convection of a material interface: results obtained at time t = 0.25 with the numerical fluxes (5.9) (REL), (5.5) (HLL), and (5.3) (GOD), on two grids with N = 100 (top) and N = 800 (bottom) elements.

  3 with ρ ∞ = 7.83 × 10 -3 kg/m 3 and T ∞ = 293K. The upstream flow is made of nitrogen and oxigen with Y N 2 = 0.79, Y O 2 = 0.21 which are uniform in the flow domain since we do not consider chemical reactions or molecular relaxation. The freestream vibration temperatures are taken at T ∞ for both species. A symmetry condition is imposed at the bottom boundary.

Figure 6 .

 6 4 displays the contours of Mach number on two different grids with the three different schemes. Neglecting chemical reactions overestimates the shock

Fig. 6 . 2 :

 62 Fig. 6.2: Shock tube filled with air: results at time t = 1.5 × 10 -4 obtained with the numerical fluxes (5.9) (REL), (5.5) (HLL), (5.3) (GOD), and the Roe solver for a monocomponent perfect gas (ROE-PG), on two grids with N = 100 (top) and N = 800 (bottom) elements.

Fig. 6 . 3 :

 63 Fig. 6.3: Examples of meshes for the 2D simulations: hypersonic flows over (a) a sphere with N = 100 elements, (b) a double cone with N = 749 elements.

Fig. 6 . 5 :

 65 Fig. 6.5: Hypersonic flow over Lobb's sphere: Mach number contours obtained with the HLL numerical flux (5.5) (red), and the HLL-PG flux considering an equivalent monocomponent perfect gas (black) on two grids with N = 20 × 20 (left) and 160 × 160 (right) elements.

Fig. 6 . 6 :

 66 Fig. 6.6: Hypersonic flow over Lobb's sphere: convergence of the shock position in the symmetry plane y = 0 under mesh refinement.

Fig. 6 . 7 :

 67 Fig. 6.7: Hypersonic flow over a double cone: Mach number contours obtained with the numerical fluxes (5.9) (REL), (5.5) (HLL), and (5.3) (GOD) on a fine grid with N = 52, 237 elements.

Fig. 6 . 8 :

 68 Fig. 6.8: Hypersonic flow over a double cone: Mach number contours obtained with the numerical flux (5.5) (red), and the HLL-PG flux considering an equivalent perfect gas (black) on a fine grid with N = 52, 237 elements.

Fig. 6 . 9 :

 69 Fig.6.9: Hypersonic flow over a double cone: pressure distribution at the wall obtained with the numerical fluxes (5.9) (REL), (5.5) (HLL), and (5.3) (GOD) on a series of five grids with N elements.

1 α=1

 1 Now, let introduce Y = (Y 1 , . . . , Y ns-1 ) (A.1) the vector of the mass fractions of the ns -1 first species. Note that the change of variables Y = Y(Y) is obviously one-to-one from (2.5) and (4.3) for w in Ω r in (4.5), so we may adopt equivalently the notations Y or Y and write r(Y) = r(Y) = rn s + ns-Yα(rα -rn s ), (A.2) with an equivalent relation for Cv t (Y) = Cv t (Y).

  -r(Y) ln τ -r(Y) γ-1 ln e p t -(γ-γ(Y))Cv t (Y) γ-1 ln e r t + l(Y) -sv(Y, ev), with l(Y) =α YαC t vα ln C t vα -Cv t (Y) ln(γ -1), with Yn s = 1 -α<ns Yα, linear in Y. Introducing the short notations ∂ k r ≡ ∂ Y k r(Y), ∂ k Cv t ≡ ∂ Y k Cv t (Y), and ∂ k γ ≡ ∂ Y k γ(Y),the Hessian H ζ (Y, τ, e p t , e r t , ev) of ζ reads with δ k,l the Kronecker symbol and ψ k = 1 if 1 ≤ k ≤ n d and ψ k = 0 if n d < k < ns, and s v k defined by (2.16b). Unless stated otherwise, the subscripts are in the range 1 ≤ k, l < ns, k corresponding to a row index and l corresponding to a column index. Likewise ∂ k r (2.9) = Cv t (Y)∂ k γ + (γ(Y) -1)∂ k Cv t , (A.4) and ∂

1 ) 2 τ= 2 rp t 2 (r t 2 (+ r k x 2 k

 122222 Cv t (Y) + r(Y) x 1)∂ k Cv t -∂ k r γ -Q(x) + x k ∂ k r -r(Y) xτ τ (Y) + x k ∂ k r -r(Y) xr e γ -1)r(Y) + x k ((γ -1)∂ k Cv t -∂ k r) -(γ -γ(Y))Cv t (Y) xs e γ -γ(Y))(γ -1)Cv t (Y) v k +r k ϑ v k ) < 0,so the four last terms are non-negative, andQ(x) = rn s ( x k ) 2 Yn s Y k -( x k ∂ k r) 2 r(Y) .

  where λ 1 and λ n s +n d +d+1 are associated to genuinely nonlinear fields and λ 2≤i≤n s +n d +d to linearly degenerate fields. The frozen sound speed reads

	c(Y, e t ) = γ(Y) γ(Y) -1 e t .	(2.12)

  with L and P defined in (4.7), is ES for the pair (η, q) in (2.13), satisfies (3.6) with Q(u -, u + , n) = Z P(u -), P(u + ), n , the minimum and maximum principles (3.7) and (3.8), and is robust, U n≥0 j∈Z ∈ Ω a .

					H ρE p t (P(u -), P(u + ), n + H ρe r t (P(u -), P(u + ), n
					n s	n d
					+	h 0 α H ρ α (P(u -), P(u + ), n +	H ρe v β (P(u -), P(u + ), n ,
					α=1	β=1
	Proof By consistency of H: h(u, u, n)	(4.23) = LH P(u), P(u), n = Lg(P(u)) •
	n	(4.6) = f (u) • n, while Lipschitz continuity and conservation are direct.
		Then, let W n j = P(U n j ) so ρζ(W n j ) = η(U n j ) and Z(W n j , W n j+1 , n) = Q(U n j , U n j+1 , n),
	and define W n+1 j	from (4.19) and U n+1 j	= LW n+1 j	. If U n j ∈ Ω a , then W n j =
	P(U n j ) ∈ Ω r since e p t = γ(Y)-1 γ-1 e t and e r t = γ-γ(Y) γ-1 e t , and W n+1 j	∈ Ω r so
	U n+1 j	= LW n+1 j	∈ Ω a by (4.6). Now, by the minimization principle (4.14) we
	have
	j η(U n+1	)	(4.14)

1 ρ

 1 , 1 ρ y are equivalent. Since E p t = e p t -1 2 v • v, the convexity of f is equivalent to the convexity of ρf 2 (Y, τ, v, E p t , e r t , ev) = ρg(Y, τ, e p t , e r t , ev) [23, chap. 2]. A.2 Proof of lemma 4.1 To prove that ρζ(w) is convex it is sufficient to prove that ζ(Y, τ, e p t , e r t , ev) is convex from lemma A.1 and, from (4.12), we rewrite ζ as ζ(Y, τ, e p t , e r t , ev) = Cv t (Y) ln (γ -γ(Y))Cv t (Y) +

  Cv t ∂ l Cv t Cv t (Y) + Cv t (Y)∂ k γ∂ l γ γ -γ(Y) γ(Y) -1 , 1 ≤ k, l < ns. (A.5) We now prove that H ζ is symmetric positive definite. Let x = (x 1≤i<ns , xτ , xr, xs, x v ∂ k Cv t ) 2 Cv t (Y) + Cv t (Y)( x k ∂ k γ) 2 γ -γ(Y) γ(Y) -1 -1)( x k ∂ k Cv t ) 2 γ -γ(Y) Cv t (Y) -2( x k ∂ k r)( x k ∂ k Cv t ) γ -γ(Y) Cv t (Y) (γ -1)∂ k Cv t -x k ∂ k r) 2 γ -γ(Y) (γ -1)Cv t (Y), so using this result together with (A.3) we obtainx H ζ x = rn s ( x k ) 2

	ns-1 k,l=1	x k ∂ 2 kl ζx l	(A.5) = ( x k (A.4) rn s ( x k ) 2 Yn s + r k x 2 k Y k + = rn s ( x k ) 2 Yn s + r k x 2 k Y k ( x k ∂ k Cv t ) 2 + Cv t (Y)
							2
				x k ∂ k r -(γ(Y) -1)∂ k Cv t
			+	γ -γ(Y) r(Y)
			=	rn s ( x k ) 2 Yn s	+	r k x 2 k Y k	+	( x k ∂ k Cv t ) 2 Cv t (Y)
			+ (γ(Y) (A.4) ( x k ∂ k r) 2 γ -γ(Y) r(Y) + = rn s ( x k ) 2 Yn s + r k x 2 k Y k ( x k Yn s + ( x k ∂ k r) 2 (γ -1)r(Y) + + r k x 2 k
	2 kl ζ = in R ns+n d +3 non zero and use the notation rn s Yn s + r k Y k δ k,l + ∂ k 1≤i≤n d ≡ ns-1 k=1 , we get	)

Acknowledgements Part of this work was funded under the Onera research project JER-OBOAM. One of the authors (F. Renac) gratefully acknowledges this support.