Equation and Conditions

Referred to literature R. Schwarze, a comparison and analyze of numerical schemes will be done. A convective flow field in a defined domain (see picture 1.1) will be applied and the scalar value of Φ is analyzed. The equation which is solved is a normal transport equation with time derivation, convective, diffusion and source terms. Due to the fact that we focus on convective term and the influence of schemes, mesh coarseness and orientation we can simplify the following equation:

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

with the following definitions:

  • The system will be solved till time derivation is zero (steady-state condition)
  • Diffusion coefficient D = 0
  • Incompressible fluid and no temperature depended density (constant)
  • No source terms

Then we can re-write the equation:

\[\nabla \bullet (\rho \textbf{U} \phi) = 0\]

Thus, this equation and the assumptions/boundary condition lead us to a physical system which is shown in picture 1.1. Without diffusion the scalar Φ must not mix with up while going through the domain. We will find a hard separation (dashed line) of the values. Now we will investigate into the mixing due to numerical schemes, numerical discretization and mesh dependency. Therefore, three different mesh refinements, mesh kinds (hexaedral, tetraedral) and in the kind of hexaedral meshes the orientation of the mesh edges are used. The solver used is a included in the standard OpenFOAM® libraries (scalarTransportFoam). Three meshes are defined and the discretization is applied to the outer edges:

  • 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 could occur that some schemes will blow up due to the time derivation and instability, although this scheme will produce good results if we choose a smaller time step. This behavior is due to the fact that some dimensionless numbers like the Courant or the Fourier number have to be in a certain range. If we go above a limit, we have in general two options a) decreasing the time step; b) numerical stabilization using different methods. This is the reason why we can observe different results by using the same time step.. That phenomenon was also discussed by Jósef Nagy from Linz. The test case is build as given below. The different meshes can be checked out below too.

Boundary Conditions

Structured mesh (transport under 0°)

  • Structured_20x20_0deg
  • Structured_40x40_0deg
  • Structured_80x80_0deg


Structured mesh (transport under 45°)

  • Structured_20x20_45deg
  • Structured_40x40_45deg
  • Structured_80x80_45deg


Unstructured mesh

  • Unstructured_20x20
  • Unstructured_40x40
  • Unstructured_80x80


Polygon mesh 

  • Polygon_20x20
  • Polygon_40x40
  • Polygon_80x80