DSMC量子動理學(xué)模型在火星再入流動中的應(yīng)用

李錦，耿湘人，陳堅強(qiáng)，江定武，李紅喆

1. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000 2. 中國空氣動力研究與發(fā)展中心 設(shè)備設(shè)計及測試技術(shù)研究所，綿陽 621000

火星探測器在進(jìn)入、下降和著陸(Entry, Descent and Landing, EDL)過程中會遭遇氣動特性方面的挑戰(zhàn)[1]。相比于地球，火星大氣非常稀薄，大約只有地球大氣密度的1%,火星表面的氣體密度大致相當(dāng)于地球30 km高度的大氣密度。與地球再入飛行器不同，火星探測器的減速高度較低，可能達(dá)不到最終的亞聲速平衡下降速度，等到速度低至可以釋放減速裝置時，探測器距離火星表面已經(jīng)非常近，留給著陸過程的時間十分有限[2]。圖1給出了火星和地球典型EDL過程的比較。一方面，稀薄的火星大氣層會導(dǎo)致探測器減速不夠；而另一方面，火星大氣層并沒有稀薄到可以忽視氣動加熱的程度。對于沒有主動控制系統(tǒng)的探測器，稀薄氣動特性的準(zhǔn)確預(yù)測將是一個重要的問題。

Moss等[2]利用直接模擬蒙特卡羅(DSMC)計算代碼G2和DAC以及自由分子流計算軟件，詳細(xì)計算了探路者號火星探測器的稀薄氣動特性，計算條件從自由分子流一直延續(xù)到連續(xù)流條件，攻角從0°～30°。結(jié)果顯示探路者號在稀薄-過渡流動大部分區(qū)域都是靜不穩(wěn)定的。Edquist等[3]根據(jù)探路者號火星車的測試分析數(shù)據(jù)并使用Navier-Stokes(N-S)方程解算器LAURA的計算結(jié)果，建立了鳳凰號火星探測器的氣動力數(shù)據(jù)庫。結(jié)果表明在高超聲速過渡流以及近極限動壓情況下，鳳凰號存在靜俯仰不穩(wěn)定性。同時，非平衡流體力學(xué)效應(yīng)還將導(dǎo)致探測器配平攻角非零，在近極限動壓情況下，配平攻角可達(dá)2.5°。

國內(nèi)呂俊明等[4]利用N-S 方程解算器分析了火星科學(xué)實驗室號進(jìn)入火星大氣時探測器周圍的流場結(jié)構(gòu)、化學(xué)非平衡效應(yīng)和氣體特性變化規(guī)律以及火星探路者號小攻角飛行的靜不穩(wěn)定性[5]。黃飛等[6]計算分析了火星稀薄大氣真實氣體效應(yīng)對氣動特性的影響。程曉麗等[7]針對火星再入環(huán)境，發(fā)展了類探路者號外形在不同火星軌道點的氣動力特性預(yù)估方法。梁杰等[8]模擬了火星科學(xué)實驗室號在火星大氣環(huán)境70 km高度、進(jìn)入速度為5.85 km/s下的高溫真實氣體效應(yīng)對氣動力、氣動熱和流場特征的影響。

量子動理學(xué)(Quantum Kinetic, QK)模型是DSMC方法中新近發(fā)展的一種化學(xué)反應(yīng)模型[9]。QK模型已經(jīng)成功吸引了許多學(xué)者的關(guān)注，在國內(nèi)相關(guān)研究工作也正在開展。陳浩等[10]評估了QK模型在氮氧離解復(fù)合反應(yīng)中的表現(xiàn)。楊超和孫泉華[11]針對高溫氣體熱化學(xué)反應(yīng)對比分析了典型化學(xué)反應(yīng)模型(TCE(Total Collision Energy)、VFD(Vibrationally-Favored Dissociation)、QK)的微觀動力學(xué)性質(zhì)。在本文作者之前的研究中[12]，利用絕熱單元熱泳平衡/非平衡測試算例，仔細(xì)評估了該模型在化學(xué)反應(yīng)速率以及化學(xué)反應(yīng)熱松弛過程的預(yù)測能力，研究結(jié)果表明，QK模型具有優(yōu)良的特性，可應(yīng)用于高超聲速化學(xué)反應(yīng)流動計算，QK模型不依賴宏觀的化學(xué)反應(yīng)速率數(shù)據(jù)，可嘗試在深空探測等難以獲取可靠實驗數(shù)據(jù)的領(lǐng)域應(yīng)用。本文在之前研究的基礎(chǔ)上，將量子動理學(xué)模型進(jìn)一步應(yīng)用于火星探測器再入流動稀薄氣動特性的數(shù)值模擬，分析了流場結(jié)構(gòu)及壁面物理量的分布情況，以進(jìn)一步評估稀薄氣體效應(yīng)和化學(xué)反應(yīng)對氣動特性的影響規(guī)律以及QK模型的性能。

1 火星大氣及化學(xué)反應(yīng)模型

火星大氣成分為95%的CO2、3%的N2、1.6% 的Ar，很少的O2、H2O等。圖2給出了火星和地球大氣層密度ρ與溫度T隨高度h變化的比較。

圖2 火星和地球大氣參數(shù)對比Fig.2 Comparison of atmosphere parameters on the Mars and the Earth

與地球化學(xué)反應(yīng)類似，火星化學(xué)反應(yīng)同樣可以采用TCE模型或者QK模型計算，只不過反應(yīng)組份要從氮氧系列的五組份擴(kuò)充到包含碳系列的八組份。CO2是線性三原子分子，有3個平動自由度、2個轉(zhuǎn)動自由度和4個振動模態(tài)[13-15]。振動模態(tài)中，第1和第2模態(tài)為退化的彎曲模態(tài)，特征振動溫度為945 K。第3模態(tài)為對稱拉伸模態(tài)，特征振動溫度為1 903 K。第4模態(tài)為反對稱拉伸模態(tài)，特征振動溫度為3 329 K。因為彎曲模態(tài)對應(yīng)的能級最低，低溫時CO2的絕大多數(shù)振動能集中在這一模態(tài)。此外，這一模態(tài)會與拉伸模態(tài)共振，導(dǎo)致振動-振動能量的快速交換[14]。

TCE模型可以采用Evans等給出的火星大氣化學(xué)反應(yīng)數(shù)據(jù)[16]。不考慮離子組份和Ar，一共有8種組份：O2、N2、O、N、NO、C、CO、CO2。反應(yīng)共計54組，其中40組離解反應(yīng)，14組置換反應(yīng)，不考慮復(fù)合反應(yīng)。

QK模型不依賴于宏觀的化學(xué)反應(yīng)速率系數(shù)，具體計算方法分析如下。

1.1 離解反應(yīng)

離解反應(yīng)AB+M→A+B+M發(fā)生的條件為

imax>Θd/Θv

(1)

式中：AB代表分子，A和B代表原子，M可以是分子，也可以是原子;Θd是分子離解特征溫度;Θv是振動特征溫度;imax是碰撞能Ec對應(yīng)的最大振動能級，即

(2)

Ec=εtr,AB,M+εv,AB

(3)

1.2 置換反應(yīng)

置換反應(yīng)A+B→C+D發(fā)生的概率為

(4)

式中：A和C是分子；B和D是原子;Ea是活化能；ω是黏性溫度指數(shù)。當(dāng)Ea/k<Θv時，分母中的求和可近似等于1。

量子動理學(xué)火星化學(xué)反應(yīng)模型與地球化學(xué)反應(yīng)模型在計算軟件內(nèi)的執(zhí)行是比較類似的，只需要注意CO2具有多個振動能級即可。

2 火星探測器再入特性計算與分析

2.1 計算條件

為了增大減阻效果，火星探測器外形基本上都是采用70°半頂角的大鈍錐前體外形，最早在海盜號上取得了成功。錐越鈍，相同表面積的阻力越大，熱流越?。诲F越尖，更加穩(wěn)定[17]。后來的探測任務(wù)大都吸取了海盜號的成功經(jīng)驗，只是大小和后體形狀存在不同。

探路者號的外形是一個半頂角為70.19°的球形鈍頭大底后接一個46.63°的錐形后體，最大直徑2.65 m，見圖3，圖中Location表示不同的位置，x和y分別是各位置點的坐標(biāo)，s為各點的弧長，從頭部開始算起,Rn為頭部半徑，Rc為肩部倒圓半徑，Rb為探測器最大半徑，Ra為底部半徑，V∞表示來流速度。探路者號的外形及其氣動數(shù)據(jù)很大部分借鑒于海盜號。由于重心偏后，在自由分子流以及大部分過渡流區(qū)域，探路者號是靜不穩(wěn)定的[1，18]。

圖3 火星探路者號外形[18]Fig.3 Configuration of Mars Pathfinder[18]

來流氣體中，兩種氣體的體積分?jǐn)?shù)為XCO2=0.953 7，XN2=0.046 3。壁面條件為完全漫反射條件，不考慮催化。計算采用變徑硬球(Variable Hard Sphere，VHS)模型，平動能和內(nèi)能之間的松弛采用Larsen-Borgnakke模型,轉(zhuǎn)動和振動松弛數(shù)分別為5和50?；瘜W(xué)反應(yīng)采用TCE模型和QK模型計算，為了顯示化學(xué)反應(yīng)的影響，還進(jìn)行了理想氣體無化學(xué)反應(yīng)的對比計算。由于流動的對稱性，為節(jié)約計算資源，只計算一半的流場。計算網(wǎng)格包含了前體和尾渦。

計算條件如表1所示，分別對應(yīng)的高度為85 km、95 km和110 km。表中，ρ∞為來流密度，T∞為來流溫度，Tw為壁面溫度，Kn∞,HS為基于硬球模型的Knudsen數(shù)。

表1 計算條件Table 1 Computational conditions

2.2 結(jié)果分析

圖4給出了使用TCE模型以及QK模型計算的駐點線上溫度分布與理想氣體結(jié)果的比較。從圖中可以看出，探測器前方有耗散的激波結(jié)構(gòu)形成，經(jīng)過激波之后，流場溫度迅速升高。流場呈現(xiàn)非常明顯的熱力學(xué)非平衡現(xiàn)象，平動溫度Tt遠(yuǎn)高于轉(zhuǎn)動溫度Tr和振動溫度Tv，振動溫度的激發(fā)非常有限。QK模型計算的平動溫度和轉(zhuǎn)動溫度的峰值都要遠(yuǎn)低于理想氣體的結(jié)果，振動溫度由于相對較低，差異不明顯，這是與離解化學(xué)反應(yīng)的吸熱特性相吻合的。但是TCE模型的計算結(jié)果則不同，反而與理想氣體的計算結(jié)果幾乎一致，沒有反映出化學(xué)反應(yīng)的影響。對兩種化學(xué)反應(yīng)模型的平衡反應(yīng)速率進(jìn)行比較發(fā)現(xiàn)，文獻(xiàn)[11]給出的TCE模型的反應(yīng)速率系數(shù)比QK模型要低好幾個數(shù)量級。因而使用這套TCE模型參數(shù)計算出來的化學(xué)反應(yīng)影響與QK模型的結(jié)果相比幾乎可以忽略。很多化學(xué)反應(yīng)的實驗數(shù)據(jù)都是在較低溫度下測量得到的，而本算例的流場溫度高達(dá)10 000 K，因此直接使用這些數(shù)據(jù)是不準(zhǔn)確的。而且針對同一化學(xué)反應(yīng)可能有許多組TCE模型參數(shù)，彼此差異很大，QK模型可以作為驗證這些參數(shù)適用性的一種手段。

圖4 95 km高度真實氣體與理想氣體在駐點線上的溫度分布比較Fig.4 Comparison of temperature distributions along stagnation line of real air and ideal air at 95 km altitude

為了便于對比驗證，本文根據(jù)QK模型的理論反應(yīng)速率擬合了一組TCE模型參數(shù)，詳見附錄A和文獻(xiàn)[19]。以下的計算結(jié)果中，TCE模型相關(guān)結(jié)果均使用根據(jù)QK模型擬合得到的參數(shù)。圖5 給出了無化學(xué)反應(yīng)(NR)和考慮化學(xué)反應(yīng)后駐點線上的溫度分布對比。從圖中可以看出，化學(xué)反應(yīng)使得平動溫度、轉(zhuǎn)動溫度的峰值都大大下降，振動溫度因為激發(fā)程度不高，所以差別并不大。考慮化學(xué)反應(yīng)的總溫度在激波位置比無化學(xué)反應(yīng)的總溫度更高，是因為離解反應(yīng)產(chǎn)生了大量的原子組分，其內(nèi)自由度均為0，使得流場統(tǒng)計得到的平均轉(zhuǎn)動和振動自由度降低，而總溫度是平動溫度、轉(zhuǎn)動溫度和振動溫度與相應(yīng)的自由度加權(quán)平均的結(jié)果。使用根據(jù)QK模型理論結(jié)果擬合得到的TCE模型參數(shù)的計算結(jié)果與QK模型直接計算結(jié)果吻合良好。

圖5 95 km高度化學(xué)反應(yīng)對駐點線上溫度分布的影響Fig.5 Effects of chemical reaction on temperature distribution along stagnation line at 95 km altitude

圖6給出了95 km高度下氣動力系數(shù)隨攻角α的變化曲線，F(xiàn)MF為自由分子流計算結(jié)果，MN為修正牛頓公式計算結(jié)果。從圖中可以看出，隨著攻角的增加，軸向力系數(shù)CA和阻力系數(shù)CD逐漸減小;升力系數(shù)CL為負(fù)，且與阻力相比很小;法向力系數(shù)CN和升力系數(shù)都隨攻角增加而增大;俯仰力矩系數(shù)Cm隨攻角變化不大，基本位于0附近。

圖6 95 km高度氣動力系數(shù)隨攻角的變化Fig.6 Variations of aerodynamic coefficients with angle of attack at 95 km altitude

圖7給出了氣動力系數(shù)隨稀薄程度的變化。曲線分別代表本文不同稀薄程度的計算結(jié)果，對應(yīng)的符號(圓點、三角和方框)代表文獻(xiàn)[18]的結(jié)果，分別代表自由分子流和修正牛頓理論的結(jié)果。從圖中可以看出，隨著Knudsen數(shù)的增大，即流動稀薄程度增加，軸向力系數(shù)、法向力系數(shù)和阻力系數(shù)都隨之增大，升力系數(shù)則相反，俯仰力矩系數(shù)甚至發(fā)生變號。所有的計算結(jié)果均位于自由分子流和修正牛頓理論結(jié)果之間。

圖7 氣動力系數(shù)隨稀薄程度的變化Fig.7 Variations of aerodynamic coefficients with degree of rarefaction

隨著稀薄程度的增加，分子自由程逐漸增大，來流分子與探測器表面的碰撞比來流分子之間以及來流分子與從探測器表面碰撞反射回的分子碰撞頻繁得多，距離較遠(yuǎn)的分子也有可能撞擊到探測器表面，因而黏性增強(qiáng)導(dǎo)致摩阻對軸向力的貢獻(xiàn)增大[20]。無量綱阻力系數(shù)增大主要是隨高度的增加動壓迅速降低造成的[6]。

雖然化學(xué)反應(yīng)對氣動力的影響較小，但是其對氣動熱的影響不容忽略。95 km高度0°攻角狀態(tài)下，不考慮化學(xué)反應(yīng)的駐點熱流約為21.7 kW；考慮化學(xué)反應(yīng)后，則降至約19.1 kW，降幅達(dá)12.0%。其余攻角下的駐點熱流降幅也有14%左右。這說明吸熱離解反應(yīng)在流場中占據(jù)主導(dǎo)地位，若考慮壁面催化反應(yīng)，熱流降幅將更顯著。駐點熱流隨攻角的變化不是單調(diào)的，不考慮化學(xué)反應(yīng)時，15°攻角下的駐點熱流值最低，考慮化學(xué)反應(yīng)后，20°攻角下的駐點熱流最低，詳見表2。

表2 95 km高度化學(xué)反應(yīng)對駐點熱流的影響

圖8給出了探測器對稱面剖面上的熱流分布，從圖中可以看出，熱流計算結(jié)果展現(xiàn)出較大的統(tǒng)計漲落；考慮化學(xué)反應(yīng)后，迎風(fēng)面熱流普遍降低；當(dāng)攻角較小時，駐點熱流值最高；隨著攻角的增大，肩部的熱流逐漸升高并超過駐點值；背風(fēng)面熱流很低，幾乎為零。

圖8 95 km高度化學(xué)反應(yīng)對對稱面剖面上熱流分布的影響Fig.8 Effects of chemical reactions on heat flux distributions on symmetry plane of body at 95 km altitude

3 結(jié) 論

本文概述了火星大氣環(huán)境以及火星探測器在進(jìn)入、下降和著陸階段所面臨的氣動問題。TCE模型和QK模型都可以用于火星大氣化學(xué)反應(yīng)模擬，只是反應(yīng)的組份、數(shù)目及類型與地球大氣化學(xué)反應(yīng)相比略有區(qū)別?；鹦谴髿夥磻?yīng)數(shù)據(jù)較少且不完整，不同的反應(yīng)速率系數(shù)得到的計算結(jié)果差別較大。

QK模型擺脫了化學(xué)反應(yīng)宏觀速率系數(shù)的限制，計算結(jié)果與內(nèi)在物理機(jī)理相符合，且不包含自由參數(shù)，結(jié)果具有唯一性。根據(jù)QK模型平衡反應(yīng)速率的理論結(jié)果還可以擬合出對應(yīng)的TCE模型參數(shù)，經(jīng)驗證使用該參數(shù)的TCE模型計算結(jié)果與QK模型直接計算的結(jié)果吻合良好。QK模型可以作為驗證TCE模型參數(shù)適用性的一種手段?；瘜W(xué)反應(yīng)對火星再入高空氣動力的影響較小，但是對氣動熱的影響不容忽略，因此在化學(xué)反應(yīng)模型的選取應(yīng)用上應(yīng)力求能貼近真實飛行狀況。

綜上，QK模型的應(yīng)用具有很大的優(yōu)勢，值得進(jìn)一步的發(fā)展。

附錄A

表A1 根據(jù)QK模型擬合的火星大氣化學(xué)反應(yīng)速率系數(shù)
Table A1 Atmosphere chemical reaction rates coefficient of Mars fitted by results of QK model

No.ReactionΛηEa/J1CO+O2→C+O+O21.108×10-8-1.3331.780×10-182CO+N2→C+O+N21.232×10-8-1.3301.780×10-183CO+O→C+O+O9.666×10-9-1.3371.780×10-184CO+N→C+O+N1.010×10-8-1.3371.780×10-185CO+NO→C+O+NO1.122×10-8-1.3361.780×10-186CO+C→C+O+C9.829×10-9-1.3371.780×10-18

續(xù)表

No.ReactionΛηEa/J7CO+CO→C+O+CO6.295×10-9-1.3291.780×10-188CO+CO2→C+O+CO21.115×10-8-1.3511.780×10-189N2+O2→N+N+O28.233×10-9-1.3261.563×10-1810N2+N2→N+N+N24.555×10-9-1.3221.563×10-1811N2+O→N+N+O7.217×10-9-1.3301.563×10-1812N2+N→N+N+N7.537×10-9-1.3301.563×10-1813N2+NO→N+N+NO8.369×10-9-1.3291.563×10-1814N2+C→N+N+C7.337×10-9-1.3301.563×10-1815N2+CO→N+N+CO9.292×10-9-1.3211.563×10-1816N2+CO2→N+N+CO28.525×10-8-1.3461.563×10-1817NO+O2→N+O+O22.789×10-9-1.2921.040×10-1818NO+N2→N+O+N23.036×10-9-1.2871.040×10-1819NO+O→N+O+O2.506×10-9-1.2971.040×10-1820NO+N→N+O+N2.619×10-9-1.2971.040×10-1821NO+NO→N+O+NO1.435×10-9-1.2951.040×10-1822NO+C→N+O+C2.554×10-9-1.2971.040×10-1823NO+CO→N+O+CO3.078×10-9-1.2851.040×10-1824NO+CO2→N+O+CO23.156×10-8-1.3191.040×10-1825CO2+O2→CO+O+O22.290×10-10-1.1654.741×10-1926CO2+N2→CO+O+N22.416×10-10-1.1574.741×10-1927CO2+O→CO+O+O2.291×10-10-1.1734.741×10-1928CO2+N→CO+O+N2.408×10-10-1.1734.741×10-1929CO2+NO→CO+O+NO2.406×10-10-1.1714.741×10-1930CO2+C→CO+O+C2.393×10-10-1.1734.741×10-1931CO2+CO→CO+O+CO2.421×10-10-1.1554.741×10-1932CO2+CO2→CO+O+CO21.454×10-10-1.2084.741×10-1933O2+O2→O+O+O23.588×10-10-1.1938.197×10-1934O2+N2→O+O+N27.742×10-10-1.1878.197×10-1935O2+O→O+O+O6.517×10-10-1.1998.197×10-1936O2+N→O+O+N6.820×10-10-1.1998.197×10-1937O2+NO→O+O+NO7.441×10-10-1.1978.197×10-1938O2+C→O+O+C6.646×10-10-1.1998.197×10-1939O2+CO→O+O+CO7.824×10-10-1.1858.197×10-1940O2+CO2→O+O+CO28.596×10-10-1.2258.197×10-1941CO+N→NO+C9.590×10-14-0.6727.404×10-1942NO+C→CO+N1.085×10-13-0.713043CO+CO→CO2+C4.370×10-14-0.6379.072×10-1944CO2+C→CO+CO1.349×10-13-0.831045CO+O→O2+C9.183×10-14-0.6729.601×10-19

續(xù)表

No.ReactionΛηEa/J46O2+C→CO+O1.019×10-13-0.711047N2+O→NO+N9.382×10-14-0.6695.238×10-1948NO+N→N2+O1.113×10-13-0.713049NO+CO→CO2+N1.008×10-13-0.6781.666×10-1950CO2+N→NO+CO1.386×10-13-0.836051NO+O→O2+N1.064×10-13-0.7132.687×10-1952O2+N→NO+O9.663×10-14-0.716053CO2+O→O2+CO1.319×10-13-0.8363.823×10-1954O2+CO→CO2+O8.753×10-14-0.6813.300×10-19