三带数值模拟报告0502

1、51121(3)采空区“三带”数值模拟采空区自然发火防治的首要技术前提是要搞清采空区自燃“三带”的分布范 围，它是采空区注氮、堵漏风等技术措施的主要技术依据。因此采空区“三带” 的划分显得尤为重要。5. 1 1121(3)工作面概况1121(3)工作面位于东一采区13-1煤层第二块段，煤层呈块状及粉末状，煤 层产状变化较大，煤层倾向40左右，倾角1342。，煤厚3.06.0m,平均 煤厚4.4m。工作面平均长130m，宽5.256.05m，高3.24.0m。上风巷标高： -448.0m，下顺槽标高：-566.0moX作面顶、底板岩性如表5-1所示。该工作面南侧有1111 (3)工作面，与111

2、1 (3)工作面下顺槽留设10m煤 柱，与1121 (1)工作面相邻，高差均在80100m左右。其中1111 (3) 已回采 完毕，1121 (1)正在回采。工作面采用区内后退走向长壁式，一次采全高综合 机械化采煤法，平均推进速度为93m/月。根据精查地质报告，本矿井恒温 带深度为30m，温度16.8C，地温梯度2.3C/100mH作面实际温度在2630C。 工作面实际供风量为1800 m3/min。该煤层具有自然发火性，自然发火期36 个月。瓦斯相对涌出量为2.52m3 /t，绝对瓦斯涌出量为8.84m3 /min，矿压为 21-24mpa。表5-1工作面顶、底板岩性岩石种类备注老顶细砂岩厚

3、度为07.2/5.6m，灰白色、细粒为主直接顶泥岩伪顶炭质泥岩厚度为00.4/0.2m，灰黑色，主要为炭质泥岩直接底泥岩老底细砂岩、中细砂岩厚度为10.018.0/16.5m，浅灰灰白色，坚硬5. 2采空区“三带”划分指标工作面正常生产时，采空区自燃“三带”处于一个动态的稳定状态。采空区 自燃“三带”主要指散热带、氧化带和窒息带。散热带内由于冒落不充分，漏风 流较大，采空区遗煤气化产生的热量不能积聚，一般不会发生自然发火。氧化带 内漏风风速适当，具有热量积聚的条件，O2体积分数又能满足氧化需求，因此最 容易发生自然发火。窒息带内由于漏风难以到达，O2体积分数往往难以满足氧化 需求，一般也不会发

4、生自然发火。“三带”是客观存在的，但如何划分具有一定的困难。由于探测手段和方法 的局限，想要定量地准确划分是难以做到的。目前，一些研究者提出确定划分“三 带”的指标有漏风风速(V)、采空区氧浓度和温升速率3种：1、采空区漏风风速(v)。从理论上说，漏风风速相对较好。因为它可以体现 氧浓度分布、氧化生热与散热的平衡关系。B.M.Maebcka(1946)等研究表明，采 空区及煤柱的漏风强度在0.10.24 m3(nin m2，寸容易自然发火；有的研究者认 为不会导致自燃的极限风速低于0.02 0.05m3： m2；封闭采空区密闭墙漏风 压差在300Pa、漏风强度在0.02 1.2m3： inin

5、- m21寸容易自燃发火。由于这些参 数都是在一定条件下取得的，因此有一定局限性，但对研究自燃问题是有一定的 参考价值。2、采空区氧浓度(c)分布。采用氧浓度指标不能划分散热带和自燃带。因为 在自燃带中氧浓度也有可能达20%以上。因为氧化初期，耗氧量是很小的。对 于划分自燃带和窒息带的指标有的研究者认为是5%6%；有的认为是8%，即 氧浓度c8%为自燃带。实际上，氧化速度随氧浓度降低而减 小，究竟氧化速度小到何值才算之是“窒息”是没有确切标淮；3、采空区遗煤温升速度(dt1OC/d为自燃带)也是一个划分“三带”的指 标，但目前应用不是很广泛。按照采空区内漏风风速指标划分可知，散热带：采空区内漏

6、风风速大于 0.24m/min的区域；氧化带：漏风风速在0.240.1m/min之间的区域；窒息带： 漏风风速小于0.1m/min的区域。考虑到1121(3)工作面的推进速度快，采空区浮煤还没有足够的氧化时间， 氧气浓度随工作面的推进，变化并不明显。因此，本章主要采用专业的流体力学 模拟软件FLUENT对1121(3)工作面采空区进行流场模拟，并根据采空区内漏风 风速划分指标来确定1121(3)工作面采空区“三带”的位置及范围。最后改变工 作面的供风量，得到工作面不同供风条件下的采空区“三带”位置及范围。5.3物理模型的建立根据潘北矿1121(3)工作面现场实际提供的相关资料，为了使模拟的结果

7、更 能够反映工作面及采空区的实际流场分布情况，把上、下顺槽，工作面和采空区 三者结合起来一起作为采场整体来进行数值模拟。考虑到现场实际的上、下顺槽及采空区的走向很长，而我们主要是关注工作 面采空区“三带”的位置范围，因此，在对“三带”模拟结果影响不大的前提下， 模型进行了简化。上、下顺槽的长度都取20m，采空区走向取150m，其他模型尺寸完全和现场实际一致。根据潘北矿1121（3）工作面采空区实际情况，建立复杂的采空区三维物理模 型，如图5-1所示，模型尺寸如表5-2所示。表5-2物理模型尺寸模型名称长/m宽/m高/m备注工作面130.05.63.6长：指倾向长采空区150.0130.03.6

8、长：指走向长上顺槽20.03.02.2下顺槽20.04.43.05.4数学模型的建立一、连续性方程的建立根据质量守衡原理直接可以得到连续性方程为:（5.4.1） TOC o 1-5 h z du6udu八A + j + k = 0 dxdydz其中：u，u，u为单元体的平均流速，它和孔隙中的平均流速V.关系为：i j kiu = v. - n（5.4.2）N孔隙率。以下数学模型中，除特别指明外，所指速度都是指单元体的平均流速。二、动量方程的建立在惯性（非加速）坐标系中i方向上的动量守恒方程为：E（p u）+（p uu ）一祝 + 当 + F（5.4.3）dt f 1 dx. f 1 jdx.

9、6x.1式中：p静压，Pa；t .应力张量，Pa；且应力张量由下式给出：七+习-:啊8.5F包含了其它的模型相关源项，如多孔介质和自定义源项；采空区内是属于多孔介质流体区域，F.可由下式得出：F= D 四 +X C 1 p vv（5.4.5）i. j . 2 j jj=1j=1其中D和C是规定的矩阵。在多孔介质单元中，动量损失对于压力梯度有贡 献，压降和流体速度（或速度方阵）成比例。对于简单的均匀多孔介质：S = v + C p v v（5.4.6）i 以 i 2 2 j j其中a是渗透系数，C2为内部阻力因子，简单的指定D和C分别为对角阵和 C2，其它项为零。通过多孔介质的层流流动中，压降和

10、速度成比例，常数C2可 以考虑为零。5.5 FLUENT软件简介FLUENT是美国FLUENT公司推出的用于模拟具有复杂外形的流体流动以及热 传导的计算机软件。它提供了很好的网格灵活性，可以使用结构化网格和非结构 化网格，例如二维三角形或四边形网格、三维四面体/六面体/金字塔形网格来 解决具有复杂外形的流动，也可以用混合型非结构网格。它还允许用户根据解的 具体情况对网格进行修改（细化或粗化）。FLUENT软件采用有限容积法，可选择一阶迎风、二阶迎风、中心差分、幕指 数以及QUICK等格式。速度和压力耦合采用SIMPLE或SIMPLEC算法，代数方程 的求解可以米用多重网格及最小残差法。FLUE

11、NT是用C语言写的，具有很大的灵活性与能力。除此之外，为了执行高 效，易于交互控制，以及灵活的适应各种机器与操作系统，Fluent使用client / server结构，因此它允许同时在用户桌面工作站和强有力的服务器上分离地 运行程序。在FLUENT中，解的计算与显示可以通过交互界面、菜单界面来完成。用户界 面是用Scheme语言及LISP dialect写成的。高级用户可以通过写菜单宏及菜单函 数自定义及优化界面。FLUENT软件包括以下几个部分：FLUENT 解法器；prePDF，用于模拟PDF燃烧过程；GAMBIT，网格生成几何图形模拟以及网格生成的预处理程序；TGrid，从已有边界网格

12、中生成体网格的附加前处理程序；Filters(Translators)，转换其它程序生成的网格，用于FLUENT计算。 可以与ANSYS、I-DEAS、NASTRAN、PATRAN等程序进行数据交换。在FLUENT中可以用GAMBIT产生所需的几何结构以及网格，也可以在已知边界 网格(由GAMBIT或者第三方CAD/CAE软件产生的)中用TGrid产生三角网格，四面 体网格或者混合网格。也可以用其他软件产生FLUENT所需要的网格，比如 ANSYS(Swanson Analysis Systems.Inc.)、I-DEAS (SDRC)，或者MSC/ARIES， MSC/PATRAN以及MSC

13、/NASTRAN。大多数CAD/CAE软件都可以产生上述格式的网格。FLUENT软件解算器对于可压缩与不可压缩流动，稳态和非稳态流动，无粘流、 层流及湍流，牛顿流体及非牛顿流体，对流换热(包括自然对流和混合对流)， 导热与对流换热耦合，辐射换热，惯性坐标系和非惯性坐标系下的流动，多运动 坐标系下的流动，化学组分混合与反应，多孔介质流动，两相流，复杂表面形状 下的自由面流动等问题都可以进行较好得模拟。因此，FLUENT软件在航天、机械、 电子、汽车工业、材料处理、建筑设计和火灾研究以及暖通空调等领域都得到了 广泛的应用。5.6采空区“三带”模拟条件模拟采用隐式求解的求解方法。模拟需要的其它已知条

14、件如下：1、 工作面的推进速度取3.1m/d；X作面供风量为Q=30.0ms。工作面2、采空区倾向上的风阻变化系数在倾斜方向上，由于采空区的上、下巷道附近顶板受实体煤的支撑，而采空 区中部的上覆岩层无支撑，所以顶板岩石冒落的下沉量也不同。一般是采空区中 部覆岩的下沉量大于其上、下边界覆岩的下沉量。因此，在倾斜方向上，一般采 空区中部较上、下巷附近的孔隙率小，风阻大。令k(y)表示风阻变化系数，则 可得：5.6.1H1采空区中部覆岩的下沉量，m；H：上、下边界覆岩的下沉量，m；ly采空区倾向长，m；根据矿井实际条件及相关资料，H1/H2取1.30。3、采空区漏风风阻的计算（5.6.2）（5.6.

15、3）由采空区漏风风阻的计算公式为：&lR = axc 1slR bx o.5c 式中：a,b经验系数，取决于顶板冒落岩石性质，其值如表5-3所示;c冒落岩石的压实系数，它与工作面推进速度（m/d）有关，按式（5.6.4） 计算：c = 1.0e M .）（5.6.4）式中：x采空区内距工作面距离，m；l采空区滤流分支的长度，m；s采空区滤流分支的截面积，m2；Vf工作面推进速度，m/d；考虑采空区倾向的风阻变化，则可得采空区漏风风阻的计算公式为:（5.6.5）（5.6.6）laxcs k （y JR bx 0.5c一-2 s 2 k （y ）由采空区走向长为150m，工作面长度为130m，工作

16、面的推进速度取：3.1m/d。（5.6.7 ）（5.6.8）a，b系数，根据实际顶板岩性：a=0.04,b=40，则可以得到:R = 0.04 X1.209 k （）R 40 x0.6：（）表5-3不同岩性冒落岩石的a、b值冒落岩石种类ab松软粘土岩，页岩0.6 1.0101中硬粘土页岩0.2 0.571 100硬粘土页岩、砂岩0.06 0.15170砂岩、石灰岩0.03 0.0535 504、模拟基本假设1）、采空区内遗煤、冒落岩石与空气等混合物视为各向同性的均匀一致的多 孔介质；2）、上、下顺槽在同一水平面，不考虑位能的影响；3）、粘性阻力系数和内部阻力系数在z方向不发生变化。5.7采空区

17、“三带”模拟的结果及分析根据模拟计算结果，得到采场在z=1.0m及z=2.0m处切面的流场立体等值线 图，如图5-3所示，并得到z=1.0m处切面的二维流场等值线图，如图5-4所示。从图5-3中可以看出，两个切面的流场等值线分布基本上一致。这主要是因 为在Z方向上的粘性阻力系数和内部阻力系数不发生变化。从图5-4中可知，工作面风速明显大于采空区滤流速度。随着采空区的走向， 滤流速度下降的较快，这是由于随着采空区越来越被压实，阻力系数逐渐增大的 缘故。采空区的走向距工作面6.08.0m，倾向距下顺槽外帮30.0m左右处有一个 低速区，另外，距工作面上隅角11.0m处也出现了一个低速区，这主要是因

18、为采 空区阻力分布不均匀及壁面作用，改变了风流的前进方向，形成旋涡所致。1.7676 D1SI6 QD4J10。漩 0.DM2 ODMO 0.DQI7 Q.顺 皿6 DDCC3 QD0Q2 口硕I图5-3 z=1.0m、2.0m流场立体等值线图图5-4 z=1.0m处二维流场等值线图以滤流速度0.1m/min和0.24m/min为三带”分界线，得到在z=1.0m及 z=2.0m处切面的“三带”立体图，如图5-5所示，并得到z=1.0m处切面的二维 “三带”位置及分布，如图5-6所示。从图5-5中可知，z=1.0m及z=2.0m两切面的“三带”分布位置及大小、形 状基本一致，这也是因为粘性阻力系

19、数和内部阻力系数不随Z方向变化的缘故。 从图5-6中可知，在这种采场条件下，氧化带的最大的宽度为40.0m左右。且在 采空区倾向不同的位置，氧化带的位置及宽度不一样，进风侧宽度最小，倾向距 下顺槽外帮30m左右氧化带宽度最大。根据推进速度为3.1m/d，则氧化带存在的最长时间为12.9d。由经验得出采空区自然发火期为36个月，所以推得采空区内部自然发火的可能性较小。图5-6 z=1.0m处二维三带分布图图5-5 z=1.0m及z=2.0m处三带立体图根据模拟结果，得到z=1.0m处切面采场二维风流矢量分布图，如图5-7所 示。从图5-7可以看出，风量从下顺槽大部分流入工作面，有小部分风量漏入采

20、 空区，且漏风量的大部分是从工作面下隅角附近流入采空区的。II1I：；,图5-7 z=1.0m处二维风流矢量图从模拟结果中得知，从下顺槽外帮起，每隔10m处的工作面风量分布，如 表5-4所示，并得到工作面的风量分布曲线，如图5-8所示。最后根据表5-4及 图5-8,得到工作面每段长为10m的漏风量分布，漏风量为正时表示向采空区漏 入，负值表示由采空区漏出，其分布柱形图如图5-9所示。表5-4工作面风量分布,m3/s距离/m102030405060708090100110120风量/m/min169017501770175017501750175017501760177017901800-601

21、23450080012340070距离段/m图5-9图工作面漏漏量分布柱形图OOOOOOOOO208642005588777776612010080-40从表5-4及图5-8中可知，工作面风量分布是不均匀的。由于采空区较大部 分漏风量是从工作面下隅角流入的，因此工作面风量在下隅角附近是最小的。从图5-9可以看出，采空区的大部分漏风量是发生在010m范围内，而 1030m区间内，风流就开始返回至工作面内。3040m区间内，又有小部分风 量漏入采空区，说明进风侧3040m处也存在漏风源，采空区漏风源不止一处。 因此，为了减少采空区的漏风量，应重点防止采空区下隅角附近的漏风，同时， 也应注意工作面其

22、它地方的漏风。以上结果表明，漏风量及漏风风流在工作面倾向上，沿垂直于工作面煤壁中 心线并不是呈对称分布的，这种现象从采场风流矢量分布图及流场等值线图中也 可以看出。这主要是因为在倾斜方向上，采空区中部比上、下顺槽附近的孔隙率 更小，风阻更大的缘故。考虑到现场综采工作面的实际条件，无法准确测出工作面的风量分布，因此 无法同现场工作面实际风量分布曲线进行比较，但同论文采空区孔隙率非均质 模型及流场分布模型2中的工作面内漏风量分布曲线做了对比，漏风量分布趋 势较为吻合。5.8工作面供风量对采空区“三带”分布的影响为了考察工作面供风量对采空区“三带”位置及范围的影响，在原采场条件不变的基础上，除现场实

23、际的供风量Q2=1800m3/min夕卜，再提出工作面供风量Q1=1200m3/min和工作面供风量为Q3=2400m3/min两种方案进行模拟。模拟结果 分别如图5-105-15所示。图5-14 Q3=2400m3/min的三带立体图图5-15 Q3=2400m3/min二维三带分布图从图5-105-15中可知，工作面供风量不同，采空区内流场分布不同，从而 导致采空区“三带”的位置及范围也不同。随着工作面的供风量越大，采空区“三 带”位置越向采空区深部延伸，而且“三带”的范围也越大，由此说明，采空区 “三带”的位置及范围与工作面的供风量有直接的关系。另外，从图5-105-15中还可以看出，尽管工作面的供风量不同，但是采空 区走向距工作面6.08.0m,倾向距下顺槽外帮30.0m左右处的旋涡始终存在。虽 然此处的“三带”范围小，但是它并不随着工作面供风量的改变而改变，而是受 采空区走向及倾向的粘性阻力系数和内部