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

0 引 言

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

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

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

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

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

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

(1)

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

(2)

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

(3)

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

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

1.1 離解反應(yīng)

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

(4)

imax>θd/θv

(5)

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

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

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

Prec=nVcoll=naTbVref

(6)

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

2 結(jié)果與討論

2.1 基于離解度的Q-K模型評(píng)估和對(duì)比分析

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

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

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

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

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

2.2 基于非平衡離解的Q-K模型評(píng)估和對(duì)比分析

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

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

(7)

式中，Tv為振動(dòng)溫度，T為平動(dòng)-轉(zhuǎn)動(dòng)溫度。針對(duì)氮?dú)怆x解反應(yīng)，Sharma[17]建議取值q=0.3。

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

[9]黃飛, 趙波, 程曉麗, 等. 返回器高空稀薄氣動(dòng)特性的真實(shí)氣體效應(yīng)研究[J]. 宇航學(xué)報(bào), 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ù)對(duì)氮氧離解反應(yīng)計(jì)算影響的分析[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 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é)非平衡流動(dòng)的DSMC方法[M]. 長沙: 國防科技大學(xué)出版社, 1999.