考慮轉捩的SST模型在沖擊射流傳熱中的應用

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

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

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

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

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

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

1 數值計算方法

1.1 基本控制方程

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

(1)

(3)

1.2 SST k-ω湍流模型

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

(4)

(5)

Gω-Yω+Dω+Sω

(6)

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

Gk=μtS2

(7)

Gk=μtSΩ

(8)

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

(9)

在轉捩模型的源項中,通過式(10)觸發該模型

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

(10)

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

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

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

盡管如此,Wienand等同時對SSTk-ω模型以及Kato-Launder模型進行了研究,發現努塞爾數第二峰值依然難以被準確地捕捉[9]。由于第二峰值與轉捩相關,因此間歇轉捩模型的引入有可能提高模型對轉捩的預測能力。考慮到沖擊射流涉及到強逆壓梯度、流動分離和渦的形成與發展,因此將上述3種模型進行結合可能會提高其對沖擊射流傳熱的預測能力。基于以上討論,將Kato-Launder模型和間歇轉捩模型與式(5)進行耦合,如下所示

(11)

(12)

(13)

1.4 計算模型及求解設置

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

圖1 計算模型和邊界條件

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

1.5 網格劃分及無關性驗證

圖2 網格無關性驗證

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

2 結果分析與討論

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

2.1 溫度場分析

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

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方向的分布

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

2.3 不同沖擊距離的影響

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

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

(a)H/B=2

(b)H/B=4

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

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

3 結 論

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

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

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

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