1

Capitolo 1

Flussi euleriani non stazionari bidimensionali per miscele di gas ideali comprimibili.

Trattazione analitica.

1.1 Flussi monodimensionali non stazionari per fluidi costituiti da un solo gas perfetto.

L’evoluzione di tali flussi è governata dal seguente sistema di equazioni di Eulero:

( ) ( 1.1.1 )

( )

Nella ( 1.1.1 ) il sistema di Eulero è stato espresso in forma di divergenza.

Con ed ( ) si sono indicati rispettivamente il vettore delle variabili conservate ed il vettore dei

flussi generalizzati, definiti in funzione delle seguenti grandezze, descriventi lo stato termodinamico

della particella di gas:

la densità

la velocità

la pressione

l’energia totale per unità di massa

l’entalpia totale per unità di massa

andando ad esplicitare le grandezze che caratterizzano rispettivamente l’energia totale e l’entalpia

totale per unità di massa si ottiene:

energia interna per unità di massa

entalpia interna per unità di massa

2

Per riportarsi alle quantità che compaiono nel sistema di Eulero espresso nella formulazione ( 1.1.1 ) ,

occorre definire le seguenti grandezze: l’energia totale per unità di volume e l’entalpia totale per

unità di volume:

(energia totale per unità di volume)

(entalpia totale per unità di volume)

Sotto l’ipotesi di gas termicamente perfetto si può introdurre la seguente equazione di stato (valida

appunto nel caso di gas perfetti):

( 1.1.2 )

Dove è la costante del gas definita come :

( 1.1.3 )

nella quale con [ ] si indica la costante universale dei gas e con

[ ] il peso molecolare del gas.

Introducendo l’ulteriore ipotesi di gas caloricamente perfetto le grandezze energia interna per unità

di massa ed entalpia interna per unità di massa possono essere espresse nella seguente forma:

( 1.1.4 )

( 1.1.5 )

dove con e si indicano rispettivamente il calore specifico a pressione costante e a volume

costante.

Si definisce con la quantità adimensionale ottenuta dal rapporto dei i calori specifici :

( 1.1.6 )

I calori specifici possono essere espressi in funzione di e di nel seguente modo:

( 1.1.7 )

( 1.1.8 )

Per completezza si riporta anche la relazione che lega tra loro i calori specifici e la costante della

miscela:

( 1.1.9 )

3

Sfruttando la definizione dell’ energia totale per unità di volume, l’equazione di stato ( 1.1.2 ) e

l’espressione ( 1.1.8 ) relativa al calore specifico a volume costante è possibile ricavare una

formulazione, che verrà chiamata “equazione di stato per esteso” , per mezzo della quale si potrà

ottenere la chiusura del sistema di Eulero ( 1.1.1 ) .

(

)

(

)

nella quale una volta sostituita l’equazione di stato:

è possibile ottenere :

(

)

(1.10)

Tale equazione, utilizzando la definizione dell’energia interna per unità di volume, è riscrivibile nella

seguente forma equivalente:

(1.11)

Concludendo, si può affermare che il sistema delle equazioni di Eulero rappresentante la dinamica di

un fluido costituito da un solo gas termicamente e caloricamente perfetto , in campo

monodimensionale non stazionario è integrabile, cioè se ne può ricavare la soluzione analitica. Le

variabili dipendenti sono , e , mentre le variabili indipendenti sono . Come espresso dalla

( 1.1.10 ) , o analogamente dalla ( 1.1.11 ), la chiusura del sistema è garantita dall’aggiunta, come dato

in ingresso, della variabile costitutiva ed adimensionale ( rapporto dei calori specifici ).

E’ da tener presente che per la determinazione di alcune grandezze derivate, quali ad esempio la

temperatura è necessario aggiungere alla un’ altra variabile costitutiva, ad esempio (o

analogamente in base alla ( 1.1.3 ) il peso molecolare ).

In alternativa, per chiudere il sistema e ricavare la temperatura, in luogo della coppia di variabili

costitutive ed è possibile assegnare la coppia relativa ai calori specifici del gas e , o una

qualsiasi altra coppia, di grandezze costitutive appena citate purché indipendenti tra loro, ed utilizzare

eventualmente le ( 1.1.7) , ( 1.1.8 ) e ( 1.1.9 ) per ricavare le grandezze necessarie.

4

1.2 Flussi monodimensionali non stazionari per un fluido costituito da una miscela di gas perfetti.

Introduzione del modello Termodinamico Standard

Nella seguente trattazione viene fatto un passo avanti rispetto al paragrafo precedente, rendendo il

discorso più generale, e cioè prendendo in analisi fluidi costituiti da miscele di gas.

Si assegnano le seguenti ipotesi:

o (1) mescolamento isotermo tra i gas costituenti la miscela

o (2) che ogni specie di gas occupi lo stesso volume

o (3) che ogni particella abbia la stessa velocità.

o (4) che ogni specie di gas costituente la miscela sia termicamente e caloricamente perfetta

Si procede ora con la definizione delle grandezze che caratterizzano la miscela di gas. In questa

trattazione ogni grandezza, definita nel precedente paragrafo per un fluido ad una sola specie, verrà

tradotta nel caso multispecie.

Densità della miscela:

se si indica con la massa della k_esima specie e con il volume occupato dalla miscela di gas, per

l’ipotesi (2), è possibile definire la densità riferita alla k_esima specie come:

( 1.2.12 )

Pertanto così definita la rappresenterebbe la densità del gas qualora nel volume sia presente la

sola k_esima specie.

Sempre in base all’ipotesi (2), cioè che i gas occupino tutti lo stesso volume si può definire la densità

della miscela come:

∑ ( 1.2.13 )

conseguentemente possono definirsi le “frazioni in massa” dette anche “concentrazioni”:

; ( 1.2.14 )

Velocità della miscela :

5

definita direttamente dall’ ipotesi (3) .

Pressione della miscela:

si definisce con la pressione riferita alla k_esima specie di gas qualora fosse l’unica ad

occupare l’intero volume del dominio di definizione spaziale.

Per la legge di Dalton:

∑ ( 1.2.15 )

Ovviamente in base all’ipotesi (1) e (4) per ogni k_esimo gas risulta valida la legge di stato:

( 1.2.16 )

dove:

( 1.2.17 )

avendo indicato con il peso molecolare della k_esima specie di gas.

energia interna per unità di volume relativa alla miscela:

dalla definizione di energia interna per unità di volume relativa alla k_esima specie:

( 1.2.18 )

andando ad eseguire la sommatoria per ogni specie presente, si ottiene:

∑ ( 1.2.19 )

Calore specifico a volume costante per la miscela:

considerando l’ipotesi (4) cioè di gas termicamente e caloricamente perfetto, l’energia interna per

unità di massa relativa alla k_esima specie è data dalla:

( 1.2.20 )

6

e sostituendo le ( 1.1.4 ) , ( 1.2.20 ) nella ( 1.2.19 ) , si ottiene :

∑

ricordando la definizione ( 1.2.14 ) sulle frazioni in massa, e considerando l’ipotesi ( 1 ) di

mescolamento isotermo delle specie di gas si ottiene infine:

∑ ( 1.2.21 )

Calore specifico a pressione costante per la miscela

Analogamente a quanto appena svolto per ricavare il calore specifico a volume costante, e quindi

sfruttando la definizione dell’energia interna per unità di volume, si procede col ricavare il calore

specifico a pressione costante utilizzando questa volta la definizione dell’entalpia interna per unità

di volume:

∑

dove è proprio l’entalpia interna per unità di volume associata alla k_esima specie.

Nel caso di gas caloricamente perfetto l’entalpia interna è esprimibile come:

infine:

∑ ( 1.2.22 )

Rapporto dei calori specifici per la miscela:

∑

∑

7

Una volta definite le grandezze che caratterizzano le variabili di stato in gioco, si passa a verificare

come anche nel caso di una miscela di gas il sistema di Eulero sia chiuso e di quale sia stavolta la

formulazione dell’ “equazione di stato per esteso”, (cioè funzione dell’energia interna per unità di

volume e di ) necessaria alla chiusura del sistema.

Partendo dalla definizione dell’ energia totale per unità di volume relativa alla k_esima specie di gas,

costituente la miscela, e sfruttando l’ipotesi (3) (sulla velocità delle particelle) si ottiene :

( 1.2.24 )

Pertanto l’energia totale relativa alla miscela è ottenibile attraverso la sommatoria:

∑

∑

∑

ed introducendo la ( 1.2.20 ) , cioè , e la ( 1.2.14 ) ossia la definizione della

concentrazione per la k_esima specie, è possibile riscrivere l’energia totale per unità di volume

come:

∑

∑

considerando a questo punto le ( 1.2.13 ) , ( 1.2.20 ) si ottiene:

Per l’equazione di stato espressa dalla ( 1.1.2 ) e per la definizione del calore specifico a volume

costante ( 1.1.8 ) in questa analisi riferito alla miscela, si ricava:

e quindi risolvendo rispetto alla pressione:

o analogamente:

( 1.2.25 )

8

cioè si è ritrovata la ( 1.1.9 ) e quindi “ l’equazione di stato per esteso” ricavata nel caso di fluido

costituito da una sola specie di gas. In questo caso però le grandezze in questione sono di miscela,

e bisogna porre attenzione nella definizione del .

Infatti se nel caso di fluido monocomponente il viene definito dalla ( 1.1.6 ) semplicemente

come rapporto tra i calori specifici del gas, in questo caso , come indicato nella ( 1.2.23 ) è

necessario conoscere anche le concentrazioni dei gas costituenti la miscela, nuova incognita del

problema.

Discorso analogo per la determinazione della costante di miscela:

∑ ( 1.2.26 )

Introduzione del modello Termodinamico Standard

Al fine di ricavare la soluzione nel caso in cui il fluido sia composto da una miscela di gas bisogna

modificare il sistema di Eulero ( 1.1.1 ) in modo tale da renderlo in grado di poter ricavare

l’evoluzione delle concentrazioni.

Un tale modello di equazioni di Eulero è noto in letteratura come “modello Termodinamico

Standard”, ed è caratterizzato dal fatto che l’equazione di conservazione della massa viene

espressa in termini delle densità parziali, e pertanto considerata tante volte quante sono le specie

di gas costituenti la miscela.

Grandezze derivate (per il modello Termodinamico Standard).

Definizione della velocità del suono nel caso di un fluido costituito da k specie di gas.

Per definizione la velocità del suono corrisponde alla velocità di propagazione dei disturbi nel caso

di flusso omentropico:

(

)

( 1.2.27 )

dalle ipotesi di mescolamento isotermico e che ogni specie di gas occupi lo stesso volume si sono

potute scrivere la legge di Dalton ( 1.2.15 ) e l’ equazione di stato per la k_esima specie ( 1.2.16 )

che per chiarezza vengono di seguito riportate:

∑

9

pertanto si può scrivere:

∑

∑ ((

)

(

)

)

per la quale le derivate parziali contenute, grazie all’equazione di stato, assumono la forma :

(

)

( 1.2.28a)

(

)

( 1.2.28b)

Le (1.2.28) sono state ottenute rispettivamente per differenziazione dell’equazione di stato

tenendo nella prima costante la temperatura e nella seconda la densità della

k_esima specie.

Pertanto:

∑

( 1.2.29 )

Ora si procede col ricavare un’ espressione di in funzione dell’energia interna per unità di

volume della miscela :

si parte differenziando la (1.2.19) e cioè l’equazione ∑

∑

∑

∑

considerando che per un gas caloricamente perfetto è possibile scrivere :

∑

∑

10

Ed andando a risolvere in termini di :

∑

∑

se inoltre si tiene presente che, ∑

:

∑

( 1.2.30 )

Inserendo quest’ultima nella ( 1.2.21 ), e considerando la ( 1.2.20 ) :

∑

( ∑

) ∑

∑(

)

A questo punto si ha la da inserire nell’espressione della velocità del suono

(

)

∑ (

)

(

)

e considerando la ( 1.2.28a)

(

) ∑ (

)

che può essere riscritta nella forma:

∑

( ∑

)

A questo punto ricordando che:

11

∑

∑

e considerando la definizione di entalpia interna per unità di massa si ottiene :

∑

Si è pertanto dimostrato come per un fluido multispecie il sistema di Eulero sia sostanzialmente lo

stesso di quello utilizzato nel caso di fluido monocomponente, a differenza del fatto che ora le

grandezze in questione sono riferite alla miscela e quindi ottenibili in funzione delle concentrazioni,

ricavabili a partire dal sottosistema delle equazioni di conservazione della massa scritte in termini di

densità parziali e sfruttando le (1.2.13) , (1.2.14).

Si riporta di seguito il sistema di equazioni che caratterizza la struttura del modello Termodinamico

Standard, d’ora in poi indicato anche come T.S. :

{

( 1.2.31 )

Tale sistema è chiuso il momento in cui si considera “l’equazione di stato per esteso”:

Pertanto al sistema ( 1.2.31 ) è necessario aggiungere le seguenti relazioni:

∑

; k=1,…,K

12

∑

∑

Così facendo si è garantita la chiusura del sistema, infatti note ora le densità dei componenti è

possibile ottenere la densità della miscela, valutare le frazioni in massa, ricavare il rapporto dei calori

specifici e quindi ottenere “l’equazione di stato per esteso”.

1.3 Flussi bidimensionali non stazionari multispecie di gas ideali:

Il sistema di Eulero ( 1.1.1 ) in campo bidimensionale nel caso si consideri il modello Termodinamico Standard

assume la seguente struttura, (formulazione quasi lineare):

{

( )

( )

nel caso in cui il fluido sia costituito da una miscela di due gas il precedente sistema diventa :

{

( )

( )

( )

13

tale sistema è inoltre rappresentabile nella seguente formulazione matriciale:

( 1.3.32 )

nella quale con si è indicato il vettore delle grandezze fondamentali :

mentre con e le matrici Jacobiane:

[ ]

[ ]

Approssimazione per mezzo di una linearizzazione locale e risoluzione del problema 2_D per

sovrapposizione degli effetti

Se si suppongono note le proprietà di una particella, posizionata in un punto generico del campo

bidimensionale, il sistema (1.3.32 ) è in grado di risolvere il problema dell’ evoluzione di tali proprietà,

nell’intorno infinitesimo del punto.

Per integrare nel tempo il sistema ( 1.3.32 ) , cioè per ricavarne la soluzione ad una distanza finita dal

punto di partenza, nel caso bidimensionale non risulta conveniente percorrere la stessa strada

utilizzata per il caso monodimensionale, applicando quindi la teoria delle caratteristiche, (argomento

del prossimo paragrafo ) direttamente al sistema (1.2.31).

Per ottenere la soluzione del sistema (1.3.32) è comodo invece inquadrare il problema direttamente

dal punto di vista numerico, essendo oltretutto questa la strada che verrà seguita nel corso dei

prossimi capitoli.

14

Ragionare dal punto di vista numerico vuol dire avere a che fare con elementi di dimensione finita,

pertanto (nell’intorno del punto) il dominio spazio temporale dovrà essere discretizzato in piccoli

elementi di ampiezza finita rispettivamente di dimensioni Δt, Δx, Δy .

Per mezzo di un’approssimazione lineare la soluzione temporale della può essere valutata in termini

del seguente sviluppo in serie di Taylor al primo ordine:

Δ

Se ora nell’intorno del punto si ipotizza che i coefficienti moltiplicativi delle derivate ossia le

componenti delle matrici Jacobiane ed siano costanti all’interno del volume Δ Δ Δ

(e pertanto nell’intervallo d’integrazione temporale tra n ed n+1), si passa da un sistema che è quasi

lineare ad un sistema di tipo lineare, e quindi la soluzione potrà essere ottenuta per sovrapposizione

degli effetti. L’approssimazione introdotta consente di riscrivere il generico problema bidimensionale

come sovrapposizione di due fronti d’onda onda piani. La sovrapposizione degli effetti è possibile in

quanto il problema è stato preventivamente linearizzato.

Una tale approssimazione risulterà ovviamente possibile anche in campo continuo ma avrà comunque

validità nell’intorno infinitesimo del punto. Sarà quindi possibile ottenere una soluzione approssimata

della ( 2.3.32 ) come:

( 1.3.33 a )

dove :

( 1.3.33 b )

( 1.3.33 c )

Le ( 1.3.33 b ) e ( 1.3.33 c ) rappresentano nel continuo due casi particolari di onde bidimensionali,

rispettivamente un’onda piana che propaga in direzione x ed una che propaga in direzione y.

Se ci si limita a considerare solo il problema dell’evoluzione di un’onda piana e quindi considerando

separatamente la ( 1.3. 33 b ) o la ( 1.3. 33 c ) , è possibile ricercare la soluzione per via analitica

utilizzando la teoria delle caratteristiche, seguendo un percorso analogo a quello che si intraprende

nel caso 1_D non stazionario.

In effetti il fronte d’onda piano costituisce un caso particolare di una generica onda bidimensionale e

la sua evoluzione è risolvibile analiticamente, sfruttando la teoria delle caratteristiche, similmente a

quanto fatto per il caso 1_D .

15

Infatti per ottenere il sistema di equazioni di Eulero necessario a descrivere l’evoluzione di un’onda

piana si parte dalla sua generale formulazione 2D la ( 1.3.32 ) e si impongono le seguenti condizioni (

in particolare per un’onda piana che propaga in direzione x) :

o siano assenti o nulle le variazioni spaziali, gradienti e/o salti in direzione y

o sia in generale diversa da zero la componente della velocità in direzione y.

Sistema di Eulero per l’evoluzione di onde piane:

Il caso particolare in analisi restituisce pertanto un’approssimazione delle onde bidimensionali,

considerandole come fronti piani ed ovviamente non percepisce eventuali curvature delle onde

stesse.

Tale modello può considerarsi oltre che come un caso particolare del 2D non stazionario, anche come

un’ estensione del modello 1-D non stazionario, visto che come verrà dimostrato, entrambi i modelli

saranno caratterizzati dalle stesse proprietà iperboliche.

Il sistema di Eulero (modello termodinamico standard) che descrive la propagazione di un’onda piana

in direzione x per un fluido costituito da una miscela di gas, è dato dalla ( 1.3.33 b ) che in forma di

divergenza diventa:

{

( 1.3.34 a)

ed analogamente per l’onda piana in direzione y :

{

16

( 1.3.34 b)

Rispetto al caso 1-D compare in aggiunta l’equazione della conservazione quantità di moto in

direzione y ( o in direzione x nel caso si consideri la ( 1.3.34 b ) ).

Nella seguente trattazione si lavorerà unicamente in termini di fronti piani che propagano in

direzione x e quindi di (1.3.32a), essendo lo studio legato alla (1.3.34 b) e quindi al fronte piano in

direzione y, perfettamente analogo.

Per l’equazione di conservazione dell’energia bisogna tener presente che nell’esprimere l’energia

totale per unità di volume, per quanto riguarda il termine relativo all’energia cinetica si dovrà

considerare anche il contributo della componente della velocità in direzione y :

energia totale per unità di massa

entalpia totale per unità di massa

Anche in questo caso il sistema ( 2.3.32 a) risulta chiuso una volta introdotte le:

∑

∑

∑

Ovviamente in questo particolare caso 2_D “l’equazione di stato per esteso” ha la seguente

espressione :

(

)

17

( 1.3.35 )

Calcolo della temperatura:

se si assegnano le costanti di miscela (oppure i calori specifici a pressione e a volume costante

utilizzando la ( 1.1.9 ) ) è possibile ricavare la costante della miscela come:

∑

la temperatura sarà ricavata conseguentemente attraverso l’equazione di stato:

1.4 Proprietà iperboliche del sistema di Eulero

In questo paragrafo verranno messe a confronto le proprietà iperboliche relative ai sistemi di Eulero

descriventi rispettivamente l’evoluzione di un flusso 1-D e di un’onda piana.

Come primo punto si procede con lo studio delle proprietà iperboliche relativamente al caso 1_D

Il sistema ( 1.2.31) nel caso 1-D a due specie espresso in formulazione di divergenza assume la forma:

{

Ed in formulazione quasi lineare, (comodo punto di partenza per affrontare lo studio con la teoria

delle caratteristiche) :

e quindi:

18

{

( 1.4.36 )

(

)

( 1.4.37)

[

]

( 1.4.38 )

è il vettore avente per componenti le grandezze di stato fondamentali

è la matrice Jacobiana, costituita cioè dai coefficienti moltiplicativi delle derivate delle variabili di

stato in direzione x.

Calcolo degli autovalori

Ottenimento dell’ equazione caratteristica:

( )

[

]

Si applica il criterio di Laplace per ricavarne il determinante:

Tale equazione caratteristica ha per radici:

19

Tali radici rappresentano gli autovalori del sistema, grandezze che descrivono nello spazio fisico (in

questo caso piano fisico, spazio_tempo) le direzioni lungo le quali, (attraverso lo studio degli

autovettori) sarà possibile ricavare le variabili d’onda (grandezze scalari rappresentanti la soluzione

del problema nello spazio delle soluzioni ).

Pertanto ora lo studio prosegue andando a ricavare quali proprietà viaggiano lungo le direzioni

caratteristiche.

Variazione delle variabili d’onda :

Dove con si rappresenta la variazione della variabile d’onda che viaggia nello spazio fisico lungo la

direzione ed l’ i_esimo autovettore sinistro

Determinazione degli autovettori sinistri:

:autovettore sinistro associato all’ autovalore

[

]

le ultime due equazioni dicono la stessa cosa,

20

scelgo il valore

e per sostituzione si ottiene:

:autovettore sinistro associato all’ autovalore

[

]

le ultime due equazioni dicono la stessa cosa,

scelgo il valore

e per sostituzione si ottiene:

, :autovettori sinistri associati all’ autovalore di molteplicità 2.

21

[

]

;

:

:

per la determinazione di si può scegliere: {

per la determinazione di si può scegliere: {

22

{

Proprietà iperboliche relative al sistema di Eulero 2-D non stazionario per un fluido costituito da due

specie di gas nel caso particolare di fronte d’onda piano:

Formulazione quasi lineare:

{

( 1.4.39 )

(

)

( 1.4.40 )

[ ]

( 1.4.41 )

Calcolo degli autovalori

Ottenimento dell’ equazione caratteristica:

23

( )

[

]

Vista la presenza nella matrice di molti zeri è conveniente al fine di ottenere il determinante applicare

due volte il criterio di Laplace:

L’equazione caratteristica ha per radici, i seguenti autovalori:

; autovalore di molteplicità 3

( 1.4.42 )

Variazione della variabile d’onda :

Dove rappresenta la variazione della variabile d’onda che viaggia lungo la direzione

Determinazione degli autovettori sinistri

:autovettore sinistro associato all’ autovalore

[ ]

24

le ultime due dicono la stessa cosa,

similmente a quanto svolto nel caso 1-D si scegli

e per sostituzione nel sistema sarà possibile ottenere:

:autovettore sinistro associato all’ autovalore

Procedendo analogamente a quanto fatto per si ottiene:

[ ]

Scegliendo

Autovettori relativi all’ autovalore di molteplicità 3:

25

[ ]

Si possono assegnare tre combinazioni di valori purché in grado di generare autovettori

indipendenti tra loro.

Per si può scegliere :

, ,

Per si può scegliere :

, ,

Per si può scegliere :

, ,

È pertanto possibile riscrivere il sistema di Eulero in termini della variazione delle variabili d’onda:

26

{

Equazioni di compatibilità:

[

]

{

Che proiettate lungo le direzioni caratteristiche, restituiscono:

lungo la direzione

27

lungo la direzione:

lungo la direzione

Pertanto è stato messo in evidenza il come le proprietà iperboliche dei casi 1-D e 2-D nel caso

particolare di onde piane siano strutturalmente analoghe.

Un importante risultato restituito dalle equazioni di compatibilità è il fatto che lungo la direzione

“traiettoria della particella” siano trasportate rispettivamente le concentrazioni e la

componente della velocità (normale alla direzione del flusso).

Riepilogando:

lungo le direzioni caratteristiche viaggiano le perturbazioni acustiche in

accordo con l’equazione di compatibilità

lungo la direzione caratteristica viaggia l’ onda entropica e vengono trasportate le

concentrazioni come indicato dalle equazioni ,

, ed avviene il trasporto della componente della velocità , come indicato da dv=0

28

1.5 Relazioni di salto

Anche in questa analisi si parte dallo studio del sistema di equazioni di bilancio che regola l’evoluzione

di un flusso monodimensionale e monocomponente per poi estendere la trattazione al caso di una

onda piana per un fluido composto da una miscela di gas.

Caso 1D per un fluido costituito da un solo gas

Si parte dal sistema ( 1.1.1 )

( )

dove i vettori delle “quantità conservate” e dei “flussi generalizzati” in campo monodimensionale

sono rispettivamente:

( )

Si introducono le seguenti notazioni:

salto della generica grandezza attraverso la discontinuità:

[ ]

media tra le grandezze attraverso la discontinuità:

Sistema delle relazioni di salto:

[ ] [ ]

[ ] [ ] ( 1.5.53 )

[ ] [ ]

Nelle quali con si indica la velocità della discontinuità.

È possibile quindi definire la velocità del fluido relativamente alla discontinuità attraverso la scrittura:

( 1.5.54 )

29

Sviluppo delle relazioni di salto:

l’equazione di salto associata all’equazione di bilancio della massa per mezzo di una trasformazione di

tipo galileiano ( ) restituisce:

[ ]

[ ] [ ]

passaggio possibile visto che e quindi si è potuto portare all’interno della parentesi

quadra la w.

L’introduzione di una trasformazione di tipo galileiana è utile anche per sviluppare l’equazione di

salto associata all’ equazione di bilancio della quantità di moto, pertanto analogamente a quanto

fatto per la relazione di salto associata all’equazione di conservazione della massa, si porta la w

all’interno della parentesi quadra e si ottiene:

[ ] [ ]

[ ] [ ]

[ ] [ ]

Se all’interno della parentesi quadra si introduce la seguente quantità nulla definita come [ ]

(e cioè la relazione di salto associata all’equazione di conservazione della massa) è possibile ottenere:

[ ] [ ]

[ ] [ ]

[ ] [ ] [ ]

E considerando ancora una volta la relazione di salto associata all’equazione di conservazione della

massa si ottiene:

[ ] [ ]

30

Equazione di salto associata all’equazione di conservazione dell’energia:

[ ]

Si ottiene quindi il seguente sistema di relazioni di salto valido appunto per il caso 1_D

monocomponente:

{

[ ] [ ]

[ ] [ ] [ ]

( 1.5.55 )

Relazioni di salto nel caso monodimensionale per una miscela di gas

Per affrontare il caso di un fluido multispecie occorre aggiungere alla trattazione appena svolta

k relazioni di salto, in termini di densità parziali, associate all’equazione di conservazione della massa.

[ ]

che può scriversi:

[ ]

[ ] [ ]

e considerando la relazione di salto associata all’equazione di conservazione della massa [ ] :

[ ]

Resta ora da affrontare lo studio relativo alla relazione di salto associata all’ equazione di

conservazione dell’ energia. Nel caso di fluido multispecie è conveniente riscrivere tale relazione di

salto in funzione di e di grandezza adimensionale definita in funzione di “gamma” nel seguente

modo:

Si parte pertanto dalla:

31

[ ] [ ]

nella quale si sostituiscono le seguenti definizioni valide rispettivamente per l’energia interna per

unità di volume, e per l’entalpia interna per unità di volume:

(avendo utilizzato nelle precedenti la definizione di )

[

] [

]

e andando avanti con i passaggi si ottiene:

[ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ]

Pertanto il sistema di relazioni di salto per un flusso monodimensionale multispecie è il seguente:

{

[ ]

[ ] [ ]

[ ] [ ] [ ] [ ] [ ]

( 1.5.56 )

Relazioni di salto nel caso in cui si considerino “onde piane”, per un fluido multispecie:

Bisogna in questo caso aggiungere al sistema (1.5.56) la relazione di salto associata alla conservazione

della quantità di moto in direzione y, e tener presente la diversa formulazione assunta dall’equazione

32

di conservazione dell’energia visto che il termine riguardante l’energia cinetica deve tener conto

anche della componente v della velocità:

o relazione di salto associata all’equazione di conservazione della quantità di moto in direzione

y:

[ ] [ ]

[ ]

e dalla definizione ( 1.5.54 ) :

[ ]

e sfruttando la relazione di salto associata all’equazione di conservazione della massa, scritta nella

forma [ ] :

[ ] [ ]

si ottiene :

[ ]

o relazione di salto associata all’equazione di conservazione dell’energia:

Si parte dalla relazione di salto associata all’energia relativa al sistema (1.5.53) e la si sviluppa nel

seguente modo:

[ ] [ ]

( 1.5.57 )

Ora si esprimere l’energia totale per unità di volume separando i contributi relativi all’energia

interna e all’energia cinetica

essendo:

33

e sostituendo nella relazione di salto ( 1.5.57 ) si ottiene:

[

] [

]

[ ] [ ] (

) [ ] (

) [ ] [

]

[

] (

) [ ] (

) [ ] [

] [

] [ ]

[ ] [ ] (

) [ ] (

) [ ] [

]

[

] (

) [ ] (

) [ ] [

] [

] [ ]

[ ] [ ] (

) [ ] [

] (

) [ ]

[

] [ ]

[ ] [ ] [ ] ((

) (

)) [

]

[

] [ ]

e per la relazione di salto della massa [ ]=[ ] ,

[ ] [ ] [

] [

] [ ]

34

( 1.5.58 )

e tenendo presente che :

[

] [ ]

[

] [ ]

è possibile riscrivere la ( 1.5.58 ) nel seguente modo :

[ ] [ ] [ ] [ ] [ ] [ ]

Ora è possibile mettere in evidenza all’interno della relazione di salto dell’energia i termini relativi alle

relazioni di salto associate alle equazioni di conservazione delle quantità di moto in direzione x ed y:

[ ] [ ] [ ] [ ] [ ] [ ] [ ]

dove si ricorda che : [ ] [ ] ; [ ] [ ] [ ] [ ] [ ]

Pertanto la relazione di salto associata all’equazione di conservazione dell’energia nel caso di un

fluido multispecie per un onda piana assume la seguente forma:

[ ] [ ] [ ]

( 1.5.59 )

Che poi è la stessa relazione di salto associata all’energia ottenuta nel caso monodimensionale per

una miscela di gas.

E quindi sviluppandola come fatto in precedenza si ottiene:

[ ] [ ] [ ] [ ]

e per la

35

[ ] [ ] [ ]

Infine:

[ ] [ ] [ ] [ ]

( 1.5.60 )

Riepilogando, il sistema delle relazioni di salto necessarie allo studio di un fronte piano per una

miscela di gas è il seguente:

{

[ ]

[ ] [ ]

[ ]

[ ] [ ] [ ] [ ]

( 1.5.61 )

Quindi il sistema di relazioni di salto da considerare per lo studio di onde piane nel caso in cui il fluido

sia composto da una miscela di due gas differenti è il seguente:

{

[ ]

[ ]

[ ] [ ]

[ ]

[ ] [ ] [ ] [ ]

( 1.5.62 )

A partire dal sistema ( 1.5.62 ) si passa ora a studiare la forma che assumono le relazioni di salto nei

casi in cui si considerino rispettivamente discontinuità di contatto ed urti.

Relazioni di salto per una discontinuità di contatto non stazionaria ,relativamente al caso particolare

di onda piana:

le condizioni da imporre sono:

[ ]

[ ]

E per la definizione dove rappresenta la velocità della discontinuità risultano

conseguentemente valide le scritture:

[ ] [ ]

36

essendo nel caso della d. di c.

dalle due relazioni di salto associate alla conservazione della massa espresse in termini di

densità parziali, si ottiene :

[ ]

[ ]

Tenendo presente che per la d. di c. :

sono possibili salti delle concentrazioni e cioè:

[ ] [ ]

Relazione di salto associata alla quantità di moto in direzione x:

[ ] [ ]

nel caso di una d. di c. risulta identicamente verificata.

Relazione di salto associata alla quantità di moto in direzione y:

[ ]

per la condizione risulta,

sono possibili salti della componente della velocità in direzione y,

cioè: [ ]

E quindi qualora si fosse in presenza di un salto non nullo della componente della velocità

trasversale alla discontinuità, la d. di c. può essere considerata anche come una “linea di

scorrimento”.

37

R azi n i sa ass cia a a ’ n ia a :

[ ] [ ] [ ] [ ]

Considerando che per la d. di c. [ ] [ ] :

[ ]

E considerando che e :

attraverso la discontinuità di contatto sono possibili salti della

quindi: [ ]

Cioè essendo la sinonimo di , l’equazione di salto dell’energia nel caso di d. di c. indica che,

attraverso la discontinuità di contatto è possibile che si verifichi un salto in termini di variazione di

composizione.

Relazioni di salto per un urto non stazionario ,relativamente al caso particolare di onda piana:

[ ]

[ ]

relazioni di salto associate alla massa espresse in termini di densità parziali:

[ ]

[ ]

[ ]

[ ]

38

attraverso un urto non sono possibili salti sulle concentrazioni, ed analogamente per la

(1.2.23) è possibile affermare che risulti nullo anche il [ ] . Si può concludere che attraverso

un urto non sono possibili variazioni di composizione.

Relazione di salto associata alla conservazione della quantità di moto in direzione x:

[ ] [ ]

[ ] [ ]

[ ]

e quindi tale relazione fornisce solamente un legame tra il salto di velocità ed il salto di pressione

a av s ’

La relazione di salto associata alla conservazione della quantità di moto in direzione y:

[ ]

[ ]

Attraverso un urto la componente della velocità tangenziale alla discontinuità rimane invariata.

La azi n i sa ass cia a a a c ns vazi n ’ n ia a :

[ ] [ ] [ ] [ ]

In base a quanto riscontrato dalle relazioni di salto associate alla massa si è dimostrato che

attraverso un urto non sono possibili variazioni di concentrazioni ( non può esserci variazione

di composizione ):

ciò vuol dire, ricordando la definizione (1.2.23) che attraverso un urto non sono possibili salti dei

calori specifici e quindi di

39

Si può pertanto affermare che attraverso un urto non sono possibili variazioni di composizione.

Ess n a av s ’ [ ] [ ] ; , la relazione di salto relativa all’energia

assume la seguente forma:

[ ] [ ] [ ]

[ ] [ ] [ ]

cioè si è ritrovata la stessa relazione di salto valida ricavata nel caso di un fluido monocomponente.

A completamento della trattazione sulle relazioni di salto, vengono introdotte tre diverse riscritture

della relazione di salto dell’energia, ottenute esprimendo l’energia interna per unità di volume

rispettivamente nelle seguenti tre forme equivalenti:

o (1) energia interna per unità di volume espressa in termini di

o (2) “ “ “ “ “ “ in termini di

o (3) “ “ “ “ “ “ in termini di

per le quali le quantità e sono delle grandezze costitutive della miscela definite nel seguente

modo:

Si parte dalla relazione di salto,

[ ] [ ] [ ]

Per la quale, di volta in volta l’energia interna per unità di volume verrà espressa in termini di

pressione di temperatura e di velocità del suono (al quadrato).

Come primo punto si passa ad esprimere l’energia interna in termini di pressione:

[ ] [ ] [ ]

40

[ ] [ ] [ ] [ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ] [ ]

(

) [ ] [ ] [ ] [ ] [ ] [ ] [ ]

(

) [ ] [ ] [ ] [ ] [ ]

Def: =

Tale relazione di salto è riscrivibile nel sistema di due equazioni che devono essere soddisfatte

contemporaneamente:

{(

) [ ]

[ ] [ ] [ ] [ ]

Caso di una discontinuità di contatto ( [ ] [ ] ) :

la seconda relazione del sistema è identicamente nulla mentre la prima diventa:

[ ]

Caso di un urto:

come anticipato in questo caso [ ] ,pertanto è la prima relazione del sottosistema ad essere

identicamente nulla mentre la seconda diventa:

[ ] [ ] [ ] [ ]

Si è pertanto evidenziato come l’equazione di conservazione dell’energia, assuma due diverse

formulazioni, rispettivamente nei casi in cui sia presente una di discontinuità di contatto o un urto :

[ ] [ ] [ ]

[ ] [ ] [ ] [ ]

41

Rispetto a quanto restituito nel caso di discontinuità di contatto si riconosce che la

[ ] [ ] [ ] ,

rappresenta a azi n i sa ass cia a a ’ q azi n i bi anci :

Ta q azi n a s n a i “trasporto di eta” ed è da evidenziare che nel caso di una

discontinuità di contatto la dovrà essere considerata costante.

La [ ] ò an ss vis a c m a c n izi n ch a ’an am n i

salti di relativamente al caso di una discontinuità di contatto. Tale condizione verrà di seguito

utilizzata anche in campo discreto al fine di regolare correttamente il trasporto di “ a” tra una

c a ’a a n cas in c i sia s n na isc n in i i c n a .

Come secondo punto si passa ad esprimere l’energia interna in termini della Temperatura:

[ ] [ ] [ ]

[ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ]

{ [ ] [ ]

[ ] [ ] [ ]

Se si ha una discontinuità di contatto che presenta un salto nullo in termini di temperatura:

[ ] [ ] [ ] ,

42

resta non nulla solo la prima, la quale essendo [ ] potrà essere riscritta come:

[ ] [ ]

Ed essendo per la discontinuità di contatto [ ]

Attraverso una discontinuità di contatto possono verificarsi salti di “csi” diversi da zero.

Come terzo punto si passa ad esprimere l’energia interna in termini della velocità del suono

(al quadrato).

[ ] [ ] [ ]

[ ] [ ] [ ]

[ ] [ ] [ ] [ ] [ ]

{ [ ] [ ]

[ ] [ ] [ ]

Se si ha una discontinuità di contatto che presenti un salto nullo in termini di velocità del suono:

[ ] [ ] [ ] ,

resta non nulla solo la prima, la quale essendo [ ] potrà essere riscritta come:

[ ] [ ]

[ ]

43

Con le equazioni ( 1.5.65, 66, 67 ) si sono trovate le condizioni esatte da imporre, in presenza di

discontinuità di contatto, rispettivamente nei seguenti casi particolari :

o [ ]

o [ ]

Tali relazioni verranno introdotte anche nel discreto, quali condizioni necessarie da imporre per

ottenere dei trasporti corretti in termini di “eta”, “csi” e “zita” . Verrà pertanto dimostrato come solo

attraverso tali condizioni risulti possibile annullare eventuali oscillazioni spurie in termini di pressione,

temperatura e velocità del suono, nel caso in cui si sia in presenza di una discontinuità di contatto.

44

Capitolo 2

Equazioni di Eulero nel discreto.

L’obiettivo di questo capitolo è la trattazione di modelli numerici in grado di simulare l’evoluzione di

flussi bidimensionali non stazionari per miscele di gas. Per una soluzione di tipo numerico si andranno

a considerare, valutazioni approssimate al posto degli integrali analitici e differenze finite al posto

delle derivate nel tempo.

Al fine d’integrare numericamente le equazioni di bilancio, tra i vari metodi agli elementi finiti, si è

scelto di utilizzare il metodo di Godunov. Tale metodo riscrive nel discreto le equazioni di

conservazione a partire dalla loro formulazione integrale, pertanto costituisce una forma generale per

risolverne il sistema di Eulero.

Il fatto di riscrivere nel discreto le equazioni di Eulero. partendo dalla loro formulazione integrale

conferisce al metodo di G. ma più in generale a tutti i metodi del tipo “shock capturing” , la capacità di

poter gestire l’evoluzione temporale di flussi in presenza di discontinuità.

In questa trattazione si procede con l’introduzione del metodo di G. del primo ordine in campo

monodimensionale per poi introdurre le modifiche necessarie ad estenderne l’applicabilità anche in

campo bidimensionale.

2.1 Metodo di Godunov

formulazione monodimensionale per un fluido costituito da una miscela di gas perfetti.

Per prima cosa si procede con l’integrare, in un intervallo finito, il sistema di Eulero (1.1.1) nel quale le

equazioni ci conservazione sono espresse in forma di divergenza:

( )

si ricorda che i vettori ed a seconda che si consideri un fluido composto da un so lo gas perfetto o

da una miscela di gas perfetti, sono definiti nel seguente modo:

caso monodimensionale per un fluido ad una sola specie

45

( )

caso monodimensionale per un fluido composto dalla miscela di due diverse specie di gas,

(modello Termodinamico Standard):

( )

Discretizzazione del dominio di definizione spazio temporale

o Il dominio d’integrazione spaziale viene suddiviso in n elementi discreti di uguale ampiezza

pari a

o Per la suddivisione del dominio d’integrazione temporale si utilizzano invece elementi la

cui ampiezza, come verrà dimostrato, dipenderà di volta in volta (ad ogni passo d’integrazione

temporale) dalla soluzione stessa, e da un opportuno coefficiente imposto dalla condizione di

stabilità.

Il volume d’integrazione finito è rappresentato dal seguente intervallo spazio temporale:

[

] [

]

A questo punto si procede con l’integrazione della (1.1.1) nell’intervallo appena definito:

∫ ∫ ( ( ) )

Si introducono la definizione di media spaziale e di media temporale delle grandezze ed ( )

rispetto agli intervalli d’ integrazione,

( (

))

pertanto l’integrale in esame assume la seguente forma:

{

} {

}

46

Si introducono ora le seguenti notazioni:

avendo considerato che per uno schema del primo ordine è corretto scrivere:

,

e quindi

Inoltre è comodo rappresentare con il rapporto tra le ampiezze degli elementi discreti :

Pertanto il sistema di Eulero ( ) integrato nell’intervallo discreto assume la seguente

formulazione,

, la quale nel caso in cui si consideri un modello del primo ordine, può essere scritta come:

( 2.1.1 )

, ed analogamente per esteso :

(

)

( 2.1.2 )

2.2 Metodo di Godunov, in campo bidimensionale:

Come discusso nel corso del precedente capitolo il problema dell’evoluzione di flussi in campo

bidimensionale è governato dal seguente sistema di Eulero

( 2.2.3 a)

L’ evoluzione nello spazio tempo del sistema ( 2.2.3 a), come anticipato nel paragrafo 1.3 , sarà

ottenuta approssimando il fenomeno come una sovrapposizione di onde piane.

47

Questo risulta possibile sulla base della linearizzazione introdotta, per la quale i coefficienti delle

matrici Jacobiane vengono considerati costanti nell’ intervallo discreto (spazio_temporale di volume

)

Cioè nel discreto la soluzione del sistema ( 2.2.3 a ) sarà vista come sovrapposizione delle soluzioni

prodotte dai due fronti d’onda piani, trattati l’uno indipendentemente dall’altro :

La precedente scrittura mette in evidenza il fatto che d’ora in poi si lavorerà separatamente lungo la

direzione x sfruttando le proprietà iperboliche relative alla matrice Jacobiana e lungo la

direzione y con le proprietà iperboliche relative alla matrice .

L’intervallo d’integrazione è rappresentato dal seguente dominio spazio temporale:

[

] [

] [

]

Il sistema ( 2.2.3 a) espresso in forma di divergenza nel caso di un onda piana assume la seguente

forma:

b

Per il quale con si indica il vettore delle quantità conservate, e con ( ) e ( ) i vettori dei flussi

generalizzati, rispettivamente in direzione x e y.

( )

( )

Se con gli indici e si indicano le coordinate discrete della generica cella di calcolo, il relativo

intervallo d’integrazione spaziale ha per estremi (

) e (

)

Si può a questo punto integrare il sistema b rispetto alla generica cella di calcolo:

48

∫

(

∫

)

(

∫

∫

)

(

∫

∫

)

Utilizzando la seguente definizione delle grandezze mediate per cella,

∫

(

∫

)

si assa a in a i sis ma is a ’in va m a i am i zza compreso tra gli

istanti di tempo discreti ed

∫ (

(

) (

))

(

)

(

)

( 2.2.4 a)

Si è quindi riscritto nel discreto quanto anticipato nel capitolo 2. La soluzione in termini di valor medio

per la generica cella al tempo è determinata dalla sommatoria tra il valor medio al tempo

e la sovrapposizione di quantità legate ai flussi all’interfaccia rispettivamente in direzione x ed y (da

ricavare per mezzo della soluzione di problemi di Riemann).

Si riporta di seguito l’espressione in forma compatta del sistema di Eulero in campo 2D riscritto nel

discreto sfruttando il metodo di Godunov

(2.2.4 b)

tenendo presente che ;

49

Scelta del :

Condizione di stabilità secondo il criterio di Courant Friedrics Levi :

si parte da una valutazione della stabilità in campo monodimensionale per poi estenderla al caso

bidimensionale:

dove rappresenta la massima pendenza delle direzioni caratteristiche valutata tra

tutte le celle di calcolo al tempo n .

o Caso 1-D ma | |

j=1,…,N

o Caso 2-D ma | |

⋃ma | |

j=1,…,M ; i=1,…,N

Con si rappresenta il coefficiente C.F.L. (Courant Friedrics Levi).

Al fine di ottenere la stabilità della soluzione si dovrà porre:

o Caso 1-D c (2.2.5 a)

o Caso 2-D c

condizione di Runghe Kutta (2.2.5 b)

Per semplicità di notazione nel seguito della trattazione discreta si utilizzerà in luogo di .

Definizione del modello Termodinamico Standard nel discreto

È possibile ora definire con completezza il modello T.S. in campo discreto:

il sistema di equazioni di bilancio assume la forma seguente:

o monodimensionale per una miscela composta da due gas:

{

(

)

(

)

(2.2.6 a)

50

o bidimensionale ( onda piana in direzione x) per una miscela composta da due gas:

{

(

)

(

)

(2.2.6 b)

o bidimensionale ottenuto come sovrapposizione di onde piane per una miscela di due gas :

{

(

) (

)

(

)

(2.2.6 c)

Nel seguito della trattazione si studierà in primis il sistema (2.2.6 a) e quindi il problema dal punto di

vista monodimensionale, infine i risultati saranno ottenuti conseguentemente anche per i sistemi

(2.2.6 b) e (2.2.6 c), con semplici considerazioni aggiuntive.

È utile introdurre le seguenti notazioni

o def: media aritmetica temporale relativamente alla cella j:

( 2.2.7 a)

o def: media aritmetica spaziale nella cella j al tempo :

( 2.2.7 b)

51

o def: derivata temporale nella cella j al tempo :

( 2.2.8 a)

o def: derivata spaziale nella cella j al tempo :

( 2.2.8 b)

Conseguentemente risultano valide anche le seguenti scritture:

( 2.2.9 a)

Nella quale per semplicità di notazione si è posto :

( 2.2.9 b)

Pertanto è possibile sviluppare le equazioni di bilancio del sistema (2.2.6 a) nel seguente modo:

equazione di bilancio della massa:

Equazione di bilancio della quantità di moto:

che sfruttando l’equazione di bilancio della massa: ( )

si può riscrivere nella forma:

52

Equazione di bilancio dell’energia totale per unità di volume:

per la quale si esprimono per esteso i termini relativi all’energia e all’entalpia totale per unità di

volume

(

) (

)

e ricordando la definizione della variabile “eta” ,

:

(

) (

)

Si sono ottenute quindi le equazioni che rappresentano il modello T.S. discreto nel caso 1D :

{

(

) (

)

(2.2.10 a )

, e similmente il modello T.S. discreto relativo al fenomeno di propagazione di un’onda piana:

{

(

) ( (

))

(2.2.10 b)

53

2.3 Studio di una discontinuità di contatto non stazionaria isolata con il modello T.S.

In questo paragrafo l’attenzione è rivolta allo studio dell’evoluzione di una discontinuità di contatto

non stazionaria per un flusso composto da una miscela di gas. Si passerà pertanto a dimostrare come

in generale il modello T.S. nella sua formulazione discreta, ottenuta con un metodo numerico di tipo

“shock capturing”, già al primo passo d’integrazione temporale, restituisca soluzioni affette da errore

in termini di pressione .

Caso 1D ci si riferisce pertanto al sistema ( 2.2.10 a ):

Una discontinuità di contatto non stazionaria impone le seguenti condizioni:

Le condizioni iniziali in termini di vettore delle variabili di stato nell’ipotesi che la d. di c. al tempo

sia collocata tra la generica cella e la cella adiacente sono le seguenti:

_esima cella a sinistra della cella (compresa)

_esima cella a destra della cella (compresa)

Si vuole ricavare la soluzione restituita dal modello T.S. espresso in base alla( 2.2.10 a ) a seguito di un

passo d’integrazione temporale. Come cella d’indagine si considera quella “etichettata” con

cioè quella alla destra della discontinuità di contatto sui dati iniziali.

Per prima cosa si vede cosa restituiscono le due equazioni di conservazione della massa espresse in

termini di densità parziali:

o Specie 1: c n izi n im s a a a i c

Si ottiene infine:

( 2.3.11 a)

54

o Specie 2: c n izi n im s a a a i c

Si ottiene infine:

(2.3.11 b)

A questo punto è possibile esprimere la soluzione in termini di densità di miscela:

ottenendo pertanto:

(2.3.11 c)

che poi rappresenta la soluzione al tempo nella cella restituita dall’equazione di bilancio:

Si è ora in grado di definire le concentrazioni strumento indispensabile per la chiusura del sistema di

Eulero discreto il momento in cui si scelga di utilizzare il modello T.S.

per definizione la concentrazione delle k-esima specie è data da:

Ed utilizzando le ( 3.1.6 a ) ,( 3.1.6 b ), è facile ricavare:

(2.3.12 a)

(2.3.12 b)

55

o Equazione di bilancio della quantità di moto:

Ricordando le notazioni introdotte ( 3.1.6 a ) e ( 3.1.6 b ), (

e

) , si

ottiene una pratica riscrittura dell’equazione di bilancio della quantità di moto :

(

)

per la condizione imposta dalla discontinuità di contatto si ha che , quindi:

(

) (

)

Che dal punto di vista fisico può essere soddisfatta se si impone,

(2.3.13)

o l’equazione di bilancio dell’energia partendo dalla definizione riportata nel sistema ( 3.1.8 a)

diventa:

(

) (

)

(

)

(

)

e tenendo presenti le (2.3.13 ) e le condizioni imposte dalla discontinuità di contatto, si

ottiene:

e per il bilancio della massa relativo alla miscela ottenuto nella

é utile introdurre la seguente quantità chiamata “sigma” , definita come:

56

e quindi per sostituzione si ottiene:

(2.3.14)

Riassumendo il sistema T.S. che regola la soluzione numerica per una miscela di gas nel caso di una

discontinuità di contatto è il seguente:

{

La ( 2.3.13 ) assicura l’assenza di oscillazioni spurie della variabile di stato rappresentante la velocità.

In generale in base alla ( 2.3.14 ) non si verifica lo stesso per l’altra variabile di stato, la pressione.

A tal fine è utile la riscrittura di quest’ultima in modo da esprimerla in termini di errore relativo:

partendo proprio dalla (2.3.14 ) ed introducendo la

si ottiene:

(

) (

)

Allo scopo di ricavare un’espressione dell’errore relativo si procede nel seguente modo:

E considerando che nel caso di una d. di c. è vera la seguente scrittura , si ottiene:

In quest’ultima si riconosce l’equazione del trasporto di :

57

(2.3.15)

È noto dalla ( 1.2.23 ) come per il metodo Termodinamico Standard l’espressione di “gamma” e quindi

conseguentemente anche di si ricavi a partire dalle concentrazioni, pertanto ricorrendo alle

(1.2.21 ), (1.2.22 ) e considerando che si ottiene:

∑ (

)

∑ ( )

espressione di eta che sostituita nella ( 3.2.13 ) :

( )

E quindi:

Per maggior praticità conviene riportarsi in termini di :

e per sostituzione si ottiene infine:

(2.3.16)

58

È la ( 3.2.16) l’equazione alla quale si fa riferimento per condurre lo studio sull’errore relativo della

pressione.

Sfruttando quest’ultima infatti è facile verificare come in generale il modello T.S. restituisca valori di

diversi da zero e di come pertanto in generale si ottenga nel caso di una d. di c. non

stazionaria:

In completo disaccordo con il fatto che la pressione nel caso della d. di c. debba conservare il suo

valore costante nel tempo in tutto il campo ( e pertanto anche a ridosso della discontinuità )

Esistono tuttavia dei casi particolari di discontinuità di contatto nei quali il modello T.S. attraverso la

( 3.2.16 ) restituisce risultati corretti in termini di pressione.

Condizioni sulla sigma:

o

tale condizione ha senso solo se ottenuta per e cioè nel caso in cui il flusso sia stazionario

o

cioè si va ad imporre un vincolo sull’ incremento tale che nella cella alla destra della

discontinuità di contatto sia presente la sola specie di gas imposta dalle condizioni iniziali di destra

ed analogo discorso a sinistra. Tale considerazione può essere tradotta facilmente in termini di

concentrazioni scrivendo:

e quindi:

( )

( )

In generale la discontinuità di contatto per un metodo di tipo “shock capturing” è catturata e non

inseguita, e nella cella j+1 (adiacente alla d. di c. ) al tempo n+1 si avrà sia che

59

Caso di un fluido monocomponente

Nel modello T.S. la è ricavabile per mezzo dei calori specifici nel seguente modo:

∑ (

)

∑ ( )

e considerando che per un fluido mono_componente:

e per sostituzione si ottiene:

Salto nullo sulla temperatura attraverso la discontinuità di contatto:

Considerando l’equazione di stato per i gas perfetti ( ), ed il fatto che per una

discontinuità di contatto =0, imporre equivale a imporre

Ovviamente è analogamente vero il viceversa:

Riepilogando, in generale la condizione che deve essere rispettata per non avere oscillazioni spurie

della pressione è la seguente:

Se questa condizione si sostituisce direttamente nella (2.3.14) si ottiene :

ed essendo in questo caso ed ovviamente anche , si ottiene semplicemente:

(2.3.17)

si è pertanto ricavato nel discreto un’equazione che rappresenta il trasporto della variabile .

60

La ( 3.2.17 )è la condizione che si deve imporre se non si vogliono avere errori in termini di oscillazioni

spurie della variabile di stato pressoria. Si osserva che tale equazione di trasporto è compatibile con la

relazione di salto ricavata alla fine del capitolo precedente.

Pertanto per non aver nel discreto errori in termini di pressione, la “eta” ( e quindi la “gamma” ) dovrà

essere ricavata sfruttando la (2.3.17) e non il criterio standard (che utilizza le concentrazioni).

Motivo per il quale non è possibile applicare il metodo standard per ricavare nel discreto una

soluzione corretta in termini di gamma, ma più in generale, in termini delle quantità costitutive di

miscela.

In realtà per il modello T.S. ,se si utilizza un metodo numerico di tipo “shock capturing”,in generale le

grandezze derivate non sono ben calcolate.

Il problema riguarda la forma assunta dalle grandezze di miscela ,ed in particolare del “gamma”

ricavato in funzione delle concentrazioni. Pertanto si può affermare che, pur verificandosi una certa

monotonicità sulle , questa non assicura la condizione necessaria per ottenere nel discreto il

corretto andamento del rapporto dei calori specifici necessario alla corretta valutazione della

pressione.

Caso di continuo per soluzioni regolari

Si considera l’equazione di conservazione della massa relativa alla k-esima specie:

Sostituendo e sviluppando:

Che per l’equazione di conservazione della massa :

e quindi:

Ora si passa a valutare le relazioni di salto:

[ ] [ ]

61

[ ] [ ]

[ ] [ ] [ ] [ ]

e considerando la relazione di salto associata all’equazione di conservazione della massa:

[ ] [ ]

Nel caso di una discontinuità di contatto è possibile scrivere in luogo di :

[ ]

quindi il salto delle concentrazioni non dipende da .

Si procede ora col mettere in evidenza l’equazione che governa il trasporto delle per il modello T.S.

in campo discreto.

sviluppando,

Si ricorda che,

,

allora è possibile scrivere :

per la quale risulta nulla

La quale non costituisce un semplice trasporto e soprattutto non rappresenta il modo utilizzato dal

modello T.S. per ricavare le concentrazioni nel discreto.

62

Discontinuità di contatto o linea di scorrimento:

Un’osservazione aggiuntiva necessaria alla descrizione di una discontinuità di contatto nel caso in cui

si consideri un’onda piana, riguarda il salto della velocità trasversale.

o se è corretto parlare di discontinuità di contatto;

o se è invece più corretto parlare di “linea di scorrimento”

2.4 Modelli discreti corretti

Considerati gli errori restituiti dal modello T.S. in termini di pressione (al primo passo d’integrazione

temporale) nasce l’esigenza di formulare nuovi modelli.

A questo punto si sfrutta quanto trattato in materia di relazioni di salto, ( paragrafo 1.5), al fine di

introdurre uno o più modelli in grado di ricavare correttamente le grandezze derivate.

Modelli Eta Conservativi:

Si è dimostrato nel paragrafo 1.5 come la relazione di salto associata all’equazione di conservazione

dell’energia, possa essere riscritta attraverso il seguente sottosistema, per il quale le due equazioni

devono essere soddisfatte contemporaneamente:

{(

) [ ]

[ ] [ ] [ ] [ ]

Tale sistema andrà ad assumere, nel caso di una discontinuità di contatto, la seguente forma

specifica:

[ ]

Si vuole pertanto ,a partire dall’equazione di bilancio dell’energia in campo discreto, ottenere una

formulazione in grado di far emergere la condizione ottenuta nell’ambito delle relazioni di salto:

63

(

) (

)

(

) ( (

) )

tale equazione di bilancio pertanto sarà soddisfatta se sono nulli separatamente i seguenti termini :

(2.4.18)

(

) ( (

))

(2.4.19)

Ottenendo il seguente sottosistema:

{

(

) ( (

))

Da osservare che il seguente sistema nei casi di urto e di discontinuità di contatto sarà caratterizzato

da una sola equazione:

o infatti nel caso di un urto la prima equazione sarà identicamente nulla;

o nel caso di una d. di c. la prima sarà proprio l’equazione di trasporto cercata

, mentre sarà la seconda ad essere identicamente nulla.

La precedente riscrittura dell’equazione di conservazione dell’energia consente di definire unitamente

alle altre equazioni di bilancio il modello denominato Eta_Conservativo_Implicito.

64

Eta_Conservativo_Implicito:

o caso 1-D ,miscela di due gas:

{

(

)

(

)

(

) ( (

))

(2.4.20 a)

o Caso 2-D “onda piana in direzione x” per la miscela di due gas:

{

(

)

(

)

(

) ( (

))

(2.4.20 b)

Il modello in questione è di tipo implicito in quanto le ultime due equazioni di bilancio rappresentanti

l’equazione di conservazione dell’energia per unità di volume, sono poste tra loro in forma implicita.

Infatti dipende da

e viceversa dipende da

.

Introducendo nella terza e nella quarta equazione del sistema (2.4.20 a) o nella quarta e nella quinta

del sistema (2.4.20 b) le seguenti approssimazioni:

e cioè andando a confondere la media spaziale delle grandezze con la media temporale, si riesce a

disaccoppiare le ultime due equazioni, ottenendo così una versione esplicita del modello E.C.I. che

prenderà appunto il nome di modello Eta_Conservativo_Esplicito:

65

Eta_Conservativo_Esplicito

caso 1-D:

{

(

)

(

)

(

) ( (

))

(2.4.21 a)

Altro modello d’interesse ottenuto anch’esso a partire dal E.C.I. è quello proposto da Abgrall e Karni,

imponendo la seguente approssimazione

Abgrall_Karni

{

(

)

(

)

(2.4.22 a)

Modello coincidente al E.C.E. ad eccezione dell’equazione di bilancio dell’energia che qui è riportata

nella sua formulazione standard.

Si nota inoltre come la struttura del modello A.K. sia sostanzialmente la stessa del modello T.S. ma in

questo caso la chiusura viene effettuata per mezzo della determinata grazie all’equazione di

trasporto:

Per una discontinuità di contatto e la assume la forma:

e quindi la condizione che si deve imporre se non si vogliono avere oscillazioni spurie di pressione.

66

Un’ importante considerazione da fare riguarda la chiusura dei sistemi appena introdotti. Come

discusso precedentemente per il modello T.S. (sia nella formulazione continua che discreta), la

chiusura avviene mediante l’equazione di stato per esteso, dove per ricavare il rapporto dei calori

specifici si utilizza la ( 1.2.23 ) ,ed è quindi necessario ricavare le concentrazioni e quindi utilizzare

tante equazioni di bilancio della massa (riferite alle densità parziali) quante sono le specie di gas.

Nel caso dei modelli “discreti corretti” appena introdotti il rapporto dei calori specifici si ricava in

funzione di ottenuta per mezzo della sua equazione di evoluzione imposta dal soddisfacimento

dell’equazione di bilancio dell’energia.

Pertanto per lo studio di fluidi costituiti da miscele di gas nei modelli “eta conservativi” , non è

necessaria la stesura di tante equazioni di conservazione della massa (in termini di densità parziali)

quante le specie di gas presenti.

Altro discorso riguarda il calcolo della temperatura, in effetti per i modelli E.C.I. ,E.C.E., ed A.K. è

necessario sviluppare una procedura a parte, che consenta di determinare la costante R di miscela.

A tale scopo è necessario aggiungere al sistema di equazioni l’equazione di trasporto del peso

molecolare in forma non conservativa:

, e conseguentemente in termini discreti:

(2.4.24)

Una volta ricavato il peso molecolare della miscela è possibile ricavare la costante R della miscela

proprio utilizzando la definizione:

Quindi la struttura assunta dai modelli “eta conservativi” è la seguente :

caso 1-D per una miscela di gas:

Eta_Conservativo_Implicito ( con determinazione della temperatura)

{

(

) ( (

))

(2.4.25)

67

Eta_Conservativo_Esplicito ( con determinazione della temperatura)

{

(

) ( (

))

( 2.4.26 )

Abgrall_Karni ( con determinazione della temperatura)

{

(2.4.27)

Caso di un’onda piana:

Eta_Conservativo_Implicito “onda piana” ( con determinazione della temperatura)

{

(

) ( (

))

(2.4.28)

Eta_Conservativo_Esplicito “onda piana” ( con determinazione della temperatura)

{

(

) ( (

))

(2.4.29)

Abgrall_Karni “onda piana” ( con determinazione della temperatura)

68

{

(2.4.30)

Studio di una discontinuità di contatto isolata con i modelli E.C.I. , E.C.E. , A.K.

Le condizioni da imporre nei sistemi ( 2.4.25 ), ( 2.4.26 ), ( 2.4.27 ) sono:

Equivalenza dei metodi E.C.I. , E.C.E. , A.K. . per una d. di c.

Si comincia col ricavare la struttura assunta dal modello Eta_Conservativo_Implicito nel caso in cui si

voglia studiare una discontinuità di contatto isolata e non stazionaria:

l’equazione di conservazione della massa diventa,

Per l’ equazione di conservazione della quantità di moto, si procede con il seguante sviluppo:

e sfruttando l’equazione di conservazione della massa , si ottiene:

ed imponendo le condizioni e

69

Infine ricordando che

e che

, si ottiene:

(

)

e quindi:

Ponendo , , nella terza equazione del sistema ( 2.4.25 ) si ottiene:

(

) (

)

e ricordando l’ equazione di conservazione della massa, , si ottiene:

La quarta equazione ,considerando le condizioni imposte dalla discontinuità di

contatto , , ,

L’ultima equazione, rappresentante il trasporto del peso molecolare resta invariata.

Il modello E.C.I. nel caso di una discontinuità di contatto isolata e non stazionaria assume la forma:

70

{

ed utilizzando la notazione:

{

Il modello E.C.E. differisce dal modello E.C.I. per la terza e la quarta equazione di bilancio,

di seguito riportate:

(

) ( (

))

e analogamente a quanto svolto per il modello E.C.I. :

imponendo , ,

la terza restituisce e di conseguenza:

Mentre la quarta per il fatto che sono nulle le seguenti variazioni , ,

potendo scrivere:

71

Si è quindi dimostrato come il modello E.C.E. coincida con il modello E.C.I. nel caso di una

discontinuità di contatto.

Modello A.K.:

anche il modello A.K. differisce dal E.C.I. per la formulazione della terza e della quarta equazione di

bilancio ma anche in quest’ultimo per una discontinuità di contatto assume la forma ( 2.4.31 )

Da osservare il come nei modelli E.C.I. ed E.C.E. la condizione sia verificata direttamente dalla

terza equazione di bilancio, indipendentemente dall’ equazione imposta per determinare il parametro

di miscela.

Nel modello A.K. invece la condizione è subordinata alla soluzione della quarta equazione di

conservazione essendo la terza essenzialmente l’equazione di conservazione dell’energia nella sua

forma standard.

Concludendo tutti e tre i modelli nel caso di una discontinuità di contatto non stazionaria sono

rappresentati dal sistema ( 2.4.31)

{

72

3.5 Modello G.A.K. :

Si sfruttano i seguenti tre modi tra loro equivalenti di esprimere l’energia interna per unità di volume:

( 2.5.32 )

Dove si sono introdotte le seguenti definizioni:

( 2.5.33 )

Si procede similmente a quanto fatto per il metodo E.C.I. , e quindi partendo dall’equazione di bilancio

dell’energia espressa nella sua formulazione standard, si andranno di volta in volta a sostituire le

(2.5.32) .

Si ricorda che l’espressione dell’equazione di conservazione dell’energia in campo discreto e nella

formulazione standard è la seguente:

Si ricorda inoltre che:

, e quindi per sostituzione:

(

) (

)

a questo punto si procede col sostituire al posto di le equivalenti quantità riportate nelle ( 2.5.32 )

73

Sostituendo :

si ottiene quanto già descritto in precedenza, e cioè il fatto che l’equazione di conservazione

dell’energia è soddisfatta se sono verificate separatamente le due equazioni del sistema:

{ (

) ( (

))

( 2.5.34 a )

,e sfruttando l’approssimazione proposta da A.K.

, la seconda equazione del sistema restituisce l’equazione del trasporto di eta

( 2.5.34 b )

Sostituendo si ottiene:

(

) (

)

Questa volta l’obiettivo è quello di ottenere l’ equazione di evoluzione di in forma

conservativa:

(

) (

)

Che risulta verificata se è verificato il sistema:

{

(

) (

)

( 2.5.35 a )

e nell’ ipotesi che sia , cioè introducendo un’approssimazione che confonde la media

spaziale di con la media temporale, la prima equazione del sottosistema diventa :

74

( 2.5.35 b )

, che rappresenta quanto cercato e cioè l’equazione di evoluzione di

Sostituendo si ottiene:

(

) (( )

)

E sviluppando anche in questo caso con l’idea di ricavare un’ equazione di evoluzione per di

tipo conservativo :

(

) (

)

Che risulta verificata se entrambe le equazioni del sistema sono verificate:

{

(

) (

)

( 2.5.36 a )

anche in questo caso ipotizzando si ottiene per la prima:

( 2.5.36 b )

Anche in questo caso un’ equazione di evoluzione.

Formulazione di tipo non conservativo (valido solo nel caso di discontinuità di contatto):

Andando a sostituire si ottiene:

(

) (

)

questa volta l’obiettivo è quello di ottenere l’ equazione di evoluzione di in forma non conservativa:

75

(

) (

)

Che risulta verificata se è verificato il sistema:

{

(

) (

)

( 2.5.37 a )

e similmente a quanto ipotizzato da A.K.

si ottiene per la prima:

( 2.5.37 b )

E quindi un’ equazione del trasporto per la

Andando a sostituire si ottiene:

(

) (

)

sviluppando:

(

) (

)

Che risulta verificata se è soddisfatto il sistema:

{

(

) (

)

e con l’ipotesi :

76

ottenendo un’equazione di trasporto per la .

Per lo studio della soluzione in termini di e di la strada corretta da intraprendere è quella di

utilizzare le equazioni in forma conservativa in quanto non avrebbe senso in generale trattare la

evoluzione temporale di queste due grandezze utilizzando una forma non conservativa e quindi per

mezzo di un semplice trasporto come fatto per la . In effetti sia la che la sono funzione oltre che

dei parametri costitutivi della miscela anche della densità ; , e quindi

un’equazione di trasporto non riuscirebbe a simulare correttamente la soluzione nel caso in cui si

abbia a che fare con problemi di Riemann. Le quantità costitutive in effetti attraverso un’ espansione

o un urto non variano, ma la densità si ( eq. di Poisson e di compatibilità per l’espansione e R.

Hugoniot per l’urto ).

Ovviamente per quanto anticipato si può facilmente dimostrare come per lo studio di una

discontinuità di contatto non stazionaria ed isolata ( assenza di espansioni e/o di urti) le due

formulazioni ( conservativa e non conservativa ) siano perfettamente equivalenti.

Si può ora passare a definire il seguente modello termodinamico, denominato G.A.K. (Abgrall Karni

generalizzato) :

{

( 2.5.39 a )

Per completezza si riporta anche la versione nella quale le equazioni di evoluzione di e di hanno

forma non conservativa, valida solo per lo studio di una discontinuità di contatto non stazionaria.

{

( 2.5.39 b )

77

Modello G.A.K. e chiusura del sistema di equazioni di Eulero

Tale modello è di notevole importanza in quanto dimostra come in realtà anche il sistema di Eulero

nel discreto per una miscela di gas se correttamente formulato sia chiuso e il come oltretutto non si

debbano aggiungere altre equazioni di trasporto per la determinazione della temperatura, come

risulta invece necessario per i modelli eta conservativi precedentemente studiati.

Con il modello G.A.K. è stato trovato un modo per chiudere il sistema di Eulero ( rispetto alle variabili

) utilizzando le sole informazioni contenute nelle equazioni di bilancio della massa, della

quantità di moto e dell’ energia totale.

Si nota che entrambi i sistemi ( 2.5.39 a ) e ( 2.5.39 b ) presentano una ridondanza legata

all’introduzione dell’equazione di evoluzione della variabile “zita” , quantità utilizzata per ricavare

correttamente la soluzione in termini di velocità del suono.

In effetti per garantire la chiusura e quindi la soluzione del sistema di Eulero in

sarebbero sufficienti le sole prime quattro equazioni e quindi:

{

Il precedente sistema infatti non presenta ridondanze, ma non consente di ottenere in generale

soluzioni corrette in termini di velocità del suono.

Il problema della ridondanza può tradursi nel fatto che il sistema ( 2.5.39 a ) “vede” due differenti

modi per definire la velocità del suono:

√

√

Tale questione non comporta problemi dal punto di vista della soluzione numerica, in quanto può

essere aggirata scegliendo a priori di ottenere la velocità del suono utilizzando unicamente la (2) , e

cioè la strada che consente di ottenere una soluzione corretta.

La gestione della ridondanza all’interno del solutore di Riemann, (necessario alla valutazione dei flussi

all’interfaccia) sarà invece argomento del prossimo capitolo (paragrafo 3.2).

78

Equazione di stato per il modello G.A.K. :

(

)

Entrambi i sistemi ( 2.5.39 a ), ( 2.5.39 b ) nel caso di una d. di c. assumono la seguente forma:

{

( 2.5.40 )

Dove .

3.6 Studio di una d. di c. non stazionaria con i modelli E.C.I. , E.C.E., A.K. .

Si procede ora con l’applicare i modelli eta conservativi al caso di una discontinuità di contatto isolata.

La soluzione sarà ricavata al primo passo d’integrazione temporale, e valutata nelle celle adiacenti la

discontinuità al fine di verificare o meno la correttezza della soluzione numerica.

Come dati in input saranno introdotte le seguenti condizioni iniziali (che unitamente alle condizioni

) rappresentano dei casi particolari di discontinuità di contatto:

79

L’interfaccia presa in considerazione sarà la , compresa tra la cella e la cella , mentre

per il passo d’integrazione temporale si considera quello compreso tra gli istanti di tempo n ed

n .

Condizioni iniziali:

celle alla sinistra della discontinuità di contatto e quindi alla sinistra dell’interfaccia :

celle alla destra della discontinuità di contatto e quindi alla destra dell’interfaccia :

Si passa ora a valutare la forma assunta dal sistema ( 2.4.31 ) dopo un passo d’integrazione

temporale.

Per la prima equazione del sistema riguardante il bilancio della massa si procede nel seguente modo:

si considera la densità riferita alla cella j+1 (a destra dell’ interfaccia j+1/2 ) al tempo n+1

(

)

nella quale si introducono le condizioni iniziali ( per n=0) :

ottenendo quindi :

che può anche essere riscritta come :

80

Per la seconda e la terza equazione , si ottiene direttamente:

Mentre la quarta e la quinta procedendo analogamente a quanto fatto per la massa

restituiscono :

Quindi il sistema ( 2.4.31 ) nel caso in cui si voglia ricavare la soluzione nella cella j+1 (a destra

della discontinuità) dopo un passo d’integrazione temporale assume la seguente forma:

{

( 2.6.41 )

Caso a) Sia nullo il salto della temperatura attraverso la discontinuità di contatto:

Considerando l’equazione di stato dei gas perfetti ed il fatto che attraverso la discontinuità di

contatto è nullo il salto di pressione :

La soluzione in termini di temperatura ad un generico istante di tempo può essere

ricavata per mezzo dell’equazione di stato:

81

e nel caso in cui la soluzione si voglia ottenere al tempo n=1

dove come imposto dal sistema ( 2.6.41 ).

La p al tempo è esprimibile (utilizzando l’equazione di stato) in funzione della quantità

e quindi:

per la quale se non si vogliono avere oscillazioni spurie di temperatura bisogna imporre:

e tale condizione sarà verificata se:

Ma il valore di è noto dalla ( 2.6.41 )e sostituendo:

e quindi si ottiene:

82

che può essere riscritta in termini di peso molecolare considerando che :

Si è ricavata l’equazione del trasporto del peso molecolare, contenuta nel sistema ( 2.6.41 ),

dimostrando quindi come i modelli “eta conservativi” presi in esame, siano in grado di simulare

correttamente una d. di c. con distribuzione uniforme della temperatura.

Caso b) Sia nullo il salto dell’energia interna per unità di massa attraverso la discontinuità di contatto:

(

) (

)

Dove ovviamente essendo

(

)

(

)

(

)

(

)

La pressione al tempo n=0 è esprimibile in funzione del rapporto :

(

)

e quindi è possibile scrivere:

(

)

per la quale se non si vogliono avere oscillazioni spurie di energia interna per unità di massa

bisogna imporre:

83

condizione verificata se si impone:

(

)

e moltiplicando a destra e a sinistra per ,

(

)

, per la quale, sostituendo l’espressione ottenuta per , è possibile ricavare:

(

)

, e dividendo a destra e a sinistra per :

( )

si è quindi ottenuta l’equazione del trasporto della “eta”, contenuta nel sistema ( 2.6.41 ).

Pertanto i modelli E.C.I. , E.C.E. ed A.K. simulano correttamente l’evoluzione di una discontinuità di

contato nel caso particolare di salto nullo sull’energia interna per unità di massa.

Caso c) Sia nullo il salto dell’entalpia interna per unità di massa attraverso la discontinuità di contatto:

(

) (

) (

)

e tenendo conto che

(

)

(

)

(

)

(

)

84

Esprimendo p in funzione dell’ entalpia al tempo :

(

)

(

)

Per non avere oscillazioni spurie in termini di entalpia interna per unità di massa si deve imporre

la condizione:

che sarà verificata se:

(

)

(

)

E per la definizione :

(

)

( ) (

)

Ma per il legame che sussiste tra ed tale scrittura equivale all’equazione del trasporto per

contenuta nel sistema (2.6.41):

Anche in questo caso particolare di discontinuità di contatto ( ) i modelli E.C.I. , E.C.E. ed A.K.

restituiscono risultati corretti.

85

Caso d) Sia nullo il salto della velocità del suono attraverso la discontinuità di contatto:

per comodità di esposizione si può studiare analogamente il problema in termini di

(

)

ed essendo

(

) (

)

(

)

(

)

La velocità del suono al quadrato nella cella j+1 al tempo n=1 si può ricavare come:

Ed esprimendo la pressione in funzione di

(

)

(

)

e per escludere errori sulla velocità del suono al tempo :

Che sarà verificata se:

(

)

(

)

Nella quale sostituendo la prima equazione del sistema (2.6.41):

(

)

86

Tale equazione non è in generale verificata nei modelli E.C.I. , E.C.E. ed A.K.

In effetti nel discreto le equazioni di trasporto per “eta” e per “gamma” e quindi le

;

risultano equivalenti solo nei seguenti casi particolari :

o caso stazionario

o

o gas monocomponente

In generale i modelli E.C.I. , E.C.E. ed A.K. non simulano correttamente l’evoluzione di una

discontinuità di contato nel caso particolare di salto nullo sulla velocità del suono.

Caso e) Sia nullo il salto del calore specifico a volume costante attraverso la discontinuità di contatto:

(

)

Calore specifico nella cella j+1 al tempo n=1:

Calore specifico nella cella j+1 al tempo n=0:

Procedendo similmente a quanto fatto per le precedenti dimostrazioni:

ed imponendo

si ottiene :

che conviene esprimere in funzione del peso molecolare,

87

( )

( )

( )

( )

(

)

ed introducendo l’equazione del trasporto per il peso molecolare,

(

)

si ottiene :

Che è l’equazione di evoluzione di “eta” contenuta nel sistema (2.6.41)

Pertanto i modelli E.C.I. , E.C.E. ed A.K. simulano correttamente l’evoluzione di una discontinuità di

contatto nel caso particolare di salto nullo del calore specifico a volume costante attraverso la

discontinuità.

Caso f) Sia nullo il salto del calore specifico a pressione costante attraverso la discontinuità di

contatto:

(

)

Il calore specifico a pressione costante all’istante di tempo n=1 nella cella j+1 è ottenibile come:

Mentre all’istante di tempo n=0:

88

e procedendo similmente al precedente caso:

ed imponendo

si ottiene:

Esprimendo quest’ultima in termini di peso molecolare:

( )

( )

(

)

(

)

che è l’equazione di evoluzione di equivalente all’ equazione di evoluzione di “eta”.

Pertanto i modelli E.C.I. , E.C.E. ed A.K. simulano correttamente l’evoluzione di una discontinuità di

contatto nel caso particolare di salto nullo del calore specifico a pressione costante attraverso la

discontinuità.

Si può concludere affermando che, i tre modelli “eta conservativi” presi in considerazione, simulano

correttamente tutti i casi particolari di discontinuità di contatto analizzati ad eccezione del caso in cui

sia nullo il salto della velocità del suono.

89

2.7 Studio di una discontinuità di contatto utilizzando il modello G.A.K.:

Il modello G.A.K. nel caso di una discontinuità di contatto non stazionaria assume la forma ( 2.5.40 )

{

Si vuole pertanto verificare se il modello simuli correttamente il problema dell’evoluzione di una

discontinuità di contatto sui dati iniziali relativamente ai casi particolari presi in considerazione nel

precedente paragrafo:

Rispetto a quanto svolto nel paragrafo 3.6 ( modelli E.C.I. , E.C.E. , A.K. ) la dimostrazione risulta

analoga nei seguenti casi:

;

Mentre per i restanti casi la dimostrazione è simile dal punto di vista dell’ approccio ma ovviamente

diversa per quanto riguarda gli strumenti utilizzati.

caso a) discontinuità di contatto non stazionaria nel caso in cui si abbia allo stato iniziale

:

(

) (

)

(

)

(

)

( 2.7.42 )

90

L’equazione di stato permette di ricavare la temperatura in funzione di (o analogamente

utilizzando la definizione (2.7.42) in funzione di ).

(

)

( 2.7.44 )

Che valutata al tempo

(

)

( 2.7.45 )

Al tempo iniziale n=0 è possibile ricavare la pressione in funzione della e della (

)

(

)

e per quanto imposto dal sistema ( 2.5.40 )

, e sostituendo nella ( 2.7.44 ) :

(

)

(

)

La condizione che si deve imporre per non avere oscillazioni spurie di temperatura attraverso la

discontinuità di contatto è la seguente :

( )

e per l’equazione di evoluzione di “csi” :

91

( )

( )

( )

si è ottenuto che la condizione nel caso della discontinuità di contatto è retta dall’equazione

di trasporto per “eta”, contenuta nel sistema (3.4.22) (modello G.A.K.) il modello G.A.K. simula

correttamente un discontinuità di contatto non stazionaria di un fluido multispecie nel caso in cui sia

nullo il salto sulla temperatura.

caso b) discontinuità di contatto non stazionaria nel caso in cui si abbia allo stato iniziale

:

(

)

che per diventa:

(

) (

) (

) (

)

(

)

(

)

ed al tempo n=1 :

mentre al tempo n=0 è possibile esprimere la pressione in funzione di ed del rapporto

( )

( )

92

per non avere in output oscillazioni spurie in termini di velocità del suono bisogna imporre :

ne consegue che:

( )

E per l’equazione dell’evoluzione di “zita”

( )

( )

( )

anche in questo caso particolare si è ottenuta l’equazione di trasporto di .

Il modello G.A.K. pertanto, a differenza di quanto avviene per i modelli E.C.I. , E.C.E. , ed A.K., è in

grado si simulare correttamente una discontinuità di contatto anche nel caso particolare di salto della

velocità del suono nullo.

caso c) discontinuità di contatto non stazionaria nel caso in cui si abbia allo stato iniziale

:

(

) (

) (

)

(

)

(

)

Il calore specifico a volume costante nella cella j+1 al tempo n=1 :

( 2.7.46 )

Mentre al tempo n=0 è esprimibile come :

( )

93

e moltiplicando quest’ultima a secondo membro della ( 3.7.28 ) si ottiene:

( )

per non avere oscillazioni spurie in termini di calore specifico a volume costante la condizione da

imporre è la

, si ottiene così :

( )

( )

( )

( )

, che rappresenta l’equazione di conservazione della massa, e quindi si è dimostrato come il modello

G.A.K. simuli correttamente una discontinuità di contatto anche in questo caso particolare.

94

Capitolo 3

Il problema di Riemann in questa trattazione ricopre un ruolo cardine in quanto da una parte

consente di ottenere la soluzione “esatta” (caso 1D o caso di onda piana), dall’altra costituisce uno

strumento necessario allo sviluppo del metodo di Godunov. In effetti come espresso nel corso del

capitolo 2 , la valutazione dei flussi all’interfaccia è subordinata alla soluzione di problemi di Riemann.

3.1 Problema di Riemann nel caso di un’onda piana.

L’algoritmo che governa il problema di Riemann nel caso in cui si voglia studiare l’evoluzione di

un’onda piana (caso particolare di un problema bidimensionale), per un fluido costituito da una

miscela di gas, è strutturalmente analogo a quello utilizzato per lo studio di un flusso

monodimensionale riferito ad una sola specie di gas.

Quanto affermato può essere dimostrato spostando l’attenzione sullo studio, dei sistemi di Eulero che

rappresentano i due problemi, in termini di proprietà iperboliche e di relazioni di salto.

Caso 1D fluido costituito da un solo gas:

{

(

) [

]

Caso di onda piana per un fluido multispecie:

{

(

) (3.1.2)

che può essere riscritta anche nel seguente modo andando a mettere in evidenza il trasporto delle

concentrazioni:

{

(

)

(3.1.3)

95

e quindi per il caso di nostro interesse (miscela di due specie di gas) :

{

(

)

(3.1.4)

La concentrazione della specie 2 è ottenibile dalla (1.2.14)

In riferimento ai sistemi (3.1.1) e (3.1.4) si passa ora ad analizzarne le proprietà iperboliche.

caso 1D :

{

caso 2D (onda piana):

{

Il fatto che entrambi i sistemi siano iperbolici ne consente una riscrittura lungo le direzioni

caratteristiche (ottenimento delle equazioni di compatibilità).

Grazie alla determinazione delle equazioni di compatibilità è ora possibile valutare con facilità le

proprietà che regolano l’espansione.

In effetti si ricorda che il fenomeno dell’espansione in campo 1D non stazionario è governato dalle:

{ q azi n i c m a ibi i n a i zi n ca a is ica

q azi n i iss n

nel caso di un’espansione sulla famiglia di linee caratteristiche

q azi n i c m a ibi i n a i zi n ca a is ica

q azi n i iss n

nel caso di cui l’espansione sia sull’altra famiglia di linee caratteristiche

96

Il caso di un’onda piana per un fluido costituito da due specie di gas differisce dal precedente per

l’aggiunta delle due equazioni di trasporto rispettivamente delle concentrazioni e della velocità

trasversale all’onda, e quindi delle:

Tali equazioni regolano il trasporto delle quantità e lungo la direzione traiettoria della particella,

pertanto e resteranno invariate a prescindere dal fatto che la particella aumenti la sua velocità

(andando da monte a valle dell’ espansione).

Relazioni di salto:

Anche in considerazione di quanto restituito dalle relazioni di salto si evince come i due sistemi siano

sostanzialmente analoghi ad eccezione del fatto che per il caso di un’onda piana si dovranno

considerare due relazioni in più, quella che governa il salto di una delle due concentrazioni e quella

che rappresenta il salto della velocità trasversale. Tale argomento è stato ampiamente discusso nel

capitolo 1, ma per maggiore chiarezza verranno di seguito riportati i risultati in termini [ ] e di [ ] :

Per la k-esima concentrazione si avranno le seguenti condizioni:

discontinuità di contatto [ ]

urto [ ]

Anche per quanto riguarda le relazioni di salto che regolano il comportamento della componente

della velocità trasversale (all’onda), il discorso è essenzialmente lo stesso:

discontinuità di contatto [ ]

urto [ ]

Sulla base di quanto restituito dallo studio delle caratteristiche iperboliche del sistema e delle

relazioni di salto, si può concludere come l’unico caso nel quale sia le concentrazioni che la velocità

trasversale possono variare sia attraverso una discontinuità di contatto.

Considerando la rappresentazione grafica sul piano fisico x,t di un generico tubo d’urto, è possibile

inquadrare per via schematica quanto appena descritto.

97

Concludendo si può affermare come il problema di Riemann relativo all’evoluzione di una

discontinuità sui dati iniziali (costanti a tratti), in riferimento al caso di un’onda piana per un fluido a

due specie di gas, sia chiuso in p,u, e analogamente al caso di un flusso monodimensionale per

un fluido monocomponente.

Algoritmo del “solutore di Riemann” caso 1D:

Come noto il sistema ondoso generato da un generico problema di Riemann, può essere del tipo :

o RCR

o RCS

o SCR

o SCS

Se si considerano i seguenti dati in input la procedura è articolata nel

seguente modo:

98

1. Si determina una stima di prima approssimazione della velocità della discontinuità di contatto:

( (

))

(

)

per la quale:

(

)

{

2. Inizio di un ciclo iterativo:

o Espansione (governata dalla Poisson e dall’equazione di compatibilità lungo la direzione

caratteristica che attraversa l’onda)

o Urto (governato dalle Rankine Hugoniot):

Espansione a sinistra (sulle linee caratteristiche u-a): (procedura attivata se

)

dati in input

dove : mentre è la stessa andando da monte a valle dell’espansione.

eq di compatibilità che attraversa le u-a:

Poisson:

(

)

99

Direzioni sul piano fisico dei fronti d’onda:

Espansione a destra (sulle linee caratteristiche ) (procedura attivata se

)

dati in input

dove : mentre è la stessa andando da monte a valle dell’espansione.

Eq di compatibilità sulle direzioni caratteristiche a che attraversano le a

Poisson:

(

)

Direzioni sul piano fisico dei fronti d’onda:

100

Urto a sinistra: (procedura attivata se

)

dati in input

dove : mentre è la stessa andando da monte a valle dell’urto per quanto imposto dalle

relazioni di salto, (non esiste variazione di composizione attraversando un urto)

Rankine Hugoniot:

√ (

)

√(

)

Direzione sul piano fisico dell’urto:

Urto a destra: (procedura attivata se

)

dati in input

dove : mentre è la stessa andando da monte a valle dell’urto per quanto imposto

dalle relazioni di salto, (non esiste variazione di composizione attraversando un urto)

Rankine Hugoniot:

√ (

)

101

√(

)

Direzione sul piano fisico dell’urto:

3. Aggiornamento della velocità della discontinuità di contatto all’ iterata k+1 (metodo delle tangenti

Newton-Rapson)

4. Convergenza della soluzione e logica di uscita:

dato un piccolo a piacere si uscirà dal ciclo iterativo se:

|

|

Altrimenti si aggiorna l’indice d’iterazione e si ritorna al punto 2.

5. Una volta ultimati i cicli iterativi la soluzione nei settori 2 e 3 è nota in a , è necessario

pertanto al fine di completare la determinazione delle variabili fondamentali ricavare la densità:

per le quali ovviamente risulta:

102

3.2 Solutore di Riemann per il modello G.A.K.

Come trattato nel capitolo 2, il nuovo modello presenta un eccesso di equazioni rispetto al numero di

variabili.

Dal punto di vista numerico il problema non crea ambiguità in quanto si sceglie a priori di ottenere la

soluzione in termini di velocità di propagazione dei disturbi sfruttando l’ultima equazione del sistema

( 2.5.39 a ) e quindi in funzione dell’energia interna per unità di volume e di “zita”.

Lo stesso ragionamento non è sufficientemente valido il momento in cui si debba passare attraverso

l’algoritmo “Solutore di Riemann”. La struttura del solutore di Riemann, in effetti è chiusa rispetto alle

variabili e considera la velocità del suono attraverso la definizione :

√

Si ricorda che per i modelli “eta conservativi”, ma anche per il G.A.K. , il rapporto dei calori specifici è

ricavabile attraverso “eta”:

Per il modello G.A.K. in virtù del fatto che esistono due modi di definire la velocità di propagazione dei

disturbi, si può procedere secondo due tipologie d’inserimento dei dati in input:

o Criterio tradizionale: vengono passate al solutore di Riemann le grandezze necessarie a

chiudere il P.R. :

o Criterio basato sul modello G.A.K. che fornisce al solutore di Riemann altre due variabili di

campo, , e quindi complessivamente:

L’eccesso di equazioni rispetto al numero delle variabili per quanto riguarda il modello G.A.K. , può

tradursi nel fatto che il modello G.A.K. è in grado di vedere la velocità di propagazione delle

perturbazioni in due modi differenti:

√

103

√

Nel primo caso d’inserimento, criterio tradizionale, si trascura a priori la definizione , e si passa

a ridefinire all’interno della procedura solutore di Riemann la velocità del suono a partire dalle

grandezze di campo e quindi come:

√

che poi ovviamente coincide con la definizione ,(quella “scartata” da G.A.K.)

Così operato si è escluso il problema della ridondanza all’interno de solutore di Riemann, e dal punto

di vista della soluzione si può procedere tradizionalmente, percorrendo l’algoritmo riportato nel

precedente paragrafo in quanto vengono utilizzate le sole variabili necessarie e sufficienti alla

chiusura del problema di Riemann.

o Secondo criterio d’inserimento, (metodo G.A.K.) :

si inserisce nel problema di Riemann la definizione (3.2.6) e si sceglie di gestire la ridondanza

all’interno del solutore stesso. La ridondanza è presente in quanto il problema di Riemann è

naturalmente chiuso in grandezze che vedono la velocità di propagazione dei disturbi

in termini della (3.2.5).

Si dimostra come la compresenza di due definizioni della velocità di propagazione dei disturbi,

imponga un’opportuna modifica della procedura risolutiva del solutore di Riemann.

Ragionando per il momento in termini di condizioni iniziali, si procede col ricavare un’espressione

della (3.2.6) in grado di far emergere il rapporto , al fine di poterla mettere a confronto con la

(3.2.5):

√

c T

T

e se si definisce:

, si ottiene l’espressione di cercata:

104

√

√

per la quale le variabili sono le stesse variabili di campo contenute nella (3.2.1).

Dal precedente studio si conclude che il rapporto è esprimibile in due differenti modi a seconda

che si scelga il criterio d’inserimento tradizionale o il criterio suggerito dal modello G.A.K. :

a

Pertanto si è ricavata la condizione che dovrà essere imposta, al termine della procedura iterativa, per

ricavare correttamente la densità:

o per il criterio d’inserimento tradizionale la densità sarà ricavata dalla:

a

o mentre per il criterio d’inserimento indicato dal modello G.A.K. si userà la:

Ora si vuole dimostrare come sia il metodo tradizionale, che quello imposto dal modello G.A.K.

restituiscono soluzioni praticamente sovrapponibili tra loro.

La dimostrazione nel caso di una discontinuità di contatto isolata verrà ottenuta algebricamente,

mentre per il caso più generale (presenza di urti o di d. di c. ) sarà ottenuta attraverso una simulazione

numerica.

o Caso di una discontinuità di contatto:

verranno coinvolte le sole procedure:

“espansione a sinistra” per una discontinuità di contatto che propaga verso sinistra

“espansione a destra” “ “ “ “ che propaga verso destra

Nel caso di una discontinuità di contatto gli stati 2 e 3 coincidono banalmente con gli stati 1 e

4.

105

Pertanto essendo non verrà percorso il ciclo iterativo, e le procedure

“espansione a sinistra” ed “espansione a destra” restituiranno direttamente:

Si procede con il calcolo della densità:

a sinistra della discontinuità di contatto:

a destra della discontinuità di contatto :

Nel caso di una discontinuità di contatto si è pertanto dimostrato come entrambi i criteri

d’inserimento dati, al solutore di Riemann, tradizionale e G.A.K., restituiscano soluzioni coincidenti in

termini di .

o Caso generale del problema di Riemann (presenza di un generico sistema ondoso):

A titolo di esempio si considera il problema relativo ad un’espansione sulla famiglia di linee

caratteristiche a . Tale fenomeno risulta governato dalle:

(

)

(

)

106

Mettendo a confronto i due tipi di inserimento è possibile dimostrare come, in questo caso, la

pressione ottenuta con l’inserimento G.A.K. sia diversa dalla ottenuta con l’inserimento

tradizionale. In effetti:

essendo si otterranno:

(

)

(

)

per le quali ovviamente

Un discorso analogo si manifesta anche per l’espansione sull’altra famiglia di linee caratteristiche ed

ovviamente anche per quanto riguarda il caso di un urto (che impone l’utilizzo delle R.H.)

Una stima di quanto le grandezze si discostino dalle (o viceversa), può essere condotta

direttamente per via numerica.

Le seguenti informazioni sono relative al test SRM per il quale per una discussione più dettagliata si

rimanda al capitolo 5. In questa sede ci si limita solamente a valutare le differenze in termini di

soluzioni, restituite dai due diversi tipi di inserimento di dati nella procedura “solutore di Riemann”.

Ai primi passi d’integrazione temporale i due criteri d’inserimento restituiscono in output valori

praticamente identici. Una certa variazione tra le grandezze diventa apprezzabile, se pur minima, al

10° passo d’integrazione temporale:

Le seguenti tabelle contengono i valori medi per cella delle grandezze di stato

relativamente alle celle investite dalla perturbazione a seguito di 10 passi di integrazione temporale, e

mettono a confronto le soluzioni restituite dai due diversi tipi di l’inserimento dei dati in input al

solutore di Riemann.

Velocità del suono [ ]:

202 203 204 205 206 207 208 209

1224.959 1239.387 1258.212 1294.726 1377.178 1509.264 1436.391 1060.715

1224.962 1239.525 1258.445 1294.434 1376.418 1509.854 1437.852 1059.738

107

Pressione in [ ]

202 203 204 205 206 207 208 209

934155.7 1020856.2 1045766.4 1038018.6 1023251.3 924669.3 472960.8 151927.3

934250.6 1022720.7 1049856.5 1041438.6 1026947.5 930202.9 474486.8 151415.9

Velocità [ ]

202 203 204 205 206 207 208 209

1860.605 1760.453 1727.974 1728.310 1721.678 1624.77 1002.78 195.580

1860.491 1758.401 1724.363 1726.597 1723.182 1631.758 1005.605 193.506

densità

202 203 204 205 206 207 208 209

0.712001 0.741758 0.689300 0.597717 0.529361 0.485573 0.355667 0.225374

0.71204 0.742613 0.690305 0.596803 0.528076 0.486044 0.355846 0.225034

Le seguenti figure rappresentano graficamente i dati riportati nella tabella relativa alla velocità

del suono:

Figura 3.1

108

Figura 3.2

A questo punto si procede con l’analizzare le soluzioni restituite dai due differenti metodi d’inserimento al

tempo t=0.0001534 sec a seguito di 266 passi d’integrazione temporale, al fine di mostrare come

entrambi i metodi, restituiscano in output, andamenti delle grandezze sovrapponibili tra loro.

109

In questa “scala” non sono apprezzabili le differenze tra le grandezze restituite dai due criteri.

110

Figura 3.3

Figura 3.4

111

Figura 3.5

Figura 3.6

112

Figura 4.7

Figura 4.8

113

Si può concludere che entrambi i metodi sono in grado di restituire in output grandezze tra loro

praticamente sovrapponibili, anche in presenza di problemi di Riemann.

114

Capitolo 5

Simulazione numerica di problemi di Riemann

Test 5.1

problema di SOD

Modello Termodinamico Standard

Per dominio di definizione spaziale si è considerato un condotto di lunghezza pari ad un metro e

diametro unitario, discretizzato in 400 elementi finiti. Allo stato iniziale, a metà della lunghezza del

condotto è posto un setto fittizio separante i seguenti stati assunti da un unico gas:

Assegnazione dei dati iniziali:

o

o

o

o

La simulazione è stata eseguita utilizzando il modello T.S. , mentre la soluzione è stata ricavata al

tempo t=0.2 sec, considerando un coefficiente di Courant pari a 0.8

115

Figura 5.1

Figura 5.2

116

Figura 5.3

Figura 5.4

117

Figura 5.5

Figura 5.6

118

Si nota come nel caso di un fluido composto da una sola specie di gas, il modello T.S. sia in grado di

simulare correttamente l’evoluzione nello spazio tempo di un generico tubo d’urto.

Test 5.2

Test SRM:

Si assume come dominio di definizione spaziale un tubo di lunghezza pari ad un metro e diametro

unitario, all’interno del quale allo stato iniziale sono presenti due gas ideali in differenti condizioni

termodinamiche, separati da un setto fittizio posto a metà lunghezza. La soluzione è stata ricavata al

tempo t=0.0001534 sec

dati iniziali:

o o o o o o

Modello T.S. : Rosso

Modello A.K. :Blu

Modello G.A.K. :Verde

Soluzione esatta : linea tratteggiata

119

120

121

122

123

124

125

Nei grafici riportati di seguito è messo in evidenza il particolare degli andamenti delle variabili pressione,

velocità, densità, temperatura ed energia interna per unità di volume nell’intorno della discontinuità di

contatto.

126

127

128

129

Test 5.3

(Test H)

dati iniziali:

o

o o o

o

Si assume come dominio di definizione spaziale un condotto avente lunghezza pari ad un metro e

diametro unitario, all’interno del quale allo stato iniziale sono presenti due gas in differenti condizioni

e separati da un setto ideale posto a metà lunghezza.

La soluzione è stata ricavata al tempo t = 0.205 sec, sfruttando un coefficiente di Courant pari a 0.8,

mentre il dominio spaziale è stato discretizzato in 400 elementi finiti

I modelli utilizzati sono, il Termodinamico Standard, l’ Abgrall-Karni e l’ Abgrall-Karni generalizzato.

Modello T.S. :Rosso

Modello A.K. :Blu

Modello A.K._G. :Verde

Soluzione esatta : linea tratteggiata

130

131

132

133

134

135