基于三维空隙率模型的采空区瓦斯运移及富集规律

2021-06-03曹明月李小芳孙浩石闫志铭

煤矿安全 2021年5期

徐超，曹明月，李小芳，孙浩石，闫志铭

（1.中国矿业大学（北京）共伴生能源精准开采北京市重点实验室，北京100083；2.中国矿业大学（北京）应急管理与安全工程学院，北京100083；3.山东科技大学矿山灾害预防控制重点实验室，山东青岛266590；4.阳泉煤业（集团）有限责任公司，山西阳泉045000）

煤与瓦斯突出等瓦斯事故严重威胁煤矿安全生产，而且矿井瓦斯事故多发生在综采工作面[1-2]。因此基于FLUENT模拟软件研究采空区内瓦斯运移及富集规律对预防矿井瓦斯事故有指导性作用。采空区是由垮落的岩石和遗煤构成的多孔介质，多孔介质的空隙率和渗透率分布对流体流动有较大影响。李树刚等通过研究垮落岩石的碎胀特性，提出了基于碎胀系数的采空区空隙率计算公式，并基于空隙率模型研究了采空区瓦斯渗流规律[3]；高建良等研究了在采空区渗透率均匀分布、分段均匀分布和非均匀连续分布下的采空区流体运移规律，其中非均匀连续分布最符合现实情况[4]；陈鹏等研究采空区及其上覆岩层的岩石碎胀性质，基于“O”型圈理论提出了更符合实际的采空区空隙率三维分布模型[5]。在参考前人研究结果的基础上，以平舒煤矿15111综采工作面为研究对象，基于“O”型圈和砌体梁理论建立采空区三维空隙率和渗透率模型，运用FLUENT软件，模拟在有无重力2种情况下的采空区瓦斯流动状态，此外还对的Blake-Kozeny公式涉及到的多孔介质平均粒子直径的取值进行了模拟，得到了在不同粒径下采空区瓦斯运移和富集规律，并对比模拟结果和实测数据，分析模拟结果的准确性。

1 采空区三维流场数学模型

煤层开采后采场中形成矩形采动空间，受二次应力影响，围岩发生垮落、断裂和变形，形成了采场上覆岩体结构的“砌体梁”模型[6]。基于“砌体梁”理论和“O”型圈理论，提出如下模型：取进风巷侧工作面与采空区的交点为坐标原点O，x轴方向为沿采空区走向方向，y轴方向为沿采空区倾向方向，z轴方向为沿采空区竖直方向。

1.1 采空区三维空隙率模型

沿采空区走向方向，破碎煤体碎胀系数引起多孔介质空隙率变化范围较大。随着距工作面越远，空隙率呈负指数关系递减[7]。因此基于以上关系，建立采空区空隙率沿x轴方向的分布函数。

由煤岩体碎胀系数定义可知采空区内破碎岩体的碎胀系数沿x轴的分布为：

式中：kp为破碎岩石的碎胀系数，无量纲；∑h为直接顶厚度，m；m为煤层开采厚度，m；W为采空区上覆岩层下沉量，m。

式中：kp1为直接顶破碎后岩块碎胀系数，无量纲；x为采空区内某一点距工作面的距离，m；l为基本顶破断后岩体长度，m。

根据破碎煤岩体空隙率与碎胀系数之间的关系，空隙率沿x轴方向的分布函数如下：

式中：n（x）为采空区底板处y=0时沿走方向方向空隙率，无量纲。

根据“O”型圈理论，可得沿工作面倾向方向上偏离y轴原点的空隙率变化系数变化函数[8]：

式中：n（y）为采空区地板处空隙率沿y轴方向变化系数，无量纲；y1为工作面倾向长度，m；y为采空区内某一点距进风巷的倾向距离，m，取值范围［0，y1］。

在xy平面上采空区空隙率变化函数为：

在竖直方向上，将关键层视为临界线，在关键层以下的层位，从垮落带到裂隙带，破碎岩体空隙率沿竖直方向呈现逐渐减小的分布特征；在关键层以上层位，由于关键层支撑作用，受采动影响较小，空隙率趋近于0。由于数值模拟模型主要位于底部层位，关键层以上层位对于数值模拟结果影响较小，因此，只关注关键层以下的区域[9]。设该区域内采空区空隙率在z方向上服从线性递减规律，得到空隙率在xyz空间上的分布函数为：

式中：a、b为待定系数，无量纲；z为采空区内某一点竖直高度，m。

式中：H为关键层高度，m。

将15111综采工作面相关参数（∑h=5.9 m、m=3.5 m、kp1=1.3、l=30 m、y1=186 m、H=50 m）代入式（7）得到采空区三维空隙率分布函数：

运用MATLAB软件，得到采空区底板z=0处空隙率变化曲面，采空区底板垮落岩石空隙率分布如图1。从图1可以看出，沿x轴方向，空隙率呈递减趋势，这是由于采空区深部垮落岩石逐渐被压实，从而导致空隙率减小；沿y轴方向，空隙率呈先减小后增大变化趋势，这是由于两道附近顶板受支护巷道支撑作用，使得两侧采空区的空隙率比采空区中部大。在采空区底板处，采空区垮落岩石的空隙率基本呈“滑梯”状。

图1 采空区底板垮落岩石空隙率分布Fig.1 Distribution of caving rock porosity in bottom plate of gob

1.2 采空区三维渗透率模型

根据Blake-Kozeny方程，得到采空区多孔介质渗透率k与空隙率n之间关系，即：

式中：k为多孔介质渗透率，10-15m2；Dm为多孔介质平均粒子直径，m；n为多孔介质空隙率，无量纲。

1.3 采空区非均质多孔介质动量损失模型

在FLUENT软件中模拟定义采空区多孔介质时，需要定义其黏性和惯性阻力系数作为模拟采空区对气体的流动阻力[10]。

式中：Si为采空区多孔介质的动量损失源，无量纲；μ为动力黏度，Pa·s；ρ为流体密度，kg/m3；Dij和Cij分别为黏性阻力和惯性阻力损失系数矩阵；vj（j=1，2，3）为流体微元在x、y、z方向上的速度分量，m/s。

在黏性阻力中黏性阻力系数C1，惯性阻力中内部惯性阻力系数为C2[11]。通过UDF导入FLUENT模拟软件中，并在其中调用。

1.4 采空区流场控制方程

根据质量及动量守恒方程，建立采空区流场控制方程。

连续性方程：

式中：ui为在i方向的速度，m/s；xi为i方向上的坐标值，m；t为时间，s；qm为质量源项，kg/（m3·s）。

动量守恒方程：

式中：uj为在j方向的速度，m/s；p为流体微团上的压力，Pa；τij为应力张量，无量纲；xj为j方向上的坐标值，m；fi为在单位质量流体微团上的i方向上体积力，N；Fi为流体微元上的体力，N[12]。

2 采空区瓦斯运移数值模拟及结果

2.1 数值模型建立

平舒煤矿15111综采工作面位于阳泉煤矿15#煤东翼一采区，所采的煤层为15#煤层，煤层厚度平均为3.8 m。该综采工作面走向长为1 130.7 m，倾斜长为186 m，预计回采期间瓦斯绝对涌出量为42.2 m3/min，采用“U”型通风系统，核定配风量为2 500 m3/min。15111综采工作面采用综合机械化一次采全高的采煤工艺。

由于实际采场的复杂性与不确定性，在建模时需要进行简化。根据平舒煤矿15111综采工作面的实际情况，工作面及进、回风巷的参数以专业人员现场实测为主，建立了三维采场物理模型，15111综采工作面三维采场物理模型参数表见表1。运用ICEM软件对15111综采工作面简化物理模型进行非结构化网格划分，并进行局部加密，划分成1 432 500个网格单元。

表1 15111综采工作面三维采场物理模型参数表Table 1 Parameter table of 3D stope physical model in 15111 fully mechanized face

工作面进风口设置为速度入口（inlet），进风风速设置为2.89 m/s；回风口设置为自由出口（outflow），工作面与采空区之间的面和垮落带与裂隙带之间的面，设为交界面（interior）[13]。在模拟过程中将采场流体视为不可压流体、采空区视为各向同性多孔介质并假定采空区瓦斯涌出源来自于采空区，并且瓦斯涌出量为定值。

此次模拟研究重力因素（有无重力）和Blake-Kozeny公式中采空区平均粒径取值（0.08、0.12、0.16 m）对采空区瓦斯运移及富集规律的影响[14]。

2.2 模拟结果

2.2.1 重力环境设置的影响

考虑重力环境设置，在模拟三维流场中很有必要，在竖直方向上，由于瓦斯密度为0.716 kg/m3，空气密度是瓦斯密度的1.81倍，瓦斯在采空区内流动时，会出现由于重力作用的上浮效应，使得采空区瓦斯流动规律有所差别，采空区瓦斯体积分数分布云图如图2。

图2 采空区瓦斯体积分数分布云图Fig.2 Cloud charts of gas volume fraction distribution

由图2可以看出，采空区深部和靠近回风巷侧瓦斯体积分数较大，可知工作面向采空区漏风主要区域为靠近进风巷侧工作面处。从图2（a）、图2（b）可以看出，由于垮落带瓦斯涌出量较大且在无重力操作环境下，瓦斯在垮落带深处积聚，造成采空区顶部瓦斯体积分数比采空区底部瓦斯体积分数小，说明在不考虑重力因素时采空区内瓦斯不会出现上浮效应。从图2（c）、图2（d）可知，上隅角瓦斯积聚效应明显，这是因为采空区内积存大量高体积分数瓦斯，漏入采空区的风流通过回风巷带出大量瓦斯气体，在重力环境下，瓦斯比空气密度低，会产生上浮效应和流体分层现象，使得瓦斯在上隅角处以下扎姿态流入回风巷，并在上隅角上部采空区形成“倒三角锥”区域。

在不考虑重力因素下采空区最大瓦斯体积分数为0.876 4，在考虑重力因素下采空区最大瓦斯体积分数为0.839 3，发现在无重力因素下采空区瓦斯体积分数较高，这是因为不考虑重力因素时风流静压变化较小，渗流速度较低，漏风量较小。通过以上分析，在不考虑重力时，瓦斯不会出现明显的上浮效应，这与实际采空区内瓦斯在其中分布存在很大差异，因此在模拟采空区瓦斯运移时考虑重力很有必要。

2.2.2 平均粒子直径的影响

根据FLUENT模拟结果，不同平均粒子直径下高（低）瓦斯体积分数占比如图3。

图3 不同平均粒子直径下高（低）瓦斯体积分数占比Fig.3 Proportions of high（low）gas volume fraction under different average particle diameters

由图3可知，当瓦斯体积分数为0～5%，平均粒径为0.16 m时占总网格数最高为38.81%，平均粒径为0.08 m时占比最低为37.73%，设置不同平均粒径值对采空区低体积分数瓦斯富集影响较小，总体上随着平均粒径值的增加，低体积分数瓦斯所占区域占比增大。当瓦斯体积分数为5%～15%，平均粒径为0.16 m时占总网格数最高为13.89%，平均粒径为0.08 m时占比最低为8.72%，当采空区瓦斯体积分数5%～15%区段内，有明火的情况下就会发生爆炸，当采空区粒径越大，也易发生瓦斯爆炸事故。在瓦斯体积分数为60%～100%之内，平均粒子直径为0.16 m时占比最低为0。由此可以看出，随着采空区平均粒径的增大，采空区高体积分数瓦斯富集区变小，低体积分数瓦斯富集区增大。采空区平均粒径增大使采空区整体渗透率增大，工作面向采空区漏风量也增大，风流从进风巷附近漏入采空区并从回风巷附近漏出，带出更多高体积分数瓦斯，从而使得采空区整体瓦斯体积分数下降。

不同平均粒子直径下采空区瓦斯体积分数分布云图如图4。不同平均粒子直径下z=5 m截面处瓦斯体积分数分布等值线如图5。

图4 不同平均粒子直径下采空区瓦斯体积分数分布云图Fig.4 Cloud charts of gas volume fraction distribution in goaf with different average particle diameters

图5 不同平均粒子直径下z=5 m截面处瓦斯体积分数分布等值线Fig.5 Contours of gas volume fraction distribution at z=5 m cross section at different mean particle diameters

由图4可知，不同平均粒径下采空区瓦斯分布差别较大。当平均粒径为0.08 m时，采空区瓦斯体积分数范围为0～1；平均粒径为0.12 m时，采空区瓦斯体积分数范围为0～0.839；平均粒径为0.16 m时，采空区瓦斯体积分数范围为0～0.569。

由图5可知，在z=5 m截面处，不同平均粒径下采空区整体瓦斯分布规律一致，在回风巷侧采空区深部形成高体积分数瓦斯富集区。为了更准确的看到z=5 m处的瓦斯运移情况，布置1条沿走向的观测线y=100 m，z=5 m，观测采空区瓦斯体积分数变化和采空区流体渗流速度变化。

不同平均粒子直径下观测线上采空区瓦斯体积分数及渗流速度变化如图6。

图6 不同平均粒子直径下观测线上采空区瓦斯体积分数及渗流速度变化Fig.6 Variation of seepage velocity and gas volume fraction at the observation line at different mean particle diameters in the area

从图6可知，越远离工作面，观测线上瓦斯体积分数越小。不同平均粒径下采空区瓦斯体积分数相差较大，当平均粒径为0.08 m时，观测线上最大瓦斯体积分数为1；当平均粒径为0.16 m时，观测线上最高瓦斯体积分数为0.463，随着平均粒径增大，观测线上采空区最大瓦斯体积分数减少了约1/2。由此可见，采空区平均粒径的取值对采空区瓦斯流动影响较大。采空区初始渗流速度较小，说明工作面漏入采空区风量只为整体风量的一小部分，不同平均粒子直径下速度变化不明显。在距离工作面0～3 m内采空区平均粒径越大，流体初始渗流速度越小，这是因为风流漏风位置主要位于靠近进风巷侧工作面内，平均粒径越大漏入采空区的风量越多，从而导致工作面中部风流越小，渗流速度越小。在距离工作面3～32.3 m内，采空区平均粒径越大，流体渗流速度越小，渗流速度衰减的越慢。在距离工作面32.3～200 m，随着风流向采空区深部流动，流体渗流速度逐渐减小，并在最深处达到0 m/s。

结合图5、图6可知，平均粒径的取值影响着采空区漏风情况，从而影响瓦斯运移及富集规律，采空区平均粒径越大，表明采空区内煤岩空隙越大，也易引起采空区瓦斯事故。因此实际工程中，使用填充法填充采空区，减少采空区煤岩裂隙，从而减少采空区漏风，以达到防止瓦斯事故的目的。

3 现场实测对比

为了验证模拟得到结果是否与实际相吻合，此次研究对比了15111工作面实际工作情况下进回风巷风量和模拟得到的进回风巷风量。不同平均粒径下进回风巷风量分布见表2。

由表2可知，不同平均粒子直径下进回风巷风量差别不大，整体规律都呈进风巷比回风巷风量大的趋势，这是由于风流有一部分漏入采空区所致，随着平均粒径的增大回风量越来越小，说明漏入采空区的风量变多。对比实测结果和模拟结果，可知当平均粒子直径为0.12 m时误差最小，为3.981 3%，说明当设定初始条件平均粒子直径为0.12 m时，模拟结果最符合实际情况，也说明了此次模拟较为准确。