• CFD, Fluid Flow, FEA, Heat/Mass Transfer
  • +91 987-11-19-383
  • amod.cfd@gmail.com

Multi-phase Flow in OpenFOAM

Multi-phase Flows and Discrete Phase Models

Multi-phase flows have wide applications in process, automotive, power generation and metal industries including phenomena like mixing, particle-laden flows, CSTR - Contunuously Stirred Tanks Reactor, Water Gas Shift Reaction (WGSR), fluidized bed, fuel injection in engines, bubble columns, mixer vessels, Lagrangian Particle Tracking (LPT). Some of the general characteristics and categories of multi-phase flow are described below before moving to actual application of OpenFOAM utilities.

 Multiphase flow regimes are typically grouped into five categories: gas-liquid (which are naturally immiscible) flows and (immiscible) liquid-liquid flows, gas-solid flows, liquid-solid flows, three or more phase flows. As can be seen, the immiscibility is a important criteria. In a multi-phase flow, one of the phase is usually continuous and the other phase(s) are dispersed in it. The adjective "Lagrangian" indicates that it relates to the phenomena of tracking a moving points ("fluid particles") - such as tracking a moving vehicle on a road. On the other hand, the adjective "Eulerian" is used to describe correlations between two fixed points in a fixed frame of reference - such as counting the type of vehicles and their speeds while passing through a fixed point on the road.

Gas-liquid flows are further grouped into many categories depending upon the distribution and shape of gas parcels. Three such types are described below.

Bubbly Flow: it represents a flow of discrete gaseous or fluid bubbles in a continuous phase.

Bubbly Flow - Gas Liquid

Slug Flow: This is characterized by flow of large gas bubbles in liquid.

Annular Flow - Gas Liquid

Annular Flow: Here one of the phase if confined to area near the wall forming an annular section.

Annular Flow - Gas Liquid

Some other types of flows are particle-laden flow such as air carrying dust particles, slurry flow where particles are transported in a liquid, hydrotransport which describes densely-distributed solid particles in a continuous liquid such as cement concrete mix. Gas assisted mixing of solid such as fluidized-bed and settling tank where particles tend to sediment near the bottom of the tank forming thick sludge are some other examples of multi-phase flows.

The two dominant method of multi-phase simulations are listed below.
multiphase Methods

This page summarizes the multi-phase cases supplied with OpenFOAM 1606+. Some of the resource from the web are listed below. IPR: the ownership and copyright of these documents lie with the author mentioned in the document.

There are around 21 utilities (solvers) in OpenFOAM dealing with multiphase cases. They provide a ready-to-use sample cases for various flow configurations. A comparison of few utilities are described in the following table. There seems to be a possibility of combining some of the utilities while keeping overall capabilities same.

Comparison of Multiphase Solvers

Excerpts from "LPT for erosion modeling in OpenFOAM: Differences between solidParticle and kinematicParcel, and how to add erosion modeling by Alejandro Lopez" - When dealing with the movement of a group of particles inside a fluid, there are basically two different ways to approach the problem. In the Eulerian-Eulerian models, the particles are treated as a continuous phase and conservation equations are solved for the particulate phase. This method is suitable for large particle concentrations, where two-way coupling between the fluid and the particulate phases as well as particle-particle collisions are important. On the other hand, in the Eulerian-Lagrangian approach, the Eulerian continuum equations are solved for the fluid phase, while Newton’s equations for motion are solved for the particulate phase in order to determine the trajectories of the particles (or groups of particles).

  Thus, if both or all the phases are to be modelled as continuous phases, Eulerian approach is used. When a discrete phase (spatially not continuous) is to be modeled, a Lagrangian frame of reference is used in which spherical particles of pre-defined size distribution is dispersed in the continuous phase (the spherical particles may represent parcels of droplets or bubbles). The fluid phase is modeled as a continuum where time-averaged Navier-Stokes equations are solved to get flow field while the discrete phase is solved by tracking a large number of particles (parcels, bubbles, or droplets) through this calculated flow field. Thus, a Lagrangain approach includes and requires:
  • calculation of discrete phase trajectory for both steady and unsteady flows
  • estimation of effects of turbulence on dispersion of particles due to vortices and eddies present in the continuous phase
  • account for heating or cooling of the discrete phase including break-up of droplets in sprays, vaporization and boiling of such droplets
  • the volume fraction of discrete phase should be < 10% even though mass fraction can be higher due to higher density of solids as compared to liquids. Thus, this method will not be appropriate for fluidized bed where volume fraction of discrete phase is > 50%. Note that twoPhaseEulerFoam should be used in such applications. The Euler-Lagangian solver coalChemistryFoam is not related to fluidized beds.

CSTR, Emulsification

  • Vigorous stirring or agitating two or more immiscible liquids results in the dispersion of one phase in the other in the form of small droplets whose size and shape depend on the equipment, the material properties such as density and surface tension and the operating conditions.
  • The mixing of two immiscible liquids is called emulsification where a binding agent called emulsifier holds the two phase together and dispersed. Icecream, milk, mayonnaise, sauces are some of the examples of emulsion.
  • Prediction of drop size distributions in turbulent liquid-liquid dispersions in static mixers thus becomes one of the most important output of any test set-up or numerical simulation methods.
  • Similarly, any phenomena dealing with droplets require a mathematical model to capture break-up, aggregation and coalesce phenomena of the droplets. One such approach is Population Balance Equations (PBE) - particle size distribution by means of a transport equation, similar to k and ε equations.
  • Since each phase is a continuum, Euler-Euler multiphase models are used to simulation emulsification process.
  • There are three different Euler-Euler multiphase models available in most of the general purpose commercial software and OpenFOAM
    1. VOF Model: The VOF Model is a surface-tracking technique designed for two or more immiscible fluids where the position of the interface between the fluids is of interest. In this model, a single set of momentum equations is shared by the fluids, and the volume fraction of each of the fluids in each computational cell is tracked throughout the domain.
    2. Mixture Model:the phases are treated as interpenetrating continua and solves momentum equation for the mixture and prescribes relative velocities to describe the dispersed phases (but it can also be used without relative velocities for the dispersed phases to model homogeneous multiphase flow). It is a good substitute for the full Eulerian multiphase model (it can perform a full multiphase simulation while solving a smaller number of variables than the Eulerian approach and it is less computationally expensive compared to the Eulerian approach. In this case the momentum, continuity, and energy equations for the mixture, the volume fraction equations for the secondary phases, and algebraic expressions for the relative velocities are solved simultaneously. An algebraic formulation is used for the slip velocity which is based on the assumption of reaching to a local equilibrium between the phases over a short spatial length scale.
    3. Eulerian Model: the phases are treated as interpenetrating continua and solves a set of momentum and continuity equations for each phase. Then the coupling is achieved through the pressure and interphase exchange coefficients.

Variables specific to multi-phase fows
  • alpha: volume fraction
  • p_rgh: Modified pressure p* = p - ρg.x obtained by removing hydrostatic pressure component from static pressure P = pρ.

DPMFoam: Discrete Phase Model - This is another multi-phase solver which includes the effect of the discrete phase particulate volume fraction on the continuous phase. This utility is recommended for dense particle flow simulation. Though there is no strict and standard definition of 'dense' particle-laden flow, a volume fraction of 10% or higher is usually considered 'dense'. The solver uses existing functionality for particle clouds and their collisions, which directly resolves particle-particle interactions. Official description is "Transient solver for the coupled transport of a single kinematic particle cloud including the effect of the volume fraction of particles on the continuous phase.". The term 'cloud' represents the overall presence of all particles, whether they're active or not. Since, this is a transient solver, it is based on PIMPLE algorithm. Turbulence models available in this utility are [a] laminar, [b] k-ε RANS and [c] LES as shown in the header file of the source code.

DPMFoam Turbulence Models

Forces acting on particles
  • Gravity - the well-known force present everywhere!
  • Buoyancy - force due to density difference of the particle and continuous phase.
  • Pressure gradient - the force due to pressure difference along the direction of motion of the particle.
  • Intertia: The acceleration or deceleration of the particle in a fluid requires also an accelerating or decelerating of a fluid surrounding the particle (important for liquid-particle flows).
  • Slip-shear lift force: particles moving in a shear layer experience a transverse lift force due to the nonuniform relative velocity over the particle and the resulting nonuniform pressure distribution
  • Slip-rotation lift force: particles, which are freely rotating in a flow, may also experience a lift force due to their rotation (Magnus force and/or Saffman force)
  • Thermophoretic force: a thermal force moves fine particles in the direction of negative temperature gradients (important for gas-particle flows)
  • The solidParticleCloud class in OpenFOAM is a class that calculates the movement of particles.
Newton Second Law for Particle Tracking
Here subscript 'p' denotes 'particle'. Various empirical correlations of drag laws for viscous forces over solid particles are:
Drag Laws Solid Particles

The object 'particleProperties' defined in dictionary file 'kinematicCloudProperties' is the additional information required for Lagrangian-Eulerian simulations. It defines particle injection rate, particle size and distribution, the forces acting of the particle and empirical model to be used. A sample of kinematicCloudProperties file with comments to describe each of the variable and its significance in the simulation can be find in the this file.

MPPICFoam (MultiPhase Particle-In-Cell method):This is a "Lagrangian solver" having LPT (Lagrangian Particle Tracking) capabilities and can be used for modelling particles in continuum. This method simulates solid phase as parcels and is used to represent collisions without resolving particle-particle interactions. This solver is identical to DPMFoam, but without the collisions between particles.. Quote from a post by Bruno Santos on cfd-online.com: ||The solvers DPMFoam and MPPICFoam are virtually identical, except that the first one uses "basicKinematicCollidingCloud" as the base cloud type and the second one uses "basicKinematicMPPICCloud"||

Phase-coupling mechanisms strongly influences the behavior of the continous and dispersed phase. There are 3 different types of couplings present in particle-laden fluid flows.

  1. One-way coupling: fluid → particles - this is appropriate only with very dilute concentration (< 0.0001% volume fraction) of particles and has no significant effect on turbulence. The flow field is calculated assuming no presence of particle (that is before the particles is injected into the flow) and tracked as they injected into the flow. The particles does not interact with any other particles during its track throughout the flow domain.
  2. Two-way coupling: fluid ↔ particles - this coupling becomes important when the volume fraction of particles are in the range 0.0001% - 0.10% and affects both dissipation and production of TKE. The flow of fluid is necessarily solved along with the movement for Lagrangian particles.
  3. Four-way coupling: fluid ↔ particles + particle collisions: the particle - particle interaction becomes important if voume fraction of particles are > 0.10%.
  4. The particles are assumed to be rigid, spherical and are thus described by diameter, density, coefficient of restitution (rebound properties) and coefficient of friction (tangential drag during collision). Hence, collisions between particles are modeled using keywords identifying properties of spring, friction slider and dash-pot in appropriate dictionaries.

Lagangian >> MPPICFoam >> injectionChannel - The computational geometry and animation of particle tracks in a mixing phenomena of two streams of particle-laden gas is described below.
MPPICFoam in Channel

In addition to the methods outlined earlier to post-process particle tracks, the "Block Extract" features can be directly applied on the kinematicCloud to extract the particles as shown below. This can be done at each time step or the final time step of the simulation.

ParaView Block Extract - DPMFoam

Particle tracking in ParaView

To show particles in ParaView from a kinematicParcleFoam simulation, choose "kinematicCloud-lagrangian" instead of internalMesh in the Mesh Parts. For particle tracking in ParaView, [a] either the utility particleTracks can be used or [b] convert the case to VTK format using foamToVTK utility. Copy the particleTrackProperties dictionary into the constant directory and execute the utility: particleTracks. This will generate necessary files to visualize the particle trajectories in ParaView. The location of the file particleTrackProperties is:

  Both the particleTracks and foamToVTK utility will create a folder named VTK in the case directory which will have futher sub-level folders. The folder VTK/lagrangian/kinematicCloud will contain a file kinematicCloud_TTT.vtk files where TTT is the time step. For particle tracking visualization, one needs to open this file. Then, to display the particles the user must filter the Lagrangian field with the Glyph filter, subsequently the scalar option should be selected in the scale mode window along with the spheres as glyphs where Gyph type options [Sphere, arrow, cone, box...] are available.

Glyph settings to plot particles as sphere is shown below

MPPICFoam - Particle Track VTK Files

MPPICFoam - Particle Track Glyph Setting

MPPICFoam - Particle Track Glyph Setting

The particles coloured by velocity magnitude plotted using VTK files are shown below.
MPPICFoam - Particle Track and Velocity Contour

Lagangian >> simpleReactingParcelFoam >> verticalChannel : a steady state Lagrangian - Eulerian solver for chemical reaction, combustion and particle clouds.

interFoam utility for Multi-phase Flows

Two-phase flows are often broadly categorised by the physical states of the constituent components and by the topology of the interfaces. Thus, a two-phase flow can be classified as gas-solid, gas-liquid, solid-liquid and liquid-liquid in the case of two immiscible liquids . Similarly, a flow can be broadly classfied topologically as separated, dispersed or transitional.


Description as per source code interFoam.c: Solver for 2 incompressible, isothermal immiscible fluids using a VOF (volume of fluid) phase-fraction based interface capturing approach. The momentum and other fluid properties are of the "mixture" and a single momentum equation is solved. Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. For two-fluid approach see twoPhaseEulerFoam. interFoam has limitations in prediction of sharp interfaces, highly self-aerated flows and air-entrainment phenomena.

The following video is an attempt to model air entrainment using multiphase flow. In the left video, you may notice how spherical bubbles are formed due air entrained by water stream entering into the bottle. The interface of the air bubble and water is sharp. This feature is not as prominent in the simulation primarily due to corse mesh. Capturing a sharp interface between gas and liquid will require a very fine mesh. Some improvements planned are:

  1. Make the domain shape and size identical: as of now they are in 2:1 ratio to keep the number of mesh counts low.
  2. Match the boundary conditions such as velocity and its profile at inlet, turbulence parameters at inlet.

multiPhase >> interFoam >> ras >> angledDuct - Flow of water in an empty angled duct: immiscible two phase flow - there are 3 three implementation with LAMINAR, LES and RAS turbulence as tutorial cases. One of the tutorial with RAS implementation is shown below where water enters into the inlet and fills the duct as it moves towards the outlet.

multiPhase >> interFoam >> ras >> weirOverflow - Over-flow a Weir: 2D - immiscible two phase flow
Volume fraction at t = 0:
interFoam Weir
interFoam Weir
interFoam Weir

The animation of change in water volume fraction with time can be downloaded from here.
multiPhase >> interFoam >> ras >> DTCHull - Flow over a ship
Geometry and boundary condition:
interFoam flow over ship

Volume fraction of water at t = 0:
interFoam flow over ship volume fraction of water
multiPhase >> interMixingFoam >> laminar >> damBreak interMixingFoam - "Solver for 3 incompressible fluids, two of which are miscible, using a VOF method to capture the interface." As per the source code, PIMPLE method is used for pressure-velocity coupling.
The settings of tutorial with official version is as follows:

Geometry and Volume fraction of water at t = 0:
interMixingFoam Dam Break alpha

Specification of miscible and immiscible phases:
interMixingFoam Phases

Volume fraction of miscible phases after 2 [s]:
interMixingFoam Volume Fraction


Cavitation refers to the formation of vapor bubbles by local boiling in regions of low pressure within the flow field of a liquid. In some respects, cavitation is similar to boiling, except that the later is generally considered to occur as a result of an increase in temperature rather than a decrease of pressure. Some of the differences between "vapour formation due to cavitation" and that due to boiling are as follows:
  1. It is not possible to raise the temperature of a liquid uniformly where as uniform change in pressure in a liquid can be easily achieved just by controlling the flow passage.
  2. Cavitation always results in self-destruction downstream the flow, causing large pressure fluctuations. Though vapour formation are also self-destructing as the vapour bubbles move away from the cause, "the surface heat source".
  3. Bubble formation involves knowledge of Heat Transfer and boundary layer behaviour of the liquid. understanding the Cavitation phenomena does not need any such information.
  4. Cavitation is almost always detrimental to the surfaces where bubbles collapse. Bubble generation and extinction phenomena in boiling is not so detrimental to the enclosing surfaces.
  5. Typical location of cavitation is low-pressure (suction side) of the centrifugal pumps. A pump is said to get incipient cavitation if the pressure head generated by the pump reduces by 10% of a "non-cavitating" performance.
multiPhase >> cavitatingFoam >> ras >> throttle - Cavitation after sudden expansion from a narrow slit: This case demonstrates use of multi-stage local mesh refinement using combination of "topoSet and refineMesh" utility inside OpenFOAM.

Basic Mesh: Level-0
Base Mesh - Cavitation

Mesh with level-3 refinement near throat

Close-up view of final mesh near throat
cavitatingFoam Refine all


This is an Eulerian-Eulerian solver for two (including compressible) fluid phases where one phase is continuous say water and the other phase is dispersed say gas bubbles or solid particles. It may or may not involve heat transfer. Fluidised bed simulations can be performed using this solver. Note that the flow of solids (the particle bed is bubbling) is being modeled here though it is gas assisted and not the granular flow of solid by its own weight.

 For two phase flows involving fluids, the thermophysicalProperties are specified by adding suffix after the dictionary 'thermophysicalProperties' such as thermophysicalProperties.air, thermophysicalProperties.water. Special variables for this solver and few settings are described below.

  • Theta.particles: Granular temperature of solids phase
  • alpha.air: Volume fraction of air (fluid) phase
  • alpha.particles: Volume fraction of solid phase
  • Note that there is a minimum inlet velocity required to ensure fluidization - that is to move the fluid particles from the bed. Hence, chosen inlet velocity should be higher than or equal to minimum fluidization velocity.
  • Gravity should be active.

Volume fraction and velocity after 1 [s]
Volume fraction and velocity - fluidized bed

Three or more phases, interface capturing capabilities configured to work with LES and not RANS. Phases can be segregated or dispersed. As per source code: "Solver for a system of many compressible fluid phases including heat-transfer". The application to 2D geometry of a mixing vessel is shown below. This is based on Multiple Reference Frame (MRF) method for sliding interfaces - set by dictionary constant/MRFProperties.
multiphaseEulerFoam - mixing vessel MRF

reactingTwoPhaseEulerFoam / reactingMultiphaseEulerFoam:

The simulation of this categories require following 5 files in 'constant' folder -

  1. chemistryProperties - required to include chemical reactions when 'chemistry' is switched on (for example refer to tutorial case reactingTwoPhaseEulerFoam >> RAS >> bubbleColumnEvaporatingReacting),
  2. environmentalProperties - specifies gravity,
  3. combustionProperties - to activate combustion models such as PaSR / XiFoam,
  4. thermophysicalProperties - specifies type of mixtures and its properties, gas phase reaction scheme, presence of inert species,
  5. phaseProperties - this dictionary file is specific to this solver and is used to describe interaction between the phases,
  6. MRFProperties - in case geometry is dealing with rotating domains such as mixer vessel.

A sample 'phaseProperties' dictionary file with comments as per tutorials, online literatures and information available in approppriate *.C and *H file can be acccessed here.

The computational domain for case reactingTwoPhaseEulerFoam >> RAS >>bubbleColumnEvaporatingReacting is shown below. This case simulates Water Gas Shift Reaction (WGSR)taking place in a bubble column reactor. The reacting fluids are mixture of air and CO as reacting gas and liquid water with dissolved air. The reaction products are CO2 and H2 along with water vapor. Like in any other Eulerian solver, the flow of phases are modeled with Eulerian-Eulerian approach along with separate species mass fraction transport equations for each species namely CO, H2O, air, CO2 and H2. This tutuorial deals with a 3D block with a 0.15 x 0.10 [m2] base and a length of 1 [m]. The reactor is filled with water till half of the overall height and have 0.1% dissolved air. Air (with mass fraction of CO = 0.10) enters from the bottom inlet at 0.1 [m/s]. The system is initially assumed to be at 400 [K]. The overall reaction involved is: 0.93CO + 0.24H2O -> 0.69CO2 + H2

Two Phase Reacting Foam


This utility deals with multi-phase flow simulation where phases need to be considered compressible. The suffix 'interFoam' suggest that it is based on Volume-Of-Fluid (VOF) approach. The tutorial case multiphase >> compressibleInterFoam >> laminar >> depthCharge2D deals with a large high-temperature (578 [K]) and pressureized (10 [bar]) air bubble trapped under a water column in a closed cavity. This can be considered a simplified version of high intensity detonation inside water. The computational domain is shown below.
compressibleInterFoam Geometry

As in any VOF simulation with mesh non-optimal mesh refinement, the interface between air and water is not as sharp as observed in real-life examples (refer to the bottle filling animation above). The variation of volume fraction of water (alpha) with time is demonstrated in the following animation. Note how water is raised till upper wall. Also note symmetry during initial stage which gets disturbed once water hits the upper wall and starts falling down. During later stage of simulation, there are zone where the volume fraction of water is in the range 0.5 which is not realistic as water and air are assumed to immiscible and the volume fracion has to be either 0 or 1.

How to convert laminar case to a turbulent case
Change are required in constant/turbulenceProperties file, 0/ directory, system/fvScheme and system/fvSolution files.
constant/turbulenceProperties file
simulationType laminar;simulationType RAS;
  RASModel     kEpsilon;
  turbulence   on;
  printCoeffs   on;
0/ folder
* Add following file and modify boundary names:
divSchemes {
  div(phi,alpha)   Gauss vanLeer;
  div(phirb,alpha) Gauss linear;
  div(rhoPhi,U)   Gauss upwind;
  div(rhoPhi,T)   Gauss upwind;
  div(rhoPhi,K)   Gauss upwind;
  div(phi,p)      Gauss upwind;

  div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
divSchemes {
  div(phi,alpha)   Gauss vanLeer;
  div(phirb,alpha) Gauss linear;
  div(rhoPhi,U)   Gauss upwind;
  div(rhoPhi,T)   Gauss upwind;
  div(rhoPhi,K)   Gauss upwind;
  div(phi,p)      Gauss upwind;
div(rhoPhi,k)   Gauss upwind;
  div(rhoPhi,epsilon)   Gauss upwind;

  div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
"(T|B).*" {
  solver      smoothSolver;
  smoother   symGaussSeidel;
  tolerance   1e-08;
  relTol       0;

"(T|k|epsilon|nuTilda|omega|B).*" {
  solver      smoothSolver;
  smoother   symGaussSeidel;
  tolerance   1e-08;
  relTol       0;


A complete list of tutorial cases related to Lagrangian Particle Tracking (LPT) and multiphase flow is summarized in following tables.

lagrangian coalChemistryFoam simplifiedSiwek
icoUncoupledKinematicParcelDyMFoam mixerVesselAMI2D
icoUncoupledKinematicParcelFoam hopperEmptying
MPPICFoam column
reactingParcelFilmFoam cylinder
reactingParcelFoam counterFlowFlame2DLTS
simpleReactingParcelFoam verticalChannel
sprayFoam aachenBomb
multiphase cavitatingFoam LES throttle
RAS throttle
compressibleInterDyMFoam RAS sloshingTank2D
compressibleInterFoam laminar depthCharge2D
compressibleMultiphaseInterFoam laminar damBreak4phase
driftFluxFoam RAS dahl
interDyMFoam laminar damBreakWithObstacle
multiphase interFoam laminar capillaryRise
LES nozzleFlow2D
RAS angledDuct
interMixingFoam laminar damBreak
interPhaseChangeDyMFoam propeller
interPhaseChangeFoam cavitatingBullet
multiphaseEulerFoam bubbleColumn
multiphaseInterDyMFoam laminar mixerVesselAMI2D
multiphaseInterFoam laminar damBreak4phase
potentialFreeSurfaceDyMFoam oscillatingBox
potentialFreeSurfaceFoam oscillatingBox
reactingMultiphaseEulerFoam laminar bubbleColumn
multiphase reactingTwoPhaseEulerFoam laminar bubbleColumn
LES bubbleColumn
RAS bubbleColumn
twoLiquidMixingFoam lockExchange
twoPhaseEulerFoam laminar bubbleColumn
LES bubbleColumn
RAS bubbleColumn

Contact us
Disclaimers and Policies

The content on CFDyna.com is being constantly refined and improvised with on-the-job experience, testing, and training. Examples might be simplified to improve insight into the physics and basic understanding. Linked pages, articles, references, and examples are constantly reviewed to reduce errors, but we cannot warrant full correctness of all content.