劉方彬, 袁軍婭
(1. 北京航空航天大學宇航學院, 北京 100083; 2. 深圳易信科技股份有限公司, 深圳 518000)
火星大氣的主要成分是CO2,大氣非常稀薄,密度僅為地球大氣的1%。相對于地球大氣,再入飛行器進入火星大氣時,雖然同樣具有較高的飛行速度,但是由于火星大氣密度較低,在高速飛行時激波后流動的熱力學非平衡特性更為突出。作為火星大氣主要成分的CO2氣體,具有和N2、O2等氣體不同的熱化學特性[1]。目前對于火星再入飛行器實際飛行時,仍存在諸多認知上的不足,包括剪切層的復雜性、尺度問題和一些現(xiàn)象的物理模型等方面,地面試驗仍是火星再入飛行器研制中必不可少的驗證手段[2]。
由于風洞條件的限制,難以完全模擬飛行環(huán)境,且在風洞試驗中有尺寸效應、邊界效應或邊界干擾、支架干擾等干擾,這就存在一個與真實飛行條件下的相似理論及相關性的問題。目前關于地球再入飛行器的風洞試驗相關性研究較多,相似準則的選擇也有多種,而關于火星探測器風洞試驗相關性研究的國內外公開文獻極少,且成熟度不高。例如苗文博等[1]則利用數(shù)值方法研究了火星再入飛行器表面熱環(huán)境,通過研究認為,在駐點附近區(qū)域近似為熱力學平衡,隨著高度逐漸增加,則熱力學非平衡效應更明顯。董維中[3]在博士論文中根據(jù)鈍錐標模ELECTRE模型風洞試驗和飛行試驗數(shù)據(jù)之間對比,認為保持總焓值不變,利用雙尺度參數(shù)可以對模型頭部進行氣動熱/力的外推。然而ELECTRE模型的試驗介質均為空氣,沒有進行介質為CO2的試驗。Bur等[4]通過研究風洞與“Pathfinder”(“探路者號”)的計算數(shù)據(jù),認為風洞條件下駐點和前體部分熱流計算峰值要高28%。Armenise等[5]則對風洞條件和飛行條件下的計算化學反應模型進行了分析和驗證。Paterna等[6]通過試驗和計算,認為火星再入飛行器在CO2氣體環(huán)境下的計算,采用催化壁條件獲得到的物理參數(shù)更為合適。
鑒于風洞試驗在火星再入飛行試驗中的重要性,需要建立風洞條件與飛行條件之間的相關性?;诒姸嗫蒲泄ぷ髡叩拈_拓性研究,本文利用了“Pathfinder”風洞試驗數(shù)據(jù)[7-8]和飛行試驗數(shù)據(jù)[9-10],通過數(shù)值方法,建立了火星再入飛行器氣動力和氣動熱的風洞參數(shù)外推方法。
本次計算采用的是“Pathfinder”火星再入飛行器的飛行數(shù)據(jù)[11]和風洞試驗數(shù)據(jù)[12],其中風洞數(shù)據(jù)來源于“MP-1”模型的風洞試驗數(shù)據(jù)[12],飛行數(shù)據(jù)來源于“Pathfinder”的飛行數(shù)據(jù)[10]。兩者前體幾何相似,縮尺比例為1|∶ 52.15,如圖1[11]和圖2[12]所示,其中圖2中1 in= 0.025 4 m。表1[12]中CASE 1和CASE 2是兩種常規(guī)高超聲速空氣風洞試驗條件及其對應的飛行狀態(tài)。CASE 3是試驗介質分別為空氣、CO2的高焓風洞試驗數(shù)據(jù)及其對應的飛行數(shù)據(jù)。在本文計算中,認為火星大氣其主要成分是由體積分數(shù)97%的CO2和3%的O2組成;空氣的主要成分由體積分數(shù)79%的N2和21%的O2組成;來流風洞的主要成分由體積分數(shù)為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∞為來流馬赫數(shù),Re∞為來流雷諾數(shù)。“97%CO2+3%N2”代表來流介質的主要組成成分是體積分數(shù)97%的CO2和3%的N2,即真實飛行條件中的火星大氣,“21%O2+79%N2”代表來流介質的主要組成成分是體積分數(shù)21%的O2和79%的N2,即空氣風洞試驗中的空氣,“100%CO2”代表來流介質的主要組成成分是體積分數(shù)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]
求解的控制方程為采用化學反應全Navier-Stokes方程求解二維軸對稱外形,并基于以下假設:流動為熱力學非平衡;忽略輻射和徹體力的影響;流動質量采用雙組元氣體模型;溫度模型采用Park雙溫模型[13];壁面采用完全催化壁,壁面單原子和離子濃度為0;常溫風洞條件下壁面溫度設置為300 K,高焓風洞條件和飛行條件下的壁溫設為2 000 K。
對于空氣的化學反應模型,采用Park 5組分化學反應模型[13],如表2所示。對于CO2的化學反應模型,則采用8組分9化學反應模型[14],如表3所示。
采用有限體積法對控制方程進行數(shù)值求解,空間格式采用二階TVD格式,黏性項采用二階中心差分,時間格式采用隱式求解。
計算網(wǎng)格采用代數(shù)關系方法生成,為了模擬附面層,壁面網(wǎng)格采用了線性函數(shù)來加密。由于熱流是以黏性為主導的物理量,高超聲速飛行器的氣動熱計算對網(wǎng)格依賴性較高。根據(jù)文獻[15],在高超聲速中,網(wǎng)格雷諾數(shù)小于8可以獲得收斂的熱流結果[15]。所以在本文計算中,所有網(wǎng)格雷諾數(shù)均小于8?!癕P-1”與“Pathfinder”網(wǎng)絡如圖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”的網(wǎng)格Fig.3 Mesh of “MP-1”
圖4 “Pathfinder”的網(wǎng)格Fig.4 Mesh of “Pathfinder”
通過與NASA Langley研究中心“Pathfinder”飛行狀態(tài)計算結果[9-12]的對比,來驗證本文所采用的熱力學模型和計算方法能否用于飛行條件下的飛行器流場計算;采用表1中CASE 3的CO2條件開展計算,將計算結果與風洞試驗結果[12]進行對比,來驗證本文所采用的熱力學模型和計算方法能否運用于風洞條件下的模型流場計算。
圖5為采用本文方法計算得到的風洞條件下的模型表面熱流與風洞試驗數(shù)據(jù)[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]中的風洞條件下的試驗數(shù)據(jù)(下文圖中的曲線標識類似)。根據(jù)參考文獻[6],在飛行條件下的飛行器流場計算中,壁面條件采用完全催化壁。從圖5可以看出,風洞條件下的熱流峰值計算q0,CFD與試驗q0,Experiment相差0.1%,試驗點上最大的誤差為8%,證明所采用的熱力學模型和計算方法在風洞條件下的計算精度可以滿足分析要求。
表4為采用本文的計算方法得到的飛行狀態(tài)駐點熱流與文獻[8]中計算結果的對比。從表4可以看出,本文的計算方法在飛行狀態(tài)時的計算結果與NASA Langley中心計算結果最大誤差為6%,表明本文的計算方法對飛行條件下的計算精度滿足分析要求。
表4 對比條件和計算結果(飛行條件)Table 4 Comparison conditions and computation results(flight condition)
由于低焓風洞中流場溫度較低,所以在“MP-1”計算中僅考慮冷壁(壁溫Tw=300 K)非催化壁條件;在高焓風洞條件下,來流溫度較高,為了比較冷壁和熱壁(Tw=2 000 K)條件對相關性的影響,在“MP-1”計算中分別采用冷壁和熱壁2 種計算條件,以便對比分析火星再入飛行器飛行狀態(tài)下冷壁和熱壁對相關性的影響。
圖6 風洞條件和飛行條件壓力系數(shù)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)中可以看出,對于低焓風洞試驗的壓力系數(shù)和無量綱壓力,在駐點附近時與飛行數(shù)據(jù)相差很??;從駐點沿壁面向肩部發(fā)展,低焓風洞試驗結果與飛行數(shù) 據(jù)之間差距逐漸增大,最大相差11.5%。壁面溫度對飛行狀態(tài)下的壓力系數(shù)和無量綱壓力沒有明顯的影響。
圖8 風洞條件和飛行條件的馬赫數(shù)對比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)可以發(fā)現(xiàn),對于高焓風洞試驗壓力系數(shù)和無量綱壓力,以及高焓CO2風洞試驗壓力系數(shù)與無量綱壓力,除了模型駐點附近區(qū)域以外,其他地方與飛行狀態(tài)下的壓力系數(shù)和無量綱壓力有明顯的差距;對于高焓空氣風洞試驗的壓力系數(shù)和無量綱壓力,從駐點沿壁面向肩部發(fā)展,與飛行條件下的壓力系數(shù)和無量綱壓力幾乎吻合。高焓空氣風洞試驗無量綱壓力和壓力系數(shù),比常溫空氣風洞試驗的無量綱壓力和壓力系數(shù)更接近于飛行條件下的無量綱壓力與壓力系數(shù)。
在低焓風洞條件與真實飛行之間出現(xiàn)以上情況的主要原因是馬赫數(shù)影響。由于常溫下空氣的比熱比與火星大氣比熱比接近,根據(jù)馬赫數(shù)無關原理,對于高超聲速飛行,當馬赫數(shù)較高時,壓力系數(shù)基本上與馬赫數(shù)無關,僅與流動介質的比熱比相關。根據(jù)表1可以發(fā)現(xiàn),不管CASE 1和CASE 2 的風洞條件還是飛行條件,來流馬赫數(shù)都較高,因此低焓風洞試驗模型的駐點處的壓力系數(shù)與飛行條件下的駐點壓力系數(shù)接近。根據(jù)圖8(a)可以知道,從模型駐點處沿著壁面向肩部發(fā)展,激波層內的馬赫數(shù)逐漸降低,模型肩部附近馬赫數(shù)較小。根據(jù)牛頓修正公式[16],當馬赫數(shù)逐漸減小時,壓力系數(shù)等參數(shù)逐漸減小,而且當比熱比越大,壓力系數(shù)的變化幅度也越大。所以當壓力系數(shù)和無量綱壓力逐漸從模型駐點沿壁面向肩部發(fā)展而逐漸減小,且空氣介質中的減小幅度大于CO2介質中的減小幅度。
分析高焓風洞試驗和真實飛行的流場,根據(jù)圖8(b)可知,兩者來流馬赫數(shù)都比較高,所以引起壓力系數(shù)和無量綱壓力變化差異的原因主要在于溫度。標準大氣壓下,CO2在5 000 K時完全分解,且CO2分解的溫度域很小。O2在2 000 K左右開始分解,4 000 K左右時完全分解,N2在3 000 K左右開始分解,9 000 K幾乎完全分解。從圖9和圖10可知,在高焓CO2風洞試驗中,模型駐點到肩部區(qū)域的激波層內溫度不超過4 000 K,模型頭部駐點附近激波層內的壓力略大于101 325 Pa。根據(jù)高超聲速理論,除了振動激勵現(xiàn)象不受壓力影響以外,如果壓力上升,其他起始溫度將會對應上升[16]。即壓力大于101 325 Pa時,氣體的分解起始溫度比101 325 Pa下的分解起始溫度要高;當壓力小于101 325 Pa時,氣體的分解起始溫度溫度比101 325 Pa下的分解起始溫度要低。所以在CASE 3中的高焓CO2風洞試驗中,CO2并未發(fā)生分解,整個流場均為CO2氣體。在CASE 3中的高焓空氣風洞試驗中,由于其來流條件與高焓CO2風洞試驗來流條件接近,激波層內的壓力與高焓CO2風洞試驗相差不大,略大于101 325 Pa;空氣主要成分是O2和N2,且O2和N2均發(fā)生了分解反應,O2甚至可能已經(jīng)完全分解,因此從模型駐點沿壁面向肩部發(fā)展的流場區(qū)域內,氣體的性質主要以單原子氣體為主。
分析圖9和圖10可知,在飛行條件下,飛行器駐點附近到肩部區(qū)域激波層內溫度明顯高于6 000 K,而激波層內壓力小于101 325 Pa,因此飛行條件下的CO2起始分解溫度小于5 000 K。由于火星大氣主要成分以CO2為主,在飛行條件下的流場中CO2會發(fā)生分解反應。此時,從模型駐點沿壁面向肩部發(fā)展的流場區(qū)域內的氣體性質主要以單原子氣體為主。
因此,在保證雙尺度ρ∞L,在零攻角情況下的駐點處,均有
Cp,windtunnel=Cp,flight
(1)
(2)
式中:下標windtunnel表示風洞條件下的物理量;flight表示飛行條件下對應飛行器的物理量;“0”表示駐點處的物理量。
圖11是高焓風洞試驗模型壁面熱流分布比較。圖11(a)是兩種不同介質情況下在不同壁溫時的熱流分布,同時給出了駐點熱流的計算結果和試驗結果??梢钥闯?,壁面條件采用高溫催化條件(Tw=2 000 K)時,計算結果與試驗結果峰值接近,最大相差0.12%。說明本文熱流計算方法合理。但是當采用冷壁條件(Tw=300 K)時,計算結果與試驗結果相差較大,其中,高焓CO2風洞試驗中,風洞模型的駐點熱流與試驗值相差18.7%;在高焓空氣風洞中,風洞模型的駐點熱流與試驗值相差7.4%。結合圖5,可以發(fā)現(xiàn)在高焓CO2風洞試驗中,冷壁條件下風洞模型的熱流值與試驗點的熱流之間的最大誤差為25%,最大誤差的幾何位置在肩部;熱壁條件下風洞模型的熱流值與試驗點的熱流之間的最大誤差為8%。在表1中,高焓CO2風洞和高焓空氣風洞中來流的溫度均比較高,模型壁面為高溫壁面。所以采用熱壁條件,得到的計算結果,比采用冷壁條件得到的計算結果,更接近于試驗值。所以,在數(shù)值計算過程中,應盡量模擬壁溫條件。
根據(jù)圖11(a)可知,在高焓空氣、高焓CO2風洞試驗中,風洞模型與飛行條件下對應飛行器的雙尺度一致時,以空氣為介質的風洞試驗熱流數(shù)據(jù)明顯大于以CO2為介質風洞試驗結果。圖11(b)是飛行條件下壁溫不同時得到的熱流分布,可以看出,冷壁條件下的熱流分布變化趨勢比較陡峭,熱壁條件下熱流分布相對平緩;在駐點附近,冷壁條件下的熱流明顯高于熱壁條件。圖11(c)是飛行條件與風洞條件下的無量綱熱流q/q0對比數(shù)據(jù),可以看出,無量綱化之后,以空氣和CO2為介質的風洞試驗的熱流分布與飛行狀態(tài)下的熱流分布仍然不匹配。
圖12為高焓風洞與飛行條件Stanton數(shù)的分布。與飛行條件相比,高焓風洞條件下的Stanton數(shù)分布明顯高于飛行條件下飛行器的Stanton數(shù) 分布。高焓CO2風洞試驗在駐點附近區(qū)域的Stanton數(shù),與飛行條件下的Stanton數(shù)吻合;從駐點沿壁面向后,兩者之間的差距逐漸增大,其中肩部區(qū)域附近最大相差27.3%。
圖11 壁面熱流分布Fig.11 Distribution of wall surface heat flow
圖12 壁面Stanton數(shù)的分布Fig.12 Distribution of wall Stanton number
根據(jù)上述分析可知:在數(shù)值過程中,應盡量模擬壁溫條件;高焓空氣、高焓CO2風洞試驗的熱流和無量綱熱流數(shù)據(jù),都不宜直接作為將風洞條件外推至飛行條件的外推參數(shù)。在風洞中,不能直接將高焓空氣風洞試驗的Stanton數(shù)作為外推參數(shù)外推到飛行條件中;在高焓CO2風洞條件下,可以在風洞模型駐點附近區(qū)域內,將風洞試驗的Stanton數(shù)作為風洞試驗模型與飛行試驗下飛行器的相關性條件,且有:
St0,windtunnel≈St0,flight
(3)
通過對類“探路者號”外形的火星再入飛行器二維模型的計算,分析了風洞條件和飛行狀態(tài)的氣動力特性和氣動熱特性,在保證風洞條件下風洞模型與飛行條件下飛行器的雙尺度一致且攻角均為零的情況下,得到:
1) 在低焓空氣風洞試驗、高焓空氣風洞試驗和高焓CO2風洞試驗條件中,可以利用模型駐點附近區(qū)域內的無量綱壓力和壓力系數(shù),將風洞條件與飛行條件相關聯(lián)起來。
2) 在高焓空氣和CO2風洞試驗條件下,不能直接將試驗熱流和無量綱熱流數(shù)據(jù)作為風洞模型與飛行條件下飛行器的關聯(lián)參數(shù)。
3) 在高焓空氣風洞試驗條件下,不能直接利用試驗所得Stanton數(shù)作為外推參數(shù)將飛行條件與風洞條件相關聯(lián)起來;在高焓CO2風洞試驗條件下,可以直接將模型駐點附近區(qū)域的Stanton數(shù)作為關聯(lián)參數(shù),將風洞條件與飛行條件相關聯(lián)起來。
4) 在數(shù)值過程中,應盡量模擬壁溫條件,從而得到符合更接近于實際物理條件的計算結果。