Università degli Studi della Basilicata

FACOLTA' DI INGEGNERIA

Corso di laurea specialistica in Ingegneria Meccanica

Tesi di laurea

in Gasdinamica

Modellazione teorica e numerica di gas ionizzati

RELATORI CANDIDATO

Prof. Aldo Bonfiglioli Raffaele Pepe Matr. 31747 – IM Dr. Antonio D’Angola

Anno Accademico 2010 – 2011

1

Indice

Introduzione ................................................................................................................................................. 2

1 Modello termochimico per i gas ad alta temperatura ................................................................ 3

1.1 Il plasma ................................................................................................................................................................................. 5

1.2 Modello termodinamico .................................................................................................................................................. 6

1.3 Modello chimico.................................................................................................................................................................. 8

1.4 Modello chimico per l’aria .............................................................................................................................................. 9

1.5 Fit delle proprietà termodinamiche ........................................................................................................................ 12

2 Le equazioni della fluidodinamica ................................................................................................. 17

2.1 Le equazioni di conservazione .................................................................................................................................. 18

2.2 Equazione di conservazione della massa .............................................................................................................. 20

2.3 Equazione di conservazione della quantità di moto ......................................................................................... 21

2.4 Equazione di conservazione dell’energia .............................................................................................................. 22

2.5 Equazioni di Navier-Stokes ......................................................................................................................................... 24

2.6 Equazioni di Eulero ........................................................................................................................................................ 26

3 Fluctuation splitting per gas reale ................................................................................................. 29

3.1 Introduzione...................................................................................................................................................................... 30

3.2 Approccio ai Volumi Finiti ........................................................................................................................................... 30

3.3 Approccio Fluctuation Splitting ................................................................................................................................ 33

3.4 Decomposizione caratteristica delle equazioni di Eulero 1D ....................................................................... 35

3.5 Linearizzazione conservativa .................................................................................................................................... 44

3.6 Generalizzazione di Vinokur-Montagnè ................................................................................................................ 56

3.7 Distribuzione del residuo ............................................................................................................................................ 61

4 Condizioni al contorno ....................................................................................................................... 65

4.1 Introduzione...................................................................................................................................................................... 66

2

4.2 Condizioni al contorno NSCBC ................................................................................................................................... 67

4.3 Metodo del flusso correttivo ...................................................................................................................................... 71

5 Fluctuation splitting 2D per gas reale ........................................................................................... 75

5.1 Introduzione...................................................................................................................................................................... 76

5.2 Decomposizione caratteristica delle equazioni di Eulero 2D per gas reale ............................................ 79

5.3 Linearizzazione conservativa delle equazioni di Eulero 2D per gas reale .............................................. 83

6 Risultati delle simulazioni ................................................................................................................ 92

6.1 Introduzione...................................................................................................................................................................... 93

6.2 Il Flusso di Rayleigh ....................................................................................................................................................... 93

6.3 Simulazioni numeriche per un flusso di Rayleigh 1D ...................................................................................... 99

6.4 Modello con sorgente variabile ...............................................................................................................................106

5.5 Simulazioni per il modello con sorgente variabile ..........................................................................................111

Conclusioni ............................................................................................................................................... 117

Bibliografia .............................................................................................................................................. 119

Appendice A: Proprietà termidinamiche dell’aria ..................................................................... 122

Appendice B: Trasformazioni Jacobiane ....................................................................................... 132

Appendice C: Codice numerico in Fortran ................................................................................... 139

3

Introduzione

Negli ultimi decenni lo studio teorico e numerico dei gas ionizzati ha avuto delle notevoli ricadute

nell’ambito delle applicazioni industriali ed aerospaziali delle tecnologie che coinvolgono aeriformi

che si trovino allo stato di plasma. A causa di ciò, si è reso necessario sviluppare modelli chimici e

termodinamici in grado di prevedere in maniera accurata le proprietà termodinamiche e di trasporto

dei gas lontani dalle condizioni per cui risulta valido il modello di gas perfetto. Nella maggior parte

delle applicazioni, e con buona approssimazione, è possibile considerare il plasma in condizioni di

equilibrio locale termodinamico e chimico: in questo modo non è necessario considerare le reazioni

chimiche in maniera esplicita.

La ionizzazione in un gas può essere ottenuta mediante l’aumento di temperatura dovuto ad una

sorgente di calore o grazie ad una scarica elettrica, per cui parte del presente lavoro di tesi è stato

dedicato allo studio del flusso unidimensionale, non viscoso, con sorgente interna di calore, noto in

letteratura come flusso di Rayleigh. A tale scopo è stato modificato un codice di simulazione scritto

in linguaggio FORTRAN che utilizza il solutore approssimato di Roe [10] per simulare flussi in-

stazionari uni-dimensionali. Il codice è stato opportunamente modificato aggiungendovi il modello

termodinamico del plasma, così come recentemente proposto da D’Angola et al. [3]. Uno degli

elementi chiave nella generalizzazione del metodo numerico è costituito dalla “linearizzazione

conservativa”, necessaria a garantire che anche le equazioni “discretizzate” conservino massa,

energia e quantità di moto. Utilizzando il codice di simulazione così modificato è stato possibile

condurre una valutazione comparativa delle capacità previsionali dei due modelli: plasma e gas

perfetto, al variare delle condizioni iniziali ed al contorno. Sono inoltre state implementate le

opportune condizioni al contorno necessarie per simulare adeguatamente flussi subsonici.

La seconda parte del lavoro di tesi è consistito nell’estensione della predetta metodologia allo

spazio bidimensionale. Il passaggio da una a due dimensioni spaziali comporta notevoli

complicazioni di carattere teorico, ancor prima che numerico. Tali problematiche sono state

affrontate mediante un attento lavoro di ricerca bibliografica che ha portato alla formulazione di un

originale modello teorico-numerico adatto alla simulazione del flusso di plasmi comprimibili su

reticoli di calcolo a maglie triangolari, mediante schemi numerici di distribuzione del residuo o

Fluctuation Splitting [8].

4

Capitolo 1

Modello termochimico per i gas ad alta

temperatura

5

1.1 Il plasma

Con il termine plasma viene indicato un gas globalmente neutro, costituito da specie sia cariche che

neutre le quali possono essere: elettroni, ioni, atomi o molecole. Lo stato di plasma può essere

ottenuto inducendo in un gas costituito da particelle neutre il fenomeno della ionizzazione per

mezzo di una scarica elettrica o di un aumento di temperatura. Un parametro molto importante per

la determinazione delle caratteristiche di un plasma è il grado di ionizzazione, ossia la frazione di

specie ionizzate rispetto alle specie totali. Un plasma che presenta un grado di ionizzazione basso

viene detto debolmente ionizzato, mentre uno che presenta un grado di ionizzazione prossimo

all’unità viene detto completamente ionizzato [1].

Nei gas costituiti da specie neutre, si assume che le molecole o gli atomi interagiscano tra di loro

solo per mezzo di collisioni binarie e che quindi seguano traiettorie rette tra un urto ed il successivo.

Nei plasmi la presenza di specie cariche rende rilevanti le interazioni a lungo raggio di tipo

coulombiano, dovute ai campi elettromagnetici, per cui una singola particella risente

contemporaneamente della presenza di tutte le altre. A causa di tale caratteristica le particelle

presenti in un plasma possono presentare dei moti estremamente complessi, caratterizzati da aspetti

collettivi [2].

Spesso il plasma viene considerato come il quarto stato della materia, è importante sottolineare però

che il passaggio di stato da gas neutro a plasma non presenta le discontinuità di fase tipiche dei

passaggi di stato canonici (gas-liquido, liquido-solido, gas-solido).

Negli ultimi decenni le applicazioni industriali del plasma hanno registrato un considerevole

sviluppo. In particolare tecnologie sviluppate in ambito metallurgico, come ad esempio il taglio e la

saldatura al plasma, o in ambito aerospaziale, come gli scudi termici per il rientro dei veicoli

Figura 1.1 Traiettorie tipiche di una particella neutra in un gas (a) e di una particella carica in un plasma (b) [2]

(a) (b)

6

spaziali o il controllo MHD per i flussi ipersonici, richiedono modelli numerici sempre più accurati

per lo studio e la simulazione del comportamento sia termodinamico che fluidodinamico dei plasmi.

Per molte di tali applicazioni è possibile assumere l’ipotesi di equilibrio locale sia chimico che

termodinamico, è possibile quindi modellare il comportamento del plasma considerando soltanto

due variabili termodinamiche, come ad esempio la pressione o la temperatura [3].

1.2 Modello termodinamico

All’aumentare della temperatura, il modello valido per i gas caloricamente perfetti risulta

inadeguato poiché a causa dell’attivazione dell’energia vibrazionale, la dipendenza dell’energia

interna dalla temperatura diventa di tipo non lineare. Alle alte temperature inoltre non è possibile

considerare il gas nemmeno termicamente perfetto, poiché la presenza di reazioni chimiche altera la

composizione del gas.

Il modo più generale di trattare un gas reale1 è quindi quello di considerarlo in non equilibrio sia

termico che chimico, tale trattazione comporta però notevoli svantaggi per quanto riguardo il costo

computazionale. L’introduzione dell’ipotesi di equilibrio termochimico locale permette quindi di

ottenere delle simulazioni dei gas ad alta temperatura sufficientemente accurate, con un costo

computazionale relativamente basso. Da considerazioni di tipo termodinamico, si ottiene che un

flusso in equilibrio termochimico locale presenta proprietà locali funzioni di due sole variabili di

stato indipendenti [4].

Nel modello utilizzato nel presente lavoro, il gas viene considerato come una miscela di CN

specie

chimiche, ognuna delle quali viene considerata come un gas termicamente perfetto; le proprietà

termodinamiche della miscela vengono quindi espresse in funzione di quelle delle singole specie.

L’energia interna per unità di massa e

è pari a:

)(,1

TeTee i

i

ii

CN

(1.1)

dove i e ie sono rispettivamente la densità e l’energia interna per unità di massa della i-esima

specie, T la temperatura e la densità della miscela definita come somma delle densità delle

singole specie:

1 Nel presente lavoro verrà indicato come “gas reale”, un gas che non è né termicamente e né caloricamente perfetto.

7

CN

i

i

1

(1.2)

L’energia interna per unità di massa viene espressa in diversi modi a seconda delle specie chimiche,

ad esempio per un gas biatomico termicamente perfetto può essere espressa come somma

dell’energia traslazionale, rotazionale, vibrazionale, elettronica e dal calore di formazione,

f

iielT

iiii he

e

RTRTe

i

,/12

5)(

(1.3)

dove il contributo vibrazionale è stato approssimato a quello di un oscillatore armonico semplice,

con kh ii pari alla temperatura vibrazionale caratteristica della molecola considerata. La

costante della specie i-esima è pari a ii mR / , in cui è la costante universale dei gas ed im è

la massa molare della specie i. Il contributo dovuto all’energia elettronica viene generalmente

trascurato. Per un gas monoatomico l’energia interna per unità di massa viene espressa come:

f

iielii heTRTe ,2

3)(

(1.4)

mentre per un elettrone libero vale la seguente espressione:

f

iii hTRTe 2

3)(

(1.5)

La pressione P può essere ottenuta combinando le equazioni di stato per le singole specie chimiche

con la legge di Dalton:

RTTRTPPCN

i

iii 1

,

(1.6)

dove

i

N

i

i

i RRRC

1

. Si ottiene quindi che l’entalpia della miscela può essere determinata

come:

TRTePeThh iii ,/,

(1.7)

8

Infine si ottiene l’equazione dell’energia interna per unità di volume:

CN

i

iii TeeT1

, .

Dalle equazioni precedenti si evince che le proprietà termodinamiche dipendono sia dalla

temperatura che dalla composizione chimica. In condizioni di equilibrio locale, sia T che i

dipendono soltanto da due variabili termodinamiche.

1.3 Modello chimico

La composizione chimica di una miscela di gas può essere determinata risolvendo un sistema non

lineare di 1CN equazioni le cui incognite sono T e i

[4]. Supponiamo che il numero di

elementi atomici, inclusi gli elettroni liberi, sia EN , allora la generica specie chimica iX può

essere espressa come segue:

E

ENijii

N

a

j

aai ZZZX,,1,

1

(1.8)

dove jZ indica l’elemento j-esimo e jia ,

il numero di elementi j presenti nella specie i. La k-esima

reazione invece può essere espressa come:

i

N

i

p

kii

N

i

r

ki XXCC

1

,

1

,

(1.9)

dove r

ki, e p

ki, sono rispettivamente i coefficienti stechiometrici dei reagenti e dei prodotti della k-

esima reazione. Dalla conservazione dei nuclei e delle cariche si ottiene la seguente relazione:

E

t

N

i

i

i

jiN

i

i

i

jiNj

m

a

m

a CC

...1,

01

,

1

,

(1.10)

9

dove l’espressione a destra rappresenta la composizione della miscela ad un istante fissato 0t .

Inoltre per ogni reazione è possibile scrivere l’equazione di equilibrio seguente:

C r

ki

C pki

N

i

ii

N

i

ii

p

k

m

m

K

1

1

,

,

, EC NNk 1

(1.11)

in cui tKK pp

kk è la costante di equilibrio della reazione k-esima alla pressione p .

Lo stato termochimico di una miscela di gas è noto dalla risoluzione del sistema di equazioni dato

dalle (1.6 - 1.11). Tali equazioni costituiscono un sistema notevolmente complesso di equazioni non

lineari accoppiate tra di loro, la cui risoluzione può essere ottenuta utilizzando diversi algoritmi. Nel

presente lavoro è stato considerato l’algoritmo presentato da D’Angola et al. [3], che consiste nella

risoluzione di un’equazione di equilibrio alla volta. La sequenza delle reazioni viene determinata a

partire dalla reazione più lontana dall’equilibrio.

1.4 Modello chimico per l’aria

Per applicare le relazioni mostrate nelle sezioni precedenti all’aria è necessario scegliere un modello

chimico per la composizione dell’aria. Essendo l’aria secca composta all’incirca per il 78.08% V/V

da azoto, per il 20.95% V/V da ossigeno e per il restante 0.97% V/V da argon ed altri gas, è

possibile considerare l’aria costituita esclusivamente da azoto ed ossigeno.

La scelta dell’accuratezza del modello chimico, ossia delle specie chimiche da trascurare, dipende

fortemente dai range di pressione e di temperatura considerati. Di seguito vengono mostrati alcuni

modelli riportati in letteratura [4]:

L’Air Mixture 1 (AM1) è composta dalle seguenti specie neutre

ΝΟΝ,Ο,,Ν,Ο 22

legate dalle seguenti reazioni chimiche:

OOO2

10

NNN2

ONON

L’Air Mixture 2 (AM2) comprende anche alcune specie ionizzate e gli elettroni

e,ΝΟΝΟ,Ν,Ο,,Ν,Ο 22

inoltre presenta la seguente reazione aggiuntiva

eNNO

L’Air Mixture 3 (AM3) è composta dalle seguenti specie

e,Ν,Ο,Ν,Ο,ΝΟΝΟ,Ν,Ο,,Ν,Ο2222

ed è caratterizzata dalle seguenti reazioni chimiche aggiuntive:

e2OOO

e2NNN

eOO

eNN

Tali miscele presentano dei range di validità limitati sia per quanto riguarda la temperatura che per

quanto riguarda la pressione. Considerando un range più ampio è necessario introdurre delle

ulteriori specie chimiche aumentando così notevolmente la complessità del modello:

Miscela presentata in [3]

e,ΝΟΝΟ,,Ο,Ν,Ο,Ν,Ο,Ν,Ν,Ο,Ν,ΟΝ,Ο,,Ν,Ο 443322

22 22

Tale miscela presenta una validità molto estesa sia per la temperatura che per la pressione (50 -

60.000 K e 0.01 – 100 atm). Nelle figure 1.2 e 1.3 viene mostrata la concentrazione delle specie

chimiche della miscela precedente alla pressione di 1 bar e al variare della temperatura. Si nota

11

come al crescere della temperatura, le concentrazioni delle specie trascurate dal modello AM3

diventano rilevanti.

Figura 1.2 Frazione molare delle specie chimiche, ad alta concentrazione, presenti nell’aria al variare della temperatura

Figura 1.3 Frazione molare delle specie chimiche, a bassa concentrazione, presenti nell’aria al variare della temperatura

12

1.5 Fit delle proprietà termodinamiche

Il fit delle grandezze termodinamiche, funzioni delle specie chimiche presenti nell’aria, in un range

esteso di temperatura e pressione (50 – 60000 K e 0.01 – 100 atm), rappresenta un problema

complesso a causa del comportamento non monotono delle variabili in funzione della temperatura.

Per tenere in conto di tale comportamento sono state introdotte le seguenti funzioni [3]:

Gaussiana

)exp(,; 2qcT ,

cTq (1.12)

iii cTT ,;

Sigmoide

)exp()exp(

)exp(,;

qq

qcT

,

cTq (1.13)

iii cTT ,;

Funzioni speciali

Figura 1.4 Frazione molare delle specie chimiche della miscela AM3 al variare della temperatura [4]

13

wT

cacT exp,; (1.14)

iiiii wcaTT ,,,;

waTwaT ,; (1.15)

iii waTT ,;

La dipendenza delle funzioni precedenti dalla pressione viene espressa considerando i parametri di

fit iiii wca ,,, , come funzioni polinomiali del logaritmo della pressione. Il generico parametro di

fit x si può esprimere quindi come

n

j

j

j px0

log (1.16)

Utilizzando le funzioni 1.12 – 1.14 è possibile esprimere le quantità termodinamiche in funzione di

due sole variabili, la temperatura T e la pressione p. L’espressione della densità 3/ mkg è la

seguente:

MRT

P (1.17)

dove M è la massa molare media molkg / espressa dalla seguente relazione:

6

1

0

j

jjacM (1.18)

Le espressioni dell’energia interna specifica e kgJ / e dell’energia interna per unità di volume

3/ mJ si ricavano a partire da quella dell’entalpia:

TaTchj

jj

j

j

j

7

1

2

1

(1.19)

phe (1.20)

14

phe (1.21)

In Appendice A vengono riportati i coefficienti utilizzati per i fit e i codici di calcolo Matlab

utilizzati.

Figura 1.5 Densità al variare della temperatura per differenti valori di pressione (0.01, 0.1, 1, 10, 100 atm)

Figura 1.6 Entalpia al variare della temperatura per differenti valori di pressione (0.01, 0.1, 1, 10, 100 atm)

15

Note le funzioni Tp, e Tp, , è stato possibile invertirle ottenendo delle tabelle in cui

la pressione e la temperatura sono state espresse in funzione della densità e dell’energia interna per

unità di volume:

,

,

,

,

TT

pp

Tp

Tp

Conoscendo ,pp in forma tabellare inoltre è stato possibile determinare numericamente le

derivate parziali della pressione rispetto alla densità , e rispetto all’energia interna per

unità di volume ,

p (1.22)

p (1.23)

L’accuratezza delle tabelle, ossia il numero di nodi della griglia, è stata scelta in maniera tale da

limitare l’errore d’interpolazione e contemporaneamente non compromettere l’efficienza

computazionale.

Figura 1.7 Pressione al variare dell’energia interna per unità di volume e della densità

16

Tali tabelle sono state validate effettuando un confronto con i fit ottenuti da Tannehill et al. [5]. In

tale articolo l’equazione di stato utilizzata è la seguente

1~1

e

heehp (1.24)

in cui eh~ . La relazione usata per il fit di ~ viene riportata di seguito:

3

8

2

7

2

6

2

54321~ ZaYZaZaYaYZaZaYaa

16151413

1211109

exp1 aYaZYaa

YZaZaYaa

(1.25)

dove

292.1log10 Y (1.26)

4.78408log10 Z (1.27)

I coefficienti dell’equazione 1.25 sono stati riportati in appendice A.

Figura 1.8 Pressione al variare dell’energia interna per unità di volume e a densità costante(1 kg/m3 ), confronto tra il fit NASA e le Tabelle

17

Capitolo 2

Le equazioni della fluidodinamica

18

2.1 Le equazioni di conservazione

Le equazioni che governano la fluidodinamica sono basate sulle seguenti leggi di conservazione

fondamentali:

1. Conservazione della Massa

2. Conservazione del Momento

3. Conservazione dell’energia

L’equazione scalare che si ottiene applicando la conservazione della massa viene detta equazione di

continuità. La conservazione del Momento non è nient’altro che la seconda legge di Newton,

pertanto essa è un’equazione vettoriale. Infine l’equazione di conservazione dell’energia consiste

nel primo principio della termodinamica. Per chiudere il sistema di equazioni derivanti dalle leggi

di conservazione è necessario considerare delle relazioni aggiuntive tra le proprietà del fluido; un

esempio tipico è l’equazione di stato che lega tra loro le proprietà termodinamiche del fluido.

Inoltre è necessario specificare le condizioni iniziali e al contorno.

Le leggi di conservazione si possono esprimere in maniera del tutto generale per grandezze sia

scalari che vettoriali [6]. Consideriamo un generico volume di controllo VC fissato nello spazio e

limitato da una superficie SC , l’equazione di conservazione della generica quantità Q in forma

euleriana si può scrivere come:

SC

S

VC

V

SCVC

ddSdQdt

nSnF (2.1)

Figura 2.1 Volume di controllo e superficie di controllo

19

Nella (2.1) F rappresenta il vettore dei flussi attraverso la superficie infinitesimale d , VS e SS

rappresentano il termine sorgente rispettivamente di volume e di superficie ed n il vettore normale

alla superficie d . Il vettore dei flussi si può scomporre in un termine convettivo pari a

QC uF

e in un termine diffusivo determinato dalla legge di Fick

qD F

dove è un coefficiente di diffusione, è la densità del fluido ed è tale che qQ .

Poiché il volume di controllo è costante nel tempo la derivata rispetto al tempo si può portare

all’interno del segno di integrale come segue

VCVC

dt

QQd

t, considerando inoltre il teorema

di Gauss è possibile esprimere gli integrali di superficie presenti nella (2.1) in integrali di volume:

SC

S

VC

V

SCVC

ddSddt

QSF (2.2)

Quest’ultima espressione permette di scrivere l’equazione di conservazione in forma differenziale:

SVSt

QSF

(2.3)

Nel caso in cui la quantità conservata sia un vettore Q , i flussi e il termine sorgente di superficie

diventano dei tensori, mentre il termine sorgente di volume diventa un vettore; si ottengono quindi

le seguenti espressioni per le equazioni di conservazione in forma integrale e differenziale:

SC

S

VCSCVC

ddddt

nSSnFQ

V (2.4)

St

SSFQ

V

(2.5)

20

I termini convettivi e diffusivi del flusso in questo caso si esprimono nel modo seguente

uQF C

qF D

È importante sottolineare che le equazioni (2.1) e (2.4) restano valide anche in presenza di

discontinuità delle grandezze che caratterizzano il flusso, come ad esempio le onde d’urto nel caso

di un flusso comprimibile. Le forme differenziali (2.3) e (2.5) delle equazioni di conservazione sono

valide solo in regioni in cui le proprietà del fluido variano in modo continuo; per cui il passaggio

dalla forma integrale a quella differenziale è lecito solo sotto l’ipotesi che la soluzione delle

equazioni di conservazione sia continua.

2.2 Equazione di conservazione della massa

L’equazione di conservazione della massa si ottiene dalla generica legge di conservazione scalare

(2.1) sostituendo a Q la massa per unità di volume, ossia la densità . Considerando il fluido

omogeneo i termini sorgente si annullano poiché non sono presenti reazioni chimiche, inoltre in tal

caso il flusso è dato unicamente dalla componente convettiva.

0

SCVC

ddt

nu

(2.6)

Attraverso il teorema della divergenza si ottiene:

0

VCVC

ddt

u

(2.7)

e quindi in forma differenziale

0

u

t (2.8)

21

Introducendo l’operatore di derivata sostanziale o convettiva definita come

u

tDt

D: si

ottiene la seguente espressione non conservativa dell’equazione di continuità

0 u

Dt

D (2.9)

2.3 Equazione di conservazione della quantità di moto

Dalla seconda legge della dinamica si ottiene che la variazione della quantità di moto all’interno del

volume di controllo è data dal flusso della quantità di moto attraverso la superficie più la somma

delle forze esterne V

iF agenti sul volume di controllo e delle forze S

jF agenti sulla superficie di

controllo. Dall’equilibrio delle forze agenti si può scrivere quindi in forma euleriana l’equazione di

conservazione in forma integrale

j

S

j

i

V

i

SCVC

ddt

FFnuuu

(2.10)

esplicitando le forze esterne si ottiene:

SCVCSCVC

ddddt

nTfnuuu

(2.11)

dove f rappresenta la generica forza esterna agente sull’elemento infinitesimo di volume d e

nT il generico vettore degli sforzi agenti sull’elemento infinitesimo di superficie d .

L’espressione del tensore degli sforzi dipende dalla natura del fluido considerato, in particolare

assumendo il fluido Newtoniano esso è dato dalla somma di un contributo dovuto alla pressione

isotropica p e da un tensore degli sforzi viscosi τ ; si può scrivere quindi che τIT p .

Sostituendo l’espressione del tensore degli sforzi nella (2.11) si ottiene:

SCSCVCSCVC

ddpdddt

nτnIfnuuu

(2.12)

22

La forma differenziale della (2.12) si ottiene applicando il teorema della divergenza come segue

VCVC

dpdt

fτuuu

(2.13)

fτuuu

p

t (2.14)

Considerando l’operatore di derivata sostanziale e l’equazione di conservazione della massa si

ottiene la seguente forma non conservativa per l’equazione della quantità di moto

fτu

pDt

D (2.15)

Nel caso di fluidi Newtoniani il tensore degli sforzi viscosi τ può essere espresso dalla seguente

relazione costitutiva:

DIuτ 2 (2.16)

in cui e sono due coefficienti di viscosità e t

uuD 2

1.

Infine adottando l’ipotesi semplificativa di Stokes secondo cui 32 si ottiene la cosiddetta

equazione di Navier-Stokes per flussi comprimibili.

DIufu

23

2p

Dt

D (2.17)

2.4 Equazione di conservazione dell’energia

L’equazione di conservazione dell’energia è una diretta conseguenza del primo principio della

termodinamica. La variazione dell’energia totale del fluido contenuto nel volume di controllo è pari

alla somma del lavoro compiuto sul sistema nell’unità di tempo L e del calore fornito al sistema

nell’unità di tempo Q . In forma euleriana si può scrivere:

23

QLdEdt

E

SCVC

nu

(2.18)

dove 2uu eE è l’energia interna totale. La potenza scambiata dal sistema con l’ambiente

esterno si può scomporre in due componenti, la prima dovuta all’azione delle forze di volume VL e

la seconda dovuta all’azione delle forze di superficie SL .

VC

V dL uf

e

SC

S dL unT

(2.19)

Analogamente è possibile scomporre la potenza termica scambiata dal sistema in due termini VQ e

SQ . Il termine VQ è un termine di generazione di calore interna legato ad esempio a delle reazioni

chimiche o a riscaldamenti di tipo ohmico dovuti a scariche elettriche, SQ è invece un termine

dovuto allo scambio termico che avviene attraverso la superficie di controllo.

VC

V dqQ

e

SC

S dQ nΦ

(2.20)

con q si è indicato il calore per unità di volume generato, mentre con Φ il flusso termico per unità

di area che attraversa la superficie di controllo. La legge di Fourier permette di legare il flusso

termico, in modo proporzionale al gradiente di temperatura T mediante la seguente relazione:

TkΦ

(2.21)

dove è la conducibilità termica. Tale coefficiente è legato in modo proporzionale alla viscosità

dinamica mediante il numero di Prandtl PrCp , per cui nel caso di flussi inviscidi il flusso

termico Φ è nullo. Sostituendo le espressioni (2.19-2.21) nella (2.18) si giunge alla forma integrale

dell’equazione di conservazione dell’energia

SCSCVCSCVC

dpdddEdt

E unIunτufnu

SC VC

dqdT n (2.22)

24

e attraverso il teorema di Reynolds si ottiene la sua forma differenziale:

qTpEt

E

uuτufu (2.23)

Analogamente a quanto fatto per le equazioni di conservazione della massa e della quantità di moto,

si può introdurre l’operatore di derivata sostanziale nella (2.23) ottenendo così la seguente forma

non conservativa per l’equazione di conservazione dell’energia:

qTpDt

DE uuτuf (2.24)

Un’ulteriore formulazione alternativa si ottiene introducendo l’entalpia totale pEH

qTHt

E

uτufu (2.25)

2.5 Equazioni di Navier-Stokes

Le equazioni (2.8), (2.14) e (2.23) vengono denominate equazioni di Navier-Stokes, esse

costituiscono il modello di descrizione della dinamica dei fluidi più generale possibile. È possibile

esprimere tali equazioni in maniera compatta nel seguente modo [6]:

SU

FF v-

t (2.26)

dove U è il vettore colonna (5x1) delle variabili conservative

u

U

E

25

F e vF sono dei tensori (5x3) che rappresentano rispettivamente i flussi inviscidi e viscosi

Iuu

u

u

p

H

F

τ

τu

0

TvF

ed S è il vettore dei termini sorgente

f

fu

0

S

q .

Le equazioni descritte dalla (2.26) costituiscono un set di equazioni alle derivate parziali non lineari

accoppiate con soluzione analitica solo in alcuni casi particolari. In particolare l’equazione di

conservazione della massa è un’equazione alle derivate parziali del primo ordine di tipo iperbolico,

mentre le equazioni di conservazione della quantità di moto e dell’energia, a causa della presenza

dei termini viscosi, sono del secondo ordine di tipo parabolico-iperbolico. Nelle regioni in cui

l’effetto dei termini viscosi è rilevante, come ad esempio in prossimità delle pareti solide, il

comportamento di tali equazioni è di tipo parabolico, mentre nelle restanti regioni è di tipo

iperbolico.

L’equazione (2.26) può essere resa adimensionale mediante un’opportuna normalizzazione delle

variabili fluidodinamiche, ciò permette di esprimere le equazioni di Navier-Stokes in funzione di

alcuni parametri adimensionali come il numero di Reynolds o il numero di Mach. Ad esempio si

può utilizzare la seguente procedura di normalizzazione [7]:

zyxL

zyx ,,1

,, *** , VL

tt* ,

V

uu

* ,

*

, 2

*

V

pp

,

T

TT * ,

2

*

V

ee ,

*

,

*

in cui l’asterisco indica le variabili adimensionalizzate, l’indice ∞ le condizioni di flusso

indisturbato ed L è la grandezza caratteristica del flusso usata nella definizione del numero di

Reynolds

26

LVRe .

Si ottiene quindi la seguente forma per le equazioni (2.21)

***

*

*

SU

F - F *

v*

t (2.27)

dove

**

**

*

*

u

U

E

**

******

f

fu

0

S

q

Iuu

u

u

****

***

**

p

H

*F

*

2

****

1

1

τ

τu

0

PrM

T

Re

*vF

(2.28)

Nelle (2.28) sono state introdotte delle variabili adimensionali aggiuntive le cui espressioni vengono

riportate di seguito

2

**** uu eE ,

*

***

pEH , L

V

ττ

*, L

V 2

*

f

f , LV

qq

3

*

pCPr ,

RT

VM

2.6 Equazioni di Eulero

Dalle equazioni (2.28) si nota che per valori del numero di Reynolds elevati i termini viscosi

diventano trascurabili, in tal caso è possibile considerare il flusso inviscido e privo di scambio

termico con l’esterno. Le equazioni che ne derivano in tale caso vengono dette equazioni di Eulero

e presentano un comportamento puramente iperbolico. L’equazione (2.26) si riduce quindi alla

seguente espressione:

27

SU

F

t (2.29)

Il tensore dei flussi F può essere scomposto considerando le sue componenti lungo gli assi

coordinati cartesiani come segue

zy HeGeFex F (2.30)

Nella (2,30) i vettori F , G ed H rappresentano le componenti del flusso lungo rispettivamente gli

assi x , y e z . Di seguito vengono esplicitati considerando le componenti del vettore velocità

twvu ,,u

uw

uv

pu

uH

u

2F ,

vw

pv

uv

uH

u

2

G ,

pw

vw

uw

uH

u

2

H (2.31)

si ottiene quindi la seguente espressione per l’equazione (2.29)

SHGFU

zyxt (2.32)

Nel caso di flusso comprimibile il campo di forza esterno dovuto all’attrazione gravitazionale non

da contributi rilevanti per cui generalmente viene trascurato, inoltre se non sono presenti sorgenti di

calore interne l’intero termine sorgente S può essere trascurato ottenendo le equazioni di Eulero

omogenee

0

zyxt

HGFU (2.33)

Nel caso unidimensionale si ottiene la seguente equazione

28

0

xt

FU (2.34)

Il sistema di equazioni definito dalla (2.29) necessita di una relazione costitutiva che definisca le

proprietà termodinamiche del fluido considerato. Tipicamente viene considerata l’equazione di stato

(EOS) del gas Tpp , accompagnata dalla definizione dell’energia interna Tpee , .

Come è stato già esposto nel paragrafo 2.1, le soluzioni delle equazioni di conservazione in forma

differenziale devono essere delle funzioni di tipo continuo, mentre nel caso della forma integrale la

soluzione può ammettere anche delle discontinuità di prima specie. Per il sistema di equazioni in

forma differenziale le soluzioni vengono dette forti, mentre per quello in forma integrale vengono

dette deboli. Nel caso in cui il campo fluidodinamico ammetta delle discontinuità, nelle varie zone

di continuità, si avranno diverse soluzioni forti. Tali soluzioni continue andranno raccordate tra di

loro in corrispondenza delle discontinuità mediante delle opportune relazioni note con il nome di

condizioni di Rankine-Hugoniot. Bisogna notare inoltre che nel caso delle equazioni di Eulero i

processi irreversibili legati alle viscosità non sono presenti, pertanto è necessario considerare delle

ulteriori condizioni che permettano di escludere le soluzioni non fisiche; tali condizioni prendono il

nome di condizioni entropiche.

29

Capitolo 3

Fluctuation Splitting per gas reale

30

3.1 Introduzione

Per la risoluzione numerica delle equazioni di conservazione iperboliche è di fondamentale

importanza la scelta del modello di discretizzazione dei termini convettivi. A seconda del tipo di

discretizzazione del dominio spaziale o del tipo di rappresentazione funzionale delle variabili

dipendenti, si ottengono diversi tipi di approcci.

L’approccio “tradizionale” consiste in una discretizzazione ai Volumi Finiti (FV) centrata sulle

celle, e su una rappresentazione discontinua delle variabili dipendenti. L’interazione tra celle

adiacenti viene quindi modellata considerando un problema di Riemann per ogni interfaccia. Una

discretizzazione spaziale centrata sui vertici piuttosto che sulle celle, ed una rappresentazione delle

variabili di tipo continuo sono legati al cosiddetto approccio Fluctuation Splitting (FS) o Residual

Distribution (RD). Tale approccio presenta notevoli vantaggi per quanto riguarda la

discretizzazione dei termini convettivi nel caso multidimensionale.

Di seguito dopo una breve descrizione degli approcci FV e FS 1D verrà presentata la

generalizzazione per gas reale del metodo FS per le equazioni di Eulero 1D e 2D.

3.2 Approccio ai Volumi Finiti

L’approccio ai Volumi Finiti si basa su una discretizzazione del tipo mostrato in figura 3.1 [8]. In

tale tipo di discretizzazione ogni volume finito ovvero ogni cella, viene individuata in modo

univoco da un nodo. Nel caso unidimensionale la cella esimai è centrata sul nodo ix e viene

individuata dai punti 2121 ; ii xx , tali per cui 22121 iii xxx e 2121 iii xxh .

Il valore assunto dalle variabili dipendenti nei nodi equivale al valore medio nella cella

corrispondente, inoltre la rappresentazione funzionale delle variabili dipendenti viene considerata di

tipo lineare e discontinua a tratti, per cui all’interfaccia sono presenti delle discontinuità.

Seguendo tale approccio, le equazioni di Eulero 1D (Eq. 2.34) possono essere approssimate

mediante la seguente formula di discretizzazione:

02121

1

xx

n

i

n

i

nn

ii FFUU (3.1)

31

in cui n

i 21F e n

i 21F rappresentano i valori dei flussi calcolati rispettivamente all’interfacce 21ix e

21ix . A causa della rappresentazione funzionale discontinua, per calcolare i flussi alle interfacce è

necessario risolvere il seguente problema di Riemann [9]

0

0

0,

0,

0

x

x

se

se

x

xxt

R

L

UU

UU

FU

(3.2)

La risoluzione del problema di Riemann nel caso di equazioni di conservazione non lineari può

essere affrontata risolvendo in maniera approssimata il problema esatto o risolvendo in maniera

esatta un problema approssimato. A quest’ultima categoria appartiene il più conosciuto ed utilizzato

dei risolutori del problema di Riemann ossia quello presentato da Roe nei primi anni ’80 [10]. In

tale metodo l’equazione di conservazione esatta viene sostituita da una equazione approssimata di

tipo lineare.

Applicando la regola della catena all’equazione (3.2) si possono riscrivere le equazioni di Eulero

nella cosiddetta forma quasi-lineare:

0

xt

UA

U (3.3)

dove A è la matrice jacobiana del vettore dei flussi

Figura 3.1 Discretizzazione spaziale FV 1D [8]

32

U

FA

(3.4)

Nel metodo di Roe la matrice jacobiana esatta viene sostituita da una matrice costante A~

, funzione

degli stati LU e RU , si ottiene così il seguente problema di Riemann approssimato.

0

00,

~

x

xx

R

L

xt

U

UU

0UAU

(3.5)

La matrice RL UUAA ,~~

deve essere scelta in maniera tale da essere rappresentativa delle

condizioni locali LU e RU , ed inoltre deve rispettare il seguente set di proprietà, denominato da

Roe come Proprietà “U”:

(i). RL UUA ,~

costituisce un’applicazione lineare dallo spazio vettoriale U allo spazio

vettoriale F

(ii). RL UUA ,~

ha autovalori reali RLii UU ,~~ , mi ,...,1 , e un corrispondente set di

autovettori linearmente indipendenti m

rrr ,...,, 21

(iii). UAUUA RL ,~

se LU e RU U

(iv). RL UUA ,~

soddisfa la seguente relazione

UUUAF RL ,~

(3.6)

La proprietà (ii) assicura che il problema approssimato conservi il carattere iperbolico del sistema

non lineare originale. La proprietà (iii) assicura la consistenza con le leggi di conservazione, mentre

la proprietà (iv) è una condizione sufficiente per assicurare che il metodo sia conservativo.

Quest’ultima proprietà permette di derivare algebricamente la matrice A~

come funzione di uno

stato medio U~

detto stato medio di Roe tale per cui si può scrivere UAA~~

. Nel caso di gas

perfetto lo stato medio può essere definito in funzione di due sole variabili, grazie alla proprietà di

omogeneità del primo grado della funzione UF . Tipicamente le due variabili che vengono scelte

sono la velocità u e l’entalpia totale H , si ottiene quindi che

33

Hu~

,~~AA (3.7)

Sostituendo la (3.7) nella (3.6) si ottiene che le medie u~ e H~

vengono definite dalle seguenti

espressioni [10]:

RL

RRLL uuu

~

(3.8)

RL

RRLL HHH

~ (3.9)

3.3 Approccio Fluctuation Splitting

L’approccio Fluctuation Splitting o Residual Distribution, contrariamente all’approccio FV si basa

su una rappresentazione continua delle variabili dipendenti e sul tipo di discretizzazione mostrata in

figura 3.2. In tale caso le celle presentano come vertici i nodi della discretizzazione, in particolare si

ottiene che la cella avente indice 21i ha per vertici i nodi localizzati nei punti ix e 1ix ed ha

quindi ampiezza pari a iii xxh 121 . Nell’ intorno di ogni nodo della griglia viene definito un

volume di controllo, detto volume di controllo mediano, individuato dai punti 21ix e 21ix e avente

ampiezza 2121 iii xxh dove: 2121 iii xxx e 2121 iii xxx [8].

Poiché si assume che la rappresentazione funzionale sia di tipo continuo, l’approccio FS non

richiede l’utilizzo del modello di Riemann per la determinazione del flusso all’intercella.

Figura 3.2 Discretizzazione spaziale FS 1D [8]

34

Alla base dell’approccio FS o RD vi sono il concetto di “fluttuazione” o “residuo” e il concetto di

“segnale” [11]. La fluttuazione viene associata a ciascuna cella ed è pari all’integrale dei flussi

lungo la frontiera della cella, rappresenta quindi il bilancio dei flussi lungo la cella. All’equilibrio si

ottiene che la fluttuazione è nulla per ogni cella, per cui si può affermare che essa rappresenta una

misura di quanto ciascuna cella si trovi lontana dall’equilibrio. Il segnale rappresenta invece la

frazione di fluttuazione che viene “inviata” ad uno dei vertici che circonda la cella considerata.

Considerando l’equazione di convezione scalare con velocità costante

0

x

w

t

w (3.10)

si ottiene che la fluttuazione relativa alla cella 21i è pari a

ii

ix

ix

ix

ixi

wwdxx

wdx

t

w

1

11

2

1 (3.11)

La distribuzione della fluttuazione, ossia la determinazione dei segnali da inviare ai nodi i e 1i ,

viene fatta in base al segno della velocità . In particolare la fluttuazione viene scomposta in una

componente positiva

21i e in una negativa

21i , tali per cui

212121 iii .

0

0

,0

,2

1

21

i

i

0

0,

,0

2

121

ii

Quindi al nodo esimoi viene inviata la fluttuazione positiva della cella 21i e quella negativa

della cella 21i , come mostrato in figura 3.3.

Figura 3.3 Distribuzione della fluttuazione per l’equazione della convezione scalare

35

Si ottiene quindi il seguente schema di discretizzazione per l’equazione (3.10)

2121

1

ii

n

i

n

i x

tww (3.12)

in cui t viene scelto in maniera tale da rispettare la nota condizione di stabilità

1

x

tCFL (3.13)

Per quanto riguarda le equazioni di Eulero è possibile applicare il metodo sviluppato per le

equazioni iperboliche scalari considerando le seguenti tre fasi distinte [12]:

Decomposizione caratteristica

Linearizzazione conservativa

Distribuzione di tipo upwind

Nella prima fase il sistema di equazioni iperboliche accoppiate fra di loro, costituito dalle equazioni

di Eulero, viene trasformato in un sistema di equazioni di convenzione scalari. La seconda fase

consiste nel determinare una linearizzazione locale del problema (3.3) tale per cui venga rispettata

la condizione (3.6). Infine è necessario determinare uno schema di distribuzione del residuo che sia

di tipo upwind.

3.4 Decomposizione caratteristica delle equazioni di Eulero 1D

Si consideri le equazioni di Eulero 1D

0

xt

FU (3.14)

dove

u

E

U

e

pu

uH

u

2

F

36

Per un gas ideale è possibile esprimere il vettore di flussi inviscidi come funzione del vettore delle

variabili conservative UFF nel seguente modo [8]:

2

1

2

3

2

1

2

3

2

1

2

3

1

23

3

21

2

1

u

uu

u

u

u

u

u

uu

u

F

(3.15)

La 3.15 è stata ottenuta considerando che per un gas perfetto valgono le seguenti relazioni

21

1 2upE

(3.16)

1

2

3

2

2

21

21

u

uu

uEp (3.17)

e

2

2

2

1

21 uE

uEE

pEH

2

1

2

3

1

2

2

1

u

u

u

uH

(3.18)

Differenziando la (3.15) rispetto al vettore delle variabili conservative, si ottiene la matrice

jacobiana dei flussi che compare nella forma quasi-lineare delle equazioni di Eulero (3.3)

1

3

2

1

2

3

2

1

3

3

1

2

1

3

3

1

3

3

2

1

32

312

3

12

31

100

u

u

u

u

u

u

u

u

u

u

u

u

u

uu

A (3.19)

37

uu

uHuu

Hu

312

3

12

1

100

2

22

(3.20)

Poiché il sistema (3.14) è di tipo iperbolico, la matrice A presenta degli autovalori reali e un set

completo di autovettori linearmente indipendenti. In particolare nel caso unidimensionale presenta i

seguenti tre autovalori:

auu ,,0 (3.21)

dove a è la velocità del suono e si esprime come

pu

Ha

21

22 (3.22)

Gli autovalori possono essere rappresentati mediante la matrice diagonale ,,0diagΛ .

A tali autovalori vengono associati degli autovettori destri e sinistri raggruppati rispettivamente

nelle matrici R ed L

11

2

1

2

a

u

a

uu

ua

Hu

a

Huaa

R (3.23)

11

2

1

2

1

2

1

2

1

11

2

1

2

1

2

1

2

1

11

2

11

2

2

222

2

uaa

uu

a

uaa

uu

a

uaaa

u

L (3.24)

38

Poiché valgono le seguenti relazioni

RΛAR (3.25)

ΛLLA (3.26)

la matrice A può essere fattorizzata nel seguente modo:

LRA (3.27)

Tale scomposizione permette di diagonalizzare le equazioni di Eulero 1D rendendole un sistema di

equazioni di convezione scalari. Infatti premoltiplicando l’equazione (3.3) per la matrice degli

autovettori sinistri si ottiene

xtxtxt

UΛL

UL

ULA

UL

UA

UL

0

xt

WΛ

W (3.28)

in cui sono state definite le variabili caratteristiche W come:

ua

p

ua

p

pa

2

1

2

1

12

ULW (3.29)

Il set di equazioni (3.28) è costituito quindi da tre equazioni scalari di convezione non accoppiate

tra loro

0

x

w

t

w kk

k

(3.30)

di cui la prima rappresenta l’onda di entropia e le altre due le onde acustiche rispettivamente veloce

e lenta. Le (3.30) possono essere esplicitate come segue:

39

011

22

x

p

axu

t

p

at

(3.31)

011

x

u

x

p

aau

t

u

t

p

a (3.32)

011

x

u

x

p

aau

t

u

t

p

a (3.33)

Considerando l’operatore di derivata sostanziale introdotto nel paragrafo 2.2 , si può scrivere

0k

k

wDt

D (3.34)

dove

xtDt

D kk

: .

Dall’espressione (3.34) si ottiene che le variabili caratteristiche restano costanti lungo delle curve

che vengono dette curve caratteristiche. Tali curve, nel piano tx , hanno tangente pari a k

,,0dt

dx (3.35)

Figura 2.2 Curve caratteristiche nel piano x-t [8]

40

Nel caso di un gas reale, o in generale nel caso di un gas avente una generica equazione di stato, per

poter calcolare la matrice jacobiana A è necessario esprimere la pressione come una funzione

generica delle variabili conservative [4].

mEp ,ˆ, U (3.36)

dove um è la quantità di moto ed EE ˆ è l’energia interna totale per unità di volume.

L’equazione (3.36) non rappresenta una relazione di tipo termodinamico; infatti nel caso di gas in

equilibrio termodinamico la pressione dipende soltanto da due variabili termodinamiche. Esprimere

la pressione mediante la (3.36) permette quindi di determinare in maniera conveniente i termini

dalla matrice jacobiana dipendenti dall’equazione di stato.

Differenziando la (3.36) si ottiene

dmEdddp mE ˆ

ˆ (3.37)

in cui le derivate parziali della pressione rispetto alle variabili conservative sono pari a:

Em

p

ˆ,

,

m

E

E

p

,

ˆˆ

,

E

mm

p

ˆ,

(3.38)

Essendo la pressione funzione di due sole variabili termodinamiche, come ad esempio la densità e

l’energia interna epp , , o la densità l’energia interna per unità di volume ,pp , segue

che

Em

mˆ

(3.39)

Considerando la funzione (3.36), il vettore dei flussi inviscidi si esprime nel seguente modo:

321

1

2

3

3212

1

3

3

2

,,

,,

,ˆ,

uuuu

u

uuuuu

u

u

mEu

uH

u

F

(3.40)

41

si può quindi determinare la matrice jacobiana per il caso di gas reale come segue

uu

uHuHu

mE

mE

2

1

100

ˆ2

ˆ

UA (3.41)

Gli autovalori della matrice A presentano la stessa forma del caso di gas perfetto2, si ottiene infatti

,,03,2,1

auu ,,0 (3.42)

dove la velocità del suono è espressa come

E

uHa ˆ22 (3.43)

Gli autovettori destri e sinistri della matrice (3.41) vengono riportati di seguito

11

1

ˆ

2

a

u

a

uu

ua

Hu

a

HaH

aa

E

R (3.44)

2 Si riportano per completezza i passaggi algebrici:

0det IA

0121 2ˆˆˆˆ uuHuuHuuEEmEmE

023 ˆ23

ˆˆ222

ˆ3

EEEmEm uHuuHuuuuu

03 2222ˆ

22ˆ

3 auuauuuuuEmEm

03 222223 auuauu

0 auauu

42

12

1

22

1

12

1

22

1

ˆˆ2ˆ

ˆˆ2ˆ

2

ˆ

2

ˆ2

2

ˆ

uaa

uHa

ua

uaa

uHa

ua

uaa

uHa

EEE

EEE

EEE

L (3.45)

Nota la matrice degli autovettori sinistri è possibile diagonalizzare le equazioni di Eulero anche per

il caso di gas reale ottenendo così l’equazione (3.28)3.

Nel caso in cui l’equazioni di Eulero non siano omogenee, ossia nel caso in cui sia presente il

termine sorgente l’equazione (3.28) si scrive come segue

SW

ΛW ˆ

xt

(3.46)

dove LSS ˆ . In tale caso le variabili caratteristiche non restano costanti lungo le curve

caratteristiche, infatti l’espressione 3.34 diventa la seguente

kk

k

swDt

Dˆ (3.47)

in cui il termine ks indica il “tasso di variazione” della variabile caratteristica kw , lungo la curva

caratteristica esimak .

3Le variabili caratteristiche nel caso di gas reale presentano la stessa forma del caso di gas perfetto.

Si riportano di seguito i passaggi algebrici

mua

Ea

uHa

ua

mua

Ea

uHa

ua

mua

Ea

uHa

EEE

EEE

EEE

12

1ˆ22

1

12

1ˆ22

1

ˆ

ˆˆ2ˆ

ˆˆ2ˆ

2

ˆ

2

ˆ2

2

ˆ

ULW

ua

p

ua

pa

p

uupa

u

uupa

u

paa

uuma

Eaa

u

uuma

Eaa

u

mEaa

mE

mE

mE

2

1

2

1

2

1

2

1

2

1

2

12

1

2

1

2

1

2

1

1

2

1

2ˆ

22

1

2

1

2ˆ

22

1

ˆ12

2

2

ˆ

ˆ

ˆ

2

2

43

Di seguito viene esplicitato il termine S sia per gas reale

22

22ˆ

ˆ

ˆ

2

ˆ

fq

a

fq

a

qa

E

E

E

S

(3.48)

che per gas perfetto

22

122

1

1

ˆ

2

fq

a

fq

a

qa

S

(3.49)

Nel caso di flussi comprimibili è possibile considerare trascurabile il campo di forze esterne,

ottenendo quindi le seguenti espressioni semplificate per la (3.48) e la (3.49)

qa

qa

qa

E

E

E

2

2ˆ

ˆ

ˆ

2

ˆ

S

(3.50)

qa

qa

qa

2

12

1

1

ˆ

2

S

(3.51)

44

3.5 Linearizzazione conservativa

Assumendo una rappresentazione della soluzione di tipo lineare a tratti e di tipo continuo, come

rappresentato in figura 3.2, è necessario determinare una linearizzazione dell’equazione (2.34) tale

per cui la fluttuazione relativa alle variabili conservate 21iΦ sia valutata in maniera esatta [13]

21

11

2

1

i

ix

ix

ix

ixi

xx

dxx

dxt

FFUΦ (3.52)

in cui il gradiente dei flussi viene valutato come segue

dxxxx

ix

ixi

1

21

1 FF (3.53)

Utilizzando la forma quasi-lineare delle equazioni di Eulero (3.3) è possibile scrivere

21

11

2

1

~

i

ix

ix

ix

ixi

xx

dxx

dxx

UA

UA

FΦ (3.54)

in cui il gradiente del vettore delle variabili conservative viene mediato nella cella come segue

dxxxx

ix

ixi

1

21

1 UU (3.55)

ed A~

è un media sulla cella della matrice jacobiana dei flussi, valutata in maniera tale che venga

rispettatala la condizione di conservazione seguente

xx

UA

F ~ iiii UUAFF 11

~ (3.56)

45

Una linearizzazione conservativa nel senso della (3.56), garantisce che considerando l’intero

dominio di calcolo, i flussi interni si elidano in maniera tale da lasciare solo i contributi al bordo [8]

0

1

0

121

~FFUUA

N

N

i

iii (3.57)

La condizione (3.57) viene definita generalmente come proprietà “telescopica” del residuo ed è

fondamentale affinché lo schema considerato sia capace di determinare correttamente le variazioni

delle variabili in corrispondenza di soluzioni discontinue come gli urti.

La matrice 21

~iA viene considerata valutando la matrice jacobiana in corrispondenza di uno stato

medio μ~ , funzione degli stati dei vertici della cella iU e 1iU ,

121 ,~~ iii UUμAA (3.58)

3.5.1 Stato medio per gas perfetto

Per un gas perfetto lo stato medio può essere ottenuto considerando il parametro vettore seguente,

introdotto da Roe [10]

u

H

z

z

z 1

3

2

1

Z (3.59)

Tale scelta è particolarmente conveniente poiché ogni componente del vettore delle variabili

conservative e del vettore dei flussi inviscidi risulta essere una funzione bilineare delle componenti

di Z come espresso dalle (3.60) e (3.61)

31

21

2

1

zz

pzz

z

u

E

U

(3.60)

pz

zz

zz

pu

uH

u

2

3

32

31

2

F (3.61)

46

in cui la pressione è una funzione quadratica delle componenti di Z :

2

12

3

21

zzzp

(3.62)

Si ottiene quindi che le matrici jacobiane ZUB e ZFC sono funzioni lineari delle

componenti di Z :

13

312

1

0

1

002

zz

zzz

z

Z

UB (3.63)

312

23

13

111

0

0

zzz

zz

zz

Z

FC (3.64)

Considerando il parametro vettore con un andamento lineare nella cella si ottiene che il gradiente

xZ è costante ed è pari a:

21

11

21

1

i

iiix

ixi xdx

xxx

ZZZZ (3.65)

Considerando le equazioni (3.63-3.65) si ottiene che i gradienti medi (3.53) e (3.55) si esprimono

come segue:

dxxxx

ix

ixi

1

21

1 UU

dxxx

dxxx

ix

ixi

ix

ixi

1

21

1

21

11

Z

UZZ

Z

U

ZBZ

ZZ

UZ

xx

(3.66)

47

dxxxx

ix

ixi

1

21

1 FF

dxxx

dxxx

ix

ixi

ix

ixi

1

21

1

21

11

Z

FZZ

Z

F

ZCZ

ZZ

FZ

xx

(3.67)

Essendo B e C lineari rispetto a Z si ottiene che lo stato medio Z è pari alla media aritmetica tra

iZ e 1iZ .

2

2

2

2

11

11

1

1

iiii

iiii

ii

ii

uu

HH

ZZZ (3.68)

Sostituendo le (3.66) e la (3.67) nella (3.56) è possibile ottenere quindi l’espressione della matrice

A~

x

xxx

U

ZBZ

AZCZF ~

ZBAZC~

1

3

2

1

2

3

2

1

2

3

1

2

1

3

1

2

2

1

2

3

1

31

312

3

12

1

100

~

z

z

z

z

z

z

z

z

z

z

z

z

z

z

z

z

ZBZCZAA (3.69)

Si può notare come la (3.69) sia funzione unicamente dei rapporti 13 zz e 12 zz i quali coincidono

con le medie definite dalle (3.8) e (3.9)

RL

RRLL uu

z

zuu

1

3~ Z (3.70)

48

RL

RRLL HH

z

zHH

1

2~Z (3.71)

Lo stato medio per un gas perfetto è quindi individuato unicamente dal vettore Huii

~,~,~

1 UUμ .

La decomposizione caratteristica del problema linearizzato in maniera conservativo può essere

ottenuto valutando in corrispondenza dello stato medio Z , le espressioni analitiche presentate nel

paragrafo 3.4

auu ~~~,,~~

0 (3.72)

1~

~~1~

~~~

~~

~~~

~

~~

2

~

~

~

~

~1

~2

a

u

a

uu

ua

Hu

a

Hu

aa

ZRR (3.73)

in cui la velocità del suono a~ viene espressa come

2

~~1~

22 u

Ha (3.74)

e viene introdotta la densità media ~

2

2

12

~

RLz

Z (3.75)

Infine la matrice degli autovettori sinistri è pari a 11 ~~ RZRZLL .

La fluttuazione può essere quindi valutata nel seguente modo

21

1

2

1

~~

i

ix

ixi

xx

dxx

WΛR

FΦ (3.76)

49

dove ~

,~

,~~

0diagΛ e xx UZLW .

3.5.2 Stato medio per gas reale

Nel caso di un gas reale che presenti una equazione di stato generica, non è verificata la proprietà di

omogeneità di primo grado della funzione UF , per cui non è possibile individuare lo stato medio

considerando soltanto due variabili termodinamiche come nel caso di gas perfetto.

Per poter estendere la linearizzazione proposta da Roe ai gas reali sono stati proposti in letteratura

diversi approcci. È possibile distinguere tra tali approcci due famiglie principali [14]: una prima,

che comprende la maggior parte dei metodi proposti [15-18], basata su una forma quasi-Jacobiana

della matrice A~

, e una seconda famiglia in cui si considera una forma strettamente Jacobiana della

matrice A~

. Nel primo caso si suppone che la matrice linearizzata sia della forma seguente

qJμAA ~~~

(3.78)

in cui lo stato medio viene assunto come segue

mEqJ Hu ~

,~

,~

,~

,~~ˆμ (3.79)

In tale caso A~

viene detta in forma quasi-Jacobiana poiché i termini ~

, E

~ e m

~ non

rappresentano le derivate della funzione Up valutate in corrispondenza di uno stato medio U~

, ma dei parametri aggiuntivi che vengono introdotti dalla procedura di linearizzazione.

Un inconveniente di tale formulazione è che le derivate medie della pressione ~

, E

~ e m

~ sono

delle incognite artificiali e quindi non conservano il loro significato termodinamico. Ciò può portare

a delle inconsistenze nel caso in cui tali quantità vengano utilizzate per la determinazione di altre

variabili termodinamiche come la velocità del suono, la quale determina la struttura della

decomposizione caratteristica del problema linearizzato [19].

I metodi che appartengono alla seconda famiglia [14] vengono costruiti considerando lo stato medio

funzione di sole tre variabili indipendenti. Generalmente si suppone che lo stato medio sia del

seguente tipo

50

UμAμAA~~~~~

JJ (3.80)

dove 1,~~

ii UUUU . In questo caso le derivate della pressione sono funzioni dello stato medio e

quindi sono delle variabili dipendenti.

U~~

U~~

ˆˆ EE U

~~mm (3.81)

Lo stato medio viene ottenuto risolvendo il sistema di equazioni dato dalla (3.56), in cui si assume

che la matrice jacobiana abbia la forma data dalla (3.80). Poiché la prima equazione di tale sistema

viene identicamente soddisfatta ( uu ) è necessario introdurre un’equazione

supplementare per poter determinare lo stato medio. Tipicamente si considera la seguente equazione

mEp mE UUU

~~ˆ~~~ˆ (3.82)

In questo modo si ottiene un sistema di tre equazioni aventi tre incognite

Em~ˆ,~,~~

U .

Per i gas ad alta temperatura molto spesso non è disponibile l’espressione analitica dell’equazione

di stato, per cui non è possibile utilizzare i metodi basati sulla formulazione strettamente Jacobiana.

Nel caso della formulazione quasi-Jacobiana, lo stato medio non è univocamente determinato

poiché il numero di incognite definite dal vettore ausiliario mEqJ Hu ~

,~

,~

,~

,~~ˆμ

è maggiore

del numero di condizioni espresse dall’equazione (3.56).

L’imposizione della condizione (3.56) quindi non basta da sola ad assicurare l’univocità dello stato

medio come avviene nel caso di gas perfetto, per cui è necessario considerare delle condizioni

aggiuntive. Assumendo per u~ e H~

la stessa forma ottenuta per gas perfetto si ottiene che i

parametri aggiuntivi risultano legati tra loro in modo lineare dalla seguente espressione

mEp mE

~ˆ~~ˆ (3.83)

Tale espressione permette di ridurre il numero dei gradi di libertà della soluzione da tre a due. Infine

si impone la condizione (3.39) tra i parametri E

~ , m

~ e u~

Em u ˆ

~~~ (3.84)

51

per imporre la consistenza dei parametri con il loro significato termodinamico e per ridurre il

numero di gradi di libertà della soluzione ad uno. In seguito verrà dimostrato che assumendo per

l’equazione di stato l’espressione generica data dalla (3.36), si giunge alla stessa definizione dello

stato medio, pur utilizzando diversi metodi.

Considerando il parametro vettore definito dalla (3.59), è possibile ottenere la generalizzazione

dell’approccio originale di Roe, mostrato nel paragrafo 3.5.1. Per poter determinare le matrici C e

B è necessario conoscere il differenziale della pressione in termini dei differenziali delle

componenti del parametro vettore [4]. Dalla definizione di entalpia totale si ricava che

dpHdEd )(ˆ (3.85)

Sostituendo la (3.85) nella (3.37) si ottiene

udHdddp mE

E

)(1

1ˆ

ˆ

(3.86)

I termini presenti nell’equazione (3.86) sono facilmente esprimibili in funzione delle componenti

del parametro vettore

112 dzzd (3.87)

2113 dzzdzzud (3.88)

3112)( dzzdzzHd (3.89)

Assumendo che sia valida la (3.83) si ottiene la seguente discretizzazione

)(~

)(~~

~1

1ˆ

ˆ

uHp mE

E

(3.90)

e quindi

52

2ˆ1311ˆ231

~~)

~~~2(~

1

1zzzzzzzzp

EmEm

E

(3.91)

in cui è stato assunto che il parametro vettore vari linearmente nella cella come nel caso di gas

perfetto, e quindi 21 ii ZZZ . Dalla (3.91) segue in maniera banale che:

13

3

ˆ

ˆ

1

ˆˆ

312

1

ˆ

0

~1

~

~1

1~

1

~~2

002

)~

,~

,~

,(

zz

zzzzz

z

E

E

EE

m

mE

ZB (3.92)

3

ˆ

ˆ

1

ˆ

ˆ

ˆ

ˆ231

23

13

ˆ

1~

2~

1~

~

~1

~~~2

0

0

)~

,~

,~

,(

zzzzz

zz

zz

E

E

E

E

E

Em

mE

ZC (3.93)

infine si ottiene l’espressione della matrice linearizzata

1

3ˆ2

1

2

3

1

3

1

2

1

3ˆ

1

2

1

31

ˆ

2~~~

~~1

~

100

~~~,

~,

~,

~

z

z

z

z

z

z

z

z

z

z

z

z

z

z

mE

mEmE

BCZAA (3.94)

uu

uHuHu

mE

mEqJ

~2~~~~

~~~~1~~~~

100~

ˆ2

ˆ

μA (3.95)

in cui u~ e H~

rappresentano le usuali medie di Roe (3.70) e (3.71).

uz

z ~

1

2 Hz

z ~

1

3

53

È possibile giungere alla stessa espressione dello stato medio ottenuta attraverso l’utilizzo del

parametro vettore, sostituendo direttamente l’espressione (3.95) nella (3.56). Per la prima

equazione si ottiene l’identità seguente:

uu (3.96)

La seconda equazione è uguale a

EuuuHuHuHEm

ˆ~~1)(

~~~~~~)( ˆ (3.97)

e può essere riscritta usando la (3.91) come segue

)(~

~~~)()(~)(~~

uuuuHHuuuHE

m

(3.98)

utilizzando la (3.84) si ottiene

)()(~)(~~uHHuuuH (3.99)

La terza equazione è uguale a

EuuupuEm

ˆ~)(

~~2~~)( ˆ

22 (3.100)

sostituendo la (3.83) si ottiene

0)()(~2~ 22 uuuu (3.101)

La (3.99) e la (3.101) presentano come soluzione le stesse espressioni (3.70) e (3.71)ottenute per il

caso di gas perfetto.

Infine è possibile determinare lo stato medio imponendo le seguenti espansioni spettrali

nell’equazione (3.56)

54

3

1

~~~~~

p

ppp rWΛRF (3.102)

3

1

~~~

p

pprWRU (3.103)

in cui p

~ sono gli autovalori della matrice A

~ espressa dalla (3.95)

u~~

1 au ~~~2 au ~~~

3 (3.104)

pr~ sono i corrispondenti autovettori destri, raccolti nelle colonne della matrice R~

seguente

auauu

uaHuaHa

HE

~~~~~

~~~~~~~

~~111

~2

R (3.105)

e p

~ sono le intensità delle onde caratteristiche, cioè le componenti del vettore β

~:

ua

p

a

ua

p

a

a

p

uuaa

p

uuaa

pa

p

ˆ~~2

1

ˆ~~2

1

~

~)(~2

1~2

~)(~2

1~2

~

~~

2

2

2

2

ULβ (3.106)

in cui 1~~ RL e 1

ˆ ii .

La velocità del suono definita nel calcolo degli autovalori è pari a

EuHa ~

)~~(

~~ 22

(3.107)

Anche in questo caso per u~ e H~

si ottengono le espressioni ottenute da Roe per i gas perfetti.

55

3.5.3 Determinazione delle derivate della pressione

Per completare la formulazione generale qui proposta, è necessario valutare le derivate medie della

pressione (3.38) in termini delle derivate di opportune variabili termodinamiche in base

all’equazione di stato scelta. Nelle condizioni di equilibrio termodinamico locale, la pressione può

essere riferita ad altre due variabili termodinamiche attraverso una generica equazione di stato della

forma

),()ˆ,,( ipEmp (3.108)

dove la variabile i può essere l’energia interna, per unità di massa e o per unità di volume ε, oppure

la temperatura T. A seconda della scelta fatta, le derivate della pressione (3.38) e la matrice

Jacobiana del flusso assumono una forma differente, influenzando così la determinazione dello stato

medio di Roe.

La valutazione delle derivate medie della pressione può essere effettuata seguendo le fasi:

Scelta dell’equazione di stato nella forma (3.108);

Da considerazioni termodinamiche, dedurre le relazioni:

Emxeup

i

pD

i

xˆ,,,,...,,,,

(3.109)

che soddisfino la (3.37);

Assumere che le relazioni (3.109) presentino come valori medi:

Emxeup

i

pD

i

xˆ,,,,...~,~,~,

~,

~~

(3.110)

Inserire le relazioni (3.110) nella (3.83).

Come esempio si consideri un gas perfetto per cui vale

2

2

1ˆ)1(m

Ep (3.111)

56

Le derivate della pressione sono facilmente calcolabili

2

2

1u

um )1( 1ˆ E

(3.112)

Supponendo che le relazioni (3.112) valgano anche per i valori medi ed inserendole nella (3.83) si

ottiene

Euuup

)1()(~)1(~

2

1 2

(3.113)

Differenziando la (3.111) e sostituendo p nella (3.113) si ricava

0)()(~2~

2

1 22

uuuu

(3.114)

che è soddisfatta poiché vale l’equazione (3.101). Le derivate sopra scritte provano che per un gas

perfetto lo stato medio di Roe è definito in modo univoco e che la formulazione generale data dalle

equazioni (3.83) e (3.95) si riduce esattamente allo schema originale per gas perfetto.

3.6 Generalizzazione di Vinokur-Montagnè

Nel presente lavoro viene considerata la linearizzazione in forma quasi-Jacobiana presentata da

Vinokur e Montagnè [15], in cui vengono scelte come variabili termodinamiche indipendenti la

densità e l’energia interna per unità di volume

),( pp (3.115)

Differenziando la (3.115) si ottiene

dddp (3.116)

dove

57

p,

p

(3.117)

Ricordando che

2

2

1ˆ mE

(3.118)

si può esprimere l’equazione di stato come segue

2

2

1ˆ,m

EpU (3.119)

Applicando la regola della derivazione a catena alla (3.119) si ottengono le derivate della pressione

nella forma data dalla (3.109)

2

21 u , um ,

E~ (3.120)

Sostituendo i valori medi delle (3.120) nell’equazione (3.83) si ottiene

Emuup ˆ~~~~~~ 2

21 (3.121)

Dalla definizione di energia totale per unità di volume 2

21ˆ uE si ottiene

22

21 ~2~~~~ umuup (3.122)

Grazie all’equazione (3.101) i termini tra parentesi quadre si elidono, per cui si ha

~~p (3.123)

La relazione (3.123) rappresenta l’equivalente discreto della (3.116), inteso nell’ambito

dell’approccio quasi-Jacobiano. Tale relazione lineare non è sufficiente a definire in maniera

univoca lo stato medio.

58

Per definire in modo univoco i valori di ~ e ~ , Vinokur e Montagnè propongono una procedura

che utilizza le informazioni date dai due stati termodinamici L e R individuati in corrispondenza dei

nodi i ed 1i . Tali stati sono individuati unicamente dalle due variabili termodinamiche e .

Integrando l’equazione (3.116) lungo uno dei possibili percorsi tra gli stati L ed R e usando la

(3.123) si arriva alle relazioni generali seguenti

1

0)](),([~ dttt (3.124.a)

1

0)](),([~ dttt (3.124.b)

in cui il parametro t è normalizzato in modo che tL=0 e tR=1.

Come mostrato da Toumi [19] lo stato medio dipende fortemente dal percorso lungo il quale viene

eseguita l’integrazione. Vinokur e Montagné propongono di integrare le equazioni (3.124) lungo il

percorso lineare definito come segue:

tt

tt

R

R (3.125)

L’esatta valutazione degli integrali (3.124) è, in generale, laboriosa ed in pratica è necessaria

un’approssimazione. Si può pensare di approssimare i valori delle derivate a quelli del punto medio

Figura 3.4 Possibili percorsi tra gli stati termodinamici L ed R

59

),(ˆMMM (3.126.a)

),(ˆMMM (3.126.b)

in cui lo stato medio è definito da 2/)( RLM e 2/)( RLM , oppure secondo il metodo

dei trapezi

2/)(ˆRL (3.127.a)

2/)(ˆRL (3.127.b)

Infine, quando i due stati L ed R sono molto distanti, si può applicare la regola di Simpson

6/)4(ˆRML (3.128.a)

6/)4(ˆRML (3.128.b)

Una volta trovati i valori approssimati e è necessario determinare i valori ~ e ~ che

verificano la (3.123) e che contemporaneamente sono il più vicino possibile ai valori approssimati.

Tale procedura può essere formulata geometricamente come una proiezione del punto )ˆ,ˆ( sulla

retta definita dall’equazione (3.123).

Figura 3.5 Proiezione dello stato medio approssimato sulla retta (3.123)

60

In particolare ~ e ~ vengono definiti come i valori di e che minimizzano la funzione

seguente

22

2 ˆ

11

ˆ

ˆ

ˆ

1,

sf (3.129)

Tale tipo di minimizzazione viene scelta per assicurare che le derivate medie della pressione

risultanti siano indipendenti dalla costante arbitraria introdotta nella definizione di [20].

Si ottengono le seguenti relazioni

ppD

psD

ˆˆ~ (3.130.a)

ppD

D

ˆ~ (3.130.b)

in cui sono stati introdotti i seguenti parametri

hs

psD

pp

ˆˆ

)()ˆ(

ˆˆ22 (3.131)

dove h è valutato applicando al prodotto h , la stessa formula di integrazione usata per e .

Quando o tendono a zero, le espressioni (3.130) non diventano singolari. Questo può

accadere nel caso in cui sia che tendono a zero; nel secondo caso anche il p è nullo per

cui il problema non si pone.

La scelta dell’energia interna per unità di volume e della densità come variabili termodinamiche

indipendenti, permette di minimizzare il numero di variabili indipendenti che devono essere

mediate nella cella. In questo modo, lo stato medio di Roe è definito come

~,~,~

,~~ HuVM μ (3.132)

Bisogna ricordare che i valori di ~ e ~ in generale non coincidono con i valori che si otterrebbero

dall’integrazione esatta delle (3.124). Comunque l’integrazione approssimata seguita dalla fase di

61

proiezione permette ottenere una coppia di valori di ~ e ~ tali per cui la condizione (3.56) sia

soddisfatta. Si riporta infine l’espressione della matrice A~

uuu

uHuHuuVM

~~2~~~~

2

1~

~~~~~1~~~

2

1~~

100

~

22

22

μA (3.133)

e della velocità del suono

22 ~

2

1~~~~ uHa (3.134)

3.7 Distribuzione del residuo

Una volta determinato il bilancio dei flussi in modo conservativo su ogni cella è necessario

distribuire la fluttuazione ai nodi secondo uno schema di tipo upwind. La decomposizione

caratteristica permette di suddividere la fluttuazione riferita alle variabili conservative Φ in tre

componenti k che rappresentano il bilancio dei flussi delle onde caratteristiche [8]. Dunque si può

scrivere

3

1 2

1

2

1

2

1

2

1

2

1

2

1

2

1

2

1

2

1~~~

k

k

i

k

i

au

i

au

i

u

i

iiii

rRφRΦ (3.135)

in cui k rappresenta la fluttuazione associata alla esimak onda caratteristica e viene

determinata considerando l’equazione (3.11) utilizzata per l’equazione di convezione scalare

kkix

ix

kk

ix

ix

kk

iwdx

x

wdx

t

w

~11

2

1 (3.136)

62

Le velocità medie ~

corrispondono agli autovalori della matrice jacobiana mediata alla Roe.

La distribuzione viene fatta quindi suddividendo la fluttuazione associata ad ogni variabile

caratteristica in una parte diretta al nodo i e una diretta al nodo 1i nel seguente modo

WΛWΛφφφ

~~212121 iii

(3.137)

in cui

,,0diagΛ e ΛΛΛ . Dalla (3.137) si ricava il seguente schema di

discretizzazione per la (3.28):

WΛWΛWW

~~1

x

tn

i

n

i (3.138)

Nel caso di flusso subsonico la fluttuazione corrispondente alle prime due variabili caratteristiche

viene inviata al nodo 1i , mentre quella corrispondente alla terza viene inviata al nodo i , come

mostrato in figura (3.6).

Per un flusso supersonico invece tutte le componenti della fluttuazione vengono inviate al nodo

1i

Figura 3.6 Distribuzione della fluttuazione nel caso di flusso localmente subsonico [8]

Figura 3.7 Distribuzione della fluttuazione nel caso di flusso localmente supersonico [8]

63

Considerando l’equazione (3.27) è possibile suddividere la matrice Jacobiana in due componenti,

una positiva ed una negativa, in maniera analoga a quanto fatto per la matrice degli autovalori Λ .

LRΛA (3.139)

dove AAA .

Dalla (3.139) è possibile suddividere la fluttuazione espressa in funzione delle variabili

conservative nella componente diretta al nodo i e uno diretta al nodo 1i

UAUAΦΦΦ

~~

2

1

2

1

2

1iii

(3.140)

si ottiene quindi lo schema di integrazione temporale espresso in funzione delle variabili

conservative

UAUAUU

~~1

x

tn

i

n

i (3.141)

Per il caso di gas reale le componenti della matrice

A si esprimono come segue:

a

u

aa

2202011

02

ˆ

122

aa E

22

10ˆ13

a

u

aa

E

a

uHu

a

Ha

220

2

221

a

u

a

Ha

E

Eˆ2

ˆ

002222

a

u

a

H

a

Hua

E

E

2

ˆ02

ˆ

2322

1

64

a

u

aaua

2

023122

1

220

ˆ

32

a

u

aa E

2

122

ˆ02

2

ˆ33

a

u

a

ua

EE

Nel caso in cui si consideri la presenza del termine sorgente, le equazioni di Eulero si

decompongono come espresso dalla (3.46)

k

kk

k

sx

w

t

wˆ

(3.142)

La fluttuazione associata ad ogni variabile caratteristica viene quindi calcolata come segue

2

1

11

2

1

~ˆ

~ˆ

i

kkix

ix

kk

ix

ix

kk

ixswdxs

x

wdx

t

w (3.143)

in cui il termine s~ˆ rappresenta il termine sorgente medio nella cella, calcolato in corrispondenza

dello stato medio qJss μ~ˆ

~ˆ . Per il termine sorgente espresso dalla (3.48) si ottiene:

qa

qa

qa

E

E

E

~2

~~2

~

~~

~

~ˆ

ˆ

ˆ

2

ˆ

S

(3.144)

Una volta determinata la fluttuazione, lo schema di distribuzione è identico a quello mostrato per il

caso di flusso privo di termine sorgente.

65

Capitolo 4

Condizioni al contorno

66

4.1 Introduzione

La scelta di condizioni al contorno appropriate rappresenta uno dei problemi chiave per la

risoluzione numerica delle equazioni di Navier-Stokes o di Eulero. Affinché la soluzione delle

equazioni di conservazione esista e sia unica è necessario scegliere condizioni al contorno per cui il

problema sia ben posto; inoltre è necessario scegliere un’opportuna discretizzazione, tale da non

pregiudicare l’accuratezza e la stabilità del metodo di risoluzione utilizzato all’interno del dominio

di calcolo [6]. Generalmente vengono distinte due categorie di condizioni al contorno: le condizioni

al contorno fisiche e le condizioni al contorno numeriche. Le condizioni al contorno fisiche sono

legate alle informazioni che propagano dall’esterno verso l’interno del dominio di calcolo, in

particolare dipendono dalle variabili associate alle caratteristiche entranti nel dominio. Le restanti

condizioni al contorno vengono denominate numeriche, in quanto sono legate alle caratteristiche

uscenti dal dominio e quindi dipendono dalla soluzione calcolata all’interno del dominio. Tali

condizioni al contorno devono essere scelte in modo tale da essere consistenti con le proprietà

fisiche del flusso, e compatibili con il metodo numerico di risoluzione.

In figura 4.1 e 4.2 vengono mostrate le curve caratteristiche entranti ed uscenti, rispettivamente per

un flusso supersonico ed uno subsonico, in tabella 4.1 viene riportato invece il numero di variabili

da imporre per un flusso unidimensionale inviscido.

Le condizioni al contorno sono state implementate seguendo due approcci diversi, uno basato

sull’ampiezza delle onde caratteristiche e l’altro sulla definizione di un flusso correttivo nelle celle

al bordo.

Figura 4.1 Caratteristiche entranti ed uscenti dal dominio per un flusso supersonico unidimensionale [6]

67

4.2 Condizioni al contorno NSCBC

Il metodo basato sulle caratteristiche è il cosiddetto metodo NSCBC (Navier-Stokes Characteristic

Boundary Conditions), introdotto da Thompson [22] e sviluppato successivamente da Poinsot e

Lele [23]. L’idea base di tale approccio è quella di esprimere le condizioni al contorno fisiche in

termini di ampiezze delle onde entranti nel dominio di calcolo. Per assicurare che il problema sia

ben posto, le onde uscenti dal dominio non devono essere specificate, ma devono essere calcolate

risolvendo sul bordo le stesse equazioni di conservazione utilizzate all’interno del dominio di

calcolo.

La determinazione delle condizioni al contorno corrette può essere ottenuta mediante delle relazioni

di compatibilità tra le variabili scelte come condizioni al contorno fisiche e l’ampiezza delle onde

entranti. Riferendosi alle variabili primitive , u e p si ottengono, per un flusso unidimensionale

inviscido, le seguenti relazioni di compatibilità dette LODI (Local One-Dimensional Inviscid) [23]:

Inflow Outflow

Subsonico 2 1

Supersonico 3 0

Figura 4.2 Caratteristiche entranti ed uscenti dal dominio per un flusso subsonico unidimensionale [6]

Tabella 4.1 Numero di condizioni al contorno fisiche

68

02

112

-0 LLLct

(4.1)

02

1

-LLct

u

(4.2)

02

1

-LLt

p (4.3)

dove c è la velocità del suono e 0L , L e -L rappresentano le ampiezze delle onde caratteristiche

associate rispettivamente alle velocità caratteristiche u , cu e cu .

Le iL vengono date dalle seguenti espressioni:

x

wu

x

p

xau

12 0L (4.4)

x

wau

x

ua

x

pau

2

L (4.5)

x

wau

x

ua

x

pau

3

L (4.6)

Utilizzando le relazioni 4.1 - 4.3 è possibile determinare delle relazioni di compatibilità valide per

altre grandezze come la temperatura, l’entropia o l’entalpia totale [23]:

012

12

-0 LLL- a

T

t

T

(4.7)

0

1

1

0LTt

s

(4.8)

0112

1

1

1

-0 LLL- MMt

h

.

(4.9)

Note le relazioni di compatibilità è possibile esprimere le condizioni al contorno fisiche in termini

di ampiezze delle onde caratteristiche, ad esempio facendo riferimento all’equazione 4.3 si ottiene

che imporre la pressione costante in un nodo al bordo equivale a richiedere che -LL , oppure

dalla 4.8 si ottiene che imporre l’entropia costante in un nodo al bordo equivale a richiedere che

0L sia nullo in tale nodo. Considerando i vettori LLL ,,0L e il vettore delle variabili

caratteristiche 321 ,w,wwW si ottiene:

69

x

WΛL (4.10)

dove è la matrice degli autovalori.

Mediante la 4.10 è possibile determinare la fluttuazione nelle celle al bordo come:

ULΛRUAF iiii

ii

i

2

1

2

1

2

1

2

1

~~~~

xi

i

i

i

ii

ii

2

1

2

1

0

2

1

2

1

2

1

2

1

~~~

L

L

L

RWΛR

Di seguito vengono riportati alcuni esempi di implementazione.

4.2.1 Outlet subsonico

Nel caso di un outlet subsonico le onde sono tutte uscenti dal dominio tranne quella associata

all’autovalore au 3 , per cui è necessario specificare una sola condizione al contorno fisica.

In particolare imponendo la pressione statica all’outlet uguale a quella del flusso indisturbato p , si

ottiene un problema ben posto, per cui viene scelta la seguente espressione per l’ampiezza dell’onda

caratteristica entrante L [23]:

ppl

cL (4.12)

in cui l è una dimensione caratteristica del dominio di calcolo e è un parametro che va scelto

opportunamente. Per valori di elevati l’outlet risulta essere molto riflessivo, in particolare per

1 l’outlet risulta essere totalmente riflessivo; per valori tendenti a zero si ottiene invece un

outlet perfettamente non-riflessivo.

70

4.2.2 Inlet subsonico

Per un inlet subsonico le onde caratteristiche entranti sono individuate unicamente dagli autovalori

positivi, per cui è necessario specificare due condizioni al contorno. Affinché il problema risulti ben

posto bisogna fare particolare attenzione alla scelta delle variabili da specificare come condizioni al

contorno fisiche. In particolare si dimostra che la scelta di u e p come condizioni al contorno

fisiche dà luogo ad un problema mal posto [6]; inoltre è stato mostrato che un flusso subsonico

stazionario, con area uguale all’inlet e all’outlet, può convergere a soluzioni non uniche, nel caso in

cui venga specificata una stessa variabile all’inlet e all’outlet [6]. La scelta di ed u come

variabili da specificare all’ingresso, invece dà luogo ad un problema ben posto [23], per cui

l’ampiezza delle caratteristiche entranti 0L e L si può esprimere come:

0

20

al

cL (4.13)

0uual

c

L (4.14)

Molto spesso per un inlet subsonico si considerano fissate la temperatura totale e la pressione totale,

in tal caso le ampiezze delle onde caratteristiche vengono espresse come:

oo TTRl

a0

0

L (4.15)

oo ppl

a0

L (4.16)

4.2.3 Outlet ed inlet supersonico

Per un outlet supersonico le onde caratteristiche sono tutte uscenti per cui non è necessario

specificare nessuna condizione al contorno fisica.

Nel caso di inlet supersonico è necessario invece inferire l’ampiezza di tutte le onde caratteristiche

mediante le equazioni di compatibilità 4.1 - 4.3. Fissando le variabili primitive , u e p si

ottengono le seguenti relazioni:

71

0

20

cl

cL (4.17)

0uucl

c

L (4.18)

0ppl

c

L (4.19)

4.3 Metodo del flusso correttivo

Nelle celle al bordo il residuo corretto può essere ottenuto aggiungendo alla fluttuazione calcolata

mpcoΦ un termine correttivo

corrΦ , tale che le condizioni al contorno fisiche siano rispettate:

corrcompΦΦΦ

corrcomp

FFF (4.20)

Nel caso in cui le condizioni al contorno vengano espresse in termini delle variabili primitive

tpu,,V , è necessario esprimere il flusso nelle celle al bordo in funzione del gradiente del

vettore V , in maniera tale da garantire che la linearizzazione sia conservativa

ZZZ

VZ

V

FZZ

Z

FZ

Z

FFF ii

cellacella

i dxdx

ddx

dx

d

(4.21)

La matrice ZV non è lineare rispetto alle componenti del vettore Z , per cui non è possibile

determinare il gradiente medio di V nella cella semplicemente come iiii xxVV 11 [13], si

ottiene quindi:

cella i

i

icellaicella xdx

dx

d

xdx

dx

d

xdx

d ZZ

Z

VZ

Z

VVV 11. (4.22)

Sostituendo nell’equazione 1.22 l’espressione del gradiente di Z e della matrice ZZV si

esplicita l’espressione del gradiente di V come segue:

72

x

pp

uux

zzzzzz

zzzzz

zz

xdx

d

i

i

ii

ii

ii

i

iii

ii

i

icella

VV ˆ

~ˆ1

1

1

2

1

1

1

1

331221

13312

1

11

(4.23)

dove 1

ˆ ii e 21 2~

ii . Il gradiente delle variabili primitive assume la stessa

espressione sia nel caso di gas perfetto che nel caso di gas reale; si ottiene infatti che:

33122111

13312

1

11

21

1

2

1

zzzzzzzz

zzzzz

zz

xdx

d

iiii

ii

i

icella

V

ii

ii

ii

i

iiiiii

ii

ii

i

pp

uux

pp

uux

1

1

1

111

1

1~1

1

1

~1

(4.24)

Per ragioni di semplicità espositiva viene introdotto il vettore

iiii pu ,~

ˆ,

V

tale per cui

VV

ii . Le densità medie ~ e vengono considerate rispetto alla cella 21 per l’inlet, e

rispetto alla cella 21n per l’outlet.

Il flusso nelle celle al bordo può essere quindi espresso in funzione delle variabili primitive:

Inlet

bcompcompbcompcompb ΔΔ 00100011011

~~~~~VVQVQVVVVQVVQVQFFF

b

01

bcompcomp

0011

~VVQFF

73

Outlet

comp

n

b

n

comp

nn

comp

n

comp

n

b

nn

b

nnn

b

nn VVQVQVVVVQVVQVQFFF

~~~~~111

comp

n

b

n

comp

nn VVQFF

~

Il termine correttivo è quindi pari a:

comp

n

b

n

corr

n VVQF

~

(4.25)

bcompcorr

001

~VVQF

. (4.26)

I termini che presentano l’apice b si riferiscono alle condizioni al contorno fisiche imposte; la

matrice Q~

è pari allo Jacobiano del flusso rispetto al vettore delle variabili primitive tpu,,V

calcolato in corrispondenza dello stato medio Z

ZV

FZQQ

~

ZV

ZZ

Z

FZ

V

FQ

~ (4.27)

per un gas ideale si ottiene:

12

12

0

2

23

uu

uuHu

u

V

FQ . (4.28)

Nel caso di una generica equazione di stato si ottiene invece:

12

1

0

2 uu

uuHuu

u

m

m

V

FQ (4.29)

74

in particolare considerando la generalizzazione di Vinokur e Montagné la matrice (4.29) si riduce

alla seguente:

12

1

2

0

2

22

uu

uuHuu

u

V

FQ

(4.30)

Bisogna sottolineare che nel caso di gas reale la matrice Q~

viene definita come segue funzione sia

Z che di ~ e ~ .

Di seguito vengono riportati alcuni esempi di implementazione per un outlet subsonico con

pressione costante, ed un inlet subsonico con densità e velocità costanti.

comp

n

comp

n

b

n

corr

n

pp

0

0~~QVVQF

(4.31)

0

~ˆ~~

00

00

001

bcalc

bcalc

bcompcorr uu

QVVQF . (4.32)

75

Capitolo 5

Fluctuation Splitting 2D per gas reale

76

5.1 Introduzione

Gli schemi di tipo FV centrati sulle celle, vengono estesi al caso multidimensionale considerando

una serie di problemi di Riemann localmente unidimensionali, lungo le direzioni normali alle facce

delle celle. Tale tipo di generalizzazione presenta notevoli svantaggi a causa della natura

unidimensionale del modello fisico utilizzato.

L’estensione multidimensionale dell’approccio FS d’altro canto permette di ottenere una

discretizzazione genuinamente multidimensionale delle equazioni di Eulero multidimensionali [11].

La naturale estensione della rappresentazione lineare continua assunta per il caso unidimensionale

(figura 3.2) è quella di considerare una griglia costituita da elementi di forma triangolare su cui le

variabili dipendenti varino in maniera lineare e continuo.

In figura 5.2 viene mostrata la cella elementare di forma triangolare e le normali che individuano le

facce. Le normali alle facce possono essere calcolate come segue

ykjkji xxyy een x

con kji ,, permutazione circolare

dove il vettore iii yx ,x individua il nodo esimoi . Le normali alle facce sono tali che

03

1

i

in .

Figura 5.1 Rappresentazione funzionale lineare delle variabili dipendenti nel caso bidimensionale

77

La superficie dell’elemento triangolare può essere calcolato come:

3

14

1

i

iiTS nx (5.1)

Considerando l’equazione scalare di convezione lineare 2 D

0

w

t

wλ (5.2)

è possibile scrivere il gradiente della variabile dipendente come:

3

12

1

i

ii

TTSTS

wS

dwwdSw nn (5.3)

In tale caso la fluttuazione sull’elemento triangolare T viene calcolata come

wSwdS T

TS

T λλ (5.4)

Sostituendo la (5.3) nella (5.4) si ottiene per la fluttuazione la seguente espressione

3

12

1

i

ii

T w nλ .

Figura 5.2 Cella elementare e normali alle facce

78

Facendo riferimento alle equazioni di Eulero 2D

0

F

t

U

(5.5)

0

yxt

GFU (5.6)

l’espressione della fluttuazione assume la seguente forma:

TSTSTS

t

T dSdSdS yx GFUΦ F (5.7)

Applicando il teorema della divergenza, la fluttuazione può essere ottenuta calcolando il seguente

integrale circuitale

TSTS

T dxdyd GFnΦ F (5.8)

Nota la fluttuazione della cella T è necessario determinarne la distribuzione ai nodi come mostrato

in figura 5.3 in maniera tale che T

k

T

j

T

i

TΦΦΦΦ .

Figura 5.3 Distribuzione della fluttuazione in 2D

79

Nel presente lavoro verrà affrontata unicamente la generalizzazione al caso 2D della

decomposizione caratteristica e della linearizzazione conservativa delle equazioni di Eulero scritte

per gas reale, verrà quindi tralasciata la trattazione relativa alla distribuzione della fluttuazione ai

nodi.

5.2 Decomposizione caratteristica delle equazioni di Eulero 2D per gas reale

Per poter determinare la struttura caratteristica delle equazioni 5.5 è necessario determinare la forma

quasi lineare seguente

0

A U

U

t (5.9)

in cui la matrice jacobiana è pari a

U

FA (5.10)

Usando la notazione vettoriale si può scrivere che il vettore delle variabili conservative e il vettore

dei flussi sono pari rispettivamente a uU ,, E e a pH uuuu ,,F [11].

Considerando la (5.5) è possibile scomporre la matrice jacobiana A nelle sue componenti

cartesiane

yx BeAe A (5.11)

dove A e B sono le matrici jacobiane delle componenti cartesiane del flusso F e G .

L’iperbolicità del sistema di equazioni (5.9) implica che la matrice

nAn A (5.12)

abbia autovalori reali ed un set completo di autovettori destri e sinistri per ogni versore

yyyx nn eeeen xx sincos , scelto in maniera arbitraria.

80

yx nn BAAn (5.13)

Per poter esprimere la matrice nA per un gas reale è necessario generalizzare al caso

bidimensionale la funzione introdotta precedentemente.

nmEp ,,ˆ, U (5.14)

Differenziando la (5.13) si ottiene

dndmEdddp nmE ˆ

ˆ (5.15)

in cui le derivate parziali della pressione rispetto alle variabili conservative sono pari a:

Enm

p

ˆ,,

,

nmE

E

p

,,

ˆ ˆ

,

nE

mm

p

,ˆ,

,

mE

nn

p

,ˆ,

(5.16)

Considerando il vettore quantità di moto nm,m si può scrivere in modo più compatto la (5.15)

come mm dEdddpE

Πˆˆ , in cui nm ,mΠ .

La matrice nA si può esprimere quindi come

Iunnnun

n

n

A

m

mn

nEn

n

t

nEn

t

uu

uHuuH

Π

Πt

ˆ

ˆ1

00

(5.17)

Dove nu nu e I è la matrice identità 2x2. Gli autovalori della matrice nA sono pari a

,,, 003,2,1 dove

a 00 , nu (5.18)

e possono essere raccolti nella matrice diagonale seguente

,,, 00diagnΛ (5.19)

81

La velocità del suono presente nell’espressione degli autovalori viene calcolata come:

E

Ha ˆ2 uu (5.20)

Agli autovalori (5.18) viene associato un set di autovettori destri rappresentati dalle colonne della

matrice seguente

nu

nu

su

nunusuRn

aa

a

H

a

HaH

aa

E

ˆ

2

01

(5.21)

e un set di autovettori sinistri rappresentati dalle righe della matrice seguente

ttEEE

ttEEE

t

tEEE

aaH

aa

aaH

aa

aaH

a

nuuunu

nuuunu

ssu

uuu

Ln

ˆˆˆ

ˆˆˆ

2

ˆ

2

ˆ

2

ˆ

2

1

22

1

2

1

22

1

0

(5.22)

Il vettore s forma con n una base ortonormale, per cui yyxy nn eeees xx cossin .

Noti gli autovettori destri e sinistri è possibile fattorizzare la matrice nA in maniera simile al caso

unidimensionale ottenendo

nnnn LΛRA (5.23)

Tale fattorizzazione permette di scomporre la matrice jacobiana nelle sue componenti positive e

negative lungo n analogamente a quanto esposto per il caso unidimensionale

82

nnnnLΛRA

(5.24)

Le componenti della matrice

nA vengono riportate in Appendice B.

Considerando delle particolari variabili, note come variabili simmetrizzanti, è possibile rendere

simmetrica la matrice jacobiana, ottenendo notevoli vantaggi dal punto di vista computazionale.

u

Ua

p

pa

2

1

~ (5.25)

Tali variabili sono legate alle variabili conservative mediante le seguenti trasformazioni

I0u

U

U m

m

1

ΠΠ

ΠΠ1

~ˆ

22

ˆ

2

aaa

aaa

Eρ

Eρ

t

t

Π

Π

(5.26)

Iuu

uu

0

U

U m

a

Ha

a

t

E

ρ

t

ˆΠ

Π

1

~

tΠ

(5.27)

Le equazioni di Eulero possono essere scritte in funzione delle variabili simmetrizzanti

premoltiplicando l’equazione (5.9) per la matrice UU ~

0~~

~~

U

UU

U

U

UA A

tt (5.28)

dove yx eBeA~~~

A e

83

u

ua

au

u

000

00

00

000

~

~~

U

UA

U

UA ,

va

v

av

v

00

000

00

000

~

~~

U

UB

U

UB (5.29)

Si ottiene quindi la seguente espressione

Inun0

nnu

0nu

nAn

a

a t

t

0

0~~

A (5.30)

La matrice nA~

presenta gli stessi autovalori della matrice nA , mentre presenta i seguenti

autovettori:

nns

Rn

0

1100

0001~

e

t

t

t

t

n

n

s

0

Ln

21

21

21

21

0

0

00

01

~ (5.31)

Le equazioni di Eulero scritte in funzione delle variabili simmetrizzanti (5.28) hanno il vantaggio di

presentare la medesima espressione sia per gas perfetto che per gas reale.

5.3 Linearizzazione conservativa delle equazioni di Eulero 2D per gas reale

Per poter discretizzare la (5.7) è necessario definire sulla cella T i gradienti medi dei vettori dei

flussi e del vettore delle variabili conservative [13]

yxyx GFGF

TSTTST

dSS

dSS

11FF (5.32)

TST

dSS

UU1

(5.33)

84

È possibile determinare il gradiente dei flussi considerando uno stato medio U~

tale per cui

UUGFΦ yx ~

AF TT

TS

T SSdS (5.34)

Considerando la seguente generalizzazione del parametro vettore di Roe al caso bidimensionale

v

u

H

z

z

z

z 1

4

3

2

1

Z (5.35)

si ottiene che per gas perfetto F , G ed U sono funzioni omogenee di grado due delle componenti

del parametro vettore. Assumendo che Z vari linearmente nella cella si ottengono le seguenti

espressioni per i gradienti medi (5.32) e (5.33)

TST

dSS

UU1

ZZZ

UZ

Z

UZ

Z

U

TSTTST

dSS

dSS

11 (5.36)

TST

dSS

FF1

ZZZ

ZZ

ZZ

FFF

TSTTST

dSS

dSS

11

(5.37)

dove ZZ è costante nella cella e Z è lo stato medio che assicura che la linearizzazione sia

conservativa. Poiché si assume che il parametro vari in modo lineare sulla cella lo stato medio Z è

pari semplicemente alla media algebrica dei valori assunti nei nodi.

85

3

3

3

3

3

332211

332211

332211

321

321

vvv

uuu

HHH

ZZZZ (5.38)

Per gas perfetto si ottiene quindi che fluttuazione nella cella è pari a

ZZZ

UZGFΦ yx

AF TT

TS

T SSdS (5.39)

Nel caso di gas reale la pressione non è una funzione di secondo grado delle componenti di Z per

cui per definire lo stato medio è necessario considerare dei parametri aggiuntivi oltre a Z .

È possibile decomporre il flusso F in due termini [21], uno quadratico rispetto ai termini di Z e

l’altro funzione della pressione:

ttp pH I0,0uuuu ,,,* FFF (5.40)

e in maniera simile è possibile decomporre il vettore delle variabili conservative

ttp pH 0,0,,0,,* uUUU (5.41)

Dalla funzione (5.14) si ottiene che la matrice jacobiana mediata alla Roe si può esprimere come

funzione del parametro vettore medio e di alcuni parametri aggiuntivi legati alle derivate parziali

della pressione rispetto alle variabili conservative

nmE~

,~

,~

,~

, ˆZAA (5.40)

Considerando la linearizzazione proposta da Vinokur e Montagné per il caso unidimensionale è

possibile legare le derivate della pressione rispetto alle variabili conservative alle derivate rispetto

alle variabili termodinamiche

86

mmU

2

1ˆ, Ep (5.41)

uu 21 , um , vn ,

E~ (5.42)

Quindi la matrice (5.40) diventa ~,~,ZAA .

Dalla (3.116), tenendo conto che pH 2uu si può scrivere

dddp (5.43)

dpdHdddp

2

uu (5.44)

21

uu dHdddp (5.45)

I termini di destra della (5.47) sono tutti funzioni di secondo grado di Z

221

2

4

2

3

21

2

1

zzdzzdzddp (5.46)

Sostituendo la (5.40) nella (5.34) si ottiene

TTT S

p

S

*

S

dSdSdS FFF

TTT S

p*

S

p

S

*

dSdSdS FF

FF

ZZZ

ZZ

(5.47)

Poiché

p

p

p

0

0

00

00

F (5.48)

87

segue che ttp p ,0F . Poiché la pressione non è una funzione omogenea di ordine due delle

componenti del parametro vettore non è possibile determinare il gradiente medio esatto

semplicemente come ZZZ

ZZ

pdS

p

SpdS

Sp

TT STST

11, analogamente al caso di

gas perfetto. Per gas reale si ottiene

ZZZ

ZZ

~,~,

1 pdS

p

Sp

TST

(5.49)

Considerando la (5.45), il gradiente della pressione medio si può scrivere come

2

~~~1

1 uu

Hp (5.50)

e quindi

4433211211

~2~~1

1zzzzzzzzzzp

(5.51)

Nel caso unidimensionale la (5.50) si riduce ad una sola equazione, per cui non è possibile

determinare in maniera univoca le derivate medie ~ e ~ . Nel caso bidimensionale invece la (5.50)

rappresenta un sistema di due equazioni nelle due incognite ~ e ~ . È possibile quindi ottenere lo

stato medio risolvendo semplicemente il seguente sistema lineare

bM

~

~

(5.52)

dove la matrice M

e il vettore b sono pari a:

y

p

y

zz

y

zz

y

zz

y

zz

y

zz

x

p

x

zz

x

zz

x

zz

x

zz

x

zz

44

3

31

22

11

1

44

3

31

22

11

1

2

2

M

,

y

px

p

pb (5.53)

88

bM1

~

~

(5.54)

È importante sottolineare che pur essendo determinate in maniera univoca, i parametri ~ e ~

ottenuti risolvendo il sistema (5.52) non rappresentano le derivate della pressione calcolate in

corrispondenza dello stato medio Z ( Z ~ e Z ~ ).

Il gradiente medio della pressione che appare nel sistema (5.52) deve essere determinato scegliendo

una formula di quadratura appropriata per l’integrale circuitale seguente4

TT STST

dpS

pdSS

p n11

(5.55)

Nel caso in cui la matrice M sia singolare, ossia quando 0det M , il sistema (5.52) non è

risolvibile. I diversi casi in cui non è possibile determinare ~ e ~ risolvendo il sistema vengono

riportati nelle figure seguenti

4 Nel caso di gas perfetto l’integrale (5.55) può essere risolto in maniera esatta considerando una formula di quadratura

di ordine superiore al secondo (ad esempio la formula di Simpson). Nel caso unidimensionale invece la (5.55) si riduce

semplicemente alla differenza p tra gli stati relativi ai due nodi che individuano la cella, per cui non è necessaria

nessuna formula di quadratura.

Figura 5.4 Rappresentazione grafica del sistema (5.52)

89

Figura 5.5 Caso a

Figura 5.6 Caso b

90

Figura 5.7 Caso c

Figura 5.8 Caso d

91

Nei casi a, b e c è impossibile determinare lo stato medio in maniera tale che vengano soddisfatte

entrambe le equazioni della (3.52), mentre il caso d è riconducibile al caso unidimensionale, in cui

l’equazione (3.123) viene sostituita dalla seguente

iii xxx

p

~~ (5.56)

Nel caso tridimensionale il sistema dato dalla (5.49) è costituito da tre equazioni in due incognite,

per cui non è risolvibile in nessun caso. È evidente quindi che l’utilizzo di una linearizzazione di

tipo quasi-Jacobiana, come quella proposta da Vinokur e Montagnè, non è estendibile in modo

sufficientemente robusto alle equazioni di Eulero multidimensionali scritte per gas reale in

equilibrio termochimico.

Nel caso di gas reale in non equilibrio chimico, è possibile utilizzare lo stesso approccio seguito per

il caso unidimensionale a causa della presenza delle derivate s~ della pressione rispetto alle densità

delle varie specie chimiche [20]-[21].

Utilizzando una linearizzazione in forma strettamente Jacobiana si potrebbe ottenere un metodo

estendibile al caso multidimensionale in maniera più robusta. In tale caso il numero di incognite

sarebbe pari al numero di equazioni sia per il caso bidimensionale che per quello tridimensionale.

92

Capitolo 6

Risultati delle simulazioni

93

6.1 Introduzione

Lo schema numerico presentato nei capitoli precedenti è stato utilizzato per la risoluzione delle

equazioni di Eulero unidimensionali con termine di produzione di calore interno q costante in

tutto il dominio. La soluzione stazionaria di tali equazione corrisponde al ben noto flusso di

Rayleigh. Lo studio del flusso con termine sorgente è stato esteso anche al caso in cui la sorgente di

calore interna non è costante ed è dovuta ad un riscaldamento di tipo ohmico, dovuto ad un

passaggio di corrente elettrica trasversale al flusso. Di seguito vengono presentati brevemente i due

modelli ed i risultati ottenuti.

6.2 Il Flusso di Rayleigh

Il flusso di Rayleigh rappresenta una buona approssimazione per sistemi in cui lo scambio termico

con l’esterno non è trascurabile come ad esempio gli scambiatori di calore, o per sistemi in cui è

presente una sorgente di calore interna come ad esempio le camere di combustione.

In tali sistemi la variazione di entropia dovuta all’attrito è trascurabile rispetto a quella dovuta allo

scambio termico o alla sorgente di calore interna [24]

fh dsds

per cui si ottiene:

fh dsdsds (6.1)

hdsds . (6.2)

Figura 6.1 Schematizzazione per il flusso di Rayleigh [24]

94

Considerando le ipotesi seguenti:

- Flusso unidimensionale

- Flusso stazionario

- Attrito trascurabile

- Assenza di lavoro esterno

- Condotto a sezione costante

è possibile scrivere le equazioni di conservazione nel seguente modo [25]:

- l’ equazione di continuità:

2211 uu (6.3)

- l’ equazione di conservazione della quantità di moto:

2

222

2

111 upup (6.4)

- l’equazione di conservazione dell’energia:

22

2

22

2

11

uhq

uh

(6.5)

dove i è la densità, iu la velocità, ip la pressione ed ih l’entalpia nella sezione i-esima.

Si ottiene quindi che la variazione di entalpia totale 0h tra due sezioni è pari alla quantità di calore

fornito o sottratto q

0012 hhq

(6.6)

Per un gas caloricamente perfetto è possibile scrivere l’entalpia totale come prodotto del calore

specifico a pressione costante pc per la temperatura totale 0T

)( 00

12TTcq p

(6.7)

95

Per cui una fornitura o una sottrazione di calore danno luogo rispettivamente ad un aumento o ad

una riduzione della temperatura totale.

Per un gas perfetto inoltre vale la seguente relazione:

pMMau 22

(6.8)

dove M è il numero di Mach e a è la velocità del suono, dalla quale è possibile ricavare i rapporti

tra le grandezze termofluidodinamiche tra due sezioni diverse in funzione dei numeri di Mach:

2

2

2

1

1

2

1

1

M

M

p

p

(6.9)

2

1

2

2

2

2

2

1

1

2

1

1

M

M

M

M

T

T

(6.10)

2

2

1

2

1

2

2

1

2

1

1

M

M

M

M

(6.11)

2

1

2

22

1

2

2

2

2

2

1

0

1

0

2

2

11

2

11

1

1

M

M

M

M

M

M

T

T

. (6.12)

Lungo il flusso la portata massica per unità di sezione u e la funzione impulso 2up restano

costanti, per cui una volta fissate tali quantità il flusso è univocamente determinato.

Utilizzando le forme differenziali delle equazioni (6.3) e (6.4), e dell’equazione di stato per gas

perfetto si ottengono le seguenti relazioni [25]:

0u

dud

(6.13)

0p

udu

p

dp

(6.14)

96

d

T

dT

p

dp

(6.15)

mediante l’equazione (6.8) è possibile riscrivere l’equazione (6.14) come segue

02 u

duM

p

dp

(6.16)

per cui si ha

T

dT

M

d

u

du21

1

.

(6.17)

Sostituendo tali espressioni nella relazione seguente:

dR

T

dTcds v

(6.18)

si ottiene

T

dTc

M

Mds v2

2

1

1

.

(6.19)

L’equazione precedente permette di rappresentare nel piano sT la curva rappresentata in figura

6.2, detta curva di Rayleigh. Il ramo superiore della curva rappresenta le condizioni di flusso

supersonico, mentre il ramo inferiore rappresenta le condizioni subsoniche. Nel punto a, poiché

0ds il numero di Mach è pari ad 1, mentre nel punto b, in cui 0dT , si ha che 1M .

97

Dalla seguente relazione di Gibbs

0dhududhdp

dhTds

0dhTdsq

(6.20)

si ottiene che in corrispondenza della condizione sonica (punto a) il differenziale dell’entalpia totale

è nullo, 00 dh , ossia in tale punto sia l’entropia che l’entalpia totale presentano un massimo. Per

cui il punto a rappresenta un punto limite per il flusso oltre il quale non è possibile fornire ulteriore

calore. Risulta evidente che uno scambio termico positivo tende ad accelerare un flusso subsonico e

viceversa a rallentare un flusso supersonico, fino al raggiungimento della condizione limite di flusso

sonico. Nel caso in cui un flusso subsonico raggiunga la condizione sonica, un’ ulteriore fornitura di

calore determinerà una riduzione della portata, e una variazione delle condizioni a monte del flusso;

in tal caso si dirà che il flusso è termicamente strozzato [24]. Se il flusso invece è supersonico, una

fornitura di calore ulteriore determinerà la formazione di un urto all’interno del condotto, come

mostrato in figura 6.3. Si dimostra che nota la quantità di calore fornito e la pressione in uscita, non

è possibile determinare in maniera univoca la posizione dell’urto lungo l’asse x.

Figura 6.2 Curva di Rayleigh [25]

98

La quantità di calore massima che può essere fornita può essere calcolata conoscendo le condizioni

del flusso entrante ed imponendo all’uscita un numero di Mach pari ad uno. Si ottiene quindi:

1

2

11

0

1

0

*2

11maxT

TMTcq p

(6.21)

1

121

1

2

11

2

1

2

1

22

12

11MM

MMTc p

2

1

22

11max

12

1

M

MTcq p

(6.22)

Nella tabella seguente viene riassunto il comportamento del flusso :

0q

Flusso subsonico Flusso supersonico

u aumenta diminuisce

p diminuisce aumenta

diminuisce aumenta

T

aumenta per 1M

aumenta diminuisce per 1M

M aumenta diminuisce

0T aumenta aumenta

Tabella 6.1 Comportamento delle variabili nel caso in cui venga fornito del calore al

flusso. Nel caso in cui venga sottratto del calore il comportamento è speculare

Figura 6.3 Flusso di Rayleigh in presenza di un urto

99

Anche la sottrazione di calore presenta un limite, in quanto in linea di principio è possibile sottrarre

al più un quantità di calore pari all’entalpia totale posseduta inizialmente dal fluido.

Contrariamente a quanto ci si potrebbe aspettare una fornitura di calore non corrisponde

necessariamente ad un aumento della temperatura del gas, infatti nel caso in cui il flusso sia

subsonico e il numero di Mach sia superiore a 1 , la temperatura subisce una diminuzione

nonostante la fornitura di gas. Ciò è spiegabile tenendo conto che in corrispondenza di 1M

l’entalpia presenta un massimo, per cui il calore fornito può dare luogo solo ad un aumento

dell’energia cinetica, e quindi della velocità.

6.3 Simulazioni numeriche per un flusso di Rayleigh 1D

Nel presente capitolo vengono presentati ed analizzati i risultati ottenuti mediante la simulazione di

un flusso di Rayleigh unidimensionale. I test sono stati effettuati sia per gas perfetto che per gas

reale allo stato di plasma, in condizioni di flusso subsonico e di flusso supersonico.

6.3.1 Risultati ottenuti per gas perfetto

Nel caso di gas perfetto è possibile determinare la soluzione analitica mediante le equazioni (6.9) -

(6.12) a patto di conoscere il numero di Mach all’ingresso e all’uscita. In figura 6.4 vengono

mostrate le condizioni al contorno e le condizioni iniziali considerate nel primo test, per un flusso

supersonico (6.4.a) e per un flusso subsonico (6.4.b).

(a) Flusso supersonico

100

La quantità di calore fornita può essere espressa in funzione del numero di Mach all’ingresso e

all’uscita del flusso, mediante le equazioni (6.7) e (6.12). Nella tabella 6.2 vengono riportati i valori

utilizzati per le simulazioni.

1

0

1

0

20

1T

TTcq p

(6.23)

Nelle figure 6.5, 6.6 e 6.7 vengono mostrati gli andamenti lungo x, delle variabili primitive, si può

notare come la soluzione numerica ottenuta con il metodo di Roe sia in totale accordo con la

soluzione analitica.

Inlet Outlet

Subsonico 2.01 M 4.02 M

Supersonico 5.11 M 2.12 M

Flusso subsonico Flusso supersonico

Figura 6.4 Condizioni iniziali e condizioni al contorno per il gas perfetto

(b) Flusso subsonico

Figura 6.5 Pressione lungo x

Tabella 6.2 Numeri di Mach utilizzati peri il primo test

101

6.3.2 Risultati ottenuti per gas reale

Per il caso di gas reale è stato effettuato un primo test, in cui la temperatura iniziale è stata posta

pari a 300K. In tali condizioni il comportamento del gas reale non presenta sostanziali differenze

rispetto a quello di un gas termicamente e caloricamente perfetto; inoltre è stata considerata la

stessa quantità di calore fornita nel caso di gas perfetto. In figura 6.8 vengono mostrate le

condizioni al contorno ed iniziali utilizzate per il flusso supersonico.

Figura 6.6 Velocità lungo x

Flusso supersonico Flusso subsonico

Figura 6.7 Densità lungo x

Flusso supersonico Flusso subsonico

102

Nelle figure seguenti vengono mostrati i risultati del primo test considerando solo il caso di flusso

supersonico.

Figura 6.8 Condizioni iniziali e al contorno per il gas reale

Figura 6.9 Andamento della pressione lungo x

Figura 6.10 Andamento della velocità lungo x

103

Come ci si aspettava, alle basse temperature, il comportamento del gas reale non si discosta molto

da quello del gas perfetto. Se il calore fornito al flusso non è elevato e la temperatura iniziale è

prossima a quella standard, è possibile considerare il gas caloricamente e termicamente perfetto con

buona approssimazione.

Il secondo test è stato effettuato considerando una temperatura iniziale del gas pari a 3000k, ed un

termine sorgente elevato. In figura 6.12 vengono mostrate le condizioni al contorno ed iniziali

utilizzate per il secondo test.

Figura 6.11 Andamento della densità lungo x

Figura 6.12 Condizioni iniziali e al contorno per il gas reale

(a) Flusso supersonico

(b) Flusso subsonico

104

In tali condizioni non sono più trascurabili gli effetti delle vibrazioni atomiche e la composizione

chimica del gas non può essere più considerata costante con la temperatura e la pressione, inoltre

per valori della temperatura superiore ai 2000 K non sono più trascurabili nemmeno gli effetti delle

dissociazioni molecolari.

Le figure seguenti mostrano i risultati ottenuti per il gas reale confrontati con quelli ottenuti per il

gas perfetto.

Figura 6.13 Andamento della pressione lungo x

Figura 6.14 Andamento della velocità lungo x

Flusso supersonico

Flusso supersonico Flusso subsonico

Flusso subsonico

105

Si può notare come a parità di calore fornito la temperatura del gas reale sia inferiore a quella del

gas perfetto, ciò è dovuto essenzialmente all’energia vibrazionale, non trascurabile alle alte

temperature. Dalla figura 6.14 si osserva che il flusso supersonico ed il flusso subsonico, subiscono

rispettivamente una decelerazione ed un’accelerazione minore rispetto al caso ideale. La minore

variazione della velocità, rispetto al caso di gas perfetto, si riflette su una minore variazione della

densità, sia per il caso di flusso subsonico che per il caso di flusso supersonico; ciò è spiegabile

tenendo in conto della costanza della portata. La costanza della funzione impulso 2up , infine

permette di spiegare la minore variazione della pressione rispetto al caso ideale.

Per quanto riguarda le condizioni al contorno, sono state implementate seguendo entrambe gli

approcci mostrati nel capitolo 4. Per il metodo NSCBC sono state effettuate delle prove per

Figura 6.15 Andamento della densità lungo x

Figura 6.16 Andamento della temperatura lungo x

Flusso subsonico Flusso supersonico

Flusso subsonico Flusso supersonico

106

analizzare l’effeto del parametro presente nelle (4.12)-(4.19). In particolare è stato studiato

l’effetto di tale parametro nel caso di un outlet subsonico sulla convergenza della soluzione. In

figura 6.17 viene riportato l’andamento del residuo massimo al variare del numero di iterazioni per

diversi valori di .

Dalla figura precedente si osserva che 0.5 risulta essere un valore ottimale per il parametro .

Nella figura seguente viene mostrato l’andamento del termine ppn

durante il transitorio per

diversi valori del parametro . Anche in questo caso si osserva che la soluzione a regime viene

raggiunta prima nel caso in cui 5.0 .

Figura 6.17 Storia di convergenza per gas reale subsonico al variare del parametro (0.2, 0.5, 0.9)

107

6.4 Modello con sorgente variabile

Il modello con sorgente interna di tipo ohmico si basa sulla schematizzazione riportata in figura

6.17. In tale modello il campo elettrico agisce in maniera trasversale al flusso e viene generato da

due elettrodi piani e paralleli connessi ad un generatore G mediante un circuito avente resistenza

cR [26]. La legge di Ohm per il circuito mostrato in figura 6.17 si può scrivere come

tIRVtV cgdis (6.24)

dove disV è la differenza di potenziale tra i due elettrodi e gV è la differenza di potenziale ai capi del

generatore G e I è la corrente che attraversa il circuito.

Figura 6.18 Andamento della differenza tra pressione nel nodo n e pressione imposta al variare di (0.2, 0.5, 0.9)

108

È possibile determinare la corrente considerando la densità di corrente J che attraversa il fluido tra

gli elettrodi funzione soltanto dell’asse x e del tempo t

eS

JdStI (6.25)

dove dS è l’elemento infinitesimo di superficie ed eS rappresenta la superficie trasversale

dell’elettrodo come mostrato in figura 6.18.

Figura 6.17 Schema del circuito

Figura 6.18 Dimensioni degli elettrodi

109

Scrivendo la densità di corrente J come prodotto della conducibilità elettrica del gas per il campo

elettrico si E

EJ (6.26)

e sostituendo dxldS e nella (6.25) si ottiene

b

a

x

x

edxElI (6.27)

È necessario sottolineare che l’equazione (6.27) è valida solo sotto l’ipotesi che il campo elettrico

sia costante e si diretto lungo la normale all’asse x ; sotto tali i potesi il campo elettrico risulta

essere espresso come

e

dis

d

VE (6.28)

Sostituendo la (6.27) e la(6.28) nella legge di Ohm (6.24) si ottiene

b

a

x

x

ecge dxElRVEd (6.29)

Poiché il campo elettrico viene assunto costante lungo x , può essere portato fuori dal segno di

integrale presente nella (6.29)

b

a

x

x

ecge dxElRVEd (6.30)

Essendo stato assunto il fluido in condizioni di equilibrio termochimico locale, segue che le sue

proprietà termodinamiche e di trasporto sono funzioni di due sole variabili come ad esempio la

pressione e la temperatura, o la densità e l’energia interna. La conducibilità elettrica in un punto x

e

110

txTtxptx ,,,, (6.31)

o

txtxtx ,,,, (6.32)

La determinazione della conducibilità elettrica per un gas ionizzato presenta notevoli difficoltà a

causa della sua dipendenza non monotona dalle variabili termodinamiche. Nel presente lavoro è

stato fatto riferimento alla procedura di calcolo presentata da D’Angola et al. [3]:

7

1

0 loglogj

jj TaT (6.32)

in cui e j rappresentano le funzioni (1.14) e (1.13), mentre i coefficienti ja vengono riportati

in appendice A.

Una volta nota la distribuzione della temperatura e della pressione nel dominio è possibile calcolare

la conducibilità elettrica e quindi il campo elettrico tra gli elettrodi.

b

a

x

xe

e

c

eg

dxtxd

lR

dVtE

,1

(6.33)

Se si suppone che nelle singole celle in cui viene discretizzato il dominio di calcolo, la conducibilità

elettrica media sia pari al valore assunto in corrispondenza dello stato medio U~~ , si ottiene

per la (6.33) la seguente discretizzazione

1

1

~11

N

i

ii

e

c

egn

xd

R

dVE

(6.34)

Tale campo elettrico dà luogo ad un riscaldamento Q all’interno del flusso per effetto Joule

tEtxJtxQJ ,, (6.35)

Considerando la (6.26) si ottiene

111

2,, tEtxtxQ (6.36)

Tale termine sorgente può essere considerato costante nelle varie celle ad ogni istante di tempo, si

ottiene quindi la seguente discretizzazione

2~ n

i

n

i EQ (6.37)

Il termine (6.36) rappresenta il termine sorgente per l’equazione dell’energia, si ottiene quindi

QHt

E

u

(6.38)

6.5 Simulazioni per il modello con sorgente variabile

Nel presente paragrafo vengono presentati ed analizzati i risultati ottenuti considerando il modello

presentato nel paragrafo precedente. In figura 6.19 vengono mostrati i valori dei parametri del

circuito e le condizioni al contorno ed iniziali del flusso, considerati per il primo test. La distanza

tra gli elettrodi è stata posta pari a mde

310

Figura 6.19 Test 1

112

La soluzione delle equazioni di Eulero raggiunge una condizione stazionaria; nelle figure seguenti

viene riportato l’andamento delle variabili termofluidodinamiche lungo x .

Figura 6.20 Test 1: Andamento della pressione lungo x

Figura 6.21 Test 1:Andamento della temperatura lungo x

113

Nelle regioni non comprese tra gli elettrodi il flusso è stato considerato adiabatico, per cui a regime

le variazioni lungo x sono concentrate soltanto tra

ax e bx . Poiché il flusso è subsonico il

riscaldamento per effetto Joule determina un aumento del numero di Mach e una diminuzione della

pressione.

Il campo elettrico è stato considerato costante lungo l’asse x e dipendente unicamente dal tempo.

Nelle figure 6.23 e 6.24 viene riportato l’andamento nel tempo del campo elettrico tra gli elettrodi e

della corrente che attraversa il circuito.

Figura 6.22 Test 1:Andamento del numero di Mach lungo

x

Figura 6.23 Test 1:Andamento del campo elettrico nel tempo

114

Il circuito raggiunge una condizione stazionaria dopo un transitorio pseudoperiodico.

È stato effettuato un secondo test considerando la temperatura e la pressione iniziali maggiori

rispetto a quelle utilizzate nel test 1.

In figura 6.25 vengono mostrate le condizioni iniziali e al contorno per il test 2.

Di seguito vengo riportate gli andamenti lungo x delle variabili termofluidodinamiche

Figura 6.24 Test 1:Andamento della corrente nel tempo

Figura 6.25 Test 2

115

Figura 6.26 Test 2: Andamento della pressione lungo x

Figura 6.27 Test 1: Andamento della temperatura lungo x

Figura 6.28 Test 2: Andamento del numero di Mach lungo x

116

In figura 6.29 e 6.30 vengono riportati glia andamenti nel tempo del campo elettrico e della corrente

elettrica

Figura 6.29 Test 2: Andamento del campo elettrico nel tempo

Figura 6.30 Test 2: Andamento del corrente nel tempo

117

Conclusioni

Nel presente lavoro è stato mostrato ed analizzato un metodo numerico per lo studio dei gas

ionizzati basato su un nuovo modello termochimico dell’aria e su una generalizzazione del metodo

Fluctuation Splitting al caso di gas reale. La procedura di calcolo delle proprietà termodinamiche

dell’aria mediante i fit proposti si è dimostrato particolarmente accurata ed efficiente nel caso

unidimensionale.

La generalizzazione del metodo Fluctuation Splitting è stata costruita a partire da una formulazione

del tutto generale proposta in letteratura [4] per il metodo di Roe, valida nell’ambito dell’approccio

ai volumi finiti. È stato mostrato come la definizione dello stato medio non è univoca nel caso di

gas reale e di come venga influenzata dal tipo di equazione di stato considerata. In particolare è

stato mostrato come l’equazione di stato proposta da Vinokur e Montagnè [15], in cui le variabili

indipendenti considerate sono la densità e l’energia interna per unità di volume, permetta di

utilizzare un minor numero di parametri per la definizione dello stato medio.

Per quanto riguarda le condizioni al contorno è stato mostrato come l’approccio basato

sull’ampiezza delle onde caratteristiche possa essere applicato al caso di gas reale in equilibrio ed

inoltre è stato proposto un metodo innovativo basato sull’utilizzo di un flusso correttivo nelle celle

al bordo.

Le simulazioni effettuate hanno mostrato come il modello di Rayleigh per gas reale si discosti

notevolmente da quello per gas perfetto nel caso in cui la temperatura e la pressione siano lontane

da quelle standard. In particolare è stato mostrato come nel caso di gas ionizzato l’utilizzo

dell’equazione di stato dei gas perfetti produca delle sovrastime per quanto riguarda le variazioni di

temperatura dovute all’apporto di calore.

Infine la linearizzazione di Vinokur e Montagnè [15], è stata generalizzata per il metodo

Fluctuation Splitting multidimensionale. È stato mostrato come l’utilizzo di tale linearizzazione, e

in generale delle linearizzazioni in forma quasi-Jacobiana, non è estendibile al caso

tridimensionale, è auspicabile quindi che in futuro venga sviluppata una linearizzazione di tipo

strettamente Jacobiana per i gas reali in equilibrio.

118

Ringraziamenti

Per la realizzazione di questo lavoro vorrei ringraziare il Prof. A. Bonfiglioli e il Dr. A. D’Angola

per la loro disponibilità e per i loro fondamentali suggerimenti, ringrazio inoltre il Dr. G. Colonna

del CNR di Bari. Un ringraziamento speciale va alla mia famiglia: i miei genitori Giuseppe e

Agnese, mio fratello Donato e mia sorella Mariateresa, per avermi sostenuto in questi anni. Infine

vorrei ringraziare tutti i miei amici e i compagni di corso.

119

Bibliografia

[1] Goldston R.; Rutherford P. H., Introduction to Plasma Physics, Institute of

Physics Publishing, Philadelphia, p. 1-3, 1995

[2] Choudhuri A. R., The Physics of Fluids and Plasmas, Cambridge University

Press, 1998

[3] D’Angola A.; Colonna G.; Gorse C.; Capitelli M., Thermodynamic and Transport

Properties in Equilibrium Air Plasmas in a Wide Pressure and Temperature Range ,

The European Physical Journal, vol. 46, 129-150, 2008

[4] Mottura L.; Vigevano L.; Zaccanti M., An evaluation of Roe’s scheme

generalizations for equilibrium real gas, Journal of Computational Physics, vol. 138,

354-399, 1997

[5] Tannehill J. C.; Mugge P. H., Improved Curve Fits for the Thermodynamic

Properties of Equilibrium Air Suitable for Numerical Computation using Time-

Dependent or Shock-Capturing Method, Nasa CR-2470, October 1974

[6] Hirsch C., Numerical Computation of Internal and External Flows Vol II ,

Second Edition, John Wiley & Sons, New York, 2007

[7] Tannehill J. C.; Anderson D. A.; Pletcher R. H. Computational Fluid Mechanics

and Heat Transfer, Second Edition, Taylor & Francis, New York, 1997

[8] Bonfiglioli A., Unstructured Grid Method in CFD with Emphasis on Fluctuation

Splitting Schemes, Dipartimento di Ingegneria e Fisica dell’Ambiente, Università degli

studi della Basilicata, 2009

[9] Toro E. F., Riemann Solvers and Numerical Methods for Fluid Dynamics ,

Springer, Berlin, 1999

120

[10] Roe P. L., Approximate Riemann solvers, parameter vectors, and difference

schemes, Journal of Computational Physics, vol. 43, 357-372, 1981

[11] Bonfiglioli A.; Deconinck H., Multidimensional Upwind Scheme for the 3D Euler

Equations on Unstructured Tetrahedral Meshes, Notes on Numerical Fluid Mechanics,

vol. 57, 141-185, Vieweg, 1996

[12] Paillère H., Multidimensional Upwind Distribution Scheme for the Euler and

Navier-Stokes Equations on Unstructured Meshes, PhD thesis, Universitè Libre de

Bruxelles, Bruxelles, 1995

[13] Deconinck H.; Roe P. L.; Struijs R., A Multidimensional Generalization of Roe’s

Flux Difference Splitter for the Euler Equations, Computer & Fluids, vol. 22, 215-222,

1993

[14] Guardone A.; Vigevano L., Roe Linearization for the van der Waals Gas , Journal

of Computational Physics, vol. 175, 50-78, 2002

[15] Vinokur M.; Montagné L., Generalized flux-vector splitting and Roe average for

an equilibrium real gas , Journal of Computational Physics, vol. 89, 276-300, 1990

[16] Glaister P., An approximate linearized Riemann solver for the Euler equations for

real gases , Journal of Computational Physics, vol. 74, 382-408, 1988

[17] Abgrall R., An extension of Roe’s upwind scheme to algebraic equilibrium real gas

models, Computer & Fluids, vol. 19, 171-182, 1991

[18] Liou M. S.; van Leer B.; Shuen J. S., Splitting of inviscid fluxes for real gases ,

Journal of Computational Physics, vol. 87, 1-24, 1990

[19] Toumi I., A weak formulation of Roe’s approximate Riemann solver , Journal of

Computational Physics, vol. 102, 360-373, 1992

[20] Liu Y.; Vinokur M., Upwind algorithms for general thermo-chemical

nonequilibrium flows, AIAA Paper 89-201, 1989

[21] Degrez G.; van der Weide E., Upwind residual distribution schemes for chemical

nonequilibrium flows, Collection of Technical Papers, vol. 2, 978-987, AIAA, 1999

121

[22] Thompson K.W., Time-dependent boundary conditions for hyperbolic systems ,

Journal of Computational Physics, vol. 68(1), 1-24, 1987

[23] Poinsot T.J.; Lele S.K., Boundary conditions for direct simulations of compressible

viscous flows, Journal of Computational Physics, vol. 101, 104-129, 1992

[24] Zucker R. D.; Biblarz O., Fundamental of Gas Dynamics, Second Edition, John Wiley

& Sons, New York, 2002

[25] Napolitano M., Flusso stazionario unidimensionale non isentropico, Gasdinamica,

Capitolo 3, Poliba

[26] Colonna G.; Capitelli M., Boltzmann and master equations for

Magnetohydrodynamics in weakly ionized gases, Journal of Thermophysics and Heat

Transfer, vol. 22 (3), 414-423, 2008

122

Appendice A

Proprietà termodinamiche dell’aria

123

A.1 Coefficienti per il calcolo della massa molare media.

function y=dens(mT,mp)

mp=log(mp);

R=8.314472; % J mol-1 K-1

catm_pascal=1.01325e5;

sa1=([0.028811 ]);

sa2=(fliplr([ -5.452539e+000 -2.762076e-002 -3.327630e-003 -2.453118e-004 -6.332107e-006

]));

sc2=(fliplr([ 8.170734e+000 5.708244e-002 1.293374e-003 ]));

sd2=(fliplr([ 6.380594e+000 1.046470e-001 8.553860e-004 -1.572857e-004 ]));

sa3=(fliplr([ -4.595514e+000 1.328152e-002 9.294096e-004 -8.243998e-005 -9.490079e-006

]));

sc3=(fliplr([ 8.805680e+000 5.468057e-002 1.121881e-003 ]));

sd3=(fliplr([ 7.080690e+000 1.142540e-001 6.869247e-004 -2.257365e-004 ]));

sa4=(fliplr([ -4.971377e+000 -1.668833e-002 -2.409638e-003 -2.840529e-004 -2.934495e-005

]));

sc4=(fliplr([ 9.525862e+000 6.639994e-002 7.836529e-004 -2.447910e-004 -2.415297e-005

]));

sd4=(fliplr([ 7.888211e+000 9.954169e-002 -1.327510e-004 -2.926560e-004 -4.717532e-005

]));

sa5=(fliplr([ -6.720756e+000 7.203127e-002 6.766486e-003 -1.019894e-003 9.196578e-005

]));

124

sc5=(fliplr([ 1.055726e+001 8.397717e-003 9.849480e-004 3.539965e-004 -4.236150e-005

]));

sd5=(fliplr([ 8.707609e+000 3.713173e-002 -1.670186e-002 -5.094908e-004 4.248200e-004

]));

sa6=(fliplr([ -6.218117e+000 -7.145834e-002 6.529894e-004 1.599394e-003 1.981881e-005

]));

sc6=(fliplr([ 1.020784e+001 2.553473e-002 -3.549988e-003 ]));

sd6=(fliplr([ 8.422438e+000 1.125955e-001 -3.204629e-003 -1.655103e-003 -2.051312e-004

]));

sa7=(fliplr([ -6.611171e+000 8.990124e-002 -5.418532e-003 ]));

sc7=(fliplr([ 1.096136e+001 2.887564e-002 -3.621097e-004 ]));

sd7=(fliplr([ 9.253817e+000 1.341329e-002 -6.004835e-003 1.860800e-003 -1.229602e-004

]));

conc=sa1-exp(polyval(sa2,mp)).*sgm(mT,exp(polyval(sc2,mp)),exp(polyval(sd2,mp)))...

-exp(polyval(sa3,mp)).*sgm(mT,exp(polyval(sc3,mp)),exp(polyval(sd3,mp)))...

-exp(polyval(sa4,mp)).*sgm(mT,exp(polyval(sc4,mp)),exp(polyval(sd4,mp)))...

-exp(polyval(sa5,mp)).*sgm(mT,exp(polyval(sc5,mp)),exp(polyval(sd5,mp)))...

-exp(polyval(sa6,mp)).*sgm(mT,exp(polyval(sc6,mp)),exp(polyval(sd6,mp)))...

-exp(polyval(sa7,mp)).*sgm(mT,exp(polyval(sc7,mp)),exp(polyval(sd7,mp)));

mpp=exp(mp)*1.01325e5;

y=conc/R.*mpp./mT;

125

A.2 Coefficienti per il calcolo dell’entalpia specifica.

function y=ei(mT,mp)

mp=log(mp);

sa0=fliplr([ 2.350912e-001 -1.120236e-003 -2.508755e-005 ]);

sa1=fliplr([ 1.542966e-005 6.556647e-007 ]);

sa2=(fliplr([ 6.587335e+000 -6.112145e-002 -9.108114e-003 -9.569561e-004 -1.128838e-004 -

8.757988e-006 ]));

sc2=(fliplr([ 8.164839e+000 5.283021e-002 4.741812e-004 -1.276598e-004 -9.877950e-006

]));

sd2=(fliplr([ 6.513247e+000 1.040239e-001 -8.104042e-004 -2.991537e-004 4.348437e-005

6.258153e-006 ]));

sa3=(fliplr([ 8.740885e+000 3.050736e-003 1.599171e-003 -2.859059e-004 -5.371695e-005

]));

sc3=(fliplr([ 8.856133e+000 5.964702e-002 1.745638e-003 2.343688e-005 -3.102821e-006

]));

sd3=(fliplr([ 6.981907e+000 1.119408e-001 4.185626e-003 -2.499247e-004 -5.209456e-005

]));

sa4=(fliplr([ 1.014496e+001 -1.833015e-002 -4.265166e-003 -8.321612e-004 -6.481810e-005

]));

126

sc4=(fliplr([ 9.593196e+000 7.089945e-002 1.640521e-003 -1.055407e-004 -1.510653e-005

]));

sd4=(fliplr([ 7.910995e+000 1.006930e-001 -1.608832e-003 -2.526731e-004 ]));

sa5=(fliplr([ 1.082665e+001 -4.777223e-002 -4.682547e-003 ]));

sc5=(fliplr([ 1.030572e+001 6.607308e-002 1.512694e-003 -5.009486e-005 -5.192881e-006

1.116840e-006 ]));

sd5=(fliplr([ 8.320951e+000 7.474585e-002 1.789257e-003 5.273341e-004 3.755570e-005

3.425485e-006 ]));

sa6=(fliplr([ 1.145937e+001 5.122940e-004 -8.805300e-003 -1.193042e-003 ]));

sc6=(fliplr([ 1.076031e+001 6.404003e-002 9.621465e-004 -1.883920e-005 ]));

sd6=(fliplr([ 8.846750e+000 1.307197e-001 -2.943134e-004 -6.425060e-004 ]));

sa7=(fliplr([ 1.172458e+001 -5.461477e-002 3.413385e-003 7.407737e-004 -1.644220e-004

-1.878043e-005 ]));

sc7=(fliplr([ 1.109244e+001 6.026294e-002 1.125935e-003 -2.170126e-005 -3.141895e-006

]));

sd7=(fliplr([ 8.942747e+000 8.687938e-002 1.554323e-002 3.584506e-005 -2.447405e-004

]));

sa8=(fliplr([ -1.011841e+001 -2.295953e+001 -1.220667e+001 -3.504472e+000 -4.373233e-001

1.127311e-002 6.598926e-003 -2.119755e-004 -1.369506e-004 -8.311253e-006]));

sc8=(fliplr([ 1.314544e+001 2.079129e+000 9.992304e-001 2.223931e-001 1.963954e-002

-1.622592e-004 -1.094608e-005 2.304744e-005 1.817656e-006 ]));

sd8=(fliplr([ -1.743314e+000 -1.807206e+001 -1.393980e+001 -5.232064e+000 -7.607736e-001

8.529592e-002 4.967101e-002 7.733746e-003 5.507513e-004 1.527569e-005]));

y0=polyval(sa0,mp).*mT+polyval(sa1,mp).*mT.^2 ...

+exp(polyval(sa2,mp)).*sgm(mT,exp(polyval(sc2,mp)),exp(polyval(sd2,mp)))...

+exp(polyval(sa3,mp)).*sgm(mT,exp(polyval(sc3,mp)),exp(polyval(sd3,mp)))...

+exp(polyval(sa4,mp)).*sgm(mT,exp(polyval(sc4,mp)),exp(polyval(sd4,mp)))...

+exp(polyval(sa5,mp)).*sgm(mT,exp(polyval(sc5,mp)),exp(polyval(sd5,mp)))...

+exp(polyval(sa6,mp)).*sgm(mT,exp(polyval(sc6,mp)),exp(polyval(sd6,mp)))...

+exp(polyval(sa7,mp)).*sgm(mT,exp(polyval(sc7,mp)),exp(polyval(sd7,mp)))...

+exp(polyval(sa8,mp)).*sgm(mT,exp(polyval(sc8,mp)),exp(polyval(sd8,mp)));

y0=y0*4.1868*1e3;

mp=exp(mp);

y=y0-mp*1.01325e5/dens(mT,mp);

127

A.3 Coefficienti per il calcolo della conducibilità elettrica.

function y=sig(mT,mp)

mp=log(mp);

R=8.314472; % J mol-1 K-1

catm_pascal=1.01325e5;

fa =(fliplr([ 1.635045e+000 4.450390e-002 -5.928863e-004 ]));

fc =(fliplr([ 5.748398e+000 6.411299e-002 ]));

fd =(fliplr([ 1.786355e+000 -1.212690e-002 -2.375673e-004 ]));

fw =(fliplr([ 1.419925e+000 -3.875497e-002 ]));

sa1 =(fliplr([ 4.493934e-002 -9.063708e-003 -2.489500e-003 ]));

sc1 =(fliplr([ 8.930803e+000 5.718843e-002 1.093759e-003 ]));

sd1 =(fliplr([ 7.014976e+000 7.625175e-002 3.011941e-004 ]));

sa2=(fliplr([ 1.593153e+000 4.137850e-002 1.430491e-002 -4.403957e-007 ]));

sc2=(fliplr([ 8.576847e+000 1.004174e-001 7.406762e-003 -1.095186e-003 ]));

sd2=(fliplr([ 9.113182e+000 -8.202725e-002 6.299430e-003 9.099828e-004 ]));

sa3=(fliplr([ 2.627897e-001 2.917558e-003 3.036205e-003 -1.926651e-004 -2.917018e-

005 ]));

sc3=(fliplr([ 1.023493e+001 6.651575e-002 1.090308e-003 -6.576415e-005 4.715318e-007

]));

128

sd3=(fliplr([ 8.039563e+000 1.435966e-001 8.862611e-003 -3.478227e-004 -3.614711e-

005 ]));

sa4=(fliplr([ 1.707216e-001 2.035164e-002 1.809127e-003 -9.630175e-005 1.781249e-005

]));

sc4=(fliplr([ 1.072380e+001 5.671452e-002 1.468210e-004 2.608196e-005 6.511719e-

006 ]));

sd4=(fliplr([ 8.556977e+000 2.227207e-001 -2.773160e-003 -1.684219e-003 1.878188e-

004 ]));

sa5=(fliplr([ 2.480007e-001 2.217818e-002 9.094614e-004 ]));

sc5=(fliplr([ 1.106431e+001 5.578774e-002 6.639655e-004 ]));

sd5=(fliplr([ 9.309043e+000 1.366510e-001 -2.599317e-003 ]));

sa6=(fliplr([ 3.636707e+000 -1.436268e-001 -2.934926e-003 ]));

sc6=(fliplr([ 1.023203e+001 8.703300e-002 5.007602e-003 ]));

sd6=(fliplr([ 1.130562e+001 -2.184155e-002 -1.865543e-004 ]));

sc7=(fliplr([ 2.946755e+001 -4.289010e+000 -3.224136e-001 9.371814e-002 ]));

sd7=(fliplr([ 2.430324e+001 -2.653523e+000 -3.309222e-001 4.769061e-002 ]));

conc=exp(...

exp(polyval(fa,mp))-exp(polyval(fc,mp)).*exp(-

(log(mT)./exp(polyval(fd,mp))).^exp(polyval(fw,mp)))+...

exp(polyval(sa1,mp)).*sgm(mT,exp(polyval(sc1,mp)),exp(polyval(sd1,mp)))...

+(polyval(sa2,mp)).*sgm(mT,exp(polyval(sc2,mp)),exp(polyval(sd2,mp)))...

-(polyval(sa3,mp)).*sgm(mT,exp(polyval(sc3,mp)),exp(polyval(sd3,mp)))...

-(polyval(sa4,mp)).*sgm(mT,exp(polyval(sc4,mp)),exp(polyval(sd4,mp)))...

-(polyval(sa5,mp)).*sgm(mT,exp(polyval(sc5,mp)),exp(polyval(sd5,mp)))...

+(polyval(sa6,mp)).*sgm(mT,exp(polyval(sc6,mp)),exp(polyval(sd6,mp)))...

-(polyval(sa1,mp)+polyval(sa2,mp)+polyval(sa6,mp)-polyval(sa3,mp)-polyval(sa4,mp)-

polyval(sa5,mp)).*sgm(mT,exp(polyval(sc7,mp)),exp(polyval(sd7,mp)))...);

y=conc;

129

A.4 Coefficienti per il calcolo della pressione [5] .

load tabella.txt

%T=100:100:25000;

%p=([1.0e-2 2.0e-2 5.0e-2 1.0e-1 2.0e-1 5.0e-1 1.0e0 2.0e0 5.0e0]);

Nrho=80;

),( pp

130

Ne=100;

%rho_g=linspace(1e-2,17.5548,Nrho);

rho_g=exp(linspace(log(1e-2),log(18),Nrho));

ep_g=linspace(2e3,3e6,Ne);

%ep_g=exp(linspace(log(8.83e3),log(2.98e6),Ne));

[ep,rr]=meshgrid(ep_g,rho_g);

Y=log10(rr/1.292);

Z=log10(ep./rr./78408.4);

%YL=[-Inf -7 -4.5 -0.5 +Inf];

YL=[-7 -4.5 -0.5 +Inf];

ZL=[-Inf .65 1.5 1.54 1.68 2.2 2.22 2.46 2.9 3.05 3.38 +Inf];

%matrice dei coefficienti

a1=zeros(Nrho,Ne);

a2=zeros(Nrho,Ne);

a3=zeros(Nrho,Ne);

a4=zeros(Nrho,Ne);

a5=zeros(Nrho,Ne);

a6=zeros(Nrho,Ne);

a7=zeros(Nrho,Ne);

a8=zeros(Nrho,Ne);

a9=zeros(Nrho,Ne);

a10=zeros(Nrho,Ne);

a11=zeros(Nrho,Ne);

a12=zeros(Nrho,Ne);

a13=zeros(Nrho,Ne);

a14=zeros(Nrho,Ne);

a15=zeros(Nrho,Ne);

a16=zeros(Nrho,Ne);

%matr=zeros(length(ZL)-1,length(YL)-1);

matr=[10 5 1

11 6 2

12 6 2

12 7 2

12 7 3

13 7 3

13 8 3

13 8 4

13 9 4

14 9 4

15 9 4];

for jy=1:length(YL)-1

for jz=1:length(ZL)-1

%c1=((YL(jy)<Y<YL(jy+1))&(ZL(jz)<Z<ZL(jz+1)));

c1=(Y>YL(jy)&Y<YL(jy+1)&Z>ZL(jz)&Z<ZL(jz+1));

131

a1(c1)=tabella(matr(jz,jy),1);

a2(c1)=tabella(matr(jz,jy),2);

a3(c1)=tabella(matr(jz,jy),3);

a4(c1)=tabella(matr(jz,jy),4);

a5(c1)=tabella(matr(jz,jy),5);

a6(c1)=tabella(matr(jz,jy),6);

a7(c1)=tabella(matr(jz,jy),7);

a8(c1)=tabella(matr(jz,jy),8);

a9(c1)=tabella(matr(jz,jy),9);

a10(c1)=tabella(matr(jz,jy),10);

a11(c1)=tabella(matr(jz,jy),11);

a12(c1)=tabella(matr(jz,jy),12);

a13(c1)=tabella(matr(jz,jy),13);

a14(c1)=tabella(matr(jz,jy),14);

a15(c1)=tabella(matr(jz,jy),15);

a16(c1)=tabella(matr(jz,jy),16);

end

end

gamma=a1+a2.*Y+a3.*Z+a4.*Y.*Z+a5.*Y.^2+a6.*Z.^2+a7.*Y.*Z.^2+a8.*Z.^3+...

(a9+a10.*Y+a11.*Z+a12.*Y.*Z)./(1+exp((a13+a14.*Y).*(Z+a15.*Y+a16)));

p=ep.*(gamma-1);

[kappa,chi] = GRADIENT(p,ep_g,rho_g);

%catm_pascal=1.01325e5;

%p=p/catm_pascal;

save densita.txt rho_g -ASCII

save energia.txt ep_g -ASCII

save pressione.txt p -ASCII

save chi.txt chi -ASCII

save kappa.txt kappa -ASCII

132

Appendice B

Trasformazioni Jacobiane

133

B.1 Variabili Conservative

Considerando il vettore delle variabili conservative tE uU ,, si ottiene

15

14

13

543112

1

000

000

000

1111

12

1

1

00002

zz

zz

zz

zΠ

Πz

Π

Πz

Π

Πz

ΠΠzz

Π

z

E

E

E

E

E

E

EE

qq Πz

Z

U

11

5

11

4

11

3

2

1

5

2

1

4

2

1

3

11

2

1

1

1000

2

01

002

001

02

1

2

1

00002

1

zz

z

zz

z

zz

z

Πz

zΠ

z

zΠ

z

z

z

Π

z

zΠ

z

z

EEEE

ρ

U

Z

dove ttlnmwvu ,,,, q , tt

zzzwvu 543 ,,,, qz e tlnm ΠΠΠ ,,qΠ .

Considerando le componenti cartesiane del flusso inviscido

uw

uv

pu

uH

u

2F ,

uw

pv

uv

vH

v

2

G ,

pw

vw

uw

wH

w

2

H

si ottengono le seguenti matrici jacobiane

134

uwuw

uvuv

ΠΠuΠΠuΠ

uΠuΠuΠHΠuHΠu

lnmEρ

lnmEρ

00

00

2

1

00100

2

U

FA

vwvw

ΠvΠΠΠvΠ

uvvu

vΠvΠHvΠΠvHΠv

lnmEρ

lnmEρ

00

2

00

1

01000

2U

GB

wΠΠΠΠwΠ

vwwv

uwwu

wΠHwΠwΠΠwHΠw

lnmEρ

lnmEρ

2

00

00

1

10000

2

U

HC

35

34

543121

23

13

000

000

111

2

12

1

1

000

000

zz

zz

zΠ

Πz

Π

Πz

Π

Πz

Π

ΠΠzΠz

Π

zz

zz

E

E

E

E

E

E

E

EE

E

qq zΠZ

F

45

543121

34

24

14

000

11

2

112

1

1

000

000

000

zz

zΠ

Πz

Π

Πz

Π

Πz

Π

ΠΠzΠz

Π

zz

zz

zz

E

E

E

E

E

E

E

EE

E

qq zΠZ

G

135

543121

45

35

25

15

1

2

1112

1

1

000

000

000

000

zΠ

Πz

Π

Πz

Π

Πz

Π

ΠΠzΠz

Π

zz

zz

zz

zz

E

E

E

E

E

E

E

EE

E

qq zΠ

Z

H

B.2 Variabili Primitive

Considerando il vettore delle variabili primitive tp uV ,, si ottiene:

11

5

11

4

11

3

543121

1

1000

01

00

001

0

11112

1

1

00002

zz

z

zz

z

zz

z

zΠ

Πz

Π

Πz

Π

Πz

Π

ΠΠzΠz

Π

z

E

E

E

E

E

E

E

EEρ

E

qq Πz

Z

V

12

1

5

12

1

4

12

1

3

543

11

2

5

2

4

2

3

1

2

1

1

0002

1

0002

1

0002

1

11

2

1

00002

1

zz

z

zz

z

zz

z

zzzzΠ

Π

Π

Π

z

zzz

z

z

z

z

E

E

E

ρ

V

Z

136

B.3 Variabili Simmetrizzanti

Considerando il vettore delle variabili simmetrizzanti definito come tapap uU ,/,/~ 2

si ottiene

1

2

5

1

2

4

1

2

3

543121

2

5

2

4

2

3

2

12121

1000

01

00

001

0

11112

1

11

11112

1

12

~

1

1

1

zz

z

zz

z

zz

z

a

z

Π

Π

a

z

Π

Π

a

z

Π

Π

a

z

Π

ΠΠzΠz

Πa

a

z

Π

Π

a

z

Π

Π

a

z

Π

Π

a

z

Π

ΠΠzΠz

Πaz

E

E

E

E

E

E

E

EEρ

E

E

E

E

E

E

E

E

EEρ

E

qq

qq

Πz

Πz

Z

U

1

5

2

1

5

14

2

1

4

1

3

2

1

3

543

1

121

1

2

2

1

1

0022

1

0022

1

0022

1

2

11

22

1

00022

1

~

1

za

z

z

z

za

z

z

z

za

z

z

z

zzzΠz

ΠzzΠa

zΠ

Π

z

z

z

a

z

z

EE

E

ρ qqqq zzzz

U

Z

1000

01

00

001

0

1

~

22222

w

v

ua

Π

a

Π

a

Π

a

Π

a

Πa

Π

a

Π

a

Π

a

Π

a

Π

lnmEρ

lnmEρ

U

U

137

00

00

00

0001

~

wa

w

va

v

ua

u

wvuHaΠ

wΠvΠuΠΠa

E

lnmρ

U

U

B.4 Componenti della matrice

nA

333231

23n

Aaa

aA 2221

131211

aa

aaa

a

u

aa

2202011

02

ˆ

122

aa E

n

ua13

22

10ˆ

aa E

a

uHu

a

Ha n

n

220

2

221

a

u

a

Ha n

E

Eˆ2

ˆ

002222

ununa23

a

u

a

H

a

Hu n

E

En ˆ02

ˆ

22

unn

ua31

a

u

au

a

n

n

2202

138

n

ua32

220

ˆ

aa

E

Inuunuu

nnA33

0

ˆ

02

ˆ

2

1

2

aaa

EE

139

Appendice C

Codice numerico in Fortran

140

Program main

program main

!

! S1D Rayleigh using 1D fluctuation splitting

!

implicit real*8 (a-h,o-z)

integer ndof, nx

parameter (ndof=3)

real *8 mdc

parameter(nmax=1024)

parameter(nx=64)

! nodal stuff

dimension x(0:nx),du(ndof,0:nx),mdc(0:nx),dx(nx), z(ndof,0:nx)

! cell based stuff

real*8 phi(ndof,nx),al(ndof,nx),dw(ndof,nx),xc(nx)

real*8 p1,d1,u1,p2

real*8 cfl

real*8 df(3),sumdf(3)

common /abc/df

common /def/sumdf

data df,sumdf/6*0.d0/

logical lstop,lflag

character*1 answ

integer iok,niter,i,j

include 'aria.inc'

! initialize data and re-define grid

!

open(1,file='pressione.txt')

do j=1,Nr

read(1,*) (p(j,i),i=1,Ne)

enddo

close(1)

open(1,file='chi.txt')

do j=1,Nr

read(1,*) (chi(j,i),i=1,Ne)

enddo

close(1)

open(1,file='kappa.txt')

do j=1,Nr

read(1,*) (kappa(j,i),i=1,Ne)

enddo

close(1)

open(1,file='densita.txt')

read(1,*) (rho(i),i=1,Nr)

close(1)

open(1,file='energia.txt')

read(1,*) (ei(i),i=1,Ne)

close(1)

it = 0

call init(z,x,xc,dx,mdc,ndof,nx,tmax,p1,d1,u1,p2)

call plot(z,x,ndof,0,nx,it)

!

141

7 write(6,*)' Preconditioning (Y/N) '

read(5,FMT="(A1)")answ

if( answ .EQ. "Y" .OR. answ .EQ. "y" )then

lflag = .TRUE.

elseif( answ .EQ. "N" .OR. answ .EQ. "n" )then

lflag = .FALSE.

else

goto 7

endif

write(6,*)' Enter CFL [0<CFL<=1] & # of ITS '

read(5,FMT=*)cfl,niter

!

time = 0.d0

lstop = .FALSE.

do 1000 it = 1,niter

!

do 1 i = 0,nx

do 1 j = 1,ndof

du(j,i) = 0.d0

1 continue

do 3 j = 1,ndof

sumdf(j) = 0.d0

3 continue

!

! loop over cells

!

dt = 1.e38

do i = 1,nx

if(lflag)then

call he(z(1,i-1),z(1,i),dx(i),dt,du(1,i-

1),du(1,i),phi(1,i),al(1,i),dw(1,i),i)

else

call splitw(z(1,i-1),z(1,i),dx(i),dt,du(1,i-

1),du(1,i),phi(1,i),al(1,i),dw(1,i),i)

endif

enddo

!

dt = cfl * dt

if((time+dt).GE.tmax)then

dt = tmax-time

lstop = .TRUE.

endif

time = time + dt

call bcs(z,du,ndof,nx,dt,p1,d1,u1,p2) ! applica condizioni al contorno

call parm2cons(z,ndof,nx) ! trasforma Z in U

call update(z,du,ndof,mdc,time,dt,nx,iok) ! aggiorna U

call cons2parm(z,ndof,nx) ! trasforma U in Z

if(lstop.OR.(iok.EQ.0))goto 2000

!

! compute error

! call fd(u,du,x,dt,nx,phi)

! write(6,*)it,phi

!

1000 continue

2000 continue

write(6,*)time,tmax

! call w(phi,xc,ndof,nx,"fluct.dat")

! call w(dw,xc,ndof,nx,"strength.dat")

! call w(al,xc,ndof,nx,"speed.dat")

write(6,*)'Dumping solution at iteration ',it

call plot(z,x,ndof,0,nx,it)

call plot2(z,x,ndof,0,nx,it)

call check(z,mdc,ndof,nx)

142

write(22,*)(df(j),j=1,ndof)

stop

9900 format(i3,6(2X,e12.7))

end

Subroutine init

subroutine init(z,x,xc,dx,mdc,ndof,nx,timeout,p1,d1,u1,p2)

include 'paramt.h'

integer nx,ndof,je

include 'aria.inc'

real*8 x(0:nx),z(ndof,0:nx),xc(nx),dx(nx),mdc(0:nx)

real*8 zz(3),flx(3),df(3),heat(8),dom(3)

integer ix

real*8 TIMEOUT,qe,dotq,p1,t1,u1,d1,e1,p2

real*8 DD,Z1L,Z2L,Z3L,HELP,DOMLEN,DIAPH

real*8 v1,v2,checksum(3),eps1

double precision peps(Ne)

c OPEN(UNIT=1,FILE='exact.ini',STATUS='UNKNOWN')

timeout = 10000.d0

toler = 1.d-8

c open(33,file='initial.dat')

OPEN(UNIT=1,FILE='dom.txt')

READ(1,*)(dom(i),i=1,3)

close(1)

DOMLEN=dom(1)

DIAPH=dom(2)

c NX=dom(3)

c

OPEN(UNIT=2,FILE='heat.txt')

READ(2,*)(heat(i),i=1,8)

close(2)

p2=heat(8)

dotq=heat(2)

p1=heat(3)

t1=heat(4)

u1=heat(5)

d1=heat(6)

e1=heat(7)

dd = domlen/real(nx)

z1l = sqrt(d1)

z3l = z1l * u1

c z2l = gogm1 * p1/d1 + 0.5d0*(u1*u1)

call Pepsilon(peps,d1)

je=0

do j=1,Ne

if( peps(j) .LT. p1 )then

je=je+1

143

endif

enddo

eps1=ei(je)+(ei(je+1)-ei(je))/(peps(je+1)-peps(je))*(p1-peps(je))

ereps=abs(e1*d1-eps1)/eps1*100

write(6,*)' e1 = ',e1*d1,' eps1= ',eps1

write(6,*)' error =', ereps,'% '

c

c z2l = eps1+p1+0.5d0*d1*(u1*u1)

z2l = e1*d1+p1+0.5d0*d1*(u1*u1)

z2l = z2l/z1l

c

c write(6,*)' Z(0) = ',z1l,z2l,z3l

zz(1)=z1l

zz(2)=z2l

zz(3)=z3l

call flux(flx(1),zz(1),ndof)

c do j =1,ndof

c df(j) = flx(j,2)-flx(j,1)

c enddo

c write(6,*)(flx(j,1),j=1,ndof)

c write(6,*)(flx(j,2),j=1,ndof)

c write(6,*)(df(j),j=1,ndof)

do 1 ix = 0,nx

x(ix) = ix*dd

if( x(ix) .LE. DOMLEN )then

z(1,ix) = z1l

z(2,ix) = z2l

z(3,ix) = z3l

endif

1 continue

c

c compute cell size and cell centers

c

do 3 ix = 1,nx

dx(ix) = x(ix)-x(ix-1)

xc(ix)=0.5d0*(x(ix)+x(ix-1))

3 continue

c

c compute median dual cell

c

do 5 ix = 1,nx-1

mdc(ix) = xc(ix+1)-xc(ix)

5 continue

mdc(0) = xc(1)-x(0)

mdc(nx) = x(nx)-xc(nx)

c

v1 = 0.d0

do 7 ix = 1,nx

v1 = v1 + dx(ix)

7 continue

v2 = 0.d0

do 9 ix = 0,nx

v2 = v2 + mdc(ix)

9 continue

c write(6,*)'Check volumes',v1,v2,'should be equal'

open(22,file='test')

call check(z,mdc,ndof,nx)

144

return

end

subroutine check(z,mdc,ndof,nx)

c

c computed the amount of mass, energy and momentum being stored

c in the system

c on entry z must be the parameter vector

c

real*8 z(ndof,0:nx),mdc(0:nx)

real*8 checksum(3)

integer ix,j

call parm2cons(z,ndof,nx)

do j = 1,ndof

checksum(j) = 0.d0

enddo

do ix = 0,nx

do j = 1,ndof

checksum(j) = checksum(j) + z(j,ix)*mdc(ix)

enddo

enddo

call cons2parm(z,ndof,nx)

write(22,*)(checksum(j),j=1,ndof)

return

end

Subroutine spitw

subroutine splitw(zleft,zright,dx,dt,fleft,fright,

&phiw,al,dw,ncell)

implicit none

integer ndof,ncell

parameter (ndof=3)

include 'paramt.h'

double precision zleft(ndof),zright(ndof)

double precision fleft(ndof),fright(ndof)

double precision uleft(ndof),uright(ndof)

double precision fl(ndof),fr(ndof),df(ndof),heat(6)

double precision dt,dx

double precision al(*),eigen(ndof,ndof),zavg(ndof)

double precision dz(ndof),dw(*),phiw(*),help(ndof),source(ndof)

double precision ravg,Havg,uavg,asqr,aavg,ra,drho,dp,du,Mach

double precision toler,temp,helpme,dru,drE,druH,dH,dr,cost

double precision dpdrho,dpde,dpdm,PIM,PIE,PIR,eps,dotq

double precision pavg,chiavg,kappaavg,deltp,esse,DD,chic,kappac

double precision rleft,epsleft,pleft,chileft,kappaleft

double precision hmedio,hleft,hright,alpha,deltau,zero

double precision rright,epsright,pright,chiright,kapparight

parameter (TOLER=1.e-12)

integer k,j,i

logical lflag

double precision sumdf(3)

common /def/ sumdf

c dt = 1.e38

c

c gradient of the parameter vector and averaged state

145

c

do 1 i = 1, ndof

zavg(i) = 0.5d0 *(zright(i)+zleft(i))

dz(i) = (zright(i)-zleft(i))/dx

1 continue

c

c Roe averaged

c

Havg = zavg(2)/zavg(1)

uavg = zavg(3)/zavg(1)

ravg = zavg(1)*zavg(1)

alpha=zleft(1)/(zleft(1)+zright(1))

hleft=zleft(2)/zleft(1)-1/2*(zleft(3)/zleft(1))

&*(zleft(3)/zleft(1))

hright=zright(2)/zright(1)-1/2*(zright(3)/zright(1))

&*(zright(3)/zright(1))

deltau=zright(3)/zright(1)-zleft(3)/zleft(1)

hmedio=alpha*hleft+(1-alpha)*hright

&+1/2*alpha*(1-alpha)*deltau*deltau

c destra

cost=zright(3)*zright(3)/2-zright(1)*zright(2)

rright = zright(1)*zright(1)

call epsi(cost,rright,epsright)

call pchikappa(rright,epsright,pright,chiright,kapparight)

c sinistra

cost=zleft(3)*zleft(3)/2-zleft(1)*zleft(2)

rleft = zleft(1)*zleft(1)

call epsi(cost,rleft,epsleft)

call pchikappa(rleft,epsleft,pleft,chileft,kappaleft)

chic=(chiright+chileft)/2

kappac=(kapparight+kappaleft)/2

deltp=(pright-pleft)-chic*(rright-rleft)-kappac*(epsright-epsleft)

esse=chic+kappac*hmedio

DD=(esse*(rright-rleft))*(esse*(rright-rleft))

&+(pright-pleft)*(pright-pleft)

if ((rright-rleft).eq.0)then

chiavg=chic

else

chiavg=(DD*chic+esse*esse*(rright-rleft)*deltp)/(DD-(pright-pleft)

&*deltp)

endif

if ((epsright-epsleft).eq.0)then

kappaavg=kappac

else

kappaavg=DD*kappac/(DD-(pright-pleft)*deltp)

endif

c zero=chiavg*(zright(1)*zright(1)-zleft(1)*zleft(1))

c &+kappaavg*(epsright-epsleft)-(pright-pleft)

c if (abs(zero).gt.toler)then

c write(*,*)ncell

c write(*,*)zero

c write(*,*)chiavg,chic

c write(*,*)kappaavg,kappac

c pause

c endif

146

c derivate della pressione

c

dpdrho=pir(uavg,chiavg,kappaavg)

dpde=pie(uavg,chiavg,kappaavg)

dpdm=pim(uavg,chiavg,kappaavg)

asqr = dpdrho+(Havg-uavg*uavg)*dpde

if( asqr .LT. 0.d0 )then

write(6,*)'Negative a^2 ',asqr,' in cell no. = ',ncell

stop

endif

aavg = sqrt(asqr)

Mach = uavg/aavg

ra = ravg*aavg

c

c gradient of the primitive

c

drho = 2.*zavg(1)*dz(1)

dp = 1/(1+dpde)*((2*zavg(1)*dpdrho+zavg(3)*dpdm+zavg(2)*dpde)

&*dz(1)+zavg(1)*dpdm*dz(3)+zavg(1)*dpde*dz(2))

du = (-zavg(3)/zavg(1)*dz(1)+dz(3))/zavg(1)

c

c characteristic jumps (gradient)

c

dw(1) = (asqr*drho-dp)

dw(2) = 0.5*(dp+ra*du)

dw(3) = 0.5*(dp-ra*du)

c source terms for the characteristic variables

open(1,file='heat.txt')

read(1,*)(heat(i),i=1,6)

close(1)

dotq=heat(2)

c

source(1) = -dpde* dotq

source(2) = 0.5d0 * dpde * dotq

source(3) = source(2)

c

c autovalori

c

al(1) = uavg

al(2) = uavg+aavg

al(3) = uavg-aavg

c

c fluttuazione delle variabili caratteristiche

c

do i = 1,ndof

phiw(i) = (-al(i)*dw(i)+source(i))*dx

enddo

c

c right eigenvector (entropy)

c

eigen(1,1) = 1.d0/asqr

eigen(2,1) = Havg/asqr-1/dpde

eigen(3,1) = uavg/asqr

c

147

c right eigenvector (fast wave)

c

eigen(1,2) = 1.d0/asqr

eigen(2,2) = (Havg/asqr+uavg/aavg)

eigen(3,2) = (uavg/asqr+1./aavg)

c

c right eigenvector (slow wave)

c

eigen(1,3) = 1.d0/asqr

eigen(2,3) = (Havg/asqr-uavg/aavg)

eigen(3,3) = (uavg/asqr-1./aavg)

c

c compute minimum dt

c

do k = 1, ndof

dt = min( dt, dx/abs(al(k)) )

enddo

c

do j = 1, ndof

help(j) = 0.0

enddo

c

do 20 k = 1, ndof

do j = 1, ndof

help(j) = help(j) + eigen(j,k) * phiw(k)

enddo

if( al(k) .GT. 0.0 )then

do j = 1, ndof

fright(j) = fright(j) + eigen(j,k) * phiw(k)

enddo

elseif( al(k) .LT. 0.0 )then

do j = 1, ndof

fleft(j) = fleft(j) + eigen(j,k) * phiw(k)

enddo

else

goto 20

endif

20 continue

do j = 1, ndof

sumdf(j) = sumdf(j) + help(j)

enddo

c

c test the decomposition:

c check that f_R-f_L = \sum ...

c the sum is stored in help

c

call flux(fl,zleft,ndof)

call flux(fr,zright,ndof)

source(1) = 0.d0

source(2) = dotq

source(3) = 0.d0

lflag = .FALSE.

do j = 1, ndof

c temp normalizes

temp = abs(0.5*(fl(j)+fr(j)))

if(temp.LE.TOLER)temp = 1.d0

df(j) = -(fr(j)-fl(j))

if(abs(help(j)-df(j))/abs(temp) .GT. TOLER)lflag= .true.

enddo

c if(lflag)then

148

c write(6,*)"cell no. ",ncell

c do j = 1, ndof

c write(6,*)j,df(j),help(j)

c enddo

c write(6,*)"dw",(dw(j),j=1,ndof)

c write(6,*)"al",(al(j),j=1,ndof)

c write(6,*)"phiw",(phiw(j),j=1,ndof),du

c write(6,*)"u,a,",uavg,aavg

c pause

c endif

12 continue

1333 format(1(E12.6,1X))

c

return

c

c temp = (2.+gm1*Mach*(Mach-2.))/(2.*gm1)

do j = 1,ndof

uleft(j) = zleft(j)

uright(j) = zright(j)

enddo

call parm2cons(uleft,ndof,0)

call parm2cons(uright,ndof,0)

c

c jumps in ru, rE, ruH, ru^2+p H

c

dr = dp/asqr

dru = ravg*(1-Mach)*du

drE = temp*dp

druH = aavg*(Mach-1)*temp*dp

c

helpme = (0.5*gm1*uavg*uavg-Havg)*dr+gam*drE-gm1*uavg*dru

helpme = helpme/ravg

c

return

c if(ncell.EQ.51)then

c write(36,*)fr(1)-fl(1),dru

c write(37,*)uright(2)-uleft(2),drE

c write(38,*)fr(2)-fl(2),druH

c write(39,*)fr(3)-fl(3),dp*(Mach-1)**2

c dH = zright(2)/zright(1)-zleft(2)/zleft(1)

c write(40,*)dH,(1.d0-Mach)*dp/ravg,helpme

c write(41,*)Mach,dp

c endif

c

return

end

!derivata della pressione rispetto all'energia

double precision function PIE(uavg,chi,kappa)

double precision uavg,chi,kappa

PIE = kappa

return

end

!derivata della pressione rispetto alla quantità di moto

double precision function PIM(uavg,chi,kappa)

double precision uavg,chi,kappa

PIM = -uavg*kappa

149

return

end

!derivata della pressione rispetto alla densità

double precision function PIR(uavg,chi,kappa)

double precision uavg,chi,kappa

PIR = chi + 0.5*uavg*uavg*kappa

return

end

Subroutine epsi

subroutine epsi(cost,rhos,zero)

!dichiarazione variabili

include 'aria.inc'

double precision toll, zero, cost,rhos

double precision press(Ne)

external funz

toll=1d-7

!cost=-3.05d5*1.17

!rhos=1.17

call Pepsilon(press,rhos)

call fzero(toll,zero,cost,press)

!write(*,*) zero

return

end

Funzione 21

2

3 2, zzzpF

REAL FUNCTION funz(x,cost,press,eii)

include 'aria.inc'

integer je

double precision x, cost,pressx

double precision press(*),eii(*)

!pressx è da calcolare con interpolazione capendo in che intervallo si trova

je=0

do j=1,Ne

if( (eii(j) .LT. x).or. (eii(j) .eq. x))then

je=je+1

endif

enddo

150

pressx=press(je)+(press(je+1)-press(je))/(eii(je+1)-eii(je))*(x-eii(je))

funz=x+cost+pressx

END FUNCTION funz

Subroutine fzero

subroutine fzero(toll,zero,cost,press)

include 'aria.inc'

double precision cost,err, toll, zero, a, b, m

double precision press(*)

external funz

a=ei(1)

b=ei(Ne-1)

if ((funz(a,cost,press,ei)*funz(b,cost,press,ei)).gt.0d0) then

write(*,*) 'intervallo non valido'

stop

endif

err=1d0

do while (err.gt.toll)

m=.5d0*(a+b)

err=b-a

if (funz(b,cost,press,ei)*funz(m,cost,press,ei).lt.0) then

a=m

else

b=m

endif

enddo

zero=a

return

end

Subroutine interp1D

subroutine interp1D(jr,je,rho,p,rhos,peps)

integer Nr,Ne,jr,je

parameter (Nr=80)

parameter (Ne=100)

double precision p(Nr,Ne),rho(Nr),rhos,peps

151

peps=p(jr,je)+(p(jr+1,je)-p(jr,je))/(rho(jr+1)-rho(jr))*(rhos-rho(jr))

return

end

Subroutine interp2D

subroutine interp2D(jr,je,rho,ei,p,rhos,es,press)

integer Nr,Ne,jr,je

parameter (Nr=80)

parameter (Ne=100)

double precision p(Nr,Ne),ei(Ne),rho(Nr)

real*8 rhos,es,press1,press2,press

press1=p(jr,je)+(p(jr+1,je)-p(jr,je))/(rho(jr+1)-rho(jr))*(rhos-rho(jr))

press2=p(jr,je+1)+(p(jr+1,je+1)-p(jr,je+1))/(rho(jr+1)-rho(jr))*(rhos-rho(jr))

press=press1+(press2-press1)/(ei(je+1)-ei(je))*(es-ei(je))

return

end

Subroutine pchikappa

subroutine pchikappa(rhos,es,press,chis,kappas)

!dichiarazione variabili

include 'aria.inc'

integer jr,je

real*8 rhos,es,press,chis,kappas

jr=0

do j=1,Nr

if( rho(j) .LT. rhos )then

jr=jr+1

endif

enddo

!var1=(rho .LT. rhos)

!jr=sum(var1)

je=0

do j=1,Ne

if( ei(j) .LT. es )then

je=je+1

endif

enddo

call interp2D(jr,je,rho,ei,p,rhos,es,press)

call interp2D(jr,je,rho,ei,chi,rhos,es,chis)

152

call interp2D(jr,je,rho,ei,kappa,rhos,es,kappas)

return

end

Subroutine pepsilon

subroutine Pepsilon(press,rhos)

!dichiarazione variabili

integer jr,je

include 'aria.inc'

double precision press(*)

double precision rhos

jr=0

do j=1,Nr

if( rho(j) .LT. rhos )then

jr=jr+1

endif

enddo

do je=1,Ne

call interp1D(jr,je,rho,p,rhos,press(je))

enddo

close(13)

return

end

Subroutine bcs

subroutine bcs(z,du,ndof,nx,dt,p1,r1,u1,p2) implicit none

integer nx,ndof

double precision sigma1, sigma2

parameter(sigma1=0.5,sigma2=0.5)

double precision dt,dx

double precision z(ndof,0:nx),du(ndof,0:nx)

double precision ux,rho,h,cs,csqr,eps,press,chi,k,M

double precision uav,rav,hav,csav,csavq,epsav,pav,chiav,kav,Mav

double precision fl(3),fr(3),df(3),sumdf(3)

double precision p1,r1,u1,p2,cc1,cc2,cos

double precision zavgN(ndof),zavg0(ndof)

common /abc/df

common /def/sumdf

integer i,j

include 'paramt.h'

include 'nozzle.com'

c

153

dx=length/nx

c

c node 0 (densità e velocità costanti)

c Roe averaged nella cella 1

do i=1,ndof

zavg0(i)=0.5d0*(z(i,0)+z(i,1))

enddo

c

rav = zavg0(1)*zavg0(1)

hav = zavg0(2)/zavg0(1)

uav = zavg0(3)/zavg0(1)

c

c Pressione, chi e k nella cella 1

cos=rav*uav*uav/2-zavg0(1)*zavg0(2)

call epsi(cos,rav,epsav)

call pchikappa(rav,epsav,pav,chiav,kav)

c

c Velocità del suono mediata nella cella 1

csavq=chiav+(hav-0.5*uav**2)*kav

if(csavq.LT.0.d0)then

write(6,*)'Problem with cs^2 ',csavq

stop

endif

csav=sqrt(csavq)

c

c Variabili primitive nel nodo 0

rho = z(1,0)*z(1,0)

h = z(2,0)/z(1,0)

ux = z(3,0)/z(1,0)

cos=rho*ux*ux/2-z(1,0)*z(2,0)

call epsi(cos,rho,eps)

call pchikappa(rho,eps,press,chi,k)

c Velocità del suono nel nodo 0

csqr=chi+(h-0.5*ux**2)*k

if(csqr.LT.0.d0)then

write(6,*)'Problem with cs^2 ',csqr

stop

endif

cs=sqrt(csqr)

c ********************************************************

c Charcteristic boundary conditions Poinsot - Lee

c

cc1=sigma1*cs*dx/length

c contributo dovuto a Lo

du(1,0)=du(1,0)-(1/csavq)*(rho-r1)*csqr*cc1

du(2,0)=du(2,0)-(hav/csavq-1/kav)*(rho-r1)*csqr*cc1

du(3,0)=du(3,0)-(uav/csavq)*(rho-r1)*csqr*cc1

c contributo dovuto a L+

du(1,0)=du(1,0)-(1/csavq)*rho*cs*(ux-u1)*cc1

du(2,0)=du(2,0)-(hav/csavq+uav/csav)*rho*cs*(ux-u1)*cc1

du(3,0)=du(3,0)-(uav/csavq+1/csav)*rho*cs*(ux-u1)*cc1

154

c ********************************************************

c Velocità e densità costanti nel nodo 0

c du(1,0)=du(1,0)+(r1-rho)*uav+(u1-ux)*rav

c du(2,0)=du(2,0)+(r1-rho)*uav*(0.5*uav**2-chiav/kav)+

c +rav*(hav+uav**2)*(u1-ux)

c du(3,0)=du(3,0)+(u1-ux)*2*rav*uav+(r1-rho)*uav**2

open(12,file='pp.dat')

write(12,FMT=400) rho-r1,ux-u1

c

c -------------------------------------------------------------------

c call flux(fr,z(1,NX),ndof)

c call flux(fl,z(1,0),ndof)

c do j =1,ndof

c df(j) = df(j) + (fr(j)-fl(j))*dt

c enddo

c write(42,FMT=400)((fr(j)-fl(j)),j=1,ndof),(sumdf(j),j=1,ndof)

c -------------------------------------------------------------------

c

c node NX (pressione costante)

c Roe averaged nella cella nx

c

do i=1,ndof

zavgN(i)=0.5d0*(z(i,nx)+z(i,nx-1))

enddo

rav = zavgN(1)*zavg0(1)

hav = zavgN(2)/zavg0(1)

uav = zavgN(3)/zavg0(1)

c Pressione, chi e k nella cella nx

cos=rav*uav*uav/2-zavgN(1)*zavgN(2)

call epsi(cos,rav,epsav)

call pchikappa(rav,epsav,pav,chiav,kav)

c

c Velocità del suono media nella cella nx

csavq=chiav+(hav-0.5*uav**2)*kav

if(csavq.LT.0.d0)then

write(6,*)'Problem with cs^2 ',csavq

stop

endif

csav=sqrt(csavq)

c Variabili primitive nel nodo Nx

rho = z(1,Nx)*z(1,Nx)

h = z(2,Nx)/z(1,Nx)

ux = z(3,Nx)/z(1,Nx)

cos=rho*ux*ux/2-z(1,nx)*z(2,nx)

call epsi(cos,rho,eps)

call pchikappa(rho,eps,press,chi,k)

press=z(1,nx)*z(2,nx)-eps-z(3,nx)*z(3,nx)/2

c velocità del suono nel nodo Nx

csqr=chi+(h-0.5*ux**2)*k

if(csqr.LT.0.d0)then

write(6,*)'Problem with cs^2 ',csqr

stop

endif

cs=sqrt(csqr)

155

c ****************************************************************

c Charcteristic boundary conditions Poinsot - Lee

M=ux/cs

Mav=uav/csav

cc2=sigma2*dx/length

c cc2=sigma2*(1-M**2)*dx/length

c contributo dovuto a L-

du(1,nx)=du(1,nx)-(1/csav)*(press-p2)*cc2

du(2,nx)=du(2,nx)-(hav/csav-uav)*(press-p2)*cc2

du(3,nx)=du(3,nx)-(uav/csav-1)*(press-p2)*cc2

c ****************************************************************

c du(2,nx)=du(2,nx)+(kav/(kav+1))*uav*(press-p2)

c du(3,nx)=du(3,nx)+(press-p2)

open(15,file='pp1.dat')

write(15,FMT=400) press-p2

400 FORMAT(6(E14.6,1X))

return

end

Subroutine parm2cons

SUBROUTINE parm2cons(ZROE,NDOF,NX) C

C $Id:$

C

IMPLICIT NONE

C

INCLUDE'paramt.h'

C

C .. Array Arguments ..

C

INTEGER NDOF,NX

LOGICAL LFLAG

DOUBLE PRECISION ZROE(NDOF,0:NX)

DOUBLE PRECISION cos,ep

C

C .. Local Scalars ..

C

INTEGER IX

DOUBLE PRECISION SUM

C

C .. External Functions ..

C

C

C .. Executable Statements ..

C

do 100 IX = 0 , NX

SUM = ZROE(3,IX)**2

cos=SUM/2-ZROE(1,IX) * ZROE(2,IX)

156

call epsi(cos,ZROE(1,IX)**2,ep)

ZROE(2,IX) = ep + SUM/2

C

ZROE(3,IX) = ZROE(1,IX) * ZROE(3,IX)

ZROE(1,IX) = ZROE(1,IX) * ZROE(1,IX)

100 CONTINUE

c

RETURN

END

SUBROUTINE cons2parm(ZROE,NDOF,NX)

C

IMPLICIT NONE

C

c

c This routine transforms the conservative vector into the

c paramter vector overwriting the array Z.

c

INCLUDE'paramt.h'

C

C .. Array Arguments ..

C

INTEGER NDOF,NX

LOGICAL LFLAG

DOUBLE PRECISION ZROE(NDOF,0:NX)

C

C .. Local Scalars ..

C

INTEGER IX,NERR,IOPT

DOUBLE PRECISION SUM, ep,press,chis,kappas

C

C .. External Functions ..

C

C

C .. Executable Statements ..

C

DO 100 IX = 0 , NX

C

IF(ZROE(1,IX) .LE. 0.d0 )THEN

WRITE(6,FMT=105)ZROE(1,IX),IX

105 FORMAT('CONS_TO_PARM',1X,'NEG. DENSITY ',F7.3,' IN NODE ',I6)

STOP

ENDIF

C

ZROE(1,IX) = SQRT(ZROE(1,IX))

ZROE(3,IX) = ZROE(3,IX) / ZROE(1,IX)

SUM = ZROE(3,IX)**2

ep=ZROE(2,IX)-SUM/2

call pchikappa(ZROE(1,IX)**2,ep,press,chis,kappas)

ZROE(2,IX) = (ep+press+SUM/2)/ZROE(1,IX)

100 CONTINUE

C

RETURN

END

157

Subroutine update

subroutine update(u,du,ndof,mdc,t,dt,nx,iok)

implicit none

integer nx,ndof,iok

double precision u(ndof,0:nx),du(ndof,0:nx),mdc(0:nx)

double precision dt,t

double precision TOLER, EPS

parameter(TOLER=1.d-12,EPS=5.d-11)

double precision errmax

!

! u is conserved variables ......

!

double precision err(3),duji,denom

integer i,j

!

! update interior nodes only

!

open(7,FILE='conv.dat')

iok=1

errmax=0.d0

do j =1,ndof

err(j) = 0.d0

enddo

! dt = 1.d38

! do 4 i = 1, nx-1

! dt = min(dt, mdc(i)/abs(u(i)) )

! 4 continue

do 1 i = 0, nx

do 1 j = 1,ndof

duji = du(j,i)

denom = u(j,i)**2

if(denom.LE.TOLER)denom = 1.d0

err(j) = err(j) + duji*duji/denom

u(j,i) = u(j,i) + dt/mdc(i)*duji

1 continue

do j = 1, ndof

err(j) = sqrt(err(j))/(nx+1)

errmax = dmax1(err(j),errmax)

enddo

write(6,FMT=120)t,dt,(err(j),j=1,ndof)

write(7,FMT=1333)t,(err(j),j=1,ndof)

if(abs(errmax).LT.EPS)iok = 0

120 FORMAT('t= ',F12.8,' dt= ',F12.8,' err= ',3(1X,E12.6))

1333 FORMAT(6(E12.6,1X))

return

end