火星再入飛行器風洞試驗與真實飛行 之間相關性的探討

劉方彬， 袁軍婭

(1. 北京航空航天大學宇航學院, 北京 100083； 2. 深圳易信科技股份有限公司, 深圳 518000)

火星大氣的主要成分是CO2，大氣非常稀薄，密度僅為地球大氣的1%。相對于地球大氣，再入飛行器進入火星大氣時，雖然同樣具有較高的飛行速度，但是由于火星大氣密度較低，在高速飛行時激波后流動的熱力學非平衡特性更為突出。作為火星大氣主要成分的CO2氣體，具有和N2、O2等氣體不同的熱化學特性[1]。目前對于火星再入飛行器實際飛行時，仍存在諸多認知上的不足，包括剪切層的復雜性、尺度問題和一些現象的物理模型等方面，地面試驗仍是火星再入飛行器研制中必不可少的驗證手段[2]。

由于風洞條件的限制，難以完全模擬飛行環境，且在風洞試驗中有尺寸效應、邊界效應或邊界干擾、支架干擾等干擾，這就存在一個與真實飛行條件下的相似理論及相關性的問題。目前關于地球再入飛行器的風洞試驗相關性研究較多，相似準則的選擇也有多種，而關于火星探測器風洞試驗相關性研究的國內外公開文獻極少，且成熟度不高。例如苗文博等[1]則利用數值方法研究了火星再入飛行器表面熱環境，通過研究認為，在駐點附近區域近似為熱力學平衡，隨著高度逐漸增加，則熱力學非平衡效應更明顯。董維中[3]在博士論文中根據鈍錐標模ELECTRE模型風洞試驗和飛行試驗數據之間對比，認為保持總焓值不變，利用雙尺度參數可以對模型頭部進行氣動熱/力的外推。然而ELECTRE模型的試驗介質均為空氣，沒有進行介質為CO2的試驗。Bur等[4]通過研究風洞與“Pathfinder”(“探路者號”)的計算數據，認為風洞條件下駐點和前體部分熱流計算峰值要高28%。Armenise等[5]則對風洞條件和飛行條件下的計算化學反應模型進行了分析和驗證。Paterna等[6]通過試驗和計算，認為火星再入飛行器在CO2氣體環境下的計算，采用催化壁條件獲得到的物理參數更為合適。

鑒于風洞試驗在火星再入飛行試驗中的重要性，需要建立風洞條件與飛行條件之間的相關性。基于眾多科研工作者的開拓性研究，本文利用了“Pathfinder”風洞試驗數據[7-8]和飛行試驗數據[9-10]，通過數值方法，建立了火星再入飛行器氣動力和氣動熱的風洞參數外推方法。

1 物理模型和計算方法

1.1 物理模型

本次計算采用的是“Pathfinder”火星再入飛行器的飛行數據[11]和風洞試驗數據[12]，其中風洞數據來源于“MP-1”模型的風洞試驗數據[12]，飛行數據來源于“Pathfinder”的飛行數據[10]。兩者前體幾何相似，縮尺比例為1|∶ 52.15，如圖1[11]和圖2[12]所示，其中圖2中1 in= 0.025 4 m。表1[12]中CASE 1和CASE 2是兩種常規高超聲速空氣風洞試驗條件及其對應的飛行狀態。CASE 3是試驗介質分別為空氣、CO2的高焓風洞試驗數據及其對應的飛行數據。在本文計算中，認為火星大氣其主要成分是由體積分數97%的CO2和3%的O2組成；空氣的主要成分由體積分數79%的N2和21%的O2組成；來流風洞的主要成分由體積分數為100%的CO2組成。

CASE 1中，風洞試驗焓值達到H0,2-H298=(0.756±0.037 8) MJ/kg。CASE 2 中，風洞試驗焓值為H0,2-H298=(0.764±0.007 64) MJ/kg。CASE 3中，介質為CO2的風洞試驗焓值為H0,2-H298=(12.25±0.26) MJ/kg，介質為空氣時的風洞試驗焓值為H0,2-H298=(14.18±0.28) MJ/kg。 CASE 1、CASE 2和CASE 3中，風洞模型的雙尺度ρ∞L與飛行器的雙尺度一致。雙尺度ρ∞L中的ρ∞為來流密度，L為飛行器或風洞模型頭部的有效直徑。表1中u∞為來流速度，T∞為來流溫度，p∞為來流壓力，Ma∞為來流馬赫數，Re∞為來流雷諾數。“97%CO2+3%N2”代表來流介質的主要組成成分是體積分數97%的CO2和3%的N2，即真實飛行條件中的火星大氣，“21%O2+79%N2”代表來流介質的主要組成成分是體積分數21%的O2和79%的N2，即空氣風洞試驗中的空氣，“100%CO2”代表來流介質的主要組成成分是體積分數100%的CO2，即CO2風洞試驗中CO2。CASE 1和CASE 2中的風洞試驗條件為低焓風洞試驗條件，CASE 3 中的風洞試驗條件為高焓風洞試驗條件。

圖1 “Pathfinder”計算模型[11]Fig.1 Calculation model of “Pathfinder”[11]

圖2 “MP-1”計算模型[12]Fig.2 Calculation model of “MP-1”[12]

表1 計算條件[12]Table 1 Calculation conditions[12]

1.2 計算方法

求解的控制方程為采用化學反應全Navier-Stokes方程求解二維軸對稱外形，并基于以下假設：流動為熱力學非平衡；忽略輻射和徹體力的影響；流動質量采用雙組元氣體模型；溫度模型采用Park雙溫模型[13]；壁面采用完全催化壁，壁面單原子和離子濃度為0；常溫風洞條件下壁面溫度設置為300 K，高焓風洞條件和飛行條件下的壁溫設為2 000 K。

對于空氣的化學反應模型，采用Park 5組分化學反應模型[13]，如表2所示。對于CO2的化學反應模型，則采用8組分9化學反應模型[14]，如表3所示。

采用有限體積法對控制方程進行數值求解，空間格式采用二階TVD格式，黏性項采用二階中心差分，時間格式采用隱式求解。

計算網格采用代數關系方法生成，為了模擬附面層，壁面網格采用了線性函數來加密。由于熱流是以黏性為主導的物理量，高超聲速飛行器的氣動熱計算對網格依賴性較高。根據文獻[15]，在高超聲速中，網格雷諾數小于8可以獲得收斂的熱流結果[15]。所以在本文計算中，所有網格雷諾數均小于8。“MP-1”與“Pathfinder”網絡如圖3、圖4所示。

表2 5組分化學反應模型[13]Table 2 Mechanism with five species chemical reactions[13]

表3 8組分9化學反應模型[14]Table 3 Mechanism with eight species and nine chemical reactions[14]

圖3 “MP-1”的網格Fig.3 Mesh of “MP-1”

圖4 “Pathfinder”的網格Fig.4 Mesh of “Pathfinder”

2 計算結果

2.1 對比驗證

通過與NASA Langley研究中心“Pathfinder”飛行狀態計算結果[9-12]的對比，來驗證本文所采用的熱力學模型和計算方法能否用于飛行條件下的飛行器流場計算；采用表1中CASE 3的CO2條件開展計算，將計算結果與風洞試驗結果[12]進行對比，來驗證本文所采用的熱力學模型和計算方法能否運用于風洞條件下的模型流場計算。

圖5為采用本文方法計算得到的風洞條件下的模型表面熱流與風洞試驗數據[12]的對比曲線。

圖5 氣動熱對比(風洞條件)Fig.5 Comparison of aerodynamics (wind tunnel condition)

圖中“4 772 m-s-wall=2 000-FCW-CO2”表示來流速度為4 772 m/s、來流介質為CO2，對應表1中的CASE 3下的風洞條件中的CO2風洞條件；壁溫Tw=2 000 K、壁面條件為完全催化壁(FCW)。x/S表示壁面投影到x軸的坐標，x軸為飛行器或風洞模型的幾何旋轉軸，坐標原點為駐點位置，下同；Qdot為壁面熱流；“Experiment”表示文獻[12]中的風洞條件下的試驗數據(下文圖中的曲線標識類似)。根據參考文獻[6]，在飛行條件下的飛行器流場計算中，壁面條件采用完全催化壁。從圖5可以看出，風洞條件下的熱流峰值計算q0,CFD與試驗q0,Experiment相差0.1%，試驗點上最大的誤差為8%，證明所采用的熱力學模型和計算方法在風洞條件下的計算精度可以滿足分析要求。

表4為采用本文的計算方法得到的飛行狀態駐點熱流與文獻[8]中計算結果的對比。從表4可以看出，本文的計算方法在飛行狀態時的計算結果與NASA Langley中心計算結果最大誤差為6%，表明本文的計算方法對飛行條件下的計算精度滿足分析要求。

表4 對比條件和計算結果(飛行條件)Table 4 Comparison conditions and computation results(flight condition)

2.2 表面壓力試驗相關性

由于低焓風洞中流場溫度較低，所以在“MP-1”計算中僅考慮冷壁(壁溫Tw=300 K)非催化壁條件；在高焓風洞條件下，來流溫度較高，為了比較冷壁和熱壁(Tw=2 000 K)條件對相關性的影響，在“MP-1”計算中分別采用冷壁和熱壁2 種計算條件，以便對比分析火星再入飛行器飛行狀態下冷壁和熱壁對相關性的影響。

圖6 風洞條件和飛行條件壓力系數Fig.6 Pressure coefficient of wind tunnel and flight conditions

圖7 風洞條件和飛行條件無量綱壓力Fig.7 Dimensionless pressure of wind tunnel and flight conditions

從圖6(a)、圖6(b)、圖7(a)和圖7(b)中可以看出，對于低焓風洞試驗的壓力系數和無量綱壓力，在駐點附近時與飛行數據相差很小；從駐點沿壁面向肩部發展，低焓風洞試驗結果與飛行數 據之間差距逐漸增大，最大相差11.5%。壁面溫度對飛行狀態下的壓力系數和無量綱壓力沒有明顯的影響。

圖8 風洞條件和飛行條件的馬赫數對比Fig.8 Mach number comparison of wind tunnel and flight conditions

圖9 風洞條件和飛行條件的溫度對比Fig.9 Temperature comparison of wind tunnel and flight

圖10 風洞條件和飛行條件的壓力對比Fig.10 Pressure comparison of wind tunnel and flight

從圖6(c)和圖7(c)可以發現，對于高焓風洞試驗壓力系數和無量綱壓力，以及高焓CO2風洞試驗壓力系數與無量綱壓力，除了模型駐點附近區域以外，其他地方與飛行狀態下的壓力系數和無量綱壓力有明顯的差距；對于高焓空氣風洞試驗的壓力系數和無量綱壓力，從駐點沿壁面向肩部發展，與飛行條件下的壓力系數和無量綱壓力幾乎吻合。高焓空氣風洞試驗無量綱壓力和壓力系數，比常溫空氣風洞試驗的無量綱壓力和壓力系數更接近于飛行條件下的無量綱壓力與壓力系數。

在低焓風洞條件與真實飛行之間出現以上情況的主要原因是馬赫數影響。由于常溫下空氣的比熱比與火星大氣比熱比接近，根據馬赫數無關原理，對于高超聲速飛行，當馬赫數較高時，壓力系數基本上與馬赫數無關，僅與流動介質的比熱比相關。根據表1可以發現，不管CASE 1和CASE 2 的風洞條件還是飛行條件，來流馬赫數都較高，因此低焓風洞試驗模型的駐點處的壓力系數與飛行條件下的駐點壓力系數接近。根據圖8(a)可以知道，從模型駐點處沿著壁面向肩部發展，激波層內的馬赫數逐漸降低，模型肩部附近馬赫數較小。根據牛頓修正公式[16]，當馬赫數逐漸減小時，壓力系數等參數逐漸減小，而且當比熱比越大，壓力系數的變化幅度也越大。所以當壓力系數和無量綱壓力逐漸從模型駐點沿壁面向肩部發展而逐漸減小，且空氣介質中的減小幅度大于CO2介質中的減小幅度。

分析高焓風洞試驗和真實飛行的流場，根據圖8(b)可知，兩者來流馬赫數都比較高，所以引起壓力系數和無量綱壓力變化差異的原因主要在于溫度。標準大氣壓下，CO2在5 000 K時完全分解，且CO2分解的溫度域很小。O2在2 000 K左右開始分解，4 000 K左右時完全分解，N2在3 000 K左右開始分解，9 000 K幾乎完全分解。從圖9和圖10可知，在高焓CO2風洞試驗中，模型駐點到肩部區域的激波層內溫度不超過4 000 K，模型頭部駐點附近激波層內的壓力略大于101 325 Pa。根據高超聲速理論，除了振動激勵現象不受壓力影響以外，如果壓力上升，其他起始溫度將會對應上升[16]。即壓力大于101 325 Pa時，氣體的分解起始溫度比101 325 Pa下的分解起始溫度要高；當壓力小于101 325 Pa時，氣體的分解起始溫度溫度比101 325 Pa下的分解起始溫度要低。所以在CASE 3中的高焓CO2風洞試驗中，CO2并未發生分解，整個流場均為CO2氣體。在CASE 3中的高焓空氣風洞試驗中，由于其來流條件與高焓CO2風洞試驗來流條件接近，激波層內的壓力與高焓CO2風洞試驗相差不大，略大于101 325 Pa；空氣主要成分是O2和N2，且O2和N2均發生了分解反應，O2甚至可能已經完全分解，因此從模型駐點沿壁面向肩部發展的流場區域內，氣體的性質主要以單原子氣體為主。

分析圖9和圖10可知，在飛行條件下，飛行器駐點附近到肩部區域激波層內溫度明顯高于6 000 K，而激波層內壓力小于101 325 Pa，因此飛行條件下的CO2起始分解溫度小于5 000 K。由于火星大氣主要成分以CO2為主，在飛行條件下的流場中CO2會發生分解反應。此時，從模型駐點沿壁面向肩部發展的流場區域內的氣體性質主要以單原子氣體為主。

因此，在保證雙尺度ρ∞L，在零攻角情況下的駐點處，均有

Cp,windtunnel=Cp,flight

(1)

(2)

式中：下標windtunnel表示風洞條件下的物理量；flight表示飛行條件下對應飛行器的物理量;“0”表示駐點處的物理量。

2.3 表面熱流試驗相關性

圖11是高焓風洞試驗模型壁面熱流分布比較。圖11(a)是兩種不同介質情況下在不同壁溫時的熱流分布，同時給出了駐點熱流的計算結果和試驗結果。可以看出，壁面條件采用高溫催化條件(Tw=2 000 K)時，計算結果與試驗結果峰值接近，最大相差0.12%。說明本文熱流計算方法合理。但是當采用冷壁條件(Tw=300 K)時，計算結果與試驗結果相差較大，其中，高焓CO2風洞試驗中，風洞模型的駐點熱流與試驗值相差18.7%；在高焓空氣風洞中，風洞模型的駐點熱流與試驗值相差7.4%。結合圖5，可以發現在高焓CO2風洞試驗中，冷壁條件下風洞模型的熱流值與試驗點的熱流之間的最大誤差為25%，最大誤差的幾何位置在肩部；熱壁條件下風洞模型的熱流值與試驗點的熱流之間的最大誤差為8%。在表1中，高焓CO2風洞和高焓空氣風洞中來流的溫度均比較高，模型壁面為高溫壁面。所以采用熱壁條件，得到的計算結果，比采用冷壁條件得到的計算結果，更接近于試驗值。所以，在數值計算過程中，應盡量模擬壁溫條件。

根據圖11(a)可知，在高焓空氣、高焓CO2風洞試驗中，風洞模型與飛行條件下對應飛行器的雙尺度一致時，以空氣為介質的風洞試驗熱流數據明顯大于以CO2為介質風洞試驗結果。圖11(b)是飛行條件下壁溫不同時得到的熱流分布，可以看出，冷壁條件下的熱流分布變化趨勢比較陡峭，熱壁條件下熱流分布相對平緩；在駐點附近，冷壁條件下的熱流明顯高于熱壁條件。圖11(c)是飛行條件與風洞條件下的無量綱熱流q/q0對比數據，可以看出，無量綱化之后，以空氣和CO2為介質的風洞試驗的熱流分布與飛行狀態下的熱流分布仍然不匹配。

圖12為高焓風洞與飛行條件Stanton數的分布。與飛行條件相比，高焓風洞條件下的Stanton數分布明顯高于飛行條件下飛行器的Stanton數 分布。高焓CO2風洞試驗在駐點附近區域的Stanton數，與飛行條件下的Stanton數吻合；從駐點沿壁面向后，兩者之間的差距逐漸增大，其中肩部區域附近最大相差27.3%。

圖11 壁面熱流分布Fig.11 Distribution of wall surface heat flow

圖12 壁面Stanton數的分布Fig.12 Distribution of wall Stanton number

根據上述分析可知：在數值過程中，應盡量模擬壁溫條件；高焓空氣、高焓CO2風洞試驗的熱流和無量綱熱流數據，都不宜直接作為將風洞條件外推至飛行條件的外推參數。在風洞中，不能直接將高焓空氣風洞試驗的Stanton數作為外推參數外推到飛行條件中；在高焓CO2風洞條件下，可以在風洞模型駐點附近區域內，將風洞試驗的Stanton數作為風洞試驗模型與飛行試驗下飛行器的相關性條件，且有：

St0,windtunnel≈St0,flight

(3)

3 結 論

通過對類“探路者號”外形的火星再入飛行器二維模型的計算，分析了風洞條件和飛行狀態的氣動力特性和氣動熱特性，在保證風洞條件下風洞模型與飛行條件下飛行器的雙尺度一致且攻角均為零的情況下，得到：

1) 在低焓空氣風洞試驗、高焓空氣風洞試驗和高焓CO2風洞試驗條件中，可以利用模型駐點附近區域內的無量綱壓力和壓力系數，將風洞條件與飛行條件相關聯起來。

2) 在高焓空氣和CO2風洞試驗條件下，不能直接將試驗熱流和無量綱熱流數據作為風洞模型與飛行條件下飛行器的關聯參數。

3) 在高焓空氣風洞試驗條件下，不能直接利用試驗所得Stanton數作為外推參數將飛行條件與風洞條件相關聯起來；在高焓CO2風洞試驗條件下，可以直接將模型駐點附近區域的Stanton數作為關聯參數，將風洞條件與飛行條件相關聯起來。

4) 在數值過程中，應盡量模擬壁溫條件，從而得到符合更接近于實際物理條件的計算結果。