INDICE

Introduzione ………………………………………………………………………….6

1 Le equazioni fluidodinamiche…………………………………………………..8

1.1 Descrizione matematica del flusso……………………………………...8

1.2 Equazioni di Eulero……………………………………………………14

1.3 Sistema completo delle equazioni di Navier-Stokes…………………..14

1.4 Turbolenza……………………………………………………………..19

1.5 Modelli di turbolenza………………..…………………………………21

1.5.1 Modelli Eddy-Viscosity………………………………………...21

1.5.2 Modello ad un equazione di Spalart- Allmaras………………….24

1.5.3 Modello RSM…………………………………………………...27

2 Il codice SU2…………………………………………………………………..28

2.1 Struttura e classi del codice…………………………………………….29

2.2 griglia di calcolo………………………………………………………..31

3 Simulazioni…………………………………………………………………….34

3.1 Simulazione con fluido non viscoso…………………………………….34

3.1.1 Tempi di calcolo e numero d‟iterazioni…………………………..46

3.2 Simulazioni con fluido viscoso……………………...………………….47

3.2.1 Griglie di calcolo………………………………………………….48

3.2.2 Tempi di calcolo e numero d‟iterazioni...…………………………67

4

4 Conclusioni……………………………………………………………….……69

Bibliografia……………………………………………………………….……71

Sitografia…………………………………………………………………........71

5

6

Introduzione

La fluidodinamica computazionale meglio conosciuta con l‟acronimo CFD

(Computational Fluid Dynamics) è la tecnica che permette di determinare una

soluzione numerica per le equazioni che governano un flusso fluido.

Il principale utilizzo della CFD è quello di risolvere le equazioni di Navier-Stokes e

le equazioni ad esse collegate. La risoluzione per via analitica di queste equazioni è

fattibile solamente in casi semplici con flussi laminari, e geometrie semplici (sfere,

lastre piane), mentre le risoluzioni di casi reali, in cui compaiono di frequente flussi

turbolenti, richiedono per forza un approccio numerico. Esistono quindi diversi

metodi per risolvere le equazioni di Navier-Stokes, e generalmente sono operazioni

dall'elevato costo computazionale.

Dagli anni settanta ad oggi c'è stata una notevole crescita nello sviluppo e

applicazione di CFD perché le potenze di calcolo necessarie ad un analisi di questo

tipo sono divenute nel corso degli anni di domino pubblico, lo sviluppo a interessato

tutti gli aspetti della dinamica dei fluidi che conduce al CFD, diventando uno

strumento di modellazione standard ampiamente utilizzato nell'industria in quanto

può essere considerato una valida alternativa alla sperimentazione fisica in galleria

del vento.

Il presente lavoro di tesi si pone l‟obiettivo di valutare le capacità del software libero

SU2, nell‟ottenere risultati affidabili nella simulazione fluidodinamica.

Per avere una valutazione completa si sono effettuate simulazioni con fluido non

viscoso che consentano una formulazione Euleriana del problema e una simulazione

con fluido viscoso, dove si sono utilizzate le equazioni di Navier-Stokes ed il

modello di turbolenza Spalart-Allmaras (SA).

In entrambi i casi l‟oggetto dello studio è stato il profilo NACA0012 con geometria

bidimensionale. Una volta ricavati i dati delle simulazioni essi sono stati confrontati

con risultati noti e attendibili derivanti da altri codici di simulazione. Il confronto non

è stato solo numerico ma si è valutata anche la fisica in gioco.

Dai risultati ottenuti, non è possibile verificare quale codice sia migliore rispetto ad

altri con lo stesso oggetto di studio, ma si è potuta constatare l‟attendibilità dei

risultati ottenuti da SU2.

7

Il software presenta diversi moduli che svolgono compiti diversi, come ad esempio

risoluzione delle equazioni, deformazione della mesh e scomposizione del dominio

di calcolo. Sono poi presenti programmi scritti in linguaggio Python, che combinano

le funzioni dei moduli sopra citati per eseguire lavori più complessi come

l‟ottimizzazione geometrica.

Nel Capitolo 1 viene proposto l‟inquadramento matematico delle equazioni del

flusso, con particolare attenzione alle equazioni turbolente e ai modelli di turbolenza.

Nel Capitolo 2 viene introdotta la struttura del codice utilizzato e la numerica

implementata nelle diverse classi del programma.

Nel Capitolo 3 si propongono simulazioni con fluido viscoso e fluido non viscoso,

inoltre si analizzano i risultati numerici dei differenti casi e si confrontano con valori

di riferimento.

8

CAPITOLO 1

1 Le equazioni della fluidodinamica

Vengono di seguito introdotte le equazioni di Navier-Stokes complete per un fluido

newtoniano. Verrà poi esposta la formulazione delle equazioni mediate alla Favre-

Reynolds, focalizzando l‟attenzione sulla chiusura di queste ultime mediante i

modelli di turbolenza. Dopo una breve introduzione sulle varie tipologie di modelli

verrà presentato con maggiore dettaglio il modello Spalart-Allmaras, disponibile in

SU2.

1.1 Descrizione matematica del flusso

La fluidodinamica è lo studio del movimento di un numero di atomi o molecole; si

suppone quindi che la densità del fluido sia elevata abbastanza da poterlo considerare

un continuo. Assunta valida l‟ipotesi di continuo, le equazioni che regolano la

dinamica dei fluidi derivano direttamente dalle seguenti leggi di conservazione[1]:

conservazione della massa,

conservazione della quantità di moto,

conservazione dell‟energia.

La conservazione di una qualsiasi quantità del fluido implica che la variazione totale

della stessa, all‟interno di un volume di controllo arbitrario, sia la somma del flusso

netto della quantità attraverso il bordo del volume, delle forze interne e delle

eventuali sorgenti e forze esterne agenti sul volume di controllo.

Tale volume di controllo Ω è definito come la regione arbitraria di fluido delimitata

dalla superficie chiusa Ω e fissa nello spazio (punto di vista euleriano). É possibile

definire anche la superficie infinitesima dS e il versore associato normale n, come

mostrato nella Figura seguente:

9

Figura 1.1

Conservazione della massa

La legge di conservazione della massa stabilisce che non si può generare o

distruggere massa all‟interno del sistema fluido. Si consideri quindi il modello di

volume di controllo definito in Figura 1.1. In questo caso la quantità conservata è la

densità ρ. La variazione temporale della massa all‟interno del volume di controllo si

può esprimere come:

∫

Per la determinazione del flusso netto di massa che attraversa il volume si consideri

un elemento infinitesimo di superficie dS, e sia υ il vettore velocità del fluido e n il

versore normale alla superficie; il flusso di massa che attraversa la superficie dS sarà:

-ρ(υ n)

dove per convenzione il versore normale è diretto verso l‟esterno del volume di

controllo, quindi il flusso è entrante quando il prodotto scalare υ · n è negativo.

Questi sono gli unici termini presenti nell‟equazione di continuità poiché non sono

10

presenti sorgenti; per definizione è possibile scrivere:

∫ ∮ ( )

(1.1)

che rappresenta la forma integrale della legge di conservazione della massa, o

equazione di continuità.

Equazione della quantità di moto

È possibile ricavare l‟equazione partendo dal secondo principio della dinamica: la

variazione di quantità di moto è causata dalle forze agenti sull‟elemento di massa.

Sempre facendo riferimento alla definizione di volume di controllo (Figura 1.1), la

quantità di moto di un elemento infinitesimo del volume sarà:

la cui variazione nel tempo in tutto il volume è :

∫

La quantità conservata è la quantità di moto per unità di volume e il flusso può

essere espresso come :

∮ ( )

con la solita convenzione sulla normale e sul segno del flusso entrante. Le forze

agenti sul fluido sono di due tipi:

Forze di volume: sono forze che agiscono direttamente sulla massa, come la

forza di gravità, di Coriolis o la forza centrifuga e sono definite come forze

per unità di volume: ρfe;

11

Forze di superficie: agiscono sul bordo del volume e derivano dalla

distribuzione di pressione o dagli sforzi dovuti all‟attrito e sono forze per

unità di superficie.

Il contributo delle forze di volume è:

∫

Le forze di superficie sono dette sforzi e possono essere divise in due parti, una

componente isotropa di pressione e una componente di sforzi viscosi:

Fs= -pI + σ (1.2)

dove I è il tensore identità e σ il tensore degli sforzi viscosi.

È ora possibile mettere insieme tutti gli elementi per scrivere l‟equazione integrale

della conservazione della quantità di moto:

∫

∮ ( )

∮

∮( )

(1.3)

Equazione dell’energia

L‟equazione dell‟energia si ricava dal primo principio della termodinamica, il quale

stabilisce che la variazione nel tempo dell‟energia totale all‟interno del volume di

controllo è dovuta al lavoro delle forze agenti sul volume nell‟unità di tempo e al

flusso netto di calore all‟interno di esso.

L‟energia totale per unità di massa è definita come la somma dell‟energia interna

specifica (e) e dell‟energia cinetica specifica:

| |

(1.4)

12

La quantità conservata è l‟energia totale e, quindi la sua variazione nel volume si può

scrivere come:

∫

Il flusso attraverso il bordo di Ω è :

∮ ( )

con la solita convenzione per la direzione del versore normale.

A differenza delle due equazioni precedenti, nell‟equazione dell‟energia è presente

anche un flusso diffusivo, proporzionale al gradiente della quantità conservata;

siccome il flusso diffusivo è definito solo per fluidi a riposo, sarà proporzionale alla

sola energia interna:

dove

è il rapporto tra i coefficienti di calore specifico è il coefficiente di

diffusività termica. Questo termine, che rappresenta il flusso di calore dovuto al

gradiente termico, può essere espresso più semplicemente nella forma della legge di

Fourier :

con k coefficiente di conducibilità termica.

Il termine di sorgente di volume è composto dal lavoro per unità di tempo delle

forze di volume e dal flusso di calore nell‟unità di tempo dovuto all‟irraggiamento:

Rimane ora solo il contributo delle forze di superficie, che interviene come lavoro

nell‟unità di tempo delle stesse e può essere espresso come:

(1.5)

13

Unendo tutti i termini precedenti si giunge all‟equazione dell‟energia:

∫

∮ ( )

∮ ( )

∮ ( )

∮ ( )

∮( )

(1.6)

Solitamente l‟equazione (1.6) viene proposta in una forma leggermente differente.

A questo proposito si introduce l‟entalpia totale H e la relazione che lega questa

variabile termodinamica all‟energia totale e alla pressione:

| |

(1.7)

Raccogliendo insieme il flusso di energia e il termine di pressione, utilizzando la

(1.7) si può giungere alla forma più comune dell‟equazione di conservazione

dell‟energia:

∫

∮ ( )

∮ ( )

∮ ( )

∮ ( )

∮( )

(1.8)

Si è quindi giunti alla formulazione completa delle equazioni di governo del moto di

un fluido.

Prima di procedere è necessaria una riflessione su tali equazioni; infatti la forma

integrale, detta anche debole, delle leggi di conservazione ha due proprietà

fondamentali:

14

se non sono presenti sorgenti di volume, la variazione delle quantità

conservate dipende esclusivamente dai flussi sul bordo ∂Ω;

le equazioni sono valide anche in presenza di discontinuità nel flusso, ad

esempio onde d‟urto o discontinuità di contatto.

É evidente che, vista la sua generalità e le proprietà sopra descritte, questa sia la

forma più semplice da implementare in un codice fluidodinamico.

1.2 Equazioni di Eulero

Le equazioni di Eulero rappresentano una particolare forma semplificata delle

equazioni di Navier-Stokes, ottenute nel caso sussista l'ipotesi semplificativa di

flusso inviscido, ovvero flusso con viscosità trascurabile. Le equazioni di Eulero

sono riportate di seguito nella forma differenziale.

{

( )

( )

[ ( )]

Per quanto riguarda le prime due equazioni del sistema, esse descrivono il bilancio

della massa (equazione di continuità) e della quantità di moto. In queste equazioni

viene inoltre trascurata ogni forma di trasmissione del calore, per cui lo studio del

flusso inviscido è puramente fluidodinamico e non termofluidodinamico.

1.3 Sistema completo delle equazioni di Navier-Stokes

Nella sezione precedente sono state ricavate le equazioni di bilancio che governano

la dinamica dei fluidi, è ora possibile riscriverle tutte in un unico sistema di

equazioni, in una forma più generale e con l‟aggiunta delle relazioni termodinamiche

per la chiusura delle equazioni. Inoltre il sistema sarà proposto sia in forma integrale

15

che in forma differenziale. A tale proposito si introducono il tensore dei flussi

convettivi Fc, i tensori dei flussi viscosi divisi in contributo delle forze viscose Fv1 e

del trasferimento di calore Fv2, sia inoltre Q il vettore dei termini sorgenti; allora il

sistema completo delle equazioni di conservazione in forma integrale è:

∫ ∮ [ ( )]

∫

( )

dove U é il vettore delle variabili conservative: {

} ( )

mentre i flussi sono definiti come: {

} ( )

{

( )} {

( ) } ( )

Infine i termini sorgente sono: {

( )} ( )

L‟equazione rappresenta un sistema di cinque equazioni nelle cinque variabili

conservative, ma sono presenti sette incognite: ρ, v, E, p, T.

É quindi necessario aggiungere altre due equazioni, che dovranno essere relazioni

termodinamiche tra le variabili di stato, per chiudere il problema e inoltre fornire una

relazione per la viscosità dinamica μ e una per la conducibilità termica k.

Si ricorre quindi alla formulazione per i gas perfetti poiché, in aerodinamica classica,

è ragionevole ritenere di lavorare con fluidi che si comportano come se fossero gas

caloricamente perfetti.

16

Per tali gas la legge di stato è la seguente:

p = ρRT (1.14)

con R la costante specifica del gas. L‟entalpia e l‟energia interna specifiche sono

definite come:

(1.15)

Combinando l‟equazione (1.7) dell‟entalpia totale con (1.14) e (1.15) è possibile

scrivere la pressione come:

( ) ( | |

) (1.16)

mentre la temperatura si ricava dall‟equazione di stato. Per il coefficiente di viscosità

si utilizza la legge di Sutherland:

(

)

(1.17)

Dove . Mentre per

quanto riguarda il coefficiente di conducibilità termica si ritiene che sia costante in

tutto il fluido, tale ipotesi è sicuramente vera per i liquidi, ma è anche usata per l‟aria,

e vale:

(1.18)

dove Pr rappresenta il numero di Prandtl, che si considera anch‟esso costante nel

fluido e per l‟aria vale Pr = 0.72. Si ridefiniscono ora le due componenti dei flussi

viscosi raccogliendo la viscosità:

{

( ) } {

( ) } (1.19)

17

dove τ rappresenta i gradienti di velocità del tensore degli sforzi, ovvero:

Definendo:

il sistema delle equazioni di governo in forma integrale si riscrive come:

∫ ∮ [ (

)] ∫

(1.20)

Ora il bilancio equazioni-incognite risulta soddisfatto a patto di introdurre le corrette

condizioni al contorno e iniziali.

Applicando il teorema di Gauss è possibile riscrivere l‟integrale di superficie (il

secondo integrale dell‟equazione (1.20)) in integrale di volume e l‟equazione (1.20)

si può riscrivere in forma differenziale:

(

) (1.21)

Per semplicità successiva si riportano le equazioni di governo complete in notazione

tensoriale:

{

( )

( )

( )

( )

( )

( )

( )

(

)

18

1.4 Turbolenza

Le equazioni precedenti descrivono completamente il moto del flusso per qualsiasi

numero di Reynolds; nel caso di flussi turbolenti però una teoria analitica semplice

non esiste, la speranza è quindi di risolvere le equazioni sfruttando la potenza di

calcolo dei moderni computer per il calcolo di quantità interessanti dei flussi

turbolenti [3]. L‟affermazione precedente risulta più chiara riflettendo su alcune

proprietà dei flussi turbolenti. In un flusso turbolento, anche con condizioni al

contorno costanti, le grandezze scalari e vettoriali che caratterizzano il flusso non

sono stazionarie e presentano oscillazioni attorno ad un valor medio. Inoltre l‟analisi

di un flusso turbolento è complicata poiché occorre seguire l‟evoluzione di fenomeni

su scale di tempo e lunghezza molto differenti tra loro.

La scelta più comune per l‟analisi numerica di un flusso turbolento è quindi la

risoluzione delle equazioni tramite calcolatore; a questo scopo esistono tre differenti

approcci: DNS, LES, RANS.

La DNS (Direct Numerical Simulation) consiste nella risoluzione delle equazioni di

Navier-Stokes con una discretizzazione spaziale e temporale tale da risolvere tutte le

scale del moto. Tale approccio è concettualmente quello più semplice e fornisce un

livello di dettaglio massimo, non vi sono sostanziali perdite di informazione. Bisogna

però notare che il costo computazionale è estremamente elevato (cresce come Re3),

motivo per il quale, fino agli anni „70, tale approccio era inapplicabile a causa della

scarsa potenza di calcolo a disposizione.

Un approccio alternativo è dato dalla LES (Large Eddy Simulation) che prevede la

rappresentazione, senza approssimazioni, delle scale più grandi della turbolenza,

mentre le scale minori sono rappresentate tramite appositi modelli numerici. In

pratica le LES risolvono le equazioni del moto per un campo di velocità filtrato

rappresentativo delle scale maggiori della turbolenza.

Dal punto di vista del costo computazionale, le LES sono una via di mezzo tra la

DNS e le RANS; tenendo presente che il costo dipende dalle scale che si vogliono

rappresentare, idealmente, a parità di costo computazionale di una DNS, si riesce ad

19

ottenere lo stesso livello di dettaglio, cosa praticamente impossibile per le RANS

vista la differenza di approccio al problema, come sarà esposto in seguito.

Come discusso precedentemente, la maggior parte del costo computazionale di una

DNS è richiesto per la rappresentazione delle scale più piccole (o numeri d‟onda

grandi se si fa riferimento allo spazio di Fourier), mentre il maggior contenuto

energetico è racchiuso nelle scale più grandi del moto; nella LES la dinamica delle

scale maggiori è rappresentata esplicitamente, l‟influenza delle scale minori è

rappresentata da semplici modelli, si evita quindi il costo enorme della

rappresentazione delle scale minori della turbolenza.

In molte applicazioni ingegneristiche non è necessario il livello di dettaglio di DNS o

LES, molto più semplicemente, interessano proprietà macroscopiche, in generale

infatti basta la conoscenza dei valori medi delle quantità fluidodinamiche coinvolte.

É quindi conveniente in questi casi utilizzare la formulazione RANS (Reynolds

Averaged Navier-Stokes equations) che prevede la risoluzione delle equazioni

mediate.

In particolare si effettua dapprima la scomposizione alla Reynolds delle variabili

fluidodinamiche in una parte media ed una fluttuante, poi si procede con la

mediazione delle equazioni di Navier-Stokes. L‟operazione di media delle equazioni

produce un termine aggiuntivo detto tensore degli sforzi di Reynolds, che diventa

un‟incognita delle equazioni mediate, per cui il sistema risultante è sotto-

determinato.

Per pareggiare il conto equazioni-incognite è necessaria l‟introduzione dei cosiddetti

modelli di turbolenza che si dividono in modelli basati sull‟ipotesi di viscosità

turbolenta e altri sulla modellazione diretta del tensore degli sforzi di Reynolds.

I vantaggi di questo approccio sono sicuramente il ridotto costo della simulazione,

dovuto al fatto che si sta praticamente cancellando tutta la parte fluttuante del moto

del fluido; considerare solo i valori medi ha però dei contro.

Infatti gli svantaggi principali delle RANS sono dati dal fatto che è limitante

considerare solo il campo medio delle variabili fluidodinamiche (si pensi ad esempio

alle strutture coerenti che con questo approccio sono del tutto cancellate), inoltre tutti

i modelli di turbolenza sono caratterizzati da un certo empirismo, il quale impedisce

20

la determinazione a priori di modelli che forniscano previsioni ragionevoli per ogni

geometria e condizione di flusso.

Di seguito viene riportato il sistema completo delle equazioni Navier-Stokes mediate

alla Reynolds-Favre[4].

{

̅

( )

̃

( ̅ ̃)

( ̅ ̃ ̃)

̅

( ̃

̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅)

( ̅ ̃)

( ̅ ̃ ̃)

( ̌

̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅)

[ ̌ ( ̃ ̃ )]

(1.22)

É necessario fornire un‟interpretazione fisica dei seguenti termini dell‟equazione:

(

̌

) - diffusione molecolare di calore

(

̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ) - trasporto turbolento di calore

( ̃ ̃ ̌ ) - diffusione molecolare

( ̌

) - lavoro fatto dal tensore degli sforzi di Raynolds-Favre

21

1.5 Modelli di turbolenza

Le incognite del sistema possono essere determinate tramite i modelli di turbolenza,

che si suddividono in due principali tipologie: modelli basati sulla viscosità

turbolenta (Eddy-viscosity models) e modelli che si basano sulla scrittura di

equazioni differenziali che forniscono direttamente le componenti del tensore,

(Reynolds-Stress Models o RSM).

1.5.1 Modelli eddy-viscosity

Il più grande contributo alla modellazione turbolenta fu fornito da Boussinesq nel

1877; la sua idea si basa sull‟ipotesi di viscosità turbolenta derivante

dall‟osservazione del fenomeno fisico della turbolenza[4]. Tale ipotesi può essere

suddivisa in due parti; la prima è detta ipotesi intrinseca, afferma che la parte

anisotropa del tensore degli sforzi di Reynolds:

̅̅ ̅̅ ̅̅ ̅̅ ̅̅

è determinata dal gradiente di velocità media ̃

.

La seconda è detta ipotesi specifica e stabilisce la relazione tra ̃

(

̃

̃

)

̃

̃

(1.23)

Dove è la viscosità turbolenta; quest‟ultima è una quantità scalare positiva che

non rappresenta caratteristiche fisiche del fluido, ma è funzione delle condizioni

locali dello stesso. In generale l‟ipotesi intrinseca della viscosità turbolenta non ha

validità generale, ma per flussi semplici, in cui i gradienti di velocità media e le

caratteristiche turbolente evolvono lentamente, l‟ipotesi risulta essere ragionevole.

Nella pratica l‟ipotesi viene spesso utilizzata anche al di fuori dei suoi limiti di

validità, in questi casi l‟approccio RANS è ancora la scelta più valida quando si è

22

interessati a quantità macroscopiche del flusso dato il basso costo computazionale di

questa tecnica.

Bisogna ora ricavare un‟approssimazione del termine di trasporto turbolento di

calore; una soluzione comunemente utilizzata è la seguente:

( ̅̅ ̅̅ ̅̅ ̅̅ ̅)

̃

Dove è il coefficiente di conducibilità termica turbolento, definito come:

in analogia con il caso laminare. In questa relazione rappresenta il numero di

Prandtl turbolento, che in generale si ritiene costante nel fluido e per l‟aria vale

= 0.9.

Applicando l‟approccio appena descritto al sistema di equazioni differenziali il

tensore degli sforzi turbolenti è sostituito con una forma del tutto analoga al caso

laminare, ovvero:

̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̃

A questo punto l‟ultimo termine del membro di destra della seconda equazione della

(1.22) può essere riscritto come:

( ̃ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ )

( ̃)

(1.24)

dove il coefficiente di viscosità dinamica è dato dalla somma tra il coefficiente

laminare e quello turbolento:

Allo stesso modo il coefficiente di conducibilità termica k si riscrive come:

(

)

23

Data la completa somiglianza alle equazioni di Navier-Stokes non mediate è

possibile concludere che la loro versione mediata si può ottenere considerando le

incognite del problema come valori medi e sostituendo il coefficiente di viscosità e

quello di conducibilità termica con la somma della loro parte laminare più quella

turbolenta, non più uniforme in tutto il campo di moto, ottenendo il seguente sistema

di equazioni mediate:

(

)

(1.25)

dove in questo caso le variabili e

sono:

Questa tecnica, almeno dal punto di vista ingegneristico, è molto comoda poiché

richiede solo la determinazione della viscosità turbolenta. Quindi l‟ipotesi di

Boussinesq è la base di tutte le tecniche del prim‟ordine di chiusura delle equazioni

turbolente. Come già detto le limitazioni sono molteplici, ma i risultati, al di fuori del

campo di validità delle ipotesi, possono essere migliorati tramite correzioni o

utilizzando modelli di viscosità turbolenta non lineare, più complessi e che in questo

lavoro non verranno trattati.

I modelli eddy-viscosity si suddividono in tre categorie:

modelli algebrici,

modelli ad una equazione differenziale,

modelli a due equazioni differenziali.

Di seguito viene descritto un particolare modello ad una equazione differenziale, il

modello Spalart-Allmaras,

24

1.5.2 Modello ad una equazione di Spalart-Allmaras

Questo modello utilizza un‟equazione di trasporto per la viscosità turbolenta ̃

l‟equazione è stata sviluppata basandosi su considerazioni empiriche, analisi

dimensionale e invarianza galileiana. Nell‟equazione sono presenti tre funzioni

empiriche e ben otto coefficienti empirici, che sono stati tarati basandosi su risultati

per il mixing-layer 2D, scie e strato limite su lamina piana. Il modello risulta

accurato se utilizzato con flussi turbolenti con gradiente avverso ed è ottimizzato per

flussi turbolenti su ali e scie ad esse connesse. Ha particolari peculiarità anche dal

punto di vista numerico, infatti e un‟equazione locale che può essere facilmente

implementata su strutture multi-blocco o su griglie non strutturate. Inoltre risulta

essere un modello robusto, converge velocemente alla soluzione stazionaria e non

richiede elevata risoluzione nelle zone vicino parete.

La forma differenziale dell‟equazione, scritta in notazione tensoriale, è la seguente:

̅

( ̃ ̃)

( ) ̃ ̌

,

*( ̃)

̃

+ ̃

̃

-

*

+ (

̃ )

(1.26)

I termini di destra dell‟equazione rappresentano rispettivamente: produzione di

viscosità turbolenta, diffusione e dissipazione, distruzione di turbolenza.

̅̅̅̅

mentre il termine di produzione si ricava dalla seguente formula:

̃ ̌

(

)

25

( )

( )( )

̃

dove S è definita da:

√ √

( ̃

̃

)

Le altre funzioni empiriche dell‟equazione sono:

(

)

( )

̃

̃

Infine le costanti sono:

( )

Per poter scrivere le equazioni in un codice a volumi finiti è necessario ricondursi

alla forma integrale, che è:

∫ ̃ ∮ ( ) ∫

(1.27)

dove è il flusso convettivo ed è definito da:

̃ ̃

solitamente discretizzato tramite schema upwind del primo ordine, mentre il

26

flusso viscoso è:

dove è il vettore degli sforzi viscosi normali:

( ̃)

̃

Infine il termine sorgente è

( ) ̃ ̌ ̃

̃

[

] ( ̃

)

Sfruttando il teorema di Gauss è possibile riscrivere il secondo integrale

dell‟equazione (1.27) in integrale di volume e ottenere l‟equazione differenziale:

̃

(1.28)

Dove rappresenta il flusso viscoso e convettivo, ovvero:

̃

̃ ̃ ̃

Nell‟equazione è presente la variabile d che rappresenta la distanza dalla parete più

vicina; tale variabile è comune nei modelli di turbolenza ed è soluzione

dell‟equazione di Eikonal:

| |

(1.29)

in questo caso la soluzione d è proprio la distanza dalla superficie . Solitamente si

impone come condizione iniziale ̃ = 3 ÷5 , la stessa condizione si impone

solitamente sul contorno di inflow; sul contorno di outflow la variabile turbolenta

27

viene semplicemente estrapolata dal dominio. Sulle pareti solide invece si impone ̃

nulla e di conseguenza: = 0.

1.5.3 Modelli RSM

Nei modelli di turbolenza RSM (Reynolds-Stress Models)[1], introdotti per la prima

volta da Launder nel 1975, vengono risolte delle equazioni di trasporto modellate per

ogni componente del tensore degli sforzi e solitamente si aggiunge anche

un‟equazione aggiuntiva per la dissipazione . Di conseguenza, l‟ipotesi di viscosità

turbolenta non è più necessaria, eliminando uno dei difetti più grandi dei modelli

visti in precedenza. Per contro, un aspetto negativo è sicuramente l‟elevato costo

computazionale legato al fatto di dover risolvere diverse equazioni differenziali

aggiuntive. L‟equazione del trasporto del tensore degli sforzi di Reynolds si ottiene

dalle equazioni di Navier-Stokes ed è:

̅

In questa equazione il termine di convezione del flusso medio e il tensore di

produzione, , sono in forma chiusa; devono invece essere forniti modelli per il

tensore di dissipazione, , per il tensore legato ai processi distributivi, e per il

flusso degli sforzi di Reynolds,.

28

CAPITOLO 2

2 IL CODICE SU2

In questo capitolo si esporranno la struttura e i moduli principali del codice Stanford

University Unstructured (SU2)[2]. Il codice è stato sviluppato dal Dipartimento di

Aeronautica e Astronautica dell‟Università di Stanford e rappresenta un programma

innovativo nell‟ambito della CFD. Infatti è un software open-source scritto in C++

per la risoluzione di problemi descritti da equazioni alle derivate parziali (PDE),

inoltre presenta diversi strumenti (scritti in linguaggio Python) che permettono di

eseguire, tra le altre cose, ottimizzazione geometrica o calcolo parallelo. Nonostante

il programma sia stato sviluppato per applicazioni di CFD e ottimizzazione

geometrica, può essere esteso a qualsiasi problema con differenti equazioni di

governo. Il fatto stesso di essere stato concepito come un codice versatile e

facilmente generalizzabile, ha obbligato gli sviluppatori a progettare il codice in

moduli separati, che rappresentano modelli fisici e procedure differenti. Sempre

basandosi su questa filosofia è stato possibile sviluppare sia metodi a volumi finiti

(FVM) che ad elementi finiti (FEM), sistemi aggiunti corrispondenti e solutori multi-

fisici capaci di combinare entrambi gli approcci. L‟utilizzo di classi permette di

identificare chiaramente le parti di codice, in modo da poter facilmente modificarle o

migliorarle lasciando inalterate le altre funzioni del software. L‟ultima versione

disponibile sul sito http://su2.stanford.edu è la 3.2[11], ma vengono frequentemente

rilasciate versioni beta per sviluppatori. Il programma include numerosi moduli per

la risoluzione di sistemi di PDE su griglie non strutturate, che sono, tra gli altri:

SU2 CFD: il modulo principale per la risoluzione delle equazioni.

SU2 DDC: il modulo che suddivide il dominio per il calcolo in parallelo.

SU2 GPC: permette il calcolo della sensitività aerodinamica sulla superficie

da ottimizzare proiettando la sensitività su un set ridotto di variabili di design.

SU2 MDC: è il modulo che deforma la mesh a seguito di una variazione di

geometria.

29

2.1 STRUTTURA E CLASSI DEL CODICE

Nella sezione successiva sarà descritta l‟architettura del software e particolare

attenzione verrà posta nella descrizione della struttura delle classi ad oggetti. Si farà

riferimento al modulo SU2 CFD, che è il modulo principale del programma, ed

quello con cui sono state effettuate tutte le simulazioni di seguito riportate , ma molte

di queste classi sono utilizzate anche negli altri moduli.

CConfig: legge dal file di configurazione (estensione .cfg) la descrizione

fisica, matematica e numerica del caso da simulare.

COutput: scrive l‟output della simulazione in un formato che viene

specificato dall‟utente all‟interno del file di configurazione (Paraview,

Tecplot, o .csv).

CIntegration: risolve le equazioni di governo richiamando la classe figlia

CMultiGridIntegration. Questa classe è richiamata sia se la simulazione

viene effettuata su griglia singola che multipla e connette le seguenti

sottoclassi: CGeometry, CSolution e CNumerics per effettuare l‟integrazione

in spazio e tempo.

Di seguito vengono descritte brevemente le sotto classi del codice.

Classe CGeometry

Questa classe legge e processa il file della mesh (in formato nativo .su2) e include

diverse classi tra cui:

CPhysicalGeometry: costruisce la mesh duale, da quella in input, sulla quale

si basa la formulazione a volumi finiti.

CMultiGridGeometry: crea griglie consecutive sempre meno fitte a partire da

quella originale, utilizzate per il multi-grid.

CPrimalGrid e CDualGrid: sono le classi che definiscono le caratteristiche

geometriche della mesh primaria e duale.

30

Classe CSolution

In questa classe sono definite le procedure per il calcolo della soluzione; ogni classe

figlia rappresenta un solutore per le diverse equazioni di governo. Ad esempio

CEulerSolution è il solutore delle equazioni di Eulero e CAdjEulerSolution è il

solutore delle equazioni aggiunte di Eulero.

Queste sottoclassi ne richiamano diverse altre in CNumerics per discretizzare i vari

termini delle equazioni di governo. In generale le sottoclassi di CSolution hanno

come classi instanziate:

CVariable: è utilizzata per raccogliere le variabili, che dipendono dal tipo di

equazioni, in ogni punto della griglia. Una classe figlia può essere ad esempio

CNSVariable che raccoglie le variabili del sistema di equazioni di Navier-

Stokes.

CSparseMatrix: raccoglie i Jacobiani dei flussi e i termini sorgente in matrici

sparse per il calcolo implicito.

Classe CNumerics

È la classe che discretizza i termini delle equazioni di governo utilizzando gli schemi

numerici specificati nel file di configurazione. Ad esempio, in un calcolo implicito,

questa classe calcola i residui e i Jacobiani in ogni nodo della mesh, utilizzando le

variabili presenti in CVariables, passa i risultati a CSolution che a sua volta chiama

le funzioni in CSparseMatrix per la risoluzione del sistema lineare.

Il codice è stato concepito in modo da rimanere il più generale possibile e può

approcciare senza particolari problemi di programmazione a casi multidisciplinari; di

conseguenza, le equazioni di Navier-Stokes e l‟equazione della turbolenza vengono

risolte separatamente. Per questo motivo sono necessarie due classi instanziate in

CSolution (CNS- Solution e CTurbSolution), così come due classi instanziate in

CVariable (CNSVariable e CTurbVariable).

31

2.2 GRIGLIA DI CALCOLO

Le griglie di calcolo o mesh native, hanno estensione .su2 e sono in formato ASCII.

Poiché il codice è scritto per lavorare con griglie non strutturate, sono richieste

informazioni sui nodi e sulla connettività degli stessi, sul tipo di elemento

tridimensionale, bidimensionale o lineare e infine da quali nodi esso è formato. Nella

prima riga del file che definisce la mesh viene dichiarata la dimensione del problema,

infatti SU2 può operare sia su mesh bidimensionali che tridimensionali. Nella parte

successiva del file sono descritti gli elementi interni nel modo seguente:

%

%Inner element connectivity

% NELEM = 48970

9 0 1 3 2 0

9 1 4 5 3 1

9 4 6 7 5 2

dove ogni riga rappresenta la descrizione di un elemento interno. La prima colonna

definisce l‟ identificatore del tipo di elemento che viene descritto nella riga (secondo

il formato VTK). Seguono poi, i nodi che compongono l‟elemento descritto e il

numero che identifica l‟elemento costruito.

32

Una volta che la connettività è stata descritta, bisogna fornire le coordinate di ogni

nodo come esposto di seguito:

%

%Node coordinates

%

NPOIN = 49434

− 1.4008731282782294e − 01 − 1.6638704286217959e – 01 0

− 1.2947872467754354e – 01 − 1.5578240873270222e − 01 1

− 1.4164312923291425e – 01 − 1.6478390506414181e – 01 2

− 1.3097929203469980e – 01 − 1.5424567685235527e – 01 3

− 1.1887380719309064e – 01 − 1.4518144390092497e − 01 4

− 1.2031914461897376e – 01 − 1.4371109496123094e − 01 5

− 1.0839453032401439e – 01 − 1.3470607285589833e − 01 6

− 1.0978529214927452e – 01 − 1.3330132039213211e – 01 7

Dopo la dichiarazione del numero di nodi attraverso la definizione della variabili

NPOIN, si passa alla definizione delle coordinate e ad ogni punto viene assegnato il

numero che lo identifica. L‟ultima parte del file è costituita dalla definizione degli

elementi dei contorni del dominio di calcolo:

NMARK = 2

MARKER TAG = AIRFOIL

MARKER ELEMS = 277

3 138 139

3 139 209

3 209 279

3 279 349

3 349 419

3 419 489

Si dichiara prima il numero di contorni nel dominio, che prendono il nome di

Markers; per ogni Marker è necessario dichiarare il nome e il numero di elementi da

33

cui è costituito. Ancora una volta gli elementi di ogni Marker vengono dichiarati

mediante un identificatore che dichiara il tipo di elemento e i punti da cui è costituito

l‟elemento di contorno.

34

CAPITOLO 3

3 SIMULAZIONI

In questa sezione verranno esposti i risultati ottenuti da due campagne di

simulazione una per flusso non viscoso e una per flusso viscoso. Dalle prove

effettuate si ricava il valore dei coefficienti di portanza CL e di resistenza CD e si

valuta l‟andamento del residuo dell‟equazione di continuità in funzione del numero

di iterazioni per giudicare se la soluzione stazionaria sia stata, o meno, raggiunta.

Il residuo deve essere di almeno di 10-6

per considerare la soluzione a convergenza e,

dunque, effettivamente stazionaria. Questi risultati saranno poi confrontati con

risultati reperibili in letteratura per effettuare una validazione del codice stesso.

Tutte le simulazioni sono state effettuate usando la versione 3.2 del codice SU2 ed

un elaboratore dotato di processore Intel (R) core(TM) i3 M 350 2.27 GHz con 4

gigabyte di memoria RAM e sistema operativo a 64 bit. I file di configurazione

sono stati presi dai casi test diffusi da SU2 e sono specifici per questo tipo di

simulazione[10].

3.1 SIMULAZIONI CON FLUIDO NON VISCOSO

L‟oggetto di questo primo studio è stato quello di effettuare una simulazione

bidimensionale su un profilo NACA0012 (figura 3.1) e di confrontare i risultati

ottenuti con dati resi disponibili da Hegedus Aerodynamics[9] derivanti da analoghe

simulazioni.

Nella seguente campagna di simulazione si è ipotizzato un fluido non viscoso, perciò

si è utilizzata la formulazione Euleriana delle equazioni di conservazione e si è

utilizzata una sola griglia di calcolo la cui frontiera esterna, come si evince dalla

figura 3.2, ha forma circolare con raggio pari a 20 corde; la lunghezza della corda del

profilo è stata posta pari alla lunghezza unitaria.

35

La mesh è composta da 10216 elementi triangolari e 5233 punti. Sul bordo d‟uscita

gli elementi sono di dimensioni molto minori rispetto al resto del campo di moto per

permettere una migliore accuratezza della simulazione, data la criticità della zona. La

figura 3.2 mostra il reticolo di calcolo in prossimità del profilo, mentre la Figura 3.3

mostra il dettaglio della griglia in prossimità del bordo di uscita.

Figura 3.1: profilo alare NACA0012

Figura 3.2 : Dimensione del domino di calcolo

Dorso Bordo d’uscita Bordo d’attacco

Ventre

36

Figura 3.3 : Profilo NACA0012 Figura 3.4 : Bordo d‟uscita

Lo studio si compone di quattro simulazioni effettuate facendo variare l‟angolo

d‟attacco e il numero di Mach del flusso all‟infinito. L‟obiettivo della simulazione è

quello di osservare e confrontare i risultati derivanti dalla simulazione con quelli resi

disponibili da Hegedus Aerodynamics[9].

Le caratteristiche del flusso all‟infinito sono le seguenti:

Pressione: 101325 Pa

Temperatura: 273.15 K

SIMULAZIONE CON ANGOLO D’ATTACCO DI 0 E MACH 0.5

Questa prima simulazione è caratterizzata da un angolo d‟attacco nullo e flusso

ovunque subsonico; i valori ricavati del coefficiente di portanza e coefficiente di

resistenza sono i seguenti:

Coefficiente di portanza CL: 0.00056

Coefficiente di resistenza CD: 0.000113

37

Le Figure 3.5 e 3.6 riportano le isolinee del numero di Mach e della pressione

relative alla simulazione effettuata per M = 0.5, AOA = 0ᵒ.

Figura 3.5: Distribuzione del Mach; M = 0.5, AOA = 0ᵒ.

Figura 3.6: Distribuzione della pressione; M = 0.5, AOA = 0ᵒ.

L‟angolo d‟attacco nullo porta ad una simmetria delle distribuzioni del modulo della

velocità e della pressione tra la parte superiore e quella inferiore del profilo,

conseguentemente, la portanza dovrebbe avere valore esattamente nullo.

La presenza di errori numerici fa sì che il coefficiente di portanza sia piccolo, ma non

esattamente nullo. Inoltre, il flusso in regime subsonico fa sì che non si generi

nessun‟onda d‟urto su tutto il profilo.

38

La soluzione convergente, converge verso lo stato stazionario in quanto l‟andamento

del residuo dell‟equazione di continuità mostrato in figura 3.7 diminuisce

progressivamente fino a -8 (in scala logaritmica).

Figura 3.7: Andamento del residuo [RES_0] in funzione del numero d‟iterazioni.

Per valutare l‟attendibilità del risultato si è effettuato il confronto con la

distribuzione del coefficiente di pressione lungo il profilo, ottenuto da SU2 e il

medesimo reso disponibile da Hegedus Aerodynamics[figura3.8].

39

Figura 3.8: Confronto tra Cp

Da quanto emerge dalla figura 3.8, le due distribuzioni sono pressoché identiche se

non per la zona in prossimità del bordo di uscita; questo è dovuto alla chiusura del

profilo fornito da Hegedus Aerodynamics che è di poco superiore all‟unita mentre

per SU2 il profilo è unitario.

SIMULAZIONE CON ANGOLO D’ATTACCO DI 0ᵒ E MACH 0.8

La simulazione presenta lo stesso angolo d‟attacco delle precedente ma una velocità

del flusso maggiore, i risultati sono i seguenti:

Coefficiente di portanza CL : 0.00117

Coefficiente di resistenza CD : 0.00857

Dalle figure 3.8 e 3.9 relative alla distribuzione del numero di Mach e della

pressione, si nota la presenza di un‟onda d‟urto sul dorso e sul ventre del profilo

alare.

Anche se il flusso indisturbato è in regime subsonico, la deviazione che la corrente

subisce, dovuta alla presenza del profilo, porta ad un‟accelerazione del flusso che

raggiunge velocità supersoniche per poi decelerare e ricomprimersi attraverso

40

un‟onda d‟urto. L‟angolo d‟attacco nullo fa sì che, anche in questo secondo caso, il

campo di moto sia simmetrico rispetto alla corda del profilo alare.

Figura 3.9: Distribuzione del Mach; M = 0.8 , AOA = 0ᵒ.

Figura 3.10: Distribuzione della pressione; M = 0.8 , AOA = 0ᵒ.

L‟andamento del residuo mostrato in figura 3.11, conferma che la soluzione è

ottenuta è stazionaria.

41

Figura 3.11 : Andamento del residuo [RES_0] rispetto al numero d‟iterazioni; M =

0.8 , AOA = 0ᵒ.

Dal confronto fra le due distribuzioni di pressione[figura 3.12] si evince come queste

siano pressoché sovrapponibili. Come raffigurato in figura 3.12 l‟andamento del

coefficiente di pressione è decrescente dal bordo d‟attacco fino a circa il 50% della

corda del profilo alare, per poi aumentare repentinamente attraverso l‟onda d‟urto.

Figura 3.12 : Confronto tra Cp; M = 0.8 , AOA = 0ᵒ.

42

SIMULAZIONE CON ANGOLO D’ATTACCO DI 1.25ᵒ E MACH 0.5

Nel terzo caso considerato, l‟angolo d‟attacco maggiore di zero genera un valore di

portanza maggiore rispetto ai casi precedenti, inoltre si perde la simmetria che

caratterizzava le simulazioni ad angolo d‟attacco nullo. I risultati ottenuti dalla

presente simulazione sono riportati di seguito:

Coefficiente di portanza CL: 0.1749

Coefficiente di resistenza CD: 0.000214

Il flusso e subsonico su tutta la superficie del profilo perciò non si crea nessuna onda

d‟onda d‟urto [figure 3.13 e 3.14]. La soluzione è stazionaria perché i valori del

residuo sono prossimi allo zero macchina [figura 3.15].

Figura 3.13: Distribuzione del Mach; M = 0.5, AOA = 1.25.

Figura 3.14: Distribuzione della pressione; M = 0.5, AOA = 1.25.

43

Figura 3.15: Andamento del residuo [RES_0] rispetto al numero d‟iterazioni

Il confronto dei coefficienti di pressione è riportato in figura 3.16, nella quale si nota

l‟uguaglianza delle due distribuzioni, ad eccezione della zona terminale dove vi è

una lieve differenza che presumibilmente è dovuta alla diversa chiusura del profilo

alare, con cui Hegedus Aerodynamics ha effettuato le simulazioni. Rispetto ai casi

precedenti, nel grafico 3.16 ci sono due curve che sono dovute alla non simmetria

della distribuzione della pressione.

Figura 3.16: Confronto tra Cp; M = 0.5, AOA = 1.25.

44

SIMULAZIONE CON ANGOLO D’ATTACCO DI 1.25ᵒ E MACH 0.8

La simulazione ha fornito i risultati seguenti:

Coefficiente di portanza CL : 0.3269

Coefficiente di resistenza CD : 0.02134

A differenza del caso precedente, relativo ad un angolo di attacco α=0, questa

configurazione di flusso, crea due onde d‟urto di diversa entità. Come si vede dalle

figure 3.17 e 3.18 nella parte inferiore del profilo (intradosso) si ha un onda d‟urto

meno intensa rispetto a quella che si crea sull‟estradosso del profilo.

Figura 3.17:Distribuzione del Mach; M = 0.8, AOA = 1.25ᵒ.

Figura 3.18: Distribuzione della pressione; M = 0.8, AOA = 1.25ᵒ.

45

Anche in questa simulazione si è giunti a valori del residuo [figura3.19] che sono

prossimi allo zero macchina. Dunque, la soluzione si può considerare stazionaria.

Figura: 3.19: Andamento del residuo [RES_0] rispetto al numero d‟iterazioni;

M = 0.8, AOA = 1.25.

Anche in questa configurazione di flusso si osserva un buon accordo tra le

distribuzioni del coefficiente di pressione Cp ottenute nelle due soluzioni messe a

confronto, se non per le zone interessate dal salto di pressione. Si nota, infatti, come

per Hegedus il salto di pressione è netto, mentre per SU2 la distribuzione dei valori

in prossimità dell‟urto è graduale, soprattutto nel ventre del profilo. Questo è

presumibilmente dovuto al minor numero di punti con cui è stata discretizzata la

corda del profilo. Infatti, Hegedus Aerodinamics presenta 4097 punti, mentre il

profilo usato per la simulazione con SU2 presenta solo 198 punti.

46

Figura 3.20: Confronto tra Cp; M = 0.8, AOA = 1.25ᵒ.

3.1.1 TEMPI DI CALCOLO

Nelle simulazioni con fluido non viscoso il tempo necessario ad effettuare

un‟iterazione è di circa un secondo, per un numero totale d‟iterazioni che va dalle

230 alle 400 iterazioni. Pertanto il tempo necessario per creare una soluzione è di

pochi minuti. I ridotti tempi di calcolo sono dovuti alla natura del problema.

La formulazione Euleriana del problema, derivante dal fluido non viscoso, è la

condizione più semplice per effettuare una simulazione e anche quella che impegna

meno il codice dal punto di vista della potenza di calcolo; infatti è necessario un

numero ridotto di iterazione per arrivare a valori del residuo dell‟equazione di

continuità molto bassi.

47

3.2 SIMULAZIONI CON FLUSSO VISCOSO

Alla base di quest‟analisi vi è l‟obiettivo di valutare le effettive capacità del codice

SU2 nel generare risultati attendibili per flussi viscosi. Per la risoluzione del

problema, si utilizzano le equazioni di Navier-Stokes e il modello di turbolenza

Spalart-Allmaras (SA). Il profilo bidimensionale scelto è ancora il NACA0012.

Si è effettuata una campagna di prove di validazione che ci ha permesso di

confrontare i risultati ottenuti dal codice SU2, con quelli forniti dalla National

Aeronautics and Space Administration (NASA)[8] la quale, ha raccolto i risultati

ottenuti utilizzando sette codici CFD indipendenti. Essi sono: CFL3D (NASA

LARC, Stati Uniti d'America), FUN3D (NASA LaRC, Stati Uniti d'America), NTS

(NTS, Russia), JOE (Stanford, USA), SUMB (Stanford, USA), TURNS (Stanford /

Maryland, USA), e GGNS (Boeing, USA). I primi sei codici hanno usato la stessa

griglia di calcolo 897 x 257, mentre il codice GGNS ha utilizzato un reticolo di

calcolo di tipo adattativo che, durante la simulazione, modifica la grandezza locale

della griglia. Tutti e sette i codici hanno fornito risultati molto simili tra di loro. Le

piccole differenze sono presumibilmente dovute a errori di discretizzazione,

differenze di convergenza iterativa e differenze di condizioni al contorno.

Le simulazioni sono state eseguite in condizioni di flusso stazionario, a parità di

condizioni del flusso all‟infinito (cioè a grande distanza dall‟oggetto). Facendo

variare unicamente l'angolo d‟ attacco.

Le condizioni ambientali e di flusso utilizzate per le diverse simulazioni sono le

seguenti:

Numero di mach 0.15

Pressione 101325 Pa

Temperatura 273,15 K

Numero di Raynolds 6 000 000

Per effettuare il confronto, nelle simulazioni con SU2 si sono ricreate le stesse

condizioni ambientali e di flusso con cui NASA ha ottenuto i propri valori.

Utilizzeremo delle griglie di calcolo caratterizzate da diversa grandezza delle celle,

48

allo scopo di analizzare la variabilità dei risultati ottenuti in funzione della

risoluzione del reticolo di calcolo.

La campagna di prove si compone di tre simulazioni con tre angoli di attacco

differenti: 0 , 10, e 15 gradi. Per ogni angolo d‟attacco si sono effettuate due

simulazioni, una per ogni tipologia di mesh.

I quantitativi d‟interesse per il confronto sono i seguenti:

Coefficiente di portanza CL ad angolo di attacco (α) di 0, 10, e 15 gradi.

Coefficiente di resistenza aerodinamica CD ad angolo di attacco (α) di 0, 10,

e 15 gradi.

Coefficiente di pressione Cp vs x / c ad angolo di attacco (α) di 0, 10, e 15

gradi.

Coefficiente di attrito della superficie Cf vs x / c ad angolo di attacco (α) di

0, 10, e 15 gradi.

3.2.1 GRIGLE DI CALCOLO

Di seguito si analizzeranno nel dettaglio le griglie di calcolo con cui sono state

effettuate le simulazioni. Esse sono state prese direttamente dal caso test disponibile

sul sito ufficiale di SU2 [TestCases/rans/naca0012] e sono le mesh 449 x 129 e 225

x 65.

MESH 449 x 129

Questa mesh, mostrata in figura 3.21 è la più grande con cui si effettuano le

simulazioni. Essa presenta la classica forma a “C” ed è formata da 57824 punti e

57344 elementi quadrilateri con la frontiera esterna del dominio di calcolo posta a

508 corde dal profilo.

Il profilo alare, mostrato figura 3.22, ha il bordo d'uscita aguzzo e non tozzo, ed è

stato discretizzato con 256 segmenti [figura 3.23].

49

Figura 3.21: Domino di calcolo

Figura 3.22: Profilo NACA0012 Figura 3.23: Bordo d‟uscita

50

MESH 225 x 65

La mesh mostrata in figura 3.24 è formata da 14576 punti e 14336 celle quadrilatere

con la frontiera esterna posta 508 corde dal profilo [figura 3.25]. In questa griglia il

bordo d‟uscita è stato discretizzato con 128 segmenti [figura 3.26], ciò determina una

minore accuratezza di calcolo rispetto alla griglia precedentemente illustrata.

Figura 3.24: Domino di calcolo

Figura 3.25 : Profilo NACA0012 Figura 3.26 : Bordo d‟uscita

51

SIMUAZIONI CON ANGOLO D’ATTACCO DI 0ᵒ

Questa è la simulazione con angolo d‟attacco nullo; le condizioni di flusso, comuni a

tutte le prove fatte, sono state descritte in precedenza.

La tabella 3.1 riporta, il valore asintotico del residuo dell‟equazione di continuità e i

valori dei coefficienti di resistenza CD ottenuti con SU2 e delle simulazioni

effettuate dalla NASA.

Si è omesso il risultato del coefficiente di portanza perché con angolo d‟attacco nullo

anche il CL è nullo.

Inoltre, la NASA per le sue simulazioni non diffonde i valori del residuo

dell‟equazione di continuità ma, afferma che tutte le soluzione sono stazionarie.

Tabella 3.1 : Risultati simulazione; AOA = 0ᵒ.

I dati riportati in Tabella 3.1 sono riportati in forma grafica in Figura 3.30

Da questa simulazione emerge un primo dato rilevante giacché la simulazione che

presenta la migliore convergenza iterativa, cioè quella che più si avvicina alla

soluzione stazionaria è quella ottenuta con la mesh 449 x 129 infatti, il residuo

scende progressivamente fino a valori di – 6 in scala logaritmica come si vede

anche dalla figura 3.27. Nella simulazione con mesh 225 x 65 il valore del

SIMULAZIONI EFFETTUATE CON SU2 α = 0ᵒ

CODICI MESH CD RESIDUO (scala logaritmica)

SU2 449 x 129 0.00811 -6.33

SU2 225 x 65 0.00846 -5.80

RIFERIMENTI NASA α = 0ᵒ

CODICI MESH CD

CFL3D 897 x 257 0.00819 Soluzione stazionaria

FUN3D 897 x 257 0.00812 Soluzione stazionaria

NTS 897 x 257 0.00813 Soluzione stazionaria

JOE 897 x 257 0.00812 Soluzione stazionaria

SUMB 897 x 257 0.00813 Soluzione stazionaria

TURNS 897 x 257 0.00830 Soluzione stazionaria

GGNS 897 x 257 0.00817 Soluzione stazionaria

52

logaritmo del residuo non arriva a -6 (anche se per poco), ma comunque non vi è

una significativa differenza tra i risultati ottenuti con la mesh 449 x 129. Le figure

3.28 e 3.29 riportano le isolinee del numero di Mach e della pressione relative alla

simulazione effettuata con griglia 449 x 129. Si nota anche la simmetria che vi è tra

la parte superiore e quella inferiore della distribuzione di velocità e pressioni, questo

è dovuto al fatto che l‟angolo d‟attacco è nullo e di conseguenza la risultante delle

forze aerodinamiche in direzione ortogonale al flusso è nulla, ovvero il valore del

coefficiente di portanza è nullo.

Figura 3.27: andamento del residuo [RES 0] rispetto al numero d‟iterazioni;

AOA = 0ᵒ

Figura 3.28: Distribuzione del numero di Mach; AOA = 0ᵒ.

53

Figura 3.29 : Distribuzione della Pressione; AOA = 0ᵒ.

Dalla tabella 3.2 e dal grafico riportato in figura 3.30 è possibile osservare la

differenza di valori tra le simulazioni della NASA (nella tabella 3.2 sono riportati

solo i valori estremi) e quelle effettuate con SU2 con le due griglie di calcolo.

Nella figura 3.30 il valore 1 in ascissa rappresenta la griglia 897 x 257 con cui la

NASA ha ottenuto i propri risultati, mentre con il valore 2 in ascissa viene indicata la

griglia più fine con cui si è effettuata la simulazione con SU2, ovvero la griglia 449 x

129 e infine con il valore 4 in ascissa viene indicata la griglia 225 x 65.

Da quanto si evince in tabella 3.2 e dalla figura 3.30, con la griglia più fine la

differenza è modesta, mentre cresce con la griglia più grossolana (225 x 65).

La NASA afferma che per i propri valori di CD , lo scarto massimo passando dal

valore più piccolo a quello più grande, è del 4%. Quindi con la griglia 449 x 129

siamo nello scarto che la NASA certifica per le proprie simulazioni.

Tabella 3.2 : Differenze CD tra valori SU2 e NASA; AOA = 0ᵒ.

CD (codici NASA) CD ( codice SU2)

SCARTO %

MIN MAX MIN MAX

0.00812 0.00830 0.00811 (mesh 449 x129) 0.001 2.28

0.00846 (mesh 225 x 65) 1.92 10.41

54

Figura 3.30 : Valori del CD in funzione della grandezza della griglia.

Per avere una valutazione completa della capacità previsionale del codice SU2

confrontiamo non solo il valore dei coefficienti di portanza e resistenza ma, anche la

distribuzione delle forze aerodinamiche lungo il profilo.

Nella figure 3.31 viene riportata la distribuzione del coefficiente di pressione Cp

lungo il profilo. I valori per SU2 sono ottenuti con mesh 449 x 129 mentre, quelli di

riferimento NASA sono stati ottenuti dal codice CFL3D con la mesh 897 x 257.

Come vediamo dal grafico 3.31, la distribuzione è pressoché identica ossia i risultati

ottenuti con SU2 ricalcano perfettamente quelli ottenuti con CLF3D.

55

Figura 3.31: Distribuzione del coefficiente di pressione lungo il profilo; AOA = 0ᵒ.

In figura 3.32 viene mostrato il confronto tra l‟andamento del coefficiente di attrito

superficiale Cf ottenuto dal codice SU2 con mesh 449 x 129 e quelli resi noti dalla

NASA , ed ottenuti con codice CFL3D e mesh 897 x 257.

Come si vede in figura , la sovrapposizione dei due andamenti non è identica come

nel caso precedente. Il primo tratto è perfettamente uguale ma dal 30% della corda

verso valle vi è aumento del Cf ottenuto con SU2 rispetto a quello di CFL3D.

Sarebbe necessario poter effettuare una simulazione con SU2 e la medesima griglia

utilizzata dalla NASA per accertare se tale differenza sia imputabile alla risoluzione

del reticolo.

Figura 3.32 : Distribuzione del coefficiente d‟attrito lungo il profilo.

56

SIMULAZIONI CON ANGOLO D’ATTACCO DI 10ᵒ

In tabella 3.3 sono riportati i valori di CD e CL e sono presenti anche i valori ottenuti

dai sette codici utilizzati per il confronto. Anche i questo caso, la NASA non fornisce

i residui dell‟equazione di continuità, ma afferma che tutte le soluzioni sono

stazionarie. In questo caso il valore del coefficiente di portanza non è nullo perché

l‟angolo d‟attacco è diverso da zero, ovvero la risultante delle forze aerodinamiche

in direzione ortogonale al flusso non è nulla.

Tabella 3.3: Risultati simulazione; AOA = 10ᵒ.

In questa simulazione, come quella precedente, la griglia di calcolo più fine, è quella

che presenta la migliore convergenza iterativa, nonostante i valori del logaritmo del

residuo differiscono per pochi decimi rispetto alla griglia più grossolana.

La figura 3.33 mostra l‟andamento del numero di Mach mentre la figura 3.34 mostra

l‟andamento della pressione. Rispetto al caso precedente non vi è più la simmetria in

quanto l‟angolo d‟attacco è maggiore di zero quindi l‟intradosso del profilo è

caratterizzata dall‟alta pressione e bassa velocità, mentre l‟estradosso si caratterizza

SIMULAZIONI EFFETTUATE CON SU2 α = 10ᵒ

CODICI MESH CD CL RESIDUO (scala

logaritmica)

SU2 449 x 129 0.01363 1.1090 -6.208

SU2 225 x 65 0.01895 1.0865 -5.866

RIFERIMENTI NASA α = 10ᵒ

CODICI MESH CD CL

CFL3D 897 x 257 0.01231 1.0909 Soluzione stazionaria

FUN3D 897 x 257 0.01242 1.0983 Soluzione stazionaria

NTS 897 x 257 0.01243 1.0891 Soluzione stazionaria

JOE 897 x 257 0.01245 1.0918 Soluzione stazionaria

SUMB 897 x 257 0.01233 1.0904 Soluzione stazionaria

TURNS 897 x 257 0.01230 1.1000 Soluzione stazionaria

GGNS 897 x 257 0.01225 1.0941 Soluzione stazionaria

57

per una bassa pressione e l‟alta velocità del fluido. Nella figura 3.37 è mostrato

l‟andamento dei residui per la mesh 449 x 129.

Figura 3.35: Distribuzione del numero di Mach; AOA = 10ᵒ.

Figura 3.36 : Distribuzione della Pressione; AOA = 10ᵒ.

58

Figura 3.37: andamento del residuo [RES 0] rispetto al numero d‟iterazioni;

AOA = 10ᵒ.

La tabella 3.4 mostra le differenze relative al coefficiente di resistenza tra i valori

NASA e quelli ottenuti con il codice SU2. Invece la tabella 3.5 identifica la

differenza tra i coefficienti di portanza. Tutti i valori di scarto sono ottenuti passando

dal valore minore a quello maggiore.

Tabella 3.4 : Differenze CD tra valori SU2 e NASA; AOA = 10ᵒ.

Nelle figure 3.38 e 3.39 si vede come all‟aumentare delle dimensioni delle celle

della griglia di calcolo, i valori di CD e CL si allontanano dalla soluzione con la mesh

più fine, identificata nel grafico 3.38, con il valore 1 in ascissa.

La determinazione del coefficiente di resistenza è, in generale più complessa rispetto

a quella del CL , infatti le differenze con i valori di riferimento sono maggiori. Ciò

dipende dal fatto che la resistenza viscosa dipende dai gradienti della velocità.

CD (valori NASA) CD ( codice SU2)

SCARTO %

MIN MAX MIN MAX

0.01225 0.01245 0.01363 (mesh 449 x129) 9.4 11.2

0.01895 (mesh 225 x 65) 51.2 54.6

59

Figura 3.38: Andamento del CD rispetto alla grandezza della mesh.

Tabella 3.5: Differenze CL tra valori SU2 e NASA; AOA = 10ᵒ.

Figura 3.39 : Andamento del CL rispetto alla grandezza del mesh.

CL (valori NASA) CL ( codice SU2)

SCARTO %

MIN MAX MIN MAX

1.0891 1.1000 1.1090 (mesh 449 x129) 1.01 1.82

1.0865 (mesh 225 x 65) 0.2 1.24

60

All‟aumentare dell‟angolo d‟attacco i valori del coefficiente di resistenza presentano

uno scarto maggiore circa 10% per la soluzione con griglia 449 x 129, rispetto al

caso in cui l‟angolo d‟attacco era nullo.

Per quanto riguarda il coefficiente di portanza entrambe le soluzioni di SU2

presentano uno scarto intorno all‟ 1% .

Nelle figure 3.40 e 3.41 vengono mostrati l‟andamento del coefficiente di pressione

e quello di attrito superficiale lungo il profilo, ottenuti con il codice SU2 e griglia

449 x 129 e con il codice CFL3D con griglia 897 x 257. Come si vede dalle suddette

figure l‟andamento dei due coefficienti è pressoché identico.

Per quanto riguarda la figura 3.41, i valori di Cf generati dal codice CFL3D sono

disponibili solo per l‟estradosso del profilo alare. La figura dimostra in questo caso i

due andamenti coincidono perfettamente lungo nell‟estradosso. Per l‟intradosso

abbiamo solo i valori ottenuti dal codice SU2.

Figura 3.40 :

Figura 3.40: Distribuzione del coefficiente di pressione lungo il profilo; AOA = 10ᵒ.

61

Figura 3.41: Distribuzione del coefficiente d‟attrito lungo il profilo; AOA = 10ᵒ.

62

SIMULAZIONE CON ANGOLO D’ATTACCO DI 15ᵒ

Questa è la simulazione con angolo d‟attacco maggiore ed è quella che mette in

evidenza le criticità del codice. Difatti, all‟aumentare dell‟angolo d‟attacco diventa

sempre più difficile raggiungere la soluzione stazionaria. In tabella 3.5 sono riportati

i risultati derivanti da questa simulazione:

Tabella 3.5: Risultati simulazione; AOA = 15ᵒ.

Dalla simulazione emergono dati indicativi sulle differenti mesh e sull‟affidabilità

che le stesse hanno nel generare risultati attendibili. Dai dati riportati nella tabella 3.5

si evince che, con la diminuzione del numero di celle, i coefficienti aerodinamici si

allontanano dalla soluzione di riferimento. Anche se i valori dei residui sono in linea

con quelli ottenuti dalle simulazioni ad angolo d‟attacco minore, i valori dei

coefficienti CD e CL delle due griglie presentano differenze sostanziali.

Pertanto l‟unica soluzione stazionaria deriva dalla simulazione fatta con la mesh 449

x 129 mostrata in figura 3.42.

SIMULAZIONI EFFETTUATE CON SU2 α = 15ᵒ

CODICI MESH CD CL RESIDUO (scala

logaritmica)

SU2 449 x 129 0.02747 1.5384 -6.102

SU2 225 x 65 0.05848 1.3946 -5.758

RIFERIMENTI NASA α = 15ᵒ

CODICI MESH CD CL

CFL3D 897 x 257 0.02124 1.5461 Soluzione stazionaria

FUN3D 897 x 257 0.02159 1.5547 Soluzione stazionaria

NTS 897 x 257 0.02105 1.5461 Soluzione stazionaria

JOE 897 x 257 0.02148 1.5490 Soluzione stazionaria

SUMB 897 x 257 0.02141 1.5446 Soluzione stazionaria

TURNS 897 x 257 0.02140 1.5642 Soluzione stazionaria

GGNS 897 x 257 0.02073 1.5576 Soluzione stazionaria

63

Figura 3.42: andamento del residuo [RES 0] rispetto al numero d‟iterazioni;

AOA = 15ᵒ.

Infine in figura 3.43 per la griglia 449 x 129 e in figura 3.44 per la griglia 225 x 65

sono riportate le distribuzioni della velocità intorno al profilo alare. Si evidenzia la

differenza delle due distribuzioni difatti, sul dorso del profilo vi è un elevata

discordanza tra i due campi di moto, mentre per il ventre del profilo le due

distribuzioni sono similari.

Figura 3.43: Mach (mesh 449 x 129) Figura 3.45: Mach(mesh 225 x 65)

64

Di seguito [Figura 3.45 e 3.46] vengono riportate le distribuzioni delle pressioni

ottenute dalle rispettive griglie di calcolo.

Figura 3.45: Pressione (mesh 449 x 129) Figura 3.46: pressione (mesh 225x 65)

Quanto accennato in precedenza viene quantificato nelle tabelle 3.7 e 3.8.

La differenza tra i valori di riferimento NASA e quelli ottenuti con SU2 è ben più

elevata dei casi precedenti, addirittura lo scarto tra il valore ottenuto con i

coefficienti di resistenza generati dalla mesh 225 x 65 è del 282.1 % che fa risultare

questa griglia insufficiente per una simulazione che presenta un angolo d‟attacco di

questa entità.

Tabella 3.8: Differenze CD tra valori SU2 e NASA; AOA = 15ᵒ.

Le Figure 3.47 e 3.48 mostrano come variano i valori dei coefficienti in rispetto alla

grandezza della griglia di calcolo.

CD (valori NASA) CD ( codice SU2)

SCARTO %

MIN MAX MIN MAX

0.02073 0.02159 0.02747 (mesh 449 x129) 27.23 32.51

0.05848 (mesh 225 x 65) 270.8 282.1

65

Figura 3.47 : Andamento del CD rispetto alle dimensioni della mesh.

Il coefficiente di portanza ha ancora uno scarto accettabile rispetto hai valori riportati

dalla NASA infatti si attesta intorno al 1% nella mesh 449 x 129. Per la mesh

meno fitta, invece, il coefficiente il CL ottenuto presenta una differenza elevata

rispetto hai valori di riferimento NASA.

Tabella 3.9 : Differenze CL tra valori SU2 e NASA; AOA = 15ᵒ.

CL (valori NASA) CL ( codice SU2)

SCARTO %

MIN MAX MIN MAX

1.5446 1.5642 1.5384 (mesh 449 x129) 0.04 1.6

1.3946 (mesh 225 x 65) 10.7 12.16

66

Figura 3.48: Andamento del CL rispetto alla dimensione della mesh.

Le figure 3.49 e 3.50 riportano il confronto tra le distribuzione del coefficiente di

pressione e di attrito superficiale, generati dalla griglia più fine, con quello di

riferimento NASA ottenuto dal codice CFL3D. Anche in questo caso le due

distribuzioni sono quasi sovrapposte. Ciò conferma la validità del codice in quanto,

con una griglia di calcolo di minore risoluzione rispetto a quella usata dalla NASA è

possibile ottenere risultati soddisfacenti nonostante il valore elevato angolo

d‟attacco.

Figura 3.49 : Distribuzione del coefficiente di pressione lungo il profilo; AOA = 15ᵒ.

67

Figura 4.50: Distribuzione del coefficiente d‟attrito lungo il profilo; AOA = 15ᵒ.

3.2.2 TEMPI DI CALCOLO E NUMERO D’ITERAZIONI

Le simulazioni effettuate con le due mesh, oltre alle evidenti differenze nel livello di

convergenza iterativa raggiunto e valori numerici dei coefficienti, presentano

differenze sostanziali in termini di velocità di calcolo e numero d‟iterazioni da

compiere per raggiungere soluzione la soluzione stazionaria.

MESH 449 x 129

Essendo questa mesh la più grande con cui si è eseguita la simulazione, è quella che

presenta tempi di calcolo maggiori rispetto all‟altra, circa 4.2 secondi per ogni

iterazione ed un numero d‟iterazioni totali anch‟esso molto grande, ovvero compreso

tra le 12000 e le 13500 iterazioni per raggiungere la soluzione stazionaria. Pertanto,

i tempi totali di simulazioni sono molto elevati; circa 13 ore per ogni prova eseguita.

68

MESH 225 x 65

Questa mesh è di dimensioni ridotte rispetto alla precedente perciò permette una

notevole diminuzione dei tempi di calcolo. Infatti, per ogni iterazione l‟elaboratore

impiega 0.8 secondi e il numero d‟iterazioni che occorrono per generare una

soluzione è di circa la meta rispetto alla mesh più grande, quindi nell‟ordine delle

6000 - 7000 iterazioni. Questo comporta un risparmio di tempo cospicuo, difatti per

eseguire una simulazione, il tempo necessario è di circa 1,5 ore.

69

CAPITOLO 4

CONCLUSIONI

Nel presente lavoro di tesi è stata condotta una campagna di validazione di un codice

di calcolo “open-source” per la simulazione fluidodinamica sviluppato presso la

Università di Stanford.

A tale scopo, sono state condotte due diverse campagne di simulazione: una

nell‟ipotesi di flusso inviscido, la seconda nell‟ipotesi di flusso turbolento. É stato

così possibile, avere una visione chiara delle reali capacità del codice.

É stato esaminato il flusso attorno ad un profilo bidimensionale NACA0012 ed i

risultati ottenuti sono stati confrontati con altre soluzioni disponibili in letteratura.

Nella formulazione Euleriana del problema, non sono emerse discrepanze

significative fra i risultati ottenuti con il codice SU2 e le soluzioni di letteratura in

mio possesso.

Questo ci consente di affermare che il codice SU2 è in grado di generare risultati

attendibili quando si affrontano problemi di flusso non viscoso. Un‟affermazione di

questo tipo era prevedibile, poiché il tipo di problema non presenta una grande

complessità computazionale.

Il vero banco di prova è dato dalla simulazione con fluido viscoso perché come

spiegato nei capitoli precedenti la risoluzione del problema implica l‟utilizzo delle

equazioni di Navier-Stokes e del modello di turbolenza ad un equazione di Spalart –

Allmaras che ha una complessità maggiore.

La simulazione con flusso viscoso ci ha permesso di capire le reali capacità del

codice e soprattutto con quale griglia di calcolo si possono ottenere risultati

affidabili.

Le soluzioni ottenute con SU2 e mesh 449 x 129 mostrano che le differenze del

coefficiente di portanza e di resistenza rispetto ai risultati ottenuti dalla NASA, con

una griglia di calcolo che presenta circa il quadruplo delle celle che presenta la mesh

449 x 129 sono minime e tuttavia crescenti con il crescere dell‟angolo d‟attacco.

La risoluzione spaziale fornita dalla griglia 225 x 65 si è rivelata sostanzialmente

insufficiente, soprattutto all‟aumentare dell‟angolo di attacco.

70

Considerando le differenze tra griglie di calcolo e i risultati ottenuti possiamo

considerare il codice SU2 un valido strumento per la simulazione fluidodinamica in

quanto i risultati da esso generati sono globalmente verificati e in buon accordo con

quelli ottenuti da altri codici, sia che si operi con fluido viscoso e fluido non

viscoso.

Sarebbe opportuno poter riprodurre con SU2 i calcoli turbolenti utilizzando la

medesima mesh utilizzata dalla NASA, ma ciò richiede mezzi di calcolo che vanno

oltre un semplice laptop.

71

BIBLIOGRAFIA

[1] J. Blazek. Computational fluid dynamics: principles and applications. 2001.

[2] Francisco Palacios, Michael R. Colonno, Aniket C. Aranake, Alejandro Campos,

Sean R. Copeland, Thomas Economon, Amrita K. Lonkar, Trent W. Lukaczyk,

Thomas W. R. Taylor, and Juan J. Alonso. Stanford university unstructured (su2):

An open-source integrated computational environment for multi-physics simulation

and design. AIAA, 2013.

[3] Stephen B. Pope. Turbulent Flows. 2011.

[4] D.C. Wilcox. Turbulence Modeling for CFD. 1993.

[5] Anderson J.D.Jr. Computational fluid dynamics. The basics with applications.

1995.

[6] J.H. Ferziger, M. Peric Computational Methods for Fluid Dynamics 3rd Ed

(Springer 2002).

[7] David. W. Zingg Fundamentals of Computational Fluid Dynamics (1999).

SITOGRAFIA

[8] Implementing turbulence models into the compressible rans equations.

http://turbmodels.larc.nasa.gov/implementrans.html.

[9] http://www.hegedusaero.com/cases.html.

[10] https://github.com/su2code/TestCases.git.

[11] http://su2.stanford.edu/download.html