OpenFOAM: turbulence, k-ε and k-ω SST RANS models and near-wall treatment

With the tutorial data (U = 10 m/s, ν = 1·10⁻⁵ m²/s) and a step height and inlet-channel height of 25.4 mm (to be checked in the mesh dictionary, section 3):

The last number is the most instructive: already at the inlet the turbulent viscosity is about 85 times the molecular one. In a RANS computation the flow behaviour therefore depends decisively on how the model produces and dissipates ν_t: this is why the choice of model shifts reattachment measurably.

Reynolds decomposition writes every quantity as mean plus fluctuation; averaging the equations produces the Reynolds-stress tensor, which must be modelled. With the Boussinesq hypothesis:

u = Ū + u'          (mean + fluctuation)
−⟨u'ᵢ u'ⱼ⟩ = ν_t (∂Ūᵢ/∂xⱼ + ∂Ūⱼ/∂xᵢ) − (2/3) k δᵢⱼ
k = ½ ⟨u'ᵢ u'ᵢ⟩

All the physics of the fluctuations is compressed into a single scalar, ν_t. Practical consequences: normal stresses come out isotropic apart from the mean gradient, curvature and rotation effects are not represented naturally, and the separation zone behind the step, where turbulence is strongly anisotropic, is exactly where two-equation models differ most.

Before preparing the cases, write in a text file inside the campaign folder which quantities you want to compare and how you will measure them: for example “reattachment length on the lower wall, defined by the sign change of the x component of wall shear stress” and “ΔP between sections x = −0.02 m and x = 0.2 m, area-weighted average”. Deciding afterwards, while looking at the results, inevitably leads to choosing the definition that confirms the initial idea.

For a DNS the ratio between the largest scale and the Kolmogorov scale grows as Re^(3/4) in each direction, so the number of 3D points grows as Re^(9/4). With Re ≈ 2.5·10⁴:

N_DNS ~ Re^(9/4) = 25400^2.25 ≈ 8·10⁹ points (order of magnitude)
N_RANS 2D (pitzDaily) ≈ 1·10⁴ cells
ratio ≈ 10⁶

Moreover DNS is necessarily three-dimensional and transient, with thousands of time steps to obtain converged statistics, whereas steady RANS converges in a few thousand iterations on a 2D mesh. LES sits in between: it resolves the large eddies, but near walls the cost again grows rapidly with Re if the boundary layer is resolved. This gap of six or more orders of magnitude is what makes RANS the everyday tool of engineering, with the limits the guide highlights.

• Steady RANS: the mesh must resolve mean-field gradients; refining it improves the solution of the RANS equations, but does not bring the result closer to real turbulence beyond what the model allows;
• URANS: represents only large-scale instabilities clearly separated from the turbulent spectrum (for example periodic shedding); turbulent fluctuations remain modelled by ν_t;
• LES: the filter size is tied to the mesh, so the mesh is part of the model; a spatial convergence study has a different meaning than in RANS and requires time statistics.

The case is named after the Pitz and Daily (1983) experiment on a turbulent mixing layer downstream of a rearward-facing step. The dimensions can be read directly from the resource dictionary:

D="$FOAM_TUTORIALS/resources/blockMesh/pitzDaily"
foamDictionary "$D" -entry convertToMeters
foamDictionary "$D" -entry vertices
foamDictionary "$D" -entry blocks

From the vertex coordinates you obtain the inlet-channel height, the step height and the expansion ratio, which should be recorded in the project notebook before any comparison. In the literature on steps with an expansion ratio around 2 and Re of order 10⁴ the measured reattachment length is typically a few step heights (of order 6–8 h), and standard k-ε is known to tend to underpredict it because it produces too much turbulent viscosity in the shear layer. Use these values as an order of magnitude to recognise an anomalous result, not as validation data: a real comparison needs the original measurements and the same inlet conditions.

T="$FOAM_TUTORIALS/incompressibleFluid/pitzDailySteady"
cat "$T/Allrun"
ls "$T" "$T/0" "$T/constant" "$T/system"
foamDictionary "$T/system/fvSolution" -entry SIMPLE
foamDictionary "$T/system/controlDict" -entry endTime
foamDictionary "$T/0/U" -entry boundaryField

The official Allrun is the source of truth for the commands: if a new version changes the sequence (for example the mesh resource path), the guide’s script must be adapted to it, not the other way round. The SIMPLE block shows the stopping criteria (residualControl) with which the solver declares convergence.

The ν_t limiter is the difference that matters most in this case: in the shear layer starting from the step edge k-ε produces high turbulent viscosity, which mixes rapidly and shortens the recirculation bubble; SST limits it and tends to produce longer bubbles. The coefficients actually used are printed at the start of log.foamRun (block kEpsilonCoeffs or kOmegaSSTCoeffs): keep them as part of the run documentation.

R="$HOME/OpenFOAM/guide09_run01"
diff -ru "$R/kEpsilon" "$R/kOmegaSST"
foamDictionary "$R/kEpsilon/constant/momentumTransport" -entry RAS/model
foamDictionary "$R/kOmegaSST/constant/momentumTransport" -entry RAS/model

Before launching, the diff must show only the model line and the GUIDE09_MODEL.txt file. Save the output in a campaign file (diff -ru ... > diff_before_run.txt): it is the documentary proof that the comparison is controlled.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Calculate indicative RANS inlet k/epsilon/omega from U, intensity and length scale."""
import argparse
from math import sqrt


def estimates(U, intensity, length, cmu=0.09):
    if not (U > 0 and 0 < intensity < 1 and length > 0 and cmu > 0):
        raise ValueError("U>0, 0<I<1, length>0 and Cmu>0 are required")
    k = 1.5 * (U * intensity)**2
    epsilon = cmu**0.75 * k**1.5 / length
    omega = sqrt(k) / (cmu**0.25 * length)
    return k, epsilon, omega


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument("--U", type=float, default=10.0, help="Mean inlet speed in m/s")
    p.add_argument("--I", type=float, default=0.05, help="Turbulence intensity as a fraction")
    p.add_argument("--L", type=float, default=0.00254, help="Turbulence length scale in m")
    a = p.parse_args()
    try:
        k, epsilon, omega = estimates(a.U, a.I, a.L)
    except ValueError as ex:
        p.error(str(ex))
    print(f"k       = {k:.8g} m2/s2")
    print(f"epsilon = {epsilon:.8g} m2/s3")
    print(f"omega   = {omega:.8g} 1/s")
    print("Indicative estimates; verify consistency with the selected model and experiment.")

if __name__ == "__main__":
    main()

$ python3 scripts/estimate_inlet.py --U 10 --I 0.05 --L 0.00254
k       = 0.375 m2/s2
epsilon = 14.855767 m2/s3
omega   = 440.17086 1/s
Indicative estimates; verify consistency with the selected model and experiment.

The three relations are built so that both models start from the same turbulent viscosity: ν_t = Cμ k²/ε = 0.09 · 0.375² / 14.856 ≈ 8.52·10⁻⁴ m²/s and ν_t = k/ω = 0.375 / 440.17 ≈ 8.52·10⁻⁴ m²/s. This is the check to make every time the inlet values are changed: if the two ν_t differ, the two models start from different conditions and the comparison is no longer controlled.

The length ℓ = 2.54 mm corresponds to 10 % of the inlet-channel height; a widespread rule of thumb is ℓ ≈ 0.07 Dₕ, which with Dₕ = 50.8 mm would give 3.6 mm. Neither is “right”: they are assumptions. Since inlet ν_t changes proportionally with I and ℓ, a sensitivity campaign on these two parameters (separate from the model comparison, as the guide requires) tells you how much reattachment depends on information usually not known precisely.

For the inlet channel, upstream of the step, Dean’s correlation for turbulent plane-channel flow can be used, Cf = 0.073 Re^(−1/4) with Re based on the channel height. The following script, added to this guide, standard library only, applies the estimate:

#!/usr/bin/env python3
"""A-priori estimate of y+ for the first cell centre in a plane channel (engineering estimate).

Uses Dean's correlation Cf = 0.073 Re^-1/4 (Re on full channel height and bulk velocity),
valid roughly for 6e3 < Re < 6e5 in smooth, fully developed channel flow.
Usage: first_cell.py --U 10 --H 0.0254 --nu 1e-5 [--yp 0.25e-3 0.5e-3 1e-3] [--target 1 30 100]
Standard library only. License: GNU GPL v3 or later.
"""
from __future__ import annotations
import argparse, math


def main() -> None:
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("--U", type=float, required=True, help="bulk velocity [m/s]")
    ap.add_argument("--H", type=float, required=True, help="full channel height [m]")
    ap.add_argument("--nu", type=float, required=True, help="kinematic viscosity [m2/s]")
    ap.add_argument("--yp", type=float, nargs="*", default=[], help="first-cell-centre distances [m]")
    ap.add_argument("--target", type=float, nargs="*", default=[1.0, 30.0, 100.0], help="target y+ values")
    a = ap.parse_args()
    Re = a.U * a.H / a.nu
    Cf = 0.073 * Re ** -0.25
    utau = a.U * math.sqrt(Cf / 2.0)
    lv = a.nu / utau
    print(f"Re_H = {Re:.4g}   Cf = {Cf:.4g}   u_tau = {utau:.4g} m/s   viscous length nu/u_tau = {lv:.4g} m")
    if not 6e3 <= Re <= 6e5:
        print("warning: Re outside the usual range of Dean's correlation")
    for y in a.yp:
        print(f"y_P = {y:.4g} m  ->  y+ = {y / lv:.4g}   (first cell height ~ {2 * y:.4g} m)")
    for t in a.target:
        print(f"target y+ = {t:g}  ->  y_P = {t * lv:.4g} m  (first cell height ~ {2 * t * lv:.4g} m)")
    print("Estimate only: separated and reattaching regions have very different wall shear; measure y+ after the run.")


if __name__ == "__main__":
    main()

$ python3 first_cell.py --U 10 --H 0.0254 --nu 1e-5 --yp 0.25e-3 0.5e-3 1e-3
Re_H = 2.54e+04   Cf = 0.005782   u_tau = 0.5377 m/s   viscous length nu/u_tau = 1.86e-05 m
y_P = 0.00025 m  ->  y+ = 13.44   (first cell height ~ 0.0005 m)
y_P = 0.0005 m  ->  y+ = 26.89   (first cell height ~ 0.001 m)
y_P = 0.001 m  ->  y+ = 53.77   (first cell height ~ 0.002 m)
target y+ = 1  ->  y_P = 1.86e-05 m  (first cell height ~ 3.72e-05 m)
target y+ = 30  ->  y_P = 0.0005579 m  (first cell height ~ 0.001116 m)
target y+ = 100  ->  y_P = 0.00186 m  (first cell height ~ 0.00372 m)

The viscous length is about 19 μm: resolving the viscous sublayer (y⁺ ≈ 1) would need a first cell of about 37 μm, log-law wall functions (y⁺ ≈ 30) about 1.1 mm. Downstream of the step the estimate does not hold: inside the bubble the wall shear passes through zero, so y⁺ tends to zero near the separation and reattachment points whatever the mesh. That is why the guide asks to look at the distribution, not at a single value.

• nutkWallFunction: in the wall cell it computes ν_t from k so that wall shear follows the log law when y⁺ exceeds a threshold (of order 11, where the linear and log laws meet); below the threshold only molecular viscosity acts;
• epsilonWallFunction and omegaWallFunction: impose ε or ω in the wall cell according to the log-region scale (ω also according to the viscous one, blending the two in recent versions);
• kqRWallFunction: zero gradient of k at the wall, consistent with the functions above.

Details (thresholds, blended formulations) change between versions: the reference is the source code of your installation, browsable in the C++ documentation cited in the sources. The operational consequence remains: y⁺ must be compared with the validity range of the wall function actually chosen.

cd "$HOME/OpenFOAM/guide09_run01/kOmegaSST"
for f in U p k epsilon omega nut; do
  printf "%-8s lowerWall: %s\n" "$f" "$(foamDictionary 0/$f -entry boundaryField/lowerWall/type -value)"
done

A table like this, saved in the case folder, immediately shows a frequent error: a field with the wrong or missing condition on a patch. Repeat it for upperWall, inlet and outlet.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Prepare isolated OpenFOAM Foundation 14 pitzDailySteady RANS cases."""
import argparse
import csv
import os
from pathlib import Path
import shutil
import sys

MODELS = {"kEpsilon": "kEpsilon", "kOmegaSST": "kOmegaSST"}


def prepare(template: Path, output: Path):
    template = template.expanduser().resolve(strict=True)
    output = output.expanduser().resolve()
    if not template.is_dir():
        raise ValueError("Tutorial template is not a directory")
    for required in ("0/U", "0/k", "0/epsilon", "0/omega", "0/nut",
                     "system/controlDict", "system/fvSchemes", "system/fvSolution",
                     "constant/momentumTransport", "constant/physicalProperties"):
        if not (template / required).is_file():
            raise ValueError(f"Missing tutorial file: {required}")
    if output.exists():
        raise FileExistsError(f"Destination exists (will not overwrite): {output}")
    if output == template or template in output.parents and output == template:
        raise ValueError("Output cannot be the template")
    output.mkdir(parents=True, exist_ok=False)
    try:
        for case_name, model in MODELS.items():
            dst = output / case_name
            shutil.copytree(template, dst, symlinks=False,
                            ignore=shutil.ignore_patterns("processor*", "log.*", "postProcessing"))
            path = dst / "constant/momentumTransport"
            text = path.read_text(encoding="utf-8")
            import re
            if len(re.findall(r"(?m)^\s*model\s+kEpsilon\s*;", text)) != 1:
                raise ValueError("Expected exactly one 'model kEpsilon;' in tutorial momentumTransport")
            text = re.sub(r"(?m)^(\s*model\s+)kEpsilon(\s*;)",
                          lambda m: m.group(1) + model + m.group(2), text, count=1)
            path.write_text(text, encoding="utf-8", newline="\n")
            (dst / "GUIDE09_MODEL.txt").write_text(
                f"model={model}\nsource={template}\n"+
                "Run independently; no experimental validation is implied.\n", encoding="utf-8")
        with (output / "cases.csv").open("w", encoding="utf-8", newline="") as f:
            writer = csv.writer(f)
            writer.writerow(["case", "model", "mesh", "status"])
            for name, model in MODELS.items():
                writer.writerow([name, model, "pitzDaily tutorial resource", "not_run"])
    except Exception:
        # Only clean up the new directory created by this function, never pre-existing work.
        shutil.rmtree(output)
        raise
    print(f"Created two independent cases at {output}")


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--out", type=Path, required=True)
    parser.add_argument("--template", type=Path,
                        default=Path(os.environ.get("FOAM_TUTORIALS", "")) / "incompressibleFluid/pitzDailySteady")
    args = parser.parse_args()
    if not os.environ.get("FOAM_TUTORIALS") and "--template" not in sys.argv:
        parser.error("FOAM_TUTORIALS is not set; load OpenFOAM 14 or provide --template")
    try:
        prepare(args.template, args.out)
    except (FileNotFoundError, FileExistsError, ValueError, OSError) as exc:
        parser.exit(2, f"Preparation failed: {exc}\n")

if __name__ == "__main__":
    main()

• checks that every file needed for the comparison exists in the tutorial before creating anything;
• requires exactly one model kEpsilon; line: should a future tutorial change, the script stops instead of producing a wrong comparison;
• on error it removes only the folder it has just created itself, never pre-existing work, because an existing destination is refused at the start;
• writes cases.csv with status not_run: a minimal campaign register to be updated by hand.

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Run one untouched pitzDailySteady-derived case, without overwriting logs/results.
set -Eeuo pipefail
if [[ $# -ne 1 ]]; then
  echo "Usage: $0 /absolute/path/to/prepared/case" >&2
  exit 2
fi
: "${FOAM_TUTORIALS:?Source OpenFOAM 14 environment before running}"
case_dir="$(realpath "$1")"
[[ -f "$case_dir/GUIDE09_MODEL.txt" ]] || { echo "Not a Guide 09 prepared case" >&2; exit 2; }
[[ -f "$FOAM_TUTORIALS/resources/blockMesh/pitzDaily" ]] || { echo "Missing OpenFOAM tutorial mesh resource" >&2; exit 2; }
for logfile in log.blockMesh log.checkMesh log.foamRun log.yPlus; do
  [[ ! -e "$case_dir/$logfile" ]] || { echo "Existing $logfile: refusing to overwrite" >&2; exit 2; }
done
[[ ! -d "$case_dir/constant/polyMesh" ]] || { echo "Existing mesh: refusing to alter case" >&2; exit 2; }
cd "$case_dir"
echo "Running case: $case_dir"
blockMesh -dict "$FOAM_TUTORIALS/resources/blockMesh/pitzDaily" > log.blockMesh 2>&1
checkMesh > log.checkMesh 2>&1
foamRun > log.foamRun 2>&1
# Failure in optional yPlus post-processing does not invalidate a finished flow solution.
if ! foamPostProcess -solver incompressibleFluid -func yPlus -latestTime > log.yPlus 2>&1; then
  echo "WARNING: yPlus post-processing failed; see log.yPlus" >&2
fi
echo "Completed flow solution for $case_dir"

The check on constant/polyMesh is an important protection: re-running blockMesh on an already computed case would replace the mesh, leaving incompatible saved fields. The y⁺ post-processing is separated from the run with if !: an error in it produces a warning but neither deletes nor invalidates the flow solution.

• in a steady run Time = ... is an iteration counter, not physical time;
• if residualControl is active, the solver stops before endTime with a message such as “SIMPLE solution converged in N iterations”: record N for each model;
• bounding k, bounding epsilon or bounding omega lines mean negative values clipped to a minimum: a few in the first iterations are common, persisting to the end they signal a mesh, condition or scheme problem;
• continuity errors sum local that do not decrease indicate a non-converged solution even if residuals look low.

grep -E "converged in|^End" log.foamRun
grep -c "^bounding" log.foamRun
grep "bounding" log.foamRun | tail -3

The two cases are independent: with at least two free cores they can be launched in two separate tmux windows (Ctrl+b then c for a new window). On a 2D mesh of about ten thousand cells each run uses one core and little memory; total time halves with no risk to the data, because the folders are separate. Check the load first with uptime or htop if the server is shared.

The yPlus function object writes the yPlus field into the last saved time folder and prints in the log a per-wall-patch summary with minimum, maximum and average. The summary is the first thing to look at:

grep -i "y+" log.yPlus
foamListTimes -latestTime
ls "$(foamListTimes -latestTime)"

An acceptable average with a minimum near zero and a very high maximum is not contradictory: in the recirculation bubble wall shear is small, under the inlet shear layer it is large. The assessment must be made by zones: inlet channel, region behind the step, reattachment zone, upper wall.

The same post-processing mechanism computes wall shear stress, which is needed to define reattachment reproducibly (section 11). Before using it, check the name of the function available in your installation:

foamPostProcess -list | grep -i -E "wallShear|yPlus"
foamPostProcess -solver incompressibleFluid -func wallShearStress -latestTime > log.wallShearStress 2>&1
tail -n 20 log.wallShearStress

The resulting field is per unit density [m²/s²] and has three components: for reattachment the component parallel to the lower wall matters, that is x in this geometry. This command also writes a separate log and does not modify the flow fields.

• in ParaView apply Extract Block and select only the lowerWall patch (with wallShearStress among the loaded fields);
• with a 2D one-cell-thick geometry the patch is a strip: Plot Data or Plot Over Line along x at constant z give the x component versus position;
• export to CSV with File → Save Data and compute the sign changes with the following script, added to this guide.

#!/usr/bin/env python3
"""Find wall-shear sign changes along a wall from a ParaView CSV export.

The CSV must contain a coordinate column (default 'Points:0') and the streamwise
wall-shear component (default 'wallShearStress:0'). Rows are sorted by coordinate;
each sign change is located by linear interpolation between neighbouring samples.
Values are per unit density [m2/s2] in incompressible OpenFOAM. Establish the sign
convention from a region where the flow is certainly attached (far downstream)
before naming a crossing 'separation' or 'reattachment'.
Usage: find_reattachment.py lowerWall.csv [--x Points:0] [--tau wallShearStress:0]
       [--xmin 0] [--h 0.0254]
Standard library only. License: GNU GPL v3 or later.
"""
from __future__ import annotations
import argparse, csv
from pathlib import Path


def main() -> None:
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("csv", type=Path)
    ap.add_argument("--x", default="Points:0", help="coordinate column")
    ap.add_argument("--tau", default="wallShearStress:0", help="streamwise wall shear column")
    ap.add_argument("--xmin", type=float, default=float("-inf"), help="ignore samples upstream of this x")
    ap.add_argument("--h", type=float, default=None, help="step height for x/h output")
    a = ap.parse_args()
    with a.csv.open(newline="", encoding="utf-8") as f:
        rd = csv.DictReader(f)
        for col in (a.x, a.tau):
            if col not in (rd.fieldnames or []):
                ap.error(f"column '{col}' not found; available: {rd.fieldnames}")
        pts = sorted((float(r[a.x]), float(r[a.tau])) for r in rd if float(r[a.x]) >= a.xmin)
    if len(pts) < 2:
        ap.error("not enough samples")
    found = 0
    for (x0, t0), (x1, t1) in zip(pts, pts[1:]):
        if t0 == 0 or t0 * t1 < 0:
            xc = x0 if t0 == 0 else x0 + (x1 - x0) * (-t0) / (t1 - t0)
            kind = "negative -> positive" if t1 > t0 else "positive -> negative"
            extra = f"   x/h = {xc / a.h:.4g}" if a.h else ""
            print(f"sign change at x = {xc:.6g} m ({kind}){extra}")
            found += 1
    print(f"{len(pts)} samples, {found} sign change(s). Report the criterion and the wall patch used.")


if __name__ == "__main__":
    main()

Tested on a synthetic profile with a known sign change at x = 0.17 m, the script returns:

$ python3 find_reattachment.py test_wall.csv --xmin 0.001 --h 0.0254
sign change at x = 0.17 m (negative -> positive)   x/h = 6.693
200 samples, 1 sign change(s). Report the criterion and the wall patch used.

More than one sign change is normal: near the step corner a small counter-rotating secondary bubble may exist. Report all points found and state which one you take as main reattachment, and why.

Besides its length, the bubble has a measurable area. In ParaView: Threshold on U_X with upper bound 0 (cells with backflow), then Integrate Variables: the integrated volume divided by the domain thickness in z gives the in-plane bubble area. Done with the same filter on both cases, it is an extra indicator, insensitive to a single sampling point. Always state the criterion (U_x < 0, not “blue zone of the colour map”).

In kinematic units the reference dynamic pressure is ½U² = 50 m²/s² whatever the fluid; the value in pascals depends entirely on the declared density. Note that ν = 1·10⁻⁵ m²/s corresponds exactly neither to air (about 1.5·10⁻⁵) nor to water (about 1·10⁻⁶): it is a tutorial value. This is why density must be declared as an assumption of the measurement, and a ΔP in pascals only makes sense together with it.

For inlet and outlet, which are patches, averages can be computed without ParaView using the predefined foamPostProcess configurations. Names and arguments depend on the version: list them and read the configuration file header before using them.

foamPostProcess -list | grep -i -E "patch|surface"
# example to adapt to what is listed
foamPostProcess -solver incompressibleFluid -func "patchAverage(name=inlet, p)" -latestTime
foamPostProcess -solver incompressibleFluid -func "patchAverage(name=outlet, p)" -latestTime
foamPostProcess -solver incompressibleFluid -func "patchFlowRate(patch=outlet)" -latestTime

The result goes to postProcessing/ with one file per function: the difference of the two averages, multiplied by the declared density, is the ΔP between inlet and outlet. For internal sections (for example at x = 0.2 m) you need a sampling plane or Slice + Integrate Variables in ParaView, with area-weighted averaging.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Collect traceable last residuals from Guide 09 solver logs (no invented CFD values)."""
import argparse
import csv
from pathlib import Path
import re

RESIDUAL = re.compile(
    r"Solving\s+for\s+([A-Za-z0-9_]+)\s*,\s*Initial residual\s*=\s*"
    r"([+\-0-9.eE]+)\s*,\s*Final residual\s*=\s*([+\-0-9.eE]+)"
)
FIELDS = ("p", "Ux", "Uy", "Uz", "k", "epsilon", "omega")

def summarize(case: Path):
    log = case / "log.foamRun"
    result = {"case": case.name, "model": "", "log_present": str(log.is_file()),
              "end_marker": "", "residual_records": 0}
    for f in FIELDS:
        result[f + "_initial_residual_last"] = ""
    info = case / "GUIDE09_MODEL.txt"
    if info.is_file():
        for line in info.read_text(encoding="utf-8").splitlines():
            if line.startswith("model="):
                result["model"] = line.split("=", 1)[1]
    if log.is_file():
        content = log.read_text(encoding="utf-8", errors="replace")
        result["end_marker"] = str(bool(re.search(r"(?m)^End\s*$", content)))
        for field, initial, final in RESIDUAL.findall(content):
            result["residual_records"] += 1
            key = field + "_initial_residual_last"
            if key in result:
                result[key] = initial
    return result

def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument("root", type=Path)
    p.add_argument("--output", type=Path, default=Path("summary_log.csv"))
    a = p.parse_args()
    root = a.root.expanduser().resolve()
    cases = [root / name for name in ("kEpsilon", "kOmegaSST")]
    records = [summarize(c) for c in cases]
    if a.output.exists():
        p.error(f"Refusing to replace existing output: {a.output}")
    with a.output.open("w", encoding="utf-8", newline="") as fh:
        writer = csv.DictWriter(fh, fieldnames=list(records[0]))
        writer.writeheader()
        writer.writerows(records)
    print(f"Wrote {a.output}; residuals alone do not establish convergence.")

if __name__ == "__main__":
    main()

measurement_template.csv

case,model,mesh_cells,yplus_wall_min,yplus_wall_mean,yplus_wall_max,pressure_upstream_kinematic_m2_s2,pressure_downstream_kinematic_m2_s2,density_kg_m3,pressure_drop_Pa,reattachment_length_m,flow_rate_m3_s,measurement_method,notes
kEpsilon,kEpsilon,,,,,,,,,,,,
kOmegaSST,kOmegaSST,,,,,,,,,,,,

The script reads the last initial residual per field: a snapshot, not a history. Note also that it refuses to overwrite an existing CSV, consistent with the rest of the package. The measurement sheet separates value columns from a measurement_method column: filling it in without the method makes the numbers useless for a future comparison.

• initial residuals have dropped several orders of magnitude and levelled off, not merely “below the threshold”;
• the quantities of interest (ΔP, reattachment, outlet flow rate) no longer change, within the required accuracy, over the last few hundred iterations: checking this needs function objects active during the run or several saved instants;
• outlet flow rate matches inlet flow rate (mass balance, Guide 2);
• no persistent bounding lines remain.

If one model meets these criteria and the other does not, the comparison must be suspended: you would be comparing a solution with an unfinished computation.

To follow ΔP and flow rate during the iterations, add function objects to controlDict in a new campaign (not in the one already prepared). Indicative structure, to adapt to the names listed by foamPostProcess -list:

functions
{
    #includeFunc patchAverage(name=inlet, p)
    #includeFunc patchAverage(name=outlet, p)
    #includeFunc patchFlowRate(patch=outlet)
}

The series in postProcessing/ make the criterion “the quantity of interest no longer changes” verifiable, which residuals alone do not guarantee.

If exercise B shows y⁺ out of range, the right correction is a new family of cases with a modified mesh, where first-cell thickness and wall function are chosen consistently: for example wall functions with y⁺ ≈ 30–100 throughout the inlet channel, or a first cell with y⁺ ≈ 1 and formulations suited to low y⁺. Mixing the two strategies in different meshes of the same convergence study makes the observed order (Guide 4) meaningless.

The indication “SST typically longer” comes from the known behaviour of the two models in separated flows and serves as a plausibility check, not as a result to obtain at all costs: if the opposite emerges, the interesting thing is understanding why.

Case: pitzDailySteady (OpenFOAM Foundation 14), mesh resources/blockMesh/pitzDaily, N cells = ...
Models: kEpsilon | kOmegaSST — only difference verified with diff -ru
Inlet: U = 10 m/s, I = 0.05, l = 2.54 mm, nu_t/nu = 85
Convergence: iterations ...; final residuals ...; mass balance ...
y+ lowerWall (min/mean/max): ... | ...   critical zones: ...
Reattachment (sign change of tau_x on lowerWall): x_r = ... m | ... m  (x_r/h = ... | ...)
Delta P (patchAverage inlet-outlet, rho = ... kg/m3): ... Pa | ... Pa
Observations / limits / assumptions: ...

• S. B. Pope, Turbulent Flows, Cambridge University Press, 2000: Reynolds decomposition, turbulence scales, eddy-viscosity models;
• D. C. Wilcox, Turbulence Modeling for CFD, DCW Industries: k-ε and k-ω models, wall functions, test cases;
• B. E. Launder, D. B. Spalding, The numerical computation of turbulent flows, Computer Methods in Applied Mechanics and Engineering 3 (1974): the standard k-ε model and wall functions;
• F. R. Menter, Two-equation eddy-viscosity turbulence models for engineering applications, AIAA Journal 32 (1994): the SST model;
• R. W. Pitz, J. W. Daily, Combustion in a turbulent mixing layer formed at a rearward-facing step, AIAA Journal 21 (1983): the experiment the case is named after.

cd "$HOME/OpenFOAM"
foamVersion > guide09_run01/OPENFOAM_VERSION.txt 2>&1
tar czf guide09_run01_$(date +%Y%m%d).tar.gz guide09_run01 guide09_run01_log_summary.csv
sha256sum guide09_run01_*.tar.gz > guide09_run01_SHA256.txt

No file is removed: the archive is a fingerprinted copy, to be kept together with the completed measurement sheet.