Close

17th February 2014

OpenFOAM 2.3.0: Numerical Methods

Polynomial Surface Normal Gradient

The surface normal gradient calculation, snGrad, is integral to pressure-velocity algorithms on unstructured meshes and the choice of method has an important impact on both accuracy and stability. 2nd-order accuracy on arbitrary unstructured meshes cannot be achieved using an snGrad calculated on a face only from its neighbouring cells (even with non-orthogonal correction). Instead, a larger stencil is required to which a higher-order polynomial may be fitted. Using the generalised stencil mechanism provided in a previous OpenFOAM this version now includes the following snGrad schemes:

  • quadraticFit, 2nd-order, using a CFC stencil;
  • linearFit, 1st-order using FEC stencil.

The quadraticFit scheme may be used on good quality meshes for runs requiring formal 2nd-order accuracy on all terms but at the cost of additional run-time and lower stability and robustness. The linearFit scheme is provided for completeness, providing similar accuracy to the standard corrected scheme, but using a more compact stencil.

Extended ddtCorr

The ddtPhiCorr calculation reduces decoupling between pressure, velocity and velocity flux in the pressure-velocity coupling algorithm for unstructured meshes. However, it is difficult to include in moving-mesh solvers due to the way in which the relative flux couples to the motion of the faces.

In version 2.3.0, ddtPhiCorr is replaced with the more general ddtCorr method which for moving-mesh solvers is based on the face velocity rather than the face flux, thus handling arbitrary motion of the face including translation and rotation, which was not handled accurately in the previous ddtPhiCorr formulation.

Source code

  • finiteVolume library
    $FOAM_SRC/finiteVolume

Ordinary Differential Equation Solvers

The semi-implicit Burlish-Stoer (SIBS) solver is recommended for stiff systems of ordinary differential equations (ODEs). It is usually reliable and efficient, but for some systems it converges poorly or fails to converge. This recently prompted a detailed investigation into ODE solver research and development to find the most promising methods for application to OpenFOAM cases, making particular use of Solving Ordinary Differential Equations Hairer et al. volume 1 and volume 2 A set of algorithms presented in those texts (and some other texts stated in the code documentation) have been implemented in OpenFOAM and tested on a range of cases. The set of ODE solvers, now provided with OpenFOAM-2.3.0, includes:

  • Euler, Euler ODE solver of order (0)1;
  • EulerSI, Semi-implicit Euler ODE solver of order (0)1;
  • Trapezoid, Trapezoidal ODE solver of order (1)2;
  • RKCK45, Cash-Karp Runge-Kutta ODE solver of order (4)5;
  • RKDP45, DormandPrince Runge-Kutta ODE solver of order (4)5;
  • RKF45, Runge-Kutta-Fehlberg ODE solver of order (4)5;
  • rodas23, L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (2)3;
  • rodas34, L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (3)4;
  • Rosenbrock12, L-stable embedded Rosenbrock ODE solver of order (1)2;
  • Rosenbrock23, L-stable embedded Rosenbrock ODE solver of order (3)4;
  • Rosenbrock34, L-stable embedded Rosenbrock ODE solver of order (3)4;
  • seulex, extrapolation-algorithm, based on the linearly implicit Euler method;
  • SIBS, semi-implicit Burlish-Stoer (deprecated).

For very stiff systems, in particular complex chemistry, the seulex method has proved to be particularly robust and efficient, followed closely by rodas23 which is more efficient when low accuracy is sufficient. For less stiff cases the Rosenbrock43 may be more efficient than either seulex or rodas23 and for non-stiff systems a range of Runge-Kutta solvers are available.

Source code

  • ODE library
    $FOAM_SRC/ODE

Symmetry Plane handling

In this version, the behaviour of the symmetryPlane boundary condition has changed, in which it now explicitly requires the underlying patch to be a single, flat plane. The change has been made to avoid accumulation of truncation errors and to guarantee parallel consistency in any dealing with point fields, e.g. mesh motion.

The behaviour from previous versions, in which it was not an absolute requirement that the patch is flat, is now available in a new symmetry condition. However, for any case dealing with point fields, it is advised to upgrade to use the symmetryPlane boundary condition.