D Grapsas

R Herbin

J.-C Latché

Y Nasseri

A staggered pressure correction numerical scheme to compute a travelling reactive interface in a partially premixed mixture

We address a turbulent deflagration model with a flow governed by the compositional Euler equations and the flame propagation represented by the transport of the characteristic function. The numerical scheme works on staggered unstructured, meshes with a time-marching algorithm solving first the chemical species mass balances and then the mass, momentum and energy balances. A pressure correction technique is used for this latter step, which solves a balance equation for the sensible enthalpy with corrective terms to ensure consistency of the total energy. The approximate solutions respect the physical bounds and satisfy a conservative weakly-consistent discrete total energy balance equation. Numerical evidence shows that they converge to the solution of the infinitely fast chemistry continuous problem when the chemical time scale tends to zero with the space and time steps.

Problem position

In this paper, we study a numerical scheme for the computation of large scale turbulent deflagrations occurring in a partially premixed atmosphere. In usual situations, such a physical phenomena is driven by the progress in the atmosphere of a shell-shaped thin zone, where the chemical reaction occurs and which thus separates the burnt area from fresh gases; this zone is called the flame brush. The onset of the chemical reaction is due to the tempera-D. Grapsas, R. Herbin and Y. Nasseri Aix-Marseille Université, CNRS, Centrale Marseille, I2M, UMR 7373, 13453 Marseille, France, e-mail: {dionysios.grapsas, raphaele.herbin,youssouf.nasseri}@univamu.fr J.-C. Latché Institut de Radioprotection et de Sûreté Nucléaire (IRSN), BP 3, 13115 Saint-Paullez-Durance cedex, France e-mail: jean-claude.latche@irsn.fr 1 ture elevation, so the displacement of the flame brush is driven by the heat transfers inside and in the vicinity of this zone. Modelling of deflagrations still remains a challenge, since the flame brush has a very complex structure (sometimes presented as fractal in the literature), due to thermo-convective instabilities or turbulence [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF][START_REF] Peters | Turbulent Combustion. Cambridge Monographs of Mechanics[END_REF]. Whatever the modelling strategy, the problem thus needs a multiscale approach, since the local flame brush structure is out of reach of the computations aimed at simulating the flow dynamics at the observation scale, i.e. the whole reactive atmosphere scale. A possible way to completely circumvent this problem is to perform an explicit computation of the flame brush location, solving a transport-like equation for a characteristic function of the burnt zone; such an approach transfers the modelling difficulty to the evaluation of the flame brush velocity (or, more precisely speaking, to the relative velocity of the flame brush with respect to the fresh gases), by an adequate closure relation, and the resulting model is generally referred to as a Turbulent Flame velocity Closure (TFC) model [START_REF] Zimont | Gas premixed combustion at high turbulence. Turbulent flame closure combustion model[END_REF]. The transport equation for the characteristic function of the burnt zone is called in this context the G-equation, its unknown being denoted by G [START_REF] Peters | Turbulent Combustion. Cambridge Monographs of Mechanics[END_REF]. Such a modelling is implemented in the in-house software P 2 REMICS (for Partially PREMIxed Combustion Solver) developed, on the basis of the software components library CALIF 3 S (for Components Adaptative Library For Fluid Flow Simulations, see [START_REF]CALIF 3 S: A software components library for the computation of reactive turbulent flows[END_REF]) at the French Institut de Radioprotection et de Sûreté Nucléaire (IRSN) for safety evaluation purposes; this is the context of the work presented in this paper.

Usually, TFC models apply to perfectly premixed flows (i.e. flows with constant initial composition), and the chemical state of the flow is governed by the value of G only: G ∈ [0, 1], for G ≥ 0.5, the mixture is supposed to be in its fresh (initial) state and G < 0.5 is supposed to correspond to the burnt state; in both cases, the composition of the gas is known (it is equal to the initial value in the fresh zones, and to the state resulting from a complete chemical reaction in the burnt zone).

However, for partially premixed turbulent flows (i.e. flows with nonconstant initial composition), the situation is more complex, since the composition of the mixture can no more be deduced from the value of G. An extension for this situation, in the inviscid case, is proposed in [START_REF] Beccantini | The reactive Riemann problem for thermally perfect gases at all combustion regimes[END_REF]. The line followed to formulate this model is to write transport equations for the chemical species initially present in the flow, as if no chemical reaction occured, and then to compute the actual composition in the burnt zone (i.e. the part of the physical space where G < 0.5) as the chemical equilibrium composition, thus supposing an infinitely fast reaction. This model is referred to in the following as the "asymptotic model", and is recalled in the first part of Section 2.

We propose here an alternate extension, which consists in keeping the classical reactive formulation of the chemical species mass balance, but evaluating the reaction term as a function of G: it is set to zero in the fresh zone (G ≥ 0.5), and to a finite (but possibly large) value in the burnt zone (G < 0.5). This model is referred to as the "relaxed model"; it is in fact more general, as it may be readily extended to cope with diffusion terms, while the "asymptotic model" cannot (to this purpose, a balance for the actual mass fractions is necessary). We then build a numerical scheme, based on a staggered discretization of the unknowns, for the solution of the relaxed model; this algorithm is of fractional step type, and employs a pressure correction technique for hydrodynamics. The balance energy solved by the scheme is the so-called (non conservative) sensible enthalpy balance, with corrective terms in order to ensure the weak consistency (in the Lax-Wendroff sense) of the scheme. It enjoys the same stability properties as the continuous model: positivity of the density and, thanks to the choice of the enthalpy balance, the internal energy, conservation of the total energy, chemical species mass fractions lying in the interval [0, 1]. In addition, it is shown to be in fact conservative: indeed, its solutions satisfy a discrete conservative total energy balance whose time and space discretization is non-standard, but weakly consistent with its continuous counterpart. This algorithm is an extension to the reactive case of the numerical scheme for compressible Navier-Stokes equations described and tested in [START_REF] Grapsas | An unconditionally stable finite element-finite volume pressure correction scheme for the compressible Navier-Stokes equations[END_REF].

As the reaction term gets stiffer, the relaxed model should boil down to the asymptotic one, for which a closed form of the solution of Riemann problems is available. Numerical tests are performed which show that this is indeed the case. In addition, we observe that the accuracy of the scheme (for this kind of application) is highly dependent on the numerical diffusion introduced by the scheme in the mass balance equation for the chemical species, comparing the results for three approximations of the convection operator in these equations: the standard upwind scheme, a MUSCL-like scheme introduced in [START_REF] Piar | A formally second order cell centered scheme for convection-diffusion equations on general grids[END_REF] and a first order scheme designed to reduce diffusion proposed in [START_REF] Després | Contact discontinuity capturing scheme for linear advection and compressible gas dynamics[END_REF].

The presentation is structured as follows. We first introduce the asymptotic and the relaxed models in Section 2. Then we give an overview of the content of this paper in Section 3, writing the scheme in the time semi-discrete setting and stating its stability and consistency property. The fully discrete setting is given in two steps, first describing the space discretization (Section 4) and then the scheme itself (Section 5). The conservativity of the scheme is shown in Section 6. Finally, numerical experiments are presented in Section 7.

The physical models

We begin with the description of the asymptotic model introduced in [START_REF] Beccantini | The reactive Riemann problem for thermally perfect gases at all combustion regimes[END_REF] and then turn to the relaxed model proposed in the present work.

The asymptotic model -For the sake of simplicity, only four chemical species are supposed to be present in the flow, namely the fuel (denoted by F ), the oxydant (O), the product (P ) of the reaction, and a neutral gas (N ). A one-step irreversible total chemical reaction is considered, which is written:

ν F F + ν O O + N → ν P P + N,
where ν F , ν O and ν P are the molar stoichiometric coefficients of the reaction. We denote by I the set of the subscripts used to refer to the chemical species in the flow, so I = {F, O, N, P } and the set of mass fractions of the chemical species in the flow reads {y i , i ∈ I} (i.e. {y F , y O , y N , y P }). We now define the auxiliary unknowns {ỹ i , i ∈ I} as the result of the (inert) transport by the flow of the initial state, which means that the {ỹ i , i ∈ I} are the solutions to the following system of equation:

∂ t (ρỹ i ) + div(ρỹ i u) = 0, ỹi (x, 0) = y i,0 (x), for i ∈ I, (1)
where ρ stands for the fluid density, u for the velocity, and y i,0 (x) is the initial mass fraction of the chemical species i in the flow. These equations are supposed to be posed over a bounded domain

Ω of R d , d ∈ {1, 2, 3}
and a finite time interval (0, T ). The initial conditions are supposed to verify i∈I y i,0 = 1 everywhere in Ω, and this property is assumed to be valid for any t ∈ (0, T ), which is equivalent with the mixture mass balance, given below. The characteristic function G is supposed to obey the following equation:

∂ t (ρG) + div(ρGu) + ρ u u f |∇G| = 0, (2)
associated to the initial conditions G = 0 at the location where the flame starts and G = 1 elsewhere. The quantity ρ u is a constant density, which, from a physical point of view, stands for a characteristic value for the unburnt gases density. The chemical mass fractions are now computed as:

if G > 0.5,

y i = ỹi for i ∈ I, if G ≤ 0.5, y F = ν F W F z+ , y O = ν O W O z-, y N = ỹN , with z = 1 ν F W F ỹF - 1 ν O W O ỹO . (3)
In these relation, z+ and z-stand for the positive and negative part of z, respectively, i.e. z+ = max(z, 0) and z-= -min(z, 0), and, for i ∈ I, W i is the molar mass of the chemical species i. The physical meaning of Relation (3) is that the chemical reaction is supposed to be infinitely fast, and thus that the flow composition is stuck to the chemical equilibrium composition in the so-called burnt zone, which explains why the model is qualified as "asymptotic". The product mass fraction is given by y

P = 1-(y F +y O +y N ).
The flow is governed by the Euler equations:

∂ t ρ + div(ρu) = 0, (4a) ∂ t (ρu i ) + div(ρu i u) + ∂ i p = 0, i = 1, d, (4b)
∂ t (ρE) + div(ρEu) + div(pu) = 0, (4c)

p = (γ -1) ρe s , E = 1 2 |u| 2 + e, e = e s + i∈I y i ∆h 0 f,i , (4d)
where p stands for the pressure, E for the total energy, e for the internal energy, e s for the so-called sensible internal energy and, for i ∈ I, ∆h 0 f,i is the formation enthalpy of the chemical species i. The equation of state (4d) supposes that the fluid is a perfect mixture of ideal gases, with the same iso-pressure to iso-volume specific heat ratio γ > 1. This set of equations is complemented by homogeneous Neumann boundary conditions for the velocity:

u • n = 0 a.e. on ∂Ω, (5)
where ∂Ω stands for the boundary of Ω and n its outward normal vector.

The "relaxed" model -This model retains the original form of the governing equations for reactive flows: a transport/reaction equation is written for each of the chemical species mass fractions; the value of G controls the reaction rate ω, which is set to zero when G ≥ 0.5, and takes non-zero (and possibly large) values otherwise. The unknowns {y i , i ∈ I} are thus now solution to the following balance equations:

∂ t (ρy i ) + div(ρy i u) = ωi , ỹi (x, 0) = y i,0 (x) for i ∈ I, (6)
where the reactive term ωi is given by: ωi

= 1 ε ζ i ν i W i ω, with ω = η(y F , y O ) (G -0.5) - and η(y F , y O ) = min( y F ν F W F , y O ν O W O
), [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF] with

ζ F = ζ O = -1, ζ P = 1 and ζ N = 0. Note that, since ν F W F + ν O W 0 = ν P W P ,
we have i∈I ωi = 0, which, summing on i ∈ I the species mass balances, allows to recover the equivalence between the mass balance and the fact that i∈I y i = 1. The factor η(y F , y O ) is a cut-off function, which prevents the chemical species mass fractions from taking negative values (and, consequently, values greater than 1, since their sum is equal to 1).

The rest of the model is left unchanged.

3 General description of the scheme and main results

Time semi-discrete algorithm

Instead of the total energy balance equation, the scheme solves a balance equation for the sensible enthalpy h s = e s + p/ρ, which is formally derived as follows. The first step is to establish the kinetic energy balance formally and subtract from (4c) to obtain a balance equation for the internal energy. Thanks to the mass balance equation, for any regular function ψ

∂ t (ρψ) + div(ρψu) = ρ∂ t ψ + ρu • ∇ψ.
Using twice this identity and then the momentum balance equation, we have for 1 ≤ i ≤ d:

1 2 ∂ t (ρu 2 i ) + 1 2 div(ρu 2 i u) = ρ u i ∂ t u i + ρu i u • ∇u i = u i ∂ t (ρu i ) + div(ρu i u) = -u i ∂ i p,
and, summing for i = 1 to d, we obtain the kinetic energy balance:

1 2 ∂ t (ρ|u| 2 ) + 1 2 div(ρ|u| 2 u) = u • ∂ t (ρu) + div(ρu ⊗ u) = -u • ∇p.
Substituting the expression of the total energy in (4c), yields

∂ t (ρe) + div(ρeu) + 1 2 ∂ t (ρ|u| 2 ) + 1 2 div(ρ|u| 2 ) + u • ∇p + pdiv(u) = 0,
which, using the kinetic energy balance, gives the total internal energy balance:

∂ t (ρe) + div(ρeu) + pdiv(u) = 0. ( 8
)
Using the linearity of the mass balance of the chemical species i, for any i ∈ I, we derive the reactive energy balance:

∂ t ρ i∈I ∆h 0 f,i y i + div ρ i∈I ∆h 0 f,i y i u = i∈I ∆h 0 f,i ωi = -ωθ . (9)
Subtracting ( 9) from (8) yields the sensible internal energy balance:

∂ t (ρe s ) + div(ρe s u) + pdiv(u) = ωθ . (10)
Finally, using the relation between the sensible energy and the sensible enthalpy, we obtain the sensible enthalpy balance:

∂ t (ρh s ) + div(ρh s u) -∂ t p -u • ∇p = ωθ . (11)
The numerical resolution of the mathematical model is realized by a fractional step algorithm, which implements a pressure correction technique for hydrodynamics in order to separate the resolution of the momentum balance from the other equations of the Euler system. Supposing that the time interval (0, T ) is split in N sub-intervals, of constant length δt = T /N , the semi-discrete algorithm is given by: Reactive step:

G n+1 : 1 δt (ρ n G n+1 -ρ n-1 G n ) + div(ρ n G k u n ) + ρ u u f |∇G k | = 0, (12a) Y n+1 N : 1 δt (ρ n y n+1 N -ρ n-1 y n N ) + div(ρ n y k N u n ) = 0. ( 12b
)
z n+1 : 1 δt (ρ n z n+1 -ρ n-1 z n ) + div(ρ n z k u n ) = 0. ( 12c
)
Y n+1 F : 1 δt (ρ n y n+1 F -ρ n-1 y n F ) + div(ρ n y k F u n ) = - 1 ε ν F W F ω(y n+1 F , z n+1 ), (12d)
Y n+1 P : y n+1 F + y n+1 O + y n+1 N + y n+1 P = 1. (12e) (12f)
Euler step: ũn+1 :

1 δt (ρ n ũn+1 i -ρ n-1 u n i ) + div(ρ n ũn+1 i u n ) + ρ n ρ n-1 1/2 ∂ i p n = 0, i = 1, . . . , d, (12g)
u n+1 ρ n+1 h n+1 s p n+1 : 1 δt ρ n (u n+1 i -ũn+1 i ) + ∂ i p n+1 - ρ n ρ n-1 1/2 ∂ i p n = 0, i = 1, . . . , d, 1 δt (ρ n+1 -ρ n ) + div(ρ n+1 u n+1 ) = 0, 1 δt (ρ n+1 h n+1 s -ρ n h n s ) + div(ρ n+1 h n+1 s u n+1 ) - 1 δt (p n+1 -p n ) -u n+1 • ∇p n+1 = ωn+1 θ + S n+1 , p n+1 = γ -1 γ ρ n+1 h n+1 s . (12h)
Equations (12a)-(12h) are solved successively, and the unknown for each equation is specified before each equation. In the convection term of the equations of the reactive step, the index k may take the value n (so the scheme is explicit) or n + 1 (so the scheme is implicit). The unknown z is an affine combination of y F and y O , defined so that the reactive term cancels:

z = 1 ν F W F y F - 1 ν O W O y O . (13)
Thus the value of y n+1 O is deduced from y n+1 F and z n+1 , which allows to express ω in (12d) as a function of y n+1 F and z n+1 , instead of y n+1 F and y n+1 O as suggested by Relation [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]. In addition, we have:

η(y n+1 F , y n+1 O ) = min( y n+1 F ν F W F , y n+1 O ν O W O ) = 1 ν F W F y n+1 F if z n+1 ≤ 0, 1 ν O W O y n+1 O = 1 ν F W F y n+1 F -z n+1 otherwise.
Hence, because of the specific form of the function η, the right hand side of (12d) boils down to an affine term, even if η vanishes when y F or y O vanishes, and the scheme is fully implicit in time with respect to the reaction term. This is the motivation for the choice of the form of η. It is fundamental to remark that Equations (12b)-(12e) are equivalent to the following system:

1 δt (ρ n y n+1 i -ρ n-1 y n i ) + div(ρ n y k i u n ) = 1 ε ζ i ν i W i ω(y n+1 F , y n+1 O ), i ∈ I, (14)
where we recall that

ζ F = ζ O = -1, ζ P = 1 and ζ N = 0.
Indeed, dividing the fuel mass balance equation (12d) by ν F W F , substracting Equation (12c) and finally multiplying by ν O W O yields the desired mass balance equation for the oxydant chemical species. Finally, we suppose that the product mass balance holds:

1 δt (ρ n y n+1 P -ρ n-1 y n P ) + div(ρ n y k P u n ) = 1 ε ν P W P ω(y n+1 F , y n+1 O ). ( 15
)
Since the sum of the chemical reaction terms vanishes, we have for Σ = y F + y O + y P + y N , summing all the chemical species mass balances,

1 δt (ρ n Σ n+1 -ρ n-1 Σ n ) + div(ρ n Σ k u n ) = 0, (16)
and this equation may equivalently replace the product mass balance equation [START_REF] Piar | A formally second order cell centered scheme for convection-diffusion equations on general grids[END_REF]. Thanks to the mixture balance, we see that, provided that Σ n satisfies Σ n = 1 everywhere in Ω, the solution to Equation ( 16) is Σ n+1 = 1 everywhere in Ω. Since the initialization yields Σ 0 = 1, this last equality is indeed true, and ( 15) is equivalent to (12e). Finally, note that, when the chemical step is performed, the mass balance at step n + 1 is not yet solved; hence the (unusual) backward time shift for the densities and for the mass fluxes in the equations of this step.

Equations (12g)-(12h) implement a pressure correction technique, where the correction step couples the velocity correction equation, the mass balance and the sensible enthalpy balance. This coupling ensures that the pressure and velocity are kept constant through the contact discontinuity associated to compositional non-reactive Euler equations (precisely speaking, the usual contact discontinuity, already present in 1D equations, but not slip lines); for this property to hold, it is necessary that all chemical species share the same heat capacity ratio γ. The term S n+1 K in the sensible enthalpy balance equation is a corrective term which is necessary for consistency; schematically speaking, it compensates the numerical dissipation which appears in a discrete kinetic energy balance that is obtained from the discrete momentum balance. Its expression is given in Section 5, and its derivation is explained in Section 6, where the conservativity of the scheme is discussed.

Space discretization

The space dicretization is performed by a finite volume technique, using a staggered arrangement of the unknowns (the scalar variables are approximated at the cell centers and the velocity components at the face centers), using either a MAC scheme (for structured discretizations) or the degrees of freedom of low-order non-conforming finite elements: Crouzeix-Raviart [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations[END_REF] for simplicial cells and Rannacher-Turek [START_REF] Rannacher | Simple nonconforming quadrilateral Stokes element[END_REF] for quadrangles (d = 2) or hexahedra (d = 3). For the Euler equations (i.e. Steps (12g)-(12h)), upwinding is performed by building positivity-preserving convection operators, in the spirit of the so-called Flux-Splitting methods, and only first-order upwinding is implemented. The pressure gradient is built as the transpose (with respect to the L 2 inner product) of the natural velocity divergence operator. For the balance equations for the other scalar unknowns, the time discretization is implicit when first-order upwinding is used in the convection operator (in other words, k = n+ 1 in (12a)-(12d)) or explicit (k = n in (12a)-(12d)) when a higher order (of MUSCL type, cf. Appendix 8) flux or an anti-diffusive flux (cf. Appendix 9) is used.

Properties of the scheme

First, the positivity of the density is ensured by construction of the discrete mass balance equation, i.e. by the use of a first order upwind scheme. In addition, the physical bounds of the mass fractions are preserved thanks to the following (rather standard) arguments: first, building a discrete convection operator which vanishes when the convected unknown is constant thanks to the discrete mass balance equation ensures a positivity-preservation property [START_REF] Larrouturou | How to preserve the mass fractions positivity when computing compressible multi-component flows[END_REF], under a CFL condition if an explicit time approximation is used; second, the discretization of the chemical reaction rate ensures either that it vanishes when the unknown of the equation vanishes (for y F and y O ), or that it is non-negative (for y P ). Consequently, mass fractions are non-negative and, since their sum is equal to 1 (see above), they are also bounded by 1.

The positivity of the sensible energy stems from two essential arguments: first, a discrete analog of the internal energy equation ( 8) may be obtained from the discrete sensible enthalpy balance, by mimicking the continuous computation; second, this discrete relation may be shown to have only positive solutions, once again thanks to the consistency of the discrete convection operator and the mass balance. This holds provided that the equation is exothermic ( ωθ ≥ 0) and thanks to the non-negativity of S n+1 (see below).

In order to calculate correct shocks, it is crucial for the scheme to be consistent with the following weak formulation of the problem:

∀φ ∈ C ∞ c (Ω × [0, T ) , T 0 Ω ρ∂ t φ + ρu • ∇φ dx dt + Ω ρ 0 (x)φ(x, 0)dx = 0, T 0 Ω ρu i ∂ t φ + (ρuu i ) • ∇φ + p∂ i φ dx dt + Ω ρ 0 (x)(u i ) 0 (x)φ(x, 0)dx = 0, 1 ≤ i ≤ d, T 0 Ω ρE∂ t φ + (ρE + p)u • ∇φ dx dt + Ω ρ 0 (x)E 0 (x)φ(x, 0)dx = 0, T 0 Ω ρy i ∂ t φ + ρy i u • ∇φ dx dt + T 0 Ω ρ 0 (x)y i,0 (x)φ(x, 0)dx = - T 0 Ω ωi φ dx dt, 1 ≤ i ≤ d, p = (γ -1)ρe s . (17)
Remark that this system features the total energy balance equation and not the sensible enthalpy balance equation, which is actually solved here. However, we show in Section 6 that the solutions of the scheme satisfy a discrete total energy balance, with a time and space dicretization which is unusual but allows however to prove the consistency in the Lax-Wendroff sense. Finally, the integral of the total energy over the domain is conserved, which yields a stability result for the scheme (irrespectively of the time and space step, for this relation; recall however that the overall stability of the scheme needs a CFL condition if an explicit version of the convection operator for chemical species is used).

Meshes and unknowns

Let the computational domain Ω be an open polygonal subset of R d , 1 ≤ d ≤ 3, with boundary ∂Ω and let M be a decomposition of Ω, supposed to be regular in the usual sense of the finite element literature (e.g. [START_REF] Ciarlet | Basic error estimates for elliptic problems[END_REF]). The cells may be:

-for a general domain Ω, either convex quadrilaterals (d = 2) or hexahedra (d = 3) or simplices, both type of cells being possibly combined in a same mesh for two-dimensional problems, -for a domain the boundaries of which are hyperplanes normal to a coordinate axis, rectangles (d = 2) or rectangular parallelepipeds (d = 3) (the faces of which, of course, are then also necessarily normal to a coordinate axis). By E and E(K) we denote the set of all (d -1)-faces σ of the mesh and of the element K ∈ M respectively. The set of faces included in the boundary of Ω is denoted by E ext and the set of internal edges (i.e. E \ E ext ) is denoted by E int ; a face σ ∈ E int separating the cells K and L is denoted by σ = K|L. The outward normal vector to a face σ of K is denoted by n K,σ . For K ∈ M and σ ∈ E, we denote by |K| the measure of K and by |σ| the (d -1)-measure of the face σ. The size of the mesh is denoted by h:

h = max diam(K), K ∈ M . For 1 ≤ i ≤ d, we denote by E (i) ⊂ E and E (i)
ext ⊂ E ext the subset of the faces of E and E ext respectively which are perpendicular to the i th unit vector of the canonical basis of R d .

The space discretization is staggered, using either the Marker-And Cell (MAC) scheme [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF][START_REF] Harlow | A numerical fluid dynamics calculation method for all flow speeds[END_REF], or nonconforming low-order finite element approximations, namely the Rannacher and Turek (RT) element [START_REF] Rannacher | Simple nonconforming quadrilateral Stokes element[END_REF] for quadrilateral or hexahedric meshes, or the lowest degree Crouzeix-Raviart (CR) element [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations[END_REF] for simplicial meshes.

For all these space discretizations, the degrees of freedom for the pressure, the density, the enthalpy, the mixture, fuel and neutral gas mass fractions and the flame indicator are associated to the cells of the mesh M and are denoted by:

p K , ρ K , h K , y F,K , y N,K , z K , G K , K ∈ M .
Let us then turn to the degrees of freedom for the velocity (i.e. the discrete velocity unknowns).

-Rannacher-Turek or Crouzeix-Raviart discretizations -The degrees of freedom for the velocity components are located at the center of the faces of the mesh, and we choose the version of the element where they represent the average of the velocity through a face. The set of degrees of freedom reads:

{u σ , σ ∈ E}, of components {u σ,i , σ ∈ E, 1 ≤ i ≤ d}.
-MAC discretization -The degrees of freedom for the i th component of the velocity are defined at the centre of the faces of E (i) , so the whole set of discrete velocity unknowns reads:

u σ,i , σ ∈ E (i) , 1 ≤ i ≤ d .
For the definition of the schemes, we need a dual mesh which is defined as follows.

-Rannacher-Turek or Crouzeix-Raviart discretizations -For the RT or CR discretizations, the dual mesh is the same for all the velocity components. When K ∈ M is a simplex, a rectangle or a rectangular cuboid, for σ ∈ E(K), we define D K,σ as the cone with basis σ and with vertex the mass center of K (see Figure 1). We thus obtain a partition of K in m sub-volumes, where m is the number of faces of the mesh, each sub-volume having the same measure

|D K,σ | = |K|/m.
We extend this definition to general quadrangles and hexahedra, by supposing that we have built a partition still of equal-volume subcells, and with the same connectivities; note that this is of course always possible, but that such a volume D K,σ may be no longer a cone; indeed, if K is far from a parallelogram, it may not be possible to build a cone having σ as basis, the opposite vertex lying in K and a volume equal to |K|/m (note that these dual cells do not need to be constructed in the implementation of the scheme, only their volume is needed). The volume D K,σ is referred to as the half-diamond cell associated to K and σ. In order to be able to write a unique expression of the discrete equations for both MAC and CR/RT schemes, we introduce the set of faces E (i) S associated with the degrees of freedom of each component of the velocity (S stands for "scheme"):

D σ D σ ′ σ ′ = K |M K L M |σ | σ = K |L ε = D σ |D σ ′
E (i) S = E (i) \ E (i) ext for the MAC scheme, E \ E (i)
ext for the CR or RT schemes. Similarly, we unify the notation for the set of dual faces for both schemes by defining:

Ẽ(i) S = Ẽ(i) \ Ẽ(i)
ext for the MAC scheme, Ẽ \ Ẽ(i) ext for the CR or RT schemes, where the symbol ˜refers to the dual mesh; for instance, Ẽ(i) is thus the set of faces of the dual mesh associated with the i th component of the velocity, and Ẽ(i) ext stands for the subset of these dual faces included in the boundary. Note that, for the MAC scheme, the faces of Ẽ(i) are perpendicular to a unit vector of the canonical basis of R d , but not necessarily to the i th one.

The scheme

In this section, we give the fully discrete form of the scheme. Even if it corresponds to the reverse order with respect to the semi-discrete scheme given in [START_REF] Herbin | Kinetic energy control in the MAC discretisation of the compressible Navier-Stokes equations[END_REF], we begin with the hydrodynamics (Section 5.1) and then turn to the mass balance step for chemical species and the transport of the characteristic function for the burnt zone (Section 5.2). This choice is due to the fact that the definition of the convection operators for scalar variables necessitates to introduce the discretization of the mixture mass balance equation first.

Euler step

For 0 ≤ n < N , the step n + 1 of the algorithm for the resolution of the Euler equations reads:

Pressure gradient scaling step -Solve for ( ∇p) n+1 :

∀σ ∈ E, ( ∇p) n+1 σ = ρ n Dσ ρ n-1 Dσ 1/2 (∇p) n σ . ( 18a
)
Prediction step -Solve for ũn+1 :

For 1 ≤ i ≤ d, ∀σ ∈ E (i) S , 1 δt (ρ n Dσ ũn+1 σ,i -ρ n-1 Dσ u n σ,i ) + div σ (ρ n ũn+1 i u n ) + ( ∇p) n+1 σ,i = 0. ( 18b
)
Correction step -Solve for ρ n+1 , p n+1 and u n+1 :

For 1 ≤ i ≤ d, ∀σ ∈ E (i) S , 1 δt ρ n Dσ (u n+1 σ,i -ũn+1 σ,i ) + (∇p) n+1 σ,i -( ∇p) n+1 σ,i = 0, (18c)
∀K ∈ M, 1 δt (ρ n+1 K -ρ n K ) + div K (ρu) n+1 = 0, (18d)
∀K ∈ M, 1 δt ρ n+1 K (h s ) n+1 K -ρ n K (h s ) n K + div K (ρh s u) n+1 - 1 δt (p n+1 K -p n K ) -u • ∇p n+1 K = ( ωθ ) n+1 K + S n+1 K , (18e)
∀K ∈ M, p n+1 K = γ -1 γ (h s ) n+1 K ρ n+1 K . (18f)
The initial approximations for ρ -1 , h 0 s and u 0 are given by the mean values of the initial conditions over the primal and dual cells:

∀K ∈ M, ρ -1 K = 1 |K| K ρ 0 (x) dx and (h s ) 0 K = 1 |K| K (h s ) 0 (x) dx, ∀σ ∈ E (i) S , 1 ≤ i ≤ d, u 0 σ,i = 1 |D σ | Dσ (u 0 (x)) i dx.
Then, ρ 0 is computed by the mass balance equation (18d) and p 0 is computed by the equation of state (18f).

We now define each of the discrete operators featured in System [START_REF] Zimont | Gas premixed combustion at high turbulence. Turbulent flame closure combustion model[END_REF].

Mass balance equation

Equation ( 18d) is a finite volume discretisation of the mass balance (4a) over the primal mesh. For a discrete density field ρ and a discrete velocity field u, the discrete divergence is defined by:

div K (ρu) = 1 |K| σ∈E(K) F K,σ , F K,σ = |σ| ρ σ u K,σ ,
where u K,σ is an approximation of the normal velocity to the face σ outward K. The definition of this latter quantity depends on the discretization: in the MAC case, u K,σ = u σ,i e (i) • n K,σ for a face σ of K perpendicular to e (i) , with e (i) the i-th vector of the orthonormal basis of R d , and, in the CR and RT cases, u K,σ = u σ • n K,σ for any face σ of K. The density at the face σ = K|L is approximated by the upwind technique, so ρ σ = ρ K if u K,σ ≥ 0 and ρ σ = ρ L otherwise. Since we assume that the normal velocity vanishes on the boundary faces, the definition is complete.

Convection operators associated to the primal mesh

We may now give the general form of the discrete convection operator of any discrete field z defined on the primal cell:

div K (ρzu) = 1 |K| σ∈E(K) F K,σ z σ , (19)
where z σ is an approximation of the unknown z at the face σ.

Momentum balance equation and pressure gradient scaling

We now turn to the discrete momentum balance (18b). For the MAC discretization, but also for the RT and CR discretizations, the time derivative and convection terms are approximated in (18b) by a finite volume technique over the dual cells, so the convection term reads:

div σ (ρũ i u) = div σ ũi (ρu) = 1 |D σ | ε∈ Ẽ(Dσ) F σ,ε ũε,i ,
where F σ,ε stands for a mass flux through the dual face ε, and ũε,i is a centered approximation of the i th component of the velocity ũ on ε. The density in the dual cell ρ Dσ is obtained by a weighted average of the density in the neighbour cells:

|D σ | ρ Dσ = |D K,σ | ρ K + |D L,σ | ρ L for σ = K|L ∈ E int ,
and ρ Dσ = ρ K for an external face of a cell K. The mass fluxes (F σ,ε ) ε∈E(Dσ ) are evaluated as linear combinations, with constant coefficients, of the primal mass fluxes at the neighbouring faces, in such a way that the following discrete mass balance over the dual cells is implied by the discrete mass balance (18d):

∀σ ∈ E and n ∈ N, |D σ | δt (ρ n+1 Dσ -ρ n Dσ ) + ε∈E(Dσ) F n+1 σ,ε = 0. ( 20
)
This relation is critical to derive a discrete kinetic energy balance (see Section 6 below). The computation of the dual mass fluxes is such that the flux through a dual face lying on the boundary, which is then also a primal face, is the same as the primal flux, that is zero. For the expression of these fluxes, we refer to [START_REF] Gastaldo | Staggered discretizations, pressure correction schemes and all speed barotropic flows[END_REF][START_REF] Herbin | On some implicit and semi-implicit staggered schemes for the shallow water and Euler equations[END_REF][START_REF] Herbin | Kinetic energy control in the MAC discretisation of the compressible Navier-Stokes equations[END_REF]. Since the mass balance is not yet solved at the velocity prediction stage, they have to be built from the mass balance at the previous time step: hence the backward time shift for the densities in the time-derivative term.

The term (∇p) σ,i stands for the i-th component of the discrete pressure gradient at the face σ. This gradient operator is built as the transpose of the discrete operator for the divergence of the velocity, i.e. in such a way that the following duality relation with respect to the L 2 inner product holds:

K∈M |K|p K div K (u) + d i=1 σ∈E (i) S |D σ |u σ,i (∇p) σ,i = 0.
This leads to the following expression:

∀σ = K|L ∈ E int , (∇p) σ,i = |σ| |D σ | (p L -p K ) n K,σ • e (i) .
The scaling of the pressure gradient (18a) is necessary for the solution to the scheme to satisfy a local discrete (finite volume) kinetic energy balance [8, Lemma 4.1].

Sensible enthalpy equation

The convection term for the sensible enthalpy takes the form (19), with an implicit and upwind (with respect to the mass flux F K,σ ) approximation of the unknown at the face. In addition, this equation is discretized in such a way that the present enthalpy formulation is strictly equivalent to the internal energy formulation of the energy balance equation used in [START_REF] Grapsas | An unconditionally stable finite element-finite volume pressure correction scheme for the compressible Navier-Stokes equations[END_REF]. Consequently, the term -u • ∇p K reads:

-u • ∇p K = 1 |K| σ∈E(K) |σ| u K,σ (p K -p σ ),
where p σ is the upwind approximation of p at the face σ with respect to u K,σ . The reaction heat, ( ωθ ) K , is written in the following way:

( ωθ ) K = - Ns i=1 ∆h 0 f,i ( ωi ) K = ν F W F ∆h 0 f,F +ν O W O ∆h 0 f,O -ν P W P ∆h 0 f,P ωK .
The definition of ωK is given in Section 5.2, and the definition of the corrective term S n+1 K is given in Section 6 (see Equation (31) and Remark 3 below).

Chemistry step

For 0 ≤ n < N , the step n + 1 for the solution of the transport of the characteristic function of the burnt zone and the chemical species mass balance equations reads:

Computation of the burnt zone characteristic function -Solve for G n+1 :

∀K ∈ M, 1 δt (ρ n K G n+1 K -ρ n-1 K G n K ) + div K (ρ n G k u n ) + (ρ u u f |∇G k |) K = 0. (21a)
Computation of the variable z -Solve for z n+1 :

∀K ∈ M, 1 δt (ρ n K z n+1 K -ρ n-1 K z n K ) + div K (ρ n z k u n ) = 0. ( 21b
)
Neutral gas mass fraction computation -Solve for y n+1 N :

∀K ∈ M, 1 δt ρ n K (y N ) n+1 K -ρ n-1 K (y N ) n K + div K (ρ n y k N u n ) = 0. ( 21c
)
Fuel mass fraction computation -Solve for y n+1 F :

∀K ∈ M, 1 δt ρ n K (y F ) n+1 K -ρ n-1 K (y F ) n K + div K (ρ n y k F u n ) = - 1 ε ν F W F ωn+1 K . ( 21d
)
Product mass fraction computation -Compute y n+1 P given by:

∀K ∈ M, (y P ) n+1 K = 1 -(y F ) n+1 K -(y O ) n+1 K -(y N ) n+1 K . (21e)
The initial value of the chemical variables is the mean value of the initial condition over the primal cells:

∀K ∈ M, G 0 K = 1 |K| K G 0 (x) dx, z 0 K = 1 |K| K z 0 (x) dx, (y i ) 0 K = 1 |K| K (y i ) 0 (x) dx, i = N, F,
where the reduced variable z is the linear combination of y F and y O given by Equation [START_REF] Larrouturou | How to preserve the mass fractions positivity when computing compressible multi-component flows[END_REF].

In Equations (21a)-(21d), the discretization of the convection terms is performed by a discrete operator of the form (19). Several choices are possible (and compared in numerical tests) for the evaluation of the value at the face: either an implicit scheme (i.e. k = n + 1) with a first-order upwind space discretization, either an explicit scheme (i.e. k = n) with a MUSCL or an anti-diffusive space approximation. These latter discretizations are described in Section 8 and Section 9 of the appendix, respectively.

According to the developments of Section 3, the chemical reaction term reads ωn+1

K = η((y F ) n+1 K , z n+1 K ) (G n+1 K -0.5) -with η((y F ) n+1 K , z n+1 K ) = 1 ν F W F (y F ) n+1 K if z n+1 ≤ 0, 1 ν F W F (y F ) n+1 K -z n+1 K otherwise,
and the chemical species mass fractions satisfy the following system, which is equivalent to (21b)-(21e):

1 δt ρ n K (y i ) n+1 K -ρ n-1 K (y i ) n K + div K (ρ n y k i u n ) = 1 ε ζ i ν i W i ωn+1 , for i ∈ I and K ∈ M. ( 22
)
At the continuous level, the last term of equation (21a) may be written

ρ u u f |∇G| = a • ∇G = div(G a) -G div(a), with a = ρ u u f ∇G |∇G| ,
and we use the same decomposition at the discrete level:

|K| (ρ u u f |∇G|) K = σ∈E(K) |σ| (G k σ -G k K ) a n σ • n K,σ ,
where G k σ may be given by one of the three above-mentioned schemes, namely an implicit upwind (with respect to a n • n K,σ ) scheme, an explicit MUSCL or an explicit anti-diffusive scheme. The flame velocity on σ, a n σ , is evaluated as

a n σ = ρ u u f (∇G) n σ |(∇G) n σ |
, where the gradient of G on σ = K|L is computed as:

(∇G) σ = 1 |K ∪ L| τ ∈E(K) |τ | Ĝτ n K,τ + τ ∈E(L) |τ | Ĝτ n L,τ ,
where Ĝτ is a second order approximation of G at the center of the face τ .

Scheme conservativity

Let the discrete sensible internal energy be defined by

p n K = (γ -1) ρ n K (e s ) n K for K ∈ M and 0 ≤ n ≤ N . In view of the equation of state (18f), this definition implies ρ n K (h s ) n K = ρ n K (e s ) n K + p n K , for K ∈ M and 0 ≤ n ≤ N .
The following lemma states that the discrete solutions satisfy a local internal energy balance.

Lemma 1 (Discrete internal energy balance).

A solution to (18)-( 21) satisfies the following equality, for any K ∈ M and 0 ≤ n < N :

1 δt (ρe) n+1 K -(ρe) n K + div K (ρeu) n+1 + p n+1 K div K (u) n+1 = S n+1 K , ( 23
)
where

(ρe) n+1 K = ρ n+1 K (e s ) n+1 K + ρ n K i∈I ∆h 0 f,i (y i ) n+1 K , div K (ρeu) n+1 = div K (ρe s ) n+1 u n+1 + ρ n i∈I ∆h 0 f,i y k i u n .
Proof. We begin by deriving a local sensible internal energy balance, starting from the sensible enthalpy balance (18e) and mimicking the previously given formal passage between these two equations at the continuous level (i.e. the passage from Equation [START_REF] Herbin | On some implicit and semi-implicit staggered schemes for the shallow water and Euler equations[END_REF] to Equation ( 10)). To this purpose, let us write (18e) as T 1 + T 2 = T 3 with

T 1 = 1 δt ρ n+1 K (h s ) n+1 K -ρ n K (h s ) n K - 1 δt (p n+1 K -p n K ), T 2 = div K (ρh s u) n+1 -u • ∇p n+1 K , T 3 = ( ωθ ) n+1 K + S n+1 K . Using ρ ℓ K (h s ) ℓ K = ρ ℓ K (e s ) ℓ K + p ℓ K for ℓ = n and ℓ = n + 1, we easily get T 1 = 1 δt ρ n+1 K (e s ) n+1 K -ρ n K (e s ) n K .
The term T 2 reads:

|K| T 2 = σ∈E(K) |σ| ρ n+1 σ (h s ) n+1 σ -p n+1 σ + p n+1 K u n+1 K,σ . If u n+1 K,σ > 0, by definition, ρ n+1 σ (h s ) n+1 σ = ρ n+1 K (h s ) n+1
ρ n+1 σ (h s ) n+1 σ -p n+1 σ = ρ n+1 σ (e s ) n+1 σ , so, finally |K| T 2 = σ∈E(K) F n+1 K,σ (e s ) n+1 σ + p n+1 K σ∈E(K) |σ| u n+1 K,σ .
We thus get the following sensible internal energy balance:

|K| δt ρ n+1 K (e s ) n+1 K -ρ n K (e s ) n K + σ∈E(K) F n+1 K,σ (e s ) n+1 σ + p n+1 K σ∈E(K) |σ| u n+1 K,σ = |K| ( ωθ ) n+1 K + S n+1 K , (24)
or, using the discrete differential operator formalism,

1 δt ρ n+1 K (e s ) n+1 K -ρ n K (e s ) n K + div K (ρe s u) n+1 + p n+1 K div K u n+1 = ( ωθ ) n+1 K + S n+1 K . ( 25
)
We now derive from this relation a discrete (sensible and chemical) internal energy balance. Multiplying the mass fraction balance equations by the corresponding formation enthalpy (∆h 0 f,i ) i∈I and summing over i ∈ I yields:

1 δ t i∈I ∆h 0 f,i ρ n K (y i ) n+1 K -ρ n+1 K (y i ) n K + σ∈E(K) F n K,σ i∈I ∆h 0 f,i (y i ) k σ = i∈I ∆h 0 f,i ( ωi ) n+1 K = ( ωθ ) n+1 K .
Adding this relation to (24) yields the balance equation ( 23).

Remark 1 (Positivity of the sensible internal energy). Equation ( 25) implies that the sensible internal energy remains positive, provided that the right-hand side is non-negative, which is true if ωθ ≥ 0, i.e. if the chemical reaction is exothermic. The proof of this property may be found in [START_REF] Grapsas | An unconditionally stable finite element-finite volume pressure correction scheme for the compressible Navier-Stokes equations[END_REF]Lemma 4.3], and relies on two arguments: first, the convection operator may be recast as a discrete positivity-preserving transport operator thanks to the mass balance, and, second, the pressure p n+1 K vanishes when (e s ) n+1 K vanishes, by the equation of state.

The following local discrete kinetic energy balance holds on the dual mesh (see [START_REF] Grapsas | An unconditionally stable finite element-finite volume pressure correction scheme for the compressible Navier-Stokes equations[END_REF]Lemma 4.1] for a proof).

Lemma 2 (Discrete kinetic energy balance on the dual mesh).

A solution to (18)-( 21) satisfies the following equality, for

1 ≤ i ≤ d, σ ∈ E (i) S and 0 ≤ n < N : |D σ | δt (e k ) n+1 σ,i -(e k ) n σ,i + ε∈ Ẽ(Dσ) F n σ,ε (e k ) n+1 ε,i + |D σ |(∇p) n+1 σ,i u n+1 σ,i = -R n+1 σ,i , (26)
where

(e k ) n+1 σ,i = 1 2 ρ n Dσ (u n+1 σ,i ) 2 + δt 2 |D σ | 2ρ n Dσ (∇p) n+1 σ,i 2
,

(e k ) n+1 ε,i = 1 2 ũn+1 σ,i ũn+1 σ ′ ,i , R n+1 σ,i = |D σ | ρ n-1 Dσ 2δt (ũ n+1 σ,i -u n σ,i ) 2 .
We now derive a kinetic energy balance equation on the primal cells from Relation (26). For the sake of clarity, we make a separate exposition for the Rannacher-Turek case and the MAC case. The case of simplicial discretizations, with the degrees of freedom of the Crouzeix-Raviart element, is an easy extension of the Rannacher-Turek case.

The Rannacher-Turek case

Since the dual meshes are the same for all the velocity components in this case, we may sum up Equation (26) over i = 1, . . . d to obtain, for σ ∈ E and 0 ≤ n < N :

|D σ | δt (e k ) n+1 σ -(e k ) n σ + ε∈ Ẽ(Dσ) F n σ,ε (e k ) n+1 ε + |D σ |(∇p) n+1 σ • u n+1 σ = -R n+1 σ , (27)
K ε1 -Fσ ,ε 1 ε2 -F σ , ε2 ε3 F σ , ε3 ε 4 F σ , ε 4 σ GK,σ = 1 2 (-Fσ,ε 1 (e k )ε 1 -Fσ,ε 2 (e k )ε 2 +Fσ,ε 3 (e k )ε 3 + Fσ,ε 4 (e k )ε 4 )
Fig. 2: From fluxes at dual faces to fluxes at primal faces, for the Rannacher-Turek discretization.

with

(e k ) ℓ σ = d i=1 (e k ) ℓ σ,i , for ℓ = n or ℓ = n + 1, (e k ) n+1 ε = d i=1 (e k ) n+1 ε,i and R n+1 σ = d i=1 R n+1 σ,i .
For K ∈ M, let us define a kinetic energy associated to K and the flux G n+1 K,σ as follows (see Figure 2):

(e k ) ℓ K = 1 2 |K| σ∈E(K) |D σ | (e k ) ℓ σ , ℓ = n or ℓ = n + 1, G n+1 K,σ = - 1 2 ε∈E(Dσ),ε⊂K F n σ,ε (e k ) n+1 ε + 1 2 ε∈E(Dσ ),ε ⊂K F n σ,ε (e k ) n+1 ε .
We easily check that the fluxes G n+1 K,σ are conservative, in the sense that, for

σ = K|L, G n+1 K,σ = -G n+1 L,σ .
Let us now divide Equation ( 27) by 2 and sum over the faces of K. A reordering of the summations, using the conservativity of the mass fluxes through the dual edges and the expression of the discrete pressure gradient, yields:

|K| δt (e k ) n+1 K -(e k ) n K + σ∈E(K) G n+1 K,σ + 1 2 σ=K|L |σ| (p n+1 L -p n+1 K ) u n+1 K,σ = -R n+1 K , with R n+1 K = 1 2 σ∈E(K) R n+1 σ . (28)
The MAC case

Let 1 ≤ i ≤ d, let K ∈ M,
let us denote by σ and σ ′ the two faces of E (i) (K), and let us define:

(e k ) ℓ K,i = 1 2 |K| |D σ | (e k ) ℓ σ,i + |D σ | (e k ) ℓ σ,i , for ℓ = n or ℓ = n + 1.
Case of primal faces parallel to the dual faces. Let τ = σ or τ = σ ′ , let ε 1 and ε 2 be the two faces of D τ perpendicular to e (i) , and let ε 2 be the one included in K (see Figure 3). Then we define

G n+1 K,τ,i = 1 2 F τ,ε1 (e k ) n+1 ε1,i -F τ,ε2 (e k ) n+1 ε2,i . K σ σ ′ ε1 Fσ,ε 1 ε2 -Fσ,ε 2 GK,σ,1 = 1 2 Fσ,ε 1 (e k )ε 1 ,1 -Fσ,ε 2 (e k )ε 2 ,1
Fig. 3: From fluxes at dual faces to fluxes at primal faces, for the MAC discretization, primal faces parallel to the dual edges, first component of the velocity.

Case of primal faces orthogonal to the dual faces. For τ ∈ E(K) \ {σ, σ ′ }, let ε and ε ′ be such that τ ⊂ (ε ∪ ε′ ) with ε a face of D σ and ε ′ a face of D σ ′ (see Figure 4).

K τ σ σ ′ ε Fσ,ε ε ′ F σ ′ ,ε ′ GK,τ,1 = 1 2 Fσ,ε(e k )ε,1 + F σ ′ ,ε ′ (e k ) ε ′ ,1
Fig. 4: From fluxes at dual faces to fluxes at primal faces, for the MAC discretization, primal faces orthogonal to the dual edges, first component of the velocity.

Then we define

G n+1 K,τ,i = 1 2 F σ,ε (e k ) n+1 ε,i + F σ ′ ,ε ′ (e k ) n+1 ε ′ ,i .
Summing Equation (26) written for σ and for σ ′ and dividing the result by 2 yields:

|K| δt (e k ) n+1 K,i -(e k ) n K,i + σ∈E(K) G n+1 K,σ,i + 1 2 σ∈E (i) (K) σ=K|L |σ| (p n+1 L -p n+1 K ) u n+1 K,σ = - 1 2 R n+1 σ,i + R n+1 σ,i . (29)
Now let (e k ) ℓ K = d i=1 (e k ) ℓ K,i , for ℓ = n or ℓ = n + 1, and

G n+1 K,σ = d i=1 G n+1 K,σ,i , for σ ∈ E(K).
Since only one equation is written for a given face σ of the mesh (for the velocity component i with i such that the normal vector to σ is parallel to e (i) ), we may define in the MAC case R n+1 σ = R n+1 σ,i . Summing Equation (29) over the space dimension, we finally get

|K| δt (e k ) n+1 K -(e k ) n K + σ∈E(K) G n+1 K,σ + 1 2 σ=K|L |σ| (p n+1 L -p n+1 K ) u n+1 K,σ = -R n+1 K , with R n+1 K = 1 2 σ∈E(K) R n+1 σ , (30)
which is formally the same equation as Relation (28) (although with a different definition of all the terms in the equation except the pressure gradient).

Remark 2 (On the definition of the cell kinetic energy). Note that, both in the Rannacher-Turek and the MAC case, the cell kinetic energy is not a convex combination of the face kinetic energies, since, on a non-uniform mesh, the equalities

|K| = 1 2 σ∈E(K) |D σ | (Rannacher Turek case) and |K| = 1 2 σ∈E (i) (K) |D σ | (MAC case) do not hold in general.
Consequently, the cell kinetic energy may oscillate from cell to cell while the face kinetic energy does not. Nevertheless, the discrete time derivative of the cell kinetic energy is consistent in the Lax-Wendroff sense, because, despite of these oscillations, the cell kinetic energy still converges weakly if the velocity converges. Equations ( 28) and (30) suggest a choice for the term S n+1 K , the purpose of which is to compensate the numerical dissipation terms appearing in the kinetic energy balance:

S n+1 K = R n+1 K , for K ∈ M and 0 ≤ n < N. (31)
This expression yields a conservative scheme, in the sense that the discrete solutions satisfy a discrete total energy balance without any remainder term (see Equation (4c) below); as a consequence, the scheme can be proven to be consistent in the Lax-Wendroff sense. However, different definitions are possible (and this latitude may be useful in explicit variants of the scheme, to ensure the positivity of S n+1 K see Remark 3 below). We are now in position to state a total energy balance for the scheme.

Theorem 1 (Discrete total energy and stability of the scheme). A solution to (18)-(21) satisfies the following discrete total energy balance, for any K ∈ M and 0 ≤ n < N :

1 δt (ρE) n+1 K -(ρE) n K + div K ((ρE + p)u) n+1 = 0, (32)
where

(ρE) ℓ K = (e k ) ℓ K + ρ ℓ K (e s ) ℓ K + ρ l-1 K i∈I ∆h 0 f,i (y i ) ℓ K , for ℓ = n and ℓ = n + 1, div K ((ρE + p) u) n+1 = div K (ρe s ) n+1 u n+1 + ρ n i∈I ∆h 0 f,i y k i u n + 1 |K| σ∈E(K) G n+1 K,σ + 1 |K| σ=K|L |σ| p n+1 K + p n+1 L 2 u n+1 K,σ .
Let us suppose that e 0 s , ρ 0 and ρ -1 are positive. Then, a solution to (18)-( 21) satisfies ρ n > 0, e n s > 0 and the following stability result:

E n = E 0 ,
where, for 0 ≤ n ≤ N ,

E n = K∈M |K|(ρe) n K + 1 2 d i=1 σ∈E (i) S |D σ | ρ n-1 Dσ (u n σ,i ) 2 + δt 2 2 σ∈Eint |D σ | ρ n-1 Dσ |(∇p) n σ | 2 .
Proof. The discrete total energy balance equation ( 32) is obtained by summing the internal energy balance (23) and the kinetic energy balance, i.e. Equation (28) in the Rannacher-Turek case and Equation (30) for the MAC scheme, and remarking that the numerical dissipation terms in the kinetic energy balance R n+1 K exactly compensate with the corrective terms S n+1 K in the internal energy balance. Then the stability result is obtained by summation over the time steps.

Remark 3 (Consistency of the scheme). The consistency in the Lax-Wendroff sense follows from the conservativity of the scheme (for all balance equations) so, in particular, from the fact that the discrete solutions satisfy the discrete total energy balance (32), thanks to standard (but technical) arguments. Note however that the consistency of the scheme does not require a strict conservativity, and in particular, variants for the choice (31) of the compensation term in the sensible enthalpy balance are possible; indeed, what is really needed is only that the difference between the dissipation in the kinetic energy balance and its compensation tend to zero in a distributional sense. In practice, this allows a different redistribution of the face residuals to the neighbour primal cells, and this can help to preserve the non-negativity of the compensation term for explicit versions of the scheme.

Numerical tests

At the continuous level, the boundedness of the chemical mass fractions formally implies that, when ε → 0, the relaxed model converges to the asymp-totic one. Indeed, integrating any of the reactive species mass balance equations with respect to time and space, we observe that || ω|| L 1 (Ω×(0,T )) tends to zero as ε, and thus two separate zones appear: a zone characterized by G < 0.5 where the reaction is complete, and a zone corresponding to G ≥ 0.5, where no reaction has occured.

A closed form of the solution of the Riemann problem for the asymptotic model is available [START_REF] Beccantini | The reactive Riemann problem for thermally perfect gases at all combustion regimes[END_REF]. In order to perform numerical tests, a Riemann problem with initial conditions such that the analytic solution has the profile presented in Figure 5 is chosen. Moreover, the selected configuration imposes zero amplitude for the contact discontinuity and the left non linear wave, thus the solution consists of three different constant states: W * R , W * * and W R . The right state corresponds to a stoichiometric mixture of hydrogen and air (so the molar fractions of Hydrogen, Oxygen and Nitrogen are 2/7, 1/7 and 4/7 respectively) at rest, at the pressure p = 9.9 10 4 Pa and the temperature T = 283 • K. The velocity is supposed to be zero in the left state, which is sufficient to determine the solution. Physically, speaking, supposing that the initial discontinuity lies at x = 0, this situation corresponds to the left part of a (symmetrical) constant velocity plane deflagration starting at x = 0. The flame velocity is u f = 63 m/s and the formation enthalpies are zero except for the product (i.e. steam), with ∆h 0 f,O = -13.255 10 6 J (Kg K) -1 . The quantity ρ u is the analytical density in the intermediate state (so the total velocity of the flame brush is equal to the sum of u f and the material velocity on the right side of the reactive shock, see [START_REF] Beccantini | The reactive Riemann problem for thermally perfect gases at all combustion regimes[END_REF]). The computation is initialized by the analytical solution at t = 0.002 and the final time is t = 0.005. The computational domain is the interval (0, 4.5).

The numerical tests performed aim at checking the convergence of the scheme to such a solution, which in fact may result from two different properties: the convergence of the relaxed model to the asymptotic model when ε tends to zero, and the convergence of the scheme towards a numerical solution when the time and space steps tend to zero. To this purpose, we choose ε proportional to the space step and make it tend to zero, with a constant CFL number. We test the scheme behaviour with three different discretizations of the convection operator in the chemical mass species balances: the standard upwind scheme, a MUSCL-like discretization which is an extension to variable density flows of the scheme proposed in [START_REF] Piar | A formally second order cell centered scheme for convection-diffusion equations on general grids[END_REF] and is described in Appendix 8, and a first-order anti-diffusive scheme which is an adaptation to our setting of the scheme proposed in [START_REF] Després | Contact discontinuity capturing scheme for linear advection and compressible gas dynamics[END_REF]; we detail it in Appendix 9 for the sake of completeness.

Results obtained at t = 0.005 with the upwind scheme, the MUSCL-like scheme and the anti-diffusive scheme, for increasingly refined meshes, are shown on Figure 6, Figure 7 and Figure 8 respectively, together with the analytical solution. The expected convergence is indeed observed but, with the upwind discretization, the rate of convergence is poor. This seems to be due to the interaction between the numerical diffusion of the upwind scheme, which artificially introduces unburnt reactive chemical species into the burnt zone, and the stiffness of the reaction term. As expected in such a case, the results are significantly improved by the use of a less diffusive scheme for the chemical species balance equations. Indeed, passing from the upwind to the MUSCL-like and to the anti-diffusive discretization improves the accuracy of the scheme, as may be observed in Figure 9, where the results obtained by the three discretizations for a regular mesh composed of 500 cells are plotted together with the continuous solution. This observation is comforted by the measure, in L 1 -norm, of the difference between the discrete and continuous solutions, see Table 1. For every mesh and variable, the anti-diffusive scheme is the most accurate and the upwind one the least. The calculated order of convergence is close to 0.5 for the upwind scheme, and to 1 for the MUSCL-like and anti-diffusive schemes. 8 Appendix A: The MUSCL scheme

The MUSCL discretization of the convection operators of the chemical species balance and G-equation closely follows the technique proposed in [START_REF] Piar | A formally second order cell centered scheme for convection-diffusion equations on general grids[END_REF]. To present this discretization, we consider the following system of equations:

∂ t ρ + div(ρu) = 0, ∂ t (ρu) + div(ρuy) = 0.

We suppose for short that this system is complemented by impermeability boundary conditions, i.e. that the normal velocity, both at the continuous and the discrete level, vanishes on the boundary of the computational domain.

The discretization of the above system reads:

∀K ∈ M, ρ n+1 K -ρ n K δt + 1 |K| σ∈E(K)
F n+1 K,σ = 0,

ρ n+1 K y n+1 K -ρ n K y n K δt + 1 |K| σ∈E(K)
F n+1 K,σ y n σ = 0.

For any σ ∈ E, the procedure consists in three steps:

-calculate a tentative value for y σ as a linear interpolate of nearby values, -calculate an interval for y σ which guarantees the stability of the scheme, -project the tentative value y σ on this stability interval.

For the tentative value of y σ , let us choose some real coefficients (α σ K ) K∈M such that

x σ = K∈M α σ K x K , K∈M α σ K = 1.
The coefficients used in this interpolation are chosen in such a way that as few as possible cells, to be picked up in the closest cells to σ, take part. For example, for σ = K|L and if x K , x σ , x L are aligned, only two non-zero coefficients exist in the family (α σ K ) K∈M , namely α σ K and α σ L . Then, these coefficients are used to calculate the tentative value of y σ by

y σ = K∈M α σ K y K .
The construction of the stability interval must be such that the following property holds:

y n+1 K = y n K 1 - δt ρ n+1 K |K| σ∈E(K) β σ K F n+1 K,σ + δt |K| σ∈E(K) y n M σ K β σ K F n+1 K,σ ,
which concludes the proof under the hypothesis that CFL ≤ 1.

We now need to reformulate (33) in order to construct the stability interval. Let σ ∈ E, let us denote by V -and V + the upstream and downstream cell separated by σ, and by V σ (V -) and V σ (V + ) two sets of neighbouring cells of V -and V + respectively, and let us suppose: Remark 5. For σ ∈ E, if V -∈ V σ (V + ), the upstream choice y n σ = y n V -always satisfies the conditions (H1)-(H2), and is the only one to satisfy them if we choose ζ -= ζ + = 0.

(H1) -∃M ∈ V σ (V + ) s.t. u n σ ∈ |[u n M , u n M + ζ + 2 (u n V + -u n M )]|, (H2) -∃M ∈ V σ (V -) s.t. u n σ ∈ |[u n V -, u n V -+ ζ - 2 (u n V --u n M )]|,
Remark 6 (1D case). Let us take the example of an interface σ separating K i and K i+1 in a 1D case (see Figure 10 for the notations), with a uniform meshing and a positive advection velocity, so that V -= K i and V + = K i+1 . In 1D, a natural choice is V σ (K i ) = {K i-1 } and V σ (K i+1 ) = {K i }. On Figure 10, we sketch: on the left, the admissible interval given by (H1) with ζ + = 1 (green) and ζ + = 2 (orange); on the right, the admissible interval given by (H2) with ζ -= 1 (green) and ζ -= 2 (orange). The parameters ζ - and ζ + may be seen as limiting the admissible slope between (x i , y n i ) and (x σ , y n σ ) (with x i the abscissa of the mass centre of K i and x σ the abscissa of σ), with respect to a left and right slope, respectively. For ζ -= ζ + = 1, one recognizes the usual minmod limiter (e.g. [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]Chapter III]). Note that, since, on the example depicted on Figure 10, the discrete function y n has an extremum in K i , the combination of the conditions (H1) and (H2) imposes that, as usual, the only admissible value for y n σ is the upwind one.

9 Appendix B: an anti-diffusive scheme

The scheme proposed in [START_REF] Després | Contact discontinuity capturing scheme for linear advection and compressible gas dynamics[END_REF] by of B. Després and F. Lagoutière for the constant velocity advection problem presents some interesting properties in one space dimension (and may be extended to structured multi-dimensional meshes using alternate directions techniques); in particular, it notably limits the numerical diffusion. We extend here this scheme to work with unstructured meshes for which the "opposite cell to a face" (in the sense introduced in the previous section) may be defined and with a variable density. With the same notations as in the previous section, for σ ∈ E int , σ = K|L with F n+1 K,σ ≥ 0, -the tentative value for y σ is chosen as the downwind value, i.e. y n σ = y n L , -Then we project y n σ on the interval

I σ = y n K , y n K + 1 -ν ν (y K -y M ) , ν = |F n+1 K,σ | δt ρ n+1 K |K| ,
where M ∈ M is the control volume which stands at the opposite side of K with respect to L. The original scheme presented in [START_REF] Després | Contact discontinuity capturing scheme for linear advection and compressible gas dynamics[END_REF] is recovered by this formulation for the one-dimensional constant velocity convection equation. In addition, by arguments similar to those of the previous section, the discretization proposed here may be shown to satisfy a discrete maximum principle.

  For σ ∈ E int , σ = K|L, we now define the diamond cell D σ associated to σ by D σ = D K,σ ∪ D L,σ ; for an external face σ ∈ E ext ∩ E(K), D σ is just the same volume as D K,σ . -MAC discretization -For the MAC scheme, the dual mesh depends on the component of the velocity. For each component, the MAC dual mesh only differs from the RT or CR dual mesh by the choice of the half-diamond cell, which, for K ∈ M and σ ∈ E(K), is now the rectangle or rectangular parallelepiped of basis σ and of measure |D K,σ | = |K|/2. We denote by |D σ | the measure of the dual cell D σ , and by ε = D σ |D σ ′ the dual face separating two diamond cells D σ and D σ ′ .

Fig. 1 :

 1 Fig. 1: Primal and dual meshes for the Rannacher-Turek and Crouzeix-Raviart elements.

K

  to the assumptions on the boundary conditions, σ is an internal face and, denoting by L the adjacent cell to K such that σ = K|L, both cases, denoting by (e s ) n+1 σ the upwind choice for (e s ) n+1 at the face σ, we get

Fig. 5 :

 5 Fig. 5: The analytic solution of the numerical test configuration.

Fig. 6 :

 6 Fig.6: Upwind scheme -From top left to bottom right, fuel mass fraction, G, velocity, pressure, temperature and density at t = 0.005, as a function of the space variable.

  Fig.7: MUSCL scheme -From top left to bottom right, fuel mass fraction, G, velocity, pressure, temperature and density at t = 0.005, as a function of the space variable.

Fig. 8 :Fig. 9 :

 89 Fig.8: Anti-diffusive scheme -From top left to bottom right, fuel mass fraction, G, velocity, pressure, temperature and density at t = 0.005, as a function of the space variable.

  where for a, b ∈ R, we denote by |[a, b]| the interval {αa+(1-α)b, α ∈ [0, 1]}, and ζ + and ζ -are two numerical parameters lying in the interval [0, 2].Conditions (H1)-(H2) and (33) are linked in the following way: let K ∈ M and σ ∈ E(K). If F n K,σ ≤ 0, i.e. K is the downstream cell for σ, denoted above by V + , since ζ + ∈ [0, 2], condition (H1) yields that there existsM ∈ M such that u n σ ∈ |[u n K , u n M ]|, which is (33). Otherwise, i.e. if F n K,σ ≥ 0 and K is the upstream cell for σ, denoted above by V -, condition (H2) yields that there existsM ∈ M such that y n σ ∈ |[y n K , 2y n K -y n M ]|, so y n σ -y n K ∈ |[0, y n K -y n M ]|, which is once again (33).

Table 1 :

 1 7: MUSCL scheme -From top left to bottom right, fuel mass fraction, G, velocity, pressure, temperature and density at t = 0.005, as a function of the space variable.h ||p -p ex || L 1 × 10 -4 ||uu ex || L 1 × 10 -2 ||ρ -ρ ex || L 1 × 10 L 1 normof the error between the discrete and continuous solutions for the various schemes -Black : upwind scheme, blue: MUSCL scheme, orange: anti-diffusive scheme; h 0 = 4.5/250 is the size of the least refined mesh.

	h 0 16.5 7.26 4.59	2.17 1.56	1.07	7.69 3.71 2.74
	h 0 2 h 0 4 h 0 8 h 0 20 h 0 40	12.5 3.88 2.43 9.66 2.05 1.38 7.58 1.17 0.708 0.958 0.263 0.175 3.63 0.691 0.476 1.64 0.787 0.579 6.16 2.23 1.65 1.23 0.471 0.371 4.73 1.26 0.913 5.78 0.673 0.375 0.728 0.160 0.103 2.77 0.382 0.267 4.31 0.414 0.194 0.543 0.0786 0.0458 2.03 0.201 0.134

∀K ∈ M, ∀σ ∈ E(K) ∩ E int , ∃β σ K ∈ [0, 1] and M σ K ∈ M such that

Indeed, under this latter hypothesis and a CFL condition, the scheme preserves the initial bounds of y.

Remark 4. Note that, in Assumption (33), only internal faces are considered, since the fluxes through external faces are supposed to vanish. However, the present discussion may easily be generalized to cope with convection fluxes entering the domain.

Definition 1. The so-called CFL number reads for any 0 ≤ n ≤ N :

Proof. The discrete mass balance equation yields:

Replacing this expression of ρ n K in the discrete balance equation of y and using the relations provided by (33), we obtain:

This relation yields

Fig. 10: Conditions (H1) and (H2) in 1D.

Fig. 11: Notations for the definition of the limitation process. In orange, control volumes of the set V σ (V -) for σ = V -|V + , with a constant advection field F: upwind cells (a) or opposite cells (b).

Finally, we need to specify the choice of the sets V σ (V -) and V σ (V + ). Here, we just set V σ (V + ) = {V -}; such a choice guarantees that at least the upstream choice is in the intersection of the intervals defined by (H1) and (H2), as explained in Remark 6. The set V σ (V -) may be defined in two different ways (cf. Figure 11): -as the "upstream cells" to V -, i.e.

V σ (V -) = {L ∈ M, L shares a face σ with V -and F V -,σ ≤ 0}, -when this makes sense (i.e. with a mesh obtained by Q 1 mappings from the (0, 1) d reference element), the opposite cells to σ in V -are chosen. Note that for a structured mesh, this choice allows to recover the usual minmod limiter.