基于非線性耦合本構關系模型的尖化前緣氣動加熱影響研究1)

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

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

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

引言

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

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

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

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

1 非線性耦合本構關系模型

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

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

經過Eu 的高階矩封閉以及Eu 和Myong 絕熱假設簡化[8],可以得到非線性耦合本構關系模型表達式為

為方便比較,NS 方程線性本構模型也在這給出

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

2 算例驗證

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

表1 流動狀態Table 1 Flow state

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

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

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

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

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

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

3 尖化前緣氣動加熱計算

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

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

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

風洞試驗在中國科學院力學研究所高溫氣體動力學國家重點實驗室JF8A 和JF10 高超聲速風洞中進行.試驗模型的來流條件見表2,Ma=8.0 的飛行條件在JF8A 中試驗[31],JF10 風洞中進行Ma=10.0 飛行條件試驗[32],后續計算模型的流動參數同風洞試驗參數保持一致.

表2 流場狀態參數Table 2 Flow field status parameters

3.1 33 km 處高超聲速尖化前緣連續流動效應分析

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

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

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

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

表3 網格無關性對比Table 3 Mesh independent contrast

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

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

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

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

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

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

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

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

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

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

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

60 km 等效高度時,前緣半徑R=0.5 mm 模型尖化前緣局部努森數為0.23,屬于過渡流域,計算條件為U=3162.08 m/s,Ma=10.0,T=248.76 K,Tw=300 K,P=23.56 Pa.計算仿真模型的各項尺寸同試驗模型保持一致.為了檢驗非線性本構模型自身對稀薄非平衡現象的捕獲能力和準確性,本文在計算過程中壁面采用無滑移等溫壁.

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

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

表4 網格無關性對比Table 4 Mesh independent contrast

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

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

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

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

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

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

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

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

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

4 結論

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

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

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

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