內(nèi)置變徑障礙物組管道內(nèi)爆轟波傳播規(guī)律的研究*

丁建旭，杜群貴，王新華，吳雨蒙

(1.廣州特種機電設(shè)備檢測研究院 國家防爆設(shè)備質(zhì)量監(jiān)督檢驗中心，廣東 廣州 510760；2. 華南理工大學(xué) 機械與汽車工程學(xué)院，廣東 廣州 510006；3. 廣東工業(yè)大學(xué) 機電工程學(xué)院，廣東 廣州510006)

0 引言

典型?；返谋?zāi)害經(jīng)常涉及爆轟波與復(fù)雜障礙物的強耦合作用，造成如反射和繞射等強非線性現(xiàn)象的發(fā)生，形成梯度較大的復(fù)雜爆炸環(huán)境，對人員及建筑設(shè)備造成嚴重的傷害和破壞。因此，研究爆轟波與管道內(nèi)復(fù)雜障礙物的耦合作用機制，對工業(yè)重大爆炸災(zāi)害事故的預(yù)防及快速響應(yīng)能力的提升均有著重要的現(xiàn)實意義。

隨著描述爆轟的數(shù)學(xué)物理模型的不斷完善，采用數(shù)值模擬方法研究爆轟波傳播動力學(xué)行為已成為新的發(fā)展方向，受到了許多學(xué)者的關(guān)注[1-3]。張寶亮等[4]通過變化障礙物的高度及間距，數(shù)值研究了矩形障礙物對甲烷/氧氣混合物的氣相爆轟波傳播過程的影響；隨著數(shù)值技術(shù)的不斷發(fā)展，Wang等[5]采用高精度WENO格式研究了一組矩形障礙物對爆轟波傳播規(guī)律的影響；真實的爆炸環(huán)境往往是復(fù)雜的，文獻[6-7]在無約束空間內(nèi)研究了一簇方形和圓柱形障礙物對火焰加速以及爆燃轉(zhuǎn)爆轟過程的影響機制；文獻[8]采用商業(yè)計算軟件對爆炸沖擊波在建筑群中的傳播規(guī)律進行了模擬，表明數(shù)值模擬技術(shù)在解決實際問題時具有一定的優(yōu)越性。

從以上的討論可以看出，許多學(xué)者為了簡化模型和研究的方便，基本上采用了規(guī)則的如矩形等障礙物作為火焰?zhèn)鞑サ臄_動對象。但是，在實際工程應(yīng)用中所涉及的障礙物往往更加復(fù)雜。管道中變尺寸圓柱障礙物組對爆轟波傳播過程的影響機制，目前相關(guān)的研究鮮有報道，為此，本文采用數(shù)值手段開展爆轟波在內(nèi)置一組變徑圓柱障礙物管道中的傳播規(guī)律研究，以揭示預(yù)混氣體爆轟波在復(fù)雜障礙物組作用下發(fā)生的復(fù)雜流場改變及爆炸特征參數(shù)的時空分布規(guī)律。

1 控制方程組及數(shù)值方法

采用兩步化學(xué)反應(yīng)[9-10]描述爆轟波與一組變徑圓柱障礙物的相互作用過程，相應(yīng)無量綱守恒型控制方程組的具體形式為：

(1)

式中：

U=[ρρuρvρEραρβ]T

F=[ρuρu2+pρuvρu(e+p/ρ)ρuαρuβ]T

G=[ρvρuvρv2+pρv(e+p/ρ)ρvαρvβ]T

S=[0 0 0ρωαρωβ]T

p=ρT

(2)

(3)

(4)

(5)

式中：q為放熱量；R為氣體常數(shù)；r為密度；T為溫度；P為壓強；U和v分別表示x和y方向的速度分量；E1和E2為活化能；k1和k2為反應(yīng)速率常數(shù)；r是比熱率；a和b為控制反應(yīng)進行的變量，當(dāng)0

控制方程組的通量采用五階精度的WENO格式[11]進行離散，時間推進采用三階精度的TVD Runge-Kutta格式[12]。根據(jù)CLF條件確定時間步長：

(6)

式中：CFL取0.6；dx和dy為網(wǎng)格尺寸；c為聲速。

爆轟波在內(nèi)置復(fù)雜障礙物管道中的傳播過程，特別在邊界附近的數(shù)值流場強烈地依賴于邊界的數(shù)值處理。因此，文中固壁邊界均采用高精度ILW方法[13-15]進行處理。

2 計算工況及數(shù)值結(jié)果討論

2.1 工況設(shè)置

在實際化工生產(chǎn)中，圓形輸運管道的管徑及其布局方式往往是復(fù)雜的。若發(fā)生爆炸事故，不同直徑管道的布局方式對爆炸場的影響規(guī)律值得深入探討，可為有限約束空間內(nèi)輸運管道的布局優(yōu)化以及對局部結(jié)構(gòu)強化抗爆安全提供指導(dǎo)。矩形無量綱管道的長度為90，寬為15，內(nèi)置5個可變半徑的圓柱形障礙物組。在計算中障礙物的個數(shù)保持不變，均為5個。工況1的圓柱形障礙物半徑均為3，工況2中的障礙物半徑從1逐漸增加到5，工況3的圓柱形障礙物半徑的變化規(guī)律與工況2正好相反，由大到小變化。在各個圓柱形障礙物的5個位置設(shè)置監(jiān)測點Mijk，記錄壓力的變化，其中下標(biāo)i代表工況，j代表從左到右的第幾個圓柱形障礙物，k代表每個圓柱形障礙物的5個監(jiān)測點，從1到5代表點(xk0+rkcos(qk),yk0+rksin(qk))，點(xk0,yk0)代表的是5個障礙物的圓心，分別為(15.0,7.5)，(30.0,7.5)，(45.0,7.5)，(60.0,7.5) 和 (75.0,7.5)，3種工況圓心位置均相同，rk為各個圓柱的半徑，qk從180°以45°減幅逐漸過渡到0°，共5個角度。采用ZND條件作為入流條件，右邊為自由出流邊界，管道上、下表面以及5個圓柱障礙物組的表面均采用固壁邊界條件，采用高精度ILW邊界處理方法進行處理，保證固壁邊界附近與管道內(nèi)部爆轟流場均為五階精度，提高數(shù)值分辨率。無量綱的初始條件為(ρ,u,v,p,α,β)=(1,0,0,0.16,1,1). 網(wǎng)格大小為dx=dy=0.017，保證了數(shù)值方法的收斂性[13]。

2.2 數(shù)值結(jié)果及討論

為了節(jié)約篇幅，文中僅以圓柱障礙物半徑逐漸增加的工況2和圓柱障礙物半徑逐漸減小的工況3的密度梯度云圖為例，考察爆炸場內(nèi)復(fù)雜波系的演化過程，工況1條件下的爆轟波與5個等半徑圓柱障礙物作用的密度梯度云圖不再贅述。

圖1顯示5個障礙物半徑從1漸增到5爆轟波傳播過程中的密度梯度變化。從圖1(a)可以看出，由于障礙物1的半徑較小，對平面爆轟波陣面的破壞區(qū)域有限，形成的反射波經(jīng)過長時間的演變，抵達管道上下壁面時強度損失較多。經(jīng)上下壁面碰撞形成的反射波強度也較弱，繞射波陣面從三波點到對稱軸兩側(cè)的密度梯度線開口逐漸增大，如圖1(b)所示。從圖1(c)可以看出，在波陣面抵達圓柱形障礙物2之前，波陣面上存在的兩對相向運動的三波點即將發(fā)生碰撞。如圖1(d)所示，向下運動的三波點遭到破壞，而向上運動的三波點增強，且在圓柱形障礙物2左側(cè)發(fā)生了反射現(xiàn)象。圖1(e)表明，向上運動的三波點經(jīng)管道壁面反射后強度增強，向著對稱軸方向運動，此時繞射波發(fā)生碰撞，同時，繞射陣面成弧形彎曲，但沒有明顯的局部爆炸產(chǎn)生。隨著反射波的持續(xù)相互作用以及壓縮波對未反應(yīng)氣體的壓縮活化，在對稱軸附近的流場內(nèi)部發(fā)生了局部爆炸，釋放大量的化學(xué)能，提升了局部的壓力和溫度，如圖1(f)所示。從圖1(g)中可以看出，局部爆炸形成的射流頭部加速傳播，并與障礙物3左側(cè)首先發(fā)生碰撞并反射，形成反射波。射流頭部橫向運動的橫波沿著波陣面繞障礙物固壁面運動，增強了圓柱壁面附近波陣面的傳播速度，但是，隨著傳播的繼續(xù)，三波點不斷向管道壁面移動，導(dǎo)致稀疏效應(yīng)逐漸顯現(xiàn)，三波點以下陣面結(jié)構(gòu)呈弧形彎曲狀，如圖1(h)所示。隨后，出現(xiàn)了與圖1(f)類似的爆轟波陣面結(jié)構(gòu)，如圖1(i)所示。隨著爆轟波傳播路徑上障礙半徑的逐漸增大，繞射爆轟波陣面呈弧形相交結(jié)構(gòu)，在相交軸線附近不再出現(xiàn)內(nèi)凹結(jié)構(gòu)，如圖1(k)和圖1(l)，密度梯度強度分布出現(xiàn)變化，主要集中在局部爆炸射流的頭部，如圖1(l)所示，此時波陣面沒有出現(xiàn)明顯的三波點結(jié)構(gòu)。

圖1 工況2條件下爆轟波傳播過程的密度梯度Fig.1 Density gradient contours in the propagation of detonation wave for case 2

圖2表示在五個障礙物半徑從5逐漸遞減到1的條件下爆轟波傳播過程中的密度梯度變化云圖。與工況2類似，爆轟波與圓柱障礙物組相互作用產(chǎn)生的連續(xù)反射和繞射是爆轟波動力學(xué)演化的主要機制。但是由于圓柱障礙物半徑的逐漸減小，從圖2(l)中可以看出，在爆轟波與障礙物5作用之前，波陣面上已經(jīng)出現(xiàn)了多個三波點結(jié)構(gòu)。從圖2(m)可以看出，障礙物5右側(cè)靠近軸線附近的三波點個數(shù)增多，其他位置三波點能量增強。由于爆轟波在水平傳播方向上自由通道漸寬，反射和繞射強度減小，隨著爆轟波的繼續(xù)傳播，波陣面的三波點大量出現(xiàn)，如圖2(n)所示，釋放的能量維持爆轟波繼續(xù)向下游傳播。

圖2 工況3條件下爆轟波傳播過程的密度梯度Fig.2 Density gradient contours in the propagation of detonation wave for case 3

通過以上分析可以發(fā)現(xiàn)，不同的變徑障礙物組形成了不同工況條件下阻礙物阻塞比漸變趨勢的變化，造成了爆轟波在水平傳播方向上約束程度的差異，以及障礙物對爆轟波反射和繞射強度的改變。這些耦合效應(yīng)使得較大阻塞比管道段圓柱障礙物以及管道壁面反射和繞射碰撞產(chǎn)生的局部爆炸是爆轟波維持傳播的重要原因，而在較小阻塞比管道段，波陣面演化產(chǎn)生的多個三波點提供的大量的化學(xué)能是此段管道中爆轟波自持傳播的主要機制。

圖3(a)～(c)分別表示在工況1～3條件下5個障礙物的5個監(jiān)測點的最大壓力。圖中p*是無量綱參考壓力[13]。在工況1條件下，爆轟波傳播通道中存在5個等半徑圓柱障礙物。從圖3(a)中可以看出，在爆轟波反射區(qū)監(jiān)測點M111和M112的最大壓力均為2.1，由于繞射區(qū)的稀疏效應(yīng)，使得監(jiān)測點M114和M115的最大壓力均有所降低，約為1.3，在反射結(jié)束與繞射開始的臨界區(qū)的監(jiān)測點M113的最大壓力為1.6，處于兩者之間，表明反射對爆轟波的最大壓力提升具有正反饋效應(yīng)，而繞射會使最大壓力顯著地降低。隨著爆轟波的繼續(xù)傳播，與更多的圓柱形障礙物繼續(xù)發(fā)生相互作用，使得爆轟流場內(nèi)部波系更加復(fù)雜，同一個圓柱障礙物壁面處的監(jiān)測點的最大壓力不再具有類似的分布規(guī)律。但是各個圓柱形障礙物的正反射區(qū)域的監(jiān)測點的最大壓力顯著高于其他監(jiān)測點，特別是監(jiān)測點M141的最大壓力躍升至10.0，為所有監(jiān)測點的時空最大壓力，且M111的最大壓力值遠低于其他該工況下正反射處監(jiān)測點所獲得的最大壓力。在繞射區(qū)的45°方向上監(jiān)測點M114，M124，M134，M144和M154的最大壓力基本沒有顯著升高或降低。對于爆轟波傳播通道中圓柱形障礙物半徑漸增的工況2，如圖3(b)所示，在第一個半徑為1的障礙物壁面上獲得不同作用區(qū)域的最大壓力出現(xiàn)明顯的變化，正反射的監(jiān)測點M211的最大壓力為8.0，遠大于該圓柱障礙物壁面其他位置監(jiān)測點獲得的最大壓力。在工況2條件下，所有障礙物正反射區(qū)域的監(jiān)測點的最大壓力遠大于對應(yīng)圓柱障礙物其他監(jiān)測點的最大壓力，監(jiān)測點M241處的最大壓力為時空最大壓力，達到11.4，但是最后一個圓柱障礙物處監(jiān)測點M251的最大壓力僅為6.2，小于其他任何障礙物正反射區(qū)監(jiān)測點的最大壓力。在各個圓柱形障礙物的繞射區(qū)的45°方向上監(jiān)測點的最大壓力雖然有降低的趨勢，但不是非常顯著，位于其他相同監(jiān)測點位置但屬不同圓柱障礙物的最大壓力均有一定程度的改變。對于圓柱形障礙物半徑逐漸減小的工況3，如圖3(c)所示，第一個圓柱形障礙物的5個監(jiān)測點的最大壓力歷史位于監(jiān)測點M312，達到5.7，正反射區(qū)監(jiān)測點M311的最大壓力甚至小于監(jiān)測點M315，僅為2.9。在各個圓柱形障礙物的繞射區(qū)45°方向上監(jiān)測點的最大壓力變化不顯著，但是對于各個圓柱形障礙物對應(yīng)的在135°方向上的監(jiān)測點出現(xiàn)了較為明顯的差異。

圖3 圓柱形障礙物監(jiān)測點的最大壓力Fig.3 The maximum pressures at the monitoring points around cylindrical obstacles

從以上的分析可以看出，變徑圓柱障礙物組對爆轟波傳播具有重要的影響，特別在障礙物壁面處，總體上正反射區(qū)獲得的最大壓力均較大(除監(jiān)測點M311)，其他區(qū)域特別是45°方向上監(jiān)測點的最大壓力基本沒有出現(xiàn)顯著的起伏變化，維持在較低值。就監(jiān)測點處最大壓力的空間分布而言，3種工況的最大壓力的最大值均出現(xiàn)在第4個圓柱障礙物的正反射區(qū)的監(jiān)測點處，且工況2在該點處獲得的最大壓力均大于障礙物組半徑維持不變或者逐漸減小的2種工況。

圖4為3種工況條件下中心對稱軸線上最大壓力隨時間的變化曲線，圖中t*是無量綱參考時間[13]。爆轟波繞過障礙物形成的稀疏波從管道上下兩個方向朝管道對稱中心匯聚并發(fā)生碰撞，在局部區(qū)域釋放大量的化學(xué)能，導(dǎo)致溫度和壓力的提升，同時，反射波在圓柱障礙物之間以及與管道上下壁面的不斷反射增加了軸線附近的流場復(fù)雜度，從而反映于軸線上最大壓力的周期性變化。對比三條曲線可以發(fā)現(xiàn)，對于相同障礙物半徑的工況1，最大壓力的第3個峰值為16.7，均大于其他工況的最大壓力的最大峰值，工況3的最大壓力的最大峰值與其他工況最大壓力的最大峰值相比最小。工況2的最大壓力的最大峰值比工況1和工況3均提前出現(xiàn)。

圖4 對稱軸線上最大壓力歷史Fig.4 The maximum pressure history along the central axis of tube

以上分析表明，變半徑障礙物組的不同半徑變化設(shè)置對中心軸線上的最大壓力的最大峰值具有重要的影響，等半徑障礙物組的工況1條件下獲得的最大壓力的最大峰值比障礙物半徑逐漸增加的工況2或遞減的工況3的最大壓力的最大峰值都大。

3 結(jié)論

1)較大阻塞比管道段圓柱障礙物以及管道壁面反射和繞射碰撞產(chǎn)生的局部爆炸是爆轟波維持傳播的重要原因，而在較小阻塞比管道段，波陣面演化產(chǎn)生的多個三波點提供的大量的化學(xué)能是此段管道中爆轟波自持傳播的主要機制。

2)在障礙物壁面處，正反射區(qū)的最大壓力較大，45°方向上的最大壓力維持在較低值。3種工況的最大壓力的最大峰值均出現(xiàn)在第4個圓柱障礙物的正反射區(qū)，且逐漸增大障礙物半徑的工況2的監(jiān)測點最大壓力的最大峰值為所有工況最大。

3)不同障礙物半徑變化規(guī)律對中心軸線上的最大壓力的最大峰值具有重要影響，等半徑障礙物工況下獲得的最大壓力的最大峰值比障礙物半徑逐漸遞減或增加的最大壓力的最大峰值都大。

[1] 王公忠, 張建華, 李登科, 等. 障礙物對預(yù)混火焰特性影響的大渦數(shù)值模擬[J]. 爆炸與沖擊, 2017, 37(1):69-75.

WANG Gongzhong, ZHANG Jianhua, LI Dengke, et al. Large eddy simulation of impacted obstacles’effects on premixed flame’s characteristics [J]. Explosion and Shock Waves, 2017, 37(1):69-75.

[2] 徐景德, 楊庚宇. 置障條件下的礦井瓦斯爆炸傳播過程的數(shù)值模擬研究[J]. 煤炭學(xué)報, 2004, 29(1):53-56.

XU Jingde, YANG Gengyu. Numerical simulation of the barricade encouraging effect in the process of gas explosion propagation [J]. Journal of China Coal Society, 2004, 29(1):53-56.

[3] 陳志華, 葉經(jīng)方, 范寶春, 等. 方形管內(nèi)楔形障礙物對火焰結(jié)構(gòu)與傳播的影響[J]. 爆炸與沖擊, 2006, 26(3):208-213.

CHEN Zhihua, YE Jingfang, FAN Baochun, et al. Effects of a wedge obstacle on flame propagation and its structure [J]. Explosion and Shock Waves, 2006, 26(3):208-213.

[4] 王可強, 蘇經(jīng)宇, 王志濤. 爆炸沖擊波在建筑群眾傳播規(guī)律的數(shù)值模擬研究[J]. 中國安全科學(xué)學(xué)報, 2007, 17(10):121-127.

WANG Keqiang, SU Jingyu, WANG Zhitao. Numerical investigation of propagation rules of blast shock wave in building cluster [J]. China Safety Science Journal, 2007, 17(10):121-127.

[5] WANG C, MA T B, LU J. Influence of obstacle disturbance in a duct on explosion characteristics of coal gas [J]. Physics, Mechanics & Astronomy, 2010, 53(2):269-278.

[6] OGAWA T, ORAN E S, GAMEZO V N. Numerical study on flame acceleration and DDT in an inclined array of cylinders using an AMR technique [J]. Computers & Fluids, 2016(85):63-70.

[7] OGAWA T, GAMEZO V N, ORAN E S. Flame acceleration and transition to detonation in an array of square obstacles [J]. Journal of Loss Prevention in the Process Industries, 2013(26):355-362.

[8] 張寶亮, 丁玨, 王慶濤, 等. 約束空間可燃氣體燃燒爆轟特性的數(shù)值研究[J]. 中國安全生產(chǎn)科學(xué)技術(shù), 2012,8(8):23-27.

ZHANG Baoliang, DING Jue, WANG Qingtao, et al. Numerical study on the combustion and detonation characteristics of combustible gas in constraint space [J]. Journal of Safety Science and Technology, 2012,8(8):23-27.

[9] 王剛, 張德良, 劉凱欣. 氣相爆轟波數(shù)值模擬中化學(xué)反應(yīng)模型研究[J]. 高壓物理學(xué)報, 2008(4): 350-356.

WANG Gang, ZHANG Deliang, LIU Kaixin. Study on chemical reaction models in gaseous detonation numerical simulation [J]. Chinese Journal of High Pressure Physics, 2008(4):350-356.

[10] 王剛, 張德良, 劉凱欣. 對數(shù)值模擬氣相爆轟中二階段化學(xué)反應(yīng)模型再研究[J]. 空氣動力學(xué)報，2009, 27(3):363-368.

WANG Gang, ZHANG Deliang, LIU Kaixin. Re-study on the two-step chemical reaction model in numerical simulation of gaseous detonation [J]. Acta Aerodynamica Sinica, 2009, 27(3):363-368.

[11] SHU C W, OSHER S. Efficient implementation of essentially non-oscillatory shock-capturing schemes [J]. Journal of Computational Physics, 1988, 77(2):439-471.

[12] HAN W H, GAO Y, LAW C K. Flame acceleration and deflagration-to-detonation transition in micro and macro channels: an integrated mechanistic study [J]. Combustion and Flame, 2017(176): 285-298.

[13] WANG C, DING J X, TAN S R, et al. High order numerical simulation of detonation wave propagation through complex obstacles with the inverse Lax-Wendroff treatment [J]. Communications in Computational Physics, 2015, 18(5):1264-1281.

[14] TAN S, WANG C, SHU C W, et al. Efficient implementation of high order inverse Lax-Wendroff boundary treatment for conservation laws [J]. Journal of Computational Physics, 2012(231):2510-2517.

[15] LI T, SHU C W, ZHANG M. Stability analysis of the inverse Lax-Wendroff boundary treatment for high order central difference schemes for diffusion equations [J]. Journal of Scientific Computing, 2017(70):576-607.