Appunti di Fluidodinamica Numerica

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