剛性源項處理方法在爆震模擬中的應用對比*

鄧 利，馬 虎，武曉松，周長省

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

爆震波是激波和燃燒波的強烈耦合，涉及爆震燃燒的流動存在多個時間尺度，容易出現嚴重的剛性源項，使數值求解流動方程存在困難。

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

Leveque等[6]利用時間分裂法和MacCormack預測校正方法求解了含剛性源項的雙曲型方程，認為分裂算法在空間分辨率足夠的情況下，計算效果比耦合算法稍好。Oran等[7]對剛性源項方程的積分方法做了全面的論述。一步法雖然在求解剛性方程組時穩定性要求嚴格，但由于只需已知變量值，求解相對方便，并針對特征值為正且相當大的物理剛性問題時，使用一步法是比較有利的。Young等[8]采用逼近法與預測校正思想來處理化學反應源項，對剛性和普通常微分方程組都具有較好的計算性能，且方便與時間步長估算方法[9]配合使用，求解效率高，但容易出現質量分數不守恒的問題，需要在流場推進時不斷修正。Mott等[10-11]在逼近法的基礎上發展了αQSS方法，通過迭代校正組分值來保證質量守恒，比逼近法更加穩定和精確，并針對碳氫燃料的燃燒和爆震模擬取得了令人滿意的效果。

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

1 計算方法介紹

1.1 控制方程

激波誘導燃燒為激波和燃燒波耦合的結果，可忽略黏性的影響，因此使用有限速率化學反應的歐拉方程描述為:

(1)

1.2 數值方法

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

(2)

式中：算符Lf為流動求解算子，Lr為化學反應的求解算子，Δt為流場推進的時間步長，上標n為求解時間步數。

一步法為常微分方程初值問題的經典求解方法，求解變量n+1時刻的值時，只需n時刻的信息。解耦后的化學源項方程如下：

(3)

由于Strang分裂對算子求解精度的要求，選取2階Runge-Kutta法作為化學反應求解算子。

在利用逼近法求解化學反應源項時，將組分i的凈生成率寫成如下形式：

(4)

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

(5)

(6)

αQSS方法在逼近法的基礎上進一步改進其迭代形式，在求解時將組分凈生成率寫為：

(7)

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

(8)

整理式(8)得到：

(9)

(10)

式中：pΔt→∞時，α→1，表示無窮快的反應；而pΔt→0時，α→1/2，表示無窮慢的反應。αQSS方法中利用預測校正的方法獲得一定的精度并保證質量分數守恒，而在計算中需要不斷對校正步進行迭代，直到滿足迭代標準，預測和校正步的詳細實施過程見文獻[10-11]。

使用點隱方法時，將源項寫成n+1時間步的值，通過線性化消除源項剛性的問題，實現方式為：

(11)

在與流動耦合，使用具有TVD性質的二階Runge-Kutta法推進時，其形式可寫為:

(12)

式中：I為單位質量矩陣；R為空間方程離散后的余項。源項雅各比矩陣推導見文獻[2,16]。

2 源項處理方法的理論分析

2.1 源項處理方法的穩定性分析

考慮模型方程：

(13)

式中：a為常數，s(u)為反應源項。

為了簡便，采用一階迎風差分格式離散式(13)的對流項，得到如下形式:

(14)

式(14)中的源項可根據處理方法寫為不同的形式。

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

(15)

G=1-λΔt

(16)

其中，λ的具體形式為:

(17)

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

(18)

從式(17)可以看出，當τ(c)非常小時，其倒數遠大于其他項的值，式(17)變為λ≈1/τ(c)。因此式(18)可簡化為：

Δt≤2τ(c)

(19)

可見，在使用一步法求解剛性問題嚴重的化學反應流動時，化學源項積分的時間步長需要小于或等于兩倍最小特征反應時間。

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

(20)

該形式與梯形公式相同，其穩定性分析見文獻[1]，結果表明，逼近法只要時間步長滿足對流項穩定條件即可。

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

(21)

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

αQSS方法的穩定性分析可見文獻[10]，在滿足對流項穩定的條件下是絕對穩定的。

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

2.2 對化學反應流動模擬的適用性分析

化學反應流動存在多個時間尺度，但可能不完全符合剛性問題的數學定義[17]，因為化學反應源項的雅各比矩陣特征值實部可能不全為負值[7]。此時，特征值為正且相當大的情況表現為物理上的剛性問題，求解時必須保證時間步長足夠小以捕捉到物理特征的變化。如果源項雅各比矩陣特征值實部全為負值，則對應的物理特征表現為組分快速變化到反應平衡狀態。

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

將組分i的化學反應凈生成項寫為式(7)形式，并將式(5)寫成式(9)形式，得到：

(22)

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

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

(23)

(24)

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

(25)

(26)

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

3 數值計算

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

3.1 計算結果可靠性分析

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

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

表1 實驗數據與理論模型結果對比Table 1 Comparison of experimental and theoretical data

圖2(a)表示使用一步法時積分步長為6倍最小化學反應特征時間的計算結果，與圖2(b)使用兩倍化學反應特征時間的計算結果對比可發現，在燃燒陣面之后，H2O2變化曲線存在非物理振蕩現象，與穩定性分析的預測結果一致，即積分時間步長超過穩定性要求的最大時間步長后，計算結果將出現不穩定或者發散現象，說明在使用一步法積分化學反應剛性源項時，要想獲得準確的組分質量分數變化曲線，時間步長必須滿足穩定性要求。積分時間步長取值滿足穩定要求時，幾種源項處理方法的結果相同，如圖3(a)～(c)所示，因此，在計算化學反應流動時，逼近法、αQSS和點隱在穩定性方面要優于一步法。

觀察圖2(a)，發現振蕩并沒有隨時間和空間變化而加劇，這是因為化學反應具有強非線性與隨時間變化的特點，源項的雅各比矩陣特征值將在反應過程中不斷變化，使得引入的小擾動在增長模態和衰

減模態之間反復轉變，因此擾動不會進一步增強。另外，從表1可知，在使用一步法計算時，積分時間步長不滿足穩定性要求時，也能獲得正確的燃燒振蕩頻率。這是由于振蕩頻率主要由化學反應放熱量決定的，H2O2組分質量分數很小，對于整個化學反應放熱量的影響很微弱，因此其質量分數的小波動對于振蕩頻率幾乎沒有影響。此外，圖3中組分變化曲線表明這幾種源項處理方法對本算例中化學反應特征的捕捉效果基本相同，其原因在于化學反應是一個動態的非線性系統，某組分在這一時刻呈快速的增長模態，在下一時刻變為快速衰減的模態或者趨于平衡；并且本算例的剛性系數不大，化學反應特征在單個流動時間步內的變化也很小，因此各個源項處理方法對于該算例的捕捉效果基本相同。

3.2 計算效率分析

在計算硬件條件相同且CFL數為0.85的情況下，以來流馬赫數4.79且完全發展的振蕩燃燒流場作為初始計算流場，比較幾種源項處理方法的計算效率。表2所示為流場推進物理時間473 μs時，不同源項處理方法所用的CPU時間，其中一步法的化學推進時間步長為6倍最小反應特征時間。

表2 不同源項處理方法消耗CPU時間Table 2 The CPU time in differentstiff source term methods

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

4 結 論

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

[1] 劉君,周松柏,徐春光.超聲速流動中燃燒現象的數值模擬方法及應用[M].長沙:國防科技大學出版社,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] 劉瑜.化學非平衡流的計算方法研究及其在激波誘導燃燒現象模擬中的應用[D].長沙:國防科技大學,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] 劉君,張涵信,高樹椿.一種新型的計算化學非平衡流動的解耦方法[J].國防科技大學學報,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] 劉世杰,林志勇,孫明波,等.采用不同化學反應源項處理方法的胞格爆轟數值研究[J].國防科技大學學報,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].長沙:國防科技大學,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] 劉世杰,孫明波,林志勇,等.鈍頭體激波誘導振蕩燃燒現象的數值模擬[J].力學學報,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.