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

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

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

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

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

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

1.1 控制方程

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

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

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

YK是第K個(gè)組分的質(zhì)量分?jǐn)?shù)，滿足

能量E定義為：

式中：hs為顯焓。

化學(xué)反應(yīng)速率用Arrhenius 公式來(lái)計(jì)算

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

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

爆震室直徑為2 cm，充滿初始?jí)毫橐粋€(gè)標(biāo)準(zhǔn)大氣壓、初始溫度為300 K 的氫氣和氧氣混合氣體，其中氣體組分的物質(zhì)的量之比為m(H2):m(O2)=2:1。爆震室左端點(diǎn)火，示意圖如圖1 所示。

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

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

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

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

1.3 數(shù)值計(jì)算方法

1.3.1 爆轟的超聲速流場(chǎng)

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

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

1.3.2 化學(xué)反應(yīng)動(dòng)力學(xué)

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

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

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

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

2.1 敏感性分析法

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

總體敏感性分析表達(dá)式為

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

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

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

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

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

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

3.1 爆轟波參數(shù)

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

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

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

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

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

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

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

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

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

3.3 爆轟場(chǎng)中單個(gè)胞格的傳播性質(zhì)

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

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

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

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

4 惰性氣體Ar 的稀釋作用

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

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

5 爆震室初始?xì)庀鄩簭?qiáng)的影響

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

圖7 初始?jí)簭?qiáng)對(duì)爆轟參數(shù)的影響Fig.7 Influence of initial pressure on detonation parameters

表4 初始?jí)簭?qiáng)對(duì)爆轟波參數(shù)的影響Table 4 Influence of initial pressure on detonation parameters

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

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

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

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

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