Troubleshooting Guide

interFoam VOF Simulation: Setup and Common Errors

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.

By CFDpilot · Updated April 2026

1. What interFoam does

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.

2. The alpha.water field

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.).

3. Using setFields to initialise the interface

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.

4. Physical properties — surface tension and densities

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.

5. Gravity vector

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.

6. Courant number for the interface — maxAlphaCo

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;

7. fvSolution for interFoam

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.

8. p_rgh boundary conditions

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.

9. fvSchemes for interFoam

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.

10. Common errors and fixes

Alpha field blowing up (alpha > 1 or alpha < 0)

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.

Wrong boundary conditions for alpha

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.

MULES instability — oscillations at interface

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.

p_rgh reference cell error

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.

Missing g file

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.

Interface too diffuse after many time steps

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.

Diagnose your interFoam setup instantly — free

Upload your multiphase case and CFDpilot checks your alpha field, surface tension settings, and VOF scheme configuration.

Diagnose my case →
Official documentation

Frequently Asked Questions

What is the difference between p and p_rgh in interFoam?

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.

What value of maxAlphaCo should I use in interFoam?

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.

How do I initialise the alpha.water field for a complex geometry?

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.

Why does my interFoam simulation show alpha values outside 0 and 1?

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.

What surface tension value should I use for water-air in interFoam?

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).

Can I use interFoam with turbulence models in OpenFOAM?

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.