Convergence of the MAC scheme for the incompressible Navier-Stokes equations

TL;DR

This paper proves the convergence of the MAC scheme for incompressible Navier-Stokes equations on non-uniform grids, establishing existence, compactness, and that the limit solutions are weak solutions of the continuous problem.

Contribution

It provides a rigorous convergence proof of the MAC scheme without regularity assumptions on solutions, extending previous results to non-uniform grids.

Findings

A priori estimates on scheme solutions

Existence of discrete solutions

Convergence to weak solutions of Navier-Stokes

Abstract

We prove in this paper the convergence of the Marker and cell (MAC) scheme for the dis-cretization of the steady-state and unsteady-state incompressible Navier-Stokes equations in primitive variables on non-uniform Cartesian grids, without any regularity assumption on the solution. A priori estimates on solutions to the scheme are proven ; they yield the existence of discrete solutions and the compactness of sequences of solutions obtained with family of meshes the space step of which tends to zero. We then establish that the limit is a weak solution to the continuous problem.

Equations

$$\mathrm{div}\bar{\mathbf{u}}=0$$

$$-\Delta \bar{\mathbf{u}}+(\bar{\mathbf{u}}\cdot \nabla )\bar{\mathbf{u}}+\nabla \bar{p}=\bar{\mathbf{f}}$$

$$\bar{\mathbf{u}}=0$$

$$\text{mbox}Find(\bar{\mathbf{u}},\bar{p})\in H_0^1(\Omega {)}^d\times L_0^2(\Omega )\text{mbox}suchthat,\forall (\mathbf{v},q)\in H_0^1(\Omega {)}^d\times L_0^2(\Omega ),$$

$${\int }_{\Omega }\nabla \bar{\mathbf{u}}:\nabla \mathbf{v}\mathrm{d}\mathbf{x}+{\int }_{\Omega }((\bar{\mathbf{u}}\cdot \nabla )\bar{\mathbf{u}})\cdot \mathbf{v}\mathrm{d}\mathbf{x}-{\int }_{\Omega }\bar{p}\mathrm{div}\mathbf{v}\mathrm{d}\mathbf{x}={\int }_{\Omega }\bar{\mathbf{f}}\cdot \mathbf{v}\mathrm{d}\mathbf{x},$$

$${\int }_{\Omega }q\mathrm{div}\bar{\mathbf{u}}\mathrm{d}\mathbf{x}=0,$$

$$\mathrm{div}\bar{\mathbf{u}}=0$$

$${\partial }_t\bar{\mathbf{u}}-\Delta \bar{\mathbf{u}}+(\bar{\mathbf{u}}\cdot \nabla )\bar{\mathbf{u}}+\nabla \bar{p}=\bar{\mathbf{f}}$$

$$\bar{\mathbf{u}}=0$$

$$\bar{\mathbf{u}}(\mathbf{x},0)={\mathbf{u}}_0$$

$${\boldsymbol{E}}(\Omega)=\bigl{\{}{\boldsymbol{u}}\in H^{1}_{0}(\Omega)^{d}~{};\ {\mathrm{div}}{\boldsymbol{u}}=0\mbox{ a.e. in }\Omega\bigr{\}}.$$

$$\begin{array}[]{l}\displaystyle\mbox{Find }{\boldsymbol{u}}\in L^{2}(0,T;{\boldsymbol{E}}(\Omega))\cap L^{\infty}(0,T;L^{2}(\Omega)^{d})\mbox{ such that, }\forall{\boldsymbol{v}}\in L^{2}(0,T;{\boldsymbol{E}}(\Omega))\cap C_{c}^{\infty}(\Omega\times[0,T)),\\[4.30554pt] \displaystyle\hskip 30.1388pt-\int_{0}^{T}\int_{\Omega}\bar{\boldsymbol{u}}(\boldsymbol{x},t)\cdot\partial_{t}{\boldsymbol{v}}(\boldsymbol{x},t)\ \mathrm{d}\boldsymbol{x}\ \mathrm{d}t-\int_{\Omega}{{\boldsymbol{u}}_{0}}(\boldsymbol{x})\cdot{\boldsymbol{v}}(\boldsymbol{x},0)\ \mathrm{d}\boldsymbol{x}+\int_{0}^{T}\int_{\Omega}{\boldsymbol{\nabla}}\bar{\boldsymbol{u}}(\boldsymbol{x},t):{\boldsymbol{\nabla}}{\boldsymbol{v}}(\boldsymbol{x},t)\ \mathrm{d}\boldsymbol{x}\ \mathrm{d}t\\[4.30554pt] \hfill\displaystyle+\int_{0}^{T}\int_{\Omega}((\bar{\boldsymbol{u}}\cdot{\boldsymbol{\nabla}})\bar{\boldsymbol{u}})(\boldsymbol{x},t)\cdot{\boldsymbol{v}}(\boldsymbol{x},t)\ \mathrm{d}\boldsymbol{x}\ \mathrm{d}t=\int_{0}^{T}\int_{\Omega}\bar{\boldsymbol{f}}(\boldsymbol{x},t)\cdot{\boldsymbol{v}}(\boldsymbol{x},t)\ \mathrm{d}\boldsymbol{x}\ \mathrm{d}t.\end{array}$$

$$\Omega=\cup_{{\sigma}\in{{\mathcal{E}}^{(i)}}}D_{\sigma},\quad i\in\bigl{[}|1,d|\bigr{]},$$

$$d_{ϵ}=\{\begin{array}{ll}d({\mathbf{x}}_{\sigma },{\mathbf{x}}_{{\sigma }^{\prime }})& \text{mbox}ifϵ=\sigma |{\sigma }^{\prime }\in {\tilde{\mathcal{E}}}_{\mathrm{int}}^{(i)},\\ d({\mathbf{x}}_{\sigma },ϵ)& \text{mbox}ifϵ\in {\tilde{\mathcal{E}}}_{\mathrm{ext}}^{(i)},\end{array}$$

$$\displaystyle h_{\mathcal{M}}=\max\bigl{\{}{\mathrm{diam}}(K),K\in{\mathcal{M}}\bigr{\}},$$

$$\displaystyle\eta_{\mathcal{M}}=\max\Bigl{\{}\frac{|{\sigma}|}{|{\sigma}^{\prime}|},\ {\sigma}\in{{\mathcal{E}}^{(i)}},\ {\sigma}^{\prime}\in{{\mathcal{E}}^{(j)}},\ i,j\in\bigl{[}|1,d|\bigr{]},\ i\not=j\Bigr{\}},$$

$$H_{{\mathcal{E}}^{(i)},0}=\Bigl{\{}u\in H_{{\mathcal{E}}^{(i)}},\ u(\boldsymbol{x})=0\ \forall\boldsymbol{x}\in D_{\sigma},\ {\sigma}\in{\mathcal{E}}_{\mathrm{ext}}^{(i)}\Bigr{\}}.$$

$$u_{i}=\sum_{{\sigma}\in{{\mathcal{E}}^{(i)}}}u_{\sigma}{1\!\!1}_{D_{\sigma}},\ i\in\bigl{[}|1,d|\bigr{]}\quad\mbox{and}\quad p=\sum_{K\in{\mathcal{M}}}p_{K}{1\!\!1}_{K}.$$

$$\begin{array}[]{l|l}-\Delta_{{\mathcal{E}}^{(i)}}:&\quad H_{{\mathcal{E}}^{(i)},0}\longrightarrow H_{{\mathcal{E}}^{(i)},0}\\ &\displaystyle\quad u_{i}\longmapsto-\Delta_{{\mathcal{E}}^{(i)}}u_{i}=-\sum_{{\sigma}\in{{\mathcal{E}}^{(i)}}}(\Delta u)_{\sigma}{1\!\!1}_{D_{\sigma}},\mbox{ with }-(\Delta u)_{\sigma}=\frac{1}{|D_{\sigma}|}\sum_{\epsilon\in{\widetilde{{\mathcal{E}}}}(D_{\sigma})}\phi_{{\sigma},\epsilon}(u_{i}),\\ &\displaystyle\quad\phi_{{\sigma},\epsilon}(u_{i})=\begin{cases}\displaystyle\frac{|\epsilon|}{d_{\epsilon}}(u_{\sigma}-u_{{\sigma}^{\prime}}),\text{ if }\ \epsilon={\sigma}|{\sigma}^{\prime}\in{{\widetilde{{\mathcal{E}}}}^{(i)}_{{\rm int}}},\\[8.61108pt] \displaystyle\frac{|\epsilon|}{d_{\epsilon}}u_{\sigma},\text{ if }\epsilon\in{{\widetilde{{\mathcal{E}}}}^{(i)}_{{\rm ext}}}\cap{\widetilde{{\mathcal{E}}}}(D_{\sigma}),\end{cases}\end{array}$$

$${\varphi }_{\sigma ,ϵ}(u_i)=-{\varphi }_{{\sigma }^{\prime },ϵ}(u_i),\forall ϵ=\sigma |{\sigma }^{\prime }\in {\tilde{\mathcal{E}}}_{\mathrm{int}}^{(i)}.$$

$$\begin{array}[]{l|l}-\Delta_{\mathcal{E}}:&{\boldsymbol{H}}_{{\mathcal{E}},0}\longrightarrow{\boldsymbol{H}}_{{\mathcal{E}},0}\\[4.30554pt] &{\boldsymbol{u}}\mapsto-\Delta_{\mathcal{E}}{\boldsymbol{u}}=(-\Delta_{{\mathcal{E}}^{(1)}}u_{1},\ldots,-\Delta_{{\mathcal{E}}^{(d)}}u_{d}).\end{array}$$

$$\forall({\boldsymbol{u}},{\boldsymbol{v}})\in{\boldsymbol{H}}_{{\mathcal{E}},0}\times{\boldsymbol{H}}_{{\mathcal{E}},0},\ \int_{\Omega}-\Delta_{\mathcal{E}}{\boldsymbol{u}}\cdot{\boldsymbol{v}}\ \mathrm{d}\boldsymbol{x}=[{\boldsymbol{u}},{\boldsymbol{v}}]_{1,{\mathcal{E}},0}=\sum_{i=1}^{d}[u_{i},v_{i}]_{1,{{\mathcal{E}}^{(i)}},0},\\ \mbox{with, for }i\in\bigl{[}|1,d|\bigr{]},\ [u_{i},v_{i}]_{1,{{\mathcal{E}}^{(i)}},0}=\sum_{\begin{subarray}{c}\epsilon\in{{\widetilde{{\mathcal{E}}}}^{(i)}_{{\rm int}}}\\ \epsilon={\sigma}|{\sigma}^{\prime}\end{subarray}}\frac{|\epsilon|}{d_{\epsilon}}\ (u_{\sigma}-u_{{\sigma}^{\prime}})\ (v_{\sigma}-v_{{\sigma}^{\prime}})+\sum_{\begin{subarray}{c}\epsilon\in{{\widetilde{{\mathcal{E}}}}^{(i)}_{{\rm ext}}}\\ \epsilon\in{\widetilde{{\mathcal{E}}}}(D_{\sigma})\end{subarray}}\frac{|\epsilon|}{d_{\epsilon}}\ u_{\sigma}\ v_{\sigma}.$$

$$\left|\begin{array}[]{l}H_{{\mathcal{E}}^{(i)},0}\times H_{{\mathcal{E}}^{(i)},0}\to\xR\\[4.30554pt] (u_{i},v_{i})\mapsto[u_{i},v_{i}]_{1,{{\mathcal{E}}^{(i)}},0}\end{array}\right.\quad\mbox{and}\quad\left|\begin{array}[]{l}{\boldsymbol{H}}_{{\mathcal{E}},0}\times{\boldsymbol{H}}_{{\mathcal{E}},0}\to\xR\\[4.30554pt] ({\boldsymbol{u}},{\boldsymbol{v}})\mapsto[{\boldsymbol{u}},{\boldsymbol{v}}]_{1,{\mathcal{E}},0}\end{array}\right.$$

$$\displaystyle\|u_{i}\|^{2}_{1,{{\mathcal{E}}^{(i)}},0}=[u_{i},u_{i}]_{1,{{\mathcal{E}}^{(i)}},0}=\sum_{\begin{subarray}{c}\epsilon\in{{\widetilde{{\mathcal{E}}}}^{(i)}_{{\rm int}}}\\ \epsilon={\sigma}|{\sigma}^{\prime}\end{subarray}}\frac{|\epsilon|}{d_{\epsilon}}\ (u_{\sigma}-u_{{\sigma}^{\prime}})^{2}+\sum_{\begin{subarray}{c}\epsilon\in{{\widetilde{{\mathcal{E}}}}^{(i)}_{{\rm ext}}}\\ \epsilon\in{\widetilde{{\mathcal{E}}}}(D_{\sigma})\end{subarray}}\frac{|\epsilon|}{d_{\epsilon}}\ u_{\sigma}^{2},\quad\mbox{for }i\in\bigl{[}|1,d|\bigr{]},$$

$$\Vert \mathbf{u}{\Vert }_{1,\mathcal{E},0}^2=[\mathbf{u},\mathbf{u}{]}_{1,\mathcal{E},0}=\sum _{i=1}^d\Vert u_i{\Vert }_{1,{\mathcal{E}}^{(i)},0}^2.$$

$$(D_{\epsilon})_{\epsilon\in{{\widetilde{{\mathcal{E}}}}^{(i)}},\epsilon\perp{\boldsymbol{e}}^{(j)}},\mbox{ with }\left|\begin{array}[]{ll}D_{\epsilon}=\epsilon\times[\boldsymbol{x}_{\sigma}\,\boldsymbol{x}_{{\sigma}^{\prime}}]&\mbox{ if }\epsilon\mbox{ lies inside }\Omega,\ \epsilon={\sigma}|{\sigma}^{\prime},\\[4.30554pt] D_{\epsilon}=\epsilon\times[\boldsymbol{x}_{\sigma}\,\boldsymbol{x}_{{\sigma},\epsilon}]&\mbox{ if }\epsilon\mbox{ lies on }\partial\Omega,\ \epsilon\in{\widetilde{{\mathcal{E}}}}(D_{\sigma}),\end{array}\right.$$

$$\begin{array}[]{lll}\mbox{ if }\epsilon\mbox{ lies inside }\Omega,\ \epsilon=\overrightarrow{{\sigma}|{\sigma}^{\prime}},&(\eth_{j}u_{i})_{D_{\epsilon}}=\dfrac{u_{{\sigma}^{\prime}}-u_{\sigma}}{d_{\epsilon}},\\ \mbox{ if }\epsilon\mbox{ lies on }\partial\Omega,\ \epsilon\in{\widetilde{{\mathcal{E}}}}(D_{\sigma}),&(\eth_{j}u_{i})_{D_{\epsilon}}=\dfrac{-u_{\sigma}}{d_{\epsilon}}\ \overrightarrow{\boldsymbol{x}_{\sigma}\boldsymbol{x}_{{\sigma},\epsilon}}\cdot{\boldsymbol{e}}^{(j)},\end{array}$$

$${\nabla }_{{\tilde{\mathcal{E}}}^{(i)}}{u}_i=({ð}_1{u}_i,\dots ,{ð}_d{u}_i)\text{mbox}with{ð}_j{u}_i=\sum _{\begin{array}{c}ϵ\in {\tilde{\mathcal{E}}}^{(i)}\\ ϵ\perp {\mathbf{e}}^{(j)}\end{array}}({ð}_j{u}_i{)}_{D_{ϵ}}{11}_{D_{ϵ}}.$$

$$\int_{\Omega}{\boldsymbol{\nabla}}_{{\widetilde{{\mathcal{E}}}}^{(i)}}\psi\cdot{\boldsymbol{\nabla}}_{{\widetilde{{\mathcal{E}}}}^{(i)}}\chi\ \mathrm{d}\boldsymbol{x}=[\psi,\chi]_{1,{{\mathcal{E}}^{(i)}},0},\mbox{ for }\psi,\chi\in H_{{\mathcal{E}}^{(i)},0},\mbox{ and }i\in\bigl{[}|1,d|\bigr{]}.$$

$${\nabla }_{\tilde{\mathcal{E}}}\mathbf{u}=({\nabla }_{{\tilde{\mathcal{E}}}^{(1)}}{u}_1,\dots ,{\nabla }_{{\tilde{\mathcal{E}}}^{(d)}}{u}_d),$$

$${\int }_{\Omega }{\nabla }_{\tilde{\mathcal{E}}}\mathbf{u}:{\nabla }_{\tilde{\mathcal{E}}}\mathbf{v}\mathrm{d}\mathbf{x}=[\mathbf{u},\mathbf{v}{]}_{1,\mathcal{E},0}.$$

$$\Pi_{\mathcal{E}}{\boldsymbol{\varphi}}=(\Pi_{{\mathcal{E}}^{(1)}}\varphi_{1},\cdots,\Pi_{{\mathcal{E}}^{(d)}}\varphi_{d})\in H_{{\mathcal{E}}^{(1)},0}\times\cdots\times H_{{\mathcal{E}}^{(d)},0},\mbox{ where }\\ \bigl{|}(\Pi_{{\mathcal{E}}^{(i)}}\varphi_{i})_{\sigma}-\varphi_{i}(\boldsymbol{x}_{\sigma})\bigr{|}\leq C_{{\boldsymbol{\varphi}}}\ h_{\mathcal{M}}^{2},\mbox{ for }{\sigma}\in{{\mathcal{E}}^{(i)}},\ i\in\bigl{[}|1,d|\bigr{]}.$$

$$|\eth_{j}\Pi_{{\mathcal{E}}^{(i)}}\varphi_{i}(\boldsymbol{x})-\partial_{j}\varphi_{i}(\boldsymbol{x})|\leq C_{{\boldsymbol{\varphi}},\eta}\ h_{\mathcal{M}}\mbox{ for a.e. }\boldsymbol{x}\in\Omega\mbox{ and for }i,j\in\bigl{[}|1,d|\bigr{]}.$$

$$\displaystyle\begin{array}[]{l| l}\displaystyle{\mathrm{div}}_{\mathcal{M}}:&\quad{\boldsymbol{H}}_{{\mathcal{E}},0}\longrightarrow L_{{\mathcal{M}},0}\\[4.30554pt] &\displaystyle\quad{\boldsymbol{u}}\longmapsto{\mathrm{div}}_{\mathcal{M}}\,{\boldsymbol{u}}=\sum_{K\in{\mathcal{M}}}\frac{1}{|K|}\sum_{{\sigma}\in{\mathcal{E}}(K)}|{\sigma}|\ u_{K,{\sigma}}\ {1\!\!1}_{K},\end{array}$$

$$\displaystyle\mbox{with }u_{K,{\sigma}}=u_{\sigma}\boldsymbol{n}_{K,{\sigma}}\cdot{\boldsymbol{e}}^{(i)}\mbox{ for }{\sigma}\in{{\mathcal{E}}^{(i)}}\cap{\mathcal{E}}(K),\ i\in\bigl{[}|1,d|\bigr{]}.$$

Full text

Keywords Finite-volume methods, MAC scheme, incompressible Navier-Stokes.

Communicated by Douglas N. Arnold.

1. 1 Introduction

2. 2 Space discretization

3. 3 The steady case

4. 3.1 The scheme

5. 3.2 Stability and existence of a solution

6. 3.3 Convergence analysis

7. 4 The time-dependent case

8. 4.1 Time discretization

9. 4.2 Estimates on discrete solutions and existence

10. 4.3 Convergence analysis

11. 4.4 Case of the unsteady Stokes equations

12. 5 Appendix: Discrete functional analysis

1. Introduction

Let $\Omega$ be an open bounded domain of $\xR^{d}$ with $d=2$ or $d=3$. The steady-state incompressible Navier-Stokes equations read:

[TABLE]

where $\bar{\boldsymbol{u}}$ stands for the (vector-valued) velocity of the flow, $\bar{p}$ for the pressure and $\bar{\boldsymbol{f}}$ is a given field of $L^{2}(\Omega)^{d}$, and where, for two given vector fields ${\boldsymbol{v}}=(v_{1},\ldots,v_{d})$ and ${\boldsymbol{w}}=(w_{1},\ldots,w_{d})$, the quantity $({\boldsymbol{v}}\cdot{\boldsymbol{\nabla}}){\boldsymbol{w}}$ is a vector field whose components are $(({\boldsymbol{v}}\cdot{\boldsymbol{\nabla}}){\boldsymbol{w}})_{i}=\sum_{k=1}^{d}v_{k}\partial_{k}w_{i}$, i\in\bigl{[}|1,d|\bigr{]}. A weak formulation of Problem (1) reads:

[TABLE]

where $L^{2}_{0}(\Omega)$ stands for the subspace of $L^{2}(\Omega)$ of zero mean-valued functions.

The time-dependent Navier-Stokes equations are also considered:

[TABLE]

This problem is posed for $(\boldsymbol{x},t)$ in $\Omega\times(0,T)$ where $T\in\xR_{+}^{*}$; the right-hand side $\bar{\boldsymbol{f}}$ is now a given vector field of $L^{2}(\Omega\times(0,T))^{d}$ and the initial datum ${\boldsymbol{u}}_{0}$ belongs to the space ${\boldsymbol{E}}(\Omega)$ of divergence-free functions, defined by:

[TABLE]

A weak formulation of the transient problem (3) reads (see e.g. [3]):

[TABLE]

The Marker-And-Cell (MAC) scheme, introduced in the middle of the sixties [21], is one of the most popular methods [28, 33] for the approximation of the Navier-Stokes equations in the engineering framework, because of its simplicity, its efficiency and its remarkable mathematical properties. The aim of this paper is to show, under minimal regularity assumptions on the solution, that sequences of approximate solutions obtained by the discretization of problem (1)(resp. (3)) by the MAC scheme converge to a solution of (2)(resp. (4)) as the mesh size (resp. the mesh size and the time step) tends (resp. tend) to 0.

For the linear problems, the first error analysis seems to be that of [29] in the case of the time-dependent Stokes equations on uniform square grids. The mathematical analysis of the scheme was performed for the steady-state Stokes equations in [26] for uniform rectangular meshes with $H^{2}$-regularity assumption on the pressure. Error estimates for the MAC scheme applied to the Stokes equations have been obtained by viewing the MAC scheme as a mixed finite element method [19, 20] or a divergence conforming DG method [22]. Error estimates for rectangular meshes were also obtained for the related covolume method, see [6] and references therein. Using the tools that were developed for the finite volume theory [11, 12], an order 1 error estimate for non-uniform meshes was obtained in [1], with order 2 convergence for uniform meshes, under the usual regularity assumptions ($H^{2}$ for the velocities, $H^{1}$ for the pressure). It was recently shown in [24] that under higher regularity assumptions ($C^{4}$ for the velocities and $C^{3}$ for the pressure) and an additional convergence assumption on the pressure, superconvergence is obtained for non uniform meshes. Note also that the convergence of the MAC scheme for the Stokes equations with a right-hand side in $H^{-1}(\Omega)$ was proven in [2].

Mathematical studies of the MAC scheme for the nonlinear Navier-Stokes equations are scarcer. A pioneering work was that of [27] for the steady-state Navier-Stokes equations and for uniform rectangular grids. More recently, a variant of the MAC scheme was defined on locally refined grids and the convergence proof was performed for both the steady-state and time dependent cases in two or three space dimensions [4]. A MAC-like scheme was also studied for the stationary Stokes and Navier-Stokes equations on two-dimensional Delaunay-Voronoï grids [9]. For the Stokes equations on uniform grids, the scheme given in [4] coincides with the usual MAC scheme that is classically used in CFD codes. However, for the Navier-Stokes equations, the nonlinear convection term is discretized in [4] and [9] in a manner reminiscent of what is sometimes done in the finite element framework (see e.g. [32]), which no longer coincides with the usual MAC scheme, even on uniform rectantular grids; this discretization entails a larger stencil, and numerical experiments [5] seem to show that it is not as efficient as the classical MAC scheme.

Our purpose here is to analyse the genuine MAC scheme for the steady-state and transient Navier-Stokes equations in primitive variables on a non-uniform rectangular mesh in two or three dimensions, and, as in [4], without any assumption on the data nor on the regularity of the the solutions. The convergence of a subsequence of approximate solutions to a weak solution of the Navier-Stokes equations is proved for both the steady and unsteady case, which yields as a by product the existence of a weak solution, well known since the work of J. Leray [23]. In the case where uniqueness of the solution is known, the whole sequence of approximate solutions can be shown to converge, see remarks 3.14 and 4.4.

This paper is organized as follows. In Section 2, the MAC space grid and the discrete operators are introduced. In particular, the velocity convection operator is approximated so as to be compatible with a discrete continuity equation on the dual cells ; this discretization coincides with the usual discretization on uniform meshes [28], contrary to the scheme of [4]. The MAC scheme for the steady state Navier-Stokes equations and its weak formulation are introduced in Section 3. Velocity and pressure estimates are then obtained, which lead to the compactness of sequences of approximate solutions. Any prospective limit is shown to be a weak solution of the continuous problem. Section 4 is devoted to the unsteady Navier-Stokes equations. An essential feature of the studied scheme is that the (discrete) kinetic energy remains controlled. We show the compactness of approximate sequences of solutions thanks to a discrete Aubin-Simon argument, and again conclude that any limit of the approximate velocities is a weak solution of the Navier-Stokes equations, thanks to a passage to the limit in the scheme. In the case of the unsteady Stokes equations, some additional estimates yield the compactness of sequences of approximate pressures; this entails that the approximate pressure converges to a weak solution of the Stokes equations as the mesh size and time steps tend to 0.

2. Space discretization

Let $\Omega$ be a connected subset of $\xR^{d}$ consisting in a union of rectangles ($d=2$) or orthogonal parallelepipeds ($d=3$); without loss of generality, the edges (or faces) of these rectangles (or parallelepipeds) are assumed to be orthogonal to the canonical basis vectors, denoted by $({\boldsymbol{e}}^{(1)},\ldots,{\boldsymbol{e}}^{(d)})$.

Definition 2.1** (MAC grid).**

A discretization of $\Omega$ with a MAC grid, denoted by $\mathcal{D}$, is defined by $\mathcal{D}=({\mathcal{M}},{\mathcal{E}})$, where:

• –

${\mathcal{M}}$ stands for the primal grid, and consists in a conforming structured partition of $\Omega$ in possibly non uniform rectangles ($d=2$) or rectangular parallelepipeds ($d=3$). A generic cell of this grid is denoted by $K$, and its mass center by $\boldsymbol{x}_{K}$. The pressure is associated to this mesh, and ${\mathcal{M}}$ is also sometimes referred to as ”the pressure mesh”.

• –

The set of all faces of the mesh is denoted by ${\mathcal{E}}$; we have ${\mathcal{E}}={\mathcal{E}}_{\mathrm{int}}\cup{\mathcal{E}}_{\mathrm{ext}}$, where ${\mathcal{E}}_{\mathrm{int}}$ (resp. ${\mathcal{E}}_{\mathrm{ext}}$) are the edges of ${\mathcal{E}}$ that lie in the interior (resp. on the boundary) of the domain. The set of faces that are orthogonal to ${\boldsymbol{e}}^{(i)}$ is denoted by ${{\mathcal{E}}^{(i)}}$, for i\in\bigl{[}|1,d|\bigr{]}. We then have ${{\mathcal{E}}^{(i)}}={\mathcal{E}}_{\mathrm{int}}^{(i)}\cup{\mathcal{E}}_{\mathrm{ext}}^{(i)}$, where ${\mathcal{E}}_{\mathrm{int}}^{(i)}$ (resp. ${\mathcal{E}}_{\mathrm{ext}}^{(i)}$) are the edges of ${{\mathcal{E}}^{(i)}}$ that lie in the interior (resp. on the boundary) of the domain.

For ${\sigma}\in{\mathcal{E}}_{\mathrm{int}}$, we write ${\sigma}=K|L$ if ${\sigma}=\partial K\cap\partial L$. A dual cell $D_{\sigma}$ associated to a face ${\sigma}\in{\mathcal{E}}$ is defined as follows:

if ${\sigma}=K|L\in{\mathcal{E}}_{\mathrm{int}}$ then $D_{\sigma}=D_{K,{\sigma}}\cup D_{L,{\sigma}}$, where $D_{K,{\sigma}}$ (resp. $D_{L,{\sigma}}$) is the half-part of $K$ (resp. $L$) adjacent to ${\sigma}$ (see Fig. 1 for the two-dimensional case);

• -

if ${\sigma}\in{\mathcal{E}}_{\mathrm{ext}}$ is adjacent to the cell $K$, then $D_{\sigma}=D_{K,{\sigma}}$.

We obtain $d$ partitions of the computational domain $\Omega$ as follows:

[TABLE]

and the $i^{th}$ of these partitions is called $i^{th}$ dual mesh, and is associated to the $i^{th}$ velocity component, in a sense which is clarified below. The set of the faces of the $i^{th}$ dual mesh is denoted by ${{\widetilde{{\mathcal{E}}}}^{(i)}}$ (note that these faces may be orthogonal to any vector of the basis of $\xR^{d}$ and not only ${\boldsymbol{e}}^{(i)}$) and is decomposed into the internal and boundary edges: ${{\widetilde{{\mathcal{E}}}}^{(i)}}={{\widetilde{{\mathcal{E}}}}^{(i)}_{{\rm int}}}\cup{{\widetilde{{\mathcal{E}}}}^{(i)}_{{\rm ext}}}$. The dual face separating two duals cells $D_{\sigma}$ and $D_{{\sigma}^{\prime}}$ is denoted by $\epsilon={\sigma}|{\sigma}^{\prime}$.

To define the scheme, we need some additional notations. The set of faces of a primal cell $K$ and a dual cell $D_{\sigma}$ are denoted by ${\mathcal{E}}(K)$ and ${\widetilde{{\mathcal{E}}}}(D_{\sigma})$ respectively. For ${\sigma}\in{\mathcal{E}}$, we denote by $\boldsymbol{x}_{\sigma}$ the mass center of ${\sigma}$. The vector $\boldsymbol{n}_{K,{\sigma}}$ stands for the unit normal vector to ${\sigma}$ outward $K$. In some case, we need to specify the orientation of a geometrical quantity with respect to the axis:

a primal cell $K$ will be denoted $K=[\overrightarrow{{\sigma}{\sigma}^{\prime}}]$ if ${\sigma},{\sigma}^{\prime}\in{{\mathcal{E}}^{(i)}}\cap{\mathcal{E}}(K)$ for some i\in\bigl{[}|1,d|\bigr{]} are such that $(\boldsymbol{x}_{{\sigma}^{\prime}}-\boldsymbol{x}_{\sigma})\cdot{\boldsymbol{e}}^{(i)}>0$;

• -

we write ${\sigma}=\overrightarrow{K|L}$ if ${\sigma}\in{{\mathcal{E}}^{(i)}}$ and $\overrightarrow{\boldsymbol{x}_{K}\boldsymbol{x}_{L}}\cdot{\boldsymbol{e}}^{(i)}>0$ for some i\in\bigl{[}|1,d|\bigr{]};

• -

the dual face $\epsilon$ separating $D_{\sigma}$ and $D_{{\sigma}^{\prime}}$ is written $\epsilon=\overrightarrow{{\sigma}|{\sigma}^{\prime}}$ if $\overrightarrow{\boldsymbol{x}_{\sigma}\boldsymbol{x}_{{\sigma}^{\prime}}}\cdot{\boldsymbol{e}}^{(i)}>0$ for some i\in\bigl{[}|1,d|\bigr{]}.

For the definition of the discrete momentum diffusion operator, we associate to any dual face $\epsilon$ a distance $d_{\epsilon}$ as sketched on Figure 1. For a dual face $\epsilon\in{\widetilde{{\mathcal{E}}}}(D_{\sigma})$, ${\sigma}\in{{\mathcal{E}}^{(i)}}$, i\in\bigl{[}|1,d|\bigr{]}, the distance $d_{\epsilon}$ is defined by:

[TABLE]

where $d(\cdot,\cdot)$ denotes the Euclidean distance in $\xR^{d}$.

The size $h_{\mathcal{M}}$ and the regularity $\eta_{\mathcal{M}}$ of the mesh are defined by:

[TABLE]

where $|\cdot|$ stands for the $(d-1)$-dimensional measure of a subset of $\xR^{d-1}$ (in the sequel, it is also used to denote the $d$-dimensional measure of a subset of $\xR^{d}$).

The discrete velocity unknowns are associated to the velocity cells and are denoted by $(u_{\sigma})_{{\sigma}\in{{\mathcal{E}}^{(i)}}}$, i\in\bigl{[}|1,d|\bigr{]}, while the discrete pressure unknowns are associated to the primal cells and are denoted by $(p_{K})_{K\in{\mathcal{M}}}$. The discrete pressure space $L_{\mathcal{M}}$ is defined as the set of piecewise constant functions over each of the grid cells $K$ of ${\mathcal{M}}$, and the discrete $i^{th}$ velocity space $H_{{\mathcal{E}}^{(i)}}$ as the set of piecewise constant functions over each of the grid cells $D_{\sigma},\ {\sigma}\in{{\mathcal{E}}^{(i)}}$. The set of functions of $L_{\mathcal{M}}$ with zero mean value is denoted by $L_{{\mathcal{M}},0}$. As in the continuous case, the Dirichlet boundary conditions are (partly) incorporated into the definition of the velocity spaces, by means of the introduction of the spaces H_{{\mathcal{E}}^{(i)},0}\subset H_{{\mathcal{E}}^{(i)}},\ i\in\bigl{[}|1,d|\bigr{]}, defined as follows:

[TABLE]

We then set ${\boldsymbol{H}}_{{\mathcal{E}},0}=\prod_{i=1}^{d}H_{{\mathcal{E}}^{(i)},0}$. Defining the characteristic function ${1\!\!1}_{A}$ of any subset $A\subset\Omega$ by ${1\!\!1}_{A}(\boldsymbol{x})=1$ if $\boldsymbol{x}\in A$ and ${1\!\!1}_{A}(\boldsymbol{x})=0$ otherwise, the $d$ components of a function ${\boldsymbol{u}}\in{\boldsymbol{H}}_{{\mathcal{E}},0}$ and a function $p\in L_{\mathcal{M}}$ may then be written:

[TABLE]

Let us now introduce the discrete operators which are used to write the numerical scheme.

Discrete Laplace operator – For i\in\bigl{[}|1,d|\bigr{]}, the $i^{th}$ component of the discrete Laplace operator is defined by:

[TABLE]

where $d_{\epsilon}$ is defined by (5). The numerical diffusion flux is conservative:

[TABLE]

The discrete Laplace operator of the full velocity vector is defined by

[TABLE]

Let us now recall the definition of the discrete $H^{1}_{0}$-inner product [11]: the $H^{1}_{0}$-inner product between ${\boldsymbol{u}}\in{\boldsymbol{H}}_{{\mathcal{E}},0}$ and ${\boldsymbol{v}}\in{\boldsymbol{H}}_{{\mathcal{E}},0}$ is obtained by taking, for each dual cell, the inner product of the discrete Laplace operator applied to ${\boldsymbol{u}}$ by the test function ${\boldsymbol{v}}$ and integrating over the computational domain. A simple reordering of the sums (which may be seen as a discrete integration by parts) yields, thanks to the conservativity of the diffusion flux (9):

[TABLE]

The bilinear forms

[TABLE]

are inner products on $H_{{\mathcal{E}}^{(i)},0}$ and ${\boldsymbol{H}}_{{\mathcal{E}},0}$ respectively, which induce the following scalar and vector discrete $H^{1}_{0}$ norms:

[TABLE]

This inner product may also be formulated as the $L^{2}$-inner product of discrete gradients. To this purpose, we introduce $d\times d$ new partitions of the domain $\Omega$, where the $(i,j)^{th}$ partition consists in an union of rectangles ($d=2$) or orthogonal parallelepipeds ($d=3$) associated to the dual faces orthogonal to ${\boldsymbol{e}}^{(j)}$ of the dual mesh ${{\widetilde{{\mathcal{E}}}}^{(i)}}$ for the $i^{th}$ component of the velocity. This $(i,j)^{th}$ partition reads:

[TABLE]

where $\boldsymbol{x}_{{\sigma},\epsilon}$ is defined as the orthogonal projection of $\boldsymbol{x}_{\sigma}$ on $\epsilon$ (which is also, in two space dimensions, the vertex of ${\sigma}$ lying on $\epsilon$). The discrete derivative $\eth_{j}u_{i}$ is defined on the $(i,j)^{th}$ partition and reads:

[TABLE]

with $d_{\epsilon}$ defined by (5). These definitions are illustrated on Figure 2. Note that some of these partitions are the same: the $(i,j)^{th}$ partition coincide with the $(j,i)^{th}$ and $(i,i)^{th}$ partitions are the same for i\in\bigl{[}|1,d|\bigr{]}. In addition, these latters also coincide with the primal mesh: for any sub-volume $D_{\epsilon}$ of such a partition, there is $K\in{\mathcal{M}}$ such that $D_{\epsilon}=K$, and we may thus write equivalently $(\eth_{i}u_{i})_{D_{\epsilon}}$ or $(\eth_{i}u_{i})_{K}$. We choose this latter notation in the definition of the discrete divergence below for the sake of consistency, since, if we adopt a variational point of view for the description of the scheme, the discrete velocity divergence has to belong (and indeed does belong) to the space of discrete pressures (see Sections 3 and 4 below for a varitional form of the scheme, in the steady and time-dependent case, respectively). The discrete discrete gradient of each velocity component $u_{i}$ may now be defined as:

[TABLE]

With this definition, it is easily seen that

[TABLE]

If we extend this definition to the velocity vector by

[TABLE]

we get

[TABLE]

This operator satisfies the following consistency result.

Lemma 2.2** (Consistency of the discrete partial derivatives of the velocity).**

Let $\Pi_{\mathcal{E}}$ be an interpolation operator from $C_{c}^{\infty}(\Omega)^{d}$ to ${\boldsymbol{H}}_{{\mathcal{E}},0}$ such that, for any ${\boldsymbol{\varphi}}=(\varphi_{1},\cdots,\varphi_{d})\in C_{c}^{\infty}(\Omega)^{d}$, there exists $C_{\boldsymbol{\varphi}}\geq 0$ depending only on ${\boldsymbol{\varphi}}$ such that

[TABLE]

Let $\eta_{\mathcal{M}}$ be the parameter measuring the regularity of the mesh defined by (6). Then there exists $C_{{\boldsymbol{\varphi}},\eta_{\mathcal{M}}}\geq 0$, only depending in a non-decreasing way on $\eta_{\mathcal{M}}$, such that

[TABLE]

As a consequence, if $({\mathcal{M}}_{m},{\mathcal{E}}_{m})_{m\in\xN}$ is a sequence of MAC grids whose regularity is bounded and whose size tends to 0 as $m$ tends to $+\infty$, then ${\boldsymbol{\nabla}}_{{\widetilde{{\mathcal{E}}}}_{m}}\,(\Pi_{{\mathcal{E}}_{m}}{\boldsymbol{\varphi}})\to{\boldsymbol{\nabla}}{\boldsymbol{\varphi}}$ uniformly as $m\to+\infty$.

Discrete divergence and gradient operators – The discrete divergence operator ${\mathrm{div}}_{\mathcal{M}}$ is defined by:

[TABLE]

Note that the numerical flux is conservative, \ie

[TABLE]

We can now define the discrete divergence-free velocity space:

[TABLE]

The discrete divergence of ${\boldsymbol{u}}=(u_{1},\ldots,u_{d})\in{\boldsymbol{H}}_{{\mathcal{E}},0}$ may also be written as

[TABLE]

where the discrete derivative $(\eth_{i}u_{i})_{K}$ is defined by Relation (13).

The gradient (which applies to the pressure) in the discrete momentum balance equation is built as the dual operator of the discrete divergence, and reads:

[TABLE]

where $\eth_{i}p\in H_{{\mathcal{E}}^{(i)},0}$ is the discrete derivative of $p$ in the $i$-th direction, defined by:

[TABLE]

Note that, in fact, the discrete gradient of a function of $L_{\mathcal{M}}$ should only be defined on the internal faces, and does not need to be defined on the external faces; it is chosen to be in ${\boldsymbol{H}}_{{\mathcal{E}},0}$ (that is zero on the external faces) for the sake of simplicity. Again, the definition of the discrete derivatives of the pressure on the MAC grid is consistent in the sense made precise in the following lemma.

Lemma 2.3** (Discrete gradient consistency).**

Let $\Pi_{\mathcal{M}}$ be an interpolation operator from $C_{c}^{\infty}(\Omega)$ to $L_{\mathcal{M}}$ such that, for any $\psi\in C_{c}^{\infty}(\Omega)$, there exists $C_{\psi}\geq 0$ depending only on $\psi$ such that

[TABLE]

then there exists $C_{\psi,\eta_{\mathcal{M}}}\geq 0$ depending only on $\psi$ and, in a non-decreasing way, on $\eta_{\mathcal{M}}$, such that

[TABLE]

Lemma 2.4** (Discrete ${\mathrm{div}}-{\boldsymbol{\nabla}}$ duality).**

Let $\ q\in L_{\mathcal{M}}$ and ${\boldsymbol{v}}\in{\boldsymbol{H}}_{{\mathcal{E}},0}$ then:

[TABLE]

Proof.

Let $\ q\in L_{\mathcal{M}}$ and ${\boldsymbol{v}}\in{\boldsymbol{H}}_{{\mathcal{E}},0}$. By the definition (19) of the discrete divergence operator and thanks to the conservativity (21) of the flux:

[TABLE]

Therefore, by the definition (23) of the discrete derivative of $q$,

[TABLE]

which concludes the proof. ∎

Discrete convection operator – Let us consider the momentum equation (1b) for the $i^{th}$ component of the velocity, and integrate it on a dual cell $D_{\sigma}$, ${\sigma}\in{{\mathcal{E}}^{(i)}}$. By the Stokes formula, we then need to discretize $\sum_{\epsilon\in{\widetilde{{\mathcal{E}}}}(D_{\sigma})}\int_{\epsilon}u_{i}\,{\boldsymbol{u}}\cdot\boldsymbol{n}_{{\sigma},\epsilon}\ \mathrm{d}\gamma(\boldsymbol{x}),$ where $\boldsymbol{n}_{{\sigma},\epsilon}$ denotes the unit normal vector to $\epsilon$ outward $D_{\sigma}$ and $\ \mathrm{d}\gamma(\boldsymbol{x})$ denotes the $d-1$-dimensional Lebesgue measure. For $\epsilon={\sigma}|{\sigma}^{\prime}$, the convection flux $\int_{\epsilon}u_{i}{\boldsymbol{u}}\cdot\boldsymbol{n}_{{\sigma},\epsilon}\ \mathrm{d}\gamma(\boldsymbol{x})$ is approximated by $|\epsilon|\,u_{{\sigma},\epsilon}\,u_{\epsilon}^{\ast}$; usually, $u_{\epsilon}^{\ast}$ is chosen as the mean value of the two unknowns $u_{\sigma}$ and $u_{{\sigma}^{\prime}}$. In some situations (high Reynolds number for instance), an upwind choice may be preferred. The two possible choices that will be considered for $u_{\epsilon}$ are thus:

[TABLE]

The quantity $|\epsilon|\,u_{{\sigma},\epsilon}$ is the numerical mass flux through $\epsilon$ outward $D_{\sigma}$; it must be chosen carefully to obtain the $L^{2}$-stability of the scheme. More precisely, a discrete counterpart of ${\mathrm{div}}{\boldsymbol{u}}=0$ should be satisfied also on the dual cells. To define $u_{{\sigma},\epsilon}$ on internal dual edges, we distinguish two cases (see Figure 3):

First case – The vector ${\boldsymbol{e}}^{(i)}$ is normal to $\epsilon$, and $\epsilon$ is included in a primal cell $K$, with ${{\mathcal{E}}^{(i)}}(K)=\{{\sigma},{\sigma}^{\prime}\}$. Then the mass flux through $\epsilon={\sigma}|{\sigma}^{\prime}$ is given by:

[TABLE]

Note that, in this relation, all the measures of the face are the same, so this definition equivalently reads $u_{{\sigma},\epsilon}=(-u_{K,{\sigma}}+u_{K,{\sigma}^{\prime}})/2$.

• -

Second case – The vector ${\boldsymbol{e}}^{(i)}$ is tangent to $\epsilon$, and $\epsilon$ is the union of the halves of two primal faces ${\tau}$ and ${\tau}^{\prime}$ such that ${\sigma}=K|L$, ${\tau}\in{\mathcal{E}}(K)$ and ${\tau}^{\prime}\in{\mathcal{E}}(L)$. The mass flux through $\epsilon$ is then given by:

[TABLE]

Again, the numerical flux on a dual face is conservative:

[TABLE]

Moreover, if ${\mathrm{div}}_{\mathcal{M}}\,{\boldsymbol{u}}=0$, the following discrete free divergence condition holds on the dual cells:

[TABLE]

On the external dual faces associated to free degrees of freedom (which means that we are in the second of the above cases), this definition yields $u_{{\sigma},\epsilon}=0$, which is consistent with the boundary condition (1c).

The $i$-th component $C_{{\mathcal{E}}^{(i)}}({\boldsymbol{u}})$ of the non linear convection operator is defined by:

[TABLE]

where $v_{\epsilon}^{\ast}$ is chosen centred or upwind, as defined in (26). The full discrete convection operator ${\boldsymbol{C}}_{\mathcal{E}}({\boldsymbol{u}}),\ {\boldsymbol{H}}_{{\mathcal{E}},0}\longrightarrow{\boldsymbol{H}}_{{\mathcal{E}},0}$ is defined by

[TABLE]

3. The steady case

3.1. The scheme

With the notations introduced in the previous sections, the MAC scheme for the discretization of the steady Navier-Stokes equations (1) on a MAC grid $({\mathcal{M}},{\mathcal{E}})$ reads:

[TABLE]

The discrete right-hand side of the momentum balance equation reads ${\boldsymbol{f}}=\mathcal{P}_{\mathcal{E}}\bar{\boldsymbol{f}}$, where $\mathcal{P}_{\mathcal{E}}$ is the cell mean-value operator defined by $\mathcal{P}_{\mathcal{E}}{\boldsymbol{v}}=(\mathcal{P}_{{\mathcal{E}}^{(1)}}v_{1},\cdots,\mathcal{P}_{{\mathcal{E}}^{(d)}}v_{d})\in H_{{\mathcal{E}}^{(1)},0}\times\cdots\times H_{{\mathcal{E}}^{(d)},0}$ and, for i\in\bigl{[}|1,d|\bigr{]},

[TABLE]

Let us define the weak form $b_{\mathcal{E}}$ of the nonlinear convection term:

[TABLE]

We can now introduce a weak formulation of the scheme, which reads:

[TABLE]

This formulation is equivalent to the strong form (32).

Remark 3.1* (Convergence of the MAC scheme for the Stokes problem and the gradient schemes theory).*

Omitting the convection terms in (35), we obtain a weak formulation of the MAC scheme for the linear Stokes problem. Moreover, formulating the discrete $H^{1}$-inner product as the integral over $\Omega$ of dot products of discrete gradients, the MAC scheme can be interpreted as a gradient scheme in the sense introduced in [13] (see [15] and [8] for more details on the generalization of this formulation to other schemes). Thanks to this result, the (strong) convergence of the velocity and of its discrete gradient to the exact velocity and its gradient can be shown, and thus also the strong convergence of the pressure.

3.2. Stability and existence of a solution

To prove the scheme stability, it is convenient to first reformulate the trilinear form associated to the velocity convection term. To this purpose, we introduce a reconstruction of the velocity components on the partitions which where used for the definition of the discrete velocity gradient. This leads to define $d\times d$ class of reconstruction operators, denoted by $\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(i,j)}$, with $\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(i,j)}$ acting on the $i^{th}$ component of the velocity and providing a reconstruction of this field on the partition of $\Omega$ associated to its $j^{th}$ partial derivative.

Definition 3.2** (Velocity reconstructions).**

Let $({\mathcal{M}},{\mathcal{E}})$ be a given MAC mesh, and let i,j\in\bigl{[}|1,d|\bigr{]}. Let $\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(i,j)}$ be a reconstruction operator defined as follows:

[TABLE]

where $(\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(i,j)}v)_{D_{\epsilon}}$ is a convex combination of the (one of two) discrete values of the $i^{th}$ component of the velocity lying on faces of $D_{\epsilon}$ (see Figure 4).

Such a reconstruction operator satisfies the following stability result.

Lemma 3.3** (Stability of the velocity reconstruction operators).**

Let $({\mathcal{M}},{\mathcal{E}})$ be a given MAC mesh, i,j\in\bigl{[}|1,d|\bigr{]}, and $\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(i,j)}$ be a reconstruction operator, in the sense of Definition 3.2. Then, for $p\in[1,+\infty)$, there exists $C_{\eta_{\mathcal{M}}}\geq 0$, depending only on $p$ and on the parameter $\eta_{\mathcal{M}}$ characterizing the regularity of the mesh defined by (7), and non-decreasing with respect to $\eta_{\mathcal{M}}$, such that, for any $v\in H_{{\mathcal{E}}^{(i)},0}$,

[TABLE]

Proof.

Let $p\in[1,+\infty)$, i,j\in\bigl{[}|1,d|\bigr{]} and $v\in H_{{\mathcal{E}}^{(i)},0}$. We have:

[TABLE]

Since $|\alpha_{\epsilon}|\leq 1$ and $(a+b)^{p}\leq 2^{p-1}(a^{p}+b^{p})$, for $a,b\in[0,+\infty)$, we get:

[TABLE]

Reordering the sums, we obtain that:

[TABLE]

where the volume $V_{\sigma}$ is the sum of the two volumes $D_{\epsilon}$ such that $\epsilon$ is a face of $D_{\sigma}$. It may now be easily checked that there exists $C_{\eta_{\mathcal{M}}}$ depending only on the the parameter $\eta_{\mathcal{M}}$ and non-decreasing with respect to this parameter such that $|V_{\sigma}|\leq C_{\eta_{\mathcal{M}}}|D_{\sigma}|$, which concludes the proof. ∎

The discretization of the velocity convection term in the $i^{th}$ momentum balance equation may be seen as a discrete counterpart of ${\mathrm{div}}(v_{i}\,{\boldsymbol{u}})$, where $v_{i}$ is the convected component of the velocity field (in the scheme, $v_{i}=u_{i}$). Multipling this expression by $w_{i}$ and inegrating over $\Omega$ yields a continuous counterpart $b^{(i)}({\boldsymbol{u}},v_{i},w_{i})$ of $b_{{\mathcal{E}}^{(i)}}({\boldsymbol{u}},v_{i},w_{i})$ which reads:

[TABLE]

An integration by parts (supposing that $w_{i}$ vanishes on the boundary) yields:

[TABLE]

The following lemma states a discrete equivalent of this relation.

Lemma 3.4** (Reformulation of $b_{\mathcal{E}}$).**

Let $({\mathcal{M}},{\mathcal{E}})$ be a given MAC mesh, i\in\bigl{[}|1,d|\bigr{]}, and $({\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}})\in{\boldsymbol{H}}_{{\mathcal{E}},0}\times{\boldsymbol{H}}_{{\mathcal{E}},0}\times{\boldsymbol{H}}_{{\mathcal{E}},0}$. Let $b_{{\mathcal{E}}^{(i)}}({\boldsymbol{u}},v_{i},w_{i})$ be given by (34). Then there exists two reconstruction operators in the sense of Definition 3.2, denoted by $(\mathcal{R}_{{\widetilde{{\mathcal{E}}}}_{n}}^{(i,j)})^{u}$ and $(\mathcal{R}_{{\widetilde{{\mathcal{E}}}}_{n}}^{(j,i)})^{v}$, such that:

[TABLE]

Proof.

Let $({\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}})\in{\boldsymbol{E}}_{\mathcal{E}}\times{\boldsymbol{H}}_{{\mathcal{E}},0}^{2}$. By definition,

[TABLE]

Reordering the sums, we get by conservativity:

[TABLE]

The sum is over the whole set of dual faces ${{\widetilde{{\mathcal{E}}}}^{(i)}}$, so over the $d$ partitions involved in the definition of the discrete gradient of $w_{i}$. In addition, without loss of generality, we may suppose that we have chosen for $\epsilon$ the orientation such that $\epsilon=\overrightarrow{{\sigma}|{\sigma}^{\prime}}$. Hence, we get, by definition (13),

[TABLE]

where $j$ is the index such that $\epsilon$ is normal no ${\boldsymbol{e}}^{(j)}$. For the centered version of the convection operators, $v^{\ast}_{\epsilon}=(v_{\sigma}+v_{{\sigma}^{\prime}})/2$; in the upwind case, it is equal to either $v_{\sigma}$ or $v_{{\sigma}^{\prime}}$ (Relation (26)). In both cases, it is a convex combination of the two discrete values of $v_{i}$ lying on the faces of $D_{\epsilon}$; there exists thus an operator $\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(i,j)}$ (still with the same meaning for $j$), in the sense of Definition 3.2 such that $v^{\ast}_{\epsilon}=(\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(i,j)}v_{i})_{D_{\epsilon}}$. Finally, from the definition of the convection operator and with the chosen orientation for $\epsilon$, $u_{{\sigma},\epsilon}$ is a convex combination of the two values of $u_{j}$ lying on the faces on $D_{\epsilon}$: either the mean value given by (27), if $j=i$, either the convex combination of (28), if $j\neq i$. In addition, $D_{\epsilon}$ is a volume used in the definition of the $j^{th}$ discrete partial derivative of the $i^{th}$ component, and thus also a volume used in the definition of the $i^{th}$ discrete partial derivative of the $j^{th}$ component (both partitions are the same). So there exists one reconstruction operator $\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(j,i)}$ such that $u_{{\sigma},\epsilon}=(\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(j,i)}u_{j})_{D_{\epsilon}}$, which concludes the proof. ∎

Lemma 3.5** (Estimates on $b_{\mathcal{E}}$).**

Let $({\mathcal{M}},{\mathcal{E}})$ be a MAC grid and let $b_{\mathcal{E}}$ be defined by (34). There exists $C_{\eta_{\mathcal{M}}}>0$, depending only and in non-decreasing way on the regularity parameter $\eta_{\mathcal{M}}$ of the mesh defined by (7), such that:

[TABLE]

and

[TABLE]

Proof.

Let i\in\bigl{[}|1,d|\bigr{]}. Thanks to Lemma 3.4, there exists two reconstruction operators in the sense of Definition 3.2, denoted by $(\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(i,j)})^{v}$ and $(\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(j,i)})^{u}$ such that:

[TABLE]

Thanks to Hölder’s inequality, we get, for j\in\bigl{[}|1,d|\bigr{]}:

[TABLE]

which, in view of Lemma 3.3 and the identity (15), concludes the proof of Estimate (36). We then deduce (37) by the discrete Sobolev inequality [11, Lemma 3.5] which allows to control the $L^{4}$-norm by the discrete $H^{1}$-norm. ∎

Let us now prove that $b_{\mathcal{E}}$ is skew-symmetrical with respect to the last two variables. At the continuous level, this result is obtained as follows. For i\in\bigl{[}|1,d|\bigr{]}, on one side, we have by integration by parts:

[TABLE]

On the other side, since ${\mathrm{div}}\,{\boldsymbol{u}}=0$, ${\mathrm{div}}(v_{i}\,{\boldsymbol{u}})={\boldsymbol{u}}\cdot{\boldsymbol{\nabla}}v_{i}$, so:

[TABLE]

which yields the conclusion. The following lemma states a discrete analogue of this property.

Lemma 3.6** ($b_{\mathcal{E}}$ is skew-symmetrical).**

Let $({\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}})\in{\boldsymbol{E}}_{\mathcal{E}}\times{\boldsymbol{H}}_{{\mathcal{E}},0}\times{\boldsymbol{H}}_{{\mathcal{E}},0}$, and let i\in\bigl{[}|1,d|\bigr{]}. Assume the centred choice for $v_{\epsilon}$ in the expression of $b_{{\mathcal{E}}^{(i)}}$; then

[TABLE]

and therefore,

[TABLE]

Assume now the upwind choice for $v_{\epsilon}$ in the expression of $b_{{\mathcal{E}}^{(i)}}$; then,

[TABLE]

Proof.

We mimick the computation performed in the continuous case. At the discrete level and for the centred formulation of the convection term, we have, by a simple reordering of the sum:

[TABLE]

This relation is just obtained by conservativity of the mass flux, by a process which may be seen as a discrete integration by parts, and we have seen that it may be written as a discrete analogue of (38) (Lemma 3.4). On the other hand, thanks to (30) (\ie the discrete analogue of ${\mathrm{div}}\,{\boldsymbol{u}}=0$), we have, for any face ${\sigma}\in{\mathcal{E}}_{\mathrm{int}}^{(i)}$:

[TABLE]

Hence:

[TABLE]

This concludes the proof of (39) and (40). In the upwind case, we have, for i\in\bigl{[}|1,d|\bigr{]}, $b_{{\mathcal{E}}^{(i)}}({\boldsymbol{u}},v_{i},v_{i})=T_{1}+T_{2}$ with:

[TABLE]

From (39), we know that $T_{1}=0$. By definition, for $\epsilon={\sigma}|{\sigma}^{\prime}$,

[TABLE]

Thus, reordering the sums:

[TABLE]

∎

In order to obtain an a priori estimate on the pressure, we introduce a so-called Fortin interpolation operator, \ie a continuous operator from $H^{1}_{0}(\Omega)^{d}$ to ${\boldsymbol{H}}_{\mathcal{E}}$ (equipped with the discrete $H^{1}$-norm) which preserves the divergence. The following lemma is given in [18, Theorem 1, case $q=2$], and we re-state it here with our notations for the sake of clarity.

Lemma 3.7** (Fortin interpolation operator).**

Let $({\mathcal{M}},{\mathcal{E}})$ be a MAC grid of $\Omega$. For ${\boldsymbol{v}}\in H^{1}_{0}(\Omega)^{d}$, we define $\widetilde{\mathcal{P}}_{\mathcal{E}}{\boldsymbol{v}}$ by $\widetilde{\mathcal{P}}_{\mathcal{E}}{\boldsymbol{v}}=(\widetilde{\mathcal{P}}_{{\mathcal{E}}^{(1)}}v_{1},\cdots,\widetilde{\mathcal{P}}_{{\mathcal{E}}^{(d)}}v_{d})\in{\boldsymbol{H}}_{\mathcal{E}}$, where, for i\in\bigl{[}|1,d|\bigr{]},

[TABLE]

For $q\in L^{2}(\Omega)$, we define $\mathcal{P}_{\mathcal{M}}q\in L_{\mathcal{M}}$ by:

[TABLE]

Let ${\boldsymbol{\varphi}}\in(H_{0}^{1}(\Omega))^{d}$. Then:

[TABLE]

In particular, if ${\mathrm{div}}{\boldsymbol{\varphi}}=0$, then ${\mathrm{div}}_{\mathcal{M}}\,(\widetilde{\mathcal{P}}_{\mathcal{E}}{\boldsymbol{\varphi}})=0$. In addition, there exists a real number $C_{\eta_{\mathcal{M}}}$, depending only on $\Omega$ and, in a non-decrasing way, on $\eta_{\mathcal{M}}$ defined by (7), such that:

[TABLE]

Theorem 3.8** (Existence and estimates).**

There exists a solution to (35), and there exists $C_{\eta_{\mathcal{M}}}>0$ depending only on $\Omega$ and, in a non-decreasing way, on the parameter $\eta_{\mathcal{M}}$ characterizing the regularity of the mesh, such that any solution of (35) satisfies the following stability estimate:

[TABLE]

Proof.

Let us start by an a priori estimate on the approximate velocity. Assume that $({\boldsymbol{u}},p)\in{\boldsymbol{H}}_{{\mathcal{E}},0}\times L_{{\mathcal{M}},0}$ satisfies (32); taking ${\boldsymbol{v}}={\boldsymbol{u}}$ in (35a) we get that:

[TABLE]

Since ${\mathrm{div}}_{\mathcal{M}}\,{\boldsymbol{u}}=0$ and $b_{\mathcal{E}}({\boldsymbol{u}},{\boldsymbol{u}},{\boldsymbol{u}})=0$, this yields that

[TABLE]

thanks to the fact that $\|{\boldsymbol{f}}\|_{(L^{2}(\Omega))^{d}}\leq\|\bar{\boldsymbol{f}}\|_{(L^{2}(\Omega))^{d}}$ and to the discrete Poincaré inequality [11, Lemma 9.1].

An a priori estimate on the pressure is obtained by remarking as in [30] that the MAC scheme is inf-sup stable, which is a consequence of the existence of a Fortin operator. Indeed, since $p\in L_{0}^{2}(\Omega)$, there exists ${\boldsymbol{\varphi}}\in H_{0}^{1}(\Omega)^{d}$ such that ${\mathrm{div}}{\boldsymbol{\varphi}}=p$ a.e. in $\Omega$ and

[TABLE]

where $C_{2}$ depends only on $\Omega$ [25]. Taking ${\boldsymbol{v}}=\widetilde{\mathcal{P}}_{\mathcal{E}}{\boldsymbol{\varphi}}$ (defined by (42)) as test function in (35a), we obtain thanks to Lemma 3.7 that

[TABLE]

Thanks to the estimate (37) on $b_{\mathcal{E}}$ and the Cauchy-Schwarz inequality we get:

[TABLE]

where the real number $C_{2}$ is a non-decreasing function of $\eta_{\mathcal{M}}$. This yields

[TABLE]

with $C_{\eta_{\mathcal{M}}}$ non-decreasing with respect to $\eta_{\mathcal{M}}$, thanks to (45), (48) and to the estimate (47).

Let us now prove the existence of a solution to (35). Consider the continuous mapping:

[TABLE]

where $(\hat{{\boldsymbol{u}}},\hat{p})\in{\boldsymbol{H}}_{{\mathcal{E}},0}\times L_{{\mathcal{M}},0}$ is such that:

[TABLE]

It is easily checked that $F$ is well defined, since the values of \hat{u}_{i},\ i\in\bigl{[}|1,d|\bigr{]}, and $\hat{p}$ are readily obtained by setting, for i\in\bigl{[}|1,d|\bigr{]} and ${\sigma}\in{\mathcal{E}}_{\mathrm{int}}^{(i)}$, $v_{i}={1\!\!1}_{D_{\sigma}}$, $v_{j}=0,\ j\neq i$ in (50a) and $q={1\!\!1}_{K}$ in (50b). We also note that the constraint $\hat{p}\in L_{{\mathcal{M}},0}$ is satisfied, thanks to the boundary conditions on ${\boldsymbol{u}}$ (choose $q=1$ in (50b)). The mapping $F$ is continuous; moreover, if $({\boldsymbol{u}},p)\in{\boldsymbol{H}}_{{\mathcal{E}},0}\times L_{{\mathcal{M}},0}$ is such that $F({\boldsymbol{u}},p,\zeta)=(0,0)$, then for any $({\boldsymbol{v}},q)\in{\boldsymbol{H}}_{{\mathcal{E}},0}\times L_{\mathcal{M}}$,

[TABLE]

The arguments used in the above estimates on possible solutions of (35) may be used in a similar way to show that such a pair $({\boldsymbol{u}},p)$ is bounded independently of $\zeta$. Since $F({\boldsymbol{u}},p,0)=0$ is a bijective affine function by the stability of the linear Stokes problem (see [2]), the existence of at least one solution $({\boldsymbol{u}},p)$ to the equation $F({\boldsymbol{u}},p,1)=0$, which is exactly (35), follows by a topological degree argument (see [7] for the theory, [10] for the first application to a nonlinear scheme and [14, Theorem 4.3] for an easy formulation of the result which can be used here). ∎

3.3. Convergence analysis

Lemma 3.9** (Convergence of the velocity reconstructions).**

*Let $({\mathcal{M}}_{n},{\mathcal{E}}_{n})_{n\in\xN}$ be a sequence of MAC meshes such that $h_{{\mathcal{M}}_{n}}\to 0$ as ${n\to+\infty}$ ; assume that there exists $\eta>0$ such that $\eta_{{\mathcal{M}}_{n}}\leq\eta$ for any $n\in\xN$ (with $\eta_{{\mathcal{M}}_{n}}$ defined by (7)). Let i,j\in\bigl{[}|1,d|\bigr{]}, let $\bar{v}\in L^{2}(\Omega)$, and let $(v_{n})_{n\in\xN}$ be such that $v_{n}\in H_{{{\mathcal{E}}^{(i)}}_{n},0}$ and $v_{n}$ converges to $\bar{v}$ as ${n\to+\infty}$ in $L^{2}(\Omega)$. Let $\mathcal{R}_{{\widetilde{{\mathcal{E}}}}_{n}}^{(i,j)}$ be a velocity reconstruction operator, in the sense of Definition (3.2).

Then $\mathcal{R}_{{\widetilde{{\mathcal{E}}}}_{n}}^{(i,j)}v_{n}\to\bar{v}$ in $L^{2}(\Omega)$ as ${n\to+\infty}$.*

Proof.

Let i,j\in\bigl{[}|1,d|\bigr{]}. Denoting $\mathcal{R}_{{\widetilde{{\mathcal{E}}}}_{n}}^{(i,j)}$ by $\mathcal{R}_{n}$ and $\mathcal{P}_{{\mathcal{E}}^{(i)}_{n}}$ (defined by (33)) by $\mathcal{P}_{n}$ for short, we have, for any $\varphi\in C_{c}^{\infty}(\Omega)$:

[TABLE]

Since $\mathcal{R}_{n}v_{n}=\mathcal{R}_{n}\circ\mathcal{P}_{n}v_{n}$, and thanks to the fact that $\|\mathcal{R}_{n}\|_{L^{2}(\Omega)}$ is bounded (see Lemma 3.3) and that $\mathcal{P}_{n}$ is an $L^{2}$-orthogonal projection, we get that there exists $C\geq 0$ such that

[TABLE]

Let $\varepsilon>0$. Let us choose $\varphi\in C_{c}^{\infty}(\Omega)$ such that $\|\mathcal{\varphi}-\bar{v}\|_{L^{2}(\Omega)}\leq\varepsilon/(C+1)$. There exists $n_{1}$ such that $C\|v_{n}-\bar{v}\|_{L^{2}(\Omega)}\leq\varepsilon,\,\forall n\geq n_{1}$, and there exists $n_{2}$ such that $\|\mathcal{R}_{n}\circ\mathcal{P}_{n}\varphi-\varphi\|_{L^{2}(\Omega)}\leq\varepsilon,\,\forall n\geq n_{2}$. Therefore $\|\mathcal{R}_{n}v_{n}-\bar{v}\|_{L^{2}(\Omega)}\leq 3\varepsilon$ for $n\geq\max(n_{1},n_{2})$, which concludes the proof. ∎

Lemma 3.10** (Weak consistency of the nonlinear convection term).**

Let $({\mathcal{M}}_{n},{\mathcal{E}}_{n})_{n\in\xN}$ be a sequence of meshes such that $h_{{\mathcal{M}}_{n}}\to 0$ as $n\to+\infty$ ; assume that there exists $\eta>0$ such that $\eta_{{\mathcal{M}}_{n}}\leq\eta$ for any $n\in\xN$ (with $\eta_{{\mathcal{M}}_{n}}$ defined by (7)). Let $({\boldsymbol{v}}_{n})_{n\in\xN}$ and $({\boldsymbol{w}}_{n})_{n\in\xN}$ be two sequences of functions such that

${\boldsymbol{v}}_{n}\in{\boldsymbol{H}}_{{\mathcal{E}}_{n},0}$* and ${\boldsymbol{w}}_{n}\in{\boldsymbol{H}}_{{\mathcal{E}}_{n},0}$, for $n\in\xN$,*

• -

the sequences $({\boldsymbol{v}}_{n})_{n\in\xN}$ and $({\boldsymbol{w}}_{n})_{n\in\xN}$ converge in $L^{2}(\Omega)^{d}$ to $\bar{\boldsymbol{v}}$ and $\bar{\boldsymbol{w}}$ respectively.

*Let $(\Pi_{{\mathcal{E}}_{n}})_{n\in\xN}$ be a family of interpolation operators satisfying (16) and let ${\boldsymbol{\varphi}}\in C_{c}^{\infty}(\Omega)^{d}$.

Then*

[TABLE]

Proof.

We have $b_{{\mathcal{E}}_{n}}({\boldsymbol{v}}_{n},{\boldsymbol{w}}_{n},\Pi_{{\mathcal{E}}_{n}}{\boldsymbol{\varphi}})=\sum_{i=1}^{d}b_{{\mathcal{E}}^{(i)}}({\boldsymbol{v}},w_{i},\Pi_{{\mathcal{E}}^{(i)}}\varphi_{i})$, where we have omitted the sub- and superscripts $n$ for the sake of clarity in the right-hand side of the equality, with, thanks to Lemma 3.4:

[TABLE]

where $(\mathcal{R}_{{\widetilde{{\mathcal{E}}}}_{n}}^{(i,j)})^{v}$ and $(\mathcal{R}_{{\widetilde{{\mathcal{E}}}}_{n}}^{(j,i)})^{w}$ are two reconstruction operators, in the sense of Definition 3.2. Thanks to the convergence properties of the reconstruction operators (Lemma 3.9) and the strong consistency of the discrete partial derivatives of the velocity (Lemma 2.2), we obtain:

[TABLE]

which concludes the proof. ∎

Lemma 3.11** (Weak consistency of the nonlinear convection term, continued).**

Let $({\mathcal{M}}_{n},{\mathcal{E}}_{n})_{n\in\xN}$ be a sequence of meshes such that $h_{{\mathcal{M}}_{n}}\to 0$ as $n\to+\infty$ ; assume that there exists $\eta>0$ such that $\eta_{{\mathcal{M}}_{n}}\leq\eta$ for any $n\in\xN$ (with $\eta_{{\mathcal{M}}_{n}}$ defined by (7)). Let $({\boldsymbol{u}}_{n})_{n\in\xN}$, $({\boldsymbol{v}}_{n})_{n\in\xN}$ and $({\boldsymbol{w}}_{n})_{n\in\xN}$ be three sequences of functions such that

$({\boldsymbol{u}}_{n},{\boldsymbol{v}}_{n},{\boldsymbol{w}}_{n})\in{\boldsymbol{H}}_{{\mathcal{E}}_{n},0}^{3}$, for $n\in\xN$,

• -

the sequences $({\boldsymbol{u}}_{n})_{n\in\xN}$ and $({\boldsymbol{v}}_{n})_{n\in\xN}$ converge in $L^{p}(\Omega)^{d}$, $1\leq p<6$, to $\bar{\boldsymbol{u}}$ and $\bar{\boldsymbol{v}}$ respectively,

• -

the sequence $({\boldsymbol{w}}_{n})_{n\in\xN}$ converge in $L^{p}(\Omega)^{d}$, $1\leq p<6$, to $\bar{\boldsymbol{w}}\in H^{1}_{0}(\Omega)^{d}$, and $({\boldsymbol{\nabla}}_{{\widetilde{{\mathcal{E}}}}_{n}}{\boldsymbol{w}}_{n})_{n\in\xN}$ converges to ${\boldsymbol{\nabla}}\bar{\boldsymbol{w}}$ weakly in $L^{2}(\Omega)^{d\times d}$.

Then

[TABLE]

Proof.

Once again, we use the reformulation of the form $b_{\mathcal{E}}$, provided by Lemma 3.4. Omitting sub- and superscripts $n$ for short, we have:

[TABLE]

where $(\mathcal{R}_{{\widetilde{{\mathcal{E}}}}_{n}}^{(i,j)})^{u}$ and $(\mathcal{R}_{{\widetilde{{\mathcal{E}}}}_{n}}^{(j,i)})^{v}$ are two reconstruction operators, in the sense of Definition 3.2. Thanks to the stability and convergence properties of the reconstruction operators (Lemma 3.3 and 3.9), the sequences $((\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(j,i)})^{u}u_{j,n})_{n\in\xN}$ and $((\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(i,j)})^{v}v_{n,i})_{n\in\xN}$ are uniformly bounded in $L^{p}(\Omega)^{d}$, for $1\leq p<6$ and i,j\in\bigl{[}|1,d|\bigr{]}, and converge in $L^{2}(\Omega)^{d}$ to $\bar{\boldsymbol{u}}$ and $\bar{\boldsymbol{v}}$, respectively. Hence, these sequences also converge in $L^{2}(\Omega)^{d}$, $1\leq p<6$, and the result follows thanks to the weak convergence in $L^{2}(\Omega)^{d\times d}$ of the partial derivatives. ∎

Lemma 3.12** (A discrete integration by parts formula).**

Let $({\mathcal{M}},{\mathcal{E}})$ be a given MAC mesh, and i,j\in\bigl{[}|1,d|\bigr{]}. Let $u$ and $v$ be two functions of $H_{{\mathcal{E}}^{(i)},0}$. Then there exists a reconstruction operator, in the sense of Definition 3.2, such that:

[TABLE]

Proof.

Let i,j\in\bigl{[}|1,d|\bigr{]} and $(u,v)\in H_{{\mathcal{E}}^{(i)},0}$. We have, by conservativity:

[TABLE]

where $\eta_{\epsilon}=\pm 1$, depending on the relative locations of ${\sigma}$ and ${\sigma}$. For any real number $\alpha_{\epsilon}\in[0,1]$, we have:

[TABLE]

We thus have $T_{1}+T_{2}$=0, with:

[TABLE]

When $i=j$, all the dual faces are included in the domain (so the last sum vanishes). In addition, a dual face $\epsilon$ is included in a cell of the primal mesh, say $K$, and $D_{\epsilon}=K$; we choose in this case $\alpha_{\epsilon}=|D_{K,{\sigma}^{\prime}}|/|K|$ and, by definition of the half-diamond cells, $1-\alpha=|D_{K,{\sigma}}|/|K|$. With this choice, we obtain:

[TABLE]

with, for $\epsilon={\sigma}|{\sigma}^{\prime}$,

[TABLE]

Let us choose now consider the case $i\neq j$. In this case, we choose $\alpha=|D_{{\sigma}^{\prime}}|/(2\ |D_{\epsilon}|)$, so $1-\alpha=|D_{\sigma}|/(2\ |D_{\epsilon}|$, by definition of $D_{\epsilon}$, and we get (52) with:

[TABLE]

∎

We are now in position to state and prove the convergence of the scheme.

Theorem 3.13** (Convergence of the scheme, steady case).**

Let $({\mathcal{M}}_{n},{\mathcal{E}}_{n})_{n\in\xN}$ be a sequence of meshes such that $h_{{\mathcal{M}}_{n}}\to 0$ as $n\to+\infty$ ; assume that there exists $\eta>0$ such that $\eta_{{\mathcal{M}}_{n}}\leq\eta$ for any $n\in\xN$ (with $\eta_{{\mathcal{M}}_{n}}$ defined by (7)). Let $({\boldsymbol{u}}_{n},p_{n})$ be a solution to the MAC scheme (32) or its weak form (35), for ${\mathcal{M}}={\mathcal{M}}_{n}$. Then there exists $\bar{\boldsymbol{u}}\in H^{1}_{0}(\Omega)^{d}$ and $\bar{p}\in L^{2}(\Omega)$ such that, up to a subsequence:

the sequence $({\boldsymbol{u}}_{n})_{n\in\xN}$ converges to $\bar{\boldsymbol{u}}$ in $L^{2}(\Omega)^{d}$,

• -

the sequence $({\boldsymbol{\nabla}}_{{\widetilde{{\mathcal{E}}}}_{n}}{\boldsymbol{u}}_{n})_{n\in\xN}$ converges to ${\boldsymbol{\nabla}}\bar{\boldsymbol{u}}$ in $L^{2}(\Omega)^{d\times d}$,

• -

the sequence $(p_{n})_{n\in\xN}$ converges to $\bar{p}$ in $L^{2}(\Omega)$,

• -

$(\bar{\boldsymbol{u}},\bar{p})$* is a solution to the weak formulation of the steady Navier-Stokes equations (2).*

Proof.

Thanks to the estimate (47) on the velocity, applying the classical estimate on the translates [11, Theorem 14.2] we obtain the existence of a subsequence of approximate solutions $({\boldsymbol{u}}_{n})_{n\in\xN}$ which converges to some $\bar{\boldsymbol{u}}\in L^{2}(\Omega)^{d}$. From the estimates on the translates, we also get the regularity of the limit, that is $\bar{\boldsymbol{u}}\in H_{0}^{1}(\Omega)^{d}$. The estimate (49) on the pressure then yields the weak convergence of a subsequence of $(p_{n})_{n\in\xN}$ to some $\bar{p}$ in $L^{2}(\Omega)$. Let us then pass to the limit in the scheme in order to prove its (weak) consistency.

Passing to the limit in the mass balance equation – Let $\psi\in C_{c}^{\infty}(\Omega)$. Taking $\psi_{n}=\Pi_{{\mathcal{M}}_{n}}\psi$, the pointwise interpolate defined by (24), as test function in (35b) and using (25), we get that:

[TABLE]

Therefore, thanks to Lemma 2.3,

[TABLE]

and therefore $\bar{\boldsymbol{u}}$ satisfies (35b).

Passing to the limit in the momentum balance equation – Let ${\boldsymbol{\varphi}}=(\varphi_{1},\cdots,\varphi_{d})\in C_{c}^{\infty}(\Omega)^{d}$, and take ${\boldsymbol{\varphi}}_{n}=\Pi_{{\mathcal{E}}_{n}}{\boldsymbol{\varphi}}=(\varphi_{n,1},\cdots,\varphi_{n,d})\in{\boldsymbol{H}}_{{\mathcal{E}}_{n},0}$ as test function in (35a). This yields:

[TABLE]

Thanks to the weak $L^{2}$-convergence of $p_{n}$ to $\overline{p}$ and to the uniform convergence of $\mathcal{P}_{{\mathcal{E}}_{n}}\bar{\boldsymbol{f}}$ to $\bar{\boldsymbol{f}}$ and of ${\mathrm{div}}_{{\mathcal{M}}_{n}}\,{\boldsymbol{\varphi}}_{n}$ to ${\mathrm{div}}{\boldsymbol{\varphi}}$ (see Lemma 2.2) as $n\to+\infty$, we have

[TABLE]

From [11, Proof of Theorem 9.1], thanks to the $L^{2}$-convergence of ${\boldsymbol{u}}_{n}$ to $\bar{\boldsymbol{u}}$, we get that, for i\in\bigl{[}|1,d|\bigr{]},

[TABLE]

Therefore,

[TABLE]

By Lemma 3.10, we have

[TABLE]

Passing to the limit as $n\to+\infty$ in (53) thus yields that $\bar{\boldsymbol{u}}$ and $\bar{p}$ satisfy (2).

Strong convergence of ${\boldsymbol{\nabla}}_{{\widetilde{{\mathcal{E}}}}_{n}}{\boldsymbol{u}}_{n}$ to ${\boldsymbol{\nabla}}\bar{\boldsymbol{u}}$ in $L^{2}(\Omega)^{d\times d}$ – The sequence $({\boldsymbol{\nabla}}_{{\widetilde{{\mathcal{E}}}}_{n}}{\boldsymbol{u}}_{n})_{n\in\xN}$ is bounded in $L^{2}(\Omega)^{d\times d}$ and therefore, there exists $\xi\in L^{2}(\Omega)^{d\times d}$ and a subsequence still denoted by $({\boldsymbol{\nabla}}_{{\widetilde{{\mathcal{E}}}}_{n}}{\boldsymbol{u}}_{n})_{{n\in\xN}}$ converging to $\xi$ weakly in $L^{2}(\Omega)^{d\times d}$. Let i,j\in\bigl{[}|1,d|\bigr{]}, and let $\varphi$ be a function of $C^{\infty}_{c}(\Omega)$. We denote by $\varphi_{n}$ the interpolate of $\varphi$ by the projection operator $\Pi_{{\mathcal{E}}^{(i)}_{n}}$ associated to the $i^{th}$ component of the velocity. By Lemma 3.12, we know that there exists a reconstruction operator $\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(i,j)}$, in the sense of Definition 3.2, such that:

[TABLE]

By the strong convergence of $\varphi_{n}$ to $\varphi$, of $\mathcal{R}_{\widetilde{{\mathcal{E}}}}^{(i,j)}u_{n,i}$ to $\bar{u}_{i}$ and of $\eth_{j}\varphi_{n}$ to $\partial_{j}\varphi$, passing to the limit in the above relation, we get:

[TABLE]

Integrating by parts in the right-hand side thanks to the regularity of $\bar{\boldsymbol{u}}$, we obtain:

[TABLE]

Hence, by density, $\xi={\boldsymbol{\nabla}}\bar{\boldsymbol{u}}$. Taking ${\boldsymbol{\varphi}}_{n}={\boldsymbol{u}}_{n}$ in (53) yields:

[TABLE]

Passing to the limit as $n\to+\infty$ we get that:

[TABLE]

which implies the strong convergence of the discrete gradient of the velocity.

Strong convergence of the pressure – Let ${\boldsymbol{\varphi}}_{n}\in H_{0}^{1}(\Omega)^{d}$ be such that ${\mathrm{div}}{\boldsymbol{\varphi}}_{n}=p_{n}$ a.e. in $\Omega$ and $\|{\boldsymbol{\varphi}}_{n}\|_{H_{0}^{1}(\Omega)^{d}}\leq C\ \|p_{n}\|_{L^{2}(\Omega)},$ where $C$ depends only on $\Omega$. Let ${\boldsymbol{\psi}}_{n}=\widetilde{\mathcal{P}}_{{\mathcal{E}}_{n}}{\boldsymbol{\varphi}}_{n}$; thanks to Lemma 3.7, we have $\|{\boldsymbol{\psi}}_{n}\|_{1,{\mathcal{E}}_{n},0}\leq C\ C_{\eta_{n}}\,\|p_{n}\|_{L^{2}(\Omega)}$, and since $p_{n}\in L_{{\mathcal{M}}_{n}}$, we get that ${\mathrm{div}}_{{\mathcal{M}}_{n}}\,{\boldsymbol{\psi}}_{n}=p_{n}$. Therefore, taking ${\boldsymbol{\psi}}_{n}=\widetilde{\mathcal{P}}_{{\mathcal{E}}_{n}}{\boldsymbol{\varphi}}_{n}$ as test function in (35a), we obtain:

[TABLE]

From the bound on $\|{\boldsymbol{\psi}}_{n}\|_{1,{\mathcal{E}},0}$, we know that ${\boldsymbol{\psi}}_{n}$ converges to some ${\boldsymbol{\psi}}\in H^{1}_{0}(\Omega)^{d}$ in $L^{2}(\Omega)^{d}$ and, by the same arguments as for the identification of $\xi$ with ${\boldsymbol{\nabla}}\bar{\boldsymbol{u}}$, that ${\boldsymbol{\nabla}}_{{\widetilde{{\mathcal{E}}}}_{n}}{\boldsymbol{\psi}}_{n}\to{\boldsymbol{\nabla}}{\boldsymbol{\psi}}$ weakly in $L^{2}(\Omega)^{d\times d}$ as $n\to+\infty$. In addition, we also have that ${\mathrm{div}}{\boldsymbol{\varphi}}=p$ a.e. in $\Omega$. By Lemma 3.11, $b_{\mathcal{E}}({\boldsymbol{u}}_{n},{\boldsymbol{u}}_{n},{\boldsymbol{\psi}}_{n})$ converges to $b(\bar{\boldsymbol{u}},\bar{\boldsymbol{u}},{\boldsymbol{\psi}})$. Passing to the limit as $n\to+\infty$, we thus get that

[TABLE]

Since $(\bar{\boldsymbol{u}},\bar{p})$ satisfies (2), this implies that $\|p_{n}\|_{L^{2}(\Omega)}\to\|\bar{p}\|_{L^{2}(\Omega)}$, which in turn yields that $p_{n}\to\bar{p}$ in $L^{2}(\Omega)$ as $n\to+\infty$. ∎

Remark 3.14* (Uniqueness of the continuous solution and convergence of the whole sequence).*

In the case where uniqueness of the solution is known, then a classical argument can be used to show that the whole sequence converges ; this is for instance the case for small data, see e.g. [31, Theorem 1.3] or [3, Theorem V.3.5].

4. The time-dependent case

4.1. Time discretization

Let us now turn to the time discretization of the problem (3). We consider a MAC grid $({\mathcal{M}},{\mathcal{E}})$ of $\Omega$ in the sense of Definition 2.1, and a partition $0=t_{0}<t_{1}<\cdots<t_{N}=T$ of the time interval $(0,T)$, and, for the sake of simplicity, a constant time step $\delta t=t_{n+1}-t_{n}$; hence $t_{n}=n\,\delta t$, for n\in\bigl{[}|0,N|\bigr{]}. Let \{u_{\sigma}^{n+1},\ {\sigma}\in{\mathcal{E}},\ n\in\bigl{[}|0,N-1|\bigr{]}\} and \{p_{K}^{n+1},\ K\in{\mathcal{M}},\ n\in\bigl{[}|0,N-1|\bigr{]}\} be sets of discrete velocity and pressure unknowns. For n\in\bigl{[}|1,N|\bigr{]}, we first define the corresponding piecewise constant space-dependent functions ${\boldsymbol{u}}=(u_{1}^{n},\ldots,u_{d}^{n})$ and $p^{n}$ by:

[TABLE]

We enforce that $u_{\sigma}^{n}=0$ for ${\sigma}\in{\mathcal{E}}_{\mathrm{ext}}$ and n\in\bigl{[}|1,N|\bigr{]} (so $u_{i}^{n}\in H_{{\mathcal{E}}^{(i)},0}$ and the sum in the relation above may be restricted to ${\mathcal{E}}_{\mathrm{int}}^{(i)}$), and we set ${\boldsymbol{u}}^{n}=(u_{1}^{n},\ldots,u_{d}^{n})\in{\boldsymbol{H}}_{{\mathcal{E}},0}$. Then, we define the discrete (time- and space-dependent) velocities and pressures functions by:

[TABLE]

where ${1\!\!1}_{]t_{n},t_{n+1}]}$ is the characteristic function of the interval $]t_{n},t_{n+1}]$. For i\in\bigl{[}|1,d|\bigr{]}, we denote by $X_{{\mathcal{E}},\delta t}^{(i)}$ the set of such piecewise constant functions on time intervals and dual cells for the $i^{th}$ velocity component approximation, we set $\boldsymbol{X}_{{\mathcal{E}},\delta t}=\prod_{i=1}^{d}X_{{\mathcal{E}},\delta t}^{(i)}$, and we denote by $Y_{{\mathcal{M}},\delta t}$ the space of piecewise constant functions on time intervals and primal cells for the pressure approximation. Setting

[TABLE]

we define the discrete time derivative $\eth_{t}{\boldsymbol{u}}\in\boldsymbol{X}_{{\mathcal{E}},\delta t}$ by:

[TABLE]

Finally, we define the discrete right-hand side by:

[TABLE]

With these notations, the time-implicit MAC scheme for the transient Navier-Stokes reads:

[TABLE]

Step $n$, n\in\bigl{[}|0,N-1|\bigr{]}, of the scheme (55) admits the following weak formulation:

[TABLE]

The equivalence between this relation and (55b)-(55d) (in the sense that (56) implies the existence of a discrete pressure field such that (55b)-(55d) is satisfied) is a consequence of the stability of the MAC scheme for the Stokes problem (\ie the fact that this scheme satisfies a discrete inf-sup condition).

4.2. Estimates on discrete solutions and existence

Let us define the two following discrete norms for functions of space and time:

[TABLE]

Lemma 4.1** (Existence and first estimates on the velocity).**

There exists at least a solution ${\boldsymbol{u}}\in\boldsymbol{X}_{{\mathcal{M}},\delta t}$ satisfying (55). Furthermore, there exists $C>0$ depending only on ${{\boldsymbol{u}}_{0}}$ and $\bar{\boldsymbol{f}}$ such that, for any function ${\boldsymbol{u}}\in\boldsymbol{X}_{{\mathcal{M}},\delta t}$ satisfying (55), the following estimates hold:

[TABLE]

Proof.

We prove the a priori estimates (57) and (58). The existence of a solution then follows by a topological degree argument, as for the stationary case.

Let M\in\bigl{[}|0,N-1|\bigr{]}; taking ${\boldsymbol{v}}={\boldsymbol{u}}^{n+1}$ in (56), multiplying by $\delta t$ and summing the result over n\in\bigl{[}|0,M|\bigr{]}, we obtain thanks to Lemma 3.6 and to the Cauchy-Schwarz inequality:

[TABLE]

Using the fact that for all $a,b\in\xR$, $2a(a-b)=(a-b)^{2}+a^{2}-b^{2}$ for the first term of the left-hand side and the discrete Poincaré and Young inequalities for the right-hand side, we get that

[TABLE]

where $C_{P}>0$ depends only on $\Omega$. On one hand, this inequality yields the $L^{\infty}$-estimate (58); on the other hand, taking $M=N-1$, we get the $L^{2}$-estimate (57). ∎

Next we turn to an estimate on the discrete time derivative. To this end, we introduce the following discrete dual norms on ${\boldsymbol{H}}_{{\mathcal{E}},0}$ and $\boldsymbol{X}_{{\mathcal{E}},\delta t}$:

[TABLE]

Lemma 4.2** (Estimate on the dual norm of the velocity discrete time derivative).**

Let ${\boldsymbol{u}}\in\boldsymbol{X}_{{\mathcal{E}},\delta t}$ be a solution to (55). Then there exists $C>0$ depending only on ${{\boldsymbol{u}}_{0}}$, $\Omega$, $\bar{\boldsymbol{f}}$ and, in a non-decreasing way, on $\eta_{\mathcal{M}}$, such that:

[TABLE]

Proof.

Taking ${\boldsymbol{v}}\in{\boldsymbol{E}}_{\mathcal{E}}$ such that $\|{\boldsymbol{v}}\|_{1,{\mathcal{E}},0}\leq 1$ as test function in (56), we have, for n\in\bigl{[}|0,N-1|\bigr{]}:

[TABLE]

By Lemma 3.6 and thanks to the estimate (36), we have:

[TABLE]

Using the Cauchy-Schwarz inequality, we note that:

[TABLE]

Therefore, thanks to the estimate (58) of Lemma 4.1 and to the discrete Poincaré inequality, there exists $\widetilde{C}_{\eta_{\mathcal{M}}}>0$ depending only on $\Omega$ and on the regularity of the mesh, such that:

[TABLE]

Hence,

[TABLE]

Multiplying this latter inequality by $\delta t$ and summing for n\in\bigl{[}|0,N-1|\bigr{]}, we get:

[TABLE]

We conclude by the discrete Sobolev inequality [11, Lemma 3.5] and thanks to the $L^{2}(0,T;{\boldsymbol{H}}_{{\mathcal{E}},0})$-estimate (57) of ${\boldsymbol{u}}$. ∎

4.3. Convergence analysis

Theorem 4.3** (Convergence of the scheme, time-dependent case).**

Let $(\delta t_{m})_{m\in\xN}$ and $({\mathcal{M}}_{m},{\mathcal{E}}_{m})_{m\in\xN}$ be a sequence of time steps and MAC grids (in the sense of Definition 2.1) such that $\delta t_{m}\rightarrow 0$ and $h_{{\mathcal{M}}_{m}}\to 0$ as $m\to+\infty$. Assume that there exists $\eta>0$ such that $\eta_{{\mathcal{M}}_{m}}\leq\eta$ for any $m\in\xN$ (with $\eta_{{\mathcal{M}}_{m}}$ defined by (7)). Let ${\boldsymbol{u}}_{m}$ be a solution to (56) for $\delta t=\delta t_{m}$ and $({\mathcal{M}},{\mathcal{E}})=({\mathcal{M}}_{m},{\mathcal{E}}_{m})$. Then there exists $\bar{\boldsymbol{u}}\in L^{2}(0,T;{\boldsymbol{E}}(\Omega))$ such that, up to a subsequence:

the sequence $({\boldsymbol{u}}_{m})_{m\in\xN}$ converges to $\bar{\boldsymbol{u}}$ in $L^{4/3}(0,T;L^{2}(\Omega)^{d})$,

• -

$\bar{\boldsymbol{u}}$* is a solution to the weak formulation (4).*

• -

$\partial_{t}\bar{\boldsymbol{u}}\in L^{4/3}(0,T;E^{\prime}(\Omega))$.

Proof.

We proceed in four steps.

First step: compactness in $L^{4/3}(0,T;L^{2}(\Omega)^{d})$ – The first step consists in applying the discrete Aubin-Simon theorem 5.3 in order to obtain the existence of a subsequence of $({\boldsymbol{u}}_{m})_{{m\in\xN}}$ which converges to $\bar{\boldsymbol{u}}$ in $L^{4/3}((0,T);L^{2}(\Omega)^{d})$. In our setting, we apply Theorem 5.3 with $p=4/3$; the Banach space $B$ is $L^{2}(\Omega)^{d}$, and the spaces $X_{m}$ and $Y_{m}$ consist in the space ${\boldsymbol{H}}_{{\mathcal{E}}_{m},0}$ endowed with the norms defined respectively by Relations (12) and (59). By [11, Theorem 14.2] and the Kolmogorov compactness theorem (see e.g. [11, Theorem 14.1]), we obtain that $(X_{m},Y_{m})_{m\in\xN}$ is compactly embedded in $B$ in the sense of Definition 5.1. Let us then show that the sequence $({X_{m}},Y_{m})_{{m\in\xN}}$ is compact-continuous in $L^{2}(\Omega)^{d}$ in the sense of Definition 5.2. Let ${\boldsymbol{v}}_{m}\in{\boldsymbol{H}}_{{\mathcal{E}}_{m},0}$ such that $(\|{\boldsymbol{v}}_{m}\|_{1,{\mathcal{E}}_{m},0})_{{m\in\xN}}$ is bounded and $\|{\boldsymbol{v}}_{m}\|_{{\boldsymbol{E}}^{\prime}_{m}}\rightarrow 0$ as $m\to+\infty$. Assume that ${\boldsymbol{v}}_{m}\rightarrow{\boldsymbol{v}}$ in $(L^{2}(\Omega))^{d}$; by definition (59) of the dual norm, we have:

[TABLE]

Passing to the limit in this inequality as $m\to\infty$, we get that ${\boldsymbol{v}}=0$, so that the sequence $({X_{m}},Y_{m})_{{m\in\xN}}$ is compact-continuous in $L^{2}(\Omega)^{d}$. We now check the three assumptions (H1), (H2) and (H3) of Theorem 5.3: by Lemma 4.1, the sequence $\|{\boldsymbol{u}}_{m}\|_{L^{1}(0,T;{\boldsymbol{H}}_{{\mathcal{E}},0})}$ is bounded, and thanks to the discrete Poincaré inequality, the sequence $({\boldsymbol{u}}_{m})_{m\in\xN}$ is also bounded in $L^{4/3}(0,T;(L^{2}(\Omega)^{d}))$; furthermore, the sequence $\|\eth_{t}{\boldsymbol{u}}_{m}\|_{L^{4/3}(0,T;{\boldsymbol{E}}^{\prime}_{\mathcal{E}})}$ is bounded by Lemma 4.2. Hence, Theorem 5.3 applies and there exists $\bar{\boldsymbol{u}}\in L^{4/3}(0,T;L^{2}(\Omega)^{d})$ such that, up to a subsequence,

[TABLE]

Step 2: Convergence in $L^{2}(0,T;L^{2}(\Omega)^{d})$ – Thanks to Lemma 4.1, the sequence $({\boldsymbol{u}}_{m})_{m\in\xN}$ is bounded in $L^{\infty}(0,T,L^{2}(\Omega)^{d})$, and therefore, there exists $\hat{{\boldsymbol{u}}}\in L^{\infty}(0,T;L^{2}(\Omega)^{d})$ and a subsequence $({\boldsymbol{u}}_{\phi(m)})_{m\in\xN}$ converging to $\hat{{\boldsymbol{u}}}$ $\star$-weakly in $L^{\infty}(0,T;L^{2}(\Omega)^{d})$. Since ${\boldsymbol{u}}_{\phi(m)}\rightarrow\bar{{\boldsymbol{u}}}$ in $L^{4/3}(0,T;L^{2}(\Omega)^{d})$, the uniqueness of the limit in the sense of distributions implies that $\bar{\boldsymbol{u}}=\hat{{\boldsymbol{u}}}$ so that $\bar{\boldsymbol{u}}\in L^{\infty}(0,T;L^{2}(\Omega)^{d})$. By a classical interpolation result on $L^{p}(0,T)$ spaces, we have:

[TABLE]

which implies that ${\boldsymbol{u}}_{m}$ converges towards $\bar{\boldsymbol{u}}$ in $L^{2}(0,T;L^{2}(\Omega)^{d})$ as $m$ tends to infinity.

Step 3: Weak consistency of the scheme – The notion of weak consistency that we use here is the Lax-Wendroff notion: we show that if a sequence of approximate solutions of the scheme converges to some limit, then this limit is a weak solution to the original problem. Let us then show that $\bar{\boldsymbol{u}}$ satisfies (4). Let ${\boldsymbol{\varphi}}\in C_{c}^{\infty}(\Omega\times[0,T))^{d}$, such that ${\mathrm{div}}{\boldsymbol{\varphi}}=0$. By Lemma 3.7, we have ${\mathrm{div}}_{{\mathcal{M}}_{m}}\,\widetilde{\mathcal{P}}_{{\mathcal{E}}_{m}}{\boldsymbol{\varphi}}(\cdot,t_{n})=0$, and so we can take ${\boldsymbol{\varphi}}_{m}^{n}=\widetilde{\mathcal{P}}_{{\mathcal{E}}_{m}}{\boldsymbol{\varphi}}(\cdot,t_{n})\in{\boldsymbol{E}}_{\mathcal{E}}$ as test function in (56) ; multiplying by $\delta t_{m}$ and summing for $n=\{0,\ldots,N_{m}-1\}$ (with $N_{m}\delta t_{m}=T$), we then get:

[TABLE]

where the subscript $m$ in ${\boldsymbol{f}}_{m}^{n+1}$ is here to recall that the discrete right-hand side is an interpolation of the continuous one, which depends on the mesh and time step. The first term of the left-hand side reads $T_{1,m}=\sum_{i=1}^{d}T_{1,m,i}$ with:

[TABLE]

We know that $u_{m,i}\to\bar{u}_{i}$ in $L^{2}(0,T;L^{2}(\Omega))$ as $m\to+\infty$. By definition, the discrete partial derivative $\eth_{t}\varphi_{m,i}$ converges uniformly to $\partial_{t}\varphi_{i}$ as $m\to+\infty$. Moreover, $\widetilde{\mathcal{P}}^{(i)}_{{\mathcal{E}}_{m}}\bar{u}_{0,i}$ converges to $\bar{u}_{0,i}$ in $L^{q}(\Omega)$ for all $q$ in $[1,2]$, and $\varphi_{m,{\sigma}}^{0}$ converges to $\varphi_{i}(\cdot,0)$ in $L^{q}(\Omega)$ for all $q$ in $[1,\infty]$. Hence:

[TABLE]

Let us then study the second term of the left-hand side. We have:

[TABLE]

By the same arguments as in the stationary case, we get that

[TABLE]

Moreover, thanks to the regularity of ${\boldsymbol{\varphi}}$,

[TABLE]

where $C_{\varphi}$ only depends on ${\boldsymbol{\varphi}}$. We thus obtain that

[TABLE]

Similarly, we have:

[TABLE]

so that

[TABLE]

The convection term is dealt with by remarking that an easy adaptation of Lemma 3.10 to the time-dependent framework implies that

[TABLE]

Therefore, $\bar{\boldsymbol{u}}$ is indeed a solution of (4).

Step 4: Regularity of the limit – Thanks to [11, Theorems 14.1 and 14.2], the sequence of normed vector spaces $({\boldsymbol{H}}_{{\mathcal{E}}_{m},0},\|\cdot\|_{1,{\mathcal{E}}_{m},0})_{m\in\xN}$ is $L^{2}(\Omega)^{d}$-limit-included in $H_{0}^{1}(\Omega)^{d}$ in the sense of Definition 5.4. We have ${\boldsymbol{u}}_{m}\rightarrow\bar{\boldsymbol{u}}$ in $L^{2}(0,T,L^{2}(\Omega)^{d})$ as $m\to+\infty$ and $(\|{\boldsymbol{u}}_{m}\|_{L^{2}(0,T;{\boldsymbol{H}}_{{\mathcal{E}}_{m},0})})_{{m\in\xN}}$ is bounded thanks to Lemma 4.1. Therefore Theorem 5.5 applies, so that $\bar{\boldsymbol{u}}\in L^{2}(0,T;H^{1}_{0}(\Omega)^{d})$; then, adapting the proof that ${\mathrm{div}}\bar{\boldsymbol{u}}=0$ of the stationary case (see the proof of Theorem 3.13), we get that $\bar{\boldsymbol{u}}\in L^{2}(0,T;{\boldsymbol{E}}(\Omega))$.

Let us finally show that $\partial_{t}\bar{\boldsymbol{u}}\in L^{4/3}(0,T;E^{\prime}(\Omega))$. Let ${\boldsymbol{\varphi}}\in C_{c}^{\infty}(\Omega\times(0,T))$ such that ${\mathrm{div}}{\boldsymbol{\varphi}}=0$. Let ${\boldsymbol{\varphi}}_{m}\in\boldsymbol{X}_{{\mathcal{E}}_{m},\delta t_{m}}$ be defined by

[TABLE]

Note that, for n\in\bigl{[}|0,N-1|\bigr{]}, ${\boldsymbol{\varphi}}_{m}^{n+1}$ is discretely divergence-free, \ie ${\boldsymbol{\varphi}}_{m}^{n+1}\in{\boldsymbol{E}}_{{\mathcal{E}}_{m}}$. Thanks to Lemma 4.2, there exists $C\geq 0$ depending only on ${\boldsymbol{u}}_{0}$, $\Omega$, $\eta$ and $\bar{\boldsymbol{f}}$ such that:

[TABLE]

By Lemma 3.7, there exists $C_{2}$ depending only on $\eta$ and $\Omega$, such that $\|{\boldsymbol{\varphi}}_{m}\|_{L^{4}(0,T;{\boldsymbol{H}}_{{\mathcal{E}},0})}\leq C_{2}\|{\boldsymbol{\varphi}}\|_{L^{4}(0,T;{\boldsymbol{E}}(\Omega))},$ where ${\boldsymbol{E}}(\Omega)$ is endowed with the $H^{1}_{0}$ norm. Hence, passing to the limit as $m\to+\infty$ in a similar way as for $T_{1,m}$ in Step 3, we get that

[TABLE]

We then obtain that $\partial_{t}\bar{\boldsymbol{u}}\in L^{4/3}(0,T;{\boldsymbol{E}}^{\prime}(\Omega))$ by density (see [31, Theorem 1.6] for the density of divergence-free regular functions in divergence-free functions of $H^{1}_{0}(\Omega)^{d}$). ∎

Remark 4.4* (Uniqueness and convergence of the whole sequence).*

In the case where uniqueness of the solution is known, then again the whole sequence converges ; this is for instance the case for $d=2$, see e.g. [31, Theorem 3.2], under a small data assumption [31, Theorem 3.7] or under a short time assumption [31, Theorem 3.11].

4.4. Case of the unsteady Stokes equations

In the case of the unsteady Stokes equations, that is Problem (3) where the nonlinear convection term in (3b) is omitted, stronger estimates can be obtained, which entail the weak convergence of the pressure. To obtain these bounds, the assumption that ${{\boldsymbol{u}}_{0}}\in H^{1}(\Omega)^{d}$ and that ${\mathrm{div}}{{\boldsymbol{u}}_{0}}=0$ plays a central role.

Let us consider the following weak formulation of the unsteady Stokes problem:

[TABLE]

Note that this formulation does not use divergence-free test functions as in (4), so the pressure still appears.

The scheme – We look for an approximation $({\boldsymbol{u}},p)\in\boldsymbol{X}_{{\mathcal{E}},\delta t}\times Y_{{\mathcal{M}},\delta t}$ of $({\boldsymbol{u}},p)$ solution to the problem (61); we consider the time-implicit MAC scheme which reads:

[TABLE]

Note that the choice of the discretization of the initial condition in (62a), together with the assumption ${\mathrm{div}}{\boldsymbol{u}}_{0}=0$, implies that ${\mathrm{div}}_{\mathcal{M}}\,{\boldsymbol{u}}^{0}=0$; this fact is important for the obtention of the estimates. A weak formulation of (62b)–(62d) reads:

[TABLE]

The estimates of Lemma 4.1 on the approximate solutions obtained in the case of the Navier-Stokes equations are of course still valid. However we get stronger estimates on $\eth_{t}{\boldsymbol{u}}$ and on $p$, as we proceed to show.

Lemma 4.5** (Estimates on the velocity and its discrete time derivative).**

Let ${\boldsymbol{u}}\in\boldsymbol{X}_{{\mathcal{E}},\delta t}$ be a solution to (62); then there exists $C>0$ depending only on ${{\boldsymbol{u}}_{0}}$, $\Omega$, $\bar{\boldsymbol{f}}$ and, in a non-decreasing way, on $\eta_{\mathcal{M}}$, such that:

[TABLE]

Proof.

Let ${\boldsymbol{u}}^{n+1}\in{\boldsymbol{E}}_{\mathcal{E}}$ be a solution to (63). Taking ${\boldsymbol{v}}=\eth_{t}{\boldsymbol{u}}^{n+1}$ as test function, we get:

[TABLE]

By linearity of the discrete time derivative and the discrete divergence operators, and thanks to (62d), we get that ${\mathrm{div}}_{\mathcal{M}}\,(\eth_{t}{\boldsymbol{u}}^{n+1})=\eth_{t}({\mathrm{div}}_{\mathcal{M}}\,{\boldsymbol{u}}^{n+1})=0$. Multiplying (66) by $\delta t$ and summing the result over n\in\bigl{[}|0,M|\bigr{]}, for M\in\bigl{[}|0,N-1|\bigr{]}, we obtain $T_{1}+T_{2}=T_{3}$ where

[TABLE]

We have, by linearity of the discrete gradient operator:

[TABLE]

By continuity of the Fortin operator, we have in addition that $|{\boldsymbol{u}}^{0}\|_{1,{\mathcal{E}},0}\leq C\ \|{\boldsymbol{u}}_{0}\|_{H^{1}(\Omega)^{d}}$, with $C$ depending only on $\Omega$ and (in a non-decreasing way) on $\eta_{\mathcal{M}}$. Let us now turn to $T_{3}$. By the Cauchy-Schwarz and the Young inequalities, we obtain:

[TABLE]

and the Cauchy-Schwarz inequality, together with the definition of ${\boldsymbol{f}}$, yields for the first term at the right-hand side:

[TABLE]

Gathering the above inequalities, we get that:

[TABLE]

This in turn yields the $L^{2}$-estimate (64) (taking $M=N-1$) on the discrete time derivative of the velocity, and the $L^{\infty}(H^{1})$-estimate (65) on the velocity itself. ∎

Lemma 4.6** (Estimate on the pressure).**

Let $({\boldsymbol{u}},p)\in\boldsymbol{X}_{{\mathcal{M}},\delta t}\times Y_{{\mathcal{M}},\delta t}$ be a solution to (62). There exists $C\geq 0$ depending only on $\Omega$, $\bar{\boldsymbol{f}}$ and, in a non-decreasing way, on $\eta_{\mathcal{M}}$, such that:

[TABLE]

Proof.

We follow the same strategy as in the proof of the pressure estimate in Proposition 3.8. Therefore, let ${\boldsymbol{\varphi}}\in H^{1}_{0}(\Omega)^{d}$ be such that ${\mathrm{div}}{\boldsymbol{\varphi}}=p^{n+1}$ and $\|{\boldsymbol{\nabla}}{\boldsymbol{\varphi}}\|_{L^{2}(\Omega)^{d\times d}}\leq C\ \|p^{n+1}\|_{L^{2}(\Omega)}$, with $C$ depending only on $\Omega$. Taking ${\boldsymbol{v}}=\widetilde{\mathcal{P}}_{\mathcal{E}}{\boldsymbol{\varphi}}$ as test function in (63), we obtain, thanks to (44):

[TABLE]

Thanks to the Cauchy-Schwarz and Poincaré inequalities and to the continuity of the Fortin operator $\widetilde{\mathcal{P}}_{\mathcal{E}}$, we then get that there exists $C_{\eta_{\mathcal{M}}}$ depending on $\Omega$ and on the regularity of the mesh such that

[TABLE]

Summing this relation over n\in\bigl{[}|0,N-1|\bigr{]} and multiplying by $\delta t$ yields the result thanks to (57) and (64). ∎

Theorem 4.7** (Convergence of the scheme, time-dependent Stokes problem).**

Let $(\delta t)_{m\in\xN}$ and $({\mathcal{M}}_{m},{\mathcal{E}}_{m})_{m\in\xN}$ be a sequence of time steps and meshes such that $(\delta t)_{m}\to 0$ and $h_{{\mathcal{M}}_{m}}\to 0$ as $m\to+\infty$; assume that there exists $\eta>0$ such that $\eta_{{\mathcal{M}}_{m}}\leq\eta$ for any $m\in\xN$ (with $\eta_{{\mathcal{M}}_{m}}$ defined by (7)). Let $({\boldsymbol{u}}_{m},p_{m})$ be a solution to (62) for $\delta t=\delta t_{m}$ and ${\mathcal{M}}={\mathcal{M}}_{m}$. Then there exists $(\bar{\boldsymbol{u}},\bar{p})\in L^{2}(0,T;{\boldsymbol{E}}(\Omega))\times L^{2}(0,T;L^{2}(\Omega))$ such that, up to a subsequence:

the sequence $({\boldsymbol{u}}_{m})_{m\in\xN}$ converges to $\bar{\boldsymbol{u}}$ in $L^{2}(0,T;L^{2}(\Omega)^{d})$,

• -

the sequence $(p_{m})_{m\in\xN}$ weakly converges to $\bar{p}$ in $\in L^{2}(0,T;L^{2}(\Omega))$,

• -

$(\bar{\boldsymbol{u}},\bar{p})$* is a solution to the weak formulation (61).*

Proof.

The convergence of the sequence of discrete solutions of the velocity follow from Theorem 4.3 and the weak convergence of the sequence of discrete solutions of the pressure in $L^{2}(0,T;L^{2}(\Omega))$ follow from the estimate (68). Let us then show that $(\bar{\boldsymbol{u}},\bar{p})$ satisfies (61). Let ${\boldsymbol{\varphi}}\in C_{c}^{\infty}(\Omega\times[0,T))^{d}$. Taking ${\boldsymbol{\varphi}}_{m}^{n}=\widetilde{\mathcal{P}}_{{\mathcal{E}}_{m}}{\boldsymbol{\varphi}}(\cdot,t_{n})\in{\boldsymbol{H}}_{{\mathcal{E}}_{m},0}$ as test function in (63), multiplying by $\delta t_{m}$ and summing for n\in\bigl{[}|0,N_{m}-1|\bigr{]} (with $N_{m}\delta t_{m}=T$), we obtain:

[TABLE]

Let us deal with the pressure term (all other terms of the equation can be dealt with as in the proof of Theorem 4.3). We have, by the divergence preservation property of the Fortin operator:

[TABLE]

Hence, thanks to the regularity of ${\boldsymbol{\varphi}}$ (\ie the fact that $|{\boldsymbol{\varphi}}(\boldsymbol{x},t)-{\boldsymbol{\varphi}}(\boldsymbol{x},t_{n})|\leq C_{\varphi}\,\delta t_{m}$ for $\boldsymbol{x}\in\Omega$ and $t\in(t_{n},t_{n+1})$) and the weak convergence of $p_{m}$ to $\bar{p}$,

[TABLE]

∎

5. Appendix: Discrete functional analysis

Definition 5.1** (Compactly embedded sequence of spaces).**

Let $B$ be a Banach space; a sequence $(X_{m})_{{m\in\xN}}$ of Banach spaces included in $B$ is compactly embedded in $B$ if any sequence $(u_{m})_{m\in\xN}$ satisfying:

• •

$u_{m}\in X_{m}$ ($\forall m\in\xN$),

• •

the sequence $(\|u_{m}\|_{X_{m}})_{{m\in\xN}}$ is bounded,

is relatively compact in $B$.

Definition 5.2** (Compact-continuous sequence of spaces).**

Let B be a Banach space, and let $(X_{m})_{{m\in\xN}}$ and $(Y_{m})_{{m\in\xN}}$ be sequences of Banach spaces such that $X_{m}\subset B$ for $m\in\xN$. The sequence $(X_{m},Y_{m})_{{m\in\xN}}$ is compact-continuous in $B$ if the following conditions are satified:

• $\bullet$

The sequence $(X_{m})_{{m\in\xN}}$ is compactly embedded in $B$ (see Definition 5.1),

• $\bullet$

$X_{m}\subset Y_{m}$ (for all ${m\in\xN}$),

• $\bullet$

if the sequence $(u_{m})_{{m\in\xN}}$ is such that $u_{m}\in X_{m}$ (for all ${m\in\xN}$), $(\|u_{m}\|_{X_{m}})_{{m\in\xN}}$ is bounded and $\|u_{m}\|_{Y_{m}}\rightarrow 0$ as $m\to+\infty$, then any subsequence of $(u_{m})_{{m\in\xN}}$ converging in $B$ converges to [math] (in $B$).

The following theorem is proved [4] and is a generalization of a previous work carried out in [16].

Theorem 5.3** (Aubin-Simon Theorem with a sequence of subspaces and a discrete derivative.).**

Let $1\leq p<\infty$, let B be a Banach space, and let $(X_{m})_{{m\in\xN}}$ and $(Y_{m})_{{m\in\xN}}$ be sequences of Banach spaces such that $X_{m}\subset B$ for $m\in\xN$. We assume that the sequence $(X_{m},Y_{m})_{{m\in\xN}}$ is compact-continuous in $B$. Let $T>0$ and $(u^{(m)})_{{m\in\xN}}$ be a sequence of $L^{p}(0,T;B)$ satisfying the following conditions:

1. $\bullet$

(H1) the sequence $(u^{(m)})_{{m\in\xN}}$ is bounded in $L^{p}(0,T;B)$. 2. $\bullet$

(H2) the sequence $(\|u^{(m)}\|_{L^{1}(0,T;X_{m})})_{{m\in\xN}}$ is bounded. 3. $\bullet$

(H3) the sequence $(\|\eth_{t}u^{(m)}\|_{L^{p}(0,T;Y_{m})})_{{m\in\xN}}$ is bounded.

Then there exists $u\in L^{p}(0,T;B)$ such that, up to a subsequence, $u^{(m)}\rightarrow u$ in $L^{p}(0,T;B)$.

Definition 5.4** ($B$-limit-included).**

Let $B$ be a Banach space, $(X_{m})_{m\in\xN}$ be a sequence of Banach spaces included in $B$ and $X$ be a Banach space included in $B$. The sequence $(X_{m})_{m\in\xN}$ is $B$-limit-included in $X$ if there exists $C\in\xR$ such that if $u$ is the limit in $B$ of a subsequence of a sequence $(u_{m})_{{m\in\xN}}$ verifying $u_{m}\in X_{m}$ and $\|u_{m}\|_{X_{m}}\leq 1$, then $u\in X$ and $\|u\|_{X}\leq C$.

The regularity of a possible limit of approximate solutions may be proved thanks to the theorem which we recall below [17, Theorem B1].

Theorem 5.5** (Regularity of the limit).**

Let $1\leq p<\infty$ and $T>0$. Let $B$ be a Banach space, $(X_{m})_{m\in\xN}$ be a sequence of Banach spaces included in $B$ and $B$-limit-included in $X$ (where $X$ is a Banach space included in $B$). Let $T>0$ and, for $m\in\xN$, Let $u_{m}\in L^{p}(0,T;X_{m})$. We assume that the sequence $(\|u_{m}\|_{L^{p}(0,T;X_{m})})_{m\in\xN}$ is bounded and that $u_{m}\rightarrow u$ a.e. as $m\rightarrow\infty$. Then $u\in L^{p}(0,T;X)$.