, 乔磊, 张玉良, 孙建华, 陈华天津大学 土木建筑工程学院, 水利工程仿真与安全国家重点实验室, 天津 300350
收稿日期：2021-01-16
基金项目：国家重点研发计划项目(2018YFC0406804);国家自然科学基金资助项目(11772221)
作者简介：赵高峰(1982-), 男, 教授。E-mail: gaofeng.zhao@tju.edu.cn
摘要：超大规模非均质模型的生成速度制约着离散数值方法在复杂工程实践中的应用，为了解决超大规模模型非均质岩石的快速建立问题，该文基于四维弹簧模型（four-dimensional lattice spring model，4D-LSM），提出了偏心四维弹簧模型（eccentric 4D-LSM，ECC4D），开展了Poisson比、弹性模量和法向变形量的参数分析研究。通过弹性参数和破坏参数的研究论证了模型的可行性。结果表明：采用质心随机偏移的方式生成了非规则偏心四维模型ECC4D，通过特定算法可以轻松实现大规模计算，既能在一定程度上表征岩石材料的非规则性，也能快速获取任意几何模型。该方法在本地工作站上已经成功实现了超过1亿颗粒的大坝模型计算。
关键词：四维弹簧模型(4D-LSM)偏心高性能计算岩石破坏
Eccentric four-dimensional lattice spring model for heterogeneous rock fracturing
ZHAO Gaofeng
, QIAO Lei, ZHANG Yuliang, SUN Jianhua, CHEN HuaState Key Laboratory of Hydraulic Engineering Simulation and Safety, School of Civil Engineering, Tianjin University, Tianjin 300350, China
Abstract: The time needed to generate ultra-large heterogeneous models limits the application of discrete numerical methods in rock engineering. This paper presents an eccentric four-dimensional spring model (eccentric 4D-LSM, ECC4D) based on a 4D-LSM (four-dimensional lattice spring model) for parameter analyses of the effects of Poisson ratio, the elastic modulus and the ultimate normal deformation. The model feasibility is demonstrated through a study of the elastic parameters and failure parameters. The irregular, eccentric four-dimensional model ECC4D is generated using a random offset of the center of mass that can be used for large-scale calculations. ECC4D characterizes the irregularities of the rock materials while quickly generating an arbitrary geometry model. This method has been used to prepare a dam model with more than 100 million particles on a workstation.
Key words: four-dimensional lattice spring model (4D-LSM)eccentrichigh performance computingrock fracturing
岩石的物理力学性质研究对于工程实践具有重要意义，面对工程建设规模日益扩大、地质条件愈加复杂、载荷和边界条件复杂多变的现实情况下，常规的力学分析方法难以获得精确解，计算力学的发展极大地解决了该难题，不仅加快了工程问题的求解速度，还具有成本低、可重复性强、数据全面等优点，工程中已经有较多的应用[1-4]。随着计算机性能的提高，岩石力学数值计算方法得到了极大的发展，一些比较成熟的数值方法如离散元法[5]、有限元法[6]、非连续变形分析方法[7-8]、数值流形法[9]等已经被广泛应用到岩石工程中。同时，一些以岩石力学为出发点的新的数值方法也相继出现。格点弹簧模型(lattice spring model, LSM) 最早由Hrennikoff[10]于1941年提出，但是由于Poisson比限制问题，LSM发展较为缓慢。离散弹簧模型(distinct lattice spring model, DLSM) 是Zhao等[11]提出的一种不同的LSM模型。为了解决Poisson比限制问题，DLSM引入了多体剪切弹簧来计算多体非中心相互作用[12]，该模型已被成功应用于岩石与煤的动态破坏研究[13-15]，但是计算效率仍然受到较大的制约。为了解决计算效率问题，研究者基于经典LSM理论引入了第四维相互作用，不仅突破了Poisson比的限制还可以求解更加复杂的固体动力学问题，称之为四维弹簧模型(4D-LSM)[16]，该模型也成功应用于岩石的动态破坏研究[17-18]。然而岩石的非均质性和各向异性决定了其物理力学性质的复杂性，因此通过构建随机模型反映岩石的非均质性和各项异性是必要的。本文在4D-LSM的基础上将质心随机偏移，快速建立具有非均匀性特点的计算模型，称为偏心四维弹簧模型(eccentric 4D-LSM, ECC4D)。模拟岩石破坏需要大量的数值单元才能反应真实情况，除了计算过程需要大量的计算资源，模型的前处理也一直是离散介质数值方法的薄弱环节[19-21]。本文提出的快速构建超大规模非规则模型对模拟真实条件下的岩石力学行为具有积极意义。
1 偏心四维弹簧模型1.1 数学模型构建偏心四维弹簧模型的构建基于4D-LSM。4D-LSM中颗粒质心和形心是重合的，质心之间由弹簧连接；而ECC4D中，通过随机偏移质心的方式使质心产生移动，质心与形心的位置不一定具有重合关系，如此构建成ECC4D随机模型，二维模型示意图如图 1所示。按照质心随机偏移方式的不同，ECC4D模型可分为2种，记为Type A和Type B。
 |
| 图 1 4D-LSM与ECC4D二维模型区别 |
| 图选项 |
Type A的质心偏移方式如下：在球形颗粒中，质心在以形心为中心、边长为L的立方体中随机移动，该立方体形心与颗粒形心一致，如图 2所示。设颗粒i的形心坐标分别为xiC、yiC和ziC，颗粒半径为R，立方体的边长为L，定义偏心系数Ra=(L/2)/R，则新的质心坐标xiECC、yiECC和ziECC可表示为
| $\left\{\begin{array}{l}x_{i}^{\mathrm{ECC}}=x_{i}^{\mathrm{C}}+(L / 2) \cdot \operatorname{Rand}(-1, 1)= \\x_{i}^{\mathrm{C}}+R \cdot R_{\mathrm{a}} \cdot \operatorname{Rand}(-1, 1), \\y_{i}^{\mathrm{ECC}}=y_{i}^{\mathrm{C}}+(L / 2) \cdot \operatorname{Rand}(-1, 1)= \\y_{i}^{\mathrm{C}}+R \cdot R_{\mathrm{a}} \cdot \operatorname{Rand}(-1, 1), \\y_{i}^{\mathrm{ECC}}=y_{i}^{\mathrm{C}}+(L / 2) \cdot \operatorname{Rand}(-1, 1)= \\y_{i}^{\mathrm{C}}+R \cdot R_{\mathrm{a}} \cdot \operatorname{Rand}(-1, 1)\end{array}\right.$ | (1) |
 |
| 图 2 ECC4D的2种生成方式 |
| 图选项 |
其中，Rand()为随机函数，可随机生成某区间内的一个数，上述3个Rand(-1, 1)均相互独立。
Type B的质心偏移方式如下：在球形颗粒中，质心只能在以形心为中心，半径为r的球面上随机移动，该球面所在球体的形心与颗粒形心一致，如图 1和图 2所示。设颗粒i的形心坐标分别为xiC、yiC和ziC，颗粒半径为R，限制质心移动范围的球体半径为r，定义偏心系数Ra=r/R，则新的质心坐标xiECC、yiECC和ziECC可表示为
| $\left\{\begin{array}{l}x_{i}^{\mathrm{ECC}}=x_{i}^{\mathrm{C}}+r \cdot \sin \theta \cdot \cos \varphi= \\x_{i}^{\mathrm{C}}+R \cdot R_{\mathrm{a}} \cdot \sin \theta \cdot \cos \varphi, \\y_{i}^{\mathrm{ECC}}=y_{i}^{\mathrm{C}}+r \cdot \sin \theta \cdot \sin \varphi= \\y_{i}^{\mathrm{C}}+R \cdot R_{\mathrm{a}} \cdot \sin \theta \cdot \sin \varphi, \\z_{i}^{\mathrm{ECC}}=z_{i}^{\mathrm{C}}+r \cdot \cos \theta=z_{i}^{\mathrm{C}}+R \cdot R_{\mathrm{a}} \cdot \cos \theta\end{array}\right.$ | (2) |
| $\left\{\begin{array}{l}\theta=\arccos (\operatorname{Rand}(0, 1)) \\\varphi=2 \pi \cdot \text { Rand }(-1, 1)\end{array}\right.$ | (3) |
 |
| 图 3 4D-LSM和ECC4D三维模型区别 |
| 图选项 |
2 模型验证分析2.1 偏心系数在提出ECC4D模型时引入了偏心系数Ra，为了研究Ra对几何构型的影响，分别对不同Ra取值下弹簧键的方向分布和长度分布进行统计。考虑Ra=0.1~0.7，间隔为0.2，结果如图 4和5所示。对于弹簧方向分布规律，Type A和Type B这2种模型的变化规律是一致的。随着Ra的增加，弹簧键在球面上以某些点为中心逐渐地朝整个球面扩散；对于弹簧长度分布，2种模型的规律也趋于一致，从间断的双峰分布逐渐过渡到连续的双峰分布最终到单峰分布。图 6表示Ra为0.1、0.5、0.9这3种偏心系数下2种模型的弹簧连接键的结构图，可以看到模型的棱角逐渐减小而趋于“圆润”。通过引入Ra，弹簧键方向和长度从原来的“确定量”逐步转变成了具有一定统计规律的“随机量”，即ECC4D在一定程度上表征了岩石材料的几何非均匀特征。后文若未加说明，则默认的研究对象是Type A方式生成的ECC4D。
 |
| 图 4 ECC4D的2种生成方式弹簧键方向分布随Ra的变化图 |
| 图选项 |
 |
| 图 5 ECC4D的2种生成方式弹簧键长度分布随Ra的变化图 |
| 图选项 |
 |
| 图 6 ECC4D不同的Ra对应的弹簧键结构图 |
| 图选项 |
2.2 Poisson比Poisson比作为材料重要的力学参数，是参数研究的重点对象。在4D-LSM中，通过试验获取的Poisson比μ和数值计算得到的Poisson比μs是相等的，而ECC4D在构建模型的过程中随机偏移了质心，其μs和μ的关系是否吻合，还需要进一步研究验证。
相对于4D-LSM，ECC4D引入了偏心系数Ra=r/R，其中r是限制质心移动生成范围的球体的半径(或者是立方体边长的一半)，R表示颗粒半径。因此，Ra的取值范围为(0, 1)，Ra=0表示质心没有偏移，Ra=1表示质心位于颗粒表面的极限位置。在物理世界，2个物体的相互作用大多数都是随着距离的增加而减小，当距离增加到一定程度，这种相互作用就可以忽略不计了。在ECC4D中同样考虑这种规律，通过引入另一个参数Rc，取值范围同为(0, 1)，称之为截断系数。其具体意义如下：记ECC4D中弹簧键的最大长度为Lmax，那么，当长度超过Rc×Lmax时，弹簧键将被截断。当Rc=1时，Rc×Lmax=Lmax的弹簧键长已经是最大弹簧键长了，因此没有弹簧键被截断；当Rc=0时，Rc×Lmax=0，此时所有的弹簧键都将被截断。
建立一个立方体单轴压缩模型(见图 7)进行Poisson比参数分析，通过测点监测横向和竖向上的变形并计算数值模拟所得μs，选取3种截断系数Rc进行模拟，输入的μ为0.2，研究Ra与μs的关系。如图 8，当截断系数Rc=1.0时，即弹簧键均未被截断，可以看作未引入截断系数。随着Ra的增加，重现μs整体呈现单调下降趋势。当引入截断系数后，曲线的变化趋势不再呈现单调关系，存在峰值点，并且当Rc=0.6时，μs已经突破了μ。其数值模型也发生了相应的变化，连接键从有序逐渐过渡到随机。因此，可以推测ECC4D模型重现的Poisson比可突破0.25的限制。
 |
| 图 7 立方体单轴压缩试验的计算模型 |
| 图选项 |
 |
| 图 8 Ra与μs的关系(μ的输入值为0.2) |
| 图选项 |
为了进一步探究μs峰值变化情况，考虑不同Ra条件下Rc和μs的变化关系。如图 8所示，一方面曲线整体均为先上升后下降的规律，部分曲线已经突破了Poisson比的限制，超过了0.25，最大的重现Poisson比已经达到了0.36。另一方面，不同的μ对应的曲线变化规律是一致的，峰值点的位置与Ra是一一对应的。但是，直接截断弹簧键的处理方式也存在一些问题，图 8中Rc小于0.35的部分没有曲线变化数据，这是由于在偏心系数Ra为0.7和0.8的条件下，截断系数过小会产生计算错误，因此曲线的变化形式存在一定的缺陷。初步分析原因可能由于ECC4D是基于4D-LSM演变而来，在4D-LSM中，颗粒所连接的弹簧数量有限并且长度固定；ECC4D首先通过随机生成质心的方式使得键长在某一范围内成为了随机量，而引入Rc使得“过长”的弹簧被截断，但这也导致某些颗粒的所有弹簧键有可能都因为“过长”而被截断，致使该颗粒处于“无连接状态”，进而使模型本身处于不稳定状态并造成计算错误。
为了解决上述问题，将弹簧键直接截断的简单处理方式变为“弱化”。当过长的弹簧需要截断时，将其弹簧刚度按比例大幅度降低，弱化比例记为γ=kc/k，其中kc为弱化后的弹簧刚度，k为原弹簧刚度。图 9表明输入Poisson比的不同对曲线变化规律没有影响，下文输入Poisson比均采取μ=0.25进行研究。引入弱化系数后，重新探究Rc和μs的关系，结果如图 10和11所示。当Ra=0.0~0.5时，重现Poisson比与截断系数的变化规律表现为先降低后升高，不再是先升高后降低的变化规律。并且重现Poisson比对截断系数的敏感性较大，致使Rc以0.1为间距变化时无法有效地找到其最小值点。当Ra=0.6~0.9时，重现Poisson比和截断系数的变化规律表现为先升高后降低，曲线变化形式较为平缓，对于不同的Ra，μs峰值位置都是Rc=0.5，并且峰值随着Ra的增大有一定的增幅。
 |
| 图 9 Rc与μs的关系 |
| 图选项 |
 |
| 图 10 Rc与μs的关系(Ra=0.0~0.6，γ=0.1，μ=0.25) |
| 图选项 |
 |
| 图 11 Rc与μs的关系(Ra=0.6~0.9，γ=0.1，μ=0.25) |
| 图选项 |
μs的峰值也是一个重要的研究参数，记为μs-max。图 11中可以看出，曲线呈现一定的规律性，即在Rc=0.5处出现峰值且Ra=0.9时峰值最大。因此，相关参数分别取为Ra=0.9、Rc=0.5、μ=0.25，研究弱化系数和重现Poisson比的峰值之间的变化规律，同时也对重现的弹性模量Es进行分析，并且考虑了Type A和Type B两种生成方式的影响，其结果如图 12所示。Es对生成方式的影响以Es/E表示，Es是重现的弹性模量，E是输入的弹性模量。由图 12可知，随着弱化系数的增加，重现Poisson比呈现下降的趋势，重现弹性模量呈现上升的趋势。其原因是γ的值越大，代表弹簧键弱化的程度越小，即弹簧刚度越大，进而使Poisson比减小而弹性模量增大；Type A和Type B的2条μs-max曲线几乎平行，Type A的μs-max稍大，两者的最大值都达到了0.3左右；Type A和Type B的2条Es曲线在γ较小时差别较大，随着γ的增大，两者的差距逐渐减小并且有重合的趋势，同时，由于部分弹簧键被弱化，Es/E的比值始终小于1。
 |
| 图 12 γ对μs-max和Es的影响(Ra =0.9，Rc=0.5，μ=0.25) |
| 图选项 |
由于ECC4D的模型生成是随机的。因此，其弹性参数必然存在一定的波动性，但是ECC4D的目的主要是为了快速的构建非均匀模型并用于大规模计算。因此，需要对ECC4D弹性参数是否能够随着模型规模的扩大而收敛进行研究。模拟依然采用立方体单轴压缩模型，但是改变模型的物理尺寸，边长10 ~ 80 mm，以10 mm递增。其他参数分别为：Ra=0.9，Rc=0.5，γ=0.001，μ=0.25，参数相同的模型随机生成10次并计算出μs和Es，结果分别如图 13和14所示。图 13表明，随着模型规模的增大，μs-max具有很好的收敛性，收敛值大约为0.3；图 14中，Es同样也有收敛的趋势，但是当模型规模较小时，先有一段上升的趋势。
 |
| 图 13 μs-max随模型规模的变化规律(Ra=0.9，Rc=0.5，μ=0.25) |
| 图选项 |
 |
| 图 14 Es随模型规模的变化规律(Ra=0.9，Rc=0.5，μ=0.25) |
| 图选项 |
2.3 破坏参数当采取质心随机偏移后，ECC4D模型获取的破坏参数与4D-LSM是否存在差异，需要进一步讨论。4D-LSM中弹簧键的变形定义为U，当U达到最大的极限值U*时弹簧键即断开，U*是4D-LSM的唯一破坏参数，其值可以用公式U*=Dσt/E进行估算，其中：D为颗粒直径，σt为抗剪强度，E为弹性模量。同样采用单轴拉伸数值实验来研究ECC4D中破坏参数UECC*的变化规律，单轴拉伸模型仍然采用立方体，边长20 mm，与图 7的单轴压缩模型相同，但是荷载方向改为竖直向上，颗粒直径为1 mm，Ra=0.9，Rc=0.5，弱化比例γ=0.001。首先，考虑不同的Ra、Rc和γ，输入的破坏参数为UECC*=0.001~0.001 6，间隔0.000 1，分别进行单轴拉伸模拟，根据监测的顶面边界受力峰值、峰值点变形计算出模型的弹性模型E和抗拉强度σt，通过公式U* =Dσt / E计算出U*，图 15给出了UECC*和U*的关系图。
 |
| 图 15 UECC*与U*关系 |
| 图选项 |
结果表明，UECC*和U*成正比关系，记为UECC*=αU*，其中α为比例系数。在图 16中，存在2个α值，分别对应Rc=0.5和Rc=0.8的情况，与Ra、γ无关。由以上分析可知，比例系数α与Rc唯一对应，因此，选取参数Ra=0.9、γ=0.001、UECC*=0.001为例，以Rc为自变量研究Rc与α的关系，Rc选取0.45~1.0，间隔0.05，结果如图 16所示。结果表明，当Rc较小时，α随Rc的增长非常快，在Rc=0.6的时候α达到峰值点，峰值约为1.5，之后缓慢下降，而当Rc>0.8之后，α则一直保持在1.3附近。
 |
| 图 16 α与Rc关系 |
| 图选项 |
3 动态破坏算例以二维双孔爆炸问题为例，通过简易荷载模型，分析不同参数对破坏过程的影响。模型为长方形，尺寸为500 mm×250 mm，直径10 mm的2个爆孔位于水平中心线上，中心间距100 mm，如图 17所示。
 |
| 图 17 二维双孔爆炸模型几何尺寸及荷载条件 |
| 图选项 |
荷载持续时间为25 μs，最大荷载为107 kN，在实际计算中爆破荷载按照圆环方式加载，每个颗粒的荷载按照圆环总体荷载除以颗粒数目进行分配。考虑6种不同的情形，即4D-LSM计算、仅引入Ra的ECC4D计算(Ra分别为0.2、0.8)以及引入3个参数(Ra、Rc、γ)的ECC4D计算，其爆炸的破坏过程分别如图 18和19所示。
 |
| 图 18 4D-LSM双孔爆炸破坏过程 |
| 图选项 |
 |
| 图 19 不同偏心系数Ra的爆孔破坏图 |
| 图选项 |
图 18是4D-LSM的模拟结果，不同颜色在DLSM中表示不同状态，直观地表明不同阶段模型发生的变形或者破坏等行为。蓝色代表初始未破坏区域，红色代表拉伸破坏状态，绿色表示拉伸破坏后再闭合状态。破坏首先以裂缝的形式发展，以爆孔为中心呈45°交叉状延伸。接着两孔之间的裂纹转为水平方向发展并相互连接。之后两侧的裂缝分散出一条水平方向的分支，整个破坏形态具有很好的对称性。图 19是只引入Ra的计算结果，整体来看，裂缝起始方向也呈45°交叉状延伸。
Ra较小时的破坏过程差异不大，双孔之间的裂缝继续发展之后贯通，整体具有一定的对称性。图 19a和19b中Ra分别为0.2和0.8，相比Ra比较小的情况，裂缝起始发展的方向开始增多，以爆孔为中心呈放射状，各条裂缝的发展速度大致均等。同时，后期上下两端的中心位置开始出现新的破坏带，整体有一定的对称性。图 19c中引入了3个参数Ra、Rc和γ的ECC4D的计算结果，此时，裂缝的起始方向依然以爆孔为中心呈45°发展。中期开始出现一些其他方向的裂缝，但是2个爆孔之间的裂缝不具有对称性，并且后期也没有在上下两端出现新的破坏带，更能吻合实际爆破过程中裂纹的非对称性扩展分布。
4 大规模前处理方法离散数值方法将研究对象划分成离散单元，利用Newton第二定律和Hooke定律自下而上对物理现象进行重现。很大程度上，模拟结果的可靠度取决于数值计算单元的规模。然而，当模型规模逐渐变大时，数值计算的前处理和后处理的数据量也将异常庞大，因此需要合理地对数据进行压缩、读取、调用。在4D-LSM中，常规处理方式是将模型的全部颗粒文件信息储存起来，向内存中输入海量数据，这就极大地影响了计算效率和读取效率。经过改变存储方式，利用三维有限元网格来表征计算模型，将计算模型离散成若干个简单三维形体并携带相应的材料信息传输到内存中，该有限元模型仅用来表征计算模型的几何信息，不参与数值计算。
如图 20所示，四维弹簧模型的原始文件不再是以往的庞大的颗粒信息文件，而仅仅是一个由有限元软件生成的几何模型文件，计算机首先将数据量非常小的几何文件快速地读入内存，然后通过相应的算法在内存中生成海量颗粒数据，最后进行求解。另一方面，针对百亿颗粒的模型，若保存所有求解结果将造成结果文件过大。可仅对特定的颗粒进行保存，即选择性地保存用户感兴趣的指定颗粒，而不是整个模型数据。例如，如果只对100亿颗粒的立方体的表面颗粒进行保存，其数量仅为2 700万，则后处理工作可以在本地台式工作站上完成。此外，还可以对指定的切片区域进行保存，这样可以进一步减少输出文件所需的存储空间。
 |
| 图 20 计算流程对比 |
| 图选项 |
基于该思路，在本地工作站上已经实现了超过1亿规模的颗粒数值算例试验，本地CPU服务器测试平台为艮泰，CPU的核心数20，线程数40，内存容量256 GB，内存带宽68.3 GB/s。由于ECC4D的基本颗粒是基于规则颗粒生成规则，因此可以直接通过判断规则颗粒的球心与四面体单元的空间关系生成对应的计算模型。如图 21所示的超1亿颗粒的ECC4D计算模型显示了该方法的可行性，通过提取模型表面颗粒进行后处理分析也是可行的，表面颗粒总数约为160万，在个人工作站的处理能力范围内。利用有限元网格建立4D-LSM模型所消耗的时间与颗粒的数量成正比，并且划分的网格越多，计算时间越长，对于上述算例划分的网格，每生成100万的颗粒大约耗时160 s，1亿颗粒模型的生成花费不到5 h，效率比较可观。
 |
| 图 21 大规模前处理方法应用实例(计算颗粒数为104 866 784) |
| 图选项 |
5 结论针对快速构建超大规模非均匀模型困难的现状，本文提出一种非规则偏心四维模型(ECC4D)，以立方体单轴压拉数值实验为例，通过研究弹性和破坏参数的变化规律证明了ECC4D模型的正确性，并通过一个具体的双孔爆破数值算例研究岩石的破坏过程。同时提出了一种解决超大规模模型的生成和处理方法。主要得到了如下结论：
1) ECC4D模型中，弹性参数与输入值参数具有明确的规律性，文中通过引入3个新参数，偏心参数、截断系数和弱化比例来控制模型的非均匀性，并研究了这些参数对模型弹性参数重现的影响规律。
2) ECC4D模型的随机性导致重现的弹性参数是离散的，但是随着模型规模的增大，μs和Es都有收敛的趋势。
3) 破坏参数的研究表明，ECC4D与4D-LSM的破坏参数之间存在一个比例系数α的关系，且α的取值与Rc唯一对应。
4) 超大规模前处理方法在本地工作站上可实现超过1亿离散单元的模型计算，为后续工程应用奠定了基础。
参考文献
| [1] | 赵尚毅, 郑颖人, 邓卫东. 用有限元强度折减法进行节理岩质边坡稳定性分析[J]. 岩石力学与工程学报, 2003, 22(2): 254-260. ZHAO S Y, ZHENG Y R, DENG W D. Stability analysis on jointed rock slope by strength reduction FEM[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(2): 254-260. DOI:10.3321/j.issn:1000-6915.2003.02.018 (in Chinese) |
| [2] | 李世海, 高波, 燕琳. 三峡永久船闸高边坡开挖三维离散元数值模拟[J]. 岩土力学, 2002, 23(3): 272-277. LI S H, GAO B, YAN L. 3D simulation of the excavation of high steep slope of Three-Gorges permanent lock by distinct element method[J]. Rock and Soil Mechanics, 2002, 23(3): 272-277. DOI:10.3969/j.issn.1000-7598.2002.03.004 (in Chinese) |
| [3] | LI X F, ZHANG Q B, LI H B, et al. Grain-based discrete element method (GB-DEM) modelling of multi-scale fracturing in rocks under dynamic loading[J]. Rock Mechanics and Rock Engineering, 2018, 51(12): 3785-3817. DOI:10.1007/s00603-018-1566-2 |
| [4] | 唐春安, 张永彬. 岩体间隔破裂机制及演化规律初探[J]. 岩石力学与工程学报, 2008, 27(7): 1362-1369. TANG C A, ZHANG Y B. Discussion on mechanism and evolution laws of fracture spacing in rock mass[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(7): 1362-1369. DOI:10.3321/j.issn:1000-6915.2008.07.008 (in Chinese) |
| [5] | 王卫华, 李夕兵. 离散元法及其在岩土工程中的应用综述[J]. 岩土工程技术, 2005, 19(4): 177-181. WANG W H, LI X B. A review on fundamentals of distinct element method and its application in geotechnical engineering[J]. Geotechnical Engineering Technique, 2005, 19(4): 177-181. DOI:10.3969/j.issn.1007-2993.2005.04.005 (in Chinese) |
| [6] | 郑颖人, 赵尚毅. 有限元强度折减法在土坡与岩坡中的应用[J]. 岩石力学与工程学报, 2004, 23(19): 3381-3388. ZHENG Y R, ZHAO S Y. Application of strength reduction fem in soil and rock slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(19): 3381-3388. DOI:10.3321/j.issn:1000-6915.2004.19.029 (in Chinese) |
| [7] | 焦玉勇, 张秀丽, 刘泉声, 等. 用非连续变形分析方法模拟岩石裂纹扩展[J]. 岩石力学与工程学报, 2007, 26(4): 682-691. JIAO Y Y, ZHANG X L, LIU Q S, et al. Simulation of rock crack propagation using discontinuous deformation analysis method[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(4): 682-691. DOI:10.3321/j.issn:1000-6915.2007.04.004 (in Chinese) |
| [8] | SHI G H, GOODMAN R E. Two dimensional discontinuous deformation analysis[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1985, 9(6): 541-556. DOI:10.1002/nag.1610090604 |
| [9] | MA G W, AN X M, HE L. The numerical manifold method: A review[J]. International Journal of Computational Methods, 2010, 7(1): 1-32. DOI:10.1142/S0219876210002040 |
| [10] | HRENNIKOFF A. Solution of problems of elasticity by the framework method[J]. Journal of Applied Mechanics, 1941(8): 619-715. |
| [11] | ZHAO G F, FANG J N, ZHAO J. A 3D distinct lattice spring model for elasticity and dynamic failure[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2011, 35(8): 859-885. |
| [12] | ZHAO G F, FANG J, ZHAO J N. A 3D distinct lattice spring model for elasticity and dynamic failure[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2011, 35(8): 859-885. |
| [13] | ZHAO G F, XIA K W. A study of mode-I self-similar dynamic crack propagation using a lattice spring model[J]. Computers and Geotechnics, 2018, 96: 215-225. DOI:10.1016/j.compgeo.2017.11.001 |
| [14] | WANG Y B, YANG R S. Study of the dynamic fracture characteristics of coal with a bedding structure based on the NSCB impact test[J]. Engineering Fracture Mechanics, 2017, 184: 319-338. |
| [15] | ZHAO Y X, ZHAO G F, JIANG Y D, et al. Effects of bedding on the dynamic indirect tensile strength of coal: Laboratory experiments and numerical simulation[J]. International Journal of Coal Geology, 2014, 132: 81-93. |
| [16] | ZHAO G F. Developing a four-dimensional lattice spring model for mechanical responses of solids[J]. Computation Methods in Applied Mechanics and Engineering, 2017, 315: 881-905. |
| [17] | 张奔, 赵高峰. 基于四维离散数值方法的岩石圆环试样动态破坏及耗能规律[J]. 土木与环境工程学报, 2019, 41(2): 20-28. ZHANG B, ZHAO G F. Dynamic failure and energy dissipation of rock ring specimen based on 4D lattice spring model[J]. Journal of Civil and Environmental Engineering, 2019, 41(2): 20-28. (in Chinese) |
| [18] | HU X D, ZHAO G F, DENG X F, et al. Application of the four-dimensional lattice spring model for blasting wave propagation around the underground rock cavern[J]. Tunnelling and Underground Space Technology, 2018, 82: 135-147. |
| [19] | 郑文刚, 刘凯欣. 离散元法工程计算软件的前后处理系统[J]. 计算机工程与科学, 2000, 22(6): 14-15, 23. ZHENG W G, LIU K X. A preprocessing and postprocessing system for an engineering computation software of discrete element methods[J]. Computer Engineering and Science, 2000, 22(6): 14-15, 23. (in Chinese) |
| [20] | 赵高峰, 陈华. 基于阿里云的四维弹簧模型并行运算性能[J]. 土木与环境工程学报, 2019, 41(3): 1-10. ZHAO G F, CHEN H. Performance of the parallel four-dimensional lattice spring model using Alibaba cloud[J]. Journal of Civil and Environmental Engineering, 2019, 41(3): 1-10. (in Chinese) |
| [21] | LIU C, POLLARD D D, GU K, et al. Mechanism of formation of wiggly compaction bands in porous sandstone: 2 numerical simulation using discrete element method[J]. Journal of Geophysical Research: Solid Earth, 2015, 120(12): 8153-8168. |
