楊俊沅 李旭東 曾舒華 趙文文,2) 張 賦 陳偉芳

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

? (北京航天長征飛行器研究所,北京 100076)

引言

由于臨近空間高超聲速飛行器有著更快的飛行速度和高空機動能力,逐漸成為了世界各航空航天大國追逐的熱點.對高超聲速飛行器而言,升阻比越高,機動性越好[1].由于尖化前緣這一構(gòu)型能夠為飛行器提供更大的升阻比,高超聲速飛行器的設(shè)計普遍選擇尖前緣外形.然而,相比于鈍頭設(shè)計,高超聲速來流在飛行器的尖化前緣處會產(chǎn)生更劇烈的氣動加熱現(xiàn)象,引發(fā)諸如氣動熱載荷大、熱防護難等難題,給快速準(zhǔn)確的熱環(huán)境預(yù)示提出了更高的要求和挑戰(zhàn).尖化前緣熱環(huán)境的準(zhǔn)確刻畫是尖化前緣熱防護研究的前提[2].黃飛等[3]分別采用納維-斯托克斯(Navier-Stokes,NS)方程和蒙特卡洛直接數(shù)值模擬(direct simulation Monte Carlo,DSMC)方法計算不同努森數(shù)條件下尖前緣氣動特性,發(fā)現(xiàn)駐點熱流密度受局部稀薄效應(yīng)影響較大.姜貴慶等[4]發(fā)現(xiàn)尖化前緣表面熱流密度在離駐點2～3 個分子自由程長度時下降至駐點熱流密度的1/3,而表面溫度此處僅下降至駐點溫度的90%.表明尖化前緣的局部稀薄效應(yīng)對駐點熱流密度的精確預(yù)示產(chǎn)生了較大影響.工程問題中,連續(xù)流條件下駐點熱流密度的工程快速預(yù)測通常可采用Fay-Riddell 公式[5],熱流密度和前緣半徑的關(guān)系可以表示為

式中,qs為尖化前緣的駐點熱流密度,RN表示飛行器的前緣半徑.當(dāng)前緣半徑過小時,駐點處易出現(xiàn)顯著的局部稀薄氣體效應(yīng),基于連續(xù)介質(zhì)假設(shè)和層流邊界層假設(shè)的Fay-Riddell 公式也會失效.理論上前緣半徑趨近于0 時,駐點熱流密度將趨近于無窮大,但實際上熱流密度會趨近于自由分子流的極限[6].

在含局部稀薄效應(yīng)的連續(xù)-稀薄耦合流動區(qū)域,NS 方程線性本構(gòu)模型不足以準(zhǔn)確描述稀薄非平衡和多尺度非線性的流動機理.高超聲速飛行器的飛行環(huán)境介于稠密大氣層和稀薄大氣層之間,飛行環(huán)境不再單一,在尖前緣可能出現(xiàn)顯著的局部稀薄氣體效應(yīng),疊加滑移邊界條件的NS 方程結(jié)果仍需通過實驗數(shù)據(jù)或其他多尺度高精度計算模型進行校準(zhǔn)[7].近年來,針對連續(xù)-稀薄跨流域流動現(xiàn)象的描述,Myong 等[8]在廣義流體動力學(xué)方程基礎(chǔ)上發(fā)展了非線性耦合本構(gòu)關(guān)系(nonlinear coupled constitutive relations,NCCR)模型,并通過經(jīng)典非平衡流動問題驗證了NCCR 模型準(zhǔn)確描述高速多尺度和稀薄非平衡流動現(xiàn)象的能力[9-11].隨后,NCCR 模型得到了國內(nèi)外眾多學(xué)者的深入研究和理論驗證[12-14],在計算框架[15]、求解算法[16-18]、邊界條件[19-20]、熱化學(xué)非平衡流動[21-23]、非結(jié)構(gòu)網(wǎng)格求解[24-25]等方面進展顯著.然而,NCCR 模型在高速連續(xù)/稀薄流實驗驗證方面的研究工作寥寥[26],并且缺少專門針對局部稀薄條件下的尖化前緣氣動加熱機理的影響研究.

本文首先通過圓柱繞流算例驗證NCCR 模型的非平衡多尺度模擬能力;再通過對比等效高度33 km和60 km 條件下NCCR 模型計算的流場溫度和壁面熱流分布與實驗值的偏差,檢驗NCCR 模型在尖化前緣構(gòu)型中準(zhǔn)確描述局部稀薄非平衡流動和物面氣動熱的性能.

1 非線性耦合本構(gòu)關(guān)系模型

非線性耦合本構(gòu)關(guān)系模型是Eu 在廣義流體動力學(xué)方程的基礎(chǔ)上對高階非守恒量的時間項和對流項進行簡化處理得到的[27-28].守恒量和非守恒量的方程可寫為

式中,ρ,u,p,E,T分別為密度、速度、壓力、能量和溫度,Π,?,Q,I,cp分別表示剪切應(yīng)力張量、附加體積應(yīng)力、熱流、單位矩陣和氣體的定壓比熱;η,ηb,λ為剪切黏性系數(shù)、體積黏性系數(shù)和熱傳導(dǎo)系數(shù);是氣體的比熱比,q(κ)是非線性耗散項,非線性耗散項q(κ)的引入不僅增強了方程組的非線性度,同時使非守恒變量都耦合在一起

經(jīng)過Eu 的高階矩封閉以及Eu 和Myong 絕熱假設(shè)簡化[8],可以得到非線性耦合本構(gòu)關(guān)系模型表達式為

為方便比較,NS 方程線性本構(gòu)模型也在這給出

通過對比可以看出,NCCR 模型具有更強的非線性,從而在非平衡流動問題描述中可能具備更高的模擬精度和準(zhǔn)確性.

2 算例驗證

為驗證非線性耦合本構(gòu)關(guān)系模型求解連續(xù)-稀薄跨流域流動問題,本文選取了以壓縮流動為主的高超聲速圓柱繞流算例對NCCR 模型理論體系進行驗證.圓柱算例半徑為6 inch (1 inch=2.54 cm),計算域網(wǎng)格采用結(jié)構(gòu)網(wǎng)格,壁面和徑向方向網(wǎng)格數(shù)目分別為100 和100.氣體黏性使用逆冪律分子模型求解,參考溫度Tref=273K其逆冪律指數(shù)為s=0.74,NCCR常數(shù)為1.02029,氣體常數(shù)為296.8 m2/(s2·K),普朗特常數(shù)Pr=0.72,氣體比熱比 γ=1.4,參考黏性ηref=1.656×10-5N·s/m2.數(shù)值模擬的計算狀態(tài)均與文獻[29]中DSMC 的計算狀態(tài)一致[29].具體來流參數(shù)見表1.

表1 流動狀態(tài)Table 1 Flow state

圖1 給出了Kn分別為0.002,0.01 和0.05 時NCCR 模型、NS 方程和DSMC 模擬得到的駐點線對比曲線.可以看出在連續(xù)流條件下,NCCR 模型計算結(jié)果同NS 方程數(shù)據(jù)保持一致,且兩種方法計算結(jié)果同DSMC 數(shù)據(jù)也保持一致.在非連續(xù)流區(qū)域,NCCR 模型計算結(jié)果同DSMC 數(shù)據(jù)更為吻合、誤差更小,NS 方程結(jié)果與DSMC 模擬結(jié)果的偏差隨流動稀薄程度的增加而增加,流動越稀薄,NS 方程的準(zhǔn)確性越差.通過此算例,驗證了NCCR 模型在稀薄流動領(lǐng)域能獲得比NS 方程更為準(zhǔn)確、更符合物理的數(shù)值解,同時也證明了NS 方程線性本構(gòu)模型在高努森數(shù)稀薄流動領(lǐng)域的不適用性.

圖1 不同努森數(shù)條件下圓柱沿駐點線溫度分布Fig.1 Temperature distribution of cylinders in different Kn along the stagnation line

圖2 給出了Kn=0.002 連續(xù)流條件下圓柱模型壁面摩擦、熱流和壓力系數(shù)分布,NCCR 模型、NS 方程及DSMC 數(shù)據(jù)基本一致.證明在連續(xù)流區(qū)域NCCR 模型計算結(jié)果能夠回歸到NS 方程結(jié)果上.圖3～圖4 為Kn=0.01,0.05 時圓柱模型壁面摩擦、熱流、壓力系數(shù)的DSMC 數(shù)據(jù)、NCCR 模型和NS 方程的結(jié)果對比.從連續(xù)流過渡到滑移流的過程中可以看出,NCCR 模型同NS 方程計算結(jié)果出現(xiàn)差距,且隨著Kn數(shù)的增大,NCCR 模型同NS 方程之間的差距逐漸增大,NCCR 模型計算結(jié)果比NS 方程更貼合DSMC 數(shù)據(jù)結(jié)果.隨著來流稀薄程度的增加,壁面熱流系數(shù)和壁面壓力系數(shù)的變化較小,壁面摩擦系數(shù)的變化最為明顯.此時NCCR 模型與NS 方程計算結(jié)果的差距進一步擴大,從結(jié)果上看,NCCR 模型仍然比NS 方程更為貼近DSMC 的仿真數(shù)據(jù).NCCR 模型的計算結(jié)果同DSMC 的仿真數(shù)據(jù)之間仍有部分差異,這可能是由于NCCR 模型為降低收斂難度而對某些項簡化導(dǎo)致了計算結(jié)果精度的降低,在后續(xù)工作中將進一步開展研究.

圖2 Kn=0.002 時物面系數(shù)分布Fig.2 Surface coefficient distribution when Kn=0.002

圖3 Kn=0.01 時物面系數(shù)分布Fig.3 Surface coefficient distribution when Kn=0.01

圖4 Kn=0.05 時物面系數(shù)分布Fig.4 Surface coefficient distribution when Kn=0.05

3 尖化前緣氣動加熱計算

在連續(xù)流域,前緣的駐點熱流密度可以通過工程中常用的Fay-Riddell 公式進行預(yù)測,但在尖化前緣駐點區(qū)域發(fā)生局部稀薄效應(yīng)時,Fay-Riddell 公式預(yù)測值將偏離真值[30].本文選擇前緣半徑R=0.5 mm 的飛行器,尖化前緣外形如圖5 所示.通過分別對比連續(xù)流和滑移流域NCCR 模型、NS 方程及尖化前緣風(fēng)洞試驗結(jié)果,分析飛行器的尖化前緣區(qū)域是否出現(xiàn)局部稀薄流動,并探究NCCR 模型在計算局部稀薄流動方面的準(zhǔn)確性.

圖5 模型示意圖Fig.5 Model diagram

模型長度為640 mm,尾部寬度為616.78 mm,前緣半徑0.5 mm,尾部高度113.35 mm.

風(fēng)洞試驗在中國科學(xué)院力學(xué)研究所高溫氣體動力學(xué)國家重點實驗室JF8A 和JF10 高超聲速風(fēng)洞中進行.試驗?zāi)Ｐ偷膩砹鳁l件見表2,Ma=8.0 的飛行條件在JF8A 中試驗[31],JF10 風(fēng)洞中進行Ma=10.0 飛行條件試驗[32],后續(xù)計算模型的流動參數(shù)同風(fēng)洞試驗參數(shù)保持一致.

表2 流場狀態(tài)參數(shù)Table 2 Flow field status parameters

3.1 33 km 處高超聲速尖化前緣連續(xù)流動效應(yīng)分析

等效高度33 km 飛行條件下,前緣半徑R=0.5 mm 模型的尖化前緣局部Kn數(shù)為0.007,處于連續(xù)流區(qū)域,計算條件為U=1190.25 m/s,Ma=8.0,T=55.07 K,Tw=300 K,P=798.95 Pa.計算仿真模型各項尺寸同試驗?zāi)Ｐ捅３忠恢?圖6 對比了NCCR模型與NS 方程計算的沿駐點線溫度云圖,兩種計算方法的結(jié)果基本吻合、數(shù)值模擬結(jié)果保持一致.

圖6 駐點附近的溫度對比(Ma=8)Fig.6 Comparison of the temperature near the stagnation point (Ma=8)

首先,為了消除計算網(wǎng)格對計算結(jié)果的影響,采用3 套不同分辨率的網(wǎng)格進行網(wǎng)格收斂性研究.相關(guān)網(wǎng)格信息由表3 給出.圖7 給出的表面熱流系數(shù)分布情況可以看出,細網(wǎng)格2 和網(wǎng)格3 計算結(jié)果較為吻合,粗網(wǎng)格結(jié)果表現(xiàn)出一定差異.為最大化提升計算效率并降低收斂因素,本算例采用網(wǎng)格2 作為模擬計算模型.

圖7 不同網(wǎng)格分辨率下Z=0 m 處壁面熱流系數(shù)分布Fig.7 Distribution of heat flux coefficient at Z=0 m under different grid resolutions

表3 網(wǎng)格無關(guān)性對比Table 3 Mesh independent contrast

壁面熱流系數(shù)由下式計算所得,其中q為表面熱流密度

圖8 對比了前緣半徑R=0.5 mm 模型對稱面中心線上NCCR 模型和實驗結(jié)果的壁面熱流系數(shù),在連續(xù)流區(qū)域,根據(jù)公式計算駐點熱流系數(shù)實驗值為0.026 98,NS 模型計算的駐點熱流系數(shù)為0.026 4,偏差為 2.15%,NCCR 模型計算的駐點熱流系數(shù)為0.026 49,偏差為1.81%,Fay-Riddell 公式計算的駐點熱流系數(shù)值約為 0.028 14,偏差為 4.3%.3 種計算方法同實驗值偏差均保持在5%以內(nèi),證明在33 km高度,Ma=8.0 來流條件下飛行器尖化前緣區(qū)域并未出現(xiàn)局部稀薄氣體效應(yīng),且在連續(xù)流域NCCR 模型計算結(jié)果同實驗值基本一致.

圖8 Z=0 處模型的迎風(fēng)面與背風(fēng)面物面熱流系數(shù)分布(Ma=8)Fig.8 The distribution of surface heat flux coefficient at Z=0 m(Ma=8)in windward side and leeward side

因模型下表面為迎風(fēng)面,上表面為背風(fēng)面,故尖化前緣構(gòu)型的下表面壁面熱流系數(shù)大于上表面的壁面熱流系數(shù).從圖8 中可以看出背風(fēng)面的壁面熱流實驗值與數(shù)值模擬值保持一致.迎風(fēng)面的前半段部分,實驗值與NCCR 模型計算值吻合,但在迎風(fēng)面尾部存在偏差,偏差產(chǎn)生原因可能是迎風(fēng)面尾部壁面發(fā)生轉(zhuǎn)捩出現(xiàn)湍流,引起壁面與氣體分子間的額外能量輸運,導(dǎo)致壁面熱流系數(shù)偏高,而本文NCCR 模型在計算時均使用層流模型,因此可能在湍流區(qū)NCCR 模型與實驗值之間存在偏差.

圖9 是前緣半徑R=0.5 mm 模型Ma=8.0 條件下,單側(cè)尖化前緣熱流系數(shù)分布曲線,考慮到風(fēng)洞試驗的誤差,對比時選擇實驗值的 ± 10% 作為偏差帶的上下限.在尖化前緣區(qū)域,各點的熱流系數(shù)值大致在試驗值偏差帶內(nèi),同實驗值的吻合性較好.

圖9 尖化前緣區(qū)域壁面熱流系數(shù)分布(Ma=8)Fig.9 Heat flux coefficient distribution on sharpened leading edge(Ma=8)

圖10 是模型橫切對比圖,對比Ma=8.0 來流條件下X=0.185 m 處上表面的壁面熱流系數(shù).NS 方程、NCCR 模型計算結(jié)果同實驗值均保持一致.

圖10 X=0.185 m 處物面熱流系數(shù)分布(Ma=8)Fig.10 The distribution of surface heat flux coefficient at X=0.185 m(Ma=8)

圖11 對比距模型對稱面Z=0.04 m 處縱切面的物面熱流系數(shù)分布.模擬曲線同實驗值基本保持一致.此切面上駐點熱流系數(shù)為0.014 26,NCCR 模型計算的駐點熱流系數(shù)為0.015 47,偏差為8.48%,NS 模型計算的駐點熱流系數(shù)為0.015 48,偏差為8.55%.

圖11 Z=0.04 m 處物面熱流系數(shù)分布(Ma=8)Fig.11 The distribution of surface heat flux coefficient at Z=0.04 m(Ma=8)

3.2 60 km 處高超聲速尖化前緣局部稀薄氣體流動效應(yīng)分析

60 km 等效高度時,前緣半徑R=0.5 mm 模型尖化前緣局部努森數(shù)為0.23,屬于過渡流域,計算條件為U=3162.08 m/s,Ma=10.0,T=248.76 K,Tw=300 K,P=23.56 Pa.計算仿真模型的各項尺寸同試驗?zāi)Ｐ捅３忠恢?為了檢驗非線性本構(gòu)模型自身對稀薄非平衡現(xiàn)象的捕獲能力和準(zhǔn)確性,本文在計算過程中壁面采用無滑移等溫壁.

為了消除網(wǎng)格分布對計算結(jié)果的影響,采用3 套不同分辨率的網(wǎng)格進行網(wǎng)格收斂性研究.網(wǎng)格信息由表4 給出.圖12 給出的駐點熱流系數(shù)分布情況可以看出,網(wǎng)格2 和網(wǎng)格3 計算結(jié)果較為吻合,網(wǎng)格1 結(jié)果表現(xiàn)出一定差異.網(wǎng)格1 計算駐點壁面熱流系數(shù)為0.630 56,網(wǎng)格2 計算駐點壁面熱流系數(shù)為0.633 06,網(wǎng)格3 計算值為0.634 55,為最大化提升計算效率并降低收斂因素,本算例采用網(wǎng)格2 作為模擬計算模型.

圖12 不同網(wǎng)格分辨率下Z=0 m 壁面熱流系數(shù)分布Fig.12 Distribution of heat flux coefficient at Z=0 m under different grid resolutions

表4 網(wǎng)格無關(guān)性對比Table 4 Mesh independent contrast

圖13 對比了NCCR 模型和NS 方程計算的沿駐點線流場溫度分布曲線,NCCR 模型與NS 方程計算的激波位置、激波層內(nèi)溫度分布和溫度峰值存在差異.在60 km 高度時,尖化前緣駐點存在局部稀薄氣體效應(yīng),NCCR 模型不僅能修正存在局部稀薄效應(yīng)的駐點熱流密度,而且計算沿駐點線溫度分布結(jié)果更接近真實物理解.

圖13 兩種方法計算沿駐點線溫度分布(Ma=10)Fig.13 Temperature distributions along the stagnation line (Ma=10)

圖14 對比了NCCR 模型與采用無滑移邊界條件的NS 方程分別計算前緣半徑R=0.5 mm 飛行器的駐點前溫度云圖,因尖化前緣處于過渡流區(qū)域,NS 方程計算結(jié)果在此區(qū)域內(nèi)已不再準(zhǔn)確,對比發(fā)現(xiàn)NCCR 模型計算的駐點前區(qū)域,激波、激波脫體區(qū)域、邊界層三者之間的界限逐漸模糊和重合,NCCR模型計算的激波層內(nèi)溫度峰值要比NS 方程的峰值要低,激波層內(nèi)部溫度變化比NS 方程更加均勻.許多已發(fā)表的研究指出NCCR 模型能夠準(zhǔn)確描述過渡流區(qū)的高超聲速流動[16,33],證明NCCR 模型在連續(xù)/稀薄耦合流動求解中的潛力和工程價值.

圖14 駐點附近的溫度對比(Ma=10)Fig.14 Comparison of the temperature near the stagnation point(Ma=10)

圖15 對比了前緣半徑R=0.5 mm 模型的NCCR模型和實驗結(jié)果的壁面熱流系數(shù),高超聲速飛行器的熱防護系統(tǒng)設(shè)計初期會將局部高熱流密度作為設(shè)計重點,而通常情況下,駐點熱流密度屬于模型的熱流密度峰值,因此我們重點關(guān)注駐點熱流密度的偏差情況.在過渡流區(qū)域,駐點實驗熱流系數(shù)為0.474 88,NS 方程計算駐點熱流系數(shù)為0.633 06,與實驗值偏差為33.31%;NCCR 模型計算的駐點熱流系數(shù)為0.530 76,偏差為11.77%;而Fay-Riddell 公式計算的駐點熱流系數(shù)值約為0.615 11,偏差為29.5%.在非駐點區(qū)域局部稀薄效應(yīng)影響較小,NCCR 模型計算結(jié)果同實驗值基本保持一致.在圖14 中激波層內(nèi)溫度分布及溫度峰值的計算上NCCR 模型的結(jié)果更準(zhǔn)確和可靠.

圖15 Z=0 m 處迎風(fēng)面與背風(fēng)面物面熱流系數(shù)分布(Ma=10)Fig.15 The distribution of surface heat flux coefficient at Z=0 m(Ma=10)in windward side and leeward side

通過圖15 壁面熱流系數(shù)的分布可以看出: 除駐點附近熱流系數(shù)外,NCCR 模型與NS 方程計算迎風(fēng)面與背風(fēng)面熱流系數(shù)與實驗值基本一致.證明壁面除前緣外的部分局部稀薄氣體效應(yīng)較弱,NS 方程的計算結(jié)果能夠滿足精度要求.

圖16 對比了距模型對稱面Z=0.02 m 處尖化前緣模型縱切面的物面熱流系數(shù)分布.存在局部稀薄效應(yīng)的尖化前緣區(qū)域中,NS 方程計算駐點熱流系數(shù)為0.491 43,NCCR 模型計算駐點熱流系數(shù)為0.414 54,NS 方程計算結(jié)果同NCCR 模型偏差為18.55%.而非前緣區(qū)域流動處于連續(xù)流狀態(tài),NS 方程同NCCR模型結(jié)果與實驗值基本保持一致.

圖16 Z=0.02 m 處物面熱流系數(shù)分布(Ma=10)Fig.16 The distribution of surface heat flux coefficient at Z=0.02 m(Ma=10)

4 結(jié)論

本文針對尖化前緣構(gòu)型在高超聲速連續(xù)/稀薄流中的氣動加熱問題開展數(shù)值計算研究和風(fēng)洞實驗驗證.根據(jù)實驗來流條件,借助NCCR 模型對尖化前緣構(gòu)型在Ma=8,10 的0°攻角狀態(tài)下的物面熱流系數(shù)進行了數(shù)值計算,并與風(fēng)洞試驗結(jié)果進行對比驗證.通過對比物面熱流系數(shù)和沿駐點線溫度分布,檢驗了NCCR 模型相比于NS 方程在計算三維高速局部稀薄流動中的改善水平.本文研究結(jié)論總結(jié)如下.

(1)飛行器尖化前緣處于低空連續(xù)流條件下,駐點區(qū)域局部稀薄氣體效應(yīng)并不明顯,NCCR 模型與NS 方程計算結(jié)果同實驗值保持一致,偏差均保持在5%以內(nèi).

(2)飛行器尖化前緣處于中高空來流條件時,尖化前緣的局部稀薄氣體效應(yīng)隨高度增加逐漸顯著,NCCR 模型不僅能夠修正局部稀薄流動區(qū)域的壁面熱流系數(shù)值,而且近壁面流場內(nèi)參數(shù)計算的結(jié)果更接近真實物理解.

(3)在等效高度60 km,來流馬赫數(shù)Ma=10.0的條件下,NS 方程計算駐點熱流系數(shù)偏差為33.31%,NCCR 模型計算駐點熱流系數(shù)偏差為11.77%;在相同計算條件下,不考慮邊界條件帶來的影響,NCCR模型計算駐點熱流系數(shù)在誤差允許范圍內(nèi)同試驗值更為接近,體現(xiàn)出非線性本構(gòu)方程求解稀薄區(qū)域流動的優(yōu)勢.