OpenFOAM® v1806: New and updated solvers and physics


This release includes a major new suite of functionality that builds upon the Volume Of Fluid (VOF) method by adding new heat and mass transfer capabilities, with an example application that can be used to approximate laser melting

The animation above describes a moving laser heat source applied to an array of aluminium spheres, created using the new ParaView Catalyst adapter

icoReactingMultiphaseInterFoam solver

icoReactingMultiphaseInterFoam is a multi-phase, multi-component incompressible solver based on a Volume Of Fluid (VOF) method with per-phase choice of thermodynamics model (sharing pressure and temperature). In addition it supports mass and heat transfer across phases and optional two-way radiation interaction via the Discrete Transfer Radiation Model (DTRM) (see below), including optional reflection at phase interfaces.

The basic set up is inside the new phaseProperties dictionary which defines

  • the names of the individual phases,
  • per phase the model to evolve it,
  • the mass transfer model for the combination of two phases,
  • optional: a porosity model for the combination of two phases,
  • optional: surface tension model for the combination of two phases.

Phase model types

Options for the model for each phase include:

  • pureMovingPhaseModel : pure moving phase, e.g. a fluid
  • pureStaticSolidPhaseModel : pure static phase, e.g. a solid
  • multiComponentMovingPhaseModel : multi-component moving phase, e.g. a multi-component fluid

An example of an evaporating liquid simulation could for example specify the two phases, one the pure liquid and the other a mixture of air and evaporated liquid:

phases  (liquid gas);

    type            pureMovingPhaseModel;

    type            multiComponentMovingPhaseModel;

Each phase has its own thermophysical model described by a thermophysicalProperties.phaseName dictionary, where the user can select from the available run-time selectable thermophysical model packages.

Mass transfer models

The mass transfer models are selectable for each pair of phases. The current choices are:

  • Lee : solid melting and solidification
  • kineticGasEvaporation : condensation and evaporation

i) Solid melting/solidification model

Mass transfer between phases is driven by the Lee model, taking the form:

pict\relax \special {t4ht=

where C  \relax \special {t4ht= is a model constant. If C >  0  \relax \special {t4ht=, i.e. melting

pict\relax \special {t4ht=

For C <  0  \relax \special {t4ht=, i.e. solidification:

pict\relax \special {t4ht=

Based on the reference:

  • W. H. Lee. ”A Pressure Iteration Scheme for Two-Phase Modeling”. Technical Report LA-UR 79-975. Los Alamos Scientific Laboratory, Los Alamos, New Mexico. 1979.

A typical specification might be:

    (solid to liquid)
        type            Lee;
        C               130;
        Tactivate       659;

ii) The condensation/evaporation model. (kineticGasEvaporation)

This model is based on the kinetic gas theory, considering the Hertz Knudsen formula, which gives the evaporation-condensation flux based on the kinetic theory for a flat interface:

pict\relax \special {t4ht=

where F  \relax \special {t4ht= represents the mass flux rate [Kg/s/m2], M  \relax \special {t4ht= the molecular weight, T
  activate  \relax \special {t4ht= the activation temperature, C  \relax \special {t4ht= the accommodation coefficient, R  \relax \special {t4ht= the universal gas constant, psat  \relax \special {t4ht= the saturation pressure and p  \relax \special {t4ht= vapour partial pressure

The Clapeyron-Clausius equation relates the pressure to the temperature for the saturation condition:

pict\relax \special {t4ht=

Where L  \relax \special {t4ht= is the latent heat, nuv  \relax \special {t4ht= is the inverse of the vapor density and nul  \relax \special {t4ht= the inverse of the liquid density

Using the above relations:

pict\relax \special {t4ht=

This assumes liquid and vapour are in equilibrium, whereby the accommodation coefficients are equivalent for the interface. This relation is known as the Hertz-Knudsen-Schrage.

Based on the reference:

  • Van P. Carey, Liquid-Vapor Phase Change Phenomena, ISBN 0-89116836, 1992, pp. 112-121.

A typical specification might be:

    (liquid to gas)
        type            kineticGasEvaporation;
        species         vapour.gas;
        C               0.1;
        alphaMin        0.0;
        alphaMax        0.2;
        Tactivate       373;

where :

  • C  \relax \special {t4ht= is the model constant. If C  > 0  \relax \special {t4ht= the model acts as evaporation, if C <  0  \relax \special {t4ht= is condensing.
  • alphaMin and alphaMax are the range of the liquid on which the model is applied on the interface.
  • Tactivate is the saturation/condensation temperature
  • species is the name of the component of the gas being transferred. If the ’gas’ phase is pure (no components) then the entry species does not need to be specified.

Inter-phase porosity models

These models can add an artificial momentum source on/next to the interface between the two phases. A good example is to influence the behaviour when solidifying/melting.

The only porous interface method is VollerPrakash, based on the references:

  • V.R. Voller and C. Prakash, A fixed grid numerical modelling methodology for convection-diffusion mushy phase-change problems, Int. J. Heat Mass Transfer 30(8):17091719, 1987.
  • C.R. Swaminathan. and V.R. Voller, A general enthalpy model for modeling solidification processes, Metallurgical Transactions 23B:651664, 1992.

It can be specified as:

    (solid and liquid)
        type            VollerPrakash;
        solidPhase      alpha.solid;
        Cu              1e10;

where solidPhase is the phase representing the solid phase and Cu is a model constant.

Inter-phase surface tension models

The surface tension can be specified on a phase pair basis. The only current choice is the constant model:

    (gas and liquid)
        type            constant;
        sigma           0.07;

Source code

DTRM radiation model

A new radiation model is introduced to enable simulation of collimated radiation flux. It uses the Discrete Transfer Radiation Model (DTRM) to enable ful two-way interaction between radiation and participating media. It is currently specialised to represent a collimated laser beam passing through a VOF system.

Laser specification

The DTRM computes radiation effects through a set of rays, passing through the mesh, i.e. equivalent to Lagrangian particle tracking where the laser represents the particle injector. The specification for the laser includes:

  • Laser power, shape, distribution
  • Laser position, orientation
  • Laser DTRM resolution

The laserPower (W) is a Function1 Type of entry to allow time based variation. The laser power distribution can be uniform, Gaussian or manual. When using the mode manual the total laser power output is :

pict\relax \special {t4ht=

allowing control over the time and space distribution.

The focalLaserPosition describes the current centre of the circular injection disk, laserDirection describes the direction of the rays. Both are of type Function1 to enable time dependency.

The DTRM resolution is specified as a polar division of a 2D circular disk:

nTheta      70;
nr          30;
where nTheta and nr are the number of particles in the theta and radial directions respectively.

Phase interaction (absorbtion, emissivity)

The model is developed to work in a VOF framework, representing different phases.The absorption and emission is calculated from the individual phases’ absorption, emission through a simple mixture model:

pict\relax \special {t4ht=

pict\relax \special {t4ht=

where α
 N  \relax \special {t4ht= is the phase volume fractions

Inter-phase interaction (reflection)

The only reflection model supplied is the FresnelLaser model. This model is based on :

  • Implementation of real-time multiple reflection and Fresnel absorption of FresnelLaser beam in keyhole. J. Phys. D: Appl. Phys. 39 (2006) 5372-5378.

The reflection model is specified per pair of phases:

(gas and liquid)
    type            FresnelLaser;
    epsilon         0.25;

This model will attempt to find the interface between the two phases and the particles coming from the gas to the liquid will be reflected back to the gas and transmitted on to the liquid.

It is possible to specify several reflection models for different interfaces. This will increase considerably the execution time as new particles are inserted for each reflection until almost all energy is absorbed or lost through the outlet patches.

The interface is located using a user-supplied threshold of the phase fraction taking a default value of 0.5  \relax \special {t4ht=.

Consolidation of moving mesh solvers

Moving mesh functionality has been incorporated into many of the static mesh solver applications from earlier releases:

Old solver New solver
pimpleDyMFoam pimpleFoam
rhoPimpleDyMFoam rhoPimpleFoam
interDyMFoam interFoam
multiphaseInterDyMFoam multiphaseInterFoam

These updates were performed by

Extended interIsoFoam solver for dynamic meshes

The interIsoFoam solver and its core isoAdvector library have been extended to work with dynamic meshes.

In particular, the code now works with the dynamicRefineFvMesh class for automatic mesh refinement (AMR). The damBreakWithObstacle case using this functionality has been added to the interIsoFoam tutorials:

AMR functionality is controlled in the dynamicMeshDict e.g.


Also the code now works with dynamicMotionSolverMesh the solid body mesh motion solver. The sloshingTank2D case using this functionality has been added to the interIsoFoam tutorials:

Source code

Known issues
Warnings from Foam::isoCutFace::cutPoints and Foam::isoCutFace::quadAreaCoeffs may appear in the log file due to occurrences of multiple cut points along a mesh edge. This will be investigated further, but does not appear to cause visual effects on the solution quality.
Note: Although interIsoFoam can now also be run with other moving mesh classes, the solver has not yet been modified to work consistently with cells of changing shape and volume. The aim is to include this in v1812.

This extension was provided by Johan Roenby

Angular velocity restraint for solid-body simulations

The Rigid-body-dynamics (RBD) framework, introduced in v1606+, has been extended with the capability to specify the motion of an object relative to its parent. This specification is through a new restraint, prescribedRotation, that specifies the relative (angular) velocity with respect to its parent. This velocity is of the generic Function1 type so can be e.g. a single constant, an in-line table or read from a file csvFile, all with interpolation in time. In below example the propeller rotates with an increasing rotational velocity.

    type                    prescribedRotation;
    body                    propeller;
    referenceOrientation    (1 0 0 0 1 0 0 0 1);
    axis                    (1 0 0);
    omega                   table
        (0 (0 0 0))
        (1 (16 0 0))

Known issues
The restraint uses an iterative procedure to obtain the momentum for the desired rotational velocity so might have problems starting up - hence it is a good idea to ramp the velocity. The other recommendation is to more than the default single nIter overall iteration.

New coded mesh motion

A new displacement type motion solver using the coded framework allows a full specification of the motion through user coding. In the constant/dynamicMeshDict one can now specify coded as the motionSolver:

dynamicFvMesh   dynamicMotionSolverFvMesh;
motionSolver    coded;
name            myMotion;

    const Time& tm = mesh().time();
    const pointField& p0 = points0();

    // Inflate w.r.t. origin
    return p0*tm.value();

In above example the solver will be compiled into sub-directory dynamicCode/myMotion and automatically updated, loaded executed. The animation from the example tutorial below shows a twisting motion implemented using the coded motion solver.

Source code

New wave-mangrove interaction

Wave interaction with mangroves can now be modelled using the new momentum and turbulence finite volume options, based on the reference:

  • Tsunami wave interaction with mangrove forests: A 3-D numerical approach. Maza, M, Lara, J.L., & Losada, I.J. [2015], Coastal Engineering [Vol.98, pp. 33-54]


Source code
These boundary conditions were supplied by the Environmental Hydraulics Institute IHCantabria
The code has been updated by OpenCFD and added to the $FOAM_SRC/waveModels library available to the interFoam family of solvers