OpenFOAM® v1706: New and improved numerics

New overset mesh functionality

OpenFOAM v1706 includes the first release of a major new set of functionality to enable users to perform overset mesh calculations. The approach performs cell-to-cell mapping between multiple, disconnected mesh regions, to generate a composite domain applicable to both static and dynamic analyses.

The mapping is used to determine which parts of the mesh are:

• solved, calculated;
• interpolated from solution cells, interpolated; and
• not used, holes

The result is a straightforward method to perform complex mesh interactions without the need to employ, e.g. challenging topological changes. An example showing the background mesh (blue) with two sub-meshes that represent rotors is shown below:

Here, the overset analysis takes this 3-part mesh and determines the cell types. For the background mesh, the result is given by:

where the colours correspond to:

• blue: calculated
• white: interpolated
• red: blocked

Adapting solvers to include overset capabilities is straightforward for any dynamic mesh solver (”DyM”) incorporating Finite Volume discretisation. Currently options include:

• overLaplacianDyMFoam : laplacianFoam (moving mesh)
• overSimpleFoam : simpleFoam (steady-mesh)
• overPimpleDyMFoam : pimpleDyMFoam (moving mesh)
• overRhoPimpleDyMFoam : rhoPimpleDyMFoam (moving mesh)
• overInterDyMFoam : interDyMFoam (moving mesh)

Case set-up

Updating a moving mesh case for overset requires the following changes:

• The case must have at least one patch of type overset. This should ideally be the first patch. As a constraint type, it will automatically assume a boundary condition type overset by including the default set-up file in the boundaryField section of the field files:

  #includeEtc "caseDicts/setConstraintTypes"
  
• Any cells on the outside of the mesh (part) that require interpolation should be on this patch.
• Set the mesh type to dynamicOversetFvMesh in the dynamicMeshDict (located in the constant directory)

  dynamicFvMesh       dynamicOversetFvMesh;
  
  solver              displacementLaplacian;
  
  displacementLaplacianCoeffs
  {
      diffusivity     uniform 1;
  }
  
• Generate a 0/zoneID volScalarField indicating the mesh zone (not to be confused with cellZones). This field can usually be generated using standard OpenFOAM tools, e.g.topoSet, setFields. It should start at 0 and be consecutively numbered. There is no limit to the number of zones.
• In the fvSchemes dictionary (in the system directory) select the relevant overset interpolation method and which additional variables should be included. By default, all solved-for variables and some solver-specific fields use overset interpolation.

  oversetInterpolation
  {
      method          inverseDistance;
  }
  
  oversetInterpolationRequired
  {
      nut;
  }
  
• If necessary adapt the wall distance calculation to a continuous method, e.g.Poisson, advectionDiffusion. The default meshWave method does not walk through the overset interpolation.
• Update the solver to an asymmetric variant in the fvSolution dictionary (the matrix from interpolation is asymmetric). All asymmetric solvers except from GAMG are supported. In practice a good choice is the smoothSolver with symGaussSeidel for transport equations (U, k, etc.) and PBiCGStab with DILU preconditioner for elliptic equations (p, yPsi).

Running

All overset solvers can run in parallel using the normal tools, e.g.decomposePar and mpirun. No special set up is required; however, overset interpolation benefits when the participating cells reside on the same processor as the interpolation becomes more implicit. Additional output is provided Compared to the normal dynamic mesh solvers related to the chosen oversetInterpolation method:

This output shows that there are two mesh ’zones’, i.e. described by the 0/zoneID field. Of the cells, some are identified as unreachable (hole), some are interpolated, and the majority are solved (calculated).

Additional fields are written for post-processing:

• zoneID : a copy of the input field
• cellTypes : result of the hole detection analysis:
  • 0 : cell is calculated
  • 1 : cell is interpolated
  • 2 : cell is hole/inactive

The cellTypes field can be used to subset only the active set of cells — in ParaFoam: Filter->Threshold

Further information

Source code

$FOAM_SOLVERS/basic/overLaplacianDyMFoam
$FOAM_SOLVERS/incompressible/simpleFoam/overSimpleFoam
$FOAM_SOLVERS/incompressible/pimpleFoam/pimpleDyMFoam
$FOAM_SOLVERS/multiphase/interFoam/overInterDyMFoam
$FOAM_SOLVERS/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam

Examples

$FOAM_TUTORIALS/basic/overLaplacianDyMFoam/heatTransfer
$FOAM_TUTORIALS/incompressible/overSimpleFoam/aeroFoil
$FOAM_TUTORIALS/incompressible/overPimpleDyMFoam/cylinder
$FOAM_TUTORIALS/incompressible/overPimpleDyMFoam/simpleRotor
$FOAM_TUTORIALS/incompressible/overPimpleDyMFoam/twoSimpleRotors
$FOAM_TUTORIALS/multiphase/overInterDyMFoam/floatingBody

New interface capturing functionality - iso-advector

The new IsoAdvector algorithm captures the fluid interface in multiphase calculations. It is currently avaiable in a variant of the interFoam solver, named interIsoFoam, offering an alternative to the MULES interface capturing method. IsoAdvector employs a geometric Volume-of-Fluid (VOF) method for the advection of a sharp interface between two immiscible, incompressible fluids, applicable to both structured and unstructured meshes.

Tests have shown that the method generally offers more accurate interface advection and a sharper interface representation than the existing MULES method. The following videos compare the results obtained for the interIsoFoam and interFoam solvers.

Pure advection tests

Note that the time spent in the advection step is significantly reduced. Depending on the case and the set-up, the sharper interface can cause both reduction or increase in overall simulation time due to its effects on the pressure-velocity system.

Usage

IsoAdvector controls are set in the phase fraction solver section of the fvSolution dictionary, e.g. for the field alpha.water. Many parameters have default settings which are appropriate for most cases. The complete set of options includes:

Note: to operate the interIsoFoam and interFoam in pure advection mode (without solving the pressure, velocity, phase fraction system) the frozenFlow entry should be activated in the PIMPLE sub-dictionary of the fvSolution file, i.e.

Attribution

• The isoAdvector method was developed by Dr. Johan Roenby, DHI, Associate Prof. Henrik Bredmose at DTU Wind Energy and Prof. Hrvoje Jasak at University of Zagreb, Department Faculty of Mechanical Engineering and Naval Architecture.
• The work was funded by a Sapere Aude postdoc from The Danish Council for Independent Research, Technology and Production Sciences (Grant-ID: DFF - 1337-00118B - FTP). Co-funding is also provided by the GTS grant to DHI from the Danish Agency for Science, Technology and Innovation.
• We also acknowledge the significant code contributors by Dr. Vuko Vukcevik and Dr. Andrew Heather.

Limitations

• IsoAdvector is explicit in nature and therefore operates best with a Courant number less than 1. However, sub-cycling can be used to run with larger time steps.
• Moving mesh capabilities will be added in the near future
• Due to the sharper interface, surface tension dominated flows may benefit from a smoothed version of the phase fraction field for improved robustness of the curvature calculation.

Reference

Roenby J, Bredmose H, Jasak H. 2016 A computational method for sharp interface advection. R. Soc. open sci. 3: 160405. http://dx.doi.org/10.1098/rsos.160405

Source code

$FOAM_SRC/finiteVolume/fvMatrices/solvers/isoAdvector
$FOAM_SOLVERS/multiphase/interIsoFoam

Examples

$FOAM_TUTORIALS/multiphase/interIsoFoam/damBreak
$FOAM_TUTORIALS/multiphase/interIsoFoam/discInReversedVortexFlow
$FOAM_TUTORIALS/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation
$FOAM_TUTORIALS/multiphase/interIsoFoam/sphereInReversedVortexFlow
$FOAM_TUTORIALS/multiphase/interIsoFoam/standingWave

Improved second order restart

The restart behaviour in older versions of OpenFOAM often exhibited noticeable instabilities, particularly for compressible flows, and using 2nd order time integration.

This version includes many incremental updates to aid smooth restart behaviour, leading to significant improvements in restart robustness.