大面积采空区煤自燃环境下气体运移规律研究

2022-09-16张鹏宇陈晓坤魏高明

煤矿安全 2022年9期

张鹏宇，陈晓坤，赵亮，魏高明，王奇

（1.西安科技大学安全科学与工程学院，陕西西安 710054；2.西安科技大学煤火灾害防治陕西省重点实验室，陕西西安 710054；3.兖矿能源集团公司鲍店煤矿，山东邹城 273500）

随着现代化矿井开采强度与机械化程度不断提升，矿井采空区面积和范围随之增加，形成大面积采空区[1]。在工作面正常推采期间，采掘活动以及顶板周期来压使得沿空侧预留煤柱压酥破碎，形成大量裂缝，导致采空区之间相互贯通，为煤的氧化升温提供漏风通道[2]。随着大面积采空区持续供氧，采空区遗煤长时间被氧化，增大了采空区煤自燃危险性。采空区煤自燃产生如CH4、CO2、CO、C2H6、C2H4、C3H8、H2、SO2、H2S 等有害气体，会在漏风通道内运移并涌出[3-5]。因此，采空区煤自燃环境作用下气体非控运移有可能引发瓦斯、煤尘爆炸和矿井通风系统紊乱等一系列严重的次生灾害[6-7]。采空区煤-气-热共生问题日益严峻，大面积采空区流场特性和有害气体运移机制亟待进一步深入研究。

国内外学者对采空区流场特征进行了大量研究，主要采用现场观测、相似模拟实验和数值模拟方法[8-10]。由于采空区条件复杂，人员无法进入其内部，采空区隐蔽灾害信息判断和气体流场探测极为困难[11]。物理相似模拟实验虽然在一定程度上解决了现场对于采空区垮落环境下束管监测系统监测点位破坏和数据局限性等问题，但在实验中不可避免产生的误差，无法在全域条件下研究采空区时空演化规律。而数值模拟可以弥补这一缺陷，并进一步为物理相似模拟实验的研究提供验证依据，为现场工程应用的可行性提供指导[12-13]。许多研究人员已经建立了多物理场耦合模型，旨在分析采空区中的气体渗流[14-15]、气体浓度场[16-17]、以及温度场[18-19]变化规律，这些研究为采空区模拟复杂流场和气体运移提供了必要的理论基础。由于已建立数值模型通常基于单一工作面采空区的漏风规律进行研究，缺乏考虑采空区连通后气体流场的运移规律；此外，忽略了由于采空区内部气体密度不同，煤自燃温度导致采空区压力以及气体渗流和扩散能力的改变。为此，针对大面积采空区建立了一个涉及多孔介质的气体运输、扩散以及能量传递之间复杂相互作用的多物理场耦合模型。考虑在重力影响下，煤自燃温度会使大面积采空区局部气体运动变得更加复杂，根据大面积采空区压力、温度和气体浓度分布规律，对采空区气体流场形成过程进行分析，为大面积采空区的漏风控制和煤自燃防治提供科学指导。

1 大面积采空区多物理场耦合模型

1.1 基本守恒方程

流体运动必须遵循3 个规则，即质量守恒、动量守恒和能量守恒[20]。

质量守恒方程又称为连续性方程，是指单位时间内微元体中流体质量的增加等于同一时间间隔内流入该微元体的净质量，如式（1）：

式中：ρ 为混合气体的密度，kg/m3；ε 为采空区孔隙度；u 为气体渗流速度，m/s；Sm为采空区遗煤耗氧源项，kg/（m3·s）；t 为时间，s，。

动量守恒方程描述了流体系统所受外力与流体系统动量变化之间的关系，即微元体中流体动量的增加率等于作用在微元体上各种力之和，如式（2）[21]：

式中：p 为气体静压，Pa；τ 为黏性应力张量；ρgn为重力体积力，kN/m3；α 为多孔介质渗透率，m2；μ 为气体动力黏度，Pa·s。

能量守恒方程建立在物质运动变化过程中控制体内能量转换的等量关系，决定了采空区高温点的温度分布和变化[22]。采空区气体在运移过程中，温度的变化同时也会引起采空区气体动能和内能改变，如式（3）：

式中：ρf、ρs为气体和固体介质的密度，kg/m3；Cpf、Cps为气体和固体介质的定压比热容，J/（kg·K）；kf、ks为气体和固体介质的导热系数，W/（m·K）；T 为温度，K；hi为组分i 的焓，kJ/mol；Ji为组分i 的扩散通量，kg/（m2·s）；Sf为遗煤的耗氧放热源项，W/m3。

煤氧化反应速率r 与氧气浓度和环境温度有关，可根据Arrhenius 方程定义为式（4）：

式中：A 为指前因子，s-1；E 反应活化能，kJ/mol；R 为普适气体常数，kJ/（mol·K）；CO2为氧气浓度，kmol/m3；n 为反应级数。

指前因子A 和反应活化能E 可根据煤低温氧化实验获得，利用Fluent 二次开发工具UDF 编写遗煤耗氧及放热源项。

1.2 粒子输送方程

粒子输送方程是指任一瞬时系统内物理量（质量、动量和能量）随时间的变化率等于该时间控制体内物理量的变化率与通过控制体表面净通量之和，是流体质点物理量在流场中在时间空间上输运传送变化过程[23]。粒子输送方程如式（5）：

式中：ci为气体组分i 的质量分数；Qi为组分i的增减源项，kg/（m3·s）。

1.3 大面积采空区孔隙率和渗透率方程

在采空区形成过程中，采空区上覆岩层垮落破碎，垮落岩石填充采空区形成多孔介质空间。在采空区边界煤岩层的支撑作用和地层的垮落压实作用下，上覆岩层的变形呈现出一定垮落规律，用碎胀系数Kp表征煤岩块堆积与压实状况，采空区走向上碎胀系数Kp（x）如式（6）：

采空区内部的孔隙度分布与岩石的破碎程度相关，当岩石破碎程度较小时，岩块尺寸较大，孔隙度也变大；在采空区中心区域，上覆岩层在垮塌过程中破裂较为剧烈，在上覆岩层的自重压实作用下，岩块尺寸变小，孔隙度也越小。因此，有研究表明采空区孔隙度分布与边界距离呈现“O”型分布形状[24]。

渗透率是在压差作用下，介质中允许流体通过的能力，是用于描述流体在介质内传导液体能力的参数。其中，渗透率越大，流体通过多孔介质的速度越快。采空区孔隙率ε 与碎胀系数Kp如式（7），采空区渗透率α 根据Blake-Kozeny 公式计算得到[25]：

式中：Dp为垮落带煤岩块平均直径，m。

2 大面积采空区数值模型

2.1 工程实例

鲍店煤矿位于山东省济宁市，隶属兖州煤业股份有限公司，鲍店煤矿煤样的基础参数为：水分1.910%，灰分12.790%，挥发分13.040%，固定碳52.760%，碳元素72.660%，氢元素4.274%，氮元素1.276%，临界温度70 °C。该地区煤样为气煤，易发生自燃。5316 工作面邻近多个采空区，煤层埋藏深、通风系统和地质条件复杂、采空区范围大、内部遗煤多、漏风通道多。5316 工作面推采过程中邻近采空区漏风流场极为复杂，存在煤自燃隐患。以U 型通风5316 工作面采空区和邻近五采区的5309、5310、5311 和5312 串联采空区为物理模型进行3D 数值模拟研究。大面积采空区模型如图1。

图1 大面积采空区模型Fig. 1 Large area goaf model

图1（a）是工作面的二维示意图，工作面回采率约90%，两侧留有5 m 煤柱。5309、5310、5311 和5312 采空区已封闭，5316 为生产面。CFD 模型如图1（b），大面积采空区采用六面体网格，为提高计算精度，对边界层网格进行加密处理。

2.2 模型验证

5309 和5310 工作面采空区预留钻孔，通过气相色谱监测采空区内部气体体积分数，反映采空区气体渗流情况，大面积采空区氧气体积分数的实测值与模拟值的对比验证如图2。

图2 大面积采空区氧气体积分数的实测值与模拟值的对比验证（风速v=1.5 m/s，高度切片Z=1 m）Fig. 2 Verification of the measured and simulated values of oxygen concentration in large area goaf（ventilation flux v=1.5 m/s, extracted slices Z=1 m）

数值模拟结果与现场实测O2体积分数变化规律基本一致，由于煤柱漏风，随着埋深的增加，采空区内O2体积分数迅速下降。风流从进风口到回风口的过程，5309 和5310 工作面O2体积分数也有所降低。对于5316 工作面采空区，当进风侧采空区距工作面72 m 时，从散热区到氧化升温区，O2体积分数降低到18%，距工作面150 m 时，从氧化区到窒息区，O2体积分数降低至8%，氧化带宽为78 m；回风侧距工作面82 m 时，O2体积分数降低至18%，从散热区到氧化升温区，距工作面109 m 时，O2体积分数降低8%，从氧化区到窒息区，氧化带宽度为27 m。

3 结果与讨论

3.1 煤自燃环境下大面积采空区温度场演化规律

大面积采空区温度随时间的演化云图如图3，大面积采空区温度随高度变化云图如图4，大面积采空区温度分布曲线如图5。

图3 大面积采空区温度随时间的演化云图（风速v=1.5 m/s，高度切片z=1 m）Fig.3 Cloud images of temperature evolution with time in large area goaf（ventilation flux v=1.5 m/s,extracted slices z=1 m）

图4 大面积采空区温度随高度变化云图（风速v=1.5 m/s，第30 d）Fig.4 Cloud images of temperature variation with height in large area goaf（ventilation flux v=1.5 m/s,30th day）

图5 大面积采空区温度分布曲线（风速v=1.5 m/s，第30 d，切片y=320 m）Fig.5 Temperature distribution curves of large area goaf（ventilation flux v=1.5 m/s, 30th day,extracted slices y=320 m）

由图3 可知：随着时间的推移，采空区高温区域向四周扩散，并随风流运动方向迁移，通过温度的等值线可以看出温度向外扩散和蔓延的幅度较小，与环境温度（300 K）形成较大的温度梯度；由于风流流经隅角，在隅角处积聚，风流对于热量耗散作用较小，高温点逐渐趋向隅角处，形成高温点向隅角偏移现象。

由图4 和图5 可知：受传导和对流作用影响，温度在采空区竖直方向上进行传递，高温区域面积没有明显改变；在切片位置为z=1、5、10 m 处的最高温度分别为346、339、336 K，温度上升梯度分别为1.75 K/m 和0.6 K/m；在离开遗煤区域后，距火源距离越远，温度上升梯度越小，这是由于采空区属于多孔介质空间，内部的岩石和空气的导热能力较差，导致热量不断积聚，从而形成高温点并能够持续的维持其影响范围。

3.2 煤自燃环境下大面积采空区气体运移规律

压力是促使大面积采空区内部气体运移的直接动力，掌握采空区氧气体积分数分布和流动规律是防治自然发火的关键技术基础。大面积采空区不同高度下静压分布云图如图6，大面积采空区静压分布曲线如图7，大面积采空区氧气浓度分布云图如图8，大面积采空区氧气浓度分布曲线如图9，不同风速条件下大面积采空区氧气云图如图10。

图6 大面积采空区不同高度下静压分布云图（风速v=1.5 m/s，第30 d）Fig.6 Cloud images of static pressure distribution at different heights in large area goaf（ventilation flux v=1.5 m/s, 30th day）

图8 大面积采空区氧气体积分数分布云图（风速v=1.5 m/s，第30 d）Fig.8 Cloud images of oxygen concentration distribution in large area goaf（ventilation flux v=1.5 m/s, 30th day）

图9 大面积采空区氧气体积分数分布曲线（风速v=1.5 m/s，第30 d）Fig.9 Oxygen concentration distribution curves of large area goaf（ventilation flux v=1.5 m/s, 30th day）

图10 不同风速条件下大面积采空区氧气云图（切片高度z=10 m）Fig.10 Oxygen cloud images of large area goaf under different wind speeds（extracted slices z=10 m）

由图6 和图7 可知：①无煤自燃的情况下，通风作用会在进风侧形成较大压力，在回风侧形成负压。靠近进风侧采空区深部压力随着高度增加而增大，而靠近回风侧采空区压力随高度增大而急剧减小，这是由于进风侧风流随漏风通道进入采空区，向采空区上部运移，导致采空区上部压力增大，最终风流汇聚在回风侧流出，致使压力骤降；②煤自燃产生高温会使气体膨胀，密度减小，从而在采空区底部形成负压区域，相应的受浮力作用影响，气体受热膨胀向上运动，在采空区顶部积聚形成高压区域。

由图8 和图9 可知：①氧气进入采空区，受孔隙率分布影响，在向采空区深部扩散过程中呈现出不同的浓度分布，且随着采空区竖直高度越高氧气体积分数越低，这是由于氧气密度大于空气，受重力影响，导致氧气在采空区底部积聚；②在煤自燃作用下，可以看到温度产生的火风压作用使局部风流方向发生改变，使得采空区煤自燃区域不断卷吸周边氧气，产生强对流作用涌入采空区上部区域并在隅角处积聚。

加强通风会导致采空区漏风强度增大，进而影响采空区气体的运移和积聚。

由图10 可知：随着进风风速的增加，进风侧的压力明显增大，采空区漏风更加严重；因此，进风风速越大，氧化升温带越宽，并且明显从工作面向采空区深部移动；随着进风风速从0.5 m/s 增加到2.0 m/s，5309 采空区进风侧氧化升温带从距工作面55 m移动到65 m；由于氧化升温带向深部转移，采空区煤自燃的风险会更加严重。