Gasdinamica 2 Progetto del corso -...

POLITECNICO DI MILANO

Facolta di Ingegneria Industriale

Corso di Laurea in Ingegneria Aeronautica

Gasdinamica 2Progetto del corso

Silvia Adami 682281Angelo Colbertaldo 681338

Dario Isola 682855

Anno Accademico 2006/07

Indice

Introduzione v

1 Il contesto teorico 11.1 Ipotesi globali . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Fluidodinamica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Termodinamica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.4 Termo-fluidodinamica . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2 I modelli matematici 72.1 Modello quasi–monodimensionale . . . . . . . . . . . . . . . . . . . . . . 72.2 Metodo di Sauer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2.1 Caso assisimmetrico . . . . . . . . . . . . . . . . . . . . . . . . . 132.3 Il metodo di Hall-Kliegel assisimmetrico . . . . . . . . . . . . . . . . . . 142.4 Metodo delle caratteristiche . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.4.1 Caso assisimmetrico . . . . . . . . . . . . . . . . . . . . . . . . . 18

3 Algoritmica 213.1 Metodo quasi–monodimensionale . . . . . . . . . . . . . . . . . . . . . . 21

3.1.1 Soluzione subsonica: Pa ≥ P1 . . . . . . . . . . . . . . . . . . . . 223.1.2 Ugello adattato: Pa = P2 . . . . . . . . . . . . . . . . . . . . . . 233.1.3 Urto normale nel divergente: P3 ≤ Pa < P1 . . . . . . . . . . . . 233.1.4 Ugello sovraespanso: P2 < Pa < P3 . . . . . . . . . . . . . . . . . 243.1.5 Ugello sottoespanso: Pa < P2 . . . . . . . . . . . . . . . . . . . . 24

3.2 Metodo delle caratteristiche . . . . . . . . . . . . . . . . . . . . . . . . . 24

4 Risultati 334.1 Ugelli cartesiani . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 344.2 Ugelli assisimmetrici . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 514.3 Osservazioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

A Integrazione sulla sezione 77

B Ugello con curve di Bezier 79

iii

iv INDICE

Introduzione

Il problema piu classico della gasdinamica, relativo a flussi interni, e quello di risolverela corrente che si sviluppa all’interno di un ugello convergente-divergente e in letteraturae possibile trovare diversi metodi analitici e numerici che consentano di conoscere effet-tivamente le caratteristiche di questi flussi. Con l’intento di semplificare il modello dipartenza e impostare un’algoritmo numerico che ne permetta l’ effettiva risoluzione, talimetodi introducono approssimazioni e ipotesi dalla cui effettiva realizzazione dipendea priori la validita del metodo. Obiettivo di questo laboratorio e quello di analizzaree verificare le ipotesi introdotte dai vari modelli e, una volta implementati, mettere aconfronto i risultati ottenuti per mezzo di quantita che siano effettivamente comparabilie indicative di qualita.

I metodi che verranno analizzati, per risolvere la corrente comprimibile all’internodi un ugello convergente-divergente, sono: quasi monodimensionale, di Sauer e dellecaratteristiche.

Nel capitolo 1 si chiarisce il contesto che fa da sfondo a tutti i capitoli successivi, siadal punto di vista fluidodinamico sia termodinamico con la scrittura delle equazioni diEulero comprimibile. Nel capitolo 2 queste sono specializzate per i casi relativi, a par-tire dalla semplificazione spaziale: 1-D, 2-D cartesiano o assisimmetrico. Nel capitolo 3sono invece presentati gli algoritmi effettivi di risoluzione, con particolare attenzione almetodo delle caratteristiche, l’unico che effettivamente presenta reali complicazioni im-plementative. Nel capitolo 4 sono invece messi a confronto i risultati prodotti, andandoin particolare a osservare l’andamento del numero di Mach lungo il condotto anche alvariare della geometria locale dell’ugello.

Questo lavoro e scritto in LATEX.

v

vi INTRODUZIONE

Capitolo 1

Il contesto teorico

In questo capitolo si descrive in maniera generale, anche se brevemente, il modello ma-tematico che permette di rappresentare la dinamica di un flusso all’interno di un ugelloconvergente divergente. In questo tipo di condizione reale si assiste ad una interazionetra le grandezze meccaniche e quelle termodinamiche proprie del fluido, pertanto unadescrizione puramente fluidodinamica non e sufficiente, ma deve essere necessariamenteaccoppiata ad una che consideri anche la termodinamica dei gas.

Come ben noto, le equazioni che permettono di descrivere in maniera completa ilcomportamento di un gas sono quelle di Navier–Stokes comprimibili. Esse sono equa-zioni iperboliche in tempo e ellittiche nello spazio e tengono conto sia degli scambienergetici tra quantita meccaniche e termiche, sia della presenza all’interno dei fluididella viscosita, con la nascita di sforzi di taglio. Volendo descrivere in maniera sempli-ficata il comportamento di un gas, e possibile trascurare tali sforzi, almeno in alcunezone del campo di moto (spesso quelle lontane dalle pareti), andando cosı a scrivere leequazioni di Eulero comprimibili. Questo tipo di semplificazione comporta la perditadel termine di derivata seconda, lasciando intatto quello advettivo e, conseguentemente,alla modifica della natura delle equazioni, che diventano paraboliche in tempo e lega-no il loro comportamento nello spazio al numero di Mach. Questa modifica puo nellamaggioranza dei casi, portare a problemi piu complicati di quello di partenza.

Storicamente si sono affermati metodi analitici e numerici che permettessero di ri-solvere le equazioni di Eulero, tramite l’introduzione di ipotesi specifiche per alcuni pro-blemi della gasdinamica. Di seguito verranno introdotte le ipotesi generali che permet-tono di scrivere le equazioni di Eulero specifiche per un ugello convergente-divergente;nel capitolo successivo saranno invece descritti i metodi risolutivi e le considerazionispecifiche.

1.1 Ipotesi globali

Vengono di seguito riportate le ipotesi globali valide per tutti i metodi presentati:

Ipotesi del continuo: il fluido-gas e rappresentato tramite il modello del continuo.L’equilibrio meccanico e termodinamico e riferito all’elemento infinitesimo di con-

1

2 Capitolo 1. Il contesto teorico

tinuo, che costituisce un sistema aperto a scambi di massa, quantita di moto,energia, entropia.

Flusso stazionario: nessuna dipendenza dal tempo delle variabili che descrivono lameccanica del fluido.

Fluido in equilibrio termodinamico locale: i processi termodinamici sono infinita-mente lenti (quasistatici), rappresentando il passaggio tra continui stati di equi-librio. Cio implica anche che le due quantita intensive che descrivono il gas, peresempio P e T , sono sempre definite.

Termodinamica piu rapida della fluidodinamica: i processi fluidodinamici sonoinfinitamente piu lenti di quelli termodinamici che, all’interno dell’elemento infini-tesimo di fluido, non causano alcun transitorio, ma raggiungono istantaneamentel’equilibrio. Il fluido e sempre in equilibrio termodinamico locale.

Serbatoi di capacita infinita: l’ugello convergente–divergente (De-Laval) e connessoad un serbatoio di massa, energia e quantita di moto di capacita infinita, siaall’ingresso sia all’uscita (ambiente).

Fluido non viscoso: i blobs fondamentali di fluido sono in grado di trasferire solosforzi normali alla loro superficie e gli sforzi di taglio sono trascurabili in tutto ilcampo di moto; si ipotizza direttamente ν = 0.

1.2 Fluidodinamica

La meccanica del continuo permette di determinare il moto di un solido o di un fluidotramite le equazioni di conservazione. Le leggi di conservazione principali sono: massa,quantita di moto ed energia. Da sole sono sufficienti a descrivere i sistemi meccanici,governati da processi reversibili. Riferite all’elemento infinitesimo di fluido, sotto leipotesi presentate, assumono la forma:

∇ · (ρu) = 0

∇ · (ρu⊗ u + P I) = 0

∇ · ((ρet + P )u)

= 0

(1.1)

con et = e + 12 |u|2 e completato dalle necessarie condizioni al contorno.

La (1.1) forma un sistema di 5 equazioni per 6 incognite: ρ(x), e(x) , P (x) e u(x).Se si decidesse di scindere la fluidodinamica dalla termodinamica basterebbe rendere ladensita una costante del problema, cosı da separare le prime due equazioni dalla terza,ottenendo le equazioni di Eulero incomprimibile. Per trattare i due aspetti del problemain maniera accoppiata, e necessario introdurre gli elementi base della termodinamica eaggiungere almeno una equazione al sistema, per ottenere un numero di equazioni paria quello delle incognite.

1.3. Termodinamica 3

1.3 Termodinamica

Per poter scrivere un sistema termo-meccanico completo, e necessario modellare le tra-sformazioni di energia introducendo il concetto di irreversibilita, ovvero l’esistenza diuna direzione obbligata in alcuni processi fisici. Un sistema termodinamico puo esseredescritto in maniera esaustiva dalla cosı detta equazione fondamentale: e una equazionescalare, a due varaibili indipendenti e relativa a quantita estensive. Due sue possibiliformulazioni, specifiche alla massa, sono quella energetica e = e(s, v) e quella entropicas = s(e, v). Quest’ultima gode delle seguenti proprieta:

• omogenea di ordine 1 rispetto a e e v,

• strettamente monotona rispetto a e,

• superadditiva.

La superadditiva e l’omogeneita implicano la convessita

s(λe1 + µe2, λv1 + µv2) ≥ λs(e1, v1) + µs(e2, v2) con λ + µ = 1

che permette di rappresentare le trasformazioni fisiche irreversibili. La conoscenza del-l’equazione fondamentale, unita ad un’altra relazione detta equazione di Eulero, forniscetutte le informazioni necessarie a rappresentare la termodinamica della materia.

La scelta delle variabili indipendenti presenta un certo grado di arbitrarieta. Appli-cando la trasformazione di Legendre, e infatti possibile ottenere diverse rappresentazionidell’equazione fondamentale, in modo da manipolarla al fine di ottenere la forma chemeglio si adatta al problema preso in esame, senza perdita alcuna di informazione. Eper questo motivo che e possibile incontrare diverse formulazioni: energetica, entropica,entalpica, energia libera di Gibbs, funzione di Helmoltz . . .

Come inserito nelle ipotesi iniziali, la termodinamica e sempre all’equilibrio. Questaaffermazioni si riflette immeditatamente sulla nostra possibilita di conoscere le proprietadel fluido. Solo all’equilibrio e infatti possibile definire le seguenti variabili intensive:

T = T (s, v) =(

∂e

∂s

)

v

, P = P (s, v) = −(

∂e

∂v

)

s

.

La conoscenza di una sola delle equazioni di stato fra T (s, v) e P (s, v), che puo essereil risultato di una analisi empirica o di un modello di gas, ci permette di conoscere s o ein funzione delle sole variabili T, v o P, v. Questo risultato, purtroppo, non e sufficientea descrivere il sistema termodinamico in maniera completa, poiche, nell’operazione diderivazione che permette di definire le equazioni di stato, si incappa necessariamente inuna perdita di informazione rispetto all’equazione fondamentale che, per poter essererappresentata in maniera completamente equivalente, ci obbliga a conoscere entrambele equazioni di stato.


Il modello piu semplice di gas e quello di gas ideale e politropico (cv = cost), le cuiequazioni di stato sono:

P (T, ρ) = ρRT (1.2)

e(T ) =R

γ − 1(T − T0), (1.3)

le quali ci permettono di conoscere l’equazione fondamentale a meno di una costante:

s(e, ρ) = s0 + Rln

[(e

e0

) 1γ−1 ρ0

ρ

]. (1.4)

I modelli termo-fluidodinamici che verranno descritti nei capitoli seguenti necessita-no della definizione di un modello di gas, e possono rappresentare gas reali, di Van DerWaals, ecc. . . La costruzione del modello di gas ideale e politropico comporta notevolivantaggi in termini pratici, infatti molte delle equazioni possono essere semplificate ealcune risolte addirittura in forma chiusa. Questa scelta presenta anche dei ben preci-si limiti di valitida fisica che schemi piu complessi e generali non possiedono, ma talenecessita non e dovuta alla scarsa potenza della formulazione matematica, che e perfet-tamente in grado di descrivere il comportamento di un qualsiasi gas reale, ma e figliadei limiti dell’intelligenza (sebbene ne simboleggi contemporaneamente la straordinariacapacita esplorativa).

1.4 Termo-fluidodinamica

Per accoppiare i due campi principali che dirigono la dinamica dei gas, e necessarioaggiungere le equazioni di stato, e.g. la (1.2) e la (1.3), oppure direttamente l’equazionefondamentale, e.g (1.4), alla (1.1).

∇ · (ρu) = 0

∇ · (ρu⊗ u + P I) = 0

∇ · ((ρet + P )u)

= 0

P = P (e, ρ)

T = T (e, ρ)

(1.5)

1.4. Termo-fluidodinamica 5

nel caso di gas ideale politropico, introducendo l’entalpia totale ht = h + Pρ + 1

2 |u|2,otteniamo:

∇ · (ρu) = 0

∇ · (ρu⊗ u + P I) = 0

∇ · (ρhtu) = 0

P = ρRT

h = cpT

(1.6)

Nel presente lavoro, si utilizzara sempre l’ipotesi di gas ideale politropico, le cui sempli-ficazioni rientrano a diversi livelli nei modelli matematici presentati.

Capitolo 2

I modelli matematici

In questo capitolo vengono introdotti gli schemi oggetto di confronto di questo labora-torio, vengono enunciate le ipotesi che sono alla base e scritte le equazioni risolte daimetodi.

2.1 Modello quasi–monodimensionale

Questo metodo si basa su una serie di ipotesi molto restrittive, che lo rendono il piu sem-plice tra i metodi esposti in questo lavoro, ma anche quello che meglio si presta a trattarele varie condizioni di funzionamento a cui si trova ad operare un ugello convergente–divergente. Infatti, dato che la natura delle equazioni che costituiscono il metodo noncambia al variare del numero di Mach, come nel caso 2-D, con l’approssimazione diflusso quasi–monodimensionale siamo in grado di rappresentare tutte e tre le condizionisubsoniche, soniche e supersoniche; questa possibilita lo rende l’unico metodo, tra quelliproposti in questo lavoro, in grado di descrivere il flusso lungo tutto l’ugello, di tratta-re urti normali all’interno del convergente e di considerare anche flussi completamentesubsonici.

Le ipotesi che vengono introdotte a priori sono:

• approssimazione quasi–monodimensionale della geometria:

- A(x) ∈ C0

-∣∣∣ 1√

AdAdx

∣∣∣ ¿ 1,

- hr ¿ 1, con h altezza del condotto ed r raggio di curvatura;

• assenza di attriti mondimensionali (τw = 0);

• flusso adiabatico (q = 0);

Tali ipotesi permetto di semplificare notevolmente il problema in esame. La prima ipo-tesi permette di ritenere costanti le grandezze termodinamiche e la velocita sulla sezione,consentendo quindi di passare da un sistema tridimensionale ad uno monodimensionale,

7

8 Capitolo 2. I modelli matematici

dove cio che conta sono gli integrali delle grandezze o piuttosto il valore medio di taligrandezze sulla sezione, definiamo H

.=∣∣∣ 1√

AdAdx

∣∣∣. La seconda e la terza ipotesi permet-tono di ritenere il flusso omoentropico ed omoentalpico, escludendo la possibilita dellapresenza di urti che verranno trattati in maniera speciale. Le grandezze sono dunquetutte continue.

Il sistema di equazioni che governa la dinamica della corrente diventa quindi:

d

dx(ρuA) = 0

d

dx

[(ρu2 + P )A

]= P

dA

dxd

dx

(h + 1

2u2)

= 0

equazione di stato

(2.1)

che, specializzato per un gas ideale politropico, diventa

d

dx(ρuA) = 0

d

dx

[(ρu2 + P )A

]= P

dA

dxd

dx

(h + 1

2u2)

= 0

P = ρRT

h = cpT

(2.2)

un sistema chiuso che puo essere risolto analiticamente tramite l’introduzione del numerodi Mach e la scrittura dell relazioni isoentropiche.

Trattandosi di un flusso adiabatico non viscoso, la variazione di velocita all’inter-no dell’ugello dipende solamente dalla variazione di area; sfruttando la conservazionedell’entalpia totale e possibile esprimere la velocita in funzione del numero di Mach, equindi, combinando le due espressioni, ottenere un legame tra il numero di Mach e lageometria dell’ugello:

A(M)A∗

=1M

[2

γ + 1

(1 +

γ − 12

M2

)] γ+12(γ−1)

, (2.3)

dove A∗ rappresenta l’area della sezione critica, cioe quella in cui si raggiungono condi-zioni soniche (M = 1).Grazie all’omoentropicita del flusso e quindi possibile, conoscendo l’andamento del nu-mero di Mach lungo l’ugello, ottenere le grandezze termodinamiche dalla conservazione

2.2. Metodo di Sauer 9

delle grandezze totali, note in corrispondenza del serbatoio:

c20

c2= 1 +

γ − 12

M2,

T0

T= 1 +

γ − 12

M2,

ρ0

ρ=

(1 +

γ − 12

M2

) 1γ−1

,

P0

P=

(1 +

γ − 12

M2

) γγ−1

.

(2.4)

L’estrema semplicita del modello utilizzato consente di trattare anche condizioni difunzionamento che presentano delle discontinutita nell’andamento delle grandezze ter-modinamiche e quindi studiare anche correnti supersoniche in presenza di urti normali.Infatti, riscrivendo il sistema (2.1) in forma integrale e ritenendo le correnti a monte eda valle dell’urto omoentropiche, e possibile ottenere le relazioni di salto delle grandezzetermodinamiche in funzione del numero di Mach a monte dell’urto. Utilizzando quindil’equazione di continuita e le relazioni (2.3) e (2.4) per la corrente a valle dell’urto e pos-sibile ricavare il numero di Mach a monte e quindi, sempre usando la (2.3), la posizionedell’urto.

2.2 Metodo di Sauer

Il secondo metodo che viene preso in esame in questo lavoro e il metodo di Sauer per lasoluzione delle equazioni del potenziale comprimibile. Sotto le seguenti ipotesi:

• quantita che descrivono il flusso continue,

• flusso omoentalpico,

• flusso omoentropico,

• la funzione che descrive l’ugello y = η(x) deve, in gola, soddisfare la seguente:

η(xt) ∈ C2.

Il sistema delle equazioni di Eulero 2D (1.5) puo essere scritto nella seguente forma:

12u(x, y) · ∇|u(x, y)|2 − c2(|u|)∇u(x, y) = 0

∇× u(x, y) = 0

s(x, y) = s

ht(x, y) = ht

c2 = c2(|u|)

(2.5)


dove l’ultima rappresenta sostanzialmente l’equazione di stato che modella il comporta-mento del gas. Per un gas ideale politropico essa e del tipo:

c2 = (γ − 1)(ht − |u|22

)

e si ricava utilizzando la (1.2) e la (1.3). La seconda equazione impone l’irrotazionalitadel campo di velocita e poiche il dominio considerato e semplicemente connesso, questopermette di introdurre una funzione potenziale di velocita:

φ(x, y) : u = ∇φ. (2.6)

Decidendo di trattare esclusivamente il caso bidimensionale cartesiano, e possibile ri-scrivere il sistema (2.5) come:

(φ2x − c2)φxx + 2φxφyφxy + (φ2

y − c2)φyy = 0,

c2 = c2(|u|).(2.7)

attraverso l’uso del potenziale di velocita e tenendo gia conto della omoentropicita edella omoentalpicita. Un ulteriore semplificazione e resa possibile mediante l’utilizzodel potenziale di perturbazione ϕ(x, y). Introducendo le velocita di perturbazione u e v,definite come: {

U = (1 + u)U∞V = vU∞

(2.8)

sotto le ipotesi

• u ¿ 1,

• v ¿ 1,

e possibile definire ϕ nel seguente modo:

ϕ :

{u = ϕx

v = ϕy

(2.9)

Il potenziale completo φ puo essere riscritto in funzione del potenziale di perturbazionecome:

φ = U∞x + ϕ. (2.10)

Utilizzando il gas ideale e la sua espressione della celerita del suono, trascurando itermini quadratici sotto l’ipotesi di piccole perturbazioni e considerando numeri di Machdi riferimento M∞ non troppo grandi (condizioni non ipersoniche), otteniamo un’unicaequazione non lineare:

(1−M2

∞ + (γ + 1)M2∞ϕx

)ϕxx + 2M2

∞ϕyϕxy +((γ + 1)M2

∞ϕx − 1)ϕyy = 0 . (2.11)

Un’analisi degli ordini di grandezza dei termini ci permette di considerare quelliimportanti e di trascurare gli altri, cosı da ottenere un’espressione semplice per tutti etre i regimi di moto:


Subsonico e supersonico:(M2

∞ − 1)ϕxx − ϕyy

Transonico:(M2

∞ − 1)ϕxx − ϕyy + M2∞(γ + 1)ϕxϕxx = 0, (2.12)

Se siamo interessati a cio che succede in prossimita della gola, dobbiamo considerarel’ultimo caso, dove la (2.12) e ancora un’equazione non lineare alle derivare parziali. Perapplicare le piccole perturbazioni nel caso della gola di un ugello, si prende come statodi riferimento la condizione sonica (M∞ = 1, U∞ = c∗), da cui si ottiene un equazioneancora non lineare per il potenziale comprimibile:

(γ + 1)ϕxϕxx − ϕyy = 0. (2.13)

Figura 2.1: Geometria dell’ugello e sistema di riferimento

Il metodo introdotto da Sauer permette di trovare la soluzione di questa equazioneipotizzando a priori la forma assunta dal potenziale, infatti si richiede che sia del tipo:

φ(x, y) =∞∑

i=0

f2i(x)y2i = f0(x) + f2(x)y2 + f4(x)y4 + . . . (2.14)

dove e evidente che l’utilizzo di termini esclusivamente pari e dovuto a ragioni di sim-metria rispetto all’asse. Troncando quindi la serie ai termini moltiplicati per y4 esostituendo questo sviluppo nell’equazione (2.13) si ottiene la seguente espressione:

{(γ + 1)f ′0f′′0 − 2f2}+ {(γ + 1)

[f ′0f

′′2 + f ′2f

′′0

]− 12f4}y2

+ {(γ + 1)[f ′0f

′′4 + f ′4f

′′0 + f ′2f

′′2

]− 30f6}y4 = 0.(2.15)


Poiche l’equazione (2.15) deve essere soddisfatta per qualsiasi y, il coefficiente di ciascunapotenza di y deve essere nullo:

f2 =γ + 1

2f ′0f

′′0 ,

f4 =γ + 1

2[f ′0f

′′2 + f ′′0 f ′2

].

(2.16)

In tal modo, conoscendo f ′0(x) e possibile, risolvendo in cascata il sistema, ottenere tuttele funzioni f2i(x) e quindi determinare il potenziale di perturbazione. Sauer pone

f ′0(x) = αx. (2.17)

Sostituendo le funzioni cosı trovate in (2.15) si ottengono le espressioni delle velocita:

u(x, y) = ϕx = αx +γ + 1

2α2y2, (2.18)

v(x, y) = ϕy =γ + 1

2α2xy +

(γ + 1)2

6α3y3. (2.19)

A questo punto e possibile determinare la linea sonica (dove M = 1), notando cheessendo questa la condizione attorno a cui e stata linearizzata l’equazione del potenziale,e sufficiente sostituire la condizione di perturbazione nulla (u = 0) in (2.18):

x = −γ + 12

αy2. (2.20)

Per determinare l’origine del sistema di riferimento si impone la condizione di velocitanormale nulla a parete nel punto in corrispondenza della gola:

ε = −γ + 16

αy2t , (2.21)

dove il parametro α puo essere determinato mediante considerazioni geometriche e risultapari a:

α =√

1γ + 1

1rtyt

. (2.22)

La soluzione delle ultime due equazioni algebriche permette di determinare l’intero cam-po di moto bidimensionale in un intorno della gola e, come si puo ben osservare, essedipendono solo da rt e yt; possiamo qunidi affermare che, per il metodo di Sauer, l’interageometria di un ugello e esclusivamente rappresentata dal raggio di curvatura in golae dalla altezza della gola. E interessante capire, a questo punto, quali sono i valori diquesti due parametri sufficienti a garantire il verificarsi delle ipotesi che sono alla basedel concetto di piccole perturbazioni: ovvero quali ugelli possiamo sperare di descriveree quanto lontano dalla condizione di riferimento e possibile spingersi. Generalmente ilmetodo di Sauer viene considerato corretto solo per valori di rt/yt > 2, il che limita nonpoco le sue possibilita di utilizzo.


2.2.1 Caso assisimmetrico

Introducendo un sistema di coordinate cilindriche uz, ur, uθ e definendo le ipotesi ulte-riori di:

• assialsimmetria: ∂∂θ = 0,

• swirl nullo: uθ = 0,

possiamo riscrivere la prima equazione del sistema (1.5) nel seguente modo:

(u2r − c2(|u|))∂ur

∂r+ 2uruz

∂uz

∂r+ (u2

z − c2(|u|))∂uz

∂z− c2(|u|)ur

r= 0. (2.23)

Sotto tali ipotesi il problema diventa bidimensionale, percio, per comodita di scrittura,e possibile introdurre il seguente cambio di variabili: z = x e r = y.Operando questa sostituzione e rifacendosi all’espressione (2.6) del potenziale di velocita,e possibile riscrivere l’equazione (2.23) come:

(φ2x − c2)φxx + 2φxφyφxy + (φ2

y − c2)φyy − c2 φy

y= 0. (2.24)

In modo analogo a quanto gia descritto per il caso cartesiano, utilizzando il potenziale diperturbazione e linearizzando attorno alla condizione di riferimento M∞ = 1, la (2.24)diventa:

(γ + 1)ϕxϕxx − ϕyy − ϕy

y= 0. (2.25)

Sostituendo la forma del potenziale ipotizzata da Sauer ed imponendo l’annullamentodei coefficienti si ottiene:

f2 =γ + 1

4f ′0f

′′0 ,

f4 =γ + 116

[f ′0f

′′2 + f ′′0 f ′2

],

f6 =γ + 136

[f ′0f

′′4 + f ′4f

′′0 + f ′2f

′′2

],

(2.26)

che, con la scelta di f ′0(x) = αx operata da Sauer, diventano:

f2 =γ + 1

4α2x,

f4 =(γ + 1)2

64α3,

f6 = 0.

(2.27)

Le espressioni di f2 e f4 sostituite nell’equazione (2.15) permettono di trovare le velocitau e v :

u(x, y) = ϕx = αx +γ + 1

4α2y2, (2.28)

v(x, y) = ϕy =γ + 1

2α2xy +

(γ + 1)2

16α3y3. (2.29)


Analogamente a quanto gia fatto si ottiene poi:

x = −(γ + 1)4

αy2, ε = −(γ + 1)8

αy2t . (2.30)

2.3 Il metodo di Hall-Kliegel assisimmetrico

L’analisi transonica di Sauert e la piu semplice tra quelle proposte per risolvere il pro-blema del flusso in un ugello convergente-divergente. Come detto in precedenza, talemetodo e adeguato solo per valori di rt/yt > 2. In questa sezione sono velocementeanalizzati alcuni metodi che hanno una accuratezza superiore e sono validi anche pervalori di rt/yt < 2.

Il metodo di Hall e basato sull’impiego di una espansione in serie per le componentidella velocita in potenze inverse del parametro di espansione R. Le seguenti forme diespansione in serie sono assunte, con x misurato a partire da xt:

u =u1(r, z)

R+

u2(r, z)R2

+u3(r, z)

R3+ . . . (2.31)

v =

√γ + 12R

[v1(r, z)

R+

v2(r, z)R2

+v3(r, z)

R3+ . . .

](2.32)

dove

R =rt

ytr =

y

ytz =

x

yt

√2R

γ + 1.

Le funzioni u1, v1, etc.. devono essere determinate in modo che le equazioni di governoe le condizioni al contorno siano soddisfatte. Le espansioni in serie sono sostituite all’in-terno dell’equazione non lineare completa (2.7), espressa naturalmente per la velocitae non per il potenziale, e tutti i termini relativi alla stessa potenza di R sono raccolti.La prima approssimazione e ottenuta considerando il coefficiente 1/R, la seconda 1/R2

e cosı via. La prima porta ad una soluzione identica a quella ottenuta da Sauer; sesi considerano invece termini di ordine superiore l’accuratezza aumenta e il metodo evalido per valori di R > 1. I seguenti risultati sono stati ottenuti da Hall per il casoassisimmetrico.

u1 =12r2 − 1

4+ z

v1 =14r3 − 1

4r + rz

u2 =2γ + 9

24r4 − 4γ + 15

24r2 +

10γ + 57288

+ z

(r2 − 5

8

)− 2γ − 3

6z2

v2 =γ + 3

9r5 − 20γ + 63

96r3 +

28γ + 93288

r + z

(2γ + 9

6r3 − 4γ + 15

12r

)+ rz2

2.4. Metodo delle caratteristiche 15

Il metodo dovuto a Kliegel e una modifica di quello dovuto ad Hall. Esso rimpiazzail parametro di espansione R nelle equazioni (2.31) e (2.32) con R + 1. La serie chesi ottiene e, a questo punto, convergente per valori di R < 1. La soluzione di Kliegelprende la forma

u =u1

R + 1+

u2 + u1

(R + 1)2+

u3 + 2u2 + u1

(R + 1)3+ . . . (2.33)

v =

√γ + 1

2(R + 1)

[v1

R + 1+

32v2 + v1

(R + 1)2+

v3 + 52v2 + 15

8 v1

(R + 1)3+ . . .

](2.34)

dove u1, v1, ecc... sono quelli calcolati per il metodo di Hall.

2.4 Metodo delle caratteristiche

Il metodo delle caratteristiche permette di risolvere l’equazioni di Eulero 2D (2.5), sottoipotesi del tutto generali che sono:

• η(x) ∈ C1

• grandezze che descrivono il flusso continue,

• flusso omoentalpico, omoentropico e supersonico.

Per un sistema cartesiano ortogonale, indicando con u e v le componenti di velocita, ilsistema puo essere scritto come:

(u2 − c2)ux + (2uv)vy + (v2 − c2)vy = 0uy − vx = 0c2 = c2(|u|)

(2.35)

Come gia detto in precedenza, la natura matematica delle equazioni di Eulero, varia alvariare del rapporto tra velocita del gas e celerita del suono e si puo dimostrare che incaso di corrente supersonica, il sistema di equazioni (2.35) e iperbolico.Il metodo delle caratteristiche permette di risolvere equazioni alle derivate parziali ditipo iperbolico, sotto l’ipotesi fondamentale della continuita della soluzione. Il puntocardine del metodo e la possibilita di trasformare un sistema di n equazioni alla derivateparziali in uno formato da 2n equazioni alle derivate ordinarie del tutto equivalente, lacui risoluzione risulta, in termini teorici, piu semplice.

E pertanto necessario considerare il caso completamente supersonico, e quindi re-stringere il processo risolutivo alla parte d’ugello interessata dal flusso a M(x) > 1,generalmente il divergente. Per questo motivo imponiamo:

M2(x) =u2(x) + v2(x)

c2(x)> 1


ovunque nel campo di moto.Un sistema iperbolico, come quello che stiamo considerando, e caratterizzato dalla

esistenza di linee particolari, dette caratteristiche, lungo le quali e sufficiente integrareuna equazione del prim’ordine alle derivate ordinarie, generalmente non lineare, perdeterminare la soluzione del problema PDE. In un sistema 2 × 2 come il (2.35), leequazioni che permettono di risolvere il campo di moto sono due e devono essere risoltecontemporaneamente. Per poter determinare la posizione delle linee caratteristiche, chee incognita, e necessaria l’introduzione di due ulteriori equazioni, dette di compatibilita,che sono generalmente anch’esse non lineari.

Il prezzo piu alto da pagare, per aver effettuato questo tipo di operazione, e larinuncia alla possibilita di rappresentare soluzioni discontinue e quindi urti, all’internodel campo di moto. L’unica concessione che il metodo fa in questo senso, e la possibiltadi capire quando un’urto si crea, segnato dall’accavallarsi delle caratteristiche: il cheimplica anche la fine del campo di validita del metodo stesso.

Definiamo, a questo punto, le linee caratteristiche come i luoghi del piano chesoddisfano l’equazione:

y = Y (x)

e, in maniera analoga, introduciamo le funzioni maiuscole, che sott’intendono il fattoche la generica quantita g(x, y) e valutata lungo una linea caratteristica, quindi G =g(x, Y (x)) = G(x), per cui:

U(x) = u(x, Y (x)),V (x) = v(x, Y (x)),C(x) = c(U(x), V (x)).

Introduciamo anche le derivate composte di tali quantita:

dU

dx=

∂u

∂x+

∂u

∂y

dY

dx= ux + uy

dY

dx,

dV

dx=

∂v

∂x+

∂v

∂y

dY

dx= vx + vy

dY

dx.

A questo punto, partendo dalla (2.35), tenendo conto del vincolo di irrotazionalita eutilizzando le precedenti, si ottiene la seguente equazione per la derivata parziale di urispetto a y:

uy =(U2 − C2)

dU

dx

dY

dx+ (V 2 − C2)

dV

dx

(U2 − C2)(

dY

dx

)2

+ 2UVdY

dx+ (V 2 − C2)

. (2.36)

Come si vede, si puo calcolare uy ovunque tranne che lungo le linee caratteristiche.L’annullarsi del numeratore consente di determinare l’equazione caratteristica, mentre


ponendo a zero il denominatore permette di ottenere dell’equazione di compatibilita,un’equazione quadratica che ha quindi due soluzioni del tipo:

dY

dx

±=

UV ± C√

U2 + V 2 − C2

U2 − C2= λ±(U, V ). (2.37)

Come si vede nella (2.37) la pendenza delle linee caratteristiche dipende esclusivamentedalla soluzione del sistema (2.35) e mai direttamente dalla posizione; questo fatto emolto importante, dal punto di vista metodologico, per poter risolvere numericamentele equazioni di Eulero. E anche evidente come, in ogni punto del campo di moto, esistanodue pendenze diverse, corrispondenti a due linee caratteristiche differenti dette C+ e C−,le quali identificano rispettivamente anche le due famiglie di appartenenza.

Dividendo l’equazione di compatibilita perdU

dxe poi esplicitando per questo termine,

otteniamo la seguente equazione differenziale ordinaria:

dV

dU=

U2 − C2

C2 − V 2λ(U, V ) (2.38)

che non dipende direttamente dalla posizione, ma solamente dalle componenti dellavelocita U e V . E quindi possibile pensare di integrarla direttamente, dato che e nellaforma dy

dx = f(x, y). Facendo tale operazione di divisione, e in realta stato operato uncambio di variabile, scrivendo V in funzione di U , definendo di fatto una nuova funzioneV = V (U). Per cui la (2.38), valutata per le due diverse soluzioni di λ, diventa:

dV

dU

±=

[U2 − C2

C2 − V 2λ(U, V )

]±

ottenendo quindi due versioni dell’equazione precedente, a seconda di quale caratteristicasi osserva.

Il sistema di due equazioni (2.35), corrisponde a un sistema di quattro equazioni, acui va aggiunta per entrambi l’equazione di stato del gas.

λ±(U, V ) =UV ± C

√U2 + V 2 − C2

U2 − C2

dV

dU

±=

[U2 − C2

C2 − V 2λ(U, V )

]±

C2 = c2(U, V )

(2.39)

A tale sistema, naturalmente, vanno aggiunte le condizioni al contorno necessarie arendere ben posto il problema di Cauchy. Dal punto di vista del metodo delle carat-teristiche, anche se non risulta molto chiaro dalla (2.39), e necessario specificare undato iniziale che permetta l’avvio del metodo. Poiche le caratteristiche rappresentanole linee spaziali lungo le quali si propaga l’informazione, si intuisce che aver definito ildato iniziale lungo una di tale linee costituisce un problema mal posto. Alla linea del


dato iniziale devono essere affiancate le altre condizioni al contorno del problema, chevengono dette di inflow-outflow. Possono essere formulate come condizioni di non pe-netrazione, come nel caso della parete dell’ugello, oppure imponendo la pressione, comesuccede all’espulsione in atmosfera dei gas.

»»»9

6

¾x

y

z

6»»9

Rθ

θ

Figura 2.2: Geometria assisimmetrica.

2.4.1 Caso assisimmetrico

Per lo studio della corrente che si sviluppa all’interno di un ugello con geometria assi-simmetrica, e necessario riformulare il sistema (2.5) in coordinate cilindriche (r, z, θ). Ilvettore velocita espresso in questo nuovo sistema di coordinate diventa quindi:

u(R, z, θ) = uR(R, z, θ) R(θ) + uz(R, z, θ) z + uθ(R, z, θ) θ(θ). (2.40)

La divergenza in coordinate cilindriche assume la seguente forma:

∇·u =1R

∂

∂R(RuR) +

∂uz

∂z+

1R

∂uθ

∂θ, (2.41)

mentre il termine convettivo si scrive come:

u ·∇ = uR∂

∂R+

uθ

R

∂

∂θ+ uz

∂

∂z. (2.42)

Infine l’espressione dell’operatore rotore che compare nell’equazione di continuita e:

∇×u =[

1R

∂uz

∂θ− ∂uθ

∂z,

∂uR

∂z− ∂uz

∂R,

1R

∂

∂R(Ruθ)− 1

R

∂uR

∂θ

]T

. (2.43)

Sotto le uteriori ipotesi di:

• simmetria assiale: ∂∂θ = 0,

• nessuna componente di swirl: uθ = 0,


il sistema (2.5) diventa un sistema di 2 equazioni in 2 incognite come nel caso cartesiano;ponendo z = x e R = y e chiamando ancora u e v le 2 componenti della velocita, possiamoriscriverlo come:

(u2 − c2)ux + (2uv)vy + (v2 − c2)vy − c2 v

y= 0

uy − vx = 0(2.44)

Tale sistema e formalmente identico al sistema (2.35), se si eccettua l’aggiunta del termi-ne c2 v

y , potenzialmente singolare per y = 0. E quindi possibile procedere in modo analo-go a quanto gia visto nel caso cartesiano, ottenendo le seguenti equazioni caratteristicheed equazioni di compatibilita:

[dY

dx

]±= λ±(U, V ) =

UV ± C√

U2 + V 2 − C2

U2 − C2, (2.45)

(U2 − C2)λ±(U, V )dU

dx+ (V 2 − U2)

dV

dx− C2λ±(U, V )

V

Y= 0. (2.46)

Mentre l’equazione carratteristica (2.45) e identica alla (2.37), l’equazione di compati-bilita (2.46) differisce dalla (2.38) per la comparsa del termine singolare C2λ(U, V )V

Y .Rimane infine il problema delle condizioni sull’asse. Al fine di rispettare il vincolo

di assisimmetria e necessario che il flusso di tutte le quantita scalari attraverso l’asse sianullo: ∮

Γasse

f(w) · ndγ = 0, (2.47)

dove, nel caso delle equazioni di Eulero complete:

f(w) ≡

ρu ρvρu2 + P ρuv

ρuv ρv2 + Pu(Et + P ) v(Et + P )

.

Per garantire l’annullamento dell’integrale nella condizione (2.47), e sufficiente imporrel’annullamento sull’asse della funzione integranda:

f(w) · n|asse = 0,

che nel caso in esame equivale a:

fρ(w) · n = [ρu ρv] [0 − 1]T = 0.

Pertanto il vincolo di assisimmetria diventa:

v|asse = 0. (2.48)

Capitolo 3

Algoritmica

In questo capitolo verranno descritti gli algoritmi di soluzione del metodo quasi-1D edel metodo delle caratteristiche.

3.1 Metodo quasi–monodimensionale

Per studiare la corrente che fluisce in un ugello il primo passo e fissare alcuni parametrifondamentali, quali:

- γ: natura del gas,

- P0, T0: condizioni di serbatoio,

- Ae/Ag: rapporto tra area di efflusso e area di gola,

- Pa/P0: rapporto tra pressione ambiente e pressione totale.

Con l’ipotesi che il rapporto Ai/Ag tra la sezione d’ingresso e quella di gola sia suffi-cientemente grande, e possibile ritenere che le condizioni sulla sezione d’ingresso sianopraticamente coincidenti con quelle del serbatoio.Denotando con Pe la pressione sulla sezione di efflusso dell’ugello, si possono presentaretre diversi casi:

1. Pe = P1 : la corrente raggiunge condizioni soniche in gola e torna ad essere sub-sonica nel divergente.Per il calcolo della pressione P1 si sfrutta inanzitutto la relazione (2.3) per correntistazionarie isoentropiche che lega il Mach e la geometria dell’ugello:

Ae

Ag=

1Me

[2

γ + 1

(1 +

γ − 12

M2e

)] γ+12(γ−1)

. (3.1)

Questa relazione fornisce due valori di Mach all’efflusso distinti, uno subsonico euno supersonico: ovviamente in questo caso sara di interesse soltanto la soluzione

21

22 Capitolo 3. Algoritmica

subsonica, denominata Me,1. Noto il Mach all’efflusso, e possibile calcolare P1

sfruttando direttamente la relazione isoentropica (2.4) per le pressioni:

P1 =P0(

1 + γ−12 M2

e,1

) γγ−1

. (3.2)

2. Pe = P2 : la corrente raggiunge condizioni soniche in gola e diventa supersonicanel divergente.Il procedimento per il calcolo della pressione P2 e lo stesso di quello adottato nelcaso precedente, con l’accorgimento di sfruttare la soluzione supersonica (Me,2)che risulta dalla (3.1).

3. Pe = P3 : formazione di un urto in corrispondenza della sezione di efflusso.In questo caso e necessario ricorrere alle relazioni d’urto: ricavate le grandezzesupersoniche a monte dell’urto (indicate col pedice 2) come illustrato al puntoprecedente, e possibile trovare la pressione a valle (indicata col pedice 3):

P3 = P2

2γM2e,2 − (γ − 1)(γ + 1)

. (3.3)

Note P1, P2, P3, e possibile studiare la corrente nell’ugello al variare della pressioneambiente Pa.

3.1.1 Soluzione subsonica: Pa ≥ P1

In questo caso Pe = Pa e il flusso nell’ugello e completamente isoentropico subsonico:con tali condizioni al contorno non si raggiungono condizioni soniche in gola e dunquel’area critica risulta essere minore dell’area di gola.Il procedimento da seguire per la risoluzione della corrente nell’ugello e il seguente:

1. calcolo del Mach sulla sezione di efflusso attraverso la relazione isoentropica per lepressioni:

P0

Pa=

(1 +

γ − 12

M2e

) γγ−1

⇒ Me =

√√√√ 2γ − 1

[(P0

Pa

) γ−1γ

− 1

]; (3.4)

2. determinazione dell’area critica A∗:

A∗ =Ae

1Me

[2

γ+1

(1 + γ−1

2 M2e

)] γ+12(γ−1)

; (3.5)

3. calcolo dell’andamento del Mach nel condotto attraverso l’espressione implicita:

A(M) = A∗1M

[2

γ + 1

(1 +

γ − 12

M2

)] γ+12(γ−1)

; (3.6)

3.1. Metodo quasi–monodimensionale 23

4. determinazione dell’andamento della pressione sempre attraverso la relazione isoen-tropica:

P

P0=

(1 +

γ − 12

M2

)− γγ−1

. (3.7)

3.1.2 Ugello adattato: Pa = P2

In questo caso si ha ancora Pe = Pa e il flusso e ovunque isoentropico, ma nel divergentediventa supersonico. L’area critica coincide con la sezione di gola (A∗ = Ag) e quindi larisoluzione della corrente risulta particolarmente facile:

1. calcolo dell’andamento del Mach nel condotto attraverso l’espressione implicita:

A(M) = Ag1M

[2

γ + 1

(1 +

γ − 12

M2

)] γ+12(γ−1)

; (3.8)

2. determinazione dell’andamento della pressione attraverso la relazione isoentropica:

P

P0=

(1 +

γ − 12

M2

)− γγ−1

. (3.9)

3.1.3 Urto normale nel divergente: P3 ≤ Pa < P1

La corrente e ovunque isoentropica a meno della sezione in corrispondenza della qualesi forma l’onda d’urto: l’intensita dell’urto dipende dalla sua posizione, e la posizionedipende a sua volta dal rapporto Pa/P0. Essendo un fenomeno adiabatico, attraversol’urto si conserva la temperatura totale T0 ma non la pressione totale P0.I passaggi per risolvere questo tipo di flusso sono:

1. calcolo del Mach sulla sezione di efflusso:

Me =

√−1 +

√1 + β2(γ − 1)γ − 1

(3.10)

dove

β =P0

Pa

Ag

Ae

(2

γ + 1

) γ+12(γ−1)

; (3.11)

2. determinazione della pressione totale a valle dell’urto:

P0,valle = Pa

(1 +

γ − 12

M2e

) γγ−1

; (3.12)

3. calcolo del Mach a monte dell’urto per mezzo dell’espressione implicita:

P0,valle

P0=

[1 +

2γ

γ + 1(M2

monte − 1)] −1

γ−1[

(γ + 1)M2monte

(γ − 1)M2monte + 2

] γγ−1

; (3.13)


4. risoluzione della corrente isoentropica a monte dove M ∈ (0,Mmonte);

5. risoluzione della corrente isoentropica a valle dove M ∈ (Mvalle,Me);

6. determinazione della posizione dell’urto attraverso il calcolo della sezione a cuicorrisponde M = Mmonte:

Amonte =Ag

1Mmonte

[2

γ+1

(1 + γ−1

2 M2monte

)] γ+12(γ−1)

. (3.14)

3.1.4 Ugello sovraespanso: P2 < Pa < P3

Poiche Pe < Pa, avviene una ricompressione all’esterno dell’ugello tramite onde d’urtooblique.

3.1.5 Ugello sottoespanso: Pa < P2

In questo caso Pe > Pa e si forma un ventaglio di espansione all’esterno dell’ugello.

3.2 Metodo delle caratteristiche

Γ

•

•

•P+

P−

I

Figura 3.1: Unit Process

L’algoritmo su cui si basa la procedura di soluzione mediante il metodo delle carat-teristiche prende il nome di Unit process. Con riferimento alla Figura 3.1, tale algoritmopermette di ottenere la soluzione del sistema (2.39) nel punto I di intersezione dellecaratteristiche C+ e C−, che hanno origine rispettivamente nei punti P+ e P−, che sitrovano lungo la linea Γ dove la soluzione e nota:

V +(UP+) = VP+ ,

V −(UP−) = VP− .


A questo punto e possibile integrare separatamente le 2 equazioni di compatibilita,ottenendo:

∫

V +

P+

dV + = −∫

U+

P+

(U+)2 − (C+)2

(V +)2 − (C+)2λ+(U+, V +) dU+,

∫

V −P−

dV − = −∫

U−P−

(U−)2 − (C−)2

(V −)2 − (C−)2λ−(U−, V −) dU−,

(3.15)

dove gli estremi di integrazione sono incogniti. Tuttavia, essendo il punto I il punto diintersezione delle caratteristiche, in tal punto le soluzioni coincidono:

V +I = V −

I ,

U+I = U−

I .

E quindi possibile integrare le equazioni (3.15) tra P± e I:

∫ V +I

V +

P+

dV + = −∫ U+

I

U+

P+

(U+)2 − (C+)2

(V +)2 − (C+)2λ+(U+, V +) dU+,

∫ V −I

V −P−

dV − = −∫ U−I

U−P−

(U−)2 − (C−)2

(V −)2 − (C−)2λ−(U−, V −) dU−,

(3.16)

ottenendo il seguente sistema di 4 equazioni nelle 4 incognite U+I , U−

I , V +I e V −

I :

dV

dU

+

= β+(U+, V +)

dV

dU

−= β−(U−, V −)

V +I = V −

I

U+I = U−

I

V +(U+P+) = v+

P+

V −(U−P−) = v−

P−

dove le ultime due rappresentano le condizioni al contorno, poiche vP+ e vP− sono note,e β± = (U±)2−(C±)2

(V ±)2−(C±)2λ±(U±, V ±). E’ possibile poi ottenere un sistema di 2 equazioni


nelle incognite UI e VI , considerando il fatto che U+I = U−

I = UI :

VI =∫ UI

UP+

β+(U+, V +) dU+

VI =∫ UI

UP−β−(U−, V −) dU−

V +(U+P+) = v+

P+

V −(U−P−) = v−

P−

(3.17)

La soluzione di questo sistema, fortemente non lineare, non richiede la conoscenza delladislocazione spaziale del punto di intersezione. Cio significa che e sempre possibilerisolverlo indipendentemente dalle equazioni caratteristiche. Una volta risolto questosistema, per determinare la posizione del punto di intersezione delle caratteristiche enecessario risolvere l’equazione caratteristica

dY

dx

+

= λ+(U+(x), V +(x)

)

dY

dx

−= λ−

(U−(x), V −(x)

)

Y +(xI) = yI = Y −(xI)

Y +(xP+) = yP+

Y −(xP−) = yP−

(3.18)

dove le ultime due equazioni rappresentano le condizioni al contorno, mentre la terzae un vincolo, che costituisce una sorta di terza condizione al contorno. Come si puoben vedere, questo sistema e decisamente problematico, in quanto e possibile conoscerei valori di U± e V ± solo agli estremi, cosı come mostrato in Figura 3.1 e dall’ equazione(3.17). Questo rende necessario l’utilizzo di metodi di integrazione numerica che nonrichiedano la valutazione della funzione λ±

(U±(x), V ±(x)

)nei punti interni, ma solo

agli estremi, che risultano quindi poco accurati. Un esempio di tali metodi e quello deitrapezi, che opera una linearizzazione delle equazioni:

YI − YP+ =λ+

(UP+ , VP+

)+ λ+

(UI , VI

)

2(xI − xP+)

YI − YP− =λ−

(UP− , VP−

)+ λ−

(UI , VI

)

2(xI − xP−)

che rappresenta un sistema di 2 equazioni nelle 2 incognite YI e xI .


Figura 3.2: Rete caratteristica

Applicando questa procedura a partire da un certo numero di punti iniziali giacientisu una linea che non sia una caratteristica, in cui si suppone noto il vettore velocita, siottiene la rete caretteristica riportata in figura 3.2. Se n e il numero di punti inizialiall’i-esimo passo, lo Unit Process permette di calcolare la soluzione solamente in n − ipunti. La procedura sembra quindi destinata ad esaurirsi dopo n− 1 iterazioni.

•PC+

y(x)•I

Figura 3.3: Unit Process al contorno.

Affinche ad ogni step il numero di punti in cui calcoliamo la soluzione non si ridu-ca, e necessario imporre delle condizioni al contorno. In corrispondenza del contornosolido, la direzione del vettore velocita deve essere parallela alla pendenza del contorno.Come riportato in figura 3.3, solamente una caratteristica interseca il contorno soli-do, quindi sono disponibili solamente una equazione caratteristica ed una equazione dicompatibilita. Tuttavia, sul contorno abbiamo a disposizione 2 ulteriori equazioni:

y = η(x) (3.19)dη

dx=

v

u. (3.20)


dove l’equazione (3.19) descrive la geometria dell’ugello, mentre la (3.20) rappresenta lacondizione di non compenetrabilita sul contorno solido. Otteniamo cosı, in riferimentoal caso mostrato in figura 3.3, il seguente sistema di equazioni:

YI = η(xI)

UIdη

dx= VI .

YI − YP+ =λ+

(UP+ , VP+

)+ λ+

(UI , VI

)

2(xI − xP+)

VI =∫ UI

UP+

β+(U+, V +) dU+

(3.21)

La formulazione del metodo delle caratteristiche presentata e la piu generale possibileed e facilmente estendibile al caso assisimmetrico. Infatti il caso assisimmetrico risulta,come gia visto, formalmente identico al caso cartesiano. E pertanto conveniente, dalpunto di vista implementativo, adottare la seguente formulazione:

dY

dx= λ(U, V )

(U2 − C2)λ(U, V )dU

dx+ (V 2 − C2)

dV

dx− δ C2λ(U, V )

V

Y= 0

(3.22)

dove {δ = 0 nel caso cartesianoδ = 1 nel caso assisimmetrico

Inoltre i sistemi non lineari ottenuti mediante tale formulazione possono essere facil-mente risolti mediante un metodo predictor-corrector. In particolare si e utilizzato unmetodo predictor-corretor di Eulero modificato, che ha un’accuratezza del secondo ordinenell’integrazione di equazioni differenziali ordinarie. L’algoritmo di correzione si basasul metodo dei coefficienti mediati: i valori dei coefficienti dell’equazione differenzialesono calcolati come media tra i valori nel punto iniziale e nel punto di intersezione dellecaratteristiche. Data l’equazione differenziale ordinaria:

dy

dx= f(x, y),

il metodo consiste in un primo passo di predizione, in cui viene stimato in manieragrossolana il valore dell’incognita:

P: y0I = yP + f(xP , yP )∆x,

seguito da passi di correzione in cui la stima cosı effettuata viene progressivamenteraffinata:

C: ymI = yP + f

(xP +

∆x

2,yP + ym−1

I

2

)∆x.


Tale schema applicato al caso in esame assume la seguente forma:

P:

[U2

P+ − C2P+

]λ+

P+

U0I − UP+

x0I − xP+

+[V 2

P+ − C2P+

] V 0I − VP+

x0I − xP+

− δ C2P+λ+

P+

VP+

yP+

= 0

[U2

P− − C2P−

]λ−

P−U0

I − UP−

x0I − xP−

+[V 2

P− − C2P−

] V 0I − VP−

x0I − xP−

− δ C2P−λ−

P−VP−

yP−= 0

y0I − yP+

x0I − xP+

= λ+(UP+ , VP+)

y0I − yP−

x0I − xP−

= λ−(UP− , VP−)

C:

[U2⊕ − C2

⊕]λ+⊕

UmI − UP+

xmI − xP+

+[V 2⊕ − C2

⊕] V m

I − VP+

xmI − xP+

− δ C2⊕λ+

⊕V⊕y⊕

= 0

[U2ª − C2

ª]λ−ª

UmI − UP−

xmI − xP−

+[V 2ª − C2

ª] V m

I − VP−

xmI − xP−

− δ C2ªλ−ª

Vªyª

= 0

ymI − yP+

xmI − xP+

= λ+(U⊕, V⊕)

ymI − yP−

xmI − xP−

= λ−(Uª, Vª)

dove

U⊕ =U

(m−1)I + UP+

2Uª =

U(m−1)I + UP−

2

V⊕ =V

(m−1)I + VP+

2Vª =

V(m−1)I + VP−

2

y⊕ =y

(m−1)I + yP+

2yª =

y(m−1)I + yP−

2

C⊕ = c(U⊕, V⊕) Cª = c(Uª, Vª)

λ+⊕ = λ+(U⊕, V⊕) =

U⊕V⊕ + C⊕√

U2⊕ + V 2

⊕ − C2⊕

U2⊕ − C2

⊕

λ−ª = λ−(Uª, Vª) =UªVª − Cª

√U2ª + V 2

ª − C2ª

U2ª − C2

ª

Il passo di corrector viene ripetuto fino ad ottenere l’accuratezza desiderata, utilizzandoper l’arresto un criterio sul residuo:

∣∣∣∣∣( · )m−1

I − ( · )mI

( · )mI

∣∣∣∣∣ < ε.


Come si puo vedere, i problemi legati al termine singolare sono stati risolti scegliendo ditrascurare il termine singolare nel passo di predictor. In questo modo, se la discretizzazio-ne dei punti sulla linea iniziale e tale da evitare di collocare un nodo in corrispondenzadell’asse, il termine non singolare non viene mai valutato in corrispondenza dell’assestesso.

Per chiarezza riportiamo il sistema risolutivo per il passo di corrector:

−λ+(U⊕, V⊕) 1 0 0

−λ−(Uª, Vª) 1 0 0

−δ C2⊕λ+

⊕ 0(U2⊕ − C2

⊕)λ+⊕

(V 2⊕ − C2

⊕)

−δ C2ªλ−ª 0

(U2ª − C2

ª)λ−ª

(V 2ª − C2

ª)

xmI

ymI

UmI

V mI

=

=

yP+ − λ+(U⊕, V⊕)xP+

yP− − λ−(Uª, Vª)xP−

−δ C2⊕λ+

⊕V⊕y⊕

xP+ +(U2⊕ − C2

⊕)λ+⊕UP+ +

(V 2⊕ − C2

⊕)VP+

−δ C2ªλ−ª

Vªyª

xP− +(U2ª − C2

ª)λ−ªUP− +

(V 2ª − C2

ª)VP−

Come si puo vedere, in tale sistema le prime 2 equazioni sono indipendenti dalle ultime.E pertanto possibile pensare dividere questo sistema 4×4 in due sistemi 2×2 da risolverein cascata.Rimane a questo punto da considerare il problema dei contorni. Trattando il problemadel flusso all’interno di ugelli che presentano delle simmetrie, vi sono due tipologie dicontorni da considerare: contorni solidi e assi di simmetria. Come gia visto trattandoil vincolo sull’asse per le caratteristiche nel caso assisimmetrico, queste contorni sonoformalmente diversi. Tuttavia, dal punto di vista implementativo risulta convenienteconsiderarli entrambi come contorni solidi, su cui la pendenza del vettore velocita eimposta. Trattando quindi i punti estremi, in cui si ha la condizione vista in Figura 3.3,il metodo di predictor-corrector assume la seguente forma:

P:

[U2

P+ − C2P+

]λ+

P+

U0I − UP+

x0I − xP+

+[V 2

P+ − C2P+

] V 0I − VP+

x0I − xP+

− δ C2P+λ+

P+

VP+

yP+

= 0

y0I − yP+

x0I − xP+

= λ+(UP+ , VP+)

y0I = η(x0

I)

U0I

dη

dx= V 0

I


C:

[U2⊕ − C2

⊕]λ+⊕

UmI − UP+

xmI − xP+

+[V 2⊕ − C2

⊕] V m

I − VP+

xmI − xP+

− δ C2⊕λ+

⊕V⊕y⊕

= 0

ymI − yP+

xmI − xP+

= λ+(U⊕, V⊕)

ymI = η(xm

I )

UmI

dη

dx= V m

I

Poiche la funzione che descrive l’ugello e nota, e possibile sostituirla direttamente nellaequazione di compatibilia, ottenendo cosı una equazione non lineare per la sola incognitaxm

I :η(xm

I )− λ+(U⊕, V⊕)xmI = yP+ − λ+(U⊕, V⊕)xP+

che una volta risolta mediante opportuni algoritmi, come bisezione o Newton, permettela conoscenza anche di ym

I , e consente di valutare la condizione al contorno dηdx(xm

I ).Per inizializzare il metodo delle caratteristiche, la linea dei dati iniziali viene calcolatamediante il metodo di Sauer (oppure il metodo di Kliegel nel caso assisimmetrico),utilizzando il luogo dei punti in cui la pendenza del vettore velocita e nulla (quindi lalinea V = 0) e una discretizzazione variabile, al fine di consentire un addensamento deidati iniziali nelle zone di interesse.

Capitolo 4

Risultati

Figura 4.1: L’ugello di riferminento cartesiano.

Figura 4.2: L’ugello di riferminento assisimmetrico.

In questo capitolo verranno confrontati i risultati ottenuti calcolando il flusso al-l’interno di diversi ugelli convergenti-divergenti utilizzando i metodi precedentemente

33

34 Capitolo 4. Risultati

Figura 4.3: L’ugello di riferminento rt/yt = 1.2.

esposti. La geometria, che chiamerermo di riferimento, e disegnata utilizzando due cur-ve di Bezier ed e riportata in Figura 4.3. A partire da questa, modificando i coefficientiche descrivono la curva, sono state ottenute tutte le geometrie riportate in Tabella 4. Inparticolare, gli ugelli 1, 2, 3 e 4 sono ottenuti a partire dal caso di riferimento modifi-cando l’area di gola; mentre gli ugelli 5, 6 e 7 variando solamente la geometria interna emantenendo costante l’area di gola. Infine nell’ugello 8 e stata modificata anche l’area diefflusso per ottenere un caso estremo su cui testare la robustezza dei metodi. Gli ugellida 9 a 13 sono quelli utilizzati nel caso assisimmetrico e sono stati ottenuti, partendosempre dalla stessa geometria di riferimento, variando il raggio di gola.

Per rendere confrontabili le soluzioni calcolate nei diversi casi, vengono ovviamenteutilizzate le stesse condizioni di ristagno nel serbatoio; in particolare si utilizza unatemperatura totale di 300K, e un rapporto tra le pressioni p0/pa tale da evitare laformazione di urti normali nel tratto divergente.

4.1 Ugelli cartesiani

Il primo ugello analizzato e caratterizzato da una geometria molto regolare, come si puovedere osservando l’andamento del coefficiente H riportato in Figura 4.4, che risultamolto inferiore all’unita lungo tutto il condotto. Cio garantisce il pieno rispetto delleipotesi sulla geometria viste nel Capitolo 2 per il metodo quasi-monodimensionale erende quindi l’ugello adatto al primo confronto tra tutti e tre i metodi presi in esame.

Come visto nel capitolo 2, nel metodo di Sauer non viene mai presa in considerazionela geometria dell’ugello, ma solamente l’altezza della sezione di gola ed il suo raggio dicurvatura, ed i risultati ottenuti possono essere ritenuti accurati solamente se il para-metro di espansione R risulta maggiore di 2. Tuttavia, in questo caso e nel caso degliugelli 2 e 3, il parametro R risulta minore di 2, per cui i risultati calcolati con Sauer nonsaranno accurati. Pero, poiche il numero di Mach aumenta abbastanza lentamente, lasoluzione che si avra e abbastanza prossima a quella del quasi-1D da poter ritenere lasoluzione di Sauer valida in gola. Inoltre si puo notare fin da subito come Sauer risenta

4.1. Ugelli cartesiani 35

Ugello 1 cartesianoxi yi xt yt xe ye Ai/Ag Ae/Ag L rt/yt T0

0 0.9 1 0.7 3 0.9 1.29 1.29 1 1.7 300Ugello 2 cartesiano

xi yi xt yt xe ye Ai/Ag Ae/Ag L rt/yt T0





0 1.3 2 0.3 5 0.9 4.3 3 1 12.1 300Ugello 5-6-7 cartesiano


0 1.3 2 0.1 5 0.9 13 9 1 30.3 300Ugello 8 cartesiano


0 1.3 2 0.1 5 3 13 30 1 11.7 300Ugello 9 assisimmetrico


0 0.9 1 0.7 3 0.9 1.65 1.65 / 1.7 300Ugello 10 assisimmetrico


0 0.9 1 0.5 3 0.9 3.2 3.2 / 1.2 300Ugello 11 assisimmetrico


0 0.9 1 0.3 3 0.9 9 9 / 1.3 300Ugello 12 assisimmetrico


0 0.9 1 0.1 3 0.9 81 81 / 3 300Ugello 13 assisimmetrico


0 0.9 1 0.05 3 0.9 324 324 / 5.6 300

Tabella 4.1: Tabella degli ugelli. Lunghezze dimensionali espresse in [m].


pesantemente delle ipotesi semplificative che ne stanno alla base: essendo le equazioniottenute mediante una linearizzazione nell’intorno delle condizioni soniche, lo schemarisulta accurato solamente in prossimita della gola e diventa completamente inadeguatoallontanandosi da essa, in particolar modo nel tratto divergente. Risultando tuttaviaestremamente accurato nell’intorno delle condizioni soniche, rappresenta sicuramenteuna buona soluzione per il calcolo della linea dei dati iniziali necessaria per l’avvio delmetodo delle caratteristiche che, come visto nel capitolo 2, e applicabile solamente acorrenti supersoniche e permette quindi la soluzione del solo flusso che si sviluppa al-l’interno del divergente. Il risultato che si ottiene e pressoche identico a quello ottenutomediante il quasi-1D, come riportato nelle Figure 4.5 e 4.8. Al diminuire del parametroR, la soluzione di Sauer rimane valida in un intorno sempre piu piccolo delle condizionisoniche.

Com’e naturale aspettarsi, gli schemi quasi-1D e delle caratteristiche forniscono an-cora gli stessi risultati, essendo ancora verificate le ipotesi per entrambi i metodi (inparticolare l’ipotesi H < 1 del quasi-1D).

E invece molto interessante notare come nei casi degli ugelli 3 e 4 (Figure 4.10 e 4.13),dove il valore del coefficiente H e addirittura maggiore di 1, il metodo quasi-1D continuia rimanere valido, restituendo risultati in perfetto accordo con quelli del metodo dellecaratteristiche e scostandovisi di poco solamente in prossimita dell’efflusso. In Figura

0 0.5 1 1.5 2 2.5 30

0.2

0.4

0.6

0.8

1

x/xg

| A−

1/2 d

A/d

x |

Figura 4.4:∣∣∣ 1√

AdAdx

∣∣∣ lungo tutto il condotto per l’ugello cartesiano 1.

4.5 e riportato l’andamendo del numero di Mach medio sulla sezione lungo il condottoottenuto nei diversi casi.


0 0.5 1 1.5 2 2.5 3

0.6

0.8

1

1.2

1.4

1.6

x/xg

M

Quasi 1DSauerCarat.

Figura 4.5: Andamento del numero di Mach lungo il condotto per l’ugello cartesiano 1.

0 0.5 1 1.5 2 2.5 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x/xg

P/P

0

Quasi 1DSauerCarat.

Figura 4.6: Andamento del rapporto P/P0 lungo il condotto per l’ugello cartesiano 1.


0 0.5 1 1.5 2 2.5 30

0.2

0.4

0.6

0.8

1

x/xg

| A−

1/2 d

A/d

x |

Figura 4.7:∣∣∣ 1√

AdAdx


0 0.5 1 1.5 2 2.5 30.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

x/xg

M

Quasi 1DSauerCarat.

Figura 4.8: Andamento del numero di Mach lungo il condotto per l’ugello cartesiano 2.


0 0.5 1 1.5 2 2.5 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x/xg

P/P

0

Quasi 1DSauerCarat.


0 0.5 1 1.5 2 2.5 30

0.2

0.4

0.6

0.8

1

1.2

x/xg

| A−

1/2 d

A/d

x |

Figura 4.10:∣∣∣ 1√

AdAdx



0 0.5 1 1.5 2 2.5 3

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.2

x/xg

M

Quasi 1DSauerCarat.

Figura 4.11: Andamento del numero di Mach lungo il condotto per l’ugello cartesiano3.

0 0.5 1 1.5 2 2.5 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x/xg

P/P

0

Quasi 1DSauerCarat.



0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.5

1

1.5

x/xg

| A−

1/2 d

A/d

x |

Figura 4.13:∣∣∣ 1√

AdAdx


0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.5

1

1.5

2

2.5

x/xg

M

Quasi 1DSauerCarat.



0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x/xg

P/P

0

Quasi 1DSauerCarat.



Il metodo comincia a diventare meno accurato solamente a partire dall’ugello 5, per ilquale il coefficiente di riferimento diventa maggiore di 1.5 nel tratto convergente e moltoprossimo a 1 nel tratto divergente (Figura 4.16). Nonostante il mancato rispetto delleipotesi di validita, la soluzione si mantiene comunque valida, anche se poco accurata,finche il rapporto Ae/Ag si mantiene piccolo. E importante osservare, pero, che i casistudiati in questo lavoro sono solamente una parte dei casi che tale metodo e in gradodi trattare e non bisogna dimenticare che il soddisfacimento delle ipotesi di valita puodiventare un requisito molto piu stringente nello studio di casi piu complessi.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.5

1

1.5

2

x/xg

| A−

1/2 d

A/d

x |

Figura 4.16:∣∣∣ 1√

AdAdx


0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.5

1

1.5

2

2.5

3

3.5

4

x/xg

M

Quasi 1DSauerCarat.



0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x/xg

P/P

0

Quasi 1DSauerCarat.



0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.5

1

1.5

2

x/xg

| A−

1/2 d

A/d

x |

Figura 4.19:∣∣∣ 1√

AdAdx


Gli ugelli 6 e 7 presentano delle forti variazioni della geometria dell’ugello sia neltratto convergente sia in quello divergente, con valori di H maggiori di 1 sia a montesia a valle della gola. Tali irregolarita nella geometria rendono inaffidabile il metodoquasi-1D, venendo a cadere completamente una delle ipotesi in base alla quale era statoottenuto. La forti pendenze delle pareti nel tratto divergente mettono in crisi anche ilmetodo delle caratteristiche. Infatti, l’irregolarita del divergente rappresentata dai picchinei grafici 4.19 e 4.22 fa sı che le caratteristiche finiscano per intersecarsi in prossimitadella sezione di efflusso, generando quello che sarebbe un urto obliquo. Viene quindi acadere una delle ipotesi fondamentali, cioe la continuita delle grandezze che descrivonoil flusso. Come si puo vedere nelle Figure 4.21 e 4.23, l’andamento del numero di Machcalcolato con le caratteristiche diventa quindi oscillante ed il metodo non e piu valido.Da notare infine come anche per questi casi, pur limitatamente alla zone di gola, ilmetodo di Sauer continui a funzionare fornendo una soluzione accurata, essendo perquesti ugelli il parametro di espansione maggiore di 2.


0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.5

1

1.5

2

2.5

3

3.5

4

x/xg

M

Quasi 1DSauerCarat.


0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x/xg

P/P

0

Quasi 1DSauerCarat.



0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.5

1

1.5

2

2.5

3

x

| A−

1/2 d

A/d

x |

Figura 4.22:∣∣∣ 1√

AdAdx


0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.5

1

1.5

2

2.5

3

3.5

4

x

M

Quasi 1DSauerCarat.



0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x

P/P

0

Quasi 1DSauerCarat.


0 1 2 3 4 5−3

−2

−1

0

1

2

3

Figura 4.25: Ugello cartesiano No. 8.


I 3 metodi vengono infine testati sull’ugello 8, riportato in Figura 4.25, che rappre-senta un caso limite, caratterizzato da bursche variazioni di pendenza in zone relativead aree decisamente piccole, specialmente nella zona di gola (Figura 4.26). Come sipuo vedere in Figura 4.23, per quest’ugello Sauer e quasi-1D diventano completamen-te inadeguati, venendo meno le ipotesi sulle quali sono stati ricavati. Rimane invecemolto accurato il metodo delle caratteristiche se si esclude il piccolo salto alla fine deldivergente causato da un’onda d’urto obliqua in prossimita della sezione di efflusso.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.5

1

1.5

2

2.5

x

| A−

1/2 d

A/d

x |

Figura 4.26:∣∣∣ 1√

AdAdx



0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

x

M

Quasi 1DSauerCarat.


0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x

P/P

0

Quasi 1DSauerCarat.


4.2. Ugelli assisimmetrici 51

4.2 Ugelli assisimmetrici

Trattando il caso di ugelli a simmetria assiale, nell’analisi transonica al metodo di Sauersi affiancano alcuni metodi che hanno un’accuratezza maggiore e hanno un campo divalidita in funzione del parametro di espansione piu ampio, come visto nel Capitolo 2.Prima di presentare il confronto tra i diversi metodi nei casi di ugelli assisimmetricivengono presi in esame il metodo di Sauer e quello di Kliegel nello studio del flussotransonico in prossimita della sezione di gola. Come si puo vedere nelle Figure 4.30 e4.31, le isolinee del metodo di Kliegel descrivono un flusso piu uniforme ed e inoltrepossibile notare come la transizione del flusso da subsonico a supersonico risulti piugraduale rendendolo, come si vedra, accurato in un intorno delle condizioni soniche piuampio rispetto al metodo di Sauer. Nelle Figure 4.32 e 4.29 vengono comparate le lineesoniche e le linee a v = 0 ottenute con entrambi i metodi. Poiche la ε valutata conSauer risulta maggiore di quella valutata con Kliegel, entrambe le linee nel caso di Sauerrisultano piu avanzate; osservando Figura 4.29, si nota come per valori di yt piccoli i duemetodi abbiano lo stesso valore di ε, il che implica anche linee v = 0 abbastanza simili,ma all’aumentare della sezione di gola questi divergano.

yt [m] 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1R 3 1.71 1.33 1.2 1.2 1.33 1.71 3

Tabella 4.2: Paramentro di espansione R in funzione dell’altezza di gola.

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.91

1.05

1.1

1.15

yt

epsi

lon

Sauer

Kliegel

Figura 4.29: ε in funzione dell’altezza di gola yt.


0.6 0.7 0.8 0.9 1 1.1 1.2 1.30

0.1

0.2

0.3

0.4

0.5

0.6

0.26949

0.42513

0.58077

0.73642

0.73642

0.89206

0.89206

1.0477

1.04771.2033

1.2033

1.359

1.359

1.5146

1.5146

1.6703

1.6703

1.8259

1.8259

1.9816

2.1372

2.2928

2.4485

2.6041

2.7598

2.9154

3.0711

3.2267

3.3823

3.538

θ=0M=1

Figura 4.30: Isolinee di numero di Mach del metodo di Sauer.

0.6 0.7 0.8 0.9 1 1.1 1.2 1.30

0.1

0.2

0.3

0.4

0.5

0.6

0.510720.64231

0.64231

0.77389

0.90547

1.037

1.037

1.1686

1.3002

1.4318

1.4318

1.5634

1.6951.8265

1.9581

2.0897

θ=0M=1

Figura 4.31: Isolinee di numero di Mach del metodo di Kliegel.


0.85 0.9 0.95 1 1.05 1.1 1.150

0.1

0.2

0.3

0.4

0.5

v=0 Sauer

v=0 Kliegel

M=1 Sauer

M=1 Kliegel

Figura 4.32: Linee θ = 0 e M = 1 comparate.


Analizzando gli ugelli assisimmetrici e possibile fare considerazioni simili al casocartesiano. Nell’ugello 9, in cui il valore di H si mantiene prossimo a 1, i metodiquasi-1D e quello delle caratteristiche restituiscono lo stesso valore del numero di Machmedio sulla sezione lungo tutto il condotto (Figura 4.34). I metodi di Sauer e di Kliegelsono validi solamente per condizioni transoniche e risultano pertanto accurati solo inprossimita della gola.

Avendo i primi 3 ugelli assisimmetrici un rapporto di espansione minore e 2 e unayt maggiore di 0.1, a partire dall’ugello 10 (Figura 4.36) si assiste ad un netta riduzionedel dominio di validita del metodo di Sauer rispetto a quello di Kliegel, in particolarmodo nel tratto convergente, come riportato in Figura 4.37. Vale ancora quanto fattonotare nel caso cartesiano, cioe che, anche non rispettando a priori i vincoli geometri-ci monodimensionali, il quasi-1D continua a rimanere in accordo con il metodo dellecaratteristiche.

0 0.5 1 1.5 2 2.5 30

0.2

0.4

0.6

0.8

1

1.2

x

| A−

1/2 d

A/d

x |

Figura 4.33:∣∣∣ 1√

AdAdx

∣∣∣ lungo tutto il condotto per l’ugello assisimmetrico 9.


0 0.5 1 1.5 2 2.5 30

0.5

1

1.5

2

x

M

Quasi 1DSauerKliegelCarat.

Figura 4.34: Andamento del numero di Mach lungo il condotto per l’ugelloassisimmetrico 9.

0 0.5 1 1.5 2 2.5 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x

P/P

0


Figura 4.35: Andamento del rapporto P/P0 lungo il condotto per l’ugello assisimmetrico9.


0 0.5 1 1.5 2 2.5 30

0.5

1

1.5

2

2.5

x

| A−

1/2 d

A/d

x |

Figura 4.36:∣∣∣ 1√

AdAdx


0 0.5 1 1.5 2 2.5 30

0.5

1

1.5

2

2.5

3

x

M




0 0.5 1 1.5 2 2.5 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x

P/P

0



1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 31

1.5

2

2.5

3

3.5

x

M

Quasi 1DCar−MaxCar−Min

Figura 4.39: Confronto con i valori massimi e minimi di Mach nell’ugello No. 10.


0 0.5 1 1.5 2 2.5 30

0.5

1

1.5

2

2.5

3

3.5

x

| A−

1/2 d

A/d

x |

Figura 4.40:∣∣∣ 1√

AdAdx


0 0.5 1 1.5 2 2.5 30

0.5

1

1.5

2

2.5

3

3.5

4

x

M




0 0.5 1 1.5 2 2.5 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x

P/P

0



1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 31

1.5

2

2.5

3

3.5

4

4.5

x

M




L’aumento del parametro H provoca una perdita di accuratezza del quasi-1D, special-mente a partire dall’ugello 12, in cui il parametro di riferimento aumenta nel divergentefino a valori superiori a 2. In tali casi l’unico metodo che rimande valido nel trattodivergente resta quello delle caratteristiche. Infatti, come si puo vedere nelle figure 4.43e 4.47, un forte aumento dell’irregolarita della geometria comporta una forte escursionedei valori del numero di Mach lungo la sezione, particolarmente in prossimita dell’efflus-so, rendendo estremamente inaccurato il modello monodimensionale rispetto a quellobidimensionale.

A partite dall’ugello 12, il parametro di espansione R raggiunge valori maggiori di 2,ed i metodi di Sauer e Kliegel resituiscono essenzialmente gli stessi risultati nell’intornodell’aerea di gola.

Gli schemi precedenti vanno ancora piu in crisi nel caso dell’ugello 13, che rappre-senta l’ugello piu irregolare tra quelli considerati nel caso assisimmetrico e che presentauna fortissima escursione del numero di Mach lungo la sezione d’efflusso (figura 4.51).Per questo ugello il campo di funzionamento dei metodi di Sauer e Kliegel torna adessere lo stesso ed a ridursi ulteriormente rispetto ai casi precedenti. Il quasi-1D perdeulteriormente di accuratezza, e si discosta dal metodo delle caratteristiche lungo tuttoil condotto divergente.

0 1 2 30

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

x

| A−

1/2 d

A/d

x |

Figura 4.44:∣∣∣ 1√

AdAdx



0 0.5 1 1.5 2 2.5 30

1

2

3

4

5

6

7

x

M



0 0.5 1 1.5 2 2.5 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x

P/P

0




1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 31

2

3

4

5

6

7

x

M



0 1 2 30

1

2

3

4

5

x

| A−

1/2 d

A/d

x |

Figura 4.48:∣∣∣ 1√

AdAdx



0 0.5 1 1.5 2 2.5 30

1

2

3

4

5

6

7

8

9

x

M



0 0.5 1 1.5 2 2.5 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x

P/P

0




1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 31

2

3

4

5

6

7

8

9

x

M



0 0.5 1 1.5 2 2.5 30

5

10

15

Hmax

err

[%]

Assisimmetrico

Cartesiano

Figura 4.52: L’errore massimo percentuale tra caratteristiche e quasi 1D in relazione aM , in funzione del massimo valore di H.

4.3. Osservazioni 65

4.3 Osservazioni

La natura molto differente dei metodi trattati ha richiesto, come gia accennato, moltaattenzione nella scelta dei casi su cui effettuare i confronti e sulle grandezze effettiva-mente confrontabili. Per quanto riguarda il primo caso, si e dovuta porre particolarecura nella scelta delle condizioni nel serbatoio e all’efflusso, per fare evitare la forma-zione di urti normali nel tratto divergente, che avrebbero reso inutilizzabili Sauer e lecaratteristiche. Analogamente, dovendo confrontare metodi monodimensionali con me-todi bidimensionali, si e scelto di eseguire il confronto su grandezze significativa perentrambi i casi, che sono il numero di Mach e il rapporto P/P0 mediati sulla sezione.

I confronti effettuati non devono pero trarre in inganno circa la reale potenza deimetodi presentati: nonostante l’ottimo accordo tra quasi-monodimensionale e caratte-ristiche sui valori medi nella maggior parte degli ugelli trattati, quest’ultimo e moltopiu potente in quanto e l’unico che permette la completa risoluzione di tutto il campodi moto bidimensionale (almeno nel tratto divergente). In questo senso le figure 4.39,4.43, 4.47 e 4.51 sono un indice semplice ma chiaro di quanto il valore medio del nu-mero di Mach, possa non rispecchiare le reali condizioni di funzionamento dell’ugello.Per sottolineare ulteriormente le potenzialita del metodo delle caratteristiche, rispettoagli altri visti in questo lavoro, vengono riportate le soluzioni complete sia del campodi velocita (in termini di numero di Mach e angolo di deflessione θ) sia del campo dipressione (in termini di P/P0) per entrambe le tipologie di ugello studiate, cartesiani eassisimmetrici.


11.

21.

41.

61.

82

2.2

2.4

2.6

2.8

30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.91

Car

atte

ristic

he

Figura 4.53: Le caratteristiche dell’ugello No. 3.


11.

21.

41.

61.

82

2.2

2.4

2.6

2.8

30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.91

Figura 4.54: Le caratteristiche dell’ugello No. 10.


Figura 4.55: Numero di Mach all’interno dell’ugello No. 3.


Figura 4.56: Numero di Mach all’interno dell’ugello No. 10.


Figura 4.57: P/P0 all’interno dell’ugello No. 3.


Figura 4.58: P/P0 all’interno dell’ugello No. 10.


Figura 4.59: Angolo di deflessione del flusso theta all’interno dell’ugello No. 3.


Figura 4.60: Angolo di deflessione del flusso theta all’interno dell’ugello No. 10.

Conclusioni

In questo lavoro sono stati confrontati diversi metodi per la soluzione del flusso all’internodi un ugello convergente-divergente. I diversi schemi analizzati non mirano a risolvereil sistema termo-fluidodinamico completo delle equazioni di Navier–Stokes comprimibili3D che e l’unico di grado di descrivere compiutamente il flusso in esame, ma una versioneapparentemente semplificata costituita dalle equazioni di Eulero Comprimibili. Alla basedelle diverse formulazioni risiedono quindi alcune ipotesi semplificative che sono stateenunciate in maniera chiara e completa nella parte introduttiva del lavoro e, in alcunicasi, sono state sottoposte anche a verifica: in particolare sono state sottolineate leeffettive conseguenze che un loro mancato rispetto puo comportare dal punto di vistadell’accuratezza della soluzione.

Uno dei risultati piu interessanti e stato ottenuto in relazione al metodo quasi-monodimensionale la cui potenza si e rivelata decisamente grande. Infatti, pur essendoottenuto a partire da una serie di ipotesi molto restrittive, questo riesce a descrive-re molto bene l’andamento delle grandezze medie lungo il condotto e il confronto conun metodo “esatto”, come quello delle caratteristiche, ci consente di poter affermareche, anche rilassando alcuni vincoli geometrici richiesti a priori nella formulazione quasimonodimensionale, questo schema restituisce effettivamente le quantita integrali sul-la sezione (sottolineando ancora una volta che tali considerazioni sono limitate ai casianalizzati e che rappresentano solo una parte di quelli trattabili con tale metodo).

Stupisce meno il risultato relativo al metodo delle caratteristiche che, seppur limita-tamente ai flussi supersonici, e effettivamente quello piu generale e completo in relazioneal problema di un flusso comprimibile in un ugello. E chiaro che la mancanza di un ef-fettivo confronto non ci permette di dire molto sulla soluzione ottenuta, che dovrebbecomunque essere quella esatta di un gas inviscido supersonico. Positiva si e mostrata adogni modo la versatilita, dal punto di vista implementativo, della formulazione utilizzata,consentendo l’immediata estensione del caso cartesiano a quello assisimmetrico.

Infine e stato possibile riscontre l’estrema limitatezza dei metodi di Sauer e Kliegel:sia in relazione alla geometria dell’ugello, o meglio al parametro di espansione R, siaalla espansione reale del gas, piu grandi sono i gradienti di velocita peggiore e l’accordocon gli altri metodi. Questi risultati sono perfettamente in sintonia con il concetto dilinearizzazione che e alla base di tutti i metodi i quali hanno, dal lato pratico, massimautilita nel calcolo della linea dei dati iniziali necessaria per l’avvio del metodo dellecaratteristiche.

75

Appendice A

Integrazione sulla sezione

Integrazione di una quantita q sulla sezione. Nel caso quasi1D quello che otteniamo egia Q. Nel caso di flusso 2D abbiamo:

Q =(∫

Aq dA

)/ (∫

AdA

)

che nel caso cartesiano, diventa:

Q(x) =

(∫ L/2

−L/2

∫ η(x)

−η(x)q(x, y) dy dx

)/(∫ L/2

−L/2

∫ η(x)

−η(x)dy dx

)

=1

η(x)

∫ η(x)

0q(x, y) dy

nel caso assisimmetrico, in assenza di swirl, invece:

Q(x) =

(∫ 2π

0

∫ η(x)

0q(z, R)R dR dθ

)/(∫ 2π

0

∫ η(x)

0R dR dθ

)

=2η2

∫ η

0q(x,R)R dR

Quasi 1D 2D cartesiano 2D assisimmetrico

Q1η

∫ η

0q(x, y) dy

2η2

∫ η

0q(x,R)R dR

Tabella A.1: Valore medio Q sulla sezione della quantita q.

77

78 Capitolo A. Integrazione sulla sezione

Appendice B

Ugello con curve di Bezier

Si descrive il contorno di un ugello tramite una funzione del tipo y = η(x), che e definitaattraverso l’uso di due funzioni, ovvero:

η(x) =

{ηc(x) : x ≤ xt

ηd(x) : x > xt

Le due funzioni ηc e ηd devono naturalmente soddisfare alcuni vincoli, che sono:

1. ηc(xi) = yi

2. ηc(xt) = yt

3. ηd(xe) = ye

4. ηd(xt) = yt

5. η′c(xt) = 0

6. η′d(xt) = 0

7. η′′c (xt) = η′′d(xt)

8. η′c(xi) = 0

9. η′d(xe) = 0

dove naturalmente le condizioni 8 e 9 non sono strettamente necessarie, ma ci consentonodi ottenere un ugello piu regolare. La condizione numero 9 e invece necessaria peravere un raggio di curvatura continuo in relazione al metodo di Sauer. Questo tipo dicondizioni sono assolutamente generali e ora andiamo a introdurre le curve di Bezier.

x(s) =N∑

i=1

XiBi s ∈ [0, 1] (B.1)

79

80 Capitolo B. Ugello con curve di Bezier

con N il grado del polinomio di Bezier,Xi = [Xi, Yi] i gradi di liberta associati, e Bi lefunzioni dette di mescolamento, definite come segue per un polinomio di grado 4:

B1 = (1− s)3

B2 = 3s(1− s)2

B3 = 3s2(1− s)

B4 = s3

e per uno di grado 5

B1 = (1− s)4

B2 = 4s(1− s)3

B3 = 6s2(1− s)2

B4 = 4s3(1− s)

B5 = s4

Figura B.1: L’ugello di riferminento, sono segnati anche i punti X.

I punti X1 e XN sono punti particolari, poiche rappresentano i punti di ancoraggiodella funzione, ovvero x(0) = X1 e x(1) = XN . Poiche vogliamo descrivere l’ugello infunzione di x, e possibile confondere s con x/xend, e analogamente la derivata dη

ds conxend

dηdx , dove xend sar rispettivamente xt per ηc e xe per ηd.

Il numero dei vincoli da imporre e 9, a cui ne vanno aggiunti altri 4 per potergarantire che s = x/xend ∈ [0, 1], e la particolarita di questi 9 ci obbliga alla scelta di unpolinomio di grado N = 9 o meglio alla scelta di uno di grado N = 4 e l’altro N = 5. Ilrisultato di questa operazione e l’ottenimento di una funzione η(x) che soddisfa tutti irequisiti e che pero ci lascia 5 ulteriori gradi di libera per poter manipolare la sua forma, ecorrispondono alla componente xi dei punti Xi necessari per costruire le curve di Bezier.

81

Traslando questi cinque punti, che corrispondono a quelli interni, lungo l’orizzontale epossibile deformare l’ugello a piacere. Va posto l’accento sul fatto che non sempre sisi ottiene una forma accettabile: e infatti possibile che si generi piu di un punto conderivata prima nulla all’interno dell’ugello, con quindi piu di un’area di gola, oppure sipossono ottenere forme che non definiscono una funzione η(x), ma rappresentano casidecisamente rari e comunque non interessanti ai fini di cio che concerne gli ugelli.

82 Capitolo B. Ugello con curve di Bezier

Bibliografia

[1] J. D. Hoffman, M. J. Zucrow, Gas Dynamics Vol I, Robert E. KriegerPublishing Company, (1977).

[2] J. D. Hoffman, M. J. Zucrow, Gas Dynamics Vol II, Robert E. KriegerPublishing Company, (1977).

[3] F. Sabetta, Gasdinamica, Edizioni ingegneria 2000.

[4] J. R. Kliegel, V. Quan, Convergent-Divergent Nozzle Flows, N.A.S.A.-MannedSpacecraft Center, (1966).

[5] J. B. Callen, Thermodynamics and an introduction to thermostatistic, aaa,(1982).

83

Gasdinamica 2 Progetto del corso -...

Documents

Transcript of Gasdinamica 2 Progetto del corso -...