OpenFOAM Multiphase Flows — VOF, Eulerian-Eulerian, and DPM

#+CATEGORY: openfoam # :PROPERTIES: # :ID: uuid-openfoam-multiphase # :END:

OpenFOAM Multiphase Flows :: VOF, Eulerian-Eulerian, and Eulerian-Lagrangian

This note documents OpenFOAM's multiphase modelling capabilities: Volume of Fluid (VOF) for free-surface flows, Eulerian-Eulerian for interpenetrating continua, and Eulerian-Lagrangian (DPM) for discrete particle tracking. These methods map to theory in Conservation Laws (phase conservation) and Incompressible Flow (multiphase incompressible limits).

Multiphase flow is inherently complex: each phase has its own velocity field (in Eulerian-Eulerian), own thermodynamics, and possibly its own turbulence model. The interface-tracking methods (VOF) are simpler but limited to two immiscible phases with a sharp interface. The Eulerian-Eulerian approach is more general but requires closure models for inter-phase momentum, heat, and mass transfer. The Eulerian-Lagrangian approach is the most flexible for dilute particle flows but does not scale well to dense suspensions.

Each approach is a different answer to the same question: how do I represent multiple phases in the conservation equations?

**Table of Contents**

• Volume of Fluid (VOF) Method :: interFoam, MULES, interDyMFoam
• Eulerian-Eulerian (Multi-Fluid) Method :: multiphaseEulerFoam
• Eulerian-Lagrangian (DPM) Method :: DPMFoam, MPPICFoam
• Alpha BCs :: Interfacial boundary conditions
• Surface Tension (CSF) Model :: Continuum Surface Force
• Multiphase Solver Selection :: When to use which approach
• See Also
• References

Volume of Fluid (VOF) Method :: interFoam, MULES, interDyMFoam

The Volume of Fluid method tracks a single volume fraction field =\alpha_i/ for each phase i. For two-phase flow with phases 1 and 2:

= \alpha_1 + \alpha_2 = 1 /

where =\alpha_1 = 1/ indicates pure phase 1, =\alpha_1 = 0/ indicates pure phase 2, and =0 < \alpha_1 < 1/ indicates the interface region (one cell wide). The density and viscosity are computed as volume-fraction-weighted averages:

= \rho = \alpha_1 \rho_1 + (1 - \alpha_1) \rho_2 / = \mu = \alpha_1 \mu_1 + (1 - \alpha_1) \mu_2 /

The phase transport equation is:

= \frac{\partial (\alpha_1 \rho_1)}{\partial t} + \nabla \cdot (\vec{v} \alpha_1 \rho_1) = 0 /

which simplifies (for constant density) to:

= \frac{\partial \alpha_1}{\partial t} + \nabla \cdot (\vec{v} \alpha_1) = 0 /

interFoam solves this equation coupled with the incompressible Navier-Stokes equations. The key challenge is interface compression: without compression, numerical diffusion smears the interface over many cells (turning a sharp interface into a diffusive transition zone). OpenFOAM addresses this through MULES (Multidimensional Universally Limited Encoder).

**MULES (Multidimensional Universally Limited Encoder)**

MULES is a bounded, conservative scheme for advecting the volume fraction field. It enforces = \alpha \in [0, 1]/ at each timestep by applying a multidimensional limiter:

= \alpha_f^{n+1} = \alpha_f^{corrected} + \alpha_{correction} = \alpha_f^{corrected} + \alpha_{max} \min(\alpha_{available}, \alpha_{needed}) /

where the correction is limited by the available volume (the amount that can be removed from the upwind cell) and the needed volume (the amount required at the downwind cell to fill the advected volume). MULES guarantees that the total volume of each phase is conserved globally (the sum of alpha over all cells is constant, up to roundoff).

// fvSchemes for interFoam
divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwind Grad U;
    div(phi,alpha)  MULES interIsoAdvector;  // VOF with MULES + interface compression
    // or: div(phi,alpha)  MULES;           // MULES without compression
}

The interIsoAdvector scheme is the compressive variant — it adds an artificial compression velocity at the interface to keep it sharp:

= \vec{v}_c = \psi |\vec{v}| \hat{n} /

where =\psi/ controls the compression rate and =\hat{n} = \nabla \alpha / |\nabla \alpha|/ is the interface normal (computed from the alpha gradient). The stronger the compression, the sharper the interface — but the higher the risk of oscillations (alpha overshoot > 1 or undershoot < 0). MULES limits these oscillations to maintain boundedness.

**interDyMFoam — VOF with Dynamic Mesh**

For free-surface flows over moving bodies (ship waves, slamming, wave-structure interaction), interDyMFoam combines VOF with OpenFOAM's dynamic mesh framework. The mesh deforms via the ALE (Arbitrary Lagrangian-Eulerian) formulation, and the mesh motion must be solved at each timestep simultaneously with the flow equations:

/* dynamicMeshDict */
dynamicFvMesh       dynamicRefineFvMesh;  // or: solidBodyMeshMotion, svfMotion

solidBodyMeshMotion
{
    solidBodyMotionFunction linearRotation;
    linearCoeffs
    {
        origin (0 0 0);
        axis (0 0 1);
        omega 1.0;    // angular velocity
    }
}

The linearRotation function rotates the mesh at constant angular velocity (for rotating domains). solidBodyMotion functions available in OpenFOAM include: linearTranslation, linearRotation, sixDoF (6 degree-of-freedom rigid body motion), and user-defined motion via a UDF.

The Arbitrary Mesh Interface (AMI) is the OpenFOAM implementation of Fluent's moving mesh interface and Star-CCM+'s sliding mesh — it couples information between rotating and stationary domains. See OpenFOAM Dynamic Meshes for the full dynamic mesh framework.

Eulerian-Eulerian (Multi-Fluid) Method :: multiphaseEulerFoam

The Eulerian-Eulerian approach treats each phase as an interpenetrating continuum — each phase has its own velocity field, its own pressure (or shares a common pressure), and its own conservation equations. The phases are coupled through drag, heat transfer, mass transfer, and turbulence exchange terms.

// Phase 1 (continuum 1, e.g., gas bubbles)
\frac{\partial}{\partial t}(\alpha_1 \rho_1) + \nabla \cdot (\alpha_1 \rho_1 \vec{v}_1) = 0

\frac{\partial}{\partial t}(\alpha_1 \rho_1 \vec{v}_1) + \nabla \cdot (\alpha_1 \rho_1 \vec{v}_1 \vec{v}_1) = -\alpha_1 \nabla p + \nabla \cdot \vec{\tau}_1 + \vec{F}_{1,exchange} + \alpha_1 \rho_1 \vec{g}

// Phase 2 (continuum 2, e.g., liquid)
\frac{\partial}{\partial t}(\alpha_2 \rho_2) + \nabla \cdot (\alpha_2 \rho_2 \vec{v}_2) = 0

\frac{\partial}{\partial t}(\alpha_2 \rho_2 \vec{v}_2) + \nabla \cdot (\alpha_2 \rho_2 \vec{v}_2 \vec{v}_2) = -\alpha_2 \nabla p + \nabla \cdot \vec{\tau}_2 + \vec{F}_{2,exchange} + \alpha_2 \rho_2 \vec{g}

// Volume fraction constraint
\alpha_1 + \alpha_2 + ... + \alpha_N = 1

The interphase exchange terms =\vec{F}_{i,exchange}/ include:

multiphaseEulerFoam is OpenFOAM's default Eulerian-Eulerian solver. It supports up to 9 phases with configurable interphase exchange models. Each phase can have its own turbulence model (RANS, LES, laminar), its own thermophysical properties, and its own equation of state.

**bubbleFoam** — The simplest Eulerian-Eulerian solver. Tracks only the gas volume fraction (= \alpha_g/) and liquid is assumed continuous and stationary:

= \frac{\partial \alpha_g}{\partial t} + \nabla \cdot (\vec{v}_g \alpha_g) = 0 /

with =\vec{v}_g = \vec{v}_l - v_{slip}/, where =v_{slip}/ is the slip velocity from a drift-flux model. This is the open-source answer to Fluent's Eulerian model for bubbly flow.

**twoLiquidMixingFoam** — Two immiscible liquids with no mass transfer. The density and viscosity vary according to the local volume fraction (as in VOF), but the velocity field is shared between phases (single-velocity formulation). This is a simplified Eulerian-Eulerian approach for miscible liquids that share the same velocity.

**fourPhaseEulerFoam / sixPhaseEulerFoam** — Multiphase Eulerian-Eulerian with multiple dispersed phases in a continuous phase. Used for fluidised beds, particle-laden flows, and spray combustion. Each dispersed phase has its own size distribution, its own velocity field, and its own turbulence model.

Eulerian-Lagrangian (DPM) Method :: DPMFoam, MPPICFoam

The Eulerian-Lagrangian approach treats the continuous phase (gas or liquid) as a continuum (Eulerian) and the dispersed phase (particles, droplets, bubbles) as discrete entities (Lagrangian). Each particle is tracked individually:

// Particle motion (Newton's second law)
\frac{d\vec{v}_p}{dt} = \underbrace{\frac{3 \rho C_D Re_p}{4 \rho_p d_p}}_{\text{Drag}} (\vec{v} - \vec{v}_p) + \underbrace{\frac{\vec{g}(\rho_p - \rho)}{\rho_p}}_{\text{Gravity}} + \underbrace{\vec{F}_{other}}_{\text{Lift, Saffman, Virtual mass, etc.}}

// Particle energy
\frac{dT_p}{dt} = \frac{Nu k_g}{\rho_p c_p d_p} (T_g - T_p) + \text{radiation} + \text{phase-change}

// Particle mass
\frac{dm_p}{dt} = -\pi d_p Sh D_v \rho_g (Y_{v,s} - Y_{v,\infty}) \quad \text{(evaporation)}

**DPMFoam** — The standard OpenFOAM DPM solver. Tracks particles through the continuous phase flow field. The continuous phase is solved first (with a standard incompressible or compressible solver), then the particle trajectories are computed as a post-processing step (one-way coupling). For two-way coupling, the particle source terms are fed back into the continuous phase equations:

// Two-way coupling: particle source terms in the continuous phase
\vec{S}_U = \sum_p \frac{\vec{F}_{drag, p}}{V_{cell}} \quad \text{(momentum source)}
\vec{S}_E = \sum_p \frac{\dot{Q}_p}{V_{cell}} \quad \text{(energy source)}
\vec{S}_m = \sum_p \frac{\dot{m}_p}{V_{cell}} \quad \text{(mass source)}

The source terms are computed at each timestep from the particle forces/heat transfer/mass transfer, then fed back into the continuous phase equations. This is two-way coupling (the particles affect the continuous phase, and the continuous phase affects the particles). Four-way coupling adds particle-particle collisions.

**MPPICFoam** (Multiphase Particle-in-Cell) — Handles dense particle suspensions with explicit particle-particle collision modeling. Instead of tracking individual particles, MPPICFoam tracks parcel trajectories where each parcel represents many physical particles. The collision model uses a kinetic-theory-based approach for dense granular flows:

=MPPICFoam/ is the open-source answer to Fluent's Dense Discrete Phase Model (DDPM) with granular kinetics. It is used for fluidised beds, pneumatic conveying, and dense particle transport.

**kineticTheoryFunctionObject** — Provides the granular temperature (=T_g = \frac{1}{3} \overline{|\vec{v}'_p|^2}/) and the related granular pressure and viscosity for dense suspensions. The granular kinetic theory is the OpenFOAM implementation of the Enskog-corrected kinetic theory for dense granular flows (Gidaspow et al., 1992).

Alpha BCs :: Interfacial Boundary Conditions

The alphaContactAngle BCs control the contact angle — the angle at which the free surface meets a solid boundary. This is critical for capturing contact-line dynamics in wetting/drying problems:

// AlphaBC: constant contact angle
alphaContactAngle
{
    type            alphaContactAngle;
    theta0          90;            // equilibrium contact angle (degrees)
    limit           0;             // 0 = constant; 1 = dynamic contact angle
    // for dynamic contact angle:
    // theta0          90;
    // thetaS          60;          // static contact angle
    // thetaA          120;         // advancing contact angle
    // CaMax           0.1;         // maximum capillary number
}

// AlphaBC: inletOutlet (for alpha at inlet)
alphaInlet
{
    type            inletOutlet;
    alpha           alphaInlet;    // phase fraction at inlet (1 = full phase A)
    values          uniform 0;     // zeroGradient for outflow
}

Contact angle directly affects capillary forces in the CSF model (see below). If you simulate droplet impact, capillary flow, or meniscus formation, the contact angle must be specified accurately. For water on clean glass: theta ~ 0° (complete wetting). For water on Teflon: theta ~ 110° (strongly non-wetting).

**dynamicAlphaContactAngle** — The dynamic contact angle model computes theta as a function of the capillary number =Ca = \mu U / \sigma/:

= \theta(Ca) = \theta_0 + A Ca^n /

where =A/ is a coefficient and =n/ is an exponent (typically 0.5–1.0). This accounts for the fact that the effective contact angle changes with contact-line velocity — advancing contact angle (theta_a) > static angle > receding contact angle (theta_r).

Surface Tension (CSF) Model :: Continuum Surface Force

The Continuum Surface Force (CSF) model, originally from Brackbill et al. (1992), models surface tension as a volumetric force applied at the interface:

= \vec{F}_\sigma = \sigma \kappa \nabla \alpha / \rho_{avg} /

where =\sigma/ is the surface tension coefficient (=N/m), =\kappa = -\nabla \cdot \hat{n} = -\nabla \cdot (\nabla \alpha / |\nabla \alpha|)/ is the interface curvature (negative of the divergence of the unit normal), and =\rho_{avg} = \frac{1}{2}(\rho_1 + \rho_2)/ is the average density.

The curvature calculation is critical: a noisy =\kappa/ field generates spurious currents (spurious velocities that circulate along the interface, also known as "parasitic currents"). These can be 10–20% of the free-stream velocity for poorly-resolved interfaces. Recommendations:

The interFoam solver includes surface tension as a built-in capability (controlled by the =sigma/ field). For cases where surface tension is negligible (large-scale ship hydrodynamics), it can be set to zero for computational efficiency.

InterFoam can model the water-air interface with \sigma = 0.072 N/m, or the oil-water interface with \sigma = 0.040 N/m.

Multiphase Solver Selection :: When to Use Which Approach

The choice between VOF, Eulerian-Eulerian, and DPM is a function of the physical regime:

For free-surface ship hydrodynamics (which is the naval application from the training material), interDyMFoam (VOF + dynamic mesh) or LTSInterFoam (VOF with local time stepping in the streamwise direction) are the standard choices. The LTS (Local Time Stepping) option allows larger timesteps in the streamwise direction while the normal-to-free-surface direction is constrained by the Courant limit — this is critical for ship resistance simulations where the domain is much longer than it is deep.

See Also :: Related Notes

. Incompressible Flow — multiphase incompressible limit . Boundary Conditions — alphaContactAngle, alphaInletOutlet . OpenFOAM Dynamic Meshes — interDyMFoam, AMI, dynamicFvMesh . OpenFOAM Solver Selection — multiphase solver decision tree . OpenFOAM Compression Schemes — MULES, interface compression . OpenFOAM Combustion — spray combustion (Eulerian-Lagrangian + DPM) . ANSYS Multiphase — VOF, Eulerian, DPM comparison with ANSYS

References

. Hnes, J.R.F., Brackbill, J.U., & Samuel, B. (1992). "A continuum method for modeling surface tension." Journal of Computational Physics, 100(2), 335-354. (CSF model.) . Weller, H.G. (1999). "A new approach to N-S equation-based multiphase flow simulation." 5th World Conference on Experimental CFD Proceedings. (interFoam derivation.) . OpenFOAM User Guide. ESI OpenCFD. (Multiphase solver documentation.) . Gidaspow, D., Nagendra, R., & Sinclair, C.W. (1992). "Hydrodynamics of fluidization using review data: modeling gas-solid flow regimes." AIChE Journal, 38(5), 685-711. (Granular kinetic theory.) . Munson, B.R., Young, D.F., Okiishi, T.H., & Huebsch, W.W. (2012). Fundamentals of Fluid Mechanics. Chapter on Multiphase flow.