Finite Volume Method

The finite volume method is a numerical method for the discretization of conservation equations, i.e. of special, partial differential equations based on a conservation law.

If the finite volume method is used correctly, the conservation properties of these equations are preserved, which is why one speaks of a conservative discretization method. This is the main reason why the finite volume method has established itself in numerical fluid mechanics . But it is also used in numerical structural mechanics and electrical engineering.

In this method, the computation area is discretized by finite volumes, which can have any polygonal or polyhedral shape. This is why even complex geometries can be easily networked. The solution of the searched variable takes place in the center of these finite volumes.

Derivation

A conservation law is given by an equation of the kind

${\displaystyle {\frac {\partial }{\partial t}}u(x,t)+\nabla \cdot f(u(x,t))=g(u(x,t))}$

given in a field , where denotes the Nabla operator , which here stands for the divergence . The common case considered here is . The derivation for equations with other terms is analogous. First, the area is divided into a finite (finite) number of volume elements (see grid cells ). The law of conservation applies in each of these cells. If each of the cells fulfills the conditions of the Gaussian integral theorem , e.g. Lipschitz continuity of the edge, then the integration over a cell and conversion of the integral of the divergence into a surface integral results : ${\displaystyle \Omega \subset \mathbb {R} ^{d}\times [0,\infty )}$ ${\displaystyle \nabla }$ ${\displaystyle g(u)=0}$ ${\displaystyle \Omega _{i}}$

${\displaystyle \int _{\Omega _{i}}{\frac {\partial }{\partial t}}u\,\mathrm {d} \Omega +\int _{\partial \Omega _{i}}f(u)\cdot n\,\mathrm {d} S=0}$.

The change of a received quantity (e.g. the energy) in a cell can only happen by flowing out or flowing in (in this case of energy) over the edge of the cell. In each cell, the mean value of the maintenance quantity in this cell is calculated and, in the event that the cell does not change over time, an equation is obtained which describes the change in the mean values ​​in the cells over time: ${\displaystyle u_{i}={\frac {1}{|\Omega _{i}|}}\int _{\Omega _{i}}u\,\mathrm {d} \Omega }$

${\displaystyle {\frac {\partial }{\partial t}}u_{i}=-{\frac {1}{|\Omega _{i}|}}\int _{\partial \Omega _{i}}f(u)\cdot n\,\mathrm {d} S}$.

In numerical processes, cells with polygonal borders are usually selected, so that the integral over the edge can be represented as the sum of surface integrals over simple structures (straight edges in the two-dimensional case).

Solving the equation

A second-order Gaussian quadrature is generally used to calculate the surface integrals . After averaging the values ​​in the individual cells, the problem arises that the numerical solution is discontinuous along the grid edges. However, the situation at the edge can be seen as a Riemann problem . The use of an approximate Riemann solver then allows the calculation of the flows. The consistency of the Riemann solver is required here, which in this case means, on the one hand, continuity or even Lipschitz continuity , as well as the condition that it supplies the physical flow with identical data from both cells.

The system of ordinary differential equations to be created only provides this if an entropy condition is added. Because the purely mathematical consideration of the discontinuity at the cell edge allows not only the correct solution for the Riemann problem by means of a compression shock, but also an unphysical dilution shock. The entropy condition, however, excludes the dilution surge. The Riemann problem is then solved approximately (e.g. Osher ) or iteratively-exactly ( Godunov ) using numerical methods for ordinary differential equations (taking into account the entropy condition) .

1D example

The following convection equation is considered:

${\displaystyle {\frac {\partial u}{\partial t}}=-v{\frac {\partial u}{\partial x}}\,,\qquad t\geq 0\,,\qquad x\in \left[x_{\mathrm {min} },x_{\mathrm {max} }\right]}$

V describes the convection speed. This partial differential equation is to be discretized with the help of the finite volume method along the position coordinate. For this purpose, the location coordinate is first broken down into n discrete sections

${\displaystyle x_{\mathrm {min} }=x_{\frac {1}{2}}<x_{\frac {3}{2}}<\dots <x_{n-{\frac {1}{2}}}<x_{n+{\frac {1}{2}}}=x_{\mathrm {max} }}$

The control volumes , their midpoints and the sizes of the control volumes are defined by ${\displaystyle \Omega _{i}}$ ${\displaystyle x_{i}}$ ${\displaystyle \Delta x_{i}}$

${\displaystyle \Omega _{i}=\left[x_{i-{\frac {1}{2}}},x_{i+{\frac {1}{2}}}\right]}$

${\displaystyle x_{i}={\frac {x_{i-{\frac {1}{2}}}+x_{i+{\frac {1}{2}}}}{2}}}$

${\displaystyle \Delta x_{i}=x_{i+{\frac {1}{2}}}-x_{i-{\frac {1}{2}}}}$

Here, integer indices relate to the center of a control volume and non-integer indices relate to the edge of a control volume.

Now every term of the differential equation to be discretized is integrated over any internal control volume. For the accumulation term follows

${\displaystyle \int _{x_{i-{\frac {1}{2}}}}^{x_{i+{\frac {1}{2}}}}{\frac {\partial u}{\partial t}}\,\mathrm {d} x={\frac {\mathrm {d} }{\mathrm {d} t}}\int _{x_{i-{\frac {1}{2}}}}^{x_{i+{\frac {1}{2}}}}u\,\mathrm {d} x={\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}\Delta x_{i}}$

According to the Leibniz rule for parameter integrals, the integration and differentiation may be swapped, provided that the control volume cannot be changed over time. The discretized function value is then approximated with the aid of the mean value theorem of the integral calculus as an integral mean value of the actual function profile in the control volume. ${\displaystyle u_{i}}$

Then follows for the convective term with the help of the Gaussian integral theorem

${\displaystyle -v\int _{x_{i-{\frac {1}{2}}}}^{x_{i+{\frac {1}{2}}}}{\frac {\partial u}{\partial x}}\,\mathrm {d} x=-v\left(\left.u\right|_{x_{i+{\frac {1}{2}}}}-\left.u\right|_{x_{i-{\frac {1}{2}}}}\right)}$

Here the function values ​​on the edges of the control volumes must be approximated as a function of the known function values. A simple method for this is the so-called 1st-order UPWIND method, which states that the function value on the edge is approximated by the next known function value located upstream, i.e. H. ${\displaystyle \left.u\right|_{x_{i+{\frac {1}{2}}}}}$ ${\displaystyle u_{i}}$

${\displaystyle \left.u\right|_{x_{i+{\frac {1}{2}}}}={\begin{cases}u_{i}&:v>0\\u_{i+1}&:v<0\end{cases}}}$

If the two terms are put together again, one obtains the set of ordinary differential equations

${\displaystyle {\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}=-{\frac {v}{\Delta x_{i}}}\left(u_{i}-u_{i-1}\right)}$

which can be solved for ordinary differential equations using any method.

Higher Order Procedure

The method described so far is only first order due to the averaging of the values ​​in each cell. A higher order is achieved by using higher order polynomials in the cells. I.e. a distribution (constant, linear, parabolic, etc.) is assumed that receives the integral value.

The main difficulty here is that shock waves or shocks in the solution can lead to oscillations. To avoid this, total variation diminishing methods ( TVD methods ) are used, which do not increase the total variation and thus do not allow any new extremes (since polynomials interpolate discontinuities with overshoots). The most important classes of processes here are the flux limiter processes and the ENO processes (or WENO).

Convergence theory

Finite volume methods can be understood as special finite element methods for elliptic equations , in which one selects piece-wise constant or piece-wise linear shape functions that live on the cells and not on the grid points. This allows a convergence analysis with the help of the comprehensive theory there.

For parabolic or hyperbolic equations such as the Euler or Navier-Stokes equations , however, the mathematical convergence theory is less advanced.

Hyperbolic conservation equations

In the case of hyperbolic problems, shocks occur in particular , which make analysis considerably more difficult. Another difficulty here is that the solution of the equations is usually not unique. The set of Lax-Wendroff provides that a finite-volume method at convergence actually against a weak solution converges the equation. Entropy conditions or numerical viscosity are then used to show that this is actually the physically sensible one.

The convergence of a numerical approximation to an exact solution is defined by means of the global error : ${\displaystyle u_{\Delta t}}$ ${\displaystyle u}$ ${\displaystyle E_{\Delta t}(x,t):=u_{\Delta t}(x,t)-u(x,t)}$

${\displaystyle \|E_{\Delta t}(\cdot ,t)\|\longrightarrow 0\quad {\text{für}}\;\Delta t\rightarrow 0,\quad t\in \mathbb {R} ^{+}.}$

Another statement that applies to all finite volume methods is the necessary CFL condition that the numerical dependency area must contain the actual dependency area. Otherwise the process is unstable.

Convergence theory is particularly difficult for multidimensional equations. In the one-dimensional domain there are also results for methods of higher order which are based on the fact that the space of functions with limited variation yields compact sets in L ₁ .

consistency

A numerical procedure is a rule for the construction of the numerical approximation for the next time step: ${\displaystyle H_{\Delta t}}$

${\displaystyle u_{j}^{n+1}=H_{\Delta t}(U^{n},j),}$

where the time index, the location index and the approximate solution at the previous point in time of all location indices represent: ${\displaystyle n}$ ${\displaystyle j}$ ${\displaystyle U^{n}}$

${\displaystyle t=t_{n},\quad x=x_{j},\quad U^{n}=\left(u_{1}^{n},u_{2}^{n},\cdots \right).}$

The numerical procedure can also be defined continuously : ${\displaystyle H_{\Delta t}}$

${\displaystyle u(x,t+\Delta t)=H_{\Delta t}\left(u(\cdot ,t),x\right).}$

The local truncation error of the method is defined as: ${\displaystyle L_{\Delta t}}$ ${\displaystyle H_{\Delta t}}$

${\displaystyle L_{\Delta t}(x,t):={\frac {1}{\Delta t}}\left(u(x,t+\Delta t)-H_{\Delta t}\left(u(\cdot ,t),x\right)\right),}$

where the exact solution is currently assumed. The local truncation error is now used to define the consistency (numerics) of the procedure : ${\displaystyle u}$ ${\displaystyle t}$ ${\displaystyle H_{\Delta t}}$ ${\displaystyle L_{\Delta t}}$

${\displaystyle \|L_{\Delta t}(\cdot ,t)\|\longrightarrow 0\quad {\text{mit}}\;\Delta t\rightarrow 0,\quad t\in \mathbb {R} ^{+}.}$

The procedure has order of consistency , if there is one , so that: ${\displaystyle H_{\Delta t}}$ ${\displaystyle n\in \mathbb {N} }$ ${\displaystyle C>0}$

${\displaystyle \|L_{\Delta t}(\cdot ,t)\|\leq C(\Delta t)^{n}\quad {\text{mit}}\;\Delta t\rightarrow 0,\quad t\in \mathbb {R} ^{+}}$

gives.

software

The most popular commercial program packages for numerical flow simulation using FVM are Fluent and CFX from Ansys Inc and Star-CCM + from Siemens . Different codes are used in the aerospace industry, including codes developed by NASA , the flo codes by Antony Jameson , and at Airbus SE the codes ELSA and the TAU code of the German Aerospace Center . OpenFOAM is a software package released under the GNU General Public License .

history

The basic theoretical and practical ideas were developed for space travel from the 1950s , in particular by the Russian Godunow and the Hungarian Lax . The first finite volume methods come from Richard S. Varga (1962) and Preissmann (1961). The term finite volume method was then introduced independently of one another in 1971 by McDonald and in 1972 by MacCormack and Paullay for the solution of the time-dependent two-dimensional Euler equations.

The idea of ​​the approximate Riemann solver first appeared in the 1980s when Roe , Osher , van Leer and others also independently presented such methods.

literature

• Charles Hirsch: Numerical computation of internal and external flows. 2 volumes. Wiley & Sons, Chichester et al. 1988–1990 ( Wiley series in numerical methods in engineering ).
• Randall J. LeVeque: Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, Cambridge 2002, ISBN 0-521-00924-3 ( Cambridge Texts in Applied Mathematics ).