F Réjasse

O Rapaud

G Trolliard

A Maître

Experimental investigation and thermodynamic modeling of the C-Me-O systems (where Me = Hf, Zr)

Introduction

• Zirconium and hafnium carbide (ZrC x and HfC x ) are considered as Ultra High Temperature Ceramics (UHTC) which can be easily oxidized under flowing dioxygen. • Oxygen atoms can substitute carbon atoms in the crystalline lattice of carbides involving the formation of oxycarbide phases ZrC x O y and HfC x O y .

• The thermodynamic modeling of the C-Me-O ternary phase diagrams requires thermodynamic and diagrammatic data especially to determine the solubility of the oxygen into the carbide sublattice.

Reaction

••

  . David et al., J. Eur. Ceram. Soc. 33 (2013) 165-179. b -F. Réjasse et al., J. David, RSC Adv. 5 (2015) 45341-45350. c -A. Fernández Guillermet, J. Alloys Compd. 217 (1995) 69-89. d -P.Y. Chevalier et al., Calphad. 28 (2004) 15-40. e -H. Bittermann, P. Rogl, J. Phase Equilibria. 18 (1997) 344-356. f -D. Shin et al., Calphad. 30 (2006) 375-386. Structural approach by XRD and TEM  understanding of the reaction mechanisms involved for different advancements ξreaction of MeO (g) and CO (g) within carbon-rich areas Evolution of the MeC x O 1-x phase Boudouard equilibrium Le Bail refinements The increase of the carbon amount in the crystalline structure implies the increase of the MeC x O 1-x lattice parameter. Purely reconstructive reaction  Impossible to obtain nanosized carbides by this synthesis way Summary • The amount of initial reactants is adjusted to produce several oxycarbide stoichiometries. • Heat treatment in a graphite furnace under flowing argon at 1650°C, 1750°C and 1850°C during 8 hours (ξ = 1). Between 1650°C and 1750 °C : no evolution of the solubility of O • At the same temperature, the stability domain of HfC x O 1-x solid solution is narrower than ZrC x The carbothermal reaction is a purely reconstructive reaction which operates through solid-gas equilibriums. • This reaction allowed to determine the stability domains of MeC x O y solid solutions. On the Zr-C-O system, the extent of the domain increases with the temperature. At the same temperature, the domain of HfC x O y is narrower than ZrC x O y . • These experimental data were used to assess the two ternary phase diagrams and to describe the evolution of the solid solutions with the temperature


Giacomo Dimarco
email: giacomo.dimarco@unife.it

Rapha Ël

Loub Ère

Marie-H Él

Ène Vignal
email: mhvignal@math.univ-toulouse.fr

STUDY OF A NEW ASYMPTOTIC PRESERVING SCHEME FOR THE EULER SYSTEM IN THE LOW MACH NUMBER LIMIT

Keywords: Low Mach number limit, Asymptotic preserving schemes, Euler system, stability analysis

This article deals with the discretization of the compressible Euler system for all Mach numbers regimes. For highly subsonic flows, since acoustic waves are very fast compared to the velocity of the fluid, the gas can be considered as incompressible. From the numerical point of view, when the Mach number tends to zero, the classical Godunov type schemes present two main drawbacks: they lose consistency and they suffer of severe numerical constraints for stability to be guaranteed since the time step must follow the acoustic waves speed. In this work, we propose and analyze a new unconditionally stable an consistent scheme for all Mach number flows, from compressible to incompressible regimes, stability being only related to the flow speed. A stability analysis and several one and two dimensional simulations confirm that the proposed method possesses the sought characteristics.

1. Introduction. Almost all fluids can be said to be compressible. However, there are many situations in which the changes in density are so small to be considered negligible. We refer to these situations saying that the fluid is in an incompressible regime. From the mathematical point of view, the difference between compressible and incompressible situations is that, in the second case, the equation for the conservation of mass is replaced by the constraint that the divergence of the velocity should be zero. This is due to the fact that when the Mach number tends to zero, the pressure waves can be considered to travel at infinite speed. From the theoretical point of view, researchers try to fill the gap between those two different descriptions by determining in which sense compressible equations tend to incompressible ones [START_REF] Asano | On the incompressible limit of the compressible euler equation[END_REF][START_REF] Lions | Incompressible limit for a viscous compressible fluid[END_REF][START_REF] Klainerman | Majda Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF][START_REF] Klainerman | Compressible and incompressible fluids[END_REF][START_REF] Schochet | The compressible Euler equations in a bounded domain: existence of solutions and the incompressible limit[END_REF]. In this article we are interested in the numerical solution of the Euler system when used to describe fluid flows where the Mach number strongly varies. This causes the gas to pass from compressible to almost incompressible situations and consequently it causes most of the numerical methods build for solving compressible Euler equations to fail. In fact, when the Mach number tends to zero, it is well known that classical Godunov type schemes do not work anymore. Indeed, they lose consistency in the incompressible limit. This means that when close to the limit, the accuracy of theses schemes is not sufficient to describe the flow. Many efforts have been done in the recent past in order to correct this main drawback of Godunov schemes, for instance by using preconditioning methods [START_REF] Turkel | Preconditioned methods for solving the incompressible and low speed compressible equations[END_REF] or by splitting and correcting the pressure on the collocated meshes [START_REF] Ph | A projection method for low speed flows[END_REF], [START_REF] Dellacherie | Preliminary results for the study of the Godunov Scheme Applied to the Linear Wave Equation with Porosity at Low Mach Number[END_REF][START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF], [START_REF] Grenier | An accurate low-Mach scheme for a compressible two-fluid model applied to free-surface flows[END_REF], [START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit : II. Godunov type schemes[END_REF][START_REF] Guillard | On the behavior of upwind schemes in the low Mach limit[END_REF][START_REF] Paillère | Comparison of low mach number models for natural convection problems[END_REF], [START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[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] Klein | Asymptotic adaptive methods for multi-scale problems in fluid mechanics[END_REF], [START_REF] Munz | The extension of incompressible flow solvers to the weakly compressible regime[END_REF][START_REF] Munz | Linearized acoustic perturbation equations for low Mach number flow with variable density and temperature[END_REF][START_REF] Park | Multiple pressure variables methods for fluid flow at all Mach numbers[END_REF], or instead by using staggered grids like in the famous MAC scheme, see for instance [START_REF] Bonner | Compressible subsonic flow on a staggered grid[END_REF], [START_REF] Harlow | A numerical fluid dynamics calculation method for all flow speeds[END_REF], [START_REF] Van Der Heul | A conservative pressure-correction method for flow at all speeds[END_REF], [START_REF] Issa | The computation of compressible and incompressible flow of fluid with a free surface[END_REF], [START_REF] Kheriji | Pressure correction staggered schemes for barotropic one-phase and two-phase flows[END_REF], [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF]. Unfortunately, even if these approaches permit to bypass the consistency problem of Godunov methods, they all need to resolve the scale of the acoustic waves in the fluid in order to remain stable. This means that they suffer from a restrictive CFL (Courant-Frierichs-Levy) condition which is inversely proportional to the Mach number value. In this work, we derive a method which is able to overcome both drawbacks of standard schemes: the lack of consistency in the limit and the strong stability requirements for small Mach numbers. For sake of clarity, we start considering the isentropic/isothermal gas dynamics case and we extend the results to the case of the full Euler system. The method derived and studied in this work belongs to the class of schemes called asymptotic preserving (AP). Such type of schemes have been developed in [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] and in [START_REF] Cordier | An Asymptotic-Preserving all-speed scheme for the Euler and Navier Stokes equations[END_REF][START_REF] Degond | All speed scheme for the low mach number limit of the Isentropic Euler equation[END_REF] for the specific problem related to the compressible-incompressible passage. However, in [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF] the method is based on a Lagrange projection method and on a splitting procedure that allows to decouple the acoustic and the transport phenomenons. While in [START_REF] Degond | All speed scheme for the low mach number limit of the Isentropic Euler equation[END_REF][START_REF] Cordier | An Asymptotic-Preserving all-speed scheme for the Euler and Navier Stokes equations[END_REF], the authors split the pressure through the introduction of a numerical parameter which must be tuned depending on the problem in order for the scheme to work. Here, we use an alternative approach. Following pioneer works of [START_REF] Degond | An Asymptotic Preserving Scheme for the Euler equations in a strong magnetic field[END_REF] and of [START_REF] Cordier | An Asymptotic-Preserving all-speed scheme for the Euler and Navier Stokes equations[END_REF][START_REF] Degond | All speed scheme for the low mach number limit of the Isentropic Euler equation[END_REF], we propose a discretization based on an explicit/implicit methodology of the compressible Euler equations and a Godunov type scheme through the use of a simple Rusanov solver (see [START_REF] Eymard | Finite volume methods, Handbook of numerical analysis[END_REF], [START_REF] Rusanov | Calculation of interaction of non-steady shock waves with obstacles[END_REF]). Thanks to a stability analysis on the linearized isentropic system, we show how the numerical viscosity must be chosen in order to have a method which is consistent with the incompressible limit and stable for the chosen semi-discretization in time. Even if the method proposed in this work is only first order accurate in time and space, the proved result paves the way to high order accurate schemes both in time and space through the use of similar explicit/implicit time discretizations and of high order polynomial reconstructions for the space derivatives which employ the same numerical diffusion needed to assure the stability of the method proposed here. We show with several numerical tests in one and two dimensions that our scheme behaves as expected for different regimes ranging from very low to high Mach number flows. The last part of the paper is sacred to the extension of the scheme proposed for the isentropic/isothermal case to the full Euler equations. Numerical evidences show that the approach is fruitful and the scheme performs well on all speed flows.

The article is organized as follows. In Section 2, we present the isentropic/ isothermal system of Euler equations and its low Mach number limit. Then, in Section 3, we introduce a classical explicit semi-discretization in time and our new asymptotic preserving scheme. We perform in Section 4 a stability analysis on the linearized system around a constant solution. This analysis shows how the viscosity of the scheme must be chosen in order to ensure its uniform stability in the low Mach number limit. We present the full discretization in Section 5 and the numerical results in Section 6. Several test problems which prove the accuracy and the strong efficiency of our new method are discussed in this part. Finally, in Section 7, we extend our method to the case of the full Euler system and we test this new method on different problems. This extension is not straightforward due to non linear coupling between equations. A concluding Section ends the paper.

2. The continuous isentropic/isothermal model and its low Mach number limit. We start our discussion by considering an isentropic or isothermal flow in a bounded polygonal (or polyhedral) domain Ω ∈ R d , d = 1, 2 or 3. The extension to the full Euler system will be presented in Section 7. We denote by x ∈ Ω the space variable and by t ≥ 0 the time. Following an usual rescaling and change of variables [START_REF] Klainerman | Majda Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF][START_REF] Metivier | The incompressible limit of the non-isentropic Euler equations[END_REF], the isentropic/isothermal Euler system can be recast in a quasilinear hyperbolic system depending on the squared Mach number ε as

∂ t ρ + ∇ • (ρ U ) = 0, (2.1a) ∂ t (ρ U ) + ∇ • (ρ U ⊗ U ) + 1 ε ∇p(ρ) = 0, (2.1b)
where ρ = ρ(t, x) ∈ [0, +∞[ is the density of the fluid, U = U (t, x) ∈ R d its velocity and p(ρ) = ρ γ the pressure with γ = 1 for isothermal fluids and γ > 1 for isentropic ones. The above system is supplemented by initial and boundary conditions as for instance

U (t, x) • ν(x) = 0, for all x ∈ ∂Ω and t > 0, (2.2a) (ρ, U )(0, x) = (ρ 0 (x), U 0 (x)), for all x ∈ Ω, (2.2b)
where ν denotes the unit normal to the boundary ∂Ω, outward to Ω. We suppose that the flow to be subsonic. This means that only one boundary condition is needed at each point of the boundary. As observed in [START_REF] Klainerman | Majda Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF] well prepared or equivalently consistent initial data are necessary for system (2.1). Thus, we assume that

ρ 0 (x) = ρ 0 + ε ρ0 (x) ≥ 0, U 0 (x) = U 0 (x) + ε Ũ0 (x), ∇ • U 0 (x) = 0, (2.3)
where we stress that ρ 0 > 0 is a constant. Solutions to the isentropic/isothermal Euler system are known to exist for time intervals independent of the small parameter ε. For this problem, the low-Mach number limit is also proved rigorously to exist [START_REF] Alazard | Incompressible limit of the non isentropic Euler equations with the solid wall boundary conditions[END_REF][START_REF] Klainerman | Majda Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF][START_REF] Klainerman | Compressible and incompressible fluids[END_REF][START_REF] Metivier | The incompressible limit of the non-isentropic Euler equations[END_REF][START_REF] Schochet | The compressible Euler equations in a bounded domain: existence of solutions and the incompressible limit[END_REF]. In the next Section, we recall the formal low-Mach number limit ε → 0 which will be used in the sequel.

2.1. Low-Mach number limit for the isentropic/isothermal Euler system. When ε → 0, the momentum equation (2.1b), yields ∇ρ(t, x) = 0 and so ρ(t, x) = ρ(t), for all t ≥ 0 and x ∈ Ω. Inserting, this result in the mass equation (2.1a), integrating on Ω and using the Green formula, gives

|Ω| ρ (t) + ρ(t) ∂Ω U • ν = 0,
where |Ω| denotes the measure of Ω. Then, using the initial and boundary conditions (2.2), (2.3) and assuming that

π 1 = lim ε→0 1 ε (p(ρ) -p(ρ 0 )) < +∞,
we obtain, the incompressible isentropic Euler system

ρ = ρ 0 , (2.4a)
∇ • U = 0, (2.4b)
ρ 0 ∂ t U + ρ 0 ∇ • (U ⊗ U ) + ∇π 1 = 0, (2.4c)
with the well-prepared initial condition U (x, 0) = U 0 (x) such that ∇•U 0 (x) = 0. Note that in system (2.4) the order one correction of the pressure, π 1 , is determined thanks to the incompressibility constraint ∇ • U = 0. It is possible to derive an explicit equation for this pressure subtracting the time derivative of the incompressibility constraint to the divergence of the momentum equation. Performing these operations, we obtain

-∆π 1 = ρ 0 ∇ 2 : (U ⊗ U ), (2.5)
where ∇ 2 and : are respectively the tensor of second order derivatives and the contracted product of two tensors. Observe that equations (2.4a), (2.5), (2.4c) are equivalent to the incompressible isentropic Euler system (2.4), since (2.5) and (2.4c) yields ∂ t ∇ • U = 0 and the well prepared initial condition gives the incompressibility constraint ∇ • U = 0. We now derive the pressure wave equation which is responsible of the fast waves appearing in the Euler system. This is obtained by performing the same manipulations on system (2.1) as those done on the limit system for determining the explicit equation for the pressure π 1 . Then, we take the divergence of the momentum equation and we subtract it to the time derivative of the mass equation, this leads to

∂ 2 tt ρ - 1 ε ∆p(ρ) = ∇ 2 : (ρ U ⊗ U ). (2.6)
It is important to note that an explicit discretization of this pressure wave equation will be constrained by a stability condition of the form ∆t = O( √ ε) ∆x. This means that for being unconditionally stable with respect to the Mach number a scheme should yield an implicit discretization of the previous pressure wave equation.

3. Classical and AP semi-discretizations for the isentropic system. In the previous Section, we have seen that the low Mach number limit yields numerical difficulties for the time discretization since the pressure waves are very fast compared to the scale of the fluid motion. This may lead to think that, in order to construct a scheme which is consistent with the incompressible limit and which avoids too severe time step restrictions to remain stable, one needs to focus on the time semi-discretization of (2.1). We will see in the next section that, indeed, this is not sufficient, because the way in which the space derivatives are discretized is crucial in order to get stability and consistency with the limit equations. In particular, we will see that the numerical viscosity must be properly chosen in order to guarantee stability.

Classical explicit Euler scheme.

A standard way to discretize in time the Euler system relies on the use of explicit time integrator schemes. Among all the different possibilities, the first and simpler first order accurate scheme reads

ρ n+1 -ρ n ∆t n + ∇ • (ρU ) n = 0, (3.1a) (ρU ) n+1 -(ρU ) n ∆t n + ∇ • (ρU ⊗ U ) n + 1 ε ∇p(ρ n ) = 0. (3.1b)
The stability of such discretization, which is the less restrictive among all the possible explicit time discretizations and which will be frequently used next, is ensured under a CFL condition of type

∆t n ≤ ∆x |u i ± p (ρ)/ε| = O( √ ε) -→ ε→0 0, (3.2)
where ∆x measures the size of the space mesh. This leads to obvious stability issues because the time step ∆t n becomes tremendously small in the low-Mach number limit. Moreover, in this limit, the consistency of the discretization is not ensured because system (3.1) is reduced to ∇p(ρ n ) = 0, which is nothing but a constraint on the initial condition which may or may not be assured. Instead one would like to relate the consistency of the scheme to the model and not on the initial datum. Let us now recover the discretized pressure wave equation we get from the explicit first order time discretization. Taking the difference between equation (3.1a) at step n + 2 and n + 1, we get

ρ n+2 -ρ n+1 ∆t n+1 - ρ n+1 -ρ n ∆t n + ∇ • (ρU ) n+1 -∇ • (ρU ) n = 0.
Now, subtracting the divergence of equation (3.1b), yields

ρ n+2 -ρ n+1 ∆t n+1 - ρ n+1 -ρ n ∆t n - ∆t n ε ∆p(ρ n ) = ∆t n ∇ 2 : (ρ U ⊗ U ) n ,
which is indeed an explicit discretization of the pressure equation (2.6), constrained by a time step of order √ ε for stability to be assured.

3.2.

A new asymptotic preserving scheme for the isentropic/isothermal Euler equations. The new proposed scheme is based on an explicit/implicit time discretization strategy like in [START_REF] Degond | An Asymptotic Preserving Scheme for the Euler equations in a strong magnetic field[END_REF]. It reads

ρ n+1 -ρ n ∆t n + ∇ • (ρU ) n+1 = 0 , (3.3a) (ρU ) n+1 -(ρU ) n ∆t n + ∇ • (ρU ⊗ U ) n + 1 ε ∇p(ρ n+1 ) = 0 , (3.3b)
where the difference with respect to the full explicit scheme is that now the pressure term and the density flux are taken as implicit. Let us note that in [START_REF] Degond | All speed scheme for the low mach number limit of the Isentropic Euler equation[END_REF] a similar algorithm is proposed. However, in this previous work, the term 1/ε ∇p(ρ n+1 ) is substituted by α∇p(ρ n ) + (1/ε -α) ∇p(ρ n+1 ) where α > 0 is a numerical parameter which must be correctly fixed to ensure stability and which depends on the considered problem. Here, instead we choose α = 0. We will see that this is enough to ensure stability and preservation of the asymptotic state provided the numerical viscosity of the scheme to be properly chosen. We immediately see that this scheme is consistent for all ε ≥ 0 with systems (2.1) and(2.4). Indeed, the limit ε → 0 formally gives ∇p(ρ n+1 ) = 0, ⇒ ∇ρ n+1 = 0, which means that the scheme projects the solution over the asymptotic incompressible limit even if the initial datum is not consistent with this limit. In fact, if we operate as in the continuous case, integrating the mass equation on the domain and using the boundary condition U n+1 • ν = 0 on ∂Ω, we obtain ρ n+1 = ρ n and using the initial condition we get ρ n+1 = ρ 0 , for all n ≥ 0. Inserting this result in the mass equation, we recover the incompressibility constraint ∇ • U n+1 = 0 and then, a semi-discretization in time of the limit system (2.4)

ρ n+1 = ρ 0 , ∇ • U n+1 = 0, ρ 0 U n+1 -U n ∆t n + ρ 0 ∇ • (U ⊗ U ) n + ∇π n+1 1 = 0, where π n+1 1 = lim ε→0 1 ε p(ρ n+1
) -p(ρ 0 ) . Moreover, taking the discrete time derivative of the mass equation (i.e. the difference between equations (3.1a) at step n + 1 and n) and subtracting the divergence of the momentum equation, we obtain the following discretization for the pressure wave equation

ρ n+1 -ρ n ∆t n - ρ n -ρ n-1 ∆t n-1 - ∆t n ε ∆p(ρ n+1 ) = ∆t n ∇ 2 : (ρ U ⊗ U ) n ,
which is an unconstrained implicit discretization of the pressure wave equation (2.6).

A key observation is that the scheme even if it looks fully implicit, it can be solved partly explicitly by decoupling the mass and the momentum equations. To this aim, taking the divergence of the momentum equation (3.3b), we get

∇ • (ρU ) n+1 = ∇ • (ρU ) n -∆t n ∇ 2 : (ρU ⊗ U ) n - ∆t n ε ∆p(ρ n+1 ),
which further substituted into (3.3a) gives

ρ n+1 -ρ n ∆t n + ∇ • (ρU ) n -∆t n ∇ 2 : (ρU ⊗ U ) n - ∆t n ε ∆p(ρ n+1 ) = 0 , (3.4a) (ρU ) n+1 -(ρU ) n ∆t n + ∇ • (ρU ⊗ U ) n + 1 ε ∇p(ρ n+1 ) = 0 . (3.4b)
Thus, one can start by solving the mass equation (3.4a) which is a non linear elliptic equation determining the density at step t n+1 . Then, the momentum equation (3.4b) provides the momentum at time t n+1 explicitly. This is of great importance when dealing with multidimensional equations since the computational cost is strongly reduced compared to a standard fully implicit discretization.

4. Stability analysis for the linearized system. In this Section, we show that for the low Mach limit, it is not sufficient to choose a good time discretization to get a stable scheme. Indeed space and time discretizations are related and the choice of the space discretization is crucial in order to ensure uniform stability. We present a stability analysis for the scheme (3.3) applied to the one dimensional linearized isentropic system. In particular, we prove L 2 and L ∞ estimates depending on the choice of the viscosity coefficient of the implicit part of the flux in the momentum equation. Let us start by rewriting system (2.1) on Ω = (0, 1), denoting q = ρ U ∈ R and linearizing the system around the constant state (ρ 0 , q 0 = ρ 0 U 0 ) ∈ (0, +∞) × R . This leads to

∂ t W + A ε ∂ x W = 0,
where W = (ρ, q), and A ε = A e + A ε i , with A e and A ε i being the explicit and implicit sub-matrices related to the flux function. Denoting by c 0 = p (ρ 0 ) the sound speed (for ε = 1), these matrices are given by

A e = 0 0 -U 2 0 2 U 0 , and
A ε i = 0 1 c 2 0 /ε 0 .
The above system is supplemented by periodic boundary conditions and initial conditions which read

ρ(t = 0, x) = ρ 0 + ε ρ0 (x), q(t = 0, x) = q 0 + √ ε q0 (x).
In the following space and time meshes are chosen uniform for clarity but a more general choice will not be restrictive for our analysis. Thus, given ∆t and ∆x, the finite volume discretization of the linearized system is written as

W n+1 j -W n j ∆t + F e (W n j , W n j+1 ) -F e (W n j-1 , W n j ) ∆x + F i (W n+1 j , W n+1 j+1 ) -F i (W n+1 j-1 , W n+1 j ) ∆x =0 (4.1)
while the initial conditions are W 0 j = ρ 0 + ε ρ0,j q 0 + ε q0,j . In the following the fluxes are both discretized, the explicit and implicit part separately, with the modified Lax-Friedrichs flux. A similar analysis with the Godunov solver will give analogous results. The fluxes are then

F e (W n j , W n j+1 ) = A e W n j + W n j+1 2 -D e W n j+1 -W n j , (4.2a)
F i (W n+1 j , W n+1 j+1 ) = A ε i W n+1 j + W n+1 j+1 2 -D ε i W n+1 j+1 -W n+1 j , (4.2b)
where D e ≥ 0 and D ε i ≥ 0 are positive constants which measure the numerical viscosity and which values will be specified later. Assuming periodic boundary conditions, we also set

W n 0 = W n L , W n L+1 = W n L , for all n ≥ 0. ( 4.3)
We are now ready to prove the following result Theorem 4.1.

Let ε > 0 and (ρ 0 , U 0 ) ∈]0, +∞[×R such that U 0 /(c 0 / √ ε) ∈ [-1, 1] where c 0 = p (ρ 0 ) (linearization around a subsonic state). If the time step satisfies ∆t ≤ ∆x 2 |U 0 | . (4.4)
and

D e = |U 0 |, then scheme (4.1)-(4.3) is 1. L 2 stable, if the viscosity of the implicit part is fixed to D ε i = 0.
This means that the following estimate holds

L j=1 ∆x |α n+1 ±,j | 2 ≤ 1 2 L j=1 ∆x |α n 0,j | 2 + 1 2 L j=1 ∆x |α n 2,j | 2 , for all n ≥ 0. (4.5)
2. L ∞ stable, if the viscosity of the implicit part is fixed to

D ε i = c0 2 √
ε . This means that the following estimate holds

L max j=1 |α n+1 ±,j |≤ 1 2 L max k=1 |α n 0,k |+ 1 2 L max k=1 |α n 2,k |, for all n ≥ 0. (4.6)
Where (α n+1 -,j , α n+1 +,j ) are the coordinates in the right eigenvector basis of A ε i while (α n+1 0,j , α n+1 2,j ) are the coordinates in the right eigenvector basis of A e of W n+1 j = (ρ n+1 j , q n+1 j ) and L the number of space cells. Proof: We first need to introduce the following quantities:

λ ε ± = ± c 0 / √ ε the eigen- values of A ε i , r ε ± = (1, ±c 0 / √ ε), l ε ± = (±c 0 / √ ε, 1) the associated right and left eigen- vectors, α n+1 -,j = ρ n+1 j /2 - √ ε q n+1 j /(2 c 0 ) and α n+1 +,j = ρ n+1 j /2 + √ ε q n+1 j /
(2 c 0 ) the coordinates in the right eigenvectors basis. Similarly we introduce, λ 0 = 0 and λ 2 = 2 U 0 the eigenvalues of A e , r 0 = µ(2, U 0 ) and r 2 = (0, U 0 ) with µ a positive constant defined next the associated right eigenvectors and α n 0,j = ρ n j /(2µ) and α n 2,j = q n j -ρ n j U 0 /2 the coordinates in the right eigenvector basis (r 0 , r 2 ). The scheme (4.1), (4.2) can be rewritten separating the explicit and the implicit terms as Z n+1 j = B n j , where

Z n+1 j = W n+1 j + ∆t A ε i W n+1 j+1 -W n+1 j-1 2 ∆x - D ε i ∆t ∆x W n+1 j+1 -2 W n+1 j + W n+1 j-1 , and
B n j = W n j -∆t A e W n j+1 -W n j-1 2 ∆x + D e ∆t ∆x W n j+1 -2 W n j + W n j-1 .
Passing now to the right eigenvectors basis (r ε -, r + ε) for the implicit part leads to

Z n+1 j = α n+1 -,j 1 + 2 D ε i ∆t ∆x - ∆t ∆x α n+1 -,j+1 c 0 2 √ ε + D ε i + ∆t ∆x α n+1 -,j-1 c 0 2 √ ε -D ε i r ε - + α n+1 +,j 1 + 2 D ε i ∆t ∆x + ∆t ∆x α n+1 +,j+1 c 0 2 √ ε -D ε i - ∆t ∆x α n+1 +,j-1 c 0 2 √ ε + D ε i r ε + ,
while choosing D ε i = 0 and summing over the cells, leads to

L j=1 ∆x Z n+1 j •l ε - r ε -•l ε - α n+1 -,j = L j=1 ∆x |α n+1 -,j | 2 + c0 ∆t 2 √ ε ∆x -α n+1 -,L+1 α n+1 -,L + α n+1 -,0 α n+1 -,1 = L j=1 ∆x |α n+1 -,j | 2 ,
since the boundary conditions are periodic and α n+1 -,L+1 α n+1 -,L = α n+1 -,1 α n+1 -,0 . A similar inequality holds also true for |α n+1 +,j |. Now, performing the same passage to the right eigenvalues basis for B n j gives

B n j = α n 0,j 1 - 2 D e ∆t ∆x + D e ∆t ∆x α n 0,j+1 + D e ∆t ∆x α n 0,j-1 r 0 + α n 2,j 1 - 2 D e ∆t ∆x + ∆t 2 ∆x α n 2,j+1 (2 D e -2 U 0 ) + ∆t 2 ∆x α n 2,j-1 (2 D e + 2 U 0 ) r 2 .
Observe that for each term to be a convex combination it is enough to assure that

L j=1 ∆x |α n+1 ±,j | 2 ≤ 1 2 L j=1 ∆x |α n+1 ±,j | 2 r 0 • l ε ± r ε ± • l ε ± + r 2 • l ε ± r ε ± • l ε ± + 1 2 r 0 • l ε ± r ε ± • l ε ±   L j=1 ∆x |α n 0,j | 2 1 - 2 |U 0 | ∆t ∆x + L+1 j=2 ∆x |α n 0,j | 2 |U 0 | ∆t ∆x + L-1 j=0 ∆x |α n 0,j | 2 |U 0 | ∆t ∆x   + 1 2 r 2 • l ε ± r ε ± • l ε ±   L j=1 ∆x |α n 2,j | 2 1 - 2 |U 0 | ∆t ∆x + L+1 j=2 ∆x |α n 2,j | 2 (|U 0 |-U 0 ) + L-1 j=0 ∆x |α n 2,j | 2 (|U 0 |+U 0 )   .
Finally, since the regime is subsonic and thus

U 0 /(c 0 / √ ε) ∈ [-1, 1] we have r 0 • l ε ± r ε ± • l ε ± = µ 1 ± U 0 2 c 0 √ ε ∈ [µ/2, 3 µ/2], r 2 • l ε ± r ε ± • l ε ± = ± U 0 2 c 0 √ ε ∈ [-1/2, 1/2]. (4.7)
which choosing µ = 1/3 leads to

r 0 • l ε ± r ε ± • l ε ± + r 2 • l ε ± r ε ± • l ε ± ≤ 1.
Using, in conclusion, the periodic boundary conditions we obtain (4.5) which proves the first part of Theorem. Now, let us turn to the second part of the theorem, setting

D ε i = c 0 /(2 √ ε) gives Z n+1 j = α n+1 -,j + c 0 ∆t √ ε ∆x α n+1 -,j -α n+1 -,j+1 r ε - + α n+1 +,j + c 0 ∆t √ ε ∆x α n+1 +,j -α n+1 +,j-1 r ε + .
Denoting by

j 0 ∈ {1, • • • , L} the integer such that |α n+1 -,j0 |= max L j=1 |α n+1 -,j |= max L+1 j=0 |α n+1 -,j |, since α n+1
-,j0 and α n+1 -,j0 -α n+1 -,j0+1 have the same sign, we get for all j ∈ {1, • • • , L}

|α n+1 -,j |≤ |α n+1 -,j0 |≤ α n+1 -,j0 + c 0 ∆t √ ε ∆x α n+1 -,j0 -α n+1 -,j0+1 = |Z n+1 j0 • l -| |r -• l -| ,
where, as before, a similar inequality holds for |α n+1 +,j |. Using the expression for B n j , Z n+1 j = B n j and (4.7), the previous inequalities and remarking that

|B n j • l 0 | |l 0 • r 0 | ≤ L+1 max k=0 |α n 0,k |= L max k=1 |α n 0,k |, |B n j • l 2 | |l 2 • r 2 | ≤ L max k=1 |α n 2,k |,
we get (4.6) which concludes the proof.

5. The full Cartesian discretization for the non linear two-dimensional isentropic Euler system. In this Section, we focus on the space discretization of the two-dimensional system (2.1) discretized in time by scheme (3.3). This scheme can be recast in the uncoupled system (3.4) which has clearly advantages in terms of computational cost. This leads to consider two alternatives for the full spacetime discretization both permitting to maintain the decoupled structure of the semidiscrete scheme proposed. The first possibility is to discretize the space derivative in equations (3.3) and successively reformulate the system as done in the semi-discrete case. Alternatively, we can directly discretize equations (3.4) which naturally gives an uncoupled fully discretized scheme. Even if the first choice seems more natural, it typically yields to space discretizations with larger stencils being the full scheme obtained by two consecutive applications of the discrete Lax-Friedrichs flux functions. Instead, the second choice permits to get less diffusive schemes since the order two terms ∇ 2 : and ∆ can be directly discretized by using compact stencils. In this work, we pursue the second possibility. To this aim, we rewrite system (3.4) as

W n+1 -W n ∆t n +∇•F e (W n )+∇•F i (W n+1/2 )+    ∆p(ρ) n+1 ε + ∇ 2 : (ρ U ⊗ U ) n 0 0    = 0, where W = (ρ, ρ u, ρ v), F e (W ) = (0, ρ U ⊗ U ), F i (W n+1/2 ) = ((ρ U ) n , p n+1 /ε Id) and U = (u, v).
We consider now a rectangular domain Ω and we define by (x, y) ∈ Ω the components of the space variables. We discretize on a non uniform Cartesian mesh with L x and L y the number of cells in the two spatial directions. For all t ≥ 0 the unknown W (x, t) is approximated by the piecewise constant approximation W n k,j = (ρ n k,j , (ρ u) n k,j , (ρ v) n k,j ) with mesh steps respectively ∆x k = x k+1/2 -x k-1/2 and ∆y j = y j+1/2 -y j-1/2 with k = 1, • • • L x and j = 1, • • • L y . In this setting, the full discretized scheme reads

∆x k ∆y j W n+1 k,j -W n k,j ∆t n + (F e ) n k+1/2,j -(F e ) n k-1/2,j + (F e ) n k,j+1/2 -(F e ) n k,j-1/2 +(F i ) n+1/2 k+1/2,j -(F i ) n+1/2 k-1/2,j + (F i ) n+1/2 k,j+1/2 -(F i ) n+1/2 k,j-1/2 -∆t n ∆x k ∆y j (∆p(ρ)) n+1 k,j ε + (∇ 2 : (ρ U ⊗ U )) n k,j , 0, 0 = 0, for all k = 1, • • • , L x and all j = 1, • • • , L y .
In the above equation, the explicit fluxes considered are the Rusanov or modified Lax-Friedrichs solvers. They read

(F e ) n k+1/2,j = ∆y j F e (W n k,j ) + F e (W n k+1,j ) 2 -(D e ) n k+1/2,j W n k+1,j -W n k,j
, where (D e ) n k+1/2,j is the viscosity coefficient defined thanks to the eigenvalues of the Jacobian matrix associated to the explicit flux F e

(D e ) n k+1/2,j =    max(|u n k+1,j |, |u n k,j |), for the Rusanov solver, max k=1,•••,Lx j=1,•••,Ly |u n k,j |,
for the modified Lax-Friedrichs solver.

Similarly in the y direction we have

(F e ) n k,j+1/2 = ∆x k F e (W n k,j ) + F e (W n k,j+1 ) 2 -(D e ) n k,j+1/2 W n k,j+1 -W n k,j , with (D e ) n k,j+1/2 =    max(|v n k,j+1 |, |v n k,j |), for the Rusanov solver, max k=1,•••,Lx j=1,•••,Ly |v n k,j |,
for the modified Lax-Friedrichs solver.

The implicit fluxes are discretized by the same type of numerical fluxes used for the explicit part where, however, the choice of the numerical diffusion is done accordingly to the results of the previous section. They read, in the x direction

(F i ) n+1/2 k+1/2,j = ∆y j F i (W n+1/2 k,j ) + F i (W n+1/2 k+1,j ) 2 -(D i ) n k+1/2,j W n+1 k+1,j -W n+1 k,j
, where the viscosity coefficient can be either zero (D i ) n k+1/2,j = 0 (we expect in this case an L 2 stable scheme) or defined through the use of the eigenvalues of the Jacobian matrix associated to the implicit flux F i (we expect in this case an L ∞ stable scheme)

(D i ) n k+1/2,j = 1 2        max( p (ρ n k+1,j )/ε, p (ρ n k,j )/ε), for the Rusanov solver, max k=1,•••,Lx j=1,•••,Ly p (ρ n k,j )/ε, for the modified L-F solver.
Similarly, in the y direction

(F i ) n+1/2 k,j+1/2 = ∆x k F i (W n+1/2 k,j ) + F i (W n+1/2 k,j+1 ) 2 -(D i ) n k,j+1/2 W n+1 k,j+1 -W n+1 k,j , with (D i ) n k,j+1/2 = 0 or (D i ) n k,j+1/2 = 1 2        max( p (ρ n k,j+1 )/ε, p (ρ n k,j )/ε), for the Rusanov solver, max k=1,•••,Lx j=1,•••,Ly p (ρ n k,j )/ε, for the modified L-F solver.
Finally, concerning the discretization of the elliptic operator for the pressure, we set

(∆p(ρ)) n+1 k,j = 1 ∆x k p(ρ n+1 k+1,j ) -p(ρ n+1 k,j ) ∆x k+1/2 - p(ρ n+1 k,j ) -p(ρ n+1 k-1,j ) ∆x k-1/2 + 1 ∆y j p(ρ n+1 k,j+1 ) -p(ρ n+1 k,j ) ∆y j+1/2 - p(ρ n+1 k,j ) -p(ρ n+1 k,j-1 ) ∆y j-1/2 ,
while the last operator is discretized by two consecutive applications of a central difference scheme

(∇ 2 : (ρ U ⊗ U )) n k,j = 1 ∆x k (ρ u 2 ) n k+1,j -(ρ u 2 ) n k,j ∆x k+1/2 - (ρ u 2 ) n k,j -(ρ u 2 ) n k-1,j ∆x k-1/12 + 1 ∆x k ∆y j (ρ u v) n+1 k+1,j+1 -(ρ u v) n+1 k+1,j-1 2 - (ρ u v) n+1 k-1,j+1 -(ρ u v) n+1 k-1,j-1 2 + 1 ∆y j (ρ v 2 ) n k,j+1 -(ρ v 2 ) n k,j ∆y j+1/2 - (ρ v 2 ) n k,j -(ρ v 2 ) n k,j-1 ∆u j-1/12 ,
where ∆x k+1/2 = (∆x k+1 + ∆x k )/2 for all k = 1, • • • , L x -1, ∆x 1/2 = ∆x 1 /2 and ∆x Lx+1/2 = ∆x Lx /2. Similarly, ∆y j+1/2 = (∆y j+1 + ∆y j )/2 for all j = 1, • • • , L y -1, ∆y 1/2 = ∆y 1 /2 and ∆y Lx+1/2 = ∆y Ly /2.

6. Numerical results. In this part, we present several numerical test cases to illustrate the main features of the new method proposed, the differences and the improvements with respect to classical explicit discretization schemes as the one discussed in (3.1). The explicit scheme employed for comparisons uses the same Rusanov numerical flux of our AP method where, however, standard numerical viscosity depending on the eigenvalues of the compressible Euler equation is used. We refer to this method to as the classical scheme (CL). We test our new method against the CL method on two one dimensional problems and on a two dimensional one. For every test studied well prepared initial and boundary conditions are employed [START_REF] Degond | All speed scheme for the low mach number limit of the Isentropic Euler equation[END_REF][START_REF] Klainerman | Majda Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF]. The analysis of numerical schemes for non well prepared initial data are postponed to a future work. The CL and the AP, introduced in (3.3), schemes are run for different choices of the Mach number ranging from compressible to incompressible flows. Simulations in compressible regimes permits to observe the behavior of our scheme in situations in which explicit schemes can be employed with acceptable time step sizes. Simulations in incompressible regimes permits to observe the capability of our method to capture the limit solution without the burden of solving all micro-structures due to the presence of fast waves in the fluid. For this last situation, we also study the possibility for the AP scheme to retrieve the fast scale behavior when the time step is artificially reduced. For all problems considered, we compare the numerical results, the number of time steps needed to get the final solution and the evolution of the time step in time for the CL and the AP schemes. This permits to extrapolate the computational costs related to the two methods.

6.1. Riemann problem from Degond-Tang [START_REF] Degond | All speed scheme for the low mach number limit of the Isentropic Euler equation[END_REF]. We consider the same case test studied in [START_REF] Degond | All speed scheme for the low mach number limit of the Isentropic Euler equation[END_REF]. The initial data are

ρ 0 (x) =            1 if x ∈ [0, 0.2] 1 + ε if x ∈]0.2, 0.3] 1 if x ∈]0.3, 0.7] 1 -ε if x ∈]0.7, 0.8[ 1 if x ∈ [0.8, 1] (ρu) 0 (x) =            1 -ε/2 if x ∈ [0, 0.2] 1 if x ∈]0.2, 0.3] 1 + ε/2 if x ∈]0.3, 0.7] 1 if x ∈]0.7, 0.8[ 1 -ε/2 if x ∈ [0.8, 1].
We choose an equation of state with γ = 2, which gives p(ρ) = ρ 2 for the pressure. The final time is fixed to T f = 0.051 and 300 space cells are employed. In Figure 6.1, top panels, we report the solutions for the density on the left and the momentum on the right for the CL (red line) and the AP (blue line) schemes when the numerical viscosity is fixed equal to zero for the implicit terms in AP, i.e. D i = 0. A reference solution is also present which has been computed by the CL scheme with 3000 spatial cells. We observe that our scheme is able to reproduce the same results of the CL scheme with a slightly larger numerical diffusion due to its semi-implicit character. In addition, we clearly see some small overshoots after the shock on the left part of the domain (around x 0.14). These overshoots not only remain bounded but their intensity diminishes in time. They are linked to the choice of D i = 0. This result is expected since our stability result affirms that we can only assure L 2 stability with D i = 0 and not L ∞ stability. Figure 6.1 bottom panels presents the results of the same simulation where now the numerical diffusion of the implicit terms is not zero anymore, instead it has been fixed proportional to the Mach number as demanded to guarantee L ∞ stability. As before, we report the results for the density on the left and the momentum on the right. As expected, since the diffusion now is larger some additional dissipation is observed in the AP scheme but the spurious overshoots have completely disappeared, the scheme is L ∞ stable. The price to pay in order to have a scheme which assures a stronger stability results is that the momentum equation requires the resolution of a linear system since the viscosity terms are implicit which may become costly in higher spatial dimensions.

In Figure 6.1 are also reported the values of the time steps ∆t n for the the CL and the AP schemes (with and without extra implicit dissipation). We observe that for this relatively large Mach numbers, our scheme already demands smaller time steps compared to explicit discretizations. On the other hand, the time steps restrictions for the L 2 and the L ∞ schemes are almost the same. In Figure 6.2, we report the behaviors of our scheme when the mesh is refined. In particular, we show simulation results for 300, 900 and 1500 cells. We observe that the numerical solution does converge towards the reference solution as expected. We focus now on a low Mach regime. In Figure 6.3 are presented the results for the two schemes, CL and AP, when a 300 cell mesh is employed. Blue line for the AP, red line for CL, reference solution black line. The reference solution has been computed by using the CL scheme with 3000 cells. On the left top panel the density is reported while on the top right panel the momentum. In this case, the Mach number is 10 -4 which gives ε = 0.01. The asymptotic stable scheme projects the solution to the limit incompressible solution faster than the explicit scheme. This is exactly the scope of the method, we want to avoid to capture the micro-structures of the solution, instead the goal is to capture the limit solution at reasonable costs. To this aim, we report the size of time steps in the bottom panel of Figure 6.3. The time steps of the AP scheme are around seventy time smaller than the ones required by an explicit method. This shows that the AP scheme can employ time steps which are independent with respect the small parameter ε while time steps of explicit schemes remain close to the Mach number values. In order to prove that the AP scheme converges to the correct solution we report in Figure 6.4 the same results as in Figure 6.3 but with a small time step for the AP scheme. The AP scheme employs 604 time steps to capture this new solution shown in the Figure while the CL scheme 684 time steps. From the results we can conclude that the AP scheme is able to compute the correct rate of We finally measure in table 6 As expected, we observe that in the case of large ε the cost of our AP scheme is more important that the classical one, even if the costs are still comparable. This is due to the time needed to invert the linear system which characterize the AP scheme. However, when ε becomes small enough then the AP scheme only needs 11 cycles to get to the final time whereas the number of iterations for the CL scheme increases indefinitely. The CPU time needed to update one AP cell in one time step is about constant of the order 4.265 × 10 -6 s whereas for the CL scheme is obviously less expensive and of the order of 2.5 × 10 -7 with our implementation.

Time varying Mach number test.

In this part we consider a Riemann problem where the Mach number changes with time. In these cases, if one is interested only to the main features described by the limit incompressible solution and not to the fast scale dynamics which ends in the limit problem, the use of an Asymptotic Preserving scheme is indicated since it permits to bypass the solution of the fast waves without losing the asymptotic behaviors. In Figure 6.5, on the left panel, is shown the value of ε(t) as a function of time in logarithmic scale. The schemes, CL and AP, use both 600 spatial cells to furnish a solution. On the right panel is shown the time steps needed by the two methods as a function of time in logarithmic scale. As expected while the explicit scheme has a time step which follows the variation of ε(t), the AP scheme is able to maintain a relative constant time step throughout the simulation. The total number of time steps needed by our scheme is 163 while an explicit method needs 24897 steps in order to remain stable. In Figure 6 Up to time t = 0.05 the AP scheme still captures the main structures of the solution but with additional dissipation (first two rows). However, we see that a larger time step can be used compared to CL. When ε jumps to small values (third row) then the AP scheme bypasses the solution of the small scale dynamics and starts to project the solution to the limit incompressible solution but still it remains able to capture the main phenomena. Finally, when ε decreases to 0, our scheme captures the limit solution which is 1, while the CL scheme still tries to solve the small scale paying the price of extremely small time steps.

Double rarefaction test problem.

In this test we measure the errors produced by our scheme. To this aim, we have computed an exact solution of the isentropic Euler equations composed by two rarefaction waves emanating from the following initial condition

ρ 0 (x) = 1 + ε if x ∈ [0, 0.5], 1 if x ∈]0.5, 1.0], (ρu) 0 (x)
= (1 + ε)(1 - √ ε) if x ∈ [0, 0.5], 1 + √ ε if x ∈]0.5, 1.0].
We choose an equation of state with γ = 2, which gives p(ρ) = ρ 2 for the pressure. The structure of the exact solution is constituted of left and right moving rarefaction waves separated by a constant state. We choose three different values for the squared Mach number, i.e. ε = 0.99, 0.1, 0.001. For these three values we compute the solution up to final times 0.1, 0.05 and 0.007 respectively for the AP and the CL schemes employing 1000 spatial cells on the domain (0, 1). The results are provided in Figure 6.7 against the exact solution. The ratio between the number of cycles needed by the two schemes are respectively 984/798, 579/264, 641/30 for the three cases that is about 1. of ε = 0.99, 0.1, 0.001 as a function of the ratio between the time step used and the maximum time step allowed for stability to be guaranteed (CFL in the Figure ). We clearly observe that while the error for AP scheme is reduced when the time step decreases, this is not the case for the explicit scheme for any value of the Mach number. Fig. 6.8. Double rarefaction test case. 1000 cell simulations. L 2 norms of the error for the density for the CL and the AP schemes as a function of the ratio between the time step used and the maximum time step allowed for ε = 0.99, 0.1, 0.001. Red line explicit scheme, blue line AP scheme.

2D isentropic cylindrical explosion problem.

In this section we focus on a genuinely 2D cylindrical problem with γ = 1. The computational domain is set to Ω = [0, 1] 2 paved with L x × L y uniform cells. The initial density is given by

ρ k,j = 1 + ε if r 2 k,j ≤ 1/4 1 else ,
where r k,j is the distance from the cell center. The velocity is set as u k,j = -α k,j x i /(r k,j ρ k,j ), v k,j = -α k,j y j /(r k,j ρ k,j ), where the coefficients α k,j are given by α k,j = max(0, 1-r k,j )(1-e -16r 2 k,j ) and we set u k,j = v k,j = 0 if r k,j ≤ 10 -15 , see figure 6.9 for a representation of the initial data. Periodic boundary conditions are set everywhere. We take 200 cells in both directions. We compute the solution for different values of the squared Mach number ε by using our AP scheme and the CL scheme for ε = 1 and ε = 10 -4 . The AP scheme uses D i = 0 which means that we consider the L 2 stable scheme. In Figures 6.10 field for the CL and the AP schemes (top and bottom panels respectively) for ε = 1.0 at t = 0.1 (short time), t = 0.24 (focusing of waves) and t final = 0.5 the final time. From these Figures, we can observe that the AP scheme can reproduce the numerical results of a classical explicit scheme in the case ε = 1. Concerning the time steps, the explicit scheme demands 144 iterations while the AP scheme needs only 64 iterations to get to the final solution. However the AP scheme scheme being more expensive due to its implicit character, the ratio of CPU time is still of order 3 in favor of the explicit method. Finally, in Figures 6.12 are presented the density variable and the velocity field for the CL and AP schemes (top and bottom panels respectively) for ε = 10 -4 , i.e. close to the incompressible limit for a final time t final = 0.05. In this case, the classical explicit scheme demands 1109 time steps while the AP scheme only needs 3 time steps to get to the final solution. From the computational point of view, the AP scheme is 5 time less expensive in our implementation. From the images, we see that while the AP scheme has captured the limit solution the explicit scheme is still resolving small waves dynamic. The difference in time steps between the two schemes being of the order of the Mach number √ ε. 7. Extension to the full Euler system. In this Section we study the extension of our scheme to the case of the full Euler system. This system in rescaled variables reads

∂ t ρ + ∇ • (ρ U ) = 0, (7.1a)
∂ t (ρ U ) + ∇ • (ρ U ⊗ U ) + 1 ε ∇p = 0, (7.1b)
∂ t E + ∇ • ((E + p) U ) = 0, (7.1c)
supplemented with the following equation of state

p = (γ -1) E -ε ρ|U | 2 2 , (7.1d)
where γ ≥ 1 is the ratio of specific heats. Like for the isentropic case, appropriate initial and boundary values should be provided U (t, x) • ν(x) = 0, for all x ∈ ∂Ω and t > 0, (

(ρ, U, p)(0, x) = (ρ 0 (x), U 0 (x), p 0 (x)), for all x ∈ Ω, (7.1f) where the initial condition is well prepared to the low-Mach number regime p 0 (x) = p 0 +ε p0 (x) ≥ 0, and

U 0 (x) = U 0 (x)+ √ ε Ũ0 (x), ∇•U 0 (x) = 0, (7.1g)
with p 0 > 0 a given constant. In this case, the low Mach number limit gives the incompressible Euler system, see [START_REF] Metivier | The incompressible limit of the non-isentropic Euler equations[END_REF], [START_REF] Alazard | Incompressible limit of the non isentropic Euler equations with the solid wall boundary conditions[END_REF] ∂

t ρ + ∇ • (ρ U ) = 0, ∇ • U = 0, ∂ t (ρ U ) + ∇ • (ρ U ⊗ U ) + ∇π 1 = 0,
with the pressure given by p(t, x) = (γ-1) E(t, x) = p 0 , where π 1 = lim ε→0 1 ε (p -p 0 ) < +∞ is implicitly defined by the constraint ∇ • U = 0 and explicitly given by the equation -∆π 1 = ρ 0 ∇ 2 : (U ⊗ U ).

7.1. The semi AP discretization in time. The scheme for the full Euler system follows the isentropic case strategy. The direct extension to this previous scheme to the full Euler case reads

ρ n+1 -ρ n ∆t + ∇ • q n = 0, q n+1 -q n ∆t + ∇ • q ⊗ q ρ n - γ -1 2 ∇ |q| 2 ρ n + γ -1 ε E n+1 = 0, E n+1 -E n ∆t + ∇ • γ E q ρ n+1 -∇ • (γ -1) ε 2 |q| 2 q ρ 2 n = 0, p n+1 = (γ -1) E n+1 -ε (ρ u 2 ) n 2 .
We observe that in [START_REF] Cordier | An Asymptotic-Preserving all-speed scheme for the Euler and Navier Stokes equations[END_REF] a similar algorithm is proposed but the term 1/ε ∇E n+1 is changed by α∇E n + (1/ε -α) ∇E n+1 where α > 0 is a numerical parameter which must be correctly fixed to ensure the uniform stability. Here, like in the isentropic case, we choose α = 0 and we carefully choose the numerical space viscosity of the scheme in order to be stable. However, unfortunately, this approach has the main drawback that the Jacobian matrix associated to the explicit flux may have complex eigenvalues. In order to avoid this problem, we choose instead the following non conservative discretization

ρ n+1 -ρ n ∆t + ∇ • q n+1 = 0, q n+1 -q n ∆t + ∇ • q ⊗ q ρ n - γ -1 2 ∇ |q| 2 ρ n + γ -1 ε ∇E n+1 = 0, E n+1 -E n ∆t + ∇ • γ E n q n+1 ρ n -∇ • (γ -1) ε 2 |q| 2 q ρ 2 n = 0, p n+1 = (γ -1) E n+1 -ε (ρ u 2 ) n 2 .
Now, this semi-discretization can be rewritten in a compact form as

W n+1 -W n ∆t + A e (W n , W n+1 )∂ x W n + A i (W n , W n+1 ) ∂ y W n+1 = 0,
where A e and A i have real eigenvalues. This scheme, as the one proposed for the isentropic system, has an uncoupled formulation. This can be obtained by inserting the value of q n+1 given by the momentum equation into the energy equation. This gives a linear parabolic equation for the energy at time n + 1 which furnishes the energy E n+1 . Then, using the knowledge of the energy, one can update the momentum equation to get q n+1 and finally the mass equation gives ρ n+1 . This decoupled formulation is of fundamental importance when dealing with multidimensional problems since it avoids the inversion of very large linear systems.

Numerical tests.

In this last part, we present two numerical test cases to illustrate the behavior of this new scheme for the full Euler system in the one dimensional setting. We use the same Rusanov numerical fluxes employed for the isentropic case and a CFL condition related to the eigenvalues of the explicit Jacobian matrix A e . This gives a time step limitation of type

∆t n ≤ ∆x max |(3 -γ) U n |, γ |U n | . 7.
2.1. Riemann like problem. In this problem, the computational domain is set to Ω = (0, 1) with periodic boundary conditions on both sides. Density and pressure are uniforms and set to 1 with γ = 7/5. The velocity is defined by u(x) = 1 -ε/2 if x ≤ 0.2 or x ≥ 0.8, 1 + ε/2 if x ≤ 0.75 or x ≥ 0.25 and 1 elsewhere. The total energy E is deduced from the equation of state when ε is fixed as p

= (γ -1) E -ε 2 ρu 2 2
. The computations are stopped at the final time t f inal = 0.05 and the square of the Mach number ε is chosen equal to 0.99 and 10 -4 . In Figure 7.1 are depicted the pressure and velocity variables for the explicit classical (first order explicit Euler scheme) and AP schemes when 300 cells are employed. A reference solution is also present which has been obtained with the CL scheme and 10000 cells. From these Figures, we observe that both schemes produce rather equivalent numerical solutions although the AP scheme is a little bit more dissipative. The AP scheme demands 40 cycles while the CL one 45. The CPU time is about 0.08 for the AP scheme and 0.02 for the CL. The right panel of Figure 7.1 presents the time steps as a function of time in logarithmic scale for both schemes showing that, as expected, ∆t is of the same order in the case ε = 0.99. In Figure 7 Riemann problem for the one dimensional Euler system for ε = 0.99 at t final = 0.05 for 300 cells. Left panel show the pressure profile, middle panel show the velocity profile. CL scheme is depicted in red while the AP scheme is depicted in blue. Reference solution is in black. result as in Figure 7.1 with a Mach number squared of ε = 10 -4 . The AP scheme demands only 27 cycles to furnish a solution while the CL 2177 time steps. The CPU time is about 1.02 for the CL scheme and 0.07 for the AP one. Moreover, we can observe in the panels that the AP scheme has captured the limit incompressible solution while the explicit scheme is still resolving the small scale dynamic. The time step employed by the AP scheme is of the order ∆t = 1.9×10 -3 , while the CL scheme is restricted to ∆t = 2.3 × 10 -5 leading to a ratio of √ ε = 10 -2 for the two time steps. We finally consider the solution of the Riemann problem in which our scheme is used with a restricted time step, i.e. ∆t to 10 -5 . In this case, we observe that we are able to retrieve the small scale structures of the solution lost in the case of larger time steps. The results are reported in Figure 7.3 where a comparison between the explicit results and those of the AP scheme is shown. In this case the AP scheme demands 5000 cycles for a CPU time of the order of 6.2.

Colliding acoustic pulses.

The last test we consider is a weakly compressible problem from [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] Arun | An asymptotic preserving all Mach number scheme for the Euler equations of gas dynamics[END_REF]. The computational domain is Ω = (-L, L) with Riemann problem for the one dimensional Euler system for ε = 10 -4 at t final = 0.05 for 300 cells. Left panel show the pressure profile, middle panel show the velocity profile. CL scheme is depicted in red while the AP scheme is depicted in blue. Reference solution is in black. In Figure 7.4 we report the results in term of density, velocity and energy at an intermediate time t = 0.815 and at the final time t = 1.63 along with the initial data for the CL and the AP schemes using 440 cells and ε = 1/11. We can observe that the AP scheme is a little more dissipative but capture the main physical phenomena. The CL scheme needs 1030 cycles to reach the final time while the AP scheme only demands 224 time steps leading to a ratio of about 4.6.

8. Conclusions. The purpose of this article is the development and the analysis of a new Asymptotic Preserving scheme which is able to deal with low Mach number fluid flows bypassing the resolution of the fast waves dynamics. We have proven the good behaviors of our scheme by performing a stability analysis for the isentropic Euler equations. We have successively extended the proposed scheme to the case of the full Euler system. The numerical results show that the scheme possesses the demanded properties. It is L 2 stable or L ∞ stable depending on the choice of the numerical diffusion and consistent with the low Mach number limit. Close to incompressible regimes this new scheme permits to compute solutions of low Mach number flows at a much lower computational costs. This reduction being much more important when dealing with multidimensional problems. However, results also indicate that the scheme in intermediate regimes is more diffusive than explicit methods. For this reason, we plan to develop in future works high order schemes in time and space which will permit to get more accurate results without losing the nice properties of the scheme here presented.

1 - 2

 12 D e ∆t ∆x ≥ 0, and 2 D e ≥ 2 |U 0 |. Choosing finally D e = |U 0 |, leads to condition (4.4) while summing over the cells Z n+1 j = B n j and using the Cauchy-Schwarz and the Young inequalities, we obtain

Fig. 6 . 1 .Fig. 6 . 2 .

 6162 Fig. 6.1. Solution of the Riemann problem at t = 0.051 for 300 cells and ε = 0.99. Left panel density ρ, middle panel momentum (ρu), right panel time steps ∆t as a function of time. Reference solution in black line, classical explicit scheme CL in red and AP scheme in blue. Top panels: AP scheme without extra implicit dissipation (D i = 0). Bottom panels: AP scheme with implicit dissipation (D i = 0).

Fig. 6 . 4 .

 64 Fig. 6.4. Solution of the Riemann problem at t = 0.008 for ε = 0.01 and 300 cells for the AP and the CL schemes. Left panel density ρ, right panel momentum (ρu). Reference solution in black line, CL red line and AP scheme blue line. The time step of the AP scheme is reduced to the order of the time step of the explicit scheme.

Fig. 6 . 5 .

 65 Fig. 6.5. Riemann problem with time varying Mach number at t = 0.06 for 600 cells. Left panel shows the value of ε(t) as a function of time. Right panel shows time steps as a function of time. The scales are logarithmic. The number of time steps is 163 for the AP scheme and 24897 for the classical explicit scheme.

Fig. 6 . 6 .

 66 Fig. 6.6. Time varying Mach number test. CL (red line) and AP (blue line), 600 cells. Snapshots of the density at different times. From left to right: eight different intermediate times corresponding to different (t) values (straight vertical green line on the embedded panels).

Fig. 6 . 7 .

 67 Fig. 6.7. Double rarefaction test case. Density as a function of x at final time, 1000 spatial cells. Left panel ε = 0.99, central panel ε = 0.1, right panel ε = 0.001. Red line classical explicit scheme, blue line AP scheme, black line exact solution.

  -6.11 are presented the density variable and the velocity

. 9 -Fig. 6 . 9 .

 969 Fig. 6.9. 2D cylindrical problem for ε = 1. Initial density left panel, (ρu) and (ρv) as a function of the distance from the center middle panel, two dimensional velocity field right panel.

Fig. 6 .

 6 Fig. 6.10. 2D cylindrical problem for ε = 1 at t = 0.1 left panel, t = 0.24 middle panel and t final = 0.5 right panel. Top images show the density profile for the CL scheme. Bottom images show the density profile for the AP scheme.

Fig. 6 . 11 .

 611 Fig. 6.11. 2D cylindrical problem for ε = 1 at t = 0.1 left panel, t = 0.24 middle panel and t final = 0.5 right panel. Top images show the velocity field for the CL scheme. Bottom images show the velocity field for the AP scheme.

Fig. 6 . 12 .

 612 Fig.6.12. 2D cylindrical problem for ε = 10 -4 at t final = 0.05. Top panel density for the CL and the AP schemes. Bottom panel velocity field.

Fig. 7 .

 7 Fig. 7.1.Riemann problem for the one dimensional Euler system for ε = 0.99 at t final = 0.05 for 300 cells. Left panel show the pressure profile, middle panel show the velocity profile. CL scheme is depicted in red while the AP scheme is depicted in blue. Reference solution is in black.

Fig. 7 .

 7 Fig. 7.2.Riemann problem for the one dimensional Euler system for ε = 10 -4 at t final = 0.05 for 300 cells. Left panel show the pressure profile, middle panel show the velocity profile. CL scheme is depicted in red while the AP scheme is depicted in blue. Reference solution is in black.

Fig. 7 . 3 .

 73 Fig. 7.3.Riemann problem for the one dimensional Euler system for ε = 10 -4 at t final = 0.05 for 300 cells. Left panel show the pressure profile, right panel show the velocity profile. CL scheme is depicted in red while the AP scheme is depicted in blue. Reference solution is in black. The AP scheme employs a fixed time step ∆t = 10 -5 whereas the Classical scheme employs a stability restricted time step like in the previous figure.

L = 2 ε

 2 . The initial setup consists of two colliding acoustic pulses in a weakly compressible regimeρ(x, 0) = ρ 0 + ε 2 ρ 1 α(x), u(x, 0) = -u 0 2 sign(x) ερ 1 α(x), p(x, 0) = p 0 + ε 2 p 1 α(x),where α(x) = 1 -cos 2πx L , ρ 0 = 0.955, ρ 1 = 2, u 0 = 2 √ γ, p 0 = 1 and p 1 = 2γ.

Fig. 7 . 4 .

 74 Fig. 7.4. Colliding acoustic pulses problem for the one dimensional Euler system for ε = 1/11 at t = 0.815 (top panels) and t = 1.63 (bottom panels) for 440 cells. Left panel show the density profile, middle panel shows the velocity profile, right panel show the energy profile. CL scheme is depicted in red while the AP scheme is depicted in blue. Initial condition is in black.

  Solution of the Riemann problem at t = 0.008 for ε = 0.01 and 300 cells for the AP and CL schemes. Left panel density ρ, middle panel momentum (ρu), right panel time steps ∆t as a function of time. Reference solution in black line, CL scheme red line and AP scheme blue line.convergence to the limit solution if the time step is reduced and made closer to the time step of the explicit method.

	1.00005					1.008
											Classical 300 cells			AP sch.
											AP 300 cells			Classical sch.
	1.00004										Reference
						1.006						0.001
	1.00003
						1.004														8.30e-4
	1.00002
						1.002
	1.00001
	1.00000					1
												0.0001
	0.99999
						0.998
	0.99998																			1.17e-5
						0.996
	0.99997
						0.994
	0.99996	Classical 300 cells
		AP 300 cells
		Reference
	0.99995					0.992						1e-05
	0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1	0	0.001	0.002	0.003	0.004	0.005	0.006	0.007	0.008	0.009
		0.99996 0.99997 0.99998 0.99999 1.00000 1.00001 1.00002 1.00003 1.00005 Fig. 6.3. 0.99995 1.00004	Classical 300 cells AP 300 cells -CFL=0.008 Reference					0.992 0.994 0.996 0.998 1 1.002 1.004 1.006 1.008					Classical 300 cells AP 300 cells -CFL=0.008 Reference
			0	0.2	0.4		0.6	0.8	1	0	0.2	0.4		0.6		0.8		1

  .1 the number of time steps and the CPU times needed to reach the final time for different values of the Mach number, from compressible to incompressible situations.

		Mach #	# time steps CL AP	Ratio CL/AP	CPU time CL AP	Ratio CL/AP
		ε = 1	150	245	0.6	0.076s	0.12 s	0.63
		ε = 0.9	51	94	0.5	0.066 s	0.0083 s	0.8
		ε = 0.75	24	39	0.6	0.035 s	0.0060 s	0.58
		ε = 0.5	18	20	0.9	0.031 s	0.0033 s	0.94
	300 cells	ε = 0.25 ε = 10 -1 ε = 10 -2 ε = 10 -3 ε = 10 -4	21 28 74 220 684	16 14 11 10 11	1.3 2.0 6.7 22 62	0.031 s 0.039 s 0.078 s 0.085 s 0.092 s	0.0025 s 0.0033 s 0.0057 s 0.0055 s 0.0054 s	1.2 1.2 1.4 1.6 1.7
		ε = 10 -5	2152	11	196	0.29 s	0.0051 s	5.7
		ε = 10 -6	6797	11	618	0.83 s	0.0027 s	30.7
		ε = 10 -7	21472	11	1952	2.58 s	0.0024 s	108
		ε = 10 -8	67688	10	6768	8.10 s	0.0026 s	311.5

Table 6 . 1

 61 Number of time steps and CPU time needed for the CL and the AP schemes to reach t = 0.008 for the Riemann problem for different values of ε.

  23, 2.19, 21.4. The L 2 norms of the errors are displayed in Figure6.8 for the three values

	2					1.1						1.001
					Classical CFL=0.5						Classical CFL=0.5						Classical CFL=0.5
					AP CFL=0.5						AP CFL=0.5						AP CFL=0.5
					Reference						Reference						Reference
	1.8					1.08						1.0008
	1.6					1.06						1.0006
	1.4					1.04						1.0004
	1.2					1.02						1.0002
	1					1						1
	0.8					0.98						0.9998
	0.6					0.96						0.9996
	0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1

  .2 are reported the same

		Ref. CL 10000				1.6		Ref. CL 10000
		CL 300 cells -45 cycles						CL 300 cells -45 cycles
		AP 300 cells -40 cycles						AP 300 cells -40 cycles
	1.6					1.4
	1.4					1.2
	1.2					1
	1					0.8
	0.8					0.6
	0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1

Acknowledgments. -R. Loubère has been partially funded by the ANR under the JCJC project "ALE INC(ubator) 3D". G. Dimarco and M.H. Vignal were supported by the French ANR project MOONRISE. The material of this work has been gathered during SHARK-FV conference (Sharing Higherorder Advanced Know-how on Finite Volume www.math.univ-toulouse.fr/SHARK-FV/ held in 2014.