Christophe Chalons
email: christophe.chalons@uvsq.fr

Mathieu Girardin
email: mathieu.girardin@cea.fr

Samuel Kokh
email: samuel.kokh@cea.fr

An all-regime Lagrange-Projection like scheme for the gas dynamics equations on unstructured meshes

We propose an all regime Lagrange-Projection like numerical scheme for the gas dynamics equations. By all regime, we mean that the numerical scheme is able to compute accurate approximate solutions with an under-resolved discretization, i.e. a mesh size and time step much bigger than the Mach number M . The key idea is to decouple acoustic and transport phenomenon and then alter the numerical flux in the acoustic approximation to obtain a uniform truncation error in term of M . This modified scheme is conservative and endowed with good stability properties with respect to the positivity of the density and the internal energy. A discrete entropy inequality under a condition on the modification is obtained thanks to a reinterpretation of the modified scheme in the Harten Lax and van Leer formalism. A natural extension to multi-dimensional problems discretized over unstructured mesh is proposed. Then a simple and efficient semi implicit scheme is also proposed. The resulting scheme is stable under a CFL condition driven by the (slow) material waves and not by the (fast) acoustic waves and so verifies the all regime property. Numerical evidences are proposed and show the ability of the scheme to deal with tests where the flow regime may vary from low to high Mach values.

Introduction

In this paper, we consider the system of gas dynamics in two space dimension in situations when the flow regime may vary in terms of Mach number M across the computational domain. We propose a collocated Finite Volume method that addresses two important issues.

The first issue concerns the lack of accuracy in the low Mach regime of Godunov-type schemes. While these methods performs well at capturing shocks, they may generate spurious numerical diffusion when they are used for simulating low Mach flows over relatively coarse mesh, i.e. mesh size much bigger the Mach number. Improvements of Godunov-type schemes more generally of collocated methods have been proposed by many authors like [START_REF] Turkel | Preconditioned methods for solving the incompressible and low speed compressible equations[END_REF][START_REF] Guillard | On the behavior of upwind schemes in the low Mach limit[END_REF][START_REF] Klein | Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: One-dimensional flow[END_REF][START_REF] Ph | A projection method for low speed flows[END_REF][START_REF] Degond | Mach-number uniform asymptotic-preserving gauge schemes for compressible flows[END_REF][START_REF] Cordier | An Asymptotic-Preserving all-speed scheme for the Euler and NavierStokes equations[END_REF][START_REF] Thornber | An improved reconstruction method for compressible flows with low Mach number features[END_REF][START_REF] Toumi | Approximate Riemann solvers and flux vector splitting schemes for two-phase flow[END_REF][START_REF] Paillère | Comparison of low mach number models for natural convection problems[END_REF][START_REF] Liou | A sequel to AUSM, part ii: AUSM+-up for all speeds[END_REF][START_REF] Dauvergne | Renormalization of the numerical diffusion for an upwind finite volume method. application to the simulation of Kelvin-Helmholtz instability. Finite volumes for complex applications[END_REF][START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible euler system at low Mach number[END_REF][START_REF] Degond | All speed method for the Euler equation in the low mach number limit[END_REF][START_REF] Liu | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF]. The analysis of these authors may rely on different arguments like the analysis of the viscosity matrix [START_REF] Turkel | Preconditioned methods for solving the incompressible and low speed compressible equations[END_REF], an asymptotic expansion in terms of Mach number [START_REF] Guillard | On the behavior of upwind schemes in the low Mach limit[END_REF], a detailed study in [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible euler system at low Mach number[END_REF] that seek for invariance properties of the numerical scheme transposing the framework of Schochet [START_REF] Schochet | Fast singular limits of hyperbolic PDEs[END_REF] to the discrete setting, and also an analysis based on the so-called Asymptotic Preserving property [START_REF] Jin | Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms[END_REF] in [START_REF] Liu | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF]. Nevertheless the resulting cure usually boils down to reduce the numerical diffusion in the momentum equation for low Mach number values.

The second problem we address deals with subsonic flow when the fluid velocity is slow and the acoustic waves are not driving phenomenons. In this case, the Courant-Friedrichs-Lewy (CFL) condition on the time step for explicit Godunov-type methods that involves the (fast) acoustic wave velocity may lead to very small time steps choices and thus costly computations. It seems natural to seek for numerical schemes that enable the use of a large time steps that are not constrained by the sound velocity. This question has been examined by several authors like [START_REF] Klein | Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: One-dimensional flow[END_REF][START_REF] Cordier | An Asymptotic-Preserving all-speed scheme for the Euler and NavierStokes equations[END_REF][START_REF] Degond | All speed method for the Euler equation in the low mach number limit[END_REF][START_REF] Degond | Mach-number uniform asymptotic-preserving gauge schemes for compressible flows[END_REF][START_REF] Liu | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF] (see also [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF][START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF]) who derived mixed implicit-explicit strategies that allows to choose the time step independently of the Mach Number.

Numerical schemes that can tackle both issues, namely: accuracy for mesh sizes that do not depend on the Mach number and also stability for time steps that are not constrained by the Mach value are usually referred to as all regime, like the methods proposed by [START_REF] Klein | Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: One-dimensional flow[END_REF][START_REF] Cordier | An Asymptotic-Preserving all-speed scheme for the Euler and NavierStokes equations[END_REF][START_REF] Degond | All speed method for the Euler equation in the low mach number limit[END_REF][START_REF] Degond | Mach-number uniform asymptotic-preserving gauge schemes for compressible flows[END_REF][START_REF] Liu | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF].

In the present work, we first propose an operator splitting strategy that allows to decouple the acoustic and the transport phenomenons. The approximation algorithm is split into two steps: an acoustic step and a transport step. For one-dimensional problems, this strategy is equivalent to an explicit Lagrange-Projection [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Després | Lois de conservations Eulériennes, Lagrangiennes et méthodes numériques[END_REF] method, however the present splitting does not involve any moving Lagrangian mesh and can be naturally expressed for multi-dimensional problems. Following simple lines inspired by [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible euler system at low Mach number[END_REF][START_REF] Dellacherie | Construction of modified Godunov type schemes accurate at any Mach number for the compressible Euler system[END_REF] we investigate the dependence of the truncation error with respect to the Mach number. Let us mention that our study does not involve a Taylor expansion in the vicinity of the zero-Mach limit, nor a near-divergence free condition for the velocity field. Although this analysis is by no mean a thorough explanation of the low Mach regime behavior of our solver, it is enough to suggest simple means to obtain a truncation error with a uniform dependence on the Mach number for M < 1. The cure simply relies on modifying the pressure terms in the flux of the acoustic operator that is coherent with the correction proposed by [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible euler system at low Mach number[END_REF][START_REF] Dellacherie | Construction of modified Godunov type schemes accurate at any Mach number for the compressible Euler system[END_REF][START_REF] Toumi | Approximate Riemann solvers and flux vector splitting schemes for two-phase flow[END_REF][START_REF] Paillère | Comparison of low mach number models for natural convection problems[END_REF][START_REF] Dauvergne | Renormalization of the numerical diffusion for an upwind finite volume method. application to the simulation of Kelvin-Helmholtz instability. Finite volumes for complex applications[END_REF]. Although this modified scheme is based on a modified flux definition, one can shows that it can also be rephrased as a simple approximate Riemann solver in the sense of Harten, Lax and van Leer [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF] that is consistent with the integral form of the gas dynamics equation. This scheme is endowed with good stability properties under a CFL condition that involves the Mach number as the time step is still constrained by the sound velocity.

We propose to circumvent this time-step restriction by implementing a mixed implicit-explicit method following the ideas developped by [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF] for one-dimensional problems using a genuine Lagrange-Projection framework. This idea was also used in [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] and consists in using an implicit update for the acoustic step and an explicit march in time for the transport step. This enables stability under a CFL condition that only involves the (slow) material waves without the (fast) acoustic waves. Finally, let us mention that the overall procedure is a conservative discretization that relies on a Suliciu relaxation approach [START_REF] Suliciu | On the thermodynamics of rate-type fluids and phase transitions. i. rate-type fluids[END_REF] that allows to cope with compressible fluids equipped with very general Equation of State (EOS).

The paper is structured as follows: we first present the operator splitting considering only one-dimensional problems. Then we study the behavior of the scheme in the low Mach regime. This allows to lead to an explicit corrected scheme for the sole acoustic step that preserves the accuracy of the scheme at low Mach. Interestingly, we show that this flux-based corrected method may be expressed thanks to an approximate Riemann solver for the acoustic step. Next and thanks to this property, we investigate the ability of the corrected scheme to satisfy to a discrete entropy inequality. Afterwards, we present the extension of the operator splitting method to unstructured meshes either with a semi-implicit or full-explicit march in time. Finally we present numerical results involving low Mach and multi-regime flows.

Governing equations

We are interested in the two-dimensional gas dynamics equations

     ∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = 0, ∂ t (ρE) + ∇ • [(ρE + p)u] = 0, (1a) (1b) (1c)
where ρ, u = (u 1 , u 2 ) t , E denote respectively the density, the velocity vector and the total energy of the fluid. Let e = E -|u| 2 2 be the specific internal energy of the fluid and s its specific entropy. We note τ = 1/ρ and we suppose given an Equation of State (EOS) through the mapping (τ, s) → e EOS which satisfies to the usual Weyl assumptions [START_REF] Weyl | Shock waves in arbitrary fluids[END_REF] ∂ τ e EOS < 0, ∂ s e EOS > 0, ∂ τ τ e EOS > 0

∂ ss e EOS > 0, ∂ τ τ e EOS ∂ ss e EOS > (∂ τ s e EOS ) 2 , ∂ τ τ τ e EOS < 0. ( 2
)
The entropy s = s EOS (τ, e) verifies e = e EOS (τ, s) thanks to (2) and we can define the pressure p = -∂ τ e EOS and the sound velocity c = τ ∂ τ τ e EOS . The above assumptions imply that (τ, s) → e EOS and (τ, e) → -s EOS are strictly convex. Using a slight abuse of notation, we shall also consider p as a function of (τ, e) and note p = p EOS (τ, e).

3 Acoustic/transport operator splitting strategy for the one-dimensional problem

In this section we will consider for the sake of simplicity one-dimensional problems and propose a two-step approximation strategy based on an operator splitting. The aim of this splitting is to decouple acoustic and transport phenomena. Using this guideline we will propose an explicit numerical solver. We shall propose two simple extensions of this method to two-dimensional problems discretized over unstructured grids using either an explicit or a semi-implicit time update in section 5.4. Before going any further, we introduce classical notations for the one-dimensional setting: let ∆t > 0 and ∆x > 0 be respectively the time and space steps. We define the Eulerian mesh interfaces x j+1/2 = j∆x for j ∈ Z, and the intermediate times t n = n∆t for n ∈ N. If b is a fluid parameter, in the sequel, we will note b n j (resp. b n+1 j

) the approximate value b respectively within the

j th cell [x j-1/2 , x j+1/2 ) at instant t = t n (resp. t = t n+1 ).
For one-dimensional problems, (1) supplemented with a passive scalar variable v (that will account for the transverse velocity in two-dimensional problems) reads

           ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p) = 0, ∂ t (ρv) + ∂ x (ρuv) = 0, ∂ t (ρE) + ∂ x [(ρE + p)u] = 0. (3a) (3b) (3c) (3d)
Our discretization strategy of (3) consists in approximating successively the solutions of the following systems ( 4) and [START_REF] Ph | A projection method for low speed flows[END_REF] where

         ∂ t ρ + ρ∂ x u = 0, ∂ t (ρu) + ρu∂ x u + ∂ x p = 0, ∂ t (ρv) + ρv∂ x u = 0, ∂ t (ρE) + ρE∂ x u + ∂ x (pu) = 0, (4a) (4b) (4c) (4d) and          ∂ t ρ + u∂ x ρ = 0, ∂ t (ρu) + u∂ x (ρu) = 0, ∂ t (ρv) + u∂ x (ρv) = 0, ∂ t (ρE) + u∂ x (ρE) = 0. (5a) (5b) (5c) (5d)
In the sequel, system (4) and ( 5) will be respectively referred to as the acoustic system and the transport system.

Given a fluid state (ρ, ρu, ρv, ρE) n j , j ∈ Z at instant t n , this splitting algorithm can be decomposed as follows.

1. Update the fluid state (ρ, ρu, ρv, ρE) n j to the value (ρ, ρu, ρv, ρE) n+1- j by approximating the solution of (4); 2. Update the fluid state (ρ, ρu, ρv, ρE) n+1- j to the value (ρ, ρu, ρv, ρE) n+1 j by approximating the solution of (5).

Properties and approximation of the one-dimensional acoustic system

First, we notice that the acoustic system (4) reads equivalently

∂ t τ -τ ∂ x u = 0, ∂ t u + τ ∂ x p = 0, ∂ t v = 0, ∂ t E + τ ∂ x (pu) = 0. ( 6
)
The acoustic system ( 6) is a quasilinear system that can be simply checked to be strictly hyperbolic. Indeed, the Jacobian of the system (6) has three eigenvalues (λ 1 , λ 2 , λ 3 ) = (-c, 0, +c). The waves associated with λ 1 and λ 3 are genuinely nonlinear waves while the wave of velocity λ 2 = 0 is a stationnary contact discontinuity.

In order to derive an update process from (ρ, ρu, ρv, ρE) n j to (ρ, ρu, ρv, ρE) n+1- j , we will perform several approximations. We notice that for a smooth solution [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF] we also have ∂ t p + τ (ρc) 2 ∂ x u = 0 and we thus choose to perform a Suliciu-type approximation of (6) for t ∈ [t n , t n + ∆t) by introducing a surrogate pressure Π and considering the relaxed system

               ∂ t τ -τ ∂ x u = 0, ∂ t u + τ ∂ x Π = 0, ∂ t v = 0, ∂ t E + τ ∂ x (Πu) = 0, ∂ t Π + τ a 2 ∂ x u = ν(Π -p), (7a) (7b) (7c) (7d) (7e)
where a > 0 is parameter whose choice will be specified later. In the regime ν → +∞ we formally recover [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF]. In our numerical solver context, we classically mimic the ν → +∞ regime enforcing at each time step Π n j = p EOS (τ n j , e n j ) and then solving [START_REF] Cordier | An Asymptotic-Preserving all-speed scheme for the Euler and NavierStokes equations[END_REF] with ν = 0. At last, for t ∈ [t n , t n + ∆t) we choose to approximate τ (x, t)∂ x by τ (x, t n )∂ x in [START_REF] Cordier | An Asymptotic-Preserving all-speed scheme for the Euler and NavierStokes equations[END_REF]. If one introduces the mass variable m defined by dm = ρ(x, t n )dx our approximation of (6) (up to an abuse of notation) can be expressed in the following fully conservative form

∂ t W + ∂ m F(W) = 0, (8)
where W = (τ, u, v, E, Π) T and F(W) = (-u, Π, 0, Πu, a 2 u) T . Let us remark that ( 8) is consistent with a Suliciu relaxation of the gas dynamics equation written in Lagrangian coordinates using a mass variable formulation. The solution of the Riemann problem associated with ( 8) can be derived explicitly (see section C). This allows to write an exact Godunov solver for (8) that turns out to be an approximate Riemann solver for ( 6) following the Harten-Lax-van Leer formalism (see section B and [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF]). It provides us with the update formula

         W n+1- j = W n j - ∆t ∆x F j+1/2 -F j-1/2 , F j+1/2 = F(W n j , W n j+1 ), F(W L , W R ) = (-u * , Π * , 0, Π * u * , a 2 u * ) T , (9a) (9b) (9c) where        u * = (u R + u L ) 2 - 1 2a (Π R -Π L ), Π * = (Π R + Π L ) 2 - a 2 (u R -u L ). (10a) (10b)
The update of the conservative variables is obtained by setting

ρ n+1- j = 1/τ n+1- j , (ρu) n+1- j = ρ n+1- j × u n+1- j , (ρv) n+1- j = ρ n+1- j × v n+1- j and (ρE) n+1- j = ρ n+1- j × E n+1- j
. This can be summed up by the following update formulas

                           L j ρ n+1- j = ρ n j , L j (ρu) n+1- j = (ρu) n j - ∆t ∆x (Π * j+1/2 -Π * j-1/2 ), L j (ρv) n+1- j = (ρv) n j , L j (ρE) n+1- j = (ρE) n j - ∆t ∆x (Π * j+1/2 u * j+1/2 -Π * j-1/2 u * j-1/2 ), L j = 1 + ∆t ∆x (u * j+1/2 -u * j-1/2 ). ( 11a
) (11b) (11c) (11d) (11e)
Let us remark that (9) also proposes an update relation for Π. However in this case Π is just a disposable intermediate value whose role only consists in providing a formula for the interface pressure terms and the udpated value Π n+1- k will be discarded. Indeed in this explicit scheme, Π is updated after each time step by the equilibrium formula Π n j = p EOS (τ n j , e n j ). However, this will no longer be the case for semi-implicit strategy as we will see in section 5.4.

Let us finally note that the relaxation scheme ( 9) is equivalent to the acoustic scheme [START_REF] Després | Inégalité entropique pour un solveur conservatif du système de la dynamique des gaz en coordonnées de lagrange[END_REF]. In order to avoid numerical instabilities, the parameter a must complies with the subcharacteristic condition a > max ρc, [START_REF] Després | Inégalité entropique pour un solveur conservatif du système de la dynamique des gaz en coordonnées de lagrange[END_REF] for all possible values of ρc when considering a solution of the equilibrium system [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF]. In practice we will choose a value a LR for each interface by setting

a LR = K max(ρ n L c n L , ρ n R c n R ), (13)
where K ≥ 1, LR = j + 1/2, L = j and R = j + 1. We refer the reader to [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF][START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF][START_REF] Chalons | Navier-stokes equations with several independant pressure laws and explicit predictor-corrector schemes[END_REF][START_REF] Després | Lois de conservations Eulériennes, Lagrangiennes et méthodes numériques[END_REF] and the reference therein for more details.

Properties and approximation of the one-dimensional transport system

The transport system equation discretization is quite simple. Indeed, system ( 5) is a quasi-linear hyperbolic system that only involves the transport of the conservative variables with the velocity u. We choose to approximate the solution of (5) thanks to a standard upwind Finite-Volume approximation for ϕ ∈ {ρ, ρu, ρv, ρE}

ϕ n+1 j = ϕ n+1- j - ∆t ∆x u * j+1/2 ϕ n+1- j+1/2 -u * j-1/2 ϕ n+1- j-1/2 + ∆t ∆x ϕ n+1- j u * j+1/2 -u * j-1/2 , (14)
where

ϕ n+1- j+1/2 = ϕ n+1- j , if u * j+1/2 ≥ 0, ϕ n+1- j+1 , if u * j+1/2 < 0.
Let us finally remark that ( 14) can be recast into

ϕ n+1 j = ϕ n+1- j L j + ∆t ∆x u * j+1/2 ϕ n+1- j+1/2 -u * j-1/2 ϕ n+1- j-1/2 . ( 15
)

Properties of the operator splitting scheme

We present here a few properties of the operator splitting scheme defined by ( 9) and [START_REF] Dauvergne | Renormalization of the numerical diffusion for an upwind finite volume method. application to the simulation of Kelvin-Helmholtz instability. Finite volumes for complex applications[END_REF]. Let us first remark that this algorithm performs the same update as a classical Lagrange-Remap (or equivalently Lagrange-Projection) algorithm for one-dimensional problems (see appendix A) although the design of our algorithm does not involve a moving mesh for following the variables in a Lagrangian reference frame. This feature will be the key element of the multi-dimensional extension of the present scheme. It is also interesting to mention that the operator splitting strategy also provided a mean of treating the waves of the gas dynamics system (3) separately: the acoustic step only involves acoustic waves while freezing the transport waves. The transport step only deals with the contact discontinuity of the material transport. Let us mention that a similar operator splitting was used in [?]. The overall update from variable at instant t n to variables at instant t n+1 is fully conservative with respect to ρ, ρu, ρv and ρE. Indeed, we have

                       ρ n+1 j = ρ n j + ∆t ∆x u * j+1/2 ρ n+1- j+1/2 -u * j-1/2 ρ n+1- j-1/2 , (ρu) n+1 j = (ρu) n j + ∆t ∆x u * j+1/2 (ρu) n+1- j+1/2 + Π * j+1/2 -u * j-1/2 (ρu) n+1- j-1/2 -Π * j-1/2 , (ρv) n+1 j = (ρv) n j + ∆t ∆x u * j+1/2 (ρv) n+1- j+1/2 -u * j-1/2 (ρv) n+1- j-1/2 -, (ρE) n+1 j = (ρE) n j + ∆t ∆x u * j+1/2 (ρE) n+1- j+1/2 + Π * j+1/2 u * j+1/2 -u * j-1/2 (ρE) n+1- j-1/2 -Π * j-1/2 u * j-1/2 . ( 16a
) (16b) (16c) (16d)
The scheme ( 9)- [START_REF] Dellacherie | Construction of modified Godunov type schemes accurate at any Mach number for the compressible Euler system[END_REF] for the acoustic step is stable under the Courant-Friedrichs-Lewy (CFL) condition

∆t ∆x max j∈Z max(τ n j , τ n j+1 )a j+1/2 ≤ 1 2 . ( 17
)
If one notes b ± = b±|b| 2 , then a classical result states that the CFL condition associated with the transport scheme [START_REF] Dauvergne | Renormalization of the numerical diffusion for an upwind finite volume method. application to the simulation of Kelvin-Helmholtz instability. Finite volumes for complex applications[END_REF]

reads ∆t max j∈Z (u * j-1 2 ) + -(u * j+ 1 2 ) -< ∆x. ( 18
)
Entropy-related stability properties of the scheme will be examined in section 5.3. One can also remark that both the acoustic steps and the transport steps are achieved thanks to genuine Godunov solvers applied to simplified subsystems. [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] Behavior of the scheme with respect to the Mach regime

We are now interested in the behavior of the numerical scheme with respect to the variations of the Mach regime. In order to characterize the Mach regime of the flow, we consider a classical rescaling of the equations (3): let us introduce the following non-dimensional quantities:

x = x L , t = t T , ρ = ρ ρ 0 , ũ = u u 0 , ṽ = v v 0 , ẽ = e e 0 , p = p p 0 , c = c c 0 . ( 19
)
The parameters L, T , u 0 = v 0 = L T , ρ 0 , e 0 = p 0 ρ 0 , p 0 and c 0 = p0 ρ0 denote respectively a characteristic length, time, velocity, density, internal energy, pressure and sound speed. If M = u0 c0 is the so-called Mach-number then system (3) reads

             ∂ t ρ + ∂ x( ρũ) = 0, ∂ t( ρũ) + ∂ x( ρũ 2 ) + 1 M 2 ∂ x p = 0, ∂ t( ρṽ) + ∂ x( ρũṽ) = 0, ∂ t( ρ Ẽ) + ∂ x[( ρ Ẽ + p)ũ] = 0, (20a) (20b) (20c) (20d)
where Ẽ = ẽ + 1 2 M 2 ũ2 . For a given small value of the Mach number, we distinguish two cases: • the term ∂ x p remains of magnitude O(M 2 ). Then the variations of ρũ are of order 1 which implies that all the tilde variables will remain of order 1. We shall refer this case as the low Mach regime;

• the term ∂ x p does not remain of magnitude O(M 2 ). Then the variations of ρũ will reach a magnitude O(1/M ) or O(1/M 2 ). These large magnitude variations of the momentum will induce a growth of the Mach number and thus a change of Mach regime.

Before going any further, let us underline that in the present approach we do not intend to study the behavior of the rescaled system (3) in the limit regime M → 0. This delicate question has been widely investigated over the past years and is still a rich field of research [START_REF] Guillard | On the behavior of upwind schemes in the low Mach limit[END_REF][START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible euler system at low Mach number[END_REF][START_REF] Liu | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF]. We focus here on a simpler task that consists in examining the consistency of a rescaled approximate solution provided by the splitting operator algorithm with the solution of (20) in the low Mach regime. The framework we will place ourselves in does not require sophisticated hypotheses and may deal with the evaluation of a local behavior of the solution (a few neighbouring cells in the discrete setting). More precisely, if one considers smooth solutions of ( 20) and considers the truncation error of the rescaled numerical scheme in the sense of Finite Difference, how does it depends on M in the low Mach regime?

Introducing the rescaling defined earlier into [START_REF] Dellacherie | Construction of modified Godunov type schemes accurate at any Mach number for the compressible Euler system[END_REF] we get

ũ * j+1/2 = 1 2 (ũ n j + ũn j+1 ) - 1 2ã j+1/2 M ( Πn j+1 -Πn j ), Π * j+1/2 = 1 2 ( Πn j + Πn j+1 ) - ãj+1/2 M 2 (ũ n j+1 -ũn j ),
for [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible euler system at low Mach number[END_REF] we have

                             Lj ρn+1- j = ρn j , Lj (ρũ) n+1- j = (ρũ) n j - ∆ t M 2 ∆x ( Π * j+1/2 -Π * j-1/2 ), Lj (ρṽ) n+1- j = (ρṽ) n j , Lj (ρ Ẽ) n+1- j = (ρ Ẽ) n j - ∆ t ∆x ( Π * j+1/2 ũ * j+1/2 -Π * j-1/2 ũ * j-1/2 ), Lj = L j = 1 + ∆ t ∆x (ũ * j+1/2 -ũ * j-1/2 ), (21a) (21b) (21c) (21d) (21e)
and finally if φ ∈ {ρ, ρũ, ρṽ, ρ Ẽ} the rescaling of ( 15), reads

1 ∆ t ( φn+1 j -Lj φn+1- j ) + 1 ∆x φn+1- j+1/2 ũ * j+1/2 -ũ * j-1/2 φn+1- j-1/2 = 0. ( 22
)
Note that the CFL restriction of the acoustic step reads now

∆ t ∆x max(τ n j , τ n j+1 )ã n j+1/2 ≤ M 2 , (23)
while the CFL restriction associated with the transport step is

(ũ * j-1/2 ) + -(ũ * j+1/2 ) -∆ t ∆x ≤ 1. ( 24
)
In order to evaluate the truncation error (in the Finite Difference sense) in the low Mach regime, we use the classical tool of equivalent equations. Let (x, t) → b be a parameter of (rescaled) functions that describe a smooth flow. With a classical slight abuse of notation, we consider that φ(x j , t n ) = φn j so that we can substitute these functions into the discrete update formula when φ ∈ {ρ, ũ, ṽ, Ẽ, Π}. We suppose that we are in low Mach regime, namely ∂ x p = O(M 2 ). This hypothesis yields that Πn j+1 = Πn j + O(M 2 ∆x) for the discrete unknowns. We have the following result.

Proposition 1. In the low Mach regime, the rescaled discretization of the acoustic step is consistent with

∂ t τ -τ ∂ x ũ = O(∆ t) + O(M ∆x), ∂ t ũ + τ M 2 ∂ x p = O(∆ t) + O ∆x M , ∂ t ṽ = O(∆ t), ∂ t Ẽ + τ ∂ x( pũ) = O(∆ t) + O(M ∆x).
The rescaled discretization of the transport step is consistent with

∂ t φ + ũ∂ x φ = O(∆ t) + O(∆x) + O(M ∆x),
and the equivalent equation verified by the rescaled scheme reads

               ∂ t ρ + ∂ x( ρũ) = O(∆ t) + O(∆x) + O(M ∆x), ∂ t( ρũ) + ∂ x( ρũ 2 ) + 1 M 2 ∂ x p = O(∆ t) + O(∆x) + O(M ∆x) + O ∆x M , ∂ t( ρṽ) + ∂ x( ρũṽ) = O(∆ t) + O(∆x) + O(M ∆x), ∂ t( ρ Ẽ) + ∂ x[( ρ Ẽ + p)ũ] = O(∆ t) + O(∆x) + O(M ∆x). ( 25a
) (25b) (25c) (25d)
Proof. There exists three smooth functions A, B and C of magnitude 1 with respect to M such that

ũ * j+1/2 = ũn j+1 + ũn j 2 + M ∆xA(x j+1/2 , t n ) + O(M ∆x 2 ), Π * j+1/2 = Πn j+1 + Πn j 2 + M ∆xB(x i+1/2 , t n ) + O(M ∆x 2 ), Π * j+1/2 ũ * j+1/2 = (ũ n j+1 + ũn j )( Πn j+1 + Πn j ) 4 + M ∆xC(x i+1/2 , t n ) + O(M ∆x 2 ).
Injecting the above relation into [START_REF] Klein | Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: One-dimensional flow[END_REF] we get

                               Lj ρn+1- j = ρn j , Lj (ρũ) n+1- j = (ρũ) n j - ∆ t M 2 Πn j+1 -Πn j-1 2∆x + O ∆x∆ t M , Lj (ρ j ṽ) n+1-= (ρṽ) n j , L j (ρ Ẽ) n+1- j = (ρ Ẽ) n j -∆ t (ũ n j+1 + ũn j )( Πn j+1 + Πn j ) 4∆x - (ũ n j-1 + ũn j )( Πn j-1 + Πn j ) 4∆x + O(M ∆x∆ t), Lj = 1 + ∆ t ũn j+1 -ũn j-1 2∆x + O(M ∆x∆ t). (26a) (26b) (26c) (26d) (26e) This yields                        Lj ρn+1- j = ρn j , Lj (ρũ) n+1- j = (ρũ) n j - ∆ t M 2 ∂ x p + O ∆x∆ t M + O ∆x 2 ∆ t , Lj (ρ j ṽ) n+1-= (ρṽ) n j , L j (ρ Ẽ) n+1- j = (ρ Ẽ) n j -∆ t∂ x( pũ) + O(M ∆x∆ t) + O(∆x 2 ∆ t), Lj = 1 + ∆ t∂ x ũ + O(M ∆x∆ t) + O(∆x 2 ∆ t). (27a) (27b) (27c) (27d) (27e)
Remark 1. For smooth solutions in the low Mach regime, we have

∂ xx p = O(M 2 ). We used this relation to obtain the term O ∆x 2 ∆ t in (27b).
Let us remark that ( 27) is indeed consistent at order 1 with respect to ∆x with

∂ t τ -τ ∂ x ũ = O(∆ t) + O(M ∆x), ∂ t ũ + τ M 2 ∂ x p = O(∆ t) + O ∆x M , ∂ t ṽ = O(∆ t), ∂ t Ẽ + τ ∂ x( pũ) = O(∆ t) + O(M ∆x).
Now we turn to the transport step. Accounting for the low Mach hypothesis, [START_REF] Lax | Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[END_REF] becomes

1 ∆ t ( φn+1 j -Lj φn+1- j ) + 1 2∆x φn+1- j+1/2 (ũ n j+1 + ũn j ) -φn+1- j-1/2 (ũ n j + ũn j-1 ) = O(M ∆x), hence 1 ∆ t ( φn+1 j -Lj φn+1- j ) + ∂ x( φũ) = O(∆x) + O(M ∆x), (28)
which is consistent with

∂ t φ + ũ∂ x φ = O(∆ t) + O(∆x) + O(M ∆x).
Finally, using ( 27) into [START_REF] Thornber | An improved reconstruction method for compressible flows with low Mach number features[END_REF] we finally obtain the desired result.

Remark 2. It is important to note that the analysis we proposed in this section cannot be considered as an exhaustive explanation for the behavior of the numerical scheme in the Low Mach regime. It just merely provides magnitude estimate of the truncation error. Considering the same lines with additional hypotheses: ρ, ũ, ṽ, Ẽ are solution of the rescaled gas dynamics equations in the low Mach regime with well-prepared conditions [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible euler system at low Mach number[END_REF], then one can show that the O(∆x/M ) term in (25b) does vanish [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible euler system at low Mach number[END_REF] for one-dimensional problems set over the whole real line. The analysis is delicate and depends on many hypotheses: for two-dimensional problems same results can be obtained for discretization over a triangular mesh with periodic boundary conditions. However, this no longer works for two-dimensional Cartesian meshes where the classical Godunov-type solvers perform poorly with periodic boundary conditions. More general boundary conditions require a specific study for each case [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible euler system at low Mach number[END_REF][START_REF] Dellacherie | Construction of modified Godunov type schemes accurate at any Mach number for the compressible Euler system[END_REF].

Low Mach correction

The equivalent equation (25) satisfied by the rescaled scheme is clearly not satisfactory because of the term O( ∆x M ) which behaves bad when M ∆x. This suggests to modify the scheme accordingly.

Correction of the low Mach behavior: a simple flux modification

In the light of the previous asymptotic analysis, we propose to leave the projection step unchanged and rather focus on the acoustic step of the scheme. In the acoustic step, we suggest to simply replace Π * j+1/2 by

Π * ,θ j+1/2 = 1 2 (Π n j + Π n j+1 ) -θ j+1/2 a j+1/2 2 (u n j+1 -u n j ). ( 29
)
The associated dimensionless flux reads

Π * ,θ j+1/2 = 1 2 ( Πn j + Πn j+1 ) -θ j+1/2 ãj+1/2 M 2 (ũ n j+1 -ũn j ). ( 30
)
This yields the following modified scheme for the acoustic step.

         W n+1- j = W n j - ∆t ∆x F j+1/2 -F j-1/2 , F j+1/2 = F θ (W n j , W n j+1 ), F θ (W L , W R ) = (-u * , Π * ,θ , 0, Π * ,θ u * , a 2 u * ) T . (31a) (31b) (31c)
Let us underline that this modification solely alters the non-centered terms of the pressure flux. In other words this does not modify the ultimate consistency of Π * ,θ j+1/2 with the pressure value, it does impact the numerical dissipation involved with the discretization of the pressure terms. This approach complies with several previous works that have been investigating the approximation of the low Mach regime like [START_REF] Toumi | Approximate Riemann solvers and flux vector splitting schemes for two-phase flow[END_REF][START_REF] Paillère | Comparison of low mach number models for natural convection problems[END_REF][START_REF] Dauvergne | Renormalization of the numerical diffusion for an upwind finite volume method. application to the simulation of Kelvin-Helmholtz instability. Finite volumes for complex applications[END_REF]. While such modification is usually delicate with regards to the stability of the numerical scheme, we will nevertheless see that the resulting modified numerical scheme is still endowed with stability properties (see section 5.3). In the sequel, in order to perform an equivalent equation analysis with the modified pressure flux, we consider a smooth function x → θ such that θ j+1/2 = θ(x j+1/2 ). We have the following consistency properties for the numerical scheme with the modified pressure flux Π * ,θ j+1/2 . Proposition 2. In the low Mach regime, the rescaled discretization (31) of the acoustic step is consistent with

∂ t τ -τ ∂ x ũ = O(∆ t) + O(M ∆x), ∂ t ũ + τ M 2 ∂ x p = O(∆ t) + O θ∆x M , ∂ t ṽ = O(∆ t), ∂ t Ẽ + τ ∂ x( pũ) = O(∆ t) + O(M ∆x) + O(M θ∆x).
The rescaled discretization of the transport step is consistent with

∂ t φ + ũ∂ x φ = +O(∆ t) + O(∆x) + O(M ∆x),
and the equivalent equation verified by the rescaled scheme reads

               ∂ t ρ + ∂ x( ρũ) = O(∆ t) + O(∆x) + O(M ∆x), ∂ t( ρũ) + ∂ x( ρũ 2 ) 1 M 2 + ∂ x p = O(∆ t) + O(∆x) + O θ∆x M , ∂ t( ρṽ) + ∂ x( ρũṽ) = O(∆ t) + O(∆x) + O(M ∆x). ∂ t( ρ Ẽ) + ∂ x[( ρ Ẽ + p)ũ] = O(∆ t) + O(∆x) + O(M ∆x) + O(M θ∆x). (32a) (32b) (32c) (32d)
As a consequence, provided that we impose the asymptotic behavior θ j+1/2 = O(M ), the truncation error is uniform with respect to M .

Proof. Following similar lines as in the proof of proposition 1 and using the same notations, there exists three smooth functions A, B, C and D of magnitude 1 with respect to M such that ũ * j+1/2 = ũn j+1 + ũn

j 2 + M ∆xA(x j+1/2 , t n ) + O(M ∆x 2 ), Π * ,θ j+1/2 = pn j+1 + pn j 2 + θ j+1/2 M ∆xB(x i+1/2 , t n ) + O(M ∆x 2 ), Π * ,θ j+1/2 ũ * j+1/2 = (ũ n j+1 + ũn j )(p n j+1 + pn j ) 4 + M ∆xC(x i+1/2 , t n ) + M θ j+1/2 ∆xD(x i+1/2 , t n ) + O(M ∆x 2 ).
The rest of the analysis follows the same line as the proof of proposition 1. Using [START_REF] Klein | Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: One-dimensional flow[END_REF] we get

                       Lj ρn+1- j = ρn j , Lj (ρũ) n+1- j = (ρũ) n j - ∆ t M 2 ∂ x p + O θ∆x∆ t M + O ∆x 2 ∆ t , Lj (ρ j ṽ) n+1-= (ρṽ) n j , L j (ρ Ẽ) n+1- j = (ρ Ẽ) n j -∆ t∂ x( pũ) + O(M ∆x∆ t) + O(M θ∆x∆ t) + O(∆x 2 ∆ t), Lj = 1 + ∆ t∂ x ũ + O(M ∆x∆ t) + O(∆x 2 ∆ t). (33a) (33b) (33c) (33d) (33e) and (22) yields again 1 ∆ t ( φn+1 j -Lj φn+1- j ) + ∂ x( φũ) = O(∆x) + O(M ∆x). ( 34
)
Relations ( 33) and (34) provides the desired results.

Remark 3. In the light of the truncation error that appears in [START_REF] Schochet | Fast singular limits of hyperbolic PDEs[END_REF], one can see that it is not necessary to involve a correction for the energy flux term in (31c). It would be possible to consider a numerical scheme with the definition (10a) for the velocity at the interface, the modified pressure (29) for interface pressure terms and Π * u * for the energy flux.

Approximate Riemann solver for the modified acoustic scheme

The modified numerical scheme (31) for the acoustic step belongs to the category of flux-based solver. Indeed, this solver relies on an update formula (31a) that involves the modified flux (31c). We will prove in this section that this modified flux solver can also be obtained thanks to an approximate Riemann solver in the sense of Harten, Lax and van Leer [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF], see also Annex B for a quick refresh on this, that is consistent with the integral form of (8). This formalism is useful to establish stability properties. We have the following proposition.

Proposition 3. There exists a simple approximate Riemann solver that is an approximation of the Riemann problem associated with the relaxed acoustic problem (8) and whose associated flux matches the flux of the modified acoustic solver. More precisely, there exists a self-similar function

W θ RP m t ; W L , W R = (τ, u, v, E, Π) m t ; W L , W R =          W L , if m/t < -a, W * ,θ L , if -a ≤ m/t < 0, W * ,θ R , if 0 ≤ m/t < +a, W R , if +a ≤ m/t. (35) such that F θ (W R , W L ) = F(W L ) - 0 -∞ [W θ RP (ξ; W L , W R ) -W L ] dξ = F(W R ) + +∞ 0 [W θ RP (ξ; W L , W R ) -W R ] dξ. ( 36
)
The states

W * ,θ L = (τ * ,θ L , u * ,θ L , v * ,θ L , Π * ,θ L ) T and W * ,θ R = (τ * ,θ R , u * ,θ R , v * ,θ R , Π * ,θ R ) T are given by τ * ,θ L = τ L + 1 a (u * -u L ), τ * ,θ R = τ R + 1 a (u R -u * ), u * ,θ L = u * + 1 2 (θ -1)(u R -u L ), u * ,θ R = u * + 1 2 (1 -θ)(u R -u L ), v * ,θ L = v L , v * ,θ R = v R , E * ,θ L = E L + 1 a (Π L u L -Π * ,θ u * ), E * ,θ R = E R + 1 a (Π * ,θ u * -Π R u R ) Π * ,θ L = Π * , Π * ,θ R = Π * . ( 37a
) (37b) (37c) (37d) (37e)
Proof. Suppose that W θ RP is consistent with the integral form of the relaxed acoustic problem (8) then for a given W L and W R we have

F(W R ) -F(W L ) = -a(W * ,θ L -W L ) + a(W R -W * ,θ R ), which reads W * ,θ R + W * ,θ L = W R + W L - 1 a (F(W R ) -F(W L )). ( 38
)
If the resulting flux of this approximate Riemann solver is F θ (W L , W R ) then ( 36) is verified and yields

2F θ (W L , W R ) = F(W R ) + F(W L ) -a(W * ,θ L -W L ) -a(W R -W * ,θ R ) or equivalently W * ,θ R -W * ,θ L = W R -W L + 1 a 2F θ (W L , W R ) -F(W L ) -F(W R ) . (39)
Both ( 38) and (39) provide

W * ,θ L = W L - 1 a (F θ (W L , W R ) -F(W L )), W * ,θ R = W R + 1 a (F θ (W L , W R ) -F(W R )).
This yields the desired results.

Using this approximate Riemann solver, we can deduce that the modified acoustic solver (31) is stable under the same CFL conditions (17) that does not depend on the Mach number M . Moreover, when θ = 1 the self-similar function W θ RP defined in proposition 3 degenerates to the exact solution of the Riemann problem associated with relaxed acoustic system [START_REF] Degond | All speed method for the Euler equation in the low mach number limit[END_REF].

Finally, if one takes into account the equilibrium projection step of the relaxation strategy into the approximate Riemann solver of proposition 3, we have Π L = p EOS (τ L , e L ), and Π R = p EOS (τ R , e R ). Under this assumption, it is easy to check that the first coordinates (τ, u, v, E) of the self similar function W θ RP are consistent with the integral form of the acoustic system (6).

Properties of the modified operator splitting scheme

We start this section by examining the ability of the modified operator splitting scheme to satisfies a discrete entropy inequality. In the sequel, I(b, b ) ⊂ R will denote the interval whose bounds are b ∈ R and b ∈ R. We consider the following slightly more restrictive subcharacteristic condition

τ * L > 0, -∂ τ p EOS (τ, s L ) ≤ a 2 , ∀τ ∈ I(τ L , τ * L ), τ * R > 0, -∂ τ p EOS (τ, s R ) ≤ a 2 , ∀τ ∈ I(τ R , τ * R ), (40)
and we start with the two following technical results. we also refer the reader to Annex B for a quick refresh on this topic.

Lemma 1. Consider the solution of Riemann problem for the relaxed acoustic system [START_REF] Degond | All speed method for the Euler equation in the low mach number limit[END_REF]. Suppose that (40)

is verified. Let s k = s EOS (τ k , e k ), k = L, R, we have e * k -e EOS (τ * k , s k ) - p EOS (τ * k , s k ) -Π * 2 2a 2 ≥ 0. ( 41
)
Proof. We consider the case k = R and set for τ

∈ I(τ R , τ * R ) φ(τ ) = e EOS (τ, s R ) - p EOS (τ, s R ) 2 2a 2 -e EOS (τ * R , s R ) + p EOS (τ * R , s R ) 2 2a 2 + p EOS (τ * R , s R ) τ + p EOS (τ, s R ) a 2 -τ * R - p EOS (τ * R , s R ) a 2 .
We have

φ (τ ) = p EOS (τ, s R ) -p EOS (τ * R , s R ) 1 -ρ 2 c 2 (τ, s R )/a 2 . If τ R > τ > τ * R (resp. τ R < τ < τ * R ) the Weyl assumptions (2) provides p EOS (τ, s R ) -p EOS (τ * R , s R ) < 0 (resp. p EOS (τ, s R ) -p EOS (τ * R , s R ) > 0
) and together with hypothesis (40) this yields φ (τ ) ≥ 0 (resp. φ (τ ) ≤ 0). As φ(τ * R ) = 0 we obtain that

φ(τ R ) > φ(τ * R ) = 0 for τ ∈ I(τ R , τ * R ). Using the Riemann invariant jump relation (e * R -Π * 2a 2 ) = (e R -Π R 2a 2 ), one obtains 0 < φ(τ R ) = e * R -e EOS (τ * R , s R ) -1 2a 2 (p EOS (τ * R , s R ) -Π * ) 2 .
The same lines applies for the case k = L.

= E * ,θ k -(u * ,θ k ) 2 /2 for k = L, R then we have e * ,θ k -e EOS (τ * ,θ k , s k ) - 1 2a 2 p EOS (τ * ,θ k , s k ) -Π * 2 + (1 -θ) 2 (u R -u L ) 2 8 ≥ 0, k = L, R. (42)
Proof. One has u * ,θ

R = u * + (1 -θ)(u R -u L )/2 and Π * ,θ = Π * + (1 -θ)a(u R -u L )/2 and together with (37) one obtains e * ,θ R = e * R -(1 -θ) 2 (u R -u L ) 2 /8.
Injecting this relation into (41) and noticing that τ * ,θ R = τ * R provides the desired result for k = R. The case k = L is obtained with the same lines.

It is now clear that the inequalities

- 1 2a 2 p EOS (τ * ,θ k , s R ) -Π * 2 + (1 -θ) 2 (u R -u L ) 2 8 ≤ 0, k = L, R (43)
can help us equip the modified numerical scheme with a discrete entropy inequality.

Proposition 4. Let s * ,θ k = s EOS (τ * ,θ k , e * ,θ k ) for k = L, R. If assumption (43) is verified, we have 0 ≤ -a(s * ,θ L -s L ) + a(s R -s * ,θ R ). ( 44
)
Inequality (44) implies that the modified scheme (31) for the acoustic step is consistent with the integral form of the entropy inequality

∂ t s(τ, e) ≤ 0. ( 45
)
Moreover, the explicit modified scheme (31) is equipped with a discrete entropy inequality. Indeed there exists a numerical flux function q n j+1/2 = q(W n j , W n j+1 ) that is consistent with 0 when ∆t and ∆x tend to 0 such that

s(τ n+1- j , e n+1- j ) -s(τ n j , e n j ) + τ n j ∆t ∆x (q n j+1/2 -q n j-1/2 ) ≤ 0. (46)
Proof. Let k = L, R, under hypothesis (43), we have that e * ,θ

k ≥ e EOS (τ * ,θ k , s k ). According to (2) → s EOS (τ * ,θ k , ) is increasing, thus s EOS (τ * ,θ k , e * ,θ k ) = s * ,θ k ≥ s EOS (τ * ,θ k , e EOS (τ * ,θ k , s k )) = s k .
Inequality (44) follows trivially. Relation (44) expresses the consistency with the integral form of (45) and it provides the entropy inequality (46) (see [1, Chap. 2] and Annex B).

We can now state the following entropic property for the full modified operator splitting explicit scheme composed by (31) and ( 14). Proposition 5. If the assumptions (43), ( 17) and (18) are verified, then the explicit scheme defined by (31) and (14) verifies the following discrete entropy inequality

ρ n+1 j s(τ n+1 j , e n+1 j ) -ρ n j s(τ n j , e n j ) + ∆t ∆x g n j+1/2 -g n j-1/2 ≤ 0, (47)
where the numerical entropy flux is defined by

g n j+1/2 = (u * j+1/2 ) + ρ n+1- j s(τ n+1- j , e n+1- j ) + (u * j+1/2 ) -ρ n+1- j s(τ n+1- j+1 , e n+1- j+1 ) + q n j+1/2 . ( 48
)
Proof. Let φ ∈ (ρ, ρu, ρv, ρE), under the CFL assumption (18) the transport scheme ( 14) expresses φ n+1 j as a convex combination of φ n+1- i , i = j -1, j, j + 1, indeed one has

φ n+1 j = ∆t ∆x (u * j-1/2 ) + φ n+1- j-1 + 1 - ∆t ∆x ((u * j+1/2 ) --(u * j-1/2 ) + ) φ n+1- j + ∆t ∆x (u * j+1/2 ) -φ n+1- j+1 .
As the mapping (ρ, ρu, ρv, ρE) → -(ρs)(τ, e) is a strictly convex function (see for example [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]) we obtain that

-(ρs)(τ n+1 j , e n+1 j ) ≤ - ∆t ∆x (u * j-1/2 ) + (ρs)(τ n+1- j-1 , e n+1- j-1 ) -1 - ∆t ∆x ((u * j+1/2 ) --(u * j-1/2 ) + ) (ρs)(τ n+1- j , e n+1- j ) - ∆t ∆x (u * j+1/2 ) -(ρs)(τ n+1- j-1 , e n+1- j+1 ).
Using relation (46) one obtains (47).

We now sum up the main properties of the modified operator splitting scheme.

Theorem 1. Suppose that (17), (18) (12) are satisfied, the explicit scheme defined by (31) and (14) verifies 1. the scheme is conservative with respect to the density ρ, the momentum ρu and total energy ρE, 2. the density ρ n j is positive for all j and n > 0 provided that ρ 0 j is positive for all j, 3. if θ = O(M ), then the truncation error of the numerical scheme is uniform with respect to M < 1, 4. if (43) is verified then the numerical scheme is equipped with a discrete entropy inequality, 5. if (43) is verified then e n j > 0 for all j ∈ Z and all n ∈ N.

It is clear from (32b) that the choice θ = O(M ) is natural for the modified scheme to have an equivalent equation which is satisfactory when M ∆x (uniform consistency w.r.t. M ). At this stage θ = O(M ) is not made precise, see section 6 below. Let us now discuss the new condition which is related to the correction θ.

Behavior of condition (43) in the low-Mach regime for a perfect gas equation of state

We have just seen that the scheme is entropic provided that (43) is satisfied. In this section, we study the compatibility in the low Mach regime between the condition (43) that is required to obtain a discrete entropy inequality and the condition θ = O(M ) that is required to have uniform consistency with respect to M (see section 5). If |u R -u L | = 0, any value of θ ∈ R verifies condition (43), we can then assume that |u R -u L | > 0. We consider the case of a Perfect Gas EOS defined by p EOS (ρ, e) = (γ -1)ρe, where γ is the specific heat ratio. First, let us recall that τ * ,θ

k = τ * k and Π k = p EOS (τ k , s k ), k = R, L. For k = R, relation (43) reads |1 -θ| ≤ 2 a |p EOS (τ * R , s R ) -Π * | |u R -u L | . (49)
Let us remark that the right hand side of this inequality does not depend on θ. The Perfect Gas assumption provides that p EOS (τ * R , s R ) = Π R (τ R /τ * R ) γ , therefore thanks to the definition of Π * we get

p EOS (τ * R , s R ) -Π * = Π R (τ R /τ * R ) γ - Π L + Π R 2 + a 2 (u R -u L ). ( 50
)
The definition of τ * ,θ R = τ * R by (37) using the dimensionless parameters defined by [START_REF] Liu | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF] gives τ * R = τR + ( ΠR -ΠL )/(2ã 2 ) + M (ũ R -ũL )/(2ã). If one now supposes that the flow is locally in the low Mach regime, then we have

∂ x Π = O(M 2 ), therefore Π R -Π L = O(M 2 ∆x). Thus we obtain τ * R τR = 1 + M ũR -ũL 2ãτ R + O(M 2 ∆x).
Injecting the above relation into (50), we obtain

p EOS (τ * R , s R ) -Π * p 0 = - M 2 1 - γ ΠR ã2 τR ã(ũ R -ũL ) + O(M 2 ∆x).
Using the fact that γp EOS (τ R , s R ) = γΠ R = ρ R (c R ) 2 for a Perfect Gas in the previous relation allows to recast (49) into

|1 -θ| ≤ 1 - ρR cR ã 2 + O M ∆x |ũ R -ũL | . ( 51
)
Let us recall that by definition : ã = K max(ρ R cR , ρL cL ) with K ≥ 1. Suppose without loss of generality that ρR cR = max(ρ R cR , ρL cL ) then (ρ R cR ) 2 /ã 2 = 1/K and the condition (51) becomes Behavior when M → 0. When M → 0 the above inequality yields that θ ≥ (1/K) 2 if one wants to enforce uniform consistency with respect to M by setting θ = O(M ). tThis leads to a contradiction. As a conclusion, a correction scheme with θ = O(M ) does not provide an entropic scheme in the asymptotic limit M → 0. On the contrary, θ = 1 which correspond to the classic unmodified scheme is still entropic. Nevertheless, it is reasonable to consider that in the limit M → 0, the solution of the gas dynamics equation is smooth and therefore the consistency with an entropy criterion is a less critical matter.

|1 -θ| ≤ 1 - 1 K 2 + O M ∆x |ũ R -ũL | . Ω k Ω j n jk N S

Extension to several space dimensions with unstructured grids

Without loss of generality, we suppose that Ω ⊂ R 2 is a polygonal domain that is covered by a set of N polygonal cells (Ω j ) 1≤j≤N . Let Γ be a face of a cell Ω j , 1 ≤ j ≤ N . If Γ ⊂ ∂Ω, we suppose that there exists a single k > N that will help to index ghost values for boundary conditions and we shall note Γ = Γ jk . If Γ ∩ ∂Ω = ∅, we suppose that the mesh is admissible in the sense that there exists a single 1 ≤ k ≤ N such that Γ = Ω j ∩ Ω k . Moreover, for 1 ≤ j ≤ N and 1 ≤ k ≤ N we suppose that Ω i ∩ Ω j can either be empty, a vertex or a single face of the mesh. If Γ jk be the face of a cell Ω j then n jk will denote the unit vector normal to Γ jk pointing out of Ω j . We define N (j) the set of indices k such that Γ jk is a face of Ω j . Let E = {(j, k) | 1 ≤ j, k ≤ N, k ∈ N (j)} and

E ext = {(j, k) | 1 ≤ j ≤ N, k ∈ N (j), Γ jk ⊂ ∂Ω}.
In sequel x = (x 1 , x 2 ) ∈ R 2 will denote the space variable.

We will now present a natural extension of our discretization strategy for the case of multi-dimensional problems with unstructured grids. Within this framework, the classical Lagrange-Remap algorithm involves tracking a genuine multi-dimensional moving mesh. This task is a very delicate matter as the mesh may be dramatically distorted during the simulation. We will here present a much simpler approach that relies on the alternative guideline proposed in section 3. A similar approach was used to derive an explicit scheme for two-component interface problems in [?].

Consider the operator splitting of (1) into the following systems

     ∂ t ρ + ρdiv(u) = 0, ∂ t (ρu) + ρudiv(u) + ∇p = 0, ∂ t (ρE) + ρEdiv(u) + div(P u) = 0, ( 52a
) (52b) (52c) and      ∂ t ρ + (u • ∇)ρ = 0, ∂ t (ρu) + (u • ∇)ρu = 0, ∂ t (ρE) + (u • ∇)ρE = 0. ( 53a
) (53b) (53c)
Before going any further, let us note that we obtain similar properties as for the systems ( 4) and (5). Indeed system (52) is a quasilinear hyperbolic system that involves the two nonlinear acoustic waves of velocity ±c and two null velocity contact discontinuities waves. System (52) only involves acoustic phenomena while freezing the material transport, while (53) is pure multi-dimensional transport system at the material velocity u.

We adopt the same strategy as in section 3: given a fluid state (ρ, ρu, ρE) n j ,

• update the fluid state to the value (ρ, ρu, ρE) n+1- j by approximating the solution of (52),

• update the fluid state to the value (ρ, ρu, ρE) n+1 j by approximating the solution of (53).

Approximation of the acoustic system (52) System (52) can be expressed

∂ t τ -τ (x, t)div(u) = 0, ∂ t u + τ (x, t)∇p = 0, E t + τ (x, t)div(pu) = 0.
Using the same lines as in section 3 we consider a Suliciu-type relaxation approximation

∂ t τ -τ (x, t)div(u) = 0, ∂ t u + τ (x, t)∇Π = 0, E t + τ (x, t)div(pu) = 0, Π t + τ (x, t)a 2 div(u) = ν(p -Π),
in the regime ν → +∞. Once again, for t ∈ [t n , t n + ∆t), this task is achieved by setting Π(x, t n ) = p(x, t n ) and then solving the relaxation system for ν = 0. We approximate again τ (x, t)∂ xr by τ (x, t n )∂ xr for r = 1, 2 when t ∈ [t n , t n + ∆t). In the regime λ = 0 our approximation of (52) becomes

∂ t τ -τ (x, t n )div(u) = 0, ∂ t u + τ (x, t n )∇Π = 0, E t + τ (x, t n )div(pu) = 0, Π t + τ (x, t n )a 2 div(u) = 0. ( 54
)
If b is a flow parameter and b n j is an approximation of 1 |Ωj | Ωj b(x, t n ) dx, we solve (54) thanks to the following classical Finite-Volume method

                             u n+1- j = u n j -τ n j ∆t k∈N (j) σ jk Π * ,θ jk n jk , Π n+1- j = Π n j -τ n j ∆t k∈N (j) σ jk (a jk ) 2 u * jk , τ n+1- j = τ n j + τ n j ∆t k∈N (j) σ jk u * jk , E n+1- j = E n j -τ n j ∆t k∈N (j) σ jk Π * ,θ jk u * jk , (55a) (55b) (55c) (55d)
where

σ jk = |Γ jk |/|Ω j |.
The three scalar quantities a jk , Π * ,θ jk and u * jk that respectively represent an average sound velocity, a pressure and the normal velocity at the face Γ jk . In order to define these quantities, we classicaly take advantage of the fact that (54) is rotational invariant. This allows to associate in the referential of each face Γ jk a Suliciu-type relaxation approximation of a one-dimensional Riemann problem in the frame of the face. Noting ∈ {n, n+1-}, this leads us to set

a jk ≥ max[(ρc) n j , (ρc) n k ], u * jk = 1 2 n T jk (u j + u k ) - 1 2a jk (Π k -Π j ), Π * ,θ jk = 1 2 (Π j + Π k ) - a jk θ jk 2 n T jk (u k -u j ). ( 56a
) (56b) (56c)
When = n the solver is explicit and when = n + 1-, the solver is implicit.

Approximation of the transport system (53)

In order to approximate the solution of (53), we simply use an upwind Finite-Volume scheme. Let ϕ ∈ {ρ, ρu 1 , ρu 2 , ρE}, we set

ϕ n+1 j = ϕ n+1- j -∆t k∈N (j) σ jk u * jk ϕ n+1- jk + ∆tϕ n+1- j k∈N (j) σ jk u * jk , (57)
where ϕ n+1- jk is defined by the upwind choice with respect to the sign of u * jk , namely

ϕ n+1- jk = ϕ n+1- j , if u * jk > 0, ϕ n+1- k , if u * jk ≤ 0.
Proposition 6. The overall numerical scheme composed by the discretization steps (55a)-(55d) and (57) is conservative with respect to the variable ρ, ρu and ρE, for both the implicit solver and the explicit solver. The update of these variables from t n to t n+1 reads

ρ n+1 j -ρ n j + ∆t k∈N (j) σ jk ρ n+1- jk u * jk = 0, (ρu) n+1 j -(ρu) n j + ∆t k∈N (j) σ jk (ρu) n+1- jk u * jk + Π * ,θ jk n jk = 0, (ρE) n+1 j -(ρE) n j + ∆t k∈N (j) σ jk (ρE) n+1- jk + Π * ,θ jk u * jk = 0. (58a) (58b) (58c)
The semi-implicit solver obtained for = n + 1-can be decomposed along the following steps: the acoustic step first involves solving the linear system (55a)-(55b) for computing the acoustic velocity u n+1- j and pressure term Π n+1- j . The acoustic step is completed by the update of τ n+1- j and E n+1- j thanks to the explicit procedures (55c) and (55d). The last stage of the semi-implicit solver is achieved thanks to the explicit transport scheme (57).

We want now to investigate further the implicit system involved with the semi-implicit method for the specific case of wall-boundary conditions that we implement by imposing ghost values Π n+1- k and n T jk u n+1- k for a boundary face Γ jk ⊂ ∂Ω, where 1

≤ j ≤ N and k ∈ N (j), k > N with Π n+1- k = Π n+1- j , n T jk u n+1- k = -n T jk u n+1- j . (59)
We have the following proposition.

Proposition 7. We consider the case of the semi-implicit solver with implementation of wall boundary conditions (59) and a uniform choice of a, i.e. a jk = a for all 1 ≤ j ≤ N and k ∈ N (j). If τ n j > 0 for all 1 ≤ j ≤ N , then the linear system (55a)-(55b) always possesses a single solution for any ∆t > 0 and θ jk > 0.

Proof. For the sake of readability, we shall note here u n+1- j = u j and Π n+1- j = Π j . The finite-dimension linear system (55a)-(55b) reads

             |Ω j |u j + τ n j ∆t k∈N (j) |Γ jk | 1 2 (Π j + Π k ) - aθ jk 2 n T jk (u k -u j ) n jk = |Ω j |u n j , |Ω j |Π j + τ n j ∆t k∈N (j) |Γ jk |a 2 1 2 n T jk (u j + u k ) - 1 2a (Π k -Π j ) = |Ω j |Π n j . (60a) (60b)
This system admits a unique solution if and only if u j = 0, Π j = 0, 1 ≤ j ≤ N is the only solution of the particular case obtained for u n j = 0, Π n j = 0, 1 ≤ j ≤ N . Thus, let us now suppose that the right members of (60) are null, we proceed using an energy estimate type proof. Left multiply (60a) by 2u T j τ n j ∆t and sum over j,

we obtain 0 = N j=1 2|Ω j ||u j | 2 τ n j ∆t + N j=1 k∈N (j) |Γ jk |(Π j + Π k )(u T j n jk ) - N j=1 k∈N (j) |Γ jk |aθ jk (u T j n jk )(u k -u j ) T n jk . (61)
Accounting for the fact that k∈N (j)

|Γ jk |n jk = 0, the second term of (61) verifies

N j=1 k∈N (j) |Γ jk |(Π j + Π k )(u T j n jk ) = N j=1 k∈N (j) |Γ jk |Π k u T j n jk .
Using boundary conditions (59), the third term of (61) reads

N j=1 k∈N (j) |Γ jk |aθ jk (u T j n jk )(u k -u j ) T n jk = (j,k)∈E |Γ jk |aθ jk (u T j n jk )(u k -u j ) T n jk + (u T k n kj )(u j -u k ) T n kj + (j,k)∈E ext |Γ jk |aθ jk (u T j n jk )(u k -u j ) T n jk = - (j,k)∈E |Γ jk |aθ jk (u k -u j ) T n jk 2 -2 (j,k)∈E ext |Γ jk |aθ jk [(u T j n jk )] 2 .
Finally we see that (61) is equivalent to

0 = N j=1 2|Ω j ||u j | 2 τ n j ∆t + N j=1 k∈N (j) |Γ jk |Π k u T j n jk + (j,k)∈E |Γ jk |aθ jk (u k -u j ) T n jk 2 + 2 (j,k)∈E ext |Γ jk |aθ jk [(u T j n jk )] 2 . ( 62
)
Let us turn to the pressure equation (60b), we multiply by 2Π j τ n j a 2 ∆t and sum over all 1 ≤ j ≤ N , this yields

0 = N j=1 2|Ω j |Π 2 j τ n j a 2 ∆t + N j=1 k∈N (j) |Γ jk |n T jk (u j + u k )Π j - N j=1 k∈N (j) 1 a |Γ jk |(Π k -Π j )Π j . (63)
Using once again

k∈N (j)
|Γ jk |n jk = 0, we have for the second term of (63) that

N j=1 k∈N (j) |Γ jk |n T jk (u j + u k )Π j = N j=1 k∈N (j) |Γ jk |n T jk u k Π j .
Accounting for (59), the third term of (63) verifies

N j=1 k∈N (j) |Γ jk | a (Π k -Π j )Π j = 1 a (j,k)∈E |Γ jk | (Π k -Π j )Π j -(Π j -Π k )Π k + 1 a (j,k)∈E ext |Γ jk |(Π k -Π j )Π j = - 1 a (j,k)∈E |Γ jk |(Π k -Π j ) 2 .
Then, we see that (63) also reads

0 = N j=1 2|Ω j |Π 2 j τ n j a 2 ∆t + N j=1 k∈N (j) |Γ jk |n T jk u k Π j + 1 a (j,k)∈E |Γ jk |(Π k -Π j ) 2 . ( 64
)
We now remark that

N j=1 k∈N (j) |Γ jk |n T jk u k Π j = (j,k)∈E |Γ jk |(n T jk u k Π j + n T kj u j Π k ) + (j,k)∈E ext |Γ jk |n T jk u k Π j = (j,k)∈E |Γ jk |n T jk (u k Π j -u j Π k ) - (j,k)∈E ext |Γ jk |n T jk u j Π j ,
and also that

N j=1 k∈N (j) |Γ jk |Π k u T j n jk = (j,k)∈E |Γ jk |(Π k u T j n jk + Π j u T k n kj ) + (j,k)∈E ext |Γ jk |Π k u T j n jk = (j,k)∈E |Γ jk |n T jk (Π k u j -Π j u k ) + (j,k)∈E ext |Γ jk |Π j n T jk u j . Therefore N j=1 k∈N (j) |Γ jk |Π j u T k n jk + N j=1 k∈N (j) |Γ jk |Π k u T j n jk = 0.
Thus, summing (62) and (64), we obtain

0 = N j=1 2|Ω j | τ n j ∆t |u j | 2 + Π 2 j a 2 + (j,k)∈E |Γ jk | aθ jk (u k -u j ) T n jk 2 + (Π k -Π j ) 2 a +2 (j,k)∈E ext |Γ jk | aθ jk [(u T j n jk )] 2 .
This implies that |u j | = Π j = 0 for all 1 ≤ j ≤ N .

Remark 4. It is possible to derive a similar proof for the case of periodic boundary conditions.

We now examine the stability of the multi-dimensional operator splitting strategy (55), ( 56) and (57). The acoustic step (55) in the explicit cases = n is stable under the CFL condition ∆t max

1≤j≤N τ n j max k∈N (j) σ jk a jk ≤ 1 2 . ( 65
)
For both the explicit scheme = n and semi-implicit scheme = n + 1-, the transport step (57) is stable under the CFL condition

∆t max 1≤j≤N   k∈N (j) σ jk (n T jk u * ,θ jk )   ≤ 1. (66)
When one uses the semi-implicit scheme = n + 1-, the condition (66) becomes implicit as the computation of u * ,θ jk depends on a given ∆t. In our simulations with the semi-implicit scheme, we chose to compute ∆t thanks to the CFL condition (66) with the value u * ,θ jk given by the fully explicit scheme = n. It is then possible to check a posteriori that this ∆t value matches (66).

We gather thereafter the properties of the explicit and semi-implicit multi-dimensional schemes.

Theorem 2. Suppose that (65), (66) and (12) are satisfied. The explicit scheme defined by (55) and (57) with = n verifies 1. the scheme is conservative with respect to the density ρ, the momentum ρu and total energy ρE, 2. the density ρ n j is positive for all j and n > 0 provided that ρ 0 j is positive for all j, 3. if θ = O(M ), then the truncation error of the numerical scheme is uniform with respect to M < 1, 4. if (43) is verified then the numerical scheme is equipped with a discrete entropy inequality, 5. if (43) is verified then e n j > 0 for all j ∈ Z and all n ∈ N.

Theorem 3. Suppose that (66) and (12) are satisfied. The semi-implicit scheme defined by (55) and (57) with = n + 1-verifies 1. the scheme is conservative with respect to the density ρ, the momentum ρu and total energy ρE, 2. the density ρ n j is positive for all j and n > 0 provided that ρ 0 j is positive for all j, 3. if θ = O(M ), then the truncation error of the numerical scheme is uniform with respect to M < 1.

Let us note that the implicit treatment of the acoustic step leads to a CFL restriction (66) based only on (slow) material waves.

Backward facing step

We consider now the case of an inviscid flow passing a backward facing step as derived from [START_REF] Cordier | An Asymptotic-Preserving all-speed scheme for the Euler and NavierStokes equations[END_REF]. The computational domain is Ω = [0, 18] × [0, 2] (0, 4) × (0, 1). The initial condition is given by

ρ 0 (x 1 , x 2 ) = 10, u 0 (x 1 , x 2 ) = 1, p 0 (x 1 , x 2 ) = 10 5 , v 0 (x 1 , x 2 ) = 0.
We

Compressible flow examples

In this section, we assess the ability of our operator splitting scheme to handle cases where the flow may not remain uniformly in the same Mach regime over the whole computational domain Ω. We will see that even with a centred pressure discretization (which corresponds to the choice θ = 0), the solution remains stable but may be less precise in area where the Mach number is of order 1. The semi-implicit scheme becomes slower than the explicit scheme when the Mach number is of order 1 as the benefit from using a material CFL (66) condition instead of an acoustic CFL (65) becomes less beneficial but requires solving a linear system.

1D Sod shock tube

We consider a variant of the classical Sod shock tube [START_REF] Sod | Numerical methods in fluid dynamics. Initial and initial-boundary value problems[END_REF], that consists in solving the one-dimensional Riemann problem over Ω = [0, 1] defined by the initial conditions (ρ, u, P ) = (1.0, 0.0, 10 5 ) for x < 0.5 and (ρ, u, P ) = (0.1, 0.0, 10 4 ) for x > 0.5. We impose Neumann boundary conditions during the test. The domain is discretized over a 1000-cell grid. This resulting Mach number verifies 0 < M < 0.95, so that we have both low Mach and order 1 Mach values. We plot the solution at t = 3.1 × 10 -4 s.

Figure 7 displays the results obtained with EX(θ) and IMEX(θ) for θ = 1 and θ = 0. We use as reference solution an approximation computed with EX(θ = 1) using a 10 000-cell mesh. All schemes show a good agreement with the reference solution. The schemes EX(θ = 0) and IMEX(θ = 0) schemes are slightly less diffused than the EX(θ = 1) and IMEX(θ = 1) schemes results. Let us underline that despite part of the solutions clearly do not belong to the low Mach regime since M 0.95, the schemes EX(θ = 0) and IMEX(θ = 0) are stable and provide good numerical results while involving a centered pressure discretization with θ ij = 0.

2D-Riemann problem

We consider a 2D Riemann problem that consists of 4 shock waves [START_REF] Lax | Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[END_REF]. We consider the domain Ω

= [0, 1] 2 . The initial condition is (ρ, u 1 , u 2 , P )(x 1 , x 2 , t = 0) =          (0.1380, 1.206, 1.206, 0.029),
for x 1 < 0.5, x 2 < 0.5 (0.5323, 0.000, 1.206, 0.300), for x 1 > 0.5, x 2 < 0.5 (0.5323, 1.206, 0.000, 0.300), for x 1 < 0.5, x 2 > 0.5 (1.5000, 0.000, 0.000, 1.500), for x 1 > 0.5, x 2 > 0.5

We impose Neumann boundary conditions. This configuration leads to a Mach number that ranges from 10 -5 to 3.15, i.e. according to the regions of the computation domain, the flow belongs to the low Mach regime or the order 1 Mach regime. We consider as a reference solution the approximation obtained with EX(θ = 1) for a 200 × 200-cell Cartesian mesh. Figures 8, 9 10 and 11 display the result at t = 0.4 s. We observe in figure 8 and figure 9 that EX(θ = 0) and IMEX(θ = 0) schemes are stable for this test case with both low Mach and order 1 Mach number values regions. Both figure show that the wave pattern at the center of the domain shape is better captured with coarse meshes when one uses the corrected schemes (θ = 0). A 1D cut along the axis y = x as depicted in figure 10, also corroborates this observation: the approximation obtained with EX(θ = 0) and IMEX(θ = 0) schemes are closer to the 200 × 200-cell reference solution thanks to the numerical diffusion reduction. Nonetheless, we observe on a 1D cut along the x = 0.75 axis in figure 11 a spurious overshot for both density and pressure located at the shock front with EX(θ = 0) and IMEX(θ = 0). This suggests that a small value of θ allows to improve the precision of the scheme by reducing the numerical diffusion but it may cause overshoots if the value of θ becomes too small relatively to the local behavior of flow. In all our numerical experiments the scheme seems to remain stable for any value of θ ∈ (0, 1). Let us note that even if the pressure gradient is given a centred treatment (θ = 0), the transport step introduce some numerical diffusion (independent of M ) that stabilize the scheme see (32).

In table 3 we observe that the choice of θ does not impact the number of time steps and CPU time. For this case, while the number of time steps is slightly reduced by about 30%, the semi-implicit schemes are much slower due to the time required for solving the linear system involved with the schemes.

As a partial conclusion of this section, we can observe that for tests that strongly involve the compressibility of the fluid both semi-implicit and explicit schemes seem to be very robust, independantly of the choice of θ within [0, 1]. However, if the low Mach correction is too important, i.e. the value of θ is too close to 0 we witnessed a deterioration of the numerical approximation with the appearance of overshoots in the vicinity of shock fronts. Then some numerical criterion may be constructed with good properties, θ ij = min

|u * ij | max(c n i ,c n j )
, 1 for instance. We also observed that the benefit in terms of CPU time of the semi-implicit scheme vanishes when the Mach number becomes of order 1.

The implementation of the criterion on θ to recover a discrete entropy inequality does not allow to recover a good low Mach behaviour as it was expected from its low Mach analysis. Finding from a theoretical point of view a criterion on θ that allows to recover a good low Mach behaviour and avoid spurious phenomenon that may occur if θ is too small for a given configuration is still an open problem. Remark 6. The robustness of the scheme with respect to the modification θ ≥ 0 seems to be linked to the Lagrange-Projection decomposition approach. Indeed, numerical evidences not presented here show that a modified relaxation scheme written in Eulerian coordinates is unstable outside of the low-Mach regime for value of θ that are too small. For instance, The Sod shock test tube require θ > 0.63 to be stable.

Conclusion

We proposed a conservative operator splitting based Lagrange-Projection like numerical strategy for approximating the gas dynamics that decouples acoustic and transport phenomenons. The operator splitting scheme is positive for the density, the internal energy and entropic under classical CFL conditions. For one-dimensional problem, this procedure is equivalent to a Lagrange-Projection discretization. We presented an analysis of the way the truncation error depends on the Mach number for one-dimensional problems. In the low Mach regime, the truncation error of the scheme showed to be non-uniform with respect to the Mach number M . This allowed us to modify the operator splitting scheme in order to recover a uniform truncation error in term of M by altering the numerical flux in the acoustic approximation. We showed that this modification can be obtained thanks to a simple approximate Riemann solver that is consistent with the integral form of the PDEs. This modified operator splitting scheme is conservative and endowed with good stability properties with respect to the positivity of the density, the internal energy under classical acoustic CFL conditions that depend on M . The resulting scheme allows to deal with tests where the flow regime may vary from low to high Mach values.

We showed that this splitting strategy has a natural extension to multi-dimensional problems discretized over unstructured meshes. A simple and efficient semi-implicit scheme that is stable under CFL conditions based on the material velocity is also proposed and leads to an all-regime numerical scheme, following the ideas paved by [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF] for one-dimensional problems.

Future developments include extensions to high-order methods and approximation of other systems for the simulation of multi-material flows.

From [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF], if ∆x = 1 2 (∆x L + ∆x R ) with ∆x L > 0, ∆x R > 0 and ∆t > 0 are respectively space and time steps that verify the CFL condition

max 1≤k≤l |λ k | ∆t min(∆x L , ∆x R ) ≤ 1 2 , (76)
such an approximate Riemann solver is said to be consistent with the integral form of (73) over the interval

[ -∆x L 2 , ∆x R 2 ] × [0, ∆t] if [ -∆x L 2 , ∆x R
2 ]×[0,∆t] [∂ t U RP + ∂ x G(U RP )] dxdt = 0, in other words if

G(U R ) -G(U L ) = l k=1 λ k (U k+1 -U k ). ( 77
)
Regarding the consistency with the entropy inequality (74), the simple approximate Riemann solver is said to be consistent with the integral form of (74) if and only if under the CFL condition (76) we have

q(U R ) -q(U L ) ≤ l k=1 λ k η(U k+1 ) -η(U k ) . (78)
Hereafter and using classic notations, (∆x j ) j∈Z and ∆t represent the space steps and constant time step of the mesh under consideration to define the approximate solutions. More precisely, we define the mesh interfaces x j+1/2 = x j-1/2 + ∆x j for j ∈ Z, the intermediate times t n = n∆t for n ∈ N, and we note U n j the approximate value of U at time t n and on the cell [x j-1/2 , x j+1/2 ). For n = 0 and j ∈ Z, we set U 0 j = 1 ∆x x j+1/2

x j-1/2 U 0 (x) dx where U 0 (x) is the initial condition. Then, the explicit in time Godunov-type scheme reads

     U n+1 j = U n j - ∆t ∆x j (G n j+ 1 2 -G n j-1 2
),

G n j+ 1 2 = G(U n j , U n j+1 ), (79a) (79b)
with

G(U L , U R ) = 1 2 G(U L ) + G(U R ) - l k=1 |λ k |(U k+1 -U k ) . (80)
Moreover, if the simple approximate Riemann solver is consistent with the entropy inequality (74), then the numerical scheme defined by (79) satisfies the following discrete entropy inequality

     η(U n+1 j ) ≤ η(U n j ) - ∆t ∆x j (q n j+ 1 2 -q n j-1 2
),

q n j+ 1 2 = q(U n j , U n j+1 ), with q(U L , U R ) = 1 2 q(U L ) + q(U R ) - l k=1 |λ k | S(U k+1 ) -S(U k ) . (82)
The CFL condition associated with this (explicit in time) Godunov-type scheme reads

max 1≤k≤l | λ k (U n j , U n j+1 ) | ∆t min(∆x j , ∆x j+1 ) ≤ 1 2 ,
for all j. Again, we refer to [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF] for more details. To conclude this paragraph, let us observe that the numerical flux G(U L , U R ) and the entropy numerical flux q(U L , U R ) are clearly consistent in the classical sense, namely G(U, U) = G(U) and q(U, U) = q(U) provided that the intermediate states of the approximate Riemann solver are such that U k = U for all k = 1, ..., l as soon as U L = U R = U.

Lemma 2 .

 2 Let θ ∈ R, and e * ,θ k

Figure 1 :

 1 Figure1: the face Γ jk = Ω j ∩ Ω k defined the segment (N S) has a unit normal vector n jk oriented from Ω j to Ω k .

Table 1 :Figure 2 :

 12 Figure 2: Vortex in a box test case with a Cartesian mesh. Profile at time t = 0.125 s of the velocity magnitude for (a) EX(θ = 1), (b) EX(θ = O(M )) with a 50 × 50-cell Cartesian mesh, (c) velocity magnitude obtained with EX(θ = 1) using a 400 × 400 Cartesian mesh and (d) Mach number obtained with EX(θ = 1) using a 400 × 400 Cartesian mesh.

Figure 3 :

 3 Figure 3: Vortex in a box test case with a Cartesian mesh. Profile at time t = 0.125 s of the velocity magnitude for the IMEX(θ = 1) scheme (left) and the IMEX(θ = O(M )) scheme (right) on a 50 × 50-cell Cartesian mesh.

  impose an inlet boundary condition at {0} × [1, 2] and an outlet boundary condition at {12} × [0, 2]. Wall boundary conditions are set on other boundaries. This configuration leads to a low Mach flow with the order of magnitude 10 -3 ≤ M ≤ 10 -2 . All tests are performed with a 220 × 20 Cartesian space grid.Figure4and 6 display the flow profile at t = 50 s, we observe that EX(θ = 1) and IMEX(θ = 1) schemes do not capture the vortex of the fluid in the low Mach velocity region. On the contrary, thanks to the low Mach correction EX(θ = O(M )) and IMEX(θ = O(M )) schemes are both able to capture this vortex with a coarse Cartesian Mesh. In term of CPU cost, measure are presented in table 2. We observe that the IMEX(θ = O(M )) scheme is 139.29 times faster than the EX(θ = O(M )) scheme thanks to the implicit treatment of the acoustic step, thanks to the use of material velocity CFL condition (66).

Figure 4 :

 4 Figure 4: Backward facing step test case. Profile at t = 50 s of the velocity magnitude and stream lines for the EX(θ = 1) scheme (top) and the EX(θ = O(M )) scheme (bottom) on a 220 × 20-cell Cartesian mesh.

Figure 5 :

 5 Figure 5: Backward facing step test case. Mapping at t = 50 s of the Mach number values for EX(θ = O(M )) scheme with a 220 × 20-cell Cartesian mesh.

Figure 6 :

 6 Figure 6: Backward facing step test case. Mapping at t = 50 s of the velocity magnitude and stream lines obtained with the IMEX(θ = 1) scheme (top) and the IMEX(θ = O(M )) scheme (bottom) using a 220 × 20-cell Cartesian mesh.

Figure 7 :

 7 Figure 7: 1D Sod shock tube test case. Profile at t = 3.1×10 -4 s of the density (top left), velocity magnitude (top right), pressure (bottom left) and Mach number (bottom right) for the EX(θ = 1), EX(θ = 0), IMEX(θ = 1), IMEX(θ = 0) using a 1000-cell grid, together with reference solution.

Figure 8 :

 8 Figure 8: 2D Riemann problem with a Cartesian mesh. Profile at t = 0.4 s of the velocity magnitude for (a) EX(θ = 1), (b) EX(θ = 0) with a 50 × 50-cell mesh, (c) velocity magnitude and (d) Mach number with EX(θ = 1) using a 200 × 200 mesh.

Figure 9 :

 9 Figure 9: 2D Riemann problem test case. Profile at t = 0.4 s of the velocity magnitude for IMEX(θ = 1) (left) and IMEX(θ = 0) (right) on a 50 × 50-cell Cartesian mesh.

Figure 10 :

 10 Figure 10: 2D Riemann problem test case. Cut profile along y = x at t = 0.4 s of the density, velocity magnitude, pressure and Mach number for EX(θ = 1), EX(θ = 0), IMEX(θ = 1) and IMEX(θ = 0) using a 50 × 50 mesh together with the 200 × 200-cell reference solution.

Figure 11 :

 11 Figure 11: 2D Riemann problem test case. Cut profile along x = 0.75 at t = 0.4 s of the density, velocity magnitude, pressure and Mach number for EX(θ = 1), EX(θ = 0), IMEX(θ = 1) and IMEX(θ = 0) using a 50 × 50 mesh together with the 200 × 200-cell reference solution.

Table 2 :

 2 Backward facing step test case. Comparison of the number of iterations and CPU time of EX(θ = 1), EX(θ = O(M )), IMEX(θ = 1) and IMEX(θ = O(M )) schemes to obtain solutions of figure 4 and figure 6.

	Numerical scheme	EX(θ = 1)	EX(θ = O(M ))	IMEX(θ = 1)	IMEX(θ = O(M ))
	Number of time steps	5 258 803	5 258 784	4 384	4 970
	CPU time (s)	62 805.6 (17h 27min) 69 764.7 (19h 22min) 418.37 (6min 58s) 500.83 (8min 20s)

Acknowledgement

The authors would like to thank Pierre-Arnaud Raviart for sharing his personal notes and for the many rich discussions that brought an invaluable contribution to this work.

Numerical results

In this section, we present numerical results computed thanks to the general operator splitting strategy (55), ( 56) and (57) with the following schemes:

• EX(θ = 1): the explicit operator splitting scheme obtained for θ jk = 1 and = n,

• EX(θ = O(M )) : the explicit modified operator splitting scheme obtained with the low Mach correction θ jk = min |u * jk |/max(c n j , c n k ), 1 and = n,

• EX(θ = 0) : the explicit modified operator splitting scheme with centered pressure gradient θ jk = 0 and = n,

• IMEX(θ = 1) : the semi-implicit operator splitting scheme with θ jk = 1 and = n + 1-,

• IMEX(θ = O(M )) : the modified semi-implicit operator splitting scheme with = n + 1-and a low Mach correction θ jk defined as in the case of EX(θ = O(M )),

• IMEX(θ = 0) : the modified semi-implicit operator splitting scheme with a centered pressure gradient θ ij = 0 and = n + 1-.

Remark 5. The choice of the modification θ jk = min

, 1 corresponds to a low Mach correction. Indeed, this choice is non-dimensional, in (0, 1), such that θ = O(M ) in the low Mach regime and θ = 1 for large Mach numbers. In this latter case, we then recover the classical scheme without modification.

In the sequel, we shall consider that the fluid follows a perfect gas equation of state p = (γ -1)ρe with a specific heat ratio γ = 1.4. We will test schemes on both low Mach and order 1 Mach number test cases.

Low Mach number examples

In this section we will consider low Mach tests and try to examine two questions: the accuracy gain for simulations on coarse grid in the low Mach regime thanks to the proposed correction, then the benefit of using a semi-implicit strategy in term of CPU time.

Vortex in a Box

We consider a test performed in [START_REF] Ph | A projection method for low speed flows[END_REF]. The computational domain is Ω = [0, 1] 2 with an initial condition is given by

No-slip boundary conditions are imposed on the domain boundaries. The Mach number for the resulting flows is of order 0.026, so that we are in the low Mach regime. Results are displayed in table 1 and figures 2 and 3. We first use the schemes EX(θ = 1) with a 400 × 400-cell and a 50 × 50-cell mesh. As expected the scheme performs poorly on the coarse mesh and the gain of accuracy is obvious when one refines the mesh: a mesh size of order M is required, but it comes at a much higher price in terms of CPU time as we can see on table 1. The EX(θ = O(M )) scheme gives good results even with the coarse 50 × 50-cell grid. With the low Mach correction scheme the connection between the accuracy of the solution and the mesh size does not seem to be constrained by M . Therefore, for a given target accuracy on a relatively coarse mesh, this numerical scheme is also much cheaper in term of CPU time.

Let us now turn to the semi-implicit strategies where the time step was chosen in agreement with the material CFL condition (66). While the IMEX(θ = 1) is not CPU intensive on a coarse mesh the results are very altered by the numerical diffusion. The IMEX(θ = O(M )) scheme performs fast and allows to recover numerical results that are as good as EX(θ = O(M )). As with the EX(θ = M ) scheme the accuracy seems much less constrained by the Mach number when it comes to choosing the time step and the mesh size. As we can see in table 1, the IMEX(θ = O(M )) scheme is 3.34 times faster than the EX(θ = O(M )).

A Classical Lagrange-Projection for one-dimensional gas dynamics

In this section we briefly recall the classical Lagrange-Projection (or Lagrange-Remap) procedure for deriving a Finite Volume discretization within a one-dimensional framework. For a detailed description we refer the reader to [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Després | Lois de conservations Eulériennes, Lagrangiennes et méthodes numériques[END_REF]. Let (X, t)R × [t n , t n + ∆t) → χ be the mapping defined by

The pair (X, t) is usually referred to as the Lagrangian system of coordinates: a particle of fluid at the position X at instant t = t n will be located at x = χ(X, t), t ∈ [t n , t n + ∆t]. If (x, t) → b is a mapping that provides an Eulerian representation of a parameter b, one defines a Lagrangian representation of b as the function (X, t) → b Lag by setting b Lag (X, t) = b(χ(X, t), t). The system (3) is equivalent to

It is common to introduce a mass coordinate m defined by dm = ρ(X, t n )dX in order to obtain the equivalent conservation laws (with a slight abuse of notation)

Straightforward calculations show that (68) (which is nothing but ( 6)) is hyperbolic over the phase space is linearly degenerate. It is important to note that the material transport phenomenons are frozen in system (68) which explains why the characteristics speeds of the system only involve the sound velocity c.

Before going any further, we introduce classical notations: let ∆t > 0 and ∆x > 0 be respectively the time and space steps. We define the Eulerian mesh interfaces x j+1/2 = j∆x for j ∈ Z, and the intermediate times ) the approximate value b respectively within the j th Eulerian cell [x j-1/2 , x j+1/2 ) at instant t = t n (resp. t = t n+1 ). We need to introduce a moving Lagrangian mesh (with respect to the Eulerian mesh) whose cell j at instant t n is [x j-1/2 , x j+1/2 ) and at instant t = t n+1 is [x * j-1/2 , x * j+1/2 ). The value of the parameter b at instant t n (resp. t = t n+1 ) in the Lagrangian cell j is noted b Lag j (resp. b n+1- j ). Given a fluid state (ρ, ρu, ρv, ρE) n j , j ∈ Z at instant t n , the Lagrange-Projection strategy proposes the following update procedure.

1. Build the discrete Lagrangian fluid state at instant t n by setting (V Lag ) j = (τ n j , u n j , v n j , E n j );

Update the Lagrangian fluid state into the value (V

) by approximating the solution of (68); 3. Build the updated value (ρ, ρu, ρv, ρE) n+1 j by remapping the Lagrangian state (V Lag ) n+1- j onto the Eulerian mesh.

The Lagrangian step (t n → t n+1-) We propose to approximate the solution of (68) using the acoustic scheme [START_REF] Després | Inégalité entropique pour un solveur conservatif du système de la dynamique des gaz en coordonnées de lagrange[END_REF][START_REF] Després | Lois de conservations Eulériennes, Lagrangiennes et méthodes numériques[END_REF]. This leads to

where the interfaces terms are defined by

The acoustic scheme (69) with (70) provides the same update of the flow variable as the scheme ( 9) with [START_REF] Dellacherie | Construction of modified Godunov type schemes accurate at any Mach number for the compressible Euler system[END_REF]. Let us mention that a direct proof of stability for the acoustic scheme is available in [START_REF] Després | Inégalité entropique pour un solveur conservatif du système de la dynamique des gaz en coordonnées de lagrange[END_REF] under the CFL criterion [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF].

The projection (or remapping) step (t n+1-→ t n+1 ) The aim of this step is to project the solution obtained at the end of the Lagrangian step onto the Eulerian cells

), a standard way to achieve to goal consists in: first, approximating the position of the Lagrangian mesh interfaces at instant t n+1 by setting x * j+1/2 = x j+1/2 + u * j+1/2 ∆t; second reaveraging the conservative variable unknowns over the Eulerian mesh by setting [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]

Noting ∆x * j = x * j+1/2 -x * j-1/2 and ε(j, n) = -sign(u * j+1/2 )1/2 one obtains the update formula

The update formula (72) matches the classic upwind scheme. Consequently this is the same numerical scheme as [START_REF] Dauvergne | Renormalization of the numerical diffusion for an upwind finite volume method. application to the simulation of Kelvin-Helmholtz instability. Finite volumes for complex applications[END_REF].

B Approximate Riemann solvers: Harten Lax and van Leer formalism

We briefly recall the Harten, Lax and van Leer formalism associated with the numerical approximation of the solutions (x, t) ∈ R × [0, +∞) → U ∈ R m of the general hyperbolic system of conservation laws

by means of the so-called approximate Riemann solvers and Godunov-type methods, where G : R m → R m is a smooth function. System (73) is supplemented with the validity of an entropy inequality

where U → (η, q) is a strictly convex entropy-entropy flux pair (see [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF]).

Solving the Riemann problem amounts to find the solution of (73) with the following piecewise constant initial data

for any given U L and U R in the phase space. It is well-known that the exact Riemann solution U(x/t; U L , U R ) is self-similar, i.e. depends only on the ration x/t. In order to approximate this solution, we consider a (selfsimilar) simple approximate Riemann solver U RP ( x t ; U L , U R ) made of l + 1 intermediate states U k separated by discontinuities propagating with velocities

C Riemann problem for the relaxation approximation of the acoustic system

We consider the Suliciu relaxation approximation of the Lagrangian gas dynamics equations expressed using a mass coordinate. The system reads

where a is a constant that verifies the subcharacteristic condition a > max (ρc) in order to prevent instabilities (see for instance [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF] for a rigorous proof). It is easy to prove that the convective part of ( 83) is strictly hyperbolic with three eigenvalues given by -a, 0 and a which correspond to linearizations of the exact eigenvalues -ρc, 0 and ρc for system (67). Interestingly, the characteristic fields are linearly degenerate, which allows to solve analytically the Riemann problem associated with (83) with λ = 0. More precisely, the exact Riemann solution

The intermediate states are easily recovered from the following formulas

Then, setting U L = (τ, u, v, E) T , the classical scheme can be understood in the Harten, Lax and van Leer formalism by considering the following approximate Riemann solver W m t ; U L , U R obtained by simply extracting the first four components from W m t ; U L , U R , in which we take Π at equilibrium, namely Π L = p L , Π R = p R .

More precisely, we have

where the intermediate states are given by (85), together with Π L = p L and Π R = p R .