Appunti di Fluidodinamica Numerica

Introduzione alla fluidodinamica numerica per applicazioni industriali

Francesco Papi

22 giugno 2018

Appunti del corso di Fluidodinamica Numerica per Applicazioni Industriali tenuto dal prof. Roberto Pacciani e dal prof. Antonio Andreini, Università degli Studi di Firenze, anno accademico 2017/2018.

Indice

• I Modelli Matematici per la CFD 5
• 1 Introduzione 5
• 2 Equazioni di Governo 5
• 3 Proprietà matematiche di equazioni alle derivate parziali 10
• 4 Teoria delle Caratteristiche 11
• II Metodi di soluzione alle Differenze Finite 14
• 1 Introduzione 14
• 2 Metodi in avanti, centrali, all’indietro 14
• 3 Proprietà di uno schema Numerico 16
• 4 Analisi di stabilità di uno schema 17
• 5 Condizione CFL e Stabilità 19
• 6 Metodi Upwind e Metodi Centrati con dissipazione 21
• III Griglie Computazionali 22
• 1 Griglie Strutturate 22
• 2 Griglie non Strutturate 24
• 3 Staggered Grid 26
• IV Metodi di soluzione ai Volumi Finiti 28
• 1 Introduzione 28
• 2 Metodi per il calcolo di Volumi, Aree e Flussi 29
• 3 Schemi centrati con dissipazione artificiale 31
• 4 Metodi Upwind 33
• 5 Roe’s Upwind 35
• 6 Schemi Upwind Multidirezionali 36
• V Tecniche Numeriche di Calcolo 40
• 1 Metodi di Runge-Kutta 40
• 2 Tecniche numeriche per accelerare la Convergenza 42
• 3 Metodi di pressure-correction 45
• 4 Solutori Segregati e solutori Coupled 49
• 5 Trattamento numerico dei termini sorgente 50
• 6 Cenni a metodi di soluzione non-stazionari 50
• VI Modelli di Turbolenza 52
• 1 Fenomenologia 52
• 2 Energy Cascade 52
• 3 Equazioni di NS mediate alla Reynolds 53
• 4 Modelli di Turbolenza 54
• 5 Modelli Eddy Viscosity 54
• 6 Stress Transport Models 61
• 7 Metodi di Scale Resolving per la Turbolenza 61
• 8 Simulazioni LES 63
• 9 Sub Grid Scale Models 64
• 10 Griglie di calcolo LES 64
• 11 Altri tipi di simulazioni: DES, SST 66
• VII Condizioni al Bordo per equazioni di NS 67
• 1 Condizioni Fisiche 67
• 2 Imposizione numerica delle condizioni al contorno 68
• 3 Condizioni al bordo per simulazioni LES 70
• VIII Soluzione numerica di flussi reattivi 72
• 1 Fisica del problema 72
• 2 Combustione 72
• 3 Soluzione numerica di flussi reattivi turbolenti 74
• IX Soluzione numerica di flussi multi-fase 79
• 1 Introduzione 79
• 2 Atomizzazione di un liquido in un areiforme 79
• 3 Strumenti Numerici 81

Parte I: Modelli matematici per la CFD

1 Introduzione

Obiettivo della Fluidodinamica Numerica è la soluzione di problemi fluidodinamici tramite la discretizzazione e successiva soluzione numerica delle equazioni alle derivate parziali che governano il fenomeno. Nel presente lavoro ci si concentrerà sulla soluzione numerica di:

• Equazioni di Eulero
• Equazioni di Navier Stokes

Saranno dunque trascurati modelli di flusso a potenziale o flusso irrotazionale in quanto basati su pesanti approssimazioni e divenuti marginali con l’aumentare della potenza di calcolo a disposizione del progettista. Importante specificare che il punto di vista utilizzato è tipicamente quello Euleriano. Infatti, qualora il dominio in cui si muove il flusso sia noto, circostanza peraltro molto comune nelle applicazioni Industriali, questo modo di procedere risulta assai più semplice di un approccio Lagrangiano.

2 Equazioni di governo

Il modello fisico che permette di trattare un problema fluidodinamico in maniera più generale possibile è costituito dalle equazioni Navier Stokes. Da un punto di vista Euleriano vi sono essenzialmente due modi di scrivere le Equazioni di Governo:

• Formulazione Conservativa: Risulta nella scrittura di equazioni di conservazione di massa, momento ed energia
• Formulazione Non Conservativa: Risulta nella scrittura di Equazioni di Trasporto per le variabili sopra citate

Si noti che le due formulazioni sono equivalenti da un punto di vista Matematico, non lo sono invece da un punto di vista Numerico.

ρ, ρv, ρE Variabili Conservative

Si definiscano i ρ, p, v Variabili Primitive

Si definiscano i

Equazioni di Navier Stokes in formulazione conservativa

Equazioni di Governo in forma Integrale:

Z Z ∂ · ρ dV + ρv n dA = 0 (1)
∂t V A

Z Z Z Z Z ∂ · −ρv dV + n) dA = dV dA + dA ρv(v ρg pn τ n (2)
∂t V A V A A

Z Z Z Z Z ∂ · · ·ρE dV + ρH(v n) dA = v dV + n) dA + dA ρg τ (v qn (3)
∂t V A V A A

In cui E ed H sono rispettivamente l’Energia e l’entalpia specifiche, che per un gas perfetto:

2 2p v v1 + = c T + E = (4)
v−γ 1 ρ 2 2

2 2γ p v vH = + = c T + (5)
p−γ 1 ρ 2 2

Dove l’ultimo passaggio può esser facilmente svolto ricordando la Relazione di Mayer e l’equazione dei gas perfetti.

Si può facilmente passare da una formulazione integrale ad una formulazione differenziale ricordando il teorema della divergenza e trasformando dunque gli integrali di superficie in integrali di volume. Le equazioni definite sopra diventano allora:

Z Z Z ∂p ∂ ∇(ρv ∇(ρv· ·ρ dV + +n) dV = n) dV (6)
∂t ∂tV V V

Z Z Z Z Z ∂ ∇ · ⊗ − ∇p ∇τρv dV + (ρv v) dV = ρg dV dV + dV (7)
∂t V V V V V

Z Z Z Z Z ∂ ∇(ρvH) · ∇(τ ∇qρE dV + dV = ρg v dV + v) dV + dV (8)
∂t V V V V V

Imponendo l’arbitarietà del volume di integrazione si passa infine ad una formulazione differenziale:

∂ρ ∇(ρv ·+ n) (9)
∂t

∂ ∇ · ⊗ − ∇p ∇τρv + (ρv v) = ρg + (10)
∂t

∂ ∇(ρvH) · ∇(τ ∇qρE + = ρg v + v) + (11)
∂t

Gli argomenti delle divergenze sono dei Flussi. Si ricorda che l’operatore ∇ · (. . . ) abbassa l’ordine dell’entità in argomento. Si passa dunque ad esempio da un tensore 3x3 ad un vettore o da un vettore ad uno scalare.

Equazioni di Navier Stokes in formulazione non conservativa

Interessante vedere adesso come si può passare da una formulazione conservativa ad una non conservativa. Per fare questo si usa l’equazione di continuità esplicitandone le derivate:

∂ρv ∂ρ ∂ρ ∂v∂ρ i i+ = + v + ρ =0 (12)
i∂t ∂x ∂t ∂x ∂xi i i

ρ: Svolgendo la derivata totale di

dρ ∂ρ ∂ρ ∂x ∂ρ ∂ρi= + = + v (13)
idt ∂t ∂x ∂t ∂t ∂xi i

Da cui sostituendo in (12) si ottiene:

dρ ∂vi+ ρ =0 (14)
dt ∂xi

Allo stesso modo svolgendo le derivate per l’euquazione (10), con riferimento alla sola parte sinistra dell’equazione in componenti:

∂v∂ρv ∂ρv v ∂p ∂v ∂vii i j i i+ρ+ = v [ + ρ ] + ρv (15)
i i∂t ∂x ∂t ∂x ∂t ∂xj i j
{z }| =0

Si ricordi inoltre che:

∂v dv∂v ii + ρv = ρρ (16)
i∂t ∂x dtj

Operando in maniera analoga anche per l’equazione dell’energia si arriva a scrivere il sistema delle equazioni di governo in forma non conservativa come:

dρ ∂v = 0+ ρ idt ∂x i ∂τ∂pdv ijρ = ρg ++ ρi idt ∂x ∂xi j∂τ∂E 1∂H ij+ u = u + g u i i i i∂t ∂x ρ ∂xi j

Le formulazioni conservativa e non conservativa non sono equivalenti a livello numerico. Sopratutto nel caso comprimibile la rappresentazione non conservativa non è adatta a tenere traccia di tutta la fenomenologia, in particolare le discontinuità d’urto.

τ Deve esser specificata una formulazione per gli sforzi di taglio. L’assunzione più comune e soddisfacente per la maggior parte dei problemi di trattazione comune è quella di considerare il fluido Newtoniano. Il tensore degli sforzi di taglio assume dunque la forma:

∂v2 ∂vi j− + µ(τ = µ∇vδ + ) (17)
ij ij3 ∂x ∂xj j

Gli sforzi di taglio sono dunque proporzionali al gradiente di velocità. La costante di proporzionalità prende il nome di viscosità dinamica. Nel caso di fluido incomprimibile sparisce la parte contenente l’operatore ∇.

Equazioni di Navier Stokes in forma vettoriale

Le equazioni di governo possono essere espresse anche in forma maggiormente sintetica con una notazione vettoriale. Sia il vettore delle variabili dipendenti, conservative, raggruppate in un vettore:

TU = ρ[1, u, v, w, E] (18)

Siano definiti inoltre i Vettori:

2 TF = ρ[u, u + p, uv, uw, uH] (19)

2 TG = ρ[v, vu, v + p, vw, vH] (20)

2 TH = ρ[w, wu, wv, w + p, wH] (21)

TS = [0, g , g , g , g u + g v + g w] (22)
x y z x y z

F = [0, τ , τ , τ , τ u, τ v, τ w] (23)
v xx yx zx xx xy zx

G = [0, τ , τ , τ , τ u, τ v, τ w] (24)
v xy yy zy xy yy zy

H = [0, τ , τ , τ , τ u, τ v, τ w] (25)
v xx yz zz xz yz zz

Allora le Equazioni di Navier Stokes posso essere scritte in forma sintetica per un sistema di riferimento Cartesiano:

∂U ∂F ∂G ∂H ∂F ∂G ∂Hv v v+ + + = S + + + (26)
∂t ∂x ∂y ∂z ∂x ∂y ∂z

Equazioni di Navier Stokes in un sdr curvilineo

Per ottenere una scrittura più generale delle equazioni è necessario passare da un sistema di riferimento cartesiano, trattato fino a questo momento, ad un sistema di riferimento curvilineo. Tale passaggio è necessario sopratutto nel caso in cui si usi una metodologia di discretizzazione alle differenze finite, in quanto questo facilita il processo di discretizzazione. In un metodo risolutivo ai volumi finiti questo passaggio non è strettamente necessario, tuttavia molti moderni solutori fanno uso di quantità definite in sistemi di riferimento locali, per cui questo risulta conveniente anche in questi casi. Un esempio in cui generalizzare il sistema di riferimento risulta utile è lo studio di un condotto rotorico di turbomacchina in condizioni stazionarie, dove l’intero volume di controllo ruota attorno all’asse della machina. Sia un generico sistema di riferimento in coordinate curvilinee composto dalle seguenti componenti:

 = (x, y, z, t)
 η = η(x, y, z, t)
ζ = ζ(x, y, z, t)
 τ = t

Si noti che le coordinate curvilinee dipendono dal tempo, il sistema di riferimento infatti si muove anche se il fenomeno è stazionario. = ∂/∂x, = x y ∂/∂y, . . . J calcolano la lunghezza tra gli spazi. Lo Jacobiano correla il volume tra gli spazi. Si ottiene le equazioni di governo in coordinate generalizzate per sistema di riferimento fisso:

un −1 −1 −1 −1U F G H∂J ∂J ∂J ∂J −1 S++ + + = J∂t ∂ ∂η ∂ζ −1 −1 −1F G H∂J ∂J ∂Jv v v+ + + (27)
∂ ∂η ∂ζ

U, F, G, H I termini sono ottenuti:

∂ ∂ ∂U F= U = (F + G + H ) (28)
∂x ∂y ∂z∂η ∂η ∂η ∂ζ ∂ζ ∂ζ

G H= (F + G + H ) = (F + G + H ) (29)
∂x ∂y ∂z ∂x ∂y ∂z

Supponiamo di porci in un sistema di riferimento curvilineo rotante a velocità ω attorno all’asse x del sistema di riferimento cartesiano. Allora posso scrivere:

z(τ ) = r sin(ωt) y(τ ) = r cos(ωt) (30)

−ωz −ωyx = 0 y = z = (31)
τ τ τ

Allora riscrivendo i termini della velocità ricavati sopra e facendone la derivata parziale in si ottiene:

∂ ∂y ∂z−Ω= (32)
∂t ∂t ∂t

∂y ∂z= v = w: τ e che E scambiando le derivate in t e in e ricordando che

∂v −Ωw= (33)
∂τ z

Ragionamento analogo può esser fatto anche per τ. Allora l’introduzione del nuovo vettore tempo si riflette sul vettore delle incognite, la cui variazione nel tempo diviene:

∂U ∂U ∂U ∂U ∂U= + + η + ζ + I (34)
t t t∂t ∂τ ∂ ∂η ∂ζ

Con I vettore delle Forze di Inerzia:

T−ρωw, −ρωv,I = [0, 0, 0] (35)

Riassemblando tutti i termini si ottiene una formulazione formalmente quasi identica alla formulazione vista per il sistema di riferimento fisso. Compare in più un termine relativo alle forze di inerzia. Si ottiene una scrittura delle equazioni di Navier Stokes in coordinate curvilinee per un sistema di riferimento mobile:

−1 −1 −1 −1U F G H∂J ∂J ∂J ∂J −1 −1+ + + = J S + J I+
∂t ∂ ∂η ∂ζ −1 −1−1 G HF ∂J ∂J∂J v v v+ ++ (36)
∂ ∂η ∂ζ

Le variabili che non risentono del cambiamento di stema di riferimento sono ovviamente la pressione, gli sforzi di taglio e la densità. Il vettore delle incognite contiene le componenti cartesiane del vettore soluzione. Non è necessario allora fare un cambio di coordinate per leggere gli output. Una scrittura di questo tipo permette di trattare sia problemi steady che problemi unsteady in modo del tutto generale. Partire da una formulazione generale è la chiave di volta per costruire un sistema che permetta di risolvere varie casistiche. Inoltre partire da una formulazione di questo tipo permette di espandere le casistiche che il solutore è in grado di maneggiare in caso questo dovesse esser necessario.

Equazioni di Eulero

Le Equazioni di Eulero governano il moto di fluidi non viscosi e non conduttivi. Sono ricavate a partire dalle più generali equazioni di Navier Stokes annullando i termini viscosi. Spesso inoltre si trascurano anche eventuali termini sorgente. In questo modo le equazioni perdono i termini del secondo ordine e diventano equazioni convettive del primo ordine. Questa assunzione dal punto non semplifica la soluzione numerica del problema quanto ci si potrebbe aspettare; i termini difficili da trattare in modo consistente sono infatti i termini del primo ordine e non i termini viscosi. Applicando le semplificazioni sopra citate alla (36) si ottengono le equazioni di Eulero:

−1 −1 −1−1 F G HU ∂J ∂J ∂J∂J + + + =0 (37)
∂t ∂ ∂η ∂ζ

Comunemente le equazioni di Eulero vengono linearizzate. In questa forma è possibile studiarne le proprietà matematiche in modo da poter specificare adeguate condizioni al contorno e definire vettori di flusso numerici. Analizzando una componente della (37), si ottiene:

−1−1 FF ∂J ∂U ∂U∂J = = A(U ) (38)
∂ ∂U ∂ ∂

Con ragionamenti analoghi si ottiene:

−1−1 ∂J H∂J G C(U ) =B(U ) = (39)
∂U ∂U

A, B, C Le matrici sono chiamate matrici jacobiane non viscose.

• 6A, B, C = const. Se ottengo Equazioni di Eulero Quasi-Lineari
• A, B, C = const. Se ottengo Equazioni di Eulero Lineari

Allora raggruppando tutti i termini sopra discussi si ottengono le equazioni di Eulero Linearizzate:

−1∂J U ∂U ∂U ∂U+ A + B + C =0 (40)
∂t ∂ ∂η ∂ζ

3 Proprietà matematiche di equazioni alle derivate parziali

La maggior parte dei problemi fisici di nostro interesse può esser classificata come segue:

• Problemi di Equilibrio: sono descritti da equazioni alle derivate parziali ellittiche. Sono tipicamente problemi di tipo stazionario la cui soluzione dipende dalle condizioni al contorno sul dominio. Esempi di questo tipo di problemo sono Ricavare la distribuzione di temperatura in un solido oppure problemi di flusso a potenziale.
• Problemi di Propagazione: sono descritti da equazioni alle derivate parziali Iperboliche o Paraboliche. Sono tipicamente problemi di tipo evolutivo. Se sono assegnate esclusivamente condizioni iniziali per risolvere il problema si parla di problemi Iperbolici, se sono assegnate sia condizioni al bordo che iniziali si parla di problemi parabolici.

Si consideri una generica equazione alle derivate parziali del secondo ordine:

22 2∂ u∂ u ∂ u+ b + terminidiordineinf eriore = sa + c (41)
2 2∂x ∂x∂y ∂y

2 −4ac∆ = b

Con riferimento al discriminante, è possibile riepilogare quanto detto nella tabella seguente:

4 Teoria delle caratteristiche

L’equazione di Eulero Linearizzata (40) può essere riformulata come un sistema di propagazione d’onda. La condizione di Iperbolicità del sistema (40) permette di scrivere la soluzione del sistema nella foma:

k·s−wtU = Û e (42)

La relazione di cui sopra esprime dunq