T型管內流動氣體中爆轟繞射過程的數值模擬*

潘振華，范寶春，歸明月

(1.江蘇大學能源與動力工程學院，江蘇 鎮江 230027；2.南京理工大學瞬態物理國家重點實驗室，江蘇 南京 210094)

爆轟波傳入突擴通道時，如翻越障礙物或進入分叉管等，將發生繞射。爆轟繞射現象是爆轟波傳播過程中經常出現的現象，與許多工業過程有關。爆轟繞射的演化過程涉及許多基本問題，如激波陣面和化學反應陣面的解耦、熄火、激波反射與相互作用、局部爆炸和二次起爆等，因此爆轟繞射現象引起廣泛關注。爆轟繞射可以發生在靜止介質中，也可以發生在流動介質中，兩者的演化過程不盡相同。但迄今，已經進行的爆轟繞射的研究，皆限于靜止系統。

I.B.Zeldovich等[1]通過實驗研究了爆轟繞射，提出爆轟管臨界管徑的概念，即其他參數一定的條件下，存在臨界管徑，當爆轟管直徑小于該值時，爆轟進入突擴空腔后會熄滅。此后，V.V.Mitrofanov等[2]發現，臨界管徑與爆轟胞格寬度有關。R.I.Soloukin等[3]通過實驗得到了爆轟繞射時激波與波后化學反應區分離以及熄滅后再生的紋影照片。S.B.Murray等[4]提出了爆轟波繞射后，熄滅再生的2種機制：(1) 反射誘導起爆機制；(2) 局部爆炸引起的自起爆機制。A.Smolinska等[5]通過實驗和計算，得到爆轟繞射后的清晰胞格結構。F.Pintgen等[6]通過實驗觀察得到了爆轟波在2種不同混合氣中的繞射過程，分別研究了爆轟波繞射后的3種傳播狀態。C.M.Guo等[7]和C.J.Wang等[8]對爆轟波在分叉管中的傳播過程進行了實驗和數值研究，其中涉及到了爆轟波衰減、馬赫反射、規則反射向馬赫反射的轉變以及爆轟波的重新起爆過程。

與靜止介質不同，流動介質中傳播的爆轟可分為迎風和順風2類。爆轟波繞射進入水平管后，激波面與高速來流作用，一側為迎風傳播，另一側為順風傳播。M.Y.Gui等[9]基于Euler方程數值模擬了尖劈誘導的斜爆轟的胞格結構，得到了迎風傳播的爆轟陣面隨時間和空間位置的變化特征。T.H.Yi等[10]通過數值模擬，討論了一維流動系統中的爆轟波，得出了迎風面的傳播速度小于順風面。潘振華等[11]通過數值模擬，討論了二維流動系統中的爆轟波的傳播特性，得到了在流動介質中，爆轟流場不對稱，爆轟壓力強度在上游、橫向和下游方向上依次遞減。K.Ishii等[12]通過實驗得到了爆轟波分別在600和910 m/s的高速預混氣流中傳播后留下的胞格軌跡，結果表明：向上游傳播的迎風爆轟波，爆速高于CJ值，胞格被壓縮；向下游傳播的順風爆轟波，爆速降低，胞格被拉長。

流動介質中的爆轟傳播問題，迄今尚無系統研究，特別是涉及到爆轟波的繞射、再起爆以及穩定傳播的問題。本文中以H2、O2、Ar體積比為2∶1∶1的預混氣體為研究對象，基于帶化學反應的二維Euler方程和五階精度的WENO(weighted essentially non-oscillatory)格式，對K.Ishii等的實驗中[12]爆轟波由T型管的垂直管進入通有超高速氣流的水平管中時產生的爆轟繞射及其后續演化過程進行數值研究，與相應靜止介質中的爆轟繞射進行比較，分析流動系統中爆轟繞射的流場變化特征和三波點煙箔軌跡的變化特征，探索爆轟再生的動力學機理。

1 數值研究方法

1.1 控制方程

假設混合氣體為理想氣體，忽略擴散、粘性和熱傳導。在貼體坐標系中，帶基元化學反應的多組分二維Euler方程為：

(1)

(2)

(3)

(4)

(5)

可燃預混氣體為體積比為2∶1∶1的H2、O2、Ar混合氣體，基元反應模型采用9組元和48個化學反應的詳細化學反應機理[13]，反應組元分別為H、O、H2、OH、H2O、O2、HO2、H2O2、Ar。

1.2 計算方法

計算過程中，為處理爆轟化學反應帶來的剛性問題，采用附加半隱的龍格-庫塔法(additive semi-implicit Runge-Kutta methods，ASIRK)[14]，該方法在時間上具有二階精度，且有很好的穩定性。空間項采用五階精度的WENO格式[15]計算。化學反應源項采用基于Gear格式的LSODE程序[16]計算。

1.3 計算模型

圖1分別為靜止系統和流動系統的計算域示意圖。圖1(a)為靜止系統，長31 cm、寬3.2 cm的水平燃燒室中心位置與長1.8 cm、寬3 cm的垂直爆轟管相連。圖1(b)為流動系統，燃燒室長31 cm，寬3.2 cm，爆轟管長1.8 cm、寬3 cm，兩者相連，爆轟管中心軸線在水平管9.6 cm位置上。整個計算區域采用均勻化網格，網格尺寸為0.1 mm×0.1 mm。其中燃燒室使用的網格數為3 100×320，爆轟管使用的網格數為180×300。

圖1 計算域示意圖Fig.1 The computational domain

對初壓為10.6 kPa、初溫為300 K的預混氣體進行前期數值計算，在垂直爆轟管中形成穩定自持傳播的爆轟波，管道內橫波數為7，對應胞格數穩定在3.5，以此初壓下自持胞格爆轟作為初始化計算。在其他區域中，仍充滿以等當量比的H2、O2為燃料、Ar作為稀釋劑的預混氣體，三者體積比為2:1:1。在靜止系統中，水平管內的預混氣體流速為零。在流動系統中，左側為進氣端，右側為出口端，預混氣以1 km/s的速度在燃燒室內從左向右流動。設靜止系統和流動系統2個算例初始壓力為10.6 kPa，溫度為300 K。計算中，均采用量綱一量，量綱一參考值為pf=15.4 kPa，Tf=300 K，L0=0.1 m。邊界處理中，CDJ、KEF和BA用無催化、絕熱固壁邊界條件；邊界AC和BF用外推邊界條件；JK用插值外流邊界條件。

2 數值驗證

為了驗證計算的準確性，對E.S.Oran等[17]的算例進行數值模擬，與E.S.Oran等[17]的計算結果進行對比，并將計算得到的爆轟波速與C.A.Eckett[18]的實驗結果進行比較。

圖2 直管中的計算胞格Fig.2 Numerical cells structure in straight tube

為了直接與文獻[17]計算結果進行對比，數值驗證采用與文獻[17]相同的計算域和初始條件，即在直管中，充滿H2、O2、Ar的預混氣體，三者體積比為2∶1∶7，初始壓力和溫度分別為6.67 kPa和298 K。圖2為本文數值模擬得到的胞格圖像，表1為本文計算結果與文獻[17]的對比，其中：Δx、Δy分別為單個網格橫向和縱向的尺寸，va為爆轟波的平均速度，L和W分別為胞格橫向和縱向的尺寸。計算得到的胞格尺寸為約54 mm×30 mm，與E.S.Oran在相同條件下得到的胞格尺寸54 mm×31 mm[17]基本一致；胞格縱橫比為0.556，與文獻[18]中實驗得到的縱橫比0.54符合較好。本文中爆轟波的傳播速度為1 589 m/s，與Gordon-McBride程序算出的理論爆轟波CJ爆速(1 617 m/s)相比誤差在小于2%，與文獻[18]實驗得到的爆轟波速度(1 550 m/s)基本一致。因此，本文中的數值模擬能夠重復文獻[17]的計算，并且計算結果與文獻[18]的實驗結果相符。

表1 數值驗證結果與文獻[17]計算結果的對比Table 1 Comparison between current simulation and referance [17]

3 結果與討論

3.1 爆轟繞射的波陣面結構

圖3 t=16.2 μs時靜止系統爆轟波紋影圖Fig.3 Numerical schlieren of detonation in a quiescent system while t=16.2 μs

圖3為t=16.2 μs時，靜止系統中爆轟繞射流場的計算紋影圖。其中右圖為左圖的局部放大顯示。此時，爆轟陣面尚未抵達水平管的上壁面。右側局部放大圖中，實線為激波陣面，虛線為化學反應陣面。由圖可知，在未受擾動區域，激波與反應陣面耦合很好，平面形狀的爆轟波以固有的速度向前發展。但在被擾動區域，陣面彎曲，激波已與反應陣面解耦，兩者之間存在一個經前導激波預壓的高溫高壓未燃氣體區域。

圖4為t=16.2 μs時，流動系統中爆轟繞射流場的計算紋影圖。其中左圖為右圖的局部放大顯示。對于右側的順風繞射，被擾動的波陣面上，激波與反應陣面解耦，大致情形與靜止系統的相仿。但對于左側的迎風繞射，從垂直管中傳出的爆轟，其波后爆轟產物對水平管內的定向來流而言，具有與尖劈類似的壓縮作用，這使得迎風繞射的爆轟陣面具有斜爆轟的某些特點。穩定的斜爆轟由一系列子波構成，其精細結構如圖5所示[9-11]。由圖5可知，激波陣面的彎曲，使陣面分為S1和S2兩部分，并在三波點P處碰撞，形成橫波T1，并誘發化學反應，也進而使T1和S2在三波點附近的部分成為爆轟波D1，T1稱為橫向爆轟波。與正爆轟的精細結構相仿，此類子波結構(精細結構)使得斜爆轟得以穩定自持。圖4左側的局部放大圖中同樣出現了穩定斜爆轟所具有的子波結構(如三波點P1和P2)。

圖4 t=16.2 μs時流動系統中的計算紋影圖Fig.4 Numerical schlieren of detonation in a flowing system while t=16.2 μs

圖5 迎風面上波系結構示意圖Fig.5 Schematic diagram of wave structure on windward side

3.2 爆轟繞射的再生機理

當爆轟波抵達水平管的上壁面時，將在壁面反射。圖6給出了流通系統中爆轟波在上壁面反射后不同時刻流場的計算紋影、壓力分布、溫度分布和H組分的質量分數分布圖。

圖6 流動系統中爆轟波的再生及波系結構的演變Fig.6 Reinitiation event of detonation and evolution of wave structure in a flowing system

激波在壁面的反射包括2部分：未受擾動的平面爆轟在壁面的正反射R1和隨后的被擾動的彎曲爆轟在壁面的馬赫反射R2。反射后，在壁面附近形成向兩端傳播的馬赫干M和向下方傳播的橫向反射波R(圖6(a)和圖6(b)紋影圖中的水平白線)。反射波波后為高溫高壓區域(見圖6的壓力和溫度圖)。由于馬赫干的強度較高，故可誘導爆轟波。而橫向傳播的反射波，在傳播過程中，有部分陣面掃過因激波陣面與反應區解耦而形成的預壓未燃區，從而使之點火，形成橫向傳播的爆轟波T。因此，整個反射波陣面是爆轟波和激波組成的復合波陣面，兩端為彎曲的爆轟波T，中間由激波R連接。在反射波陣面向下傳播的同時，左側類斜爆轟陣面上的三波點P3和P4也作相對運動，直至碰撞(如圖6(a)～(c)所示)，這與斜爆轟精細結構的情形一樣。左側的反射波(爆轟子波)在與相對穩定的類斜爆轟陣面的相互作用過程中向左傳播并向下擴展，陣面趨于平整并與水平方向垂直。而右側反射波(強度較弱的爆轟子波)在與不斷衰減的激波的相互作用過程中向右傳播并向下擴展，其強度也不斷衰減，從而形成傾斜的波陣面。在形成斜爆轟的過程中，激波陣面與火焰陣面未很好耦合(如圖6(a)～(c)所示)。由圖6～9還可看出，雖然迎風爆轟的強度大于順風爆轟，但順風爆轟的傳波速度大于迎風爆轟。此后，順風斜爆轟在下壁面反射，由于反射爆轟的強度大于入射爆轟，故隨著反射爆轟陣面的擴展，爆轟得以加強(如圖6(d)所示)。經過上下壁面的若干次反射，右側爆轟波也可逐漸發展為與水平方向垂直的正爆轟。至此，爆轟波經歷了繞射導致的熄火和二次起爆，最終成為沿水平管，迎風和順風向兩端傳播的正爆轟波。靜止系統中爆轟波繞射過程與流動系統順風側的情形類似。

3.3 爆轟繞射過程的胞格結構

從精細流場角度講，穩定傳播的爆轟陣面的強度是不均勻的。該陣面是由一系列馬赫碰撞事件構成的，三波點上的壓力最高。三波點(即高壓點)的軌跡通常具有圖2所示的網狀圖形，稱為爆轟胞格。圖7為靜止系統中爆轟在T型管中繞射時胞格結構演變圖。爆轟在垂直管中穩定傳播時，形成圖7中a1區的規則胞格。爆轟進入水平管繞射時，部分陣面因擾動而熄火，b1區域為爆轟波未受擾動區域，胞格形狀不變，與a1區一致。c1區為“爆轟熄火”區，胞格消失。b1區與c1區的邊界為傾斜直線，傾斜角等于擾動角。繞射爆轟抵達上壁面后，在壁面反射，反射波傳播區域中，僅有強度最大的橫向爆轟留下傳播軌跡，即圖中d1區域。反射波在下壁面再次反射，形成向兩端傳播，強度不大的斜爆轟波。波陣面上，馬赫碰撞的三波點單向運動，形成e1區平行的三波點軌跡線。反射波在上下壁面反復反射，強度不斷增加，最終會重新成長為穩定傳播的平面爆轟，胞格也重新形成f1區中規則的網狀結構。

圖8為流動系統中爆轟繞射的胞格結構演變圖。受水平管中定向流動的影響，未受擾動區b2向右側偏移，胞格也發生傾斜。在右側(順風側)，胞格圖像與圖7中大致相同，但各區的面積被拉長。此外，橫向爆轟的軌跡出現兩次(圖8中e2區域)，而靜止系統只有一次，這說明靜止系統中，繞射爆轟的增長速度更快。在左側(迎風側)，爆轟繞射時，基本未出現熄滅陣面(熄火僅出現在拐角附近很小的c2區域內)，繞射陣面上三波點以及發源于此的橫向爆轟的運行軌跡如圖中d2區所示。反射波抵達下壁面時，已基本成為正爆轟，左側g2區胞格尺寸被顯著壓縮，右側g2區胞格尺寸被顯著拉長。

圖8 流動系統中的胞格結構演變Fig.8 Cellular structure in a flowing system

4 結 論

對比了靜止系統和流動系統中爆轟波從繞射到再生的整個過程，得到如下結論。(1) 爆轟波在流動系統中發生繞射，順風傳播的爆轟波趨于熄火；在迎風側，三波點附近局部區域內激波陣面始終與波后化學反應區耦合在一起，形成斜爆轟結構，使爆轟得以維持。(2) 由于兩側爆轟傳播速度不同，迎風側的胞格圖像明顯被壓縮，順風側被拉長。(3) 在爆轟波的傳播過程中，波陣面上橫向爆轟波的產生對爆轟波的起爆過程起到了關鍵作用。

[1] Zeldovich I B, Kogarko S M, Simonov N N.An experimental investigation of spherical detonation in gases[J].Soviet Physics: Technical Physics, 1956,1(8):1689-1713.

[2] Mitrofanov V V, Soloukhin R I.The diff'raction of mutifront detonation waves[J].Soviet Physics: Doklady, 1965,9(12):1055-1058.

[3] Soloukhin R I, Ragland K W.Ignition processes in expanding detonations[J].Combustion and Flame, 1965,13(3):295-302.

[4] Murray S B, Lee J H.On the transformation of planar detonation to cylindrical detonation[J].Combustion and Flame, 1983,52:269-289.

[5] Smolinska A, Khasainov B, Virot V, ea al.Detonation diffraction from tube to space via frontal obstacle[C]∥Proceedings of the Fourth European Combustion Meeting.Vienna, Austria, 2009.

[6] Pintgen F, Shepherd J E.Detonation diffraction in gases[J].Combustion and Flame, 2009,156(3):665-677.

[7] Guo C M, Wang C J, Xu S L, ea al.Cellular pattern evolution in gaseous detonation diffraction in a 90°-branched channel[J].Combustion and Flame, 2007,148(3):89-99.

[8] Wang C J, Xu S L, Guo C M.Gaseous detonation propagation in a bifurcated tube[J].Journal of Fluid Mechanics, 2008,599:81-110.

[9] Gui M Y, Fan B C.Wavelet structure of wedge-induced oblique detonation waves[J].Combustion Science and Technology, 2012,184(10/11):1456-1470.

[10] Yi T H, Wilson D R, Lu F K.Numerical study of unsteady detonation wave propagation in a supersonic combustion chamber[R].Arlington, Texas, USA: University of Texas at Arlington, 2004.

[11] 潘振華，范寶春，歸明月，等.流動系統中爆轟波傳播特性的數值模擬[J].爆炸與沖擊,2010，30(6)：593-597.Pan Zhen-hua, Fan Bao-chun, Gui Ming-yue, et al.Numerical study of detonation wave propagation in a flow system[J].Explosion and Shock Waves, 2010,30(6):593-597.

[12] Ishii K, Kataoka H, Kojima T.Initiation and propagation of detonation waves in combustible high speed flows[J].Proceedings of Combustion Institute, 2009,32(2):2323-2330.

[13] Oran E S, Young T R, Boris J P, et al.Weak and strong ignition I: Numerical simulation of shock tube experiments[J].Combustion and Flame, 1982,48:135-148.

[14] Zhong X L.Additive semi-implicit Runge-Kutta methods for computing high-speed nonequilibrium reactive flows[J].Journal of Computational Physics, 1996,128(1):19-31.

[15] Jiang G S, Shu C W.Efficient implementation of weighted ENO schemes[J].Journal of Computational Physics, 1996,126(2):202-228.

[16] Radhakrishnan K, Hindmarsh A C.Description and use of LSODE, the livermore solver for ordinary differential equations, UCRL-ID-113855[R].Washington, USA: NASA, 1993.

[17] Oran E S, Weber J W, Stefaniw E I, et al.Numerical study of a two-dimensional H2-O2-Ar detonation using a detailed chemical reaction model[J].Combustion and Flame, 1998,113(1/2):147-163.

[18] Eckett C A.Numerical and analytical studies of the dynamic of gaseous detonation[D].CA, USA: California Institute of Technology, 2001.