OpenFOAM® v1912: New and updated solvers and physics


New interface-tracking dynamic mesh

The moving mesh interface tracking library has been integrated into this release from the foam-extend project. Porting to OpenFOAM v1912 was performed by Zeljko Tukovic (FMENA, Zagreb) in collaboration with OpenCFD Ltd.

Basic functionality is implemented through the interfaceTrackingFvMesh class. Kinematic and dynamic conditions at the free-surface are imposed by custom boundary conditions for velocity and pressure. Interface tracking calculations are performed using standard pimpleFoam solver.

The example below shows the concentration and shape evolution of a contaminated, perturbed droplet.

Source code
Z. Tukovic and H. Jasak. A moving mesh finite volume interface tracking method for surface tension dominated interfacial fluid flow. Computers & fluids. 55 (2012) 70-84.
Code written and contributed by Zeljko Tukovic

New Six Degree of Freedom (DoF) features

New linearSpringDamper restraint model

The Six DoF library includes a new linearSpringDamper model.

This is a spring-damper system that acts as a ”soft” rope restraint when the distance between the anchor and refAttachmentPt exceeds the restLength.

When the rope is slack there is no force applied to body. The specification in dynamicMeshDict is as follows:

        sixDoFRigidBodyMotionRestraint  linearSpringDamper;
        refAttachmentPt     (3.7 4.2 4.5);
        anchor              table
            (0         (3.7 4.2 6.5))
            (2         (4.7 4.2 6.5))
            (6         (5.7 4.2 4.5))
            (8         (6.7 4.2 4.5))

        psi                 1;
        wn                  8.28;
        numberOfChains      4;
        restLength          2;


  • refAttachmentPt is the attachment point to the solid object.
  • anchor is the anchor of the rope, specified as a Function1 time dependant entry.
  • psi, ψ  \relax \special {t4ht= is the damping ratio.
  • wn, ωn  \relax \special {t4ht= is the un damped natural frequency of the system.
  • numberOfChains, n  \relax \special {t4ht= is the number of chains attached to the body.
  • restLength is the rest length of the rope.

The damping and stiffness of the system is calculated as:

pict\relax \special {t4ht=

Note: psi and wn should be set up for an over-damped system to model the behaviour of a soft rope.

Source code

New soft wall restraints

The new soft wall model is a damper-linear-spring restraint that acts as a ”soft” wall when the distance between anchor and refAttachmentPt in the wall-normal direction becomes negative.

When the rope is slack there is no force applied to body. The specification in dynamicMeshDict is as follows:

        type                    softWall;
        body                    floatingObject;
        anchor                  (0.5 0.5 0.7);
        refAttachmentPt         (0.5 0.5 0.58);
        wallNormal              (0 0 -1);
        psi                     2.0;
        C                       0.01;


  • refAttachmentPt is the attachment point to the solid object
  • anchor is the anchor of the rope.
  • psi is the damping ratio
  • wallNormal is the normal wall direction
  • C is a model constant

The damping and stiffness of the system is calculated as:

pict\relax \special {t4ht=

Source code

New multi-motion capabilities for overset dynamic meshes

Overset dynamic meshes can now handle multiple motion types. This allows users to have different objects moving under the action of different motion solvers.

As an example see:

In the dynamicMeshDict dictionary the dynamicOversetFvMesh has two motion solvers:

  • solidBody using a new drivenLinearMotion (see below)
  • sixDoFRigidBodyMotion

In the tutorial two motion solvers are applied to different cell sets.

The solver named cube is the floatingObject which uses the motion solver sixDoFRigidBodyMotion. A new keyword, CofGvelocity writes the centre of gravity velocity of the object for use in the drivenLinearMotion solver.

The new drivenLinearMotion motion solver is a rigid body type that moves points according the entry CofGvelocity. In this tutorial it is applied to the c0 cell set which corresponds to the background mesh in the overset calculation.

The drivenLinearMotion enables users to set a background mesh that ’follows’ a given object in the ’inset’ region of the overset calculation, this helps to reducing mesh size when the expected displacement of the object is large.


New collimated radiation models

New collimated solar beam model for fvDOM

The fvDOM radiation model now includes support for external beams from a solar external load to provide a straightforward means to include solar load in the domain. The effect is added to the rays in the fvDOM model, activated by setting the new entry useExternalBeam to true in the radiationProperties dictionary.

The fvDOM ray solid angles are rotated to match the best ray direction according to the angular discretization.

To obtain the energy spectra of the source the following entry is used:

spectralDistribution (2 1)

Note: The number of bands in the spectralDistribution should match the number of bands of the absorptivity model. i.e.multiBandAbsorption in the boundaryRadiationProperties:

type       opaqueDiffusive;
    type            multiBandAbsorption;
    absorptivity    (0.3 0.7);
    emissivity      (0.3 0.7);

The radiative boundary condition excample specifies an opaqueDiffusive wall and a multiBandAbsorption model using two bands.

The collimated model should NOT be used in conjunction with the solarLoad radiation model.

The collimated model uses the solar calculator to derive the sun direction and energy flux. The standard models can be used as for the solarLoad radiation model.

Source code

New multi band zone absorption-emission model

The new multiBandZoneAbsorptionEmission model enables the use of non-uniform absorptivity and emissivity fields in the domain. For example, for two bands model:

    absorptivity  (0.01 0.01);
    emissivity    (0.01 0.01);

        absorptivity  ("trees" (20 20));
        emissivity    ("trees" (20 20));

This example sets an absorptivity value of (0.01 0.01) in the background mesh and (20 20) in the zone named trees.

Source code

New multiphase turbulence stabilization

Turbulence kinetic energy is over-predicted in incompressible VOF solvers at the phase interface and throughout the water column in nearly-potential flow regions beneath surface waves.

This behaviour is corrected using the new multiphaseStabilizedTurbulence fvOption that applies corrections to the turbulence kinetic energy equation and turbulence viscosity field.

It is specified using:

    type            multiphaseStabilizedTurbulence;
    active          yes;

        // Optional coefficients
        lambda2         0.1;   // A value of 0 sets the nut correction to 0
        Cmu             0.09;  // from k-epsilon model
        C               1.51;  // model coefficient from k-omega model
        alpha           1.36;  // 1/Prt

Source code
  • Devolder, B., Rauwoens, P., and Troch, P. (2017). Application of a buoyancy-modified k-w SST turbulence model to simulate wave run-up around a monopile subjected to regular waves using OpenFOAM. Coastal Engineering, 125, 81-94.
  • Larsen, B.E. and Fuhrman, D.R. (2018). On the over-production of turbulence beneath surface waves in Reynolds-averaged Navier-Stokes models J. Fluid Mech, 853, 419-460
  • Implemented by OpenCFD
  • Thanks go to the Turbulence Technical Committee, and the useful discussions with and code testing by Bjarke Eltard-Larsen and David Fuhrman (Technical University of Denmark)

New k-epsilon-phit-f turbulence model

A new three transport-equation linear-eddy-viscosity RANS turbulence closure model with an elliptic relaxation function, named kEpsilonPhitF, has been implemented based on the work of (Laurence et al., 2005).

kEpsilonPhitF is a kEpsilon-based model which originated from (Durbin, 1995)’s v2-f methodology. The main benefit of the v2-f methodology was reportedly to provide a mid-way between the linear eddy-viscosity models and Reynolds stress equation models whilst keeping the minimalism of an eddy-viscosity model formulation. This manifests itself in its inherent and simulation-independent reproduction capability for the parabolic decay of the wall-normal velocity scale towards the wall without using any wall-distance or low-Reynolds number type modifications.

However, the majority of v2-f model variants proved to be numerically stiff for segregated solution algorithms due to the coupled formulations of v2 and f fields, particularly on wall boundaries.

Despite this drawback, the v2-f models have been observed to successfully reproduce the properties of a broad set of generic flows, such as flows exhibiting various non-equilibrium/non-local effects; involving strong pressure gradients, separations, impingements; effects of fluctuating pressure on turbulence fields, including heat transfer and more.

The v2-f variant (i.e. OpenFOAM’s v2f model) due to (Lien and Kalitzin, 2001) reformulated the original v2-f model to enable segregated computations; however, a number of shortcomings regarding the model fidelity were reported in the literature.

To overcome the shortcomings of the v2-f methodology, the v2-f approach was re-evaluated by (Hanjali et al., 2004), and (Laurence et al., 2005) at the same time by transforming   2
v   \relax \special {t4ht= scale into its equivalent non-dimensional form, i.e.       2
ϕ =  v ∕k  \relax \special {t4ht=, to reduce the numerical stiffness. These variants were believed to provide numerical robustness, and insensitivity to grid anomalies while retaining the theoretical model fidelity of the original v2-f model.

Please note that the kEpsilonPhitF model is in the beta stage. Further internal changes to the model may or may not be performed in the future. In addition, the v2f model available in OpenFOAM will be deprecated in favour of kEpsilonPhitF or another v2-f model variant.

The implementation of kEpsilonPhitF model was verified by means of the following canonical flows: the smooth-wall plane channel flows (                                           6
Re τ = 180,395,590, 800,5186; Re = 80 × 10   \relax \special {t4ht=), zero-pressure gradient turbulent flat plate, two-dimensional backward-facing step, and two dimensional hill in a channel.

As an illustration, three verification plots were shown in the following:

Pressure coefficient profile from the two-dimensional hill flow
Pressure coefficient profile from the two-dimensional hill flow

Wall-unit normalised streamwise flow speed profiles
Wall-unit normalised streamwise flow speed profiles from the smooth-wall plane channel flow at Reτ  \relax \special {t4ht==5186.

Wall-unit normalised Reynolds stress profiles
Wall-unit normalised Reynolds stress profiles from the smooth-wall plane channel flow at Reτ  \relax \special {t4ht==5186.

More detailed information regarding the model equations, implementation, and usage can be found in the (kEpsilonPhitF section in the Extended Code Guide).

Source code
OpenCFD would like to thank Dr. Mirza Popovac (Austrian Institute of Technology) for useful discussions and his helpful suggestions.

Updated avalanche modelling module

With winter approaching, the avalanche modelling module received some updates to the entrainment models, including an implementation of the RAMMS model (Christen et al, 2010) and stability improvements for the Medina model.

Source code
The avalanche code is developed by Matti Rauter

New equation of state for liquids

A reciprocal polynomial equation of state for liquids and solids, rPolynomial, has been integrated from, having the following form:

pict\relax \special {t4ht=

The new option is specified in the thermophysicalProperties dictionary as:

    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState rPolynomial; // new entry
    specie          specie;
    energy          sensibleInternalEnergy;

        molWeight   18.0;
        // Coefficients for the reciprocal polynomial equation of state
        C (0.001278 -2.1055e-06 3.9689e-09 4.3772e-13 -2.0225e-16);
        Cp          4195;
        Hf          0;
        mu          3.645e-4;
        Pr          2.289;

Source code
The model was integrated from