Charlotte Perrin
email: charlotte.perrin@univ-amu.fr

Khaled Saleh
email: saleh@math.univ-lyon1.fr

Numerical staggered schemes for the free-congested Navier-Stokes equations

Keywords: Compressible Navier-Stokes equations, singular limit, staggered discretization, finite volume -finite element method, pressure correction scheme 35Q30, 65N12, 76M10, 76M12

In this paper we analyze two numerical staggered schemes, one fully implicit, the other semiimplicit, for the compressible Navier-Stokes equations with a singular pressure law; system which degenerates towards the free-congested Navier-Stokes equations with the density constraint 0 ≤ ρ ≤ 1 as the intensity of the pressure, ε, tends to 0. The two schemes are shown to "preserve" the limit ε → 0, in the following sense: via discrete energy and pressure estimates we control uniformly with respect to ε the discrete density, velocity and pressure; and, as ε → 0, the limit pressure is shown to be activated only in the congested region where the discrete limit density is maximal, equal to 1. Finally, the semi-implicit scheme is implemented via the software CALIF 3 S on 1D and 2D test cases and the numerical solutions are shown to capture well the transition between the congested regions {ρ = 1} and the free regions {ρ < 1}.

Introduction

The purpose of this paper is to present and analyze new numerical schemes for the following singular compressible Navier-Stokes equations: The viscous stress tensor reads: div (τ (u ε )) = µ∆u ε + (µ + λ)∇(div u ε ) with µ > 0, µ + λ > 0.

∂ t ρ ε + div(ρ ε u ε ) = 0, ∂ t (ρ ε u ε ) + div(ρ ε u ε ⊗ u ε ) + ∇p ε (ρ ε ) -div (τ (u ε )) = f , (ρ ε , (ρ ε u ε )) |t=0 = (ρ ε 0 , q ε 0 ) (1.1a) (1.1b) (1.
Finally, the pressure law depends on a parameter ε > 0 and is given by

p ε (ρ) = ε ρ γ (1 -ρ) β , γ, β > 1.
(1.3)

Due to the singularity, the pressure (1.3) plays the role of a barrier, it prevents the density to exceed the critical value 1: 0 ≤ ρ ε (t, x) < 1 a.e. (t, x) ∈ (0, T ) × Ω, (1.4) provided that the constraint is satisfied initially. From the modeling point of view, p ε can be seen as the average of microscopic repulsion forces that prevent the overlapping of the microscopic components of the fluid and the parameter ε represents the average intensity of these forces. As this latter decreases and tends to 0, the range of action of the repulsive forces shrinks. Heuristically, as ε → 0, the pressure p ε tends to a limit potential p which is activated only when the density is maximal, equal to 1. Hence, the following exclusion constraint holds: (1 -ρ)p = 0. From another point of view, observing that in the congested (or saturated) regions where ρ = 1 one must have div u ≥ 0 (otherwise the density could exceed 1 according to the mass conservation equation), the limit potential p can be seen as a Lagrange multiplier associated to the velocity constraint. Therefore, in the asymptotics ε → 0, solutions to the problem (1.1) should converge towards solutions of the following hybrid system:

∂ t ρ + div(ρu) = 0, ∂ t (ρu) + div(ρu ⊗ u) + ∇p -div (τ (u)) = f , 0 ≤ ρ ≤ 1, (1 -ρ)p = 0, p ≥ 0,
(ρ, (ρu)) |t=0 = (ρ 0 , q 0 ),

(1.5d) which couples a free (compressible, pressureless) dynamics in {ρ < 1} with a congested (incompressible) dynamics in {ρ = 1}, and will be consequently referred as the free-congested Navier-Stokes equations.

This limit system has important and various applications in physics, geophysics and biology. It has been first derived by Bouchut et al. in [START_REF] Bouchut | A hierarchy of models for two-phase flows[END_REF] from biphasic liquid/gas equations through the formal asymptotic r = ρgas ρ liq → 0. In another domain, several recent works have been dedicated to wave-structure interactions and floating body dynamics. In particular the reader is referred to [START_REF] Lannes | On the dynamics of floating structures[END_REF], [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF] (and the references given therein) which present a formulation of the dynamics like a free-congested model similar to (1.5). Finally, congestion models are used in the modelling of collective motions, see for instance [START_REF] Degond | Numerical simulations of the euler system with congestion constraint[END_REF], [START_REF] Maury | Prise en compte de la congestion dans les modeles de mouvements de foules[END_REF]- [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF]. Note that the approximation (1.1) is particularly relevant for collective motions, since the singular pressure p ε (1.3) models then the social (short-range) repulsive forces between the agents.

The singular limit ε → 0 and the convergence of (weak) solutions of (1.1)-(1.3) towards (weak) solutions of the free-congested system (1.5) has been proved in [START_REF] Perrin | Free/congested two-phase model from weak solutions to multidimensional compressible navier-stokes equations[END_REF]. As explained in Section 2 below, the convergence result relies on energy and pressure estimates combined with refined weak compactness arguments specific to the compressible Navier-Stokes equations. In particular, it can be rigorously shown that the incompressibility constraint divu = 0 is satisfied inside the congested domain {ρ = 1}. We can thus connect the limit ε → 0 with the famous low Mach number limit which classically formalizes the transition between the compressible regime and the incompressible regime. The important difference between the two singular limits lies in fact that the convergence towards the incompressible regime in (1.1)-(1.3) is not only local but also depends on the solution itself, which is not the case for the low Mach limit where the speed the sound in the material is supposed to tend to +∞ (relatively to the characteristic speed of the flow) uniformly in the whole domain.

From the numerical standpoint, it seems that only few numerical schemes have been proposed for the limit system (1.5), the difficulty associated to that free boundary problem being that no information is available on the interface (no explicit dynamics, no transmission conditions imposed). Let us mention the study [START_REF] Bouchut | A hierarchy of models for two-phase flows[END_REF] which presents in one dimension a numerical scheme based on a "sticky blocks" approximation of the dynamics of (1.5), and [START_REF] Maury | Pressureless euler equations with maximal density constraint: a timesplitting scheme[END_REF], [START_REF] Perrin | One-dimensional granular system with memory effects[END_REF] (restricted to 1D) which rely on a projection of the dynamics onto a set of admissible trajectories. Nevertheless, these three papers can only deal with inviscid flows µ = λ = 0.

The studies [START_REF] Degond | Numerical simulations of the euler system with congestion constraint[END_REF] and [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF] adopt a "pseudo-compressibility" strategy by approximating the compressibleincompressible system (1.5) (with zero viscosity) by a fully compressible system. For Degond et al. [START_REF] Degond | Numerical simulations of the euler system with congestion constraint[END_REF] the approximated system is the inviscid counterpart of (1.1), while in [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF] Godlewski et al. relax the constraint by replacing the incompressible pressure p by p λ (ρ) = (ρ-1)+ λ and let λ → 0. In both cases, ε → 0 or λ → 0, one faces numerical difficulties similar to the low Mach number limit. First, one observes that the CFL condition associated to a fully explicit scheme degenerates as ε → 0 (resp. λ → 0) in the dense regions. Indeed, the CFL, which is associated to the speed of propagation of acoustic waves c ε = p ′ ε (ρ ε ), blows up as ε → 0 and ρ ε → 1 (resp. c λ ∝ λ -1 for ρ λ > 1). On the other hand, the schemes for the compressible case usually use a collocated arrangement of the unknowns, and the diffusion operator appearing in the pressure wave equation is unstable, since collocated approximations do not satisfy a form of the so-called discrete inf-sup condition. Hence, we expect that an additional stabilization mechanism will be required as ε → 0 (or λ → 0) in the dense regions whose dynamics is nearly incompressible. To remedy to these problems, one can either adapt the classical compressible methods via preconditionning techniques [START_REF] Turkel | Preconditionning techniques in computational fluid dynamics[END_REF], or extend the incompressible methods to compressible regimes by adapting pressure-correction schemes [START_REF] Harlow | A numerical fluid dynamics calculation method for all flow speeds[END_REF]. The authors of [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equations[END_REF] propose a finite volume (Rusanov) scheme where the pressure is split into a bounded part, which is treated explicitly in time, and an implicit part which carries the singularity close to 1. This technique allows to circumvent the degeneracy of the CFL condition in the dense regimes. An analysis of the L 1 -error is performed on the numerical solutions for several 1D inviscid test cases based on Riemann problems. Physical viscosity is taken into account in the more recent study [START_REF] Degond | Transport of congestion in two-phase compressible/incompressible flows[END_REF] (which also deals with a variable congestion constraint), keeping a fully finite volume space discretization. In the study [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF], the authors implement a Centered-Potential Regularization scheme (see also [START_REF] Parisot | Centered-potential regularization for the advection upstream splitting method[END_REF]) which ensures the consistency the incompressible regime but also the well-balanced property (in the wave-structure interaction context, topography may be taken into account and numerical schemes are required to preserve the lake at rest solution). Nevertheless, the analysis of these schemes is quite far from the analysis of the underlying continuous equations performed in [START_REF] Perrin | Free/congested two-phase model from weak solutions to multidimensional compressible navier-stokes equations[END_REF], in particular no uniform control of the discrte pressure (wrt to ε or λ) seems available in [START_REF] Degond | Numerical simulations of the euler system with congestion constraint[END_REF], [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF].

In this paper, we consider a numerical scheme which is implemented in the industrial software CALIF 3 S [2] developed by the french Institut de Radioprotection et de Sûreté Nucléaire (IRSN, a research center devoted to nuclear safety). This scheme falls in the class of staggered discretizations in the sense that the scalar variables (density, pressure) are associated with the cells of a primal mesh M while the vectorial variables (velocity, external force) are associated with the set E of faces of the primal mesh. Such decoupling, associated here with a Crouzeix-Raviart or Rannacher-Turek finite element discretization [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations[END_REF] (but other non-conforming finite elements are possible, such as the Rannacher-Turek discretization [START_REF] Rannacher | Simple nonconforming quadrilateral Stokes element[END_REF] on quadrangles (2D) or hexahedra (3D)) of the viscous stress tensor, provides a discrete pressure estimate, thanks to the so-called discrete inf-sup stability condition (see for instance [START_REF] Girault | Finite element methods for Navier-Stokes equations: theory and algorithms[END_REF]). This condition, which is also satisfied by the MAC scheme (see [START_REF] Harlow | Numerical calculation of almost incompressible flow[END_REF], [START_REF] Harlow | A numerical fluid dynamics calculation method for all flow speeds[END_REF], [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF]) on structured grids, ensures the unconditional stability of the scheme in almost incompressible regimes (for instance in the low Mach regime, see [START_REF] Gallouët | An unconditionally stable pressure correction scheme for compressible barotropic Navier-Stokes equations[END_REF] and [START_REF] Herbin | Low mach number limit of some staggered schemes for compressible barotropic flows[END_REF]). A finite element discretization for the viscous stress tensor is here coupled with finite volume discretizations of the convective terms which allow, thanks to standard techniques, to obtain discrete convection operators satisfying maximum principles (e.g. [START_REF] Larrouturou | How to preserve the mass fractions positivity when computing compressible multi-component flows[END_REF]). The discrete mass convection operator is a standard finite volume operator defined on the cells of the primal mesh M while the discrete momentum convection operator is also a finite volume operator written on dual cells, i.e. cells centered at the location of the velocity unknowns, namely the faces E. A difficulty implied by such staggered discretization lies in the fact that, as in the continuous case, the derivation of the energy inequality needs that a mass balance equation be satisfied on the same (dual) cells, while the mass balance in the scheme is naturally written on the primal cells. A procedure is therefore needed to define the density on the dual mesh cells and the mass fluxes through the dual faces from the primal cell density and the primal faces mass fluxes, which ensures a discrete mass balance on dual cells (see Section 3). We address in this paper, two different time discretizations: a purely implicit one, a pressure correction scheme. Both are unconditionally stable, and, associated to the previously described space discretization, they both provide discrete energy and pressure estimates that allow us to extract converging subsequences as ε → 0 and recover a discrete counterpart of (1.5). Although the pressure correction scheme is the only implemented scheme in practice, we present an analysis of the fully implicit scheme to emphasize the correspondence between the estimates associated to the discrete equations and the ones of the continuous problem.

The paper is organized as follows. First, we summarize in Section 2 the main ideas of the analysis of (1.1) and the singular limit ε → 0 towards (1.5). Next, in Section 3 we introduce the tools and notions associated to the spatial (staggered) discretization that is common to both schemes. The two schemes are then presented and analyzed separately in Section 4 (implicit scheme) and Section 5 (pressure correction scheme). Numerical 1D and 2D simulations associated to the pressure correction scheme are presented in the final Section 6.

Continuous setting

The aim of the section is to recall the theoretical results obtained in [START_REF] Perrin | Free/congested two-phase model from weak solutions to multidimensional compressible navier-stokes equations[END_REF] concerning the existence of solutions to (1.1) and their asymptotic behavior as ε → 0. We also provide the main tools used in the proof of these results, tools that will be adapted in the next sections to the discrete setting.

First of all, let us introduce the free energy H ε associated to the pressure p ε which is such that

H ′ ε (r)r -H ε (r) = p ε (r). (2.1)
Our choice of pressure (1.3), with γ > 1, ensures that the function r → H ε (r) is convex. We assume that f ∈ L 2 (R + × Ω). Initially, we assume that (ρ 0 ε ,

q 0 ε = ρ 0 ε u 0 ε ) is such that: 0 ≤ ρ ε 0 < 1, m(ρ 0 ε ) := 1 |Ω| Ω ρ ε 0 (x) dx ≤ ρ ⋆ < 1 a.e. in Ω, (2.2)
q 0 ε X {ρ 0 ε =0} = 0 a.e. in Ω, Ω |q ε 0 | 2 2ρ ε 0 + H ε (ρ ε 0 ) ≤ E 0 , (2.3)
for some ρ ⋆ , E 0 > 0 independent of ε.

Theorem 2.1 (Perrin & Zatorska [START_REF] Perrin | Free/congested two-phase model from weak solutions to multidimensional compressible navier-stokes equations[END_REF]). Assume the pressure law (1.3) with β > 3. Assume initially (2.2)-(2.3).

• For any ε > 0, there exists a global weak solution (ρ ε , u ε ) to (1.1)-(1.2) with finite energy;

• For ε → 0, there exists a subsequence (ρ ε , u ε , p ε (ρ ε )) ε converging to (ρ, u, p) a global weak solution to (1.5)-(1.2). More precisely:

ρ ε → ρ strongly in L q (0, T ) × Ω for all 1 ≤ q < +∞, u ε → u weakly in L 2 (0, T ; H 1 0 (Ω)), p ε (ρ ε ) → p weakly in M + (0, T ) × Ω .
Remark 2.1. As shown by Feireisl et al. in [START_REF] Feireisl | On pde analysis of flows of quasi-incompressible fluids[END_REF], the condition β > 3 can be relaxed to β > 5 2 . This condition is required for the construction of (ρ ε , u ε ) at ε fixed, but not used for the singular limit ε → 0. We will see in the next sections that no such condition, except β > 1, is needed in the discrete setting. We would expect to need this assumption if we wanted show the convergence of the discrete numerical solutions towards global weak solutions of (1.1) as the mesh size tends to 0, but this convergence is beyond the scope of this study. Sketch of the proof. For the construction of global weak solutions (ρ ε , u ε ) at ε > 0 fixed, the first step is the definition of an approximate, non-singular, problem. For that purpose, we introduce a cut-off parameter δ ∈ (0, 1) and the approximate pressure p ε,δ such that

p ε,δ (ρ) = p ε (ρ) if 0 ≤ ρ ≤ 1 -δ, εδ -β ρ γ if ρ > 1 -δ.
(2.4)

At δ > 0 fixed, the existence of global weak solutions (ρ ε,δ , u ε,δ ) is ensured by the classical theory developed by Lions and Feireisl (see for instance the book [START_REF] Novotný | Introduction to the mathematical theory of compressible flow[END_REF]). The solutions satisfy the energy inequality:

sup t∈[0,T ] Ω ρ ε,δ |u ε,δ | 2 2 + H ε,δ (ρ ε,δ ) dx + T 0 Ω µ 2 |∇u ε,δ | 2 + (λ + µ)(divu ε,δ ) 2 dx dt ≤ C 1 + f 2 L 2 ([0,T ]×Ω) ,
for some positive constant C depending on the initial datum and the viscosity, but independent of ε and δ. In particular we control the internal energy H ε,δ associated with the truncated pressure p ε,δ :

H ε,δ (ρ ε,δ ) L ∞ ([0,T ];L 1 (Ω)) ≤ C for some C > 0 independent of ε, δ.
Classical arguments based on the weak compactness of the effective flux enable to pass to the limit δ → 0 in the (weak formulation of) equations by ensuring the strong convergence of the sequence (ρ ε,δ ) δ>0 towards a limit ρ ε . The sequence (p ε,δ (ρ ε,δ )) δ>0 is proved to be equi-integrable under the condition β > 3 and it converges (strongly) in L 1 ((0, T ) × Ω) towards p ε (ρ ε ). Finally, coming back to the control of the energy

H ε,δ , one proves that meas {ρ ε,δ ≥ 1 -δ} ≤ C(ε)δ β-1 .
As a consequence of this inequality, one ensures the maximal density constraint on the limit ρ ε : 0 ≤ ρ ε < 1 a.e. in (0, T )×Ω. The second step of the proof consists in passing to the limit with respect to the parameter ε. Thanks to the energy inequality:

sup t∈[0,T ] Ω ρ ε |u ε | 2 2 + H ε (ρ ε ) dx + T 0 Ω µ 2 |∇u ε | 2 + (λ + µ)(divu ε ) 2 dx dt ≤ C 1 + f 2 L 2 ([0,T ]×Ω) , (2.5)
we control (u ε ) ε>0 in L 2 (0, T ; H 1 (Ω)), while the sequence of densities (ρ ε ) ε>0 is automatically controlled in L ∞ (0, T ) × Ω thanks to the maximal density constraint. The control of the pressure p ε (ρ ε ) follows from the use of the Bogovskii operator which is basically an inverse of the divergence operator (see [START_REF] Novotný | Introduction to the mathematical theory of compressible flow[END_REF] Chapter 3) Using

v ε = B(ρ ε -m(ρ ε )) ∈ W 1,q 0 (Ω), 1 ≤ q < +∞
, as a test function in the weak formulation of (1.1b), one gets

- T 0 Ω ρ ε u ε • B(∂ t ρ ε ) dx dt - T 0 Ω ρ ε u ε ⊗ u ε : ∇v ε dx dt + µ T 0 Ω ∇u ε : ∇v ε dx dt + (λ + µ) T 0 Ω div u ε div v ε dx dt + T 0 Ω f • v ε dx dt = T 0 Ω p ε (ρ ε ) ρ ε -m(ρ ε ) dx dt.
Thanks to the continuity properties of the Bogovskii operator ( [START_REF] Novotný | Introduction to the mathematical theory of compressible flow[END_REF], Lemma 3.17), the energy inequality and the control of the density ensure that the left-hand side of the above equality is bounded and therefore:

T 0 Ω p ε (ρ ε ) ρ ε -m(ρ ε ) dx dt ≤ C for some C > 0 independent of ε.
Splitting this integral into two parts according to the value of ρ ε and using assumption (2.2), we have: and(ρ, u, p) is a solution in the sense of distributions of the free-congested Navier-Stokes equations (1.5a)-(1.5b). Then, one can show the strong convergence of (ρ ε ) ε>0 in L q (0, T ) × Ω ∀ 1 ≤ q < +∞ using the concepts of effective flux and renormalized solutions. We do not detail further this point which can be found in classical books [START_REF] Novotný | Introduction to the mathematical theory of compressible flow[END_REF], [START_REF] Feireisl | Convergence of a mixed finite element-finite volume scheme for the isentropic navier-stokes system via dissipative measure-valued solutions[END_REF]. Eventually, the justification of the exclusion relation (1 -ρ)p = 0 is achieved in two steps. On the one hand, the strong convergence of (ρ ε ) ε>0 yields:

C ≥ T 0 Ω p ε (ρ ε ) ρ ε -m(ρ ε ) X {ρε≥(1+ρ ⋆ )/2} dx dt ≥ 1 + ρ ⋆ 2 -m(ρ 0 ε ) T 0 Ω p ε (ρ ε ) X {ρε≥(1+ρ ⋆ )/2} dx dt ≥ 1 -ρ ⋆ 2 T 0 Ω p ε (ρ ε ) X {ρε≥(1+ρ ⋆ )/2} dx dt. This achieves to show that (p ε (ρ ε )) ε is controlled in L 1 (0, T )×Ω . There exist ρ, u and p such that ρ ε → ρ weakly-* in L ∞ ((0, T ) × Ω), u ε → u weakly in L 2 (0, T ; H 1 0 (Ω)), p ε (ρ ε ) → p weakly in M + (0, T ) × Ω ,
(1 -ρ ε )p ε (ρ ε ) -→ (1 -ρ)p in D ′
(see [START_REF] Perrin | Free/congested two-phase model from weak solutions to multidimensional compressible navier-stokes equations[END_REF] for the sense given to the product ρp). On the other hand, one has

(1 -ρ ε )p ε (ρ ε ) = ε 1 β ρ α β ε (p ε (ρ ε )) β-1 β → 0 strongly in L β β-1 ((0, T ) × Ω), Hence, the limit (ρ, p) satisfies 0 ≤ ρ ≤ 1 a.e., (1 -ρ)p = 0 in D ′ , p ∈ M + ((0, T ) × Ω).

Discrete setting

Meshes and unknowns Let Ω be an open bounded connected subset of R d , d = 2 or 3. We assume that Ω is polygonal if d = 2 and polyhedral if d = 3.

Definition 3.1 (Staggered mesh).

A staggered discretization of Ω, denoted by D, is given by a pair D = (M, E), where: • M, the primal mesh, is a finite family composed of non empty simplices forming a partition of Ω : Ω = ∪ K∈M K. We denote by E the set of faces of the mesh, and we suppose that two neighboring cells share a whole face: for all σ ∈ E, either σ ⊂ ∂Ω or there exists (K, L) ∈ M 2 with K = L such that K ∩ L = σ; we denote in the latter case σ = K|L. We denote

n K,σ K L M D σ D σ ′ σ = K | L σ ′ = K | M ǫ = D σ |D σ ′
E ext = {σ ∈ E, σ ⊂ ∂Ω} and E int = E \ E ext .
For K ∈ M, E(K) stands for the set of faces of K. The unit vector normal to σ ∈ E(K) outward K is denoted by n K,σ .

• We define a dual mesh associated with the faces σ ∈ E as follows. For K ∈ M and σ ∈ E(K), we define D K,σ as the cone with basis σ and with vertex the mass center of K (see Figure 1). We 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 Eventually, we define the size of the discretization as h M = max

|D K,σ | = |K|/m. For σ ∈ E int , σ = K|L, we now define the diamond cell D σ associated with σ by D σ = D K,σ ∪ D L,σ . For σ ∈ E ext ∩ E(K),
K∈M h K .
Relying on Definition 3.1, we now define a staggered space discretization.

Definition 3.2. Let D = (M, E) be a staggered discretization of Ω. We denote L M (Ω) the space of scalar functions that are constant on each primal cell K ∈ M. For ρ ∈ L M (Ω) and K ∈ M, we denote ρ K the constant value of ρ on K. We denote L M,0 (Ω) the subspace of L M (Ω) composed of zero average functions over Ω. We denote H M (Ω) the space of functions u such that u |K ∈ P 1 (K) for all K ∈ M and such that:

1 |σ| σ [u] σ dσ(x) = 0, ∀σ ∈ E int , (3.1)
where

[u] σ is the jump of u through σ which is defined on σ = K|L by [u] σ = u |L -u |K . We define H M,0 (Ω) ⊂ H M (Ω)
the subspace of H M (Ω) composed of functions the degrees of freedom of which are zero over ∂Ω, i.e. the functions u ∈ H M (Ω) such that 1 |σ| σ u dσ(x) = 0 for all σ ∈ E ext . Finally, we denote H M (Ω) := H M (Ω) d and H M,0 (Ω) := H M,0 (Ω) d . For a discrete velocity field u ∈ H M (Ω) and σ ∈ E, the degree of freedom associated with σ is given by:

u σ = 1 |σ| σ u dσ(x). (3.2)

Discrete operators

Mass convection operator Given discrete density and velocity fields ρ ∈ L M (Ω) and u ∈ H M,0 (Ω), the discretization of the mass convection term is given by:

div M (ρu)(x) = K∈M 1 |K| σ∈E(K) F K,σ (ρ, u) X K (x), (3.3)
where X K is the characteristic function of the subset K of Ω. The quantity F K,σ (ρ, u) 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,σ (ρ, u) = |σ| ρ σ u σ • n K,σ , ∀σ ∈ E int , σ = K|L. (3.4)
The density at the face σ = K|L is approximated by the upwind technique, i.e.

ρ σ = ρ K if u σ • n K,σ ≥ 0 ρ L otherwise. (3.5)
Velocity convection operator We now describe the approximation of the convection operator ∂ t (ρu)+ div(ρu ⊗ u). The approximation of the time derivative part ∂ t (ρu) is naturally discretized at the dual cells D σ , σ ∈ E int and an approximation ρ Dσ of the density on these dual cells D σ is thus needed. Given a density field ρ ∈ L M (Ω), this approximation is built as follows:

|D σ | ρ Dσ = |D K,σ | ρ K + |D L,σ | ρ L , ∀σ ∈ E int , σ = K|L. (3.6)
Given discrete a density ρ ∈ L M (Ω) and velocity fields u, v ∈ H M,0 (Ω), the discretization of the convection term is given by:

div E (ρu ⊗ v)(x) = σ∈Eint 1 |D σ | ǫ∈ Ẽ(Dσ ) F σ,ǫ (ρ, u) v ǫ X Dσ (x).
(3.7)

F σ,ǫ (ρ, u) is the mass flux across the edge ǫ of the dual cell D σ . Its value is zero if ǫ ∈ Ẽext . Otherwise, it is defined as a linear combination, with constant coefficients, of the primal mass fluxes at the neighboring faces. For K ∈ M and σ ∈ E(K), let ξ σ K be given by:

ξ σ K = |D K,σ | |K| = 1 d + 1 , (3.8)
so that σ∈E(K) ξ σ K = 1. The dual mass fluxes F σ,ǫ (ρ, u) are then computed to as to satisfy the following three conditions: (H1) For all primal cell K in M, the set (F σ,ǫ (ρ, u)) ǫ⊂K of dual fluxes included in K solves the following linear system

F K,σ (ρ, u) + ǫ∈ Ẽ(Dσ ), ǫ⊂K F σ,ǫ (ρ, u) = ξ σ K σ ′ ∈E(K) F K,σ ′ (ρ, u), σ ∈ E(K). (3.9)
(H2) The dual fluxes are conservative, i.e. for any dual face

ǫ = D σ |D σ ′ , we have F σ,ǫ (ρ, u) = -F σ ′ ,ǫ (ρ, u).
(H3) The dual fluxes are bounded with respect to the primal fluxes (F K,σ (ρ, u)) σ∈E(K) , in the sense that

|F σ,ǫ (ρ, u)| ≤ max {|F K,σ ′ (ρ, u)|, σ ′ ∈ E(K)} , (3.10) for K ∈ M, σ ∈ E(K), ǫ ∈ Ẽ(D σ ) with ǫ ⊂ K.
This convection operator is built so that a discrete mass conservation equation can be shown to be satisfied on the cells of the dual mesh:

|D σ | δt (ρ n+1 Dσ -ρ n Dσ ) + ǫ∈ Ẽ(Dσ ) F σ,ǫ (ρ n+1 , u n+1 ) = 0, ∀σ ∈ E int . (3.11)
To complete the definition of the convective flux, we must give the expression of the velocity v ǫ at the dual face:

v ǫ = 1 2 (v σ + v σ ′ ), for ǫ = D σ |D ′ σ .
Diffusion operator -Let us define the shape functions associated with the Crouzeix-Raviart finite element. These are the functions (ζ σ ) σ∈E where for all σ ∈ E, ζ σ is the element of H M (Ω) which satisfies:

1 |σ ′ | σ ′ ζ σ dσ(x) = 1, if σ ′ = σ, 0, if σ ′ = σ.
Given a discrete velocity field u ∈ H M,0 (Ω), the discretization of the diffusion terms is given by:

div E (τ (u)) = µ∆ E u + (µ + λ)(∇ • div) E u with ∆ E u(x) = σ∈Eint 1 |D σ | K∈M K ∇u .∇ζ σ dx X Dσ (x), (∇ • div) E u(x) = σ∈Eint 1 |D σ | K∈M K div u∇ζ σ dx X Dσ (x).
(3.12)

Truncated pressure and gradient operator -Given a discrete pressure field p ∈ L M (Ω), the pressure gradient term is discretized as follows:

∇ E p(x) = σ∈Eint σ=K|L |σ| |D σ | (p L -p K )) n K,σ X Dσ (x). (3.13)
We also introduce a truncated pressure p ε,M defined as follows:

p ε,M (r) =      ε r γ (1 -r) β if r ≤ 1 -δ ε,M , h -αβ M r γ if r > 1 -δ ε,M , where δ ε,M = h α M ε 1 β-1 and α(β -1) > 3. (3.14)
As said in the introduction, the cut-off of the pressure is a technical detail that ensures the existence of a numerical solution to the schemes. In practice, there is no problem to keep the original pressure in the software CALIF 3 S. We will also see later in the analysis that the cut-off density 1 -δ ε,M is never reached provided the initial density is less than 1 -δ ε,M .

Discrete norms and properties

We start with defining the piecewise smooth first order differential operators associated with the Crouzeix-Raviart non-conforming finite element representation of velocities u ∈ H M (Ω) :

∇ M u(x) = K∈M ∇u(x)X K (x), div M u(x) = K∈M div u(x)X K (x). (3.15)
We define an interpolation operator Π E which associates a piecewise constant function over the cells of the dual mesh to any function u ∈ H M (Ω) as follows:

Π E u(x) = σ∈Eint u σ X Dσ (x).
(3.16)

The mapping u → Π E u is a one-to-one mapping which is continuous with respect to the L q -norm, for all q ∈ [1, +∞]: there exists a constant C = C(q, d, D) such that:

Π E u L q (Ω) ≤ C u L q (Ω) .
Lemma 3.1. Let D = (M, E) be a staggered discretization of Ω. The following discrete integration by parts formulas are satisfied for all (p, u) ∈ L M (Ω) × H M,0 (Ω). One has for all v ∈ H M,0 (Ω):

-

Ω ∆ E u • Π E v dx = Ω ∇ M u : ∇ M v dx, - Ω (∇ • div) E u • Π E v dx = Ω div M u div M v dx, Ω ∇ E (p) • Π E v dx = - Ω p div M v dx.
(3.17)

We then define for q ∈ [1, ∞) the broken Sobolev W 1,q semi-norm . 1,q,M associated with the Crouzeix-Raviart finite element representation of the discrete velocities. For any u ∈ H M (Ω) it is given by:

u q 1,q,M = Ω |∇ M u| q dx.
An important property is the following discrete Sobolev embedding: for all q ∈ [1, +∞) if d = 2 and for all q ∈ [1, 6] if d = 3, there exists C = C(q, d, D) > 0 such that:

u L q (Ω) ≤ C u 1,2,M , ∀u ∈ H M,0 (Ω).
Note that the semi-norm u 1,2,M is in fact a norm on the space H M,0 (Ω).

Presentation of the scheme

Let δt > 0 be a constant time step. The approximate solution (ρ n , u n ) ∈ L M (Ω) × H M,0 (Ω) at time t n = nδt for 1 ≤ n ≤ N = ⌊T /δt⌋ is computed by induction through the following implicit scheme.

Knowing (ρ n , u n ) ∈ L M (Ω) × H M,0 (Ω), solve for ρ n+1 ∈ L M (Ω) and u n+1 ∈ H M,0 (Ω): 1 δt (ρ n+1 K -ρ n K ) + div(ρ n+1 u n+1 ) K = 0, ∀K ∈ M, 1 δt ρ n+1 Dσ u n+1 σ -ρ n Dσ u n σ + div(ρ n+1 u n+1 ⊗ u n+1 ) σ -div(τ (u n+1 )) σ + (∇p n+1 ε,M ) σ = f n+1 σ , ∀σ ∈ E int (4.1a) (4.1b)
As in the continuous case (see (2.1)), we can define an internal energy H ε,M ∈ C 1 (R + ) associated to this truncated pressure, it is such that:

H ′ ε,M (r)r -H ε,M (r) = p ε,M (r) ∀r ∈ [0, 1). (4.2)
Written differently,

H ′′ ε,M (r) = p ′ ε,M (r)/r for any r = 1 -δ ε,M . Since p ′ ε,M (r) > 0 for all r = 1 -δ ε,M , H ε,M is convex.
Initialization of the scheme As in Section 2, we shall assume that the continuous initial datum is such that ρ ε 0 ∈ (0, 1) a.e. and

ρ ε 0 u ε 0 L 2 (Ω) + H ε (ρ ε 0 ) L 1 (Ω) ≤ C 0 , (4.3)
for some C 0 independent of ε and the mesh. The initial approximations are given by the average of the initial density ρ ε 0 on the primal cells and the initial velocity u ε 0 on the dual cells:

ρ 0 K = 1 |K| K ρ ε 0 (x) dx, ∀ K ∈ M, u 0 σ = 1 |D σ | Dσ u ε 0 (x) dx, ∀ σ ∈ E int . (4.4)
We also define

f n σ = 1 δt 1 |D σ |

Estimates

We present in this subsection how to derive uniform estimates w.r.t the parameter ε. To lighten the presentation, we will just detail the estimates that are specific to our singular pressure case, other estimates that are identical to the classical Navier-Stokes equations and the low Mach limit are just stated and we refer the reader to [START_REF] Herbin | Low mach number limit of some staggered schemes for compressible barotropic flows[END_REF] for their proof.

The first result is a reformulation of the discrete mass balance (4.1a). It is simply obtained by multiplying (4.1a) H ′ ε,M (ρ n+1 K ) and by using the definition of the flux (3.4). Lemma 4.1 (Discrete renormalization identity). Any solution to the numerical implicit scheme (4.1) satisfies for all K ∈ M and 0 ≤ n ≤ N -1 the following identity

1 δt H ε,M (ρ n+1 K ) -H ε,M (ρ n K ) + div(H ε,M (ρ n+1 )u n+1 ) K + R n+1 K = -p ε,M (ρ n+1 K ) div(u n+1 ) K (4.6)
with:

R n+1 K = 1 δt H ε,M (ρ n K ) -H ε,M (ρ n+1 K ) -H ′ ε,M (ρ n+1 K )(ρ n K -ρ n+1 K ) + 1 |K| σ∈E(K) |σ|r n+1 K,σ u n+1 σ • n K,σ , r n+1 K,σ = (H ε,M ) ′ (ρ n+1 K )(ρ n+1 σ -ρ n+1 K ) + H ε,M (ρ n+1 K ) -H ε,M (ρ n+1 σ
).

(4.7)

Moreover, since H ε,M is convex we have 0 ≤ K∈M |K|R n+1 K .

We now are able to state a discrete kinetic energy balance obtained by multiplying the momentum equation (4.1b) with the discrete velocity u n+1 σ . Proposition 4.2 (Discrete kinetic energy balance).

1 2δt ρ n+1 Dσ |u n+1 σ | 2 -ρ n Dσ |u n σ | 2 + 1 2|D σ | ǫ=Dσ|D σ ′ F σ,ǫ (ρ n+1 , u n+1 )u n+1 σ • u n+1 σ ′ -div(τ (u n+1 )) σ • u n+1 σ + (∇p n+1 ε ) σ • u n+1 σ + R n+1 σ = f n+1 σ • u n+1 σ , (4.8) where R n+1 σ = 1 2δt ρ n Dσ |u n+1 σ -u n σ | 2 .
Proof. Let us take the scalar product of the momentum equation (4.1b) with u n+1 σ , we get:

1 δt ρ n+1 Dσ u n+1 σ -ρ n Dσ u n σ + div(ρ n+1 u n+1 ⊗ u n+1 ) σ • u n+1 σ -div(τ (u n+1 )) σ • u n+1 σ + (∇p n+1 ε ) σ • u n+1 σ = f n+1 σ • u n+1 σ
Let us rewrite the first term as follows:

1 δt ρ n+1 Dσ u n+1 σ -ρ n Dσ u n σ + div(ρ n+1 u n+1 ⊗ u n+1 ) σ • u n+1 σ =      1 δt ρ n+1 Dσ u n+1 σ -ρ n Dσ u n σ + 1 2|D σ | ǫ∈ Ẽ(Dσ ) ǫ=Dσ|D σ ′ F σ,ǫ (ρ n+1 , u n+1 ) (u n+1 σ + u n+1 σ ′ )      • u n+1 σ = 1 2δt ρ n+1 Dσ |u n+1 σ | 2 -ρ n Dσ |u n σ | 2 + 1 2|D σ | ǫ=Dσ|D σ ′ F σ,ǫ (ρ n+1 , u n+1 ) u n+1 σ • u n+1 σ ′ + 1 2δt ρ n Dσ |u n+1 σ -u n σ | 2 +   1 δt (ρ n+1 Dσ -ρ n Dσ ) + |D σ | -1 ǫ∈E(Dσ) F σ,ǫ (ρ n+1 , u n+1 )   |u n+1 σ | 2 2
where we have used the identity 2 (ρ|a|

2 -ρ * a • b) = ρ|a| 2 -ρ * |b| 2 + ρ * |a -b| 2 + (ρ -ρ * )|a| 2 for ρ = ρ n+1 Dσ , ρ * = ρ n Dσ , a = u n+1 σ and b = u n σ .
Eventually, we observe that the last term is equal to 0 thanks to the discrete mass conservation satisfied on the dual cells (3.11) Next, we get a total (local in time) energy inequality by combining discrete kinetic energy balance (4.8) over the faces σ ∈ E int and the discrete renormalized equation (4.6) summed over K ∈ M. We also prove the existence of a discrete solution to (4.1). Lemma 4.3 (Local-in-time energy inequality, existence of the numerical solutions). Let ε > 0 and assume that the initial data satisfy (4.3). Then, there exists a solution (ρ n ε , u n ε ) 0≤n≤N to the scheme (4.1), such that ρ n ε > 0 for 0 ≤ n ≤ N , and the following inequality holds for 0 ≤ n ≤ N -1:

1 2 σ∈Eint |D σ | ρ n+1 Dσ |u n+1 σ | 2 -ρ n Dσ |u n σ | 2 + K∈M |K| H ε (ρ n+1 K ) -H ε (ρ n K ) + µ 2 δt u n+1 2 1,2,M + R n+1 ≤ Cδt f n+1 2 L 2 (Ω) , (4.9)
where

R n+1 = σ∈Eint |D σ |R n+1 σ + K∈M |K|R n+1 K ≥ 0.
Sketch of proof. We sum the discrete kinetic energy balance (4.8) over the faces σ ∈ E int and add the discrete renormalized equation (4.6) summed over K ∈ M. The desired inequality is then obtained as in [START_REF] Herbin | Low mach number limit of some staggered schemes for compressible barotropic flows[END_REF] thanks to the discrete grad-div duality (3.19), the conservativity of the dual fluxes (H2), and the coercivity of the diffusion operator. The positivity of the density is a consequence of the properties of the upwind choice (3.4) for ρ [11, Lemma 2.1]. Given discrete density and velocity fields (ρ n , u n ), the proof of existence of a solution (ρ n+1 , u n+1 ) to the implicit scheme at the time step n + 1 is the same as that given in [START_REF] Herbin | correction staggered schemes for barotropic one-phase and two-phase flows[END_REF]Proposition 3.3] in the case of the Euler equations. It may be inferred by the Brouwer fixed point theorem, by an easy adaptation of the proof of [START_REF] Eymard | Convergence of the MAC Scheme for the Compressible Stokes Equations[END_REF]Proposition 5.2]. 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 L 1 -norm. Then, by a norm equivalence argument, the truncated pressure p ε,M is controlled in any norm. Next, for a given density ρ, the discrete global kinetic energy inequality (i.e. (4.8) summed over the control volumes) provides a control on the velocity. Therefore, computing ρ from the mass balance for fixed u, then p ε,M from ρ by the equation of state and finally u from the momentum balance equation with fixed ρ, yields an iteration in a bounded convex subset of a finite dimensional space. Remark 4.1. It is important to note that the truncation of the pressure is necessary in the argument above to infer a control of the pressure from the control of the density ρ. The control of the discrete truncated pressure is then obviously not uniform with respect to ε, but this is not an issue here since we just want to ensure at this stage the existence of a solution. Suitable bounds are derived below, see in particular Lemma 4.6.

Integrating over time the previous inequality, we get the final energy inequality: Lemma 4.4 (Global discrete energy inequality). Let ε > 0 and assume that the initial data (ρ ε 0 , u ε 0 ) has finite energy. There exists C > 0 independent of ε and of the mesh such that, for ε small enough and for all 1 ≤ n ≤ N :

1 2 σ∈Eint |D σ | ρ n Dσ |u n σ | 2 + µ 2 n k=1 δt ũk 2 1,2,M + K∈M |K| H ε (ρ n K ) ≤ C. (4.10)
Proof. Summing (4.9) over n yields the expected inequality (4.10

) with C = 1 2 σ∈Eint |D σ | ρ 0 Dσ |u 0 σ | 2 + K∈M |K| H ε (ρ 0 K ) + f 2 L 2
, which is bounded uniformly with respect to ε by hypothesis (4.3).

The next lemma gives the crucial upper bound on the density.

Lemma 4.5 (Maximal constraint on the density). Assume that α(β -1) > 3. There exists a constant C M > 0 independent of ε (but not on the mesh), such that

ρ n K ≤ 1 -C M ε 1 β-1 ∀ K ∈ M, n ∈ N. (4.11)
Proof. The bound follows from the control of the energy H ε,M (ρ n ) provided by the discrete energy estimate (4.10). Indeed, let us rewrite H ε,M as

H ε,M (ρ) = ρe ε,M (ρ) with e ε,M (ρ) = ρ 0 p ε,M (s) s 2 ds.
To simplify, let us assume that γ is an integer. We have

Ω H ε,M (ρ n ε ) dx ≥ Ω ρ n ε e ε,M (ρ n ε )X {ρ n ε >1-δε,M} dx ≥ Ω ρ n ε e ε,M (1 -δ ε,M )X {ρ n ε >1-δε,M} dx = ε Ω ρ n ε 1-δ 0 s γ-2 (1 -s) β ds X {ρ n ε >1-δε,M} dx ≥ ε Ω ρ n ε 1 δ (1 -s) γ-2 s β ds X {ρ n ε >1-δε,M} dx = ε Ω ρ n ε γ-2 k=0 (-1) k γ -2 k 1 δ s k-β ds X {ρ n ε >1-δε,M} dx.
Hence,

Ω H ε,M (ρ n ε ) dx ≥ δ 1-β -1 β -1 ε Ω ρ n ε X {ρ n ε >1-δε,M} dx + ε Ω ρ n ε γ-2 k=1 (-1) k γ -2 k 1 δ s k-β ds X {ρ n ε >1-δε,M} dx ≥ C 1 (γ, β) ε δ 1-β ε,M meas{ρ n ε > 1 -δ ε,M } -C 2 (γ, β)
For some positive constants C 1,2 independent of ε. As a result of the global energy inequality (4.10), we deduce that meas

{ρ n ε > 1 -δ ε,M } ≤ Cε -1 δ β-1 ε,M ≤ Ch α(β-1)
for some constant C independent of h and ε. Since we have assumed α(β -1) > 3, we ensure that meas

{ρ n ε > 1 -δ ε,M } < |K| for all K ∈ M or, in other words, ρ n K ≤ 1 -δ ε,M = 1 -h α ε 1 β-1 for all K ∈ M.
The previous bound on the density ensures that ρ n K < 1. Unfortunately this bound does not yield any uniform (with respect to ε) control of the pressure since

p ε,M (1 -δ ε,M ) ≤ Cε -1 β-1 .
Proposition 4.6 (Control of the pressure). Assume initially (4.3) as well as (2.2). For 1 ≤ n ≤ N the following inequality holds:

p n ε,M L 1 (Ω) ≤ C M , (4.13)
where the real number C M > 0 depends on the mesh size but not on ε.

The proof of the lemma relies on the following result which ensures the existence of a discrete equivalent of the Bogovskii operator (cf Section 2). Lemma 4.7 (L q Discrete inf-sup property). Let D = (M, E) be a staggered discretization of Ω. Then, there exists a linear operator B M : L M,0 (Ω) -→ H M,0 (Ω), depending only on Ω and on the discretization such that the following properties hold:

(i) For all p ∈ L M,0 (Ω), Ω r div M (B M p) dx = Ω r p dx, ∀r ∈ L M (Ω).
(ii) For all q ∈ (1, +∞), there exists C = C(q, d, Ω, D), such that

B M p 1,q,M ≤ C p L q (Ω) .
Proof of Proposition 4.6. Let us set v n+1 = B M (ρ n+1 -m(ρ n+1 )) ∈ H M,0 (Ω) where m(ρ n+1 ) denotes the average of ρ n+1 . We multiply the discrete momentum equation (4.1b) by v n+1 σ , multiply by |D σ | and sum over σ ∈ E int to get:

- σ∈Eint |D σ |(∇p ε,M ) n+1 σ v n+1 σ = 1 δt σ∈Eint |D σ |(ρ n+1 Dσ u n+1 σ -ρ n Dσ u n σ )v n+1 σ - σ∈Eint |D σ |div(τ (u n+1 )) σ • v n+1 σ .
Thanks to the bounds previously derived on ρ and u, we bound the right-hand side uniformly with respect to ε. On the left-hand side, thanks to (3.19) we have by Lemma 3.1

- σ∈Eint |D σ |(∇p ε,M ) n+1 σ v n+1 σ = Ω p n+1 ε,M div M v n+1 dx,
and therefore

Ω p n+1 ε,M div M v n+1 dx = Ω p n+1 ε,M (ρ n+1 -m(ρ n+1 )) dx ≤ C M .
Now, we can split the pressure integral into two parts:

Ω p n+1 ε,M • (ρ n+1 -m(ρ n+1 )) dx = Ω p n+1 ε,M • (ρ n+1 -m(ρ n+1 ))X {ρ n+1 > 1 2 (1+ρ * )} dx + Ω p n+1 ε,M • (ρ n+1 -m(ρ n+1 ))X {ρ n+1 ≤ 1 2 (1+ρ * )} dx.
In the domain {ρ n+1 ≤ 1 2 (1+ρ * )}, the pressure is uniformly bounded with respect to ε so the corresponding integral is controlled uniformly with respect to ε. Hence, using the assumption (2.2):

C M ≥ Ω p n+1 ε,M • (ρ n+1 -m(ρ n+1 ))X {ρ n+1 > 1 2 (1+ρ * )} dx ≥ 1 + ρ * -2m(ρ n+1 ) 2 Ω p n+1 ε,M X {ρ n+1 > 1 2 (1+ρ * )} dx ≥ 1 -ρ * 2 Ω p n+1 ε,M X {ρ n+1 > 1 2 (1+ρ * )} dx
which yields the control of the pressure in L 1 (Ω).

Limit ε → 0

Thanks to the previous discrete estimates, we infer that the sequences (ρ n ε,M ) ε<ε0 , (p ε (ρ n ε,M )) ε<ε0 and (u n ε,M ) ε<ε0 are bounded respectively in L ∞ (Ω), L 1 (Ω) and H M,0 (Ω). Since all the discrete norms are equivalent in finite dimension, we deduce that the previous sequences are bounded in any discrete norm. We infer that there exist limits ρ, p, u such that the following convergences (up to subsequences) hold

ρ n ε,M → ρ n , p ε (ρ n ε,M ) → p n , u n ε,M → u n ,
in any discrete norm. As a consequence of (4.11), it is easy to see that we ensure the maximal density constraint on the limit density

0 ≤ ρ n K ≤ 1 ∀ K ∈ M, n = 0, . . . , N. (4.14)
We also have p n ≥ 0. Moreover, since on each cell we have

1 -(ρ n ε,M ) K p ε,M ((ρ n ε,M ) K ) = ε (ρ n ε,M ) γ K (1 -(ρ n ε,M ) K ) β-1 = ε 1 β (ρ n ε,M ) γ β K p ε ((ρ n ε,M ) K ) β-1 β
and we get then as ε → 0

1 -ρ n K p n K = 0 where p n K = lim ε→0 p ε (ρ n ε,M ) K .
Remark 4.2. In the continuous case, it was proved in [START_REF] Lions | On a free boundary barotropic model[END_REF] that divu = 0 in {ρ = 1}, i.e. the dynamics of the congested domain is incompressible. We could a priori hope a discrete equivalent of this incompressibility condition but the only (trivial) observation that we can make is the following: assume that on the cell K, ρ n K , ρ n+1 K , ρ n+1 σ are all equal to 1, then by (4.1a) (divu) K = 0 (see [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equations[END_REF] where a similar observation is made).

Pressure correction scheme

In practice, the implementation of the previous implicit scheme (4.1) requires to find the solution of a fully non-linear coupled system. This task is difficult in a real computational context due to the computational cost and lack of robustness. For that reason, we also perform the analysis of the limit ε → 0 for the pressure correction scheme implemented in CALIF 3 S [2]. The scheme is obtained thanks to a partial decoupling of the discrete equations. It consists (after a rescaling step for the pressure gradient) in first 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, and then 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. The scheme reads as follows:

Knowing (ρ n-1 , ρ n , u n ) ∈ L M (Ω) × L M (Ω) × H M,0 (
Ω), we compute (ρ n+1 , u n+1 ) as follows:

Pressure gradient scaling step:

(∇p ε ) n σ = ρ n Dσ ρ n-1 Dσ 1/2 (∇p ε ) n σ ∀σ ∈ E int (5.1)
Prediction step: Solve for ũn+1 ∈ H M,0 (Ω):

1 δt ρ n Dσ ũn+1 σ -ρ n-1 Dσ u n σ + 1 |D σ | ǫ∈ Ẽ(Dσ ) F σ,ǫ (ρ n , u n ) ũn+1 ǫ +(∇p ε ) n σ -divτ (ũ n+1 ) σ = f n+1 σ ∀σ ∈ E int (5.2)
Correction step: Solve for ρ n+1 ∈ L M (Ω) and u n+1 ∈ H M,0 (Ω):

1 δt ρ n+1 K -ρ n K + 1 |K| σ∈E(K) F K,σ (ρ n+1 , u n+1 ) = 0 ∀K ∈ M, 1 δt ρ n Dσ u n+1 σ -ũn+1 σ + (∇p ε ) n+1 σ -(∇p ε ) n σ = 0 ∀σ ∈ E int . (5.3a) (5.3b)
Initial data and initialization of the scheme Definition 5.1. The initial data

(ρ ε 0 , u ε 0 ) is said to be well-prepared if 0 < ρ ε 0 < 1, u ε 0 ∈ H 1 0 (Ω) d for all ε > 0 satisfy (2.2)-(2.
3) and if, in addition,

u ε 0 H 1 (Ω) d + p ε (ρ ε 0 ) L 1 (Ω) + (1 -ρ ε 0 ) -1 [div u ε 0 ] -L 2 (Ω) ≤ C. (5.4)
Remark 5.1. • Note that the initial control of the pressure requires

1 -ρ 0 ε ≥ C 0 ε 1 β for some C 0 > 0 independent of ε, (5.5)
which is actually more stringent than the control (2.3) of the energy H ε (ρ 0 ε ). • The condition on the divergence of the initial velocity ensures that the initial compression (i.e.

[divu] -) is very small in the dense regions where ρ ε 0 is close to 1. The source term f n σ is defined as in (4.5). The initial velocities u ε 0 are in H 1 0 (Ω) d . In particular, their trace is thus well defined as an L 2 function on any smooth (or Lipschitz-continuous) hypersurface of Ω. The initialization of the pressure correction scheme is performed as follows. First, ρ 0 and u 0 are given by the average of the initial conditions ρ ε 0 and u ε 0 respectively on the primal cells and on the faces of the primal cells:

ρ -1 K = 1 |K| K ρ ε 0 (x) dx ∀K ∈ M, u 0 σ = 1 |σ| σ u ε 0 (x) dσ(x) ∀σ ∈ E int . (5.6)
Finally, we compute ρ 0 by solving the mass balance equation (5.3a) for n = -1:

1 δt ρ 0 K -ρ -1 K + 1 |K| σ∈E(K) F K,σ (ρ 0 , u 0 ) = 0 ∀K ∈ M. (5.7)
Basically the steps for the derivation of the discrete energy estimate and the convergence as ε → 0 are completely analogous to the previous implicit scheme. Nevertheless, due to the initialization of scheme, we need to ensure to properties on the the first iteration ρ 0 . Lemma 5.1 (Bounds on ρ 0 ). For any ε ≤ ε 0 , there exists ρ 0 , A 0 ε > 0 such that A 0 ε ≤ 1 -Cε 1/β for some C > 0 independent of ε, and

ρ 0 ≤ ρ 0 K ≤ A 0 ε < 1 ∀ K ∈ M. (5.8)
Proof.

• Lower bound. Let us set

ρ 0 = min K∈M ρ -1 K 1 + δt max 0, max K∈M div(u 0 ) K .
and assume that min K∈M ρ 0 K = ρ 0 K < ρ 0 . Equation (5.7) rewrites as follows:

1 δt ρ 0 K -ρ -1 K + ρ 0 K (div u 0 ) K + 1 | K| σ∈E( K) |σ|(ρ 0 σ -ρ 0 K )u 0 σ • n K,σ = 0,
where, by definition of ρ 0 σ (upwind) and

ρ 0 K = min K∈M ρ 0 K we have (ρ 0 σ -ρ 0 K )u 0 σ • n K,σ ≤ 0.
As a consequence, we get the inequality

1 δt ρ 0 K -ρ -1 K + ρ 0 K (div u 0 ) K ≥ 0.
If (div u 0 ) K < 0, we obtain ρ 0 K ≥ ρ -1 K while for (div u 0 ) K ≥ 0:

ρ 0 K ≥ ρ -1 K 1 + δt(div u 0 ) K ,
In both cases, we have ρ 0 K > ρ 0 and thus a contradiction.

• Upper bound. Let K ∈ M such that max K∈M ρ 0 K = ρ 0 K . We have

ρ 0 K = ρ -1 K -δtρ 0 K div(u 0 ) K - δt | K| σ∈E( K) |σ|(ρ 0 σ -ρ 0 K )u 0 σ • n K,σ
Similarly to the previous case we have:

1 δt ρ 0 K -ρ -1 K + ρ 0 K (div u 0 ) K ≤ 0.
If (div u 0 ) K > 0, we directly get thanks to (5.5) that ρ 0 K ≤ ρ -1 K ≤ 1-C 0 ε 1/β while for (div u 0 ) K ≤ 0, we have:

ρ 0 K ≤ ρ -1 K 1 -δt(div u 0 ) K .
Now, using assumption (5.4) we have max

K [div(u 0 )] -≤ C M min K (1 -ρ -1 K ) ≤ C M ε 1/β
and therefore, for some C > 0 independent of ε,

max ρ 0 K ≤ max ρ -1 K 1 + Cδt min(1 -ρ -1 K ) ≤ 1 -Cε 1/β .
Once we have the control of ρ 0 , we can deduce with a similar argument the positivity of the density at any time step:

ρ n K > 0 ∀K ∈ M, n = 1, . . . , N.
(5.9)

Identically to the implicit scheme, we have a discrete renormalization identity (4.6). Next, we derive the kinetic energy estimate by taking the scalar product of the velocity prediction equation (5.2) with the corresponding velocity unknown ũn+1

σ : 1 2δt ρ n Dσ |u n+1 σ | 2 -ρ n-1 Dσ |u n σ | 2 + 1 2|D σ | F σ,ǫ (ρ n , u n )ũ n+1 σ • ũn+1 σ ′ -div(τ (ũ n+1 )) σ • ũn+1 σ + (∇p n+1 ε ) σ • u n+1 σ + δt 2 (∇p n+1 ε ) σ 2 ρ n Dσ - (∇p n ε ) σ 2 ρ n-1 Dσ + R n+1 σ = δtf n+1 σ • ũn+1 σ , (5.10)
where

R n+1 σ = 1 2δt ρ n-1 Dσ |ũ n+1 σ -u n σ | 2 .
Note that we used the rescaling step (5.1) to obtain this final expression (see details of the calculations in [START_REF] Herbin | Low mach number limit of some staggered schemes for compressible barotropic flows[END_REF]).

Combining the renormalized equation with (5.10), we infer the following result.

Lemma 5.2 (Local-in-time discrete energy inequality, existence of a solution). There exists a solution (ρ n , u n ) 0≤n≤N to the scheme such that ρ n > 0 for 1 ≤ n ≤ N and the following inequality holds for 0 ≤ n ≤ N -1: The existence of a solution (ρ n+1 , u n+1 ) to the correction step (5.3a)-(5.3b) follows from the Brouwer fixed point theorem, by an easy adaptation of the proof of [START_REF] Eymard | Convergence of the MAC Scheme for the Compressible Stokes Equations[END_REF]Proposition 5.2].

1 2 σ∈Eint |D σ | ρ n Dσ |u n+1 σ | 2 -ρ n-1 Dσ |u n σ | 2 + K∈M |K| H ε (ρ n+1 K ) -H ε (ρ n K ) + µ 2 δt ũn+1 2 1,2,M + δt 2 2 σ∈Eint |D σ | (∇p n+1 ε ) σ 2 ρ n Dσ - (∇p n ε ) σ 2 ρ n-1 Dσ + R n+1 ≤ Cδt f n+1 2 L 2 (Ω) , (
Let us now turn to the global entropy inequality; it was first proven in [10, Theorem 3.8], we prove here that the bound is independent of ε if the initial data is well-prepared.

Lemma 5.3 (Global discrete energy inequality). Let ε > 0 and assume that the initial data (ρ ε 0 , u ε 0 ) is well-prepared. By Lemma 5.2, there exists a solution (ρ n , u n ) n≤N to the scheme (5.2)-(5.3). In addition, there exists C > 0 independent of ε such that, for ε small enough and for all 1 ≤ n ≤ N :

1 2 σ∈Eint |D σ | ρ n-1 Dσ |u n σ | 2 + µ 2 n k=1 δt ũk 2 1,2,M + K∈M |K| H ε (ρ n K ) + δt 2 2 σ∈Eint |D σ | ρ n-1 Dσ |(∇p n ε ) σ | 2 ≤ C. (5.12)
Proof. Summing (5.11) over n yields the expected inequality (5.12) with

C = 1 2 σ∈Eint |D σ | ρ -1 Dσ |u 0 σ | 2 + K∈M |K| H ε (ρ 0 K ) + δt 2 2 σ∈Eint |D σ | ρ -1 Dσ |(∇p 0 ε ) σ | 2 + f 2 L 2 .
Let us prove that if the initial data is well-prepared in the sense of (5.4), then for ε small enough, C is uniformly bounded independently of ε. Since 0 < r 1 ≤ ρ -1 K < 1 -Cε 1/β for all K ∈ M, so is ρ -1 Dσ for all σ ∈ E int . Hence, since u ε 0 is uniformly bounded in H 1 (Ω) d by (5.4), a classical trace inequality yields the boundedness of the first term. By initial conditions , one has |ρ 0 K -1| ≥ Cε 1/β for all K ∈ M. Hence, the second term vanishes as ε → 0. The third term is also uniformly bounded with respect to ε since the initial data is supposed to be well-prepared.

Control of the density and estimate of the pressure As for the implicit case, the control of the energy (5.12) provides a control on the density close to 1: there exists a constant C > 0 independent of ε (but not on the mesh), such that

ρ n K ≤ 1 -Cε 1 β-1 ∀ K ∈ M, n ∈ N.
(5.13) Proposition 5.4 (Control of the pressure). One has, for 1 ≤ n ≤ N :

∇p n ε L 2 (Ω) + p n ε L 1 (Ω) ≤ C M,δt , (5.14)
where the real number C M,δt depends on the mesh and the time step δt but not on ε.

Proof of Proposition 5.4. The control of the pressure gradient directly follows from the estimate (5.12) by introducing (as done for the implicit scheme) v n+1 = B M (ρ n+1 -m(ρ n+1 )) ∈ H M,0 (Ω). We multiply the velocity correction step (5.3b) by v σ , multiply by |D σ | and sum over σ ∈ E int to get:

- σ∈Eint |D σ |(∇p ε ) n+1 σ v n+1 σ = 1 δt σ∈Eint |D σ |ρ n Dσ (u n+1 σ -ũn+1 σ )v n+1 σ - σ∈Eint |D σ |(∇p ε ) n σ v n+1 σ
Thanks to the bounds previously derived on ρ, ∇p ε , u and ũ, we bound the right-hand side uniformly with respect to ε. On the left-hand side we have thanks to (3.19)

- σ∈Eint |D σ |(∇p ε ) n σ v n+1 σ = - Ω ∇ E (p n ε ) • Π E v n+1 dx = Ω p n ε div M v n+1 dx.
and therefore

Ω p n+1 ε div M v n+1 dx = Ω p n+1 ε (ρ n+1 -m(ρ n+1 )) dx ≤ C M,δt .
We finally conclude as in the implicit case by splitting the pressure integral into two parts.

We finally pass to the limit exactly as for the full implicit scheme.

Numerical simulations

The numerical test-cases presented here have been run with the pressure-correction scheme programmed in the open source software CALIF 3 S. This software was originally designed for the compressible Navier-Stokes system. We adapted the code to the soft congestion model by setting the barotropic pressure law to (1.3). The inviscid case µ = λ = 0 can be handled by the software CALIF 3 S by setting the viscosity to the order of magnitude of a numerical diffusion, namely µ = h where h is the space step and λ = -2 3 µ. For all the test-cases, we set γ = β = 2.

1D validation test cases

We first present Riemann type test-cases taken from [START_REF] Degond | Numerical simulations of the euler system with congestion constraint[END_REF] for the 1D inviscid soft congestion system. In both test cases we set ε = 10 -4 .

(P1): A double shock wave The initial density and momentum (q = ρu) for the first test-case test-case are given by: (P 1) :

(ρ, q)(x, 0) = (0.7, 0.8)

x ∈ [0, 0.5) (0.7, -0.8) x ∈ (0.5, 1]

In (P 1), the two opposite velocities lead to a rise of the density and the appearance of a congested area in the middle region. The exact solution thus consists in two symmetric left and right going shocks with a congested (up to ≈ ε) intermediate state in between. In order to check that the pressure correction scheme is indeed unconditionally stable with respect to the time step, the scheme has been run for various CFL numbers varying from CF L = 0.27 to CF L = 43.9 and we obtain similar results. However, it is to be noticed that the larger is the CFL number, the worsed conditioned is the Newton method solving the correction step. For CFL with values larger that 44, the Newton solver needs more than 30 steps to converge which causes the program to stop.

A mesh refinement process has also been implemented in order to check numerically the convergence of the method. For this purpose, we compute the discrete L 1 -error between the approximate solution and the exact one at the final time T = N δt = 0.05, normalized by the discrete L 1 -norm of the exact solution:

E(δx) = j |φ N j -φ ex (x j , T )|δx j |φ ex (x j , T )|δx ,
where φ is any of the non-conservative variables (ρ, u, p). The calculations have been implemented on several meshes composed of 200 × 2 n uniform cells with n = 0, 1, .., 7 (knowing that the domain size is L = 1). The time step is adapted in order to fix the CFL number to 1.73 in all computations. In Fig. 3, the error E(δx) at the final time T = 0.05, is plotted against δx in a log -log scale. For the three most refined meshes, we observe the convergence of the method at the expected first order.

(P2): A double rarefaction wave with vacuum The initial density and momentum (q = ρu) for the second test-case are given by:

(P 2) : (ρ, q)(x, 0) = (0.7, -0.8) x ∈ [0, 0.5) (0.7, 0.8) x ∈ (0.5, 1]
The solution of the problem (P 2) consists in a left-going rarefaction wave and a symmetric right-going rarefaction wave giving rise to a vacuum state in between. Fig. 4 displays the exact and numerical solutions.

2D Simulation test-cases

2D collision between blocks This is a 2D inviscid test-case taken from [START_REF] Degond | Numerical simulations of the euler system with congestion constraint[END_REF] which illustrates the collision of two congested regions. The computational domain is Ω = (0, 2 . The initially congested regions are A = [ 1 6 , 5 12 ] × [ 1 3 , 7 12 ] and B = [ 7 12 , 5 6 ] × [ 5 12 , 2 3 ]. The initial density and momentum are given by:

ρ = 0.6 + 0.2 × 1 A∪B , q = 1 0 1 A + -1 0 1 B .
The computation is run on a non structured grid composed of 41648 triangles. We impose a no slip boundary condition on the velocity. The time step is δt = 5.10 -4 . The CPU time of the computation is 13 min for 400 time iterations. Fig. 5 displays the color map for the density at times t = 0.005 and t = 0.2. Fig. 6 displays cut lines of the density and velocity in the x direction on the line y = 0.5 at timet = 0.005. We obtain similar results as in [START_REF] Degond | Numerical simulations of the euler system with congestion constraint[END_REF].

2D convergent corridor This 2D test-case is inspired by a test case in [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF] featuring the global motion of crowd through a convergent corridor. The computational domain described in polar coordinates is the 22 We model the presence of an "exit" (which is actually closed) at the bottom end (r = 0.1) of the cone by adding the following source term:

f = ∇ D -u
where D is the euclidean distance with the origin of the cone (0, 0). This source term tends to force the flow velocity towards the "exit" and the further (x, y) is from the exist, the larger is the velocity. In practice, the part of the source term depending on the velocity unknown is incorporated as an additional implicit reaction term in the prediction step (5.2) of the numerical scheme.

To assess the behaviour of the pressure correction scheme on the viscous model (1.1), we take µ = 0.1 and λ = -2 3 µ in this test case. On the boundary of Ω, we impose a non penetration condition on the velocity u |∂Ω • n = 0. Hence, we expect a congested region to appear at the bottom of the computational domain and to propagate towards the top of the domain.

The computation is run on a non structured grid composed of 12020 triangles. The time step is δt = 10 -2 . The CPU time of the computation is 23 min for 1000 time iterations. Fig. 7 displays the color map for the density at various times. Fig. 8 displays cut lines of the density and velocity with respect to y along the line x = 0 at time t = 10. We observe that the radial velocity is constant (close to zero) in the congested region. In the vicinity of the upper boundary (r = 1) where the system is close to vacuum, the velocity is close to zero which is coherent with (2.3).

Conclusion

We have performed the analysis of a (spatially) staggered and semi-implicit (in time) numerical scheme for the simulation of free-congested viscous and inviscid compressible flows, where congestion phenomena are modelled with a singular pressure p ε (ρ) = ερ γ /(1 -ρ) β enforcing a maximum constraint on the density: ρ < 1. The numerical scheme falls in the framework of pressure correction schemes, it was originally designed for the Navier-Stokes equations and programmed in the open source software CALIF 3 S.

Thanks to an analysis which is similar to that of the continuous model, we prove that both the fully implicit and semi-implicit schemes satisfy a control of the discrete pressure which is independent of the mesh size and of the parameter ε. We also prove that the expected maximum constraint on the density is satisfied. In the limit ε → 0, it is known that the model converges towards a hybrid model combining a pressureless non congested compressible flow (ρ < 1) and an incompressible congested (ρ = 1) flow. At the discrete level, we prove that (up to extraction of a subsequence), the limit pressure only activates in the congested area as we have (1 -ρ)p = 0 at the limit in all the mesh cells. An interesting question is to go further and prove that similarly to the continuous model, the solutions of the numerical scheme for the free-congested model converge towards the solutions of a numerical scheme for the limit hybrid model. This will be studied in a forthcoming work, where we analyse a numerical scheme for the limit hybrid model.

The convergence of the numerical scheme as the time step goes to zero is still an open problem for this time-dependent system. However we refer to [START_REF] Perrin | A convergent FV-FE scheme for the stationary compressible Navier-Stokes equations[END_REF] where such a result was proven for the stationary compressible Navier-Stokes system.

In order to validate the numerical scheme, 1D Riemann type test-cases for the inviscid model have been performed and the discrete solutions have been shown to compare well with the exact solutions. Finally we have performed 2D test-cases which are interesting for the targeted applications, namely the simulation of crowd motion.

  1c) in the domain (0, T ) × Ω for T > 0 and Ω an open bounded subset of R d , d = 2 or 3, with Lipschitz boundary. The quantities ρ ≥ 0 and u = (u 1 , .., u d ) T are respectively the density and velocity of the fluid, while f = (f 1 , .., f d ) T is an external force. System (1.1) is complemented with homogeneous Dirichlet boundary conditions on the velocity: u| ∂Ω = 0. (1.2)

Figure 1 :

 1 Figure 1: Notations for primal cells (delimited with bold lines) and dual cells (in grey).

  we define D σ = D K,σ . We denote by Ẽ(D σ ) the set of faces of D σ , and by ǫ = D σ |D σ ′ the face separating two diamond cells D σ and D σ ′ . We denote by Ẽint the set of dual faces included in the domain and by Ẽext the set of dual faces lying on the boundary ∂Ω.

5 . 11 )

 511 where R n+1 = σ∈Eint |D σ |R n+1 σ + K∈M |K|R n+1 K ≥ 0, with R n+1 σ defined by (5.10) and R n+1 K defined by (4.7).

Fig. 2 Figure 2 :

 22 Figure 2: (P1): space variations of (ρ, u) at time T = 0.05, δx = 5.10 -3 , δt = 5.10 -4 , CF L = 0.27. Solid line : exact solution, red line with points : pressure correction scheme.

21

 21

Figure 3 :Figure 4 :

 34 Figure 3: (P1): L 1 -Error with respect to δx for the pressure correction scheme. CF L = 1.73.

Figure 5 :

 5 Figure 5: 2D collision: color map of the density and (rescaled) velocity vector field at times t = 0.005 (left) and t = 0.2 (right). Number of cells 41648, δt = 5.10 -4 .

Figure 6 :

 6 Figure 6: 2D collision: cut of the density (left) and velocity (right) in the x direction on the line y = 0.5 at timet = 0.005.

Figure 7 :

 7 Figure 7: 2D convergent corridor: color map of the density at times t = 0.0 (top left), t = 0.7 (top right), t = 1.7 (bottom left) and t = 10 (bottom right). Number of cells 12020, δt = 10 -2 .

Figure 8 :

 8 Figure 8: 2D convergent corridor: cut of the density (left) and velocity (right) with respect to y on the line x = 0 at time t = 10.

We begin with the analysis of a fully implicit scheme. As we shall see, the estimates associated with this time discretization is really close to the analysis of the continuous equations.

Acknowledgments

This work was supported by the SingFlows and CRISIS projects, grants ANR-18-CE40-0027 and ANR-20-CE40-0020-01 of the French National Research Agency (ANR). The authors would like to thank J.-C. Latché, F. Babik and G. Boyer of the Institut de Radioprotection et de Sûreté Nucléaire, and B. Fabrège from the Institut Camille Jordan for their help regarding installation and implementation issues with the software CALIF 3 S.