Una bella immagine e tre domande
Uno studente ti mostra una splendida mappa del flusso attorno a un’auto. Colori morbidi, linee di flusso plausibili, un coefficiente di resistenza con tre decimali. Fai tre domande: quanto valeva y+ sulla carrozzeria, il bilancio di massa tornava, e di quanto è cambiata la resistenza raffinando la mesh? Nessuna ha una risposta.
Quello studente ha imparato un software di CFD. Non ha ancora imparato la CFD.
Per imparare la CFD da zero, impara a distinguere una risposta giusta da una risposta sicura di sé. In pratica sono quattro fasi, in quest’ordine: la poca matematica che spiega perché le soluzioni numeriche sbagliano, un piccolo solutore che scrivi tu e rompi apposta, una scala di casi di riferimento con risposte pubblicate, e solo allora un software commerciale o open source. Questo articolo descrive il percorso, con un esempio di calcolo che mostra in dieci righe di algebra la lezione più importante della materia, e un piano settimana per settimana alla fine.
Che cosa significa imparare la CFD?
Tre cose, e il software è una quarta che conta meno di tutte.
- La fisica. Che cosa dicono le equazioni di conservazione, che cosa fa ogni termine e quali termini dominano nel tuo problema. Un numero di Reynolds, uno di Péclet e uno di Mach dicono di un flusso, prima della mesh, più di quanto dica la maggior parte delle mappe a colori dopo il calcolo.
- La numerica. Come un’equazione continua diventa algebra, e quale errore introduce. Ogni risposta ai volumi finiti è sbagliata di un errore di discretizzazione, e l’abilità sta nel conoscerne l’entità e la natura.
- Il giudizio. Decidere che cosa modellare, verificare che i numeri siano convergenti, bilanciati e poco sensibili alla mesh, e confrontarli con qualcosa di indipendente. Sono la verifica e la validazione, ed è per questo che pagano datori di lavoro e clienti.
- Il software è l’interfaccia verso tutte e tre. Si impara in poche settimane quando il resto è a posto, e comunque cambia a ogni rilascio.
La matematica che serve, e perché
Meno di quanto molti temano, e più di quanto ammetta la maggior parte dei tutorial.
Serve il calcolo vettoriale (divergenza, gradiente, teorema della divergenza), perché il metodo dei volumi finiti è il teorema della divergenza applicato a una piccola cella. Serve riconoscere le tre classi di equazioni alle derivate parziali, perché si comportano in modo diverso: i problemi ellittici (diffusione stazionaria, pressione) sentono tutti i contorni insieme, quelli parabolici (diffusione transitoria) avanzano nel tempo, e quelli iperbolici (convezione pura, onde, urti) trasportano l’informazione lungo le caratteristiche e puniscono gli schemi che ignorano la direzione del flusso.
E servono le serie di Taylor, perché è da lì che nasce l’errore di discretizzazione. Sviluppa una funzione regolare attorno a un punto e fai le differenze:
Qui D_c è la differenza centrata e D_u la differenza upwind. La differenza centrata è del secondo ordine: dimezza la cella e l’errore principale scende a un quarto. La differenza upwind è del primo ordine, e il suo termine di errore principale è proporzionale a una derivata seconda. Tieni a mente questo punto. È tutta la sezione successiva.
Scrivi un piccolo solutore e rompilo
Prima di aprire un’interfaccia commerciale, scrivi un solutore di una quarantina di righe. Imparerai di più guardandolo sbagliare che in un mese di tutorial.
Un punto di partenza gratuito e ben fatto è CFD Python, i «12 passi verso Navier-Stokes» di Lorena Barba e Gilbert Forsyth: notebook Jupyter che vanno dalla convezione lineare 1D alla diffusione, all’equazione di Burgers, alle equazioni di Laplace e di Poisson fino al flusso 2D in cavità e in canale. Il contenuto è sotto licenza CC-BY 4.0 e il codice sotto BSD-3, quindi puoi usarli e adattarli liberamente.
Poi fai l’unico esperimento che ogni ingegnere CFD dovrebbe aver fatto con le proprie mani.
Esempio di calcolo: il numero di Péclet di cella
Prendiamo la convezione-diffusione stazionaria e monodimensionale di uno scalare φ su 0 ≤ x ≤ L, con φ(0) = 0 e φ(L) = 1. L’equazione e la sua soluzione esatta sono:
Con Pe = 50 la soluzione è quasi nulla su tutto il dominio e sale bruscamente in uno strato sottile all’estremità destra. In x = 0,9 L il valore esatto è circa exp(−5) ≈ 0,007.
Discretizza con N celle uniformi. La grandezza che governa tutto è il numero di Péclet di cella, Pe_Δ = ρuΔx / Γ = Pe / N. Su una griglia uniforme ogni schema si riduce a una ricorrenza le cui soluzioni crescono da nodo a nodo di un rapporto fisso, e quel rapporto dice che cosa fa lo schema:
La prima riga è il rapporto esatto tra nodi vicini, r_c il rapporto prodotto dallo schema centrato e r_u quello prodotto dallo schema upwind del primo ordine.
Su 10 celle, Pe_Δ = 5:
| Schema | Rapporto tra nodi vicini | Che cosa vedi |
|---|---|---|
| Esatta | exp(5) ≈ 148 | Uno strato ripido nell’ultima cella |
| Centrato | (1 + 2.5)/(1 − 2.5) ≈ −2.33 | Il segno cambia da nodo a nodo: oscillazioni e valori fuori dall’intervallo da 0 a 1 |
| Upwind | 1 + 5 = 6 | Liscio, monotono e con uno strato fin troppo spesso |
Su 40 celle, Pe_Δ = 1,25: lo schema centrato dà un rapporto di circa 4,3 contro 3,5 esatto, liscio e vicino; l’upwind dà 2,25, ancora visibilmente spalmato.
Ne seguono due lezioni, e sono le più importanti di questo articolo.
Lo schema centrato oscilla quando Pe_Δ > 2. Il rapporto diventa negativo e la soluzione presenta oscillazioni che non sono fisica. È il risultato classico del manuale di Patankar, ed è il motivo per cui i solutori commerciali non usano la differenza centrata pura per la convezione nella maggior parte dei flussi.
L’upwind non oscilla mai e sbaglia in un modo che sembra giusto. L’errore di troncamento dello sviluppo di Taylor sopra, (Δx/2)φ″, agisce esattamente come una diffusione aggiuntiva. Lo schema si comporta come se la diffusività fosse:
Su 10 celle fa 3,5 Γ. La tua bella soluzione upwind è, al primo ordine, la soluzione di un problema con un numero di Péclet di circa 50/3,5 ≈ 14, non 50. Niente sullo schermo te lo dice. Le isolinee sono pulite, i residui scendono splendidamente, e la risposta appartiene a un altro flusso.
Per questo preferisco che un principiante veda prima una soluzione centrata che oscilla. Le oscillazioni sono oneste.
Ogni solutore moderno offre schemi limitati di ordine superiore (upwind del secondo ordine con limitatori, QUICK, famiglie TVD) che stanno tra i due. Capirai che cosa fanno queste impostazioni solo dopo aver visto fallire i due estremi.
Controlla l’ordine di accuratezza
Ora raffina. Calcola l’errore rispetto alla soluzione esatta su due griglie con rapporto di raffinamento r. Con e₁ l’errore sulla griglia grossolana ed e₂ quello sulla griglia fine, l’ordine osservato è:
Esempio di calcolo, numeri illustrativi. Uno schema del secondo ordine dà errori di 0,020 e 0,0051 su griglie con passo Δx e Δx/2, quindi p = ln(3,92)/ln(2) ≈ 1,97. Corretto. Se ottieni 1,0 da uno schema che dovrebbe essere del secondo ordine, qualcosa da qualche parte è del primo ordine: una condizione al contorno, un gradiente, un’interpolazione sulle facce. Trovarlo è la migliore formazione numerica che esista. È lo stesso calcolo che sta dentro l’indice di convergenza di griglia, spiegato passo per passo nel nostro articolo sull’indipendenza dalla mesh.
Sali una scala di casi con risposte pubblicate
Quando il tuo codice si comporta bene, passa a casi di cui qualcun altro ha pubblicato la risposta. Ogni gradino insegna una nuova difficoltà, e ognuno ha un controllo che si supera o no.
| Gradino | Caso | Che cosa insegna | Riferimento | Il controllo |
|---|---|---|---|---|
| 1 | Convezione-diffusione 1D | Schemi, numero di Péclet, ordine di accuratezza | La soluzione esatta sopra | L’ordine osservato corrisponde allo schema |
| 2 | Cavità con coperchio mobile, Re 100 e 1000 | Accoppiamento pressione-velocità in 2D | Ghia, Ghia e Shin (1982), velocità sulle linee mediane | I tuoi profili passano sui loro punti tabulati |
| 3 | Flusso laminare in tubo | Lunghezza di ingresso, gradienti di parete | Hagen-Poiseuille: f = 64/Re | Il fattore d’attrito tende a 64/Re raffinando la mesh |
| 4 | Gradino all’indietro, laminare | Separazione e riattacco | Armaly et al. (1983), misure | Lunghezza di riattacco in funzione del numero di Reynolds |
| 5 | Cilindro in un canale, DFG 2D-2 | Flusso transitorio, passo temporale | FeatFlow: C_D,max = 3,2200, St = 0,30188 | Compare il distacco di vortici ed entrambi i valori si avvicinano |
| 6 | Canale turbolento | Trattamento di parete e limiti dei RANS | DNS a Re_τ da 180 a 5200 (Lee e Moser, 2015; dati dell’Oden Institute) | Velocità media in unità di parete rispetto al profilo DNS |
| 7 | Tubo riscaldato, turbolento | Coefficienti di scambio termico | La correlazione di Gnielinski | Numero di Nusselt entro la dispersione della correlazione |
Il gradino 3 merita un’osservazione. La lunghezza di ingresso laminare è circa 0,05 Re D, quindi a Re = 1000 il profilo ha bisogno di circa 50 diametri per svilupparsi. Un tubo corto dà un fattore d’attrito ben sopra 64/Re, ed è il flusso in sviluppo che fa quello che deve. Imparare a distinguere le due cose è proprio lo scopo.
Al gradino 6 la maggior parte degli autodidatti scopre che un modello di turbolenza è un modello. Confronta il profilo di velocità media in unità di parete con la DNS, una volta con la prima cella a y+ ≈ 1 e una volta con le funzioni di parete a y+ ≈ 50. L’articolo su k-omega SST e k-epsilon di questo blog spiega perché i due si comportano in modo diverso.
Il gradino 5 è anche un buon test per chiunque, o qualunque cosa, dichiari di impostare casi CFD al posto tuo; l’articolo sugli agenti di IA per la CFD lo usa proprio così.
Poi impara il software
A questo punto i menu avranno senso, perché ognuno è una decisione che hai già preso a mano: lo schema di convezione, l’accoppiamento pressione-velocità, il trattamento di parete, il passo temporale.
Impara bene un software prima di provarne diversi. Per imparare bastano due strade gratuite. OpenFOAM è open source e gira su Linux, o su Windows tramite WSL. Ansys Student è gratuito per studenti e docenti, gira su Windows e consente fino a 1 milione di celle per problemi di fluidi su un massimo di quattro core, solo per uso didattico. Un milione di celle basta per ogni gradino della scala sopra, e non è un caso.
Ripeti i gradini 2, 3 e 5 nel software. Se la tua cavità scritta a mano e quella commerciale concordano tra loro e con Ghia, le hai verificate entrambe, e hai imparato le impostazioni del software per confronto e non per fiducia.
Libri da avere
- Versteeg e Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Il primo manuale più chiaro sul metodo dei volumi finiti usato dai codici commerciali.
- Patankar, Numerical Heat Transfer and Fluid Flow. Breve, datato e ancora la migliore spiegazione del perché gli schemi si comportano come si comportano.
- Ferziger, Perić e Street, Computational Methods for Fluid Dynamics. Il riferimento per la numerica una volta superate le basi.
- Moukalled, Mangani e Darwish, The Finite Volume Method in Computational Fluid Dynamics. Il ponte tra la teoria e l’implementazione in OpenFOAM.
- Pope, Turbulent Flows. Per quando «quale modello di turbolenza?» diventa una domanda vera.
Un piano di 16 settimane
La mia proposta per chi ha una laurea in ingegneria e circa otto o dieci ore a settimana. Ogni blocco termina con qualcosa che puoi mostrare.
| Settimane | Lavoro | Risultato |
|---|---|---|
| 1 e 2 | Serie di Taylor, classi di EDP; passi 1 a 4 di CFD Python | Un codice di diffusione 1D con la sua soluzione esatta |
| 3 e 4 | Convezione-diffusione; l’esperimento del numero di Péclet; ordine di accuratezza | Un grafico delle soluzioni esatta, centrata e upwind a Pe_Δ = 5 e 1,25 |
| 5 e 6 | Passi 9 a 11 di CFD Python; cavità con coperchio mobile | La tua cavità confrontata con Ghia et al. a Re 100 |
| 7 e 8 | Installare un software; ripetere la cavità; flusso laminare in tubo | Fattore d’attrito confrontato con 64/Re, con raffinamento della mesh |
| 9 e 10 | Mesh: strati limite, metriche di qualità; uno studio su tre mesh | Un indice di convergenza di griglia per il tubo |
| 11 e 12 | Turbolenza e pareti; canale turbolento | u+ in funzione di y+ confrontato con la DNS, per due trattamenti di parete |
| 13 e 14 | Flusso transitorio; cilindro DFG | Numero di Strouhal e resistenza massima confrontati con il riferimento |
| 15 e 16 | Un piccolo progetto tuo | Una breve relazione con convergenza, bilanci, studio di mesh e un confronto con dati |
L’ultimo risultato è quello che ti fa assumere. Una relazione che dice «la perdita di carico è 1.245 Pa con un’incertezza numerica dell’1,6%, e concorda con la correlazione entro il 4%» vale più di qualsiasi numero di mappe a colori.
Errori comuni, come appaiono sullo schermo
- Una soluzione liscia mai confrontata con nulla. L’ha prodotta l’upwind. Può essere sbagliata di un fattore, come mostra l’esempio.
- Oscillazioni «curate» passando al primo ordine. Le oscillazioni dicevano che la mesh era troppo grossolana per lo schema. Il primo ordine nasconde il messaggio e tiene l’errore.
- Un tutorial ripetuto clic per clic. Gira, e non sai dire perché ogni impostazione sia quella che è. Cambia un’impostazione alla volta e prevedi l’effetto prima di lanciare.
- Un fattore d’attrito «sbagliato» del 30% in un tubo corto. È la lunghezza di ingresso, e il solutore ha ragione.
- Un caso turbolento senza controllo di y+. Il trattamento di parete non è mai stato nel suo campo di validità, e il modello non faceva ciò che pensi.
- I residui come unica prova di convergenza. Monitora la grandezza che riporterai e controlla il bilancio di massa; l’articolo sulla convergenza di questo blog spiega perché.
Che cosa farei
Se hai una laurea in ingegneria e qualche ora a settimana, segui l’ordine sopra: matematica che spiega l’errore, un piccolo solutore, casi di riferimento, poi il software. Non saltare l’esperimento del Péclet. Richiede una serata e cambia il modo in cui leggerai ogni mappa a colori per il resto della carriera.
Se usi già un software commerciale al lavoro senza aver avuto una formazione numerica, fai lo stesso percorso in forma compressa: l’esperimento del Péclet, la cavità contro Ghia, uno studio su tre mesh di un caso che conosci. Tre fine settimana. Comincerai a vedere quali dei tuoi vecchi risultati firmeresti ancora.
Se vuoi un aiuto strutturato, il nostro corso Modellazione CFD di flussi di liquidi e gas segue lo stesso ordine: impostazione, fisica, simulazione, verifica e interpretazione, e il corso Scambio termico nelle applicazioni ingegneristiche copre il gradino 7 e oltre. In ogni caso tieni la regola: nessun risultato lascia la tua scrivania finché non è stato confrontato con qualcosa che non hai calcolato tu.
Fonti e approfondimenti
- Barba and Forsyth: CFD Python, 12 steps to Navier-Stokes
- FeatFlow, TU Dortmund: DFG flow around cylinder benchmark 2D-2 (Re = 100)
- Oden Institute, University of Texas at Austin: turbulent channel DNS data
- Ansys Student: free student software and its limits
- Ghia, Ghia and Shin (1982), Journal of Computational Physics: high-Re solutions for incompressible flow
Servizio correlato
Revisione dei modelli e studi su misura
Capire quali conclusioni sono supportate dal modello, dove resta incertezza e cosa richiede verifica.
Scopri i servizi