劉君S，董海波，劉瑜

1.大連理工大學(xué) 航空航天學(xué)院，大連 116024 2.航天工程大學(xué) 航天裝備系，北京 101416

化學(xué)非平衡流動(dòng)與量熱完全氣體相比，控制方程多出了化學(xué)反應(yīng)源項(xiàng)和組元質(zhì)量守恒方程組。由于化學(xué)反應(yīng)速率對(duì)溫度、密度等熱力學(xué)狀態(tài)函數(shù)極其敏感，導(dǎo)致不同基元反應(yīng)和同一基元反應(yīng)的正逆向反應(yīng)之間的特征時(shí)間尺度相差較大，計(jì)算過(guò)程中存在剛性問(wèn)題。在20世紀(jì)70年代左右，為解決航天器再入過(guò)程的熱障和黑障問(wèn)題，需要模擬空氣介質(zhì)的化學(xué)非平衡流動(dòng)，為此提出了點(diǎn)隱算法和全隱算法(統(tǒng)稱(chēng)隱式算法)，目前依然是處理剛性問(wèn)題的主要算法?，F(xiàn)代精細(xì)化學(xué)研究得到的反應(yīng)機(jī)理越來(lái)越復(fù)雜，對(duì)于碳?xì)淙剂?，常包含成百上千種組元和基元反應(yīng)。復(fù)雜反應(yīng)機(jī)理引入的組元方程數(shù)目遠(yuǎn)大于三維量熱完全氣體流動(dòng)方程的5個(gè)求解變量，例如對(duì)于生物柴油燃燒過(guò)程，文獻(xiàn)[1]給出了1 034種組元、4 236個(gè)基元反應(yīng)的機(jī)理模型，采用現(xiàn)有計(jì)算方法準(zhǔn)確模擬這種化學(xué)非平衡流動(dòng)，控制方程至少需要1 038個(gè)偏微分方程構(gòu)成方程組。采用隱式算法求解剛性問(wèn)題涉及到矩陣求逆，由此帶來(lái)的巨大計(jì)算量限制了模擬手段在實(shí)際三維工程問(wèn)題中的應(yīng)用。

為減少求解變量，國(guó)內(nèi)外提出很多算法。早期采用忽略鏈?zhǔn)椒磻?yīng)中間過(guò)程的總包反應(yīng)，例如，1997年梁劍寒和王承堯[2]采用文獻(xiàn)[3]提出的兩步總包反應(yīng)模擬超燃沖壓發(fā)動(dòng)機(jī)內(nèi)部的氫氧燃燒過(guò)程。隨著質(zhì)譜儀等測(cè)試手段的進(jìn)步，對(duì)碳?xì)淙剂戏磻?yīng)機(jī)理研究逐步深入，在此基礎(chǔ)上提出反應(yīng)機(jī)理較為簡(jiǎn)明的替代反應(yīng)模型，例如針對(duì)煤油，黃生洪等[4]對(duì)國(guó)外提出的兩套C12H23替代模型進(jìn)行考核，修正為17種組元、30個(gè)基元反應(yīng)；楊順華等[5]根據(jù)天津大學(xué)測(cè)量結(jié)果提出10種組元、12個(gè)基元反應(yīng)的n-C10H22替代模型。不同替代反應(yīng)模型對(duì)燃燒最終產(chǎn)物影響較小，但是低溫段存在明顯差異，對(duì)于煤油燃燒這樣復(fù)雜的化學(xué)反應(yīng)，尋找高低溫段特性接近的替代反應(yīng)模型非常困難。目前，超燃沖壓發(fā)動(dòng)機(jī)的點(diǎn)火延遲過(guò)程和內(nèi)燃機(jī)減少NOx排放的低溫燃燒技術(shù)作為化學(xué)非平衡流動(dòng)模擬的主要研究背景，新的發(fā)展趨勢(shì)是基于原始復(fù)雜模型采用系統(tǒng)辨識(shí)、優(yōu)化算法等數(shù)學(xué)方法，在保持燃燒過(guò)程主要特征的條件下盡可能合并基元反應(yīng)從而減少組元數(shù)目，這一方法稱(chēng)之為骨架反應(yīng)機(jī)理[1]。對(duì)上述生物柴油，文獻(xiàn)[1]研究結(jié)論是可以簡(jiǎn)化到49種組元、163個(gè)基元反應(yīng)，即使這樣采用點(diǎn)隱算法或全隱算法的計(jì)算量仍然巨大，為此，近些年部分研究者開(kāi)始重新審視解耦算法。

解耦算法又稱(chēng)為時(shí)間分裂法、分?jǐn)?shù)步法[6]。解耦算法將復(fù)雜的多物理問(wèn)題按照物理意義分解為若干個(gè)子問(wèn)題，然后分別對(duì)子問(wèn)題進(jìn)行求解。分解后的子問(wèn)題具有各自的物理特性，針對(duì)子問(wèn)題的計(jì)算，通?？梢圆捎孟鄳?yīng)的成熟算法進(jìn)行求解，這樣相對(duì)于完全耦合算法而言，便可以減少問(wèn)題求解的難度；另外，復(fù)雜問(wèn)題往往存在跨度很大的時(shí)間或空間尺度，通過(guò)解耦算法分解得到的問(wèn)題包含相同的物理尺度，因而對(duì)于物理尺度的處理會(huì)更加靈活，從而可以減少計(jì)算時(shí)間。由于解耦算法具有以上優(yōu)點(diǎn)，其應(yīng)用范圍十分廣泛，比如生物學(xué)[7-8]、金融學(xué)[9]、物理學(xué)[10-11]等眾多領(lǐng)域都有涉及。由于化學(xué)反應(yīng)流動(dòng)是燃燒過(guò)程和流動(dòng)過(guò)程相互耦合的復(fù)雜物理現(xiàn)象，是典型的非線性多尺度問(wèn)題，采用解耦算法對(duì)化學(xué)非平衡流動(dòng)進(jìn)行模擬是一個(gè)重要的研究手段。解耦算法應(yīng)用于化學(xué)非平衡流動(dòng)的模擬，即通過(guò)在空間點(diǎn)上把控制方程分解為流動(dòng)和化學(xué)反應(yīng)兩部分，各自采用適合的算法。實(shí)際上最早的化學(xué)非平衡流動(dòng)模擬就是采用解耦算法[12]，源項(xiàng)計(jì)算采用顯式格式，需要上千計(jì)算步才能與流動(dòng)時(shí)間推進(jìn)步長(zhǎng)相匹配，由于計(jì)算效率低下且不穩(wěn)定而逐漸被淘汰。雖然采用解耦算法可以大大降低化學(xué)反應(yīng)流動(dòng)數(shù)值模擬的難度，但是由于熱物性參數(shù)顯著依賴(lài)于溫度的變化，因而在流動(dòng)計(jì)算中，需要大量的工作才能將已有算法移植到化學(xué)非平衡流模擬中，而且此過(guò)程往往較為繁瑣，特別是對(duì)于復(fù)雜的高階精度格式更是如此。為了解決這一問(wèn)題，劉君在其博士論文中提出了一種新穎的非平衡流動(dòng)解耦算法[13]，通過(guò)對(duì)能量方程進(jìn)行分解，將適用于量熱完全氣體的流動(dòng)計(jì)算格式直接推廣到化學(xué)非平衡流動(dòng)的模擬，提高算法改進(jìn)的效率。隨后這種解耦算法不斷完善[14-22]，在化學(xué)反應(yīng)求解算法上，從最初采用穩(wěn)定性較差的梯形算法[19]，替換為穩(wěn)定性較好的剛性常微分變系數(shù)方法(VODE)和擬穩(wěn)態(tài)逼近方法(α-QSS)求解器[20]；汪洪波等[23]將這種解耦算法推廣到雙時(shí)間步法；針對(duì)不同網(wǎng)格結(jié)構(gòu)，首先在結(jié)構(gòu)網(wǎng)格上實(shí)現(xiàn)了解耦算法[13-14,24]，隨后為了適用于復(fù)雜計(jì)算域，又推廣到非結(jié)構(gòu)網(wǎng)格[20-21]，空間離散的精度從2階推廣到了5階[22]。該解耦算法也被應(yīng)用于爆轟波模擬[25-27]、固體火箭發(fā)動(dòng)機(jī)尾噴流動(dòng)模擬[24]、橫向噴流模擬[28-29]、沖壓發(fā)動(dòng)機(jī)流場(chǎng)數(shù)值模擬[18,24]、旋轉(zhuǎn)爆震波的模擬[30]和彈箭底排裝置尾部化學(xué)非平衡流模擬[31]等。經(jīng)過(guò)20多年的理論研究和實(shí)踐應(yīng)用，驗(yàn)證該解耦算法是成功的，為了區(qū)別于傳統(tǒng)的化學(xué)非平衡流動(dòng)解耦算法，本文將其命名為“劉君化學(xué)非平衡流動(dòng)解耦算法”，簡(jiǎn)稱(chēng)為“劉君算法”。

在非平衡流的數(shù)值模擬中化學(xué)反應(yīng)的求解工作量巨大，采用常規(guī)的求解器，化學(xué)反應(yīng)常微分方程組的計(jì)算量為O(N2)～O(N3),其中N為反應(yīng)組元數(shù)目。對(duì)于包含上千組元的詳細(xì)反應(yīng)機(jī)理，化學(xué)反應(yīng)的計(jì)算占據(jù)整個(gè)計(jì)算量的絕大部分。為此需要發(fā)展高效的化學(xué)反應(yīng)模擬方法。Ren等[41]通過(guò)動(dòng)態(tài)自適應(yīng)化學(xué)方法，在每一計(jì)算步得到詳細(xì)的骨架機(jī)理，在整個(gè)計(jì)算精度未受影響的情況下，計(jì)算效率提高了3～5倍。Gou等[47]提出了動(dòng)態(tài)多時(shí)間尺度(Multi-TimeScale, MTS)方法和混合多時(shí)間尺度(Hybrid Multi-TimeScale, HMTS)方法來(lái)提高化學(xué)反應(yīng)的計(jì)算效率，Sun等[48]將Gou等的動(dòng)態(tài)多時(shí)間尺度方法與動(dòng)態(tài)自適應(yīng)化學(xué)方法結(jié)合起來(lái)形成HMTS/CO-DAC(Hybrid Multi-TimeScale/COrrelated Dynamic Adaptive Chemistry)方法，顯著提高了計(jì)算效率。值得關(guān)注的是目前出現(xiàn)了一種非解耦算法，自適應(yīng)混合積分(Adaptive Hybrid Integration, AHI)方法[49]及其最新改進(jìn)型算法(AHI-S)[50]。這種算法在每一步識(shí)別出快反應(yīng)組元和慢反應(yīng)組元，對(duì)于快反應(yīng)組元采用隱式方法求解，對(duì)于慢反應(yīng)組元采用顯式方法求解。這種方法可將化學(xué)反應(yīng)求解的計(jì)算量降低到O(N)量級(jí)。

1 劉君解耦算法

1.1 出發(fā)方程

包含有n種組元的非平衡流守恒型Euler方程為

(1)

式中：F、G和H分別為對(duì)流通量在x、y和z方向的分量，具體表達(dá)形式參見(jiàn)文獻(xiàn)[51]，這里僅給出后面推導(dǎo)用到的守恒變量和源項(xiàng)表達(dá)式，即

(2)

(3)

E=0.5ρ(u2+v2+w2)+

(4)

第i個(gè)組元定壓比熱采用的擬合公式為

cpi=Ri(ai+biT+ciT2+diT3+eiT4)

(5)

式中：Ri為組元的氣體常數(shù)；系數(shù)ai、bi、ci、di、ei為常數(shù)，具體取值參見(jiàn)文獻(xiàn)[52]。

熱焓的表達(dá)式為

(6)

(7)

(8)

然后用E′代替E,并且對(duì)式(1)進(jìn)行修改得到新的守恒變量和化學(xué)反應(yīng)源項(xiàng)表達(dá)式為

U0=[ρρuρvρwE′ρ1ρ2…ρn-1]T

(9)

(10)

為了說(shuō)明方程重構(gòu)所表達(dá)的物理意義，進(jìn)一步引入符號(hào)γ′,把E′整理為

(11)

(12)

式中：等效比熱比γ′為溫度的函數(shù)，計(jì)算過(guò)程中采用二分法迭代得到。

從表達(dá)形式上看，E′與量熱完全氣體的總內(nèi)能完全相同，稱(chēng)為等效內(nèi)能；γ′與量熱完全氣體的比熱比相同，稱(chēng)為等效比熱比。E′和γ′不符合氣體動(dòng)力學(xué)中化學(xué)反應(yīng)混合氣體的內(nèi)能和比熱比定義，僅看作物理意義不明確的中間變量。

所有解耦算法在計(jì)算源項(xiàng)時(shí)隱含假設(shè)氣體微團(tuán)總密度和速度保持不變，即

(13)

(14)

代入源項(xiàng)方程E′后得到

(15)

引入等效內(nèi)能E′和等效比熱比γ′以后，將守恒變量U0進(jìn)一步分解為

U1=[ρρuρvρwE′]T

(16)

U2=[ρ1ρ2…ρn-1]T

(17)

采用時(shí)間分裂法把式(1)分解為

(18)

(19)

(20)

U3=[E′ρ1ρ2…ρn-1]T

(21)

(22)

式(18)～式(20)構(gòu)成劉君算法的出發(fā)方程，它和傳統(tǒng)解耦算法的差異主要有兩點(diǎn)：① 源項(xiàng)方程組包括等效內(nèi)能E′;② 流動(dòng)方程的原始守恒變量分解為U0=[U1U2]。

值得說(shuō)明的是，2003年Deiterding[32]也采用等效內(nèi)能E′構(gòu)造了解耦算法，但論文沒(méi)有給出明確的物理解釋?zhuān)瑫r(shí)沒(méi)有下面的二次解耦過(guò)程。

1.2 算法實(shí)現(xiàn)

點(diǎn)隱算法計(jì)算量隨求解變量非線性增加是由源項(xiàng)線化矩陣D=?S/?U∈R(n+4)×(n+4)引起的。全隱算法除D矩陣外，至少還涉及到對(duì)流項(xiàng)的系數(shù)矩陣，以式(1)中x方向?yàn)槔?，?duì)流項(xiàng)通量的系數(shù)矩陣和維數(shù)為A=?F/?U∈R(n+4)×(n+4)。因此計(jì)算量比點(diǎn)隱算法更大。全隱算法隨組元變化不但需要改變D和A的維數(shù)，其表達(dá)形式也需要重新推導(dǎo)。2000年國(guó)外學(xué)者Choi等[53-54]比較特定算例以后得出的結(jié)論表明點(diǎn)隱算法的計(jì)算量與方程數(shù)目呈2次方規(guī)律增加，全隱算法為3次方或更高。

線化矩陣在解耦算法中不是必須的，采用剛性常微分方程計(jì)算程序CHEMEQ、LSODE求解化學(xué)反應(yīng)方程組時(shí)無(wú)需計(jì)算Jacobian矩陣，但在構(gòu)造流動(dòng)方程計(jì)算格式的過(guò)程中大多會(huì)用到系數(shù)矩陣。用于超聲速流動(dòng)的計(jì)算格式為了體現(xiàn)擾動(dòng)沿特征線傳播特性，在構(gòu)造通量分裂格式時(shí)需要系數(shù)矩陣的特征值。傳統(tǒng)解耦算法沒(méi)有二次解耦，組元變化后需要重新推導(dǎo)表達(dá)形式并修改計(jì)算程序。

以上算法的控制方程與反應(yīng)機(jī)理相關(guān)，不利于解耦算法的工程化應(yīng)用。守恒變量U1∈R5對(duì)應(yīng)的方程形式與量熱完全氣體完全相同，無(wú)論顯式格式或隱式格式離散以后得到的能夠體現(xiàn)格式精度和特點(diǎn)的右端項(xiàng)與反應(yīng)機(jī)理無(wú)關(guān)，編程時(shí)僅需將量熱完全氣體的比熱比從常數(shù)修改為等效比熱比變量，直接采用原來(lái)的各種算法。根據(jù)解耦算法的機(jī)理分析，計(jì)算流動(dòng)方程時(shí)流體微團(tuán)內(nèi)部組元變換規(guī)律和總密度一樣，因此，守恒變量U2中各組元的通量分裂格式均采用U1中密度ρ的表達(dá)式，這樣形成的右端項(xiàng)與組元數(shù)目相關(guān)，但是計(jì)算表達(dá)式是統(tǒng)一的，不需重新推導(dǎo)公式和修改程序，易于實(shí)現(xiàn)軟件工程化。

用于隱格式求解時(shí)，守恒變量U1系數(shù)矩陣A1=?F1/?U1∈R5×5與量熱完全氣體完全相同，守恒變量U2系數(shù)矩陣A2=?F2/?U2∈R(n-1)(n-1)為對(duì)角矩陣。如果U1和U2耦合計(jì)算，總變量U0的系數(shù)矩陣由A1和A2按照塊對(duì)角構(gòu)成，形式簡(jiǎn)單，組元數(shù)目增減也不需要重新編程。U1和U2也可以不耦合，在后面新進(jìn)展中介紹。

按照以上解耦算法建立流體和反應(yīng)求解模塊，化學(xué)反應(yīng)方程的求解涉及到剛性問(wèn)題，因此反應(yīng)源項(xiàng)方程表示為L(zhǎng)s,而流動(dòng)算子為L(zhǎng)f，為保證時(shí)間二階精度，總算子采用的形式為

(23)

傳統(tǒng)解耦算法在源項(xiàng)方程時(shí)間推進(jìn)以后僅改變組元，反應(yīng)過(guò)程引起的壓力等流場(chǎng)參數(shù)變化需要等到流動(dòng)方程求解以后才能體現(xiàn)，式(3)中源項(xiàng)的第5項(xiàng)對(duì)應(yīng)為能量項(xiàng)，其值為0，說(shuō)明在進(jìn)行化學(xué)反應(yīng)計(jì)算中壓力等流場(chǎng)參數(shù)未更新，只改變了組元的質(zhì)量分?jǐn)?shù)。相對(duì)于傳統(tǒng)解耦算法，劉君算法同樣保持氣體微團(tuán)總內(nèi)能E不變，式(10)中源項(xiàng)的第5項(xiàng)對(duì)應(yīng)為能量項(xiàng)，其表達(dá)形式涉及到壓力、溫度等物性參數(shù)的求解，通過(guò)對(duì)源項(xiàng)方程的計(jì)算既改變組元的質(zhì)量分?jǐn)?shù)，同時(shí)也計(jì)算了化學(xué)能轉(zhuǎn)化為壓力能的熱力學(xué)過(guò)程。采用圖1可以說(shuō)明兩種算法的差異，圖中使用不同形狀代表流體微團(tuán)的狀態(tài)參數(shù)，Ci表示不同狀態(tài)下組元的質(zhì)量分?jǐn)?shù)矩陣。圖1(a)代表傳統(tǒng)解耦算法，流體微團(tuán)通過(guò)反應(yīng)算子Ls計(jì)算源項(xiàng)方程得到不同組元以后，改變了流體微團(tuán)內(nèi)部的質(zhì)量分?jǐn)?shù)，而壓力、總密度未發(fā)生改變，然后通過(guò)流動(dòng)算子Lf的求解確定流體微團(tuán)熱力學(xué)狀態(tài)參數(shù)；從圖1(b)中可以清晰地看出劉君算法的物理意義，通過(guò)反應(yīng)源項(xiàng)算子Ls求解源項(xiàng)方程，此過(guò)程模擬了流體微團(tuán)在當(dāng)?shù)亟^熱(E不變)、定容(ρ不變)的熱力學(xué)系統(tǒng)內(nèi)發(fā)生的化學(xué)反應(yīng)過(guò)程，狀態(tài)參數(shù)和組元變化耦合求解，通過(guò)Lf求解流動(dòng)方程，模擬組元質(zhì)量分?jǐn)?shù)yi保持不變的流體微團(tuán)沿流線運(yùn)動(dòng)(凍結(jié)流)，包括密度在內(nèi)的狀態(tài)參數(shù)隨位置變化。

圖1 傳統(tǒng)解耦算法和劉君算法的原理差異Fig.1 Difference between classical uncoupled algorithm and modified algorithm

2 發(fā)展歷程

劉君[13]于1993年在其博士學(xué)位論文中提出了一種解耦算法，其研究團(tuán)隊(duì)20多年來(lái)不斷對(duì)算法進(jìn)行改進(jìn)，在介紹最新進(jìn)展前先簡(jiǎn)要回顧發(fā)展歷史。

2.1 有限差分法

最早的劉君算法采用時(shí)空二階的NND2M(Nonoscillatory Nonfree Dissipative second order Modified)單步差分格式求解流動(dòng)偏微分方程組，對(duì)流項(xiàng)通量采用VanLeer格式，反應(yīng)常微分方程組采用梯形公式求解，由于源項(xiàng)計(jì)算耗時(shí)很少，總算子采用式(23)。通過(guò)模擬彈道靶試驗(yàn)[55]來(lái)檢驗(yàn)算法的精度，計(jì)算發(fā)現(xiàn)彈丸飛行速度大于氣體爆轟速度時(shí)“超爆轟”模態(tài)效果很好，駐點(diǎn)線上Von Neumann Spike現(xiàn)象的位置和峰值與國(guó)內(nèi)外文獻(xiàn)結(jié)果完全一致，表明解耦算法的空間精度達(dá)到預(yù)期，但是針對(duì)彈丸飛行速度接近爆轟速度的“跨爆轟”模態(tài)始終沒(méi)有模擬出振蕩現(xiàn)象，當(dāng)時(shí)懷疑算法時(shí)間精度不滿足要求，一直未在學(xué)術(shù)期刊發(fā)表。直到2000年看到文獻(xiàn)[53-54]才發(fā)現(xiàn)網(wǎng)格尺度和反應(yīng)機(jī)理的重要性，最終模擬出燃燒陣面進(jìn)行周期性振蕩的穩(wěn)定燃燒流場(chǎng)[15]，進(jìn)一步分別對(duì)氫氣和空氣燃燒的3種反應(yīng)動(dòng)力學(xué)模型進(jìn)行比較研究[17]，得出的結(jié)論認(rèn)為組元HO2和H2O2及其相關(guān)基元反應(yīng)是模擬振蕩燃燒的關(guān)鍵因素，采用復(fù)雜反應(yīng)機(jī)理能夠準(zhǔn)確模擬混合氣體的點(diǎn)火延遲特性。在此基礎(chǔ)上模擬了高速?gòu)椡柙跉錃?空氣預(yù)混氣體、甲烷/空氣預(yù)混氣體中形成的流場(chǎng)[18]，通過(guò)對(duì)比不同反應(yīng)動(dòng)力學(xué)模型的計(jì)算耗時(shí)后得到的結(jié)論認(rèn)為劉君解耦算法的計(jì)算量與求解變量呈線性規(guī)律變化。為了考察計(jì)算效率，采用40萬(wàn)網(wǎng)格和19種組元、66個(gè)基元反應(yīng)的甲烷/空氣燃燒模型，在主頻為1 GHz、內(nèi)存為128 MB的Winbook6630筆記本電腦上完成了超燃沖壓加速器三維流場(chǎng)模擬。其后，采用該算法針對(duì)固體火箭發(fā)動(dòng)機(jī)裝藥裂紋、超聲速飛行器噴流控制系統(tǒng)和液體火箭發(fā)動(dòng)機(jī)火箭底部阻力等包含有非平衡流動(dòng)的實(shí)際工程問(wèn)題開(kāi)展研究并取得實(shí)效[56-59]。

通過(guò)對(duì)工程應(yīng)用的經(jīng)驗(yàn)總結(jié)，發(fā)現(xiàn)采用梯形公式求解存在剛性問(wèn)題的非線性常微分方程組是造成收斂魯棒性較差等問(wèn)題的主要原因。經(jīng)過(guò)對(duì)算法進(jìn)行改進(jìn)，采用VODE和α-QSS兩種求解器替換梯形公式解決了上述問(wèn)題[60]，但是計(jì)算效率明顯降低，尤其是在總算子采用式(23)模擬“超爆轟”模態(tài)時(shí)，Ls算子的CPU時(shí)間超過(guò)Lf算子。在其他計(jì)算條件相同的情況下，流場(chǎng)穩(wěn)定以后比較每個(gè)時(shí)間步內(nèi)Ls所占CPU時(shí)間比例，采用1/10流動(dòng)時(shí)間步長(zhǎng)的梯形公式為8%、VODE為84%、α-QSS為53%[24]。

2.2 有限體積法

對(duì)式(1)進(jìn)行空間積分，通過(guò)高斯公式把體積分變化為面積分，得到有限體積法的出發(fā)方程為

(24)

化學(xué)反應(yīng)過(guò)程采用Arrhenius模型，源項(xiàng)S表示為組元密度、溫度的指數(shù)型函數(shù)?？臻g一階精度有限體積法假設(shè)物理量在控制體內(nèi)為常數(shù)，源項(xiàng)的體積分可以直接得到,其表達(dá)式為

(25)

(26)

時(shí)間分裂法早期被用來(lái)解決模擬復(fù)雜工程問(wèn)題時(shí)所遇到的計(jì)算機(jī)內(nèi)存不足難題，其理論基礎(chǔ)是在空間點(diǎn)上采用時(shí)間方向迭代降低空間算子維數(shù)，以增加CPU時(shí)間為代價(jià)減少內(nèi)存開(kāi)銷(xiāo)。采用解耦算法在空間點(diǎn)上進(jìn)行時(shí)間分裂法的同時(shí)可以解決有限體積法的源項(xiàng)精度問(wèn)題。時(shí)間分裂以后流動(dòng)偏微分方程不包含化學(xué)反應(yīng)源項(xiàng)，因此不需要計(jì)算源項(xiàng)函數(shù)的體積分；源項(xiàng)函數(shù)只在化學(xué)反應(yīng)的常微分方程中出現(xiàn)，與空間離散無(wú)關(guān)；如果偏微分方程和常微分方程均采用時(shí)間二階精度算法，決定整個(gè)化學(xué)非平衡流動(dòng)計(jì)算空間精度的是流動(dòng)方程的空間離散精度。基于解耦算法的有限體積法理論上具有時(shí)空二階精度。實(shí)際上利用有限體積法空間平均特性還可以建立不需要求解U2方程組的優(yōu)化算法。

3 研究新進(jìn)展

通過(guò)對(duì)國(guó)內(nèi)外研究現(xiàn)狀調(diào)研，認(rèn)為提高計(jì)算效率是解決化學(xué)非平衡流動(dòng)模擬在實(shí)際工程應(yīng)用中發(fā)揮效益的關(guān)鍵，近期對(duì)解耦算法理論和工程應(yīng)用進(jìn)行系統(tǒng)整理后發(fā)現(xiàn)，計(jì)算效率還有改進(jìn)空間，提出如下優(yōu)化算法和處理技術(shù)。

3.1 流動(dòng)方程優(yōu)化算法

結(jié)合劉君解耦算法機(jī)理和有限體積法空間平均特性，提出如圖2所示的優(yōu)化算法：假設(shè)在Δt時(shí)間步內(nèi)，流體微團(tuán)從控制體單元k流出到單元3，從相鄰單元1和2流入。由于U1和U2方程解耦，可以單獨(dú)計(jì)算U1的方程，記錄進(jìn)出邊界的質(zhì)量流量，例如，單元1外法向速度為un,流量為m1=(ρunSΔt)1,其中S為面元面積，流場(chǎng)從tn時(shí)刻推進(jìn)到tn+1時(shí)刻。優(yōu)化解耦算法的流動(dòng)方程本質(zhì)是模擬凍結(jié)流，隨流體微團(tuán)進(jìn)出控制體過(guò)程中沒(méi)有發(fā)生化學(xué)反應(yīng)，組元質(zhì)量分?jǐn)?shù)yi保持不變，根據(jù)質(zhì)量守恒，Δt時(shí)間內(nèi)控制體組元i的密度平均值為

(27)

優(yōu)化算法的優(yōu)勢(shì)在于方程求解過(guò)程中不需要計(jì)算U2的偏微分方程組就可以得到U2。

圖2 優(yōu)化算法機(jī)理示意圖Fig.2 Schematic of optimization algorithm mechanism

3.2 驗(yàn)證算例

采用250×200的計(jì)算網(wǎng)格、修正的Jachimowski(9種組元、19個(gè)基元反應(yīng))化學(xué)反應(yīng)機(jī)理模擬經(jīng)典的激波誘導(dǎo)燃燒試驗(yàn)，試驗(yàn)來(lái)流條件壓力p∞為42 650 Pa，溫度T∞為293 K。模擬彈丸飛行馬赫數(shù)Ma=6.46的超爆轟模態(tài)得到的密度云圖中激波陣面和燃燒陣面位置與試驗(yàn)比較如圖3所示，試驗(yàn)采用的彈丸半徑R為7.5 mm，圖中顯示的橫縱坐標(biāo)均為無(wú)量綱坐標(biāo)?；瘜W(xué)反應(yīng)的誘導(dǎo)時(shí)間主要取決于流場(chǎng)中的溫度，即溫度越高化學(xué)反應(yīng)誘導(dǎo)時(shí)間越短?；旌蠚怏w經(jīng)過(guò)彈丸頭部附近正激波后，由于激波強(qiáng)度大、溫度高，導(dǎo)致此區(qū)域附近誘導(dǎo)時(shí)間最短。而誘導(dǎo)區(qū)厚度是由誘導(dǎo)時(shí)間決定的，在駐點(diǎn)流線上激波陣面與燃燒陣面幾乎重合在一起，垂直于弓形激波的速度分量等于爆轟速度，文獻(xiàn)[63]也對(duì)此現(xiàn)象有詳細(xì)的解釋。通過(guò)對(duì)比發(fā)現(xiàn)數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果是一致的。遠(yuǎn)離駐點(diǎn)流線，激波強(qiáng)度和流場(chǎng)溫度開(kāi)始降低，所以誘導(dǎo)時(shí)間增加，導(dǎo)致激波陣面與燃燒陣面逐漸分離。在出口附近，模擬結(jié)果與試驗(yàn)結(jié)果之間存在一定的誤差，主要原因是Jachimowski反應(yīng)機(jī)理中放熱效果更加顯著，斜激波強(qiáng)度高于試驗(yàn)值，使得模擬結(jié)果中陣面的位置均高于試驗(yàn)結(jié)果。

圖3 計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比Fig.3 Comparison between calculated and test results

沿駐點(diǎn)線壓力、溫度分布和文獻(xiàn)[24]的比較結(jié)果如圖4所示。ZND(Zeldovich Neumann Doering)爆轟理論[64]認(rèn)為：前導(dǎo)激波將可燃?xì)怏w壓縮加熱到其自燃狀態(tài)，經(jīng)過(guò)誘導(dǎo)區(qū)后開(kāi)始發(fā)生化學(xué)反應(yīng)并迅速釋放化學(xué)能，誘導(dǎo)反應(yīng)區(qū)后溫度、壓強(qiáng)、密度等狀態(tài)量急劇變化，借助化學(xué)反應(yīng)釋放的能量使燃燒氣體膨脹并支撐前導(dǎo)激波以一定馬赫數(shù)傳播。類(lèi)似于傳統(tǒng)的燃燒現(xiàn)象，化學(xué)反應(yīng)區(qū)是一個(gè)溫度升高、壓力基本持平的過(guò)程，燃燒氣體的熱力學(xué)參數(shù)在放熱區(qū)的末端達(dá)到CJ(Chapman Jouget)爆轟狀態(tài)。ZND模型考慮了爆轟燃燒的非平衡過(guò)程，并提出了能夠解釋爆轟波自持傳播的化學(xué)反應(yīng)過(guò)程與前導(dǎo)激波相互作用的機(jī)制。本文采用的基元化學(xué)反應(yīng)模型中，化學(xué)反應(yīng)速率遵循Arrhenius規(guī)律，符合ZND模型的基本假設(shè)，從模擬結(jié)果中顯示出Von Neumann Spike現(xiàn)象，壓力峰值是由于波后的劇烈燃燒和能量釋放，難以維持等壓燃燒造成的。顯然，如果能量釋放率越高，那么Von Neumann Spike現(xiàn)象將更加明顯。從計(jì)算結(jié)果可以看出本文使用的方法很好地模擬出Von Neumann Spike現(xiàn)象，與文獻(xiàn)[24]的結(jié)果也基本一致，說(shuō)明優(yōu)化算法具有很好的精度。

圖5所示為優(yōu)化算法與劉君算法得到的密度云圖對(duì)比。 圖5分為上下兩部分，分別采用兩種方法的收斂解進(jìn)行結(jié)果對(duì)比分析。以上模擬結(jié)果驗(yàn)證了優(yōu)化算法對(duì)于定常流場(chǎng)具有較好的空間分辨率。

圖4 沿駐點(diǎn)線壓力、溫度分布的計(jì)算結(jié)果與文獻(xiàn)[24]試驗(yàn)結(jié)果的比較Fig.4 Comparison of pressure and temperature distribution obtained by calculation with test results of Ref.[24] along stagnation line

圖5 超爆轟密度云圖Fig.5 Density contours of super detonation

圖6 駐點(diǎn)線上密度等值線的時(shí)間歷程(Ma=4.79)Fig.6 History of density contours on stagnation line (Ma=4.79)

為了驗(yàn)證優(yōu)化算法的準(zhǔn)確性，對(duì)Lehr試驗(yàn)[55]中激波誘導(dǎo)周期性振蕩燃燒算例進(jìn)行數(shù)值模擬。與試驗(yàn)來(lái)流條件相同，采用H2/Air混合氣體在壓力p∞為42 650 Pa、溫度T∞為293 K的條件下進(jìn)行計(jì)算，彈丸的飛行馬赫數(shù)分別為Ma=4.79和Ma=4.48。為得到較好的計(jì)算結(jié)果，使用350×300的計(jì)算網(wǎng)格，采用文獻(xiàn)[62]提出的通量格式和限制器函數(shù)分別對(duì)兩種條件下跨爆轟算例進(jìn)行模擬，圖6和圖7分別為采用優(yōu)化算法模擬跨爆轟模態(tài)得到的駐點(diǎn)線上密度的時(shí)間歷程等值線。模擬結(jié)果顯示出明顯的周期性振蕩燃燒現(xiàn)象并且在表1中列出了與其他結(jié)果[53,55]的數(shù)據(jù)對(duì)比。本文的模擬結(jié)果與其他數(shù)據(jù)誤差很小，驗(yàn)證了本文模擬結(jié)果的準(zhǔn)確性。

圖7 駐點(diǎn)線上密度等值線的時(shí)間歷程(Ma=4.48)Fig.7 History of density contours on stagnation line (Ma=4.48)

表1 振蕩頻率對(duì)比Table 1 Comparison of oscillation frequency

數(shù)據(jù)來(lái)源頻率/kHzMa=4.79Ma=4.48Lehr試驗(yàn)[55]712425文獻(xiàn)[53]707426本文結(jié)果706422

3.3 耦合過(guò)程的空間優(yōu)化

優(yōu)化算法求解過(guò)程中不計(jì)算U2,計(jì)算效率僅提高14.6%，主要原因是采用VODE和α-QSS求解器在改善常微分方程組求解穩(wěn)定性的同時(shí)增加了計(jì)算源項(xiàng)的CPU時(shí)間。注意到在很多工程問(wèn)題中化學(xué)反應(yīng)僅發(fā)生在流場(chǎng)局部區(qū)域，在未反應(yīng)或反應(yīng)達(dá)到平衡的空間區(qū)域組元質(zhì)量生成率σi≈0,沒(méi)有必要計(jì)算源項(xiàng)。在t0時(shí)刻整個(gè)預(yù)混燃?xì)獾牧鲌?chǎng)溫度保持相同，采用流動(dòng)算子Lf對(duì)流場(chǎng)參數(shù)進(jìn)行更新以后，激波后區(qū)域的溫度會(huì)逐漸上升，溫度越高化學(xué)反應(yīng)進(jìn)行得越劇烈，計(jì)算得到的組元生成率越大。由于激波前區(qū)域仍保持為初始來(lái)流參數(shù)，不需要進(jìn)行化學(xué)反應(yīng)計(jì)算，而激波后區(qū)域溫度上升，化學(xué)反應(yīng)更加劇烈，因此可以把初始時(shí)刻進(jìn)行的化學(xué)反應(yīng)作為基準(zhǔn)來(lái)對(duì)全場(chǎng)區(qū)域進(jìn)行標(biāo)識(shí)?；谝陨险J(rèn)識(shí)，利用解耦算法的特點(diǎn)，在t0時(shí)刻采用反應(yīng)算子Ls對(duì)所有單元進(jìn)行一次化學(xué)反應(yīng)計(jì)算，同時(shí)記錄下所有單元中組元生成率絕對(duì)值的最大值為

(28)

以式(28)為判據(jù)準(zhǔn)則在隨后的計(jì)算中對(duì)所有控制體單元的質(zhì)量生成率進(jìn)行判斷，當(dāng)滿足判斷條件

(29)

將控制體及其相鄰控制體進(jìn)行標(biāo)識(shí)，整理流場(chǎng)完成分類(lèi)后僅在標(biāo)識(shí)過(guò)的控制體內(nèi)調(diào)用Ls算子。

針對(duì)馬赫數(shù)Ma=6.46的超爆轟模態(tài)，通過(guò)t0時(shí)刻計(jì)算得到的準(zhǔn)則參數(shù)ε1=2.66×10-13kg/s進(jìn)行辨識(shí)，CPU時(shí)間隨計(jì)算步數(shù)變化曲線如圖8所示，通過(guò)標(biāo)識(shí)得到的化學(xué)反應(yīng)計(jì)算區(qū)域如圖9中陰影部分所示，大約占總體計(jì)算區(qū)域的2/3左右，圖9(a)中對(duì)局部計(jì)算網(wǎng)格區(qū)域進(jìn)行放大處理，圖中加粗黑線為進(jìn)行質(zhì)量生成率判別后進(jìn)行標(biāo)記的反應(yīng)區(qū)分界面。表2列出不同計(jì)算步數(shù)時(shí)CPU時(shí)間對(duì)比，通過(guò)采用質(zhì)量生成率判據(jù)，化學(xué)反應(yīng)算子計(jì)算效率明顯提高，總CPU時(shí)間也隨之降低?？臻g優(yōu)化算法同樣適用于跨爆轟模態(tài)，圖9(b)流場(chǎng)中兩條曲線分別代表兩個(gè)時(shí)刻的化學(xué)反應(yīng)區(qū)域，由于化學(xué)反應(yīng)區(qū)域隨時(shí)間發(fā)生變化，沒(méi)有比較計(jì)算效率。

圖8 CPU時(shí)間隨計(jì)算步數(shù)的變化曲線Fig.8 Variation curves of CPU time with calculated steps

圖9 添加判據(jù)后進(jìn)行化學(xué)反應(yīng)區(qū)域示意圖Fig.9 Schematic of chemical reaction region after addition of criterion

表2 總CPU時(shí)間對(duì)比Table 2 Comparison of total CPU time

計(jì)算步數(shù)總CPU時(shí)間/s無(wú)判據(jù)準(zhǔn)則ε有判據(jù)準(zhǔn)則ε1000023222.515121.22000047096.932968.14000095130.471515.6

4 結(jié)束語(yǔ)

在回顧1993年劉君提出的化學(xué)非平衡流動(dòng)解耦算法及其發(fā)展歷史之后，采用有限體積法空間平均特性對(duì)劉君解耦算法進(jìn)行了改進(jìn)并提出耦合過(guò)程中的空間優(yōu)化算子，模擬結(jié)果顯示優(yōu)化算法可以顯著提高化學(xué)非平衡流模擬的計(jì)算效率。在此基礎(chǔ)上，今后研究的重點(diǎn)將轉(zhuǎn)向集成國(guó)外新成果來(lái)提高化學(xué)反應(yīng)算子的計(jì)算效率。

[1] 常亞超. 基于解耦法的柴油和生物柴油表征燃料骨架反應(yīng)機(jī)理研究[D]. 大連： 大連理工大學(xué), 2016: 1-13.

CHANG Y C. Investigation of skeletal chemical mechanisms for diesel and biodiesel surrogate fuel based on decoupling methodology[D]. Dalian: Dalian University of Technology, 2016: 1-13 (in Chinese).

[2] 梁劍寒, 王承堯. 超燃沖壓發(fā)動(dòng)機(jī)燃燒室三維化學(xué)反應(yīng)流場(chǎng)的數(shù)值模擬[J]. 推進(jìn)技術(shù), 1997, 18(4): 1-4.

LIANG J H, WANG C Y. 3D numerical simulation of supersonic reacting flows in scramjet combustion[J]. Journal of Propulsion Technology, 1997, 18(4): 1-4 (in Chinese).

[3] ROGERS R C, CHINITZ W. Using a global hydrogen-air combustion model in turbulent reacting flow calculations[J]. AIAA Journal, 1983, 21(4): 586-592.

[4] 黃生洪, 徐勝利, 劉小勇. 煤油超燃沖壓發(fā)動(dòng)機(jī)兩相流場(chǎng)數(shù)值模擬(I)數(shù)值校驗(yàn)及總體流場(chǎng)特征[J].推進(jìn)技術(shù), 2004, 25(6): 484-490.

HUANG S H, XU S L, LIU X Y. Numerical investigation on two phase flow of akerosene-fueled scramjet with 3D cavity (I) Numerical calibration and characteristics of general flow[J]. Journal of Propulsion Technology, 2004, 25(6): 484-490 (in Chinese).

[5] 楊順華， 樂(lè)嘉陵， 趙慧勇， 等. 煤油超燃沖壓發(fā)動(dòng)機(jī)三維大規(guī)模并行數(shù)值模擬[J]. 計(jì)算物理， 2009， 26(4)： 534-540.

YANG S H, LE J L, ZHAO H Y, et al. Three-dimensional massively parallel numerical simulation of kerosene-fueled scramjet[J]. Chinese Journal of Computational Physics, 2009, 26(4): 534-540 (in Chinese).

[6] YANENKO N N. The method of fractional steps[M]. Heidelberg: Springer-verlag, 1971: 42-50.

[7] DESCOMBES S, DUMONT T. Numerical simulation of a stroke: Computational problems and methodology[J]. Progress in Biophysics & Molecular Biology, 2008, 97(1): 40-53.

[8] DUMONT T, DUARTE M, DESCOMBES S, et al. Simulation of human ischemic stroke in realistic 3D geometry[J]. Communications in Nonlinear Science & Numerical Simulation, 2013, 18(6): 1539-1557.

[9] CUI S, KURGANOV A, MEDOVIKOV A. Particle methods for PDEs arising in financial modeling[J]. Applied Numerical Mathematics, 2015, 93: 123-139.

[10] DESCOMBES S, THALHAMMER M. The Lie-Trotter splitting for nonlinear evolutionary problems with critical parameters: A compact local error representation and application to nonlinear schr?dinger equations in the semiclassical regime[J]. IMA Journal of Numerical Analysis, 2013, 33(2): 722-745.

[11] KULIKOV I, VOROBYOV E. Using the PPML approach for constructing a low-dissipation, operator-splitting scheme for numerical simulations of hydrodynamic flows[J]. Journal of Computational Physics, 2016, 317: 318-346.

[12] CAMBIER J, ADELMAN H, MENEES G P. Numerical simulations of an oblique detonation wave engine[J]. Journal of Propulsion and Power, 1990, 6(3): 315-323.

[13] 劉君. 超音速完全氣體和H2/O2燃燒非平衡氣體的復(fù)雜噴流流場(chǎng)數(shù)值模擬[D]. 綿陽(yáng)： 中國(guó)空氣動(dòng)力研究與發(fā)展中心, 1993: 56-65.

LIU J. Numerical simulation of the complex jet flow field of supersonic gas and H2/O2combustion[D]. Mianyang: China Aerodynamics Research and Development Center, 1993: 56-65 (in Chinese).

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

LIU J, ZHANG H X, GAO S C. A new uncoupled method for numerical simulation of non-equilibrium flow[J]. Journal of National University of Defense Technology, 2000, 22(5): 19-22 (in Chinese).

[15] 劉君. 沖壓加速器非平衡流動(dòng)數(shù)值模擬[J]. 彈道學(xué)報(bào), 2002, 14(4): 31-35.

LIU J. Numerical study on the non-equilibrium flow of ram accelerator in the combustive mixture gas[J]. Journal of Ballistics, 2002, 14(4): 31-35 (in Chinese).

[16] 劉君. 非平衡流計(jì)算方法及其模擬激波誘導(dǎo)振蕩燃燒[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2003, 21(1): 53-58.

LIU J. A new non-equilibrium numerical method and simulation of oscillating shock-induced combustion[J]. Acta Aerodynamica Sinica, 2003, 21(1): 53-58 (in Chinese).

[17] 劉君. 化學(xué)動(dòng)力學(xué)模型對(duì)H2/Air超燃模擬的影響[J]. 推進(jìn)技術(shù), 2003, 24(1): 67-70.

LIU J. Numerical study on chemical mechanism in supersonic H2/Air mixture gas flow[J]. Journal of Propulsion Technology, 2003, 24(1): 67-70 (in Chinese).

[18] 劉君. 甲烷/空氣超聲速燃燒流動(dòng)數(shù)值模擬[J]. 推進(jìn)技術(shù), 2003, 24(4): 296-299.

LIU J. Numerical study for CH4/air supersonic combustion[J]. Journal of Propulsion Technology, 2003, 24(4): 296-299 (in Chinese).

[19] 劉君, 劉瑜, 周松柏. 基于新型解耦算法的激波誘導(dǎo)燃燒過(guò)程數(shù)值模擬[J]. 力學(xué)學(xué)報(bào), 2010, 42(3): 572-578.

LIU J, LIU Y, ZHOU S B. Simulation of shock induced combustion based on a novel uncoupled method[J]. Chinese Journal of Theoretical and Applied Mechanics, 2010, 42(3): 572-578 (in Chinese).

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

LIU Y. Investigations into numerical methods of chemical nonequilibrium flow and its application to simulation of shock-induced combustion phenomena[D]. Changsha： National University of Defense Technology, 2008: 20-24 (in Chinese).

[21] 劉瑜, 劉君, 白曉征. 基于新型非結(jié)構(gòu)有限體積解耦算法的激波誘導(dǎo)燃燒數(shù)值模擬[J]. 國(guó)防科技大學(xué)學(xué)報(bào), 2011, 33(6): 139-144.

LIU Y, LIU J, BAI X Z. Numerical simulation of shock-induced combustion with a new uncoupled algorithm in university finite volume method[J]. Journal of National University of Defense Technology, 2011, 33(6): 139-144 (in Chinese).

[22] 劉瑜, 劉君, 唐玲艷, 等. 一種求解化學(xué)非平衡流動(dòng)的新型解耦算法[J]. 力學(xué)學(xué)報(bào), 2015, 47(1): 82-94.

LIU Y, LIU J, TANG L Y, et al. A novel uncoupled algorithm for solving chemical nonequilibrium flows[J]. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(1): 82-94 (in Chinese).

[23] 汪洪波, 孫明波, 梁劍寒, 等. 含雙時(shí)間步法的化學(xué)非平衡流解耦算法[J]. 計(jì)算物理, 2010, 27(5): 685-691.

WANG H B, SUN M B, LIANG J H, et al. An uncoupled method with dual time-step for chemical nonequilibrium flows[J]. Chinese Journal of Computational Physics, 2010, 27(5): 685-691 (in Chinese).

[24] 周松柏. 超聲速內(nèi)外流干擾的數(shù)值方法研究及其實(shí)驗(yàn)驗(yàn)證與應(yīng)用[D]. 長(zhǎng)沙： 國(guó)防科學(xué)技術(shù)大學(xué), 2009: 76-92.

ZHOU S B. Research on numerical method, experimental verification and application in simulation jet interaction in supersonic external flow[D]. Changsha： National University of Defense Technology, 2009: 76-92 (in Chinese).

[25] 孫明波, 梁劍寒, 王振國(guó). 非平衡流解耦方法及其計(jì)算激波誘導(dǎo)燃燒的應(yīng)用驗(yàn)證[J]. 航空動(dòng)力學(xué)報(bào), 2008, 23(11): 2055-2061.

SUN M B, LIANG J H, WANG Z G. Validation of an uncoupled solver of non-equilibrium flow for shock-induced combustion[J]. Journal of Aerospace Power, 2008, 23(11): 2055-2061 (in Chinese).

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

LIU S J, SUN M B, LIN Z Y, et al. Numerical research on blunt body shock-induced oscillating combustion phenomena[J]. Chinese Journal of Theoretical and Applied Mechanics, 2010, 42(4): 597-606 (in Chinese).

[27] 韓旭, 周進(jìn), 林志勇, 等. 突躍型與平滑型斜爆震波起爆機(jī)制數(shù)值模擬[J]. 航空動(dòng)力學(xué)報(bào), 2012, 27(12): 2674-2680.

HAN X, ZHOU J, LIN Z Y, et al. Initiation mechanism investigation of saltation and smoothness oblique detonation waves[J]. Journal of Aerospace Power, 2012, 27(12): 2674-2680 (in Chinese).

[28] 劉君, 楊彥廣. 帶有橫噴控制的導(dǎo)彈流場(chǎng)數(shù)值模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2004, 22(3): 309-312.

LIU J, YANG Y G. Numerical simulation of lateral jet control of a hypersonic missile[J]. Acta Aerodynamica Sinica, 2004, 22(3): 309-312 (in Chinese).

[29] 楊彥廣, 劉君, 唐志共. 橫向噴流干擾中的真實(shí)氣體效應(yīng)研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2006, 24(1): 28-33.

YANG Y G, LIU J, TANG Z G. A study of real gas effects on lateral jet interaction[J]. Acta Aerodynamica Sinica, 2006, 24(1): 28-33 (in Chinese).

[30] 劉世杰, 林志勇, 孫明波, 等. 旋轉(zhuǎn)爆震波發(fā)動(dòng)機(jī)二維數(shù)值模擬[J]. 推進(jìn)技術(shù), 2010, 31(5): 634-640.

LIU S J, LIN Z Y, SUN M B, et al. Two-dimensional numerical simulation of rotating detonation wave engine[J]. Journal of Propulsion Technology, 2010, 31(5): 634-640 (in Chinese).

[31] 余文杰, 余永剛. 底排裝置尾部化學(xué)非平衡流的數(shù)值模擬[J]. 推進(jìn)技術(shù), 2015, 36(11): 1610-1615.

YU W J, YU Y G. Numerical simulation of base flow with chemical non-equilibrium for base bleed equipment[J]. Journal of Propulsion Technology, 2015, 36(11): 1610-1615 (in Chinese).

[32] DEITERDING R. Parallel adaptive simulation of multidimensional detonation structures[D]. Cottbus: Brandenburg University of Technology, 2003: 18-27.

[33] ZIEGLER J L, DEITERDING R, SHEPHERD J E, et al. An adaptive high-order hybrid scheme for compressive, viscous flows with detailed chemistry[J]. Journal of Computational Physics, 2011, 230(20): 7598-7630.

[34] BILLET G, RYAN J. A runge kutta discontinuous Galerkin approach to solve reactive flows: The hyperbolic operator[J]. Journal of Computational Physics, 2011, 230(4): 1064-1083.

[35] BILLET G, RYAN J, BORREL M. A Runge Kutta discontinuous Galerkin approach to solve reactive flows on conforming hybrid grids: The parabolic and source operators[J]. Computers and Fluids, 2014, 95: 98-115.

[36] LV Y, IHME M. Discontinuous Galerkin method for multicomponent chemically reacting flows and combustion[J]. Journal of Computational Physics, 2014, 270: 105-137.

[37] GLOWINSKI R, OSHER S J, YIN W. Splitting methods in communication, imaging, science, and engineering[M]. Heidelberg: Springer-Verlag, 2016: 627-641.

[38] TROTTER H F. On the product of semi-groups of operators[J]. Proceedings of the American Mathematical Society, 1959, 10(4): 545-551.

[39] STRANG G. On the construction and comparison of difference schemes[J]. SIAM Journal on Numerical Analysis, 1968, 5(3): 506-517.

[40] LU P. A consistent-splitting approach to computing stiff steady-state reacting flows with adaptive chemistry[J]. Combustion Theory and Modelling, 2003, 7(2): 383-399.

[41] REN Z, XU C, LU T, et al. Dynamic adaptive chemistry with operator splitting schemes for reactive flow simulations[J]. Journal of Computational Physics, 2014, 263: 19-36.

[42] DUARTE M. Adaptive numerical methods in time and space for the simulation of multi-scale reaction fronts[D]. Paris: Ecole Centrale Paris, 2011: 104-128.

[43] DUARTE M, DESCOMBES S, TENAUD C, et al. Time-space adaptive numerical methods for the simulation of combustion fronts[J]. Combustion and Flame, 2013, 160(6): 1083-1101.

[44] HELLANDER A, LAWSON M J, DRAWERT B, et al. Local error estimates for adaptive simulation of the reaction-diffusion master equation via operator splitting[J]. Journal of Computational Physics, 2014, 266: 89-100.

[45] SPORTISSE B. An analysis of operator splitting techniques in the stiff case[J]. Journal of Computational Physics, 2000, 161(1): 140-168.

[46] SANTILLANA M, ZHANG L, YANTOSCA R. Estimating numerical errors due to operator splitting in global atmospheric chemistry models: Transport and chemistry[J]. Journal of Computational Physics, 2015, 305: 909-917.

[47] GOU X, SUN W, CHEN Z, et al. A dynamic multi-timescale method for combustion modeling with detailed and reduced chemical kinetic mechanisms[J]. Combustion and Flame, 2010, 157(6): 1111-1121.

[48] SUN W, GOU X, ELASRAG H A, et al. Multi-timescale and correlated dynamic adaptive chemistry modeling of ignition and flame propagation using a real jet fuel surrogate model[J]. Combustion and Flame, 2015, 162(4): 1530-1539.

[49] YANG G, LIU Y, REN Z, et al. A dynamic adaptive method for hybrid integration of stiff chemistry[J]. Combustion and Flame, 2015, 162(2): 287-295.

[50] XU C, GAO Y, REN Z, et al. A sparse stiff chemistry solver based on dynamic adaptive integration for efficient combustion simulations[J]. Combustion and Flame, 2016, 172: 183-193.

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

LIU J, ZHOU S B, XU C G. The numerical simulation method of supersonic flow and its application[M]. Changsha: National University of Defense Technology Press, 2008: 61-63 (in Chinese).

[52] GUPTA R N, YOS J M, THOMPSON R A. A review of reaction rates and thermodynamic and transport properties for the 11-species air model for chemical and thermal non-equilibrium calculations to 30000 K[J]. Philosophy, 1990, 92(6): 32-34.

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

[54] CHOI J Y, JEUNG I S, YOON Y. Computational fluid dynamics algorithms for unsteady shock-induced combustion, Part 2: Comparison[J]. AIAA Journal, 2000, 38(7): 1188-1195.

[55] LEHR H F. Experiments on shock-induced combustion[J]. Acta Astronautica, 1972(1): 589-597.

[56] 楊彥廣, 劉君. 高超聲速主流中橫向噴流干擾非定常特性研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2004, 22(3): 295-302.

YANG Y G, LIU J. Unsteady characteristic research of lateral jet in hypersonic external flow[J]. Acta Aerodynamica Sinica, 2004, 22(3): 295-302 (in Chinese).

[57] 徐春光, 劉君. 某型導(dǎo)彈尾噴流形狀的數(shù)值模擬[J]. 推進(jìn)技術(shù), 2003, 24(2): 141-143.

XU C G, LIU J. Numerical simulation on the plume flow structures for a missile[J]. Journal of Propulsion Technology, 2003, 24(2): 141-143 (in Chinese).

[58] 周松柏, 劉君, 郭正, 等. 激波誘導(dǎo)異質(zhì)氣體界面失穩(wěn)的數(shù)值模擬[J]. 國(guó)防科技大學(xué)學(xué)報(bào), 2009, 31(1): 1-4.

ZHOU S B, LIU J, GUO Z, et al. Numerical simulation for shock induced interfacial instabilities of heterogeneous gases[J]. Journal of National University of Defense Technology, 2009, 31(1): 1-4 (in Chinese).

[59] 崔小強(qiáng), 白曉征, 周偉, 等. 點(diǎn)火沖擊波與藥柱裂紋流固耦合作用機(jī)理數(shù)值模擬研究[J]. 固體火箭技術(shù), 2010, 33(1): 9-12.

CUI X Q, BAI X Z, ZHOU W, et al. Numerical simulation of the fluid-structure interaction mechanism between ignition shock wave and elastic cracks of grain[J]. Journal of Solid Rocket Technology, 2010, 33(1): 9-12 (in Chinese).

[60] 劉瑜. 量熱完全氣體/化學(xué)非平衡氣體流動(dòng)ALE有限體積計(jì)算方法研究[D]. 長(zhǎng)沙： 國(guó)防科學(xué)技術(shù)大學(xué), 2014: 128-134.

LIU Y. Numerical research of ALE finite volume method for calorically perfect gas chemical nonequilibrium gas flow[D]. Changsha: National University of Defense Technology, 2014: 128-134 (in Chinese).

[61] VANLEER B. Towards the ultimate conservative difference scheme: V. A second-order sequel to Godunov’s method[J]. Journal of Computational Physics, 1979, 32(1): 101-136.

[62] ZHANG F, LIU J, CHEN B, et al. Evaluation of rotated upwind schemes for contact discontinuity and strong shock[J]. Computers and Fluids, 2016, 134(1): 11-22.

[63] WILSON G J, MACCORMACK R W. Modeling supersonic combustion using a fully implicit numerical method[J]. AIAA Journal, 1992, 30(4): 1008-1015.

[64] ZELDOVICH Y B. On the theory of the propagation of detonation in gaseous systems[J]. Journal of Experimental and Theoretical Physics, 1940, 10(5): 542-568.