岩石材料在细观上由矿物、胶结材料和微孔隙等构成，且构成方式及矿物形状、尺寸、组分等在空间上具有一定的随机性，由此造成了岩石材料的非均匀性或非均质性，并导致岩石力学参数和变形参数在空间上也具有一定的随机性，如弹性模量、粘聚力和抗拉强度等，直接影响到岩石试样的力学特性和破坏行为，也影响隧洞开挖后围岩的破裂碎胀变形特征^[1]。

准确描述实验室岩样或隧洞围岩每个矿物颗粒及胶结材料的力学行为和变形特征是不必要的，同时也是难以做到的，因此，可采用数值模拟方法，一般分为连续性方法(如有限元法FEM^[2-3]、有限差分法FLAC^[4])、非连续方法(如颗粒元法PFC^[5]、离散单元法UDEC^[6-7])和耦合方法(如有限元-离散元耦合方法FDEM^[8-9])；而对非均质性的表征则可分为矿物晶体模型法、数字图像处理法和随机参数赋值法。

对于矿物晶体模型法(GBM)，Potyondy等首次构建了PFC单轴压缩数值模型^[10]和巴西劈裂模型^[11]，即PFC-GBM模型，ZHOU等^[5]和LI等^[12]对该模型进行了改进；PENG等^[13]和LIU等^[14]采用PFC-GBM模型研究了矿物颗粒圆滑度对岩样单轴压缩力学行为和破裂形态的影响；谭鑫等^[6]基于Voronoi块体随机离散方法建立了含层理构造的非均质片麻岩巴西劈裂UDEC数值模型；LAN等^[7]基于GBM建立了非均质岩样单轴压缩UDEC数值模型；ABDELAZIZ等^[15]基于GBM建立了FDEM非均质岩样单轴压缩模型。

对于数字图像处理法(DIP)，陈沙等^[16]详细介绍了DIP技术原理及实现方法，并将数字图像与FLAC相结合，研究了非均质花岗岩巴西劈裂裂纹扩展过程；徐文杰等^[17]和丁秀丽等^[18]基于DIP建立了土石混合体真实细观结构PFC^2D数值计算模型；刘占新等^[19]采用扫描电镜图像处理建立了包含孔裂隙、矿物质和煤基质的非均质单轴压缩煤样PFC数值模型；朱泽奇等^[20]基于DIP建立了FLAC非均质岩样单轴压缩数值模型；张世瑞等^[8]结合GBM和DIP，构建了花岗岩细观FDEM数值模型；TAN等^[21]采用DIP建立了非均质花岗岩巴西劈裂UDEC数值模型。

对于随机参数赋值法，唐春安等^[22]提出了岩石强度参数服从Weibull分布的RPFA数值分析方法；LIU等^[23]基于RFPA数值方法和有限元法，提出了一种新的数值方法(R-T^2D)以模拟岩石的非均质特征；王旭一^[24]基于Weibull分布非均质线性平行粘结模型(WLPB)和Weibull分布非均质光滑节理模型(WSJ)构建了PFC非均质层状岩体数值模型；夏海城等^[1]基于弹性模量服从正态分布的假设建立了非均质岩样DDA数值模型；罗荣等^[25-26]提出了岩石矿物细胞元随机性参数赋值方法；冯增朝等^[27]对岩石细胞元进行随机参数赋值时，假定细胞元的强度、弹性模量等的变化是成比例的；MANOUCHEHRIAN等^[3]采用有限元软件ABAQUS研究了考虑非均质特征下的岩石单轴压缩力学特征及破坏行为，其非均质性通过随机参数赋值法来体现，包括弹性模量、粘聚力和内摩擦角；程昊等^[4]提出了考虑岩石弹性模量非均匀分布的FLAC^3D数值模型，文志杰等^[28]对上述模型做了进一步修正，即考虑了岩石弹性模量的峰后损伤；李昂等^[29]在非均质单轴压缩FLAC^3D数值模型中考虑了弹性模量和粘聚力服从Weibull分布，同时也考虑了抗拉强度和粘聚力随塑性应变参量的变化；王学滨等^[30]基于Weibull分布建立了FLAC隧洞非均质围岩数值模型；杜修力等^[31]采用随机参数赋值法建立了考虑非均质特征的混凝土单轴压缩和直接拉伸有限元数值模型。

总结现有研究发现，还存在如下不足之处：

1) 对于矿物晶体模型法，一般采用Voronoi块体表征多边形矿物颗粒，并在Voronoi块体内部划分更细小的块体/颗粒，如PFC圆球颗粒^{[5, 10, 12-14]}、UDEC多边形块体^{[7, 32]}和FDEM三角形单元块体^{[8, 15, 33]}等；其优势在于能够模拟多边形矿物颗粒及岩样的穿晶断裂和沿晶断裂；然而，Voronoi多边形块体形状相对规则、尺寸相对均匀、边界平直，与真实矿物空间分布和形状、尺寸、矿物边界不符；在数值模拟中，需要将Voronoi块体划分足够小，随即在Voronoi块体内部再划分更为细小的颗粒/块体会造成颗粒/块体数目激增，降低计算效率，因而难以用于工程尺度的模拟，仅在实验室尺度得到较多的研究，因为对于工程尺度计算而言，将Voronoi块体划分至矿物尺度、并在Voronoi块体内部再划分计算网格是不现实的；当然，对于GBM而言，Voronoi多边形块体并非唯一的表征方法，亦可采用其他形状颗粒，如圆形^[34]和其他不规则形状的颗粒/块体^[8]等，但亦存在上述研究尺度局限性的问题。

2) 对于数字图像处理法，其原理是通过拍照或扫描获取岩样断面的矿物空间分布形态和颗粒/块体几何形态，随后对图像进行灰度处理获得数字图像，沿着不同矿物的边界进行建模，真实刻画岩样的实际断面形态、不同矿物的空间分布形态及矿物颗粒自身的几何形态^[35]，已应用于诸多数值模拟方法，如有限元法^[36-37]、UDEC^[21]、PFC^[17-18]、FLAC^{[20, 38]}和FDEM^[39]等。其优势在于能够真实重构岩样剖面的矿物组成和几何形态。然而也面临诸多问题，如仅能对具体断面进行二维建模，难以建立三维模型；若仅对某一个二维剖面进行建模分析，难以反映真实岩样的破裂特征。为克服上述不足，更为先进的X射线CT扫描和逐层切片成像技术已被应用起来^[40-41]，但存在成本高昂、仅适用于小尺度岩样等不足；此外，面对诸如隧洞工程尺度，难以通过拍照扫描获得深入围岩内部的岩体非均质信息，因此也大都仅在实验室尺度得到了研究。

3) 对于随机参数赋值法，首先根据均质材料的方法建立数值模型，随后对该模型的参数进行随机赋值，典型如唐春安等提出RFPA法^[22]。其他数值方法也可采用随机参数赋值法，如DDA ^[1]、UDEC^[42]、PFC^[43]、有限元法^{[3, 34]}、FLAC^{[4, 28, 44]}、岩石矿物细胞元法^[25-27]、虚内键模型法^[45]、元胞自动机法^[46]和光滑粒子流体动力学法^[47]等。岩石材料的非均质性一般体现在弹性模量、泊松比、密度、抗拉强度、粘聚力和内摩擦角上，其概率密度函数一般服从Weibull分布，并根据概率密度函数参数的不同，体现不同的非均质度。随机参数赋值法的优势在于不局限于对特定的岩样进行分析，能够得到更为普适性的结论和规律，且因为采用均质岩体的形式进行建模，仅在赋参时候将力学参数和变形参数根据Weibull分布的方式进行赋值，因此也适用于诸如隧洞工程尺度的模拟^{[24, 30, 48-50]}。此外，该方法在建立三维数值模型方面亦是便捷高效的^[48]，因为建模方法与均质岩样无异，其非均质性仅体现在赋参方面。然而当前大都采用单一的连续或非连续数值方法进行随机参数赋值，而耦合方法的研究尚未见诸报道。

因此，本文将采用基于随机参数赋值法的FDEM数值模拟手段开展非均质隧洞围岩破裂碎胀大变形机理研究，且围岩力学参数概率密度服从Weibull分布。相比于现有非均质隧洞围岩变形破裂数值模拟方法(如FLAC3D^[30]、有限元法^[50]、RFPA^[48-49]和PFC^[24]等)，FDEM能够模拟隧洞开挖过程中完整围岩从弹塑性连续变形至断裂失效全过程，即围岩裂隙网络孕育萌生、起裂、交叉贯通和扩展过程，再现具有真实开度的宏观裂隙，且能够模拟破断后块体间的相互接触挤压效应^[51]。实际上，FDEM已被广泛应用于均质岩体^[51-52]和各向异性岩体^[53-54]隧洞开挖的模拟，在非均质岩体变形和破坏机制方面，已有学者基于FDEM采用GBM^{[8, 15, 33]}和DIP^[39]进行了研究，但仍存在前文所述不足，即仅适用于实验室小尺度岩样的模拟，尚无采用随机参数赋值法对隧洞非均质围岩破裂碎胀大变形机制的研究。

1   基于Weibull分布的非均质岩体FDEM数值模拟原理

1.1   FDEM原理概述

在FDEM中，将材料划分为三角形单元和四边形节理单元，如图1(a)所示。有限元特征通过三角形单元来体现：对每个三角形单元建立独立的应力求解矩阵，根据当前时步与上一时步的坐标差求解应变增量并消除刚体运动，根据应变增量求解三角形变形柯西应力；离散元特征通过四边形节理单元的破裂和三角形单元接触来体现：根据四边形节理单元节点间的相对位移值，判定其处于何种状态(弹性变形或塑性屈服)，当达到极限位移时，节理单元破裂并退出粘结应力计算流程，产生一条微裂隙，其两侧的三角形单元由粘结转换为块体接触，并进行接触检索和接触应力计算。对于宏观裂隙而言，断裂过程区(FPZ)的存在能够真实模拟裂隙尖端的应力状态，避免了应力的异常集中，如图1(c)所示。三角形单元变形柯西应力计算、四边形节理单元粘结应力计算、三角形单元接触检索、接触应力计算和节点坐标更新可在诸多文献[55-56]中找到，本节主要介绍非均质特性在FDEM中的表征方法。须指出，当四边形节理单元发生压剪破坏时，为了实现三角形单元由粘结关系平稳过渡至接触关系、消除应力异常震荡现象，可采用笔者提出的新型法向接触刚度算法^[55]，且四边形节理单元本构模型采用如下改进形式以防节理单元发生压剪破坏时残余剪应力降至为0^[55](图1(b))：

式中：σ_n、τ分别为法向和切向应力(σ_n>0代表拉伸应力，σ_n<0代表压缩应力)；o、s分别为拉伸位移和剪切位移；o_p、s_p分别为法向和切向峰值位移；o_t、s_t分别为法向和切向极限位移；c为粘聚力；φ_i为内摩擦角；f_t为抗拉强度；z为峰后软化函数^[56]。

图  1  FDEM数值模拟基本原理

Figure  1.  Basic principles of FDEM numerical simulation

下载: 全尺寸图片 幻灯片

1.2   非均质性在FDEM中的表征原理

唐春安等^[22]指出，岩石材料在细观上其变形参数(主要为弹性模量)和强度参数(主要为粘聚力和抗拉强度)服从Weibull分布，并据此开发了岩石破裂过程分析软件RFPA。本文亦据此假设开发考虑岩石非均质特征的FDEM数值模拟软件。尽管随机参数赋值方法在诸多文献中有描述^[26]，但在FDEM中尚未得到研究，而FDEM有其特殊性，如弹性模量赋值在三角形单元上、而粘聚力和抗拉强度赋值在四边形节理单元上，因此，有必要对随机参数在FDEM中的赋值方法进行说明。

在任一数值模型中，岩石弹性模量服从Weibull分布，其概率密度函数如下^[22]：

式中：E为任意三角形单元的弹性模量；E₀为弹性模量期望值；m为非均质度，m值越大，岩样越均质，当m趋于无穷大时，即为均质岩样。式(3)所对应的累积分布函数为^[26]：

不同m值下的岩石弹性模量(本文弹性模量期望值E₀=25 GPa)概率密度函数和累积分布函数如图2所示。

图  2  不同m值弹性模量概率密度和累积分布

Figure  2.  Probability density and cumulative distribution of elastic modulus with different m values

下载: 全尺寸图片 幻灯片

根据式(4)可知，任意三角形单元的弹性模量E为：

F(E)为区间(0, 1)的均匀分布函数。因此，在FDEM中，三角形单元弹性模量的赋值流程概括为：采用C语言rand函数产生并获得0~1间的随机数，即为F(x)函数值；随后得到任意一三角形的弹性模量值E，并对所有三角形单元循环得到相应的弹性模量值；因为式(4)、式(5)是根据式(3)而来，因此采用式(5)得到的三角形单元弹性模量分布必然服从式(3)所示的形式。不同m值得到的单轴压缩岩样弹性模量分布如图3所示，其随机分布特征与真实的花岗岩随机性相近^[57]。

弹性模量随机赋值方法一方面能够确保得到的具有不同弹性模量的三角形单元数目服从Weibull概率密度函数形式，如图4所示；另一方面，在空间上，采用Gmsh软件划分的非结构化三角形网格在编号上具有随机性，即相邻三角形单元在编号上不连续，如图5所示，并且C语言rand函数能够产生任意随机数，因此，当对所有的三角形单元进行C语言for循环时，能够确保具有不同弹性模量的三角形单元在空间分布上具有随机性，如图3所示。

图  3  不同m值下单轴压缩岩样弹性模量分布云图

Figure  3.  Nephogram of elastic modulus distribution of uniaxial compression rock samples with different m values

下载: 全尺寸图片 幻灯片

图  4  不同m值下具有不同弹性模量值的三角形单元数目分布(以m=1.5和m=4为例)

Figure  4.  Number distribution of triangular elements with different elastic moduli with different m values (taking m=1.5 and m=4 as examples)

下载: 全尺寸图片 幻灯片

冯增朝等^[27]指出，矿物颗粒的粘聚力及抗拉强度与弹性模量成正比关系，即弹性模量越大的颗粒，其粘聚力和抗拉强度亦越大，并服从如下关系：

式中：E_i、c_i和f_{t, i}为颗粒/块体i的弹性模量、粘聚力和抗拉强度；E_j、c_j、f_{t, j}则为颗粒/块体j相应的参数。对于FDEM而言，弹性模量体现在三角形单元上，而粘聚力和抗拉强度则体现在四边形节理单元上。因此，如图5所示，若将1706号三角形单元记为i单元、5286号三角形单元记为j单元，它们所夹的四边形节理单元记为ij单元，则ij单元的粘聚力和抗拉强度、以及所对应的I型和II型断裂能分别为：

式中：c_ij、f_t,ij、G_I,ij和G_II,ij分别为四边形节理单元ij的粘聚力、抗拉强度、I型断裂能和II型断裂能；c₀、f_t,0、G_I,0和G_II,0分别为粘聚力、抗拉强度、I型断裂能和II型断裂能的期望值。对其他所有四边形节理单元粘聚力、抗拉强度、I型断裂能和II型断裂能亦采用该种方法进行赋值。

图  5  三角形单元编号(试样尺寸和网格尺寸与图4相同)

Figure  5.  Triangle element number (sample size and grid size are the same as those in Fig. 4)

下载: 全尺寸图片 幻灯片

2   参数标定与验证

2.1   参数标定

对三角形单元赋值的前提是需确定弹性模量、粘聚力、抗拉强度、I型和II型断裂能等的期望值，即在均质、各向同性条件下的对应参数值。采用单轴和巴西劈裂模拟标定均质材料的输入参数，并采用三轴压缩模拟验证被标参数的可靠性。数值模型如图6所示，单轴压缩或三轴压缩模型为标准试样尺寸，即为50 mm×100 mm，网格尺寸h=1.3 mm，该网格尺寸符合笔者所标定的网格尺寸范围^[58]，通过上、下各为0.05 m/s的加载速率获得实际为0.1 m/s的加载速率，该加载速率符合笔者^[58]和TATONE等^[59]提出的范围，加载板与岩样间的摩擦系数f_r=0.1，亦符合笔者所提出的取值范围^[58]。巴西劈裂为直径50 mm的圆盘，如图6(b)所示，其他条件与单轴压缩相同。

图  6  单轴/三轴压缩(σ₃=0时为单轴压缩)和巴西劈裂数值模型　/mm

Figure  6.  Numerical models of uniaxial/triaxial compression (when σ₃=0, it is a uniaxial compression model) and Brazilian splitting

下载: 全尺寸图片 幻灯片

采用表1输入参数得到如图7所示的模拟结果(图7中红色为拉伸裂隙、黄色为剪切裂隙)，对于单轴压缩而言，试件产生2条明显的X型共轭剪切带，且与水平面夹角约为62°，这与理论值非常接近(理论破断角为$\dfrac{{{\varphi _i}}}{2} + \dfrac{\pi }{4}$=62.5°^[60])；此外，单轴抗压强度模拟结果为51.49 MPa，亦与理论值非常接近，因为当采用摩尔-库仑强度准则时，单轴抗压强度σ_c理论值为^[61]：

表  1  FDEM模拟输入参数

Table  1.  Input parameters of FDEM simulation

参数	值		参数	值
计算时步Δt/s	0.8×10−9		I型断裂能GI/(J·m−2)	7000
密度ρ/(kg·m−3)	2700		II型断裂能GII/(J·m−2)	27 000
弹性模量E/GPa	25		内摩擦角φi/(°)	35
泊松比ν	0.25		节理单元罚值Pf/GPa	625
室内试验模拟最小 网格尺寸h/mm	1.3		滑动摩擦角φr/(°)	35
粘滞阻尼[56]μ/(kg·(m·s−1))	$2h\sqrt {E\rho }$		法向接触刚度Pn/GPa	90.5
抗拉强度ft/MPa	5.0		切向接触刚度Pt/GPa	62.5
粘聚力c/MPa	13

下载: 导出CSV 

| 显示表格

图  7  单轴压缩和巴西劈裂标定试验模拟结果

Figure  7.  Simulation results of uniaxial compression and Brazilian splitting

下载: 全尺寸图片 幻灯片

巴西劈裂是一种间接测量抗拉强度的方法，抗拉强度σ_t为^[58]：

式中：P为加载板荷载；D为试件直径。巴西劈裂模拟结果如图7(b)所示，首先，试样产生一条中央拉伸裂隙，且从试件中心起裂，符合巴西劈裂间接测量抗拉强度的起裂要求；其次，抗拉强度模拟结果5.09 MPa，与输入值非常接近。须指出，单轴抗压强度和巴西劈裂抗拉强度模拟值均略微偏大，主要是由于模拟加载速率远大于室内试验速率造成的^[62]，此外，网格拓扑结构也会对模拟强度产生轻微影响。

2.2   参数验证

图7模拟结果表明表1的输入参数是可靠的，下面将采用三轴压缩模拟试验再次验证。数值模型如图6(a)所示，围压σ₃分别设为−2 MPa (即围压为拉伸应力)、0 MPa (即单轴压缩)、2 MPa、4 MPa、6 MPa、8 MPa和10 MPa，得到如图8所示的模拟结果。三轴压缩理论强度σ₁为^[61]：

由图8(a)可知，在不同围压下，三轴抗压强度模拟值和理论值均非常吻合。此外，根据围压σ₃和强度σ₁的关系可绘制出如图8(b)所示的摩尔圆，根据摩尔圆得到的粘聚力c和内摩擦角φ_i均与输入值相吻合，如图8(b)所示，因此，表1的输入参数是可靠的。

图  8  三轴压缩模拟结果

Figure  8.  Simulation results of triaxial compression

下载: 全尺寸图片 幻灯片

3   非均质岩样三轴压缩和单轴压缩模拟

3.1   三轴压缩

采用图7(a)所示三轴压缩数值模型进行非均质岩样FDEM数值模拟，计算时步为0.5×10⁻⁹ s，采用1.2节所述原理对岩样弹性模量、粘聚力、抗拉强度及相应的I型、II型断裂能进行赋值，m值分别设为1.5、2、3、4、6、8、10和12，不同m值下的试样弹性模量云图如图4所示，围压分别设定为0 MPa、2 MPa、4 MPa、6 MPa、8 MPa和10 MPa，得到如图9~图10所示模拟结果，可知：

1) 与均质岩样相似(图8(a))，对于具有不同均质度(即m值不同)的非均质岩样，随着围压σ₃的增大，三轴抗压强度σ₁亦呈线性增加，如图9(a)所示。

2) 然而，在相同围压下，单轴抗压强度(当σ₃=0时)或三轴抗压强度均随着m的增大呈指数增长，如图9(b)所示。

3) 与图8(b)类似，绘制不同m值下的摩尔圆，根据摩尔圆得到不同m值对应的等效粘聚力${c^\prime }$和等效内摩擦角${\varphi _{\rm i}^\prime }$，如图9(c)所示，可知，与单轴抗压强度或三轴抗压强度相似，等效粘聚力${c^\prime }$和等效内摩擦角${\varphi _{\rm i}^\prime}$也随着m值的增大呈指数增长，且不断向均质岩样的粘聚力和内摩擦角逼近，即向各自的期望值逼近。

4) 随着非均质度m的增大，不同岩样呈现出相似的破坏特征，即均大致沿着各自的理论破坏角$\dfrac{{{\varphi _{\rm i}^\prime} }}{2} + \dfrac{\pi }{4}$发生剪切破坏，并伴随少量拉伸裂隙，如图10所示。

图  9  不同m值下三轴压缩模拟结果

Figure  9.  Simulation results of triaxial compression with different m values

下载: 全尺寸图片 幻灯片

图  10  不同m值下的三轴压缩岩样破坏模式(以围压σ₃=2 MPa为例)

Figure  10.  Failure modes of rock samples under triaxial compression with different m values (taking σ₃=2 MPa as examples)

下载: 全尺寸图片 幻灯片

3.2   单轴压缩

尽管3.1节已对不同非均质度岩样的三轴压缩力学特征和破坏行为进行了研究，相比于三轴压缩，单轴压缩状态下的非均质岩样力学特征和破坏行为有其独特性，本节将详细分析。不同非均质度下岩样单轴压缩模拟结果如图11~图12所示，可知：

1) 随着m值的增大，除了单轴抗压强度σ_c外，弹性模量E也呈指数函数增大，且不断向均质岩样逼近，如图11(a)所示，其中弹性模量被定义为单轴抗压强度1/2时的割线模量。

2) 对于应力-应变全过程曲线而言，不同m值所对应的非均质岩样除了峰前斜率(即弹性模量)和峰值抗压强度有所区别外，峰后斜率亦有很大不同，如图11(b)所示。以m=1.5和m=12为例，当m=1.5时，应力-应变曲线达到峰值后缓慢下降，表现出明显的峰后塑性和残余抗压能力；而当m=12时，峰后曲线则迅速跌落，表现出脆性特征，尤其均质岩样更为明显；这与唐春安等^[22]采用RFPA得到的模拟结果极为相似。

图  11  不同m值下非均质岩样单轴压缩模拟结果

Figure  11.  Simulation results of uniaxial compression of heterogeneous rock samples with different m values

下载: 全尺寸图片 幻灯片

不同m值呈现出不同的峰后特性可用图12所示的非均质岩样单轴压缩破裂过程解释，对于单轴压缩状态下的试样而言，不同弹性模量的单元应变值是相同的，故而E值大的单元其变形应力也大，反之亦然。然而，尽管E值小的单元其变形应力小，但是由于该单元的抗剪强度和抗拉强度亦较低(据式(7))，因此，裂纹起裂和扩展过程将沿着弹模较小的三角形单元边界，也即沿着强度较小的四边形节理单元进行；此外，由于FDEM采用了摩尔-库仑强度准则，因此，主剪切裂隙带与水平面约呈$\dfrac{{{\varphi _{\rm i}^\prime} }}{2} + \dfrac{\pi }{4}$角度，如图12(a)所示，对于m=1.5的非均质岩样，主剪切裂隙带沿着与水平面呈$\dfrac{{{\varphi _{\rm i}^\prime} }}{2} + \dfrac{\pi }{4}$角度扩展，且在扩展过程中沿着强度较小的四边形节理单元进行，并由此造成主剪切裂隙带宽度较大、剪切面不平整、不光滑，这是由于每个单元的弹性模量、应力状态和抗拉/抗剪强度不尽相同，且分布极不均匀造成的；对于m=12的非均质岩样而言，主剪切裂隙仍与水平面约呈$\dfrac{{{\varphi _{\rm i}^\prime} }}{2} + \dfrac{\pi }{4}$角度，但主剪切裂隙带宽度较小、剪切面较平整、较光滑，如图12(b)所示，这是因为单元弹性模量、应力状态和抗拉/抗剪强度分布较均匀所致。

图  12  不同m值下非均质岩样单轴压缩破裂过程(以m=1.5和m=12为例)

Figure  12.  Fracture process of heterogeneous rock samples under uniaxial compression with different m values (taking m=1.5 and m=12 as examples)

下载: 全尺寸图片 幻灯片

3.3   巴西劈裂

采用图6(b)所示的巴西劈裂数值模型研究不同非均质度下岩样抗拉强度特征，得到如图13所示的模拟结果，可知：巴西劈裂抗拉强度σ_t亦随着m值的增大呈指数升高，这与单轴/三轴抗压强度(图9(b))、弹性模量(图11(a))、等效粘聚力和等效内摩擦角(图9(c))的变化规律相似。

图  13  巴西劈裂抗拉强度σ_t与m值的关系

Figure  13.  Relationship between tensile strength σ_t and m value

下载: 全尺寸图片 幻灯片

3.4   随机参数赋值法可靠性分析

由3.1节~3.3节得到了等效粘聚力、等效内摩擦角、弹性模量、单轴抗压强度和抗拉强度与非均质度m的关系，模拟结果与已有研究采用其他数值模拟方法得到的结果^[22]相吻合，然而，FDEM数值方法有其特殊性，其中之一便是对网格尺寸的依赖性，因为裂隙扩展只能沿着三角形单元边界发生。因此，对于任意FDEM数值模型而言，应尽可能将网格划分足够小以便为裂隙尖端提供足够多的下一步扩展路径，也为更准确反映如图1(c)所示裂隙尖端的应力状态。然而，受计算效率限制，网格尺寸无法划分至无限小，当然也是不必要的，基于此，笔者提出了均质岩体实验室尺度和工程尺度模拟网格尺寸上限值取值方法^[58]，然而该网格尺寸是否也适用于非均质岩体呢？也即当网格尺寸发生变化时，基于随机参数赋值法的非均质岩样模拟结果是否仍保持稳定需要进一步研究。本节将研究单轴压缩模拟结果与不同m值和不同网格尺寸h的关系。

以m=2、3、4、6、10为例，将网格尺寸h分别设为1.3 mm、1.6 mm、1.9 mm、2.2 mm、2.5 mm、2.8 mm、3.1 mm、3.4 mm、3.7 mm和4.0 mm，试样尺寸和加载速率与图6(a)相同，得到如图14所示的模拟结果。一方面，对于不同非均质度的岩样(即m值不同)，单轴抗压强度σ_c在h=1.3 mm、1.6 mm和1.9 mm时均保持稳定，而当网格尺寸大于2.2 mm (对于m=2、3、4的非均质度岩样)或大于2.8 mm (对于m=6、10的非均质度岩样)时，σ_c值将产生震荡、不再保持稳定值，如图14(a)所示，即，即使对于极不均质的岩样(如m=2时)，当网格尺寸足够小时，其模拟强度也能保持稳定，更不必说对于较为均质的岩样，或即，在均质岩体中所适用的网格尺寸范围能够直接用于非均质岩体的模拟；另一方面，在适用网格尺寸范围内，不同m值所模拟得到的岩样均大致沿着各自的破坏角($\dfrac{{{\varphi _{\rm i}^\prime} }}{2} + \dfrac{\pi }{4}$)发生主剪切破坏并伴随少量拉伸裂隙，如图14(b)所示。因此，本节结果表明，随机参数赋值法模拟非均质岩样力学特性和破坏特征是可靠的，其适用网格尺寸范围可采用对应的均质岩样标定得到的网格尺寸。

图  14  不同m值、不同网格尺寸单轴压缩模拟结果

Figure  14.  Simulation results of uniaxial compression with different m values and different element sizes

下载: 全尺寸图片 幻灯片

4   非均质隧洞围岩破裂碎胀大变形模拟

4.1   数值模型及隧洞开挖模拟方法

隧洞开挖数值模型如图15所示，以圆形隧洞为例，直径4 m，网格细化区直径32 m；在远场区内采用较大的网格尺寸，逐渐过渡至模型边界，其网格尺寸为15 m，模型直径为100 m。以m=3 为例，非均质围岩弹性模量分布云图如图15(b)所示。

图  15  隧洞开挖模拟数值模型

Figure  15.  Numerical model of tunnel excavation

下载: 全尺寸图片 幻灯片

研究表明^[63]，掌子面会对前后方岩体产生径向支撑作用，并由掌子面前方一定距离内向掌子面后方逐渐递减。因此，在二维平面应变条件下，为了模拟这种径向支撑效应的逐渐减弱过程，可采用核心材料软化法^[64]，即逐步软化隧洞材料的弹性模量及与弹性模量有关的粘滞阻尼($\mu = 2h\sqrt {E\rho }$)，其具体模拟过程和软化路径的选取可见笔者的研究^[64]。

隧洞开挖的模拟可分为2个步骤，即地应力施加阶段和隧洞材料软化卸荷阶段。在地应力施加阶段，根据所需地应力，采用式(11)得到所有节点的节点力，随后将节点力反向施加至对应节点上。

式中：f_x0、f_y0、f_x1、f_y1、f_x2、f_y2分别为节点0、1、2在x轴和y轴方向的节点力；x₀、y₀、x₁、y₁、x₂、y₂分别为节点0、1、2的坐标；σ_xx、σ_yy、τ_xy、τ_yx分别为x方向和y方向正应力及剪应力，在水平和垂直地应力条件下：σ_xx=σ_h、σ_yy=σ_v、τ_xy=τ_yx=0。

在地应力施加阶段，模型边界自由，仅采用三角形单元、不插入四边形节理单元；施加地应力后，模型产生巨大动能，在粘滞阻尼和临界迟滞阻尼^[65]作用下，迅速达到平衡状态，即模型总动能达到极小值(如小于0.01 kJ^[64])，地应力施加完毕。地应力施加完毕后，插入四边形节理单元、固定边界以保持住得到的地应力。由于四边形节理单元的插入，相邻三角形单元间会发生微小的相互嵌入现象产生较大动能，在粘滞阻尼和临界迟滞阻尼作用下，模型再次达到平衡状态(总动能小于0.01 kJ)，随后即可开始隧洞材料软化卸荷模拟。采用表1的输入参数，在开挖模拟阶段计算时步Δt=3×10⁻⁸ s，而在地应力加载阶段计算时步Δt=2×10⁻⁷ s，这是因为在开挖模拟阶段，插入了四边形节理单元，且围岩破裂后，存在三角形单元间的接触应力计算，需采用较小的计算时步。

4.2   非均质围岩破裂碎胀大变形

施加38 MPa的静水应力场(即侧压系数λ=1)，以m=3为例，得到如图16~图16模拟结果。在约49万步地应力加载完毕，模型达到初次平衡状态；在约82万步达到二次平衡状态。随着核心材料的逐步软化直至隧洞材料完全移除，围岩裂隙扩展过程如图16(a)~图16(f)所示，可知，围岩主要以共轭剪切裂隙形式向围岩深处扩展，并伴随少量拉伸裂隙，这种裂隙形态与顾金才等^[66]在室内隧洞开挖模型试验得到的极为相似。实际上，顾金才等^[66]采用水泥砂浆拌和而成的相似模型，其围岩也具有非均质特征。此外，尽管围岩是非均质性的，但在静水压力条件下，围岩最终损伤破裂区形态依然呈圆形，与均质围岩^[52]的损伤破裂区形态是相似的，半径约为6.6 m，如图16(f)所示。

图  16  非均质隧洞围岩裂隙扩展过程及围岩位移场(m=3)

Figure  16.  Fracture development process and displacement field of heterogeneous surrounding rock mass (m=3)

下载: 全尺寸图片 幻灯片

图  17  非均质隧洞围岩应力场(m=3)

Figure  17.  Stress field of heterogeneous surrounding rock mass (m=3)

下载: 全尺寸图片 幻灯片

本文将隧洞材料完全移除，也即模拟无支护隧洞，因此被共轭剪切裂隙和拉伸裂隙切割的碎裂块体沿着主剪切带发生滑移剪胀扩容，即碎裂块体沿着主剪切带滑移时发生相互不啮合的现象：一方面使得块体间产生宏观裂隙，使其体积膨胀；另一方面，块体在滑移过程中挤压更靠近隧洞周界的浅部碎裂块体，使得浅部块体向隧洞空间发生平移和翻转大运动，造成隧洞断面积的缩小，是为破裂碎胀大变形，如图16(g)所示。在本节工况下，隧道表面最大位移量约达0.6 m，如图16(h)所示。

隧洞围岩径向应力和切向应力场如图17所示，分别为^[64]：

式中：σ_r、σ_θ为极坐标系下的径向应力和切向应力；θ为极坐标与x轴正方向的夹角，以逆时针为正。

尽管对模型施加的地应力为38 MPa，由于弹性模量的非均匀性，单元最大初始地应力可达60 MPa。隧洞开挖后围岩径向应力场分布与均质围岩^[64]相似，即在隧洞周边径向应力降至为0，并向围岩深处逐渐增大至初始地应力；对于切向应力而言，亦在隧洞周界降至为0、在裂隙尖端达到集中状态，最大集中应力约为120 MPa，即为2倍初始地应力。正是由于切向集中应力的作用，围岩发生了共轭剪切破坏。随着浅部岩体的破裂，无法承担巨大切向集中应力致使其向深部完整岩体转移，剪切裂隙也随之向深处扩展，直至在深处随着径向应力(即围压)的增大，岩体强度增高，与切向集中应力达到极限平衡状态，裂隙停止扩展。因此，隧洞围岩裂隙的扩展和隧洞表面变形具有显著的时间和空间效应，即在不同时间段和距掌子面不同位置，围岩裂隙形态、应力场和位移场具有很大差异。

将水平地应力σ_h设为38 MPa、垂直地应力σ_v设为27.1 MPa，即地应力侧压系数λ=1.4，得到如图18所示模拟结果，可知：① 当λ=1.4时，围岩仍以X型共轭剪切裂隙为主，并伴随少量拉伸破裂；② 与λ=1.0相比，当λ=1.4时，围岩仅在顶、底板发生对称性破裂，而两帮围岩保持完整，此外，围岩产生显著的V型松散剥落坑，在V型坑内的碎裂岩体从围岩体中分离，而V型坑外侧仅产生共轭剪切裂隙，未与围岩体发生显著分离，这与诸多实际工程中^[67]观测到的破裂现象极为相似；③ 再者，当λ=1.4时，围岩最大裂隙扩展半径约为4.5 m，小于λ=1.0时的裂隙扩展半径，这是因为围岩产生V型裂坑后易形成稳定结构，最大切向集中应力在裂隙尖端达到平整状态，形成平衡拱结构。

图  18  侧压系数λ=1.4模拟结果

Figure  18.  Simulation results with lateral pressure coefficient λ=1.4

下载: 全尺寸图片 幻灯片

4.3   非均质度m值的影响

将m值分别设为3、4、6、8、10和12，地应力和模型条件与前述相同，不同m值的模拟结果如图19~图20所示(m=3的模拟结果可见4.2节)。从图20可知，对于不同的m值，围岩裂隙形态是相似的，均主要以共轭剪切裂隙的形式发生破裂，并伴随少量拉伸裂隙，且最终的围岩裂隙网络形态均大致呈圆形；即使对于均质围岩，其裂隙形态变化大不，如图19(f)所示。然而，随着m值的增大，即围岩逐渐向均质性靠近，呈现如下规律：

图  19  不同m值的非均质围岩和均质围岩裂隙网络

Figure  19.  Fracture network of heterogeneous surrounding rock mass with different m values and homogeneous rock mass

下载: 全尺寸图片 幻灯片

图  20  不同m值的围岩破裂率和裂隙扩展范围

Figure  20.  Failure degree and fracture propagation range of surrounding rock mass with different m values

下载: 全尺寸图片 幻灯片

1) 定义围岩破裂程度η为破裂四边形节理单元数目与围岩总节理单元数目的比值，即：

式中：n_failed为破裂节理单元数目；N为围岩总节理单元数目。随着m值的增大，围岩破裂程度η呈指数降低，如图20所示。

2) 此外，定义围岩裂隙扩展范围L为围岩损伤破裂区半径与隧洞半径的差值，如图16(f)所示，可知，随着m值的增大，裂隙最大扩展范围L亦呈指数降低，但相比于围岩破裂程度η，裂隙扩展范围随m值增大的降低速率更为缓和，如图20所示，这是因为随着m值的增大，除了裂隙扩展范围变小外，碎裂岩块的尺寸亦有所增大，使得围岩破裂程度下降得更快。

4.4   地应力量值的影响

将静水地应力分别设为34 MPa、36 MPa、38 MPa (4.2节)、40 MPa和42 MPa，非均质度m设为3，其他条件与4.2节相同，得到如图21所示模拟结果，可知：

1) 随着地应力增大，围岩总体裂隙形态基本保持不变，即主要以共轭剪切裂隙形式发生破坏并伴随少量拉伸裂隙，最终裂隙网络形态呈圆形、且变形后的隧洞断面仍大致呈圆形；剪切裂隙沿着主剪切带发生滑移剪胀扩容，且碎裂块体间相互接触推挤向隧洞空间发生挤压碎胀大变形，如图21(a)、图21(b)所示。

2) 此外，随着地应力增大，围岩破裂率呈线性增大，而裂隙扩展范围呈指数函数增大，如图21(c)所示；围岩破裂率增长快于裂隙扩展范围，这是因为随着地应力增大，裂隙网络既向围岩深处扩展，同时浅部块体也碎裂成更为细小的块体。

图  21  不同静水地应力下模拟结果(m=3)

Figure  21.  Simulation results with different in situ stresses (m=3)

下载: 全尺寸图片 幻灯片

5   结论

本文采用有限元-离散元耦合数值模拟方法FDEM提出了基于Weibull分布的非均质岩样随机参数赋值方法，随后研究了非均质岩体在单轴压缩、三轴压缩和巴西劈裂下的力学特征和变形行为，最后对隧洞开挖后非均质围岩的破裂碎胀变形机制进行了研究，结果表明：

(1) 采用基于Weibull分布的随机参数赋值法在统计学上能够表征岩样的非均质特征，采用C语言rand函数、结合Gmsh软件建立的三角形网格和本文提出的FDEM单个三角形单元弹性模量赋值公式，可建立不同弹性模量的三角形单元数目服从Weibull分布、空间分布具有随机性的FDEM数值模型，且四边形节理单元强度(主要为粘聚力和抗拉强度)亦服从该分布；这种方法既适用于实验室岩样尺度，也适用于隧洞开挖等工程尺度。

(2) 单轴压缩、三轴压缩和巴西劈裂模拟结果表明，随着非均质度m值的增大，即逐渐向均质性岩样靠近，单轴抗压强度、三轴抗压强度和巴西劈裂抗拉强度均呈指数增大并向均质岩样的强度值靠近；此外，通过三轴压缩模拟得到的岩样等效粘聚力${c^\prime }$和等效内摩擦角${\varphi _{\rm i}^\prime}$也随着m值的增大呈指数增长，且不断向均质岩样的粘聚力和内摩擦角逼近；在单轴压缩条件下，岩样主剪切裂隙带与水平面呈$\dfrac{{{\varphi _{\rm i}^\prime} }}{2} + \dfrac{\pi }{4}$角度。

(3) 隧洞开挖模拟结果表明，对于非均质围岩，裂隙网络形态与均质围岩保持相似、且m值的影响不大，即均在集中切向应力作用下发生共轭剪切破裂并伴随少量拉伸断裂；被剪切和拉伸裂隙切割的碎裂岩块沿着主剪切带发生滑移剪胀扩容，一方面使得岩块间不再啮合完整，产生宏观空隙，造成碎裂岩块体积膨胀，另一方面深部岩块的滑移剪胀推挤浅部岩块向隧洞空间产生平移和翻转大运动，造成隧洞断面积的减小，是为破裂碎胀大变形。非均质度m值的影响体现在，随着m值的增大，围岩破坏率η和裂隙扩展范围L呈指数降低，且前者降低速率更快，因为随着m的增大，除了裂隙扩展范围有所减缓外，岩块尺寸亦在增大。

(4) 随机参数赋值法的关键是如何获取合适的m值，根据模拟结果可知，在单轴压缩状态下，不同m值呈现显著不同的峰后脆延特性，因此可根据单轴压缩应力-应变曲线峰后特性并配合峰值抗压强度获取合适的m值。此外，本文研究对象为诸如花岗岩等非均质岩石，而类似鹅卵石-砂土等隧道的非均质地层^[68]则不在本文考虑范围之内。