Simulazione canale turbolento1 Il problema e la sua adimensionalizzazionePer flusso in un canale si intende il flusso tra due pareti piane parallele edillimitate.

Il problema in esame è presentato in fig.1. Si indica con h la semialtezzadel canale (ovvero con 2h la distanza tra le due pareti). Con x, y e z sonoriportate, rispettivamente le direzioni longitudinale (cioè la direzione del motomedio), quella normale alle pareti e quella trasversale. Si simula soltanto untratto del canale, ovvero un parallelepipedo lungo Lx, alto 2h e largo Lz. Lxed Lz devono essere abbastanza grandi da non disturbare l’evoluzione delle piùgrandi strutture turbolente. Per questo si scelgono:

Lx = 2πhLz = πh

Figura 1: Dominio di calcolo.

Come condizioni al contorno si ipotizza che il flusso sia omogeneo sia indirezione x che in direzione z, ammettendo inoltre che le fluttuazioni sianoperiodiche: {

u′i (0, y, z, t) = u′i (Lx, y, z, t) ∀y, z, tu′i (x, y, 0, t) = u′i (x, y, Lz, t) ∀x, y, t

mentre sulle pareti si impone la condizione di aderenza, con ui = 0.

1

Bisogna quindi risolvere le equazioni di Navier-Stokes (DNS) o le equazio-ni filtrate della LES con queste condizioni nel parallelepipedo {0 ≤ x ≤ Lx} ×{−h ≤ y ≤ h} × {0 ≤ z ≤ Lz}:

∂ui∂xi

= 0

∂ui∂t

+∂

∂xj(uiuj) = −1

ρ

∂p

∂xi+ ν∇2ui

Se si considera la media di queste equazioni, tenendo conto che la velocitàmedia è solo in direzione x e che il flusso risulta omogeneo in x e z (le derivatein x e z delle grandezze mediate sono nulle), si ha:

∂U

∂t= −1

ρ

∂p

∂x+ ν

∂2

∂y2U − ∂

∂yu′v′

con U = U (y).Derivando rispetto ad x si ottiene:

∂2p

∂x2= 0 =⇒ ∂p

∂x= costante = −G =⇒ p = −Gx+ costante

con G modulo del gradiente di pressione medio.Applicando la decomposizione alla Reynolds, la pressione può essere scritta

come somma di una componente media (quella appena trovata) e una fluttuantep′:

p = Gx+ p′

Sostituendo nell’equazione vettoriale della quantità di moto, si ha:

∂ui∂t

+∂

∂xj(uiuj) = +

1

ρGδi1 −

1

ρ

∂p′

∂xi+ ν∇2ui

nella quale il primo termine del secondo membro è un termine costante indirezione x che agisce come una forzante.

In questo modo si impone il gradiente di pressione medio (ossia la differenzadi pressione media su un tratto fissato) e si calcola il flusso risultante.

2

Adimensionalizzazione• Per le lunghezze si usa la semialtezza del canale h.

• Per le velocità, poiché si impone il valore di G, conviene utilizzare unascala di velocità dedotta da G. Considerando che [G] = [ρ] [u]

2/ [L], ne

deriva che:

[u] =

([G] [L]

[ρ]

)1/2

=⇒ uτ =

(Gh

ρ

)1/2

con uτ velocità di attrito a parete (in quanto Gh è lo sforzo a parete τw).

Infatti, se si integra l’equazione per il moto medio tra −h ed h si ottiene:

0 = −G · 2h+ µ∂U

∂y

∣∣∣∣h−h

Ma lo sforzo medio a parete è proprio τw = µ∂U/∂y, quindi:

0 = −2Gh+ 2τw =⇒ τw = Gh

Intuitivamente, applicando l’equazione di equilibrio delle forze in direzionex si ottiene, fig. :

∆p · 2h · Lz = 2τw · Lx · Lz =⇒ τw =∆p

Lx· h = Gh

Figura 2: Bilancio delle forze.

3

Le grandezze adimensionali sono quindi:

xi =xih

ui =uiuτ≡ u+i

t =t

h/uτ=uτ t

h

p′ =p′

ρu2τ=

p′

Gh

Introducendole nelle equazioni di Navier-Stokes si ha:∂ui∂xi

= 0

∂ui

∂t+

∂

∂xj(uiuj) = +δ1i −

∂p′

∂xi+

1

Rτ∇2ui

con Rτ = uτh/ν numero di Reynolds definito dalle grandezze di parete.Il dominio, in coordinate adimensionali, è pari a: 0 ≤ x1 ≡ x ≤ 2π

−1 ≤ x2 ≡ y ≤ 10 ≤ x3 ≡ x ≤ π

4

2 Calcolo delle statistiche del flussoIl flusso è:

• mediamente stazionario;

• omogeneo in x ed in z;

• disomogeneo in y.

Tutte le grandezze mediate sono quindi funzione della sola variabile y (direzionenormale alla parete). Possono essere calcolate operando delle medie spaziali indirezione x e z e/o delle medie temporali:

f (y) =1

T

´ T0

[1

LxLz

´ Lz

0

´ Lx

0f (x, y, z, t) dxdz

]dt ∀f

Quindi, se si discretizza con Nx punti in x ed Nz punti in z, e si hanno adisposizione i risulitati in NT istanti temporali, si ottiene:

f (y) =1

NT

1

NxNz

∑i,k,s (xi, y, zk, ts)

i = 1, ..., Nxk = 1, ..., Nzs = 1, ..., Ns

Le statistiche dipendono dalla distanza dalla parete, quindi i piani ad essaparalleli sono statisticamente simili.

Inoltre, visto che il flusso è mediamente stazionario, tutti gli istanti sono allastessa densità di probabilità.

Spettri unidimensionaliCi sono due direzioni omogenee x e z in cui è possibile valutare gli spettro. Laprocedura che segue il codice utilizzato è la seguente:

1. le variabili interessate vengono antitrasformate;

2. si estraggono uno o più segmenti in una delle due direzioni (x o z);

3. si calcola lo spettro unidimensionale sul segmento |uk|2;

4. si media tra tutti i segmenti paralleli a pari distanza dalla parete;

5. si media anche tra più istanti temporali (in quanto il flusso è mediamentestazionario).

5

3 Metodo numericoLe direzione x e z hanno condizioni di periodicità, quindi è lecito applicare ilmetodo di Fourier-Galerkin.In direzione y invece si utilizza un metodo alle differenze finite, con una formulacentrata al IV ordine.

In questo modo si valuta ˆui (kx, y, kz, t), trasformata di Fourier della varia-bile ui (x, y, z, t) nelle direzioni x e z soltanto.Di seguito sono proposte le equazioni di bilancio della quantità di moto adimen-sionati in forma semidiscreta:

∂u∂t

+ ikxu2 +∂

∂yuv + ikz wu = −ikxp′ +

1

Rτ∇2u+ 1

∂v∂t

+ ikxuv +∂

∂yv2 + ikz vw = −∂p

′

∂y+

1

Rτ∇2v

∂ w∂t

+ ikx uw +∂

∂yvw + ikz w2 = −ikz p′ +

1

Rτ∇2 w

con ∇2 pari a:

∇2 = −(k2x + k2z

)+

∂2

∂y2

L’avanzamento nel tempo avviene mediante il metodo dei passi frazionari.Nel caso più semplice (Eulero), dato il valore u(n) al passo n, si calcola la velocitàprovvisoria u∗ come:

u∗ = u(n) + ∆t[C(n) +D(n)

]dove C e D costituiscono i termini convettivo e diffusivo.

Si aggiunge il gradiente di pressione:

u(n+1) = u∗ −∆t · ∇p

e si impone che la divergenza della velocità sia nulla, ossia ∇ · u(n+1) = ∇ ·u(n) = 0. è quindi possibile valutare la velocità al passo successivo, risolvendol’equazione che segue:

∇2p =1

∆t∇ · u∗

6

4 Variabili di pareteLe variabili di parete sono:

u+i =uiuτ

y+ =yPuτν

con yP distanza dalla parete.Le velocità calcolate nella simulazione sono già espresse in unità di parete,

a differenza delle coordinate y, adimensionalizzate con h.Se si fa riferimento alla parete inferiore si ha yP = y + h, da cui:

y+ =(y + h)uτ

ν=uτh

ν

(1 +

y

h

)= Rτ (1 + y)

Per la parete superiore, da yP = h− y, si ricava analogamente:

y+ =(h− y)uτ

ν=uτh

ν

(1− y

h

)= Rτ (1− y)

nelle quali si è indicato con y la coordinata y adimensionale del codice.

7

5 Approfondimento metodi DNS - LES

5.1 DNSLa Direct Numerical Simulation è la simulazione numerica in cui le equazioni diNavier-Stokes descritte in precedenza sono risolte numericamente.

Tutte le scale spaziali della turbolenza devono essere risolte nella grigliacomputazionale, dalle scale dissipative più piccole (scala di Kolmogorov, η), finoalla scala integrale L, associata ai moti contenenti la maggiorparte dell’energiacinetica.

Tale scala integrale deve essere contenuta nel dominio computazionale, talecioè che sia Nh = L e, al tempo stesso, la scala di Kolmogorov deve essere taleda avere h ≤ η, dove N è il numero di punti della griglia e h l’incremento diogni intervallo.

5.2 LESIl metodo Large Eddy Simulation si basa sulla teoria di Kolmogorov, il cheimplica la dipendenza delle grandi strutture turbolente dal contorno e il com-portamento universale delle strutture più piccole.

Si risolvono le equazioni fino ad una certa scala della cascata di Richardson,fino ai vortici di grandi dimensioni, tramite l’uso di un filtro spaziale.

Le piccole strutture, in quanto filtrate, vengono modellate con i modelli diturbolenza di sottogriglia.Le variabili filtrate (nello spazio) sono definite come:

〈f〉 (x, t) =

ˆf (x + η)G (η) dη

In particolare, per le componenti di velocità si ottiene:

〈u〉 (x, t) =

ˆu (x + η)G (η) dη

dove G è la funzione filtro, che può essere omogenea o non omogenea ed èsempre funzione dell’ampiezza δ del filtro.

In fig.3-4 sono proposti un filtro “a scatola” (box filter) e un filtro di tipoGaussiano. Il primo ha il vantaggio di avere un supporto compatto; il secondo,a fronte dell’assenza di discontinuità (quindi la trasformata di Fourier è sempredefinita) presenza lo svantaggio di mantenere valori non nulli, per la variabilefiltrata, anche a distanze infinite.

8

Figura 3: Filtro a scatola (box filter).

Figura 4: Filtro Gaussiano.

Applicando il filtro spaziale alle equazioni di Navier-Stokes, si ottengono leseguenti equazioni filtrate:

∂ 〈ui〉∂xi

= 0

∂

∂t〈ui〉+

∂

∂xj〈uiuj〉 = −1

ρ

∂

∂xi〈p〉+ ν∇2 〈ui〉+

∂

∂xjτSGSij = 0

con il termine che contiene gli sforzi di sottogriglia indicato con τSGSij (modellodi Smagorinsky):

τSGSij = 〈ui〉 〈uj〉 − 〈uiuj〉L’effetto che questi sforzi devono riprodurre è quello di sottrarre alla parte

risolta della corrente fluida quell’energia cinetica che i vortici più piccoli, inquanto filtrati, non sono in grado di dissipare.

La parte isotropa del tensore τij ha gli effetti di una pressione, la parteanisotropa viene legata al tensore velocità di deformazione Sij , tramite l’eddyviscosity νδ:

9

τAij = 2νδ 〈Sij〉

dove:〈Sij〉 =

1

2

[∂ 〈ui〉∂xj

+∂ 〈uj〉∂xi

]Uno dei modelli di sottogriglia più utilizzati è appunto quello di Smagorinsky,

nel quale si impone che l’eddy viscosity sia proporzionale alla scala del filtrospaziale δ, e alla scala di velocità delle fluttuazioni δ |Sij | al livello δ:

νδ = (Csδ)2 |Sij |

Il problema di questo modello è quello di determinare il coefficiente Cs, in quan-to l’eddy viscosity è funzione del campo di moto e non delle proprietà del fluido.

5.2.1 Problemi

Filtro non uniforme Il filtro descritto sopra non è uniforme, quindi si ha:

∂

∂xi〈f〉 6=

⟨∂f

∂xi

⟩Infatti, ad esempio per un filtro a scatola, nel caso unidimensionale si ha:

〈f〉δ (x) =1

δ

ˆ δ/2

−δ/2f (x+ η) dη =

1

δ

ˆ 1/2

−1/2f (x+ δξ) dξ

Dalla derivazione nello spazio si ricava:

d

dx〈f〉δ =

d

dx

ˆ 1/2

−1/2f (x+ δξ) dξ =

=

ˆ 1/2

−1/2

d

dx[f (x+ δξ)] dξ =

ˆ 1/2

−1/2f ′ (x+ δξ) (1 + δ′ξ) dξ =

=

ˆ 1/2

−1/2f ′ (x+ δξ) dξ︸ ︷︷ ︸ +δ′

ˆ 1/2

−1/2f ′ (x+ δξ) dξ︸ ︷︷ ︸

〈f ′〉 ∂

∂δ〈f〉

nella quale l’ultima uguaglianza deriva da:

∂

∂δ〈f〉δ =

∂

∂δ

ˆ 1/2

−1/2f (x+ δξ) dξ =

ˆ 1/2

−1/2f ′ (x+ δξ) ξdξ

10

In conclusione si può scrivere:⟨df

dx

⟩δ

=d

dx〈f〉δ −

dδ

dx

∂

∂δ〈f〉︸ ︷︷ ︸

errore di commutazione

Nel caso generale, senza l’ipotesi di unidimensionalità:⟨∂f

∂xi

⟩δ

=∂

∂xi〈f〉δ −

∂δ

∂xi

∂

∂δ〈f〉

Bisogna approssimare e discretizzare questi termini.Uno dei metodi utilizzabili è quello delle differenze finite. Nel caso di metodocentrato del primo ordine, si ha:

∂

∂δ〈f〉δ =

〈f〉δ+h − 〈f〉δ−h2h

+O(h2)

Considerando δ = h:

∂

∂δ〈f〉δ =

〈f〉2δ − 〈f〉02δ

+O(h2)

Il problema diventa così quello di ricavare i valori di 〈f〉2δ e 〈f〉0 = f a partiredalla 〈f〉δ.

Si valuta lo sviluppo in serie, nel caso unidimensionale di filtro a scatola:

f (x+ δξ) = f (x) + δξf ′ (x) +1

2δ2ξ2f ′′ (x) +

1

3!δ3ξ3f ′′′ (x) +O

(δ4)

e si ricava:

〈f〉δ = f (x)+δ

ˆ 1/2

−1/2ξdξf ′ (x)+

1

2δ2ˆ 1/2

−1/2ξ2dξf ′′ (x)+

1

3!δ3ˆ 1/2

−1/2ξ3dξf ′′′ (x)+O

(δ4)

=

= f (x) +1

24δ2f ′′ (x) +O

(δ4)

da cui:〈f〉0 = f = 〈f〉δ −

1

24δ2f ′′ +O

(δ4)

Analogamente si ottiene:

〈f〉2δ = f (x) +1

6δ2f ′′ (x) +O

(δ4)

11

Applicando il filtro di ampiezza 2δ alla funzione filtrata (ma con filtro di am-piezza δ):

〈〈f〉δ〉2δ = 〈f (x)〉2δ +1

24

⟨δ2f ′′ (x)

⟩2δ

+O(δ4)

= 〈f〉2δ +1

24δ2f ′′ +O

(δ4)

che permette la seguente:

〈f〉2δ = 〈〈f〉δ〉2δ −1

24δ2f ′′ +O

(δ4)

In questo modo si ottiene:

〈f〉2δ − 〈f〉0 = 〈〈f〉δ〉2δ − 〈f〉δ +O(δ4)

ed infine l’espressione della derivata:

∂

∂δ〈f〉 =

〈〈f〉δ〉2δ − 〈f〉δ2δ

+O(δ3)

da cui un errore di commutazione dell’ordine di O(δ4).

Sovrastima degli sforzi a parete del modetto di Smagorinsky Il valoredella eddy viscosity è troppo elevato in prossimità della parete, quindi si devecorreggere il modello.

La correzione di Van Driest consiste nel moltiplicare la νSGSδ per una funzionecorrettiva il cui valore dipende dalla distanza da parete, il cui andamento èpresentato in fig.5 e proposto di seguito:

fV D,1(y+)

= 1− e−y+

A+

con A+ ' 25.In generale la eddy viscosity si esprime come:

vδ = vSGSδ · fV D(y+)

In alternativa con la funzione correttiva può essere del tipo:

fV D,2(y+)

= 1− e−

y+

A+

2

con la quale si ha una transizione meno graduale da 0 ad 1 e si dimezza il bufferlayer.

12

Figura 5: Correzione di Van Driest.

Condizioni al contorno In prossimità della parete, quindi per y → 0, occorrediscriminare tra due casi, in funzione dell’ampiezza del filtro applicato δ.

1. Nel caso in cui δ → 0 allora si ottiene 〈u〉 → u; qui è possibile applicarele condizioni di aderenza alle equazioni di Navier-Stokes, ovviamente conuna griglia fitta.

2. Nel caso in cui δ 9 0 non è possibile applicare le condizioni di aderenza.

Si valuta il valore della coordinata di parete minima ymin, al di sotto della qualesi deve intervenire inserendo delle condizioni al contorno opportune per simulareil problema reale.

• Per ymin grande si inseriscono condizioni basate sul modello fisico delcomportamento del flusso dalla parete (y = 0) alla ymin; ad esempio nelmodello di Schumacher ymin è considerato appartenente alla log region1 esi ipotizza che ∂ 〈u〉 /∂y ∝ τ .

• Per ymin piccolo (vicino a parete, cioè 0 ≤ y+min ≤ 5) si inseriscono dellecondizioni basati sugli sviluppi in serie:

u (0) = u (ymin)− ymindu

dy(ymin) +

1

2y2min

d2u

dy2(ymin) + ... = 0

Con 〈u〉 = u + αδ2u′′ + O(δ4)e alcuni semplici passaggi algebrici, qui

omessi, si trova la condizione al contorno per 〈u〉:

〈u〉 = ymin∂

∂y〈u〉+

1

2y2min

∂2

∂y2〈u〉+

δ2

2∇2 〈u〉

1Si rimanda al capitolo successivo per ulteriori chiarimenti.

13

6 Codici e programmi

Credenziali di accessoI file necessari alla simulazione del canale turbolento si trovano sul server 130.192.25.149.

L’accesso al server in modalità SSH può essere effettuare con le seguentichiavi di accesso:

utente: esercitazioni_snft

password: davide

Quindi si digiti: ssh esercitazioni_snft@130.192.25.149.

Percorsi• ~/CanaleTurbolento/Codici → contiene i file sorgente del codice

• ~/CanaleTurbolento/DatiIniziali→ contiene i dati iniziali per le simulazioniDNS e LES.

Lista eseguibiliIl compilatore dei codici che si trovano nelle cartelle è sempre gfortran.Tutte le istruzioni sono riportate nel Makefile all’interno delle directory ~/Pro-cessing/ e ~/Post-Processing/ situate nella cartella ~/Codici/. Si hanno a dispo-sizione i file riportati di seguito.

Processing

• dnsles.exe → soluzione delle equazioni di Navier-Stokes; si genera digitan-do make dnsles.exe nella cartella ~/Codici/Processing/. Ha bisogno del fileparam.h reperibile nella relativa cartella ~/DatiIniziali/. Si consiglia di spo-stare l’eseguibile nella cartella ~/Simulazione/ unitamente al file param0correlato e quindi avviare la simulazione2.

Post-processing

• mediat180media.exe → calcola le componenti medie di velocità lungo ladirezione x; il file *.umedia generato si compone di quattro colonne (y, y+,U , Us)3;

2Il file param0 contiene i parametri fisici del problema, mentre il filse param.h contiene iparametri della discretizzazione.

3Con Us si indica la componente di velocità specchiata.

14

• mediat180up.exe → calcola le componenti fluttuanti della velocità e glisforzi di Reynolds; il file d’uscita *.up è diviso in otto colonne (y, y+, u′,v′, w′, u′v′, u′w′, v′w′), mentre il file *.upm riporta i valori delle stessevariabili, ma specchiate rispetto alla metà del canale;

• mediat180v.exe → calcola le componenti di vorticità; il file *.c3 generatopresenta i valori di y, y+ e della vorticità nelle tre direzioni ωx, ωy e ωz;per il file *.c3m vale quanto scritto sopra.

I codici sorgenti che generano gli eseguibili di cui sopra si trovano nella cartella~/Codici/Post-Processing/. Per la compilazione, occorre digitare ./compilazioneall’interno della cartella. I file media*.exe essi devono quindi essere spostatinella cartella contenente i risultati della simulazione (directory ~/Simulazione/),unitamente al file param0.

Visualizzazione Nella cartella ~/Codici/Visualizzazione/ si trovano i codicinecessari alla separazione nelle tre direzioni delle componenti di velocità mediee fluttuanti dal generico file *.txt risultato della simulazione al generico istantetemporale. Occorre digitare il comando ./compila all’interno della directory (ènecessaria la presenza del file param.h), quindi spostare l’eseguibile conv.x nellacartella ~/Simulazione/.

I file in uscita, in formato .bin possono essere convertiti con il programmabin2vtk e visualizzati tramite VisIt.

Parametri fisici I parametri fisici della simulazione possono essere scelti emodificati dall’utente nel file param0. Esso contiene:

Rτ numero di Reynolds d’attrito;

dt passo temporale;

ndt numero totale di passi;

car indica il suffisso con cui verranno salvati i file della simulazione;

nomedati nome del file con i dati iniziali (D010.txt per la DNS e z180.t0 perla LES).

Inoltre è possibile modificare le dimensioni del dominio in x, y e z, sceglierela posizione della frontiera ymin, e se utilizzare o meno la correzione di VanDriest e/o valutare l’errore di commutazione dato dalla disuniformità del filtroapplicato.

Per maggiori chiarimenti si rimanda alle note contenute nelle varie cartelle.

15