Appunti di Fluidodinamica Numerica

Introduzione alla Fluidodinamica Numerica per

Applicazioni Industriali

Francesco Papi

22 giugno 2018

Appunti del corso di Fludidodinamica Numerica per Applicazioni Industriali tenuto

dal prof. Roberto Pacciani e dal prof. Antonio Andreini, Università degli Studi di

Firenze, anno accademico 2017/2018

1

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

2

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

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

4

Parte I

Modelli Matematici per la CFD

1 Introduzione

Obiettivo della Fluidodinamica Numerica è la soluzione di problemi fluidodi-

namici tramite la discretizzazione e successiva soluzione numerica delle equa-

zioni 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

domino 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

Goveno:

• Formulazione Conservativa: Risulta nella scrittura di equazioni di conser-

vazione 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 Mate-

matico, non lo sono invece da un punto di vista Numerico.

ρ, ρv , ρE Variabili Conservative

Si definiscano i

ρ, p, v Variabili Primitive

Si definiscano i

Euqazioni di Navier Stokes in Fomulazione 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

5

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 2

p v v

1 + = c T +

E = (4)

v

−

γ 1 ρ 2 2

2 2

γ p v v

H = + = 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 formulazio-

ne 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 ∂t

V 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 Interes-

sante 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 ∂x

i i i

6

ρ:

Svolgendo la derivata totale di

dρ ∂ρ ∂ρ ∂x ∂ρ ∂ρ

i

= + = + v (13)

i

dt ∂t ∂x ∂t ∂t ∂x

i i

da cui sostituendo in (12) si ottiene:

dρ ∂v

i

+ ρ =0 (14)

dt ∂x

i

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 ∂v

i

i i j i i

+ρ

+ = v [ + ρ ] + ρv (15)

i i

∂t ∂x ∂t ∂x ∂t ∂x

j i j

{z }

| =0

Si ricordi inoltre che: ∂v dv

∂v i

i + ρv = ρ

ρ (16)

i

∂t ∂x dt

j

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

+ ρ i

dt ∂x

 i

 ∂τ

∂p

dv ij

ρ = ρg +

+ ρ

i i

dt ∂x ∂x

i j

∂τ

∂E 1

∂H

 ij

+ u = u + g u

 i i i i

∂t ∂x ρ ∂x

i j

Le formulazioni conservativa e non conservativa non sono equivalenti a livello

numerico. Sopratutto nel caso comprimibile la rappresentazione non conserva-

tiva 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’assun-

ij

zione 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: ∂v

2 ∂v

i j

− + µ(

τ = µ∇vδ + ) (17)

ij ij

3 ∂x ∂x

j 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 notazio-

ne vettoriale. Sia il vettore delle variabili dipendenti (conservative) raggruppate

in un vettore: T

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

7

Siano definiti inoltre i Vettori: 2 T

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

2 T

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

2 T

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

T

S = [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 ∂H

v v v

+ + + = S + + + (26)

∂t ∂x ∂y ∂z ∂x ∂y ∂z

Equazioni di Navier Stokes in un sdr curvilineo Per ottenere una scrit-

tura più generale delle equazioni è necessario passare da un sistema di riferi-

mento 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 rife-

= ∂/∂x, =

rimento infatti si muove anche se il fenomeno è stazionario. 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 −1

U F G H

∂J ∂J ∂J ∂J −1 S+

+ + + = J

∂t ∂ ∂η ∂ζ −1 −1 −1

F G H

∂J ∂J ∂J

v v v

+ + + (27)

∂ ∂η ∂ζ

8

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 −ωy

x = 0 y = z = (31)

τ τ τ

Allora riscrivendo i termini della velocità ricavati sopra e facendone la deri-

t

vata 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 ∂t ∂t

∂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 equa-

sistema di riferimento

zioni di Navier Stokes in coordinate curvilinee per un

mobile:

−1 −1 −1 −1

U F G H

∂J ∂J ∂J ∂J −1 −1

+ + + = J S + J I+

∂t ∂ ∂η ∂ζ −1 −1

−1 G H

F ∂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.

9

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 an-

che 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 quan-

to 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 H

U ∂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 F

F ∂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.

• 6

A, 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 deri-

vate 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.

10

• 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:

2

2 2

∂ u

∂ u ∂ u

+ b + terminidiordineinf eriore = s

a + c (41)

2 2

∂x ∂x∂y ∂y 2 −4ac

∆ = b

Con riferimento al discriminante , è possibile riepilogare quanto

detto nella tabella seguente: ∆ PDE Type

< 0 Ellittica

=0 Parabolica

> 0 Iperbolica

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−wt

U = Û e (42)

La relazione di cui sopra esprime dunq