VOF Method (Volume of Fluid) for Multiphase Flow CFD

Exactly. VOF (Volume of Fluid) shines for problems with free surfaces or fluid-fluid interfaces: ocean wave simulation, liquid sloshing in LNG tanks, fuel sloshing in aircraft wing tanks, droplet impact on hydrophobic surfaces, dam break, fuel injection sprays, and ship hydrodynamics. The method represents "how much of each fluid is in each cell" through a volume fraction $\alpha$, implicitly tracking the interface. Hirt and Nichols proposed VOF in 1981 at Los Alamos, and it remains one of the most widely used multiphase CFD approaches today. The key strength is exact mass (volume) conservation — a property many competing methods struggle to maintain.

Yes — cells that straddle the interface have $0 < \alpha < 1$. The interface is smeared over 2–3 cells wide in the numerical representation. How sharply this smearing is controlled depends on the interface reconstruction scheme. The Geo-Reconstruct (PLIC) method reconstructs a planar interface segment within each interface cell, keeping the interface very sharp. Simpler methods like HRIC are faster but can allow more smearing. In ship hydrodynamics or tank sloshing where interface sharpness affects wave pattern accuracy, PLIC is the professional standard.

The core equation is the advection equation for volume fraction $\alpha$:

$\alpha = 1$ for fluid 1 (e.g., water), $\alpha = 0$ for fluid 2 (e.g., air). One set of Navier-Stokes equations is solved for the shared velocity field, with mixture density and viscosity interpolated:

OpenFOAM's interFoam solver adds an artificial compression term to counteract numerical diffusion:

The compression velocity $\mathbf{u}_r$ is directed normal to the interface with magnitude proportional to the local velocity. This term is only active in the interface region ($0 < \alpha < 1$) and artificially compresses it toward the sharp interface limit without violating mass conservation.

Yes — the Continuum Surface Force (CSF) model by Brackbill et al. (1992) converts the surface tension discontinuity at the interface into a continuous volumetric body force:

$\sigma$ [N/m] is the surface tension coefficient (water-air: 0.072 N/m at 20°C), $\kappa = -\nabla \cdot \hat{n}$ is the interface curvature, $\hat{n} = \nabla\alpha / |\nabla\alpha|$ is the interface unit normal. This body force appears in the momentum equation. A known weakness of CSF is "spurious (parasitic) currents" — nonphysical velocity circulations around a stationary droplet due to curvature discretization error. These artificial velocities can be 1–100× the physical flow velocity for very fine mesh, causing solver instability. Using the Height Function method for curvature calculation (standard in modern codes) reduces spurious currents by 1–2 orders of magnitude compared to gradient-based curvature.

Each method has different strengths:

VOF is the dominant industrial choice because mass conservation is non-negotiable for long-time simulations. Level Set gives better curvature accuracy (useful for droplet shape studies) but requires reinitialization and mass correction schemes. Phase Field (Cahn-Hilliard) naturally handles coalescence and breakup physics but is expensive. CLSVOF combines the best of both VOF and Level Set and is used in research codes for high-accuracy droplet dynamics.

The interface reconstruction scheme is the key control. Here's a comparison:

The interface Courant number must be controlled — if the interface moves more than one cell per time step, $\alpha$ cannot be transported accurately and you get interface smearing or numerical instability. The target is:

For explicit schemes this constraint is strict, making high-velocity interface problems (e.g., high-pressure fuel injection sprays at 500 m/s) extremely expensive. For a 10 μm mesh and 500 m/s velocity, $\Delta t \le 10$ ns — billions of time steps for a 10 ms injection event. Fluent's implicit VOF scheme relaxes the Courant limit to Co ≤ 2–4, but at some cost to interface sharpness. Adaptive time step control (ATSC) in OpenFOAM automatically adjusts the time step to maintain Co ≤ 0.5 using the maxCo and maxAlphaCo settings. There is also a surface tension Courant number constraint for surface-tension-dominated flows:

This constraint can be more restrictive than the convective Courant number for fine meshes with high surface tension (small droplets in water).

Tank sloshing is a classic VOF application. Here's a complete checklist:

1. Initial condition: Patch the $\alpha$ field to match the initial fill level. Use Fluent's Patch command or OpenFOAM's setFields utility with a box region defining the liquid zone.
2. Gravity vector: Don't forget to set the gravity direction accurately ($g = 9.81$ m/s² in the correct direction). For oscillating or tilting loads (ship motion, road vehicle on a curve), implement the time-varying acceleration via UDF or expressions. For LNG carriers, the excitation is the wave-induced ship motion — typically 6-DOF rigid body motion.
3. Mesh refinement at the free surface: Ensure 3–5 cells across the expected wave height at the quiescent free surface level. AMR (adaptive mesh refinement) that tracks the interface automatically is highly efficient for long transient runs.
4. Time step: Verify $Co_\alpha \le 0.5$ throughout. Use adaptive time stepping. Typical $\Delta t$ for a 10 m tank at 3 m/s peak velocity: 0.001–0.01 s.
5. Volume conservation check: Monitor total liquid volume periodically — it should remain within 1% of the initial value throughout. Larger errors indicate time step or scheme issues.
6. Pressure loads: Record wall pressure time histories at multiple points — especially near the corners and liquid impact zones. Peak sloshing pressure (impact pressure) for structural design is the key output. Extreme value statistics require long simulations (1000+ cycles).

Droplet impact is one of the most demanding VOF problems because it involves topology changes (droplet merges with the pool), high curvature at the impact point, and potential air entrainment. Key settings: (1) Use PLIC (Geo-Reconstruct) for the sharpest interface capture. (2) Set the contact angle at the solid surface correctly — advancing angle for wetted surfaces (glass: ~20°), larger for hydrophobic surfaces (lotus leaf: ~150°). (3) Use an extremely fine mesh around the impact zone — at least 10 cells across the droplet diameter. For a 100 μm ink droplet, that means 10 μm cells over a 200–300 μm region. (4) The surface tension Courant constraint will likely dominate — check it explicitly before setting the time step.

Here's a comparison of the major platforms:

Several exciting directions are active:

• CLSVOF (VOF + Level Set): Hybrid approach combining VOF's mass conservation with Level Set's smooth interface representation for better curvature accuracy. Used in research codes for high-fidelity droplet and bubble simulations.
• AMR with VOF: Automatic local mesh refinement near the interface reduces computational cost by 10–100× compared to uniform fine mesh. STAR-CCM+ and OpenFOAM's dynamicRefineFvMesh enable this.
• ML-assisted interface reconstruction: Neural networks trained on high-resolution VOF data reconstruct sharp interfaces from coarse-mesh $\alpha$ fields — effectively a super-resolution approach for VOF.
• VOF for additive manufacturing: Modeling the melt pool in laser powder bed fusion (PBF) and DED processes requires VOF with coupled heat transfer, phase change, and Marangoni (thermocapillary) convection. An extremely active application area since 2020.
• Lattice Boltzmann VOF: Combining LBM's mesoscale physics with VOF's interface tracking for pore-scale multiphase flow in porous media — important for battery electrodes and fuel cell diffusion layers.

The main failure modes in order of frequency:

• Divergence from large Courant number: Check $Co_\alpha$ in every time step. If using fixed time step and interface speed is higher than expected, the solver will diverge. Switch to adaptive time stepping.
• Interface smearing: Interface diffuses from 2–3 cells to 5–10 cells width over time. Causes: high Courant number, upwind-dominated scheme, or insufficient mesh resolution. Switch from HRIC to PLIC, reduce time step, or refine the mesh.
• Spurious currents around droplets/bubbles: Artificial velocities appear near the interface even in the absence of external flow. Reduce mesh or switch to Height Function curvature. Limit surface tension coefficient as a diagnostic — if spurious currents disappear with $\sigma = 0$, it's a CSF discretization issue.
• Mass non-conservation: Monitor total liquid volume every 100 time steps. If volume changes by more than 1%, check for improper outlet boundary conditions or a time step too large for the Courant constraint.
• $\alpha$ values outside [0,1]: MULES (Multidimensional Universal Limiter with Explicit Solution) in OpenFOAM bounds $\alpha$ to [0,1]. In Fluent, the bounded scheme is automatic. If values outside this range appear, check for overly coarse mesh or excessive time step.