探月返回器稀薄氣體熱化學(xué)非平衡特性數(shù)值模擬

梁 杰,李志輝,杜波強(qiáng),方 明

(中國(guó)空氣動(dòng)力研究與發(fā)展中心超高速所，綿陽(yáng)621000)

探月返回器稀薄氣體熱化學(xué)非平衡特性數(shù)值模擬

梁 杰,李志輝?,杜波強(qiáng),方 明

(中國(guó)空氣動(dòng)力研究與發(fā)展中心超高速所，綿陽(yáng)621000)

通過(guò)發(fā)展直角/非結(jié)構(gòu)混合網(wǎng)格下，適于復(fù)雜返回器含內(nèi)能激發(fā)的五組元熱化學(xué)非平衡稀薄氣體直接模擬蒙特卡洛(DSMC)方法，數(shù)值模擬了鈍體返回器稀薄流區(qū)7.5 km/s和10.6 km/s再入速度下的非平衡流場(chǎng)特性和對(duì)物面氣動(dòng)力/熱分布的影響。計(jì)算結(jié)果表明在月球返回速度下稀薄過(guò)渡流區(qū)存在著強(qiáng)烈的熱化學(xué)非平衡，近連續(xù)滑移區(qū)僅在激波層區(qū)域存在著較大的平動(dòng)、轉(zhuǎn)動(dòng)和振動(dòng)非平衡度，而在高稀薄流區(qū)，熱力學(xué)非平衡遍布繞流物體四周，包括整個(gè)壓縮區(qū)和尾跡區(qū)。高溫真實(shí)氣體效應(yīng)對(duì)表面熱環(huán)境影響顯著，對(duì)表面壓力和摩擦力的影響相對(duì)較弱。計(jì)算結(jié)果顯示探月返回器再入繞流過(guò)程100 km飛行高度仍需要考慮熱化學(xué)非平衡的影響。

探月返回器；稀薄流域；熱化學(xué)非平衡；DSMC方法；數(shù)值模擬

1 引言

月球探測(cè)返回艙是月球樣品的承載和運(yùn)輸工具，以近第二宇宙速度、半彈道跳躍式再入大氣層，其間經(jīng)歷自由分子流、稀薄過(guò)渡流和連續(xù)流等多個(gè)流動(dòng)區(qū)域，其再入飛行過(guò)程的靜動(dòng)態(tài)氣動(dòng)特性和熱環(huán)境都具有新的特點(diǎn)。當(dāng)返回艙以極高的速度再入大氣層時(shí)，由于激波和粘性的作用會(huì)產(chǎn)生嚴(yán)重的氣動(dòng)加熱，使返回器表面及其周?chē)鲌?chǎng)高溫氣體發(fā)生復(fù)雜的物理、化學(xué)變化，這個(gè)過(guò)程伴隨著氣體分子的離解和電離，同時(shí)也會(huì)出現(xiàn)分子的振動(dòng)激發(fā)、束縛電子的激發(fā)和輻射熱，壁面的化學(xué)反應(yīng)也可能發(fā)生[1-2]。在高空稀薄氣體流域，由于氣體的密度低，粒子之間的碰撞頻率降低，超高速返回器繞流流場(chǎng)中氣體分子的平動(dòng)溫度、轉(zhuǎn)動(dòng)溫度和振動(dòng)溫度處于高度的非平衡狀態(tài)，激波后高溫氣體發(fā)生化學(xué)反應(yīng)的弛豫時(shí)間遠(yuǎn)大于流動(dòng)的特征時(shí)間，因而在返回器的繞流中會(huì)出現(xiàn)嚴(yán)重的熱化學(xué)非平衡現(xiàn)象[3-4]。這些氣體流動(dòng)現(xiàn)象會(huì)造成流動(dòng)結(jié)構(gòu)的顯著變化，并對(duì)返回器的氣動(dòng)特性產(chǎn)生重要影響。詳細(xì)了解稀薄過(guò)渡流區(qū)的熱化學(xué)非平衡特性以及對(duì)返回器表面特性的影響對(duì)高超聲速返回器的設(shè)計(jì)非常重要，但地面試驗(yàn)設(shè)備很難模擬低密度環(huán)境下的高焓真實(shí)氣體效應(yīng)并且試驗(yàn)費(fèi)用昂貴。隨著高性能計(jì)算機(jī)的發(fā)展以及不斷更新的物理化學(xué)模型，使大規(guī)模并行計(jì)算成為研究稀薄過(guò)渡流區(qū)熱化學(xué)非平衡效應(yīng)的主要手段，基于氣體分子速度分布函數(shù)為Maxwell平衡態(tài)分布概念的連續(xù)流氣體動(dòng)力學(xué)計(jì)算理論，如求解NS方程，已經(jīng)不再適用于稀薄氣體流域的模擬，取而代之的是基于分子氣體動(dòng)力學(xué)的直接模擬蒙特卡洛(DSMC)方法[5]，通過(guò)模擬分子在計(jì)算區(qū)域內(nèi)運(yùn)動(dòng)、相互之間碰撞以及與固體表面之間的相互作用來(lái)實(shí)現(xiàn)真實(shí)氣體流動(dòng)的模擬。本文采用流場(chǎng)直角/表面非結(jié)構(gòu)混合網(wǎng)格結(jié)構(gòu)的DSMC方法(已在文獻(xiàn)[6，7]驗(yàn)證了該算法可靠性)，發(fā)展適于復(fù)雜返回器含內(nèi)能激發(fā)的五組元熱化學(xué)非平衡模擬技術(shù)，研究分析鈍體返回器在稀薄過(guò)渡流區(qū)不同再入速度下的非平衡流場(chǎng)特性和對(duì)物面氣動(dòng)力/熱參數(shù)分布的影響規(guī)律。

2 模擬熱化學(xué)非平衡流動(dòng)的DSMC方法

2.1 依賴(lài)于溫度的轉(zhuǎn)動(dòng)和振動(dòng)松弛模型

式中，ζt、ζr和ζv分別是分子平動(dòng)、轉(zhuǎn)動(dòng)和振動(dòng)自由度。直到二十世紀(jì)九十年代DSMC方法中的轉(zhuǎn)動(dòng)和振動(dòng)碰撞數(shù)都取為常數(shù)(Zr＝5，Zv＝50)[10]。根據(jù)分子運(yùn)動(dòng)論可知，同一密度條件下的氣體流動(dòng)，溫度越高，分子間的碰撞越頻繁，松弛的過(guò)程也就越快，因此松弛速率總是與溫度有關(guān)的函數(shù)。Parker[11]給出了連續(xù)流氣體分子轉(zhuǎn)動(dòng)松弛所需要的碰撞數(shù)ZCr與氣體溫度之間的關(guān)系表達(dá)式如式(2):

式中，T?是氣體分子作用勢(shì)的特征溫度，(Zr)∞是實(shí)驗(yàn)測(cè)定的極限值。

連續(xù)流的振動(dòng)松弛碰撞數(shù)定義為式(3):

式中，τc為分子平均碰撞時(shí)間，τMW為振動(dòng)松馳時(shí)間，τP為高溫修正項(xiàng)。

由Landau-Teller理論結(jié)果給出式(4):

其中，p＝nkT。

Millikan和White[12]在高溫激波管內(nèi)通過(guò)干涉儀觀測(cè)了氣體分子振動(dòng)松馳過(guò)程，在104K的溫度范圍內(nèi)，給出了如式(5)所示擬合參數(shù):

其中，mr為分子折合質(zhì)量，θv是振動(dòng)特征溫度。Park[13]給出了高溫修正的τP為式(6):

式中，n是數(shù)密度，σv是分子振動(dòng)碰撞截面(＝10－20m2)，是分子平均熱運(yùn)動(dòng)速度。

2.2 化學(xué)反應(yīng)模型及非彈性碰撞模擬策略

通常反應(yīng)速率常數(shù)可以寫(xiě)成式(7)所示的Arrhenius方程形式:

式中，Ea是反應(yīng)中需要的活化能，k是Boltzmann常數(shù)，Λ和η是常數(shù)，由試驗(yàn)定出。

化學(xué)反應(yīng)幾率可以推導(dǎo)為式(8)所示表達(dá)式[5]:

式中，ε是對(duì)稱(chēng)因子，對(duì)于不同類(lèi)分子ε＝1，對(duì)同類(lèi)分子ε＝2。為碰撞分子的平均內(nèi)自由度，ωAB為碰撞對(duì)的粘性溫度指數(shù)。當(dāng)碰撞中的總碰撞能量Ec＞Ea時(shí)，反應(yīng)截面σR與總碰撞截面σT的比值就是反應(yīng)發(fā)生的概率，這種化學(xué)反應(yīng)模型也稱(chēng)為T(mén)CE(Total Collision Energy)模型。

在DSMC模擬過(guò)程中為簡(jiǎn)化碰撞算法，假設(shè)非彈性碰撞的化學(xué)反應(yīng)過(guò)程和內(nèi)能交換是相互獨(dú)立的事件。由于碰撞分子發(fā)生化學(xué)反應(yīng)的幾率在所有非彈性碰撞事件中是最低的，所以首先判斷其是否發(fā)生化學(xué)反應(yīng)。離解反應(yīng)、交換反應(yīng)等不同的化學(xué)反應(yīng)過(guò)程在滿(mǎn)足反應(yīng)幾率的條件下都有可能發(fā)生。在無(wú)法獲知反應(yīng)分子遵循何種反應(yīng)軌跡生成新的反應(yīng)生成物時(shí)，可以假定每種反應(yīng)過(guò)程都是獨(dú)立的。對(duì)所有可能發(fā)生化學(xué)反應(yīng)的幾率求和得到總的反應(yīng)幾率，如果Pr＞ Rf，碰撞分子發(fā)生化學(xué)反應(yīng)[14]。具體發(fā)生的是何種反應(yīng)類(lèi)型，則同樣根據(jù)幾率Pir/Pr與隨機(jī)數(shù)Rf的比較來(lái)選取，然后按照相應(yīng)的化學(xué)反應(yīng)類(lèi)型產(chǎn)生新的生成物(分子或原子)。

對(duì)于不發(fā)生化學(xué)反應(yīng)的碰撞分子，即Rf＞Pr時(shí)，則根據(jù)分子振動(dòng)抽樣概率，判斷是否發(fā)生振動(dòng)自由度的激發(fā)；對(duì)于未發(fā)生振動(dòng)激發(fā)的分子，再根據(jù)分子轉(zhuǎn)動(dòng)抽樣概率，判斷是否發(fā)生轉(zhuǎn)動(dòng)自由度的激發(fā)；對(duì)于發(fā)生內(nèi)自由度激發(fā)的分子，按相應(yīng)的能量交換模型進(jìn)行能量的再分配，否則分子按彈性碰撞處理。

本文模擬了五組元(O2，N2，O，N，NO)氣體的離解反應(yīng)和交換反應(yīng)，未考慮電離的發(fā)生。

2.3 不同溫度模式的非平衡度定義

為研究航天器再入高超聲速繞流稀薄氣體環(huán)境下，分子內(nèi)能激發(fā)致平動(dòng)、轉(zhuǎn)動(dòng)和振動(dòng)非平衡的程度，本文定義變量η作為不同溫度模式非平衡程度的度量參數(shù)。對(duì)于平動(dòng)溫度的非平衡度用三個(gè)分量來(lái)表示，其定義為式(9):

其中，Tt是分子的平動(dòng)溫度，(Tt)i，i＝x，y，z分別是三個(gè)方向的速度分量定義的平動(dòng)溫度。分別定義轉(zhuǎn)動(dòng)溫度非平衡度和振動(dòng)溫度非平衡度為式(10)、(11):

其中，Tr和Tv分別是轉(zhuǎn)動(dòng)溫度、振動(dòng)溫度。

3 計(jì)算結(jié)果分析

本文選取類(lèi)飛船返回艙的探月返回器外形作為研究對(duì)象，來(lái)流速度分別為 7.5 km/s和10.6 km/s，計(jì)算的高度范圍為80～120 km(Kn數(shù)為0.0025～2.44)，攻角為－20°，壁面溫度設(shè)置為500 K，采用完全漫反射模型，因不考慮側(cè)滑角，僅計(jì)算對(duì)稱(chēng)的半個(gè)流場(chǎng)，計(jì)算中考慮了分子的轉(zhuǎn)動(dòng)、振動(dòng)能量激發(fā)以及五組元空氣的離解反應(yīng)和交換反應(yīng)。

圖1、圖2分別是來(lái)流速度7.5 km/s條件下80 km和100 km高度返回器繞流流場(chǎng)的平動(dòng)溫度、轉(zhuǎn)動(dòng)溫度和振動(dòng)溫度沿中心流線的變化曲線。遠(yuǎn)前方氣流在返回器的前方被激波壓縮加熱，吸收的能量首先供應(yīng)給平動(dòng)能，使平動(dòng)溫度迅速上升；隨著分子碰撞松弛影響的進(jìn)一步加劇，氣體分子碰撞致內(nèi)部能量迅速增加，分子轉(zhuǎn)動(dòng)和振動(dòng)能態(tài)也相繼激發(fā)。受壁面溫度條件限制，最后平動(dòng)、轉(zhuǎn)動(dòng)和振動(dòng)溫度在物面附近被強(qiáng)制平衡到壁面溫度。由此可見(jiàn)，在激波和激波壓縮層區(qū)域，平動(dòng)溫度、轉(zhuǎn)動(dòng)溫度和振動(dòng)溫度之間較大程度的偏離，反映出氣體流動(dòng)呈現(xiàn)出強(qiáng)烈的熱力學(xué)非平衡特征，這也是稀薄氣體流域所特有的流動(dòng)現(xiàn)象。隨著氣體稀薄程度的增大，來(lái)流氣體密度降低，分子之間的碰撞頻率減弱，分子轉(zhuǎn)動(dòng)和振動(dòng)自由度的激發(fā)程度明顯下降，引起轉(zhuǎn)動(dòng)溫度和振動(dòng)溫度的迅速降低，導(dǎo)致平動(dòng)溫度與轉(zhuǎn)動(dòng)溫度、振動(dòng)溫度之間差異增大、非平衡度增加。

圖1 返回器再入80 km流場(chǎng)沿中心流線溫度分布Fig.1 Temperature distribution around reentry capsule along central streamline at 80 km height

圖2 返回器再入100 km流場(chǎng)沿中心流線溫度分布Fig.2 Temperature distribution around reentry capsule along central streamline at 100 km height

從前面定義的不同溫度模式的非平衡度，可以了解整個(gè)流場(chǎng)的熱力學(xué)非平衡分布情況。圖3、圖4分別給出了來(lái)流速度7.5 km/s返回器再入80 km和100 km高度繞流流場(chǎng)的轉(zhuǎn)動(dòng)溫度和振動(dòng)溫度非平衡度等值線分布。80 km的Kn數(shù)為0.0025，流動(dòng)屬于近連續(xù)滑移流區(qū)，圖中顯示出頭部激波的壓縮性非常強(qiáng)，激波層很薄，在流場(chǎng)的大部分區(qū)域都處于平衡態(tài)，熱力學(xué)非平衡主要分布在激波層附近區(qū)域。返回器鈍體前方激波壓縮溫度迅速上升引起嚴(yán)重的轉(zhuǎn)動(dòng)和振動(dòng)非平衡，在后體激波層區(qū)域振動(dòng)非平衡度仍然較高，并且分布的區(qū)域明顯寬于轉(zhuǎn)動(dòng)非平衡度，主要是由于振動(dòng)能激發(fā)相對(duì)較弱且滯后于轉(zhuǎn)動(dòng)能，振動(dòng)溫度與平動(dòng)溫度的偏離較大造成的。在100 km高度時(shí)已看不到明顯的激波，而在返回器前部出現(xiàn)較強(qiáng)的氣體壓縮強(qiáng)擾動(dòng)變化過(guò)程，邊界層與強(qiáng)擾動(dòng)激波壓縮層相互融合，出現(xiàn)厚厚的強(qiáng)擾動(dòng)激波過(guò)渡層，并且是隨著繞流物體后部稀薄程度增加，激波層厚度也隨之增大，變成了一個(gè)強(qiáng)擾動(dòng)壓縮過(guò)渡帶。圖4顯示出在整個(gè)壓縮區(qū)和尾跡區(qū)均存在較嚴(yán)重的轉(zhuǎn)動(dòng)和振動(dòng)非平衡。同樣，在80 km高度激波層區(qū)域由不同方向速度分量定義的平動(dòng)溫度本身也存在著較大的非平衡。

圖3 返回器再入80 km流場(chǎng)轉(zhuǎn)動(dòng)（左）和振動(dòng)（右）溫度非平衡度等值線分布Fig.3 Nonequilibrium degree contours of rotational（left）and vibrational（right）temperature at 80 km height

圖4 返回器再入100 km流場(chǎng)轉(zhuǎn)動(dòng)（左）和振動(dòng)（右）溫度非平衡度等值線分布Fig.4 Nonequilibrium degree contours of rotational（left）and vibrational（right）temperature at 100 km height

圖5 是x方向和y方向的平動(dòng)溫度非平衡度分布，在頭部脫體激波壓縮最強(qiáng)的地方也是x方向速度變化最劇烈的位置，x方向的平動(dòng)溫度非平衡最強(qiáng)，而在沿返回器身部的脫體激波層區(qū)域，則是y方向的平動(dòng)溫度非平衡度最大。z方向的平動(dòng)溫度非平衡度在激波層區(qū)域從0.5～0.7之間變化。同時(shí)也看到，在后體尾跡膨脹區(qū)也有局部的平動(dòng)非平衡。隨著飛行高度的增加，平動(dòng)溫度非平衡區(qū)域不斷擴(kuò)大，到返回器再入100 km高度繞流流場(chǎng)，在整個(gè)壓縮區(qū)域和尾跡區(qū)也都存在著較嚴(yán)重的平動(dòng)溫度非平衡。

圖5 返回器再入80 km流場(chǎng)平動(dòng)溫度非平衡度等值線分布Fig.5 Nonequilibrium degree contours of translational temperature around reentry capsule at 80 km height

返回器在以超高速再入飛行時(shí)，雙原子分子之間較高的碰撞能量不但會(huì)引起轉(zhuǎn)動(dòng)和振動(dòng)自由度的激發(fā)，形成較強(qiáng)的熱力學(xué)非平衡流動(dòng)，還會(huì)發(fā)生一定程度的化學(xué)反應(yīng)，生成新的化學(xué)組元。圖6給出了返回器再入80 km高度、來(lái)流速度分別為7.5 km/s和10.6 km/s繞流流場(chǎng)中化學(xué)組元摩爾分?jǐn)?shù)沿中心流線的變化情況，圖中顯示出返回器繞流流場(chǎng)存在較強(qiáng)的化學(xué)非平衡流動(dòng)特征。在激波層區(qū)域，氣體分子的高密度產(chǎn)生高碰撞率，使離解活化能較低的O2分子離解速率遠(yuǎn)大于N2的離解速率，在氣流到達(dá)物面前，O2全部離解完畢。

N2分子在10.6 km/s的來(lái)流速度下大部分離解，使N原子的摩爾濃度接近了N2分子發(fā)生離解前的初始摩爾濃度，如圖7分別展示了N2分子和N原子摩爾分?jǐn)?shù)流場(chǎng)等值線分布，激波后的返回艙周?chē)植贾罅康碾x解后的N原子，顯示出強(qiáng)烈的化學(xué)非平衡效應(yīng)。

圖6 返回器再入80 km流場(chǎng)沿中心流線組元摩爾分?jǐn)?shù)分布Fig.6 Mole fraction distribution of gas species around reentry capsule along central streamline at 80 km height

圖7 返回器再入80 km和10.6 km/s速度下流場(chǎng)N2分子和N原子摩爾分?jǐn)?shù)分布Fig.7 Mole fraction distribution of N2molecule and N atom around reentry capsule at 80 km height with 10.6 km/s velocity

對(duì)比圖6(a)中N2分子在7.5 km/s來(lái)流速度下離解反應(yīng)要弱得多，表明同一飛行高度，返回器以近第二宇宙速度10.6 km/s再入比第一宇宙速度7.5 km/s再入致流場(chǎng)呈現(xiàn)更加嚴(yán)重的熱化學(xué)非平衡流動(dòng)特征。在五個(gè)化學(xué)組元中NO的含量是最低的，并且是先上升再下降接近于零，說(shuō)明發(fā)生交換反應(yīng)生成NO的幾率并不高，同時(shí)也表明生成的NO在靠近物面前也基本離解完畢。圖8繪出返回器以10.6 km/s再入100 km高度飛行繞流流場(chǎng)沿中心流線組元摩爾分?jǐn)?shù)分布曲線，與圖6(b)對(duì)比看出，隨著飛行高度由80 km增加到100 km，稀薄效應(yīng)變得很?chē)?yán)重，致使化學(xué)反應(yīng)的強(qiáng)度不斷減弱，但在 10.6 km/s的來(lái)流速度下，100 km高度仍有較強(qiáng)的化學(xué)反應(yīng)發(fā)生。這與通常情況，返回器以第一宇宙速度7.5 km/s飛行繞流在100 km以上幾乎不再有化學(xué)反應(yīng)發(fā)生，形成鮮明的對(duì)比，預(yù)示未來(lái)我國(guó)發(fā)展的探月返回器再入繞流過(guò)程100 km以上就要考慮化學(xué)反應(yīng)的發(fā)生為宜。因此不同的來(lái)流速度和來(lái)流密度條件下，發(fā)生離解反應(yīng)和交換反應(yīng)的程度以及生成新組元的分布規(guī)律有明顯不同，對(duì)比圖6、圖8各組元在流場(chǎng)中的分布，同樣發(fā)現(xiàn)來(lái)流初始條件對(duì)流場(chǎng)的化學(xué)非平衡特性以及組元分布有較大影響。

圖8 返回器再入100 km流場(chǎng)沿中心流線組元摩爾分?jǐn)?shù)分布（V＝10.6 km/s）Fig.8 Mole fraction distribution of gas species along central streamline at 100 km height（V＝10.6 km/s）

圖9～11展示了不同高度、來(lái)流速度7.5 km/s時(shí)高溫氣體熱化學(xué)非平衡效應(yīng)對(duì)返回器子午面φ＝0°，180°沿x－方向不同位置氣動(dòng)力/熱參數(shù)分布的影響，圖中實(shí)線是計(jì)算中僅考慮N2和O2混合氣體分子的轉(zhuǎn)動(dòng)能激發(fā)，定義為T(mén)R氣體。而反應(yīng)氣體則是考慮了N2和O2混合氣體分子的轉(zhuǎn)動(dòng)能和振動(dòng)能激發(fā)以及化學(xué)反應(yīng)。圖9(a)顯示出在80 km高度時(shí)TR氣體計(jì)算的熱流明顯高于反應(yīng)氣體的熱流，在駐點(diǎn)區(qū)域大一倍左右，說(shuō)明高溫氣體熱化學(xué)非平衡效應(yīng)對(duì)返回器表面熱環(huán)境有顯著的影響，主要是考慮熱化學(xué)反應(yīng)真實(shí)氣體分子能量有很大一部分因激發(fā)分子振動(dòng)自由度與發(fā)生離解反應(yīng)吸熱，導(dǎo)致氣體溫度降低，物面溫度梯度與熱流下降。隨著飛行高度的增加，無(wú)論是分子的振動(dòng)激發(fā)還是化學(xué)反應(yīng)都不斷減弱，到100 km高度時(shí)兩種氣體計(jì)算結(jié)果表明，表面熱流幾乎是重合的(圖9(b))。另外，熱流分布也表明在－20°攻角下，最大熱流位于肩部而不是駐點(diǎn)?；瘜W(xué)反應(yīng)對(duì)表面壓力和摩檫力系數(shù)分布的影響，相比熱流的影響要弱得多?；瘜W(xué)反應(yīng)使得氣體壓縮和氣體膨脹更劇烈，在前體脫體激波層化學(xué)反應(yīng)致氣體壓縮性更強(qiáng)，波后壓力更大，因此反應(yīng)氣體的頭部激波更靠近物體，頭部區(qū)域的壓力要高于TR氣體，但在身部膨脹區(qū)，熱化學(xué)非平衡效應(yīng)加劇了膨脹，使壓力減小，以致TR氣體計(jì)算的壓力要稍高些，見(jiàn)圖10。對(duì)于摩檫力系數(shù)來(lái)說(shuō)，無(wú)論頭部還是身部都是TR氣體計(jì)算值高于反應(yīng)氣體，見(jiàn)圖11，主要是熱化學(xué)反應(yīng)真實(shí)氣體分子能量因振動(dòng)能激發(fā)與離解吸熱而致能量降低，跨越脫體激波層氣流速度較僅考慮轉(zhuǎn)動(dòng)能激發(fā)的TR氣體波后速度更小，以致物面切向速度變化梯度低于TR氣體所致，這種情況在φ＝0°的迎風(fēng)子午面因高超聲速再入脫體激波明晰強(qiáng)烈表現(xiàn)更為明顯，而在φ＝180°背風(fēng)子午面因脫體激波較弱使得考慮與不考慮熱化學(xué)反應(yīng)真實(shí)氣體效應(yīng)，兩種情況計(jì)算值差別不大，圖11中摩擦力系數(shù)在φ＝0°、180°兩個(gè)子午面的變化曲線顯示出熱化學(xué)非平衡效應(yīng)對(duì)返回器表面摩檫力系數(shù)造成影響這一變化規(guī)律。因此，頭部壓力增大對(duì)俯仰力矩的貢獻(xiàn)與身部摩擦力減小的貢獻(xiàn)相抵消，造成高溫氣體熱化學(xué)非平衡效應(yīng)對(duì)這種返回器配平攻角的影響較小。為了定量化分析上述兩種氣體模型計(jì)算得到的俯仰力矩系數(shù)變化特點(diǎn)，圖12藍(lán)色虛點(diǎn)線與藍(lán)色實(shí)線繪出兩種氣體模型計(jì)算的該返回器以攻角－20°再入120～80 km高度所受的俯仰力矩系數(shù)比較，看出兩者相差微弱，表明高溫?zé)峄瘜W(xué)非平衡效應(yīng)對(duì)返回器俯仰力矩系數(shù)影響甚微。

圖12、圖13繪出了返回器以來(lái)流速度7.5 km/s攻角－20°再入120～80 km高度的氣動(dòng)力系數(shù)隨高度變化規(guī)律。隨著飛行高度的增加，阻力系數(shù)是不斷增大的，升力系數(shù)和升阻比則迅速下降，壓心前移，引起俯仰力矩系數(shù)(負(fù)值為低頭力矩)在數(shù)值上不斷增大。如果要達(dá)到配平狀態(tài)，則需要更大的負(fù)向攻角產(chǎn)生正的抬頭力矩才能實(shí)現(xiàn)，也就是說(shuō)隨著飛行高度增加，配平攻角是增大的。在高空稀薄流域，壓力對(duì)氣動(dòng)力系數(shù)的貢獻(xiàn)隨高度的變化比較小，但摩擦力的貢獻(xiàn)則是隨高度的增加而增大，這也是引起配平攻角隨飛行高度有較大變化的主要原因。

圖9 化學(xué)反應(yīng)對(duì)返回器0°和180°子午面熱流影響比較Fig.9 Effect of chemical reactions on surface heat flux around reentry capsule along φ ＝0°，180°

圖10 化學(xué)反應(yīng)對(duì)返回器0°和180°子午面壓力影響比較Fig.10 Effect of chemical reactions on surface pressure around reentry capsule along φ＝0°，180°

圖11 化學(xué)反應(yīng)對(duì)返回器0°和180°子午面摩檫力影響比較Fig.11 Effect of chemical reactions on surface friction around reentry capsule along φ＝0°，180°

圖12 俯仰力矩和壓心系數(shù)隨高度變化Fig.12 Pitching moment and center of pressure coefficient vs flight altitude

圖13 升阻特性隨高度變化Fig.13 Lift-drag feature vs flight altitude

4 結(jié)論

本文給出了雙原子分子轉(zhuǎn)動(dòng)和振動(dòng)松弛碰撞數(shù)計(jì)算公式、TCE化學(xué)反應(yīng)模型，建立了模擬分子非彈性碰撞先后次序的算法，采用基于流場(chǎng)直角與物面非結(jié)構(gòu)混合網(wǎng)格的DSMC方法模擬了鈍體返回器不同再入速度、不同飛行高度下的熱化學(xué)非平衡流動(dòng)現(xiàn)象和氣動(dòng)力/熱變化特性，得出以下幾點(diǎn)結(jié)論:

1)在近連續(xù)滑移區(qū)，流場(chǎng)大部分區(qū)域都處于平衡態(tài)，只有在較窄的激波層區(qū)域，平動(dòng)溫度、轉(zhuǎn)動(dòng)溫度和振動(dòng)溫度之間存在著較嚴(yán)重的非平衡；隨著稀薄度的不斷增大，激波層和邊界層逐漸融合，熱力學(xué)非平衡區(qū)域擴(kuò)大至整個(gè)壓縮區(qū)和尾跡區(qū)；

2)在返回器再入80 km高度，氣體分子發(fā)生劇烈的化學(xué)反應(yīng)，流場(chǎng)呈現(xiàn)化學(xué)非平衡狀態(tài)。隨高度增加，化學(xué)反應(yīng)逐漸減弱，但探月返回器若以近第二宇宙速度從月球再入返回，在離地面100 km高空的稀薄條件下也會(huì)產(chǎn)生強(qiáng)烈的離解反應(yīng)。

3)高溫真實(shí)氣體效應(yīng)對(duì)低高度的表面熱環(huán)境影響較大，但對(duì)表面壓力和摩擦力分布影響相對(duì)較小，隨著飛行高度增加，影響呈不斷減弱的趨勢(shì)。

4)稀薄氣體效應(yīng)對(duì)探月返回器的氣動(dòng)特性影響顯著，使升阻比降低，壓心前移，配平攻角增大。

（

）

[1]Boyd I D.Modeling of associative ionization reactions in hypersonic rarefied flows[J].Physics of Fluids，2007，19(9): 3-14.

[2]Ozawa T，Levin D A，Nompelis I，et al.Particle and continuum method comparison of a high-altitude，extreme-Machnumber reentry flow[J].Journal of Thermophysics and Heat Transfer，2010，24(2):225-240.

[3]Lofthouse A J，Scalabriny L C，Boyd I D.Hypersonic aerothermodynamics analysis across nonequilibrium regimes using continuum and particle methods[R].AIAA 2007-3903，2007.

[4]Wen C Y，Massimi H S，Chen Y S，et al.Numerical simulations of nonequilibrium flows over rounded models at reentry speeds[R].AIAA 2012-5906，2012.

[5]Bird G A.Molecular gas dynamics and the direct simulation of gas flows[M].Oxford:Clarendon Press，1994:123-256.

[6]梁杰，閻超，杜波強(qiáng).基于兩級(jí)直角網(wǎng)格結(jié)構(gòu)的三維DSMC算法研究[J].空氣動(dòng)力學(xué)學(xué)報(bào)，2010，28(4):466-471.Liang Jie，Yan Chao，Du Boqiang.An algorithm study of three-dimensional DSMC simulation based on two-level Cartesian coordinates grid structure[J].Acta Aerodynamica Sinica，2010，28(4):466-471.(in Chinese)

[7]梁杰，閻超，楊彥廣，等.過(guò)渡區(qū)側(cè)向噴流干擾的并行DSMC數(shù)值模擬研究[J].宇航學(xué)報(bào)，2011，32(5):1012-1018.Liang Jie，Yan Chao，Yang Yanguang，et al.Parallel DSMC simulation of lateral jet interaction in rarefied transitional region[J].Journal of Astronautics，2011，32(5):1012-1018.(in Chinese)

[8]沈青.稀薄氣體動(dòng)力學(xué)[M].北京:國(guó)防工業(yè)出版社，2003:246-247.Shen Qing.Rarefied Gas Dynamics[M].Beijing:National Defense Industry Press，2003:246-247.(in Chinese)

[9]Lumpkin F E，Hass B L，Boyd I D.Resolution of differences between collision number definition in particle and continuum simulations[J].Physics of Fluids A:Fluid Dynamics(1989-1993)，1991，3(9):2282-2284.

[10]Boyd I D.Rotational and vibrational nonequilibrium effects in rarefied hypersonic flow[J].Journal of Thermophysics and Heat Transfer，1990，4(4):478-484.

[11]Parker J G.Rotational and vibrational relaxation in diatomic gases[J].Physics of Fluids，1959，2(4):449-462.

[12]Millikan R C，White D R.Systematics of vibrational relaxation[J].The Journal of Chemical Physics，1963，39(12): 3209-3213.

[13]Park C.Problems of rate chemistry in the flight regimes of aeroassisted orbital transfer vehicles[C]//Progress in Astronautics and Aeronautics，AIAA，New York，1985，96:511-537.

[14]Gimelshein N E，Gimelshein S F，Levin D A，et al.Reconsideration of DSMC models for internal energy transfer and chemical reactions[C]//23rd International Symposium of Rarefied Gas Dynamics，American Institute of Physics，2003: 349-357.

Numerical Simulation of Rarefied Gas Thermochemical Nonequilibrium when Lunar Exploration Vehicle Re-entering into Atmosphere

LIANG Jie，LI Zhihui?，DU Boqiang，F(xiàn)ANG Ming
(Hypervelocity Aerodynamics Institute，China Aerodynamics Research and Development Center，Mianyang 621000，China)

The Direct Simulation Monte Carlo(DSMC)method was developed to simulate rarefied gas flows around complex bodies with internal energy excitation and thermochemical nonequilibrium effect of five air species，using a hybrid grid technique of Cartesian coordinate meshes and surface unstructured triangular cells.The nonequilibrium flowfield characteristics and their effects on surface parameters distribution around blunt re-entry capsule in rarefied flow regime were numerically simulated under the re-entering velocities of 7.5km/s and 10.6km/s，respectively.The computational results showed that intensive thermochemical nonequilibrium occurred in the rarefied transitional region under lunar return velocity.In near continuum slip flow regime，the great degrees of translational，rotational and vibrational nonequilibrium only exited in the shock wave layer.While thermodynamic nonequilibrium distributed in the whole compress and wake regions at high rarefied flow regime.The high temperature real gas effect affected the surface thermal environment significantly and impacted relatively weak on the surface pressure and friction.It is indicated that the effect of thermochemical nonequilibrium must be considered up to 100km during the re-entering into the atmosphere of the lunar exploration vehicle.

lunar exploration reentry vehicle；rarefied flow region；thermochemical nonequilibrium；DSMC method；numerical simulation

V211.3

A

1674-5825(2015)03-0295-08

2014-09-11；

2015-03-15

國(guó)家自然科學(xué)基金(91016027、11325212)；973計(jì)劃(2014CB744100)；國(guó)防基礎(chǔ)科研基金(51313030104)

梁杰(1966－)，男，博士，研究員，研究方向?yàn)橄”怏w動(dòng)力學(xué)的數(shù)值方法及應(yīng)用。E-mail:liangjie29501＠163.com

李志輝(1968－)，男，博士，研究員，研究方向?yàn)榭缌饔蚩諝鈩?dòng)力學(xué)。E-mail:zhli0097＠x263.net