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

(2.12)

where is a square matrix, is the column vector of dependent variable
and is the source vector. The description of and 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.

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. , and
, 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.

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 . Most spatial
derivative terms are then converted to integrals over the cell surface bounding
the volume using Gauss’s theorem

(2.13)

where is the surface area vector, can represent any tensor field and the star
notation is used to represent any tensor product, i.e. inner, outer and cross
and the respective derivatives: divergence , gradient and .
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

laplacian(phi)

laplacian(Gamma, phi)

Time derivative

Imp/Exp

ddt(phi)

ddt(rho,phi)

Second time derivative

Imp/Exp

d2dt2(rho, phi)

Convection

Imp/Exp

div(psi,scheme)*

div(psi, phi, word)*

div(psi, phi)

Divergence

Exp

div(chi)

Gradient

Exp

grad(chi)

gGrad(phi)

lsGrad(phi)

snGrad(phi)

snGradCorrection(phi)

Grad-grad squared

Exp

sqrGradGrad(phi)

Curl

Exp

curl(phi)

Source

Imp

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:

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

(2.14)

The face gradient discretisation is implicit when the length vector between the
centre of the cell of interest and the centre of a neighbouring cell is
orthogonal to the face plane, i.e. parallel to :

(2.15)

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.

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

(2.16)

The face field can be evaluated using a variety of schemes:

Central differencing (CD)

is second-order accurate but unbounded

(2.17)

where where is the distance between and cell
centre and is the distance between cell centres and
.

Upwind differencing (UD)

determines from the direction of flow and is
bounded at the expense of accuracy

(2.18)

Blended differencing (BD)

schemes combine UD and CD in an attempt to
preserve boundedness with reasonable accuracy,

(2.19)

OpenFOAM has several implementations of the Gamma differencing scheme to
select the blending coefficient [?] but it offers other well-known schemes
such as van Leer, SUPERBEE, MINMODetc..

The first time derivative is integrated over a control volume as
follows:

(2.20)

The term is discretised by simple differencing in time using:

new values

at the time step we are solving for;

old values

that were stored from the previous time step;

old-old values

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:

(2.21)

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:

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:

(2.24)

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

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 timeSchemeGauss or directly using the fvc::gGrad function. The discretisation is
performed using the standard method of applying Gauss’s theorem to the
volume integral:

(2.25)

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:

a value at point can be extrapolated to neighbouring point
using the gradient at ;

the extrapolated value at can be compared to the actual value
at , the difference being the error;

if we now minimise the sum of the square of weighted errors at
all neighbours of with the respect to the gradient, then the
gradient should be a good approximation.

Least squares is invoked using the fvc::grad function with timeSchemeleastSquares or directly using the fvc::lsGrad function. The
discretisation is performed as by first calculating the tensor at every
point by summing over neighbours :

(2.26)

where is the vector from to and the weighting function
. The gradient is then evaluated as:

(2.27)

Surface normal gradient

The gradient normal to a surface can
be evaluated at cell faces using the scheme

(2.28)

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**].

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 is .

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

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 , we would define phi and
f as volScalarField and then do

solve(fvm::laplacian(phi) == f)

Implicit

An implicit source term is integrated over a control volume and
linearised by

(2.29)

Implicit/Explicit

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 is and the source term is
.

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. . It returns a volField<Type>.

Surface sum

fvc::surfaceSum performs
a summation of surface<Type>Field face values bounding each cell, i.e.
returning a volField<Type>.

Average

fvc::average produces an area weighted average of
surface<Type>Field face values, i.e. , and returns a
volField<Type>.

Reconstruct

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.