DSMC量子動理學化學反應模型的數值模擬

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

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

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

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

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

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

1 量子動理學模型

1.1 離解反應

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

imax>Θd/Θv

(1)

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

(2)

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

(3)

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

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

(4)

(5)

(6)

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

1.2 置換反應

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

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

(7)

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

解析的反應速率表達式為

(8)

(9)

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

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

2 絕熱單元熱泳測試

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

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

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

(10)

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

(11)

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

2.1 平衡反應測試

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

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

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

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

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

2.2 非平衡反應測試

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

工況2～4對全部19組空氣反應同時測試。組分均為21%的氧氣和79%的氮氣，僅初始溫度和數密度不同。初始的單元溫度分別是30 000、20 000、10 000 K。圖4給出了不同工況下組分濃度與單元溫度隨時間的變化。隨著溫度的降低，氮氣的離解變弱，平衡后的組分濃度分別約為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

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

3 數值算例

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

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

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

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

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

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

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

4 結 論

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

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

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

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

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