interFoam is OpenFOAM's solver for two-phase immiscible flows using the Volume of Fluid (VOF) method. It tracks the interface between two fluids (typically water and air) using a phase fraction field alpha.water. Setting it up correctly requires careful attention to the alpha field initialisation, the Courant number limit for the interface, surface tension, and the MULES compression scheme. This guide covers all of these with real dictionary examples and explains every common error.
interFoam solves the Navier-Stokes equations for a mixture of two fluids, where each cell has a phase fraction alpha between 0 and 1. A cell with alpha = 1 is entirely fluid 1 (water), alpha = 0 is entirely fluid 2 (air), and values in between represent the interface region.
The solver uses MULES (Multidimensional Universal Limiter for Explicit Solution) to keep alpha bounded between 0 and 1 and to compress the interface to prevent it from smearing out over many cells. The interface typically spans 2–3 cells regardless of mesh size.
interFoam is appropriate for: dam break simulations, sloshing, wave propagation, coastal engineering, free-surface channel flows, and filling or draining tanks. It is not designed for flows with phase change (boiling, evaporation) — use a dedicated solver like icoReactingMultiphaseInterFoam for those cases.
The phase fraction field lives at 0/alpha.water. It must be initialised to define which regions contain water (1) and which contain air (0). For a dam-break case where the left half of the domain is water:
// 0/alpha.water
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object alpha.water;
}
dimensions [0 0 0 0 0 0 0];
internalField uniform 0; // entire domain starts as air
boundaryField
{
inlet
{
type fixedValue;
value uniform 1; // pure water inlet
}
outlet
{
type inletOutlet;
inletValue uniform 0;
value uniform 0;
}
walls
{
type zeroGradient;
}
atmosphere
{
type inletOutlet;
inletValue uniform 0;
value uniform 0;
}
}
For complex initial interface geometries, use setFields with a setFieldsDict in the system/ directory to set alpha = 1 in a geometric region (box, sphere, cylinder, etc.).
setFields is the standard tool for initialising the alpha field when the initial interface is not trivial. It applies field values to geometric regions defined in system/setFieldsDict.
// system/setFieldsDict — initialise water in a rectangular region
defaultFieldValues
(
volScalarFieldValue alpha.water 0 // default: all air
);
regions
(
boxToCell
{
box (0 0 0) (0.5 0.3 1); // water fills left half up to y=0.3
fieldValues
(
volScalarFieldValue alpha.water 1
);
}
);
Run setFields after blockMesh and before interFoam. For STL-defined initial interfaces, replace boxToCell with surfaceToCell and specify the STL file that defines the interface boundary.
In OpenFOAM v2106+ the fluid properties are defined in constant/physicalProperties. In older versions, use constant/transportProperties:
// constant/physicalProperties (v2106+)
phases (water air);
water
{
transportModel Newtonian;
nu 1e-06; // kinematic viscosity [m²/s]
rho 1000; // density [kg/m³]
}
air
{
transportModel Newtonian;
nu 1.48e-05;
rho 1;
}
sigma 0.07; // surface tension coefficient [N/m], water-air at 20°C
The density ratio for water-air (1000:1) is numerically challenging. If your simulation diverges early, try reducing the density ratio temporarily (e.g. rho_water = 100) to diagnose whether it is a stability or setup issue.
The gravity constant is defined in constant/g. This is separate from the physical properties and is required for any simulation involving buoyancy or free-surface flows. Getting the direction wrong is a very common mistake.
// constant/g
FoamFile { object g; }
dimensions [0 1 -2 0 0 0 0];
value (0 -9.81 0); // gravity in -y direction (vertical domain)
If your domain has z as the vertical axis, use (0 0 -9.81). If you forget the g file, interFoam will run without buoyancy, and the free surface will not behave correctly.
interFoam requires two Courant number limits in system/controlDict: the standard maxCo (for the velocity field) and maxAlphaCo (for the alpha advection, i.e. the interface movement). The interface Courant number should be kept below 0.5 to maintain a sharp interface.
// system/controlDict
application interFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 5;
deltaT 0.001;
writeControl adjustableRunTime;
writeInterval 0.1;
adjustTimeStep yes;
maxCo 1;
maxAlphaCo 0.5; // critical for interface sharpness
maxDeltaT 0.05;
interFoam uses the PIMPLE loop internally (like pimpleFoam). The fvSolution file needs entries for the pressure solver, the velocity solver, and the PIMPLE settings. The alpha equation is handled internally and does not appear in fvSolution:
// system/fvSolution
solvers
{
"alpha.water.*"
{
nAlphaCorr 2;
nAlphaSubCycles 1;
cAlpha 1; // interface compression factor (1=standard)
icAlpha 0;
MULESCorr yes;
nLimiterIter 3;
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-8;
relTol 0;
}
p_rgh
{
solver PCG;
preconditioner DIC;
tolerance 1e-7;
relTol 0.05;
}
p_rghFinal
{
$p_rgh;
relTol 0;
}
U
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-6;
relTol 0;
}
}
PIMPLE
{
momentumPredictor no;
nOuterCorrectors 1;
nCorrectors 3;
nNonOrthogonalCorrectors 0;
}
Note that interFoam uses p_rgh (modified pressure minus hydrostatic head), not p directly. Your boundary conditions for pressure should reference p_rgh, not p.
The boundary conditions for p_rgh differ from those for p in single-phase solvers. The correct pairings for typical interFoam cases are:
// 0/p_rgh — typical interFoam setup
boundaryField
{
inlet
{
type fixedFluxPressure; // consistent with velocity flux at inlet
value uniform 0;
}
outlet
{
type fixedValue;
value uniform 0; // gauge pressure = 0 at outlet
}
walls
{
type fixedFluxPressure;
value uniform 0;
}
atmosphere
{
type totalPressure;
p0 uniform 0; // open atmosphere: total pressure = 0 gauge
value uniform 0;
}
}
The atmosphere patch (top boundary in sloshing or wave tank cases) uses totalPressure with p0 = 0 to represent the open air boundary. The corresponding U condition at the atmosphere patch is pressureInletOutletVelocity.
The alpha advection scheme is critical for interface sharpness. The interfaceCompression scheme is the standard choice. Avoid using central differencing or linear schemes for the alpha equation — they cause severe smearing.
// system/fvSchemes — relevant entries for interFoam
divSchemes
{
default none;
div(rhoPhi,U) Gauss linearUpwind grad(U);
div(phi,alpha) Gauss interfaceCompression vanLeer 1;
div(phi,k) Gauss limitedLinear 1;
div(phi,omega) Gauss limitedLinear 1;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}
The compression factor (second argument after vanLeer) controls how aggressively the interface is sharpened. A value of 1 is the standard. Values above 1 (up to 2) give a sharper interface but can cause instability. Values below 1 produce a more diffuse interface.
MULES is supposed to bound alpha between 0 and 1, but it can fail if the time step is too large. The interface Courant number being above 1 is the most common cause. Fix: ensure maxAlphaCo 0.5 is set in controlDict and that adjustTimeStep yes is active.
At open boundaries (outlets, atmosphere patches), alpha must use inletOutlet — not zeroGradient. If backflow occurs at an outlet with zeroGradient on alpha, unphysical values can be advected back in. Use inletOutlet with inletValue uniform 0 for air-side outlets.
If the interface develops checkerboard oscillations or ripples that grow over time, increase nAlphaSubCycles from 1 to 2, and increase nLimiterIter to 5. Also check that cAlpha is set to 1 (the default). Setting cAlpha 0 disables compression and leads to a very diffuse interface.
interFoam may print a warning about the pressure reference cell if the domain has no outlet with a fixed pressure. Add a pRefCell 0; pRefValue 0; entry in the PIMPLE block of fvSolution to fix this.
If constant/g does not exist, interFoam runs without gravity. The simulation appears to start normally but free-surface shapes are wrong. Always create the g file with the correct orientation.
The interface width in VOF naturally grows over time due to numerical diffusion. To limit this, increase the compression factor in the div(phi,alpha) scheme from 1 to 1.5, and increase nAlphaSubCycles to 2. Also check your mesh resolution at the interface — at least 3–5 cells across the interface region is needed for a reasonable result.
Upload your multiphase case and CFDpilot checks your alpha field, surface tension settings, and VOF scheme configuration.
Diagnose my case →p_rgh is the modified pressure equal to p minus the hydrostatic head (rho * g * h). interFoam solves for p_rgh rather than p directly because it eliminates the hydrostatic pressure gradient from the momentum equation. All boundary conditions in 0/ should reference p_rgh, not p.
Set maxAlphaCo to 0.5 in system/controlDict with adjustTimeStep yes. Values above 1 cause MULES to fail at bounding alpha between 0 and 1, leading to alpha blow-up. For cases with very sharp interfaces or high-density ratios (water-air), reduce maxAlphaCo to 0.2 for extra stability.
Use setFields with a setFieldsDict in the system/ directory. Define defaultFieldValues with alpha.water = 0 (all air), then add region overrides using boxToCell, sphereToCell, or surfaceToCell to set alpha.water = 1 in the water region. Run setFields before running interFoam.
Alpha values outside [0,1] indicate MULES has failed to bound the phase fraction. The most common cause is maxAlphaCo being too high. Fix: set maxAlphaCo 0.5 and adjustTimeStep yes in controlDict. Also verify that nAlphaCorr is at least 2 and MULESCorr is set to yes in the alpha.water solver block of fvSolution.
The surface tension coefficient sigma for water-air at 20°C is 0.07 N/m. Set this in constant/physicalProperties (v2106+) or constant/transportProperties (older versions). At 60°C, sigma drops to approximately 0.066 N/m. Surface tension only matters significantly when the Weber number is low (We < 10).
Yes. interFoam supports RANS and LES turbulence models through the same momentumTransport (or turbulenceProperties) file as single-phase solvers. For free-surface flows with turbulence, kOmegaSST is recommended. Configure the turbulence fields (k, omega/epsilon, nut) in the 0/ directory as you would for any other solver.