rhoPimpleFoam: Compressible Flow Setup in OpenFOAM

rhoPimpleFoam is OpenFOAM's transient compressible flow solver. It handles density-based coupling between the momentum, energy, and continuity equations using the PIMPLE algorithm. Setting it up correctly requires understanding thermophysical models, the role of pressure and temperature fields, and how Mach number affects solver settings. This guide covers the complete setup from solver selection to fvSchemes.

By CFDpilot · Updated April 2026

1. Choosing the right compressible solver

OpenFOAM provides several compressible solvers. Choosing the wrong one for your Mach number regime is the most common mistake:

• rhoPimpleFoam: Transient, density-based, pressure-implicit. Best for Ma < 0.3 to low supersonic with time-accuracy. Uses PIMPLE loop. Good for flows with moving boundaries or transient effects.
• rhoSimpleFoam: Steady-state equivalent. Use for steady compressible flows at low to moderate Mach numbers. Faster per iteration than rhoPimpleFoam for steady problems.
• rhoCentralFoam: Density-based with central-upwind (Kurganov-Tadmor) scheme. Designed for high-speed flows (Ma > 0.5), shock-capturing. Less suitable for low-Ma subsonic flows due to compressibility noise.

Rule of thumb: if Ma < 0.3, consider whether you truly need compressibility. rhoSimpleFoam or rhoPimpleFoam are appropriate. If Ma > 0.8 with shocks, use rhoCentralFoam.

The PIMPLE algorithm in rhoPimpleFoam combines the outer SIMPLE relaxation (for large time steps) with the inner PISO pressure-velocity correction (for stability). The number of outer PIMPLE iterations is set with nOuterCorrectors. With nOuterCorrectors 1, rhoPimpleFoam behaves as a pure PISO solver and requires Co < 1. With nOuterCorrectors > 1, larger Courant numbers are stable but accuracy is reduced.

2. thermophysicalProperties setup

The constant/thermophysicalProperties file defines the equation of state and transport properties. For an ideal gas with constant specific heats:

// constant/thermophysicalProperties
thermoType
{
    type            hePsiThermo;
    mixture         pureMixture;
    transport       sutherland;
    thermo          janaf;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}

mixture
{
    specie
    {
        nMoles      1;
        molWeight   28.96;  // air, g/mol
    }
    thermodynamics
    {
        Cp          1005;   // J/(kg·K)
        Hf          0;
    }
    transport
    {
        As          1.458e-6;  // Sutherland coefficient
        Ts          110.4;    // Sutherland temperature
    }
}

For simple constant-viscosity air, replace sutherland with const and janaf with hConst, then specify mu and Pr directly.

The Sutherland law models viscosity as a function of temperature: mu(T) = As * T^(3/2) / (T + Ts). For air, As = 1.458e-6 kg/(m·s·K^0.5) and Ts = 110.4 K, giving mu ≈ 1.81e-5 Pa·s at 293 K. Use Sutherland transport for any simulation where the temperature range exceeds ±50 K from the reference condition.

// Simplified constant-property thermophysical model
thermoType
{
    type            hePsiThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}
mixture
{
    specie    { nMoles 1; molWeight 28.96; }
    thermodynamics { Cp 1005; Hf 0; }
    transport { mu 1.81e-5; Pr 0.713; }
}

3. Pressure field: p vs p_rgh

rhoPimpleFoam solves for thermodynamic pressure p in Pascals, not kinematic pressure. Do not confuse this with the p used in incompressible solvers which has units of m²/s².

For flows with significant buoyancy (e.g. natural convection or hot-gas plumes), rhoPimpleFoam can be configured to use p_rgh (modified pressure subtracting the hydrostatic component). Most forced-flow compressible cases use p directly.

// 0/p — thermodynamic pressure in Pa
dimensions      [1 -1 -2 0 0 0 0];
internalField   uniform 101325;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 110000;  // 1.087 atm
    }
    outlet
    {
        type            fixedValue;
        value           uniform 101325;
    }
    walls
    {
        type            zeroGradient;
    }
}

For compressible outlet boundary conditions, avoid zeroGradient for pressure when the flow can be locally supersonic at the outlet. Use waveTransmissive to reduce acoustic reflections in time-accurate cases:

// Outlet BC for subsonic compressible flow (wave-transmissive)
    outlet
    {
        type            waveTransmissive;
        field           p;
        psi             thermo:psi;
        gamma           1.4;
        fieldInf        101325;
        lInf            1.0;
        value           uniform 101325;
    }

4. Temperature field initial conditions

The temperature field T must be initialised to physically realistic values. Setting T = 0 will produce negative absolute temperatures, causing immediate divergence through the equation of state.

// 0/T
dimensions      [0 0 0 1 0 0 0];
internalField   uniform 300;  // Kelvin, not Celsius

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 300;
    }
    outlet
    {
        type            zeroGradient;
    }
    walls
    {
        type            zeroGradient;  // adiabatic wall
    }
}

Common error: specifying temperature in Celsius (e.g. uniform 27 instead of uniform 300). OpenFOAM always works in Kelvin for temperature fields used in thermophysical calculations.

For heated wall boundary conditions, use fixedValue with the wall temperature in Kelvin, or externalWallHeatFluxTemperature when the heat flux (W/m²) is known:

// Wall with specified heat flux (W/m2)
    heatedSurface
    {
        type            externalWallHeatFluxTemperature;
        mode            power;
        q               uniform 5000;
        value           uniform 300;
    }

5. Mach number and numerical scheme selection

The Mach number regime determines the appropriate convection schemes in fvSchemes:

• Ma < 0.3 (subsonic): Use Gauss linearUpwind for momentum. Standard PIMPLE settings.
• 0.3 < Ma < 0.8 (transonic approach): Switch to Gauss upwind or Gauss linearUpwind for all transported quantities. Reduce maxCo to 0.5.
• Ma > 0.8 (transonic/supersonic): rhoPimpleFoam becomes less suitable. Consider rhoCentralFoam with its specialised central-upwind scheme.

// system/fvSchemes for subsonic rhoPimpleFoam
divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwind grad(U);
    div(phi,e)      Gauss linearUpwind grad(e);
    div(phi,K)      Gauss linearUpwind grad(K);
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div((muEff*dev2(T(grad(U))))) Gauss linear;
}

The div(phi,K) term represents convection of kinetic energy K = 0.5*|U|². In compressible solvers, both internal energy and kinetic energy are transported. If you use sensibleInternalEnergy instead of sensibleEnthalpy, replace div(phi,e) with div(phi,Ekp).

6. fvSolution settings for rhoPimpleFoam

The fvSolution file must include solvers for p, U, e (or h), rho, and turbulence fields. The PIMPLE block controls the outer iteration loop:

// system/fvSolution for rhoPimpleFoam
solvers
{
    rho         { solver diagonal; }
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-7;
        relTol          0.05;
    }
    pFinal      { $p; relTol 0; }
    U
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0.1;
    }
    UFinal      { $U; relTol 0; }
    e
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0.1;
    }
    eFinal      { $e; relTol 0; }
}

PIMPLE
{
    nOuterCorrectors        3;
    nCorrectors             2;
    nNonOrthogonalCorrectors 1;
    rhoMin                  0.1;   // minimum density guard
    rhoMax                  10.0;  // maximum density guard
}

The rhoMin and rhoMax bounds clip density during the initial transient to prevent negative values. For air at near-atmospheric conditions, rhoMin 0.1 and rhoMax 10 covers a wide temperature range. The rho solver uses diagonal because density is computed directly from the equation of state (p = rho * R * T) rather than solved iteratively.

7. Common errors and fixes

Negative temperature: Usually caused by wrong initial T (Celsius instead of Kelvin), p initialised at 0, or too large a time step. Fix: check 0/T and 0/p dimensions and values, reduce deltaT.

Pressure units in m²/s² instead of Pa: Incompressible p has dimensions [0 2 -2 0 0 0 0]. If you copy a p file from an incompressible case into rhoPimpleFoam, the solver will crash or give nonsensical density values. Always use [1 -1 -2 0 0 0 0] for compressible p.

Negative density: Caused by a combination of wrong pressure/temperature initialisation and aggressive time stepping. Run the first few steps with maxCo 0.1 and ramp up gradually.

// controlDict — safe starting settings
adjustTimeStep  yes;
maxCo           0.3;  // conservative for compressible startup
maxDeltaT       1e-4;

Floating point exception during first time step: If the solver crashes before writing any output, check 0/p dimensions ([1 -1 -2 0 0 0 0], not [0 2 -2 0 0 0 0]), 0/T in Kelvin, and 0/U velocity magnitude. Also verify that thermophysicalProperties uses molWeight in g/mol — using kg/mol produces density 1000x too high and causes immediate crash.

Very slow convergence or stalled pressure residual: Increase nNonOrthogonalCorrectors to 2 or 3 if mesh non-orthogonality exceeds 50 degrees. Verify the outlet boundary condition allows mass to leave — a fixedValue outlet for U combined with zeroGradient for p can create an over-constrained system.