# Discretization of the viscous dissipation term with the MAC scheme

Fabrice Babik, Raphaele Herbin, Walid Kheriji, Jean-Claude Latché

► To cite this version:

Fabrice Babik, Raphaele Herbin, Walid Kheriji, Jean-Claude Latché. Discretization of the viscous dissipation term with the MAC scheme. FVCA 6, Internal Symposium, Jun 2011, Prague, Czech Republic. pp.571-579, 10.1007/978-3-642-20671-9_60. hal-00576268

HAL Id: hal-00576268

https://hal.science/hal-00576268v1

Submitted on 13 Aug 2011

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Discretization of the viscous dissipation term with the MAC scheme

F. Babik*, R. Herbin†, W. Kheriji*, J.-C. Latché*

*Institut de Radioprotection et de Sûreté Nucléaire (IRSN)

†Université de Provence

[fabrice.babik, walid.kheriji, jean-claude.latche]@irsn.fr, herbin@cmi.univ-mrs.fr

## Abstract

We propose a discretization for the MAC scheme of the viscous dissipation term $\tau(\mathbf{u}) : \nabla \mathbf{u}$ (where $\tau(\mathbf{u})$ stands for the shear stress tensor associated to the velocity field $\mathbf{u}$), which is suitable for the approximation of this term in a conservation equation for a scalar variable. This discretization enjoys the property that the integral over the computational domain $\Omega$ of the (discrete) dissipation term is equal to what is obtained when taking the inner product of the (discrete) momentum balance equation by $\mathbf{u}$ and integrating over $\Omega$. As a consequence, it may be used as an ingredient to obtain an unconditionally stable scheme for the compressible Navier-Stokes equations. It is also shown, in some model cases, to ensure the strong convergence in $L^1$ of the dissipation term.

**Key words** : MAC scheme, compressible Navier-Stokes equations, RANS turbulence models.

## 1 Introduction

Let us consider the compressible Navier-Stokes equations, which may be written as:

$$ \partial_t \rho + \text{div}(\rho \mathbf{u}) = 0, \qquad (1a) $$

$$ \partial_t(\rho \mathbf{u}) + \text{div}(\rho \mathbf{u} \otimes \mathbf{u}) + \nabla p - \text{div}(\tau(\mathbf{u})) = 0, \qquad (1b) $$

$$ \partial_t(\rho e) + \text{div}(\rho e \mathbf{u}) + p \text{div} \mathbf{u} + \text{div}(\mathbf{q}) = \tau(\mathbf{u}) : \nabla \mathbf{u}, \qquad (1c) $$

$$ \rho = \varphi(p, e), \qquad (1d) $$

where $t$ stands for the time, $\rho$, $\mathbf{u}$, $p$ and $e$ are the density, velocity, pressure and internal energy in the flow, $\tau(\mathbf{u})$ stands for the shear stress tensor, $\mathbf{q}$ for the energy diffusion flux, and the function $\varphi$ is the equation of state. This system of equations is posed over $\Omega \times (0, T)$, where $\Omega$ is a domain of $\mathbb{R}^d$, $d \le 3$. It must be supplemented by a closure relation for $\tau(\mathbf{u})$ and for $\mathbf{q}$, assumed to be:

$$ \tau(\mathbf{u}) = \mu(\nabla \mathbf{u} + \nabla^t \mathbf{u}) - \frac{2\mu}{3} \text{div} \mathbf{u} I, \quad \mathbf{q} = -\lambda \nabla e, \qquad (2) $$
where $\mu$ and $\lambda$ stand for two (possibly depending on $x$) positive parameters.

Let us suppose, for the sake of simplicity, that $\mathbf{u}$ is prescribed to zero on the whole boundary, and that the system is adiabatic, i.e. $\mathbf{q} \cdot \mathbf{n} = 0$ on $\partial\Omega$. Then, formally, taking the inner product of (1b) with $\mathbf{u}$ and integrating over $\Omega$, integrating (1c) over $\Omega$, and, finally, summing both relations yields the stability estimate:

$$ \frac{d}{dt} \int_{\Omega} \left[ \frac{1}{2} \rho |\mathbf{u}|^2 + \rho e \right] d\mathbf{x} \le 0. \qquad (3) $$

If we suppose that the equation of state may be set under the form $p = f(\rho, e)$ with $f(\cdot, 0) = 0$ and $f(0, \cdot) = 0$, Equation (1c) implies that $e$ remains positive (still at least formally), and so (3) yields a control on the unknown. Mimicking this computation at the discrete level necessitates to check some arguments, among them:

(i) to have available a discrete counterpart to the relation:

$$ \int_{\Omega} \left[ \partial_t(\rho \mathbf{u}) + \operatorname{div}(\rho \mathbf{u} \otimes \mathbf{u}) \right] \cdot \mathbf{u} \, d\mathbf{x} = \frac{d}{dt} \int_{\Omega} \frac{1}{2} \rho |\mathbf{u}|^2 \, d\mathbf{x}. $$

(ii) to identify the integral of the dissipation term at the right-hand side of the discrete counterpart of (1c) with the (discrete) $L^2$ inner product between the velocity and the diffusion term in the discrete momentum balance equation (1b).

(iii) to be able to prove that the right-hand side of (1c) is non-negative, in order to preserve the positivity of the internal energy.

The point (i) is extensively discussed in [6] (see also [7]), and is not treated here. Indeed, we focus here on a discretization technique which allows to obtain (ii) and (iii) with the usual Marker and Cell (MAC) discretization [4, 5], and which is implemented in the ISIS free software developed at IRSN [9] on the basis of the software component library PELICANS [11]. We complete the presentation by showing how (ii) may also be used, in some model problems, to prove the convergence in $L^1$ of the dissipation term.

# 2 Discretization of the dissipation term

## 2.1 The two-dimensional case

Let us begin with a two-dimensional case. The first step is to propose a discretization for the diffusion term in the momentum equation. We begin with the $x$-component of the velocity, for which we write a balance equation on $K_{i-\frac{1}{2},j}^x = (x_{i-1}, x_i) \times (y_{j-\frac{1}{2}}, y_{j+\frac{1}{2}})$ (see Figures 1 and 2 for the notations). Integrating the $x$ component of the momentum balance equation over $K_{i-\frac{1}{2},j}^x$, we get for the diffusion term:

$$ \bar{T}_{i-\frac{1}{2},j}^{\text{dif}} = - \left[ \int_{K_{i-\frac{1}{2},j}^x} \operatorname{div}[\tau(\mathbf{u})] d\mathbf{x} \right] \cdot \mathbf{e}^{(x)} = - \left[ \int_{\partial K_{i-\frac{1}{2},j}^x} \tau(\mathbf{u}) \, n \, d\gamma \right] \cdot \mathbf{e}^{(x)}, \quad (4) $$

where $\mathbf{e}^{(x)}$ stands for the first vector of the canonical basis of $\mathbb{R}^2$. We denote by $\sigma_{i,j}^x$ the right face of $K_{i-\frac{1}{2},j}^x$, i.e. $\sigma_{i,j}^x = \{x_i\} \times (y_{j-\frac{1}{2}}, y_{j+\frac{1}{2}})$. Splitting the boundary integral in (4), the part of $\bar{T}_{i-\frac{1}{2},j}^{\text{dif}}$ associated to $\sigma_{i,j}^x$, also referred to as the viscous flux through $\sigma_{i,j}^x$, reads:

$$ - \left[ \int_{\sigma_{i,j}^x} \tau(\mathbf{u}) \, n \, d\gamma \right] \cdot \mathbf{e}^{(x)} = -2 \int_{\sigma_{i,j}^x} \mu \, \partial_x \mathbf{u}^x \, d\gamma + \frac{2}{3} \int_{\sigma_{i,j}^x} \mu (\partial_x \mathbf{u}^x + \partial_y \mathbf{u}^y) \, d\gamma, $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Dual cell for the x-component of the velocity</figcaption></figure>

and the usual finite difference technique yields the following approximation for this term:

$$ -\frac{4}{3} \int_{\sigma_{i,j}^x} \mu \partial_x \mathbf{u}^x \mathrm{d}\gamma + \frac{2}{3} \int_{\sigma_{i,j}^x} \mu \partial_y \mathbf{u}^y \mathrm{d}\gamma \\ \approx -\frac{4}{3} \mu_{i,j} \frac{h_j^y}{h_i^x} (\mathbf{u}_{i+\frac{1}{2},j}^x - \mathbf{u}_{i-\frac{1}{2},j}^x) + \frac{2}{3} \mu_{i,j} \frac{h_j^y}{h_i^y} (\mathbf{u}_{i,j+\frac{1}{2}}^y - \mathbf{u}_{i,j-\frac{1}{2}}^y), \quad (5) $$

where $\mu_{i,j}$ is an approximation of the viscosity at the face $\sigma_{i,j}^x$. Similarly, let $\sigma_{i-\frac{1}{2},j+\frac{1}{2}}^x = (x_{i-1}, x_i) \times \{y_{j+\frac{1}{2}}\}$ be the top edge of the cell. Then:

$$ - \left[ \int_{\sigma_{i-\frac{1}{2},j+\frac{1}{2}}^x} \tau(\mathbf{u}) \, n \, d\gamma \right] \cdot e^{(x)} = - \int_{\sigma_{i-\frac{1}{2},j+\frac{1}{2}}^x} \mu (\partial_y \mathbf{u}^x + \partial_x \mathbf{u}^y) \, d\gamma \\ \approx -\mu_{i-\frac{1}{2},j+\frac{1}{2}} \left[ \frac{h_{i-\frac{1}{2}}^x}{h_{j+\frac{1}{2}}^y} (\mathbf{u}_{i-\frac{1}{2},j+1}^x - \mathbf{u}_{i-\frac{1}{2},j}^x) + \frac{h_{i-\frac{1}{2}}^x}{h_{i-\frac{1}{2}}^x} (\mathbf{u}_{i,j+\frac{1}{2}}^y - \mathbf{u}_{i-1,j+\frac{1}{2}}^y) \right], $$

where $\mu_{i-\frac{1}{2},j+\frac{1}{2}}$ stands for an approximation of the viscosity at the edge $\sigma_{i-\frac{1}{2},j+\frac{1}{2}}^x$.

Let us now multiply each discrete equation for $\mathbf{u}^x$ by the corresponding degree of freedom of a velocity field $\mathbf{v}$ (i.e. the balance over $K_{i-\frac{1}{2},j}^x$ by $\mathbf{v}_{i-\frac{1}{2},j}^x$) and sum over $i$ and $j$. The viscous flux at the face $\sigma_{i,j}^x$ appears twice in the sum, once multiplied by $\mathbf{v}_{i-\frac{1}{2},j}^x$ and the second one by $-\mathbf{v}_{i+\frac{1}{2},j}^x$, and the corresponding term reads:

$$ T_{i,j}^{\text{dis}}(\mathbf{u}, \mathbf{v}) = \mu_{i,j} \left[ -\frac{4}{3} \frac{h_j^y}{h_i^x} (\mathbf{u}_{i+\frac{1}{2},j}^x - \mathbf{u}_{i-\frac{1}{2},j}^x) + \frac{2}{3} \frac{h_j^y}{h_i^y} (\mathbf{u}_{i,j+\frac{1}{2}}^y - \mathbf{u}_{i,j-\frac{1}{2}}^y) \right] (\mathbf{v}_{i-\frac{1}{2},j}^x - \mathbf{v}_{i+\frac{1}{2},j}^x) \\ = \mu_{i,j} h_j^y h_i^x \left[ \frac{4}{3} \frac{\mathbf{u}_{i+\frac{1}{2},j}^x - \mathbf{u}_{i-\frac{1}{2},j}^x}{h_i^x} - \frac{2}{3} \frac{\mathbf{u}_{i,j+\frac{1}{2}}^y - \mathbf{u}_{i,j-\frac{1}{2}}^y}{h_j^y} \right] \frac{\mathbf{v}_{i+\frac{1}{2},j}^x - \mathbf{v}_{i-\frac{1}{2},j}^x}{h_i^x}. \quad (6) $$
Similarly, the term associated to $\sigma_{i-\frac{1}{2},j+\frac{1}{2}}^x$ appears multiplied by $v_{i-\frac{1}{2},j}^x$ and by $-v_{i-\frac{1}{2},j+1}^x$, and we get:

$$ T_{i-\frac{1}{2},j+\frac{1}{2}}^{\text{dis}}(\mathbf{u}, \mathbf{v}) = \mu_{i-\frac{1}{2},j+\frac{1}{2}} h_{i-\frac{1}{2},j+\frac{1}{2}}^x \left[ \frac{\mathbf{u}_{i-\frac{1}{2},j+1}^x - \mathbf{u}_{i-\frac{1}{2},j}^x}{h_{j+\frac{1}{2}}^y} + \frac{\mathbf{u}_{i,j+\frac{1}{2}}^y - \mathbf{u}_{i-1,j+\frac{1}{2}}^y}{h_{i-\frac{1}{2}}^x} \right] \frac{\mathbf{v}_{i-\frac{1}{2},j+1}^x - \mathbf{v}_{i-\frac{1}{2},j}^x}{h_{j+\frac{1}{2}}^y}. \quad (7) $$

Let us now define the discrete gradient of the velocity as follows:

- The derivatives involved in the divergence, $\partial_x^M \mathbf{u}^x$ and $\partial_y^M \mathbf{u}^y$, are defined over the primal cells by:

$$ \partial_x^M \mathbf{u}^x(\mathbf{x}) = \frac{\mathbf{u}_{i+\frac{1}{2},j}^x - \mathbf{u}_{i-\frac{1}{2},j}^x}{h_i^x}, \quad \partial_y^M \mathbf{u}^y(\mathbf{x}) = \frac{\mathbf{u}_{i,j+\frac{1}{2}}^y - \mathbf{u}_{i,j-\frac{1}{2}}^y}{h_j^y}, \quad \forall \mathbf{x} \in K_{i,j}. \quad (8) $$

- For the other derivatives, we introduce another mesh which is vertex-centred, and we denote by $K^{xy}$ the generic cell of this new mesh, with $K_{i+\frac{1}{2},j+\frac{1}{2}}^{xy} = (x_i, x_{i+1}) \times (y_j, y_{j+1})$. Then, $\forall \mathbf{x} \in K_{i+\frac{1}{2},j+\frac{1}{2}}^{xy}$:

$$ \partial_y^M \mathbf{u}^x(\mathbf{x}) = \frac{\mathbf{u}_{i+\frac{1}{2},j+1}^x - \mathbf{u}_{i+\frac{1}{2},j}^x}{h_{j+\frac{1}{2}}^y}, \quad \partial_x^M \mathbf{u}^y(\mathbf{x}) = \frac{\mathbf{u}_{i+1,j+\frac{1}{2}}^y - \mathbf{u}_{i,j+\frac{1}{2}}^y}{h_{i+\frac{1}{2}}^x}. \quad (9) $$

With this definition, we get:

$$ T_{i,j}^{\text{dis}}(\mathbf{u}, \mathbf{v}) = \mu_{i,j} \int_{K_{i,j}} \left[ \frac{4}{3} \partial_x^M \mathbf{u}^x - \frac{2}{3} \partial_y^M \mathbf{u}^y \right] \partial_x^M \mathbf{v}^x \, d\mathbf{x}, $$

and:

$$ T_{i-\frac{1}{2},j+\frac{1}{2}}^{\text{dis}}(\mathbf{u}, \mathbf{v}) = \mu_{i-\frac{1}{2},j+\frac{1}{2}} \int_{K_{i-\frac{1}{2},j+\frac{1}{2}}^{xy}} (\partial_y^M \mathbf{u}^x + \partial_x^M \mathbf{u}^y) \partial_y^M \mathbf{v}^x \, d\mathbf{x}. $$

Let us now perform the same operations for the y-component of the velocity. Doing so, we are lead to introduce an approximation of the viscosity at the edge $\sigma_{i-\frac{1}{2},j+\frac{1}{2}}^y = \{x_{i-\frac{1}{2}}\} \times (y_j, y_{j+1})$ (see Figure 2). Let us suppose that we take the same approximation as on $\sigma_{i-\frac{1}{2},j+\frac{1}{2}}^x$. Then, the same argument yields that multiplying each discrete equation for $\mathbf{u}^x$ and for $\mathbf{u}^y$ by the corresponding degree of freedom of a velocity field $\mathbf{v}$, we obtain a dissipation term which reads:

$$ T^{\text{dis}}(\mathbf{u}, \mathbf{v}) = \int_{\Omega} \tau^{\mathcal{M}}(\mathbf{u}) : \nabla^{\mathcal{M}} \mathbf{v} \, d\mathbf{x}, \quad (10) $$

where $\nabla^{\mathcal{M}}$ is the discrete gradient defined by (8)-(9) and $\tau^{\mathcal{M}}$ the discrete tensor:

$$ \tau^{\mathcal{M}}(\mathbf{u}) = \begin{bmatrix} 2\mu \partial_x^{\mathcal{M}} \mathbf{u}_x & \mu^{xy} (\partial_y^{\mathcal{M}} \mathbf{u}_x + \partial_x^{\mathcal{M}} \mathbf{u}_y) \\ \mu^{xy} (\partial_y^{\mathcal{M}} \mathbf{u}_x + \partial_x^{\mathcal{M}} \mathbf{u}_y) & 2\mu \partial_y^{\mathcal{M}} \mathbf{u}_y \end{bmatrix} - \frac{2}{3} \mu (\partial_x^{\mathcal{M}} \mathbf{u}_x + \partial_y^{\mathcal{M}} \mathbf{u}_y) I, \quad (11) $$

where $\mu$ is the viscosity defined on the primal mesh by $\mu(\mathbf{x}) = \mu_{i,j}$, $\forall \mathbf{x} \in K_{i,j}$ and $\mu^{xy}$ is the viscosity defined on the vertex-centred mesh, by $\mu(\mathbf{x}) = \mu_{i+\frac{1}{2},j+\frac{1}{2}}$, $\forall \mathbf{x} \in K_{i+\frac{1}{2},j+\frac{1}{2}}^{xy}$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Dual cell for the y-component of the velocity</figcaption></figure>

Now the form (10) suggests a natural to discretize the viscous dissipation term in the internal energy balance in order for the consistency property (ii) to hold. Indeed, if we simply set on each primal cell $K_{i,j}$:

$$ (\tau(\mathbf{u}) : \nabla \mathbf{u})_{i,j} = \frac{1}{|K_{i,j}|} \int_{K_{i,j}} \tau^{\mathcal{M}}(\mathbf{u}) : \nabla^{\mathcal{M}} \mathbf{u} \, d\mathbf{x}, \quad (12) $$

then, thanks to (10), the property (ii) which reads:

$$ T^{\text{dis}}(\mathbf{u}, \mathbf{u}) = \sum_{i,j} |K_{i,j}| (\tau(\mathbf{u}) : \nabla \mathbf{u})_{i,j}. $$

holds. Furthermore, we get from Definition (11) that $\tau^{\mathcal{M}}(\mathbf{u})(\mathbf{x})$ is a symmetrical tensor, for any $i, j$ and $\mathbf{x} \in K_{i,j}$, and therefore an elementary algebraic argument yields:

$$ \begin{aligned} (\tau(\mathbf{u}) : \nabla \mathbf{u})_{i,j} &= \frac{1}{|K_{i,j}|} \int_{K_{i,j}} \tau^{\mathcal{M}}(\mathbf{u}) : \nabla^{\mathcal{M}} \mathbf{u} \, d\mathbf{x} \\ &= \frac{1}{2|K_{i,j}|} \int_{K_{i,j}} \tau^{\mathcal{M}}(\mathbf{u}) : [\nabla^{\mathcal{M}}\mathbf{u} + (\nabla^{\mathcal{M}}\mathbf{u})^t] \, d\mathbf{x} \geq 0. \end{aligned} $$

*Remark 1 (Approximation of the viscosity)* Note that, for the symmetry of $\tau^{\mathcal{M}}(\mathbf{u})$ to hold, the choice of the same viscosity at the edges $\sigma_{i-\frac{1}{2},j+\frac{1}{2}}^x$ and $\sigma_{i-\frac{1}{2},j+\frac{1}{2}}^y$ is crucial even though other choices may appear natural. Assuming for instance the viscosity to be a function of an additional variable defined on the primal mesh, the following construction seems reasonable:

1. define a constant value for $\mu$ on each primal cell,

2. associate a value of $\mu$ to the primal edges, by taking the average between the value at the adjacent cells,

3. finally, split the integral of the shear stress over $\sigma_{i-\frac{1}{2},j+\frac{1}{2}}^x$ in two parts, one for the part included in the (top) boundary of $K_{i-1,j}$ and the second one in the boundary of $K_{i,j}$.
<figure><img src="image_4.png" /><figcaption>Figure 3: The xy-staggered cell $K_{i+\frac{1}{2}, j+\frac{1}{2}, k}^{xy}$, used in the definition of $\partial_y^M u^x$, $\partial_x^M u^y$, and $\tau^M(u)_{x,y} = \tau^M(u)_{y,x}$.</figcaption></figure>

Then the viscosities on $\sigma_{i-\frac{1}{2}, j+\frac{1}{2}}^x$ and $\sigma_{i-\frac{1}{2}, j+\frac{1}{2}}^y$ coincide only for uniform meshes, and, in the general case, the symmetry of $\tau^M(u)$ is lost.

## 2.2 Extension to the three-dimensional case

Extending the computations of the preceding section to three space dimensions yields the following construction.

- First, define three new meshes, which are "edge-centred": $K_{i+\frac{1}{2}, j+\frac{1}{2}, k}^{xy} = (x_i, x_{i+1}) \times (y_i, y_{i+1}) \times (z_{k-\frac{1}{2}}, z_{k+\frac{1}{2}})$ is staggered from the primal mesh $K_{i,j,k}$ in the $x$ and $y$ direction (see Figure 3), $K_{i+\frac{1}{2}, j+\frac{1}{2}}^{xz}$ in the $x$ and $z$ direction, and $K_{i,j+\frac{1}{2}, k+\frac{1}{2}}^{yz}$ in the $y$ and $z$ direction.

- The partial derivatives of the velocity components are then defined as piecewise constant functions, the value of which is obtained by natural finite differences:

- for $\partial_x^M u^x$, $\partial_y^M u^y$ and $\partial_z^M u^z$, on the primal mesh,

- for $\partial_y^M u^x$ and $\partial_x^M u^y$ on the cells ($K_{i+\frac{1}{2}, j+\frac{1}{2}, k}^{xy}$),

- for $\partial_z^M u^x$ and $\partial_x^M u^z$ on the cells ($K_{i+\frac{1}{2}, j, k+\frac{1}{2}}^{xz}$),

- for $\partial_y^M u^z$ and $\partial_z^M u^y$ on the cells ($K_{i,j+\frac{1}{2}, k+\frac{1}{2}}^{yz}$).

- Then, define four families of values for the viscosity field, $\mu$, $\mu^{xy}$, $\mu^{xz}$ and $\mu^{yz}$, associated to the primal and the three edge-centred meshes respectively.

- The shear stress tensor is obtained by the extension of (11) to $d=3$.

- And, finally, the dissipation term is given by (12).
### 3 A strong convergence result

We conclude this paper by showing how the consistency property (ii) may be used, in some particular cases, to obtain the strong convergence of the dissipation term, and then pass to the limit in a coupled equation having the dissipation term as right-hand side. To this purpose, let us just address the model problem:

$$ -\Delta \underline{u} = f \text{ in } \Omega = (0,1) \times (0,1), \quad \underline{u} = 0 \text{ on } \partial\Omega, \qquad (13) $$

with $\underline{u}$ and $f$ two scalar functions, $f \in L^2(\Omega)$. Let us suppose that this problem is discretized by the usual finite volume technique, with the uniform MAC mesh associated to the $x$-component of the velocity. We define a discrete function as a piecewise constant function, vanishing on the left and right sides of the domain (so on the left and right stripes of staggered (half-)meshes adjacent to these boundaries), and we define the discrete $H^1$-norm of a discrete function $v$ by:

$$ \|v\|_1^2 = \int_{\Omega} (\partial_x^M v)^2 + (\partial_y^M v)^2 \, d\mathbf{x}. $$

Let $(M^{(n)})_{n \in \mathbb{N}}$ be a sequence of such meshes, with a step $h^n$ tending to zero, and let $(u^{(n)})_{n \in \mathbb{N}}$ be the corresponding sequence of discrete solutions. Then, with the variational technique employed in the preceding section (i.e. multiplying each discrete equation by the corresponding unknown and summing), we get, with the usual discretization of the right-hand side:

$$ \|u^{(n)}\|_1^2 = \int_{\Omega} (\partial_x^M u^{(n)})^2 + (\partial_y^M u^{(n)})^2 \, d\mathbf{x} = \int_{\Omega} f u^{(n)} \, d\mathbf{x}. \qquad (14) $$

Since the discrete $H^1$-norm controls the $L^2$-norm (i.e. a discrete Poincaré inequality holds [2]), this yields a uniform bound for the sequence $(u^{(n)})_{n \in \mathbb{N}}$ in discrete $H^1$-norm. Hence the sequence $(u^{(n)})_{n \in \mathbb{N}}$ converges in $L^2(\Omega)$ to a function $\bar{u} \in H_0^1(\Omega)$, possibly up to the extraction of a subsequence [2], and the discrete derivatives $(\partial_x^M u^{(n)})_{n \in \mathbb{N}}$ and $(\partial_y^M u^{(n)})_{n \in \mathbb{N}}$ weakly converge in $L^2(\Omega)$ to $\partial_x \bar{u}$ and $\partial_y \bar{u}$ respectively. This allows to pass to the limit in the scheme, and we obtain that $\bar{u}$ satisfies the continuous equation (13) (so, since the solution to (13) is unique, the whole sequence converges). Thus, taking $\bar{u}$ as a test function in the variational form of (13):

$$ \int_{\Omega} (\partial_x \bar{u})^2 + (\partial_y \bar{u})^2 \, d\mathbf{x} = \int_{\Omega} f \bar{u} \, d\mathbf{x}. $$

But, passing to the limit in (14), we get:

$$ \lim_{n \to \infty} \int_{\Omega} (\partial_x^M u^{(n)})^2 + (\partial_y^M u^{(n)})^2 \, d\mathbf{x} = \lim_{n \to \infty} \int_{\Omega} f u^{(n)} \, d\mathbf{x} = \int_{\Omega} f \bar{u} \, d\mathbf{x}, $$

which, comparing to the preceding relation, yields:

$$ \lim_{n \to \infty} \int_{\Omega} (\partial_x^M u^{(n)})^2 + (\partial_y^M u^{(n)})^2 \, d\mathbf{x} = \int_{\Omega} (\partial_x \bar{u})^2 + (\partial_y \bar{u})^2 \, d\mathbf{x}. $$

Since the discrete gradient weakly converges and its norm converges to the norm of the limit, the discrete gradient strongly converges in $L^2(\Omega)^2$ to the gradient of the solution. Let us now imagine that Equation (13) is coupled to a balance equation for another variable, the right-hand side of which is $|\nabla \underline{u}|^2$; this situation occurs in several physical situations, as the modelling of
Joule effect [1], or RANS turbulence models [10, 3]. The discretization (12) of the dissipation term in the cell K, which reads here:

$$ (|\nabla u^{(n)}|^2)_K = \frac{1}{|K|} \int_K (\partial_x^M u^{(n)})^2 + (\partial_y^M u^{(n)})^2 d\mathbf{x}, $$

thus yields a convergent right-hand side, in the sense that, for any regular function $\varphi \in C_c^\infty(\Omega)$, we have:

$$ \lim_{n \to \infty} \sum_K \int_K (|\nabla u^{(n)}|^2)_K \varphi d\mathbf{x} = \int_{\Omega} |\nabla \underline{u}|^2 \varphi d\mathbf{x}. $$

(A declination of) this argument has been used to prove the convergence of numerical schemes in [1, 10, 3].

## References

[1] A. Bradji, R. Herbin: Discretization of the coupled heat and electrical diffusion problems by the finite element and the finite volume methods. IMA Journal of Numerical Analysis, **28**, 469–495 (2008).

[2] R. Eymard, T. Gallouët, R. Herbin, Finite Volume Methods. Handbook of Numerical Analysis, Volume VII, 713–1020, North Holland (2000)

[3] T. Gallouët, J.-C. Latché: Compactness of discrete approximate solutions to parabolic PDEs. Submitted (2011).

[4] F.H. Harlow, J.E. Welsh: Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Physics of Fluids, **8**, 2182–2189 (1965).

[5] F.H. Harlow, A.A. Amsden: A numerical fluid dynamics calculation method for all flow speeds. Journal of Computational Physics, **8**, 197–213 (1971).

[6] L. Gastaldo, R. Herbin, W. Kheriji, C. Lapuerta, J.-C. Latché: Staggered discretizations, pressure correction schemes and all speed barotropic flows. Finite Volumes for Complex Applications VI (FVCA VI), Prague, Czech Republic, June 2011.

[7] R. Herbin, J.-C. Latché: A kinetic energy control in the MAC discretization of compressible Navier-Stokes equations. International Journal of Finite Volumes **2** (2010).

[8] R. Herbin, W. Kheriji, J.-C. Latché: An unconditionally stable Finite Element-Finite Volume pressure correction scheme for compressible Navier-Stokes equations. In preparation (2011).

[9] ISIS: a CFD computer code for the simulation of reactive turbulent flows, https://gforge.irsn.fr/gf/project/isis.

[10] A. Larcher, J.-C. Latché: Convergence analysis of a finite element-finite volume scheme for a RANS turbulence model. Submitted (2011).

[11] PELICANS: Collaborative Development Environment. https://gforge.irsn.fr/gf/project/pelicans.