FLUIDODINAMICA COMPUTAZIONALE - METODI NUMERICI PER LE EQUAZIONI DI NAVIER-STOKES

Ingegneria Meccanica e Gestionale

· METODI NUMERICI PER LE EQUAZIONI DI NAVIER-STOKES

· Valerio Marra

· INTRODUZIONE

· EQUAZIONI DI NAVIER-STOKES

· METODO DI PROIEZIONE DI CHORIN

· CODIZIONI AL CONTORNO

· DISCRETIZZAZIONE DEL PROBLEMA ED EQUAZIONE DI POISSON

· TRATTAMENTO DEL TERMINE CONVETTIVO

· TEST: SHEAR LAYER E INSTABILITà DI KELVIN-HELMHOLTZ

· ANALISI DELLE PERFORMANCE DEI METODI PRESENTATI

INTRODUZIONE

La tecnica numerica usata per risolvere il sistema di equazioni di Navier-Stokes per un flusso non stazionario e incomprimibile è basata sul metodo di proiezione sviluppato da Chorin. L'equazione di conservazione della quantità di moto viene risolta avanzandola temporalmente con i soli termini convettivo, diffusivo ed eventuali forzanti esterni per ottenere un campo di velocità “provvisorio”. Successivamente questo campo viene proiettato sullo spazio dei campi vettoriali solenoidali risolvendo un'equazione di Poisson per la pressione.

EQUAZIONI DI NAVIER-STOKES

Il moto dei fluidi Newtoniani, incomprimibili e isotermi è governato dal sistema di equazioni differenziali alle derivate parziali di Navier-Stokes, che rappresentano le equazioni fondamentali di continuità e di conservazione della quantità di moto:

(1)

(2)

dove t è il tempo, r la densità, u=(u, v) il vettore velocità del fluido, p la pressione idrostatica, m la viscosità cinematica e f_e rappresenta il vettore delle (eventuali) forze volumetriche esterne (e.g. la forza gravitazionale modellata come f_e=g). Si noti che in questa formulazione r e m sono ritenute costanti e che le incognite sono costituite dai campi di velocità (vettoriale) e di pressione (scalare).

Il metodo di soluzione del sistema di equazioni (1) e (2) proposto e descritto nel seguito è basato su un metodo di proiezione esplicito.

METODO DI PROIEZIONE DI CHORIN

La tecnica utilizzata per risolvere il sistema di equazioni di Navier-Stokes (1) e (2) è derivata dal metodo di proiezione di Chorin che è basato sul seguente teorema di decomposizione (di cui sarà omessa la dimostrazione).

Teorema di decomposizione

Sia D una regione del piano con contorno continuo ¶D. Un qualsiasi campo vettoriale u definito su D può essere unicamente decomposto nella forma:

(3a)

dove u è un campo vettoriale solenoidale con u×n=0 sul contorno ¶D e p è una funzione scalare (Figura 1).

Figura 1

Il teorema precedente permette di definire un operatore P, detto di proiezione ortogonale, che applicato a u* lo mappa sulla sua componente solenoidale u. L’operatore P risulta essere lineare e con le seguenti proprietà:

(3b)

(3c)

(3d)

Il metodo di proiezione consiste nell’applicare all’equazione (2), in assenza di forzanti esterni, l’operatore P:

(4)

Se u è solenoidale e ha componente normale nulla sul contorno, per continuità temporale lo stesso vale per ¶u/¶t e quindi P[¶u/¶t]=¶u/¶t. Quindi si ottiene:

(5)

Sebbene nÑ²u sia solenoidale, non necessariamente la sua componente normale è nulla sul contorno e quindi non vale P[nÑ²u]=nÑ²u. La forma (5) dell’equazione di conservazione della quantità di moto elimina la pressione ed esprime ¶u/¶t solo in termini della u. La pressione può essere quindi vista come la “parte gradiente” del termine -(u×Ñ)u+nÑ²u. La formulazione (5) oltre che d’interesse teorico lo è anche pratico per la definizione dell’algoritmo numerico che sarà presentato nel seguito.

CODIZIONI AL CONTORNO

Le precedenti equazioni devono essere corredate da opportune condizioni al contorno. L’equazione di Eulero (che si ottiene dalla (2) in assenza di forzanti esterni e viscosità nulla) valida per i fluidi ideali

(6)

è del primo ordine nello spazio e richiede una sola condizione al contorno. Su eventuali pareti solide e ferme si impone u×n=0 (free slip): il fluido non può penetrare la parete ma è libero di scorrervi tangenzialmente. Per le equazioni di Navier-Stokes il termine aggiuntivo nÑ²u aumenta l’ordine delle derivate spaziali da uno a due. Sia per ragioni sperimentali che matematiche questa situazione è accompagnata dall’incremento del numero di condizioni al contorno. Su una parete solida e ferma dobbiamo annullare anche la componente tangenziale della velocità e quindi u=0 (no-slip). Si fa notare che per i fluidi reali (m ¹ 0) il termine viscoso è responsabile della conversione continua di energia meccanica in energia interna (calore).

DISCRETIZZAZIONE DEL PROBLEMA ED EQUAZIONE DI POISSON

L’equazione di conservazione della quantità di moto (2) è risolta utilizzando il metodo delle differenze finite (FDM – Finite Difference Method) su un grigliato di tipo MAC (staggered MArker and Cell) con uno schema di discretizzazione temporale di tipo split esplicito. L’incomprimibilità del flusso fluido – i.e. il soddisfacimento della condizione di solenoidalità del campo di velocità (1) – è assicurata da un’equazione di Poisson per la pressione risolta utilizzando una tecnica multigriglia (MGM – Multi-Grid Method).

La pressione e le componenti della velocità sono discretizzate su una griglia cartesiana uniforme (Dx=Dy=h), la distribuzione delle variabili sulle celle è come detto di tipo MAC con le componenti della velocità definite sui bordi delle celle (Figura 2). I volumi di controllo W_u e W_v delle componenti della velocità sono mostrati in Figura 3 e sono utilizzati per la risoluzione del sistema di equazioni (1) e (2).

Figura 2 Figura 3

La tecnica di risoluzione può essere schematizzata nei tre passi seguenti:

1) Si costruisce una soluzione “provvisoria” u* dell’equazione (2) al tempo t_n+1=t_n+Dt (dove Dt è il passo temporale scelto) utilizzando il campo di velocità noto all’istante t_n:

(7)

L’operatore differenziale Ñ^h è discretizzato mediante differenze finite centrate mentre la discretizzazione del termine convettivo -Ñ^h×Õ^h è descritta in dettaglio nel seguito. L’equazione (7) è avanzata temporalmente usando una tecnica di tipo split esplicita: il campo di velocità è prima aggiornato con il tensore viscoso mÑ^huⁿ e successivamente con il termine convettivo.

2) Risolviamo per p l’equazione di Poisson:

(8)

La soluzione di questa equazione è ardua poiché il campo di velocità può esibire delle discontinuità. Tecniche di risoluzione del tipo multigriglia si sono rivelate particolarmente adatte alla soluzione di questo tipo di problema.

3) Proiettiamo il campo “provvisorio” u* sul campo uⁿ⁺¹ solenoidale:

(9)

dove p è soluzione della precedente equazione di Poisson (8) e garantisce che:

(10)

TRATTAMENTO DEL TERMINE CONVETTIVO

Il temine convettivo:

(11)

è discretizzabile utilizzando diverse tecniche, ognuna delle quali tesa a catturare diverse peculiarità fisico-matematiche del campo di velocità, per un fluido ideale o meno, soluzione di una determinata configurazione iniziale e al contorno.

Metodo delle differenze finite centrate

La componente del termine convettivo in direzione x è discretizzata secondo il seguente schema (Figura 4):

(12)

Analogamente si procede per la componente in direzione y.

Figura 4

Metodo upwind del primo ordine

Per la componente del termine convettivo in direzione x si deriva il seguente schema di discretizzazione:

(13)

dove:

(14)

Analogamente per la componente in direzione y. L’idea alla base del metodo upwind è la modifica della direzione del trasferimento di informazione a seconda del segno della velocità d’onda locale u(v) in direzione x(y).

Metodo di Godunov

Le tecniche di discretizzazione del termine convettivo fin qui proposte si rivelano adeguate alla soluzione di problemi con flussi continui. Tuttavia nella simulazione di flussi discontinui, comuni in applicazioni ingegneristiche di tipo avanzato, è necessario utilizzare opportuni strumenti matematici.

Nella trattazione di tale tipo di flussi, affinché la soluzione debole del problema sia unica, si deve introdurre la “condizione entropica”, cosí che al tendere a zero della viscosità la soluzione “viscosa” delle equazioni del moto (equazioni di Navier-Stokes) tenda alla soluzione non viscosa (equazione di Eulero) per le stesse condizioni iniziali. Nei flussi comprimibili e viscosi caratterizzati da un elevato numero di Reynolds possono esistere delle regioni limitate con forti gradienti delle variabili fisiche. Se il fluido è ideale (i.e. non viscoso) in queste regioni ci saranno delle onde di shock, vale a dire delle vere e proprie discontinuità fisico-matematiche in movimento. Quando la viscosità è piccola ma non nulla, la regione caratterizzata da gradienti elevati può in generale essere piú piccola del passo di griglia e quindi la soluzione numerica viscosa si comporta come una soluzione non viscosa. Di conseguenza, lo schema numerico che si intende utilizzare per risolvere le equazioni del moto per piccoli valori della viscosità deve essere capace di produrre soluzioni contenenti discontinuità. Per questa ragione, i metodi generalmente usati sono una diretta estensione dei metodi previsti per la soluzione dell'equazione non viscosa.

Nella ricerca della soluzione delle equazioni di Navier-Stokes secondo condizioni iniziali e al contorno che possono portare alla formazione di flussi caratterizzati da elevati gradienti, la discretizzazione del campo di velocità u porta alla situazione discussa sopra (anche per flussi con m >> 1) con la conseguente trattazione da parte dello schema numerico adottato di valori delle variabili non continui sul contorno dei volumi di controllo adiacenti.

Il metodo di Godunov, specializzato alla discretizzazione del termine convettivo non lineare dell'equazione di conservazione della quantità di moto, è una tecnica ai volumi finiti sviluppata nell'ambito dei fluidi comprimibili che si è mostrata particolarmente utile nel calcolo di flussi ideali estremamente complicati e con elevati gradienti.

L’idea che ne sta’ alla base è l’uso di soluzioni localizzate di problemi fisici monodimensionali per stimare il comportamento di flussi multidimensionali. Questa tecnica fornisce un trattamento robusto del termine convettivo non lineare anche a numeri di Reynolds elevati, può essere del secondo ordine spazialmente per flussi continui ed è stabile per condizioni iniziali discontinue (anche nel limite di viscosità nulla). Come nel metodo upwind, anche il metodo di Godunov modifica la direzione del trasferimento di informazione a seconda del segno della velocità d’onda locale, ma in più dissipa gli shock di rarefazione (risultato dei sopra discussi fenomeni) aggiungendo un termine viscoso artificiale.

Anche se la trattazione del metodo di Godunov esula dallo scopo di questo articolo poiché richiederebbe una discussione approfondita delle proprietà delle leggi di conservazione per i sistemi iperbolici, se ne è dato cenno al fine fornire una visione completa delle difficoltà incontrate e delle sfide che deve affrontare la fluidodinamica computazionale (CFD – Computational Fluid Dynamics).

Nel seguito verranno presentati i risultati di una simulazione numerica effettuata utilizzando nell’ordine: il metodo alle differenze finite centrate, il metodo upwind e quello di Godunov del secondo ordine (schema di tipo “predictor/corrector”).

TEST: SHEAR LAYER E INSTABILITà DI KELVIN-HELMHOLTZ

In questo test si risolve l’equazione di Eulero con il campo di velocità iniziale così definito:

(15)

Questo profilo corrisponde a una struttura costituita da due getti in controcorrente che si ripetono in direzione y con velocità lungo x, dove il parametro r controlla l'estensione della zona di transizione tra i due getti. La componente del moto lungo y è una piccola perturbazione sinusoidale di ampiezza d che innesca l'instabilità di Kelvin-Helmholtz. Il dominio d'integrazione è il quadrato unitario con condizioni contorno periodiche sia in orizzontale che verticale. In Figura 5 sono riportate le isolinee del capo di vorticità per shear layer su grigliato 512x512 al tempo t=0.4 s.

Figura 5

La zona di transizione (shear layer) tra i due getti evolve in una struttura periodica di grandi vortici che manifestano nel tempo delle strutture sempre più fini che si avvolgono a spirale (Figura 6). La simmetria delle condizioni iniziali da’ luogo a una soluzione che deve mantenere una simmetria rilevabile dalla soluzione numerica e che costituisce un buon test per la qualità dell'algoritmo numerico assunto per il termine convettivo. La figura 6 riporta le isolinee del capo di vorticità per shear layer su grigliato 512x512 al tempo t=1.8 s.

Figura 6

Nelle figure riportate di seguito sono mostrate, attraverso le isolinee del campo di vorticità (a) e il campo vettoriale di velocità (b) all’istante t=0.99 s, i risultati ottenuti su di un grigliato con 256² punti con il Metodo delle Differenze Finite Centrate (Figure 7), con il Metodo upwind del primo ordine (Figure 8) e con il Metodo di Godunov (Figure 9).

Figura 7a Figura 7b

Figura 8a Figura 8b

Figura 9a Figura 9b

ANALISI DELLE PERFORMANCE DEI METODI PRESENTATI

Come noto, il metodo alle differenze finite centrate non è adatto alla discretizzazione dell'equazione di Eulero essendo leggermente instabile poiché deficitario di un termine viscoso, per l’appunto, stabilizzante. Questo fatto è chiaramente rilevabile dalle figure che mostrano delle oscillazioni selvagge del campo di vorticità e un campo di velocità senza un preciso significato fisico in gran parte del dominio (Figure 7).

Si evince dall'analisi delle Figure 8 delle Figure 9 che nelle zone del dominio nelle quali il campo di velocità tende a zero il metodo upwind non riesce a diffondere correttamente i non realistici shock di rarefazione dovuti alle discontinuità nel campo di velocità. Questo fatto da luogo alla creazione di vortici spuri che non permettono una corretta rappresentazione della dinamica. Come anticipato, il metodo di Godunov a fronte di un costo computazionale solo leggermente più elevato, ma di uno sforzo teorico ben più rilevante, permetta una descrizione delle dinamica dell’evoluzione del campo di velocità maggiormente dettagliato e con un significato fisico maggiormente aderente alla realtà.