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:
Results from cell type analysis
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)
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:
- 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)
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.
- 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).
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:
Overset analysis : nCells : 660
calculated : 532
interpolated : 98
hole : 30
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
- Source code
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.
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:
// Only cells with surfCellTol < alpha < 1 - surfCellTol are treated as
// surface cells
// Number of times the ad-hoc bounding step should try to correct
// unboundedness. Strictly volume conserving (provided that sum(phi) = 0 for
// a cell).
// Optional: cells with alpha < snapAlphaTol are snapped to 0 and
// cells with 1 - alpha < snapAlphaTol are snapped to 1
// Optional: clip remaining unboundedness
// Optional: PLIC-like surface reconstruction using a (smoothed) grad(alpha)
// vector for the surface orientation inside surface cells
// Optional: write isofaces to <case>/isoFaces/isoFaces_[timeIndex].obj
// Optional: write a surface cellSet at each time step
// Optional: write a cellSet of cells which required heuristic bounding
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.
- 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.
- 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.
- 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
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.