颗粒形状对堆石料力学特性影响的离散元分析

王蕴嘉，宋二祥，张千里

(1. 中国铁道科学研究院集团有限公司铁道建筑研究所，北京 100081；2. 清华大学土木工程系，北京 100084)

堆石料广泛应用于土石坝、机场高填方路基和铁路公路填方路基等填筑工程中[1-2]。堆石料取材便利、来源广泛，不同工程中采用的堆石料在颗粒强度、颗粒形状等方面存在较大差异。大量工程经验表明，颗粒形状是影响堆石料力学性能的重要因素之一，值得进行深入研究。

近年来，很多学者采用室内试验方法研究颗粒形状对堆石料、砂等散粒体材料力学特性的影响[3-13]。Chuhan 等[3]对砂土进行侧限压缩试验，结果表明砂颗粒棱角越鲜明、越易发生破碎，材料的屈服应力越低。Varadarajan 等[4]分别对砂卵石堆石料和爆破堆石料进行相似级配试验，发现二者表现出相反的试验规律，进而认为颗粒形状对材料破碎程度有较大影响。Afshar 等[5]采用电镜扫描的方式统计了大量岩石颗粒的形状特征并对其进行单颗粒压碎试验，发现颗粒球度越高越难破碎。邹德高等[6]对单一粒组堆石料进行三轴试验，研究颗粒颗粒形状对颗粒破碎率的影响，发现堆石料颗粒破碎率随颗粒球度的增加而减小。刘钢等[7]对大量碎石颗粒进行统计分析，提出了颗粒几何尺寸的量化表征参数，认为颗粒棱角与纹理是影响堆石料孔隙率的重要因素。Lashkari等[8]对六种不同形状的砂进行直剪试验，分析了颗粒形状对试样临界状态的影响，发现形状越不规则，临界状态摩擦角越高。此外，田继荣等学者[9-12]也进行了大量试验研究，分析颗粒形状对散粒体力学特性的影响。现阶段室内试验测得的大多是材料宏观力学表现，难以捕捉颗粒细观表现，不利于进行细观机理分析，且天然材料形状不规则难以控制单一变量，定量研究较困难。

堆石料是典型的颗粒材料，越来越多的学者采用多尺度分析方法研究其宏、细观力学特性[13-14]。离散单元法由Cundall 和Strac[15]提出，该方法将介质视为一系列离散的独立运动单元，可方便的处理非连续介质力学问题，常被用于研究颗粒形状对散粒材料力学特性影响及其细观机理[16-24]。Ng[16]对长短轴比值不同的椭球颗粒试样进行三轴剪切试验模拟，发现椭球颗粒长短轴之比越大、材料抗剪强度越小。边学成等[17]采用自行开发的三维块体离散元程序建立道砟颗粒模型，模拟分析道砟颗粒料在直剪试验中的细观力学特性。王蕴嘉和宋二祥[18]采用PFC3D 模拟研究球度对砂卵石堆积特性及强度的影响，并分析其细观机理。Xu 等[19]采用三维离散元模拟分析颗粒形状对试样力学特性的影响，发现颗粒形状越不规则，材料剪切强度越大，且颗粒形状对临界状态线的影响不可忽略。Katagiri 等学者[20-24]也分别就考虑颗粒形状的数值模型构造和形状影响机理开展研究。

与一般土体相比，堆石料粒径大、颗粒破碎明显，研究中应考虑颗粒破碎的影响[25]。基于此，常晓林等[26]采用SGDD 模型，分析颗粒形状及粒间摩擦角的影响，并建立粒间摩擦角与非线性强度参数的对应关系。Zhou 等[27-28]采用FEM与DEM 相结合的模拟方法研究颗粒形状的影响，结果表明颗粒形状越狭长、颗粒破碎越显著。Wang 等[29]模拟分析了颗粒形状对堆石料蠕变特性的影响，发现堆石料颗粒棱角越尖锐、蠕变速率越大。Zhang 等[30]对不同颗粒形状试样进行模拟试验，发现不规则颗粒试样会产生更多的破碎，且颗粒破碎与能量输入的关系受颗粒形状影响。可破碎材料的力学性能与颗粒强度密切相关，目前对综合考虑颗粒形状和颗粒强度影响的研究尚不充分，且尚未建立颗粒形状参量与试样宏观力学性能参数的对应关系。

本文以爆破堆石料为研究对象，采用离散元软件PFC3D 根据真实颗粒形状建立不同球度的单颗粒模型，分别生成不破碎和可破碎试样并进行三轴剪切试验模拟，分析颗粒球度对堆石料力学特性的影响及细观机理。

1 数值模型及试验方案

1.1 单颗粒模型

Barrett[31]对已有的岩石颗粒形状相关研究进行综述，提出颗粒形状一般可从宏观、中观、细观三个尺度来描述。宏观尺度一般指颗粒的整体形状，例如正方形或三角形；中观尺度指颗粒局部变化的剧烈程度，例如角部圆滑或尖锐；细观尺度是指颗粒表面形态，例如光滑或粗糙。本文暂不考虑堆石料颗粒的中观、细观形态，仅从宏观尺度出发讨论颗粒形状对堆石料力学特性的影响。参考文献[18]，采用球度S作为描述颗粒宏观形态的计算参数。

爆破堆石料由于爆破开采的原因，颗粒往往棱角分明，这使其在受力过程中应力集中更明显，颗粒破碎现象十分显著。对于爆破堆石料而言，颗粒破碎是影响材料力学特性的重要因素，模拟时需考虑颗粒破碎。基于此，采用可模拟颗粒破碎的颗粒簇单元来构建单颗粒模型。

如图1 所示，依据真实颗粒形状建立不同球度S的爆破堆石料单颗粒模型。单颗粒模型由半径相同的球颗粒粘结而成，球颗粒按密排六方结构排列。综合考虑单颗粒模型的精细程度和计算效率，分别采用13 颗粒、23 颗粒和33 颗粒的颗粒簇作为球度分别为0.878、0.617 和0.446 的单颗粒模型，不同模型的圆度及表面粗糙度均相同。

图1 真实颗粒与单颗粒模型Fig. 1 Rocks and particle model

单颗粒模型内的球颗粒互相不重叠，颗粒真实密度即为计算密度，无需修正。球颗粒之间采用平行粘结模型，当接触所受的拉应力或剪应力超出平行粘结的法向或切向强度时，粘结破坏，此后接触仍采用线性接触模型和滑移模型来模拟。模型的颗粒力学参数如表1 所示。

表1 颗粒力学参数Table 1 Granular mechanical parameters

1.2 单元试验模型

颗粒流模拟试验中，数值试样尺寸过小可能会导致结果不可靠[32]。参考已有研究，综合考虑计算结果可靠性和计算效率，试样边长与颗粒平均直径之比在5～9，最大与最小颗粒粒径之比为1.25。数值试样及试样尺寸参数分别如图2 和表2所示。采用文献[18]提出的制样方法进行制样，制样摩擦系数取0.1，试样相对密实度为0.75。

表2 数值试样参数Table 2 Parameters of numerical models

图2 数值试样(S=0.878)Fig. 2 Numerical sample (S=0.878)

1.3 试验方案

为讨论颗粒破碎的影响，分别生成不破碎爆破料试样和可破碎爆破料试样。两组试样的接触刚度、摩擦系数等颗粒力学参数均相同，仅粘结强度不同。不破碎试样用于模拟颗粒强度较高的情况，可认为材料在常用应力水平范围内几乎不发生破碎，颗粒破碎的影响可忽略不计，此时粘结强度均值取2000 MPa；可破碎试样用于模拟材料强度较低的情况，颗粒破碎显著，其影响不可忽略，此时粘结强度按照表1 取值。

分别对两组试样进行三轴剪切试验模拟，围压取值为400 kPa、600 kPa、800 kPa 和1000 kPa。颗粒球度取值分别为0.878、0.617 和0.446。

2 模拟结果及分析

2.1 三轴剪切试验结果

2.1.1 应力-应变曲线

图3 和图4 为不同围压、不同球度时不破碎试样的应力-应变曲线。由图3 可知，当试样颗粒形状不变时，偏应力随围压增大而增大，试样剪胀性随围压增大而减小，这与已有的试验结果[2]是一致的。由图4 可知，围压为800 kPa 时，峰值偏应力及残余偏应力均随S的减小而增大，试样剪胀性随S的减小而减小。这与已有模拟研究[26]得到的规律相符。

图3 围压对不破碎试样三轴剪切试验表现的影响(S=0.878)Fig. 3 The effect of confining pressure on triaxial compression behavior of unbreakable specimens (S=0.878)

图4 球度对不破碎试样三轴剪切试验表现的影响(800 kPa)Fig. 4 The effect of S on triaxial compression behavior of unbreakable specimens (800 kPa)

图5 和图6 为不同围压、不同球度时可破碎试样的应力-应变曲线。由图可知，围压为400 kPa时，峰值偏应力受球度影响较小，试样剪胀性随S的减小而减小；围压增加至800 kPa 时，峰值摩擦角随S的减小而减小，试样剪胀性随S减小而减小，峰值摩擦角变化规律与400 kPa 围压时并不相同。

图5 球度对可破碎试样三轴剪切试验表现的影响(400 kPa)Fig. 5 The effect of S on triaxial compression behavior of breakable specimens (400 kPa)

图6 球度对可破碎试样三轴剪切试验表现的影响(800 kPa)Fig. 6 The effect of S on triaxial compression behavior of breakable specimens (800 kPa)

对比不破碎试样和可破碎试样的三轴剪切试验结果可知，颗粒球度对试样强度的影响规律与颗粒破碎即颗粒强度有关，颗粒强度高、围压低时，试样强度随S的增大而减小，颗粒强度低、围压高时，试样强度随S的增大而增大；试样剪胀性均随S的增大而增大，且围压和S均相同时，可破碎试样剪胀性小于不破碎试样。

2.1.2 峰值摩擦角

图7 分别为不破碎和可破碎试样的峰值摩擦角。由图可知，不考虑颗粒破碎时(颗粒强度高)，试样峰值摩擦角随围压增大稍有降低，随S增大而减小；考虑颗粒破碎时(颗粒强度低)，峰值摩擦角随围压增大而迅速减小，且S越小峰值摩擦角减小的程度越大。

图7 峰值摩擦角Fig. 7 Peak friction angle

粗粒料峰值摩擦角由两部分组成：一是颗粒间滑动时的滑动摩擦角；二是颗粒间由于咬合产生的咬合摩擦，也可理解为剪胀角。不考虑颗粒破碎时，剪胀角随围压增大而稍有降低(图7(a))；考虑颗粒破碎时，围压增大，颗粒破碎量增加，咬合摩擦即剪胀角显著减小，使得峰值摩擦角随围压增大而显著减小(图7(b))。

S的影响可从以下两方面分析：首先，S越小颗粒不规则程度越高，试样内颗粒咬合、内锁程度越高，会提高试样强度；其次，S越小越容易发生颗粒破碎，尤其是平均压应力水平高时，这会降低试样强度。试样强度的提高或降低取决于两种作用叠加后的效果，因此，在讨论颗粒形状影响时，需综合考虑颗粒破碎的影响。基于以上分析，S对堆石料强度的影响可依据平均压应力水平的不同大致分为三个区域：当平均压应力水平低时，S的影响以增大试样内颗粒咬合、内锁为主，试样强度随S的减小而增大；当平均压应力水平高时，S的影响以增加试样内颗粒破碎为主，试样强度随S的减小而减小；当平均压应力水平适中时，S减小时增大颗粒咬合和增加颗粒破碎的作用相当。

峰值摩擦角和围压在半对数坐标系内有较好的线性关系，围压越高，颗粒破碎越多，咬合摩擦越小。本文在Duncan 等[33]建议的非线性强度公式基础上进行改进，提出可反映球度S影响的峰值摩擦角表达式：

其中：

式中：pa用于平衡量纲，取值一个标准大气压；Sref为参考球度； φ0,ref和 Δφref为S=Sref时的 φ0和Δφ取值；a和b为反映球度影响的拟合参数。

本算例中，取Sref=0.878 ，不同球度时 φ0和Δφ取值如图8 所示，此时各参数取值如表3 所示。

表3 峰值摩擦角公式参数Table 3 Parameters of the peak friction angle formula

图8 S 对峰值摩擦角公式参数的影响Fig. 8 The effect of S on parameters of the peak friction angle formula

需注意的是，半对数坐标系下，试样峰值摩擦角随围压增大而减小，近似线性关系，但减小的斜率绝对值并不是固定的，可大致分为两阶段：当围压低时，颗粒几乎不破碎，峰值摩擦角-围压曲线的斜率绝对值较小；当围压高时，颗粒大量破碎，峰值摩擦角-围压曲线的斜率绝对值较大。因此，若只采用高围压时的峰值摩擦角拟合强度曲线，其曲线斜率绝对值较大，会高估低围压时的试样强度。

2.1.3 初始弹性模量

图9 为S对三轴剪切试验中初始弹性模量E0的影响。由于剪切初始阶段颗粒破碎较少，对初始弹性模量的影响可忽略不计，此处仅给出不破碎试样的初始弹性模量。由图可知，E0随围压增大而增大，且随S增大而增大。

图9 初始弹性模量Fig. 9 Initial tangent modulus

根据计算结果对Janbu 公式[34]进行修正，提出可反映球度S影响的初始弹性模量表达式：

其中：

式中， α1、 α2和n为拟合参数。

对图9 中数据进行拟合，K取值如图10 所示，拟合参数取值如表4 所示。

图10 S 对参数K 的影响Fig. 10 The effect of S on K

表4 初始弹性模量公式参数Table 4 Parameters of the initial tangent modulus formula

2.2 细观组构分析

颗粒材料的细观组构是指颗粒在空间中的分布、定向，孔隙在空间中的分布，颗粒间的相互作用，以及这些因素在空间中的分布情况。细观组构及其演化规律常用于解释和预测颗粒材料的宏观力学行为。本部分针对不破碎和可破碎试样，采用不同的细观组构进行分析。

2.2.1 不破碎试样

采用接触力平均值、拉应力平均值和接触各向异性指标来描述和分析不破碎材料的细观表现[13]。接触力平均值是指试样内所有颗粒间接触的力的平均值，描述试样内颗粒整体受力水平；拉应力平均值是指颗粒内所有接触的拉应力的平均值，描述颗粒内拉应力水平；接触各向异性指标 φ描述试样内颗粒间接触的各向异性程度，可采用下式计算：

式中：C是试样内颗粒间的接触个数；mi1、mi2和mi3分别为第i个接触法向向量的第一、第二和第三主分量。

图11 为S=0.878 时不同围压下的接触力平均值、拉应力平均值和接触各向异性指标曲线。由图11(a)和图11(b)可知，接触力平均值随轴向应变的增加而逐渐增大，且当围压从400 kPa 增大至1000 kPa 时，相同轴向应变时接触力平均值随围压增大而增大；拉应力平均值变化规律与接触力平均值基本相同。以上结果说明，围压越高，试样内接触力水平和颗粒内拉应力水平都越高。

图11 围压对细观组构参量的影响(S=0.878)Fig. 11 The effect of confining pressure on meso-fabric parameters (S=0.878)

由图11(c)可知，接触各向异性指标随轴向应变增大而增大并最终保持不变；围压越低，接触各向异性指标随轴向应变增加而增大的速率越快，这与图3(a)中的应力应变曲线是相对应的，围压越低，试样达到峰值强度时所需的轴向应变越小。

图12 为围压800 kPa 时不同球度下的接触力平均值、拉应力平均值和接触各向异性指标曲线。由图12(a)和图12(b)可知，接触力平均值和拉应力平均值均随着轴向应变的增加而逐渐增大，且当轴向应变相同时，随S减小而增大。以上结果说明，相同应力水平时，颗粒S越小、颗粒间接触力和颗粒内拉应力越大，即颗粒越容易破碎。图12(c)为接触各向异性指标曲线，结合图4(a)分析可知，剪切初期(轴向应变小于5%时)接触各向异性指标增大速率与试样达到峰值强度所需的轴向应变是相关的，球度越大，接触各向异性指标增大速率越大，试样达到峰值强度所需的轴向应变越小。

图12 球度对细观组构参量的影响(800 kPa)Fig. 12 The effect of S on meso-fabric parameters (800 kPa)

2.2.2 可破碎试样

对于可破碎试样，重点研究围压和球度对颗粒破碎的影响，采用相对破碎率表征试样中颗粒破碎程度。相对破碎率Br由Hardin[35]提出，反映破碎前后试样级配曲线的变化，可采用下式计算：

式中：Bt为试验前后试样破碎势之差；Bp为试样初始破碎势。

图13 为S=0.878 时不同围压下的相对破碎率曲线。由图可知，相对破碎率曲线可分为两阶段描述，当轴向应变小于界限应变时相对破碎率为0，当轴向应变大于界限应变时，相对破碎率随应变增大而增大；当轴向应变相同时，相对破碎率随围压增大而增大。

图13 围压对相对破碎率的影响(S=0.878)Fig. 13 The effect of confining pressure on relative breakage rate (S=0.878)

图14 为围压800 kPa 时不同球度的相对破碎率曲线。由图可知，相对破碎率关系曲线同样可分为二阶段描述，当轴向应变小于界限应变时，颗粒几乎不破碎，相对破碎率近似为0，当轴向应变大于界限应变时，相对破碎率随应变增大而增大；当轴向应变相同时，相对破碎率随S减小而增大。以上计算结果说明，颗粒球度越小，形状越不规则，相同应力水平时越易破碎，这与2.2.1节分析结果是一致的。

图14 球度对相对破碎率的影响(800 kPa)Fig. 14 The effect of S on relative breakage rate (800 kPa)

图15 是不同围压和球度的相对破碎率。由图可知，S相同时，相对破碎率随围压增大而增加，近似线性关系；围压相同时，相对破碎率随S增大而减小，近似成半对数线性关系。分析表明，围压越大，S越小，颗粒破碎越多，相对破碎率越大，这与邹德高等[6]的试验结果是一致的，也与2.1.2 节峰值摩擦角计算结果相吻合。

图15 围压和球度对相对破碎率的影响Fig. 15 The effect of confining pressure and S on relative breakage rate

计算不破碎试样与破碎试样的峰值摩擦角差值，并将其与相对破碎率对应，如图16 所示。由图可知，二者近似线性关系，即颗粒破碎越显著，峰值摩擦角差值越大，且基本不受球度的影响。

图16 相对破碎率和峰值摩擦角差值的关系Fig. 16 Relation between relative breakage rate and difference of peak friction angle

3 结论

本文以爆破堆石料为研究对象，采用颗粒流软件PFC3D 模拟分析了颗粒球度对堆石料力学特性的影响及其细观机制，得到以下结论：

(1) 堆石料峰值摩擦角同时受颗粒强度、围压以及颗粒球度的影响。颗粒强度高、围压低时，峰值摩擦角随球度增大而减小；颗粒强度低、围压高时，峰值摩擦角随球度增大而增大。

(2) 颗粒球度对堆石料强度影响体现在颗粒咬合和颗粒破碎两方面。颗粒球度越小、形状越不规则，颗粒间咬合、内锁程度越高，增大试样强度；而颗粒球度越小，颗粒内拉应力水平越高、颗粒破碎越多，降低试样强度。颗粒球度对堆石料强度的影响是二者综合作用的结果。

(3) 堆石料初始弹性模量受颗粒破碎影响较小，随围压与增大而增大，随颗粒球度的增大而增大，但与颗粒是否破碎关系不大。