Cfd Per Valvole Di Sicurezza

Dipartimento di Meccanica ed Aeronautica TR-ISPESL-010505

Linee guida per l'applicazione della fluidodinamica computazionale (CFD) per la caratterizzazione dell'efflusso nei dispositivi di protezione e regolazione di attrezzature in pressione. Realizzazione di un modello su geometrie assegnate, al fine di estendere i risultati ottenuti in prove sperimentali a diverse tipologie di fluidi, campi di pressione e condizioni di esercizio. [B120/DOM/03]

Paolo Bernardini Enrico Sciubba Giugno 2005

1 - Descrizione del progetto e piano delle attività.

1. Descrizione del progetto e piano delle attività. 1.1. Descrizione del progetto. Lavoro di ricerca svolto dalla Cattedra di Turbomacchine del Dipartimento di Meccanica e Aeronautica dell’Università di Roma “La Sapienza” per conto di ISPESL nell’ambito dell’AREA TEMATICA DI RICERCA n. 9: Affidabilità, manutenzione e conservazione. Verifica dell'efficacia di applicazioni tecnologiche e normative di prevenzione infortuni su impianti, macchine, strutture, attrezzature e sistemi. Armonizzazione delle prassi tecnologiche dell'Istituto con i criteri imposti dalla evoluzione normativa nazionale ed internazionale. Motivazioni: L'esigenza di disporre di metodi che consentano di estendere i risultati ottenuti con le prove sperimentali a diverse tipologie di fluidi e campi di pressione, limitando i costi ed i tempi che dette prove comportano, richiede lo sviluppo di un modello fluido-termodinamico di tipo numerico. Obiettivi. L'attività di ricerca si propone la messa a punto di un modello di calcolo fluido-termodinamico che, applicato a geometrie reali di dispositivi di protezione e regolazione, permetta di simulare l'efflusso e prevedere il comportamento del dispositivo in esame a condizioni non sperimentate.

1.2. Piano delle attività e scadenze. Durata del progetto: annuale Assegnazione: inizio 2003. Inizio dei lavori: giugno 2004. Termine per la consegna: giugno 2005.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 2

2 – Scelta delle geometrie.

2. Scelta delle geometrie rappresentative di apparati reali (“geometrie tipo”). Dati sperimentali.

2.1. Geometrie tipo. L’ISPESL ha messo a disposizione diverse geometrie di valvole di sicurezza a molla di diversi produttori italiani del tipo a scarico libero e a scarico convogliato. Si sottolinea che tutte le geometrie a disposizione caratterizzano tipologie di valvole ampiamente diffuse in commercio e normalmente utilizzate nelle apparecchiature industriali. Di seguito sono illustrati alcuni disegni schematici delle due tipologie insieme ad alcune foto.

Valvole a scarico libero (con e senza cappellotto di protezione):

Foto 2.1.1 – Valvole a scarico libero

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 3

2 – Scelta delle geometrie.

Foto 2.1.2

Figura 2.1. 1 – Sezione della valvola in Foto 2.1.2

Foto 2.1.3 Figura 2.1. 2 – Sezione della valvola in Foto 2.1.3

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 4

2 – Scelta delle geometrie.

Figura 2.1. 3 – Vista e sezione della valvola in Foto 2.1.4

Foto 2.1.4 Foto 2.1.5 Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 5

2 – Scelta delle geometrie.

Figura 2.1. 4 – Vista e sezione della valvola in Foto 2.1.6

Foto 2.1.6 Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 6

2 – Scelta delle geometrie.

Valvole a scarico convogliato:

Foto 2.1.7 – Valvole a scarico convogliato

Foto 2.1.8

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Figura 2.1. 5 – Sezione della valvola in Foto 2.1.8

Roma, 16/06/2005

Pagina 7

2 – Scelta delle geometrie.

Dopo una valutazione preliminare, in accordo con l’ISPESL, è stata scelta la geometria della Figura 2.1. 3 e Foto 2.1.4 e Foto 2.1.5. Si tratta di una valvola (denominata Modello D-14) a carico diretto, a sede conica, con molla, e a scarico libero in atmosfera. Presenta due feritoie di scarico diametralmente opposte. La scelta della suddetta geometria è dettata principalmente dalla natura del presente lavoro: essendo questo un progetto pilota, è finalizzato a definire la base per una metodologia innovativa di valutazione delle prestazioni di questo tipo di dispositivi. Si è individuata quindi la tipologia più diffusa dal punto di vista prettamente geometrico e che offre un esempio chiaro e completo nella realizzazione della griglia tridimensionale di calcolo e nell’analisi numerica e fisica dei risultati. Infatti dall’analisi delle diverse geometrie, prevedendo per quanto possibile il moto dell’aria attraverso esse, si è stabilito che nelle altre valvole, l’efflusso sarebbe stato più semplice nel caso delle valvole del tipo in Foto 2.1.6, o con interazioni fluidodinamiche tra fluido e pareti più accentuate come nel caso delle valvole a scarico convogliato della Foto 2.1.7. Nella valvola presente si prevedeva (considerando l’efflusso da una delle due feritoie) un flusso in uscita verosimilmente composto da un unico getto radiale idealmente divisibile in due flussi principali controrotanti simmetrici rispetto ad un piano passante per l’asse della valvola e per la mezzeria della feritoia. Questo lavoro ha tra i suoi scopi quello di valutare se valvole analoghe potranno essere “simulate” con lo stesso approccio numerico. Sottoliniamo la necessità di una simulazione fluidodinamica tridimensionale, non essendo l’efflusso riconducibile ad una trattazione bi-dimensionale per via del suo sviluppo naturale nelle tre dimensioni spaziali. I risultati (vedi capitoli 7 e 8) confermano la validità di questa scelta, e mostrano la natura altamente tridimensionale dell’efflusso soprattutto all’interno della valvola. 2.2. Dati sperimentali. Oltre alle geometrie delle valvole suddette, l’ISPESL ha fornito le relazioni sulle prove sperimentali atte a rilevare le caratteristiche di funzionamento e di portata delle valvole di sicurezza. Dalla relazione, ed in particolare dai rapporti di prova, si sono ricavati i valori assunti dalle grandezze fisiche nelle condizioni di prova. I dati utili per le simulazioni sono i seguenti: 1. pressione relativa nel serbatoio di prova; 2. temperatura dell’aria nel serbatoio; 3. pressione ambiente (atmosferica); 4. alzata1 dell’otturatore nelle condizioni delle “prove di qualifica”2; 5. coefficiente di efflusso; da questo e dalla portata teorica si ricava la portata in massa (non necessaria per il calcolo ma indispensabile per il confronto con i valori calcolati).

1

Nota bene: tutte le simulazioni (quindi per tutti i valori della pressione relativa nel serbatoio) sono state effettuate in relazione ai dati sperimentali delle “prove di qualifica” effettuate sulla valvola avente l’otturatore con l’alzata bloccata a 5mm (rispetto alla sede). 2 Come specificato nella Circolare 15 novembre 1979, n.38468 – Raccolta E, le “prove di qualifica” servono alla determinazione della portata di efflusso all’alzata h, ai fini di accertare il coefficiente di efflusso del tipo di valvola. Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 8

2 – Scelta delle geometrie.

I valori non presenti nei rapporti di prova e necessari per le simulazioni numeriche sono stati scelti arbitrariamente e sulla base di casi simulati in studi precedenti e di ciò che mette a disposizione la letteratura scientifica. Essi sono la temperatura dell’aria (nell’ambiente esterno al serbatoio), le temperature delle superfici metalliche della valvola e del serbatoio e i valori di grandezze fisiche correlate alla turbolenza dell’efflusso. Come vedremo in seguito, si sono ipotizzati valori in base allo scambio termico presunto nelle condizioni di prova, e si sono attribuiti valori di temperatura sulle pareti suddette pari o molto vicini a quelle dell’aria nel serbatoio. Si sono imposti valori decrescenti (gradienti dell’ordine di qualche grado/metro) sulle superfici esterne della valvola. Le temperature delle superfici interne della valvola sono invece stati assunti uguali a quella del serbatoio. Anche la rugosità superficiale delle pareti solide a contatto con il fluido è stata scelta arbitrariamente, ipotizzando la tipologia di lavorazione meccanica utilizzata per l’ottenimento dei particolari costruttivi della valvola. Al capitolo 6, paragrafo 6.4 Condizioni al contorno, si riportano i valori attribuiti alle grandezze correlate alla turbolenza.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 9

3 - Le equazioni di governo del moto di un fluido.

3. Le equazioni di governo del moto di un fluido.

3.1. Introduzione Il moto di un mezzo continuo è governato dai principi della meccanica e della termodinamica classica. Esso è rappresentato mediante le equazioni esprimenti le leggi della conservazione della massa, della quantità di moto e dell’energia. Nell’applicazione di questi principi ci si avvale della descrizione Euleriana del moto in un sistema di riferimento Galileiano (assoluto). Il punto di vista Euleriano fissa una data posizione x,y,z ed osserva, al trascorrere del tempo, quel che accade in tale posizione. Le variabili indipendenti sono quindi le x,y,z (in notazione compatta xi ) ed il tempo t. Dato che il presente studio avviene in condizioni stazionarie la variabile tempo non comparirà esplicitamente. Le proprietà caratteristiche del mezzo fluido sono considerate quindi come funzioni dello spazio nel sistema di riferimento. Il mezzo fluido è ritenuto continuo: questa assunzione implica che esistono le derivate di tutte le variabili dipendenti; in altre parole, proprietà locali come la densità e la velocità sono definite come medie su elementi “grandi” se comparati con la struttura microscopica del fluido, ma abbastanza “piccoli” in confronto alla scala dei fenomeni macroscopici. Ciò permette di descriverli con l’uso del calcolo differenziale. Laddove necessario, per tenere conto della natura reale e turbolenta dell’efflusso il fluido è considerato viscoso e si introducono modelli ad una o più equazioni per simulare le interazioni turbolente. Nei paragrafi seguenti vengono presentate le diverse forme delle equazioni di conservazione nella loro applicazione ai fluidi incompressibili e compressibili.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 10

3 - Le equazioni di governo del moto di un fluido.

3.2. Equazione di continuità o di conservazione della massa. Il principio di conservazione della massa, nel caso del moto di un fluido, si esprime dicendo che resta invariata nel tempo la massa contenuta in un volume che si muove insieme al fluido. Si scrive allora: dM d dt = dt

∫

ρdV = 0

(3.2.1)

V (t )

Applichiamo ora il teorema del trasporto di Reynolds. Quest’ultimo permette di determinare la derivata temporale dell’integrale di una certa quantità A(xi,t) esteso ad un volume arbitrario che si muove con il fluido. Nel caso generale si scrive: d dt

∫

AdV =

V (t )

∫

V (t )

∂A ∂t dV +

∫

Au ⋅ nds

(3.2.2)

S (t )

Bisogna ora applicare all’integrale superficiale della (3.2.2) il teorema della divergenza:

∫

a ⋅ udS =

S

∫

∇⋅ adV

(3.2.3)

V

Il teorema della divergenza stabilisce che il flusso di un vettore uscente da una superficie chiusa è uguale all’integrale della divergenza del vettore stesso, esteso al volume racchiuso. Abbiamo così: d dt

∫

AdV =

V (t )

∂A { ∫ ∂t + ∇ ⋅ ( Au)}dV

(3.2.4)

V (t )

Osservando che: ∇ ⋅ ( Au) = (u ⋅∇) A + A∇ ⋅ u

(3.2.5)

ed utilizzando l’espressione della derivata sostanziale: dA ∂A dt = ∂t + ( u ⋅ ∇) A

(3.2.6)

si ottiene: d dt

∫

AdV =

V (t )

dA { ∫ dt + A∇ ⋅ u}dV

(3.2.7)

V (t )

La derivata sostanziale di una grandezza fisica esprime la variazione totale nel tempo della grandezza stessa percepita da un osservatore solidale al moto della particella fluida. Nell’espressione (3.2.6) possiamo individuare due termini a secondo membro; il primo, detto derivata locale, esprime a livello fisico la variazione nel tempo della grandezza in un punto fissato; il secondo termine, detto derivata convettiva, equivale alla variazione temporale dovuta al movimento dell’elemento fluido da un punto all’altro di un campo fluido dove le proprietà del flusso sono diverse nello spazio. Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 11

3 - Le equazioni di governo del moto di un fluido.

Riprendendo la (3.2.1) nella forma (3.2.4 ) o nella forma (3.2.7) e, inserendo per la generica A(xi,t) la densità ρ, risulta che devono essere nulli gli integrali a secondo membro “qualunque sia il volume di integrazione V(t)”; devono essere quindi identicamente nulli gli integrandi:

∂ρ

∂t + ∇ ⋅ ( ρu) = 0

(3.2.8)

dρ + ρ∇ ⋅ u = 0 dt

(3.2.9)

oppure in forma equivalente:

La (3.2.8), o la (3.2.9), costituisce la “equazione di continuità” per un fluido. La forma in cui è espressa la (3.2.8), viene chiamata “conservativa”; la forma della (3.2.9) è detta “non conservativa”. Le (3.2.8) e (3.2.9) possono essere scritte nelle rispettive forme esplicite:  ∂u ∂v ∂w  ∂ρ ∂ρ ∂ρ ∂ρ  = 0 + ρ  + + +u +v +w ∂z ∂t ∂x ∂y  ∂x ∂y ∂z 

(3.2.8a)

dρ ∂u ∂v ∂w +ρ +ρ +ρ =0 dt ∂x ∂y ∂z

(3.2.9a)

Se il moto è stazionario, ossia se le proprietà locali non variano nel tempo, la (3.2.8) si riduce a: ∇ ⋅ ( ρu) = 0

(3.2.10)

Le (3.2.8) e la (3.2.9) sono valide nel caso generale di fluido compressibile e si equivalgono solo se la densità ρ non è funzione della posizione. dρ Altro caso particolare è quello di un fluido incompressibile, per il quale è dt = 0 ; la (3.2.9) fornisce in questo caso: ∇⋅u = 0

⇔

∂ui ∂xi = 0

(3.2.11)

1 dV La (3.2.11), essendo V dt = divu = ∇ ⋅ u , ci dice che la velocità di variazione relativa del volume della particella fluida è nulla. Il fluido è quindi incompressibile.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 12

3 - Le equazioni di governo del moto di un fluido.

3.3. Equazione della conservazione della quantità di moto L’equazione della quantità di moto applicata alla massa contenuta in un certo volume di fluido, che si muove con esso, si scrive:

dQ dt = Fe

(3.3.1)

dove Fe è la risultante delle forze esterne di massa (ad es. la gravità, la forza centrifuga) e di quelle di superficie, come gli sforzi dovuti al fluido esterno e agenti sulla superficie S di contorno. La (3.3.1) diventa allora: d dt

∫

ρudV =

V (t )

∫

ρ f dV +

V (t )

∫

T ⋅ ndS

(3.3.2)

S (t )

dove f è la forza di volume che si esercita per unità di massa e T è il tensore degli sforzi. Tale tensore permette di descrivere gli sforzi attorno ad un punto nelle varie direzioni possibili: esso è un tensore a nove componenti scalari (a tre componenti vettoriali rispetto alle tre direzioni degli assi coordinati prescelti). Rimane da sottolineare che T è un tensore simmetrico per cui Tij = Tji , cosicché le nove componenti si riducono a sei quantità indipendenti. Occorre ora esprimere il principio di conservazione della quantità di moto per un volume che non varia nel tempo, cioè un volume fisso nello spazio (formulazione Euleriana). Ciò può essere semplicemente realizzato esprimendo le derivate temporali degli integrali sul volume materiale V ( t ) mediante il teorema del trasporto di Reynolds. Trasformiamo così il primo membro e il secondo termine del secondo membro della (3.3.2) in modo che vi appaiano soltanto integrali di volume, come è il primo termine a secondo membro. Se poi applichiamo il teorema della divergenza:

∫

∫

T ⋅ ndS = divTdV

S

(3.3.3)

V

all’integrale superficiale avremo:

∫

  du ρ dt − ρ f − ∇ ⋅ T dV = 0  V

(3.3.4)

Dovendo valere la (3.3.4) per qualsiasi volume di integrazione, l’integrando deve essere identicamente nullo; quindi: du ρ dt = ρ f + ∇ ⋅ T

(3.3.5)

Passando dalla espressione vettoriale a quella delle componenti ed introducendo la simbologia della somma introdotta da Einstein, o notazione indiciale1: 1

In ogni prodotto di termini, l’indice ripetuto implica una somma rispetto allo stesso indice per i valori

1,2,3. Quello non ripetuto può assumere uno qualsiasi dei valori 1,2,3.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 13

3 - Le equazioni di governo del moto di un fluido.

du ∂ ρ dti = ρfi + ∂x Tij j

(3.3.6)

dove u i è la componente i-esima della velocità istantanea, ρ è la densità, Tij è il tensore degli sforzi e fi è la componente i-esima della forza di volume per unità di massa. Dalle ipotesi fatte il fluido è di tipo Newtoniano e a comportamento isotropo. Per i “fluidi Newtoniani” le componenti del tensore degli sforzi sono funzioni lineari delle componenti delle velocità di deformazione. Per le componenti del tensore degli sforzi si dimostra che è: Tij = 2µεij per (i ≠ j)

(3.3.7)

Tjj = ( − p + λ ∇ ⋅ u) + 2µε jj (senza somma su j)

(3.3.8)

Dalla (3.3.8) risulta che nel caso generale i tre sforzi normali sono diversi tra loro, e la loro media non coincide con la pressione. Nelle (3.3.7) e (3.3.8) µ è la viscosità dinamica che riteniamo costante, ed εij sono le componenti del tensore delle deformazioni. La relazione tra sforzi e velocità di deformazione si scrive: 2   Tij = 2 µε ij −  p + µ ∇ ⋅ u δ ij 3  

(3.3.9)

dove p è la pressione termodinamica e δ ij è il delta di Kronecker: 1 per i = j δij =  0 per i ≠ j

Nella (3.3.8) il termine λ∇ ⋅ u descrive l’effetto della viscosità dovuto alla variazione di volume di una particella fluida. I due coefficienti di viscosità λ e µ sono legati tra loro dalla relazione:

2 µ′ = λ + 3 µ

(3.3.10)

0 < µ ′ << µ

(3.3.11)

2 λ = −3µ

(3.3.12)

Essendo per i gas poliatomici:

possiamo scrivere:

La quantità µ ′ , detta “bulk viscosity” o viscosità di massa, descrive la differenza esistente, e dovuta alla viscosità, tra sforzo normale medio e pressione in un fluido in espansione. In altre parole, supponendo di avere una massa di gas viscoso che si espande rapidamente (se ∇ ⋅ u > 0 ), il suo comportamento coincide con quello di un gas non viscoso a pressione p′ inferiore. Si dimostra, infatti, che è: p′ = p − µ ′ div u

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

(3.3.13)

Roma, 14/06/2005

Pagina 14

3 - Le equazioni di governo del moto di un fluido.

Questo significa che un gas viscoso, sottoposto a pressione decrescente nel tempo, si espande meno rapidamente di un gas non viscoso sottoposto alla stessa legge temporale delle pressioni esterne e parità delle altre condizioni. La (3.3.9) è detta “equazione costitutiva”. Come si vede, per i fluidi isotropi Newtoniani, essa dipende dai due coefficienti di viscosità λ e µ e, in virtù della (3.3.12), dal solo coefficiente µ. A questo punto se introduciamo nelle (3.3.6) l’equazione costitutiva (3.3.9) e le componenti del “tensore di deformazione”: 1  ∂u ∂u j  εij = 2  ∂x i + ∂x   j i

(3.3.14)

otteniamo le “equazioni di Navier-Stokes” in forma differenziale:

du ∂p ∂ ∂   ∂u ∂u j   ρ dti = ρfi − ∂x + ∂x ( λ ∇ ⋅ u) + ∂x µ ∂x i + ∂x   i i j  i    j 

(3.3.15)

L’ultimo termine della (3.3.15) può essere scritto: ∂µ  ∂u ∂u j  ∂ µ∇2ui + µ ∂x ∇ ⋅ u + ∂x  ∂xi + ∂x  i j j i

(3.3.16)

Nel caso in cui µ e λ possono essere ritenuti costanti la (3.3.15) diventa così:  ∂u du ∂u  ∂p ∂ ρ dti = ρ ∂ti + u j ∂x i  = ρfi − ∂x + µ∇2ui + ( λ + µ) ∂x (∇ ⋅ u)  j i i

(3.3.17)

La (3.3.17) può essere scritta in forma vettoriale ricordando però che la somma ( λ + µ ) è, 1 per via della (3.3.12), pari a 3 µ : ρ

µ du = ρ f − ∇ p + µ∇2 u + ∇(∇ ⋅ u) dt 3

(3.3.18)

Nel caso particolare di fluido incompressibile si aggiunge la seguente:

∇⋅u = 0 e la (3.3.18) diventa: ρ

du = ρ f − ∇ p + µ∇ 2 u dt

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

(3.3.19)

Pagina 15

3 - Le equazioni di governo del moto di un fluido.

3.4. Equazione dell’energia Il primo principio della termodinamica, definita l’equivalenza tra energia e calore, stabilisce il principio della conservazione dell’energia:

dE

= L −Q

(3.4.1)

dt Nella (3.4.1) la dE è l’energia totale del sistema (la particella fluida):

E = ∫ ρ edV

(3.4.2)

V

essendo “e” l’energia totale locale per unità di massa. Quest’ultima è pari alla somma dell’energia interna U (molecolare), posseduta in virtù del moto di agitazione “termica” attorno u2 al baricentro, e della energia cinetica macroscopica 2 per unità di massa, dovuta al moto di insieme della particella fluida. L è il lavoro fatto nell’unità di tempo, sul sistema, dalle forze esterne di volume e di superficie (sforzi). La Q equivale al calore che fluisce, per unità di tempo, verso l’esterno del sistema. Il principio della conservazione dell’energia può scriversi: d dt

∫

V

ρEdV =

∫

ρ f ⋅ udV +

V

∫(

)

T ⋅ u − K ⋅ ndS

S

(3.4.3)

dove il lavoro L è così diviso in un lavoro della forze di volume L′ ed un lavoro degli sforzi superficiali L′′ . Il primo termine del secondo integrale a secondo membro esprime il lavoro degli sforzi viscosi; infatti:

(T ⋅ u) ⋅ ndS = (T ⋅ ndS) ⋅ u = T dS ⋅ u = d R ⋅ u

(3.4.4)

n

dove d R è la risultante della forza totale agente sull’elemento di superficie dS ad opera delle molecole situate sulla faccia opposta dell’areola dS ; quindi il lavoro totale che le forze d R compiono per unità di tempo è proprio pari al prodotto d R ⋅ u . Il secondo termine del secondo integrale a secondo membro della (3.4.3) equivale all’energia trasmessa, per unità di tempo, attraverso la superficie dS e legata alla rapidità di variazione spaziale della temperatura (gradienti termici) nella direzione normale n ; esprime quindi il flusso termico. Infatti la:

[

]

− ( K ⋅ n)dS = ( k ∇T ) ⋅ n dS

(3.4.5)

esprime la legge di Fourier, dove il segno meno puntualizza il fatto che tale flusso di calore è positivo nel verso in cui le temperature sono decrescenti. E’ da notare che la conducibilità termica k del fluido è stata espressa mediante uno scalare.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 16

3 - Le equazioni di governo del moto di un fluido.

Applicando ora alla (3.4.3) il teorema del trasporto di Reynolds al primo membro e il teorema della divergenza all’integrale superficiale, si ottiene l’equazione della conservazione dell’energia: d dt

∫

ρEdV =

V

∫

 d ( ρE ) + ρE ∇ ⋅ udV  dt  V

(3.4.6)

ossia d dt

∫

ρEdV =

V

∫

ρ

V

∫

 dρ  E  + ρ∇ ⋅ udV dt  V 

dE dV + dt

(3.4.7)

Ma l’ultimo integrale è nullo in conseguenza della continuità (3.2.8), e utilizzando il secondo membro della (3.4.3) si ottiene:

(

∫

ρ dE − ρ f ⋅ u − ∇ ⋅ T ⋅ u − K  dt V

)dV = 0

(3.4.8)

Da questa, dovendo essere nullo l’integrando e poiché è del tutto arbitrario il volume di integrazione, si scrive:

( )

dE ρ dt = ρ f ⋅ u + ∇ ⋅ T ⋅ u − ∇ ⋅ K

(3.4.9)

Volendo esprimere l’equazione della conservazione dell’energia in termini di variazione dell’energia interna, bisogna combinare la (3.4.6) con la equazione relativa ai soli scambi di energia meccanica, ottenibile dal lavoro delle forze che figurano nell’equazione della quantità di moto (3.3.5): du ρ dt = ρ f + ∇ ⋅ T

Da questa infatti si scrive:

( )

du d  u2  ρ dt ⋅ u ≡ ρ dt  2  = ρ f ⋅ u + ∇ ⋅ T ⋅ u  

(3.4.10)

u2 e sottraendo la (3.4.7) dalla (3.4.6), tenuto conto della e = U + 2 , si ha allora:

( ) ( )

dU ρ dt = ∇ ⋅ T ⋅ u − ∇ ⋅ T ⋅ u − ∇ ⋅ ( − k ∇T )

(3.4.11)

A questo punto, sostituendo al tensore degli sforzi l’espressione in funzione delle velocità di deformazione, si ottiene:  p Dρ  DU ρ Dt = εijTij − ∇ ⋅ ( k ∇T ) =  ρ Dt + µΦ + ∇ ⋅ ( k ∇T )  

(3.4.12)

∑∑

(3.4.13)

dove è: Φ = 2εijεij = 2

i

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

j

εij2 −

2 3

∑

εii2

i

Roma, 14/06/2005

Pagina 17

3 - Le equazioni di governo del moto di un fluido.

La (3.4.13), tenendo conto della (3.3.12), diventa (esplicitata) per il caso bidimensionale:

 ∂u  2  ∂v  2  ∂w  2 1  ∂u ∂v ∂w  2  2  ∂u ∂v ∂w  2    −  + + Φ = 2   +   +   +  + +  ∂x   ∂y   ∂z  2  ∂y ∂x ∂z   3  ∂x ∂y ∂z 

(3.4.14)

Dρ La (3.4.12), in virtù della equazione di continuità Dt + ρ∇ ⋅ u = 0 , diventa:

dU ρ dt = − p∇ ⋅ u + µΦ + ∇ ⋅ ( k ∇T )

(3.4.15)

e scritta in forma esplicita per un caso bidimensionale:

 ∂U ∂U ∂U  + +u +v ρ  ∂x ∂y   ∂t

 ∂ 2T ∂ 2T   ∂u ∂v  p +  = k  2 + 2  + µΦ ∂y   ∂x ∂y   ∂x

(3.4.15a)

La quantità Φ è detta “funzione di dissipazione”: essa si può dimostrare essere definita positiva, e rappresenta la frazione di energia meccanica che si trasforma in energia termica e che viene dissipata sotto questa forma. La Φ è quindi pari al tasso di generazione di calore dovuto alla dissipazione viscosa. La (3.4.15) può essere scritta anche nella forma seguente: µ dh 1 dp 1 = + ∇ ⋅ k ∇ ⋅ T + ( ) dt ρ dt ρ ρΦ

(3.4.16)

p dρ dp d  p εij Tij = ρ dt + µΦ = dt − ρ dt  ρ  + µΦ  

(3.4.17)

essendo:

e avendo introdotto l’entalpia del fluido: p h =U + ρ

(3.4.18)

Ancora, la (3.4.15) può essere scritta anche in termini di entropia introducendo la definizione di questa grandezza nell’espressione differenziale del primo principio della termodinamica: dp dh = ρ + Tds

(3.4.19)

ed ottenendo così, confrontando la (3.4.16) con la (3.4.19):

ds 1 T dt = ρ ∇ ⋅ ( k ∇T ) + νΦ

(3.4.20)

µ dove è stata introdotta la viscosità cinematica ν = ρ e dh, dp, ds sono rispettivamente le variazioni di entalpia, pressione ed entropia di una particella fluida. Il secondo membro della (3.4.20) dà i contributi della conduzione termica e della dissipazione per attrito alla variazione temporale dell’entropia. In particolare, si nota come l’ultimo termine (dissipativo) dà un contributo definito positivo alla variazione di entropia.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 18

3 - Le equazioni di governo del moto di un fluido.

3.5. Fluidi incompressibili Se si assume che il fluido sia incompressibile, si ottengono considerevoli semplificazioni. Di solito tale assunzione è ritenuta ragionevole quando la velocità caratteristica del fluido è minore di 0.3 volte la velocità del suono alle condizioni di temperatura e pressione a cui si trova il fluido stesso. Per un fluido incompressibile vale la seguente relazione (dalla (3.2.8)): ∇⋅u = 0

⇔

∂ui =0 ∂xi

(3.5.1)

In questo caso l’espressione della conservazione dell’energia (3.4.15), con l’introduzione della: dU = ρc pdT

(3.5.2)

diventa: ρc p

dT = ∇ ⋅ (k ∇T ) + µΦ dt

(3.5.3)

con cp pari al calore specifico a pressione costante. 3.6. Fluidi compressibili Per fluidi compressibili l’equazione della conservazione dell’energia assume la forma seguente: ρcv

dT = − p∇ ⋅ u + µΦ + ∇ ⋅ (k ∇T ) dt

(3.6.1)

ottenuta direttamente dalla (3.4.15) essendo cv il calore specifico a volume costante. Per un gas ideale vale la relazione: c p − cv = R

(3.6.2)

dove R è la costante del gas in questione. 3.7. Equazione di stato Per completare l’impostazione del problema è necessario aggiungere un’ulteriore relazione tra le variabili indipendenti, l’equazione di stato, e cioè una relazione termodinamico-fisica tra densità, temperatura, pressione. Viene anche richiesto un parametro caratteristico del fluido in questione mediante l’immissione del valore della costante del gas. Esiste però una differenza nella forma dell’equazione di stato dipendente dalla considerazione o dall’esclusione della compressibilità del fluido.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 19

3 - Le equazioni di governo del moto di un fluido.

3.8. Fluidi incompressibili Per i fluidi incompressibili l’equazione di stato assume semplicemente la forma seguente: ρ = ρ0

(3.8.1)

cioè, la pressione non compare nella equazione e la densità rimane costante. E’ opportuno sottolineare come l’ipotesi di incomprimibilità non sia equivalente a quella di densità costante; la densità, infatti, è in generale funzione della temperatura oltre che della pressione. 3.9. Fluidi compressibili L’equazione di stato per un gas ideale e omogeneo si scrive: pv =

m m R0 R0T ⇒ p = T ⇒ p = ρRT M v M

(3.9.1)

dove M è il peso molecolare del gas ed R0 è la costante universale dei gas. Nelle simulazioni (vedi capitolo 6) eseguite è stata considerata la comprimibilità, ma ne sono state trascurate le variazioni con la temperatura.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 20

3 - Le equazioni di governo del moto di un fluido.

3.10. Condizioni al contorno e condizioni iniziali L’equazione di continuità (3.2.8), l’equazione della quantità di moto (3.3.15), equivalente a tre equazioni scalari, e l’equazione dell’energia (3.6.1), costituiscono insieme all’equazione di stato un sistema di sei equazioni differenziali nelle sei incognite u, v, w, ρ , p, T . Esse possono essere risolte, analiticamente o numericamente, con opportune condizioni al contorno. In realtà la soluzione analitica delle equazioni di Navier-Stokes (in genere è così chiamato tutto il sistema di equazioni) presenta in generale difficoltà insormontabili, principalmente per il fatto che le equazioni stesse sono alle derivate parziali e non lineari. Da questo risulta l’importanza e la necessità dell’impiego di solutori numerici che possono sfruttare la velocità di calcolo e la precisione del calcolatore elettronico. La risoluzione del sistema necessita della specifica delle cosiddette condizioni al contorno, e delle condizioni iniziali. Sono proprio queste condizioni che decidono le soluzioni da ottenere dalle equazioni di governo. Su ciascuna “linea” o “superficie” di confine del dominio di calcolo Ω è necessario specificare appropriate condizioni. Dato che le equazioni di Navier-Stokes stazionarie sono di tipo ellittico, se il fluido è incompressibile o compressibile con valori del numero di Mach minori dell’unità, la loro soluzione necessita della specificazione di due condizioni al contorno per ogni coordinata, ed una condizione iniziale (tranne che per la pressione per la quale il valore si ricava, senza l’imposizione di condizioni esplicite, a meno di una costante). L’ellitticità delle equazioni implica il fatto che una variazione del valore di una condizione al contorno in un punto qualsiasi del bordo modifica istantaneamente la soluzione in tutto il dominio di calcolo.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 21

4 – La fisica del fenomeno dell’efflusso in una valvola.

4. La fisica del fenomeno dell’efflusso in una valvola.

4.1. Generalità La geometria ed il principio di funzionamento della tipologia delle valvole oggetto di questo lavoro sono molto semplici. Più complessa risulta invece la previsione e l’analisi approfondita dell’efflusso. Ci riferiremo alla valvola oggetto delle simulazioni. La parte mobile della valvola è l’otturatore (vedi Figura 4.1.1). La molla tarata esercita la sua forza elastica sull’otturatore che è premuto, con l’interposizione di una Gruppo guarnizione di tenuta, contro la sede. In tal modo, si evita stelo-molla la fuoriuscita del gas o vapore contenuto all’interno del e serbatoio da controllare. Il fluido esercita la sua pressione vite di regolazione della taratura sulla superficie della guarnizione limitata internamente dal ciglio della sede. Sull’area di questa superficie si calcola la spinta necessaria della molla per la corretta taratura (ovviamente non si può prescindere dalla taratura al banco). Otturatore con Quando si raggiunge o si supera la pressione di guarnizione di taratura (anche detta impropriamente “di sfioro”), tenuta l’otturatore inizia a sollevarsi. Tanto maggiore è la sovrapressione (rispetto alla pressione di taratura) all’interno del serbatoio, tanto maggiore è la spinta sull’otturatore e l’altezza raggiunta dall’otturatore rispetto Corpo con: • sede integrata alla sede. •attacco filettato In funzione della taratura della molla, delle portate •condotto d’entrata da smaltire e dell’andamento temporale della pressione •feritoie di scarico nel serbatoio (dipendente dalle condizioni termodinamiche al suo interno e dalla sua capacità), la legge di alzata dell’otturatore può differire in modo sostanziale. La forma dell’otturatore, ed in particolare il suo profilo laddove il fluido Figura 4.1.1 viene da esso guidato al momento del distacco, può comportare differenze nella spinta risultante. Questa sua caratteristica, insieme al suo sviluppo diametrale, gioca un ruolo importante soprattutto alla richiusura del gruppo otturatore e nel fenomeno del blowdown (differenza che si riscontra fra il valore della pressione relativa all’apertura e alla chiusura). In effetti in condizioni dinamiche, il fluido uscente dalla sede esercita una spinta su di una sezione maggiore. Si deve inoltre evidenziare come soprattutto all’apertura e/o alla richiusura, si possono manifestare fenomeni di chatter o flutter che consistono in oscillazioni più o meno marcate dell’otturatore che possono comportare anche urti continuati contro la sede. Tali fenomeni sono direttamente collegati alla natura turbolenta dell’efflusso e ad un campo fluidodinamico nel quale le variazioni e le oscillazioni di velocità e pressione locali comportano squilibri ed instabilità. Immaginiamo infatti una valvola avente un otturatore molto ampio ed una pressione nel serbatoio che cresca molto lentamente. Alla pressione di taratura l’otturatore si solleva, ma il deflusso del fluido causa una veloce riduzione della pressione e l’otturatore si abbassa quindi

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 22

4 – La fisica del fenomeno dell’efflusso in una valvola.

verso la sede. La pressione necessaria può essere invece mantenuta solo in presenza di un’adeguata portata. Se ciò non accade si verificheranno oscillazioni della pressione che comporteranno oscillazioni nel grado di apertura dell’otturatore. Quindi gli effetti dinamici si vanno a sovrapporre localmente a livello di variazioni di pressione generando spinte risultanti variabili nel tempo. In questo caso la turbolenza non è la causa dell’instabilità dell’otturatore. Si deve notare inoltre che in funzione del valore dell’alzata variano le forme geometriche delle sezioni a disposizione del fluido (in particolare la legge della variazione delle sezioni). Ciò induce forti variazioni nel regime fluido-termodinamico e nella fenomenologia dell’efflusso. Il percorso del fluido all’interno della valvola dipende dalle geometrie in gioco: dalla forma dell’otturatore (la presenza su quest’ultimo di bordini che deviano il fluido); dalla forma della sede e di tutta la camera interna; dalla forma delle feritoie di uscita (e da quella del convogliatore se presente); da tutto ciò che è nelle immediate vicinanze della valvola ed in particolare nella traiettoria del getto (o dei getti) di uscita. Come è noto, le perdite a livello energetico valutabili in termini di perdita di carico sono imputabili alle irreversibilità naturali dell’efflusso e dipendono dalla viscosità del fluido, dalla sua comprimibilità e dai fenomeni ad essa legati in funzione delle velocità e della pressione. L’interazione con le superfici solide, il percorso da esse stabilito e le perdite per attrito viscoso, portano ad introdurre il concetto di rendimento fluidodinamico della valvola, rendimento valutabile in funzione del coefficiente di efflusso: esso consiste nel rapporto fra portata misurata (al valore di pressione prestabilito) e la portata teorica. 4.2. La valvola vista come un ugello Teoricamente la valvola in oggetto può essere considerata come un ugello. È ovvio che la sua geometria complessa e poco lineare la porti ad essere un ugello “mal disegnato”. Il fluido è costretto a defluire da un ambiente “grande” passando attraverso un restringimento per poi uscire nell’ambiente esterno, previa il passaggio in una zona in cui le sezioni di passaggio aumentano. Quindi in base alle condizioni termodinamiche dell’ambiente esterno ed interno si avranno fenomenologie analoghe a quelle di un ugello del tipo De Laval. Accenniamo ora le basi teoriche dell’efflusso di un gas da un ugello convergente vedi Figura 4.2.1.

Figura 4.2.1

Prendiamo in considerazione un gas in un grande contenitore alla pressione p1, temperatura T1 e velocità v1 trascurabile.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 23

4 – La fisica del fenomeno dell’efflusso in una valvola.

La pressione nell’ambiente di uscita è p2. Dato che p2 è minore di p1, il gas “sente” una differenza di pressione che lo accelera. Il gas si muove quindi dall’ambiente a pressione p1 all’ambiente a pressione p2 attraverso l’ugello. L'energia potenziale del gas sotto forma di pressione si converte gradualmente in energia cinetica, quindi il gas acquista velocità. Aumentando la pressione a monte il gas è soggetto ad un’ulteriore accelerazione ed aumenta la velocità in gola (sezione di area minima). Varia certamente anche la sua densità (la pressione si riduce ed il fluido si espande). Per un certo valore di pressione p1, la differenza in pressione tra p1 e p2 è sufficiente per accelerare il gas alla velocità del suono locale a = kRTt alle condizioni in gola (pt, Tt), un limite cinetico che non può essere oltrepassato, nel senso che un ulteriore piccolo aumento della pressione p1 (energia potenziale) non può aumentare la velocità (energia cinetica) in gola. All’uscita dell’ugello il gas risulta non più “costretto” dalle pareti solide, ed espande irreversibilmente (in condizioni reali) alla pressione circostante p2, decelerando rapidamente. Se a questo punto si abbassa la pressione a valle, le condizioni del gas in gola non cambiano e non varia la portata in massa. Si verifica la stessa situazione se la pressione a valle viene aumentata, finché non si raggiunge nuovamente la pressione di gola: non si ha nessuna influenza sulle condizioni di gola e la portata rimane costante. Questa situazione è nota come flusso in condizioni di “choking” (condizione per cui si raggiunge la velocità limite in gola) in cui, per fissate condizioni a monte, si raggiunge una portata in massa costante indipendentemente dalle condizioni dell’ambiente a valle. Ma cosa accade se la pressione di gola è aumentata ulteriormente? Dato che solamente una frazione dell’energia potenziale (sotto forma di pressione) è necessaria per accelerare il gas alla velocità del suono (con trasformazione in energia cinetica), il surplus di energia si manifesta nell’aumento della pressione in gola. Aumentano così la velocità locale del suono e la densità del gas che si traducono nell’aumento della portata in massa. Si osserva quindi che dalle condizioni di “choking”, variazioni della pressione nell’ambiente di uscita non hanno effetto sulla portata, al contrario delle variazioni della pressione in gola. La relazione matematica per il calcolo della portata in condizioni di choking, anche dette “critiche” è la seguente: k +1

m& = At p1 *

k +1

k  2  k −1  2  k −1   = At ρ1a1   RT1  k + 1   k +1

(4.2.1)

dove: • m& * è la portata critica (corrispondente alle condizioni del fluido a monte); • At è l’area della sezione di gola; p1 , ρ1 , T1 sono rispettivamente la pressione, la densità e la temperatura del fluido nell’ambiente a monte; • a1 = kRT1 è la velocità del suono locale alle condizioni termodinamiche del fluido a monte; • R è la costante del gas; • k è il rapporto delle capacità specifiche di calore a pressione e volume costante del fluido. Da questa equazione si vede come la portata massica critica, fissata la sezione di gola (fisso At), è una funzione esclusivamente della pressione di ingresso e della temperatura e della natura del gas evoluente.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 24

4 – La fisica del fenomeno dell’efflusso in una valvola.

La portata massica critica è realizzata se si è nelle condizioni in cui è rispettata la seguente: k

p2  2  k −1 ≤  p1  k + 1 

(4.2.2)

dove p2 è la pressione di uscita dall’ugello. Da questa equazione si vede che il rapporto di pressione critico (p2/p1) è funzione esclusivamente della natura del gas. Per l’aria (ed in generale per i gas biatomici a condizioni ambiente) quando il rapporto delle pressioni assolute è 0.528, cioè quando la pressione assoluta (p2) a valle è 52.8% della pressione assoluta (p1) a monte. Per flusso di aria attraverso un orifizio con temperatura dell’aria di gola di 20°C la velocità di “bloccaggio” (sonica) risulta circa 345 m/s (circa 1239 km/h). Quando la velocità dell’aria raggiunge la velocità sonica, gli ulteriori aumenti in p1 non causano ulteriore aumento della velocità di gola ma la densità aumenta; e dato che la portata in massa è anche funzione della densità, la portata aumenta linearmente con la pressione a monte. Anche se la velocità dell’aria attraverso l'orifizio è limitata dalla velocità del suono, la portata massica continua ad aumentare con la pressione assoluta p1 (vedi Figura 4.2.2 nella quale si riporta l’andamento qualitativo della portata in massa in funzione delle variazioni della pressione di scarico p2 con, a parametro la pressione a monte p1.

Figura 4.2.2

Nella Figura 4.2.3 invece si mostra l’andamento qualitativo della portata in massa attraverso l’ugello al variare della pressione a monte (con pressione di scarico p2 costante).

Figura 4.2.3

Come abbiamo affermato in precedenza però, la valvola può essere assimilata ad un ugello convergente-divergente. Ciò comporta la possibilità dell’instaurarsi di condizioni per cui Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 25

4 – La fisica del fenomeno dell’efflusso in una valvola.

l’efflusso può raggiungere almeno localmente velocità supersoniche: parleremo di efflusso transonico. Dalla letteratura e dalla sperimentazione si evince che la presenza del tratto divergente permette il raggiungimento di tali velocità. Tutto dipende dalle condizioni fluidotermodinamiche a monte e a valle. Sempre in funzione della differenza di pressione tra monte e valle si sperimenteranno fenomeni quali onde d’urto attraverso le quali il flusso da supersonico viene riportato bruscamente a velocità subsoniche. Gli spessori di questi fronti sono dell’ordine delle decine di micron e quindi tali bruschi salti sono accompagnati ad un repentino aumento di entropia connesso con fenomeni dissipativi ed irreversibili in virtù della viscosità del fluido nonché dal forte gradiente di temperatura che si genera localmente. Si può affermare che la posizione (a valle della sezione di gola) in corrispondenza della quale si produce l’onda d’urto dipende, per ogni valore di p1, dalla pressione di scarico p2, e che in generale essa si sposta verso valle al decrescere di quest’ultima. Si deve notare come la trattazione teorica relativa al classico ugello De Laval può essere ritenuta valida per il caso della valvola in oggetto con la riserva sul reale orientamento (e forma) delle onde d’urto che è molto influenzato dalle sedi degli attriti. Infatti l’interazione fra onda d’urto e strato limite e la geometria “tormentata” delle valvole introducono generazione di onde che si influenzano vicendevolmente. Inoltre al valore minimo delle sezioni di passaggio (sezione alla quale corrisponde la cosiddetta “area critica”) non corrisponde necessariamente la velocità sonica. Infatti in un flusso reale tridimensionale le sezioni minime di efflusso non corrispondono necessariamente con quelle minime geometriche proprio a causa della realtà dell’efflusso ed ai fenomeni di compressibilità. Da tutti questi fattori deriva la difficoltà della previsione di questi fenomeni e dell’evolversi dell’efflusso.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 26

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

5. Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

5.1. Dominio di calcolo: come isolare una porzione del sistema fisico completo. Per chiarezza e praticità applicheremo la seguente trattazione teorica ed applicativa alle geometrie della valvola in studio. Per la realizzazione delle simulazioni fluidodinamiche (calcolo fluido-termodinamico), si applicano metodi matematici di risoluzione numerica delle equazioni di Navier-Stokes nel dominio fluido con determinate condizioni al contorno ed iniziali. Tali metodi offrono la possibilità di discretizzare anche le equazioni differenziali alle derivate parziali, tipiche dei problemi della fluidodinamica. Si ottiene così una “soluzione “ con una approssimazione molto buona. Esistono varie metodologie di risoluzione ma tutte, come base della loro teoria, introducono il concetto di “discretizzazione”. Esso consiste nella rappresentazione (mappatura) del dominio spaziale continuo (quindi contenente infiniti punti) come un dominio discreto avente un limitato e quindi finito numero di punti e/o volumi, definiti usualmente nodi e celle rispettivamente. Il dominio di calcolo è una porzione spaziale, nelle tre dimensioni, del sistema fisico completo. Inizialmente si deve definire tale porzione con criteri legati alla fisica del fenomeno ed al procedimento numerico di calcolo che si prevede di applicare. 5.2. Superfici di contorno: pareti e sezioni del volume fluido completo. Il dominio è limitato da pareti solide (quelle interne ed esterne della valvola e del serbatoio etc.), e da superfici corrispondenti a sezioni del volume fluido intorno alla valvola (nel nostro caso: aria a condizioni ambiente e aria in pressione all’interno del serbatoio). Nel caso in questione, non considerando fenomeni di conduzione termica nel solido (valvola e serbatoio) perché irrilevanti, il dominio non è altro che il volume “fluido” nel quale avviene l’efflusso dell’aria. Quindi la definizione dei limiti spaziali di tale dominio è imposta principalmente dalla geometria e dalla configurazione della valvola nelle condizioni di prova (ad esempio si deve tenere conto dell’alzata effettiva dell’otturatore che, come è stato indicato al paragrafo 2.2 Dati sperimentali, in tutte le prove di qualifica è stata fissata a 5mm) ma, come si spiega oltre, anche dalle sue eventuali simmetrie. Un’assunzione riguardo alla geometria è stata fatta riguardo all’installazione della valvola sul serbatoio di prova. La scelta arbitraria è stata quella di ipotizzare un collegamento diretto della valvola sulla superficie orizzontale del serbatoio, senza interposizione di condotti appositi e/o altri tipi di adattatori. L’altra assunzione si basa sull’aver ipotizzato una superficie superiore esterna del serbatoio piana orizzontale ed estesa fino alle superfici verticali costituenti i confini del dominio di calcolo. Un’altra approssimazione consiste nell’aver semplificato (senza modifiche alle dimensioni generali reali) con superfici semplici (cilindriche e piane) le superfici esterne della valvola in realtà più complesse.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 27

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

Si sono usati raccordi a spigolo vivo per tutte le linee di intersezione tra le varie superfici senza tenere conto dei raggi di raccordo reali e delle tolleranze dimensionali di lavorazione. Tali assunzioni sono state fatte valutando percentualmente poco influente il loro effetto ai fini del calcolo soprattutto per le superfici esterne che il fluido lambisce a velocità molto basse. In Figura 5.2.1 si riportano le sezioni significative con le quote di interesse della valvola in studio.

Figura 5.2.1

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 28

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

In Figura 5.2.2 si riporta il disegno della valvola montata sul serbatoio. Si è assunto il montaggio della valvola direttamente su di un serbatoio con superficie piana. Questa assunzione è stata fatta per semplicità, prevedendone un’influenza trascurabile sulle simulazioni.

Figura 5.2.2

5.3. Individuazione di piani di simmetria. Nel caso della valvola in studio si sono individuati due piani ortogonali di simmetria passanti per l’asse geometrico (vedi Figura 5.3.1 – Piani di simmetria). Essi fissano quattro porzioni spaziali nelle quali si prevedeva verosimilmente che l’efflusso, a meno delle fluttuazioni dovute alla turbolenza, mostrasse le stesse caratteristiche fluido-termodinamiche mediate nel tempo. L’efflusso nelle quattro porzioni si ottiene “specchiando” quello ottenuto per una porzione, rispetto ai piani di simmetria definiti. Il vantaggio dell’utilizzo dei due piani di simmetria è quello di poter utilizzare le risorse di calcolo in una porzione fisica (un quarto della totale) molto ridotta rispetto a quella complessiva, senza perdere efficacia. Infatti i limiti di velocità (CPU) e di capacità (dimensioni della memoria RAM) del calcolatore limitano il numero delle celle da poter generare per ottenere risultati in tempi ragionevoli. Quindi la possibilità di ridurre lo spazio fisico da grigliare permette un miglior utilizzo del numero di celle a disposizione.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 29

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

I nomi delle entità che compongono le superfici sono: simx, simx:064 sul piano X-Z, e simy, simy:066 sul piano YZ (vedi Figura 5.3.2- piano di simmetra X-Z, Figura 5.3.3- piano di simmetra Y-Z, Figura 5.3.4– superfici sui piani di simmetria, superfici esterne della valvola e superfici sul mantello superiore esterno del serbatoio e Figura 5.3.5– particolare delle superfici sui piani di simmetria, superfici esterne della valvola e superfici sul mantello superiore esterno del serbatoio).

Figura 5.3.1 – Piani di simmetria

Figura 5.3.2- piano di simmetra X-Z

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Figura 5.3.3- piano di simmetra Y-Z

Roma, 14/06/2005

Pagina 30

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

simx – simx0:64

Valvola simy – simy0:66

serbout1-2-3-4

Figura 5.3.4– superfici sui piani di simmetria, superfici esterne della valvola e superfici sul mantello superiore esterno del serbatoio

Valvola simx – simx0:64

simy – simy0:66

serbout1-2-3

Figura 5.3.5– particolare delle superfici sui piani di simmetria, superfici esterne della valvola e superfici sul mantello superiore esterno del serbatoio

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 31

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

5.4. Definizione delle superfici di contorno La scelta della posizione, della forma e della distanza delle superfici di contorno “fluide” del dominio dipendono dalla natura attesa dell’efflusso, dalle necessità di costruzione della griglia e della sua semplificazione. Limiti e/o vincoli di tipo numerico (tipo di condizione al contorno da assegnare), geometria, fisica del problema, scelta della tipologia della “griglia di calcolo” da applicare ed esperienza guidano quindi in questa fase. Nel nostro caso le superfici di confine sono state definite nel modo seguente: 5.4.1. Sezioni “fluide”. 1) La superficie verticale parallela alla feritoia, e verso la quale il getto ha direzione predominante, è stata posta a circa cento diametri1 di distanza. Nome dell’entità: pressout3. 2) Le altre superfici verticali al di sopra del serbatoio, sono state poste ad una distanza minima pari a circa 10 diametri. Nomi delle entità: pressout1, pressout2. Allontanandosi dalla valvola, in direzione x la superficie verticale pressout2 si allarga in direzione y per tenere conto del previsto allargamento del getto (vedi Figura 5.4.1 – Sezioni “fluide”). La superficie verticale all’interno del serbatoio ha forma cilindrica ed è posta a circa sette pressoutup1-2-3-4 pressout1 pressoutup5 Pressin1-2-3 pressout3 pressout2

Figura 5.4.1 – Sezioni “fluide”

diametri dalla sezione di entrata nella valvola. Nome dell’entità: pressin3. 3) Le superfici orizzontali parallele al piano del serbatoio e sopra ad esso, sono state poste ad una distanza minima pari a circa sette diametri dal serbatoio stesso. La quota corrisponde a quella massima della valvola. Nomi delle entità che compongono la 1

E’ conveniente, specie in studi “prototipali” di questo tipo, adimensionalizzare le geometria rispetto ad una “lunghezza caratteristica” del problema. In questo caso la fenomenologia nota porta ad assumere il diametro D dell’orifizio di ingresso della valvola come lunghezza caratteristica. Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 32

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

superficie: pressoutup1, pressoutup2, pressoutup3, pressoutup4. Allontanandosi dalla valvola, in direzione x (vedi Figura 5.4.1 – Sezioni “fluide” e Figura 5.4.2 – Sviluppo del dominio di calcolo nella direzione x) la superficie superiore ha una pendenza positiva sempre per tenere conto dell’allargamento del getto. Nome dell’entità: pressoutup5. La superficie orizzontale sulla quale giace la sezione fluida all’interno del serbatoio è stata posta a circa sette diametri dalla sezione di entrata nella valvola. Nomi delle entità che compongono la superficie: pressin1, pressin2. Figura 5.4.2 – Sviluppo del dominio di calcolo nella direzione x

Figura 5.4.3 – Particolare del dominio di calcolo intorno alla valvola

Figura 5.4.4 – Particolare del dominio di calcolo in prossimità dell’otturatore e della feritoia di uscita

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 33

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

Per tutte le superfici le distanze sopra dette sono da interpretarsi come minime accettabili. Infatti si deve evitare la coincidenza di tali superfici con zone in cui il flusso deve ancora svilupparsi, zone di ricircolo, di riattacco, in generale zone in cui si possa prevedere l’esistenza di alti gradienti delle grandezze fisiche. 5.4.2. Pareti solide. 4) Parete interna del serbatoio. Nome dell’entità: serbin. 5) Parete esterna del serbatoio. Nomi delle entità che compongono la superficie: serbout1, serbout2, serbout3, serbout4. 6) Condotto interno della valvola e sede conica. Nomi delle entità componenti: condotto, sede1in, sede2in, sede1out, sede2out, sedevert1, sedevert2. 7) Superficie dell’otturatore (bordino compreso). Nomi delle entità componenti: orizzout5, orizzout6, vertwall6, vertwall7, inclinup1, inclinup2, bordup1, bordup2, bordvertup1, bordvertup2, upwall1, upwall2, upwall3, upwall4. 8) Superfici orizzontali all’interno della valvola. Nomi delle entità: upwall5, upwall6, orizzout1, orizzout2, orizzout3, orizzout4. 9) Superfici verticali all’interno della valvola. Nomi delle entità: vertwall1, vertwall2, vertwall3, vertwall4, vertwall5, vertwall8, vertwall9. 10) Superfici esterne della valvola. Nomi delle entità: ext1, ext2, ext3, ext4, ext5, ext6, extup. Individuato il dominio di calcolo si è effettuata la generazione della griglia (o mesh) di calcolo suddividendolo in un elevato numero di volumi e sottovolumi (discretizzazione spaziale). 5.5. Generazione della griglia di calcolo. Nella prima fase si procede quindi alla generazione della geometria 3D di tipo “wireframe” con l’ausilio di un software CAD (o direttamente con programmi per la modellazione 3D). In questa fase si delimita quindi il dominio di calcolo. Il disegno della geometria si può limitare alla generazione dei punti, estremi dei segmenti o circonferenze o archi che la formano. In questa fase le entità sono disegnate in funzione del tipo di griglia che si vuole creare e dagli eventuali sottovolumi in cui la si prevede di dividere. Per la creazione di quest’ultimi si è fatto uso del software commerciale GAMBIT® release 1.3.0 dopo avervi importato le entità geometriche. GAMBIT® è un pre-processore integrato che permette, tramite interfaccia grafica, di costruire la geometria e di generare e controllare la qualità della griglia di modelli per simulazioni termofluidodinamiche ed altre applicazioni scientifiche. Altra sua caratteristica è la possibilità dell’assegnazione del tipo di “boundary” alle diverse entità geometriche e la completa compatibilità con il codice di calcolo usato per le simulazioni: FLUENT 5.4.8. 5.6. Tipologia di griglia impiegata. Nel caso in esame è stata generata una griglia multiblocco con griglie strutturate e non strutturate. Si è fatto anche uso di griglie dette “non conformi”, (vedi Figura 5.6.11 - Griglia) caratterizzate dalla mancanza di coincidenza tra i nodi giacenti su due superfici affacciate e Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 34

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

coincidenti, ma appartenenti a sottodomini separati. FLUENT elabora con interpolazione i flussi delle grandezze attraverso tali “superfici di interfaccia”. Per tutta la griglia sono stati impiegati elementi piani quadrilateri a quattro nodi ed esaedrici ad otto nodi. GAMBIT® permette un controllo avanzato delle celle vicino alle pareti mediante l’attivazione dei “boundary layers” (vedi Figura 5.6.11 - Griglia). Si possono controllare gli strati di celle vicino alle pareti impostando il numero e la legge di variazione dello spessore di quest’ultime. I “boundary layers” possono servire anche per effettuare transizioni all’interno della griglia per ridurre uniformemente il numero dei nodi (e delle celle) allontanandosi dal primo strato di controllo. Tale primo strato può anche non essere adiacente ad una parete (vedi Figura 5.6.11 - Griglia). La qualità e l’efficienza della soluzione numerica sono fortemente dipendenti dalla griglia usata nel modello di simulazione. La posizione delle celle della griglia esercita un’influenza sostanziale sulla stabilità e la convergenza della soluzione numerica: per esempio, se le celle non sono sufficientemente concentrate nelle zone con alti gradienti delle grandezze fisiche del flusso, il solutore non può “catturare” tali gradienti. Quindi nelle regioni dove si prevedono elevati gradienti, è necessario che la griglia sia abbastanza fine da minimizzare la variazione delle variabili del flusso da cella a cella. Nella costruzione della griglia occorre raggiungere un compromesso tra la densità dei nodi, la velocità di calcolo e l’accuratezza della soluzione. Generalmente più punti vengono assegnati in una griglia e più accurata sarà la soluzione, dal momento che diminuisce la spaziatura della griglia nella discretizzazione del dominio di calcolo, ma anche più lento sarà il processo di calcolo poiché aumentano le celle nelle quali vengono calcolate le variabili del flusso. Quindi è necessario concentrare i nodi nelle zone con elevati gradienti cercando di ottenere una griglia sufficientemente “raffinata” per la precisione della soluzione ma anche abbastanza “larga” per la velocità di calcolo. Esistono due tipi di griglie: la griglia strutturata (vedi Figura 5.6.1 - Griglia strutturata) nella quale i nodi sono posizionati all’intersezione di due o tre famiglie di linee (che formano un sistema di coordinate curvilinee) e la griglia non strutturata (vedi Figura 5.6.2 - Griglia non strutturata) nella quale i nodi non possono essere identificati con linee coordinate.

Figura 5.6.1 - Griglia strutturata

Figura 5.6.2 - Griglia non strutturata

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 35

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

Nelle griglie strutturate la trasformazione da spazio fisico a spazio di calcolo è introdotta per “mappare” un sistema di coordinate non rettangolare nello spazio fisico in un sistema rettangolare nello spazio di calcolo. Le griglie strutturate permettono di controllare in modo più preciso la distribuzione dei nodi e la forma delle celle. Nel caso dei volumi finiti due caratteristiche vengono considerate nel determinare la forma delle celle: la “skewness” che è una misura della differenza tra la forma della cella e la forma di una cella equilatera di pari volume (una griglia quadrilatera ottimale ha tutte le celle con angoli prossimi a 90°) e lo “aspect ratio” che è una misura dell’allungamento delle celle; una variazione rapida nel volume di celle adiacenti può portare ad errori di troncamento elevati. Questi due parametri, influenzando la soluzione, indicano l’effetto che ha la fattura della mesh sul modello di discretizzazione. Nel caso di elementi che presentano un alto fattore di skewness (Figura 5.6.3 - Esempi di non ortogonalità-a), il calcolo del gradiente tra la faccia di separazione e la congiungente i centroidi a destra e sinistra della faccia stessa restituisce un valore che non è la derivata rispetto alla normale della superficie, bensì rispetto alla sopraccitata congiungente. In questo caso è necessario modificare le equazioni discretizzate introducendo dei termini correttivi che in generale rallentano il calcolo. Un secondo tipo di non ortogonalità è rappresentato in Figura 5.6.3 Esempi di non ortogonalità-b (non coincidenza tra il centro della faccia e l’intersezione della stessa con la congiungente i centroidi)

Figura 5.6.3 - Esempi di non ortogonalità

Nel caso in cui le celle siano particolarmente allungate in una direzione (l’aspect ratio è alto) può accadere che la massa che attraversa l’elemento non sia valutata correttamente e il flusso che attraversa il lato minore sia valutato in maniera sbagliata. In questo caso la soluzione si dice disaccoppiata e generalmente sono necessari un numero maggiore di iterazioni per raggiungere la soluzione. Si deve comunque valutare il peso relativo dei due parametri posti a confronto. Riportiamo un esempio riguardo alla griglia (non ottimizzata) in prossimità del bordino dell’otturatore e della sede sul piano di simmetria x-z. Nelle Figura 5.6.4 – equiangle skew-bordino, Figura 5.6.5 – aspect ratio-bordino, Figura 5.6.6 – equiangle skew-sede, Figura 5.6.7 – aspect ratio-sede sono evidenziate con i colori della scala del rosso le celle con aspect ratio o equiangle skewness elevati. Si osserva che mentre lo aspect ratio assume valori bassi in tutta la zona, la mappa corrispondente per il secondo parametro (equiangle skew) evidenzia elementi da tenere sotto controllo.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 36

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

Figura 5.6.4 – equiangle skew-bordino

Figura 5.6.6 – equiangle skew-sede

Figura 5.6.5 – aspect ratio-bordino

Figura 5.6.7 – aspect ratio-sede

L’utilizzo delle griglie strutturate presenta maggiori difficoltà nel generare la mesh poiché il dominio di calcolo deve essere diviso in sottodomini per avere una forma che sia “trattabile” e la distribuzione dei nodi sulle linee di confine dei sottovolumi deve essere tale da avere una corrispondenza ‘nodo-nodo’ in ogni sottodominio (a meno dell’uso di mesh non conformi). Facendo riferimento alla Figura 5.6.8- Griglia strutturata, i segmenti 1 e 2 devono avere lo stesso numero di nodi, così come i segmenti 3 e 4.

Figura 5.6.8- Griglia strutturata

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 37

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

Le griglie bidimensionali non strutturate possono essere generate attraverso elementi di due forme, triangolare e quadrata. Le griglie tridimensionali non strutturate possono essere generate attraverso tetraedri, piramidi, prismi a base triangolare e quadrata (esaedri) in funzione della tipologia degli elementi bidimensionali. Le griglie costruite con elementi triangolari sono più facili da generare intorno a geometrie complesse, ma presentano problemi per la velocità di calcolo. Infatti per una data distribuzione di nodi il numero di elementi triangolari è maggiore di quello dei corrispondenti elementi quadrangolari, quindi l’utilizzo di elementi triangolari è penalizzante per la velocità di calcolo. Le griglie non strutturate permettono di avere una diradazione dei nodi nelle zone di minor interesse, non essendo soggette ad una corrispondenza uno a uno dei nodi, ma permettono un minore controllo nella generazione della griglia. In questa simulazione sono state utilizzate entrambe le tipologie di griglia. Ad esempio all’interno della valvola, dove i gradienti sono elevati, sono stati creati sottovolumi dotati di griglie strutturate controllate negli strati adiacenti alle pareti. In alcuni di essi (vedi quelli adiacenti contemporaneamente alla sede ed al bordino circonferenziale dell’otturatore) le griglie di alcune facce sono non strutturate perché di più semplice generazione e perché permettono infittimenti localizzati grazie al posizionamento opportuno dei nodi sui segmenti costituenti i lati delle facce stesse; senza i vincoli delle strutturate. Nella Figura 5.6.9 – Particolare della griglia all’interno della valvola, si osservano gli infittimenti proprio vicino alla sede e al bordino. Il sottodominio in questione è stato grigliato adottando il “Cooper scheme”. Vedi il seguente paragrafo 5.7 per un maggiore approfondimento su tale schema. Altre griglie non strutturate sono state usate nel caso di facce triangolari (vedi Figura 5.6.10 – Particolare della griglia all’interno della valvola, in prossimità del bordino dell’otturatore) sulla superficie dell’otturatore) o laddove non è richiesto un gran controllo delle celle di calcolo (vedi Figura 5.6.9 – Particolare della griglia all’interno della valvola e Figura 5.6.11 - Griglia). Anche in questi casi si è adottato il “Cooper scheme”. Ovviamente nulla vieta un approccio diverso alla generazione della griglia. Nelle zone con gradienti molto elevati, per esempio al bordo della sede dove il flusso da assiale diventa radiale con elevati gradienti dei parametri di flusso, la griglia strutturata avrebbe consentito celle più regolari sia nella forma che nella distribuzione, ma sarebbe stato necessario introdurre un maggior numero di sottovolumi.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 38

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

Figura 5.6.9 – Particolare della griglia all’interno della valvola

Figura 5.6.10 – Particolare della griglia all’interno della valvola, in prossimità del bordino dell’otturatore

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 39

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

Esempio di griglia “non strutturata” (Cooper scheme)

Esempio di griglia “strutturata”

Esempio di utilizzo dei “boundary layers” con elementi di transizione

Griglie “non conformi” con griglie d’interfaccia

Esempi di griglia “non strutturata” (Cooper scheme)

Figura 5.6.11 - Griglia

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 40

5 - Realizzazione delle meshes “caratteristiche”: la discretizzazione del continuo.

5.7. Brevi cenni descrittivi sul “Cooper scheme”

In alternativa alla creazione delle meshes, a partire dal volume è possibile generare la griglia in un volume utilizzando una faccia come sorgente. Nella letteratura tecnica questo schema è chiamato “di Cooper” (Cooper scheme); la mesh di una delle facce viene presa come modello di struttura da ripetere nel volume all’interno del quale essa viene estrusa. Quando si applica questo schema, l’intero volume è trattato come se fosse costituito da uno o più cilindri logici, quindi con le due facce superiore ed inferiore e la superficie laterale. La faccia superiore ed inferiore sono le facce sorgente, mentre tutte le altre facce del volume generico non utilizzate come sorgente sono viste come la superficie laterale del cilindro (Figura 5.7.1). Il metodo per quanto sia estremamente potente è difficilmente generalizzabile in applicazioni complesse, per cui spesso necessitano trattamenti correttivi che includono anche la creazioni di sottodomini.

Figura 5.7.1

La struttura della faccia sorgente oltre ad essere estrusa, può anche essere ruotata di 360° intorno ad un asse; anche così si ottiene una mesh in 3D data dalla ripetizione di una struttura base. La creazione di una mesh non è un processo irreversibile: in qualsiasi momento deve essere possibile cancellare una (parte della) mesh. Esistono però delle gerarchie: non è possibile eliminare una mesh se essa serve per la creazione di una mesh di ordine superiore. Se si è meshato il cubo, non è possibile ad esempio eliminare la mesh di una delle facce, mentre è possibile eliminare la mesh del cubo e lasciare inalterate quelle delle facce, che solo a questo punto possono essere cancellate.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 41

6 – Le simulazioni numeriche. Generalità.

6. Le simulazioni numeriche. Generalità. Il software utilizzato per lo studio effettuato è FLUENT 5.4.8. Il metodo di calcolo utilizzato è il metodo ai volumi finiti che consiste in: 1. Divisione del dominio di calcolo in volumi di controllo discreti utilizzando una griglia di calcolo. 2. Integrazione delle equazioni di governo del flusso su ogni volume di controllo per determinare le equazioni algebriche per le variabili incognite (velocità, pressione, temperatura, etc.); l’integrazione porta ad equazioni discrete che comportano comunque la conservazione di ogni grandezza sul singolo volume di controllo. 3. Linearizzazione delle equazioni discretizzate e soluzione del sistema di equazioni per produrre valori aggiornati delle variabili. Per la chiusura del sistema delle equazioni di Navier-Stokes (che esprimono l’equazione di bilancio della quantità di moto nella quale si è introdotta l’equazione costitutiva del fluido), si affiancano l’equazione della continuità, l’equazione dell’energia e l’equazione di stato per il fluido sotto l’ipotesi di continuità di quest’ultimo; questa assunzione implica che esistono le derivate di tutte le variabili dipendenti. In altre parole, proprietà locali come la densità e la velocità sono definite come medie su elementi “grandi” se comparati con la struttura microscopica del fluido, ma abbastanza “piccoli” in confronto alla scala dei fenomeni macroscopici. Ciò permette di descriverli con l’uso del calcolo differenziale. Per tenere conto della natura reale e turbolenta dell’efflusso il fluido è considerato viscoso e si introduce il modello k-ε RNG per simulare le interazioni turbolente. Un moto turbolento non è mai realmente stazionario e quando si parla di moto turbolento stazionario si intende riferirsi alla invariabilità nel tempo della velocità mediata su di un tempo sufficientemente lungo. Si è optato per tale modello per i seguenti motivi: § È uno dei più utilizzati in campo industriale; § I suoi punti di forza, come pure i suoi limiti, sono ben documentati; § La variante “realizable k-ε” del modello di base “Standard k-ε” fornisce prestazioni superiori per: 1. flussi in presenza di getti piani o circolari; 2. strati di celle sulle superfici di confine in cui esistono alti gradienti positivi (nel senso del moto) della pressione e separazione; 3. rotazione e ricircoli; 4. linee di flusso ad alta curvatura § Fornisce un’accuratezza accettabile per una vasta tipologia di flussi turbolenti in applicazioni industriali che comportano anche scambi di calore; § Siamo in possesso di un’ampia casistica di simulazioni basate sull’applicazione di tale modello ad efflussi caratterizzati dal numero di Reynolds e da effetti dovuti alla compressibilità anche molto diversi. Tutte queste ragioni sono state infine confermate dai risultati che saranno mostrati nei paragrafi successivi.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 42

6 – Le simulazioni numeriche. Generalità. 6.1. Impostazioni del calcolo. Riportiamo nel seguente elenco le caratteristiche e le impostazioni del calcolo effettuato: Simulazioni in stato stazionario (il campo di moto e termico sono indipendenti dal tempo). Le proprietà caratteristiche del mezzo fluido sono considerate quindi come funzioni dello spazio nel sistema di riferimento. q Geometria fissa 3D (tutte le entità geometriche non cambiano di posizione e/o forma nel tempo). q Griglie di base: griglia multiblocco con uso di griglie strutturate e non strutturate e “non conformal grids” con griglie di interfaccia. q Griglie ottimizzate nei primi due o tre strati vicino alle pareti e/o in funzione della y+ (“solution-adaptive refinement”: l’ottimizzazione delle griglie è effettuata sulla griglia base in funzione dei risultati ottenuti dal calcolo preliminare sulla griglia base stessa)1. q 170.000 ÷ 420000 celle di calcolo (le griglie con maggior numero di elementi sono quelle “ottimizzate”). q Solutore: accoppiato, esplicito con discretizzazione “first order upwind”. Significato: § “accoppiato”: tutte le equazioni sono risolte simultaneamente (cioè “accoppiate”); § “esplicito”: metodo di linearizzazione delle equazioni di governo. Per una data variabile, il valore incognito in ogni cella è calcolato usando una relazione che include solo i valori esistenti. Inoltre ogni incognita appare in una sola equazione nel sistema e le equazioni per i valori incogniti in ogni cella possono essere risolti uno alla volta per fornire le quantità non note. § “first order upwind”: discretizzazione del primo ordine in avanti. FLUENT memorizza i valori discreti delle grandezze scalari nei centri cella. Comunque, i valori di tali grandezze sulle facce delle celle servono per i termini convettivi delle equazioni di governo discretizzate e devono essere interpolati dai valori nei centri cella. Ciò si può fare con uno schema “upwind” nel quale tali valori sono derivati da quelli delle celle adiacenti a monte (nel senso dell’efflusso), da cui il nome dello schema. Il solutore impiegato è adatto alla risoluzione di flussi compressibili e ad alta velocità. Dato che le equazioni di governo sono non lineari ed accoppiate, parecchi cicli iterativi devono essere completati prima di raggiungere la convergenza. Ogni iterazione consiste nei seguenti passaggi: 1. Le proprietà del fluido sono aggiornate basandosi sulla soluzione corrente (per il primo ciclo si parte dalla soluzione di inizializzazione di cui parleremo in seguito). 2. Le equazioni di Navier-Stokes e dell’energia sono risolte simultaneamente. 3. Le equazioni per le grandezze scalari, quali quelle per la turbolenza, sono risolte usando i valori aggiornati delle altre variabili. 4. Controllo della convergenza per il sistema di equazioni. 5. Il ciclo viene ripetuto sino al raggiungimento della “convergenza”, che viene valutata sia in base ai valori assunti dai “residui” sia in base alla storia di convergenza. Anche l’analisi del bilancio della massa entrante e uscente con il procedere del procedimento iterativo (si è valutata la sua stabilizzazione) ha permesso di stabilire la convergenza. Anche la stabilizzazione dei valori q

1

Per il significato di y+ vedi i paragrafi 6.2 e 6.3

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 43

6 – Le simulazioni numeriche. Generalità. assunti dalle grandezze fisiche in zone diverse e significative rappresenta un criterio per la valutazione della convergenza e della validità del calcolo stesso. q Fluido: aria. Si tiene conto della sua compressibilità e dei fenomeni ad essa collegati. L’aria è trattata come un gas ideale con viscosità dinamica µ pari a µ = 1.7894*10-5Kg/m*s. Quindi la densità ρ viene calcolata ad ogni iterazione con l’equazione dei gas ideali: ( p + p) ρ= o RT dove po è la pressione ambiente, p è la pressione statica locale relativa alla pressione ambiente, R che vale 287 J/(kg*K), è la costante del gas (aria nel nostro caso) e T è la temperatura locale del fluido espressa in kelvin. q Pressione operativa: Fluent permette di stabilire una pressione detta “operativa” per evitare errori di troncamento in simulazioni a basso numero di Mach. Il codice sottrae la pressione operativa dalla pressione assoluta e usa nei calcoli la pressione relativa. Nel caso in oggetto (anche se si raggiungono valori del numero di Mach pari 2.5) la pressione operativa è stata fatta coincidere con la pressione ambiente. La relazione che lega pressione assoluta pabs, pressione operativa (ambiente) po e pressione relativa prel è: p abs = p o + p rel Tutte le pressioni inserite e calcolate nelle simulazioni in Fluent sono in termini di pressione relativa prel, per cui, d’ora in avanti, il suffisso ‘rel’ sarà omesso. q Modello di turbolenza: k-ε RNG (a due equazioni) con opzione attivata per “swirl dominated flow”. Le RANS (Reynolds-averaged Navier-Stokes), equazioni di NavierStokes mediate nel tempo, rappresentano le equazioni di trasporto per le sole grandezze medie del flusso. Tutte le scale della turbolenza devono essere modellate e si devono introdurre termini aggiuntivi nelle equazioni di governo. Per la “chiusura” del sistema si introducono equazioni aggiuntive. L’approccio Reynolds-averaged è generalmente adottato per problemi ingegneristici. Il modello k-ε introduce due equazioni di trasporto insieme alla “energia cinetica turbolenta” ed al “tasso di dissipazione dell’energia cinetica turbolenta”. La soluzione delle due equazioni aggiuntive permette il calcolo della viscosità turbolenta per le equazioni RANS. Esso è un modello semi empirico; la sua robustezza, l’economia in termini di risorse di calcolo ed una ragionevole accuratezza per una vasta casistica di flussi turbolenti rende conto del suo ampio uso in molteplici applicazioni. Il modello k-ε RNG introduce termini aggiuntivi nelle sue equazioni che migliorano l’accuratezza per flussi con alti effetti di “swirl” ed altamente ”stirati”. q Funzioni di parete: standard. Funzioni che permettono la risoluzione del campo fluidotermodinamico in prossimità delle pareti. Affinchè tali funzioni siano ben applicate si devono controllare gli spessori degli strati di celle adiacenti alle pareti. Esistono prescrizioni diverse in base al tipo di funzioni di parete applicate. Il controllo si esegue in funzione dei valori assunti dalla grandezza y+. 6.2. Cenni sulla turbolenza. Nel campo di applicazione delle valvole in oggetto, i flussi dei fluidi sono quasi sempre turbolenti.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 44

6 – Le simulazioni numeriche. Generalità. Questo significa che il moto del fluido è stocastico, non stazionario e tridimensionale. Altre caratteristiche dei flussi turbolenti sono le seguenti: • la diffusività che causa rapidi mescolamenti e alti tassi di trasporto della quantità di moto, di calore e di massa; • la dissipazione: gli sforzi tangenziali viscosi producono lavoro di deformazione con relativi incrementi di energia interna del fluido a spese della energia cinetica della turbolenza. Questa ha bisogno di un continuo apporto di energia per compensare le perdite viscose: se non c’è apporto di energia la turbolenza decade rapidamente. Per queste ragioni, il moto turbolento ed i fenomeni di trasferimento di calore e di massa con esso associati sono estremamente difficili da descrivere e da predire in maniera teorica. Se si considera il moto di un fluido in un condotto rettilineo o a bassa curvatura, si osserva che, finché la velocità si mantiene sufficientemente bassa, ogni particella fluida si muove con velocità uniforme descrivendo una traiettoria parallela all’asse del condotto. Le forze viscose rallentano le particelle vicino alle pareti. Nel contempo le particelle a diverse altezze, anche se con velocità diverse, scorrono “ordinatamente” in strati contigui. In questo caso si è in condizioni di regime laminare. All’aumentare del numero di Reynolds2 si osserva che, in corrispondenza ad un certo Re critico, scompare improvvisamente questo comportamento ordinato, caratteristico del moto laminare. Si verifica cioè la (brusca) transizione al regime turbolento, nel quale le traiettorie descritte dalle diverse particelle non sono più rettilinee e la velocità in un dato punto è diversa da istante a istante. Al moto medio nella direzione dell’asse del condotto si sovrappone un moto fluttuante in direzione trasversale e longitudinale, che genera un trasferimento di quantità di moto orizzontale da uno strato all’altro e dà luogo ad un mescolamento molto più energico. In conseguenza di ciò, nel moto turbolento si ha una maggiore uniformità trasversale della velocità (Figura 6.2.1). Nel caso di trasporto turbolento quindi, le dimensioni delle particelle e la distanza che esse percorrono prima di “cedere” le loro proprietà (quantità di moto, energia) dipendono oltre che dalla natura del fluido e dal suo stato termodinamico anche dal tipo di moto.

Figura 6.2.1

La dimensione degli agglomerati fluidi, che si formano e si disintegrano continuamente e che per brevi periodi si muovono quasi come un unico corpo, determina la scala della turbolenza e per uno stesso fluido può essere sostanzialmente diversa a seconda della velocità, delle condizioni al contorno (pareti, ecc.) e della “storia” del fluido. Le proprietà di trasporto turbolento non sono quindi una caratteristica del fluido ma variano da un tipo di moto ad un altro. In un moto turbolento il trasporto macroscopico predomina nettamente sul trasporto microscopico nel quale le proprietà di trasporto (µ,k) dipendono dalle caratteristiche delle molecole e dal loro libero cammino medio. Quest’ultima considerazione implica che ai fini degli

2

Definizione del numero di Reynolds:

Re =

ρLu ; dove r è la densità, L è una lunghezza caratteristica, u indica la µ

velocità media e µ è la viscosità dinamica. Tale parametro adimensionale quantifica il rapporto tra forze viscose e forze d’inerzia.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 45

6 – Le simulazioni numeriche. Generalità. scambi energetici, tutto avviene come se la viscosità o la conduttività termica fossero enormemente maggiori (anche più di 10000 volte). Pertanto nel moto turbolento gli sforzi resistenti e la conduzione termica sono molto maggiori che nel moto laminare. Il fenomeno della turbolenza comporta quindi il continuo rimescolamento della massa fluida, con cessione di energia da parte di particelle di fluido animate da velocità medie più elevate a quelle più lente, nonché il verificarsi di fenomeni di trasporto e dissipativi di rilevante entità. Nel caso dei flussi ad elevato numero di Reynolds le forze viscose sono piccole rispetto a quelle d’inerzia. Questo caso è frequente in quanto i due principali fluidi che si incontrano nella pratica, l’aria e l’acqua, hanno una viscosità molto piccola e quindi il caso di Re molto elevati si verifica anche per velocità moderate. Ciò non significa che, nelle equazioni di Navier-Stokes, si possano trascurare semplicemente i termini viscosi, e ciò vale in particolare nel caso del flusso in prossimità di una parete solida. Fisicamente ciò significa piuttosto che non è possibile trascurare le forze viscose vicino alla parete poiché, essendo proprio queste forze la causa del “frenamento” della corrente esterna, esse devono essere, in questa zona, dello stesso ordine di grandezza delle forze d’inerzia. Poiché le forze viscose sono del tipo ∂u τ =µ ∂y (relazione di Newton anche attribuita a Petroff e Maxwell), esse possono essere grandi, pur essendo la viscosità molto piccola, se è molto grande il gradiente di velocità normale alla parete, cioè se è molto piccolo lo spessore nel quale la velocità passa dal valore nullo alla parete al valore finito della velocità della corrente. Esiste cioè, in prossimità della parete, una zona detta strato limite (vedi il paragrafo seguente 6.3), che è tanto più sottile quanto più grande è il numero di Reynolds e nella quale non possono essere trascurate le forze viscose. Però quando il numero di Reynolds aumenta, anche il fluido interno allo strato limite mostra una notevole transizione dal regime laminare a quello turbolento. L’intero campo del fluido intorno ad un corpo immerso in una corrente, e in particolare, le forze esercitate su esso, sono fortemente dipendenti dal fatto che il fluido nello strato limite sia laminare o turbolento. La transizione in uno strato limite su un corpo solido in una corrente dipende da molti parametri, i più importanti dei quali sono la distribuzione della pressione nel flusso esterno, la natura della parete (scabrosità) e la natura dei disturbi nel fluido libero (intensità della turbolenza). Un fenomeno che spesso si verifica consiste nel distacco dello strato limite. In questo caso si verificano delle circostanze (dipendenti molto spesso dalle forti curvature dei profili delle pareti) per le quali il fluido, che è stato rallentato nello strato limite, viene trasportato man mano nella corrente esterna. Ciò accade più facilmente nel caso in cui la pressione lungo la parete cresca nella direzione del moto. In queste condizioni il fluido, che per effetto dell’attrito nello strato limite ha perso energia cinetica, non riesce a vincere le forze di pressione che si oppongono al moto ed a penetrare nella zona a pressione crescente. Il fluido in vicinanza della parete tende quindi ad arrestarsi mentre nella zona a valle, per effetto del gradiente di pressione, si genera un flusso che si muove in direzione opposta a quella della corrente principale e che costringe lo strato limite a staccarsi dalla parete (Figura 6.2.2). Il punto di separazione S rappresenta il confine fra la zona di flusso diretto e quella di flusso invertito ed è definito dalla condizione  ∂u    = 0  ∂y  y = 0 come si può osservare da un esame dei profili di velocità nelle due zone (Figura 6.2.3).

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 46

6 – Le simulazioni numeriche. Generalità.

Figura 6.2.2

Figura 6.2.3

La separazione dello strato limite, che avviene sempre in maniera più o meno marcata nella zona posteriore di un corpo immerso in una corrente fluida, dà luogo alla cosiddetta scia che altera sensibilmente la struttura del flusso esterno e dà origine a perdite energetiche aggiuntive (Figura 6.2.4).

Figura 6.2.4

6.3. Effetto delle pareti e significato del parametro dimensionale y+. La maggior parte degli studi sono effettuati su flussi entro canali con pareti solide. La presenza di tali superfici comporta una considerevole differenza nel comportamento del flusso e nella formazione di turbolenza se paragonati a un flusso turbolento libero. Il flusso vicino alle pareti è divisibile idealmente in tre strati: • il sottostrato laminare (inner layer) nel quale sono preponderanti gli sforzi viscosi e il profilo di velocità è lineare; • lo strato esterno (outer layer) nel quale sono preponderanti gli sforzi turbolenti; • lo strato di “sovrapposizione” (overlap layer o “log law” region) nel quale il profilo della velocità media mostra un andamento logaritmico (da cui il nome di “log law” cioè legge logaritmica). L’effetto di una parete viene spiegato definendo un numero di Reynolds locale,

Re =

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

ρUy µ

(6.3.1)

Roma, 16/06/2005

Pagina 47

6 – Le simulazioni numeriche. Generalità. dove U è la velocità, µ è la viscosità dinamica del fluido, ρ è la sua densità e y rappresenta la distanza normale dalla parete. Le forze di inerzia prevalgono nel flusso lontano dalle pareti. Se la distanza diminuisce, anche il numero di Reynolds diminuisce, e, poco prima che la distanza y si annulli, c’è un insieme di valori di y per il quale il valore del numero di Reynolds locale è nell’ordine dell’unità. A questa distanza dalle pareti, e più vicino, le forze viscose sono dello stesso ordine di grandezza, o maggiori, delle forze di inerzia. Le differenti regioni nel flusso vicino alle pareti sono definite sulla base di un parametro adimensionale y+

y+ =

τw ρ

yp v

(6.3.2)

dove yp è la distanza dalla parete del punto più vicino nel quale viene calcolata la velocità, (nel caso dei “volumi finiti” questo coincide con il centro cella ),ν è la viscosità cinematica del fluido e τw è lo sforzo di taglio alla parete. τw = µ

∂U ∂y

(6.3.3) y =0

y+ è simile ad un numero di Reynolds locale, quindi il suo valore determina l’importanza relativa del contributo viscoso e del contributo turbolento nello sforzo di taglio. Per valori di y+ minori di 50 esiste un forte contributo della viscosità molecolare allo sforzo di taglio, mentre per valori di y+ maggiori di 30450 questo contributo diventa trascurabile. La distanza dalla parete delle celle ad essa adiacenti deve essere determinata considerando il campo oltre il quale è valida la “log-law”. Questa è valida per y+>30÷60. Un valore della y+ vicino al limite inferiore (y+ circa uguale a 30) è consigliato (da manuale di Fluent). Un altro parametro dimensionale utile nella descrizione delle regioni di un flusso entro canali con pareti è u+

u+ =

u uτ

(6.3.4)

dove uτ è la ‘friction velocity’ definita dalla relazione

uτ =

τw ρ

(6.3.5)

Lo sforzo di taglio τ(y) è la somma dello sforzo viscoso e dello sforzo di Reynolds − ρ uv . Alla parete, la condizione è che la velocità sia nulla, quindi lo sforzo di Reynolds si annulla, e lo sforzo alla parete τw è dovuto interamente al contributo viscoso. In Figura 6.3.1 Profili degli sforzi viscosi e degli sforzi di Reynolds sono presentati due grafici che mostrano i profili degli sforzi viscosi e degli sforzi di Reynolds in un canale di raggio δ per due tipi di flussi, entrambi turbolenti ma con numero di Reynolds pari a 5600 (linea tratteggiata) e 13750 (linea continua).

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 48

6 – Le simulazioni numeriche. Generalità.

Re=13750

Re=13750

Re=5600

Re=5600

Figura 6.3.1 - Profili degli sforzi viscosi e degli sforzi di Reynolds

6.3.1. Il sottostrato laminare. Nella zona immediatamente adiacente ad una parete, gli sforzi sono totalmente dovuti agli sforzi viscosi. Questo strato è chiamato sottostrato laminare, o lineare, ed è molto sottile, con un basso valore di y+ (y+<5). Lo sforzo di taglio è praticamente costante e uguale allo sforzo di taglio alla parete τw. In questa zona esiste una relazione lineare tra i due parametri dimensionali u+ e y+

u+ = y+

(6.3.6)

Nella Figura 6.3.2 - Profilo della velocità nel sottostrato laminare, la linea continua rappresenta il profilo del parametro u+ ottenuto con una simulazione numerica diretta (DNS). Lo scostamento per valori di y+ inferiori a 5 della relazione (6.2.6) è trascurabile.

Figura 6.3.2 - Profilo della velocità nel sottostrato laminare

6.3.2. Lo strato di ‘sovrapposizione’. Nella regione “log law region” in cui il valore di y+ è maggiore di trenta (y+>30) gli effetti della viscosità molecolare sono trascurabili. La relazione tra u+ e y+ è logaritmica, da cui il nome “log-law” velocity profile

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 49

6 – Le simulazioni numeriche. Generalità.

u+ =

1 ln y + + B k

(6.3.7)

dove k=0.41, B=5.2, il valore delle quali è stato trovato sperimentalmente. Questa è la legge di parete dovuta a Von Karman (1930). B è inversamente proporzionale alla rugosità delle pareti. In Figura 6.3.3 - Profilo della velocità nel sottostrato laminare e nella regione ‘log-law’ è mostrato il confronto tra la legge logaritmica e la simulazione numerica diretta. Si può notare la perfetta rispondenza tra le due curve per valori di y+>30.

Figura 6.3.3 - Profilo della velocità nel sottostrato laminare e nella regione ‘log-law’

La regione tra il sottostrato laminare e la regione ‘log-law’ è detta ‘buffer layer’. Questa è la regione di transizione tra la parte del flusso dominato dalle forze viscose e la parte del flusso dominato dalle forze di turbolenza. I valori di y+ sono compresi tra 5 e 30. 6.3.3. Lo strato esterno. La zona esterna, con valori di y+ maggiori di 50, è una zona nella quale gli effetti della viscosità molecolare sono quasi interamente trascurabili rispetto a quelli della “viscosità turbolenta”. Il profilo delle velocità è governato da una legge determinata da Von Karman (1930) e sviluppata da Millikan (1938) detta ‘velocity-defect law’. Essa devia leggermente dalla “log law” soprattutto nei boundary layers non in equilibrio con gradienti di pressione:

U max − U 1  y = − ln  + D ut k δ 

(6.3.8)

dove D è una costante che dipende dal moto. In Figura 6.3.4 - Schema riassuntivo delle regioni del flusso è proposto uno schema riassuntivo delle regioni che contraddistinguono un flusso entro un canale con pareti solide, in termini di y+ e di y/δ dove y rappresenta la distanza dalla parete mentre δ è la grandezza caratteristica del canale.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 50

6 – Le simulazioni numeriche. Generalità.

Figura 6.3.4 - Schema riassuntivo delle regioni del flusso [adattato da Pope (2000)]

6.4. Condizioni al contorno. Le condizioni al contorno assegnate ai confini del dominio fluido sono le seguenti: “Pressure-Inlet”: è stata assegnata per ogni simulazione nelle sezioni “fluide” (pressin1, pressin2, pressin3) all’interno del serbatoio. Si devono specificare: 1. la sovrapressione (gauge total pressure) rispetto alla pressione di riferimento (pressione ambiente come precedentemente spiegato); 2. la temperatura totale (coincidente con la statica per il fluido quasi in quiete nelle sezioni scelte poste a debita distanza dalla sezione di ingresso alla valvola); 3. intensità della turbolenza “turbulence intensity” che è il rapporto tra il valore medio delle fluttuazioni turbolente e la velocità media (è stato assunto un valore pari a 0.1% per i motivi sopra detti); 4. “turbulent viscosity ratio” che è un parametro adimensionale dato dal rapporto tra la viscosità turbolenta µt e la viscosità dinamica del fluido µ µ t.v.r = t µ per flussi all’interno di condotti questo parametro assume valori compresi tra 1 e 10. Nel nostro caso è stato assegnato un valore pari a 10 perché la fenomenologia nota mostra una forte prevalenza degli effetti turbolenti. • “Pressure-Outlet”: all’uscita del dominio di calcolo (sulle facce pressout1, pressout2, pressout3, pressoutup1, pressoutup2, pressoutup3, pressoutup4, pressoutup5). E’ stato inserito il valore della pressione ambiente nelle sezioni di “uscita”. • “Wall”: per gli elementi del dominio che rappresentano le pareti, cioè corpo valvola (interno ed esterno), otturatore con guarnizione di tenuta. Per queste superfici la condizione di “Wall” equivale ad imporre le condizioni delle componenti parallele alla parete della velocità entrambe nulle (no-slip). Dato che si sono imposte distribuzioni uniformi delle temperature sulle diverse pareti, e non considerando trasferimenti di calore nel materiale solido della valvola e del serbatoio, si è definita soltanto la rugosità superficiale tenendo conto della tipologia delle lavorazioni meccaniche tradizionali (tornitura, fresatura, foratura) applicate per la costruzione. • “Simmetry”: per le sezioni dei piani di simmetria sopra indicate. FLUENT attribuisce flusso nullo a tutte le grandezze di trasporto attraverso un piano di simmetria. Non esiste cioè flusso convettivo: la componente normale (al piano) della velocità è nulla. Non esiste flusso diffusivo: i gradienti (normali al piano) di tutte le variabili sono nulli. •

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 51

6 – Le simulazioni numeriche. Generalità. Le condizioni di simmetria possono essere così riassunte: 1. velocità normale al piano di simmetria pari a zero; 2. gradienti di tutte le variabili (grandezze fisiche) normali al piano di simmetria pari a zero. Dato che non esiste shear stress (sforzo di taglio) sul piano di simmetria, quest’ultimo può essere interpretato come una parete senza attrito (slip wall) se utilizzato nel calcolo di tipo “viscoso”.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 52

6 – Le simulazioni numeriche. Generalità.

6.5. Condizioni iniziali. Per risolvere numericamente il problema, si devono specificare le “condizioni iniziali” (c.i. in seguito). Queste coincidono con una inizializzazione del campo fluido-termodinamico. In pratica si devono specificare in Fluent i valori puntuali, nel dominio di calcolo, assunti dalle grandezze fisiche e corrispondenti ad un valore di primo tentativo (initial guess). Questo compito è demandato completamente all’operatore. In realtà, la scelta delle condizioni segue criteri e regole basate sulle necessità derivanti dal calcolo da produrre. Si può infatti usare la soluzione di un flusso con caratteristiche più “semplici” nella stessa geometria. Ad esempio si può partire da una simulazione “stazionaria” o nella quale è negata la viscosità del fluido. In sintesi, è buona norma e talvolta necessario che la condizione iniziale si avvicini il più possibile alla soluzione (approssimata). Questa affermazione implica che il problema sia stato “ben posto”. Con questo si intende che la scelta dei limiti fisici del dominio di calcolo, la griglia di calcolo, le caratteristiche del fluido e del solido e le condizioni al contorno specificate siano rispondenti alla fisica del fenomeno. La necessità di una corretta inizializzazione deriva da implicazioni di carattere numerico che coinvolgono direttamente la “convergenza” del calcolo iterativo effettuato. Infatti una condizione iniziale “mal posta” può facilmente non condurre alla convergenza. Quindi il nostro sistema “discretizzato” non fornisce soluzione. Nel presente caso di efflusso da una valvola in condizioni di “flusso critico” è stata posta cura nell’inizializzazione del campo fluidodinamico assegnando valori iniziali di pressione, velocità (modulo, direzione e verso mediante le tre componenti spaziali), temperatura, intensità della turbolenza e “turbulent viscosity ratio” operando per zone; cioè dividendo idealmente il dominio fluido in zone diverse (regions). Questa caratteristica di FLUENT permette una operazione di “patching” con successive sovrapposizioni di condizioni iniziali. Si inizia infatti con c.i. identiche per tutto il dominio fluido. Successivamente si effettuano le sovrapposizioni nelle zone preventivamente individuate e contrassegnate (marked) con altre c.i. che sostituiscono quelle ivi precedentemente assegnate. Si riesce in tal modo ad assegnare gradienti (ovviamente a gradini) per tutte le grandezze. I valori di primo tentativo sono stati scelti in base a simulazioni precedenti effettuate su griglie molto rade (∼100000 celle) e con il solutore “segregato”. Ciò ha comportato un calcolo preventivo con l’efflusso subsonico per via delle caratteristiche del solutore sopra detto, incapace di trattare flussi supersonici o localmente supersonici. I valori poi assegnati sono stati maggiorati per tenere conto dei salti di pressione maggiori corrispondenti alle prove in condizioni di “efflusso critico”. Mostriamo ora come è stato applicato l’operazione di “patching” nel caso in studio facendo riferimento alla Figura 6.5.1 - Quote per il “patching” dove si indicano i valori delle coordinate relative alle zone scelte per l’applicazione delle diverse condizioni iniziali. Come primo input si sono imposte le seguenti condizioni sull’intero dominio (vedi Figura 6.5.2 - I set C.I.): § velocità nulla; § pressione (relativa) nulla: § temperatura: 300K § energia cinetica turbolenta: 1m2/s2; § tasso di dissipazione dell’energia cinetica turbolenta: 1m2/s3. I valori assegnati a queste due ultime grandezze non influenzano sensibilmente la storia iniziale di convergenza; non si dimostrano critici come pressione e velocità.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 53

6 – Le simulazioni numeriche. Generalità. Secondo set di c.i. applicate alla regione cilindrica avente raggio pari a quello del mantello esterno della valvola e coincidente con la feritoia e limitato dall’otturatore (vedi Figura 6.5.3 - II set C.I.): § componente vx positiva (modulo assunto da simulazioni precedenti come sopra esposto); § pressione relativa positiva (valore leggermente inferiore a quella del serbatoio assunto da simulazioni precedenti). E’ scelto in funzione del valore da attribuire alla zona del condotto verticale della valvola. Terzo set di c.i. applicate alla regione cilindrica avente raggio pari a quello del ciglio della sede e limitato alla quota z della sede stessa (vedi Figura 6.5.4 - III set C.I.): § componente vx nulla (annulla il valore positivo precedentemente assegnato); § componente vz positiva (modulo assunto da simulazioni precedenti come sopra esposto); § pressione relativa positiva intermedia tra quella imposta nel II set di condizioni e quella del serbatoio (valore assunto da simulazioni precedenti come sopra esposto); Quarto set di c.i. applicate alla regione corrispondente alla zona del serbatoio (vedi Figura 6.5.5 - IV set C.I.): § pressione relativa positiva (pari alla pressione di prova nel serbatoio). Una volta effettuati tutti i cicli iterativi e raggiunta la convergenza, il risultato ottenuto può essere usato per l’inizializzazione di simulazioni che ad esempio comportano l’imposizione di una pressione superiore nel serbatoio. Per talune simulazioni si è ricorso a questa metodologia per abbreviare il tempo necessario all’esecuzione di altre simulazioni. Ovviamente sono state modificate le altre condizioni al contorno laddove necessario.

Figura 6.5.1 - Quote per il “patching”

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 54

6 – Le simulazioni numeriche. Generalità.

Figura 6.5.2 - I set C.I.

Figura 6.5.3 - II set C.I.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 55

6 – Le simulazioni numeriche. Generalità.

Figura 6.5.4 - III set C.I.

Figura 6.5.5 - IV set C.I.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 56

7 – Il procedimento di calcolo.

7. Il procedimento di calcolo. 7.1. Le simulazioni preliminari. Come affermato nel paragrafo 6.5, per l’inizializzazione del problema numerico si sono utilizzati risultati ottenuti da simulazioni preliminari. In questa fase, definita Fase-1, si è cercato di effettuare il maggior numero di semplificazioni possibili nella geometria del modello e di variazioni mirate nelle impostazioni del codice di calcolo per valutare e verificare l’importanza di diversi parametri. Per gradi si sono introdotte le variazioni geometriche che hanno portato il modello alle forme delle simulazioni definitive con le ridotte approssimazioni già esposte nel paragrafo 5.2. In linea di massima sono state effettuate simulazioni preliminari con la stessa sequenza di variazioni nelle impostazioni del codice al variare delle geometrie. Le informazioni ottenute dalle simulazioni precedenti hanno fornito le indicazioni per le successive. Non riporteremo le sequenza cronologica completa, come pure non accompagneranno tale documento tutti i files generati. Si riporteranno (e si accluderanno) soltanto quelli significativi. Sottolineiamo come alla luce dei risultati che verranno in seguito mostrati, una parte consistente di questi “tentativi” abbia non solo un ruolo per così dire “accademico”, ma costituisca anche e soprattutto un’investigazione necessaria alla validazione del procedimento generale per l’effettuazione di vere e proprie “prove numeriche” attendibili ed economicamente accettabili. Le principali approssimazioni geometriche sopra dette sono state le seguenti: 1. le dimensioni e le forme del dominio fluido sono state oggetto di diverse modifiche per la valutazione della loro importanza agli effetti dell’accuratezza e della convergenza (i limiti del dominio sono stati allontanati gradualmente dalla valvola); 2. in diverse simulazioni non si è esteso il dominio di calcolo alla zona del serbatoio, ma si è disegnato invece un semplice condotto diritto di sezione circolare di diametro pari a quello nominale della valvola. (vedi Figura 7.1.2 relativa alla geometria del Caso 5-1---147917; 3. il piattello è stato “disegnato” piano e senza il bordino: la superficie inferiore del piattello è alla stessa quota del bordo superiore delle feritoie di uscita (vedi Figura 7.1.6, Figura 7.1.7 e Figura 7.1.8. Le impostazioni del solutore e dei “modelli” matematici sono state le seguenti: 1. adozione del solutore “segregato” (con tutti i limiti riguardo all’impossibilità di trattare flussi supersonici e/o localmente supersonici) per una valutazione preliminare del campo fluidodinamico; 2. tentativo di utilizzare lo schema di discretizzazione (approssimazione) del secondo ordine. Quando il flusso è allineato con la griglia (ad esempio un flusso laminare in un condotto rettangolare modellato con griglia ad elementi quadrilateri o esaedrici) la discretizzazione “upwind” del primo ordine può essere sufficiente. Quando il flusso non è allineato con la griglia (cioè quando il flusso attraversa la griglia obliquamente), tale discretizzazione è affetta da un errore numerico (diffusione numerica). Per le griglie ad elementi triangolari e tetraedrici, dato che il flusso non è allineato mai con la griglia, si ottengono risultati generalmente più accurati Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 57

7 – Il procedimento di calcolo. usando la discretizzazione del secondo ordine. Anche per griglie ad elementi quadrilateri o esaedrici in genere si ottengono migliori risultati con il secondo ordine, specialmente per flussi complessi. Mentre di solito il primo ordine produce la migliore convergenza, fornisce per contro risultati meno accurati, specialmente su griglie ad elementi triangolari o tetraedrici. Per la maggior parte dei casi si può usare lo schema del secondo ordine dall'inizio del calcolo. A volte può essere conveniente partire con uno schema del primo ordine, per poi passare allo schema del secondo ordine dopo alcune iterazioni. Per esempio, per flussi ad alto numero di Mach aventi una soluzione iniziale diversa da quella finale attesa, si possono effettuare diverse iterazioni con lo schema del primo ordine e poi si attiva lo schema del secondo ordine continuando il calcolo fino alla convergenza. Per flussi semplici “allineati” con la griglia (ad esempio, flusso laminare in un condotto rettangolare modellato con elementi quadrilateri o esaedrici), la diffusione numerica sarà naturalmente bassa. Si può usare lo schema del primo ordine senza perdita significativa dell'accuratezza. Nel caso si incontrino difficoltà di convergenza con lo schema al secondo si può passare al primo. Nel presente lavoro lo schema al secondo ordine è stato abbandonato dopo diversi tentativi proprio per problemi di convergenza. Si sono comunque ottenuti risultati accettabili con lo schema del primo ordine. Si è ipotizzato che sebbene le diverse griglie utilizzate siano adeguate alla soluzione del problema, non sono sufficientemente “raffinate” (il numero degli elementi nelle zone in cui esistono alti gradienti delle grandezze fisiche non è abbastanza elevato) per la discretizzazione del secondo ordine. 3. scelta iniziale del modello di turbolenza k-ε standard, un modello semi empirico basato su equazioni di trasporto sia per l’energia cinetica turbolenta (k) che sul tasso di dissipazione (ε) di quest’ultima. L’equazione di trasporto per k coincide con la sua espressione matematica esatta. Nell’equazione per la ε compaiono invece dei coefficienti in alcuni termini. Questi valori di default per Fluent sono stati determinati da esperimenti con aria ed acqua in “tipici” flussi turbolenti che includono flussi omogenei e flussi isotropici con turbolenza indotta da griglie. Si è sperimentato che tali coefficienti permettono una buona approssimazione anche nel caso di flussi confinati entro pareti solide ed anche in flussi esterni. Essi possono essere modificati da un utente esperto. Il modello suddetto è valido solo per flussi completamente turbolenti e laddove l’effetto della viscosità molecolare può essere trascurato. Notiamo come tale modello sia molto diffuso nell’ambito delle simulazioni a carattere industriale e che si ha un’ampia casistica di applicazioni. Sia la sua accuratezza che la sua “robustezza”, cioè la sua capacità intrinseca di “risolvere” problemi appartenenti a tipologie abbastanza diverse, sono documentati in modo esauriente. Riassumendo, si è operato scegliendo i limiti del dominio di calcolo, la necessaria rappresentazione geometrica della valvola, le impostazioni del solutore e il modello di turbolenza. Tale iter preliminare è stato adottato per le simulazioni della prova sperimentale denominata 2/a a 497000Pa di pressione (relativa) nel serbatoio. Dalle relazioni sulle prove sperimentali atte a rilevare le caratteristiche di funzionamento e di portata sono stati desunti i dati delle grandezze fisiche che costituiscono le condizioni al contorno, come esposto nel paragrafo 6.4.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 58

7 – Il procedimento di calcolo. Fase-1

7.1.1. Fase-1 - Caso 11-1---174010. Nelle prime simulazioni la geometria del dominio è quella mostrata in Figura 7.1.1 L’otturatore è piatto come nei casi successivi Caso 5-1---147917 e Caso 5-2---160896 (vedi Figura 7.1.6). Il nome del “caso” è 11-1---174010. Le cartelle presenti nel DVD sono le seguenti: 1) 88100Pa-turb-segregated; 2) 88300Pa-turb-segregated-NO-ADAPT; 3) 88400Pa-turb-segregated-NO-ADAPT. Il termine “caso”, che ricorrerà in seguito, fa riferimento al termine inglese “case” con cui in Fluent sono memorizzati i files (con estensione “.cas”) delle simulazioni. In realtà nei files.cas sono memorizzate tutte le impostazioni del solutore e la geometria (importata)1 del “caso” in studio. I dati numerici delle simulazioni sono memorizzati nei rispettivi files con estensione “.dat” sempre generati da Fluent.

Figura 7.1.1

I due volumi, quello del serbatoio e quello dell’ambiente esterno hanno, nell’ordine, le seguenti dimensioni: 1) serbatoio: raggio pari a 0.35m e altezza pari a 0.35m; 2) volume ambiente esterno: raggio pari a 0.6m (circa 42 diametri) e altezza pari a 0.5m. 1

I files delle geometrie generati da Gambit hanno l’estensione “.dbs”; tali files convertiti per la successiva importazione in Fluent hanno estensione “.msh”. Questi possono già includere la definizione delle diverse entità in funzione delle condizioni al contorno che vi si debbono assegnare.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 59

7 – Il procedimento di calcolo. Fase-1 Su questa geometria si sono effettuate simulazioni con il solutore segregato raggiungendo un numero di Mach massimo pari all’unità per una pressione nel serbatoio pari a 88400Pa. Il campo fluido-termodinamico è però qualitativamente non accettabile come anche la storia di convergenza. L’inizializzazione è stata effettuata imponendo valori delle incognite (velocità, pressione, temperatura, k ed ε) costanti in tutto il dominio. 7.1.2. Fase-1 - Caso 5-1---147917. La geometria seguente relativa al caso 5-1---147917 è quella mostrata nella Figura 7.1.2. In questo “caso” è stato abolito il serbatoio per facilitare la convergenza. Si riduce così l’onere di calcolo nella zona d’imbocco nel condotto nella cui sezione d’entrata s’impone più semplicemente il valore della pressione del serbatoio. Il condotto e il volume dell’ambiente esterno hanno, nell’ordine, le seguenti dimensioni: 1) condotto: diametro pari a 14mm e lunghezza pari a 50mm; 2) volume ambiente esterno: lunghezza massima pari a 0.31m (circa 22 diametri), larghezza pari a 0.16m e altezza pari a 0.5m. E’ stata raggiunta la pressione nel serbatoio pari a 497000Pa in virtù dell’utilizzo del solutore accoppiato e si è ottenuto un campo fluido-termodinamico che fornisce importanti indicazioni numeriche per le simulazioni seguenti e soprattutto per l’inizializzazione del problema. Le cartelle presenti nel DVD sono le seguenti: 1) 5-1-100Pa-segregated 2) 5-1-45000Pa-segregated; 3) 5-1-497000Pa-coupled. I valori ottenuti nel calcolo con il solutore segregato hanno fornito l’inizializzazione per il solutore accoppiato.

Figura 7.1.2

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 60

7 – Il procedimento di calcolo. Fase-1 7.1.3. Fase-1 - Caso 5-2---160896. La geometria seguente relativa al caso 5-2---160896 è quella mostrata nella Figura 7.1.3. Dalle indicazioni ottenute dai casi precedenti riguardo all’evoluzione della storia della convergenza, è stato “reintrodotto il serbatoio” affinché il campo fluidodinamico all’imbocco del canale sia risolto dal calcolatore. Questo per simulare più fedelmente le condizioni reali e per non alterare con valori imposti di pressione non corretti il campo fluidodinamico. Il volume del serbatoio e il volume dell’ambiente esterno hanno, nell’ordine, le seguenti dimensioni: 1) serbatoio: raggio pari a 0.1m e altezza pari a 0.1m; 2) volume ambiente esterno: lunghezza massima pari a 0.31m (circa 22 diametri), larghezza pari a 0.15m e altezza pari a 0.11m.

Figura 7.1.3

Si mostrano di seguito alcune figure della geometria della valvola del caso 5-2---160896. I particolari messi in evidenza sono identici a quelli dei casi precedenti a meno di differenze sulle griglie.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 61

7 – Il procedimento di calcolo. Fase-1

Figura 7.1.4

Figura 7.1.5

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 62

7 – Il procedimento di calcolo. Fase-1

Figura 7.1.6

Figura 7.1.7

Figura 7.1.8

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 63

7 – Il procedimento di calcolo. Fase-1 Da questa simulazione si sono ottenuti risultati quantitativamente e qualitativamente accettabili. Infatti la storia di convergenza del caso 5-2---160896 – 497000Pa mostra, nel complesso, (vedi Figura 7.1.9) un andamento regolare e gradatamente decrescente.

Figura 7.1.9

Anche i campi di velocità (vedi Figura 7.1.10), Mach (vedi Figura 7.1.11), pressione (vedi Figura 7.1.12) totale e lo sviluppo dei getti (vedi Figura 7.1.13) sono accettabili. Per la geometria considerata la fisica dell’efflusso appare nel complesso correttamente calcolata (simulata).

Figura 7.1.10

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 64

7 – Il procedimento di calcolo. Fase-1

Figura 7.1.11

Figura 7.1.12

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 65

7 – Il procedimento di calcolo. Fase-1

Figura 7.1.13

Si sono quindi effettuate le simulazioni per gli altri valori di pressione: 919000Pa e 1218000Pa corrispondenti rispettivamente ai valori delle prove sperimentali 2/b e 2/c. Nella Figura 7.1.14 si riportano le relative storie di convergenza. I restart a 919000Pa

II restart a 1218000Pa

Figura 7.1.14

Per questi due valori di pressione il campo fluido-termodinamico del caso precedente, alla pressione più bassa, ha fornito l’inizializzazione del calcolo. Si modificano semplicemente le condizioni al contorno e si riavvia (si effettua il cosiddetto “restart”) il ciclo di iterazioni successivo. Nella Figura 7.1.14 ai picchi nella storia di convergenza corrispondono le fasi iniziali (al restart) del calcolo nel quale, avendo imposto nuovi e più alti valori della pressione sulle

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 66

7 – Il procedimento di calcolo. Fase-1 superfici “fluide” interne del serbatoio, si vengono a determinare salti bruschi delle grandezze fisiche. Con il procedere del calcolo il campo fluido-termodinamico si “assesta” e le curve degradano con andamento mediamente monotono. Dal confronto con i valori delle portate reali si sono osservati gli scostamenti riportati in Tabella 1 Tabella 1

Pressione 497000Pa 919000Pa 1218000Pa

Caso 5-2---160896 Portata misurata Portata calcolata (kg/h) (kg/h) 588 566 1027 963 1297 1244

Scarto percentuale (%) -3.74 -6.23 -4.08

Le cartelle presenti nel DVD sono le seguenti: 1) 5-2---160896-coupled-497000Pa; 2) 5-2---160896-coupled-919000Pa; 3) 5-2---160896-coupled-1218000Pa. Nota: è stata utilizzata la medesima griglia per tutti i tre valori di pressione e il valore dell’alzata (rispetto alla sede) dell’otturatore è sempre pari a 5mm.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 67

7 – Il procedimento di calcolo. Fase-2

7.2. Fase-2 - Affinamento della geometria del modello. Alla luce dei risultati preliminari si è deciso di apportare variazioni alla geometria del modello, ed in particolare alla zona dell’otturatore ed alla posizione della superficie di confine perpendicolare al getto. Nella Fase 2 l’otturatore è stato rappresentato con le sue forme reali (non più piano) e la superficie di confine suddetta è stata spostata a circa cento diametri dall’asse valvola (circa 1.4m) in direzione x. Si tralasciano diverse simulazioni intermedie in cui è stata valutata l’opportunità della scelta dello schema di discretizzazione del primo ordine. La denominazione del caso della geometria in questione è Caso 6-3---171604. Nelle figure seguenti si mostrano le geometrie del dominio, della valvola e dei particolari della zona dell’otturatore in prossimità delle feritoie di scarico. Si nota la forma più complessa dell’otturatore.

Figura 7.2.1

Figura 7.2.3 Figura 7.2.2

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 68

7 – Il procedimento di calcolo. Fase-2

Figura 7.2.4

Figura 7.2.5

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 69

7 – Il procedimento di calcolo. Fase-2

Figura 7.2.6

Figura 7.2.7

Figura 7.2.8

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 70

7 – Il procedimento di calcolo. Fase-2 Anche da queste simulazioni si sono ottenuti risultati quantitativamente e qualitativamente accettabili. Infatti la storia di convergenza del caso 6-3---171604 – 497000Pa mostra (vedi Figura 7.2.9) un andamento regolare e gradatamente decrescente.

Figura 7.2.9

Si fa notare che l’inizializzazione del campo è stata effettuata come spiegato al paragrafo 6.4, con l’ausilio del “patching”. I valori assegnati alle grandezze fisiche sono stati desunti dai valori ottenuti dalle simulazioni del caso 5-2---160896 per lo stesso valore di pressione nel serbatoio. Nelle diverse “regioni” sono stati assegnati valori corrispondenti circa pari a quelli medi riscontrati. Lo sviluppo dei getti (vedi Figura 7.2.10), i campi di velocità (vedi Figura 7.2.11), Mach (vedi Figura 7.2.12), pressione totale (vedi Figura 7.2.13) ed energia cinetica turbolenta (vedi Figura 7.1.14) sono accettabili. Per la geometria considerata la fisica dell’efflusso appare, anche in questa occasione, correttamente calcolata (simulata). Al paragrafo 7.3, si approfondirà lo studio della fenomenologia dell’efflusso e si valuteranno le differenze con l’efflusso delle simulazioni della Fase-1.

Figura 7.2.10

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 71

7 – Il procedimento di calcolo. Fase-2

Figura 7.2.11

Figura 7.2.12

Figura 7.2.13

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 72

7 – Il procedimento di calcolo. Fase-2

Figura 7.2.14

Si nota come la presenza del bordino dell’otturatore e della camera anulare interna superiore alle feritoie modifichi il campo rispetto alle simulazioni precedenti. Si sono poi effettuate le simulazioni per gli altri valori di pressione (919000Pa e 1218000Pa). Anche per questi due valori di pressione il campo fluido-termodinamico del caso precedente, a pressione inferiore, ha fornito l’inizializzazione del calcolo. Dal confronto con i valori delle portate reali si sono osservati gli scostamenti riportati in Tabella 2. Tabella 2

Pressione 497000Pa 919000Pa 1218000Pa

Caso 6-3---171604 Portata misurata Portata calcolata (kg/h) (kg/h) 588 581.7 1027 994.1 1297 1000.5

Scarto percentuale (%) -1.12 -3.51 -1.17

Rispetto alle simulazioni della Fase-1, relative al Caso 5-2---160896, si è ottenuta una prima riduzione degli scarti dalle prove sperimentali. Vedremo nel paragrafo seguente come gli affinamenti successivi riducano ulteriormente tali scarti. Le cartelle presenti nel DVD sono le seguenti: 1) 6-3---171604-497000Pa; 2) 6-3---171604-919000Pa; 3) 6-3---171604-1218000Pa.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 73

7 – Il procedimento di calcolo. Fase-2 7.3. Simulazioni definitive. Analisi numerica e fisica dei risultati. Per migliorare l’accuratezza delle simulazioni si è proceduto all’analisi della griglia di calcolo specialmente in prossimità delle pareti. In particolar modo si è valutata la necessità di operare infittimenti (o diradamenti) nelle zone critiche, cioè nelle zone affette da alti gradienti delle grandezze fisiche ed in particolar modo della velocità del flusso. Come spiegato ai paragrafi 6.2 e 6.3, lo spessore degli strati di celle adiacenti alla pareti influenza la risoluzione del campo fluidodinamico nello strato limite. In base al tipo di “near-wall treatment” cioè in base al tipo di approccio usato per trattare la zona vicino alle pareti, diversi sono i requisiti per la griglia. Per le simulazioni svolte, il campo fluido-termodinamico nella “inner region” cioè nel ”sottostrato limite” e nel “buffer layer”, non è risolto in modo diretto, bensì vengono adottate funzioni semi empiriche chiamate “wall functions”. Esse funzionano da “ponte” tra la parete e la regione completamente turbolenta. L’uso delle “wall functions” ovvia al bisogno di modificare i modelli di turbolenza per tenere conto della presenza delle pareti. Il metodo per controllare e modificare se necessario lo spessore degli strati suddetti consiste nell’analisi del campo fluido-termodinamico e nella valutazione della y+. Si è proceduto quindi all’analisi dei risultati del Caso 6-3---171604 – 497000Pa come esposto nelle figure di seguito riportate, cercando di evidenziare i punti critici nell’evoluzione spaziale del flusso. Di seguito riportiamo le figure (dalla Figura 7.3.1 alla Figura 7.3.7) relative alla distribuzione del numero di Mach su di un fascio di piani passanti per l’asse della valvola, iniziando dal piano di simmetria x-z fino al piano denominato plane-75 passante vicino alla parete verticale della feritoia.

Figura 7.3.1

Figura 7.3.2

Figura 7.3.3

Figura 7.3.4

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 74

7 – Il procedimento di calcolo. Fase-2

Figura 7.3.5

Figura 7.3.6

Figura 7.3.7

Dalle figure si notano le zone in cui il fluido è sottoposto alle maggiori accelerazioni. Le pressioni in gioco e la geometria interna “tormentata” obbligano il fluido ad un percorso complesso (meglio evidenziato nelle figure seguenti, le quali visualizzano le linee di flusso colorate in base alle diverse grandezze fisiche). Si può osservare (Figura 7.3.8) come il getto si allarghi (con simmetria rispetto al piano xz) allontanandosi dalla valvola (esso si espande in virtù dei gradienti di pressione).

Figura 7.3.8

Internamente alla valvola il flusso si suddivide in diversi rami principali alcuni dei quali evidenziano un alto effetto di “swirl”. Con riferimento alla geometria simulata (un quarto della valvola reale) le linee di flusso centrali del condotto si dividono all’uscita del bordino dell’otturatore in due flussi contro rotanti (vedi Figura 7.3.9 e Figura 7.3.10), in special modo in prossimità del piano di simmetria y-z. Il primo ramo (con riferimento alle figure suddette ed alle Figura 7.3.11, Figura 7.3.12, Figura 7.3.13, e Figura 7.3.14) rotante in senso antiorario al di sotto dell’otturatore, si dirige

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 75

7 – Il procedimento di calcolo. Fase-2 verso la feritoia di uscita con moto approssimativamente elicoidale nella regione anulare inferiore compresa tra la superficie interna del corpo valvola, l’otturatore e l’esterno della sede. Regione anulare superiore

Regione anulare inferiore

Figura 7.3.9

Regione anulare superiore

Regione anulare inferiore

Figura 7.3.10

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 76

7 – Il procedimento di calcolo. Fase-2

Figura 7.3.11

Figura 7.3.12

Figura 7.3.13

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 77

7 – Il procedimento di calcolo. Fase-2

Figura 7.3.14

Il secondo flusso, ruota in senso orario (con rotazione mediamente più ampia ed allungata) all’interno della regione anulare al di sopra del bordino compresa tra la superficie interna del corpo valvola e quelle superiori dell’otturatore (vedi Figura 7.3.9, Figura 7.3.10, e le “trasparenze” in Figura 7.3.15, Figura 7.3.16, e Figura 7.3.17.

Figura 7.3.15

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 78

7 – Il procedimento di calcolo. Fase-2

Figura 7.3.16

Figura 7.3.17

Il primo vortice, nell’avvicinamento alla feritoia, viene stirato, distorto e ridotto nelle sue dimensioni medie radiali dal flusso uscente dalla regione anulare superiore con conseguente accelerazione (come si evince dalla Figura 7.3.3 Figura 7.3.4, 7.3.5 e 7.3.6, nella zona intermedia tra i piani di simmetria e visibile grazie alle Figura 7.3.13 e Figura 7.3.14). Inoltre, in prossimità della feritoia, il flusso uscente dalla regione anulare superiore schiaccia e devia verso il basso il flusso principale già deviato dal bordino. Ne risulta una forte accelerazione con il raggiungimento dei valori più elevati del numero di Mach nella zona esterna alla feritoia. Tale accelerazione contribuisce allo “stiramento” del vortice sopra descritto. Le linee di flusso adiacenti alla superficie del condotto che oltrepassano la sede sono deviate e ruotate verso il basso nelle vicinanze della sede stessa. Sempre con riguardo alle linee di flusso che lambiscono la superficie esterna della sede, si nota che in prossimità del piano di simmetria x-z, esse si sollevano, decelerando, in direzione verticale su tale piano. Ciò dipende dal flusso proveniente dall’altra metà (simmetrica) della valvola. Intorno alla sede, davanti alla feritoia i due flussi convergono verso il piano di simmetria. Da tale convergenza (e relativa decelerazione, come visibile nelle Figura 7.3.1 Figura 7.3.2 e Figura 7.3.17) risulta la Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 79

7 – Il procedimento di calcolo. Fase-2 deviazione verso l’alto e l’immediata rotazione nella direzione (radiale) ad opera del flusso principale uscente (vedi anche Figura 7.3.18).

Figura 7.3.18

Anche la “costrizione” dovuta al flusso ora descritto è una concausa dell’accelerazione del flusso principale. L’area efficace (reale) di efflusso risulta mediamente inferiore alla geometrica a disposizione. Dalle Figura 7.3.1, Figura 7.3.2, Figura 7.3.6, Figura 7.3.10 e Figura 7.3.12 si osserva come, in corrispondenza della feritoia, nella zona interna della sede dove il condotto si raccorda con il ciglio (con raccordo circolare), il flusso sia accelerato vigorosamente. Si nota una vasta zona nella quale il numero di Mach locale è superiore all’unità. Proprio la presenza della camera anulare superiore, e la presenza del bordino modificano nettamente il flusso rispetto alla geometria con otturatore piatto usato nella Fase-1. Infatti la camera anulare permette una via secondaria per una portata aggiuntiva che aiuta ed è aiutata ad essere smaltita dal flusso deviato dall’otturatore. Inoltre il bordino risulta dalla presenza di un incavo cilindrico inferiormente all’otturatore. Tale incavo permette un’area efficace di efflusso maggiore. Per confronto mostriamo ora due figure relative al Caso 5-2---160896 a 497000Pa della Fase-1.

Figura 7.3.19

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 80

7 – Il procedimento di calcolo. Fase-2

Figura 7.3.20

Sull’esterno del corpo valvola il campo fluidodinamico non presenta caratteristiche particolari. Le pareti verticali delle feritoie generano ricircoli esterni di entità limitata. Il richiamo di aria dall’ambiente esterno (e che lambisce anche l’esterno del corpo valvola) contribuisce alla riduzione di quest’ultimi (vedi Figura 7.3.21).

Figura 7.3.21

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 81

7 – Il procedimento di calcolo. Fase-2

7.4. Affinamento delle griglie di calcolo. Di seguito riportiamo le figure che rappresentano la distribuzione della y+ sulle superfici nella zona critica dell’interno della valvola ed in prossimità delle feritoie di scarico. Le due figure hanno la sola differenza nella scala attribuita alla colorazione. Infatti nella prima si è imposto un valore massimo pari a 150 (tutte le zone nel quale la y+ ha un valore maggiore sono colorate con lo stesso tono di colore rosso) per meglio evidenziare le zone più critiche (in questo modo sono esaltate anche le zone con valori più bassi). Nella seconda figura sono rappresentati tutti i valori calcolati lasciando la scalatura automatica del colore.

Figura 7.4.1

Figura 7.4.2

Come si osserva, diverse zone (superfici) denunciano valori dell’y+ circa pari o superiori a 100÷150. Come descritto al paragrafo 6.3, sarebbe auspicabile riscontrare valori di y+>30÷60. Un valore della y+ vicino al limite inferiore (y+ circa uguale a 30) è consigliato (da manuale di Fluent). In realtà sono ancora accettabili valori come quelli sopra trovati ma non nelle zone più critiche. Per ridurre la y+ si è proceduto all’affinamento (o adapting) della mesh di calcolo. Tale operazione è stata effettuata inizialmente sulla mesh del Caso 6-3---171604 e per il valore di pressione pari a 497000Pa. Per l’affinamento è stato scelto di raddoppiare gli elementi costituenti

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 82

7 – Il procedimento di calcolo. Fase-2 i primi due strati adiacenti alle pareti interne ed esterne alla valvola. Di seguito riportiamo le figure che mostrano particolari delle griglie prima e dopo l’adapt.

Figura 7.4.3

Figura 7.4.4

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 83

7 – Il procedimento di calcolo. Fase-2

Figura 7.4.5

Figura 7.4.6

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 84

7 – Il procedimento di calcolo. Fase-2

Figura 7.4.7

Figura 7.4.8

Come si intuisce, l’operazione di affinamento ha comportato l’aumento del numero degli elementi della griglia. Si è passati da 171604 a 384019 celle. Alla simulazione effettuata con l’inizializzazione del campo con la stessa procedura di “patching” del Caso 6-3---171604 – 497000Pa, è stato assegnato il nome: 6-3-boundadapt2celle---384019 – 497000Pa. Osserviamo ora le variazioni ottenute nella distribuzione della y+ nelle zone critiche. Anche in questo caso si faccia attenzione alle diverse scale usate per la rappresentazione a colori.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 85

7 – Il procedimento di calcolo. Fase-2

Figura 7.4.9

Figura 7.4.10

Figura 7.4.11

Figura 7.4.12

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 86

7 – Il procedimento di calcolo. Fase-2 Per confronto con le precedenti Figura 7.4.1, e Figura 7.4.2 si verifica un abbassamento generalizzato dei valori nelle aree affette da alte velocità del flusso come sopra descritto (vedi per esempio la zona interna in prossimità della sede. L’affinamento ha comportato un miglioramento dell’accuratezza dato che lo scarto della portata calcolata da quella misurata è sceso da –1.12% a –0.51%. Si verifica inoltre come valori anche non prossimi al limite atteso di 30÷60 della y+ non comportino, per le zone non critiche, elevati problemi numerici. Per il Caso 6-3---171604-919000Pa avendo ottenuto nelle precedenti simulazioni lo scarto maggiore si è scelto di effettuare l’affinamento non sui primi due, bensì sui primi tre strati vicino alle pareti. Il caso corrispondente ha nome: Caso 6-3-boundadapt2e3---415029 – 919000Pa. Anche in quest’occasione l’affinamento ha comportato un miglioramento dell’accuratezza dato che lo scarto della portata calcolata da quella misurata è sceso da –3.51% a –2.57%. Per l’inizializzazione del campo fluido-termodinamico si è operato con l’operazione di patching. Per il Caso 6-3-boundadapt2e3---415029 – 1218000Pa invece, si è operato un “restart” dal Caso 6-3-boundadapt2e3---415029 – 919000Pa usando perciò la stesa griglia del caso precedente. Il miglioramento della griglia ha comportato una riduzione dello scarto tra portata misurata e calcolata e si è passati da –1.17% a –0.30%. Riassumiamo nella seguente Tabella 3 i risultati definitivi dopo l’affinamento della griglia per i casi corrispondenti alle “prove di qualifica” (prove 2/a a 4970000Pa, 2/b a 919000Pa e 2/c a 1218000Pa) della valvola in studio. Nota: • il numero di Reynolds è stato calcolato considerando le medie dei valori della densità e della velocità dell’aria nel condotto di entrata e come lunghezza caratteristica il suo diametro; • nella tabella si riportano anche i valori del numero di Mach massimo raggiunto nell’efflusso. Tabella 3

Pressione

Numero di Portata misurata Portata calcolata Mach max Reynolds (kg/h) (kg/h)

Scarto percentuale (%)

6-3-boundadapt2celle---384019 497000Pa

814500

919000Pa

1393800

1218000Pa

1800600

1.967 588 585.2 6-3-boundadapt2e3---415029 2.249 1027 1000.5 6-3-boundadapt2e3---415029 2.508 1297 1292.7

-0.51 -2.57 -0.30

Le cartelle presenti nel DVD sono le seguenti: 1) 6-3-boundadapt2celle---384019-497000Pa; 2) 6-3-boundadapt2e3---415029-919000Pa; 3) 6-3-boundadapt2e3---415029-1218000Pa.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 87

7 – Il procedimento di calcolo. Fase-2

7.5. Le “prove intermedie”. Per valutare il comportamento del codice di calcolo per valori della pressione nel serbatoio per i quali non erano a disposizione prove sperimentali si è proceduto all’effettuazione di due simulazioni per valori della pressione intermedi tra quelli precedentemente descritti. Più precisamente si sono scelti i seguenti valori: 700000Pa e 1050000Pa. Le condizioni al contorno sono state scelte ragionevolmente in base ai valori presenti nei rapporti di prova in nostro possesso, ipotizzando valori di temperatura nel serbatoio poco diversi. 7.5.1. Prima prova intermedia: 700000Pa. Per il caso a 700000Pa (6-3-boundadapt2celle---384019 – 70000Pa) si è semplicemente effettuato il restart dal caso, a pressione inferiore pari a 497000Pa, 6-3-boundadapt2celle--384019 – 497000Pa. Il valore della portata calcolata è 783.5kg/h per un numero di Reynolds pari a 1090000 ed un Mach massimo di 2.115. 7.5.2. Seconda prova intermedia: 1050000Pa. Per questa prova ( Caso 6-3-boundadapt2e3---415029 – 1050000Pa) si è effettuato il restart dal caso, a pressione inferiore pari a 919000Pa, 6-3-boundadapt2e3---415029 – 919000Pa. Il valore della portata calcolata è 1129.2kg/h per un numero di Reynolds pari a 1565000 ed un Mach massimo di 2.348. Si riportano le rispettive curve dei residui.

Figura 7.5.1

Figura 7.5.2

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 88

7 – Il procedimento di calcolo. Fase-2 7.6. Le “prove bassa pressione”. Per valutare il comportamento del codice di calcolo per valori della pressione nel serbatoio abbastanza diversi dal campo in cui si sono svolte le prove di qualifica del modello della valvola, si è proceduto all’effettuazione delle simulazioni delle prove a “bassa pressione” cioè per pressioni inferiori a 3bar (300000Pa). Più precisamente si sono effettuate “prove numeriche” a 50000Pa, 76000Pa, 101000Pa, 129000Pa, 147000Pa, 176000Pa, 217000Pa, 261000Pa e 309000Pa. Si è scelta una sola griglia per tutte le simulazioni, le quali partendo da quella a 50000Pa sono state ottenute con dei “restart”. La griglia usata è stata quella ottimizzata per le prove ad alta pressione a 919000Pa e 1218000Pa, cioè la 6-3-boundadapt2e3---415029. Questa scelta ha voluto essere una forzatura proprio per valutare l’attendibilità di una griglia affinata per un campo di pressioni abbastanza diverso. Come si vedrà nel prossimo capitolo, i risultati confermano l’importanza della griglia in relazione alla fisica dell’efflusso. Dal confronto con i valori delle portate reali si sono osservati gli scostamenti riportati in Tabella 4. Tabella 4

Pressione 50000Pa 76000Pa 100000Pa 129000Pa 147000Pa 176000Pa 217000Pa 261000Pa 309000Pa

6-3-boundadapt2e3---415029 Numero di Portata misurata Portata calcolata Scarto percentuale Mach max Reynolds (kg/h) (kg/h) (%) 176000 0.826 115.6 124.8 +7.4 221000 1.14 149.2 158.7 +6.0 259200 1.3 180.7 186.2 +3.0 302700 1.418 212.1 218.5 +2.9 330800 1.484 230.7 237.8 +3.0 372900 1.532 262.3 268.2 +2.2 430300 1.576 305.8 309.4 +1.2 493000 1.646 351.6 354.6 +0.8 558800 1.774 401.2 401.6 +0.1

Dai risultati numerici si osserva che: § gli scarti sono tutti positivi, con un’inversione di tendenza rispetto a quelli relativi alle prove precedenti; § esiste un salto tra i primi due valori ed i seguenti. Nei primi due casi, come vedremo meglio al capitolo seguente, si è in condizioni di salto di pressione tali da non instaurare il “choking” (il regime dell’effusso è subsonico). Le condizioni fluidotermodinamiche delle prove successive invece rientrano nella situazione di flusso “bloccato”; § gli scarti, dalla pressione di 100000Pa alla prova a 309000Pa, si riducono. Le spiegazioni dei risultati e della loro “distribuzione” possono avere origine da due fattori distinti e/o concomitanti: uso della stessa griglia per tutte le simulazioni e correlazione con la fisica dell’efflusso; prove a bassa pressione in cui è stata limitata l’alzata dell’otturatore a 5mm ma questo non è stato bloccato per quel valore di alzata.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 89

7 – Il procedimento di calcolo. Fase-2 Per quanto riguarda il primo punto possiamo affermare che, se gli scarti dipendono dalla griglia usata, questa è stata ottimizzata per una tipologia di flusso che è giocoforza diversa per i casi a bassa pressione (ed in particolare per i due in regime subsonico o transonico). La conferma deriverebbe dall’evidenza del miglioramento (riduzione) dello scarto all’aumentare della pressione di prova nel serbatoio. Tali risultati confermano l’importanza fondamentale della griglia di calcolo. Riguardo al secondo punto, per i valori più bassi di pressione, l’efflusso delle prove numeriche potrebbe essere falsato dall’aver imposto l’alzata dell’otturatore pari a 5mm quando invece essa era minore nelle condizioni di prova, sopravalutando la capacità di portata della valvola. Mostriamo nelle due figure seguenti i campi del Mach per i due casi in regime transonico.

Figura 7.6.1

Figura 7.6.2

NOTA: Per brevità, alleghiamo nel DVD diversi files in formato “.doc” nei quali riportiamo figure significative relative agli efflussi per i diversi valori di pressione nel serbatoio. In particolare essi aiutano a visualizzare l’evoluzione spaziale dell’efflusso. Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba.

Roma, 16/06/2005

Pagina 90

8 – Confronto con i risultati sperimentali.

8. Confronto con i risultati sperimentali. 8.1. Riferimenti normativi. Nel presente capitolo si farà riferimento a formule per la caratterizzazione del regime di funzionamento delle valvole e alla terminologia presente nella Raccolta E – Edizione gennaio 1979 e nella INTERNATIONAL STANDARD ISO 4126-1 – Safety Valves. In particolare ci si riferisce a: § paragrafi E.1.D.1, E.1.D.2 e E.1.D.3 del fascicolo “E.1-Accessori di sicurezza e controllo” della suddetta Raccolta; § paragrafo “7-Determination of safety valve performance” della suddetta ISO 41261. Tutti i valori numerici delle grandezze fisiche e delle portate relativi alle prove sperimentali sono stati estratti dalla “Relazione sulle prove sperimentali atte a rilevare le caratteristiche di funzionamento e di portata della valvola di sicurezza” messa a disposizione da ISPESL, e dai rapporti di prova relativi alla valvola (denominata Modello D-14) in studio. Nel seguito useremo diverse denominazioni per le serie di prove sperimentali e numeriche in funzione del campo di pressione e/o di altre caratteristiche salienti. Con il termine “prove ad alta pressione” indicheremo le prove 2/a a 497000Pa, 2/b a 919000Pa e 2/c a 1218000Pa. Le “prove intermedie” sono le due prove “numeriche” a 700000Pa e 1050000Pa, sempre nel campo ad “alta pressione”. Con il termine “prove a bassa pressione” indicheremo le nove prove a pressione inferiore a 300000Pa. Nelle norme sopra indicate si stabiliscono definizioni inerenti alle condizioni di funzionamento di un generico dispositivo per il controllo della pressione e/o dell’efflusso di gas o vapore tra due ambienti a pressione differente. Tale dispositivo generico come anche la valvola in studio è assimilabile ad un orifizio (con area di efflusso a geometria variabile per via del movimento dell’otturatore). Per questi dispositivi valgono le considerazioni di natura fisica fatte al Capitolo.4. In particolare si definiscono condizioni di “flusso o salto critico” quelle per cui risulta: k

pb  2  k −1 ≤  p  k +1

(8.1)

dove: § p è la pressione (in bar assoluti) all'entrata della valvola alla quale si misura la portata (essa è data dalla somma della pressione di taratura e della sovrapressione); § pb è la contropressione in bar assoluti a valle del dispositivo e k è l'esponente dell'equazione dell'espansione isentropica. Come già esposto in precedenza, per l’aria è k=1,4 e quindi si ha flusso critico se è:

pb ≤ 0.528 . p

(8.2)

La formula per il calcolo della portata, in queste condizioni è:

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 91

8 – Confronto con i risultati sperimentali.

qmg = 0.2883 ⋅ C

p ⋅S . ν

(8.3)

In quest’ultima: § qmg indica la portata massica teorica espressa in kg/h; § C=2.7 (è funzione del coefficiente k); § p è espresso in bar; § ? è il volume specifico alle condizioni di pressione e temperatura del serbatoio durante l’efflusso. Esso è espresso in m3/kg e, dalla legge di stato dei gas perfetti, è:

ν=

RT . p

(8.4)

Dalle formule precedenti si verifica che in tutte le prove effettuate, eccetto le prime due con i valori di pressione relativa nel serbatoio a 50000Pa e 76000Pa, la valvola è in condizioni di salto critico. Infatti, con riferimento alla prova a 50000Pa, valore di prova più basso in assoluto, si ha: § p=(0.5+0.99)bar=1.49bar; § pb=0.99bar; e quindi è pb/ p=0.664>0.528 (condizione di salto non critico). Per la prova a 76000Pa si ha: § p=(0.76+0.99)bar=1.75bar; § pb=0.99bar; Risulta perciò pb/ p=0.566>0.528 (condizione di salto non critico). Al paragrafo seguente esporremo i valori delle grandezze significative nelle condizioni delle prove ad alta pressione ed i relativi valori dei coefficienti di efflusso e delle portate (reali).

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 92

8 – Confronto con i risultati sperimentali.

8.2. Coefficienti di efflusso e portate reali. Riportiamo ora sotto forma (sintetica) di tabelle i valori delle grandezze significative nelle condizioni delle “prove ad alta pressione”. Prova 2/a Prova 2/b Prova 2/c Pressione (relativa) nel serbatoio espressa in bar e Pa:

4,97 497000

9,19 919000

12,18 1218000

bar Pa

Pressione (assoluta) nel serbatoio espressa in bar e Pa: (assoluta=relativa+pressione atmosferica)

5,96 596000

10,18 1018000

13,17 1317000

bar Pa

Pressione atmosferica in bar e Pa:

0,99 99000

0,99 99000

0,99 99000

bar Pa

287 KJ/(kmol K)

R del gas (aria): 308

Temperatura nel serbatoio in K: 2

Portata teorica in kg/h ed in kg/s: Nota: calcolate con la formula (8.3)

311

S= 153,86

Area della sezione trasversale del condotto in mm : Volume specifico ?

310

K mm

2

0,148

0,087

0,068

kg/m

759

1293

1670

kg/h

0,211

0,359

0,464

kg/s

3

Dati dichiarati da ISPESL Prova 2/a Prova 2/b Prova 2/c Coefficienti di efflusso:

0,773

0,793

0,775

Portata teorica in kg/h ed in kg/s:

761 0,211

1295 0,360

1673 0,465

588

1027

1297

kg/h

0,163

0,285

0,360

kg/s

Portata reale in kg/h ed in kg/s: (calcolata moltiplicando la portata teorica per il rispettivo coefficiente di efflusso)

kg/h kg/s

0,781 Coefficiente di efflusso del modello (vedi ALLEGATO E): (valore medio della serie di prove) -5% 0 5% Scostamento percentuale accettato: ±5% 0,742 0,781 0,820 Valori corrispondenti del coefficiente di efflusso: Nella tabella seguente si riportano i valori delle portate, in base al coefficiente di efflusso del modello e ai suoi scostamenti percentuali accettati. Essi forniscono il campo di tolleranza per le prove della valvola in studio.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 93

8 – Confronto con i risultati sperimentali.

-5% 0,742 0,781 +5% 0,820

4,97 565 594 624

Pressione 9,19 961 1011 1062

12,18 1241 1307 1372

bar Kg/h Kg/h Kg/h

Valori delle portate calcolati dalle portate teoriche e dai coefficienti di efflusso indicati a sinistra.

Nota: il calcolo fluidodinamico è stato effettuato su di un quarto della valvola scegliendo due piani di simmetria tra essi perpendicolari. La portata calcolata equivale quindi ad un quarto della totale. Per questa ragione nella tabella seguente si sono riportati i valori della portata calcolata indicati come “portata/4”. E’ stato poi indicato lo scarto percentuale della portata totale calcolata rispetto alla totale misurata.

Valori calcolati delle portate e corrispondenti coefficienti di efflusso. Calcolo dello scarto percentuale dalle portate misurate. Caso

6-3---171604-497000Pa (171604 celle-186677 nodi) Prima dell’adapt

Prova 2/a - 4,97bar (497000Pa) portata totale scarto portata/4 calcolata calcolata percentuale 0,040 0,162 kg/s 581,67 kg/h -1,12% Coefficiente di efflusso corrispondente 0,764 “ portata/4

6-3-boundadapt2celle---384019 -497000Pa (384019 celle-460684 nodi) Dopo l’adapt

Caso (restart dai 4.97bar) 6-3---171604-919000Pa (171604 celle-186677 nodi) Prima dell’adapt

scarto percentuale

portata totale

0,041

0,163 kg/s 585,24 kg/h Coefficiente di efflusso corrispondente 0,769

Dopo l’adapt

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

“

Prova 2/b - 9,19bar (919000Pa) portata totale scarto portata/4 calcolata calcolata percentuale 0,069 0,276 kg/s 994,08 kg/h -3,20% Coefficiente di efflusso corrispondente 0,768 “ portata totale scarto calcolata percentuale 0,069 0,278 kg/s 1000,55 kg/h -2,57% Coefficiente di efflusso corrispondente 0,773 “

portata/4 calcolata 6-3-boundadapt2e3---415029 -919000Pa (415029 celle-489410 nodi)

-0,51%

Roma, 16/06/2005

Pagina 94

8 – Confronto con i risultati sperimentali.

Caso

6-3---171604-1218000Pa (171604 celle-186677 nodi) Prima dell’adapt

Prova 2/c - 12,18bar (1218000Pa) portata totale scarto portata/4 calcolata calcolata percentuale 0,089 0,356 kg/s 1281,36 kg/h -1,17% Coefficiente di efflusso corrispondente 0,766 “ portata totale scarto calcolata percentuale 0,090 0,359 kg/s 1292,75 kg/h -0,30% Coefficiente di efflusso corrispondente 0,773 “

portata/4 calcolata 6-3-boundadapt2e3---415029 -1218000Pa (415029 celle-489410 nodi) Dopo l’adapt

Valori calcolati delle portate delle prove “numeriche” intermedie. Caso

6-3-boundadapt2celle---384019 -70000Pa (384019 celle-460684 nodi)

Caso

Prima prova intermedia a 7bar (70000Pa) portata totale portata/4 calcolata calcolata 0.054 0.218 kg/s 783,46 kg/h

Seconda prova intermedia a 10.5bar (1050000Pa) portata/4 calcolata

6-3-boundadapt2e3---415029 -1050000Pa (384019 celle-460684 nodi)

0,078

portata totale calcolata 0,314 1129,20

kg/s kg/h

Dalle tabelle si osserva come tutte le portate calcolate siano inferiori rispetto alle misurate. Gli scarti massimi si hanno per le griglie prima dell’operazione di affinamento. Lo scarto massimo assoluto si verifica per la prova a 919000Pa. Sono stati effettuati diversi altri tentativi numerici, ma lo scarto è rimasto sempre dello stesso ordine di grandezza. E’ stata anche effettuata una prova di sensibilità all’alzata dell’otturatore nei limiti dell’accuratezza dichiarata del sensore (trasduttore) di alzata dell’otturatore. Non si sono ottenute variazioni significative. Per brevità non esponiamo i risultati, ma alleghiamo i files corrispondenti nella cartella Fase-3-919000-Pa. Riportiamo nelle pagine seguenti i grafici dei risultati descritti fino a questo punto. Tali grafici evidenziano, accettando lo scarto percentuale del ±5% rispetto al “coefficiente di efflusso del modello”, che tutte le “prove numeriche” rientrano nel campo di tolleranza accettabile. Si disegnano i grafici delle prove sperimentali, delle prove prima dell’adapt (spezzate alle quali è stato dato il nome di “Valori minimi calcolati”), di quelle successive all’operazione di affinamento (spezzate alle quali è stato dato il nome di “Valori massimi calcolati”) e delle “prove intermedie”. Si nota come queste ultime siano perfettamente allineate sui segmenti che uniscono le prove ad esse precedenti e seguenti. Da questo risultato, si può anche ipotizzare una deviazione del valore misurato relativo alla pressione di 919000Pa dal valore reale atteso, addirittura in sede di prova. Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 95

8 – Confronto con i risultati sperimentali.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 96

8 – Confronto con i risultati sperimentali.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 97

8 – Confronto con i risultati sperimentali.

Nella tabella seguente si riportano invece i valori della portata calcolata, di quella misurata e dello scarto percentuale relativamente alle prove a “bassa pressione”. Misurati

Calcolati

Scarto percentuale portata totale (%) (kg/h)

Prova

pressione(Pa)

portata/4 (kg/s)

portata totale(kg/h)

portata/4 (kg/s)

6/a

50000

0,0080

115,6

0.0087

124,8

7,4

6/b

76000

0,0104

149,2

0.0110

158,7

6,0

6/c

100000

0,0126

180,7

0.0129

186,2

3,0

6/d

129000

0,0147

212,1

0,0152

218,5

2,9

6/e

147000

0,0160

230,7

0,0165

237,8

3,0

6/f

176000

0,0182

262,3

0,0186

268,2

2,2

6/g

217000

0,0212

305,8

0,0215

309,4

1,2

6/h

261000

0,0244

351,6

0,0246

354,6

0,8

6I

309000

0,0279

401,2

0,0279

401,6

0,1

Ed alla pagina seguente si riportano tali risultati sotto forma di grafico. Per il commento di questi risultati vedi il Capitolo7 al paragrafo7.6. Le “prove a bassa pressione”.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 98

8 – Confronto con i risultati sperimentali.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 99

9 - Conclusioni e proposta di una procedura per le “prove numeriche”.

9. Conclusioni e proposta di una procedura per le “prove numeriche”. Alla luce dei risultati esposti nei precedenti capitoli, possiamo affermare che la metodologia scelta per l’effettuazione delle “prove numeriche” si rivela corretta ed attendibile nei limiti del campo di tolleranza accettabile definito nell’ALLEGATO E. Una conferma ulteriore della corretta impostazione del problema numerico deriva dalle storie di convergenza delle simulazioni definitive sia prima che dopo l’affinamento delle griglie. In effetti trattandosi di un efflusso complesso tridimensionale, supersonico (localmente), si raggiungono valori accettabili dei residui con un numero ragionevolmente basso di cicli iterativi (4000÷6000 iterazioni, come evidente dalle relative figure nel capitolo precedente). Inoltre, come già esposto nel capitolo precedente, anche le modalità di decremento delle curve di convergenza sono buone. Un limite che tali curve mettono in evidenza in alcune simulazioni, è l’impossibilità di raggiungere limiti più bassi dei valori dei residui anche se gli scarti percentuali dei risultati sono molto bassi. In effetti si osserva nella figura seguente (relativa alla simulazione con griglia affinata a pressione pari a 1218000Pa) che raggiunto un determinato valore le curve dei residui iniziano ad oscillare.

L’oscillazione si stabilizza mantenendosi limitata in ampiezza. Questo evidenzia limiti di carattere numerico connessi alla discretizzazione del continuo: esiste un limite quantitativo e/o qualitativo nella griglia di calcolo. Probabilmente (non sono stati effettuati tentativi in tal senso) una griglia più fitta o di diversa fattura e la probabile connessa possibilità dell’utilizzo dello schema del secondo ordine potrebbero portare miglioramenti. Il presente lavoro fornisce indicazioni importanti per l’esecuzione di simulazioni su valvole analoghe, stabilendo un iter procedurale ben definito. In realtà il software Fluent si presta all’esecuzione di simulazioni basate su griglie dotate di elementi di forma diversa e/o su diverse modelli per la soluzione della turbolenza, ma alla base del calcolo sussistono sempre le stesse basi fisico-matematiche.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 100

9 - Conclusioni e proposta di una procedura per le “prove numeriche”.

Come precedentemente spiegato, il “set-up” del problema e le scelte riguardo agli aspetti “numerici” sono stati guidati dal know how acquisito in casi fluido-termodinamici svolti in precedenza ed anche molto diversi dal “caso fisico” in studio. Ovviamente nell’applicazione ad un caso analogo non sarebbero necessari tutti i tentativi guida che si sono rivelati invece necessari per la “messa a punto” del problema. Le numerose simulazioni numeriche effettuate sono state ritenute necessarie per le seguenti ragioni: 1. Verificare lo scarto percentuale tra i valori misurati ed i valori calcolati; 2. Verificare l’importanza ed il “peso a livello numerico” della qualità della griglia di calcolo e della minore o maggiore discretizzazione del dominio di calcolo; 3. Verificare la sensibilità a diverse opzioni a disposizione in virtù del solutore numerico adottato; 4. Stabilire procedure generalizzate per la geometria in studio in funzione della fisica dell’efflusso alle condizioni di prova. 5. Confermare l’affidabilità del calcolo in condizioni diverse da quelle, certificate, di prova. Gli scarti assoluti delle portate calcolate dalle misurate sono evidentemente molto ridotti per le “prove ad alta pressione”. Nel caso delle prove a “bassa pressione”, come già spiegato, sussistono dei ragionevoli dubbi sulle cause dell’aumento degli scarti che aumentano per i casi a pressione più bassa. Le motivazioni addotte nel Capitolo 7 al paragrafo 7.6 e cioè: q stessa griglia affinata, corrispondente alla prova a pressione 6-3-boundadapt2e3---415029 1218000Pa, usata per tutte le prove a bassa pressione ovvero; q incertezza sulle modalità delle prove sperimentali (piattello bloccato o alzata limitata?), spiegano il fenomeno, e confermano la presenza di un margine di miglioramento che esula però dallo scopo di questo lavoro.

Per chiarezza e per completezza si riassumono ora i passi significativi per una proposta di procedura per l’esecuzione di “prove numeriche” in problemi analoghi. Specifichiamo che non riteniamo corretto prescindere dall’esecuzione di prove sperimentali, ma è ragionevole ridurne il numero. Infatti esse devono avere lo scopo di porre dei punti di riferimento per le “prove numeriche”, le quali possono “coprire” i campi o le configurazioni (caso in cui si abbia bisogno di prove sulla stessa valvola dotata di accessori quali ad esempio cappellotti parapioggia) per cui non si hanno dati sperimentali.

q

Geometria Nella fase preliminare si studia la geometria della valvola e dell’apparato sperimentale usato. Specialmente per le valvole a scarico convogliato e dove si effettuano prove con contropressione (backpressure) impressa, risulta fondamentale la conoscenza delle geometrie (e delle condizioni fisiche) a valle del dispositivo. Si è confermata l’importanza della fedeltà del disegno delle geometrie dei particolari specialmente nell’interno valvola e comunque nelle zone “critiche” (vedi passaggio dalla geometria con otturatore piatto all’otturatore con guarnizione, bordino e camera anulare). Si raccomanda inoltre un rilievo dal vero del modello reale provato (ove possibile). Ciò per assicurare la corrispondenza tra disegno costruttivo e valvola e verificare le tolleranze di lavorazione.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 101

9 - Conclusioni e proposta di una procedura per le “prove numeriche”.

q

q

q

q

q q

Ricerca di “periodicità” e/o simmetrie Si deve valutare la possibilità (come avvenuto nel caso in studio) di poter “simulare” una porzione limitata del dispositivo (e di ciò che lo circonda). Infatti la presenza di “periodicità” e/o “simmetrie” può ridurre considerevolmente il volume fisico da simulare. Ciò comporta vantaggi per l’ottimizzazione delle risorse di calcolo e soprattutto per l’accuratezza della soluzione. Infatti, come esposto nel Capitolo 5, si possono destinare le risorse di calcolo a disposizione ad un volume limitato, e quindi operare le discretizzazione del “continuo” con maggiore libertà e possibilità. Ciò equivale a poter destinare un maggior numero di nodi (e celle) per le zone critiche, così com’esposto nei capitoli precedenti. Inoltre, si ha anche il vantaggio di poter ottenere una griglia finale in cui il rapporto tra i volumi degli elementi più piccoli e quelli più grandi (in genere destinati alle zone più lontane dalla zona critica dei getti), risulterà meno ridotto a tutto vantaggio della perdita di accuratezza per problemi “numerici”. Le superfici di contorno Si stabiliscono i limiti (i confini) del dominio di calcolo. Come visto in precedenza la superficie di contorno ortogonale al getto è stata posta a circa cento diametri dall’asse valvola. In generale la distanza dalla valvola delle superfici di contorno, corrispondenti a sezioni del dominio fluido, dipende dalle caratteristiche del campo fluido-termodinamico atteso. Queste superfici si pongono alle distanze per cui non sussistono alti gradienti delle grandezze fisiche. Quindi, esse possono essere poste anche in prossimità della valvola (come nel caso della superficie di contorno superiore orizzontale del caso studiato). La griglia di calcolo Come esposto nel Capitolo 5 e come indicato nei Manuali del Fluent, esistono molteplici possibilità per questa operazione, a seconda della scelta della geometria degli elementi e dalla “filosofia” adottata per la costruzione della griglia a livello macroscopico. Volendo seguire l’esempio della valvola in studio, si deve dividere il dominio di calcolo in sotto domini da grigliare. Si deve prestare attenzione però alla necessità di “matching”, cioè di corrispondenza tra le facce dei sotto domini adiacenti una volta grigliati. Le posizioni dei sotto domini devono essere scelte nell’ottica della successiva operazione di “meshing” degli stessi. Essi devono aiutare ad aggirare i problemi derivanti da geometrie complesse nelle tre dimensioni e ad ottimizzare la griglia nelle zone in questione (in special modo in prossimità delle pareti ed eventualmente con l’ausilio dei boundary layers). Si genera quindi la “griglia di base”: griglia multiblocco con uso di griglie strutturate e non strutturate ed eventuali “non conformal grids” con griglie di interfaccia. Impostazioni del calcolo Simulazioni in stato stazionario (il campo di moto e termico sono indipendenti dal tempo). Le proprietà caratteristiche del mezzo fluido sono considerate quindi come funzioni dello spazio nel sistema di riferimento: si attiva l’opzione “steady”. Geometria fissa 3D (tutte le entità geometriche non cambiano di posizione e/o forma nel tempo). Solutore: accoppiato, esplicito con discretizzazione “first order upwind”. Tale solutore è adatto alla risoluzione di flussi compressibili e ad alta velocità. Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 102

9 - Conclusioni e proposta di una procedura per le “prove numeriche”.

q

q

q

Definizione del fluido evoluente. Si tiene conto della sua compressibilità e dei fenomeni ad essa collegati (in generale si attinge al database del Fluent, altrimenti si introducono le caratteristiche fisiche del fluido aggiungendole alle librerie esistenti). Pressione operativa: Fluent permette di stabilire una pressione detta “operativa”. La relazione che lega pressione assoluta pabs, pressione operativa po e pressione relativa prel è: p abs = p o + p rel Tutte le pressioni inserite e calcolate nelle simulazioni in Fluent sono in termini di pressione relativa prel. Nel caso in studio la pressione operative è stata posta coincidente con quella “ambiente”. Modello di turbolenza: k-ε RNG (a due equazioni) con opzione attivata per “swirl dominated flow”. Condizioni al contorno Le condizioni al contorno assegnate ai confini del dominio fluido sono le seguenti: • “Pressure-Inlet”: nelle sezioni “fluide” all’interno del serbatoio o comunque nelle sezioni di ingresso. Si devono specificare: 1. la sovrapressione (gauge total pressure) rispetto alla pressione di riferimento (pressione ambiente come precedentemente spiegato); 2. la temperatura totale (coincidente con la statica per il fluido quasi in quiete nelle sezioni scelte poste a debita distanza dalla sezione di ingresso alla valvola); 3. intensità della turbolenza “turbulence intensity” che è il rapporto tra il valore medio delle fluttuazioni turbolente e la velocità media (Fluent suggerisce valori di default); 4. “turbulent viscosity ratio” per flussi all’interno di condotti questo parametro assume valori compresi tra 1 e 10. • “Pressure-Outlet”: all’uscita del dominio di calcolo, sulle sezioni di contorno del dominio fluido. Si deve specificare il valore della pressione relativa ambiente nelle sezioni di “uscita”. Nel caso in studio data la coincidenza con la pressione operativa, la pressure-outlet ha valore nullo. • “Wall”: per gli elementi del dominio che rappresentano le pareti, cioè corpo valvola (interno ed esterno), otturatore con guarnizione di tenuta. Si impongono le distribuzioni delle temperature (o dei flussi termici) sulle diverse pareti e si definisce soltanto la rugosità superficiale tenendo conto della tipologia delle lavorazioni meccaniche tradizionali (tornitura, fresatura, foratura) applicate per la costruzione. • “Simmetry”: per le sezioni dei piani di simmetria sopra indicate. FLUENT attribuisce flusso nullo a tutte le grandezze di trasporto attraverso un piano di simmetria. Non esiste cioè flusso convettivo: la componente normale (al piano) della velocità è nulla. Non esiste flusso diffusivo: i gradienti (normali al piano) di tutte le variabili sono nulli.

q

Condizioni iniziali L’inizializzazione del campo ha notevole importanza. Può essere necessario usare una soluzione corrispondente ad un valore di pressione inferiore. Nel caso in studio abbiamo inizializzato il campo operando tramite “patching”. Come discusso nel Capitolo 6, esso consiste nel suddividere idealmente il dominio in zone diverse (regions). Questa caratteristica di FLUENT permette una operazione di “patching” con successive sovrapposizioni di condizioni iniziali. Si inizia infatti con c.i. identiche per tutto il dominio fluido.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 103

9 - Conclusioni e proposta di una procedura per le “prove numeriche”.

Successivamente si effettuano le sovrapposizioni nelle zone preventivamente individuate e contrassegnate (marked) con altre c.i. che sostituiscono quelle ivi precedentemente assegnate. Si riesce in tal modo ad assegnare gradienti (ovviamente a gradini) per tutte le grandezze. I valori di primo tentativo sono stati scelti in base a simulazioni precedenti effettuate su griglie molto rade (∼100000 celle) e con il solutore “segregato”. Ciò ha comportato un calcolo preventivo con l’efflusso subsonico per via delle caratteristiche del solutore sopra detto incapace di trattare flussi supersonici o localmente supersonici. I valori poi assegnati sono stati maggiorati per tenere conto dei salti di pressione maggiori corrispondenti alle prove in condizioni di “efflusso critico”.

q

q q

q

q

q

Valutazione della convergenza e dei primi risultati Una volta effettuati i cicli iterativi e raggiunta la convergenza (perché raggiunti tutti i limiti inferiori impostati precedentemente per le incognite) si valutano i diagrammi relativi alle incognite e le modalità di decremento. Essa è stata valutata sia in base ai valori assunti dai “residui” sia in base alla storia di convergenza. Anche l’analisi del bilancio della massa entrante e uscente con il procedere del procedimento iterativo (si è valutata la sua stabilizzazione) ha permesso di stabilire la convergenza. In molti casi, la stabilizzazione dei valori assunti dalle grandezze fisiche in zone diverse e significative fornisce criterio per la valutazione della convergenza e della validità del calcolo stesso. Nel caso le storie di convergenza siano accettabili, si può passare alla successiva fase di affinamento della griglia. Nel caso non siano accettabili (per via della forma dei diagrammi non regolarmente decrescenti o perché si rivela necessario un numero di cicli troppo elevato) o addirittura non si raggiunga la convergenza, si devono ripercorrere i passi precedenti operando alcune modifiche. Ad esempio il problema può derivare da una griglia non ben costruita; anche l’inizializzazione può essere stata non efficace. Si verifica la correttezza della simulazione dal punto di vista fisico. Si controlla che sia stata riprodotta la fenomenologia dell’efflusso. Si valuta lo scarto con i valori sperimentali (laddove presenti). Se ne valuta l’accettabilità. La griglia “affinata” Come esposto nei precedenti capitoli si generano griglie ottimizzate nei primi due o tre strati vicino alle pareti e/o in funzione della y+ (“solution-adaptive refinement”: l’ottimizzazione delle griglie è effettuata sulla griglia base in funzione dei risultati ottenuti dal calcolo preliminare sulla griglia base stessa). Il calcolo definitivo In base agli affinamenti effettuati sulla griglia si ottiene una nuova soluzione che si confronterà ovviamente con la precedente. Si valuterà lo scarto con i risultati sperimentali e lo scarto percentuale rispetto alle prima simulazione. Le simulazioni fuori “campo” Le simulazioni nei campi di funzionamento per cui non si hanno prove sperimentali usano, per la loro inizializzazione, le simulazioni effettuate per valori di pressione inferiore. Per

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 104

9 - Conclusioni e proposta di una procedura per le “prove numeriche”.

condizioni molto diverse da quest’ultime può essere necessario inizializzare come esposto nei passi precedenti addirittura rigenerando una griglia ad hoc. Per concludere, possiamo affermare che lo studio effettuato risponde alle aspettative prefissate ed evidenzia la validità dell’uso “intelligente” e sistematico delle simulazioni CFD (che costituiscono d’altronde uno strumento di calcolo oramai diffusamente utilizzato in molteplici campi dell’ingegneria) per la progettazione e l’ottimizzazione di una generica attività di Progettazione Ingegneristica sotto il profilo tecnico ed economico. Qui, in particolare, si è dimostrato come l’applicazione sistematica di una procedura di “prove numeriche” su una geometria di valvola di sovrapressione possa venire integrata con una molto ristretta serie di prove sperimentali per fornire un risultato complessivo equivalente, sia in termini di affidabilità che di completezza, a quello ottenibile mediante una serie molto più estesa –e quindi molto più costosa in termini di tempo e di risorse- di prove sperimentali. Le prove numeriche hanno fra l’altro il vantaggio di poter essere usate per incrementare la conoscenza specifica del progettista, ponendolo di fronte a fenomeni che siano sì a lui noti, ma che vengono quantizzati in modo nuovo dal codice numerico, o addirittura a fenomeni imprevisti che possano giustificare certe caratteristiche delle prestazioni della valvola. In questa ottica, si può senz’altro affermare che la metodologia è valida sia nel particolare caso esaminato che in generale, e meriterebbe di essere estesa ad altri tipi di verifiche di collaudo.

RINGRAZIAMENTI Desideriamo ringraziare gli ingegneri Corrado Delle Site, Fausto Di Tosto ed Emanuela Franchi per la disponibilità e la cordiale collaborazione con cui ci hanno facilitato enormemente lo svolgimento del lavoro qui presentato.

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 16/06/2005

Pagina 105

Bibliografia

RIFERIMENTI BIBLIOGRAFICI

Celik, I.B. Introductory Turbulence Modeling (1999)

Cunsolo, D. Dispense delle lezioni del Corso di Aerodinamica. Università degli Studi di Roma “La Sapienza”.

FLUENT 5 User’s guide (1998)

Mattiussi, C. The Finite Volume, Finite Element and Finite Difference Methods as Numerical Methods for Physical Fields Problems (2000) – Swiss Federal Institute of Technology CH1015 Lausanne Switzerland Millikan, C. B. A critical discussion of turbulent flows in channel and circular tubes (1938) – Proc. 5th Int. Congr. Applied Mechanics, New York, pp386-392

Peyret, R.; Taylor, T.D., Computational Methods for Fluid Flow. Springer-Verlag 1984.

Pope, S.B. Turbulent Flows (2000) - Cambridge University Press

Quori, F., Aerodinamica (1998)

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 106

Indice

INDICE 1.

Descrizione del progetto e piano delle attività........................................................... 2 1.1

Descrizione del progetto................................................................................................. 2

1.2

Piano delle attività e scadenze........................................................................................ 2

2

Scelta delle geometrie rappresentative di apparati reali (“geometrie tipo”). Dati

sperimentali.................................................................................................................................. 3 2.1

Geometrie tipo................................................................................................................ 3

2.2

Dati sperimentali............................................................................................................ 8

3

Le equazioni di governo del moto di un fluido.......................................................... 10 3.1

Introduzione................................................................................................................... 10

3.2

Equazione di continuità o di conservazione della massa............................................... 11

3.3

Equazione della conservazione della quantità di moto.................................................. 13

3.4

Equazione dell’energia.................................................................................................. 16

3.5

Fluidi incompressibili.................................................................................................... 19

3.6

Fluidi compressibili....................................................................................................... 19

3.7

Equazione di stato.......................................................................................................... 19

3.8

Fluidi incompressibili.................................................................................................... 20

3.9

Fluidi compressibili....................................................................................................... 20

3.10 Condizioni al contorno e condizioni iniziali................................................................. 21 4

La fisica del fenomeno dell’efflusso in una valvola.................................................. 22 4.1

Generalità...................................................................................................................... 22

4.2

La valvola vista come un ugello................................................................................... 23

5

Realizzazione delle meshes “caratteristiche”: la discretizzazione del

continuo........…………………………………………………………………………………... 27 5.1

Dominio di calcolo: come isolare una porzione del sistema fisico completo............... 27

5.2

Superfici di contorno: pareti e sezioni del volume fluido completo............................. 27

5.3

Individuazione di piani di simmetria............................................................................. 29

5.4

Definizione delle superfici di contorno......................................................................... 32 5.4.1 Sezioni “fluide”.................................................................................................... 32 5.4.2 Pareti solide.......................................................................................................... 34

5.5

Generazione della griglia di calcolo.............................................................................. 34

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 107

Indice

5.6

Tipologia di griglia impiegata....................................................................................... 34

5.7

Brevi cenni descrittivi sul “Cooper scheme”................................................................ 41

6

Le simulazioni numeriche. Generalità...................................................................... 42 6.1

Impostazioni del calcolo............................................................................................... 43

6.2

Cenni sulla turbolenza................................................................................................... 44

6.3

Effetto delle pareti e significato del parametro dimensionale y+................................. 47 6.3.1 Il sottostrato laminare.......................................................................................... 49 6.3.2 Lo strato di ‘sovrapposizione’............................................................................. 49 6.3.3 Lo strato esterno.................................................................................................. 50

6.4

Condizioni al contorno................................................................................................. 51

6.5

Condizioni iniziali........................................................................................................ 53

7

Il procedimento di calcolo.......................................................................................... 57 7.1

Le simulazioni preliminari............................................................................................ 57 7.1.1 Fase-1 - Caso 11-1---174010.............................................................................. 59 7.1.2 Fase-1 - Caso 5-1---147917................................................................................ 60 7.1.3 Fase-1 - Caso 5-2---160896................................................................................ 61

7.2

Fase-2 - Affinamento della geometria del modello..................................................... 68

7.3

Simulazioni definitive. Analisi numerica e fisica dei risultati..................................... 74

7.4

Affinamento delle griglie di calcolo............................................................................ 82

7.5

Le “prove intermedie”.................................................................................................

88

7.5.1 Prima prova intermedia: 700000Pa....................................................................

88

7.5.2 Seconda prova intermedia: 1050000Pa..............................................................

88

7.6 8

9

Le “prove a bassa pressione”....................................................................................... 89 Confronto con i risultati sperimentali...................................................................... 91

8.1

Riferimenti normativi...................................................................…............................ 91

8.2

Coefficienti di efflusso e portate reali...............................................................…....... 93 Conclusioni e proposta di una procedura per le “prove numeriche”................... 100 Bibliografia................................................................................................................. 106 Indice.......................................................................................................................... 107

Dipartimento di Meccanica ed Aeronautica Autori: P. Bernardini, E. Sciubba

Roma, 14/06/2005

Pagina 108