Numerical Schemes of OpenFOAM®

For numerical analysis, one deals with partial differential equations. These equations are solved in time and space. Performing numerical simulations by using OpenFOAM®, one has to choose a discretization scheme for each partial differential term from a large list of numerical schemes. However, each discretization method has its numerical behavior. Thus, depending on the scheme in use, higher accuracy in terms of quantity representation (physically more correct) or numerical stability can be achieved. The following scheme analyzer allows one to investigate a wide range of different numerical schemes. The analysis is based on a simple convection equation. Tobias tried to make the analyzer flexible so that people can compare different meshes, schemes, and also OpenFOAM versions.

» Download the numerical case ...


— Tobias investigated into new scripting languages to develop the below presented analyzer —

The numerical schemes were analyzed using a simple convection equation. The equation, its derivation, and the conditions and the numerical meshes, are presented below. Previously, Tobias also added 2D slices for each numerical scheme, which is not added right now. This feature might come in the future if the community thinks it is worth extending the analyzer concerning this topic.

Analysis (1)
Comparison (2)






Equation and Conditions

Referred to literature R. Schwarze, a comparison and analysis of different numerical schemes were performed and can be investigated with the analyzer given above. A convective flow field in a defined domain (see picture 1.1) will be set-up, and the value of the scalar Φ is analyzed. The equation which is solved is an ordinary conserved transport equation including the time derivation, convective, diffusion, and source term. Since the investigation focuses onto numerical schemes, the equation is going to be simplified. The general conservation equation is given by:

\[\frac{\partial \rho \phi}{\partial t} + \nabla \bullet (\rho \textbf{U} \phi) - \nabla \bullet (D \nabla \phi) = S_\phi\]


while the quantities represent the general ones in computational fluid dynamics. Further information can be taken from the book »Mathematics, Numerics, Derivations and OpenFOAM®«. The following definitions and assumptions are applied to simplify the equation:

  • Diffusion coefficient D = 0, thus, the diffusion term vanishes
  • No source terms
  • Incompressible fluid and, therefore, no temperature-dependent density (constant)
  • The system will be solved until a steady-state condition is reached, the time derivation is included

 Hence, one can re-write the equation to:

\[\frac{\partial \phi}{\partial t} + \nabla \bullet (\rho \textbf{U} \phi) = 0\]


This equation is used in combination with the boundary conditions shown in picture 1.1. Without diffusion, the scalar Φ must keep a clear separation between the blue and yellow color; it should not mix as indicated by the dashed line. Using this setup, one can investigate into numerical schemes, while different phenomena can be shown, such as numerical diffusion, instability, and accuracy. Furthermore, the schemes are depended on the used numerical mesh. Hence, different meshes types were investigated as well as three different mesh-densities. Therefore, the analysis considered hexahedral, tetrahedral, and polyhedral meshes. For the hexahedral mesh, a 0° and 45° analysis was performed. The angle refers to the angle between the flow vector and the normal surface vector. The investigation was performed by using a standard OpenFOAM® solver named »scalarTransportFoam«. The mesh densities are equal to:

    20 x 20 cells
    40 x 40 cells
    80 x 80 cells

Important notice: All simulations are done with the same control settings. Hence, it can happen that some schemes will blow up due to the time derivation and numerical instability, although this scheme will produce good results if we choose a smaller time step. This behavior is related to some dimensionless numbers such as the Courant or Fourier number. These numbers always have to be in a certain »valid« range. If one goes above the validity limit, one commonly has two options to take this into consideration and allow a stable calculation: a) decreasing the time-step; b) numerical stabilization using different methods. This is the reason why one might observe different results for different schemes by using the same time-step. In other words, the analysis performed by Tobias Holzmann is consistent in the numerical set-up. However, for some schemes, this set-up does not fit to the numerical needs and thus the scheme does not provide accurate results. Jósef Nagy also discussed this phenomenon. The test case is identical to the one given on the right side. The different meshes are given below.

Numerical Schemes Case Setup

Structured Mesh (Hexahedral)

(transport under 0°)

Structured Mesh 20 x 20
Structured Mesh 40 x 40
Structured Mesh 80 x 80

Structured Mesh (Hexahedral)

(transport under 45°)

Structured Mesh 20 x 20
Structured Mesh 40 x 40
Structured Mesh 80 x 80

Unstructured Mesh (Tetrahedral)

Unstructured Mesh 20 x 20
Untructured Mesh 60 x 60
Unstructured Mesh 80 x 80

Unstructured Mesh (Polyhedral)

Polygon Mesh 20 x 20
Polygon Mesh 40 x 40