0.   引　　言

水岩相互作用是地质学和岩土工程领域的研究热点之一。水岩作用过程中涉及许多物理和化学反应，主要包括润滑、沉淀、氧化还原、离子交换等^[1–4]。许多地质灾害和工程安全问题都与水岩作用有关，如边坡稳定性^[5]、水库大坝稳定性^[6]、岩爆^[7]、岩溶塌陷^[8]、矿井突水^[9]、巷道失稳等^[10]。水对岩体力学和渗流行为的影响是分析和解释上述问题的基础。近年来，通过大量试验研究发现，水的存在削弱了岩石的力学参数(如弹性模量、抗压强度、黏聚力和抗拉强度等)以及脆性，也改变了岩石破坏后的块体分布^[11-12]。一方面，水的存在增加了土壤和岩体中的孔隙压力，从而降低了介质中骨架颗粒的有效应力，进而改变了岩石的物理力学参数。另一方面，水可以溶解矿物成分和颗粒之间的黏结物，从而增加孔隙，弱化黏聚力。如果水中含有一些亲水性矿物成分，水岩作用将更加严重^[12-13]。

煤系地层中有许多含水层，煤矿开采活动对这些含水层造成扰动。因此，水不可避免地会流入煤矿开采影响区域，威胁煤矿安全开采，甚至引发突水事故。除此之外，淋水区域将进一步弱化煤岩体，造成巷道、煤柱等失稳事故频发^[14]。随着测试技术，微观表征技术的发展，关于煤岩体水岩作用力学性质以及水岩作用微观机理的研究已经非常丰富^[11]。但现阶段的研究主要集中在实验室尺度，岩体强度的跨尺度效应与边界环境的差异，导致很多研究成果无法直接应用于工程现场问题。因此，如何实现实验室成果的扩尺度工程应用至关重要。随着计算机软硬件的发展，基于实验室试验结果的数值模拟正成为工程尺度煤矿开采问题分析的主要研究方法。然而，以往的研究主要集中在煤矿开采原岩与采动应力、开采工艺与参数、岩层赋存特征、地质构造等因素对围岩稳定性的影响^[15–18]，很少涉及煤矿开采水浸软化等因素的影响。除此之外，在工程尺度流固耦合模拟方面，也大多只关注渗流场和力学场之间的相互耦合，很少考虑水对煤岩体的弱化作用^[19-20]。或者为了计算方便，在考虑水浸弱化作用时不进行渗流场的分析。比如POULSEN等^[21]根据煤体干燥和饱和时的力学性质研究水浸对煤柱稳定性影响，但在模拟中没有考虑含水率以及水渗的影响。白庆升等^[22]则根据水浸时间与煤岩体强度的关系模拟富水巷道顶板变形机理，这一方面没有考虑水渗的影响，另一方面模拟时间很难和实际时间相对应。

基于此，笔者基于常用的离散元模拟方法，提出了一种适合工程尺度的水浸软化流固耦合模拟方法，并讨论了其在有限元模拟方法实现中的区别。在此基础上，模拟分析了充水断裂构造区巷道围岩弱化特征并验证了模拟方法的可靠性。最后，基于提出的模拟方法，讨论了其适用性和优缺点，为后续研究提供了基础。

1.   数值模型介绍

数值模拟是煤矿开采水岩作用分析的一种必不可少的方法，尤其是在工程尺度问题中。根据水岩作用在煤矿开采中研究目的不同，模拟可分为连续岩体模拟和裂隙岩体模拟。连续岩体模拟主要集中在巷道围岩及煤柱失稳等问题。裂隙岩体模拟包括断层突水、导水裂隙发育等问题。实际上，2种模拟方法在一定假定条件下均可以应用。基于此，笔者给出了基于离散元的水浸软化流固耦合模拟计算方法，并进一步扩展至有限元模拟方法。

1.1   基于离散元的水浸弱化流固耦合模型

离散元数值模拟方法因其能够很好的模拟裂隙发育特征而广泛应用于水压致裂、实验室加载渗流、导水裂隙发育渗流等方面的流固耦合模拟^[23–25]。在离散元水浸弱化模拟过程中，主要考虑渗透率的更新和煤岩体强度弱化2个方面。

1.1.1   渗透率计算模型

以煤矿开采中最常用的离散元模拟软件UDEC为例进行离散元流固耦合分析。UDEC离散元数值模拟软件能够对不渗透块体间的流体流动进行分析。该软件能够进行全流固耦合分析，其中节理的渗透率取决于固体变形，而流体压力则反过来影响力学计算。在渗流计算过程中，渗流规律符合常规的立方体定律，渗流速率q可以根据上下游压力差ΔP、节理长度l、节理开度a和流体动力黏度μ进行计算；而节理开度a可以由节理法向刚度k_n、法向应力σ_n以及初始节理开度a₀^[23]计算获得，具体见式(1)和式(2)。

根据式(2)可以看出，随着应力的增加，节理开度呈线性减小，为了确保节理开度大于0，需要设置残余节理开度a_res，具体如图1所示。

图  1  节理开度与法向应力的关系

Figure  1.  Relationship between joint opening and normal stress

下载: 全尺寸图片 幻灯片

节理开度随应力的变化不断变化，通过节点的流量也随之变化，孔隙压力P可以根据该区域内所有裂隙流量Q、流体的体积模量K_w、运算时步Δt和原始区域的压力P₀计算：

其中，V_m=(V+V₀)/2；ΔV=V–V₀，V及V₀分别为新旧区域的面积。运用式(1)—式(3)就可以进行相应的流固耦合计算，模拟过程如图2所示。

图  2  块体离散元流固耦合模拟方法

Figure  2.  Main fluid-solid coupling ways of discrete element method

下载: 全尺寸图片 幻灯片

由式(2)及对应的图1可以看出，要保证在高应力阶段节理开度始终为正，则需要设置相应的残余节理开度(a_res)。而如果煤层埋深较大，地应力也会很大，按照上面的方法，煤样节理开度近似接近保留在残余节理开度，流固耦合的结果会被影响。得到这种结果是因为计算过程中认为k_n不变，这与试验结果不符。许多专家通过试验等研究发现，由于有效应力的改变，节理的可压缩性系数(c_f)也会发生变化，进而造成节理刚度也在不断变化，所以在模拟过程中应该根据有效应力不断更新节理刚度，需要增加节理压缩系数的改变系数α_f ^[26-27]，具体为

假设开始阶段应力为0，则式（5）可化简为

由式(6)可以看出，随着有效应力的增加节理刚度同样增加，能够避免在高应力处节理开度一直维持在残余状态的情况。式(6)中的节理初始可压缩系数c_f0和改变系数α_f可以通过数值反演实验室应力–渗流试验确定^[28]。需要说明的是，相对于固定节理刚度模拟方法，变节理刚度因为需要实时更新节理刚度，使得模型的计算量会一定程度的增加。因此，在模型应力变化不大的工程问题分析过程中仍可以采用固定节理刚度模拟方法。

1.1.2   煤岩体弱化计算模型

在UDEC模拟过程中，为了实现块体划分的均匀性和随机性，一般采用Voronoi多边形块体模型。二维Voronoi模型由单元块体(Block)和块体间的节理(Joint)组成。在UDEC数值模拟中，假定块体是不渗水的，水只能通过块体之间的裂缝渗流。在采掘扰动过程中，围岩应力环境发生变化，块体将出现旋转、移动、变形等运动形式，而节理开度将发生变化，如图3所示。根据上述流固耦合模拟方法可知岩层内地下水的流动取决于节理的开度以及含水层的水压。由于块体离散元模拟过程中煤岩体内的水主要存储与煤岩体节理裂隙中，则煤岩体内含水率的变化与节理裂隙的开度有关。

图  3  含水率计算方法

Figure  3.  Water content calculation method

下载: 全尺寸图片 幻灯片

由于Voronoi模型为不规则块体划分，为了有效统计煤岩含水率变化规律，可以将研究区域的煤岩体划分成1 m×1 m的单元体矩形网格，单元体含水率为

其中，ρ₁为水的密度；ρ₂为单元体内块体的密度；$\displaystyle\sum V $为单元体内块体的体积。由于煤岩体存在饱和含水率，因此，在建模过程中需要根据式(7)和饱和含水率确定单元体内煤岩体初始节理开度和长度。

在含水率获得之后就可以根据含水率更新煤岩体强度。在UDEC–Voronoi模型中岩石的力学行为由3组参数控制，分别为强度控制参数(节理和块体的内摩擦角、抗拉强度和黏聚力)、变形控制参数(块体的体积模量和剪切模量、接触的法向及切向刚度)、渗流控制参数(初始及残余开度和节理渗透率)。

为了准确获得数值模型中煤岩体在水浸弱化作用下的参数变化，可以将数值结果与煤岩体在不同含水率条件下实验室数据进行反演。当块体单元尺寸足够小时，数值结果(如峰值抗压强度和抗拉强度)几乎与单元尺寸和网格划分无关。为此，采用块体单元边长为0.002 m的数值试件，对不同含水率煤体和岩体分别进行单轴压缩和劈裂试验。块体弹性模量E可以根据应力应变曲线获得，并通过块体质心之间的距离与弹性模量的关系推导出节理弹性参数。利用单轴压缩的应力–应变曲线以及试样的破坏形式，通过试错优化算法，确定煤岩体的黏聚力和内摩擦角。然后将数值结果与所有含水率的试验数据进行比较，直到模拟结果与试验数据匹配，具体如图4所示。

图  4  在不同含水率条件下煤岩体的弱化情况

Figure  4.  Weakening situation of coal and rock masses under different water content conditions

下载: 全尺寸图片 幻灯片

综合变节理刚度渗透率更新模型和水浸软化模型可以实现离散元的水浸软化流固耦合模拟，具体模拟程序如图5所示。与以往流固耦合主要的区别在于：① 在模型构建过程中需要根据煤岩体饱和含水率和渗透率设置裂隙尺寸；② 在强度弱化方面区分节理和裂隙，在节理破坏前，其强度根据含水率进行更新，破坏后则根据强度折减法进行强度更新；③ 在渗透率更新方面，同样区分节理和裂隙，并引入了变节理裂隙刚度渗透率更新方法。

图  5  离散元水浸弱化流固耦合模拟程序

Figure  5.  Program diagram of fluid-solid coupling water leaching weakening of coal and rock mass

下载: 全尺寸图片 幻灯片

1.2   基于有限元的水浸弱化流固耦合模型

相比于离散元模拟方法，有限元模拟方法不存在节理裂隙，渗流发生在单元体内。因此，有限元模拟方法在计算渗透率和含水率方面与离散元模拟方法不一致。

在渗透率更新方面，与离散元类似需要首先判断块体是否发生屈服。在块体屈服前，需要考虑渗流块体的各向异性设置横向和纵向渗透率，煤体的可以参考文献[29]获得的煤样轴向和径向渗透率与应力的关系。在块体屈服后，其渗透率将发生大幅度增加，需要采用裂隙煤岩体应力–渗透率计算关系，具体可以参考文献[30]给出的贯穿裂隙渗透率应力关系。

在含水率更新方面，有限元同样只更新弹性块体的含水率，对于塑性块体的强度弱化则采用应变软化模型更新(图5)。弹性块体的饱和含水率可以通过孔隙率进行更新，而不同应力状态下的孔隙率大小则可以根据更新后的渗透率进行计算。需要说明的是，煤岩体在水浸后并不是立马达到饱和含水率，如样品尺度的煤岩样一般在50 h左右能接近饱和含水率^[31–33]。因此，在模拟过程中需要考虑饱和前和饱和后2种状态，其饱和周期可以通过流体渗流速度设定确定。在获得煤岩体最终含水率后，可以通过文献[34-35]中实测获得的煤岩体强度参数与含水率的关系更新数值模型块体强度参数。

2.   水浸弱化流固耦合模型的工程应用

虽然本文提出的水浸弱化流固耦合模拟方法能够同时适用于离散元模拟方法和有限元模拟方法，但考虑到篇幅限制，且主要对离散元模拟方法进行了详细介绍，因此，主要进行离散元水浸弱化流固耦合模拟方法的应用于验证。

2.1   工程概况

P8201为北辛窑煤矿12采区首采工作面，P8201走向长度为1 546 m，倾向长度为190 m，此工作面采煤方法选择的是综合机械化放顶煤采煤法，其中，采高为3.9 m，放煤高度为3.1 m，工作面布置如图6所示。工作面直接顶为6.1 m厚的中粒砂岩，基本顶是均厚12 m的细粒砂岩，均厚15 m的粉砂岩和细粒砂岩分布在底板。2号煤层上方为二叠系砂岩和奥灰含水层，其中，二叠系砂岩含水层距2号煤层16～20 m，奥灰含水层距2号煤层45～55 m，两层含水层局部导通，顶板承压水水压实测值为1.8～2.1 MPa。

图  6  P8201工作面布置

Figure  6.  P8201 working face layout plan

下载: 全尺寸图片 幻灯片

回风和运输巷道均沿煤层底板掘进，平行布置于工作面两侧。巷道断面规格为5.6 m×3.9 m(宽×高)，断面面积为21.84 m²。巷道顶板支护使用锚杆锚索，巷道两侧使用钢丝网和锚杆支护的方式。P8201工作面内包括14个断层，而断层落差大于10 m的断层有4条，如图6所示。其中正断层F2、F6分别位于运输和回风巷道掘进的路径上，其落差高度分别为15、14 m。由于断层沟通顶板含水层，导致巷道淋水不断，围岩变形量大。为了研究充水断裂构造围岩弱化特征，采用本文提出的离散元水浸弱化流固耦合模拟方法进行模拟。

2.2   巷道数值模型构建

根据北辛窑工作面P8201巷道建立2D Voronoi模型来模拟巷道围岩水浸弱化渐进破坏的过程以及力学行为的变化，如图7所示。由于采用二维模拟，选择三维掘进过程中巷道最不稳定的截面进行模拟。模型长为100 m，高为90.9 m。为了研究巷道围岩的破坏过程，将研究区域设置边长为0.4 m的不规则块体，同时考虑到模型的计算效率，将研究区域以外的模型划分为矩形块体。根据现场4处原岩应力测量结果设置模型垂直应力5 MPa，水平应力为4 MPa，模型侧面限制水平位移，模型底部限制垂直位移。考虑篇幅，模型模拟参数确定和修正参考文献[36]。

图  7  充水断裂构造巷道数值模型与边界条件

Figure  7.  Numerical models and boundary conditions

下载: 全尺寸图片 幻灯片

2.3   水浸弱化渐进损伤特征

为分析水浸弱化对巷道围岩稳定性的影响程度，根据图5所示的水浸弱化条件下巷道围岩稳定性分析的数值模拟方法，即回风巷道开挖后进行200次循环(每次循环渗流和力学各计算100步)，模拟过程中实时监测模型各区域孔隙压力、围岩损伤量、应力位移分布情况和围岩水浸弱化参数。图8为不同迭代次数结束后，巷道围岩不同位置处孔隙压力、参数弱化和损伤变化情况，具体巷道围岩监测布置如图8c所示。图8c中的损伤量是含水率计算单元范围内裂隙发育长度占节理总长度的比值。

图  8  充水断裂构造巷道围岩损伤弱化及孔隙水压演化特征

Figure  8.  Characteristics of pore water pressure evolution and surrounding rock damage weakening in water filled fracture structure roadway

下载: 全尺寸图片 幻灯片

在UDEC模拟中，巷道围岩的含水率和渗透率仅与节理的开度有关。由图8可以看出，由于断层的存在，造成巷道围岩稳定性差，巷道围岩破坏严重，A₁、A₄和A₅的损伤量在巷道开挖后急剧增加。位于断层内的A₁监测点产生张拉裂隙，损伤量达到51.5%，力学参数下降为残余强度。此时，含水层水并未渗流至巷道围岩(图9a)，除A₁以外，其他监测点力学参数并未弱化。随着流固耦合的继续运算，含水层水沿着断层裂隙面向巷道入渗(图9b)，测点A₁和A₂的水压开始升高。这导致A₂测点处围岩力学参数出现弱化(图8b)，进一步导致该区域损伤量急剧增加至53.4%。随着断层内岩体不断水浸弱化，裂隙进一步发育，含水层渗水区域扩展至A₄和A₅(图9c、图9d)，进一步造成上述2处范围内岩体黏聚力和内摩擦角下降，并最终达到残余强度(图8b)，损伤量大幅度提升。可以看出，水浸弱化会造成巷道围岩参数弱化–裂隙发育–加速水渗–参数弱化的渐进损伤特征。除此之外，位于断层上方的A₃监测处由于裂隙未贯通到该区域，此处水压一直为0(图8a和图9d)，进而该处力学参数并未弱化(图8b)，这从侧面反映出本文提出的水浸弱化流固耦合模型的可行性。

图  9  水浸弱化条件下巷道围岩孔隙压力变化过程

Figure  9.  Change process of pore pressure in surrounding rock of roadway under water immersion weakening conditions

下载: 全尺寸图片 幻灯片

2.4   模拟可靠性分析

为了进一步验证水浸弱化流固耦合模拟的可靠性，图10给出了P8201巷道顶板下沉量与数值模拟的对比结果。由于数值模拟很难有效模拟岩层内块体之间的最大间距，导致模拟结果略小于现场测得的结果。但由图10可以看出，顶板沉降的数值模拟结果与现场结果基本一致，证明了水浸弱化模拟方法研究充水断裂构造巷道渐进破坏特征的适用性。值得注意的是，UDEC计算时间和现场监测时间并不相等，两者之间存在线性相关性。

图  10  充水断裂构造巷道顶板下沉量的数值模拟和现场实测对比

Figure  10.  Comparison of numerical roof subsidence and field result in water filled fracture structure roadway

下载: 全尺寸图片 幻灯片

除了巷道变形量，巷道裂隙发育情况也是验证水浸弱化模拟可靠性的重要方面。图11为充水断裂区域巷道倾斜钻孔和垂直钻孔窥视结果。在窥视过程中，充水断裂构造区域成孔率差，钻孔周边围岩裂隙发育，破碎严重(图11a)。距顶板4 m处的垂直孔内可以看到断层面内为不稳定的松动小石块，层内离层清晰可见。随着钻孔高度的增加，砂岩地层中出现明显的水浸弱化面，这与图9中的孔隙压力分布基本一致，并在孔壁上产生垂直和水平裂隙(图11b)。对比图11c数值模拟获得的最终裂隙发育可以看出，数值模拟结果与现场裂隙发育范围基本一致，进一步验证了本文模拟方法的可靠性。

图  11  充水断裂构造区巷道围岩损伤实测与数值模拟结果

Figure  11.  Measurement and simulation results of surrounding rock damage in water filled fracture structure roadway

下载: 全尺寸图片 幻灯片

3.   讨论与展望

煤岩体水浸弱化作用是煤矿开采围岩失稳的重要因素，大量实验室研究表明水岩作用将大幅度弱化煤岩体强度，并给出了煤岩体力学参数与含水率的量化关系。如何将实验室获得的量化结果应用于现场工程实践至关重要。基于此，笔者提出了考虑水浸软化的流固耦合模拟方法，该方法主要解决两方面的问题：一是渗透率根据采动应力和煤岩体屈服状态的实时更新；二是煤岩体强度根据含水率的实时更新。在此基础上，模拟分析了充水断裂构造区域巷道水浸软化特征，结合现场实测进一步验证了本文模型的可行性。除了应用于煤矿开采相关水浸弱化工程问题的分析，笔者提出的流固耦合模拟方法同样适用于其他水浸工程问题。如降雨影响下的边坡稳定性，水库大坝的稳定性，富水地层隧道硐室稳定性分析等。但笔者提出的水浸软化流固耦合模拟方法仍然存在需要改进的地方。

1)提出的水浸软化的流固耦合模拟方法采用的是间接流固耦合模拟方法。在固定的计算步骤之后，渗流和力学计算被分离并相互循环。这种方法的优点是该计算相对简洁高效，计算量小，适用于工程尺度的数值模拟。但实际情况下，渗流和力学计算是同时进行的，这使得在间接流固耦合计算时设置合理的间隔步数非常重要但又难以确定。但对于工程尺度，特别是煤矿开采这种采动损伤特别大、尺度特别大的场景内，很难实现全耦合计算。

2)在实现水浸弱化时采用的是含水率和煤岩体强度的相关关系。但实际上，除了含水率，干湿循环次数、浸水时间、渗透水压等都会影响煤岩体强度。且对于煤矿开采而言，长期水浸对于煤岩体强度的弱化更贴近现场实际情况。但是利用水浸时间进行煤岩体强度弱化存在2个方面的问题。一方面，现阶段关于长期水浸(1 a以上)的试验结果相对较少，且很少有煤岩体强度与水浸时间相关量化关系。另一方面，数值模拟中关于时间的标定并不容易，现有关于数值模拟时间的标定主要采用现场实测结果与模拟结果进行对比。但数值模拟计算时间受到模型大小，网格划分影响，很难与实际时间匹配。

3)提出的水浸软化流固耦合模拟方法有助于实验室成果的工程尺度应用。但煤岩体强度与渗流特征本身就具有尺度效应。实验室获得弱化系数是否与工程弱化系数一致需要进一步确认。特别是煤岩体渗流特征，随着尺度的增加，煤岩体连续面也会增加，其渗流能力将逐渐增加。因此，实验室获得的渗透率更新模型和力学参数弱化模型需要进行尺度校验。除此之外，由于煤岩体赋存环境的差异性造成煤岩体力学和渗流性能具有工程特殊性。对于特定煤矿开采工程问题，这就需要针对性进行水浸弱化试验和渗透率应力敏感性试验，进而获得模拟过程中力学参数和渗透率更新公式。

4)提出的考虑水浸弱化流固耦合模拟方法可以适用于离散元和有限元模拟方法，能够很好的适用于现阶段煤矿开采常用的UDEC、PFC、FLAC3D等数值模拟软件。随着模拟方法的发展，很多软件能够实现有限元与离散元的耦合模拟(有限离散元法)，不需要像离散元事先划分裂隙网格，单元体能够根据破坏准则随机产生裂隙。因此，后续工作可以将现有考虑水浸弱化的数值模拟方法进一步扩展到有限离散元方法。

4.   结　　论

1)针对煤矿开采水岩作用弱化工程问题，提出了水浸弱化流固耦合模拟方法。该方法主要分为渗透率更新模型和水浸弱化模型，可以实现采动水浸过程中煤岩体渗透率和强度的实时更新。针对解决问题的不同，笔者将提出的水浸弱化流固耦合模拟方法可以同时适用于离散元模拟方法和有限元模拟方法。

2)离散元水浸弱化流固耦合模拟方法主要用于节理裂隙水渗弱化模拟。模型中采用变节理裂隙刚度方法进行渗透率更新以适应煤矿开采应力扰动大的特点。在节理破坏前，煤岩体强度弱化则根据含水率与煤岩体强度的关系进行弱化，并利用节理开度计算煤岩体实时含水率。在节理破坏形成裂隙后，则采用强度折减法更新裂隙强度。

3)有限元水浸弱化流固耦合模拟方法在渗透率更新中同时考虑了应力、屈服状态和煤岩体的各向异性。在块体破坏前，煤岩体含水率主要通过孔隙率与渗透率的关系进行更新，之后进一步更新煤岩体力学参数。在块体破坏发生塑性变形后，则采用应变软化模型更新块体强度。

4)运用本文提出的水浸弱化流固耦合模拟方法模拟分析了充水断裂构造区域巷道围岩稳定性。实现了充水断裂构造区域巷道围岩渐进损伤破坏和参数弱化，并利用现场巷道顶板下沉量和围岩损伤裂隙发育范围实测结果验证了模拟的可靠性。

5)主要采用间接流固耦合模拟方法，虽然能够提高计算效率但降低了模拟精度。煤岩体弱化主要是根据含水率进行更新，但煤矿开采也要考虑水浸时间的影响，这需要解决模拟时间和实际时间匹配性问题。最后，渗流和弱化模型均基于实验室实测获得，后期需要考虑实验室尺度和工程尺度之间的相关性。