岩石的物理力学性质研究对于工程实践具有重要意义，面对工程建设规模日益扩大、地质条件愈加复杂、载荷和边界条件复杂多变的现实情况下，常规的力学分析方法难以获得精确解，计算力学的发展极大地解决了该难题，不仅加快了工程问题的求解速度，还具有成本低、可重复性强、数据全面等优点，工程中已经有较多的应用^[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。

Type A的质心偏移方式如下：在球形颗粒中，质心在以形心为中心、边长为L的立方体中随机移动，该立方体形心与颗粒形心一致，如图 2所示。设颗粒i的形心坐标分别为x_i^C、y_i^C和z_i^C，颗粒半径为R，立方体的边长为L，定义偏心系数R_a=(L/2)/R，则新的质心坐标x_i^ECC、y_i^ECC和z_i^ECC可表示为

$ \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)

其中，Rand()为随机函数，可随机生成某区间内的一个数，上述3个Rand(-1, 1)均相互独立。

Type B的质心偏移方式如下：在球形颗粒中，质心只能在以形心为中心，半径为r的球面上随机移动，该球面所在球体的形心与颗粒形心一致，如图 1和图 2所示。设颗粒i的形心坐标分别为x_i^C、y_i^C和z_i^C，颗粒半径为R，限制质心移动范围的球体半径为r，定义偏心系数R_a=r/R，则新的质心坐标x_i^ECC、y_i^ECC和z_i^ECC可表示为

$ \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)

其中：θ和φ是随机产生的球面坐标分量，分别表示向量与z轴和x轴的夹角，具体产生的方式为

$ \left\{\begin{array}{l} \theta=\arccos (\operatorname{Rand}(0, 1)) \\ \varphi=2 \pi \cdot \text { Rand }(-1, 1) \end{array}\right. $

(3)

1.2 几何模型构建

模型的几何结构主要由弹簧和颗粒构成，颗粒具有不变性，因此弹簧的长度和方向决定了三维模型的最终构成。在4D-LSM中弹簧的方向和长度具有固定的排列规律，在ECC4D中，由于质心随机偏移，造成了弹簧长度和弹簧方向的随机性，没有固定的比例。关于4D-LSM中弹簧键的方向性，模型晶格的弹簧键固定，同样的晶格按照一定的规则组合成整体模型。如果将模型中所有弹簧键的一端移至颗粒的形心处，并将这一端作为球面坐标系下向量的起点。那么，所有这些弹簧向量的方向必然是特定的。而ECC4D由于质心在一定范围内随机移动过，因此其弹簧向量的方向就呈现出了更多的可能性，排列分布更加具有随机性。对于弹簧键的长度，与弹簧向量的方向同理。4D-LSM中弹簧键的长度也是有固定规律的，一般有3种长度的弹簧键，即正弹簧、斜弹簧和对角弹簧，并且3种弹簧键的数量比例是一定的。因此，整个模型弹簧键的平均长度也是固定的，而ECC4D中弹簧键则长短不一，并且长、短弹簧键的数量也没有固定的比例，弹簧的平均键长是不固定的。图 3给出了4D-LSM和ECC4D的三维模型，可以清晰地观察2种模型的区别。

2 模型验证分析 2.1 偏心系数

在提出ECC4D模型时引入了偏心系数R_a，为了研究R_a对几何构型的影响，分别对不同R_a取值下弹簧键的方向分布和长度分布进行统计。考虑R_a=0.1~0.7，间隔为0.2，结果如图 4和5所示。对于弹簧方向分布规律，Type A和Type B这2种模型的变化规律是一致的。随着R_a的增加，弹簧键在球面上以某些点为中心逐渐地朝整个球面扩散；对于弹簧长度分布，2种模型的规律也趋于一致，从间断的双峰分布逐渐过渡到连续的双峰分布最终到单峰分布。图 6表示R_a为0.1、0.5、0.9这3种偏心系数下2种模型的弹簧连接键的结构图，可以看到模型的棱角逐渐减小而趋于“圆润”。通过引入R_a，弹簧键方向和长度从原来的“确定量”逐步转变成了具有一定统计规律的“随机量”，即ECC4D在一定程度上表征了岩石材料的几何非均匀特征。后文若未加说明，则默认的研究对象是Type A方式生成的ECC4D。

2.2 Poisson比

Poisson比作为材料重要的力学参数，是参数研究的重点对象。在4D-LSM中，通过试验获取的Poisson比μ和数值计算得到的Poisson比μ_s是相等的，而ECC4D在构建模型的过程中随机偏移了质心，其μ_s和μ的关系是否吻合，还需要进一步研究验证。

相对于4D-LSM，ECC4D引入了偏心系数R_a=r/R，其中r是限制质心移动生成范围的球体的半径(或者是立方体边长的一半)，R表示颗粒半径。因此，R_a的取值范围为(0, 1)，R_a=0表示质心没有偏移，R_a=1表示质心位于颗粒表面的极限位置。在物理世界，2个物体的相互作用大多数都是随着距离的增加而减小，当距离增加到一定程度，这种相互作用就可以忽略不计了。在ECC4D中同样考虑这种规律，通过引入另一个参数R_c，取值范围同为(0, 1)，称之为截断系数。其具体意义如下：记ECC4D中弹簧键的最大长度为L_max，那么，当长度超过R_c×L_max时，弹簧键将被截断。当R_c=1时，R_c×L_max=L_max的弹簧键长已经是最大弹簧键长了，因此没有弹簧键被截断；当R_c=0时，R_c×L_max=0，此时所有的弹簧键都将被截断。

建立一个立方体单轴压缩模型(见图 7)进行Poisson比参数分析，通过测点监测横向和竖向上的变形并计算数值模拟所得μ_s，选取3种截断系数R_c进行模拟，输入的μ为0.2，研究R_a与μ_s的关系。如图 8，当截断系数R_c=1.0时，即弹簧键均未被截断，可以看作未引入截断系数。随着R_a的增加，重现μ_s整体呈现单调下降趋势。当引入截断系数后，曲线的变化趋势不再呈现单调关系，存在峰值点，并且当R_c=0.6时，μ_s已经突破了μ。其数值模型也发生了相应的变化，连接键从有序逐渐过渡到随机。因此，可以推测ECC4D模型重现的Poisson比可突破0.25的限制。

为了进一步探究μ_s峰值变化情况，考虑不同R_a条件下R_c和μ_s的变化关系。如图 8所示，一方面曲线整体均为先上升后下降的规律，部分曲线已经突破了Poisson比的限制，超过了0.25，最大的重现Poisson比已经达到了0.36。另一方面，不同的μ对应的曲线变化规律是一致的，峰值点的位置与R_a是一一对应的。但是，直接截断弹簧键的处理方式也存在一些问题，图 8中R_c小于0.35的部分没有曲线变化数据，这是由于在偏心系数R_a为0.7和0.8的条件下，截断系数过小会产生计算错误，因此曲线的变化形式存在一定的缺陷。初步分析原因可能由于ECC4D是基于4D-LSM演变而来，在4D-LSM中，颗粒所连接的弹簧数量有限并且长度固定；ECC4D首先通过随机生成质心的方式使得键长在某一范围内成为了随机量，而引入R_c使得“过长”的弹簧被截断，但这也导致某些颗粒的所有弹簧键有可能都因为“过长”而被截断，致使该颗粒处于“无连接状态”，进而使模型本身处于不稳定状态并造成计算错误。

为了解决上述问题，将弹簧键直接截断的简单处理方式变为“弱化”。当过长的弹簧需要截断时，将其弹簧刚度按比例大幅度降低，弱化比例记为γ=k_c/k，其中k_c为弱化后的弹簧刚度，k为原弹簧刚度。图 9表明输入Poisson比的不同对曲线变化规律没有影响，下文输入Poisson比均采取μ=0.25进行研究。引入弱化系数后，重新探究R_c和μ_s的关系，结果如图 10和11所示。当R_a=0.0~0.5时，重现Poisson比与截断系数的变化规律表现为先降低后升高，不再是先升高后降低的变化规律。并且重现Poisson比对截断系数的敏感性较大，致使R_c以0.1为间距变化时无法有效地找到其最小值点。当R_a=0.6~0.9时，重现Poisson比和截断系数的变化规律表现为先升高后降低，曲线变化形式较为平缓，对于不同的R_a，μ_s峰值位置都是R_c=0.5，并且峰值随着R_a的增大有一定的增幅。

μ_s的峰值也是一个重要的研究参数，记为μ_s-max。图 11中可以看出，曲线呈现一定的规律性，即在R_c=0.5处出现峰值且R_a=0.9时峰值最大。因此，相关参数分别取为R_a=0.9、R_c=0.5、μ=0.25，研究弱化系数和重现Poisson比的峰值之间的变化规律，同时也对重现的弹性模量E_s进行分析，并且考虑了Type A和Type B两种生成方式的影响，其结果如图 12所示。E_s对生成方式的影响以E_s/E表示，E_s是重现的弹性模量，E是输入的弹性模量。由图 12可知，随着弱化系数的增加，重现Poisson比呈现下降的趋势，重现弹性模量呈现上升的趋势。其原因是γ的值越大，代表弹簧键弱化的程度越小，即弹簧刚度越大，进而使Poisson比减小而弹性模量增大；Type A和Type B的2条μ_s-max曲线几乎平行，Type A的μ_s-max稍大，两者的最大值都达到了0.3左右；Type A和Type B的2条E_s曲线在γ较小时差别较大，随着γ的增大，两者的差距逐渐减小并且有重合的趋势，同时，由于部分弹簧键被弱化，E_s/E的比值始终小于1。

由于ECC4D的模型生成是随机的。因此，其弹性参数必然存在一定的波动性，但是ECC4D的目的主要是为了快速的构建非均匀模型并用于大规模计算。因此，需要对ECC4D弹性参数是否能够随着模型规模的扩大而收敛进行研究。模拟依然采用立方体单轴压缩模型，但是改变模型的物理尺寸，边长10 ~ 80 mm，以10 mm递增。其他参数分别为：R_a=0.9，R_c=0.5，γ=0.001，μ=0.25，参数相同的模型随机生成10次并计算出μ_s和E_s，结果分别如图 13和14所示。图 13表明，随着模型规模的增大，μ_s-max具有很好的收敛性，收敛值大约为0.3；图 14中，E_s同样也有收敛的趋势，但是当模型规模较小时，先有一段上升的趋势。

2.3 破坏参数

当采取质心随机偏移后，ECC4D模型获取的破坏参数与4D-LSM是否存在差异，需要进一步讨论。4D-LSM中弹簧键的变形定义为U，当U达到最大的极限值U^*时弹簧键即断开，U^*是4D-LSM的唯一破坏参数，其值可以用公式U^*=Dσ_t/E进行估算，其中：D为颗粒直径，σ_t为抗剪强度，E为弹性模量。同样采用单轴拉伸数值实验来研究ECC4D中破坏参数U_ECC^*的变化规律，单轴拉伸模型仍然采用立方体，边长20 mm，与图 7的单轴压缩模型相同，但是荷载方向改为竖直向上，颗粒直径为1 mm，R_a=0.9，R_c=0.5，弱化比例γ=0.001。首先，考虑不同的R_a、R_c和γ，输入的破坏参数为U_ECC^*=0.001~0.001 6，间隔0.000 1，分别进行单轴拉伸模拟，根据监测的顶面边界受力峰值、峰值点变形计算出模型的弹性模型E和抗拉强度σ_t，通过公式U^* =Dσ_t / E计算出U^*，图 15给出了U_ECC^*和U^*的关系图。

结果表明，U_ECC^*和U^*成正比关系，记为U_ECC^*=αU^*，其中α为比例系数。在图 16中，存在2个α值，分别对应R_c=0.5和R_c=0.8的情况，与R_a、γ无关。由以上分析可知，比例系数α与R_c唯一对应，因此，选取参数R_a=0.9、γ=0.001、U_ECC^*=0.001为例，以R_c为自变量研究R_c与α的关系，R_c选取0.45~1.0，间隔0.05，结果如图 16所示。结果表明，当R_c较小时，α随R_c的增长非常快，在R_c=0.6的时候α达到峰值点，峰值约为1.5，之后缓慢下降，而当R_c>0.8之后，α则一直保持在1.3附近。

3 动态破坏算例

以二维双孔爆炸问题为例，通过简易荷载模型，分析不同参数对破坏过程的影响。模型为长方形，尺寸为500 mm×250 mm，直径10 mm的2个爆孔位于水平中心线上，中心间距100 mm，如图 17所示。

荷载持续时间为25 μs，最大荷载为10⁷ kN，在实际计算中爆破荷载按照圆环方式加载，每个颗粒的荷载按照圆环总体荷载除以颗粒数目进行分配。考虑6种不同的情形，即4D-LSM计算、仅引入R_a的ECC4D计算(R_a分别为0.2、0.8)以及引入3个参数(R_a、R_c、γ)的ECC4D计算，其爆炸的破坏过程分别如图 18和19所示。

图 18是4D-LSM的模拟结果，不同颜色在DLSM中表示不同状态，直观地表明不同阶段模型发生的变形或者破坏等行为。蓝色代表初始未破坏区域，红色代表拉伸破坏状态，绿色表示拉伸破坏后再闭合状态。破坏首先以裂缝的形式发展，以爆孔为中心呈45°交叉状延伸。接着两孔之间的裂纹转为水平方向发展并相互连接。之后两侧的裂缝分散出一条水平方向的分支，整个破坏形态具有很好的对称性。图 19是只引入R_a的计算结果，整体来看，裂缝起始方向也呈45°交叉状延伸。

R_a较小时的破坏过程差异不大，双孔之间的裂缝继续发展之后贯通，整体具有一定的对称性。图 19a和19b中R_a分别为0.2和0.8，相比R_a比较小的情况，裂缝起始发展的方向开始增多，以爆孔为中心呈放射状，各条裂缝的发展速度大致均等。同时，后期上下两端的中心位置开始出现新的破坏带，整体有一定的对称性。图 19c中引入了3个参数R_a、R_c和γ的ECC4D的计算结果，此时，裂缝的起始方向依然以爆孔为中心呈45°发展。中期开始出现一些其他方向的裂缝，但是2个爆孔之间的裂缝不具有对称性，并且后期也没有在上下两端出现新的破坏带，更能吻合实际爆破过程中裂纹的非对称性扩展分布。

4 大规模前处理方法

离散数值方法将研究对象划分成离散单元，利用Newton第二定律和Hooke定律自下而上对物理现象进行重现。很大程度上，模拟结果的可靠度取决于数值计算单元的规模。然而，当模型规模逐渐变大时，数值计算的前处理和后处理的数据量也将异常庞大，因此需要合理地对数据进行压缩、读取、调用。在4D-LSM中，常规处理方式是将模型的全部颗粒文件信息储存起来，向内存中输入海量数据，这就极大地影响了计算效率和读取效率。经过改变存储方式，利用三维有限元网格来表征计算模型，将计算模型离散成若干个简单三维形体并携带相应的材料信息传输到内存中，该有限元模型仅用来表征计算模型的几何信息，不参与数值计算。

如图 20所示，四维弹簧模型的原始文件不再是以往的庞大的颗粒信息文件，而仅仅是一个由有限元软件生成的几何模型文件，计算机首先将数据量非常小的几何文件快速地读入内存，然后通过相应的算法在内存中生成海量颗粒数据，最后进行求解。另一方面，针对百亿颗粒的模型，若保存所有求解结果将造成结果文件过大。可仅对特定的颗粒进行保存，即选择性地保存用户感兴趣的指定颗粒，而不是整个模型数据。例如，如果只对100亿颗粒的立方体的表面颗粒进行保存，其数量仅为2 700万，则后处理工作可以在本地台式工作站上完成。此外，还可以对指定的切片区域进行保存，这样可以进一步减少输出文件所需的存储空间。

基于该思路，在本地工作站上已经实现了超过1亿规模的颗粒数值算例试验，本地CPU服务器测试平台为艮泰，CPU的核心数20，线程数40，内存容量256 GB，内存带宽68.3 GB/s。由于ECC4D的基本颗粒是基于规则颗粒生成规则，因此可以直接通过判断规则颗粒的球心与四面体单元的空间关系生成对应的计算模型。如图 21所示的超1亿颗粒的ECC4D计算模型显示了该方法的可行性，通过提取模型表面颗粒进行后处理分析也是可行的，表面颗粒总数约为160万，在个人工作站的处理能力范围内。利用有限元网格建立4D-LSM模型所消耗的时间与颗粒的数量成正比，并且划分的网格越多，计算时间越长，对于上述算例划分的网格，每生成100万的颗粒大约耗时160 s，1亿颗粒模型的生成花费不到5 h，效率比较可观。

5 结论

针对快速构建超大规模非均匀模型困难的现状，本文提出一种非规则偏心四维模型(ECC4D)，以立方体单轴压拉数值实验为例，通过研究弹性和破坏参数的变化规律证明了ECC4D模型的正确性，并通过一个具体的双孔爆破数值算例研究岩石的破坏过程。同时提出了一种解决超大规模模型的生成和处理方法。主要得到了如下结论：

1) ECC4D模型中，弹性参数与输入值参数具有明确的规律性，文中通过引入3个新参数，偏心参数、截断系数和弱化比例来控制模型的非均匀性，并研究了这些参数对模型弹性参数重现的影响规律。

2) ECC4D模型的随机性导致重现的弹性参数是离散的，但是随着模型规模的增大，μ_s和E_s都有收敛的趋势。

3) 破坏参数的研究表明，ECC4D与4D-LSM的破坏参数之间存在一个比例系数α的关系，且α的取值与R_c唯一对应。

4) 超大规模前处理方法在本地工作站上可实现超过1亿离散单元的模型计算，为后续工程应用奠定了基础。