内置变径障碍物组管道内爆轰波传播规律的研究*

丁建旭，杜群贵，王新华，吴雨蒙

(1.广州特种机电设备检测研究院 国家防爆设备质量监督检验中心，广东 广州 510760；2. 华南理工大学 机械与汽车工程学院，广东 广州 510006；3. 广东工业大学 机电工程学院，广东 广州510006)

0 引言

典型危化品的爆炸灾害经常涉及爆轰波与复杂障碍物的强耦合作用，造成如反射和绕射等强非线性现象的发生，形成梯度较大的复杂爆炸环境，对人员及建筑设备造成严重的伤害和破坏。因此，研究爆轰波与管道内复杂障碍物的耦合作用机制，对工业重大爆炸灾害事故的预防及快速响应能力的提升均有着重要的现实意义。

随着描述爆轰的数学物理模型的不断完善，采用数值模拟方法研究爆轰波传播动力学行为已成为新的发展方向，受到了许多学者的关注[1-3]。张宝亮等[4]通过变化障碍物的高度及间距，数值研究了矩形障碍物对甲烷/氧气混合物的气相爆轰波传播过程的影响；随着数值技术的不断发展，Wang等[5]采用高精度WENO格式研究了一组矩形障碍物对爆轰波传播规律的影响；真实的爆炸环境往往是复杂的，文献[6-7]在无约束空间内研究了一簇方形和圆柱形障碍物对火焰加速以及爆燃转爆轰过程的影响机制；文献[8]采用商业计算软件对爆炸冲击波在建筑群中的传播规律进行了模拟，表明数值模拟技术在解决实际问题时具有一定的优越性。

从以上的讨论可以看出，许多学者为了简化模型和研究的方便，基本上采用了规则的如矩形等障碍物作为火焰传播的扰动对象。但是，在实际工程应用中所涉及的障碍物往往更加复杂。管道中变尺寸圆柱障碍物组对爆轰波传播过程的影响机制，目前相关的研究鲜有报道，为此，本文采用数值手段开展爆轰波在内置一组变径圆柱障碍物管道中的传播规律研究，以揭示预混气体爆轰波在复杂障碍物组作用下发生的复杂流场改变及爆炸特征参数的时空分布规律。

1 控制方程组及数值方法

采用两步化学反应[9-10]描述爆轰波与一组变径圆柱障碍物的相互作用过程，相应无量纲守恒型控制方程组的具体形式为：

(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为气体常数；r为密度；T为温度；P为压强；U和v分别表示x和y方向的速度分量；E1和E2为活化能；k1和k2为反应速率常数；r是比热率；a和b为控制反应进行的变量，当0

控制方程组的通量采用五阶精度的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.