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].
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(,;
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
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
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.
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.
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)
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)
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.
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
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
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
Π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
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