1.1   数值模型

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

图 1.  PFC^2D模型和平行粘结键

Figure 1.  The PFC ^2D model and parallel-bond bond

下载: 全尺寸图片 幻灯片

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

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

表 1.  平行粘结模型参数

Table 1.  Parameters of the parallel-bond model

细观参数	假定值	标定值
E*/GPa	30	41
k*	2.0	2.6
/GPa	30	20
	2.0	2.6
/MPa	40	45
/MPa	50	85
/（°）	30	30
μ	0.70	0.56

 | Show Table

DownLoad: CSV

图 2.  细观参数对宏观参数的影响

Figure 2.  Effects of meso-parameters on macro-parameters

下载: 全尺寸图片 幻灯片

1.3   参数的标定与验证

本部分采用“试错法”对参数进行标定，其思路为：先确定各细观参数对宏观参数的影响，运用得到的变化规律，反复调试细观参数，直到模拟出的力学响应与试验的力学响应接近为止。通过调整平行粘结有效模量 ${\bar E^ * }$ 和有效模量E^*对弹性模量E进行标定，然后通过调整刚度比k^*和平行粘结刚度比 ${\bar k^ * }$ 对泊松比ν进行标定，并调整粘结强度 $\bar c$ 和抗拉强度 ${\bar \sigma _{\rm{c}}}$ 对峰值强度σ_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.  试验结果与模拟结果的比较

Figure 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.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。裂隙①的倾角α₁=45°，位置不变，端点分别为A和B。裂隙②的倾角α₂在0°～180°变化，端点分别记为C和D。

图 4.  预制双裂隙及相应的岩样模型

Figure 4.  Prefabricated double fractures and the corresponding sample model

下载: 全尺寸图片 幻灯片

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

图 5.  不同围压下裂隙角度α₂对宏观参数的影响

Figure 5.  Effects of fracture angle α₂ on macroscopic parameters under different confining pressures

下载: 全尺寸图片 幻灯片

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

图 6.  不同裂隙角度α₂下最终裂纹扩展形态（σ₃=5 MPa）

Figure 6.  Final crack growth pattern under different fracture angles α₂（σ₃=5 MPa）

下载: 全尺寸图片 幻灯片

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

图 7.  不同裂隙角度α₂下裂纹演化过程（σ₃=5 MPa）

Figure 7.  Crack evolution processes under different fracture angle α₂（σ₃=5 MPa）

下载: 全尺寸图片 幻灯片

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

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

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

2.2   力学特征微观表现

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

图 8.  α₂=0º时不同围压下岩样力链破坏情况

Figure 8.  Force chains failure states of the samples under different confining pressures when α₂=0º

下载: 全尺寸图片 幻灯片

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

图 9.  α₂=0°时裂隙②周围力链

Figure 9.  Force chains around crack ② when α₂=0°

下载: 全尺寸图片 幻灯片

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

图 10.  围压5 MPa作用下的裂隙②周围力链

Figure 10.  Force chains around the fissure ② under the confining pressure of 5 MPa

下载: 全尺寸图片 幻灯片