發動機爆震室氣相化學反應機理及爆轟波性質的數值研究

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

(上海大學上海市應用數學和力學研究所，上海 200072)

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

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

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

1 爆轟波形成和發展的控制方程和數值計算方法

1.1 控制方程

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

式中：u是守恒變量，Fi是對流通量，Q為源項，分別為

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

YK是第K個組分的質量分數，滿足

能量E定義為：

式中：hs為顯焓。

化學反應速率用Arrhenius 公式來計算

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

1.2 計算區域的網格及初始邊界條件

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

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

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

表1 網格尺寸對爆轟參數的影響Table 1 Influence of mesh size on detonation parameters

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

1.3 數值計算方法

1.3.1 爆轟的超聲速流場

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

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

1.3.2 化學反應動力學

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

表2 化學反應和反應參數Table 2 Chemical reactions and reaction parameters

2 九步基元反應模型

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

2.1 敏感性分析法

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

總體敏感性分析表達式為

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

2.2 基元反應的敏感性分析

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

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

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

3 二十步基元反應模型

3.1 爆轟波參數

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

表3 基元反應模型模擬爆轟波參數的對比Table 3 Comparison of simulated detonation parameters in elementary reaction models

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

3.2 爆轟波瞬態結構

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

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

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

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

三波點碰撞后的爆轟流場(見圖3(c))。由于三波點的二次反應產生大量的能量，未反應氣體和產物被擠壓形成向前的射流。由于不穩定性，射流在其頭部形成蘑菇狀的結構。先前的入射激波轉變為馬赫干，原先的馬赫干退化成入射激波。強度較大的馬赫干與反應區發生強耦合作用。

3.3 爆轟場中單個胞格的傳播性質

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

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

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

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

4 惰性氣體Ar 的稀釋作用

為了研究惰性氣體Ar 的稀釋作用，在相同氫氧混合氣體中添加質量分數為70%的Ar。計算得到的穩定爆轟波爆速，溫度和壓強分別為1 639.676 m/s, 2 933.018 K，1.958 MPa。將未添加與添加了惰性氣體的穩定爆轟參數進行對比，結果如圖6 所示?？梢钥吹剑合♂尯蟮谋Z波參數，如爆轟壓力和溫度峰值分別降低0.047 MPa 和749.738 k。分析其原因，主要是添加了惰性氣體Ar 后，延緩了反應的釋放，爆轟場溫度變化顯著。而壓強變化不大，主要是惰性氣體并未參與化學反應所致。

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

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

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

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

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

6 結束語

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

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

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

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