采动影响下工作面覆岩渗透率变化研究

孙祺钰，祝 彦，赵晓夏

(1.国家能源集团神东煤炭集团有限责任公司，陕西 神木 719315； 2.中煤科工集团重庆研究院有限公司，重庆 400039)

工作面上覆岩层具有不易接触和隐蔽的特点，在采动影响下地表的沉陷和覆岩的移动、变形会经历非常复杂的过程，针对这一复杂过程的理论研究和预计预测一直是学者们高度关注的重点。刘宝琛等[1]提出的随机介质理论逐渐为大家所接受，进而发展成为我国常用的概率积分法，为地表移动预计提供了重要支撑；煤炭科学研究院北京开采研究所[2]针对不同倾角煤层的岩层移动和地表沉陷进行了深入研究，将煤层分为水平煤层、缓倾斜煤层、急倾斜煤层，在研究过程中提出了导水裂隙的概念，并得出冒落带与导水裂隙带高度的计算公式；钱鸣高等[3-4]创新性地提出了关于上覆岩层下沉的“关键层”理论，该理论是岩层移动、顶板周期性矿压、地表沉陷等多方面理论的有机结合，综合考虑了上覆岩层岩性对岩层移动的影响，突破了前人开采沉陷理论中的局限性；谢和平等[5]利用相似材料模拟煤层上覆岩层节理的分布特征，研究表明节理的存在使得煤层覆岩的破坏加剧并控制着采动覆岩裂隙的发育状况；郭文兵等[6]对厚煤层综采放顶煤开采引起的地表沉陷情况进行了实际观测，并分析了关于厚煤层综采放顶煤开采条件下的地表下沉规律。这些理论及研究成果对我国上覆岩层移动变形的预测起到了非常重要的作用。

在采动影响下采煤工作面上覆岩层的应力、应变的变化和煤岩层渗透率的演化规律，对工作面瓦斯的抽采和治理有着重要的指导意义，国内外学者针对这一课题做了大量的研究并取得了许多研究成果。彭信山[7]系统分析了急倾斜近距离下保护层开采岩层移动变形、裂隙动态演化及应力变化规律，建了保护层开采煤岩变形与卸压瓦斯流动气固耦合模型；陆卫东等[8]建立了岩体渗透率的计算模型，确立了“面应变—孔隙率—渗透率”模型，用以表征煤岩体渗透率的变化，并对其可靠性进行了验证；马建宏[9]基于煤的孔—裂隙双重介质物理模型，建立了煤层有效应力、吸附膨胀效应和外载应力共同控制下的煤体渗透率模型，得到了保护层开采过程中被保护层瓦斯渗透率随煤层应力变化的规律；周银波等[10]研究发现在下伏被保护层双重采动影响下，采场覆岩在短时间内会发生大范围运动，周期性地形成二次采动影响区，采场内原本富集稳定的瓦斯涌入二次采动影响区，造成该区域在短时间内瓦斯富集并发生周期性变化。

1 采煤工作面上覆岩层开采沉陷预计模型

目前，国内外相关学者针对地表移动与变形的预计进行了大量的研究，并获得了较多的成果。地表沉陷预计方法主要有[11-13]：基于现场观测经验的预计方法(包括图表法、剖面函数法等)、基于影响函数的预计方法(以Knothe影响函数为典型代表)和理论模型法(包括随机介质理论、基于弹性理论或塑性理论的有限单元法、离散元法和边界元法等)。由于影响函数方法所采用的计算参数都是通过现场实测得到的，并且该方法对于任意的工作面开采程度、开采形状、地表任意点的位移的预计都具有较好的适用性，已成为国内外较为成功并被广泛应用的矿山地表沉陷预计的主流方法之一。

影响函数方法是对地表沉陷进行预计的一种有效的方法[14-17]、一种由经验方法向理论模型过渡的方法，其理论基础是分布函数，故称之为影响函数方法。该方法可描述为：在地下矿层之中每开采任意一个单元点，都会导致其周围一定范围内的地表产生沉陷，且越远离该单元点，地表沉陷量越小。地表某一点的最终沉陷则是矿层中所有开采单元点对该地表点所造成影响的叠加结果，地表沉陷示意图如图1所示。该方法的关键在于选择合适的影响函数，用以准确地描述地表沉陷及对所选择的影响函数在地下开采区域内进行正确积分。

图1 由矿层中某一单元点的开采所引起的地表沉陷示意图

在国内外相关学者提出的众多不同的影响函数方法中[18-21]，波兰学者Knothe推导的Knothe函数是较为成功并被广泛应用的一种影响函数，该影响函数能全面地揭示地下煤层开采单元点对其周围一定范围内地表的沉陷所造成的影响。Knothe影响函数在二维情形下的表达式如下：

(1)

Smax=ma

(2)

式中：Smax为采矿地质条件下地表所能达到的最大下沉值；R为采动影响半径；x′为地表预计点与地下开采单元点之间的水平距离；m为煤层开采厚度；a为下沉系数。

基于影响函数法的基本原理，极易将二维情形下的沉陷预计公式推广到一般三维情形下。三维情形下的地表下沉预计的影响函数表达式如下：

(3)

由二维情形下的地表水平移动影响函数表达式推广得到三维情形下的沿x、y轴方向的地表水平移动影响函数表达式，见式(4)～(5)：

(4)

(5)

式中:y′为地表预计点与地下开采单元点之间的垂直距离；h为采深。

显然，二维情形下的影响函数积分计算方法在三维情形下依然适用。采用三维情形下影响函数计算模型及其积分方法对矩形采空区的地表沉陷进行预计。引入xOy全局坐标系，其原点位于采空区左下角，x轴位于工作面走向方向，y轴位于工作面倾向方向。矩形采空区的长度和宽度分别为L和W。全局坐标系与局部坐标系如图2所示。

图2 全局坐标系与局部坐标系示意图

在预计点(x,y)处的最终下沉值S(x,y)是三维情形下的地表下沉影响函数(见式(3))在图2中积分区域A内的积分值。由于拐点偏移距的存在，积分区域并非实际采空区的边界而是其向内收缩拐点偏移距为d后的区域。为了简化数学模型，引入了一个局部坐标系x′O′y′，该坐标系的原点位于地表预计点处。因此，三维情形下地表最终下沉数学表达式为：

(6)

式中：d1、d2、d3、d4为x、y轴方向上的拐点偏移距；A为积分面积。

同理，在积分区域A内分别对式(4)和式(5)进行积分，可得到三维情形下的沿x、y轴方向的地表最终水平移动预计的数学表达式：

Ux(x,y)=∬fux(x′,y′)dA

(7)

Uy(x,y)=∬fuy(x′,y′)dA

(8)

式中：fux为沿x轴方向的地表水平移动影响函数；fuy为沿y轴方向的地表水平移动影响函数；Ux为沿x轴方向地表最终水平移动距离；Uy为沿y轴方向的地表最终水平移动预计距离。

根据覆岩沉陷预计的基本原理及覆岩最终沉陷预计参数a(h)、R(h)和d(h)，可得到三维情形下的覆岩最终移动预计模型，见式(9)～(11)：

(9)

(10)

(11)

式中：S(x,y,h)为覆岩在(x,y,h)处最终下沉距离；a(h)为覆岩下沉系数；R(h)为采动覆岩主要影响半径；d(h)为覆岩拐点偏移距；x为覆岩预计点在全局坐标系中的x轴坐标；y为覆岩预计点在全局坐标系中的y轴坐标；x′为局部坐标系中煤层开采单元体与覆岩预计点间在x′轴方向上的距离；y′为局部坐标系中煤层开采单元体与覆岩预计点间在y′轴方向上的距离。

垂直应变是覆岩下沉值对-h的一阶导数，即:

(12)

x轴方向和y轴方向上的水平应变是其相应方向上水平移动的一阶导数，即:

(13)

(14)

三维情形下，采动覆岩的总应变可看作是体积应变，即覆岩在变形过程中，单位体积覆岩体的体积改变。根据体积应变的计算公式，可得三维情形下覆岩的总应变数学表达式：

εt(x,y,h)=εx(x,y,h)+εy(x,y,h)+

εz(x,y,h)

(15)

式(9)～(15)即为三维情形下采空区上覆岩层的静态沉陷预计模型。

覆岩下沉系数ai是在沉陷预测中最重要的参数，其与岩层的厚度和坚硬岩石占比有关，其经验计算公式如下：

ai=1.003 2(h/n)-0.009e0.000 05(35-ηi)

(16)

式中：n为岩层总数；ηi为第i层中坚硬岩石所占百分比；i为岩层层数。

要确定某一层的主要影响半径，就把这一层看作h/n的悬梁臂。左侧悬梁通过弹性固定端受到垂直限制，右上悬梁受到前一层Smax(zi-1)上表面的最大可能沉陷的限制。覆岩荷载大小qi也是影响岩层主要影响半径Ri的重要因素。Ri必须进行修正，以适合在地下岩层中应用，修正后的公式如下：

(17)

(18)

式中:Ri为岩层主要影响半径；Smax(zi-1)为第i-1 岩层最大下沉量；E为岩石地层的弹性模量，可以估算成软岩层弹性模量的0.49倍；Qi为第i个岩层的岩层因子；qi为第i个岩层的上覆荷载；γ为覆盖层密度。

(19)

第i层的偏移量可由预测地面沉陷的经验公式计算得到，见式(20)：

(20)

式中di为第i层的偏移量。

2 试验工作面覆岩渗透率变化规律研究

2.1 试验工作面背景

阳泉五矿位于山西省沁水煤田的东北部，井田内含煤10余层，煤层标高为456～535 m，煤层厚度为6 m左右，煤层倾角为3°～8°。8136试验工作面平均标高为517 m，工作面走向长度为1 000 m，倾向长度为280 m，采高为6.0 m，煤层平均倾角为5°，近似为水平煤层。

2.2 覆岩层移动计算结果及分析

采用三维情形下的工作面上覆岩层静态沉陷预计模型并通过Matlab软件对阳泉五矿8136试验工作面覆岩的移动与变形情况进行预计、分析。以三维情形下煤层上方100 m范围内的覆岩总应变分布规律为例进行详细分析。

煤层上方40、50、60、70、80、90 m高度处覆岩总应变等值线分布如图3所示，不难发现在煤层上方100 m高度范围内覆岩存在高总应变值。

(a)斜45°视角

(b)倾向视角

(c)走向视角图3 不同观察视角下煤层上方40、50、60、70、80、90 m 高度处覆岩总应变等值线分布云图

由图3可知，在工作面倾向距离煤层上方40～45 m的高度处，拉伸总应变区的宽度为12.5 m，最大拉伸应变大于0.200;压缩总应变区的宽度为10.0 m，略小于拉伸总应变区的宽度，最大压缩应变为-0.024。距离煤层上方90～95 m的高度处，拉伸总应变区的宽度为20.0 m，最大拉伸应变为0.066；压缩总应变区的宽度为35.0 m，最大压缩应变为 -0.010。

在工作面走向距离煤层上方40～45 m的高度处，拉伸总应变区的宽度为7.0 m，最大拉伸应变为0.020；压缩总应变区的宽度为7.0 m，最大压缩应变为-0.065。距离煤层上方90～95 m高度处，拉伸总应变区宽度为20.0 m，最大拉伸应变为0.066；压缩总应变区的宽度为40.0 m，最大压缩应变为-0.025。

由图3可知，覆岩总应变区包括3个区域：①拉应变圈(由下至上，倾向上宽度和走向上宽度均逐渐增大，但倾向宽度略大于走向宽度。煤层上方90 m高度处的拉应变圈在倾向上可覆盖至进、回风巷内侧0～32 m附近，在走向上可覆盖至开切眼、终采线内侧0～30 m附近；另外，拉应变圈在采空区4个边角处为弱拉应变区，高拉应变区的倾向长度和走向长度由下至上均逐渐减小)；②压应变圈(由下至上，走向上宽度和倾向上宽度也均逐渐增大，但倾向上宽度略小于走向上宽度。煤层上方90 m高度处压应变圈在倾向上可覆盖至进、回风巷内侧 32～50 m附近，在走向上可覆盖至开切眼、终采线内侧30～53 m附近)；③中部卸压区(位于采空区中部，且由下至上基本呈逐渐变窄趋势，处于较弱的拉伸状态)。

总应变能够综合反映工作面覆岩在采动影响下的移动与变形情况，可指导判断覆岩中裂缝、断裂或离层的存在位置和发育范围。另外，总应变值的正负反映采动影响下岩石体的膨胀与收缩，正值代表岩石内孔隙的强度，表明覆岩内透气性的增加。根据采空区覆岩总应变的分布规律，将覆岩裂隙划分为3个区域：①环形拉伸裂隙发育区(煤岩体处于拉应变圈内，其充分膨胀变形，采动裂隙发育，透气性增高，是瓦斯主要的运移与富集区)；②环形压缩裂隙闭合区(煤岩体处于压应变圈内，裂隙受到挤压而收缩、封闭，透气性差，瓦斯运移困难)；③中部压实区(煤岩体处于卸压区内，采动裂隙在上覆岩层压力的作用下基本被压实，拉应变值较小甚至趋近于0，透气性较差，瓦斯运移较困难)。因此，14#、13#、12#、11#煤层及赋存瓦斯的K3和K4灰岩围岩(高度为43～45 m、46～48 m、55～58 m)在采空区四周均存在着高环形拉伸裂隙发育区，采动裂隙发育程度高，瓦斯解吸并向邻近裂隙发育的围岩运移。

2.3 覆岩层渗透率变化计算结果及分析

由采动渗透率变化可知，在煤岩体原始孔隙率已知的情况下，采动渗透率变化仅为总应变的函数[10]。通过调研相关文献和经验对比，确定了阳泉五矿8136试验工作面上覆岩层的原始孔隙率。利用前述数学模型，对煤层上方100 m高度范围内覆岩渗透率变化分布规律与瓦斯渗透运移路径进行分析。根据覆岩岩性和厚度等，合理地将覆岩进行了分组处理，并求得各分组平均孔隙率，如表1和表2 所示。

表1 8136试验工作面覆岩孔隙率

表2 8136试验工作面覆岩各分组平均孔隙率

煤层上方40、50、60、70、80、90 m高度处覆岩渗透率变化等值线分布云图如图4所示。

图4 煤层上方40、50、60、70、80、90 m高度处 覆岩渗透率变化等值线分布云图

由图4可知，采动覆岩渗透率变化在采空区四周呈环形分布，并且进、回风巷内侧的渗透率变化要明显大于开切眼、终采线内侧的渗透率变化。根据等值线的分布特征，可以将采动覆岩在采空区四周的渗透率变化分为3个区域，即环形渗透率增高区、环形渗透率降低区和渗透率恢复区，并且这3个区域与总应变的拉应变圈、压应变圈和中部卸压区的位置、范围基本相同。距离煤层上方50 m处，渗透率增加到原始渗透率的13倍；距离煤层上方60 m处，渗透率增加到原始渗透率的37倍多；此后随着距离煤层高度的增加，渗透率逐渐减小，距离煤层上方高度90 m处，渗透率仅增加到原始渗透率的7倍。

采动影响下工作面倾向、走向主断面覆岩渗透率变化局部等值线云图如图5～6所示。

图5 采动影响下工作面倾向主断面覆岩渗透率变化局部 等值线图(y=-3～45 m，h=30～170 m)

图6 采动影响下工作面走向主断面覆岩渗透率变化局部 等值线图(x=-10～45 m，h=30～170 m)

由图5～6进一步对煤层上方60 m高度附近覆岩采动渗透率变化进行对比分析，可以确定覆岩高渗透率变化区位于煤层上方59～67 m内，在该范围内采动渗透率比原始渗透率可高出22～37倍。在h=59 m高度处覆岩的采动渗透率变化最大，是原始渗透率的37倍多，推断是由于此高度处总应变值比其他高度处总应变值大所致。这一区域的瓦斯运移活动最为活跃，采动条件下高浓度瓦斯将在此大量聚积，会形成高位环形瓦斯富集区。

3 结论

1)以影响函数法为基础，开发了覆岩静态三维沉陷预计模型，提出了三维总应变的概念及其数学表达式。推导得出总应变与孔隙率、渗透率变化之间的表征关系式，在总应变基础上对采动覆岩孔隙率和渗透率变化进行了预计，分析表明二者的分布情况基本与覆岩总应变分布情况具有一致性。

2)阳泉五矿8136试验工作面在最终沉陷状态下，进、回风巷内侧的渗透率变化大于开切眼、终采线内侧的渗透率变化，并根据采动覆岩渗透率变化在采空区四周呈环形分布的特征将其划分为3个区域：环形渗透率增高区、环形渗透率降低区、渗透率恢复区。

3)分析得出8136试验工作面最终沉陷状态下高位环形瓦斯富集区的具体位置和范围，以及其采动过程中高位瓦斯富集区的动态演化规律。