火星大氣非平衡條件下駐點熱流計算公式的探討

劉方彬，袁軍婭

（北京航空航天大學(xué)宇航學(xué)院，北京100191）

1 引言

火星因為與地球具有很多類似的環(huán)境，成為人類空間探索的熱點?；鹦翘綔y器再入時，由于飛行速度很高，再入熱環(huán)境非常惡劣，因此氣動熱的計算十分重要，但是其熱環(huán)境的準(zhǔn)確計算目前存在較大的不確定性。然而，對于火星探測器的外形和軌道的快速確定及飛行器的優(yōu)化設(shè)計而言，駐點區(qū)域熱流的快速準(zhǔn)確計算又十分重要。目前可用于火星再入飛行器零攻角駐點熱流的計算公式有 Fay-Riddell公式［1］、Sutton-Grave 公式［2］、 Page 公 式［3］、 NASA 經(jīng) 驗 擬 合 公 式［4］、Edquist公式［5］及 Sutton-Grave 修正公式［6］等。 但是，這些公式均未考慮熱力學(xué)非平衡效應(yīng)，而火星再入飛行器再入過程中，駐點區(qū)域經(jīng)歷明顯的熱力學(xué)非平衡現(xiàn)象［7］，以上公式的適用性和準(zhǔn)確性值得探究。

對此，本文先討論典型的零攻角駐點熱流公式——Fay-Riddell公式在計算火星再入飛行器駐點熱流的局限性；然后針對Mars Pathfinder（詳見文獻7）外形和涵蓋飛行軌道的大范圍計算狀態(tài)，采用二維軸對稱熱化學(xué)非平衡數(shù)值計算方法對完全催化壁條件下的熱環(huán)境開展數(shù)值仿真；利用仿真得到的駐點熱流數(shù)據(jù)進行工程擬合，以得到適用于火星探測器再入飛行過程中連續(xù)流階段熱力學(xué)非平衡條件下零攻角駐點熱流的計算公式，并利用Mars Pathfinder飛行試驗數(shù)據(jù)對公式進行驗證。

2 Fay-Riddell公式的局限性

1958 年，F(xiàn)ay等［1］利用多羅德尼津-曼格勒變換和相似性假定，將高溫氣體邊界層偏微分方程轉(zhuǎn)化為常微分方程，和氣體熱力學(xué)特性、輸運特性有關(guān)的無因此參數(shù)假設(shè)為一系列常數(shù)，Pr＝0.71，Le＝1.0 ～2.0， ρsμs／ρwμw＝ 0.17 ～ 1.0。 利用薩瑟蘭黏性定律，在總焓 hs＝（1549 ～24158）kJ／kg，壁面溫度Tw＝ （300 ～ 3000）K 范圍內(nèi)，對駐點邊界層方程進行了數(shù)值計算，得出了式（1）所示完全催化壁條件下的駐點熱流估算公式［8］：

式中：qs為駐點熱流；Pr為普朗特常數(shù)；Le為路易斯數(shù)；ρ為密度，μ為動力粘性，h為焓值；下標(biāo)“s”和“w”分別代表駐點附近和壁面附近的物理量；hD為空氣離解焓；（due／dx）為駐點外緣速度梯度，在高超聲速流動情況下，可用式（2）所示修正牛頓流動公式計算［3］：

式中：下標(biāo)“∞ ”代表來流物理量，RN為駐點等效半徑；ps為駐點附近壓力。Fay-Riddell公式對于火星再入飛行器熱力學(xué)非平衡計算，有著很大的局限性：

1）Fay-Riddlell公式推導(dǎo)時，對于路易斯數(shù)等物理量的選取的時候，采用的是空氣的相關(guān)參數(shù)（詳細(xì)推導(dǎo)過程見文獻［1］），而火星大氣的主要成分是CO2，CO2的物理性質(zhì)和空氣的物理性質(zhì)從粘性、密度等物理量上都有差別［9］。

2）Fay-Riddell公式在理論推導(dǎo)過程中未考慮熱力學(xué)非平衡的影響，在路易斯數(shù)和 ρsμs／ρwμw擬合時采用的數(shù)值計算也未考慮到熱力學(xué)非平衡的影響，僅有化學(xué)非平衡的影響。而火星大氣環(huán)境中，火星大氣密度和壓力都比地球小兩個量級，溫度也比地球平均溫度低，且CO2本身的氣體分子振動特征也比空氣低。因此，在火星環(huán)境下的高超速流動呈現(xiàn)低雷諾數(shù)、高馬赫數(shù)的特點，熱力學(xué)非平衡效應(yīng)明顯，在計算中需要考慮到CO2的振動溫度［10］。因此，在火星大氣條件下，駐點熱流計算采用熱力學(xué)非平衡計算比熱化學(xué)非平衡計算計算更合適［11］。

采用Fay-Riddell公式對Mars Pathfinder飛行器典型飛行狀態(tài)［8，11］進行計算，結(jié)果如圖 1、表 1所示，可以看出，F(xiàn)ay-Riddle公式在高空熱力學(xué)非平衡效應(yīng)非常明顯的情況下計算誤差較大，在56 km時誤差為33%，85 km時誤差達(dá)到40%。

圖1 Sutton-Grave公式、Fay-Riddell公式和飛行試驗對比Fig.1 Comparison of Sutton-Grave formula fay-Riddell formula and flight experiment

因此，為了更好的估算火星大氣條件下再入飛行器的駐點熱流，許多研究人員做了不同的假設(shè)，采用了不同的方法，其中Sutton等人在1971年根據(jù)激波管試驗，擬合出了式（3）所示CO2環(huán)境下的 Sutton-Grave 公式［2］：

但是Sutton在擬合公式過程中，相應(yīng)的物理量來源于激波管風(fēng)洞試驗，而風(fēng)洞條件與真實大氣條件下的飛行所得到的參數(shù)存在較大差異。采用該公式對飛行條件下駐點熱流進行計算，如圖1、表2所示，在高空和低空時均存在較大誤差，56km時誤差為24%，21 km時誤差為75%。

表1 Fay-Riddell公式計算對比結(jié)果Table 1 Comparison of Fay-Riddell formula calculation with the results in literature

表2 Sutton-Grave公式計算對比結(jié)果Table 2 Comparison of Sutton-Grave Formulation calculation with the results in literature

由以上結(jié)果可以看出，應(yīng)對火星大氣條件下的再入飛行器的駐點熱流公式進行數(shù)值仿真和理論分析，以得到適用于飛行條件下的零攻角駐點熱流估算公式。

3 數(shù)值計算方法

3.1 控制方程和數(shù)值方法

考慮零攻角情況，控制方程采用軸對稱的化學(xué)振動非平衡的流動的無量綱化方程如式（4）：

其中：

這里i＝1，2，…，Ns，Ns為組分個數(shù)；t、ρ、p、u、v分別為時間、密度、壓力、相對于坐標(biāo)x和y的速度分量；ρi、fi和Di分別為組分i的密度、質(zhì)量分?jǐn)?shù)和擴散系數(shù)；Ms、hs、evs為組分s的摩爾質(zhì)量、焓值和振動能（J·mol－1）； E、H0、Ev分別為總能量、焓、和振動能；qv和q分別為振動熱流和總熱流；Re為雷諾數(shù)。部分量有下列關(guān)系：H0＝E＋e為總內(nèi)能；μ∞和μ黏性系數(shù)的量綱和無量綱后的值；T和Tv為輸運溫度和振動溫度；ktr和kv為輸運和振動的熱傳導(dǎo)系數(shù)。對于化學(xué)振動非平衡，狀態(tài)方程如式（5）：

本次數(shù)值離散采用簡單隱式全耦合TVD格式的計算方法。壁面條件采用滿足無滑移條件（u＝v＝0），等溫壁條件（T ＝Tv＝Tw），壓力法向?qū)?shù)為外邊界采用激波捕捉方法，所有變量均為來流值。對稱軸滿足和v＝0。出口原始變量滿足線性關(guān)系。以上所有和均采用二階差分離散。

計算網(wǎng)格采用代數(shù)方法生成。為了是網(wǎng)格在物面附近加密以模擬附面層，采用了指數(shù)壓縮函數(shù)來加密網(wǎng)格。

3.2 火星熱化學(xué)動力學(xué)模型

火星大氣由97%的CO2和3%的N2組成?；瘜W(xué)反應(yīng)模型考慮v-d耦合效應(yīng)，假設(shè)速率常數(shù)為控制溫度Tf和Tb的函數(shù)，正逆反應(yīng)速率常數(shù)有式（6）～（7）所示關(guān)系［2］：

這里 i＝1，2，……，Nr、Nr為反應(yīng)數(shù)；A 為指前因子，B為溫度系數(shù)，C活化能，f為正向反應(yīng)，b為逆向反應(yīng)；Af，i，Bf，i，C，fi見文獻［12］，逆反應(yīng)常數(shù)kb，i通過平衡狀態(tài)獲得。本次化學(xué)反應(yīng)為CO2的8組分9化學(xué)反應(yīng)模型，如表3所示，化學(xué)反應(yīng)速率見文獻［11］。對于第i個反應(yīng)對于第i個反應(yīng)的生成化學(xué)源項如式（8）［2］：

第 j種組分的生成化學(xué)源項如式（9）［2］：

表3 8組分9化學(xué)反應(yīng)［13］Table 3 Eight Species Nine Reactions Mechanism［13］

對于CO2，由于有三個振動模態(tài)，且每個模態(tài)對于不同的振動特征溫度和特征時間。為了減少控制方程數(shù)量，在振動模型中，振動非平衡效應(yīng)只考慮一個振動溫度。則振動能量如式（10）～（11）：

這里θvs是組分s的振動特征問。振動非平衡能量源項Wv為兩項之和，一項是輸運能與振動能的交換項，另一項是化學(xué)反應(yīng)中分子理解的振動能量的變化項，Wv的計算公式為式（12）［3］：

v-d耦合效應(yīng)在本文中采用簡單的Park耦合模型。在該模型中，用控制溫度的概念來描述振動非平衡相應(yīng)對離解的影響，即離解速率用控制溫度計算如式（13） ～ （14）：

這里的α一般取值范圍為0.5≤α≤0.7［3］，本次計算取值為常用的經(jīng)驗值，即0.5。壁面附近化學(xué)反應(yīng)與流場中化學(xué)反應(yīng)相同，其中壁面處考慮完全催化（fi＝0，fi為單原子）。

3.3 數(shù)值方法驗證

采用Mars Pathfinder的外形進行驗證，氣體為97%的CO2和3%的N2，驗證標(biāo)準(zhǔn)為駐點熱流。對于熱流數(shù)值計算，網(wǎng)格的影響主要體現(xiàn)在法向網(wǎng)格間距上，這里網(wǎng)格無關(guān)性分析著重考察法向網(wǎng)格間距的影響。網(wǎng)格如圖2所示，計算條件為驗證算例工況：Ma＝32，T∞＝ 169 K，Tw＝2000 K，ρ∞＝ 2.8×10－4kg·m3。 計算結(jié)果如表 4。

圖2 計算網(wǎng)格Fig.2 Computational grid

表4 網(wǎng)格無關(guān)性驗證Table 4 Validation of grid independence

根據(jù)表4可以發(fā)現(xiàn)，熱流計算對于最小網(wǎng)格間距要求極其嚴(yán)格。因此，在每次更換狀態(tài)計算時，均需要按照表4方法驗證網(wǎng)格無關(guān)性，直至計算駐點熱流幾乎不再變化，以保證計算結(jié)果的準(zhǔn)確性。

4 駐點熱流計算公式擬合

4.1 擬合公式

本章利用數(shù)值計算得到的駐點熱流數(shù)據(jù)，重新進行公式擬合，以得到適性更廣、準(zhǔn)確性更好的駐點熱流計算公式。

借鑒 Fay-Riddell公式［1］、Sutton-Grave 公式［7］、Sutton-Gave 修正公式［7］等駐點熱流計算公式及其推導(dǎo)過程，可以得出駐點與各個變量之間的關(guān)系如式（15）：

式中： c1，c2，c3為常數(shù)。 結(jié)合 Sutton-Grave 公式及修正公式，認(rèn)為Sutton-Gave公式的適用范圍較窄的原因在于將壁面溫度影響忽略，因此考慮壁面因素的影響，擬采用的公式形式為式（16）：

其中hs、hw和h∞分別為駐點焓值、壁面焓值和來流焓值。數(shù)值計算的工況以Mars Pathfinder飛行器的軌道數(shù)據(jù)為基礎(chǔ)，如表5所示［14］。

表5 Pathfinder典型飛行參數(shù)［14］Table 5 Typical flight data of Mars Pathfinder［14］

4.2 駐點熱流與速度的關(guān)系

為了建立熱流與速度之間的關(guān)系，選取的速度的變化范圍為 3000 m／s、4000 m／s、5000 m／s、5333 m／s、6000 m／s、6041 m／s、6500 m／s、6596 m／s、7000 m／s、8000 m／s和9000 m／s；其余遠(yuǎn)場物理量采用表5中狀態(tài)3和狀態(tài)4的物理參量。數(shù)據(jù)擬合時采用線性回歸進行擬合，可以得到速度隨駐點熱流變化曲線如圖3。其中縱軸C＝logqs－k為線性回歸時擬合到的斜率，來流物理參量為表5狀態(tài)4時進行的熱流計算和擬合。

圖3 速度-駐點熱流變化Fig.3 Relationship between stagnation heat flux＆velocity

根據(jù)圖3，可以發(fā)現(xiàn)駐點熱流和速度的2.89次方成線性關(guān)系如式（17）：

4.3 駐點熱流與半徑的關(guān)系

為了建立熱流與半徑之間的變化關(guān)系，來流參量采用表5中狀態(tài)3的物理量。形狀為圖1中的形狀，半徑從0.01R到1.4R之間選取95個值，其中R＝0.6638 m。數(shù)據(jù)擬合時采用線性回歸進行擬合，可以得到速度與駐點熱流之間的變化關(guān)系。其中圖 4縱軸 C1＝ logqs－為線性回歸時擬合得到的斜率。根據(jù)圖4，可以得到駐點熱流與半徑之間的關(guān)系如式（18）：

4.4 駐點熱流與密度的關(guān)系

密度選取表5中的狀態(tài)3、狀態(tài)4和狀態(tài)5的高度進行變化。得到的熱流與密度之間的變化關(guān)系如圖5所示。數(shù)據(jù)擬合時采用線性回歸。其中圖5中縱軸為線性回歸時擬合到的斜率。

圖4 熱流-半徑變化Fig.4 Relationship between stagnation heat flux＆radius

圖5 熱流-密度關(guān)系Fig.5 Relationship between heat flux＆density

根據(jù)圖5，可以發(fā)現(xiàn)駐點熱流與密度之間的關(guān)系如式（19）所示：

4.5 熱力學(xué)非平衡下零攻角駐點熱流公式

建立熱流值與速度、密度和和半徑的關(guān)系后，進一步采用線性回歸進行擬合，如圖6所示，可以得到熱流的擬合公式。圖6中橫軸 C3＝縱軸

圖6 線性關(guān)系式Fig.6 Linear formula

根據(jù)圖6可知，熱流密度與來流速度、來流密度和等效半徑存在式（20）所示一定關(guān)系：

式中：故根據(jù)前面計算，得到零攻角駐點熱流公式為式（21）：

4.6 驗證擬合公式的壁溫修正

壁溫分別取值 1000 K、1500 K、1600 K、2000 K、2100 K、2500 K、3000 K、3500 K、4000 K、4500 K和5000 K，其余參數(shù)選取表5中的狀態(tài)3，驗證公式（21）中壁溫修正的可行性。得到的熱流與密度之間的變化關(guān)系如圖7所示。圖7中qs為駐點熱流，單位 MW／m2，橫坐標(biāo) C5為

進行線性回歸求得趨勢圖，可以得到 qs＝1.8688×C5，而公式（21）在表5中狀態(tài)3下熱流與壁溫修正的系數(shù)為106＝1.86448，相差2.3%。因此，可以接受公式（21）的壁面修正。

同時也可以發(fā)現(xiàn)，壁面溫度的取值對于駐點熱流的影響也很突出。因此，在駐點熱流估算中，也應(yīng)該考慮到壁面溫度對于駐點熱流的影響。

圖7 壁面溫度修正驗證Fig.7 Verification of wall temperature correction

4.7 驗證擬合公式

采用擬合得到的計算公式（21）對表1中的工況進行計算，與文獻給出的飛行試驗數(shù)據(jù)［14］進行對比，見表6。

表6 擬合公式計算對比結(jié)果Table 6 Comparison of fitting formula computation with flight results

從表6可以看出，在雙溫?zé)崃W(xué)非平衡高溫完全催化壁條件Tw≥1000 K情況下，公式（21）與文獻試驗數(shù)據(jù)符合的很好，最大誤差為6.25%，可以認(rèn)為公式（21）能夠很好的運用于火星探測器再入過程中的高溫完全催化壁面的計算，而且公式（21）的計算高度范圍也比Fay-Riddell公式和SuttonGrave公式要寬許多。

5 結(jié)論

本文通過數(shù)值仿真計算和工程擬合得到的駐點熱流公式相對于Fay-Riddell公式，更適合于火星大氣條件下的零攻角駐點熱流計算：

1）該公式在擬合過程中使用的數(shù)據(jù)均是在熱力學(xué)非平衡狀態(tài)下的數(shù)據(jù)，比Fay-Riddell公式更適用于火星條件下的熱力學(xué)非平衡計算過程。

2）該公式采用遠(yuǎn)場焓作為基準(zhǔn)進行壁面修正，比Fay-Riddell的適用高度范圍更為寬廣。

3）采用Mars Pathfinder的飛行試驗數(shù)據(jù)對比表明，該公式的準(zhǔn)確性更高。

4）由于壁面溫度對于駐點熱流的計算影響很大，因此，需要在駐點熱流公式中考慮到壁面溫度對于駐點熱流的影響。因此，該公式也比Sutton-Grave公式的適用高度范圍更寬。

該公式目前仍有不足之處，主要是擬合數(shù)據(jù)有限，因為目前可查到的有真實飛行試驗數(shù)據(jù)的零攻角再入的火星探測器僅有 Pathfinder。隨著后續(xù)火星探測器的數(shù)據(jù)的不斷增加，擬合公式還可以進一步驗證和細(xì)化。