预制裂隙岩样宏细观力学行为颗粒流数值模拟

王星辰，王志亮，黄佑鹏，贾帅龙

（合肥工业大学土木与水利工程学院，安徽 合肥 230009）

颗粒流程序（PFC）已成为研究岩石力学特性重要工具，在岩土工程领域得到了广泛应用，如模拟岩石力学性能演化、边坡的变形失稳[1]以及桩-土相互作用等。通过PFC标定出的岩石细观参数能较好地表征岩石压缩行为[2]，并可进一步研究含预制裂隙岩石的力学特性，如Yang等[3]采用平行粘结接触模型对预制双裂隙红砂岩进行了单轴压缩模拟，结果表明裂纹扩展模式取决于裂隙的倾角。另有学者[4−7]采用PFC模拟了含预制裂隙的各种岩石的双轴压缩，并对其力学特性和裂纹演化过程和机理进行了分析。此外，王振永[8]运用PFC分析了不同裂隙参数（单裂隙角度和双裂隙延桥倾角）对岩体试样强度和裂纹扩展模式的影响。

综上可知，凭借围压作用下的平行粘结接触模型细观参数标定，来开展预制双裂隙岩样力学行为模拟方面的研究并不多见。因此，本文采用花岗岩常规三轴压缩下的力学特性标定花岗岩的细观参数，通过得到的参数模拟预制双裂隙花岗岩的力学特性，最后分析预制双裂隙花岗岩力学特性及裂纹发展的演化过程，力求得出有参考价值的结论。

1 建立模型及标定细观参数

1.1 数值模型

PFC2D通过离散单元方法模拟单位厚度圆形颗粒的运动和相互作用，从基本颗粒角度来考虑材料的力学特性，反映出裂纹的形成和扩展过程。因此，可采用PFC2D模拟岩石的宏观力学特性[9]。以华山花岗岩为模拟研究对象，其平均密度为2 600 kg/m3。试验采用直径为50 mm、高100 mm的圆柱体试样。为方便与试验结果对比，在PFC2D中建立宽50 mm、高100 mm的二维岩石颗粒模型来进行模拟研究（图1）。该数值模型共随机生成19 791个半径在0.20～0.332 mm之间的圆形颗粒，最大粒径与最小粒径之比为Rmax/Rmin=1.66，颗粒密度为2 600 kg/m3，半径因子的默认值为1.0，设置模型的孔隙率为1%。花岗岩属于硬岩，采用平行粘结模型能够更好地模拟此类材料的压缩力学特性[10]，模型中产生51 459个平行粘结。标定中固定围压σ3为15 MPa，设置加载到轴压降低为峰值强度的70%时停止。

图1 PFC2D模型和平行粘结键Fig.1 The PFC 2D model and parallel-bond bond

1.2 平行粘结模型细观参数分析

压缩作用下岩石的宏观力学响应主要受如下细观参数影响：有效模量E*、刚度比k*、平行粘结有效模量、平行粘结刚度比、抗拉强度、粘结强度、摩擦角和摩擦系数μ。而宏观参数主要包括弹性模量E、峰值强度σp、泊松比ν、内摩擦角φ和黏聚力c。细观参数的初始假定值如表1所示。为简化分析过程，选取对宏观参数影响较大的细观参数进行讨论。细观变形参数包括E*、k*、和。参照表1中的设定值，采用控制变量法，分别设置E*、为10，15，20，25，30，35 GPa，k*、为1，1.5，2，2.5，3，3.5。图2（a）为变形细观参数对宏观参数的影响，可见随着E*和E¯∗的增加，E均呈近似线性增加，且在讨论范围内，其增加幅度分别为58.1%和71.1%。随着k¯∗的增加，E的降低幅度逐渐减小。而k*对E的影响最小。采用相同方法，将、和作为强度细观参数，分别设置、为10～100 MPa，为10°～70°。从图2（b）中可知，σp受和的影响较大，σp的增加幅度分别为106.6%、72.3%，且随的增加σp增加幅度较小并渐趋稳定。

图2 细观参数对宏观参数的影响Fig.2 Effects of meso-parameters on macro-parameters

表1 平行粘结模型参数Table 1 Parameters of the parallel-bond model

1.3 参数的标定与验证

本部分采用“试错法”对参数进行标定，其思路为：先确定各细观参数对宏观参数的影响，运用得到的变化规律，反复调试细观参数，直到模拟出的力学响应与试验的力学响应接近为止。通过调整平行粘结有效模量和有效模量E*对弹性模量E进行标定，然后通过调整刚度比k*和平行粘结刚度比对泊松比ν进行标定，并调整粘结强度和抗拉强度对峰值强度σp进行标定，最后调整围压，运用莫尔-库伦准则来验证内摩擦角φ和黏聚力c。所得到的细观参数列于表1。

图3为采用上述参数模拟和试验得到的应力-应变曲线和破坏模式的对比，其中破坏效果图中红色部分代表裂纹。可以看出，模拟破坏结果中形成1条沿对角线的贯穿裂纹，与试验结果吻合较好。试验和模拟峰值强度分别是272.5 MPa和270.5 MPa，弹性模量分别是58.6 GPa和56.3 GPa，结果比较接近，表明标定的参数具有可靠性。同理，可进行围压为25 MPa和35 MPa的模拟试验，得到高围压下峰值强度。25 MPa围压下试验和模拟的峰值强度分别为320.0 MPa和313.5 MPa，35 MPa围压下试验和模拟的峰值强度为370.0 MPa和375.5 MPa。采用库伦准则来获得黏聚力和内摩擦角：

图3 试验结果与模拟结果的比较Fig.3 Comparison of the experimental and simulated results

式中：c—黏聚力；

φ—内摩擦角。

采用式（1）分别拟合试验和模拟数据。求得试验值c=43.9 MPa，φ=41.7°，模拟值c=41.2 MPa，φ=42.8°，内摩擦角φ和黏聚力误差分别为6.15%和2.64%，表明模拟效果较好。

综上所述，通过标定的细观参数是可行的，能够很好地反映华山花岗岩常规三轴压缩的力学特性，故采用标定的材料参数来模拟带围压的含预制双裂隙华山花岗岩的力学特性与破坏特征。

2 双裂隙岩样双轴压缩模拟

2.1 岩样宏观力学特性与破坏特征

为了研究含缺陷岩样的力学特性，可将不同形式的孔洞、裂隙预制在岩样中，如张科等[11]研究了含圆形灌浆孔洞岩样的力学特性和裂纹扩展模式。Shi等[12]分析了预制裂隙的不均匀性对岩石的影响。本文研究不同围压下（5，15，25 MPa）含双裂隙花岗岩的力学特性，参照文[13]中红砂岩试样模型和预制裂隙样式，设计本文模型和预制裂隙。如图4所示，模型高为100 mm，宽为50 mm，裂隙均为非闭合裂隙，不考虑摩擦阻力[14]，岩桥长度b=10 mm，岩桥倾角β=45°，裂隙长度a=18 mm，宽度d=1 mm。裂隙①的倾角α1=45°，位置不变，端点分别为A和B。裂隙②的倾角α2在0°～180°变化，端点分别记为C和D。

图4 预制双裂隙及相应的岩样模型Fig.4 Prefabricated double fractures and the corresponding sample model

图5为不同围压下峰值强度与弹性模量随裂隙角度α2的变化情况，在各围压下，随着角度α2的增加，σp与E均是先增加后减小，且在α2=90°时达到了峰值。同时，随着围压的增大，各角度下σp与E均逐渐增加。这主要与预制裂隙②在水平方向的投影长度有关。随着角度的增加，2条预制裂隙水平方向的投影长度先减小后增加；当α2到达90°时，可近似看作只有预制裂隙①。易见水平投影长度越长，对岩样的影响越大，岩样的力学性能也就相对较弱。

图5 不同围压下裂隙角度α2对宏观参数的影响Fig.5 Effects of fracture angle α2 on macroscopic parameters under different confining pressures

图6为5 MPa围压条件下预制双裂隙岩样在不同裂隙角度α2下的最终裂纹扩展形态。当α2<90°时，A和D端点裂纹均向试样边缘发展，B和C端点的裂纹搭接并沿主应力方向发展；当α2=90°时，在裂隙②的长度方向上出现大量的裂纹，且有1条裂纹发展到试样的左下角，B端点和预制裂隙②长度方向裂纹搭接形成贯通。当α2=120°时，出现B和D端裂隙贯通，且主要裂隙的宽度增大。当α2>120°时，B和C端点裂隙贯通现象减弱，而A和D端点裂隙贯通，在试样左侧形成以C端点为主要起裂点的裂纹。另一方面，围压对裂隙拓展形态的影响不甚明显，主要体现出围压越大裂隙越宽的特点，此处不再赘述。当α2为0°和180°时，裂隙②中部会产生垂直拉伸裂纹，当α2为120°和150°时，B和D端点裂隙贯通。

图6 不同裂隙角度α2下最终裂纹扩展形态（σ3=5 MPa）Fig.6 Final crack growth pattern under different fracture angles α2（σ3=5 MPa）

为观察不同裂隙角度α2的试样的裂纹发展情况，针对围压5 MPa模拟组开展讨论。限于篇幅，仅对0°、60°、90°、150°进行分析（图7）。试样起裂强度定义为峰值强度σp时裂纹数目的1%所对应的主应力大小[15]。不同加载应力水平下裂纹扩展情况如图7所示，可见各裂隙角度α2下试样均存在裂纹萌生、增长和稳定3个主要阶段：

图7 不同裂隙角度α2下裂纹演化过程（σ3=5 MPa）Fig.7 Crack evolution processes under different fracture angle α2（σ3=5 MPa）

（1）σ1≤0.7σp为裂纹萌生阶段。所有岩样的裂隙尖端开始萌发出微小的翼裂纹，当α2较大时，裂隙附近产生几条明显的拉伸裂纹。特别是在α2=0°时，裂隙②中间段产生了微小的拉伸裂纹。当α2=90°时，裂隙②端点没有萌发裂纹。

（2）0.7σp<σ1≤σp为裂纹增长阶段。裂纹数量迅速增多，先前形成的微裂纹逐渐加宽加长，几条主裂纹开始形成规模，2条预制裂隙通过裂纹接近贯通。当α2=0°时，裂纹主要在岩样上半部分发展；当α2=90°时，裂纹在预制裂隙②在长度方向产生。α2较大时裂纹总面积相对较小，且裂纹宽度也相对较窄，而α2较小时裂隙端点裂纹向岩样边缘拓展现象更明显。

（3）σ1>σp为裂纹形成稳定阶段。各岩样的主要裂纹完全形成，并贯通至岩样边缘，2条预制裂隙通过裂纹完全贯通，且随着加载所有裂纹长度与宽度几乎没有变化。最终的裂纹形态表现为：当α2<90°时，形成左、中、右3条分离的大裂纹；当α2>90°时，形成左、右2条大裂纹。α2=0°时，裂纹区最紧凑，而α2=90°时裂纹区最分散。

2.2 力学特征微观表现

图8为力链在不同围压下的最终破坏情况。图中红色代表拉伸力链，而黑色代表受压力链。当外力大于力链的承载能力时，力链断裂使得颗粒间的作用力消失，从而致使裂纹萌发。在预制裂隙的端点处有受压力链的集中，且随着围压的增加红色部分逐渐减小。在预制裂隙周围力链相对稀薄，低围压下力链的空白区域比较高围压情况下面积相对较大。力链越密集，则承载能力越大，达到破坏而产生的裂纹越多。因此，岩石的峰值强度会随着围压的增大而增大。

图8 α2=0°时不同围压下岩样力链破坏情况Fig.8 Force chains failure states of the samples under different confining pressures when α2=0°

图9为不同围压下预制裂隙②周围力链。可知，在预制裂隙②长度方向上有许多红色的拉伸裂纹，且方向多与预制裂隙平行，表明此处有拉应力的集中，这就是在破坏的过程中裂纹产生并且裂纹方向（力链断裂方向）是沿着主压力方向的原因。同时，在预制裂隙长度方向的力链较为稀疏，在端点处有密集的压缩力链。随围压增大，预制裂隙②长度方向上的红色逐渐变淡，表明拉伸力链逐渐减少。可见，在高围压下，拉应力集中并不明显，不会出现沿着主压力方向产生拉伸裂纹。α2=180°时，预制裂隙②也是如此。

图9 α2=0°时裂隙②周围力链Fig.9 Force chains around crack ② when α2=0°

图10为加载前预制裂隙②周围力链。当α2=90°时，预制裂隙②周围并没有产生力链空白区，这种情况下对岩石产生的影响最小，其力学特性与没有预制裂隙②的时候最为接近。在预制裂隙②的C端点与预制裂隙①的B端点之间并没有拉伸力链，而在B端点与预制裂隙②中部之间区域分布较多，所以在两预制裂隙的贯通上与α2<90°时有所不同。在α2=150°时，两预制裂隙之间力链复杂，B端点与D端点之间以及A端点和D端点之间均有较多拉伸力链，在B、C端点间仅有少量的拉伸裂纹。从裂纹的发展过程可知，在B和C端点之间产生的裂纹是C端点萌发的翼裂纹产生的，最终与C端点贯通，这与其他的两预制裂隙贯通的机制并不一样。所以，α2=150°时，有3条使两预制裂隙贯通的裂纹。

图10 围压5 MPa作用下的裂隙②周围力链Fig.10 Force chains around the fissure ② under the confining pressure of 5 MPa

3 结论

（1）平行粘结接触模型的强度细观参数对花岗岩峰值强度的影响较大，而变形细观参数对弹性模量和泊松比的影响较大。采用“试错法”标定出了一组细观参数，并基于该参数模拟的花岗岩破坏形貌与试验结果吻合较好。

（2）围压作用下预制双裂隙花岗岩的峰值强度和弹性模量均随α2增大呈先增大后减小的趋势，当α2=90°时，两者均达到最大值。裂纹演化过程包括萌生、扩展和稳定3个阶段；α2的不同导致裂隙的贯通方式和贯通裂纹数量均有所差异。

（3）力链承载力随围压增大而增大，表现出宏观的峰值强度越大。拉伸力链对水平裂隙上垂直裂纹的发展和预制裂隙的贯通方式影响较大。围压越高，则水平裂隙周围的拉伸力链越稀疏，产生拉伸裂纹较少。在不同α2下，裂隙端点间及端点与裂隙中部间会产生不同形式的拉伸裂纹。