OpenFOAM 14, parte 5: fvSchemes, dalle equazioni alla discretizzazione

Il metodo dei volumi finiti parte dalla forma integrale delle equazioni di conservazione. Si integra l’equazione su ogni cella (volume di controllo) V_P e si usa il teorema di Gauss (della divergenza) per trasformare gli integrali di volume delle divergenze in integrali di superficie, cioè in somme sulle facce della cella. Per la grandezza generica a dell’equazione di trasporto della sezione 2:

∫V ∂a/∂t dV  +  ∮S (U a)·dS  =  ∮S (Γ ∇a)·dS  +  ∫V S dV

(∂a/∂t)_P V_P  +  Σ_f F_f a_f  =  Σ_f Γ_f |S_f| (∇a)_f·n_f  +  S_P V_P

Qui F_f = S_f·U_f è il flusso volumetrico attraverso la faccia f (il campo phi), a_f il valore della grandezza sulla faccia, Γ_f il coefficiente di diffusione interpolato sulla faccia e (∇a)_f·n_f il gradiente normale alla faccia. Ognuno di questi ingredienti richiede una regola di approssimazione: è esattamente ciò che si sceglie in fvSchemes. Quando ogni a_f e ogni gradiente di faccia sono espressi come combinazione dei valori ai centri delle celle vicine, l’equazione di ogni cella diventa una relazione lineare:

a_P · a_P_value  +  Σ_N a_N · a_N_value  =  b_P

con un coefficiente diagonale a_P, un coefficiente a_N per ogni cella vicina che condivide una faccia e un termine noto b_P che raccoglie sorgenti, contributi espliciti e valori al contorno. Mettendo insieme tutte le celle si ottiene una matrice sparsa: in una mesh esaedrica ogni riga ha al massimo 7 elementi non nulli (la cella e le sue 6 vicine). fvSchemes decide come si costruiscono i coefficienti; fvSolution (Parte 6) decide come si risolve il sistema che ne risulta. Nel codice sorgente la distinzione si vede nei due spazi dei nomi: fvm:: (finite volume method) produce contributi impliciti alla matrice, fvc:: (finite volume calculus) valuta un termine esplicitamente con i valori correnti; entrambi leggono lo schema dallo stesso sotto-dizionario di fvSchemes.

L’obiettivo «trovare la documentazione ufficiale v14» ha tre strumenti concreti, tutti disponibili in una shell OpenFOAM 14. Il primo è la sezione 4.5 della User Guide v14, organizzata per sotto-dizionario (4.5.1 schemi temporali, 4.5.2 gradienti, 4.5.3 divergenze, 4.5.4 gradienti normali, 4.5.5 laplaciani, 4.5.6 interpolazione). Il secondo sono i tutorial installati in $FOAM_TUTORIALS, che mostrano come gli sviluppatori usano gli schemi in casi reali: foamSearch estrae da tutti i file con un certo nome il valore di una voce e ne stampa l’elenco ordinato senza duplicati; con l’opzione -count (vedi foamSearch -help) aggiunge quante volte compare ciascun valore. Il terzo è foamToC, che elenca le classi selezionabili a run-time compilate nell’installazione: se un nome di schema non compare, non è disponibile nella tua versione.

# Which scheme names are actually used, and how often, in the v14 tutorials
foamSearch -count "$FOAM_TUTORIALS" fvSchemes "divSchemes/div(phi,U)"

# Which tutorials use a given keyword (plain grep over the installed tree)
grep -rl "linearUpwind" "$FOAM_TUTORIALS"/incompressibleFluid/*/system/fvSchemes

# Which run-time selectable classes contain a name (v14 table-of-contents tool)
foamToC -search linearUpwind

Esiste anche un modo «sperimentale» e del tutto sicuro, da usare solo in una copia del caso: scrivere di proposito un nome inesistente (per esempio div(phi,U) Gauss xyz;). All’avvio OpenFOAM si ferma con un errore del tipo «Unknown discretisation scheme xyz» seguito dall’elenco Valid schemes are: è l’elenco autentico degli schemi compilati per quel tipo di operatore e di campo. È il criterio finale: ciò che il tuo eseguibile accetta, non ciò che dice una pagina web di un altro ramo di OpenFOAM.

fvSchemes non viene letto «tutto e subito» come un modulo di configurazione generico. All’avvio, foamRun legge controlDict, carica il modulo indicato da solver, costruisce la mesh e i campi; il dizionario fvSchemes viene aperto insieme alla mesh, ma ogni singola voce viene cercata solo nel momento in cui un operatore la richiede, cioè quando il modulo assembla per la prima volta il termine corrispondente. Due conseguenze pratiche: un errore su una voce mancante può comparire dopo diverse righe di log apparentemente normali (alla costruzione della prima equazione, non alla lettura del file); e una voce scritta ma mai richiesta non produce né errori né effetti.

Il dizionario è registrato come modificabile durante l’esecuzione: OpenFOAM controlla durante il calcolo se i dizionari di system/ sono stati modificati (un comportamento globale, regolato nel controlDict di sistema in $FOAM_ETC); se si salva una modifica a fvSchemes mentre il calcolo gira, il file viene riletto all’inizio di un passo successivo e il log segnala la rilettura. È comodo per esperimenti, ma rompe la tracciabilità: in un confronto controllato (sezione 15) non si modifica mai un caso in esecuzione. Nei casi multi-regione ogni regione ha il proprio file, system/<regione>/fvSchemes.

Un controllo rapido che il file sia sintatticamente valido, senza avviare il solutore, è chiedere a foamDictionary di stamparlo interamente: foamDictionary system/fvSchemes. Se c’è una parentesi o un punto e virgola mancante, l’errore compare subito con il numero di riga; se il file viene stampato, la sintassi è corretta (la correttezza dei nomi degli schemi, invece, si scopre solo all’avvio del modulo).

Per il canale laminare il modulo incompressibleFluid risolve l’equazione della quantità di moto e un’equazione per la pressione (cinematica, p/ρ, come visto nella Parte 3). In forma continua la quantità di moto è:

∂U/∂t + ∇·(U U) − ∇·[ν_eff (∇U + (∇U)^T − (2/3)(∇·U) I)] = −∇p

Il termine viscoso viene diviso in due parti: un laplaciano implicito ∇·(ν_eff∇U), che contiene la parte dominante e va nella matrice, e una divergenza esplicita del resto, ∇·[ν_eff dev2((∇U)^T)], dove dev2(A) = A − (2/3) tr(A) I. Per un fluido incomprimibile la seconda parte è quasi nulla, ma il modulo la assembla comunque. Ecco come i termini diventano nomi di operatori e quale voce di fvSchemes li serve, con il dizionario didattico della sezione 3:

// incompressibleFluid, laminar: momentum-equation operators and their names
ddt(U)                              -> ddtSchemes/default
div(phi,U)                          -> divSchemes/div(phi,U)
laplacian(nuEff,U)                  -> laplacianSchemes/default
div((nuEff*dev2(T(grad(U)))))       -> divSchemes/div((nuEff*dev2(T(grad(U)))))
grad(p)                             -> gradSchemes/default
grad(U)  (needed by linearUpwind)   -> gradSchemes/default or gradSchemes/grad(U)

Nel regime laminare ν_eff coincide con la viscosità ν di constant/physicalProperties (Parte 7), ma il nome dell’operatore resta nuEff perché è il modello di trasporto della quantità di moto a costruirlo. L’equazione della pressione aggiunge un laplaciano di p (servito da laplacianSchemes/default), interpolazioni sulle facce (servite da interpolationSchemes/default) e gradienti normali (serviti da snGradSchemes/default). La divergenza dei flussi sulle facce nell’equazione di continuità, invece, non richiede alcuna voce: la divergenza di un campo già definito sulle facce è una semplice somma, senza interpolazione da scegliere. È il motivo per cui un caso laminare incompressibleFluid funziona con solo due voci esplicite in divSchemes.

Prendiamo tre celle in fila, 1 – 2 – 3, con le facce interne A (tra 1 e 2) e B (tra 2 e 3), una faccia di ingresso a sinistra della cella 1 e una di uscita a destra della cella 3. In OpenFOAM ogni faccia interna appartiene a una cella owner (quella con indice minore) e a una cella neighbour; il vettore area S_f punta dall’owner al neighbour, e il flusso phi della faccia è un unico numero memorizzato una sola volta. Quando si somma attorno alla cella owner il flusso entra con segno +, attorno alla neighbour con segno −. Supponiamo phi_in = 2,5 × 10⁻⁶ m³/s (la portata del canale a 0,05 m/s), phi_A = 2,5 × 10⁻⁶, phi_B = 2,5 × 10⁻⁶, phi_out = 2,5 × 10⁻⁶ m³/s.

Il bilancio netto in uscita da ogni cella vale: cella 1, phi_A − phi_in = 0; cella 2, phi_B − phi_A = 0; cella 3, phi_out − phi_B = 0. Sommando le tre celle, i flussi interni A e B compaiono una volta con il segno + e una volta con il segno − e si cancellano («somma telescopica»): resta phi_out − phi_in. È la conservazione discreta: ciò che il dominio guadagna o perde dipende soltanto dai flussi sul contorno, qualunque sia lo schema scelto per calcolare i valori sulle facce, purché ogni faccia abbia un solo valore di flusso. Con un metodo alle differenze finite su griglia non strutturata questa proprietà non è automatica.

Conservazione non significa accuratezza: anche uno schema molto diffusivo come upwind è perfettamente conservativo. Le scelte in fvSchemes cambiano quale valore si attribuisce alla faccia, non il fatto che lo stesso valore venga usato da entrambe le celle.

Nella figura 2 il termine sorgente S non ha un sotto-dizionario. Una sorgente è un termine di volume, non una derivata: nella forma discreta diventa semplicemente S_PV_P, senza bisogno di interpolazione. In OpenFOAM 14 le sorgenti aggiuntive (forze di volume imposte, sorgenti di massa o di calore, zone porose) sono definite in constant/fvModels, e i vincoli (limiti di velocità o temperatura, per esempio) in system/fvConstraints. Una sorgente che dipende dalla variabile stessa viene di solito linearizzata come S = S_u + S_p·a: la parte S_p (se negativa) va sulla diagonale della matrice e migliora la stabilità, la parte S_u va nel termine noto.

Anche il gradiente di pressione nella quantità di moto, pur essendo una derivata, non si trova «per nome» nell’equazione generica: è un gradiente come gli altri e usa gradSchemes (in pratica, nell’algoritmo pressione–velocità, la parte che conta è il gradiente sulle facce, coerente con il laplaciano della pressione). Ricordalo quando leggi i log: il fatto che un termine non compaia nella figura non significa che non abbia uno schema.

Per un confronto con un file reale, ecco il contenuto (senza l’intestazione) di tutorials/incompressibleFluid/airFoil2D/system/fvSchemes di OpenFOAM-14, uno dei due esempi citati nell’appendice. È un caso stazionario con il modello di turbolenza Spalart–Allmaras: si notano steadyState, il prefisso bounded, una seconda voce di avvezione per la variabile di turbolenza nuTilda e il sotto-dizionario wallDist, richiesto dal modello.

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;

    div(phi,U)      bounded Gauss linearUpwind grad(U);
    div(phi,nuTilda) bounded Gauss linearUpwind grad(nuTilda);
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method          meshWave;
}

Le differenze rispetto al file didattico sono tutte spiegabili con le equazioni attive: il regime stazionario cambia ddtSchemes e motiva bounded (sezione 10), il modello di turbolenza aggiunge un’equazione di trasporto (quindi una divergenza in più) e la distanza dalle pareti (quindi wallDist, sezione 13). Il resto è identico.

I dizionari OpenFOAM possono contenere macro ($nome), direttive #include e chiavi espresse come espressioni regolari. Il file che si vede nell’editor può quindi differire da quello che il programma interpreta. foamDictionary mostra la versione interpretata:

foamDictionary -expand system/fvSchemes                    # whole file, macros and #include resolved
foamDictionary -keywords -entry divSchemes system/fvSchemes  # only the keywords of one sub-dictionary
foamDictionary -entry "divSchemes/div(phi,U)" -value system/fvSchemes

Le virgolette attorno a divSchemes/div(phi,U) sono necessarie: le parentesi hanno un significato speciale per la shell. Il risultato atteso per l’ultimo comando, con il file didattico, è la sola riga Gauss linearUpwind grad(U). Se invece compare un errore del tipo «keyword … is undefined», la voce non esiste con quel nome esatto: spesso per uno spazio di troppo, come in div(phi, U), che per OpenFOAM è una chiave diversa.

Quando il modulo costruisce un operatore, genera una stringa con il suo nome a partire dai nomi dei campi coinvolti: fvm::div(phi, U) produce div(phi,U), fvc::grad(p) produce grad(p), fvm::laplacian(nuEff, U) produce laplacian(nuEff,U). Poi cerca quella stringa nel sotto-dizionario della famiglia:

1. se esiste una chiave identica (o una chiave-espressione regolare che la riconosce), usa quel valore;
2. altrimenti, se esiste default con un valore diverso da none, usa il default;
3. altrimenti (default assente o none) si ferma con un errore fatale di input che nomina la chiave cercata e il dizionario in cui l’ha cercata.

Il confronto è letterale: div(phi,U) e div(phi, U) sono chiavi diverse, così come grad(U) e grad(u). Il messaggio d’errore tipico, per una divergenza mancante con default none, ha questa forma (percorsi e numeri di riga dipendono dal caso):

--> FOAM FATAL IO ERROR:
keyword div((nuEff*dev2(T(grad(U))))) is undefined in dictionary ".../system/fvSchemes/divSchemes"

file: .../system/fvSchemes/divSchemes from line 32 to line 36.

Il messaggio è la prova più affidabile di quale operatore il modulo richieda davvero: basta copiare la chiave dal log nel file e assegnarle uno schema adatto.

Con default none i file dei casi turbolenti diventerebbero lunghi e ripetitivi. I tutorial Foundation usano due strumenti del linguaggio dei dizionari (Parte 1) per evitarlo. Il primo è la macro: si definisce una voce con un nome qualsiasi e la si richiama con $. Il secondo è la chiave espressione regolare, scritta tra virgolette, che vale per tutti i nomi che riconosce. Un esempio dello schema tipico:

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);
    turbulence      bounded Gauss limitedLinear 1;
    div(phi,k)      $turbulence;
    div(phi,omega)  $turbulence;
    "div\(phi,(k|epsilon)\)"  bounded Gauss upwind;   // regular-expression alternative
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

La voce turbulence non è un operatore: nessun modulo la richiede, quindi viene ignorata come tale e serve solo come «variabile». In caso di conflitto, una chiave letterale ha la precedenza su un’espressione regolare (nell’esempio, div(phi,k) userebbe $turbulence, non upwind); in un file reale non si mescolano le due tecniche per la stessa chiave, proprio per evitare ambiguità. foamDictionary -expand mostra il risultato finale.

Indicando con n il livello temporale noto e con n+1 quello da calcolare, la derivata di cella (∂a/∂t)_P viene approssimata così:

Euler          (a^(n+1) − a^n) / Δt                                   error O(Δt)
backward       (3 a^(n+1) − 4 a^n + a^(n−1)) / (2 Δt)                 error O(Δt²), uses two old levels
CrankNicolson  trapezoidal average of the RHS between n and n+1,
               blended with Euler by the coefficient ψ                error O(Δt²) for ψ = 1
steadyState    0

«Implicito» significa che il resto dell’equazione (avvezione, diffusione) è valutato al nuovo livello n+1: il sistema va risolto, ma non esiste un limite di stabilità rigido sul passo come per gli schemi espliciti. Il passo resta però limitato dall’accuratezza e, nell’algoritmo PIMPLE, dalla convergenza dell’accoppiamento pressione–velocità. backward ha bisogno di due livelli passati: al primo passo, quando aⁿ⁻¹ non esiste, OpenFOAM usa automaticamente una formula del primo ordine. Per CrankNicolson la User Guide v14 spiega che il coefficiente ψ vale 1 per il Crank–Nicolson puro e 0 per Euler, e che in genere si usa ψ = 0,9 per stabilizzare lo schema: la sintassi è default CrankNicolson 0.9;. Se si omette il numero, la lettura fallisce: è il senso di «il parametro è obbligatorio» nella tabella.

Esempio numerico. Prendiamo il decadimento da/dt = −a/τ con τ = 0,01 s, a(0) = 1 e il passo della Parte 4, Δt = 0,002 s (Δt/τ = 0,2). Dopo un passo la soluzione esatta vale e^−0,2 = 0,8187. Euler implicito dà 1/(1 + 0,2) = 0,8333 (errore 0,0146, 1,8 %); Crank–Nicolson puro dà (1 − 0,1)/(1 + 0,1) = 0,8182 (errore 0,0005, 0,06 %). Dimezzando Δt l’errore di Euler circa si dimezza, quello di Crank–Nicolson si riduce di circa quattro volte: è il significato concreto di «primo» e «secondo ordine».

Nel canale il numero di Courant è piccolo: con U = 0,05 m/s, Δt = 0,002 s e Δx = 2,5 mm, Co = UΔt/Δx = 0,04. Per un caso così l’errore temporale di Euler è modesto; lo diventa meno nei transitori rapidi, nelle scie e nei flussi periodici, dove si preferiscono backward o CrankNicolson 0.9.

Nella forma discreta, il termine temporale implicito aggiunge V_P/Δt alla diagonale di ogni cella: è un contributo che rende la matrice più «diagonalmente dominante» e quindi più facile da risolvere e più stabile da iterare. Con steadyState questo contributo sparisce. Per non perdere stabilità, il solutore stazionario sostituisce il ruolo del passo temporale con il sotto-rilassamento, impostato in fvSolution (relaxationFactors, Parte 6), e usa un algoritmo di accoppiamento pressione–velocità diverso (SIMPLE invece di PIMPLE con più correttori).

Inoltre il «tempo» cambia significato: in un calcolo stazionario deltaT 1; e endTime 2000; indicano 2000 iterazioni, non 2000 secondi, e le directory dei risultati sono numerate per iterazione. Gli stati intermedi non sono istanti fisici: un calcolo stazionario che oscilla non descrive un flusso oscillante, segnala soltanto che la soluzione stazionaria non è (ancora) raggiunta o non esiste. Viceversa, un transitorio lungo con Euler converge alla stessa soluzione stazionaria solo se questa esiste ed è stabile, e i suoi residui misurano la variazione tra passi, non la convergenza di un’iterazione stazionaria.

Con localEuler ogni cella avanza con il proprio passo temporale, calcolato per rispettare un numero di Courant massimo locale: le celle grandi, dove il flusso è lento, fanno passi lunghi; quelle piccole fanno passi corti. L’evoluzione non è più sincronizzata tra le celle, quindi gli stati intermedi non hanno alcun significato fisico, ma la soluzione converge verso lo stato stazionario spesso più rapidamente e in modo più robusto di un transitorio con passo globale. La User Guide v14 lo descrive come schema pseudo-transitorio implicito del primo ordine, per accelerare la convergenza stazionaria.

Non tutti i moduli supportano LTS e i parametri che lo controllano (per esempio un Courant massimo e coefficienti di smorzamento del campo del passo locale) stanno nel dizionario dell’algoritmo in fvSolution, non in fvSchemes. Il modo sicuro per scoprire sintassi e moduli compatibili nella tua installazione è cercare i tutorial che lo usano: foamSearch "$FOAM_TUTORIALS" fvSchemes ddtSchemes/default mostra se localEuler compare; poi grep -rl localEuler "$FOAM_TUTORIALS" indica in quali casi, e il loro fvSolution mostra i parametri associati.

Il teorema di Gauss applicato al gradiente dice che l’integrale di ∇a sul volume di una cella è uguale all’integrale di a·n sulla sua superficie. Approssimando il valore medio del gradiente con il suo valore al centro e l’integrale di superficie con una somma sulle facce:

(∇a)_P ≈ (1 / V_P) Σ_f S_f a_f          S_f = face area vector (outward, magnitude = face area)
a_f = w a_P + (1 − w) a_N               linear interpolation, w = weight of the owner cell

In OpenFOAM il peso lineare è calcolato con le distanze dei due centri dal piano della faccia, misurate lungo la normale: w = d_fN/(d_Pf + d_fN). Esempio: centro P a x = 0, centro N a x = 3 mm, faccia a x = 1 mm (cella P più piccola). Allora d_Pf = 1 mm, d_fN = 2 mm, w = 2/3: la faccia, più vicina a P, prende 2/3 del valore di P. Con a_P = 1 e a_N = 4 risulta a_f = (2/3)·1 + (1/3)·4 = 2, che è esattamente il valore di una funzione lineare a(x) = 1 + x (in mm) nel punto x = 1 mm. Il nome dello schema viene da qui: è esatto per variazioni lineari, quindi accurato al secondo ordine.

Sulle facce di contorno non c’è una cella N: il valore di faccia è quello della condizione al contorno della patch (Parte 3). È il punto in cui le condizioni di 0/U e 0/p entrano negli operatori discreti: un fixedValue fornisce a_f direttamente, un zeroGradient pone a_f = a_P.

I passi 3 e 4 del modello mentale hanno una struttura molto concreta in OpenFOAM. La matrice del sistema (classe fvMatrix, basata su lduMatrix) memorizza una diagonale (un valore per cella) e due vettori di coefficienti fuori diagonale, «upper» e «lower», uno per faccia interna. Ogni faccia interna collega la cella owner alla neighbour e porta il suo coefficiente in entrambe le righe. Per esempio, per la diffusione tra P e N con coefficiente di faccia Γ_f|S_f|/|d|:

row P:   diag_P   −= Γ_f |S_f| / |d|      coefficient of N in row P = + Γ_f |S_f| / |d|
row N:   diag_N   −= Γ_f |S_f| / |d|      coefficient of P in row N = + Γ_f |S_f| / |d|

Il numero di coefficienti fuori diagonale è quindi pari al doppio del numero di facce interne, e la struttura di sparsità della matrice è la topologia della mesh. Per il canale della Parte 2 (40 × 10 × 5 = 2000 celle, 5350 facce interne) la matrice ha 2000 elementi diagonali e 10700 fuori diagonale, su 4 milioni di posizioni possibili: lo 0,3 %. È questa sparsità che rende possibili i solutori iterativi della Parte 6.

Gauss linear usa solo le facce della cella e i valori interpolati su di esse; è economico, conservativo nella forma ed esatto per campi lineari su mesh regolari. Su mesh molto distorte (forte asimmetria tra le celle, cioè skewness, o celle con facce molto diverse) l’interpolazione lineare non avviene nel centro della faccia e l’accuratezza degrada. leastSquares sceglie invece il gradiente che minimizza, in senso dei minimi quadrati pesati, la differenza tra a_N − a_P e (∇a)_P·d_PN su tutte le celle vicine: è esatto per campi lineari su qualsiasi mesh e quindi più robusto in presenza di skewness, ma tende a essere più sensibile su celle molto allungate e costa di più.

La User Guide v14 indica Gauss linear come scelta usuale e menziona leastSquares e Gauss cubic come alternative usate raramente. Per il canale esaedrico regolare le due scelte danno lo stesso gradiente a meno dell’arrotondamento nelle celle interne: un buon test di sanità, non un esperimento interessante.

Il limitatore cellLimited calcola prima il gradiente con lo schema che lo segue (qui Gauss linear), poi controlla ogni cella: per ogni faccia estrapola il valore a_P + (∇a)_P·(x_f − x_P) e verifica che non esca dall’intervallo [min, max] dei valori della cella e delle sue vicine. Se esce, moltiplica l’intero gradiente della cella per il fattore (≤ 1) che riporta l’estrapolazione peggiore dentro l’intervallo. Il coefficiente finale, tra 0 e 1, regola quanto margine si concede: 1 impone il limite in modo stretto (la «limitazione forte» della tabella), 0 equivale a non limitare.

Esempio. In una fila di celle con passo 1 mm i valori sono 0, 1, 1 (W, P, E). Il gradiente centrale in P è (1 − 0)/2 = 0,5 per mm. Estrapolato alla faccia est (0,5 mm) dà 1 + 0,25 = 1,25: supera il massimo locale 1. Il limitatore riduce il gradiente a 0 in quella cella (l’unico fattore che non supera 1 a est), e il valore di faccia diventa 1. Senza limitatore, uno schema come linearUpwind che usa questo gradiente creerebbe un nuovo massimo, cioè un’oscillazione non fisica.

La User Guide v14 mostra cellLimited Gauss linear 1 per i gradienti di velocità e di turbolenza (grad(U), grad(k), grad(epsilon)) per mantenere la limitatezza su mesh di qualità scadente. Il codice sorgente di OpenFOAM contiene anche varianti come cellMDLimited (limitazione direzionale, meno diffusiva) e faceLimited; non sono trattate nella guida, quindi prima di usarle verifica nell’installazione che esistano (per esempio con l’elenco «Valid schemes» descritto nella sezione iniziale) e cerca un tutorial che le impieghi.

L’argomento di linearUpwind (sezione 8) non è una formula ma il nome di una voce di gradSchemes: linearUpwind grad(U) significa «per la ricostruzione usa il gradiente definito da gradSchemes/grad(U), o dal default se manca». Ne segue che un override come grad(U) cellLimited Gauss linear 1; agisce anche sulla ricostruzione convettiva della velocità. Alcuni tutorial Foundation sfruttano questo meccanismo definendo in gradSchemes una voce con nome libero, per esempio limited cellLimited Gauss linear 1;, e richiamandola in divSchemes come linearUpwind limited: così la velocità usa un gradiente limitato solo nella ricostruzione convettiva, mentre gli altri usi di grad(U) restano non limitati. Prima di copiare questo schema, controlla con foamSearch come lo scrive il tutorial della tua versione.

phi è un campo di superficie (surfaceScalarField): un numero per faccia, F_f = S_f·U_f, con dimensioni [m³/s] nei moduli incomprimibili. Nei moduli comprimibili lo stesso nome indica il flusso di massa ρS_f·U_f [kg/s], ed è per questo che la stessa stringa div(phi,U) può avere un significato fisico leggermente diverso a seconda del modulo. Il modulo crea phi all’avvio dalla velocità iniziale (o lo legge, se esiste un file phi nella directory del tempo iniziale) e lo aggiorna dopo ogni correzione della pressione in modo che rispetti la continuità discreta; quando scrive i risultati compare un file phi in ogni directory temporale.

Numeri del canale. Una faccia di ingresso della mesh della Parte 2 misura 1 mm × 1 mm = 1,0 × 10⁻⁶ m²; con U = 0,05 m/s il suo flusso è 5,0 × 10⁻⁸ m³/s, e le 50 facce di ingresso sommano 2,5 × 10⁻⁶ m³/s, la portata del canale. Il contributo convettivo di una faccia all’equazione della quantità di moto è F_f·U_f: un flusso di quantità di moto (per unità di densità) [m⁴/s²]. Lo schema in div(phi,U) decide solo U_f; F_f viene dall’accoppiamento pressione–velocità.

Il flusso si può esaminare come qualsiasi altro campo: foamDictionary -entry internalField 0.02/phi | head mostra i primi valori (uno per faccia interna, nell’ordine delle facce), e la somma sulle patch è calcolabile con il function object patchFlowRate, che si può lanciare a posteriori, per esempio con foamPostProcess -func "patchFlowRate(patch=outlet)" -latestTime (verifica gli argomenti esatti con foamPostProcess -list e con il file in $FOAM_ETC/caseDicts/functions).

Nella voce div((nuEff*dev2(T(grad(U))))) il campo di cui si calcola la divergenza è un tensore già noto, valutato esplicitamente con i valori correnti: ν_eff per la parte deviatorica della trasposta del gradiente di velocità. Non c’è un «verso di trasporto» né un flusso che porta informazione da monte: si tratta di una grandezza diffusiva, per cui l’interpolazione centrale linear, del secondo ordine, è quella naturale. Gli schemi upwind non avrebbero significato: lo stesso upwind richiede un flusso per sapere quale cella sta a monte, e qui non ce n’è uno.

Per l’algebra: ∇·[ν(∇U + (∇U)^T)] = ∇·(ν∇U) + ∇·[ν(∇U)^T]. Con ν costante e ∇·U = 0 il secondo termine si annulla esattamente (è ν∇(∇·U)), ma nella forma discreta, con ν_eff variabile (turbolenza) o continuità non ancora soddisfatta, non è zero. La User Guide v14 indica che i termini non avvettivi in divSchemes usano tipicamente Gauss linear. Nei modelli non lineari compaiono altre voci analoghe (per esempio una divergenza di sforzi non lineari): stesso ragionamento, stesso schema.

Mesh uniforme monodimensionale, passo Δx, flusso positivo (da sinistra a destra). Tre celle consecutive W, P, N con valori a_W = 0, a_P = 1, a_N = 3; cerchiamo il valore sulla faccia f tra P e N. Poiché il flusso va da P a N, la cella a monte (upwind) è P.

Il rapporto r confronta la pendenza a monte (a_P − a_W = 1) con quella attraverso la faccia (a_N − a_P = 2). Dove il profilo è regolare (r vicino a 1 o, come qui, non troppo piccolo) limitedLinear si comporta come linear; dove la pendenza cambia bruscamente o di segno (r ≤ 0, cioè un massimo o minimo locale) torna a upwind, che non può creare nuovi estremi. Il coefficiente 1 è, secondo la User Guide v14, la limitazione più forte; valori più piccoli (verso 0) limitano meno. La formula esatta del limitatore va presa dal codice sorgente; l’esempio serve a capirne il comportamento, non a sostituirlo.

Le quattro risposte diverse (1; 1,75; 1,94; 2) mostrano perché la scelta conta: il valore di faccia moltiplica il flusso F_f e quindi cambia la quantità di moto trasportata. Il teorema di Godunov spiega perché non esiste uno schema perfetto: nessuno schema lineare di ordine superiore al primo può essere monotono (privo di nuove oscillazioni). Gli schemi limitati aggirano il teorema essendo non lineari: il loro peso dipende dalla soluzione stessa.

Uno sviluppo in serie di Taylor mostra che lo schema upwind risolve, invece dell’equazione originale, un’equazione con una diffusione aggiuntiva pari a circa Γ_num ≈ |U|Δx/2 nella direzione del flusso. Per il canale (U = 0,05 m/s, Δx = 2,5 mm) Γ_num ≈ 6,3 × 10⁻⁵ m²/s, circa 60 volte la viscosità fisica ν = 1,0 × 10⁻⁶ m²/s. Lo stesso confronto si esprime con il numero di Péclet di cella, Pe = UΔx/ν = 125: quando Pe ≫ 2, la diffusione numerica dell’upwind domina su quella fisica.

Questo non significa che i risultati del canale saranno sbagliati di un fattore 60. La diffusione numerica agisce dove il campo trasportato varia lungo la direzione del flusso; nel tratto sviluppato di un canale il profilo di velocità quasi non cambia lungo x, e il termine convettivo è piccolo. Le differenze tra upwind e linearUpwind sono quindi attese soprattutto nella regione d’ingresso, dove il profilo uniforme si trasforma nel profilo parabolico, e nei transitori iniziali. È una previsione da verificare nel laboratorio della sezione 15, non un risultato già garantito.

• Varianti V (limitedLinearV 1, linearUpwindV grad(U)): per campi vettoriali, usano un unico limitatore per tutte le componenti, calcolato nella direzione di massima variazione, invece di limitare ogni componente separatamente. La User Guide v14 ne riporta esempi per div(phi,U), anche con il prefisso bounded.
• vanLeer, vanAlbada: altri schemi limitati di tipo TVD, usati soprattutto per scalari.
• Varianti 01 (per esempio limitedLinear01 1): limitano più strettamente tra 0 e 1, per scalari che per definizione stanno in quell’intervallo (frazioni, variabili di progresso).
• multivariateSelection: applica un limitatore comune a più equazioni accoppiate, per esempio frazioni di massa delle specie ed entalpia, perché restino coerenti tra loro.

Sono tutti strumenti per problemi specifici; per il canale laminare con la sola velocità i cinque schemi della tabella bastano ampiamente. La regola dell’appendice resta valida: si parte da un tutorial nativo v14 con fisica simile e si annota ogni deviazione.

Nella forma discreta, il termine div(U) di una cella è la somma dei flussi uscenti divisa per il volume: (Σ_f F_f)/V_P. L’operatore bounded Gauss … costruisce prima la convezione con lo schema indicato e poi sottrae implicitamente a_P·Σ_f F_f, cioè il termine a·div(U) dell’identità: in pratica aggiunge −Σ_fF_f alla diagonale della cella. Se la continuità è esatta, Σ_f F_f = 0 e non cambia nulla; se non lo è, il termine compensa lo squilibrio di massa della cella.

Perché serve. Immagina una cella in cui, a metà convergenza, entra più flusso di quanto ne esca: Σ_f F_f = −1 × 10⁻⁹ m³/s. Nella forma conservativa pura la cella «accumula» la grandezza trasportata in proporzione allo squilibrio, come se ci fosse una sorgente fittizia; nelle iterazioni stazionarie, dove non c’è il termine V/Δt a smorzare, questa sorgente spuria può spingere un campo fuori dai suoi limiti fisici (per esempio una grandezza di turbolenza negativa). Sottraendo a·div(U) l’operatore diventa la forma non conservativa U·∇a, che non può creare nuovi estremi a causa dell’errore di continuità. La User Guide v14 afferma che questo aiuta a mantenere la limitatezza della variabile risolta e favorisce la convergenza, soprattutto nei calcoli stazionari.

Nei calcoli transitori, con PIMPLE convergito a ogni passo, la continuità è soddisfatta fino alla tolleranza del solutore di pressione, il termine è trascurabile e la forma conservativa è preferibile: per questo la riga transitoria della sezione non ha il prefisso. Il prefisso si combina con qualsiasi schema convettivo (bounded Gauss upwind, bounded Gauss limitedLinear 1, bounded Gauss linearUpwindV grad(U)…) e non sostituisce il limitatore: sono due correzioni indipendenti.

Integrando ∇·(Γ∇a) sulla cella e applicando Gauss si ottiene una somma di flussi diffusivi di faccia: Σ_f Γ_f |S_f| (∂a/∂n)_f. Le tre parole della voce corrispondono esattamente ai tre fattori: Gauss è la somma, linear fornisce Γ_f, corrected fornisce (∂a/∂n)_f. Su una mesh ortogonale il gradiente normale è semplicemente (a_N − a_P)/|d|, con |d| la distanza tra i centri, e la parte implicita va interamente nella matrice.

Esempio. Faccia tra due celle del canale sovrapposte in y (mesh della Parte 2: |S_f| = 2,5 mm × 1 mm = 2,5 × 10⁻⁶ m², |d| = 1 mm). Con ν = 1,0 × 10⁻⁶ m²/s il coefficiente della faccia vale ν|S_f|/|d| = 1,0 × 10⁻⁶ · 2,5 × 10⁻⁶/1,0 × 10⁻³ = 2,5 × 10⁻⁹ m³/s. Se le velocità assiali delle due celle sono 0,060 e 0,066 m/s, il flusso diffusivo di quantità di moto attraverso la faccia è 2,5 × 10⁻⁹ · 0,006 = 1,5 × 10⁻¹¹ m⁴/s²: è questo scambio, sommato su tutte le facce, che produce la caduta di pressione del canale.

Secondo la User Guide v14, Gauss è l’unica scelta disponibile per i laplaciani, nella forma Gauss <schema di interpolazione> <schema snGrad>: per questo la voce ha sempre tre parti (quattro se lo schema snGrad ha un coefficiente, come in limited corrected 0.5). Gli operatori laplaciani del canale sono laplacian(nuEff,U) nella quantità di moto e un laplaciano della pressione il cui nome contiene il coefficiente costruito dall’algoritmo pressione–velocità; con un default non serve conoscerlo, ma se si usa default none il messaggio d’errore lo riporta esattamente.

Nel canale ν è uniforme e l’interpolazione del coefficiente è irrilevante: qualsiasi media di due valori uguali dà lo stesso numero. Diventa importante quando Γ varia molto da una cella all’altra: viscosità turbolenta vicino alle pareti, conducibilità termica all’interfaccia tra due materiali in un caso multi-regione, permeabilità in un mezzo poroso. In un problema di conduzione in serie tra due strati, il flusso corretto attraverso l’interfaccia è governato dalla media armonica delle conducibilità, non da quella aritmetica: con Γ_P = 1 e Γ_N = 100 la media lineare dà ≈ 50,5, quella armonica 2/(1/1 + 1/100) ≈ 1,98, una differenza di 25 volte nel flusso di faccia. OpenFOAM dispone di uno schema di interpolazione harmonic per questi casi; prima di usarlo, verifica in un tutorial v14 pertinente che sia impiegato con quella sintassi.

Sia d il vettore da P a N e S_f il vettore area della faccia. Se d è parallelo a S_f (mesh ortogonale), il gradiente normale è esattamente (a_N − a_P)/|d| al secondo ordine. Se non lo è, OpenFOAM scompone S_f in due parti: una parallela a d, Δ, che si tratta implicitamente con la differenza tra i due centri, e una correzione k = S_f − Δ, che si moltiplica per il gradiente interpolato sulla faccia (∇a)_f e si tratta esplicitamente:

S_f · (∇a)_f  ≈  |Δ| (a_N − a_P) / |d|   +   k · (∇a)_f
                  implicit (matrix)           explicit (non-orthogonal correction)

over-relaxed decomposition:  Δ = d (S_f · S_f) / (d · S_f)

La parte esplicita usa il gradiente della soluzione corrente, quindi è esatta solo se la si ricalcola dopo aver risolto: è il motivo dei correttori non ortogonali (nNonOrthogonalCorrectors in fvSolution, Parte 6), che ripetono la soluzione dell’equazione della pressione aggiornando la correzione. Su una mesh ortogonale k = 0 e la correzione non fa nulla: corrected e orthogonal danno lo stesso risultato, e i correttori aggiuntivi sono solo tempo di calcolo sprecato.

Il coefficiente ψ di limited (tra 0 e 1) limita la correzione esplicita in rapporto alla parte ortogonale: la correzione non può superare ψ/(1 − ψ) volte la parte implicita. Con ψ = 0,5 la correzione può al massimo eguagliare la parte ortogonale; con ψ = 0,33 può valerne al massimo circa la metà; ψ = 0 equivale a uncorrected e ψ = 1 a corrected. Limitare significa scambiare un po’ di accuratezza su celle molto distorte con la stabilità della correzione esplicita, che su facce quasi tangenti al vettore d può diventare dominante.

checkMesh riporta l’angolo di non ortogonalità, cioè l’angolo tra d e S_f, nella forma Mesh non-orthogonality Max: … average: …, e segnala come errore le facce oltre la sua soglia predefinita di 70°. Per il canale blockMesh a blocco singolo e rettilineo, l’atteso è Max: 0 average: 0: tutte le facce sono perfettamente ortogonali. Per una mesh non strutturata generata con Gmsh (il caso del canale con Gmsh e gmshToFoam) i valori sono tipicamente diversi da zero e vanno letti prima di scegliere snGradSchemes e il numero di correttori non ortogonali.

checkMesh > log.checkMesh 2>&1
grep -i "non-orthogonality" log.checkMesh
grep -i "skewness" log.checkMesh

• Max < 5° circa: corrected e orthogonal/uncorrected sono praticamente equivalenti; corrected resta la scelta più sicura.
• Max fino a circa 70°: corrected, eventualmente con uno o due correttori non ortogonali.
• Max oltre circa 75°: considerare limited corrected 0.5 o 0.33 (e la stessa politica in laplacianSchemes), ma prima chiedersi se la mesh vada migliorata.
• Max oltre circa 85°: la convergenza diventa difficile; il rimedio sta nella mesh, non negli schemi.

Le soglie sono quelle indicate dalla User Guide v14 e da checkMesh; sono indicazioni orientative, non regole rigide. La politica di correzione in snGradSchemes e quella nell’ultima parola di laplacianSchemes dovrebbero essere coerenti tra loro.

Le interpolazioni «già specificate da una voce più specializzata» sono quelle interne agli altri schemi: il linear di Gauss linear in gradSchemes, lo schema convettivo in divSchemes, l’interpolazione del coefficiente in laplacianSchemes. Restano le interpolazioni «nude» che il modulo chiede esplicitamente, la più importante delle quali nel modulo incompressibleFluid è la costruzione del flusso sulle facce a partire dalla velocità predetta, nel passo di correzione della pressione. Anche il calcolo iniziale di phi dalla velocità iniziale usa un’interpolazione cella-faccia.

Cambiare questo default in un caso incompressibleFluid modifica quindi l’accoppiamento pressione–velocità, non la convezione: è un esperimento sconsigliato senza un motivo preciso. La User Guide v14 menziona come eccezione l’uso di cubic in un tutorial di analisi delle sollecitazioni.

La distanza dalla parete y è un campo di volume che alcuni modelli di turbolenza usano direttamente: Spalart–Allmaras (il modello di airFoil2D), k-ω SST (nelle funzioni di raccordo), i modelli LES con smorzamento di parete. Nel caso laminare del canale nessun modello la richiede, quindi il sotto-dizionario può mancare senza conseguenze. Se invece un modello la richiede e wallDist manca, il calcolo si ferma all’avvio con un errore che segnala la voce mancante in system/fvSchemes: la correzione è aggiungerla, non cambiare modello.

meshWave calcola la distanza propagando l’informazione dalle facce di parete attraverso la mesh, cella dopo cella: è veloce e robusto, ed è il metodo dei tutorial. Esistono altri metodi selezionabili a run-time, per esempio basati sulla soluzione di un’equazione di Poisson, più regolari su mesh molto distorte; l’elenco autentico per la tua installazione si ottiene con l’errore «unknown method» o con foamToC -search meshWave.

La derivata seconda nel tempo, ∂²a/∂t², compare nelle equazioni della dinamica dei solidi (per esempio lo spostamento nel modulo di analisi delle sollecitazioni). Secondo la User Guide v14 per d2dt2Schemes è disponibile solo Euler. Per i fluidi non serve mai.

Ogni equazione di trasporto aggiunta porta con sé almeno una divergenza con il proprio nome: un modello k-ε aggiunge div(phi,k) e div(phi,epsilon), un modello k-ω div(phi,omega), l’equazione dell’energia una divergenza della variabile energetica scelta dal modulo, ogni specie una divergenza della sua frazione di massa. Nei moduli comprimibili phi diventa un flusso di massa e possono comparire termini come la divergenza dell’energia cinetica o del lavoro di pressione. Il modo affidabile per sapere i nomi esatti è leggere un tutorial nativo dello stesso modulo, oppure lasciare che default none li faccia comparire uno per uno negli errori.

foamVersion deve stampare la versione 14 (per esempio OpenFOAM-14): se stampa altro, o se il comando non esiste, l’ambiente non è attivo e i comandi successivi useranno un’altra installazione o falliranno. echo "$FOAM_TUTORIALS" stampa il percorso della directory dei tutorial dell’installazione attiva; se la riga è vuota, la variabile non è definita e l’ambiente non è stato caricato. foamSearch attraversa tutti i file fvSchemes sotto quel percorso, estrae il valore della voce richiesta e stampa l’elenco ordinato dei valori distinti, uno per riga, nella forma default Euler;, default steadyState; e così via. La lunghezza e il contenuto esatti dell’elenco dipendono dall’installazione: non c’è un output «giusto» da confrontare, c’è un inventario da leggere.

Per trasformare l’inventario in informazione, conviene risalire ai casi che usano un certo valore e confrontare il loro regime. Per esempio, per sapere quali tutorial usano bounded nella convezione della velocità e se sono stazionari:

foamSearch -count "$FOAM_TUTORIALS" fvSchemes "divSchemes/div(phi,U)"
grep -l "bounded Gauss linearUpwind" "$FOAM_TUTORIALS"/incompressibleFluid/*/system/fvSchemes
for f in "$FOAM_TUTORIALS"/incompressibleFluid/*/system/fvSchemes; do
  printf "%s : %s\n" "$f" "$(foamDictionary -entry ddtSchemes/default -value "$f")"
done

Ci si aspetta di vedere che bounded compare soprattutto nei casi con steadyState: è la regola della sezione 10 osservata nei dati, non imparata a memoria.

Nella versione pubblicata nel repository OpenFOAM-14 (link nell’appendice) il file del cilindro contiene, al momento della stesura: ddtSchemes/default steadyState, con un commento che indica Euler per l’uso transitorio; gradSchemes con default Gauss linear e una versione limitata con cellLimited per i gradienti di velocità e delle grandezze di turbolenza (U, k, omega, epsilon); divSchemes con default none, convezione bounded della velocità basata su linearUpwind con un gradiente limitato, convezione bounded con limitedLinear per le grandezze di turbolenza (riutilizzata con una macro) e altre voci pensate per moduli comprimibili, come divergenze dell’energia e della pressione; laplacianSchemes Gauss linear corrected, interpolationSchemes linear, snGradSchemes corrected e wallDist con meshWave. Apri il file installato per controllare la sintassi esatta.

Il fatto che un file contenga voci per moduli che non stai usando non è un errore: le voci mai richieste vengono ignorate (sezione 1). È una scelta degli sviluppatori per rendere il file riusabile.

• div(phi,U) contro divergenza degli sforzi. La prima è il trasporto convettivo della quantità di moto da parte del flusso di faccia phi: richiede uno schema che sappia da che parte arriva il flusso (upwind, linearUpwind…). La seconda è la parte esplicita della divergenza dello sforzo viscoso, ∇·[ν_eff dev2((∇U)^T)]: non c’è trasporto, il tensore è noto, e lo schema naturale è Gauss linear.
• Perché linearUpwind ha un gradiente. Ricostruisce il valore di faccia estrapolando dal centro della cella a monte: a_f = a_P + (∇a)_P·(x_f − x_P). Il gradiente serve a questa estrapolazione, e l’argomento dice quale voce di gradSchemes usare per calcolarlo.
• Transitorio o stazionario. La natura del problema si legge in fvSchemes: ddtSchemes/default (o una voce ddt(U) specifica) decide se la derivata temporale esiste (Euler, backward…) o è annullata (steadyState). controlDict fissa deltaT e endTime, ma con steadyState questi contano iterazioni, non secondi; fvSolution sceglie l’algoritmo coerente (PIMPLE o SIMPLE). La risposta completa è quindi: la voce ddtSchemes di fvSchemes, sostenuta da impostazioni coerenti negli altri due dizionari.
• Nota del passo 06. default none elimina solo il ripiego: dice «non scegliere uno schema per una divergenza che non ho nominato». L’avvezione della velocità è calcolata eccome, con lo schema assegnato esplicitamente a div(phi,U); se quella voce mancasse, il calcolo non partirebbe affatto invece di proseguire senza avvezione.

Il blocco della sezione mostra l’inizio del flusso di lavoro: due copie e la lettura della voce attuale. Ecco una versione completa e verificabile. foamDictionary -set modifica la voce nel file senza aprire un editor e senza toccare il resto; diff dimostra che le due copie differiscono solo in quella riga; il checksum fissa lo stato esatto dei file nel registro dell’esperimento.

# 0. Start from a clean, complete channel case (mesh, 0/, constant/, system/)
cd /path/to/your/run
cp -a channelCase channel_upwind
cp -a channelCase channel_linearUpwind

# 1. Change exactly one entry in each copy
foamDictionary -entry "divSchemes/div(phi,U)" -set "Gauss upwind" \
  channel_upwind/system/fvSchemes
foamDictionary -entry "divSchemes/div(phi,U)" -set "Gauss linearUpwind grad(U)" \
  channel_linearUpwind/system/fvSchemes

# 2. Prove that only that line differs
diff channel_upwind/system/fvSchemes channel_linearUpwind/system/fvSchemes
diff -r -q channel_upwind/constant channel_linearUpwind/constant
diff -r -q channel_upwind/0 channel_linearUpwind/0
sha256sum channel_*/system/fvSchemes > fvSchemes.sha256

# 3. Run each case in its own directory, keeping its log
(cd channel_upwind && foamRun > log.foamRun 2>&1)
(cd channel_linearUpwind && foamRun > log.foamRun 2>&1)
tail -n 5 channel_*/log.foamRun

Risultato atteso. L’output di diff sui due fvSchemes deve mostrare una sola riga cambiata, quella di div(phi,U) (foamDictionary -set può riformattare spaziature e commenti: per questo il confronto si fa tra le due copie modificate, che sono state trattate allo stesso modo). I due diff -r -q su constant e 0 non devono stampare nulla. Le ultime righe dei due log devono contenere End, segno che il calcolo è arrivato a endTime.

Attenzione alla ripartenza. Se channelCase contiene già directory dei tempi e controlDict usa startFrom latestTime (Parte 4), le copie ripartono dall’ultimo risultato del caso originale, calcolato con lo schema originale: il confronto sarebbe falsato. Prima di lanciare, controlla con ls che in ciascuna copia ci sia solo 0/, oppure imposta esplicitamente startFrom startTime in entrambe (la stessa modifica in tutte e due, annotata nel registro).

I punti 03 e 04 del registro richiedono misure fatte con lo stesso strumento, nelle stesse posizioni e agli stessi tempi. In OpenFOAM 14 lo strumento naturale sono i function object eseguiti a posteriori con foamPostProcess: graphUniform campiona i campi lungo un segmento con punti equispaziati e patchAverage calcola la media di un campo su una patch. Con il riferimento p = 0 all’uscita (Parte 3), la media di p sull’ingresso è direttamente la caduta di pressione cinematica.

# Same line, same time, both cases (adapt the coordinates to your mesh, in metres)
for c in channel_upwind channel_linearUpwind; do
  (cd "$c" && foamPostProcess -time 0.02 \
     -func "graphUniform(start=(0.08 0 0.0025), end=(0.08 0.01 0.0025), nPoints=50, fields=(U))")
  (cd "$c" && foamPostProcess -time 0.02 -func "patchAverage(name=inlet, fields=(p))")
done

Le coordinate dell’esempio presuppongono il canale 100 × 10 × 5 mm con origine in un angolo, espresso in metri: un segmento verticale a x = 80 mm, a metà della larghezza. Controlla l’estensione reale con checkMesh (riga Overall domain bounding box) e gli argomenti dei function object con foamPostProcess -list e con i file in $FOAM_ETC/caseDicts/functions. I risultati vengono scritti in postProcessing/ dentro ciascun caso, con lo stesso nome di file: confrontali con un grafico o con paste.

Con le impostazioni della Parte 4 (Δt = 0,002 s, endTime 0,02 s) il calcolo copre soltanto i primi 20 ms dopo un avvio impulsivo. In quel tempo lo strato limite cresce di circa √(νt) = √(1,0 × 10⁻⁶ · 0,02) ≈ 0,14 mm: il profilo è ancora quasi uniforme, con strati sottili vicino alle pareti. Il tempo di sviluppo viscoso sulla semi-larghezza più piccola (2,5 mm) è dell’ordine di (2,5 × 10⁻³)²/ν ≈ 6 s, e il tempo di attraversamento del canale L/U = 2 s. Un confronto a t = 0,02 s misura quindi soprattutto la risposta all’avvio; per confrontare il flusso sviluppato servono tempi di qualche secondo (migliaia di passi, comunque economici su 2000 celle) oppure la versione stazionaria dell’esercizio C.

Qualitativamente, ci si aspetta che le differenze tra upwind e linearUpwind si concentrino dove il campo varia lungo il flusso, cioè nella regione d’ingresso, e che siano piccole nel tratto sviluppato (sezione 9). Il riferimento di Poiseuille del primo caso pratico e lo studio di convergenza della mesh sono gli strumenti per decidere quale dei due risultati sia più vicino alla soluzione esatta: il confronto A/B da solo dice solo quanto la scelta numerica conta in questo caso.

I messaggi di OpenFOAM hanno una struttura stabile: una riga --> FOAM FATAL ERROR o --> FOAM FATAL IO ERROR (quest’ultima quando il problema è nella lettura di un file), il testo dell’errore, poi il file e le righe coinvolte e la funzione del codice che l’ha generato. Le forme tipiche, con percorsi e numeri di riga che variano da caso a caso:

# first fatal message, with context
grep -n -m 1 -A 15 "FOAM FATAL" log.foamRun

# which scheme keyword was missing or unknown
grep -n -i -E "undefined|Unknown discretisation scheme|Valid schemes" log.foamRun

# what the file really contains, after macros and #include
foamDictionary -expand system/fvSchemes

Leggi sempre il primo errore: quelli successivi sono spesso conseguenze. E ricorda che la chiave riportata è la prova di ciò che il modulo richiede; se non coincide con ciò che ti aspettavi (per esempio compare un termine di turbolenza in un caso che credevi laminare), il problema non è in fvSchemes ma nel modello selezionato in constant/.

Campo instabile. L’ordine dei controlli della tabella non è casuale: la maggior parte delle divergenze nasce da mesh, condizioni al contorno o passo temporale, non dallo schema. In pratica: checkMesh (non ortogonalità, skewness, rapporto d’aspetto); le condizioni di 0/U e 0/p (Parte 3), in particolare riflusso all’uscita; il numero di Courant nel log (Courant Number mean: … max: …), che in un caso transitorio con PIMPLE dovrebbe restare nell’ordine dell’unità o meno; il bilancio dei flussi (portata in ingresso uguale a quella in uscita); solo alla fine la limitatezza degli schemi, passando per esempio da linear a linearUpwind o aggiungendo cellLimited al gradiente. Annotare il primo passo in cui un residuo o il Courant esplode indica spesso la regione e la causa.

Risultati invariati. Le cause tipiche sono banali e frequenti: hai modificato il fvSchemes di un’altra copia (controlla pwd e il percorso nel comando); il calcolo è ripartito da un tempo già calcolato (startFrom latestTime); stai confrontando la directory 0 o un tempo diverso; la voce modificata non viene mai richiesta (per esempio hai cambiato div(phi,k) in un caso laminare); oppure la modifica è reale ma il suo effetto è piccolo in quel caso, come nel tratto sviluppato del canale. Il checksum e la data dei file (ls -l --time-style=full-iso system/fvSchemes log.foamRun) dicono se il log è successivo alla modifica.

«Chiedere a OpenFOAM di verificare» ha un significato operativo per ciascun esercizio. Per A, in una copia del caso: imposta default none in una famiglia alla volta e osserva quale chiave compare nell’errore all’avvio: è la prova di quale sotto-dizionario serve ciascun operatore (ripristina subito dopo). Per B: nel caso del canale, con quel dizionario, la seconda divergenza div((nuEff*dev2(T(grad(U))))) non è specificata, quindi la previsione corretta è un errore fatale di voce mancante; esegui e confronta. Per C: foamDictionary -entry ddtSchemes/default -value system/fvSchemes deve restituire steadyState, e il log di un calcolo stazionario mostra il «tempo» che avanza di 1 a ogni iterazione. Per D: i comandi di ricerca sono nel riquadro delle soluzioni.

Una configurazione candidata coerente con la guida v14 e con il tutorial airFoil2D cambia due righe rispetto al file della sezione 3: lo schema temporale e il prefisso bounded sulla convezione. Tutto il resto resta identico, perché le altre famiglie non dipendono dal regime temporale.

ddtSchemes
{
    default     steadyState;
}
gradSchemes
{
    default     Gauss linear;
}
divSchemes
{
    default     none;
    div(phi,U)  bounded Gauss linearUpwind grad(U);
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
    default     Gauss linear corrected;
}
interpolationSchemes
{
    default     linear;
}
snGradSchemes
{
    default     corrected;
}

• fvSolution: serve il dizionario dell’algoritmo stazionario (SIMPLE, con residualControl come criterio di arresto) e fattori di sotto-rilassamento in relaxationFactors, perché manca il termine V/Δt che stabilizzava il transitorio (Parte 6).
• controlDict: deltaT 1;, endTime pari al numero massimo di iterazioni (per esempio 2000) e writeInterval espresso in iterazioni; il modulo resta incompressibleFluid (Parte 4).
• Criterio di convergenza: residui iniziali scesi sotto la soglia scelta e grandezze di interesse (caduta di pressione, profilo a valle) che non cambiano più tra iterazioni; il solo calo dei residui non basta.
• Confronto con il transitorio: la soluzione stazionaria va confrontata con il transitorio sviluppato (qualche secondo fisico, sezione 15), non con lo stato a 0,02 s.

Il primo passo è l’inventario, il secondo è risalire ai file e al loro regime:

foamSearch "$FOAM_TUTORIALS" fvSchemes "divSchemes/div(phi,U)"
grep -l "div(phi,U)" "$FOAM_TUTORIALS"/incompressibleFluid/*/system/fvSchemes
foamDictionary -entry "divSchemes/div(phi,U)" -value \
  "$FOAM_TUTORIALS"/incompressibleFluid/airFoil2D/system/fvSchemes
foamDictionary -entry "ddtSchemes/default" -value \
  "$FOAM_TUTORIALS"/incompressibleFluid/airFoil2D/system/fvSchemes

Un esempio verificato nel repository OpenFOAM-14: in incompressibleFluid/airFoil2D/system/fvSchemes la voce è div(phi,U) bounded Gauss linearUpwind grad(U); con ddtSchemes/default steadyState. Il caso è stazionario e turbolento (Spalart–Allmaras): il prefisso bounded segue la regola della sezione 10, e lo schema linearUpwind del secondo ordine è una scelta comune per la quantità di moto. Il secondo esempio va scelto tra i tutorial che foamSearch elenca nella tua installazione, idealmente uno transitorio (ddtSchemes/default pari a Euler o backward): ci si aspetta che la sua voce div(phi,U) sia senza bounded, e che lo schema rifletta le esigenze del caso (per esempio uno schema poco diffusivo come LUST in un calcolo che risolve strutture turbolente nel tempo). La risposta corretta documenta file, voce e regime di entrambi e collega le differenze alle equazioni e al regime, senza dichiarare uno dei due «migliore».

La User Guide descrive gli schemi più usati; l’elenco completo e il comportamento esatto stanno nel codice. Nell’installazione v14 gli schemi ai volumi finiti sono sotto $FOAM_SRC/finiteVolume/interpolation/surfaceInterpolation/schemes (schemi di interpolazione e convettivi, per esempio linearUpwind, LUST, limitedLinear), $FOAM_SRC/finiteVolume/finiteVolume/gradSchemes (gradienti e limitatori), $FOAM_SRC/finiteVolume/finiteVolume/snGradSchemes, $FOAM_SRC/finiteVolume/finiteVolume/ddtSchemes, $FOAM_SRC/finiteVolume/finiteVolume/convectionSchemes (dove si trova anche la variante bounded) e $FOAM_SRC/finiteVolume/finiteVolume/laplacianSchemes. L’intestazione di ogni file .H contiene una descrizione e, spesso, la sintassi d’uso. Un modo rapido per orientarsi:

ls "$FOAM_SRC"/finiteVolume/finiteVolume/ddtSchemes
find "$FOAM_SRC"/finiteVolume -name "linearUpwind*.H"
sed -n '1,80p' "$(find "$FOAM_SRC"/finiteVolume -name linearUpwind.H | head -1)"

Se la struttura delle directory della tua installazione differisce, find "$FOAM_SRC" -name "<nome>.H" trova comunque il file. È lo stesso principio di tutta la parte: la prova sta nell’installazione v14 che userai, non in una fonte di un altro ramo.