0.   引言

自然界中的岩石因受到长期地质构造作用和复杂物理地质环境的影响,表现为非均质、不连续等复杂特性^[1]。在岩土工程的计算分析中,常采用离散元方法分析岩体的复杂力学特性,随着近年来离散元方法的快速发展,其中的颗粒流方法已经实现了对脆性岩石材料微裂纹行为的模拟分析^[2]。通过对岩石细观及微观的观察,发现多数岩石都是由不规则的矿物颗粒紧密结合在一起组成,矿物颗粒的形状对岩石宏观力学行为存在一定的影响。所以,提出一种能够符合不同类型岩石力学特性的颗粒流模型,并实现构建大尺度颗粒模型的数值方法,对深入研究岩体的细观破裂特征,分析深部洞室围岩宏观破裂尺度效应规律等方面具有重要意义。

近年来,学者们主要从裂隙扩展的角度开展了岩石细观破裂特征的颗粒流研究,如胡波等^[3]采用颗粒流程序对试样的三轴压缩及蠕变进行了数值模拟,研究了的含裂隙硬岩长期蠕变变形行为及蠕变模型；邓树新等^[4]利用试验设计法针对硬质岩体颗粒离散元数值计算中的细观参数标定问题进行了研究；Bahrani等^[5]提出SDA法估算微裂纹岩石的强度以及含不同种类岩石试样的极限强度裂纹密度；Shi等^[6]采用离散元方法,研究了矿物配比和微裂隙分布对岩石宏观力学性能的影响。丛宇等^[7]基于颗粒流原理并结合大理岩室内加、卸荷试验确定了适用于岩石类材料的细观参数,提出了裂纹分布多样化的本质原因。

在颗粒流分析岩石宏观破裂过程中,岩体颗粒流模型具有一定的尺度效应,如Stavrou等^[8]采用离散元方法进行岩体的尺度分析,提出了块石强度与其体积和条件的关系；孙超等^[9]通过不同高宽比均质岩样在不同端部约束条件下的离散元单轴压缩数值试验,探讨了无微裂隙均质岩石尺寸效应产生的微观机制；谢璨等^[10]开展了不同试件尺寸及裂隙倾角的光弹性单轴压缩试验,分析了岩石裂隙扩展失稳尺寸效应的影响；Peng等^[11]对粒子效应进行了数值研究,分析了粒度配比对材料的变形和强度行为、微裂纹扩展过程和诱发微裂纹形态的影响；陈庆发等^[12]利用块体化程度指标,从岩体空间结构角度,揭示了岩体几何与力学尺寸效应的关联性；梁昌玉等^[13]通过对中低应变率加载条件下岩石尺寸效应能量特征进行物理试验,探究了引起强度尺寸效应的本质动因。

在颗粒流分析中,如何获得准确的岩石拉压比细观参数,如何构建大尺度颗粒流模型等方面,仍然是科研工作者面临的主要问题。本文基于颗粒流方法,提出可变半径比例的Clump模型构建方法,计算分析模型细观参数及细观结构特征对模拟岩石试件拉、压特性的影响,采用室内试验方法分析不同类型裂隙花岗岩细观破裂特征,并构建适用于花岗岩力学特性的Clump颗粒流结构,提出构建大尺度颗粒流模型方法,分析深部围岩局部破裂的宏观尺度效应。研究成果对深入研究硬岩的岩爆、水压致裂及锚固失效机理等具有重要参考意义。

1.   PFC模拟岩石力学行为的方法

颗粒流程序PFC (particle flow code)是通过离散单元方法来模拟圆形颗粒的运动及其相互作用,从基本颗粒结构角度考虑基本力学特性,在分析材料从线弹性阶段至断裂破坏的大变形过程中,能直观反映裂纹的形成、扩展及贯通过程^[2],因此,颗粒流方法已经成为模拟岩石类材料宏细观力学行为的主要方法。

PFC中提供了两种基本的黏结模型：接触黏结模型（CBM）和平行黏结模型（PBM）。在接触黏结模型中,黏结剪切破坏后,只要颗粒保持接触,剪切强度就跌落至由摩擦决定的恒定值,接触刚度仍存在。然而,一旦黏结张拉破坏,不论颗粒是否接触,黏结的拉伸强度为零,如图1所示。因此,在接触黏结模型中,如果保持颗粒接触,黏结破坏可能不会显著影响宏观刚度,这对于岩石来说是不可能的。

图  1  接触黏结模型（CBM）法向与切向力学响应

Figure  1.  Normal and tangential mechanical responses of CBM

下载: 全尺寸图片 幻灯片

平行黏结模型由线性模型与黏结模型组成,当荷载过大导致黏结破坏时,平行黏结模型就等价于线性模型^[14],如图2所示,k_n为法向刚度,k_s为切向刚度,g_s为黏结激活间隙,${\overline{\sigma }}_{\text{c}}$为黏结抗拉强度,$\overline{c}$为黏结内聚力,$\overline{\varphi }$为摩擦角,μ为摩擦系数。平行黏结模型的刚度由接触刚度和黏结刚度共同决定,黏结断裂会立即导致刚度降低,这不仅会影响所黏结颗粒的刚度,还会影响整个模型的宏观刚度,从这个意义上讲,平行黏结模型更加适用于岩石类材料。因此,本文选用平行黏结模型来模拟花岗岩的宏细观破裂特征。

图  2  平行黏结模型(PBM)示意图

Figure  2.  Diagrams of parallel bond model (PBM)

下载: 全尺寸图片 幻灯片

2.   Clump模型及参数

2.1   Clump模型

Cho等^[15]首先提出了Clump平行黏结模型（CPBM）,随机将多个颗粒黏结组合形成一个Clump结构体,团聚的颗粒可以像单个具有不规则形状但作为刚体移动的颗粒一样运行。Clump中颗粒的旋转速度是固定的,只有Clump结构体本身才能有旋转速度,如图3所示。Potyondy等^[16]的研究表明,使用普通PFC模型会与试验实测值有较大的偏差,而采用Clump建模方法,可以限制过多的颗粒旋转,能够较好地模拟颗粒的加载和破坏过程。

图  3  颗粒模型中的Ball与Clump旋转机制

Figure  3.  Ball and Clump rotation mechanism in particle model

下载: 全尺寸图片 幻灯片

2.2   Clump模型构建方法

目前建立Clump模型的方法主要有3种：Cho等^[15]通过限制一个Clump中颗粒的数量,形成不规则颗粒模型；Yoon等^[17]通过定义各个Clump外轮廓圆形的位置,将中心位于同一圆形中的颗粒组成一个Clump；Shi等^[18]先建立模型外轮廓,通过一系列pebble相互重叠以形成Clump。

本文基于颗粒流的基本原理,提出一种可变半径比例的Clump结构,该结构的建模方法如下：

（1）在给定范围的墙体（Wall）内生成半径比例为R_min∶R_max=1∶x（x >1）的大颗粒（Ball_max）,如图4（a）所示；然后通过编写Fish语言将这些颗粒的坐标点及半径等信息存储到数据库中,再将这些大颗粒（Ball_max）删除。

图  4  Clump结构建模方法示意图

Figure  4.  Diagram of modeling method for Clump structure

下载: 全尺寸图片 幻灯片

（2）在给定范围的墙体（wall）内重新生成r_min∶r_max=1∶y（r <R,y >1）的小颗粒（Ball_min）,如图4（b）所示。

（3）从数据库中调出大颗粒的坐标及半径数据,并对墙体（Wall）内的小颗粒进行标记,如图4（c）所示,即将每个大颗粒圆形范围内的小颗粒分别划分成组（Group）,而未被标记分组的颗粒,通过Fish语言将该颗粒所有接触的未分组的颗粒分成一组,所建立的这些组数就是后续建立的Clump数。

（4）通过Fish语言将每一组小颗粒转化成Clump模板,然后按照这些模板建立Clump结构,最后将这些小颗粒（Ball_min）删除,最终的结构模型如图4（d）所示。

2.3   Clump模型参数

基于本文提出的可变半径比例Clump模型,通过颗粒流程序PFC^2D,进行单轴压缩与巴西圆盘劈裂试验数值模拟,探究Clump模型与Ball模型的粒径尺寸、细观参数对模拟岩石试件宏观力学特性的影响,找到能够更好地匹配室内试验参数的规律。其中单轴压缩试验颗粒流试件尺寸为长100 mm、宽50 mm,巴西圆盘劈裂试验颗粒流试件尺寸为直径50 mm,模型内小颗粒（Ball_min）的r_min=0.26mm,r_min∶r_max=1∶1.5。

（1）粒径尺寸的影响

以表1中的模型力学参数为基础,保持模型细观力学参数及大颗粒的粒径比R_min∶R_max=1∶1.5不变,改变Clump与Ball的最小粒径R_min,变化区间为0.3～1.0 mm,研究Clump模型与Ball模型粒径尺寸变化对模拟岩石试件宏观力学特性的影响。

表  1  模型力学参数

Table  1.  Mechanical parameters of model

参数	Ball模型/ Clump模型
颗粒刚度比kn/ks	2.0
颗粒摩擦系数µ	0.2
黏结模量E*/GPa	2.0
黏结刚度比$k_{\text{n}}^{*}/k_{\text{s}}^{*}$	2.0
黏结抗拉强度${\sigma }_{\text{b}}$/MPa	36
黏结内聚力$c_{\text{b}}$/MPa	27
摩擦角ƒ/(°)	32

下载: 导出CSV 

| 显示表格

拉压强度变化如图5（a）所示,由于Clump结构限制了其内颗粒的旋转,导致抗压强度比Ball模型大,随着Clump半径的增大,Clump结构内的颗粒增多,Clump结构体之间的边界摩擦效应随之增强,因此Clump模型试件的单轴抗压强度与抗拉强度均不断增加,而Ball模型试件的单轴抗压强度随着Ball半径的增大无明显变化,抗拉强度呈增大趋势。

图  5  Clump/Ball粒径对试件宏观力学特性的影响

Figure  5.  Effects of Clump/Ball particle size on macro-mechanical properties of specimens

下载: 全尺寸图片 幻灯片

张拉裂纹比例及单轴抗压强度（UCS）与抗拉强度（TS）之比如图5（b）所示,通过计算可知,Clump模型试件的张拉裂纹比例随着半径的增大而逐渐减小,这是因为张拉裂纹沿轴向扩展时,Clump结构特性导致应力路径发生改变,转为剪切破坏。粒径尺寸变化对Ball模型试件微裂纹扩展路径的影响并不明显,张拉裂纹比例趋于平稳,基本保持不变。Clump模型试件的拉压比随半径的增大呈显著增大趋势,而Ball模型试件的拉压比受颗粒半径的影响很小。

（2）模型力学参数的影响

选取半径比例R_min∶R_max=1∶1.5的Ball模型与Clump模型,以表1中模型细观力学参数为基础,逐渐调整其中一个参数,保持其它参数不变,在一定参数变化范围内,研究黏结模量E^*、黏结刚度比$k_{\text{n}}^{*}/k_{\text{s}}^{*}$、黏结内聚力$c_{\text{b}}$与黏结抗拉强度${\sigma }_{\text{b}}$对单轴抗压强度与抗拉强度之比UCS/TS的影响。

图6为Ball模型模拟结果,随着各参数的增大,拉压比逐渐趋于稳定。通过对Ball模型参数变化的计算可知,各参数对Ball模型的拉压比有一定的影响,其中黏结模量E^*的影响最大。但是,在参数达到一定值时,对拉压比影响减小,尤其是在黏结刚度比$k_{\text{n}}^{*}/k_{\text{s}}^{*}$达到1.6后影响较小。

图  6  细观参数对Ball模型UCS/TS的影响

Figure  6.  Effects of meso parameters on UCS/TS of Ball model

下载: 全尺寸图片 幻灯片

图7为Clump模型模拟结果,黏结内聚力$c_{\text{b}}$与黏结抗拉强度${\sigma }_{\text{b}}$对Clump模型的拉压比影响较明显,其中黏结抗拉强度${\sigma }_{\text{b}}$的影响最大,随着黏结抗拉强度${\sigma }_{\text{b}}$的增大,拉压比呈减小趋势。黏结刚度比$k_{\text{n}}^{*}/k_{\text{s}}^{*}$对Clump模型的拉压比无较大影响,对模型的破裂模式有一定影响。

图  7  细观参数对Clump模型UCS/TS的影响

Figure  7.  Effects of meso parameters on UCS/TS of Clump model

下载: 全尺寸图片 幻灯片

通过对比参数分析可知,各参数在一定程度上对岩石的宏观力学特性有影响,但是通过对比Ball模型与Clump模型可知,在一定参数范围内,拉压比随参数的变化并不明显,可见,在构建不同岩石材料的颗粒流模型时,参数能起到一定的微调作用,而颗粒的形状尺寸是影响岩石拉压力学特性的主要因素。

通过以上方法建立的Clump结构,可以任意调节Clump结构的大小比例,也可以调节Clump结构内Pabble的半径比例,通过这些结构尺寸比例的调节,来调整颗粒流模型的宏观结构力学特性,进而较为快捷准确地匹配不同的岩石材料。这一方法不仅能够有效克服普通颗粒模型与试验实测值有较大偏差的问题,而且能够提高数值模拟的精度与效率,更真实准确地反映出岩石的宏细观破裂特征。

3.   裂隙花岗岩细观破裂特征

3.1   裂隙花岗岩力学特性试验

试验所选用的花岗岩岩样取自山东新城金矿1020 m水平工作面,岩块质地均匀,总体完整性好。采集的花岗岩材料经试验室精加工,制备成ϕ50 mm×100 mm的圆柱试件及ϕ50 mm×25 mm的圆盘试件,其中裂隙试件在完整试件基础上在试件中心点水力切割不同倾角的贯穿裂隙,试件的加工精度满足国际岩石力学学会建议试验规范要求,岩石单轴压缩及巴西劈裂试验在朝阳仪器有限公司生产的GAW-2000微机控制电液伺服刚性压力试验机上进行。完整试件及裂隙花岗岩试件如图8所示,其中a为12.5 mm,b为10 mm。试验工况及力学参数计算结果如表2所示。

图  8  完整及含裂隙花岗岩试件

Figure  8.  Complete and fractured granite specimens

下载: 全尺寸图片 幻灯片

表  2  试验工况及力学参数

Table  2.  Test conditions and mechanical parameters

单轴压缩试验			巴西圆盘劈裂试验
试件编号	裂隙角度β/(°)	单轴抗压强度σc/MPa	试件编号	裂隙角度γ/(°)	抗拉强度σt/MPa
无裂隙	—	138.73
A	35	55.14	1	—	11.85
	45	60.91	2	0	4.97
B	35	48.03	3	20	3.69
	45	60.68	4	40	4.07
C	35	47.87	5	60	2.27
	45	56.58	6	90	3.11

下载: 导出CSV 

| 显示表格

3.2   颗粒流模型参数标定

采用可变半径比例Clump模型,通过颗粒流PFC^2D软件进行单轴压缩试验及巴西圆盘劈裂试验的数值模拟。数值计算模型尺寸与室内试验试件尺寸一致,Clump模型如图9所示,其中单轴Clump模型的Pebble数为14616,巴西Clump模型的Pebble数为4986。参数标定过程中,基于前文粒径尺寸及力学参数对拉压强度的影响分析,配置花岗岩颗粒流模型,根据花岗岩拉压强度经验值确定Clump R_min为0.6 mm,考虑黏结内聚力c_b与黏结抗拉强度${\sigma }_{\text{b}}$对Clump模型拉压比影响,以及黏结刚度比$k_{\text{n}}^{*}/k_{\text{s}}^{*}$对破裂模式的影响,对表1中的细观力学参数进行微调,确定模型力学参数如表3所示。

图  9  Clump模型示意图

Figure  9.  Diagram of Clump model

下载: 全尺寸图片 幻灯片

表  3  模型力学参数

Table  3.  Mechanical parameters of model

参数	取值	参数	取值
Ball Rmin/mm	0.26	黏结模量E*/GPa	2.0
Ball粒径比Rmin∶Rmax	1.5	黏结刚度比$k_{\text{n}}^{*}/k_{\text{s}}^{*}$	2.5
颗粒刚度比kn/ks	2.0	黏结抗拉强度${\sigma }_{\text{b}}$/MPa	32
颗粒摩擦系数µ	0.2	黏结内聚力$c_{\text{b}}$/MPa	24
Clump Rmin/mm	0.6	黏结摩擦系数µ*	0.2
Clump粒径比Rmin∶Rmax	1.5	摩擦角ƒ/(°)	32

下载: 导出CSV 

| 显示表格

3.3   单轴压缩试验结果及分析

单轴压缩数值模拟与室内试验的轴向应力–应变曲线如图10所示,最终破裂模式如表4所示。

图  10  花岗岩单轴压缩试验与模拟应力–应变曲线对比

Figure  10.  Comparison of uniaxial compression tests on granite and simulated stress-strain curves

下载: 全尺寸图片 幻灯片

表  4  花岗岩单轴压缩试验与模拟破裂特征

Table  4.  Uniaxial compression tests and numerical simulation fracture characteristics of granite

试件类型	无	A		B		C
裂隙角度β/(°)	—	35	45	35	45	35	45
颗粒流计算结果
室内试验结果

下载: 导出CSV 

| 显示表格

由图10（a）可知,无裂隙花岗岩试件单轴压缩数值模拟与室内试验的应力–应变曲线吻合较好,室内试验单轴抗压强度为135.86 MPa,数值模拟单轴抗压强度为141.39 MPa,相对误差为4.07%。同时,在峰值强度之后,轴向应力迅速跌落,试件失稳破坏,表现出典型的脆性特征,由表4可知,数值模拟的破裂模式为沿加载轴线的劈裂破坏,与室内试验结果吻合。

由图10（b）,（c）与（d）可知,35°裂隙试件单轴压缩数值模拟与室内试验的应力–应变曲线吻合较好,其中B-35°试件的数值计算曲线与室内试验曲线的应力峰值吻合,而峰值点对应的应变值相差0.05%；45°裂隙试件数值模拟应力–应变曲线与室内试验应力–应变曲线基本吻合。由表4可知,35°及45°的裂隙试件数值模拟的破裂模式与室内试验结果吻合。

3.4   巴西圆盘劈裂试验结果分析

巴西圆盘劈裂数值模拟与室内试验的拉应力–应变曲线如图11所示,最终破裂模式如表5所示。

图  11  花岗岩巴西劈裂试验与模拟应力–应变曲线对比

Figure  11.  Comparison of Brazilian splitting tests on granite and simulated stress-strain curves

下载: 全尺寸图片 幻灯片

表  5  花岗岩巴西劈裂试验与模拟破裂特征

Table  5.  Brazilian splitting tests and numerical simulation fracture characteristics of granite

裂隙角度γ/(°)	无	0	20	40	60	90
颗粒流计算结果
室内试验结果

下载: 导出CSV 

| 显示表格

由图11（a）可知,无裂隙花岗岩圆盘试件数值模拟应力–应变曲线与室内试验应力–应变曲线吻合较好,室内试验抗拉强度为11.85 MPa,数值模拟抗拉强度为10.37 MPa,相对误差为12.49%。同时,数值模拟结果的破裂模式为沿加载轴线的拉伸破坏,由表5可知,其最终破裂模式与室内试验结果吻合较好。

由图11（a）,（b）与表5可知,20°,40°及60°倾角花岗岩圆盘试件数值计算曲线与室内试验曲线基本吻合,破裂模式吻合较好；0°倾角试件数值计算曲线与室内试验曲线相比略有偏差,破裂模式较为相似；90°倾角试件数值计算曲线与室内试验曲线的应力峰值基本吻合,而峰值点对应的应变值相差0.06%,破裂模式基本吻合。

4.   高应力围岩宏观破裂尺度效应

4.1   数值计算模型

不同尺寸及边界的岩体,其力学性质存在一定差别,形成了岩体力学性质及破裂特征的尺度效应。在构建大尺度颗粒流模型方面,目前仍然面临较大困难,尤其是大尺度Clump模型的构建,按照常规方法会浪费大量的时间。基于前文所提出的可变半径比例Clump模型构造方法,以及所选取的模型细观力学参数,首先建立100 mm×100 mm小尺度的Ball-Brick,然后转化为Clump模型,最后通过Copy功能建立不同尺度的颗粒流模型。值得注意的是,在构建小尺度Ball-Brick模型时,模型的尺度不宜过大,应控制颗粒在1万个内,否则在转化时,会出现边界不能正常接触的情况。

颗粒流模型的边界条件及尺寸如图12所示,模型以边长为L的正方形为边界,在保证模型颗粒密度相同的条件下,L分别取0.5,2.0,8.0 m,圆形洞室的半径取0.3L,侧压力系数$\lambda$=0.5,原岩应力$\sigma$=100 MPa,模型中的颗粒细观结构尺寸按本文校对的花岗岩模型参数选取,3种颗粒模型的Clump个数分别为1.2万,19.1万和304.1万,分析不同尺度的高应力洞室围岩局部破裂特征及尺度效应。

图  12  数值计算模型示意图

Figure  12.  Diagram of numerical model

下载: 全尺寸图片 幻灯片

4.2   洞室围岩局部破裂尺度效应分析

洞室的开挖是一个卸载过程,基于所构建的不同尺度的颗粒流模型,从破裂模式及裂纹扩展这两个方面,研究洞室围岩的宏观破裂尺度效应。围岩模拟开挖后运行10000步的计算结果如图13所示。

图  13  不同尺度围岩破裂区及裂纹扩展区

Figure  13.  Fracture zones and crack propagation zones of surrounding rock under different scales

下载: 全尺寸图片 幻灯片

图13（a）为模型尺寸L = 0.5 m的洞室开挖后围岩破裂区及裂纹扩展区分布特征,洞室开挖后,围岩沿水平方向出现了三角形的破碎区域,该破裂区域的边界清晰,厚度约为0.25L,并且已经逐渐脱离洞室轮廓线向内部移动,围岩破裂区内岩体破碎,但是没有出现明显的剪切带,裂纹扩展区的分布也主要集中在三角区域内。

图13（b）为模型尺寸L = 2.0 m的洞室开挖后围岩破裂区及裂纹扩展区分布特征,开挖后围岩出现了不规则的破碎区域,破裂边界较为清晰,围岩破裂区内完整岩块占破裂区的大部分,并且沿着开挖边界出现了明显的剪切带,裂纹扩展区主要集中在0.13L范围内,相比图13（a）的破裂区域,该模型的破裂范围比例明显减小。

图13（c）为模型尺寸L = 8.0 m的洞室开挖后围岩破裂区及裂纹扩展区分布特征,在保证相同颗粒密度而增大模型尺度条件下,围岩的破裂特征与图13（a）,（b）不同,围岩破裂深度较小,在0.12L范围内,沿着围岩垂直方向的扩展范围较大,破裂形式主要是以剪切滑移破坏和板裂破坏为主,裂纹扩展范围比例明显减小,围岩体现出了明显的高应力局部破裂特征。

5.   结论

（1）通过对比Ball和Clump模型在不同细观力学参数及不同粒径尺寸变化情况下的拉压比可知,Ball和Clump模型的拉压比受细观力学参数的影响较小,Ball模型受粒径尺寸的影响同样较小,而Clump模型受粒径尺寸及比例的影响较大,在构建其它类型岩石的颗粒模型时,可通过调整Clump颗粒结构来体现岩石材料不同的拉压比力学特性。

（2）对比分析裂隙花岗岩试件数值模拟与室内试验的拉、压强度曲线及破裂模式,结果表明二者吻合度较高,本文所提出的Clump建模方法与Ball建模方法相比,可以较好的控制拉压比,得到更符合工程实际的计算结果。

（3）当围岩尺度较小时,围岩的破裂区主要以局部区域破碎为主,而当围岩尺度逐渐增大时,围岩表现出了明显的剪切滑移及板裂破坏特征。可见,基于花岗岩细观颗粒模型的深部围岩破坏特征具有明显的尺度效应。在分析工程问题时,应建立符合工程尺度的颗粒模型,才能在硬岩破裂及支护失效机理等研究方面获得合理的计算结果。