17th February 2014

OpenFOAM 2.3.0: Multiphase Modelling

Predictor-Corrector Semi-Implicit MULES

The success of the volume of fluid (VoF) solvers for multiphase flow in OpenFOAM is underpinned by the development of the multi-dimensionsal limiter for explicit solution (MULES) as a very effective method of guaranteeing boundedness of scalar fields, in particular phase/mass-fractions. 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 sub-cycling.

However, in this version of OpenFOAM, a new semi-implicit 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 semi-implicit approach is overall substantially faster than the explicit method with its very strict limit on time-step.

The semi-implicit 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 semi-implicit version of MULES can be switched on through the MULESCorr switch, see example below. For large scale VoF simulations, e.g. ship-keeping simulations, the semi-implicit 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).

solvers
{
    "alpha.water.*"
    {
        nAlphaCorr      1;
        nAlphaSubCycles 1;
        cAlpha          1;
        icAlpha         0.25;

        MULESCorr       yes;      // Switches on semi-implicit MULES
        nLimiterIter    8;        // Number of MULES iterations over the limiter
        alphaApplyPrevCorr  true;fluidised bedmotorbike
        ...

hull_00s
hull_80s

Reference:
el Moctar, O., Shigunov, V., Zorn, T., Duisburg Test Case: Post-Panamax Container Ship for Benchmarking, Journal of Ship Technology Research, 59:50-65, 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.

phases (air water);

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 phase-fraction 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 phase-fraction 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, phase-incompressible and phase-compressible options including support for laminar, RAS and LES simulations with wall-functions and other specialised boundary conditions.

The new architecture is demonstrated through a two-phase turbulence modelling system built for the general twoPhaseEulerFoam solver, with specialised two-phase turbulence models derived from the standard models, but including additional dispersed-phase turbulence generation terms. For RAS modelling, this includes:

the continuousGasKEpsilon model;
the LaheyKEpsilon model from Lahey R.T, Nucl. Eng. & Design 235:1043-1060, 2005;
the mixtureKEpsilon model which takes a basic structure from Behzadi A. et al., Chemical Engineering Science 59:759-770, 2004 with a model for bubble-generated turbulence from Lahey.

The LES modelling includes:

the continuousGasKEpsilon model;
the SmagorinskyZhang model from Zhang D. et al., Chemical Engineering Science 61:7593-7608, 2006;
the NicenoKEqn model from Niceno B. et al., Chemical Engineering Science 63:3923-3931, 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 sub-models, 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, two-phase RAS turbulence models
$FOAM_SRC/TurbulenceModels/phaseIncompressible/RAS
Incompressible, two-phase LES turbulence models
$FOAM_SRC/TurbulenceModels/phaseIncompressible/LES
snappyHexMesh

Examples

bubble column, RAS
$FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn
bubble column, LES
$FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/LES/bubbleColumn

Generalised Two-Phase 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;
heat-transfer;
generalised two-phase turbulence modelling, including specialised dispersed-phase generation terms);
generalised run-time 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*

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 sub-models. Many of the models were hard-coded, or had hard-coded coefficients, but the user can now specify the models and coefficients through case input files at run-time. 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 near-wall 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:

virtualMass
(
    (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*