连续加载应力下真实裂缝流场和渗透率演化规律数值研究

梁运培，张怀军，王礼春，秦朝中，田键，陈强，史博文

（1.重庆大学煤矿灾害动力学与控制全国重点实验室，重庆 400044；2.重庆大学资源与安全学院，重庆 400044；3.天津大学表层地球系统科学研究院，天津 300072）

水力压裂技术已被大规模运用到非常规油气、深部地热资源等地下能源开采中，且取得了显著成效[1-3]。通过水力压裂作业，原本致密的地层在高水压下产生大量人工裂缝，最终形成主裂缝、次裂缝以及天然裂缝相互串通的复杂裂缝网络体系[4-9]。由于尺寸匹配、颗粒悬浮等问题，压裂液中的一些支撑剂无法进入次裂缝以及天然裂缝中进行有效支撑[10]，这类未被有效支撑的裂缝极易在地层应力扰动下出现不同程度的闭合，从而基本上失去裂缝的导流能力。因此，研究裂缝在应力扰动下的渗透率演化特征成了揭示裂缝导流能力损伤机制的前提。

国内外学者通过开展裂缝样品应力（加载与卸载）实验，测试了不同应力作用下的裂缝渗透率变化特点[11-15]，这些信息对分析物质在裂缝内部的运移行为至关重要。但采用实验方法难以直接探究应力作用下裂缝内部流场、流动通道与速度场的演化规律，因此，越来越多的学者开始采用数值模拟方法揭示裂缝在应力作用下的导流能力变化机制。学者们通过数值生成的裂缝，发现应力扰动下裂缝的导流能力主要与裂缝的平均开度、表面粗糙度以及接触面积等因素有关，并在此基础上建立起表征裂缝渗透率的相关经验模型[16-20]。然而，通过数值随机生成的裂缝通常在裂缝几何形态上满足平均化特征，与人工裂缝和部分天然裂缝的几何形态表现出的强非均质性存在显著差别。因此，通过直接数值模拟研究非均质裂缝渗透率演化行为存在适用性受限问题。

本研究首先利用巴西劈裂法制备真实人工裂缝岩样，提取裂缝表面形貌基本特征参数后，再通过数值方法拼接裂缝表面形成具有非均质裂缝特征信息的单条裂缝结构，随后，研究裂缝在有效应力下的渗透率演化特征，分析裂缝流场、接触面积、空间相关长度等参数的响应行为。最后，将研究模拟结果与传统经典经验模型进行验证，证实了采用真实裂缝获得的非均质裂缝与直接数值生成的裂缝在渗透率演化行为以及渗透率模型匹配关系上存在明显差别。

1 单裂缝开度分布及非均质性

1.1 岩石裂缝提取

裂缝在页岩气藏和致密砂岩气藏中普遍存在，对天然气的采出发挥了重要作用。因此，选取页岩和致密砂岩作为研究对象，同时选取花岗岩作为对比样，分析3种典型岩石中裂缝的单相渗流行为。首先，对3 类岩石钻取直径为50 mm、长度为100 mm 的柱塞岩样用于开展力学性能测试分析，获取弹性模量和泊松比，为建立接触力学模型提供基准参数。随后，钻取直径为25 mm、长度为50 mm 的3 类岩石标准柱塞岩样用于人工造缝，获取真实裂缝的几何形貌等信息。

采用MTS815 岩石力学试验设备测量50 mm×100 mm样品的多组岩石力学参数，平均结果如表1所示。测试时，通过轴向位移控制加载方式，页岩和致密砂岩的加载速率为0.1 mm/min，考虑花岗岩强度更高、变形更慢，对花岗岩的加载速率设定为0.05 mm/min。

表1 不同种类岩石基本力学参数测试平均结果Table 1 Average results of basic mechanicalparameters of different kinds of rocks

分别选取25 mm×50 mm 的页岩、致密砂岩和花岗岩样品各2块，采用巴西劈裂法对样品进行人工造缝（图1）。造缝仪（型号AG-250kN IS）加载方式为轴向位移控制，试件的加载速率为0.1 m/min。为了使岩样均匀受力，造缝时将岩样水平放在凹槽中，造缝仪压头与岩样表面保持接触平整，使岩样受压形成张拉裂缝，随后，采用Cronos Dual 三维光学扫描仪对人工裂缝进行扫描（设备分辨率为300 万像素，扫描精度为0.02 mm），获取裂缝表面的空间位置等形貌参数信息。

图1 基于巴西劈裂法的人工裂缝样品制作实物Fig.1 Photos of artificial fracture sample via Brazilian splitting method

将每个裂缝样品上下2个表面的形貌数据，通过位置坐标进行拼接得到整合裂缝，整合裂缝的拼接需满足2 个条件：①在垂直方向上，选择上下2 个裂缝面的最高点和最低点，确保其在同一条竖直线上呈对准状态；②在水平方向上，当上下2 个裂缝表面拼接后，2个裂缝面轮廓线能够重合。通过插值方法对拼接后的裂缝表面进行长度x方向和宽度y方向网格单元划分，网格单元尺寸为l1（长度）×lw（宽度）=0.16 mm×0.16 mm。在扫描提取裂缝时，由于出现了部分边界损失，将这部分损失的边界切除后，得到所有裂缝中尺寸最小的裂缝长度为46.88 mm，宽度为21.92 mm。为了方便计算和对比，将其余扫描裂缝的尺寸也处理为长度为46.88 mm，宽度为21.92 mm，裂缝拼接完成后，裂缝的开度可由下列公式得出：

式中：b为裂缝开度大小，单位mm；Z1、Z2分别为裂缝同一基准线上的上表面位置高程和下表面位置高程。在有效应力下，如果裂缝面发生相互接触，则b=0。

1.2 裂缝开度分布特征

在开展数值研究前，对6条裂缝初始状态下的裂缝开度进行分析，结果如图2 和表2 所示。开度场云图（图2）表明，3 种样品裂缝开度分布的空间相关性与岩石类型有关[16,21]。其中，花岗岩裂缝G1 和G2 开度分布呈现出左侧开度小、右侧开度大的特点；致密砂岩裂缝TS1和TS2开度整体呈现出与水平y方向为斜角45°分布特征；页岩裂缝S1 和S2 开度则表现为水平y方向底部开度大。此外，裂缝开度大小概率密度分布统计显示，花岗岩裂缝的开度分布范围相对更宽，平均开度更大，表明裂缝开度分布相对不均匀，开度高低起伏较明显；页岩裂缝和致密砂岩裂缝开度范围较为集中，表明开度大小分布相对集中。表2 给出了3 种岩石裂缝开度分布在x和y方向上的相关长度，表现为花岗岩裂缝开度空间分布相对集中，致密砂岩裂缝开度空间分布范围相对更广，页岩裂缝开度空间分布则主要沿着裂缝y方向分布且y方向上裂缝开度的相关长度与宽度相等。

图2 不同裂缝样品开度场分布云图及相应的裂缝开度概率密度分布Fig.2 Results of fracture aperture distribution and corresponding aperture probability density distribution

表2 不同裂缝样品开度初始参数统计Table 2 Statistical analysis results of apertures of different fracture samples

2 裂缝接触力学及渗流模型

2.1 接触力学模型

学者们通常采用接触力学模型数值研究裂缝在有效应力载荷下的变形规律，该模型遵循HOPKINS[22]、PYRAK-NOLTE等[23]、PETROVITCH等[24]、KLING等[25]提出的方法，且通过了大量实验验证。研究根据接触力学模型分析裂缝应力敏感，在进行裂缝变形分析时，将裂缝面内的微凸起简化为圆柱体，并假设在裂缝的上表面和下表面分别规定一块基准平板，将有效应力均匀施加在该平板上，而在有效应力连续加载过程中，当裂缝渗透率为0时，停止分析，从而得到裂缝在不同有效应力下开度、粗糙度和相关长度等变化规律。根据研究方法，裂缝表面某一圆柱凸起i的变形计算公式为[20,24]：

式中：Wi为粗糙圆柱凸起i的总变形；wii为圆柱凸起i受作用力后的法向变形；wij是圆柱凸起i由邻近凸起j引起的变形。

在圆柱i施加的力为fi，则圆柱i的法向变形wii为：

式中：θ为角度，单位（°）；I1为圆柱i在法向应力作用下的法向变形函数，单位m；r0为圆柱的半径，单位m；r为圆柱i中心的距离，单位m。

其中：

式中：fi为法向应力，单位N；E为弹性模量，单位GPa；v为泊松比。

wij可以用下式表达：

式中：r1和r2是圆柱j沿r方向到圆柱i中心的距离，单位m；θ1和θ2为圆弧和圆柱体j之间的交点；I2为圆柱i由临近凸起j在挤压作用下的横向变形函数，单位m。

其中：

除了2个裂缝面的变形外，凸起长度在法向应力fi的作用下同样被压缩，则：

被压缩后，上下裂缝面距离D和凸起变形Wi之和需等于凸起变形后的长度，则：

式中：D为上下裂缝面距离，单位m。

2.2 单相渗流模型

单相渗流模型将裂缝内的流动空间划分为多个长方体，如图3a所示，每一个长方体内的流量遵循质量守恒定律。在不可压缩假设下，对每个开度单元采用体积守恒定律，通过Hagen–Poiseuille 方程计算裂缝内流体压力分布[26-27]（图3）：

图3 裂缝开度网格示意图Fig.3 Schematic diagram of fracture aperture grid

由于考虑裂缝表面之间的起伏，加入修正因子，将式（9）修正为：

式中：l′ij为i到j的倾斜距离，单位m；li和lj分别为i、j网格长度的一半，单位m。

对每个网格i上的水力传导系数gi通过局部立方定律得到：

式中：gi为对每个网格i上的水力传导系数，单位m3/Pa；bi为i点的开度；Ai为i的横截面积，单位m2；μ为黏度，单位Pa·s，研究中为单相水的流动，黏度设置为1×10-3Pas；lw为网格宽度，单位m。

式（11）为确定每个开度网格上的水力传导系数，然而流体总是从一个网格i流进另一个网格j，因此，不能只用一个网格的水力传导系数计算压力，需要对其做调和平均：

将计算得到的gij代入公式（10），同时将裂缝入口压力设置为10 Pa，出口压力设置为1 Pa，通过求解得到裂缝中流体的压力分布情况，将压力分布带入式（10）得出裂缝内流体流量的分布，并根据压力分布计算裂缝入口流量Q。

通过立方定律[28]计算出水力开度bh：

式中：bh为水力开度，单位m；Q为裂缝入口流量，单位m3/s；∆p为入口与出口之间的压力差，单位Pa；w为裂缝宽度，单位m；L为裂缝长度，单位m。

根据平板模型，裂缝渗透率可以用下式表示[19,28]：

式中：k为裂缝渗透率，单位m2。

3 结果与讨论

3.1 不同有效应力下渗透率变化行为

由于同类裂缝样品在不同有效应力下流场和渗透率演化模拟结果趋势一致，因此，针对3 种裂缝样品各选取其中1块作为代表样品，分析裂缝内流体流场随有效应力变化的响应特征，所有流动方向为从左往右。模拟结果表明（图4）：由于裂缝开度的强非均质性，在有效应力作用下花岗岩和致密砂岩流场演化方向与流体流动方向呈现为垂直交叉特点，而页岩裂缝的流场变化方向与流体流动方向平行，这种由于裂缝开度非均质性和沿流动方向上相关长度的影响对渗透率的控制作用，将在进一步研究中讨论。

图4 有效应力加载过程中不同裂缝内流场演化结果Fig.4 Results of flow field distribution under loading effective stress

研究分析了3 个裂缝样品的渗透率随有效应力变化的模拟结果（图5）：其中，花岗岩裂缝样品G1 的初始渗透率为5.64×103µm2，模拟结束后渗透率降至6.58×10-5µm2，降幅为8 个数量级；致密砂岩裂缝样品TS1 初始渗透率为1.51×103µm2，模拟结束后渗透率降至3.39×10-3µm2，降幅为6 个数量级；页岩裂缝样品S1 的初始渗透率为1.27×103µm2，模拟结束后渗透率降至3.55×101µm2，降幅在2 个数量级。结合图2可以看出：花岗岩和致密砂岩裂缝样品的渗透率随有效应力的增加下降程度更为明显，而页岩裂缝的渗透率变化则相对较小，主要是页岩S1 裂缝下边缘开度较大，在有效应力下一直难以压实，这部分未被压实的开度对整个裂缝的导流能力贡献大，所以渗透率变化较缓慢。

图5 有效应力下不同岩石裂缝渗透率演化结果Fig.5 Evolution of fracture permeability with increasing loading effective stress

3.2 有效应力作用下裂缝特征参数响应

3.2.1 裂缝开度

研究分析了有效应力变化下3 种裂缝样品的裂缝开度演化过程（图6）：裂缝开度演化与裂缝开度的非均质性有关，花岗岩G1 在有效应力不断增加过程中，裂缝右侧开度降幅有限，而裂缝左侧更容易发生闭合；致密砂岩TS1 在有效应力作用下，裂缝内整体开度变化相对均匀；页岩S1的开度变化则表现为从y方向上的一侧往另一侧逐渐压合的过渡行为。结果表明，有效应力加载过程中，裂缝开度变化行为与图4中的流场分布演化呈现出良好的对应关系。

图6 不同有效应力下裂缝开度演化结果Fig.6 Evolution of fracture aperture distribution under different effective stress

3.2.2 空间相关长度

根据裂缝在x和y方向上的实际缝长和缝宽大小，将裂缝开度分布与x和y方向上的空间相关长度做归一化处理，得到归一化后裂缝开度的空间相关长度介于0~1，其中，0 代表没有空间相关长度，1 代表x或y方向上的裂缝开度分布贯穿整个裂缝。分析在初始状态下，只有页岩S1 裂缝开度的相关长度在y方向上贯穿裂缝，且所有裂缝的相关长度在x和y方向上呈现显著差异，说明3 类裂缝开度表现出各向异性[29]。当施加的有效应力小于0.5 MPa 时，所有裂缝的空间相关长度变化较小；随着有效应力超过0.5 MPa，所有裂缝的空间相关长度开始随有效应力的增加呈下降趋势，并最终趋于平缓。分析结果表明（图7），致密砂岩裂缝的相关长度和页岩裂缝沿y方向的相关长度降幅明显，而花岗岩裂缝的空间相关长度则相对缓慢。由表1可知，花岗岩的弹性模量为26 GPa，高于致密砂岩和页岩，且花岗岩裂缝开度分布不均匀，因而有效应力加载过程中，花岗岩裂缝整体上不容易被压缩，其裂缝右侧还有大面积的开度未被压缩（图6），因此，花岗岩的相关长度减小缓慢。

图7 裂缝相关长度随有效应力变化的演化结果Fig.7 Changes of fracture correlation length with increasing effective stress

3.2.3 接触面积

分析有效应力加载过程中裂缝岩样接触面积的变化特征，结果表明（图8）：随着有效应力的增加，所有裂缝接触面积占比逐渐增大且呈现近似线性关系，但不同裂缝接触面积增长速率差异越发明显。在相同的有效应力下，由于花岗岩的强度高和开度分布不均匀，花岗岩裂缝的接触面积远小于致密砂岩和页岩裂缝的接触面积，而致密砂岩TS1裂缝开度的相关长度较小且岩样的强度更低，因此，在有效应力下致密砂岩裂缝接触面积增长的速率最快。

图8 有效应力作用下接触面积与渗透率变化结果Fig.8 Changes of contact area and fracture permeability under effect of effective stress

通过分析接触面积占比与渗透率之间的关系，对比图8a 和图8b 可知：裂缝接触面积的增加与裂缝渗透率的减小并未呈现出良好的正相关性，其中，页岩裂缝有效应力加载后的最终接触面积增加幅度更明显，但其渗透率降幅却相对最小。结合流场分布和裂缝开度演化可知，相比接触面积本身大小，裂缝接触面积在方向上的变化对渗透率的控制更加关键。因此，针对强非均质性开度分布的裂缝，仅仅利用裂缝接触面积变化分析有效应力作用对裂缝渗透率的影响存在局限性，需要考察裂缝接触面积在空间方向上的变化特性，进而综合分析对渗透率的控制作用。针对这种现象，有学者提出了无效渗流区域概念，即裂缝网格中除了接触的这部分区域变成无效流动区域外，还将裂缝网格中局部流量低于裂缝整体平均流量1%的区域也划分为无效流动区域，这2 类区域对裂缝整体流动的贡献忽略不计[16]，反之，对裂缝渗流有贡献的区域定义为有效渗流区域。因此，建议在今后的工作中，采用有效渗流区域作为关键参数表征裂缝导流能力随有效应力变化的影响。

3.3 研究案例与经验公式匹配性

为了验证传统模型是否可以预测通过巴西劈裂法形成的非均质性裂缝的渗透率演化行为，采取常用的ZIMMERMAN 等[19]提出的经典经验公式（简称ZB模型）进行分析，该模型表达如下：

ZB模型得到的计算结果与研究数值计算结果对比，结果表明（图9）：随着有效应力的增大，由于相对花岗岩G1 和页岩S1 裂缝样品，致密砂岩裂缝TS1 开度的非均质性程度更低且相关长度更小，能在有效应力下基本满足均匀接触，因而能够与经典模型较好匹配，因此，除了致密砂岩TS1的数据拟合较好，其余通过经验公式计算的结果与本研究数据结果有较大差异。另外，虽然经验公式计算结构与研究实际结果存在较大误差，但是可以看出两者结果在相对低有效应力下的趋势是一致的，说明在低有效应力作用下，裂缝开度的非均质性对渗透率演化的影响较小。因此，传统的经验公式在分析预测不满足平均化特征的非均质裂缝应力敏感下的渗透率演化时，存在明显的不适性。

图9 研究数据与经验公式计算本文数据结果对比Fig.9 Comparison between results calculated by our method and ZB model

4 结论

本研究中分别利用花岗岩、致密砂岩和页岩样品，采取巴西劈裂法制成真实裂缝岩样，通过三维光学扫描仪获取裂缝上下表面微观结构特征参数后，数值拼接形成整合裂缝，并开展了连续应力加载过程中裂缝流场和渗透率演化行为研究，取得了以下主要结论和认识：

1）采用巴西劈裂法形成的裂缝具有较强的非均质性，且不同种类的岩样在裂缝开度空间分布上表现出不同的非均质性，这种表现会影响裂缝在有效应力下的渗流演化规律。

2）裂缝开度、接触面积和空间相关长度是影响裂缝渗透率的重要参数。针对均质裂缝，有效应力加载过程中，裂缝渗透率降低趋势与该参数线性关系明显，尤其是与接触面积密切相关。随着裂缝非均质性结构的加强，这种线性关系逐渐失效，表现为即使有些裂缝开度没有闭合，但未能形成渗流通道，渗透率仍然降低。

3）传统的经验公式在研究裂缝应力敏感时，裂缝开度分布满足平均化条件，因而拟合程度较好。但在拟合本文利用巴西劈裂法形成的真实裂缝时，拟合程度随有效应力增加而偏差严重。因此，对于存在非均质性特性或裂缝开度并没有整体满足平均化尺度时，采用传统经验公式研究裂缝导流能力变化会存在显著的不适应性。