2016, Vol. 39
2. 四川省地质调查院, 成都 610081;
3. 北京理工大学 机电工程学院, 北京 100081
2. Geology Investigation Institute of Sichuan Province, Chengdu 610081, P.R.China;
3. School of Mechatronic Engineering, Beijing Institute of Technology, Beijing 100081, P.R.China
岩石通常可以被认为是一种非线性、非均匀性、各向异性的脆性材料。在细观上，岩石内部具有不同的弹性模量、泊松比、破坏强度等，从而具有非均匀特性[1]。这种细观上的非均匀特性将影响到岩石材料的宏观力学响应和破坏形式，因为岩石材料的宏观力学特征是其内部细观力学特性的综合反映。在数值模拟中若要更好的再现岩石的复杂力学行为，就需要引入非均匀性的概念。
统计学方法已经被广泛用于研究脆性材料的非均匀特性[2]。如Danzer等[3]的实验研究结果表明，陶瓷的破坏统计符合Weibull分布，Gulino等[4]采用Weibull分布函数研究了石墨纤维的非均匀性，Van Mier等[5]在颗粒复合材料的断裂机制研究中则分别采用了高斯和Weibull分布。一些模拟软件或方法已经可以模拟岩石的非均匀性，比如基于FEM(finite element method, FEM)方法的RFPA(rock failure proeess analysis, RFPA)软件[6]，以及DEM(distinct element method, DEM)[7]和DDA(discontinuous deformation analysis, DDA)[8]方法等。Tang等[6]在RFPA中提出了有关岩石细观单元强度满足Weibull分布的假设，结果表明，假定岩石材料单元的力学性质满足Weibull分布。朱万成等[9]则在RFPA的基础上，利用细观非均匀单元对混凝土的断裂进行了模拟分析，并提出采用最大拉应变破坏准则进行开裂判断。Liu等[10-11]则进一步基于RFPA模拟了岩石的非均匀性，分析了非均匀指数对岩石开裂过程以及开裂强度的影响。此外，Chen等[7]在DEM方法中引入了非均匀模型，并对岩石的力学行为进行了有效的分析；王士民等[12]采用演化元胞自动机方法，数值模拟了Weibull坡度参数表示的非均匀性对脆性岩石材料破坏形式的影响；Fang等[13]则通过采用Weibull函数描述岩石强度的非均匀特性，用FLAC模拟了矿柱系统的破坏过程；唐欣薇等[14-15]还将常规Weibull分布模型进行了改进，引入空间相关尺度因子，分别对岩石试件以及混凝土的细观力学模型进行了分析。Jiao等[8]在DDA方法中也考虑了岩石的非均匀性，但仅是对弹性模量和泊松比等弹性参数进行了非均匀化，并讨论其对岩石开裂的影响，而没有考虑强度参数的非均匀特性。
笔者采用非连续变形分析(DDA)方法，基于Weibull分布函数引入非均匀模型，考虑弹性参数(弹性模量、泊松比)以及强度参数(抗拉强度、内聚力、内摩擦角)的非均匀特性，计算模拟了巴西圆盘的劈裂实验，基于模拟结果得到了一些新的结论。
1 数值模拟模型 1.1 非均匀模型笔者所用的非连续变形分析(DDA)方法是对离散可变形块体系统的静态和动态力学行为进行计算的不连续数值计算方法。由于DDA方法中采用的是线弹性本构关系，细观单元以及宏观模拟结果的应力应变关系均为线性。在DDA方法中引入非均匀模型，在采用线性细观单元的条件下，得到呈现出非线性特性的宏观模拟结果，从而更好地模拟岩石的物理力学特性。根据Weibull双参数函数的分布公式
| $p\left( x \right) = \left\{ {\begin{array}{*{20}{c}} {\frac{{m{x^{m- 1}}}}{{x_0^m}}\exp \left[{-{{\left( {\frac{x}{{{x_0}}}} \right)}^m}} \right], x \ge 0}\\ {0, x < 0} \end{array}} \right\}, $ | (1) |
对式(1)在0~x区间上进行积分获得累积分布函数
| $Q\left( x \right) = \int_0^x {P\left( x \right){\rm{d}}x = 1- \exp \left[{-{{\left( {\frac{x}{{{x_0}}}} \right)}^m}} \right], } $ | (2) |
式(2)中令Q(x)=y，简化后可以得到
| $x = {x_0}{\left[{-\ln \left( {1-y} \right)} \right]^{\frac{1}{m}}}, $ | (3) |
式(3)中：y为(0, 1)之间的随机数，由随机数函数产生；m是非均匀指数控制参数；x0是均匀模型材料参数；x是非均匀化的模型材料参数。只需要将弹性模量E，泊松比μ，内摩擦角θ，内聚力c，抗拉强度σt分别作为x0代入，即可获得这些参数非均匀化后的x。这样就能对均匀模型材料参数进行非均匀化，并应用于程序计算中。
1.2 巴西圆盘劈裂模型采用引入非均匀模型的DDA程序对巴西圆盘的劈裂实验进行数值模拟。巴西圆盘的几何模型如图 1所示，圆盘半径为50 mm，上下两端各削一个10°的平台，以减小端部的应力集中。在圆盘上下两端均采用准静态线性位移加载，加载速率为2 mm/s。为了测量圆盘的宏观应力应变关系，在圆盘上端安放2块加载板，上板布置加载点，下板中布置一个测点。在采用一阶位移函数的DDA方法中，单个块体具有常应力，因此，测点的Y方向应力能够反映圆盘在竖直方向上的平均应力。圆盘采用三角形单元进行子块体单元划分，共划分子块体单元4 633个。在径向压载作用下，圆盘的开裂破坏采用子块体开裂模拟算法[16]进行模拟。
|
图 1 巴西圆盘几何模型 Figure 1 Geometrical model of a Brazilian disk |
假定圆盘为线弹性岩石材料，其弹性模量E=20 GPa，泊松比μ=0.2，密度ρ=0.1 kg/m3。在DDA方法中，静动载条件下均考虑材料的惯性力加载项，这里取较小的材料密度，目的是在静载条件下，减小惯性力加载，加快试件中的应力平衡分布，从而提高计算效率[15]。岩石的强度参数取为:抗拉强度σt=5 MPa，内聚力c=20 MPa，内摩擦角θ=30°。DDA模拟中最大位移比、时间步长和弹簧刚度分别取为5×10-4、2×10-6 s和100 GPa。
2 模拟结果及分析 2.1 非均匀巴西圆盘和均匀巴西圆盘的模拟结果为了研究非均匀性对模拟结果的影响，进行了一系列的巴西圆盘劈裂实验的DDA模拟。算例采取不同的非均匀指数(m=1.5, 2.5, 5, 200)进行了4组非均匀巴西圆盘劈裂数值模拟。模型中同时考虑了弹性参数(弹性模量、泊松比)以及强度参数(抗拉强度、内聚力、内摩擦角)的非均匀性(定义为“参数条件1”)，并将模拟结果与均质圆盘的数值模拟结果进行了对比。图 2是均匀巴西圆盘的模拟结果。图 3给出了不同非均匀指数条件下的非均匀巴西圆盘的DDA模拟结果。根据图 2、图 3可以发现，随着非均匀指数值的增大，即圆盘非均匀性的减弱，非均匀巴西圆盘的模拟结果越来越接近均匀巴西圆盘的模拟结果。具体而言，当m=1.5时，裂纹产生的随机性强，裂纹较分散，圆盘的非均匀性明显，在圆盘中心形成了一条较宽的裂纹带；随着非均匀指数m的增大，裂纹分布的随机性明显减弱，次生裂纹明显减少，主裂纹位于圆盘径向的程度更加明显。当m=200时，非均匀巴西圆盘的模拟结果与均匀巴西圆盘的模拟结果几乎一致。
|
图 2 均匀巴西圆盘的DDA模拟结果 Figure 2 DDA Simulation result of the homogeneous Brazilian disk |
|
图 3 不同非均匀指数的非均匀巴西圆盘的DDA模拟结果 Figure 3 DDA Simulation results of heterogeneous Brazilian disks with different heterogeneity index (参数条件1：同时考虑弹性、强度参数的非均匀性) (Condition 1: considering heterogeneities of the elastic and strength parameters) |
图 4是不同非均匀指数巴西圆盘的宏观应力位移曲线，其中Y方向应力是指通过测点测得的圆盘Y方向的平均应力。结果表明，对于较小的m值，由于其非均匀性更强，存在着很弱的块体单元，在加载的早期就出现了开裂；而较大的m值，非均匀性较弱，每个块体单元的材料强度差异并不大，从而表现为整体材料的平均强度较高，在加载过程中出现开裂也较慢。另一方面，较小m值的应力位移曲线呈现出非线性，是因为这种情况下材料强度梯度很大，在加载初期，那些强度较弱的单元率先发生了破坏，产生局部裂纹；随着加载步的进行，这些微小的裂纹逐渐贯通、汇聚形成宏观裂纹。正是这些微小裂纹的逐渐贯通过程使得整个圆盘的力学响应出现了非线性，这比较符合岩石材料真正的力学特性。
|
图 4 不同非均匀程度的巴西圆盘的应力位移曲线 Figure 4 Stress-displacement curves of Brazilian disks under different heterogeneity degree (参数条件1：同时考虑弹性、强度参数的非均匀性) (Condition 1: considering heterogeneities of the elastic and strength parameters) |
在DDA计算中，虽然细观单元是线性的，但是由于模型材料具有非均匀性，因而呈现出宏观的非线性。DDA模拟反映了岩石的非线性、非均匀性的特点。在非均匀的模型中，m值的取值对材料的非均匀性和宏观强度影响很大。m值越小非均匀性越明显，岩石的宏观强度越低，即非均匀性的引入使得岩石表现出更低的宏观强度特性。
2.2 强度参数的非均匀性对模拟结果的影响在DDA方法中，由于以往的工作仅针对弹性参数(弹性模量、泊松比)的非均匀性进行了研究[5]，笔者则将强度参数(抗拉强度、内聚力、内摩擦角)也同时进行了非均匀化，认为这样更加接近于真实的岩石非均匀性特性。图 5给出了仅考虑弹性参数非均匀性(定义为“参数条件2”)的巴西圆盘劈裂破坏模拟结果。与图 3所示的同时考虑弹性、强度参数非均匀性(参数条件1)的巴西圆盘劈裂破坏模拟结果相比可以看出，当m值相同时，参数条件1得到的裂纹更为分散，这种现象在m值越小的情况下越为显著，而随着m值的增大，2种参数条件下的模拟结果都逐渐接近均匀试件的破坏结果，因此当m值较大时2种参数条件得到的破坏结果之间的差异变得不明显。
|
图 5 不同非均匀指数的非均匀巴西圆盘的DDA模拟结果 Figure 5 DDA simulation results of heterogeneous Brazilian disks with different heterogeneity index (参数条件2：仅考虑弹性参数的非均匀性) (Conditon 2: only considering heterogeneity of the elastic parameters) |
图 6给出了m=1.5和m=200圆盘的宏观应力位移曲线。可以看出，当m值相同时，相比于参数条件2，参数条件1情况下圆盘中的裂纹在更早的加载阶段产生，圆盘表现出的宏观强度更低，宏观非线性也更明显。同时，m值越大，2种参数条件下得到的宏观强度相差越小，即非均匀性带来的影响越不明显。
|
图 6 考虑不同参数的非均匀性的巴西圆盘的应力位移曲线 Figure 6 Stress-displacement curves of Brazilian disks considering heterogeneities of different parameters |
通过引入基于Weibull函数得到的非均匀模型到DDA方法中，并对非均匀巴西圆盘劈裂实验的数值模拟结果进行分析，可以得出如下结论：
1) DDA方法的功能得到了改进，该算法实现了对岩石材料的非均匀特性的模拟功能；
2) 非均匀模型的引入，实现了采用细观的线性单元得到岩石的宏观非线性力学响应。非均匀指数m值的大小对岩石材料的非均匀性和宏观强度有重要影响，m值越小意味着非均匀性越明显，岩石材料的宏观强度越低，宏观应力应变关系(即宏观应力位移关系，如图 4所示)的非线性也越明显；
3) 在对弹性模量，泊松比进行非均匀化的基础上，添加对强度参数(抗拉强度，内聚力，内摩擦角)的非均匀化，可以使得模拟结果更接近岩石的真实非均匀和非线性特性。
| [1] |
张永兴, 贺永年. 岩石力学[M]. 北京: 中国建筑工业出版社, 2008 .
ZHANG Yongxing, HE Yongnian. Rock mechanics[M]. Beijing: China Architecture & Building Press, 2008 . (in Chinese) |
| [2] | Tang S B, Zhang H, Tang C A, et al. Numerical model for the cracking behavior of heterogeneous brittle solids subjected to thermal shock[J]. International Journal of Solids and Structures , 2016, 80 : 520–531. DOI:10.1016/j.ijsolstr.2015.10.012 |
| [3] | Danzer R, Supancic P, Pascual J, et al. Fracture statistics of ceramics-weibull statistics and deviations from weibull statistics[J]. Engineering Fracture Mechanics , 2007, 74 (18) : 2919–2932. DOI:10.1016/j.engfracmech.2006.05.028 |
| [4] | Gulino R, Phoenix S L. Weibull strength statistics for graphite fibres measured from the break progression in a model graphite/glass/epoxy microcomposite[J]. Journal of Materials Science , 1991, 26 : 3107–3118. DOI:10.1007/BF01124848 |
| [5] | VanMier J G, VanVliet M R A, Wang T K. Fracture mechanisms in particle composites: statistical aspects in lattice type analysis[J]. Mechanics of Materials , 2002, 34 (11) : 705–724. DOI:10.1016/S0167-6636(02)00170-9 |
| [6] | Tang C A, Yang W T, Fu Y F, et al. A new approach to numerical method of modelling geological processes and rock engineering problems-continuum to discontinuum and linearity to nonlinearity[J]. Engineering Geology , 1998, 49 : 207–214. DOI:10.1016/S0013-7952(97)00051-3 |
| [7] | Chen W, Konietzky H. Simulation of heterogeneity, creep, damage and lifetime for loaded brittle rocks[J]. Tectonophysics , 2014, 633 : 164–175. DOI:10.1016/j.tecto.2014.06.033 |
| [8] | Jiao Y Y, Zhang X L, Zhao J. Two-dimensional DDA contact constitutive model for simulating rock fragment-ation[J]. Engineering Mechanics , 2012, 138 : 199–209. DOI:10.1061/(ASCE)EM.1943-7889.0000319 |
| [9] |
朱万成, 赵启林, 唐春安, 等.
混凝土断裂过程的力学模型与数值模拟[J]. 力学进展 , 2002, 32 (4) : 579–598.
ZHU Wancheng, ZHAO Qilin, TANG Chunan, et al. Mechanical model and numerical simulation of fracture process of concrete[J]. Advances in Mechanics , 2002, 32 (4) : 579–598. (in Chinese) |
| [10] | Liu H Y, Roquete M, Kou S Q, et al. Characterization of rock heterogeneity and numerical verification[J]. Engineering Geology , 2004, 72 : 89–119. DOI:10.1016/j.enggeo.2003.06.004 |
| [11] | Liu H Y. Numerical modelling of the rock fragmentation process by mechanical tools[D]. Sweden: Lulea University of Technology, 2004. |
| [12] |
王士民, 朱合华, 冯夏庭, 等.
细观非均匀性对脆性岩石材料宏观破坏形式的影响[J]. 岩土力学 , 2006 (2) : 224–227.
WANG Shimin, ZHU Hehua, FENG Xiating, et al. Influence of heterogeneity on macroscopical crack form of the brittle rock[J]. Rock and Soil Mechanics , 2006 (2) : 224–227. (in Chinese) |
| [13] | Fang Z, Harrison J P. Development of a local degradation approach to the modelling of brittle fracture in heterogeneous rocks[J]. International Journal of Rock Mechanics & Mining Sciences , 2002, 39 (4) : 443–457. |
| [14] | Tang X W, Zhou Y D, Zhang C H, et al. Study on the heterogeneity of concrete and its failure behavior using the equivalent probabilistic model[J]. Journal of Materials in Civil Engineering , 2011, 23 (4) : 402–413. DOI:10.1061/(ASCE)MT.1943-5533.0000179 |
| [15] |
唐欣薇, 周元德.
基于改进Weibull分布模型的岩石单轴压缩变形破裂规则研究[J]. 岩石力学与工程学报 , 2014, 33 (1) : 2802–2807.
TANG Xinwei, ZHOU Yuande. Study of rock uniaxial compression deformation and failure character based on an improved Weibull distribution model[J]. Chinese Journal of Rock Mechanics and Engineering , 2014, 33 (1) : 2802–2807. (in Chinese) |
| [16] |
倪克松, 甯尤军.
DDA子块体开裂模拟算法的优化与验证[J]. 地下空间与工程学报 , 2014, 10 (5) : 1017–1022.
NI Kesong, NING Youjun. The optimization and validation of sub-block method for modeling rock fracturing within discontinuous deformation analysis framework[J]. Chinese Journal of Underground Space and Engineering , 2014, 10 (5) : 1017–1022. (in Chinese) |