OpenFOAM: trasferimento di calore, convezione naturale e numero di Nusselt

Con i dati del tutorial (aria come gas perfetto, M = 28,96 g/mol, p = 10⁵ Pa, c_p = 1004,4 J/(kg K), μ = 1,831·10⁻⁵ Pa s, Pr = 0,705) e le proprietà valutate alla temperatura media fra le pareti, lo script cavity_estimates.py (riquadro successivo) fornisce:

Il rapporto d’aspetto H/W = 28,7 e Ra_W di ordine 10⁶ collocano il problema nel campo di validità della correlazione di Elsherbiny, Raithby e Hollands per cavità verticali (5 ≤ H/W ≤ 110, Ra_W fino a circa 2·10⁷). L’intestazione dei file sperimentali distribuiti con il tutorial riporta Ra = 8,6·10⁵: la piccola differenza rispetto al valore calcolato dipende dalle proprietà usate per l’aria. La scheda del caso nella banca dati ERCOFTAC (KBwiki, UFR 4-10) riporta per l’esperimento un Nusselt medio di 5,85 con ΔT = 19,6 °C (Ra = 0,86·10⁶) e di 7,57 con ΔT = 39,9 °C (Ra = 1,43·10⁶): la correlazione è entro il 4 % e, con Nu = 5,85, la potenza attesa è circa 44,6 W. Se la potenza calcolata da OpenFOAM si discosta di un fattore 2 da questi valori, prima di ogni interpretazione fisica si cercano errori di segno, di area o di convergenza.

Calcola proprietà al valore medio, numeri di Rayleigh, potenza di sola conduzione, Nusselt da correlazione e un limite superiore dello scambio radiativo fra due lastre parallele grigie (riquadro della sezione 6). Usa solo la libreria standard di Python e funziona su Windows e su Linux:

python cavity_estimates.py
python cavity_estimates.py --Th 317.75
python cavity_estimates.py --emissivity 0.05 0.9

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Order-of-magnitude estimates for the OpenFOAM 14 buoyantCavity case.

Computes ideal-gas properties at the mean temperature, Rayleigh numbers based on
gap width and height, the conduction-only heat rate, the Elsherbiny-Raithby-
Hollands (1982) correlation for vertical cavities and the radiative exchange
between two parallel grey plates. Nothing here is a CFD result: use the numbers
only to judge whether simulated wall powers are plausible.
"""
import argparse
import math

R_UNIVERSAL = 8.314462618  # J/(mol K)
SIGMA = 5.670374419e-8     # W/(m2 K4)


def erh_nusselt(ra, aspect):
    """Elsherbiny, Raithby & Hollands (1982), vertical cavity, isothermal walls.

    Stated range: 5 <= H/W <= 110 and Ra_W up to about 2e7. Returns Nu_W.
    """
    nu1 = 0.0605 * ra ** (1.0 / 3.0)
    nu2 = (1.0 + (0.104 * ra ** 0.293 / (1.0 + (6310.0 / ra) ** 1.36)) ** 3) ** (1.0 / 3.0)
    nu3 = 0.242 * (ra / aspect) ** 0.272
    return max(nu1, nu2, nu3), (nu1, nu2, nu3)


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument('--Th', type=float, default=307.75, help='hot wall temperature [K]')
    p.add_argument('--Tc', type=float, default=288.15, help='cold wall temperature [K]')
    p.add_argument('--W', type=float, default=0.076, help='gap width [m]')
    p.add_argument('--H', type=float, default=2.18, help='cavity height [m]')
    p.add_argument('--D', type=float, default=0.52, help='cavity depth [m]')
    p.add_argument('--p', type=float, default=1e5, help='pressure [Pa]')
    p.add_argument('--M', type=float, default=28.96, help='molar mass [g/mol]')
    p.add_argument('--Cp', type=float, default=1004.4, help='specific heat [J/(kg K)]')
    p.add_argument('--mu', type=float, default=1.831e-5, help='dynamic viscosity [Pa s]')
    p.add_argument('--Pr', type=float, default=0.705, help='Prandtl number [-]')
    p.add_argument('--g', type=float, default=9.81, help='gravity [m/s2]')
    p.add_argument('--emissivity', type=float, nargs='*', default=[0.05, 0.1, 0.3, 0.9],
                   help='wall emissivities for the radiative bound (both walls equal)')
    a = p.parse_args()
    if not a.Th > a.Tc > 0:
        p.error('Expected Th > Tc > 0 K')

    dT = a.Th - a.Tc
    Tm = 0.5 * (a.Th + a.Tc)
    Rs = R_UNIVERSAL / (a.M * 1e-3)
    rho_m, rho_h, rho_c = (a.p / (Rs * T) for T in (Tm, a.Th, a.Tc))
    k = a.mu * a.Cp / a.Pr
    nu = a.mu / rho_m
    alpha = k / (rho_m * a.Cp)
    beta = 1.0 / Tm
    ra_w = a.g * beta * dT * a.W ** 3 / (nu * alpha)
    ra_h = ra_w * (a.H / a.W) ** 3
    aspect = a.H / a.W
    area = a.H * a.D
    q_cond = k * dT * area / a.W
    nu_erh, parts = erh_nusselt(ra_w, aspect)
    u0 = math.sqrt(a.g * beta * dT * a.W)

    print('Mean temperature            Tm   = %.2f K' % Tm)
    print('Density hot / mean / cold         = %.4f / %.4f / %.4f kg/m3' % (rho_h, rho_m, rho_c))
    print('Relative density difference       = %.2f %%' % (100 * (rho_c - rho_h) / rho_m))
    print('Conductivity k = mu Cp / Pr       = %.5f W/(m K)' % k)
    print('Kinematic viscosity nu (at Tm)    = %.4e m2/s' % nu)
    print('Thermal diffusivity alpha (at Tm) = %.4e m2/s' % alpha)
    print('Expansion coefficient beta = 1/Tm = %.4e 1/K' % beta)
    print('Aspect ratio H/W                  = %.1f' % aspect)
    print('Ra based on W                     = %.3e' % ra_w)
    print('Ra based on H                     = %.3e' % ra_h)
    print('Buoyancy velocity sqrt(g beta dT W) = %.3f m/s' % u0)
    print('Wall area H x D                   = %.4f m2' % area)
    print('Conduction-only power (Nu_W = 1)  = %.2f W' % q_cond)
    print('ERH 1982: Nu1, Nu2, Nu3           = %.2f, %.2f, %.2f' % parts)
    print('ERH 1982: Nu_W = max(...)         = %.2f  ->  Q ~ %.1f W' % (nu_erh, nu_erh * q_cond))
    if not (5 <= aspect <= 110 and ra_w <= 2e7):
        print('WARNING: outside the stated range of the ERH correlation')
    q_bb = SIGMA * (a.Th ** 4 - a.Tc ** 4)
    for eps in a.emissivity:
        if not 0 < eps <= 1:
            p.error('emissivity must be in (0, 1]')
        q_rad = q_bb / (2.0 / eps - 1.0)
        print('Radiation, eps = %.2f: q = %.1f W/m2, Q = %.1f W (infinite parallel plates)' % (eps, q_rad, q_rad * area))


if __name__ == '__main__':
    main()

Le formule della correlazione sono riportate nella funzione erh_nusselt: Nu₁ = 0,0605 Ra^1/3, Nu₂ = [1 + (0,104 Ra^0,293/(1 + (6310/Ra)^1,36))³]^1/3, Nu₃ = 0,242 (Ra/A)^0,272, Nu = max(Nu₁, Nu₂, Nu₃). Sono correlazioni sperimentali con un’incertezza tipica dell’ordine del 10 %: servono per giudicare la plausibilità, non per validare.

Scrivi in un file di testo nella cartella dello studio le grandezze che confronterai e la loro definizione esatta: «Q_hot = valore Q della patch hot nel file di wallHeatFlux all’ultima iterazione, positivo se entra nel fluido», «Nu_W con L = 0,076 m, k = 0,026086 W/(m K), ΔT = T_hot − T_cold», «convergenza: variazione di Q_hot inferiore allo 0,1 % nelle ultime 100 iterazioni». Fissare definizioni e criteri prima dei risultati evita di scegliere a posteriori quella che conferma l’ipotesi.

Nel modulo fluid (file applications/modules/fluid/thermophysicalPredictor.C del sorgente della versione 14) l’incognita è l’entalpia sensibile h, non la temperatura. L’equazione assemblata è, in forma simbolica:

ddt(rho, h) + div(phi, h) + ddt(rho, K) + div(phi, K) - dp/dt + divq(h)
    = rho*(U & g) + fvModels sources

K = |U|^2/2            energia cinetica per unità di massa
phi                    flusso di massa sulle facce [kg/s]
divq(h)                flusso termico molecolare + turbolento

Con thermo hConst l’entalpia sensibile è h = c_p (T − T_std), con T_std = 298,15 K: dopo ogni soluzione T viene ricavata da h in thermo.correct(). I termini in K e il lavoro di pressione sono trascurabili qui (velocità di 0,2 m/s, variazioni di pressione di pochi pascal), ma sono presenti: per questo fvSchemes deve contenere div(phi,K). Il termine rho*(U & g) compare solo se il caso ha il campo p_rgh, cioè se è attiva la classe di galleggiamento.

La classe buoyancy di OpenFOAM 14 (src/finiteVolume/cfdTools/general/buoyancy) definisce:

p_rgh = p - rho*gh - pRef
gh    = (g & x) - ghRef,      ghRef = -|g| hRef   (hRef opzionale, predefinito 0)

Nel tutorial constant/pRef vale 10⁵ Pa e 0/p parte da 10⁵ Pa, quindi p_rgh parte da zero. Con g = (0 −9,81 0) la differenza idrostatica fra fondo e sommità della cavità è ρ g H ≈ 1,17 · 9,81 · 2,18 ≈ 25 Pa, mentre la pressione dinamica ½ρU² con U ≈ 0,15 m/s vale circa 0,013 Pa: tre ordini di grandezza di meno. Risolvere su p_rgh toglie dall’incognita la parte idrostatica dominante, che altrimenti renderebbe imprecisi gradienti piccoli. Conseguenza pratica: per leggere pressioni fisiche usare p; per confrontare zone alla stessa quota si può usare p_rgh; mai confrontare p_rgh di casi con pRef o hRef diversi.

Il tutorial usa equationOfState perfectGas: ρ = p M/(R T) con la pressione locale. L’approssimazione di Boussinesq (densità costante salvo nel termine di galleggiamento, ρ ≈ ρ₀[1 − β(T − T₀)]) è considerata accurata quando βΔT ≲ 0,1. Qui βΔT = 0,066 nel caso base e 0,098 nella variante: il caso più caldo è al limite. Con il gas perfetto non serve scegliere β: è implicito (β = 1/T per un gas perfetto). Lasciare il modello del tutorial; se si vuole studiare l’effetto, si confronta in una copia nuova con un modello di densità diverso, documentandolo come variante di modello.

Le prime righe di logs/foamRun.log riportano i modelli effettivamente costruiti. Prima di guardare i risultati verifica:

cd "$HOME/OpenFOAM/guide10_thermal_study/baseline_19p6K"
grep -n "Selecting" logs/foamRun.log | head -20

Devono comparire, fra le altre, la termodinamica heRhoThermo con perfectGas, il modello kOmegaSST e la riga sul modello di trasporto termico turbolento (riquadro della sezione 4). Una riga mancante o diversa indica che il caso non è quello che si crede di avere.

Dal system/blockMeshDict del tutorial (unità in mm): x da 0 a 76, y da 0 a 2180, z da −260 a 260; un solo blocco (35 150 15) con simpleGrading (1 1 1). La parete cold è il piano x = 0, la parete hot il piano x = 76 mm; frontAndBack sono i piani z = ±260 mm, topAndBottom y = 0 e y = 2180 mm. Tutte le patch sono di tipo wall.

I dati sperimentali distribuiti con il tutorial (validation/exptData/mt_z0_50_lo.dat) mostrano a metà altezza un salto di circa 4 K entro i primi 3 mm dalla parete fredda: lo strato limite termico ripido è risolto da una o due celle. Stimando lo sforzo di parete dalla prima velocità misurata (−0,081 m/s a 1,2 mm) si ottiene u_τ ≈ 0,03 m/s e quindi y⁺ ≈ 2 per il centro della prima cella: la parete è nel sottostrato viscoso, dove le funzioni di parete del tutorial (nutUWallFunction, omegaWallFunction) lavorano in regime a basso Reynolds. La conclusione operativa è che il flusso termico di parete, e quindi Nu, è sensibile alla risoluzione in x: lo studio di mesh della sezione 13 deve infittire soprattutto quella direzione.

Oltre ai file del caso, la cartella del tutorial contiene system/sample, che estrae T e U lungo nove linee orizzontali a metà profondità (z = 0) alle quote y/H = 0,1 … 0,9, e validation/exptData, con i profili di temperatura media (mt_z0_*_lo.dat, x in mm, T in °C) e di velocità verticale media (mv_z0_*_lo.dat, x in mm, V in m/s) misurati da Betts e Bokhari. Il file README precisa che il caso considera la minore delle due differenze di temperatura dell’esperimento, 19,6 °C. L’Allrun ufficiale esegue:

runApplication blockMesh
runApplication foamRun
runApplication -s sample foamPostProcess -latestTime -func sample
runApplication validation/createGraphs

Lo script prepare_cases.py copia anche system/sample e validation/, quindi entrambe le copie possono essere confrontate con le misure (riquadri della sezione 10). La variante a 29,6 K invece non ha dati sperimentali corrispondenti: per essa sono possibili soltanto verifiche interne (bilancio, mesh, correlazioni).

Dopo la preparazione, un confronto ricorsivo deve mostrare per il caso base nessuna differenza e per la variante una sola riga cambiata:

cd "$HOME/OpenFOAM/guide10_thermal_study"
diff -ru "$FOAM_TUTORIALS/fluid/buoyantCavity" baseline_19p6K && echo "baseline identical"
diff -ru "$FOAM_TUTORIALS/fluid/buoyantCavity" hotter_29p6K

Eseguendo lo script sui file del repository ufficiale di OpenFOAM 14 la seconda differenza è esattamente value uniform 307.75; // 34.6 degC → value uniform 317.75; // 34.6 degC. Nota il commento: lo script cambia il valore ma lascia il commento originale, che nella variante diventa sbagliato (317,75 K corrispondono a 44,6 °C). Correggi il commento a mano nella copia oppure annotalo nel registro dello studio, per non confondere chi leggerà il file.

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}

Il file constant/physicalProperties del tutorial contiene anche hf 0 (entalpia di formazione, irrilevante senza reazioni). Qualsiasi variazione di questi termini è una variazione di modello, non di parametro, e va trattata in una copia separata.

Il flusso termico turbolento è modellato con una diffusività turbolenta α_t = ρ ν_t/Pr_t (dimensioni [kg/(m s)], da cui l’indicazione «dimensioni dinamiche» della tabella). Nel tutorial 0/alphat usa su tutte le pareti compressible::alphatWallFunction con Prt 0.85. Nel volume, invece, il modello di trasporto termico si sceglie nel file constant/thermophysicalTransport, che il tutorial non contiene. Leggendo il sorgente di OpenFOAM 14 (RASThermophysicalTransportModel.C e unityLewisEddyDiffusivity.C) si vede che, in assenza del file, viene selezionato unityLewisEddyDiffusivity con Pr_t predefinito pari a 1.

Il log lo conferma con la riga «Selecting default RAS thermophysical transport model unityLewisEddyDiffusivity». Il caso usa quindi Pr_t = 0,85 alle pareti e 1 nel volume. Per renderlo esplicito, o per studiarne la sensibilità, in una copia nuova si crea:

// constant/thermophysicalTransport (FoamFile header: object thermophysicalTransport;)
RAS
{
    model           unityLewisEddyDiffusivity;
    Prt             0.85;
}

e si verifica nel log la riga «Selecting RAS thermophysical transport model unityLewisEddyDiffusivity» seguita dal valore stampato. La differenza fra le due scelte è una misura dell’incertezza di modello sul flusso termico: va riportata, non scelta in base all’accordo con un riferimento.

Due conseguenze pratiche. Primo: se il calcolo arriva a 1000 iterazioni senza soddisfare residualControl, non è «convergito», è stato fermato; la cosa si riconosce dal log (riquadro della sezione 8). Secondo: con purgeWrite 3 non resta la storia dei campi, quindi l’andamento di Q_hot durante le iterazioni si ottiene solo con una function object attiva durante il calcolo (riquadro della sezione 8), non a posteriori.

Il valore 0,02608 della guida è 0,026086 troncato: nei calcoli di Nu usa sempre lo stesso valore in tutti i casi e dichiaralo. Poiché κ = μ c_p/Pr con tre costanti, κ non dipende dalla temperatura; ν e α invece cambiano con T attraverso la densità. Per l’aria reale fra 288 e 318 K la conducibilità aumenta di circa il 7 %: il modello a proprietà costanti è una scelta del tutorial, da citare fra le ipotesi.

Ogni campo di OpenFOAM dichiara le dimensioni nella riga dimensions. Nel tutorial si leggono nomi simbolici; la corrispondenza con le unità SI è:

Se una variante richiede di scrivere temperature in °C o pressioni relative, convertile prima: OpenFOAM non converte, e una temperatura di 34,6 scritta al posto di 307,75 viene accettata senza errori come 34,6 K.

Imporre T su una parete equivale a supporre una parete di capacità termica infinita, mantenuta isoterma qualunque sia il flusso che riceve. È una buona idealizzazione quando la parete reale è un buon conduttore termostatato (piastre metalliche con circuito di fluido, come negli esperimenti su cavità) e il suo numero di Biot verso il fluido è piccolo. Il flusso termico è allora un risultato del calcolo, distribuito in modo non uniforme sull’altezza: elevato dove lo strato limite inizia (in basso sulla parete calda, in alto sulla fredda) e minore dove è già sviluppato.

Il caso di base non include l’irraggiamento. Per due lastre parallele grigie con la stessa emissività ε il flusso netto è q = σ(T_h⁴ − T_c⁴)/(2/ε − 1). Con le temperature del tutorial cavity_estimates.py dà:

Le stime sono per lastre infinite e trascurano l’assorbimento dell’aria, ma dicono una cosa chiara: con superfici non lucidate lo scambio radiativo fra le pareti è dello stesso ordine o maggiore di quello convettivo. Nella simulazione non entra nel bilancio del fluido (l’aria è trasparente), ma modifica le temperature delle pareti adiabatiche e quindi la stratificazione. Prima di confrontare Nu con un esperimento verifica come la fonte ha trattato l’irraggiamento e quale emissività avevano le superfici.

Nel sorgente di OpenFOAM 14 la condizione externalTemperature accetta, in alternativa: potenza totale Q [W], flusso q [W/m²] oppure coefficiente h [W/(m² K)] con temperatura ambiente Ta [K]; se oltre a Q o q si danno anche h e Ta, i flussi si sommano. Sono opzionali emissivity (irraggiamento verso l’ambiente) e relaxation. Ta, Q e q sono funzioni del tempo (per esempio constant 300), mentre h è un campo sulla patch: nei tutorial della versione 14 è sempre scritto come h uniform 8;. Se la riga h constant 10; dell’esempio produce un errore di lettura, usa:

wall
{
    type            externalTemperature;
    Ta              constant 300;
    h               uniform 10;
    value           uniform 310;
}

Esempi completi e funzionanti sono nei tutorial ufficiali fluid/roomHeating e multiRegion/CHT/hotWire (h = 8) e multiRegion/CHT/vacuumFurnace (h = 100): leggerli con grep -A6 externalTemperature $FOAM_TUTORIALS/fluid/roomHeating/0/T.

Le pareti frontAndBack e topAndBottom hanno zeroGradient su T: nessun flusso le attraversa. In un esperimento reale le pareti «isolate» scambiano sempre un poco di calore, conducono lungo il proprio spessore e irraggiano verso le pareti termiche. Le differenze fra calcolo e misura vicino al fondo e alla sommità della cavità (y/H = 0,1 e 0,9) dipendono spesso da queste condizioni più che dal modello di turbolenza. Nel confronto con i dati della sezione 10 è utile separare le quote centrali da quelle estreme.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Create non-destructive OpenFOAM 14 buoyantCavity thermal study cases."""
import argparse
import os
from pathlib import Path
import re
import shutil
import sys


def update_hot_temperature(text: str, new_temp: float) -> str:
    # This parser is intentionally constrained to the verified OF14 tutorial input.
    pattern = r'(\bhot\s*\{[^{}]*?\bvalue\s+uniform\s+)([0-9.+eE-]+)(\s*;)'
    matches = list(re.finditer(pattern, text, re.DOTALL))
    if len(matches) != 1 or float(matches[0].group(2)) != 307.75:
        raise ValueError('Unexpected hot patch: tutorial version differs from verified OF14 source')
    return re.sub(pattern, lambda m: m.group(1) + format(new_temp, '.2f') + m.group(3), text, count=1, flags=re.DOTALL)


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--source', type=Path, default=Path(os.environ.get('FOAM_TUTORIALS', '')) / 'fluid/buoyantCavity', help='Official OF14 buoyantCavity source directory')
    parser.add_argument('--dest', type=Path, required=True, help='New, nonexisting study directory')
    args = parser.parse_args()
    source = args.source.expanduser().resolve()
    dest = args.dest.expanduser().resolve()
    if not (source/'0/T').is_file() or not (source/'system/controlDict').is_file():
        parser.error('Official tutorial source not found: load OpenFOAM 14 shell first')
    if dest.exists():
        parser.error(f'Destination already exists; refusing to overwrite: {dest}')
    temp_data = (source/'0/T').read_text(encoding='utf-8')
    if 'uniform 288.15' not in temp_data:
        parser.error('Unexpected cold-wall temperature in source tutorial')
    # Validate before creating anything.
    baseline_text = update_hot_temperature(temp_data, 307.75)
    hot_text = update_hot_temperature(temp_data, 317.75)
    dest.mkdir(parents=True)
    for name, content in (('baseline_19p6K', baseline_text), ('hotter_29p6K', hot_text)):
        target = dest/name
        shutil.copytree(source, target, ignore=shutil.ignore_patterns('processor*', 'postProcessing', 'log.*'))
        (target/'0/T').write_text(content, encoding='utf-8')
    print(f'Created two independent cases in {dest}')
    print('Baseline: Th=307.75 K, Tc=288.15 K, dT=19.60 K')
    print('Hotter:   Th=317.75 K, Tc=288.15 K, dT=29.60 K')


if __name__ == '__main__':
    try:
        main()
    except (ValueError, OSError) as exc:
        sys.exit(f'ERROR: {exc}')

• L’espressione regolare cerca, dentro il blocco hot, la riga value uniform …;; se le occorrenze non sono esattamente una o il valore non è 307,75, si ferma: protegge da tutorial di versioni diverse.
• Le due versioni del file 0/T vengono costruite e validate prima di creare qualsiasi cartella.
• shutil.ignore_patterns('processor*', 'postProcessing', 'log.*') evita di copiare risultati di eventuali esecuzioni precedenti nella cartella sorgente.
• Se la destinazione esiste, parser.error termina con codice 2 senza toccare nulla.

Le function object di monitoraggio (sezione 8), il file di Pr_t (sezione 4) o una mesh più fitta (sezione 13) vanno inseriti in copie aggiuntive, non nelle due di riferimento. Lo script accetta qualunque destinazione nuova:

python3 scripts/prepare_cases.py \
  --source "$FOAM_TUTORIALS/fluid/buoyantCavity" \
  --dest "$HOME/OpenFOAM/guide10_thermal_study_monitor"
# then edit only baseline_19p6K inside the new folder, before running it

Annota nel registro dello studio quale copia contiene quale modifica e verificala con diff -ru rispetto al tutorial.

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Run an untouched study case exactly once; preserve all existing outputs.
set -euo pipefail
if [[ $# -ne 1 ]]; then
  echo 'Usage: run_case.sh /absolute/path/to/new/case' >&2
  exit 2
fi
CASE_DIR=$(realpath -- "$1")
[[ -f "$CASE_DIR/system/controlDict" ]] || { echo 'Invalid case' >&2; exit 2; }
command -v foamRun >/dev/null || { echo 'Load the OpenFOAM 14 environment first' >&2; exit 2; }
if [[ -e "$CASE_DIR/logs" || -e "$CASE_DIR/constant/polyMesh" ]]; then
  echo 'Case already has logs or a mesh. Refusing to rerun or overwrite it.' >&2
  exit 3
fi
mkdir "$CASE_DIR/logs"
cd "$CASE_DIR"
blockMesh > logs/blockMesh.log 2>&1
checkMesh > logs/checkMesh.log 2>&1
foamRun > logs/foamRun.log 2>&1
# Optional post-process. Its failure does not change the solver success state.
if command -v foamPostProcess >/dev/null; then
  if ! foamPostProcess -latestTime -func wallHeatFlux > logs/wallHeatFlux.log 2>&1; then
    echo 'Warning: wallHeatFlux post-processing failed; review logs/wallHeatFlux.log' >&2
  fi
fi
echo "Solver completed. Logs are in $CASE_DIR/logs"

• set -euo pipefail: lo script si ferma al primo comando fallito, su variabili non definite e su errori in una pipeline; quindi se checkMesh fallisce, foamRun non parte.
• Il controllo su logs e constant/polyMesh impedisce di rieseguire un caso già lanciato (codice di uscita 3).
• Il post-processing wallHeatFlux è dentro un if !: il suo eventuale fallimento produce solo un avviso e non annulla il successo del solver.
• Lo script non esegue sample: per il confronto con i dati sperimentali vedi la sezione 10.

cd "$HOME/OpenFOAM/guide10_thermal_study/baseline_19p6K"
grep -c "^Time = " logs/foamRun.log            # iterazioni eseguite
grep -n -i "converged" logs/foamRun.log | tail -3
grep "Solving for h" logs/foamRun.log | tail -3  # residui di h
tail -n 5 logs/foamRun.log                      # deve terminare con End

Se il numero di iterazioni è 1000 e non compare alcun messaggio di convergenza, il calcolo è stato fermato da endTime. Per estrarre l’andamento dei residui senza disturbare la cartella logs/ del caso, esegui foamLog da una cartella separata, perché crea a sua volta una sottocartella logs nella directory corrente:

mkdir "$HOME/OpenFOAM/guide10_thermal_study/foamLog_baseline"
cd "$HOME/OpenFOAM/guide10_thermal_study/foamLog_baseline"
foamLog ../baseline_19p6K/logs/foamRun.log
ls logs/   # h_0, p_rgh_0, Ux_0 ... una colonna per iterazione

OpenFOAM 14 legge automaticamente il file system/functions, se presente, oltre all’eventuale voce functions di controlDict. In una copia nuova, prima di run_case.sh, salva questo file come system/functions:

/*--------------------------------*- C++ -*----------------------------------*\
  Guide 10 - monitoring functions for a NEW copy of buoyantCavity.
  Save as <case>/system/functions before running the solver.
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      functions;
}

// Integrated wall heat flow Q [W] on the hot and cold walls every 10 iterations.
// Positive values: heat entering the fluid. No field is written (writeFields false).
wallHeatFluxMonitor
{
    type            wallHeatFlux;
    libs            ("libfieldFunctionObjects.so");
    patches         (hot cold);
    executeControl  timeStep;
    executeInterval 10;
    writeControl    timeStep;
    writeInterval   10;
    writeFields     false;
    log             yes;
}

// Initial residuals of the solved fields, one line per iteration.
#includeFunc residuals(p_rgh, U, h, k, omega)

// Maximum velocity magnitude: a cheap indicator of the circulation strength.
#includeFunc cellMaxMag(name=UmaxMag, U)

// ************************************************************************* //

La sintassi è quella dei tutorial ufficiali della versione 14 (fluid/concentricCylinders, multiRegion/CHT/wallBoiling); writeFields false evita di scrivere il campo wallHeatFlux a ogni esecuzione. Il risultato è un file postProcessing/wallHeatFluxMonitor/0/wallHeatFlux.dat con una riga per patch ogni 10 iterazioni, leggibile con whf_extract.py (sezione 9): l’andamento di Q_hot e Q_cold è il criterio di convergenza fisicamente significativo. Nel log, all’avvio, deve comparire la creazione delle tre function object; se una manca, controlla l’intestazione del file.

Dal file src/functionObjects/field/wallHeatFlux/wallHeatFlux.C di OpenFOAM 14:

• il campo vale meno il flusso termico del modello di trasporto sulla faccia di bordo, wallHeatFlux = −q, dove q è positivo nel verso della normale uscente dal dominio. Ne segue che wallHeatFlux > 0 significa calore che entra nel fluido attraverso la parete;
• se esiste un campo radiativo qr viene sottratto anch’esso, quindi il risultato include l’irraggiamento quando è modellato;
• per ogni patch vengono scritti min, max, la potenza integrata Q = Σ |S<sub>f</sub>| q [W] e la media d’area q = Q/Σ|S<sub>f</sub>| [W/m²], sia nel log sia nel file postProcessing/<nome>/<tempo>/wallHeatFlux.dat;
• senza voce patches vengono elaborate tutte le patch di tipo wall.

Quindi l’integrazione della sezione 9 è già fatta da OpenFOAM, sulla mesh effettiva e con le aree esatte delle facce: Q_hot è il valore Q della patch hot, Q_cold è il valore Q della patch cold cambiato di segno. Le patch adiabatiche, con zeroGradient, devono dare Q = 0. ParaView resta utile come verifica indipendente e per la distribuzione locale.

Nel sorgente di foamPostProcess (versione 14) i modelli fisici vengono costruiti solo se si indica il solutore con -solver; altrimenti l’utilità legge soltanto i campi. wallHeatFlux ha bisogno del modello di trasporto termico e, se non lo trova, termina con «Unable to find thermophysicalTransportModel». È proprio il caso in cui lo script run_case.sh stampa l’avviso senza compromettere il calcolo. Controlla e, se serve, ripeti con il solutore esplicito e un log con nome nuovo:

cd "$HOME/OpenFOAM/guide10_thermal_study/baseline_19p6K"
grep -i "error\|unable" logs/wallHeatFlux.log
foamPostProcess -solver fluid -latestTime -func wallHeatFlux > logs/wallHeatFlux_solver.log 2>&1
grep "for patch" logs/wallHeatFlux_solver.log

Nell’uscita cerca le righe «min, max, Q [W], q [W/m^2] for patch hot = …» e «… for patch cold = …»: il terzo numero è la potenza integrata. Lo stesso vale per yPlus e per ogni function object che richiede il modello di turbolenza o di trasporto termico; la function object sample invece legge solo T e U e funziona anche senza -solver, come nell’Allrun ufficiale.

Legge uno o più file wallHeatFlux.dat, converte i segni nella convenzione della guida, mostra la somma con segno su tutte le patch (a convergenza deve tendere a zero), la deriva di Q_hot nelle ultime registrazioni e, a richiesta, scrive la riga di misura per thermal_balance.py:

cd "$HOME/OpenFOAM/guide10_thermal_study/baseline_19p6K"
find postProcessing -name wallHeatFlux.dat
python3 whf_extract.py postProcessing/wallHeatFlux/*/wallHeatFlux.dat \
  --case baseline_19p6K --Th 307.75 --Tc 288.15 \
  --output ../baseline_measurements.csv
python3 scripts/thermal_balance.py ../baseline_measurements.csv ../baseline_results.csv

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Read OpenFOAM 14 wallHeatFlux.dat files and prepare a thermal_balance.py row.

The wallHeatFlux function object writes, for every processed time and wall
patch: Time, patch, min [W/m2], max [W/m2], Q [W], q [W/m2]. In OpenFOAM 14 the
field is minus the boundary heat flux along the outward normal, so a positive Q
means heat entering the fluid through that patch. This tool converts the raw
signed values to the positive Q_hot / Q_cold convention of Guide 10, reports the
sum over all wall patches and, when several times are present, the drift of Q
over the last records. It never overwrites an existing output file.
"""
import argparse
import csv
from collections import OrderedDict
from pathlib import Path
import sys


def read_dat(paths):
    data = OrderedDict()  # time -> {patch: Q}
    for path in paths:
        with open(path, encoding='utf-8') as f:
            for number, line in enumerate(f, 1):
                s = line.strip()
                if not s or s.startswith('#'):
                    continue
                parts = s.split()
                if len(parts) < 6:
                    raise ValueError(f'{path}:{number}: expected 6 columns, found {len(parts)}')
                t, patch, q_total = float(parts[0]), parts[1], float(parts[4])
                data.setdefault(t, {})[patch] = q_total
    if not data:
        raise ValueError('No data rows found')
    return OrderedDict(sorted(data.items()))


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument('dat', nargs='+', type=Path, help='wallHeatFlux.dat file(s)')
    p.add_argument('--case', required=True, help='case label for the output row')
    p.add_argument('--hot', default='hot', help='hot patch name')
    p.add_argument('--cold', default='cold', help='cold patch name')
    p.add_argument('--Th', type=float, required=True, help='hot wall temperature [K]')
    p.add_argument('--Tc', type=float, required=True, help='cold wall temperature [K]')
    p.add_argument('--area-hot', type=float, default=1.1336, help='hot wall area [m2]')
    p.add_argument('--area-cold', type=float, default=1.1336, help='cold wall area [m2]')
    p.add_argument('--L', type=float, default=0.076, help='characteristic length [m]')
    p.add_argument('--k', type=float, default=0.026086, help='fluid conductivity [W/(m K)]')
    p.add_argument('--last', type=int, default=10, help='records used for the drift check')
    p.add_argument('--output', type=Path, help='NEW csv for thermal_balance.py')
    a = p.parse_args()

    data = read_dat(a.dat)
    t_last, row = next(reversed(data.items()))
    for name in (a.hot, a.cold):
        if name not in row:
            p.error(f'patch {name} not found at time {t_last}; available: {", ".join(row)}')
    q_hot, q_cold = row[a.hot], -row[a.cold]
    print(f'Last time: {t_last:g}  ({len(data)} time records)')
    for name, value in row.items():
        print(f'  raw Q[{name}] = {value:+.6g} W')
    print(f'Q_hot  = {q_hot:.6g} W   (heat entering the fluid)')
    print(f'Q_cold = {q_cold:.6g} W   (heat leaving the fluid)')
    print(f'Sum of signed Q over all wall patches = {sum(row.values()):+.6g} W')
    if q_hot <= 0 or q_cold <= 0:
        print('WARNING: unexpected sign; check patch names and the sign convention before continuing')
    if len(data) > 1:
        tail = list(data.items())[-max(2, a.last):]
        series = [r[a.hot] for _, r in tail if a.hot in r]
        ref = max(abs(v) for v in series)
        drift = (max(series) - min(series)) / ref * 100 if ref else float('nan')
        print(f'Drift of Q[{a.hot}] over the last {len(series)} records: {drift:.3g} %')
    if a.output:
        if a.output.exists():
            p.error('Output already exists; refusing to overwrite it')
        with a.output.open('x', newline='', encoding='utf-8') as f:
            w = csv.writer(f)
            w.writerow(['case', 'Th_K', 'Tc_K', 'Ahot_m2', 'Acold_m2', 'Qhot_W', 'Qcold_W', 'L_m', 'k_W_mK'])
            w.writerow([a.case, a.Th, a.Tc, a.area_hot, a.area_cold, f'{q_hot:.6g}', f'{q_cold:.6g}', a.L, a.k])
        print('Measurement row written to', a.output)


if __name__ == '__main__':
    try:
        main()
    except (OSError, ValueError) as exc:
        sys.exit(f'ERROR: {exc}')

Provato su un file costruito con il formato scritto dal sorgente (colonne Time, patch, min, max, Q, q): con Q(hot) = 36,5 W e Q(cold) = −35,8 W restituisce Q_hot = 36,5 W, Q_cold = 35,8 W, somma +0,7 W, e la riga prodotta dà in thermal_balance.py Nu_hot = 4,786 con k = 0,026086. Rifiuta un file di uscita esistente e una patch inesistente.

Sulla faccia di parete il flusso è calcolato come κ_eff ∂T/∂n con κ_eff = κ + c_p α_t,w, dove α_t,w è dato da alphatWallFunction con Pr_t = 0,85, e il gradiente normale è valutato fra il centro della prima cella e la parete. Se y⁺ della prima cella è dell’ordine di 1–2, ν_t alla parete è quasi nullo e il flusso è essenzialmente molecolare; se la mesh è più grossolana, la funzione di parete sostituisce la parte mancante del profilo con una legge universale, che per la convezione naturale è meno affidabile che per gli strati limite forzati. Il valore di Nu va sempre riportato insieme a y⁺ della prima cella e alla mesh usata.

cd "$HOME/OpenFOAM/guide10_thermal_study/baseline_19p6K"
foamPostProcess -solver fluid -latestTime -func yPlus > logs/yPlus.log 2>&1
grep -A4 "patch hot\|patch cold" logs/yPlus.log

La function object yPlus stampa minimo, massimo e media per ogni patch di parete e scrive il campo yPlus visualizzabile in ParaView lungo l’altezza. Il nome del log è nuovo, quindi non sovrascrive nulla; se esiste già, scegline un altro.

• Seleziona nel lettore OpenFOAM solo la patch (nell’elenco «Mesh Regions» le patch compaiono con il prefisso patch/), oppure isolala con Extract Block: integrare l’intero volume produce un numero senza significato.
• Integra l’array di celle wallHeatFlux: il lettore crea per default anche valori ai punti interpolati, che sul bordo non coincidono con quelli di faccia.
• Lascia disattivata l’opzione «Divide Cell Data By Volume»: con l’opzione attiva il filtro restituisce una media, non un integrale.
• Il filtro riporta anche Area: per la parete calda deve valere 1,1336 m². Un valore diverso rivela una selezione sbagliata.
• Il valore integrato deve coincidere, a meno di arrotondamenti, con il Q stampato da wallHeatFlux per la stessa patch e lo stesso tempo.

Dopo aver generato i profili con la function object sample del tutorial, lo script confronta T e velocità verticale con i dati di Betts e Bokhari interpolando il calcolo nelle ascisse sperimentali e riportando scarto medio, RMS e massimo per ogni quota:

cd "$HOME/OpenFOAM/guide10_thermal_study/baseline_19p6K"
foamPostProcess -latestTime -func sample > logs/sample.log 2>&1
ls postProcessing/sample/*/
python3 compare_expt.py . --output ../baseline_vs_expt.csv

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Compare buoyantCavity sample profiles with the Betts & Bokhari data.

Reads postProcessing/sample/<latest time>/y0.<i>.xy written by
'foamPostProcess -latestTime -func sample' (columns: x T Ux Uy Uz, as assumed by
the tutorial's validation/createGraphs script) and the experimental files
validation/exptData/mt_z0_<i>0_lo.dat (x [mm], T [degC]) and
mv_z0_<i>0_lo.dat (x [mm], Uy [m/s]). The simulated profile is linearly
interpolated at each experimental abscissa. Prints mean, RMS and maximum
deviations per height and optionally writes them to a NEW CSV file.
"""
import argparse
import csv
import math
from pathlib import Path
import sys


def read_columns(path):
    rows = []
    with open(path, encoding='utf-8') as f:
        for line in f:
            s = line.strip()
            if not s or s.startswith('#'):
                continue
            rows.append([float(v) for v in s.split()])
    if not rows:
        raise ValueError(f'No numeric data in {path}')
    return rows


def interp(xs, ys, x):
    if x < xs[0] or x > xs[-1]:
        return None
    for i in range(1, len(xs)):
        if x <= xs[i]:
            x0, x1, y0, y1 = xs[i - 1], xs[i], ys[i - 1], ys[i]
            return y0 if x1 == x0 else y0 + (y1 - y0) * (x - x0) / (x1 - x0)
    return ys[-1]


def stats(diffs):
    n = len(diffs)
    mean = sum(diffs) / n
    rms = math.sqrt(sum(d * d for d in diffs) / n)
    worst = max(diffs, key=abs)
    return n, mean, rms, worst


def latest_time_dir(root):
    dirs = []
    for d in root.iterdir():
        try:
            dirs.append((float(d.name), d))
        except ValueError:
            continue
    if not dirs:
        raise ValueError(f'No time directories in {root}')
    return max(dirs)[1]


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument('case', type=Path, help='case directory containing postProcessing/sample')
    p.add_argument('--expt', type=Path, help='experimental data directory (default: <case>/validation/exptData)')
    p.add_argument('--output', type=Path, help='NEW csv summary file')
    a = p.parse_args()
    expt = a.expt or a.case / 'validation' / 'exptData'
    sdir = latest_time_dir(a.case / 'postProcessing' / 'sample')
    print('Simulation profiles from', sdir)
    if a.output and a.output.exists():
        p.error('Output already exists; refusing to overwrite it')
    results = []
    for i in range(1, 10):
        sim_file = sdir / f'y0.{i}.xy'
        if not sim_file.is_file():
            continue
        sim = sorted(read_columns(sim_file))
        xs = [r[0] for r in sim]
        for kind, col, fname, conv, unit in (('T', 1, f'mt_z0_{i}0_lo.dat', 273.15, 'K'),
                                              ('Uy', 3, f'mv_z0_{i}0_lo.dat', 0.0, 'm/s')):
            efile = expt / fname
            if not efile.is_file():
                continue
            ys = [r[col] for r in sim]
            diffs = []
            for row in read_columns(efile):
                value = interp(xs, ys, row[0] / 1000.0)
                if value is not None:
                    diffs.append(value - (row[1] + conv))
            if not diffs:
                continue
            n, mean, rms, worst = stats(diffs)
            results.append(dict(y_over_H=f'0.{i}', quantity=kind, unit=unit, points=n,
                                mean_dev=f'{mean:.4g}', rms_dev=f'{rms:.4g}', max_abs_dev=f'{worst:.4g}'))
            print(f'y/H=0.{i}  {kind:2s}: n={n:2d}  mean={mean:+.4g} {unit}  rms={rms:.4g} {unit}  worst={worst:+.4g} {unit}')
    if not results:
        p.error('No matching simulation/experiment pairs found')
    if a.output:
        with a.output.open('x', newline='', encoding='utf-8') as f:
            w = csv.DictWriter(f, fieldnames=list(results[0]))
            w.writeheader()
            w.writerows(results)
        print('Summary written to', a.output)
    print('Deviation = simulation - experiment. The data come from the tutorial; cite the original publication.')


if __name__ == '__main__':
    try:
        main()
    except (OSError, ValueError) as exc:
        sys.exit(f'ERROR: {exc}')

Lo script assume il formato a colonne «x T Ux Uy Uz» dei file y0.N.xy, lo stesso usato da validation/createGraphs del tutorial (colonna 4 per Uy): se ls mostra nomi o colonne diversi, adegua gli indici prima di usarlo. È stato provato su profili costruiti ad arte e sui file sperimentali reali del tutorial. Una deviazione di pochi decimi di kelvin al centro e di 1–2 K vicino alle pareti è tipica dei calcoli RANS su questo caso; scarti di molti kelvin indicano un problema.

• Per confrontare a vista con i dati sperimentali usa Plot Over Line da (0, y, 0) a (0,076, y, 0) con y = 1,09 m (y/H = 0,5): sono le stesse coordinate della linea y0.5 di system/sample.
• Imposta la stessa scala di colore (per esempio 288–318 K) per entrambi i casi e salvala come preset: confronti fra immagini con scale diverse non hanno valore.
• Per la stratificazione, un Plot Over Line verticale al centro (x = 0,038 m, z = 0) mostra l’andamento di T con l’altezza: nelle cavità alte il nucleo è stratificato in modo approssimativamente lineare.
• Salva lo stato con File → Save State (.pvsm) accanto ai CSV esportati.

Un controllo indipendente si ottiene dai dati del tutorial: a metà altezza vicino alla parete fredda la temperatura misurata passa da 15,0 °C (parete) a 15,62 °C a 0,33 mm, poi 15,94 °C a 0,58 mm e 16,45 °C a 0,84 mm. Il gradiente vicino alla parete è quindi fra 1,3 e 2 K/mm, cioè un Nu locale fra circa 5 e 7,5. È coerente con Nu_W ≈ 5,65 della correlazione e con il valore sintetico 4,79 della figura: l’ordine di grandezza è giusto, i numeri precisi vanno ottenuti dal calcolo. Il valore medio misurato da Betts e Bokhari per questo ΔT, riportato dalla scheda ERCOFTAC, è Nu = 5,85: è il riferimento sperimentale con cui confrontare Nu_hot e Nu_cold, dopo aver verificato nella fonte che la definizione (lunghezza W, conducibilità, ΔT) coincida con quella della guida.

Per Nu = Q L / (A κ ΔT), con grandezze indipendenti, l’incertezza relativa è:

(dNu/Nu)^2 = (dQ/Q)^2 + (dL/L)^2 + (dA/A)^2 + (dk/k)^2 + (d(DeltaT)/DeltaT)^2

In una simulazione L, A e ΔT sono esatti per costruzione; κ è un dato di ingresso (incertezza di modello, non numerica); il termine dominante è dQ/Q, da stimare come combinazione di: scarto fra Q_hot e Q_cold, deriva iterativa delle ultime iterazioni e differenza fra mesh (sezione 13). Per il valore centrale usa la media (Q_hot + Q_cold)/2 e riporta metà della differenza come una delle componenti dell’incertezza.

In ParaView, sulla patch hot, usa un Calculator sui dati di cella con l’espressione wallHeatFlux*0.076/(0.026086*19.6) e chiamala Nu_local; poi Plot Data o Plot Over Line lungo y a z = 0. La curva attesa è alta in basso, dove lo strato limite caldo nasce, e decresce verso l’alto; per la parete fredda l’andamento è speculare (e i valori sono negativi per la convenzione di segno). La media d’area di Nu_local deve coincidere con Nu_hot calcolato da Q.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Evaluate signed wall-power magnitudes, balance error, and cavity Nusselt number."""
import argparse
import csv
from pathlib import Path
import sys

REQUIRED = ('case', 'Th_K', 'Tc_K', 'Ahot_m2', 'Acold_m2', 'Qhot_W', 'Qcold_W', 'L_m', 'k_W_mK')


def evaluate(row):
    label = row['case'].strip()
    if not label:
        raise ValueError('Empty case name')
    values = {}
    for col in REQUIRED[1:]:
        values[col] = float(row[col])
    th, tc = values['Th_K'], values['Tc_K']
    ah, ac = values['Ahot_m2'], values['Acold_m2']
    qh, qc = values['Qhot_W'], values['Qcold_W']
    gap, k = values['L_m'], values['k_W_mK']
    if not th > tc or min(ah, ac, qh, qc, gap, k) <= 0:
        raise ValueError('Expected Th > Tc and strictly positive areas, powers, gap, conductivity')
    denom = k * (th - tc)
    return dict(case=label, deltaT_K=f'{th-tc:.6g}', Nu_hot=f'{qh*gap/ah/denom:.6g}',
                Nu_cold=f'{qc*gap/ac/denom:.6g}',
                relative_imbalance_percent=f'{100*abs(qh-qc)/max(qh,qc):.6g}',
                Qhot_W=f'{qh:.6g}', Qcold_W=f'{qc:.6g}')


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('input_csv', type=Path)
    parser.add_argument('output_csv', type=Path)
    args = parser.parse_args()
    if args.output_csv.exists():
        parser.error('Output already exists; refusing to overwrite it')
    with args.input_csv.open(newline='', encoding='utf-8-sig') as f:
        reader = csv.DictReader(f)
        if not set(REQUIRED).issubset(set(reader.fieldnames or [])):
            parser.error('CSV missing required columns: ' + ', '.join(REQUIRED))
        rows = [evaluate(row) for row in reader]
    if not rows:
        parser.error('Empty data table')
    with args.output_csv.open('x', newline='', encoding='utf-8') as f:
        writer = csv.DictWriter(f, fieldnames=list(rows[0]))
        writer.writeheader()
        writer.writerows(rows)
    for row in rows:
        print(f"{row['case']}: Nu_hot={row['Nu_hot']}, Nu_cold={row['Nu_cold']}, balance error={row['relative_imbalance_percent']}%")
    print('Results saved:', args.output_csv)


if __name__ == '__main__':
    try:
        main()
    except (OSError, ValueError, KeyError) as exc:
        sys.exit(f'ERROR: {exc}')

• Colonne obbligatorie: case, Th_K, Tc_K, Ahot_m2, Acold_m2, Qhot_W, Qcold_W, L_m, k_W_mK; il BOM di Excel è gestito da utf-8-sig.
• Rifiuta T_h ≤ T_c e qualsiasi area, potenza, distanza o conducibilità non positiva: per questo Q_cold va inserito già positivo.
• Il file di uscita è aperto in modalità 'x': se esiste, l’operazione fallisce invece di sovrascrivere.
• Eseguito sul file sintetico del pacchetto produce un file identico a examples/derived_SYNTHETIC.csv (Nu_hot 4,78722 e 5,29766, scarti 1,91781 % e 0,983607 %); sul modello vuoto termina con «Empty data table», come indicato nella guida.

examples/measurement_TEMPLATE.csv

case,Th_K,Tc_K,Ahot_m2,Acold_m2,Qhot_W,Qcold_W,L_m,k_W_mK

examples/example_MEASUREMENTS_SYNTHETIC.csv (dati inventati)

case,Th_K,Tc_K,Ahot_m2,Acold_m2,Qhot_W,Qcold_W,L_m,k_W_mK
baseline_SYNTHETIC,307.75,288.15,1.1336,1.1336,36.5,35.8,0.076,0.02608
hotter_SYNTHETIC,317.75,288.15,1.1336,1.1336,61.0,60.4,0.076,0.02608

examples/derived_SYNTHETIC.csv (risultato atteso del calcolatore)

case,deltaT_K,Nu_hot,Nu_cold,relative_imbalance_percent,Qhot_W,Qcold_W
baseline_SYNTHETIC,19.6,4.78722,4.69541,1.91781,36.5,35.8
hotter_SYNTHETIC,29.6,5.29766,5.24555,0.983607,61,60.4

python scripts/thermal_balance.py examples/example_MEASUREMENTS_SYNTHETIC.csv examples/test_results.csv
fc examples\test_results.csv examples\derived_SYNTHETIC.csv        # Windows
diff examples/test_results.csv examples/derived_SYNTHETIC.csv      # Linux

Se il confronto non riporta differenze, il calcolatore e l’interprete Python funzionano come previsto. Usa poi una cartella separata per i risultati reali (per esempio results/), così i file sintetici non si mescolano mai con le misure.

In una copia nuova si può infittire solo vicino alle pareti termiche con una gradazione a più tratti nel blockMeshDict. Esempio, in cui ogni metà della larghezza ha metà delle celle e un rapporto 4 fra cella centrale e cella di parete:

blocks
(
    hex (0 1 2 3 4 5 6 7) (50 150 15)
    simpleGrading
    (
        ((0.5 0.5 4) (0.5 0.5 0.25))   // x: refined at both walls
        1
        1
    )
);

Ogni terna è (frazione della lunghezza, frazione delle celle, rapporto fra ultima e prima cella del tratto). Prima di ogni calcolo controlla checkMesh e la dimensione della prima cella. Una sequenza di tre mesh con raffinamento costante (per esempio 35, 50 e 70 celle in x, gradazione invariata) permette di applicare il metodo di Richardson e il GCI della Guida 4 a Q_hot. Mantieni costante il numero di celle in y e z nella prima serie: cambiare tutto insieme confonde gli effetti.

A Ra_H ≈ 2·10¹⁰ il flusso nella cavità alta può presentare instabilità e fluttuazioni a bassa frequenza. Un calcolo stazionario RANS può allora oscillare indefinitamente fra stati vicini: residui fermi su un plateau e Q_hot che oscilla di qualche percento. In questo caso non si stringono i fattori di rilassamento fino a «forzare» la convergenza; si riporta l’ampiezza dell’oscillazione come incertezza iterativa oppure, in una copia nuova, si passa a un calcolo transitorio (ddtSchemes Euler o backward, passo controllato dal numero di Courant come nella Guida 8) e si mediano le grandezze su un intervallo statisticamente stazionario.

Cambia un fattore per volta, ognuno in una copia propria, e riporta i risultati in una tabella unica. Con l’esponente 0,293 della correlazione, passando da ΔT = 19,6 a 29,6 K ci si aspetta Nu_W più alto di circa il 13 %: un aumento molto diverso va spiegato.

Nell’estensione a un condotto riscaldato esistono soluzioni di riferimento che permettono una vera verifica del codice, non solo una plausibilità:

• flusso laminare completamente sviluppato in tubo circolare: Nu_D = 3,66 a temperatura di parete costante e 4,36 a flusso di parete costante;
• canale piano fra due lastre: Nu_Dh = 7,54 (temperatura costante) e 8,23 (flusso costante) su entrambe le pareti;
• lunghezza d’ingresso termica laminare ≈ 0,05 Re Pr D_h: il dominio deve essere più lungo per raggiungere il valore sviluppato;
• in turbolento, correlazioni come Gnielinski hanno incertezze del 10–20 % e servono solo come confronto di ordine di grandezza.

Il bilancio ṁ c_p (T_out − T_in) va calcolato con la temperatura media pesata dalla portata, non con la media d’area: in OpenFOAM 14 si ottiene con una surfaceFieldValue con operation average e weightField phi (procedura descritta nella Guida 11).

Si elencano tutti con ls $FOAM_TUTORIALS/fluid $FOAM_TUTORIALS/multiRegion/CHT. Il file system/functions di ciascuno è un buon modello per le proprie misure.

cd "$HOME/OpenFOAM"
foamVersion > guide10_thermal_study/OPENFOAM_VERSION.txt 2>&1
tar czf guide10_thermal_study_$(date +%Y%m%d).tar.gz guide10_thermal_study
sha256sum guide10_thermal_study_*.tar.gz > guide10_thermal_study_SHA256.txt

L’archivio è una copia con impronta: nessun file viene rimosso. Conserva insieme i CSV di misura, i risultati di thermal_balance.py, il confronto sperimentale e gli stati ParaView.

• Q delle patch adiabatiche pari a zero e somma con segno di tutte le patch prossima a zero.
• Riga «Selecting … thermophysical transport model» letta e Pr_t dichiarato (parete 0,85, volume 1 se non diversamente impostato).
• Numero di iterazioni e motivo dell’arresto (convergenza o endTime) registrati.
• y⁺ della prima cella sulle pareti termiche riportato.
• Q confrontato con la stima da correlazione (≈ 43 W nel caso base) e scostamento spiegato.
• Profili confrontati con i dati sperimentali del tutorial, con scarti per quota.
• Effetto dell’irraggiamento discusso, almeno come ordine di grandezza.
• Commento della variante corretto o annotato in hotter_29p6K/0/T.

• P. L. Betts, I. H. Bokhari, Experiments on turbulent natural convection in an enclosed tall cavity, International Journal of Heat and Fluid Flow 21(6) (2000) 675–683, doi:10.1016/S0142-727X(00)00033-3: l’esperimento del tutorial (citato nel suo README);
• Scheda del caso nella banca dati ERCOFTAC: KBwiki UFR 4-10, natural convection in a tall cavity (geometria, Ra, Nusselt medi misurati, dati scaricabili);
• S. M. ElSherbiny, G. D. Raithby, K. G. T. Hollands, Heat transfer by natural convection across vertical and inclined air layers, Journal of Heat Transfer 104(1) (1982) 96–102, ASME Digital Collection: la correlazione usata nelle stime;
• I. B. Celik, U. Ghia, P. J. Roache, C. J. Freitas, H. Coleman, P. E. Raad, Procedure for estimation and reporting of uncertainty due to discretization in CFD applications, Journal of Fluids Engineering 130 (2008) 078001, testo pubblicato da ASME: il GCI usato nello studio di mesh;
• T. L. Bergman, A. S. Lavine, F. P. Incropera, D. P. DeWitt, Fundamentals of Heat and Mass Transfer, 8th ed., Wiley, 2017 (ISBN 978-1-119-35388-1), Wiley: convezione naturale in cavità, irraggiamento fra superfici, convezione forzata in condotti;
• A. Bejan, Convection Heat Transfer, 4th ed., Wiley, 2013 (ISBN 978-0-470-90037-6): analisi di scala degli strati limite di convezione naturale;
• F. R. Menter, Two-equation eddy-viscosity turbulence models for engineering applications, AIAA Journal 32(8) (1994) 1598–1605: il modello k-ω SST del tutorial (Guida 9).

Il README del tutorial indica anche la pagina del database ERCOFTAC/Manchester (caso 79) da cui provengono i dati sperimentali.

Ogni affermazione sul comportamento del codice si può verificare direttamente nel repository ufficiale della Foundation:

• tutorials/fluid/buoyantCavity — il caso completo, con validation/exptData e system/sample
• applications/modules/fluid/thermophysicalPredictor.C — equazione dell’energia del modulo fluid
• src/finiteVolume/cfdTools/general/buoyancy/buoyancy.H — definizione di p_rgh, gh, pRef e hRef
• src/functionObjects/field/wallHeatFlux/wallHeatFlux.C — segno del flusso e potenza integrata Q
• RASThermophysicalTransportModel.C + unityLewisEddyDiffusivity.C — modello predefinito e Pr_t = 1
• externalTemperatureFvPatchScalarField.H — opzioni Q, q, h, Ta, emissivity
• src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H — entalpia sensibile e T_ref
• foamPostProcess.C — l’opzione -solver e la costruzione dei modelli
• etc/caseDicts/functions — le configurazioni richiamate da #includeFunc e -func