OpenFOAM: condizioni al contorno e modelli fisici in un caso completo

Per un fluido newtoniano a densità costante, in forma cinematica (divisa per ρ):

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

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

Due conseguenze pratiche. Primo: il solver non conosce ρ, quindi né la pressione né le forze escono in unità SI finché non le si moltiplica per una densità dichiarata. Secondo: nel caso incomprimibile la pressione non è una variabile termodinamica ma il moltiplicatore che impone ∇·U = 0; per questo è definita a meno di una costante, ed è il motivo per cui serve almeno un valore di riferimento (qui p = 0 all’uscita). Il termine div((nuEff*dev2(T(grad(U))))) presente in fvSchemes corrisponde alla parte trasposta dello sforzo viscoso.

Per le equazioni di Navier–Stokes incomprimibili, su ogni porzione di contorno si prescrive in genere o la velocità (ingressi, pareti) o una condizione «di trazione», che in pratica OpenFOAM approssima imponendo la pressione e un gradiente nullo di velocità (uscite). La combinazione corretta per un condotto è quindi: velocità o portata imposta all’ingresso con gradiente nullo di pressione; pressione imposta all’uscita con gradiente nullo di velocità; aderenza sulle pareti. Imporre sia velocità sia pressione sulla stessa patch sovradetermina il problema; non imporre la pressione da nessuna parte lo lascia indeterminato.

Nei condotti la transizione al moto turbolento avviene, in condizioni ordinarie, intorno a Re ≈ 2000–2300 basato sul diametro idraulico; con ingressi molto disturbati può anticipare, con ingressi molto curati può ritardare di molto. Con Re ≈ 133 il caso è laminare con ampio margine. Il secondo numero da guardare è il rapporto fra lunghezza del canale e lunghezza d’ingresso, L/L_e ≈ 0,20/0,044 ≈ 4,5: il moto si sviluppa e resta sviluppato per circa i tre quarti del canale, condizione necessaria per confrontarlo con la soluzione analitica.

Il controllo più rapido è leggere l’intestazione di un log qualsiasi e confrontare i nomi dei file con quelli del tutorial installato. Un dizionario in più non sempre produce errore: viene semplicemente ignorato, e questo è il caso più insidioso perché si crede di aver impostato qualcosa che il solver non legge.

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

La conversione non si fa a memoria: si copia il tutorial equivalente della versione 14 e vi si trasferiscono mesh, campi e valori del caso originale, controllando ogni dizionario con foamDictionary -keywords. Il log del primo lancio va letto per intero alla ricerca di avvisi su voci ignorate.

Per moto laminare completamente sviluppato in un condotto rettangolare esiste una soluzione in serie (Shah e London). Con sezione 10 mm × 5 mm (rapporto 0,5) il prodotto fra fattore d’attrito di Darcy e Reynolds vale f_D·Re ≈ 62,19 (Fanning 15,55). Da qui:

Con ingresso uniforme la Δp calcolata sarà maggiore del valore sviluppato: il tratto d’ingresso aggiunge una caduta concentrata dovuta all’accelerazione del nucleo e agli sforzi di parete più alti. Per un confronto pulito con la teoria si misura il gradiente fra due sezioni nel tratto sviluppato, per esempio x = 0,10 m e x = 0,18 m, non fra ingresso e uscita. Questo è un vero controllo di verifica, molto più informativo del solo «STATUS OK».

# gradiente sviluppato e Δp attesi (Python, solo libreria standard)
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

Le celle misurano 2,5 mm in x e 0,625 mm in y e z: rapporto d’aspetto 4, allineato con la direzione del moto, dove i gradienti sono deboli. Sul lato corto (5 mm) ci sono 8 celle, sul lato lungo (10 mm) 16: per un profilo quasi parabolico è una risoluzione sufficiente a cogliere la forma, ma un confronto quantitativo con la velocità massima (1,99·U) o con il gradiente di pressione richiede lo studio di convergenza della Guida 4, per esempio con 40×8×4, 80×16×8 e 160×32×16 celle (rapporto r = 2).

La tentazione di risparmiare celle con un modello bidimensionale qui porta a errori grossolani, perché la profondità (5 mm) è minore dell’altezza (10 mm) e le pareti laterali dominano l’attrito. Gradiente di pressione sviluppato a parità di velocità media (0,02 m/s, μ = 10⁻³ Pa·s):

Un modello 2D è lecito solo quando la dimensione trascurata è molto maggiore delle altre (rapporto d’aspetto della sezione di almeno 8–10) o quando il problema è davvero piano per simmetria. È una decisione di modellazione da giustificare con numeri, non una scorciatoia di calcolo.

Il tipo geometrico scritto in constant/polyMesh/boundary vincola le condizioni ammesse. Le patch empty, symmetryPlane, wedge e cyclic sono vincolate: nei campi il tipo deve essere lo stesso, altrimenti il solver si ferma. Il tipo wall non impone nulla ai campi di velocità e pressione, ma è necessario a funzioni che cercano le pareti: distanza dalla parete, funzioni di parete dei modelli RAS, calcolo di y⁺ e di wallShearStress. Una parete dichiarata patch funziona in laminare ma rompe silenziosamente un caso turbolento.

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

I nomi esatti e le voci richieste da ciascuna condizione si controllano nella documentazione API della propria versione o con foamInfo (Parte 8 del percorso).

foamInfo flowRateInletVelocity
foamInfo pressureInletOutletVelocity
foamInfo totalPressure

Per un ingresso che varia nel tempo (avvio graduale, pulsazione) si usa una condizione che legge una funzione del tempo, per esempio una tabella interpolata; per condizioni di tipo Robin, che combinano valore e gradiente, esiste mixed. Esempio di rampa d’ingresso in 10 s, da usare solo in un caso transitorio:

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

La sintassi delle funzioni del tempo (table, tableFile, sine, …) si controlla con foamInfo uniformFixedValue. In un caso stazionario una rampa non ha senso fisico: il «tempo» conta le iterazioni.

Se nessuna patch impone la pressione (per esempio un dominio chiuso, o un condotto con ingresso e uscita entrambi a velocità imposta), il livello di pressione è indeterminato e va fissato in una cella. Nel dizionario dell’algoritmo si aggiungono pRefCell (o pRefPoint) e pRefValue. In un dominio chiuso la portata netta imposta deve essere esattamente nulla, altrimenti l’equazione della pressione non ha soluzione.

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;
    }
}

Le dimensioni seguono l’ordine SI di OpenFOAM [massa lunghezza tempo temperatura quantità corrente intensità]: [0 1 -1 0 0 0 0] sono m/s, [0 2 -2 0 0 0 0] m²/s². OpenFOAM verifica la coerenza dimensionale di ogni operazione: una pressione scritta in Pa [1 -1 -2 0 0 0 0] in un solver cinematico produce un errore di dimensioni, non un risultato sbagliato silenzioso.

Partire da velocità nulla è lecito ma lento: le prime iterazioni devono «inventare» un campo che rispetti la continuità. Due alternative pratiche. La prima è inizializzare il campo interno con la velocità media, internalField uniform (0.02 0 0);, che rispetta già la portata. La seconda è risolvere un flusso potenziale che soddisfa ∇·U = 0 con le stesse condizioni al contorno e usarlo come campo iniziale:

potentialFoam -help
potentialFoam -writep > logs/log.potentialFoam 2>&1   # scrive U (e p) in 0/

potentialFoam sovrascrive i campi in 0/: va eseguito sulla copia di lavoro, mai sul template, e il fatto va annotato. Nei casi stazionari la soluzione finale non deve dipendere dall’inizializzazione; se dipende, il problema ha più soluzioni o non è convergente.

L’ingresso uniforme del template è una scelta di comodo: genera un tratto di sviluppo di circa 44 mm in cui gradiente di pressione e sforzo di parete sono maggiori del valore sviluppato. Se l’obiettivo è il comportamento del canale in una rete, si preferisce imporre il profilo sviluppato oppure allungare il dominio a monte e misurare solo a valle. Un’opzione pratica, nelle versioni che la offrono, è il parametro profile di flowRateInletVelocity (profili di strato limite laminare o turbolento); in ogni caso la scelta va dichiarata nella scheda del caso.

/* 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 ();

Le parti 5 e 6 del percorso spiegano in dettaglio schemi e solutori. Per sapere quali voci la propria installazione usa nei casi stazionari, confrontare con il tutorial pitzDailySteady: è il riferimento più sicuro per le parole chiave della versione 14.

• constant/momentumTransport: simulationType RAS; e un sottodizionario RAS { model kOmegaSST; turbulence on; printCoeffs on; }.
• Nuovi campi 0/k, 0/omega, 0/nut, con condizioni d’ingresso coerenti e funzioni di parete (kqRWallFunction, omegaWallFunction, nutkWallFunction o equivalenti a bassi Reynolds).
• Schemi per div(phi,k) e div(phi,omega), solutori e rilassamento per k e omega.
• Una mesh di parete coerente con il trattamento scelto (Guida 6, y⁺).

Valori d’ingresso tipici: k = 1,5·(I·U)² con intensità I di qualche percento e ω = √k / (C_μ^¼·ℓ) con C_μ = 0,09 e scala ℓ legata al diametro idraulico (per esempio 0,07·Dₕ). A Re ≈ 133 nessuno di questi è giustificato: un modello di turbolenza in un flusso laminare aggiunge viscosità fittizia e falsa la Δp. Il punto di partenza corretto, quando servirà, è il tutorial RAS della propria installazione.

Con fluidi come sangue, polimeri o fanghi la viscosità dipende dalla velocità di deformazione. In Foundation 14 il modello si sceglie in constant/physicalProperties con viscosityModel; i nomi disponibili e i coefficienti richiesti si ricavano dai tutorial installati, non da esempi di altre versioni:

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

Con viscosità variabile la soluzione analitica di questa guida non vale più, il profilo si appiattisce (fluidi pseudoplastici) e il numero di Reynolds va definito con una viscosità di riferimento dichiarata. È un passo da fare dopo aver verificato il caso newtoniano.

Con i fattori classici p 0,3 e U 0,7 l’algoritmo SIMPLE converge in modo robusto ma lento. La variante «consistente» (SIMPLEC), dove disponibile, permette fattori più alti (per esempio U 0,9 e pressione senza rilassamento) e riduce sensibilmente il numero di iterazioni. Il criterio per cambiare è empirico ma misurabile: confrontare il numero di iterazioni necessario per portare i residui sotto soglia e verificare che la soluzione finale (Δp, Q) non cambi. Le parole chiave della versione 14 si controllano 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 crea una directory con suffisso casuale: due lanci nello stesso secondo non collidono mai.
• La trap onexit EXIT scrive STATUS.txt in ogni caso, anche quando set -e interrompe lo script: il file dice in quale fase e con quale codice si è fermato.
• Le varianti si ottengono sostituendo un solo file nella copia (0/U o constant/physicalProperties): il template non viene mai toccato.
• Codici di uscita: 2 per modalità sconosciuta, 127 se manca un eseguibile OpenFOAM; altrimenti quello dell’utility fallita, riportato 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.')

Si esegue anche su Windows, prima del trasferimento: verifica classe e nome di ogni dizionario, bilanciamento di parentesi, presenza delle tre patch nei campi, dimensioni di U e p e il solver selezionato, e stampa i valori di riferimento (area, portata, Reynolds). Non sostituisce l’esecuzione: un dizionario sintatticamente corretto può essere fisicamente sbagliato.

python scripts/preflight.py

Le tre esperienze sono indipendenti e brevi: si possono lanciare in sequenza in una sessione tmux, registrando un riepilogo finale. Ogni variante crea la propria directory, quindi un fallimento non ferma le altre se si evita set -e nel ciclo:

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

Con flowRateInletVelocity e profilo uniforme, la velocità imposta è Q/A = 1·10⁻⁶ / 5·10⁻⁵ = 0,02 m/s: identica all’esperienza A. Il confronto A–B è quindi un test di coerenza: le due soluzioni devono coincidere entro la tolleranza dei solutori. Una differenza sensibile indica un errore di area, di segno o di unità.

Raddoppiando ν a velocità fissa, nel tratto sviluppato il gradiente di pressione raddoppia (circa 14 → 28 Pa/m), perché f_D = 62,19/Re e Re si dimezza; la lunghezza d’ingresso si dimezza (44 → 22 mm); il profilo sviluppato e quindi la velocità massima restano gli stessi, perché in laminare la forma del profilo non dipende da ν. La Δp totale fra ingresso e uscita raddoppia meno che proporzionalmente, perché il contributo del tratto d’ingresso cresce meno della parte viscosa.

Se il calcolo mostra un rapporto 2,00 fra i gradienti misurati nel tratto sviluppato, il caso è verificato; se mostra un rapporto 2 sulla Δp totale, probabilmente si stanno confrontando grandezze diverse o la soluzione non è convergente.

• D — Posizione dell’uscita. Allungare il canale a 0,30 m (vertici e numero di celle in x in proporzione, 120 × 16 × 8) mantenendo tutto il resto: il gradiente misurato fra 0,10 e 0,18 m non deve cambiare. Se cambia, l’uscita influenzava la soluzione.
• E — Raffinamento. Ripetere A con 40 × 8 × 4 e 160 × 32 × 16 celle (r = 2) e applicare la procedura della Guida 4 al gradiente di pressione sviluppato. Il valore estrapolato va confrontato con 14,0 Pa/m: è una verifica completa, numerica e analitica insieme. La mesh più fine ha 81 920 celle, ancora leggera.

Il caso esegue 700 iterazioni senza un criterio d’arresto. In un uso serio si aggiunge un controllo sui residui, così che il calcolo si fermi quando i residui iniziali di tutte le equazioni scendono sotto soglie dichiarate. Nella versione 14 la voce si inserisce nel dizionario dell’algoritmo; la sintassi esatta si verifica nel tutorial pitzDailySteady:

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

Il residuo da solo non basta: si affianca la storia delle grandezze di interesse. Qui la scelta naturale è aggiungere misure di pressione media su due sezioni del tratto sviluppato e verificare che la loro differenza resti costante nelle ultime centinaia di iterazioni. Un esempio di function object da aggiungere a functions nella copia di lavoro:

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; }
    }
}

Le parole chiave delle superfici di campionamento sono cambiate fra le versioni: prima di usarlo confrontare con foamInfo surfaces o con i file in $FOAM_ETC/caseDicts/functions. In alternativa le medie di sezione si calcolano in ParaView (Guida 3).

A ogni iterazione il log riporta una riga «time step continuity errors» con tre numeri: errore locale medio (sum local), globale (global) e cumulativo. Sono errori sulla divergenza del flusso dopo la correzione di pressione, normalizzati: devono diventare molto piccoli e il «global» deve oscillare intorno a zero senza deriva. Una deriva del cumulativo indica una perdita sistematica di massa, tipicamente da tolleranze troppo larghe della pressione.

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 da (0,15; 0; 0,0025) a (0,15; 0,01; 0,0025) campiona il profilo lungo l’altezza a metà profondità, nel tratto sviluppato. Lo script seguente (solo libreria standard) calcola la soluzione in serie negli stessi 21 punti; l’RMSE fra i due profili, normalizzato su U, misura l’errore di discretizzazione locale (Guida 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}")

Il profilo lungo il lato lungo è più piatto al centro di una parabola, perché a metà profondità il gradiente è dominato dalle pareti vicine (lato corto). Una parabola di Poiseuille fra lastre piane darebbe un massimo di 1,5·U: un risultato vicino a 1,5 indicherebbe un errore di geometria o un caso inavvertitamente bidimensionale.

Un errore fatale riporta il messaggio, il file e spesso la riga del dizionario coinvolto, e la funzione del codice che l’ha generato. Il messaggio utile è quello che precede «From function»; le righe dello stack che seguono servono agli sviluppatori. Gli errori più frequenti in questa fase: «Cannot find patchField entry for walls» (nome patch diverso), «keyword nu is undefined» (dizionario di un’altra release), «Different dimensions for (p + …)» (unità sbagliate), «Unknown patchField type» (condizione scritta male o non disponibile nella versione).

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

• Leggere in quale equazione esplode il residuo per prima (Ux o p) e a quale iterazione.
• Rieseguire checkMesh e confrontare con la Guida 6: una mesh appena modificata è il primo sospetto.
• Ridurre i fattori di rilassamento (per esempio p 0,2, U 0,5) e passare temporaneamente a bounded Gauss upwind per div(phi,U): se così converge, il problema è numerico; se diverge comunque, è nelle condizioni al contorno.
• Ripristinare le impostazioni originali partendo dalla soluzione convergente (startFrom latestTime): lo schema del secondo ordine deve essere quello del risultato finale.

Lo stesso schema della Guida 5 vale qui: un archivio per variante, senza i campi intermedi che si possono rigenerare, con le impronte dei file e la descrizione dell’ambiente.

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())

• Per ogni istante tiene un solo valore; se lo stesso tempo compare in più file (ripartenze), prevale il file modificato per ultimo.
• Usa l’ultimo istante comune alle due misure, quindi non confronta mai tempi diversi.
• Lo squilibrio è relativo: |Qin + Qout| / max(|Qin|, |Qout|); se un dato manca la cella resta vuota.

Con il caso nominale ci si aspetta Q_in ≈ −1·10⁻⁶ e Q_out ≈ +1·10⁻⁶ m³/s; lo squilibrio relativo di un caso incomprimibile convergente è dell’ordine della tolleranza della pressione, molto sotto 10⁻⁶.

Nel tratto sviluppato il bilancio di forze su un volume di lunghezza L impone Δp·A = τ̄_w·P·L, con P perimetro bagnato (30 mm). È un secondo controllo indipendente: con 14 Pa/m e A/P = Dₕ/4 = 1,667 mm lo sforzo medio atteso è τ̄_w ≈ 0,023 Pa. In OpenFOAM lo sforzo a parete si calcola con il function object wallShearStress (cinematico, da moltiplicare per ρ) e si media sulla patch in ParaView o con surfaceFieldValue.

foamPostProcess -help
foamPostProcess -func wallShearStress -latestTime   # sulla copia di lavoro, dopo il calcolo

A e B impongono lo stesso campo d’ingresso, quindi la differenza fra le loro soluzioni misura solo il rumore iterativo. Un criterio ingegneristico ragionevole: la differenza sui valori di interesse (Δp fra le due sezioni, Q) deve essere almeno un ordine di grandezza più piccola di qualsiasi effetto che si vuole attribuire alla fisica, per esempio la differenza A–C. Se A e B differiscono per più dell’1 %, la convergenza iterativa è insufficiente e il confronto A–C non è ancora interpretabile.

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,,,,,,,,

Valori attesi per la verifica: dpdx_developed ≈ 14,0 / 14,0 / 28,0 Pa/m; umax ≈ 0,0398 m/s in tutti e tre i casi. Le colonne vuote si riempiono con i dati di collect_fluxes.py e di ParaView; le unità sono nel nome della colonna, così il file resta leggibile anche fuori contesto.

• R. K. Shah, A. L. London, Laminar Flow Forced Convection in Ducts, Academic Press, 1978. Soluzioni e correlazioni per condotti di sezione qualsiasi, incluso f·Re dei rettangoli.
• F. M. White, Viscous Fluid Flow, McGraw-Hill. Moto laminare nei condotti, lunghezza d’ingresso, strato limite.
• S. B. Pope, Turbulent Flows, Cambridge University Press, 2000. Fondamenti dei modelli RAS e LES citati nella sezione 6.
• Tutorial della propria installazione: incompressibleFluid/pitzDailySteady, incompressibleFluid/cylinder, da confrontare voce per voce con il caso della guida.