级配与颗粒形状对复杂堆积体路基填料剪切性能影响的离散元模拟研究

华文俊¹，肖源杰^{1, 2}，王萌¹，郑俊星³，王小明¹，陈晓斌^{1, 2}

(1. 中南大学 土木工程学院，湖南 长沙，410075；

2. 中南大学 重载铁路工程结构教育部重点实验室，湖南 长沙，410075；

3. 爱荷华州立大学 土木建造与环境工程系，爱荷华州 埃姆斯市，50011)

摘要：为探究复杂堆积体填料用于路基填筑的可行性，基于计算机视觉实现堆积体颗粒的三维外形重构，获取颗粒的几何形状参数(棱角指数和长细比)、粒径及其分布；通过开展不同竖向应力水平下的大型直剪试验，研究级配和颗粒形状对堆积体填料剪切性能的影响规律，进而构建基于真实颗粒形状的堆积体填料直剪试验的离散元模型，研究颗粒运动、接触网络、组构异性和配位数等微细观特征的演化规律，建立颗粒微细观参数与堆积体填料宏观剪切行为之间的关联机制。研究结果表明：棱角指数和长细比大的堆积体填料具有更高的抗剪强度，中值粒径与抗剪强度呈明显的正相关；粗颗粒含量多、级配良好的堆积体填料，颗粒间的相互作用随棱角指数的增大而增强，颗粒间的相互错动随颗粒长细比的增大而加剧；残余强度与峰值强度之比随棱角指数的增大、长细比的减小而增大；压密后的密实程度随棱角指数的减小、长细比的增大而增大，初始变形模量和剪胀量也更大。

关键词：隧道施工弃渣；离散单元法；级配；颗粒形状；直剪试验

中图分类号：U416.1                    文献标志码：A

文章编号：1672-7207（2021）07-2332-17

Discrete element modeling(DEM) study on effect of gradation and morphology on shear strength behavior of rock debris as embankment fill materials

HUA　Wenjun¹, XIAO　Yuanjie^{1, 2}, WANG　Meng¹, ZHENG　Junxing³, WANG　Xiaoming¹, CHEN　Xiaobin^{1, 2}

(1. School of Civil Engineering, Central South University, Changsha 410075, China;

2. Key Laboratory of Engineering Structure of Heavy Railway, Central South University, Changsha 410075, China;

3. Department of Civil, Construction and Environmental Engineering, Iowa State University, Ames 50011, USA)

Abstract: In order to explore the feasibility of rock debris used for subgrade filling, three-dimensional shape reconstruction of particles was realized based on computer vision, and the geometric shape parameters including angularity index and flat and elongated ratio, particle size and its distribution were obtained. The effects of gradation and particle morphology on the shear strength behavior of rock debris were studied from large-scale direct shear tests with different levels of normal stress. Further, the discrete element model(DEM) of laboratory direct shear tests was established by using element shape library digitized from real particles and calibrated from laboratory test results. The micro-mechanical characteristics including particle movement, inter-particle contact force network, fabric evolution, and coordination number during the shearing process were simulated and linked to the macro-scale shear behavior. The results show that particles of rock debris with larger angularity index and flat and elongated ratio exhibit greater shear strength. A significant positive correlation was found between the median diameter and shear strength. Specimens of rock debris with a larger median diameter and wider size distribution show stronger inter-particle interaction with the increase of angularity index and greater tendency of particle sliding with the increase of flat and elongated ratio. The ratio of residual shear stress to peak shear stress increases with increasing angularity index and decreasing flat and elongated ratio. Specimens of rock debris with smaller angularity index and larger flat and elongated ratio show better compaction, greater levels of initial deformation modulus, and shear-induced volume dilatancy.

Key words: rock debris from tunneling; discrete element method; gradation; particle morphology; direct shear tests

在交通强国国家战略的驱动下，公路工程建设逐渐向地质条件复杂的西部山区转移，西部山区高差显著，公路建设线路中隧道占比大，方量巨大的隧道施工弃渣可形成厚度巨大的复杂堆积体；此外，高山冲洪积物等也可形成深厚复杂堆积体。受显著地形高差及选线限制，部分线路被迫建在复杂堆积体上，高速公路对路基提出了较高的平顺性和稳定性要求，如何将这些复杂堆积体填料用于公路粒料基层的填筑以及如何确保复杂堆积体上修筑的路面结构的服役性能是亟待研究的课题。粒料基层作为路面结构的必要组成部分，可充分传递和扩散车辆荷载引起的应力，起到保护下层软弱土基的作用。大量尺寸和形状各异的复杂堆积体填料颗粒随机分布在粒料基层内部，使其表现出明显的非连续性、非均匀性和各向异性。在外荷载作用下，堆积体填料颗粒之间的摩擦、嵌挤和咬合等作用十分复杂，从微细观角度深入研究复杂堆积体填料的剪切特性有助于揭示其内在机制和劣化机理。直接剪切试验相较于其他复杂的试验方法而言操作简便且成本更低，故广泛用于确定粒料基层填料等散粒体材料的抗剪强度特性^[1]。

传统的基于连续介质力学理论的有限单元法无法反映粒料基层在颗粒尺度上的细观力学特性，且难以考虑颗粒级配、尺寸和形状等参数对宏观力学性能的影响。近年来，离散单元法^[2-3]已成为求解散粒体内部应力与变形问题的有效工具，广泛用于模拟室内试验和揭示散粒体颗粒系统的力学特性。采用这一非连续介质分析方法可以很好地反映粒料基层的散体特性，从微细观尺度模拟颗粒之间的相互作用，从而研究粒料基层填料受力变形的内在机制。

复杂堆积体填料颗粒通常具有形状不规则、棱角尖锐和表面纹理粗糙等形态特征，已有研究表明颗粒的级配、棱角度和长细比会显著影响散粒体材料的强度、模量和变形特性。徐旸等^[4]提出考虑细颗粒含量的级配分形公式，通过模拟不同级配道砟的直剪试验发现级配分维度和细颗粒含量与道砟抗剪强度之间存在较强的负相关性。MAHMUD等^[5]采用离散元法研究了特定级配范围内不同的粗细比水平对道砟抗剪强度的影响。王鹏程等^[6]利用PFC3D研究了颗粒粒径对碎石集料抗剪强度的影响，发现随着粒径的增大，集料的抗剪强度增大，内摩擦角减小。赵仕威等^[7]利用开源离散元程序YADE研究了颗粒棱角度在直剪试验模拟中的影响，发现竖向荷载越大，颗粒棱角度的影响越明显。MISHRA等^[8]研究了颗粒形状对道砟抗剪强度的影响，发现增大颗粒棱角度和长细比均可显著增强道砟的抗剪强度。BIAN等^[9]利用自行开发的块体离散元程序BLOCK3D对不同形状道砟的直剪试验进行了模拟，发现试样的剪切强度和体积变形随着颗粒长细比和棱角度的增大而增大。然而，上述研究多采用随机生成的黏结球簇或理想凸多面体模拟颗粒外形，未能与真实颗粒外形建立联系或过高估计颗粒的棱角性；单个数值试样仅由一种特定形状的颗粒组成，与实际情况不符，仍需进一步从微细观角度研究剪切作用下不同级配堆积体填料的宏观力学性能机理。

为揭示堆积体内单个颗粒的性质与复杂堆积体整体力学行为之间的关联机制，本文作者采用计算机视觉成像技术获取单个颗粒的几何形态和尺寸^[10]，构建与真实颗粒具有相同几何特征的离散单元，结合室内大型直剪试验和离散元模拟手段研究不同竖向荷载作用下不同级配和颗粒形状的复杂堆积体填料的剪切力学行为，着重分析复杂堆积体填料在剪切过程中颗粒运动和接触状态等微细观特征的演化规律及其与剪切强度的关联机制。

1  堆积体填料颗粒几何形态分析

本文采用计算机视觉技术从3个互相垂直的方向对堆积体填料粗颗粒进行拍照^[11-14]，精确地获取堆积体内颗粒的真实三维外形(图1(a))，经过阈值化处理后得到二值图像，如图1(b)所示。采用RAO等^[13]提出的针对单个颗粒棱角度和长细比的计算方法，将投影轮廓简化为n边形(图1(b))，计算n边形相邻顶点对应内角的变化值，并将其位于18个长度为10°的等值区间内的期望值作为棱角度：

图1　基于计算机视觉的堆积体颗粒形状分析与三维重建示意图^[15]

Fig. 1　Schematic diagram of particle shape analysis and three-dimensional reconstruction based on computer vision^[15]

(1)

式中：为区间的左端点；为n边形相邻内角的变化值，；为位于区间内的概率。

按照式(2)计算3个视图下棱角度的加权平均值，即为颗粒的棱角指数：

(2)

式中：为视图编号；为投影的面积。

和分别为颗粒投影在2个互相垂直方向上的最大和最小尺寸，按照式(3)计算3个视图下颗粒投影的最大尺寸与最小尺寸的比值(最小包围盒的长宽比)，即为颗粒的长细比：

(3)

式中：为视图中颗粒投影的最大尺寸；为视图中颗粒投影的最小尺寸。

本文从级配编号为2的堆积体填料试样中随机挑选78个粗颗粒，采用前述计算机视觉方法对大粒径颗粒或者利用集料图像分析仪对中小粒径颗粒的形态特征进行量化分析，得到该组堆积体填料试样的平均棱角指数和平均长细比分别为410和2.54。选取6种具有代表性形状的堆积体颗粒，计算其形态特征参数，并得到3个视图下的颗粒投影；继而在三维空间中，将简化为n边形的投影轮廓沿投影的反方向进行拉伸，经布尔运算得到重建的三维颗粒^[15]，如图1(b)所示。

2  堆积体填料室内大型直剪试验

本文所选取的试验材料为洗净的火成岩材质堆积体碎石，相对密度为2.32。在进行室内大型直剪试验之前，人为剔除4组堆积体填料试样中的大粒径颗粒以便消除尺寸效应的影响后，再进行筛分，结果如图2所示。

图2　试样的颗粒级配曲线

Fig. 2　Gradation curves of specimens

为研究不同级配对堆积体填料剪切性能的影响并建立定量的关联，本文采用式(4)所示的Logistic函数来定量描述各组堆积体填料试样的级配特征。

(4)

式中：为堆积体填料颗粒的累计过筛质量分数；为衡量曲线分布宽度的形状参数，越大表示曲线越窄，粒径分布越均匀；为中值粒径，mm。

4组级配曲线拟合参数，的取值和拟合优度见表1。

表1　不同级配曲线的拟合参数

Table 1　Parameter fitted for different gradation curves

直剪试验装置由长×宽×高分别为305 mm×305 mm×76 mm和356 mm×305 mm×152 mm的上、下剪切盒组成^{[1, 16-18]}。试验开始前，将大约20 kg的堆积体填料分3层装入剪切盒内，每层均振动击实至密实状态，控制试样的总高度为203 mm，其中下剪切盒试样高152 mm，上剪切盒试样高51 mm，初始孔隙率为0.35。试验过程中，竖向应力经由压缩气体充气气囊均匀地施加在加载板上，加载板随上剪切盒沿水平方向运动，剪切速率保持在12.2 mm/min，下剪切盒固定不动，当剪应力几乎不变时停止测试。

3  直剪试验的离散元模拟

为研究级配和颗粒形状对堆积体填料颗粒剪切性能的影响，并更好地理解其在剪切作用下宏观力学响应的微细观机制，本文的离散元模型采用clump单元模拟堆积体颗粒的真实三维形状，通过生成12面墙体来模拟实际的剪切盒，且暂未考虑颗粒破碎的影响。参考MAHMUD等^[5]和MISHRA等^[8]的研究，颗粒之间、颗粒与墙体之间的接触采用线性刚度模型进行模拟，基于室内直剪试验数据并通过试错法对离散元模型中堆积体颗粒的微细观接触参数进行标定，标定后的微细观力学参数见表2。

表2　离散元模型中的微细观力学参数

Table 2　Micro-mechanical parameters used in DEM simulation

LIM等^[19]指出当单个簇颗粒所包含的球体数量大于8个时，就能反映颗粒间的多点接触和咬合关系。赵云哲^[20]认为当簇颗粒外包络网格填充率大于90%时，再增加球形单元的数目对计算结果影响不大。因此，本文的离散元模型采用以上文献标准控制堆积体颗粒三维形状的模拟精度，从堆积体颗粒三维模型库^[21]中选取了编号为Ⅱ，Ⅳ，Ⅴ，Ⅷ，Ⅸ和Ⅺ的典型颗粒，采用bubblepack算法^[22]对三维几何边界进行填充，得到表3所示的簇颗粒单元。对比真实颗粒的不规则外形与离散元模型的簇颗粒单元可以看出，各模型库下的簇颗粒单元可以较好地模拟对应颗粒的复杂外形。

表3　离散元模型所采用的簇颗粒单元

Table 3　Clump elements used in DEM simulation

为保证离散元模型的计算效率，将颗粒的最小粒径设置为9.5 mm，按照试验级配和初始空隙率生成堆积体颗粒，采用落雨法^{[20, 23]}分层填装堆积体颗粒，得到不同级配的堆积体填料的虚拟直剪试样如图3(a)所示。为研究不同颗粒形状对堆积体填料剪切性能的影响，本文设计4组不同棱角指数和长细比组合的堆积体填料虚拟试样(图3(b))，其中各编号簇颗粒单元的体积分数由多目标规划算法求得，即形状1：I_A=410，R_FE=2.54；形状2：I_A=591，R_FE=2.54；形状3：I_A=591，R_FE=3.66；形状4：I_A=410，R_FE=2.54。

图3　堆积体填料的离散元虚拟直剪试样

Fig. 3　Virtual direct shear specimens for DEM simulations

为实现竖向应力的施加，引入数值伺服机制控制加载板的运动，通过不断调整顶部约束速度使得约束应力保持不变。当加载板上的竖向应力水平达到目标值158 kPa时，测得压实后试样的空隙率如表4所示。从表4可以看出：在同一竖向应力水平作用下，棱角指数、长细比以及级配曲线形状参数b较小的试样更容易被压密。

表4　压实前后试样的空隙率

Table 4　Porosity of specimens before and after compaction

对上剪切盒的五面墙体和两面翼墙沿长度方向施加恒定的剪切速率，而下剪切盒的五面墙体保持静止以实现水平加载，当水平位移达到对应于15%的剪应变水平时停止加载。通过分析直剪过程中的剪切应力和体积变形可评价堆积体填料试样的剪切力学性能，试样的体积应变由加载板的竖向位移确定，剪切面上的切向力为

(5)

式中：为堆积体填料颗粒与上剪切盒内壁的接触力沿剪切方向的分量。由于上剪切盒由五面墙体组成，故n=5。

剪切应力为

(6)

式中：为下剪切盒的长度；为剪切速率；为剪切过程时间。

4  结果与分析

4.1　离散元数值模拟与室内试验对比分析

对不同竖向应力条件下的直剪试验进行离散元模拟，将模拟结果与试验结果进行对比。当离散元模型所设置的微细观参数可以吻合竖向应力为158 kPa时的试验结果，并较准确地预测竖向应力为248 kPa和317 kPa的试验结果时，认为所选取的离散元模型参数具有合理性，模型能准确地模拟堆积体填料在剪切作用下的力学响应。

图4所示为不同竖向应力水平条件下，级配编号为1的堆积体填料试样的剪切应力-剪切应变曲线的模拟与试验结果。从图4可以看出：在3种不同的竖向应力水平下，离散元模拟得到的剪切应力-剪切应变曲线与试验结果较为吻合，呈现出先硬化后稍微软化的特点。这说明试样在剪切前已经压实稳定，即初始试样较为密实；试样的剪切强度随着竖向应力的增加而不断增大，当竖向应力分别为158，248和317 kPa时，离散元模拟得到的抗剪强度分别为200，309和415 kPa，与试验结果相近，表明随着竖向应力的增加，试样在压实过程中达到的密实程度更高，其初始变形模量和峰值剪切应力也越大。由于室内直剪试验过程中竖向应力不断变化，而离散元模拟采用伺服机制使得竖向应力基本不变，故模拟结果与试验结果存在一定差异。

图4　不同竖向应力水平下剪切应力-剪切应变曲线的室内试验与离散元模拟结果对比

Fig. 4　Comparison of laboratory-measured and DEM-simulated shear stress-shear strain curves under different normal stress levels

4.2　级配对堆积体填料剪切性能的宏观影响

采用校准后的离散元模型对158 kPa竖向应力作用下的其余3种不同级配的堆积体填料试样进行模拟，分别得到剪切应力-剪切应变曲线和竖向应变-剪切应变曲线如图5和图6所示。

图5　不同试样级配下剪切应力-剪切应变曲线的离散元模拟结果

Fig. 5　DEM simulation results of shear stress-shear strain curves for different specimen gradations

图6　不同试样级配下竖向应变-剪切应变曲线的离散元模拟结果

Fig. 6　DEM simulation results of vertical strain-shear strain curves for different specimen gradations

从图5可以看出：级配曲线形状参数b相近的试样具有相似的应力-应变曲线；试样的抗剪强度与其中值粒径d₅₀呈明显的正相关，与级配曲线形状参数b则无明显相关性；试样达到峰值应力状态时对应的剪切应变随中值粒径d₅₀的增大而增大，与级配曲线形状参数b无明显相关性。这说明较高的细颗粒含量会明显削弱堆积体填料的抗剪强度，粗颗粒含量高的堆积体填料抵抗剪切作用和塑性流动的能力更强。

从图6可以看出：堆积体填料试样在剪切初始阶段的竖向应变保持为0，说明试样在剪切前已压实稳定；随着剪切应变的增加，试样的竖向应变逐渐增大并最终表现出明显的剪胀现象。级配4的堆积体填料试样在剪切作用下先体缩后体胀，原因可能是当剪切应变较小时，大量的细颗粒填充了粗颗粒间因相互错动、重新排列而形成的空隙，使得颗粒堆装结构更稳定从而剪胀量减小。此外，堆积体填料的最大剪胀量与中值粒径d₅₀呈明显的正相关，与级配曲线形状参数b呈较强的负相关，说明粒径较小且分布均匀的堆积体填料抵抗体积变形的能力更强。

4.3　级配影响堆积体填料剪切性能的微细观机制

以竖向应力为158 kPa的试验工况模拟结果为例，从颗粒运动状态(位移、转角)和接触特征(力链网络、法向接触力组构、配位数)2个方面对堆积体填料在直剪过程中的微细观力学特性进行分析。

图7所示为在158 kPa竖向应力作用下，当剪切应变为15%时，剪切盒内不同试样级配下的颗粒位移矢量在x-z平面内的分布情况(包括沿剪切方向的水平位移和不同高度处竖向位移的分布情况)。图7中箭头的长度与位移成正比关系，箭头的颜色对应位移的范围，箭头指向位移矢量的方向。4组不同试样级配下的颗粒平均位移分别为12.3，10.7，10.7和10.3 mm。可以看出，填料颗粒在剪切作用下的体积变形具有不均匀性：上剪切盒内的颗粒位移较大，且以向右上方的位移为主，与试样的剪胀相对应；下剪切盒内的颗粒位移较小，且以向右的位移为主，右侧颗粒有向下运动的趋势，与试样的剪缩相对应；上、下剪切盒内颗粒位移的差异凸显出试样内部的弧形滑动面，颗粒位移与中值粒径d₅₀成正比关系，与级配曲线形状参数b无明显关系。

图7　不同试样级配下颗粒的位移矢量

Fig. 7　Displacement vectors of particles for different specimen gradations

图8所示为在158 kPa竖向应力作用下，当剪切应变为15%时，剪切盒内不同试样级配下的颗粒转角在x-z平面内的分布情况(包括不同角度范围的转角分布情况)。图8中气泡的直径与转角成正比，气泡的颜色对应旋转角度的范围。可以看出，不同级配的填料颗粒绕Y轴转角的分布情况大致相同，上剪切盒内位于剪切面附近的颗粒其转角较大，与图7所示的位移矢量分布情况相似。由于颗粒与墙体之间的摩擦因数较小，靠近剪切盒壁面的颗粒更易发生转动。随着中值粒径d₅₀的增大，细颗粒含量降低，转角在0°~10°范围内的颗粒所占比例也随之降低。在颗粒形状相同的条件下，由于大粒径颗粒的转动惯量较大，故粗颗粒含量较高的填料试样不容易发生颗粒转动，平均转角较小。随着级配曲线形状参数b的增大，细颗粒与粗颗粒含量同时降低，转角在0°~10°范围内的颗粒占比随之降低，而转角在20°~30°范围内的颗粒所占比例明显升高，故填料颗粒的平均转角增大。这说明在级配不良的堆积体填料内部，颗粒之间的相互咬合嵌挤作用较弱，颗粒更容易发生转动。

图8　不同试样级配下颗粒绕Y轴的转角

Fig. 8　Particle rotation angle around Y-axis for different specimen gradations

图9所示为在158 kPa竖向应力作用下，当剪切应变分别为0，5%，10%和15%时，不同试样级配下的颗粒接触网络在x-z平面内的分布情况。其中，为粒间接触力的平均值；为滑动接触百分比。图9中行反映接触网络中的强/弱接触，接触力大于平均接触力的接触为强接触，反之为弱接触；图9中行反映接触网络中的滑动/非滑动接触，切向接触力大于法向接触力与摩擦因数乘积的接触为滑动接触，反之为非滑动接触。图9中接触力与力链的长度成正比关系。可以看出，粒间接触力的分布呈现出极大的不均匀性，级配曲线形状参数b相近的填料试样具有相似的接触网络；对于级配良好的填料试样，强力链主要集中在上剪切盒左部至下剪切盒右部的对角线区域；颗粒之间的平均接触力随中值粒径d₅₀和级配曲线形状参数b的增大而增大，说明在粗颗粒含量高且粒径分布均匀的试样内部，颗粒之间具有较强的相互作用；滑动接触百分比随中值粒径d₅₀的增大而增大，随级配曲线形状参数b的增大而减小。从图7、图8和图9可知：对于粗颗粒含量较高且粒径分布不均匀的堆积体填料，颗粒之间的相互错动更加剧烈，试样更容易发生整体剪切破坏。

图9　不同试样级配下的颗粒接触网络

Fig. 9　Inter-particle contact force network for different specimen gradations

为反映颗粒体系组构(例如接触网络及其空间分布)的各向异性，可采用ROTHENBURG等^[24-25]建议的傅里叶函数对法向接触力的分布进行描述：

(7)

式中：为接触法向落在该角度区间的法向接触力的平均值；为所有粒间法向接触力的平均值；为傅里叶拟合系数，反映法向接触力各向异性程度；为法向接触力各向异性主方向。

图10所示为在158 kPa竖向应力作用下，当剪切应变分别为0，5%，10%和15%时，不同试样级配下的粒间法向接触力在x-z平面内的分布情况。图10中极角为接触法向与X轴的夹角，极径为接触法向落在各角度区间内的法向接触力的平均值，虚线为式(7)中傅里叶函数的拟合曲线，箭头指向为法向接触力各向异性的主方向。从图10可以看出：不同剪切应变下粒间法向接触力近似呈“花生状”分布；未发生剪切时，不同试样级配下的粒间法向接触力各项异性的主方向均近似为90°，沿竖直方向的法向接触力明显大于沿水平方向的法向接触力；随着剪切应变的增加，粒间法向接触力不断增大，法向接触力方向由竖直方向偏转至水平方向。剪切前后4组不同试样级配下的法向接触力各项异性系数的变化值分别为0.11，-0.02，-0.11和-0.17。可见，级配曲线形状参数b越大，法向接触力各向异性主方向的偏转角度越小，法向接触力方向的各向异性变化程度越大；粗颗粒含量高(中值粒径d₅₀大)且级配良好(级配曲线形状参数b小)的颗粒体系内部接触状态良好，在外界荷载作用下法向接触力各向异性主方向的偏转角度更大。

图10　不同试样级配下粒间法向接触力的分布

Fig. 10　Distribution of inter-particle normal contact force for different specimen gradations

配位数是颗粒料内部每个颗粒与周围颗粒的接触数目，受空隙率、级配和颗粒形状等因素的影响，可用于反映颗粒空间分布状态与颗粒间接触密度等微细观特征。平均配位数可由式(8)计算得到：

(8)

式中：为试样的平均配位数；为颗粒与颗粒之间的接触数目；为颗粒总数。

图11所示为当竖向应力为158 kPa时，不同试样级配下的平均配位数随剪切应变的变化曲线。从图11可以看出：随着剪切应变的增加，试样内颗粒的配位数不断减小，且存在一定的波动，说明剪切过程中剪切面积的减小与试样体积的膨胀使得颗粒之间的挤压作用变弱，颗粒之间的接触密度减小；剪切前后4组不同试样级配下的平均配位数的变化值分别为0.17，0.58，0.84和0.42。从图5、图6和图11可知：对于级配曲线形状参数b较大的试样，其残余强度与峰值强度之比更大，剪胀量更小，颗粒骨架越稳定。

图11　不同试样级配下的配位数

Fig. 11　Coordination number for different specimen gradations

4.4　颗粒形状对堆积体填料剪切性能的宏观影响

通过离散元模拟得到竖向应力为158 kPa时，不同颗粒形状下的剪切应力-剪切应变曲线和竖向应变-剪切应变曲线分别如图12和图13所示。由图12可以看出：对于棱角指数较小的堆积体填料试样，其剪切应力-剪切应变曲线具有明显的应变软化特征；试样的抗剪强度与颗粒的棱角指数和长细比呈现出较强的正相关性；试样达到峰值应力状态时对应的剪切应变随棱角指数和长细比的增大而增大，这说明棱角指数和长细比较大的堆积体填料抵抗塑性流动的能力更强。

图12　不同颗粒形状下剪切应力-剪切应变曲线的离散元模拟结果

Fig. 12　DEM simulation results of shear stress-shear strain curves for different particle shape categories

图13　不同颗粒形状下竖向应变-剪切应变曲线的离散元模拟结果

Fig. 13　DEM simulation results of vertical strain-shear strain curves for different particle shape categories

由图13可以看出：试样的最大剪胀量随着棱角指数的增大而减小，说明颗粒棱角能显著增强颗粒间的“互锁效应”；试样的最大剪胀量随颗粒长细比的增大而增大，其原因可能是剪切作用下长细比大的颗粒相较水平方向的滑动(摩擦消耗的能量更多)更倾向于沿竖直方向运动(竖向位移明显增大)。

4.5　颗粒形状影响堆积体填料剪切性能的微细观机制

图14所示为在158 kPa竖向应力作用下，当剪切应变为15%时，剪切盒内不同颗粒形状下的颗粒位移矢量在x-z平面内的分布情况(包括剪切方向水平位移和不同高度处竖向位移的分布情况)。4组不同颗粒形状下的颗粒平均位移分别为10.7，7.8，10.3和11.8 mm。从图14可以看出：颗粒位移随棱角指数的增大而减小，说明颗粒棱角可以增强颗粒间的相互咬合嵌挤作用；颗粒位移随长细比的增大而增大，滑动面(竖向位移为0的等值面)的相对位置也随之向下移动，说明长细比大的颗粒在剪切作用下更容易沿竖直方向运动。由于本文未能考虑颗粒破碎的影响，考虑到长细比越大的颗粒在剪切作用下也越容易破碎或断裂，故上述规律仍有待进一步验证。

图15所示为在158 kPa竖向应力作用下，当剪切应变为15%时，剪切盒内不同颗粒形状下的颗粒转角在x-z平面内的分布情况。从图15可见：不同颗粒形状下颗粒绕Y轴转角的 分布情况大致相同，颗粒的平均转角随棱角指数的增大而减小；随着长细比的增大，颗粒的平均转角仅略微减小，转角在0°~10°范围内的颗粒占比却明显升高，这说明棱角指数对颗粒转动的影响更大，而增大长细比可以有效抑制颗粒的剧烈转动。又因为在垂直平面内，长细比大的颗粒发生单位角度的转动对试样高度的影响更大，故颗粒长细比越大，试样的最大剪胀量也越大。

图14　不同颗粒形状下的颗粒位移矢量

Fig. 14　Displacement vectors of particles for different specimen gradations

图15　不同颗粒形状下颗粒绕Y轴的转角

Fig. 15　Particle rotation angle around Y-axis for different particle shape categories

图16所示为在158 kPa竖向应力作用下，当剪切应变分别为0，5%，10%和15%时，不同颗粒形状下的颗粒接触网络在x-z平面内的分布情况。从图16可知：当棱角指数和长细比均较小时，粒间接触力的分布相对均匀；随着棱角指数的增大，强力链主要集中在上剪切盒左部至下剪切盒右部的对角线区域，颗粒间的相互作用更加强烈，使填料试样具有更高的抗剪强度和抗变形能力，而棱角指数与滑动接触百分比无明显关系；随着颗粒长细比的增大，滑动接触网络由上剪切盒左侧向下剪切盒右侧发展，滑动接触百分比明显增大。从图14、图15和图16可知：对于长细比大的堆积体填料，颗粒之间的相互错动更加剧烈，试样更容易发生整体剪切破坏。

图16　不同颗粒形状下的颗粒接触网络

Fig. 16　Inter-particle contact force network for different particle shape categories

图17所示为在158 kPa竖向应力作用下，当剪切应变分别为0，5%，10%和15%时，颗粒间法向接触力在x-z平面内的分布情况。剪切前后4组不同颗粒形状下的粒间法向接触力各项异性系数的变化值分别为-0.02，0.06，0.16和-0.18。从图17可以看出：对于棱角指数较大的试样，粒间法向接触力各项异性主方向的偏转角度更大，法向接触力方向的各项异性程度变大；对于颗粒长细比较大的试样，法向接触力各项异性主方向的偏转角度更小，法向接触力方向的各项异性变化程度更大。这说明棱角度较大和长细比较小的颗粒体系内部接触状态良好，在外界荷载作用下法向接触力各向异性主方向的偏转角度更大。

图17　不同颗粒形状下粒间法向接触力的分布

Fig. 17　Distribution of inter-particle normal contact force for different particle shape categories

图18所示为当竖向应力为158 kPa时，不同颗粒形状下的平均配位数随剪切应变的变化曲线。剪切前后4组不同颗粒形状下的平均配位数的变化值分别为0.58，0.38，0.55和0.72。从图12、图13和图18可知：对于棱角指数较大和颗粒长细比较小的试样，其残余强度与峰值强度之比更大，剪切初始阶段的剪胀量更小，颗粒骨架更稳定；对于棱角指数较小或长细比较大的试样，其压实后的颗粒配位数更大，密实程度更高，初始变形模量和最大剪胀量也更大。

图18　不同颗粒形状下的配位数

Fig. 18　Coordination number for different particle shape categories

5  结论

1) 堆积体填料试样的抗剪强度随中值粒径、颗粒长细比和棱角指数的增大而增大；试样的剪胀量随级配曲线形状参数和棱角指数的增大而减小，随中值粒径和颗粒长细比的增大而增大；粗颗粒含量高、棱角指数和长细比大的堆积体填料抗塑性流动能力更强。

2) 堆积体填料的颗粒位移随中值粒径和棱角指数的增大而减小，随颗粒长细比的增大而增大，增大粗颗粒含量、棱角指数和长细比以及控制级配良好可有效抑制填料内部的颗粒转动。

3) 随着级配曲线形状参数、中值粒径和棱角指数的增大，强力链的分布更为密集，颗粒之间的相互作用更加强烈；随着中值粒径和颗粒长细比的增大以及级配曲线形状参数的减小，滑动接触百分比明显增大，颗粒之间的相互错动更加强烈，堆积体填料更容易发生整体剪切破坏。

4) 粗颗粒含量高、级配良好、棱角指数较大和长细比较小的堆积体填料试样内部的粒间接触状态良好，残余强度与峰值强度之比更大；棱角指数较小和颗粒长细比较大的试样压实后的密实程度更高，初始变形模量和剪胀量也更大。

参考文献：

[1] WNEK M A, TUTUMLUER E, MOAVENI M, et al. Investigation of aggregate properties influencing railroad ballast performance[J]. Transportation Research Record, 2013, 2374(1): 180-189.

[2] CUNDALL P A. A computer model for simulating progressive, large-scale movement in blocky rock system[C]// Proceedings of the International Symposium on Rock Mechanics. Rotterdam, Holland: AA Balkama, 1971: 8-12.

[3] CUNDALL P A, STRACK O D L. A discrete numerical model for granular assemblies[J]. Géotechnique, 1979, 29(1): 47-65.

[4] 徐旸, 高亮, 王红, 等. 道砟级配的分形方法及对道床剪切性能影响[J]. 铁道学报, 2016, 38(12): 94-101.

XU Yang, GAO Liang, WANG Hong, et al. Study of fractal method and influence of ballast gradation on ballast bed shear behavior[J]. Journal of the China Railway Society, 2016, 38(12): 94-101.

[5] MAHMUD S M N, MISHRA D. Development of 'coarse to fine ratio' as a new gradation parameter to study ballast packing[C]// Proceedings of 2018 Joint Rail Conference. Pittsburgh, PA, USA : The American Society of Mechanical Engineers, 2018: 50978.

[6] 王鹏程, 刘建坤, 李旭, 等. 直剪条件下碎石集料力学特性研究[J]. 铁道学报, 2014, 36(2): 103-108.

WANG Pengcheng, LIU Jiankun, LI Xu, et al. Mechanical behavior of crushed stone aggregate subjected to direct shear conditions[J]. Journal of the China Railway Society, 2014, 36(2): 103-108.

[7] 赵仕威, 周小文, 刘文辉, 等. 考虑颗粒棱角影响的直剪试验的离散元模拟[J]. 岩土力学, 2015, 36(S1): 602-608.

ZHAO Shiwei, ZHOU Xiaowen, LIU Wenhui, et al. DEM simulation of direct shear tests considering particle angularity effect[J]. Rock and Soil Mechanics, 2015, 36(S1): 602-608.

[8] MISHRA D, MAHMUD S M N. Effect of particle size and shape characteristics on ballast shear strength: a numerical study using the direct shear test[C]// Proceedings of 2017 Joint Rail Conference. Philadelphia, PA, USA: The American Society of Mechanical Engineers, 2017: 50718.

[9] BIAN Xuecheng, LI Wei, QIAN Yu, et al. Micromechanical particle interactions in railway ballast through DEM simulations of direct shear tests[J]. International Journal of Geomechanics, 2019, 19(5): 04019031.

[10] TUTUMLUER E, HUANG H, HASHASH Y M A, et al. Imaging based discrete element modeling of granular assemblies[C]// 9th International Conference on Multiscale and Functionally Graded Materials. Oahu, HI, USA: American Institute of Physics, 2008: 544-549.

[11] RAO C, TUTUMLUER E. Determination of volume of aggregates: new image-analysis approach[J]. Transportation Research Record, 2000, 1721(1): 73-80.

[12] PETERSEN D R, LINK R E, RAO C, et al. Coarse aggregate shape and size properties using a new image analyzer[J]. Journal of Testing and Evaluation, 2001, 29(5): 461-471.

[13] RAO C, TUTUMLUER E, KIM I T. Quantification of coarse aggregate angularity based on image analysis[J]. Transportation Research Record, 2002, 1787(1): 117-124.

[14] PAN Tongyan, TUTUMLUER E, ANOCHIE-BOATENG J. Aggregate morphology affecting resilient behavior of unbound granular materials[J]. Transportation Research Record, 2006, 1952(1): 12-20.

[15] HUANG Hai, LUO Jiayi, TUTUMLUER E, et al. Size and shape determination of riprap and large-sized aggregates using field imaging[R]. Urbana: Illinois Center for Transportation, 2020: 17-27.

[16] DOMBROW W, HUANG Hai, TUTUMLUER E. Comparison of coal dust fouled railroad ballast behavior-granite vs. limestone[C]// 8th International Conference on the Bearing Capacity of Roads, Railways and Airfields. Champaign, IL, USA: Taylor & Francis Group, 2009: 1349-1357.

[17] BOLER H. On the shear strength of polyurethane coated railroad ballast[D]. Urbana: University of Illinois at Urbana-Champaign, 2012: 16-20.

[18] DERSCH M S, TUTUMLUER E, MARENO A P, et al. Characterization of polyurethane coated railroad ballast behavior[C]// Transportation Research Board 90th Annual Meeting. Washington D C, USA: Transportation Research Board, 2011: 1-17.

[19] LIM W L, MCDOWELL G R. Discrete element modelling of railway ballast[J]. Granular Matter, 2005, 7(1): 19-29.

[20] 赵云哲. 基于道床非线性特性的桥上有砟轨道系统动力学特性研究[D]. 北京: 北京交通大学, 2018: 24-30.

ZHAO Yunzhe. Research of dynamic characteristic of ballast track on bridge based on nonlinear property of ballast bed[D]. Beijing: Beijing Jiaotong University, 2014: 24-30.

[21] HUANG Hai. Discrete element modeling of railroad ballast using imaging based aggregate morphology characterization[D]. Urbana: University of Illinois at Urbana-Champaign, 2010: 39-51.

[22] TAGHAVI R. Automatic clump generation based on mid-surface[C]// 2nd International FLAC/DEM Symposium. Melbourne, Australia, 2011: 791-797.

[23] 王子杰. 基于离散单元法的道砟力学特性研究[D]. 北京: 北京交通大学, 2015: 31-47.

WANG Zijie. Railway ballast mechanical mechanism research based on discrete element method[D]. Beijing: Beijing Jiaotong University, 2015: 31-47.

[24] ROTHENBURG L, BATHURST R J. Analytical study of induced anisotropy in idealized granular materials[J]. Géotechnique, 1989, 39(4): 601-614.

[25] BATHURST R J, ROTHENBURG L. Observations on stress-force-fabric relationships in idealized granular materials[J]. Mechanics of Materials, 1990, 9(1): 65-80.

(编辑  杨幼平)

收稿日期： 2020 -11 -10； 修回日期： 2021 -01 -20

基金项目(Foundation item)：国家重点研发计划项目(2019YFC1904704)；国家自然科学基金资助项目(51878673，U1734208，U1934209)；中国铁道科学研究院科技研究开发计划重大课题(2019YJ026) (Project(2019YFC1904704) supported by the National Key Research Development Program of China; Projects(51878673, U1734208, U1934209) supported by the National Natural Science Foundation of China; Project(2019YJ026) supported by the Key Research Development Program of Chinese Academy of Railway Sciences)

通信作者：肖源杰，博士，副教授，从事路基工程研究；E-mail：yjxiao@csu.edu.cn

DOI: 10.11817/j.issn.1672-7207.2021.07.020

引用格式： 华文俊, 肖源杰, 王萌, 等. 级配与颗粒形状对复杂堆积体路基填料剪切性能影响的离散元模拟研究[J]. 中南大学学报(自然科学版), 2021, 52(7): 2332-2348.

Citation: HUA Wenjun, XIAO Yuanjie, WANG Meng, et al. Discrete element modeling(DEM) study on effect of gradation and morphology on shear strength behavior of rock debris as embankment fill materials[J]. Journal of Central South University(Science and Technology), 2021, 52(7): 2332-2348.