鄧 利，馬 虎，武曉松，周長(zhǎng)省

(南京理工大學(xué)機(jī)械工程學(xué)院,江蘇 南京 210094)

爆震波是激波和燃燒波的強(qiáng)烈耦合，涉及爆震燃燒的流動(dòng)存在多個(gè)時(shí)間尺度，容易出現(xiàn)嚴(yán)重的剛性源項(xiàng)，使數(shù)值求解流動(dòng)方程存在困難。

目前處理剛性源項(xiàng)問(wèn)題的方法主要有分裂求解和耦合求解兩類。采用顯式耦合方法求解時(shí)，時(shí)間步長(zhǎng)需要足夠小以滿足穩(wěn)定性要求，否則可能導(dǎo)致計(jì)算不穩(wěn)定甚至出現(xiàn)非物理解[1]。Bussing等[2-3]針對(duì)反應(yīng)流中出現(xiàn)的多尺度問(wèn)題，提出了點(diǎn)隱算法，將反應(yīng)與流動(dòng)的積分時(shí)間步長(zhǎng)統(tǒng)一到偽時(shí)間步長(zhǎng)上，消除了源項(xiàng)的剛性。Zhong[4]推廣了Rosenbrock[5]從半隱出發(fā)求解剛性常微分方程組的方法，采用不同的方法分別求解方程的剛性項(xiàng)和非剛性項(xiàng)，增強(qiáng)了計(jì)算穩(wěn)定性，提高了時(shí)間方向的精度，并在化學(xué)反應(yīng)流的模擬中得到了成功應(yīng)用。隱式處理剛性源項(xiàng)時(shí)，雖然很好地解決了剛性問(wèn)題，但由于涉及矩陣求逆運(yùn)算，隨著反應(yīng)組分增多，計(jì)算量急劇增大。

Leveque等[6]利用時(shí)間分裂法和MacCormack預(yù)測(cè)校正方法求解了含剛性源項(xiàng)的雙曲型方程，認(rèn)為分裂算法在空間分辨率足夠的情況下，計(jì)算效果比耦合算法稍好。Oran等[7]對(duì)剛性源項(xiàng)方程的積分方法做了全面的論述。一步法雖然在求解剛性方程組時(shí)穩(wěn)定性要求嚴(yán)格，但由于只需已知變量值，求解相對(duì)方便，并針對(duì)特征值為正且相當(dāng)大的物理剛性問(wèn)題時(shí)，使用一步法是比較有利的。Young等[8]采用逼近法與預(yù)測(cè)校正思想來(lái)處理化學(xué)反應(yīng)源項(xiàng)，對(duì)剛性和普通常微分方程組都具有較好的計(jì)算性能，且方便與時(shí)間步長(zhǎng)估算方法[9]配合使用，求解效率高，但容易出現(xiàn)質(zhì)量分?jǐn)?shù)不守恒的問(wèn)題，需要在流場(chǎng)推進(jìn)時(shí)不斷修正。Mott等[10-11]在逼近法的基礎(chǔ)上發(fā)展了αQSS方法，通過(guò)迭代校正組分值來(lái)保證質(zhì)量守恒，比逼近法更加穩(wěn)定和精確，并針對(duì)碳?xì)淙剂系娜紵捅鹉M取得了令人滿意的效果。

劉瑜[12]對(duì)比了VODE(variable-coefficient ordinary differential equation solver)、梯形公式[13]和αQSS方法處理化學(xué)反應(yīng)剛性源項(xiàng)的能力，認(rèn)為在計(jì)算精度基本相同的情況下，梯形公式和αQSS方法效率更高。劉世杰等[14]利用梯形公式和αQSS方法對(duì)爆震波胞格進(jìn)行了模擬，認(rèn)為使用αQSS方法時(shí)計(jì)算的化學(xué)反應(yīng)更完全，放熱量更大。目前對(duì)剛性源項(xiàng)處理方法在化學(xué)反應(yīng)流動(dòng)模擬中的適用性，以及各方法之間的聯(lián)系方面還缺少一定的認(rèn)識(shí)。本文中在對(duì)各源項(xiàng)處理方法的穩(wěn)定性分析基礎(chǔ)上，分析各種方法之間的關(guān)系以及對(duì)化學(xué)反應(yīng)源項(xiàng)的處理能力，并通過(guò)具體算例比較不同處理方法的計(jì)算效率。

1 計(jì)算方法介紹

1.1 控制方程

激波誘導(dǎo)燃燒為激波和燃燒波耦合的結(jié)果，可忽略黏性的影響，因此使用有限速率化學(xué)反應(yīng)的歐拉方程描述為:

(1)

1.2 數(shù)值方法

具有二階時(shí)間精度的Strang分裂格式表示為：

(2)

式中：算符Lf為流動(dòng)求解算子，Lr為化學(xué)反應(yīng)的求解算子，Δt為流場(chǎng)推進(jìn)的時(shí)間步長(zhǎng)，上標(biāo)n為求解時(shí)間步數(shù)。

一步法為常微分方程初值問(wèn)題的經(jīng)典求解方法，求解變量n+1時(shí)刻的值時(shí)，只需n時(shí)刻的信息。解耦后的化學(xué)源項(xiàng)方程如下：

(3)

由于Strang分裂對(duì)算子求解精度的要求，選取2階Runge-Kutta法作為化學(xué)反應(yīng)求解算子。

在利用逼近法求解化學(xué)反應(yīng)源項(xiàng)時(shí)，將組分i的凈生成率寫(xiě)成如下形式：

(4)

式中：qi為組分i在反應(yīng)中的生成項(xiàng)，τi為組分i的消耗特征時(shí)間。則剛性方程組中組分i在n+1時(shí)刻的密度為:

(5)

(6)

αQSS方法在逼近法的基礎(chǔ)上進(jìn)一步改進(jìn)其迭代形式，在求解時(shí)將組分凈生成率寫(xiě)為：

(7)

式中：pi為特征消耗時(shí)間的倒數(shù)，為了形式簡(jiǎn)便略去下標(biāo)，如果qi與pi在時(shí)間段(t,t+Δt)之內(nèi)為常數(shù)時(shí)，則可積分式(7)得到精確解：

(8)

整理式(8)得到：

(9)

(10)

式中：pΔt→∞時(shí)，α→1，表示無(wú)窮快的反應(yīng)；而pΔt→0時(shí)，α→1/2，表示無(wú)窮慢的反應(yīng)。αQSS方法中利用預(yù)測(cè)校正的方法獲得一定的精度并保證質(zhì)量分?jǐn)?shù)守恒，而在計(jì)算中需要不斷對(duì)校正步進(jìn)行迭代，直到滿足迭代標(biāo)準(zhǔn)，預(yù)測(cè)和校正步的詳細(xì)實(shí)施過(guò)程見(jiàn)文獻(xiàn)[10-11]。

使用點(diǎn)隱方法時(shí)，將源項(xiàng)寫(xiě)成n+1時(shí)間步的值，通過(guò)線性化消除源項(xiàng)剛性的問(wèn)題，實(shí)現(xiàn)方式為：

(11)

在與流動(dòng)耦合，使用具有TVD性質(zhì)的二階Runge-Kutta法推進(jìn)時(shí)，其形式可寫(xiě)為:

(12)

式中：I為單位質(zhì)量矩陣；R為空間方程離散后的余項(xiàng)。源項(xiàng)雅各比矩陣推導(dǎo)見(jiàn)文獻(xiàn)[2,16]。

2 源項(xiàng)處理方法的理論分析

2.1 源項(xiàng)處理方法的穩(wěn)定性分析

考慮模型方程：

(13)

式中：a為常數(shù)，s(u)為反應(yīng)源項(xiàng)。

為了簡(jiǎn)便，采用一階迎風(fēng)差分格式離散式(13)的對(duì)流項(xiàng)，得到如下形式:

(14)

式(14)中的源項(xiàng)可根據(jù)處理方法寫(xiě)為不同的形式。

考慮一步法，將式(14)改寫(xiě)為:

(15)

G=1-λΔt

(16)

其中，λ的具體形式為:

(17)

為滿足求解穩(wěn)定性，要求誤差放大因子G的模小于或等于1，解出滿足穩(wěn)定要求的Δt，其形式為：

(18)

從式(17)可以看出，當(dāng)τ(c)非常小時(shí)，其倒數(shù)遠(yuǎn)大于其他項(xiàng)的值，式(17)變?yōu)棣恕?/τ(c)。因此式(18)可簡(jiǎn)化為：

Δt≤2τ(c)

(19)

可見(jiàn)，在使用一步法求解剛性問(wèn)題嚴(yán)重的化學(xué)反應(yīng)流動(dòng)時(shí)，化學(xué)源項(xiàng)積分的時(shí)間步長(zhǎng)需要小于或等于兩倍最小特征反應(yīng)時(shí)間。

逼近法的分析類似，假設(shè)在Δt內(nèi)q和τ為常數(shù)，則可將式(14)改寫(xiě)為：

(20)

該形式與梯形公式相同，其穩(wěn)定性分析見(jiàn)文獻(xiàn)[1]，結(jié)果表明，逼近法只要時(shí)間步長(zhǎng)滿足對(duì)流項(xiàng)穩(wěn)定條件即可。

點(diǎn)隱方法的分析與一步法類似，只需將式(15)中源項(xiàng)的變量換為n+1時(shí)刻的值，解出誤差放大因子的模為:

(21)

式中：aΔt/Δx為CFL(Courant Friedrichs Lewy)數(shù)，而Δt/τ(c)為剛性系數(shù)，可知只要滿足對(duì)流項(xiàng)穩(wěn)定條件，即CFL數(shù)小于1，誤差放大因子的模是恒小于1的。

αQSS方法的穩(wěn)定性分析可見(jiàn)文獻(xiàn)[10]，在滿足對(duì)流項(xiàng)穩(wěn)定的條件下是絕對(duì)穩(wěn)定的。

穩(wěn)定性分析的結(jié)果表明:隱式方法在求解數(shù)學(xué)意義上的剛性問(wèn)題時(shí)具有卓越的性能；一步法為了滿足求解的穩(wěn)定性，需要很小的時(shí)間步長(zhǎng)；逼近法和αQSS方法中的特征消耗時(shí)間τi滿足在求解時(shí)間步長(zhǎng)內(nèi)為常數(shù)的情況下，也具有很好的求解性能。

2.2 對(duì)化學(xué)反應(yīng)流動(dòng)模擬的適用性分析

化學(xué)反應(yīng)流動(dòng)存在多個(gè)時(shí)間尺度，但可能不完全符合剛性問(wèn)題的數(shù)學(xué)定義[17]，因?yàn)榛瘜W(xué)反應(yīng)源項(xiàng)的雅各比矩陣特征值實(shí)部可能不全為負(fù)值[7]。此時(shí)，特征值為正且相當(dāng)大的情況表現(xiàn)為物理上的剛性問(wèn)題，求解時(shí)必須保證時(shí)間步長(zhǎng)足夠小以捕捉到物理特征的變化。如果源項(xiàng)雅各比矩陣特征值實(shí)部全為負(fù)值，則對(duì)應(yīng)的物理特征表現(xiàn)為組分快速變化到反應(yīng)平衡狀態(tài)。

基于以上分析，以點(diǎn)隱為首的隱式方法在求解數(shù)學(xué)意義上的剛性問(wèn)題時(shí)具有優(yōu)秀的性能，但針對(duì)物理意義上的剛性問(wèn)題時(shí)，由于點(diǎn)隱對(duì)化學(xué)反應(yīng)源項(xiàng)積分沒(méi)有時(shí)間步長(zhǎng)限制，可能會(huì)由于時(shí)間步長(zhǎng)大而導(dǎo)致分辨不出物理特征變化的細(xì)節(jié)；另外，流場(chǎng)求解時(shí)網(wǎng)格數(shù)較多，使用隱式方法時(shí)，在每個(gè)網(wǎng)格單元都要進(jìn)行矩陣的求逆運(yùn)算，運(yùn)算量隨變量個(gè)數(shù)二次方增加[1]，導(dǎo)致流場(chǎng)推進(jìn)速度很慢。一步法受穩(wěn)定性要求的嚴(yán)格限制，在剛性很?chē)?yán)重的時(shí)候，時(shí)間步長(zhǎng)取值造成的計(jì)算量需求將難以接受，但在求解物理和數(shù)學(xué)上都為剛性的問(wèn)題時(shí)還是具有一定的優(yōu)勢(shì)，至少能夠保證捕捉到變化細(xì)節(jié)所需要的時(shí)間步長(zhǎng)。逼近法和αQSS方法在求解時(shí)簡(jiǎn)化了物理問(wèn)題，將組分反應(yīng)源項(xiàng)考慮為解耦形式，從而降低了求解的難度，但造成了在計(jì)算時(shí)容易出現(xiàn)質(zhì)量不守恒的問(wèn)題，使用逼近法時(shí)需要不斷修正，而αQSS方法則利用組分迭代校正到一定精度來(lái)滿足質(zhì)量守恒。盡管如此，αQSS在積分化學(xué)反應(yīng)源項(xiàng)時(shí)還是具有良好的計(jì)算性能，且對(duì)于化學(xué)反應(yīng)的變化適應(yīng)性很好，具體分析如下。

將組分i的化學(xué)反應(yīng)凈生成項(xiàng)寫(xiě)為式(7)形式，并將式(5)寫(xiě)成式(9)形式，得到：

(22)

逼近法對(duì)應(yīng)αQSS方法中pΔt→0的情況，而梯形公式與逼近法本質(zhì)上相同，只是缺少校正步。

對(duì)式(7)分別采用隱式和顯式進(jìn)行數(shù)值求解，得到如下兩種形式:

(23)

(24)

分別將式(23)、(24)改寫(xiě)成式(9)形式，得到：

(25)

(26)

由αQSS中參數(shù)α的定義，式(25)為解耦的隱式方法，對(duì)應(yīng)pΔt→∞的情況，適合快速衰減的反應(yīng)，與點(diǎn)隱方法類似。而式(26)則為一步法，對(duì)應(yīng)pΔt→-∞的情況，表明組分消耗特征時(shí)間為負(fù)值，而αQSS方法不存在這種情況。因此，由式(22)～(26)可以看出，如果組分反應(yīng)可表示為解耦的形式，如式(7)，則αQSS方法幾乎涵蓋了其他幾種方法。但是，αQSS方法中的參數(shù)α變化只能表征組分衰減到平衡狀態(tài)快慢的速度，即對(duì)應(yīng)特征值小于零的情況，在生成項(xiàng)q遠(yuǎn)大于消耗項(xiàng)pρ時(shí)，只依賴參數(shù)α變化并不能捕捉到快速增長(zhǎng)的化學(xué)反應(yīng)特征，此時(shí)αQSS方法通過(guò)對(duì)比生成項(xiàng)和消耗項(xiàng)的大小，自動(dòng)調(diào)整時(shí)間步長(zhǎng)以適應(yīng)化學(xué)反應(yīng)的變化[11]，其性質(zhì)與一步法相同。從以上分析可以得到，對(duì)應(yīng)物理模態(tài)快速衰減的剛性問(wèn)題時(shí)，αQSS方法通過(guò)參數(shù)α的變化來(lái)適應(yīng)不同的衰減速度；而對(duì)于物理模態(tài)快速增長(zhǎng)的剛性問(wèn)題時(shí)，αQSS方法通過(guò)改變積分時(shí)間步長(zhǎng)來(lái)適應(yīng)增長(zhǎng)速度，因此對(duì)化學(xué)反應(yīng)變化的適應(yīng)性強(qiáng)。

3 數(shù)值計(jì)算

Lehr[18]進(jìn)行了一系列彈道靶實(shí)驗(yàn)，將球形彈丸以接近爆震的速度射入一定條件下的氫氣和空氣預(yù)混氣體中，觀察到彈丸頭部不穩(wěn)定的燃燒現(xiàn)象，并拍攝了不同入射速度下彈丸頭部振蕩燃燒現(xiàn)象的陰影圖，圖1(a)、(b)分別對(duì)應(yīng)來(lái)流馬赫數(shù)為4.48和4.79的工況。McVey等[19]提出了解釋激波誘導(dǎo)振蕩燃燒的機(jī)理，其后Choi[20]對(duì)該現(xiàn)象做了詳盡的數(shù)值模擬，認(rèn)為激波誘導(dǎo)周期振蕩燃燒現(xiàn)象是驗(yàn)證非平衡化學(xué)流計(jì)算程序的典型算例。本文中使用該算例，從計(jì)算結(jié)果、計(jì)算效率等方面比較了源項(xiàng)處理方法的性能。

3.1 計(jì)算結(jié)果可靠性分析

計(jì)算模型和網(wǎng)格選取與文獻(xiàn)[21]相同，取球頭前方2.5 mm，上方10 mm，網(wǎng)格點(diǎn)數(shù)為251×201。計(jì)算工況與Lehr[18]實(shí)驗(yàn)條件一致，來(lái)流溫度293 K，壓力42 665 Pa，預(yù)混氣體來(lái)流速度分別為1 804 m/s(Ma=4.48) 和1 931 m/s (Ma=4.79)，該實(shí)驗(yàn)條件下對(duì)應(yīng)的爆震馬赫數(shù)約為4.845，CFL數(shù)取為0.85。通過(guò)記錄球頭軸線駐點(diǎn)處壓力隨時(shí)間變化的數(shù)據(jù)，對(duì)其進(jìn)行快速傅里葉變換得到振蕩頻率，對(duì)比計(jì)算振蕩頻率和實(shí)驗(yàn)所得振蕩頻率，驗(yàn)證化學(xué)反應(yīng)模塊的正確性。為了避免由于取值區(qū)間小而造成的振蕩頻率計(jì)算的波動(dòng)，用于計(jì)算主頻的壓力時(shí)間曲線選取至少100個(gè)振蕩周期。

為驗(yàn)證源項(xiàng)處理方法計(jì)算結(jié)果的正確性，對(duì)來(lái)流馬赫數(shù)分別為4.48和4.79的工況進(jìn)行了計(jì)算，其軸線駐點(diǎn)壓力振蕩頻率見(jiàn)表1。從表1可以看出，來(lái)流馬赫數(shù)分別為4.48和4.79時(shí)，不同源項(xiàng)處理方法得到的振蕩頻率與實(shí)驗(yàn)值相對(duì)誤差都較小。在使用251×201的網(wǎng)格尺度且來(lái)流馬赫數(shù)4.79時(shí)的計(jì)算條件下，計(jì)算得到的組分最小反應(yīng)特征時(shí)間約為6×10-11s，流動(dòng)特征時(shí)間為3.5×10-9s，剛性系數(shù)Rs約為58，從穩(wěn)定性分析可知，采用的化學(xué)反應(yīng)時(shí)間步長(zhǎng)應(yīng)小于等于2倍特征反應(yīng)時(shí)間，但文獻(xiàn)[2]表明，在實(shí)際計(jì)算中化學(xué)反應(yīng)推進(jìn)步長(zhǎng)可取為10倍最小反應(yīng)特征時(shí)間。表1所列的一步法中，化學(xué)反應(yīng)時(shí)間步長(zhǎng)為6倍最小時(shí)間尺度，計(jì)算得到的振蕩頻率與實(shí)驗(yàn)值吻合，但從穩(wěn)定性方面考慮，計(jì)算過(guò)程中可能存在一定的振蕩。在Jachimowski化學(xué)反應(yīng)模型中，組分H2O2的化學(xué)反應(yīng)特征時(shí)間最小，因此出現(xiàn)振蕩的可能性最大，圖2～3所示為對(duì)應(yīng)表1中Ma=4.79時(shí)，不同源項(xiàng)處理方法計(jì)算得到的球頭軸線上組分變化曲線，其中圖2為不同積分時(shí)間步長(zhǎng)的一步法計(jì)算結(jié)果，圖3(a)、(b)、(c)分別為αQSS、點(diǎn)隱和逼近法的計(jì)算結(jié)果。由于各組分之間的質(zhì)量分?jǐn)?shù)相差較大，為了對(duì)比的直觀性，對(duì)圖中顯示的各組分質(zhì)量分?jǐn)?shù)進(jìn)行歸一化處理，歸一化之后，H2的歸一化質(zhì)量分?jǐn)?shù)為0.03，H2O的歸一化質(zhì)量分?jǐn)?shù)采用0.25，H2O2的歸一化質(zhì)量分?jǐn)?shù)為0.000 5，而N2的歸一化質(zhì)量分?jǐn)?shù)為0.76。

表1 實(shí)驗(yàn)數(shù)據(jù)與理論模型結(jié)果對(duì)比Table 1 Comparison of experimental and theoretical data

圖2(a)表示使用一步法時(shí)積分步長(zhǎng)為6倍最小化學(xué)反應(yīng)特征時(shí)間的計(jì)算結(jié)果，與圖2(b)使用兩倍化學(xué)反應(yīng)特征時(shí)間的計(jì)算結(jié)果對(duì)比可發(fā)現(xiàn)，在燃燒陣面之后，H2O2變化曲線存在非物理振蕩現(xiàn)象，與穩(wěn)定性分析的預(yù)測(cè)結(jié)果一致，即積分時(shí)間步長(zhǎng)超過(guò)穩(wěn)定性要求的最大時(shí)間步長(zhǎng)后，計(jì)算結(jié)果將出現(xiàn)不穩(wěn)定或者發(fā)散現(xiàn)象，說(shuō)明在使用一步法積分化學(xué)反應(yīng)剛性源項(xiàng)時(shí)，要想獲得準(zhǔn)確的組分質(zhì)量分?jǐn)?shù)變化曲線，時(shí)間步長(zhǎng)必須滿足穩(wěn)定性要求。積分時(shí)間步長(zhǎng)取值滿足穩(wěn)定要求時(shí)，幾種源項(xiàng)處理方法的結(jié)果相同，如圖3(a)～(c)所示，因此，在計(jì)算化學(xué)反應(yīng)流動(dòng)時(shí)，逼近法、αQSS和點(diǎn)隱在穩(wěn)定性方面要優(yōu)于一步法。

觀察圖2(a)，發(fā)現(xiàn)振蕩并沒(méi)有隨時(shí)間和空間變化而加劇，這是因?yàn)榛瘜W(xué)反應(yīng)具有強(qiáng)非線性與隨時(shí)間變化的特點(diǎn)，源項(xiàng)的雅各比矩陣特征值將在反應(yīng)過(guò)程中不斷變化，使得引入的小擾動(dòng)在增長(zhǎng)模態(tài)和衰

減模態(tài)之間反復(fù)轉(zhuǎn)變，因此擾動(dòng)不會(huì)進(jìn)一步增強(qiáng)。另外，從表1可知，在使用一步法計(jì)算時(shí)，積分時(shí)間步長(zhǎng)不滿足穩(wěn)定性要求時(shí)，也能獲得正確的燃燒振蕩頻率。這是由于振蕩頻率主要由化學(xué)反應(yīng)放熱量決定的，H2O2組分質(zhì)量分?jǐn)?shù)很小，對(duì)于整個(gè)化學(xué)反應(yīng)放熱量的影響很微弱，因此其質(zhì)量分?jǐn)?shù)的小波動(dòng)對(duì)于振蕩頻率幾乎沒(méi)有影響。此外，圖3中組分變化曲線表明這幾種源項(xiàng)處理方法對(duì)本算例中化學(xué)反應(yīng)特征的捕捉效果基本相同，其原因在于化學(xué)反應(yīng)是一個(gè)動(dòng)態(tài)的非線性系統(tǒng)，某組分在這一時(shí)刻呈快速的增長(zhǎng)模態(tài)，在下一時(shí)刻變?yōu)榭焖偎p的模態(tài)或者趨于平衡；并且本算例的剛性系數(shù)不大，化學(xué)反應(yīng)特征在單個(gè)流動(dòng)時(shí)間步內(nèi)的變化也很小，因此各個(gè)源項(xiàng)處理方法對(duì)于該算例的捕捉效果基本相同。

3.2 計(jì)算效率分析

在計(jì)算硬件條件相同且CFL數(shù)為0.85的情況下，以來(lái)流馬赫數(shù)4.79且完全發(fā)展的振蕩燃燒流場(chǎng)作為初始計(jì)算流場(chǎng)，比較幾種源項(xiàng)處理方法的計(jì)算效率。表2所示為流場(chǎng)推進(jìn)物理時(shí)間473 μs時(shí)，不同源項(xiàng)處理方法所用的CPU時(shí)間，其中一步法的化學(xué)推進(jìn)時(shí)間步長(zhǎng)為6倍最小反應(yīng)特征時(shí)間。

表2 不同源項(xiàng)處理方法消耗CPU時(shí)間Table 2 The CPU time in differentstiff source term methods

從表2中可以看到，點(diǎn)隱與一步法消耗的CPU時(shí)間最多，采用源項(xiàng)完全求逆的點(diǎn)隱所需時(shí)間約為αQSS的2倍，點(diǎn)隱方法需要計(jì)算源項(xiàng)的雅各比矩陣，并涉及復(fù)雜的矩陣求逆運(yùn)算，求解時(shí)間隨組分的增加呈指數(shù)增長(zhǎng)，在使用基元反應(yīng)模型求解爆震燃燒問(wèn)題時(shí)，計(jì)算量很大；而使用一步法時(shí)，時(shí)間步長(zhǎng)需取足夠小以滿足穩(wěn)定性要求，極大降低了計(jì)算效率；αQSS計(jì)算時(shí)間與組分求解迭代時(shí)指定的誤差ε有關(guān)，ε越小，所需計(jì)算時(shí)間越長(zhǎng)；逼近法相對(duì)αQSS來(lái)說(shuō)，缺少迭代過(guò)程，因此所用時(shí)間最少。由此可見(jiàn)，在計(jì)算結(jié)果精度基本相同的情況下，逼近法和αQSS的計(jì)算效率較其他兩種處理方法高。

4 結(jié) 論

本文中從穩(wěn)定性分析和數(shù)值計(jì)算方面對(duì)比一步法、逼近法、αQSS和點(diǎn)隱方法計(jì)算化學(xué)反應(yīng)剛性源項(xiàng)問(wèn)題的性能，得到結(jié)論如下：(1)一步法在積分剛性源項(xiàng)時(shí)，積分步長(zhǎng)需小于或等于2倍最小時(shí)間尺度，而其他3種源項(xiàng)計(jì)算方法對(duì)時(shí)間步長(zhǎng)取值沒(méi)有影響；(2)一步法、逼近法實(shí)質(zhì)上是αQSS的特例，αQSS方法可根據(jù)化學(xué)反應(yīng)特征自動(dòng)調(diào)整系數(shù)α和時(shí)間步長(zhǎng)，適用范圍廣；(3)一步法在積分化學(xué)反應(yīng)剛性源項(xiàng)時(shí)，可適當(dāng)放寬穩(wěn)定性要求的時(shí)間步長(zhǎng)限制；(4)在計(jì)算結(jié)果相近以及和精度基本相同的前提下，逼近法和αQSS在求解剛性源項(xiàng)時(shí)具有較高的效率，并且αQSS能夠適應(yīng)化學(xué)反應(yīng)的變化，是采用分裂法求解剛性源項(xiàng)時(shí)綜合性能最好的方法。

[1] 劉君,周松柏,徐春光.超聲速流動(dòng)中燃燒現(xiàn)象的數(shù)值模擬方法及應(yīng)用[M].長(zhǎng)沙:國(guó)防科技大學(xué)出版社,2008:76-79.

[2] BUSSING T R A. A finite volume method for the Navier-Stokes equations with finite rate chemistry[D]. Cambridge: Massachusetts Institute of Technology, 1985.

[3] BUSSING T R A, Murman E M. Finite-volume method for the calculation of compressible chemically reacting flows[J]. AIAA Journal, 1988,26(9):1070-1078.

[4] ZHONG X. Additive semi-implicit Runge-Kutta methods for computing high-speed non-equilibrium reactive flows[J]. Journal of Computational Physics, 1996,128(1):19-31.

[5] ROSENBROCK H H. Some general implicit processes for the numerical solution of differential equations[J]. The Computer Journal, 1963,5(4):329-330.

[6] LEVEQUE R J, YEE H C. A study of numerical methods for hyperbolic conservation laws with stiff source terms[J]. Journal of Computational Physics, 1990,86(1):187-210.

[7] ORAN E S, BORIS J P. Numerical simulation of reactive flow[M]. Cambridge: Cambridge University Press, 2005:114-158.

[8] YOUNG T R, BORIS J P. A numerical technique for solving stiff ordinary differential equations associated with the chemical kinetics of reactive-flow problems[J]. The Journal of Physical Chemistry, 1977,81(25):2424-2427.

[9] CHIANG T, HOFFMANN K. Determination of computational time step for chemically reacting flows[C]∥Proceedings of AIAA 20th Fluid Dynamics, Plasma Dynamics and Laser Conference. Buffalo, New York, USA, 1989.

[10] MOTT D R, ORAN E S, VAN L B. A quasi-steady-state solver for the stiff ordinary differential equations of reaction kinetics[J]. bJournal of Computational Physics, 2000,164(2):407-428.

[11] MOTT D R, ORAN E S. CHEMEQ2: A solver for the stiff ordinary differential equations of chemical kinetics[R]. Naval Research Lab, Washington D C, 2001.

[12] 劉瑜.化學(xué)非平衡流的計(jì)算方法研究及其在激波誘導(dǎo)燃燒現(xiàn)象模擬中的應(yīng)用[D].長(zhǎng)沙:國(guó)防科技大學(xué),2008.

LIU Yu. Investigations into numerical methods of chemical non-equilibrium flow and its application to simulation of shock-induced combustion phenomenon[D]. Changsha: National University of Defense Technology, 2008.

[13] 劉君,張涵信,高樹(shù)椿.一種新型的計(jì)算化學(xué)非平衡流動(dòng)的解耦方法[J].國(guó)防科技大學(xué)學(xué)報(bào),2000,22(5):19-23.

LIU Jun, ZHANG Hanxin, GAO Shuchun. A new uncoupled method for numerical simulation of non-equilibrium flow[J]. Journal of National University of Defense Technology, 2000,22(5):19-23.

[14] 劉世杰,林志勇,孫明波,等.采用不同化學(xué)反應(yīng)源項(xiàng)處理方法的胞格爆轟數(shù)值研究[J].國(guó)防科技大學(xué)學(xué)報(bào),2010,32(5):1-6.

LIU Shijie, LIN Zhiyong, SUN Mingbo, et al. Numerical simulation of cellular detonation using different chemical reaction source term methods[J]. Journal of National University of Defense Technology, 2010,32(5):1-6.

[15] JACHIMOWSKI C J. An analytical study of the hydrogen-air reaction mechanism with application to scramjet combustion[R]. NASA TechnicaI Paper 2791, 1988.

[16] 潘沙.高超聲速氣動(dòng)熱數(shù)值模擬方法及大規(guī)模并行計(jì)算研究[D].長(zhǎng)沙:國(guó)防科技大學(xué),2010.

PAN Sha. Hypersonic aerothermal numerical simulation method and massive parallel computation research[D]. Changsha: National University of Defense Technology, 2010.

[17] TORO E F. Riemann solvers and numerical methods for fluid dynamics: A practical introduction[M]. Springer Science & Business Media, 2009:531-542.

[18] LEHR H F. Experiments on shock-induced combustion[J]. Astronautica Acta, 1972,17:589-597.

[19] MCVEY J B, TOONG T Y. Mechanism of instabilities of exothermic hypersonic blunt-body flows[J]. Combustion Science and Technology, 1971,3(2):63-76.

[20] CHOI J Y, JEUNG I S, YOON Y. Computational fluid dynamics algorithms for unsteady shock-induced combustion: Part 1: Validation[J]. AIAA Journal, 2000,38(7):1179-1187.

[21] 劉世杰,孫明波,林志勇,等.鈍頭體激波誘導(dǎo)振蕩燃燒現(xiàn)象的數(shù)值模擬[J].力學(xué)學(xué)報(bào),2010,42(4):597-606.

LIU Shijie J, SUN Mingbo, LIN Zhiyong, et al. Numerical research on blunt body shock-induced oscillation combustion phenomenon[J]. Chinese Journal of Theoretical and Applied Mechanics, 2010,42(4):597-606.