综放采空区瓦斯与煤自燃多场耦合模拟研究*

贾廷贵，娄和壮，刘 剑，曲国娜

(1.辽宁工程技术大学 安全科学与工程学院，辽宁 阜新 123000；(2.内蒙古科技大学 矿业与煤炭学院，内蒙古 包头 014010)

0 引言

我国煤矿受瓦斯与煤交织影响发生共生灾害威胁越发凸显，这种共生灾害发展成为煤矿特重大事故导火索的致灾模式，正逐步蚕食我国煤矿安全生产长效机制，遏制了煤炭工业持续健康发展。针对瓦斯与煤自燃双重威胁，文献[1-6]对瓦斯与煤自燃复合致灾条件与灾变演化进程进行了深入研究，针对采空区煤氧共存环境，得出了复合致灾的显著特性；文献[7-9]对瓦斯与煤自燃复合致灾协同防治进行了充分研究，取得一定成果。本文结合工作面受采动影响引起的煤体变形、裂隙瓦斯混合气体流动、气体组分渗流与热量传输等多物理过程，建立瓦斯抽采诱导煤自燃的热-流-化模型多场耦合数学方程，分析了瓦斯与煤自燃多场耦合致灾特性，为瓦斯与煤自燃的协同防治提供指导。

1 理论模型

1.1 气体流动控制方程

假设瓦斯与空气混合气体为不可压缩理想气体，对于煤层瓦斯，其吸附与解吸满足朗格缪尔等温吸附方程[10]，则有质量守恒：

(1)

m1=βp1φ

(2)

(3)

根据非达西福希海默尔方程[11]有瓦斯与空气混合气体在采空区流动动量守恒：

(4)

(5)

式中：m为气体质量，kg；m1为空气质量，kg；m2为瓦斯气体质量，kg；ρ为气体密度，kg/m3；ρs为煤的密度，kg/m3；ρ2a为标准状况下瓦斯气体密度，kg/m3；β为压缩因子，kg/(m3·Pa)；VL为朗格缪尔常数，m3/kg；p1为空气压力，Pa；p2为瓦斯气体压力，Pa；q为气体流速，m/s；μ为气体动力黏度，kg/m·s；g为重力加速度，m/s2；S为源汇项；φ为孔隙率，%；Kp为渗透率，m2。

1.2 气体组分质量守恒方程

瓦斯混合气体满足质量传输守恒：

(6)

(7)

(8)

式中：c为气体浓度，mol/m3；D为气体分子扩散系数，m2/s；Ic为源汇项；c(O2)为O2气体摩尔浓度，mol/(m3·s)；c0为参考O2气体摩尔浓度，mol/m3；γ0为参考O2气体消耗速率，mol/(m3·s)；α为氧化温度指数，℃-1；H为采高，m；d0为采空区某点到工作面距离，m；wa为采空区局部瓦斯涌出强度，wa=2.3×10-4mol/(m2·s)；wb为采空区下部均匀瓦斯涌出强度，wb=1.5×10-4mol/(m2·s)；w1为本煤层瓦斯涌出强度，w1=7.0×10-4mol/(m2·s)；w2为上部瓦斯涌出强度，w2=3.2×10-4mol/(m2·s)；λ为瓦斯释放衰减系数，λ=0.076 d-1。

1.3 渗透率变化方程

采空区的孔隙率与渗透率逐渐变化，影响着采空区瓦斯混合气体流动[12]，在实际采空区中，某点孔隙率和渗透率与该点的空间坐标相关，推导得到计算式为：

(9)

(10)

式中：Kp为渗透率，m2；dp为煤颗粒直径，m；L为工作面倾向长度，m；x,y,z为空间坐标。

1.4 热量传输守恒方程

假设煤自燃过程无矸石与混合气体热传递，忽略采空区水蒸气热传递[13]，则有采空区固相气相热量传输守恒：

(11)

(12)

(13)

(14)

式中：ρg为气体密度，kg/m3；cp(s)为固相比热容，J/(kg·℃)；cp(g)为气相比热容，J/(kg·℃)；κs为固相导热系数，W/(m·℃)；κg为气相导热系数，W/(m·℃)；IT(s)为源汇项；hsg为固相气相界面导热系数，W/(m·℃)；asg为采空区多孔介质比表面积，m-1；Ts为固相温度，℃；Tg为气相温度，℃；ΔH为煤氧化反应热，J/mol；Ic(O2)为O2源汇项。

上述控制方程定义了采空区瓦斯与空气混合气体流动、气体组分对流扩散、受采动影响下采空区煤体渗透率变化与热量传输等多物理过程数学模型，其相互之间耦合关系见图1。

图1 综放采空区瓦斯与煤自燃多场耦合作用Fig.1 Multi-field coupling effect of gas and coal spontaneous combustion in goaf of fully mechanized caving face

2 瓦斯与煤自燃多场耦合数值模拟

2.1 物理模型

寸草塔二矿31102综放工作面位于31号煤一盘区，工作面推进长度1 179 m，工作面倾向长度186 m，煤层倾角为1°，平均采高3 m，放煤高度3 m。31102综放工作面采用“U”型通风，风量2 255 m3/min。31102综放工作面为自燃煤层，煤自燃等级为Ⅰ类。建立数值模拟模型，工作面倾向长度×宽度×高度为186 m×5 m×6 m，采空区走向长度×宽度×高度为200 m×186 m×60 m，运输巷与回风巷长度×宽度×高度为80 m×4 m×3 m，工作面后方未采的实体煤长度×宽度×高度为80 m×186 m×11.5 m，未采区域煤层以上均设置为岩层。

根据式(15)～(16)确定的冒落带高度H1与裂缝带高度H2分别为12，37 m。

(15)

(16)

为减少模拟计算量，仅对单个高位钻孔进行研究，在回风侧设置1号钻场，钻孔仰角α1=10°，方位角β1=30°，钻孔长度L=140 m，开孔位置距离底板2 m，终孔位置距离回风巷巷帮70 m，距离底板25.3 m。几何模型见图2(a)，网格划分结果见图2(b)，网格处理得到45 835个域单元、8 659个边界单元和1 526个边单元。

2.2 边界条件与初始条件

基于寸草塔二矿31102综放工作面相关物理参数进行模拟研究，相关物理参数见表1。

设置进风速度v_inlet=2.5 m/s，设置进风巷与回风巷两端间压差110 Pa，原煤瓦斯压力0.4 MPa。高位钻孔设置为质量流出边界，抽采流量60 m3/min，抽采负压2 MPa。设定初始条件：c(O2)t=0=0 mol/m3，Ts=Tg=22 ℃，设定采空区与工作面交界为浓度控制边界与温度控制边界，其他为无通量边界。

图2 高位钻孔抽采瓦斯模型Fig.2 Model of high-level borehole gas drainage

物理参数符号数值煤的密度/(kg·m-3)ρ01 320瓦斯朗格缪尔压力常数/MPaPL6.02瓦斯朗格缪尔常数/(m3·kg-1)VL0.015瓦斯气体密度/(kg·m-3)ρ2a0.716气体动力黏度/(kg·(m·s)-1)μ0.3气体分子扩散系数/(m2·s-1)D2×10-5固相比热容/(J·(kg·℃)-1)cp(s)1 002气相比热容/(J·(kg·℃)-1)cp(g)1 012固相导热系数/(W·(m·℃)-1)κs0.2气相导热系数/(W·(m·℃)-1)κg0.026煤氧化反应热/(kJ·mol-1)ΔH400煤颗粒直径/mdp0.04

3 数值模拟结果分析

模拟得到高位钻孔仰角10°，方位角30°，钻孔长度140 m，抽采流量60 m3/min条件下采空区瓦斯浓度分布见图3，采空区冒落带与裂缝带最高瓦斯浓度变化见图4。模拟高位钻孔瓦斯抽采初始阶段，图3(a)抽采时间24 h时，靠近回采工作面的采空区瓦斯浓度得到有效降低，沿采空区走向，进风侧瓦斯浓度比回风侧降低幅度大，即在距离工作面同等长度条件下，进风侧瓦斯浓度比回风侧低。采空区深处瓦斯浓度依旧较高，高位钻孔抽采瓦斯使得采空区漏风从进风侧流入，回风侧流出，漏风风流流经采空区降低了靠近工作面的采空区浓度，但漏风风流存在沿程损失，采空区深处瓦斯仍将聚集。模拟得到抽采时间24 h冒落带与裂缝带瓦斯平均浓度分别为63.97%，59.55%。

图3 采空区瓦斯浓度分布Fig.3 Distribution of gas concentration in goaf

图4 采空区瓦斯平均浓度变化Fig.4 Change of average gas concentration in goaf

随着模拟抽采时间达到120 h，见图3(b)，沿采空区走向回风侧存在较高浓度瓦斯，相比较抽采时间24 h时瓦斯浓度明显延伸到采空区深处，同时进风侧仅存较低浓度的瓦斯，可见高位钻孔抽采瓦斯效果明显。模拟得到抽采时间达到120 h的冒落带与裂缝带瓦斯平均浓度分别为25.59%，18.14%。在模拟范围内，冒落带与裂缝带瓦斯平均浓度分别降低了68%，77%，将高位钻孔布置在裂缝带，可有效提升瓦斯抽采效率。

抽采时间24 h与120 h回风侧瓦斯浓度分别为0.86%，0.18%，见图5，瓦斯抽采初始阶段回风侧瓦斯浓度较高，采空区漏风携带着高浓度瓦斯涌向回风隅角，但随着抽采的进行，回风侧瓦斯浓度在可控范围内。由于受工作面采动影响，采空区上覆荷载造成煤岩体发生移动断裂垮落，当工作面形成一定范围时，采空区中部垮落煤岩体逐渐被压实，但压实区周围即巷帮位置仍将保持堆积状态[14-16]，采空区的这种非均质性为采空区漏风提供有效渠道。漏风风流由进风侧流入，在抽采负压作用下，采空区高浓度瓦斯向回风侧诱导流出，在回风侧形成瓦斯汇聚区域，随着抽采时间的进行，采空区外瓦斯浓度降低。

图5 回风侧瓦斯浓度变化Fig.5 Change of gas concentration on return air side

通过在31102综放工作面现场实践，设置高位钻场参数为上述2.1节模拟参数，采集1号高位钻场实测回风巷数据，对比模拟结果发现二者具有较强一致性。通过监测抽采期间工作面上隅角瓦斯浓度在120 h稳定在0.21%左右，保证了工作面的安全回采。

选择模型z=1.2 m平面进行采空区O2浓度与温度研究，模拟得到抽采时间120 h的采空区O2浓度分布、温度分布见图6。采空区进风侧O2浓度18%时宽度为66.7 m，8%时宽度为128.9 m，以O2浓度8%～18%为氧化带宽度，得到进风侧氧化带宽度62.2 m，回风侧氧化带宽度31.5 m。最大氧化带边界位于进风侧，即距离工作面同等距离条件下，进风侧O2浓度高于回风侧，这是由于进风侧漏风导致的。模拟得到抽采时间120 h的采空区最高温度为46.2 ℃，靠近进风侧，其距离工作面59.4 m，距离进风侧巷帮22.3 m。沿采空区走向方向上温度先增加后减少，形成高温区域，温度45 ℃及以上范围距离工作面41～77 m左右，距离进风侧巷帮约11～51 m。此区域范围基本覆盖氧化带，由于氧化带内O2浓度适宜煤氧复合反应的发生，随着瓦斯抽采的进行，煤体蓄热进程加快，采空区煤自燃危险增加。

图6 采空区氧气分布、温度分布Fig.6 Distribution of oxygen and temperature in goaf

通过在距离工作面59 m与60 m，距离进风侧巷帮20 m处埋设2个温度测点，捕捉到抽采时间120 h的测点温度为47.3 ℃与47.0 ℃，与模拟结果误差为2.38%与1.73%。

通过改变抽采流量，得到不同抽采流量条件下的回风侧瓦斯浓度变化、进风侧氧化带宽度与采空区最高温度见图7。在模拟抽采时间范围内，抽采流量分别为40，60，80，100 m3/min条件下回风侧瓦斯浓度逐渐降低，抽采时间120 h时瓦斯浓度分别为0.34%，0.18%，0.10%，0.14%，即随着抽采强度的增加，回风侧瓦斯浓度先降低后上升。

图7 模拟抽采流量影响抽采效果Fig.7 Influence of simulated drainage flow on drainage effect

抽采流量40，60，80，100 m3/min条件下抽采时间120 h的进风侧氧化带宽度分别为51.8，62.2，71.0，75.9 m，采空区最高温度为40.1，46.2，51.8，61.6 ℃。随着抽采强度的增加进风侧氧化带宽度增大，煤自燃范围增大，采空区最高温度逐渐上升，采空区煤体自燃可能性增强。加大抽采流量增大采空区漏风，扩大了采空区氧化带宽度，煤氧复合放热造成煤体温度升高，产生有毒有害气体在风压的推动下沿风流方向流动，涌出工作面。当瓦斯气体中混入CO气体会扩大瓦斯爆炸浓度上下限[17-19]，促使瓦斯爆炸，因此在提升高效抽采瓦斯的同时会降低采空区安全性。

4 结论

1)建立了高位钻孔瓦斯抽采诱导采空区煤自燃的热-流-化模型多场耦合数学模型，采用COMSOL软件模拟了寸草塔二矿31102综放采空区高位钻孔抽采瓦斯诱导煤自燃过程，对比模拟结果与现场测试结果吻合程度较高，模型对研究高位钻孔瓦斯抽采诱导采空区煤自燃具有可靠性。

2)高位钻孔布置在裂缝带有效地提升瓦斯抽采效率，在瓦斯抽采参数如钻孔仰角10°，方位角30°，钻孔长度140 m，抽采流量60 m3/min条件下抽采时间120 h的采空区进风侧氧化带宽度由于漏风的存在大于回风侧，采空区最高温度靠近进风侧为46.2 ℃。

3)综放采空区高位钻孔抽采瓦斯强度与煤自燃相互交织相互影响，加大抽采强度降低了回风巷瓦斯浓度，增大了进风侧氧化带宽度，提高了采空区最高温度，降低采空区安全性。瓦斯与煤自燃二者之间的矛盾在于，采场得到高效抽采瓦斯的同时，会造成采空区氧化带边界向深处蔓延，扩大煤自燃高温区域，漏风携氧充分的参与煤氧复合反应，采空区最高温度逐渐上升，煤自燃风险增大。