松散顶煤自燃流—固—热—化耦合数值模拟

张红芬，曹 伟，王宝夫，高尔新

(1.华北科技学院 建筑工程学院,北京 东燕郊 065201;2.中国矿业大学(北京) 力学与建筑工程学院,北京 100083)

0 引言

煤自燃是煤矿的主要灾害之一，是所有产煤国遭遇的棘手难题。以往井下局部自燃一般发生在采空区或破裂的煤柱内。随着综采放顶煤工艺的广泛实施，掘进巷道松散顶煤的高自燃危险性逐渐引起关注[1-4]。在煤自燃研究不断深入的过程中，人们逐渐认识到：它是由能量传输、气体流动、力学效应、氧气传输多个物理场的耦合作用推动的复杂物理化学过程[5-6]。2008年，杨胜强等指出了漏风与自燃之间的耦合关系[7]；2009年，余明高等对顶煤高温区位置和气体分布进行了分析，研究过程局限于二维稳态[8]；2012年，朱红青等求解了高冒区内漏风风速和氧气浓度的稳态和温度的动态分布[9]；2012年，文虎等进行了大断面全煤巷高冒区煤自然发火三维模拟[10]，得到了不同时间氧浓度场和温度场的分布，但没有考虑地质力学效应。总体来看，已有的研究中，多数忽略升温引起的力学效应，而进行三维动态研究的较少。在多场耦合问题中，任何物理过程或边界条件的改变都会影响其他物理过程[11]，可能造成对灾害的预测不准，贻误防治时机。本文结合工程实际建立巷道顶部松散煤体自燃三维模型，考虑力学效应、气体流动、传质场、温度场的耦合作用，获得物理场的变化过程，总结煤自燃过程高温区域、高位移区域位置和变化规律，为巷道顶部松散煤体自燃时空预测及防治提供依据。

1 工程概况

神华宁煤集团枣泉煤矿主采煤层为2#煤层，埋藏深度为200～300 m，煤层平均厚度为7.88 m。其11201工作面采用综放开采方式，工作面巷道为高3.56 m、宽5 m的矩形。巷道风温26～28 ℃，通风量1300～1500 m3/min。受矿压影响，顶煤破碎、离层，松散煤层厚度约5 m。顶棚网上堆积大量松散浮煤，多次险些发生火情。

2 松散顶煤自燃流—固—热—化耦合三维模型

2.1 几何模型

根据枣泉矿11201工作面的实际情况建立巷道顶部松散煤体自燃三维模型如图1所示，上部半圆柱代表巷道顶部松散煤体，下部中间为掘进后的巷道，左右两侧为煤壁，计算区长100 m。巷道风流从ZOY面开始沿X轴正向流动。模型ZOX面两边是对称的。网格划分先在横截面上进行，然后在巷道长度方向扫掠，划分结果如图1所示。

图1 计算模型与网格划分(m)

2.2 数学模型

假定煤是均匀的各向同性多孔介质，忽略瓦斯及水分影响。假设氧气主要通过对流和扩散，热量主要通过热传导和热对流进行传递，考虑煤的热力学效应。根据煤体的热、质传输过程中能量守恒原理和质量守恒原理，得到煤自然发火数学模型[5,12]：

流场方程：

(1)

(2)

氧气浓度场方程：

(3)

温度场方程：

(4)

(ρCp)whole=(1-n)ρccc+nρCp

(5)

kwhole=(1-n)kc+nkg

(6)

变形场方程：

(7)

式中，ρ为空气密度，取1.05 kg/m3；u为松散煤体内的流速，m/s；p为压力，Pa；T为煤的温度，K；F为体积力，N/m3；c为氧气浓度，mol/m3；D为氧气扩散系数，m2/s；R为耗氧速率，mol/(m3·s)；ρc为煤的密度，取为1350 kg/m3；CP、cc分别为空气、煤的热容，二者分别取为1005 J/(kg·K)，1540 J/(kg·K)；kc、kg分别为煤、空气的导热系数，分别取1.2 W/(m·K)，0.029 W/(m·K)；Q为放热强度，W/m3；E为煤的杨氏模量，1500 MPa；v为煤的泊松比，取0.4；K为体积模量，取为4900 MPa；u′为位移，m；αT为热膨胀系数，取8.7×10-6K-1；n为孔隙率，巷道顶部松散煤体取0.4，煤壁取0.15。模型中耗氧速率、扩散系数、放热强度使用煤自然发火试验结果[12]：

(8)

(9)

(10)

在上述模型中，流场受温度和压力的影响，氧气浓度场受流速和温度的影响，温度场受流速和氧气浓度的影响，变形场受流速和温度的影响，耗氧速率是温度和浓度的函数，氧气扩散系数是温度和压力的函数，放热强度是温度和浓度的函数，各物理场的耦合作用决定煤自燃发展。

2.3 边界条件

在进行巷道顶部松散煤体自燃模拟时边界条件设定为：巷道顶部松散煤体初始温度为围岩原始温度293 K；巷道风流初始温度为正常工作时巷道风温的平均值300 K；新鲜风流中氧气体积浓度为21%；巷道通风量为1300 m3/min，设进风口风速为1.2 m/s，风流出口按自由流设置；松散煤体上部边界为绝热壁，巷道壁为对流换热边界；巷道侧壁和顶部为自由面，其余为固定约束。

3 松散顶煤自燃流—固—热—化耦合数值模拟结果

3.1 温度场模拟结果

根据上述模型和边界条件，利用多物理场仿真有限元软件COMSOL进行巷道顶部松散煤体自燃三维流—固—热—化耦合数值模拟，借助PARDISO求解器获得了煤体自热升温过程各物理场的分布。图2、图3、图4为不同时刻横截面上的温度分布图，顶煤出现大范围破碎10天后，在煤体中形成了314 K的区域，第17天在煤体中形成了388 K的区域，已接近其临界温度393 K，第21天时煤体局部温度高达510 K。可见，松散顶煤经历了以临界温度为界由缓慢升温到急剧升温的过程。因此，在煤体发生破碎后应密切关注温度变化，尽力在临界温度前采取防治措施，临界温度后升温迅速，难以控制。根据文献[13]中的现场观测数据，巷道顶部出现大范围破碎现象10天后，巷道表面温度升高到30 ℃左右，顶煤中下部温度达到40 ℃左右，数值模拟结果与现场实测数据吻合。

图2 第10天横截面上的温度分布

图3 第17天横截面上的温度分布

图4 第21天横截面上的温度分布

由图2、图3、图4可见，在横截面上，松散顶煤中贴近巷道部分温度最低，顶部附近次之，高温区域位于顶煤中下部，呈现明显的三带分布。自热升温初期，高温区域位于距巷道顶部1.5～3 m的位置，随着时间的推移煤体高温区逐渐向下移动，即向进风侧移动。

图5、图6为纵截面上的温度分布图，可见高温区域位于距离巷道进风口大约5～20 m的范围，随着温度的升高，高温区域逐渐向进风口移动。

图5 第17天纵截面上的温度分布

图6 第21天纵截面上的温度分布

3.2 位移场模拟结果

图7为数值模拟计算得到的第21天时巷道顶部松散煤体位移分布，在煤氧复合自热升温过程中，松散顶煤中下部向巷道内位移，不断挤入巷道。第10天时最大位移为1.3 cm，第17天时最大位移为4.6 cm，第21天时最大位移达到了6.3 cm，随着温度升高，位移数值不断增大。松散顶煤升温在横截面上形成以巷道顶为中心向周围扩展的位移分布，位移由中心向外逐渐降低，高位移区域在松散煤体下部，巷道顶面上的位移最大。在巷道顶面上，高位移区域位于距离进风口5～20m的范围，和高温区位置相同。受松散顶煤向巷道内位移的影响，与松散顶煤接触的一部分煤壁发生了一定程度的位移，第21天时煤壁内存在平均位移3 cm左右的区域。因此在关注松散顶煤高温高位移问题同时，也要注意煤壁位移可能带来的破坏。

图7 第21天的位移分布

在煤自热升温过程中，松散顶煤高温区域位于松散煤体中下部，松散煤体最下部温度最低，使得自燃现象难以被及时发现。位移场则不同，位移最大区域位于松散煤体最底部，容易观察和发现异常。因此在进风口附近的巷道顶部布置位移观测点有助于及时发现自燃现象，尽早采取处理措施。

4 结论

(1) 数值模拟计算结果较好地预测了松散顶煤自燃过程，与现场实测数据吻合较好。

(2) 巷道顶部煤体破碎后，经历了以临界温度为界由缓慢升温到急剧升温的过程。温度场呈三带分布：贴近巷道部分温度最低，顶部附近次之，高温区位于松散顶煤中下部。氧化初期高温区在横向截面上位于距巷道顶部1.5～3m的位置，在纵向截面上位于距进风口5～20m的位置。随温度升高，高温区在横向和纵向截面都逐渐向进风侧移动。

(3) 在升温过程中，巷道顶部松散煤体中下部向巷道内位移，位移值随温度升高而增大。在横截面上形成以巷道顶为中心向周围扩展的位移分布，由中心向外位移逐渐降低，松散煤体下部为高位移区。在巷道顶面上，高位移区域位于距离进风口5～20 m的范围。

(4) 防治建议：进风口附近巷道顶部松散煤体中下部为煤自燃防治重点区域；巷道顶面上的位移最大，便于观察和发现异常，可在巷道顶部进风口附近布置位移观测点，以便及时发现自燃现象。