Cap. 10. Modelli di turbolenza per equazioni mediate alla Reynolds
by user
Comments
Transcript
Cap. 10. Modelli di turbolenza per equazioni mediate alla Reynolds
Cap. 10. Modelli di turbolenza per equazioni mediate alla Reynolds 10.1 Introduzione Richiamiamo le equazioni di Navier-Stokes che descrivono una corrente fluida: ∂ρ ∂(ρvj ) + =0 ∂t ∂xj ∂p ∂ ∂(ρvi ) ∂(ρvi vj ) + =− + ∂t ∂xj ∂xi ∂xi (10.1) ∂vi µ ∂xj (10.2) Dalla teoria di Kolmogorov (vedi capitolo 5) sappiamo che il rapporto tra la scala dissipativa ld e la scala integrale L è dell’ordine 3 ld ≃ Re− 4 . L Pertanto una simulazione numerica completa di una corrente turbolenta richiederebbe una 9 griglia con N 3 ≃ Re 4 punti di collocazione, cosa attualmente non affrontabile per elevati numeri di Reynolds. Una soluzione del problema è stata proposta da Sir Osborne Reynolds (1895) applicando l’operazione di media alle equazioni di Navier-Stokes. Secondo questo approccio una quantità f viene decomposta in una componente media hf i e una fluttuante f ′ f = hf i + f ′ . La media alla Reynolds può assumere diversi significati, a seconda del tipo di corrente a cui viene applicata: • turbolenza omogenea hf i = hf (t)i 1 V →∞ V hf i = lim f (x, t)d3 x V Media spaziale • turbolenza statisticamente stazionaria 1 hf i = lim T →∞ T hf i = hf (x)i Z Media temporale 107 Z 0 T f (x, t)dt CAPITOLO 10. MODELLI DI TURBOLENZA PER EQUAZIONI MEDIATE 108 ALLA REYNOLDS • turbolenza generica N 1 X f(α) (x, t) N →∞ N α=1 hf i = lim hf i = hf (x, t)i Media di insieme su N eventi. In fluidodinamica viene comunemente accettata l’estensione del teorema di ergodicità per cui le medie temporali e spaziali possono essere ritenute medie d’insieme. Date f = hf i + f ′ e g = hgi + g ′ valgono le seguenti regole ′ ′ f = g =0 hhf ii = hf i hf i f ′ = hgi g′ = hf i g′ = hgi f ′ = 0 f 2 = hf i2 + f ′2 hf gi = hf i hgi + f ′ g ′ ∂ hf i ∂f = ∂xi ∂xi ∂f ∂ hf i = ∂t ∂t Applichiamo alle equazioni di Navier-Stokes la decomposizione alla Reynolds al campo di velocità e pressione: vi = Ui + vi′ p = hpi + p′ . dove abbiamo posto hvi i = Ui . Si sottolinea che nel caso di fluido a densità ρ costante hρvi i = ρUi mentre nel caso di fluido a densità variabile viene utilizzata una media alla Favre hρvi i = hρi hρvi i = hρi ṽi hρi Nel seguito considereremo ρ costante. Mediando il termine convettivo hρvi vj i = ρ (Ui + vi′ )(Uj + vj′ ) = ρ(Ui Uj + vi′ vj′ ) e definendo τijR = − vi′ vj′ 10.1. INTRODUZIONE 109 come le componenti del tensore degli sforzi di Reynolds, otteniamo le equazioni RANS (Reynolds-averaged Navier-Stokes) per il moto medio ∂Uj =0 ∂xj (10.3) ∂Ui ∂(Ui Uj ) 1 ∂ hpi ∂ + =− + ∂t ∂xj ρ ∂xi ∂xj ∂Ui ν ∂xj + ∂τijR (10.4) ∂xj dove ν = µ/ρ è la viscosità cinematica. Per la presenza di 6 nuove incognite rappresentate dalle componenti del tensore degli sforzi di Reynolds, le equazioni (10.3,10.4) non cosituiscono un sistema chiuso che permetta di determinare i campi di velocità hvi e pressione hpi medie. Nelle equazioni RANS è quindi necessaria l’assunzione di un modello che leghi in modo fisicamente consistente il tensore degli sforzi di Reynolds alla storia globale del campo di velocità medio τijR = Fij [ v(x′ , t′ ) ; x, t], x′ ∈ V, t′ ∈ (−∞, t). (10.5) Al fine di cercare di risolvere il problema della chiusura delle equazioni RANS, ricaviamo ora delle equazioni per il moto turbolento. Sottraendo alla (10.2) la (10.4) si ottiene l’equazione per la componente fluttuante del campo di velocità ∂v ′ ∂Ui ∂v ′ 1 ∂p′ ∂ ∂vi′ + Uj i + vj′ i + vj′ =− + ∂t ∂xj ∂xj ∂xj ρ ∂xi ∂xj ν ∂vi′ ∂xj − ∂τijR ∂xj (10.6) e analogamente si ottiene l’equazione di continuità per le fluttuazioni ∂vj′ =0 ∂xj (10.7) Si faccia la somma dei prodotti della j-esima componente della (10.6) moltiplicata per vi′ e della i-esima componente della (10.6) moltiplicata per vj′ , e si operi la media; si ottiene cosı̀ l’equazione di trasporto per gli sforzi di Reynolds: ∂τijR ∂τijR + Uk ∂t ∂xk ′ ∂vj′ ∂vi ∂Uj ∂Ui ′ = −τik − τjk − p + + ∂xk ∂xk ∂xj ∂xi ′ ′ ′ ′ 1 ′ ′ ∂vi ∂vj′ ∂ 1 ′ ′ 2ν vi vj vk + p vi δjk + p vj δik + + ∂xk ∂xk ∂xk ρ ρ R ∂τjk ∂τ R (10.8) −vj′ ik − vi′ ∂xk ∂xk Nell’equazione (10.8) sono comparsi numerosi termini derivati dalla media di prodotti di quantità fluttuanti che costituiscono 22 nuove incognite. Pertanto la soluzione delle equazioni per gli sforzi di Reynolds rende ancora più complesso il problema della chiusura. Proprio a causa della non-linearità delle equazioni di Navier-Stokes, scrivere equazioni per momenti di ordine sempre piú elevato introduce sempre più incognite senza poter raggiungere un bilancio con le equazioni. CAPITOLO 10. MODELLI DI TURBOLENZA PER EQUAZIONI MEDIATE 110 ALLA REYNOLDS 10.2 Modelli a viscosità turbolenta Per risolvere il problema della chiusura, Boussinesq (1877) per primo propose un modello per gli sforzi turbolenti nell’ipotesi di turbolenza sviluppata ad alti numeri di Reynolds, dove convezione e diffusione sono trascurabili, mentre produzione e dissipazione degli sforzi di Reynolds sono in equilibrio. Boussinesq assunse l’ipotesi di similarità tra il moto delle strutture turbolente e quello molecolare. Per meglio comprendere le motivazioni di Boussinesq analizziamo il moto molecolare. y U(y) Q 2l m t xy t xy x P Figura 10.1: Consideriamo un flusso bidimensionale come in figura (10.1) in cui la velocità macroscopica è U = U (y)i dove i è il versore dell’asse x. Ricordando che il moto molecolare è casuale, decomponiamo la velocità delle molecole in una componente macroscopica U corrispondente al profilo di figura e una componente u′′ corrispondente al moto casuale u = U + u′′ Consideriamo ora il flusso di quantità di moto longitudinale specifica Qx = ρu attraverso una superficie infinitesima dS del piano y = 0 dfx = ρuvdS = ρ(U + u′′ )v ′′ dS dove v è la componente della velocità normale alla parete (la velocità normale media viene supposta nulla). Effettuando la media d’insieme si ottiene dFy = ρ u′′ v ′′ dS Per la definizione di sforzo secondo la convenzione della normale entrante 10.2. MODELLI A VISCOSITÀ TURBOLENTA 111 dFy dS e considerando la parte anisotropa del tensore degli sforzi τij = σij − pδij risulta σxy = − τxy = −ρ u′′ v ′′ È evidente la somiglianza tra questo tensore ed il tensore degli sforzi di Reynolds τijR (10.9) D E = − u′i u′j dove le fluttuazioni turbolente di velocità u′ a livello macroscopico sostituiscono le fluttuazioni di velocità molecolare u′′ . y φ vt x Figura 10.2: Nell’ipotesi di gas perfetto, le molecole si muovono con uguale probabilità in tutte le direzioni (distribuzione Maxwelliana) con una velocità di agitazione termica vt pari circa a 4/3 la velocità del suono nell’aria. La componente secondo y della velocità molecolare è vt cosφ (vedi fig.10.2) e integrando su una semisfera si ottiene la sua media nella direzione y positiva R 2π R π/2 1 dφvt cos φ sin φ = vt /2. Quindi se n è il numero di molecole per unità di vm = 2π 0 dψ 0 volume, nvt /4 è la velocità media per unità di area delle molecole che si muovono in direzione y positiva (metà delle molecole andranno verso l’alto, metà verso il basso). Una molecola di massa m partendo dal punto P , attraversando il piano y = 0 acquista una differenza di quantità di moto ∆q = m (U (0) − U (−lm )) poichè u′′ ≪ U il flusso di quantità di moto risulta 1 1 dU ∆Q+ = ρvt (U (0)) − U (−lm )) ≃ ρvt lm 4 4 dy dove ρ = nm, e U (−lm ) è stato sostituito dall’espansione in serie di Taylor al 1o ordine U (−lm ) = U (0) − lm dU dy Allo stesso modo il flusso di quantità di moto verso y negative risulta ∆Q− = 1 dU 1 ρvt (U (0) − U (lm )) ≃ − ρvt lm 4 4 dy Essendo lo sforzo τxy legato al flusso totale di quantità di moto CAPITOLO 10. MODELLI DI TURBOLENZA PER EQUAZIONI MEDIATE 112 ALLA REYNOLDS τxy = ∆Q+ − ∆Q− = dU dU 1 ρvt lm =µ 2 dy dy dove si definisce la viscosità molecolare 1 µ = ρvt lm 2 Seguendo quindi l’assunzione di Boussinesq possiamo studiare il moto delle strutture turbolente in modo analogo a quanto fatto per il moto molecolare, riuscendo cosı̀ a legare il tensore degli sforzi τij al moto medio U . Sostituiamo alle quantità molecolari quelle turbolente e sottolineiamo le ipotesi di validità: • le scale caratteristiche della turbolenza devono essere molto più piccole delle scale del moto medio. Poniamo | dU/ dy| L= 2 | d U/ dy 2 | la lunghezza a grande scala, questa condizione si traduce nel fatto che il numero di Knudsen per la turbolenza deve essere molto piccolo: Kn = lm ≪1 L • la turbolenza deve essere isotropa (distribuzione Maxwelliana) anche in presenza di shear. L’ipotesi di alta collisionalità vale se tcoll. ≪ tevol. dove tcoll. ≃ lm /vt è il tempo caratteristico delle collisioni e tevol. = | dU/1 dy| è il tempo caratteristico di evoluzione del moto medio. Introdotto il numero di Mach Ma = L| dU/ dy| vt si ottiene la condizione M a Kn ≪ 1 Arriviamo cosı̀ ad un’espressione per gli sforzi turbolenti di taglio d hui dy e generalizzando per la parte anisotropa del tensore degli sforzi turbolenti ∂Ui ∂Uj 1 + τij − τkk δij = νT 3 ∂xj ∂xi τxy = νT (10.10) dove νT è la viscosità turbolenta, caratteristica locale della corrente turbolenta. Un altro ragionamento per stabilire la validità dell’approssimazione di viscosità turbolenta è il seguente: la relazione (10.10) può essere vista come uno sviluppo del tensore τij in serie di potenze del tensore velocità di deformazione troncato ai termini lineari. Essa è valida a condizione che il rapporto tra la parte lineare (anisotropa) e quella costante (isotropa) sia piccolo, in modo che lo sviluppo troncato al primo ordine sa arrivato a convergenza e non ci sia bisogno di tenere in considerazione i termini quadratici. Nel caso molecolare |2νSij | λvth U/L λ U ∼ = KnM a ≪ 1 = 2 P/ρ L vth vth 10.3. MODELLI A 0 EQUAZIONI 113 dove λ è il libero cammino medio. Questo numero è generalmente molto piccolo (10−10 nel caso dell’aria). Allora l’espansione della parte deviatorica del tensore degli sforzi in funzione dei gradienti di velocità può fermarsi al termine lineare e il termine lineare più generale possibile compatibile con le simmetrie è 2νSij . In un fluido turbolento il rapporto vale |τija | 2 3k che in situazioni realistiche spesso non è trascurabile, anzi raggiunge valori dell’ordine di 0.5. Ciò mette in luce il fatto che l’ipotesi di viscosità turbolenta per un fluido turbolento è un’approssimazione abbastanza grossolana. Per modellare la viscosità turbolenta sono state presentate numerose proposte. Nei modelli a 0 equazioni, la viscosità turbolenta viene modellata sulla base di assunzioni empiriche. In questo modo vengono risolte le equazioni RANS per il solo campo mediato. Nei modelli a 1 o 2 equazioni invece, la viscosità turbolenta viene determinata in funzione di quantità turbolente (per es. l’energia turbolenta) per le quali è necessario risolvere delle equazioni di trasporto-diffusione. In questo modo la viscosità turbolenta assume il compito di trasferire informazioni dal campo turbolento a quello medio. Sottolineiamo che l’ipotesi di Boussinesq impone un allineamento tra il tensore degli sforzi di Reynolds ed il tensore velocità di deformazione. Quindi i modelli basati su questa ipotesi cadono in difetto in correnti • con rapide variazioni del tensore velocità di deformazione medio, • su superfici curve, • in condotti con flussi secondari, • rotatorie, • tridimensionali, • con separazione dello strato limite. 10.3 Modelli a 0 equazioni y u(y) x Figura 10.3: CAPITOLO 10. MODELLI DI TURBOLENZA PER EQUAZIONI MEDIATE 114 ALLA REYNOLDS Il primo contributo significativo è stato quello di Prandtl (1925) con il modello detto mixing length: in uno strato limite 2D, su parete con curvatura trascurabile (fig.10.3), si può ritenere che i gradienti longitudinali del moto medio siano associati principalmente al gradiente di pressione nella medesima direzione e non agli sforzi di Reynolds; i gradienti di velocità media in direzione longitudinale sono trascurabili rispetto a quelli normali alla parete, quindi gli sforzi di Reynolds si riducono alla sola componente ∂ hui − u′ v ′ = ν T ∂y (10.11) In base a considerazioni dimensionali νT = lV dove l e V sono rispettivamente una lunghezza ed una velocità tipiche del moto turbolento. Lo shear tende ad aumentare il mescolamento turbolento e quindi è ragionevole assumere una velocità caratteristica proporzionale al gradiente di velocità: ∂ hui V = lm ∂y e νT = 2 ∂ hui lm ∂y dove lm è la mixing length, che non è una caratteristica del fluido ma dipenderà dal tipo di corrente. Pertanto la sua determinazione, di origine empirica, non è universale. Vediamo qualche esempio. Lo strato limite In uno strato limite turbolento si possono distinguere tre regioni: lo strato viscoso, la regione logaritmica e la regione di scia. Dal punto di vista matematico esistono due regioni, il substrato viscoso e la zona di scia, mentre lo strato logaritmico è la regione in cui le due soluzioni si sovrappongono. Il modello mixing length con lunghezza di mixing lm = κy, dove κ è la costante di von Kàrmàn, è capace di riprodurre lo strato logaritmico: consideriamo le equazioni per il moto medio in uno strato limite incomprimibile a pressione costante ∂ hui ∂ hvi + =0 ∂x ∂y ∂ hui ∂ ∂ hui ∂ hui ′ ′ + hvi = − uv hui ν ∂x ∂y ∂y ∂y Nello strato logaritmico i termini convettivi sono trascurabili, quindi lo sforzo totale può essere considerato costante ∂ hui ′ ′ ∂ hui ν − u v = cost ≃ ν = τw = u2τ ∂y ∂y w Inoltre, poichè nello strato logaritmico gli sforzi viscosi sono trascurabili rispetto a quelli turbolenti, si ha: τxy = − u′ v ′ ≃ u2τ ρ 10.3. MODELLI A 0 EQUAZIONI 115 Applicando il modello mixing length otteniamo 2 lm ∂ hui ∂y 2 ≃ u2τ (10.12) Assumendo l’ipotesi di Prandtl lm = κy l’equazione (10.12) fornisce: hu(y)i = uτ ln y + costante κ Sostituendo le quantità adimensionali u+ = u , uτ y+ = uτ y ν si ottiene la classica legge di parete + 1 u = ln y + + C κ dove la costante di von Kàrmàn viene assunta pari a κ ≃ 0.41, e C ≃ 5. Si è cosı̀ verificato che il modello mixing length è in accordo con la legge di parete. Il sottostrato viscoso non è in regime turbolento, e il modello di Prandtl non è direttamente applicabile in esso. Van Driest (1956) propose una modifica per y → 0 valida anche nel sottostrato laminare i h + + (10.13) lm = κy 1 − e−y /Ao dove la costante A+ o = 26. Notiamo come vicino alla parete lo sforzo di taglio abbia un 2 (∂v /∂y)2 ∼ y 2 , mentre nel modello originale comportamento migliore: lm ∼ y 2 e τxy ∼ 2lm x era lm ∼ y e τxy ∼ y, e nella realtà τxy ∼ vx vy ∼ y 3 . Corrsin e Kistler (1954) e Klebanoff (1956), al fine di tenere conto dell’elevata intermittenza riscontrata nella regione più esterna dello strato limite, hanno proposto di moltiplicare la viscosità turbolenta per la funzione y 6 −1 Fk (y, δ) = 1 + 5.5 δ (10.14) dove δ è lo spessore di strato limite. Infine Cebeci e Smith (1967) e Baldwin e Lomax (1978) combinano le due varianti viste sopra; secondo Cebeci e Smith: ( νT i se y ≤ ym νT = νT o se y ≥ ym dove ym è il valore più piccolo di y per il quale νT i = νT o . νT i e νT o sono definite da: νT i = lm " ∂ hui ∂y 2 + ∂ hvi ∂x 2 #1/2 CAPITOLO 10. MODELLI DI TURBOLENZA PER EQUAZIONI MEDIATE 116 ALLA REYNOLDS con lm dp/ dx A = 26 1 + y ρu2τ i h + + = κy 1 − e−y /A , + −1/2 questa definizione di A+ permette di trattare situazioni con gradienti di pressione longitudinali. νT o = αUe δV∗ Fk (y, δ) con α = 0.0168, Ue la velocità longitudinale esterna, Fk la funzione di Klebanoff (10.14) e Z δ u ∗ δV = 1− dy Ue 0 Altri esempi In correnti autosimilari come getti, scie o strati di mescolamento, si pone frequentemente νT = U0 (x)δ(x) ReT dove U0 (x) e δ(x) sono velocità e lunghezze di riferimento. ReT è una costante, il numero di Reynolds turbolento, determinata empiricamente per il tipo di flusso considerato: • getto cilindrico ReT ∼ 35 • getto piano ReT ∼ 31 • strato di mescolamento ReT ∼ 60 % 110 • scia piana ReT ∼ 13 % 19 • scia assialsimmetrica ReT ∼ 2 % 22 I difetti più evidenti di questo tipo di modelli sono che la costante del modello dipende fortemente dal tipo di flusso considerato, e che il modello non è invariante per trasformazioni di Galileo. 10.4 Modelli a 1 equazione La viscosità turbolenta viene ipotizzata come funzione di una quantità turbolenta da determinare risolvendo un’equazione differenziale. La prima proposta fu di Prandtl e Kolmogorov, durante la seconda guerra mondiale, √ che hanno assunto come velocità caratteristica la radice dell’energia cinetica turbolenta k, quindi √ νT = l k (10.15) L’equazione di evoluzione per l’energia cinetica turbolenta può essere ottenuta ponendo i = j nell’eq.(10.8) e sommando per i = 1, .., 3: ∂Ui ∂k ∂Ti ∂k + Uk = τijR −ǫ− ∂t ∂xk ∂xj ∂xi (10.16) 10.5. MODELLI A 2 EQUAZIONI dove 117 ∂ ′ ′ 1 ′ ′2 ′ ′ ν ∂ v ′2 −ν vv vi v + p vi − Ti = 2 2 ∂xi ∂xj i j i τijR ∂U ∂xj è il termine di produzione associato al lavoro degli sforzi turbolenti, responsabile del trasferimento di energia turbolenta dal moto D medioE a quello turbolento. Il tasso di dissipa′ S′ zione di energia cinetica turbolenta ǫ = 2ν Sij ij deriva dal lavoro degli sforzi turbolenti che dissipano energia cinetica turbolenta con un corrispondente aumento di energia interna. i Il termine ∂T ∂xi rappresenta il trasporto turbolento e viene modellato come un processo di diffusione. Si pone: νT ∂k 1 ′ ′ ′ 1 ′ ′ vvv + p vj = − 2 i i j ρ σk ∂xj dove σD k è un E coefficiente da determinare su base empirica; inoltre si trascura il termine ∂ ′ ′ −ν ∂xj vi vj . L’equazione (10.16) risulta quindi cosı̀ semplificata τijR ∂Ui ∂k ∂k ∂ ∂k + Uj = −ǫ+ (ν + νT /σk ) ∂t ∂xj ρ ∂xj ∂xj ∂xj (10.17) In base a considerazioni dimensionali (Prandtl 1945) dovrà essere ǫ = CD k3/2 /l dove la lunghezza caratteristica l è da Ddeterminare. Si sottolinea che in generale νT è il E ′ ′ rapporto tra una quantità turbolenta (− vi vj ) e una del moto medio (hSij i), quindi soltanto per turbolenza in equilibrio le scale turbolente e del moto medio sono proporzionali e possono essere usate indifferentemente per νT , altrimenti è necessaria una combinazione di scale, ed il predetto modello non è valido. Successivamente Nee e Kovasnay (1968), Baldwin e Bar (1990), Spalart e Allmaras (1992) hanno risolto un’equazione di trasporto-diffusione per νT . Questi modelli hanno il difetto di richiedere un’espressione algebrica per la lunghezza caratteristica (fortemente dipendente dal tipo di corrente) da determinare su base empirica. 10.5 Modelli a 2 equazioni Kolmogorov per primo (1942) suggerı̀ il fatto che la turbolenza può essere descritta adeguatamente da due grandezze indipendenti, determinabili quindi mediante equazioni differenziali. Il modello k − ǫ Il modello k − ǫ è il modello a due equazioni più diffuso. La prima proposta risale a Chou (1945), ma il primo contributo significativo è dovuto a Jones e Launder (1972) con il modello denominato k − ǫ standard. Alcune variazioni nella determinazione dei coefficienti di chiusura sono dovuti a Launder e Sharma (1974). Si parte dalle equazioni (10.10) e (10.16) con ν T = Cµ k2 ǫ CAPITOLO 10. MODELLI DI TURBOLENZA PER EQUAZIONI MEDIATE 118 ALLA REYNOLDS a cui bisogna aggiungere un’equazione per il tasso di dissipazione di energia turbolenta ǫ. Per ricavarla si può partire dall’equazione (10.6); scriviamola nella forma ∂vi′ = Ni (U , v ′ , p′ ) ∂t ′ : Formiamo l’equazione di evoluzione per Sij ′ ∂Sij ∂ ∂vj′ ∂ ∂vi′ ∂Nj ∂Ni = + = + 2 ∂t ∂xi ∂t ∂xj ∂t ∂xi ∂xj Ricordando la definizione del tasso di dissipazione ′ ′ ǫ = 2ν Sij Sij arriviamo a scrivere l’equazione di evoluzione per ǫ: ∂Nj ∂ǫ ∂Ni ′ = 4ν Sij + ∂t ∂xi ∂xj con numerosi calcoli si può dunque scrivere un’equazione differenziale per ∂ǫ ∂t . Questa equazione però è poco utile perché i termini da modellare sono troppo complessi. Si preferisce dunque scrivere un’equazione ottenuta sulla base di considerazioni fisiche e di un’analisi dimensionale. Come per l’equazione di evoluzione per k, essa conterrà un termine di produzione, uno di dissipazione e uno di trasporto turbolento. i Il termine di produzione analogamente a quello per k deve contenere la quantità τijR ∂U ∂xj i e per ragioni dimensionali esso dovrà essere Cǫ1 kǫ τijR ∂U ∂xj con Cǫ1 costante adimensionale. Il termine di dissipazione si ottiene con considerazioni dimensionali: esso dipenderà da ǫ e da 2 k e l’unica possiblità è −Cǫ2 ǫk ; infine il termine di trasporto turbolento viene scritto nella forma di un termine diffusivo con contributi di viscosità molecolare e viscosità turbolenta. Il modello k − ǫ standard risulta cosı̀ formulato: ∂k ∂k ∂ νT ∂k R ∂Ui + Uj ) = τij −ǫ+ (ν + ∂t ∂xj ∂xj ∂xj σk ∂xj ∂ǫ ǫ R ∂Ui ǫ2 ∂ νT ∂ǫ ∂ǫ + Uj = Cǫ1 τij − Cǫ2 + ) (ν + (10.18) ∂t ∂xj k ∂xj k ∂xj σǫ ∂xj νT = Cµ k2 /ǫ, Cǫ1 = 1.44, Cǫ2 = 1.92, Cµ = 0.09, σK = 1., σǫ = 1.3 Modello k − ω Questo modello è stato ricavato in origine da Kolmogorov (1942) ed è stato molto sviluppato in seguito. Il significato di ω cambia a seconda degli autori: Saffman (1970) utilizzò un’equazione per ω 2 , inteso come il quadrato medio della vorticità associata agli eddies energetici. Launder e Spalding (1972) attribuirono ad ω il significato di RMS delle fluttuazioni di vorticità, e quindi a ω 2 il doppio dell’enstrofia. Wilcox (1988) e Speziale (1990) definirono ω = ǫ/k cioè il rapporto fra il tasso di dissipazione e l’energia cinetica turbolenta. 10.5. MODELLI A 2 EQUAZIONI 119 Nella formulazione originale Kolmogorov scrisse l’equazione di evoluzione (10.17) per l’energia cinetica turbolenta, e la seguente equazione di evoluzione per ω: ∂ω ∂ ∂ω ∂ω 2 + Uj = −βω + σνT ∂t ∂xj ∂xj ∂xj Questa equazione è stata ricavata da Kolmogorov con considerazioni simili a quelle che hanno condotto all’equazione (10.18) per ǫ. In essa si riconosce un terminedi dissipazione −βω 2 e ∂ω un termine di trasporto turbolento scritto in forma di diffusione ∂x∂ j σνT ∂x . Sono assenti j invece il termine di produzione e quello di diffusione molecolare. Negli sviluppi successivi del modello entrambi questi termini sono stati inclusi. Il termine di produzione deve contenere ω R ∂Ui i il prodotto τijR ∂U ∂xj e per ragioni dimensionali l’unica possibilità è: α k τij ∂xj . Si ottiene cosı̀ il modello k − ω standard: ∂k ∂k ∂ ∂k R ∂Ui ∗ ∗ + Uj = τij − β kω + (ν + σ νT ) ∂t ∂xj ∂xj ∂xj ∂xj ∂ω ∂ω ω ∂Ui ∂ ∂ω + Uj = α τijR − βω 2 + (ν + σνT ) ∂t ∂xj k ∂xj ∂xj ∂xj νT = k , ω α= 13 , 25 β= 9 , 125 β∗ = 9 , 100 σ = σ∗ = 1 2 La variante più interessante del modello è stata proposta da Wilcox (1993) per correnti con shear ed estendibile anche a correnti comprimibili: ∂k ∂ ∂k ∂k R ∂Ui ∗ ∗ + Uj = τij − β kω + (ν + σ νT ) (10.19) ∂t ∂xj ∂xj ∂xj ∂xj ∂ω ω R ∂Ui ∂ ∂ω ∂ω 2 + Uj = α τij − βω + (ν + σνT ) (10.20) ∂t ∂xj k ∂xj ∂xj ∂xj νT = α= β0 = 13 , 25 9 , 125 k ω (10.21) 1 1 β ∗ = β0∗ fβ ∗ , σ = , σ ∗ = , 2 2 Ωij Ωij Ski 1 + 70χω , β∗ = 9 fβ = , χω = 0 1 + 80χω (β0∗ ω)3 100 ( 1, χK ≤ 0 fβ ∗ = 1+680χ2K , χK ≥ 0 1+400χ2 β = β0 fβ , K χk = 1 ∂k ∂ω , ω 3 ∂xj ∂xj ǫ = β ∗ ωk, l = k1/2 /ω Si sottolinea che in correnti 2D χω = 0 ed il termine sorgente si annulla. Inoltre questo modello k − ω può essere utilizzato anche nello strato limite, dove ω è alta, i fattori χK e χω diventano piccoli riducendo i termini sorgenti. I coefficienti e le funzioni fβ ∗ , fβ , α, β0 sono stati calibrati per scie, mixing layers e getti, ma il modello ha dato buoni risultati anche in correnti di parete. CAPITOLO 10. MODELLI DI TURBOLENZA PER EQUAZIONI MEDIATE 120 ALLA REYNOLDS Valutazioni L’applicazione di modelli a 1 equazione è poco conveniente perchè non dà risultati sensibilmente migliori rispetto ai modelli algebrici. Inoltre i modelli a 1 equazione richiedono che sia definita per ogni tipo di corrente la lunghezza caratteristica l. L’unico loro vantaggio è la semplice implementabilità rispetto ai modelli a 2 equazionni. Il modello k − ǫ, nonostante sia il più diffuso tra quelli a due equazioni, • è difficile da implementare, • deve essere accuratamente calibrato per ogni nuova applicazione, • richiede correzioni nello strato limite, • è assolutamente inadeguato a simulare correnti con gradiente di pressione contrario. Il modello k − ω presenta invece alcuni vantaggi: • non richiede correzioni nello strato limite, • riproduce adeguatamente anche correnti con gradiente di pressione contrario ed in transizione. • è stato applicato anche a correnti separate e con ricircolo senza modifiche rispetto al modello base. Questo modello non è però adatto a simulare l’interazione onda d’urto-strato limite e correnti su pareti curve. 10.6. CALIBRAZIONE DEI COEFFICIENTI DEL MODELLO K − ω 10.6 121 Calibrazione dei coefficienti del modello k − ω a cura di Andrea Mola Allo scopo di meglio comprendere come vengano sviluppati i modelli di turbolenza, studieremo come siano scelti i coefficienti di chiusura per un modello a due equazioni. In particolare, si farà riferimento al modello k − ω di Wilcox (1988), e sarà dunque seguita la procedura esposta da Wilcox stesso nel capitolo 4 del libro T̈urbulence modeling for CFD.̈ Partiamo dunque dalla scrittura delle equazioni per k ed ω relative a tale modello, che sono rispettivamente ∂k ∂k + vj ∂t ∂xj ∂ω ∂ω + vj ∂t ∂xj con ∂v i ∂ ∂k − β ∗ kω + (ν + σ ∗ νT ) ∂xj ∂xj ∂xj ∂v i ω ∂ ∂ω = α τijR − βω 2 + (ν + σνT ) k ∂xj ∂xj ∂xj = τijR (10.22) (10.23) 1 τijR = νT Sij − kδij , νT = k/ω, 3 e dove Sij è la parte simmetrica dello jacobiano del campo di velocità media. Ci sono dunque cinque coefficienti di chiusura da determinare: α, β ∗ , β, σ ∗ , σ. Solitamente, i valori di questi coefficienti di chiusura vengono scelti in modo da riprodurre risultati sperimentali per la legge di decadimento temporale della turbolenza omogenea isotropa, e per alcuni flussi di riferimento, che nella maggior parte dei casi sono gli strati limite all’equilibrio, ed i cosı̀ detti free shear flows (getti, strati di mescolamento, scie). Il motivo per cui gli strati limite ed i free shear flows sono i principali benchmark per i modelli di turbolenza, si basa principalmente sulla supposizione che se un modello mostra di dare buoni risultati nella simulazione di fenomeni viscosi sia in prossimità di pareti solide, sia lontano da queste, allora si può ritenere in grado, almeno dal punto di vista qualitativo, di simulare praticamente tutti i flussi di interesse industriale. Nel nostro caso dunque, i valori dei coefficienti di chiusura saranno ottenuti studiando il comportamento del modello di turbolenza in esame nel caso di turbolenza omogenea isotropa, e di varie regioni dello strato limite turbolento. In seguito, analizzeremo le soluzioni numeriche ottenute impiegando questi coefficienti, per i freee shear flows. Osserveremo infine come i risultati per quest’ultima categoria di flussi possano essere migliorati tramite delle correzioni ai coefficienti di chiusura. 10.6.1 Turbolenza omogenea isotropa Nel caso di turbolenza omogenea isotropa, nelle equazioni per k ed ω si annullano tutti i gradienti spaziali. Si ottiene dunque il sistema di equazioni differenziali ordinarie nella sola variabile temporale t ∂k = −β ∗ kω ∂t ∂ω = −βω 2 ∂t CAPITOLO 10. MODELLI DI TURBOLENZA PER EQUAZIONI MEDIATE 122 ALLA REYNOLDS dal quale si può ricavare per via analitica l’andamento temporale dell’energia turbolenta ∗ k∼t − ββ . Poiché la legge di decadimento osservata sperimentalmente è k ∼ t−n , in cui l’esponente n = 1.25 ± 0.06, possiamo scrivere β∗ 5 = , β 4 (10.24) prima relazione derivata da osservazioni sperimentali, tra i nostri coefficienti di chiusura. 10.6.2 Strato limite Strato logaritmico Come ben noto, in questa regione dello strato limite, la componente orizzontale della velocità media v x segue la legge vx 1 (10.25) = log y + + A V∗ κ dove κ = 0.41 è la costante di Von Kàrmàn, ed A = 5.2. Scriviamo ora le equazioni di moto e quelle delle variabili turbolente per lo strato logaritmico: in tale regione si suppone che la velocità sia abbastanza piccola da consentire di trascurare i termini non lineari nelle equazioni RANS, e che la viscosità molecolare sia trascurabile rispetto a quella turbolenta. Inoltre sfruttiamo le ipotesi relative allo strato limite stazionario su lastra piana orizzontale: la componente verticale della velocità è trascurabile rispetto a quella longitudinale alla lastra, mentre le derivate in direzione orizzontale sono trascurabili rispetto a quelle in direzione verticale. Otteniamo dunque ∂v x ∂ νT 0 = ∂y ∂y ∂v x 2 ∂ ∂k ∗ ∗ 0 = νT − β ωk + σ νT (10.26) ∂y ∂y ∂y 2 ∂v x ∂ ∂ω 0 = α − βω 2 + σ νT ∂y ∂y ∂y Poiché si desidera che la velocità segua la legge logaritmica (10.25) analizzando la prima equazione di questo sistema si ricava che νT = k/ω ∼ y, mentre dalla seconda avremo kω ∼ 1 , y se supponiamo che k ∼ y n , ω ∼ y m , otteniamo da queste due equazioni il semplice sistema lineare n + m = −1 n−m = 1 dal quale concludiamo che 10.6. CALIBRAZIONE DEI COEFFICIENTI DEL MODELLO K − ω k ∼ cost., ω∼ 1 . y 123 (10.27) Con analoghe considerazioni sulle costanti, possiamo dedurre dunque che ∂v x V∗ = , ∂y κy V∗ 2 k = √ ∗, β ω=√ V∗ . β ∗ κy (10.28) Inserendo queste nel sistema (10.26), si osserva come la prima e la seconda equazione siano identicamente soddisfatte, mentre dalla terza equazione si ricava una seconda relazione tra i coefficienti di chiusura, derivata dalla richiesta di un andamento di velocità logaritmico con coefficiente κ nel log layer, α= β∗ σκ2 − √ ∗. β β (10.29) Dall’osservazione invece della prima equazione del sistema (10.26), si deduce che lo sforzo tangenziale τxy è costante nel log layer, e pari a τxy = V∗2 = k p β∗. Risultati sperimentali confermano questa circostanza, e fissano ad un valore di circa 3/10 il coefficiente di proporzionalità tra k e τxy . Quindi scegliamo β∗ = 9 . 100 Fissato il valore di β ∗ , dalla (10.24), risulta automaticamente che 9 . 125 Analizzando dunque il comportamento del nostro modello di turbolenza nel log layer, siamo riusciti ad ottenere due ulteriori relazioni tra i coefficienti di chiusura. Come vedremo, la richiesta che il modello di turbolenza sia in grado di riprodurre i risultati sperimentali anche nel substrato viscoso e nel defect layer, ci consentirà di fissare i coefficienti residui. Tuttavia, non sarà possibile giungere in questi casi a delle relazioni esplicite simili a quelle trovate in precedenza, ma ci si dovrà accontentare di valori che di volta in volta consentono alle soluzioni numeriche di avvicinarsi il più possibile alle osservazioni sperimentali. β= Substrato viscoso A causa della presenza delle due equazioni relative al modello di turbolenza, non è possibile ottenere per questa regione dello strato limite, una soluzione autosimilare, simile a quella ottenuta nel caso laminare, o nel caso di modelli di turbolenza algebrici (mixing length). Tuttavia, per potersi ancora ricondurre ad un sistema di equazioni differenziali ordinarie nella variabile y + = V∗ y/ν, si lavora su una perturbazione del problema originale, asintotica a quest’ultimo per Reδ∗ → ∞, dove Reδ∗ è il numero di Reynolds basato sullo spessore di spostamento Z ∞ vx ∗ dy. δ = 1− U 0 CAPITOLO 10. MODELLI DI TURBOLENZA PER EQUAZIONI MEDIATE 124 ALLA REYNOLDS Figura 10.4: Substrato viscoso: confronto tra risultati sperimentali e la soluzione numerica per substrato ottenuta usando il modello k − ω Le soluzioni numeriche cosı̀ ottenute possono essere confrontate con i risultati sperimentali. È importante che la soluzione numerica tenda alla soluzione analitica dello strato logaritmico per y + → ∞. In particolare, dalla soluzione nel substrato viscoso dipenderà il valore calcolato numericamente An della costante A nell’equazione 10.25. Data la soluzione numerica vn (y + ), avremo 1 + + An = lim vn (y ) − log y . κ y + →∞ Wilcox osserva come il comportamento della soluzione numerica nel substrato viscoso dipenda quasi esclusivamente dal coefficiente σ, e che un valore pari ad 1/2 consente un’ottima riproduzione dei dati sperimentali. A tal proposito, riportiamo nella Tabella 10.1 il valore numerico di A per vari modelli di turbolenza, e mostriamo in Figura 10.4 il grafico di un confronto tra risultati sperimentali e numerici. Modello di turbolenza k−ω Wilcox (1988) k−ǫ Launder–Sharma (1974) k − ω 2 Launder–Spalding (1972) k−ω Kolmogorov (1942) An 5.1 -2.2 5.7 3.1 Tabella 10.1: Costanti An calcolate con diversi modelli a due equazioni Si nota immediatamente come ci siano alcuni modelli (soprattutto il k − ǫ) che non si comportano adeguatamente nel substrato viscoso, cosa che ne rende l’utilizzo sconveniente in 10.6. CALIBRAZIONE DEI COEFFICIENTI DEL MODELLO K − ω 125 tale regione. Questo è uno dei motivi per cui vengono utilizzate le cosı̀ dette wall functions, che consentono di non simulare il flusso all’interno di questa regione, imponendo delle condizioni derivate dalle (10.28) su una parete fittizia che si trovi ad un’opportuna distanza da quella reale. Tornando dunque ai nostri coefficienti di chiusura, avendo scelto σ = 1/2 e sfruttando la (10.29) otteniamo 13 . 25 Abbiamo quindi determinato il valore di quattro dei cinque coefficienti di chiusura. Per determinare il valore di σ ∗ analizzeremo il comportamento del modello di turbolenza nel defect layer. α= Defect Layer Questa zona dello strato limite riveste particolare importanza, essendo la regione in cui si fa maggiormente sentire la presenza di un eventuale gradiente di pressione esterno dp/dx. Anche in questo caso non è possibile giungere ad una soluzione autosimilare, se non considerando una perturbazione asintotica del problema, questa volta ottenuta per V∗ /U → 0 (condizione verificata ad alti Rex ). Per il problema perturbato si è addirittura in grado di trovare una soluzione in forma chiusa, nell’ipotesi che sia piccola la variabile di similarità η = y/∆(x), con ∆(x) = U δ∗ /V∗ . In questo caso, la legge di difetto trovata è 1 U − vn = − log η + v0 − v1 η log η V∗ κ (10.30) dove e √ (β/αβ ∗ )(σ ∗ κ/2 β ∗ ) √ v1 = βT (1 − β/αβ ∗ )(σ ∗ κ2 /2 β ∗ ) δ∗ dp . τw dx Uno sguardo all’espressione della costante v1 , mostra come quest’ultima dipenda linearmente da βT : l’ultimo termine nella (10.30) dunque, provocherà un’inflessione nelle curve del difetto di velocità, quando βT 6= 0. Poiché da un’ulteriore osservazione dell’espressione di v1 , ci si rende conto che l’entità di tale inflessione è pressochè poroporzionale a σ ∗ , ultimo coefficiente di chiusura incognito, fisseremo quest’ultimo in modo che i risultati numerici ottenuti applicando il modello di turbolenza in esame, riproducano il più fedelmente possibile quelli sperimentali. Con questa motivazione, l’ultimo coefficiente di chiusura del modello k − ω di Wilcox viene fissato ad un valore σ ∗ = 1/2. In Figura 10.5 viene mostrato il confronto tra risultati sperimentali e numerici (ottenuti tramite diversi modelli di turbolenza) nel defect layer, per βT = 0 (sinistra) e βT = 8.7 (destra). Si nota come il modello k − ω ottenga risultati soddisfacenti sia nel caso di gradiente di pressione esterna nullo, che nel caso di gradiente di pressione esterna significativamente βT = CAPITOLO 10. MODELLI DI TURBOLENZA PER EQUAZIONI MEDIATE 126 ALLA REYNOLDS Figura 10.5: Defect layer: confronto tra risultati sperimentali e numerici per il difetto di velocità (U − v x )/V∗ ; — k − ω; - - k − ω 2 ; –·– k − ǫ diverso da zero. Al contrario gli altri modelli di turbolenza citati non sono in grado di riprodurre correttamente i risultati degli esperimenti per βT 6= 0. Riassumendo, i valori dei coefficienti di chiusura per il modello k − ω di Wilcox saranno α= 10.6.3 13 , 25 β∗ = 9 , 100 β= 9 , 125 σ∗ = 1 , 2 σ= 1 . 2 Free shear flows Vedremo ora come il modello analizzato si comporta nella simulazione di getti (Figura 10.6), scie, e strati di mescolamento. Per farlo, presenteremo risultati numerici ottenuti usando il modello k − ω, con quelli relativi ad altri modelli di turbolenza, e con risultati sperimentali. Per la simulazione numerica di questo tipo di flussi, si ricorre ancora una volta all’autosimilarità, in modo da ricondursi ad un sistema di equazioni differenziali ordinarie nella variabile autosimilare η = y/x, le cui incognite sono la velocità del fluido e le grandezze turbolente caratteristiche di ciascun modello di turbolenza. In Tabella 10.2 sono riportate le tangenti degli angoli di semiapertura per getti (in Figura 10.6 è mostrato l’angolo α di semiapertura del getto), scie, e strati di mescolamento. È possibile notare immediatamente come il modello k − ω fornisca in questo caso risultati sensibilmente meno accurati del k − ǫ. Più precisamente, il modello k − ω predice valori di viscosità turbolenta troppo elevati, che portano ad un’eccessiva diffusione. Questo è il motivo per cui getti, scie e strati di mescolamento calcolati con il modello k − ω sono più larghi di quelli calcolati con il modello k − ǫ, e rilevati sperimentalmente. Per spiegare questa differenza, si può√considerare il modello k−ǫ ed applicare all’equazione per ǫ un cambiamento di variabili ǫ = β ∗ kω, in modo da trovare un’equivalente equazione 10.6. CALIBRAZIONE DEI COEFFICIENTI DEL MODELLO K − ω 127 y α x, v x Figura 10.6: Free shear flows: l’esempio del getto Modello di turbolenza k−ω Wilcox (1988) k − ǫ Launder–Sharma (1974) Risultati sperimentali Getto Piano 0.135 0.108 0.100-0.110 Getto Ax. 0.369 0.120 0.086-0.096 Scia 0.496 0.256 0.365 Strato Mesc. 0.141 0.098 0.115 Tabella 10.2: tan α calcolata e misurata, per i free shear flows per ω. Dal confronto con l’equazione per ω del modello k − ω, si osserva la comparsa del termine aggiuntivo (detto termine di crossflow) 2σ 1 ∂k ∂ω ω ∂xj ∂xj nel secondo membro dell’equazione (10.23). L’interessante caratteristica di questo termine, è che aumenta la dissipazione, diminuendo cosı̀ k, e di conseguenza νT . Molti ricercatori (Wilcox, Menter, Speziale e Peng tra gli altri) hanno dunque tentato di includere questo termine nei propri modelli, con risultati tuttavia spesso insoddisfacenti, sia perchè nella maggior parte dei casi il modello di turbolenza ottenuto produce troppa dissipazione, sia perchè il termine di crossflow può dare qualche problema all’interno dello strato limite, dove il modello k−ω originale invece dava risultati soddisfacenti. Inoltre questo termine è di difficile trattamento numerico poiché modifica il termine convettivo nell’equazione per ω ∂ω ∂ω 2σ ∂k −→ . vj vj − ∂xj ω ∂xj ∂xj CAPITOLO 10. MODELLI DI TURBOLENZA PER EQUAZIONI MEDIATE 128 ALLA REYNOLDS Essendo infatti ω molto piccolo lontano da parete, la nuova velocità può risultare molto elevata costringendo ad usare time step molto ridotti per l’integrazone temporale. Per aggirare questi problemi, e basandosi proprio sull’osservazione che il termine di crossflow risulta particolarmente elevato lontano da parete, ma praticamete nullo nello strato limite, Wilcox (1993) ha proposto la seguente correzione per il coefficiente di chiusura β ∗ β ∗ = β0∗ fβ ∗ (χk ) con β0∗ = 9 , 100 680 + χk 2 fβ ∗ (χk ) = 400 + χk 2 1 χk > 0 , χk < 0 χk = 1 ∂k ∂ω . ω 3 ∂xj ∂xj Come possiamo osservare, il coefficiente di chiusura β ∗ , non è più una costante, ma dipende ora dalla funzione delle varabili turbolente χk , adimensionalizzazione del termine di crossflow analizzato precedentemente. Va dunque notato come questa correzione, cosı̀ come il termine di crossflow, ha effetti apprezzabili esclusivamente lontano da pareti solide, cioè nelle sole regioni in cui il comportamento del modello andava corretto. Infine si fa notare come questa correzione compaia nell’equazione per k anziché in quella per ω, per garantire un controllo più diretto sull’effetto desiderato (νT è proporzionale a k). In Tab. 10.3 riportiamo i risultati ottenuti applicando al modello k − ω la correzione appena presentata. Getto Piano 0.101 Getto Ax. 0.139 Tabella 10.3: tan α per il modello k − ω con il termine β corretto Come auspicato, i valori calcolati per le ampiezze dei getti, si riducono e divengono compatibili con quelli misurati. Tuttavia (e questo è un difetto comune a molti modelli a due equazioni, quali il k − ǫ) sebbene negli esperimenti si osservi come l’ampiezza del getto piano sia maggiore di quella del getto assialsimmetrico, le soluzioni numeriche mostrano un comportamento opposto. Per ovviare a questo problema, conosciuto come anomalia del getto piano/getto assialsimmetrico, Pope (1970) ha proposto di introdurre modifiche che aumentino il termine di produzione di dissipazione nel caso di flusso sia tridimensionale. In questo modo, l’aumento della dissipazione, contribuirebbe a diminuire k, e qundi νT per tali flussi. Infatti, un meccanismo con cui l’energia viene passata dalle grandi alle piccole scale, dove agisce la dissipazione, è lo stiramento dei vortici sulla scala del moto medio. Questo meccanismo è tipicamente tridimensionale, e dunque non entra in azione nel caso di moto piano. Viene cosı̀ introdotta una correzione al coefficiente β di produzione per la dissipazione ω della forma β = β0 fβ (χω ) dove 9 , β0 = 125 1 + 70χω fβ (χω ) = , 1 + 80χω Ωij Ωjk Ski . χω = (β ∗ ω)3 10.6. CALIBRAZIONE DEI COEFFICIENTI DEL MODELLO K − ω 129 Essendo Ωij e Sij la parte antisimmetrica e simmetrica rispettivamente della matrice jacobiana del campo di velocità medio, la funzione χω è nulla nel caso di moto medio piano. Infatti in questo caso avremo un tensore antisimmetrico del tipo 0 −a 0 Ωij = a 0 0 0 0 0 e quindi Ωij Ωjk −a2 0 0 = 0 −a2 0 0 0 0 −→ Ωij Ωjk Ski = −2a 2 ∂v x ∂v y + ∂x ∂y =0 nel caso incomprimibile. Si noti infine come il termine χω , proporzionale ad ω −3 , risulti pressoché nullo nello strato limite, dove la correzione non deve agire. Riassumiamo dunque il modello k − ω cosı̀ ottenuto: ∂k ∂k + vj ∂t ∂xj ∂ω ∂ω + vj ∂t ∂xj ∂v i ∂ ∂k − β ∗ kω + (ν + νT ) ∂xj ∂xj ∂xj ∂v i ω ∂ ∂ω = α τijR − βω 2 + (ν + νT ) , k ∂xj ∂xj ∂xj = τijR 1 τijR = νT Sij − kδij , 3 α= β0∗ = 13 , 25 9 , 100 β ∗ = β0∗ fβ ∗ (χk ), 1 σ∗ = , 2 β = β0 fβ (χω ), 680 + χk 2 fβ ∗ (χk ) = 400 + χk 2 1 9 , β0 = 125 νT = k/ω, χk > 0 , χk = χk < 0 1 + 70χω fβ (χω ) = , 1 + 80χω 1 σ= , 2 1 ∂k ∂ω , ω 3 ∂xj ∂xj Ωij Ωjk Ski . χω = (β ∗ ω)3 In Tabella 10.4 sono infine riportati i risultati numerici ottenuti impiegando questo modello k − ω di Wilcox (1993) nella simulazione dei free shear flows. Modello di turbolenza k−ω Wilcox (1993) Risultati sperimentali Getto Piano 0.101 0.100-0.110 Getto Ax. 0.088 0.086-0.096 Scia 0.339 0.365 Strato Mesc. 0.105 0.115 Tabella 10.4: tan α calcolata e misurata, per i free shear flows È possibile quindi concludere che modello k − ω in cui sono presenti le correzioni analizzate, si comporta molto bene nella simulazione dei free shear flows. Poiché poi i termini correttivi sono stati costruiti in modo da non avere significativi effetti nella vicinanza di pareti CAPITOLO 10. MODELLI DI TURBOLENZA PER EQUAZIONI MEDIATE 130 ALLA REYNOLDS solide, questo modello conserverà le buone prestazioni per la simulazione degli strati limite all’equilibrio, analizzata nella sezione 10.6.2. Si è dunque ottenuto un modello di turbolenza che presenta un buon comportamento per tutti i flussi di benchmark considerati.