Low Mach number limit of a pressure correction MAC scheme for compressible barotropic flows

TL;DR

This paper analyzes how a pressure correction MAC scheme for compressible barotropic flows behaves as the Mach number approaches zero, showing convergence to an incompressible flow scheme under well-prepared initial data.

Contribution

It proves the convergence of a pressure correction MAC scheme for compressible flows to the incompressible limit as Mach number tends to zero, under certain initial conditions.

Findings

Numerical scheme converges to incompressible solution as Mach number approaches zero.

Convergence holds for well-prepared initial data.

Provides theoretical validation for the scheme's asymptotic behavior.

Abstract

We study the incompressible limit of a pressure correction MAC scheme [3] for the unstationary compressible barotropic Navier-Stokes equations. Provided the initial data are well-prepared, the solution of the numerical scheme converges, as the Mach number tends to zero, towards the solution of the classical pressure correction inf-sup stable MAC scheme for the incompressible Navier-Stokes equations.

Equations

$${\partial }_t{\rho }^{\epsilon }+\mathrm{div}({\rho }^{\epsilon }{\mathbf{u}}^{\epsilon })=0,$$

$${\partial }_t({\rho }^{\epsilon }{\mathbf{u}}^{\epsilon })+\mathrm{div}({\rho }^{\epsilon }{\mathbf{u}}^{\epsilon }\otimes {\mathbf{u}}^{\epsilon })-\mathrm{div}(\mathbf{\tau }({\mathbf{u}}^{\epsilon }))+\frac{1}{{\epsilon }^2}\nabla \wp ({\rho }^{\epsilon })=0,$$

$${\mathbf{u}}^{\epsilon }{|}_{\partial \Omega }=0,{\rho }^{\epsilon }{|}_{t=0}={\rho }_0^{\epsilon },{\mathbf{u}}^{\epsilon }{|}_{t=0}={\mathbf{u}}_0^{\epsilon },$$

$$\mathrm{div}(\mathbf{\tau }(\mathbf{u}))=\mu \Delta \mathbf{u}+(\mu +\lambda )\nabla (\mathrm{div}\mathbf{u}),$$

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

$$\bar{\rho }{\partial }_t\bar{\mathbf{u}}+\bar{\rho }\mathrm{div}(\bar{\mathbf{u}}\otimes \bar{\mathbf{u}})-\mu \mathbf{\Delta }\bar{\mathbf{u}}+\nabla \pi =0,$$

$$||{\mathbf{u}}_0^{\epsilon }|{|}_{{\mathrm{H}}^1(\Omega {)}^d}+\frac{1}{\epsilon }||\mathrm{div}{\mathbf{u}}_0^{\epsilon }|{|}_{{\mathrm{L}}^2(\Omega )}+\frac{1}{{\epsilon }^2}||{\rho }_0^{\epsilon }-1|{|}_{{\mathrm{L}}^{\infty }(\Omega )}\le C.$$

$$\displaystyle\displaystyle\hskip 21.52771pt\mbox{For }1\leq i\leq d,\ \forall\sigma\in{\mathcal{E}}_{{\rm int}}^{(i)},\quad(\overline{{\boldsymbol{\nabla}}p})^{n}_{\sigma,i}=\Bigl{(}\frac{\rho^{n}_{\scalebox{0.6}{$D_{\sigma}$}}}{\rho^{n-1}_{\scalebox{0.6}{$D_{\sigma}$}}}\Bigr{)}^{1/2}({\boldsymbol{\nabla}}p^{n})_{\sigma,i}.$$

$${\tilde{\mathbf{u}}}^{n+1}$$

$$\text{mbox}For1\le i\le d,\forall \sigma \in {\mathcal{E}}_{\mathrm{int}}^{(i)},$$

$$\displaystyle\hskip 21.52771pt\dfrac{1}{\delta t}\bigl{(}\rho^{n}_{\scalebox{0.6}{$D_{\sigma}$}}\tilde{u}^{n+1}_{\sigma,i}-\rho^{n-1}_{\scalebox{0.6}{$D_{\sigma}$}}u_{\sigma,i}^{n}\bigr{)}\!+\!{\rm div}(\rho^{n}\tilde{u}^{n+1}_{i}{\boldsymbol{u}}^{n})_{\sigma}\!-\!{\rm div}{\boldsymbol{\tau}}(\tilde{\boldsymbol{u}}^{n+1})_{\sigma,i}\!+\!\dfrac{1}{\varepsilon^{2}}(\overline{{\boldsymbol{\nabla}}p})^{n}_{\sigma,i}\!=0.$$

$${\rho }^{n+1}$$

$$\text{mbox}For1\le i\le d,\forall \sigma \in {\mathcal{E}}_{\mathrm{int}}^{(i)},$$

$$\displaystyle\displaystyle\phantom{\hskip 21.52771pt\forall K\in{\mathcal{M}},\quad}\dfrac{1}{\delta t}\ \rho^{n}_{\scalebox{0.6}{$D_{\sigma}$}}\ (u^{n+1}_{\sigma,i}-\tilde{u}_{\sigma,i}^{n+1})+\dfrac{1}{\varepsilon^{2}}\,({\boldsymbol{\nabla}}p^{n+1})_{\sigma,i}-\dfrac{1}{\varepsilon^{2}}\,(\overline{{\boldsymbol{\nabla}}p})_{\sigma,i}^{n}=0,$$

$$\forall K\in \mathcal{M},\frac{1}{\delta t}({\rho }_K^{n+1}-{\rho }_K^n)+\mathrm{div}({\rho }^{n+1}{\mathbf{u}}^{n+1}{)}_K=0,$$

$$\forall K\in \mathcal{M},p_K^{n+1}=\wp ({\rho }^{n+1}),$$

$$\mathrm{div}(\rho \mathbf{u}{)}_K=\frac{1}{|K|}\sum _{\sigma \in \mathcal{E}(K)}{F}_{K,\sigma }(\rho ,\mathbf{u}),K\in \mathcal{M},$$

$$F_{K,\sigma }(\rho ,\mathbf{u})=|\sigma |{\rho }_{\sigma }{u}_{K,\sigma },\sigma \in {\mathcal{E}}_{\mathrm{int}},\sigma =K|L,$$

$$(\nabla p{)}_{\sigma ,i}=\frac{|\sigma |}{|D_{\sigma }|}(\wp ({\rho }_L)-\wp ({\rho }_K)){\mathbf{n}}_{K,\sigma }\cdot {\mathbf{e}}^{(i)},1\le i\le d,\sigma \in {\mathcal{E}}_{\mathrm{int}}^{(i)},\sigma =K|L.$$

$$\sum _{K\in \mathcal{M}}|K|\wp ({\rho }_K)(\mathrm{div}\mathbf{u}{)}_K+\sum _{i=1}^d\sum _{\sigma \in {\mathcal{E}}_{\mathrm{int}}^{(i)}}|D_{\sigma }|u_{\sigma ,i}(\nabla p{)}_{\sigma ,i}=0.$$

$$||\mathbf{u}|{|}_{1,\mathcal{M}}=1\text{ and }\sum _{K\in \mathcal{M}}|K|p_K(\mathrm{div}\mathbf{u}{)}_K\ge \beta ||p-\frac{1}{|\Omega |}{\int }_{\Omega }p\mathrm{d}\mathbf{x}|{|}_{L^2(\Omega )},$$

$$|D_{\sigma}|\,\rho_{\scalebox{0.6}{$D_{\sigma}$}}=|D_{K,\sigma}|\,\rho_{K}+|D_{L,\sigma}|\,\rho_{L},\qquad\sigma\in{\mathcal{E}}_{{\rm int}},\,\sigma=K|L,$$

$$\mathrm{div}(\rho v_i\mathbf{u}{)}_{\sigma }=\sum _{\epsilon \in \bar{\mathcal{E}}(D_{\sigma })}{F}_{\sigma ,\epsilon }(\rho ,\mathbf{u})v_{i,\epsilon },\sigma \in {\mathcal{E}}_{\mathrm{int}}^{(i)},1\le i\le d.$$

$$F_{\sigma,\varepsilon}(\rho,{\boldsymbol{u}})=\frac{1}{2}\bigl{(}F_{K,\sigma}(\rho,{\boldsymbol{u}})\ {\boldsymbol{n}}_{D_{\sigma},\varepsilon}\cdot{\boldsymbol{n}}_{K,\sigma}+F_{K,\sigma^{\prime}}(\rho,{\boldsymbol{u}})\ {\boldsymbol{n}}_{D_{\sigma},\varepsilon}\cdot{\boldsymbol{n}}_{K,\sigma^{\prime}}\bigr{)}.$$

$$F_{\sigma,\varepsilon}(\rho,{\boldsymbol{u}})=\frac{1}{2}\ \bigl{(}F_{K,\tau}(\rho,{\boldsymbol{u}})+F_{L,\tau^{\prime}}(\rho,{\boldsymbol{u}})\bigr{)}.$$

$$\frac{|D_{\sigma}|}{\delta t}(\rho^{n+1}_{\scalebox{0.6}{$D_{\sigma}$}}-\rho^{n}_{\scalebox{0.6}{$D_{\sigma}$}})+\sum_{\varepsilon\in\bar{\mathcal{E}}(D_{\sigma})}F_{\sigma,\varepsilon}(\rho^{n+1},{\boldsymbol{u}}^{n+1})=0,\qquad\sigma\in{\mathcal{E}}_{{\rm int}}.$$

$$-\sum _{i=1}^d\sum _{\mathcal{E}\in {\mathcal{E}}_{\mathrm{int}}^{(i)}}|D_{\sigma }|u_{\sigma ,i}\mathrm{div}\mathbf{\tau }(\mathbf{u}{)}_{\sigma ,i}\ge \mu ||\mathbf{u}|{|}_{1,\mathcal{M}}^2.$$

$$\forall K\in \mathcal{M},{\rho }_K^0=\frac{1}{|K|}{\int }_K{\rho }_0^{\epsilon }(\mathbf{x})\mathrm{d}\mathbf{x}\text{mbox}and\forall \sigma \in {\mathcal{E}}_{\mathrm{int}}^{(i)},1\le i\le d,u_{\sigma ,i}^0=\frac{1}{|\sigma |}{\int }_{\sigma }{\mathbf{u}}_0^{\epsilon }(\mathbf{x})\cdot {\mathbf{e}}^{(i)}\mathrm{d}\mathbf{x},$$

$$\frac{1}{{\epsilon }^2}\underset{K\in \mathcal{M}}{\max }|{\rho }_K^0-1|+\frac{1}{{\epsilon }^2}\underset{1\le i\le d}{\max }\underset{\sigma \in {\mathcal{E}}_{\mathrm{int}}^{(i)}}{\max }|(\nabla p{)}_{\sigma ,i}^0|+\frac{1}{\epsilon }\underset{K\in \mathcal{M}}{\max }|{\rho }_K^{-1}-1|\le C.$$

$$\forall K\in \mathcal{M},{\rho }_K^{-1}-1=\underset{=\mathcal{O}({\epsilon }^2)}{\underbrace{{\rho }_K^0-1}}+\underset{=\mathcal{O}(\epsilon )}{\underbrace{\delta t{\rho }_K^0(\mathrm{div}{\mathbf{u}}^0{)}_K}}+\underset{=\mathcal{O}({\epsilon }^2)}{\underbrace{\delta t\sum _{\sigma \in \mathcal{E}(K)}\frac{|\sigma |}{|K|}({\rho }_{\sigma }^0-{\rho }_K^0){\mathbf{u}}_{K,\sigma }^0}}.$$

$$\dfrac{1}{2\delta t}\Bigl{(}\rho_{\scalebox{0.6}{$D_{\sigma}$}}^{n}\,|u_{\sigma,i}^{n+1}|^{2}-\rho_{\scalebox{0.6}{$D_{\sigma}$}}^{n-1}\,|u_{\sigma,i}^{n}|^{2}\Bigr{)}+\frac{1}{2|D_{\sigma}|}\sum_{\begin{subarray}{c}\varepsilon\in\bar{\mathcal{E}}(D_{\sigma})\\ \varepsilon=D_{\sigma}|D_{\sigma^{\prime}}\end{subarray}}F_{\sigma,\varepsilon}(\rho^{n},{\boldsymbol{u}}^{n})\,\tilde{u}_{\sigma,i}^{n+1}\,\tilde{u}_{\sigma^{\prime},i}^{n+1}\\ -{\rm div}{\boldsymbol{\tau}}(\tilde{\boldsymbol{u}}^{n+1})_{\sigma,i}\,\tilde{u}_{\sigma,i}^{n+1}+\frac{1}{\varepsilon^{2}}({\boldsymbol{\nabla}}p)_{\sigma,i}^{n+1}\,u_{\sigma,i}^{n+1}+\frac{\delta t}{\varepsilon^{4}}\Bigl{(}\frac{|({\boldsymbol{\nabla}}p)^{n+1}_{\sigma,i}|^{2}}{2\,\rho^{n}_{\scalebox{0.6}{$D_{\sigma}$}}}-\frac{|({\boldsymbol{\nabla}}p)^{n}_{\sigma,i}|^{2}}{2\,\rho^{n-1}_{\scalebox{0.6}{$D_{\sigma}$}}}\Bigr{)}\\ +R_{\sigma,i}^{n+1}=0,\qquad\mbox{with }R_{\sigma,i}^{n+1}=\dfrac{1}{2\delta t}\rho_{\scalebox{0.6}{$D_{\sigma}$}}^{n-1}(\tilde{u}_{\sigma,i}^{n+1}-u_{\sigma,i}^{n})^{2}.$$

$$\dfrac{1}{\delta t}\Bigl{(}\Pi_{\gamma}(\rho^{n+1}_{K})-\Pi_{\gamma}(\rho^{n}_{K})\Bigr{)}+{\rm div}\bigl{(}b_{\gamma}(\rho^{n+1}){\boldsymbol{u}}^{n+1}-b_{\gamma}^{\prime}(1)\rho^{n+1}{\boldsymbol{u}}^{n+1}\bigr{)}_{K}+p_{K}^{n+1}\,{\rm div}({\boldsymbol{u}}^{n+1})_{K}+R_{K}^{n+1}=0,$$

$$\frac{1}{2}\!\sum_{i=1}^{d}\!\!\sum_{\sigma\in{\mathcal{E}}_{{\rm int}}^{(i)}}\!|D_{\sigma}|\!\Big{(}\rho^{n}_{\scalebox{0.6}{$D_{\sigma}$}}\ |u^{n+1}_{\sigma,i}|^{2}\!-\rho^{n-1}_{\scalebox{0.6}{$D_{\sigma}$}}\ |u^{n}_{\sigma,i}|^{2}\Big{)}\!+\!\frac{1}{\varepsilon^{2}}\!\!\sum_{K\in{\mathcal{M}}}\!\!|K|\!\Big{(}\Pi_{\gamma}(\rho^{n+1}_{K})-\Pi_{\gamma}(\rho^{n}_{K})\Big{)}\\ \!+\!\mu\delta t||\tilde{\boldsymbol{u}}^{n+1}||_{1,{\mathcal{M}}}^{2}\!+\!\frac{1}{\varepsilon^{4}}\!\sum_{i=1}^{d}\!\!\sum_{\sigma\in{\mathcal{E}}_{{\rm int}}^{(i)}}\!\!|D_{\sigma}|\delta t^{2}\Bigl{(}\!\frac{|({\boldsymbol{\nabla}}p)^{n+1}_{\sigma,i}|^{2}}{2\,\rho^{n}_{\scalebox{0.6}{$D_{\sigma}$}}}-\frac{|({\boldsymbol{\nabla}}p)^{n}_{\sigma,i}|^{2}}{2\,\rho^{n-1}_{\scalebox{0.6}{$D_{\sigma}$}}}\!\Bigr{)}\!\!+\!\mathcal{R}^{n+1}\!\!\leq 0$$

$$\displaystyle\bullet\ \mbox{For all $\gamma\geq 1$ there exists $C_{\gamma}$ such that:}\ \Pi_{\gamma}(\rho)\leq C_{\gamma}\,|\rho-1|^{2},\ \forall\rho\in(0,2).$$

$$\displaystyle\bullet\ \mbox{If $\gamma\geq 2$ then}\ \Pi_{\gamma}(\rho)\geq|\rho-1|^{2},\ \forall\rho>0.$$

$$\displaystyle\begin{array}[]{l}\bullet\ \mbox{If $\gamma\in[1,2)$ then for all $R\in(2,+\infty)$, there exists $C_{\gamma,R}$ such that:}\\ \hskip 56.9055pt\begin{array}[]{ll}\Pi_{\gamma}(\rho)\geq C_{\gamma,R}|\rho-1|^{2},&\ \forall\rho\in(0,R),\\ \Pi_{\gamma}(\rho)\geq C_{\gamma,R}|\rho-1|^{\gamma},&\ \forall\rho\in[R,\infty).\end{array}\end{array}$$

$$\hskip 30.1388pt\frac{1}{2}\sum_{i=1}^{d}\sum_{\sigma\in{\mathcal{E}}_{{\rm int}}^{(i)}}|D_{\sigma}|\rho^{n-1}_{\scalebox{0.6}{$D_{\sigma}$}}\ |u^{n}_{\sigma,i}|^{2}+\mu\ \sum_{k=1}^{n}\delta t\ ||\tilde{\boldsymbol{u}}^{k}||_{1,{\mathcal{M}}}^{2}\\ +\frac{1}{\varepsilon^{2}}\sum_{K\in{\mathcal{M}}}|K|\,\Pi_{\gamma}(\rho^{n}_{K})+\frac{1}{\varepsilon^{4}}\sum_{i=1}^{d}\sum_{\sigma\in{\mathcal{E}}_{{\rm int}}^{(i)}}\frac{|D_{\sigma}|\ \delta t^{2}}{2\,\rho^{n-1}_{\scalebox{0.6}{$D_{\sigma}$}}}\ |({\boldsymbol{\nabla}}p)^{n}_{\sigma,i}|^{2}\leq C_{0}.\hskip 30.1388pt$$

Full text

1. Introduction

Let $\Omega$ be parallelepiped of $\mathbb{R}^{d}$, with $d\in\{2,3\}$ and $T>0$. The unsteady barotropic compressible Navier-Stokes equations, parametrized by the Mach number $\varepsilon$, read for $({\boldsymbol{x}},t)\in\Omega\times(0,T)$:

[TABLE]

where $\rho^{\varepsilon}>0$ and ${\boldsymbol{u}}^{\varepsilon}=(u_{1}^{\varepsilon},..,u_{d}^{\varepsilon})^{T}$ are the density and velocity of the fluid. The pressure satisfies the ideal gas law $\wp(\rho^{\varepsilon})=(\rho^{\varepsilon})^{\gamma}$, with $\gamma\geq 1$, and

[TABLE]

where the real numbers $\mu$ and $\lambda$ satisfy $\mu>0$ and $\mu+\lambda>0$. The smooth solutions of (1) are known to satisfy a kinetic energy balance and a renormalization identity. In addition, under assumption on the initial data, it may be inferred from these estimates that the density $\rho^{\varepsilon}$ tends to a constant $\bar{\rho}$, and the velocity tends, in a sense to be defined, to a solution $\bar{u}$ of the incompressible Navier-Stokes equations [4]:

[TABLE]

where $\pi$ is the formal limit of $(\wp(\rho^{\varepsilon})-\wp(\bar{\rho}))/\varepsilon^{2}$.

In this paper, we reproduce this theory for a pressure correction scheme, based on the Marker-And-Cell (MAC) space discretization: we first derive discrete analogues of the kinetic energy and renormalization identities, then establish from these relations that approximate solutions of (1) converge, as $\varepsilon\to 0$, towards the solution of the classical projection scheme for the incompressible Navier-Stokes equations (2).

For this asymptotic analysis, we assume that the initial data is “well prepared”: $\rho_{0}^{\varepsilon}>0$, $\rho_{0}^{\varepsilon}\in{\rm L}^{\infty}(\Omega)$, ${\boldsymbol{u}}_{0}^{\varepsilon}\in{\rm H}_{0}^{1}(\Omega)^{d}$ and, taking without loss of generality $\bar{\rho}=1$, there exists $C$ independent of $\varepsilon$ such that:

[TABLE]

Consequently, $\rho_{0}^{\varepsilon}$ tends to $1$ when $\varepsilon\to 0$; moreover, we suppose that ${\boldsymbol{u}}_{0}^{\varepsilon}$ converges in ${\rm L}^{2}(\Omega)^{d}$ towards a function $\bar{\boldsymbol{u}}_{0}\in{\rm L}^{2}(\Omega)^{d}$ (the uniform boundedness of the sequence in the ${\rm H}^{1}(\Omega)^{d}$ norm already implies this convergence up to a subsequence).

2. The numerical scheme

Let ${\mathcal{M}}$ be a MAC mesh (see e.g. [1] and Figure 1 for the notations). The discrete density unknowns are associated with the cells of the mesh ${\mathcal{M}}$, and are denoted by \big{\{}\rho_{K},\ K\in{\mathcal{M}}\big{\}}. We denote by $\mathcal{E}$ the set of the faces of the mesh, and by $\mathcal{E}^{(i)}$ the subset of the faces orthogonal to the $i$-th vector of the canonical basis of $\mathbb{R}^{d}$. The discrete $i^{th}$ component of the velocity is located at the centre of the faces $\sigma\in\mathcal{E}^{(i)}$, so the whole set of discrete velocity unknowns reads \big{\{}u_{\sigma,i},\ \sigma\in\mathcal{E}^{(i)},1\leq i\leq d\big{\}}. We define ${\mathcal{E}}_{{\rm ext}}=\{\sigma\in\mathcal{E},\sigma\subset\partial\Omega\}$, ${\mathcal{E}}_{{\rm int}}=\mathcal{E}\setminus{\mathcal{E}}_{{\rm ext}}$, ${\mathcal{E}}_{{\rm int}}^{(i)}={\mathcal{E}}_{{\rm int}}\cap\mathcal{E}^{(i)}$ and ${\mathcal{E}}_{{\rm ext}}^{(i)}={\mathcal{E}}_{{\rm ext}}\cap\mathcal{E}^{(i)}$. The boundary conditions (1c) are taken into account by setting $u_{\sigma,i}=0$ for all $\sigma\in{\mathcal{E}}_{{\rm ext}}^{(i)}$, $1\leq i\leq d$. Let $\delta t>0$ be a constant time step. The approximate solution $(\rho^{n},{\boldsymbol{u}}^{n})$ at time $t_{n}=n\delta t$ for $1\leq n\leq N={\lfloor T/\delta t\rfloor}$ is computed as follows: knowing $\{\rho_{K}^{n-1},\rho_{K}^{n},K\in{\mathcal{M}}\}\subset\mathbb{R}$ and $(u_{\sigma,i}^{n})_{\sigma\in{\mathcal{E}}_{{\rm int}}^{(i)},1\leq i\leq d}\subset\mathbb{R}$, find $(\rho_{K}^{n+1})_{K\in{\mathcal{M}}}\subset\mathbb{R}$ and $(u_{\sigma,i}^{n+1})_{\sigma\in{\mathcal{E}}_{{\rm int}}^{(i)},1\leq i\leq d}\subset\mathbb{R}$ by the following algorithm:

[TABLE]

where the discrete densities and space operators are defined below (see also [3, 2]).

Mass convection flux – Given a discrete density field $\rho=\{\rho_{K},\,K\in{\mathcal{M}}\}$, and a velocity field ${\boldsymbol{u}}=\{u_{\sigma,i},\,\sigma\in\mathcal{E}^{(i)},\,1\leq i\leq d\}$, the convection term in (4d) reads:

[TABLE]

where $F_{K,\sigma}(\rho,{\boldsymbol{u}})$ stands for the mass flux across $\sigma$ outward $K$. This flux is set to 0 on external faces to account for the homogeneous Dirichlet boundary conditions; it is given on internal faces by:

[TABLE]

where $u_{K,\sigma}=u_{\sigma,i}\ {\boldsymbol{n}}_{K,\sigma}\cdot{\boldsymbol{e}}^{(i)}$, with ${\boldsymbol{e}}^{(i)}$ the $i$-th vector of the orthonormal basis of $\mathbb{R}^{d}$. The density at the face $\sigma=K|L$ is approximated by the upwind technique, i.e. $\rho_{\sigma}=\rho_{K}$ if $u_{K,\sigma}\geq 0$ and $\rho_{\sigma}=\rho_{L}$ otherwise.

Pressure gradient term – In (4a) and (4c), the term $({\boldsymbol{\nabla}}p)_{\sigma,i}$ stands for the $i^{th}$ component of the discrete pressure gradient at the face $\sigma$. Given a discrete density field $\rho=\{\rho_{K},\,K\in{\mathcal{M}}\}$, this term is defined as:

[TABLE]

Defining for all $K\in{\mathcal{M}}$, $({\rm div}{\boldsymbol{u}})_{K}={\rm div}(1\times{\boldsymbol{u}})_{K}$ (see (5)), the following discrete duality relation holds for all discrete density and velocity fields $(\rho,{\boldsymbol{u}})$:

[TABLE]

The MAC scheme is inf-sup stable: there exists $\beta>0$, depending only on $\Omega$ and the regularity of the mesh, such that, for all $p=\{p_{K},\,K\in{\mathcal{M}}\}$, there exists ${\boldsymbol{u}}=\{u_{\sigma,i},\,\sigma\in\mathcal{E}^{(i)},\,1\leq i\leq d\}$ satisfying homogeneous Dirichlet boundary conditions with:

[TABLE]

where $||{\boldsymbol{u}}||_{1,{\mathcal{M}}}$ is the usual discrete ${\rm H}^{1}$-norm of ${\boldsymbol{u}}$ (see [1]).

Velocity convection operator – Given a density field $\rho=\{\rho_{K},\,K\in{\mathcal{M}}\}$, and two velocity fields ${\boldsymbol{u}}=\{u_{\sigma,i},\,\sigma\in\mathcal{E}^{(i)},\,1\leq i\leq d\}$ and ${\boldsymbol{v}}=\{v_{\sigma,i},\,\sigma\in\mathcal{E}^{(i)},\,1\leq i\leq d\}$, we build for each $\sigma\in{\mathcal{E}}_{{\rm int}}$ the following quantities:

• •

an approximation of the density on the dual cell \rho_{\scalebox{0.6}{D_{\sigma}}} defined as:

[TABLE]

• •

a discrete divergence for the convection on the dual cell $D_{\sigma}$:

[TABLE]

For $i\in\{1,..,d\}$, and $\sigma\in{\mathcal{E}}_{{\rm int}}^{(i)}$, $\sigma=K|L$,

If the vector ${\boldsymbol{e}}^{(i)}$ is normal to $\varepsilon$, $\varepsilon$ is included in a primal cell $K$, and we denote by $\sigma^{\prime}$ the second face of $K$ which, in addition to $\sigma$, is normal to ${\boldsymbol{e}}^{(i)}$. We thus have $\varepsilon=D_{\sigma}|D_{\sigma^{\prime}}$. Then the mass flux through $\varepsilon$ is given by:

[TABLE]

• -

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

[TABLE]

With this definition, the dual fluxes are locally conservative through dual faces $\varepsilon=D_{\sigma}|D_{\sigma^{\prime}}$ (i.e. $F_{\sigma,\varepsilon}(\rho,{\boldsymbol{u}})=-F_{\sigma^{\prime},\varepsilon}(\rho,{\boldsymbol{u}})$), and vanish through a dual face included in the boundary of $\Omega$. For this reason, the values $v_{\varepsilon,i}$ are only needed at the internal dual faces, and are chosen centered, i.e., for $\varepsilon=D_{\sigma}|D_{\sigma^{\prime}}$, $v_{\varepsilon,i}=(v_{\sigma,i}+v_{\sigma^{\prime},i})/2$.

As a result, a finite volume discretization of the mass balance (1a) holds over the internal dual cells. Indeed, if $\rho^{n+1}=\{\rho_{K}^{n+1},\,K\in{\mathcal{M}}\}$, $\rho^{n}=\{\rho_{K}^{n},\,K\in{\mathcal{M}}\}$ and ${\boldsymbol{u}}^{n+1}=\{u_{\sigma,i}^{n+1},\,\sigma\in\mathcal{E}^{(i)},\,1\leq i\leq d\}$ are density and velocity fields satisfying (4d), then, the dual quantities \{\rho_{\scalebox{0.6}{D_{\sigma}}}^{n+1},\rho_{\scalebox{0.6}{D_{\sigma}}}^{n},\,\sigma\in{\mathcal{E}}_{{\rm int}}\} and the dual fluxes $\{F_{\sigma,\varepsilon}(\rho^{n+1},{\boldsymbol{u}}^{n+1}),\,\sigma\in{\mathcal{E}}_{{\rm int}},\,\varepsilon\in\bar{\mathcal{E}}(D_{\sigma})\}$ satisfy a finite volume discretization of the mass balance (1a) over the internal dual cells:

[TABLE]

Diffusion term – The discrete diffusion term in (4b) is defined in [2] and is coercive in the following sense: for every discrete velocity field ${\boldsymbol{u}}$ satisfying the homogeneous Dirichlet boundary conditions, one has:

[TABLE]

The initialization of the scheme (4) is performed by setting

[TABLE]

and computing $\rho^{-1}$ by solving the backward mass balance equation (4d) for $n=-1$ where the unknown is $\rho^{-1}$ and not $\rho^{0}$. This allows to perform the first prediction step with \{\rho_{\scalebox{0.6}{D_{\sigma}}}^{0},\rho_{\scalebox{0.6}{D_{\sigma}}}^{-1},\,\sigma\in{\mathcal{E}}_{{\rm int}}\} and the dual mass fluxes $\{F_{\sigma,\varepsilon}(\rho^{0},{\boldsymbol{u}}^{0}),\,\sigma\in{\mathcal{E}}_{{\rm int}},\,\varepsilon\in\bar{\mathcal{E}}(D_{\sigma})\}$ satisfying the mass balance (12). Moreover, since $\rho_{0}^{\varepsilon}>0$, one clearly has $\rho_{K}^{0}>0$ for all $K\in{\mathcal{M}}$ and therefore \rho_{\scalebox{0.6}{D_{\sigma}}}^{0}>0 for all $\sigma\in{\mathcal{E}}_{{\rm int}}$. The positivity of $\rho^{-1}$ is a consequence of the following Lemma.

Lemma 2.1**.**

If $(\rho_{0}^{\varepsilon},{\boldsymbol{u}}_{0}^{\varepsilon})$ satisfies (3), then there exists $C$, depending on the mesh but independent of $\varepsilon$ such that:

[TABLE]

Proof.

We sketch the proof. The boundedness of the first two terms is a straightforward consequence of (3). For the third term we remark that, again by (3):

[TABLE]

∎

3. Asymptotic analysis of the zero Mach limit

By the results of [3], there exists a solution $(\rho^{n},{\boldsymbol{u}}^{n})_{0\leq n\leq N}$ to the scheme (4) and any solution satisfies the following relations:

•  a discrete kinetic energy balance:  for all $\sigma\in{\mathcal{E}}_{{\rm int}}^{(i)}$, $1\leq i\leq d$, $0\leq n\leq N-1$:

[TABLE]

•  a discrete renormalization identity:  for all $K\in{\mathcal{M}}$, $0\leq n\leq N-1$:

[TABLE]

with $R_{K}^{n+1}\geq 0$, where the function $b_{\gamma}$ is defined by $b_{\gamma}(\rho)=\rho\log\rho$ if $\gamma=1$, $b_{\gamma}(\rho)=\rho^{\gamma}/(\gamma-1)$ if $\gamma>1$ and satisfies $\rho b_{\gamma}^{\prime}(\rho)-b_{\gamma}(\rho)=\rho^{\gamma}=\wp(\rho)$ for all $\rho>0$, and $\Pi_{\gamma}(\rho)=b_{\gamma}(\rho)-b_{\gamma}(1)-b_{\gamma}^{\prime}(1)(\rho-1)$.

Summing (15) and (16) over the primal cells from one side, and over the dual cells and the components on the other side, and invoking the grad-div duality relation (8), we obtain a local-in-time discrete entropy inequality, for $0\leq n\leq N-1$:

[TABLE]

where $\displaystyle\mathcal{R}^{n+1}=\sum_{i=1}^{d}\sum_{\sigma\in{\mathcal{E}}_{{\rm int}}^{(i)}}R^{n+1}_{\sigma,i}+\frac{1}{\varepsilon^{2}}\sum_{K\in{\mathcal{M}}}R_{K}^{n+1}\geq 0$.

The function $\Pi_{\gamma}$ has some important properties:

[TABLE]

Lemma 3.1** (Global discrete entropy inequality).**

Under assumption (3), there exists $C_{0}>0$ independent of $\varepsilon$ such that the solution $(\rho^{n},{\boldsymbol{u}}^{n})_{0\leq n\leq N}$ to the scheme (4) satisfies, for $\varepsilon$ small enough, and for $1\leq n\leq N$:

[TABLE]

Proof.

Summing (17) over $n$ yields the inequality (19) with

[TABLE]

By (14), for $\varepsilon$ small enough, one has $\rho_{K}^{-1}\leq 2$ for all $K\in{\mathcal{M}}$ and therefore \rho_{\scalebox{0.6}{D_{\sigma}}}^{-1}\leq 2 for all $\sigma\in{\mathcal{E}}_{{\rm int}}^{(i)}$ and $1\leq i\leq d$. Hence, since ${\boldsymbol{u}}_{0}^{\varepsilon}$ is uniformly bounded in ${\rm H}^{1}(\Omega)^{d}$ by (3), a classical trace inequality yields the boundedness of the first term. Again by (14), one has $|\rho_{K}^{0}-1|\leq C\varepsilon^{2}$ for all $K\in{\mathcal{M}}$. Hence, by (18a), the second term vanishes as $\varepsilon\to 0$. The third term is also uniformly bounded with respect to $\varepsilon$ thanks to (14). ∎

Lemma 3.2** (Control of the pressure).**

Assume that $(\rho_{0}^{\varepsilon},{\boldsymbol{u}}_{0}^{\varepsilon})$ satisfies (3) and let $(\rho^{n},{\boldsymbol{u}}^{n})_{0\leq n\leq N}$ satisfy (4). Let $p^{n}=\wp(\rho^{n})$ and define ${\delta p}^{n}=\{{\delta p}_{K}^{n},\,K\in{\mathcal{M}}\}$ where ${\delta p}_{K}^{n}=(p_{K}^{n}-|\Omega|^{-1}\int_{\Omega}p^{n}\,\mathrm{d}{\boldsymbol{x}})/\varepsilon^{2}$. Then, one has, for all $1\leq n\leq N$:

[TABLE]

where $C_{{\mathcal{M}},\delta t}\geq 0$ depends on the mesh and $\delta t$ but not on $\varepsilon$, and $||\cdot||$ stands for any norm on the space of discrete functions.

Proof.

By (19), the discrete pressure gradient is controlled in ${\rm L}^{\infty}$ by $C_{{\mathcal{M}},\delta t}\,\varepsilon^{2}$, so that ${\boldsymbol{\nabla}}(\delta p^{n})$ is bounded in any norm independently of $\varepsilon$. Using the discrete $({\rm H}^{-1})^{d}$-norm (see e.g. [1]), invoking the gradient divergence duality (8) and the inf-sup stability of the scheme, $||{\boldsymbol{\nabla}}(\delta p^{n})||_{-1,{\mathcal{M}}}\leq C_{{\mathcal{M}},\delta t}$ implies that $||{\delta p}^{n}||_{{\rm L}^{2}}\leq\beta^{-1}C_{{\mathcal{M}},\delta t}$. ∎

Theorem 3.3** (Incompressible limit of the MAC pressure correction scheme).**

Let $(\varepsilon^{(m)})_{m\in\mathbb{N}}$ be a sequence of positive real numbers tending to zero, and let $(\rho^{(m)},{\boldsymbol{u}}^{(m)})$ be a corresponding sequence of solutions of the scheme (4). Then the sequence $(\rho^{(m)})_{m\in\mathbb{N}}$ converges to the constant function $\rho=1$ when $m$ tends to $+\infty$ in ${\rm L}^{\infty}((0,T),{\rm L}^{q}(\Omega))$, for all $q\in[1,\min(\gamma,2)]$.

In addition, the sequence $({\boldsymbol{u}}^{(m)},{\delta p}^{(m)})_{m\in\mathbb{N}}$ tends, in any discrete norm, to the solution $({\boldsymbol{u}},{\delta p})$ of the usual MAC pressure correction scheme for the incompressible Navier-Stokes equations, which reads:

[TABLE]

Proof.

By (18b) and the global entropy estimate (19), one has for $\gamma\geq 2$,

[TABLE]

For $1\leq\gamma\leq 2$, invoking (18g) and estimate (19), we obtain for all $t\in(0,T)$ and for all $R\in(2,+\infty)$:

[TABLE]

which proves the convergence of $(\rho^{(m)})_{m\in\mathbb{N}}$ to the constant function $\rho=1$ as $m\to+\infty$ in ${\rm L}^{\infty}((0,T),{\rm L}^{q}(\Omega))$ for all $q\in[1,\min(\gamma,2)]$. Using again (19), the sequence $({\boldsymbol{u}}^{(m)})_{m\in\mathbb{N}}$ is bounded in any discrete norm and the same holds for the sequence $({\delta p}^{(m)})_{m\in\mathbb{N}}$ by Lemma 3.2. By the Bolzano-Weiertrass theorem and a norm equivalence argument, there exists a subsequence of $({\boldsymbol{u}}^{(m)},{\delta p}^{(m)})_{m\in\mathbb{N}}$ which tends, in any discrete norm, to a limit $({\boldsymbol{u}},{\delta p})$. Passing to the limit cell-by-cell in (4), one obtains that $({\boldsymbol{u}},{\delta p})$ is a solution to (21). Since this solution is unique, the whole sequence converges, which concludes the proof. ∎