0 引 言
我国煤田水文地质条件十分复杂，煤炭开采受到多种形式的水害威胁[1-4]，从以往突水资料统计来看，80%的突水事故与断层有关[5-7]。
目前，许多专家、学者对断层突水问题进行了研究，代表性的如黎良杰等[8]把断层分为张开型与闭合型，分析得出张开型断层的突水机理是断层两盘在承压水作用下产生了张开，承压水沿张开裂隙突出，闭合型断层的突水机理主要是断层两盘按板的规律破坏或断层两盘关键层接触部产生强度失稳的结论；卜万奎等[9]针对贯穿煤层顶底板的大型断层采动活化及突水的力学特征，建立了采场底板断层活化机理分析的力学模型；浦海等[10]通过数值模拟的手段，模拟了在有无断层的2种情况下松散含水层下煤层开采围岩及断层的应力、塑性区等分布规律，掲示了含断层构造顶板突水通道形成的机制；VICTOR Bense等[11]通过数字图像分析的方法，研究了断层带周围物质导水性的空间分布、变形机制及水力学特性；武强等[12]阐述了断层滞后突水的弱化机理和主要控制因素，提出了煤层底板断裂构造突水时间“弱化效应”的概念；以上研究成果对于解释断层突水这一灾害现象有着重要意义，但目前关于断层突水问题的研究大都集中在断层作为导水通道沟通顶底板含水层的领域，对于断层既作为导水通道又作为唯一突水水源的突水类型研究较少，因地层中没有充水含水层的存在，这种类型的断层突水便具有极强的不可预见性和突发性，危害极大，同时由于水文地质条件具有的特殊性，断层作为突水水源的突水机制与其他类型的断层突水有着明显差异。
针对上述问题，笔者以弘财煤矿1301回风巷发生的突水灾害为案例，通过数值模拟方法再现了掘进巷道断层突水的动态过程，分析了突水过程中围岩应力场、渗流场和位移场的多场灾变耦合演化规律，提炼了有效的突水前兆信息特征，相关研究可以为富水型断层突水现象的解释及监测预警提供科学依据。
1 工程概况
弘财煤矿矿井设计生产能力为30万t/a，位于贵州省水城县勺米乡，距六盘水市西南约24 km，总体地形南西低、北东高，相对高差399 m，以中低山为主。区内地层由新到老分别为第四系、下三叠统飞仙关组、上二叠统龙潭组、峨眉山玄武岩组，主采煤层位于上二叠统龙潭组，煤厚2～4 m，倾角8°～10°。研究区地质构造位于扬子准台地威宁北西向构造变形区格目底向斜西翼，整体为单斜构造，倾向北东。研究区地表属于珠江水系巴拉河流域，地表水系不发育，区内主要依靠大气降水补给，受到地势影响，一般就近补给，就近排泄，区内碎屑岩构造裂隙水均以大气降水为主要补给来源，地下水动态随季节变化明显,一般每年5月地下水流量、水位开始回升，6—9月为最高值，期间出现1～3次峰值，10—12月进入平水期，水位、流量开始逐渐递减，4月降为最低值。
2016年7月弘财煤矿1301回风巷掘进至438 m时发生突水灾害，总涌水量约4 200 m3，持续时间4 h，灾害发生后经过勘查，发现1301回风巷掘进工作面揭露的一条强富水正断层是造成突水事故的主要原因，断层倾向216°，倾角69°，断距10 m。断层所切割地层岩性主要为黏土岩、细砂岩、泥质粉砂岩及煤层，断层两盘岩体裂隙发育程度和含水性均较差，大气降雨在断层上覆地表洼地汇集后沿断层破碎带长期入渗造成断层带强富水，在本次突水事故中断层本身即是导水通道也是储水空间。
2 富水型断层突水流固耦合分析
2.1 流固耦合数学模型
假设岩体视为多孔介质，流体在岩体中的流动满足 Darcy定律，同时满足Biot流固耦合方程，即
(1)
式中：λ、G为Lame常数；p为孔隙水压；εv为体应变； xj为j向的坐标；uj为j向的位移；fj为j向的体积力；K为渗透系数；
为拉普拉斯算子；S为弹性释水系数；∂p/∂xj反映了渗流对岩体骨架的影响，∂εv/∂t则反映了岩体的应变对渗流的影响。
式(1)中并没有包含岩体变形对渗透系数及孔隙率的描述，在以往的研究中,为了简化模型，普遍把这些参数看做是常数,但岩石在受力状态时，其渗透系数和孔隙率会随着应力-应变的变化而变化[13]，为了能够更准确地进行突水模拟，笔者以文献[14-15]提供的渗透性和孔隙率控制方程为依据，通过编写程序把式(2)嵌入计算模型中，从而可以实现巷道掘进过程中围岩渗透系数和孔隙率随岩体变形同步变化。即
(2)
式中：φ为孔隙率；φ0为初始孔隙率；φr为孔隙率极限值；
为平均有效应力(拉应力为正)；α1为应力影响系数，其值可以取5.0×10-8Pa；α为Biot有效应力系数；σ1、σ2、σ3分别为3个主应力；K0为岩土体的初始渗透系数；Δε为体应变增量。
2.2 数值计算模型
以弘财煤矿1301回风巷的掘进工作面的实际地质条件为背景，建立Mohr-Coulomb本构模型(图1)，模型地层倾角10°，断层倾角60°，断层带宽度10 m，掘进工作面位于煤层上覆砂岩地层，分19次推进，前11次每次推进4 m，后8次靠近断层破碎带，每次推进2 m，在巷道顶底板及断层周围设置有位移监测点和最大剪应变监测点。
模型力学边界条件为底部采用全约束边界条件，模型顶部采用自由边界条件，模型左、右采用x方向固定边界。渗流边界条件为断层破碎带左、右及下部为隔水边界，断层破碎带顶部根据突水发生前洼地实测水位高度，设置4 m给定水头，断层破碎带内部根据p=ρgh估算初始静水压力 (p即孔隙水压，ρ为液体密度，g为重力加速度，取9.8 m2/kg，h 为断层破碎带内任一位置到洼地积水面的高度)，断层破碎带渗透率初始值为1.2×10-10 m2，模型各岩层力学参数取值见表1。
图1 数值计算模型
Fig.1 Numerical calculation model
表1 岩体力学计算参数
Table 1 Mechanical parameters of rock mass
岩层体积模量/MPa剪切模量/MPa泊松比内摩擦角/(°)黏聚力/kPa抗拉强度/kPa密度/(kg·m-3)砂质黏土岩4 7303 5500.20479001 1002 650泥质粉砂岩2 4001 7300.21455006402 510粉砂岩4 9504 3500.16501 0009302 600煤层1 6707700.30324504202 400细砂岩4 8004 2200.17489801 2002 500断层1 5007800.28305005502 200
2.3 数值计算结果与分析
2.3.1 围岩应力场
掘进工作面至断层不同距离时围岩应力场的变化模拟结果如图2所示，由于掘进扰动引起应力的重分布，巷道两端出现应力集中区范围不断增大，断层破碎带内应力比周围围岩低，破碎带内为相对卸压区，且靠近断层端掘进工作面的应力集中程度小于巷道另一端。随着工作面的推进，掘进面附近围岩应力出现不断变化，从而引起拉、剪破坏，这与巷道顶板监测点最大剪应变的变化是相对应的(图3)，即最大剪应变先增大后减小，后在体系调整下已发生变形破坏的塑性围岩最大剪应变再次发生增大。隔水层的应力变化特征基本可概化为：掘进距离L≤24 m的缓慢增大区(≤20 000时步)，掘进距离24 m<L≤52 m的快速增大区(20 000-30 000时步)，掘进距离L>52 m(>30 000时步)的迅速减小区。
图2 不同掘进距离时围岩应力场分布
Fig.2 Distribution of stress field of surrounding rock at different excavation distances
图3 巷道顶板监测点最大剪应变变化
Fig.3 Maximum shear strain variation curves at roadway roof monitoring points
2.3.2 围岩位移场
工作面掘进过程中围岩垂向位移的变化模拟结果如图4所示，由图4可知，由于掘进扰动造成位移场的重新分布，巷道顶板出现垂向负位移集中，底板出现正位移集中，且随着工作面推进，不断增大；当掘进工作面远离断层时，断层破碎带内垂向负位移比周围围岩低，带内为相对减小区，而工作面越接近断层破碎带，上述减小现象越不明显。
各监测点变化曲线如图5和图6所示，随着工作面的推进，掘进工作面顶板及断层左盘的垂向位移总体趋势呈现增大，位移变化特征基本可概化为掘进距离L≤24 m的平缓增大区域，掘进距离24 m<L≤52 m的快速化增大区域，掘进距离L>52 m的急剧增大区域。而断层右盘的各个监测点，由于与掘进工作面相隔断层破碎带，因而其受到的影响有一定减小。
图4 不同掘进距离时围岩位移场分布
Fig.4 Distribution of displacement field of surrounding rock at different driving distances
2.3.3 围岩塑性区
图7是掘进工作面距离断层不同距离时围岩塑性区的变化模拟结果，由图可知：随着工作面的推进，靠近断层一端的围岩塑性区范围不断增大，当掘进距离L=48 m时掘进工作面与断层破碎带过渡区域产生了孤立的损伤破坏区；当掘进距离L=52 m时掘进工作面与断层破碎带过渡区域的孤立损伤破坏区进一步扩大；当掘进距离L=56 m时断层破碎带与掘进工作面之间的损伤破坏区已相连，此时极有可能发生突水事故。
2.3.4 断层带水压
图8是掘进工作面距离断层不同距离时断层带内水压的变化模拟结果，从图8可知，当掘进距离L=24 m时，掘进扰动尚未对断层产生明显影响，当掘进距离L=48 m时，富水断层下部的水压显著增加，当掘进距离增加到L=52 m时，水压继续缓慢增加，当距离增加到L=56 m时，断层带内的水压迅速下降，此时，断层与工作面之间的受损区相连， 断层内的水压变化规律与隔水层的应力变化规律一致。
图5 巷道顶板监测点位移变化
Fig.5 Vertical displacement variation curves of roof monitoring point ofroadway
图6 断层两盘监测点位移变化
Fig.6 Vertical displacement variation curves of two walls of fault monitoring points
图7 不同掘进距离时围岩塑性区分布
Fig.7 Distribution of plastic zone of surrounding rock at different driving distances
图8 不同掘进距离时水压分布
Fig.8 Distribution of water pressure at different excavation distances
3 富水型断层突水过程流场分析
3.1 非线性渗流耦合模型
目前的研究普遍认为破碎岩石的渗流一般不服从Darcy定律而服从Forchheimer非线性渗流方程[16-18]，研究区断层破碎带胶结程度差，颗粒粒径和孔隙均较大，突水通道形成后，断层破碎带流速明显增大，渗流非线性特征明显，适用于Forchheimer方程描述，而突水发生后，地下水在巷道内的流动为紊流，则满足黏性牛顿流体流动的Navier-Stokes 方程。本次研究把地下水在这2个区域内的流动联合在一起，应用多物理场耦合模拟软件Comsol Multiphysics进行强富水断层突水过程中的渗流场的径流规律分析。数值模型中流体密度为1 000 kg/m3，动力黏度为1×10-3 N·s/m2，断层破碎带顶部为定水头边界，与巷道交界处附近水压为0.52 MPa，其余边界为无流动边界，模型其他参数见表1。地下水在不同阶段的动量方程如下：
1)断层破碎带非达西Forchheimer流方程式为
(3)
2)巷道Navier-Stokes紊流方程式为
(4)
式中：u为流体流速；η为动力黏度系数；F为流体阻力(与重力和流体可压缩性有关)；I为单位矩阵；Q为源(汇)强度；φ为孔隙度；βF为Forchherimer曳力参数,取
且
3.2 数值计算与分析
地下水从断层破碎带的非达西Forchheimer流到巷道Navier-Stokes 紊流过程中的压力分布如图9所示，压力在断层破碎带与突水通道交界面附近产生低值，主要是由于水流在压力差下涌入巷道，把压力转化为动能，造成此处压力急剧下降。
图9 水流路径上压力分布
Fig.9 Distribution of water pressure on flow path
图1中所示各关键监测点的流速和压力变化曲线图如图10所示。由图10a可知流动过程中速度总体趋于增大，自断层破碎带顶部到断层破碎带上部小幅增大，速度大幅增加于突水通道内和巷道内。在Forchherimer流区域，断层顶部到断层上部流速出现小幅增加，从4.2 mm/s增大到5.5 mm/s，但在断层与突水通道交界处流速出现明显减小并伴随波动，从5.5 mm/s下降到4.0 mm/s，原因是在该处流向发生改变，部分地下水的动能消耗于克服水分子之间的摩擦。地下水进入突水通道后流速又急剧增大，从4.0 mm/s增大到8.3 mm/s，在突水通道出口处再次因流向变化而减小，从8.3 mm/s下降到5.0 mm/s，进入巷道后，地下水呈现紊流，流速增大至1.0 mm/s，并保持稳定。从图10b可以看到沿流程压力不断减小，从断层顶部到突水通道位置压力迅速减小，从0.52 MPa下降到0.16 MPa，当水流进入突水通道后，压力的减小速度趋缓，从0.15 MPa下降到0.10 MPa，当水流进入巷道后，压力保持稳定。
(a)各监测点流速
(b)压力曲线
a—断层顶部；b—断层上部；c—突水通道入口；d—突水通道内；e—突水通道出口；f—巷道内
图10 各监测点的流速与压力
Fig.10 Flow velocty Water pressure at each monitoring point
4 结 论
1)以弘财煤矿巷道突水为案例，利用数值模拟方法再现了巷道掘进过程中富水型断层突水从发育、发展到发生的过程中围岩应力场、位移场和渗流场的变化规律，计算表明，富水型断层突水是在开挖扰动和断层破碎带内水压共同作用下的结果。
2)开挖扰动初期，巷道及其附近围岩发生损伤，隔水层内部开始萌生裂隙，围岩的应力、位移以及渗流等多场信息均平稳增长；当掘进工作面接近富水断层时，损伤破坏区进一步扩大，隔水层及其附近围岩内部裂隙开始局部扩展、水压跟踪传递，隔水层和断层带的多场信息数据进入快速增长阶段；随着工作面的持续推进，形成大范围的损伤破坏区，隔水层的裂隙通道逐渐贯通，当工作面推进至临界值时，隔水层围岩的多场信息发生突跳，围岩瞬间失稳突水。整体上来看隔水层岩体内的多场信息在突水过程中可以分为3个阶段，即采掘工作面远离断层时的初始变化阶段，接近断层时的快速变化阶段和突变阶段。
3)突水发生过程中，地下水沿断层破碎带自上向下涌入巷道，径流过程中压力持续减小，流速总体趋势增大，但在进入突水通道和巷道过程中流向均发生改变，此时流速明显降低，随后再次增大，并在进入巷道后流速和压力逐渐趋于稳定。
4)通过对多场信息变化规律的内在联系及其对富水型断层突水预测的有效性分析，发现多场变化信息在富水型断层突水前具有相对明显的前兆特征：在隔水层破断诱发突水前，围岩应力持续增长后发生回跌；围岩位移平缓增长后急剧增大，同时出现断层破碎带内垂向负位移减小的现象；在接近断层后断层带内水压明显增大，但在濒临突水前出现减小。以上信息可以作为富水型断层突水灾害监测与预警的前兆信息源，为矿井突水灾害的防治提供科学依据。
[1] 武 强.我国矿井水防控与资源化利用的研究进展、问题和展望[J].煤炭学报，2014，39(5)：795-805.
WU Qiang.Progress,problems and prospects of prevention and control technology of mine water and reutilization in China [J].Journal of China Coal Society,2014,39(5):795-805.
[2] 钱鸣高，石平五.矿山压力与岩层控制[M].徐州：中国矿业大学出版社，2003.
[3] 李 博，武 强.煤层底板突水变权脆弱性评价模型参数灵敏度分析[J].采矿与安全工程学报，2015，32(6)：911-917.
LI Bo,WU Qiang.An analysis of parameters sensitivity for vulnerability assessment of groundwater inrush during mining from underlying aquifers based on variable weight model [J].Journal of Mining & Safety Engineering,2015,32(6):911-917.
[4] 李 博，武 强.煤层底板突水危险性变权评价理论及其工程应用[J].应用基础与工程科学学报，2017,25(3):500-508.
LI Bo,WU Qiang.Risk evaluation of coal floor water inrush based on variable weight theory and its application [J].Journal of Basic Science and Engineering,2017,25(3):500-508 .
[5] 张金才，张玉卓，刘天泉.岩体渗流与煤层底板突水[M].北京：地质出版社，1997.
[6] SANTOS C F,BIENIAWSKI Z T .Floor design in underground coalmines[J].Rock Mech Rock Eng,1967,22(4):249-271.
[7] LI Bo,WU Qiang,DUAN Xianqian,et al.Risk analysis model of water inrush through the seam floor based on set pair analysis [J].Mine Water & the Environment,2017.DOI：10.1007/s10230-017-0498-5.
[8] 黎良杰，钱鸣高，李树刚.断层突水机理分析[J].煤炭学报，1996，21(2)：119-123.
LI liangjie,QIAN minggao,LI shugang.Mechanism of waterinrush through fault [J].Journal of China Coal Society,1996,21(2):119-123.
[9] 卜万奎，茅献彪.断层倾角对断层活化及底板突水的影响研究[J].岩石力学与工程学报，2009，28(2)：386-394.
BU Wankui,MAO Xianbiao.Research on effect of fault dip on fault activation and water inrush of coal floor[J].Chinese Journal of Rock Mechanics and Engineering,2009，28(2)：386-394.
[10] 浦 海，张 柬.断层影响下顶板突水特性的数值模拟研究[J].采矿与安全工程学报，2010,27(3)：421-424.
PU Hai,ZHANG Jian.Numerical simulation of water inrush from coal roof affected by fault [J].Journal of Mining & Safety Engineering,2010,27(3):421-424.
[11] VICTOR Bense,BERG E H Van den,RONALD Van Balen.Deformation mechanisms and hydraulic properties of fault zones in unconsolidated sediments: the Roer Valley Rift System,The Netherlands [J].Hydrogeology Journal,2003,11(3):319-332.
[12] 武 强，朱 斌，李建民，等.断裂带煤矿井巷滞后突水机理数值模拟[J].中国矿业大学学报，2008,37(6)：421-424.
WU Qiang,ZHU Bin,LI Jianmin,et al.Numerical simulation of lagging water-inrush mechanism of rock road ways near fault zone [J].Journal of China University of Mining & Technology,2008,37(6):421-424.
[13] 彭苏萍，孟召平，王 虎，等.不同围压下砂岩孔渗规律试验研究[J].岩石力学与工程学报，2003，22(5):742-746.
PENG Suping,MENG Zhaoping,WANG Hu,et al.Experimental study on pore ratio and permeability of sandstone under different confining pressures [J].Chinese Journal of Rock Mechanics and Engineering,2003,22(5):742-746.
[14] ELSWORTH D,BAI M.Flow-deformation response of dual-porosity media [J].Journal of Geotechnical Engineering,1992,118(1) :107-124.
[15] RUTQVIST J,TSANG C F.A study of cap rock hydro-mechanical changes associated with CO2-injection into a brine formation [J].Environmental Geology,2002,42:296-305.
[16] MA D,BAI H B,MIAO X X,et al.Compaction and seepage properties of crushed limestone particle mixture: an experimental investigation for Ordovician karst collapse pillar groundwater inrush [J].Environmental Earth Sciences,2016，75(1):1-14.
[17] PRADIPKUMAR G N,VENKATARAMAN P.Non-Darcy converging flow through coarse granular media [J].Journal of the Institution of Engineers (India) : Civil Engineering Division,1995,76: 6-11.
[18] 师文豪，杨天鸿，刘洪磊，等.矿山岩体破坏突水非达西流模型及数值求解[J].岩石力学与工程学报， 2016，35(3)，446-455.
SHI Wenhao,YANG Tianhong,LIU Honglei,et al.Non-Darcy flow model and numerical simulation for water-inrush in fractured rock mass [J].Chinese Journal of Rock Mechanics and Engineering,2016,35(3),446-455.
