陳 浩, 李林穎, 張 斌,2,*, 劉 洪(. 上海交通大學(xué) 航空航天學(xué)院, 上海 200240; 2. 上海交通大學(xué) 數(shù)學(xué)系, 上海 200240)

0 引 言

高超聲速再入飛行通常發(fā)生在稀薄氣體環(huán)境中，由于激波層和高空低密度的特征，在激波下游會產(chǎn)生顯著的熱化學(xué)非平衡松弛現(xiàn)象[1]。流場中的化學(xué)反應(yīng)對飛行器氣動性能和表面?zhèn)鳠嵊兄@著影響，多原子振動激發(fā)和離解反應(yīng)的發(fā)生對高超聲速飛行器來說是很重要的特征，低溫下的振動激發(fā)和有限的氧氣離解仍然很重要[2]。

高速稀薄氣體環(huán)境下，直接模擬蒙特卡洛(DSMC)方法[3]在模擬熱化學(xué)非平衡方面具有很好的優(yōu)勢，該方法是目前從自由分子流區(qū)到稀薄過渡流區(qū)用于數(shù)值分析氣體流動問題的主要模擬工具。DSMC方法的化學(xué)反應(yīng)模型發(fā)展有著近40年的歷史，總能碰撞(TCE, total collision energy)模型是最為廣泛使用的一種[4]。該模型基于碰撞理論，將標準Arrhenius形式的連續(xù)速率方程轉(zhuǎn)變?yōu)榉磻?yīng)碰撞截面，因此需要每一個化學(xué)反應(yīng)的Arrhenius速率系數(shù)[3-4]。雖然DSMC化學(xué)反應(yīng)模型發(fā)展比較緩慢并且種類不多，但在精細化學(xué)反應(yīng)模擬方面的研究很多。這一類模型直接采用QCT(Quasi classical trajectory)計算結(jié)果[5]，或者采用態(tài)-態(tài)碰撞截面實驗數(shù)據(jù)[6]來描述碰撞截面。Bird[7]于2011年提出了一種全新的化學(xué)反應(yīng)模型：量子-動理學(xué)(Q-K, quantum-kinetic)模型。該方法基于量子振動模型，幾乎所有的化學(xué)反應(yīng)概率都基于碰撞分子的振動能態(tài)而描述，不需要采用任何的宏觀速率方程，也不存在氣體平衡態(tài)的假設(shè)[2,4,7]。

TCE模型不能細致的描述內(nèi)能(振動能、轉(zhuǎn)動能)對反應(yīng)過程的貢獻，針對不同反應(yīng)采用不同的碰撞總能來滿足對連續(xù)反應(yīng)速率的準確模擬[2,4,8-9]。而Q-K模型將分子振動能態(tài)與化學(xué)反應(yīng)直接相關(guān)聯(lián)，根據(jù)平衡原理得到精準的模擬結(jié)果。本文關(guān)注的便是兩種模型在這一類細致模擬的問題上會有怎樣的不同，因此本文旨在評估Q-K模型模擬平衡態(tài)、非平衡離解反應(yīng)速率的能力，完成Q-K模型和TCE模型在離解度、非平衡反應(yīng)速率和振動非平衡效應(yīng)方面的對比研究工作。

下一節(jié)簡略地介紹了Q-K模型的特點，以及離解、復(fù)合反應(yīng)的在DSMC程序中的實施策略。結(jié)果與討論部分分別在平衡態(tài)離解度和非平衡離解問題上對Q-K模型進行了數(shù)值評估，同時對比分析了Q-K和TCE模型在模擬非平衡離解速率和離解速率中的溫度非平衡效應(yīng)兩方面的差異，并與經(jīng)典Park雙溫度模型進行了對比分析。

1 Q-K化學(xué)反應(yīng)模型

振動能態(tài)在Q-K模型中具有重要的作用，離解和置換反應(yīng)的發(fā)生都是基于碰撞分子對的振動能態(tài)而描述的，采用量子化碰撞溫度能夠?qū)崿F(xiàn)振動能平衡均分的結(jié)果[10]。Q-K模型采用量子振動和簡諧振子模型，處理振動能松弛時采用“量子碰撞溫度”求解松弛碰撞數(shù)。量子碰撞溫度可寫為：

(1)

式中，imax為碰撞分子對所具有的最大能級，ω為氣體黏性系數(shù)中的溫度指數(shù)。采用量子碰撞溫度Tcoll定義振動松弛碰撞數(shù)如下[11]：

(2)

其中，Zref為參考溫度Tref下的振動松弛碰撞數(shù)，該松弛碰撞數(shù)可通過TCE模型中求解振動碰撞數(shù)的方式獲得[2]：

(3)

式(2)和(3)中的參考溫度Tref通常用振動特征溫度θv表示[2]。

本文所考慮的化學(xué)反應(yīng)類型包括離解反應(yīng)和復(fù)合反應(yīng)兩種，Q-K模型中兩種化學(xué)反應(yīng)的處理方式簡要表達如下小節(jié)所示，詳細的Q-K理論可查閱文獻[7]。

1.1 離解反應(yīng)

離解反應(yīng)的發(fā)生條件可描述為：一次碰撞的能量足夠大以至于能夠達到離解閥值所對應(yīng)的振動能級?？紤]一個離解反應(yīng)AB+T=A+B+T，AB為離解對象，T為碰撞對象，A和B為分子AB的離解產(chǎn)物。碰撞后的最大振動能級imax由方程(4)確定：

(4)

imax>θd/θv

(5)

其中，θd為分子AB的特征離解溫度。

1.2 復(fù)合反應(yīng)

A和B發(fā)生復(fù)合反應(yīng)的條件為：有另外一個分子或者原子進入到碰撞體積Vcoll中，這里的體積為“三元碰撞體積”。因此，復(fù)合反應(yīng)A+B+T=AB+T發(fā)生的概率可描述為：

Prec=nVcoll=naTbVref

(6)

式中，n為分子數(shù)密度，Vref為參考碰撞體積，其值為半徑等于三個碰撞分子(A、B和T)半徑之和的球體的體積；系數(shù)aTb中的溫度一般采用宏觀溫度，則T=Tcell/θv，Tcell為單位網(wǎng)格內(nèi)的宏觀溫度；系數(shù)aTb中的參數(shù)a、b根據(jù)理論平衡常數(shù)進行調(diào)整，推薦值如表1所示。

2 結(jié)果與討論

2.1 基于離解度的Q-K模型評估和對比分析

文獻[7]對Q-K理論直接推導(dǎo)的反應(yīng)速率進行了詳細的驗證，Body教授[2]通過零維熱浴和圓柱擾流數(shù)值模擬將Q-K模型的DSMC采樣結(jié)果與實驗數(shù)據(jù)[12-15]、傳統(tǒng)TCE模型[14]進行了詳盡的對比，證明了Q-K模型在DSMC程序中的可行性，因此本文不再贅述平衡態(tài)反應(yīng)速率的對比結(jié)果。工程上關(guān)注更多的是化學(xué)反應(yīng)模型在模擬氣體組分和反應(yīng)后流場溫度上的精準度，因此我們對標準分子數(shù)密度狀態(tài)下(n=2.687×1025m-3)的空氣進行熱浴實驗?zāi)M，對比Q-K和TCE化學(xué)反應(yīng)模型在模擬雙原子離解度上的精準性。

表1 三元碰撞體積參數(shù)Table 1 Parameters in ternary collision volume

空氣組分為78%的氮氣和22%的氧氣，采用1×104個模擬分子進行熱浴模擬。注意的是，該模擬中只考慮氮氣和氧氣的離解、復(fù)合反應(yīng)，因此不會有NO分子的生成。兩種化學(xué)反應(yīng)模型的對比結(jié)果如圖1所示，圖中每一個點代表當(dāng)前平衡溫度下所對應(yīng)的離解度。圓形和三角形標記分別為Q-K和TCE模型的計算結(jié)果，線條為離解度理論分析值[16]，紅色表示氮氣，藍色表示氧氣。我們發(fā)現(xiàn)，Q-K和TCE模型所得到的結(jié)果幾乎完全一致。相比于氮氣，氧氣表現(xiàn)出極易容易離解的特征。當(dāng)溫度約為5000 K左右時，氮氣開始離解，而氧氣離解度則高達90%。此外，如圖1橢圓形區(qū)域所示，Q-K模型表現(xiàn)出比TCE模型更精確的仿真能力，TCE模型得到的離解度過低，當(dāng)平衡溫度約為6000 K時，離解度相對誤差高達75%。

圖1 不同溫度下的氮氧離解度Fig.1 The degree of dissociation for nitrogen and oxygen

根據(jù)兩種化學(xué)反應(yīng)模型的理論和實現(xiàn)過程，可理性地得到TCE模型在低離解度時產(chǎn)生大誤差的原因。TCE離解反應(yīng)的發(fā)生是建立在振動松弛條件之下的，即離解反應(yīng)的觸發(fā)是由振動能唯一控制的。然而空氣中的化學(xué)反應(yīng)并非由振動能唯一控制[16]，在分子振動能還沒有高度激發(fā)的狀態(tài)下，平動能和轉(zhuǎn)動能對化學(xué)反應(yīng)的影響將變得不可忽略。Q-K模型能得到精準離解度的原因在于，Q-K離解反應(yīng)的觸發(fā)是獨立于振動松弛而設(shè)置的，并且由方程(4)可看出，離解反應(yīng)發(fā)生條件中的最大振動能級imax考慮到了相對平動能對反應(yīng)的貢獻。

2.2 基于非平衡離解的Q-K模型評估和對比分析

通常來說，非平衡效應(yīng)有兩種。一種是速度分布不再滿足Maxwell分布，另一種是內(nèi)能松弛引起的溫度非平衡。本節(jié)討論的非平衡便是由平動溫度、轉(zhuǎn)動溫度和振動溫度所體現(xiàn)的溫度非平衡效應(yīng)。

考慮三種非平衡的情況：1、平動-轉(zhuǎn)動溫度保持在6000 K，振動溫度逐漸上升；2、平動-轉(zhuǎn)動溫度保持在10000 K，振動溫度逐漸上升；3、平動-轉(zhuǎn)動溫度保持在20000 K，振動溫度逐漸升高。在DSMC程序中很容易實現(xiàn)這種溫度非平衡：振動能在碰撞過程中不再與平動能、轉(zhuǎn)動能進行能量交換，而是按照初始給定的振動溫度給予振動能Boltzmann分布。作為比照的參考，本文采用Park[12]發(fā)展的雙溫度模型從理論上描述前述三種非平衡條件下的非平衡反應(yīng)速率。在Park的雙溫度模型中，采用一個平動-轉(zhuǎn)動溫度和振動溫度耦合的幾何平均溫度來替代平衡態(tài)Arrhenius化學(xué)反應(yīng)速率方程中的宏觀溫度，以此來表征振動溫度對化學(xué)反應(yīng)速率的影響，該幾何平均溫度Ta定義如下：

(7)

式中，Tv為振動溫度，T為平動-轉(zhuǎn)動溫度。針對氮氣離解反應(yīng)，Sharma[17]建議取值q=0.3。

針對上述的人為的給定的三種非平衡離解反應(yīng)過程，對Q-K模型、TCE模型和Park雙溫度模型得到的非平衡離解反應(yīng)速率進行了對比。比較結(jié)果如圖2所示，線條表示Park非平衡離解反應(yīng)速率、標記塊為Q-K模型和TCE模型模擬所得到的非平衡速率，具體標記區(qū)別如圖2中圖例所示。

如圖2所示，TCE模型所描述的非平衡離解速率更接近Park理論，而Q-K模型的結(jié)果比Park非平衡速率要大近1～2個量級左右。還可以發(fā)現(xiàn)，隨著平動-轉(zhuǎn)動溫度的增加，非平衡反應(yīng)速率對振動溫度的敏感性變得很弱，其中TCE模型的這種趨勢更為明顯。當(dāng)平動-轉(zhuǎn)動溫度降低時，振動溫度的升高能有效的提高化學(xué)反應(yīng)速率，其中Q-K模型的這種趨勢更為明顯；當(dāng)振動溫度很低時，兩種化學(xué)反應(yīng)模型DSMC模擬結(jié)果與Park雙溫度模型之間的誤差很大，其原因在于：Park雙溫度模型中的非平衡反應(yīng)速率對幾何平均溫度Ta的依賴導(dǎo)致了反應(yīng)速率的快速下降，以及DSMC程序中采用了相對平動能來滿足反應(yīng)所需的激活能。

圖2 Q-K模型和TCE模型的非平衡離解速率Fig.2 The nonequilibrium dissociation rate using Q-K and TCE model

另外，本文還對Q-K模型和TCE模型描述溫度非平衡效應(yīng)的能力進行了對比研究。專著[18]中建議采用比值Kf(T,Tvib)/Kf(T,T)來考察振動非平衡效應(yīng)對離解反應(yīng)速率的影響，本文也采用該比值進行對比分析。隨著振動溫度與平動-轉(zhuǎn)動溫度之間的溫差越大，這種溫度非平衡效應(yīng)也就越強。為得到較大程度的振動非平衡效應(yīng)，本文取定振動溫度為Tvib=6000 K，而平動-轉(zhuǎn)動溫度Ttran-rot從6000 K逐漸上升至20000 K。同樣，采用Park雙溫度模型作為理論參照，并使用兩種q值，一種為Park推薦值[12]：q=0.5，另一種為Sharma推薦值[17]：q=0.3。

圖3 溫度非平衡對離解速率的影響Fig.3 The effect of vibrational nonequilibrium on the dissociation rate

Q-K模型、TCE模型以及Park雙溫度模型的比較結(jié)果如圖3所示，縱坐標為比值Kf(T,Tvib)/Kf(T,T)的對數(shù)結(jié)果，橫坐標為平動-轉(zhuǎn)動溫度。首先，由圖左上角所示，四條曲線都相交于一點(6,0)，該點意味著：當(dāng)平動-轉(zhuǎn)動溫度與振動溫度均等于6000 K時，比值Kf(T,Tvib)/Kf(T,T)為1，即此時并沒有非平衡效應(yīng)。隨著平動-轉(zhuǎn)動溫度的升高，曲線逐漸下降并逐漸趨于穩(wěn)定，這表明振動非平衡效應(yīng)先增強，然后會趨于一定的程度，該趨勢與圖2中所描述的“當(dāng)平動-轉(zhuǎn)動溫度很高時，振動溫度對非平衡反應(yīng)速率的影響變得很小”這一特征相吻合。更重要的是，圖3中四條曲線的對比表明：在描述離解速率中的溫度非平衡效應(yīng)時，Q-K模型與Park雙溫度模型(q=0.5)一致吻合，而TCE模型卻過大地預(yù)估了這種非平衡效應(yīng)。

3 結(jié) 論

基于工程上比較關(guān)心的大氣高超聲速流動中的氮氧離解反應(yīng)，本文對近期提出的量子-動理學(xué)(Q-K, quantum-kinetic)模型[7]在模擬平衡態(tài)離解度、非平衡離解反應(yīng)速率以及溫度非平衡方面進行了評估，并基于廣泛使用的TCE模型進行詳細的對比研究，得到如下結(jié)論：

1) 在模擬平衡態(tài)離解度方面，Q-K模型與理論值精準吻合，而TCE模型在低離解度時具有較大誤差。該結(jié)果表明，Q-K模型中將離解判斷和振動松弛獨立而設(shè)，并采用能級判斷的方式是滿足物理真實的；

2) 在模擬非平衡離解速率方面，兩模型均與理論描述大致吻合，相比而言TCE結(jié)果與Park雙溫度模型更為接近，而Q-K非平衡離解速率則比Park理論大1～2個量級；

3) 在模擬離解速率中的溫度非平衡效應(yīng)方面，Q-K模型與Park雙溫度模型高度吻合，而TCE模型卻過大地預(yù)估了這種非平衡效應(yīng)。

[1]Park C. Review of chemical-kinetic problems future NASA missions, I: Earth entries[J]. Journal of thermophysics and heat transfer, 1993, 7(3): 385-398

[2]Scanlon T J, Boyd I D, et al. Open-source direct simulation Monte Carlo chemistry modeling for hypersonic flows[J]. AIAA Journal, 2015, 53(6): 1670-1680

[3]Bird G A. Molecular gas dynamics and the direct simulation of gas flows[M]. Oxford: Oxford Univ. Press, 1994

[4]Wysong I, Gimelshein S, et al. Comparison of direct simulation Monte Carlo chemistry and vibrational models applied to oxygen shock measurements[J]. Physics of Fluids, 2014, 26(043101): 1-15

[5]Levin D A, Gimelshein S F. Modeling of OH vibrational distributions using molecular dynamics with the direct simulation Monte Carlo method[C]//34thThermophysics Conference, Doi: 10.2514/6.2000-2432

[6]Kim J G, Boyd I D. Monte Carlo simulation of nitrogen dissociation based on state-resolved cross sections[J]. Physics of Fluids, 2014, 26(012006): 1-20

[7]Bird G A. The Q-K model for gas-phase chemical reaction rates[J]. Physics of Fluids, 2011, 23(106101): 1-13

[8]梁杰, 李志輝, 杜波強, 等. 探月返回器稀薄氣體熱化學(xué)非平衡特征性數(shù)值模擬[J]. 載人航天, 2015, 21(3): 295-302

[9]黃飛, 趙波, 程曉麗, 等. 返回器高空稀薄氣動特性的真實氣體效應(yīng)研究[J]. 宇航學(xué)報, 2014, 35(3): 283-290

[10]Bergemann F, Boyd I D. New discrete vibrational energy model for the direct simulation Monte Carlo method[J]. Progress in Astronautics and Aeronautics, 1994, 158: 174-183

[11]Bird G A. A comparison of collision energy-based and temperature-based procedures in DSMC[C]//Rarefied Gas Dynamics: 26thInternational Symposium, 2008, 1084: 245-250

[12]Park C. Nonequilibrium hypersonic aerothermodynamics[M]. New York: A wiley-interscience publication, 1989

[13]Gupta R N, Thompson R A, et al. A review of reaction rates

and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30000 K[R]. Hampton, NASA Langley Research Center, 1990

[14]Carlson A B, Bird G A. Implementation of a vibrationally linked chemical reaction model for DSMC[R]. Hampton, NASA Langley Research Center, 1994

[15]Gallis M A, Bond R B. Assessment of collision-energy-based models for atmospheric species reactions in hypersonic flows[J]. Journal of the Thermophysics and Heat Transfer, 2010, 24(2): 241-253

[16]楊浩森, 張斌, 劉洪, 等. DSMC熱化學(xué)模擬分子數(shù)對氮氧離解反應(yīng)計算影響的分析[J]. 空氣動力學(xué)學(xué)報, 2014, 32(4): 511-517

[17]Sharma S P, Huo W M, Park C. Rate parameters for coupled vibrational-dissociation in a generalized SSH approximation[J]. Journal of Thermophysics and Heat Transfer, 1988, 6: 9-21

[18]吳其芬, 陳偉芳. 高溫稀薄氣體熱化學(xué)非平衡流動的DSMC方法[M]. 長沙: 國防科技大學(xué)出版社, 1999.