DSMC量子動(dòng)理學(xué)化學(xué)反應(yīng)模型的數(shù)值模擬

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

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

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

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

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

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

1 量子動(dòng)理學(xué)模型

1.1 離解反應(yīng)

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

imax>Θd/Θv

(1)

式中：Θd為分子的離解特征溫度；Θv為振動(dòng)特征溫度；imax為碰撞能Ec對(duì)應(yīng)的最大振動(dòng)能級(jí)。

(2)

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

(3)

式中：為截?cái)嗳≌?hào)；εtr,AB,M為碰撞對(duì)的平動(dòng)能；εv,AB為分子的振動(dòng)能;k為Boltzmann常數(shù)。

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

(4)

(5)

(6)

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

1.2 置換反應(yīng)

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

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

(7)

式中：Ea為反應(yīng)活化能，當(dāng)Ea/k小于Θv時(shí)，分母中的求和可近似等于1。

解析的反應(yīng)速率表達(dá)式為

(8)

(9)

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

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

2 絕熱單元熱泳測(cè)試

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

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

反應(yīng)速率的理論計(jì)算公式為[20]

(10)

反應(yīng)概率的理論計(jì)算公式為

(11)

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

2.1 平衡反應(yīng)測(cè)試

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

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

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

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

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

2.2 非平衡反應(yīng)測(cè)試

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

工況2～4對(duì)全部19組空氣反應(yīng)同時(shí)測(cè)試。組分均為21%的氧氣和79%的氮?dú)猓瑑H初始溫度和數(shù)密度不同。初始的單元溫度分別是30 000、20 000、10 000 K。圖4給出了不同工況下組分濃度與單元溫度隨時(shí)間的變化。隨著溫度的降低，氮?dú)獾碾x解變?nèi)酰胶夂蟮慕M分濃度分別約為20%、50%和76%。氧氣離解比較充分，但離解的速度也隨溫度的降低而變慢。3個(gè)工況均以吸熱反應(yīng)為主，單元溫度隨時(shí)間不斷下降。除工況3中氮?dú)夥肿雍偷拥慕M分濃度與理論解有一定偏差外，QK模型預(yù)測(cè)的組分濃度及單元溫度隨時(shí)間的變化與理論解基本吻合。

表1 非平衡反應(yīng)測(cè)試條件Table 1 Test conditions of non equilibrium reactions

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

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

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

3 數(shù)值算例

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

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

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

圖5 不同模型溫度場(chǎng)的對(duì)比Fig.5 Comparison of temperature contours by different models

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

圖6 駐點(diǎn)線上溫度對(duì)比Fig.6 Comparison of temperatures along stagnation line

圖7 駐點(diǎn)線上組分?jǐn)?shù)密度對(duì)比Fig.7 Comparison of number densities along stagnation line

4 結(jié) 論

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

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

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

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

下一步，將應(yīng)用QK模型進(jìn)行火星探測(cè)器稀薄氣動(dòng)特性的計(jì)算。