高超聲速飛行器復(fù)雜外形轉(zhuǎn)捩預(yù)測

李 齊, 董 穎, 趙雅甜, 趙 瑞

(1. 北京空間飛行器總體設(shè)計部, 北京 100094; 2. 北京理工大學(xué), 北京 100081; 3. 中南大學(xué), 湖南長沙 410083)

引 言

在高超聲速飛行器的飛行高度、 速度和Reyno-lds數(shù)范圍內(nèi)非常容易出現(xiàn)邊界層轉(zhuǎn)捩[1-2], 轉(zhuǎn)捩后壁面摩阻和熱流可增大數(shù)倍之多, 嚴(yán)重影響飛行器的氣動性能和熱防護(hù)系統(tǒng)[3-4], 造成壁面燒蝕嚴(yán)重、 顫振加劇、 姿態(tài)控制困難等后果[2]. 因此, 研究高超聲速邊界層轉(zhuǎn)捩問題, 準(zhǔn)確預(yù)測邊界層的轉(zhuǎn)捩起始位置, 對飛行器的氣動設(shè)計和熱防護(hù)設(shè)計至關(guān)重要.

近幾十年來人們對邊界層轉(zhuǎn)捩預(yù)測開展了全面的探索和研究[5]. 當(dāng)前, 面向高超聲速復(fù)雜外形飛行器的邊界層轉(zhuǎn)捩預(yù)測手段, 主要包括實驗研究、 基于穩(wěn)定性理論的eN方法、 轉(zhuǎn)捩準(zhǔn)則以及模式理論. 實驗研究具有直觀性、 真實性的特點, 面向邊界層轉(zhuǎn)捩的實驗研究主要包括飛行試驗和風(fēng)洞試驗兩大類. 其中具有代表性的飛行試驗工作如: 美國的HyBoLT(Hypersonic Boundary Layer Transition)計劃[6]、 法國的LEA項目[7]、 歐洲的EXPERT(European Experimental Reentry Testbed)項目[8]等. 2015年12月完成的MF-1飛行是我國高超聲速模型飛行試驗歷史上首次以空氣動力學(xué)基礎(chǔ)研究為目的的飛行試驗[9]. 飛行試驗花費極其高昂, 試驗的設(shè)計周期也較長. 為了復(fù)現(xiàn)高空飛行條件下邊界層轉(zhuǎn)捩過程, 實現(xiàn)轉(zhuǎn)捩預(yù)測結(jié)果的天地一致性, 靜音風(fēng)洞應(yīng)運而生. Yao等[10]在Ma=6靜音風(fēng)洞中對平板三角翼的邊界層轉(zhuǎn)捩情況進(jìn)行了研究, 分析了攻角和Reynolds數(shù)對邊界層轉(zhuǎn)捩的影響, 結(jié)果表明, 側(cè)線邊界層的轉(zhuǎn)捩比中線邊界層的轉(zhuǎn)捩早, 大攻角促進(jìn)邊界層的轉(zhuǎn)捩. 隨著Reynolds數(shù)增加, 邊界層的轉(zhuǎn)捩位置向上游移動. Lu等[11]在Ma=6高超聲速低噪聲風(fēng)洞中對45°后掠平板邊界層轉(zhuǎn)捩過程進(jìn)行了研究, 得到了邊界層從層流到湍流的時空演變特征. 劉向宏等[12]總結(jié)了有關(guān)高超聲速邊界層轉(zhuǎn)捩實驗研究.

eN方法基于穩(wěn)定性理論得到主要擾動的增長率, 通過經(jīng)驗參數(shù)N判斷轉(zhuǎn)捩, 其估算轉(zhuǎn)捩位置和范圍的技術(shù)水平取決于對所涉及的物理機(jī)理的了解程度[13]. 轉(zhuǎn)捩準(zhǔn)則是根據(jù)大量地面實驗和飛行試驗數(shù)據(jù)擬合得到的有效經(jīng)驗公式[14]. 歷史上, 針對各類湍流問題突出的飛行器均發(fā)展了相應(yīng)的轉(zhuǎn)捩判據(jù), 如Potter & Whitfield轉(zhuǎn)捩判據(jù)(1962), Van Driest-Blumer轉(zhuǎn)捩判據(jù)(1967), PANT轉(zhuǎn)捩判據(jù)(1975), Boudreau轉(zhuǎn)捩判據(jù)(1981),Reθ/Me判據(jù)(1997), BLT判據(jù)[14]等. 孔維萱等[15]提出“層流+轉(zhuǎn)捩準(zhǔn)則模型”來預(yù)測流向不穩(wěn)定失穩(wěn)誘導(dǎo)形成的邊界層自然轉(zhuǎn)捩, 周玲等[16]在此基礎(chǔ)上增加了橫流轉(zhuǎn)捩準(zhǔn)則判據(jù).

轉(zhuǎn)捩模式是計算復(fù)雜外形轉(zhuǎn)捩的主要方法, 包括基于當(dāng)?shù)刈兞康摩?Reθt模式[17-18], 適用于高超聲速邊界層轉(zhuǎn)捩預(yù)測的k-ω-γ轉(zhuǎn)捩模式[19-20]等.γ-Reθt模式已在低速邊界層轉(zhuǎn)捩流動模擬中獲得了成功應(yīng)用, 近年來有學(xué)者開始將其拓展至高超聲速邊界層轉(zhuǎn)捩模擬[3]. Krause等[21]考慮來流湍流度構(gòu)造了Flength和Reθc關(guān)系式, 較好地模擬了Ma=8.3的雙楔壁面壓力分布. Zhao等[22-23]評估對比了3種經(jīng)典轉(zhuǎn)捩模式(γ-Reθ,k-ω-γ和kT-kL-ω轉(zhuǎn)捩模式)在高超聲速復(fù)雜外形轉(zhuǎn)捩預(yù)測中的性能表現(xiàn). Zhou等[24-26]對k-ω-γ轉(zhuǎn)捩模式進(jìn)行了改進(jìn). 由于改進(jìn)的k-ω-γ轉(zhuǎn)捩模式采用網(wǎng)格預(yù)處理技術(shù), 可以快速大規(guī)模并行獲取具有非當(dāng)?shù)刈兞繉傩缘倪吔鐚訁?shù). 同時, 改進(jìn)的k-ω-γ轉(zhuǎn)捩模式考慮了橫流轉(zhuǎn)捩, 因此具有推廣至工程復(fù)雜外形應(yīng)用的潛力[27]. 曾宏剛等[28]針對高超聲速飛行器前體及壓縮面組合外形, 在Ma=6下采用k-ω-γ轉(zhuǎn)捩模式預(yù)測轉(zhuǎn)捩起始位置, 結(jié)果與風(fēng)洞實驗較為吻合. Zhao等進(jìn)一步對k-ω-γ轉(zhuǎn)捩模式經(jīng)驗參數(shù)[29]和來流擾動[30]導(dǎo)致的預(yù)測不確定性開展了不確定的量化分析工作. 本文轉(zhuǎn)捩模式預(yù)測中即采用k-ω-γ轉(zhuǎn)捩模式.

目前, 采用轉(zhuǎn)捩模式預(yù)測工程復(fù)雜外形, 尤其是高超聲速流動中復(fù)雜外形的轉(zhuǎn)捩研究較少. 本文基于改進(jìn)的k-ω-γ轉(zhuǎn)捩模式對高超聲速飛行器復(fù)雜主體外形進(jìn)行了流場分析和轉(zhuǎn)捩預(yù)測, 分析了Reynolds數(shù)和攻角對邊界層轉(zhuǎn)捩的影響, 為類似的工程應(yīng)用提供參考.

1 數(shù)值方法

1.1 計算模型

該飛行器采用上反翼布局, 進(jìn)氣方式為背部進(jìn)氣, 飛行器總長45 m, 外形布局關(guān)鍵幾何參數(shù)詳見文獻(xiàn)［31-32］. 由于計算工況無側(cè)滑角, 因此采用半模計算. 計算網(wǎng)格為多區(qū)對接結(jié)構(gòu)網(wǎng)格. 飛行器半?？偩W(wǎng)格量6.7×106, 如圖 1所示. 壁面第1層網(wǎng)格高度為0.001 mm, 保證y+≤1.

(a) Global grid topology

(b) Symmetry plane grids

(c) Wall grids圖1 飛行器全機(jī)計算網(wǎng)格示意圖Fig. 1 Computational zone and grid distribution of the configuration

1.2 轉(zhuǎn)捩模式

本研究使用的k-ω-γ轉(zhuǎn)捩模式以SST湍流模式為基礎(chǔ), 同時引入γ輸運方程, 并進(jìn)行了Reynolds數(shù)和橫流轉(zhuǎn)捩的預(yù)測修正. 該k-ω-γ模式的總體框架為

式中,Pγ和Dγ分別為γ輸運方程的生成項和耗散項, 基于量綱分析構(gòu)造, 具體表達(dá)式為

本文選取有限體積方法求解控制方程. 層流N-S方程和湍流、 轉(zhuǎn)捩輸運方程采用統(tǒng)一的離散方式: 無黏通量采用差分分裂的Roe格式, 網(wǎng)格界面上采用MUSCL(monotone upstream-centered schemes for conservation laws)格式插值獲得高階精度; 黏性項采用2階中心差分格式進(jìn)行離散; 時間推進(jìn)采用隱式LU-SGS(lower-upper symmetric Gauss-Seidel)格式. Zhou等[26]將改進(jìn)的k-ω-γ模式用于預(yù)測X-51A飛行器前體的邊界層轉(zhuǎn)捩, 得到了和風(fēng)洞實驗一致的結(jié)果. 同時與Reθ/Me轉(zhuǎn)捩準(zhǔn)則的預(yù)測結(jié)果進(jìn)行了對比, 結(jié)果較為吻合, 因此本文不再贅述算法驗證, 計算方法相關(guān)驗證細(xì)節(jié)可參考文獻(xiàn)[26].

計算的來流條件為: 飛行高度分別為H=25， 35 km, 飛行Mach數(shù)為Ma=6, 攻角分別為α=0°， 5°, 壁面溫度為Tw=300 K, 不同飛行高度對應(yīng)的單位Reynolds數(shù)分別為Re=4.9×106， 1.0×106/m, 以研究不同Reynolds數(shù)和攻角對飛行器邊界層轉(zhuǎn)捩的影響.

2 流場分析

圖2 對稱面和流向截面等Mach線分布(4個截面位置分別為x/L=0.1,0.3,0.6,0.9, L為機(jī)身全長)Fig. 2 Mach number distribution on the symmetry plane and flow cross sections(The four sections are at x/L=0.1,0.3,0.6, 0.9, respectively, L is the total length of aircraft body)

圖3 流向截面無量綱壓力云圖(4個截面位置分別為x/L=0.1,0.3,0.6,0.9) Fig. 3 Normalized pressure contours at four sections: x/L=0.1, 0.3, 0.6, 0.9

圖4 橫流Reynolds數(shù)Recf和表面極限流線分布Fig. 4 Recf and limit streamline distribution on the vehicle surface

3 轉(zhuǎn)捩預(yù)測結(jié)果分析

采用轉(zhuǎn)捩模式對飛行器在不同Reynolds數(shù)(高度)、 攻角條件下的邊界層轉(zhuǎn)捩位置進(jìn)行預(yù)測. 本文采用的轉(zhuǎn)捩模式考慮了第一模態(tài)、 第二 Mack模態(tài)以及橫流模態(tài), 用于預(yù)測流向不穩(wěn)定和橫流不穩(wěn)定誘導(dǎo)的轉(zhuǎn)捩. 通過飛行器壁面間歇因子的變化判斷飛行器全機(jī)邊界層自然轉(zhuǎn)捩起始位置的分布. 重點分析了飛行器背風(fēng)面、 迎風(fēng)面和側(cè)面轉(zhuǎn)捩區(qū)域的分布; 同時對頭部和舵面兩處局部結(jié)構(gòu)的轉(zhuǎn)捩起始位置進(jìn)行了詳細(xì)分析.

3.1 Reynolds數(shù)對轉(zhuǎn)捩的影響

圖 5為0° 攻角時不同Reynolds數(shù)條件下飛行器全機(jī)壁面間歇因子分布云圖. 轉(zhuǎn)捩模式計算將間歇因子γ=0.2～1的范圍定義為轉(zhuǎn)捩區(qū). 整體而言, 隨著飛行高度的增加, 來流密度減小, 自由來流Reynolds數(shù)逐漸降低, 轉(zhuǎn)捩位置明顯滯后. 注意到, 圖 6中背風(fēng)面中心線附近轉(zhuǎn)捩位置較展向其他位置提前, 呈現(xiàn)明顯的“凸”字型轉(zhuǎn)捩陣型, 對比圖 4表面極限流線分布可以看到, 相應(yīng)位置處發(fā)生了流動分離, 因此“凸”字型轉(zhuǎn)捩陣型是由于進(jìn)氣道壓縮面逆壓梯度誘導(dǎo)產(chǎn)生小的流動分離, 流動分離誘導(dǎo)了提前轉(zhuǎn)捩. 進(jìn)氣道噴管區(qū)域由于膨脹加速效應(yīng), 在順壓梯度的作用下使得流動再層流化. 機(jī)翼前緣橫流效應(yīng)并不顯著(見圖 4), 該處轉(zhuǎn)捩主要由流向不穩(wěn)定模態(tài)主導(dǎo)(見圖 7).

(a) Re=4.9×106/m

(b) Re=1.0×106/m圖5 不同Reynolds數(shù)下壁面間歇因子分布云圖 (α=0°)Fig. 5 γ distribution on the vehicle wall at different Reynolds numbers(α=0°)

(a) Re=4.9×106/m

(b) Re=1.0×106/m圖6 不同Reynolds數(shù)下頭部壁面間歇因子分布云圖(α=0°)Fig. 6 γ distribution on the vehicle head wall at different Reynolds numbers(α=0°)

圖7 不同Reynolds數(shù)下舵面間歇因子分布云圖(α=0°)Fig. 7 γ distribution on the rudder wall at different Reynolds numbers(α=0°)

橫流Reynolds數(shù)Recf常用來判斷流場中橫流區(qū)域分布及其強(qiáng)度.Recf的表達(dá)式為Recf=δ0.1wmax/ve, 其中δ0.1是當(dāng)?shù)貦M流速度滿足w/wmax=0.1的壁面法向距離,wmax為邊界層內(nèi)最大橫流速度,ve為邊界層外緣分子運動黏性系數(shù). 橫流Reynolds數(shù)越大, 說明橫流效應(yīng)越強(qiáng). 一般認(rèn)為, 對于可壓縮流動,Recf大于其臨界值200～400即可認(rèn)為橫流誘導(dǎo)轉(zhuǎn)捩發(fā)生. 結(jié)合圖 4可知, 迎風(fēng)面和背風(fēng)面大面積的轉(zhuǎn)捩均由橫流不穩(wěn)定主導(dǎo).

3.2 攻角對轉(zhuǎn)捩的影響

圖 8為5° 攻角,Re=1.0×106/m下飛行器全機(jī)壁面間歇因子分布云圖. 與圖 5(b)對比可知, 隨著攻角增大, 頭部以及機(jī)翼轉(zhuǎn)捩位置向上游移動, 背風(fēng)面噴管區(qū)域?qū)恿髅娣e增加. 圖 9為5° 攻角下頭部和舵面的局部細(xì)節(jié)圖, 經(jīng)定量比較可以看出攻角增加導(dǎo)致頭部迎風(fēng)面和背風(fēng)面轉(zhuǎn)捩起始位置均前移. 舵面轉(zhuǎn)捩受攻角影響不大, 這是因為舵面邊界層轉(zhuǎn)捩由流向不穩(wěn)定主導(dǎo), 對于本文外形, 攻角對流向不穩(wěn)定影響較小. 圖10給出了不同攻角下的對稱面無量綱壓力云圖, 隨著攻角增加, 迎風(fēng)面前緣斜激波更強(qiáng), 背風(fēng)面進(jìn)氣道壓縮激波減弱, 噴管區(qū)域膨脹效應(yīng)增強(qiáng). 飛行器迎風(fēng)面壓力增大, 背風(fēng)面壓力減小, 頭部附近背風(fēng)面與側(cè)邊的展向壓力梯度增大, 橫流效應(yīng)增強(qiáng). 如圖 11所示, 背風(fēng)面流動匯聚線向上游移動, 匯聚線下游的橫向流動更加顯著, 導(dǎo)致轉(zhuǎn)捩提前. 0°攻角下由流動分離誘導(dǎo)產(chǎn)生的“凸”字型轉(zhuǎn)捩型線不再存在, 對比圖 11可知, 5°攻角下流動分離發(fā)生在轉(zhuǎn)捩之后, 故“凸”字型消失. 迎風(fēng)面的轉(zhuǎn)捩起始位置前移則是由攻角的增大導(dǎo)致頭部激波增強(qiáng), 波后密度增加, 如圖 12所示, 邊界層外緣Reynolds數(shù)增大, 導(dǎo)致轉(zhuǎn)捩提前.

圖8 壁面間歇因子分布云圖(α=5°, Re=1.0×106/m)Fig. 8 γ distribution on the vehicle wall(α=5°, Re=1.0×106/m)

(a) Head

(b) Rudder圖9 局部壁面間歇因子的分布云圖(α=5°, Re=1.0×106/m)Fig. 9 γ distribution on the local wall(α=5°, Re=1.0×106/m)

(a) α=0°

(b) α=5°圖10 不同攻角下對稱面無量綱壓力云圖(Re=1.0×106/m)Fig. 10 Pressure distribution on symmetry plane at different angles of attack(Re=1.0×106/m)

圖11 橫流Reynolds數(shù)Recf和表面極限流線(Re=1.0×106/m, α=5°)Fig. 11 Recf and limit streamline distribution on the vehicle wall(Re=1.0×106/m, α=5°)

(a) α=0°

(b) α=5°圖12 不同攻角下對稱面Reynolds數(shù)云圖(Re=1.0×106/m)Fig. 12 Reynolds number distribution on symmetry plane at different angles of attack(Re=1.0×106/m)

4 結(jié)論

本文以高超聲速飛行器復(fù)雜外形為研究對象, 采用改進(jìn)的k-ω-γ轉(zhuǎn)捩模式對主體外形進(jìn)行了流場分析和邊界層轉(zhuǎn)捩預(yù)測, 并分析了Reynolds數(shù)和攻角對邊界層轉(zhuǎn)捩的影響. 研究得到的主要結(jié)論如下:

橫流是影響飛行器大面積轉(zhuǎn)捩的主要因素. 復(fù)雜飛行器外形由于側(cè)緣和機(jī)翼尖前緣激波后壓力較大, 形成展向壓力梯度, 使得迎風(fēng)面和背風(fēng)面大部分區(qū)域轉(zhuǎn)捩主要由橫流主導(dǎo). 隨著高度增加, 來流Reynolds數(shù)減小, 迎風(fēng)面和背風(fēng)面的轉(zhuǎn)捩起始位置均向下游移動. 攻角增大導(dǎo)致頭部激波增強(qiáng), 波后迎風(fēng)面密度增加, 邊界層外緣Reynolds數(shù)增大, 導(dǎo)致轉(zhuǎn)捩提前. 背風(fēng)面尾噴管區(qū)域存在膨脹加速效應(yīng), 在順壓梯度的作用下存在再層流化, 攻角增大后膨脹效應(yīng)增強(qiáng), 層流區(qū)域增加. 0°攻角下由于背風(fēng)側(cè)進(jìn)氣道壓縮面逆壓梯度誘導(dǎo)產(chǎn)生小的流動分離, 流動分離誘導(dǎo)了提前轉(zhuǎn)捩, 形成“凸”字型轉(zhuǎn)捩型線; 5°攻角下背風(fēng)面橫流效應(yīng)增強(qiáng), 流動分離發(fā)生在轉(zhuǎn)捩之后, “凸”字型不再存在.