An extension of the MAC scheme to locally refined meshes: convergence analysis for the full tensor time-dependent Navier–Stokes equations

A variational formulation of the standard marker-and-cell scheme for the approximation of the Navier–Stokes problem yields an extension of the scheme to general 2D and 3D domains and more general meshes. An original discretization of the trilinear form of the nonlinear convection term is proposed; it is designed so as to vanish for discrete divergence free functions. This property allows us to give a mathematical proof of the convergence of the resulting approximate solutions, for the nonlinear Navier–Stokes equations in both steady-state and time-dependent regimes, without any small data condition. Numerical examples (analytical steady and time-dependent ones, inclined driven cavity) confirm the robustness and the accuracy of this method.


Introduction
The marker-and-cell (MAC) scheme, introduced in [14] is one of the most popular methods [20,26] for the approximation of the Navier-Stokes equations in the engineering framework, because of its simplicity and of its remarkable mathematical properties. In its standard implementation on a rectangular parallelepipedic mesh, the discrete unknowns are the normal component of the velocity at the faces of the mesh and the pressure at the center of the grid blocks. The mass conservation is discretely satisfied by a finite volume scheme in this grid. The discrete conservation of the component k = 1, . . . , d (where d = 2 or 3 is the dimension of space) of the momentum is written as a finite volume scheme in the staggered grid M (k) , made of rectangular parallelepipedic grid blocks, centered at the center of the faces whose normal vector is the basis vector e (k) .
The first error analysis seems to be that of [21] in the case of the time-dependent Stokes equations on uniform square grids. The mathematical analysis of the scheme was also performed for the Stokes equations in [17] for uniform rectangular meshes, and generalized to non uniform rectangular meshes and irregular source terms in [1]. Error estimates may be obtained by viewing the MAC scheme as a mixed finite element method of the vorticity formulation [12], or by a mixed method in primitive variables, with the pressures approximated by Q 1 finite elements [13]. Along the same lines, it is proven in [15] that a divergence conforming DG scheme based on the lowest order Raviart-Thomas space on rectangular meshes is algebraically equivalent to the MAC scheme. Error estimates for rectangular meshes were also obtained for the related covolume method, see [4] and references therein.
Mathematical studies of the MAC scheme for the non linear Navier-Stokes equations are scarcer. To our knowledge, the only convergence study is that of [18] for the steady-state Navier-Stokes equations and for uniform rectangular grids.
Extensions of the MAC scheme on general unstructured grids are not easy to derive, and most of the mathematical analyses that we described are restricted to rectangular cells. The covolume approach [19] may be seen as one way to generalize the MAC scheme on unstructured grids. MAC schemes on triangular meshes are proposed and tested in [25]. A variational MAC scheme based on acute triangles is proposed in [6], and the convergence analysis is completed.
The aim of this paper is to provide the complete mathematical analysis of an extension of the MAC scheme on possibly non conforming meshes, allowing local refinement. This extended MAC scheme is a slight modification of a MAC-like scheme which was presented in [3]. The modification concerns the discretization of the momentum equation, which was performed on dual Voronoï cells in [3], while it is performed by a linear finite element method on a Delaunay triangulation built from the centers of the set of edges (in 2D), or faces (in 3D) where each component of the velocity is defined. This modification was found necessary in the mathematical analysis of the scheme for the steady-state and time-dependent nonlinear Navier-Stokes equations, in order to prove the convergence of the scheme. It also allows to easily handle the case of the full tensor viscosity. The convergence of a collocated finite volume scheme was proven in [9] for the steady state and time-dependent Navier-Stokes equations, and we shall use some of the tools therein in the present analysis.
The outline of the paper is the following. In Sect. 2, we give the weak form of the steady-state and time-dependent Navier-Stokes equations. In Sect. 3 we first write a discrete variational formulation of the standard MAC scheme on the linear Stokes problem, and use this variational formulation to extend the MAC scheme to more complex geometries. The extension to the steady-state Navier-Stokes equation is presented in Sect. 4 where we also prove the convergence of the approximate solutions to a weak solution. The proof of convergence relies on a proper choice of the convection term, which takes the mass conservation into account. We then consider the case of the time-dependent problem in Sect. 5, and again prove the convergence of the method in this case, thanks to a discrete Aubin-Simon type result. In Sect. 6, the efficiency of the extended scheme is illustrated by numerical examples on a non-rectangular domain, using a local refinement along the boundary of the domain. The first example is a steady-state problem on a circular domain, for which the exact solution is known so that we may assess the numerical order of convergence. We then consider the inclined driven cavity problem, which implies the discretization of a non rectangular domain.
The comparison with the literature shows a very good accuracy. We finally consider the time-dependent Green-Taylor vortex problem, that we approximate on a circular domain.

Weak formulation of the Navier-Stokes equations
Let be an open bounded set of R d , where d denotes the space dimension. In the remainder of this paper, we assume the domain to be polygonal (d = 2) or polyhedral (d = 3), in the sense that ∂ is a finite union of subsets of hyperplanes of R d . Let β ∈ [0, +∞) denote the Reynolds number.
We first consider the steady state case with a time-independent forcing term f ∈ L 2 ( ) in the momentum equation. Then, a weak solution to the steady-state Navier-Stokes equations with homogeneous Dirichlet boundary conditions is a vector function u with components (u (i) ) i=1,...,d , such that We consider two cases for the stress tensor S(u, v), namely: -the usual simplified form for the incompressible Stokes or Navier-Stokes problem, which reads -the full stress tensor where μ > 0 and 3λ + 2μ ≥ 0 (these values could depend on the space variable through a coupled variable, such as the temperature) and The coefficient β is strictly positive in the general (nonlinear) case and is set to 0 to obtain the linear Stokes problem.
We then consider the time-dependent case, for which we consider a finite time T of study of the flow, a time-dependent forcing term f ∈ L 2 ( × (0, T )) d in the momentum equation, and an initial condition u ini ∈ L 2 ( ) d . Then, a weak solution to the time-dependent Navier-Stokes equations with homogeneous Dirichlet boundary conditions is a vector function u with components (u (i) ) i=1,...,d , such that Remark 2.1 In the case where ∂ is C 2 , it is proven in e.g. [24] that any weak solution u of (4) satisfies ∂ t u ∈ L 4/d (0, T ; E( ) ) in the following classical sense: One may prove the same regularity for a polygonal boundary by passing to the limit on the numerical scheme, but this is out of the scope of the present paper.

A variational MAC scheme
Our extension of the MAC scheme to non conforming meshes is based on a discrete variational formulation. Hence in this section, we begin by considering a variational formulation of the standard MAC scheme for the approximation of the Stokes problem, that is (1) with β = 0. We then extend this variational scheme to non conforming meshes.
3.  Let N and M be two positive integers, and let M be the set of pressure grid cells (with the notations given in Fig. 1: The above notations may easily be extended to the case d = 3. For the space dimension d equal to 2 or 3, we denote by E = ∪ d i=1 E (i) the set of the edges or faces of the mesh, where E (i) is the set of edges associated to the i-th component of the velocity. In order to define the normal velocity flux from one cell to a neighbouring one, we introduce, for any pair σ, σ ∈ E (k) , k = 1 or 2, the transmissivity τ where σ , and d(x σ , x σ ) denotes the distance between the points x σ and x σ . For instance, for a vertical edge (1) , one has: For any K ∈ M, we denote by E K the subset of E int containing all edges (or faces) of K which are internal, and by E (k) (k) . For an internal edge σ ∈ E int separating two cells K and L, we shall write σ ∈ K |L = L|K (at this stage, we could write σ = K |L but in the generalized MAC scheme, interfaces will be allowed to contain more than one edge or face). We denote by (e (k) ) k=1,...,d the canonical orthonormal basis of R d and, for σ ∈ K |L, with K , L ∈ M and by n K ,σ the unit normal vector to σ outward to K . We then write σ ∈ − − → K |L in the case where σ ∈ K |L ⊂ E (k) for some k = 1, . . . , d and n K ,σ · e (k) = 1. We finally represent by D = (M, E) the collection of all the space discretization data.
Let us consider here the simplest case of the Stokes equations, i.e. β = 0, with the stress tensor given by (2). The standard MAC scheme may then be written, only detailing the case of interior edges: where v K ,σ = v σ n K ,σ · e (k) . The discrete mass conservation Eq. (7a) can then be simply written as: We now introduce a variational formulation of the discrete momentum Eq. (7b) by recalling (see e.g. [7,Chapter 3] that if the points x σ are the nodes of a Delaunay triangulation, then τ (1) σ,σ = ∇ξ (1) σ (x) · ∇ξ (1) σ (x)dx, Fig. 2 Triangular mesh for the P 1 finite element approximation u (1) of u (1) where the functions (ξ (1) σ ) σ ∈E (1) are the P 1 finite element basis functions defined on the Delaunay triangulation, such as the one depicted in Fig. 2 (in which the nodes are the points x σ , for all σ ∈ E (1) ).
The functions (ξ (2) σ ) σ ∈E (2) are defined in a similar way, and we may then define We then define an inner product ·, · E on H E ( ) (this holds for both cases (2) and (3)) by which we expect to approximate the inner product S(u, v)(x)dx of the continuous problem. We then obtain, multiplying (7) by v σ and summing on k = 1, 2 and σ ∈ E (k) , A discrete variational formulation of the standard MAC scheme (7) is therefore: K ∈M |K | p K = 0 and (9) and (12) hold. (13) Note that, in the case of a 2D triangular mesh obtained by diagonal splitting of rectangles, the analogy of a five-point finite difference and the P 1 conforming finite element is classically due to the fact that the gradients of the P 1 basis functions associated to the two vertices of a diagonal of the mesh are orthogonal. This property cannot be generalized when splitting rectangular parallelepipedic grid blocks into tetrahedra in 3D.

The extended MAC scheme for non conforming meshes
Let us now turn to an extension of the above depicted standard MAC scheme to more general meshes (including local refinement as in Fig. 6; note that these local refinements may be used to follow the contours of a general non rectangular domain, as in Figs. 3 and 4. We consider 2D or 3D meshes of , which are such that all internal edges (2D) or faces (3D) (from now on, we only use the word "face", in 2D or 3D) have their normal vector parallel to one of the basis vector e (k) of the space R d , for some k = 1, . . . , d.
In other words, all internal edges must be aligned with one of the reference axes. Note that on the other hand, the external faces, that is the faces of the mesh lying on the boundary ∂ need not be aligned with the axes: they are only assumed to be planar. In fact a tilted cavity such as the one depicted in Fig. 3 is easily meshed with such a grid. Curved boundaries may also be meshed with such grids, by using local refinement close to the boundaries, such as in Fig. 5.
We denote by M the set of pressure cells. Examples of resulting pressure grids are depicted in the left part of Figs. 3 and 5. For all K ∈ M, we again denote by E K the set of all internal faces of K (therefore, the faces of K which are on ∂ are not elements of E K ), and we define the set E as the union over K ∈ M of all the sets E K . It is assumed that a given σ ∈ E int is entirely included in an interface between two cells, say K and L; we shall write σ ∈ K |L. Note that K |L is allowed to contain several faces of the mesh; this may for instance happen in adaptive mesh refinement (and de-refinement) procedures. We then introduce the set E (k) as the subset of E which contains all the internal faces whose normal is parallel to the basis vector e (k) . For any σ ∈ E, x σ denotes the center of gravity of σ .
In order to get a discrete variational MAC-like scheme, we consider, for any k = 1, . . . , d, the set of internal points V (k) int = (x σ ) σ ∈E (k) and a given family of external points V (k) ext containing at least all the vertices of . The only additional difficulty on the non conforming grid is the discretization of the diffusion term, which we simply discretize by a linear finite element approximation. To this purpose, we introduce a Delaunay triangulation T (k) of whose vertices are V ext . Such a triangulation of is a set of simplices: triangles in 2D, tetrahedra in 3D. Each simplex has d + 1 vertices which belong to The triangulation is assumed to satisfy the Delaunay property, which means that the interior of the circumcircle (in 2D) or of the circumsphere (in 3D) of any simplex T ∈ T (k) does not contain any element of V ext . Examples of Delaunay triangulations constructed from the edge mid-points are illustrated in Figs. 3 and 5 (middle and right parts).
We then denote, for any σ ∈ E (k) , the function ξ (k) σ , which is continuous, piecewise P 1 on any T ∈ T (k) , and whose value is 1 at the point x σ and 0 at points x σ for any , defining the Voronoï cells as follows:

See Figs. 4 and 6 for examples of the superposition of the Delaunay triangulations and
the Voronoï mesh construction. We finally denote by D the collection of all the space discretization data.
Remark 3.1 Note that in the case of a uniform rectangular mesh, the Voronoï cells thus defined are equal to the velocity cells defined in the previous section. However, this is no longer true if a non uniform mesh is used, even in the conforming case; indeed, in this latter case, the Voronoï cells V (k) σ are again rectangles, but they are not equal to the (rectangular) velocity cells of the classical MAC scheme. In the case of hanging nodes, they are no longer rectangular, as can be seen in Fig. 6, where we depict the Delaunay grid and the Voronoï cells for the horizontal and vertical velocities.
We may again define H M ( ) as the set of piecewise functions constant on the pressure cells K ∈ M, the set H (k) E ( ) of piecewise constant functions on the dual grid cells V σ , for σ ∈ E (k) ; this discrete set is the space of functions meant to approximate the k-th component of the velocity. We then denote by we denote by v = ( v (k) ) k=1,...,d , and we define the norm The extended MAC scheme for the Stokes equations (β = 0) is again (8)-(13), applying Definition (11) for ·, · E .

Discretization of the nonlinear convection term
In order to write this generalized scheme for the Navier-Stokes equations, we only need to give a discretization of the nonlinear term (u(x) · ∇)u(x) · v(x)dx. To this purpose, we introduce a discrete trilinear form b E which aims at discretizing the trilinear form b defined over ( We begin by defining some interpolation operators between H E ( ) and (H M ( )) d .
and M v ∈ H M ( ) = (H M ( )) d as the piecewise constant function equal to ( K v) on cell K . We then define E v ∈ H E ( ) as the following piecewise constant function on the Voronoï cells:

Definition 4.2 (Discrete gradient and convection term).
For where, for any v ∈ H M ( ), It is well-known that the continuous trilinear form satisfies Similarly, we have the following result for the discrete trilinear form. Lemma 4.1 (Properties of the trilinear form). Let u, v, w ∈ H E ( ) and let b E be defined by (18).
where the discrete divergence operator div M is defined by (8); Reordering the summation over the edges of each cell K yields (19).
The extended MAC scheme for the Navier-Stokes equation then reads: (21c) Remark 4.1 (Link with the classical MAC scheme). In the case of the Stokes problem on a conforming rectangular grid, the scheme presented here is equivalent to the classical MAC scheme. In the case of the nonlinear Navier-Stokes equations, the scheme is not quite the same, since the centred MAC scheme has a 3 point stencil in each direction whereas the scheme presented here has a 5 point stencil. This is due to the fact that we chose to take E w rather than w in the expression of the trilinear form (18), in order to obtain the property (20), which is quite useful to obtain the estimates on u.

Mathematical analysis
Let D be a possibly non conforming mesh, such as introduced in Sect. 3.2. We define h D as the maximum value of the diameters of all T ∈ T (k) and all V The regularity θ D of D is defined as the minimum value of: 3. all ratios |V Let us then introduce the following interpolation operators.

Definition 4.3 (Interpolation operators). For
For ϕ ∈ H 1 0 ( ), let P M ϕ ∈ H M ( ) be defined by its piecewise constant values on the cells K ∈ M: Following the ideas of [22] applied to the piecewise P 1 interpolation P E v of v, one may check that the interpolation operator P E satisfies the following property.

Lemma 4.3 (Properties of the discrete divergence).
Let v ∈ (H 1 0 ( )) d , let P E be the interpolation operator defined by (22).and let p ∈ H M . Then Moreover, for any u ∈ H E , and for any q ∈ H M ( ), one has: Proof The first result is an obvious consequence of the relation 8 defining the discrete divergence operator and of the fact that p is piecewise constant on the cells K ∈ M. The second result follows from the fact that The next lemma, which is a technical result allowing to bound, under adequate geometrical conditions, the difference between various interpolates by the gradient of a P1 reconstruction, is used several times in the estimates and convergence proofs. It is proved in the appendix.

Lemma 4.4 (Comparison between interpolates). Let
⊂ R d , with d ∈ N * be an open polygonal connected open set, such that there exists a H 1 conforming simplicial mesh T over . We denote by T the elements of M, by V the set of all vertices, and by V T the set of all vertices of T and by T T ⊂ T the set of the simplices sharing a face with T . For s ∈ V, we denote by ξ s the P 1 finite element basis function associated to the node x s . We assume that are given some non-negative functions ψ s , for all s ∈ V, such that For (u s ) s∈V ∈ R V , let us denote u = s∈V u s ξ s and u = s∈V u s ψ s , and, for any T ∈ T and x ∈ T , let us denote δ(x) = diam(T ). Then there exists C 3 > 0 only depending on C 2 , such that so that in particular, We also need several times the following lemma, also proved in the appendix, allowing the comparison between discrete norms.
The following lemma gives some estimates on the trilinear form which are used to obtain some estimates on the solutions of the schemes, for both the steady-state treated in this section and the time-dependent case in Sect. 5.  (19). Then there exist C 4 and C 5 , only depending on any θ ≤ θ D and on , such that Moreover, there exists C 6 , only depending on and θ , such that Proof From the definition (18) of the trilinear form b E denoting byṽ ( j) (resp.w ( j) ) the j-th component of M v (resp. E w), we get thanks to the Cauchy-Schwarz inequality, we get that: Let us then remark that the components of K w+ L w 2 are barycentric combinations of neighbouring terms w σ ; hence there exists C 7 , only depending on and θ , such that In order to bound the term ∇ E ( M v) (L 2 ( )) d×d , we first define the functions We finally define δ (k) σ . We then get that In order to apply Lemma 4.4, we first remark that, thanks to the regularity hypotheses of the mesh, there exists some C 8 > 0, only depending on θ , such that δ (k) only depend on the discrete unknowns v σ , we may introduce the functions ψ k, −,σ such that which are the piecewise constant functions defined by: The functions ψ k, −,σ satisfy the hypothesis (27) of Lemma 4.4; furthermore, thanks to the regularity of the mesh, the mesh dependent bound C 2 defined by (28) remains bounded by a function of θ . We may thus apply Lemma 4.4, which yields (29). This implies Similarly, we have the same inequality with (k) − v ( ) , which shows the existence of C 9 > 0, only depending on θ , such that thus concluding the left inequality of (32). Now, thanks to the equivalence of the norms proved in Lemma 4.5 in the appendix (inequality (31)), there exists C 10 ∈ R + , only depending on any θ ≤ θ D , such that Applying the standard Sobolev inequality u (k) where C sob only depends on and d, we conclude that there exists C 11 , only depending on and on any θ ≤ θ D , such that Therefore, we conclude the right inequality of (32). Let us then prove (33); since div M u = 0, we have u, v, u), which proves, using (32), that Using the Cauchy-Schwarz inequality, we have: We again apply Lemma 4.5: there exists C 12 ∈ R + , only depending on θ and , such that Applying the standard Sobolev inequality, we get that u L 6 ( ) ≤ C sob ∇ u L 2 ( ) d×d , where C sob ∈ R + depends only on and d; therefore, by the definition (15) of the norm · E , we get: which concludes the proof of Lemma 4.6.  (21). Then there exists C 13 , only depending on (and on λ and μ in the case (3), this dependency is no longer mentioned in this paper), such that u satisfies the following estimates: Moreover, there exists C 14 only depending on β, θ and such that As a direct consequence, there exists one and only one solution (if β = 0) and at least one solution (if β > 0) to the scheme (21).
Proof We let v = u in (21). Let us first remark that we have p(x)div M u(x)dx = 0 and that by the definition (19) Hence we obtain where we denote by C 15 > 0 the constant involved in the Korn inequality, which holds for both cases (2) and (3). We then have, thanks to the Cauchy-Schwarz inequality, Let u be defined by (10). Applying Lemma 4.5, we get that there exists C 16 , only depending on θ , such that Thanks to the Poincaré inequality, we may write Therefore we get (34). Let us now turn to the proof of (35). We use the following property, first due to Nečas [16], see also [2]: there exists v ∈ H 1 0 ( ) d such that p(x) = divv(x) for a.e. x ∈ and there exists C 17 , only depending on , such that Let P E be the interpolation operator defined by (22). By Lemma 4.3, we have: We then get, using P E v ∈ H E ( ) as test function in Scheme (21), Thanks to the Cauchy-Schwarz inequality, we get that there exists C 18 (only depending on λ and μ in the case (3), this dependency is no longer recalled), such that Using inequality (24), we get that there exists C 19 , only depending on θ and , such that We may also write, thanks to the Poincaré inequality and to (37), Therefore, using (34) and (32), we get the existence of C 20 such that which concludes the proof of (35).
In the case where β = 0, the inequalities (34) and (35) suffice to prove that the square linear system issued from Scheme (21) (replacing one equation of (9) by p(x)dx = 0) is invertible. If β > 0, using the topological degree argument as done in [9], we conclude to the existence of at least one solution to Scheme (21).
The next lemma states the weak convergence of the discrete gradient and the regularity of a possible limit; it is a straightforward adaptation of [8,Lemma 5.7]. In particular, we obtain the H 1 0 regularity of the limit, thanks to the fact that the Dirichlet boundary conditions are taken into in the definition of the gradient. In [7, Chapters 2 and 3], the fact that u ∈ H 1 0 is obtained by extending the discrete functions by 0 and passing to the limit on the estimates on the translations in the L 2 norm: these latter estimates are obtained under an orthogonality condition on the mesh, thanks to the fact that the Dirichlet boundary conditions are taken into account in the discrete inner product and norm defined on the space of discrete functions.
In [8,Lemma 5.7], the fact that u ∈ H 1 0 is obtained again by extending the discrete functions by 0, but now with a passage to the limit on an approximate gradient, for the same type of meshes than those considered here. The proof of Lemma 4.8 is contained in step 3 of the proof of Lemma 5.7 of [8], we give it here for the sake of completeness.
the family (v D ) D∈F converges weakly in L 2 ( ) to v ∈ L 2 ( ) as h D → 0. (40) Then v ∈ H 1 0 ( ), and the family (∇ E v D ) D∈F converges to ∇v weakly in L 2 ( ) as h D tends to 0.
where P E is the interpolation operator defined by (22), and |R E | ≤ C ψ h D ∇ E v D (L 2 ( )) 2 → 0 as h D → 0, thanks to Assumption (39). Now, by (26) and (25) of Lemma 4.3, Therefore, by Assumption (40), which shows that ∇ E v D tends to ∇v in the distribution sense. Since the sequence (∇ E ) D∈F is bounded in L 2 ( ), we have that ∇v ∈ L 2 ( ), and prolonging v D by 0 outside of , we get that  (21). Then there exists a weak solution (u, p) of (1) such that, up to a subsequence, u E converges in L 2 ( ) d to u and p M converges in L 2 ( ) d to p as h D → 0. Moreover, if β = 0, the whole family converges to the unique weak solution of (1) as h D → 0.
Proof Let (D ( ) ) ∈N be a sequence of elements of F, such that h D ( ) tends to 0 as → ∞ and such that there exists θ > 0 with θ D ( ) ≥ θ , for all m ∈ N. Let (u E ( ) , p M ( ) ) denote a solution of (21) for the discretization D = (M , E , δt ). From estimate (34), using Rellich's theorem, we deduce that there exists u ∈ H 1 0 ( ) d (with u = (u (k) ) k=1,...,d ) and a subsequence of (D ( ) ) m∈N , again denoted by (D ( ) ) ∈N , such that u (k) M ( ) ∈ H 1 0 ( ) weakly converges in H 1 0 ( ) d (therefore strongly in L 2 ( ) d ) to u (k) as → ∞ for k = 1, . . . , d. Using Lemma 4.4 and the estimate (34), we get that u E ( ) also converges in L 2 ( ) d to u as → ∞. Thanks to the estimate (35) on the pressure, we may then consider a subsequence of (D ( ) ) ∈N , again denoted by (D ( ) ) ∈N , such that p M ( ) weakly converges in L 2 ( ) to some function p ∈ L 2 ( ) as → ∞, with p(x)dx = 0. We now have to prove that (u, p) is a weak solution of (1), which we do by passing to the limit on the weak form of the scheme.
Let us first show that u ∈ E( ), i.e. divu(x) = 0 for a.e. x ∈ . Let ϕ ∈ C ∞ c ( ); multipying (21b) by P M ϕ and integrating over yields, thanks to Lemma 4.3: Thanks to the regularity of the mesh, there exists C 21 depending only on ϕ and θ such that ∇ E P M ϕ (L 2 ( )) d ≤ C 21 . Therefore, by Lemma 4.8, ∇ E P M ϕ tends to ∇ϕ weakly in (L 2 ( )) d as h D tends to 0. Hence, passing to the limit in the above equation, we get that ∇ϕ(x) · u(x)dx = 0 for any ϕ ∈ C ∞ c ( ), which shows that divu(x) = 0 for a.e. x ∈ .
Let ϕ ∈ C ∞ c ( ) d . Let us now show (1) with v = ϕ. We take (again omitting some indices ( ), thus denoting D = D ( ) ) v E = P E ϕ ∈ H E ( ), as defined by (22), as test function in Scheme (21). We get that Using Lemma 4.3 and the fact that p D ( ) weakly converges in L 2 ( ) to p ∈ L 2 ( ), we obtain that Thanks to the weak and strong convergence properties of the different sequences, we get that We also have, thanks to the definition (22) of P E ϕ, Let us now turn to the study of the limit of b E ( ) (u D ( ) , u D ( ) , P E ( ) ϕ). By the definition (19) of b E , again dropping some indices ( ), we have: Thanks to the regularity of the mesh, the family ∇ E ( D u E ) is bounded in (L 2 ( )) d . Moreover, using Lemma 4.4 and the estimate (34), we get that there exists C 22 , only depending on θ , such that which completes the proof that (1) holds for all v ∈ H 1 0 ( ) d by density.

Definition of the scheme
Let D be the collection of all the space discretization data as given in Sect. 4.1. The discrete approximation of the time-dependent Navier-Stokes equations (4) is then given by the following relations. The initial condition is approximated by and, for a given time step δt > 0, let us denote by t (n) = nδt for n ≥ 0. The scheme is defined by and, for a given α ∈ [ 1 2 , 1], Note that we use the notation u (n+ 1 2 ) for any value α ∈ [ 1 2 , 1], and not only α = 1 2 . We then denote and we define the discrete time derivative Remark 5.1 (Time discretization) We consider a constant time step only for the sake of clarity of notations. The mathematical analysis is still valid with a variable time step. Note that the above scheme corresponds to a Crank-Nicolson-like scheme for α = 1 2 , and to an implicit scheme for α = 1.

Mathematical analysis
Our goal is to use, in our proof of convergence of the scheme, an adaptation of a previous result [11]. The differences between the theorem given below and [11,Theorem 3.4] are the following; 1. Theorem 5.1 is a compactness result in L 1 (0, T ; B) where B is some Banach space, while [11,Theorem 3.4] is a compactness result in L p (0, T ; B) for p ∈ [1, +∞). We chose to restrict to L 1 for two reasons: • The proof, which is given in the appendix, is much simpler in the L 1 case, especially in the case of variable time steps. • In our framework,we can easily get the result in L 2 by using compactness in L 1 and the discrete estimates (47) and (48), using the "L p − L q compactness" property and Sobolev inequalities, see [10, Lemma 3.1]. 2. The assumption (h2) that we require here is weaker than the assumption (h2) of [11,Theorem 3.4], which does not require that ( w n X n ) ∈N be bounded. A quick look at the proofs of [11,Lemma 3.1 and Theorem 3.4] shows that in fact they only require our hypothesis (h2) (which is easier to verify in our framework). 3. We give here a version with variable time steps. 4. The definition (46) of u (·, t) corresponds to an α−scheme whereas that of [11,Theorem 3.4] corresponds to an implicit scheme, i.e. α is equal to 1.
Let δ u be the "discrete time derivative", defined by: Let · X and · Y be two norms on B . We denote by X the space B endowed with the norm · X and by Y the space B endowed with the norm · Y . We assume that (h1) For any sequence (w ) ∈N such that w ∈ B and ( w X ) ∈N is bounded, then, up to a subsequence, there exists w ∈ B such that w → w in B as → +∞. Then there exists u ∈ L 1 (0, T ; B) such that, up to a subsequence, u → u in L 1 (0, T ; B) as → +∞.
Remark 5.2 (The original Aubin-Simon theorem) Note that in fact, our weaker assumption (h2) has a continuous equivalent, and indeed, we may weaken one of the hypotheses of the classical Aubin-Simon theorem [23]; indeed, the hypothesis "B continuously embedded in Y " , which is equivalent to the following hypothesis: If w n → w in B and w n Y → 0, then w = 0, may be replaced by the weaker assumption "if (w n ) n∈N is bounded in X , converge in B to w and w n Y → 0, then w = 0".  23 , only depending on f , u ini , T and , such that the following discrete L 2 (H 1 ) and L ∞ (L 2 ) estimates hold: and u (n) Moreover we have the following L 1 (L 2 ) estimate on the pressure: where C 24 only depends on D, β, δt, f , u ini and . Therefore there exists at least one solution to (43)-(44) if β > 0 and exactly one solution if β = 0.
Proof Let us first remark that (48), for n = 0, is due to the Cauchy-Schwarz inequality applied to (41). We let u = u (n+ 1 2 ) in (44). The nonlinear term again vanishes, which leads, using the Young, Poincaré and Korn inequalities, to: Summing for n = 0, . . . , N − 1 lead to both (47) and (48). Following the same ideas as in the proof of Lemma 4.7, we get the existence of C 25 , only depending on θ , β and such that Summing the above inequality on n = 0, . . . , N T − 1, using the Cauchy Schwarz inequality for the terms δt N T −1 E and inequality (48) for the term (49). The existence of at least one solution to the scheme is then again deduced from the use of the topological degree results, as in [9], as well as the uniqueness of the solution if β = 0.
Since our aim is to apply Theorem 5.1, we define a second norm on H E ( ) in the next definition.
Let us denote by Y E the space H E ( ) equipped with the norm · ,E . We define a continuous embedding of where P E v is defined by (22). Note that, thanks to (24) we get that where C 1 only depends on θ and .
In view of applying Theorem 5.1, let us study the dual norm of the discrete time derivative.
Since u (n+ 1 2 ) satisfies (43) we get thanks to the estimate (33) of Lemma 4.6 and to the estimate (48) that: where C 27 , only depends on f , u ini , T , and θ . We then get that there exists C 28 , only depending on f , u ini , T , β, and θ , such that Remarking that u (n+ 1 2 ) we get that there exists C 29 , only depending on f , u ini , T , β, and θ , such that δ D,δt u(t) 4 3 * ,E ≤ C 29 Integrating for t ∈ (0, T ) and using the discrete L 2 (H 1 ) estimate on the velocity (47), we get that there exists C 26 , only depending on f , u ini , T , , β and θ , such that (53) holds, which concludes the proof.  1]. Then there exists a weak solution u of (4) such that, up to a subsequence, u E,δt converges in L 2 ( × (0, T )) d to u as h D → 0 and δt → 0, which shows that Moreover, if β = 0, the whole family converges to the unique weak solution of (4) as h D → 0 and δt → 0.
Proof Let D = (M , E , δt ) ∈N be a sequence of elements of F, such that h D and δt tend to 0 as → ∞. From thereon, we denote u M ,δt by u for short.
Step 1 Proof that hypotheses (h1-h4) of Theorem 5.1 hold, and consequences. In our setting, the space B of Theorem 5.1 is L 2 ( ) d . We take B = {w ∈ H E ; div D w = 0}. The norm · X is the norm · E , and the norm · Y is defined in Definition 5.1.
Let (w ) ∈N be a sequence of functions of B such that w E ≤ C for some C ∈ R + . Then by definition of the norm w E and thanks to Rellich's theorem, the sequence ( w ) E converges in L 2 ( ) d to some w ∈ L 2 ( ) d ; therefore, by the inequality (30) of 4.4 given in the Appendix, the sequence (w ) E also converges to w in L 2 ( ) d . Thus, assumption (h1) of Theorem 5.1 is satisfied.
Let us then show that assumption (h2) is also satisfied. Let (w ) ∈N be a sequence of functions of L 2 ( ) d such that w ∈ B and w E ≤ C for some C ∈ R + , and such that there exists w ∈ B with w → w in B and w ,E → 0 as → ∞. By definition of the norm · ,E , we have which shows that w = 0.
Step 2 Proof that u is a weak solution of (4). As in [9], we easily get that u ∈ L 2 (0, T ; H 1 0 ( ) d ) and we prove the convergence in L 2 ( × (0, T )) d to u of u . Lemma 4.4 implies the same convergence property for u . We get that divu(x, t) = 0 in the same manner as in the steady case. The proof that u satisfies (4) follows the same steps as the Navier-Stokes steady case, using a divergence free test function. The convergence of the additional time term to the corresponding continuous limit is then classical.
Step 3 Proof that (54) holds. Thanks to (53), extracting a subsequence such that δ D ,δt u weakly converges in the Sobolev space L 4/3 (0, T ; E( ) ), we may pass to the limit in (51). We thus get The preceding relation allows to conclude that and therefore that (54) holds.

An analytical steady problem
We consider a problem where the continuous solution of the Navier-Stokes equations (1) in the case (2) with β = 1 is given by: in a circle with centre (0.5, 0.5) and radius 0.45. We consider four meshes for the mass conservation M j , j = 0, . . . , 3, defined in the following way: Left part of Fig. (7) shows the errors log 10 (e 2 (u 1 )) and log 10 (e 2 ( p)) with respect to log 10 (1/ card(M j )) for j = 0, . . . , 3. On right part of Fig. 7 are plotted the stream lines for the finest mesh. The velocity components and the pressure are shown in Fig. (8) for two meshes.
Although the velocity fields are accurately computed on the coarsest mesh, the pressure fields show oscillations where neighbouring control volumes have contrasted sizes. However, these oscillations are decreasing while refining the mesh.

An inclined driven cavity
We consider a 30 • inclined driven cavity, which corresponds to non homogeneous Dirichlet boundary conditions for the velocity on the horizontal upper boundary of a parallelogram (see Fig. 9) in the case (2). We consider the case where the Reynolds number is = 1,000, and we discretize the domain using refinement strategies. We thus show the flexibility of this extended MAC scheme in this case. We obtain the result shown in Fig. 10. The results are compared to the literature [5] in Table 1. The results show the expected precision. The solution is given, in the case (2), by We plot in Fig. 11 the L 2 errors evaluated at time t = 0.1 as a function of the space step (−0.5 × log 10 (#M)). For each computation, the time step is set such that the temporal error is negligible with respect to the spatial one. The sequence of time and space steps are given in Table 2: the ratio of two successive time steps is 4 whereas it is only 2 for the space step. We observe that the approximate pressure does not seem to converge in the case where α = 1/2, which does not occur in the case α = 1. If we plot the values 1 2 ( p (n− 1 2 ) + p (n+ 1 2 ) ) instead of p (n+ 1 2 ) , we observe that this postprocessed pressure numerically converges. All the converging curves show an order 2. In order to check the behaviour of the pressure in the case α = 1/2, we show, in Figs. 12 and 13, the pressures fields, for two meshes, the left and right figures showing two consecutive times steps, and the middle one showing the average value between these ones. This shows that, although the scheme obtained in the case α = 1/2 may lead to oscillations in pressures, it is simple to get back converging pressures.
In order to check the convergence order with respect to the time step, we have chosen the sequence described by Table 3. The errors are plotted in Fig. 14 (for the case α = 1/2, we have compared the exact pressure with the average pressure    1 2 ) ) instead of p (n+ 1 2 ) ). This shows an order 2 for α = 1/2, and 1 for α = 1.

Conclusion
The extension of the MAC scheme presented here seems to show interesting properties for practical computations. The grids are easy to construct, and complex domains may easily be discretized with the possibility of local refinement. This could lead to the use of this scheme in industrial and large scale problems. Work is in progress to design a discrete nonlinear convection term which would resume to the classical MAC scheme in the case of a uniform rectangular grid, while retaining the same properties in order to obtain the mathematical convergence of the discrete solutions.
Proof of Lemma 4.5 Proof The left inequality of (31) is an immediate consequence of the convexity of the function x → |x| p , which provides s∈V u s ξ s (x) p ≤ s∈V |u s | p ξ s (x).
Let us turn to the proof of the right inequality of (31). Let T ∈ T and (a 1 , . . . , a d+1 ) ∈ R d+1 be given, and let, for any x ∈ T , ξ 1 (x), . . . , ξ d+1 (x) be the d + 1 barycentric coordinates of x with respect to the vertices (recall that the P 1 basis functions are defined by such coordinates in each simplex). We then denote, for a given ε ∈ (0, 1/2), which will be chosen later, by Since, on T i , we have, for j = i, ξ j (x) ≤ ε < 1 − ε, we get that all the T i , i = 1, . . . , d + 1, are disjoint. We then have We may write, for x ∈ T d+1 (ε), and therefore, thanks to the Young inequality Let us remark that the measure of T i (ε) is equal to ε d |T |. We then get We then choose ε such that We then have and therefore Denoting C(d, p) > 0 the quantity defined by we have then proved that which implies (31).