颗粒材料如堆石体、道砟、砂等被广泛应用于大型工程建设中，这些天然材料的工程力学特性不仅受颗粒强度、颗粒间的摩擦系数、级配、孔隙率等物理属性的影响^[1]，其不规则的颗粒形状也是影响其宏细观力学行为的关键因素之一.目前，通过室内试验研究颗粒形状对颗粒材料力学特性的影响工作量大，且试验条件控制难度大，因此，数值模拟成为研究这一影响的有效方法.

在有限元-离散元耦合分析方法(Combined Finite-Discrete Element Method, FDEM)中，可以采用随机多边形(或多面体)颗粒来模拟颗粒形状^[2-6]，然而，当颗粒形状拓扑结构比较复杂(如凹多边形)时，FDEM接触检索方法将无法准确识别接触颗粒.在离散元(Discrete Element Method, DEM)模拟中，颗粒的形状多为二维圆盘或三维圆球^[7]，无法考虑颗粒之间的嵌入咬合作用，数值模拟结果与实际颗粒的力学性质存在一定差异.为此，有学者将形状简单的圆盘或者圆球通过一定的方式组合在一起形成“颗粒簇”或“超级颗粒”(clump单元)来反映颗粒的不规则形状^[8-10]，clump单元内部的圆盘或圆球可以任意重叠，内部的颗粒之间不产生接触力，整个clump单元作为一个完整的刚性体而存在，其边缘的颗粒与其外部颗粒产生相互作用.此外，在DEM方法中也有以椭圆或者椭球为基本颗粒的，通过改变椭圆的长短轴之比来模拟不同颗粒形状^[11-12].上述DEM研究只能模拟一些形状简单且规则的颗粒，对于岩土工程中的复杂形状颗粒往往需要借助图像获取技术才能实现真实颗粒形状的再现.Wang等^[13]和金磊等^[14]采用CT扫描技术获得颗粒表面信息后，在离散元模拟中用“超级颗粒”重构了真实颗粒形状，但CT扫描方法能提供的颗粒形状很有限，无法反映真实岩土颗粒形状的随机性和多样性.

目前，通过室内试验研究颗粒形状对颗粒集合体宏观力学影响的研究较少；采用抗转动模型^[15-16]等效考虑颗粒形状，增加了细观参数且无法充分揭示颗粒形状对细观力学响应的影响.为此，本文提出了复杂形状颗粒DEM模拟方法，首先在有限元软件ANSYS中对生成的随机多边形进行网格划分，随后在PFC^2D中用圆盘替换网格，并将多个圆盘交接在一起模拟非规则形状的颗粒(clump单元)，实现了包括凹多边形在内的不规则随机多边形颗粒模拟；通过双轴压缩试验，获得了不同形状颗粒试样的宏观力学响应，揭示了颗粒形状对宏观力学响应的细观影响机理.

1 随机多边形生成算法 1.1 多边形颗粒的生成

复杂多边形颗粒由边数、极角、圆半径表征确定，如图 1所示.

具体生成方法包括如下6个步骤：

步骤1，在指定区域内按照给定级配曲线生成随机圆形颗粒.此时，颗粒圆心位置(x₀，y₀)和半径r已知.

步骤2，指定最大边数和最小边数，生成多边形随机边数：

式中：n为生成随机多边形的边数；n_min为最小边数；int[]为取整函数；n_max为最大边数；rand为0~1之间的随机数.

步骤3，生成第k条边对应的随机角度:

式中：θ_k为第k条边的对角；δ为小于1的常数.上述n个随机角θ_k之和一般不闭合，需要进行修正：

式中：为修正后的第k条边对角.

步骤4，获取多边形顶点坐标：

式中：λ为小于1的正常数；grand为服从Gaussian正态分布-1~1之间的随机数；θ′_k为修正后的极角；β为颗粒随机方位角；当λ=0时生成的随机多边形为凸多边形，λ≠0时，可以生成凹多边形.

步骤5，根据以上复杂多边形的边数、极角、圆半径等几何信息，在有限元软件ANSYS中生成多边形面，并对各个面进行网格划分.

步骤6，通过在有限元软件ANSYS中获取的网格形心位置及网格面积等信息，在离散元软件PFC^2D中用等面积同形心圆盘替换各个网格，并将同一多边形内的圆盘采用clump命令组构成为“超级颗粒”，以此模拟复杂随机多边形.

图 2、图 3分别为λ=0和λ=0.2时生成的随机多边形.当λ=0时，所生成的随机多边形为圆内接凸多边形，随着边数的增加，随机多边形的形状越接近圆.当λ≠0时，本文提出的随机多边形生成方法可以生成凸多边形和凹多边形，当边数较少时，所生成的多边形为凸多边形但与圆无内接关系，随着边数的增加，颗粒形状越来越复杂.

图 4为工程中实际岩石颗粒，对比图 2~4，本文提出的复杂形状颗粒模拟方法可以通过参数控制使得生成的颗粒形状更加真实.

1.2 颗粒表面粗糙度敏感性研究

对生成的多边形颗粒进行有限元网格划分，以5边形和10边形颗粒为例，图 5和图 6所示为DEM模拟的颗粒形状.

由图 5和6可知，随着小球数目的增加(网格加密)，颗粒形状的边界被模拟得越真实.当离散颗粒所采用的网格尺寸越小，所得颗粒边界的凹凸程度(粗糙度)越小但所需的小颗粒数目越大.为了寻找到一个合理计算规模又能控制颗粒粗糙度对宏观力学响应的影响，本文对随机凸多边形采用不同尺寸网格进行离散，网格尺寸以多边形每条边划分的份数为控制指标，最后生成3组不同颗粒粗糙度的同一随机多边形离散元计算模型(颗粒粒径线性分布)，计算参数见表 1.

根据表 1中参数生成不同颗粒粗糙度的模型，计算结果见图 7.由图可知，颗粒集合体的体积变形与粗糙度关系不明显且差异较小，粗糙度对颗粒集合体的峰值偏应力和残余偏应力的影响较小.因此，可以认为采用划分密度大于每条边3个单元的网格离散后所得到的随机多边形离散元计算结果对颗粒表面粗糙度不敏感.

1.3 颗粒形状的描述

Barrett^[17]提出将颗粒形状量化指标分为颗粒形态、圆度以及表面粗糙度.颗粒表面的粗糙程度通常隐含在颗粒的摩擦系数中，因此只需引入两个颗粒形状指标，即圆度和凹凸度.参考孔亮等^[18]对颗粒形状的表征方法，对于二维情况，圆度和凹凸度定义为

式中：S_2D和F分别为颗粒的圆度和凹凸度；A_f和A′_f分别为颗粒面积和颗粒最大内接标准椭圆面积；A_s为与多边形颗粒等周长圆的面积.为了更便于描述颗粒形状，用颗粒的圆度和凹凸度的加权平均获得表征颗粒外轮廓特征的形状系数，即形状系数K=α₁S_2D+α₂F，其中，α₁+α₂=1，α₁和α₂的取值与颗粒的实际形状密切相关，且与颗粒组的受力变形特性也存在联系.在变形过程中，假设颗粒的自转分量比颗粒的平移分量要大，颗粒的圆度影响更大；若平移的分量大，则凹凸度的影响更大^[18].本文取形状系数为圆形度和凹凸度的平均值，即K=(S_2D+F_2D)/2.

颗粒圆度和凹凸度的理论最大值为1.0(纯圆颗粒)，图 8为不同规则多边形的形状指标.由图 8可知：随着多边形边数的增加，颗粒的圆度、凹凸度及形状系数K均趋近于1.0.

2 颗粒圆度对试样力学响应研究 2.1 试样的制备

为了分析颗粒圆度对颗粒力学特性的影响，本文建立了6个不同规则多边形形状的颗粒集合体(颗粒边数从5变化到10)，通过双轴压缩试验研究颗粒形状对宏观力学响应的影响.试样尺寸为120 mm×60 mm，采用的计算参数见表 1.首先生成了875个圆形颗粒，采用前述方法，取λ=0，控制n的大小，得到表征各个复杂多边形的边数和顶点坐标信息.基于以上信息在有限元软件ANSYS中生成相应的多边形，并对其进行网格控制划分，最后将网格数据导入PFC^2D中，用等形心等面积替换法生成不同形状的DEM颗粒.为了使不同形状的颗粒集合体力学响应具有可比性，不同颗粒集合体均在1.0 MPa围压和颗粒摩擦系数μ=0下等向固结，再将摩擦系数设置为0.5后进行加载试验，即保证不同颗粒集合体在加载前均处于最密实状态^[19].图 9所示为加载前不同形状颗粒集合体的接触分布图.

图 9中蓝色线条表示颗粒之间接触数大于1(即边-边接触)，红色线条表示颗粒之间只有1个接触关系(即点-点接触或者点-边接触).可以看出，在最密实状态下颗粒间的孔隙填充与颗粒边数(顶点数)有关，当颗粒边数越少时，规则颗粒边-边接触所占的比例越小，其内部的“拱效应”越不明显，形成的“无效颗粒”越少，颗粒填充越密实，孔隙率n越小.

2.2 宏观结果分析

为了研究颗粒形状对试样应力变形影响，采用峰值摩擦角φ和剪胀角ψ进行定量描述，按照下式计算：

式中：σ₁为大主应力；σ₃为施加的围压；ε_v为体积应变；ε_a为轴向应变；ε_50%为对应于50%峰值偏应力点的轴向应变，对于没有明显峰值的应力应变曲线，峰值偏应力取15%轴向应变对应的偏应力.

图 10为不同试样在双轴压缩下的偏应力及体积变形曲线.从图中可知，不同形状的颗粒集合体在加载初期的变形模量区别不大，体积变形均表现出较小的加载初期剪缩和加载后期剪胀，但强度并不随颗粒边数的增加而单调递减，当规则颗粒边数为5~10时，颗粒集合体的峰值强度和体积变形特性差异较小.

图 11为不同形状试样峰值摩擦角与形状系数K关系曲线，由图可知试样峰值摩擦角与颗粒形状系数存在很好的线性关系.图 12为试样颗粒平均旋转速度演化图，不同形状颗粒的平均旋转速度演化规律与应力应变曲线规律一致，随着颗粒形状边数的增加，平均旋转速度总体上呈增大趋势，但颗粒边数为5~10时，相互区别不大.

2.3 细观机理分析

颗粒旋转主要有2种模式：绕形心自转和绕顶点转动，如图 13.不同边数的颗粒外接圆直径相同，在最密实状态下不同形状颗粒绕形心自转同一个角度θ时，颗粒顶点转动的轨迹与外接圆重合，颗粒边数越少，颗粒每一个顶点与圆心连线转过的面积总和与颗粒面积之差越大.当颗粒转动一个圆心角时，颗粒顶点与圆心连线转过的面积之和为外接圆的面积，颗粒边数越多，颗粒面积越接近其外接圆面积，即颗粒边数越多，颗粒转动所需的额外空间越少.以5边形和6边形为例，5边形颗粒转动一个圆心角时，试样体积增大10A₁，6边形试样体积增大12A₂，显然有10A₁＞12A₂，因此颗粒绕形心自转时，颗粒边数越少局部剪胀越大，当颗粒边数相差越大时，剪胀变形差异越大.由于点-点或点-边接触比较不稳定，因此这种接触容易绕颗粒接触点转动，且在加载过程中，颗粒顶点数越多，发生点-点或点-边接触的概率越大.当颗粒绕顶点转动时，偶数边颗粒的旋转半径R₂等于其外接圆的直径，奇数边颗粒旋转半径R₁ < R₂，因此偶数边颗粒绕顶点转动时，试样剪胀变形更大.

图 14为加载过程中不同试样边-边接触比例关系演化图，从图可知在最密实状态下颗粒接触以边-边接触为主，随着加载进行，边-边接触比例逐渐降低.图 10中的峰值强度和体积变形并没有出现与颗粒形状单调变化的关系，这可以用颗粒旋转加以解释：加载过程中颗粒绕自身旋转，导致试样局部体积剪胀，颗粒边数越少剪胀越严重，颗粒“摆脱”周围约束所需能量越大，因此试样峰值强度越大；当颗粒边数相差很小时，颗粒顶点数的增加使得点-点和点-边接触的概率增加，颗粒发生绕顶点旋转的可能性增大，因此偶数边颗粒与奇数边颗粒的体积变形差异较小.同时，与绕形心旋转相比，颗粒绕顶点旋转时，颗粒转动惯量增大，所需外力矩增大，因此偶数边颗粒试样的峰值强度未明显降低.

下式是由Oda^[20]提出的定量描述颗粒接触法向的组构张量，组构张量主值φ₁₁和φ₂₂能够很好地表征颗粒接触形成和失效程度，在二维情况下有φ₁₁+φ₂₂=1.该式为

式中：n为接触面法向的单位矢量；Θ描述了全局坐标系中单位矢量n法向；E(Θ)是接触法向的概率密度函数.

图 15为加载结束时颗粒接触分布图，由图可知，颗粒边数越少，拱效应越明显，拱内的空隙越大.图 16为试样组构张量主值φ₁₁演化图(下角标1表示为围压方向，下角标2表示为加载方向)，由图可知在加载过程中，由于颗粒滑动和颗粒转动(即颗粒位置调整)，加载方向上不断有新接触形成而围压方向上的接触不断失效，且颗粒边数越少，试样中颗粒接触失效与重构越激烈.图 15和图 16的结果印证了上述对图 10的解释，即颗粒边数越少颗粒位置调整越强烈，其形成的力链稳定结构(拱结构力链)越多，试样整体强度最高且剪胀最明显.

3 颗粒凹凸度对试样力学响应研究

在上节875个圆形颗粒中，随机生成5~10边的随机多边形，分别取λ=0，0.1，0.2，0.3生成4组不同凹凸度的试样，并在1.0MPa围压和颗粒摩擦系数μ=0下等向固结，再将摩擦系数设置为0.5后进行加载试验.

图 17为不同试样在双轴压缩下的偏应力及体积变形曲线.从图中可知，当试样颗粒全部为凸颗粒(λ=0)时，试样的峰值强度最低，且变形模量也最小.当λ≠0时，试样的变形模量比λ=0时大，但不同λ的试样初始变形模量差异很小；随着λ值的增大，试样的峰值强度先增大后减小，说明随着颗粒凹凸度变化幅度的增加，试样在最密实状态下的咬合作用明显增强，但这种咬合作用在λ=0.2时最强.当λ从0增大到0.2时，试样的剪胀变形增大，而后开始变小，这与峰值偏应力结果一致.说明在最密实状态下，当颗粒形状更加随机多样时，颗粒间的排列更加密实，咬合作用更强，颗粒的旋转需要更多的外部能量和旋转空间；同时，颗粒形状随机多样性与颗粒间的咬合作用并无单调关系，当λ=0.2时咬合作用最强，峰值强度最高，变形模量最大.

4 结论

本文提出了复杂形状颗粒DEM模拟方法，通过双轴压缩试验研究了颗粒形状对强度变形的影响，并从细观层次上解释其影响机制，主要结论如下：

1) 提出了复杂形状颗粒DEM模拟方法，在PFC^2D中实现了包括凹多边形在内的复杂形状颗粒模拟，通过参数控制使得生成的颗粒形状更加真实.

2) 在最密实状态下，颗粒间的孔隙填充与颗粒边数有关，颗粒边数越少，孔隙率越小.这主要是由于颗粒的边数越少，周长越短，力链中的稳定接触(边-边接触)所占比例越少，试样内部的拱效应越不明显，形成的“无效颗粒”越少，颗粒填充越密实.

3) 在最密实状态下，不同形状颗粒集合体加载初期的变形模量很接近，偏应力和体积变形规律相似，峰值摩擦角与形状系数K呈线性关系，但规则颗粒边数为5~10时，颗粒集合体的峰值强度和体积变形特性差异较小.

4) 在双轴压缩过程中，颗粒边数越少的试样，颗粒位置调整(颗粒滑动和颗粒转动)越剧烈.颗粒形状对颗粒集合体宏观力学响应的影响可以从细观层次上解释为：颗粒边数越少，颗粒绕形心转动所需外界功越大，所导致的体积剪胀也越明显；颗粒绕接触点旋转时，偶数边颗粒所需外界功和导致的体积剪胀比奇数边颗粒试样大，因此随着颗粒边数的降低，试样峰值强度和剪胀变形越大.边数为5~10的试样宏观响应差异很小这一现象可通过相应的边-边接触比例关系演化结果以及接触法相组构张量演化结果得到解释.

5) 在最密实状态下，当试样颗粒全部为凸颗粒时，试样的峰值强度和变形模量最小.随着颗粒凹凸度变化幅度的增加，颗粒间的排列更加密实，咬合作用更强，颗粒的旋转需要更多的外部能量和旋转空间，但颗粒形状随机多样性与颗粒间的咬合作用并无单调关系，当λ=0.2时咬合作用最强.

6) 本文所提出的这一复杂形状颗粒模拟方法中，为使颗粒形状的边界被模拟得更真实，离散颗粒所采用的网格需划分得更细，即颗粒数目更大，计算工作量会更大.但是以目前的计算机硬件水平，在二维条件下的计算速度是可以接受的.如何优化此方法，减少颗粒数目，将是下一步的研究工作.