Y 型通风采煤工作面瓦斯爆炸传播规律模拟研究*

2023-09-15刘佳佳聂子硕

爆炸与冲击 2023年8期

刘佳佳，张扬，张翔，聂子硕

（1.河南理工大学安全科学与工程学院，河南焦作 454003；2.安徽理工大学深部煤矿采动响应与灾害防控国家重点实验室，安徽淮南 232001；3.河南理工大学煤炭安全生产与清洁高效利用省部共建协同创新中心，河南焦作 454003；4.河南理工大学瓦斯地质与瓦斯治理国家重点实验室培育基地，河南焦作 454003）

近年来，我国煤矿生产状况得到显著改善，但煤矿重大事故依然时有发生。其中，瓦斯爆炸事故是我国煤矿最严重的事故之一，尤其是瓦斯爆炸所产生的动力灾害主要是高压冲击波和高温气体，都会对井下工作者、井下通风设施甚至矿井风硐造成毁灭性的破坏，造成非常恶劣的影响[1]。

学者们针对不同连接形式的瓦斯巷道开展了瓦斯爆炸研究，发现瓦斯巷道结构不同时，瓦斯爆炸传播规律也各有不同。Emami 等[2]、Blanchard 等[3]和孟显华等[4]证明了L 型巷道的弯曲结构能够提高火焰速度和超压。Yuan 等[5]对衍射爆轰的再起爆现象进行了研究，发现其实质是弯曲结构内外壁相互作用引起的爆轰能量的再分配。在此基础上，Cao 等[6]、Zhang 等[7]和Qiu 等[8]发现不同角度的弯曲结构对瓦斯爆炸超压传播规律的影响具有显著差异；Liu 等[9]发现，W 型巷道内，瓦斯爆炸超压衰减速率是一个由低到高再到低的过程；马恒等[10]证实了H 型巷道内爆炸冲击波会改变联络巷内气流的流动状态，爆炸后气流由双涡旋模式向单涡旋模式变化；高建良等[11]和高智慧等[12]研究了角联网络结构中的瓦斯爆炸传播特性，证明了角联网络对冲击波超压的衰减作用大于并联网络对超压的衰减作用；景国勋等[13]在T 型巷道末端利用不同约束条件进行了瓦斯爆炸试验，发现瓦斯最大爆炸压力和管道中火焰锋面传播速度存在明显差异；丁浩等[14]、Yan 等[15]、周辉等[16]和时本军等[17]均利用管道附带的特殊腔体结构进行了瓦斯爆炸研究，发现空腔结构具有消焰和消波性能，可以作为一种新型瓦斯抑爆装置。

除了研究巷道结构本身对瓦斯爆炸传播规律的影响，学者们同样对爆炸气体的组分和含量进行了大量的研究。Wang 等[18]研究了不同氢气添加比的氢气/甲烷/空气混合物的爆炸传播规律，发现爆炸超压随氢气添加比的增加而增大；Li 等[19]开展了12 种不同体积分数的甲烷/空气混合物的预混爆炸试验，发现当气体体积分数为9.5%时，超压峰值比其他体积分数至少高15%；Gao 等[20]发现瓦斯浓度的增大使超压峰值呈线性变化，且瓦斯浓度和瓦斯体积对超压峰值和超压平均上升速率均有显著影响；罗振敏等[21]、梁国栋[22]和贾泉升等[23]均利用20 L 爆炸特性测试系统研究了不同体积分数瓦斯的爆炸压力传播规律，证明爆炸后产生的H2和CO 会加速瓦斯的爆炸，且与初始瓦斯体积分数呈一次函数关系。

综上所述，Y 型通风对瓦斯涌出量较大的工作面解决回风流中瓦斯浓度过高和上隅角瓦斯积聚有明显效果，已被越来越多的矿井工作面采用[24]。但是，目前尚未见关于Y 型通风采煤工作面的瓦斯爆炸研究，尤其缺少针对该方面的数值模拟研究。本文中，利用Fluent 软件结合大涡模拟（large eddy simulation）数学模型，研究Y 型通风采煤工作面瓦斯爆炸冲击波传播规律，以期为采煤工作面瓦斯爆炸灾害防治提供理论基础和技术支撑，预防工作面瓦斯爆炸事故发生。

1 模型建立

1.1 物理模型

1.1.1 基本假设

瓦斯爆炸是一个复杂快速的化学反应过程，为了便于模拟计算，使模拟结果更加可靠，对模拟过程作如下假设[25]：

(1) 预混气体和燃烧产物为满足状态方程的理想气体；(2) 气体混合物的比热容遵循混合规则，各组分的比热容为温度的函数；(3) 物理模型壁面为刚性绝热壁面，不会产生相对位移；(4) 瓦斯爆炸反应为单向不可逆反应；(5) 物理模型中瓦斯充填区为常态的瓦斯和空气均匀混合气体，且为理想状态。

1.1.2 物理模型及网格划分

利用ANSYS Workbench 软件，结合余吾煤矿N2105 工作面实际情况，建立如图1 所示的Y 型通风工作面三维物理模型。其中胶带顺槽和进风顺槽均为进风巷道，以回风巷道和工作面交界中心处为原点，工作面走向为Y轴，回风巷道为X轴正方向，垂直于工作面方向为Z轴，进回风巷道截面尺寸为5.0 m×3.5 m，工作面倾斜长度为285.0 m，采煤工作面宽5.0 m，采高3.5 m。利用ANSYS 软件自带的Mesh 功能对模型进行结构化网格划分。

图1 余吾煤矿N2105 工作面的三维物理模型Fig.1 A 3D physical model of N2105 working face in Yuwu Coal Mine

1.2 数学模型

大涡模拟采用非稳态的N-S 方程，对大尺度涡的湍流脉动进行计算，对小尺度的脉动进行过滤，通过亚网格模型对小尺度脉动产生的影响进行模型假设，经过过滤后的大涡模拟的控制方程[26-28]如下。连续性方程为：

1.3 边界条件与初始条件

模拟环境设为常温常压，设定出口的组分为氧气21.2% ，氮气78.8% ，压力出口的湍流强度均设为5%，壁面的粗糙常数为0.5，不考虑壁面的粗糙度厚度。

采用体积分数为9.5%的CH4-空气混合气体。瓦斯填充区尺寸为6.0 m ×5.0 m×3.5 m，瓦斯填充区域内（除点火区外）的初始条件为：初始温度T0=300 K，气体中各组分的质量分数分别为w(CH4)=0.055，w(O2)=0.220，w(H2O)=0，w(CO2)=0，w(N2)=0.725。一般空气区的初始操作条件为：初始温度T0=300 K，气体中各组分的质量分数分别为w(CH4)=0，w(O2)=0.22，w(H2O)=0，w(CO2)=0，w(N2)=0.78。采用Patch 高温区实现点火，假设氧化反应完全，初始条件为：w(CH4)=0，w(O2)=0.220 00，w(H2O)=0.119 25，w(CO2)=0.145 60，w(N2)=0.515 15。

1.4 数学模型验证

开展模拟研究之前，需要对数学模型进行验证，以确保数学模型的可靠性。根据洪溢都等[29]开展的管道爆炸试验，建立如图2所示的三维物理模型并进行网格划分，管道模型长为5.0 m，横截面尺寸为8 cm×8 cm。自管道左端2.0 m 处开始设置压力监测点，每隔0.5 m 设置一个监测点，并采用Patch 高温区方式实现点火。

图2 模拟管道示意图Fig.2 Schematic diagram of simulated pipeline

表1 给出了爆炸超压的数值模拟结果与实验数据的对比。由表1 可知，数值模拟结果与实验结果之间的最大相对误差为11.3%，最小相对误差仅为1.7%，除3.5 m 处测点的相对误差较大外，其他测点的相对误差均在10%以内。可知模拟结果与实验结果吻合较好，验证了数学模型的可靠性。

表1 爆炸超压模拟结果与实验结果的对比Table 1 Comparison between simulation and experimental results of explosion overpressure

2 瓦斯爆炸模拟关键参数优化

2.1 网格尺寸参数优化

为了使模拟结果更加可靠，需要对模拟关键参数进行优化，如：网格尺寸、迭代步长和点火温度等。参数确定过程共选取4 个监测点，监测点位置分布如图3 所示（后续参数优化监测点布置与其相同），监测点A～D的坐标分别为(12.5，0，1.75)、(0，-7.5，1.75)、(27.5，0，1.75)和(0，-27.5，1.75)。

图3 参数优化监测点示意图Fig.3 Parameter optimization monitoring points

利用ANSYS Workbench 中的Mesh 模块将整体全尺寸模型分别划分为0.8、0.6、0.5、0.4 和0.2 m 的不同尺寸网格，如表2 所示。

表2 网格分布Table 2 Grid distribution

通过对比分析不同网格尺寸情况下各测点超压及超压峰值到达时间的变化规律，确定合理的网格尺寸。各测点超压峰值和超压峰值到达时间的变化规律分别如图4～5 所示。

图4 超压峰值随网格尺寸的变化Fig.4 Variations of overpressure peak with mesh size

由图4 可知，随着网格尺寸从0.8 m 减小至0.4 m，测点A、B和D处的超压峰值均不断增大；网格尺寸由0.5 m 减小至0.4 m 时，测点C处的超压峰值有所下降，但幅度并不明显；当网格尺寸从0.4 m 继续减小至0.2 m 时，测点A～D处的超压峰值基本保持不变。由图5 可知，随着网格尺寸的不断减小，测点A～D处的爆炸超压峰值所到达的时间随网格尺寸的减小而缩短；当网格尺寸从0.4 m 减小至0.2 m 时，爆炸超压峰值所到达的时间变化幅度较小，网格尺寸的影响基本可以忽略不计。综合考虑数值模拟可靠性和计算机配置的影响可知，在不改变模拟结果的基础上，太小的网格尺寸会大幅降低模拟计算速度。因此，最终确定最合理的网格尺寸为0.4 m。

图5 超压峰值到达时间随网格尺寸的变化Fig.5 Variations of overpressure peak arrival time with mesh size

2.2 迭代步长参数优化

在确定网格尺寸为0.4 m 的基础上，分别选取1.00、0.50、0.25、0.10 和0.05 ms 等5 种不同的迭代步长，通过对比分析不同迭代步长情况下各测点超压及超压峰值到达时间的变化规律，确定合理的迭代步长。各测点超压峰值和超压峰值到达时间的变化规律分别如图6～7 所示。

图6 超压峰值随迭代步长的变化Fig.6 Variation of overpressure peak with time step

由图6 可知，迭代步长从1.00 ms 缩短至0.10 ms 时，测点A～D处的超压峰值均不断变化，当迭代步长从0.10 ms 继续缩短至0.05 ms 时，测点A、B和D处的超压峰值基本保持不变，测点C处的超压峰值有所浮动，但相对误差在5%以内；由图7 可知，随着迭代步长的不断缩短，测点A～D处的爆炸超压峰值到达的时间随迭代步长的缩短而提前，当迭代步长从0.10 ms 缩短至0.05 ms 时，爆炸超压峰值到达的时间基本保持不变。综合考虑数值模拟可靠性和计算机配置的影响，最终确定合理的迭代步长为0.10 ms。

图7 超压峰值到达时间随迭代步长的变化Fig.7 Variation of overpressure peak arrival time with time step

2.3 点火温度参数优化

在确定网格尺寸为0.4 m，迭代步长为0.10 ms 的基础上，分别选取1 500、1 600、1 800、2 000 和2 500 K 等5 种不同的点火温度，通过对比分析不同点火温度情况下各测点超压及超压峰值到达时间的变化规律，确定合理的点火温度。各测点超压峰值和超压峰值到达的时间变化规律分别如图8～9 所示。

图8 超压峰值随点火温度的变化Fig.8 Variation of overpressure peak with ignition temperature

由图8 可知，点火温度从2 500 K 降至1 800 K 时，测点A处超压峰值不断降低，测点B～D处超压峰值有所变化，但变化幅度不大；当点火温度从1 800 K 继续降低至1 500 K 时，测点A～D处超压峰值变化幅度很小，甚至基本保持不变。这是因为测点A距离爆炸源过近，温度越高，活化能也越高，爆炸威力相比其他测点更大。由图9 可知，随着点火温度不断降低，测点A～D处爆炸超压峰值到达的时间基本保持不变。综合考虑数值模拟可靠性和计算机配置的影响，最终确定合理的点火温度为1 800 K。

图9 超压峰值到达时间随点火温度的变化Fig.9 Variation of overpressure peak arrival time with ignition temperature

3 结果分析与讨论

3.1 瓦斯爆炸超压分析

3.1.1 巷道瓦斯爆炸超压

模拟过程中共选取32 个监测点，监测点位置分布为：以回风巷道、胶带顺槽和工作面三者连接处的中心点为原点，回风巷道为X轴方向，工作面为Y轴方向，爆炸源中心点距离原点2.5 m 处。原点设置为监测点1，在回风巷道内，水平距离原点17.5 m 处设置监测点2，并每隔10.0 m 分别设置监测点3～5；在胶带顺槽内，水平距离原点7.5 m 处设置监测点6，并每隔10.0 m 分别设置监测点7～10；在工作面内，垂直距离原点2.5、7.5、17.5 和27.5 m 处分别设置监测点11～14，之后每隔20.0 m 分别设置监测点15～26，在与监测点26 距离17.5 m 处设置监测点27；在进风顺槽内，水平距离原点7.5 m 处设置监测点28，之后每隔10.0 m 分别设置监测点29～32。

图10 为Y 型通风采煤工作面的瓦斯爆炸超压（Δp）变化云图。图11～14 为Y 型通风采煤工作面瓦斯爆炸超压和超压到达时间。

图10 瓦斯爆炸超压传播云图Fig.10 Contours of gas explosion overpressure propagation

图11 回风巷道瓦斯爆炸传播规律Fig.11 Gas explosion propagation law in return airway

由图10 和进风顺槽、胶带顺槽、回风巷道和工作面瓦斯爆炸超压的模拟结果可知，瓦斯爆炸后，回风巷道内的瓦斯爆炸超压峰值明显高于胶带顺槽和工作面，距离爆源15.0 m 处，回风巷道的瓦斯爆炸超压峰值为451.874 kPa，胶带顺槽和工作面的瓦斯爆炸超压峰值分别为228.776 和227.886 kPa，主要原因是在瓦斯爆炸区内存在分叉巷道，由于分叉巷道具有分压的作用，因此胶带顺槽和工作面内的瓦斯爆炸超压峰值明显低于回风巷道。

由图11 可知，回风巷道内超压峰值与爆源之间的距离（d）符合指数函数关系Δp=142.761+1061.871 e-d/12.187，决定系数R2=0.990，超压峰值到达时间（Δt）与爆源之间的距离符合线性函数关系Δt=1.699d-6.995，决定系数R2=0.976，拟合度均良好。

回风巷道内，瓦斯爆炸超压峰值从水平距离爆炸源15.0 m 处的451.874 kPa 降至距离爆炸源45.0 m处的174.506 kPa，共需51.30 ms。在瓦斯爆炸传播初期，瓦斯爆炸超压传播速率较快，且超压衰减率随时间变化逐渐降低。究其原因，是瓦斯爆炸初期，火焰燃烧剧烈，大量瓦斯被引爆，反应剧烈，气体呈紊流状态，加速了冲击波向出口处传播，冲击波经爆源拐角处壁面反射，反射波与前驱冲击波叠加，使得压力下降相对缓慢。

由图12 可知，胶带顺槽内超压峰值与爆源之间的距离符合指数函数关系Δp=115.212+284.816 e-d/19.335，决定系数R2=0.990，超压峰值到达时间与爆源之间的距离符合线性函数关系Δt=1.709d+6.443，决定系数R2=0.997，拟合度均良好。

图12 胶带顺槽瓦斯爆炸传播规律Fig.12 Gas explosion propagation law along belt fluting

胶带顺槽内，瓦斯爆炸超压峰值从水平距离爆炸源2.5 m 处的370.103 kPa 下降至50.0 m 处的132.264 kPa，共需82.00 ms，胶带顺槽内超压峰值衰减速率整体上小于回风巷道，这是因为瓦斯主要积聚在回风巷道，当瓦斯被引爆后，回风巷道附近反应更加剧烈，产生的能量更高，冲击波传播速度更快，由于冲击波经过拐角处会发生反射，能量损耗更大，导致胶带顺槽和工作面的冲击波传播速度低于回风巷道，超压衰减率相对较小，超压峰值时间间隔更长。

由图13 可知，工作面内超压峰值与爆源之间的距离符合指数函数关系Δp=61.369+249.108 e-d/33.43，决定系数R2=0.933，超压峰值到达时间与爆源之间的距离符合线性函数关系Δt=2.241d-7.959，决定系数R2=0.976，拟合程度均良好。

图13 工作面瓦斯爆炸传播规律Fig.13 Gas explosion propagation law of working face

工作面内，瓦斯爆炸超压峰值从垂直距离爆源0 m 处的370.103 kPa 降至285.0 m 处的73.104 kPa，共需644.40 ms。随着冲击波的传播，前50.0 m 工作面超压峰值衰减率较大，甚至在距离爆源17.5 m 处超压峰值异常增大（异常点1）；在工作面距离爆源285.0 m 处，同样出现了瓦斯爆炸超压大幅升高的现象（异常点2），而超压峰值时间间隔却逐渐延长。究其原因，是在拐角处，瓦斯爆炸超压扰动空气流动形成湍流，冲击波经壁面反射叠加，反射波的能量高于冲击波传播损耗的能量，造成超压峰值不降反升。因此，矿井瓦斯爆炸防治需要特别注意具有明显角度变化的巷道。

由图14 可知，进风顺槽内超压峰值与爆源之间的距离符合指数函数关系Δp=42.697+50.262 e-d/4.955，决定系数R2=0.987，超压峰值到达时间与爆源之间的距离符合线性函数关系Δt=2.230d+650.807，决定系数R2=0.999，拟合度均良好。

图14 进风顺槽瓦斯爆炸传播规律Fig.14 Gas explosion propagation law in the inlet channel

因为拐角区域的冲击波反射叠加作用，传播至进风顺槽的超压峰值明显增大，从水平距离爆源2.5 m 处的73.104 kPa 降至50.0 m 处的47.495 kPa，共需105.00 ms。瓦斯爆炸超压峰值在761.50 ms 时异常增大，但超压峰值时间间隔由25.00 ms 缩短为20.00 ms。究其原因，是冲击波由工作面向进风顺槽传播经过90 度拐角时，冲击波发生反射，随后在50.0 m 处反射波追上前驱冲击波，造成超压峰值增大（异常点3）。

3.1.2 巷道分叉对瓦斯爆炸超压影响分析

为了更好地研究巷道分叉对瓦斯爆炸冲击波的影响，并分析工作面和进风顺槽内异常点1～3 的形成原因，图15～16 给出了胶带顺槽和工作面距离原点50.0 m 范围内瓦斯爆炸超压峰值和超压峰值时间的对比，图17 给出了35.00、145.00 和700.00 ms 时在不同拐角处的瓦斯爆炸超压矢量图。

图15 爆炸超压峰值对比Fig.15 Comparison of gas explosion overpressure peaks

由图15～16 可知，距巷道分叉口50.0 m 范围内，胶带顺槽内瓦斯爆炸超压峰值整体大于工作面，第1 个拐角中心处的超压峰值为370.103 kPa，当超压峰值由第1 个拐角点向胶带顺槽内传播时，距离巷道分叉口7.5 m 处，超压峰值降至273.899 kPa，超压衰减率为25.99%，向工作面传播时，距离巷道分叉口7.5 m 处，超压峰值降至218.233 kPa，超压衰减率为41.03%。因此，与爆源拐角相连的工作面内的超压衰减率明显大于与爆源直接连通的胶带顺槽，瓦斯爆炸时与分叉巷道相比，与爆源直接连通的巷道更具危险性。

瓦斯爆炸35.30 ms 时，距巷道分叉口17.5 m 处（异常点1），工作面瓦斯爆炸超压峰值上升至227.886 kPa，与距巷道分叉口7.5 m 处超压峰值相比，增幅为4.42%，胶带顺槽内达到超压峰值220.372 kPa 所需时间为40.40 ms。究其原因，是在35.00 ms 时，工作面在该处附近形成湍流（见图17(a)），湍流的形成造成空气流动紊乱，加速了冲击波的传播。这也解释了图16 中相同距离处工作面的超压峰值到达时间比胶带顺槽短的原因；由于冲击波在拐角处经过反射重叠，反而使得工作面瓦斯爆炸超压峰值有所上升，这也是图15 中出现拐点的原因。

图16 爆炸超压峰值到达时间对比Fig.16 Comparison of overpressure peak arrival time of gas explosion

图17 瓦斯爆炸超压矢量图Fig.17 Vector diagrams of gas explosion overpressure

分析图17 可知，在瓦斯爆炸传播初期，爆炸冲击波锋面最先经过工作面右侧，如图10(a)所示，由于存在拐角，所以爆炸冲击波在工作面右侧扰动空气流动，此时该处最先出现湍流区，此后，爆炸冲击波锋面逐渐向工作面和胶带顺槽传播，而工作面左侧和胶带顺槽之间也存在拐角，因此湍流区逐渐由工作面右侧向左侧移动。这说明随着瓦斯爆炸冲击波传播，湍流区并不是一成不变的，但是湍流区主要在分岔口附近形成；瓦斯爆炸700.00 ms 时（见图13），瓦斯爆炸冲击波由工作面传播至进风顺槽时，超压峰值有所上升（异常点2），这是因为冲击波经壁面反射叠加，在第2 个拐角处产生湍流区，使该处超压峰值由50.597 kPa 上升至73.104 kPa，超压增长率为44.48%；在湍流扰动以及壁面反射的双重作用下，反射波与前驱冲击波重叠，反射波能量大于沿程损耗的能量，造成异常点3 的超压峰值增大。

3.2 瓦斯爆炸火焰传播特性分析

通过分析瓦斯爆炸温度变化云图对瓦斯爆炸火焰传播过程进行分析。图18 给出了Y 型通风采煤工作面瓦斯爆炸火焰温度随时间的变化云图。

图18 火焰的温度云图Fig.18 Contours of flame temperature variation

由图18 可知，在受限空间平直巷道内，瓦斯爆炸火焰传播主要经过如下4 个阶段。

第1 阶段，球形自由扩散阶段。瓦斯爆炸8.00 ms 时，火源开始向周围的可燃气体快速传播，火焰以球形形状燃烧并向周围环境自由扩散，已经燃烧的气体会释放出大量的热量，周围未燃烧的气体在高温作用下也会向周围环境自由扩散。

第2 阶段，蘑菇状加速传播阶段。瓦斯爆炸11.00 ms 时，火焰在壁面处回弹，壁面两端火焰传播相对中部较快，随后壁面限制火焰传播，火焰呈蘑菇状；50.00 ms 时，火焰碰撞壁面回弹形成湍流，火焰阵面被拉伸，两侧火焰传播加快，火焰阵面逐渐形成一个平面。

第3 阶段，郁金香形扭曲传播阶段。瓦斯爆炸85.00 ms 时，由于湍流作用，两侧火焰传播速度高于巷道中部，火焰前锋发生反转，火焰阵面向内凹陷，形成典型的郁金香形火焰。

第4 阶段，火焰退散阶段。瓦斯爆炸2 008.00 ms 时，回风顺槽内火焰基本熄灭，12 078.00 ms 时，胶带顺槽内火焰基本熄灭。

由上述讨论可知，无角度变化的平直巷道内，火焰从传播至最远距离处开始消散，而存在角度变化的工作面内火焰传播规律明显不同。瓦斯爆炸11.00 ms 时，火焰顺着工作面右端逐渐向工作面内蔓延，工作面左侧短暂形成一个安全区，并无火焰进入；11.00～85.00 ms，火焰逐渐进入工作面左侧，在工作面内整体向前延伸，但右侧火焰传播距离整体大于左侧火焰传播距离，火焰整体呈坡状，且工作面火焰传播至一定距离处开始呈消散趋势；约2 008.00 ms，工作面火焰沿着工作面继续向前传播；12 078.00 ms时，胶带顺槽和回风顺槽内火焰基本都已熄灭，但是工作面火焰继续向前传播，且随着时间推移，工作面火焰从火焰进口处开始熄灭。