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

(6)
式中:CFL取0.6;dx和dy為網格尺寸;c為聲速。
爆轟波在內置復雜障礙物管道中的傳播過程,特別在邊界附近的數值流場強烈地依賴于邊界的數值處理。因此,文中固壁邊界均采用高精度ILW方法[13-15]進行處理。
2 計算工況及數值結果討論
2.1 工況設置
在實際化工生產中,圓形輸運管道的管徑及其布局方式往往是復雜的。若發生爆炸事故,不同直徑管道的布局方式對爆炸場的影響規律值得深入探討,可為有限約束空間內輸運管道的布局優化以及對局部結構強化抗爆安全提供指導。矩形無量綱管道的長度為90,寬為15,內置5個可變半徑的圓柱形障礙物組。在計算中障礙物的個數保持不變,均為5個。工況1的圓柱形障礙物半徑均為3,工況2中的障礙物半徑從1逐漸增加到5,工況3的圓柱形障礙物半徑的變化規律與工況2正好相反,由大到小變化。在各個圓柱形障礙物的5個位置設置監測點Mijk,記錄壓力的變化,其中下標i代表工況,j代表從左到右的第幾個圓柱形障礙物,k代表每個圓柱形障礙物的5個監測點,從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邊界處理方法進行處理,保證固壁邊界附近與管道內部爆轟流場均為五階精度,提高數值分辨率。無量綱的初始條件為(ρ,u,v,p,α,β)=(1,0,0,0.16,1,1). 網格大小為dx=dy=0.017,保證了數值方法的收斂性[13]。
2.2 數值結果及討論
為了節約篇幅,文中僅以圓柱障礙物半徑逐漸增加的工況2和圓柱障礙物半徑逐漸減小的工況3的密度梯度云圖為例,考察爆炸場內復雜波系的演化過程,工況1條件下的爆轟波與5個等半徑圓柱障礙物作用的密度梯度云圖不再贅述。
圖1顯示5個障礙物半徑從1漸增到5爆轟波傳播過程中的密度梯度變化。從圖1(a)可以看出,由于障礙物1的半徑較小,對平面爆轟波陣面的破壞區域有限,形成的反射波經過長時間的演變,抵達管道上下壁面時強度損失較多。經上下壁面碰撞形成的反射波強度也較弱,繞射波陣面從三波點到對稱軸兩側的密度梯度線開口逐漸增大,如圖1(b)所示。從圖1(c)可以看出,在波陣面抵達圓柱形障礙物2之前,波陣面上存在的兩對相向運動的三波點即將發生碰撞。如圖1(d)所示,向下運動的三波點遭到破壞,而向上運動的三波點增強,且在圓柱形障礙物2左側發生了反射現象。圖1(e)表明,向上運動的三波點經管道壁面反射后強度增強,向著對稱軸方向運動,此時繞射波發生碰撞,同時,繞射陣面成弧形彎曲,但沒有明顯的局部爆炸產生。隨著反射波的持續相互作用以及壓縮波對未反應氣體的壓縮活化,在對稱軸附近的流場內部發生了局部爆炸,釋放大量的化學能,提升了局部的壓力和溫度,如圖1(f)所示。從圖1(g)中可以看出,局部爆炸形成的射流頭部加速傳播,并與障礙物3左側首先發生碰撞并反射,形成反射波。射流頭部橫向運動的橫波沿著波陣面繞障礙物固壁面運動,增強了圓柱壁面附近波陣面的傳播速度,但是,隨著傳播的繼續,三波點不斷向管道壁面移動,導致稀疏效應逐漸顯現,三波點以下陣面結構呈弧形彎曲狀,如圖1(h)所示。隨后,出現了與圖1(f)類似的爆轟波陣面結構,如圖1(i)所示。隨著爆轟波傳播路徑上障礙半徑的逐漸增大,繞射爆轟波陣面呈弧形相交結構,在相交軸線附近不再出現內凹結構,如圖1(k)和圖1(l),密度梯度強度分布出現變化,主要集中在局部爆炸射流的頭部,如圖1(l)所示,此時波陣面沒有出現明顯的三波點結構。

圖1 工況2條件下爆轟波傳播過程的密度梯度Fig.1 Density gradient contours in the propagation of detonation wave for case 2
圖2表示在五個障礙物半徑從5逐漸遞減到1的條件下爆轟波傳播過程中的密度梯度變化云圖。與工況2類似,爆轟波與圓柱障礙物組相互作用產生的連續反射和繞射是爆轟波動力學演化的主要機制。但是由于圓柱障礙物半徑的逐漸減小,從圖2(l)中可以看出,在爆轟波與障礙物5作用之前,波陣面上已經出現了多個三波點結構。從圖2(m)可以看出,障礙物5右側靠近軸線附近的三波點個數增多,其他位置三波點能量增強。由于爆轟波在水平傳播方向上自由通道漸寬,反射和繞射強度減小,隨著爆轟波的繼續傳播,波陣面的三波點大量出現,如圖2(n)所示,釋放的能量維持爆轟波繼續向下游傳播。

圖2 工況3條件下爆轟波傳播過程的密度梯度Fig.2 Density gradient contours in the propagation of detonation wave for case 3
通過以上分析可以發現,不同的變徑障礙物組形成了不同工況條件下阻礙物阻塞比漸變趨勢的變化,造成了爆轟波在水平傳播方向上約束程度的差異,以及障礙物對爆轟波反射和繞射強度的改變。這些耦合效應使得較大阻塞比管道段圓柱障礙物以及管道壁面反射和繞射碰撞產生的局部爆炸是爆轟波維持傳播的重要原因,而在較小阻塞比管道段,波陣面演化產生的多個三波點提供的大量的化學能是此段管道中爆轟波自持傳播的主要機制。
圖3(a)~(c)分別表示在工況1~3條件下5個障礙物的5個監測點的最大壓力。圖中p*是無量綱參考壓力[13]。在工況1條件下,爆轟波傳播通道中存在5個等半徑圓柱障礙物。從圖3(a)中可以看出,在爆轟波反射區監測點M111和M112的最大壓力均為2.1,由于繞射區的稀疏效應,使得監測點M114和M115的最大壓力均有所降低,約為1.3,在反射結束與繞射開始的臨界區的監測點M113的最大壓力為1.6,處于兩者之間,表明反射對爆轟波的最大壓力提升具有正反饋效應,而繞射會使最大壓力顯著地降低。隨著爆轟波的繼續傳播,與更多的圓柱形障礙物繼續發生相互作用,使得爆轟流場內部波系更加復雜,同一個圓柱障礙物壁面處的監測點的最大壓力不再具有類似的分布規律。但是各個圓柱形障礙物的正反射區域的監測點的最大壓力顯著高于其他監測點,特別是監測點M141的最大壓力躍升至10.0,為所有監測點的時空最大壓力,且M111的最大壓力值遠低于其他該工況下正反射處監測點所獲得的最大壓力。在繞射區的45°方向上監測點M114,M124,M134,M144和M154的最大壓力基本沒有顯著升高或降低。對于爆轟波傳播通道中圓柱形障礙物半徑漸增的工況2,如圖3(b)所示,在第一個半徑為1的障礙物壁面上獲得不同作用區域的最大壓力出現明顯的變化,正反射的監測點M211的最大壓力為8.0,遠大于該圓柱障礙物壁面其他位置監測點獲得的最大壓力。在工況2條件下,所有障礙物正反射區域的監測點的最大壓力遠大于對應圓柱障礙物其他監測點的最大壓力,監測點M241處的最大壓力為時空最大壓力,達到11.4,但是最后一個圓柱障礙物處監測點M251的最大壓力僅為6.2,小于其他任何障礙物正反射區監測點的最大壓力。在各個圓柱形障礙物的繞射區的45°方向上監測點的最大壓力雖然有降低的趨勢,但不是非常顯著,位于其他相同監測點位置但屬不同圓柱障礙物的最大壓力均有一定程度的改變。對于圓柱形障礙物半徑逐漸減小的工況3,如圖3(c)所示,第一個圓柱形障礙物的5個監測點的最大壓力歷史位于監測點M312,達到5.7,正反射區監測點M311的最大壓力甚至小于監測點M315,僅為2.9。在各個圓柱形障礙物的繞射區45°方向上監測點的最大壓力變化不顯著,但是對于各個圓柱形障礙物對應的在135°方向上的監測點出現了較為明顯的差異。

圖3 圓柱形障礙物監測點的最大壓力Fig.3 The maximum pressures at the monitoring points around cylindrical obstacles
從以上的分析可以看出,變徑圓柱障礙物組對爆轟波傳播具有重要的影響,特別在障礙物壁面處,總體上正反射區獲得的最大壓力均較大(除監測點M311),其他區域特別是45°方向上監測點的最大壓力基本沒有出現顯著的起伏變化,維持在較低值。就監測點處最大壓力的空間分布而言,3種工況的最大壓力的最大值均出現在第4個圓柱障礙物的正反射區的監測點處,且工況2在該點處獲得的最大壓力均大于障礙物組半徑維持不變或者逐漸減小的2種工況。
圖4為3種工況條件下中心對稱軸線上最大壓力隨時間的變化曲線,圖中t*是無量綱參考時間[13]。爆轟波繞過障礙物形成的稀疏波從管道上下兩個方向朝管道對稱中心匯聚并發生碰撞,在局部區域釋放大量的化學能,導致溫度和壓力的提升,同時,反射波在圓柱障礙物之間以及與管道上下壁面的不斷反射增加了軸線附近的流場復雜度,從而反映于軸線上最大壓力的周期性變化。對比三條曲線可以發現,對于相同障礙物半徑的工況1,最大壓力的第3個峰值為16.7,均大于其他工況的最大壓力的最大峰值,工況3的最大壓力的最大峰值與其他工況最大壓力的最大峰值相比最小。工況2的最大壓力的最大峰值比工況1和工況3均提前出現。

圖4 對稱軸線上最大壓力歷史Fig.4 The maximum pressure history along the central axis of tube
以上分析表明,變半徑障礙物組的不同半徑變化設置對中心軸線上的最大壓力的最大峰值具有重要的影響,等半徑障礙物組的工況1條件下獲得的最大壓力的最大峰值比障礙物半徑逐漸增加的工況2或遞減的工況3的最大壓力的最大峰值都大。
3 結論
1)較大阻塞比管道段圓柱障礙物以及管道壁面反射和繞射碰撞產生的局部爆炸是爆轟波維持傳播的重要原因,而在較小阻塞比管道段,波陣面演化產生的多個三波點提供的大量的化學能是此段管道中爆轟波自持傳播的主要機制。
2)在障礙物壁面處,正反射區的最大壓力較大,45°方向上的最大壓力維持在較低值。3種工況的最大壓力的最大峰值均出現在第4個圓柱障礙物的正反射區,且逐漸增大障礙物半徑的工況2的監測點最大壓力的最大峰值為所有工況最大。
3)不同障礙物半徑變化規律對中心軸線上的最大壓力的最大峰值具有重要影響,等半徑障礙物工況下獲得的最大壓力的最大峰值比障礙物半徑逐漸遞減或增加的最大壓力的最大峰值都大。
[1] 王公忠, 張建華, 李登科, 等. 障礙物對預混火焰特性影響的大渦數值模擬[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] 徐景德, 楊庚宇. 置障條件下的礦井瓦斯爆炸傳播過程的數值模擬研究[J]. 煤炭學報, 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]. 爆炸與沖擊, 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]. 中國安全科學學報, 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] 張寶亮, 丁玨, 王慶濤, 等. 約束空間可燃氣體燃燒爆轟特性的數值研究[J]. 中國安全生產科學技術, 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] 王剛, 張德良, 劉凱欣. 氣相爆轟波數值模擬中化學反應模型研究[J]. 高壓物理學報, 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] 王剛, 張德良, 劉凱欣. 對數值模擬氣相爆轟中二階段化學反應模型再研究[J]. 空氣動力學報,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.