1. General Finite Volume Formulation

1.1. Introduction

We consider a system of partial differential equations under the following conservative form :

(1.1)\[\frac{\partial U}{\partial t}+\vec{\nabla}.\vec{F}(U)=0\]

define on a domain \(\mathcal{V}\) with a contour \(\mathcal{C}_W\).

For Euler’s equations :

\(U = U(\vec{OM}, t) = \left(\begin{matrix} \rho \vec{u} \\ \rho E\\ \rho \end{matrix} \right)\)

and

\(\vec{F} = \left( \begin{matrix} \rho \vec{u} \vec{u} +p{\overline{\overline{I}}}_d \\ \rho \vec{u} \left(E+\frac{p}{\rho}\right) \\ \rho\vec{u} \\ \end{matrix} \right)\)

We seek a weak solution of this problem: \(\int_{\mathcal{V}}{\frac{\partial U}{\partial t}+\vec{\nabla}.\vec{F}(U)}dV=0\)

This equation is transformed thanks to Green-Ostrogradski theorem :

(1.2)\[\int_{\mathcal{V}}{\frac{\partial U}{\partial t} }dV + \int_{\mathcal{C}_W}{ \vec{F}(U) . \vec{n}}dC=0\]

1.2. Local equations

1.2.1. Local equation in elemental volume

The equation Eq.1.2 is plit over N elemental volumes (or cells) \(\mathcal{V}_I\) :

\(\sum_{I =1,N}{\int_{\mathcal{V}_I}{\frac{\partial U}{\partial t}dV}+\int_{\mathcal{C}_{I}}{\vec{F}\left(U\right).{\vec{n}}_IdC}}=0\)

with N the total number of cells.

To ensure this result one can write : \(\int_{\mathcal{V}_I}{\frac{\partial U}{\partial t}dV}+\int_{\mathcal{C}_{I}}{\vec{F}\left(U\right).{\vec{n}}_IdC}=0\) for each elementary volume \(\mathcal{V}_I\).

One must take care of the boundaries of each cell that may be a part of the global contour \(\mathcal{C}_W\) that is to say a part of the boundary. This is expressed, in the local equation :

\(\int_{\mathcal{V}_I}{\frac{\partial U}{\partial t}dV} + \int_{\mathcal{C}_{int,I}}{\vec{F}\left(U\right).{\vec{n}}_IdC} + \int_{\mathcal{C}_{W,I}}{\vec{F}\left(U\right).{\vec{n}}_IdC} =0\)

The term \(\mathcal{C}_{int,I}\) is for internal surfaces of the domain \(\mathcal{V}\) while \(\mathcal{C}_{W,I}\) is for boundary surface. Obviously, \(\mathcal{C}_{int,I} \cup \mathcal{C}_{W,I} = \mathcal{C}_{I}\) is the contour of the elementary cell \(\mathcal{V}_{I}\).

_images/elementalCells.jpg

Figure 1.1 : A whole domain \(\mathcal{V}\) with its contour \(\mathcal{C}_W\) and its elemntal cells. In orange are the internal contours of the cells. The cell I has some internal contour \(\mathcal{C}_int,I\) and a boundary contour (in red) \(\mathcal{C}_W,I\). Notice that the cell M has only internal contours.

Hint

One can notice \(\sum_{I}{\mathcal{C}_{W,I}} = \mathcal{C}_{W}\).

The finite volume method consists in numerically solve the equation define for each elemental volume :

(1.3)\[\int_{\mathcal{V}_I}{\frac{\partial U}{\partial t}dV} + \int_{\mathcal{C}_{I}}{\vec{F}\left(U\right).{\vec{n}}_IdC} =0\]

1.2.2. Properties

The summation of Eq.1.3 over all the elemantary cells must lead to Eq.1.1 .

\(\sum_{I}\int_{\mathcal{V}_I}{\frac{\partial U}{\partial t}dV}+\sum_{I}\int_{\mathcal{C}_{W,I}}{\vec{F}\left(U\right).{\vec{n}}_IdC}+\sum_{I}\int_{\mathcal{C}_{int,I}}{\vec{F}\left(U\right).{\vec{n}}_IdC}=0\)

The two first terms recovered Eq.1.1. One must verify that the third integral with internal fluxes is equal to zero.

(1.4)\[\sum_{I}\int_{\mathcal{C}_{\mathcal{int},I}}{\vec{F}\left(U\right).{\vec{n}}_IdC}=0\]

Let’s consider the cell number I. Its associated term in this integral reads :

\(\int_{\mathcal{C}_{\mathcal{int},I}}{\vec{F}\left(U\right).{\vec{n}}_IdC}\).

The integral over \(\mathcal{C}_{int,I}\) is split in a sum over partial countours \(\mathcal{C}_{int,I,k}\), labelled k :

\(\sum_{k}\int_{\mathcal{C}_{\mathcal{int},I,k}}{\vec{F}\left(U\right).{\vec{n}}_{I,k}dC}=0\)

(For example, the cell depicted in Figure 2.3, which has no boundary part, k is varying for 1 to 4).

Let’s consider the cell J, a neighbour of cell I. In Eq.1.4 let’s isolate the summation over both cells I and J :

\(\sum_{k}\int_{\mathcal{C}_{\mathcal{int},I,k}}{\vec{F}\left(U\right).{\vec{n}}_{I,k}dC} + \sum_{k}\int_{\mathcal{C}_{\mathcal{int},J,k}}{\vec{F}\left(U\right).{\vec{n}}_{J,k}dC}\)

The contact between both cells I and J is done through a part of a common contour, corresponding to the index K. This part of contour is \(\mathcal{C}_{int,I,K}\) or \(\mathcal{C}_{int,J,K}\). In the last equation, the contribution of this shared face is :

\(\int_{\mathcal{C}_{\mathcal{int},I,K}}{\vec{F}\left(U\right).{\vec{n}}_{I,K}dC} + \int_{\mathcal{C}_{\mathcal{int},J,K}}{\vec{F}\left(U\right).{\vec{n}}_{J,K}dC}\).

_images/CellAndContour.jpg

Figure 1.2 : Cells I and J in a domain \(\mathcal{V}\). In dashed green lines, the internal contours and in solid blue lines the boundary contours of each cell. Here, the common face of cell I and J is the face number 2 in cell I and the face number 4 in cell J. The numbering will be detailed in Mesh

To ensure the property of Eq.1.4, one must have:

\(\int_{\mathcal{C}_{\mathcal{int},I,K}}{\vec{F}\left(U\right).{\vec{n}}_{I,K}dC} + \int_{\mathcal{C}_{\mathcal{int},J,K}}{\vec{F}\left(U\right).{\vec{n}}_{J,K}dC} = 0\).

One can obvisouly remark that their normal vectors verify : \(\vec{n}_{I,K} = - \vec{n}_{J,K}\). Then, because this is the same contour K, one must verify :

\(\int_{\mathcal{C}_{\mathcal{int},I,K} \equiv \mathcal{C}_{\mathcal{int},J,K} }{ \left\{ \vec{F}_I\left(U\right).{\vec{n}}_{I,K}+\vec{F}_J\left(U\right).{\vec{n}}_{J,K} \right\} dC}=0\).

(1.5)\[\boxed{\vec{F}_I\left(U\right).{\vec{n}}_{I,K} + \vec{F}_J\left(U\right).{\vec{n}}_{J,K}=0}\]

Important

On the same face of two neighbouring cells the fluxes vanish.

1.3. Discrete Formulation

The discretization is the process leading to a numerical represention of the problem of Eq.1.3 in the computationnal domain in a form suitable for computing. Let’s define :

\(U_i^n=\frac{1}{\mathcal{V}_i}\int_{\mathcal{V}_i}{U(\vec{OM},\ t)dV}=\frac{1}{\mathcal{V}_i}\int_{\mathcal{V}_i}UdV\)

Because \(\mathcal{V}_i\) is a fixed volume, the temporal term of Eq.1.3 is :

\(\int_{\mathcal{V}_i}{\frac{\partial U}{\partial t}}dV=\frac{\partial}{\partial t} \int_{\mathcal{V}_i}{U}dV=\frac{\partial U_i^n \mathcal{V}_i }{\partial t}=\frac{\partial U_i^n }{\partial t} \mathcal{V}_i\)

One can approximate this derivative by a forwards 1rst order formulae : \(\frac{\partial U_i^n}{\partial t}V_i\approx\frac{U_i^{n+1}-U_i^n}{\Delta t}V_i\)

In this relation, \(\Delta t\) is the timestep that must be carefully chosen to avoid unstable behavior of the solution.

The integral over the contour of each cell : \(\int_{\mathcal{C}_I}{\vec{F}\left(U\right).{\vec{n}}_IdC}=\sum_{k}\int_{\mathcal{C}_{I,k}}{\vec{F}\left(U\right).{\vec{n}}_{I,k}dC}\) is now approximate by :

\(\int_{\mathcal{C}_I}{\vec{F}\left(U\right).{\vec{n}}_IdC}=\sum_{k}{{\vec{F}}_{I,k}.{\vec{n}}_{I,k}\int_{\mathcal{C}_{I,k}} d C=\sum_{k}{{\vec{F}}_{I,k}.{\vec{n}}_{I,k}A_{I,k}}}\)

where \(A_{I,k}\) is the area of the k-th part of the contour. It means that all the flow quantities are supposed constant along the k- th part of the contour.

\(\vec{F}_{I,k}\) is the vector of the fluxes evaluated on this surface. Numerically, the exact or an approximate solution of a Riemann problem provides an excellent solution. Then one must locally solve the following discrete equation:

(1.6)\[U_i^{n+1}=U_i^n-\frac{\Delta t}{V_i}\sum_{k=1, N_{face}}{{\vec{F}}_{i,k}.{\vec{n}}_{i,k}A_{I,k}}\]

Important

The fluxes \(\vec{F}_{i,k}\) depend on the conservative vector \(U_i\) but also on the neighbour of cell i at face k, that is to say i+1 or i-1 according the face considered (and in 1D). One then must read in Eq.1.6:

(1.7)\[{\vec{F}}_{i,k} = \vec{F}_{i,k}(U_i^n , U_{i\pm1}^n )\]