OpenFOAM® v2006: New and updated solvers and physics


Improved thermal baffle

The compressible::thermalBaffle temperature boundary condition is used to create a thermally-aware solid region (baffle) adjacent to the fluid region. It provides all the information needed to solve the energy equation in a solid and couple to the fluid.

In previous versions, the thermalBaffleModel could not be used on an external fluid patches. In this release users can now extrude from a patch, creating a conjugate heat transfer (CHT) case using any of the heat transfer solvers.

The condition is specified in the temperature, T field boundary, and requires:

  • solid thermoType model
  • solid mixture type
  • radiation model
  • mesh extrusion Model

The following figure shows the main fluid and solid regions. Here, the solid region was extruded from the top patch of the fluid.


type               compressible::thermalBaffle;

// Solid baffle region name
region              baffle3DRegion;

Tnbr               T;
kappaMethod        fluidThermo;

// Is this baffle internal
internal            true;

// Solid thermo
    type            heSolidThermo;
    mixture         pureMixture;
    transport       constIso;
    thermo          hConst;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleEnthalpy;

// Solid mixture
        molWeight       20;
        kappa           0.01;
        Hf              0;
        Cp              15;
        rho             80;

// Radiation Model
    radiationModel  opaqueSolid;
    absorptionEmissionModel none;
    scatterModel    none;
// Extrude model
extrudeModel        linearNormal;
nLayers             50;
expansionRatio      1;
columnCells         false;
    thickness           0.1;

Note the new keyword entry internal that controls whether the baffle is extruded from this patch is internal or external to the fluid region. When the internal option is false, a coupling solid/T is required named bottom which must also be of the type compressible::thermalBaffle.

Source code

Improved handling for particle injection bounds

Users now have the flexibility to state whether particle injections located outside of the domain should simply be ignored as opposed to generating a fatal error. This can be useful, e.g. when using experimental datasets to initialise particle locations where the domain geometry may not fully correspond to the computational set-up.

The option is included when specifying the particle injection model, e.g.

    type            reactingMultiphaseLookupTableInjection;

    // New entry to allow out of bounds
    ignoreOutOfBounds yes;

If activated the number of injection locations identified as out-of-bounds is reported in the solver log.

Source code

New Weber number cloud function object

The new WeberNumber cloud function object computes and writes the particle Weber number to disk at simulation write times. The field is written as a ’standard’ cloud field in <time>/lagrangian/<cloudName>/We available for further post-processing, e.g. using ParaView.

The cloud function object is specified in the constant/<cloudName>Properties<i> file as:

        type    WeberNumber;

The following example is taken from the aachenBomb tutorial case where the particles are coloured according to the Weber number and the carrier by the veocity magnitude.


Source code

Improved Curle function object

The Curle function object released in OpenFOAM v1706 now includes the near field term and has been refactored to compute the acoustic pressure for a set of points based on a user-defined list or surface.

When using point mode, the entries include a list of observer positions, e.g.

    type            Curle;
    libs            (fieldFunctionObjects);
    patches         (surface1 surface2);
    c0              330;

    input           points;
        (0 0 0)
        (1 0 0)
        (x y z)

For surface mode the face centres of the input surface imply the observer locations. The output can then be written back to a surface format or a plain text file:

    type            Curle;
    libs            ("");
    patches         (surface1 surface2);
    c0              330;

    input           surface;
    surface         "inputSurface.obj"

    // Output - either points or surface
    output          points;
    //output          surface;
    //surfaceType     ensight;

The image below shows an example from the vortexShed tutorial case, where the surface geometry is shown in white on the left (overlaying a velocity magnitude prediction); the right image shows the acoustic pressure computed on the input surface.


Both point and surface ouputs are supported by the noise utility, and can be used to e.g. highlight the vortex shedding frequency, here shown to be approximately 0.1 Hz.


Source code

Improved multiphase phase change modelling

The phase change incompressible solvers interCondensatingEvaporatingFoam and icoReactingMultiphaseInterFoam have been improved:

  • New interfaceHeatResistance phase change Model
  • Semi-implicit treatment of the temperature equation
  • Use of the isoCutCell method for interface area density
  • Use of dynamic mesh

The new interfaceHeatResistance phase change model uses a spread source distribution technique (see References):

The phase change source is spread across the sharp interface in order to smooth the mass transfer rate. The mass transfer rate is calculated as:

pict\relax \special {t4ht=


  • Ai  \relax \special {t4ht= is the interface area density
  • Tsat  \relax \special {t4ht= saturation temperature
  • L  \relax \special {t4ht= latent heat

In this model the location of the isosurface can be controlled by the entry isoAlpha.

Additionally, specifically to the solver icoReactingMultiphaseInterFoam, new options are now available:

In order to give more flexibility to the phase change models new options are now available to include or neglect volumetric change due to phase change (includeVolChange)

In the solver icoReactingMultiphaseInterFoam the mass exchange model kineticGasEvaporation was made more robust.

The usage of these models can be seen on the following links:

Source code
Hardt, S., Wondra, F. (2008). Evaporation model for interfacial flows based on a continuum- field representation of the source term Journal of Computational Physics 227 (2008), 5871-5895
The interfaceHeatResistance implementation was a contribution by Henning Scheufler.

New and improved isoAdvector-based modelling

In this version the isoAdvector method has been updated and new sub-models and methods added:

  • New reconstruction methods: isoAlpha, plicRDF and gradAlpha
  • Compressible formulation of the interIsoFoam solver, compressibleInterIsoFoam
  • An improved alpha interpolation scheme for adaptive mesh refinement
  • A post processing function for piecewise linear interface construction (PLIC) interfaces based on sampledSurfaces

The most significant updates concern the newly developed reconstruction scheme (see References).

Tests have shown that the isoAlpha reconstruction method does not perform well on unstructured meshes. The newly developed method, plicRDF, improves the normal orientation calculation and is a better choice when using unstructured meshes.

The new top-level solver compressibleInterIsoFoam provides support for the isoAdvector method for compressible flows. This solver and its incompressible counterpart interIsoFoam both support Automated Mesh Refinement (AMR).


The new reconstructionScheme keyword is specified in the fvSolution file:

    nAlphaCorr      1;
    nAlphaSubCycles 1;
    cAlpha          1;

    reconstructionScheme plicRDF;
    vof2IsoTol      1e-8;
    surfCellTol     1e-6;
    nAlphaBounds    3;
    snapTol         1e-12;
    clip            true;

The following video shows a comparison of the results obtained for the depthCharge3D case for the new compressibleInterIsoFoam using the isoAdvector method and compressibleInterFoam using MULES.

Source code
Henning Scheufler and Johan Roenby. ”Accurate and efficient surface reconstruction from volume fraction data on general meshes.” Journal of Computational Physics (2019), doi: 10.1016/
This implementation was carried out by Henning Scheufler (DLR) and Johan Roenby (DHI).
The integration was performed by OpenCFD Ltd with contributions from the authors.

New pipelined Conjugate Gradient solvers

When running in parallel any conjugate gradient (CG) solver requires all processors to use the same search direction. This requires a global reduction (MPI_Allreduce). On large numbers of cores this might become a scaling bottleneck.

To get around this two new linear solvers have been introduced which ’pipeline’ the CG loop to

  • replace three global reductions of a scalar with a single one for three scalars (reducing the overhead)
  • allow these global reductions to overlap with ’halo-swaps’ (communication to neighbouring processors only)

On large numbers of cores this might help the scaling. The disadvantage of the methodology is that it

  • increases the number of (local) operations
  • increases the memory needed

Pipelined Preconditioned Conjugate Gradient solver (PPCG)

The first of the new solvers recasts the existing PCG solver in the pipelined framework. It should behave exactly the same as the existing PCG solver, apart from truncation error differences. For simple cases the residual and number of iterations should be identical.

Pipelined Preconditioned Conjugate Residual solver (PPCR)

This is a modification of the PPCG solver which bases the elimination direction on the preconditioned variables. It shows slightly better residual behaviour at the first iterations but at the cost of less communication overlap since now the precondition has to be finished before starting the global reductions. For details see above reference. This solver is interesting to use for solving the coarsest level of the GAMG solver with a loose tolerance, e.g. from the incompressible/pisoFoam/LES/motorBike tutorial:

    relTol          0;

    // Explicit specify solver for coarse-level correction to override
    // solution tolerance
        solver          PPCR;
        preconditioner  DIC;
        relTol          0.05;

Purely looking at the coarsest level solution (overall number of sweeps needed to solve the coarsest level):

Coarse level number of sweeps

The PPCR solver requires less number of sweeps even though occasionally it requires more top-level GAMG cycles so more invocations of the coarse level solver. Overall statistics:

Solver Coarse level sweeps Run time (s)
(P)PCG 7539 405
PPCR 6324 401


  • simulation has only been run once so the timings are not representative
  • simulation was run on 8 cores of the same node - no interconnect used
  • no simulation statistics were compared - however no changes were made to the case setup apart from the solver choice.

Currently these settings are being tested on more representative cases (32, 64 milllion cells, >1000 cores). Timings so far show an improvement in scalability; we will report our findings later on separately from these release notes.

Source code
Ghysels, P., and Vanroose, W. (2014). Hiding global synchronization latency in the preconditioned conjugate gradient algorithm. Parallel Computing, 40(7), 224-238.

New cloud function object to remove particles

The new removeParcels cloud function object enables users to remove particles that reach predefined faceZones. This can be beneficial when the particles are no longer of interest and would otherwise be tracked until they exited a domain boundary or eliminated by some other means, consuming valuable resources, e.g. for rain simulation perhaps only particles that reach the air filter are of interest; having passed the vehicle they could be discarded.

The removeParcels object can be specified as follows:

        type            removeParcels;
        log             yes;
        resetOnWrite    no;
        resetOnStart    no;
        faceZones       (cycLeftcycRight);

When activated the object reports number and mass of parcels that have been removed per faceZone, e.g.

removeParcels1 output:
    faceZonecycLeft: removed 994 parcels with mass 0.009817407917
    faceZonecycRight: removed 0 parcels with mass 0
The same information is also saved to the file

Source code

Improved lumped point movement for large structures

This release extends the external coupling of structures with OpenFOAM. The lumped point method is devised to handle Fluid-Structure-Interactions for large structures with an FEA model that is much, much coarser than the CFD mesh and movements that are much slower than the fluid domain. For this, the given geometry

Bridge geometry

is defined in by its ”rigging”, which are motion controllers defined in terms of lumped points and connectivity. They can be considered the ”skeleton” for the movement.

Bridge rigging

The motion controllers are associated with mesh patches to define the interpolation influences.

Bridge zones

Moving the lumped points on the rigging causes the patch points to move, which are used in the boundary conditions for the mesh motion solver. Note: both rigging displacement and bridge displacement are shown at exaggerated scales.

Source code

New and updated avalanche and mud slide tools

This release includes many updates to the avalanche community submodule:

  • New physical model (solver) faParkerFukushimaFoam (based on Parker et al. 1986; for the simulation of turbidity currents (in collaboration with Kate Heerema, supported by
  • Improved integration of geospatial data formats enables simple case setup, georeferenced results and post processing in Geographic Information Systems, e.g. QGIS)
  • Improved utility releaseAreaMapping to directly import geospatial vector data (ESRI(R) shapefiles) and geospatial raster data (ESRI(R) girds) to finite area fields
  • New utility gridToSTL for the generation of geometries (STL-files) from geospatial terrain/bathymetry data
  • New function object shapefileWrite for the export of finite area fields as geospatial vector data (ESRI(R) shapefiles)
  • New function object gridfileWrite for the export of finite area fields as geospatial raster data (ESRI(R) gird)

Exported and georeferenced results of the Monterey Canyon in QGIS

Source code
Matthias Rauter, Norwegian Geotechnical Institute