OpenFOAM: trasferimento termico coniugato (CHT) fra fluido e solidi

I tre file in constant/geometry (heatedDuct.stl, fluidToMetal.stl, metalToHeater.stl) descrivono una geometria semplice, con l’asse del condotto lungo z:

L’interfaccia fluido-metallo è il perimetro interno del canale (4 × 40 mm × 150 mm = 0,024 m²); quella metallo-riscaldatore è la faccia superiore del riscaldatore (35 × 150 mm = 5,25·10⁻³ m²). Il calore generato nel riscaldatore deve quindi attraversare prima una sezione piccola e poi distribuirsi lungo la parete metallica: la conduzione «circonferenziale» nel metallo è la chiave fisica del caso. Con la sorgente del tutorial la potenza nominale è Q̇_gen = 10⁷ × 5,25·10⁻⁵ = 525 W; quella reale dipende dal volume misurato sulla mesh (sezione 7).

Per un’interfaccia ideale fra fluido (f) e solido (s), con normale n uscente dal fluido:

T_f = T_s                                continuità della temperatura
kappa_f dT_f/dn = kappa_s dT_s/dn        continuità del flusso normale

In una discretizzazione ai volumi finiti con le due celle adiacenti a distanza δ_f e δ_s dalla faccia, queste due condizioni danno una temperatura d’interfaccia pesata dalle conduttanze:

T_i = (kappa_f/delta_f * T_f + kappa_s/delta_s * T_s) / (kappa_f/delta_f + kappa_s/delta_s)
q   = (T_s - T_f) / (delta_f/kappa_f + delta_s/kappa_s)

È la «media armonica» delle conducibilità di Patankar: il flusso è limitato dalla resistenza maggiore. Con acqua (κ = 0,61) e alluminio (κ = 200) e celle uguali, la temperatura d’interfaccia è praticamente quella del metallo: la resistenza è tutta nel fluido. Per questo la qualità della mesh dal lato del fluido decide l’accuratezza del flusso scambiato.

A differenza della Guida 10, i campi di heatedDuct dichiarano le dimensioni come vettore di sette esponenti [massa lunghezza tempo temperatura quantità corrente intensità luminosa]:

Una sorgente in W/m³ ha dimensioni [1 −1 −3 0 0 0 0], una potenza in W [1 2 −3 0 0 0 0]: confrontare gli esponenti è il modo più rapido per accorgersi di aver scambiato q‴ e Q̇.

I commenti nei file del tutorial associano il metallo all’alluminio e il riscaldatore al ferro; le proprietà sono costanti e arrotondate. La cautela della guida resta valida: modificando solo κ non si ottiene un altro materiale reale (riquadro della sezione 5).

La conseguenza è netta: a t = 20 s l’acqua entrata dall’ingresso ha percorso solo 20 mm, e quella che esce è acqua presente nel condotto fin dall’inizio, riscaldata solo vicino alle pareti. La temperatura media di miscelamento all’uscita è quindi ancora molto vicina a 300 K e quasi tutta la potenza generata è accumulata nei solidi e nell’acqua vicina alle pareti. Il caso del tutorial è un transitorio iniziale, non uno stato stazionario: il bilancio della sezione 10 deve includere l’accumulo, esattamente come la guida raccomanda.

Con g = (0 0 0) e densità costante il moto non dipende dalla temperatura: il problema termico è disaccoppiato da quello fluidodinamico. In una prova reale con acqua riscaldata dal basso a 1 mm/s la situazione sarebbe opposta. Con β ≈ 2,76·10⁻⁴ 1/K e una differenza di 10 K su 40 mm, Gr = gβΔT D³/ν² ≈ 1,9·10⁶ e il numero di Richardson Gr/Re² ≈ 1000: la convezione naturale dominerebbe largamente quella forzata, con celle convettive trasversali. Il tutorial è quindi un dimostratore numerico del CHT, non un modello di un apparato reale. Per studiare la convezione mista servono una copia nuova con g diversa da zero, un’equazione di stato che faccia variare la densità e una mesh adeguata alle celle convettive.

cd "$FOAM_TUTORIALS/multiRegion/CHT/heatedDuct"
foamDictionary system/controlDict -entry regionSolvers
foamDictionary system/controlDict -entry endTime
foamDictionary constant/metal/physicalProperties -entry mixture/transport/kappa
foamDictionary constant/heater/fvModels -entry energySource/q
foamDictionary constant/fluid/g -entry value

foamDictionary legge e stampa una voce, con il percorso separato da «/». Usato senza opzioni di scrittura non modifica nulla: è il modo più sicuro per controllare i parametri di tutorial e copie prima di ogni esecuzione.

• blockMesh: un blocco unico 20 × 24 × 60 = 28 800 celle cubiche da 2,5 mm che contiene tutta la geometria (x 0–50, y −10–50, z 0–150 mm, con 0,001 mm di margine).
• snappyHexMesh: raffina al livello 1 (celle da 1,25 mm) vicino alle superfici esterne e alle due interfacce, aggancia la mesh alle superfici STL, non aggiunge strati (addLayers false). Le superfici fluidToMetal e metalToHeater creano le faceZone omonime e, con mode insidePoint, le cellZone metal e heater; il punto insidePoint (0.025 0.025 0.075) sta nel fluido.
• splitMeshRegions -cellZones all -defaultRegion fluid: separa una mesh per ogni cellZone, assegna le celle rimanenti alla regione fluid e crea le patch d’interfaccia con nomi <regione>_to_<regione>.

Con celle da 1,25 mm la parete metallica di 5 mm è descritta da circa 4 celle nello spessore e il riscaldatore da circa 8. Per un primo studio va bene; per misure di flusso all’interfaccia è la risoluzione dal lato acqua quella da verificare (riquadro della sezione 1).

cd "$HOME/OpenFOAM/CHT_Guida11/reference"
grep -A8 "fluid_to_metal" constant/fluid/polyMesh/boundary
grep -A8 "metal_to_fluid" constant/metal/polyMesh/boundary
grep -A8 "metal_to_heater" constant/metal/polyMesh/boundary
grep -A8 "heater_to_metal" constant/heater/polyMesh/boundary

Ogni patch d’interfaccia deve essere di un tipo «mapped» che indica la regione e la patch confinanti, e il numero di facce (nFaces) delle due patch di una coppia deve essere uguale: le regioni derivano dalla stessa mesh e le facce coincidono. Numeri diversi o una regione confinante errata spiegano temperature discontinue all’interfaccia (sezione 12).

L’Allrun ufficiale decompone in 4 sottodomini con method hierarchical; n (1 1 4); cioè taglia il condotto in quattro tronchi lungo z. Per usarlo in una copia nuova, dopo la mesh:

decomposePar -allRegions > guide11_logs/decomposePar.log 2>&1
mpirun -np 4 foamMultiRun -parallel > guide11_logs/solve_parallel.log 2>&1
reconstructPar -allRegions > guide11_logs/reconstructPar.log 2>&1

Con 28 800 celle di partenza il caso è piccolo e il guadagno parallelo è modesto; serve soprattutto come prova della procedura prima di mesh più fini. Non lanciare il percorso parallelo nella stessa copia in cui è stato eseguito solve_case.sh: i risultati seriali e ricostruiti finirebbero nelle stesse cartelle temporali.

Nel sorgente della versione 14 il modulo solid (applications/modules/solid/solid.C) risolve per ogni regione solida l’energia interna sensibile e:

fvm::ddt(rho, e) + thermophysicalTransport->divq(e) == fvModels().source(rho, e)

Con eConst e rhoConst è esattamente ρ C_v ∂T/∂t = ∇·(κ∇T) + q‴, dove q‴ arriva dal heatSource di fvModels. Il modulo fluid risolve l’entalpia sensibile h con convezione, energia cinetica e lavoro di pressione (Guida 10); con g = 0 manca il termine di galleggiamento. Le due famiglie di equazioni sono accoppiate solo attraverso le condizioni coupledTemperature sulle patch d’interfaccia.

foamMultiRun risolve le regioni una dopo l’altra in un ciclo PIMPLE esterno comune. Il passo è scelto dal minimo dei maxDeltaT dei moduli: il fluido lo limita con maxCo 1, mentre il modulo solid della versione 14 non ha un limite sul numero di diffusione e usa solo maxDeltaT di controlDict, se presente. Con U = 1 mm/s e celle da 1,25 mm, il Courant permette passi fino a circa 1 s. Il numero di Fourier di cella nell’alluminio diventa allora Fo = αΔt/Δx² ≈ 8,2·10⁻⁵ × 1/(1,25·10⁻³)² ≈ 50.

Lo schema di Eulero implicito resta stabile, ma con Fo ≫ 1 la storia temporale nei solidi è poco accurata, e con un solo passaggio del ciclo esterno l’accoppiamento fra regioni è ritardato di un passo. Per un transitorio da misurare (sezione 10) in una copia nuova si limitano il passo e si aumentano le correzioni esterne:

// system/controlDict
maxDeltaT       0.05;

// system/fvSolution (top level, not system/<region>/fvSolution)
PIMPLE
{
    nOuterCorrectors 3;
}

Il blocco PIMPLE di primo livello è quello che nel tutorial è vuoto; i tutorial CHT ufficiali della versione 14 vi impostano nOuterCorrectors. La verifica è sempre la stessa: dimezzare maxDeltaT e controllare che T_max del riscaldatore e l’energia accumulata non cambino oltre la tolleranza scelta.

Il valore di h è una stima d’ordine (κ_acqua diviso lo spessore di penetrazione di circa 1,7 mm). La resistenza convettiva domina: ridurre κ del metallo da 200 a 50 cambia poco la resistenza in serie. Cambia invece molto la distribuzione: le pareti laterali e superiore funzionano come alette alimentate dal basso. Con il parametro d’aletta m = √(h/(κ t)) e t = 5 mm, la parete laterale (45 mm) ha efficienza ≈ 0,82 con κ = 200 e ≈ 0,55 con κ = 50 (cht_estimates.py). Ipotesi da verificare: con il metallo meno conduttivo la parete inferiore e il riscaldatore diventano più caldi, le pareti laterali e superiore più fredde, e la potenza ceduta all’acqua si concentra sul fondo del canale.

Raccoglie in un solo programma tutte le stime di questa guida: volumi e capacità termiche, potenza, portata, Re, Pr, Pe, tempi caratteristici, lunghezze d’ingresso, tempi di diffusione, efficienza d’aletta e numero di Richardson ipotetico. I volumi nominali vanno sostituiti con quelli misurati (--V-heater):

python cht_estimates.py
python cht_estimates.py --kappa-metal 50
python cht_estimates.py --q 5e6 --V-heater 5.31e-5

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Hand estimates for the OpenFOAM 14 multiRegion/CHT/heatedDuct case.

Nominal geometry from the tutorial STL files: square duct 50 x 50 mm outside,
40 x 40 mm fluid channel (5 mm metal wall), heater 35 x 10 mm under the duct,
all 150 mm long. Properties from constant/*/physicalProperties. The numbers are
order-of-magnitude checks for the CFD results, not results themselves; the
volumes must be replaced by the values measured on the actual mesh
(checkMesh -region <name>, 'Total volume').
"""
import argparse
import math


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument('--q', type=float, default=1e7, help='heater source density [W/m3]')
    p.add_argument('--kappa-metal', type=float, default=200.0, help='metal conductivity [W/(m K)]')
    p.add_argument('--U', type=float, default=1e-3, help='inlet velocity [m/s]')
    p.add_argument('--V-heater', type=float, help='measured heater volume [m3] (default: nominal)')
    p.add_argument('--h', type=float, default=350.0, help='assumed fluid-side film coefficient [W/(m2 K)] for the fin estimate')
    p.add_argument('--endTime', type=float, default=20.0, help='simulated time [s]')
    a = p.parse_args()

    L = 0.150
    fluid = dict(rho=1000.0, c=4181.0, mu=959e-6, Pr=6.62, V=0.040 * 0.040 * L)
    fluid['k'] = fluid['mu'] * fluid['c'] / fluid['Pr']
    metal = dict(rho=2700.0, c=900.0, k=a.kappa_metal, V=(0.050 ** 2 - 0.040 ** 2) * L)
    heater = dict(rho=8000.0, c=450.0, k=80.0, V=a.V_heater or 0.035 * 0.010 * L)

    area_in = 0.040 * 0.040
    dh = 0.040
    nu = fluid['mu'] / fluid['rho']
    mdot = fluid['rho'] * a.U * area_in
    q_gen = a.q * heater['V']
    re = a.U * dh / nu
    pe = re * fluid['Pr']

    print('Region   V [m3]      rho*c*V [J/K]   alpha [m2/s]   effusivity [W s^0.5/(m2 K)]')
    total_c = 0.0
    for name, r in (('fluid', fluid), ('metal', metal), ('heater', heater)):
        cap = r['rho'] * r['c'] * r['V']
        total_c += cap
        alpha = r['k'] / (r['rho'] * r['c'])
        eff = math.sqrt(r['k'] * r['rho'] * r['c'])
        r['alpha'] = alpha
        print(f'{name:7s}  {r["V"]:.4e}  {cap:12.1f}    {alpha:.3e}      {eff:9.0f}')
    print(f'Water conductivity k = mu cp / Pr = {fluid["k"]:.4f} W/(m K)')
    print(f'Heater power Q_gen = q * V_heater = {q_gen:.1f} W')
    print(f'Mass flow rate = rho U A = {mdot:.4e} kg/s')
    print(f'Re (D_h = 40 mm) = {re:.1f},  Pr = {fluid["Pr"]},  Pe = {pe:.0f}')
    print(f'Flow-through time L/U = {L / a.U:.0f} s  (simulated time {a.endTime:g} s)')
    print(f'Hydrodynamic entry length ~ 0.05 Re D_h = {0.05 * re * dh:.3f} m')
    print(f'Thermal entry length ~ 0.05 Re Pr D_h = {0.05 * pe * dh:.3f} m  (duct length {L} m)')
    print(f'Steady bulk temperature rise Q_gen/(mdot cp) = {q_gen / (mdot * fluid["c"]):.1f} K')
    print(f'Lumped time constant C_total/(mdot cp) = {total_c / (mdot * fluid["c"]):.0f} s')
    print(f'Energy generated up to endTime = {q_gen * a.endTime:.0f} J')
    print(f'Adiabatic heater heating rate q/(rho c) = {a.q / (heater["rho"] * heater["c"]):.2f} K/s')
    for name, r, d in (('metal wall (5 mm)', metal, 0.005), ('heater (10 mm)', heater, 0.010), ('water (5 mm)', fluid, 0.005)):
        print(f'Diffusion time d^2/alpha, {name}: {d * d / r["alpha"]:.2f} s')
    print(f'Thermal penetration in water after endTime sqrt(alpha t) = {1000 * math.sqrt(fluid["alpha"] * a.endTime):.2f} mm')
    m = math.sqrt(a.h / (metal['k'] * 0.005))
    for label, length in (('side wall (45 mm)', 0.045), ('side + top half (70 mm)', 0.070)):
        ml = m * length
        print(f'Fin estimate, {label}: mL = {ml:.2f}, efficiency tanh(mL)/mL = {math.tanh(ml) / ml:.2f}')
    beta, dT, g = 2.76e-4, 10.0, 9.81
    gr = g * beta * dT * dh ** 3 / nu ** 2
    print(f'If gravity were active (beta = {beta} 1/K, dT = {dT} K): Gr = {gr:.2e}, Ri = Gr/Re^2 = {gr / re ** 2:.0f}')


if __name__ == '__main__':
    main()

Uscita per il caso di riferimento (estratto):

Region   V [m3]      rho*c*V [J/K]   alpha [m2/s]   effusivity [W s^0.5/(m2 K)]
fluid    2.4000e-04        1003.4    1.449e-07           1591
metal    1.3500e-04         328.1    8.230e-05          22045
heater   5.2500e-05         189.0    2.222e-05          16971
Heater power Q_gen = q * V_heater = 525.0 W
Mass flow rate = rho U A = 1.6000e-03 kg/s
Flow-through time L/U = 150 s  (simulated time 20 s)
Steady bulk temperature rise Q_gen/(mdot cp) = 78.5 K
Lumped time constant C_total/(mdot cp) = 227 s
Adiabatic heater heating rate q/(rho c) = 2.78 K/s
Diffusion time d^2/alpha, metal wall (5 mm): 0.30 s
Diffusion time d^2/alpha, heater (10 mm): 4.50 s
Diffusion time d^2/alpha, water (5 mm): 172.58 s

Per i solidi la capacità termica si legge sotto Cv, non Cp: con densità costante le due coincidono. Nel fluido l’entalpia di hConst è c_p(T − T_ref) più il contributo p/ρ dell’equazione di stato rhoConst, con T_ref = 298,15 K se non specificato (sorgente hConstThermo.C); il termine p/ρ è trascurabile qui perché le pressioni relative sono minime.

Quando due corpi a temperature diverse vengono messi a contatto, la temperatura dell’interfaccia è la media delle due pesata con le effusività. Il metallo ha un’effusività circa 14 volte quella dell’acqua: nei transitori impone la propria temperatura all’interfaccia e l’acqua si adegua. La diffusività dice invece quanto rapidamente il calore si propaga all’interno: nell’alluminio 570 volte più in fretta che nell’acqua, ed è per questo che la parete si uniforma in frazioni di secondo mentre l’acqua resta fredda a pochi millimetri dalla parete.

Nel sorgente della versione 14 il modello heatSource (src/fvModels/general/heatSource) richiede o la potenza totale Q [W] o la densità di potenza q [W/m³], entrambe funzioni del tempo, applicate alla cellZone indicata (all significa tutte le celle della regione). Con q la potenza effettiva dipende dal volume della mesh; con Q è fissata indipendentemente dalla mesh. In uno studio di mesh conviene quindi una copia con:

energySource
{
    type            heatSource;
    cellZone        all;
    Q               525;
}

così le differenze fra mesh non includono una variazione della potenza iniettata. Va scritta una sola delle due voci: con entrambe il modello non sa quale usare.

Se interessa un metallo reale meno conduttivo, si cambiano insieme tutte le proprietà. Per un acciaio al carbonio valori tipici da manuale sono κ ≈ 45–50 W/(m K), ρ ≈ 7850 kg/m³, c ≈ 460–490 J/(kg K). In una copia nuova:

// constant/metal/physicalProperties (mixture section)
equationOfState { rho 7850; }
transport       { kappa 50; }
thermodynamics  { hf 0; Cv 470; }

La capacità termica volumica sale da 2,43·10⁶ a 3,69·10⁶ J/(m³ K): nel transitorio il metallo si scalda più lentamente a parità di potenza. Confrontare questo caso con lowConductivity separa l’effetto della conducibilità da quello dell’accumulo, che il caso «ipotetico» mescola.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Create isolated CHT tutorial variants without changing existing work."""
import argparse
import os
from pathlib import Path
import re
import shutil
import sys

CASES = {
    "reference": None,
    "lowConductivity": ("constant/metal/physicalProperties", "kappa", "50"),
    "lowHeatSource": ("constant/heater/fvModels", "q", "5e6"),
}


def replace_dictionary_scalar(file_path: Path, key: str, value: str) -> None:
    text = file_path.read_text(encoding="utf-8")
    pattern = rf"(?m)^(\s*{re.escape(key)}\s+)([0-9eE.+-]+)(\s*;)"
    matches = list(re.finditer(pattern, text))
    if len(matches) != 1:
        raise ValueError(f"Expected exactly one {key} scalar entry in {file_path}, found {len(matches)}")
    updated = re.sub(pattern, lambda match: match.group(1) + value + match.group(3), text, count=1)
    file_path.write_text(updated, encoding="utf-8")


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--output", required=True, type=Path, help="Directory containing NEW case folders")
    parser.add_argument("--tutorial", type=Path, help="Optional explicit path to the official heatedDuct case")
    args = parser.parse_args()
    tutorials = os.environ.get("FOAM_TUTORIALS", "")
    tutorial = args.tutorial or (Path(tutorials) / "multiRegion/CHT/heatedDuct" if tutorials else None)
    if not tutorial or not tutorial.is_dir():
        parser.error("Official heatedDuct tutorial not found; source OpenFOAM and set FOAM_TUTORIALS or use --tutorial")
    expected = ["system/controlDict", "constant/metal/physicalProperties", "constant/heater/fvModels", "0/fluid/T"]
    if any(not (tutorial / name).is_file() for name in expected):
        parser.error("The template is not the expected OpenFOAM Foundation 14 heatedDuct case")
    out = args.output.expanduser().resolve()
    if any((out / case).exists() for case in CASES):
        parser.error("One or more destination case directories already exist; refusing to overwrite")
    out.mkdir(parents=True, exist_ok=True)
    created = []
    try:
        for name, change in CASES.items():
            destination = out / name
            shutil.copytree(tutorial, destination, symlinks=False)
            created.append(destination)
            if change:
                relative, key, value = change
                replace_dictionary_scalar(destination / relative, key, value)
            (destination / "GUIDE11_CASE.txt").write_text(
                "Guide 11 (OpenFOAM Foundation 14)\n"
                f"Case name: {name}\n"
                f"Official source: {tutorial.resolve()}\n"
                f"Modification: {change or 'none'}\n"
                "Generated from a copy; upstream tutorial was not modified.\n",
                encoding="utf-8",
            )
            print(f"Prepared {destination}")
    except Exception as exc:
        # Preserve all created data to support forensic review; never clean up silently.
        print(f"Preparation stopped: {exc}. Existing copies were left in place.", file=sys.stderr)
        return 1
    return 0


if __name__ == "__main__":
    sys.exit(main())

• replace_dictionary_scalar cerca una riga che inizi con la chiave (kappa o q) seguita da un numero e da «;» e pretende che ce ne sia esattamente una; altrimenti si ferma. Nei file ufficiali kappa 200; e q 1e7; compaiono una sola volta.
• Ogni copia riceve un file GUIDE11_CASE.txt con origine e modifica: gli script di mesh e soluzione lo pretendono, quindi non possono essere lanciati per errore sul tutorial installato.
• Se una qualsiasi delle tre destinazioni esiste, non viene creato nulla. In caso di errore a metà, le copie già create restano per l’analisi: lo script non cancella mai.
• Eseguito sui file del repository ufficiale di OpenFOAM 14, diff -r reference lowConductivity mostra soltanto kappa 200; → kappa 50; e il file descrittivo; diff -r reference lowHeatSource soltanto q 1e7; → q 5e6;. Una seconda esecuzione con la stessa destinazione termina con codice 2.

Il comando source /opt/OpenFOAM-14/etc/bashrc va adattato: con i pacchetti Ubuntu distribuiti da openfoam.org il percorso tipico è /opt/openfoam14/etc/bashrc (minuscolo, senza trattino); con una compilazione da sorgente è nella cartella dell’utente. Per scoprirlo e verificarlo:

ls -d /opt/openfoam* /opt/OpenFOAM* 2>/dev/null
source /opt/openfoam14/etc/bashrc     # adjust to what ls shows
foamVersion
echo "$FOAM_TUTORIALS"
ls "$FOAM_TUTORIALS/multiRegion/CHT/heatedDuct"

Lo script di preparazione va lanciato indicando il suo percorso (python3 scripts/prepare_cases.py dalla cartella del pacchetto, come nel README) oppure entrando prima nella cartella scripts, come nel comando della guida.

cd ~/guide11_kit
file scripts/*.sh          # must NOT report "with CRLF line terminators"
bash -n scripts/mesh_case.sh && bash -n scripts/solve_case.sh && echo "syntax OK"
python3 -m py_compile scripts/prepare_cases.py scripts/energy_balance.py && echo "python OK"

Se file segnala CRLF, in Notepad++ usa Modifica → Conversione fine riga → Unix (LF), salva e ritrasferisci con FileZilla in modalità binaria. bash -n controlla la sintassi senza eseguire nulla. py_compile crea una cartella __pycache__ accanto agli script: è innocua.

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Mesh a fresh OpenFOAM Foundation 14 heatedDuct case and retain all logs.
set -euo pipefail
if [[ $# -ne 1 ]]; then
    echo "Usage: bash mesh_case.sh /absolute/path/to/case" >&2
    exit 2
fi
case_dir="$(realpath -- "$1")"
[[ -f "$case_dir/GUIDE11_CASE.txt" ]] || { echo "Not a prepared Guide 11 case: $case_dir" >&2; exit 2; }
command -v blockMesh >/dev/null || { echo "Source the OpenFOAM 14 environment first" >&2; exit 2; }
command -v snappyHexMesh >/dev/null || exit 2
command -v splitMeshRegions >/dev/null || exit 2
command -v checkMesh >/dev/null || exit 2
cd "$case_dir"
if [[ -d constant/fluid/polyMesh || -d constant/metal/polyMesh || -d constant/heater/polyMesh ]]; then
    echo "Regional meshes already exist; refusing to rebuild the case in place" >&2
    exit 1
fi
mkdir -p guide11_logs
run_logged() {
    local output="$1"; shift
    if [[ -e "$output" ]]; then
        echo "Existing log $output; refusing to overwrite" >&2
        return 1
    fi
    printf 'Running:'; printf ' %q' "$@"; printf '\n'
    "$@" > >(tee "$output") 2> >(tee -a "$output" >&2)
}
run_logged guide11_logs/mesh_blockMesh.log blockMesh
run_logged guide11_logs/mesh_snappyHexMesh.log snappyHexMesh
run_logged guide11_logs/mesh_splitMeshRegions.log splitMeshRegions -cellZones all -defaultRegion fluid
for region in fluid metal heater; do
    if [[ ! -f "constant/$region/polyMesh/boundary" ]]; then
        echo "No region mesh boundary file for $region" >&2
        exit 1
    fi
    run_logged "guide11_logs/checkMesh_${region}.log" checkMesh -region "$region"
done
printf 'Regional meshes created and checked in %s\n' "$case_dir"

• Il controllo su GUIDE11_CASE.txt garantisce che il caso sia una copia preparata; quello sulle tre polyMesh impedisce di rifare la mesh sopra una esistente.
• run_logged stampa il comando, scrive l’uscita a video e nel log con tee (anche gli errori, con tee -a) e si rifiuta di riusare un log esistente.
• Dopo splitMeshRegions verifica che esista constant/<regione>/polyMesh/boundary per ciascuna regione prima di lanciare checkMesh -region.

cd "$HOME/OpenFOAM/CHT_Guida11/reference"
grep -E "cells:|Total volume" guide11_logs/checkMesh_fluid.log guide11_logs/checkMesh_metal.log guide11_logs/checkMesh_heater.log

checkMesh stampa il numero di celle e, nel controllo geometrico, il volume totale della regione. Atteso per il riscaldatore: circa 5,25·10⁻⁵ m³. La potenza effettiva è Q̇_gen = 10⁷ × V_heater: uno scarto del 2 % nel volume è uno scarto del 2 % nella potenza, da usare nel bilancio al posto dei 525 W nominali. Annota anche i volumi di fluid e metal: servono per le capacità termiche della sezione 10.

• «Mesh OK» alla fine di ciascuno dei tre log; in caso contrario leggere quali controlli falliscono.
• Non-ortogonalità massima: le celle tagliate da snappyHexMesh vicino alle interfacce sono le più distorte; oltre circa 70° servono correttori non ortogonali.
• Elenco delle patch: nel fluido fluidInlet, fluidOutlet, fluid_to_metal; nel metallo metalInlet, metalOutlet, metalExternal, metal_to_fluid, metal_to_heater; nel riscaldatore heaterInlet, heaterOutlet, heaterExternal, heater_to_metal.
• Patch d’ingresso, d’uscita o d’interfaccia con zero facce, oppure regioni con un numero di celle diverso fra le tre copie preparate (che usano la stessa geometria e lo stesso dizionario di mesh): indicano un problema di mesh, non di fisica.

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Run a meshed CHT tutorial case once; do not silently restart old solutions.
set -euo pipefail
if [[ $# -ne 1 ]]; then
    echo "Usage: bash solve_case.sh /absolute/path/to/case" >&2
    exit 2
fi
case_dir="$(realpath -- "$1")"
[[ -f "$case_dir/GUIDE11_CASE.txt" ]] || { echo "Not a prepared Guide 11 case" >&2; exit 2; }
command -v foamMultiRun >/dev/null || { echo "Source the OpenFOAM 14 environment first" >&2; exit 2; }
cd "$case_dir"
for region in fluid metal heater; do
    [[ -f "constant/$region/polyMesh/boundary" ]] || { echo "Missing $region mesh" >&2; exit 1; }
done
mkdir -p guide11_logs
log="guide11_logs/solve_foamMultiRun.log"
if [[ -e "$log" ]]; then
    echo "Solver log already exists; refusing to overwrite or restart" >&2
    exit 1
fi
# Use one serial solver process; the official Allrun is a separate MPI workflow.
foamMultiRun > >(tee "$log") 2> >(tee -a "$log" >&2)
printf 'Solver finished; log: %s/%s\n' "$case_dir" "$log"

Lo script controlla il file descrittivo, la disponibilità di foamMultiRun e le tre mesh regionali, poi esegue il calcolo seriale una sola volta, salvando uscita ed errori in guide11_logs/solve_foamMultiRun.log. Il commento finale ricorda che l’Allrun ufficiale è un percorso MPI distinto.

Le grandezze del bilancio si misurano meglio con OpenFOAM stesso, a ogni passo e sulla mesh effettiva. In una copia preparata e già con la mesh, prima di solve_case.sh, salva questo file come system/functions (OpenFOAM 14 lo legge automaticamente):

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

// Mass flow rate [kg/s]: in the fluid module phi is a mass flux.
#includeFunc patchFlowRate(name=mdotIn, region=fluid, patch=fluidInlet)
#includeFunc patchFlowRate(name=mdotOut, region=fluid, patch=fluidOutlet)

// Signed enthalpy flow sum(phi*h) [W]: negative on the inlet, positive on the outlet.
// Q_fluid = hFlowOut + hFlowIn
hFlowIn
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    log             yes;
    writeFields     false;
    region          fluid;
    patch           fluidInlet;
    operation       sum;
    weightField     phi;
    fields          (h);
}

hFlowOut
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    log             yes;
    writeFields     false;
    region          fluid;
    patch           fluidOutlet;
    operation       sum;
    weightField     phi;
    fields          (h);
}

// Mass-flow-weighted (bulk) temperatures [K]: sum(phi*T)/sum(phi).
TbulkIn
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    log             yes;
    writeFields     false;
    region          fluid;
    patch           fluidInlet;
    operation       average;
    weightField     phi;
    fields          (T);
}

TbulkOut
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    log             yes;
    writeFields     false;
    region          fluid;
    patch           fluidOutlet;
    operation       average;
    weightField     phi;
    fields          (T);
}

// Volume-averaged temperatures [K]: stored energy E = rho*c*V*Tavg + const.
#includeFunc volAverage(name=TavgFluid, region=fluid, T)
#includeFunc volAverage(name=TavgMetal, region=metal, T)
#includeFunc volAverage(name=TavgHeater, region=heater, T)

// Hottest cell of the heater [K].
#includeFunc cellMax(name=TmaxHeater, region=heater, T)

// Heat flow through every wall patch of each region [W], positive = into the region.
// Interface check: Q(fluid_to_metal) = -Q(metal_to_fluid), Q(metal_to_heater) = -Q(heater_to_metal).
#includeFunc wallHeatFlux(name=qFluid, region=fluid, writeFields=false, log=yes)
#includeFunc wallHeatFlux(name=qMetal, region=metal, writeFields=false, log=yes)
#includeFunc wallHeatFlux(name=qHeater, region=heater, writeFields=false, log=yes)

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

• patchFlowRate somma il flusso phi sulla patch: nel modulo fluid è un flusso di massa, quindi il risultato è ṁ in kg/s (negativo all’ingresso).
• Con operation sum e weightField phi la surfaceFieldValue calcola Σ φ_f h_f, il flusso di entalpia con segno. Q̇_fluid = hFlowOut + hFlowIn è la potenza netta portata via dall’acqua, già corretta per eventuali riflussi.
• Con operation average e weightField phi si ottiene Σ φT/Σ φ: la temperatura media pesata dalla portata della sezione 9.
• wallHeatFlux per regione dà la potenza su ogni patch di parete, positiva se entra nella regione: Q(fluid_to_metal) è il calore che entra nell’acqua, e deve valere −Q(metal_to_fluid).

Operazioni e pesi sono verificati nel sorgente di surfaceFieldValue della versione 14; la sintassi riprende i tutorial ufficiali wallBoiling, hotWire e vacuumFurnace. All’avvio il log deve elencare la creazione di ogni function object e, a ogni passo, i valori. I file finiscono sotto postProcessing/: find postProcessing -name "*.dat" li elenca tutti.

L="$HOME/OpenFOAM/CHT_Guida11/reference/guide11_logs/solve_foamMultiRun.log"
grep -c "^Time = " "$L"                    # passi eseguiti
grep "^deltaT = " "$L" | tail -3            # passo attuale
grep "Courant Number" "$L" | tail -3
grep "Solving for h" "$L" | tail -3        # fluido
grep "Solving for e" "$L" | tail -6        # metallo e riscaldatore
grep -i -E "warning|error|nan" "$L" | head

Il passo deve crescere da 10⁻³ s fino al limite imposto dal Courant e poi stabilizzarsi. «Solving for e» compare due volte per passo, una per regione solida. Il tempo di calcolo per secondo simulato, letto da ExecutionTime, permette di stimare la durata di una copia con endTime più lungo prima di lanciarla.

• Crea il file vuoto di apertura nella cartella del caso: touch reference.foam. Aprendolo da ParaView collegato a pvserver il lettore OpenFOAM mostra in «Mesh Regions» le mesh interne e le patch di ogni regione, con il nome della regione come prefisso.
• Seleziona le tre mesh interne e premi Apply; con Extract Block separa le regioni per colorarle in modo indipendente. Il campo si chiama T in tutte e tre: dopo l’estrazione è chiaro a quale regione appartiene.
• Il caso è scritto in binario (writeFormat binary): ParaView lo legge senza problemi, ma i file non si possono ispezionare con un editor di testo.
• Per un calcolo eseguito in parallelo e non ricostruito, imposta «Case Type» su Decomposed Case.

Con Plot Over Line da (0,025, −0,010, 0,075) a (0,025, 0,045, 0,075) si attraversano, a metà lunghezza, il riscaldatore, il fondo metallico e l’acqua. La continuità del flusso impone che le pendenze ai due lati di un’interfaccia stiano nel rapporto inverso delle conducibilità: all’interfaccia metallo-acqua il gradiente nell’acqua deve essere circa 330 volte quello nell’alluminio (200/0,606). Nel grafico la temperatura appare quasi piatta nel metallo e ripida nell’acqua: se il salto fra le due regioni non è continuo, l’accoppiamento non funziona.

La formula della guida coincide con quanto calcolano le function object TbulkIn e TbulkOut (sezione 8): sulle facce di bordo φ_f = ρ (U·n) A_f, e Σ φ_fT_f/Σ φ_f è la versione discreta di ∫ρT(U·n)dA/∫ρ(U·n)dA. Per un tempo già scritto, senza rilanciare il calcolo, si può usare foamPostProcess sulla sola regione fluida, con il solutore esplicito perché phi sia disponibile:

cd "$HOME/OpenFOAM/CHT_Guida11/reference"
foamPostProcess -region fluid -solver fluid -latestTime \
  -func "patchFlowRate(patch=fluidOutlet)" > guide11_logs/post_mdot_outlet.log 2>&1
grep -A3 "patchFlowRate" guide11_logs/post_mdot_outlet.log

Il modulo -solver è necessario anche per wallHeatFlux (Guida 10): senza, foamPostProcess legge i campi ma non costruisce il modello di trasporto termico. Per le regioni solide si usa -region metal -solver solid.

Due grandezze spesso confuse vanno quindi tenute distinte: il calore che attraversa la parete verso l’acqua (misurato da wallHeatFlux) e il calore che l’acqua porta fuori dal dominio (ṁ c_pΔT). La loro differenza è l’energia accumulata nell’acqua. Solo in regime stazionario coincidono, e coincidono con Q̇_gen.

Lo stazionario richiede alcune costanti di tempo globali (≈ 227 s) più il tempo di attraversamento (150 s): un endTime di 1500–2000 s in una copia nuova. In quelle condizioni il bilancio darebbe ΔT̄ = Q̇_gen/(ṁ c_p) = 525/(1,6·10⁻³ × 4181) ≈ 78,5 K, cioè acqua in uscita a circa 378 K: oltre l’ebollizione a pressione atmosferica, con pareti ancora più calde. Il modello a densità costante e monofase non può rappresentarlo. Per uno studio stazionario fisicamente sensato si riduce la potenza (per esempio q 1e6, ΔT̄ ≈ 7,9 K) o si aumenta la velocità d’ingresso, sempre in una copia nuova e documentando la scelta.

Nota sulla figura sintetica: la curva «uscita fluido» sale a 312 K in 20 s, mentre le stime indicano una temperatura di miscelamento all’uscita quasi invariata in quell’intervallo. La figura illustra la forma di un transitorio, non una previsione per questo caso.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Evaluate signed CHT energy balance from measured CSV rows, not CFD files."""
import argparse
import csv
from pathlib import Path
import sys

REQUIRED = ["case", "time_s", "m_dot_kg_s", "T_in_K", "T_out_K", "cp_J_kgK", "q_gen_W", "q_ext_W", "dE_dt_W"]


def parse_float(row, key, number):
    try:
        value = float(row[key])
    except (ValueError, KeyError, TypeError) as exc:
        raise ValueError(f"Row {number}: invalid {key}") from exc
    if not -float("inf") < value < float("inf"):
        raise ValueError(f"Row {number}: {key} is not finite")
    return value


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("input_csv", type=Path)
    parser.add_argument("--output", type=Path, help="New output CSV; existing files are never overwritten")
    args = parser.parse_args()
    target = args.output or args.input_csv.with_name(args.input_csv.stem + "_balance.csv")
    if target.resolve() == args.input_csv.resolve() or target.exists():
        parser.error("Output is the input or already exists: no overwriting")
    with args.input_csv.open(encoding="utf-8-sig", newline="") as source:
        reader = csv.DictReader(source)
        if not reader.fieldnames or any(key not in reader.fieldnames for key in REQUIRED):
            parser.error("Required columns: " + ", ".join(REQUIRED))
        entries = []
        for number, row in enumerate(reader, 2):
            data = {key: parse_float(row, key, number) for key in REQUIRED[1:]}
            q_fluid = data["m_dot_kg_s"] * data["cp_J_kgK"] * (data["T_out_K"] - data["T_in_K"])
            residual = data["q_gen_W"] - q_fluid - data["q_ext_W"] - data["dE_dt_W"]
            denom = max(abs(data["q_gen_W"]), 1e-12)
            entries.append(dict(row, q_fluid_W=f"{q_fluid:.12g}", residual_W=f"{residual:.12g}", residual_percent_source=f"{100*residual/denom:.9g}"))
    with target.open("x", encoding="utf-8", newline="") as out:
        writer = csv.DictWriter(out, fieldnames=reader.fieldnames + ["q_fluid_W", "residual_W", "residual_percent_source"])
        writer.writeheader()
        writer.writerows(entries)
    print(f"Wrote {len(entries)} rows to {target}")
    print("Convention: q_gen - q_fluid - q_ext - dE/dt = residual; positive q_ext is heat leaving the system.")
    print("This arithmetic does not prove that the CFD measurements are correct.")
    return 0


if __name__ == "__main__":
    sys.exit(main())

• Colonne obbligatorie: case, time_s, m_dot_kg_s, T_in_K, T_out_K, cp_J_kgK, q_gen_W, q_ext_W, dE_dt_W. Per ogni riga calcola q_fluid = ṁ c_p(T_out − T_in) e il residuo q_gen − q_fluid − q_ext − dE/dt, anche in percentuale della sorgente.
• Senza --output scrive accanto all’ingresso un file con suffisso _balance.csv, rifiutando di sovrascriverlo.
• Convenzione stampata a fine esecuzione: q_ext positivo è calore che esce dal sistema.
• Eseguito sul CSV sintetico del pacchetto produce un file identico a example_measured_values_SYNTHETIC_balance.csv. Il residuo di 0,02 W della riga reference deriva dall’arrotondamento di T_out nel file di esempio, non dal calcolo.

esercitazione/measurement_template.csv

case,time_s,m_dot_kg_s,T_in_K,T_out_K,cp_J_kgK,q_gen_W,q_ext_W,dE_dt_W

esercitazione/example_measured_values_SYNTHETIC.csv (valori inventati)

case,time_s,m_dot_kg_s,T_in_K,T_out_K,cp_J_kgK,q_gen_W,q_ext_W,dE_dt_W
reference,10,0.004,300,309.565893326,4181,1000,100,740
lowConductivity,10,0.004,300,305.979430758,4181,1000,100,800
lowHeatSource,10,0.004,300,305.979430758,4181,500,50,350

esercitazione/example_measured_values_SYNTHETIC_balance.csv (uscita attesa)

case,time_s,m_dot_kg_s,T_in_K,T_out_K,cp_J_kgK,q_gen_W,q_ext_W,dE_dt_W,q_fluid_W,residual_W,residual_percent_source
reference,10,0.004,300,309.565893326,4181,1000,100,740,159.979999984,0.0200000159764,0.0020000016
lowConductivity,10,0.004,300,305.979430758,4181,1000,100,800,99.9999999968,3.20835624734e-09,3.20835625e-10
lowHeatSource,10,0.004,300,305.979430758,4181,500,50,350,99.9999999968,3.20835624734e-09,6.41671249e-10

esercitazione/README_IT_EN.md

# Guida 11 — Materiale di esercitazione / Exercise kit

**IT** — Questa cartella contiene tre strumenti: `prepare_cases.py` (copia il tutorial ufficiale e prepara tre casi indipendenti), `mesh_case.sh` (costruisce e controlla le mesh regionali) e `solve_case.sh` (esegue la soluzione seriale con `foamMultiRun`). `energy_balance.py` elabora esclusivamente misure registrate in un CSV; non legge automaticamente i campi OpenFOAM. I file sorgente e tutti i commenti nel codice sono in inglese. Gli script non sovrascrivono log, simulazioni o file di output già esistenti. Occorre avere OpenFOAM Foundation 14 attivato nel terminale Linux.

**EN** — This directory contains `prepare_cases.py` (copies the upstream tutorial and prepares three independent cases), `mesh_case.sh` (generates and checks the region meshes), and `solve_case.sh` (runs the serial solution using `foamMultiRun`). `energy_balance.py` processes values that **you measured and entered** into CSV; it does not automatically extract fields from OpenFOAM. The scripts do not overwrite prior cases, solver logs or outputs. Source OpenFOAM Foundation 14 in the Linux shell before running the scripts.

```bash
source /opt/OpenFOAM-14/etc/bashrc
python3 scripts/prepare_cases.py --output "$HOME/OpenFOAM/CHT_Guide11"
bash scripts/mesh_case.sh "$HOME/OpenFOAM/CHT_Guide11/reference"
bash scripts/solve_case.sh "$HOME/OpenFOAM/CHT_Guide11/reference"
python3 scripts/energy_balance.py esercitazione/example_measured_values_SYNTHETIC.csv
```

**Warning:** the output CSV generated by the last command is based on **synthetic test values**, not on a completed CFD simulation. The real heater power depends on the *actual* heater region volume and mesh. The transient energy balance must include internal energy storage and external losses where applicable. The included `measurement_template.csv` is intentionally empty and requires measured values.

**Official case:** https://github.com/OpenFOAM/OpenFOAM-14/tree/master/tutorials/multiRegion/CHT/heatedDuct
**Licence for new scripts:** GNU GPL v3 or later. The upstream OpenFOAM case is **not redistributed** in this bundle; it is copied from the user's installed OpenFOAM tutorial directory.

Calcola il bilancio transitorio passo per passo a partire dai file delle function object della sezione 8. Con densità e calori specifici costanti l’energia di una regione è E = ρcV·T̄_vol + costante, quindi la derivata si ottiene dalle temperature medie con differenze centrate. Le capacità C = ρcV si calcolano con i volumi misurati da checkMesh:

cd "$HOME/OpenFOAM/CHT_Guida11/reference"
find postProcessing -name "*.dat"
# example with nominal volumes: fluid 1000*4181*2.4e-4, metal 2700*900*1.35e-4, heater 8000*450*5.25e-5
python3 cht_balance_series.py --qgen 525 \
  --hin  <path>/hFlowIn/.../surfaceFieldValue.dat \
  --hout <path>/hFlowOut/.../surfaceFieldValue.dat \
  --region fluid:<path>/TavgFluid/.../volFieldValue.dat:1003.4 \
  --region metal:<path>/TavgMetal/.../volFieldValue.dat:328.1 \
  --region heater:<path>/TavgHeater/.../volFieldValue.dat:189.0 \
  --every 20 --output ../reference_balance_series.csv

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Transient energy balance of a CHT run from OpenFOAM function-object files.

Inputs are the .dat files written by the monitoring functions of Guide 11:
  --hin / --hout   surfaceFieldValue 'sum' of h weighted by phi on the fluid
                   inlet and outlet [W] (signed: negative on the inlet);
  --region         NAME:FILE:C  volume-averaged temperature of a region and its
                   heat capacity C = rho*c*V [J/K] (V measured with checkMesh).
With constant density and specific heat the stored energy of a region is
E = C * T_avg + const, so dE/dt follows from finite differences in time.
Balance: residual = Q_gen - Q_fluid - Q_ext - dE/dt, Q_fluid = hout + hin.
The first column of each file is time, the selected column (default 2) the
value. A NEW output CSV is written; existing files are never overwritten.
"""
import argparse
import csv
from pathlib import Path
import sys


def read_series(path, column):
    t, v = [], []
    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.replace('(', ' ').replace(')', ' ').split()
            try:
                t.append(float(parts[0]))
                v.append(float(parts[column - 1]))
            except (IndexError, ValueError) as exc:
                raise ValueError(f'{path}:{number}: cannot read column {column}') from exc
    if len(t) < 3:
        raise ValueError(f'{path}: at least three records are needed')
    pairs = sorted(dict(zip(t, v)).items())
    return [p[0] for p in pairs], [p[1] for p in pairs]


def interp(ts, vs, x):
    if x <= ts[0]:
        return vs[0]
    if x >= ts[-1]:
        return vs[-1]
    lo, hi = 0, len(ts) - 1
    while hi - lo > 1:
        mid = (lo + hi) // 2
        if ts[mid] <= x:
            lo = mid
        else:
            hi = mid
    w = (x - ts[lo]) / (ts[hi] - ts[lo])
    return vs[lo] + w * (vs[hi] - vs[lo])


def main():
    p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument('--qgen', type=float, required=True, help='heater power [W] = q * V_heater')
    p.add_argument('--qext', type=float, default=0.0, help='heat leaving through external walls [W]')
    p.add_argument('--hin', type=Path, required=True)
    p.add_argument('--hout', type=Path, required=True)
    p.add_argument('--region', action='append', required=True, help='NAME:FILE:C, repeat for each region')
    p.add_argument('--column', type=int, default=2, help='value column in the .dat files')
    p.add_argument('--every', type=int, default=1, help='keep one output row every N time records')
    p.add_argument('--output', type=Path, required=True)
    a = p.parse_args()
    if a.output.exists():
        p.error('Output already exists; refusing to overwrite it')
    hin = read_series(a.hin, a.column)
    hout = read_series(a.hout, a.column)
    regions = []
    for spec in a.region:
        try:
            name, file_name, cap = spec.rsplit(':', 2)
            regions.append((name, read_series(Path(file_name), a.column), float(cap)))
        except ValueError as exc:
            p.error(f'bad --region {spec}: {exc}')
    times = regions[0][1][0]
    energy = [sum(c * interp(s[0], s[1], t) for _, s, c in regions) for t in times]
    rows = []
    for i in range(1, len(times) - 1, max(1, a.every)):
        t = times[i]
        dedt = (energy[i + 1] - energy[i - 1]) / (times[i + 1] - times[i - 1])
        q_fluid = interp(*hout, t) + interp(*hin, t)
        residual = a.qgen - q_fluid - a.qext - dedt
        rows.append(dict(time_s=f'{t:.6g}', Q_gen_W=f'{a.qgen:.6g}', Q_fluid_W=f'{q_fluid:.6g}',
                         dE_dt_W=f'{dedt:.6g}', residual_W=f'{residual:.6g}',
                         residual_percent=f'{100 * residual / a.qgen:.4g}' if a.qgen else 'nan'))
    with a.output.open('x', newline='', encoding='utf-8') as f:
        w = csv.DictWriter(f, fieldnames=list(rows[0]))
        w.writeheader()
        w.writerows(rows)
    last = rows[-1]
    print(f'{len(rows)} rows written to {a.output}')
    print(f"Last row t = {last['time_s']} s: Q_fluid = {last['Q_fluid_W']} W, dE/dt = {last['dE_dt_W']} W, "
          f"residual = {last['residual_W']} W ({last['residual_percent']} % of Q_gen)")
    print('Central differences need closely spaced records; check the residual for sensitivity to --every.')


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

Provato su serie analitiche costruite in modo che la somma dei termini sia nulla: il residuo calcolato resta sotto 10⁻³ W su 525 W. Sui dati reali un residuo di pochi percento nei primi passi è normale (passi piccoli, differenze centrate su un transitorio rapido); un residuo che non si riduce con passi più fini indica una misura sbagliata o un problema di conservazione all’interfaccia, da controllare con i valori di wallHeatFlux delle due regioni. Sostituisci <path> con i percorsi elencati da find.

In questo tutorial tutte le proprietà sono costanti, la gravità è nulla e la densità non dipende da T: il campo di moto è identico nei tre casi e l’equazione dell’energia è lineare nella differenza T − 300 K, con sorgente proporzionale a q‴. Anche i passi temporali coincidono, perché dipendono solo dalla velocità. Ne segue una previsione esatta: in lowHeatSource ogni incremento di temperatura deve essere esattamente metà di quello di reference, in ogni cella e a ogni istante.

(T_lowHeatSource - 300) / (T_reference - 300) = 0.5

Lo scostamento misurato, per esempio su T_max del riscaldatore e sulle temperature medie delle function object, dà direttamente l’ordine dell’errore delle tolleranze dei solutori lineari. La guida osserva giustamente che le temperature non variano linearmente se il campo convettivo o le proprietà cambiano: qui, per costruzione, non cambiano, e il test sfrutta proprio questo. Se lo scostamento supera qualche millesimo, controlla che le due copie differiscano davvero solo per q (diff -r).

Compila una riga per ogni caso allo stesso istante (t = 20 s) e indica sempre: istante, tolleranze, passo massimo, volumi da checkMesh. Per l’esercizio C aggiungi il rapporto (T − 300)/(T_ref − 300) e per l’esercizio B le temperature medie della parete inferiore e di quella superiore, che misurano l’effetto «aletta».

In caso di divergenza o di risultati incoerenti, crea una copia nuova per ogni tentativo e annota la modifica nel file GUIDE11_CASE.txt (aggiungendo una riga, senza cancellare le precedenti). Le prove tipiche, nell’ordine: maxDeltaT più piccolo; nOuterCorrectors maggiore; correttori non ortogonali nelle regioni con celle distorte; mesh più fine. Cambiare più parametri insieme rende impossibile capire quale abbia agito.

• Volume del riscaldatore misurato e Q̇_gen calcolato con quel volume.
• Conservazione all’interfaccia: Q(fluid_to_metal) + Q(metal_to_fluid) ≈ 0 e Q(metal_to_heater) + Q(heater_to_metal) ≈ 0.
• Distinzione dichiarata fra calore entrato nell’acqua e calore portato fuori dall’acqua.
• Test di linearità dell’esercizio C superato, con lo scostamento riportato.
• Sensibilità al passo temporale (maxDeltaT dimezzato) documentata.
• Limiti del modello dichiarati: gravità nulla, proprietà costanti, nessuna ebollizione, contatto ideale.

• S. V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, 1980, Routledge: discretizzazione ai volumi finiti, media armonica delle conducibilità alle interfacce;
• T. L. Bergman, A. S. Lavine, F. P. Incropera, D. P. DeWitt, Fundamentals of Heat and Mass Transfer, 8th ed., Wiley, 2017, Wiley: resistenze termiche, alette, transitori conduttivi, convezione in condotti;
• R. K. Shah, A. L. London, Laminar Flow Forced Convection in Ducts, Academic Press, 1978, Elsevier: soluzioni di riferimento per condotti a sezione quadrata e lunghezze d’ingresso;
• J. H. Ferziger, M. Perić, R. L. Street, Computational Methods for Fluid Dynamics, 4th ed., Springer, 2020, doi:10.1007/978-3-319-99693-6: accoppiamento segregato e passo temporale.

La cartella $FOAM_TUTORIALS/multiRegion/CHT della versione 14 contiene, oltre a heatedDuct: VoFcoolingCylinder2D, circuitBoardCooling, coolingCylinder2D, coolingSphere, engine2Valve2D, heatExchanger, hotWire, misalignedDuct, multiphaseCoolingCylinder2D, notchedRoller, reverseBurner, shellAndTubeHeatExchanger, thermalBending, vacuumFurnace, wallBoiling. I loro file system/functions e system/fvSolution sono esempi verificati di misure e controlli dell’accoppiamento.

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

• tutorials/multiRegion/CHT/heatedDuct — il caso completo, con STL, dizionari di mesh e Allrun
• applications/solvers/foamMultiRun/foamMultiRun.C — ciclo PIMPLE multiregione
• applications/modules/solid/solid.C — equazione dell’energia nei solidi e maxDeltaT
• applications/modules/fluid/thermophysicalPredictor.C — equazione dell’energia nel fluido
• coupledThermophysicalTransportModels/coupledTemperature — la condizione d’interfaccia
• src/fvModels/general/heatSource/heatSource.H — sorgente con Q oppure q
• surfaceFieldValueTemplates.C — operazioni sum e average con weightField
• src/functionObjects/field/wallHeatFlux/wallHeatFlux.C — segno e potenza integrata
• foamPostProcess.C — opzioni -region e -solver
• etc/caseDicts/functions — le configurazioni richiamate da #includeFunc

Utilità per chi pubblica i frammenti HTML della guida su un proprio sito: sostituisce il segnaposto __MEDIA_BASE_URL__ con l’indirizzo delle immagini in una copia nuova del file, senza mai sovrascrivere.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Replace WordPress media placeholders in a NEW HTML copy."""
import argparse
from pathlib import Path

p = argparse.ArgumentParser(description=__doc__)
p.add_argument('html', type=Path)
p.add_argument('--base-url', required=True, help='Absolute URL of the WordPress upload folder, without final slash')
p.add_argument('--output', type=Path, required=True)
a=p.parse_args()
if not a.base_url.startswith(('https://','http://')):
    p.error('Base URL must start with https:// or http://')
if a.output.resolve()==a.html.resolve() or a.output.exists():
    p.error('Output cannot overwrite an existing file')
s=a.html.read_text(encoding='utf-8')
if '__MEDIA_BASE_URL__' not in s:
    p.error('No media URL placeholder was found')
a.output.write_text(s.replace('__MEDIA_BASE_URL__', a.base_url.rstrip('/')),encoding='utf-8')
print('Created',a.output)

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

L’archivio comprende le tre copie con log, GUIDE11_CASE.txt, file system/functions e postProcessing. Conserva accanto la scheda di misura, i CSV di bilancio e gli stati ParaView. Nessun file viene rimosso.