OpenFOAM® v2.3.0: Multiphase Modelling
PredictorCorrector SemiImplicit MULES
The success of the volume of fluid (VoF) solvers for multiphase flow in OpenFOAM is underpinned by the development of the multidimensionsal limiter for explicit solution (MULES) as a very effective method of guaranteeing boundedness of scalar fields, in particular phase/massfractions. MULES was introduced in version 1.4 of OpenFOAM and subsequently extended in version 2.1.0 to systems of 3 or more phases and further improved version 2.1.1, with better boundedness and consistency. The effectiveness of MULES comes at a cost, however, since the method is fundamentally explicit which introduces a strict Courant number limit, and hence time step limit, when running the solvers, which is only partially mitigated by the introduction of time step subcycling.
However, in this version of OpenFOAM, a new semiimplicit variant of MULES is introduced which combines operator splitting with application of the MULES limiter to an explicit correction rather than to the complete flux. It first executes an implicit predictor step, based on purely bounded numerical operators, e.g. Euler implicit in time, upwind for convection, etc., before constructing an explicit correction on which the MULES limiter is applied. This approach maintains boundedness and stability at an arbitrarily large Courant number. Accuracy considerations generally dictate that the correction is updated and applied frequently, but the semiimplicit approach is overall substantially faster than the explicit method with its very strict limit on timestep.
The semiimplicit method is supported in all of the interFoam family of incompressible VoF solvers including interPhaseChangeFoam. The controls for MULES, located in the fvSolution file, are now moved out of the PIMPLE controls into the solver controls for the phase fraction field, e.g. alpha.water. The semiimplicit version of MULES can be switched on through the MULESCorr switch, see example below. For large scale VoF simulations, e.g. shipkeeping simulations, the semiimplicit method reduces run times significantly, e.g. in the ship hull case shown below with initial (top) and final (bottom) images of the interface. The simulation of the Duisburg Test Case (DTC) involves motion of the ship (heave and pitch) using the mesh and solid body motion capability in OpenFOAM; careful inspection of the images reveals some pitching (bow down).
{
"alpha.water.*"
{
nAlphaCorr 1;
nAlphaSubCycles 1;
cAlpha 1;
icAlpha 0.25;
MULESCorr yes; // Switches on semiimplicit MULES
nLimiterIter 8; // Number of MULES iterations over the limiter
alphaApplyPrevCorr true;
...
Reference:
el Moctar, O., Shigunov, V., Zorn, T., Duisburg Test Case: PostPanamax
Container Ship for Benchmarking, Journal of Ship Technology Research, 59:5065,
2012.
Phase Property Naming
The naming of phase properties in OpenFOAM has evolved over time to accommodate increasing physics modelling in multiphase flows. The range of physics has reached a point where there is a need to formalise a naming convention that conforms to the policies of code consistency and generality, for sustainable development of the software.
In this version, therefore, multiphase solvers in OpenFOAM use the following conventions:
 alpha.<phase_name> denotes “phase fraction of phase <phase_name>”, e.g. alpha.water is phase fraction of phase water
 <property>.<phase name> denotes “property <property> of phase <phase_name>”, e.g. turbulenceProperties.water are the turbulence properties of phase water
This supports arbitrary numbers of phases in a natural way and separates the specification of the phase physical properties from the ordering of phases in the solvers. The properties of each phase are now specified in separate files, make it much easier to create standard repositories of physical properties and to share data between cases. For example, in the bubbleColumn example for the twoPhaseEulerFoam solver, the constant directory contains the files: phaseProperties; thermophysicalProperties.air; thermophysicalProperties.water; turbulenceProperties.air; and, turbulenceProperties.water. The phaseProperties file contains a phases entry, e.g. see below, that lists the names of the phases, which determines the names of the files for which the solver looks for phase properties.
Multiphase Turbulence
OpenFOAM can simulate a wide range of physical systems, such as multiphase and compressible flows, heat transfer, etc. which requires a compatible range of turbulence modelling. The issue of compressibility has been managed for many years using two distinct turbulence modelling frameworks, one for constant density flows and another for variable density flows. However, neither framework is appropriate for multiphase systems, in conservative form, for which the phasefraction must be included into all transport and source terms of the turbulence property equations. Code is largely duplicated between the two frameworks, which is inconsistent with the OpenFOAM code development policy to minimise code duplication to promote code maintainablity and sustainability. Extension of the current code architecture to multiphase flows would increase the number of hierarchies from two to four, one for each combination of phasefraction and density representation.
Therefore, as part of the development of multiphase turbulence modelling, a new single, general templated turbulence modelling framework was created that covers all four of the combinations. The framework can be found in $FOAM_SRC/TurbulenceModel (upper case “T”) and includes the incompressible, compressible, phaseincompressible and phasecompressible options including support for laminar, RAS and LES simulations with wallfunctions and other specialised boundary conditions.
The new architecture is demonstrated through a twophase turbulence modelling system built for the general twoPhaseEulerFoam solver, with specialised twophase turbulence models derived from the standard models, but including additional dispersedphase turbulence generation terms. For RAS modelling, this includes:
 the continuousGasKEpsilon model;
 the LaheyKEpsilon model from Lahey R.T, Nucl. Eng. & Design 235:10431060, 2005;
 the mixtureKEpsilon model which takes a basic structure from Behzadi A. et al., Chemical Engineering Science 59:759770, 2004 with a model for bubblegenerated turbulence from Lahey.
The LES modelling includes:
 the continuousGasKEpsilon model;
 the SmagorinskyZhang model from Zhang D. et al., Chemical Engineering Science 61:75937608, 2006;
 the NicenoKEqn model from Niceno B. et al., Chemical Engineering Science 63:39233931, 2008.
This new framework is very powerful and supports all of the turbulence modelling requirements needed so far. It will be enhanced and extended in future OpenFOAM releases to include a wide range of models and submodels, with the expectation to replace the current dual hierarchies of turbulence models, to aid code maintainability and sustainability.
 Source code
 New phase and density templated turbulence model libraries

$FOAM_SRC/TurbulenceModels
Incompressible, twophase RAS turbulence models 
$FOAM_SRC/TurbulenceModels/phaseIncompressible/RAS
Incompressible, twophase LES turbulence models 
$FOAM_SRC/TurbulenceModels/phaseIncompressible/LES
 Examples
 bubble column, RAS 
$FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn
bubble column, LES 
$FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/LES/bubbleColumn
Generalised TwoPhase Solver
The twoPhaseEulerFoam has been consolidated to include all of the functionality from compressibleTwoPhaseEulerFoam, which itself has been deprecated. The solver now includes:
 compressibility effects;
 heattransfer;
 generalised twophase turbulence modelling, including specialised dispersedphase generation terms);
 generalised runtime model selection;
 new models, mainly for flows containing gas bubbles.
 Source code
 twoPhaseEulerFoam solver 
$FOAM_SOLVERS/multiphase/twoPhaseEulerFoam
 Examples
 fluidised bed 
$FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/*/fluidisedBed
bubble column 
$FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/*/bubbleColumn
mixer vessel 
$FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/*/mixerVessel*
Phase Interaction Modelling
The phase interaction modelling, e.g. used in twoPhaseEulerFoam, has been built into a new library of submodels. Many of the models were hardcoded, or had hardcoded coefficients, but the user can now specify the models and coefficients through case input files at runtime. The models are:
 drag models, including noDrag, TomiyamaAnalytic, TomiyamaCorrelated and segregated;
 virtual mass models including noVirtualMass and Lamb;
 heat transfer models, including noHeatTransfer
 lift models, representing the transverse force experienced by a particle of a dispersed phase in a shearing flow, including noLift andTomiyamaLift;
 turbulent dispersion models, to represent the effect of turbulence in the continuous phase in dispersing particles or bubbles, including noTurbulentDispersion, constantTurbulentDispersionCoefficient and Gosman
 wall lubrication models, that represent the tendency of the continuous phase to preferentially occupy nearwall regions, due to the effect of surface tension, including noWallLubrication and Antal.
The effect of phase inversion, where the suspension and droplet phases can switch, has been considered for all the implemented models. For each model type, three models can be specified: one for each phase in a dispersed state, and one for the segregated where neither phase is truly dispersed. For example, for air and water, virtual mass may be specified as follows:
(
(air in water) // air bubbles dispersed in water
{
type Lamb;
}
(water in air) // water droplets dispersed in air
{
type constantCoefficient;
Cvm 0.5;
}
(air and water) // no obvious dispersed phase
{
type none;
}
);
Blending methods are employed to mix the effect of the three models. They are specified per model type, and a default can also be set. The methods available are available:
 noBlending, where the continuous phase is known throughout the flow and can be specified;
 linear, uses a linear basis between specified maximum dispersed phase fractions;
 hyperbolic uses a hyperbolic basis, which is smoother than the linear method, but more expensive to evaluate.
 Source code
 compressibleEulerianInterfacialModels library 
$FOAM_SOLVERS/multiphase/twoPhaseEulerFoam/interfacialModels
 Examples
 fluidised bed 
$FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/*/fluidisedBed
bubble column 
$FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/*/bubbleColumn
mixer vessel 
$FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/*/mixerVessel*