随着能源需求的不断增加，砾岩油气藏受到越来越多的关注，特别是在新疆准噶尔盆地玛湖凹陷，其砾岩储层发育突出，现已发现石油地质储量 20× 10⁸ t ^[1]。砾岩储层普遍埋藏较深，具有物性较差、孔隙度和渗透率较低、岩性复杂、强非均质性等特点，自然产能低^[2-5]，需要通过“井工厂+体积压裂”组合工艺形成复杂的裂缝网络以提高经济效益^[6-8]。然而，由于砾岩储层中存在随机嵌入的砾石，载荷作用下砾石表面出现应力集中现象，显著影响局部力学响应^[9]。高强度砾石引起的应力干扰使原始裂缝的扩展发生偏转，形成更复杂的裂缝^[10-14]。

对于常规砂岩储层，裂缝一般沿着最大主应力方向扩展，呈双翼状；对于砾岩储层，强非均质性导致强度和应力分布不均匀，裂缝形成和扩展存在高度不确定性。邹雨时等对人造砾岩进行了室内真三轴压裂实验，发现界面性质对裂缝的生长路径影响很大，观察到裂缝与砾石之间的多种相互作用模式，包括止裂、穿砾、绕砾和穿砾且绕砾等^[15-17]。LIU 等通过真三轴压裂实验结合计算机断层扫描和声发射 （AE）监测砾岩中裂缝的形成，发现在高应力差环境下，裂缝极大可能穿过砾石，形成较为简单平整的裂缝面^[18]。李连崇等^[13]、LI 等^[19]、张子麟等^[20] 通过岩石破裂过程分析系统（rock failure process analysis sys‐ tem，RFPA）进行水力裂缝扩展模拟，发现裂缝遇到砾石后发生局部偏转，导致裂缝路径曲折复杂；但在较高的水平应力差影响下，裂缝沿着最大主应力的方向扩展，在宏观上形成了较为平直的裂缝。 ZHANG 等利用颗粒流程序（ particle flow code， PFC）研究了地应力差异、砾石强度和粒度分布的影响，得出了相似的结论^[21]。余东合等^[23]、RUI 等^[10] 通过自编渗流-应力-损伤（flow-stress-damage，FSD）模型研究了基质渗透率、压裂液黏度和排量对裂缝遇到砾石时的扩展形态和裂缝总长度的影响。

中外学者对于砾岩储层中压裂裂缝的扩展行为机理研究主要集中在砾石宏观分布的影响上，对于单个砾石研究较少，导致裂缝遇单个砾石的扩展规律认识不清。而实际上，裂缝在砾岩储层中扩展是由多个遇单砾行为相互协同的结果，首先需要对裂缝遇单个砾石的扩展机理进行分析。此外，水力裂缝遇砾扩展行为涉及渗流-应力-损伤的多场耦合，裂缝形态是三者耦合作用的最终结果，裂缝遇砾扩展过程中，渗流场、应力场以及裂缝形态的变化规律仍不清楚。因此，在细观尺度上将砾岩视为砾石、基质和胶结面组成的三相复合材料，基于全局黏聚力单元建立二维渗流-应力-损伤全耦合模型，模拟多场耦合下裂缝遇单个砾石的扩展行为；通过对多点、线进行应力、流体压力监测，分析裂缝扩展过程中，渗流场、应力场以及裂缝形态的变化规律，研究裂缝遇砾扩展机理，旨在为砾岩储层压裂设计提供一定的理论依据。

1 基于全局黏聚力单元的裂缝扩展模型

1.1 数学模型

建立了在流体压力作用下岩石介质中的裂缝扩展模型（图1）。y轴为一竖直井筒，左侧为储层压裂裂缝，右侧为压裂裂缝模拟单元。裂缝单元内充满高压流体，在流体压力作用下，裂缝单元不同部位有着不同的响应：蓝色网格为已经破裂的区域（图1b），裂缝顶部的网格点A为即将破裂区（图1a）。在黏聚力单元的裂缝扩展模型中，用于计算的理论主要包括两部分：孔隙与裂缝中的流体流动方程以及黏聚力单元的牵引分离准则。

对于黏度为 μ 的牛顿流体，裂缝内的压裂液流动可以近似为平行板之间的流动，裂缝内流体流速 q_f可以由压力梯度和局部裂缝宽度确定^[23]：

裂缝内流体满足质量守恒，其连续方程为^[24]：

裂缝中的流体滤失表示为单元中流体的法向流动^[25]，即图1b中黏聚力单元的上表面和下表面的流体滤失，其计算方程为：

根据达西定律，多孔介质中的流体扩散满足^[26]：

根据有效应力原理引入 biot 系数 α，总应力 σ_i，j 可由有效应力σ′ _i，j表示为：

基于虚功原理，岩石多孔介质中变形平衡方程可写为^[26]：

单位体积 dV 中流体总质量的变化率与穿过单位面积 dS 的流体质量流速相等，可得出某时刻 t 流体质量连续性方程^[27]：

在平面应变条件下，岩石发生Ⅰ型断裂的临界断裂能G_IC可由岩石Ⅰ型断裂韧性K_IC、弹性模量E以及泊松比v确定^[28]：

采用二次名义应力准则来判断裂缝起裂^[27-28]：

单元发生损伤后，应力分量可由损伤系数 D 获得^[29]：

加载过程单元开始损坏，D 从 0 线性增加至 1。黏聚力单元完全破坏发生的位移δ_f的表达式为^[27-28]：

单位未失效前，其应力-应变曲线关系满足^[27-28]：

d_z可表示为^[27-28]：

单元失效后的软化机制可表示为^[27-28]：

单元完全破坏发生的位移δ_f满足^[27-28]：

单元整体损伤系数D的表达式为^[27]：

δ_m可表示为：

当单元发生剪切滑移，利用最大拉应力和摩尔库伦理论进行判定^[27]：

1.2 全局黏聚力单元模型

基于 Abaqus 平台通过 Fortran 编程语言编写插件，实现在任意 2 个相邻四边形块体单元之间嵌入黏聚力单元，得到流-固耦合的全局黏聚力单元模型 （图2）。通过该模型可以将有限元法和离散元法耦合起来，模拟砾岩中流体驱动裂缝的萌生、扩展、勾通和分叉。模型中多个块体单元之间通过黏聚力单元连接，通过设定不同参数来表征砾岩的非均质性，包括高强度砾石、低强度基质以及胶结弱面等。块体单元之间的断裂是通过黏聚力单元的破坏失效实现的，为裂缝扩展提供了潜在路径。

1.3 模型验证

MA 等采用边长为 30 cm 的立方体人造砾岩进行室内真三轴压裂实验，研究了砾石分布对压裂裂缝扩展的影响^[16]。为了验证模型的可靠性，依据 MA 等给出的实验参数^[16]，建立 2D 砾岩模型，通过确定数值模拟中裂缝遇到砾石的扩展形态和扩展规律是否与物理实验一致，来验证模型的准确性。建立模型尺寸为0.3 m×0.3 m，包含48 269个实体单元和 98 188 个黏聚力单元；颗粒半径为 5～20 mm，共 223块砾石随机嵌入试样中。模拟输入的最小主应力为 5 MPa，主应力差为 10 MPa；流体注入排量为 20 mL/min；基质、砾石和胶结面物理参数如表1 所示。

对比裂缝在砾岩中扩展形态的数值模拟与室内实验结果（图3）发现，数值模拟得到的结果无论在宏观还是细观都与物理实验有较好的一致性。宏观上，数值模拟和室内实验的砾岩试样压裂裂缝形态主要受应力状态和砾石的控制。裂缝遇到砾石后发生局部偏转，导致裂缝路径曲折复杂。但在较高的水平应力差影响下，裂缝沿着最大主应力的方向扩展，在宏观上形成了较为平直的裂缝。细观上，数值模拟中裂缝遇到砾石出现了止裂、穿砾、绕砾、裂缝分叉等扩展行为，这与MA等的物理实验结果一致。此外，MA 等还观察到压裂后的砾岩岩体中存在多个裂缝带以及裂缝内部脱落的砾石颗粒^[16]，这些现象都能够通过模型模拟出来。模拟中水力裂缝扩展形态与室内实验结果具有很好的一致性。基于全局黏聚力单元的裂缝扩展模型不仅可以准确模拟砾岩岩体中水力裂缝的扩展，还可以捕捉扩展过程中裂缝与砾石之间的微观相互作用。综上，利用全局黏聚力单元的裂缝扩展模型来研究水力裂缝遇砾行为是可行的。

2 多场耦合作用下裂缝遇砾扩展行为

2.1 算例分析

用于模拟裂缝遇砾扩展行为的压裂区域是长为 0.24 m、宽为 0.12 m 的长方形（图4b）。为避免边界效应，根据 Saint-Venant 原理，将压裂区域嵌入在边长为1.2 m的二维网格正方形模型中心区域，即满足模型尺寸为压裂区域的5倍大小（图4a）。由于模型尺寸远大于目标压裂区域，水力压裂对固定边界应力分布的干扰很小。模型由 193 776 个块体单元和 386 992 个黏聚力单元组成，最小单元尺寸为 0.000 4 m，最大为 0.02 m。模型设定的岩石基本物理参数、施加应力条件以及压裂液性质如表2所示。注入点位于压裂区域中心，压裂液注入速度为0.072 m³ /min；施加水平方向最大主应力为 21 MPa，最小主应力为15 MPa。

压裂液注入时长 t 分别为 0.432 2，0.543 5， 0.680 4 s的裂缝形态、流体流速场、流体压力场以及应力场分布（图5）表明裂缝遇砾发生了穿砾行为。根据t为0.432 2 s的流体流速场、流体压力场以及应力场可以看到：当裂缝扩展至砾石，由于砾石具有低渗透高强度的特性，裂缝尖端附近的应力场与渗流场明显受砾石遮挡的影响。为了更好地研究裂缝扩展至砾石过程中，应力场和渗流场的相互耦合作用关系，分析两者变化规律，对监测点P在不同注入时长的水平方向应力S₁₁、竖直方向应力S₂₂以及流体压力p_p进行提取；以监测点P₀作为无砾石对照组，如图4c所示。

图6为不同注入时长下监测线 L和 L₀的应力和流体压力分布，以及监测点 P 和 P₀处的变化率。需要注意的是，根据弹性力学约定此处拉应力为正，压应力为负。可以看到，随着裂缝不断向砾石靠近，P 和 P₀处的水平方向应力 S₁₁与竖直方向应力 S₂₂都在减小，而流体压力 p_p在不断增加。通过与无砾石的对照组比较发现：在高弹性模量、低渗透砾石的阻挡作用下，裂缝前端区域的水平与竖直方向应力都有减小的趋势，而流体压力不断增加。出现这种现象是由于压裂液在多孔介质渗透过程中，由于砾石的低渗透性，引起基质岩体中流体流速明显受阻，裂缝预扩展区域发生憋压现象，流体压力不断增加。根据有效应力原理，在原有地应力场不变的前提下，流体压力的增加导致有效应力的减小。由图6d 可以得到，流体压力 p_p 相较于对照组平均增加了 9.8%， S₁₁ 平均减小了 1.85%，S₂₂ 平均减小了 1.73%。相较于原地应力场，现地应力场存在一定的应力释放现象，裂缝尖端压裂环境发生了改变^[24]。

提取了作用在砾石表面即胶结面的应力分布，其中包括S₁₁，S₂₂以及τ_s（剪应力）。受到裂缝向两侧的挤压作用以及砾石憋压效应，由应力叠加后作用在砾石表面的应力分布（图7）可以看到，S₁₁ 在砾石表面圆心角 φ 为 0°～22°和 338°～360°（图4c）范围内变化较大；S₂₂在φ为0°～60°和300°～360°范围内都有较为明显的应力释放现象。除此之外还发现，随着裂缝前端不断地向砾石扩展，砾石表面逐渐出现了剪应力集中的现象（图7c）。在砾石表面 φ 为 11°和 349°的位置，剪应力集中的现象最为明显，以注入时长为 0.43 s 为例，作用在 φ 为 11°和 349°处的剪应力最大，达到 0.5 MPa。当裂缝逐渐靠近砾石，剪应力集中现象越发明显，增加胶结面弱面破坏的可能性。

为了研究裂缝扩展过程中应力场对裂缝破坏损伤作用机理，对裂缝遇砾行为进行局部放大分析（图8）。裂缝尖端的静水压力较高，在砾石入口处产生了分支裂缝，但裂缝仍穿透砾石，主裂缝不偏斜（图8a）。在裂缝扩展至砾石表面，即注入时长为0.543 5 s时，裂缝张开导致上下单元发生位移，对胶结面形成剪应力载荷（图8b），胶结面单元发生剪切损伤，造成No.285928单元破坏并失效（图8c）；由于施加在胶结面上的剪切能量不足，胶结面上相邻单元No.140146未能破坏，裂缝无法继续（图8d）；砾石由于受到的张应力达到其极限抗张强度，发生张性破坏（图8e）。综上，裂缝遇砾的扩展行为是一个涉及渗流-应力-裂缝损伤的多场耦合作用复杂现象。

2.2 渗透率对裂缝遇砾行为的影响

岩石渗透率不仅影响岩石孔隙中的流动特性，还影响孔隙压力变化和压裂液漏失系数。相同排量下，低渗透油藏压裂液漏失系数低于高渗透油藏，导致缝内静水压力高，压裂液对裂缝尖端的载荷作用增大。为研究渗透率对裂缝遇砾行为的影响，建立了不同砾石渗透率的遇砾模型。为了控制变量，维持基质渗透率K_m = 0.1 mD不变，通过改变砾石与基质渗透率比（K_g /K_m = 0.1，0.2，0.4，0.6，0.8，1.0）的方式，改变砾石渗透率。

提取了注入时长为0.432 2 s时刻作用在监测线 L 上的应力与流体压力，并以无砾石实验组监测线 L₀作为对照组，得到不同 K_g /K_m下监测点 P 点处应力和流体压力的变化率（图9）。从图9d 中可以看到，当砾石渗透率最小时，即 K_g /K_m = 0.1 时，砾石对流体流动的阻挡作用最明显，裂缝在砾石表面的憋压作用显著。此时，相较于无砾石对照组，在相同位置处的流体压力 p_p 提高了 13%，而 S₁₁ 与 S₂₂ 分别减小了 10.1% 和 12.8%。这是由于砾石相较于基质的低渗透性，导致裂缝尖端存在憋压现象，相较于无砾石对照组，流体压力明显升高。根据有效应力原理，流体压力的增加导致了该处有效应力减小的趋势。

从不同渗透率条件下作用在砾石胶结面的剪应力分布（图10a）可以看到，受砾石渗透率影响，砾石表面的剪应力分布差异明显，在砾石表面φ为33°和 327°处，出现剪应力集中现象。当K_g /K_m = 1.0时，即砾石渗透率最高，砾石表面应力集中处 τ_s 为 0.48 MPa；随着砾石渗透率降低，当 K_g /K_m = 0.1 时，裂缝尖端的憋压作用最显著，相同位置处的 τ_s 为 0.61 MPa，提升幅度达28.3%。

由不同渗透率下的裂缝遇砾行为（图10b）可见，砾石渗透率对裂缝遇到砾石的扩展行为影响很大。随着砾石渗透率的降低，裂缝遇砾行为由穿砾向穿砾-绕砾转变，即砾石渗透率越低，裂缝偏斜的可能性越大。这与 RUI的研究结论一致，不同的是 RUI在研究中改变了基质的渗透性^[10]。

2.3 弹性模量对裂缝遇砾行为的影响

在其他参数不变的条件下，维持基质弹性模量 E_m = 36 GPa 不变，改变砾石弹性模量 E_g，满足 E_g /E_m = 1.0，1.2，1.4，1.6，1.8，2.0，2.2，2.4，2.6，2.8，3.0。

由监测线L和L₀上的应力与流体压力分布以及变化率（图11）可以看到，受砾石高弹性模量特性影响，裂缝尖端靠近砾石位置（P 点）的 S₁₁相较于无砾石对照组减小了 0.65%，并且随着砾石弹性模量的增加，S₁₁ 减小幅度增大，最大减小了 5.8%（图11d）。和 S₁₁相比，S₂₂受砾石弹性模量差异的影响较小，这说明砾石对主应力的屏蔽作用与裂缝的扩展方向有关。由于砾石本身相较于基质具有低渗透特性，裂缝尖端区域同样出现了憋压现象，造成 P 点的流体压力明显大于无砾石对照组。虽然实验中没有对砾石渗透性进行改变，但可以看到不同的实验组中监测点P的流体压力在不同弹性模量条件下仍有一定差异（图11c）：随着砾石弹性模量的增加，P 点的流体压力有增大的趋势，由4%增至6.8%（图11d）。造成这种现象可能是应力-渗流耦合作用造成的：由于砾石的弹性模量增加，其抗形变能力也有一定提升； 高压流体作为裂缝扩展的能量供给，在砾石前端聚集而出现流体压力增大。

从图12a可以看到，不同弹性模量的情况下，作用在砾石表面的剪应力分布差异明显，并在砾石表面 φ 为 33°和 327°处，出现剪应力集中现象。当 E_g /E_m = 1.0时，砾石弹性模量最低，胶结面出现剪应力集中，τ_s 达 0.51 MPa；随着砾石弹性模量增至 E_g /E_m = 3.0，相同位置处的 τ_s为 0.82 MPa，提升幅度达60.7%。

由不同弹性模量的裂缝遇砾行为（图12b）可以看到，砾石的弹性模量差异对裂缝遇砾扩展行为影响很大。随着弹性模量的增加，裂缝遇砾行为由穿砾向穿砾-绕砾再向绕砾转变，即砾石弹性模量越高，裂缝偏斜的可能性越大。裂缝在弹性模量较高的砾石附近易发生转向、扭结或产生分支，导致主裂缝局部弯曲，这与许多学者的研究结论一致^[10，16-18]。

3 结论

裂缝遇砾扩展行为是一个涉及渗流-应力-损伤多场耦合作用的复杂现象。基于全局黏聚力单元的裂缝扩展模型可以准确模拟砾岩岩体中水力裂缝的宏观扩展，捕捉扩展过程中裂缝与砾石之间的细观相互作用。

砾石的渗透性和弹性模量对裂缝遇砾行为的作用机理不同，主要分为由于砾石低渗透特性造成憋压现象以及由于砾石高弹性模量特性造成应力屏蔽现象两方面。

砾石低渗透特性导致多孔介质中流体传导受阻，造成憋压现象引起裂缝前端流体压力增加，有效应力减小，同时导致砾石表面出现剪应力集中。在相同条件下，砾石渗透率降低，胶结面剪应力增加，裂缝遇砾行为将由穿砾向穿砾-绕砾转变。

砾石高弹性模量特性导致裂缝前端的应力场受到屏蔽，引起最大水平主应力大幅度减小；由于压力-渗流耦合作用使得流体压力有较小幅度的增加，砾石表面出现剪应力集中。在相同条件下，砾石的弹性模量增大，一方面提高砾石的抗变形能力，另一方面加剧了作用在胶结弱面的剪应力集中，造成裂缝遇砾行为将由穿砾向穿砾-绕砾再向绕砾转变。

符号解释

c_l ——滤失系数，m/（Pa·s）；

d_z ——单元初始厚度，m；

D——单元整体损伤系数，无量纲；

E——弹性模量，MPa；

E_g——砾石弹性模量，MPa；

E_m——基质弹性模量，GPa；

f——单元体积力，N/m³；

G_c——混合模式断裂能，N/m；

G_IC——岩石发生Ⅰ型断裂的临界断裂能，N/m；

I——单位矩阵；

K——渗透率，mD；

K_g——砾石渗透率，mD；

K_m——基质渗透率，mD；

K_IC——I型断裂韧性，Pa·m^0.5 ；

L——监测线长度，m；

n_w——岩石孔隙度；

p_b——下表面流体压力，Pa；

p_f ——裂缝内流体压力，Pa；

∇p_f ——压力梯度，Pa/m；

p_m——裂缝附近地层孔隙压力，Pa；

p_n ′ ——裂缝尖端周向应力，Pa；

p_p——流体压力，Pa；

p_t ——上表面流体压力，Pa；

p_w——地层孔隙压力，Pa；

q_f ——裂缝内流体流速，m² /s；

q_l ——流体滤失，m² /s；

q_m——多孔介质中的流体流速，m² /s；

r——极坐标半径，m；

S——面力作用的单元面积，m²；

S₁₁——水平方向应力，Pa；

S₂₂——竖直方向应力，Pa；

t——注入时长，s；

t_l ——液体作用在单元面积上的法向力，Pa；

T_max——抗张或抗剪强度，Pa；

v——泊松比，无量纲；

v_w——穿过单位面积的流体流速，m² /s；

V——单位体积，m³；

w——裂缝开度，m；

x_i ——裂缝扩展长度，m；

α——biot系数，无量纲；

β——极坐标角，（°）；

δ——单元位移，m；

δ₀——单元初始损伤发生的位移，m；

δ_f ——单元完全破坏发生的位移，m；

δ_m——单元有效位移，m；

δ⁰ _m——单元初始损伤发生的有效位移，m；

δ^f _m——单元完全破坏发生的有效位移，m；

δ_n，δ_s和δ_t ——n，s和t方向的位移分量，m；

δ_v——虚速度场，m/s；

δ_ε——虚应变率，s^-1；

φ——圆心角，（°）；

η——摩擦系数；

μ——流体动力黏度，Pa·s；

ρ_w——流体密度，kg/m³；

σ——应力分量，Pa；

$\stackrel{-}{\sigma }$——单位未损伤时所对应的极限应力，Pa；

σ_i，j ——作用于多孔介质内的总应力，Pa；

σ′ _i，j ——作用于多孔介质内的有效应力，Pa；

σ_n，σ_s和 σ_t ——作用在单元上在 n，s和 t方向的牵引力， Pa；

σ⁰ _n，σ⁰ _s 和σ⁰ _t——单元在n，s和t方向的极限强度，Pa；

σ_βx——沿β-x方向应力，Pa；

σ_h——最小主应力，MPa；

σ_H——最大主应力，MPa；

σ_v——垂向应力，MPa；

τ_max——单元接触面极限剪应力，Pa；

τ′ ——摩擦剪应力，Pa；

τ_s ——作用在胶结面上的剪应力，Pa；

$⟨⟩$——麦考利括号，表示法向压缩状态不会引发损坏。

参考文献