Convergence of a finite volume scheme for the incompressible fluids

TL;DR

This paper proves the convergence of a finite volume scheme for 2D incompressible Navier-Stokes equations using a triangular mesh and a projection method, building on previous stability results.

Contribution

It establishes the convergence of a specific finite volume scheme for incompressible fluids, extending prior stability analysis.

Findings

Scheme is stable and convergent for 2D incompressible Navier-Stokes

Uses triangular mesh with piecewise constant velocity and affine pressure

Builds on previous stability proof to show convergence

Abstract

We consider a finite volume scheme for the two-dimensional incompressible Navier-Stokes equations. We use a triangular mesh. The unknowns for the velocity and pressure are respectively piecewise constant and affine. We use a projection method to deal with the incompressibility constraint. In a former paper, the stability of the scheme has been proven. We infer from it its convergence.

Equations

$${\mathbf{u}}_t-\frac{1}{Re}\mathbf{\Delta u}+(\mathbf{u}\cdot \text{boldmath}\nabla \text{unboldmath})\mathbf{u}+\nabla p=\mathbf{f},$$

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

$$\mathbf{u}{|}_{\partial \Omega }=0,\mathbf{u}{|}_{t=0}={\mathbf{u}}_0.$$

$$\forall \sigma \in {\mathcal{E}}_h,d({\mathbf{x}}_{K_{\sigma }},\sigma )\ge C|\sigma |and|\sigma |\ge Ch.$$

$${\tau }_{\sigma }\ge C$$

$$|K|=\frac{1}{2}|\sigma |h_{K,\sigma }\ge \frac{1}{2}|\sigma |d({\mathbf{x}}_K,{\mathbf{x}}_{\sigma })\ge C{h}^2.$$

$$P_0=\{q\in L^2;\forall K\in {\mathcal{T}}_h,q{|}_{K}is\; a\; constant\},{\mathbf{P}}_0=(P_0{)}^2.$$

$$\Vert {\mathbf{v}}_h{\Vert }_h={\left(\sum _{\sigma \in {\mathcal{E}}_h^{int}}{\tau }_{\sigma }|{\mathbf{v}}_{L_{\sigma }}-{\mathbf{v}}_{K_{\sigma }}{|}^2+\sum _{\sigma \in {\mathcal{E}}_h^{ext}}{\tau }_{\sigma }|{\mathbf{v}}_{K_{\sigma }}{|}^2\right)}^{1/2}.$$

$$|{\mathbf{v}}_h|\le C\Vert {\mathbf{v}}_h{\Vert }_h.$$

$$\Vert {\mathbf{v}}_h{\Vert }_{-1,h}=\underset{{\text{boldmath}\psi \text{unboldmath}}_h\in {\mathbf{P}}_0}{\sup }\frac{({\mathbf{v}}_h,{\text{boldmath}\psi \text{unboldmath}}_h)}{\Vert {\text{boldmath}\psi \text{unboldmath}}_h{\Vert }_h}.$$

$$\forall K\in {\mathcal{T}}_h,({\Pi }_{{\mathbf{P}}_0}\mathbf{w}){|}_K=\frac{1}{|K|}{\int }_K\mathbf{w}(\mathbf{x})d\mathbf{x}.$$

$$\forall K\in {\mathcal{T}}_h,{\tilde{\Pi }}_{{\mathbf{P}}_0}\mathbf{v}{|}_K=\mathbf{v}({\mathbf{x}}_K).$$

$$|\mathbf{v}-{\tilde{\Pi }}_{{\mathbf{P}}_0}\mathbf{v}|\le Ch\Vert \mathbf{v}{\Vert }_2$$

$$P_1^d$$

$$P_1^{nc}$$

$${\nabla }_h{q}_h{|}_K=\frac{1}{|K|}{\int }_K\nabla q_{h}d\mathbf{x}.$$

$$\forall \sigma \in {\mathcal{E}}_h,{\int }_{\sigma }({\Pi }_{P_1^{nc}}q)d\sigma ={\int }_{\sigma }qd\sigma .$$

$$\left|\nabla q-{\nabla }_h({\Pi }_{P_1^{nc}}q)\right|\le Ch\Vert q{\Vert }_2,|\mathbf{v}-{\Pi }_{{\mathbf{P}}_1^{nc}}\mathbf{v}|\le Ch\Vert \mathbf{v}{\Vert }_1,$$

$$\mathbf{R}{\mathbf{T}}_0^{\mathbf{d}}$$

$$\mathbf{R}{\mathbf{T}}_0$$

$$\forall K\in {\mathcal{T}}_h,\forall \sigma \in {\mathcal{E}}_K,({\Pi }_{\mathbf{R}{\mathbf{T}}_0}\mathbf{v}\cdot {\mathbf{n}}_{K,\sigma }{)}_{\sigma }=\frac{1}{|\sigma |}{\int }_{\sigma }\mathbf{v}d\sigma .$$

$$|\mathbf{v}-{\Pi }_{\mathbf{R}{\mathbf{T}}_0}\mathbf{v}|\le Ch\Vert \mathbf{v}{\Vert }_1.$$

$$\forall \sigma \in {\mathcal{E}}_h^{int},({div}_h{\mathbf{v}}_h)({\mathbf{x}}_{\sigma })=\frac{3|\sigma |}{|K_{\sigma }|+|L_{\sigma }|}({\mathbf{v}}_{L_{\sigma }}-{\mathbf{v}}_{K_{\sigma }})\cdot {\mathbf{n}}_{K,\sigma };$$

$$\forall \sigma \in {\mathcal{E}}_h^{ext},({div}_h{\mathbf{v}}_h)({\mathbf{x}}_{\sigma })=-\frac{3|\sigma |}{|K_{\sigma }|+|L_{\sigma }|}{\mathbf{v}}_{K_{\sigma }}\cdot {\mathbf{n}}_{K,\sigma }.$$

$$\forall q_h\in P_1^{nc},{\Delta }_h{q}_h={div}_h({\nabla }_h{q}_h).$$

$${\tilde{\mathbf{\Delta }}}_h{\mathbf{v}}_h{|}_K=\frac{1}{|K|}\sum _{\sigma \in {\mathcal{E}}_K\cap {\mathcal{E}}_h^{int}}{\tau }_{\sigma }({\mathbf{v}}_{L_{\sigma }}-{\mathbf{v}}_{K_{\sigma }})-\frac{1}{|K|}\sum _{\sigma \in {\mathcal{E}}_K\cap {\mathcal{E}}_h^{ext}}{\tau }_{\sigma }{\mathbf{v}}_{K_{\sigma }}.$$

$${\widetilde{\hbox{\rm\bf b}}}_{h}({\bf u}_{h},{\bf v}_{h})\big{|}_{K}=\frac{1}{|K|}\sum_{\sigma\in{\mathcal{E}}_{K}\cap{\mathcal{E}}^{int}_{h}}\,|\sigma|\,\Big{(}({\bf u}\cdot{\mathbf{n}}_{K,\sigma})^{+}_{\sigma}\,{\bf v}_{K}+({\bf u}\cdot{\mathbf{n}}_{K,\sigma})^{-}_{\sigma}\,{\bf v}_{L_{\sigma}}\Big{)}.$$

$${\hbox{\rm b}}_{h}({\bf u}_{h},{\bf v}_{h},{\bf w}_{h})=\sum_{K\in{\mathcal{T}}_{h}}|K|\,{\bf w}_{K}\cdot{\widetilde{\hbox{\rm\bf b}}}_{h}({\bf u}_{h},{\bf v}_{h})\big{|}_{K}.$$

$${\mathbf{u}}_h^0\in {\mathbf{P}}_0\cap \mathbf{R}{\mathbf{T}}_0,{\mathbf{u}}_h^1\in {\mathbf{P}}_0\cap \mathbf{R}{\mathbf{T}}_0,p_h^1\in P_1^{nc}\cap L_0^2.$$

$$\frac{3{\tilde{\mathbf{u}}}_h^{n+1}-4{\mathbf{u}}_h^n+{\mathbf{u}}_h^{n-1}}{2k}$$

$$+{\tilde{\mathbf{b}}}_h(2{\mathbf{u}}_h^n-{\mathbf{u}}_h^{n-1},{\tilde{\mathbf{u}}}_h^{n+1})+{\nabla }_h{p}_h^n={\mathbf{f}}_h^{n+1},$$

$${\Delta }_h(p_h^{n+1}-p_h^n)=\frac{3}{2k}{div}_h{\tilde{\mathbf{u}}}_h^{n+1},$$

$${\mathbf{u}}_h^{n+1}={\tilde{\mathbf{u}}}_h^{n+1}-\frac{2k}{3}{\nabla }_h(p_h^{n+1}-p_h^n).$$

Full text

††journal: Numerische Mathematik

1 Introduction

We consider the flow of an incompressible fluid in a open bounded polyhedral set $\Omega\subset\mathbb{R}^{2}$ during the time interval $[0,T]$. The velocity field $\mathbf{u}:\Omega\times[0,T]\to\mathbb{R}^{2}$ and the pressure field $p:\Omega\times[0,T]\to\mathbb{R}$ satisfy the Navier-Stokes equations

[TABLE]

with the boundary and initial condition

[TABLE]

The terms $\mathbf{\Delta}{\bf u}$ and ({\bf u}\cdot{\hbox{\boldmath\nabla \unboldmath\!\!}}){\bf u} are associated with the physical phenomena of diffusion and convection, respectively. The Reynolds number Re measures the influence of convection in the flow. For equations (1)–(2), finite element and finite difference methods are well known and mathematical studies are available (see giraultr for example). For finite volume schemes, numerous computations have been conducted (kimchoi and boicaya for example). However, few mathematical results are available in this case. Let us cite Eymard and Herbin herb3 and Eymard, Latché and Herbin eymard . In order to deal with the incompressibility constraint, these works use a penalization method. Another way is to use the projection methods which have been introduced by Chorin chorin and Temam temam . This is the case in Faure faure where the mesh is made of squares. In Zimmermann zimm1 the mesh is made of triangles, which allows more complex geometries. In the present paper the mesh is also made of triangles, but we consider a different discretisation for the pressure. It leads to a linear system with a better-conditioned matrix. The layout of the article is the following. We first introduce (section 2.1) some notations and hypotheses on the mesh. We define (section 2.2) the spaces we use to approximate the velocity and pressure. We define also (section 2.3) the operators we use to approximate the differential operators in (1)–(2). By combining this with a projection method, we build the scheme in section 3. In order to provide a mathematical analysis, we state in section 4 that the differential operators in (1)–(2) and their discrete counterparts share similar properties. In particular, the discrete operators for the gradient and the divergence are adjoint. We then prove in section 5 the convergence of the scheme.

We conclude with some notations. We denote by $\chi_{I}$ the characteristic function of an interval $I\subset\mathbb{R}$. We denote by ${\mathcal{C}}^{\infty}_{0}={\mathcal{C}}^{\infty}_{0}(\Omega)$ the set of the functions with a compact support in $\Omega$. The spaces $(L^{2},|.|)$ and $(L^{\infty},\|.\|_{\infty})$ are the usual Lebesgue spaces and we set $L^{2}_{0}=\{q\in L^{2}\,;\int_{\Omega}q\,d{\bf x}=0\}$. Their vectorial counterparts are $({\bf L}^{2},|.|)$ and $({\bf L}^{\infty},\|.\|_{\infty})$ with ${\bf L}^{2}=(L^{2})^{2}$ and ${\bf L}^{\infty}=(L^{\infty})$. For $k\in\mathbb{N}^{*}$, $(H^{k},\|\cdot\|_{k})$ is the usual Sobolev space. Its vectorial counterpart is $({\bf H}^{k},\|.\|_{k})$ with ${\bf H}^{k}=(H^{k})^{2}$. For $k=1$, the functions of ${\bf H}^{1}$ with a null trace on the boundary form the space ${\bf H}^{1}_{0}$. Also, we set {\hbox{\boldmath\nabla \unboldmath\!\!}}{\bf u}=(\nabla u_{1},\nabla u_{2})^{T} if ${\bf u}=(u_{1},u_{2})\in{\bf H}^{1}$. If $\mathbf{X}\subset{\bf L}^{2}$ is a Banach space, we define ${\mathcal{C}}(0,T;\mathbf{X})$ (resp. $L^{2}(0,T;\mathbf{X})$) as the set of the applications $\mathbf{g}:[0,T]\to\mathbf{X}$ such that $t\to|\mathbf{g}(t)|$ is continous (resp. square integrable). The norm $\|.\|_{{\mathcal{C}}(0,T;\mathbf{X})}$ is defined by $\|\mathbf{g}\|_{{\mathcal{C}}(0,T;\mathbf{X})}=\sup_{s\in[0,T]}|\mathbf{g}(s)|$. Finally in all calculations, $C$ is a generic positive constant, depending only on $\Omega$, ${\bf u}_{0}$ and ${\bf f}$.

2 Discrete setting

First, we introduce the spaces and operators needed to build the mesh.

2.1 The mesh

Let ${\mathcal{T}}_{h}$ be a triangular mesh of $\Omega$: $\overline{\Omega}=\cup_{K\in{\mathcal{T}}_{h}}K$. For each triangle $K\in{\mathcal{T}}_{h}$, we denote by $|K|$ its area and ${\mathcal{E}}_{K}$ the set of his edges. If $\sigma\in{\mathcal{E}}_{K}$, ${\mathbf{n}}_{K,\sigma}$ is the unit vector normal to $\sigma$ pointing outwards of $K$.

The set of edges of the mesh is ${\mathcal{E}}_{h}=\cup_{K\in{\mathcal{T}}_{h}}{\mathcal{E}}_{K}$. The length of an edge $\sigma\in{\mathcal{E}}_{h}$ is $|\sigma|$. The set of edges inside $\Omega$ (resp. on the boundary) is ${\mathcal{E}}^{int}_{h}$ (resp. ${\mathcal{E}}^{ext}_{h}$): ${\mathcal{E}}_{h}={\mathcal{E}}^{int}_{h}\cup{\mathcal{E}}^{ext}_{h}$. If $\sigma\in{\mathcal{E}}^{int}_{h}$, $K_{\sigma}$ and $L_{\sigma}$ are the triangles sharing $\sigma$ as an edge. If $\sigma\in{\mathcal{E}}^{ext}_{h}$, only the triangle $K_{\sigma}$ inside $\Omega$ is defined.

We denote by ${\bf x}_{K}$ the circumcenter of a triangle $K$. We assume that the measure of all interior angles of the triangles of the mesh are below $\frac{\pi}{2}$, so that ${\bf x}_{K}\in K$. If $\sigma\in{\mathcal{E}}^{int}_{h}$ (resp. $\sigma\in{\mathcal{E}}^{ext}_{h}$ ) we set $d_{\sigma}=d({\bf x}_{K_{\sigma}},{\bf x}_{L_{\sigma}})$ (resp. $d_{\sigma}=d({\bf x}_{\sigma},{\bf x}_{K_{\sigma}})$). We define for all edge $\sigma\in{\mathcal{E}}_{h}$: $\tau_{\sigma}=\frac{|\sigma|}{d_{\sigma}}$. The maximum circumradius of the triangles of the mesh is $h$. We assume that there exists $C>0$ such that

[TABLE]

It implies that there exists a constant $C>0$ such that for all edge $\sigma\in{\mathcal{E}}_{h}$

[TABLE]

and for all triangles $K\in{\mathcal{T}}_{h}$ we have (with $\sigma\in{\mathcal{E}}_{K}$ and $h_{K,\sigma}$ the matching altitude)

[TABLE]

2.2 The discrete spaces

We first define

[TABLE]

For the sake of concision, we set for all $q_{h}\in P_{0}$ (resp. ${\bf v}_{h}\in{\bf P}_{0}$) and all triangle $K\in{\mathcal{T}}_{h}$: $q_{K}=q_{h}|_{K}$ (resp. ${\bf v}_{K}={\bf v}_{h}|_{K}$). Although ${\bf P}_{0}\not\subset{\bf H}^{1}$