rhoCentralFoam: Negative Temperature and Floating Point Exceptions

rhoCentralFoam is OpenFOAM's density-based, shock-capturing solver — and one of the most temperamental. A run that looks fine for a few iterations suddenly dies with a "negative temperature" message or a floating point exception. This guide explains why that happens, how the solver's numerics actually work, and the concrete settings — Courant number, reconstruction schemes, boundary conditions, the limitTemperature safety net, and source-term pitfalls including porous zones — that get it stable.

By CFDpilot · Updated June 2026

1. Why rhoCentralFoam is fragile

Unlike the SIMPLE/PIMPLE pressure-based solvers, rhoCentralFoam is density-based and explicit. It uses a central-upwind (Kurganov-Tadmor / Kurganov-Noelle-Petrova) scheme to capture shocks without a Riemann solver. Two consequences follow, and they explain almost every crash:

• The time step is limited by the acoustic Courant number — based on |U| + c (flow speed plus speed of sound), not just |U|. In supersonic flow that is a tight limit.
• The stabilisation lives in the reconstruction limiters, not in divSchemes. The convective fluxes are built from limited interpolations of the face values (the reconstruct(...) entries) and the fluxScheme. This is where you control how dissipative — and how stable — the solver is.

The solver recovers temperature from the total energy field. If the energy in a cell is pushed below the physical floor — by an overshoot at a shock, a too-large time step, a reflecting boundary, or an inconsistent source term — temperature goes negative and the run dies on a floating point exception.

2. Courant number and time control

The single most common cause of a sudden crash is an excessive Courant number. For shocked or supersonic flow:

// system/controlDict
adjustTimeStep  yes;
maxCo           0.2;     // 0.1 or lower if it still crashes
maxDeltaT       1e-6;

If the case dies in the first handful of iterations, the start-up transient is the problem: the initial field is far from physical and the fluxes are violent. Start from a lower maxCo (0.05–0.1), let it settle, and ramp up.

3. Schemes that stabilise it

For rhoCentralFoam the convective robustness is set by the fluxScheme and the reconstruct(...) limiters — not by divSchemes.

// system/fvSchemes
fluxScheme      Kurganov;     // or Tadmor

ddtSchemes      { default Euler; }
gradSchemes     { default Gauss linear; }

divSchemes
{
    default         none;
    div(tauMC)      Gauss linear;
}

interpolationSchemes
{
    default             linear;
    reconstruct(rho)    Minmod;     // most dissipative / most robust
    reconstruct(U)      MinmodV;
    reconstruct(T)      Minmod;
}

Key points:

• Use Minmod first. Minmod is the most dissipative limiter and the safest for strong shocks. vanLeer is sharper but overshoots more easily; SuperBee is the least dissipative. Get a stable run on Minmod, then relax toward vanLeer if you need crisper shocks.
• Don't carry over pressure-based divSchemes. Entries like div(phi,e), div(phi,K) and div(phi,U) belong to segregated solvers (rhoPimpleFoam); stock rhoCentralFoam only needs div(tauMC). If you pasted them in from another case they are at best inert and at worst (on a modified solver) an unbounded linearUpwind on energy that overshoots at shocks.

4. Boundary conditions for high-speed flow

Bad boundary conditions reflect waves back into the domain, and the reflection is often where the temperature spikes. Two to watch:

• waveTransmissive outlet needs a correct lInf (relaxation length) and fieldInf (far-field value), plus gamma. A poorly tuned one reflects pressure waves and drives T excursions upstream.
• Supersonic inlets: if the flow is fully supersonic at the inlet, every variable must be fixed; totalPressure/totalTemperature are for subsonic/transonic inlets (e.g. the stagnation chamber of a nozzle). Make sure the inlet condition matches the local Mach number.

5. The limitTemperature safety net

While you track down the cause, you can stop the run from dying on the exception by clipping temperature with an fvOption:

// system/fvOptions
limitT
{
    type            limitTemperature;
    active          true;
    selectionMode   all;       // or cellZone <name>
    min             150;       // Kelvin (must be > 0)
    max             1500;
}

This is a diagnostic crutch, not a fix. If the temperature is being pinned at min every iteration, the physics or numerics are still wrong — but the run survives long enough for you to see where it is going wrong (usually a shock, a boundary, or a source-term cell zone).

6. Source terms and porous zones

Adding an fvOptions source to rhoCentralFoam is where many negative-temperature crashes start. The most common case is a porous zone via explicitPorositySource / DarcyForchheimer:

• Porosity is a momentum sink only. The DarcyForchheimer term removes momentum, but the standard implementation does not return the dissipated kinetic energy to the energy equation as heat — the porosity model is not energy-consistent for compressible flow. At low speed this is negligible; in a supersonic nozzle the kinetic energy removed is large, so the energy bookkeeping no longer balances and the temperature collapses at the porous interface.
• The porosity sources were built for pressure-based solvers (rhoSimpleFoam, rhoPimpleFoam). On rhoCentralFoam's explicit density-based update, a stiff momentum sink is fragile on its own, before the energy issue.
• A huge isotropic resistance makes it worse. A Darcy coefficient like d = (1e10 1e10 1e10) is effectively a solid block applied in every direction including streamwise — slamming supersonic flow into a near-wall.

What works: for compressible flow with porous media, move to rhoPimpleFoam, where the fvOptions porosity is properly integrated and the energy coupling is handled. Calibrate d and f against a known pressure drop, make the resistance anisotropic (little or no streamwise resistance for a flow-aligned perforated plate), and ramp it up over time rather than switching a huge value on at t = 0. And always confirm the case is stable with the source turned off (active false) first — that isolates whether the crash is the source or the base nozzle setup.

FAQ

Why does rhoCentralFoam give a negative temperature error?

The solver recovers temperature from total energy; when fluxes overshoot at a shock, the Courant number is too high, a boundary reflects, or a source term removes energy inconsistently, the internal energy in a cell goes below zero and temperature with it. Make the numerics more dissipative (Minmod), lower the time step, and check the BCs and source terms.

How do I fix a floating point exception in rhoCentralFoam?

Lower maxCo to 0.2 or below, switch reconstruction limiters to Minmod, confirm fluxScheme is Kurganov or Tadmor, check supersonic boundary conditions, and add a limitTemperature fvOption as a safety net while you find the real cause.

What Courant number should I use for rhoCentralFoam?

Because it is explicit and density-based, the limit is the acoustic Courant number (|U| + c). Keep maxCo at 0.2 or below for shocked/supersonic flow, lower with source terms or strong gradients, and use adjustTimeStep with a sensible maxDeltaT.

Which reconstruction scheme is most stable in rhoCentralFoam?

Minmod — it is the most dissipative limiter and the safest for strong shocks. vanLeer is sharper but overshoots more; SuperBee is the least dissipative. Start with Minmod, then relax if you need sharper resolution.

Can I use a porous zone with rhoCentralFoam?

Not robustly. Porosity (DarcyForchheimer) is a momentum sink whose dissipated energy is not returned as heat, so it drives temperature negative in high-speed flow, and the fvOptions porosity sources were built for pressure-based solvers. Use rhoPimpleFoam instead, calibrate the resistance, and ramp it up.

What is limitTemperature in OpenFOAM?

An fvOption that clips the temperature field between a minimum and maximum over a region. Useful to keep a run alive past a negative-temperature crash, but a diagnostic crutch, not a fix — if temperature sits at the limit constantly, the underlying problem remains.