李錦，耿湘人，陳堅強，江定武，李紅喆

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

高超聲速飛行器，經(jīng)常在稀薄氣體環(huán)境中飛行，面臨極端速度和高度[1]。當鈍頭體高超聲速飛行器以近乎軌道速度再入大氣層時，其頭部激波之后將出現(xiàn)超過20 000 K的高溫。這樣的流動條件會激發(fā)分子的振動模態(tài)，并導致離解反應的發(fā)生。置換反應和原子的復合反應將在冷的催化壁面附近以及氣體繞過飛行器肩部的膨脹流動中發(fā)生[2]。高超聲速稀薄流動中的化學反應過程將顯著改變飛行器周圍流場的特征，影響飛行器表面的溫度、壓力和熱流以及飛行器的氣動性能[2]。如果不考慮化學反應，就無法精確描述飛行器周圍的流動物理。

利用實驗手段研究高超聲速稀薄流動中的化學反應效應非常昂貴且難以實施。而分子模型非常適合用來研究化學反應氣體流動，因為氣相化學反應的基本理論主要與分子層次的過程相關(guān)。直接模擬蒙特卡羅(Direct Simulation Monte Carlo,DSMC)方法可用來模擬這些復雜流動，化學反應模塊能夠在分子層次插入[3]。DSMC方法中的化學反應模型有很多種，最常用的是Bird教授提出的總碰撞能(Total Collision Energy, TCE)模型[3]。不過TCE模型強烈依賴于宏觀的化學反應速率系數(shù)實驗數(shù)據(jù)。由于高溫區(qū)域測量數(shù)據(jù)的缺失和不確定性，高超聲速流動模擬中常采用的Arrhenius化學反應速率系數(shù)或是從測量的低溫平衡化學反應速率系數(shù)中插值而來，或是從輻射數(shù)據(jù)估算而來。對有的反應，不存在實驗測量數(shù)據(jù)，只能從類似的反應中通過近似估算[4]。而且在TCE模型中，從宏觀數(shù)據(jù)到微觀反應概率的推導使用了平衡理論，這也使得它可能對分布函數(shù)顯著偏離平衡形式的氣體反應不夠精確[5]。為了擺脫這一束縛，發(fā)展更高效的方法，基于量子振動模型，Bird教授又提出了新的量子動理學(Quantum Kinetic,QK)模型[5]，用于從微觀層次直接計算化學反應氣體流動。QK模型不再需要宏觀的化學反應速率系數(shù)數(shù)據(jù)，只依賴于碰撞分子的性質(zhì)，包括可用的碰撞能、離解能和量子化的振動能級[6]。該模型將化學反應過程與碰撞分子的振動模態(tài)能量關(guān)聯(lián)起來。這對于深空探測等難以獲取可靠實驗數(shù)據(jù)的領(lǐng)域可能是一個值得嘗試的選擇。

QK模型已經(jīng)成功吸引了許多國外研究者對DSMC方法的關(guān)注[6-9]，但是在國內(nèi)相關(guān)研究工作還處于起步階段。陳浩等[10]評估了QK模型在氮氧離解復合反應中的表現(xiàn)。近年來國內(nèi)在DSMC方法方面的研究主要集中在時間步長、并行算法[11-13]、碰撞模型[14]、統(tǒng)計誤差[15]以及應用[16-19]等方面。

因此，本文在自研DSMC方法計算程序RariHV中建立QK化學反應模型，通過絕熱單元熱泳測試和高超聲速繞圓柱流動對QK模型的性能進行詳細的計算評估。RariHV是中國空氣動力研究與發(fā)展中心計算空氣動力研究所自主研發(fā)的一款基于DSMC方法的計算軟件，主要用于高超聲速稀薄氣體流動的氣動特性數(shù)值模擬，支持三維復雜外形大規(guī)模并行計算[14]。

1 量子動理學模型

1.1 離解反應

QK模型中離解反應發(fā)生的條件是，如果碰撞能量足夠高，超過了離解所需能量，那么反應就會發(fā)生。離解反應AB+M→A+B+M(AB代表分子；A和B代表原子；M可以是分子，也可以是原子)發(fā)生的條件為

imax>Θd/Θv

(1)

式中：Θd為分子的離解特征溫度；Θv為振動特征溫度；imax為碰撞能Ec對應的最大振動能級。

(2)

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

(3)

式中：為截斷取整符號；εtr,AB,M為碰撞對的平動能；εv,AB為分子的振動能;k為Boltzmann常數(shù)。

QK模型的一個特點是盡管其執(zhí)行過程并不需要氣體處于平衡態(tài)，但是如果假設(shè)氣體處于平衡態(tài)的話，那么可以推導出相應的反應速率kf(T)的解析解。對于可變硬球(Variable Hard Sphere, VHS)模型而言，有

(4)

(5)

(6)

式中：ε為對稱因子，如果M分子與AB分子相同，則ε=2，反之ε=1；T為溫度；rref為在參考溫度Tref時的分子半徑；mr為簡并質(zhì)量；ω為黏性溫度指數(shù)；Q(a,x)=Γ(a,x)/Γ(a)為不完全Gamma函數(shù)；zv(T)=1/[1-exp(-Θv/T)]為振動配分函數(shù)。

1.2 置換反應

置換反應的現(xiàn)象學過程類似于離解模型。置換反應A+B→C+D(A和C為分子，B和D為原子)發(fā)生的概率為

Pr=(1-Ea/Ec)3/2-ω

(7)

式中：Ea為反應活化能，當Ea/k小于Θv時，分母中的求和可近似等于1。

解析的反應速率表達式為

(8)

(9)

本文只計算離解和置換反應。關(guān)于復合反應可以參考文獻[5]。TCE模型的介紹參見Bird教授的專著[3]，這里不再贅述。

TCE模型將總碰撞能量作為確定反應概率的主要控制參數(shù)，沒有考慮化學反應對能量依賴的選擇性。與TCE模型不同，QK模型將碰前的碰撞對平動能以及所考慮分子的振動能作為碰撞能量參與碰后振動模態(tài)能量分配，強調(diào)了振動能量在促進反應方面的重要性。這與廣義Larsen-Borgnakke模型中的平動-振動能量分配是一致的[7]。

2 絕熱單元熱泳測試

DSMC方法中碰撞計算只在單元內(nèi)進行，而化學反應模塊是碰撞模塊的子模塊，因此對化學反應模塊的測試可以只考慮一個網(wǎng)格單元。這樣做最大程度避免了幾何外形、邊界條件、網(wǎng)格分區(qū)等復雜性所帶來的干擾，有利于分析查找問題。分子在單元內(nèi)自由運動、碰撞，這樣的測試稱為熱泳。具體的測試方法分為2種：一種是平衡測試，另一種是非平衡測試。

所謂平衡測試，就是指并不真正發(fā)生化學反應，碰撞對并不真的進行離解、復合或者置換等反應，而是僅僅在程序中對發(fā)生碰撞和反應的數(shù)目進行計數(shù)。反應概率的模擬值可以根據(jù)模擬中發(fā)生化學反應的碰撞數(shù)與總碰撞數(shù)的比值來計算。反應速率則根據(jù)理論公式計算。單元內(nèi)能量是不變的，溫度也不會改變。反應概率和反應速率存在理論平均值，可以與模擬值進行對比。

反應速率的理論計算公式為[20]

(10)

反應概率的理論計算公式為

(11)

化學反應的速率系數(shù)數(shù)據(jù)采用5組分(O2、N2、O、N、NO)19反應空氣模型[1]，包括15組離解反應和4組置換反應。需要指出的是，文獻中的這些反應速率系數(shù)并非實驗數(shù)據(jù)，而是根據(jù)QK模型的理論計算結(jié)果擬合得到的。采用這套數(shù)據(jù)是為了方便QK模型和TCE模型的比較。

2.1 平衡反應測試

平衡反應測試在一個邊長為10-5m的絕熱正方體網(wǎng)格單元內(nèi)進行。單元的6個面均為鏡面反射條件。單元內(nèi)的初始密度ρ可取為ρ/ρd=107，ρd為特征離解密度，它是溫度的函數(shù)，但是其變化范圍較小，計算結(jié)果對其不敏感，通?？梢匀槌?shù)，如1.30×105kg/m3。單元內(nèi)布置50 000個 模擬分子，時間步長取為1.52×10-9s。碰撞模型選擇VHS模型。內(nèi)能和平動能之間的松弛采用Larsen-Borgnakke模型，轉(zhuǎn)動和振動松弛碰撞數(shù)都取為1。

圖1給出了氧氣和氮氣離解反應的平衡反應速率系數(shù)的理論值與模擬值的對比情況。溫度范圍為5 000～13 000 K，只考慮正向反應。從圖中可以看出，反應速率系數(shù)隨著溫度的升高而增加，模擬值與理論值基本吻合。模擬值相比理論值略低，特別是在溫度為5 000 K時。這是因為QK模型中離解反應是直接根據(jù)碰撞能量以及離解能量的比較來進行的，并不像TCE模型那樣計算反應概率，因此，模擬值與理論值無法精確匹配。氮氣的離解反應速率系數(shù)低于氧氣，這是由于氮氣離解溫度高導致的。圖2給出了2組置換反應的平衡反應速率系數(shù)隨溫度的變化情況，模擬值與理論值在所有溫度下都高度吻合。本文對其余15組反應也進行了相應的測試，結(jié)果都吻合良好，限于文章篇幅，這里沒有給出。

圖1 離解反應的平衡反應速率系數(shù)Fig.1 Equilibrium reaction rate coefficient for dissociation reactions

圖2 置換反應的平衡反應速率系數(shù)Fig.2 Equilibrium reaction rate coefficient for exchange reactions

平衡反應的測試結(jié)果表明，在一定精度范圍內(nèi)，QK模型可以給出準確的化學反應速率系數(shù)。

2.2 非平衡反應測試

非平衡反應測試采用Haas和Mcdonald[21]提出的方法。該方法可以求得反應過程中單元內(nèi)各化學組分的濃度變化以及宏觀溫度隨時間變化的解析解，以便與模擬解對比。表1是測試條件，其余條件與平衡反應測試中相同。

工況2～4對全部19組空氣反應同時測試。組分均為21%的氧氣和79%的氮氣，僅初始溫度和數(shù)密度不同。初始的單元溫度分別是30 000、20 000、10 000 K。圖4給出了不同工況下組分濃度與單元溫度隨時間的變化。隨著溫度的降低，氮氣的離解變?nèi)?，平衡后的組分濃度分別約為20%、50%和76%。氧氣離解比較充分，但離解的速度也隨溫度的降低而變慢。3個工況均以吸熱反應為主，單元溫度隨時間不斷下降。除工況3中氮氣分子和氮原子的組分濃度與理論解有一定偏差外，QK模型預測的組分濃度及單元溫度隨時間的變化與理論解基本吻合。

表1 非平衡反應測試條件Table 1 Test conditions of non equilibrium reactions

圖3 組分濃度和單元溫度隨時間的變化(工況1)Fig.3 Temporal evolution of species concentrations and temperatures (Case 1)

圖4 不同工況下組分濃度和單元溫度隨時間的變化Fig.4 Temporal evolution of species concentrations and temperatures under different cases

非平衡反應測試的結(jié)果表明，在一定精度范圍內(nèi)，QK模型可以準確預測化學反應熱松弛過程。

3 數(shù)值算例

為了進一步測試QK模型的應用能力，本文模擬了高超聲速繞圓柱流動。圓柱直徑為2 m。來流由78.85%的氮氣和21.15%的氧氣組成，來流數(shù)密度為1.43×1020m-3。來流克努森數(shù)約為0.004。來流馬赫數(shù)為24.85，迎角為0°。來流溫度為187 K，圓柱表面溫度為1 000 K，完全漫反射，不考慮催化。更詳細的信息參見文獻[1]。

圖5給出了不同模型預測的溫度場的比較。圖5(a)是無化學反應和QK模型結(jié)果的比較，圖5(b)是TCE模型和QK模型結(jié)果的比較。從圖中可以看出，TCE模型和QK模型計算的流場溫度分布吻合良好。與無化學反應的計算結(jié)果相比，盡管來流克努森數(shù)較小，屬于近連續(xù)流，但由于化學反應，流場最高溫度降低了約3 800 K。

圖6給出了駐點線上的平動、轉(zhuǎn)動及振動溫度分布。從圖中可以看出，溫度非平衡效應十分明顯。圖中各線型(直線、虛線及點劃線)分別代表本文模擬值，與之相對應的各符號(方框、菱形及圓形)為文獻[1]結(jié)果。弓形激波位于圓柱前約25 cm處。氣體經(jīng)過激波后溫度迅速上升，然后在近壁面處下降。平動、轉(zhuǎn)動和振動溫度的峰值分別為約20 000、13 000、9 000 K。從圖中可以看出，除了轉(zhuǎn)動溫度的峰值略低之外，本文計算結(jié)果與文獻[1]結(jié)果吻合良好。

圖5 不同模型溫度場的對比Fig.5 Comparison of temperature contours by different models

圖7給出了駐點線上5組分的數(shù)密度分布情況。圖中各線型為本文模擬值，與之相對應的各符號為文獻[1]結(jié)果。從圖中可以看出，QK模型可以準確捕捉到化學反應流場中的組分分布。

圖6 駐點線上溫度對比Fig.6 Comparison of temperatures along stagnation line

圖7 駐點線上組分數(shù)密度對比Fig.7 Comparison of number densities along stagnation line

4 結(jié) 論

本文在自研DSMC軟件中實現(xiàn)了新的量子動理學化學反應模型，并對其進行了計算精度的測試，得到以下初步結(jié)論：

1) QK模型在絕熱單元平衡和非平衡熱泳中表現(xiàn)良好，模擬結(jié)果與理論解吻合。這表明，QK模型可以給出準確的化學反應平衡速率并預測化學反應非平衡過程。

2) 本文將QK模型應用于高超聲速繞圓柱計算中，所得計算結(jié)果無論是宏觀的駐點線上溫度，還是微觀的組分數(shù)密度分布等，都與文獻[1]結(jié)果吻合良好。這表明QK模型具有計算真實高超聲速化學反應流動的能力。

從本文初步的評估來看，QK模型具有優(yōu)良的特性，可應用于高超聲速化學反應流動計算。QK模型不依賴宏觀的化學反應速率數(shù)據(jù)，未來有希望應用于深空探測等領(lǐng)域，當然其性能也還需要進一步的考察和檢驗。

下一步，將應用QK模型進行火星探測器稀薄氣動特性的計算。