天然岩体中存在大量原生裂隙、孔隙，这些结构缺陷在地应力的作用下成为应力集中-释放区域，不仅大幅降低工程岩体的稳定性，也严重影响岩体的渗透特性。裂隙岩体的渗流场受区域地应力环境的影响很大，而渗流场的改变反过来又对地应力场产生作用。因此，“应力-渗流”耦合作用是裂隙岩体力学中一个重要的特征[1]。
乌鲁木齐矿区赋存着大量的急倾斜煤层，工作面上覆采动覆层中多存在由季节性水源(如冬季积雪融化)所形成的承压水体[2]。“裂隙-渗流”耦合下煤岩体结构性畸变导致岩体强度劣化，采动覆层中形成宏观裂隙网络，为承压水体向工作面运移提供流动通道;特别是承压水体渗流过程中常携带大量泥土、沙石，造成工作面突水溃沙现象时有发生;此外，岩体在经过较长时间水体浸泡后，其力学性能急剧劣化，采动覆层出现大尺度突发性断裂，大体积冒落和大面积片帮频发并诱发冲击地压，更甚者会出现井下泥石流等动力灾害，对工作面安全生产带来巨大隐患。
乌鲁木齐矿区地质力学环境极其复杂，矿区内急倾斜煤层的地层赋存环境在世界范围内都极具特殊性，区域地质力学特征以强挤压应力为主，“区域地应力-裂隙”耦合作用的存在进一步增加煤岩体的强度劣化程度;此外，采用连续介质的方法求解煤岩体渗流相关问题时，裂隙宽度的微小改变将引起裂隙渗流量的巨大变化，而裂隙隙宽受裂隙面上的应力控制，所以通过裂隙面的渗流量与应力环境紧密相关。另一方面，渗流水压也可以影响裂隙面上的有效应力。综上分析可知，为了准确分析岩体裂隙中的渗流问题，在进行岩体渗流研究时应考虑岩体结构、区域地应力与渗流相互耦合作用下的煤岩体强度劣化机理，揭示考虑区域地应力特征的裂隙煤岩流固耦合机理具有重要的现实意义。
已有研究表明[3-8]:水-岩体相互作用在力学领域为“流固耦合作用”，其中研究水-岩体相互作用的多场耦合理论模型与数值模拟方法等为主要研究内容之一，主要代表性学者如ODA[9]、唐春安[10]、仵彦卿[11]、王媛[12]等;在地球科学领域，水-岩体相互作用研究主要涉及水对岩体物理力学性质、破裂特性等方面，代表性学者如冯夏庭[13]、黄润秋[14]、唐辉明[15]、汤连生[16]等。研究成果主要集中在岩体水力学的流固耦合模型和水-岩体相互作用的力学性能及水物理化学作用对力学性质的影响。而针对我国西北生态脆弱区特定地质条件与生态地质环境下区域生态-水-岩体空间赋存结构特征、水-岩体相互作用下区域地应力特征-结构-渗透率变化规律及采矿诱发的动态响应机制，研究成果还相对较少，有待进一步深入研究。
基于以上认识，本文通过位移反演方法对区域地应力特征进行理论分析，建立研究煤岩内部裂隙时空演化规律的力学实验系统，设计了区域地应力作用下裂隙煤岩流固耦合计算模型，初步揭示裂隙煤岩流固耦合特性。
1 耦合计算模型的设计
1.1 区域地应力的反演计算
煤岩体地应力特征-结构-渗透率变化规律受多种因素影响，为顺利进行应力-裂隙-渗流作用下煤岩强度耦合特性数值实验，需要完成的任务主要包括:① 反演计算、得出乌鲁木齐矿区区域地应力变化规律;② 揭示煤岩体在实际工程扰动环境影响下裂隙场的时空演化特征;③ 在上述工作的基础上，基于耦合计算程序，完成耦合数值模型的构建、计算及结果分析。此外，为保证实验结果的合理性与一致性，依照ISRM制定的相关岩石测试标准，力学实验中所采用的试样与耦合数值模型均为边长60 mm的立方体试样且模型所承受载荷均按照文中所述的区域地应力变化规律施加。
文中进行区域地应力的反演计算位置为乌鲁木齐矿区乌东煤矿+500水平45号煤层综采工作面南巷。巷道为煤巷、断面为矩形，巷道断面尺寸(宽×高)为4.4 m×3.4 m。根据前期研究成果[17]:急倾斜煤层中开挖巷道相当于在煤岩体开挖边界上作用有等效结点力。其中，等效结点力与地应力之间满足体积积分函数关系:
(1)
其中,
为作用在煤岩体开挖边界上的等效结点力;σ为待求的区域地应力
为煤岩体开挖后整个开挖边界所产生的应变矩阵，取决于开挖煤岩体的几何尺寸。在反演计算过程中，南巷被简化成二维平面问题，即σ={σx,σy,τxy}T，地应力中的σz方向平行于巷道走向方向，不列入计算范畴。因而，平面问题中式(1)可简化为
(2)
式中，Bi(i=1,2,3)为B的分列函数。此外，在巷道开挖区域内，结点力与结点位移量u的量化关系为
(3)
式中，K为整个煤岩体开挖区域的刚度矩阵。在煤岩体开挖边界上P与相等，而开挖煤岩体内部其它结点上外力均为0。联立式(2)与式(3)得:
Ku=σxB1+σyB2+τxyB3
(4)
式中，Bi为巷道开挖边界的应变矩阵的分量。假设巷道位移量测量点与有限元网格结点相重合，结点位移量u可分为已知的测量点位移量u1与未测量的位移量u2，如式(5)所示。
u={u1,u2}T
(5)
联立式(4)和(5)，并消去u2，可得
(6)
式(6)中含有3个未知待求量:σx/E，σy/E与τxy/E，其中E为岩体的弹性模量为常数。现场布置煤岩体移动量及断面收敛量的量测点。于南巷150，200及250 m处各布置1组围岩移动传感器，共计3组。各断面内按顺时针方向埋设6个GYW-300型围岩移动传感器。巷道开挖后立即开始各量测点位移值及收敛值量测。现场量测数值代入式(6)可直接求出相应的开挖煤岩体边界上地应力数值，计算结果见表1。综上所述，得出乌东煤矿+500水平45号煤层综采工作面南巷地应力反演计算结果为:垂直地应力、最小水平主应力及切应力分别为7.057，8.085，0.057 MPa。
表1 区域地应力反演计算结果
Table 1 Relevant results on back analysis of regional geostress
为了校核区域地应力反演计算的准确性，现场采用KJ743型煤矿地应力监测系统同步采集地应力数据。地应力监测传感器分别位于南巷180，220，250，290，330，360，400及430 m处，每组传感器均包括深、浅基点传感器各1个。其中，深、浅基点的埋设深度分别为10.0,5.0 m。
图1为地应力传感器地应力数值的实时监测结果，其中图1(a)为地应力传感器于2014-11-11至2015-02-02所采集垂直应力数据。深、浅基点传感器的地应力数值变化基本保持一致并在监测时段内出现垂直地应力的小幅度波动，均可视为工作面开采扰动应力的影响。实测垂直地应力数值在6.8～7.0 MPa，略小于反演计算结果。图1(b)，(c)为地应力传感器获取的最大、最小水平主应力分布。可以看出最小水平主应力数值为7.6～8.4 MPa，与反演计算结果相近。340～400 m内小范围最小水平主应力出现应力数值的小幅度增大，可判定为开采扰动影响所致。最大水平主应力分布较均匀，基本维持在10 MPa左右，400～480 m靠近南巷处，最大水平主应力增大，约为13 MPa。此外，工作面220～290 m范围内最大水平主应力出现局部性突增，250 m处数值最大，达到28.8 MPa，最大水平主应力实测结果也可为后续耦合数值计算科学设计提供参考依据。
1.2 裂隙场时空演化特征
目前，开展有关急倾斜煤层中煤岩体裂隙场时空演化特征的力学实验较少，主要是由于多数用于结构监测设备无法与加载设备配套使用。鉴于此，自主开发配套于DR扫描的岩石力学实验装置。图2(a)为煤岩体在加-卸载条件下裂隙场时空演化特征力学实验的设计方案及监测原理，所采用的DR扫描设备为DSM-80型X射线机(图2(b))，加-卸载条件下可实时获取煤岩体内部裂隙的二维多媒体图像。
图1 乌东煤矿+500水平区域地应力实测结果
Fig.1 Results on in-situ monitoring of regional geostress
加-卸载条件下煤体裂隙场时空演化特征的DR测试的实验条件为80 kVp，320 mA，80 ms，共进行3阶段加-卸载，分别为3.48-1.77，1.77-4.27-3.95和3.97-9.68-0 MPa。如图3(a)～(d)所示，当峰值载荷(σp)为3.48 MPa时，试样沿载荷作用方向发生裂隙的萌生现象，试样右下角处裂隙扩展并形成主裂隙，且卸载后主裂隙开度相比未受载时增加显著;当σp为4.27 MPa时，主裂隙已贯通试样，贯通方向为煤体试样宏观断裂面方向，由于裂隙持续性扩展，煤体右下角出现剪切、滑落，且卸载后该区域煤体率先破坏;当σp达到9.68 MPa时，煤体试样发生共轭型破裂，煤体中部残留煤柱出现应力集中现象，对比发现残留煤柱的X射线衰减程度较为显著，表明煤柱被进一步压实。
岩体的DR测试的实验条件为150 kVp，320 mA，50 ms，如图3(e)～(h)所示，未受载状态岩体试样内部无明显裂隙等原生缺陷。第1次加载至10.70 MPa(阶段1)时且卸载岩体试样后，分析DR扫描图像得出:大量裂隙萌生，新生裂隙均由试样边缘向中部扩展，且裂隙开度较小;当σp为14.65 MPa时(阶段2)，试样卸载后裂隙进一步扩展，垂直于外载方向，多条裂隙贯通并形成主裂隙，多条倾斜裂隙与主裂隙贯通。当σp达到20.12 MPa时，卸载岩体后横向宏观主裂隙已成形，并与倾斜裂隙形成X型裂隙场。对MATLAB软件进行二次开发，二值化处理DR扫描图像。MATLAB程序二值化处理过程中，DR扫描的图片透光区定义为1，非透光区为0，图3中也同时展示二值化处理后的结果。结果矢量化处理后将为后续耦合数值模型构建提供裂隙场量化参数。
图2 煤岩体裂隙场时空演化特征力学实验
Fig.2 Mechanical testing for spatio-temporal characteristics of fracture field of coal-rock mass
图3 煤岩实验过程外部情况、DR扫描及二值化处理结果
Fig.3 External settings of coal-rock mass and relative results of both DR scanning and binarization processing
1.3 耦合计算模型的构建
选用COMSOL Multiphysics软件构建考虑区域地应力特征的裂隙煤岩流固耦合特性的计算模型。COMSOL Multiphysics以有限单元算法(FEM)作为其理论基础，借助偏微分方程(单一物理场)及偏微分方程组科学地展现多物理场下的煤岩体耦合作用效应及强度劣化过程。
将煤岩体视为多相体，分析水体渗流作用下裂隙煤岩的孔隙率、渗透率、孔隙压力、有效应力与水流速度等参数变化规律。针对应力作用下裂隙中水体流动问题，通常采用平均化和整体“连续”属性来描述水体流动连续方程，达西定律可被用作控制方程并引入平均流速与总流量等整体平均化变量来描述水体的连续方程。但是，针对耦合作用下水体渗流问题，整体平均化将大大削弱计算结果的精确度。因此，考虑煤岩结构的差异性及应力时空演化特征，假定裂隙空间水体的密度恒定、温度恒定且不可压缩，引入Navier-Stokes方程为
-▽p+▽μ(▽ν+▽νT)=0
(7)
其中，p为煤岩体注水压力值，MPa;μ为水体黏性系数;ν为水体流动速度，m/s。参照过往的实验结果，水体流动速度一般在9.6×10-3～9.0×10-2 m/s之间。水体流速足够快并在煤岩体内部出现层流现象，裂隙与水体间出现滑移，以至无法忽略剪切能量耗散，需对式(7)进行修订，考虑采用岩体渗流的Brinkman方程适合描述煤岩体中水体快速流动，通过输入黏性、剪切参数描述裂隙的动能耗散。耦合计算模型所采用的Brinkman运动基本方程为
(8)
式中，ρs,ρ分别为煤岩与水体的密度;vs为运移速度，这里煤岩体可视为定位模型，因而vs=0;k为煤岩体的初始渗透率。
图4给出了耦合计算模型的设计思路，其几何尺寸为0.06 m×0.06 m，矩形子区域为本算例的计算面域。模型顶部施加载荷，煤体施加载荷值分别为3.48，4.27和9.68 MPa，岩体施加载荷值分别为10.70，14.65与20.12 MPa。模型底部为固定约束，模型左右边界为对称边界，无围压并限定位移量。水体由模型顶部进入、底部流出且顶部进口水体压力为已知量，煤体模型注水压力为1.5，3.5，4.5和6.5 MPa;相应地，岩体中注水压力分别为2.0，5.0，8.0和10.0 MPa。水体流动出口压力为水体压力，模拟过程中不考虑重力所引起的体积力。
图4 应力-裂隙-渗流耦合计算模型设计思路
Fig.4 Schematic of coupled calculation model for coal-rock mass under mechano-fracturing-hydro settings
2 实验结果及分析
2.1 平均孔隙率与渗透率
在综合考虑力学、结构和渗透等影响因素的基础上，分析孔隙率与渗透率、水体压力、水流速度、渗流流线方向及出口速度等参数变化机制，初步揭示研究区域地质与开采环境下水-岩体相互作用下区域地应力特征-结构-渗透率时空变化规律。图5给出了煤岩体的孔隙率与渗透率的分布。应力-裂隙-渗流耦合作用下煤岩平均孔隙率(φave)和渗透率(kave)计算式为
(9)
(10)
图5 煤岩体孔隙率与渗透率的分布
Fig.5 Distribution traits of both porosity and permeability
其中，φ0,k0分别为初始未加载条件下，煤岩体初始孔隙率与渗透率。其中，煤岩体参数值分别为7.35%与1.87×10-13 m2、3.26%与4.85×10-12 m2;abs为计算设定的极小值，为2.2×10-16;im(x,y)为图像函数。表2给出了煤、岩体的平均孔隙率和平均渗透率的计算结果。随着外载荷的持续增加，煤岩体体孔隙率与渗透率的变化规律基本一致:即孔隙率与渗透率均随外载荷的增加表现为持续性增加，最大渗透率可达到4.52×10-11 m2。外部载荷持续增加将导致煤岩内部微观裂隙进一步扩展贯通并出现宏观破裂。这是煤岩体的孔隙率与渗透率持续性增加的本质原因。
表2 煤岩体平均孔隙率和渗透率的计算结果
Table 2 Calculation results of porosity and permeability
2.2 水体压力分布特征
图6给出了不同荷载下煤岩体水体压力表面分布等值线。外荷载相同条件下等值线形态相似。但是，注水压力愈大，相应的水力梯度越大。水体压力最大值分别约为6.3和9.9 MPa。岩体水流的层流特点显著，注水、出水位置均出现应力集中现象，且外荷越大，应力集中位置越靠近模型出水位置区域;水体在通过岩体中部主要裂隙通道前，水体压力近似均布，水力梯度降低均匀;通过主裂隙通道后，水力梯度出现突降;相应地，煤体愈靠近进水区域，水力梯度越小。随y方向几何尺寸增加，水力梯度降低幅度显著增加。
2.3 水体流速分布特征
准确描述水体流速最大区域迁移过程可作为判定煤岩体裂隙场时空演化特征的重要方法之一。图7为煤岩体水体流速表面分布等值线。水体流速值与模型注水压力值、外部荷载值呈现正相关。煤体的水体流速极大值一般在模型右侧次生裂隙通道及模型两侧底角出水区域附近;相应地，岩体的水体流速极大值出现于模型中部主裂隙通道。煤体在完成第3阶段加-卸载前，水体流速最大区域位于模型的右侧，说明水体主要经该区域的裂隙通道导出模型;而煤体加-卸载后，水体流速最大值区域迁移至模型左侧，新生的裂隙通道流出水体量明显增多。岩体首次加-卸载后，最大流速区域仅在中部主裂隙通道附近，其它区域水体流速极小，可忽略不计，表明岩体内部并未形成整体性裂隙网;随着外部荷载值的持续性增加，模型的平均流速增加显著，流速最大区域逐渐向左侧底角出水区域转移。
2.4 渗流流线分布特征
图8为煤岩体水体渗流流线分布特征与水体流动方向示意。其中，图中黑色实线为水体流线，红色与蓝色箭头均为水体的流动方向，特别地，红色箭头为水体流动主方向。图中可以看出:注水区域的水体流动方向具有较好的一致性。煤体内渗流流线分布极为不规则，水体流动随机性大。岩体模型中流线分布情况刚好相反，沿x轴方向基本均布;外载的持续增加导致主裂隙以下区域的流线分布随机性增大。此外，水体流动主方向对应于流速最大区域;岩体裂隙通道附近出现大量水体流动聚集现象，中部主裂隙导致出水区域的水体流出量减少显著。随着外载的增加，岩体两侧底角的新生裂隙通道愈发明显，流动主方向开始向该区域发生迁移。
2.5 出水区域水体流速分量
图9中绘制了耦合计算模型出口区域y方向水体流速分量随x坐标变化的分布情况。模型出水区域均匀安设10个出水流速分量监测点。由于水体流向为y轴负方向，故流速分量标记为负值。从图中可以看出:模型出口流速分量的最大值均位于模型两侧;煤、岩体出口区域流速分量的最大值均随着外部荷载的增大而增大，相应的最大值分别为0.038与0.022 m/s。此外，不同注水压力下，模型出水区域的流速分量变化趋势基本一致，且注水压力越大、相应的出水压力越大。岩体的出水区域水体流速分量出现多个极值点。
3 结 论
(1)通过地应力的现场监测，验证了位移反演方法的准确性，即研究区域垂直地应力与最小水平主应力分别为7.057及8.085 MPa，且垂直地应力与最小水平主应力数值稳定，受开采扰动影响，最大水平主应力出现局部性突增。
(2)给出了不同载荷下内部裂隙场DR扫描图像，为耦合计算模型提供了量化裂隙场数据。发现随外载增大煤岩内主裂隙持续扩展且裂隙场分别为共轭型与X型。
(3)初步得出区域地应力影响下裂隙煤岩流固耦合特性:外载增大将导致煤岩内部微裂隙进一步相互贯通扩展并出现宏观破裂，为裂隙煤岩孔隙率与渗透率持续性增加的本质原因，且流速最大区域迁移过程可表征裂隙煤岩次生裂隙分布特点;此外，主裂隙内水体流速最大且随注水压力的增大水体流速不断地增加，较大的流速一般发生在具有较大水力梯度的细小裂隙位置处。同时，出口流速特征表现为出口流速分量的最大值均位于煤岩两侧并随着外载的增大而增加。
图6 煤岩体水体压力表面分布等值线(单位:MPa)
Fig.6 Contour of superficial water pressure of coal-rock mass(Unit:MPa)
图7 煤岩体表面水体流速分布等值线
Fig.7 Contour of superficial water velocity of coal-rock mass
图8 渗流流线分布特征与水体流动方向
Fig.8 Distribution of streamline and direction of water flowing
图9 出口区域y方向流速分量随x坐标变化规律
Fig.9 Variation of velocity of outlet area at y direction with x value changing
参考文献(References):
[1] 谢和平，高峰，鞠杨.深部岩体力学研究与探索[J].岩石力学与工程学报，2015，34(11):2161-2178.
XIE Heping，GAO Feng，JU Yang.Research and development of rock mechanics in deep ground engineering[J].Chinese Journal of Rock Mechanics and Engineering，2015，34(11):2161-2178.
[2] 张东升，李文平，来兴平，等.我国西北煤炭开采中的水资源保护基础理论研究进展[J].煤炭学报，2017，42(1):36-42.
ZHANG Dongsheng，LI Wenping，LAI Xingping，et al.Development on basic theory of water protection during coal mining in northwest of China[J].Journal of China Coal Society，2017，42(1):36-42.
[3] 朱万成，魏晨慧，田军，等.岩石损伤过程中的热-流-力耦合模型及其应用初探[J].岩土力学，2009，29(12):3851-3857.
ZHU Wancheng，WEI Chenhui，TIAN Jun，et al.Coupled thermal-hydraulic-mechanical model during rock damage and its preliminary application[J].Rock and Soil Mechanics，2009，29(12):3851-3857.
[4] 张玉卓，张金才.裂隙岩体渗流与应力耦合的试验研究[J].岩土力学，1997，18(4):59-62.
ZHANG Yuzhuo，ZHANG Jincai.Experimental study of the seepage flow-stress coupling in fractured rock masses[J].Rock and Soil Mechanics，1997，18(4):59-62.
[5] W Allen HATHEWAY,Yuji KANAORI,et al.Encompassing hydrogeology,environmental geology and the applied geosciences[J].Engineering Geology,2005,81(2):99-130.
[6] SCESI L,GATTINONI P.Roughness control on hydraulic conductivity in fractured rocks[J].Hydrogeology Journal,2007,15(2):201-211.
[7] ZHAO M,ZHOU H,CHEN D.Investigation and application on gas drive development in ultra-low permeability reservoirs[J].Journal of Hydrodynamics,2008,20(2):254-260.
[8] JAFARI M K,PELLET F,BOULON M.Experimental study of mechanical behavior of rock joints under cyclic loading[J].Rock Mechanics and Rock Engineering,2004,37(1):3-7.
[9] ODA M.An equivalent continuum model for coupled stress and fluid flow analysis in jointed rock masses[J].Water Resources Research,1986,22(13):1845-1856.
[10] 唐春安，马天辉，李连崇，等.高放废料地质处置中多场耦合作用下的岩石破裂问题[J].岩石力学与工程学报，2007，26(S2):3932-3938.
TANG Chun’an，MA Tianhui，LI Lianchong，et al.Rock failure issues in geological disposal of high-level radioactive sastes under multi-field coupling function[J].Chinese Journal of Rock Mechanics and Engineering，2007，26(S2):3932-3938.
[11] 仵彦卿.裂隙岩体应力与渗流关系研究[J].水文地质工程地质，1995，22(6):30-35.
WU Yanqing.Study of seepage law with stress in jointed rock masses[J].Hydrogeology & Engineering Geology，1995，22(6):30-35.
[12] 王媛，速宝玉，徐志英.裂隙岩体渗流模型综述[J].水科学进展，1996，7(3):276-282.
WANG Yuan,SU Baoyu,XU Zhiying.Comment on the models of seepage flow in fractured rock masses[J].Advances in Water Science，1996，7(3):276-282.
[13] 冯夏庭，赖户政宏.化学环境侵蚀下的岩石破裂特性—第一部分:实验研究[J].岩石力学与工程学报，2000，19(4):403-407.
FENG Xiating，MASAHIRO Seto.Rock fracturing behaviors under chemical corrosion-part I:Experimental study[J].Chinese Journal of Rock Mechanics and Engineering，2000，19(4):403-407.
[14] 黄润秋，徐德敏，付小敏，等.岩石高压渗透试验装置的研制与开发[J].岩石力学与工程学报，2008，27(10):1981-1991.
HUANG Runqiu，XU Demin，FU Xiaomin，et al.Development and research of high pressure permeability testing system for rocks[J].Chinese Journal of Rock Mechanics and Engineering，2008，27(10):1981-1991.
[15] TANG Huiming,JIA Hongbiao,HU Xinli,et al.Characteristics of landslides induced by the great wenchuan earthquake[J].Journal of Earth Science,2010,21(1):104-113.
[16] 汤连生，张鹏程，王思敬.水-岩化学作用之岩石断裂力学效应的试验研究[J].岩石力学与工程学报，2002，21(6):822-827.
TANG Liansheng，ZHANG Pengcheng，WANG Sijing.Testing study on effects of chemical action of aqueous solution on crack propagation in rock[J].Chinese Journal of Rock Mechanics and Engineering，2002，21(6):822-827.
[17] 单鹏飞.急倾斜煤岩体强度劣化与动力灾害防治基础研究[D].西安:西安科技大学，2015:101-102.
SHAN Pengfei.Fundamental research on coal-rock mass strength deterioration and its dynamic hazard prevention and control in mining steeply inclined coal seams[D].Xi’an:Xi’an University of Science and Technology，2015:101-102.>