不同边界对花岗岩三轴试验影响的三维离散元数值研究

王治林，郑明明，夏敏，熊亮，吴祖锐，王凯

（成都理工大学地质灾害防治与地质环境保护国家重点实验室，四川 成都 610059）

0 引言

花岗岩是由不同种类较大晶体颗粒与非晶体颗粒等构成的非连续、非均质且力学行为较为复杂的岩石，且在诸多掘进工程中极为常见，杨达等［1］认为使用气动潜孔锤并且当球齿以30 J冲击功钻进花岗岩时，可有效降低能量损耗与提高碎岩效率。汤凤林等［2］研究了钻进速度与规程参数的优化关系，提出了优化正常钻进的措施。张程等［3］通过室内实验研究了超声波碎岩技术下花岗岩径向位移与内部损伤的关系，并认为超声波碎岩技术能提高在花岗岩等硬岩中的钻进效率，以上学者要么从钻进设备要么从工艺方法等方面研究如何提高钻进效率，而未考虑到花岗岩本身的力学特性，目前许多学者利用数值方法研究岩石的力学特性，而在数值方法中，离散元法对于探究岩石在变形破坏的细观机理方面具有巨大的优势［4-6］。

许多学者［7-9］使用颗粒流软件（PFC3D）分析岩石在不同围压下、加载速率下的力学特性时其边界均使用刚性边界，而常规室内三轴试验中，试样基本都采用橡胶膜或乳胶膜密封，刚性边界不符合实际边界情况，因此，金磊等［10］根据细分概念，提出三维组合墙法模拟三轴边界，其在加载过程中体现出与室内试样相似的横向膨胀现象，能良好地模拟试样的变形特征。Xu等［11］在双轴实验时使用一定厚度的光滑颗粒集合体来模拟柔性薄膜，发现边界条件一定程度上会影响剪切带的形成和形态。Cheung等［12］，Thmos等［13］在双轴实验中对比了柔性边界和刚性边界对试验模拟结果的影响，发现刚性边界下试样变形与室内试样变形形态不符，而柔性边界下试样能够体现出与室内试样相同的横向膨胀现象。Qu等［14］在PFC中提出用膜颗粒微观参数表示柔性膜边界物理属性的方法，并与大量三轴实验结果对比分析，指出使用柔性膜边界是当下较为适合的离散元边界模拟方法。但对于室内试验所用的橡胶膜或乳胶膜，无论是柔性边界还是三维组合墙在力学本质上都难以模拟出其弹塑性特征，而对于弹塑性物体的模拟，有限元软件有着极其显著的优势。

对此，基于最新的有限元-离散元（FDM-EDM）耦合技术，使用PFC分别生成“刚性”、“柔性膜”、“Shell”单元三种边界下的花岗岩三轴试验模型，进而分析“刚性”、“柔性膜”、“Shell”单元3种边界对花岗岩常规三轴模拟试验的影响。

1 花岗岩模型

1.1 花岗岩细观结构

花岗岩中的矿物颗粒是影响花岗岩力学性质的重要因素之一，也是岩石非均质性的重要体现，花岗岩中矿物成分主要为长石、石英、黑云母以及其它矿物等，根据张涛等［15］的研究，本文花岗岩内部矿物成分含量为：长石（57.6%）、石英（30.4%）、黑云母（10.1%）、其它矿物（2.0%）。本文中圆柱形岩样半径为2.1 cm，高度9 cm，岩样中颗粒数量为14799，粒径范围为0.028～0.056 cm，粒径分布服从均匀分布。

1.2 花岗岩模型的生成

颗粒流程序中不能直接生成三维等效晶质模型（Grain-Based-Model），须借助其软件中自带的“块体（Rblock）”单元或者“几何体（geometry）”单元来间接生成三维等效晶质模型。张涛等［15］是通过3个不同阶段来生成三维等效晶质模型：（1）构建试样初始颗粒模型；（2）矿物晶体分组，在初始颗粒模型外部构建一个外接“几何体”单元，将外接几何体单元转化为“块体”单元，根据矿物成分对“块体”单元进行分组，再将分组后的“块体”单元集合再次转化“几何体”单元集合；（3）利用“几何体”单元进行填充，最终生成三维等效晶质模型。

上述过程最终实现了花岗岩三维等效晶质模型的建立，但其生成步骤十分繁琐，且“块体”单元与“几何体”单元之间需要多次转化，降低了计算效率，为此本文在其基础上通过直接使用“Rblock”对初始颗粒模型进行划分，从而简化了生成步骤，提高了计算效率。花岗岩模型生成过程中的3个阶段如图1所示。

阶段一：确定试样尺寸，构建初始颗粒模型，如图1（a）所示。

阶段二：“块体”单元分组。在初始颗粒模型外部构建一个外接“几何体单元”，将外接“几何体”单元再次转化为“块体”单元，对“块体”单元进行分组，如图1（b）、（c）所示。

阶段三：生成三维等效晶质模型。通过判断“块体”单元与初始模型颗粒位置的关系对初始模型颗粒进行分组，从而生成三维等效晶质模型，如图1（d）、（e）所示。

图1 花岗岩模型生成示意Fig.1 Schematic diagram of granite model generation

1.3 模型微观参数确定

颗粒流程序中自带多种接触模型，其中平行粘结模型（Parallel Bond Model）适合用来模拟岩石颗粒之间的接触［16-19］，因此本文中岩石颗粒之间的接触均采用平行黏结模型。此外本文将岩石内部颗粒之间的接触称为晶体内接触和晶体间接触。根据Ding等［20］以往的花岗岩常规三轴离散元模拟研究以及目前效果较好的“试错法”，本文中微观接触参数如表1所示。

表1 平行粘结模型微观参数Table 1 Micro-parameters of the parallel bonding model

2 三轴模拟试验建立

边界的建立和伺服是三轴模拟试验的关键步骤，前面已介绍过如何建立花岗岩模型的方法，因此这里主要介绍本文中花岗岩三轴压缩试验的3种边界。

2.1 “刚性”边界

PFC可以直接生成常规三轴模拟试验中的刚性边界，刚性边界属于刚体，受力后不会产生任何变形，但通过刚性伺服函数可以获得径向方向上的一个平移自由度，细观上表现为刚性边界顶点的径向位移，宏观上表现为刚性边界的整体放大与缩小［21］，由于毛海涛等［21］对刚性伺服机制的深入解释，本文中便不再详细解释刚性伺服机制，“刚性”边界如图2（a）所示。

2.2 “柔性膜”边界

蒋成龙等［22］详细介绍了柔性膜的生成方法和伺服理论。本文利用其方法生成常规三轴实验所用的柔性膜边界，基于柔性膜伺服理论，使柔性膜边界下的花岗岩模型达到目标围压。柔性膜颗粒在空间上具有6个自由度，可以任意变形，但本质上属于离散元边界，Qu等［14］经过理论推导、试验验证找到了柔性膜颗粒间的接触参数与柔性膜宏观参数（弹性模量、泊松比）的关系，基此，本文使用Qu等［14］发现的关系通过柔性膜宏观参数确定柔性膜颗粒的接触参数，本文中柔性膜边界的弹性模量和泊松比分别为7.8 MPa、0.46。“柔性膜”边界如图2（b）所示。

图2 3种不同边界条件示意图Fig.2 Schematic diagram of three different boundary conditions

2.3 “Shell”单元边界

室内试验所用的橡皮膜在力学本质上属于弹塑性物体，而在模拟弹塑性物体时，有限元软件比离散元有巨大优势［23］，有限元法中的“Shell”单元由带有3个节点的三角面组成，其中每个节点都具有6个自由度（3个平移自由度、3个旋转自由度），如图3所示：其变形分为膜变形以及梁变形，对于各向同性的“Shell”单元，其力学控制方程符合式（1）与式（2），其中E和v为“Shell”单元的弹性模量和泊松比，本文中E=7.8 MPa，v=0.46。同时其显著优点是：“Shell”单元边界的围压可以通过Structure Shell Apply命令精准施加且围压在加载过程中保持恒定不变，因此引入FLAC中的“Shell”单元作为三轴模拟试验的边界。本文中“Shell”单元边界由2304个三角形“Shell”单元均匀组成，如图2（c）所示。

图3 Shell单元示意图Fig.3 Example of Shell element

2.4 花岗岩三轴加载

为在离散元中研究边界对花岗岩三轴压缩试验的影响，设计9种岩样在3种不同边界下的模拟实验，共计27组模拟试验，结果见表2。

表2 模拟岩样标号及试验参数Table 2 The number of samples and experimental parameters

3 不同边界下模拟结果的对比分析

3.1 3种边界下的岩样变形差异

为对比不同情况下的岩样变形破坏情况，同时鉴于文章篇幅，本文仅列出了X61岩样在3种边界下的变形破坏情况，如图4（a）、（b）、（c）所示：岩样在柔性膜边界和“Shell”边界下都发生了明显的横向膨胀现象，这也与图4（d）中甘霖等［24］所做的花岗岩三轴试验中花岗岩发生横向膨胀的现象基本吻合，然而岩样在刚性边界下却没有发生横向膨胀的现象，这是一定程度上由于刚性伺服机制，为了保持伺服围压稳定在目标区间内，整个圆柱形刚性边界进行整体同步运动，因此在围压相同、轴应变（ε1）相同的情况下，相对于Shell边界而言，刚性边界下岩样的径向应变（ε3）大于“Shell”单元边界下岩样的径向应变（ε3），而根据MARTIN［25］所提出的体积应变公式（3），刚性边界下的岩样体积应变（εv）应大于Shell单元边界下的岩样体积应变（εv）。

图4 试样变形形态对比Fig.4 Comparison of deformation of different samples

为对上述猜想进行验证，列出模拟岩样的体积应变曲线，鉴于文章篇幅，本文仅列出X61岩样在3种边界情况下的体积应变曲线，如图5所示。从图5中可以看出，在同一轴应变时，X61岩样在刚性边界下的体积应变大于在Shell单元边界下的体积应变，从而验证了上诉猜想的合理性。

图5 3种边界下X61试样的应力-体积应变曲线Fig.5 Stress vs volume strain curves of X61 specimens under three boundaries

3.2 岩样应力-应变曲线

三轴试验中岩样的应力-应变曲线是反映不同边界对花岗岩三轴试验影响的主要形式，图6为27组岩样三轴模拟试验的应力-应变曲线图。

从图6中可知：不同边界下岩样的三轴压缩曲线的应力-应变曲线的变化趋势大体上是一致的，曲线在峰值应力前近似呈线性变化，岩样进入屈服阶段后，曲线斜率逐渐变小，峰后阶段岩样强度均出现明显的下降。相同加载速率下，岩样在3种边界条件下的峰值应力随围压的增加而增长，其中岩样在柔性膜边界下的峰值应力增加幅度最为剧烈。相同围压条件下，岩样在刚性边界下的峰值强度受加载速度的增加影响最大，在柔性膜边界下的峰值强度受加载速率的增加影响最小。同一加载条件下，岩样在柔性膜边界下的残余强度大于在其它2种边界下的残余强度，这主要存在2个原因：一是由于根据柔性膜伺服理论计算出的柔性膜表面积略大于实际的柔性膜边界表面积，造成岩样在柔性膜边界下所受的实际围压大于设定的目标围压，尤其是目标围压越大的情况下，岩样所受的实际围压与目标围压差距越大；二是由柔性膜边界中的膜颗粒本身具有质量，对岩样形成了端部约束［14］。因此相对比其他两种边界下的岩样，岩样在柔性膜边界下的峰值应力易受围压影响且残余强度普遍大于其他两种边界下岩样的残余强度。同时发现在同一加载条件下，岩样在刚性边界下的峰值应力大于在Shell边界下的峰值应力，这一定程度上是由于刚性伺服函数的连续运行之间存在时间间隔以及刚性伺服函数确定围压值时所存在的误差造成的［21］。

3.3 不同边界对岩样力学强度影响

根据图6中的数据绘制出3种边界下岩样峰值应力、弹性模量随加载速率变化的分布曲线，如图7、8所示。

图6 3种边界下不同围压及加载速率下试样应力-应变曲线Fig.6 Axial and radial stress vs strain curves under three boundaries at different confining pressure and loading velocity

从图7中可以看出3种边界下，随着围压的增加，不同加载速率下的岩样峰值应力均大幅增加，各力学参数之间与围压之间满足线性正相关关系。Shell边界下的岩样在加载速率为0.5 mm/s，围压为6 MPa时峰值应力为222 MPa，当围压增加到12、24 MPa时，峰值应力增加到289、415 MPa，增长了约30%、86%。而当加载速率增加到1.0、1.5 mm/s时，峰值应力增加到239、253 MPa，增长了约7.6%、13.9%。柔性膜边界下的岩样在速率为0.5 mm/s，围压为6 MPa时，其峰值应力为232 MPa，当围压增加到12 、24 MPa时，柔性膜边界下岩样的峰值应力增加到352、575 MPa，增长了约51%、147.8%，而当加载速率增加到1.0、1.5 mm/s时，峰值应力增加到242、254 MPa，增长了约4.3%、9.5%。刚性边界下的岩样在速率为0.5 mm/s，围压为6 MPa时，其峰值应力为270 MPa，当围压增加到12、24 MPa时，刚性边界下岩样的峰值应力增加到331、447 MPa，增长了约22.6%、65.6%，而当加载速率增加到1.0、1.5 mm/s时，峰值应力增加到293、318 MPa，增长了约8.5%、17.8%。对比发现，对于3种边界下岩样的峰值应力而言，由围压效应造成的影响大于由加载效应造成的影响，其中柔性膜边界下岩样的峰值应力受所受围压影响最大，刚性边界下岩样的峰值应力受加载速率影响最大。

图7 三种边界下岩样峰值应力与加载速率关系Fig.7 Relationship between peak stress and loading rate of rock samples under three boundaries

从图8中可以看出，3种边界下岩样的弹性模量都随着围压或者加载速率的增加而增加，但增长幅度不大，其中Shell边界下岩样的弹性模量受加载速率影响较大，柔性膜边界和刚性边界下岩样的弹性模量受围压影响较大。

图8 三种边界下岩样弹性模量与加载速率关系Fig.8 Relationship between elastic modulus and loading rate of rock samples under three boundaries

3.4 岩样力链分布

荷载在试样通过颗粒间的接触力链传递，分析接触力链可从细观角度上对比分析不同边界对岩样三轴压缩实验的影响，因此绘制出3种边界下岩样在轴应变（7%）的接触力链图，如图9所示，其中线条粗细代表接触力链强度大小。

从图9可知，岩样破坏后在不同边界下的力链强度、力链位置分布均存在差异。对于同一加载条件下的岩样，刚性边界下岩样在边界接触处存在一条强度远大于其他接触力链的接触力链，这说明此时岩样在刚性边界下已经发生了应力集中现象，而Shell边界下岩样则未出现此情况。而对于Shell边界下岩样，其接触力链强度较为均匀，较少出现个别接触力链强度远大于其他接触力链强度的现象。

4 结论

本文使用颗粒流程序研究了刚性、柔性膜、“Shell”3种边界对模拟花岗岩常规三轴试验的影响。针对3种边界下花岗岩的宏观破坏形态和应力应变关系、力链分布等模拟结果进行了分析。通过本文研究，得出结论如下：

（1）边界影响岩样的加载效应和围压效应。岩样在不同边界下，其受到加载速率和围压的影响程度不同，从上面的三轴模拟实验中可以看出，相同围压下，刚性边界下岩样的峰值应力易受加载速率的影响，相同加载速率下，柔性膜边界下岩样的峰值应力受围压影响最大。

（2）边界对岩样的破坏变形有显著影响。从图4、图9中可知，柔性膜和Shell边界下的岩样会发生横向不均匀膨胀破坏，与花岗岩室内实验现象相符，而刚性边界下的岩样整体上呈现整体压缩破坏，不太符合花岗岩室内实验现象。

（3）刚性边界和柔性膜边界在三轴加载过程中无法对岩样施加精准且稳定的围压。其中刚性边界由于刚性伺服函数的连续运行之间存在时间间隔以及刚性伺服函数确定实际围压值时所存在的细小误差，导致岩样在刚性边界下无法获得稳定的围压。而柔性膜边界在三轴加载过程中根据柔性膜伺服理论计算出的柔性膜表面积略大于实际的柔性膜边界表面积，导致岩样在柔性膜边界下所受的实际围压大于目标围压，所以岩样在柔性膜边界下的峰值应力易受目标围压影响且残余强度普遍大于其他2种边界下岩样的残余强度。

（4）刚性边界对岩样在三轴加载过程中应力集中现象的产生具有影响。从图9中可知，岩样在刚性边界下发生应力集中的概率较大，这是由于对于Shell边界和柔性膜边界而言，二者在力学性质上具有6个自由度，而刚性边界由于其伺服机制，其在力学性质只具有一个自由度，因此相对于其他2种边界，岩样在刚性边界下发生应力集中的概率较大。

图9 3种边界下岩样接触力链图（轴应变为7%）Fig.9 Contact force chain of rock samples under three boundaries（ε1=7%）

综合对比刚性边界、柔性膜边界和Shell边界，发现在Shell边界下，岩样围压施加简单精准且在加载过程中其所受围压恒定不变，应力集中发生概率较小，因此在花岗岩离散元常规三轴模拟时，Shell边界是合适的选择。