蔡林峰, 李昊歌, 趙文文, 陳偉芳

(浙江大學 航空航天學院, 浙江 杭州 310027)

0 引 言

對于高超聲速飛行器而言，當流動發(fā)生轉(zhuǎn)捩時飛行器表面熱流和摩阻急劇增加。因此準確預估轉(zhuǎn)捩起始點位置、轉(zhuǎn)捩區(qū)長度，對于高超聲速飛行器工程設(shè)計，尤其是熱防護設(shè)計,具有重要意義，是目前航空航天工程中的關(guān)鍵問題。

采用轉(zhuǎn)捩模型的RANS方法計算高效，是工程實踐中開展流動轉(zhuǎn)捩研究較為經(jīng)濟可行的技術(shù)手段，其中以Walters和Cokljat[1]、王亮和符松[2-3]以及Lantry和Menter[4]等人的研究最具代表性。Lantry等[4]提出的γ-Reθ t轉(zhuǎn)捩模型完全基于當?shù)鼗兞?，能夠和現(xiàn)代CFD程序相兼容，具備三維復雜外形轉(zhuǎn)捩流動模擬能力，在工程中得到了廣泛的應(yīng)用。針對機翼、平板邊界層流動的眾多研究[4-6]表明，γ-Reθ t轉(zhuǎn)捩模型能夠準確捕捉各種類型的低流速邊界層轉(zhuǎn)捩，但對高超聲速邊界層流動轉(zhuǎn)捩現(xiàn)象預測能力有限。為此，Kaynak等[7]基于來流馬赫數(shù)對γ-Reθ t轉(zhuǎn)捩模型經(jīng)驗關(guān)系式進行了可壓縮修正，并針對超聲速平板開展了數(shù)值驗證工作。Gary等[8]對馬赫數(shù)8尖錐進行了數(shù)值計算，研究了來流馬赫數(shù)對γ-Reθ t模型中壓力梯度參數(shù)的影響。國內(nèi)孔維萱等[9]采用γ-Reθ t轉(zhuǎn)捩模型對馬赫數(shù)3～7范圍內(nèi)的高超聲速尖錐流動進行了數(shù)值計算，驗證了該模型具備高超聲速邊界層流動轉(zhuǎn)捩模擬預測能力。張曉東等[10]給出了一種基于當?shù)伛R赫數(shù)的γ-Reθ t轉(zhuǎn)捩模型可壓縮修正關(guān)系式，并在高超聲速雙楔和平板算例中的得到了驗證。張毅峰等[11]通過壓力梯度參數(shù)和普朗特數(shù)對γ-Reθ t轉(zhuǎn)捩模型進行了修正，并在多個高超聲速尖錐轉(zhuǎn)捩流動的模擬中都得到了不錯的結(jié)果。鄭赟等[12]參考高超聲速風洞實驗數(shù)據(jù)對γ-Reθ t轉(zhuǎn)捩模型提出了改進，通過數(shù)值計算分析了原始γ-Reθ t轉(zhuǎn)捩模型對高超聲速邊界層轉(zhuǎn)捩流動預測結(jié)果與實驗值不吻合的原因。周玲等[13]則在符松和王亮[2-3]提出的轉(zhuǎn)捩模式上進行了壁面溫度影響修正，并針對高超聲速飛行器前體邊界層強制轉(zhuǎn)捩開展了數(shù)值研究。

本文作者所在課題組[14]在γ-Reθ t轉(zhuǎn)捩模型的基礎(chǔ)上，通過去掉轉(zhuǎn)捩動量厚度雷諾數(shù)輸運方程，提出了簡化三方程轉(zhuǎn)捩模型。該模型經(jīng)驗關(guān)系式簡單，對低速機翼、高超聲速平板等邊界層流動有良好的模擬預測能力。但在對高超聲速雙楔模型的數(shù)值模擬中發(fā)現(xiàn)，壁面熱流計算結(jié)果與實驗值存在一定差距。因此，本文將在該模型的基礎(chǔ)上引入湍流模型和轉(zhuǎn)捩模型的可壓縮修正方法，研究可壓縮修正對簡化三方程轉(zhuǎn)捩模型預測高超聲速流動轉(zhuǎn)捩的影響，為該模型的進一步改進提供參考。

1 數(shù)值計算方法

1.1 控制方程

控制方程由三維可壓縮Navier-Stokes方程、k-ωSST湍流模型方程和間歇因子輸運方程組成。方程基于結(jié)構(gòu)網(wǎng)格有限體積方法求解，時間離散采用計算效率更高的LU-SGS隱式方法，空間對流通量采用AUSMPW+格式結(jié)合MUSCL差值方法進行求解。

1.2 簡化三方程轉(zhuǎn)捩模型

在原始γ-Reθ t轉(zhuǎn)捩模型中，動量厚度雷諾數(shù)輸運方程的作用是將邊界層外的信息傳播至邊界層內(nèi)，并將轉(zhuǎn)捩判據(jù)、經(jīng)驗關(guān)系式和間歇因子聯(lián)系起來。Coder和Maughmer[15]的研究表明，通過定義一個當?shù)鼗膲毫μ荻葏?shù)Hc=yΩ/U來構(gòu)造新的經(jīng)驗關(guān)系式，可以去掉轉(zhuǎn)捩動量厚度雷諾數(shù)輸運方程并取得較為合理的模擬結(jié)果。類似的，定義參數(shù)Tw=RTΩ/ω于表征剪切力和壓力梯度大小，從而構(gòu)建一個新的臨界動量厚度雷諾數(shù)經(jīng)驗關(guān)系式：

Reθ c=900.0-810.0×

(1)

同時，定義轉(zhuǎn)捩區(qū)域控制函數(shù)為：

Flength=max(0.1,30.0lnTu∞+89.97)

(2)

式(1)和式(2)中系數(shù)由平板轉(zhuǎn)捩流動擬合得到。綜上，耦合間歇因子γ的輸運方程和k-ωSST湍流模型方程可得完整的簡化三方程轉(zhuǎn)捩模型為：

(3)

上式中Pγ和Eγ分別表示間歇因子輸運方程的生成項和破壞項，其具體形式有：

Pγ=Flengthca1ρS(γFonset)0.5(1-ce1γ)

(4)

Eγ=ca2γρΩFturb(1-ce2γ)

(5)

最終通過間歇因子γ修正湍動能輸運方程來耦合簡化轉(zhuǎn)捩模型與k-ωSST湍流模型，可得修正后的湍動能輸運方程生成項和破壞項分別為：

(6)

(7)

為有效模擬流動分離，間歇因子γ的修正公式如下：

Freattach,2)Fθ t

γeff=max(γ,γsep)

(8)

相比于原始的γ-Reθ t轉(zhuǎn)捩模型，簡化三方程轉(zhuǎn)捩模型(算例中用k-ω-γ表示)減少了一個輸運方程，且經(jīng)驗關(guān)系式更為簡單，無需迭代求解。后續(xù)將對該模型進行可壓縮修正，進而開展高超聲速轉(zhuǎn)捩流動的數(shù)值研究。

2 壓縮性修正

2.1 湍流模型壓縮性修正

針對高超聲速湍流流動，國內(nèi)外學者提出了多種不同形式的k-ωSST湍流模型可壓縮修正方法(算例中湍流模型可壓縮修正用TUCC表示)。本文將采用Wilcox等人[16]提出的可壓縮修正方法對湍流模型進行修正，從而提高k-ωSST湍流模型對高超聲速湍流流動的模擬能力。修正形式如下：

β*=β*[1+ξ*F(Mt)]

(9)

β=β-β*ξ*F(Mt)

(10)

其中：

ξ*=2.0,Mt0=0.25

2.2 簡化轉(zhuǎn)捩模型壓縮性修正

簡化轉(zhuǎn)捩模型的臨界動量厚度雷諾數(shù)經(jīng)驗公式是根據(jù)平板實驗數(shù)據(jù)擬合得到的，在應(yīng)用于高超聲速可壓縮流轉(zhuǎn)捩問題時需要進行適當?shù)目蓧嚎s修正。因此，文獻[10] 給出了一種基于來流馬赫數(shù)的可壓縮修正方法(算例中記為TRCC)。該方法定義了一個來流馬赫數(shù)函數(shù)來對臨界動量厚度雷諾數(shù)進行修正，具體修正形式如下：

F(Ma∞)=(1+0.38)0.5

(11)

3 數(shù)值計算結(jié)果與分析

3.1 高超聲速平板流動

算例取自文獻[17]，高超聲速平板長1.5 m，計算來流條件見表1。計算網(wǎng)格224×150，壁面第一層網(wǎng)格y+小于1。

表1 來流條件

圖1至圖3分別給出了三種不同計算來流條件下高超聲速平板壁面斯坦頓數(shù)分布。可以看到，完全層流和完全湍流模型計算結(jié)果與實驗值相差較大，無法捕捉流動轉(zhuǎn)捩現(xiàn)象。增加簡化三方程轉(zhuǎn)捩模型后，計算結(jié)果能夠?qū)崿F(xiàn)層流到湍流的轉(zhuǎn)捩過渡，并且轉(zhuǎn)捩起始位置、轉(zhuǎn)捩區(qū)長度、轉(zhuǎn)捩前后平板壁面斯坦頓數(shù)都能夠與實驗數(shù)據(jù)以及文獻[18]中的計算結(jié)果相吻合?？梢姾喕匠剔D(zhuǎn)捩模型對于簡單平板高超聲速轉(zhuǎn)捩流動具有不錯的模擬能力。此外，從圖1至圖3中可以看到，添加可壓縮修正后，雖然對高超聲速平板轉(zhuǎn)捩位置等無明顯影響，但壁面斯坦頓數(shù)在轉(zhuǎn)捩后湍流區(qū)域有所下降，相比實驗結(jié)果偏低，后文將通過復雜幾何模型轉(zhuǎn)捩流動的數(shù)值計算對該模型及其可壓縮修正在高超聲速轉(zhuǎn)捩流動中的應(yīng)用開展進一步的研究。

圖1 計算條件1時的高超聲速平板壁面斯坦頓數(shù)分布圖

圖2 計算條件2時的高超聲速平板壁面斯坦頓數(shù)分布圖

圖3 計算條件3時的高超聲速平板壁面斯坦頓數(shù)分布圖

3.2 高超聲速雙楔流動

雙楔模型外形如圖4所示，前緣鈍化半徑0.5 mm，其中兩個坡的角度分別為9°和20.5°。來流馬赫數(shù)為8.1，靜壓為520 Pa，靜溫為106 K，來流單位雷諾數(shù)為3.8×106/m，來流湍流度0.9%，壁面溫度為300 K。計算網(wǎng)格見圖5，壁面第一層網(wǎng)格y+小于1。

圖4 雙楔模型幾何外形

圖5 雙楔模型計算網(wǎng)格示意圖

圖6是采用原始簡化三方程轉(zhuǎn)捩模型計算得到的流場馬赫數(shù)云圖，流場結(jié)構(gòu)與實驗結(jié)果基本一致。拐角處可以看到明顯的流動分離和再附，分離激波和再附激波在拐角后交匯，之后再與頭部誘導的弓形脫體激波相交。

圖6 雙楔模型馬赫數(shù)分布云圖

雙楔模型壁面斯坦頓數(shù)、壓力系數(shù)分布圖見圖7和圖8。從圖中可以看出簡化三方程轉(zhuǎn)捩模型相比完全層流模型和完全湍流模型,計算結(jié)果與實驗值更為接近。且從計算結(jié)果可以看出，在拐角處流動分離轉(zhuǎn)捩為湍流，壁面斯坦頓數(shù)明顯增大，說明流動轉(zhuǎn)捩為湍流后壁面熱流急劇增大。但是，對比實驗結(jié)果也可以看到,流動發(fā)生轉(zhuǎn)捩后，原始簡化三方程轉(zhuǎn)捩模型計算得到的壁面斯坦頓數(shù)相比實驗值明顯偏高，與實驗值存在一定誤差。對比圖7中添加湍流模型可壓縮修正后的結(jié)果可以看到，流動轉(zhuǎn)捩后壁面斯坦頓數(shù)計算結(jié)果相比實驗值偏高，主要是由于原始的k-ωSST兩方程湍流模型不適用于高超聲速湍流流動造成的。對k-ωSST兩方程湍流模型增加可壓縮修正后湍流區(qū)熱流計算結(jié)果得到了明顯改進。但圖8中湍流模型可壓縮修正添加前后壓力系數(shù)分布曲線無明顯變化，可以看出，湍流模型可壓縮修正對壁面壓力系數(shù)分布無影響?？梢姡牧髂Ｐ偷目蓧嚎s修正方法能夠有效降低雙楔模型湍流區(qū)壁面熱流的模擬結(jié)果，使之與實驗結(jié)果更為吻合。相應(yīng)的，添加簡化轉(zhuǎn)捩模型可壓縮修正后，模型能夠更好地描述流動的壓縮性，計算得到的轉(zhuǎn)捩區(qū)域增大，雙楔模型壁面斯坦頓數(shù)、壓力系數(shù)與實驗值更為吻合，具有不錯的修正效果。

圖7 雙楔模型壁面斯坦頓數(shù)分布圖

圖8 雙楔模型壁面壓力系數(shù)分布圖

3.3 高超聲速尖錐流動

高超聲速尖錐算例取自文獻[19]。如圖9所示，尖錐長2.35 m，半錐角為7°。計算來流馬赫數(shù)7.15，靜壓7722 Pa，靜溫214 K，來流單位雷諾數(shù)9.64×106/m，湍流度0.06%，壁面溫度300 K。采用半模計算以減少網(wǎng)格量，尖錐對稱面上網(wǎng)格如圖10所示，其中壁面第一層網(wǎng)格間距為1×10-6m。

圖9 尖錐外形

圖10 尖錐對稱面網(wǎng)格示意圖

由圖11可以看到，相比于全層流模型和全湍流模型，簡化三方程轉(zhuǎn)捩模型能夠準確捕捉流動轉(zhuǎn)捩現(xiàn)象，在0.67 m位置處流動發(fā)生轉(zhuǎn)捩，與實驗結(jié)果相吻合。此外，從圖11中可以看到湍流模型和轉(zhuǎn)捩模型的可壓縮修正對轉(zhuǎn)捩位置沒有影響。其中湍流模型可壓縮修正主要用于確保k-ωSST湍流模型對高超聲速湍流流動模擬的準確性，降低了流區(qū)的熱流計算值。而轉(zhuǎn)捩模型可壓縮修正主要用于確保對轉(zhuǎn)捩區(qū)流動壓縮性的準確描述，從而保證轉(zhuǎn)捩區(qū)長度的準確計算，最終復現(xiàn)了實驗結(jié)果。

圖11 尖錐壁面熱流分布

4 結(jié) 論

通過對湍流模型和簡化轉(zhuǎn)捩模型進行可壓縮修正，研究了簡化三方程轉(zhuǎn)捩模型在高超聲速轉(zhuǎn)捩流動中的應(yīng)用。對比高超聲速平板、雙楔、尖錐的數(shù)值計算與實驗結(jié)果，得出以下結(jié)論：

1)簡化三方程轉(zhuǎn)捩模型對于高超聲速轉(zhuǎn)捩流動具有模擬預測能力，能夠較好地復現(xiàn)流場結(jié)構(gòu)，準確捕捉轉(zhuǎn)捩起始位置，但對壁面熱流和轉(zhuǎn)捩區(qū)長度的模擬結(jié)果與實驗數(shù)據(jù)存在差距，需要引入合適的可壓縮修正方法來提高模型計算精度。

2)湍流模型可壓縮修正方法主要用于修正湍流區(qū)的壁面熱流模擬結(jié)果，而簡化轉(zhuǎn)捩模型臨界動量厚度雷諾數(shù)經(jīng)驗關(guān)系式對轉(zhuǎn)捩區(qū)長度影響較大，添加可壓縮修正能夠改進模型對高超聲速流動轉(zhuǎn)捩的模擬能力。

3)簡化三方程轉(zhuǎn)捩模型經(jīng)驗關(guān)系式由平板轉(zhuǎn)捩數(shù)據(jù)擬合得到，對于雙楔、尖錐模型需要添加修正后才能夠更為準確地預測壁面熱流、轉(zhuǎn)捩區(qū)長度。計算結(jié)果與計算外形相關(guān)，發(fā)展更具普適性的高超速邊界層轉(zhuǎn)捩預測方法仍然需要更深入的探索研究。