0.   引言

强度与变形属于岩石的重要力学性质,与工程的稳定性有着密切的联系。例如高放射废物地质处置库工程,由于花岗岩具有渗透性低、导热性好、力学强度高与变形小等特点,是处置库工程的理想围岩之一^[1]。

晶质岩如花岗岩、大理岩等由于矿物晶体的颗粒大小、晶体形状、晶体取向以及各矿物晶体之间的相互关系等造成了晶质岩细观结构的非均质性,由此导致了岩石在外部荷载作用下内部应力应变分布的不均匀性。进而造成了室内试验中诸如单轴压缩试验、巴西劈裂试验等测得的强度与变形参数存在一定的离散性^[2]。例如花岗岩中石英、钾长石、黑云母的力学性质（强度、刚度）差异显著,各试验岩样中矿物含量和分布均不同,这将引起强度的离散。因此,岩石内部的细观结构控制着岩石的物理力学性质。

目前,学者主要通过数值模拟手段定量研究岩石的非均质性对强度与变形的影响。如Tang等^[3]通过RFPA软件中岩石单元的力学参数服从Weibull分布来定量表征岩石的非均质性,结果表明岩石力学参数的非均质性对岩石的强度和破坏模式均有较大影响。对岩石单元的力学参数赋予一定的统计分布规律,以此表示岩石的非均质性,但是在试验中获取特定岩石的细观单元的强度分布存在很大困难。

另一种表征非均质性的方法就是细观结构的非均质性。如采用UDEC软件,研究发现细观结构的非均质性对模型的应力分布、起裂应力等有很大影响^[4]。Potyondy在颗粒流软件（particle flow code, PFC）^[5]中引入等效晶质模型（grain-based model, GBM）^[6],以模拟晶质岩的细观结构。通过该等效晶质模型,许多学者对细观结构与岩石物理力学性质之间的关系进行了深入的研究,证明GBM模型对于模拟晶质岩细观结构的可行性和正确性^[7-9]。

岩石强度准则表征岩石在临界损伤状态下的应力与强度之间的关系,被用于研究岩石在复杂应力条件下的屈服条件和破坏规律,已经被广泛应用于边坡、水电站、深埋隧道等各种工程中,一直是岩土工程研究的难点和热点问题。例如线性的Mohr–Coulomb（M–C）强度准则^[10]和非线性的Hoek–Brown (H–B)强度准则^[11],在岩土工程领域应用最多,研究最为广泛。例如Martin等^[12]对Lac du Bonnet（LdB）花岗岩开展循环加卸载试验,以Mohr–Coulomb准则为基础,基于试验结果建立了岩石内聚力和内摩擦角与损伤变量的关系。在Hoek–Brown准则的研究进展中,Zhang等^[13]借助PFC软件,研究了颗粒微观参数如球形颗粒的大小分布和接触强度对Hoek–Brown准则中参数m_i的影响,结果表明参数m_i与微观结构的互锁程度密切相关。但是由于岩石本身细观结构的非均质性,对强度准则中各参数均有一定影响,目前对此研究尚存在不足,有待进一步深入研究。

在花岗岩体内进行深埋隧道或地下硐室的开挖施工,其高应力环境经常导致花岗岩发生岩爆、板裂等形式的脆性破坏,对施工造成严重影响^[14]。研究花岗岩的脆性破坏和发生机制具有重要意义。研究细观结构的非均质性对花岗岩脆性指标的影响能够为工程设计与建设提供一定程度的指导。

本文采用PFC^2D软件中的GBM模型建模方法,以LdB花岗岩^[15]为研究对象,开展了一系列不同围压下的常规三轴压缩试验。通过获取不同细观结构模型的应力–应变曲线,分析了细观结构的非均质性对三轴压缩强度的影响,并将非均质性和Mohr–Coulomb强度准则、Hoek–Brown强度准则建立联系,定量研究非均质性和围压共同作用下对模型峰值强度的影响。最后简要分析了非均质性和花岗岩脆性指标之间的关系,以更深入理解细观结构的非均质性对花岗岩的强度与变形之间的影响规律。

1.   岩石GBM模型

1.1   GBM模型简介

岩石的细观结构具有一定的非均质性,进而影响岩石的物理力学性质,建立的数值模型应当能够在一定程度上反映出岩石的细观结构。Zhou等^[16]采用平直节理接触模型（flat-joint contact model）和光滑节理模型（smooth-joint contact model）相组合,模拟结果证明该方法可以很好地再现花岗岩实际破坏过程。因此本文采用平直节理接触模型和光滑节理接触模型分别模拟晶体内部和晶体边界的方法,研究花岗岩细观结构的非均质性与三轴压缩强度之间的关系。

笔者在文献中[7]已证明所用细观结构能够在一定程度上反映出LdB花岗岩的晶体粒径分布,因此本文继续使用所建立的细观结构,最终模型如图1所示。图1（a）所示模型中,高为125 mm,宽为50 mm,高宽比为2.5,模型的高宽比与Martin室内试验^[15]所用LdB花岗岩岩样高宽比保持一致。

图  1  GBM模型示意图^[7]

Figure  1.  Schematic diagram of GBM model ^[7]

下载: 全尺寸图片 幻灯片

平直节理接触模型将圆形颗粒改造成多边形颗粒,平直节理接触模型的力学行为分为黏结和非黏结两个部分,使得颗粒破坏后的旋转被抑制。因此,相较于平行黏结模型而言,运用平直节理接触模型可以获得更为合理的单轴抗压强度与直接拉伸强度之比,如图2（a）所示。更加详细的解释参见PFC 5.0手册^[17]。

图  2  GBM模型接触模型示意图^[17]

Figure  2.  Schematic diagram of GBM model contact ^[17]

下载: 全尺寸图片 幻灯片

光滑节理模型模拟节理或者晶体边界的力学行为而不受颗粒接触方向的影响,光滑节理经过的相邻颗粒的平直节理模型会被光滑节理接触替代,其方向平行于节理方向,两相邻颗粒可以存在部分重叠或相对滑动,从而避免了沿颗粒表面绕行的行为。如图2（b）所示。更加详细的解释参见PFC 5.0手册^[17]。

1.2   GBM模型参数标定

PFC^2D软件中所使用的细观参数不能通过室内试验直接测量得出,一般通过模拟单轴压缩试验得出的单轴压缩强度、弹性模量和泊松比等与室内试验所得岩石宏观力学参数进行对比,采用不断试错的方法得出。在GBM等效晶质模型中不能直接运用各单晶的力学参数,但是各晶体颗粒间接触模型的细观参数应当大致符合一定的规律：即弹性模量由石英、斜长石、钾长石、黑云母递减,抗拉强度同样符合此规律。为了计算的简便,模型中将所有的晶体边界的细观参数设为一致的,本文最终标定的参数如表1所示。本文所用模型和室内试验结果对比如表2和图3所示。

表  1  GBM模型细观参数表

Table  1.  Meso-scopic parameters of GBM model

种类	颗粒基本参数					晶体内部细观参数（平直节理接触模型）					晶体边界细观参数（光滑节理接触模型）
	最小粒径/mm	半径比	密度/(g·cm-3)	刚度比	模量/GPa	模量/GPa	刚度比	抗拉强度/MPa	黏聚力/MPa	内摩擦角/(°)	法向刚度/(GPa·m-1)	切向刚度/(GPa·m-1)	抗拉强度/MPa	黏聚力/MPa	内摩擦角/(°)
石英	0.3	1.66	2.65	3.5	168	168	3.5	170	170	50	40000	10000	20	20	20
钾长石	0.3	1.66	2.56	1.5	140	140	1.5	145	145	35
斜长石	0.3	1.66	2.63	1.5	140	140	1.5	146	146	35
黑云母	0.3	1.66	3.05	1.0	100	100	1.0	85	85	30

下载: 导出CSV 

| 显示表格

表  2  室内试验^[15]和本文模型宏观力学参数对比表

Table  2.  Comparison of macroscopic mechanical parameters between experimental data^[15] and numerical simulation results

试验类型	单轴抗压强度/MPa	弹性模量/GPa	泊松比	抗拉强度/MPa
室内	200±22	69±5.8	0.26±0.04	7.4±1.04
本文	245.8	66.0	0.24	10.2

下载: 导出CSV 

| 显示表格

图  3  室内试验结果^[15]和本文模型对比

Figure  3.  Comparison between experimental data^[15] and numerical simulation results

下载: 全尺寸图片 幻灯片

2.   非均质性表征

本文采用Liu等^[9]提出的非均质性的测量方法,该方法考虑了晶体颗粒最小粒径、平均粒径和最大粒径,能更好地表征模型的非均质性。

首先计算整个模型的平均晶体粒径d：

式中,d为模型平均晶体粒径,${\overline{d}}_i$,w_i分别为第i种晶体颗粒的平均粒径和含量,共j种晶体颗粒,本文取j=4。

然后计算整个模型的非均质性因子H：

式中,H为非均质性因子；d_i1,d_i2,d_i3分别为第i种晶体颗粒的最小粒径、平均粒径、最大粒径。通过采用MATLAB R2014和犀牛Rhino软件共构建了5个晶体粒径分布不同的模型,再将模型导入PFC^2D 5.0软件中构建GBM等效晶质模型。最终生成的5个模型通过计算后非均质性因子由小到大依次为模型A（H=0.11）、模型B（H=0.29）、模型C（H=0.35）、模型D（H=0.53）、模型E（H=0.65）,见图4^[7]。模型具体的粒径信息详见文献[7]。5个模型中含有4种相同的矿物（钾长石、斜长石、石英、黑云母）,相同的矿物含量。其中,模型C在第1节中被用来标定微观参数。

图  4  不用非均质性因子模型示意图（青色代表黑云母；红色代表钾长石；蓝色代表斜长石；绿色代表石英）^[7]

Figure  4.  Specimens with different heterogeneity indexes^[7]

下载: 全尺寸图片 幻灯片

3.   结果分析

在PFC^2D软件中,通过对左右侧墙进行伺服,使得模型的围压维持稳定,再通过赋予上下墙体一定的速度实现对模型的轴向加载。通过对不同非均质性模型进行三轴压缩模拟,获取应力–应变曲线,进而分析非均质性与强度及相应的强度准则之间的关系。

3.1   应力–应变曲线

下面以5个模型在围压分别为0,5,20,40 MPa时的典型应力–应变曲线为例,分析细观结构的非均质性对应力–应变曲线的影响,如图5所示。

图  5  不同围压下各模型应力–应变曲线

Figure  5.  Stress-strain curves of various models under different confining pressures

下载: 全尺寸图片 幻灯片

由图5可知,不同非均质性因子模型在同一围压下的应力–应变曲线不同,证明细观结构的非均质性对岩样受压加载情况下的应力–应变曲线具有显著的影响。不同围压下各模型的应力–应变曲线大致有压密阶段、弹性变形阶段、峰前屈服阶段和破坏阶段。在PFC常规模拟中,压密段不明显是因为没有考虑自然界中岩样内部存在张开裂隙,因此在受压情况下应力–应变曲线不会有明显下凹的压密段^[18]。在单轴压缩情况下（如图5（a）所示）,各应力–应变曲线显示出明显的脆性特征,轴向应力到达峰值强度以后,峰后应力迅速跌落。随着围压的增加（如图5（c）,（d）所示）,各应力–应变曲线峰前屈服段长度逐渐增加。

3.2   非均质性对变形参数的影响

不同非均质性因子模型单轴压缩条件下的弹性模量E和泊松比$\nu$变化如图6所示。

图  6  不同非均质模型弹性模量和泊松比

Figure  6.  Elastic moduli and Poisson's ratios of different heterogeneous models

下载: 全尺寸图片 幻灯片

弹性模量和泊松比表示岩石材料的变形与所受荷载之间的关系,反映其弹性性质。由图6可知,随着非均质因子增加,模型由均质变为非均质,弹性模量逐渐减小,泊松比增加,说明细观结构的非均质性对模型的变形影响很大。越均质的模型（如模型A）,其晶体颗粒粒径分布越均匀,岩石达到峰值强度前,在相同的轴向应力增量下,产生的轴向应变比非均质的模型小,即弹性模量大,抵抗变形的能力强。在同一应力水平下,均质性程度高的模型横向变形更小,导致泊松比更小。随着非均质性因子的增加（如模型E）,模型抵抗变形的能力减弱,弹性减小,各晶体颗粒大小不一,且存在粒径偏大的晶体（如模型E中的钾长石）,导致模型在单轴加载情况下横向变形增加,最终导致泊松比增加。

3.3   非均质性对强度的影响

5个不同均质性因子的模型在不同围压的峰值强度如图7所示,直接拉伸强度如图8所示。

图  7  不同非均质模型压缩强度变化

Figure  7.  Variation of compressive strength of different heterogeneous models

下载: 全尺寸图片 幻灯片

图  8  直接拉伸强度与非均质性因子关系

Figure  8.  Relationship between direct tensile strength and heterogeneity factor

下载: 全尺寸图片 幻灯片

由图7可知,各围压下,随着非均质性因子的增加,模型抗压强度减小,如在围压30 MPa下,非均质性因子H=0.11时,抗压强度为672.23 MPa,H=0.65时,抗压强度为465.97 MPa,较前者减小206.26 MPa,减小了30.7%。由图8可知,随着非均质性因子的增加,模型更加非均质,其抗拉强度更小。5个模型的直接拉伸强度在9.8～10.5 MPa。在室内试验中同样发现细粒花岗岩（较均质）强度比粗粒花岗岩（较不均质）强度高^[19]。由此看见,花岗岩的细观结构对抗压强度具有显著影响。

花岗岩中不同矿物间的力学差异所造成的力学参数的非均质性和由于矿物晶体粒径不同所造成的结构非均质性都会引起岩样内部强度的非均匀性,岩样在外部荷载的作用下会产生不均衡受力。例如当岩石受压应力作用时,矿物晶体粒径尺寸和形状的几何变化会引起较大的拉应力场。因此,由结构非均质性引起的局部拉应力对拉裂纹的萌生、扩展和贯通起主导作用。而矿物结构非均质性越小,岩石内部应力分布越均匀,导致了更高的抗压强度。花岗岩中黑云母力学参数较弱,且粒径较其他矿物更小,在外部荷载下会首先发生变形,而石英的强度较大,在同一应力水平下变形较黑云母较小,黑云母变形后会将应力传递给周围矿物,当遇到力学性质存在差异的矿物时（如石英或者长石）,就会产生拉应力或者剪应力。当非均质性越大时,产生的拉应力或者剪应力就会更大,进而为微裂纹的萌生提供了更加容易的地点；非均质性越大也导致晶体边界更长,为微裂纹的扩展提供了更加连续的路径。

Mohr–Coulomb强度准则和Hoek–Brown强度准则被广泛应用于岩土工程的设计与分析中。Mohr–Coulomb强度准则物理意义明确,认为岩石破坏的最大剪切应力由内聚力c和正应力与内摩擦角$\phi$的正切值共同确定。其最大主应力（抗压强度）与围压的表达式为

式中,${\sigma }_1$为最大主应力,${\sigma }_3$为最小主应力,Q为单轴抗压强度,K为围压对强度的影响系数。内摩擦角$\phi$和内聚力c的关系为

Hoek-Brown强度准则是Hoek等^[11]通过对野外大量现场岩体试验成果的统计分析而提出的,可以反映岩石和岩体固有的非线性破坏特点,其表达式为

式中 ${\sigma }_1$,${\sigma }_3$,${\sigma }_{\text{c}}$分别为最大主应力、最小主应力、岩石单轴抗压强度；m_i为量纲为“1”的经验参数,反映岩石的软硬程度。

由于Mohr–Coulomb强度准则中考虑抗拉强度值进行拟合的话,会造成较大误差,所以本文在进行Mohr–Coulomb强度准则拟合时,不考虑模型的抗拉强度。采用Mohr–Coulomb强度准则和Hoek–Brown强度准则对不同围压下各模型的抗压强度进行拟合,所得结果分别如图9,10所示。拟合参数如表3所示。

图  9  Mohr-Coulomb强度准则拟合结果

Figure  9.  Fitting results of Mohr-Coulomb strength criterion

下载: 全尺寸图片 幻灯片

图  10  Hoek-Brown强度准则拟合结果

Figure  10.  Fitting results of Hoek-Brown strength criterion

下载: 全尺寸图片 幻灯片

表  3  强度准则拟合参数

Table  3.  Fitting parameters of strength criterion

模型	非均质因子H	M-C准则			H-B准则
		Q/MPa	K	R2	${\sigma }_{\text{c}}$/MPa	$m_i$	R2
A	0.11	269.9	11.6	0.98	233.0	45.6	0.98
B	0.29	263.7	10.6	0.98	230.2	39.5	0.98
C	0.35	242.3	8.2	0.98	219.1	25.8	0.98
D	0.53	234.2	7.9	0.99	216.5	24.3	0.98
E	0.65	234.0	7.3	0.99	216.6	21.9	0.99

下载: 导出CSV 

| 显示表格

由图9,10和表3可知,采用Mohr-Coulomb强度准则和Hoek–Brown强度准则对各非均质性因子的模型抗压强度拟合效果都较好。由于本文所模拟的花岗岩在较低围压下的强度变化,拟合效果均较为理想,已有的研究表明,当围压增大时,Mohr–Coulomb强度准则和Hoek–Brown强度准则会存在较大偏差,对抗压强度的预测值偏高^[20]。

由图9（b）可知,随着非均质因子的增加,模型由均质转变为非均质,围压影响系数K逐渐减小,证明围压对抗压强度的影响逐渐降低。当岩石越均质时,围压的增加对强度的提高效果显著,而随着岩石变得非均质时,效果减弱。由图10（b）可知,非均质性因子对Hoek–Brown强度准则中的材料参数m_i影响显著。当非均质因子增加时,材料参数m_i减小,证明模型变得更加易于破坏,强度更低。

依据式（4）,（5）可以确定各模型的内聚力和内摩擦角,其结果如图11所示。

图  11  不同非均质模型内聚力和内摩擦角

Figure  11.  Cohesions and internal friction angles of different heterogeneous models

下载: 全尺寸图片 幻灯片

由图11可知,随着非均质因子的增加,模型的内聚力增加,内摩擦角减小。Mohr–Coulomb强度准则认为岩石一点处的强度由内聚力和摩擦力所组成,但是由于岩石本身具有非均质性,各处不可能同时达到承载极限而破坏^[2]。内聚力主要由晶体颗粒之间的咬合作用形成,当非均质性增加时,模型的粒径分布更加不均匀,大小颗粒之间的咬合作用增加,导致内聚力增加。内摩擦角受晶体颗粒之间的相互嵌入、联锁形成的咬合力和各晶体边界的摩擦作用力影响。本文所建立的5个非均质模型中,采用光滑节理接触模型模拟花岗岩晶体边界,当模型非均质性增加时,容易在晶体边界处产生应力集中现象,使得边界处的接触变得更加容易断裂,从而导致模型的内摩擦角减小。

将M–C强度准则和H–B强度准则中的参数和非均质性因子进行拟合（见图9（b）,10（b））,然后代入到强度准则中,可以构建围压和非均质性因子对抗压强度的影响。其最终表达式如下：

式（7）为M–C强度准则表达式,式（8）为H–B强度准则表达式。

两者的图像如图12所示,当围压越高、非均质因子越小、模型越均质时,其抗压强度越高。表明细观结构的非均质性对岩石在受压情况下的抗压强度影响显著。

图  12  M–C强度准则和H–B强度准则与非均质因子关系示意图

Figure  12.  Diagram　of relationship among M–C intensity criterion, H–B strength criterion and heterogeneity factor

下载: 全尺寸图片 幻灯片

3.4   非均质性对岩石脆性的影响

为了更好地理解非均质性对花岗岩的脆性破坏的影响,下面简单采取3个脆性指标对不同非均质因子的模型进行脆性评价。脆性指标B_I1^[21],B_I2^[22]和B_I3^[22]计算表达式为

式中 E,$\nu$分别为单轴压缩情况下的弹性模量和泊松比；$\phi$为内摩擦角；${\sigma }_{\text{c}},{\sigma }_{\text{t}}$分别为单轴压缩强度和直接拉伸强度。3个脆性指标越大时,表明模型脆性越强。脆性指标随非均质因子H的变化如图13所示。

图  13  脆性指标随非均质因子变化关系

Figure  13.  Relationship between brittleness index and heterogeneity factor

下载: 全尺寸图片 幻灯片

脆性破坏是由于细观结构上的非均质性,导致岩石在加载情况下出现非均匀应力,产生局部破坏,并逐渐形成多维破裂面的过程。岩石的矿物组成、晶体粒径的尺寸分布等对岩石的脆性或者塑性行为具有显著影响。细粒花岗岩比粗粒花岗岩具有更高的脆性,同时等粒结构的岩石比不等粒结构的岩石强度更高,脆性更强^[23]。花岗岩颗粒粒径较小时,一般来说花岗岩更加均质。由图13可知,随着非均质因子的增加,脆性指标B_I1,B_I2,B_I3均呈线性减小,且和非均质因子的线性拟合效果较好（R²分别为0.84,0.95,0.92）。本文所建立的5个模型中,模型A属于非均质因子小、等粒结构；模型E属于非均质因子大、似斑状不等粒结构。因此模型A的脆性指标比模型E的脆性指标均更大,表明模型A的脆性更强。本文模拟结果和室内试验^[19,23]的结果吻合较好。当模型更加均质时,微裂纹通常沿着垂直的晶体边界萌生扩展,最后形成垂直的宏观裂纹,模型以劈裂破坏为主；当模型更加非均质时,因含有较大粒径的长石类矿物,通常具有较长的晶体边界,使得微裂纹更加容易沿着晶体边界扩展,形成剪切破坏^[8]。因此在评价岩石的脆性时,把细观结构所造成的非均质性考虑进去可能会更加合理。

4.   结论

本文通过GBM模型建模方法,开展了一系列不同围压下的花岗岩常规三轴压缩试验研究。揭示了晶体粒径分布所造成的细观结构上的非均质性对强度和变形参数、Mohr–Coulomb强度准则和Hoek–Brown强度准则中各参数的关系以及对岩石脆性的影响。主要得到以下3点结论。

（1）细观结构的非均质性对岩样受压加载情况下的应力–应变曲线具有显著的影响。随着非均质因子增加,模型的弹性模量逐渐减小,泊松比增加。

（2）各围压下,随着非均质性因子的增加,模型抗压强度减小。当非均质性越大时,为微裂纹的萌生提供了更加容易的地点,同时为微裂纹的扩展提供了更加连续的路径。同时模型的内聚力增加、内摩擦角减小。对同一类花岗岩而言,Mohr–Coulomb强度准则和Hoek-Brown强度准则中各参数的取值并非恒定,而是一定程度上依赖于内部的细观结构造成的非均质性。

（3）岩石细观结构上的矿物粒径及粒径分布的非均质特征对岩石脆性影响明显。细观结构的非均质性对岩石的脆性指标影响较大,随着非均质性因子的增加,岩石脆性指标减小,脆性降低。