CFD Online URL
[Sponsors]
Home > Wiki > Finite volume

Finite volume

From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
(Basic Discretisation schemes)
m (Reverted edits by Reverse22 (talk) to last revision by Praveen)
 
(47 intermediate revisions not shown)
Line 1: Line 1:
-
== Discretisation Schemes for convective terms in General Transport Equation. Finite-Volume Formulation, structured grids  ==
+
The Finite Volume Method (FVM) is one of the most versatile discretization techniques used in CFD. Based on the control volume formulation of analytical fluid dynamics, the first step in the FVM is to divide the domain into a number of control volumes (aka cells, elements) where the variable of interest is located at the centroid of the control volume. The next step is to integrate the differential form of the governing equations (very similar to the control volume approach) over each control volume. Interpolation profiles are then assumed in order to describe the variation of the concerned variable between cell centroids. The resulting equation is called the discretized or discretization equation. In this manner, the discretization equation expresses the conservation principle for the variable inside the control volume.
 +
The most compelling feature of the FVM is that the resulting solution satisfies the conservation of quantities such as mass, momentum, energy, and species. This is '''exactly''' satisfied for any control volume as well as for the whole computational domain and '''for any number of control volumes'''. Even a coarse grid solution exhibits exact integral balances.<br>
-
== Introduction ==
+
FVM is the ideal method for computing discontinuous solutions arising in compressible flows. Any discontinuity must satisfy the [[Rankine-Hugoniot jump condition]] which is a consequence of conservation. Since finite volume methods are conservative they automatically satisfy the jump conditions and hence give physically correct weak solutions.
-
Here is describing the discretization schemes of the convective terms in the finite-volume equations. The accuracy, numerical stability and the boundness of the solution depends on the numerical scheme used for these terms. The central issue is the specification of an appropriate relationship between the convected variable, stored at the cell centre and its value at each of the cell faces.
+
FVM is also preferred while solving partial differential equations containing discontinuous coefficients.
-
== Basic Equations of CFD ==
+
==Integral form of conservation law==
-
All the conservation equations can be written in the same generic differential form:
+
Consider a single partial differential equation in conservation form
-
<table width="100%"><tr><td>
+
:<math>
-
:<math>  
+
\frac{\partial u}{\partial t} + \nabla \cdot f(u) = 0
-
  \frac {\partial( \rho \phi )} {\partial t} + \frac{\partial}{\partial x_{i}} \left( \rho U \phi - \Gamma_{\phi} \frac{\partial\phi}{\partial x_{i}}\right)=S_{\phi}   
+
</math>
</math>
-
</td><td width="5%">(1)</td></tr></table>
 
-
Equation (1) is integrated over a control volume and the following discretised equation for <math>\boldsymbol{\phi}</math> is produced:
+
where <math>u</math> is a conserved quantity and
-
<table width="100%"><tr><td>
+
:<math>
-
:<math>  
+
f(u) = [ f_1(u), \ldots, f_d(u)]
-
J_{h}- J_{l} + J_{n}- J_{s} + J_{e}- J_{w} + D_{h} - D_{l} + D_{n} - D_{s} + D_{e} - D_{s} = S_{p}
+
</math>
</math>
-
</td><td width="5%">(2)</td></tr></table>
 
-
  where <math>\boldsymbol{S_{p}}</math> is the source term for the control volume <math>\boldsymbol{P}</math>, and <math>\boldsymbol{J_{f}}</math> and <math>\boldsymbol{D_{f}}</math> represent, respectively, the convective and diffusive fluxes of <math>\boldsymbol{\phi}</math> across the control-volume face <math>\boldsymbol{f}</math>
+
is the flux vector. Integrating this over any volume <math>V</math> in <math>R^d</math> we get the '''integral form''' of the conservation law
-
<math>\boldsymbol{(f=h,l,n,s,e,w)}</math> 
+
-
The convective fluxes through the cell faces are calculated as:
+
:<math>
-
 
+
\frac{\partial}{\partial t} \int_V u dx + \oint_{\partial V} f_i n_i ds = 0
-
<table width="100%"><tr><td>
+
-
:<math>  
+
-
  J_{f}=C_{f}\phi_{f}
+
</math>
</math>
-
</td><td width="5%">(1)</td></tr></table>
 
-
     
 
-
  where <math>C_{f}</math> is the mass flow rate across the cell face <math>f</math>. The convected variable <math>\phi_{f}</math> associated with this mass flow rate is usually stored at the cell centres, and thus some form of interpolation assumption must be made in order to determine its value at each cell face. The interpolation procedure employed for this operation is the subject of the various schemes proposed in the literature and the accuracy, stability and boundedness of the solution depends on the procedure used.
 
-
 
 
-
In general, the value of <math>\boldsymbol{\phi_{f}}</math> can be explicity formulated in terms of its neighbouring nodal values by a functional relationship of the form:
 
-
<table width="100%"><tr><td>
+
where the [[Greens theorem | divergence theorem]] has been used to convert the divergence term to a surface integral over the surface <math>\partial V</math> and <math>(n_1, \ldots, n_d)</math> is the unit outward normal to <math>\partial V</math>. The above equation is the starting point for the finite volume method.
-
:<math>  
+
-
  \phi_{f}=P  \left( \phi_{nb} \right)
+
-
</math>
+
-
</td><td width="5%">(1)</td></tr></table>
+
-
 
+
-
where <math>\boldsymbol{\phi_{nb}}</math> denotes the neighbouring-node <math>\boldsymbol{\phi}</math>values.
+
-
Combining equations (\ref{eq3}) through (\ref{eq4a}), the discretised equation becomes:
+
-
<table width="100%"><tr><td>
+
==Finite volume formulation==
-
:<math>
+
-
  \left\{ D_{h} + C_{h} \left[ P \left( \phi_{nb} \right) \right]_{h} \right\} -
+
-
        \left\{ D_{l} + C_{l} \left[ P \left( \phi_{nb} \right) \right]_{l} \right\} + 
+
-
        \left\{ D_{n} + C_{n} \left[ P \left( \phi_{nb} \right) \right]_{n} \right\} -  
+
The computational domain <math>V_h</math> is divided into non-overlapping cells or finite volumes
-
        \left\{ D_{s} + C_{s} \left[ P \left( \phi_{nb} \right) \right]_{s} \right\} + 
+
-
        \left\{ D_{e} + C_{e} \left[ P \left( \phi_{nb} \right) \right]_{e} \right\} -
+
:<math>V_r, r=1,\ldots, N, \quad V_h = \cup_r V_r
-
        \left\{ D_{w} + C_{w} \left[ P \left( \phi_{nb} \right) \right]_{w} \right\} = S_{p} 
+
</math>
</math>
-
</td><td width="5%">(1)</td></tr></table>
 
-
== Convection Schemes ==
+
Usually these cells are polygons (triangles, quadrilaterals) in 2-D and polyhedra (tetrahedron, hexahedron, prims, etc) in 3-D.
-
All the convection schemes involve a stencil of cells in which the values of $\phi$ will be used to construct the face value $\phi_{f}$
+
Introduce the cell average value over the cell <math>V_r</math>
-
[[Image:picture_01.jpg]].
+
:<math>
 +
u_r = \frac{1}{|V_r|} \int_{V_r} u(x) dx
 +
</math>
-
+
which is the basic unknown quantity in the finite volume method. Let
-
\begin{figure}[h]
+
-
\centering
+
-
\includegraphics{D:/TeXWorkPlace/Wikki_Article/Picture_01.jpg}
+
-
\caption{Stensil of cells}
+
-
+
-
\label{fig:Picture_01}
+
-
\end{figure}
+
-
Where flow is from left to right, and $f$ is the face in question.
+
:<math>
-
 
+
N(r) = \{ \mbox{Set of cells which share a common face with } V_r \}
-
== Basic Discretisation schemes ==
+
-
 
+
-
=== Central Differencing Scheme (CDS)===
+
-
 
+
-
The most natural assumption for the cell-face value of the convected variable <math>\phi_{f}</math> would appear to be the CDS, which calculates the cell-face value from:
+
-
 
+
-
<table width="100%"><tr><td>
+
-
:<math>
+
-
  \phi_{f}=0.5 \left( \phi_{c} + \phi_{d} \right)
+
</math>
</math>
-
</td><td width="5%">(1)</td></tr></table>
 
-
 
 
-
This scheme is 2nd-order accurate, but is unbounded so that unphysical oscillations appear in regions of strong convection and also in the presence of discontinuities such as shocks. The CDS may be used directly in very low Reynolds-number flows where diffusive effects dominate over convection.
 
-
=== Upwind Differencing Scheme (UDS)===
+
Then the integral conservation law for cell <math>V_r</math> is
-
The UDS assumes that the convected variable at the cell fase <math>f</math> is the same as the upwind cell-centre value:
+
:<math>
-
 
+
|V_r| \frac{du_r}{dt} + \sum_{s \in N(r)} \int_{V_r \cap V_s} f_i n_i ds = 0
-
<table width="100%"><tr><td>
+
-
:<math>  
+
-
  \boldsymbol{\phi_{f}\phi_{c} }  
+
</math>
</math>
-
</td><td width="5%">(1)</td></tr></table>
 
-
   
 
-
The UDS is unconditionally bounded and highly stable, but as noted earlier it is only 1st-order accurate in terms of truncation error and may produce severe numerical diffusion. The scheme is therefore highly diffusive when the flow direction is skewed relative to the grid lines.
 
-
=== Hybrid Differencing Scheme (HDS) ===
+
It now remains to approximate the flux integral. This can be achieved using gaussian quadrature. Taking <math>p</math> gaussian points we get
-
The HDS of Spalding [1972] switches the discretisation of the convection terms between CDS and UDS according to the local cell Peclet number as follows:
+
:<math>
-
 
+
\int_{V_r \cap V_s} f_i n_i ds = \Delta s_{rs} \sum_{m=1}^p \omega_m F^m_{rs}
-
<table width="100%"><tr><td>
+
-
:<math>  
+
-
    \phi_{f}=0.5 \left( \phi_{c} + \phi_{d} \right) for Pe \triangleleft 2
+
-
 
+
-
 
+
-
\phi_{f}=  \phi_{c}        for Pe \triangleright 2
+
</math>
</math>
-
</td><td width="5%">(1)</td></tr></table>
 
-
     
 
-
The cell Peclet number is defined as:
 
-
 
-
<table width="100%"><tr><td>
 
-
:<math>
 
-
  Pe= \rho \left| U_{f} \right| A_{f}/D_{f}
 
-
</math>
 
-
</td><td width="5%">(1)</td></tr></table>
 
-
 
-
in which $A_{f}$ and $D_{f}$ are respectively, the cell-face area and physical diffusion coefficient. When $Pe\triangleright 2$,CDS calculations tends to become unstable so that theHDS reverts to the UDS. Physical diffusion is ignored when $Pe\triangleright 2$ .
 
-
 
-
The HDS scheme is marginally more accurate than the UDS, because the2nd-order CDS will be used in regions of low Peclet number.
 
-
 
-
== High Resolution Schemes (HRS) ==
 
-
 
-
=== Classification of High Resolution Schemes ===
 
-
 
-
HRS can be classified as \textit{linear} or \textit{non-linear}, where \textit{linear} means their coefficients are not direct functions of the convected variable when applied to a linear convection equation. It is important to recognise that linear convection schemes of 2nd-order accuracy or higher may suffer from unboudedness, and are not unconditionally stable.
 
-
 
-
\textit{Non-linear} schemes analyse the solution within the stencil and adapt the discretisation to avoid any unwanted behavior, such as unboundedness (see Waterson [1994]). These two types of schemes may be presented in a unified way by use of the \textsl{Flux-Limiter} formulation (Waterson and Deconinck [1995]), which calculates the cell-face value of the convected variable from:
 
-
 
-
  \begin{equation} 
 
-
    \phi_{f}= \phi_{c} + 0.5 B \left( r \right) \left( \phi_{c}-\phi_{u} \right)   
 
-
  \label{eq9}
 
-
  \end{equation}
 
-
     
 
-
  where $B \left( r \right)$ is termed a limiter function and the gradient ration $r$ is defined as:
 
-
 
-
    \begin{equation} 
 
-
    r= \left( \phi_{d} - \phi_{c}  \right) / \left( \phi_{c} - \phi_{u}  \right)
 
-
  \label{eq10}
 
-
  \end{equation}
 
-
 
 
-
The generalisation of this approach to handle non-uniform meshes has been given by Waterson [1994]
 
-
 
-
From equation (\ref{eq9}) it can be seen that $B=1$ gives the UDS and $B=r$ gives the CDS.
 
-
 
-
=== Numerical Implementation of HRS ===
 
-
 
-
The HRS schemes can be introdused into equation (\ref{eq4b}) by using the deffered correction procedure of Rubin and Khosla [1982]. This procedure express the cell-face value $\phi_{f}$ by:
 
-
 
-
    \begin{equation} 
 
-
    \phi_{f}=\phi_{f}\left(U \right) + \phi^{'}_{f}
 
-
    \label{eq10a}
 
-
  \end{equation}
 
-
 
 
-
where $\phi^{'}_{f}$ is a higher-order correction which represents the difference between the UDS face value $\phi_{f}\left(U \right)$ and the higher-order scheme value $\phi_{f}\left(H \right)$, i.e.
 
-
 
-
    \begin{equation} 
 
-
    \phi^{'}_{f}= \phi_{f}\left(H \right) +  \phi_{f}\left(U \right)
 
-
    \label{eq10b}
 
-
  \end{equation}
 
-
 
-
If equation (\ref{eq10a}) is substituted into equation (\ref{eq4b}), the resulting discretised equation is:
 
-
    \begin{equation} 
 
-
\left\{ D_{h} + C_{h} \phi_{h} \left( U \right) \right\} -
 
-
\left\{ D_{l} + C_{l} \phi_{l} \left( U \right) \right\} +
 
-
 
-
\left\{ D_{n} + C_{n} \phi_{n} \left( U \right) \right\} -
 
-
\left\{ D_{s} + C_{s} \phi_{s} \left( U \right) \right\} +
 
-
 
-
\left\{ D_{e} + C_{e} \phi_{e} \left( U \right) \right\} -
 
-
\left\{ D_{w} + C_{w} \phi_{w} \left( U \right) \right\} = S_{p} + B_{p}
 
-
    \label{eq11}
 
-
  \end{equation}
 
-
 
-
\\ .
 
-
\\ where $B_{p}$ is the deferred-correction source terms, given by:
 
-
 
-
    \begin{equation} 
 
-
  B_{p} = C_{l}\phi^{'}_{l} - C_{h}\phi^{'}_{h} +
 
-
  C_{s}\phi^{'}_{s} - C_{n}\phi^{'}_{n} +
 
-
  C_{w}\phi^{'}_{w} - C_{e}\phi^{'}_{e}
 
-
    \label{eq12}
 
-
  \end{equation}
 
-
 
 
-
  This treatment leads to a diagonally dominant coefficient matrix since it is formed using the UDS.
 
-
 
 
-
  The final form of the discretised equation:
 
-
 
 
-
  \begin{equation} 
 
-
  a_{P}\phi_{P}= a_{N}\phi_{N} + a_{S}\phi_{S} +
 
-
  a_{E}\phi_{E} + a_{W}\phi_{W} +
 
-
  a_{H}\phi_{H} + a_{L}\phi_{L} + a_{T}\phi_{T} + S_{p} +  B_{p}
 
-
  \label{eq2}
 
-
  \end{equation} 
 
-
 
 
-
Subscrit $P$ represents the current computational cell; $N$, $S$, $E$, $W$, $H$, $L$ represent the six neighbouring cells and $T$ represents the previous timestep (transistent cases only)
 
-
 
-
The coefficients contain the appropriate contributions from the transient, convective and diffusive terms in  (\ref{eq1})
 
-
 
-
== Normalised Variables ==
 
-
 
-
== Total Variation Diminishing (TVD) ==
 
-
 
-
== Convection Boundedness Criterion (CBC) ==
 
-
 
-
== Schemes ==
 
-
 
-
QUICK
 
-
 
-
LUS
 
-
 
-
Fromm
 
-
 
-
CUS
 
-
 
-
------------------------------
 
-
 
-
 
-
van Leer
 
-
 
-
OSPRE
 
-
 
-
Superbee
 
-
 
-
MINMOD
 
-
 
-
ISNAS
 
-
 
-
MUSCL
 
-
 
-
UMIST
 
-
 
-
SOUCUP
 
-
 
-
HLPA
 
-
 
-
SMART
 
-
SMARTER
+
where <math>\omega_m</math> are the gaussian weights and <math>F</math> is an approximation to <math>f_i n_i</math>.
-
LPPA
+
== External links ==
 +
* [http://www.imtek.uni-freiburg.de/simulation/mathematica/imsReferencePointers/FVM_introDocu.html The Finite Volume Method (FVM) - An introduction] by Oliver Rübenkönig of Albert Ludwigs University of Freiburg, available under the GNU Free Document License|GFDL.
-
SHARP
 
-
CHARM
+
<i> Return to [[Numerical methods | Numerical Methods]] </i>
-
VONOS
+
{{stub}}

Latest revision as of 09:15, 3 January 2012

The Finite Volume Method (FVM) is one of the most versatile discretization techniques used in CFD. Based on the control volume formulation of analytical fluid dynamics, the first step in the FVM is to divide the domain into a number of control volumes (aka cells, elements) where the variable of interest is located at the centroid of the control volume. The next step is to integrate the differential form of the governing equations (very similar to the control volume approach) over each control volume. Interpolation profiles are then assumed in order to describe the variation of the concerned variable between cell centroids. The resulting equation is called the discretized or discretization equation. In this manner, the discretization equation expresses the conservation principle for the variable inside the control volume.

The most compelling feature of the FVM is that the resulting solution satisfies the conservation of quantities such as mass, momentum, energy, and species. This is exactly satisfied for any control volume as well as for the whole computational domain and for any number of control volumes. Even a coarse grid solution exhibits exact integral balances.

FVM is the ideal method for computing discontinuous solutions arising in compressible flows. Any discontinuity must satisfy the Rankine-Hugoniot jump condition which is a consequence of conservation. Since finite volume methods are conservative they automatically satisfy the jump conditions and hence give physically correct weak solutions.

FVM is also preferred while solving partial differential equations containing discontinuous coefficients.

Integral form of conservation law

Consider a single partial differential equation in conservation form


\frac{\partial u}{\partial t} + \nabla \cdot f(u) = 0

where u is a conserved quantity and


f(u) = [ f_1(u), \ldots, f_d(u)]

is the flux vector. Integrating this over any volume V in R^d we get the integral form of the conservation law


\frac{\partial}{\partial t} \int_V u dx + \oint_{\partial V} f_i n_i ds = 0

where the divergence theorem has been used to convert the divergence term to a surface integral over the surface \partial V and (n_1, \ldots, n_d) is the unit outward normal to \partial V. The above equation is the starting point for the finite volume method.

Finite volume formulation

The computational domain V_h is divided into non-overlapping cells or finite volumes

V_r, r=1,\ldots, N, \quad V_h = \cup_r V_r

Usually these cells are polygons (triangles, quadrilaterals) in 2-D and polyhedra (tetrahedron, hexahedron, prims, etc) in 3-D.

Introduce the cell average value over the cell V_r


u_r = \frac{1}{|V_r|} \int_{V_r} u(x) dx

which is the basic unknown quantity in the finite volume method. Let


N(r) = \{ \mbox{Set of cells which share a common face with } V_r \}

Then the integral conservation law for cell V_r is


|V_r| \frac{du_r}{dt} + \sum_{s \in N(r)} \int_{V_r \cap V_s} f_i n_i ds = 0

It now remains to approximate the flux integral. This can be achieved using gaussian quadrature. Taking p gaussian points we get


\int_{V_r \cap V_s} f_i n_i ds = \Delta s_{rs} \sum_{m=1}^p \omega_m F^m_{rs}

where \omega_m are the gaussian weights and F is an approximation to f_i n_i.

External links


Return to Numerical Methods


My wiki