Nabla Logo programmer
  Advanced Search
  Send us a comment/query
  Back to main web site
  OpenFOAM guides
  - User guide
  - Programmer’s guide
  Changes from OpenFOAM 2.2 to 2.3
  Trademarks in the guides
  ©2000-2007 Nabla Ltd.

2.4 Equation discretisation

Equation discretisation converts the PDEs into a set of algebraic equations that are commonly expressed in matrix form as:

[A] [x] = [b]
\special {t4ht=

where [A]  \special {t4ht= is a square matrix, [x]  \special {t4ht= is the column vector of dependent variable and [b]  \special {t4ht= is the source vector. The description of [x]  \special {t4ht= and [b]  \special {t4ht= as ‘vectors’ comes from matrix terminology rather than being a precise description of what they truly are: a list of values defined at locations in the geometry, i.e. a geometricField<Type>, or more specifically a volField<Type> when using FV discretisation.

[A]  \special {t4ht= is a list of coefficients of a set of algebraic equations, and cannot be described as a geometricField<Type>. It is therefore given a class of its own: fvMatrix. fvMatrix<Type> is created through discretisation of a geometric<Type>Field and therefore inherits the <Type>. It supports many of the standard algebraic matrix operations of addition +, subtraction - and multiplication *.

Each term in a PDE is represented individually in OpenFOAM code using the classes of static functions finiteVolumeMethod and finiteVolumeCalculus, abbreviated by a typedef to fvm and fvc respectively. fvm and fvc contain static functions, representing differential operators, e.g.  \~/ 2   \special {t4ht=,  \~/  • \special {t4ht= and @/@t  \special {t4ht=, that discretise geometricField<Type>s. The purpose of defining these functions within two classes, fvm and fvc, rather than one, is to distinguish:

  • functions of fvm that calculate implicit derivatives of and return an fvMatrix<Type>
  • some functions of fvc that calculate explicit derivatives and other explicit calculations, returning a geometricField<Type>.

Figure 2.6 shows a geometricField<Type> defined on a mesh with 2 boundary patches and illustrates the explicit operations merely transform one field to another and drawn in 2D for simplicity.

                       geometricField<Type >

                          volField<Type >
                          surfaceField <Type >
                          pointField<Type >

                                            finiteVolumeCalculus  (fvc)
finiteVolumeMethod   (fvm)
      (Implicit)                              Other explicit operations

                                             geometricField<Type >
      fvMatrix<Type  >
                                               volField<Type >
                                               surfaceField<Type >
                                               pointField <Type >
\special {t4ht=

Figure 2.6: A geometricField<Type> and its operators

Table 2.2 lists the main functions that are available in fvm and fvc to discretise terms that may be found in a PDE. FV discretisation of each term is formulated by first integrating the term over a cell volume V  \special {t4ht=. Most spatial derivative terms are then converted to integrals over the cell surface S  \special {t4ht= bounding the volume using Gauss’s theorem

 integral               integral 

    \~/  *f dV  =    dS *f
 V               S
\special {t4ht=

where S  \special {t4ht= is the surface area vector, f  \special {t4ht= can represent any tensor field and the star notation *  \special {t4ht= is used to represent any tensor product, i.e. inner, outer and cross and the respective derivatives: divergence  \~/  •f  \special {t4ht=, gradient  \~/ f  \special {t4ht= and  \~/   f  \special {t4ht=. Volume and surface integrals are then linearised using appropriate schemes which are described for each term in the following Sections. Some terms are always discretised using one scheme, a selection of schemes is offered in OpenFOAM for the discretisation of other terms. The choice of scheme is either made by a direct specification within the code or it can be read from an input file at job run-time and stored within an fvSchemes class object.

Term description Implicit / Text fvm::/fvc:: functions
Explicit expression

Laplacian Imp/Exp   2
 \~/  f  \special {t4ht= laplacian(phi)
 \~/  •G \~/ f  \special {t4ht= laplacian(Gamma, phi)

Time derivative Imp/Exp @f
@t  \special {t4ht= ddt(phi)
 @t  \special {t4ht= ddt(rho,phi)

Second time derivative Imp/Exp @ (  @f )
--- r---
@t    @t \special {t4ht= d2dt2(rho, phi)

Convection Imp/Exp  \~/  •(y)  \special {t4ht= div(psi,scheme)*
 \~/  • (yf)  \special {t4ht= div(psi, phi, word)*
div(psi, phi)

Divergence Exp  \~/  •x  \special {t4ht= div(chi)

Gradient Exp  \~/ x  \special {t4ht= grad(chi)
 \~/ f  \special {t4ht= gGrad(phi)

Grad-grad squared Exp | \~/ \ ~/ f |2   \special {t4ht= sqrGradGrad(phi)

Curl Exp  \~/   f  \special {t4ht= curl(phi)

Source Imp rf  \special {t4ht= Sp(rho,phi)
Imp/Exp SuSp(rho,phi)

fvm::SuSp source is discretised implicit or explicit depending on the sign of rho.
An explicit source can be introduced simply as a vol<Type>Field, e.g. rho*phi.
Function arguments can be of the following classes:
phi: vol<Type>Field
Gamma: scalar volScalarField, surfaceScalarField, volTensorField, surfaceTensorField.
rho: scalar, volScalarField
psi: surfaceScalarField.
chi: surface<Type>Field, vol<Type>Field.

Table 2.2: Discretisation of PDE terms in OpenFOAM

2.4.1 The Laplacian term

The Laplacian term is integrated over a control volume and linearised as follows:

 integral                    integral                sum 
     \~/  •(G \~/ f) dV =     dS •(G \~/ f) =     GfSf •( \~/ f)f
  V                  S               f
\special {t4ht=

The face gradient discretisation is implicit when the length vector d  \special {t4ht= between the centre of the cell of interest P  \special {t4ht= and the centre of a neighbouring cell N  \special {t4ht= is orthogonal to the face plane, i.e. parallel to Sf  \special {t4ht=:

                 fN  - fP
Sf •( \~/ f)f = |Sf|---------
\special {t4ht=

In the case of non-orthogonal meshes, an additional explicit term is introduced [?] which is evaluated by interpolating cell centre gradients, themselves calculated by central differencing cell centre values.

2.4.2 The convection term

The convection term is integrated over a control volume and linearised as follows:

 integral                   integral                 sum                   sum 
    \~/  • (rUf) dV  =    dS •(rUf)  =     Sf •(rU)f ff =     F ff
 V                   S               f                  f
\special {t4ht=

The face field ff  \special {t4ht= can be evaluated using a variety of schemes:

Central differencing (CD)
is second-order accurate but unbounded
ff = fxfP  + (1-  fx)fN
      \special {t4ht=

where      ---- ----
fx  =_  f N /P N  \special {t4ht= where ----
f N  \special {t4ht= is the distance between f  \special {t4ht= and cell centre N  \special {t4ht= and ----
P N  \special {t4ht= is the distance between cell centres P  \special {t4ht= and N  \special {t4ht=.

Upwind differencing (UD)
determines ff  \special {t4ht= from the direction of flow and is bounded at the expense of accuracy
ff =   fP   for F > 0
       fN   for F < 0
      \special {t4ht=

Blended differencing (BD)
schemes combine UD and CD in an attempt to preserve boundedness with reasonable accuracy,
ff =  (1 - g) (ff)UD +  g (ff)CD
      \special {t4ht=

OpenFOAM has several implementations of the Gamma differencing scheme to select the blending coefficient g  \special {t4ht= [?] but it offers other well-known schemes such as van Leer, SUPERBEE, MINMOD etc..

2.4.3 First time derivative

The first time derivative @/@t  \special {t4ht= is integrated over a control volume as follows:

-@-   rf dV
@t  V
\special {t4ht=

The term is discretised by simple differencing in time using:

new values
fn   =_  f(t + Dt)  \special {t4ht= at the time step we are solving for;
old values
f    =_  f(t)  \special {t4ht= that were stored from the previous time step;
old-old values
foo  =_  f(t-  Dt)  \special {t4ht= stored from a time step previous to the last.

One of two discretisation schemes can be declared using the timeScheme keyword in the appropriate input file, described in detail in 4.4 of the User Guide.

Euler implicit
scheme, timeScheme EulerImplicit, that is first order accurate in time:
@   integral           (rPfP V )n-  (rPfP V )o
---   rf dV  = -----------------------
@t  V                    Dt
      \special {t4ht=

Backward differencing
scheme, timeScheme BackwardDifferencing, that is second order accurate in time by storing the old-old values and therefore with a larger overhead in data storage than EulerImplicit:
    integral                     n             o            oo
@--   rf dV =  3-(rP-fP-V)----4-(rPfP-V)--+-(rPfP-V-)--
@t  V                           2Dt
      \special {t4ht=

2.4.4 Second time derivative

The second time derivative is integrated over a control volume and linearised as follows:

    integral                      n             o           oo
-@-   r@f- dV =  (rPfP-V-)----2(rP-fPV-)-+--(rP-fP-V)--
@t  V   @t                        Dt2
\special {t4ht=

It is first order accurate in time.

2.4.5 Divergence

The divergence term described in this Section is strictly an explicit term that is distinguished from the convection term of Section  2.4.2, i.e. in that it is not the divergence of the product of a velocity and dependent variable. The term is integrated over a control volume and linearised as follows:

 integral               integral 
     \~/  •f dV =     dS •f =     Sf •ff
  V             S          f
\special {t4ht=

The fvc::div function can take as its argument either a surface<Type>Field, in which case ff  \special {t4ht= is specified directly, or a vol<Type>Field which is interpolated to the face by central differencing as described in Section  2.4.10:

2.4.6 Gradient

The gradient term is an explicit term that can be evaluated in a variety of ways. The scheme can be evaluated either by selecting the particular grad function relevant to the discretisation scheme, e.g. fvc::gGrad, fvc::lsGrad etc., or by using the fvc::grad function combined with the appropriate timeScheme keyword in an input file

Gauss integration
is invoked using the fvc::grad function with timeScheme Gauss or directly using the fvc::gGrad function. The discretisation is performed using the standard method of applying Gauss’s theorem to the volume integral:
 integral             integral          sum 
    \~/ f dV  =    dS f =     Sfff
 V            S         f
      \special {t4ht=

As with the fvc::div function, the Gaussian integration fvc::grad function can take either a surfaceField<Type> or a volField<Type> as an argument.

Least squares method
is based on the following idea:
  1. a value at point P  \special {t4ht= can be extrapolated to neighbouring point N  \special {t4ht= using the gradient at P  \special {t4ht=;
  2. the extrapolated value at N  \special {t4ht= can be compared to the actual value at N  \special {t4ht=, the difference being the error;
  3. if we now minimise the sum of the square of weighted errors at all neighbours of P  \special {t4ht= with the respect to the gradient, then the gradient should be a good approximation.

Least squares is invoked using the fvc::grad function with timeScheme leastSquares or directly using the fvc::lsGrad function. The discretisation is performed as by first calculating the tensor G  \special {t4ht= at every point P  \special {t4ht= by summing over neighbours N  \special {t4ht=:

      sum     2
G  =     w Ndd
      \special {t4ht=

where d  \special {t4ht= is the vector from P  \special {t4ht= to N  \special {t4ht= and the weighting function wN  = 1/|d| \special {t4ht=. The gradient is then evaluated as:

          sum     2  -1
( \~/ f)P =     w NG    •d (fN  - fP )
      \special {t4ht=

Surface normal gradient
The gradient normal to a surface nf • ( \~/ f)f  \special {t4ht= can be evaluated at cell faces using the scheme
         fN  - fP
( \~/ f)f = ---------
      \special {t4ht=

This gradient is called by the function fvc::snGrad and returns a surfaceField<Type>. The scheme is directly analogous to that evaluated for the Laplacian discretisation scheme in Section  2.4.1, and in the same manner, a correction can be introduced to improve the accuracy of this face gradient in the case of non-orthogonal meshes. This correction is called using the function fvc::snGradCorrection [Check**].

2.4.7 Grad-grad squared

The grad-grad squared term is evaluated by: taking the gradient of the field; taking the gradient of the resulting gradient field; and then calculating the magnitude squared of the result. The mathematical expression for grad-grad squared of f  \special {t4ht= is          2
| \~/  ( \~/ f) |   \special {t4ht=.

2.4.8 Curl

The curl is evaluated from the gradient term described in Section  2.4.6. First, the gradient is discretised and then the curl is evaluated using the relationship from Equation 2.7, repeated here for convenience

 \~/   f =  2 *(skew  \~/ f)  \special {t4ht=

2.4.9 Source terms

Source terms can be specified in 3 ways

Every explicit term is a volField<Type>. Hence, an explicit source term can be incorporated into an equation simply as a field of values. For example if we wished to solve Poisson’s equation   2
 \~/  f = f  \special {t4ht=, we would define phi and f as volScalarField and then do

        solve(fvm::laplacian(phi) == f)
An implicit source term is integrated over a control volume and linearised by
   rf dV  = r  V f
 V           P  P P
      \special {t4ht=

The implicit source term changes the coefficient of the diagonal of the matrix. Depending on the sign of the coefficient and matrix terms, this will either increase or decrease diagonal dominance of the matrix. Decreasing the diagonal dominance could cause instability during iterative solution of the matrix equation. Therefore OpenFOAM provides a mixed source discretisation procedure that is implicit when the coefficients that are greater than zero, and explicit for the coefficients less than zero. In mathematical terms the matrix coefficient for node P  \special {t4ht= is VP max(rP ,0)  \special {t4ht= and the source term is V f   min(r  ,0)
 P  P      P  \special {t4ht=.

2.4.10 Other explicit discretisation schemes

There are some other discretisation procedures that convert volField<Type>s into surface<Type>Fields and visa versa.

Surface integral
fvc::surfaceIntegrate performs a summation of surface<Type>Field face values bounding each cell and dividing by the cell volume, i.e.   sum 
(  f ff)/VP  \special {t4ht=. It returns a volField<Type>.
Surface sum
fvc::surfaceSum performs a summation of surface<Type>Field face values bounding each cell, i.e.  sum 
   f ff  \special {t4ht= returning a volField<Type>.
fvc::average produces an area weighted average of surface<Type>Field face values, i.e.  sum           sum 
(  f Sf ff)/  f Sf  \special {t4ht=, and returns a volField<Type>.
Face interpolate
The geometric<Type>Field function faceInterpolate() interpolates volField<Type> cell centre values to cell faces using central differencing, returning a surface<Type>Field.