非均质性对岩石宏观力学特性的影响机制

夏海城，邬爱清，卢 波，徐栋栋

(长江科学院 水利部岩土力学与工程重点实验室，武汉 430010)

1 研究背景

岩石材料是由矿物颗粒、胶结材料和空隙等组成的具有一定缺陷的不均匀物质，非均质性是岩石材料的本质特征。这种细观结构的不均匀性导致岩石材料的弹性模量、强度、泊松比等宏观力学参量具有一定的随机性，同时也直接影响受载岩石的破裂演化过程。因此研究非均质性对岩石材料力学性质的影响机制具有重要意义。

国内外学者对岩石非均质问题展开了广泛的研究，并取得诸多有价值的成果。岩石材料物理参数分布的不均匀性是其非均质性的重要方面，胡小荣等[1]应用地质统计学方法对有限单元体力学参数的赋值进行了研究；陈永强[2]采用不同的统计分布模拟材料的非均匀性，得到了不同特性的材料破坏过程应力-应变曲线；康健等[3]研究了非均质随机分布的热膨胀系数和热传导系数对岩石热破裂的影响；冯增朝等[4-5]研究了岩石细胞单元特性及其非均质分布对岩石破坏全过程曲线的影响，岩石亚细观组构的非均质性及其分布规律是决定岩石单轴压缩全过程曲线的重要因素；张征等[6]分析了岩土参数不确定性的主要原因，并研究了岩土参数空间变异性分析的原理、方法和步骤。唐春安等[7-8]、Kaiser等[9]用Weibull分布函数模拟岩石材料参数的非均质性，研究了二维岩石模型脆性破坏的特点，提出“RFPA(Rock Failure Process Analysis)方法”。梁正召[10]开发了RFPA3D程序，模拟了岩石的破坏过程。尤名庆等[11]讨论了岩石非均质性和强度尺寸效应之间的关系；罗荣等[12-13]研究了基于材料组分的非均质岩土体数值建模方法及非均质度对岩石力学特性的影响；柯长仁等[14]将非均质性引入虚内键模型，模拟了岩石式样动态破坏过程；Cil和Buscarnera[15]利用Weibull 分布的最弱环节理论，分析了不同粒径尺寸对单轴最大抗压强度的影响。总的来说，在岩石非均质的模拟方面的研究都采用基于某种统计分布对每个岩石单元赋予各不相同的物理参数，使细观单元的力学性质各不相同，以模拟岩石宏观上的非均质性。

非均质岩石材料力学特性的影响因素不仅在于岩石细观材料参数之间的差异，这种差异形成的力学薄弱环节在空间的分布特征直接控制着受载岩石的破裂演化过程，从而影响岩石的强度和应力-应变全过程曲线的特征。上述研究主要关注因岩石成分不同导致的细观材料参数的差异对岩石宏观力学性状的影响。事实上，因成岩机制和成岩过程的差异，岩石的成分、结构和构造在空间上呈现不均匀分布特征；岩石的弹性模型、强度等宏观物理力学参量的离散性则是这种细观不均匀性的综合体现。既有研究成果未充分关注细观不均匀性的空间分布特征对岩石宏观力学特性的影响。本文采用具有伺服控制机制的非连续变形分析(Discontinuous Deformation Analysis,DDA)方法，根据岩石颗粒呈不规则多边形的特点，采用Voronoi结构来模拟岩石矿物颗粒，构建用于DDA分析的岩石数值试样，通过对不同的Voronoi块体单元赋予不同的物理力学参数的方式来实现对非均质性的模拟。通过开展岩石单轴压缩数值实验，在验证非均质度对岩石强度影响的基础上，探究了岩石单元空间分布的影响作用，补充了主要结论。

2 基于DDA的细观数值实验方法及数值模型构建

2.1 DDA基本理论

DDA方法是石根华博士于20世纪80年代提出的平行于有限单元法的数值模拟方法，其单元或块体形状可以是任何凸状形或凹状形的多边形。[16]相比有限单元法只能用标准形状单元具有一定的优越性，在DDA方法中，任意块体的位移Di可用6个位移变量表示，即

Di=(u0v0r0εxεyγxy) ；

(1)

式中：(u0v0)是块体内特殊点(x0y0)的刚体位移；角r0是块体绕转动中心(u0v0)的转动角，用弧度给出；(εxεyγxy)是该块体的法向和切向应变。

块体中任意点(x，y)处的位移可由变形变量[Di]表示为

(2)

[Ti]=

(3)

式中[Ti]为位移转换矩阵。

各连接块体依靠接触和位移约束形成块体系统。n个块体的联立方程为

式中Kij是6×6子矩阵，与材料特性有关；[Di]和[Fi]是6×1子矩阵，此处Di代表块体i的变形变量，Fi为块体i上的6个变形变量的荷载。对于平衡方程(式(4))，由力和应力产生的总势能Π最小化来推导，对Π求导，得出的总势能最小。子矩阵Kij为

(5)

子矩阵Fi是

(6)

式中d表示块体的位移变量。

在对式(2)和式(4)求解后，可计算出块体变形后的位移与应变。由于DDA法中块体无嵌入和张拉的条件，需在接触位置加上或去掉刚性弹簧，修正总体平衡方程，并再次求解直至所有边界满足无嵌入和张拉，迭代完成。

杨懋偲[17]在DDA源代码的基础上进行了二次开发，初步实现了对电液伺服控制机制的模拟，得到了符合室内常规实验结果的应力-应变全过程曲线。

本文将使用具有伺服控制机制的DDA程序开展非均质度对岩石材料宏观力学性质的影响机制研究。

2.2 正态分布函数及非均质度

岩石材料物理参数的不均匀分布是引起岩石非线性变形的重要因素，在细观数值实验中，将岩石试样划分为多个Voronoi单元，每个Voronoi单元具有不同的物理力学参数，这种细观单元物理力学参数的不均匀分布可以较好地模拟出岩石非均质性的特点。本文采用正态分布函数模拟矿物颗粒物理参数分布的随机性，并定义了量化岩石非均质度大小的方法。通过调整岩石试样非均质度和期望的大小，比较分析了岩石力学性能的变化规律，并探究岩石力学性能的影响因素。

构成岩石的矿物颗粒物理力学参数尽管呈现出强烈的随机性特点，但仍然满足一定的统计学分布规律。根据以往的研究成果，采用不同的参数分布模型开展数值实验得到的破坏过程曲线也不一样。目前，常用于模拟岩石非均质性的随机函数模型有Weibull分布模型和正态分布模型，其中正态分布模型得到的曲线更加符合岩石脆性破坏的特点[2]。本文采用正态分布函数对岩石非均质进行表征。正态分布函数的概率密度函数如式(7)，其中μ为期望，反映了非均质岩石材料的综合性能；σ为方差，反映了岩石材料弹性模量的偏离值。

(7)

将服从正态分布的随机数赋给岩石单元的弹性模量，以模拟岩石材料的非均质性。正态分布函数中，标准差反映的岩石弹性模量偏离值是度量岩石非均质性的重要参数。为准确定义岩石的非均质度，有效量化非均质度的大小，本文将岩石样本中弹性模量的标准差与期望之比Cv定义为变异系数来度量岩石非均质度的大小，即

(8)

式中：变异系数Cv无量纲，反映了岩石单元随机物理参数的偏离度，用来度量岩石的非均质性；Ed为岩石单元弹性模量的标准差；Em为岩石单元弹性模量的平均值。

生成服从正态分布的随机数组时，通过调整标准差和期望的大小来控制变异系数。显然在Cv=0时，岩石为完全均质材料，各细胞单元具有完全相同的物理参数。随着Cv增大，岩石材料弹性模量的非均质度增大，单元弹性模量离散度增强，岩石均匀性变差。Cv越小，岩石单元的物理参数越趋向于期望，岩石整体的均匀性越好。

2.3 非均质材料数值模型与参数赋值

为聚焦岩石材料非均质度对其力学性能的影响，排除非均匀度以外单元分布结构的干扰，本文构建2组不同弹性模量期望的非均质岩石数值模型(如表1)，使用正态分布函数产生的随机数组对岩石单元弹性模量进行赋值，每种变异系数生成6个不同的随机数组(也即形成6个数值试样)，满足相同的期望和标准差，并将6个试样进行相同条件的数值实验。不同变异系数的岩石样本单元弹性模量统计数目直方图如图1所示。

表1 正态分布模型参数及非均质度Table 1 Parameters of normal distribution modeland degree of heterogeneity

图1 Voronoi单元统计数目直方图Fig.1 Histogram of statistical number of Voronoielement of elastic modulus

图2 非均质岩石材料模型Fig.2 Heterogeneous rockmaterial model

图2为本文所采用尺寸为50 cm×100 cm的平面数值模型，为真实模拟岩石细胞单元的形状结构，采用Voronoi多边形细观单元形成的块体系统进行模拟，每个Voronoi单元的面积约为5 cm2，共计975个单元。

本文中数值实验为单轴压缩实验，对刚性加载板采用施加固定位移的方法对岩石模型进行压缩，每个时步的加载位移为5×10-8。在单轴压缩实验中岩石的非均质性主要是由岩石Voronoi单元之间物理参数不均匀引起的，为排除其他物理参数干扰的影响，本文所有岩石Voronoi单元泊松比为0.3，Voronoi单元之间的界面强度设置为统一值。

2.4 非均质岩石材料破坏过程曲线

本文中数值加载模型主要可分为刚性加载板、非均质岩石模型和刚性支架3部分，使用混合高阶DDA数值实验平台，通过对刚性加载板施加竖向位移进行伺服控制单轴压缩数值实验。

图3 全过程应力-应变曲线Fig.3 Complete stress-strain curve

图3为岩石样本应力-应变全过程关系曲线，岩石样本轴向应变通过加载板的位移计算，加载板在轴向的位移量在数值上等于岩石试样的变形量。岩石块体轴向应力随着加载板的不断加载而增大，并在到达极限荷载后承载力迅速下降，下降到某一值后保留一定的残余强度，符合真实的脆性岩石的承载规律。显然岩石变形破坏曲线中应力最大的点也是岩石极限承载力的点，加载板施加的荷载力最大时，岩石应力达到最大。岩石的宏观弹性模量为达到峰值应力50%时，应力与对应应变的切线模量。

图4 均质岩石应力-应变全过程曲线Fig.4 Complete stress-strain curve of homo-geneous rock

3 岩石材料强度非均质度对其力学性能的影响

岩石样本为完全均质时，应力-应变全过程曲线如图4所示，此时所有Voronoi单元的弹性模量均为30 GPa。以此为参照，研究材料非均质度对岩石宏观弹性模量和极限承载力的影响，对每组变异系数的岩石试样均开展数值模拟单轴压缩实验，取岩石样本全过程破坏曲线较为平均的一组对比。

图5为弹性模量期望为30、50 GPa时不同变异系数下的应力-应变曲线。由图5可知，随着变异系数的增大，岩石的极限承载力逐渐下降，峰后曲线具有明显的差异。对于岩石的宏观弹性模量，初始加载阶段不同变异系数的岩石样本宏观弹性模量几乎相同，随着加载的不断进行，变异系数较大的岩石宏观弹性模量先出现下降，总体来看非均质度对岩石样本宏观弹性模量的影响是较弱的。

图5 不同非均质度模型应力-应变曲线Fig.5 Stress-strain curves of rock models ofdifferent degrees of heterogeneity

为了有效研究非均质度对岩石极限承载力的影响，消除相同变异系数下不同单元分布对变形曲线的干扰，聚焦岩石极限应力随非均质度的变化趋势，对每个变异系数下6组弹性模量单元的数值实验结果取平均值，得到平均极限应力，并做出变异系数-极限应力的折线图，结果如图6所示。可以看出，服从正态分布的岩石单元弹性模量，期望值分别为30 GPa和50 GPa时，岩石极限承载力均随着材料变异系数的增大而减小，两者具有较好的线性关系。

图6 变异系数与岩石极限应力的关系Fig.6 Relation between coefficient of variation andultimate stress of rock

值得注意的是，满足相同变异系数和弹性模量期望的不同岩石样本的应力-应变曲线并不完全相同(见图7)，在峰值和峰后阶段具有较大差异。在给岩石单元的弹性模量赋值时，服从某种正态分布的单元集合具有多种分布的可能性，这导致了在相同变异系数的岩石样本中，每个样本内部单元的空间排列方式都是不同的，其岩石的破裂路径也因此存在差异。所以，岩石材料的极限应力不仅受到材料物理参数非均质度的影响，样本中岩石单元的排列方式即细观组构也对岩石力学特性存在重要作用。

图7 弹性模量期望50 GPa和变异系数0.2时不同组构的应力-应变曲线Fig.7 Stress-strain curves of rock material of differentfabrics with a modulus of elasticity of 50 GPa anda coefficient of variation of 0.2

4 材料空间组构对岩石力学性质的影响

“组构”是描述岩石中几何性与物理性要素的重要概念，岩石内所有统计单元以某种排列方式遍布整个区域，其排列方式形成某种组构。理论上讲，满足同一期望和标准差的弹性模量单元集合有无数个，都能形成不同的组构。岩石材料的非均质性不仅在于其物理参数的不均匀，也在于材料组构的随机性排列。如图7所示，不同组构的岩石样本在单轴压缩荷载时有不同的应力-应变曲线。

岩石内部不同矿物颗粒的弹性模量和承载能力存在差异。显然，如果2个接触单元的弹性模量和承载力接近，2单元的应力和应变应基本保持一致，稳定与失稳也会同步；如果2相邻单元弹性模量差异较大，2单元的应变大小会存在较大差异，接触面成为容易发生破坏的“力学薄弱环节”。在岩石样本中，因不同矿物成分力学性状差异形成的力学薄弱环节，因矿物组构不同而呈现随机分布的特征，也即在相同的非均质度下，不同的数值试样“力学薄弱环节”分布也不尽相同。这些薄弱环节往往是岩石压缩过程中率先开裂和破坏的面，其过程如图8所示。

图8 非均质岩石单元接触面剪切过程示意图Fig.8 Schematic diagram of the shearing process ofthe contact surface of heterogeneous rock element

由图8可知，率先失去强度的单元会转移原来承担的部分荷载，周围单元承担的荷载则相应地增大，增幅取决于单元自身与失稳单元弹性模量的差值。很显然，相邻单元之间的弹性模量差异越大，应力梯度则越大，越容易发生开裂破坏。因此，具有相同非均质度而细观组构不同的数值试样，力学薄弱环节分布的不同会造成不同的应力分异格局，从而导致岩石破裂路径不同。

图9 单元弹性模量分布柱状图Fig.9 Histogram of elastic modulus distribution

图10 非均质度为0.5时不同组构的岩石样本破裂过程Fig.10 Rupture process of rock specimens of differentfabrics with a degree of heterogeneity of 0.5

图9展示了岩石非均质度为0.5时Voronoi单元的弹性模量在岩石样本中的空间分布，每个柱体的高度值代表了Voronoi单元弹性模量大小，显然非均质岩石单元的“高度值”是大小不一的。非均质度为0.5、岩石弹性模量期望为30 GPa时，2个不同组构的样本裂纹萌生与扩展路径不同(如图10所示)。该现象反映出同等非均质度不同组构的样本单元分布结构会对岩石破裂路径会产生干扰，使裂纹扩展过程受到力学薄弱环节分布结构的影响，进而影响岩石的强度力学特性。弹性模量较小的单元接触面易成为岩石中的薄弱环节，笔者筛选出弹性模量<8 GPa的单元位置，并绘出其分布图(如图11)，将其与岩石样本破裂时裂纹分布图(如图12)对比，两者具有较好的对应关系。因此，岩石样本的极限承载力与岩石弹性模量单元的空间分布的组构特征存在密切联系。

图11 岩石样本弹性模量<8 GPa的单元分布Fig.11 Element distribution of rock specimen withelastic modulus smaller than 8 GPa

图12 岩石破坏时裂隙密集区域分布Fig.12 Fracture distribution when rock is destroyed

5 结 论

(1)本文使用DDA方法对岩石样本进行数值建模与计算，将岩石样本划分为Voronoi结构单元，用服从正态分布的随机数赋值给单元的弹性模量，将刚性加载板采用固定加载位移的方法进行单轴压缩实验，最终得到非均值脆性岩石材料的单轴压缩应力-应变全过程曲线。

(2)标准差是一个数据集合离散程度的反映，用岩石单元弹性模量的标准差与期望之比(变异系数)描述岩石强度的非均质性是合理的：比值越大，表明岩石强度的非均质度越大，反之岩石强度的非均质度越小。

(3)岩石材料的变异系数和岩石极限承载力之间存在较好的线性关系：随着变异系数的增大，岩石极限承载力减小。岩石的宏观弹性模量同变异系数之间没有较为明显的线性变化特征。

(4)物理参数的非均质性并非岩石强度的唯一影响因素，岩石的组构方式也对岩石强度存在重要影响作用。样本中薄弱环节不同的分布结构会改变岩石裂纹的萌生与扩展路径，从而导致相同非均质度的岩石力学强度性质有所不同。

(5)岩石试样受荷载作用由变形到破坏，是连续小变形逐渐发展成宏观断裂的非连续大变形的过程。传统基于连续介质的数值方法(如FEM)受到自身理论框架的限制，在处理大量颗粒间复杂接触行为时，模拟效果与计算准确度不如DDA方法。DDA方法由于具有完备的运动学理论和严格的接触算法，可以更加正确地描述颗粒之间的接触力学行为。但是DDA方法仍存在一些不足，如对某些物理参数取值较为敏感，大规模计算迭代收敛慢的问题。