謝 東,王學德

(南京理工大學 能源與動力工程學院,江蘇 南京 210094)

作為地球軌道外的第一近鄰,火星的地質條件等與地球較為接近[1],是人類地外行星探測活動中的重要對象。探測器進入火星、減速過程穿越火星稀薄大氣的歷時很短,但由于探測器的速度極高,需要經受嚴酷的氣動加熱過程,因此該過程是探測任務中最危險也是最重要的環(huán)節(jié)[2]。

火星的大氣成分與地球顯著不同,主要由95.7%的CO2、2.7%的N2和其它稀有氣體組成[3],其地面平均大氣密度、大氣壓力約為地球的1%。美國早期的火星探測器飛行任務模擬主要采用以CO2大氣環(huán)境修正的地面實驗方法,但是其費用高昂、試驗難度大。隨著CFD技術的快速發(fā)展,美國從探路者號探測器(MPF)之后開展的氣動評估,幾乎全部基于計算流體力學和數(shù)值模擬[4]。目前國內主要采用求解連續(xù)Navier-Stokes方程和模擬稀薄氣體流動的直接模擬蒙塔卡洛方法(DSMC)進行火星大氣高溫真實氣體效應的研究,其中化學反應模型大多采用5組分17反應[5]或者8組分54反應[6]。

探測器進入火星大氣層需要經歷以CO2氣體為主導的高超聲速化學反應流動過程。由于CO2具有多個振動模態(tài),熱力學模型比地球大氣更復雜,且CO2相較于其它雙原子分子更容易離解,并在高溫下發(fā)生置換、電離等復雜的化學反應,這使得火星大氣環(huán)境下的熱化學非平衡效應比地球再入過程更為嚴重。

因此本文旨在探究CO2分子在高溫下的振動激發(fā)模式、化學反應模型以及在DSMC框架下的實現(xiàn)方式。構建了針對CO2氣體環(huán)境下5組分17反應的熱化學非平衡流動模型,編制了計算程序,并進行了驗證,證明了該物理化學模型與計算方法的可行性。

1 CO2的振動激發(fā)

1.1 CO2分子的振動模型

CO2分子為直線型三原子分子,常溫下具有3個平動自由度和2個轉動自由度。CO2分子具有4種振動模態(tài),其激發(fā)程度由溫度決定。如表1所示,其中第1、2模態(tài)為退化的彎曲模態(tài),特征溫度均為960 K;第3模態(tài)為對稱拉伸模態(tài),特征溫度為1 920 K;第4模態(tài)為反對稱拉伸模態(tài),特征溫度為3 380 K;CO2的離解特征溫度為66 000 K。低溫下CO2的大部分振動能都儲存在彎曲模態(tài);在高溫下,振動模式會被部分或完全激發(fā),所以CO2分子的實際振動模式是這幾種模態(tài)的組合。

表1 CO2分子的振動模態(tài)Table 1 Vibration modes of CO2molecules

1.2 CO2分子的振動能級與振動能

經典模型下,BORGNAKE和LARSEN[7]提出的L-B方法可以通過給每個分子分配連續(xù)振動能級,進而求得有效振動溫度。然而此方法涉及基于碰撞總能的隱式溫度,需要迭代求解。為此在DSMC方法中一般根據L-B模型給振動模態(tài)m分配離散的振動能級im:

im=?-ln(fR)T/Θm」

(1)

式中:?」表示取整,fR為(0,1)的隨機數(shù),T為分子熱力學溫度,Θm為對應振動模態(tài)的振動特征溫度,對應的振動能為

εm=imkΘm

(2)

式中:k為玻爾茲曼常數(shù)。

由于實際的熱化學反應流場溫度可達到10 000 K,分子的振動激發(fā)十分激烈,分子的原子間距增大,因此本文考慮了振動能的非諧性效應,從MORSE勢推導得到非諧振子模型[8]下各振動能級對應的特征溫度:

(3)

式中:Θd為離解特征溫度,η為非諧性因子,取值范圍為[0,1]。當η=0時,上式對應簡單諧振子模型。本文取η=0.90,則各模態(tài)的最高振動能級為

(4)

式中:Edis為離解能閾值。

圖1為當η=0.90,離解溫度為Θd/4時,CO2分子各模態(tài)的振動能級與相應能級的振動特征溫度。由于CO2的光譜圖并未顯示出很高的振動能級,本文假設各振動模態(tài)相互獨立,結合內能按自由度均分原則[9],將各模態(tài)振動能上限均設為Edis/4,也即各模態(tài)的離解特征溫度為Θd/4。這一方法既保證了振動能級不會太高而超出實際,也考慮到了非諧性效應,各能級的能量間隔隨著特征溫度升高而變小。

圖1 振動能級與對應的振動溫度Fig.1 Vibration energy level and corresponding vibration temperature

1.3 CO2分子的振動弛豫

分子間的碰撞分為彈性碰撞和非彈性碰撞,非彈性碰撞會導致總能量在分子的平動模態(tài)和內部模態(tài)之間重新分配。在碰撞過程中一般會引入振動松弛因子Zv和轉動松弛因子ZR來調整彈性碰撞和非彈性碰撞的比例,實現(xiàn)控制內能變化的速率,確保各參數(shù)不會失真的目的,各能量模式松弛概率滿足:

(5)

式中:v為碰撞頻率,τv為振動松弛時間。

CO2分子的轉動松弛因子可取常數(shù)ZR=6.5,其它雙原子分子為5[8]。而振動松弛因子需要根據流場的局部環(huán)境參數(shù)決定。MILLKAN-WHITE[10]在高溫激波管內通過干涉儀觀測了分子振動松弛過程,擬合了在10 000 K溫度范圍內,振動松弛時間τv與壓強p和溫度T的關系式:

(6)

式中:μr為分子碰撞對的折合質量。

為保證上式在高溫下的松弛時間切合實際,PARK提出了高溫修正因子τs[11]。而在DSMC方法中,振動松弛因子不允許小于轉動松弛因子,因此本文設置了當Zv

此外由于CO2分子具有多個振動模態(tài),為了簡化多振動弛豫時間的處理,本文依據CAMAC[12]的實驗結果,將CO2分子的4種振動模態(tài)以相同的速率松弛:

(7)

圖2為CO2分子3種模態(tài)的弛豫時間和壓力的乘積與溫度的關系圖。由于彎曲模態(tài)的振動溫度最低,因此低溫下大部分的振動能量都以這種模態(tài)存在,這導致低溫下彎曲模態(tài)的弛豫速率和CAMAC的弛豫速率相差很小;而在高溫下,3種模態(tài)的弛豫時間將會趨于一致,它們與CAMAC的弛豫速率相差也不太大。

圖2 CO2各模態(tài)振動弛豫時間Fig.2 Relaxation time of various modes of CO2 vibration

2 碰撞后內能分配方法

本文根據L-B模型,將各內能模式逐一與平動能進行重新分配,這一過程中通過調整彈性碰撞與非彈性碰撞的松弛因子,來控制弛豫過程使得結果與實驗相符合。下面以兩個CO2分子碰撞后內能的分配為例,描述具體實現(xiàn)方法。

設ET為兩個分子的碰前相對平動能;EVL和EVM分別是分子L和M的碰前總振動能,EVLm和EVMm為振動模態(tài)m(m=1,2,3,4)的振動能;ERL和ERM為碰前轉動能;Mm為1.2節(jié)中各模態(tài)的最大振動能級。碰撞后依次考慮分子L和M的能量分配:

先考慮分子L(振動模態(tài)m=1,2,3,4):

①首先判斷當Zv>fR時,才會發(fā)生能量分配,否則跳到步驟⑨;

②從m=1開始,設定初始碰撞總能量為EC=ET+EVLm;

③按式(3)確定振動能為EC時的振動能級,向下取整為imax;

④由im=fR(imax+0.999…),求得此振動模態(tài)的新能級和振動能E′>VLm;

⑤若im>Mm,返回步驟④;

⑥,根據分子對黏性-溫度指數(shù)ωLM,計算概率P=(1-E′>VLm/EC)(1.5-ωLM)若P

⑦將剩余平動能ET返回第一步分配給其余模態(tài)的振動能,ET=EC-E′>VLm;

⑧所有振動模態(tài)能量分配結束,按照上述步驟邏輯進行T-R能量分配;

3 化學反應理論及實現(xiàn)

本文以純CO2氣體作為初始組分,高溫下,氣體分子將會發(fā)生離解反應、置換反應和復合反應等,主要涉及了CO2、CO、O2、C、O五種組分,共計17種基元反應,化學反應式如下,各反應的反應常數(shù)參見文獻[13]。

CO2+GCO+O+G
CO+GC+O+G
O2+GO+O+G
CO2+OO2+CO
CO+OO2+C

(8)

式中:G代表第三體分子或者催化粒子。

3.1 離解反應

對于簡單諧振子模型,一般認為當分子的碰后振動能級大于離解能所對應的振動能級時,離解發(fā)生,即:

imax>Θd/Θv

(9)

式中:imax是碰撞能EC(平動能+振動能)對應的最高振動能級。

由于本文采用的非諧振子模型,所以imax需要根據式(3)進行截斷處理,并且由于CO2分子有4種振動模態(tài),需要按照第2節(jié)中的邏輯,當總振動能εv>kΘd時,離解才發(fā)生。

3.2 復合反應

復合反應可以看作離解反應的逆反應,反應過程是一個三體碰撞過程。復合反應發(fā)生的概率,依據平衡碰撞理論[14],可簡化為

(10)

式中:σR為化學反應截面,σT為碰撞截面,nM為第三體分子的數(shù)密度,?和γ為不同反應的常數(shù)值[15]。利用網格累計概率法判斷,當網格單元中的累加概率大于1時反應發(fā)生。

3.3 置換反應

置換反應是典型的雙分子反應,其可根據反應能量的變化方向,分為吸熱置換反應和放熱置換反應兩類。本文對吸熱的置換反應采用TCE模型[16],最終簡化得到反應發(fā)生的概率表達式為

(11)

式中:Ea為反應的活化能,Pf為吸熱的置換反應發(fā)生概率。

對于放熱的置換反應,由于其本身的活化能較低,反應過程又向周圍散發(fā)熱量,反應發(fā)生的概率也相對較大。本文簡化得到的化學反應抽樣幾率公式[17]為

(12)

4 計算方法驗證

4.1 計算條件

本文選取火星探路者號(MPF)為研究對象,其外形是一個頂角為140.38°的球形鈍頭大底,后接一個46.63°的錐形,探測器最大直徑為2.65 m。其參數(shù)和幾何外形如圖3所示。來流參數(shù)見表2,參考文獻[18],選取該探測器在65 km處來流為模擬狀態(tài)。

圖3 MPF探測器外形Fig.3 Shape of MPF detector

表2 來流參數(shù)Table 2 flow Parameters of incoming flow

來流氣體為純CO2氣體,壁面條件為完全漫反射,不考慮催化。計算網格采用具有更好的貼體性的非結構網格[19]。分子模型采用變徑硬球(Variable Hard Sphere,VHS)模型。 化學反應模型采用TCE模型,轉動松弛數(shù)設為5,振動松弛數(shù)依據1.3小節(jié)設置。

4.2 結果分析

計算結果如圖4所示。圖4(a)、(b)為流場速度分布云圖和壓力云圖,MPF探測器迎風面出現(xiàn)了一道強烈的弓形激波面,氣體的壓縮效應明顯。圖4(c)～(f)分別為流場總溫、平動溫度、轉動溫度和振動溫度的分布云圖。激波后,氣體被強烈壓縮,分子數(shù)密度增加導致分子間的碰撞頻率加快,激波附近出現(xiàn)局部高溫區(qū),且平動、轉動和振動溫度差異很大,流場中出現(xiàn)了顯著的熱力學非平衡現(xiàn)象。

圖4 探測器流場分布云圖Fig.4 Cloud map of detector flow field distribution

圖5為沿駐點線附近的4種溫度分布圖。可以看出,由于遠端來流中分子密度低、溫度低,流場處于平衡態(tài);平動溫度的峰值遠高于轉動、振動溫度,轉動溫度峰值高于振動溫度,且各溫度的峰值逐步靠近壁面,這是由于平動松弛速率高于內能松弛速率,轉動弛豫速率高于振動弛豫速率導致的;隨著流場靠近壁面,分子密度變大,能量交換速率加快,3種能量模式的溫度逐步趨于一致。這4種溫度的變化趨勢及大小關系初步驗證了本文將CO2離解能均分到4種振動模態(tài),從而處理分子振動激發(fā)、化學反應及能量計算方法的正確性。

圖5 駐點線溫度Fig.5 Temperature of stagnation point line

圖6和圖7分別給出了MPF探測器物面壓力和迎風面熱流密度沿探測器上半壁面(從前駐點沿上物面至后駐點)的變化圖,橫坐標S表示壁面沿前駐點的距離。由圖中可以看出,本文的計算結果與文獻[18]中給出的結果基本一致,在前駐點迎風位置和探測器肩部拐角處的壓力、熱流密度與參考文獻幾乎一致,且在探測器尾部出現(xiàn)了0.1 Pa左右的近乎“真空”區(qū)域,符合實際情況。

圖6 物面壓力Fig.6 Surface pressure

圖7 迎風面熱流密度Fig.7 Heat flux density on the windward side

圖8為沿駐點線上5種組分的摩爾分數(shù)變化圖。CO2的摩爾分數(shù)沿駐點線逐漸減小,物面駐點位置的含量僅為0.1%,離解較為徹底。CO分子和O原子隨著CO2的離解,摩爾分數(shù)隨之升高,趨勢較一致。但隨著分子來到激波高溫區(qū),CO發(fā)生離解,使得O原子摩爾分數(shù)繼續(xù)升高,CO分子的摩爾分數(shù)隨之降低,C原子得以產生。當分子靠近壁面時,由于冷壁效應(壁面的溫度遠低于流場溫度),CO分子的復合概率增加,C原子和O原子在壁面附近生成CO,致使CO的摩爾分數(shù)重新回升,C原子和O原子的比例隨之下降。整個過程中,由于高溫,很難發(fā)生O2的復合反應,所以探測器前端O2含量幾乎為0。本文中各組分比例的變化規(guī)律與理論相符,且O原子的CO分子的最大摩爾比例分別為0.48和0.41,與文獻[18]給出的結果相近,驗證了本文計算方法的正確性。

圖8 駐點線各組分摩爾百分數(shù)Fig.8 Molar percentage of each component on the stationary point line

5 不同參數(shù)對化學反應的影響

5.1 密度影響

為了進一步模擬火星探測器進入過程中面臨的真實氣動環(huán)境的影響,選取文獻[18]中85～55 km海拔高度下的來流計算條件,探究密度變化對流場特性的影響,來流參數(shù)見表3。

表3 不同高度計算條件Table 3 Calculation conditions for different heights

圖9為進入火星大氣過程中,不同海拔高度下沿探測器前駐點線的流場總溫變化趨勢圖。各個高度下的總溫變化趨勢一致,高溫區(qū)域都集中在激波與物面之間的過渡區(qū)域。隨著海拔下降,從75～55 km高度處,總溫峰值溫度逐漸下降,而85 km處總溫分布卻高于75 km處,這是因為總溫是平動溫度、轉動溫度和振動溫度與各溫度自由度乘積之和的平均值。

圖9 不同高度下流場沿駐點線總溫分布Fig.9 The total temperature distribution of the flow field along the stationary line at different heights

圖10為不同高度駐點線上CO2分子的摩爾分數(shù),其中55 km與65 km海拔下,CO2分子在激波內部幾乎離解完全,在到達駐點位置時,含量在1%以下,并且由于55 km海拔處的分子壓縮效應更強,所以該駐點線上CO2分子開始離解的位置比65 km海拔處更靠近物面。海拔75 km和85 km處激波后方的CO2分子已經無法完全離解,75 km海拔駐點處CO2分子的摩爾分數(shù)為15%左右,這些未離解的處于高振動激發(fā)態(tài)的CO2分子隨流場運動到探測器肩部以后,將導致探測器尾部的流場總溫升高。

圖10 不同高度下駐點線CO2摩爾分數(shù)Fig.10 CO2 molar fraction at different altitudes of the stationary line

圖11 不同高度下物面特性參數(shù)Fig.11 Surface characteristics parameters at different heights

5.2 速度影響

火星大氣中的音速僅為240 m/s,選取65 km海拔處的火星大氣條件,來流馬赫數(shù)取25、27、29、31,計算模型仍為MPF探測器,其余參數(shù)見表3。

圖12為不同馬赫數(shù)下駐點線上CO2分子的摩爾比例。隨著馬赫數(shù)的增大,流場中的能量增加,提供了更多的能量供CO2分子發(fā)生離解,所以當來流數(shù)密度相同時,馬赫數(shù)越大,流場中發(fā)生的化學反應越激烈,CO2分子的離解率越高。所以在激波區(qū)域相同的橫坐標位置,馬赫數(shù)越高的流場,CO2分子的摩爾比例越低。圖13為不同馬赫數(shù)下物面壓力與迎風面熱流密度變化曲線圖。

圖12 不同馬赫數(shù)下駐點線上CO2的摩爾比例Fig.12 Molar ratio of CO2 molecules on the stationary line at different Mach numbers

圖13 不同馬赫數(shù)下物面特性參數(shù)Fig.13 Different Mach number object surface characteristics parameters

隨著馬赫數(shù)增加,物面壓力和熱流密度均增大,這與流場總溫的變化趨勢一致。但是在肩部往后的背風區(qū)域,4種流場中的物面壓力及熱流密度迅速減小且趨于一致。從駐點處的熱流峰值來看,隨著馬赫數(shù)的增加,物面熱流密度的增幅幾乎保持一致。圖13(a)中,4種流場的物面壓力系數(shù)幾乎沒有差別,這符合牛頓流體在超音速下的流動特性;從圖13(d)物面熱流密度系數(shù)變化圖來看,除了頭部駐點附近的熱流密度系數(shù)有些許偏差外,4種馬赫數(shù)流場的迎風面熱流密度系數(shù)的變化趨勢以及數(shù)值大小均保持一致,這說明馬赫數(shù)對物面的熱流密度分布情況以及分子與物面的傳能效率影響較小,僅僅對物面熱流密度影響較大。

6 結束語

本文基于非結構網格技術,利用DSMC方法,研究了CO2分子各振動模態(tài)的振動能級分配及振動能計算方法,對MPF火星探測器超高聲速進入CO2稀薄氣體環(huán)境的高溫真實氣體效應開展數(shù)值模擬。

①該方法可以有效模擬高超聲速條件下的CO2真實氣體效應。

②隨著探測器進入過程中海拔高度的下降,化學反應對流場的影響程度逐漸增強,流場中的壓力和物面熱流密度都有所升高,流場總溫隨著海拔高度的降低,呈現(xiàn)出了先增加后減小的變化趨勢。

③隨著來流馬赫數(shù)的增大,化學反應的發(fā)生越來越劇烈,化學反應加快了能量模式間的能量傳遞。探測器表面的壓力及熱流密度隨來流速度的增大而增大,來流速度并不會影響物面熱流密度的分布及分子與壁面的傳能效率,僅僅影響熱流密度的數(shù)值大小。