Transporting a passive scalar — a tracer, a species, a dye — inside only one phase of an interFoam VOF simulation is a deceptively hard problem. A naive scalar equation either leaks the scalar into the other phase or leaves it stranded above a falling free surface. This guide explains why that happens, how the phaseScalarTransport function object solves it, why you post-process alphaS rather than the scalar itself, and the limits you will hit if you need two-way coupling or MULES-sharp gradients.
In a two-phase VOF case (water + air in interFoam), you often want a scalar that lives only in the water — a dissolved species, a tracer for mixing time, a dye for dead-volume detection. If you just add a standard scalar transport equation, two failures show up:
Both come from the same root cause: the flux used to transport the scalar does not track the phase, and the scalar field itself is ambiguous where the phase fraction vanishes.
OpenFOAM Foundation (v7 and later) ships a function object, phaseScalarTransport, built exactly for this: it transports a scalar confined to a single phase. It writes two fields:
s.water — the per-unit-phase scalar it solves for;alphaS.water — the mixture-total field alpha.water * s.water, written for post-processing.The key design point: in interFoam, alphaPhi is not part of the solution, so the function object needs a way to build a phase-consistent flux. If you give it a pressure field, it solves a small pressure-like equation to construct a flux that exactly matches the time derivative of the phase volume. That is what makes the scalar follow the surface down instead of getting stuck.
Add the function object to the functions block in system/controlDict (or an included file):
functions
{
sTransport
{
type phaseScalarTransport;
libs ("libsolverFunctionObjects.so");
field s.water;
p p_rgh; // builds the phase flux in VOF
D 1e-9; // molecular diffusivity (optional)
nCorr 1;
writeControl writeTime;
}
}Two things must be in place or nothing is solved:
0/s.water field with boundary conditions — initialise the scalar where you want it (e.g. via setFields).p p_rgh entry — without it the function object has no flux to convect with and the scalar only diffuses.You do not add a separate div scheme by hand for stock use; the function object looks up its own schemes. If you run a hand-written equation instead, add div(alphaPhi.water,s.water) to fvSchemes.
This is the single most misunderstood point in every thread on the topic. Always read alphaS.water, not s.water.
s.water is the per-unit-phase concentration. Where the phase fraction alpha.water goes to zero, it is effectively undefined, so it smears and appears to leak into the air — this is not a bug, it is just an ill-posed field outside the phase. The mixture-total alphaS = alpha.water * s.water is the quantity the function object actually conserves. In a dam-break test, the mass of alphaS.water is conserved to within about 1%, while the mass of s.water alone can drift by 100-400%.
The same function object works with two-fluid Euler-Euler solvers, but the flux is supplied differently — there alphaPhi is part of the solution, so you point the function object at it instead of giving a pressure field:
sTransport
{
type phaseScalarTransport;
libs ("libsolverFunctionObjects.so");
field s.water;
alphaPhi alphaRhoPhi.water;
rho thermo:rho.water;
}0/s.water, or missing p p_rgh (interFoam) / alphaPhi (Euler). Abort right after start-up and read the top of the log; the function object prints exactly what it cannot find.p p_rgh. Note the ESI scalarTransport phase option had a known bug (issue GL#2053) where the constrained scalar only diffused; the Foundation phaseScalarTransport is the reliable route.bounded01) and use an upwind-biased scheme for the scalar flux before refining.A function object runs after the solver step, so phaseScalarTransport is strictly one-way: it cannot feed alphaS back into the momentum or density equations. If you need the density (or viscosity) to depend on the scalar — a real two-way coupling — you must move the scalar equation into the solver and build a custom solver. Trying to look up alphaS.water from the registry in a stock solver will not even compile, because the field is created at run time by the function object.
The function object uses a standard bounded div scheme, not MULES, and does not expose interface compression. There is an inherent tension: transporting the scalar with the raw compressed flux alphaPhiUn gives a sharp advected gradient but traps the scalar above a falling surface; the conservative alphaS formulation keeps it confined but loses the sharpening. To get both you have to apply MULES limiting to the alphaS flux yourself in a custom solver. For boundedness alone, bounded01 true plus a compressive scheme (interfaceCompression / vanLeer01) on the scalar flux gets part of the way.
Use the phaseScalarTransport function object (OpenFOAM Foundation v7+). Add it to controlDict, give it the scalar field and p p_rgh, and provide a 0/s.water file. It builds a phase-consistent flux so the scalar follows the phase.
alphaS.water. s.water is per-unit-phase and undefined where alpha → 0; alphaS = alpha * s is the conserved mixture quantity (mass conserved to ~1% vs 100-400% drift for s alone).
No phase flux. In interFoam add p p_rgh so the function object can build one. The ESI scalarTransport phase option had a bug (GL#2053) with the same symptom — use the Foundation phaseScalarTransport.
The transport flux does not track the phase volume change. phaseScalarTransport builds a flux from p_rgh that matches the phase volume derivative, so the scalar descends with the surface.
Not with the function object — it is one-way. For two-way coupling, move the scalar equation into a custom solver. alphaS.water does not exist at compile time in a stock solver.
CFDpilot is an AI agent for OpenFOAM — it reads your case from the terminal, checks your phaseScalarTransport setup, the missing p_rgh or 0/s field, and your schemes, and tells you why the scalar is leaking or stuck.
See how CFDpilot works →