颗粒棱度指标的改进及其对剪切特性的影响

袁 斌，霍宇翔，巨能攀，宋浩燃

(地质灾害防治与地质环境保护国家重点实验室(成都理工大学)，四川 成都 610059)

颗粒材料作为建筑材料广泛应用于交通土建工程中，而颗粒形状作为颗粒材料一个及其重要的属性参数，其对颗粒宏细观力学性质的影响一直是颗粒材料研究的一个热点问题。为了更好地了解颗粒形状对颗粒材料力学性质的影响，学者们进行了大量的室内试验和数值模拟试验研究。常晓林等[1]利用三维变形体离散元单元法研究了颗粒的长短径比在宏观和细观层面对堆石体力学性能的影响；王蕴嘉等[2]利用PFC3d研究了堆石料颗粒形状对堆积密度及强度的影响；张成功等[3]利用颗粒料软件研究了颗粒形状对颗粒材料圆柱塌落的影响；王鹏程等[4-5]利用PFC2d对碎石集料双轴压缩试验进行了数值模拟，研究了颗粒形状对不良级配碎石集料剪切特性的影响；邹德高等[6]利用三维激光扫描技术量化了堆石料三维形状特征，并研究了其对颗粒破碎的影响。已有研究成果表明，不同颗粒形状有着不同的历史成因，颗粒之间的嵌挤和摩擦受颗粒形状的影响，使整体强度发生变化[7-10]。对于颗粒形态的量化，国内外学者很早就开展了一系列研究。Barrett[11]提出颗粒形态特征包括颗粒形状、棱角性以及颗粒表面纹理等三个独立的特性，这三个特性分别代表颗粒在大、中、小三个尺寸上的变化规律。棱角度作为一种中尺度形态表征参数，最早由Wadell[12]于1932年提出的圆度(Roundness)进行评价并被广泛应用于工程实践中。Masad[13]等基于等效椭球提出了棱角度AI(Angular index)以评价颗粒表面棱角起伏程度，Zhou[14]等基于球谐函数对该等效椭球体的计算方法进行了改进。

但已有研究对于颗粒形状的描述大多还处于二维平面层面，且形状较为规则，缺乏精细模型的三维模拟。本文通过PFC3d软件建立了直接剪切试验条件，利用软件内置的clump算法生成了不规则碎石颗粒精细三维模型，同时基于文献[15]对于颗粒形状的研究，改进了棱度指数的概念，并将其作为形状评价指标，定量研究了颗粒形状对颗粒材料剪切特性的影响。

1 单颗粒三维模型及三维形状指标

1.1 颗粒三维模型的建立

为了获取原始颗粒的三维形状信息，利用CT扫描获得颗粒逐层二维CT影像[16-18]，再将CT影像导入MATLAB中，利用对影像进行二值化、边界检测处理后生成的像素信息矩阵提取颗粒轮廓点的相对坐标，从而得到颗粒表面的点云集。将点云导入三维建模软件建立颗粒三维模型，将其轮廓进行Delaunay 三角化后，利用PFC3d中的clump 算法对其进行填充得到单颗颗粒(图1)。

图1 试验颗粒模型Fig.1 Particle model of the experiment

1.2 颗粒三维形状指标

棱角度AI由Masad[13]等基于颗粒等效椭球提出，其计算公式如下：

(1)

式中：t——检测步；

rp——颗粒在球坐标上的极半径；

rEE——等效椭球体在极坐标的极半径。

该等效椭球体为体积等效椭球体，即椭球体体积等于颗粒体积，并且该椭球体长短径比与原颗粒一致。

本研究考虑到在假定颗粒不会破损的情况下，颗粒体的运动主要为平移和旋转，其旋转受外力以及其本身的转动惯量张量控制，因此颗粒等效椭球体的重心体积主惯性力矩应与颗粒一致，即在其长轴、中轴、短轴方向上的转动惯量应与颗粒体一致。对此改进了颗粒等效椭球体的计算公式。

以模型形心为原点建立空间直角坐标系，假设模型轮廓由n个三角形面组成，连接原点和每个轮廓点，则可将模型分为n个四面体。

(2)

图2 第i个四面体Fig.2 The tetrahedron i

由于计算转动张量需要用三重积分：

(3)

但四面体的边界在三维直角坐标系中表达太复杂，所以导入参数方程对其进行简化：

(4)

其边界条件Ω为：

(5)

将式(4)代入式(3)进行三重积分换元：

(6)

(7)

式中：detAi——矩阵Ai的特征值。

将式(7)代入式(6)得到换元后的三重积分公式：

(8)

假设颗粒为均质材料且密度为ρ，则第i个四面体对原点的惯性张量Bi计算式如下：

(9)

式(9)中：

颗粒惯性张量B：

(10)

假设过原点任意轴的方向向量(单位向量)为(α，β，γ)，则该轴的转动惯量I表达式为：

(11)

利用式(11)可求出过原点任意旋转轴的转动惯量，最大转动惯量旋转轴方向对应等效椭球体短轴方向，最小转动惯量旋转轴方向对应等效椭球体长轴方向，中轴方向同时垂直于长轴和短轴，可用长轴方向向量与短轴方向向量叉乘得出，并将中轴方向向量代入式(11)求出相应转动惯量。

等效椭球体长半轴长a、中半轴长b、短半轴长c可由以下方程组求解：

(12)

式中：Imin，Imax——颗粒绕过形心任意转轴的最小、最大转动惯量；

Imid——颗粒绕等效椭球体中轴的转动惯量；

M——椭球体的质量。

计算等效椭球体各轴长及方向后，将直角坐标系转化为球坐标系，采用式(1)计算颗粒棱度AI。

为了方便表述，将利用转动张量矩阵得到的等效椭球体简称为张量椭球体，将利用体积等效得到的等效椭球体简称为体积椭球体。

用2种方法分别求出了颗粒b的等效椭球体，并作出三维图形进行对比(图3)。分别计算了2个等效椭球体与实体颗粒的体积之比(V/VT)、表面积之比(S/ST)、长短径比之比(EI/EIT)、重心体积主惯性力矩之比(I1/I1T、I2/I2T、I3/I3T)，结果如表1所示。

图3 颗粒b等效椭球体对比图Fig.3 Comparison diagram of the equivalent ellipsoid of the particle b(注：自左向右依次为三维图形的俯视图、左视图、正视图以及三维视图，其中黑色线框为颗粒b轮廓面，蓝色线框为张量椭球体，红色线框为球谐函数椭球体。)

表1 两种等效椭球体参数比较

张量椭球体因与颗粒重心体积主惯性力矩等效，所以其I1/I1T、I2/I2T、I3/I3T值都为1，颗粒与其张量椭球体体积的差距在4%以内，表面积和长短径比的差距在20%以内。体积等效椭球体因与颗粒体积等效且长短径比也相同，所以其V/VT、EI/EIT值都为1，颗粒与其体积等效椭球体表面积的差距也在20%以内，略小于张量椭球体，3个方向上惯性力矩的差距都超过了10%，最大差距超过30%。

综上所述，2种等效椭球体与实体颗粒体积的差距、表面积的差距相差不大，虽然体积等效椭球体的长短径比与实体颗粒一致，但张量椭球体更加符合颗粒的力学行为。

图 4为5种不同形状的颗粒及其棱度指数，可以看出颗粒AI值越小，颗粒就越接近于一个椭球体，其轮廓也越规则。表2为5种颗粒各形状指标，其中SI为球度即颗粒等体积球体的表面积与颗粒实际表面积的比值，EI为伸长率即颗粒短轴长与长轴长的比值，e为扁平率即颗粒短轴长与中轴长的比值。

图4 不同形状颗粒及其棱度指数Fig.4 Different particle shapes and its AI

表2 5种颗粒各形状指标

2 PFC3d直剪实验数值模拟

2.1 参数标定

为了标定数值模拟的参数，进行了粒径为 1～2 cm的砂岩碎石在围压为0.1，0.2，0.3 MPa 条件下的室内中型直接剪切试验，砂岩的抗压强度为 30 MPa。试样的高度为16 cm、宽度为14 cm、长度为16 cm。试样密度为2 650 kg/m3，孔隙比约为0.65。

同时，利用PFC3d[19]建立了碎石集料的中型直接剪切试验条件，用10个墙体模拟直剪盒。为了尽可能模拟颗粒间的镶嵌咬合作用，在室内试验试样中选取10个形状较为典型的颗粒进行扫描建模，并导入PFC 中利用clump算法对其进行填充后作为生成颗粒的刚性簇模板，再由这10个模板等量地随机生成粒径范围为1～2 cm形状不规则的碎石，最终的数值模拟模型如图5所示。

图5 混合颗粒离散元模型Fig.5 Discrete element model of the mixed particle

模拟过程中颗粒之间的接触模型采用线性接触刚度模型，需要定义接触的切法向刚度及摩擦系数。表3为本文所选取碎石颗粒的主要计算参数。为了与室内试验进行对比，所有数值模型的终止剪切应变均为12%，剪切速率也与室内试验保持一致。

图6为室内试验以及数值模拟计算得到的应力-应变对比曲线，与室内试验结果对比，可以看出数值模型的曲线要比室内试验曲线略高。在室内试验过程

表3 颗粒计算参数

图6 室内试验和数值计算得到的应力-应变对比曲线Fig.6 Comparison of stress-strain of the laboratory tests and numerical calculations

中，粒间摩擦会产生细小粉尘降低粒间摩擦系数以及发生颗粒破碎的现象，而这一过程在数值模型中并没有得到体现。因此，数值模型得到的碎石剪切强度要比室内试验得到的结果略高。总体上来说，数值模型试验计算结果基本可以描述曲线的主要形态，验证了模型及计算参数选取的合理性。

2.2 单一颗粒形状直剪试验数值模拟

在标定完参数后，对图 4所示的 5种颗粒逐一进行单一颗粒形状直剪试验数值模拟，采用的接触模型及相关参数、直剪盒尺寸、颗粒粒径、剪切速率与参数标定时一致。

图7为数值模拟试验得到的在0.3 MPa围压下不同颗粒碎石集料的应力-应变曲线，可以看出，不同形状颗粒粒组的应力-应变曲线在形态上有很大的不同。从整体上看，曲线为非线性，初始剪切阶段线性较为明显，而随着剪切应变的增加，颗粒集料剪应力在应变达到10% 后均没有太大的变化。

图7 不同颗粒形状下应力-应变数值模拟试验曲线Fig.7 Stress-strain of the numerical simulation test under different particle shapes

根据摩尔-库伦准则整理得到不同形状颗粒的抗剪强度线(图 8)，抗剪强度随围压的增加显著增加。颗粒内摩擦角即抗剪强度线的斜率随颗粒形状的变化也明显不同。

图8 不同颗粒形状下围压与抗剪强度关系Fig.8 Relationship between the confining pressure and shear strength under different particle shapes

图9为 0.3 MPa围压下5个颗粒各个形状参数与抗剪强度的关系，从图9中可看出，颗粒棱度AI较其他3个形状指标与抗剪强度有更好的对应关系。

图9 不同形状指标与抗剪强度关系Fig.9 Relationship between the different shape indexes and shear strength

在验证了模型的正确性后，对另外11种块石颗粒进行建模，分别计算了颗粒的棱度指标。再分别导入PFC3d生成clump模板进行直剪试验后，用试验结果对比颗粒棱度进行分析。

在围压相同的情况下，抗剪强度随棱度的增大而增大，在棱度大于0.3后，抗剪强度增加的幅度逐渐减小(图 10)，可用对数函数进行拟合：

(13)

图10 棱度指数对抗剪强度的影响Fig.10 Effect of AI on the shear strength

计算出每个块体的内摩擦角后，与颗粒棱度指标进行对比，内摩擦角随棱度的增加而增加，增加幅度随棱度的增加而减小(图 11)，其关系也可用对数函数进行拟合，拟合关系式为：

φ=7.621lnAI+63.653

(14)

3 结论

(1)提出了用块体转动张量确定块体的等效椭球体的方法，将其与球谐函数确定方法进行比较，论证了其合理性，并在此基础上改进了颗粒三维棱度指标的计算公式。

(2)由碎石颗粒直剪试验应力-应变曲线可得出，随着应变的不断增加，剪应力不断增加，但增加的速率逐渐趋近于零，有明显的极值。曲线的形态和极值跟颗粒形状有很大的关系。

(3)级配不良碎石集料的抗剪强度主要受其颗粒棱度SI的影响，随着颗粒棱度的增加，其抗剪强度和内摩擦角都明显增加，其关系可以用对数函数进行拟合，内摩擦角与颗粒棱度自然对数拟合直线的斜率为7.621、截距为63.653。在棱度较低的时候，抗剪强度和内摩擦角增加较快，当棱度较大时，抗剪强度指标增加幅度较小。

(4)本研究中的室内试验和数值模拟试验样品颗粒级配不良，对于级配良好颗粒形状的评价还应结合其他指标深入研究。