OpenFOAM® v2006: New and improved numerics

OpenFOAM® v2006: New and improved numerics


Improved Arbitrary Mesh Interface (AMI)

The AMI code has undergone significant refactoring and extended to include optional topological changes.

AMI code refactoring

The AMI core has been updated to enable greater control and easier code maintenance, whilst maintaining user-level backwards compatibility. AMI controls are now exposed and can be specified in the constant/polyMesh/boundary file, e.g.

    type            cyclicAMI;
    inGroups        1(cyclicAMI);
    nFaces          2880;
    startFace       276288;
    matchTolerance  0.0001;
    transform       noOrdering;
    neighbourPatch  AMI2;

    // Exposed entries
    AMIMethod       faceAreaWeightAMI;
    restartUncoveredSourceFace yes;

Additional improvements:

  • The directAMI method has been deprecated; the new nearestFaceAMI method offers equivalent functionality at a lower cost.
  • The face matching procedure used by the faceAreaWeightAMI method can now ’walk’ through baffles, providing more complete addressing and weights.

AMI with topological changes

Note: this is beta-level code that has been reasonably tested and shown to be robust; please expect to see new functionality and further improvements in future releases!

This release adds topological change capabilities to the cyclicAMI and cyclicACMI patches, whereby the AMI’s many-to-many addressing is used to construct 1-to-1 face matching between faces on either side of the coupled patches, using the approach described by Aguerre et al. (see References).

Topological changes are activated when selecting the new motion solver in the constant/dynamicFvMeshDict file:

dynamicFvMesh   dynamicMotionSolverFvMeshAMI;
Note: this may be absorbed into the ’standard’ dynamicMotionSolverFvMesh model in future releases.

The effect is to ensure conservation across the AMI patches, and to remove numerical noise in the pressure leading to improved force and moment predictions:


Known issues

  • restart is not currently supported for topological change.

Source code
H.J. Aguerre, S. Marquez Damian J.M. Gimenez, N.M.Nigro, Conservative handling of arbitrary non-conformal interfaces using an efficient supermesh, Journal of Computational Physics 335(15)21-49. 2017.
Many thanks to Horacio Aguerre and Santiago Marquez for bringing the methodology to our attention and for many useful discussions and testing throughout these developments.

New OpenQBMM community module

The OpenQBMM module comprises quadrature-based moment methods for population balance modeling and multiphase flows. OpenQBMM implements methods to describe polydisperse multiphase flows, from the simplest size evolution of non-inertial particles in a carrier fluid to more complex cases with bubbles in a gas-liquid system and inertial particles in a gas-particle flow. These capabilities are implemented in a suite of solvers, among which:

  • pbeFoam solves a population balance equation in a single control volume. This solver is useful to test kernel functions in the population balance equation and to solve spatially homogeneous problems. Example cases can be found in OpenQBMM/validation/pbeFoam, which shows the validation cases discussed in E. Madadi-Kandjani, A. Passalacqua, An extended quadrature-based moment method with log-normal kernel density functions, Chemical Engineering Science. 131 (2015) 323339.
  • pbeTransportFoam allows a frozen flow field to be used to solve a population balance equation with pre-imposed flow motion. A validation case is available in OpenQBMM/validation/pbeTransportFoam/serraTaylorCouette and discussed in A. Passalacqua, F. Laurent, E. Madadi-Kandjani, J.C. Heylmun, R.O. Fox, An open-source quadrature-based population balance solver for OpenFOAM, Chemical Engineering Science. 176 (2018) 306318.
  • buoyantPbePimpleFoam allows a transient flow with the population balance equation to be modeled.
  • polydisperseBubbleFoam is a specialized solver for gas-liquid flows with the evolution of the bubble size due to coalescence and breakup. Bubbles can have a velocity distribution, also accounting for polycelerity, with bubbles with different sizes in the same control volume allowed to have different velocities. Example cases are available in OpenQBMM/validation/polydisperseBubbleFoam while the implementation is discussed in the article J.C. Heylmun, B. Kong, A. Passalacqua, R.O. Fox, A quadrature-based moment method for polydisperse bubbly flows, Computer Physics Communications. 244 (2019) 187204.
  • denseAGFoam implements an anisotropic Gaussian model for dense gas-particle flows. Implementation details are discussed in B. Kong, R.O. Fox, A solution algorithm for fluid-particle flows across all flow regimes, Journal of Computational Physics. 344 (2017) 575594.
  • The solvers in the velocityDistribitionTransport implement quadrature-based moment methods for velocity distributions. These methods are suitable to describe dispersed flows with non-equilibrium velocity distributions, such as crossing jets of particles or droplets. The diluteVdfTransportFoam solver implements the quadrature-based velocity distribution transport algorithm for a disperse phase, not coupled to a carrier fluid. One-way coupling between the disperse phase and the carrier fluid is implemented oneWayCoupledVdfTransportFoam. In both solvers, the particle size can evolve in space and time, and particles can interact through collisions.

Example of numerical simulation of a bubble column with the evolution of the bubble size and reconstruction of the bubble size distribution. The simulation was performed with polydisperseBubbleFoam.

Source code
Validation cases
Alberto Passalacqua and the OpenQBMM community

New dense-matrix eigendecomposition solver: EigenMatrix

The new EigenMatrix solver performs eigendecomposition on diagonalisable, dense, nonsymmetric, and real square matrices.

EigenMatrix decomposes matrices into their canonical forms, whereby a matrix can be represented in terms of its eigenvalues and eigenvectors.

The eigenvalue equation, i.e. eigenvalue problem, is written:

pict\relax \special {t4ht=


  • A  \relax \special {t4ht= : a diagonalisable square matrix of dimension m-by-m
  • v  \relax \special {t4ht= : a (non-zero) vector of dimension m (right eigenvector)
  • λ  \relax \special {t4ht= : a scalar corresponding to v (eigenvalue)

If A  \relax \special {t4ht= is symmetric, the following relation is satisfied:

pict\relax \special {t4ht=


  • D  \relax \special {t4ht= : Diagonal real eigenvalue matrix
  • v  \relax \special {t4ht= : Orthogonal eigenvector matrix

If A  \relax \special {t4ht= is not symmetric, D  \relax \special {t4ht= becomes a block diagonal matrix wherein the real eigenvalues are present on the diagonal within 1-by-1 blocks, and complex eigenvalues within 2-by-2 blocks, i.e. λ + iμ  \relax \special {t4ht= with [λ,μ;− μ, λ]  \relax \special {t4ht=.

The columns of v  \relax \special {t4ht= represent eigenvectors corresponding to eigenvalues, satisfying the eigenvalue equation. Even though eigenvalues of a matrix are unique, eigenvectors of the matrix are not. For the same eigenvalue, the corresponding eigenvector can be real or complex with non-unique entries. In addition, the validity of the equation               T
A  = v ∗ D ∗ v  \relax \special {t4ht= depends on the condition number of v  \relax \special {t4ht=, which can be ill-conditioned, or singular for invalid equations.

The minimum working example of this solver can be seen below:

// A is a m-by-m SquareMatrix<scalar>.
const EigenMatrix<scalar> EM(A);
const DiagonalMatrix<scalar>& realEigenvalues = EM.EValsRe();
const DiagonalMatrix<doubleScalar>& imagEigenvalues = EM.EValsIm();
const SquareMatrix<scalar> eigenvectors(EM.EVecs());
const SquareMatrix<complex> complexEigenvectors(EM.complexEVecs());

Source code
  • This implementation is an integration of the OpenQBMM eigenSolver class (2019) without any changes to its internal mechanisms. Therefore, no differences between EigenMatrix and eigenSolver (2019) classes should be expected in terms of input-process-output operations.
  • The OpenQBMM eigenSolver class derives almost completely from the TNT/JAMA implementation, a public-domain library developed by NIST and MathWorks from 1998 to 2012, available at (Retrieved June 6, 2020). Their implementation was based upon EISPACK.
  • OpenCFD would like to thank the contributors to OpenQBMM, particularly Dr. Alberto Passalacqua (Iowa State University) for the initial implementation of the functionalities.

New and improved adjoint optimisation tools

The adjoint library is enhanced with new functionality enabling the minimisation of an objective function qualitatively quantifying noise through the integral of the squared turbulent viscosity over a portion of the computational domain [1]. This objective function has successfully been used in the past for minimising noise using relatively cheap RANS simulations, in cases where optimisation using high fidelity DES transient simulations is too computationally demanding. Selected figures pertaining to the design of the side mirror of a passenger car showcase the use of this objective function (see references).

The volume used to define the objective function and parameterisation of the side mirror

Turbulent viscosity corresponding to the baseline and optimised mirror

Additionally, the sensitivity derivatives contributions emerging from the differentiation of the positions of the face centres included in the rotatingWallVelocity primal boundary condition are now taken into consideration by introducing its adjoint counterpart, i.e. adjointRotatingWallVelocity. These are an indispensable part of computing the correct sensitivity derivatives when designing rotating components such as the rim of a bicycle wheel.

Shape evolution of the rim of a bicycle wheel targeting the weighted minimisation of the drag and side forces. The geometry is a courtesy of A. Kitselis & T. Theodoridis O.E.

Finally, executing writeActiveDesignVariables before performing an automated shape optimisation loop using adjointOptimisationFoam is no longer necessary, simplifying the optimisation workflow.

E.M. Papoutsis-Kiachagias, N. Magoulas, J. Mueller, C. Othmer, K.C. Giannakoglou: ’Noise Reduction in Car Aerodynamics using a Surrogate Objective Function and the Continuous Adjoint Method with Wall Functions’, Computers & Fluids, 122:223-232, 2015
Source code
The software was developed by PCOpt/NTUA and FOSS GP with contributions from
  • Dr Evangelos Papoutsis-Kiachagias,
  • Professor Kyriakos Giannakoglou,
  • Dr Andrew Heather.
The code has been integrated jointly by OpenCFD and NTUA.
User guide
A PDF guide prepared by NTUA is available here. The user guide has been updated with all features pertaining to an automated shape optimisation loop, including the volumetric B-Splines morpher introduced in OpenFOAM v1912.

New distance-to-patch calculation method

A new variant of meshWave distance-to-patch calculation, i.e. directionalMeshWave has been introduced.

For a given point within a given mesh, the existing meshWave method provides the orthogonal distance to a specified patch. This can be problematic for meshes with very steep terrain, e.g. a hill of 90 deg.

The new directionalMeshWave method ignores the distance component in a user-specified direction, e.g. to enable the distance-to-patch to be computed in the wall-normal direction only.

A comparative example between meshWave and directionalMeshWave output in terms of the wall distance can be seen below:

Wall distance field by meshWave

Wall distance field by directionalMeshWave

The minimum working example of this distance-to-patch calculation method by using system/fvSchemes.wallDist can be seen below:

    method      directionalMeshWave;
    n           (0 0 1);

Source code
Extended code guide

New external solver module

This release provides a new module : external-solver that contains an OpenFOAM interface to the PETSc linear solver library. Building this interface consists of two steps:

  • building the PETSc library
  • building the OpenFOAM interface

One way of building the PETSc library is inside the OpenFOAM ThirdParty directory. This will keep it separate from any system-wide PETSc libraries. A typical build order would be:

# Go into the ThirdParty directory
# Try and see if the PETSc source are there
# If not it will print download instructions, e.g.:
tar xzf petsc-lite-3.13.2.tar.gz
# Again try building:
This will install any PETSc include files and libraries somewhere under the ThirdParty tree. Make a note of the include/ and lib/ directory locations.

The second part of the build is the OpenFOAM solver part. Get hold of the latest external-solver module:

git submodule init
git submodule update
This should create a modules/external-solver tree.

cd modules/external-solver
# Build
./Allwmake -prefix=openfoam
# Update any OpenFOAM settings (not actually necessary for PETSc)
The Allwmake will build the OpenFOAM solver interface to PETSc. Check that the compilation line picks up the correct include directory (specified through -I) and the link line the correct link directory (specified through -L). The link line should end in something like

-o platforms/linux64GccDPInt32Opt/lib/
This library can now be loaded through the system/controlDict:

libs (petscFoam);
Note, however, that you require the underlying PETSc libraries in the LD_LIBRARY_PATH. For our ThirdParty build, this is typically done with the command:

eval $(foamEtcFile -sh -config petsc -- -force)
You can verify that it will be loadable with the following command:

foamHasLibrary -verbose petscFoam
The external-solver module comes with a few tutorials, e.g. the simpleFoam pitzDaily tutorial with PETSc as the linear solver (specified through the system/fvSolution):

Time = 1

Initializing PETSc
PETSc-bicg:  Solving for Ux, Initial residual = 1, Final residual = 0.0165669, No Iterations 1
PETSc-bicg:  Solving for Uy, Initial residual = 1, Final residual = 0.0143104, No Iterations 1
PETSc-cg:  Solving for p, Initial residual = 1, Final residual = 0.0981241, No Iterations 154
time step continuity errors : sum local = 1.71696, global = -0.010057, cumulative = -0.010057
PETSc-bicg:  Solving for epsilon, Initial residual = 0.199669, Final residual = 0.00505847, No Iterations 2
bounding epsilon, min: -1.45735 max: 1080.25 average: 48.855
PETSc-bicg:  Solving for k, Initial residual = 1, Final residual = 0.0131958, No Iterations 3
ExecutionTime = 0.47 s  ClockTime = 1 s

If the petsc4Foam solver is to be built with another version of PETSc just adapt the settings in the (local copy of) etc/ and make sure to source those.
Source code
Developed under Governance: created by the HPC technical committee, involving Mark Olesen, Simone Bna, Stefano Zampini.