發(fā)動機爆震室氣相化學反應(yīng)機理及爆轟波性質(zhì)的數(shù)值研究

饒飛雄， 丁 玨， 翁培奮， 李孝偉

(上海大學上海市應(yīng)用數(shù)學和力學研究所，上海 200072)

脈沖爆震發(fā)動機(pulse detonation engine, PDE)具有設(shè)計簡單、熱效率高、應(yīng)用范圍廣等優(yōu)點，未來有替代傳統(tǒng)發(fā)動機而成為新型動力裝置的趨勢，因此越來越受到國際研究者的高度重視[1]。

爆燃波通過火焰區(qū)的質(zhì)量擴散和熱擴散使前方反應(yīng)物著火，從而實現(xiàn)火焰?zhèn)鞑2]。而爆轟波則是一道超聲速壓縮激波，它通過前導激波掃過混合物時的絕熱壓縮加熱點燃混合物。因反應(yīng)過程接近定容燃燒，爆轟波產(chǎn)生的能量與效率遠高于普通燃燒[3]。通常，研究者采用一步、兩步[4-6]反應(yīng)模型模擬氣相爆轟波結(jié)構(gòu)。張薇等[7]采用單步化學反應(yīng)模型研究氣相爆轟波胞格尺寸和點火延遲時間關(guān)系；張德良等[8]采用簡化二階段化學反應(yīng)模型模擬爆燃轉(zhuǎn)爆轟和爆轟波在斜面上的馬赫反射；劉云峰等[9]采用二步化學反應(yīng)模型分析爆轟波的三波點和發(fā)展過程。由于一步、兩步反應(yīng)模型簡化了真實爆轟過程鏈式反應(yīng)的復雜性和影響，難以給出模型中相關(guān)參數(shù)，且反應(yīng)產(chǎn)物組成不完全導致爆轟參數(shù)產(chǎn)生較大的偏差，進而影響波陣面結(jié)構(gòu)及其相關(guān)參數(shù)。因此，越來越多的人采用較能正確反映爆轟波化學反應(yīng)本質(zhì)的基元反應(yīng)。Oran等[10]采用基元反應(yīng)模型對二維爆轟波進行模擬來研究能量的釋放規(guī)律；Fukuda 等[11]采用詳細反應(yīng)機理模型研究了壁面條件對爆燃轉(zhuǎn)爆轟的影響；董剛等[12]基于氫氣與空氣詳細化學反應(yīng)模擬爆轟波，分析了爆轟反應(yīng)后的組分濃度和溫度分布；王昌建等[13]采用基元反應(yīng)數(shù)值模擬研究了初始擾動對爆轟波形成的影響；張寶亮等[14]基于有限體積法分析了障礙物與爆轟波的相互作用規(guī)律。

本工作基于敏感性分析，將二十步反應(yīng)簡化為九步反應(yīng)模型，并將詳細基元反應(yīng)模型應(yīng)用于爆轟波數(shù)值模擬，分析爆轟波結(jié)構(gòu)、傳播性質(zhì)及胞格。同時，本工作研究了Ar稀釋和初始壓強的變化對爆轟波參數(shù)的影響, 得出其影響爆轟波特性的相關(guān)規(guī)律，為深入研究爆震發(fā)動機非均相爆轟特性奠定了理論基礎(chǔ)。

1 爆轟波形成和發(fā)展的控制方程和數(shù)值計算方法

1.1 控制方程

爆轟波包括超聲速激波以及緊隨的化學反應(yīng)區(qū)組成的高速燃燒波。對于爆震室可燃氣體爆轟及爆轟波傳播過程，采用二維Euler 控制方程：

式中：u是守恒變量，F(xiàn)i是對流通量，Q為源項，分別為

其中δ是克羅內(nèi)克函數(shù)，ρ是密度，u1、u2是速度的分量，p是壓力，為能量方程的源項，是第K個組分的質(zhì)量生成率，N是化學組分數(shù)目。

YK是第K個組分的質(zhì)量分數(shù)，滿足

能量E定義為：

式中：hs為顯焓。

化學反應(yīng)速率用Arrhenius 公式來計算

式中：Ea為化學反應(yīng)活化能，A為指前因子，R為通用氣體常數(shù)，T為流場溫度。

1.2 計算區(qū)域的網(wǎng)格及初始邊界條件

爆震室直徑為2 cm，充滿初始壓力為一個標準大氣壓、初始溫度為300 K 的氫氣和氧氣混合氣體，其中氣體組分的物質(zhì)的量之比為m(H2):m(O2)=2:1。爆震室左端點火，示意圖如圖1 所示。

圖1 爆轟室示意圖Fig.1 Sketch map of detonation chamber

爆轟的反應(yīng)速率是普通燃燒的103～108倍[3], 為了準確地模擬爆轟波結(jié)構(gòu)，網(wǎng)格尺寸需要滿足精度要求。因此，基于CJ 理論[15]得出的爆轟相關(guān)參數(shù)，本工作采用4 種不同的網(wǎng)格尺寸數(shù)值模擬爆轟波形成和發(fā)展過程，并與實驗和爆轟的CJ 理論值進行比較，結(jié)果見表1。

表1 網(wǎng)格尺寸對爆轟參數(shù)的影響Table 1 Influence of mesh size on detonation parameters

表1 中數(shù)據(jù)顯示：隨著計算網(wǎng)格尺寸的減小，爆轟波傳播速度較接近實驗值，但計算量成倍增大；爆轟壓強變化不大，因有效小數(shù)的限制，壓強在1.911 MPa 附近變化。因此，在滿足計算精度的基礎(chǔ)上考慮計算量，本工作采用0.10 mm 的計算網(wǎng)格。

1.3 數(shù)值計算方法

1.3.1 爆轟的超聲速流場

在超音速流動中，基于密度求解器能在求解Riemann 問題時更能彰顯其優(yōu)越性。針對網(wǎng)格i+交界面，本工作計算對流通量采用Roe 格式[16]：

c為音速，U=nxu+nyv為界面法向速度，nx和ny為界面法向向量在坐標軸上的分量。另外，為了保證計算結(jié)果收斂與穩(wěn)定性，時間步采用10-8s。

1.3.2 化學反應(yīng)動力學

對于激波與緊隨化學反應(yīng)區(qū)的爆轟波，需要求解反應(yīng)區(qū)中多組分的化學反應(yīng)。本工作考慮m(H2):m(O2)=1:0 的8 組元二十步基元反應(yīng)H2/O2模型[17]，具體參數(shù)如表2 所示。

表2 化學反應(yīng)和反應(yīng)參數(shù)Table 2 Chemical reactions and reaction parameters

2 九步基元反應(yīng)模型

帶有反應(yīng)流的數(shù)值模擬中采用詳細基元反應(yīng)機理，計算量增大，同時基于詳細反應(yīng)機理的反應(yīng)流的計算常常出現(xiàn)剛性問題，給數(shù)值模擬帶來了較大的困難。因此，本工作根據(jù)敏感性分析方法[18]，分析反應(yīng)機理中參數(shù)變化對計算結(jié)果的影響程度，消除詳細機理中對系統(tǒng)主要參數(shù)值影響較小的基元反應(yīng)和組分。

2.1 敏感性分析法

敏感性分析分為局部敏感性分析和總體敏感性分析。局部敏感性主要用于研究基元反應(yīng)的動力學特性。局部敏感性分析正交表達式為

總體敏感性分析表達式為

式中：Bcj和Bfj分別為第j個反應(yīng)對N種組分生成速率和濃度的總體正交敏感性系數(shù)。本工作主要從整體研究部分基元反應(yīng)對主要計算結(jié)果影響程度，因此采用總體敏感性分析。

2.2 基元反應(yīng)的敏感性分析

爆震室預混物中m(H2):m(O2)=2:1，對二十步基元反應(yīng)進行敏感性分析，結(jié)果如圖2 所示。敏感性系數(shù)為正數(shù)代表此基元反應(yīng)對此組分起促進作用，負數(shù)代表此基元反應(yīng)對此組分進行消耗。圖2 顯示：基元反應(yīng)3 是影響整個系統(tǒng)反應(yīng)速率的關(guān)鍵基元反應(yīng)；對整個反應(yīng)系統(tǒng)影響較大的基元反應(yīng)序號為1～5、7、12～13 和16。

圖2 總體敏感性分析Fig.2 Overall sensitivity analysis

將通過敏感性分析所得出的九步基元反應(yīng)用于爆轟波數(shù)值模擬，計算所得的爆速、溫度與壓強與實驗值的誤差分別為0.236%、2.76%和4.57%，而爆轟波傳播1 μs 所需計算時間平均為50.556 min。

3 二十步基元反應(yīng)模型

3.1 爆轟波參數(shù)

表3 為網(wǎng)格為0.1 mm 二十步模型和九步模型數(shù)值計算結(jié)果與實驗、CJ 理論的比較。由表3 可知：九步簡化模型和二十步基元反應(yīng)模型計算所得的爆速、溫度和壓力與實驗值、理論結(jié)果較為一致。二十步模型數(shù)值計算所得的爆速、溫度與壓強與實驗值的誤差分別為0.2%、2.78%和2.74%，表明了本工作所采用的模型和數(shù)值計算方法的合理、正確性。

表3 基元反應(yīng)模型模擬爆轟波參數(shù)的對比Table 3 Comparison of simulated detonation parameters in elementary reaction models

相比于九步簡化模型，二十步模型數(shù)值計算時間略長。爆轟波傳播1 μs 所需計算時間平均為58.25 min。

3.2 爆轟波瞬態(tài)結(jié)構(gòu)

圖3 爆轟場壓力、溫度及OH質(zhì)量分數(shù)的分布云圖Fig.3 Distribution of detonation pressure, temperature and OH mass fraction

圖3給出了爆轟波傳播中，三波點在碰撞過程壓力、溫度和OH 離子質(zhì)量分數(shù)的分布情況，其中D 為未反應(yīng)氣囊，TP 為三波點，Hotspot 為熱點。數(shù)值模擬顯示，爆轟波傳播中具有典型的三波點結(jié)構(gòu)：入射激波(I)，橫波(TW)，馬赫干(M)。

兩個三波點碰撞前的爆轟流場(見圖3(a))。非平面的誘導激波波陣面是由入射激波和馬赫干交替組成，入射激波強度比馬赫干弱。OH 質(zhì)量分數(shù)場顯示：化學反應(yīng)區(qū)與馬赫干耦合在一起向前傳播，而與入射激波之間有一段明顯的距離，表現(xiàn)出化學反應(yīng)區(qū)與入射激波解耦的性質(zhì)。馬赫干、入射激波和橫波交匯處形成三波點，且爆轟場形成了多個未反應(yīng)的氣囊。

兩個三波點碰撞時的爆轟流場(見圖3(b))。三波點相互碰撞，形成熱點，導致該位置的溫度和壓力急劇升高。而此處的未反應(yīng)氣囊中，化學反應(yīng)速率急劇加快，釋放大量的能量，即為二次反應(yīng)。二次反應(yīng)釋放能量對爆轟波自持傳播起了重要作用。

三波點碰撞后的爆轟流場(見圖3(c))。由于三波點的二次反應(yīng)產(chǎn)生大量的能量，未反應(yīng)氣體和產(chǎn)物被擠壓形成向前的射流。由于不穩(wěn)定性，射流在其頭部形成蘑菇狀的結(jié)構(gòu)。先前的入射激波轉(zhuǎn)變?yōu)轳R赫干，原先的馬赫干退化成入射激波。強度較大的馬赫干與反應(yīng)區(qū)發(fā)生強耦合作用。

3.3 爆轟場中單個胞格的傳播性質(zhì)

爆轟波傳播中呈現(xiàn)典型的三波點結(jié)構(gòu)：馬赫干、入射激波、橫波。橫波與三波點周期性運動形成胞格結(jié)構(gòu)。圖4 給出了41.7 us 時刻爆轟場單個胞格的結(jié)構(gòu)圖。圖中顯示，爆轟波自持傳播并產(chǎn)生了周期性的爆轟胞格，在模型寬度方向上爆轟胞格不完全對稱。

圖4 單個胞格圖(單位: m)Fig.4 Single detonation cell pattern (unit: m)

圖5 為沿著單個胞格中心線的壓力分布，其中展示了3 個完整胞格結(jié)構(gòu)。三波點相互碰撞導致化學反應(yīng)速率增加，該位置的壓力和溫度升高，隨后中心線上壓力逐漸衰減，馬赫干衰退為入射激波，直至下一次三波點碰撞，壓力又達到峰值。計算初始點火位置形成的膨脹波，經(jīng)過壁面反射形成了如圖4 的反射波，波的相互作用導致流場中x軸上的0.053 m 至0.062 6 m之間區(qū)域壓力變化平緩，形成壓力平臺。隨著反射激波的發(fā)展，壓力逐漸衰減。

圖5 沿著單個胞格中心線的壓力圖Fig.5 Pressure profile along cell centerline versus distance

4 惰性氣體Ar 的稀釋作用

為了研究惰性氣體Ar 的稀釋作用，在相同氫氧混合氣體中添加質(zhì)量分數(shù)為70%的Ar。計算得到的穩(wěn)定爆轟波爆速，溫度和壓強分別為1 639.676 m/s, 2 933.018 K，1.958 MPa。將未添加與添加了惰性氣體的穩(wěn)定爆轟參數(shù)進行對比，結(jié)果如圖6 所示。可以看到：稀釋后的爆轟波參數(shù)，如爆轟壓力和溫度峰值分別降低0.047 MPa 和749.738 k。分析其原因，主要是添加了惰性氣體Ar 后，延緩了反應(yīng)的釋放，爆轟場溫度變化顯著。而壓強變化不大，主要是惰性氣體并未參與化學反應(yīng)所致。

圖6 Ar 稀釋對爆轟波參數(shù)的影響Fig.6 Influence of Ar dilution on detonation parameters

5 爆震室初始氣相壓強的影響

在真實環(huán)境中，初始氣相壓強不恒定，因此研究不同初始壓強對爆轟波參數(shù)的影響十分必要?？紤]氣體初始壓強分別為101 325、80 000 和60 000 Pa 的情況，圖7 為物質(zhì)的量之比為1:1 時的氫氧混合物爆轟參數(shù)的分布，數(shù)值計算結(jié)果如表4 所示?？梢钥吹剑罕?、溫度和壓強隨著初始壓強的降低而降低，即初始氣相壓強每下降1 Pa，爆轟波峰值壓力平均下降約18.51 Pa；而初始壓強對流場溫度變化影響不太明顯。

圖7 初始壓強對爆轟參數(shù)的影響Fig.7 Influence of initial pressure on detonation parameters

表4 初始壓強對爆轟波參數(shù)的影響Table 4 Influence of initial pressure on detonation parameters

6 結(jié)束語

脈沖爆震發(fā)動機爆震室氣相爆轟形成與發(fā)展是氣體流動和化學反應(yīng)同時發(fā)生且強烈耦合的過程。本工作建立九步簡化模型和二十步基元反應(yīng)模型，對氣相爆轟過程開展數(shù)值模擬，并將數(shù)值結(jié)果與實驗數(shù)據(jù)、理論計算結(jié)果進行對比和分析。

(1) 二十步基元反應(yīng)較為準確地反映流場中爆轟波精細結(jié)構(gòu)：馬赫干、入射激波、橫波.在爆轟波傳播中，入射激波與反應(yīng)區(qū)是分離的，而馬赫干與化學反應(yīng)區(qū)緊緊耦合在一起；三波點的碰撞導致化學反應(yīng)速率增加釋放的能量，為爆轟波持續(xù)傳播提供重要的支持。橫波和三波點周期性運動形成了爆轟波的胞格結(jié)構(gòu)。

(2) 化學反應(yīng)動力學影響著爆轟波的結(jié)構(gòu)和爆轟參數(shù)。通過總體敏感性分析，對二十步基元反應(yīng)簡化得出九步基元反應(yīng)。并將兩種基元反應(yīng)模型應(yīng)用于爆轟波模擬中，得出：基于二十步基元反應(yīng)模擬計算的爆速、溫度與壓強與實驗值的誤差較小，分別為0.2%、2.78%和2.74%。九步簡化模型的優(yōu)勢是計算相比詳細基元反應(yīng)，爆轟波傳播1 μs 所需的計算時間可節(jié)省13.21%。

(3) 爆轟波傳播過程中，惰性氣體Ar 延緩了化學反應(yīng)能量的釋放，使得爆轟波傳播速度和流場溫度降低。此外，氣相初始壓強的減少使得單位體積氣相燃料釋放的熱量降低，因此爆速、溫度和壓強也隨之降低。計算結(jié)果顯示，初始壓力每下降1 Pa，則爆壓峰值下降約18.51 Pa。Ar 稀釋和初始壓強的變化對爆轟波參數(shù)的影響和所得相關(guān)規(guī)律，為深入研究爆震發(fā)動機非均相爆轟特性奠定了理論基礎(chǔ)。