Low Mach number limit of some staggered schemes for compressible barotropic flows
R. Herbin J. -C Latché Khaled Saleh

Homa: A Receiver-Driven Low-Latency

Transport Protocol Using Network Priorities (Complete Version)

    
Behnam Montazeri, Yilong Li, Mohammad Alizadeh†,
and John Ousterhout

Stanford University, †MIT


This paper is an extended version of the paper on Homa that was published in ACM
SIGCOMM 2018.

Material had to be removed from Sections 5.1 and 5.2 to meet the
SIGCOMM page restrictions; this

version restores the missing material.

    

================================================================================================================================================================================================================================================================================================================================================


 December 30, 2023



In this paper, we study the behaviour of some staggered discretization based numerical schemes for the barotropic Navier-Stokes equations at low Mach number.
Three time discretizations are considered: the implicit-in-time scheme and two non-iterative pressure correction schemes.
The two latter schemes differ by the discretization of the convection term: linearly implicit for the first one, so that the resulting scheme is unconditionally stable, and explicit for the second one, so the scheme is stable under a CFL condition involving the material velocity only.
We rigorously prove that these three variants are asymptotic preserving in the following sense: for a given mesh and a given time step, a sequence of solutions obtained with a sequence of vanishing Mach numbers tend to a solution of a standard scheme for incompressible flows.
This convergence result is obtained by mimicking the proof of convergence of the solutions of the (continuous) barotropic Navier-Stokes equations to that of the incompressible Navier-Stokes equation as the Mach number vanishes.
Numerical results performed with a hand-built analytical solution show the behaviour that is expected from the analysis.



Compressible Navier-Stokes equations, low Mach number flows, finite volumes, Crouzeix-Raviart scheme, Rannacher-Turek scheme, finite elements, staggered discretizations.


 35Q31,65N12,76M10,76M12 








§ INTRODUCTION
 

We consider the non-dimensionalized system of time-dependent barotropic compressible Navier-Stokes equations, parametrized by the Mach number denoted thereafter by , and posed for (,t)∈Ω× (0,T): 


    ∂_t ρ^ + ( ρ^ ^) = 0,
    
       ∂_t (ρ^ ^) + (ρ^ ^⊗^) 
     - ((^)) + 1/^2 (ρ^) = 0,
 
where T is a finite positive real number, and Ω is an open bounded connected subset of ^d, with d∈{2,3}, which is polygonal if d=2 and polyhedral if d=3.
The quantities ρ^>0 and ^=(u_1^,..,u_d^)^T are the density and velocity of the fluid.

In the isentropic case, the pressure satisfies the ideal gas law (ρ^) = (ρ^)^γ, with γ≥ 1, the heat capacity ratio, a coefficient which is specific to the considered fluid.
However, more general barotropic cases can be considered provided the equation of state   is a C^1 increasing convex function  such that '(1) >0, see Remark <ref>.

Equation (<ref>) expresses the local conservation of the mass of the fluid while equation (<ref>) expresses the local balance between momentum and forces. 
We consider Newtonian fluids so that the shear stress tensor (^) satisfies:

    (()) = μΔ+ (μ+λ) ( ),

where μ and λ are two parameters with μ>0 and μ+λ>0.
System (<ref>) is complemented with the following boundary and initial conditions:

    ρ^|_t=0 = ρ_0^,         ^|_t=0=_0^,         ^|_Ω=0.

At the continuous level, when  tends to zero, the density ρ^ tends to a constant and the velocity tends, in a sense to be defined, to a solution of the incompressible Navier-Stokes equations <cit.>.
Heuristically, the momentum equation (<ref>) indicates that ρ^ behaves like ρ̅(t) + O(^2/γ) where ρ̅(t) is a function only dependent on the time variable.
Integrating the mass conservation equation (<ref>) over Ω and using the homogeneous Dirichlet boundary condition (<ref>) (a homogeneous Neumann condition would be sufficient) then implies that ρ̅ is actually a constant.
For such a result to hold, some assumptions need to be made on the initial data; in particular, the initial density ρ_0^ must be assumed to be close to ρ̅ in a certain sense.
These assumptions will be specified below.
Setting ρ̅=1 (here and throughout the paper) without loss of generality, passing to the limit in the mass conservation equation (<ref>) and in the momentum balance (<ref>), the limit velocity  is formally seen to solve the system of incompressible Navier-Stokes equations:

    ( ) = 0, 
       _t + (⊗)- μ+ π = 0,

where π is the formal limit of ((ρ^)-1)/^2.
This formal computation was justified by rigorous studies <cit.>; see also <cit.> for some of the first mathematical analyses on low Mach number limits and  <cit.> for some of the numerous related works.

For the low Mach limit of numerical schemes for (<ref>), the issue is not so clear.
Indeed, in general, schemes designed for computing compressible flows do not boil down, when → 0, to standard schemes for incompressible flows, for essentially two reasons.
First, the numerical dissipation introduced to stabilize the scheme depends on the celerity of acoustic waves, which blows up when → 0; a reasonable approximation of the incompressible solution thus may need a very small space step, depending not only on the regularity of the continuous solution but also on .
Second, these schemes are usually explicit in time (with a sophisticated derivation of fluxes, for instance through solutions of Riemann problems at interfaces, which lead to a nonlinear expression with respect to the unknowns), which is not compatible with a wave celerity blowing up in the incompressible limit; to cope with the low Mach number situation, an implicitation of some terms in the equations, usually the pressure gradient in the momentum balance equation and the mass fluxes divergence in the mass balance, is thus necessary, and this makes implicit-in-time discrete analogues of the wave equation for the pressure appear.
Unfortunately, the schemes for the compressible case usually use a collocated arrangement of the unknowns (specially if one intends to compute the numerical fluxes on the basis of the solution of a Riemann problem at faces, which is well suited to a cell by cell piecewise constant approximation), and the diffusion operator appearing in this wave equation (obtained by a discrete composition of the pressure gradient and the velocity divergence operators) is unstable, since collocated approximations do not satisfy a form of the so-called discrete inf-sup condition.
At the incompressible limit, the scheme will thus need an additional stabilization mechanism <cit.>. 
These phenomena have been widely studied and corrections have been proposed <cit.>.
To obtain a scheme accurate for all Mach number flows, an alternative route consists in starting from technologies initially designed for  the incompressible Navier-Stokes equations and extending them to compressible flows.
This approach may be traced back to the late sixties, when first attempts were done to build "all flow velocity" schemes <cit.>; these algorithms may be seen as an extension to the compressible case of the celebrated MAC scheme, introduced some years before <cit.>.
These seminal papers have been the starting point for the development of numerous schemes falling in the class of pressure correction algorithms (see  <cit.> for a presentation in the incompressible case), possibly iterative, in the spirit of the SIMPLE method, some of them based on staggered finite volume space discretizations <cit.>; a bibliography extended to the schemes using other space discretizations may be found in <cit.>.

In this paper, we address, besides a purely implicit scheme, variants of schemes falling in this latter class, namely non-iterative pressure correction schemes based on staggered discretizations.
These schemes have been developed in the last ten years, first for barotropic Euler and Navier-Stokes equations <cit.> and then for the non-barotropic case <cit.>, and have been shown both theoretically and numerically to be consistent and accurate for the Navier-Stokes and Euler equations and for Mach numbers in the range of unity (including shock solutions in the inviscid cases).
In addition, some numerical experiments <cit.> suggest that, when the Mach number tends to zero, the numerical solution tends to the solution of a standard scheme for incompressible or, in non-isothermal situations, quasi-incompressible flows (in the sense of the classical asymptotic model for low Mach numbers <cit.>).
Our aim here is to rigorously prove that, in the barotropic case, these schemes are indeed asymptotic preserving: for a given discretization ( mesh and time step), when the Mach number tends to zero, the solution tends to the solution of a standard (stable and accurate) scheme for incompressible flows.
To this purpose, we reproduce at the discrete level the analysis performed in the continuous case in <cit.> (of course, with heavy simplifications on compactness arguments, especially concerning the compactness of the sequence of discrete velocities, since we work in a finite dimensional setting); to our knowledge, this is the first presentation of such a proof.
We draw the reader's attention of the fact that the results are mainly presented in the case of the Navier-Stokes equations, but we also explain how they can be easily extended to the inviscid case of the Euler equations where μ=λ=0.  

We address three different time-discretizations: first, a fully implicit scheme, because the convergence proof in this case is simpler and necessitates less restrictive assumptions on the initial data; then we turn to two variants which are more efficient in practice, namely two pressure correction schemes, which differ by the discretization of the convection term in the momentum balance equation, linearly implicit for the first one (so that the corresponding scheme is unconditionally stable) and fully explicit for the second one (so that the corresponding scheme is stable under a CFL condition based on the material velocity).
The paper is organized as follows.
First, for the reader's convenience, we recall the (part of) the continuous analysis which is mimicked at the discrete level (Section <ref>).
Then we define the meshes and unknowns used by the schemes (Section <ref>) and the space discretization of the operators involved in the balance equations (Section <ref>).
The next three sections are devoted to the convergence analysis for the time-discretization variants: the implicit scheme (Section <ref>) and the two pressure correction schemes (Sections <ref> and <ref>).
The last section presents some numerical results for the pressure correction algorithm at various Mach numbers on a problem which is built upon a hand-built analytical solution.






§ INCOMPRESSIBLE LIMIT IN THE CONTINUOUS SETTING


The convergence when  tends to zero of a weak solution (ρ^,^) to the initial value problem (<ref>)-(<ref>)  is the purpose of various papers published in the late 90' <cit.>.
In this section, we wish to recall some of the key arguments that are used in these works to pass to the limit on the global weak solutions of (<ref>)-(<ref>) as the Mach number vanishes.

The results proven in the above-mentioned papers on the convergence of (ρ^,^) towards a weak solution of the incompressible Navier-Stokes equations follow a two-step argument. 
The first step consists in deriving a priori bounds on the quantities ρ^-1 and ^ which are uniform with respect to . These bounds imply the strong convergence of ρ^ towards ρ̅=1 in ^∞((0,T);^γ(Ω)), and up to the extraction of a subsequence, the weak convergence in ^2((0,T); ^1_0(Ω)) of ^ towards some function .  
The second step consists in passing to the limit in the weak formulation of problem (<ref>)-(<ref>) thanks to these convergence properties. 
The main difficulty in this step is the passage to the limit for the term (ρ^^⊗^) with only a weak convergence of the velocity. 

Subsequently, we describe the main arguments to obtain the estimates on ρ^-1 and ^, which we will mimick at the discrete level in order to prove the asymptotic preserving feature of the staggered schemes in the low Mach number limit.
However, we do not need a careful study of the nonlinear term since at the discrete level, i.e. for a fixed mesh, all norms are equivalent and the boundedness of (^)_>0 is enough to obtain convergence in any finite dimensional norm up to the extraction of a subsequence, and then to pass to the limit on the numerical scheme.





 §.§ A priori estimates


We begin by recalling some key identities satisfied by the smooth solutions of (<ref>), which are then incorporated in the definition of weak solutions.


Let ψ_γ be the function defined for ρ>0 as ψ_γ(ρ)=ρlogρ if γ=1, and ψ_γ(ρ)=ρ^γ/(γ-1) if γ>1 and define Π_γ(ρ)=ψ_γ(ρ)-ψ_γ(1)-ψ_γ'(1)(ρ-1). The smooth solutions of (<ref>)-(<ref>) satisfy the following identities:

 
  * A kinetic energy balance:
  
    _t (1/2ρ^  |^|^2) + (1/2ρ^  |^|^2 ^) 
      - ((^)) ·^
      + 1/^2 (ρ^)·^=0.

 
  * A renormalization identity
 
    _t ψ_γ(ρ^) + (ψ_γ(ρ^) ^) + (ρ^) ^ = 0.

 
  * A "positive" renormalization identity:
 
    _t Π_γ(ρ^) + (ψ_γ(ρ^) ^ - ψ_γ'(1)ρ^^) + (ρ^) ^ = 0.

  
  * An entropy identity:
    
    _t (1/2ρ^  |^|^2)  
        + 1/^2 _t Π_γ(ρ^) 
        +(  (1/2ρ^  |^|^2 
        + 1/^2 ψ_γ(ρ^) 
        
    - 1/^2ψ_γ'(1)ρ^ 
        + 1/^2(ρ^)  ) ^ ) 
        - ((^)) ·^ = 0.





The proof of (<ref>) and (<ref>) are classical. Multiplying the mass conservation equation (<ref>) by -ψ_γ'(1) and summing with (<ref>) yields (<ref>). 
Summing (<ref>) and ^-2×(<ref>) yields (<ref>).



We call (<ref>) a "positive" renormalization identity because the function Π_γ is positive thanks to the convexity of ψ_γ.


Integrating (<ref>) over Ω×(0,t) and recalling the homogeneous Dirichlet boundary conditions on the velocity, we obtain the following estimate on the solution pair (ρ^,^).
For all t∈ (0,T):

    1/2∫_Ωρ^(t)  |^(t)|^2 +  1/^2∫_ΩΠ_γ(ρ^(t)) + μ∫_0^t^(s)^2_^2(Ω)^d× d ds 
     + (μ+λ)  ∫_0^t   (^(s))_^2(Ω)^2  ds
      =  1/2∫_Ωρ_0^  |_0^|^2 +  1/^2∫_ΩΠ_γ(ρ_0^).






 §.§ Asymptotic behavior of the density and velocity in the zero Mach limit


The global entropy estimate (<ref>) is thus proven for any strong solution of the boundary and initial value problem (<ref>)-(<ref>). In the following, we always assume the existence of a weak solution (ρ^,^) to problem (<ref>)-(<ref>) satisfying the following inequality:

    1/2∫_Ωρ^(t)  |^(t)|^2 +  1/^2∫_ΩΠ_γ(ρ^(t)) + μ∫_0^t^(s)^2_^2(Ω)^d× d ds 
     ≤ 1/2∫_Ωρ_0^  |_0^|^2 +  1/^2∫_ΩΠ_γ(ρ_0^).

In particular, we assume that all the integrals involved in (<ref>) are convergent.

Estimate (<ref>) allows to prove the convergence of ρ^ towards 1 in ^∞((0,T);^γ(Ω)) provided that the right hand side, which depends on the initial conditions is uniformly bounded with respect to .
We show in the following that it is indeed the case provided some assumptions on the initial data.

We begin by proving the following lemma, which states crucial properties of the function Π_γ.


The function Π_γ has the following lower bounds:

 
    [                                                   ∙; [ Π_γ(ρ) ≥ C_γ,δ |ρ-1|^γ,     ∀ρ>0  with |ρ-1|≥δ, ] ]
       ∙  Π_γ(ρ) ≥ |ρ-1|^2,  ∀ρ>0. 
       [                                                                                                    ∙; [ Π_γ(ρ) ≥ C_γ,R |ρ-1|^2,              ∀ρ∈ (0,R),; Π_γ(ρ) ≥ C_γ,R |ρ-1|^γ,              ∀ρ∈ [R,∞). ] ]


Moreover, the function Π_γ has the following upper bound (for small densities): For all γ≥ 1 there exists C_γ such that:

    Π_γ(ρ) ≤ C_γ |ρ-1|^2, ∀ρ∈ (0,2).




For γ=1, we have Π_1=ρlogρ-ρ. Hence Π_1∼ρlogρ for large values of ρ, which implies (<ref>). Similarly, for γ>1, we have Π_γ(ρ)=ψ_γ(ρ)-ψ_γ(1)-ψ_γ'(1)(ρ-1)=(γ-1)^-1(ρ^γ-1-γ(ρ-1)), thus Π_γ(ρ)∼ (γ-1)^-1ρ^γ for large values of ρ which proves (<ref>).
A second order Taylor expansion of ψ_γ yields, for all γ≥ 1:
 
    Π_γ(ρ) =  |ρ-1|^2  γ∫_0^1 (1+s(ρ-1))^γ-2(1-s) ds,     for all  ρ∈(0,+∞).

The case γ≥ 2 is straightforward and we obtain |ρ-1|^2≤Π_γ(ρ) for all ρ∈(0,+∞) and Π_γ(ρ)≤ C_γ|ρ-1|^2 for all ρ≤ 2 with C_γ = γ∫_0^1 (1+s)^γ-2(1-s) ds. For 1≤γ < 2, we easily get  Π_γ(ρ)≤ |ρ-1|^2 for all ρ∈(0,+∞) and the lower bound is obtained by separating the case ρ < R and ρ≥ R. We obtain the expected lower bound (<ref>) with 

    C_γ,R=γ∫_0^1  1-s/(1+s(R-1))^2-γ ds.




	Lemma <ref> considers the isentropic case (ρ) = ρ^γ. 
	However, the results of sections <ref>-<ref> hold in a more general barotropic case.
	Indeed, let  be a continuous and derivable function. 
	Defining the functions ψ(ρ)=ρ∫(s)/s^2 ds and Π(ρ)=ψ(ρ)-ψ(1)-ψ'(1)(ρ-1), easy computations show that the renormalization identities (<ref>) and (<ref>) are still valid. Moreover,



	a straightforward calculation shows that 
	
    Π(ρ) = (ρ-1)^2 ∫_0^1 '(s(ρ-1)+1)/s(ρ-1)+1 (1-s)  ds.

	 Hence, if  is an non-decreasing C^1 function such that '(1) >0,  the above integral is positive and it is a continuous function of ρ; thus, under these assumptions on , there exists C_,R and C_,R∈_+ such that 
	
    C_,R (ρ-1)^2  ≤Π_γ(ρ) ≤C_,R (ρ-1)^2  | ρ | ≤ R.

Note that in the discrete setting considered below, this estimate is sufficient since the discrete density ρ is bounded uniformly with respect to the Mach number. 
Indeed, the conservative discretisation of the mass balance (see Section <ref>) yields that the discrete density ρ satisfies ∫_Ωρ(x,t)  dx = ∫_Ωρ_0(x)  dx, so that on a given mesh, ρ is bounded by 1/|K|∫_Ωρ_0(x)  dx  where |K| is the measure of the smallest cell. 

The results of sections <ref>-<ref> are thus still valid in the barotropic case, for a non-decreasing C^1 function   such that '(1) >0.


Let us now assume that the initial data is "ill-prepared", in the following sense: ρ_0^∈^∞(Ω) with  ρ_0^ >0 for a.e. ∈Ω, _0^∈^2(Ω)^d and there exists C independent of  such that: 

    _0^_^2(Ω)^d +  1/ ρ^_0- 1_^∞(Ω) ≤ C.

This bound implies that ρ_0^ tends to 1 in ^∞(Ω) when → 0; moreover, we suppose that _0^ weakly converges in ^2(Ω)^d towards a function _0∈^2(Ω)^d. The initial data is said to be "ill-prepared" since ρ^_0-1 behaves like  and not like ^2/γ (or ^2) as suggested by the momentum equation, and the initial velocity is not required to be close to a divergence free velocity.

Under assumption (<ref>), an easy consequence of the upper bound (<ref>) on Π_γ is that the right-hand side of the estimate (<ref>) is bounded independently of the Mach number .
We then obtain that for any weak solution (ρ^,^) to (<ref>)-(<ref>) that satisfies the global entropy estimate (<ref>), the velocity ^ is bounded in ^2((0,T);_0^1(Ω)^d) uniformly with respect to .

A further consequence of the lower bound (<ref>) on Π_γ is the convergence of ρ^ towards 1 as →0 in ^∞((0,T);^γ(Ω)) as stated in the following proposition.



Let, (ρ_0^,_0^) be a family of ill-prepared initial data and let (ρ^,^) be a corresponding family of weak solutions of (<ref>)-(<ref>) that satisfy the global entropy estimate (<ref>). Then, ρ^ converges towards 1 as →0 in ^∞((0,T);^γ(Ω)).



By (<ref>) and estimate (<ref>), we have for all δ>0, and t>0:
 
    ρ^(t)-1_^γ(Ω)^γ≤ |Ω| δ^γ + ∫_Ω |ρ^(t)-1|^γ _{|ρ^-1|≥δ}≤ |Ω| δ^γ + C ^2/C_γ,δ.

where, for a given set A, _A denotes the characteristic function of A.
Hence,

    lim sup_→ 0ρ^-1_^∞((0,T);^γ(Ω))≤ |Ω|^1/γ δ

for all δ>0, which concludes the proof.


The following proposition provides a rate of convergence in ^∞((0,T);^q(Ω)) of ρ^ towards 1 for q∈[1,min(2,γ)].

 
Let, (ρ_0^,_0^) be a family of ill-prepared initial data and let (ρ^,^) be a corresponding family of weak solutions of (<ref>)-(<ref>) that satisfy the global entropy estimate (<ref>).
Then, the following estimates hold.


  ∙ If γ≥ 2, then there exists C>0 such that, for  small enough:

    ρ^-1_^∞((0,T);^2(Ω))≤ C.


  ∙ If 1≤γ < 2, then for  small enough, for all R∈(2,+∞), there exists C_R>0 such that:

    ^-1(ρ^-1)_{ρ^<R }_^∞((0,T);^2(Ω))
    +  ^-2/γ(ρ^-1)_{ρ^≥ R }_^∞((0,T);^γ(Ω))≤ C_R.



 As a consequence, for all q∈[1,min(2,γ)], there exists C>0 such that for  small enough:
 
    ρ^-1_^∞((0,T);^q(Ω))≤ C.




As already stated, thanks to (<ref>) and using the upper bound on Π_γ(ρ) for small values of ρ, for  small enough, the right hand side of (<ref>) is bounded by some constant C_0, independent of .
We now use the lower bounds on Π_γ(ρ).
For γ≥ 2, by Lemma <ref> combined with estimate (<ref>), we have for all t∈(0,T):

    ρ^(t)-1_^2(Ω)^2 ≤∫_ΩΠ_γ(ρ^(t))≤ C_0 ^2.

For 1≤γ < 2, invoking once again Lemma <ref> and estimate (<ref>), we obtain for all t∈(0,T) and for all R∈(2,+∞):

    [                                                            (i)  (ρ^(t)-1)_{ρ^(t)≤ R}_^2(Ω)^2 ≤ 1/C_γ,R∫_ΩΠ_γ(ρ^(t))  ≤ C  ^2,;                                                           (ii) (ρ^(t)-1)_{ρ^(t)≥ R }_^γ(Ω)^γ≤ 1/C_γ,R∫_ΩΠ_γ(ρ^(t))  ≤  C  ^2, ]

which concludes the proof.






§ MESHES AND UNKNOWNS
 

 
A staggered discretization of Ω, denoted by , is given by a couple =(,), where:


  * , the primal mesh, is a finite family composed of non empty triangles and convex quadrilaterals for d=2 or non empty tetrahedra and convex hexahedra for d=3. The primal mesh  is assumed to form a partition of Ω : Ω= ∪_K ∈K.
For any K∈, let K  = K∖ K be the boundary of K, which is the union of cell faces.
We denote by  the set of faces of the mesh, and we suppose that two neighbouring cells share a whole face: for all ∈, either ⊂Ω or there exists (K,L)∈^2 with K ≠ L such that K∩L  =; we denote in the latter case = K|L.
We denote by  and  the set of external and internal faces: ={∈, ⊂Ω} and =∖.
For K ∈, (K) stands for the set of faces of K. The unit vector normal to ∈(K) outward K is denoted by _K,.
In the following, the notation |K| or || stands indifferently for the d-dimensional or the (d-1)-dimensional measure of the subset K of ^d or  of ^d-1 respectively.

  * We define a dual mesh associated with the faces ∈ as follows.
When K∈ is a simplex, a rectangle or a cuboid, for ∈(K), we define D_K, as the cone with basis  and with vertex the mass center of K (see Figure <ref>).
We thus obtain a partition of K in m sub-volumes, where m is the number of faces of K, each sub-volume having the same measure | D_K,|= |K|/m.
We extend this definition to general quadrangles and hexahedra, by supposing that we have built a partition still of equal-volume sub-cells, and with the same connectivities.
The volume D_K, is referred to as the half-diamond cell associated with K and .
For ∈, =K|L, we now define the diamond cell D_ associated with  by D_=D_K,∪ D_L,.
We denote by (D_) the set of faces of D_, and by =D_|D_' the face separating two diamond cells D_ and D_'.
As for the primal mesh, we denote by  the set of dual faces included in the domain and by  the set of dual faces lying on the boundary Ω.
In this latter case, there exists ∈ such that =.



Relying on this definition, we now define a staggered space discretization.
The degrees of freedom for the density ( the discrete density unknowns) are associated with the cells of the mesh , and are denoted by:

    {ρ_K, K ∈},

while the degrees of freedom for the velocity are located at the center of the faces of the mesh  and are therefore associated with the cells of the dual mesh D_, ∈ (as in the low-degree nonconforming finite-element discretizations proposed in <cit.>).
The Dirichlet boundary conditions are taken into account by setting the velocity unknowns associated with an external face to zero, so the set of discrete velocity unknowns reads:

    {_∈^d, ∈}.


We associate functions with the discrete unknowns of the schemes described hereinafter. To this purpose, we define the following sets of discrete functions of the space variable.


Let =(,) be a staggered discretization of Ω as defined in Definition <ref>.

 
  * We denote by _(Ω)⊂^∞(Ω) the space of scalar functions which are piecewise constant on each primal mesh cell K∈.
For all w∈_(Ω) and for all K∈, we denote by w_K the constant value of w in K, so the function w reads:

    w()= ∑_K ∈ w_K 𝒳_K()     ∈Ω,

where 𝒳_K stands for the characteristic function of K.


  * We denote by _(Ω)⊂^∞(Ω) the space of scalar functions which are piecewise constant on each diamond cell of the dual mesh D_, ∈.
For all u ∈_(Ω) and for all ∈, we denote by u_ the constant value of u in D_, so the function u reads:

    u()= ∑_∈ u_ 𝒳_D_()     ∈Ω,

where 𝒳_D_() stands for the characteristic function of D_. We denote by _(Ω)=_(Ω)^d the space of vector valued (in ^d) functions that are constant on each diamond cell D_. Finally,
we denote _,0(Ω)={ u∈_(Ω), u_=0  for all ∈} and _,0(Ω)=_,0(Ω)^d .










§ SPACE DISCRETIZATION
  

This section is devoted to the construction of the discrete space differential operators that approximate the differential operators in (<ref>).
As already said, the discretization is staggered.
The discrete operators involved in the discretization of the mass equation (<ref>) are thus associated with the cells of the primal mesh K, K ∈, while the discrete operators involved in the discretization of the momentum equation (<ref>) are associated with the cells of the dual mesh D_, ∈.





 §.§ Mass convection flux
 
The discretization of the convection term (ρ) in the mass conservation equation is defined as follows. Given a discrete density field ρ∈_(Ω) and a velocity field ∈_,0(Ω), it is a piecewise constant function on each primal cell K∈ given by:

    (ρ)_K = 1/|K|∑_∈(K) F_K,(ρ,),     ∀ K ∈.

The quantity F_K,(ρ,) stands for the mass flux across  outward K.
By the impermeability boundary conditions, it vanishes on external faces and is given on internal faces by:

    F_K,(ρ,)= || ρ_ _·_K,,     ∀∈, =K|L.

The density at the face =K|L is approximated by the upwind technique,  ρ_=ρ_K if _·_K,≥ 0 and ρ_=ρ_L otherwise.





 §.§ Velocity convection operator
 

We now describe the approximation of the convection operator ∂_t(ρ) + (ρ⊗) appearing in the momentum balance equation.
The approximation of the time derivative part ∂_t(ρ) is naturally discretized at the dual cells D_, ∈ and an approximation ρ_ of the density on these dual cells D_ is thus needed.
Given a density field ρ∈_(Ω), this approximation is built as follows:

    |D_| ρ_ = |D_K,| ρ_K + |D_L,| ρ_L,      ∀∈, =K|L.


In order to consider different time discretizations, we build a space discretization of a more general divergence part of the velocity convection operator, (ρ⊗).
Given a discrete density field ρ∈_(Ω), and two discrete velocity fields ∈_,0(Ω) and ∈_,0(Ω), this approximation is built as follows:

    (ρ⊗)_= 1/|D_|∑_∈(D_)(ρ,) _,     ∀∈.

(ρ,) is the mass flux across the edge  of the dual cell D_. Its value is zero if ∈. Otherwise, it is defined as a linear combination, with constant coefficients, of the primal mass fluxes at the neighboring faces. 
For K ∈ and ∈(K), let ξ_K^ be given by:

    ξ_K^=|D_K,|/|K|,

so that ∑_∈(K)ξ_K^=1.
With the definition of the dual mesh adopted here, the value of the coefficients ξ_K^ only depend on the type of the cell K (simplicial or quandrangular/hexahedral).
For the quadrangular and hexahedral elements, we have ξ_K^=1/(2d) and, for the simplicial elements, ξ_K^=1/(d+1).
Then the mass fluxes through the inner dual faces are supposed to satisfy the following properties.


  (H1) The discrete mass balance over the half-diamond cells is satisfied, in the following sense.
For all primal cell K in , the set ((ρ,))_⊂ K of dual fluxes included in K solves the following linear system

    F_K,(ρ,) + ∑_∈(D_σ), ⊂ K F_,(ρ,)=
    ξ_K^∑_' ∈(K) F_K,'(ρ,),   ∈(K).


  (H2) The dual fluxes are conservative,  for any dual face =D_|D_', we have F_,(ρ,)=-F_',(ρ,).

  (H3) The dual fluxes are bounded with respect to the primal fluxes (F_K,(ρ,))_∈(K), in the sense that

    |F_,ϵ(ρ,)| ≤ max { |F_K,'(ρ,)|, ' ∈(K) },

for K ∈, ∈(K), ϵ∈(D_σ) with ⊂ K.

The system of equations (<ref>) only depends on the type of the cell K (since it only depends on the coefficient ξ_K^, which is chosen as the inverse of the number of cell faces, and sub-cell connectivities) but has an infinite number of solutions, which makes necessary to impose in addition the constraint (<ref>); however, assumptions (H1)-(H3) are sufficient for the subsequent developments, in the sense that any choice for the expression of the fluxes satisfying these assumptions yields stable and consistent schemes (see <cit.>).

To complete the definition of the convective flux, we just have now to give the expression of the velocity _ at the dual face. As already said, a dual face lying on the boundary is also a primal face, and the flux across that face is zero.
Therefore, the values _ are only needed at the internal dual faces; we choose them to be centered:

    _ = 1/2 ( _ + _'),     =D_|D_'.






 §.§ Diffusion term


The space discretization of the diffusion term (()) in the momentum balance equation relies on the Crouzeix-Raviart element for the simplicial cells K and on the parametric Rannacher-Turek (or rotated bilinear) element for quadrangular or hexahedral cells  (see <cit.>).
Let ℙ_K be the affine transformation between the reference unit simplex and the simplicial cell K.
The space of discrete functions over a simplicial cell K is

    P_1(K)= { f∘ℙ_K^-1, with  f∈ span { 1, (x_i)_i=1,…,d}}.

Let ℚ_K be the standard Q_1 mapping between the reference unit cuboid and the quadrangular or hexahedral cell K.
The space of discrete functions over such a cell is 

    Q̃_1(K)= { f∘ℚ_K^-1, with f∈ span { 1, (x_i)_i=1,…,d, ( x_i^2- x_i+1^2)_i=1,…,d-1}}.

The shape functions are the functions ζ_, ∈ such that for all K∈, ζ_|_K∈ P_1(K) if K is a simplicial cell and ζ_|_K∈Q̃_1(K) if K is a quadrangular or hexahedral cell and which satisfy the following two conditions: 


    (i)      ∫_ [ζ]_() = 0,   where [ζ]_() = lim_→
    ∈ Lζ()-
    lim_→
    ∈ Kζ(),
     ∀∈,  ∀∈,  =K|L. 
     (ii)      1|'|∫_'ζ_() = δ_^',
      ∀, ' ∈,


with δ_^'=1 if =' and δ_^'=0 otherwise.
Note that condition (i) is consistent with a location of the velocity degrees of freedom at the faces.

Now, with a discrete velocity field ∈_,0(Ω), one classically associates, in the finite element context, the function ()=∑_∈_ζ_(). 
The discretization of the diffusion term is a piecewise constant function on each diamond cell D_, the value of which reads

    (())_ = 
     - μ ( 1/|D_|∑_K ∈∫_K   . ζ_ ) 
    - (μ + λ)  (1/|D_|∑_K ∈∫_K () ζ_ ).

The identification between ∈_,0(Ω) and  allows to introduce the broken Sobolev ^1 semi-norm ._, given for any ∈_(Ω) by:

    _^2=∑_K∈∫_K :.

The semi-norm _ is in fact a norm on the space _,0(Ω), thanks to a classical discrete Poincaré inequality.

As in the continuous setting, it is easily seen that the bilinear form derived from the discretization of the diffusion term controls the discrete ^1-norm of the velocity as stated in the following lemma:
 
For every discrete velocity field ∈_,0(Ω), one has:

    ∑_∈ |D_| _·( -  (())_ ) ≥μ _^2.







 §.§ Pressure gradient term


The discretization of the pressure gradient term (ρ) is a piecewise constant function on each diamond cell D_, the value of which is denoted ( p)_.
For ρ∈_(Ω), this term is defined as:

    ( p)_=||/|D_| ((ρ_L)-(ρ_K)) _K,,     ∀=K|L ∈.

This pressure gradient is only defined at internal faces since, thanks to the impermeability boundary conditions, no momentum balance equation is written at the external faces. 

The following discrete duality relation holds for all ρ∈_(Ω) and ∈_,0(Ω):

    ∑_K ∈ |K| (ρ_K)  ()_K
    +∑_∈ |D_| _· ( p)_
    =0,

where we have set for all K∈, ()_K = |K|^-1∑_∈(K)|| _·_K, (consistently with (<ref>) for ρ≡ 1).

We finish this section with the following lemma which states that the staggered approximation is inf-sup stable.
 
There exists β>0, depending only on Ω and on the mesh, such that for all p={ p_K, K∈}∈_(Ω), there exists ∈_,0 satisfying:

    _=1  and ∑_K∈ |K| p_K  ( )_K ≥β p-m(p)_L^2(Ω),

where m(p)=|Ω|^-1∑_K∈|K|p_K is the mean value of p over Ω.


The inf-sup property is crucial when passing to the limit → 0 in the various schemes presented thereafter.
It indeed provides an ^2 control on the discrete (zero-mean) pressure through the control of its gradient.

In fact, the actual inf-sup stability condition states that the constant β only depends on the regularity of the mesh (in a sense to be defined), and not on the space step; this property, which is satisfied by low-order staggered discretizations, is inherited by the limit incompressible scheme and guarantees its stability and the fact that error estimates do not blow up when the mesh is refined.
In this paper, since we work on a fixed discretization, the dependency of β with respect to the mesh does not need to be precisely stated.    





§ ASYMPTOTIC ANALYSIS OF THE ZERO MACH LIMIT FOR AN IMPLICIT SCHEME


We begin with the analysis of the zero Mach limit for a fully implicit scheme.
Let δ t>0 be a constant time step. 
The approximate solution (ρ^n,^n)∈_(Ω)×_,0(Ω) at time t_n=nδ t for 1≤ n≤ N=T/δ t is computed by induction through the following implicit scheme. 

Knowing (ρ^n,^n)∈_(Ω)×_,0(Ω), solve for ρ^n+1∈_(Ω) and ^n+1∈_,0(Ω):


    1 δ t(ρ^n+1_K-ρ^n_K) + (ρ^n+1^n+1)_K = 0,    ∀ K ∈,
    
        1 δ t(ρ^n+1_^n+1_-ρ^n__^n )
    + (ρ^n+1^n+1⊗^n+1)_
        -  ((^n+1))_
    + 1^2( p^n+1)_
    =0,    ∀∈.







 §.§ Initialization of the scheme


The initial approximations are given by the average of the initial density ρ_0^ on the primal cells and the initial velocity _0^ on the dual cells:

    [ ρ_K^0 = 1/|K|∫_K ρ_0^,                ∀  K ∈,;   _^0 = 1/|D_|∫_D__0^,                   ∀ ∈. ]


Thereafter, we prove that for every >0, there exists a solution (ρ^,^) to the implicit scheme (<ref>)-(<ref>) and that for a fixed discretization, i.e. for a fixed mesh and a fixed time step δ t, the solution  (ρ^,^) converges as → 0 towards the solution of an implicit scheme for the incompressible Navier-Stokes equations, which is stable thanks to the inf-sup condition.

Assumption on the initial data – For the convergence study performed in this section, it is sufficient to assume that the initial data is ill-prepared, in the sense of Inequality (<ref>).





 §.§ A priori estimates
 

We begin with a first lemma which states that the velocity convection operator defined in Section <ref> is built so that if a discrete mass conservation equation is satisfied on the cells of the primal mesh (as in (<ref>)) - which is consistent with the staggered discretization - then a discrete mass conservation equation is also satisfied on each cell of the dual mesh (see <cit.> for a proof).



Let two density fields ρ^n, ρ^n+1∈_(Ω) and a velocity field ^n+1∈_,0(Ω) satisfying the discrete mass conservation equation (<ref>) on every cell of the primal mesh be given.
Then, the dual densities {ρ_^n, ρ_^n+1,  ∈} and the dual fluxes {(ρ^n+1,^n+1),  ∈,  ∈(D_) } satisfy a finite volume discretization of the mass balance (<ref>) over the internal dual cells:

    |D_|/δ t (ρ^n+1_-ρ^n_)
    + ∑_∈(D_)(ρ^n+1,^n+1)=0,     ∀∈.



The following result states that any solution of the implicit scheme satisfies a discrete counterpart to the kinetic energy balance (<ref>); its derivation relies on the previous relation, namely the dual mass balance (<ref>).

 
Any solution to the implicit scheme (<ref>) satisfies the following equality, for all ∈, and 0≤ n ≤ N-1:

    1 2δ t ( ρ_^n+1 |_^n+1|^2 - ρ_^n |_^n|^2  )
    + 1/2|D_|∑_=D_|D_'(ρ^n+1,^n+1)  _^n+1·_'^n+1
    
    - ((^n+1))_·_^n+1 
    + 1/^2( p^n+1)_·_^n+1 + R_^n+1=0,

where R_^n+1= 12δ tρ_^n|_^n+1-_^n|^2.



Let us take the scalar product of the discrete momentum balance equation (<ref>) by the corresponding velocity unknown _^n+1, which gives the relation T_^ conv - ((^n+1))_·_^n+1 + 1/^2( p^n+1)_·_^n+1=0, with:

    T_^ conv =
    (
     1 δ t(ρ^n+1_^n+1_-ρ^n__^n )
    +  1 2|D_| ∑_∈(D_) 
    =(ρ^n+1,^n+1) (^n+1_+^n+1_')
    )·_^n+1.

Now, using the identity 2  (ρ||^2-ρ^*·) = ρ||^2 - ρ^*||^2 + ρ^* |-|^2 + (ρ-ρ^*)||^2 with ρ=ρ_^n+1, ρ^*=ρ_^n, =_^n+1 and =_^n, we obtain

    T_^ conv   
    =  1 2δ t( ρ^n+1_ |^n+1_|^2-ρ^n_ |_^n|^2 )
    +  1 2 |D_| ∑_=(ρ^n+1,^n+1) ^n+1_·^n+1_'
       
    +  1 2δ t ρ^n_ | ^n+1_ - ^n_ |^2
    + (1/δ t (ρ^n+1_-ρ^n_)
    + |D_|^-1∑_∈(D_)(ρ^n+1,^n+1) ) |_^n+1|^2/2.

The last term is equal to zero thanks to (<ref>), which concludes the proof.


We then prove that any solution of the implicit scheme satisfies a discrete counterpart of the renormalization identities (<ref>) and (<ref>) satisfied by any smooth solution of (<ref>).
To state this result, we need to extend the notation for divergence operators on the primal cells as follows.
For K ∈ and a smooth function φ,

    (φ(ρ) )_K = 1/|K|∑_∈(K) || φ(ρ_) _·_K,,

where ρ_ stands for the upwind value of the density at the face.

 
Define the function ψ_γ as 
ψ_γ(ρ)= ρ∫^ρ_0(s)/s^2 ds, which yields ψ_γ(ρ)=ρlogρ if γ=1 and ψ_γ(ρ)=ρ^γ/(γ-1) if γ>1, and define Π_γ(ρ)=ψ_γ(ρ)-ψ_γ(1)-ψ_γ'(1)(ρ-1). Then, any solution to the implicit scheme (<ref>) satisfies the following two identities, for all K∈ and 0≤ n ≤ N-1:

    1 δ t ( ψ_γ(ρ^n+1_K)-ψ_γ(ρ^n_K)  )
    + (ψ_γ(ρ^n+1) ^n+1)_K
    + (ρ_K^n+1)  (^n+1)_K + R_K^n+1 = 0,
    
        1 δ t ( Π_γ(ρ^n+1_K)-Π_γ(ρ^n_K)  ) 
    + ((ψ_γ(ρ^n+1) - ψ_γ'(1) ρ^n+1) ^n+1)_K
    
       
    +  (ρ_K^n+1)  ( ^n+1)_K + R_K^n+1 = 0,

with :

    R_K^n+1 = 1/2 δ t ψ_γ”(ρ̅_K^n+1/2)  (ρ_K^n+1-ρ_K^n)^2
    +1/2 |K|∑_ = K|L ||   (_^n+1·_K,)^-  ψ_γ”(ρ̅_^n+1)   (ρ_L^n+1-ρ_K^n+1)^2,

where ρ̅_K^n+1/2∈[min(ρ_K^n+1,ρ_K^n),max(ρ_K^n+1,ρ_K^n)], ρ̅_^n+1∈[min(ρ_^n+1,ρ_K^n+1),max(ρ_^n+1,ρ_K^n+1)] for all ∈(K), and for a∈, a^-≥ 0 is defined by a^-=-min(a,0). Since ψ_γ is a convex function, R_K^n+1 is non-negative.



The proof of (<ref>) is the same as that of <cit.>.
The equality  (<ref>) is obtained after multiplying the discrete mass equation (<ref>) by -ψ_γ'(1) and summing with (<ref>).


As a consequence of Lemmas <ref> and <ref>, the implicit scheme satisfies a discrete counterpart of the local-in-time entropy identity (<ref>).

 
Let >0 and assume that the initial density ρ_0^ is positive. Then, there exists a solution (ρ^n,^n)_0≤ n≤ N to the scheme (<ref>), such that ρ^n >0 for 0≤ n≤ N, and the following inequality holds for 0≤ n≤ N-1:

    1/2∑_∈
    |D_| ( ρ^n+1_ |^n+1_|^2-ρ^n_ |_^n|^2 )
    + 1/^2∑_K ∈ |K|  (Π_γ(ρ^n+1_K) - Π_γ(ρ^n_K))
     +  μ δ t ^n+1_1,^2
     + ℛ^n+1≤ 0,

where ℛ^n+1 = ∑_∈ R^n+1_ + ^-2∑_K ∈ R_K^n+1≥ 0.



The positivity of the density is a consequence of the properties of the upwind choice (<ref>) for ρ <cit.>.
Let us then sum equation (<ref>) over the faces ∈, ^-2× (<ref>) over K ∈, and, finally, the two obtained relations.
Since the discrete gradient and divergence operators are dual with respect to the ^2 inner product (see (<ref>)), noting that the conservative dual fluxes vanish in the summation and that the diffusion term is coercive (see Lemma <ref>), we get (<ref>).

Given discrete density and velocity fields (ρ^n,^n), the existence of a solution (ρ^n+1,^n+1) to the implicit scheme at the time step n+1 may be inferred by the Brouwer fixed point theorem, by an easy adaptation of the proof of <cit.>.
This proof relies on the following set of mesh-dependent estimates: the conservativity of the mass balance discretization, together with the fact that the density is positive, yields an estimate for ρ in the ^1-norm, and so, by a norm equivalence argument, of the pressure in any norm; then, for a given density ρ, the discrete global kinetic energy inequality ( (<ref>) summed over the control volumes) provides a control on the velocity.
Therefore, computing ρ from the mass balance for fixed , then p from ρ by the equation of state (ρ), and finally  from the momentum balance equation with fixed ρ, yields an iteration in a bounded convex subset of a finite dimensional space.


We may now prove that the solution of the implicit scheme satisfies a discrete counterpart of the global estimate (<ref>).


Let >0 and assume that the initial data (ρ_0^,_0^) is ill-prepared in the sense of (<ref>). By Lemma <ref>, there exists a solution (ρ^n,^n)_0≤ n≤ N to the scheme (<ref>). In addition, there exists C>0 independent of  such that, for  small enough and for all 1≤ n≤ N:

    1/2∑_∈ |D_|  ρ^n_ |^n_|^2
    +μ ∑_k=1^n δ t ^k_1,^2
    +1/^2∑_K∈|K|  Π_γ(ρ^n_K) ≤ C.




Multiplying equation (<ref>) by δ t and summing over the time steps yields for 1≤ n≤ N:

    1/2∑_∈ |D_|  ρ^n_ |^n_|^2
    +μ ∑_k=1^n δ t ^k_1,^2
    +1/^2∑_K∈|K|  Π_γ(ρ^n_K)+ℛ^n
    
    ≤ 1/2∑_∈ |D_|  ρ^0_ |^0_|^2
    +1/^2∑_K∈|K|  Π_γ(ρ^0_K),

with ℛ^n= ∑_k=0^n-1 (∑_∈ R^k+1_ + ^-2∑_K ∈ R_K^k+1 ) ≥ 0.
Let us prove that the right hand side of (<ref>) is uniformly bounded for all  small enough. By (<ref>) and (<ref>), for  small enough, one has ρ_K^0 ≤ 2 for all K∈ and therefore  ρ^0_≤ 2 for all ∈, since the dual densities are convex combinations of the primal density unknowns. Hence, one has:

    1/2∑_∈ |D_|  ρ^0_ |^0_|^2 ≤∑_∈ |D_|^-1 | ∫_D__0^ |^2 ≤_0^_^2(Ω)^d^2,

which by (<ref>) is uniformly bounded with respect to . Then, using the upper bound on Π_γ(ρ) for small values of ρ (see Lemma <ref>), we can see that the second term of the right hand side of (<ref>) is uniformly bounded with respect to , for  small enough.


 
Let >0 and assume that the initial density ρ_0^ is positive and that the initial data (ρ_0^,_0^) is ill-prepared in the sense of (<ref>). Then, there exists a solution (ρ^n,^n)_0≤ n≤ N to the scheme (<ref>). Let  p^n=(ρ^n) and define ^n = {_K^n, K∈} where _K^n = (p_K^n-m(p^n))/^2 with m(p^n) the mean value of p^n over Ω (i.e. m(p^n) = |Ω|^-1∑_K∈ |K| p_K^n). Then, one has, for all 1≤ n ≤ N:

    ^n≤ C_,δ t,

where  the real number C_,δ t depends on the mesh and the time step but not on , and · stands for any norm on the space of discrete functions.



The staggered discretization satisfies the inf-sup condition (see Lemma <ref>), which implies that there exists a positive real number β, depending only on Ω and on the mesh, such that for ^n= (p^n-m(p^n))/^2, there exists a discrete velocity field , with _=1 and satisfying:

    β ^n_L^2(Ω)≤∑_K∈ |K| p_K^n  ( )_K.

Hence by the gradient-divergence duality property (<ref>), taking the scalar product of (<ref>) with |D_| _ and summing over  in , we get β^n_L^2(Ω)≤ T^n_1+T^n_2+T^n_3 with:

    T^n_1 = 1 δ t ∑_∈ |D_| (ρ^n_^n_-ρ^n-1__^n-1)·_,
    
       
    T^n_2 =∑_∈ |D_|  (ρ^n ^n ⊗^n)_·_,
    
       
    T^n_3 =- ∑_∈ |D_|  ((^n))_·_.

The proof of (<ref>) follows if one is able to control each of these terms T^n_1, T^n_2 and T^n_3 independently of . Since the mesh and the time step δ t are fixed, a norm equivalence argument in a finite dimensional space yields the existence of a positive function C_,δ t depending on the discretization, which is non-decreasing in each of its variables, such that: 

    |T^n_1+T^n_2+T^n_3| ≤ C_,δ t (ρ^n_^1,ρ^n-1_^1,^n_, ^n-1_, _).

We have  _=1 and by the estimate (<ref>), ^n_ and ^n-1_ are controlled independently of . In addition ρ^n and ρ^n-1 are controlled in ^1 by conservativity of the discrete mass balance equation.






 §.§ Incompressible limit of the implicit scheme


We may now state the main result of this section which is the convergence, up to a subsequence, of the solution to the compressible scheme (<ref>) towards the solution of an implicit inf-sup stable incompressible scheme when the Mach number tends to zero.



Let ()_m∈ be a sequence of positive real numbers tending to zero, and let (ρ,)_m ∈ be a corresponding sequence of solutions of the scheme (<ref>).
Let us assume that the initial data (ρ_0^,_0^) is ill-prepared,  satisfies Relation (<ref>) for m ∈.
Then the sequence (ρ)_m ∈ tends to the constant function ρ =1 when m tends to + ∞ in ^∞((0,T),^γ(Ω)). Moreover, for all q∈[1,min(2,γ)], there exists C>0 such that:

    ρ-1_^∞((0,T);^q(Ω))≤ C,     for m large enough.


In addition, the sequences ()_m ∈ and ()_m ∈ are bounded in any discrete norm which may depend on the fixed discretization.  If a subsequence of (,)_m ∈ tends, in any discrete norm, to a limit (,), then (,) is a solution to the standard (Rannacher-Turek or Crouzeix-Raviart) implicit scheme for the incompressible Navier-Stokes equations: 

Knowing  ^n∈_(Ω) and ^n∈_,0(Ω), solve for ^n+1∈_(Ω) and ^n+1∈_,0(Ω):


    (^n+1)_K = 0,    ∀ K ∈,
    
        1 δ t( ^n+1_ - _^n )
    + ( ^n+1⊗^n+1)_
    - ((^n+1))_
    +(^n+1)_
    =0,    ∀∈.


At the limit, the scheme uses as initial condition only the ^2-projection of _0 (the limit of _0^) on the space E_(Ω) of the discrete divergence-free functions:
E_(Ω)={∈_,0(Ω), ()_K=0, ∀ K∈}.



The proof of the convergence of (ρ)_m ∈ towards ρ =1 when m tends to + ∞ is the same as in the continuous setting. It follows from the combination of the lower bounds on Π_γ(ρ) proven in Lemma <ref> and the discrete global entropy estimate (<ref>).

Using again estimate (<ref>), one can see that the sequence ()_m ∈ is bounded in any discrete norm and the same holds for the sequence ()_m ∈ by Lemma <ref>.
By the Bolzano-Weierstrass Theorem and a norm equivalence argument for finite dimensional spaces, there exists a subsequence of  (,)_m ∈ which tends, in any discrete norm, to a limit ( ,).
Passing to the limit cell-by-cell in the implicit scheme (<ref>), one obtains that (, ) is a solution of the standard implicit scheme (<ref>) for the incompressible Navier-Stokes equations.

If ^1 is a solution of the first time-step of (<ref>) then ^1∈ E_(Ω), and for all ∈ E_(Ω), taking the scalar product of (<ref>) with |D_|_ and summing over ∈, one gets that ^1 satisfies:

    1/δ t ( ^1, )_^2 
    + ∑_∈ |D_| (ρ^1^1⊗^1)_·_
    + ∑_K∈∫_K (^1):()= 1/δ t ( ^0, )_^2.

This system is known to be the part of the algebraic system associated with one time step of the scheme determining the velocity (even if the uniqueness of this unknown is guaranteed only for small time steps), in the sense that the velocity may be computed from (<ref>), the remaining equations yielding the pressure.
Observing that in the right hand side of (<ref>), ^0 can be replaced by its ^2-projection 𝒫(^0) onto E_(Ω) (since (^0-𝒫(^0), )_^2=0 for all ∈ E_(Ω) by definition of 𝒫(^0)), we obtain that only the divergence-free part of the initial velocity is seen by the implicit scheme.




Estimate (<ref>) provides a control on the sequence of discrete velocities ()_m∈ in a discrete ^2((0,T);^1_0(Ω)^d)-norm provided that μ>0. However, even for the Euler case where μ=0, one can derive a uniform bound on ()_m∈ in any discrete norm. Indeed, since the mesh is fixed, ρ→ 1 as m→+∞ means that the discrete density unknowns tend to 1 in every cell. Hence, by the kinetic energy part of (<ref>), there exists a minimum density ρ_ min,>0 depending on the fixed mesh such that _^2(Ω)^d≤√(2C/ρ_ min,) for all m∈. Therefore, the result of Theorem <ref> is still valid for the Euler equations: the solution to the compressible scheme (<ref>) (with μ=λ=0 and homogeneous Neumann boundary conditions on the velocity) converges, when the Mach number tends to zero, towards the solution of an implicit inf-sup stable scheme for the incompressible Euler equations.
 




§ ASYMPTOTIC ANALYSIS OF THE ZERO MACH LIMIT FOR A PRESSURE CORRECTION SCHEME
 

Since the scheme (<ref>) is fully implicit, the implementation of the algorithm implies to find the solution of a fully non-linear coupled system which is difficult in a real computational context due to the computational cost and lack of robustness.
This is why we also perform the analysis of the low Mach number limit for a semi-implicit scheme which is implemented in the software (CALIF^3S <cit.>). This scheme is obtained thanks to a partial decoupling of the discrete equations, and falls in the family of pressure correction schemes.
It consists (after a rescaling step for the pressure gradient) in two main steps. A prediction step, where a tentative velocity field is obtained by solving a linearized momentum balance in which the mass convection flux and the pressure gradient are explicit. A correction step where a nonlinear problem on the pressure is solved, and the velocity is updated in such a way to recover the discrete mass conservation equation. 

Let δ t>0 be a constant time step. 
The approximate solution at time t_n=nδ t for 1≤ n≤ N=T/δ t is denoted (ρ^n,^n)∈_(Ω)×_,0(Ω).

Knowing (ρ^n-1,ρ^n,^n)∈_(Ω)×_(Ω)×_,0(Ω), the considered algorithm consists in computing ρ^n+1∈_(Ω) and ^n+1∈_,0(Ω) through the following steps:


    ( p)^n_ = (ρ^n_/ρ^n-1_)^1/2 ( p^n)_,    ∀∈.
    
       
    
    
          1 δ t (ρ^n_^n+1_-ρ^n-1__^n )
    + (ρ^n  ^n ⊗^n+1 )_
    -  ((^n+1))_
    + 1^2 ( p)^n_
    =0,    ∀∈.
    
    
    
       
          1 δ t ρ^n_ (^n+1_-_^n+1)
    + 1^2 ( p^n+1)_ -  1 ^2 ( p)_^n
    =0,    ∀∈.
    
    
          1 δ t(ρ^n+1_K-ρ^n_K) + (ρ^n+1^n+1)_K = 0,    ∀ K ∈.







 §.§ Initial data and initialization of the scheme


As mentioned in the introduction, for the pressure correction scheme, it is not sufficient to assume "ill-prepared" initial data. In this section, the initial data (ρ_0^,_0^) are assumed to be well-prepared in the sense of the following definition:



The initial data (ρ_0^,_0^) is said to be well-prepared if ρ_0^ >0, ρ_0^∈^∞(Ω), _0^∈_0^1(Ω)^d for all >0 and if there exists C independent of  such that: 

    _0^_^1(Ω)^d + 1/  _0^_^2(Ω) + 1/^2 ρ^_0- 1_^∞(Ω)≤  C.

Consequently, ρ_0^ tends to 1 when → 0; moreover, we suppose that _0^ converges in ^2(Ω)^d towards a function _0∈^2(Ω)^d (the uniform boundedness of the sequence in the ^1(Ω)^d norm already implies this convergence up to a subsequence).


The initial velocities _0^ are in _0^1(Ω)^d. In particular, their trace is thus well defined as an ^2 function on any smooth (or Lipschitz-continuous) hypersurface of Ω. The initialization of the pressure correction scheme (<ref>) is performed as follows. First, ρ^0 and ^0 are given by the average of the initial conditions ρ_0^ and _0^ respectively on the primal cells and on the faces of the primal cells:

    [ ρ_K^0 = 1/|K|∫_K ρ_0^,                ∀  K ∈,;       _^0 = 1/||∫__0^,                   ∀ ∈. ]

Finally, we compute ρ^-1 by solving the mass balance equation (<ref>) for n=-1, where the unknown is ρ^-1 and not ρ^0. 
This procedure allows to perform the first prediction step with (ρ^-1_)_∈, (ρ^0_)_∈ and the dual mass fluxes satisfying the mass balance :
 
    |D_|/δ t (ρ^0_-ρ^-1_)
    + ∑_∈(D_)(ρ^0,^0)=0,     ∀∈.


In this section, we prove that for every >0, there exists a solution (ρ^,^) to the pressure correction scheme (<ref>) with the aforementioned initialization, and that for a fixed discretization, i.e. for a fixed mesh and a fixed time step δ t, the solution  (ρ^,^) converges as → 0 towards the solution of a pressure correction scheme for the incompressible Navier-Stokes equations, which is inf-sup stable. 

A first property must be checked so that the scheme is well defined, which is the positivity of the density at the fictitious time step n=-1. Indeed,  since ρ_0^ is assumed to be positive on Ω and by definition of the discrete density ρ^0 (<ref>), one clearly has ρ_K^0 > 0 for all K∈. This positivity property is not clear for ρ^-1. It is a consequence of the following result.

 
If the initial conditions (ρ_0^,_0^) are assumed to be well-prepared in the sense of Definition <ref>, then there exists a constant C independent of  such that:

    1/^2 max_K∈|ρ_K^0-1|  +  1/^2 max_∈ |( p^0)_|  + 1/ max_K∈|ρ_K^-1-1| ≤ C.




Since the initial conditions (ρ_0^,_0^) are well-prepared, the discrete density at the time step n=0 satisfies for all K∈:

    |ρ_K^0-1| ≤1/|K|∫_K |ρ_0^-1| ≤ ρ_0^-1_^∞(K)≤ C ^2.

Then, using the definition of the pressure gradient, we also obtain |( p^0)_| ≤ C ^2 for all ∈, where the constant C may depend on the (fixed) mesh and on =γ (or on '(1) in the general barotropic case) but is independent of . Let us then verify that for  small enough the discrete density at time step n=-1 is close to 1. This density is computed in such a manner that, for all K∈:

    ρ_K^-1 -1    = ρ_K^0 -1 + δ t ∑_∈(K)||/|K| ρ_^0  _^0·_K,
        = ρ_K^0 -1 + δ t ρ_K^0 ∑_∈(K)||/|K|  _^0·_K,  + δ t ∑_∈(K)||/|K|  (ρ_^0-ρ_K^0)  _^0·_K,
        = ρ_K^0 -1 + δ t ρ_K^0 1/|K|∑_∈(K)∫__0^·_K,
                      + δ t ∑_∈(K)||/|K|  (ρ_^0-ρ_K^0)  _^0·_K,
        = ρ_K^0 -1 + δ t ρ_K^0 1/|K| ∫_K  _0^  + δ t ∑_∈(K)||/|K|  (ρ_^0-ρ_K^0)  _^0·_K,.

We have already proven that the first term and thus the third term (using, for this latter, a trace inequality to bound the face velocity) in the right-hand side of the above equality are bounded by C ^2 for some constant C independent of  (but dependent on the mesh).
The second term is bounded by C for some constant C since the divergence of the initial velocity satisfies (<ref>).
The expected bound on |ρ_K^-1-1| follows.



All the developments below are still valid under the slightly less restrictive condition on the divergence of the initial velocities:  

    _0^_^2(Ω)→ 0   as →0.

Indeed, under this relaxed condition, we still have max_K∈|ρ_K^-1-1|→ 0 as → 0.






 §.§ A priori estimates


Let us now derive the estimates satisfied by the solutions of the pressure correction scheme.
As for the implicit scheme, any solution of the pressure correction scheme satisfies a discrete counterpart to the kinetic energy balance (<ref>).

 
Any solution to the pressure correction scheme (<ref>) satisfies the following equality, for all ∈ and 0≤ n ≤ N-1:

    1 2δ t ( ρ_^n |_^n+1|^2 - ρ_^n-1 |_^n|^2  ) + 1/2|D_|∑_=D_|D_'(ρ^n,^n)  _^n+1·_'^n+1
    - ((^n+1))_·_^n+1
    
    + 1/^2( p^n+1)_·_^n+1  
    + δ t/2 ^4 (  | ( p^n+1)_ |^2/ρ^n_   - | ( p^n)_ |^2/ρ^n-1_)+ R_^n+1=0,

where R_^n+1= 12δ tρ_^n-1|_^n+1-_^n|^2.



Let us take the scalar product of the velocity prediction equation (<ref>) with the corresponding velocity unknown _^n+1.
Thanks to the time shift in the dual densities and the dual mass fluxes in the velocity prediction equation (<ref>), together with the dual mass balance (at the previous time step) 

    1/δ t (ρ^n_-ρ^n-1_)+ |D_|^-1∑_∈(D_)(ρ^n,^n) = 0,
        ∀∈,

we obtain, by a similar computation to that in the proof of Lemma <ref>:

    12 δ t( ρ^n_ |ũ^n+1_|^2-ρ^n-1_ |_^n|^2 )
    + 12 |D_| ∑_=(ρ^n,^n) ^n+1_·^n+1_'
     - ((^n+1))_·_^n+1 + 1/^2( p)_^n·_^n+1 + R_^n+1
    =0,

where R_^n+1= 12δ tρ_^n-1|_^n+1-_^n|^2.
Dividing the velocity correction equation (<ref>) by (ρ^n_δ t )^1/2, we obtain:

    (ρ^n_/δ t)^1/2 _^n+1
    + ( δ t /ρ^n_)^1/2 1/^2  ( p^n+1)_
    = (ρ^n_/δ t)^1/2 _^n+1
    + (δ t/ρ^n_)^1/2 1/^2  (  p)_^n.

Squaring this relation and summing it with (<ref>) yields

    1 2δ t ( ρ_^n |_^n+1|^2 - ρ_^n-1 |_^n|^2  ) + 1/2|D_|∑_=D_|D_'(ρ^n,^n)  _^n+1·_'^n+1
      - ((^n+1))_·_^n+1 + 1/^2( p^n+1)_·_^n+1 + R_^n+1+ 1/^4 P_^n+1=0,

where P_^n+1 = δ t2 ρ_^n (  | ( p^n+1)_ |^2-  | ( p)^n_ |^2 ).
Recalling that the rescaled pressure is defined as 

    ( p)^n_ = (ρ^n_/ρ^n-1_)^1/2 ( p^n)_

for all ∈, we obtain (<ref>).


The discrete renormalization identities are again valid for the pressure correction algorithm, thanks to the fact that the mass balance (<ref>) is satisfied. The proof is identical to that of Lemma <ref> given for the implicit scheme.

 
A solution to the system (<ref>) satisfies for all K∈ and 0≤ n ≤ N-1 the two following identities:

    1 δ t ( ψ_γ(ρ^n+1_K)-ψ_γ(ρ^n_K)  ) + (ψ_γ(ρ^n+1) ^n+1)_K 
    +  (ρ_K^n+1)  ( ^n+1)_K + R_K^n+1 = 0, 
    
        1 δ t ( Π_γ(ρ^n+1_K)-Π_γ(ρ^n_K)  )
    + ((ψ_γ(ρ^n+1) - ψ_γ'(1) ρ^n+1) ^n+1)_K 
    
       
    + (ρ_K^n+1)  ( ^n+1)_K + R_K^n+1 = 0,

where R_K^n+1 has the same expression as in Lemma <ref>.


 
Let >0 and assume that the initial data (ρ_0^,_0^) is well-prepared in the sense of Definition <ref>,  satisfies (<ref>).
Then, for  small enough to ensure that ρ^-1 is positive, there exists a solution (ρ^n,^n)_0≤ n≤ N to the scheme (<ref>) such that ρ^n >0 for 1≤ n≤ N and the following inequality holds for 1≤ n≤ N-1:

    1/2∑_∈ |D_|  ( ρ^n_ |^n+1_|^2 -  ρ^n-1_ |^n_|^2  )  +  1/^2∑_K∈|K|  ( Π_γ(ρ^n+1_K) - Π_γ(ρ^n_K)  )
    
    
    + μ δ t ^n+1_1,^2 
    + δ t^2/2^4∑_∈|D_|  (  | ( p^n+1)_ |^2/ρ^n_   - | ( p^n)_ |^2/ρ^n-1_)  
    + ℛ^n+1≤ 0,

where ℛ^n+1 = ∑_∈ R^n+1_ + ^-2∑_K ∈ R_K^n+1≥ 0.



For  small enough, by Lemma <ref>, both ρ^-1 and ρ^0 are positive.
The positivity of ρ^n, for 1 ≤ n ≤ N, is a consequence of the properties of the upwind choice (<ref>) for ρ in the mass balance of the correction step (<ref>).

After multiplication by |D_|, we sum the kinetic energy balance equation (<ref>) over the faces, and after multiplication by ^-2|K|, we sum the relative entropy balance (<ref>) over the primal cells, and finally sum the two obtained relations.
Since the discrete gradient and divergence operators are dual with respect to the ^2 inner product (see (<ref>)), noting that the conservative fluxes vanish in the summation and that the diffusion term is coercive (see Lemma <ref>), we get (<ref>).

Given discrete density and velocity fields (ρ^n-1,ρ^n,^n), the predicted velocity ^n+1 is the solution of the linear system (<ref>).
Multiplying (<ref>) by |D_|, summing over the faces and using Young's inequality, yields the following inequality, which is valid for all α >0:

    1/2δ t√(ρ^n)^n+1_^2^2 
    - α/2^2^n+1_^2^2
    + μ ^n+1_1,^2
    ≤1/2δ t√(ρ^n-1)^n_^2^2 + 1/2α^2( p)^n_^2^2.

In this inequality, ρ^n and ρ^n-1 are the piecewise constant fields equal to ρ_^n and ρ_^n-1 on each dual face D_. Choosing α>0 small enough yields the coercivity of the linear system (<ref>) and therefore the existence of a unique solution ^n+1 to the prediction step. The existence of a solution (ρ^n+1,^n+1) to the correction step (<ref>)-(<ref>) follows from the Brouwer fixed point theorem, by an easy adaptation of the proof of <cit.>.


Let us now turn to the global entropy inequality.


Let >0 and assume that the initial data (ρ_0^,_0^) is well-prepared in the sense of Definition <ref>,  satisfies (<ref>).
By Lemma <ref>, there exists a solution (ρ^n,^n)_0≤ n≤ N to the scheme (<ref>). 
In addition, there exists C>0 independent of  such that, for  small enough and for all 1≤ n≤ N:

    1/2∑_∈ |D_|  ρ^n-1_ |^n_|^2
    + μ ∑_k=1^n δ t ^k_1,^2
    + 1/^2∑_K∈|K|  Π_γ(ρ^n_K)
    + δ t^2/2^4∑_∈|D_|/ρ^n-1_ |( p^n)_|^2
    ≤ C.




Summing (<ref>) over n yields the expected inequality (<ref>) with

    C= 1/2∑_∈ |D_|  ρ^-1_ |^0_|^2  +1/^2∑_K∈|K|  Π_γ(ρ^0_K)
    + δ t^2/2^4∑_∈|D_|/ρ^-1_ |( p^0)_|^2.

Let us prove that if the initial data is well-prepared in the sense of (<ref>), then for  small enough, C is uniformly bounded independently of . 
By Lemma <ref>, ρ_K^-1 is bounded for all K∈ for  small enough and therefore so is ρ_^-1 for all ∈. 
Hence, since _0^ is uniformly bounded in ^1(Ω)^d by (<ref>), a classical trace inequality yields the boundedness of the first term. By (<ref>), one has |ρ_K^0-1|≤ C^2 for all K∈. 
Hence, by (<ref>), the second term vanishes as → 0. The third term is also uniformly bounded with respect to  thanks to (<ref>).


 
Let >0 and assume that the initial data (ρ_0^,_0^) is well-prepared in the sense of Definition <ref>,  satisfies (<ref>).
By Lemma <ref>, for  small enough, there exists a solution (ρ^n,^n)_0≤ n≤ N to the scheme (<ref>).
Let p^n=(ρ^n) and define ^n = {_K^n, K∈} where _K^n = (p_K^n-m(p^n))/^2 with m(p^n) the mean value of p^n over Ω (i.e. m(p^n) = |Ω|^-1∑_K∈ |K| p_K^n).
Then, one has, for 1 ≤ n ≤ N:

    ^n≤ C_,δ t,

where  the real number C_,δ t depends on the mesh and the time step but not on , and · stands for any norm on the space of discrete functions.



According to (<ref>), the discrete pressure gradient is controlled in ^∞ by C_,δ t ^2 where C_,δ t is a mesh-dependent constant independent of . Hence, by a finite-dimensional argument (δ p^n) is controlled in any norm by some constant independent of .
In particular, on has 
(δ p^n)_-1,≤ C_,δ t, where . _-1, is the discrete (^-1)^d-norm defined for ∈_,0(Ω) by:

    _-1, = sup{∑_∈ |D_| _·_,
      with ∈_,0(Ω)  such that _1,≤ 1}.

Invoking the gradient divergence duality (<ref>) and the inf-sup stability of the scheme (see Lemma <ref>), (δ p^n)_-1,≤ C_,δ t implies that ^n_^2≤β^-1 C_,δ t.






 §.§ Incompressible limit of the pressure correction scheme


We may now state the main result of this section which is the convergence, up to a subsequence, of the solution to the compressible scheme (<ref>) towards the solution of a pressure correction inf-sup stable scheme when the Mach number tends to zero.

 


Let ()_m∈ be a sequence of positive real numbers tending to zero, and let (ρ,) be a corresponding sequence of solutions of the scheme (<ref>).
Let us assume that the initial data (ρ_0^,_0^) is well-prepared in the sense of Definition <ref>,  satisfies (<ref>) for m ∈.
Then the sequence (ρ)_m ∈ tends to the constant function ρ =1 when m tends to + ∞ in ^∞((0,T),^γ(Ω)). Moreover, for all q∈[1,min(2,γ)], there exists C>0 such that:

    ρ-1_^∞((0,T);^q(Ω))≤ C,     for m large enough.

In addition, the sequence (,)_m ∈ tends, in any discrete norm, to the solution (,) of the usual (Rannacher-Turek or Crouzeix-Raviart) pressure correction scheme for the incompressible Navier-Stokes equations, which reads:

0.8
Knowing  ^n∈_(Ω) and ^n∈_,0(Ω), compute ^n+1∈_(Ω) and ^n+1∈_,0(Ω) through the following steps:



    1 δ t (^n+1_-_^n )
    + (^n ⊗^n+1 )_
    -  ((^n+1))_
    + (^n)_
    =0,    ∀∈.
    
    
    
       
          1 δ t  (^n+1_-_^n+1)
    + (^n+1)_- (^n)_
    =0,    ∀∈.
    
    
         (^n+1)_K = 0,    ∀ K ∈.





The proof of the convergence of (ρ)_m ∈ towards ρ =1 when m tends to + ∞ is the same as in the continuous case. It follows from the combination of the lower bounds on Π_γ(ρ) proven in Lemma <ref> and the global entropy estimate (<ref>).

Since ρ→ 1 as m→+∞, the discrete density unknowns tend to 1 in every cell. Hence, by the kinetic energy part of (<ref>), there exists a minimum density ρ_ min,>0 depending on the fixed mesh such that _^2(Ω)^d≤√(2C/ρ_ min,) for all m∈. Moreover, the sequence ()_m ∈ is also bounded by Lemma <ref>.
By the Bolzano-Weierstrass theorem and a norm equivalence argument, there exists a subsequence of (,)_m ∈ which tends, in any discrete norm, to a limit (,).
Passing to the limit cell-by-cell in (<ref>), one obtains that (,) is a solution to (<ref>).
Since this solution is unique, the whole sequence converges, which concludes the proof.



As for the implicit scheme, one can extend the result of Theorem <ref> to the inviscid case where μ=λ=0. Indeed, the bound on the sequence of velocities is obtained thanks to the control by the estimate (<ref>) of the kinetic energy part of the global entropy, invoking the convergence for the densities towards 1.
 





§ ASYMPTOTIC ANALYSIS OF THE ZERO MACH LIMIT FOR A SEMI-IMPLICIT SCHEME
 

Another interesting semi-implicit scheme for low viscosity flows (typically a viscosity μ which is of the same order of magnitude as the space step h_) is a scheme, where in the momentum equation, only the pressure is treated in an implicit way (which is mandatory for stability reasons). 

Let δ t>0 be a constant time step. 
The approximate solution (ρ^n,^n)∈_(Ω)×_,0(Ω) at time t_n=nδ t for 1≤ n≤ N=T/δ t is computed through the following algorithm:

0.9
Knowing (ρ^n,^n)∈_(Ω)×_,0(Ω), solve for ρ^n+1∈_(Ω) and ^n+1∈_,0(Ω):



    1 δ t(ρ^n+1_K-ρ^n_K) + (ρ^n+1^n+1)_K = 0,    ∀ K ∈,
    
        1 δ t (ρ^n_^n+1_-ρ^n-1__^n )
    + (ρ^n ^n ⊗^n)_^ up
    -  ((^n))_
    + 1^2 ( p^n+1)_
    =0,    ∀∈.



Computing ^n+1_ from the second equation and inserting in the first one yields a nonlinear diffusion-convection-reaction problem on the density ρ^n+1 or the pressure p^n+1.
This algorithm may also equivalently be set under the form of a pressure-correction scheme:


    1 δ t (ρ^n_^n+1_-ρ^n-1__^n )
    + (ρ^n  ^n ⊗^n)_^ up
    -  ((^n))_
    =0,    ∀∈.
    
    
    
       
          1 δ t ρ^n_ (^n+1_-_^n+1)
    + 1^2 ( p^n+1)_
    =0,    ∀∈.
    
    
          1 δ t(ρ^n+1_K-ρ^n_K) + (ρ^n+1^n+1)_K = 0,    ∀ K ∈.



The ^2-stability of this scheme is expected to be ensured under a CFL restriction on the time step of the form δ t ≤ c(h_/||+h_^2/μ), provided that an upwind space discretization be used in the convection term of the momentum balance equation.
Hence, for the semi-implicit scheme (<ref>), this convection term is defined as:

    (ρ⊗)_^ up = ∑_∈(D_)(ρ,) _,     ∀∈,

where the dual mass fluxes (ρ,) are defined as previously in Section <ref>, and where the approximation of the convected velocity on a dual face =D_|D_' is computed with the upwind technique: _ =_ if (ρ,)≥ 0, _ = _' otherwise.

The initial data (ρ_0^,_0^) is assumed to satisfy ρ_0^ >0, ρ_0^∈^∞(Ω), _0^∈_0^1(Ω)^d and 

    _0^_^1(Ω)^d + 1/  _0^_^2(Ω) + 1/ ρ^_0- 1_^∞(Ω)≤  C,

for a real number C independent of . 
Moreover, we suppose that _0^ converges in ^2(Ω)^d towards a function _0∈^2(Ω)^d.
Note that these hypotheses are less restrictive than assuming well-prepared initial data since the density is here assumed to be close to 1 with an  rate (and not ^2).

The initialization of the scheme (<ref>) is performed similarly to that of the pressure-correction scheme.
First, ρ^0 and ^0 are given by the average of the initial conditions ρ_0^ and _0^ respectively on the primal cells and on the faces of the primal cells; then we compute ρ^-1 by solving the mass balance equation (<ref>) for n=-1, where the unknown is ρ^-1 and not ρ^0.
As in the case of the pressure correction scheme, one may prove that the discrete densities ρ_K^-1, K∈, are close to 1 for  small enough and therefore positive and bounded.





 §.§ A priori estimates


The discrete renormalization identities are again valid for the scheme (<ref>) , thanks to the fact that the mass balance (<ref>) is satisfied. The proof is identical to that of Lemmas <ref> and <ref>.

 
A solution to the system (<ref>) satisfies for all K∈ and 0≤ n ≤ N-1 the two following identities:

    1 δ t ( ψ_γ(ρ^n+1_K)-ψ_γ(ρ^n_K)  ) + (ψ_γ(ρ^n+1) ^n+1)_K 
    +  (ρ_K^n+1)  ( ^n+1)_K + R_K^n+1 = 0, 
    
        1 δ t ( Π_γ(ρ^n+1_K)-Π_γ(ρ^n_K)  )
    + ((ψ_γ(ρ^n+1) - ψ_γ'(1) ρ^n+1) ^n+1)_K 
    
       
    + (ρ_K^n+1)  ( ^n+1)_K + R_K^n+1 = 0,

where R_K^n+1 has the same expression as in Lemma <ref>.


Let us now turn to the discrete kinetic energy balance which leads to the ^2 stability of the scheme.
We begin with the following proposition which characterizes the rigidity matrix associated with the discretization of the diffusion term in the momentum equation.

 
Let A be the d♯× d♯ rigidity matrix associated with the finite element discretization of the diffusion term.  The matrix A is the block matrix A=(A_,')_,'∈, where for ,'∈, A_,' is the d× d matrix defined by:

    A_,' = ∑_K∈ ( μ∫_K ζ_·ζ_ '
    + (μ+λ) ∫_K ζ_⊗ζ_' ),

and A is a symmetric positive-definite matrix.



Let us identify =(_)_∈ and =(_)_∈ with vectors in ^d♯ and denote (.,.) the canonical scalar product in  ^d♯. The result follows from the elementary computation:

    (A,) 
        = ∑_∈∑_'∈A_,'_'·_
        = ∑_∈ |D_| ( -  (())_ ) ·_
        = ∑_K∈μ∫_K  : + (μ+λ) ∫_K ( )( ).

The expression (<ref>) of the block-entries of the matrix A then follows by developping the expression of  and  on the basis of the shape functions and identifying the terms.


The following lemma states the ^2-stability of an explicit upwind scheme for the convection-diffusion equation (<ref>), under some CFL restriction on the time step.

 
Any solution to (<ref>), and thus to the explicit scheme (<ref>), satisfies the following inequality, for all  0≤ n ≤ N-1:

    1 2δ t∑_∈ |D_|  ( ρ_^n |_^n+1|^2 - ρ_^n-1 |_^n|^2 ) 
    +  R_^n+1≤ 0,

where, denoting ϱ(A) the spectral radius of the matrix A, the remainder term  R_^n+1 is given by:

    R_^n+1= ∑_∈ ( 
    ρ_^n |D_|/2 δ t 
    -  12 ∑_∈(D_) ( (ρ^n,^n) )^-
    -  14 ϱ(A)
     ) |_^n+1-_^n|^2.




Taking the scalar product of the momentum equation (<ref>) with |D_| _^n+1, we obtain T^ conv_ + T^ diff_=0, with:

    T^ conv_ =  ( |D_|δ t (ρ^n_^n+1_-ρ^n-1__^n )
    + ∑_∈(D_)(ρ^n,^n) _^n  ) ·_^n+1,
    
       
    T^ diff_ = - |D_|   ((^n))_·_^n+1.

For the convection term, the dual density unknowns and mass fluxes are chosen so as to have:

    |D_|/δ t (ρ^n_-ρ^n-1_)+ ∑_∈(D_)(ρ^n,^n) = 0,     ∀∈,

which necessitates a time-shift in order to exploit the mass balance at the previous time step.
Hence, following the proof of <cit.>, we obtain:

    T^ conv_ 
    = |D_|2δ t ( ρ_^n |_^n+1|^2 - ρ_^n-1 |_^n|^2 ) 
    + 1/2∑_∈(D_)(ρ^n,^n) |_^n|^2 
    + |D_|2δ tρ_^n |_^n+1-_^n|^2
    
    
    - 1/2∑_∈(D_)(ρ^n,^n) |_^n-_^n|^2 
    +  ∑_∈(D_)(ρ^n,^n) (_^n-_^n) · (_^n+1-_^n).

Thanks to the upwind choice for _^n, the term _^n-_^n vanishes whenever the dual flux (ρ^n,^n) is non-negative. Applying Young's inequality to the product in the last term yields:

    T^ conv_≥|D_|2δ t ( ρ_^n |_^n+1|^2 - ρ_^n-1 |_^n|^2 ) 
    + 1/2∑_∈(D_)(ρ^n,^n) |_^n|^2 
    + |D_|2δ tρ_^n |_^n+1-_^n|^2
    
    - 1/2∑_∈(D_)(ρ^n,^n)^- |_^n+1-_^n|^2.

For the diffusion term, we observe that:

    ∑_∈ T^ diff_ = (A^n,^n+1) = (A^n,^n+1-^n) + (A^n,^n).

The Cauchy-Schwarz inequality for the scalar product associated with the real positive symmetric matrix A yields

    (A^n,^n+1-^n) ≥
    - (A^n,^n)^1/2 (A(^n+1-^n),^n+1-^n)^1/2.

Applying Young's inequality yields:

    ∑_∈ T^ diff_ ≥  - 1/4 (A(^n+1-^n),^n+1-^n)
     ≥  - 1/4ϱ(A)  ∑_∈ |_^n+1-_^n|^2.

By the conservativity of the dual fluxes, taking the sum of (<ref>) over ∈ and summing the result with (<ref>) yields the expected inequality (<ref>).


We are now in position to state the following local-in-time kinetic energy balance.

 
Any solution to the explicit scheme (<ref>) satisfies the following inequality, for all  0≤ n ≤ N-1:

    1 2δ t∑_∈ |D_|  ( ρ_^n |_^n+1|^2 - ρ_^n-1 |_^n|^2 ) 
    + 1/^2∑_∈ |D_| ( p^n+1)_·_^n+1
    
    + δ t/2 ^4∑_∈|D_|/ρ^n_ |( p^n+1)_|^2
    + R_^n+1≤ 0,

where R_^n+1 is defined in Lemma <ref>.



As for the first pressure correction scheme, we write the velocity correction equation as:

    (ρ^n_/δ t)^1/2 _^n+1
    + ( δ t /ρ^n_)^1/2 1/^2  ( p^n+1)_
    = (ρ^n_/δ t)^1/2 _^n+1,

square this relation and sum it with (<ref>), which yields the desired inequality.


We may now state the following result.



Let >0 and assume that the initial data satisfies (<ref>).
Then, for  small enough to ensure that ρ^-1 is positive, there exists a solution (ρ^n,^n)_0≤ n≤ N to the scheme (<ref>), and for 1≤ n≤ N, the density ρ^n is positive.
Moreover, assuming that the time-step satisfies the following CFL restriction:

    δ t ≤min_∈ 4 ρ_^n |D_|/ 2 ∑_∈(D_) ( (ρ^n,^n) )^- + ϱ(A),

the following inequality holds:

    1/2∑_∈ |D_|  ( ρ^n_ |^n+1_|^2 -  ρ^n-1_ |^n_|^2  )
    + 1/^2∑_K∈ |K|  ( Π_γ(ρ^n+1_K) - Π_γ(ρ^n_K)  )
    
    
    + δ t^2/2 ^4∑_∈|D_|/ρ^n_ |( p^n+1)_|^2
    + ℛ^n+1≤ 0,

where ℛ^n+1 = R^n+1_ + ^-2∑_K ∈ R_K^n+1≥ 0.



The positivity of the density is a consequence of the properties of the upwind choice (<ref>) for ρ in the mass balance of the correction step (<ref>).
After multiplication by ^-2|K|, we sum the renormalization identity (<ref>) over the primal cells, and sum the obtained relation with the kinetic energy balance (<ref>).
Since the discrete gradient and divergence operators are dual with respect to the ^2 inner product (see (<ref>)), we get (<ref>). Under the CFL condition (<ref>), the remainder term R^n+1_ defined in (<ref>) is non-negative.
The existence of a solution (ρ^n+1,^n+1) to the scheme (<ref>) follows form the Brouwer fixed point theorem, by an easy adaptation of the proof of <cit.>.


The restriction (<ref>) on the time-step is a convective-diffusive CFL condition. In particular, it is satisfied if the time step simultaneously satisfies the following two conditions:

    δ t ≤min_∈ ρ_^n |D_|/∑_∈(D_) ( (ρ^n,^n) )^-  and  δ t ≤min_∈ 2 ρ_^n |D_|/ϱ(A).

The first condition is a convective CFL restriction associated with the velocity of the fluid, which is consistent with the explicit upwind discretization of the momentum convection term.
As for the second condition, it is also a classical diffusive CFL restriction associate with an explicit treatment of the diffusion term.
Indeed, if the Lamé coefficients μ and λ are of the same order of magnitude, usual regularity assumptions on the mesh imply that: 

    ϱ(A) ∝μ  h_^d-2,

where h_ is a characteristic length of the mesh cells.
Therefore, if the viscosity μ is small, typically μ≈ h_, the diffusive CFL restriction is comparable to the convective CFL restriction.


Let >0 and assume that the initial data (ρ_0^,_0^) satisfies (<ref>).
By Lemma <ref>, for  small enough, there exists a solution (ρ^n,^n)_0≤ n≤ N to the scheme (<ref>). Moreover, if the time-step satisfies the CFL restriction (<ref>) for all 1≤ n≤ N , then there exists C>0 independent of  such that, for  small enough and for all 1≤ n≤ N :

    1/2∑_∈ |D_|  ρ^n-1_ |^n_|^2
    + 1/^2∑_K∈|K|  Π_γ(ρ^n_K)
    + δ t^2/2 ^4∑_k=1^n ∑_∈|D_|/ρ^k-1_ |( p^k)_|^2
    ≤ C.




Multiplying equation (<ref>) by δ t and summing over the time steps yields for 1≤ n≤ N:

    1/2∑_∈ |D_|  ρ^n-1_ |^n_|^2
    + 1/^2∑_K∈|K|  Π_γ(ρ^n_K)
    + δ t^2/2 ^4∑_k=1^n ∑_∈|D_|/ρ^k-1_ |( p^k)_|^2
    + ℛ^n
    
    ≤ 1/2∑_∈ |D_|  ρ^-1_ |^0_|^2
    +1/^2∑_K∈|K|  Π_γ(ρ^0_K),

with ℛ^n= ∑_k=0^n-1 (R^k+1_ + ^-2∑_K ∈ R_K^k+1 ) ≥ 0.

Let us prove that the right hand side of (<ref>) is uniformly bounded for all  small enough.
As in Lemma <ref>, one may prove that under (<ref>), ρ_K^-1 is bounded for all K∈ for  small enough and therefore so is ρ_^-1 for all ∈. 
Hence, since _0^ is uniformly bounded in ^1(Ω)^d by (<ref>), a classical trace inequality yields the boundedness of the first term. By (<ref>), one has |ρ_K^0-1|≤ C for all K∈. 
Hence, by (<ref>), the second term is also uniformly bounded with respect to .


From now on, we need to assume that we use a time step independent of the Mach number and satisfying the constraint (<ref>) (which involves the velocity field, which itself depends on ).
The existence of such a time step is proven in Appendix <ref>.

Under this assumption, by the same arguments as for the pressure correction scheme, we get that the dynamic pressure is controlled independently of ; this is stated in the following lemma.



Let >0 and assume that the initial data (ρ_0^,_0^) satisfies (<ref>).
Then, there exists a solution (ρ^n,^n)_0≤ n≤ N to the scheme (<ref>).
Let  p^n=(ρ^n) and define ^n = {_K^n, K∈} where _K^n = (p_K^n-m(p^n))/^2 with m(p^n) the mean value of p^n over Ω.
If the time step δ t satisfies the CFL condition (<ref>) for all 1≤ n≤ N independently of the Mach number , then, one has, for all 1≤ n ≤ N:

    ^n≤ C_,δ t,

where  the real number C_,δ t depends on the mesh and the time step but not on , and · stands for any norm on the space of discrete functions.






 §.§ Incompressible limit of the scheme


By the same proof as for the previous pressure correction scheme, we have the following convergence result. An analogous result is also valid in the inviscid case μ=λ=0.

 

Let ()_m∈ be a sequence of positive real numbers tending to zero.
Let the associated sequence of initial data (ρ_0^,_0^) satisfy (<ref>), and the time step satisfy the CFL condition (<ref>) for all 1≤ n≤ N independently of m.
Let (ρ,) be a corresponding sequence of solutions of the scheme (<ref>).
Then the sequence (ρ)_m ∈ tends to the constant function ρ =1 when m tends to + ∞ in ^∞((0,T),^γ(Ω)). Moreover, for all q∈[1,min(2,γ)], there exists C>0 such that:

    ρ-1_^∞((0,T);^q(Ω))≤ C,     for m large enough.


In addition, the sequence (, δ p)_m ∈ tends, in any discrete norm, to a limit (,) which is the solution to the following inf-sup stable semi-implicit scheme for the incompressible Navier-Stokes equations: 

Knowing  ^n∈_(Ω) and ^n∈_,0(Ω), solve for ^n+1∈_(Ω) and ^n+1∈_,0(Ω) :


    (^n+1)_K = 0,    ∀ K ∈,
    
        1 δ t (^n+1_-_^n )
    + (^n ⊗^n)_^ up
    -  ((^n))_
    + (^n+1)_
    =0,    ∀∈.








§ NUMERICAL TESTS
 

We assess the convergence of the scheme on a test case built for this purpose.
Analytical solutions of the 2D barotropic Euler and Navier-Stokes equations are obtained through the following steps: we first derive a compactly supported H^2(^2) solution of the stationary barotropic Euler equations consisting in a standing vortex; then, we obtain a time-dependent solution of the Euler equations by adding a constant velocity motion; finally, we pass to the Navier-Stokes equations by compensating the viscous forces at the right-hand side.
The velocity field of the standing vortex is sought under the form:

    = f(ξ) [ -x_2;  x_1 ],  ξ = x_1^2 + x_2^2.

A simple derivation of this expression yields:

    (·) = -f(ξ)^2 [ x_1; x_2 ].

Let us seek the pressure under the form p̂=g(ξ), so:

    p̂ = 2  g'(ξ) [ x_1; x_2 ].

Assuming a perfect gas pressure law, the density is given by ρ̂= g(ξ)^1/γ.
By construction, these functions satisfy (ρ̂)=0. Moreover, one has ρ̂(·) +p̂ = 0 if, and only if -g(ξ)^1/γ f(ξ)^2 + 2  g'(ξ) =0 for all ξ.
We thus obtain a solution of the stationary Euler equation if g takes the following expression:

    g=(γ-1/2γ (F +c_M) )^γ/(γ-1),

where F is such that F'=f^2 and c_M is a positive real number.
For the present numerical study, we choose γ=3 and f(ξ)=10 ξ^2 (1-ξ)^2 if ξ∈ (0,1), f=0 otherwise, which indeed yields an H^2(^2) velocity field.
The associated expression of F is:

    F(ξ)=100 (1/5ξ^5 -2/3ξ^6 + 6/7ξ^7 - 1/2ξ^8 +1/9ξ^9) ξ∈ (0,1),
     F(ξ) = F(1) = 10/63.

The problem is made unstationary by a time translation: given a constant vector field , the density ρ and the velocity   are deduced from the steady state solution ρ̂ and  by ρ(,t)=ρ̂(- t) and (,t)=(- t)+.
The center of the vortex is initially located at _0=(0,0)^t, the translation velocity  is set to =(1,1)^t, the computational domain is Ω=(-1.2, 2.8)^2 and the computation is run on the time interval (0,0.8).
We perform several computations keeping the velocity constant (and therefore of order 1, according to the expression of f) and varying the Mach number by changing the constant c_M and therefore also the pressure level and the speed of sound (given by c^2=γ  p ^(γ-1)/γ).
Choosing 1 as the reference value for the velocity and the speed of sound outside the vortex as the reference speed of sound, c_M and the Mach number are connected as given in Table <ref>:



Computations are run with the open-source CALIF^3S software developed at IRSN <cit.>, with the pressure correction algorithm described in Section <ref>.
The mesh is a 500× 500 uniform grid, and the time step is set at the same value as the space step,  δ t = 0.008, for all the computations, which corresponds to a CFL number with respect to the material velocity close to 1.5.



  
Euler equations – In the Euler case, since the convection term in the momentum balance is approximated with a centered discretization, the computations are stabilized by taking into account an artificial viscosity given by:

    μ_a = ρ_ext  v_max  h/10,

where ρ_ext stands for the density outside the vortex (which depends on the computation), v_max =1.4 is an approximation of the maximal value of the components of the velocity and h is the space step.
This viscosity is in the order of a fifth of the upwinding-induced numerical viscosity.
We plot on Figure <ref> the second component of the velocity obtained at t=0.8 along the line x_2=0.8, which crosses the center of the vortex.
The results are almost independent of the Mach number (in fact, the curves are superimposed on the figure).
To check the differences, a zoom of the curves near the minimum value of the velocity is shown on Figure <ref>; the maximum of the differences is close to 0.0016, while the amplitude of the analytical velocity variation is equal to 1.





In addition, the expression of the pressure shows that F is much lower than c_M, so a Taylor development shows that the quantity δ p defined by

    δ p =(p-p_ext)/c,

with p_ext the pressure outside the vortex and c the speed of sound given in the above table, should be approximatively independent of the Mach number.
This quantity is plotted on Figure <ref>, which shows that it is indeed the case (the observed discrepancy when the Mach number is close to 1 may be attributed to the fact that F/c_M takes in this case its greatest value, and the Taylor development is less accurate).

Finally, the L^1 norm of the difference between the numerical velocity and the piecewise constant function obtained by taking, on each diamond cell, the value of the continuous solution at the cell mass center is, at t=0.8: 0.192 for a Mach number Ma close to 1, 0.189 for Ma=0.1 and 0.187 for the other values of the Mach number.
For the pressure, the same discrete L^1 norm of the difference between the numerical and analytical solutions scales as the magnitude of the pressure variations, which, in turn, as said before, scales as the speed of sound c.
The ratio between this norm and c reads: 0.0168, 0.0223, 0.0229, 0.0227 and 0.0264, for the tested Mach numbers from 1 to 0.0001.
The value slightly greater obtained for Ma=0.0001 is probably due to the fact that the pressure is so high in this case that the algebraic solvers become less accurate.





  
Navier-Stokes equations – The same analytical solution is used for the Navier-Stokes equations, with a viscosity now given by:

    μ = ρ_ext/50,

so a Reynolds number in the range of 50 (equal to 50 if the characteristic velocity range is set to 1, equal to 75 if its is set to v_max=1.5).
The corresponding viscous term is compensated by a source term at the right-hand side of the momentum balance equation.
The same curves as for the Euler case are shown on Figures <ref>-<ref>.
The numerical error for the velocity (as defined before) is almost independent of the Mach number: it always falls in the interval (0.099,0.1).
For the pressure, the same scaled quantity as before reads: 0.0108, 0.0138, 0.0147, 0.0146 0.0183, for the tested Mach numbers from 1 to 0.0001.
The conclusions are thus the same, up to the minor difference that the results in the Navier-Stokes case are slightly more accurate.













§ EXISTENCE OF A TIME STEP SATISFYING THE CFL CONDITION (<REF>) FOR ALL MACH NUMBERS
 

In this appendix, we prove that for the semi-implicit scheme (<ref>), it is possible to define a time-step δ t which satisfies the CFL condition (<ref>) independently of the Mach number .
This allows to prove the discrete global entropy estimate (<ref>) independently of the Mach number.
More precisely, we prove that under a convective CFL condition which is more restrictive than expected (δ t≤ C h_^1+d/2 instead of δ t ≤ C h_ with C independent of ), the discrete global entropy estimate (<ref>) holds true.

Let θ_ be a measure of the regularity of the mesh in the classical finite element sense (see for instance <cit.>).
Then for all K∈ and ,'∈(K), we have |'|/|D_|≤ C_1(θ_)   h_^-1 for some nondecreasing function  C_1. Moreover, thanks to classical inverse inequalities, there exists a nondecreasing function C_2 such that  _^∞(Ω)^d≤ C_2(θ_) h_^-d/2_^2(Ω)^d for all ∈_(Ω).

Since the initial data satisfies (<ref>), the initial total energy is uniformly bounded with respect to .
We denote by C_0 a uniform upper bound (which depends on the constant C in (<ref>) and on the spatial discretization ), which thus satisfies for all >0:

    1/2∑_∈ |D_|  ρ^-1_ |^0_|^2
    + 1/^2∑_K∈|K|  Π_γ(ρ^0_K) ≤ C_0.

We define:

    C(d,θ_,C_0) := 1/2√(2)d C_1(θ_) C_2(θ_) C_0^1/2.

We have the following result.


Let η be a fixed small parameter in (0,1). If the time step δ t satisfies the following CFL condition, which is independent of the Mach number:

    δ t ≤ (1-η) min ( C(d,θ_,C_0) h_^1+d/2 ,
     min_∈ 2 |D_|/ϱ(A) ),

then there exists  (depending on γ, the spatial discretization , the constant C in (<ref>) and on η) such that for all ∈(0,):

    max_n=0,..,N {1/2∑_∈ |D_|  ρ^n-1_ |^n_|^2
    +1/^2∑_K∈|K|  Π_γ(ρ^n_K) }≤ C_0.




As a by-product of the proof, we can see that, if the time step δ t satisfies the CFL condition (<ref>), then it satisfies the classical CFL condition (<ref>) independently of the Mach number , and this also allows to prove (<ref>).
Actually the proof given hereunder consists in proving (<ref>) and (<ref>) simultaneously, thanks to an induction process.



Let δ t satisfy the CFL condition (<ref>).
Let (α^n_)_n=0,..,N and (β^n_)_n=0,..,N be the two sequences defined by:

    α^n_ := min ( 
     min_∈ ρ_^n/  |D_|^-1∑_∈(D_) ( (ρ^n,^n) )^- , min_∈ 2 ρ_^n /  |D_|^-1ϱ(A) ),
       β^n_ := 1/2∑_∈ |D_|  ρ^n-1_ |^n_|^2+1/^2∑_K∈|K|  Π_γ(ρ^n_K).

By Lemma <ref>, we know that for all 0≤ n ≤ N-1 and for all >0, δ t ≤α^n_⟹β^n+1_≤β^n_. Let us prove that for all 0≤ n ≤ N-1:

    (∃ ^n >0  such that β^n+1_≤β^n_≤ .. ≤β^0_≤ C_0   for all ∈(0,^n)) 
    
    ⟹     ( ∃ ^n+1>0  with ^n+1≤^n  such that δ t ≤α^n+1_ for all ∈(0,^n+1) ).


Let us assume that the left-hand side of the above implication holds true for some n such that 0 ≤ n ≤ N-1.
Then, for every =K|L∈, we have the following inequalities:

    |D_|^-1∑_∈(D_) ( (ρ^n+1,^n+1) )^-   ≤ 2d max_'∈(K)∪(L) |D_|^-1  |F_K,'(ρ^n+1,^n+1)| 
       ≤ 2d max_'∈(K)∪(L)|'|/|D_|  |ρ^n+1_'|   |^n+1_'|
       ≤  2d  C_1(θ_) h_^-1 (max_K∈|ρ^n+1_K|  )  (max_∈ |^n+1_| ).

The first inequality in (<ref>) follows from hypothesis (H3) in (<ref>). 
We then remark that max_∈ |^n+1_|=^n+1_^∞(Ω)^d≤ C_2(θ_) h_^-d/2^n+1_^2(Ω)^d with

    ^n+1_^2(Ω)^d   ≤ ( 2/min_K∈ρ_K^n)^1/2 (   1/2∑_∈ |D_|  ρ^n_ |^n+1_|^2
    +1/^2∑_K∈|K|  Π_γ(ρ^n+1_K) )^1/2.

Hence, since the left-hand side of (<ref>) is assumed to hold true, we obtain that

    ^n+1_^2(Ω)^d≤ ( 2/min_K∈ρ_K^n)^1/2  C_0^1/2.

Injecting in (<ref>), we obtain:

    min_∈ ρ_^n+1/|D_|^-1∑_∈(D_)( (ρ^n+1,^n+1) )^-≥ 
    C(d,θ_,C_0)   h_^1+d/2 min_K∈ρ_K^n+1/max_K∈ρ^n+1_K (min_K∈ρ_K^n)^1/2.

Invoking once again the left hand side of (<ref>), we have β_^n+1≤β_^n ≤ C_0 for all ∈(0,^n), which implies that 
Π_γ(ρ^n+1)_^1(Ω)≤ C_0 ^2 and Π_γ(ρ^n)_^1(Ω)≤ C_0 ^2 for all ∈(0,^n). By the results of Lemma <ref>, this implies that for all K∈, ρ_K^n→ 1 and ρ_K^n+1→ 1 as →0. Hence, there exists 0<^n+1≤^n such that for all ∈(0,^n+1):

    min_∈ ρ_^n+1/|D_|^-1∑_∈(D_)( (ρ^n+1,^n+1) )^-≥ (1-η) C(d,θ_,C_0)   h_^1+d/2.

Obviously, ^n+1 depends on η.
It also depends on γ and on the spatial discretization (namely θ_ and h_) since one has to bound the ^∞-norm of ρ^n-1 (resp. of ρ^n+1-1) by the ^min(2,γ)-norm of ρ^n-1  (resp. of ρ^n+1-1) thanks to an inverse inequality, combined with the estimates of Lemma <ref>.
Upon diminishing ^n+1, we can also prove that for all ∈(0,^n+1):

    min_∈ 2 ρ_^n+1/  |D_|^-1ϱ(A)≥(min_∈ 2  |D_|/ϱ(A))  (min_K∈ρ_K^n+1) ≥
    (1-η)  min_∈ 2  |D_|/ϱ(A).

Combining (<ref>) and (<ref>), we obtain that δ t ≤α^n+1_ for all ∈(0,^n+1).
Hence, by a straightforward induction process, the proposition is proved  with := ^N-1 provided that there exists ^0 such that β^1_≤β_^0 for all ∈(0,^0).
It is sufficient to prove that for some ^0>0,  δ t ≤α_^0 for all ∈(0,^0). 
Following similar steps as above, this is easily proved with ^0 only depending on the constant C in (<ref>) and on the discretization.




abbrv

