TERMOFLUIDODINAMICA - TERMOFLUIDODINAMICA BIFASE

Ingegneria Meccanica e Gestionale

· TERMOFLUIDODINAMICA BIFASE

· Valerio Marra

· LA TENSIONE SUPERFICIALE NELLE EQUAZIONI DI NAVIER-STOKES

· IL METODO VOF/PLIC

· FORMULAZIONE AD UN FLUIDO

· SOLUZIONE NUMERICA DELLE EQUAZIONI NELL’APPROSSIMAZIONE DI BOUSSINESQ PER FLUSSI BIFASE

· RISULTATI NUMERICI: INSTABILITà DI RAYLEIGH-TAYLOR

SOMMARIO

Nello studio dei flussi bifase è di cruciale importanza la descrizione e la dinamica dell’evoluzione dell’interfaccia materiale di separazione. La presenza di tale elemento comporta: l’aggiunta di un termine nell’equazione della quantità di moto che tiene in conto l’azione della tensione superficiale e l’introduzione di nuove equazioni e tecniche numeriche per il tracking dell’interfaccia.

Nel seguito verranno descritte le modifiche che saranno apportate al sistema di equazioni già presentato in Termofluidodinamica Computazionale, Ing. V. Marra al fine di soddisfare le esigenze sopra citate. Inoltre sarà fornita una breve introduzione ad uno dei metodi maggiormente affermati e sviluppati dai ricercatori operanti nel campo della fluidodinamica bifase: il Metodo VOF/PLIC. Si veda l’articolo Metodo VOF/PLIC per Flussi Bifase, Ing. V. Marra.

LA TENSIONE SUPERFICIALE NELLE EQUAZIONI DI NAVIER-STOKES

Consideriamo due fluidi immiscibili a contatto e in quiete e prendiamo un punto sull'interfaccia di separazione. Tracciamo il piano tangente per questo punto e due piani perpendicolari tra di loro e anche al piano tangente. Questi due piani intersecano l'interfaccia lungo due linee, siano R₁ e R₂ il loro raggio di curvatura nel punto in questione (per convenzione attribuiamo il segno positivo se il raggio di curvatura è diretto verso il fluido di riferimento 1). Si può dimostrare (Teorema Egregium di Gauss) che la somma 1/R₁+1/R₂ è un invariante rispetto alla scelta dei due piani ortogonali. In particolare, quando uno dei due raggi di curvatura è massimo l'altro sarà minimo: tali valori vengono detti raggi principali di curvatura. L'equilibrio soddisfa l'equazione di Laplace:

(1)

dove s è la tensione superficiale che dipende dalla natura dei fluidi e dallo stato termodinamico del sistema (temperatura). In condizioni dinamiche con tensione superficiale costante esiste una forza normale all'interfaccia data dall'espressione:

(2)

dove k=(1/R₁+1/R₂)/2 è la curvatura media, n è il versore normale alla superficie e:

(3)

è la funzione delta di Kronecker. Si può dimostrare che il termine singolare (2) può essere scritto come la divergenza di un tensore G detto tensore capillare:

(4)

dove G=-s(I-nÄn)d_S e ha la proprietà di essere tangente all'interfaccia in ogni suo punto. L'equazione di Navier-Stokes con il contributo della tensione superficiale diventa (cfr. equazione (2) in Metodi Numerici per le Equazioni di Navier-Stokes, Ing. V. Marra):

(5)

IL METODO VOF/PLIC

Il metodo descritto, in una sua versione semplificata, per il tracking dell’interfaccia di separazione è il VOF/PLIC (Volume Of Fluid/Piecewise Linear Interface Calculation). Nel caso di due fluidi immiscibili, la loro posizione è specificata per mezzo della funzione frazione di volume C che assume i seguenti valori:

(6)

Numericamente la funzione C rappresenta la frazione volumetrica di ogni singola cella del grigliato che è occupata dal fluido di riferimento (Figura 1a). In questo metodo, assegnata la frazione di volume C della fase di riferimento e una stima del vettore n normale all’interfaccia in ogni cella computazionale (Figura 1b), la generica curva rappresentante l’interfaccia di separazione è ricostruita all’interno della cella mediante un segmento di retta orientato secondo n che la taglia in due parti, le cui rispettive frazioni di volume corrispondono per il fluido di riferimento a C e per quello complementare a 1-C.

Figura 1a Figura 1b

Il vettore normale n viene legato a C nel seguente modo:

(7)

In assenza di cambiamenti di fase C segue semplicemente il moto del fluido (i.e. si comporta come uno scalare passivo) e quindi soddisfa la seguente equazione di avvezione:

(8)

dove si anche può pensare alla (8) anche come all’equazione di conservazione della fase di riferimento 1. Nel caso in cui le variazioni di C siano “brusche” la (8) va intesa in senso “debole”.

I vari coefficienti coinvolti nel sistema di equazioni da risolvere per descrivere il moto del fluido sono assegnati alle singole fasi e quindi possono essere espressi in funzione di C. Considerando ad esempio la densità si può scrivere:

(9)

dove r₁ e r₂ sono, rispettivamente, la densità della fase 1 e della fase 2.

FORMULAZIONE AD UN FLUIDO

La soluzione delle equazioni della termofluidodinamica nell’approssimazione di Boussinesq per un sistema bifase si ottiene risolvendo per ciascuna fase il seguente sistema equazioni (cfr. per la loro derivazione e descrizione l’articolo Termofluidodinamica Computazionale, Ing. V. Marra):

(10)

(11)

(12)

corredato delle opportune condizioni al contorno e iniziali. Inoltre le due soluzioni ottenute per ciascuna fase devono soddisfare delle opportune condizioni di “saldatura” in corrispondenza dell’interfaccia di separazione materiale. Ciò comporta l’imposizione di ulteriori vincoli da soddisfare che complicano ulteriormente il modello sia dal punto di vista teorico che computazionale.

L’utilizzo del metodo VOF/PLIC che comporta una nuova definizione delle diverse grandezze materiali presenti in (13) e (14), secondo lo schema presentato in (11), e la soluzione delle equazioni aggiuntive (9) e (10), rispettivamente, per n e C consente la soluzione del seguente unico set di equazioni:

(10)

(11’)

(12’)

(13)

(7)

(8)

dove, si può dimostrare, le condizioni di “saldatura” sopra menzionate sono implicitamente soddisfatte. In tale formulazione ad un fluido il campo di velocità u¢ è condiviso dalle fasi: laddove il valore di C è pari a 1 oppure a 0 si riduce al classico sistema di equazioni della termofluidodinamica, rispettivamente, per la fase di riferimento e quella complementare. Per valori 0<C<1 (i.e. nelle celle tagliate dall’interfaccia) vengono calcolati: un valore opportunamente mediato (sulla cella computazionale) delle proprietà materiali e il contributo dell’azione della tensione superficiale nel bilancio della quantità di moto. Tali valori contribuiscono alla determinazione del campo di velocità in presenza dell’interfaccia.

SOLUZIONE NUMERICA DELLE EQUAZIONI DELLA TERMOFLUIDODINAMICA NELL’APPROSSIMAZIONE DI BOUSSINESQ PER FLUSSI BIFASE

Lo schema numerico di soluzione relativo alle equazioni (10), (11¢), (12¢), (13), (7) e (8), si esplica nei seguenti quattro passi:

1. inizializzazione dei campi u¢, p, C e T; calcolo della (13) e della (7).

2. soluzione della (10) e della (11¢), considerando il campo di temperatura non come incognito ma come assegnato, secondo il metodo descritto in Metodi Numerici per le Equazioni di Navier-Stokes, Ing. V. Marra. Il campo vettoriale u¢ è aggiornato.

3. soluzione della (12¢) secondo uno schema numerico, anch’esso alle differenze finite, nel quale il campo di velocità è considerato non come incognito ma come assegnato dal passo di soluzione precedente (cfr. Termofluidodinamica Computazionale, Ing. V. Marra). Il campo scalare T è aggiornato.

4. soluzione della (8) con u¢ considerato non come incognito ma come assegnato dal passo 2. Il campo scalare C è aggiornato.

5. i passi 2., 3. e 4. sono ripetuti fino al soddisfacimento di un assegnato criterio di convergenza.

Le variabili incognite temperatura e frazione di volume sono, in questa approssimazione, scalari passivi trasportati dal fluido in movimento (i.e. non esercitano direttamente sulle equazioni di conservazione della massa e della quantità di moto alcuna influenza, la loro evoluzione dinamica è quindi completamente determinata da esse). Al fine di catturare questa peculiarità è indicato anche in questo caso l’uso di una tecnica upwind per il trattamento dei termini convettivi (ad eccezione che per la (8), oggetto di tecniche speciali). In questo schema di soluzione numerica, nel quale per la cattura della dinamica del campo di velocità è già utilizzato il metodo di Godunov, è sufficiente che la tecnica upwind utilizzata per la soluzione dell’equazione (12¢) sia un metodo upwind del primo ordine. Ciò consente una determinazione accurata del campo di temperatura evitando l’ulteriore complicazione dell’algoritmo di soluzione e un’aumento dello sforzo computazionale richiesto.

La pressione, la frazione volumetrica e le componenti della velocità sono discretizzate su una griglia cartesiana uniforme (Dx=Dy=h), la distribuzione delle variabili sulle celle è di tipo MAC con le componenti della velocità definite sui bordi delle celle (Figura 2). I volumi di controllo W_u e W_v delle componenti della velocità sono mostrati in Figura 2 e sono utilizzati per la risoluzione del sistema di equazioni della termofluiodinamica nell’approssimazione di Boussinesq. Anche se non rappresentata, la discretizzazione e la distribuzione della variabile temperatura T è identica a quella adottata per la pressione e la frazione volumetrica.

Figura 2

RISULTATI NUMERICI: INSTABILITà DI RAYLEIGH-TAYLOR

In natura si presentano diverse situazioni fisiche legate alla stabilità dell'interfaccia di separazione tra mezzi diversi. Una di queste è quella in cui un fluido pesante è sospeso su di uno più leggero in presenza di un campo gravitazionale diretto verso quest'ultimo. Una perturbazione sinusoidale della superficie viene amplificata e quando entra in un regime altamente non lineare, l'interfaccia assume la ben nota forma di un fungo rovesciato (Figura 3). Questa è conosciuta come instabilità di Rayleigh-Taylor.

Figura 3

Questo tipo di instabilità può essere avviata anche nell’ambito di uno stesso fluido, mantenendo ad esempio un gradiente di temperatura riscaldandone la parte inferiore. Quest’ultima diviene più “leggera” della parte di fluido sovrastante: questo è una disposizione “top-heavy” potenzialmente instabile poiché, nonostante la viscosità, il fluido mostra la tendenza a ridistribuirsi quando il gradiente di temperatura supera un certo valore critico. Il parametro adimensionale che determina per data geometria e tipo di fluido se il gradiente di temperatura mantenuto è in grado di innescare o meno l’instabilità di Rayleigh-Taylor è noto in letteratura come numero di Rayleigh R.

Nelle figure successive è mostrato lo sviluppo dell’instabilità di Rayleigh-Taylor innescata da tale gradiente di temperatura (Figura 4).

Figura 4