2.4 Equation discretisation

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: 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:

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

2.4.2 The convection term

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, MINMOD etc..

2.4.3 First time derivative

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:

(2.22)

2.4.4 Second time derivative

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

(2.23)

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:

(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:

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:

(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:

1. a value at point can be extrapolated to neighbouring point using the gradient at ;
2. the extrapolated value at can be compared to the actual value at , the difference being the error;
3. 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 timeScheme leastSquares 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**].

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

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

2.4.9 Source terms

Source terms can be specified in 3 ways

Explicit

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 .

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