考慮轉(zhuǎn)捩的SST模型在沖擊射流傳熱中的應用

黃華坤,孫鐵志,尤天慶,張桂勇,3,4,回達,宗智,3,4

(1.大連理工大學船舶工程學院,116024,遼寧大連;2.北京宇航系統(tǒng)工程研究所,100076,北京; 3.大連理工大學工業(yè)裝備與結構分析國家重點實驗室,116024,遼寧大連; 4.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心,200240,上海)

沖擊射流由于具有強烈的換熱效果,受到了各行業(yè)的廣泛關注[1]。對于垂直或短距起降(V/STOL)的軍用飛機,射流與地面之間的相互作用對飛機的設計具有重要的影響[2]。在航母上起降的V/STOL飛機不僅牽涉到飛機本身,更有可能影響到甲板結構的安全。除此之外,沖擊射流原理廣泛應用于造紙、紡織等涉及到的干燥、冷卻和加熱等工業(yè)過程。

早期的研究表明,使用k-ε模型在傳熱問題的預測上與實驗值偏差很大,而且在沖擊距離較短時無法準確地預測出努塞爾數(shù)Nu的第二峰值[2-6]。此外,前人的研究表明標準兩方程模型容易過高預測滯止點附近的湍動能[7]。Kato和Launder從滯止點附近的流體趨于無旋、渦量趨于零的角度出發(fā),提出了基于渦量和應變的湍動能產(chǎn)生項,以降低滯止點附近湍動能的生成[8]。因而,Kato-Launder模型被推薦應用于沖擊射流問題[9]。

然而,努塞爾數(shù)第二峰值依然很難被湍流模型準確地預測到。研究表明,努塞爾數(shù)第二峰值受層流到湍流轉(zhuǎn)捩的影響[10]。Uddin等通過LES模型進一步揭示了渦的形成和發(fā)展對努塞爾數(shù)第二峰值的影響[11]。Sharif等在研究中指出具備計算層流到湍流轉(zhuǎn)捩能力的模型能更好地預測努塞爾數(shù)第二峰值[12],當前廣泛使用的轉(zhuǎn)捩模型為低雷諾數(shù)湍流模型。低雷諾數(shù)轉(zhuǎn)捩模型是一種基于黏性底層特性、通過阻尼因子修改底層湍流黏度的方法來模擬轉(zhuǎn)捩的模型,因此模擬轉(zhuǎn)捩的能力有限[13]。Dutta等通過數(shù)值計算表明,使用低雷諾數(shù)SSTk-ω模型能獲得和非轉(zhuǎn)捩模型相比更準確的努塞爾數(shù)分布[14]。然而,在沖擊距離增大時,低雷諾數(shù)SSTk-ω模型的預測能力差于標準k-ε和k-ω模型,因為低雷諾數(shù)SSTk-ω模型預測了偽二次峰值。

因此,準確地模擬沖擊射流的傳熱現(xiàn)象仍然是國內(nèi)外學者重點關注的問題。本文在文獻[15]推薦的SSTk-ω模型的基礎上,加入了Kato-Launder模型以避免湍流模型過高地預測滯止點附近的湍動能。同時,為了更好地預測層流到湍流的轉(zhuǎn)捩過程,將一方程湍流間歇轉(zhuǎn)捩模型加入到求解方程中,并通過與實驗數(shù)據(jù)對比來評價改進的SSTk-ω湍流模型在射流沖擊問題中的適用性。一方程間歇轉(zhuǎn)捩模型是基于層流到湍流的轉(zhuǎn)捩機理,利用間歇因子控制湍流的生成,同時基于經(jīng)驗公式對轉(zhuǎn)捩過程進行求解而建立的模型[16]。基于建立的數(shù)值方法探究了不同沖擊高度下壓力梯度對努塞爾數(shù)第二峰值的影響,獲得的結果可豐富射流換熱問題的機理內(nèi)涵。

1 數(shù)值計算方法

1.1 基本控制方程

控制方程包括連續(xù)性方程、動量方程和能量方程,表達式如下

(1)

(3)

1.2 SST k-ω湍流模型

SSTk-ω模型是由Menter提出的一種基于標準k-ε模型和k-ω模型的混合模型[17],因此SSTk-ω模型在壁面射流區(qū)能準確地預測流動特征,同時可避免對遠場邊界條件的敏感性。該模型還添加了逆壓梯度修正項,這使得SSTk-ω模型廣泛應用于如沖擊射流等逆壓梯度和分離流問題中。除了前面提到的控制方程外,為封閉模型,在模型中增加了湍動能k和比耗散率ω,表達式如下

(4)

(5)

Gω-Yω+Dω+Sω

(6)

1.2.1Kato-Launder模型SSTk-ω等兩方程普遍存在的一個問題是產(chǎn)生了過高的湍動能[7],導致SSTk-ω模型在滯止點處預測的傳熱率過高。Kato和Launder對產(chǎn)生的湍動能項進行了如下修正[8]

Gk=μtS2

(7)

Gk=μtSΩ

(8)

1.2.2 湍流間歇轉(zhuǎn)捩模型 間歇轉(zhuǎn)捩模型從TransitionSST模型中發(fā)展而來,該模型考慮了自然轉(zhuǎn)捩、分離流轉(zhuǎn)捩[18]和橫流不穩(wěn)定[7]等多種轉(zhuǎn)捩機制,被廣泛用于三維機翼轉(zhuǎn)捩預測問題。間歇轉(zhuǎn)捩模型的輸運方程如下

(9)

在轉(zhuǎn)捩模型的源項中,通過式(10)觸發(fā)該模型

Reθc(Tu,λθ)=CTU1+CTU2exp[-CTU3TuFPG(λθ)]

(10)

式中:Reθc表示臨界動量厚度雷諾數(shù);Tu和λθ分別表示逼近自由湍流強度和壓力梯度的局部變量;FPG的定義參見文獻[7]。

1.3 改進的SST k-ω湍流模型

如前所述,SSTk-ω模型一方面結合了標準k-ω模型和k-ε模型的優(yōu)點,另一方面增加了對逆壓梯度的修正項,如式(4)所示。逆壓梯度修正被認為是SSTk-ω模型能準確預測沖擊射流傳熱和流動的原因之一[19]。然而,SSTk-ω模型中式(5)的應力項S在滯止點附近往往過大,使得湍動能產(chǎn)生項Gk偏大,導致SSTk-ω模型難以預測有渦脫落的流動[8],同時在滯止點處產(chǎn)生過高的傳熱率。對于Kato-Launde模型,從式(7)和式(8)中可以看出其降低了湍動能的產(chǎn)生,特別是在滯止區(qū)流動趨于無旋時,Ω將趨于0。

盡管如此,Wienand等同時對SSTk-ω模型以及Kato-Launder模型進行了研究,發(fā)現(xiàn)努塞爾數(shù)第二峰值依然難以被準確地捕捉[9]。由于第二峰值與轉(zhuǎn)捩相關,因此間歇轉(zhuǎn)捩模型的引入有可能提高模型對轉(zhuǎn)捩的預測能力?？紤]到?jīng)_擊射流涉及到強逆壓梯度、流動分離和渦的形成與發(fā)展,因此將上述3種模型進行結合可能會提高其對沖擊射流傳熱的預測能力?；谝陨嫌懻?將Kato-Launder模型和間歇轉(zhuǎn)捩模型與式(5)進行耦合,如下所示

(11)

(12)

(13)

1.4 計算模型及求解設置

計算模型如圖1所示,其中噴嘴的寬度為B,長度L=50B。由于計算模型關于y軸對稱,為提高計算效率,建立1/2模型進行計算。

圖1 計算模型和邊界條件

本文計算選擇基于壓力基的穩(wěn)態(tài)求解器,即SIMPLE算法。梯度離散格式采用LeastSquaresCellBased,壓力離散格式采用PRESTO,動量和湍動能等采用二階迎風格式。

1.5 網(wǎng)格劃分及無關性驗證

圖2 網(wǎng)格無關性驗證

從圖2中可以看出,4種網(wǎng)格計算得到的努塞爾數(shù)在滯止點處基本重合,其中y+=0.162和y+=0.14的網(wǎng)格在x/B=3處的相對誤差為3%,在x/B=7.4處的相對誤差為1.1%,兩種網(wǎng)格的努塞爾數(shù)分布曲線基本重合,可認為此時的解已滿足網(wǎng)格無關性要求。因此,在沒有特別聲明的情況下,后續(xù)研究選取y+=0.162作為H/B=4時的計算網(wǎng)格。

2 結果分析與討論

圖3 Re=20 000、H/B=4時數(shù)值計算結果與 實驗數(shù)據(jù)的對比

2.1 溫度場分析

圖3給出了H/B=4時采用不同模型計算得到的努塞爾數(shù)與實驗數(shù)據(jù)的分布情況。從圖中可以看出,在x/B=3的位置努塞爾數(shù)達到極小值,隨后傳熱率逐漸增大,在x/B=7處達到第二峰值。介于二次峰值和極小值之間的區(qū)域稱為轉(zhuǎn)捩區(qū)(3≤x/B≤7.4)。在滯止點處,傳熱率與湍流強度緊密相關,較高的湍流強度引起湍動能的過快增長,導致傳熱率增加。此外,來自上游的湍動能傳遞到下游,從而引起轉(zhuǎn)捩的發(fā)生。從圖3中可以看到,在滯止區(qū),圖中所列的模型基本都準確地預測了努塞爾數(shù)。這是因為在基于k-ω的模型中,通過限制渦黏系數(shù)(式(4))來限制雷諾應力過快增長的方法被廣泛地使用,保證了該區(qū)域處于較低湍流度的狀態(tài)[17],從而達到了限制湍動能的目的。同時,改進的SSTk-ω模型由于Kato-Launder模型的作用,預測的努塞爾數(shù)比RANSk-ω模型和k-ω-v2-f模型要低,但更符合實驗結果。在轉(zhuǎn)捩區(qū),RANS/LES模型、k-ω-v2-f模型和本文改進的SSTk-ω模型都準確地捕獲了轉(zhuǎn)捩點(x/B=3)。然而,在第二峰值的預測上,圖3中各模型表現(xiàn)出了明顯的差異,其中RANSk-ω模型、k-ω-v2-f模型和RANS/LES模型均無法準確地捕捉到努塞爾數(shù)的第二峰值的大小和位置,而對于改進的SSTk-ω模型,由于轉(zhuǎn)捩模型的加入,以及在滯止點處準確地預測了第一峰值,從而準確地捕捉到了轉(zhuǎn)捩點。這說明改進的SSTk-ω模型準確地預測了一個應力場,在應力作用下,湍流強度逐漸增大從而導致轉(zhuǎn)捩的發(fā)生。同時,在2≤H/B≤4的范圍內(nèi),RANSk-ω模型獲得的傳熱率高于改進的SSTk-ω模型,可能的解釋是改進的SSTk-ω模型在式(4)中增加了逆壓梯度修正。這是因為在逆壓梯度層中,湍動能的產(chǎn)生大于耗散[23],如果不加以修正,則其預測的湍動能往往過高。因此,如前面分析,這也是SSTk-ω模型被推薦用于沖擊射流傳熱預測的原因之一。在第二峰值附近,受渦破碎的影響,Kato-Launder模型中的渦量隨著渦的破碎而增大,湍動能的生成項變大,從而引起傳熱率增大。改進的SSTk-ω模型在滯止點處準確地預測了努塞爾數(shù),并且成功預測了努塞爾數(shù)第二峰值的位置及其分布;與實驗得到的努塞爾數(shù)最小值位置相比,計算得到的最小值位置在x/B=2.77處,相對誤差為7.6%,第二峰值的位置為x/B=7.30,相對誤差為1.35%。由于RANSk-ω模型對轉(zhuǎn)捩預測的能力較弱,因此基于RANSk-ω模型的RANS/LES模型和k-ω-v2-f模型均是在靠近壁面附近采用RANSk-ω模型,在邊界層外采用LES模型或者v2-f模型,因此這種在不改變黏性層的情況下對模型進行的改進,很難準確地預測努塞爾數(shù)第二峰值現(xiàn)象。

2.2 速度場分析

(a)x/B=1

(b)x/B=2

(c)x/B=5

(d)x/B=7圖4 Re=20 000、H/B=4情況下x方向的歸一化平均速度沿y方向的分布

為了盡可能地考慮實驗的不確定性因素,將數(shù)值計算結果與Ashforth等[20]和Zhe等[21]的實驗數(shù)據(jù)均作了對比。圖4給出了x方向的歸一化平均速度U/Vin沿y方向的分布情況。在滯止區(qū)附近,由于流線的聚集,產(chǎn)生了一個較大的應力場[21]和較低的湍流強度,離滯止點越遠,應力越小。因此,流動方向的法向速度梯度從一個極大值向出口方向逐漸降低,如圖4所示。由于應力的影響,湍流強度隨著x/B的增加逐漸變大,從而促使了轉(zhuǎn)捩的發(fā)生。在x/B=5和x/B=7兩個位置,可以明顯看到減速區(qū)的存在。減速區(qū)的形成受渦發(fā)展的影響[24]。在x/B=1,2的位置,從前面的分析中已知,所有的模型在滯止區(qū)附近產(chǎn)生了一個較低的湍流強度,因此均獲得了較符合實驗的結果,這也同時反映在了滯止點附近對努塞爾數(shù)的預測上。隨著湍流強度的增大,在轉(zhuǎn)捩區(qū)x/B=5的位置,改進的SSTk-ω模型表現(xiàn)得最好。在x/B=7的位置,即渦結構附近[24],改進的SSTk-ω模型過高地預測了速度分布,但是仍然比傳統(tǒng)的RANSk-ω模型和k-ω-v2-f模型更接近于實驗結果,此時只有采用RANS/LES混合模型才能得到與實驗較一致的結果。這是因為RANS/LES模型結合了LES模型在處理渦方面的優(yōu)點,因此在該位置預測的速度也最為準確。Had?iabdi等提供了另外一種可能的解釋,他們通過LES模型對二次峰值展開研究發(fā)現(xiàn),渦和壁面間的相互作用引起的流動分離造成了二次峰值的形成[25]。流動分離過程中伴隨著逆壓梯度的形成,壓力增大,導致速度減小,可能的原因是SSTk-ω模型對逆壓梯度并不是特別敏感[23],因此造成了改進的SSTk-ω模型過高地預測了速度的分布。然而,總體上看,考慮到模型的計算效率和計算精度,本文所用的湍流模型仍然有一定的優(yōu)勢。

2.3 不同沖擊距離的影響

如前所述,努塞爾數(shù)第二峰值歸因于層流到湍流之間的轉(zhuǎn)捩。然而,Had?iabdi等認為流動的分離才是努塞爾數(shù)第二峰值產(chǎn)生的原因[25]。結合以上兩種觀點,說明SSTk-ω模型和轉(zhuǎn)捩模型相結合,能較好地預測努塞爾數(shù)第二峰值的特征及其分布,但是對于速度場的預測,在第二峰值附近,預測的速度依然偏大。為了進一步研究流動分離對努塞爾數(shù)第二峰值的影響,在前面的基礎上,選取H/B=2,4,9.2進行計算。對應的雷諾數(shù)為20 000。圖5顯示了不同沖擊距離下努塞爾數(shù)的分布情況。從圖中可以看出,在H/B=9.2時,努塞爾數(shù)第二峰值消失,在H/B≤4時,努塞爾數(shù)第二峰值較為明顯。

圖5 不同沖擊距離下的努塞爾數(shù)分布

(a)H/B=2

(b)H/B=4

(c)H/B=9.2圖6 不同沖擊距離下沿沖擊面的壓力分布

圖6顯示了不同沖擊距離下壓力沿沖擊面的分布情況,其中P表示壓力,Pmax表示滯止點處的壓力。從圖中可以看出:當H/B≤4時,在下游區(qū)域,壓力逐漸增大,預示著逆壓梯度區(qū)的存在,從而導致速度降低;當H/B=9.2時,下游壓力幾乎不變,相對應地,努塞爾數(shù)第二峰值也消失了。在滯止點附近,壓力急劇地減小,由此形成了速度加速區(qū),如圖4所示。因此,可以認為逆壓梯度和努塞爾數(shù)的第二峰值也存在著一定的聯(lián)系,當下游區(qū)的逆壓梯度消失,不存在流動分離時,努塞爾數(shù)第二峰值也會隨著消失。

3 結 論

本文基于RANS方法,提出了一種改進的SSTk-ω模型計算了H/B=4時帶有封閉板的平板沖擊射流問題。通過將計算得到的溫度分布和速度分布與實驗數(shù)據(jù)對比,在二維平板沖擊射流傳熱問題方面驗證了改進模型的有效性,在此基礎上研究了不同沖擊距離下壓力梯度對努塞爾數(shù)分布的影響。通過本文的研究工作,得到以下主要結論。

(1)改進的SSTk-ω模型與其他模型相比能夠更準確地預測努塞爾數(shù)的分布以及第二峰值的位置。

(2)在速度預測方面,與RANSk-ω和k-ω-v2-f模型過高地預測速度分布相比,帶有附加模型的SSTk-ω模型給出了與RANS/LES模型基本相當,更符合實驗的結果。

(3)當H/B≤4時,努塞爾數(shù)分布表現(xiàn)出明顯的第二峰值的特征,此時在二次峰值附近存在著逆壓梯度區(qū);隨著沖擊距離的增大,逆壓梯度區(qū)消失,努塞爾數(shù)第二峰值也隨著消失。