OpenFOAM: scambio termico coniugato, interfaccia e stabilità

25 Ottobre 2026 di giux CFD e OpenFOAM, in evidenza 33 min di lettura

Approfondimento collegato alla Guida 11 — scambio termico coniugato (CHT) del percorso Imparare OpenFOAM 14. L’articolo è stato costruito interrogando un notebook di ricerca con le fonti citate (articoli e rassegne ad accesso libero, gli estratti del sorgente di OpenFOAM 14 e la Guida 11), verificando poi ogni numero sulle fonti originali e sul repository della versione 14; il modello monodimensionale della sezione 8 è stato eseguito e i risultati sono riportati senza modifiche. Il video in testa riassume lo stesso percorso.

Video di approfondimento: scambio termico coniugato in OpenFOAM 14, interfaccia e stabilità (generato con Gemini Notebook dalle fonti citate).

Di che cosa parla. Questo approfondimento nasce dalla Guida 11 e risponde a una domanda precisa: che cosa succede, numericamente, sull’interfaccia fra un fluido e un solido quando OpenFOAM 14 risolve le regioni una dopo l’altra, e come si verifica che il risultato non dipenda dal modo in cui sono accoppiate. Si parte dal problema coniugato nella letteratura classica, si legge riga per riga la condizione coupledTemperature e il ciclo di foamMultiRun, si ricostruisce la teoria della stabilità dei metodi partizionati (Giles, Verstraete e Scholl, Errera, CHAMP), si confrontano monolitico e partizionato con i dati di Tandis, Cardiff e Ashrafizadeh (2025) e si verifica tutto con un modello monodimensionale eseguibile. Ogni numero è ricavato dalle fonti citate o calcolato con lo script riportato per intero; dove le fonti non dicono qualcosa, lo dichiariamo.

Interfaccia, pesi e convergenza dell'accoppiamento — Interfaccia, pesi e convergenza dell’accoppiamento

1. Il problema coniugato e quando serve davvero

Nella maggior parte dei calcoli di scambio termico la parete è descritta da una condizione al contorno: temperatura imposta, flusso imposto, oppure un coefficiente di scambio h con una temperatura di riferimento (condizione del terzo tipo). Il problema coniugato rinuncia a tutte e tre: conduzione nel solido e convezione nel fluido si risolvono insieme e temperatura e flusso all’interfaccia diventano risultati del calcolo.

Il termine nasce con T. L. Perelman (1961) ed è sviluppato da A. V. Luikov (1974). La rassegna di A. Dorfman e Z. Renner (2009), ad accesso libero, passa in rassegna quasi duecento pubblicazioni e indica con esempi quantitativi quando il coefficiente h costante smette di funzionare:

resistenze confrontabili. Nel caso di una lastra metallica lambita da metallo liquido, il metodo «usuale» con h di parete isoterma sbaglia del 15–20 %; in moto turbolento l’errore massimo riportato è di circa il 7 %. Il parametro che misura il peso relativo delle resistenze è un numero di Biot di parete, Bi* = h*Δ/λw, e per due fluidi separati da una lastra il Biot complessivo è 1/BiΣ = 1/Bi*₁ + 1/Bi*₂;
temperatura di parete variabile lungo il flusso. Il flusso termico dipende dalla storia a monte dello strato limite, che un h locale non contiene;
alette. Con N² = h̄H²/(λwΔ) minore di 2 metodo usuale e coniugato concordano; al crescere di N² l’errore sull’efficienza dell’aletta arriva al 60–70 % e per N² ≥ 1,9 sul lato posteriore compare l’inversione del flusso termico;
transitori. La soluzione quasi stazionaria dà errori rilevanti soprattutto per Fo < 0,5; nell’esempio della rassegna la soluzione analitica con trasformata di Laplace concorda con quella numerica solo fino a Fo < 0,35.

Il tutorial heatedDuct della Guida 11 appartiene alla prima e alla seconda categoria: una parete di alluminio molto conduttiva, riscaldata solo dal fondo, distribuisce il calore sul perimetro di un canale d’acqua in moto laminare. Un h uniforme imposto a priori darebbe una distribuzione di temperatura sbagliata già in partenza.

2. Le condizioni d’interfaccia: dal continuo alla faccia discreta

Per un’interfaccia ideale (contatto perfetto, senza sorgenti superficiali) valgono due condizioni:

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 la faccia d’interfaccia ha una cella P da una parte, a distanza δP, e una cella N dall’altra, a distanza δN. Patankar (1980, capitolo 4, conducibilità all’interfaccia) mostra che la conducibilità di faccia va costruita come resistenze in serie, non come media aritmetica delle conducibilità:

q = (T_N - T_P) / (delta_P/kappa_P + delta_N/kappa_N)
kappa_face = 2 kappa_P kappa_N / (kappa_P + kappa_N)        se delta_P = delta_N

La media aritmetica sovrastima il flusso quando le conducibilità sono molto diverse: con acqua (0,606 W/(m K)) e alluminio (200 W/(m K)) varrebbe circa 100 W/(m K), mentre la media armonica vale 1,21 W/(m K), il doppio dell’acqua, cioè il valore fisicamente corretto per due mezze celle in serie. Un solutore monolitico tratta l’interfaccia come una faccia interna con questa conducibilità; Tandis et al. (2025) la scrivono come ki = kPkN/((ΔP/Δ)kN + (ΔN/Δ)kP), che è la stessa espressione.

3. `coupledTemperature` riga per riga

OpenFOAM 14 risolve le regioni separatamente, quindi l’interfaccia non è una faccia interna ma una coppia di patch mapped (o mappedWall), una per regione, con la condizione coupledTemperature su entrambe. Nel metodo updateCoeffs() di coupledTemperatureFvPatchScalarField.C, per ogni faccia:

Estratto da coupledTemperatureFvPatchScalarField.C, OpenFOAM 14, righe 407–481 (GPL v3, The OpenFOAM Foundation):

tmp<scalarField> kappa;
tmp<scalarField> sumKappaTByDelta;
tmp<scalarField> sumKappaByDelta;
tmp<scalarField> qByKappa;

// q = alpha.this*sumq
getThis(kappa, sumKappaTByDelta, sumKappaByDelta, sumq, qByKappa);

// Add neighbour contributions
{
    tmp<scalarField> sumKappaTByDeltaNbr;
    tmp<scalarField> sumKappaByDeltaNbr;
    tmp<scalarField> qNbr;

    coupledTemperatureNbr.getNbr
    (
        sumKappaTByDeltaNbr,
        sumKappaByDeltaNbr,
        qNbr
    );

    // Include the effect of the optional neighbour insulation layer
    if (coupledTemperatureNbr.h_.valid())
    {
        const_cast<coupledTemperatureFvPatchScalarField&>
        (
            coupledTemperatureNbr
        ).h_->update();

        const scalarField hFactor
        (
            coupledTemperatureNbr.h_()
           /(coupledTemperatureNbr.h_() + sumKappaByDeltaNbr())
        );
        sumKappaTByDeltaNbr.ref() *= hFactor;
        sumKappaByDeltaNbr.ref() *= hFactor;
    }

    tmp<scalarField> sumKappaTByDeltaNbrMapped
    (
        mapper.fromNeighbour(sumKappaTByDeltaNbr)
    );

    tmp<scalarField> sumKappaByDeltaNbrMapped
    (
        mapper.fromNeighbour(sumKappaByDeltaNbr)
    );

    // Include the effect of the optional insulation layer
    if (h_.valid())
    {
        h_->update();

        const scalarField hFactor(h_()/(h_() + sumKappaByDeltaNbrMapped()));
        sumKappaTByDeltaNbrMapped.ref() *= hFactor;
        sumKappaByDeltaNbrMapped.ref() *= hFactor;
    }

    add(sumKappaTByDelta, sumKappaTByDeltaNbrMapped);
    add(sumKappaByDelta, sumKappaByDeltaNbrMapped);

    if (qNbr.valid())
    {
        sumq += mapper.fromNeighbour(qNbr);
    }
}

this->valueFraction() =
    sumKappaByDelta()/(kappa()*patch().deltaCoeffs() + sumKappaByDelta());

this->refValue() = (sumKappaTByDelta() + sumq)/sumKappaByDelta();

this->refGrad() = qByKappa;

mixedFvPatchScalarField::updateCoeffs();

Nel codice a è kappaEff*deltaCoeffs della cella di questa regione (κeff preso dal modello di trasporto termico della regione: per questo la condizione funziona sia nei fluidi sia nei solidi) e a_n lo stesso valore della cella dall’altra parte, portato sulla faccia dal mapper. Una condizione mixed assegna alla faccia Tface = w·refValue + (1 − w)(TP + refGradient/deltaCoeffs). Sostituendo i tre coefficienti, senza sorgenti:

valueFraction  w = a_n / (a + a_n)
refValue         = (a_n*T_n + q_corr) / a_n
refGradient      = q_s / kappa

T_face = (a*T_P + a_n*T_n) / (a + a_n)
q      = a*a_n/(a + a_n) * (T_n - T_P) = (T_n - T_P) / (1/a + 1/a_n)

È esattamente la formula di Patankar: a convergenza la discretizzazione partizionata coincide con quella monolitica. La differenza è nel tempo: Tn è il valore che la regione vicina ha quando la regione corrente viene risolta, cioè quello dell’iterazione esterna precedente oppure, con un solo correttore esterno, del passo temporale precedente. In letteratura questa è una condizione di tipo Robin (mista) su entrambi i lati, con coefficiente uguale alla conduttanza della cella vicina.

Le opzioni, lette dal sorgente:

Voce	Significato	Note dal codice
`Tnbr`	nome del campo di temperatura nella regione vicina (predefinito `T`)	se la patch vicina non è `coupledTemperature` il calcolo si ferma con un errore fatale che indica il tipo trovato
`h`	resistenza di uno strato sottile non meshato	funzione `wallLayersHeatTransferCoefficient` con `thicknessLayers` e `kappaLayers`; il codice moltiplica la conduttanza vicina per h/(h + a<sub>n</sub>), cioè aggiunge 1/h in serie, con h = 1/Σ(t<sub>i</sub>/κ<sub>i</sub>)
`qs`, `Qs`	flusso sorgente sull’interfaccia (W/m²) o potenza totale (W)	`Qs` è diviso per l’area della patch; specificarli entrambi è un errore fatale
`qr`, `qrNbr`, `qrRelaxation`	flusso radiativo di questa regione o della vicina	sommato al bilancio con rilassamento (predefinito 1); usato da `vacuumFurnace` e `reverseBurner`
nomi storici	`compressible::turbulentTemperatureCoupledBaffleMixed`, `compressible::turbulentTemperatureRadCoupledMixed`	registrati come alias e letti come `coupledTemperature`

Esempio di strato di contatto, nello stile del tutorial thermalBending (dove è presente ma commentato):

interface_to_metal
{
    type            coupledTemperature;
    h
    {
        type            wallLayersHeatTransferCoefficient;
        thicknessLayers (1e-4);       // 0,1 mm di pasta termica
        kappaLayers     (1.0);        // W/(m K)
    }
    value           $internalField;
}

Con questi valori h = 1/(10⁻⁴/1) = 10⁴ W/(m² K); una cella d’alluminio da 1 mm dalla faccia ha a = 200/0,001 = 2·10⁵ W/(m² K): lo strato domina la resistenza del lato metallo di un fattore 20.

Un dettaglio che si legge solo nel codice: in updateCoeffs() la conduttanza della cella vicina è moltiplicata sia per il fattore dello strato h della patch vicina sia per quello della patch corrente. Le resistenze definite sulle due patch si sommano quindi in serie: un solo strato fisico (una pasta termica, un gioco) va definito su una sola delle due patch, altrimenti la sua resistenza conta due volte. Il tutorial thermalBending prevede, commentato, uno strato per lato (1 µm sul dissipatore e 10 µm sulla piastra, entrambi con κ = 0,025): attivandoli, la resistenza totale è quella di 11 µm.

Due osservazioni pignole sul sorgente della versione 14. L’intestazione del file .H descrive ancora voci thicknessLayers e kappaLayers direttamente nella patch, mentre il codice le legge solo dentro il sottodizionario h. E qs è letto con dimensioni dichiarate dimPower/dimTime invece di W/m², per cui conviene scriverlo come numero puro, senza unità esplicite; la tabella «Usage» della documentazione indica correttamente W/m².

4. I pesi nel caso `heatedDuct`

Con le proprietà del tutorial (acqua κ = μcp/Pr = 959·10⁻⁶ × 4181/6,62 ≈ 0,606 W/(m K), alluminio 200, riscaldatore 80) e celle alla stessa distanza dalla faccia, i pesi sono quelli del pannello (b) della figura:

Interfaccia	Regione	valueFraction	Che cosa riceve
fluido–metallo	`fluid`	200/200,6 ≈ 0,997	praticamente la temperatura del metallo (quasi Dirichlet)
fluido–metallo	`metal`	0,606/200,6 ≈ 0,003	praticamente il flusso dal fluido (quasi Neumann)
metallo–riscaldatore	`heater`	200/280 ≈ 0,714	condizione mista
metallo–riscaldatore	`metal`	80/280 ≈ 0,286	condizione mista

Sull’interfaccia fluido-metallo la condizione si riduce da sola allo schema che la teoria (sezione 6) raccomanda: Dirichlet al fluido, Neumann al solido. Sull’interfaccia fra due solidi di conducibilità confrontabile i pesi sono intermedi e l’accoppiamento è più forte. Con distanze diverse il peso segue il rapporto an/a = (κn/δn)/(κ/δ): una cella fluida molto schiacciata contro la parete (strati prismatici) aumenta a e sposta il peso verso il fluido, senza cambiare la natura del risultato finché il rapporto delle conducibilità resta di due ordini di grandezza.

5. Il ciclo di `foamMultiRun`

foamMultiRun carica un modulo per regione (voce regionSolvers di controlDict: in heatedDuct fluid fluid; metal solid; heater solid;) e li fa avanzare insieme. In foamMultiRun.C, a ogni passo temporale e per ogni correttore esterno, ogni fase è eseguita su tutte le regioni prima di passare alla successiva:

per ogni correttore esterno (nOuterCorrectors):
    moveMesh, motionCorrector             regioni con moto
    fvModels().correct()
    prePredictor
    momentumTransportPredictor            fluidi
    thermophysicalTransportPredictor
    momentumPredictor                     fluidi
    while (pimple.correctEnergy())        nEnergyCorrectors, predefinito 1
        thermophysicalPredictor           energia di TUTTE le regioni, in ordine
    pressureCorrector                     fluidi
    momentumTransportCorrector, thermophysicalTransportCorrector

Le leve disponibili sono quindi tre, con ruoli diversi:

nEnergyCorrectors, nel PIMPLE di primo livello (system/fvSolution, letto da pimpleMultiRegionControl.C): ripete solo le equazioni dell’energia di tutte le regioni, scambiando a ogni giro le temperature d’interfaccia, senza risolvere di nuovo quantità di moto e pressione. È la leva più economica quando il moto non dipende dalla temperatura; VoFcoolingCylinder2D usa 2;
nOuterCorrectors, nello stesso blocco: ripete tutta la sequenza; serve quando la temperatura agisce sul moto (galleggiamento, proprietà variabili). engine2Valve2D e multiphaseCoolingCylinder2D usano 3;
outerCorrectorResidualControl, nel PIMPLE di ciascuna regione: interrompe i correttori esterni quando i residui iniziali scendono sotto tolerance e relTol (entrambi obbligatori: altrimenti errore fatale). Si usa con nOuterCorrectors maggiore di 1 come limite superiore. Per i calcoli stazionari residualControl arresta il calcolo: circuitBoardCooling lo usa con e 1e-4 nel solido e U, h 1e-4 nel fluido.

Il passo temporale è il minimo dei maxDeltaT() delle regioni (setDeltaT.C): nel fluido basicFluidSolver.C combina maxDeltaT e maxCo; nel solido solid.C usa solo maxDeltaT di controlDict e i limiti degli fvModels, senza limite sul numero di diffusione. Il passo può crescere al più di solver::deltaTFactor volte da un passo al successivo. In heatedDuct il blocco PIMPLE di primo livello è vuoto: un correttore esterno e un correttore dell’energia, cioè scambio d’interfaccia ritardato di un passo.

6. La teoria della stabilità: Giles (1997)

M. B. Giles ha analizzato con il metodo dei modi normali un problema modello monodimensionale: due materiali con coefficienti costanti, interfaccia in x = 0, griglie uniformi Δx₋ e Δx₊, discretizzazione ai volumi finiti. I parametri sono il numero di Fourier di cella di ciascun lato, d = κΔt/(cΔx²), con c capacità termica volumica, e il rapporto delle capacità delle celle adiacenti, r = c₊Δx₊/(c₋Δx₋). Il lato − (solido) riceve il flusso dall’altro lato, il lato + (fluido) riceve la temperatura: per questa scelta r = cfΔxf/(csΔxs).

Schema	Condizione di stabilità	Verifica numerica nell’articolo
accoppiato in un unico sistema, esplicito	d₋, d₊ ≤ 1/2 (principio del massimo)	—
accoppiato in un unico sistema, implicito	incondizionatamente stabile	—
partizionato esplicito	r < √(1 − 2d₋)/(1 − √(1 − 2d₊)) (eq. 4.17)	d₋ = d₊ = 3/8: limite r < 1; r = 0,99 stabile, r = 1,01 instabile
ibrido (solido implicito, fluido esplicito)	r < √(1 + 2d₋)/(1 − √(1 − 2d₊)) (eq. 4.30)	d₋ = 4, d₊ = 3/8: limite r < 6; r = 5,95 stabile, r = 6,05 instabile
partizionato implicito, dati d’interfaccia ritardati di un passo	per d₋, d₊ ≫ 1: r < √(d₋/d₊) (eq. 4.44)	d₋ = d₊ = 50: limite r < 1; r = 1 stabile, r = 1,1 instabile

L’instabilità, quando c’è, è un modo a dente di sega all’interfaccia. La conclusione dell’articolo è netta: la chiave della stabilità è dare al calcolo strutturale una condizione di Neumann e al fluido una condizione di Dirichlet; lo schema inverso è instabile salvo passi molto piccoli. Per gas e metalli r è piccolissimo, perché la capacità volumica del gas è tre ordini di grandezza inferiore.

Il caso partizionato implicito è quello di foamMultiRun con un solo correttore esterno. Applicato all’acqua e all’alluminio con celle uguali: r = 4,18·10⁶/2,43·10⁶ ≈ 1,72 e √(ds/df) = √(αs/αf) = √(8,23·10⁻⁵/1,449·10⁻⁷) ≈ 23,8. Lo schema raccomandato è stabile con ampio margine. Lo schema inverso (Neumann all’acqua, Dirichlet all’alluminio) scambia i ruoli: r′ = 2,43/4,18 ≈ 0,58 contro un limite √(αf/αs) ≈ 0,042, quindi instabile. La sezione 8 verifica entrambe le previsioni.

Giles non tratta le condizioni miste: considera solo Dirichlet e Neumann. Il passo successivo della letteratura è proprio questo.

7. Dopo Giles: Biot, condizioni di Robin ottimali, CHAMP

T. Verstraete, S. Scholl (2016) confrontano tre metodi partizionati consolidati, ne propongono un quarto e ricavano un criterio di stabilità in forma chiusa in cui il parametro dominante è il numero di Biot: ogni metodo ha un proprio campo di applicabilità e il nuovo è più stabile soprattutto ad alto Biot.
M.-P. Errera, S. Chemin (2013) ricavano, sempre con i modi normali, i coefficienti ottimali delle condizioni di Robin all’interfaccia; il filone prosegue con Errera e Duchaine (2016) sulla procedura Dirichlet-Robin e con Moretti, Errera, Couaillier e Feyel (2018) sull’accoppiamento debole e forte.
F. Meng, J. W. Banks, W. D. Henshaw, D. W. Schwendeman (2017), schema CHAMP: condizione di Robin generalizzata con pesi scelti ottimizzando la stabilità; stabile e del secondo ordine con una sola soluzione implicita per dominio e per passo nella maggior parte dei casi, con sotto-iterazioni solo in casi estremi (griglie molto fini e passi grandi). Gli autori riportano che supera Dirichlet-Neumann quando conducibilità e diffusività sono diverse.

La condizione di OpenFOAM è una Robin-Robin con un coefficiente specifico, la conduttanza della cella vicina: non è ottimizzata nel senso di Errera o di CHAMP, ma ha due proprietà utili. A convergenza coincide con la discretizzazione monolitica di Patankar; e, quando le conducibilità sono molto diverse, degenera automaticamente nella scelta Dirichlet-fluido/Neumann-solido che Giles raccomanda.

8. Una verifica eseguibile: il modello monodimensionale

Per vedere tutto questo in numeri, lo script seguente risolve la conduzione transitoria in due lastre adiacenti (5 mm ciascuna, 20 celle per lato, Eulero implicito), con temperatura imposta alle estremità esterne e un salto iniziale 300/400 K. Confronta la soluzione monolitica con tre accoppiamenti partizionati risolti in sequenza (prima A, poi B) con nOuter correttori per passo: la condizione mista di coupledTemperature, Dirichlet ad A e Neumann a B (DN), e il contrario (ND). Richiede solo numpy e non scrive file.

#!/usr/bin/env python3
"""1D two-material transient conduction: monolithic versus partitioned coupling.

Region A (x < 0) and region B (x > 0), uniform cells, implicit Euler in each region.
Interface treatments:
  mono   : single linear system, face conductance 1/(1/a_A + 1/a_B)  (Patankar)
  mixed  : OpenFOAM coupledTemperature style; each region sees
           T_face = (a_self*T_P + a_nbr*T_nbr)/(a_self + a_nbr) with T_nbr lagged
  DN     : Dirichlet to A (face T from B), Neumann to B (flux from A)
  ND     : the reverse: Neumann to A, Dirichlet to B
With many outer correctors every partitioned scheme that converges reproduces the
monolithic solution to round-off (check with --outer 200).
Regions are solved one after the other (A then B) inside nOuter correctors per step.
Prints the maximum deviation from the monolithic solution and the interface flux.
Only numpy is required. Nothing is written to disk.
"""
import argparse
import numpy as np


def tridiag(lo, di, up, rhs):
    n = len(di); c = np.zeros(n); d = np.zeros(n)
    c[0] = up[0] / di[0]; d[0] = rhs[0] / di[0]
    for i in range(1, n):
        m = di[i] - lo[i] * c[i - 1]
        c[i] = up[i] / m if i < n - 1 else 0.0
        d[i] = (rhs[i] - lo[i] * d[i - 1]) / m
    x = np.zeros(n); x[-1] = d[-1]
    for i in range(n - 2, -1, -1):
        x[i] = d[i] - c[i] * x[i + 1]
    return x


def region_system(T_old, k, rc, dx, dt, left, right):
    """Implicit Euler for one region. left/right = ('T', value) Dirichlet at face,
    ('q', flux into region) Neumann, or ('R', a_ext, T_ext) Robin (series conductance)."""
    n = len(T_old); a = k / dx; cap = rc * dx / dt
    lo = np.zeros(n); di = np.full(n, cap); up = np.zeros(n); rhs = cap * T_old.copy()
    for i in range(n):
        if i > 0:
            di[i] += a; lo[i] = -a
        if i < n - 1:
            di[i] += a; up[i] = -a
    for side, bc, idx in (('L', left, 0), ('R', right, n - 1)):
        ab = 2 * k / dx  # centre-to-face conductance
        if bc[0] == 'T':
            di[idx] += ab; rhs[idx] += ab * bc[1]
        elif bc[0] == 'q':
            rhs[idx] += bc[1]
        elif bc[0] == 'R':
            g = 1.0 / (1.0 / ab + 1.0 / bc[1])
            di[idx] += g; rhs[idx] += g * bc[2]
    return tridiag(lo, di, up, rhs)


def run(kA, rcA, kB, rcB, nA, nB, L, dt, steps, scheme, nOuter, T0A=300.0, T0B=400.0, TendA=300.0, TendB=400.0):
    dxA = L / nA; dxB = L / nB
    TA = np.full(nA, T0A); TB = np.full(nB, T0B)
    aA = 2 * kA / dxA; aB = 2 * kB / dxB
    hist = []
    g0 = 1.0 / (1.0 / aA + 1.0 / aB)
    qA = g0 * (T0B - T0A); qB = -qA  # initial interface flux guess (into A, into B)
    for s in range(steps):
        TA_old, TB_old = TA.copy(), TB.copy()
        if scheme == 'mono':
            n = nA + nB; T_old = np.concatenate([TA_old, TB_old])
            k = np.r_[np.full(nA, kA), np.full(nB, kB)]; dx = np.r_[np.full(nA, dxA), np.full(nB, dxB)]
            rc = np.r_[np.full(nA, rcA), np.full(nB, rcB)]
            cap = rc * dx / dt; lo = np.zeros(n); di = cap.copy(); up = np.zeros(n); rhs = cap * T_old
            for i in range(n - 1):
                g = 1.0 / (dx[i] / (2 * k[i]) + dx[i + 1] / (2 * k[i + 1]))
                di[i] += g; up[i] = -g; di[i + 1] += g; lo[i + 1] = -g
            di[0] += aA; rhs[0] += aA * TendA; di[-1] += aB; rhs[-1] += aB * TendB
            T = tridiag(lo, di, up, rhs); TA, TB = T[:nA], T[nA:]
        else:
            for _ in range(nOuter):
                if scheme == 'mixed':
                    TA = region_system(TA_old, kA, rcA, dxA, dt, ('T', TendA), ('R', aB, TB[0]))
                    TB = region_system(TB_old, kB, rcB, dxB, dt, ('R', aA, TA[-1]), ('T', TendB))
                elif scheme == 'DN':
                    # A gets the face temperature of B, B gets the flux computed by A
                    Tf = TB[0] + qB / aB          # face value seen from B (lagged)
                    TA = region_system(TA_old, kA, rcA, dxA, dt, ('T', TendA), ('T', Tf))
                    qA = aA * (Tf - TA[-1])       # flux from the face into A
                    qB = -qA                      # flux from the face into B
                    TB = region_system(TB_old, kB, rcB, dxB, dt, ('q', qB), ('T', TendB))
                elif scheme == 'ND':
                    # A gets the flux computed by B, B gets the face temperature of A
                    TA = region_system(TA_old, kA, rcA, dxA, dt, ('T', TendA), ('q', qA))
                    Tf = TA[-1] + qA / aA         # face value seen from A
                    TB = region_system(TB_old, kB, rcB, dxB, dt, ('T', Tf), ('T', TendB))
                    qB = aB * (Tf - TB[0])        # flux from the face into B
                    qA = -qB
        g = 1.0 / (1.0 / aA + 1.0 / aB)
        hist.append((np.concatenate([TA, TB]), g * (TB[0] - TA[-1])))
    return hist


def main():
    np.seterr(all='ignore')  # the unstable schemes overflow on purpose
    p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument('--pair', choices=['water-aluminium', 'steel-aluminium'], default='water-aluminium',
                   help='region A - region B (A is the first, lower-effusivity material)')
    p.add_argument('--cells', type=int, default=20, help='cells per region')
    p.add_argument('--L', type=float, default=5e-3, help='thickness of each region [m]')
    p.add_argument('--steps', type=int, default=40, help='time steps')
    p.add_argument('--dt', type=float, nargs='+', default=[1e-3, 1e-2, 1e-1, 1.0])
    p.add_argument('--outer', type=int, nargs='+', default=[1, 3, 10])
    a = p.parse_args()
    water = (0.6057, 1000 * 4181); alu = (200.0, 2700 * 900); steel = (80.0, 8000 * 450)
    A, B = (water, alu) if a.pair == 'water-aluminium' else (steel, alu)
    dx = a.L / a.cells
    eA = (A[0] * A[1]) ** 0.5; eB = (B[0] * B[1]) ** 0.5
    print(f'pair {a.pair}: kA={A[0]} rhocA={A[1]:.4g}  kB={B[0]} rhocB={B[1]:.4g}  '
          f'{a.cells}+{a.cells} cells, dx={dx:.3e} m, {a.steps} steps, initial jump 300/400 K')
    print(f'effusivities A={eA:.0f} B={eB:.0f} W s^0.5/(m2 K), sigma={min(eA, eB) / max(eA, eB):.3f}, '
          f'valueFraction seen by A = {B[0] / (A[0] + B[0]):.4f}, by B = {A[0] / (A[0] + B[0]):.4f}')
    print('max|T - T_mono| over all cells and steps [K]; in brackets the value at the last step')
    for dt in a.dt:
        dA = A[0] / A[1] * dt / dx ** 2; dB = B[0] / B[1] * dt / dx ** 2
        ref = run(A[0], A[1], B[0], B[1], a.cells, a.cells, a.L, dt, a.steps, 'mono', 1)
        print(f'dt={dt:g} s  Fourier dA={dA:.3g} dB={dB:.3g}')
        for sch in ('mixed', 'DN', 'ND'):
            line = f'   {sch:5s}'
            for no in a.outer:
                h = run(A[0], A[1], B[0], B[1], a.cells, a.cells, a.L, dt, a.steps, sch, no)
                e = [np.max(np.abs(x[0] - r[0])) for x, r in zip(h, ref)]
                emax = max(e)
                if not np.isfinite(emax) or emax > 1e4:
                    line += f'  nOuter={no:<3d} unstable          '
                else:
                    line += f'  nOuter={no:<3d} {emax:9.3g} [{e[-1]:8.2g}]'
            print(line)


if __name__ == '__main__':
    main()

Coppia acqua (A)–alluminio (B), σ ≈ 0,072:

pair water-aluminium: kA=0.6057 rhocA=4.181e+06  kB=200.0 rhocB=2.43e+06  20+20 cells, dx=2.500e-04 m, 40 steps, initial jump 300/400 K
effusivities A=1591 B=22045 W s^0.5/(m2 K), sigma=0.072, valueFraction seen by A = 0.9970, by B = 0.0030
max|T - T_mono| over all cells and steps [K]; in brackets the value at the last step
dt=0.001 s  Fourier dA=0.00232 dB=1.32
   mixed  nOuter=1      0.0154 [   0.015]  nOuter=3    6.93e-12 [ 6.9e-12]  nOuter=10   5.17e-12 [ 5.2e-12]
   DN     nOuter=1      0.0146 [   0.015]  nOuter=3    8.37e-07 [ 8.4e-07]  nOuter=10    5.4e-12 [ 5.4e-12]
   ND     nOuter=1   unstable            nOuter=3   unstable            nOuter=10  unstable          
dt=0.01 s  Fourier dA=0.0232 dB=13.2
   mixed  nOuter=1       0.158 [    0.03]  nOuter=3    1.04e-07 [   2e-08]  nOuter=10   9.09e-13 [ 6.8e-13]
   DN     nOuter=1       0.148 [   0.025]  nOuter=3    6.77e-05 [ 1.2e-05]  nOuter=10   8.53e-13 [ 8.5e-13]
   ND     nOuter=1   unstable            nOuter=3   unstable            nOuter=10  unstable          
dt=0.1 s  Fourier dA=0.232 dB=132
   mixed  nOuter=1         1.2 [   0.017]  nOuter=3     0.00032 [ 4.5e-06]  nOuter=10   1.71e-12 [ 1.5e-12]
   DN     nOuter=1        1.17 [   0.016]  nOuter=3     0.00261 [ 3.7e-05]  nOuter=10   1.48e-12 [ 6.8e-13]
   ND     nOuter=1   unstable            nOuter=3   unstable            nOuter=10  unstable          
dt=1 s  Fourier dA=2.32 dB=1.32e+03
   mixed  nOuter=1        2.22 [  0.0046]  nOuter=3     0.00978 [ 1.9e-05]  nOuter=10   5.46e-11 [ 1.2e-12]
   DN     nOuter=1        2.08 [  0.0049]  nOuter=3     0.00246 [   6e-06]  nOuter=10   1.14e-12 [ 6.8e-13]
   ND     nOuter=1   unstable            nOuter=3   unstable            nOuter=10  unstable

Coppia acciaio del riscaldatore (A)–alluminio (B), σ ≈ 0,770:

pair steel-aluminium: kA=80.0 rhocA=3.6e+06  kB=200.0 rhocB=2.43e+06  20+20 cells, dx=2.500e-04 m, 40 steps, initial jump 300/400 K
effusivities A=16971 B=22045 W s^0.5/(m2 K), sigma=0.770, valueFraction seen by A = 0.7143, by B = 0.2857
max|T - T_mono| over all cells and steps [K]; in brackets the value at the last step
dt=0.001 s  Fourier dA=0.356 dB=1.32
   mixed  nOuter=1        6.89 [     1.5]  nOuter=3       0.051 [  0.0098]  nOuter=10    1.7e-09 [ 3.3e-10]
   DN     nOuter=1        4.09 [   0.052]  nOuter=3        1.69 [   0.028]  nOuter=10     0.0773 [  0.0016]
   ND     nOuter=1   unstable            nOuter=3   unstable            nOuter=10  unstable          
dt=0.01 s  Fourier dA=3.56 dB=13.2
   mixed  nOuter=1          25 [     1.7]  nOuter=3        4.71 [    0.16]  nOuter=10     0.0137 [  0.0004]
   DN     nOuter=1        10.7 [    0.11]  nOuter=3        6.01 [   0.075]  nOuter=10       0.81 [   0.015]
   ND     nOuter=1   unstable            nOuter=3   unstable            nOuter=10  unstable          
dt=0.1 s  Fourier dA=35.6 dB=132
   mixed  nOuter=1        35.5 [    0.35]  nOuter=3        20.4 [ 0.00059]  nOuter=10       2.93 [   2e-05]
   DN     nOuter=1        12.3 [  0.0003]  nOuter=3        6.44 [ 3.5e-05]  nOuter=10      0.671 [ 4.3e-06]
   ND     nOuter=1   unstable            nOuter=3   unstable            nOuter=10  unstable          
dt=1 s  Fourier dA=356 dB=1.32e+03
   mixed  nOuter=1        30.5 [    0.25]  nOuter=3          23 [ 2.5e-05]  nOuter=10        8.5 [ 2.4e-12]
   DN     nOuter=1         5.3 [ 7.3e-09]  nOuter=3        1.02 [   1e-12]  nOuter=10    0.00698 [ 1.1e-12]
   ND     nOuter=1   unstable            nOuter=3   unstable            nOuter=10  unstable

Controllo di coerenza con molti correttori (Δt = 0,1 s):

pair water-aluminium: kA=0.6057 rhocA=4.181e+06  kB=200.0 rhocB=2.43e+06  20+20 cells, dx=2.500e-04 m, 40 steps, initial jump 300/400 K
effusivities A=1591 B=22045 W s^0.5/(m2 K), sigma=0.072, valueFraction seen by A = 0.9970, by B = 0.0030
max|T - T_mono| over all cells and steps [K]; in brackets the value at the last step
dt=0.1 s  Fourier dA=0.232 dB=132
   mixed  nOuter=50   1.71e-12 [ 1.5e-12]  nOuter=200  1.71e-12 [ 1.5e-12]
   DN     nOuter=50   8.53e-13 [ 6.8e-13]  nOuter=200  8.53e-13 [ 6.8e-13]
   ND     nOuter=50  unstable            nOuter=200 unstable          
pair steel-aluminium: kA=80.0 rhocA=3.6e+06  kB=200.0 rhocB=2.43e+06  20+20 cells, dx=2.500e-04 m, 40 steps, initial jump 300/400 K
effusivities A=16971 B=22045 W s^0.5/(m2 K), sigma=0.770, valueFraction seen by A = 0.7143, by B = 0.2857
max|T - T_mono| over all cells and steps [K]; in brackets the value at the last step
dt=0.1 s  Fourier dA=35.6 dB=132
   mixed  nOuter=50   4.45e-05 [ 2.4e-10]  nOuter=200  1.25e-12 [ 9.1e-13]
   DN     nOuter=50   1.64e-06 [ 9.4e-12]  nOuter=200  1.93e-12 [   8e-13]
   ND     nOuter=50  unstable            nOuter=200 unstable

Che cosa mostrano i numeri:

Tutti gli schemi che convergono riproducono il monolitico all’errore di arrotondamento (10⁻¹² K con 200 correttori): l’accoppiamento partizionato cambia il percorso, non la soluzione.
Lo schema ND è instabile in tutti i casi, come prevede Giles: per acqua-alluminio r′ ≈ 0,58 contro un limite ≈ 0,042; per acciaio-alluminio, con il lato acciaio che riceve il flusso, r′ = 2,43/3,60 ≈ 0,68 contro √(dacciaio/dalluminio) = √(0,356/1,32) ≈ 0,52.
Con σ piccolo (acqua-alluminio) la condizione mista e DN sono quasi indistinguibili, perché w ≈ 0,997 e 0,003 le rendono di fatto lo stesso schema; con un solo correttore lo scarto massimo, dovuto ai primi passi dopo il salto iniziale, va da 0,015 K (Δt = 1 ms) a 2,2 K (Δt = 1 s), e all’ultimo passo resta entro 0,03 K. Tre correttori bastano per arrivare a 10⁻² K anche con Δt = 1 s.
Con σ ≈ 0,77 (acciaio-alluminio) l’accoppiamento è forte: con la condizione mista e un solo correttore lo scarto massimo arriva a 25–35 K su un salto di 100 K per Δt ≥ 0,01 s. Con Δt = 0,01 s la condizione mista converge più rapidamente di DN (fattore ≈ 0,43 per correttore contro ≈ 0,75: pannello (c) della figura), mentre con passi molto grandi (Δt = 1 s, d ≫ 1) è DN a convergere più in fretta. Nessuno dei due schemi è il migliore in assoluto: è la stessa conclusione, con altri parametri, di Verstraete e Scholl e del filone di Errera.

Il modello è monodimensionale e senza convezione: non sostituisce un caso tridimensionale, ma isola esattamente il meccanismo dell’interfaccia e permette di controllare in pochi secondi l’effetto di passo, correttori e proprietà.

9. Monolitico contro partizionato in OpenFOAM: Tandis, Cardiff e Ashrafizadeh (2025)

Il riferimento più vicino a questo articolo è E. Tandis, P. Cardiff, A. Ashrafizadeh, «Analysis of Coupling Strategies for Conjugate Heat Transfer Problems», OpenFOAM® Journal 5 (2025) 38–58, ad accesso libero (CC BY-SA 4.0). Gli autori passano in rassegna i solutori CHT di foam-extend 4.0 (conjugateHeatFoam, monolitico, una regione fluida e una solida, Boussinesq; chtMultiRegionFoam, partizionato), di OpenFOAM-v2112 e di OpenFOAM-9, descrivono la condizione mista dei solutori partizionati con lo stesso peso ω = (kN/ΔN)/(kP/ΔP + kN/ΔN) della sezione 3 e propongono un solutore unico (multiChtFoam) con tre algoritmi: partizionato classico a ciclo integrato (SIPM-IP, come chtMultiRegionFoam), partizionato con ciclo separato per l’interfaccia (SIPM-SP) e monolitico multiregione (SIPM-M).

La forza dell’accoppiamento termico è misurata con il rapporto delle effusività σ = (k₁/√α₁)/(k₂/√α₂), con la regione di effusività maggiore al denominatore, quindi 0 < σ < 1; σ → 1 significa che entrambe le regioni variano molto e l’accoppiamento temperatura-temperatura è forte. Il numero di Grashof misura invece l’accoppiamento pressione-velocità-temperatura. Il caso di prova è una cavità quadrata di 1 m con una parete solida di 0,2 m (80×80 celle nel fluido, 80×20 nel solido, Δt = 2,5·10⁻⁴ s, Pr = 0,7):

Stato	k<sub>s</sub>/k<sub>f</sub>	Gr	σ	Iterazioni esterne IP / SP / M	CPU IP / SP / M [s]
A	1,6	1,43·10⁵	0,93	12,2 / 9,3 / 9,3	40,1 / 32,3 / 30,17
B	1,6	10⁴	0,93	12,7 / 5,6 / 5,6	41,7 / 22,9 / 20,79
C	1600	1,43·10⁵	4·10⁻⁴	10,3 / 9,3 / 8,7	36,5 / 30,6 / 28,8
D	1600	10⁴	4·10⁻⁴	7,5 / 5,5 / 5,4	29,0 / 21,2 / 20,43

Con il predittore della quantità di moto attivo il monolitico costa dal 21 al 50 % in meno del partizionato classico; con tre cavità e tre pareti (stessi stati) il vantaggio cresce: nello stato B 137,2 s per IP contro 69,3 s per M (15,6 contro 6,6 iterazioni esterne). Gli autori concludono che imporre le condizioni d’interfaccia con un ciclo interno dedicato o con una soluzione simultanea può accelerare il calcolo fino a due volte, e che il monolitico è superiore soprattutto con più di due regioni. In OpenFOAM 14 l’equivalente più vicino al ciclo separato è nEnergyCorrectors, che ripete solo l’energia: è la prima cosa da provare quando lo scambio d’interfaccia è il collo di bottiglia.

10. Effusività e forza dell’accoppiamento in `heatedDuct`

Le effusività e = √(κρc) delle tre regioni, già calcolate dallo script cht_estimates.py della Guida 11, sono 1591 (acqua), 22045 (alluminio) e 16971 (riscaldatore) W s0,5/(m² K). Ne segue:

Interfaccia	σ	Accoppiamento	Conseguenza pratica
acqua–alluminio	1591/22045 ≈ 0,072	debole	il ritardo di un passo pesa poco; il controllo principale è la mesh dal lato fluido, dove sta tutta la resistenza
riscaldatore–alluminio	16971/22045 ≈ 0,77	forte	è l’interfaccia più sensibile ai correttori: il bilancio Q(metal_to_heater) + Q(heater_to_metal) ≈ 0 della sezione 13 della guida va controllato per primo

Il quadro è coerente con i pesi della sezione 4 e con il modello 1D della sezione 8: dove σ è piccolo l’interfaccia si comporta come Dirichlet-Neumann e converge in pochi correttori, dove σ è vicino a 1 servono più correttori o passi più piccoli.

11. Valori di riferimento per il condotto quadrato

Per moto laminare completamente sviluppato in un condotto a sezione quadrata, Shah e London (1978) riportano, con Dh come lunghezza di riferimento:

Grandezza	Valore	Condizione termica di parete
f·Re (Fanning)	14,227 (Darcy 56,91)	—
Nu<sub>T</sub>	2,976	temperatura uniforme su perimetro e lunghezza
Nu<sub>H1</sub>	3,608	flusso assiale uniforme, temperatura uniforme sul perimetro (parete molto conduttiva)
Nu<sub>H2</sub>	3,091	flusso uniforme anche sul perimetro (parete senza conduzione circonferenziale)

Gli stessi valori arrotondati (2,98 e 3,61) sono tabulati da Muzychka e Yovanovich (2002, ad accesso libero). La distinzione fra H1 e H2 è un effetto coniugato: una parete di alluminio molto conduttiva avvicina heatedDuct alla condizione H1 anche se il calore entra solo dal fondo; abbassando κ del metallo (esercizio con l’acciaio della Guida 11) ci si sposta verso H2 e il Nusselt medio diminuisce.

Nessuno dei due valori si confronta direttamente con il tutorial: con Re ≈ 42 e Pr = 6,62 la lunghezza d’ingresso termica stimata con 0,05 Re Pr Dh è circa 0,55 m, 3,7 volte il condotto da 0,15 m; a t = 20 s il caso è ancora un transitorio, con il tempo di attraversamento L/U = 150 s e una costante di tempo globale di circa 227 s. Per un confronto quantitativo servirebbe una copia stazionaria con un condotto molto più lungo, potenza ridotta e un Nusselt definito con la temperatura di miscelamento e la temperatura media di parete.

12. Flusso congelato: lo stazionario senza risolvere il moto

Il tutorial vacuumFurnace (riscaldatore da 50 kW, irraggiamento e convezione naturale) accelera il raffreddamento di un blocco d’acciaio con tre strategie descritte nel suo README: frozen (moto calcolato una volta, poi solo energia), periodic (moto ricalcolato a intervalli) e continuous. Il suo Allrun imposta PIMPLE/flow a off in system/fluid/fvSolution con foamDictionary e porta il passo da 0,05 s a 50 s.

In heatedDuct la stessa tecnica è esatta: gravità nulla, densità e viscosità costanti rendono la quantità di moto indipendente dalla temperatura. In una copia nuova, dopo aver portato a regime il campo di moto:

foamDictionary system/fluid/fvSolution -entry PIMPLE/flow -set off
foamDictionary system/controlDict -set "adjustTimeStep=no, deltaT=1, endTime=2000"
foamMultiRun > log.foamMultiRun.frozen 2>&1

adjustTimeStep va disattivato: anche con il moto congelato, basicFluidSolver::maxDeltaT() continua a limitare il passo con maxCo calcolato sul flusso phi esistente. Il risultato utile è lo stato stazionario; la storia temporale ottenuta con passi grandi non è accurata e non va usata come transitorio.

13. Un protocollo di verifica per `heatedDuct`

Ogni prova si esegue in una copia nuova del caso, annotando la modifica, come nella Guida 11.

Riferimento. Il caso del tutorial, con le misure della guida: temperature massime per regione, potenze d’interfaccia (wallHeatFlux), temperatura di miscelamento in uscita, bilancio con accumulo.
Correttori dell’energia. nEnergyCorrectors 3 nel PIMPLE di primo livello. Se le potenze d’interfaccia cambiano oltre la tolleranza, l’accoppiamento termico non era convergente.
Correttori esterni con controllo dei residui. nOuterCorrectors 10 con outerCorrectorResidualControl in ogni regione: il log mostra quanti correttori servono davvero a ogni passo.
Passo. maxDeltaT dimezzato: cambiano le temperature dei solidi? Nel solido non c’è un limite automatico sul numero di Fourier.
Conservazione. Q(fluid_to_metal) + Q(metal_to_fluid) ≈ 0 e Q(metal_to_heater) + Q(heater_to_metal) ≈ 0 a ogni istante salvato, con la somma confrontata con le potenze in gioco: uno scarto sistematico indica un problema di mappatura delle patch o un accoppiamento non convergente.
Mesh dal lato fluido. È lì che sta la resistenza: un raffinamento vicino alla parete pesa sul flusso scambiato più di qualunque altro parametro.

// system/fvSolution (primo livello), copia "ecorr3"
PIMPLE
{
    nEnergyCorrectors  3;
}

// system/fvSolution (primo livello), copia "outer10"
PIMPLE
{
    nOuterCorrectors   10;
}
// e in system/fluid/fvSolution, system/metal/fvSolution, system/heater/fvSolution
PIMPLE
{
    outerCorrectorResidualControl
    {
        "h|e"  { tolerance 1e-6; relTol 0; }
    }
}

14. Gli altri tutorial CHT di OpenFOAM 14

Tutorial	Regioni → moduli	Tempo e correttori	Che cosa mostra
`coolingCylinder2D`	fluid → fluid, solid → solid	Euler, Δt 0,01 s fisso, 1 correttore	il caso minimo: `blockMesh`, `splitMeshRegions -cellZones`, `foamMultiRun`
`VoFcoolingCylinder2D`	fluid → compressibleVoF, solid → solid	maxCo 2, `nEnergyCorrectors 2`	CHT con un fluido bifase a interfaccia libera
`multiphaseCoolingCylinder2D`	fluid → multiphaseEuler, solid → solid	maxCo 0,25, 3 correttori	CHT con il modello euleriano multifase
`wallBoiling`	fluid → multiphaseEuler, solid → solid	maxCo 0,5, maxDeltaT 0,005	ebollizione a parete, cartella `validation`, grafici di diametro e frequenza di distacco
`circuitBoardCooling`	fluid → fluid, baffle3D → solid	stazionario, `residualControl`	scheda come baffle estruso in una regione solida, da facce interne o da patch
`coolingSphere`	da `foamSetupCHT` e `templates`	Euler, Δt 5·10⁻⁴ s	impostazione automatica delle regioni da modelli di materiale (aria, acqua, alluminio, rame)
`heatExchanger`	porous → fluid, air → fluid	stazionario	regioni sovrapposte senza interfaccia: `interRegionHeatTransfer` (Av = 200 m⁻¹, htc = 10) e `interRegionPorosityForce` Darcy-Forchheimer
`shellAndTubeHeatExchanger`	solid → solid, shell → fluid, tube → fluid	stazionario	due fluidi separati da un solido, mesh `snappyHexMesh` con e senza strati
`misalignedDuct`	fluid → fluid, solid → solid	Euler, Δt 0,01 s fisso	interfaccia fra mesh non conformi (`createNonConformalCouples`)
`hotWire`	fluid → fluid, solid → solid	stazionario	filo caldo: solido raffinato e cucito, fluido ruotato di 45° e meshato con `snappyHexMesh`, accoppiamento non conforme
`notchedRoller`	fluid → fluid, solid → solid, roller → solid	Euler, Δt 0,01 s fisso	rullo solido in rotazione (`solidBody`, 20 rpm) con accoppiamento non conforme
`engine2Valve2D`	fluid → fluid; liner, piston, cylinderHead, exhaustValve, intakeValve → solid	Euler, 3 correttori, `correctPhi yes`	motore 2D con pistone e valvole in moto, sei regioni
`reverseBurner`	gas → multicomponentFluid, solid → solid	gas Euler, solido `steadyState`, maxCo 0,2	combustione e irraggiamento (`qrNbr qr`); modalità mista stazionario/transitorio dichiarata nel log
`vacuumFurnace`	solid, heater, shield → solid; fluid → fluid	Euler, Δt cambiato da `Allrun`	irraggiamento, convezione naturale e le tre strategie di flusso congelato
`thermalBending`	heatsink → solid, plate → solidDisplacement	stazionario	sollecitazioni termiche (`thermalStress yes`, `planeStress yes`); strati `h` commentati

15. Che cosa le fonti non dicono

Nessuna analisi di stabilità tratta la condizione esatta di OpenFOAM (Robin-Robin con coefficiente uguale alla conduttanza della cella vicina) con i parametri di un caso tridimensionale con convezione. Giles considera solo Dirichlet e Neumann; Errera, Verstraete e CHAMP studiano coefficienti di Robin ottimizzati. Il modello 1D della sezione 8 è una verifica, non una dimostrazione.
L’ordine delle regioni in regionSolvers determina l’ordine di soluzione nel ciclo dell’energia. Nessuna fonte consultata ne quantifica l’effetto sulla convergenza.
Celle non ortogonali. deltaCoeffs sulle patch accoppiate e le eventuali correzioni non ortogonali alterano i pesi rispetto alla formula a celle allineate; le fonti non ne valutano l’effetto.
Tandis et al. usano fluidi a densità costante con Boussinesq in cavità chiuse: i loro tempi non si trasferiscono automaticamente a fluidi comprimibili, reattivi o multifase.
Le soluzioni di Shah e London sono asintotiche per condotti infiniti a regime: non esiste una soluzione chiusa per il transitorio coniugato tridimensionale di heatedDuct con riscaldamento asimmetrico.

16. Riferimenti

T. L. Perelman, «On conjugated problems of heat transfer», Int. J. Heat Mass Transfer 3 (1961) 293–303, doi:10.1016/0017-9310(61)90044-890044-8).
A. V. Luikov, «Conjugate convective heat transfer problems», Int. J. Heat Mass Transfer 17 (1974) 257–265.
A. Dorfman, Z. Renner, «Conjugate problems in convective heat transfer: review», Math. Probl. Eng. 2009, 927350, doi:10.1155/2009/927350, accesso aperto; A. S. Dorfman, Conjugate Problems in Convective Heat Transfer, CRC Press, 2009.
S. V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, 1980, Routledge.
F. Moukalled, L. Mangani, M. Darwish, The Finite Volume Method in Computational Fluid Dynamics, Springer, 2016, doi:10.1007/978-3-319-16874-6.
J. H. Ferziger, M. Perić, R. L. Street, Computational Methods for Fluid Dynamics, 4ª ed., Springer, 2020, doi:10.1007/978-3-319-99693-6.
M. B. Giles, «Stability analysis of numerical interface boundary conditions in fluid-structure thermal analysis», Int. J. Numer. Methods Fluids 25 (1997) 421–436, DOI1097-0363(19970830)25:4%3C421::AID-FLD557%3E3.0.CO;2-J), PDF dell’autore.
T. Verstraete, S. Scholl, «Stability analysis of partitioned methods for predicting conjugate heat transfer», Int. J. Heat Mass Transfer 101 (2016) 852–869, doi:10.1016/j.ijheatmasstransfer.2016.05.041.
M.-P. Errera, S. Chemin, «Optimal solutions of numerical interface conditions in fluid–structure thermal analysis», J. Comput. Phys. 245 (2013) 431–455, ScienceDirect; M.-P. Errera, F. Duchaine, J. Comput. Phys. 312 (2016) 218–234; R. Moretti, M.-P. Errera, V. Couaillier, F. Feyel, Int. J. Therm. Sci. 126 (2018) 23–37.
F. Meng, J. W. Banks, W. D. Henshaw, D. W. Schwendeman, «A stable and accurate partitioned algorithm for conjugate heat transfer», J. Comput. Phys. 344 (2017) 51–85, OSTI.
E. Tandis, P. Cardiff, A. Ashrafizadeh, «Analysis of Coupling Strategies for Conjugate Heat Transfer Problems», OpenFOAM® Journal 5 (2025) 38–58, doi:10.51560/ofj.v5.92, accesso aperto.
R. K. Shah, A. L. London, Laminar Flow Forced Convection in Ducts, Academic Press, 1978, Elsevier.
Y. S. Muzychka, M. M. Yovanovich, «Laminar flow friction and heat transfer in non-circular ducts and channels, Part II: thermal problem», in Compact Heat Exchangers, Edizioni ETS, 2002, 131–139, accesso aperto.
T. L. Bergman, A. S. Lavine, F. P. Incropera, D. P. DeWitt, Fundamentals of Heat and Mass Transfer, 8ª ed., Wiley, 2017.
OpenFOAM Foundation 14: coupledTemperatureFvPatchScalarField.C, wallLayersHeatTransferCoefficient, foamMultiRun.C, pimpleMultiRegionControl.C, solid.C, basicFluidSolver.C, tutorial multiRegion/CHT.

Per il caso completo, le misure di potenza e temperatura, il bilancio con accumulo e gli script di stima si veda la Guida 11; per lo scambio termico in una sola regione e la convezione naturale la Guida 10.

#CFD #CHT #coupledTemperature #foamMultiRun #metodi partizionati #OpenFOAM #trasferimento di calore #trasferimento termico coniugato