Convergence of the mac scheme for variable density flows

TL;DR

This paper proves the convergence of a semi-implicit MAC scheme applied to the time-dependent variable density Navier-Stokes equations, ensuring numerical stability and accuracy for simulating variable density flows.

Contribution

It provides a rigorous proof of convergence for the MAC scheme in the context of variable density flows, which was previously unestablished.

Findings

Convergence of the scheme is mathematically established.

The scheme is shown to be stable for variable density Navier-Stokes equations.

The results support reliable numerical simulations of variable density flows.

Abstract

We prove in this paper the convergence of an semi-implicit MAC scheme for the time-dependent variable density Navier-Stokes equations.

Equations

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

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

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

$$-{\int }_0^T{\int }_{\Omega }\bar{\rho }{\partial }_t\phi +\bar{\rho }\bar{\mathbf{u}}\cdot \nabla \phi \mathrm{d}\mathbf{x}\mathrm{d}t={\int }_{\Omega }{\rho }_0(\mathbf{x})\phi (\mathbf{x},0)\mathrm{d}\mathbf{x}.$$

$$\int_{0}^{T}\int_{\Omega}\bigr{[}-\bar{\rho}\,\bar{\boldsymbol{u}}\cdot\partial_{t}{\boldsymbol{v}}-(\bar{\rho}\,\bar{\boldsymbol{u}}\otimes\bar{\boldsymbol{u}}):\nabla{\boldsymbol{v}}+\nabla\bar{\boldsymbol{u}}:\nabla{\boldsymbol{v}}\bigl{]}\ \mathrm{d}\boldsymbol{x}\ \mathrm{d}t=\int_{\Omega}\rho_{0}{\boldsymbol{u}}_{0}\cdot{\boldsymbol{v}}(\cdot,0)\ \mathrm{d}\boldsymbol{x}+\int_{0}^{T}\int_{\Omega}{\boldsymbol{f}}\cdot{\boldsymbol{v}}\ \mathrm{d}\boldsymbol{x}\ \mathrm{d}t.$$

$$\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{\}},\quad\mbox{for }1\leq i\leq d$$

$$u_i=\sum _{n=0}^{N-1}\sum _{\sigma \in {\mathcal{E}}_{\mathrm{int}}^{(i)}}{u}_{\sigma }^{n+1}\text{mathds}{1}_{D_{\sigma }}\text{mathds}{1}_{]t_n,t_{n+1}]},$$

$$p=\sum _{n=0}^{N-1}\sum _{K\in \mathcal{M}}{p}_K^{n+1}\text{mathds}{1}_K\text{mathds}{1}_{]t_n,t_{n+1}]},\rho =\sum _{n=0}^{N-1}\sum _{K\in \mathcal{M}}{\rho }_K^{n+1}\text{mathds}{1}_K\text{mathds}{1}_{]t_n,t_{n+1}]},$$

$$\text{mbox}\mathbf{Initialization}:{\mathbf{u}}^{(0)}={\tilde{\mathcal{P}}}_{\mathcal{E}}{\mathbf{u}}_0,{\rho }^{(0)}={\mathcal{P}}_{\mathcal{M}}{\rho }_0.$$

$$\text{mbox}\mathbf{For}0\le n\le N-1,\text{mbox}solvefor{\mathbf{u}}^{n+1}\in {\mathbf{H}}_{\mathcal{E},0},{\rho }^{n+1}\in L_{\mathcal{M}}\text{mbox}and{p}^{n+1}\in L_{\mathcal{M},0}:$$

$${ð}_t{\rho }^{n+1}+{\mathrm{div}}_{\mathcal{M}}({\rho }^{n+1}{\mathbf{u}}^n)=0,$$

$${ð}_t(\rho \mathbf{u}{)}^{n+1}+{\mathbf{C}}_{\mathcal{E}}({\rho }^{n+1}{\mathbf{u}}^n){\mathbf{u}}^{n+1}-{\Delta }_{\mathcal{E}}{\mathbf{u}}^{n+1}+{\nabla }_{\mathcal{E}}{p}^{n+1}={\mathbf{f}}_{\mathcal{E}}^{n+1},$$

$${\mathrm{div}}_{\mathcal{M}}{\mathbf{u}}^{n+1}=0,$$

$$v_i⟼{\tilde{\mathcal{P}}}_{{\mathcal{E}}^{(i)}}{v}_i=\sum _{\sigma \in {\mathcal{E}}^{(i)}}{v}_{\sigma }\text{mathds}{1}_{D_{\sigma }}\text{mbox}with{v}_{\sigma }=\frac{1}{|\sigma |}{\int }_{\sigma }{v}_i\mathrm{d}\gamma (\mathbf{x}),\sigma \in {\mathcal{E}}^{(i)}.$$

$${ð}_t\rho (\mathbf{x},t)=\sum _{n=0}^{N-1}{ð}_t{\rho }^{n+1}(\mathbf{x})\text{mathds}{1}_{]t_n,t_{n+1}]}(t)\text{mbox}with{ð}_t{\rho }^{n+1}=\sum _{K\in \mathcal{M}}\frac{1}{\delta t}({\rho }_K^{n+1}-{\rho }_K^n)\text{mathds}{1}_K.$$

$${\mathrm{div}}_{\mathcal{M}}:L_{\mathcal{M}}\times {\mathbf{H}}_{\mathcal{E},0}⟶L_{\mathcal{M}},(\rho ,\mathbf{u})\mapsto {\mathrm{div}}_{\mathcal{M}}(\rho \mathbf{u})=\sum _{K\in \mathcal{M}}\frac{1}{|K|}\sum _{\sigma \in \mathcal{E}(K)}{F}_{K,\sigma }\text{mathds}{1}_K,$$

$${\nabla }_{\mathcal{E}}p=\sum _{\sigma =K|L\in {\mathcal{E}}_{\mathrm{int}}^{(i)}}(ðp{)}_{\sigma }{\mathbf{n}}_{K,\sigma }\text{mathds}{1}_{D_{\sigma }},\text{mbox}with(ðp{)}_{\sigma }=\frac{|\sigma |}{|D_{\sigma }|}(p_L-p_K).$$

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

$$F_{\sigma ,ϵ}=\frac{1}{2}(F_{K,\tau }+F_{L,{\tau }^{\prime }}).$$

$$\begin{array}[]{l|l}{\mathrm{div}}_{{\mathcal{E}}^{(i)}}(\rho,{\boldsymbol{v}}):&L_{\mathcal{M}}\times{\boldsymbol{H}}_{\!{\mathcal{E}},0}\longrightarrow L_{\mathcal{E}}\\ &\displaystyle\quad(\rho,{\boldsymbol{v}})\longmapsto{\mathrm{div}}_{{\mathcal{E}}^{(i)}}(\rho,{\boldsymbol{v}})=\sum_{{\sigma}\in{\mathcal{E}}_{\mathrm{int}}^{(i)}}{\mathrm{div}}_{D_{\sigma}}(\rho{\boldsymbol{v}})\ \mathds{1}_{D_{\sigma}},\\ &\quad\displaystyle\mbox{with }{\mathrm{div}}_{D_{\sigma}}(\rho,{\boldsymbol{v}})=\frac{1}{|D_{\sigma}|}\sum_{\epsilon\in{\widetilde{{\mathcal{E}}}}(D_{\sigma})}F_{{\sigma},\epsilon},\;\forall{\sigma}\in{\mathcal{E}}_{\mathrm{int}}^{(i)}.\end{array}$$

$$|D_{\sigma }|{\rho }_{D_{\sigma }}=|D_{K,\sigma }|{\rho }_K+|D_{L,\sigma }|{\rho }_L,\sigma \in {\mathcal{E}}_{\mathrm{int}},\sigma =K|L.$$

$$\frac{1}{\delta t}({\rho }_{D_{\sigma }}^{n+1}-{\rho }_{D_{\sigma }}^n)+{\mathrm{div}}_{D_{\sigma }}({\rho }^{n+1}{\mathbf{u}}^n)=0.$$

$$\forall \rho \in L_{\mathcal{M}},\forall \mathbf{v}\in {\mathbf{H}}_{\mathcal{E},0},\forall w\in H_{{\mathcal{E}}^{(i)},0},{\int }_{\Omega }{\mathrm{div}}_{{\mathcal{E}}^{(i)}}(\rho ,\mathbf{v})w\mathrm{d}\mathbf{x}={\int }_{\Omega }\rho \mathbf{v}\cdot {\nabla }_{{\mathcal{E}}^{(i)}}w\mathrm{d}\mathbf{x},$$

$$[(\rho \mathbf{v}{)}_{{\mathcal{E}}^{(i)}}{]}_j=\sum _{ϵ\in {\tilde{\mathcal{E}}}^{(i,j)}}{\rho }_{ϵ}{\tilde{v}}_{ϵ}\text{mathds}{1}_{D_{ϵ}},[(\nabla w{)}_{{\mathcal{E}}^{(i)}}{]}_j=\sum _{ϵ\in {\tilde{\mathcal{E}}}^{(i,j)},\sigma =\overrightarrow{\sigma |{\sigma }^{\prime }}}\frac{u_{{\sigma }^{\prime }}-u_{\sigma }}{d_{ϵ}}\text{mathds}{1}_{D_{ϵ}},$$

$$\begin{array}[]{l|l}C^{(i)}_{\mathcal{E}}(\rho,{\boldsymbol{u}}):&\quad H_{\!{\mathcal{E}}^{(i)},0}\longrightarrow H_{\!{\mathcal{E}}^{(i)},0}\\ &\displaystyle\quad v\longmapsto{C}^{(i)}_{\mathcal{E}}(\rho,{\boldsymbol{u}})v=\sum_{{\sigma}\in{{\widetilde{{\mathcal{E}}}}^{(i)}_{{\rm int}}}}\frac{1}{|D_{\sigma}|}\sum_{\begin{subarray}{c}\epsilon\in\tilde{{\mathcal{E}}}(D_{\sigma})\\ \epsilon={\sigma}|{\sigma}^{\prime}\end{subarray}}\!F_{{\sigma},\epsilon}\ \frac{v_{{\sigma}}+v_{{\sigma}^{\prime}}}{2}\ \mathds{1}_{D_{\sigma}}.\end{array}$$

$$\begin{array}[]{ll}&|{\boldsymbol{C}}_{\mathcal{E}}(\rho{\boldsymbol{u}}){\boldsymbol{v}}\cdot{\boldsymbol{w}}|\leq C_{\eta_{\mathcal{M}}}\,\|\rho\|_{L^{\infty}(\Omega)}\,\|{\boldsymbol{u}}\|_{L^{4}(\Omega)^{d}}\,\|{\boldsymbol{v}}\|_{L^{4}(\Omega)^{d}}\,\|{\boldsymbol{w}}\|_{1,{\mathcal{E}},0}\\[4.30554pt] \mbox{and}&|{\boldsymbol{C}}_{\mathcal{E}}(\rho{\boldsymbol{u}}){\boldsymbol{v}}\cdot{\boldsymbol{w}}|\leq C_{\eta_{\mathcal{M}}}\,\|\rho\|_{L^{\infty}(\Omega)}\,\|{\boldsymbol{u}}\|_{1,{\mathcal{E}},0}\,\|{\boldsymbol{v}}\|_{1,{\mathcal{E}},0}\,\|{\boldsymbol{w}}\|_{1,{\mathcal{E}},0}.\end{array}$$

$${\rho }_{\min }\le {\rho }^{n+1}\le {\rho }_{\max },$$

$$\frac{1}{2\delta t}\bigl{[}\rho_{D_{{\sigma}}}^{n+1}(u_{{\sigma}}^{n+1})^{2}-\rho_{D_{{\sigma}}}^{n}(u_{{\sigma}}^{n})^{2}\bigr{]}+\frac{1}{2|D_{{\sigma}}|}\sum_{\begin{subarray}{c}\epsilon\in\tilde{{\mathcal{E}}}(D_{{\sigma}})\\ \epsilon={\sigma}|{\sigma}^{\prime}\end{subarray}}F_{{\sigma},\epsilon}(\rho^{n+1},u^{n})u_{{\sigma}}^{n+1}u_{{\sigma}^{\prime}}^{n+1}\\ -(\Delta u)_{{\sigma}}^{n+1}u_{{\sigma}}^{n+1}+(\eth p)_{\sigma}^{n+1}u_{{\sigma}}^{n+1}-f_{{\sigma}}^{n+1}u_{{\sigma}}^{n+1}=-\dfrac{1}{2\delta t}\rho_{D_{{\sigma}}}^{n}\bigl{(}u_{{\sigma}}^{n+1}-u_{{\sigma}}^{n}\bigr{)}^{2}.$$

$$\Vert \mathbf{u}{\Vert }_{L^2({\mathbf{H}}_{\mathcal{E},0})}=\sum _{n=0}^{N-1}\delta t\Vert {\mathbf{u}}^{n+1}{\Vert }_{1,\mathcal{E},0}^2\le C\text{mbox}and\Vert \mathbf{u}{\Vert }_{L^{\infty }(L^2)}=\underset{0\le n\le N-1}{\max }\Vert {\mathbf{u}}^{n+1}{\Vert }_{L^2(\Omega {)}^d}\le C.$$

$${\int }_0^{T-\tau }{\int }_{\Omega }|\mathbf{u}(\mathbf{x},t+\tau )-\mathbf{u}(\mathbf{x},t){|}^2\mathrm{d}\mathbf{x}\mathrm{d}t\le C_{{\eta }_{\mathcal{M}},T}\frac{{\rho }_{\max }}{{\rho }_{\min }}(\Vert \mathbf{u}{\Vert }_{L^2({\mathbf{H}}_{\mathcal{E},0})}^3+1)\sqrt{\tau +\delta t}$$

$$A^{(i)}={\int }_0^{T-\tau }(A_1^{(i)}(t)+A_2^{(i)}(t))\mathrm{d}t,\text{mbox}with$$

$$A_1^{(i)}(t)={\int }_{\Omega }\left(\tilde{\rho }(\mathbf{x},t+\tau )u(\mathbf{x},t+\tau )-\tilde{\rho }(\mathbf{x},t)u(\mathbf{x},t)\right)w(\mathbf{x},t)\mathrm{d}\mathbf{x},$$

$$A_2^{(i)}(t)={\int }_{\Omega }\left(\tilde{\rho }(\mathbf{x},t)-\tilde{\rho }(\mathbf{x},t+\tau )\right)u(\mathbf{x},t+\tau )w(\mathbf{x},t)\mathrm{d}\mathbf{x}.$$

Full text

1. Introduction

Let $\Omega$ be a parallelepiped of $\mathbb{R}^{d}$, with $d\in\{2,3\}$ and $T>0$, and consider the following variable density Navier-Stokes equations posed on $\Omega\times(0,T)$:

[TABLE]

where $\bar{\rho}$, $\bar{\boldsymbol{u}}$ and $\bar{p}$ are the density, the velocity and the pressure of the flow and ${\boldsymbol{f}}\in L^{2}(0,T;L^{2}(\Omega)^{d})$. This system is complemented with initial and boundary conditions $\bar{\boldsymbol{u}}|_{\partial{\Omega}}=0,$ $\bar{\boldsymbol{u}}|_{t=0}={\boldsymbol{u}}_{0}$, $\bar{\rho}|_{t=0}=\rho_{0}$, which are such that $\rho_{0}\in L^{\infty}(\Omega)$, $0<\rho_{\min}<\rho_{0}\leq\rho_{\max}$ and ${\boldsymbol{u}}_{0}\in L^{2}(\Omega)^{d}$. A pair $(\bar{\rho},\bar{\boldsymbol{u}})$ is a weak solution of problem (1) if it satisfies the following properties:

• –

$\bar{\rho}\in\{\rho\in L^{\infty}(\Omega\times(0,T)),\ \rho>0\ a.e.\ \mbox{in}\ \Omega\times(0,T)\}$.

• –

$\bar{\boldsymbol{u}}\in\{{\boldsymbol{u}}\in L^{\infty}(0,T;L^{2}(\Omega)^{d})\cap L^{2}(0,T;H^{1}_{0}(\Omega)^{d}),\ {\mathrm{div}}\ {\boldsymbol{u}}=0\ a.e.\ \mbox{in}\ \Omega\times(0,T)\}$.

• –

For all $\varphi$ in $C_{c}^{\infty}(\Omega\times[0,T))$,

[TABLE]

• –

For all ${\boldsymbol{v}}$ in $\displaystyle\{{\boldsymbol{w}}\in C_{c}^{\infty}(\Omega\times[0,T))^{d},{\mathrm{div}}\ {\boldsymbol{w}}=0\}$,

[TABLE]

The existence of such a weak solution was proven in [9]; convergence results exist for the discontinuous Galerkin approximation [8] and for a finite volume/finite element scheme [7]. Here we prove the convergence of the MAC scheme.

2. The numerical scheme

Let ${\mathcal{M}}$ be a MAC mesh (see e.g. [4] and Figure 1 for the notations). The discrete pressure and density unknowns are associated with the cells of the mesh ${\mathcal{M}}$, and are denoted by \big{\{}\rho_{K},\ K\in{\mathcal{M}}\big{\}} and \big{\{}p_{K},\ K\in{\mathcal{M}}\big{\}}. The discrete velocity unknowns approximate the normal velocity to the mesh faces, and are denoted $(u_{\sigma})_{{\sigma}\in{{\mathcal{E}}^{(i)}}}$, i\in\bigl{[}|1,d|\bigr{]}, where ${\mathcal{E}}$ is the set of the faces of the mesh, and ${\mathcal{E}}^{(i)}$ the subset of the faces orthogonal to the $i$-th vector of the canonical basis of $\mathbb{R}^{d}$. We define ${\mathcal{E}}_{\mathrm{ext}}=\{{\sigma}\in{\mathcal{E}},{\sigma}\subset\partial\Omega\}$, ${\mathcal{E}}_{\mathrm{int}}={\mathcal{E}}\setminus{\mathcal{E}}_{\mathrm{ext}}$, ${\mathcal{E}}_{\mathrm{int}}^{(i)}={\mathcal{E}}_{\mathrm{int}}\cap{\mathcal{E}}^{(i)}$ and ${\mathcal{E}}_{\mathrm{ext}}^{(i)}={\mathcal{E}}_{\mathrm{ext}}\cap{\mathcal{E}}^{(i)}$.

The regularity of the mesh is defined by:

[TABLE]

and we denote by $h_{\mathcal{M}}$ the space step. The discrete space $L_{\mathcal{M}}$ for the scalar unknowns (i.e. the pressure and the density) is defined as the set of piecewise constant functions over each of the grid cells $K$ of ${\mathcal{M}}$, and the discrete space for the $i^{th}$ velocity component, $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:

[TABLE]

(i.e. we impose $u_{\sigma}=0$ for all ${\sigma}\in{\mathcal{E}}_{\mathrm{ext}}$). We then set ${\boldsymbol{H}}_{\!{\mathcal{E}},0}=\prod_{i=1}^{d}H_{\!{\mathcal{E}}^{(i)},0}$.

Let $0=t_{0}<t_{1}<\cdots<t_{N}=T$ be a partition of the time interval $(0,T)$, with $\delta t=t_{n+1}-t_{n}$. Let $\{u_{{\sigma}}^{n+1},{\sigma}\in{{\mathcal{E}}^{(i)}},0\leq n\leq N-1,\ 1\leq i\leq d\}$, $\{p_{K}^{n+1},K\in{\mathcal{M}},\ 0\leq n\leq N-1\}$ and $\{\rho_{K}^{n+1},K\in{\mathcal{M}},0\leq n\leq N-1\}$ be the sets of discrete velocity, pressure and density unknowns. Defining the characteristic function $\mathds{1}_{A}$ of any subset $A\subset\Omega$ by $\mathds{1}_{A}(\boldsymbol{x})=1$ if $\boldsymbol{x}\in A$ and $\mathds{1}_{A}(\boldsymbol{x})=0$ otherwise, the corresponding piecewise constant functions for the velocities are of the form:

[TABLE]

and $X_{i,{\mathcal{E}},\delta t}$ denotes the set of such piecewise constant functions on time intervals and dual cells; we then set $\boldsymbol{X}_{{\mathcal{E}},\delta t}=\prod_{i=1}^{d}X_{i,{\mathcal{E}},\delta t}$. The pressure and density discrete functions are defined by:

[TABLE]

and $Y_{{\mathcal{M}},\delta t}$ denotes the space of such piecewise constant functions. The numerical scheme reads:

[TABLE]

with the interpolators and discrete operators defined as follows.

Grid interpolators – The Fortin interpolator is defined by $\widetilde{\mathcal{P}}_{\mathcal{E}}{\boldsymbol{u}}=(\widetilde{\mathcal{P}}_{{{\mathcal{E}}^{(i)}}})_{i=1,\ldots,d}$ with $\widetilde{\mathcal{P}}_{{\mathcal{E}}^{(i)}}:H^{1}_{0}(\Omega)\longrightarrow H_{\!{\mathcal{E}}^{(i)},0}$ and

[TABLE]

For $q\in L^{2}(\Omega)$, $\mathcal{P}_{\mathcal{M}}q\in L_{\mathcal{M}}$ is defined by $\displaystyle\mathcal{P}_{\mathcal{M}}q(\boldsymbol{x})=\frac{1}{|K|}\int_{K}q\ \mathrm{d}\boldsymbol{x}$ for $\boldsymbol{x}\in K$.

Discrete time derivative – For $\rho\in Y_{{\mathcal{M}},\delta t}$, $\eth_{t}\rho\in Y_{{\mathcal{M}},\delta t}$ is defined by:

[TABLE]

Discrete divergence – Let $u_{K,\sigma}$ be defined as $u_{K,\sigma}=u_{\sigma}\,\boldsymbol{n}_{K,{\sigma}}\cdot{\boldsymbol{e}}_{i}$ for any face ${\sigma}\in{{\mathcal{E}}^{(i)}}$, $i=1,\ldots,d$. The discrete (upwind finite volume) divergence operator ${\mathrm{div}}_{\mathcal{M}}$ is defined by:

[TABLE]

with $F_{K,\sigma}=|{\sigma}|\ \rho_{\sigma}u_{K,\sigma}$ for $K\in{\mathcal{M}},\ {\sigma}=K|L\in{\mathcal{E}}(K)$, and $\rho_{\sigma}=\rho_{K}$ if $u_{K,{\sigma}}\geq 0$, $\rho_{\sigma}=\rho_{L}$ otherwise. For all $K\in{\mathcal{M}}$, we set $({\mathrm{div}}{\boldsymbol{u}})_{K}={\mathrm{div}}(1\times{\boldsymbol{u}})_{K}$.

Pressure gradient operator – The discrete pressure gradient operator is defined as the transpose of the divergence operator, so $\nabla_{\mathcal{E}}:\ L_{\mathcal{M}}\longrightarrow{\boldsymbol{H}}_{\!{\mathcal{E}},0},\ p\mapsto\nabla_{\mathcal{E}}(p)$ with:

[TABLE]

Discrete Laplace operator – The discrete diffusion operator $\Delta_{\mathcal{E}}$ is defined in [4] and is coercive in the sense that $-\int_{\Omega}\Delta_{{\mathcal{E}}}{\boldsymbol{v}}\cdot{\boldsymbol{v}}\ \mathrm{d}\boldsymbol{x}=\|{\boldsymbol{v}}\|_{1,{\mathcal{E}},0}^{2}$ for any ${\boldsymbol{v}}\in{\boldsymbol{H}}_{\!{\mathcal{E}}_{m},0}$, where $\|\cdot\|_{1,{\mathcal{E}},0}$ is the usual discrete ${\rm H}^{1}$-norm of ${\boldsymbol{u}}$ (see [4]). This inner product may also be formulated as the $L^{2}$-inner product of adequately chosen discrete gradients [4].

Discrete convection operator – The numerical convection fluxes and the approximations of $\rho$ in the momentum equation are chosen so as ensure that a discrete mass balance holds on the dual cells, in order to recover a discrete kinetic energy inequality. This idea was first introduced in [3, 1] for the Crouzeix-Raviart and Rannacher-Turek scheme, in [6] for the MAC scheme and was adapted to a DDFV scheme [5]. For $\epsilon={\sigma}|{\sigma}^{\prime}$, the convection flux $\int_{\epsilon}\rho u_{i}{\boldsymbol{u}}\cdot\boldsymbol{n}_{{\sigma},\epsilon}\ \mathrm{d}\gamma(\boldsymbol{x})$ is approximated by $F_{{\sigma},\epsilon}u_{\epsilon},$ where $u_{\epsilon}=(u_{\sigma}+u_{{\sigma}^{\prime}})/2$ and $F_{{\sigma},\epsilon}$ is the numerical mass flux through $\epsilon$ outward $D_{{\sigma}}$ defined as follows:

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

[TABLE]

• -

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$ with ${\tau}\in{\mathcal{E}}(K)$ and ${\tau}^{\prime}\in{\mathcal{E}}(L)$. Then:

[TABLE]

Remark 2.1**.**

In both cases, for $\epsilon={\sigma}|{\sigma}^{\prime}$, the mass flux $F_{{\sigma},\epsilon}$ may be written as $F_{{\sigma},\epsilon}=|\epsilon|\rho_{\epsilon}\tilde{u}_{\epsilon}$, with $\rho_{\epsilon}=(\rho_{\sigma}+\rho_{{\sigma}^{\prime}})$ and $\tilde{u}_{\epsilon}=(\rho_{\sigma}u_{\sigma}+\rho_{{\sigma}^{\prime}}u_{{\sigma}^{\prime}})/(\rho_{\sigma}+\rho_{{\sigma}^{\prime}})$ in the first case, and $\rho_{\epsilon}=(|{\tau}|\rho_{\tau}+|{\tau}^{\prime}|\rho_{{\tau}^{\prime}})/(|{\tau}|+|{\tau}^{\prime}|)$ and $\tilde{u}_{\epsilon}=(|{\tau}|\rho_{\tau}u_{\tau}+|{\tau}^{\prime}|\rho_{{\tau}^{\prime}}u_{{\tau}^{\prime}})/(|{\tau}|\rho_{\tau}+|{\tau}^{\prime}|\rho_{{\tau}^{\prime}})$ in the second case.

With this expression of the flux, we may define a discrete divergence operator on the dual cells:

[TABLE]

For the definition of the time-derivative $\eth_{t}(\rho{\boldsymbol{u}})$, an approximation of the density on the dual cell $\rho_{D_{\sigma}}$ is defined as:

[TABLE]

With the above definitions, if $(\rho,{\boldsymbol{u}})\in L_{\mathcal{M}}\times\boldsymbol{X}_{{\mathcal{E}},\delta t}$ satisfies the mass balance equation (4b), then the following mass balance on the dual cells holds:

[TABLE]

Note that a discrete duality property also holds, in the sense that, for $1\leq i\leq d$,

[TABLE]

where $(\rho{\boldsymbol{v}})_{{\mathcal{E}}^{(i)}}$ and $\nabla_{\!\!\!{{\mathcal{E}}^{(i)}}}w$ are vector valued functions of components:

[TABLE]

with $\rho_{\epsilon}$ and $\tilde{v}_{\epsilon}$ defined in Remark 2.1 and ${{\widetilde{{\mathcal{E}}}}^{(i,j)}}=\{\epsilon\in{{\widetilde{{\mathcal{E}}}}^{(i)}};\epsilon\perp{\boldsymbol{e}}^{(j)}\}$. We finally define the $i$-th component ${C}^{(i)}_{\mathcal{E}}(\rho{\boldsymbol{u}})$ of the non linear convection operator by:

[TABLE]

and the full (i.e. for all the velocity components) discrete convection operator ${{\boldsymbol{C}}}_{\mathcal{E}}(\rho,{\boldsymbol{u}}),\ {\boldsymbol{H}}_{\!{\mathcal{E}},0}\longrightarrow{\boldsymbol{H}}_{\!{\mathcal{E}},0}$ by ${{\boldsymbol{C}}}_{\mathcal{E}}(\rho,{\boldsymbol{u}}){\boldsymbol{v}}=(C^{(1)}_{\mathcal{E}}(\rho,{\boldsymbol{u}})v_{1},\ldots,C^{(d)}_{\mathcal{E}}(\rho{\boldsymbol{u}})v_{d})^{t}.$ Let ${\boldsymbol{E}}_{\mathcal{E}}$ be the subspace of ${\boldsymbol{H}}_{\!{\mathcal{E}},0}$ of divergence-free functions (with respect to the discrete divergence operator). By Hölder’s inequality and [4, Lemma 3.9], there exists $C_{\eta_{\mathcal{M}}}>0$ (depending only on $\eta_{\mathcal{M}}$) such that, $\forall\ (\rho,{\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}})\in L_{\mathcal{M}}\times{\boldsymbol{E}}_{\mathcal{E}}\times{\boldsymbol{H}}_{\!{\mathcal{E}},0}^{2}$,

[TABLE]

3. Estimates and convergence analysis

Since the velocity is divergence-free, the mass equation is a transport equation on $\rho$, so that, thanks to the upwind choice, the following estimate holds:

[TABLE]

and the $L^{2}$-norm of $\rho^{n+1}$ is lower than the $L^{2}$-norm of the initial data $\rho_{0}$, for $0\leq n\leq N-1$. In addition, thanks to (6), any solution to the scheme (4) satisfies the following discrete kinetic energy balance, for $1\leq i\leq d$, ${\sigma}\in{{\mathcal{E}}^{(i)}}$, $0\leq n\leq N-1,$

[TABLE]

From this inequality, we obtain estimates on the velocity. For ${\boldsymbol{u}}\in\boldsymbol{X}_{{\mathcal{E}},\delta t}$ satisfying (4), there exists $C>0$ depending on ${{\boldsymbol{u}}_{0}}$, $\rho_{0}$ and ${\boldsymbol{f}}$ such that,

[TABLE]

These estimates yields the existence of a unique solution to the scheme: indeed, the first equation may be solved separately for $\rho^{n+1}$ and is linear with repect to this unknown and, once $\rho^{n+1}$ is known, the last two equations are a linear generalized Oseen problem for ${\boldsymbol{u}}^{n+1}$ and $p^{n+1}$, which is uniquely solvable thanks to the inf-sup stability of the MAC discretization. The convergence of the scheme requires some time compactness. Contrary to the constant density case [4], there is no uniform estimate on the time derivative, and compactness is obtained thanks to the following lemma together with the Fréchet-Kolmogorov theorem.

Lemma 3.1** (Estimate on the time translates of the velocity).**

Let ${\boldsymbol{u}}\in X_{{\mathcal{E}},\delta t}$ and $\rho\in Y_{{\mathcal{M}},\delta t}$ and let $\tau>0$ then

[TABLE]

where $C_{\eta_{\mathcal{M}},T}>0$ only depends on $\Omega$, $T$, ${\boldsymbol{f}}$ and on the regularity of the mesh $\eta_{\mathcal{M}}$.

Proof.

In the continuous case, see e.g. [2, pages 444-452], the estimate (11) is obtained by bounding the term $\int_{0}^{T-\tau}\int_{\Omega}\left(\rho(\boldsymbol{x},t){\boldsymbol{u}}(\boldsymbol{x},t+\tau)-\rho(\boldsymbol{x},t){\boldsymbol{u}}(\boldsymbol{x},t)\right)\cdot{\boldsymbol{w}}(\boldsymbol{x},t)\ \mathrm{d}t$ with ${\boldsymbol{w}}(\boldsymbol{x},t)={\boldsymbol{u}}(\boldsymbol{x},t+\tau)-{\boldsymbol{u}}(\boldsymbol{x},t)$. However, in the context of the MAC scheme, the components of ${\boldsymbol{u}}$ are piecewise constant on different meshes so we need to treat the space indices separately. For a given $i=1,\ldots,d$, we denote by $u$ and $w$ the $i$-th component of ${\boldsymbol{u}}$ and ${\boldsymbol{w}}$, and by $\widetilde{\rho}$ the piecewise constant function defined by $\widetilde{\rho}(\boldsymbol{x},t)=\rho_{D_{\sigma}}^{n+1}$ for $(\boldsymbol{x},t)\in D_{\sigma}\times[t_{n}t_{n+1})$. We then wish to bound the terms

[TABLE]

For lack of space, we only deal here with the term $A_{2}^{(i)}(t)$. Thanks to the mass balance on the dual cells (6) and to the discrete duality formula (7) we have:

[TABLE]

Using Hölder’s inequalities and the fact that $\displaystyle\sum_{n=1}^{N-1}\delta t\ \mathds{1}_{(t,t+\tau)}(t_{n})\leq\tau+\delta t$,

[TABLE]

Now, by Hölder’s inequality,

[TABLE]

Therefore, integrating over $(0,T-\tau)$ yields that

[TABLE]

Similar computations for the term $\int_{0}^{T-\tau}A_{1}^{(i)}(t)\ \mathrm{d}t$ yield the result. ∎

Theorem 3.2** (Convergence of the scheme).**

Let $(\delta t_{m})_{m\in\mathbb{N}}$ and $({\mathcal{M}}_{m})_{m\in\mathbb{N}}$ be a sequence of time steps and MAC grids 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\mathbb{N}$. Let $(\rho_{m},{\boldsymbol{u}}_{m})$ be a solution to (4) for $\delta t=\delta t_{m}$ and ${\mathcal{M}}={\mathcal{M}}_{m}$. Then there exists $\bar{\rho}$ with $\rho_{\min}\leq\bar{\rho}\leq\rho_{\max}$ and $\bar{\boldsymbol{u}}\in L^{2}(0,T;\boldsymbol{E}(\Omega))$ such that, up to a subsequence:

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

• -

the sequence $(\rho_{m})_{m\in\mathbb{N}}$ converges to $\bar{\rho}$ in $\in L^{2}(0,T;L^{2}(\Omega))$,

• -

$(\bar{\rho},\bar{\boldsymbol{u}})$* is a solution to the weak formulation (2) and (3).*

Sketch of proof:

• –

Thanks to (8), there exists a subsequence of $(\rho_{m})_{m\in\mathbb{N}}$ star-weakly converging to some $\bar{\rho}$ in $L^{\infty}(\Omega\times(0,T))$; thanks to (10) and (11), there exists a subsequence of $(u_{m})_{m\in\mathbb{N}}$ converging to some $\bar{u}$ in $L^{2}(0,T;(L^{2}(\Omega)^{d})$.

• –

Passing to the limit in (4b) yields that $(\bar{\rho},\bar{u})$ satisfies (2).

• –

The strong convergence of the approximate densities is then obtained thanks to the $L^{2}$ estimates for $\rho$ in both the discrete and continuous case [7, Proposition 8.7].

• –

Passing to the limit in (4c) yields that $(\bar{\rho},\bar{u})$ satisfies (3).

• –

We finally obtain that $\bar{\boldsymbol{u}}\in L^{2}(0,T;{\boldsymbol{E}}(\Omega))$, where $E(\Omega)=\{{\boldsymbol{v}}\in H^{1}_{0}(\Omega)$ s.t. ${\mathrm{div}}{\boldsymbol{v}}=0\}$, as in [4, Proof of Theorem 4.3].