深空再入飛行器燒蝕粗糙表面高超聲速轉捩預測

李 齊, 趙 瑞, 陳 智, 郭 斌, 王 強

(1. 北京空間飛行器總體設計部, 北京 100094; 2. 北京理工大學, 北京 100081;3. 中國航天空氣動力技術研究院, 北京100074)

引 言

人類自開啟航天事業以來, 探索地外天體的目標從未停止, 并且越來越深遠。20世紀50～60年代, 美國與蘇聯分別提出了各自的月球探測計劃, 并先后實現了月球采樣返回和載人登月返回[1]。20世紀90年代, 美國和日本等國家又先后啟動了彗塵[2]、 太陽風[3]、 小行星[4]等探測取樣返回計劃, 并于2004—2010年期間先后實現了樣品取樣返回。2014年11月, 我國首次實施月地高速再入返回獲得圓滿成功, 掌握了第二宇宙速度再入返回技術[5], 并以此為基礎于2020年12月成功實現了月球取樣返回目標[6]。2017年12月, 時任美國總統特朗普簽署白宮1號太空政策指令, 重啟登月計劃, 目標在2025年前后重返月球, 實現新世紀載人月球探測[7]。2022年1月28日由中華人民共和國國務院新聞辦公室發布的《2021中國的航天》白皮書中提出, 未來五年要“發射小行星探測器、 完成近地小行星采樣和主帶彗星探測”, 并“深化載人登月方案論證, 組織開展關鍵技術攻關”[8]。由此可見, 月球與深空地外天體探測(含載人探測)是未來世界各國空間探測與科學研究的熱點項目。而開展月球與深空地外天體探測, 必然少不了樣品或載人返回再入飛行器。

由于深空再入飛行器再入速度大, 為提高氣動減速效率, 一般采用大鈍度球冠或球錐迎風外形[9]。而大鈍頭迎風外形與超高速來流相互作用, 加之高溫真實氣體效應, 導致鈍頭表面熱流高、 加熱量大。為有效降低氣動加熱對主體結構、 艙體設備和航天員安全性帶來的影響, 防熱結構須采用碳基燒蝕型材料。而碳基防熱材料密度低、 易燒蝕, 在高焓高熱流長時間氣動加熱下會造成防熱結構表面粗糙度急劇增加, 形成分布式粗糙元表面。扁平的迎風前體與粗糙元結構結合, 極易造成深空再入飛行器迎風面在高超聲速下出現流動失穩, 導致流動轉捩甚至演化為湍流, 使表面熱流分布發生巨大變化, 給飛行器安全帶來極大挑戰。其中的典型案例包括: Genesis返回艙防熱罩對接孔凹坑后緣流動轉捩造成的局部過熱問題[10], 以及NASA對Artemis 1號飛行試驗返回后的獵戶座飛船開展技術分析時發現了“大底隔熱板燒蝕量大于地面設計預估值”的問題[11]。為適應深空探測復雜系統與運載能力之間的匹配[12], 再入飛行器質量約束是系統設計必須滿足的關鍵條件, 因此防熱結構設計必須做到節約而高效。由此可見, 深空再入飛行器高超聲速邊界層失穩機制分析和轉捩位置的準確預測是影響深空再入返回飛行安全性的關鍵難題。

高超聲速邊界層轉捩研究一般主要面向尖前緣的高超聲速飛行器外形, 國內外學者曾在失穩特征、 轉捩機理、 感受性特征以及轉捩預測方法等方面取得了一些研究成果[13]。對于再入飛行器大鈍頭前體造成的流動轉捩及其熱影響問題, 早期在Apollo[14]、 Galileo和Pioneer Venus等[15]進/再入器上都發現過轉捩現象, 但由于條件限制等問題, 未見深入的分析研究。Edquist等[16]首次針對火星探測器MSL(Mars smart lander)的大鈍頭前體外形利用地面試驗確定的動量Reynolds數半經驗工程方法來預測流動轉捩的發生位置, 但后來與飛行試驗結果對比發現有相當大的出入[17]。此外, Horvath等[18]分別利用細長錐外形、 MSL和Genesis號外形[19], 通過開展高超聲速風洞試驗和PSE方法的穩定性分析, 研究了對應外形下不同形式的表面粗糙度結構對鈍頭再入器迎風錐表面轉捩模式的作用, 并綜合了磷光測熱結果和數值模擬結果, 確定了轉捩發展造成的氣動熱分布變化, 提出了基于工程擬合的粗糙度轉捩預測模型。

綜上所述, 目前大部分對于高超聲速邊界層轉捩的計算研究都是單獨針對宏觀外形或粗糙表面微觀外形而開展的, 且大部分是基于地面試驗數據進行轉捩準則的建立與修正。對于大鈍頭宏觀外形與分布式粗糙表面微觀外形結合下的高超聲速轉捩問題, 轉捩機制與基于轉捩模式的數值預測研究較為少見。本文擬用適用于深空再入返回的大鈍頭迎風大底外形與沙粒式分布粗糙燒蝕表面為研究對象, 分別基于高超聲速與粗糙元修正的γ-Reθ轉捩模式和k-ω-γ轉捩模式, 開展高超聲速邊界層轉捩位置預測、 機制分析和參數影響規律研究。

1 數值方法

1.1 控制方程與數值格式

考慮量熱完全氣體, 對三維非定?？蓧嚎sNavier-Stokes(N-S)方程進行Favre平均, 得到可壓縮湍流的Reynolds平均N-S方程

式中,ui,p,E,h分別為速度分量、 壓力、 總能、 焓。且

其中,T,μl,κl分別為溫度、 分子黏性系數、 分子熱交換系數。

考慮高溫真實氣體效應的流動控制方程為積分形式的多組分化學非平衡N-S方程[20], 忽略輻射以及徹體力的影響, 方程形式如下

本文針對高溫真實氣體效應, 采用了基于5組分的Dunn-Kang模型[21], 考慮了完全催化壁模型[22], 采用了AUSM+up格式[23]進行數值解算。

1.2 湍流模型

采用基于渦黏性假設的兩方程剪切應力輸運(shear-stress-transport, SST) 湍流模型, 在邊界層的黏性底層和對數律層采用k-ω模型, 在邊界層的虧損律層采用k-ω模型和k-ε模型的混合, 在自由剪切層中采用k-ε模型[24]?？刂品匠虨?/p>

其中

湍流運動黏性系數μt和動力黏性系數νt為

模型常數由兩部分混合求得

φ=F1φ1+(1-F1)φ2

混合函數F1及其他函數定義如下

1.3 轉捩模式

1.3.1γ-Reθ轉捩模式

γ-Reθ轉捩模式通過聯立求解間歇因子與轉捩動量厚度Reynolds數的輸運方程, 根據局部渦量Reynolds數和臨界動量厚度Reynolds數的比值判斷轉捩[25]。γ-Reθ關聯模型的間歇因子γ與動量厚度Reynolds數Reθ的方程分別為

Pγ=Flengthca1ρS(γFonset)0.5(1-ce1γ)
Dγ=ca2ρΩγFturb(ce2γ-1)

為適應高超聲速轉捩流動

Flength=20,ce2=20

為考慮表面粗糙度對邊界層轉捩的影響, 對轉捩模型中的經驗關系式進行修正[26]。引入等效沙粒高度ks和當地邊界層位移厚度對轉捩動量厚度Reynolds數進行修正, 建立粗糙表面條件下轉捩動量厚度Reynolds數Reθt_rough, 表達式如下

其中

fTu=max[0.9, 1.61-1.15·e-Tu]

fΛ用于描述粗糙度的形狀、 排列規律等幾何構型, 本文取1。

此外, 為考慮表面粗糙度對流動轉捩后湍流邊界層的影響, 對SST湍流模型中比耗散率ω進行表面粗糙度修正[27], 具體形式如下

其中,SR定義如下

1.3.2k-ω-γ轉捩模式

k-ω-γ轉捩模式以SST湍流模型為基礎, 由關于湍動能k、 比耗散率ω以及間歇因子γ的3個輸運方程構成, 可以在有效黏性系數中考慮非湍流脈動影響, 并借鑒Langtry等[25]和Langel等[26]構造的轉捩模型基于當地變量的優點, 構造了間歇因子γ輸運方程耦合層、 湍流計算。其總體框架為

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

式中,Fonset為轉捩起始位置函數,d為物面距離,Eu為當地流體相對壁面的平均流動動能,ν為分子運動黏性系數。有

其中,ζeff為有效長度尺度防熱結構燒蝕導致的表面粗糙, 可簡化為等效沙粒分布式粗糙度, 在k-ω-γ轉捩模式中引入粗糙度放大因子Ar輸運方程來描述壁面粗糙度對邊界層轉捩的影響機理和作用效果[28], 具體構造如下

Ar的壁面邊界條件以Sigmoid函數給定

下式中,k+是無量綱的等效沙粒粗糙度高度, 由壁面摩擦速度uτ和等效粗糙度高度ks共同決定

1.4 方法驗證

為驗證本文轉捩模式能否正確預測壁面粗糙度對邊界層轉捩的影響, 采用美國NASA蘭利實驗室[29]在20 in(1 in=25.4 mm),Ma=6風洞中采用的帶有分布式沙粒粗糙度的半球頭模型進行驗證, 如圖1所示, 半球模型的直徑為152.4 mm。

(a) Oblique view (b) Side view

(c) Front view (b) Close-up view圖1 NASA半球頭模型Fig. 1 NASA hemisphere model photographs

來流條件為: Mach數6.04, 攻角0°, 壁面溫度300 K, 單位Reynolds數2.18×107/m。選取80-mesh粗糙元結構來考察粗糙誘導轉捩模式的預測精度, 該粗糙元均方根粗糙高度RRMS=0.03 mm。采用文獻中90%粗糙度包絡曲線, 并根據Dassler等[30]的經驗公式, 有ks=4.33RRMS, 可得等效粗糙高度ks為0.13 mm。

圖2給出了分別采用兩種轉捩模式考慮與不考慮粗糙度放大因子計算得到的傳熱系數h/href分布與地面測試結果的對比。其中, 后綴為orig表示不考慮粗糙放大因子的轉捩模式計算結果, 后綴為rough表示考慮粗糙放大因子的轉捩模式計算結果, exp為地面測試結果。由圖可知, 在給定來流條件下, 不考慮粗糙元誘導時半球頭表面流動不發生轉捩, 熱流由頭部向肩部逐漸減小。而考慮粗糙元誘導后, 兩種轉捩模式計算結果均顯示s/R=0.2～0.3位置邊界層流動由層流轉捩為湍流, 熱流顯著增大。由于h/href是判別流動轉捩的重要宏觀物理量, 由圖中對比可見, 本文采用的兩種粗糙元誘導轉捩模式數值結果與實驗測得的轉捩起始位置和轉捩后最高熱流吻合良好。其中,k-ω-γ模式對應轉捩位置與實驗值吻合度更高, 而γ-Reθ模式對應轉捩后熱流與實驗值更為接近, 這與γ-Reθ模式的轉捩區模型參數由地面實驗修正而來有關。

圖2 不同轉捩模式下半球頭80-mesh粗糙模型熱流密度計算值與實驗值對比Fig. 2 Comparison of calculated heat flux by different transition modes and experimental heat flux of 80-mesh rough model

2 幾何模型與計算狀態

如圖3所示, 本次研究選取的幾何模型為球錐大鈍頭迎風大底外形, 球頭半徑370 mm, 迎風半錐角63°。采用分區多塊對接結構網格, 不考慮側滑角影響, 計算網格為半模, 總網格量1×106, 沿法向進行網格加密, 首層網格高度為0.03 mm, 以保證y+≤1。

(a) Symmetry plane (b) Object plane 圖3 幾何模型與網格結構示意Fig. 3 Geometric model and grid structure

本文計算狀態如表1所示, 取典型深空再入飛行器高超聲速狀態, 高度/Mach數組合關系分別為52 km/Ma25, 49.3 km/Ma20, 45 km/Ma13, 42 km/Ma10, 根據氣動加熱狀態分別設定壁面溫度為3 400, 3 000, 2 500, 2 100 K, 考慮了1.5, 3 mm兩種等效粗糙度高度。此外, 由于返回艙采用自旋彈道式再入飛行, 從宏觀時間來看表面相同軸向位置的氣動加熱相等, 另外考慮大鈍頭迎風面熱流均勻化與三維燒蝕傳熱的拉平效應等, 本文設置整個迎風大底表面采用相同的等效粗糙高度。

表1 再入飛行器來流條件與壁面條件

3 計算結果分析

3.1 粗糙高度對轉捩模擬影響分析

采用k-ω-γ轉捩模式對1.5, 3 mm等效粗糙高度的幾何模型迎風大底表面間歇因子γ進行模擬分析, 圖4和圖5分別顯示了0°攻角下狀態2和狀態3對應等效粗糙高度ks=1.5, 3 mm大底表面間歇因子γ的分布。k-ω-γ模式計算將間歇因子γ開始顯著增長的位置定義為轉捩起始位置, 將γ在(0.1～1)的范圍定義為轉捩區,γ>1為湍流區。由兩圖對比可見, 同樣的來流狀態下, 隨著等效粗糙高度的增加, 壁面間歇因子增長起始點逐漸由肩部向中心推進。當等效粗糙高度ks=1.5 mm時, 兩個狀態的大底表面流動均為層流; 當粗糙高度增長為3 mm時, 兩狀態的大底表面在球錐面交界處開始發生轉捩, 錐面大面積流動已發展為湍流。粗糙元高度通過增加當地邊界層厚度而提高了當地Reynolds數, 從而促使大底表面轉捩提前發生。

圖4 壁面間歇因子分布云圖(ks=1.5 mm, α=0°)Fig. 4 γ distribution on the wall(ks=1.5 mm, α=0°)

圖5 壁面間歇因子分布云圖(ks=3 mm, α=0°)Fig. 5 γ distribution on the wall(ks=3 mm, α=0°)

3.2 來流Reynolds數對轉捩模擬影響分析

采用k-ω-γ轉捩模式對不同來流Reynolds數條件下3 mm等效粗糙高度的幾何模型迎風大底表面間歇因子γ進行模擬分析, 圖6顯示了0°攻角狀態的間歇因子γ分布云圖。整體來看, 在粗糙高度3 mm表面形貌下, 當前計算的所有狀態下大底表面均存在轉捩與湍流區。隨著飛行高度的減小, 來流單位Reynolds數逐漸增大, 轉捩起始位置由肩部逐漸向頭部中心即上游移動; 轉捩區寬度隨著來流Reynolds數的增加而逐漸收縮, 而湍流區則逐漸增大。

與粗糙高度影響不同的是, 來流Reynolds數不是通過增加當地邊界層厚度而是通過直接提高當地Reynolds數水平誘導轉捩提前。與粗糙高度相比, 來流Reynolds數對轉捩起始位置與轉捩區的影響線性度更強。

(a) Re/D=4.0×105/m (b) Re/D=5.0×105/m

(c) Re/D=5.8×105/m (d) Re/D=6.0×105/m圖6 不同來流Reynolds數下壁面間歇因子分布云圖(ks=3 mm, α=0°)Fig. 6 γ distribution on the wall at different inflow Reynolds numbers(ks=3 mm, α=0°)

3.3 攻角對轉捩影響分析

圖7為狀態4、 10°攻角、 等效粗糙高度ks=3 mm條件下返回艙壁面間歇因子分布云圖。與圖6(d)對比可知, 由于攻角的存在, 轉捩起始位置提前, 駐點向迎風面移動, 轉捩區不再關于y=0對稱, 而是大部分位于迎風面。圖8為狀態4、 10°攻角、 等效粗糙高度ks=3 mm條件下層流與轉捩模式計算所得大底子午線壁面熱流分布曲線的對比。與圖7相對應, 由于攻角的存在, 子午線壁面熱流分布曲線結果沒有關于y=0對稱, 轉捩模式背風面熱流值大于迎風面, 轉捩后背風肩部熱流甚至與迎風肩部相當, 可達當地層流熱流的2倍以上。由上述結果分析可知, 由于攻角的存在, 大底背風面轉捩位置提前, 湍流區擴大, 熱流增長效應顯著增加。

在一般認識下, 當有攻角存在時, 大底迎風面密度應大于背風面密度(圖9(a)), 因而迎風面動量厚度Reynolds數應大于背風面。但據上述計算結果可知, 大底背風面反而轉捩提前。從大底迎背風面動量厚度對比(圖9(b))可知, 雖然大底迎風面的邊界層外緣密度是背風面的2～3倍, 但背風面的動量厚度為迎風面的3倍以上, 且動量厚度增加導致背風面邊界層外緣速度也遠大于迎風面, 由此導致背風面的動量厚度Reynolds數是迎風面的1.5倍以上(圖9(c))。由此可見, 有攻角情況下, 大底背風面動量厚度大大增加, 從而導致背風面動量厚度Reynolds數增大, 因此轉捩位置提前, 湍流區擴大。

圖7 壁面間歇因子分布云圖(42 km/Ma10, α=10°, ks=3 mm)Fig. 7 γ distribution on the wall(42 km/Ma10, α=10°, ks=3 mm)

圖8 大底子午線壁面熱流分布曲線(42 km/Ma10, α=10°, ks=3 mm)Fig. 8 Wall heat flow distribution along the meridian of heat shield(42 km/Ma10, α=10°, ks=3 mm)

(a) Density distribution

(b) Momentum thickness distribution

(c) Momentum thickness Reynolds number distribution圖9 有攻角下大底表面特征參數分布(42 km/Ma10, α=10°)Fig. 9 Feature parameter distribution of the heat shield with attack angle(42 km/Ma10, α=10°)

3.4 化學非平衡對轉捩模擬影響分析

采用γ-Reθ轉捩模式, 分別基于完全氣體和化學非平衡兩種氣體模型的基本流, 對不同高度狀態下3 mm等效粗糙高度的幾何模型迎風大底表面的熱流密度分布進行模擬計算, 圖10給出了狀態1和狀態2不同氣體模型對應的層流/轉捩模式下粗糙大底表面熱流密度分布曲線的對比。由圖可見, 在所計算狀態下, 完全氣體基本流計算所得熱流在球錐交界位置即開始高于層流熱流, 且隨著高度的降低和來流Reynolds數的增加, 完全氣體模型在轉捩與湍流區的熱流可達1.4倍以上的當地層流熱流; 而化學非平衡基本流所得熱流在不同高度下與層流熱流則沒有明顯變化。由此可見, 化學非平衡基本流可有效抑制高超聲速邊界層轉捩的發展及對熱流的影響。

(a) 52 km/Ma25 perfect gas

(b) 52 km/Ma25 chemical non-equilibrium

(c) 49.3 km/Ma20 perfect gas

(d) 49.3 km/Ma20 chemical non-equilibrium圖10 不同氣體模型下大底子午線壁面熱流分布曲線(α=0°, ks=3 mm)Fig. 10 Wall heat flow distribution along the meridian of heat shield with different gas models (α=10°, ks=3 mm)

4 結論

本文以典型深空再入飛行器迎風大底外形為研究對象, 以分布式粗糙元結構為表面特征, 采用基于粗糙放大因子修正的γ-Reθ與k-ω-γ轉捩模式, 對深空再入高超聲速典型狀態開展了邊界層轉捩預測模擬研究, 分析了粗糙高度、 來流Reynolds數、 攻角以及化學非平衡基本流對深空再入飛行器高超聲速邊界層轉捩位置與轉捩熱效應的影響規律。主要結論如下:

1) 對于本次研究的小尺寸大鈍頭迎風前體外形, 分布式粗糙元與來流Reynolds數的增加都會通過增大當地Reynolds數從而誘導轉捩, 使轉捩起始位置逐漸向上游發展。其中, 粗糙元誘導轉捩效應更為明顯, 非線性度更強。

2) 攻角可使得背風面動量厚度Reynolds數大大增加, 從而導致轉捩位置提前, 湍流區擴大, 背風面熱流顯著增長。

3) 化學非平衡基本流可有效抑制高超聲速邊界層轉捩的發展。

后續, 對于深空再入飛行器燒蝕粗糙表面高超聲速邊界層轉捩預測與影響分析的研究, 可圍繞粗糙度對流動穩定性的影響機制分析、 建立基于粗糙元感受的流動穩定性分析模型、 建立多種分布式粗糙元等效粗糙因子數學模型等方面來開展。