OpenFOAM: boundary conditions and physical models in a complete case

For a constant-density Newtonian fluid, in kinematic form (divided by ρ):

∇·U = 0
∂U/∂t + ∇·(U U) = −∇p + ∇·[ ν_eff (∇U + (∇U)ᵀ) ]

p = p_phys / ρ      [m²/s²]
ν_eff = ν (laminare)  oppure  ν + ν_t (RAS/LES)

Two practical consequences. First: the solver does not know ρ, so neither pressure nor forces come out in SI units until multiplied by a declared density. Second: in incompressible flow pressure is not a thermodynamic variable but the multiplier enforcing ∇·U = 0; it is therefore defined up to a constant, which is why at least one reference value is needed (here p = 0 at the outlet). The div((nuEff*dev2(T(grad(U))))) term in fvSchemes corresponds to the transposed part of the viscous stress.

For incompressible Navier–Stokes, on each portion of the boundary one generally prescribes either velocity (inlets, walls) or a “traction” condition, which OpenFOAM in practice approximates by fixing pressure and a zero velocity gradient (outlets). The proper combination for a duct is therefore: velocity or flow rate at the inlet with zero pressure gradient; fixed pressure at the outlet with zero velocity gradient; no-slip on walls. Prescribing both velocity and pressure on the same patch over-constrains the problem; prescribing pressure nowhere leaves it undetermined.

In ducts, transition to turbulence occurs under ordinary conditions around Re ≈ 2000–2300 based on the hydraulic diameter; very disturbed inlets can bring it forward, very clean ones can delay it considerably. With Re ≈ 133 the case is laminar by a wide margin. The second number to look at is the ratio of channel length to entrance length, L/L_e ≈ 0.20/0.044 ≈ 4.5: the flow develops and stays developed for about three quarters of the channel, a necessary condition for comparing it with the analytic solution.

The quickest check is to read the banner of any log and compare file names with those of the installed tutorial. An extra dictionary does not always produce an error: it is simply ignored, which is the most insidious case because you believe you set something the solver never reads.

ls "$FOAM_TUTORIALS/incompressibleFluid"
ls "$FOAM_TUTORIALS/incompressibleFluid/pitzDailySteady/constant" "$FOAM_TUTORIALS/incompressibleFluid/pitzDailySteady/system"
foamDictionary "$FOAM_TUTORIALS/incompressibleFluid/pitzDailySteady/system/fvSolution" -keywords

Conversion is not done from memory: copy the equivalent version-14 tutorial and transfer mesh, fields and values of the original case into it, checking each dictionary with foamDictionary -keywords. Read the log of the first run in full, looking for warnings about ignored entries.

For fully developed laminar flow in a rectangular duct there is a series solution (Shah and London). With a 10 mm × 5 mm section (aspect ratio 0.5) the product of Darcy friction factor and Reynolds number is f_D·Re ≈ 62.19 (Fanning 15.55). Hence:

With a uniform inlet the computed Δp will be higher than the developed value: the entrance region adds a lumped loss due to core acceleration and higher wall shear. For a clean comparison with theory measure the gradient between two sections in the developed region, for example x = 0.10 m and x = 0.18 m, not between inlet and outlet. This is a genuine verification check, far more informative than a bare “STATUS OK”.

# expected developed gradient and Δp (Python, standard library only)
fRe, Dh, U, nu, rho, L = 62.19, 2*0.01*0.005/0.015, 0.02, 1e-6, 1000, 0.08
Re = U*Dh/nu
dpdx = fRe/Re/Dh * 0.5*rho*U**2          # Pa/m
print(Re, dpdx, dpdx*L, dpdx*L/rho)     # Re, Pa/m, Pa on L, m2/s2 on L

Cells are 2.5 mm in x and 0.625 mm in y and z: aspect ratio 4, aligned with the flow direction, where gradients are weak. There are 8 cells across the short side (5 mm) and 16 across the long side (10 mm): for a nearly parabolic profile this resolves the shape, but a quantitative comparison with the maximum velocity (1.99·U) or the pressure gradient requires the convergence study of Guide 4, for example with 40×8×4, 80×16×8 and 160×32×16 cells (ratio r = 2).

The temptation to save cells with a two-dimensional model leads here to gross errors, because the depth (5 mm) is smaller than the height (10 mm) and the side walls dominate friction. Developed pressure gradient at the same mean velocity (0.02 m/s, μ = 10⁻³ Pa·s):

A 2D model is legitimate only when the neglected dimension is much larger than the others (section aspect ratio of at least 8–10) or when the problem is truly planar by symmetry. It is a modelling decision to be justified with numbers, not a computational shortcut.

The geometric type written in constant/polyMesh/boundary constrains the allowed conditions. empty, symmetryPlane, wedge and cyclic patches are constraint types: the field type must be the same, otherwise the solver stops. The wall type imposes nothing on velocity and pressure fields, but it is required by functions that look for walls: wall distance, RAS wall functions, y⁺ and wallShearStress computation. A wall declared as patch works in laminar flow but silently breaks a turbulent case.

foamDictionary constant/polyMesh/boundary -entry walls/type -value 2>/dev/null || grep -A3 walls constant/polyMesh/boundary
foamDictionary 0/U -entry boundaryField -keywords
foamDictionary 0/p -entry boundaryField -keywords

Exact names and required entries of each condition are checked in the API documentation of your version or with foamInfo (Part 8 of the path).

foamInfo flowRateInletVelocity
foamInfo pressureInletOutletVelocity
foamInfo totalPressure

For an inlet varying in time (gradual start-up, pulsation) use a condition that reads a function of time, for example an interpolated table; for Robin-type conditions combining value and gradient there is mixed. Example of a 10 s inlet ramp, for transient cases only:

inlet
{
    type            uniformFixedValue;
    uniformValue    table
    (
        (0    (0    0 0))
        (10   (0.02 0 0))
    );
}

The syntax of time functions (table, tableFile, sine, …) is checked with foamInfo uniformFixedValue. In a steady case a ramp has no physical meaning: “time” counts iterations.

If no patch fixes pressure (for example a closed domain, or a duct with both inlet and outlet at fixed velocity), the pressure level is undetermined and must be fixed in one cell. In the algorithm dictionary add pRefCell (or pRefPoint) and pRefValue. In a closed domain the net imposed flow rate must be exactly zero, otherwise the pressure equation has no solution.

PIMPLE
{
    // ...
    pRefPoint   (0.1 0.005 0.0025);
    pRefValue   0;
}

/* OpenFOAM Foundation 14 - Guide 7 (original educational case) */
FoamFile
{
    version 2.0;
    format ascii;
    class volVectorField;
    location "0";
    object U;
}

dimensions      [0 1 -1 0 0 0 0];
internalField   uniform (0 0 0);
boundaryField
{
    inlet
    {
        type fixedValue;
        value uniform (0.02 0 0);
    }
    outlet
    {
        type zeroGradient;
    }
    walls
    {
        type noSlip;
    }
}

/* OpenFOAM Foundation 14 - Guide 7 (original educational case) */
FoamFile
{
    version 2.0;
    format ascii;
    class volScalarField;
    location "0";
    object p;
}

// Incompressible kinematic pressure: p_phys/rho, in m2/s2.
dimensions      [0 2 -2 0 0 0 0];
internalField   uniform 0;
boundaryField
{
    inlet
    {
        type zeroGradient;
    }
    outlet
    {
        type fixedValue;
        value uniform 0;
    }
    walls
    {
        type zeroGradient;
    }
}

Dimensions follow OpenFOAM’s SI order [mass length time temperature amount current luminosity]: [0 1 -1 0 0 0 0] is m/s, [0 2 -2 0 0 0 0] is m²/s². OpenFOAM checks the dimensional consistency of every operation: a pressure written in Pa [1 -1 -2 0 0 0 0] in a kinematic solver produces a dimension error, not a silent wrong result.

Starting from zero velocity is legitimate but slow: the first iterations must “invent” a field satisfying continuity. Two practical alternatives. The first is to initialise the internal field with the mean velocity, internalField uniform (0.02 0 0);, which already honours the flow rate. The second is to solve a potential flow satisfying ∇·U = 0 with the same boundary conditions and use it as the initial field:

potentialFoam -help
potentialFoam -writep > logs/log.potentialFoam 2>&1   # writes U (and p) into 0/

potentialFoam overwrites the fields in 0/: run it on the working copy, never on the template, and record the fact. In steady cases the final solution must not depend on initialisation; if it does, the problem has multiple solutions or is not converged.

The template’s uniform inlet is a convenience: it creates a development region of about 44 mm in which pressure gradient and wall stress exceed the developed value. If the goal is the behaviour of the channel within a network, impose the developed profile or extend the domain upstream and measure only downstream. A practical option, in versions that offer it, is the profile parameter of flowRateInletVelocity (laminar or turbulent boundary-layer profiles); in any case the choice must be declared on the case sheet.

/* OpenFOAM Foundation 14 - Guide 7 (original educational case) */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "constant";
    object physicalProperties;
}

viscosityModel constant;
nu             1e-06;  // m2/s; example fluid

/* OpenFOAM Foundation 14 - Guide 7 (original educational case) */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "constant";
    object momentumTransport;
}

simulationType laminar;

/* OpenFOAM Foundation 14 - Guide 7 (original educational case) */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object controlDict;
}

solver           incompressibleFluid;
startFrom        startTime;
startTime        0;
stopAt           endTime;
endTime          700;
deltaT           1;           // pseudo-time iteration for steadyState
writeControl     timeStep;
writeInterval    100;
purgeWrite       0;
writeFormat      ascii;
writePrecision   9;
writeCompression off;
timeFormat       general;
timePrecision    6;
runTimeModifiable true;

functions
{
    inletFlux
    {
        type surfaceFieldValue;
        libs ("libfieldFunctionObjects.so");
        log true;
        writeControl writeTime;
        writeFields false;
        patch inlet;
        operation sum;
        fields (phi);
    }
    outletFlux
    {
        type surfaceFieldValue;
        libs ("libfieldFunctionObjects.so");
        log true;
        writeControl writeTime;
        writeFields false;
        patch outlet;
        operation sum;
        fields (phi);
    }
}

/* OpenFOAM Foundation 14 - Guide 7 (original educational case) */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object fvSchemes;
}

ddtSchemes
{
    default steadyState;
}
gradSchemes
{
    default Gauss linear;
}
divSchemes
{
    default none;
    div(phi,U) bounded Gauss linearUpwind grad(U);
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
    default Gauss linear corrected;
}
interpolationSchemes
{
    default linear;
}
snGradSchemes
{
    default corrected;
}
wallDist
{
    method meshWave;
}

/* OpenFOAM Foundation 14 - Guide 7 (original educational case) */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object fvSolution;
}

solvers
{
    p
    {
        solver GAMG;
        smoother DIC;
        tolerance 1e-8;
        relTol 0.05;
    }
    pFinal
    {
        $p;
        relTol 0;
    }
    U
    {
        solver smoothSolver;
        smoother symGaussSeidel;
        tolerance 1e-10;
        relTol 0.1;
    }
    UFinal
    {
        $U;
        relTol 0;
    }
}
PIMPLE
{
    momentumPredictor no;
    nOuterCorrectors 1;
    nNonOrthogonalCorrectors 0;
}
relaxationFactors
{
    fields { p 0.3; }
    equations { U 0.7; }
}

/* OpenFOAM Foundation 14 - Guide 7 (original educational case) */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object blockMeshDict;
}

scale 1;
vertices
(
    (0 0 0) (0.20 0 0) (0.20 0.01 0) (0 0.01 0)
    (0 0 0.005) (0.20 0 0.005) (0.20 0.01 0.005) (0 0.01 0.005)
);
blocks ( hex (0 1 2 3 4 5 6 7) (80 16 8) simpleGrading (1 1 1) );
edges ();
boundary
(
    inlet
    {
        type patch;
        faces ((0 4 7 3));
    }
    outlet
    {
        type patch;
        faces ((1 2 6 5));
    }
    walls
    {
        type wall;
        faces ((0 1 5 4) (3 7 6 2) (0 3 2 1) (4 5 6 7));
    }
);
mergePatchPairs ();

Parts 5 and 6 of the path explain schemes and solvers in detail. To know which entries your installation uses in steady cases, compare with the pitzDailySteady tutorial: it is the safest reference for version-14 keywords.

• constant/momentumTransport: simulationType RAS; and a sub-dictionary RAS { model kOmegaSST; turbulence on; printCoeffs on; }.
• New fields 0/k, 0/omega, 0/nut, with consistent inlet conditions and wall functions (kqRWallFunction, omegaWallFunction, nutkWallFunction or low-Reynolds equivalents).
• Schemes for div(phi,k) and div(phi,omega), solvers and relaxation for k and omega.
• A wall mesh consistent with the chosen treatment (Guide 6, y⁺).

Typical inlet values: k = 1.5·(I·U)² with intensity I of a few percent and ω = √k / (C_μ^¼·ℓ) with C_μ = 0.09 and a length scale ℓ linked to the hydraulic diameter (for example 0.07·Dₕ). At Re ≈ 133 none of this is justified: a turbulence model in a laminar flow adds spurious viscosity and corrupts Δp. The correct starting point, when needed, is your installation’s RAS tutorial.

With fluids such as blood, polymers or slurries the viscosity depends on the strain rate. In Foundation 14 the model is selected in constant/physicalProperties with viscosityModel; available names and required coefficients come from the installed tutorials, not from examples for other versions:

grep -rh --include=physicalProperties "viscosityModel" "$FOAM_TUTORIALS" | sort | uniq -c
grep -rl --include=physicalProperties -E "CrossPowerLaw|BirdCarreau|powerLaw|HerschelBulkley" "$FOAM_TUTORIALS" | head

With variable viscosity this guide’s analytic solution no longer holds, the profile flattens (shear-thinning fluids) and the Reynolds number must be defined with a declared reference viscosity. This step comes after the Newtonian case has been verified.

With the classic factors p 0.3 and U 0.7 the SIMPLE algorithm converges robustly but slowly. The “consistent” variant (SIMPLEC), where available, allows higher factors (for example U 0.9 and no pressure relaxation) and considerably reduces the iteration count. The criterion for changing is empirical but measurable: compare the iterations needed to bring residuals below threshold and check that the final solution (Δp, Q) does not change. Version-14 keywords are checked in pitzDailySteady.

#!/usr/bin/env bash
# GNU GPL v3 or later - Guide 7. Never changes the shipped template.
set -Eeuo pipefail
ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)"
TEMPLATE="$ROOT/examples/channel_case"
OUTPUT="$ROOT/results"
MODE="${1:-}"
case "$MODE" in
  velocity|flowrate|viscosity) ;;
  *) echo "Usage: $0 velocity|flowrate|viscosity" >&2; exit 2;;
esac
for executable in blockMesh checkMesh foamRun; do
  command -v "$executable" >/dev/null 2>&1 || { echo "Missing $executable: source OpenFOAM-14/etc/bashrc" >&2; exit 127; }
done
mkdir -p "$OUTPUT"
DEST="$(mktemp -d "$OUTPUT/${MODE}_$(date -u +%Y%m%dT%H%M%SZ)_XXXXXX")"
cp -a "$TEMPLATE/." "$DEST/"
mkdir -p "$DEST/logs"
if [[ "$MODE" == flowrate ]]; then
  cp "$ROOT/examples/alternatives/U_flowrate" "$DEST/0/U"
elif [[ "$MODE" == viscosity ]]; then
  cp "$ROOT/examples/alternatives/physicalProperties_high_nu" "$DEST/constant/physicalProperties"
fi
printf '%s\n' "Mode: $MODE" "Case: $DEST" "UTC: $(date -u +%Y-%m-%dT%H:%M:%SZ)" >"$DEST/logs/RUN_INFO.txt"
CURRENT_STEP=initialisation
onexit() {
  status=$?
  if (( status == 0 )); then
    printf 'OK\n' >"$DEST/logs/STATUS.txt"
  else
    printf 'FAILED during %s (exit code %s)\n' "$CURRENT_STEP" "$status" >"$DEST/logs/STATUS.txt"
  fi
}
trap onexit EXIT
CURRENT_STEP=blockMesh
blockMesh -case "$DEST" >"$DEST/logs/log.blockMesh" 2>&1
CURRENT_STEP=checkMesh
checkMesh -case "$DEST" -allGeometry -allTopology >"$DEST/logs/log.checkMesh" 2>&1
CURRENT_STEP=foamRun
foamRun -case "$DEST" >"$DEST/logs/log.foamRun" 2>&1
printf 'Completed in %s\n' "$DEST"

• mktemp -d … _XXXXXX creates a directory with a random suffix: two launches in the same second never collide.
• The trap onexit EXIT writes STATUS.txt in every case, even when set -e stops the script: the file says at which stage and with which code it stopped.
• Variants are obtained by replacing a single file in the copy (0/U or constant/physicalProperties): the template is never touched.
• Exit codes: 2 for an unknown mode, 127 if an OpenFOAM executable is missing; otherwise the failed utility’s code, reported in STATUS.txt.

#!/usr/bin/env python3
# GNU GPL v3 or later - static structural checks, no OpenFOAM runtime required.
from pathlib import Path
import re, sys
BASE = Path(__file__).resolve().parents[1]
CASE = BASE/'examples/channel_case'
check = [('0/U','volVectorField','U'),('0/p','volScalarField','p'),('constant/physicalProperties','dictionary','physicalProperties'),('constant/momentumTransport','dictionary','momentumTransport'),('system/blockMeshDict','dictionary','blockMeshDict'),('system/controlDict','dictionary','controlDict'),('system/fvSchemes','dictionary','fvSchemes'),('system/fvSolution','dictionary','fvSolution')]
for relative,kind,name in check:
    s=(CASE/relative).read_text()
    assert re.search(r'\bclass\s+'+kind+r'\s*;',s),relative
    assert re.search(r'\bobject\s+'+name+r'\s*;',s),relative
    assert s.count('{') == s.count('}'),relative
    assert s.count('(') == s.count(')'),relative
for fname in ('0/U','0/p'):
    s=(CASE/fname).read_text()
    for patch in ('inlet','outlet','walls'):
        assert re.search(r'\b'+patch+r'\s*\{',s),f'{fname} missing patch {patch}'
u=(CASE/'0/U').read_text()
p=(CASE/'0/p').read_text()
assert 'dimensions      [0 1 -1 0 0 0 0]' in u
assert 'dimensions      [0 2 -2 0 0 0 0]' in p
assert 'solver           incompressibleFluid' in (CASE/'system/controlDict').read_text()
A=.01*.005; U=.02; nu=1e-6; D=2*.01*.005/(.01+.005)
print(f'PASS: 8 dictionaries; patches match; 80*16*8=10240 cells; A={A:.8g} m2; Q={A*U:.8g} m3/s; Re_Dh={U*D/nu:.2f}')
print('WARNING: static check only; does not run blockMesh/checkMesh/foamRun.')

It also runs on Windows, before the transfer: it checks class and name of each dictionary, bracket balance, presence of the three patches in the fields, the dimensions of U and p and the selected solver, and prints reference values (area, flow rate, Reynolds). It does not replace execution: a syntactically correct dictionary can be physically wrong.

python scripts/preflight.py

The three exercises are independent and short: they can be launched in sequence in a tmux session, with a final summary. Each variant creates its own directory, so one failure does not stop the others if set -e is avoided in the loop:

tmux new -s guide7
for m in velocity flowrate viscosity; do ./scripts/run_case.sh "$m" || echo "$m failed"; done
for d in results/*/; do printf '%s %s\n' "$d" "$(cat "$d/logs/STATUS.txt")"; done
python3 scripts/collect_fluxes.py results > flux_summary.csv

/* OpenFOAM Foundation 14 - Guide 7 (original educational case) */
FoamFile
{
    version 2.0;
    format ascii;
    class volVectorField;
    location "0";
    object U;
}

dimensions      [0 1 -1 0 0 0 0];
internalField   uniform (0 0 0);
boundaryField
{
    inlet
    {
        type flowRateInletVelocity;
        volumetricFlowRate 1e-06; // m3/s (positive into domain)
    }
    outlet { type zeroGradient; }
    walls  { type noSlip; }
}

/* OpenFOAM Foundation 14 - Guide 7 (original educational case) */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "constant";
    object physicalProperties;
}

viscosityModel constant;
nu             2e-06;  // m2/s

With flowRateInletVelocity and a uniform profile the imposed velocity is Q/A = 1·10⁻⁶ / 5·10⁻⁵ = 0.02 m/s: identical to exercise A. The A–B comparison is therefore a consistency test: the two solutions must coincide within solver tolerance. A noticeable difference points to an area, sign or unit error.

Doubling ν at fixed velocity, in the developed region the pressure gradient doubles (about 14 → 28 Pa/m), because f_D = 62.19/Re and Re halves; the entrance length halves (44 → 22 mm); the developed profile and hence the maximum velocity stay the same, because in laminar flow the profile shape does not depend on ν. The total inlet–outlet Δp grows less than proportionally, because the entrance contribution grows less than the viscous part.

If the run shows a 2.00 ratio between gradients measured in the developed region, the case is verified; if it shows a ratio of 2 on the total Δp, you are probably comparing different quantities or the solution is not converged.

• D — Outlet position. Extend the channel to 0.30 m (vertices and x cell count in proportion, 120 × 16 × 8) keeping everything else: the gradient measured between 0.10 and 0.18 m must not change. If it does, the outlet was affecting the solution.
• E — Refinement. Repeat A with 40 × 8 × 4 and 160 × 32 × 16 cells (r = 2) and apply the Guide 4 procedure to the developed pressure gradient. Compare the extrapolated value with 14.0 Pa/m: a complete verification, numerical and analytic together. The finest mesh has 81,920 cells, still light.

The case runs 700 iterations without a stop criterion. In serious use a residual control is added, so the run stops when the initial residuals of all equations fall below declared thresholds. In version 14 the entry goes into the algorithm dictionary; check the exact syntax in the pitzDailySteady tutorial:

grep -n -A6 residualControl "$FOAM_TUTORIALS/incompressibleFluid/pitzDailySteady/system/fvSolution"

Residuals alone are not enough: the history of the quantities of interest goes alongside. Here the natural choice is to add mean-pressure measurements on two sections in the developed region and check that their difference stays constant over the last few hundred iterations. An example function object to add to functions in the working copy:

pSections
{
    type            surfaces;
    libs            ("libsampling.so");
    writeControl    writeTime;
    surfaceFormat   raw;
    fields          (p U);
    interpolationScheme cellPoint;
    surfaces
    {
        x010 { type cuttingPlane; point (0.10 0 0); normal (1 0 0); interpolate true; }
        x018 { type cuttingPlane; point (0.18 0 0); normal (1 0 0); interpolate true; }
    }
}

Sampling-surface keywords have changed between versions: before using it compare with foamInfo surfaces or the files in $FOAM_ETC/caseDicts/functions. Alternatively, section averages are computed in ParaView (Guide 3).

At every iteration the log prints a “time step continuity errors” line with three numbers: mean local error (sum local), global and cumulative. They are normalised errors on flux divergence after the pressure correction: they must become very small and the “global” one must oscillate around zero without drift. A drifting cumulative value indicates systematic mass loss, typically from too loose pressure tolerances.

grep "continuity errors" logs/log.foamRun | tail -n 3
grep -E "^Time =|Solving for (Ux|p)" logs/log.foamRun | tail -n 8

In ParaView, Plot Over Line from (0.15, 0, 0.0025) to (0.15, 0.01, 0.0025) samples the profile across the height at mid-depth, in the developed region. The following script (standard library only) computes the series solution at the same 21 points; the RMSE between the two profiles, normalised by U, measures the local discretisation error (Guide 3).

#!/usr/bin/env python3
"""Analytic fully developed laminar profile in a rectangular duct (series solution).
Prints u/U_mean along the 10 mm side at mid-depth, to compare with a ParaView Plot Over Line CSV."""
import math
H, W, U = 0.010, 0.005, 0.02          # duct height (y), depth (z), mean velocity
a, b = W / 2, H / 2                   # half-widths: a short side, b long side
s = sum(math.tanh((2*k+1)*math.pi*b/(2*a)) / (2*k+1)**5 for k in range(200))
fac = 1 - 192*a/(math.pi**5*b)*s
G = 3*U/(a*a*fac)                      # (-dp/dx)/mu, from the mean velocity

def u(yy, zz, N=60):                  # yy across the short side, zz across the long side
    out = 0.0
    for k in range(N):
        n = 2*k + 1
        arg = n*math.pi*b/(2*a)
        r = math.cosh(n*math.pi*zz/(2*a))/math.cosh(arg) if arg < 700 else math.exp(n*math.pi*(abs(zz)-b)/(2*a))
        out += (-1)**k * (1 - r) * math.cos(n*math.pi*yy/(2*a)) / n**3
    return 16*a*a/math.pi**3 * G * out

print("y_m,u_over_Umean")
for i in range(21):
    y = H*i/20                        # 0 ... 10 mm, wall to wall
    print(f"{y:.4f},{u(0.0, y - b)/U:.4f}")

The profile along the long side is flatter at the centre than a parabola, because at mid-depth the gradient is dominated by the nearby walls (short side). A Poiseuille parabola between parallel plates would give a maximum of 1.5·U: a result close to 1.5 would indicate a geometry error or an unintentionally two-dimensional case.

A fatal error reports the message, the file and often the line of the dictionary involved, and the code function that raised it. The useful message is the one before “From function”; the stack lines that follow are for developers. The most frequent errors at this stage: “Cannot find patchField entry for walls” (different patch name), “keyword nu is undefined” (dictionary from another release), “Different dimensions for (p + …)” (wrong units), “Unknown patchField type” (misspelt condition or not available in the version).

grep -n -m1 -B2 -A12 "FOAM FATAL" logs/log.foamRun
cat logs/STATUS.txt

• Read in which equation the residual blows up first (Ux or p) and at which iteration.
• Re-run checkMesh and compare with Guide 6: a freshly modified mesh is the first suspect.
• Reduce relaxation factors (for example p 0.2, U 0.5) and temporarily switch to bounded Gauss upwind for div(phi,U): if it then converges, the problem is numerical; if it still diverges, it is in the boundary conditions.
• Restore the original settings starting from the converged solution (startFrom latestTime): the second-order scheme must be the one used for the final result.

The same scheme as Guide 5 applies here: one archive per variant, without the intermediate fields that can be regenerated, with file fingerprints and an environment description.

for d in results/*/; do
  ( cd "$d" && { date -u; foamVersion; } > logs/environment.txt 2>&1
    find 0 constant system logs postProcessing -type f -print0 | sort -z | xargs -0 sha256sum > logs/files.sha256 )
  tar czf "$(basename "$d").tar.gz" -C "$d" 0 constant/physicalProperties constant/momentumTransport system logs postProcessing
done

#!/usr/bin/env python3
# GPL-3.0-or-later. Parse OpenFOAM 14 surfaceFieldValue output on Linux or Windows.
"""Usage: python scripts/collect_fluxes.py results > flux_summary.csv
Read last numeric time row from postProcessing/{inletFlux,outletFlux}/*/surfaceFieldValue.dat.
Missing values remain blank; script never fabricates simulation outputs.
"""
from pathlib import Path
from csv import writer
from sys import argv, stdout, stderr
from math import isfinite

def read_series(case, name):
    rows={}
    for f in (case/'postProcessing'/name).glob('*/surfaceFieldValue.dat'):
        with f.open(encoding='utf-8',errors='replace') as inp:
            for line in inp:
                if not line.strip() or line.lstrip().startswith('#'):
                    continue
                parts=line.split()
                if len(parts)<2:continue
                try: x=float(parts[0]);y=float(parts[1])
                except ValueError:continue
                if isfinite(x) and isfinite(y):
                    old=rows.get(x)
                    if old is None or f.stat().st_mtime >= old[0]:
                        rows[x]=(f.stat().st_mtime,y)
    return {t:value for t,(_,value) in rows.items()}

def main():
    if len(argv)!=2:
        print('usage: collect_fluxes.py <results-directory>',file=stderr);return 2
    root=Path(argv[1]);
    if not root.is_dir():
        print(f'not a directory: {root}',file=stderr);return 2
    cw=writer(stdout,lineterminator='\n')
    cw.writerow(['case','status','Q_inlet_m3s','Q_outlet_m3s','imbalance_relative'])
    for case in sorted(p for p in root.iterdir() if p.is_dir()):
        status_file=case/'logs'/'STATUS.txt'
        status=status_file.read_text(encoding='utf-8').strip() if status_file.exists() else 'UNKNOWN'
        ins=read_series(case,'inletFlux');outs=read_series(case,'outletFlux')
        common=ins.keys() & outs.keys()
        t=max(common) if common else None
        qi=ins.get(t);qo=outs.get(t)
        imb=abs(qi+qo)/max(abs(qi),abs(qo)) if qi is not None and qo is not None and max(abs(qi),abs(qo))>0 else None
        cw.writerow([case.name,status,qi if qi is not None else '',qo if qo is not None else '',imb if imb is not None else ''])
    return 0
if __name__=='__main__':raise SystemExit(main())

• It keeps one value per time; if the same time appears in several files (restarts), the most recently modified file wins.
• It uses the last time common to both measurements, so it never compares different times.
• The imbalance is relative: |Qin + Qout| / max(|Qin|, |Qout|); if a value is missing the cell stays empty.

With the nominal case one expects Q_in ≈ −1·10⁻⁶ and Q_out ≈ +1·10⁻⁶ m³/s; the relative imbalance of a converged incompressible case is of the order of the pressure tolerance, well below 10⁻⁶.

In the developed region the force balance on a volume of length L gives Δp·A = τ̄_w·P·L, with P the wetted perimeter (30 mm). It is a second independent check: with 14 Pa/m and A/P = Dₕ/4 = 1.667 mm the expected mean wall stress is τ̄_w ≈ 0.023 Pa. In OpenFOAM wall shear stress is computed with the wallShearStress function object (kinematic, multiply by ρ) and averaged over the patch in ParaView or with surfaceFieldValue.

foamPostProcess -help
foamPostProcess -func wallShearStress -latestTime   # on the working copy, after the run

A and B impose the same inlet field, so the difference between their solutions measures only iterative noise. A reasonable engineering criterion: the difference in the quantities of interest (Δp between the two sections, Q) must be at least an order of magnitude smaller than any effect you want to attribute to physics, for example the A–C difference. If A and B differ by more than 1 %, iterative convergence is insufficient and the A–C comparison is not yet interpretable.

case,inlet_bc,nu_m2_s,Re_Dh,Q_in_m3_s,Q_out_m3_s,imbalance_rel,dpdx_developed_Pa_m,dp_total_Pa,umax_m_s,iterations,notes
A_velocity,fixedValue 0.02,1e-6,133.3,,,,,,,,
B_flowrate,flowRateInletVelocity 1e-6,1e-6,133.3,,,,,,,,
C_viscosity,fixedValue 0.02,2e-6,66.7,,,,,,,,

Expected values for verification: dpdx_developed ≈ 14.0 / 14.0 / 28.0 Pa/m; umax ≈ 0.0398 m/s in all three cases. Empty columns are filled with data from collect_fluxes.py and ParaView; units are in the column name, so the file stays readable out of context.

• R. K. Shah, A. L. London, Laminar Flow Forced Convection in Ducts, Academic Press, 1978. Solutions and correlations for ducts of any section, including f·Re for rectangles.
• F. M. White, Viscous Fluid Flow, McGraw-Hill. Laminar duct flow, entrance length, boundary layer.
• S. B. Pope, Turbulent Flows, Cambridge University Press, 2000. Foundations of the RAS and LES models mentioned in section 6.
• Tutorials in your installation: incompressibleFluid/pitzDailySteady, incompressibleFluid/cylinder, to compare entry by entry with the guide’s case.