Che materia stai cercando?

Appunti di Fluidodinamica Numerica

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.
Il documento non è stato visto o approvato dai prof. Andreini o Pacciani e potrebbe pertanto contenere imprecisioni o incompletezze.

Esame di Fluidodinamica Numerica per Applicazioni Industriali docente Prof. R. Pacciani

Anteprima

ESTRATTO DOCUMENTO

nel calcolo dei termini metrici che servono per potersi ricondurre al sistema di

rifermento in coordinate locali.

Figura 6: TE grigia a C/H

Figura 7: TE griglia ad O

2 Griglie non Strutturate

Griglie di questo tipo sono più adatte ad un approccio moderno alla CFD ca-

ratterizzato spesso da domini computazionali molto estesi o complessi difficili

da meshare con griglie strutturate.

Cartesian Trimmed Grids Questo è un approccio intermedio tra una griglia

strutturata ed una non strutturata. Si utilizzano celle quadrate con aspect ratio

unitario. Si passa da una spaziatura più grande ad una più fine dividendo in

4 la cella di calcolo. Questo potrebbe creare problemi nei contorni tra una

dimensione della griglia e la successiva. I contorni vengono trattati quindi con

due metodologie:

• Approccio Poliedrico: corrisponde ad un approccio di tipo non-strutturato.

Prendendo come riferimento il quadrato rosso in figura, il solutore consi-

dera al posto di una cella quadrata un poliedro con 5 facce corrispondenti

alle rispettive facce degli altri elementi della griglia

24

• Overlapping-Grids: Si considera griglia fine e griglie più grandi supe-

rimposte una sull’altra. Il valore all’interfaccia tra le griglie è calcola-

to come interpolazione tra i valori calcolati nei livelli di griglia coinvolti

nell’interfaccia.

Un altra criticità di questa discretizzazione sono i bordi fisici della griglia che

molto spesso non sono ortogonali alle celle. Il processo tramite il quale la griglia

è adattata ad un generico contorno curvilineo è quindi detto body fitting.

Griglie non-strutturate Non vi sono restrizioni sul tipo di cella da utilizzare.

Questo permette una notevole libertà nella scelta della cella che nella pressochè

totalità dei casi non è di un solo tipo. Si genera quindi una griglia di tipo ibrido,

gli elementi più comuni sono riportati in figura 8.

Figura 8: Elementi griglia non strutturata

Un’ altra possibilità è quella di realizzare griglie poliedriche con elementi

come in figura 9. Figura 9: Elementi Poliedrici

25

Gli elementi vengono creati a partire da una mesh di tipo tetraedrico. Il van-

taggio principale di queste formulazioni è la riduzione di problemi di distorsione

delle celle, che risultano avere un aspect ratio generalmente migliore.

Mesh quality assestment Il parametro principale per valutare la qualità di

una mesh tetraedrica è la Skewness: −

Area Ideale Area Reale

Skewness = (88)

Area Ideale

Per area ideale di intende l’area di un triangolo iscritto in una circonferenza che

circoscrive l’elemento reale.

Un altro parametro molto utilizzato è l’Aspect Ratio, che dovrebbe essere

più vicino possibile ad 1. La skewness ha un effetto sulla discretizzazione delle

equazioni.

Figura 10: A sx griglia con bassa skewness. a dx griglia con alta skewness

Nel caso si utilizzi una differenza centrata per calcolare il termine di flusso

si ha: −

Z φ φ

A B

n (Γgradφ) dA = n (Γgradφ)∆A = n (Γ )∆A (89)

i i i ∆

A n

Con riferimento alla figura 10 se il segmento PA non è allineato al versore il

termine di flusso calcolato non corrisponde l flusso reale attraverso la superficie.

Il flusso calcolato ha infatti una componente parallela alla superficie che divide

le celle di calcolo, deve essere allora introdotto un termine di cross-diffusion:

− · − · −

φ φ n n∆A φ φ e e ∆A φ φ

A B i A P η i A P

n (Γ )∆A = (90)

i · ·

∆ n e ∆ n e ∆η

3 Staggered Grid

L’utilizzo di uno schema centrato per discretizzare la pressione potrebbe portare

ad errori nel calcolo dei gradienti. Questo risulta evidente se si studia un test

case come quello riportato nella pagina successiva:

26

Con una griglia di questo tipo si vanno a calcolare i gradienti:

p p

∂p N S

= =0 (91)

∂y 2δy

Il gradiente calcolato come si può vedere è nullo, in evidente disaccordo con

quanto riportato in figura sopra. Uno dei modi per risolvere questo problema è

introdurre una griglia di tipo staggered :

Figura 11: Staggered Grid

Si opera un shift dei nodi su cui vengono salvati i valori di pressione e di

alcune componenti di velocità. Con riferimento alla figura 11, vengono salvati

nei nodi i valori delle grandezze scalari quali densità, pressione e temperatura,

sulle facce delle celle vengono invece salvati i valori delle componenti in direzione

y e in direzione x della velocità, identificate dalle frecce. Questo disaccoppia-

mento risolve i problemi discussi in precedenza ed è un formalismo utilizzato

anche nella trattazione dei metodi di pressure correction nel seguito di questo

lavoro. Staggered grid è utilizzata su griglie di tipo strutturato. Applicazioni

moderne fanno uso di altri metodi quali interpolazione di Rie-Chow che evitano

la necessità di utilizzare una griglia staggered.

27

Parte IV

Metodi di soluzione ai Volumi

Finiti

1 Introduzione

In modo simile a quanto accade nel caso delle differenze finite l’approccio ai

volumi finiti presuppone di semplificare il domino di calcolo suddividendolo

con una griglia. Ogni cella della griglia determina un volumetto elementare

per il quale sono scritte e risolte le equazioni di governo. L’elemento di base

della discretizzazione non è più quindi il generico nodo della mesh come avviene

nelle differenze finite, bensì il volume della cella di calcolo stessa. Le relazioni

fisiche che descrivono il problema e da cui questa trattazione comincia sono

le equazioni di Navier-Stokes scritte in forma integrale, riportate in seguito in

forma vettoriale per completezza:

Z Z Z

∂ ·

U dV + F n dA = S dV (92)

∂t V A V

La relazione (92) viene scritta per ogni cella di calcolo. Il termine di flusso

F = F + F

comprende sia i termini diffusivi che quelli convettivi; . Il

dif f conv

metodo di discretizzazione ai volumi finiti prevede di considerare costanti tutte

le grandezze all’interno della singola cella. L’equazione (92) diviene allora nella

sua forma discetizzata: ∂U X ·

+ A (F

Ω n ) = ΩS (93)

k k k

∂t k

Gli integrali di volume sono diventati dei semplici prodotti del volume della

cella per i termini medi definiti su tutta la cella. Gli integrali di superficie che

esprimono i termini di flusso sono diventati delle sommatorie sulle facce della

cella dei flussi medi entranti o uscenti da ciascuna faccia.

Si parla di sistemi centrati se U è applicato all’intera cella. Si parta di

sistemi nodali se U è applicato al singolo nodo come nel caso delle differenze

finite. Non si ha nessun vantaggio ad usare uno schema nodale rispetto ad uno

schema centrato. Inoltre uno schema centrato consente di imporre in maniera

più semplice e intuitiva le condizioni al bordo. In uno schema centrato non si è

interessati alle coordinate della mesh che servono unicamente a calcolare volumi

ed aree. Come conseguenza si ha che una discretizzazione ai volumi finiti è meno

sensibile alla qualità della mesh.

Si vanno a numerare quindi le celle e non i nodi. Le facce tra una cella e

l’altra, dove si vanno a calcolare i termini di flusso, vengono numerate con indice

come in figura 12.

Altro pregio della discretizzazione ai volumi finiti è la conservatività dello

schema. Questa è infatti tutta a carico del termine di flusso che deve essere

uguale opposto per facce di celle adiacenti in modo da garantire la conservazione

delle varie grandezze. Nei volumi finiti i termini di flusso sono calcolati una sola

volta per ogni faccia di celle adiacenti e risultano quindi conservativi.

28

Figura 12: Numerazione facce e griglia di calcolo

Si introduce il concetto di Residuo:

∂U X ·

− n ) + ΩS

Ω = A (F (94)

k

k k

∂t k

| {z }

=Residuo

In un caso Stazionario risolvere il problema corrisponde ad annullare il Re-

siduo. Poichè si adotta un metodo di soluzione numerico il residuo non sarà

completamente annullato ma si riterrà il calcolo concluso una volta che esso

scende al di sotto di un certo valore di soglia.

2 Metodi per il calcolo di Volumi, Aree e Flussi

Calcolo di Aree e Volumi Si sviluppa l’eq. (93) in un caso 2D:

∂U

Ω + [F ∆y + G ∆x + F ∆y + . . . ] = ΩS (95)

AB AB AB AB BC BC

∂t

Se la griglia è cartesiana e dunque le celle hanno la forma di rettangoli, la

relazione si semplifica notevolmente:

∂U

Ω + [F ∆y + G∆x] = ΩS (96)

AB

∂t

Questo approccio molto vantaggioso è usato in problemi del tipo immer-

sed body. Si usa una griglia cartesiana su tutto il volume di calcolo, i volumi

adiacenti la superficie del corpo immerso vanno trattati separatamente.

Il solutore deve essere ovviamente in grado di calcolare aree e volumi nel

caso si tratti un problema 3D. Nel caso 2D allora l’area è calcolata come:

1 1

|(A − ∧ − − −

Ω= C) (B D)| = (∆x ∆y ∆xBD ∆y ) (97)

AC BD AC

2 2

Per il calcolo di aree e volumi in un caso 3D si faccia riferimento alla figura

2. L’area di una generica faccia può esser calcolata come:

1 |(A − ∧ −

A = C) (B D)| (98)

ABCD 2 29

Figura 13: Cella di calcolo VF

Per il calcolo del volume invece si scompone la cella in tetraedri. Con

riferimento all’eq. (93), dividendo per il volume entrambi i membri si ottiene:

∂U A k

X ·

+ F n = S (99)

k k

∂t Ω

k A

Si deve dunque calcolare il temine , questo è per ciascuna faccia pari ad

k

uno dei termini metrici, infatti: Ω = det J = J (100)

Introduzione al calcolo dei termini di flusso Il termine di flusso va valu-

tato con riferimento allo stato sulla faccia, che non è in generale noto in quanto

in un metodo ai volumi finiti è noto lo stato sulla cella. Si hanno in generale

due opzioni:

• = f (F , F ).

F Per applicare questo

Calcolo tramite interpolazione di 1 i i+1

i+ 2 i i + 1

metodo si dovrebbe creare due facce fittizie in e in modo da poter

calcolare i relativi flussi.

• U F = f (U ).

Calcolo tramite estrapolazione di , quindi Si è spo-

12 12 12

i+ i+ i+

stato il problema da F ad U. Questo metodo è più semplice e di maggiore

applicazione.

I risultati ottenuti per i due metodi sopra elencati non sono uguali in quanto

U è fare una

il problema non è lineare. Il modo più semplice per calcolare 1

i+ 2

U U

media algebrica tra e , questo è tuttavia nella pratica improponibile in

i i+1

quanto porta a schemi di calcolo instabili. La tipologia degli schemi numerici

che è possibile utilizzare ricalca quanto visto nel caso delle Differenze finite;

metodi upwind e schemi centrati con dissipazione artificiale.

30

3 Schemi centrati con dissipazione artificiale

Schemi centrali non considerano la direzione di propagazione delle perturbazioni.

In altre parole lo schema di calcolo rimane formalmente identico al variare della

direzione del flusso. Questo porta, come visto in precedenza, all’insorgere di

instabilità dovute al disaccoppiamento dei nodi pari e dispari che devono essere

smorzate. Questo viene fatto tramite la dissipazione artificiale, che sarà oggetto

di questa sezione.

Si consideri una equazione monodimensionale di convezione semplice discre-

tizzata ai volumi finiti con termine di dissipazione artificiale:

A

∂U k

X ·

+ (F n ) = S + D(U ) (101)

k k

∂t Ω

k

Un modello di questo tipo, seppure molto semplice, permette di fare consi-

derazioni interessanti.

Schema di Jameson Lo schema di Jameson è uno degli schemi più no-

ti. Il prof. Jameson ha proposto infatti un modo per scrivere la funzione di

D(U )

dissipazione artificiale vastamente utilizzato.

X X

2 4 2 2 2 4 4 4

D(U ) = D (U ) + D (U ) = [D + D + D + D + D + D ] (102)

k k η ζ η ζ

k k

Gli operatori del secondo ordine servono a catturare correttamente gli ur-

ti. Per una adeguata rappresentazione dell’urto infatti lo schema deve andare

localmente al primo ordine.

In generale si hanno due tecniche per evidenziare la presenza di un urto:

• Shock Fitting: L’utro non viene risolto ma introdotto dall’esterno. Si trat-

ta di usare le equazioni di governo con approccio alle caratteristiche; risol-

vendo da destra e da sinistra si riesce a definire la discontinuità dell’urto.

L’urto vero e proprio tuttavia non vinene definito. Questo metodo è molto

accurato in quanto evidenzia una vera e propria discontinuità in corrispon-

denza dell’urto pur essendo estremamente complesso matematicamente.

Le Università di Bari e di Torino hanno lavorato a lungo sullo sviluppo

di questi metodi riuscendo ad applicarli solo a profili bidimensionali. Le

applicazioni pratiche sono molto limitate.

• Shock Capturing: Gli urti sono risolti con metodi al primo ordine. La zona

di fortissimo gradiente viene spalmata su una certa estensione spaziale,

seppur molto limitata, senza introdurre una vera e propria discontinuità.

La rappresentazione risulta quindi essere meno accurata ma più affidabile

computazionalmente. Sono allo studio sistemi che permettano di contenere

la discontinuità all’interno di una singola cella di calcolo.

Il generici termini di secondo e quarto ordine nella relazione (102) possono

esser scritti come: (2) −

2 + )δ U

D = δ (Λ (103)

12

(i− )

1

i+ 2

(4) − −

4 + +

D = δ (Λ )δ δ δ U (104)

12 )

(i−

1

i+ 2

31

Λ

Il termine viene calcolato come media algebrica:

12

(i− ) 1

Λ (Λ + Λ )

= (105)

1 (i) (i+1)

(i− ) 2

2

Λ

Il termine viene scritto come:

(i)

Λ = φ r con r = Raggio Spettrale (106)

(i)

SI ricorda che il raggio spettrale è pari al valore assoluto del massimo

autovalore del problema e vale:

1 q i

h 2 2 2

r = U + a + + (107)

x y z

J

U = u + v + w (108)

x y z

Inoltre: σ σ

r r

η

φ = 1 + con 0 < σ < 1

+ (109)

r r

ζ

| {z }

=AR

φ φ

Si evidenzia nel termine l’aspect ratio della cella (AR). Il termine allora

va a togliere dissipazione dove AR è molto piccolo. Si va in questo modo a

togliere dissipazione nello strato limite, dove per rappresentare adeguatamente

i forti gradienti presenti si adottano celle molto schiacciate. Nello strato limite

la dissipazione fisica introdotta dai termini viscosi è sufficiente a stabilizzare il

metodo e non è necessario introdurre ulteriore dissipazione artificiale. Per contro

se si va a schiacciare la griglia in zone che non fanno parte dello strato limite

il metodo potrebbe diventare instabile a causa della insufficiente dissipazione

artificiale.

Gli altri termini che rimangono da determinare per chiudere il modello

valgono: 1

2 = max (ν , ν , ν , ν )

(110)

i−1,j,k i,j,k i+1,j,k i+2,j,k

1

i+ 2

2 1 1

4 2

= max (0, ) (111)

1

1

i+ i+

64 64

2 2

P + P 2P

i−1,j,k i+1,j,k i,j,k

| |

ν = (112)

i,j,k P + P + 2P

i−1,j,k i+1,j,k i,j,k

ν è un sensore per gli urti basato su un gradiente di pressione. Al numeratore

vi è una derivata al quadrato, a denominatore una media su tre celle.

2

• →

ν

Se il gradiente di pressione è dolce è piccolo. Allora è piccolo

2 2 4 4 4

− → 6 → ≈

k > 0 = 0 k Lo schema è dunque prevalentemente del

quarto ordine. 2 2 2 4

• → − → ≈ →

ν k < 0 0

Se si ha un urto è grande. Allora è grande

Lo schema è dunque prevalentemente del secondo ordine.

Lo schema è robusto e non compaiono assunzioni riguardo al tipo di fluido.

Questo oltre alla sua buona accuratezza hanno contribuito alla sua diffusione.

32

4 Metodi Upwind

Si vuole in questa sezione generalizzare ad un approccio a volumi finiti mul-

tidimensionale le considerazioni fatte in precedenza per i metodi Upwind. Si

era visto come questi metodi decentrano le derivate dei termini spaziali tenendo

conto della direzione del flusso. Si ottiene uno schema stabile del primo ordine.

La ragione di tale stabilità sta nel fatto che se lo schema rispetta la condizio-

ne CFL è possibile soddisfare le condizioni di compatibilità. A differenza degli

schemi centrati, schemi upwind considerano la direzione di propagazione del-

l’informazione, introducendo nelle equazioni discretizzate le proprietà fisiche di

propagazione delle varie grandezze.

Figura 14: Andamento delle variabili discrete sulla griglia di calcolo

Si consideri per il momento una situazione monodimensionale, utile a fare

alcune considerazioni. Si introducono con riferimento alla figura sopra gli stati

di riferimento left e right, i quali si riferiscono allo stato del fluido sull’interfaccia

destra e sinistra della generica superficie di calcolo. Poichè nei volumi finiti lo

stato fisico si riferisce all’intera cella, in prossimità della superficie che delimita

una cella dalla adiacente si hanno delle discontinuità. Se tali discontinuità fos-

sero fisiche si dovrebbe andare a risolvere un problema di Riemann, in questo

caso deve essere risolto un problema di Riemann numerico.

Flux Difference Splitting Se il problema è semplicemente monodimensio-

i i + 1 U

nale a primo ordine si prende la variabile U in o come a seconda

1

i+ 2

della direzione del flusso. In generale tuttavia il problema è più complesso. Per

ricostruire i valori delle variabili all’interfaccia si utilizza l’approccio alle carat-

teristiche. Questo è necessario affinchè lo schema sia consistente con le direzioni

di propagazione delle perturbazioni. La soluzione può esser dunque scomposta

in contributi d’onda semplice: X (k)

δU = δω r (113)

k

k

(k)

r δω

Il termine indica la direzione della perturbazione di ampiezza .

k

Seguendo lo stesso formalismo, affinchè lo schema sia conservativo deve valere:

− −

F = F + δF (114)

L

1

i+ 2

+ +

F = F δF (115)

R

1

i+ 2 −

+

δF = δF + δF (116)

33 −

+

|F | |F |,

=

Sottraendo le due espressioni di cui sopra, ricordando che si

ottiene: −

δF = F F (117)

R L

Se si studia un problema lineare come avviene nel caso si voglia risolvere le

Equazioni di Eulero linearizzate, si può scrivere:

∂F ∂F ∂U ∂U

δF = = = A = AδU (118)

∂x ∂U ∂x ∂x

Sostituendo allora dalla (113), si ottiene:

X X X

(k) (k) (k)

δF = AδU = A δω r = Aδω r = λ δω r (119)

k k k k

Inoltre: 1

1 −

+ |δF |] − |δF |]

[δF + δF [δF

=

δF = (120)

2 2

Allora: 1 1 1

− − |δF |] − |A|δU

F = F + [δF = F + AδU F = F (U ) (121)

L L L L

1

i+ 2 2 2

2 1 1

1

+ |δF |] − − |A|δU

− [δF + = F AδU F = F (U )

F = F (122)

R R R

R

1

i+ 2 2 2

2 −

+

|F | |F |,

=

Ricordando la (117) e sottraendo le espressioni riportate sopra

si ottiene: X (k)

− | |A|δU |λ |δω

F F = δ|F = = r (123)

L R k k

k

Si può allora scrivere il flusso sulla generica faccia: 1

1 1

h i

+ − |A|(U −

F F F (U ) + F (U ) U )

= + F = (124)

1 R L

L R

i+ 1 1

i+ i+

2 2 2

2 2 2

Il metodo descritto è upwind in quanto rispetta le direzioni del flusso e poichè

A = cost. il metodo è conservativo. Se il problema non è lineare si può scrivere:

δF = ÃδU (125)

| |

δ|F = Ã|δU (126)

Dove la matrice descrive un problema non lineare. Affinchè il problema

sia conservativo la matrice deve rispettare le condizioni:

• − −

F F = Ã(U , U )(U U ).

Conservatività: Per un problema lineare

R L L R R L

questo è banale in quanto si prende lo Jacobiano come matrice A. Per un

problema non lineare non si ha una condizione generale.

∂F

• →

U = U = U Ã(U ) = J =

R L ∂U

(k)

• λ r

reali e autovettori indipendenti

k 34

5 Roe’s Upwind

In questa sezione si espone uno dei metodi maggiormente utilizzati per costrui-

re la matrice A. Il metodo parte dalla definizione di un nuovo vettore delle

incognite; il Vettore delle variabili di Roe:

√ T

Z = ρ[1, u, v, w, H] (127)

Sia F che U risultano essere componenti quadratiche del vettore Z; con

l’ipotesi di gas perfetto: 2

U

→ −ρ(E −

ρH = ρE + p p = H) E = ρH = z z (128)

1 5

2

2

ρ = z ρu = z z ρv = z z ρw = z z (129)

1 2 1 2 1 3

1

Allora: 1 2 2 2

p = z z + (z + z + z ) (130)

1 5 2 3 4

2

La matrice A può essere quindi calcolata tramite la definizione di alcune

variabili medie alla Roe; per una generica variabile questo si traduce in:

δ(xy) = x̄δy + ȳδx (131)

x +x ȳ.

x̄ = e una formulazione analoga può essere assunta per

in cui L R

2 Ã

Sviluppando i calcoli si può dimostrare che la matrice è identica allo Jacobiano

espresso in funzione delle seguenti variabili mediate alla Roe:

ρ̄ = ρ + ρ (132)

12 L R

i+ √ √

(u ρ) + (u ρ)

L R

ū = (133)

√ √

12

i+ ρ + ρ

L R

... (134)

√ √

(H ρ) + (H ρ)

L R

H̄ = (135)

1

i+ ρ + ρ

2 L R

u possono esser scritte per

Espressioni analoghe a quella ricavata per 1

i+ 2

v w

e . Con tali assunzioni allora la matrice Jacobiana può esser scritta

12 1

i+ i+ 2

alla Roe: Ã(U , U ) = Ā(

Ū ) (136)

1

L R i+ 2

Allora il flusso viene scritto come:

1 1

h i − | −

F = F (U ) + F (U ) Ā|(U U ) (137)

12 L R R L

i+ 2 2

In questa espressione si riconosce una differenza centrata al primo membro

e un termine di dissipazione artificiale al secondo. Si evidenzia quindi ancora

una volta l’equivalenza tra metodi upwind e metodi centrati con dissipazione

artificiale. La matrice A deve esser scritta per ogni faccia.

Un altro schema Upwind molto stabile è AUSM (Advection Upstream Split-

ting Method), di cui tuttavia non si parlerà.

35

6 Schemi Upwind Multidirezionali

Per ottenere uno schema Upwind multidirezionale si procede come esposto per

il caso monodimensionale trattando il flusso per ognuna delle sue componenti

U U

locali di velocità. Resta da definire il valore da assegnare a e . Porre

L R

U = U U = U

e si ottiene uno schema Upwind del primo ordine. Di norma

L i R i+1

tuttavia questo non è sufficiente e si procede con schemi di ordine superiore

coinvolgendo un maggior numero di celle come sarà esposto in seguito. Vale la

pena sottolineare che per prevenire oscillazioni di tipo numerico lo schema deve

essere in grado di passare al primo ordine in presenza di discontinuità. Per fare

U U

questo si definisce ed per estrapolazione inserendo termini che si riducono

R L

o annullano in presenza di discontinuità riducendo l’ordine dello schema.

Approccio MUSCL Per approccio di tipo MUSCL si intende Monotone Up-

wind Scheme for Conservation of Loads. Questo è uno degli approcci più diffusi

U , U

per calcolare le variabili necessarie alla chiusura di uno schema Upwind

L R

in quanto è generale e flessibile. Le variabili locali sono definite come:

θ h i

− +

L −

(1 k)δ u + (1 + k)δ u

U = U + (138)

i i

i

i+1 4

θ h i

R +

− −

U = U (1 + k)δ u + (1 k)δ u (139)

i i i+2

i+1 4

δu

Si ricorda che rappresenta la differenza in avanti o all’indietro a seconda

δ

dell’apice assunto. e k sono parametri che permettono di ricondursi a diversi

schemi di estrapolazione:

Figura 15: Estrapolazione dei valori all’interfaccia nel caso di k=-1

• −1:

k = Estrapolazione lineare decentrata che da luogo a schema Upwind

del secondo ordine

• k = 0: Estrapolazione coinvolge una cella a monte e una a valle, schema

del secondo ordine

• k = 1: Estrapolazione Centrata

12

• k = : Estrapolazione QUICK. Questo approccio si è sviluppato indipen-

dentemente rispetto a MUSCL. Tuttavia i due approcci sono equivalenti

36

13

• k = : Estrapolazione Upwind del terzo ordine

Si inseriscono alcune altre modifiche per evitare oscillazioni attorno ad even-

tuali disocontinuità, dove lo schema deve passare al primo ordine per essere

stabile. Questo viene fatto tramite l’introduzione di delimitatori non-lineari:

• ENO (Essentially Non Oscillatory): Approccio che va a diminuire le celle

U

coinvolte nel calcolo di mano a mano che ci si avvicina a delle di-

12

i+

scontinuità. Approccio ENO viene usato con successo con schemi di ordine

superiore al quarto.

• TVD (Total Variation Diminishing): Approccio consolidato e diffuso, ap-

profondito in seguito.

Approcci TVD Questi metodi prevengono l’insorgere di termini oscillatori

all’interno della soluzione. Tali oscillazioni insorgono in presenza di forti gra-

dienti o discontinuità nella soluzione. Per fare questo agiscono sui meccanismi

di produzione delle oscillazioni a differenza di quanto fanno i metodi con dissipa-

zione artificiale che lasciano che le oscillazioni si manifestino per poi smorzarle.

Approcci TVD sono basati sull’ipotesi che la variazione totale di ogni soluzione

ammissibile non aumenti nel tempo. Questa è una proprietà delle leggi scalari

di conservazione; posto TV come "Total Variation" si ha che:

Z ∂u dx

TV = (140)

∂x

A livello discreto, per il passo temporale n questo si traduce in:

X

n ni+1 ni

|u − |

T V (u ) = u (141)

i

La condizione TVD può essere dunque espressa come:

n+1 n

T V (u ) T V (u ) (142)

In questo modo il valore di un massimo locale non cresce nel tempo ed il

valore di un minimo locale non decresce nel tempo. Uno schema di questo tipo

preserva la monotonicità della soluzione, nel senso che se la soluzione è monotona

n + 1.

al passo temporale n lo sarà anche al passo è possibile dimostrare che se

uno schema TVD è monotono e che uno schema che preserva la monotonicità è

TVD.

Matematicamente allora, se si vuole rendere TVD lo schema MUSCL:

( L ≥

u u

i−1

1

i+

u < u 2

i−1 i+1 R ≤

u u

i+1

1

i+ 2

Tale condizione suggerisce l’introduzione di slope limiters nell’estrapolazione

MUSCL: θ h i

− + +

L +

U = U + (1 k)δ u Ψ + (1 + k)δ Ψ u (143)

i i i

i+1 i−1/2 i+1/2

4 37

con: − −

+ +

Ψ = Ψ(r ) Ψ = Ψ(r ) (144)

i−1/2 i−1/2 i+1/2 i+1/2

1 U U

i+1 i

+ +

r = r = r = (145)

L −

i−1/2 i−1/2 −

U U

r i i−1

i+1/2

Per soddisfare la condizione TVD deve valere:

Ψ(r) − ≤

Ψ(s) 2 (146)

r

Ψ(r) Ψ(r) > 0

E, se si suppone che sia una funzione positiva, dunque se

r > 0 Ψ(r) = 0 r 0,

e se è possibile dimostrare che la condizione TVD diventi:

≤ ≤

0 Ψ(r) 2r (147)

Ψ(r) 2

Se si suppone inoltre la condizione diviene:

≤ ≤

0 Ψ(r) min(2, 2r) (148)

Si impone una ulteriore condizione che fa si che in regioni smooth, dove non

ci sono forti gradienti, il metodo sia almeno del secondo ordine:

( ≤

Se r < 1 r Ψ(r) < 1

Se r = 1 Ψ(r) = 1

Le condizioni discusse possono essere riportate in un grafico:

Ψ(r)

La zona grigia delimita i valori di per cui l’estrapolazione MUSCL è

TVD. All’interno di quest’area sono riportate due funzioni proposte come limiter

molto utilizzate: |r|

r +

V an Leer Ψ(r) = (149)

1+ r

2

r + r

V an Albada Ψ(r) = (150)

2

1 + r

Altri Limiter abbastanza comuni sono il liniter min-mod ed il limiter Super-

bee. 38

Riepilogo A seconda dell’applicazione si utilizzano sia schemi centrali con dis-

sipazione artificiale, sia schemi upwind TVD. In generale gli schemi centrali con

dissipazione artificiale sono più semplici da implementare e di valenza più gene-

rale; funzionano in campi di moto che vanno dall’incomprimibile al supersonico

e sono robusti ed efficienti. Tendono tuttavia ad essere diffusivi ed a introdurre

errori nella soluzione di scie e di strati limite esterni. Ad esempio uno schema

centrale potrebbe risolvere la zona esterna della scia di un profilo aerodinamico

in maniera non fisica, fornendo in output velocità maggiori di quella calcolata

nel free stream. Schemi TVD non consentono la creazione di nuovi massimi e

non soffrono quindi di questo problema.

Schemi Upwind sono costruibili con qualsiasi ordine di accuratezza deside-

rato. Sono tuttavia computazionalmente piò pesanti in quanto l’upwinding va

fatto per ogni passo temporale, a differenza della dissipazione artificiale. Non

sono inoltre facilemente estendibili al caso di gas reale, questo in quanto non è

ne semplice, ne univoco diagonalizzare la matrice Jacobiana nel caso si tratti

un gas reale. 39

Parte V

Tecniche Numeriche di Calcolo

Per risolvere numericamente un problema fisico, discretizzato con i metodi

discussi fino a questo punto, si hanno in generale due possibili approcci:

• Metodi Time Marching: Si considera la problematica di tipo evolutivo

anche per problemi stazionari. In questo modo le equazioni rimangono

sempre iperboliche. In forma stazionaria le equazioni sono ellittiche ne

caso subsonico e iperboliche nel caso supersonico. Si inizializza il calcolo

con una soluzione di primo tentativo, la quale, non rappresentando uno

stato di equilibrio, è di tipo non fisico. Con l’avanzare temporale della

soluzione questa approssima asintoticamente la soluzione stazionaria. Nel

caso in cui si sia interessati alla soluzione stazionaria, discretizzazione

spaziale e temporale sono completamente separate. In pratica è possibile

allungare a piacere il timestep.

• Metodi di Pressure Correction: Utilizzati sopratutto per flussi a basso

Mach, saranno trattati nel dettaglio in seguito.

1 Metodi di Runge-Kutta

I Metodi di Runge-Kutta sono una famiglia di schemi numerici sia di tipo im-

plicito che di tipo esplicito sviluppati intorno ai primi del ’900 dai matematici

Tedeschi Runge e Kutta.

Metodi di Runge-Kutta Espliciti In linea con quanto già detto in (94), si

definisce il residuo R(U) come: dU

dU X

−1 →

= R(U ) = (F + F + D) S = JR(U )

J (151)

conv dif f k k

dt dt

k n

Il metodo di Runge-Kutta permette di passare dal passo temporale al

n + 1

passo seguendo uno schema di questo tipo:

(0) n

U = U (152)

(q) (0) (q−1)

U = U + α ∆tJR[U ] (153)

q

... (154)

n+1 (k)

U = U (155)

k

R α

Questo è uno schema di ordine a k passi. Il valore di dipende sia dal

q

passo che dal numero di passi:

• α = 0, 15 α = 0, 4 α = 1

3-passi: , ,

1 2 3 √

1 1 1

• →

α = α = α = α = 1 CF L = 2 2

4-passi: , , ,

1 2 3 4

4 3 2

In uno schema di Runge-Kutta esplicito si ha stabilità condizionata, si deve

dunque sottostare alla condizione CFL. Ogni passo è costoso computazional-

mente, si deve quindi cercare un compromesso tra robustezza e stabilità. Dove

possibile si usano schemi ibridi: i termini diffusivi e convettivi sono ridiscussi

solo in alcuni passi temporali. 40 n + 1:

Metodi Impliciti Il residuo R viene valutato al passo temporale

dU X X

−1 n+1 n+1

n+1

J = R(U ) = (F + F + D) S = F S (156)

conv dif f k k

k k

dt k k

Linearizzando è possibile scrivere: n

∂F n

n+1 n n δU = F + AδU

F = F + δF = F + (157)

∂U

Allora: dU X X

−1 n

− −

J = F S AδU S (158)

k k

k

dt k k δt

Portando al primo membro e moltiplicando per e J:

X n

[I + Jδt A S ]δU = JδtR(U ) (159)

k k

L’uso diretto di questa equazione corrisponde ad uno schema centrato in-

∂F

stabile. Allora in modo Upwind, questo porta allo splitting della matrice

∂U

Jacobiana: 1 1

+ |A|] − |A|]

A = [A + A = [A (160)

2 2

Schema Upwind del primo ordine garantisce la dominanza diagonale del

sistema risolutivo. Si ottiene quindi sistema della forma:

X X −

+ n −

[I + Jδt A S ]δU = R Jδt A S δU (161)

k k N B

k k

k

N B

Il pedice fa riferimento a grandezze calcolate su celle adiacenti; neigh-

bouring cells. Fattorizzazione approssimata riduce richiesta di memoria:

− − −

+

+ + +

A = A + A + A A = A + A + A A = A + A + A (162)

η ζ η η

ζ ζ

Sostituendo nella relazione ricavata precedentemente:

X X −

+

+ + n −

[I + Jδt (A + A + A ) S ]δU = R Jδt A S δU (163)

k k k N B

η ζ k

k

Si può ricavare una espressione fattorizzata approssimata:

+

+ +

[(I + JδtδA )(I + JδtδA )(I + JδtδA )]δU = RHS (164)

η ζ

ADI Altering Direction Implicit:

+

(I + JδtδA )δ U̇ = RHS (165)

+

(I + JδtδA )δ Ü = δ U̇ (166)

η

+

(I + JδtδA )δU = δ Ü (167)

ζ

Da un sistema completamente tridimensionale ci si è ricondotti con l’espres-

sione (165) a tre sistemi monodimensionali. Ognuna delle tre equazioni risulta

in un sistema tridiagonale che può essere risolta con l’algoritmo di Thomas. Si

ha una precisione del secondo ordine nello spazio e nel tempo. Questo giustifi-

ca anche l’espressione fattorizzata approssimata, in quanto i termini trascurati

sono del secondo ordine. 41

2 Tecniche numeriche per accelerare la Conver-

genza

Riepilogo Si parla in questa sezione di tecniche numeriche per accelerare la

convergenza verso condizioni stazionarie. Si ricorda che in generale i passi da

compiere per la risoluzione di un problema fluidodinamico in modo numerico

sono nell’ordine:

• Discretizzazione Spaziale delle equazioni di governo. Questa può esser

fatta con chemi centrali con dissipazione artificiale o con schemi upwind.

• Discretizzazione nel tempo. Nel caso si utilizzi uno schema implicito si

deve sottostare la condizione CFL.

Le tecniche di accelerazione sono varie e possono essere anche combinate tra

loro. Parleremo delle più comuni.

Local Timestepping Un metodo esplicito come il metodo di Runge-Kutta

deve sottostare alla condizione CFL. Si deve assicurare che il dominio di di-

pendenza fisico ricada entro il dominio di dipendenza numerico. La velocità di

propagazione delle onde è data dagli autovalori dello Jacobiano associato ad una

r

certa direzione. Se si pone allora come raggio spettrale associato alla direzio-

1

∆t =

ne (Autovalore massimo, il massimo timestep vale . Il delta t critico

r

è distinto per i termini diffusivi e per quelli convettivi e vale rispettivamente:

1 1

∆t = ∆t = (168)

c d γµ 2

2 2 2 2 2

r + r + r k J (δ δ + δ + δ + δ + δ )

η ζ t

η η

ζ

ρP r

Per cautelarsi da instabilità si dovrà prendere il minimo tra i due. Conviene

quindi definire una grandezza di riferimento:

∆t ∆t

c d

∆t = (169)

ref ∆t + ∆t

c d

∆t

Si calcola per ogni timestep e per ogni cella di calcolo. Il local time-

ref

stepping utilizza un timestep differente per ogni cella di calcolo; in altre parole

ogni cella di calcolo utilizza il timestep più lungo possibile. Questo permette di

espellere più velocemente le perturbazioni e arrivare prima alla soluzione stazio-

naria. Durante il transitorio la soluzione sarà ovviamente non fisica, ma questo

non è un problema se si è interessati alla soluzione stazionaria e se, a maggior

ragione, la soluzione di tentativo di partenza è essa stessa non fisica. Allora il

timestep locale vale: ∆t = CF L∆t (170)

ref

Residual Smoothing Questa è una tecnica generale per aumentare stabilità

e robustezza di uno schema numerico. Si ricorda che uno schema viene detto

stabile se perturbazioni non sono amplificate da un passo temporale al succes-

sivo. Per velocizzare la convergenza di uno schema si cerca di massimizzare il

CFL. Questo introduce difficoltà nella convergenza dello schema con oscillazioni

dei residui. Residual Smoothing appiattisce i residui, diminuendo le variazioni

tra una cella e l’altra. Si fa in due modi:

42

• Esplicito

• Implicito Ṙ

Con riferimento ai metodi espliciti, sia R il residuo ed il residuo smoothato:

− − +

+ +

− − −

(1 β δ δ )(1 β δ δ )(1 β δ δ ) Ṙ = R (171)

η ζ

η η ζ ζ

+

δ δ equivale ad una differenza del secondo ordine. Questa espressione può

essere fattorizzata: −

+ (1)

(1 β δ δ )R = R (172)

+ (2) (1)

(1 β δ δ )R = R (173)

η η η

+ (3) (2)

(1 β δ δ )R = R (174)

ζ ζ ζ

β può esser scritto come: r

CF L

1 k

, [ ])

β = max(a, (175)

k ˙

4 r + r + r

CF L η ζ

˙

CF L

Come si intende il CFL con residui spianati. In termini pratici questo

risulta essere circa doppio del CFL con residui non spianati. Se applicato ad

uno schema di Runge-Kutta sostituisco al residuo il residuo spianato, per cui al

q-esimo passo del metodo si avrà:

(q) (0) (r−1)

U = U + α H∆t[

Ṙ ] (176)

q

Il modo di procedere può essere riassunto come segue: Si assemblano i flussi,

si avanza con il metodo di Runge-Kutta, si smorzano i Residui e si ripete il ciclo.

Se si ha a che fare invece con un metodo implicito, vale:

X +

[I + J∆t A S ]δU = JδtR (177)

k

k

k

Dividendo le componenti:

X X −

+ −

[I + J∆t A S ]δU = JδtR Jδt A δ δU (178)

k k N B

k k

k k

(n+1) n

δU = U U

Tuttavia è il residuo R. Allora la relazione può essere

riscritta con il residuo smoothato:

X X −

+ −

[I + J∆t A S ] Ṙ = JδtR Jδt A δ Ṙ (179)

k k N B

k k

k k .

Si può passare ad un metodo di Runge-Kutta inserendo il parametro Con

= 0 = 1

il metodo è esplicito, con il passo è implicito.

X X −

+ −

[I + J∆t A S ] Ṙ = JδtR Jδt A δ Ṙ (180)

k k N B

k k

k k

43

Strategia Multi-Build Una perturbazione di lunghezza d’onda elevata ri-

spetto ad una griglia fine se riportata su una griglia di dimensioni maggiori

appare di lunghezza d’onda più fine e viene dissipata più rapidamente. Un

metodo per accelerare la convergenza di analisi stazionarie allora è quello di

costruire più griglie sovrapposte, da sfruttare nei modi che saranno discussi in

h

seguito per espellere più velocemente le perturbazioni. Se è la spaziatura di

una griglia omogenea, ottengo griglie progressivamente più rade con spaziatura

2h, 4h, 8h, . . . . Ogni cella di una griglia contiene 4 celle della grigia più fine.

Allora nel passare da una griglia più fine ad una più rada, la soluzione sulla

griglia rada viene posta uguale alla media pesata sui volumi (o sulle aree nel

caso 2D) delle soluzioni calcolate nelle celle della griglia più fine. Le correzioni

calcolate sulle griglie rade devono esser poi riportate sulla griglia fine che è la

più accurata per il calcolo. All’interno delle strategie multi-build vi sono due

modi di procedere; Cicli a "V" e cicli a "W".

Figura 16: Ciclo a "V"

Figura 17: Ciclo a "W"

Nei cicli a "V" la soluzione è calcolata per la griglia più fine. Si trasferisce

quanto calcolato sulla griglia della dimensione più grande, si avanza la soluzione

e la si trasferisce alla griglia successiva. Questo viene ripetuto fino a che non si

raggiunge la griglia più rada; a tal punto si avanza la soluzione e si torna sulla

griglia di partenza interpolando sulle griglie più piccole la soluzione calcolata

sulla griglia più rada. Quanto descritto è riassunto graficamente in figura 16.

44

Per aumentare la stabilità dei metodi Multi-Build si ricorre a sotto-iterazioni.

Il ciclo di questo tipo più comune è il ciclo a "W" riportato in figura 17.

3 Metodi di pressure-correction

I metodi di pressure-correction sono una famiglia di algoritmi impiegati per la

soluzione delle equazioni di Navier-Stokes specialmente adatti al caso in cui si

trattino flussi incomprimibili o scarsamente comprimibili. In questi casi infatti

la densità è pressochè indipendente dalla pressione e non può quindi essere

utilizzata come grandezza di appoggio per calcolare velocità e pressione stessa.

Algoritmo SIMPLE L’algoritmo SIMPLE (Semi-Implicit Method for Pres-

sure Linked Equations) è il primo e più semplice degli algoritmi di pressure

correction e costituisce la base per la maggior parte degli algoritmi adatti a

risolvere flussi a basso mach. Lo schema di calcolo è riportato in figura 18 per

una griglia di calcolo bidimensionale.

Con un primo valore di tentativo della pressione si va a risolvere le equazioni

della quantità di moto per ottenere i valori delle velocità su una griglia staggered

come in figura 11, tramite le relazioni:

X

∗ ∗ ∗ ∗

a u = a u + (p p )A + b (181)

I,J nb i,J i,J

I,J nb I−1,J I,J

X

∗ ∗ ∗ ∗

a v = a v + (p p )A + b (182)

I,j nb I,j I,j

I,j nb I,J−1 I,J

0 ∗

p p

Definendo un correttore di pressione è possibile correggere la pressione

∗ 0 ∗ 0 ∗ 0

p = p + p u = u + u v = v + v

ponendo: . Analogamente per la velocità e .

Sostituendo tale espressione nella (181) e sottraendo alle equazioni ottenute la

(181) stessa si ottiene: X

0 0 0 0

a u = a u + (p p )A + b (183)

I,J nb i,J i,J

I,J nb I−1,J I,J

X

0 0 0 0

a v = a v + (p p )A + b (184)

I,j nb I,j I,j

I,j nb I,J−1 I,J

A questo punto si introduce la principale ipotesi di modellazione dell’algo-

ritmo SIMPLE. Essa consiste nel supporre nulli i contributi indicati dal pedice

’nb’. Tale semplificazione è accettabile in quanto la relazioni (183) non sono

fisiche, nel senso che sono relazioni scritte a supporto del metodo numerico di

soluzione che saranno identicamente nulle e dunque prive di significato una vol-

d = A /a

ta che l’algoritmo sia giunto a convergenza. Posti allora e

i,J i,J I,J

0 ∗

u = u u si ottiene l’espressione che lega le velocità al correttore di pressione:

∗ 0 0

u = u + d (p p ) (185)

i,J i,J

i,J I−1,J I,J

∗ 0 0

v = v + d (p p ) (186)

I,j I,j

I,j I,J−1 I,J−1

0

p

Il correttore viene ottenuto a partire dall’equazione di continuità discre-

tizzata: − −

[(ρuA) (ρuA) ] + [(ρvA) (ρvA) ] = 0 (187)

i+1,J i,J I,j+1 I,j

Sostituendo all’interno delle equazioni i termini u e v ricavati nella (185) si

ottiene:

0 0 0 0 0 0

a p = a p + a p + a p + a p + b

I,J I+1,J I−1,J I,J+1 I,J−1

I,J I+1,J I−1,J I,J+1 I,J−1 I,J

(188)

45 0

p

Si è ottenuta una equazione in cui l’unica incognita è la variabile . Risolvendo

tale equazione è possibile ricavare dunque il valore del correttore di pressione

da inserire nella (185) per ricavare le nuove velocità e pressioni. Ancora una

volta questa relazione non ha un valore fisico; quanto il sistema è a convergenza

sarà infatti identicamente nulla. I coefficienti presenti sono funzione di quantità

0 0

a b

note. Si riporta l’espressione di uno dei coefficienti e del coefficiente , si

ricavino gli altri per analogia:

a = (ρdA) (189)

I+1,J i+1,J

... (190)

0 ∗ ∗ ∗ ∗

− −

b = (ρu A) (ρu A) + (ρu A) (ρu A) (191)

i,J i+1,J I,j I,j+1

I,J

Le relazioni appena descritte sono combinate nei passaggi descritti in figura

18. Figura 18: Algoritmo SIMPLE

46

Per aiutare la convergenza del metodo si usa spesso un fattore di sottorilas-

samento. La pressione corretta mantiene in pratica memoria del suo valore al

timestep precedente: ∗

new −

p = (1 α )p + α p (192)

p p

α

Il valore di non supera generalmente la soglia di 0.3 e deve esser comunque

p

compreso tra 0 ed 1.

Dall’algoritmo SIMPLE sono stati derivati altri algoritmi di soluzione pres-

sure based; i più diffusi sono:

• SIMPLER: SIMPLE-Revised. L’equazione di continuità viene usata per

ricavare una reazione che ricavi direttamente la pressione per il passo

successivo dell’algoritmo, senza passare da una funzione di pressure cor-

rection. Il campo di velocità è tuttavia comunque ottenuto attraverso una

funzione di correzione della velocità.

• SIMPLEC: SIMPLE-Consistent. Concettualmente identico all’algorit-

mo SIMPLE da cui deriva, tuttavia vengono supposti pari a zero i termini

meno significativi dell’equazione della quantità di moto. Viene general-

mente preferito all’algoritmo SIMPLE standard ed è il solutore pressure

based più comunemente impiegato nei moderni codici CFD.

Algoritmi PISO Gli algoritmi PISO (Pressure Implicit with Splitting of Ope-

rators) è stato originariamente sviluppato per la soluzione di flussi comprimibili

in simulazioni transitorie. Il suo successivo sviluppo lo ha portato ad essere

lo standard nelle simulazioni time dependant di flussi a basso mach. Può es-

ser visto come una estensione dell’algoritmo SIMPLE con l’introduzione di un

ulteriore passo correttivo.

Il primo passo viene detto "predictor step" e corrisponde al risolvere le equa-

zioni discretizzate della quantità di moto con i valori di tentativo per velocità e

pressione. Viene condotto poi un primo "corrector step" e si chiude il ciclo SIM-

∗∗

∗∗ , v

u . A questo punto viene

PLE ottenendo i valori corretti delle velocità I,j

i,J

condotta una sotto-iterazione in più; si ricava un nuovo correttore di pressio-

00

p

ne in maniera analoga a quanto visto per l’algoritmo SIMPLE e tramite le

equazioni della quantità di moto ottengo un nuovo valore per pressioni e velocità:

∗∗∗ ∗∗ 00 ∗ 00 0

p = p + p = p + p + p (193)

∗∗ ∗

P −

a (u u )

nb

∗∗∗ ∗∗ 00 00

nb nb −

+ d (p p )

u = u + (194)

i,J

i,J i,J I−1,J I,J

a

i,J

∗∗ ∗

P −

a (v v )

nb 00 00

∗∗∗ ∗∗ nb nb −

v = v + + d (p p ) (195)

i,J

I,J I,j I,J−1 I,J

a

I,j

Il ciclo è riepilogato in figura 19.

La relazione (194) viene ricavata sottraendo le espressioni della quantità di

∗∗∗ ∗∗

u u

moto scritte per le quantità e ;

∗∗

P a (u )

nb

∗∗∗ ∗∗∗ ∗∗∗

nb −

u = + d (p p ) (196)

i,J

i,J I−1,J I,J

a

i,J ∗

P a (u )

nb

∗∗ ∗∗ ∗∗

nb −

u = + d (p p ) (197)

i,J

i,J I−1,J I,J

a

i,J 47

Figura 19: Algoritmo PISO

48

Ricordando infine la (193) si ottiene infine la relazione (194) che si cercava.

L’algoritmo PISO è generalmente più efficiente del SIMPLE riducendo in alcuni

casi il costo computazionale di un fattore 2. Generalmente è inoltre molto

robusto, soffre tuttavia quando l’inizializzazione della soluzione è molto lontana

dalla soluzione a convergenza e non tollera condizioni CFL troppo elevate.

Interpolazione di Rie-Chow Nella fluidodinamica numerica moderna l’uso

di griglie di tipo staggered è essenzialmente superato. Per evitare l’evoluzione

oscillatoria del calcolo e assicurare l’accoppiamento tra pressione e velocità viene

comunemente utilizzata l’interpolazione di Rie Chow. Con riferimento ad una

griglia monodimensionale come in figura:

La velocità sul bordo e viene calcolata:

u + u 1 1 1

E P − − − − −

u = + (d + d )(p p ) d (p p ) d (p p ) =

e P E P E P W E E P E

2 2 4 4

3

d ∂ p

u + u

E P 3

+ ∆x

= (198)

3

2 4 ∂x e

La velocità viene calcolata introducendo un gradiente di pressione del terzo

ordine. Tale procedura può esser facilmente adattata a situazioni multidimen-

sionali e a griglie non strutturate.

4 Solutori Segregati e solutori Coupled

I metodi pressure based per la risoluzione dei flussi a basso mach possono essere

implementati in due modi:

• Solutore Segregato: Si risolve equazione della quantità di moto e continui-

tà con uno degli algoritmi di pressure correction visti in precedenza. Si

risolvono successivamente le relazioni che determinano tutti gli altri sca-

lari quali turbolenza o altri. La semplificazione introdotta è giustificata

dal fatto che, se la densità è pressoché costante, l’accoppiamento tra le

equazioni di Navier-Stokes è blando e disaccoppiando le equazioni non si

introducono eccessive distorsioni nella fisica del problema.

• Solutore Accoppiato: Solutori Coupled risolvono il sistema di equazioni dif-

ferenziali alle derivate parziali di Navier-Stokes mantenendo accoppiamen-

to fisico tra le equazioni. Si risolve il sistema completo ad ogni iterazione

senza disaccoppiare tra loro le variabili. Questi metodi sono consigliati nel

caso di flussi ad alto mach in quanto la densità introduce un forte accop-

piamento tra le equazioni. Anche metodi pressure based possono essere

adattati ad essere risolti in maniera accoppiata.

In generale metodi Coupled impiegano un minor numero di iterazioni per

arrivare a convergenza, tuttavia richiedono più memoria ad ogni iterazione.

49

5 Trattamento numerico dei termini sorgente

I termini sorgente corrispondono nelle equazioni fisiche di partenza alle forze di

volume nel caso in cui si tratti le equazioni della quantità di moto o a generazio-

ne di energia all’interno del dominio considerato qualora si tratti dell’equazioni

dell’energia. I termini sorgente possono però esser modificati in modo da mo-

dellare effetti fisici che non si vuole studiare nel dettaglio ma dei cui effetti si

è interessati. Un flusso reattivo in cui avviene combustione ad esempio può

esser modellato con un termine sorgente nella equazione dell’energia. Il flusso

attraverso un condotto molto piccolo, e quindi spesso molto oneroso da discretiz-

zare, può esser modellato con un termine sorgente nell’equazione di continuità.

Questo consente di modellare correttamente la portata che passa attraverso il

condotto. Si potrebbe inoltre introdurre un termine sorgente nella equazione

dell’energia per modellare eventuali effetti di scambio termico.

Con riferimento alla equazione (92), considerando una generica variabile

φ,

trasportata si riprende il termine sorgente, che viene spesso discretizzato:

Z S dV = S V (199)

φ φ

V S

esplicita

Una discretizzazione presuppone che il termine sia costante.

φ

S = S (φ).

implicita

Una discretizzazione pone Una discretizzazione esplicita

φ φ

del termine sorgente è accettabile qualora esso sia pressochè costante. Se questo

termine varia in maniera significativa deve esser preferita una discretizzazione

implicita. Tipicamente viene impiegato il metodo di Picard:

S(φ) = S + S φ (200)

c p p

S S

Dove è costante mentre è la parte implicita.

c p ∗

φ

Il termine sorgente viene spesso linearizzato; ponendo come il valore della

p

generica variabile trasportata alla interazione corrente, espandendo in serie di

Taylor: ∂S

∗ −

(φ φ )

S = S + (201)

p p

∂φ

Allora ricordando la (200), è possibile ricavare:

∂S ∂S

∗ ∗

S = S φ S = (202)

c p

p

∂φ ∂φ

Il metodo di Picard linearizza il termine sorgente, assumendo il termine

S (φ).

sorgente stesso pari ad una retta tangente alla curva φ

6 Cenni a metodi di soluzione non-stazionari

Per flussi ad alto mach, risolti con solutori density based di tipo time marching,

la metodologia di soluzione rimane la stessa sia che si tratti flussi stazionari, sia

che si trattino flussi unsteady. Nel caso di un calcolo non stazionario si dovrà

prestare attenzione al time step scelto che deve garantire adeguata risoluzione

temporale con accettabili costi computazionali.

Nel caso in cui si trattino solutori pressure based vi sono delle differenze

formali nel trattare un caso stazionario rispetto ad uno non stazionario. I metodi

50

visti quali gli algoritmi SIMPLE e PISO utilizzano formulazioni implicite in un

problema non stazionario. Intuitivamente si va a ripetere la procedura iterativa

degli algoritmi SIMPLE o PISO ad ogni timestep. Quando si è raggiunta la

convergenza il solutore avanza al timestep successivo ripentendo la procedura.

Il vantaggio di un solutore PISO rispetto ad un SIMPLE in calcoli non sta-

zionari sta nel fatto che, poiché nel PISO viene introdotta una sottoiterazione

aggiuntiva, è possibile dimostrate che la sua accuratezza temporale sia del terzo

ordine per la pressione e del quarto ordine per la quantità di moto. Questo

in termini pratici significa che un solutore di questo tipo con un timestep ade-

guatamente breve non necessita di iterare ad ogni passo temporale in quanto

il valore prodotto dopo una singola iterazione è sufficientemente accurato nella

gran parte dei casi. Usare un solutore SIMPLE costringe il calcolatore ad iterare

fino a convergenza ad ogni timestep, sebbene possano essere scelti per questo

motivo dei passi temporali leggermente più lunghi.

Il tempo che si impiega per risolvere una simulazione transitoria è ovvia-

mente funzione del caso che si studia. Nel caso si vogliano studiare fenomeni

transitori di un caso studio globalmente stazionario, quale può essere lo studio

di dettaglio di moto in un condotto, il tempo totale necessario ad ottenere una

buona simulazione può esser stimato come:

∆x

∆T = n (203)

u

n è il numero di volte che il dominio deve essere attraversato, tipicamente

non meno di 4-5 volte. 51

Parte VI

Modelli di Turbolenza

1 Fenomenologia

La turbolenza è un fenomeno di moto caotico di un fluido che si verifica quando

gli effetti convettivi sono di gran lunga maggiori degli effetti viscosi. Si verifica

quindi ad alti numeri di Reynods ed è uno stato naturale del fluido. Il feno-

meno è fortemente tridimensionale e intrinsecamente non stazionario, se si fa

una simulazione bidimensionale si perdono necessariamente alcuni effetti. La

turbolenza è il metodo più efficiente di trasferire energia meccanica in calore.

Lo strato limite di un fluido è una zona caratterizzata da forte turbolenza;

risulta quindi chiaro che per trattare adeguatamente le pareti in un calcolo

CFD ci si trova ad avere a che fare con questo fenomeno. Poichè inoltre si

hanno pareti solide nella pressochè totalità dei problemi ingegneristici, si farà

spesso riferimento al trattamento a parete dei vari modelli di turbolenza.

2 Energy Cascade

Per descrivere la turbolenza è importante inquadrare l’energia associata alla

varie scale della turbolenza stessa. Questo permetterà di capire da un punto di

vista qualitativo come evolve il fenomeno. Si parla di Energy Cascade:

Figura 20: Energy Cascade

La turbolenza viene generata in corrispondenza delle scale più grandi, ca-

ratterizzate da numero d’onda minore. Qui formano vortici con diametro medio

pari alla lunghezza caratteristica del problema in studio (diametro dello swirler,

altezza palettatura etc . . . ), gli effetti viscosi tenderanno poi a dissipare e rom-

pere questi vortici fino alla completa dissipazione e annullamento della velocità

turbolenta che avviene in corrispondenza delle scale più piccole. Andando nel

dettaglio delle zone in figura 20: 52

• Zona I : I vortici grandi sono generati ricevendo energia dal flusso medio.

Questi danno origine a vortici più piccoli trasferendovi energia. La strut-

tura dei vortici è fortemente dipendente dalla geometria del domino e per

questo la turbolenza in questa zona è fortemente anisotropa.

• Zona II: Inertial Subrange. In questa zona la turbolenza è Isotropa. Si ha

il trasferimento di energia da vortici grandi a vortici di scala minore. Su

un grafico bi-logaritmico la pendenza è lineare.

• Zona III: Dissipation Range: Le scale minori sono dominate dagli effetti

viscosi e l’energia è dissipata in calore dalla viscosità molecolare. Si è

in prossimità della scala di Kolmogorov. La modellazione fisica basata

sulle ipotesi di meccanica del continuo da cui sono derivate le equazioni di

governo resta valida in quanto anche la scala di Kolmogorov resta diversi

ordini di grandezza più grande del cammino libero medio delle particelle.

La turbolenza è quindi un fenomeno evolutivo in cui si ha il graduale pas-

saggio da macro-vortici a micro-vortici. La scala di Kolmogorov è la più piccola

scala in cui si possono osservare fenomeni turbolenti. Per scale spaziali minori

la turbolenza è completamente dissipata a causa della prevalenza dei termini

viscosi. La scala di Kolmogorov è definita ocme:

3 14

ν

l = (204)

k

Alle scale spazioali si possono associare anche scale temporali. La scala di

Kolmogorov vale: 3 12

ν

τ = (205)

k

Le scale turbolente che si riesce a vedere dipendono dalla dimensione della

griglia. È possibile collegare il numero di Reynolds al numero di celle necessarie

a risolvere tutte le scale turbolente:

L L 9

3 3

→ ' '

'

N = = N = N N N N Re

Re (206)

4 4

x y z

x x

∆x l k

Tipicamente ci si trova a dover risolvere problemi con numero di Reynolds

che può facilmente superare il Milione. Questo porta ad un numero di celle

difficilmente gestibile con la potenza di calcolo disponibile. Il problema deve

allora esser semplificato. Questo viene fatto mediando nel tempo le equazioni

di Navier stokes ed introducendo modelli di turbolenza.

3 Equazioni di NS mediate alla Reynolds

u (x, t)

Presa una generica variabile questa può essere scomposta nella somma

i 0

u (x, t) = U (x, t)+u (x, t).

di un termine medio e di un termine fluttuante come i i i

Il termine medio vale: t+T

Z

1

U (x) = Ū (x) = lim u (x, t) dt (207)

i i i

T

→∞

T t

53

Se si sostituisce nella generica equazione di trasporto

∂u ∂p

ρ =

+ div(ρuu) + div(µ grad(v)) (208)

∂t ∂x

0

∂u

∂U ∂p

i 0 0

0

i

+ ] + div(ρ(U + u )(U + div(µ grad((U + u )))

ρ[ + u )) =

i i

i i

i i

∂t ∂t ∂x (209)

Facendo una media su un tempo abbastanza lungo si annullano i termini

fluttuanti e si ottiene:

∂ Ū ∂p

i 0

0 ū

ρ + div(ρ(

Ū ) + div(ρ(

ū ) = + div(µ grad( Ū ))

Ū (210)

i

i i

j j

∂t ∂x

Sviluppando i calcoli si ottiene:

2 ∂u w ∂u v

∂ Ū ∂u ∂p

i i i i

i i − −

+ div(ρ(

Ū Ū + div(µ grad( Ū ))

ρ ) = + (211)

i i

j

∂t ∂x ∂z ∂y ∂x

| {z }

=Sf orzi di Reynolds

Come si può vedere non si annullano i termini fluttuanti quadratici. Poichè

le equazioni della quantità di moto sono 3, sono introdotte 9 nuove incogni-

−ρ

τ = u ¯

u

te, riassunte nel tensore detto Tensore degli sforzi di Reynolds.

i,j i j

Poichè il tensore è simmetrico le incognite sono ridotte a 6.

Per chiudere il problema allora si deve allora modellare i termini che com-

pongono il tensore degli sforzi di Reynolds tramite equazioni algebriche o dif-

ferenziali. Questo diene fatto tramite un modello di turbolenza. Si è ricavato

un set di equazioni mediate alla Reynolds o modello RANS (Reynolds Averaged

Navier Stokes Equations).

4 Modelli di Turbolenza

I modelli di turbolenza applicati ad un set di equazioni RANS, permettono di

modellare gli effetti della turbolenza sul moto medio del fluido. Vi sono due

macro approcci per chiudere i problemi:

• Eddy Viscosity Models: Modelli possono essere Algebrici, ad 1 equazione,

a 2 equazioni di trasporto. La modellazione dei termini turbolenti passa

dalla definizione della Eddy viscosity

• Stress Transport Models: Possono essere sia di tipo algebrico che di tipo

differenziale. Non si passa dal calcolo della eddy viscosity per chiudere

il modello. Sono in generale più accurati di modelli eddy viscosity ma il

costo computazionale è talvolta molto superiore.

5 Modelli Eddy Viscosity

Questo approccio è stato per primo introdotto da Businnesq. Si basa sul fatto

che gli sforzi turbolenti siano espressi tramite espressioni del tutto analoghe a

quelle per moto laminare con l’eccezione che il coefficiente di viscosità sia una

54

proprietà del moto invece che una proprietà del fluido. La possibilità di

esprimere gli sforzi di taglio in funzione della viscosità turbolenta vale se:

l m

N umero di Knundsen = K = << 1 (212)

n L

l

Dove è il cammino libero medio delle particelle e L è la lunghezza carat-

m

teristica del problema. Questa assunzione è giustificata fisicamente dal fatto

che si sta considerando un problema di meccanica del continuo, implicitamente

quindi si assume che la porzione di fluido considerata sia continua, e dunque

molto più grande del cammino medio libero della singola particella. Questa

assunzione viene meno in casi molto particolari, quali il rientro in atmosfera di

una navicella spaziale, dove l’atmosfera è molto rarefatta.

Le componenti del tensore degli sforzi di taglio sono scritte come:

1 dU

dU = ρl v

τ = µ (213)

mix mix

xy T dy 2 dy

Il problema del calcolo del tensore degli sforzi di Reynolds è adesso spostato

sul calcolo della viscosità turbolenta. La questione può essere affrontata in vari

modi come si vedrà nel seguito:

Modelli algebrici Sono i modelli più semplici e meno onerosi computazio-

nalmente. Funzionano in modo adeguato su un limitato numero di casistiche. Il

comportamento turbolento del fluido è modellato con equazioni algebriche, senza

tener conto del trasporto delle grandezze turbolente. Riferendosi all’espressione

(213), si assume: dU

·

v = const l (214)

mix mix dy

Allora la viscosità turbolenta risulta essere funzione soltanto della lunghezza

di mescolamento, che deve essere dunque calcolata in modo adeguato.

Prima di proseguire è utile fare alcune considerazioni sullo strato limite. Con

riferimento alla figura 21, riportando in ascissa la distanza adimensionalizzata

+ +

y U

dalla parete ed in ordinata la velocità adimensionalizzata , si può notare

una distinzione in tre zone:

• Sublayer: Strato più vicino alla parete in cui le velocità sono più basse.

+

y < 5

Tipicamente

• Log Layer: Strato intermedio. Su un grafico bilogaritmico questa zona ha

+ 3

15 < y < 10

andamento lineare.

• Defect Layer: Strato più lontano dalla parete in cui la velocità è ormai

3 +

10 < y < δ

prossima a quella del free stream.

Uno dei primi modelli algebrici proposti è quello di Pradtl. Questo modello è

stato tarato sul caso di flusso su lastra piana. Sono stati poi proposti in seguito

numerosi altri modelli per migliorare l’accuratezza generale e la generalità del

modello di Pradtl. Il modello di Prandtl assume:

 2 −

Sub Layer

y

 −

l = ky Log Layer

mix  −

δ(x) Def ect Layer

 55

Figura 21: Strato Limite adimensionalizzato

Il principale limite dei modelli algebrici sta nel defect layer. Fintanto che si

vuole risolvere la zona vicina alla parete infatti è sensato assumere la lunghezza

di miscelamento proporzionale alla distanza da parete come fa Prandtl. Se

la distanza da parete è piccola, la lunghezza di miscelamento è piccola e la

viscosità turbolenta è bassa. Non è mai stata trovata una relazione soddisfacente

per calcolare la lunghezza di miscelamento nello strato limite esterno. Inoltre

modelli algebrici in generale soffrono nel caso di gradienti avversi di pressione

per flussi in forte non equilibrio come nel caso di strato limite separato.

Per poter parzialmente superare queste problematiche sono stati sviluppati

modelli algebrici a due equazioni. Il primo modello di questi tipo che è stato pro-

Clebeci-Smith,

posto è il modello di secondo cui la lunghezza di miscelamento

viene calcolata come: +

−y +

l = ky[1 e A ] Inner Layer (215)

mix ˙

µT = αρU δ F (y, δ) Outer Layer (216)

ν Kleb

0 e

Per lo strato esterno non si passa dal calcolo della lunghezza di miscelamento

ma si calcola direttamente la viscosità turbolenta. Rimangono nel caso di stra-

to limite separato difficoltà a calcolare l’altezza dello strato limite turbolento

˙

δ

adimensionalizzata rispetto alla viscosità cinematica sopratutto nel caso di

ν

F

flussi complessi o di strati limite separati. è funzione della distanza da

Kleb +

δ. A

parete y e della altezza totale dello strato limite k,α, sono coefficienti di

taratura. Per quanto riguarda lo strato limite interno il calcolo della viscosità

turbolenta è svolto a partire dalla lunghezza di miscelamento come

12

2 2

∂U ∂V

i

h

µ = ρl +( (217)

T mix ∂y ∂x

Baldwin-Lomax

Il modello algebrico a due equazioni di non necessita del

calcolo di parametri integrali come l’altezza dello strato limite ed è stato pro-

posto proprio per superare le limitazioni del modello di Clebeci-Smith. La

lunghezza di miscelamento viene calcolata in maniera analoga a quanto visto

per il modello di Clebeci-Smith. Allora la viscosità turbolenta viene calcolata

56

come: ( |ω|

ρl InnerLayer

mix

µ =

T ραF (F , y )F (y, y ) OuterLayer

wake max max Kleb max

y F F

è la distanza da parete a cui si ha . è il massimo del prodotto

max max max F

tra la lunghezza di miscelamento e il modulo della vorticità. La funzione max

può avere dei massimi locali. In particolare potrebbe avere un picco nel sub

F

layer. Allora se si inizia a cercare il picco di dalla parete verso l’esterno si

max

potrebbe prendere il massimo sbagliato e sbagliare completamente lo spessore

dello strato limite. Viceversa se si inizia dall’esterno verso l’interno e si ha

strato limite separato ancora una volta si andrebbe a prendere un massimo

locale sbagliato. Nel caso di flusso separato infatti si hanno due picchi, uno

nello strato limite separato e uno nel log layer.

Modelli algebrici vanno inoltre in difficoltà quando il campo di moto si com-

plica con, ad esempio, la presenza di scie e di interazioni tra diversi profili

aerodinamici. Ad esempio nel caso di interazione statore-rotore si avrà la scia

dello statore che interagisce con il dominio di calcolo del rotore. Se nelle zone

di scia si costruisce la viscosità turbolenta come il prodotto della distanza da

parete per la vorticità, poichè la distanza da parete può esser talvolta anche

elevata, si ottengono dei valori del tutto fuori scala.

D’altra parte i modelli algebrici hanno anche dei pregi; sono molto semplici

da applicare e molto stabili. Per applicazioni specifiche, simili ai casi su cui i

modelli sono stati tarati, o per configurazioni di moto particolarmente semplici

producono risultati attendibili con sforzi molto ridotti.

Modelli ad una equazione Complicando drasticamente la questione rispet-

to ai modelli algebrici, con l’intenzione di migliorare la predittività e il campo

di applicabilità dei modelli sono stati sviluppati modelli ad una equazione di

trasporto. L’equazione di trasporto è una equazione differenziale alle deriva-

te parziali e viene scritta tipicamente per una grandezza legata alla viscosità

turbolenta. La scelta più naturale, sulla quale molti modelli si sono orientati, è

quella di slegare la lunghezza di miscelamento dalla velocità di miscelamento. La

scelta della variabile trasportata ricade quindi sulla energia cinetica turbolenta

specifica: 1 ¯ ¯ ¯

02 02 02

(

u + v + w )

k = (218)

2 µ =

L’energia cinetica turbolenta risulta legata alla viscosità turbolenta T

2

1 0

−ρu

ρl k τ =

e alla traccia del tensore degli sforzi di Reynolds . Il mo-

2

mix ii i

dello quindi richiede per la sua chiusura che sia calcolata anche la lunghezza di

miscelamento in modo adeguato. L’equazione di trasporto allora:

∂k ∂k ∂U ∂ ∂k 1

h i

i 0 0 0 0 0

− − −

ρ + ρU = τ ρ + µ ρu u u p u (219)

j ij i j i j

∂t ∂x ∂x ∂x ∂x 2

j j j j

V ar. temp. + Dif f. = P roduzione Dissipazione + Convezione (220)

All’interno del termine convettivo si riconosco rispettivamente un termine

di Diffusione Molecolare, un termine di Trasporto Turbolento e un termine di

diffusione di pressione. In questi ultimi due termini compaiono termini qua-

dratici fluttuanti che devono quindi essere modellati per chiudere il problema.

57

La modellazione di questi termini è fatta in modo semi empirico basandosi su

analisi dimensionale.

Rispetto a modelli algebrici si introducono effetti non locali e storici del

flusso turbolento in quanto si ha una equazione di trasporto. Resta tuttavia il

limite di dover definire la lunghezza di scala in quanto la viscosità turbolenta

viene calcolata come: 1 l

µ = ρk (221)

2

T

Il problema della definizione della lunghezza di mescolamento è stato risolto

Spalart-Allmaras.

da due ricercatori Boing che hanno creato il modello di ν̃

Questo modello risolve una equazione di trasporto per una variabile diretta-

mente legata alla viscosità turbolenta:

µ = ρν̃f (222)

T v1

Il modello include 8 coefficienti di chiusura e 3 funzioni smorzanti dalle quali

f

si ricava tra le altre cose anche . Il modello è molto leggero e decisamente ro-

v1

busto. I risultati sono in generale più attendibili di un modello ad una equazione

basato sul trasporto di energia cinetica turbolenta. Nel campo dell’aerodinami-

ca esterna questo modello è diventato lo standard. Si hanno difficoltà a trattare

flussi fortemente swirlati.

Modelli a due Equazioni Con l’introduzione di un’ altra equazione di tra-

sporto non è necessario specificare a priori nessuna grandezza per chiudere il

problema. Nella stragrande maggioranza dei casi la prima variabile trasportata

è l’energia cinetica turbolenta k. Per la seconda variabile trasportata sono state

proposte varie soluzioni; le più comuni sono la dissipazione turbolenta specifica

e la dissipazione turbolenta specifica per unità di energia cinetica turbolenta

ω. Queste due variabili danno origine assieme alla energia cinetica turbolenta

− −

k k ω.

ai modelli e

Non vi è alcun fondamento teorico nel considerare che la viscosità turbolenta

si funzione di parametri quali k o o tantomeno nell’ipotesi di Bousinnesq

stessa. Pertanto anche i modelli a due equazioni non hanno validità generale

e ci si aspetta accuratezza basse per alcune casistiche di flusso. I modelli a

due equazioni sono largamente diffusi per via della loro maggiore accuratezza e

versatilità e della sempre maggiore disponibilità di potenza di calcolo, sebbene

richiedano approssiamativamente il trenta percento di risorse di calcolo in più di

un modello algebrico. A causa della loro vasta diffusione saranno allora trattati

k− k−ω

separatamente i modelli e cercando di specificare le loro caratteristiche

e principali differenze.

k .

Modelli Le equazioni di trasporto sono scritte per k ed La lunghezza

di scala non è una delle variabili trasportate ma può esser ricavata dal rapporto

.

tra k ed Le due equazioni di trasporto sono:

∂k ∂k ∂U ∂ ∂k

h i

i −

ρ + ρU = τ ρ + (µ + µ /σ ) (223)

j ij T k

∂t ∂x ∂x ∂x ∂x

j j j j

2

∂ ∂ ∂U ∂ ∂k

h i

i −

ρ + ρU = C τ C + (µ + µ /σ ) (224)

j 1 ij 2 T k

∂t ∂x k ∂x k ∂x ∂x

j j j j

58

Sono evidenti della seconda equazione i coefficienti di chiusura. Il termine

∂U

C τ modellato in modo esatto sarebbe troppo oneroso dal

di produzione i

1 ij k ∂x j

punto di vista delle risorse di calcolo. Viene quindi in genere ricavato per ana-

lisi dimensionale. La taratura del modello è stata fatta su lastra piana senza

gradienti di pressione.

k

Il modello presenta notevoli criticità:

• Non può essere integrato fino a parete. La risoluzione del viscous sublayer

deve essere fatta tramite l’introduzione di Wall-Function. La griglia di

calcolo deve essere strutturata in modo da includere l’intero substrato

laminare all’interno della prima cella di calcolo. All’interno della prima

cella di calcolo si utilizza un modello algebrico per riprodurre l’andamento

delle variabili. Questo è vantaggioso in termini di risorse di calcolo in

quanto viene fortemente ridotto il numero di celle presenti vicino alla

parete.

• Il modello non risente di gradienti avversi di pressione. In presenza di gra-

dienti avversi di pressione la modellazione è quindi spesso insoddisfacente.

• Difficoltà nel riprodurre correttamente il flusso separato. La separazione

viene posticipata o in alcuni casi non prevista. Calcolo di scie è in generale

inaccurato.

Sono stati proposti approcci low-Reynolds per migliorare la predizione che

si ottiene usando le wall-function. Quando si parla di low-Reynolds si fa riferi-

µ

Re = T .

mento al numero di Reynolds Turbolento T µ

L’approccio low-Reynolds si basa sulla sostituzione di alcune delle costanti

di calibrazione con delle funzioni di smorzamento che garantiscano il corretto

comportamento a parete. La griglia di calcolo deve essere infittita in prossimità

della parete in quanto deve ora essere risolto il substrato laminare.

k−

Uno dei pregi del modello è la sensibilità alla turbolenza nel free-stream.

Viene modellato in questo caso abbastanza bene ciò che succede nella realtà al

variare del livello di turbolenza nel flusso indisturbato.

k ω

Modelli Modelli di più recente affermazione. Cercano di risolvere le

k ω

limitazioni del modello riuscendoci solo in parte. L’equazione per fu

originariamente formulata da Kolmogorov sulla base di considerazioni dimen-

sionali e poi rivista da Wilcox. Le equazioni di trasporto sono simili a quelle

viste in precedenza:

∂k ∂k ∂U ∂ ∂k

h i

i −

ρ + ρU = τ βρkω + (µ + µ /σ ) (225)

j ij T k

∂t ∂x ∂x ∂x ∂x

j j j j

∂ω ∂ω ω ∂U ∂ ∂ω

h i

i 2

ρ + ρU = τ α βρω + (µ + µ /σ ) (226)

j ij T k

∂t ∂x k ∂x ∂x ∂x

j j j j

k ω

Un problema serio parzialmente risolto dal modello è quello delle

,

condizioni al contorno su che non è in generale nota. La condizione al contorno

ω

su è invece nota ma vale infinito se la distanza a parete è nulla e non può

quindi essere posta in maniera esatta in un solutore numerico. Se si pone una

ω

corretta condizione al contorno a parete, imponendo il più grade valore di

59

compatibile con il condizionamento del sistema, anche la versione standard del

modello, comunemente detta "high-Reynolds" è integrabile fino a parete.

"low-Reynolds"

Il modello standard è stato modificato nella versione tra-

mite l’introduzione di tre funzioni di smorzamento. Questa versione del modello

k ω è in grado di prevedere la transizione dello strato limite da laminare a

turbolento, sebbene in maniera approssimata basandosi su considerazioni empi-

riche. In altre parole la transizione viene riprodotta senza riprodurre i fenomeni

fisici che la originano. Le relazioni vengono modificate aggiungendo ai termini di

P P

produzione dei coefficienti moltiplicativi e . Attraverso la trasformazione

ω k

P P

di Blausius e sono esprimibili in funzione del numero di Reynolds:

ω k 2

∂U/∂η

α

Re 1

P = (227)

x

k β W 2

αα̇ ∂U/∂η

P = Re 1 (228)

ω x

β W

P P

Quando e cambiano segno si avrà l’inizio e la fine della transizione.

k ω

Sono noti allora i Reynolds di inizio e fine della transizione e si può dunque risa-

P P

lire alla posizione della transizione. Se cambia segno prima di il flusso non

ω k

transisce, in quanto la produzione di turbolenza è smorzata. La produzione di

energia cinetica turbolenta deve essere maggiore della sua dissipazione durante

l’innesco dello strato limite turbolento. Una volta terminata la zona di tran-

sizione i due contributi si bilanciano e la produzione è pari alla distruzione di

turbolenza. Prevedere la transizione in presenza di gradienti favorevoli di pres-

sione è molto difficile in quanto questa avviene in una zona spaziale abbastanza

ampia dello strato limite. La transizione in presenza di gradienti favorevoli è

graduale, al contrario di ciò che avviene in presenza di gradienti avversi in cui

la transizione avviene bruscamente.

k−ω

Il modello si comporta bene in presenza di gradienti avversi di pressione.

Inoltre non necessita di wall function per risolvere lo strato limite prossimo alla

parete. Questo ha per contro il difetto di richiedere maggiori risorse di calcolo

k .

rispetto al modello

k ω

Modelli SST I modelli SST (Shear Stress Transport) cercano di com-

− − −

k ω k . k ω

binare i pregi dei modelli con quelli dei modelli Il modello è

infatti più adatto a trattare il flusso vicino alla parete per i motivi discussi in

k

precedenza, il modello invece si adatta molto bene a trattare la turbolenza

nel free-stream. La formulazione del modello SST è dovuta a Menter che riscrive

ω F

l’equazione per in funzione di e definisce una funzione che assume valore

1

unitario vicino a parete e tende a zero verso il bordo dello strato limite. Allora

a partire dall’equazione (223), è possibile ricavare la (225) con l’aggiunta di un

termine di Diffusione:

∂ω ∂ω ω ∂U ∂ ∂ω

h i

i 2

−βρω

ρ +ρU = τ α + (µ+µ /σ ) +(1−F )CD (229)

j ij T k 1

∂t ∂x k ∂x ∂x ∂x

j j j j F

Il termine CD (Cross-Diffusion) permette tramite il coefficiente di passare

1

− −

k k ω

dal modello al modello e viceversa. Menter introduce inoltre un

F

termine per limitare la viscosità turbolenta. Il modello di Menter è il modello

2

più diffuso e in mote applicazioni il più accurato. Nel 2008 Wilcox ha introdotto

k ω

una modifica al modello portandolo a livelli di predittività equivalenti

60


PAGINE

83

PESO

6.82 MB

PUBBLICATO

4 mesi fa


DETTAGLI
Corso di laurea: Corso di laurea in ingegneria meccanica (FIRENZE, PRATO)
SSD:
Università: Firenze - Unifi
A.A.: 2018-2019

I contenuti di questa pagina costituiscono rielaborazioni personali del Publisher France_Papi di informazioni apprese con la frequenza delle lezioni di Fluidodinamica Numerica per Applicazioni Industriali e studio autonomo di eventuali libri di riferimento in preparazione dell'esame finale o della tesi. Non devono intendersi come materiale ufficiale dell'università Firenze - Unifi o del prof Pacciani Roberto.

Acquista con carta o conto PayPal

Scarica il file tutte le volte che vuoi

Paga con un conto PayPal per usufruire della garanzia Soddisfatto o rimborsato

Recensioni
Ti è piaciuto questo appunto? Valutalo!

Altri appunti di Corso di laurea in ingegneria meccanica (firenze, prato)

Compiti svolti di fisica tecnica
Esercitazione
Turbomacchine - Elementi di base
Dispensa
Comportamento Meccanico dei Materiali
Appunto
Sistemi di raffreddamento - Cicli termodinamici
Dispensa