0 引言

碎石集料是一种非均质且各向异性的复杂多相散体岩土材料。由于其具有就地取材方便、压密性能好、抗剪抗震性能强、变形量小等优势，广泛用于交通建设工程中。隧道开挖所产生大量的岩石弃渣被加工处理成为路用碎石，在实现交通绿色发展、低碳环保中起到重要作用。级配良好的路用碎石可作为铁路道砟、公路路基及坝体填料等^[1]。目前针对废料利用后的路用碎石集料的研究，主要集中在探讨其力学特性对公路路基、路面材料等稳定性、强度的影响，而忽略了其颗粒形态、颗粒破碎等自身特征的作用。

颗粒形态作为路用碎石集料固有特征，是影响材料破碎的一个重要因素。颗粒形状通常可以量化为3个独立指标：长宽比、圆度(或棱角度)、粗糙度，其中圆度对颗粒破碎的影响十分显著^[2-4]。不少学者通过单颗粒破碎试验、大型三轴试验等手段开展了圆度对颗粒破碎的影响。Wang等^[5]发现，具有较大局部圆度的岩石颗粒大概率以较高的破碎强度呈爆炸式破裂，具有较小局部圆度的岩石颗粒更容易发生劈裂。Todisco等^[6]发现，棱角明显的颗粒容易以较低的强度破裂。显然，局部应力倾向于集中在较尖锐的棱角处，从而使颗粒更容易破碎。尽管使用各类先进手段(如X射线断层扫描、高速显微镜摄像机或光弹性技术)可以获得大量与颗粒破碎相关的数据，以揭示颗粒形状对颗粒破裂的影响；但由于碎石集料固有的分散性和各向异性，以及试验设备的限制，对颗粒尺度上的破碎机理认识仍不充分。因此，离散元法被应用于颗粒材料力学特性的研究中。随着离散元法的发展，引入了多球颗粒模型和多面体模型，以探究颗粒形态的影响。Nie等^[7]通过球谐函数拟操作方法生成不同圆度指标颗粒模型，并基于离散元法采用多球刚体模型探究了颗粒圆度对颗粒材料宏观力学特性的影响。同时，通过细观尺度变量和组构异性演化揭示了宏观变化的细观根源。Fang等^[8-9]基于离散元法，采用黏结的Voronoi胞元分析了各类颗粒形貌指标(扁平度、伸长度、棱角度)对易碎填料破碎特征和力学性能的影响规律，并基于细观机理揭示了宏观现象发生的原因。现有文献在开展圆度对路用碎石集料的力学特性影响研究时，仅将颗粒视为刚性体，而忽略了颗粒的破碎效应。若直接通过3D扫描获取不同圆度的颗粒模型，既费时费力，又很难排除其他形状指标的影响；而若采用逆球谐函数法^[10]生成目标颗粒模型，当颗粒的圆度较小时，模型会出现重叠、交叉等畸变。

本研究依托瓮马铁路南北延长线梅家巷隧道项目，收集隧道弃渣的灰岩处理后的路用碎石样本，针对其3D扫描模型提出构建不同圆度的不规则颗粒模型的新方法。利用离散元软件PFC3D开展了路用碎石的三轴试验数值仿真研究，探究了破碎效应下颗粒圆度对碎石集料破碎特征和宏细观力学特性的影响规律。研究结果以期为隧道弃渣处置为路用碎石集料的选取提供可靠的试验和理论依据，具有重要工程应用价值和推广前景。

1 数值模拟过程 1.1 颗粒形状指标

本研究将待破碎的碎石集配颗粒通过黏聚在一起的胞元进行模拟，胞元形状即可为球形，也可以为简单多面体或Voronoi多面体。由于Voronoi多面体能保证破碎过程中颗粒质量守恒，且更高效精准地反映颗粒轮廓，更真实再现裂隙发育，因此本研究选择将不规则颗粒模型剖分为Voronoi多面体集合体，并对颗粒形状指标进行定义。颗粒的伸长度E₁和扁平度F₁可量化针状颗粒和片状颗粒，分别定义为：

(1)

(2)

式中P₁，P₂，P₃分别为OBB边界框最大、中间和最小边的长度^[11]，OBB指能够容纳一个颗粒的最小尺寸的长方形体盒。E₁值越小，F₁值越小。将颗粒圆度R_d定义为：

(3)

式中，R_i为第i个棱角对应的局部内切圆半径；n为颗粒棱角的总数；R_ins为颗粒轮廓最大内切圆半径。R_d越大，颗粒越接近球体，反之则越不规则。

1.2 颗粒模型生成

利用三维扫描技术，基于圆度对颗粒形态进行重构。首先，选取50个碎石集料，通过三维扫描获得颗粒不规则轮廓；其次，基于主成分分析法获取该颗粒模型的伸长度E_I与扁平度F_I；之后，通过P₂/E₁和P₃/E₁F₁使得颗粒在3个维度的长度相等，排除伸长度和扁平度效应；最后，采用球谐阶数为15的球谐函数法重构颗粒模型的表面轮廓，重构的每个颗粒模型由1 536个三角面片组成。此时的模型能够有效地去除表面粗糙度，并很好地识别颗粒棱角，并可排除网格数量对颗粒破碎的影响^[12-14]。

评估处理完毕的颗粒模型圆度指标，并每隔0.1颗粒圆度挑选1个圆度指标作为研究变量。碎石中不存在圆度极大的颗粒模型，该模型可通过球谐函数法或建模软件获取。最终确定的6个颗粒模型及相应的形状指标如表 1所示，所有模型仅考虑圆度效应。

1.3 等向固结与三轴试验

Fang等^[15]对粒径2 cm碎石(粒径与本研究相同)开展单颗粒破碎室内试验，获取单颗粒刚度和破碎强度。通过相应的数值模拟标定了组成单颗粒数值模型的胞元细观刚度和强度范围，发现强度在1.5~10.2 MPa能满足2 cm碎石的真实破碎强度分布特征。考虑本研究重点为颗粒形态的影响，即唯一变量为颗粒圆度，在不考虑颗粒破碎强度的变化下，颗粒间黏结参数统一设置为5 MPa。朱艳贵^[16]和龚建^[17]统计了近年来国内外文献探究颗粒材料力学行为的离散元模拟所用细观参量，确定了颗粒间接触参数的取值范围，并对取值合理性进行了详细说明。因此，本研究离散元法模拟中，颗粒参数参考文献[15-17]中的离散元模型细观参数，如表 2所示。

对颗粒进行Voronoi剖分，将各颗粒模型剖分为150个Voronoi多面体，具体的剖分过程见文献[15]。随后对颗粒模型开展等向压缩。本研究选用单一粒径为2 cm的碎石集料，基于离散元法模拟三轴压缩排水试验。将设定方案中选用的颗粒模型作为样本，在被6个刚性墙面所围成的立方体空间内，生成待破碎颗粒。将颗粒设置为相互不接触，且方向均匀随机分布。为平衡计算效率和模型尺寸效应，等向压缩固结所用颗粒数量为1 000个。为保证最终模型的各向同性，也为防止压缩过程中产生应力梯度，将颗粒与颗粒间、颗粒与墙体间的摩擦系数及模型重力场暂时设置为0。本研究主要聚焦颗粒破碎对剪切强度的影响，为防止固结过程中产生破碎行为，暂时将细观黏结参数提高1 000倍。以极低速度移动六面墙体，对试样进行等向压缩，并采用伺服控制机制达到目标围压为0.5 MPa。当目标围压与实际应力值的容差低于0.1%，且平均静态不平衡力与平均接触力的容差小于0.01%时，试样达到静力平衡态，等向压缩固结完成。模拟试验前，墙体与颗粒间摩擦系数不变，颗粒间的摩擦系数重置为0.5。此外，黏结胞元的细观强度参数恢复至目标值。

所有参数调整完毕后，进行常规三轴压缩排水试验的数值模拟。试验过程中，固定底部墙体，同时以极低速率移动顶部墙体实现对试样的轴向压缩，并保持试样四周围压不变。当剪切速率极低时才能忽略试验过程中的动能，认为是准静态加载。一般采用无量纲惯性参数判断是否准静态加载。顶墙以0.1 m/s速率压缩，此时惯性参数≤1.5×10^―3，满足准静态加载要求。为突显颗粒破碎行为对集料力学属性影响，同步进行不考虑颗粒破碎的模拟试验。

2 试验结果分析 2.1 颗粒破碎特征分析

根据颗粒破碎前后级配曲线变化，Hardin^[18]提出颗粒破碎程度衡量指标，即相对破碎率B_r，定义为：

(4)

式中，B_t为颗粒破碎总量，等于初始级配曲线、破碎后颗粒分布曲线及截止粒径三者的所围区域面积；B_p为颗粒破碎势，作为衡量颗粒在破碎过程中的能量消耗和破碎程度的指标，等于初始级配曲线、100%含量线与截止粒径三者的所围区域面积。Hardin^[18]认为颗粒粒径越小，其破碎强度越高，越不易发生破碎。当粒径小于0.074 mm时，颗粒难以再破碎，所以截止粒径取值为0.074 mm。不同颗粒圆度影响下，路用碎石集料的相对破碎率随轴向应变的演化规律如图 1所示。

由图 1可知，相对破碎率随轴向应变的增加而增大。这表明持续的轴向加载和颗粒圆度的减小均会导致更多的颗粒破碎发生。这是因为圆度更小的颗粒更不规则，会产生更为明显的互锁效应^[19]，颗粒间容易产生应力集中，从而引发更多的颗粒破碎。对图 1结果进行数据拟合，拟合公式为：

(5)

式中，k为拟合曲线的斜率；ε₁为轴向应变。轴向应变与相对破碎率呈线性增加趋势，相关系数R²均大于0.91，相关性良好。

由图 1统计得到拟合曲线斜率随圆度的变化，如图 2所示。由图可见，斜率随颗粒圆度增大，呈现单调递减趋势，且二者可较好地拟合为线性关系曲线(R²=0.96)，这表明在剪切过程的任意应变阶段，均可认为相对破碎率和圆度线性相关。同时，二者的关系也可较好地拟合为对数形式(R²=0.99)。Fang等^[9]和周剑等^[19]在研究棱角度效应时均发现了拟合曲线斜率与棱角度的类似关系，这表明颗粒圆度与棱角度在促使颗粒破碎效果方面极为相似。

2.2 宏观力学特性 2.2.1 宏观力学参数定义

路用碎石的宏观力学参数可有效反映材料的宏观力学特性，对于三轴试验而言，主要力学参数包括平均应力p和偏应力q, 计算为：

(6)

(7)

式中，σ₁为轴向应力；σ₂和σ₃分别为中间主应力和最小主应力。颗粒集合体的轴向应变ε₁和体应变ε_v分别计算为：

(8)

(9)

式中，H和V分别为给定变形状态下装配体的高度和体积；H₀和V₀分别为剪切前装配体的初始高度和体积。基于摩尔库伦强度准则，集料内摩擦角可通过式(10)计算得到：

(10)

2.2.2 剪切强度变化

为探究路用集料的临界剪切状态特性，确定了三轴试验中试样产生轴向应变ε₁=40%时应终止试验，不同颗粒圆度下应力比随轴向应变的变化如图 3所示。

由图 3可知：(1)初始加载阶段，应力比随应变增大而急剧上升，在经历较小应变后达到应力比峰值，然后随着应变持续增加而逐渐减小，最终趋于稳定，即可视为临界状态。(2)圆度越大，与峰值应力比对应的应变越小。峰后应力比的下降趋势说明试样呈现出较为明显的应变软化。原因是若不考虑颗粒破碎，则试样在剪切过程中出现明显的剪胀现象，对应的峰值应力比更大，峰后试样处于高势能状态，颗粒的翻滚和颗粒间的相互挤压，导致了峰后的应变软化。(3)相同颗粒圆度下，考虑颗粒破碎的试样峰值应力比小于不考虑破碎的试样。

峰值摩擦角(φ_p)和临界摩擦角(φ_c)随颗粒圆度的变化曲线，如图 4所示。由于当轴向应变≥30%，试样应力比值几乎恒定不变，此时可认为试样处于临界状态。因此取轴向应变≥30%范围内应力比的平均值并通过式(10)计算临界摩擦角。

无论考虑破碎与否，峰值摩擦角和临界摩擦角均随圆度的增加而降低。对峰值摩擦角和临界摩擦角进行数据拟合，可表示为：

(11)

式中，φ为峰值摩擦角或临界摩擦角；A和B为拟合参数。峰值摩擦角和临界摩擦角随圆度的增加整体呈线性衰减趋势。但是，峰值摩擦角值会随圆度的增加略有上升，随后下降；而临界剪切强度与圆度呈现出很好的线性相关性(R²≥0.88)。在其他文献中也观察到类似的趋势^[7]。对比是否考虑颗粒破碎的结果，发现颗粒破碎明显降低了碎石集料的峰值摩擦角，且圆度越小，二者的差值越大。这说明颗粒破碎对集料的剪切峰值强度产生明显影响。无论颗粒破碎与否，碎石集料的临界摩擦角接近，这说明颗粒破碎对临界剪切强度影响较小。

2.2.3 碎石集料的剪胀性变化

不同颗粒圆度下体应变随轴向应变的变化如图 5所示。由图可知：(1)无论是否考虑颗粒破碎，碎石集料的体应变均表现为先剪缩后剪胀，最终趋于相对稳定状态。最大剪胀量发生在圆度为0.5时，体应变分别达到―1.78%和―1.17%。最大剪缩量发生在圆度为0.7时，体应变分别达到19.5%和16.8%。(2)存在临界圆度为0.7，使得集料体积应变的最大剪缩量体应变随圆度先增大而后减小。(3)考虑颗粒破碎的试样体应变始终要小于不考虑颗粒破碎的试样，在较低圆度下二者相差更明显。这可解释为低圆度下，颗粒破碎严重，破碎产生的小粒径碎块随之填充在大颗粒粒径间的空隙处，减小了体应变的持续发展。

统计不同圆度下考虑颗粒破碎与否2种情况下试样的临界体应变值，如表 3所示。基于表 3绘制临界体应变随圆度的变化曲线，如图 6所示。

由图 6可知，无论颗粒破碎与否，存在临界圆度R_cd=0.7，使得当R_d＜R_cd时，试样临界体应变随圆度增加逐渐增大；当R_d≥R_cd时，试样临界体应变随圆度增加而减小。对比不考虑颗粒破碎的情况，颗粒破碎引起的体积收缩将导致引入破碎效应的碎石集料临界体应变明显小于不考虑颗粒破碎的集料。此外，更小的圆度对应颗粒破碎程度更大。考虑颗粒破碎与不考虑颗粒破碎情况下临界体应变的差值Δε_d随圆度增加而减小，如图 7所示。Δε_d与圆度的关系式可拟合为：

(12)

式中A和B分别为拟合系数。拟合曲线的相关系数R²=0.986，这说明Δε_d与圆度表现出良好的线性相关，即Δε_d随圆度呈线性下降趋势。在颗粒圆度较高时，考虑破碎的颗粒集料其临界体应变较不考虑破碎时要小，该现象是颗粒破碎导致的体积收缩与破碎后不规则碎块挤压、翻越和填充共同作用的结果。

3 细观力学特性 3.1 接触力演化

外部荷载由颗粒系统内部的接触进行传递，并形成接触力。接触力影响颗粒的破碎和运动状态，进而导致材料的力学特性变化。不同圆度影响下，试样平均接触力随轴向应变的演化规律如图 8所示。由图可知，无论颗粒破碎与否，平均切向接触力始终要大于平均法向接触力，这与文献[20-21]观察到的现象一致。此外，随着剪切的进行，平均切向和法向接触力均表现出前期快速增加，随后保持不变或轻微下降的趋势。

考虑与不考虑颗粒破碎情况下，临界状态接触力随圆度的变化如图 9所示。由图可知，当不引入颗粒破碎效应时，集料的平均法向和切向接触力随圆度的增加先快速下降，随后基本保持不变。这表明圆度越小，集料发生颗粒破碎的潜力越大。与之相比，当引入破碎效应后集料的2类接触力要明显更小。颗粒破碎是出现此类差异的原因。颗粒破碎会导致碎屑含量增加，并消除颗粒间由于互锁效应产生的应力集中现象。更小的颗粒圆度对应更为严重的颗粒破碎，由此产生大量的细颗粒碎屑。已有研究表明^[21]小粒径颗粒间的接触会产生较小的接触力。另一方面，颗粒破碎减弱了颗粒间的互锁效应，因而应力集中现象也随之减少。因此，引入破碎后，产生小接触力的可能性增大。此外，颗粒破碎会破坏原有颗粒系统的承载结构，导致颗粒间接触位置的荷载重新分布。因为颗粒破碎缓解了颗粒间的应力集中，所以颗粒间的受力分布更加均匀。基于颗粒破碎对碎屑含量、互锁效应和应力集中的影响，并且由于圆度越小的颗粒发生破碎的可能性越大，因此，随着圆度增加，考虑与不考虑颗粒破碎的集料二者接触力的差值会变得越来越小。

3.2 滑动接触演化

滑动接触百分比是评价颗粒运动状态的重要表征量。基于库伦摩擦定律，颗粒发生滑动与否可通过式(13)判断。

(13)

式中，S_c为接触应力比；f_s和f_n分别为接触切向力和接触法向力；μ_b为颗粒之间的摩擦系数。若S_c≥0.999 9，则认为颗粒间为接触滑动。滑动接触百分比为滑动的接触总数与系统内总的接触数量的比值。

剪切过程中不同颗粒圆度下滑动接触百分比随轴向应变的变化如图 10所示。由图可知，无论颗粒破碎与否，路用碎石集料的滑动接触占比在整个加载过程中的演化规律相同, 即随着轴向应变增加，滑动接触百分比首先在小应变阶段急剧上升至峰值，随后逐渐下降至相对稳定状态。由图 10(a)可知，较小圆度的颗粒其形貌更加不规则，颗粒间更易于发生互锁效应，导致颗粒系统更加稳定。因此，当不引入破碎效应时，集料的滑动接触占比随颗粒圆度增加而降低。图 10(a)与图 10(b)对比发现，引入破碎效应后，集料的滑动接触百分比在低圆度条件下相比不考虑颗粒破碎情况发生下降。推测原因是在加载过程中不考虑破碎的颗粒一般通过接触滑动进行重排，而当引入破碎效应后，若施加的应力达到或超过颗粒破碎强度，则部分颗粒存在破碎而非滑动的可能，因而其滑动接触百分比相比不考虑颗粒破碎时更低。

根据图 10得到临界状态下滑动接触占比随圆度的变化，如图 11所示。由图可知，对于不考虑颗粒破碎的集料，临界状态的滑动接触百分比随圆度增加逐渐下降。另外，考虑颗粒破碎的集料滑动接触百分比小于不考虑颗粒破碎的情况，且随圆度减小，二者差距增加。若颗粒无法破碎，则在颗粒间发生滑动重排。但对于有可能发生破碎的颗粒(颗粒的外部荷载达到颗粒破碎极限)，则更倾向于发生破碎而非滑动，导致其更低的滑动接触占比。当圆度较高时，虽然颗粒破碎会引起滑动接触占比降低，但破碎产生的碎片使得集料整体不规则度升高。根据已有研究，不规则度的升高会导致滑动接触占比的提升^[22]。因此，高圆度下，颗粒破碎与不规则形状碎片二者的共同作用，导致颗粒破碎与否对滑动接触占比的影响不明显。

4 结论

本研究基于扫描的真实颗粒形貌构建了不同圆度的不规则颗粒模型，并采用黏结Voronoi多面体方法模拟待破碎颗粒。基于离散元法等向固结6种不同圆度的碎石集料，并生成密实试样。随后，开展一系列常规三轴压缩排水试验的数值仿真研究，探究颗粒圆度对路用碎石集料破碎特征和宏细观力学特性影响。

(1) 更小的颗粒圆度对应更严重的颗粒破碎程度，相对破碎率随轴向应变增加而线性增大，且拟合的线性函数斜率与圆度线性相关。

(2) 随着圆度增加，峰值剪切强度先达到峰值，随后快速降低，临界剪切强度线性递减。颗粒破碎导致了更低的峰值剪切强度。

(3) 无论是否考虑颗粒破碎，碎石集料的体应变均表现为先剪缩，然后逐渐剪胀，最终趋于相对稳定状态。存在临界圆度为0.7，使得集料体积应变的最大剪缩量随圆度先增大而后减小。临界状态体应变先随圆度增加快速达到峰值，当圆度大于临界圆度时，又快速下降。此外，临界体应变的差值随圆度呈线性减小的趋势。

(4) 随着剪切进行，平均接触力表现为先快速上升后缓慢下降或保持不变的趋势，更小的圆度对应的接触力更大。颗粒不破碎情况下滑动接触占比随圆度增加下降，颗粒破碎情况下滑动接触占比与圆度的规律性一样。此外，颗粒破碎导致同等圆度下集料的平均接触力和滑动接触占比小于不考虑颗粒破碎的集料，且圆度越小，差异越大。