Y型通风方式瓦斯与煤自燃耦合模拟

贾廷贵，李颜兵，娄和壮，李绪萍

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

0 引言

当工作面漏风风流从采空区汇入回风顺槽时，受风流裹携作用，采空区中的瓦斯极易在工作面上隅角形成积聚，对工作面造成了极大的威胁[1-3]，改变工作面的通风方式、供风量可以有效地降低工作面上隅角瓦斯体积分数，但是增加工作面供风量会导致采空区氧化带宽度的增加，进而增加了防灭火工作的难度，合理地调整工作面供风与埋管瓦斯的抽采流量可有效解决工作面上隅角瓦斯积聚问题.寸草塔二矿31102综放工作面采用Y型通风方式，在工作面推进过程中，受到采动和采空区顶板冒落影响，工作面隅角瓦斯超限的情况时有发生，并且高瓦斯区域与氧化高温区域重合极易发生煤自燃[4-5].

周福宝[6]对瓦斯与煤自燃多场耦合进行研究并建立多元气体运移模型，分析并验证了多场交汇的致灾机理.杨胜强[7]等对采空区氧化产物与采空区温度的变化进行了研究，认为氧化产物CO、CH4生成初始温度不同，总体上呈现“滞缓效应”.白振龙[8]等对采空区中孔隙分布与顶板垮落特征进行了研究，确定矿压变化与采空区孔隙分布规律.张发亮[9]对采空区三带进行相关研究，确定了煤与氧的共存环境与采空区中易发生自燃的区域.何磊[10]等使用FLUENT对Y型工作面与U型工作面瓦斯运移特征进行对比，表明Y型工作面可以有效地降低工作面隅角瓦斯体积分数.张东坡[11]针对采空区易发生自燃区进行了注氮研究，并取得良好的防灭火效果.余陶[12]等在易自燃和大采高工作面注CO2进行防灭火工作.王省身[13]等对采空区通风方式与瓦斯运移规律进行研究，并证明Y型通风方式可以减少工作面上隅角的瓦斯积聚.周西华[14]等使用COMSOL对封闭火区的流场进行模拟.秦波涛[15]等提出多种瓦斯与煤自燃的治理方法.裴晓东[16]等在治理瓦斯时对采空区自燃进行了定性分析.

采空区遗煤自燃与瓦斯治理是典型的流-固-热-化耦合过程的结果，是采空区气体的传输-运移场、应力场、煤-氧复合化学场、煤-空气之间的热质等多场耦合共同作用.寸草塔二矿的瓦斯绝对涌出量为27.95 m3/min，相对瓦斯涌出量为0.71 m3/t.本文针对寸草塔二矿31102工作面建立瓦斯气体在采空区裂隙中流动、气体渗流与传热等多种物理过程的流-化-热多场耦合模型，分析在治理工作面瓦斯超限的情况下对采空区氧浓度场与温度场的变化规律，为后续瓦斯与煤自燃共同治理提供参考.

1 数学模型

1.1 采空区流场数学模型

在采空区中，气体的流动是层流、过渡流和紊流共存的状态，仅采用达西定律来描述采空区内气体流动规律与实际存在偏差.1967年，BACHMAT提出的多孔介质中大雷诺数情况下渗流的三维运动方程为

式中，E为渗透率，二阶张量，m2；J为压力梯度，Jh=-▽，h为静水压头，m；ν为动力粘度，m2/s；g为重力加速度，m/s2；q为瓦斯涌出量，m3/s；β为介质颗粒形状系数，取1.5；d为平均调和粒径，m；n为多孔介质的空隙度，%；n=p/γ，γ为容重，N/m3；qs为比流量，L/（m·s）.

采空区的风流为紊流时，使用不可压缩流纳维斯托克斯（Navier-Stokes）方程为

式中，ρ为空气密度，kg/m3；u为工作面风速，m/s；p为工作面压力，Pa；F为体积力，N/m3；I为源汇项，包括生成项与漏入项；▽为梯度算子.

采空区风流为过渡流时，采用布林克曼（Brinkman）方程为

式中，η为空气粘性系数，Pa·m；εp为多孔性，取1；k为热导率，H/m.

采空区中的风流为层流时采用达西定律为

式中，p为初始压力，Pa；W为源项，kg/（m·s）.

1.2 采空区温度场数学模型

热源控制区域为工作面采空区及上覆采空区遗煤存在区域，该区域既有氧化产热，又存在散热、传热、辐射等现象为

式中，dz为热源场的计算高度，即浮煤厚度，m；Cρ为流体介质比热容，J/（kg·K）；keff为有效导热系数，W/（m·K）；Q为氧化产热量，W/m3；Qρ为气体流动传输的热量，W/m3；Qoop为面外热辐射，W/m3；T为工作面温度，K.

采空区中的传热区域为工作面采空区及上覆老空区遗煤部分，传热控制方程为

式中，ρp为多孔介质密度，kg/m3；Cρ,p为多孔介质比热容，J/（kg·K）；θp为相对体积分数.

1.3 物理模型的建立与边界条件

寸草塔二矿31102工作面长度为320 m，煤层平均厚为5.7 m，平均采煤高为3.7 m，煤层倾角1°～3°，设定工作面的尺寸为4 m×320 m×4 m，采空区的模型见图1，尺寸为320 m×320 m×50 m，未采区域煤层以上均设置为岩层，进、回风巷的尺寸为 150 m×4 m×5 m，坐标原点位于右下角，工作面方向为x，工作面倾向方向为y，指向顶板方向为z，采空区的模型采用自动网格划分方法，生成节点 54 892个，单元294 874个，时间步长为1 h.计算过程中首先计算材料坐标系下的各项物理参数，再开展采空区风压、采空区漏风流场的计算，其后依次开展采空区温度场分布以及采空区浓度场等计算，每1 h计算完成保存一步计算结果，并开始进入下一个时间步长的计算，依次循环至设定求解结束时间停止数值模拟计算，模拟中的物理参数见表1.

图1 采空区网格模型 Fig.1 model of goaf grid

表1 模拟参数 Tab.1 simulation parameters

2 模拟结果分析

2.1 不同风量下瓦斯分布

寸草塔二矿31102工作面的实际供风量为1 174 m3/min，将供风比b设置为1，在供风比为0.5、0.75、1、1.25、1.5时采场瓦斯分布见图2.由图2可见，当供风比从0.5增加到1时，采空区内瓦斯体积分数下降明显；当供风比从1增加到1.25时采空区瓦斯体积分数下降较低；当工作面的供风比从1.25增加到1.5时，采空区的瓦斯体积分数出现了反弹，采空区中最高瓦斯体积分数增加了0.1%.随着工作面的供风比的增加，采空区中的高瓦斯区域逐渐向采空区后方移动.

图2 不同风量下采场瓦斯体积分数分布 Fig.2 gas volume fraction distribution of stope under different air volume

工作面倾向距离瓦斯体积分数变化见图3，随着距离进风侧距离的增加，工作面的瓦斯体积分数也逐渐增加，在距工作面回风侧20 m处瓦斯体积分数开始大幅下降，Y型通风方式有效地降低工作面隅角瓦斯体积分数.当工作面的供风比从0.5增加到1.25时，工作面最大的瓦斯体积分数从1.88%下降到0.24%；当供风比从1.25增加到1.5时，工作面的瓦斯体积分数开始上升，工作面最大瓦斯体积分数增加了0.25%.当供风比从0.5增加到1.25时，风流带走的瓦斯量比采空区涌出的瓦斯量小，导致工作面瓦斯体积分数降低；当供风比增加到1.5时，采空区的漏风量增加，采空区瓦斯的涌出量大于风流带走的瓦斯量，导致工作面瓦斯体积分数上升.

图3 工作面倾向距离瓦斯体积分数变化 Fig.3 variation of gas volume fraction in inclined distance of working face

2.2 采空区三带及高温区域分布

工作面供风比为1，工作面推进速度为4 m/d时，采空区氧化带内氧气体积分数为8%～18%，进风侧氧化带范围为117～414 m，中部氧化带范围为141～442 m，回风侧氧化带范围133～401 m.进风侧比回风侧氧化带范围宽27 m，采空区中部氧化带宽度为301 m，中部氧化带范围最宽.从氧化带宽度可以看出，由于工作面的通风方式以及采空区漏风的影响，使回风侧氧化带的增加.

在供风比为1工作面推进速度为4 m/d时，采空区温度分布见图4，工作面的位置在图4中采空区的左侧，采空区最高温度点的温度为51.4 °C，采空区高温区域靠近进风侧，距工作面的距离为 352 m.从图5中可以看出采空区温度上升速率比采空区温度下降速率小.由采空区温度分布与采空区氧化带分布情况可知，采空区高温区域靠近氧化带下边界区域，在采空区中满足遗煤氧化的条件与蓄热条件时，遗煤升温氧化需要足够的蓄热条件和一定的时间才能达到采空区最高温度.随着距工作面距离的增加，氧气的体积分数显著降低，采空区的温度也迅速降低.采空区最高温度点在氧化带偏后靠近氧化带与窒息带区域.

图4 采空区温度分布 Fig.4 temperature distribution of goaf

2.3 自然发火期的影响因素

自然发火期是煤矿开采后煤矿发生自燃的最短时间，寸草塔二矿31102工作面采空区遗煤自然发火期与推进速度与工作面风速的关系见图5.

图5 采空区自然发火期变化 Fig.5 changes in spontaneous combustion period of goaf

在工作面推进速度一定的情况下，随着工作面风速的增加，寸草塔二矿31102工作面采空区的自然发火期变短.当风速增加时，工作面的漏风量逐渐增加，将大量的氧气带入采空区，使采空区遗煤发生自燃的区域面积增加，使采空区的自然发火期的时间缩短.在工作面风速一定时，采空区推进速度从3 m/d增加到12 m/d时，采空区的自然发火期增加.在工作面风速从0.8 m/s增加到2 m/s时，采空区的自然发火期下降明显，当工作面风速继续增加时，自然发火期变化基本稳定.

2.4 埋管抽采对采空区三带的影响

为解决工作面上隅角瓦斯超限的问题，决定使用埋管抽采的方法进行瓦斯治理.随着抽采瓦斯流量的增加，采空区氧化带的范围逐渐增加，当瓦斯抽采流量从30 m3/min增加到120 m3/min时，采空区氧化带范围从278 m增加至321 m，当瓦斯的抽采流量增加到150 m3/min时，采空区氧化带的宽度增加到326 m.随着抽采流量的增加，氧化带的宽度增加的速度逐渐变缓.在抽采瓦斯之前，采空区属于“一源一汇”，当进行埋管抽采瓦斯时，采空区变成“一源两汇”，采空区回风侧的负压增加，导致采空区的漏风量增加，采空区氧化带范围增加.

3 工程应用

在供风比为1.25、工作面推进速度为4 m/d、瓦斯抽采流量为60 m3/min时，工作面倾向方向瓦斯体积分数变化见图6，通过对比可以看出模拟结果与实际结果相差不大，误差在5%以内，在工作面倾向方向瓦斯体积分数都出现了先增大后减小的趋势，工作面上隅角瓦斯体积分数降低到0.16%，有效地解决了工作面上隅角瓦斯体积分数超限的问题.

图6 工作面瓦斯体积分数变化 Fig.6 change of gas volume fraction in working face

4 结论

（1）建立流化热多场耦合模型，采用COMSOL模拟软件模拟了Y型通风工作面不同供风量条件下采场瓦斯体积分数变化与氧气浓度变化，分析了其对采空区三带的影响，采区埋管抽采瓦斯方式降低采场瓦斯体积分数，现场应用效果明显.

（2）工作面风速的降低与推进速度的增加都会使采空区的自然发火期增加.采空区温度场相比于采空区三带具有滞后性.

（3）埋管瓦斯抽采流量的增加，使采空区的氧化带宽度变宽，在工作面推进速度为4 m/d、工作面供风量为1 467 m3/min、抽采流量为60 m3/min时，工作面上隅角瓦斯体积分数降至0.16%.