随常规油气开采难度的增加，非常规油气资源的勘探开发已成为常规油气资源的有效补充和世界能源发展的必然战略选择^[1]。砂砾岩储层作为一种重要的非常规油气储层，受沉积过程的影响，储层地质构造复杂，且储层内分布了大量不规则砾石颗粒。砾石与基质力学性质的差异以及空间分布极度不规则性导致砂砾岩的力学特征复杂，表现出较强的岩石力学非均质性^[2-4]。这一非均质性使得砂砾岩的宏观破裂模式和内在破裂机理明显不同于其他岩石^[5]。砂砾岩的力学性质直接影响压裂的改造效果，可能会导致压裂过程中的裂缝延伸几何形态、支撑剂运移等产生不可预知的结果，因此弄清砂砾岩的岩石力学特征能为砂砾岩油藏压裂设计提供理论依据。

已有大量学者对岩石的力学特征进行了研究，其研究方法主要包括物理实验方法和数值模拟方法。在物理实验方面，单轴及三轴压缩实验常用于岩石力学性质测试，被广泛应用于岩石裂缝的破裂形态及其影响因素研究^[6-7]。针对大理岩、砂砾岩、页岩等多种岩石也有大量的力学实验研究，研究了围压、岩性差异、尺寸效应、层理等多方面因素对岩石破坏的影响。研究表明，岩石裂缝的破裂形态主要受到围压、层理以及岩石力学参数等多方面的影响，且不同围压对抗压强度、弹性模量等力学性质的影响相当显著^[8-13]。从能量角度研究岩石破坏过程也是常用的方法之一，已有研究表明岩石的能量耗散和吸收与其单轴抗压强度和围压等因素呈良好的线性关系^[14-15]。

由于仪器、取样等因素的限制，物理实验往往很难模拟真实的储层情况，且物理实验具有成本高、可重复性较差等缺点，因此建立等效的数值模型来开展岩石力学特征分析是常用的研究方法之一。YAN等^[16]学者利用RFPA2D研究了节理方向、岩石尺寸对岩石破坏特征的影响。李兆霖等^[17]借助ABAQUS软件对岩石进行了模拟研究，研究表明岩石的宏观破裂具有方向性，且受主应力方向的影响显著。然而，从宏观尺度上虽然可以把储层岩石看作连续介质，但从微观尺度上来看，岩石颗粒之间存在孔隙，且颗粒与颗粒之间是离散的。对于砂砾岩这种非均质性较强的储层，非连续介质法和粒子离散元法在研究岩石的细观断裂扩展过程和微观裂纹上面具有更大优势^[18-19]。刘京铄等^[20]利用离散单元法研究了单轴压缩条件下不同倾角对岩石破裂特征的影响。颗粒离散元法能模拟颗粒之间的运动、变形、破坏特征^[21-22]。CAO等^[23]通过颗粒离散元法研究了循环荷载下砂岩非线性破坏，发现砂岩在塑性和弹性阶段的损伤破坏有明显不同，塑性阶段损伤破坏更迅速。CHO等^[24]采用颗粒离散元法研究了脆性岩石在剪切过程中宏观破坏带的发展过程。CHONG等^[25]研究了岩石在巴西劈裂实验条件下的压缩破坏过程，讨论了裂隙力学性质、岩石脆性及方向等因素对破坏特征的影响。

综上，国内外学者针对砂砾岩岩石力学特征开展了大量研究，但从细观尺度开展砂砾岩细观力学特征分析研究较少。鉴于此，本文以准噶尔盆地百口泉组砂砾岩储层为研究对象，针对砂砾岩的非均质性特征，首先提出了随机不规则多边形砾石生成方法，其次基于颗粒离散元法建立了砂砾岩力学数值模型，从细观角度开展了含单砾石类型变形破坏特征分析，在此基础上对混合砾石变形破坏机理进行了分析研究。

1.   砂砾岩数值模型构建

1.1   砾石随机生成数学方法

砾石多呈棱角或次棱角状，主要为四边形、五边形、六边形、七边形和八边形。故本文采用不规则多边形来近似替代几何形状更为复杂的真实砾石。如图 1所示，假设多边形砾石轮廓共有N条边，砾石可以描述成具有N个顶点的极坐标系统，每个顶点的坐标为(θ_il_i)。为了生成任意多边形砾石，Δθ_i和l_i作为两个自由变量，并遵循随机均匀分布。随机变量l_i可表示为：

图  1  随机多边形砾石的几何形状示意

Figure  1.  Geometry diagram of random polygonal gravel

下载: 全尺寸图片 幻灯片

式中：l为半径l_i的平均值，单位mm；η₁是随机系数，介于0~1；Δl是半径l_i的增量，单位mm。

随机变量Δθ_i可由两个相邻顶点的极角增量来表示:

式中：η₂是随机系数，介于0~1；N是多边形的边(顶点)个数；δ是变量参数，一般取值为0.3。

为生成闭合多边形，需对Δθ_i进行如下调整：

某一顶点的极角可由下式求得：

在笛卡尔坐标系中，顶点坐标可由下式表示：

式中：(X₀，Y₀)是砾石的中心坐标值。

已有研究表明，砂砾岩中的许多砾石呈扁平多边形状，部分形似于椭圆形，因此，扁平状的程度可由长短轴比γ表征：

式中：L_a是砾石的长轴大小，单位mm；L_b是砾石的短轴大小，单位mm。砂砾岩砾石长短轴比的取值范围一般为1.2~2.3。

为生成扁平状的多边形砾石，可以采用扩大长轴的长度或者缩小短轴的长度来实现。此外，通过控制砾石长轴与x轴的夹角可以实现对砾石长轴方位的控制。

1.2   数值模型校正

为验证数值模型的准确性，选取了三叠系百口泉组砂砾岩储层具有代表性的三轴压缩(围压40 MPa)物理实验结果作为数值模拟验证对象。以此针对砂砾岩开展了砾石矿物成分分析，同时进行纳米压痕实验，发现该储层中主要存在4种典型砾石，分别为砾石A、砾石B、砾石C及砾石D，如图 2所示。砾石A主要呈黑白色，主要矿物成分为凝灰岩、霏细岩及流纹岩，宏观弹性模量平均为11.49 GPa，单轴抗压强度平均为75.73 MPa，断裂韧性平均为1.10 MPa·m^0.5；砾石B主要呈灰黑色，主要矿物成分为灰岩、霏细岩，宏观弹性模量平均为9.82 GPa，单轴抗压强度平均为63.23 MPa，断裂韧性平均为0.78 MPa·m^0.5；砾石C主要呈亮黑色，主要矿物成分为石英岩、凝灰岩，宏观弹性模量平均为17.01 GPa，单轴抗压强度平均为117.16 MPa，断裂韧性平均为1.85 MPa·m^0.5；砾石D主要呈深红色，主要矿物成分为石英岩、流纹岩及泥岩，宏观弹性模量平均为15.52 GPa，单轴抗压强度平均为105.97 MPa，断裂韧性平均为1.59 MPa·m^0.5。此外，砂砾岩基质的宏观弹性模量平均为9.08 GPa，单轴抗压强度平均为56.70 MPa，断裂韧性平均为0.97 MPa·m^0.5。对此，本文基于颗粒离散元法建立了含不规则砾石的非均质模型，同时采用软粘结模型(Soft Bond Model)分别针对基质、砾石A、砾石B、砾石C、砾石D进行细观力学参数标定。物理实验中的4种砾石比例约为A∶B∶C∶D = 1∶1∶1∶1，因此在数值模型中同样按照此比例生成砂砾岩数值模型。砂砾岩在40 MPa围压条件下的三轴压缩物理实验和数值模拟结果见图 3，两种模拟均产生了两条倾斜的主破坏带以及一些零散的宏观裂缝，应力应变曲线结果表明高围压下砂砾岩具有较强的延性特征，且峰值强度较为接近，误差率为6.59%。由于实际岩心与砂砾岩数值模型的砾石含量及分布存在一定差别，两种模拟中的应变特征存在一定差别，但其对裂缝的形态特征影响较小。因此，本文建立的数值模型能表征砂砾岩岩石力学破坏特征。

图  2  准噶尔盆地三叠系百口泉组砂砾岩照片

Figure  2.  Photograph of glutenite from Triassic Baikouquan Formation, Junggar Basin

下载: 全尺寸图片 幻灯片

图  3  准噶尔盆地三叠系百口泉组砂砾岩三轴(围压40 MPa)压缩模拟结果

Figure  3.  Simulation results of glutenite triaxial compression (confining pressure 40 MPa) in Triassic Baikouquan Formation, Junggar Basin

下载: 全尺寸图片 幻灯片

2.   含单砾石类型变形破坏特征

为了揭示砾石对砂砾岩破坏的影响机理，分别开展了不同围压(10、20、30、40 MPa)下含4种典型单砾石(砾石A、B、C、D)的三轴压缩数值模拟，分析岩石破坏过程中微裂缝及能量演化曲线的变化特征。

2.1   含单砾石A模型

当砂砾岩模型仅含较低强度砾石A时，不同围压下三轴压缩数值模拟得到的破坏裂缝形态如图 4a所示。当围压为10 MPa时，裂缝呈双“V”字形破坏，随围压的增加，裂缝逐渐呈现为高角度的单剪切破坏带。砾石A对裂缝扩展的屏蔽作用较小，围压增加产生的剪切微裂缝数量更多，且剪切微裂缝主要分布在宏观破坏带的中心区域。图 5a的微裂缝演化曲线表明，围压增大，岩石破坏时机越晚，剪切微裂缝数量呈线性增加。当围压从10 MPa增加到40 MPa时，破坏后可产生多达3 028条剪切微裂缝，剪切微裂缝比例由16.7%增加到33.5%。图 5a的能量曲线表明，围压越大，初始应变能越大。随着围压增加，岩石破坏所需要的应变能增加，应变能降低后保持的水平更高。当围压为40 MPa时，最大应变能达到1 981.5 kJ/m³。同时围压越大，岩石破坏后产生的滑移能呈增加的趋势，围压40 MPa的最大滑移能为204.9 kJ/m³。随着围压增加，岩石破坏所需的能量和应变都增加，这与传统认识一致。

图  4  不同围压下砂砾岩三轴压缩破坏形态

Figure  4.  Triaxial compression failure patterns of glutenite under different confining pressures

下载: 全尺寸图片 幻灯片

图  5  不同围压微裂缝数量及能量曲线的演化过程

Figure  5.  Evolution of number and energy curves of microcracks under different confining pressures

下载: 全尺寸图片 幻灯片

2.2   含单砾石B模型

当砂砾岩模型仅含低强度砾石B时，不同围压下三轴压缩数值模拟破坏后的裂缝形态如图 4b所示。结果表明砂砾岩模型只含低强度砾石B时破坏形态相对单一。在不同围压下砂砾岩的主裂缝破坏带均呈双“V”字形。主破坏带附近产生的微裂缝随围压增加而增加。砂砾岩的宏观主破坏带特征表明不同围压下低强度砾石B对裂缝扩展几乎无影响。围压增加导致岩石微裂缝起裂时机延迟，且破坏后的剪切微裂缝数量明显增加，剪切微裂缝占比由15.9%增加到32.7%，但张性微裂缝数量变化较小(图 5b)。由于低强度砾石B抵抗变形破坏的能力更弱，其破坏所需能量更少，但受基质的影响，导致其滑移能略有增加。因此，与只含砾石A模型相比，只含砾石B砂砾岩模型最大应变能有小幅度降低，但滑移能增加。

2.3   含单砾石C模型

当仅含高强度砾石C时，不同围压下三轴压缩破坏后的裂缝形态如图 4c所示。数值模拟结果表明，仅含高强度砾石C的砂砾岩模型在不同围压下均形成了贯穿试样的斜剪主裂缝，倾斜角度约45°，微裂缝的分布随围压增加分散性增强。当围压为10 MPa时，存在一定的穿砾现象，但砾石的破坏程度较轻，砾石内以张性破坏为主。随围压增加，穿砾现象越来越明显，同时高围压条件下容易在高强度砾石周围诱导产生剪切缝。仅含单砾石C时三轴压缩作用下微裂缝数量及能量演化曲线如图 5c所示，试样破坏产生的张性微裂缝数量随围压增加而降低。当围压为10 MPa时，张性微裂缝数量最多，达到7 266条。试样破坏产生的剪切微裂缝数量随围压增加而增加，剪切微裂缝占比由15.2%上升到30.9%。能量演化曲线表明，高围压下砂砾岩存在明显的二次破坏特征，模型中含高强度砾石时会提升砂砾岩的延性。

2.4   含单砾石D模型

图 4d是仅含较高强度砾石D的砂砾岩模型在不同围压下的三轴压缩数值模拟结果。结果表明，含较高强度砾石D模型在4种围压下的破坏均出现了贯穿试样的斜剪主裂缝带，其破坏特征与只含砾石C模型相似。但随围压增加，含较高强度砾石D模型出现了明显的宏观分支裂缝。图 5d表明张性微裂缝数量随围压增加而减少，剪切微裂缝数量随围压增加而增加，剪切微裂缝占比由15.8%增加到31.1%。与只含砾石C的砂砾岩模型相比，含砾石D的砂砾岩模型在围压40 MPa条件下破坏所需的应变能差异不大，滑移能略微下降。

2.5   不同单砾石类型对比分析

不同围压条件下砂砾岩三轴压缩模拟结果对比表明，随围压增大(0~40 MPa)，张性微裂缝数量总体呈现先增加后减小的趋势，在围压为10 MPa附近时达到最大，而剪切微裂缝数量持续显著增加(图 6a-b)。产生该现象的主要原因是当围压达到10 MPa后，围压的增大普遍会使得砂砾岩储层的脆性降低，而岩石塑性提高，主裂缝中心区域会产生更多的剪切微裂缝代替张性微裂缝，导致张性微裂缝数量呈一定程度的下降。随围压增加，峰值应变能及滑移能都呈近线性增加，其中应变能显著增加。随围压增加，砂砾岩的抗压强度均明显增加，近似呈线性增加(图 6c)。三轴压缩下的弹性模量明显大于单轴弹性模量，但不同围压下的三轴弹性模量相差不大(图 6d)。

图  6  不同围压砂砾岩轴向压缩模拟结果对比

Figure  6.  Comparison of axial compression simulation results of glutenite under different confining pressures

下载: 全尺寸图片 幻灯片

3.   混合砾石细观变形破坏机理

上文针对仅含单砾石砂砾岩开展了三轴压缩破坏特征分析，但实际砂砾岩中含有多种砾石，为明确混合砾石条件下砂砾岩细观变形破坏机理，下文讨论了不同围压作用下含4种典型砾石砂砾岩三轴压缩细观变形破坏机理。

3.1   10 MPa围压

选取了具有代表性的4种砾石比例的砂砾岩数值模型进行10 MPa围压下的三轴压缩数值模拟(图 7a)。模拟结果表明，当模型中不含砾石C时，宏观破坏呈现双“V”字形。当模型中含有砾石C时，宏观破坏形态复杂，裂缝的破坏路径明显受到砾石C的影响，呈现为绕开砾石C破坏的形式，即含有高强度砾石对岩石的破坏形态影响较大。在不同砾石比例条件下，岩石破坏的时机较为接近(图 8a)。当砾石组合中含有砾石C时，最终产生的微裂缝数量更多。且当高强度砾石C及D含量更高时(1∶0∶1∶1)，岩石破坏所需要的应变能最大。

图  7  不同围压下不同砾石比例砂砾岩三轴压缩破坏形态

Figure  7.  Failure patterns of glutenite with different gravel ratios under different confining pressures in triaxial compression

下载: 全尺寸图片 幻灯片

图  8  不同围压下不同砾石比例砂砾岩三轴压缩微裂缝数量及能量曲线

Figure  8.  Number and energy curves of glutenite with different gravel ratios under different confining pressures in triaxial compression

下载: 全尺寸图片 幻灯片

3.2   20 MPa围压

与围压10 MPa模型相比，围压20 MPa模型破坏后产生的微裂缝更为集中，剪切微裂缝更多(图 7b)。除砾石比例为A∶B∶C∶D=1∶1∶1∶1时破坏形态呈单一的主裂缝之外，其余砂砾岩的破坏形态较为接近，均呈双“V”字形破坏。当砾石比例为A∶B∶C∶D=0∶1∶1∶0时，即当砂砾岩含高强度砾石C的比例较高时，产生的微裂缝数量更多(图 8b)。这是由于受砾石C高强度的影响，压缩破坏时穿砾石所需能量较高，其会导致微裂缝分散和数量增多的现象。与10 MPa围压下的模型相比，20 MPa模型中的砂砾岩破坏开始显现二次破坏特征，但不明显。

3.3   30 MPa围压

与围压20 MPa模型相比，30 MPa围压下试样的破坏形态变化不大，但是微裂缝数量明显增加，其中剪切微裂缝大幅增加(图 7c)。不同砾石比例条件下，岩石发生破坏的时机相近，但产生的微裂缝数量差别较大，当砾石比例为A∶B∶C∶D=0∶1∶1∶0时，微裂缝数量更多，且峰值应变后的能量曲线波动更大，表现出明显的二次破坏特征(图 8c)。

3.4   40 MPa围压

与围压为10~30 MPa模型相比，40 MPa围压下破坏后的微裂缝数量和剪切微裂缝进一步增加，且宏观裂缝形态复杂，广泛分布于试样四周(图 7d)。与10~30 MPa围压模型相比，40 MPa围压时不同砾石比例条件下产生的微裂缝数量差别较大，岩石破坏后的能量及应变曲线差别也较大，分析认为在高围压(30 MPa和40 MPa)条件下，高强度砾石的存在会使砂砾岩的塑性明显提高，导致峰后破坏特征明显不同(图 8d)。

3.5   讨论

不同砾石组合在不同围压条件下砂砾岩三轴压缩模拟结果对比如图 9所示。在混合砾石模式下，抗压强度随围压增加而显著增加，但不同围压条件下的三轴弹性模量差别不大。不同砾石组合条件下，三轴抗压强度及弹性模量存在一定差别，但较小。当砾石比例为A∶B∶C∶D = 1∶0∶1∶1时，对应的砂砾岩抗压强度及弹性模量更大，即不含低强度砾石有利于增加岩石的强度和抵抗变形的能力。在三轴压缩作用下，张性微裂缝数量与围压没有明显的变化关系，但随着围压增加，剪切微裂缝数量持续大幅增加。不同砾石比例条件下，砂砾岩压缩产生的微裂缝数量差别较大。总体上，当含高强度砾石占比较大时，产生的张性及剪切微裂缝数量更多。

图  9  不同砾石比例砂砾岩轴向压缩模拟结果对比

Figure  9.  Comparison of axial compression simulation results of gravel conglomerate with different gravel ratios

下载: 全尺寸图片 幻灯片

4.   结论和建议

以我国西部准噶尔盆地百口泉组砂砾岩为背景，通过颗粒离散元法开展了三轴压缩数值模拟。结果表明：

(1) 当模型仅含单砾石时，剪切微裂缝随围压增加而增加，围压大小与剪切破坏程度直接相关。峰值应变能及滑移能随围压增加均呈线性增加，其中应变能增加显著。围压对砂砾岩抗压强度的影响十分显著，抗压强度随围压增加而增加。高围压下，含高强度砾石的砂砾岩在达峰值强度后易出现二次破坏的现象，表现出明显的塑性及延性特征。

(2) 不同砾石组合条件下，弹性模量存在明显差别，当砂砾岩不含低强度砾石(砾石比例为A∶B∶C∶D = 1∶0∶1∶1)时，即高强度砾石占比更多时，对应的砂砾岩弹性模量更大，抵抗变形能力越强。不同砾石比例条件下，砂砾岩压缩产生的微裂缝数量差别较大。总体上，当砂砾岩含高强度砾石时，产生的张性及剪切微裂缝数量更多。

(3) 砂砾岩变形破坏过程极为复杂，其宏观破坏带的形成发育过程很大程度受控于内部细观结构，且受围压、砾石类型(力学强度)等的影响较大。因此，在实际工程压裂施工时，在压裂前选井选层方面建议选择在中低围压、砾石强度较低且砾石含量较少的储层进行压裂，以此增加砂砾岩储层的渗流通道，从而提高采收率。