OpenFOAM 14, parte 7: proprietà fisiche e modelli

Questa parte è costruita come una sequenza di unità brevi, ognuna delle quali risponde a una sola domanda. Le sezioni 1–5 fissano il livello del materiale (che cosa contiene physicalProperties, come si definisce e si misura la viscosità, come v14 gestisce le unità, dove trovare un esempio reale). Le sezioni 6–8 introducono il modello dello sforzo (laminare, l’idea di chiusura, il numero di Reynolds come guida). Le sezioni 9–15 percorrono i rami della turbolenza (RAS con kEpsilon e kOmegaSST, i campi che richiedono, le stime all’ingresso, le wall function, i coefficienti, LES). Le sezioni 16–18 mostrano i livelli confinanti (reologia non newtoniana, disponibilità dei modelli, pacchetti termofisici). Le sezioni 19–26 riguardano il metodo: verifica delle dipendenze, due esperimenti, diagnostica, una valutazione autonoma, soluzioni e riferimenti.

Ogni compito guidato va svolto in un terminale, con l’ambiente di OpenFOAM 14 caricato, su una copia di un caso. Le risposte ai compiti numerici si trovano nella sezione 25; i riquadri di approfondimento aggiungono la teoria, i dettagli dell’implementazione in v14 e i risultati attesi che il testo principale lascia impliciti.

Un suggerimento pratico: apri due terminali. Nel primo tieni la copia del tuo canale laminare (quello del caso di Poiseuille); nel secondo una copia del tutorial ufficiale pitzDailySteady. Molte sezioni chiedono di confrontare i due affiancati, e il contrasto tra un caso laminare e un caso RAS è il principale strumento didattico di questa parte.

In v14 esiste un solo eseguibile, foamRun, che carica il modulo solutore indicato dalla voce solver in system/controlDict (in questo corso solver incompressibleFluid;). È il modulo a decidere quali file di constant/ gli servono. incompressibleFluid costruisce un modello di viscosità a partire da physicalProperties e un modello di trasporto della quantità di moto a partire da momentumTransport; un modulo comprimibile come fluid costruisce invece da physicalProperties un pacchetto termofisico e legge voci del tutto diverse.

Lo si vede accadere all’inizio del log: dopo la lettura della mesh, il modulo stampa le righe che annunciano la selezione del modello di viscosità e del modello di trasporto della quantità di moto (per un caso laminare, il modello di sforzo laminare selezionato, per default quello newtoniano). Leggere queste poche righe è il modo più rapido per confermare che OpenFOAM ha interpretato i tuoi dizionari come intendevi.

foamDictionary system/controlDict -entry solver -value
foamRun 2>&1 | tee log.foamRun
grep -i "select" log.foamRun

La tabella sopra è quindi una mappa delle responsabilità: se un calcolo fallisce o dà un risultato inatteso, il livello che possiede la scelta è il primo posto dove guardare. Una viscosità sbagliata è un problema di constant/; un campo iniziale mancante è un problema di 0/; una voce div(phi,k) mancante è un problema di system/, anche se la causa è una decisione presa in constant/momentumTransport.

Invece di aprire il file in un editor, interroga esattamente le voci che il solutore leggerà. foamDictionary analizza il file con lo stesso lettore usato da OpenFOAM, quindi un errore di sintassi emerge qui, prima di un calcolo:

foamDictionary constant/physicalProperties -entry viscosityModel -value
foamDictionary constant/physicalProperties -entry nu -value
foamDictionary constant/momentumTransport -entry simulationType -value

Risposte attese per il canale laminare: constant, 1e-05 (o il valore usato dal tuo caso di base) e laminar. La voce FoamFile/object deve essere physicalProperties: è il nome con cui il file viene registrato e deve coincidere con il nome del file.

Risposta al compito 3. Se cambia solo nu, non cambiano la mesh (constant/polyMesh), l’elenco dei campi in 0/ e le loro condizioni al contorno, momentumTransport, fvSchemes, fvSolution e controlDict. Cambia la soluzione: il termine viscoso dell’equazione della quantità di moto viene moltiplicato per un coefficiente diverso, quindi lo sviluppo del profilo di velocità, la caduta di pressione e il numero di Reynolds cambiano. Se il caso è definito da una velocità imposta all’ingresso, la portata resta la stessa e risponde la caduta di pressione; se fosse guidato da una differenza di pressione imposta, risponderebbe invece la portata.

Errore tipico. Se nu è scritto male o manca, il calcolo si ferma durante la costruzione del modello di viscosità con un errore fatale di IO che dichiara la parola chiave nu non definita nel dizionario …/constant/physicalProperties. Il messaggio indica sia la parola chiave sia il file: leggilo prima di modificare qualsiasi altra cosa.

Le note di rilascio di v14 indicano che le proprietà fisiche dovrebbero specificare le unità insieme al valore invece di ripetere un insieme di dimensioni. Un numero semplice viene letto nelle unità di base dell’insieme di unità attivo (SI per default), quindi nu 1e-05; significa 1 × 10⁻⁵ m²/s. Scrivere esplicitamente l’unità aggiunge un controllo e permette la conversione in lettura:

viscosityModel    constant;
nu                1e-05 [m^2/s];

// equivalente, il valore viene convertito in m^2/s alla lettura:
// nu             10 [cSt];

Le due forme danno lo stesso valore interno: 10 cSt = 10 × 10⁻⁶ m²/s = 1 × 10⁻⁵ m²/s. Se scrivi un’unità del tipo sbagliato (per esempio [m/s]), il controllo dimensionale della proprietà dimensionata rileva l’incoerenza. Se scrivi un’unità corretta con un ordine di grandezza sbagliato, nulla lo rileva: vedi la sezione 4.

Nelle versioni Foundation più vecchie e in molti tutorial reperibili in rete la proprietà era scritta con un insieme di dimensioni a indici, nu [0 2 -1 0 0 0 0] 1e-05;. Va letto come “dimensioni di m²/s, valore 1e-05”; in v14 lo stile preferito è il valore nudo oppure il valore seguito da un’unità.

Per un fluido newtoniano lo sforzo viscoso è τ = 2μD, dove D = ½(∇U + ∇U^T) è il tensore della velocità di deformazione (in un flusso incomprimibile la parte volumetrica si annulla). Con densità costante l’equazione della quantità di moto può essere divisa per ρ:

dU/dt + div(U U) = -grad(p/rho) + div( 2 (mu/rho) D )
                 = -grad(p_k)    + div( 2 nu D )

La densità scompare come parametro indipendente: compare solo ν = μ/ρ, e la pressione risolta è la pressione cinematica p_k = p/ρ, in m²/s². Per questo incompressibleFluid legge nu e non rho, e per questo 0/p in un caso simile è in m²/s², non in Pa. Per riportare una differenza di pressione in pascal, moltiplica per la densità che assumi per il fluido: Δp [Pa] = ρ × Δp_k [m²/s²].

Risposta al compito 3. Cambiare la densità del fluido reale, a parità di viscosità dinamica, cambia ν. Nel canale incomprimibile quel cambiamento entra solo attraverso nu; non esiste una voce separata per la densità da modificare, e il campo di pressione cinematica resta indipendente da ρ finché non lo converti per presentare i risultati.

Valori di riferimento utili (approssimati, 20 °C, 1 atm): acqua μ ≈ 1,0 × 10⁻³ Pa·s, ρ ≈ 998 kg/m³, ν ≈ 1,0 × 10⁻⁶ m²/s; aria μ ≈ 1,8 × 10⁻⁵ Pa·s, ρ ≈ 1,2 kg/m³, ν ≈ 1,5 × 10⁻⁵ m²/s. Nota che in termini cinematici l’aria è circa 15 volte più “viscosa” dell’acqua, benché la sua viscosità dinamica sia circa 55 volte più bassa.

Unità d’uso pratico. 1 P (poise) = 0,1 Pa·s, quindi 1 cP = 1 mPa·s = 10⁻³ Pa·s (dinamica). 1 St (stokes) = 10⁻⁴ m²/s, quindi 1 cSt = 10⁻⁶ m²/s (cinematica). L’acqua a circa 20 °C è vicina a 1 cP e a 1 cSt, ed è l’origine di molte confusioni: i numeri coincidono solo perché ρ ≈ 1000 kg/m³.

Internamente ogni grandezza dimensionata porta con sé un insieme di dimensioni di sette esponenti nell’ordine massa, lunghezza, tempo, temperatura, quantità di sostanza (moli), corrente, intensità luminosa. La viscosità cinematica è [0 2 -1 0 0 0 0]; in v14 la stessa cosa può essere scritta con un nome, [kinematicViscosity], o con un’espressione come [area/time] o [sqr(length)/time]. I casi di esempio di v14 hanno sostituito gli insiemi di dimensioni numerici con quelli con nome.

Un’unità dice invece come è espresso il numero. Le unità si scrivono tra parentesi quadre dopo un valore ([cm/s], [cm s^-1], [mm], [bar], [cSt], [deg], [%]) e il valore viene convertito nelle unità di base dell’insieme di unità attivo al momento della lettura. Gli insiemi di unità (SI per default; sono preconfigurati anche CGS e USCS) sono definiti nella sezione units di $FOAM_ETC/configDict, e un caso può ridefinirli in system/configDict.

foamUnits (in v14 un’applicazione, User Guide §4.7.7) elenca le dimensioni e le unità che corrispondono a ciascuna. Per pressure la guida mostra un output di questa forma:

Dimension [pressure]
+ Dimensions = [mass length^-1 time^-2]
+ Exponents = [1 -1 -2 0 0 0 0]
+ Standard Units = [kg m^-1 s^-2] [Pa]
+ Other Units = [psi] [atm] [kPa] [bar] [MPa]

Eseguendo foamUnits kinematicViscosity si ottiene il resoconto analogo: esponenti [0 2 -1 0 0 0 0], unità standard m²/s ed eventuali altre unità definite per quella dimensione nella tua configurazione.

Che cosa viene controllato e che cosa no. Le proprietà e i campi dimensionati confrontano l’unità fornita con le dimensioni codificate nel solutore; i parametri scalari semplici (un fattore di rilassamento, una tolleranza) non sono soggetti a controllo dimensionale. Quando si sommano o si uguagliano due grandezze con dimensioni diverse, OpenFOAM si ferma con un errore fatale che dichiara che le dimensioni differiscono per l’operazione (il messaggio mostra i due insiemi di dimensioni). Un ordine di grandezza sbagliato — 1e-5 battuto come 1e-6, un valore in cSt scritto senza unità — supera ogni controllo: confrontalo con una stima indipendente fatta a mano, come nella sezione 3.

La sequenza completa, lavorando su una copia (mai dentro $FOAM_TUTORIALS):

mkdir -p "$FOAM_RUN"
cp -r "$FOAM_TUTORIALS/incompressibleFluid/pitzDailySteady" "$FOAM_RUN"/
cd "$FOAM_RUN/pitzDailySteady"
ls 0 constant system
foamDictionary system/controlDict -entry solver -value
foamDictionary constant/physicalProperties
foamDictionary constant/momentumTransport

Secondo la User Guide v14 (§2.1), il tutorial esegue il modulo incompressibleFluid tramite foamRun; physicalProperties contiene viscosityModel constant; e nu 1e-05;; momentumTransport contiene simulationType RAS; con model kEpsilon; e turbulence on; (la guida mostra anche una voce viscosityModel Newtonian; nelle impostazioni RAS, la parte laminare dello sforzo, e un commento che elenca altri modelli provati su questo caso). 0/ contiene U, p, k, epsilon, nut e anche omega, così che kOmegaSST possa essere selezionato senza creare un nuovo file.

La geometria è il gradino all’indietro di Pitz–Daily, con un’altezza del gradino di 25,4 mm e una velocità d’ingresso di 10 m/s. Con ν = 1 × 10⁻⁵ m²/s il numero di Reynolds basato sull’altezza del gradino è 10 × 0,0254 / 10⁻⁵ ≈ 25 400: chiaramente turbolento, ed è per questo che il tutorial usa RAS e non laminare.

Risposta al compito 4 (una formulazione possibile). Separare il materiale dal modello dello sforzo significa poter cambiare una viscosità senza toccare la chiusura, e poter cambiare una chiusura senza riscrivere i dati del materiale. Documenta anche il caso: chiunque vede subito quale parte è una proprietà del fluido e quale parte è una scelta di modellazione.

Molto materiale in rete si riferisce a release Foundation più vecchie o al fork ESI (openfoam.com), i cui file hanno lo stesso scopo ma nomi e parole chiave diversi. La corrispondenza qui sotto è una guida per leggere questi esempi; riscrivili sempre nella forma v14 e verificali con gli strumenti installati.

Le note di rilascio di v14 dicono anche che il suffisso Coeffs è stato rimosso dai nomi dei sottodizionari dei coefficienti dei modelli (con retrocompatibilità). Se una voce copiata non viene riconosciuta, spesso OpenFOAM non fallisce: una parola chiave sconosciuta in un dizionario viene semplicemente ignorata e il modello gira con i suoi valori predefiniti. È l’esito più pericoloso, perché nulla ti avverte. Usa #print (sezione 14) per vedere che cosa ha effettivamente letto il modello.

simulationType seleziona uno di tre rami: laminar, RAS o LES. All’interno del ramo laminare c’è di nuovo una selezione di modello, fatta tramite il sottodizionario facoltativo laminar { model …; }. Quando è assente si usa il modello laminare newtoniano predefinito (nel codice sorgente Foundation questa classe predefinita si chiama tradizionalmente Stokes, dal nome della legge lineare dello sforzo viscoso; fa fede il nome stampato nel log e da foamToC nella tua installazione). Gli altri modelli laminari sono quelli newtoniani generalizzati e visco-elastici della sezione 16.

Le equazioni risolte da incompressibleFluid in questa configurazione sono la continuità e l’equazione della quantità di moto con lo sforzo 2νD e ν preso da physicalProperties. Non servono i campi k, epsilon, omega o nut, e se sono presenti in 0/ semplicemente non vengono letti.

Per elencare i modelli laminari disponibili nella tua installazione:

foamToC -table laminarincompressibleMomentumTransportModel

Esempio numerico per il riquadro “Perché è importante”. Per un flusso completamente sviluppato tra lastre piane parallele distanti H, con velocità media U_b, il gradiente di pressione cinematica è −dp_k/dx = 12 ν U_b / H². Con ν = 1 × 10⁻⁵ m²/s, U_b = 0,1 m/s, H = 0,01 m e una lunghezza di 0,1 m: Δp_k = 12 × 10⁻⁵ × 0,1 / 10⁻⁴ × 0,1 = 0,012 m²/s² (12 Pa per un fluido con ρ = 1000 kg/m³). La caduta è interamente viscosa e proporzionale a ν — ma solo nel tratto completamente sviluppato del canale (vedi la sezione 20).

Scomponi la velocità in una media e in una fluttuazione e media l’equazione della quantità di moto. Il termine convettivo non lineare produce una nuova incognita, il tensore degli sforzi di Reynolds ⟨u′u′⟩: sei componenti incognite aggiuntive senza alcuna equazione in più. È il problema della chiusura.

U = <U> + u'                         (decomposizione di Reynolds)

d<U>/dt + div(<U><U>) = -grad(<p>) + div[ 2 nu D(<U>) ] - div( <u'u'> )

-<u'u'> = 2 nut D(<U>) - (2/3) k I   (ipotesi di Boussinesq)

nuEff = nu + nut

L’ipotesi di Boussinesq modella gli sforzi di Reynolds come uno sforzo viscoso aggiuntivo con una viscosità turbolenta ν_t, che è una proprietà del flusso e non del fluido: varia nello spazio e nel tempo, può essere centinaia di volte più grande di ν nel nucleo turbolento e tende a zero a parete. L’equazione della quantità di moto ha allora la stessa forma di quella laminare con ν sostituita da ν_eff = ν + ν_t. Il compito del modello RAS è calcolare ν_t, di solito a partire da uno o due scalari trasportati aggiuntivi (k ed ε, k e ω, oppure ν̃).

In OpenFOAM questo si vede nei nomi dei campi: nut è ν_t, memorizzata come campo con le proprie condizioni al contorno; la viscosità efficace nuEff compare nei nomi di alcuni termini di fvSchemes. LES usa la stessa idea a una scala diversa: le equazioni filtrate contengono uno sforzo di sottogriglia modellato con una viscosità di sottogriglia, e i vortici grandi vengono risolti nel tempo invece di essere eliminati dalla media.

Un numero di Reynolds ha senso solo insieme alla sua definizione. Per i flussi interni la scelta usuale è la velocità media (bulk) e il diametro idraulico D_h = 4A/P (A area della sezione, P perimetro bagnato). Per un tubo circolare D_h = D; per due lastre piane infinite distanti H, D_h = 2H; per il condotto 10 mm × 5 mm della parte 2, D_h = 4 × 50 / 30 ≈ 6,67 mm. Usare l’altezza del canale invece di D_h cambia Re di un fattore due per le lastre: il Re = 2000 basato su H usato qui diventerebbe Re = 4000 basato su D_h.

I valori di riferimento classici quindi non sono intercambiabili: circa 2300 per un tubo (velocità media, diametro) come limite inferiore usuale della transizione; per il flusso in canale piano la teoria della stabilità lineare dà circa 5772 basato sulla velocità in asse e sulla semi-altezza, mentre negli esperimenti la transizione può avvenire a valori notevolmente più bassi in presenza di disturbi. I disturbi all’ingresso, la rugosità e la forma del profilo d’ingresso spostano tutti la soglia effettiva.

Completamente sviluppato è un’altra domanda. Re non dice nulla sulla distanza necessaria perché il profilo smetta di cambiare. Per il flusso interno laminare la lunghezza d’ingresso cresce con Re; una stima comune per i tubi è L_e ≈ 0,05 Re D_h. A Re = 400 e D_h = 8 mm questo dà già circa 0,16 m: un canale corto potrebbe non essere mai completamente sviluppato. Verificalo sulla soluzione confrontando i profili di velocità in due stazioni, non con Re.

Risposta al compito 3. A basso Re, con un flusso fisicamente laminare, un modello laminare non è una semplificazione ma la fisica corretta: non ci sono fluttuazioni turbolente da modellare, la soluzione può essere confrontata con una soluzione analitica (Poiseuille), e aggiungere un modello RAS introdurrebbe una viscosità turbolenta che non esiste, campi in più e incertezza in più.

Il modello k–ε standard trasporta l’energia cinetica turbolenta k (m²/s²) e il suo tasso di dissipazione ε (m²/s³) e costruisce da essi la viscosità turbolenta:

nut = Cmu * k^2 / epsilon

Dk/Dt       = div[ (nu + nut/sigmak)   grad(k) ]       + G - epsilon
Depsilon/Dt = div[ (nu + nut/sigmaEps) grad(epsilon) ] + C1*G*epsilon/k - C2*epsilon^2/k

G = nut * 2 |symm(grad(U))|^2

costanti standard: Cmu = 0.09, C1 = 1.44, C2 = 1.92, sigmak = 1.0, sigmaEps = 1.3

G è la produzione di k da parte del taglio medio. Le costanti sono quelle del modello standard; in OpenFOAM sono valori predefiniti nell’implementazione e vengono stampate con #print (sezione 14).

foamToC -table RASincompressibleMomentumTransportModel stampa i nomi registrati in quella tabella di selezione a run-time. Secondo la User Guide v14 l’elenco RAS incomprimibile comprende LRR, LamBremhorstKE, LaunderSharmaKE, LienCubicKE, LienLeschziner, RNGkEpsilon, SSG, ShihQuadraticKE, SpalartAllmaras, kEpsilon, kEpsilonLopesdaCosta, kOmega, kOmega2006, kOmegaSST, kOmegaSSTLM, kOmegaSSTSAS, kkLOmega, qZeta, realizableKE e v2f. L’autorità finale è la tua installazione. foamSearch "$FOAM_TUTORIALS" momentumTransport RAS/model invece scorre tutti i file momentumTransport dei tutorial e stampa i valori distinti della voce RAS/model: ti dice quali modelli sono effettivamente usati negli esempi ufficiali.

turbulence off. Con turbulence off; le equazioni del modello non vengono risolte, quindi k, ε e ν_t conservano i valori che hanno. È uno strumento per procedure particolari (per esempio l’inizializzazione di un campo di moto), non un modo per ottenere una simulazione laminare: per quello si usa simulationType laminar.

Ogni nuovo file deve avere le dimensioni corrette e una condizione al contorno per ogni patch della mesh:

In pitzDailySteady, come descritto dalla guida v14: k ed epsilon hanno valori interni e d’ingresso pari a 0,375 e 14,855; sulle pareti k usa kqRWallFunction, epsilon usa epsilonWallFunction e nut usa nutkWallFunction; nut è calculated su ingresso e uscita; frontAndBack è empty per tutti i campi perché il caso è bidimensionale. Leggi nella tua copia le condizioni d’uscita di k ed epsilon (lo schema usuale è una condizione di efflusso che passa a valore imposto in caso di riflusso).

Le note di rilascio di v14 riportano anche che le condizioni al contorno di turbolenza all’ingresso sono state semplificate in turbulentKineticEnergy, turbulentEpsilon e turbulentOmega, che, come suggeriscono i loro nomi e le condizioni che sostituiscono, ricavano i valori d’ingresso da ipotesi sull’intensità di turbolenza e sulla scala di lunghezza (o lunghezza di mescolamento). Verifica i loro parametri effettivi con foamInfo turbulentKineticEnergy prima di usarle.

Un calcolo k–ε richiede inoltre voci come queste (schema illustrativo; confrontalo con la tua copia di pitzDailySteady):

// system/fvSchemes - voci aggiuntive richieste da un calcolo k-epsilon (illustrativo)
divSchemes
{
    // ... voci già presenti nel caso ...
    div(phi,k)          bounded Gauss limitedLinear 1;
    div(phi,epsilon)    bounded Gauss limitedLinear 1;
}

// system/fvSolution - i campi di turbolenza richiedono un solutore lineare
solvers
{
    // ... voci esistenti per p e U ...
    "(k|epsilon)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }
}

In un caso stazionario SIMPLE i campi di turbolenza hanno di solito anche fattori di rilassamento e controlli sui residui in fvSolution: vedi la parte 6.

k–ω SST trasporta k e il tasso di dissipazione specifico ω (1/s), legato a ε da ω = ε/(C_μ k). Combina una formulazione k–ω vicino alle pareti con un comportamento di tipo k–ε lontano da esse e limita la viscosità turbolenta nei forti gradienti di pressione avversi (il limitatore “shear-stress transport”). Il suo trattamento a parete usa omegaWallFunction su omega. Spalart–Allmaras è un modello a un’equazione per una viscosità modificata ν̃ (nuTilda), da cui si calcola ν_t; le sue condizioni al contorno sono di tipo diverso (niente k, niente ε/ω).

Trovare gli esempi. foamSearch elenca i valori distinti di una voce e quante volte compaiono, il che conferma che kOmegaSST è usato; per ottenere i percorsi effettivi dei tutorial, cerca nei file, oppure chiedi a foamInfo, che elenca anche casi di esempio:

foamToC -help
foamToC -table RAS
foamToC -scalarBCs | grep -i wallfunction
grep -rl --include=momentumTransport kOmegaSST "$FOAM_TUTORIALS"

Risposta al compito 2. epsilon viene sostituito da omega; k e nut restano ma le loro condizioni al contorno vanno ricontrollate. Risposta al compito 3: i campi in 0/ e le loro condizioni al contorno (soprattutto pareti e ingresso), fvSchemes (nuovo div(phi,omega)), fvSolution (solutore, rilassamento, controllo dei residui per omega), la mesh vicino alla parete (y+) e le osservabili che confronterai.

Un avvertimento concreto dal tutorial ufficiale. La guida v14 riporta che su pitzDailySteady il calcolo kEpsilon converge in 285 iterazioni e realizableKE in 255, mentre kOmegaSST con la stessa numerica non converge e arriva fino al tempo finale (2000); converge in 285 iterazioni solo dopo aver cambiato lo schema del gradiente di U in cellLimited. Stessa mesh, stesso caso, modello diverso: anche le impostazioni numeriche hanno dovuto cambiare.

Da dove viene la formula k = 3/2 (U I)²: con fluttuazioni isotrope u′ = v′ = w′ = U I, k = ½(u′² + v′² + w′²) = 3/2 (U I)². Le altre grandezze discendono da k e da una scala di lunghezza l:

epsilon = Cmu^0.75 * k^1.5 / l
omega   = k^0.5 / (Cmu^0.25 * l)  =  epsilon / (Cmu * k)
nut     = Cmu * k^2 / epsilon
(Cmu = 0.09)

Proseguendo l’esempio precedente (U = 2,0 m/s, I = 0,05, k = 0,015 m²/s²) con un l ipotizzato di 1 mm: ε ≈ 0,302 m²/s³, ω ≈ 224 1/s, ν_t ≈ 6,7 × 10⁻⁵ m²/s, cioè circa 67 volte la viscosità dell’acqua. Una regola spesso citata per il flusso interno completamente sviluppato prende l ≈ 0,07 D_h; è una regola empirica, da giustificare come I.

Ricostruire a ritroso il tutorial ufficiale. pitzDailySteady ha U = 10 m/s e k = 0,375 m²/s²: 1,5 (10 I)² = 0,375 dà I = 0,05 (5%). Il suo ε = 14,855 m²/s³ dà l = 0,09^0,75 × 0,375^1,5 / 14,855 ≈ 0,00254 m, cioè un decimo dell’altezza del gradino di 25,4 mm; e ω = 14,855/(0,09 × 0,375) ≈ 440,2 1/s, esattamente il valore che la guida riporta per 0/omega. I valori del tutorial sono quindi coerenti con I = 5% e l = 2,54 mm.

Risposta al compito 2. I descrive le fluttuazioni del flusso in arrivo: dipende da ciò che sta a monte (un lungo tubo, un ventilatore, un nido d’ape, una contrazione) ed è raramente misurata. Ipotizzarla fissa k all’ingresso, e tramite l anche ε o ω e quindi ν_t; lo sviluppo a valle, la transizione e la separazione possono dipenderne. Risposta al compito 3: k è una proprietà del flusso (energia delle fluttuazioni di velocità, m²/s²), ν è una proprietà del fluido (m²/s): dimensioni diverse e significato fisico diverso, quindi nessuna delle due può essere ricavata dall’altra.

Vicino a una parete liscia la velocità media, scalata con la velocità d’attrito u_τ = √(τ_w/ρ), segue un profilo universale nella variabile y⁺: un sottostrato viscoso (u⁺ ≈ y⁺, all’incirca per y⁺ < 5), uno strato di raccordo (buffer layer) e uno strato logaritmico (u⁺ = (1/κ) ln(E y⁺), con κ ≈ 0,41 ed E ≈ 9,8, all’incirca da y⁺ ≈ 30 in su). Una wall function usa questa legge per fornire lo sforzo a parete e i valori vicino alla parete delle grandezze turbolente quando la prima cella è troppo grossolana per risolvere il gradiente.

In OpenFOAM la wall function è una condizione al contorno: nutkWallFunction calcola ν_t a parete a partire da k; epsilonWallFunction (o omegaWallFunction) impone la dissipazione vicino alla parete; kqRWallFunction è in sostanza una condizione a gradiente nullo per k. La guida v14 elenca le wall function per ν_t nutWallFunction, nutRoughWallFunction, nutUSpaldingWallFunction, nutkWallFunction e nutkAtmWallFunction; le costanti E, κ e C_μ possono essere ridefinite per ogni condizione al contorno. Verifica i nomi esatti nella tua installazione con foamToC -scalarBCs | grep -i wallfunction.

Esempio numerico. Se τ_w/ρ = 0,25 m²/s², allora u_τ = 0,5 m/s; con ν = 1 × 10⁻⁵ m²/s e il centro della prima cella a y = 0,5 mm, y⁺ = 0,5 × 5 × 10⁻⁴ / 10⁻⁵ = 25. Per calcolare y⁺ su un risultato, la funzione di post-elaborazione deve poter costruire il modello di trasporto della quantità di moto, quindi va eseguita attraverso il modulo solutore:

foamPostProcess -solver incompressibleFluid -func yPlus -latestTime
foamPostProcess -solver incompressibleFluid -func wallShearStress -latestTime

Risposta al compito 3. Un 0/nut scritto per kEpsilon presuppone che k esista e che il trattamento a parete sia coerente con l’y⁺ della prima cella di quella mesh; dopo il passaggio a SpalartAllmaras (niente k) o a un trattamento a basso Re che richiede y⁺ ≈ 1, le stesse voci sono sbagliate o incoerenti anche se il caso parte. Inoltre i nomi delle patch in 0/nut devono coincidere esattamente con constant/polyMesh/boundary (compito 1), altrimenti il calcolo si ferma segnalando che non esiste una voce patchField per quella patch.

#print è una direttiva dei dizionari di v14: posta alla fine di un dizionario, fa stampare a OpenFOAM il dizionario così come viene usato, comprese le voci predefinite che non hai scritto. Secondo le note di rilascio di v14 sostituisce la function object writeOptionalEntries; nelle release più vecchie l’interruttore printCoeffs di ciascun modello svolgeva un ruolo simile per i coefficienti. Nel dizionario RAS di un caso kEpsilon rivela i coefficienti predefiniti (per esempio Cmu 0.09;); la stessa direttiva dentro un dizionario PIMPLE rivela valori predefiniti come nCorrectors 1;. Esegui il caso (basta una iterazione) e leggi il dizionario stampato all’inizio del log. La guida cita anche, come alternativa, l’espansione del file con foamDictionary.

Se una modifica è giustificata, scrivila nel dizionario RAS o in un sottodizionario che porta il nome del modello (in v14 senza il vecchio suffisso Coeffs):

RAS
{
    model           kEpsilon;
    turbulence      on;

    kEpsilon                // facoltativo, solo per una modifica documentata
    {
        Cmu         0.09;
    }
}

Poi esegui di nuovo con #print e verifica che il valore stampato sia il tuo: una modifica con una parola chiave scritta male viene ignorata in silenzio. I coefficienti non sono proprietà del fluido: cambiare C_μ ricalibra il modello, e ogni risultato ottenuto con esso deve dichiararlo.

LES applica alle equazioni un filtro spaziale di ampiezza Δ. I moti più grandi di Δ vengono calcolati; l’effetto di quelli più piccoli compare come sforzo di sottogriglia, modellato per esempio con una viscosità di sottogriglia. Nel modello di Smagorinsky ν_sgs = (C_s Δ)² |S|, con |S| il modulo della velocità di deformazione risolta; WALE, kEqn, dynamicKEqn e dynamicLagrangian sono altre scelte elencate dalla guida v14, insieme ai modelli ibridi RAS/LES SpalartAllmarasDES, SpalartAllmarasDDES, SpalartAllmarasIDDES e kOmegaSSTDES, e a DeardorffDiffStress.

In OpenFOAM l’ampiezza del filtro è a sua volta un modello di delta selezionato a run-time, indicato dalla voce delta; la guida cita come esempio cubeRootVol (Δ dalla radice cubica del volume della cella), con i suoi coefficienti in un sottodizionario che ne porta il nome. Il dizionario LES richiede quindi model, turbulence, delta e, se necessario, sottodizionari con il nome del modello e del delta. Trova combinazioni funzionanti con foamSearch "$FOAM_TUTORIALS" momentumTransport LES/model e foamSearch "$FOAM_TUTORIALS" momentumTransport LES/delta, ed elenca i modelli con foamToC -table LESincompressibleMomentumTransportModel.

Perché è onerosa: LES è sempre tridimensionale e instazionaria; il passo temporale deve mantenere il numero di Courant abbastanza basso da seguire i vortici; i contorni d’ingresso devono fornire fluttuazioni realistiche; i risultati sono statistiche, accumulate su molti tempi di attraversamento (per esempio con una funzione di media dei campi); vicino alle pareti i vortici che trasportano energia diventano molto piccoli, quindi una LES risolta a parete richiede una mesh molto più fine di una mesh RAS. È per questo che il corso mantiene il canale laminare come oggetto di apprendimento.

Un modello newtoniano generalizzato mantiene la forma τ = 2ν(γ̇)D ma lascia dipendere la viscosità dalla velocità di deformazione locale γ̇. La guida v14 (§8.3) descrive queste leggi; la viscosità a velocità di deformazione nulla ν₀ è il nu di physicalProperties:

velocità di deformazione:  gammaDot = sqrt(2) * |symm(grad(U))|

BirdCarreau:     nu = nuInf + (nu0 - nuInf) * [1 + (k*gammaDot)^a]^((n - 1)/a)     (a = 2 per default)
CrossPowerLaw:   nu = nuInf + (nu0 - nuInf) / [1 + (m*gammaDot)^n]
powerLaw:        nu = k * gammaDot^(n - 1),          limitata a nuMin <= nu <= nuMax
HerschelBulkley: nu = min( nu0, (tau0 + k*gammaDot^n) / gammaDot )

Nota che k e n compaiono in più modelli con significati e dimensioni diversi (una costante di tempo in BirdCarreau, un indice di consistenza in powerLaw e HerschelBulkley): è esattamente la trappola segnalata sopra. Un ramo completo a scopo illustrativo, con numeri arbitrari che non descrivono un materiale reale:

// constant/physicalProperties: nu è la viscosità a velocità di deformazione nulla nu0
viscosityModel  constant;
nu              1e-03;

// constant/momentumTransport (valori illustrativi, non un materiale reale)
simulationType  laminar;

laminar
{
    model           generalisedNewtonian;
    viscosityModel  BirdCarreau;

    nuInf           1e-05;
    k               1;
    n               0.5;
}

Oltre ai fluidi newtoniani generalizzati, il ramo laminare contiene anche modelli visco-elastici (Maxwell, Giesekus, PTT, con parametri come nuM e lambda) e il modello lambdaThixotropic; le note di rilascio di v14 aggiungono la possibilità di specificare i modelli visco-elastici con più tempi di rilassamento e le corrispondenti viscosità. Le leggi di viscosità disponibili sono elencate da foamToC -table generalisedNewtonianViscosityModel.

La stessa voce viscosityModel esiste anche nei dizionari RAS e LES, dove descrive la parte laminare (molecolare) dello sforzo: per default è Newtonian, ed è per questo che la guida v14 mostra viscosityModel Newtonian; nelle impostazioni RAS di pitzDailySteady e osserva che potrebbe essere omessa. Reologia e chiusura della turbolenza sono quindi due scelte indipendenti che possono essere combinate.

OpenFOAM costruisce i modelli per nome durante l’esecuzione: ogni famiglia ha una tabella di selezione a run-time che associa una parola (kEpsilon) a una classe. Il nome della tabella codifica la famiglia — RAS + incompressible + MomentumTransportModel — ed è per questo che un modello comprimibile si trova in una tabella diversa (la guida mostra, per esempio, buoyantKEpsilon solo tra i modelli RAS comprimibili). foamToC -help elenca le altre opzioni dello strumento; foamToC -table RAS è la forma che la guida usa per cercare le tabelle il cui nome contiene RAS.

Se scrivi un nome che non è nella tabella, il calcolo si ferma con un errore fatale che dichiara sconosciuto il tipo, seguito dall’elenco dei tipi validi (la formulazione esatta dipende dalla build): il messaggio d’errore è esso stesso un elenco alla foamToC, filtrato dal modulo effettivamente in esecuzione. È la risposta definitiva alla domanda “questo modulo supporta questo modello?”.

foamInfo kEpsilon stampa la descrizione presa dall’intestazione del sorgente, la posizione dei file sorgente e un elenco di tutorial che usano il modello: la risposta più rapida al compito 3. foamSearch lavora sui file di un albero di directory ed è lo strumento giusto per vedere quali valori si usano in pratica.

In un modulo comprimibile o termico (per esempio fluid), physicalProperties descrive un intero pacchetto termofisico. Ogni voce di thermoType seleziona un ingrediente; le opzioni qui sotto sono quelle elencate dalla guida v14 (§8.1):

Il pacchetto ha poi bisogno dei dati per gli ingredienti scelti, raggruppati in un sottodizionario: massa molare, capacità termica e calore di formazione, viscosità dinamica e numero di Prandtl per il trasporto const. Una struttura illustrativa per l’aria:

mixture
{
    specie
    {
        molWeight       28.96;
    }
    thermodynamics
    {
        Cp              1004.5;
        Hf              0;
    }
    transport
    {
        mu              1.8e-05;
        Pr              0.7;
    }
}

Nota che qui la proprietà di trasporto è la viscosità dinamica mu (Pa·s), non nu: in un modulo comprimibile la densità è un campo, quindi la scorciatoia cinematica della sezione 3 non è più disponibile. Con transport sutherland, anche μ diventa una funzione della temperatura. Un caso comprimibile nuovo richiede inoltre un campo di temperatura 0/T, una pressione in pascal e condizioni al contorno termiche.

Il modo più rapido per vedere che cosa tocca un cambiamento di fisica è confrontarsi con un caso che lo contiene già. Copia un caso laminare e un tutorial RAS e confrontane la struttura:

ls 0 constant system
diff <(ls ../laminarCase/0) <(ls 0)
foamDictionary system/fvSchemes -entry divSchemes
foamDictionary system/fvSolution -entry solvers -keywords

-keywords stampa solo le parole chiave di un dizionario, il che mostra a colpo d’occhio quali campi hanno un solutore lineare. Rispondi per iscritto alle quattro domande della regola d’oro per ogni riga della tabella. Per laminare → kEpsilon, per esempio: (1) si aggiungono due equazioni di trasporto per k ed ε e l’equazione della quantità di moto usa ν + ν_t; (2) si aggiungono le incognite k, ε e il campo derivato ν_t; (3) valori iniziali e d’ingresso di k ed ε, wall function su ogni parete, condizioni d’uscita; (4) div(phi,k), div(phi,epsilon), i solutori e, per SIMPLE, rilassamento e controlli sui residui per i nuovi campi.

Per “Newtoniano → BirdCarreau” è importante la voce “output/validazione” della tabella: con un fluido fluidificante la soluzione analitica di Poiseuille non vale più, quindi il confronto che hai usato per validare il canale laminare va sostituito con un riferimento appropriato.

cp -a copia l’intera directory conservando date e permessi. Se channelBaseline contiene già cartelle temporali di risultati, rimuovile in ogni copia (foamListTimes -rm, eseguito dentro la copia) così che i tre calcoli partano dalla stessa 0/. Poi imposta i valori con foamDictionary e verifica:

foamDictionary channel_nu_half/constant/physicalProperties   -entry nu -set 5e-06
foamDictionary channel_nu_double/constant/physicalProperties -entry nu -set 2e-05

for c in channel_nu_base channel_nu_half channel_nu_double
do
    printf "%-20s nu = %s\n" "$c" "$(foamDictionary $c/constant/physicalProperties -entry nu -value)"
done

diff -r channel_nu_base/constant channel_nu_half/constant

L’output di diff -r deve mostrare solo la riga di nu. Qualsiasi altra differenza è un cambiamento non controllato e invalida il confronto.

Che cosa prevede la teoria, e quando. Nella regione completamente sviluppata di un canale laminare a portata fissata, il gradiente di pressione cinematica è proporzionale a ν (12 ν U_b/H² per le lastre, l’analoga legge di Hagen–Poiseuille per un tubo), quindi lì HALF e DOUBLE dovrebbero dare circa la metà e il doppio del gradiente di BASE. La differenza totale tra ingresso e uscita contiene però anche la regione di sviluppo, dove la perdita di pressione in eccesso dipende da Re: dimezzando ν, Re raddoppia e la lunghezza d’ingresso cresce, quindi il rapporto totale può discostarsi da 2. Se l’ingresso fosse una differenza di pressione imposta invece di una velocità imposta, cambierebbe invece la portata. Sono queste le ragioni dell’avvertenza sull’interpretazione.

Per prima cosa rendi tracciabile il calcolo stesso: foamRun | tee cattura solo lo standard output, quindi aggiungi 2>&1 per conservare nel log anche i messaggi d’errore, e verifica fin dove è arrivato il calcolo:

foamListTimes -rm                # rimuove i risultati precedenti, conserva 0/
foamRun 2>&1 | tee log.foamRun
grep -c "^Time =" log.foamRun    # quante iterazioni/passi temporali sono stati eseguiti
tail -n 30 log.foamRun

Poi misura sulle patch, non in celle arbitrarie. Le funzioni di post-elaborazione predefinite di v14 comprendono patchAverage (media dei campi su una patch), patchDifference (differenza tra le medie su due patch), patchFlowRate e patchFlux (flusso attraverso una patch). La portata attraverso ingresso e uscita deve essere uguale e di segno opposto: è il tuo controllo del flusso. Per la pressione, la media pesata sull’area di p su inlet meno quella su outlet è una Δp_k ben definita (moltiplica per ρ per ottenere pascal).

cd channel_nu_base
foamListTimes                    # quali cartelle temporali esistono?
foamPostProcess -func "patchFlowRate(name=inlet)"  -latestTime
foamPostProcess -func "patchFlowRate(name=outlet)" -latestTime
foamGet patchAverage             # copia il modello in system/ e leggine gli argomenti

La guida v14 scrive l’argomento della patch come name= sulla riga di comando e come patch= in un esempio con #includeFunc, quindi leggi l’elenco degli argomenti nell’intestazione del modello che foamGet copia in system/ prima di eseguire patchAverage o patchDifference. Usa lo stesso istante (-latestTime di un calcolo convergente, o lo stesso tempo finale) per tutti e tre i casi, e riporta i comandi stessi: questa è la “procedura effettiva di misura” del compito 3.

Plausibilità (compito 4). Re si dimezza quando ν raddoppia (compito 1). Confronta i rapporti Δp_half/Δp_base e Δp_double/Δp_base con 0,5 e 2: un rapporto vicino a questi valori, con il controllo del flusso soddisfatto, è fisicamente plausibile per un canale laminare a portata fissata; uno scostamento grande richiede una spiegazione (flusso in sviluppo, calcolo non convergente, mesh troppo grossolana per gli strati limite più sottili del calcolo a ν minore).

Copia il tutorial (vedi i comandi nell’approfondimento della sezione 5) e compila ogni voce con una riga effettivamente letta dalla copia, non a memoria. Uno schema delle dipendenze può essere semplice come questo schizzo testuale; la tua versione deve contenere i valori e i nomi che hai trovato davvero:

physicalProperties : viscosityModel constant; nu 1e-05
        |
        v
momentumTransport  : simulationType RAS; model kEpsilon
        |
        +--> 0/k, 0/epsilon, 0/nut   (+ 0/U, 0/p)
        |        |
        |        +--> patch di parete: kqRWallFunction, epsilonWallFunction, nutkWallFunction
        |        +--> ingresso: valori di k ed epsilon (ipotesi su I e l)
        |
        +--> fvSchemes : div(phi,k), div(phi,epsilon)
        +--> fvSolution: solutori per k, epsilon; residualControl

Risultati attesi, secondo la guida v14: viscosityModel constant, nu 1e-05; simulationType RAS, model kEpsilon, turbulence on; campi k, epsilon, nut (e omega predisposto per kOmegaSST); kqRWallFunction, epsilonWallFunction, nutkWallFunction su upperWall e lowerWall; all’ingresso k = 0,375 m²/s² ed ε = 14,855 m²/s³ (I = 5% e l ≈ 2,54 mm, come ricavato nella sezione 12); SIMPLE con residualControl per p, U e i campi di turbolenza. I solutori lineari e gli schemi di divergenza vanno letti dalla tua copia.

Perché “laminare → RAS” non è una parola sola. Il caso laminare non ha campi k, ε, ν_t, non ha wall function, né dati di turbolenza all’ingresso, né schemi o solutori per questi campi; e la sua mesh è stata progettata per un profilo laminare, non per l’y⁺ di una wall function. Cambiare solo simulationType porterebbe a un arresto al primo file mancante — o, peggio, a un calcolo con dati incoerenti.

Le prime righe dopo --> FOAM FATAL ERROR o --> FOAM FATAL IO ERROR indicano quasi sempre il file e la parola chiave. Schemi tipici, uno per riga della tabella (la formulazione esatta può variare leggermente tra le build):

--> FOAM FATAL ERROR: cannot find file ".../0/k"

--> FOAM FATAL IO ERROR: keyword div(phi,k) is undefined in dictionary ".../system/fvSchemes/divSchemes"

Unknown ... type kEpsilonn
Valid ... types are: ...   (list of the names in the selection table)

Cannot find patchField entry for walls

Leggi l’errore nel suo contesto, risalendo: le righe immediatamente sopra nel log dicono che cosa stava costruendo OpenFOAM (il modello di trasporto della quantità di moto, un campo, una matrice) quando è fallito. Annota anche il percorso del file e, per gli errori di IO, il numero di riga. Solo allora apri il file indicato. Per la “convergenza scarsa” non esiste un messaggio unico: traccia i residui (foamMonitor sul file dei residui, se ne hai uno, oppure grep sul log), guarda in quale zona del dominio la soluzione si comporta male, e cambia una cosa alla volta — prima le condizioni al contorno e la mesh, per ultimi i controlli numerici.

Tieni un breve registro di ogni tentativo (che cosa hai cambiato, che cosa è successo). Con il metodo di questa sezione il registro è la prova; con le correzioni a caso, la prova va persa.

Per elencare i casi candidati: ls "$FOAM_TUTORIALS/incompressibleFluid", poi controlla ciascuno con foamDictionary <caso>/constant/momentumTransport -entry simulationType -value. Uno scheletro da compilare, una riga per ogni voce dell’elaborato:

Identificazione caso : $FOAM_TUTORIALS/incompressibleFluid/<caso>, OpenFOAM 14, solutore incompressibleFluid
Materiale            : constant/physicalProperties -> viscosityModel ..., nu ... [m^2/s]
Modello di quantità  : simulationType ...; model ...
di moto
Incognite            : 0/U 0/p 0/... (motivo per ciascuno)
Vicino a parete      : patch ... -> nut ..., k ..., epsilon/omega ...
Numerica             : fvSchemes div(phi,...) ; fvSolution solvers "(...)"
Evidenze             : User Guide v14 sezione ..., percorso del tutorial ...

Per la riga “Evidenze”, una pagina di riferimento ufficiale è una sezione della User Guide v14 (per esempio §8.2 per il modello di turbolenza, §2.1 se hai scelto pitzDailySteady), e un esempio installato è il percorso di un tutorial con la voce che citi. Per il compito facoltativo, la coppia naturale è kEpsilon contro kOmegaSST su pitzDailySteady: epsilon contro omega, epsilonWallFunction contro omegaWallFunction, le voci aggiuntive di schemi e solutori, e la differenza documentata nel comportamento di convergenza (sezione 11) — che è un’osservazione numerica, non un’affermazione su quale modello sia più accurato.

ν. 1,8 × 10⁻⁵ / 1,2 = 1,5 × 10⁻⁵ m²/s — la viscosità cinematica dell’aria a temperatura ambiente, lo stesso valore usato da molti tutorial di OpenFOAM. Re. 0,05 × 0,008 / 10⁻⁶ = 4 × 10⁻⁴ / 10⁻⁶ = 400. k. U I = 1,5 × 0,03 = 0,045 m/s; 0,045² = 0,002025; × 1,5 = 0,0030375 m²/s².

Sezione 3, compito 3: vedi l’approfondimento “La legge newtoniana e perché il solutore incomprimibile non ha bisogno della densità”. Sezione 5: foamVersion deve indicare 14; il valore del tutorial è nu 1e-05; simulationType RAS. Sezione 11: epsilon viene sostituito da omega. Sezione 13: i nomi delle patch devono essere identici in 0/nut e in constant/polyMesh/boundary; le wall function si trovano in 0/nut, 0/k e 0/epsilon (o 0/omega). Sezione 17: gli output delle tabelle di foamToC e di foamInfo dipendono dalla tua build; almeno kEpsilon e kOmegaSST sono nella tabella RAS incomprimibile, e pitzDailySteady dimostra kEpsilon.

Esperimento A, una verifica numerica incrociata. Se il canale BASE avesse H = 10 mm, U_b = 0,1 m/s e ν = 1 × 10⁻⁵ m²/s, allora Re_H = 100; HALF dà Re_H = 200 e DOUBLE Re_H = 50. I gradienti di pressione completamente sviluppati sarebbero 0,12, 0,06 e 0,24 m²/s² al metro. Usa i valori del tuo caso: questi numeri mostrano solo l’andamento atteso.

Per qualsiasi caso incomprimibile sconosciuto, questi comandi rispondono alle cinque domande sopra in meno di un minuto, e sono una buona abitudine da portare nella parte 8:

foamDictionary system/controlDict -entry solver -value          # quale modulo
foamDictionary constant/physicalProperties                        # (1) materiale
foamDictionary constant/momentumTransport                         # (2) modello
ls 0                                                              # (3) campi
foamDictionary 0/nut -entry boundaryField                         # (4) trattamento a parete/ingresso
foamInfo kEpsilon                                                 # (5) documentazione ed esempi

Sostituisci 0/nut e kEpsilon con i campi e il modello che hai effettivamente trovato. Il punto non sono i comandi in sé ma l’ordine delle domande: materiale, modello, incognite, contorni, evidenze.