OpenFOAM: turbolenza, modelli RANS k-ε e k-ω SST e trattamento delle pareti

Con i dati del tutorial (U = 10 m/s, ν = 1·10⁻⁵ m²/s) e un’altezza del gradino e del canale d’ingresso di 25,4 mm (da verificare nel dizionario della mesh, sezione 3):

L’ultimo numero è il più istruttivo: già all’ingresso la viscosità turbolenta è circa 85 volte quella molecolare. Nel calcolo RANS il comportamento del flusso dipende quindi in modo determinante da come il modello produce e dissipa ν_t: è il motivo per cui la scelta del modello sposta il riattacco in modo misurabile.

La decomposizione di Reynolds scrive ogni grandezza come media più fluttuazione; mediando le equazioni compare il tensore degli sforzi di Reynolds, che va modellato. Con l’ipotesi di Boussinesq:

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

Tutta la fisica delle fluttuazioni è compressa in un solo scalare, ν_t. Le conseguenze pratiche: gli sforzi normali risultano isotropi a meno del gradiente medio, gli effetti di curvatura e rotazione non sono rappresentati in modo naturale, e la zona di separazione dietro il gradino, dove la turbolenza è fortemente anisotropa, è proprio il punto in cui i modelli a due equazioni differiscono di più.

Prima di preparare i casi, scrivi in un file di testo dentro la cartella della campagna le grandezze che vuoi confrontare e come le misurerai: per esempio «lunghezza di riattacco sulla parete inferiore, definita dal cambio di segno della componente x dello sforzo di parete» e «ΔP fra le sezioni x = −0,02 m e x = 0,2 m, media pesata sull’area». Decidere dopo, guardando i risultati, porta inevitabilmente a scegliere la definizione che conferma l’idea di partenza.

Per una DNS il rapporto fra la scala più grande e la scala di Kolmogorov cresce come Re^(3/4) in ogni direzione, quindi il numero di punti in 3D cresce come Re^(9/4). Con Re ≈ 2,5·10⁴:

N_DNS ~ Re^(9/4) = 25400^2.25 ≈ 8·10⁹ punti (ordine di grandezza)
N_RANS 2D (pitzDaily) ≈ 1·10⁴ celle
rapporto ≈ 10⁶

A questo si aggiunge che la DNS è necessariamente tridimensionale e transitoria, con migliaia di passi temporali per ottenere statistiche convergenti, mentre il calcolo RANS stazionario converge in qualche migliaio di iterazioni su una mesh 2D. La LES si colloca in mezzo: risolve i vortici grandi, ma vicino alle pareti il costo torna a crescere rapidamente con Re se lo strato limite viene risolto. È questo divario di sei o più ordini di grandezza che rende il RANS lo strumento quotidiano dell’ingegneria, con i limiti che la guida mette in evidenza.

• RANS stazionario: la mesh deve risolvere i gradienti del campo medio; raffinarla migliora la soluzione delle equazioni RANS, non avvicina il risultato alla turbolenza reale oltre quanto il modello consente;
• URANS: rappresenta solo le instabilità di grande scala separate nettamente dallo spettro turbolento (per esempio un distacco periodico); le fluttuazioni turbolente restano modellate da ν_t;
• LES: la dimensione del filtro è legata alla mesh, quindi la mesh fa parte del modello; uno studio di convergenza spaziale ha un significato diverso rispetto al RANS e richiede statistiche temporali.

Il caso prende il nome dall’esperimento di Pitz e Daily (1983) su uno strato di miscelazione turbolento a valle di un gradino rivolto all’indietro. Le dimensioni si leggono direttamente dal dizionario della risorsa:

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

Con le coordinate dei vertici si ricavano altezza del canale d’ingresso, altezza del gradino e rapporto di espansione, che vanno annotati nel quaderno del progetto prima di qualunque confronto. Nella letteratura sui gradini con rapporto di espansione intorno a 2 e Re dell’ordine di 10⁴ la lunghezza di riattacco misurata è tipicamente di alcune altezze di gradino (dell’ordine di 6–8 h), ed è noto che il k-ε standard tende a sottostimarla, perché produce troppa viscosità turbolenta nello strato di taglio. Usa questi valori come ordine di grandezza per riconoscere un risultato anomalo, non come dato di validazione: per un confronto vero servono le misure originali e le stesse condizioni d’ingresso.

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

L’Allrun ufficiale è la fonte di verità sui comandi: se una nuova versione cambia la sequenza (per esempio il percorso della risorsa della mesh), lo script della guida va adattato a quello, non viceversa. Il blocco SIMPLE mostra i criteri di arresto (residualControl) con cui il solver dichiara la convergenza.

Il limitatore di ν_t è la differenza che conta di più in questo caso: nello strato di taglio che parte dallo spigolo del gradino k-ε produce viscosità turbolenta elevata, che rimescola rapidamente e accorcia la bolla di ricircolo; SST la limita e tende a produrre bolle più lunghe. I coefficienti effettivamente usati vengono stampati all’inizio di log.foamRun (blocco kEpsilonCoeffs o kOmegaSSTCoeffs): conservali come parte della documentazione del calcolo.

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

Prima del lancio il diff deve mostrare soltanto la riga model e il file GUIDE09_MODEL.txt. Salva l’output in un file della campagna (diff -ru ... > diff_before_run.txt): è la prova documentale che il confronto è controllato.

#!/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.

Le tre relazioni sono costruite in modo che i due modelli partano dalla stessa viscosità turbolenta: ν_t = Cμ k²/ε = 0,09 · 0,375² / 14,856 ≈ 8,52·10⁻⁴ m²/s e ν_t = k/ω = 0,375 / 440,17 ≈ 8,52·10⁻⁴ m²/s. Questo è il controllo da fare ogni volta che si modificano i valori d’ingresso: se i due ν_t non coincidono, i due modelli partono da condizioni diverse e il confronto non è più controllato.

La lunghezza ℓ = 2,54 mm corrisponde al 10 % dell’altezza del canale d’ingresso; una regola empirica diffusa è ℓ ≈ 0,07 Dₕ, che con Dₕ = 50,8 mm darebbe 3,6 mm. Nessuna delle due è «giusta»: sono ipotesi. Poiché ν_t d’ingresso cambia proporzionalmente a I e a ℓ, una campagna di sensibilità su questi due parametri (separata dal confronto dei modelli, come chiede la guida) dice quanto il riattacco dipende da informazioni che di solito non si conoscono con precisione.

Per il canale d’ingresso, a monte del gradino, si può usare la correlazione di Dean per il canale piano turbolento, Cf = 0,073 Re^(−1/4) con Re riferito all’altezza del canale. Lo script seguente, aggiunto a questa guida, solo libreria standard, applica la stima:

#!/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)

La lunghezza viscosa è di circa 19 μm: per risolvere il sottostrato viscoso (y⁺ ≈ 1) servirebbe una prima cella di circa 37 μm, per le funzioni di parete logaritmiche (y⁺ ≈ 30) circa 1,1 mm. A valle del gradino la stima non vale: nella bolla lo sforzo di parete passa per zero, quindi y⁺ tende a zero vicino ai punti di separazione e riattacco qualunque sia la mesh. È il motivo per cui la guida chiede di osservare la distribuzione, non un valore singolo.

• nutkWallFunction: nella cella di parete calcola ν_t a partire da k in modo che lo sforzo di parete segua la legge logaritmica quando y⁺ supera un valore di soglia (dell’ordine di 11, il punto in cui legge lineare e logaritmica si incontrano); sotto la soglia lascia agire la sola viscosità molecolare;
• epsilonWallFunction e omegaWallFunction: impongono ε o ω nella cella di parete secondo la scala della regione logaritmica (ω anche secondo quella viscosa, con una miscela fra le due nelle versioni recenti);
• kqRWallFunction: gradiente nullo di k alla parete, coerente con le funzioni precedenti.

I dettagli (soglie, formulazioni miscelate) cambiano fra versioni: il riferimento è il codice sorgente della propria installazione, consultabile nella documentazione C++ citata fra le fonti. La conseguenza operativa resta: y⁺ va confrontato con il campo di validità della funzione di parete effettivamente scelta.

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

Una tabella così, salvata nella cartella del caso, mostra subito un errore frequente: un campo con la condizione sbagliata o mancante su una patch. Ripetila per upperWall, inlet e 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()

• verifica che tutti i file necessari al confronto esistano nel tutorial prima di creare qualsiasi cosa;
• esige esattamente una riga model kEpsilon;: se il tutorial di una versione futura cambiasse, lo script si ferma invece di produrre un confronto sbagliato;
• in caso di errore rimuove soltanto la cartella appena creata da sé stesso, mai lavori preesistenti, perché la destinazione esistente è rifiutata all’inizio;
• scrive cases.csv con lo stato not_run: un registro minimo della campagna da aggiornare a mano.

#!/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"

Il controllo su constant/polyMesh è una protezione importante: rilanciare blockMesh su un caso già calcolato sostituirebbe la mesh lasciando campi salvati incompatibili. Il post-processing di y⁺ è separato dal calcolo con if !: un suo errore produce un avviso ma non cancella né invalida la soluzione del flusso.

• nel calcolo stazionario Time = ... è un contatore di iterazioni, non un tempo fisico;
• se residualControl è attivo, il solver si ferma prima di endTime con un messaggio del tipo «SIMPLE solution converged in N iterations»: annota N per ciascun modello;
• righe bounding k, bounding epsilon o bounding omega indicano valori negativi riportati a un minimo: poche nelle prime iterazioni sono comuni, persistenti fino alla fine segnalano un problema di mesh, condizioni o schemi;
• errori di continuità sum local che non diminuiscono indicano una soluzione non convergente anche se i residui sembrano bassi.

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

I due casi sono indipendenti: con almeno due core liberi si possono lanciare in due finestre tmux distinte (Ctrl+b poi c per una nuova finestra). Su una mesh 2D di circa diecimila celle ciascun calcolo usa un solo core e poca memoria; il tempo totale si dimezza senza alcun rischio per i dati, perché le cartelle sono separate. Controlla prima il carico con uptime o htop se il server è condiviso.

La function object yPlus scrive il campo yPlus nella cartella dell’ultimo istante salvato e stampa nel log un riepilogo per patch di parete con minimo, massimo e media. Il riepilogo è la prima cosa da guardare:

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

Una media accettabile con un minimo vicino a zero e un massimo molto alto non è contraddittoria: nella bolla di ricircolo lo sforzo di parete è piccolo, sotto lo strato di taglio d’ingresso è grande. Il giudizio va fatto per zone: canale d’ingresso, regione dietro il gradino, zona di riattacco, parete superiore.

Lo stesso meccanismo di post-processing calcola lo sforzo di parete, che serve per definire il riattacco in modo riproducibile (sezione 11). Prima di usarlo verifica il nome della funzione disponibile nella tua installazione:

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

Il campo ottenuto è per unità di densità [m²/s²] e ha tre componenti: per il riattacco interessa la componente parallela alla parete inferiore, cioè x in questa geometria. Anche questo comando scrive un log separato e non modifica i campi del flusso.

• in ParaView applica Extract Block e seleziona la sola patch lowerWall (con wallShearStress fra i campi letti);
• se la geometria è 2D con uno strato di celle, la patch è una striscia: Plot Data o Plot Over Line lungo x a z costante forniscono la componente x in funzione della posizione;
• esporta in CSV con File → Save Data e calcola i cambi di segno con lo script seguente, aggiunto a questa guida.

#!/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()

Provato su un profilo sintetico con cambio di segno noto a x = 0,17 m, lo script restituisce:

$ 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.

Più di un cambio di segno è normale: vicino allo spigolo può esistere una piccola bolla secondaria controrotante. Riporta tutti i punti trovati e indica quale consideri riattacco principale, con il motivo.

Oltre alla lunghezza, la bolla ha un’area misurabile. In ParaView: Threshold su U_X con limite superiore 0 (celle con flusso all’indietro), poi Integrate Variables: il volume integrato diviso per lo spessore del dominio in z dà l’area della bolla nel piano. Fatto con lo stesso filtro sui due casi, è un indicatore in più, poco sensibile a un singolo punto di campionamento. Dichiara sempre il criterio (U_x < 0, non «zona blu della mappa colori»).

In unità cinematiche la pressione dinamica di riferimento è ½U² = 50 m²/s² qualunque sia il fluido; il valore in pascal dipende interamente dalla densità dichiarata. Nota che ν = 1·10⁻⁵ m²/s non corrisponde esattamente né all’aria (circa 1,5·10⁻⁵) né all’acqua (circa 1·10⁻⁶): è un valore di tutorial. Per questo la densità va dichiarata come ipotesi della misura e un ΔP in pascal ha senso solo accompagnato da essa.

Per l’ingresso e l’uscita, che sono patch, le medie si possono calcolare senza ParaView con le configurazioni predefinite di foamPostProcess. I nomi e gli argomenti dipendono dalla versione: elencali e leggi l’intestazione del file di configurazione prima di usarli.

foamPostProcess -list | grep -i -E "patch|surface"
# esempio da adattare a quanto elencato
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

Il risultato finisce in postProcessing/ con un file per funzione: la differenza delle due medie, moltiplicata per la densità dichiarata, è il ΔP fra ingresso e uscita. Per sezioni interne (per esempio a x = 0,2 m) servono invece un piano di campionamento o Slice + Integrate Variables in ParaView, con la media pesata sull’area.

#!/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,,,,,,,,,,,,

Lo script legge l’ultimo residuo iniziale per campo: è un’istantanea, non una storia. Nota anche che rifiuta di sovrascrivere un CSV esistente, coerentemente con il resto del pacchetto. La scheda di misura separa colonne di valori e colonna measurement_method: compilarla senza il metodo rende i numeri inutilizzabili per un confronto futuro.

• i residui iniziali sono scesi di alcuni ordini di grandezza e si sono stabilizzati, non solo «sotto la soglia»;
• le grandezze di interesse (ΔP, riattacco, portata in uscita) non cambiano più, entro la precisione richiesta, nelle ultime centinaia di iterazioni: per verificarlo servono function object attive durante il calcolo o più istanti salvati;
• la portata in uscita coincide con quella in ingresso (bilancio di massa, Guida 2);
• non restano righe bounding persistenti.

Se un modello non soddisfa questi criteri e l’altro sì, il confronto fra i due va sospeso: si stanno confrontando una soluzione e un calcolo non terminato.

Per seguire ΔP e portata durante le iterazioni, aggiungi le function object al controlDict in una nuova campagna (non in quella già preparata). Struttura indicativa, da adattare ai nomi elencati da foamPostProcess -list:

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

Le serie in postProcessing/ rendono verificabile il criterio «la grandezza di interesse non cambia più», che i soli residui non garantiscono.

Se l’esperienza B mostra y⁺ fuori campo, la correzione corretta è una nuova famiglia di casi con mesh modificata, in cui il rapporto fra spessore della prima cella e funzione di parete è scelto in modo coerente: per esempio funzioni di parete con y⁺ ≈ 30–100 in tutto il canale d’ingresso, oppure prima cella con y⁺ ≈ 1 e formulazioni adatte alla bassa y⁺. Mescolare le due strategie in mesh diverse dello stesso studio di convergenza rende l’ordine osservato (Guida 4) privo di significato.

L’indicazione «SST tipicamente più lunga» deriva dal comportamento noto dei due modelli nei flussi separati e serve come controllo di plausibilità, non come risultato da ottenere a ogni costo: se emerge il contrario, la cosa interessante è capirne il motivo.

Caso: pitzDailySteady (OpenFOAM Foundation 14), mesh resources/blockMesh/pitzDaily, N celle = ...
Modelli: kEpsilon | kOmegaSST — unica differenza verificata con diff -ru
Ingresso: U = 10 m/s, I = 0,05, l = 2,54 mm, nu_t/nu = 85
Convergenza: iterazioni ...; residui finali ...; bilancio di massa ...
y+ lowerWall (min/medio/max): ... | ...   zone critiche: ...
Riattacco (cambio di segno di tau_x su lowerWall): x_r = ... m | ... m  (x_r/h = ... | ...)
Delta P (patchAverage inlet-outlet, rho = ... kg/m3): ... Pa | ... Pa
Osservazioni / limiti / ipotesi: ...

• S. B. Pope, Turbulent Flows, Cambridge University Press, 2000: decomposizione di Reynolds, scale della turbolenza, modelli a viscosità turbolenta;
• D. C. Wilcox, Turbulence Modeling for CFD, DCW Industries: modelli k-ε e k-ω, funzioni di parete, casi di verifica;
• B. E. Launder, D. B. Spalding, The numerical computation of turbulent flows, Computer Methods in Applied Mechanics and Engineering 3 (1974): il modello k-ε standard e le funzioni di parete;
• F. R. Menter, Two-equation eddy-viscosity turbulence models for engineering applications, AIAA Journal 32 (1994): il modello SST;
• R. W. Pitz, J. W. Daily, Combustion in a turbulent mixing layer formed at a rearward-facing step, AIAA Journal 21 (1983): l’esperimento da cui prende nome il caso.

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

Nessun file viene rimosso: l’archivio è una copia con impronta, da conservare insieme alla scheda di misura compilata.