OpenFOAM Two-Phase VOF Methods

Overview of Volume-of-Frame Method in OpenFOAM

The Volume of Fluid (VOF) method tracks the interface between two immiscible, incompressible fluids by solving a transport equation for the volume fraction $\alpha$ of one of the phases. In OpenFOAM v14, the primary VOF solver is interFoam, which solves the incompressible Navier-Stokes equations with a single velocity field shared across both phases and surface tension modeled via the Continuum Surface Force (CSF) model.

The governing equations are the standard incompressible momentum and continuity equations, with mixture properties computed as weighted averages:

$$ \rho = \alpha\rho_1 + (1-\alpha)\rho_2 $$

$$ \mu = \alpha\mu_1 + (1-\alpha)\mu_2 $$

The continuity equation remains $\nabla\cdot\mathbf{U} = 0$, enforcing a shared velocity field across the interface. The momentum equation uses the mixture density and viscosity. The volume fraction transport equation is:

$$ \frac{\partial \alpha}{\partial t} + \nabla\cdot(\mathbf{U}\alpha) + \nabla\cdot(\mathbf{U}_r \alpha(1-\alpha)) = 0 $$

where the second convective term represents artificial compression to sharpen the interface. Without this term, numerical diffusion would smear the interface over several cell widths, destroying the sharp discontinuity that defines the VOF method.

**Compression Schemes :: Mathematical treatment of interface compression, bounded differencing, and anti-diffusion strategies for sharp interface capture.

**Multiphase Flows :: Broader treatment of multiphase modeling in OpenFOAM, including Euler-Euler, DPM, and interfacial area transport approaches.

interFoam: Core Algorithm

interFoam solves the incompressible, two-phase flow with surface tension using a combination of the PIMPLE algorithm (for pressure-velocity coupling) and the MULES scheme (for volume fraction boundedness). The algorithmic structure per PIMPLE iteration:

1. Solve momentum predictor with mixture density and viscosity. 2. Solve the alpha transport equation with MULES compression. 3. Update mixture properties ($\rho$, $\mu$, $\mathbf{g}$). 4. Solve pressure equation and correct velocity. 5. Optional outer corrector loop (N_outer iterations).

The MULES (Multi-Uniform Limiting Enhancement Scheme) solver is the critical component that ensures $\alpha$ remains bounded:

$$ 0 \leq \alpha \leq 1 $$

at every cell and every time step. MULES operates as a flux limiter: given the convective flux $\phi_\alpha$ at each face, it computes a corrected flux that satisfies the upper bound $0 \leq \alpha \leq 1$ while minimizing numerical diffusion. The compression term:

$$ \nabla\cdot(\mathbf{U}_r \alpha(1-\alpha)) $$

uses a relative velocity $\mathbf{U}_r$ directed normal to the interface, typically estimated from the surface normal $\mathbf{n} = \nabla\alpha / |\nabla\alpha|$. The factor $\alpha(1-\alpha)$ naturally vanishes in regions where a single phase dominates, concentrating compression only near the interface. This is mathematically elegant and numerically robust.

**Dynamic Meshes :: Mesh motion, arbitrary Lagrangian-Eulerian formulation, and mesh quality preservation for moving/free-surface problems.

Surface Tension: CSF Model

The Continuum Surface Force model, due Brackbill, Kothe & Zemach (1992), treats surface tension as a body force localized within a narrow band near the interface:

$$ \mathbf{F}_s = \sigma \kappa \delta_s \mathbf{n} $$

where $\sigma$ is the surface tension coefficient, $\kappa = -\nabla\cdot\mathbf{n}$ is the interface curvature, $\delta_s = |\nabla\alpha|$ is a distribution function that localizes the force to the interface, and $\mathbf{n}$ is the unit normal. In practice, the curvature is computed from the volume fraction gradient, and the force is added to the momentum equation as a source term in interFoam:

$$ \frac{\partial}{\partial t}(\rho\mathbf{U}) + \nabla\cdot(\rho\mathbf{U}\mathbf{U}) = -\nabla p + \nabla\cdot(\mu_{eff}\nabla\mathbf{U}) + \rho\mathbf{g} + \mathbf{F}_s $$

The accuracy of the CSF model depends critically on the curvature calculation. Discrete computation of second derivatives of $\alpha$ is notoriously sensitive to mesh quality and interpolation schemes. OpenFOAM v14 provides multiple options for surface curvature evaluation in the constant/properties file via the curvature entry, including divSurfaceNormal, surfaceNormalConservation, and Gauss linear upwind.

A practical note: for problems dominated by surface tension, the capillary time scale $\Delta t_{cap} = \sqrt{\rho h^3 / (2\pi\sigma)}$ (where $h$ is the cell size) constrains the time step independently of the Courant number. This is the Youngs-Lancaster criterion, and violating it even temporarily leads to interface breakup artifacts.

Dynamic Mesh VOF: interDyMFoam

interDyMFoam extends interFoam with dynamic mesh capabilities. The mesh motion equation is solved alongside the flow equations, with the mesh velocity $\mathbf{U}_S$ appearing in the convection term of the geometric conservation law (GCL):

$$ \frac{\partial V}{\partial t} + \nabla\cdot((\mathbf{U}_S - \mathbf{U})A) = 0 $$

This ensures that a uniform flow field remains uniform even under mesh deformation. The arbitrary Lagrangian-Eulerian (ALE) formulation requires careful treatment at moving interfaces.

**Arbitrary Mesh Interface (AMI):** For rotating domains with free surfaces (a stirred tank with a free surface, for instance), AMI provides non-conformal communication between rotating and stationary mesh regions. The AMI interpolation:

$$ \phi_{target} = \sum_{source} w_i \phi_{source,i} $$

where $w_i$ are flux-weighted interpolation coefficients. AMI preserves conservative transport, making it suitable for volume fraction advection across interface patches. The quality assessment — `AMI check` — reports interpolation weight sums and should ideally yield values between 0.97 and 1.03 at each face; persistent deviation indicates poor mesh alignment or extreme aspect ratio mismatches.

**Dynamic Meshes :: Detailed treatment of dynamicFvMesh, refinement regions, and AMI interface handling.

**ANSYS Multiphase :: ANSYS Fluent's VOF model, geometric reconstruction algorithm, and comparison with OpenFOAM's algebraic VOF approach.

Boundary Conditions for Volume Fraction

Boundary conditions on $\alpha$ dictate how fluid interacts with walls and openings. The most common is alphaContactAngle:

The advancing/receding hysteresis model is physically important for droplet dynamics, capillary rise, and moving contact line problems. Without it, the contact angle jumps between equilibrium values, producing artificial pinning or slip.

**alphaSchumer** is an alternative boundary condition for cases where the contact angle is not well-defined or varies with local flow conditions. It imposes a homogeneous Neumann condition ($\nabla\alpha\cdot\mathbf{n}_S = 0$), effectively allowing the interface to intersect the wall freely.

Local Time Stepping (LTS)

For problems with slow interfacial dynamics (sloshing, wave propagation over long durations), LTS can significantly reduce wall-clock time. Instead of a single global time step, each cell advances at its own maximum-allowed time step based on local Courant number, wave speed, and surface tension criteria. This is activated in constant/dynamicMeshProperties:

```cpp dynamicFvMesh dynamicLtsMesh; lts ltsSchemes; scalar cMax 0.5; scalar maxDT 1; scalar DT 1; ```

The trade-off is that LTS produces non-physical results if the solution is post-averaged over global time intervals — the temporal ordering of events across cells becomes artificially scrambled. LTS is appropriate for statistical steady states (mean wave profiles, average slosh forces) but unsuitable for resolving transient interfacial events.

Comparison with ANSYS Fluent VOF

The key difference: Fluent's Geometric Reconstruction (PLIC) preserves a piecewise-linear interface representation, offering superior geometric conservation at the cost of algorithmic complexity. OpenFOAM's algebraic VOF is more straightforward to implement and parallelize, with MULES providing robust boundedness at the expense of slightly more interface diffusion. For most engineering free-surface problems, the difference is measured in interface thickness of one or two cells, not in integral quantities.

References

• Hirt, C.W. & Nichols, B.D. (1981). Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics, 39(1), 201-225.
• Brackbill, J.U., Kothe, D.B. & Zemach, C. (1992). A continuum method for modeling surface tension. Journal of Computational Physics, 100(2), 335-354.
• Weller, H.G. (1999). Derivation, implementation and validation of computer simulation techniques used in the VOF method for interfacial flows. International Journal for Nuclear Energy, Science and Technology, special issue on VOF.
• OpenFOAM Foundation. (2023). OpenFOAM v2212/v14 User Guide, Chapter on interFoam and interDyMFoam.
• Hysing, S., Kvrano, M., Tjetteland, T., et al. (2010). A benchmark solution for the moving contact line problem. International Journal for Numerical Methods in Fluids, 61(10), 1091-1115.