动态推进U 型采场煤自燃多场耦合数值研究

2023-10-26秦剑云宋双林王永敬李小超李世豪

煤矿安全 2023年10期

秦剑云，宋双林，王永敬，李小超，党龙，朱鹏，李世豪

（1.库车县榆树岭煤矿有限责任公司，新疆库车 842099；2.中煤科工集团沈阳研究院有限公司，辽宁抚顺 113122；3.煤矿安全技术国家重点实验室，辽宁抚顺 113122）

煤炭作为我国主体能源，其安全供给直接关系到国家能源安全[1]。煤矿火灾是威胁煤矿安全生产的主要原因之一，已严重威胁到煤矿企业的安全生产，其中煤炭自燃是引发煤矿火灾的最主要原因[2-3]，采动影响下的垮落岩体和遗煤在地下形成的多孔松散采空区是瓦斯与煤自燃致灾的重灾区，井下采空区遗煤自燃占井下火灾事故的60%[4]。

针对采空区的自然发火问题，国内外学者开展了一系列的研究，包括防灭火措施、自燃规律的数值模拟及实验分析等。在数值模拟方面，前期的研究以研究风流在采场的稳态分布规律为主[5-7]，根据求解的漏风速度对自燃区域进行大致区分。目前的研究多是针对采空区自然发火的多场耦合模拟研究[8-9]，将流场、温度场、浓度场、应力场等进行耦合，从矿井尺度上研究煤自燃高温区域的发生及发展的规律[10-12]。由于采空区的物理尺寸和物理特性都会随着工作面的推进而改变，一些学者通过引入移动坐标系来处理物理边界的变化[8,13]，或者利用CFD 软件的动网格技术[14]，得到孔隙度在时间和空间上的非均质分布，体现高温区域在小范围推进距离下的动态移动规律，模拟结果更符合实际。文献[15-17]通过COMSOL Multiphysics®与MATLAB 联用，在生成新的采空区网格后，调用上一时间步的场信息作为当前时间步的初值，实现模拟的物理量在时间和空间上的连续性，较大程度呈现了采空区煤自燃升温过程，但这些研究未将工作面与采空区的流动、传热及传质问题联合求解，从而也缺少围岩温度、通风温度等因素对采空区煤自燃升温过程影响的研究。为此，考虑采空区动态演化过程，以巷道出入口作为边界，将工作面与采空区的物理场联合求解，模拟工作面动态推进下采空区煤岩裂隙场中煤自燃灾害的发生、发展及演化规律；通过对采空区遗煤氧化升温情况的预测预报，可确定采空区煤自燃的危险性和高温热源的位置，进而为采空区这一隐蔽火区煤自燃的治理提供及时有效的手段。

1 模型构建

U 型采场可分为几个区域，即煤体固相区、采空多孔介质区和采煤作业区。煤体固相区包括面间煤柱、边界煤柱、待回采实体煤、顶板岩层和底板岩层；采煤作业区包括进风巷、工作面和回风巷。采空区二维示意图如图1。

图1 采空区二维示意图Fig.1 Two-dimensional schematic diagram of the mining area

1.1 基本假设

1）因采空区长度和宽度远大于煤层开挖引起的上覆煤岩层垮落高度，将采空区视为高度上平均的二维非均质多孔介质。

2）假设气体为理想气体，并忽略水蒸气对煤自热过程的影响。

1.2 采空区孔隙率及渗透率的动态演化

随着工作面向前推进，自由垮落带的矸石、遗煤在重力和矿压作用下会被逐渐压实，垮落碎胀系数在压实过程中不断减小，其数值满足式（1）[18]：

式中：Kp(x,y)为采空区垮落碎胀系数；Kp,max为初始垮落碎胀系数；Kp,min为垮落压实后的碎胀系数； αx、 αy分别为采空区碎胀系数倾向和走向方向的衰减率，1/m；dx、dy分别为采空区任意点(x,y)到采空区某一固定围岩边界和工作面边界的距离，m； ξ1为控制“O”型圈模型分布形态的调整系数， ξ1<1。

如果工作面推进速度为u，则随着时间推进，采空区内任意一点到采空区上下边界和工作面的距离分别可表示为：

式中：u为工作面推进速度，m/s；t为工作面回采时间，d；L为工作面长度，m；x0为工作面初始位置，m。

垮落采空区的高度H和孔隙率 εb的表达式分别为：

式中：M为工作面采高，m； εb为采空区孔隙率。

根据Ergun 的实验研究[19]，采空区渗透率K可以根据孔隙率、体积平均粒径计算得到，计算公式为：

式中：dp为采空区多孔介质的平均颗粒直径，m。

1.3 采空区流动及传热方程

采空区内渗透率具有明显的非均质特征，固壁处的渗透率和采空区深处的孔隙率相差1 个数量级，研究表明在靠近工作面的采空区内具有较大的空气流速，空气处于湍流状态，而采空区内部的空气流速较低，空气处于层流或者蠕动流状态，因此使用Darcy-Brinkman 方程描述采空区内气体的流动[20]。对于采空区内空气与煤岩体之间的换热过程，采用局部热平衡假设，即假定局部流体温度与煤岩体温度一致。采空区内的流动、传热及传质方程[21-23]为：

氧气消耗量与氧气浓度和煤自燃参数有关，煤氧反应热与氧气消耗量和反应热有关，分别计算如下：

式中：RO2为氧气消耗量，mol；Qa为煤氧仅应热，J；为氧气浓度，mol/m3；A为指前因子，s-1；E为活化能，J/mol；R为气体常数，J/(mol·K)；n为反应级数； ΔQ为煤氧化学反应过程中生热量，J/mol；sV为形状因子。

1.4 工作面及巷道控制方程

煤矿井下采场工作面和巷道内区域气体一般都处于湍流状态，因此需要使用湍流模型进行模拟。选用L-VEL 湍流模型[24]，该模型基于普朗特混合长度理论，根据局部流速和与最近壁面的距离来计算湍流黏度，稳定且计算效率高，适合用于计算内部流动。流动、传热及传质方程如下：

式中： µT为有效湍流黏度，kg·m2/s；I为压力在空向上的分布；ρ为流体密度，kg/m3；Cp为空气比热容，J/(kg·K)；λf为空气导热系数，W/(m·K)；Di为组分i的扩散系数，m2/s。

1.5 围岩（煤）传热控制方程

围岩（煤）固相传热以热传导为主，可用下式表示：

式中： (ρC)s为围岩（煤）热容，J/(m3·K)； λs为导热系数，J/(m·K)。

1.6 多场耦合关系

建立的模型具有3 个方面的耦合关系。首先，在采空区内部，采空区能量方程和气体组分方程相互耦合，并且均与采空区流动方程所求解的速度场相关联。其次，工作面和巷道气体传热方程和气体组分方程均与其内部速度场相关联。最后，温度、风速、压力、气体体积分数在采空区和工作面区域的界面上数值相等，即具有连续性。

1.7 工作面推进过程中计算网格的处理

随着工作面向前推进，采空区面积不断增大，计算区域物理尺寸不断变化，因此计算网格也应随之调整。利用COMSOL Multiphysics®软件的变形几何模型，通过设定计算边界的移动速度，实现网格的自由变形，从而实现对采空区动态增长过程的模拟。

1.8 初始条件及边界条件

通过模型耦合关系的分析，明确了采空区、面间煤柱、边界煤柱、底板岩层、顶板岩层、待回采实体煤、进风巷、回风巷等地点的初始条件和边界条件，如下：

进风巷及回风巷面间煤柱固壁、工作面固壁、始采线边界煤柱固壁和终采线边界煤柱固壁满足无滑移边界条件。

2 模型验证

利用在开滦集团崔家寨煤矿E12604 工作面获得的实测数据与建立的数学模型计算结果进行对比分析。崔家寨煤矿E12604 工作面推进距离570 m，倾向长度135 m，煤层厚度3.0～3.2 m，属易自燃煤层，最短自然发火期90 d，采煤方法为综采倾斜长壁、一次采全高，采高为2.8～3.0 m。现场数据监测点距离E12604 工作面开切眼120 m，距回风巷底板2.5 m。数值模拟参数为：①入口风速vin=1.01 m/s；②空气入口温度Tin=295.85 K；③活化能E=21.6 kJ/mol；④进风流O2浓度C0=9.375 mol/m3；⑤指前因子A=0.85 s-1；⑥氧化反应热400.11 kJ/mol；⑦围岩初始温度295.85 K；⑧煤岩密度1 300 kg/m3；⑨煤岩比热容1 003 J/（kg·K）；⑩煤岩导热系数λs=0.2 W/（m·K）；⑪初始碎胀系数Kp,max=1.5；⑫压实碎胀系数Kp,min=1.15；⑬碎胀系数衰减率αx=0.268 m-1；⑭碎胀系数衰减率αy=0.0368 m-1；⑮工作面采高3 m；⑯工作面推进速度7.5 m/d；⑰煤矸石直径dp=0.04 m；⑱工作面长度238 m。崔家寨煤矿E12604 工作面采空区气体温度和O2体积分数的模拟与实测结果如图2。

图2 模拟与实测数据对比Fig.2 Comparison of simulated and measured data

从图2 可以看出：在工作面回采初期采空区温度模拟结果与实测数据比较接近，随着推进距离的增加，测点位置的模拟温度值逐渐高于现场实测值，这可能是由于实际采空区靠近壁面附近存在邻近采空区的漏风源，导致实测温度偏低；O2体积分数的模拟值与实测值的变化规律一致，工作面回采初期O2体积分数较高，随着推进距离增加，采空区内部O2逐渐被煤氧反应所消耗，体积分数不断降低。总体上，可以认为建立的模型能够较好地描述采空区煤自热过程随工作面推进的动态变化。

3 工作面推进条件下多物理场的动态演化

工作面推进过程中采空区空隙率的变化如图3。

图3 工作面推进过程中采空区空隙度的动态变化Fig.3 Dynamic change of porosity in the mining area during the advance of the working face

从图3 可以看出：采空区空隙率呈现非均匀分布特征；总体上看，距离工作面和两巷煤壁较近的区域，由于液压支架及巷道煤柱的支撑作用，空隙率较大；距离工作面和两巷煤壁较远的区域，由于矿山压力作用明显，压实程度较高，空隙率较小；采空区空隙率分布规律具有典型的“O”形圈特征；随着工作面不断向前推进，采空区范围不断扩大，采空区走向长度增加，“O”形圈尺寸不断增大，但采空区空隙率分布规律始终保持不变；工作面回采时间分别为10 d 和60 d 时，采空区的空隙率最小值分别为0.18 和0.14，而在工作面液压支架附近空隙率始终保持在0.33 左右，这表明随着工作面不断推进，采空区内逐渐压实。

工作面推进过程中采空区温度场的动态变化如图4。

图4 工作面推进过程中温度场的动态变化Fig.4 Dynamic change of temperature field during the advance of working face

从图4 可以看出：随工作面推进距离的增加，高温区域的形状有较大变化；开采初期，采空区走向长度相对于倾向长度较短，采空区整体范围内氧气体积分数较高，因此升温区域及等温线分布明显地向回风侧延展；在动态采空区煤自燃的过程中，高温区域随工作面推进向前移动，且与工作面的相对距离保持不变；随着采空区走向长度的增加和氧化升温时间延续，采空区固相煤体温度上升且高温区域逐渐扩大。然而，由于采空区固相煤体与底板之间的换热效应，高温煤体在进入采空区深部窒息区后有一定的降温过程。因此，采空区高温区域的倾向宽度沿开切眼方向逐渐减小，并且由于热量的传递需要一段时间，存在“热惯性”现象，而采空区深部风流速度较小，热量不宜散失，因此在采空区的进风侧温度场形成拖尾现象。

停采后第10 d 的温度场分布如图5。

图5 停止开采后第10 d 温度场分布Fig.5 Temperature field distribution on the 10th day after mining stopped

从图5 可以看出：停采后，采空区内部由于缺少氧气供应，煤氧反应减弱，温度逐渐降低；由于工作面附近氧气供应充足，高温区域逐渐向工作面靠近，高温中心距离工作面约100 m，最高温度64.2 ℃。因此，对于工作面推进过程中形成的高温区域，如果没有形成灾变，那么在采空区的堆积压实和围岩散热作用下温度会逐渐降低，但开采停止后高温区域会向工作面迁移。

工作面推进过程中采空区氧气浓度场的动态变化如图6。

图6 工作面推进过程中采空区氧气浓度场的动态变化Fig.6 Dynamic changes of the oxygen concentration field in the mining area during the advance of the working face

从图6 可以看出：随工作面推进距离的增加，采空区中部及深部区域固相煤体逐渐被压实，漏风阻力增大，采空区氧气浓度呈现不对称性分布；进风侧的氧气浓度分布区域较回风侧更深入，氧气浓度的等值线分布与工作面的相对距离固定，采空区氧气浓度场分布趋于稳定。

停止开采后采空区第10 d 的氧气浓度分布如图7。

图7 停止开采后第10 d 氧气浓度场分布Fig.7 Oxygen concentration field distribution on the 10th day after the mining stopped

从图7 可以看出：停止开采后，氧气分布带形状变化不大，仍然呈不对称分布；进风侧宽度大于回风侧宽度，但总的宽度有所减小，这是由于采空区内部的氧气逐渐耗尽，而采空区内部压实，风流向采空区内部输送的氧气量非常有限；工作面封闭后，进风侧较回风侧积存的氧气多，煤氧化产生的热量不能及时带到回风侧，导致进风侧温度逐渐升高。

4 敏感性分析

当围岩温度为20 ℃时，不同通风温度下采空区内最高温度如图8。

图8 不同通风温度下采空区内最高温度Fig.8 Maximum temperature in the extraction zone at different ventilation temperature conditions

从图8 可以看出：首先在不同的通风温度下，随着时间推移，采空区内最高温度先迅速增加，后缓慢下降，这是因为在采空区推进的前期，由于采空区走向长度较短，风流阻力小，风速较大，采空区内氧气供应条件较好，煤氧反应剧烈，煤氧反应产热与散热相比占主导地位；随着采空区不断向前推进，采空区走向长度逐渐增加，采空区内风流阻力逐渐变大，风速减小，采空区内氧气供应条件变差，煤氧反应相对较弱，煤氧反应热小于煤岩散热，因此温度略有下降；其次，通风温度越高，采空区内最高温度越大，但最高温度的增幅要小于通风温度的增幅，且工作面推进的前期最高温度非常接近；通风温度为10、30 ℃时，第30 d 的最高温度分别为55.2、59.4 ℃，通风温度提升了10 ℃，但是最高温度的增幅小于5℃，这是由于空气在进入采空区前会与巷道壁面进行换热，进入采空区后温度降低，并且由于煤矸石的热容量远大于空气的热容量。

不同围岩温度下采空区内最高温度的演化如图9。