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

Eric Chénier, Robert Eymard, Thierry Gallouët, Raphaele Herbin

► To cite this version:

Eric Chénier, Robert Eymard, Thierry Gallouët, Raphaele Herbin. An extension of the MAC scheme to locally refined meshes: convergence analysis for the full tensor time-dependent Navier-Stokes equations. Calcolo, 2014, 52 (1), pp.69-107. 10.1007/s10092-014-0108-x. http://dx.doi.org/10.1007/s10092-014-0108-x

HAL Id: hal-00751556

https://hal.science/hal-00751556

Submitted on 14 Nov 2012

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An extension of the MAC scheme to locally refined meshes : convergence analysis for the full tensor time-dependent Navier-Stokes equations

Eric Chénier*, Robert Eymard†, Thierry Gallouët‡ and Raphaèle Herbin§

November 13, 2012

## Abstract

A variational formulation of the standard MAC 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.

**Keywords:** MAC scheme, incompressible steady and time-dependent Navier-Stokes equations, non conforming grids.

**MSC2010:** 65N08,76D05

## 1 Introduction

The Marker-And-Cell (MAC) scheme, introduced in [13] is one of the most popular methods [19, 24] for the approximation of the Navier-Stokes equations in the engineering framework, because of its simplicity and of its remarkable mathematical properties. The discrete unknowns are the components of the velocity and the pressure on staggered grids; the mass conservation equation and momentum conservation equations are discretized in such a way that the kinetic energy remains controlled. The first error analysis seems to be that of [20] 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 [16] 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 [11], or by a mixed method in primitive variables, with the pressures approximated by $Q^1$ finite elements [12]. Along the same lines, it is proven in [14] 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 [17] 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 [18] may be seen as one way to generalize the MAC scheme on unstructured grids. MAC schemes on triangular

<footnote>* Université Paris-Est, eric.chenier@univ-mlv.fr</footnote>

<footnote>† Université Paris-Est, robert.eymard@univ-mlv.fr</footnote>

<footnote>‡ Aix-Marseille Université, thierry.gallouet@univ-amu.fr</footnote>

<footnote>§ Aix-Marseille Université, raphaele.herbin@univ-amu.fr</footnote>
meshes are proposed and tested in [23]. 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 [8] 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 Section 2, we give the weak form of the steady-state and time-dependent Navier-Stokes equations. In Section 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 Section 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 Section 5, and again prove the convergence of the method in this case, thanks to a discrete Aubin-Simon type result. In Section 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.

## 2 Weak formulation of the Navier-Stokes equations

Let $\Omega$ be an open bounded set of $\mathbb{R}^d$, where $d$ denotes the space dimension. In the remainder of this paper, we assume the domain $\Omega$ to be polygonal ($d=2$) or polyhedral ($d=3$), in the sense that $\partial\Omega$ is a finite union of subsets of hyperplanes of $\mathbb{R}^d$. Let $\beta \in [0, +\infty)$ denote the Reynolds number.

We first consider the steady state case with a time-independent forcing term $f \in L^2(\Omega)$ 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

$$ egin{cases} 
u \in E(\Omega), p \in L^2(\Omega) \text{ with } \displaystyle\int_{\Omega} p(x)dx = 0, \\ 
\displaystyle\int_{\Omega} S(u, v)(x)dx + \beta \int_{\Omega} (u(x) \cdot \nabla)u(x) \cdot v(x)dx \\ \qquad - \displaystyle\int_{\Omega} p(x)\mathrm{div}v(x)dx = \displaystyle\int_{\Omega} f(x) \cdot v(x)dx, \forall v \in H_0^1(\Omega)^d, 
\end{cases} 
\quad (1) $$

where $E(\Omega) = \{v = (v^{(i)})_{i=1,...,d} \in H_0^1(\Omega)^d, \mathrm{div}v = 0 \text{ a.e. in } \Omega\}$. 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

$$ S(u, v)(x) = \nabla u(x) : \nabla v(x) = \sum_{i=1}^{d} \nabla u^{(i)}(x) \cdot \nabla v^{(i)}(x), \text{ for a.e. } x \in \Omega, \forall u, v \in H_0^1(\Omega)^d, \quad (2) $$

- the full stress tensor

$$ S(u, v)(x) = \lambda \mathrm{div} u(x) \mathrm{div} v(x) + 2\mu \epsilon(u)(x) : \epsilon(v)(x), \text{ for a.e. } x \in \Omega, \forall u, v \in H_0^1(\Omega)^d, \quad (3) $$

where $\mu > 0$ and $3\lambda + 2\mu \ge 0$ (these values could depend on the space variable through a coupled variable, such as the temperature) and $\epsilon(u)_{i,j}(x) = \frac{1}{2}(\partial_i u^{(j)}(x) + \partial_j u^{(i)}(x))$, for $i, j = 1, ..., d$.
The coefficient $\beta$ 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 \in L^2(\Omega \times (0, T))^d$ in the momentum equation, and an initial condition $u_{ini} \in L^2(\Omega)^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

$$ 
\left\{ 
\begin{array}{l}
u \in L^2(0, T; E(\Omega)) \cap L^\infty(0, T; L^2(\Omega)^d), \\
- \int_0^T \int_\Omega u(x, t) \cdot \partial_t \varphi(x, t) \,dx\,dt - \int_0^T \int_\Omega u_{ini}(x) \cdot \varphi(x, 0) \,dx\,dt \\
+ \int_0^T \int_\Omega S(u, \varphi)(x, t) \,dx\,dt + \beta \int_0^T \int_\Omega (u(x, t) \cdot \nabla)u(x, t) \cdot \varphi(x, t) \,dx\,dt \\
= \int_0^T \int_\Omega f(x, t) \cdot \varphi(x, t) \,dx\,dt, \quad \forall \varphi \in L^2(0, T; E(\Omega)) \cap C_c^\infty(\Omega \times (-\infty, T))^d.
\end{array}
\right. 
\quad (4) 
$$

**Remark 2.1** It may be proved that any weak solution $u$ of (4) satisfies $\partial_t u \in L^{4/d}(0, T; E(\Omega'))$ in the following classical sense

$$ \int_0^T \varphi(t) \langle \partial_t u(t), v \rangle_{E(\Omega)', E(\Omega)} dt = - \int_0^T \varphi'(t) \int_{\Omega} u(x, t) v(x) dx dt, \\ \forall v \in E(\Omega), \forall \varphi \in C_c^1((0, T)). $$

# 3 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 $\beta = 0$. We then extend this variational scheme to non conforming meshes.

## 3.1 The standard MAC scheme for conforming meshes

For the sake of simplicity, let us first consider a two dimensional square domain $\Omega = [x, \bar{x}[\times]y, \bar{y}[$. The domain $\Omega$ is discretized in the MAC way, that is with a staggered arrangement of the rectangular discretization cells for the pressure and each of the velocity components, as depicted in Figure 1. Let $N$ and $M$ be two

<figure><img src="image_2.png" /><figcaption>Figure 1: Notations for the standard MAC scheme</figcaption></figure>

positive integers, and let $\mathcal{M}$ be the set of pressure grid cells (with the notations given in Figure 1:

$$ \mathcal{M} = \{ ]x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}[ \times ]y_{j-\frac{1}{2}}, y_{j+\frac{1}{2}}[ , 1 \le i \le N, 1 \le j \le M \}. $$
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 $\mathcal{E} = \cup_{i=1}^d \mathcal{E}^{(i)}$ the set of the edges or faces of the mesh, where $\mathcal{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 $\sigma, \sigma' \in \mathcal{E}^{(k)}$, $k = 1$ or 2, the transmissivity $\tau_{\sigma, \sigma'}^{(k)}$ between cell $V_{\sigma}^{(k)}$ and cell $V_{\sigma'}^{(k)}$:

$$ \tau_{\sigma, \sigma'}^{(k)} = \frac{|\partial V_{\sigma}^{(k)} \cap \partial V_{\sigma'}^{(k)}|}{d(\mathbf{x}_{\sigma}, \mathbf{x}_{\sigma'})}, \qquad (5) $$

where $|\partial V_{\sigma}^{(k)} \cap \partial V_{\sigma'}^{(k)}|$ denotes the length of the line segment which is the intersection of $\partial V_{\sigma}^{(k)}$ and $\partial V_{\sigma'}^{(k)}$, and $d(\mathbf{x}_{\sigma}, \mathbf{x}_{\sigma'})$ denotes the distance between the points $\mathbf{x}_{\sigma}$ and $\mathbf{x}_{\sigma'}$. For instance, for a vertical edge $\sigma = \{x_{i+\frac{1}{2}}\} \times y_{j-\frac{1}{2}}$, $y_{j+\frac{1}{2}} \in \mathcal{E}^{(1)}$, one has:

$$ \tau_{\sigma, \sigma'}^{(1)} = \begin{cases} \frac{y_{j+\frac{1}{2}} - y_{j-\frac{1}{2}}}{x_{i+\frac{3}{2}} - x_{i+\frac{1}{2}}} & \text{if } \sigma' = \{x_{i+\frac{3}{2}}\} \times [y_{j-\frac{1}{2}}, y_{j+\frac{1}{2}}], \\ \frac{x_{i+1} - x_i}{y_{j+1} - y_j} & \text{if } \sigma' = \{x_{i+\frac{1}{2}}\} \times [y_{j+\frac{1}{2}}, y_{j+\frac{3}{2}}]. \end{cases} \qquad (6) $$

For any $K \in \mathcal{M}$, we denote by $\mathcal{E}_K$ the subset of $\mathcal{E}_{\text{int}}$ containing all edges (or faces) of $K$ which are internal, and by $\mathcal{E}_K^{(k)} = \mathcal{E}_K \cap \mathcal{E}^{(k)}$. For an internal edge $\sigma \in \mathcal{E}_{\text{int}}$ separating two cells $K$ and $L$, we shall write $\sigma \in K|L = L|K$ (at this stage, we could write $\sigma = 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 $\mathbb{R}^d$ and, for $\sigma \in K|L$, with $K, L \in \mathcal{M}$ and by $\mathbf{n}_{K,\sigma}$ the unit normal vector to $\sigma$ outward to $K$. We then write $\sigma \in \overrightarrow{K} | L$ in the case where $\sigma \in K|L \subset \mathcal{E}^{(k)}$ for some $k = 1, ..., d$ and $\mathbf{n}_{K,\sigma} \cdot \mathbf{e}^{(k)} = 1$. We finally represent by $\mathcal{D} = (\mathcal{M}, \mathcal{E})$ the collection of all the space discretization data.

Let us consider here the simplest case of the Stokes equations, i.e. $\beta = 0$, with the stress tensor given by (2). The standard MAC scheme may then be written:

$$ \text{Find } (u_{\sigma})_{\sigma \in \mathcal{E}} \subset \mathbb{R}, (p_{K})_{K \in \mathcal{M}} \subset \mathbb{R}; \sum_{K \in \mathcal{M}} |K|p_{K} = 0, \\ \sum_{k=1}^{d} \sum_{\sigma \in \mathcal{E}_{K}^{(k)}} |\sigma| u_{\sigma} \mathbf{e}^{(k)} \cdot \mathbf{n}_{K, \sigma} = 0, \forall K \in \mathcal{M}, \qquad (7a) $$

$$ - \sum_{\sigma' \in \mathcal{E}^{(k)}} \tau_{\sigma, \sigma'}^{(k)} (u_{\sigma'} - u_{\sigma}) + |\sigma|(p_L - p_K) = \int_{V_{\sigma}^{(k)}} f^{(k)}(\mathbf{x}) d\mathbf{x}, \qquad (7b) $$

$$ \forall k = 1, \dots, d, \forall \sigma \in \mathcal{E}^{(k)}, \text{ and } K, L \in \mathcal{M} \text{ such that } \sigma \in \overrightarrow{K} | L. $$

Let us define $H_{\mathcal{M}}(\Omega)$ as the set of piecewise functions constant in $K \in \mathcal{M}$ (pressure unknown), and $H_{\mathcal{E}}^{(k)}(\Omega)$ as the set of piecewise functions which are constant in $V_{\sigma}$, for $\sigma \in \mathcal{E}^{(k)}$ (kth component of the velocity unknown). Let $H_{\mathcal{E}}(\Omega) = \{v = (v^{(k)})_{k=1,...,d}; v^{(k)} \in H_{\mathcal{E}}^{(k)}(\Omega)\}$. The piecewise constant discrete divergence of $v \in H_{\mathcal{E}}(\Omega)$ is the function $\text{div}_{\mathcal{M}}v \in H_{\mathcal{M}}(\Omega)$ defined by:

$$ \text{div}_{\mathcal{M}}v(\mathbf{x}) = \text{div}_K v = \frac{1}{|K|} \sum_{k=1}^{d} \sum_{\sigma \in \mathcal{E}_K} |\sigma| v_{K,\sigma}, \quad \text{for a.e. } \mathbf{x} \in K, \forall K \in \mathcal{M}, \forall v \in H_{\mathcal{E}}(\Omega), \qquad (8) $$

where $v_{K,\sigma} = v_{\sigma} \mathbf{n}_{K,\sigma} \cdot \mathbf{e}^{(k)}$. The discrete mass conservation equation (7a) can then be simply written as:

$$ \text{div}_{\mathcal{M}}u(\mathbf{x}) = 0, \text{ for a.e. } \mathbf{x} \in \Omega. \qquad (9) $$

We now introduce a variational formulation of the discrete momentum equation (7b) by recalling (see e.g. [7, Chapter 3] that if the points $x_{\sigma}$ are the nodes of a Delaunay triangulation, then

$$ \tau_{\sigma, \sigma'}^{(1)} = \int \nabla \xi_{\sigma}^{(1)}(\mathbf{x}) \cdot \nabla \xi_{\sigma'}^{(1)}(\mathbf{x}) d\mathbf{x}, $$

where the functions $(\xi_{\sigma}^{(1)})_{\sigma \in \mathcal{E}^{(1)}}$ are the $P^1$ finite element basis functions defined on the Delaunay triangulation, such as the one depicted in Figure 2 (in which the nodes are the points $x_{\sigma}$, for all $\sigma \in \mathcal{E}^{(1)}$). The
<figure><img src="image_3.png" /><figcaption>Figure 2: Triangular mesh for the $P^1$ finite element approximation $\hat{u}^{(1)}$ of $u^{(1)}$.</figcaption></figure>

functions $(\xi_{\sigma}^{(2)})_{\sigma \in \mathcal{E}^{(2)}}$ are defined in a similar way, and we may then define $\hat{u}^{(k)} \in H_0^1(\Omega)$ such that

$$ \hat{u}^{(k)} = \sum_{\sigma \in \mathcal{E}^{(k)}} u_{\sigma} \xi_{\sigma}^{(k)}, \quad \forall k = 1, \dots, d. \qquad (10) $$

We then define an inner product $\langle \cdot, \cdot \rangle_{\mathcal{E}}$ on $H_{\mathcal{E}}(\Omega)$ (this holds for both cases (2) and (3)) by

$$ \langle u, v \rangle_{\mathcal{E}} = \int_{\Omega} S(\hat{u}, \hat{v})(x)dx, \quad \forall u, v \in H_{\mathcal{E}}(\Omega), \qquad (11) $$

which we expect to approximate the inner product $\int_{\Omega} S(u, v)(x)dx$ of the continuous problem. We then obtain, multiplying (7) by $v_{\sigma}$ and summing on $k = 1, 2$ and $\sigma \in \mathcal{E}^{(k)}$,

$$ \langle u, v \rangle_{\mathcal{E}} - \int_{\Omega} p(x) \mathrm{div}_{\mathcal{M}} v(x) dx = \int_{\Omega} f(x) \cdot v(x) dx, \quad \forall v \in H_{\mathcal{E}}(\Omega), \qquad (12) $$

A discrete variational formulation of the standard MAC scheme (7) is therefore:

$$ \text{Find } u \in H_{\mathcal{E}}(\Omega) \text{ and } p \in H_{\mathcal{M}}(\Omega) \text{ s.t. } \sum_{K \in \mathcal{M}} |K|p_K = 0 \text{ and (9) and (12) hold.} \qquad (13) $$

## 3.2 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 Figure 6; note that these local refinements may be used to follow the contours of a general non rectangular domain, as in Figures 3 and 4.

We consider 2D or 3D meshes of $\Omega$, 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 $\mathbb{R}^d$, for some $k = 1, \dots, 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 $\partial\Omega$ 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 Figure 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 Figure 5.

We denote by $\mathcal{M}$ the set of pressure cells. Examples of resulting pressure grids are depicted in the left part of Figures 3 and 5. For all $K \in \mathcal{M}$, we again denote by $\mathcal{E}_K$ the set of all internal faces of $K$ (therefore, the faces of $K$ which are on $\partial\Omega$ are not elements of $\mathcal{E}_K$), and we define the set $\mathcal{E}$ as the union over $K \in \mathcal{M}$ of all the sets $\mathcal{E}_K$. It is assumed that a given $\sigma \in \mathcal{E}_{\text{int}}$ is entirely included in an interface between two cells, say $K$ and $L$; we shall write $\sigma \in K|L$. Note that an interface $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 $\mathcal{E}^{(k)}$ as the subset of $\mathcal{E}$ which contains all the internal faces whose normal is parallel to the basis vector $e^{(k)}$. For any $\sigma \in \mathcal{E}$, $x_{\sigma}$ denotes the center of gravity of $\sigma$.
In order to get a discrete variational MAC-like scheme, we consider, for any $k = 1, \dots, d$, the set of internal points $V_{\text{int}}^{(k)} = (x_{\sigma})_{\sigma \in \mathcal{E}^{(k)}}$ and a given family of external points $V_{\text{ext}}^{(k)}$ containing at least all the vertices of $\Omega$.

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 $\mathcal{T}^{(k)}$ of $\Omega$ whose vertices are $V_{\text{int}}^{(k)} \cup V_{\text{ext}}^{(k)}$. Such a triangulation of $\Omega$ is a set of simplices: triangles in 2D, tetrahedra in 3D. Each simplex has $d+1$ vertices which belong to $V_{\text{int}}^{(k)} \cup V_{\text{ext}}^{(k)}$. 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 \in \mathcal{T}^{(k)}$ does not contain any element of $V_{\text{int}}^{(k)} \cup V_{\text{ext}}^{(k)}$. Examples of Delaunay triangulations constructed from the edge mid-points are illustrated in Figures 3 and 5 (middle and right parts).

We then denote, for any $\sigma \in \mathcal{E}^{(k)}$, the function $\xi_{\sigma}^{(k)}$, which is continuous, piecewise $P^1$ on any $T \in \mathcal{T}^{(k)}$, and whose value is 1 at the point $x_{\sigma}$ and 0 at points $x_{\sigma'}$ for any $\sigma' \in \mathcal{E}^{(k)} \setminus \{\sigma\}$.

Let $\{V_{\sigma}^{(k)}, \sigma \in \mathcal{E}^{(k)}\}$ be the Voronoï mesh associated to the family $(x_{\sigma})_{\sigma \in \mathcal{E}^{(k)}}$, defining the Voronoï cells as follows:

$$V_{\sigma}^{(k)} = \{x \in \Omega, \mathrm{d}(x, x_{\sigma}) < \mathrm{d}(x, x_{\sigma'}), \sigma' \in \mathcal{E}^{(k)} \setminus \{\sigma\}, \forall \sigma \in \mathcal{E}^{(k)}\}$$

See Figures 4 and 6 for examples of the superposition of the Delaunay triangulations and the Voronoï mesh construction. We finally denote by $\mathcal{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_{\sigma}^{(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 Figure 6, where we depict the Delaunay grid and the Voronoï cells for the horizontal and vertical velocities.

<figure><img src="image_4.png" /><figcaption>Figure 3: The pressure (left) and velocity grids (middle: horizontal velocity, right: vertical velocity).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: The Delaunay triangulation and the Voronoï cells (left: horizontal velocity, right: vertical velocity).</figcaption></figure>

We may again define $H_{\mathcal{M}}(\Omega)$ as the set of piecewise functions constant on the pressure cells $K \in \mathcal{M}$, the set $H_{\mathcal{E}}^{(k)}(\Omega)$ of piecewise constant functions on the dual grid cells $V_{\sigma}$, for $\sigma \in \mathcal{E}^{(k)}$; this discrete set is the space of functions meant to approximate the $k$-th component of the velocity. We then denote by $H_{\mathcal{E}}(\Omega)$ the set of all $v = (v^{(k)})_{k=1,\dots,d}$ with $v^{(k)} \in H_{\mathcal{E}}^{(k)}(\Omega)$. We then define:

$$\hat{v}^{(k)} = \sum_{\sigma \in \mathcal{E}^{(k)}} v_{\sigma} \xi_{\sigma}^{(k)} \in H_{0}^{1}(\Omega), \quad (14)$$

we denote by $\hat{v} = (\hat{v}^{(k)})_{k=1,\dots,d}$, and we define the norm

$$\|v\|_{\mathcal{E}} = \|\nabla \hat{v}\|_{L^2(\Omega)^{d \times d}}, \forall v \in \mathcal{H}_{\mathcal{E}}(\Omega). \quad (15)$$

The extended MAC scheme for the Stokes equations ($\beta = 0$) is again (8)-(13), applying Definition (11) for $\langle \cdot, \cdot \rangle_{\mathcal{E}}$.
<figure><img src="image_6.png" /><figcaption>Figure 5: The pressure (left) and velocity grids (middle: horizontal velocity, right: vertical velocity).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: The Delaunay triangulation and the Voronoï cells (left: horizontal velocity, right: vertical velocity).</figcaption></figure>

# 4 The extended MAC scheme for the steady Navier-Stokes equations

## 4.1 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 $\int_{\Omega}(\mathbf{u}(\mathbf{x}) \cdot \nabla)\mathbf{u}(\mathbf{x}) \cdot \mathbf{v}(\mathbf{x})d\mathbf{x}$. To this purpose, we introduce a discrete trilinear form $b_{\varepsilon}$ which aims at discretizing the trilinear form $b$ defined over $(H_0^1(\Omega))^3$ by $b(\mathbf{u}, \mathbf{v}, \mathbf{w}) = \int_{\Omega}(\mathbf{u} \cdot \nabla \mathbf{v})\mathbf{w} \, d\mathbf{x}$. We begin by defining some interpolation operators between $H_{\varepsilon}(\Omega)$ and $(H_M(\Omega))^d$.

**Definition 4.1 (From $H_{\varepsilon}(\Omega)$ to $H_M(\Omega)$...and back)** Let $\mathcal{D}$ be a possibly non conforming mesh as depicted in Section 4. For $\mathbf{v} \in H_{\varepsilon}(\Omega)$, we define $\Pi_K \mathbf{v}$ by its components $(\Pi_K \mathbf{v})^{(k)}$:

$$ (\Pi_K \mathbf{v})^{(k)} = \frac{1}{\sum_{\sigma \in \mathcal{E}_K^{(k)}} |\sigma|} \sum_{\sigma \in \mathcal{E}_K^{(k)}} |\sigma| v_{\sigma}, \quad k = 1, \dots, d, \qquad (16) $$

and $\Pi_K \mathbf{v} \in H_M(\Omega) = (H_M(\Omega))^d$ as the piecewise constant function equal to $(\Pi_K \mathbf{v})$ on cell $K$. We then define $\Pi_{\varepsilon} \mathbf{v} \in H_{\varepsilon}(\Omega)$ as the following piecewise constant function on the Voronoï cells:

$$ \Pi_{\varepsilon} \mathbf{w} = \sum_{\substack{\sigma \in \mathcal{E}_{\text{int}} \\ \sigma \in K | L}} \frac{1}{2} (\Pi_K \mathbf{w} + \Pi_L \mathbf{w}) \mathbf{1}_{V_{\sigma}}, $$

where $\mathbf{1}_{V_{\sigma}}$ denotes the characteristic function of $V_{\sigma}$, that is $\mathbf{1}_{V_{\sigma}}(\mathbf{x}) = 1$ if $\mathbf{x} \in V_{\sigma}$ and 0 otherwise,

**Definition 4.2 (Discrete gradient and convection term)** For $\mathbf{v} \in H_M(\Omega)$, we define its discrete gradient $\nabla_{\varepsilon} \mathbf{v} \in H_{\varepsilon}(\Omega)$ by:

$$ \nabla_{\varepsilon} \mathbf{v} = (\partial_{\varepsilon}^{(1)} \mathbf{v}, \dots, \partial_{\varepsilon}^{(d)} \mathbf{v})^t, \quad \text{with} \quad \partial_{\varepsilon}^{(i)} \mathbf{v} = \sum_{\sigma \in \mathcal{E}^{(i)}} \partial_{\sigma} \mathbf{v} \mathbf{1}_{V_{\sigma}} \\ \text{and} \quad \partial_{\sigma} \mathbf{v} = \begin{cases} (v_L - v_K) \frac{|\sigma|}{|V_{\sigma}|}, & \text{if } \sigma \in \mathcal{E}_{\text{int}}, \sigma \in \overrightarrow{K|L}, \\ 0, & \text{if } \sigma \in \mathcal{E}_{\text{ext}}, \end{cases} \qquad (17) $$
where $\overrightarrow{\sigma \in K|L}$ means that $\sigma \in K|L \subset \mathcal{E}^{(k)}$ for some $k = 1, \dots, d$ and $\mathbf{n}_{K,\sigma} \cdot \mathbf{e}^{(k)} = 1$. We then define the following trilinear form $b_{\mathcal{E}}$ on $(\mathbf{H}_{\mathcal{E}}(\Omega))^{3}$ by:

$$b_{\mathcal{E}}(\mathbf{u}, \mathbf{v}, \mathbf{w}) = \int_{\Omega} (\mathbf{u} \cdot \nabla_{\mathcal{E}}) \Pi_{\mathcal{M}} \mathbf{v} \cdot \Pi_{\mathcal{E}} \mathbf{w} \, d\mathbf{x}, \text{ for } \mathbf{u}, \mathbf{v}, \mathbf{w} \in \mathbf{H}_{\mathcal{E}}(\Omega) \quad (18)$$

where, for any $\tilde{\mathbf{v}} \in \mathbf{H}_{\mathcal{M}}(\Omega)$,

$$(\mathbf{u} \cdot \nabla_{\mathcal{E}}) \tilde{\mathbf{v}}(\mathbf{x}) = \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{E}^{(i)}} u_{\sigma} \partial_{\sigma} \tilde{\mathbf{v}} \mathbf{1}_{V_{\sigma}}(\mathbf{x}).$$

It is well-known that the continuous trilinear form satisfies $b(\mathbf{u}, \mathbf{v}, \mathbf{v}) = -\frac{1}{2} \int_{\Omega} |\mathbf{v}|^2 \mathrm{div} \mathbf{u} \, d\mathbf{x}$ for all $\mathbf{u}, \mathbf{v} \in H_0^1(\Omega)^2$, and therefore, $b(\mathbf{u}, \mathbf{v}, \mathbf{v}) = 0$ if $\mathrm{div} \mathbf{u} = 0$. Similarly, we have the following result for the discrete trilinear form.

**Lemma 4.1 (Properties of the trilinear form)** Let $\mathbf{u}, \mathbf{v}, \mathbf{w} \in \mathbf{H}_{\mathcal{E}}(\Omega)$ and let $b_{\mathcal{E}}$ be defined by (18). Then:

$$b_{\mathcal{E}}(\mathbf{u}, \mathbf{v}, \mathbf{w}) = \sum_{K \in \mathcal{M}} \sum_{\substack{\sigma \in \mathcal{E}_K \\ \sigma \in \overline{K|L}}} |\sigma| u_{K,\sigma} \frac{\Pi_L \mathbf{v} - \Pi_K \mathbf{v}}{2} \cdot \Pi_K \mathbf{w}, \quad \forall \mathbf{u}, \mathbf{v}, \mathbf{w} \in \mathbf{H}_{\mathcal{E}}(\Omega), \quad (19)$$

where $u_{K,\sigma} = u_{\sigma} \mathbf{n}_{K,\sigma} \cdot \mathbf{e}^{(k)}$. Moreover,

$$b_{\mathcal{E}}(\mathbf{u}, \mathbf{v}, \mathbf{v}) = -\frac{1}{2} \int_{\Omega} |\Pi_{\mathcal{M}} \mathbf{v}|^2 \mathrm{div}_{\mathcal{M}} \mathbf{u}, \quad (20)$$

where the discrete divergence operator $\mathrm{div}_{\mathcal{M}}$ is defined by (8); therefore, if $\mathrm{div}_{\mathcal{M}} \mathbf{u} = 0$, $b_{\mathcal{E}}(\mathbf{u}, \mathbf{v}, \mathbf{v}) = 0$.

**Proof.** From (18) we get that

$$b_{\mathcal{E}}(\mathbf{u}, \mathbf{v}, \mathbf{w}) = \sum_{\substack{\sigma \in \mathcal{E} \\ \sigma \in \overline{K|L}}} u_{\sigma}(\Pi_L \mathbf{v} - \Pi_K \mathbf{v}) |\sigma| \frac{1}{2}(\Pi_K \mathbf{w} + \Pi_L \mathbf{w}). \\ = \frac{1}{2} \sum_{\substack{\sigma \in \mathcal{E} \\ \sigma \in \overline{K|L}}} u_{\sigma}(\Pi_L \mathbf{v} - \Pi_K \mathbf{v}) |\sigma| \Pi_K \mathbf{w} + \frac{1}{2} \sum_{\substack{\sigma \in \mathcal{E} \\ \sigma \in \overline{L|K}}} u_{\sigma}(\Pi_K \mathbf{v} - \Pi_L \mathbf{v}) |\sigma| \Pi_K \mathbf{w}.$$

Reordering the summation over the edges of each cell K yields (19). Taking $\mathbf{w} = \mathbf{v}$ in (19) then yields that

$$b_{\mathcal{E}}(\mathbf{u}, \mathbf{v}, \mathbf{v}) = \frac{1}{2} \sum_{K \in \mathcal{M}} \sum_{\substack{\sigma \in \mathcal{E}_K \\ \sigma \in K|L}} |\sigma| u_{K,\sigma} (\Pi_L \mathbf{v} \cdot \Pi_K \mathbf{v} - |\Pi_K \mathbf{v}|^2),$$

which yields (20) by conservativity. $\square$

The extended MAC scheme for the Navier-Stokes equation then reads:

Find $\mathbf{u} \in \mathbf{H}_{\mathcal{E}}(\Omega)$ and $p \in \mathbf{H}_{\mathcal{M}}(\Omega)$ s. t.

$$\sum_{K \in \mathcal{M}} |K| p_K = 0, \quad (21a)$$

$$\mathrm{div}_{\mathcal{M}} \mathbf{u}(\mathbf{x}) = 0, \text{ for a.e. } \mathbf{x} \in \Omega. \quad (21b)$$

$$\langle \mathbf{u}, \mathbf{v} \rangle_{\mathcal{E}} - \int_{\Omega} p(\mathbf{x}) \mathrm{div}_{\mathcal{M}} \mathbf{v}(\mathbf{x}) d\mathbf{x} + \beta b_{\mathcal{E}}(\mathbf{u}, \mathbf{u}, \mathbf{v}) = \int_{\Omega} f(\mathbf{x}) \cdot \mathbf{v}(\mathbf{x}) d\mathbf{x}, \forall \mathbf{v} \in \mathbf{H}_{\mathcal{E}}(\Omega). \quad (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 $\Pi_{\mathcal{E}} \mathbf{w}$ rather than $\mathbf{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 $\mathbf{u}$.
## 4.2 Mathematical analysis

Let $\mathcal{D}$ be a possibly non conforming mesh, such as introduced in Section 3.2. We define $h_{\mathcal{D}}$ as the maximum value of the diameters of all $T \in \mathcal{T}^{(k)}$ and all $V_{\sigma}^{(k)}$, $\sigma \in \mathcal{E}^{(k)}$, $k = 1, \dots, d$.

The regularity $\theta_{\mathcal{D}}$ of $\mathcal{D}$ is defined as the minimum value of:

1. all ratios $\text{diam}(T)/\text{diam}(T')$, for all neighbouring simplices $T, T' \in \mathcal{T}^{(k)}$, $k = 1, \dots, d$,

2. all ratios $|T|/\text{diam}(T)^d$, $T \in \mathcal{T}^{(k)}$, $k = 1, \dots, d$,

3. all ratios $|V_{\sigma}^{(k)}|/|T|$, $|T|/|V_{\sigma}^{(k)}|$, $\text{diam}(V_{\sigma}^{(k)})/\text{diam}(T)$, $\text{diam}(T)/\text{diam}(V_{\sigma}^{(k)})$, for any $T \in \mathcal{T}^{(k)}$ with vertex $x_{\sigma}^{(k)}$.

4. all ratios $|\sigma|\text{diam}(K)/|V_{\sigma}|$, $\sigma \in \mathcal{E}_{\text{int}}$, $K \in \mathcal{M}_{\sigma}$.

Let us then introduce the following interpolation operators.

**Definition 4.3 (Interpolation operators)** For $v = (v^{(1)}, \dots, v^{(d)}) \in H_0^1(\Omega)^d$, let $P_{\mathcal{E}}v \in H_{\mathcal{E}}(\Omega)$ be defined by the components $(P_{\mathcal{E}}v)^{(k)}$ of its piecewise constant values on the dual cells $V_{\sigma}^{(k)}$, for $k = 1, \dots, d$.

$$ (P_{\sigma}v)^{(k)} = \frac{1}{|\sigma|} \int_{\sigma} v^{(k)} \mathrm{d}s(x), \quad \forall \sigma \in \mathcal{E}^{(k)}. \qquad (22) $$

For $\varphi \in H_0^1(\Omega)$, let $P_{\mathcal{M}}\varphi \in H_{\mathcal{M}}(\Omega)$ be defined by its piecewise constant values on the cells $K \in \mathcal{M}$:

$$ (P_{\mathcal{M}}\varphi)_{K} = \varphi(x_{K}), \quad \forall K \in \mathcal{M}. \qquad (23) $$

Following the ideas of [21] applied to the piecewise $P^1$ interpolation $\widehat{P_{\mathcal{E}}v}$ of $v$, one may check that the interpolation operator $P_{\mathcal{E}}$ satisfies the following property.

**Lemma 4.2** Under the assumptions of Definition 4.3, there exists $C_1$, only depending on $\theta$ and $\Omega$, such that

$$ \|P_{\mathcal{E}}v\|_{\mathcal{E}} \le C_1 \|v\|_{H^1(\Omega)^d}. \qquad (24) $$

**Lemma 4.3 (Properties of the discrete divergence)** Let $v \in (H_0^1(\Omega))^d$, let $P_{\mathcal{E}}$ be the interpolation operator defined by (22). and let $p \in H_{\mathcal{M}}$. Then

$$ \int_{\Omega} p(x) \operatorname{div}_{\mathcal{M}} P_{\mathcal{E}} v(x) \mathrm{d}x = \int_{\Omega} p(x) \operatorname{div} v(x) \mathrm{d}x. \qquad (25) $$

Moreover, for any $u \in H_{\mathcal{E}}$, and for any $q \in H_{\mathcal{M}}(\Omega)$, one has:

$$ \int_{\Omega} q(x) \operatorname{div}_{\mathcal{M}} u \mathrm{d}x = - \int_{\Omega} \nabla_{\mathcal{E}} q(x) \cdot u(x) \mathrm{d}x. \qquad (26) $$

**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 \in \mathcal{M}$. The second result follows from the fact that

$$ \begin{aligned} \int_{\Omega} q(x) \operatorname{div}_{\mathcal{M}} u \mathrm{d}x &= \sum_{K \in \mathcal{M}} q_K \sum_{k=1}^{d} \sum_{\sigma \in \mathcal{E}_K^{(k)}} |\sigma| u_{\sigma} e^{(k)} \cdot n_{K, \sigma} \\ &= \sum_{\sigma \in \overline{K} | L} (q_K - q_L) |\sigma| u_{\sigma} \\ &= - \int_{\Omega} \nabla_{\mathcal{E}} q(x) \cdot u(x) \mathrm{d}x. \end{aligned} $$

□

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 Section 5.
**Lemma 4.4 (Estimates on the trilinear form)** Let $\mathcal{D}$ be a discretization in the sense given in Section 4.1 and let $b_{\mathcal{E}}$ be the trilinear form defined by (19). Then there exist $C_2$ and $C_3$, only depending on any $\theta \le \theta_{\mathcal{D}}$ and on $\Omega$, such that

$$|b_{\mathcal{E}}(u, v, w)| \le C_2 \|u\|_{(L^4(\Omega))^a} \|v\|_{\mathcal{E}} \|w\|_{(L^4(\Omega))^d} \le C_3 \|u\|_{\mathcal{E}} \|v\|_{\mathcal{E}} \|w\|_{\mathcal{E}}, \forall u, v, w \in \mathbf{H}_{\mathcal{E}}(\Omega). \quad (27)$$

Moreover, there exists $C_4$, only depending on $\Omega$ and $\theta$, such that

$$|b_{\mathcal{E}}(u, v, w)| \le C_4 \|u\|_{L^2(\Omega)}^{1/2} \|u\|_{\mathcal{E}}^{3/2} \|v\|_{\mathcal{E}}, \forall u \in \mathbf{H}_{\mathcal{E}}(\Omega); \operatorname{div}_{\mathcal{M}} u = 0, \forall v \in \mathbf{H}_{\mathcal{E}}(\Omega). \quad (28)$$

**Proof.** From the definition (18) of the trilinear form $b_{\mathcal{E}}$ denoting by $\tilde{v}^{(j)}$ (resp. $\tilde{w}^{(j)}$) the j-th component of $\Pi_{\mathcal{M}}v$ (resp. $\Pi_{\mathcal{E}}w$), we get thanks to the Cauchy–Schwarz inequality, we get that:

$$|b_{\mathcal{E}}(u, v, w)| \le \sum_{k=1}^{d} \sum_{j=1}^{d} \|u^{(k)}\|_{L^4(\Omega)} \|\partial_{\mathcal{E}}^{(k)} \tilde{v}^{(j)}\|_{L^2(\Omega)} \|\tilde{w}^{(j)}\|_{L^4(\Omega)} \\ \le d^2 \|u\|_{(L^4(\Omega))^d} \|\nabla_{\mathcal{E}}(\Pi_{\mathcal{M}}v)\|_{(L^2(\Omega))^{d \times d}} \|\Pi_{\mathcal{E}}w\|_{(L^4(\Omega))^d}.$$

Let us then remark that the components of $\frac{\Pi_{\mathcal{K}}w + \Pi_{\mathcal{L}}w}{2}$ are barycentric combinations of neighbouring terms $w_{\sigma}$; hence there exists $C_5$, only depending on $\Omega$ and $\theta$, such that

$$\|\Pi_{\mathcal{E}}w\|_{(L^4(\Omega))^d} \le C_5 \|w\|_{(L^4(\Omega))^d}.$$

In order to bound the term $\|\nabla_{\mathcal{E}}(\Pi_{\mathcal{M}}v)\|_{(L^2(\Omega))^{d \times d}}$, we first define the functions $\Pi_{-}^{(k)}v$ and $\Pi_{+}^{(k)}v$ by

$$\Pi_{-}^{(k)} v(x) = \Pi_{K} v, \forall \sigma \in \mathcal{E}^{(k)}, \text{ for a.e. } x \in V_{\sigma}^{(k)}, \text{ and } K, L \text{ such that } \sigma \in \overrightarrow{K|L},$$

and

$$\Pi_{+}^{(k)} v(x) = \Pi_{L} v, \forall \sigma \in \mathcal{E}^{(k)}, \text{ for a.e. } x \in V_{\sigma}^{(k)}, \text{ and } K, L \text{ such that } \sigma \in \overrightarrow{K|L}.$$

We finally define $\delta^{(k)}(x) = \frac{|V_{\sigma}^{(k)}|}{|\sigma|}$ for all $x \in V_{\sigma}^{(k)}$. We then get that

$$\|\partial_{\mathcal{E}}^{(k)} \Pi_{\mathcal{M}} v\|_{(L^2(\Omega))^d} = \int_{\Omega} \frac{1}{\delta^{(k)}(x)^2} \left| \Pi_{+}^{(k)} v(x) - \Pi_{-}^{(k)} v(x) \right|^2 dx \\ = \sum_{\ell=1}^{d} \int_{\Omega} \frac{1}{\delta^{(k)}(x)^2} \left( \Pi_{+}^{(k)} v^{(\ell)}(x) - \Pi_{-}^{(k)} v^{(\ell)}(x) \right)^2 dx.$$

In order to apply Lemma 7.1, we first remark that, thanks to the regularity hypotheses of the mesh, there exists some $C_6 > 0$, only depending on $\theta$, such that $\delta^{(k)}(x) \ge C_6\delta(x)$, with $\delta(x)$ defined in Lemma 7.1. Since the functions $\Pi_{-}^{(k)}v^{(\ell)}(x)$ only depend on the discrete unknowns $v_{\sigma}^{\ell}$, we may introduce the functions $\psi_{-, \sigma}^{k, \ell}$ such that

$$\Pi_{-}^{(k)} v^{(\ell)}(x) = \sum_{\sigma \in \mathcal{E}^{(\ell)}} v_{\sigma}^{\ell} \psi_{-, \sigma}^{k, \ell}(x),$$

which are the piecewise constant functions defined by:

$$\forall \sigma \in \mathcal{E}^{(k)}, \text{ with } K, L \text{ such that } \sigma \in \overrightarrow{K|L}, \\ \psi_{-, \sigma}^{k, \ell}(x) = \begin{cases} \frac{|\sigma'|}{\sum_{\sigma'' \in \mathcal{E}^{(\ell)}} |\sigma''|} & \text{if } \sigma' \in \mathcal{E}^{(\ell)} \setminus \mathcal{E}_K, \text{ for a.e. } x \in V_{\sigma}^{(k)}, \\ 0 & \text{if } \sigma' \in \mathcal{E}^{(\ell)} \setminus \mathcal{E}_K \text{ or if } x \notin V_{\sigma}^{(k)}. \end{cases}$$

The functions $\psi_{-, \sigma}^{k, \ell}$ satisfy the hypothesis (49) of Lemma 7.1; furthermore, thanks to the regularity of the mesh, the mesh dependent bound $C_{28}$ defined by (50) remains bounded by a function of $\theta$. We may thus apply Lemma 7.1, which yields (51). This implies

$$\int_{\Omega} \frac{1}{\delta^{(k)}(x)^2} \left( \Pi_{-}^{(k)} v^{(\ell)}(x) - \hat{v}^{(\ell)}(x) \right)^2 dx \le \frac{C_{29}}{(C_6)^2} \|\nabla \hat{v}^{(\ell)}\|_{L^2(\Omega)^d}^2.$$
Similarly, we have the same inequality with $\Pi_{+}^{(k)}v^{(\ell)}$ instead of $\Pi_{-}^{(k)}v^{(\ell)}$, which shows the existence of $C_7 > 0$, only depending on $\theta$, such that

$$ \|\nabla_{\mathcal{E}}(\Pi_{\mathcal{M}}v)\|_{(L^2(\Omega))^{d \times d}} \leq C_7 \|\nabla \hat{v}\|_{L^2(\Omega)^{d \times d}}, $$

thus concluding the left inequality of (27).

Now, thanks to the equivalence of the norms proved in Lemma 7.2 in the appendix (inequality (55)), there exists $C_8 \in \mathbb{R}_+$, only depending on any $\theta \le \theta_D$, such that

$$ \|u^{(k)}\|_{L^4(\Omega)} \le C_8 \|\hat{u}^{(k)}\|_{L^4(\Omega)}, $$

Applying the standard Sobolev inequality $\|\hat{u}^{(k)}\|_{L^4(\Omega)} \le C_{\text{sob}} \|\nabla \hat{u}^{(k)}\|_{L^2(\Omega)^d}$, where $C_{\text{sob}}$ only depends on $\Omega$ and $d$, we conclude that there exists $C_9$, only depending on $\Omega$ and on any $\theta \le \theta_D$, such that

$$ \|u^{(k)}\|_{L^4(\Omega)} \le C_9 \|\nabla \hat{u}\|_{L^2(\Omega)^{d \times d}}. $$

Therefore, we conclude the right inequality of (27).

Let us then prove (28); since $\text{div}_{\mathcal{M}}u = 0$, we have

$$ b_D(u, u, v) = -b_D(u, v, u), $$

which proves, using (27), that

$$ |b_D(u, u, v)| \le C_2 \|u\|_{L^4(\Omega)}^2 \|v\|_{\mathcal{E}}. $$

Using the Cauchy–Schwarz inequality, we have:

$$ \|u\|_{L^4(\Omega)}^2 \le \|u\|_{L^2(\Omega)}^{1/2} \|u\|_{L^6(\Omega)}^{3/2}, $$

We again apply Lemma 7.2: there exists $C_{10} \in \mathbb{R}_+$, only depending on $\theta$ and $\Omega$, such that

$$ \|u\|_{L^6(\Omega)} \le C_{10} \|\hat{u}\|_{L^6(\Omega)}. $$

Applying the standard Sobolev inequality, we get that $\|\hat{u}\|_{L^6(\Omega)} \le C_{\text{sob}} \|\nabla \hat{u}\|_{L^2(\Omega)^{d \times d}}$, where $C_{\text{sob}} \in \mathbb{R}_+$ depends only on $\Omega$ and $d$; therefore, by the definition (15) of the norm $\| \cdot \|_{\mathcal{E}}$, we get:

$$ |b_D(u, u, v)| \le C \|u\|_{L^2(\Omega)}^{1/2} \|u\|_{\mathcal{E}}^{3/2} \|v\|_{\mathcal{E}}. $$

which concludes the proof of Lemma 4.4.

<figure><img src="image_8.png" /></figure>

**Lemma 4.5 (Estimates on the velocity and the pressure for the Navier-Stokes problem)** Let $\mathcal{D}$ be a discretization in the sense given in Section 4.1 and let $\theta \le \theta_D$. Let $(u, p) \in H_{\mathcal{E}}(\Omega) \times H_{\mathcal{M}}(\Omega)$ be a solution to (21). Then there exists $C_{11}$, only depending on $\Omega$ (and on $\lambda$ and $\mu$ in the case (3), this dependency is no longer mentioned in this paper), such that $u$ satisfies the following estimates:

$$ \|u\|_{\mathcal{E}} \le C_{11} \|f\|_{(L^2(\Omega))^d}. \qquad (29) $$

Moreover, there exists $C_{12}$ only depending on $\beta$, $\theta$ and $\Omega$ such that

$$ \|p\|_{L^2(\Omega)} \le C_{12} \|f\|_{(L^2(\Omega))^d}. \qquad (30) $$

As a direct consequence, there exists one and only one solution (if $\beta = 0$) and at least one solution (if $\beta > 0$) to the scheme (21).

**Proof.** We let $v = u$ in (21). Let us first remark that we have $\int_{\Omega} p(x) \text{div}_{\mathcal{M}} u(x) dx = 0$ and that by the definition (19) of $b_{\mathcal{E}}$, we have

$$ b_{\mathcal{E}}(u, u, u) = - \sum_{K \in \mathcal{M}} \frac{1}{2} |\Pi_K u|^2 \sum_{k=1}^{d} \sum_{\sigma \in \mathcal{E}_K^{(k)}} |\sigma| u_{\sigma} e^{(k)} \cdot n_{K, \sigma} = 0, $$
since $\mathrm{div}_{\mathcal{M}}u(\boldsymbol{x}) = 0$. Hence we obtain

$$C_{13} \|\boldsymbol{u}\|_{\varepsilon}^{2} \leq \|S(\boldsymbol{u}, \boldsymbol{u})\|_{L^{1}(\Omega)} = \int_{\Omega} \boldsymbol{f}(\boldsymbol{x}) \cdot \boldsymbol{u}(\boldsymbol{x}) \mathrm{d}\boldsymbol{x},$$

where we denote by $C_{13} > 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,

$$\|\boldsymbol{u}\|_{\varepsilon}^{2} \leq \|\boldsymbol{f}\|_{L^{2}(\Omega)^{d}} \|\boldsymbol{u}\|_{L^{2}(\Omega)^{d}}.$$

Let $\hat{\boldsymbol{u}}$ be defined by (10). Applying Lemma 7.2, we get that there exists $C_{14}$, only depending on $\theta$, such that

$$\|\boldsymbol{u}\|_{L^2(\Omega)^d} \leq C_{14} \|\hat{\boldsymbol{u}}\|_{L^2(\Omega)^d}.$$

Thanks to the Poincaré inequality, we may write

$$\|\hat{\boldsymbol{u}}\|_{L^2(\Omega)^d} \leq \mathrm{diam}(\Omega) \|\boldsymbol{u}\|_{\varepsilon}. \qquad (31)$$

Therefore we get (29).

Let us now turn to the proof of (30). We use the following property, first due to Nečas [15], see also [2]: there exists $\boldsymbol{v} \in H_0^1(\Omega)^d$ such that $p(\boldsymbol{x}) = \mathrm{div}\boldsymbol{v}(\boldsymbol{x})$ for a.e. $\boldsymbol{x} \in \Omega$ and there exists $C_{15}$, only depending on $\Omega$, such that

$$\|\boldsymbol{v}\|_{H^1(\Omega)^d} \leq C_{15} \|p\|_{L^2(\Omega)}. \qquad (32)$$

Let $P_{\varepsilon}$ be the interpolation operator defined by (22). By Lemma 4.3, we have:

$$\int_{\Omega} p(\boldsymbol{x}) \mathrm{div}_{\mathcal{M}} P_{\varepsilon} \boldsymbol{v}(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} = \int_{\Omega} p(\boldsymbol{x}) \mathrm{div} \boldsymbol{v}(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} = \int_{\Omega} p(\boldsymbol{x})^2 \mathrm{d}\boldsymbol{x}.$$

We then get, using $P_{\varepsilon}\boldsymbol{v} \in \mathcal{H}_{\varepsilon}(\Omega)$ as test function in Scheme (21),

$$\langle \boldsymbol{u}, P_{\varepsilon} \boldsymbol{v} \rangle_{\varepsilon} - \int_{\Omega} p(\boldsymbol{x})^2 \mathrm{d}\boldsymbol{x} + \beta b_{\varepsilon}(\boldsymbol{u}, \boldsymbol{u}, P_{\varepsilon} \boldsymbol{v}) = \int_{\Omega} \boldsymbol{f}(\boldsymbol{x}) \cdot P_{\varepsilon} \boldsymbol{v}(\boldsymbol{x}) \mathrm{d}\boldsymbol{x}. \qquad (33)$$

Thanks to the Cauchy–Schwarz inequality, we get that there exists $C_{16}$ (only depending on $\lambda$ and $\mu$ in the case (3), this dependency is no longer recalled), such that

$$\langle \boldsymbol{u}, P_{\varepsilon} \boldsymbol{v} \rangle_{\varepsilon} \leq C_{16} \|\boldsymbol{u}\|_{\varepsilon} \|P_{\varepsilon} \boldsymbol{v}\|_{\varepsilon}.$$

Using Inequality (24), we get that there exists $C_{17}$, only depending on $\theta$ and $\Omega$, such that

$$\langle \boldsymbol{u}, P_{\varepsilon} \boldsymbol{v} \rangle_{\varepsilon} \leq C_{17} \|\boldsymbol{u}\|_{\varepsilon} \|\boldsymbol{v}\|_{H^1(\Omega)^d} \leq C_{17} C_{11} \|\boldsymbol{f}\|_{(L^2(\Omega))^d} C_{15} \|p\|_{L^2(\Omega)}.$$

We may also write, thanks to the Poincaré inequality and to (32),

$$- \int_{\Omega} \boldsymbol{f}(\boldsymbol{x}) \cdot P_{\varepsilon} \boldsymbol{v}(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} \leq \mathrm{diam}(\Omega) C_{15} \|\boldsymbol{f}\|_{(L^2(\Omega))^d} \|p\|_{L^2(\Omega)}.$$

Therefore, using (29) and (27), we get the existence of $C_{18}$ such that

$$\|p\|_{L^2(\Omega)}^2 \leq C_{18} \left( \|\boldsymbol{f}\|_{(L^2(\Omega))^d} + \|\boldsymbol{f}\|_{(L^2(\Omega))^d}^2 \right) \|p\|_{L^2(\Omega)}.$$

which concludes the proof of (30).

In the case where $\beta = 0$, the inequalities (29) and (30) suffice to prove that the square linear system issued from Scheme (21) (replacing one equation of (9) by $\int_{\Omega} p(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} = 0$) is invertible. If $\beta > 0$, using the topological degree argument as done in [8], we conclude to the existence of at least one solution to Scheme (21). $\square$

The next lemma shows a property of weak convergence of the discrete gradient which is used in the passage to the limit in the scheme.
**Lemma 4.6 (Weak convergence of the gradient)** Let $\mathcal{F}$ be a family of discretizations in the sense given in Section 4.1. For any $\mathcal{D} \in \mathcal{F}$, let $v_D \in H_M(\Omega)$, such that:

$$ \text{the family } (\nabla_{\mathcal{E}} v_D)_{\mathcal{D} \in \mathcal{F}} \text{ is bounded in } (L^2(\Omega))^d, \quad (34) $$

$$ \text{the family } (v_D)_{\mathcal{D} \in \mathcal{F}} \text{ converges weakly in } L^2(\Omega) \text{ to } v \in L^2(\Omega) \text{ as } h_D \to 0. \quad (35) $$

Then $v \in H_0^1(\Omega)$, and the family $(\nabla_{\mathcal{E}} v_D)_{\mathcal{D} \in \mathcal{F}}$ converges to $\nabla v$ weakly in $L^2(\Omega)$ as $h_D$ tends to 0.

**Proof.** Let $\psi = (\psi_1, \dots, \psi_d) \in (C^\infty(\mathbb{R}^d))^d$; then

$$ \int_{\Omega} \nabla_{\mathcal{E}} v_D(\mathbf{x}) \cdot \psi(\mathbf{x}) \mathrm{d}\mathbf{x} = \int_{\Omega} \nabla_{\mathcal{E}} v_D(\mathbf{x}) \cdot P_{\mathcal{E}} \psi(\mathbf{x}) \mathrm{d}s(\mathbf{x}) + R_{\mathcal{E}} $$

where $P_{\mathcal{E}}$ is the interpolation operator defined by (22), and $|R_{\mathcal{E}}| \le C_{\psi} h_D \|\nabla_{\mathcal{E}} v_D\|_{(L^2(\Omega))^2} \to 0$ as $h_D \to 0$, thanks to Assumption (34). Now, by (26) and (25) of Lemma 4.3,

$$ \int_{\Omega} \nabla_{\mathcal{E}} v_D(\mathbf{x}) \cdot P_{\mathcal{E}} \psi(\mathbf{x}) \mathrm{d}s(\mathbf{x}) = - \int_{\Omega} v_D(\mathbf{x}) \mathrm{div} \psi(\mathbf{x}) \mathrm{d}\mathbf{x}. $$

Therefore, by Assumption (35),

$$ \int_{\Omega} \nabla_{\mathcal{E}} v_D(\mathbf{x}) \cdot P_{\mathcal{E}} \psi(\mathbf{x}) \mathrm{d}s(\mathbf{x}) \to - \int_{\Omega} v(\mathbf{x}) \mathrm{div} \psi(\mathbf{x}) \mathrm{d}\mathbf{x}, $$

which shows that $\nabla_{\mathcal{E}} v_D$ tends to $\nabla v$ in the distribution sense. Since the sequence $(\nabla_{\mathcal{E}})_{\mathcal{D} \in \mathcal{F}}$ is bounded in $L^2(\Omega)$, we have that $\nabla v \in L^2(\Omega)$, and prolonging $v_D$ by 0 outside of $\Omega$, we get that $v = 0$ in $\mathbb{R}^d \setminus \Omega$. Therefore, $\nabla_{\mathcal{E}} v_D$ tends to $\nabla v \in H_0^1(\Omega)$ weakly in $(L^2(\Omega))^d$. $\square$

**Theorem 4.1 (Convergence of the scheme)** Let $\mathcal{F}$ be a family of discretizations in the sense given in Section 4.1, such that there exists $\theta > 0$ with $\theta \le \theta_D$ for all $\mathcal{D} \in \mathcal{F}$. For any $\mathcal{D} = (\mathcal{M}, \mathcal{E}) \in \mathcal{F}$, let $(u_{\mathcal{E}}, p_{\mathcal{M}}) \in H_{\mathcal{E}}(\Omega) \times H_D(\Omega)$ denote a solution to Scheme (21). Then there exists a weak solution $(u, p)$ of (1) such that, up to a subsequence, $u_{\mathcal{E}}$ converges in $L^2(\Omega)^d$ to $u$ and $p_{\mathcal{M}}$ converges in $L^2(\Omega)^d$ to $p$ as $h_D \to 0$. Moreover, if $\beta = 0$, the whole family converges to the unique weak solution of (1) as $h_D \to 0$.

**Proof.** Let $(D^{(\ell)})_{\ell \in \mathbb{N}}$ be a sequence of elements of $\mathcal{F}$, such that $h_{D^{(\ell)}}$ tends to 0 as $\ell \to \infty$ and such that there exists $\theta > 0$ with $\theta_{D^{(\ell)}} \ge \theta$, for all $m \in \mathbb{N}$. Let $(u_{\mathcal{E}^{(\ell)}}, p_{\mathcal{M}^{(\ell)}})$ denote a solution of (21) for the discretization $D_{\ell} = (\mathcal{M}_{\ell}, \mathcal{E}_{\ell}, \partial_{D_{\ell}})$. From estimate (29), using Rellich's theorem, we deduce that there exists $u \in H_0^1(\Omega)^d$ (with $u = (u^{(k)})_{k=1,\dots,d}$) and a subsequence of $(D^{(\ell)})_{m \in \mathbb{N}}$, again denoted by $(D^{(\ell)})_{\ell \in \mathbb{N}}$, such that $\tilde{u}_{\mathcal{M}^{(\ell)}} \in H_0^1(\Omega)$ weakly converges in $H_0^1(\Omega)^d$ (therefore strongly in $L^2(\Omega)^d$) to $u^{(k)}$ as $\ell \to \infty$ for $k = 1, \dots, d$. Using Lemma 7.1 and the estimate (29), we get that $u_{\mathcal{E}^{(\ell)}}$ also converges in $L^2(\Omega)^d$ to $u$ as $\ell \to \infty$. Thanks to the estimate (30) on the pressure, we may then consider a subsequence of $(D^{(\ell)})_{\ell \in \mathbb{N}}$, again denoted by $(D^{(\ell)})_{\ell \in \mathbb{N}}$, such that $p_{\mathcal{M}^{(\ell)}}$ weakly converges in $L^2(\Omega)$ to some function $p \in L^2(\Omega)$ as $\ell \to \infty$, with $\int_{\Omega} p(\mathbf{x}) \mathrm{d}\mathbf{x} = 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 \in E(\Omega)$, i.e. $\mathrm{div} u(\mathbf{x}) = 0$ for a.e. $\mathbf{x} \in \Omega$. Let $\varphi \in C_c^\infty(\Omega)$; multiplying (21b) by $P_{\mathcal{M}}\varphi$ and integrating over $\Omega$ yields, thanks to Lemma 4.3:

$$ \int_{\Omega} P_{\mathcal{M}}\varphi(\mathbf{x}) \cdot \mathrm{div}_{\mathcal{M}} u_{\mathcal{E}}(\mathbf{x}) \mathrm{d}\mathbf{x} = \int_{\Omega} \nabla_{\mathcal{E}} P_{\mathcal{M}}\varphi(\mathbf{x}) \cdot u_{\mathcal{E}}(\mathbf{x}) \mathrm{d}\mathbf{x} = 0. $$

Thanks to the regularity of the mesh, there exists $C_{19}$ depending only on $\varphi$ and $\theta$ such that $\|\nabla_{\mathcal{E}} P_{\mathcal{M}}\varphi\|_{(L^2(\Omega))^d} \le C_{19}$. Therefore, by Lemma 4.6, $\nabla_{\mathcal{E}} P_{\mathcal{M}}\varphi$ tends to $\nabla\varphi$ weakly in $(L^2(\Omega))^d$ as $h_D$ tends to 0. Hence, passing to the limit in the above equation, we get that $\int_{\Omega} \nabla\varphi(\mathbf{x}) \cdot u(\mathbf{x}) \mathrm{d}\mathbf{x} = 0$ for any $\varphi \in C_c^\infty(\Omega)$, which shows that $\mathrm{div} u(\mathbf{x}) = 0$ for a.e. $\mathbf{x} \in \Omega$.

Let $\varphi \in C_c^\infty(\Omega)^d$. Let us now show (1) with $v = \varphi$. We take (again omitting some indices $(\ell)$, thus denoting $\mathcal{D} = D^{(\ell)}$) $v_{\mathcal{E}} = P_{\mathcal{E}}\varphi \in H_{\mathcal{E}}(\Omega)$, as defined by (22), as test function in Scheme (21). We get that

$$ \langle u_D, P_{\mathcal{E}}\varphi \rangle_{\mathcal{E}} - \int_{\Omega} p_D(\mathbf{x}) \mathrm{div}_{\mathcal{M}} P_{\mathcal{E}}\varphi(\mathbf{x}) \mathrm{d}\mathbf{x} + \beta \, b_{\mathcal{E}}(u_D, u_D, P_{\mathcal{E}}\varphi) = \int_{\Omega} f(\mathbf{x}) \cdot P_{\mathcal{E}}\varphi(\mathbf{x}) \mathrm{d}\mathbf{x}. $$
Using Lemma 4.3 and the fact that $p_{D(\ell)}$ weakly converges in $L^2(\Omega)$ to $p \in L^2(\Omega)$, we obtain that

$$ \lim_{\ell \to \infty} \int_{\Omega} p_{D(\ell)}(\boldsymbol{x}) \mathrm{div}_{D(\ell)} P_{\mathcal{E}(\ell)} \varphi(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} = \int_{\Omega} p(\boldsymbol{x}) \mathrm{div} \varphi(\boldsymbol{x}) \mathrm{d}\boldsymbol{x}. $$

Thanks to the weak and strong convergence properties of the different sequences, we get that

$$ \lim_{\ell \to \infty} \langle \boldsymbol{u}_{D(\ell)}, P_{\mathcal{E}(\ell)} \varphi \rangle_{\mathcal{E}(\ell)} = \lim_{\ell \to \infty} \int_{\Omega} \mathcal{S}(\widehat{\boldsymbol{u}}_{D(\ell)}, \widehat{P}_{\mathcal{E}(\ell)} \varphi)(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} = \int_{\Omega} \mathcal{S}(\boldsymbol{u}, \varphi)(\boldsymbol{x}) \mathrm{d}\boldsymbol{x}. $$

We also have, thanks to the definition (22) of $P_{\mathcal{E}}\varphi$,

$$ \lim_{\ell \to \infty} \int_{\Omega} \boldsymbol{f}(\boldsymbol{x}) \cdot P_{\mathcal{E}(\ell)} \varphi(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} = \int_{\Omega} \boldsymbol{f}(\boldsymbol{x}) \cdot \varphi(\boldsymbol{x}) \mathrm{d}\boldsymbol{x}. $$

Let us now turn to the study of the limit of $b_{\mathcal{E}(\ell)}(\boldsymbol{u}_{D(\ell)}, \boldsymbol{u}_{D(\ell)}, P_{\mathcal{E}(\ell)}\varphi)$. By the definition (19) of $b_{\mathcal{E}}$, again dropping some indices $(\ell)$, we have:

$$ b_{\mathcal{E}}(\boldsymbol{u}_{D}, \boldsymbol{u}_{D}, P_{\mathcal{E}}\varphi) = \int_{\Omega} (\boldsymbol{u}_{\mathcal{M}}(\boldsymbol{x}) \cdot \nabla_{\mathcal{E}})(\boldsymbol{\Pi}_{D}\boldsymbol{u}_{\mathcal{E}})(\boldsymbol{x}) \cdot \boldsymbol{\Pi}_{D}(P_{\mathcal{E}}\varphi(\boldsymbol{x}))\mathrm{d}\boldsymbol{x}, $$

Thanks to the regularity of the mesh, the family $\nabla_{\mathcal{E}}(\boldsymbol{\Pi}_{D}\boldsymbol{u}_{\mathcal{E}})$ is bounded in $(L^2(\Omega))^d$. Moreover, using Lemma 7.1 and the estimate (29), we get that there exists $C_{20}$, only depending on $\theta$, such that

$$ \|\boldsymbol{\Pi}_{D}\boldsymbol{u}_{D} - \boldsymbol{u}_{D}\|_{L^2(\Omega)^d} \leq h_{D}C_{20}\|\boldsymbol{u}_{D}\|_{\mathcal{E}}, $$

which proves that $\boldsymbol{\Pi}_{D(\ell)}\boldsymbol{u}_{D(\ell)}$ converges to $\boldsymbol{u}$ in $L^2(\Omega)^d$. Therefore, thanks to Lemma 4.6, that $\nabla_{\mathcal{E}(\ell)}\boldsymbol{u}_{D(\ell)}$ converges weakly in $(L^2(\Omega))^d$ to $\nabla\boldsymbol{u}$. Since we have that $\boldsymbol{\Pi}_{D(\ell)}\varphi$ converges in $(L^\infty(\Omega))^d$ to $\varphi$, we obtain the convergence of $b_{\mathcal{E}(\ell)}(\boldsymbol{u}_{D(\ell)}, \boldsymbol{u}_{D(\ell)}, P_{\mathcal{E}(\ell)}\varphi)$ to $\int_{\Omega}(\boldsymbol{u}(\boldsymbol{x}) \cdot \nabla)\boldsymbol{u}(\boldsymbol{x}) \cdot \varphi(\boldsymbol{x})\mathrm{d}\boldsymbol{x}$ as $\ell \to \infty$, which completes the proof that (1) holds for all $\boldsymbol{v} \in H_0^1(\Omega)^d$ by density. $\square$

# 5 The extended MAC scheme for the time-dependent Navier-Stokes equations

## 5.1 Definition of the scheme

Let $\mathcal{D}$ be the collection of all the space discretization data as given in Section 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 $\boldsymbol{u}^{(0)} \in \mathcal{H}_{\mathcal{E}}(\Omega)$ defined by

$$ u_{\sigma}^{(0)} = \frac{1}{|V_{\sigma}^{(k)}|} \int_{V_{\sigma}^{(k)}} u_{\mathrm{ini}}^{(k)}(\boldsymbol{x}) \mathrm{d}\boldsymbol{x}, \quad \forall \sigma \in \mathcal{E}^{(k)}, \qquad (36) $$

and, for a given time step $\delta t > 0$, let us denote by $t^{(n)} = n\delta t$ for $n \ge 0$. The scheme is defined by

$$ \delta_{D, \delta t}^{(n+\frac{1}{2})} \boldsymbol{u}(\boldsymbol{x}) = \frac{\boldsymbol{u}^{(n+1)}(\boldsymbol{x}) - \boldsymbol{u}^{(n)}(\boldsymbol{x})}{\delta t}, \quad \text{for a.e. } \boldsymbol{x} \in \Omega, \qquad (37) $$

$$ \mathrm{div}_{\mathcal{M}} \boldsymbol{u}^{(n+\frac{1}{2})}(\boldsymbol{x}) = 0, \quad \text{for a.e. } \boldsymbol{x} \in \Omega, \qquad (38) $$

and, for a given $\alpha \in [\frac{1}{2}, 1]$, for all $n \in \mathbb{N}$, $(\boldsymbol{u}^{(n+1)}, p^{(n+\frac{1}{2})}) \in \mathcal{H}_{\mathcal{E}}(\Omega) \times \mathcal{H}_{\mathcal{M}}(\Omega)$,

$$ \int_{\Omega} p^{(n+\frac{1}{2})}(\boldsymbol{x})\mathrm{d}\boldsymbol{x} = 0, \quad \boldsymbol{u}^{(n+\frac{1}{2})} = (1-\alpha)\boldsymbol{u}^{(n)} + \alpha\boldsymbol{u}^{(n+1)}, \qquad (39) $$

$$ \int_{\Omega} \delta_{D, \delta t}^{(n+\frac{1}{2})} \boldsymbol{u}(\boldsymbol{x}) \cdot \boldsymbol{v}(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} + \langle \boldsymbol{u}^{(n+\frac{1}{2})}, \boldsymbol{v} \rangle_{\mathcal{E}} - \int_{\Omega} p^{(n+\frac{1}{2})}(\boldsymbol{x}) \mathrm{div}_{\mathcal{M}} \boldsymbol{v}(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} \\ + \beta b_{\mathcal{D}}(\boldsymbol{u}^{(n+\frac{1}{2})}, \boldsymbol{u}^{(n+\frac{1}{2})}, \boldsymbol{v}) = \frac{1}{\delta t} \int_{t^{(n)}}^{t^{(n+1)}} \int_{\Omega} \boldsymbol{f}(\boldsymbol{x}, t) \cdot \boldsymbol{v}(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} \mathrm{d}t, \quad \forall \boldsymbol{v} \in \mathcal{H}_{\mathcal{E}}(\Omega). $$
We then denote

$$
\begin{gathered}
\hat{u}_{D,\partial t}(\mathbf{x}, t) = \hat{u}^{(n+\frac{1}{2})}(\mathbf{x}), \quad \mathbf{u}_{D,\partial t}(\mathbf{x}, t) = \mathbf{u}^{(n+\frac{1}{2})}(\mathbf{x}), \quad \text{and } p_{D,\partial t}(\mathbf{x}, t) = p^{(n+\frac{1}{2})}(\mathbf{x}), \\
\text{for a.e. } (\mathbf{x}, t) \in \Omega \times (t^{(n)}, t^{(n+1)}), \quad \forall n \in \mathbb{N}, \quad (40)
\end{gathered}
$$

and we define the discrete time derivative

$$
\delta_{D,\partial t} \mathbf{u}(\mathbf{x}, t) = \delta_{D,\partial t}^{(n+\frac{1}{2})} \mathbf{u}(\mathbf{x}) \text{ for a.e. } (\mathbf{x}, t) \in \Omega \times (t^{(n)}, t^{(n+1)}), \forall n \in \mathbb{N}.
$$

**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.

**5.2 Mathematical analysis**

**Lemma 5.1 (Estimates on the velocity and the pressure)** Let $\partial t = T/N_T$ with $N_T \in \mathbb{N}^*$, let $\mathcal{D}$ be a discretization as defined in Section 4.1. Let the initial condition $\mathbf{u}^{(0)}$ be defined by (36), and, for a given $\alpha \in [\frac{1}{2}, 1]$, let $(\mathbf{u}^{(n+1)}, \mathbf{u}^{(n+\frac{1}{2})}, p^{(n+\frac{1}{2})}) \in \mathbf{H}_{\mathcal{E}}(\Omega) \times \mathbf{H}_{\mathcal{M}}(\Omega) \times \mathbf{H}_{\mathcal{N}}(\Omega)$ such that (38)-(39) holds for $n=1, \dots, N_T$. Then there exists $C_{21}$, only depending on $f$, $\mathbf{u}_{\text{ini}}$, $T$ and $\Omega$, such that the following discrete $L^2(H^1)$ and $L^\infty(L^2)$ estimates hold:

$$
\sum_{n=0}^{N_T-1} \partial_t \| \mathbf{u}^{(n+\frac{1}{2})} \|_{\mathcal{E}}^2 \le C_{21}, \quad (41)
$$

$$
\text{and } \|\mathbf{u}^{(n)}\|_{L^2(\Omega)^d} \le C_{21}, n = 0, \dots, N_T. \tag{42}
$$

Moreover we have the following $L^1(L^2)$ estimate on the pressure:

$$
|p_{D,\partial t}|_{L^1(0,T;L^2(\Omega))} \le C_{22}, \quad n=0,\dots,N, \tag{43}
$$

where $C_{22}$ only depends on $\mathcal{D}$, $\beta$, $\partial t$, $f$, $\mathbf{u}_{\text{ini}}$ and $\Omega$. Therefore there exists at least one solution to (38)-(39) if $\beta > 0$ and exactly one solution if $\beta = 0$.

**Proof.** Let us first remark that (42), for $n = 0$, is due to the Cauchy–Schwarz inequality applied to (36). We let $\mathbf{u} = \mathbf{u}^{(n+\frac{1}{2})}$ in (39). The nonlinear term again vanishes, which leads, using the Young, Poincaré and Korn inequalities, to:

$$
\begin{aligned}
& \frac{1}{2} \| \mathbf{u}^{(n+1)} \|_{L^2(\Omega)^d}^2 - \frac{1}{2} \| \mathbf{u}^{(n)} \|_{L^2(\Omega)^d}^2 + (\alpha - \frac{1}{2}) \| \mathbf{u}^{(n+1)} - \mathbf{u}^{(n)} \|_{L^2(\Omega)^d}^2 \\
& \quad + \partial_t C_{13} \| \mathbf{u}^{(n+\frac{1}{2})} \|_{\mathcal{E}}^2 \le \frac{\text{diam}(\Omega)^2}{2} \| \mathbf{f} \|_{(L^2(\Omega \times (t^{(n), t^{(n+1)})))^d)^2}^2 + \frac{\partial t}{2} \| \mathbf{u}^{(n+\frac{1}{2})} \|_{\mathcal{E}}^2.
\end{aligned}
$$

Summing for $n = 0, \dots, N - 1$ lead to both (41) and (42). Following the same ideas as in the proof of
Lemma 4.5, we get the existence of $C_{23}$, only depending on $\theta$, $\beta$ and $\Omega$ such that

$$
\begin{aligned}
\partial_t \| p^{(n+\frac{1}{2})} \|_{L^2(\Omega)} &\le C_{23} \left( \right. \\
&\qquad \partial_t \| u^{(n+\frac{1}{2})} \|_{\mathcal{E}} + \partial_t \| u^{(n+\frac{1}{2})} \|_{\mathcal{E}}^2 + \\
&\qquad \left. \| u^{(n+1)} - u^{(n)} \|_{L^2(\Omega)^d} + \sqrt{\partial_t} \| f \|_{(L^2(\Omega \times (t^{(n), t^{(n+1)})))^d)} \right).
\end{aligned}
$$

Summing the above inequality on $n = 0, \dots, N_T - 1$, using the Cauchy Schwarz inequality for the terms

$$
\partial_t \sum_{n=0}^{N_T-1} \| \mathbf{u}^{(n+\frac{1}{2})} \|_{\mathcal{E}} \text{ and } \sqrt{\partial_t} \sum_{n=0}^{N_T-1} \| \mathbf{f} \|_{(L^2(\Omega \times (t^{(n)}, t^{(n+1)}))^d}, \text{ inequality (41) for the term } \partial_t \sum_{n=0}^{N_T-1} \| \mathbf{u}^{(n+\frac{1}{2})} \|_{\mathcal{E}}^2
$$

$$
\text{and inequality (42) for the term } \sum_{n=0}^{N_T-1} \|\mathbf{u}^{(n+1)} - \mathbf{u}^{(n)}\|_{L^2(\Omega)^d} \text{ lead to (43).}
$$

The existence of at least one solution to the scheme is then again deduced from the use of the topological degree results, as in [8], as well as the uniqueness of the solution if $\beta = 0$. $\square$

Since our aim is to apply Theorem 7.1, we define a second norm on $\mathbf{H}_{\mathcal{E}}(\Omega)$ in the next definition.
**Definition 5.1 (Dual norm)** Let us define the following dual norm on $H_{\mathcal{E}}(\Omega)$:

$$ \begin{gathered} \forall \mathbf{w} \in \mathcal{H}_{\mathcal{E}}(\Omega), \\ \|\mathbf{w}\|_{\star, \mathcal{E}} = \sup \left\{ \int_{\Omega} \mathbf{w}(\mathbf{x}) \cdot \mathbf{v}(\mathbf{x}) \mathrm{d}\mathbf{x}, \mathbf{v} \in \mathcal{H}_{\mathcal{E}}(\Omega), \|\mathbf{v}\|_{\mathcal{E}} = 1 \text{ and } \mathrm{div}_{\mathcal{M}} \mathbf{v} = 0 \right\}. \end{gathered} \quad (44) $$

Let us denote by $Y_{\mathcal{E}}$ the space $H_{\mathcal{E}}(\Omega)$ equipped with the norm $\|\cdot\|_{\star, \mathcal{E}}$. We define a continuous embedding of $Y_{\mathcal{E}}$ in $E(\Omega)'$ (recall that $E(\Omega) = \{\mathbf{v} \in (H_0^1(\Omega))^2; \mathrm{div}\mathbf{v} = 0\}$) by the relation

$$ \langle \mathbf{w}, \mathbf{v} \rangle_{E(\Omega)', E(\Omega)} = \int_{\Omega} \mathbf{w}(\mathbf{x}) \cdot P_{\mathcal{E}} \mathbf{v}(\mathbf{x}) \mathrm{d}\mathbf{x}, \quad \forall \mathbf{w} \in \mathcal{H}_{\mathcal{E}}(\Omega), \forall \mathbf{v} \in E(\Omega). \quad (45) $$

where $P_{\mathcal{E}}\mathbf{v}$ is defined by (22) in Lemma. ???. Note that, thanks to (24) we get that

$$ \|\mathbf{w}\|_{E(\Omega)'} \le C_1 \|\mathbf{w}\|_{\star, \mathcal{E}}, \quad \forall \mathbf{w} \in \mathcal{H}_{\mathcal{E}}(\Omega), \quad (46) $$

where $C_1$ only depends on $\theta$ and $\Omega$.

In view of applying Theorem 7.1, let us study the dual norm of the discrete time derivative.

**Lemma 5.2 (Estimate on the dual norm of the discrete time derivative)** Let $\delta t = T/N_T$ with $N_T \in \mathbb{N}^*$, let $\mathcal{D}$ be a discretization as defined in Section 4.1 and let $\theta < \theta_{\mathcal{D}}$. Let $(\mathbf{u}^{(0)})$ be defined by (36), and, for a given $\alpha \in [\frac{1}{2}, 1]$, let $(\mathbf{u}^{(n+1)}, \mathbf{u}^{(n+\frac{1}{2})}, p^{(n+\frac{1}{2})}) \in \mathcal{H}_{\mathcal{E}}(\Omega) \times \mathcal{H}_{\mathcal{E}}(\Omega) \times H_{\mathcal{M}}(\Omega)$ such that (38)-(39) holds for $n = 1, \dots, N_T$. Then there exists $C_{24}$, only depending on $f, u_{\text{ini}}, T, \Omega, \beta$ and $\theta$, such that

$$ \|\delta_{\mathcal{D}, \delta t} \mathbf{u}\|_{L^{4/3}(0,T; Y_{\mathcal{E}})} \le C_{24}. \quad (47) $$

(Note that, thanks to (46), a similar inequality holds on $\|\delta_{\mathcal{D}, \delta t} \mathbf{u}\|_{L^{4/3}(0,T; E(\Omega)')}$.)

**Proof.** Let us take $\mathbf{v} \in \mathcal{H}_{\mathcal{E}}(\Omega)$, with $\mathrm{div}_{\mathcal{M}}\mathbf{v} = 0$, as test function in Scheme (39). Using (45) leads to

$$ \int_{\Omega} \delta_{\mathcal{D}, \delta t}^{(n+\frac{1}{2})} \mathbf{u}(\mathbf{x}) \cdot \mathbf{v}(\mathbf{x}) \mathrm{d}\mathbf{x} + \langle \mathbf{u}^{(n+\frac{1}{2})}, \mathbf{v} \rangle_{\mathcal{E}} + \beta b_{\mathcal{D}}(\mathbf{u}^{(n+\frac{1}{2})}, \mathbf{u}^{(n+\frac{1}{2})}, \mathbf{v}) = \frac{1}{\delta t} \int_{t^{(n)}}^{t^{(n+1)}} \int_{\Omega} f(\mathbf{x}, t) \cdot \mathbf{v}(\mathbf{x}) \mathrm{d}\mathbf{x} \mathrm{d}t. $$

Since $\mathbf{u}^{(n+\frac{1}{2})}$ satisfies (38) we get thanks to the estimate (28) of Lemma 4.4 and to the estimate (42) that:

$$ |b_{\mathcal{D}}(\mathbf{u}^{(n+\frac{1}{2})}, \mathbf{u}^{(n+\frac{1}{2})}, \mathbf{v})| \le C_{25} \|\mathbf{u}^{(n+\frac{1}{2})}\|_{\mathcal{E}}^{3/2} \|\mathbf{v}\|_{\mathcal{E}}, $$

where $C_{25}$, only depends on $f, u_{\text{ini}}, T, \Omega$ and $\theta$. We then get that there exists $C_{26}$, only depending on $f, u_{\text{ini}}, T, \beta, \Omega$ and $\theta$, such that

$$ \|\delta_{\mathcal{D}, \delta t} \mathbf{u}(t)\|_{\star, \mathcal{E}} \le C_{26} \left( \|\mathbf{u}^{(n+\frac{1}{2})}\|_{\mathcal{E}} + \|\mathbf{u}^{(n+\frac{1}{2})}\|_{\mathcal{E}}^{3/2} + \left\| \frac{1}{\delta t} \int_{t^{(n)}}^{t^{(n+1)}} f(\cdot, t) \mathrm{d}t \right\|_{L^2(\Omega)^d} \right). $$

Remarking that $\|\mathbf{u}^{(n+\frac{1}{2})}\|_{\mathcal{E}}^{\frac{4}{3}} \le 1 + \|\mathbf{u}^{(n+\frac{1}{2})}\|_{\mathcal{E}}^2$ and that

$$ \left\| \frac{1}{\delta t} \int_{t^{(n)}}^{t^{(n+1)}} \mathbf{f}(\cdot, t) \mathrm{d}t \right\|_{L^2(\Omega)^d}^{\frac{4}{3}} \le 1 + \left( \frac{1}{\delta t} \int_{t^{(n)}}^{t^{(n+1)}} \| \mathbf{f}(\cdot, t) \|_{L^2(\Omega)^d} \mathrm{d}t \right)^2 \\ \le 1 + \frac{1}{\delta t} \int_{t^{(n)}}^{t^{(n+1)}} \| \mathbf{f}(\cdot, t) \|_{L^2(\Omega)^d}^2 \mathrm{d}t, $$

we get that there exists $C_{27}$, only depending on $f, u_{\text{ini}}, T, \beta, \Omega$ and $\theta$, such that

$$ \|\delta_{\mathcal{D}, \delta t} \mathbf{u}(t)\|_{\star, \mathcal{E}}^{\frac{4}{3}} \le C_{27} \left( 1 + \|\mathbf{u}^{(n+\frac{1}{2})}\|_{\mathcal{E}}^2 + \frac{1}{\delta t} \int_{t^{(n)}}^{t^{(n+1)}} \|\mathbf{f}(\cdot, t)\|_{L^2(\Omega)^d}^2 \mathrm{d}t \right). $$

Integrating for $t \in (0, T)$ and using the discrete $L^2(H^1)$ estimate on the velocity (41), we get that there exists $C_{24}$, only depending on $f, u_{\text{ini}}, T, \Omega, \beta$ and $\theta$, such that (47) holds, which concludes the proof. $\square$
**Theorem 5.1 (Convergence of the scheme)** Let $\mathcal{F}$ be a family of time space discretizations $(\mathcal{D}, \partial t)$, where $\mathcal{D}$ is a discretization in the sense given in Section 4.1 and there exists $N_T \in \mathbb{N}^*$ such that $\partial t = T/N_T$. We assume that there exists $\theta > 0$ with $\theta \le \theta_{\mathcal{D}}$ for all $\mathcal{D} \in \mathcal{F}$ (see Section 4.1). For any $\mathcal{D} = (\mathcal{M}, \mathcal{E}, \partial t) \in \mathcal{F}$, let $(\mathbf{u}_{\mathcal{E},\partial t}, p_{\text{mesh},\partial t})$ denote a solution to Scheme (36)-(9)-(39) for a given $\alpha \in [\frac{1}{2}, 1]$. Then there exists a weak solution $\mathbf{u}$ of (4) such that, up to a subsequence, $\mathbf{u}_{\mathcal{E},\partial t}$ converges in $L^2(\Omega \times (0,T))^d$ to $\mathbf{u}$ as $h_{\mathcal{D}} \to 0$ and $\partial t \to 0$, which shows that

$$ \|\partial_t \mathbf{u}\|_{L^{4/3}(0,T;E(\Omega)^{\prime})} \le C_{24}. \quad (48) $$

Moreover, if $\beta = 0$, the whole family converges to the unique weak solution of (4) as $h_{\mathcal{D}} \to 0$ and $\partial t \to 0$.

**Proof.** Let $\mathcal{D}_{\ell} = (\mathcal{M}_{\ell}, \mathcal{E}_{\ell}, \partial t_{\ell})_{\ell \in \mathbb{N}}$ be a sequence of elements of $\mathcal{F}$, such that $h_{\mathcal{D}_{\ell}}$ and $\partial t_{\ell}$ tend to 0 as $\ell \to \infty$. From thereon, we denote $\mathbf{u}_{\mathcal{M}_{\ell},\partial t_{\ell}}$ by $\mathbf{u}_{\ell}$ for short.

**Step 1 Proof that hypotheses (h1-h4) of Theorem 7.1 hold, and consequences.**

In our setting, the space $B$ of Theorem 7.1 is $L^2(\Omega)^d$. We take $B_{\ell} = \{\mathbf{w} \in \mathbf{H}_{\mathcal{E}_{\ell}}; \text{div}_{\mathcal{D}_{\ell}} \mathbf{w} = 0\}$. The norm $\| \cdot \|_{X_{\ell}}$ is the norm $\| \cdot \|_{\mathcal{E}_{\ell}}$, and the norm $\| \cdot \|_{Y_{\ell}}$ is defined in Definition 5.1.

Let $(\mathbf{w}_{\ell})_{\ell \in \mathbb{N}}$ be a sequence of functions of $B_{\ell}$ such that $\|\mathbf{w}_{\ell}\|_{\mathcal{E}_{\ell}} \le C$ for some $C \in \mathbb{R}_{+}$. Then by definition of the norm $\|\mathbf{w}_{\ell}\|_{\mathcal{E}_{\ell}}$ and thanks to Rellich's theorem, the sequence $(\hat{\mathbf{w}}_{\ell})_{\mathcal{E}_{\ell}}$ converges in $L^2(\Omega)^d$ to some $\mathbf{w} \in L^2(\Omega)^d$; therefore, by the inequality (52) of 7.1 given in the Appendix, the sequence $(\mathbf{w}_{\ell})_{\mathcal{E}_{\ell}}$ also converges to $\mathbf{w}$ in $L^2(\Omega)^d$. Thus, assumption (h1) of Theorem 7.1 is satisfied.

Let us then show that assumption (h2) is also satisfied. Let $(\mathbf{w}_{\ell})_{\ell \in \mathbb{N}}$ be a sequence of functions of $L^2(\Omega)^d$ such that $\mathbf{w}_{\ell} \in B_{\ell}$ and $\|\mathbf{w}_{\ell}\|_{\mathcal{E}_{\ell}} \le C$ for some $C \in \mathbb{R}_{+}$, and such that there exists $\mathbf{w} \in B$ with $\mathbf{w}_{\ell} \to \mathbf{w}$ in $B$ and $\|\mathbf{w}_{\ell}\|_{\star, \mathcal{E}_{\ell}} \to 0$ as $\ell \to \infty$. By definition of the norm $\| \cdot \|_{\star, \mathcal{E}_{\ell}}$, we have

$$ \int_{\Omega} \mathbf{w}_{\ell}(\mathbf{x}) \cdot \mathbf{w}_{\ell}(\mathbf{x}) \mathrm{d}\mathbf{x} \le \|\mathbf{w}_{\ell}\|_{\star, \mathcal{E}_{\ell}} \|\mathbf{w}_{\ell}\|_{\mathcal{E}_{\ell}} \to 0 \text{ as } \ell \to \infty, $$

which shows that $\mathbf{w} = 0$.

From estimates (41) and (47), we get that hypotheses (h3) and (h4) of Theorem 7.1 are satisfied. Therefore, we deduce that there exists $\mathbf{u} \in L^1(0, T; L^2(\Omega)^d)$ and a subsequence of $(\mathcal{D}_{\ell}, \partial t_{\ell})_{\ell \in \mathbb{N}}$, again denoted by $(\mathcal{D}_{\ell}, \partial t_{\ell})_{\ell \in \mathbb{N}}$, such that $\mathbf{u}_{\ell}$ converges in $L^1(0, T; L^2(\Omega)^d)$ to $\mathbf{u}$ as $\ell \to \infty$.

**Step 2 Proof that $\mathbf{u}$ is a weak solution of (4).**

As in [8], we easily get that $\mathbf{u} \in L^2(0, T; H_0^1(\Omega)^d)$ and we prove the convergence in $L^2(\Omega \times (0, T))^d$ to $\mathbf{u}$ of $\hat{\mathbf{u}}_{\ell}$. Lemma 7.1 implies the same convergence property for $\mathbf{u}_{\ell}$. We get that $\text{div}\mathbf{u}(\mathbf{x}, t) = 0$ in the same manner as in the steady case. The proof that $\mathbf{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 (48) holds.**

Thanks to (47), extracting a subsequence such that $\delta_{\mathcal{D}_{\ell}, \partial t_{\ell}} \mathbf{u}_{\ell}$ weakly converges in the Sobolev space $L^{4/3}(0, T; E(\Omega)^{\prime})$, we may pass to the limit in (45). We thus get

$$ \int_0^T \varphi(t) \langle \lim_{\ell \to \infty} \delta_{\mathcal{D}_{\ell}, \partial t_{\ell}} \mathbf{u}_{\ell}(t), \mathbf{v} \rangle_{E(\Omega)^{\prime}, E(\Omega)} \mathrm{d}t = - \int_0^T \varphi'(t) \int_{\Omega} \mathbf{u}(\mathbf{x}, t) \mathbf{v}(\mathbf{x}) \mathrm{d}\mathbf{x}\mathrm{d}t, \\ \forall \mathbf{v} \in E(\Omega), \forall \varphi \in C_c^1((0, T)). $$

The preceding relation allows to conclude that

$$ \partial_t \mathbf{u} = \lim_{\ell \to \infty} \delta_{\mathcal{D}_{\ell}, \partial t_{\ell}} \mathbf{u}_{\ell}(t), $$

and therefore that (48) holds.

□
# 6 Numerical examples

## 6.1 An analytical steady problem

We consider a problem where the continuous solution of the Navier-Stokes equations (1) in the case (2) with $\beta = 1$ is given by:

$$ \bar{u}_1(x_1, x_2) = 2\pi \sin^2(\pi x_1) \cos(\pi x_2) \sin(\pi x_2) $$

$$ \bar{u}_2(x_1, x_2) = -2\pi \cos(\pi x_1) \sin(\pi x_1) \sin^2(\pi x_2) $$

$$ \bar{p}(x_1, x_2) = \sin^2(\pi x_1) \sin^2(\pi x_2) $$

in a circle with centre $(0.5, 0.5)$ and radius $0.45$. We consider four meshes for the mass conservation $\mathcal{M}_j$, $j = 0, \dots, 3$, defined in the following way:

1. a structured square $10 \times 10$ is given on the square $[0, 1] \times [0, 1]$,

2. for $i = 0, \dots, 3$, let us split in 4 control volumes each grid block whose centre $(x_1, x_2)$ satisfies

$$ \sqrt{(x_1 - 0.5)^2 + (x_2 - 0.5)^2} \geq 0.45 - 0.25/2^i, $$

3. get rid of all the control volumes $K$ with centre $(x_1, x_2)$ such that

$$ \sqrt{(x_1 - 0.5)^2 + (x_2 - 0.5)^2} > 0.45. $$

Let us denote $\text{card}(\mathcal{M}_j)$ the number of control volumes of the mesh $\mathcal{M}_j$. We get that $\text{card}(\mathcal{M}_0) = 1604$, $\text{card}(\mathcal{M}_1) = 6416$, $\text{card}(\mathcal{M}_2) = 25592$ and $\text{card}(\mathcal{M}_3) = 102324$. The $L^2$ errors of unknowns $u_1, u_2, p$, respectively denoted by $e_2(u_1), e_2(u_2), e_2(p)$, are respectively computed in the Voronoï grids associated to the velocity components and in $\mathcal{M}_j$.

Left part of Figure (7) shows the errors $\log_{10}(e_2(u_1))$ and $\log_{10}(e_2(p))$ with respect to $\log_{10}(1/\sqrt{\text{card}(\mathcal{M}_j)})$ for $j = 0, \dots, 3$. On right part of Figure 7 are plotted the stream lines for the finest mesh. The velocity

<figure><img src="image_9.png" /><figcaption>Figure 7: Left: The $L^2$ error with respect to the number of control volumes. Right: Stream lines.</figcaption></figure>

components and the pressure are shown in Figure (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.

## 6.2 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 Figure 9) in the case (2). We consider the case where the Reynolds number is equal to 1000, 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 Figure 10. The results are compared to the literature [5] in Table 1. The results show the expected precision.
<figure><img src="image_10.png" /><figcaption>Figure 8: Horizontal component of the velocity (left), vertical component of the velocity (middle), pressure (right) for $j = 0$ (top) and $j = 2$ (bottom).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: An example of mesh, and the two media axes $\xi_1$ and $\xi_2$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Velocity components, $D_2^m u^{(1)}$ (left) $u^{(2)}$ (right)</figcaption></figure>

<table><caption>Table 1: Maximum and minimum of the velocity components along the centerlines $\xi_1$ and $\xi_2$.</caption><thead><tr><th rowspan="2">#M</th><th colspan="4">Generalized MAC scheme</th><th>[5]</th></tr><tr><th>924</th><th>3698</th><th>14796</th><th>59190</th><th>102400 (= 320<sup>2</sup>)</th></tr></thead><tbody><tr><td>min (u<sup>(1)</sup>(ξ<sub>2</sub>))</td><td>-1.860 · 10<sup>-1</sup></td><td>-1.987 · 10<sup>-1</sup></td><td>-1.994 · 10<sup>-1</sup></td><td>-1.984 · 10<sup>-1</sup></td><td>-1.9772 · 10<sup>-1</sup></td></tr><tr><td>ξ<sub>2</sub></td><td>7.500 · 10<sup>-1</sup></td><td>7.826 · 10<sup>-1</sup></td><td>7.875 · 10<sup>-1</sup></td><td>7.8125 · 10<sup>-1</sup></td><td>7.8185 · 10<sup>-1</sup></td></tr><tr><td>min (u<sup>(2)</sup>(ξ<sub>1</sub>))</td><td>-2.143 · 10<sup>-2</sup></td><td>-2.251 · 10<sup>-2</sup></td><td>-2.065 · 10<sup>-2</sup></td><td>-2.010 · 10<sup>-2</sup></td><td>-1.9979 · 10<sup>-2</sup></td></tr><tr><td>ξ<sub>1</sub></td><td>3.423 · 10<sup>-1</sup></td><td>3.238 · 10<sup>-1</sup></td><td>3.146 · 10<sup>-1</sup></td><td>3.192 · 10<sup>-1</sup></td><td>3.1751 · 10<sup>-1</sup></td></tr><tr><td>max (u<sup>(2)</sup>(ξ<sub>1</sub>))</td><td>1.361 · 10<sup>-2</sup></td><td>1.361 · 10<sup>-2</sup></td><td>1.254 · 10<sup>-2</sup></td><td>1.216 · 10<sup>-2</sup></td><td>1.2123 · 10<sup>-2</sup></td></tr><tr><td>ξ<sub>1</sub></td><td>7.853 · 10<sup>-1</sup></td><td>8.038 · 10<sup>-1</sup></td><td>8.269 · 10<sup>-1</sup></td><td>8.246 · 10<sup>-1</sup></td><td>8.2656 · 10<sup>-1</sup></td></tr></tbody></table>

## 6.3 A time-dependent case

We consider the Green-Taylor case with $\beta = 10$, in the same domain $\Omega$ as the one considered in Section 6.1, which is a circle, with center $(0.5, 0.5)$ and radius $0.45$, on the time domain $[0, 0.1]$. Nonhomogeneous
Dirichlet boundary conditions are imposed. The solution is given, in the case (2), by

$$
\begin{align*}
u^{(1)}(x_1, x_2, t) &= -\cos(\pi x_1) \sin(\pi x_2) \exp(-2\pi^2 t), \\
u^{(2)}(x_1, x_2, t) &= \sin(\pi x_1) \cos(\pi x_2) \exp(-2\pi^2 t), \\
p(x_1, x_2, t) &= -\beta \frac{\cos(2\pi x_1) + \cos(2\pi x_2)}{4} \exp(-4\pi^2 t).
\end{align*}
$$

We plot in Figure 11 the $L^2$ errors evaluated at time $t = 0.1$ as a function of the space step ($-0.5 \times \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 Tab. 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 $\alpha = 1/2$, which does not occur in the case $\alpha = 1$. If we plot the values $\frac{1}{2}(p^{(n-\frac{1}{2})} + p^{(n+\frac{1}{2})})$ instead of $p^{(n+\frac{1}{2})}$, we observe that this post-processed pressure numerically converges. All the converging curves show an order 2. In order to check the behaviour of the pressure in the case $\alpha = 1/2$, we show, in Figures 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 $\alpha = 1/2$ may lead to oscillations in pressures, it is simple to get back converging pressures.

<table><caption>Table 2: Time and space steps used for spatial order of convergence</caption>
  <tr>
    <td>δt</td>
    <td>6.25 10<sup>-3</sup></td>
    <td>1.5625 10<sup>-3</sup></td>
    <td>3.90625 10<sup>-4</sup></td>
    <td>9.765625 10<sup>-5</sup></td>
  </tr>
  <tr>
    <td>#M</td>
    <td>1604</td>
    <td>6416</td>
    <td>25592</td>
    <td>102324</td>
  </tr>
</table>

<figure><img src="image_13.png" /><figcaption>Figure 11: $L^2$ errors (+: velocity ×: pressure) as a function of the size of the mesh.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12: An example of pressure solution at two consecutive time steps. $p^{(n-\frac{1}{2})}$ (left) $\frac{1}{2}(p^{(n-\frac{1}{2})} + p^{(n+\frac{1}{2})})$ (middle) and $p^{(n+\frac{1}{2})}$ with $\#M = 1604$ and $\delta t = 6.25 \times 10^{-3}$</figcaption></figure>

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 Figure 14 (for the case $\alpha = 1/2$, we have compared the exact pressure with the average pressure $\frac{1}{2}(p^{(n-\frac{1}{2})} + p^{(n+\frac{1}{2})})$ instead of $p^{(n+\frac{1}{2})}$). This shows an order 2 for $\alpha = 1/2$, and 1 for $\alpha = 1$.
<figure><img src="image_15.png" /><figcaption>Figure 13: An example of pressure solution at two consecutive time steps. $p^{(n-\frac{1}{2})}$ (left) $\frac{1}{2}(p^{(n-\frac{1}{2})} + p^{(n+\frac{1}{2})})$ (middle) and $p^{(n+\frac{1}{2})}$ with $\#M = 25592$ and $\delta t = 3.90625 \ 10^{-4}$</figcaption></figure>

<table><caption>Table 3: Time and space steps used for time order of convergence</caption><tr><td>δt</td><td>0.05</td><td>0.025</td><td>0.0125</td><td>0.00635</td><td>0.003125</td></tr><tr><td>#M</td><td>1604</td><td>6416</td><td>25592</td><td>102324</td><td>409156</td></tr></table>

<figure><img src="image_16.png" /><figcaption>Figure 14: $L^2$ errors (+: velocity ×: pressure) as a function of the time step.</figcaption></figure>

# 7 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.

# Appendix

## Interpolation results

The first lemma that we give is a technical result which allows to bound, under adequate geometrical conditions, the difference between various interpolates by the gradient of a $P1$ reconstruction.

**Lemma 7.1 (Comparison between interpolates)** Let $\Omega \subset \mathbb{R}^d$, with $d \in \mathbb{N}^*$ be an open polygonal connected open set, such that there exists a $H^1$ conforming simplicial mesh $\mathcal{T}$ over $\Omega$. We denote by $\mathcal{T}$ the elements of $\mathcal{M}$, by $\mathcal{V}$ the set of all vertices, and by $\mathcal{V}_T$ the set of all vertices of $\mathcal{T}$ and by $\mathcal{T}_T \subset \mathcal{T}$ the set of the simplices sharing a face with $\mathcal{T}$. For $s \in \mathcal{V}$, we denote by $\xi_s$ the $P^1$ finite element basis function associated to the node $x_s$. We assume that are given some non-negative functions $\psi_s$, for all $s \in \mathcal{V}$, such that

$$ \sum_{s \in \mathcal{V}} \psi_s(x) = 1, \quad \text{for a.e. } x \in \Omega. \tag{49} $$
We denote by $\mathcal{W}$ the set of all pairs of vertices $(s, s') \in \mathcal{V}^2$ such that $\int_{\Omega} \psi_{s'}(\mathbf{x})\xi_s(\mathbf{x})d\mathbf{x} > 0$, and we assume that a set $\mathcal{T}(s, s') \subset \mathcal{T}$ is defined such that it contains all the elements $T \in \mathcal{T}$ with $s \in \mathcal{V}_T$, it contains at least one $T \in \mathcal{T}$ such that $s' \in \mathcal{V}_T$ and for any pair $(T, T') \in \mathcal{T}(s, s')$, there exists a sequence of simplices $T_1, \dots, T_M$ such that $T_{i+1} \in \mathcal{T}_{T_i}$ for $i = 1, \dots, M-1$, $T = T_1$ and $T' = T_M$.

Let $C_{28} > 1$ be defined by

$$C_{28} = \max\left\{\frac{\mathrm{diam}(T)^d}{|T|}, T \in \mathcal{T}\right\} \cup \left\{\frac{\mathrm{diam}(T')}{\mathrm{diam}(T)}, T \in \mathcal{T}, T' \in \mathcal{T}_T\right\} \cup \left\{\#\mathcal{T}(s, s'), (s, s') \in \mathcal{W}\right\}. \quad (50)$$

For $(u_s)_{s \in \mathcal{V}} \in \mathbb{R}^\mathcal{V}$, let us denote $\hat{u} = \sum_{s \in \mathcal{V}} u_s \xi_s$ and $u = \sum_{s \in \mathcal{V}} u_s \psi_s$, and, for any $T \in \mathcal{T}$ and $\mathbf{x} \in T$, let us denote $\delta(\mathbf{x}) = \mathrm{diam}(T)$. Then there exists $C_{29} > 0$ only depending on $C_{28}$, such that

$$\left\| \frac{1}{\delta}(u - \hat{u}) \right\|_{L^2(\Omega)^d} \le C_{29} \left\| \nabla \hat{u} \right\|_{L^2(\Omega)^d}, \quad (51)$$

so that in particular,

$$\|u - \hat{u}\|_{L^2(\Omega)^d} \le C_{29} \max_{T \in \mathcal{T}} \mathrm{diam}(T) \|\nabla \hat{u}\|_{L^2(\Omega)^d}. \quad (52)$$

**Proof.** By definition of $\hat{u}$ and $u$, thanks to the fact that $\sum_{s \in \mathcal{V}} \xi_s(\mathbf{x}) = \sum_{s \in \mathcal{V}} \psi_s(\mathbf{x}) = 1$, and applying the Cauchy–Schwarz inequality, we get:

$$\begin{aligned} \| \frac{1}{\delta}(u - \hat{u}) \|_{L^2(\Omega)^d}^2 &= \int_{\Omega} \frac{1}{\delta(\mathbf{x})^2} \left( \sum_{s \in \mathcal{V}} u_s \xi_s(\mathbf{x}) - \sum_{s' \in \mathcal{V}} u_{s'} \psi_{s'}(\mathbf{x}) \right)^2 d\mathbf{x} \\ &= \int_{\Omega} \frac{1}{\delta(\mathbf{x})^2} \left( \sum_{s \in \mathcal{V}} \sum_{s' \in \mathcal{V}} (u_s - u_{s'}) \xi_s(\mathbf{x}) \psi_{s'}(\mathbf{x}) \right)^2 d\mathbf{x}, \\ &\le \int_{\Omega} \frac{1}{\delta(\mathbf{x})^2} \sum_{(s, s') \in \mathcal{W}} (u_s - u_{s'})^2 \xi_s(\mathbf{x}) \psi_{s'}(\mathbf{x}) d\mathbf{x}. \end{aligned}$$

We remark that the closure of $\mathcal{T}(s, s')$ contains $s, s'$ and an interior connected path from $s$ to $s'$, whose intersection with any $T \in \mathcal{T}(s, s')$ is either a nonzero vector $w(T, s, s') \in \mathbb{R}^d$ or empty (then we set $w(T, s, s') = 0$, and therefore satisfies

$$\sum_{T \in \mathcal{T}(s, s')} w(T, s, s') = s' - s.$$

We may then write, since $\hat{u}$ is affine on each $T \in \mathcal{T}$,

$$u_s - u_{s'} = \sum_{T \in \mathcal{T}(s, s')} \nabla_T \hat{u} \dot{w}(T, s, s'),$$

denoting by $\nabla_T \hat{u} = \sum_{s \in \mathcal{V}_T} u_s \nabla \xi_s(\mathbf{x})$, for all $\mathbf{x} \in T$, the constant gradient of $\hat{u}$ in the simplex $T$. We get

$$(u_s - u_{s'})^2 \le \sum_{T \in \mathcal{T}(s, s')} |w(T, s, s')|^2 \sum_{T \in \mathcal{T}(s, s')} |\nabla_T \hat{u}|^2,$$

which provides, denoting by $(L(s, s'))^2 = \sum_{T \in \mathcal{T}(s, s')} \mathrm{diam}(T)^2$,

$$(u_s - u_{s'})^2 \le (L(s, s'))^2 \sum_{T \in \mathcal{T}(s, s')} |\nabla_T \hat{u}|^2,$$

and therefore

$$\begin{aligned} \| \frac{1}{\delta}(u - \hat{u}) \|_{L^2(\Omega)^d}^2 &\le \int_{\Omega} \sum_{(s, s') \in \mathcal{W}} \frac{(L(s, s'))^2}{\delta(\mathbf{x})^2} \sum_{T \in \mathcal{T}(s, s')} |\nabla_T \hat{u}|^2 \xi_s(\mathbf{x}) \psi_{s'}(\mathbf{x}) d\mathbf{x} \\ &= \sum_{T \in \mathcal{T}} |\nabla_T \hat{u}|^2 \int_{\Omega} \sum_{(s, s') \in \mathcal{W}_T} \frac{(L(s, s'))^2}{\delta(\mathbf{x})^2} \xi_s(\mathbf{x}) \psi_{s'}(\mathbf{x}) d\mathbf{x}, \end{aligned}$$
denoting by $\mathcal{W}_T$ the subset of all $(s, s') \in \mathcal{W}$ such that $T \in \mathcal{T}(s, s')$. Let us show that there exists $C_{29} > 0$ only depending on $C_{28}$, such that, for any $T \in \mathcal{T}$,

$$ \int_{\Omega} \sum_{(s,s') \in \mathcal{W}_T} \frac{(L(s, s'))^2}{\delta(\mathbf{x})^2} \xi_s(\mathbf{x}) \psi_{s'}(\mathbf{x}) \mathrm{d}\mathbf{x} \leq C_{29}^2 |T|. \quad (53) $$

We choose a point $\mathbf{x}_T \in T$. Let $(s, s') \in \mathcal{W}_T$, let $\mathbf{x} \in \Omega$ be such that $\xi_s(\mathbf{x})\psi_{s'}(\mathbf{x}) > 0$. Since $\xi_s(\mathbf{x}) > 0$, there exists $T_1 \in \mathcal{T}$ such that $\mathbf{x} \in T_1$ and $s \in \mathcal{V}_{T_1}$. Let us observe that, for any $T' \in \mathcal{T}(s, s')$, then $\mathrm{diam}(T') \leq (C_{28})^M \mathrm{diam}(T_1)$, which implies that

$$ (L(s, s'))^2 \leq C_{28} (C_{28})^{2M} \mathrm{diam}(T_1)^2 = (C_{28})^{2M+1} \delta(\mathbf{x})^2. \quad (54) $$

Since $T_1 \in \mathcal{T}(s, s')$, there exists a family of neighbouring simplices $T_1, \dots, T_M$, with $M \leq C_{28}$ (using (50)), such that $T = T_M$. Therefore, $\mathrm{diam}(T_{M-1}) \leq C_{28}\mathrm{diam}(T)$, and, by induction, $\mathrm{diam}(T_1) \leq (C_{28})^{M-1}\mathrm{diam}(T)$. This leads to the inequality $|\mathbf{x}-\mathbf{x}_T| \leq C_{30}\mathrm{diam}(T)$, denoting by $C_{30} = (C_{28})^M/(C_{28}-1)$. Therefore, we get, using $\sum_{(s,s') \in \mathcal{W}_T} \xi_s(\mathbf{x})\psi_{s'}(\mathbf{x}) \leq 1$, that

$$ \int_{\Omega} \sum_{(s,s') \in \mathcal{W}_T} \xi_s(\mathbf{x}) \psi_{s'}(\mathbf{x}) \mathrm{d}\mathbf{x} \leq |B(\mathbf{x}_T, C_{30}\mathrm{diam}(T))|. $$

Denoting by $C_{31}$ the real, only depending on the space dimension (2 or 3), such that $|B(\mathbf{x}_T, \mathrm{diam}(T))| = C_{31}\mathrm{diam}(T)^d$, we may write

$$ |B(\mathbf{x}_T, C_{30}\mathrm{diam}(T))| = C_{31}\mathrm{diam}(T)^d C_{30}^d \leq C_{31}C_{28}C_{30}^d |T|. $$

Using (54), this concludes the proof of (53).

Thanks to $\sum_{T \in \mathcal{T}} |T| |\nabla_T \hat{u}|^2 = \|\sum_{s \in \mathcal{V}} u_s \nabla \xi_s\|_{L^2(\Omega)^d}$, we then conclude (51).

□

**Lemma 7.2 (Comparison between $L^p$ norms)** Let $\Omega \subset \mathbb{R}^d$, with $d \in \mathbb{N}^*$ be an open connected set, such that there exists a $H^1$ conforming simplicial mesh $\mathcal{T}$ over $\Omega$. We denote by $T$ the elements of $\mathcal{M}$, by $\mathcal{V}$ the set of all vertices. For $s \in \mathcal{V}$, we denote by $\xi_s$ the $P^1$ finite element basis function associated to the node $s$. Then, for any $p \in [1, +\infty)$, there exists $C(d, p) > 0$, such that

$$ \| \sum_{s \in \mathcal{V}} u_s \xi_s \|_{L^p(\Omega)^d} \le \left( \sum_{s \in \mathcal{V}} |u_s|^p \int_{\Omega} \xi_s(\mathbf{x}) \mathrm{d}\mathbf{x} \right)^{1/p} \le C(d,p) \| \sum_{s \in \mathcal{V}} u_s \xi_s \|_{L^p(\Omega)^d}, \quad \forall (u_s)_{s \in \mathcal{V}} \in \mathbb{R}^V. \quad (55) $$

**Proof.** The left inequality of (55) is an immediate consequence of the convexity of the function $x \to |x|^p$, which provides

$$ \left| \sum_{s \in \mathcal{V}} u_s \xi_s(\mathbf{x}) \right|^p \leq \sum_{s \in \mathcal{V}} |u_s|^p \xi_s(\mathbf{x}). $$

Let us turn to the proof of the right inequality of (55). Let $T \in \mathcal{T}$ and $(a_1, \dots, a_{d+1}) \in \mathbb{R}^{d+1}$ be given, and let, for any $\mathbf{x} \in T$, $\xi_1(\mathbf{x}), \dots, \xi_{d+1}(\mathbf{x})$ be the $d+1$ barycentric coordinates of $\mathbf{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 $\varepsilon \in (0, 1/2)$, which will be chosen later, by

$$ T_i(\varepsilon) = \{ \mathbf{x} \in T, \xi_i(\mathbf{x}) \geq 1 - \varepsilon \}. $$

Since, on $T_i$, we have, for $j \neq i$, $\xi_j(\mathbf{x}) \leq \varepsilon < 1 - \varepsilon$, we get that all the $T_i$, $i = 1, \dots, d+1$, are disjoint. We then have

$$ \int_T \left| \sum_{i=1}^{d+1} a_i \xi_i(\mathbf{x}) \right|^p \mathrm{d}\mathbf{x} \geq \sum_{i=1}^{d+1} \int_{T_i(\varepsilon)} \left| \sum_{i=1}^{d+1} a_i \xi_i(\mathbf{x}) \right|^p \mathrm{d}\mathbf{x}. $$

We may write, for $\mathbf{x} \in T_{d+1}(\varepsilon)$,

$$ \begin{aligned} |a_{d+1}|(1-\varepsilon) & \leq |a_{d+1}\xi_{d+1}(\mathbf{x})| \leq \left| \sum_{i=1}^{d+1} a_i \xi_i(\mathbf{x}) \right| + \sum_{i=1}^{d} |a_i \xi_i(\mathbf{x})| \\ & \leq \left| \sum_{i=1}^{d+1} a_i \xi_i(\mathbf{x}) \right| + \varepsilon \sum_{i=1}^{d} |a_i|, \end{aligned} $$
and therefore, thanks to the Young inequality

$$|a_{d+1}|^p (1 - \varepsilon)^p \le (d+1)^{p-1} \left( \left| \sum_{i=1}^{d+1} a_i \xi_i(\mathbf{x}) \right|^p + \varepsilon^p \sum_{i=1}^d |a_i|^p \right),$$

Let us remark that the measure of $T_i(\varepsilon)$ is equal to $\varepsilon^d|T|$. We then get

$$\int_T \left| \sum_{i=1}^{d+1} a_i \xi_i(\mathbf{x}) \right|^p d\mathbf{x} \ge \varepsilon^d |T| \sum_{i=1}^{d+1} \left( \frac{|a_i|^p (1-\varepsilon)^p}{(d+1)^{p-1}} - \varepsilon^p \sum_{j \ne i} |a_j|^p \right),$$

which provides

$$\int_T \left| \sum_{i=1}^{d+1} a_i \xi_i(\mathbf{x}) \right|^p d\mathbf{x} \ge \varepsilon^d |T| \left( \frac{(1-\varepsilon)^p}{(d+1)^{p-1}} - d\varepsilon^p \right) \sum_{i=1}^{d+1} |a_i|^p.$$

We then choose $\varepsilon$ such that

$$\varepsilon = \min\left(\frac{1}{3}, \frac{1}{(3d(d+1)^{p-1})^{1/p}}\right).$$

We then have

$$\frac{(1 - \varepsilon)^p}{(d + 1)^{p-1}} \ge \frac{2^p}{3^p (d + 1)^{p-1}} \ge \frac{2}{3(d + 1)^{p-1}},$$

and therefore

$$\frac{(1 - \varepsilon)^p}{(d + 1)^{p-1}} - d\varepsilon^p \ge \frac{1}{3(d + 1)^{p-1}}.$$

Denoting $C(d, p) > 0$ the quantity defined by

$$\frac{C(d, p)^p}{d+1} = \left( \min\left(\frac{1}{3}, \frac{1}{(3d(d+1)^{p-1})^{1/p}}\right) \right)^d \frac{1}{3(d+1)^{p-1}},$$

we have then proved that

$$\int_T \left| \sum_{i=1}^{d+1} a_i \xi_i(\mathbf{x}) \right|^p d\mathbf{x} \ge C(d, p)^p \frac{|T|}{d+1} \sum_{i=1}^{d+1} |a_i|^p,$$

which implies (55).

□

## Discrete Aubin-Simon result

As we mentioned earlier, the compactness result that we use in our proof of convergence of the scheme is an adaptation of a previous result [10]. The differences between the theorem given below and [10, Theorem 3.4] are the following:

1. Theorem 7.1 is a compactness result in $L^1(0, T; B)$ where $B$ is some Banach space, while [10, Theorem 3.4] is a compactness result in $L^p(0, T; B)$ for $p \in [1, +\infty)$. We chose to restrict to $L^1$ for two reasons:

* The proof 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 (41) and (42), using the “$L^p - L^q$ compactness” property and Sobolev inequalities, see [9, Lemma 3.1].

2. The assumption (h2) that we require here is weaker than the assumption (h2) of [10, Theorem 3.4], which does not require that $(\|w_n\|_{X_n})_{\ell \in \mathbb{N}}$ be bounded. A quick look at the proofs of [10, 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 (56) of $u_{\ell}(\cdot, t)$ corresponds to an $\alpha$-scheme whereas that of [10, Theorem 3.4] corresponds to an implicit scheme, i.e. $\alpha$ is equal to 1.

**Theorem 7.1 (Discrete Aubin-Simon lemma)** Let $T > 0$ and let $B$ be a Banach space. Let $(B_{\ell})_{\ell \in \mathbb{N}}$ be a sequence of finite dimensional subspaces of $B$. For any $\ell \in \mathbb{N}$, let $N_{\ell} \in \mathbb{N}^*$, let $t_{\ell}^{(0)} = 0 < t_{\ell}^{(1)} < \dots < t_{\ell}^{(N_{\ell})} = T$ and $\delta t_{\ell}^{(n)} = t_{\ell}^{(n)} - t_{\ell}^{(n-1)}$, $n = 1, \dots, N_{\ell}$. Let $\{u_{\ell}^{(n)}, n = 0, \dots, N_{\ell}\} \subset B_{\ell}$ and let $u_{\ell} \in L^1(0, T; B_{\ell})$ be defined, for a given real family $(\alpha_{\ell}^{(n)})_{n=1, \dots, N_{\ell}}$, by

$$ u_{\ell}(\cdot, t) = \tilde{u}_{\ell}^{(n)} := (1 - \alpha_{\ell}^{(n)})u_{\ell}^{(n-1)} + \alpha_{\ell}^{(n)}u_{\ell}^{(n)} \in B_{\ell}, \quad (56) \\ \text{for a.e. } t \in (t_{\ell}^{(n-1)}, t_{\ell}^{(n)}), \text{ and } n \in \{1, \dots, N_{\ell}\}. $$

Let $\delta_{\ell}u_{\ell}$ be the "discrete time derivative", defined by:

$$ \delta_{\ell}u_{\ell}(\cdot, t) = \delta_{\ell}^{(n)}u_{\ell} := \frac{1}{\delta_{\ell}^{(n)}}(u_{\ell}^{(n)} - u_{\ell}^{(n-1)}) \text{ for a.e. } t \in (t_{\ell}^{(n-1)}, t_{\ell}^{(n)}), n \in \{1, \dots, N_{\ell}\}. $$

Let $\| \cdot \|_{X_{\ell}}$ and $\| \cdot \|_{Y_{\ell}}$ be two norms on $B_{\ell}$. We denote by $X_{\ell}$ the space $B_{\ell}$ endowed with the norm $\| \cdot \|_{X_{\ell}}$ and by $Y_{\ell}$ the space $B_{\ell}$ endowed with the norm $\| \cdot \|_{Y_{\ell}}$. We assume that

(h1) For any sequence $(w_{\ell})_{\ell \in \mathbb{N}}$ such that $w_{\ell} \in B_{\ell}$ and $(\|w_{\ell}\|_{X_{\ell}})_{\ell \in \mathbb{N}}$ is bounded, then, up to a subsequence, there exists $w \in B$ such that $w_{\ell} \to w$ in $B$ as $\ell \to +\infty$.

(h2) For any sequence $(w_{\ell})_{\ell \in \mathbb{N}}$ such that $w_{\ell} \in B_{\ell}$, $(\|w_{\ell}\|_{X_{\ell}})_{\ell \in \mathbb{N}}$ is bounded, there exists $w \in B$ such that $w_{\ell} \to w$ in $B$ and $\|w_{\ell}\|_{Y_{\ell}} \to 0$ as $\ell \to +\infty$, then $w = 0$.

(h3) The family $(\alpha_{\ell}^{(n)})_{n=1, \dots, N_{\ell}, \ell \in \mathbb{N}}$ and the sequence $(\|u_{\ell}\|_{L^1(0,T;X_{\ell})})_{\ell \in \mathbb{N}}$ are bounded.

(h4) The sequence $(\|\delta_{\ell}u_{\ell}\|_{L^1(0,T;Y_{\ell})})_{\ell \in \mathbb{N}}$ is bounded.

Then there exists $u \in L^1(0, T; B)$ such that, up to a subsequence, $u_{\ell} \to u$ in $L^1(0, T; B)$ as $\ell \to +\infty$.

**Proof.** The proof is mainly inspired by that of [10, Theorem 3.4]. The first step is to apply [10, Lemma 3.2], which states that, under Hypothesis (h1), there exists $C_X > 0$ such that, for all $\ell \in \mathbb{N}$, for all $u \in B_{\ell}$, $\|u\|_{B} \le \|u\|_{X_{\ell}}$. We then remark that, thanks to (h1) and (h2), for all $\varepsilon > 0$, there exists $C(\varepsilon) > 0$, such that,

$$ \forall \ell \in \mathbb{N}, \forall u \in B_{\ell}, \quad \|u\|_{B} \le \varepsilon\|u\|_{X_{\ell}} + C(\varepsilon)\|u\|_{Y_{\ell}}, \quad (57) $$

following the proof of [10, Lemma 3.1]. We then remark that, using (h4), there exists $C_t > 0$ such that

$$ \forall \ell \in \mathbb{N}, \sum_{n=1}^{N_{\ell}} \delta_{\ell}^{(n)} \|\delta_{\ell}^{(n)} u_{\ell}\|_{Y_{\ell}} \le C_t. $$

Therefore we get from (h3) that the sequence $(\|u_{\ell}\|_{BV(0,T;Y_{\ell})})_{\ell \in \mathbb{N}}$ is bounded; indeed we can write

$$ \sum_{n=1}^{N_{\ell}-1} \|\tilde{u}_{\ell}^{(n+1)} - \tilde{u}_{\ell}^{(n)}\|_{Y_{\ell}} \le \sum_{n=1}^{N_{\ell}-1} \left( \alpha_{\ell}^{(n+1)} \delta_{\ell}^{(n+1)} \|\delta_{\ell}^{(n+1)} u_{\ell}\|_{Y_{\ell}} + (1 + \alpha_{\ell}^{(n)}) \delta_{\ell}^{(n)} \|\delta_{\ell}^{(n)} u_{\ell}\|_{Y_{\ell}} \right) \\ \le (1 + 2C_{\alpha})C_t, $$

denoting by $C_{\alpha}$ a bound of the family $(\alpha_{\ell}^{(n)})_{n=1, \dots, N_{\ell}, \ell \in \mathbb{N}}$. We then prolong $u_{\ell}$ by symmetry on $(-T, 2T)$, setting $u_{\ell}(\cdot, -t) = u_{\ell}(\cdot, t)$ and $u_{\ell}(\cdot, T + t) = u_{\ell}(\cdot, T - t)$ for a.e. $t \in (0, T)$. We then get that $u_{\ell} \in BV(-T, 2T; Y_{\ell})$, with $\|u_{\ell}\|_{BV(-T, 2T; Y_{\ell})} \le 3(1 + 2C_{\alpha})C_t$. We also get, using (h3), that the sequence $(\|u_{\ell}\|_{L^1(-T, 2T; X_{\ell})})_{\ell \in \mathbb{N}}$ is bounded as well, which leads that there exists $C_0$ independent on $\ell$ such that

$$ \int_{-T}^{2T} \|u_{\ell}(\cdot, t)\|_{X_{\ell}} dt \le C_0. $$

We have, thanks to (57), for any $\tau \in (0, T)$ and a.e. $t \in (-T, 2T - \tau)$,

$$ \|u_{\ell}(\cdot, t + \tau) - u_{\ell}(\cdot, t)\|_{B} \le \varepsilon \|u_{\ell}(\cdot, t + \tau) - u_{\ell}(\cdot, t)\|_{X_{\ell}} + C(\varepsilon) \|u_{\ell}(\cdot, t + \tau) - u_{\ell}(\cdot, t)\|_{Y_{\ell}}, $$
which provides

$$ \int_{-T}^{2T-\tau} \|u_{\ell}(\cdot, t+\tau) - u_{\ell}(\cdot, t)\|_{B} dt \le \varepsilon 2C_{0} + \tau C(\varepsilon)3(1+2C_{\alpha})C_{t}. $$

This proves that $\int_{-T}^{2T-\tau} \|u_{\ell}(\cdot, t+\tau) - u_{\ell}(\cdot, t)\|_{B} dt$ tends to 0 with $\tau$, uniformly with respect to $\ell \in \mathbb{N}$. Thanks to the multiplication of $u_{\ell}$ by a function $\psi \in C_c^{\infty}(-T, 2T)$, equal to 1 on $(0, T)$, we may apply the compactness theorem [10, Theorem 2.1] to obtain that the family $(u_{\ell})_{\ell \in \mathbb{N}}$ is relatively compact in $L^1(0, T; B)$.

<figure><img src="image_17.png" /></figure>

**Remark 7.1 (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 [22]; indeed, the hypothesis “B continuously embedded in Y”, which is equivalent to the following hypothesis: If $w_n \to w$ in B and $\|w_n\|_Y \to 0$, then $w = 0$, may be replaced by the weaker assumption “if $(w_n)_{n \in \mathbb{N}}$ is bounded in X, converge in B to w and $\|w_n\|_Y \to 0$, then $w = 0$”.

## References

[1] P. Blanc. Convergence of a finite volume scheme on a MAC mesh for the Stokes problem with right-hand side in $H^{-1}$. In *Finite volumes for complex applications IV*, pages 133–142. ISTE, London, 2005.

[2] J. H. Bramble. A proof of the inf-sup condition for the Stokes equations on Lipschitz domains. *Math. Models Methods Appl. Sci.*, 13(3):361–371, 2003. Dedicated to Jim Douglas, Jr. on the occasion of his 75th birthday.

[3] R. Chénier, R. Eymard, and R. Herbin. An extension of the mac scheme to some unstructured meshes. In *Finite volumes for complex applications VI*, volume 1, pages 253–261. Springer, London, 2011. Finite Volumes for Complex Applications VI (FVCA VI), Prague, Czech Republic, June 2011.

[4] S. H. Chou and D. Y. Kwak. A covolume method based on rotated bilinears for the generalized Stokes problem. *SIAM J. Numer. Anal.*, 35(2):494–507 (electronic), 1998.

[5] I. Demirdzic, Z. Lilek, and M. Peric. Fluid flow and heat transfer test problems for non-orthogonal grids: Bench-mark solutions. *International Journal for Numerical Methods in Fluids*, 15:329–354, 1992.

[6] R. Eymard, J. Fuhrmann, and A. Linke. MAC schemes on triangular meshes. In *Finite volumes for complex applications VI*, volume 1, pages 399–407. Springer, London, 2011. Finite Volumes for Complex Applications VI (FVCA VI), Prague, Czech Republic, June 2011.

[7] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In P. G. Ciarlet and J.-L. Lions, editors, *Techniques of Scientific Computing, Part III*, Handbook of Numerical Analysis, VII, pages 713–1020. North-Holland, Amsterdam, 2000.

[8] R. Eymard, R. Herbin, and J.-C. Latché. Convergence analysis of a colocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3d meshes. *SIAM J. Numer. Anal.*, 45(1):1–36, 2007.

[9] P. Fabrie and T. Gallouët. Modeling wells in porous media flow. *Math. Models Methods Appl. Sci.*, 10(5):673–709, 2000.

[10] T. Gallouët and J.-C. Latché. Compactness of discrete approximate solutions to parabolic pdes - application to a turbulence model. to appear, 2012.

[11] V. Girault and H. Lopez. Finite-element error estimates for the MAC scheme. *IMA J. Numer. Anal.*, 16(3):247–379, 1996.

[12] H. Han and X. Wu. A new mixed finite element formulation and the MAC method for the Stokes equations. *SIAM J. Numer. Anal.*, 35(2):560–571 (electronic), 1998.

[13] F. Harlow and J. Welch. Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface. *Physics of Fluids*, 8:2182–2189, 1965.
[14] G. Kanschat. Divergence-free discontinuous Galerkin schemes for the Stokes equations and the MAC scheme. *Internat. J. Numer. Methods Fluids*, 56(7):941–950, 2008.

[15] J. Nečas. *Les méthodes directes en théorie des équations elliptiques*. Masson et Cie, Éditeurs, Paris, 1967.

[16] R. Nicolaïdes. Analysis and convergence of the mac scheme i: The linear problem. *SIAM J. Numer. Anal.*, 29:1579–1591, 1992.

[17] R. Nicolaïdes and X. Wu. Analysis and convergence of the mac scheme ii, Navier-Stokes equations. *Math. Comp.*, 65:29–44, 1996.

[18] R. A. Nicolaides. The covolume approach to computing incompressible flows. In *Incompressible computational fluid dynamics: trends and advances*, pages 295–333. Cambridge Univ. Press, Cambridge, 2008.

[19] S. Patankar. *Numerical heat transfer and fluid flow. Series in Computational Methods in Mechanics and Thermal Sciences*, volume XIII. Washington - New York - London: Hemisphere Publishing Corporation; New York. McGraw-Hill Book Company, 1980.

[20] T. A. Porsching. Error estimates for MAC-like approximations to the linear Navier-Stokes equations. *Numer. Math.*, 29(3):291–306, 1977/78.

[21] L. R. Scott and S. Zhang. Finite element interpolation of nonsmooth functions satisfying boundary conditions. *Math. Comp.*, 54(190):483–493, 1990.

[22] J. Simon. Compact sets in the space $lp(0,t;b)$. *Annali di Matematica Pura ed Applicata*, 146:65–96, 1987.

[23] D. Vidović, A. Segal, and P. Wesseling. A superlinearly convergent Mach-uniform finite volume method for the Euler equations on staggered unstructured grids. *J. Comput. Phys.*, 217(2):277–294, 2006.

[24] P. Wesseling. *Principles of Computational Fluid Dynamics*. Springer, 2001.