L Batteux

T Gallou Ët

R Herbin

J C Latch É

P Poullet

CONVERGENCE OF THE MAC SCHEME FOR THE INCOMPRESSIBLE NAVIER-STOKES EQUATIONS WITH VARIABLE DENSITY AND VISCOSITY

Keywords: 2000 Mathematics Subject Classification. 35Q30, 65M08, 65N12, 76M12 variable density, incompressible Navier-Stokes equations, staggered finite volumes, MAC scheme, analysis, convergence. 1

The present paper addresses the convergence of the implicit MAC (for Marker-and-Cell) scheme for time-dependent Navier-Stokes equations with variable density and density-dependent viscosity and forcing term. A priori estimates on the unknowns are obtained, and thanks to a topological degree argument, they lead to the existence of a discrete solution at each time step. Then, by compactness arguments relying on these same estimates, we obtain the convergence (up to the extraction of a subsequence), when the space and time steps tend to zero, of the numerical solutions to a limit; this latter is shown to be a weak solution to the continuous problem by passing to the limit in the scheme.

Introduction

We consider in this paper the numerical approximation of the incompressible Navier-Stokes equations with variable density and viscosity, ∂ t ρ + div(ρū) = 0, (1a) ρ∂ t ū + (ū • ∇)ū -div(µ(ρ)D(ū)) + ∇p = f , (1b) div ū = 0, (1c)

in Ω × (0, T ) where T ∈ R + and Ω is an open bounded connected subset of R d , with d ∈ {2, 3}, which may be meshed by a structured grid, and therefore consists of a finite union of rectangles if d = 2 or of rectangular parallelepipeds if d = 3. The variables ρ, ū and p are respectively the density, the velocity and the pressure in the flow, and the three above equations respectively enforce the mass conservation, the momentum conservation and the incompressibility of the flow. The viscosity µ of the fluid is supposed to be a continuous function of the density ρ. The strain rate tensor D is defined as the symmetric part of the velocity gradient, i.e. D(v) = ∇v + t ∇v, for any sufficiently regular vector function v. We assume that the forcing term f either belongs to L 2 (0, T ; L 2 (Ω) d ) or is a continuous function of the density ρ. This system is supplemented with initial conditions and homogeneous Dirichlet boundary conditions: [START_REF] Batteux | Simulation of a particulate flow in 3D using volume penalization methods[END_REF] ū| t=0 = u 0 (x), ρ| t=0 = ρ 0 (x) for a.e. x ∈ Ω, and ū| ∂Ω = 0, with u 0 ∈ L 2 (Ω) d , ρ 0 ∈ L ∞ (Ω) and 0 < ρ min ≤ ρ 0 (x) ≤ ρ max a.e. in Ω.

These equations model a variable density and viscosity flow, which may be encountered in many physical situations. With f (ρ) = ρ g, with g the acceleration of the gravity, we obtain the system of equations governing the natural convection problems. System (1) also applies to the motion of mixtures of immiscible fluids having different densities and viscosities, with possible applications to Rayleigh-Taylor instabilities, tracking of interfaces between fluids in multiphase flows or droplet impact onto a solid or a surface liquid, with however the limitation that surface tension between the different fluids is not taken into account. If the viscosity in one of these fluids is set to a very large value, the system matches what we obtain when describing the fluid-structure interaction by the so-called volume penalization method [START_REF] Liu | Convergence of numerical approximations of the incompressible Navier-Stokes equations with variable density and viscosity[END_REF][START_REF] Batteux | Simulation of a particulate flow in 3D using volume penalization methods[END_REF].

Let us review the formal estimates satisfied by the solution. First, supposing that the velocity field is regular enough (and such a regularity is required in the weak formulation of the problem), a consequence (see for instance [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF]) of equations (1a) and (1c), is the following maximum principle: [START_REF] Boyer | Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models[END_REF] ρ min ≤ ρ(x, t) ≤ ρ max , for a.e. (x, t) ∈ Ω × (0, T ), which shows that the natural regularity for ρ is L ∞ (Ω×(0, T )), ρ min ≤ ρ ≤ ρ max . A classical formal identity, referred to as the kinetic energy balance, allows to derive natural estimates for the velocity ū. Let us take for simplicity f = 0 in Equation (1b), take the inner product of this relation by ū and use twice the mass balance equation (1a) to obtain:

∂ t 1 2 ρ |ū| 2 + div 1 2 ρ |ū| 2 ū -div µ(ρ) D(ū) • ū + ∇p • ū = 0.
Integrating over Ω, one gets by integration by parts, since div ū = 0 and ū|∂Ω = 0, that, for all t ∈ (0, T ): Integrating over the time interval (0, t) yields, once again for all t ∈ (0, T ):

d dt Ω 1
(4)

Ω 1 2 ρ(x, t) |ū(x, t)| 2 dx + t 0 Ω µ(ρ) |D(ū(x, t))| 2 dx dt = Ω 1 2 ρ 0 (x) |u 0 (x)| 2 dx, ∀t ∈ (0, T ).
Thanks to the assumptions on the initial data, the right-hand side of this relation is bounded. If ρ(x, t) ≥ ρ min > 0, the first term thus yields an estimate on u(x, t) L 2 (Ω) d , for t ∈ (0, T ). Let us assume that the function µ is continuous over [ρ min , ρ max ] and satisfies µ(s) ≥ µ min > 0 for s in [ρ min , ρ max ]; this hypothesis is made throughout the paper. Under this assumption, the so-called Korn lemma yields µ(ρ) |D(ū(x, t))| 2 ≥ µ min 2 |∇ū(x, t)| 2 , for (x, t) ∈ Ω × (0, T ). so that the second term of Equation (4) controls the L 2 (0, T ; H 1 0 (Ω) d ) norm. The natural regularity for ū is thus to lie in L ∞ (0, T ; L 2 (Ω) d ) ∩ L 2 (0, T ; H 1 0 (Ω) d ). Together with the estimate (3) of the density, this suggests the following weak formulation of the problem. Definition 1.1. Let ρ 0 ∈ L ∞ (Ω) be such that 0 < ρ min ≤ ρ 0 (x) ≤ ρ max for a.e.

x ∈ Ω, and let u 0 ∈ L 2 (Ω) d . Let µ be a continuous function over [ρ min , ρ max ] such that µ(s) ≥ µ min > 0 for s ∈ [ρ min , ρ max ]. Finally, let f either be a given function of L 2 (Ω × (0, T )) d or be a function f (ρ) of the density, continuous over [ρ min , ρ max ] (in which case f is bounded over the same interval). A pair (ρ, ū) is a weak solution of problem [START_REF] Batteux | Modélisation numérique d'écoulements incompressibles diphasiques : Applications à l'interaction fluide-structure[END_REF] if it satisfies the following properties:

ρ ∈ ρ ∈ L ∞ (Ω × (0, T )), ρ min ≤ ρ ≤ ρ max a.e. in Ω × (0, T ) .

-

ū ∈ v ∈ L ∞ (0, T ; L 2 (Ω) d ) ∩ L 2 (0, T ; E) with E = u ∈ H 1 0 (Ω) d , divu = 0 a.e. in Ω . -For all ϕ in C ∞ c (Ω × [0, T )), (5)
- T 0 Ω ρ (∂ t ϕ + ū • ∇ϕ) dx dt = Ω ρ 0 (x) ϕ(x, 0) dx.
-For all ϕ in w ∈ C ∞ c (Ω × [0, T )) d , divϕ = 0 a.e. in Ω × (0, T ) ,

T 0 Ω -ρ ū • ∂ t ϕ -(ρū ⊗ ū) : ∇ϕ + µ(ρ) D(ū) : D(ϕ) dx dt = Ω ρ 0 u 0 • ϕ(x, 0) dx + T 0 Ω f • ϕ dx dt. (6)
The existence a weak solution to the problem (1) as given in Definition 1.1 was proven in [START_REF] Simon | Nonhomogeneous viscous incompressible fluids: existence of velocity, density, and pressure[END_REF], without dependency on the density of the viscosity and the forcing term. This proof is extended in [START_REF] Lions | Mathematical Topics in Fluid Mechanics -Volume 1: incompressible models[END_REF]Chapter 2], dealing with the possible occurence of void zones, i.e. zones where the density vanishes, and density-dependent viscosities; a recent presentation may be found [START_REF] Boyer | Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models[END_REF]Chapter VI].

The MAC scheme [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface[END_REF] is certainly among the most popular schemes for the solution of Navier-Stokes equations, and, whenever a structured grid may be used, proves to be a very efficient choice and has been used for the discretization of several problems, see e.g. [START_REF] Samelson | Surface pressure Poisson equation formulation of the primitive equations: numerical schemes[END_REF][START_REF] Samelson | A fourth-order numerical method for the planetary geostrophic equations with inviscid geostrophic balance[END_REF][START_REF] Fathy | A fourth order difference scheme for the Maxwell equations on Yee grid[END_REF][START_REF] Chen | Convergence analysis of a fully discrete finite difference scheme for Cahn-Hilliard-Hele-Shaw equation[END_REF]; in particular, for these specific domains, it generally reaches the same accuracy as alternative inf-sup stable discretizations of same order, as low order non-conforming finite elements [START_REF] Latché | A convergent staggered scheme for the variable density incompressible Navier-Stokes equations[END_REF], with a better observed stability and a number of unknowns divided by d for the velocity. The aim of this paper is to prove the convergence of the implicit-in-time scheme based on this discretization, in the following sense: first, the scheme admits solutions; second, given a sequence of grids and time steps, with both the space and time steps tending to zero, we show by compactness arguments that the sequence of discrete solutions converges, up to the extraction of a subsequence, to a limit, and that this limit is a weak solution to the continuous problem in the sense of Definition 1.1. This process requires results from the theory of transport equations [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF], but does not require any existence nor regularity of the solutions; instead, it provides an existence result of weak solutions as a by product of the convergence result. Up to our knowledge,similar studies are scarce : in [START_REF] Liu | Convergence of numerical approximations of the incompressible Navier-Stokes equations with variable density and viscosity[END_REF], the authors deal with a discontinuous Galerkin approximation of System (1); in [START_REF] Latché | A convergent staggered scheme for the variable density incompressible Navier-Stokes equations[END_REF], the authors study a scheme combining finite volumes (for convection terms) and finite elements (for the velocity diffusion term), based on the low-order non-conforming Rannacher-Turek element. The advantage of this latter approach is two-fold: first, under the discrete divergence-free constraint, the discrete mass balance, with a suitable definition of the convection fluxes, preserves the bounds of the density; second, a careful construction of the velocity convection term allows to derive a local (i.e. not integrated with respect to space) kinetic energy balance, without the need of any compensating term. In [START_REF] Latché | A convergent staggered scheme for the variable density incompressible Navier-Stokes equations[END_REF], both the velocity and the forcing term are supposed to be independent of the problem variables. We extend here this latter study in several directions. First, the space discretisation is different. Second, we deal with density-dependent viscosity and forcing term, which needs some technical complements in order to pass to the limit in the scheme, but also implies that we cope with the formulation of the stress tensor which is more realistic from a mechanical viewpoint, namely µD(u) instead of ∆u. This latter point requires to derive at the discrete level the analogue of the so-called Korn lemma; this results is given in Appendix A. A part of the material of this paper, namely the convergence theorem with a constant viscosity and the Laplace operator form for the diffusion term, was announced in the conference paper [START_REF] Gallouët | Convergence of the MAC scheme for variable density flows[END_REF], with a sketch of proof. The proof of convergence of the MAC scheme to a weak solution for the constant density incompressible Navier-Stokes equations is addressed in [START_REF] Gallouët | Convergence of the Marker-And-Cell scheme for the incompressible Navier-Stokes equations on non-uniform grids[END_REF], and several results of this latter paper are used in the sequel.

The analysis in the present paper is restricted to the fully implicit scheme. However, an easy extension allows to deal with semi-implicit schemes, where the advective field in the mass and momentum balance equations is taken at the beginning of the time step; this process decouples the mass balance from the Navier-Stokes equations (i.e. the momentum balance and the divergence constraint), yielding a scheme that is easier to handle. However, for practical applications, pressure correction algorithms are often preferred: see e.g. [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF][START_REF] Wang | Convergence of gauge method for incompressible flow[END_REF][START_REF] Liu | Projection method iii. spatial discretization on the staggered grid[END_REF] in the constant density case and [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF][START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF] in the variable density case; their convergence analysis poses severe difficulties and none of the above-quoted works address the convergence of schemes where momentum balance and divergence free constraint are decoupled. In [START_REF] Wang | Convergence of gauge method for incompressible flow[END_REF][START_REF] Liu | Projection method iii. spatial discretization on the staggered grid[END_REF], error estimates are obtained for the discretization of the constant density incompressible Navier-Stokes equations on square grids provided the exact solution is very regular. In [START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF], first order error estimates for the variable density case are obtained on the velocity, with a Galerkin-type method for the space discretization, provided that the exact solution is regular (in a sense for which no proof of existence is available to this day) and that the density is "well approximated". A second order is also presented therein, but the convergence rate is only conjectured and assessed by numerical results. Our process here is quite different. Indeed, as already mentioned, we do not assume anye regularity assumption on the exact solution. In fact, we do not even assume that there exists a solution; this is obtained as a by-product of the proof of convergence of the scheme.

Let us present here the main steps of this proof. Consider a sequence of approximate densities and velocities (ρ (m) , u (m) ) m∈N computed on a series of grids and time steps indexed by m, such that both the space and time steps tend to 0 as m tends to infinity.

• From the maximum principle which is verified by the approximate density, one obtains the existence of a discrete solution and the weak convergence of the approximate densities to a limit ρ, up to a subsequence. • Owing to the estimates of the velocity obtained by a Korn inequality (needed because of the non constant viscosity), the weak convergence of the velocity is obtained. Estimates on the time translates of the velocity are then established, so that an Aubin-Simon type time-compactness result leads to the convergence of the velocity in L q (Ω × (0, T )), 1 ≤ q < +∞. • Passing to the limit on the discrete mass equation, we get that the limit (ρ, ū) satisfies the weak mass balance equation. • The convergence of the sequence (ρ (m) ) m∈N and consequently the convergence of the sequence (µ(ρ (m) )) m∈N in L q (Ω × (0, T )), 1 ≤ q < +∞ is obtained by remarking that, thanks to the fact that divu = 0, Ω (ρ (m) ) 2 dx is conserved. • The last step consists in passing to the limit in the momentum balance equation to show that the limit (ρ, ū) of the approximate solutions is a weak soltuion to the problem. As a by product, we get the existence of a weak solution of the problem.

The paper is organized as follows. The MAC discretization is described in Section 2, then the considered scheme is given in Section 3, together with some (discrete) continuity and stability properties of the discrete operators. Section 4 is devoted to establish a priori estimates for the discrete solutions, i.e. the classical L ∞ bound on the density an L 2 (H 1 0 ) d ∩L ∞ (L 2 ) d estimate for the velocity, and then, on this basis, to prove their existence by a topological degree argument. Finally, in Section 5, we establish the essential result of this paper, namely the convergence property stated above. In the appendix, we state and prove two results which seem interesting for their own sake: the discrete Korn lemma for the MAC scheme, and the stability and convergence of general transfer operators for piecewise constant functions from one mesh to another.

MAC discretization and discrete unknowns

Consider a domain Ω which is a connected union of disjoint rectangular domains, and a coordinate system such that the edges (respectively the faces) of these rectangles (respectively parallelepipeds) are orthogonal to one vector of the canonical basis of R d , (e (1) , . . . , e (d) ). The mesh of Ω is defined as follows (see Figure 1). Definition 2.1 (MAC grid). A mesh associated to the MAC discretization of Ω, referred to by M, is defined by: -a primal (or pressure) grid M, which consists in a conforming structured partition of Ω made of rectangles if d = 2 or rectangular parallelepipeds if d = 3. A generic element K of M is called a primal cell, and we note x K its center of mass. -E is the set of edges (d = 2) or faces (d = 3) of the primal grid; throughout the following, for short, we will use "face" to denote an element of E, whatever the space dimension may be. We split

E into E = E int ∪ E ext
, where E int (resp. E ext ) are the faces of E that lie in the interior (resp. on the boundary) of the domain. For i ∈ [ 1, d ], we denote by E (i) ⊂ E the set of the faces which are orthogonal to e (i) , and we also split this set into internal and boundary faces

E (i) = E (i) int ∪ E (i)
ext . The set of faces of a primal cell K is denoted by E(K) and, for σ ∈ E(K), we define D K,σ as the half-cell of K adjacent to σ. An internal face σ separating the primal cells K and L is denoted by σ = K|L, and we define the dual cell D σ associated to σ by D σ = D K,σ ∪ D L,σ . For a face σ ∈ E ext , we set D σ = D K,σ , with K the primal cell adjacent to K. For any face σ, x σ = (x σ,1 , . . . , x σ,d ) stands for the mass center of σ.

For i ∈ [ 1, d ], the set {D σ , σ ∈ E (i) } defines a partition of Ω, which is referred to as the i-th dual mesh.

The discretization of the momentum balance equation hinges upon the sets of faces of E (i) of the i-th dual mesh, i ∈ [ 1, d ].. We distinguish once more the internal elements of E (i) from the external elements by writing

E (i) = E (i) int ∪ E (i) ext . For ∈ E (i) int we note = σ|σ if (σ, σ ) ∈ (E (i) ) 2 is such that ∂D σ ∩ ∂D σ = . D σ σ = K|L L K + x K σ DK,σ σ σ = σ|σ Figure 1. Representation of (M, E) for d = 2.
The discretization of the momentum equation requires the introduction of the faces of the i-th dual mesh, denoted by E (i) . Given a generic dual face ∈ E (i) , we distinguish three situations (see Fig. 2): The set E (i) is decomposed into three subsets,

E (i) = E (i) int ∪ E (i) ext ∪ E (i)
rec , according to the location of the edge: -

E (i) int = { ∈ E (i) ; ⊂ Ω} , -E (i) ext = { ∈ E (i) ; ⊂ ∂Ω}, -E (i) rec = { ∈ E (i) ; = int ∪ ext with int = ∩ Ω and ext = ∩ ∂Ω}. For ∈ E (i) int (resp. ∈ E (i) rec ) we note = σ|σ if (σ, σ ) ∈ (E (i) ) 2 is such that ∂D σ ∩ ∂D σ = (resp. ∂D σ ∩ ∂D σ = int )
. Finally, we define the set of faces of the i-th dual mesh normal to e j : E (i,j) = { ∈ E (i) , ⊥ e (j) }.

Remark 2.1 (Elements of E (i) rec ).
If Ω exhibits a re-entrant corner, there will be instances where

(σ, σ ) ∈ E (i) int × E (i) ext
are such that ∂D σ ∩ ∂D σ = ∅ is normal to e j , with i = j. Such cases make up the elements of

E (i) rec and we have int = ∩ ∂D σ for ∈ E(D σ ) ∩ E (i) rec .
The space step is defined by: [START_REF] Droniou | The gradient discretisation method[END_REF] h M = max{diam(K), K ∈ M}.

Finally, the regularity of the mesh is measured by the following parameter:

(8) θ M = max |σ| |σ | , (σ, σ ) ∈ E , ∂Ω σ σ = σ|σ σ σ = σ|σ σ σ = σ|σ
Figure 2. Definition of = σ|σ in the two-dimensional case for (σ, σ ) ∈ E (1) . Left: ∈ E (1,1) . Middle:

∈ E (1,2) ∩ E (1) int . Right: ∈ E (1,2) ∩ E (1) rec .
with |•| designating the Lebesgue measure, this notation being used in the following for either the R d or R d-1 measure.

Remark 2.2 (Quasi-uniformity of the mesh). The quantity θ M defined by Relation (8) measures in fact the quasi-uniformity of the mesh. The stability and convergence proofs in this paper only require a weaker assumption, namely the fact that the ratio between the size of two adjacent cells remains bounded (which, for a sequence of more and more refined general meshes, does not prevent θ M to tend to +∞). Unfortunately, this latter requirement combined with the structured character of the mesh implies quasi-uniformity.

Remark 2.3 (Assumption (8) and reconstruction operators). Let a sequence of meshes (M (m) ) m∈N be given, and let us suppose that θ M (m) ≤ θ, ∀m ∈ N, with θ a given real number. Then the dual meshes are quasi-uniform with respect to the primal one (and, conversely, the primal mesh is quasi-uniform with respect to dual ones) in the sense specified in Remark B.1. Therefore, the stability and convergence of reconstruction operators defined in Appendix B easily hold, provided that they use a constant stencil (see once again Remark B.1). The same occurs with the "gradient meshes" which will be defined in the following, together with the discrete diffusion operator.

For the time discretization, we consider a partition of the time interval [0, T ], which we suppose uniform to alleviate the notations. We denote by δt the constant time step, the integer number N = T /δt is the number of of time steps and the time t n is defined by t n = n δt, for 0 ≤ n ≤ N .

The scalar discrete fields, namely the density and the pressure, are associated to the primal mesh, so the associated unknowns read (ρ n K ) K∈M, n∈[ 0,N ] and (p n K ) K∈M, n∈[ 0,N ] , respectively. A discrete velocity field reads u = (u 1 , . . . , u d ) t and its i-th component is associated to the i-th dual mesh. In addition, to take into account (a part of) the homogeneous Dirichlet boundary conditions, the normal velocities unknowns associated to external faces are set to zero. Hence, for i ∈ [ 1, d ], the unknowns for u i read (u n σ ) σ∈E (i) int , n∈[ 0,N ] . Note that, consequently, the component of the velocity associated to an unknown is not specified by its notation but viewed only by the orientation of the associated face.

The discrete scheme

The discrete scheme considered here reads: Let u 0 and ρ 0 be given and solve, for 0 ≤ n ≤ N -1,

ð t ρ n+1 K + div(ρu) n+1 K = 0, ∀K ∈ M, (9a) ð t (ρu i ) n+1 σ + div(ρ u i u) n+1 σ -div(µ D(u)) n+1 σ + (∇p) n+1 σ = f n+1 σ , ∀σ ∈ E (i) int , for i ∈ [ 1, d ], (9b) (divu) n+1 K = 0, ∀K ∈ M, (9c)
where (9a), (9b), (9c) are the backward-in-time finite volume discretizations of (1a), (1b) and (1c) respectively, on the primal mesh for the mass balance equation and the divergence-free constraint, and on the dual mesh associated to the i-th component of the velocity for the i-th component of the momentum balance equation; the terms of these equations are detailed below. Note that no equation is written for the momentum balance on external meshes, since the velocity is prescribed to zero on the boundary. For the latter reason, the system is singular (the sum over the cells of Equation ( 9c) is zero by conservativity, see below for the definition of (divu) n+1 K ), and must be complemented by the condition (10)

K∈M |K| p K = 0,
which states that the mean value of the pressure is zero.

Initialization of the scheme and forcing term -The discrete initial data (ρ 0 , u 0 ) is obtained by averaging (ρ 0 , u 0 ) on the primal and dual cells, respectively:

(11)

ρ 0 K = 1 |K| K ρ 0 (x) dx, ∀K ∈ M, u 0 σ = 1 |D σ | Dσ u i,0 (x) dx, ∀σ ∈ E (i) int , for i ∈ [ 1, d ].
When f is a given function, the forcing term in the momentum balance equation is also obtained by averaging f (which is assumed to lie in L 2 (Ω × (0, T )) d ):

(12)

f n+1 σ = 1 δt |D σ | tn+1 tn Dσ f i (x, t) dx dt, ∀σ ∈ E (i) int , for i ∈ [ 1, d ].
When f is a function of the density, we set

(13) f n+1 σ = f i (ρ n+1 σ ), ∀σ ∈ E (i) int , for i ∈ [ 1, d ].
where ρn+1 σ may be any reasonable convex combination of the density at the n + 1 time level in the neighbouring cells.

Mass balance equation -The definition of the discrete time partial derivative of the density is standard:

ð t ρ n+1 K = 1 δt (ρ n+1 K -ρ n K ), ∀K ∈ M, 0 ≤ n ≤ N -1.
The mass convection operator is obtained by a first-order upwind scheme. For

i ∈ [ 1, d ] and σ ∈ E (i) int , let us denote by u n+1 Kσ the quantity u Kσ = u n+1 σ e (i) • n K,σ .
Then, for K ∈ M and 0 ≤ n ≤ N -1, ( 14) div(ρu

) n+1 K = 1 |K| σ∈E(K), σ=K|L F n+1 K,σ , with F n+1 K,σ = |σ| ρ n+1 σ u n+1 Kσ and ρ n+1 σ = ρ n+1 K if u Kσ ≥ 0, ρ n+1 L otherwise.
.

In this relation, F n+1 K,σ is the (conservative) mass flux through the face σ outward K. The summation index "σ ∈ E(K), σ = K|L" implicitly assumes that we restrict the sum to the internal faces of K, which in turns implicitly assumes that the mass flux vanishes across the external faces, which is consistent with impermeability boundary conditions.

Pressure gradient and velocity divergence -The discretization of these terms is standard for the MAC scheme:

(15) (∇p) n+1 σ = |σ| |D σ | (p n+1 L -p n+1 K )(n K,σ • e i ), ∀σ = K|L ∈ E (i) int , i ∈ [ 1, d ], (divu) n+1 K
= 1 |K| σ∈E(K) |σ| u n+1 K,σ ,
where we recall that, for i ∈

[ 1, d ] and σ ∈ E(K) ∩ E (i) , u K,σ = u σ e (i) • n K,σ .
These two operators satisfy the following duality property.

Lemma 3.1 (grad-div duality). Let q and v be a discrete pressure and velocity, respectively. Then,

σ∈E |D σ | (∇q) σ v σ + K∈M |K| q K (divv) K = 0. (16)
Diffusion term -For the velocity diffusion term, we use a weak formulation of the MAC diffusion scheme which was already introduced [START_REF] Grapsas | An unconditionally stable staggered pressure correction scheme for the compressible Navier-Stokes equations[END_REF]. The first step is to define piecewise constant partial derivatives of the velocity components, based on specific partitions of the computational domain. Specifically, for i, j ∈ [ 1, d ], the discrete partial derivative of the i-th component of the velocity u i with respect to the j-th coordinate, which we denote by ð j u i , is piecewise-constant over each volume

D for ∈ E (i,j) ∩ E (i) rec , recalling that E (i,j) = { ∈ E (i)
, ⊥ e (j) }. Henceforth we associate one cell D to the elements of E rec mostly for the sake of convenience. The dual cell is defined as,

(17) D = × [x σ , x σ ], for = σ|σ ∈ E (i) int , × [x σ , x σ, ], for ∈ E(D σ ) ∩ E (i) ext , D σ, ∪ D σ , for = σ|σ ∈ E (i) rec , with D σ, = × [x σ , x σ, ], D σ , = int × [x σ , x σ, ] if (σ, σ ) ∈ E (i) int × E (i)
ext and where x σ, = (x σ, ,1 , . . . , x σ, ,d ) refers to the orthogonal projection of x σ on . Note that the set {D , ∈ E (i,j) } is a partition of Ω; a volume D of this set is called in the following a (i, j)-gradient cell. For the two-dimensional case, these volumes are sketched on Figure 3. When i = j, a (i, j)-gradient cell coincides with a primal cell.

In addition, we observe that, for i = j, the set of the (i, j)-gradient cells and the set of the (j, i)-gradient cells are the same; in two space dimensions, such a (i, j)gradient cell may be associated to a grid vertex, while, in three space dimensions, for = σ|σ , it is associated to the edge equal to σ ∩ = σ ∩ (see [START_REF] Grapsas | An unconditionally stable staggered pressure correction scheme for the compressible Navier-Stokes equations[END_REF]).

D u σ u σ (a) D u σ u σ (b) D D xσ, x σ , × × u σ u σ (c) D u σ u σ (d) Figure 3. (i, j)-gradient cells in the two-dimensional case. (a): (i, j) = (1, 1), = σ|σ ∈ E (1)
int and ⊥ e (1) , (b):

(i, j) = (2, 1), = σ|σ ∈ E (1)
int and ⊥ e (2) , (c):

(i, j) = (2, 1), = E (1) ext ∩ E(Dσ) and (i, j) = (1, 2), = E (2) ext ∩ E(D σ ), (d): (i, j) = (2, 1), = σ|σ ∈ E ( rec1).
For i, j ∈ [ 1, d ], we define ð j u i a.e. in Ω by [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface[END_REF]

ð j u i (x) = u σ -u σ d(x σ , x σ ) n • e j for x ∈ D , ∈ E (i,j) ∩ E (i) int , = σ|σ , -u σ d(x σ , x σ, ) n • e j for x ∈ D , ∈ E (i,j) ∩ E (i) ext , ∈ E(D σ ). -u σ d(x σ , x σ, ) n • e j X Dσ, (x) for x ∈ D , ∈ E (i) rec ∩ E(D σ ), σ ∈ E (i) int .
Note that the last formula implicitly take into account the homogeneous Dirichlet condition satisfied by the velocity (and it is the only consequence in the scheme formulation that tangential velocities are prescribed to zero at the boundary). For

σ ∈ E (i)
int , let the discrete velocity field ϕ σ be the discrete velocity function defined by (ϕ σ ) σ = 1 and (ϕ σ ) σ = 0 for σ ∈ E, σ = σ (so the j-th component(s) of ϕ σ are zero for j = i and the i-th component has only one degree of freedom set to 1, namely the degree of freedom corresponding to σ). To each gradient cell (which, for i = j, is both a (i, j)-gradient cell and a (j, i)-gradient cell), we associate a viscosity µ D , and we introduce d × d viscosity fields defined a.e. in Ω by

µ (i,j) (x) = µ D , for x ∈ D , with ∈ E (i,j) .
For a discrete velocity field u, we are now in position to define a discrete gradient ∇ E u and a tensor associated to the multiplication of the strain rate by the viscosity, denoted by (µ D) E (u), for a.e. x ∈ Ω:

(19) (∇ E u) i,j (x) = ð j u i (x), ((µD) E (u)) i,j = µ (i,j) (x) ð j u i (x) + ð i u j (x) 2 .
The discrete gradient tensor is associated to the discrete H 1 0 norm:

(20) |u i | 2 1,E = d j=1 Ω ð j u i (x) 2 dx, |u| 2 1,E = d i=1 |u i | 2 1,E ,
which turns out to be the usual finite volumes H 1 0 discrete norm, known to dominate the L 2 -norm by a discrete Poincaré inequality [9, Lemma 9.1]. Finally, we define the diffusion term by: ( 21)

-|D σ | div(µD(u)) n+1 σ = Ω (µD) E (u)(x) n+1 : ∇ E (ϕ σ )(x) dx,
where the time index in (µD) E (u)(x) n+1 means that the viscosity fields and the velocity in expression [START_REF] Herbin | On some implicit and semi-implicit staggered schemes for the shallow water and Euler equations[END_REF] must be taken at t n+1 . It is shown in [START_REF] Grapsas | An unconditionally stable staggered pressure correction scheme for the compressible Navier-Stokes equations[END_REF] that the definition (21) of the diffusion term is the same as the standard MAC formulation, up to the specific definition of the viscosity, i.e. that we have a formulation of the form, for σ ∈ E

(i) int : -|D σ | div(µD(u)) n+1 σ = ∈ E(Dσ) | | µ n+1 D D(u n+1 )| σ • n σ, • e (i) ,
where D(u n+1 )| σ stands for an approximation of D(u n+1 ) obtained by replacing as usual the partial derivatives by differential quotients. The quantity µ n+1 D may be approximated by applying the function µ to any reasonable approximation of ρ n+1 in D σ , for instance: The (µ D) E tensor and the diffusion term enjoys the following properties.

µ n+1 D = µ(ρ n+1 K ) if D = K,
Lemma 3.2. Let u and v be two discrete velocities. The diffusion term defined by (21) satisfies the following discrete integration by part formula:

- d i=1 σ∈E∈E (i) int |D σ | div(µD(u)) σ v σ = Ω (µ D) E (u)(x) : ∇ E (v)(x) dx.
In addition, the (µ D) E tensor given by ( 19) is symmetrical and, if, for i, j

∈ [ 1, d ], the functions µ (i,j) satisfy 0 ≤ µ min ≤ µ (i,j) (x) ≤ µ max a.e. in Ω, µ min 2 Ω ∇ E (u)(x) : ∇ E (u)(x) dx ≤ Ω (µ D) E (u)(x) : ∇ E (u)(x) dx ≤ µ max Ω ∇ E (u)(x) : ∇ E (u)(x) dx.
Proof. As usual in Galerkin methods, the discrete integration by part formula is a straightforward consequence of the weak formulation of the diffusion term. The symmetry of the (µ D) E tensor is a consequence of the fact that, for i, j ∈ [ 1, d ], the gradient cells associated to ð j u i are the same as the gradient cells associated to ð i u j , so the fields µ (i,j) and µ (j,i) are the same. Since (µ D) E is symmetrical, we have, by a standard identity for the tensors contraction,

Ω (µ D) E (u) : ∇ E (u) dx = Ω 1 2 (µ D) E (u) : (∇ E u + ∇ t E u) dx,
and, since

(µ D) E (u)(x) : (∇ E (u)(x) + ∇ t E (u)(x)) ≥ µ min 2 ∇ E (u)(x) + ∇ t E (u)(x) : ∇ E (u)(x) + ∇ t E (u)(x)
for a.e. x ∈ Ω,

Ω (µ D) E (u) : ∇ E u dx ≥ µ min 4 Ω (∇ E u + ∇ t E u) : (∇ E u + ∇ t E u) dx.
The lower bound in the conclusion of the lemma follows by the discrete Korn Lemma proven in Appendix A. For the upper bound of the left-hand side, we write:

(µ D) E (u)(x) : (∇ E (u)(x) + ∇ t E (u)(x)) ≤ µ max ∇ E (u)(x) + ∇ t E (u)(x) : ∇ E (u)(x) + ∇ t E (u)(x) for a.e. x ∈ Ω, and use Ω ∇ E (u)(x) : ∇ t E (u)(x) dx ≤ Ω ∇ E (u)(x) : ∇ E (u)(x) dx,
which is an easy consequence of the Cauchy-Schwarz inequality, using the fact that

∇ t E (u)(x) : ∇ t E (u)(x) = ∇ E (u)(x) : ∇ E (u)(x) for a.e. x ∈ Ω.
Momentum convection term -We now carry out the discretization of the term

∂ t (ρu i ) + div(ρu i u) in (1b), i ∈ [ 1, d ],
to obtain the terms denoted by ð t (ρu i ) n+1 σ + div(ρu i u) n+1 σ in (9b). As shown in [START_REF] Herbin | Short notes: Kinetic energy control in the MAC discretization of compressible Navier-Stokes equations[END_REF], the derivation of a discrete kinetic energy balance requires for these terms to take the following structure, for σ ∈ E

(i) int : ð t (ρu i ) n+1 σ + div(ρu i u) n+1 σ = 1 δt (ρ n+1 Dσ u n+1 σ -ρ n Dσ u n σ ) + 1 |D σ | ∈ E(Dσ) F n+1 σ, u n+1 ,
where he face densities ρ n+1 Dσ and ρ n Dσ and the mass fluxes through the dual faces F n+1 σ, satisfy the following mass balance over D σ :

1 δt (ρ n+1 Dσ -ρ n Dσ ) + 1 |D σ | ∈ E(Dσ) F n+1 σ, = 0.
This is ensured by the following construction. The face densities are obtained by a weighted average of the density in the neighbouring cells:

|D σ | ρ k Dσ = |D K,σ | ρ k K + |D L,σ | ρ k L , for k = n and k = n + 1.
Through included in the boundary, the dual mass flux is set to zero. For an internal dual face, its expression depends on the normal to (see Figure 4). If is parallel to σ, there exists K ∈ M and σ ∈ E (i) such that = σ|σ ⊂ K; in this case, we set

(22) F n+1 σ, = 1 2 (-F n+1 K,σ + F n+1 K,σ ).
Otherwise, denoting by K and L the primal cells adjacent to σ, we observe that there exist two faces of the primal mesh τ ∈ E(K) and τ ∈ E(L) such that is the union of the half of each of these two faces. We then define (23)

F n+1 σ, = 1 2 (F n+1 K,τ + F n+1 L,τ ).
Finally, the velocity at the dual face is approximated by a first-order upwind scheme: for an internal dual face = σ|σ , ( 24)

u n+1 = u n+1 σ if F n+1 σ, ≥ 0 u n+1 σ otherwise.
For external dual faces, since the dual mass fluxes vanish, no definition of u n+1 is required in this case.

σ σ K Fσ, -FK,σ F K,σ σ σ τ FK,τ τ F K,τ Fσ, K L
int . Left: for ⊥ e (1) . Right: for ⊥ e (2) .

Remark 3.1. In any case, using the last two relations and the associated notations, we may recast F σ, as

F σ, = | | ρ û n σ, • e (j)
, with (ρ , û ) given by:

(25) (ρ , û ) =        ρ σ + ρ σ 2 , ρ σ u σ + ρ σ u σ ρ σ + ρ σ if j = i, |τ | ρ τ + |τ | ρ τ |τ | + |τ | , |τ | ρ τ u τ + |τ | ρ τ u τ |τ | ρ τ + |τ | ρ τ otherwise.
Since, for σ = K|L, the face approximation of the density ρ σ is a convex combination of the density in the adjacent cells ρ K and ρ L (in fact, for the upwind scheme, ρ σ is equal to either ρ K or ρ L ), the density ρ is itself a convex combination of the density in the neighbouring cells (precisely speaking, the cells M ∈ M such that M ∩ D is either a (d-1)-surface or a volume, but not empty). From Equation (25), supposing that the densities are positive (which is indeed the case by construction of the scheme, see Section 4), it is clear that û is a convex combination of u σ and u σ if j = i and of u τ and u τ otherwise.

We now prove two stability results for the convection operator which will be used further for the estimation of the time-translates of the velocity. To this purpose, for i ∈ [ 1, d ], to the unknowns (u σ ) σ∈E (i) int of a i-th component of a discrete velocity field, we associate the piecewise constant function u defined for a.e. x ∈ Ω by u(x) = σ∈E (i) int u σ X Dσ (x), with X Dσ the characteristic function of D σ , and the following discrete L q (Ω) norm, q ≥ 1,

u L q (Ω) = σ∈E (i) int |D σ | |v σ | q 1/q ,
which is the standard L q -norm of the function u. To the discrete unkowns (u σ ) σ∈Eint related to a vector field u = (u 1 , . . . , u d ) t ), we associate a discrete L q (Ω) d norm by

u q L q (Ω) = d i=1 u i q L q (Ω) .
We recall that, for any i-th discrete component of a velocity field u, the usual finite-volume discrete H 1 -norm is denoted by |u i | 1,E and defined by [START_REF] Herbin | Short notes: Kinetic energy control in the MAC discretization of compressible Navier-Stokes equations[END_REF]. Finally, for a discrete density field ρ associated to (ρ K ) K∈M , we set

ρ L ∞ (Ω) = max K∈M ρ K .
We are now in position to state the following results. They consist in discrete analogues of the following continuous estimates:

Ω div(ρ u i v) w i dx ≤ ρ L ∞ (Ω) w i H 1 (Ω) v L 4 (Ω) d u i L 4 (Ω) , Ω u i w i div(ρ v) dx ≤ ρ L ∞ (Ω) u i H 1 (Ω) v H 1 (Ω) d w i H 1 (Ω) .
Lemma 3.3 (Continuity results related to the convection operator).

Let i ∈ [ 1, d ] and let (u σ ) E∈E (i) int , (v σ ) E∈Eint and (w σ ) E∈E (i) int
be three families of real numbers; the first and last ones correspond to the i-th component of discrete velocity fields u and w, and the the second one corresponds all the components of the velocity field v. Let (ρ K ) K∈M ⊂ R be given, and satisfy ρ K > 0, ∀K ∈ M. Let F σ, (ρ, v) be defined by ( 22)-( 23) and let u be defined by [START_REF] Nicolaides | Analysis and convergence of the MAC scheme. II. Navier-Stokes equations[END_REF]. Finally, let us assume that there exists θ > 0 such that the parameter θ M measuring the regularity of the mesh satisfies θ M ≤ θ. Then we have:

σ∈E (i) int w σ ∈ E(Dσ) F σ, (ρ, v) u ≤ C ρ L ∞ (Ω) u L 4 (Ω) v L 4 (Ω) d w 1,E , (26)
σ∈E (i) int u σ w σ ∈ E(Dσ) F σ, (ρ, v) ≤ C ρ L ∞ (Ω) u 1,E v L 4 (Ω) d w 1,E . ( 27
)
where the positive real number C only depends on θ, d and Ω.

Proof. Reordering the sums, we get: Thanks to Remark 3.1 and using the same notations, we get:

I = σ∈E (i) int w σ ∈ E(Dσ) F σ, (ρ, v) u = ∈ E (i) , =σ|σ F σ, (ρ, v) u (w σ -w σ ). D u σ , v σ , w σ u σ , v σ , w σ u = uσ or u σ (v1) = ρσvσ + ρ σ v σ ρσ + ρ σ (ð1w) = w σ -wσ x σ ,1 -xσ,1 D u σ , w σ u σ , w σ v τ v τ u = uσ or u σ (v2) = |τ | ρτ vτ + |τ | ρ τ v τ |τ | ρτ + |τ | ρ τ (ð2w) = w σ -wσ x σ ,2 -xσ,2
(28) |I| ≤ d j=1 ∈ E (i,j) , =σ|σ |D | ρ |(v j ) | |u | w σ -w σ x σ,j -x σ ,j with ρ ≤ ρ L ∞ (Ω)
, since ρ is a convex combination of the density in the neighbouring cells. Still by the convex combination properties stated in Remark 3.1, the functions

∈ E (i,j) u X D (x)
and

∈ E (i,j) (v j ) X D (x), for j ∈ [ 1, d ],
are reconstructions of u and v j over the mesh {D , ∈ E (i,j) } of Ω in the sense of Lemma B.1, respectively; the expression of these velocities is given on Figure 5, in the two-dimensional case and for i = 1. Invoking the stability property stated in Lemma B.1 and standard L q -estimates, we thus get

|I| ≤ C d j=1 ρ L ∞ (Ω) u L 4 (Ω) v j L 4 (Ω) d |ð j w| L 2 (Ω) ,
which yields the inequality [START_REF] Samelson | A fourth-order numerical method for the planetary geostrophic equations with inviscid geostrophic balance[END_REF].

By the same reordering of the sums as for I, we get:

σ∈E (i) int u σ w σ ∈ E(Dσ) F σ, (ρ, v) = ∈ E (i) , =σ|σ F σ, (ρ, v) (u σ w σ -u σ w σ ) = d j=1 ∈ E (i,j) , =σ|σ |D | ρ (v j ) u σ w σ -u σ w σ x σ,j -x σ ,j .
The proof then ends by developing the numerator of the fraction thanks to the identity 2(ab

-cd) = (a -c)(b + d) + (a + c)(b -d), for (a, b, c, d) ∈ R 4
, then invoking the same arguments as for the estimate on I, and finally using the fact that, since Ω is bounded, the L 4 -norm is controlled by the L 6 -norm, itself controlled by the discrete H 1 0 -norm thanks to discrete Sobolev embedding results (see [7, Lemma B14]).

Estimates and Existence results

The aim of this section is two-fold: first, we prove a priori estimates for the discrete solutions, and then we establish that the scheme is well-posed, in the sense that it admits at least a solution.

ρ(x, t) = N -1 n=0 K∈M ρ n+1 K X K (x) X (tn,t n+1] (t), p(x, t) = N -1 n=0 K∈M p n+1 K X K (x) X (tn,tn+1] (t),
where, for A ⊂ Ω and x ∈ Ω, X A (x) = 1 if x ∈ A, X K (x) = 0 otherwise and, for A ⊂ (0, T ) and t ∈ (0, T ), X A (t) = 1 if t ∈ A, X A (t) = 0 otherwise. Similarly, to (u n σ ) σ∈E (i) , n∈[ 0,N ] , we associate the discrete function representing a i-th component of a discrete velocity:

(30) u i (x, t) = N -1 n=0 σ∈E (i) u n+1 σ X Dσ (x) X (tn,tn+1] (t). Recall that u n+1 σ = 0 of σ ∈ E (i) ext .
A discrete velocity is a vector valued function of the form u = (u 1 , . . . u d ) t . For such a function ρ, p, u and for n ∈ [ 1, N ], we denote by ρ n , p n and u n = (u n 1 , . . . u n d ) t the function depending on the space variable only and corresponding to the value taken by its time-dependent counterpart over the (t n , t n-1 ) interval. We define the discrete L 2 (0, T ; H 1 0 (Ω) d ) norm of the discrete velocity and its components by:

u 2 L 2 (0,T ;H 1 E ) = d i=1 u i 2 L 2 (0,T ;H 1 E ) with u i 2 L 2 (0,T ;H 1 E ) = N n=1 δt |u n i | 2 1,E .
The set of functions (ρ, u, p) is said to be a solution to the scheme (9) (or to an equation of ( 9)) if the associated families of real numbers do so.

The following estimates for the density are classical consequences of the upwind choice for the discretization of the convection term in the mass balance equation and of the fact that the velocity is divergence-free (see [START_REF] Eymard | The finite volume method. Handbook for Numerical Analysis[END_REF]). Lemma 4.1 (Estimates of the density). Let the initial data ρ 0 satisfy 0 < ρ min ≤ ρ 0 (x) ≤ ρ max for a.e. x ∈ Ω. Let ρ and u satisfy the discrete mass balance equation (9a) and divergence constraint (9c). Then, for a.e x ∈ Ω and t ∈ (0, T ), (31)

ρ min ≤ ρ(x, t) ≤ ρ max .
In addition, ρ satisfies the following discrete entropy balance, for K ∈ M and n ∈ [ 0, N -1 ]:

(32) |K| 2 δt (ρ n+1 K ) 2 -(ρ n K ) 2 + 1 2 σ∈E(K) |σ| (ρ n+1 σ ) 2 u n+1 K,σ + R n+1 K = 0,
where the non-negative remainder term R n+1 K is given by

(33) R n+1 K = |K| 2 δt ρ n+1 K -ρ n K 2 + 1 2 σ∈E(K), σ=K|L |σ| ρ n+1 K -ρ n+1 L 2 (u n+1 K,σ ) -,
with, for a ∈ R, a -= -min(a, 0). Summing over the time steps, we obtain the so-called weak BV-estimate:

(34) N n=1 δt σ∈Eint, σ=K|L |σ| |u n K,σ | ρ n L -ρ n K 2 ≤ C,
where C ≥ 0 only depends on the L 2 -norm of the initial data ρ 0 L 2 (Ω) .

The next lemma states a discrete equivalent of the continuous L 2 (0, T ; H 1 0 (Ω) d ) and L ∞ (0, T ; L 2 (Ω) d ) estimates for the velocity. Lemma 4.2 (Discrete L 2 (0, T ; H 1 0 (Ω) d ) and L ∞ (0, T ; L 2 (Ω) d ) velocity estimates). There exists C > 0 depending only on Ω and µ min such that any solution u of the scheme (9) satisfies, for n ∈ [ 0, N -1 ] :

(35) 1 2 d i=1 σ∈E (i) int |D σ | ρ n+1 Dσ (u n+1 σ ) 2 + µ min 4 δt |u n+1 | 2 1,E ≤ 1 2 d i=1 σ∈E (i) int |D σ | ρ n Dσ (u n σ ) 2 + C f 2 L 2 (Ω×(tn,tn+1)) .
Consequently, there exists C > 0 depending only on Ω, u 0 , ρ 0 , µ min and the L 2norm of f such that any solution u of the scheme (9) satisfies:

(36) u L 2 (0,T ;H 1 E ) ≤ C and u L ∞ (0,T ;L 2 (Ω) d ) = max 0≤n≤N -1 u n+1 L 2 (Ω) d ≤ C. The L 2 -
(i) int |D σ | ρ n+1 Dσ (u n+1 σ ) 2 + µ min 2 δt |u n+1 | 2 1,E ≤ 1 2 d i=1 σ∈E (i) int |D σ | ρ n Dσ (u n σ ) 2 + δt Ω u n+1 • f n+1 dx.
Inequality (35) follows by the Cauchy-Schwarz inequality and the discrete Poincaré estimate [9, Lemma 9.1]. We then get Relations (36) by summing over the time steps and using the fact that the face densities are convex combinations of the cell densities, and thus bounded by below by ρ min .

4.2.

Existence of a solution to the scheme. The proof of the existence of a solution is obtained by a topological degree arguments (see e.g. [START_REF] Deimling | Nonlinear Functional Analysis. 1st Edition[END_REF] for the theory) and very close to the proof proposed in [START_REF] Latché | A convergent staggered scheme for the variable density incompressible Navier-Stokes equations[END_REF]. Note that the system is nonlinear and that uniqueness is not guaranteed. In [START_REF] Nicolaides | Analysis and convergence of the MAC scheme. II. Navier-Stokes equations[END_REF], uniqueness is shown for a covolume approach of the MAC scheme using the total pressure form of the constant density incompressible Navier-Stokes equations on a two-dimensional uniform grid, assuming H 2 regularity of the solution and under small data conditions. Here we are concerned with the general 2D or 3D case, without any condition on the data nor regularity of the solution.

Theorem 4.3 (Existence of a solution). For a given n ∈ [ 1, N -1 ], let us assume that the density ρ n is such that 0 < ρ min ≤ ρ n K ≤ ρ max for all K ∈ M. Then the non-linear system (9)-( 10) admits at least one solution (ρ n+1 , u n+1 , p n+1 ), and any possible solution satisfies the estimates (31) and (36).

Proof. Let N M = card(M), N E = card(E int ) and let V = R N M × R N E × R N M . We introduce the function F : [0, 1] × V → V defined by: F (λ, (ρ K ) K∈M , (u σ ) σ∈Eint , (p K ) K∈M ) = (A K ) K∈M , (B σ ) σ∈Eint , (C K ) K∈M , with A K = 1 δt (ρ K -ρ n K ) + λ 1 |K| σ∈E(K) F K,σ , K ∈ M, B σ = 1 δt (ρ Dσ u σ -ρ n Dσ u n σ ) + λ 1 |D σ | ∈ E(Dσ) F σ, u -div µ λ D(u) σ + (∇p) σ -(f λ ) σ , σ ∈ E int , C K = - 1 |K| σ∈E(K) |σ| u K,σ + 1 |K| L∈M |L| p L , K ∈ M,
where div µ λ D(u) σ and (f λ ) σ are obtained by applying their definition in the scheme with (λ ρ K + (1 -λ) ρ n K ) K∈M instead of (ρ K ) K∈M . The function F is continuous from [0, 1]×V to V and the problem F (1, •) = 0 is equivalent to System [START_REF] Eymard | The finite volume method. Handbook for Numerical Analysis[END_REF]. Indeed, the first and second lines correspond to the discrete mass and momentum balance equations and, multiplying the third line by |K| and summing over the cells, we end up, by conservativity, with card(M) L∈M |L| p L = 0, that is Equation [START_REF] Fathy | A fourth order difference scheme for the Maxwell equations on Yee grid[END_REF]. Removing this term from C K , we recover the divergence-free constraint (9c).

It is clear that F (λ, •) = 0 still implies, for any λ ∈ [0, 1], the uniform bounds (31) for the density (A K = 0, ∀K ∈ M still corresponds to the discretization of a convection equation for ρ with a velocity scaled by λ, and thus still divergence-free). In addition, thanks to the linearity of the expression of the dual mass fluxes as a function of the primal ones, the structure which allows to obtain the estimate (35) is preserved when λ varies in the [0, 1] interval, so the unknowns of the velocity are bounded independently from λ. Finally, once both the density and the velocity are bounded independently of λ and knowing that the mean value of the pressure over Ω is zero, thanks to the so-called inf-sup condition satisfied by the MAC scheme, the system B σ = 0, ∀σ ∈ E int , yields a bound for the pressure unknowns (in mesh dependent norms), the expression of which may be found in [START_REF] Batteux | Modélisation numérique d'écoulements incompressibles diphasiques : Applications à l'interaction fluide-structure[END_REF]. Therefore, there exists a real number M such that the boundary of the closed ball of radius M , B M ⊂ V , does not contain any point such that

F (λ, •) = 0, ∀λ ∈ [0, 1].
The problem F (0, •) = 0 has a unique solution: the density is fixed by the first block of equations A K = 0, ∀K ∈ M, and the system B σ = 0, ∀σ ∈ E int and C K = 0, ∀K ∈ M is the discretization of the generalized Stokes system for the velocity and the pressure, with now a fixed viscosity and right-hand side. In addition, F (0, •) is regular (note that both µ λ and f λ do not depend on ρ for λ = 0), the Jacobian of F (0, •) is triangular per blocks, and the diagonal blocks are associated to the identity for ρ and the generalized Stokes problem for u and p (see [START_REF] Latché | A convergent staggered scheme for the variable density incompressible Navier-Stokes equations[END_REF][START_REF] Batteux | Modélisation numérique d'écoulements incompressibles diphasiques : Applications à l'interaction fluide-structure[END_REF]); therefore, the Jacobian determinant is not equal to 0. The topological degree of F (0, •) with respect to 0 ∈ V and B M is thus different from 0. Since it is preserved up to λ = 1 (because, as said above, F does not vanish on the boundary of B M ), we obtain that there exists at least one solution to F (1, •) = 0, i.e. to the scheme.

Convergence of the scheme

We prove in this section the main result of this paper, namely the convergence, up to the extraction of a subsequence, of a sequence of discrete solutions obtained with a sequence of time steps (δt (m) ) m∈N and of meshes (M (m) ) m∈N with both the time and space steps tending to zero. We begin with some estimates on the time translates of the velocity (Section 5.1) which will be crucial to prove the compactness of the sequence of discrete velocities; the convergence theorem is then stated and proved in Section 5.2.

5.1.

Estimate on the time translates of the velocity. We derive in this section an estimate on the time translates of the velocity in L 2 (Ω × (0, T )) which will be used to conclude to the compactness of a sequence of solutions. Lemma 5.1. Let (ρ, u, p) be a discrete solution of [START_REF] Eymard | The finite volume method. Handbook for Numerical Analysis[END_REF]. Let θ > 0 be such that the parameter θ M measuring the regularity of the mesh satisfies θ M ≥ θ, and let τ be a positive real number lower than T . Then, the following estimate holds:

(37) T -τ 0 Ω u(x, t + τ ) -u(x, t) 2 dx dt ≤ C (τ + δt + √ τ + δt)
where C > 0 only depends on T , the L 2 -norm of f , ρ min , ρ max , µ max and θ.

t n t n+l0 t n+l0+1 τ α δt
Proof. Since ρ ≥ ρ min , a bound in L 2 (Ω × (0, T )) d of √ ρ (u(•, • + τ ) -u(•, •)) readily yields a bound in the same norm of u(•, • + τ ) -u(•, •)). For t ∈ (0, T ) and x ∈ Ω, let us write ρ(x, t) u(x, t + τ ) -u(x, t) = T 1 (x, t) -T 2 (x, t), with T 1 (x, t) = ρ(x, t + τ ) u(x, t + τ ) -ρ(x, t) u(x, t), (38)
T 2 (x, t) = ρ(x, t + τ ) -ρ(x, t) u(x, t + τ ). ( 39
)
Let τ < T , t ∈ (0, T -τ ) and let us define the integer numbers n and by n = t/δt and n + = (t + τ )/δt , where, for s ∈ R, s is the integer number such that s < s ≤ s + 1. We thus have δt ≤ τ + δt, t n < t ≤ t n+1 and t n+ < t + τ ≤ t n+ +1 . For x ∈ Ω, we have

T 1 (x, t) = ρ n+ +1 (x) u n+ +1 (x) -ρ n+1 (x) u n+1 (x) = k=1 ρ n+k+1 (x) u n+k+1 (x) -ρ n+k (x) u n+k (x) .
Using the discrete momentum balance equation (9b), we thus get, for i ∈ [ 1, d ] and

x ∈ D σ , σ ∈ E (i) , T 1,i (x, t) = +1 k=2 div(µ D(u)) n+k σ -div(ρu i u) n+k σ -(∇p) n+k σ + f n+k σ , with T 1,i the i-th component of T . Let A(t) = (A 1 (t), . . . A d (t)) t , with A i (t) defined by A i (t) = Ω T 1,i (x, t) u i (x, t + τ ) -u i (x, t) dx, so that T -τ 0 Ω ρ(x, t + τ ) u(x, t + τ ) -ρ(x, t) u(x, t) • u(x, t + τ ) -u(x, t) dx dt = d i=1 T -τ 0 A i (t) dt. We split A i (t) as A i (t) = A i,c (t) + A i,d (t) + A i,p (t) + A i,f (t) with A i,d (t) = δt +1 k=2 σ∈E (i) |D σ | div(µ D(u)) n+k σ u n+ +1 σ -u n+1 σ , ( 40
) A i,c (t) = -δt +1 k=2 σ∈E (i) |D σ | div(ρu i u) n+k σ u n+ +1 σ -u n+1 σ , ( 41
) A i,p (t) = -δt +1 k=2 σ∈E (i) |D σ | (∇p) n+k σ u n+ +1 σ -u n+1 σ , ( 42
) A i,f (t) = δt +1 k=2 σ∈E (i) |D σ | f n+k σ u n+ +1 σ -u n+1 σ . ( 43
)
Thanks to the fact that both u n+1 and u n+ +1 are divergence free, we have d i=1 A i,p (t) = 0. The following of the computation consists in using bounds for each of the other terms at a given time t, derived from the stability properties of the discrete operators or by standard estimates, and then conclude by integration over the time. To provide a guideline for this computation, let us first write an analogue at the continuous level. To this purpose, we set each of these terms under the following form Let us suppose that the function (s, t) → v(s)w(t) lies in L 1 (]0, T [ 2 ), for instance because both v and w lie in L 1 (0, T ). Then, if v and w lie in L 2 (0, T ),

A(t) = v(t)
(45) |I| = T -τ 0 v(t) τ 0 w(t + s) ds dt = τ 0 T -τ 0 v(t) w(t + s) dt ds ≤ τ 0 T -τ 0 v(t) 2 dt 1/2 T -τ 0 w(t + s) 2 dt 1/2 ds.
Each integral over the time is bounded by the integral over (0, T ), and we get,

|I| ≤ τ v L 2 (0,T ) w L 2 (0,T ) .
Diffusion term -Let us reproduce this computation at the discrete level for the diffusion term. Lemma 3.2 yields

|A i,d (t)| ≤ µ max δt +1 k=2 |u n+k i | 1,E |u n+ +1 i | 1,E + |u n+1 i | 1,E .
Let α ∈ [0, 1] be defined by α δt = t n + τ -t n+l and let 0 = τ /δt (see Figure 6). Then, for t ∈ (t n , t n+1 -αδt), = 0 while, for t ∈ (t n+1 -αδt, t n+1 ), = 0 + 1.

We thus have ( 46)

T -τ 0 A i,d (t) dt ≤ µ max (1 -α) δt 2 N -0-1 n=0 0+1 k=2 |u n+k i | 1,E |u n+ 0+1 i | 1,E + |u n+1 i | 1,E + µ max α δt 2 N -l0-2 n=0 0 +2 k=2 |u n+k i | 1,E |u n+ 0+2 i | 1,E + |u n+1 i | 1,E .
Let A 1 be the first summation at the right-hand side of this equation. Then we get, reordering the summations:

A 1 ≤ µ max (1 -α) 0+1 k=2 δt N -0-1 n=0 δt |u n+k i | 1,E |u n+ 0 +1 i | 1,E + |u n+1 i | 1,E .
Since, thanks to the Cauchy-Schwarz inequality, for k ∈ [ 2, 0 + 1 ],

N -l0-1 n=0 δt |u n+k i | 1,E |u n+ 0+1 i | 1,E + |u n+1 i | 1,E ≤ 2 N n=1 δt |u n i | 2 1,E = 2 |u i | 2 L 2 (H 1 E ) , we obtain that A 1 ≤ 2 µ max (1 -α) 0 δt |u i | 2 L 2 (H 1 E )
. The same arguments yields, for the second summation of the right-hand side of Inequality (46):

A 2 ≤ 2 µ max α ( 0 + 1) δt |u i | 2 L 2 (H 1 E )
. Finally, we get, since ( 0 + 1)δt ≤ τ + δt,

(47) T -τ 0 A i,d (t) dt ≤ A 1 + A 2 ≤ 2 µ max |u i | 2 L 2 (H 1 E ) (τ + δt).
Forcing term -For the term associated to the forcing term, we write

|A i,f (t)| = δt +1 k=2 f n+k i L 2 (Ω) u n+ +1 i L 2 (Ω) + u n+1 i L 2 (Ω) ,
and thus, by the same technique as for the diffusion term,

T -τ 0 A i,f (t) dt ≤ 2 τ u i L ∞ (0,T ;L 2 (Ω)) f i L 2 (Ω×(0,T )) . (48)
Convection term -For the convection term, returning to the notations of Equation (44), the technique to obtain the estimate (45) is ineffective.The function

(s, t) → v(s)w(t) lies in L 1 (]0, T [ 2 ), v ∈ L 2 (0, T ) but w ∈ L 2 (0, T ).
The continuous analogue of the computation that we implement at the discrete level is the following one. The integral to estimate is of the form

I = T -τ 0 v(t) t+τ t w(s) ds dt,
with v ∈ L 2 (0, T ) and w ∈ L 4/3 (0, T ). We first change the order of the integrations, to obtain

I = T τ w(s) s s-τ v(t) dt ds,
The Cauchy-Schwarz inequality yields

|I| ≤ T τ |w(s)| s s-τ v(t) 2 dt 1/2 τ 1/2 ds.
Then the integral over (s -τ, s) is bounded by the integral over (0, T ) to get

|I| ≤ τ 1/2 v L 2 (0,T ) T τ |w(s)| ds,
and the Holder inequality yields:

|I| ≤ τ 1/2 v L 2 (0,T ) T 1/4 w L 4/3 (0,T ) .
Let us return to the discrete level. Lemma 3.3 yields

A i,c (t) ≤ C +1 k=2 δt u n+k 2 L 4 (Ω) |u n+ +1 i | 1,E + |u n+1 i | 1,E ,
with C only depending on θ and ρ max . We thus have

T -τ 0 A i,c (t) dt ≤ C(A 1 + A 2 ) with A 1 = N -0 n=1 δt |u n i | 1,E min(n+ 0+1, N ) k=n+1 δt u k 2 L 4 (Ω) . A 2 = N n= 0+1 δt |u n i | 1,E n k=n-0 δt u k 2 L 4 (Ω) .
Reordering the sums, we get

A 1 ≤ N k=2 δt u k 2 L 4 (Ω) S n 1 , with S n 1 = min(k-1,N ) n=max(k-0-1,1) |u n i | 1,E .
Using the fact that ( 0 + 1) δt ≤ τ + δt, the Cauchy-Schwarz inequality yields:

S n 1 ≤ 0+1 k=1 δt 1/2 N +1 k=1 δt |u n i | 2 1,E 1/2 ≤ (τ + δt) 1/2 u i L 2 (H 1 E ) .
For the term A 2 , we get

A 2 ≤ N k=1 δt u k 2 L 4 (Ω) S n 2 , with S n 2 = min(k+ 0,N ) n=max( 0+1,k) |u n i | 1,E ,
and S n 2 satisfies the same bounds as S 1 , n. On the other hand, the Cauchy-Schwarz inequality yields, for any function

u ∈ L 6 (Ω), u L 4 (Ω) ≤ u 1/4 L 2 (Ω) u 3/4 L 6 (Ω) , and u L 6 (Ω) ≤ |u| 1,E by discrete Sobolev inequalities, so T -τ 0 A i,c (t) dt ≤ 2 C (τ + δt) 1/2 u i L 2 (H 1 E ) u 1/2 L ∞ (0,T ;L 2 (Ω)) N +1 n=3 δt |u n | 3/2 1,E .
Finally, using the Holder inequality with q = 4/3 and q = 1/4 to estimate the last sum,

T -τ 0 A i,c (t) dt ≤ 2 C T 1/4 (τ + δt) 1/2 u i L 2 (H 1 E ) u 1/2 L ∞ (0,T ;L 2 (Ω)) u i 3/2 L 2 (0,T ;H 1 E ) (49)
. Term T 2 -Let us now estimate the term issued from the quantity T 2 of equation (39). By the same developments as previously, we have to estimate the integral over the time interval of

B i (t), i ∈ [ 1, d ], with B i (t) = k=1 σ∈E (i) int |D σ | u n+ +1 σ (ρ n+k+1 Dσ -ρ n+k Dσ ) (u n+ +1 σ -u n+1 σ ),
where n and depend on t, with the same definition as previously. Thanks to the mass balance over the dual cells, we get:

B i (t) = δt +1 k=2 σ∈E (i) int |D σ | u n+ +1 σ (u n+ +1 σ -u n+1 σ ) ∈ E(Dσ) F n+k σ,
.

Thanks to Lemma 3.3,

B i (t) ≤ C |u n+ +1 i | 1,E + |u n+1 i | 1,E |u n+ +1 σ | 1,E +1 k=2 δt |u n+k | 1,E ,
with C only depending on θ and ρ max . We proceed in a way similar to the derivation of the estimate for the convection term. First, we take benefit to the fact that the sum over k (in the continuous setting, the integral over (t, t+τ )) involves a function lying in L 2 (0, T ) (without needing, here, a switch of the integration order thanks to a discrete Fubini technique) to make appear, thanks to the Cauchy-Schwarz inequality, the factor (τ + δt) 1/2 :

B i (t) ≤ C |u n+ +1 i | 1,E + |u n+1 i | 1,E |u n+ +1 σ | 1,E +1 k=2 δt 1/2 +1 k=2 |u n+k | 2 1,E 1/2 ≤ C |u n+ +1 i | 1,E + |u n+1 i | 1,E |u n+ +1 σ | 1,E (τ + δt) 1/2 u L 2 (H 1 E )
. So, integrating with respect to time yields:

(50) T -τ 0 B i (t) dt ≤ C (τ + δt) 1/2 u i 2 L 2 (H 1 E ) u L 2 (H 1 E ) .
Conclusion -The conclusion follows by summing Equations ( 47), ( 48), ( 49) and (50) over i ∈ [ 1, d ] and then summing once again the obtained relations.

5.2.

Convergence of the discrete solutions. We are now in position to state and prove the convergence result which is the aim of this paper.

Theorem 5.2 (Convergence of the discrete solutions). Let (δt (m) ) m∈N be a sequence of time steps and (M (m) ) m∈N be a sequence of MAC grids (in the sense of Definition 2.1) such that δt (m) → 0 and h M (m) → 0 as m → +∞. We assume that there exists θ > 0 controlling the regularity of every MAC mesh in the sequence: θ M (m) ≤ θ for any m ∈ N, with θ M (m) defined by [START_REF] Liu | Projection method iii. spatial discretization on the staggered grid[END_REF]. For m ∈ N, let (ρ (m) , u (m) ) be a solution to [START_REF] Eymard | The finite volume method. Handbook for Numerical Analysis[END_REF] for δt = δt (m) and M = M (m) . Then there exists ρ with ρ min ≤ ρ(x, t) ≤ ρ max for a.e. x ∈ Ω and t ∈ (0, T ) and ū ∈ L ∞ (0, T ; L 2 (Ω) d ) ∩ L 2 (0, T ; E(Ω)) satisfying, up to a subsequence: -the sequence (u (m) ) m∈N converges to ū in L 2 (Ω × (0, T )) d , -the sequence (ρ (m) ) m∈N converges to ρ in L q (Ω × (0, T )), for q < +∞, -(ρ, ū) is a solution to the weak formulation given by Definition 1.1.

Proof.

Step 1-Weak convergence of (ρ (m) ) m∈N in L q (Ω × (0, T )), q ≥ 1, and maximum principle -Thanks to the estimate (31), the sequence (ρ (m) ) m∈N satisfies ρ min ≤ ρ (m) (x, t) ≤ ρ max for a.e. x ∈ Ω and t ∈ (0, T ), ∀m ∈ N, which implies that, up to the extraction of a subsequence, this sequence weakly converges to a function ρ ∈ L ∞ (Ω × (0, T )) satisfying the same bounds.

Step 2-Compactness of (u (m) ) m∈N in L 2 (Ω × (0, T )) d , regularity of the limit and divergence free constraint -We apply Theorem C.1 with B = L 2 (Ω × (0, T )) d , X (m) the space of the discrete velocities endowed with the discrete H 1 -norm defined by ( 20) and (f (m) ) m∈N = (u (m) ) m∈N . Thanks to Estimate (36), we obtain both the boundedness of (u (m) ) m∈N in B and the second assumption of the theorem. Lemma 5.1 yields the third assumption, with η(τ

) = C (τ + √ 2 τ
), where C is the constant of the inequality (37), which only depends on T , ρ min , ρ max , µ max , the L 2 -norm of f (which is either a data or controlled by T and ρ max ) and θ, and with M (τ ) such that, for m ≥ M (τ ), δt (m) ≤ τ . Hence, there exists ū ∈ L 2 (Ω × (0, T )) d such that, up to a subsequence,

u (m) → ū in L 2 (Ω × (0, T )) d as m → +∞.
In addition, the uniform bound of the sequence (u (m) ) m∈N in L 2 (0, T ; H 1 E ) yields that ū ∈ L 2 (0, T ; H 1 0 (Ω) d ) (see [START_REF] Eymard | The finite volume method. Handbook for Numerical Analysis[END_REF]Theorem 14.2]). We now prove that ū is divergence free by passing to the limit in the discrete mass balance equation. Let ϕ ∈ C ∞ c (Ω × (0, T )) and, for m ∈ N, K ∈ M (m) and n ∈ [ 0, N (m) -1 ], let ϕ n K be the mean value of ϕ over K × (t n , t n+1 ). Multiplying the divergence constraint by δt (m) |K| ϕ n K and summing over the cells and the time steps, we get, for m ∈ N,

I (m) = N (m) -1 n=0 δt (m) K∈M (m) ϕ n K d i=1 σ∈E(K)∩E (m,i) int |σ| u σ e (i) • n K,σ = 0,
where, for short, we have denoted by

E (m,i) int the set (E (m)
int ) (i) . Reordering the sums, we get

I (m) = N (m) -1 n=0 δt (m) d i=1 σ∈E (m,i) int , σ=K|L |D σ | u σ ϕ n K -ϕ n L d σ e (i) • n K,σ = 0, with d σ = |D σ |/|σ|.
This last expression may be seen as the integral over Ω × (0, T ) of the inner product of u (m) with a piecewise constant vector-valued function, the ith component of which takes the value e (i)

•n K,σ (ϕ n K -ϕ n L )/d σ over D σ ×(t n , t n+1 ), σ ∈ E (m,i) int
and n ∈ [ 0, N (m) -1 ]. Thanks to the regularity of ϕ, the latter function converges to -∇ϕ in L ∞ (Ω × (0, T )) d , and we get, invoking the convergence of (u (m) ) m∈N , lim m→+∞

I (m) = - T 0 Ω ū • ∇ϕ dx dt = 0,
which implies div ū = 0 a.e. in Ω × (0, T ).

Step 3-Passing to the limit in the mass balance equation -We show in this step that the limit (ρ, ū) obtained in the previous steps satisfies the weak mass balance equation [START_REF] Deimling | Nonlinear Functional Analysis. 1st Edition[END_REF]. = 0 with

Let ϕ ∈ C ∞ c ([0, T ) × Ω), and m ∈ N. For K ∈ M (m) and n ∈ [ 0, N (m) ], let ϕ n K = ϕ(x K , t n ),
T (m) 1 = N (m) -1 n=0 K∈M (m) |K| (ρ n+1 K -ρ n K ) ϕ n K , T (m) 2 = N (m) -1 n=0 δt (m) K∈M (m) ϕ n K σ∈E(K) F n+1 K,σ .
Rearranging the sum in T to perform a discrete integration by parts with respect to the time variable, we get, using the fact that ϕ vanishes at t = T :

T (m) 1 = - N (m) -1 n=0 δt (m) K∈M (m) |K| ρ n+1 K ϕ n+1 K -ϕ n K δt (m) - K∈M (m) |K| ρ 0 K ϕ 0 K .
We recognize in the first term the integral with respect to space and time of ρ (m) multiplied by a discrete time-derivative of ϕ which converges to ∂ t ϕ in L 1 (Ω×(0, T )) when the space and time step tends to zero. In the second term, when the space step tends to zero, thanks to the definition [START_REF] Gallouët | Equations aux dérivées partielles[END_REF] of the initial value of the scheme, the function K∈M (m) ρ 0 K X K (x) converges to ρ 0 in L q (Ω), q < +∞, and the function K∈M (m) ϕ 0 K X K (x) converges to ϕ in L ∞ (Ω). So passing to the limit, we obtain:

lim m→∞ T (m) 1 = - T 0 Ω ρ(t, x) ∂ t ϕ(t, x) dx dt - Ω ρ 0 (x) ϕ(0, x) dx.
Let us now focus on T (m) 2

. Using the expression of the mass flux F K,σ and the boundary conditions on the velocity, we reorder the sum in T (m) 2 so as to perform a discrete integration by parts with respect to the space variable:

T (m) 2 = - d i=1 N (m) -1 n=0 δt (m) σ∈E (i) int , σ=K|L |D K,σ | ρ n+1 K + |D L,σ | ρ n+1 L u n+1 σ (∂ϕ) n σ + R (m) 2 ,
with:

(∂ϕ) n σ = |σ| |D σ | (ϕ n L -ϕ n K ) (n K,σ • e (i) ), R (m) 2 = - d i=1 N (m) -1 n=0 δt (m) σ∈E (i) int , σ=K→L |D L,σ | ρ n+1 K -ρ n+1 L u n+1 σ (∂ϕ) n σ ,
where the notation σ = K → L means that K and L are chosen in such a way that u K,σ ≥ 0. The first term in T (m) 2

is the sum over i ∈ [ 1, d ] of the integral with respect to the time and space of the product of ρ (m) u (m) i by a function ð i ϕ defined for a.e. x ∈ Ω and t ∈ (0, T ) by

ð (m) i ϕ(x, t) = N (m) -1 n=0 δt (m) σ∈E (i) int (∂ϕ) n σ X Dσ X (tn,tn+1] . Since this function converges in L ∞ (Ω × (0, T )) to ∂ i ϕ, u (m)
i converges to ūi in L 2 (Ω × (0, T )) and ρ (m) converges weakly to ρ in L 2 (Ω × (0, T )), we have

lim m→∞ T (m) 2 -R (m) 2 = - T 0 Ω ρ(t, x) ū(t, x) • ∇ϕ(t, x) dx dt.
The Cauchy-Schwarz inequality yields for R (m) 2 :

|R (m) 2 | ≤ h (m) ∇ϕ L ∞ (Ω×(0,T )) d N (m) -1 n=0 δt (m) σ∈Eint, σ=K→L |D σ | |u n+1 σ | 1/2 N (m) -1 n=0 δt (m) σ∈Eint, σ=K→L |σ| ρ n+1 L -ρ n+1 K 2 |u n+1 σ | 1/2 so R (m) 2
tends to zero thanks to the weak BV-estimate on ρ stated in Lemma 4.1 and the control of u in L 2 (Ω × (0, T )) d (and therefore in L 1 (Ω × (0, T )) d ), which concludes this step.

Step 4 -Weak convergence of the velocity discrete derivatives in L 2 (Ω × (0, T )) -Let i, j ∈ [ 1, d ]. Owing to the "discrete L 2 (H 1 0 )" bound (36), the sequence of derivatives (ð j u (m) i

) m∈N (see Definition [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface[END_REF]) is bounded in L 2 (Ω × (0, T )), and is thus weakly convergent to some limit in L 2 (Ω × (0, T )). Let ϕ ∈ C ∞ c (Ω), let us denote by {D , E (m,i,j) } the set of (i, j)-gradient cells of the mesh M (m) , and, for ∈ E (m,i,j) and n ∈ [ 0, N (m) -1 ], let ϕ n be the value of ϕ at (x , (t n + t n+1 )/2), with x the mass center of . Then,

T 0 Ω ð j u (m) i ϕ dx dt = I (m) + R (m) , with
I (m) = N (m) -1 n=0 δt (m) ∈ E (m,i,j) |D | (ð j u i ) n+1 ϕ n = N (m) -1 n=0 δt (m) ∈ E (m,i,j) , =σ|σ | | (u n+1 σ -u n+1 σ ) ϕ n ,
where, by convention, we have supposed in the last sum that σ and σ are ordered in such a way that x σ,j > x σ ,j . Let us now use the notation σ = -j → +j to mean that -j and +j are the two faces of D σ normal to e (j) and that the j-th coordinate of the points of -j is lower than the j-th coordinate of the points of +j . By a reordering of the sums, we get

I (m) = - N (m) -1 n=0 δt (m) σ∈E (m,i) , σ= -j → +j | +j | u n+1 σ (ϕ n +j -ϕ n -j ).
Let x -j ,j and x +j ,j be the j-th coordinate of the points of -j and +j respectively. Then we remark that, for σ = -j → +j , (x +j ,j -x -j ,j )

| +j | = |D σ |, and
I (m) = - N (m) -1 n=0 δt (m) σ∈E (m,i) , σ= -j → +j |D σ | u n+1 σ ϕ n +j -ϕ n -j
x +j ,j -x -j ,j .

Thanks to the regularity of ϕ, the piecewise constant function equal to (ϕ n +jϕ n +j )/(x +j ,j -x -j ,j ) over D σ , σ ∈ E (m,i) converges to ∂ j ϕ in L ∞ (Ω × (0, T )), so, thanks to the convergence of (u

(m) ) m∈N in L 2 (Ω × (0, T )) d , lim m→+∞ I (m) = - T 0 Ω ūi ∂ j ϕ dx dt.
The residual term R (m) reads

R (m) = N (m) -1 n=0 δt (m) ∈ E (m,i,j) |D | (ð j u i ) n+1 (ϕ n D -ϕ n ),
with ϕ n D the mean value of ϕ over D × (t n , t n+1 ). Thanks to the regularity of ϕ and the control of (ð j u (m) i

) m∈N in L 2 (Ω × (0, T )), R (m) tends to zero when m tends to +∞, which concludes this point.

Step 5 -Strong convergence of (ρ (m) ) m∈N in L q (Ω × (0, T )), q ≥ 1 -We follow the proof in [START_REF] Latché | A convergent staggered scheme for the variable density incompressible Navier-Stokes equations[END_REF]. On one hand, a classical property of the weak topology in L 2 (Ω × (0, T )) yields, for the weak limit ρ of the sequence (ρ (m) ) m∈N :

ρ L 2 (Ω×(0,T )) ≤ lim inf m→∞ ρ (m)
L 2 (Ω×(0,T )) .

On the other hand, summing the discrete entropy estimate (32) over the cells and the first n time steps, we have for all n in [ 0, N -1 ]:

K∈M (m) |K| (ρ n+1 K ) 2 ≤ K∈M (m) |K| (ρ 0 K ) 2 ≤ ρ 0 2 L 2 (Ω) ,
and thus, integrating over the time, ρ (m) L 2 (Ω×(0,T )) ≤ T ρ 0 L 2 (Ω) . We know that the weak limit ρ satisfies the convection equation with a divergence-free velocity field ū lying in L 2 (0, T ; H 1 0 (Ω) d ) and the theory of renormalized solutions [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF] implies that ρ L 2 (Ω×(0,T )) = T ρ 0 L 2 (Ω) . We thus have

lim m→+∞ ρ (m) L 2 (Ω×(0,T )) = ρ L 2 (Ω×(0,T )) ,
which implies the convergence of the sequence (ρ (m) ) m∈N to ρ in L 2 (Ω × (0, T )).

Finally, this convergence also holds in L q (Ω × (0, T )) for any q < +∞, since the sequence (ρ (m) ) m∈N is bounded in L ∞ (Ω × (0, T )).

Step 6 -Strong convergence of the viscosity and the forcing term -Thanks to the strong convergence of (ρ (m) ) m∈N in L q (Ω × (0, T )), q ≥ 1, there exists a subsequence still denoted (ρ (m) ) m∈N such that ρ (m) (t, x) → ρ(t, x) for a.e.

x in Ω and t in (0, T ). The Lebesgue dominated convergence theorem implies that, for any real function g continuous over [ρ min , ρ max ] (and so bounded), the sequence (g(ρ (m) )) m∈N converges to g(ρ) in L q (Ω × (0, T )). This yields the convergence of the sequence (µ(ρ (m) )) m∈N defined by

µ (m) (x, t) = N (m) -1 n=0 K∈M (m) µ(ρ n+1 K ) X K (x) X (tn,tn+1] (t).
Invoking Lemma B.2, the coefficients µ (i,j) involved in the definition of the diffusion term therefore also converge to µ(ρ) in L q (Ω × (0, T )), q ≥ 1. In the case where the forcing term is defined as a function of ρ, the same arguments prove its convergence to f (ρ) in L q (Ω × (0, T )) d , q ≥ 1.

Step 7 -Passing to the limit in the momentum balance equation

-Let ϕ ∈ C ∞ c ([0, T ) × Ω) d , such that divϕ = 0, and, for m ∈ N, n ∈ [ 0, N (m) ], i ∈ [ 1, d ] and σ ∈ E (m,i) int let us define by ϕ n σ by ϕ n σ = 1 |σ| σ ϕ(x, t n ) • e (i) dγ(x).
Multiplying the momentum balance equation (9b) by δt

(m) |D σ | ϕ n σ , summing over σ ∈ E (m,i) int , over i ∈ [ 1, d ] and n ∈ [ 0, N (m) -1 ], we get T (m) 1 + T (m) 2 + T (m) 3 + T (m) 4 + T (m) 5 = 0, with T (m) 1 = d i=1 N (m) -1 n=0 δt (m) σ∈E (m,i) int |D σ | ð t (ρu i ) n+1 σ ϕ n σ , T (m) 2 = d i=1 N (m) -1 n=0 δt (m) σ∈E (m,i) int |D σ | div(ρ u i u) n+1 σ ϕ n σ , T (m) 3 = - d i=1 N (m) -1 n=0 δt (m) σ∈E (m,i) int |D σ | div(µ D(u)) n+1 σ ϕ n σ , T (m) 4 = d i=1 N (m) -1 n=0 δt (m) σ∈E (m,i) int |D σ | (∇p) n+1 σ ϕ n σ , T (m) 5 = d i=1 N (m) -1 n=0 δt (m) σ∈E (m,i) int |D σ | f n+1 σ ϕ n σ .
Thanks to the definition of ϕ n σ and to the L 2 -duality between the discrete gradient and the discrete divergence, the term T , we obtain, using the expression of the discrete time derivative,

T (m) 1 = - d i=1 N (m) -1 n=0 σ∈E (m,i) int δt (m) |D σ | ρ n+1 Dσ u n+1 σ ϕ n+1 σ -ϕ n σ δt (m) - d i=1 σ∈E (m,i) int |D σ | ρ 0 Dσ u 0 σ ϕ 0 σ .
Thanks to Lemma B.2 and to the convergence of ((ρ (m) ) m∈N ) to ρ in L 2 (Ω×(0, T )), the piecewise constant function equal to ρ n+1 Dσ over D σ × (t n , t n+1 ), for σ ∈ E (m,i) int converges to ρ in L 2 (Ω×(0, T )) when m tends to +∞, ∀i ∈ [ 1, d ] (remember that the face density is a convex combination of the density in the two adjacent cells). The regularity of ϕ yields that the piecewise function taking the value (ϕ n+1 σ -ϕ n σ )/δt (m) on the same space-time domains converges to

∂ t ϕ i in L ∞ (Ω × (0, T )), ∀i ∈ [ 1, d ].
Finally, still by Lemma B.2, the function of the space variable taking the value

ρ 0 Dσ over D σ , for σ ∈ E (m,i) int converges to ρ 0 in L 2 (Ω), ∀i ∈ [ 1, d ],
since the initial condition for ρ does so. The convergence of the velocity in L 2 (Ω × (0, T )) d was proven in the previous step, and the convergence of the discrete initial value for the velocity to u 0 in L 2 (Ω) d is standard, as is the convergence of the function equal to

ϕ 0 σ over D σ , for σ ∈ E (m,i) int , to ϕ i (•, 0) in L ∞ (Ω), ∀i ∈ [ 1, d ]. Therefore: lim m→+∞ T (m) 1 = - T 0 Ω ū • ϕ dx dt - Ω ρ 0 u 0 • ϕ(x, 0) dx.
Let us now turn to the convection term T (m) 2

. By the same computation as in the proof of Lemma 3.3, we get, with the notations of Equation (28) and Figure 5,

T (m) 2 = d i=1 N (m) -1 n=0 δt (m) d j=1 ∈ E (m,i,j) , =σ|σ |D | ρ n+1 (û j ) n+1 u n+1 ϕ n σ -ϕ n σ x σ,j -x σ ,j
The Hölder inequality with q = 3/2 and q = 3 yields, for v ∈ L ∞ (Ω × (0, T )) (which holds for the components discrete velocities), at a.e. t ∈ (0, T ):

Ω |v(x, t)| 4/3+2 dx ≤ Ω |v(x, t)| 2 dx 2/3 Ω |v(x, t)| 6 dx 1/3 , so v 10/3 L 10/3 (Ω×(0,T )) ≤ v 2/3 L ∞ (0,T ;L 2 (Ω)) v 2 L 2 (0,T ;L 6 (Ω))
. We thus get that the sequence (u (m) ) m∈N is controlled in L 10/3 (Ω×(0, T )) d and thus converges to ū in L q (Ω × (0, T )) d , 1 ≤ q < 10/3. In addition, (ρ (m) ) m∈N converges to ρ in L q (Ω × (0, T )), q < +∞, and, thanks to the regularity of ϕ, the piecewise function

ð j ϕ (m) i equal to (ϕ n σ -ϕ n σ )/(x σ,j -x σ ,j ) over D × (t n , t n+1 ), ∈ E (m,i,j) and n ∈ [ 0, N (m) -1 ], converges to ∂ j ϕ i in L ∞ (Ω × (0, T )), for (i, j) ∈ [ 1, d ] 2 . The lemma B.2 thus yields: lim m→+∞ T (m) 2 = - d i=1 T 0 Ω ūi ū • ∇ϕ i dx dt.
By definition of the discrete diffusion term, we have:

T (m) 3 = T 0 Ω (µ D) E (m) (u (m) (x, t)) : ∇ϕ (m) (x, t) dx dt.
Thanks to the weak convergence of the discrete gradients of u (m) to ∇ū in L 2 (Ω × (0, T )) d×d , the convergence of µ(ρ (m) ) to µ(ρ) in L q (Ω × (0, T )), q ≥ 1 and the convergence of the discrete gradients of ϕ to ∇ϕ in L ∞ (Ω × (0, T )) d×d , we have: , with i = j, let u and v be two discrete fields corresponding to a discrete i-th velocity component and j-th velocity component, respectively, and let the partial derivatives of the discrete velocities be defined by [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface[END_REF]. We suppose that the normal velocities vanish on the boundary. Then we have:

(51) Ω ð j u ð i v dx = Ω ð i u ð j v dx.
Proof. Let us begin with the two-dimensional case and, to fix the ideas, let us suppose that i = 1 and j = 2. Both ð 2 u and ð 1 v are defined on the same twicestaggered cells, i.e. cells associated to the vertices of the mesh; we denote by E

(1,2) int the set of the mesh vertices which lie in the interior of the domain. Let us introduce the following local notations:

-for a ∈ E

(1,2) int , we define by u a N , u a S , v a W and v a E the four neighbouring discrete velocity unknowns (see Figure 7-(a)); -for σ ∈ E (1) , we define by v σ N W , v σ SW , v σ N E and v σ SW the four neighbouring unknowns for the field v (see Figure 7-(b)); -for K ∈ M, we define by u K W , u K E , v K N and v K S the four normal velocities on the edges of K (see Figure 7-(c)). Since the normal velocity vanishes at the boundary, either ð 2 u or ð 1 v vanishes, and the integral of Equation (51) reads:

Ω ð j u ð i v dx = a∈E (1,2) int ∪ E (i) rec (u a N -u a S ) (v a E -v a W ).
Reordering twice the sums, we get:

a∈E (1,2) int ∪ E (i) rec (u a N -u a S ) (v a E -v a W ) = σ∈E (1) int u σ (v σ SE -v σ SW ) -(v σ N E -v σ N W ) = K∈M (u K E -u K W ) (v K N -v K S ),
which concludes the proof in the two-dimensional case. These arguments readily extends to the three dimensional case, replacing (still with i = 1 and j = 2) the vertices of the mesh by the vertical edges. We are now in position to prove the following result, which is the discrete counterpart of the so-called Korn lemma.

a v a W v a E u a N u a S (a) σ v σ N W v σ N E v σ SW v σ SE (b) K v K N v K S u K W u K E (c)
Lemma A.2 (Discrete Korn Lemma). For any discrete velocity field u, the following identity holds:

D E (u) : D E (u) dx = 1 2 ∇ E u : ∇ E u dx + 1 2 div M (u) 2 dx.
Proof. By the definition of the strain rate tensor D and the simple algebraic identity T : S = 1 2 (T + T t ) : S, valid for two tensors T and S as soon as S is symmetric, we have:

D E (u) : D E (u) dx = 1 2 ∇ E u : (∇ E u + ∇ T E u) dx = 1 2 ∇ E u : ∇ E u dx + T,
where

T = 1 2 Ω n i,j=1
ð j u i ð i u j dx.

Lemma A.1 yields

T = 1 2 Ω n i,j=1
ð i u i ð j u j dx = sup{|x -y|, x ∈ Q, y ∈ ∪ P ∈N Q P }.

We assume that δ (m) tends to zero when m tends to +∞. In addition, we suppose that there exists n ∈ N and C ∈ R such that, for any m ∈ N,

card(N Q ) ≤ n, ∀Q ∈ Q (m)
, and

Q∈N P |Q| ≤ C |P |, ∀P ∈ P (m) .
Then, if a sequence (u (m) ) m∈N of piecewise functions over (P (m) ) m∈N converges in L r (Ω), r ≥ 1, to ū, then the sequence (R (m) (u (m) )) m∈N also converges to ū in L r (Ω).

Proof. Let > 0 and ϕ ∈ C ∞ c (Ω) such that ū -ϕ L r (Ω) ≤ . For P ∈ P (m) , let ūP and ϕ P be the mean value of ū and ϕ over P , respectively, and let us denote ū(m) and ϕ (m) the corresponding piecewise constant functions. We write R (m) (u (m) ) -ū L r (Ω) = T = ϕ -ū L r (Ω) .

We have u (m) -ū(m) L r (Ω) ≤ u (m) -ū L r (Ω) , since u (m) -ū(m) is the projection of u (m) -ū over the space of piecewise constant functions over the elements of P (m) ; for the same reason, ū(m) -ϕ (m) L r (Ω) ≤ ū -ϕ L r (Ω) . Since Lemma B.1 yields that the norm of the linear operator R (m) is lower than R = n (r-1)/r C 1/r , there exists M 1 so that, for m ≥ M 1 , both T compares ϕ to a piecewise constant function taking in Q ∈ Q (m) a value obtained by a convex combination of the averages of ϕ in cells of P (m) included in a ball δ (m) containing Q; therefore, thanks to the regularity of ϕ, this term tends to zero when m tends to +∞, so there exists M 2 such that, for m ≥ M 2 , T (m) 3

≤ . Combining the obtained estimates, we have R (m) (u (m) ) -ū L r (Ω) ≤ 2 (R + 1) for m ≥ max(M 1 , M 2 ), which concludes the proof.

The following remark clarifies the way these results are used in the present paper.

Remark B.1. Let us assume that the sequence (Q (m) ) m∈N is quasi-uniform with respect to the sequence (P (m) ) m∈N in the sense that there exists θ > 0 such that: max{|Q|, Q ∈ Q (m) } ≤ θ min{|P |, P ∈ P (m) }, ∀m ∈ N.

In addition, let us suppose that (52) holds and that there exists n ∈ N such that, for all m ∈ N, card(N Q ) ≤ n, ∀Q ∈ Q (m) and card(N P ) ≤ n, ∀P ∈ P (m) . and that Assumption (53) holds. This situation corresponds to the case where the mesh step for (P (m) ) m∈N and (Q (m) ) m∈N tends to zero and, as usual, the reconstruction operator is defined with a constant stencil.

Then the assumptions of Lemmas B.1 and B.2 are satisfied.

  2 ρ(x, t) |ū(x, t)| 2 dx + Ω µ(ρ) D(ū(x, t)) : D(ū(x, t)) dx = 0.

  n+1 K ) otherwise. Note that, thanks to the definition of D , |D | = K∈M, K∩D =∅ |K| 4 in the latter case.

Figure 4 .

 4 Figure 4. In the two-dimensional case, primal mass fluxes involved in the computation of the dual mass flux associated to = σ|σ ∈ E

Figure 5 .

 5 Figure 5. Velocities and velocity partial derivative over D involved in the right-hand side of Equation (28), in the twodimensional case, for i = 1. Left: j = 1. Right: j = 2.

4. 1 .

 1 Definition and estimates of the discrete solutions. To any family of real numbers (ρ n K ) K∈M, n∈[ 0,N ] and (p n K ) K∈M, n∈[ 0,N ] , we associate the discrete function ρ and p representing the density and the pressure, respectively, defined by (29)

Figure 6 .

 6 Figure 6. Notations for the velocity time translates estimates.

  ds, where the quantities v and w stand for (possibly discrete) norms of the solution.

  let us multiply the discrete mass balance (9a) by δt (m) |K| ϕ n K and sum over the time steps, to obtain T

f

  ) D(ū) : ∇ϕ dx dt. Finally, f is either fixed or converges to f (ρ), so, invoking once again the regularity of the test function: • ϕ dx dt. This concludes the convergence proof. Appendix A. The discrete Korn lemma Lemma A.1. Let i, j ∈ [ 1, d ] 2

Figure 7 .

 7 Figure 7. Local notations for the proof of the discrete Korn lemma.

ð

  j u j dx, which concludes the proof.Lemma B.2.Let Ω be a domain of R d , d ∈[START_REF] Batteux | Modélisation numérique d'écoulements incompressibles diphasiques : Applications à l'interaction fluide-structure[END_REF][START_REF] Boyer | Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models[END_REF], and let (P (m) ) m∈N and (Q (m) ) m∈N be two sequences of finite partitions of Ω. For m ∈ N, let the coefficients (γ P,Q ) (P,Q)∈(P (m) ×Q (m) ) satisfy the assumption (52), and let R (m) be the associated reconstruction operator. For P ∈ P and Q ∈ Q, we denote by N P and N Q the sets N P = {Q ∈ Q s.t. γ P,Q = 0} and N Q = {P ∈ P s.t. γ P,Q = 0}. Let δ (m) be defined by(53) δ (m) = max Q∈Q (m)

1 = 2 = 3 =

 123 R (m) (u (m) ) -R (m) (ū (m) ) L r (Ω) R (m) (ū (m) ) -R (m) (ϕ (m) ) L r (Ω) R (m) (ϕ (m) ) -ϕ L r (Ω) , T (m) 4

) 4 ≤

 4 lower than R . In addition, by definition of ϕ, T (m. Finally, the term term T (m) 3

Appendix B. On transfer operators from one mesh to another one This section gathers some stability and convergence results for transfer operators defining a piecewise constant function on a second mesh as a function of a piecewise constant function on a first one. Since the results of this section are quite general, we state them using generic notations; in this paper, they apply to transfers (or "reconstructions", in the sense defined below) between the primal mesh, a dual mesh or a mesh associated to the discrete velocity gradients. We first begin with a L p stability result. [START_REF] Boyer | Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models[END_REF], and let P and Q be two finite partitions of Ω. Let (u P ) P ∈P and (γ P,Q ) (P,Q)∈(P×Q) be two family of real numbers, and let us define the functions u and v, Ω → R, by

For P ∈ P and Q ∈ Q, we denote by N P and N Q the sets

and

Proof. Using the inequality (

which concludes the proof.

Let us now define what we mean by a reconstruction operator, then state and prove its convergence properties. [START_REF] Boyer | Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models[END_REF], and let P and Q be two finite partitions of Ω. Let (γ P,Q ) (P,Q)∈(P×Q) be a family of real numbers such that:

Let R the operator associating to a piecewise constant functions over the cells of P the following piecewise constant function over the cells of Q:

The operator R is referred to in the following as the reconstruction operator associated to the coefficients (γ P,Q ) (P,Q)∈(P×Q) .

Appendix C. A compactness result

Let us begin with a definition.

Definition C.1 (Compactly embedded sequence). Let B be a Banach space and (X (m) ) m∈N be a sequence of Banach spaces included in B. The sequence (X (m) ) m∈N is said to be compactly embedded in B if any sequence (u (m) ) m∈N satisfying:

Then the following compactness result is a corollary of [7, Proposition C5] (see also [START_REF] Gallouët | Equations aux dérivées partielles[END_REF]Chapter 4]).

Theorem C.1 (Time compactness with a sequence of subspaces). Let 1 ≤ p < ∞ and T > 0. Let B be a Banach space and (X (m) ) m∈N be a sequence of Banach spaces compactly embedded in B. Let (f (m) ) m∈N be a sequence of L p (0, T ; B) satisfying the following conditions:

(1) The sequence (f (m) ) m∈N is bounded in L p ((0, T ); B).

(2) The sequence ( f (m) L 1 (0,T ;X (m) ) ) m∈N is bounded. (3) There exists a non-decreasing function η from (0, T ) to R + such that lim τ →0 η(τ ) = 0 and, ∀τ ∈ (0, T ) we have:

for all m ≥ M (τ ), where M (τ ) is an integer possibly depending on τ (often, in applications, tending to +∞ when τ tends to zero).

Then, (f (m) ) m∈N is relatively compact in L p (0, T ; B).

Proof. Since f (m) ∈ L p (0, T ; B), one has

B dt → 0 as m → +∞, so that, by the third assumption of the theorem,

Therefore, [7, Proposition C5] applies. Email address: pascal.poullet@univ-antilles.fr