土石坝因其建造速度快、对基础适应能力强等优点已经成为应用最广的坝型之一。目前,在建或者拟建的200～300 m级土石坝多位于地震频发的区域，其抗震安全问题至关重要。堆石料作为土石坝的主要支撑材料其动力特性研究得到了广泛关注。

一般来说，堆石料的最大动模量随着围压和固结应力比的增加而增大^[1-3]。相同循环荷载作用下，母岩的岩性越差，颗粒破碎率越大^[4]。低围压下堆石料的动力特性受级配优劣的影响并不显著，但随着围压增加，级配优的试样动力特性指标更优。对浑圆度较好的砂砾料，级配优劣对颗粒破碎率的影响较小，相应地，对动力变形特性的影响逐渐减弱^[5]。施加动载前的颗粒破碎情况会影响后续动荷载作用过程中的颗粒破碎量，且颗粒破碎是引起试样体积收缩的主要因素^[6]。由此可知颗粒破碎是影响堆石料动力特性的重要方面。

离散单元法是从细观角度分析岩土体静动力特性的有效手段^[7-9]。刘汉龙等^[10]采用2维离散元方法对筑坝反滤料低围压下的模量阻尼特性进行模拟，提出Hertz–Mindlin非线性接触模型可以较好地反映材料的非线性特征。杨贵等^[11]根据中小尺寸堆石料动模量变化规律采用离散元模拟预测大尺寸试样的模量值。通过引入破碎准则研究了循环荷载作用下围压和固结应力比对颗粒破碎的影响^[12]。但是目前进行堆石料动力特性模拟时均采用圆盘或圆球模拟堆石颗粒，未考虑颗粒形状的不规则性且考虑颗粒破碎对动模量影响的数值模拟研究较少。

以糯扎渡Ⅱ区粗堆石料为研究对象，采用多个小球组合模拟不规则的颗粒形状，引入颗粒破碎准则，通过拟合不同围压下动三轴骨架曲线确定一套细观参数。在此基础上研究了颗粒破碎和孔隙率对堆石料动弹性模量的影响规律。

堆石颗粒采用多个直径相等的小球组合而成。首先，生成长宽比范围为1.2～1.8之间的椭球体，在椭球面上随机选取50个点作为多面体的顶点^[13]。随后，按最密六方排列生成一系列小球，小球完全覆盖多面体的区域。最后，通过判断小球与多面体的位置关系，删掉处于多面体外部的小球，剩下的小球组合得到不规则形状的颗粒模板。部分颗粒模板如图1（a）所示。

图1(Fig. 1)

图1 数值模拟的不规则颗粒与级配曲线 Fig. 1 Typical particles with irregular shape in simulations and grain size distribution

试验材料为糯扎渡Ⅱ区粗堆石料，岩性为花岗岩^[14]。为保证数值模拟的可行性和计算效率，研究者多采用截断级配的方式^[15-16]。张宜等^[17]指出堆石体的强度和变形特性随截断粒径的增大而单调变化，且截断粒径取15 mm时其影响可以接受。截断后的模拟级配与室内试验级配如图1（b）所示。采用PFC内置的Hertz接触模型来模拟堆石料的非线性特征^[10]。制样过程中通过给定不同的摩擦系数可以得到不同孔隙率的试样^[18]。模拟时不考虑墙体与颗粒之间的摩擦，制样过程中不考虑破碎且颗粒间的摩擦系数设为0.15。制样时,首先在给定空间范围内生成满足级配要求的球形颗粒。然后将球形颗粒替换为不规则颗粒，替换完成后缓慢移动墙体将试样压缩至30 kPa围压下固结稳定。最后将摩擦系数设为其真实值。生成的数值试样尺寸为300 mm×300 mm×300 mm，由3 519个颗粒共84 310个小球组成，初始孔隙率为0.55。该孔隙率大于室内试样的孔隙率主要是因为模拟时对级配进行截断处理导致没有足够的小颗粒填充孔隙^[19]。Huang等^[20]模拟循环荷载下砂土的力学特性时设置局部阻尼系数为0.1，并通过敏感性分析发现振动角频率ω小于10π rad/s时模拟结果几乎不受影响。局部阻尼系数和振动角频率均参考其取值。

动三轴试验围压为900、1500和2 200 kPa，固结比K_c=σ₁/σ₃为1.5，加载工况如表1所示。数值模拟时先将试样固结到指定围压，稳定后施加偏应力直到试样的固结比为1.5。动力试验采用应力控制，加载波型为正弦波。一个试样施加一个动应力σ_d，每一级动应力条件下振动3圈，取第3圈计算动弹性模量。不排水循环加载条件的模拟与刘汉龙等学者的处理方法一致^{[10, 20-21]}。循环加载过程中固定上侧墙体，每循环一步均对墙体的速度进行调整。先根据当前时间和时步大小计算下一循环步结束对应的时间t，计算得到动应力的目标值为σ_dsin(ωt)。通过伺服机制指定底部墙体的运动速度使试样动应力满足要求。再根据底部墙体的速度计算周围四面墙体的运动方向与速度，确保整个振动过程中试样的体积不变。

表1(Tab. 1)

表1 模拟的试验方案 Tab. 1 Test schemes for dynamic simulation

表1 模拟的试验方案 Tab. 1 Test schemes for dynamic simulation 围压/kPa 固结比 动应力/kPa 900 1.5 200, 400, 600, 900 1 500 1.5 200, 400, 600, 900, 1200, 1400 2 200 1.5 200, 400, 600, 900, 1200, 1400, 1800, 2000

当颗粒受到的最大接触力或等效应力达到极限状态时颗粒会发生破碎^[22-23]。Russell等^[24]提出遭受多个外荷载的颗粒，其破碎主要受最大接触力控制。在其工作的基础上，Ciantia等^[25]提出一个通用高效的破坏准则。如图2所示，当最大法向接触力满足如下关系式时，颗粒发生破碎。

图2(Fig. 2)

图2 接触力作用于球体表面 Fig. 2 Contact force acts to sphere surface

式中：F_max、d分别为作用在颗粒上的最大接触力和颗粒直径；θ₀决定了接触面积的大小，根据赫兹接触理论确定^[25]；σ_lim为极限强度，与材料的泊松比、单轴抗压和抗拉强度相关。

极限强度σ_lim假设符合Weibull分布，且极限强度是粒径相关的：

式中：m为Weibull模量，表征颗粒强度分布的离散程度；d₀为参考粒径，本文选为38 mm；σ_lim,0为粒径d₀的颗粒残存概率37%对应的极限强度；η为颗粒极限强度分布的尺寸效应参数。

图3为颗粒在加载过程中可能发生的破坏过程示意图。

图3(Fig. 3)

图3 颗粒的可能破碎过程示意图 Fig. 3 Possible crushing process of a particle

对任意颗粒，当其受力满足公式（1）时颗粒发生破碎，该颗粒沿最大接触力方向劈裂为2个碎片。当碎片颗粒再次达到破碎条件时劈裂为2个新的碎片，直到碎片颗粒仅由1个基本球组成时停止。新生成的碎片颗粒属性与原颗粒一致，强度根据碎片颗粒的等效直径重新指定^[25-26]。

细观模拟需要确定的参数主要分为接触模型相关参数和颗粒强度参数两类。接触模型参数为剪切模量G、泊松比v、摩擦系数μ。颗粒强度参数包括Weibull模量m、尺寸效应参数η和极限强度σ_lim,0。细观参数的确定流程是在参考前人细观参数选取规则和成果的基础上，确定细观参数的大致范围；然后，对不同的参数进行敏感性分析，了解各个参数对动力特性的影响；最后，通过微调得到一套可以反映研究对象变形特性的参数。根据杨贵等^[11]进行的参数敏感性分析可知，细观泊松比v在小于0.3的情况下，对试样的动模量影响不大，参考其取值，细观泊松比v取0.3。文献[27]确定花岗岩颗粒强度相关的分布参数m和η。这样需要确定的细观参数为剪切模量G、摩擦系数μ和极限强度σ_lim,0。

首先，对剪切模量G、摩擦系数μ和极限强度σ_lim,0进行敏感性分析，3个参数对骨架曲线的影响如图4所示。

图4(Fig. 4)

图4 参数敏感性分析 Fig. 4 Parameter sensitivity analysis

剪切模量直接影响试样的变形能力。由图4（a）可知：剪切模量越大，相同动应变对应的动应力越大；但是随着剪切模量增加，骨架曲线随动应变的增长趋势变缓。图4（b）为摩擦系数对骨架曲线的影响。由图4（b）可知：摩擦系数越小，颗粒之间越容易发生相对滑动，在相同荷载作用的条件下试样产生的变形越大；动应变相同时，动应力随着摩擦系数的增加而增大。图4（c）为极限强度对骨架曲线的影响。由图4（c）可知：动应变较小时，不同极限强度对应的动应力差异不大；随着动应变增加，极限强度越大的试样其动应力也越大。

根据参数敏感性分析结果，通过与2 200 kPa围压下的骨架曲线对比，可以大致确定细观参数的取值范围。通过进一步试算，并与2 200 kPa围压下的骨架曲线对比，可以得到一套合理的细观参数，如表2所示。

表2(Tab. 2)

表2 细观模拟参数 Tab. 2 Micro-parameters for simulation

表2 细观模拟参数 Tab. 2 Micro-parameters for simulation 模型参数 值 模型参数 值 密度/(kg·m–3) 2 630 d0/mm 38 G/GPa 22 m 3.41 v 0.3 η 0.42 μ 0.58 σlim,0/MPa 1 600

为进一步验证参数的合理性，采用表2中的细观参数，模拟900和1 500 kPa围压下的骨架曲线，如图5所示。模拟的动应力–动应变骨架曲线与室内试验结果具有较好的一致性，说明该模型可以基本再现糯扎渡Ⅱ区粗堆石料的动力变形特性。

图5(Fig. 5)

图5 模拟和试验骨架曲线对比 Fig. 5 Comparison of simulation and experimental skeleton curves

图6为考虑颗粒破碎和不考虑颗粒破碎条件下，动应变随循环振次的变化规律。由图6可知：试样在循环振动过程中，均发生了比较明显的塑性变形，这主要是因为初始偏应力的存在使得动应变曲线偏向压缩侧，累积塑性变形均为压应变，同一围压下，动应变幅值和塑性变形均随动应力的增加而增大；动荷载相同时，围压大的试样对应的动应变幅值小；相同应力状态下，考虑颗粒破碎的试样动应变幅值较不考虑颗粒破碎的试样大，相应的塑性变形也越大，说明颗粒破碎导致试样整体刚度降低，使试样产生更多不可恢复的变形。

图6(Fig. 6)

图6 动应变随循环振次的变化关系 Fig. 6 Curves of dynamic strain–cycle number

颗粒破碎程度采用Marsal破碎率衡量。整理振动过程中不同工况下试样的颗粒破碎率，如图7所示。由图7可知：围压越大，固结阶段颗粒破碎越严重；颗粒破碎率在初始加载的1/2周期内显著增加，随后增长减缓，说明循环荷载作用下颗粒破碎主要发生在第1圈的前1/2周期内；围压相同时，动应力越大，颗粒破碎率也越大。

图7(Fig. 7)

图7 循环加载过程中的颗粒破碎率 Fig. 7 Particle crushing during cyclic loading process

对任意颗粒，配位数为与该颗粒接触的周围颗粒数量。Thornton^[28]指出，配位数小于2的悬浮颗粒对集合体的力学特性贡献较小，可以忽略。试样的有效配位数定义为：

式中， ${N_{{\rm{total}}}} $ 、 ${N_{{\rm{eff}}}} $ 、N₀和N₁分别为总的颗粒数量、去除悬浮颗粒后的有效颗粒数、配位数为0和配位数为1的颗粒数，C为颗粒接触总数。

图8（a）为循环振动过程中不同动应力条件下，有效配位数的变化规律。由图8（a）可知：不考虑颗粒破碎条件下，动应力相同时围压越大，试样的有效配位数越多；围压相同时，动应力越大，试样的有效配位数越小；循环过程中有效配位数会缓慢减小，且颗粒破碎也会导致有效配位数的减少；配位数的频率为具有相同配位数的颗粒数量与有效颗粒数之比。整理围压2 200 kPa、动应力为2 000 kPa时，振动前后的配位数的频率分布，如图8（b）所示，图8（b）中振动前指试样固结完成施加动荷载前的状态。由图8（b）可知，考虑破碎和不考虑破碎条件下配位数的频率分布图相似；振动前，考虑破碎和不考虑破碎，试样的峰值配位数均为5.0；振动后，位于峰值配位数左侧的频率增加，右侧降低，且发生颗粒破碎的试样峰值配位数降低到4.0。说明振动和颗粒破碎均会导致试样中力学不稳定颗粒的增加，使试样内部的结构性变差，且颗粒破碎量越大，对试样的结构性影响越显著，最终导致试样的刚度和承载能力降低。

图8(Fig. 8)

图8 有效配位数的演化曲线与配位数频率分布 Fig. 8 Evolution curves of effective coordination number and frequency distribution of coordination number of samples before and after cyclic loading

为了研究孔隙率对动力特性的影响，将制样时颗粒的摩擦系数设置为0，采用同样的制样方法制得最密实的试样，孔隙率为0.52。施加与前文孔隙率为0.55的试样相同的振动工况，记录振动过程中的颗粒破碎率和配位数等参数。

图9为孔隙率对动应变、颗粒破碎率和配位数分布的影响。

图9(Fig. 9)

图9 孔隙率对动应变、颗粒破碎和配位数分布的影响 Fig. 9 Influence of porosity on dynamic strain, particle crushing and the distribution of coordination number

由图9（a）可知，孔隙率对动力变形特性影响显著。在同样的围压和动应力条件下，孔隙率大的试样具有更大的动应变。振动过程中的颗粒破碎发展如图9（b）所示，孔隙率小的试样在相同应力状态下，颗粒破碎率较孔隙率大的试样小。

整理动应力为2 000 kPa振动前后配位数的频率分布图如图9（c）所示。由图9（c）可知：孔隙率为0.55和0.52的试样频率峰值对应的配位数均为5；孔隙率为0.55的试样峰值左侧的频率均高于孔隙率为0.52的试样，峰值右侧频率均比孔隙率为0.52的试样低；振动后频率分布曲线的峰值有所下降，峰值左侧上抬，右侧下降，且孔隙率为0.55的试样变化趋势更显著。说明同一级配的同种材料，密实试样内部的平均接触数更多，在相同受力状态下不易发生颗粒破碎，力学性能更好。

同一围压下的动应变与动弹性模量可以近似用直线关系式（4）拟合^[29]：

式中，a和b分别为拟合直线的截距和斜率，1/a为最大动弹性模量E_dmax。斜率b越大，表明动弹性模量随着动应变的增加衰减越快。

最大动弹性模量与试样的应力状态相关，可表示为：

式中：σ_m为平均有效固结主应力，σ_m=（σ₁+2σ₃）/3；P_a为大气压力；k和n为试验常数。

整理不同围压下的动弹性模量与动应变值，绘制动弹性模量的倒数与动应变的关系，采用式（4）拟合。图10汇总了不同情况下，斜率b的分布规律。由图10可知：斜率b均随着平均有效固结主应力的增加而减小，说明动弹性模量的衰减速率随有效固结主应力的增加变缓；平均有效固结主应力相同时，考虑颗粒破碎的试样较不考虑颗粒破碎的试样具有更大的b值，动弹性模量衰减速率增加。不同孔隙率试样对应的斜率如图10（b）所示。由图10（b）可知，相同围压条件下小孔隙试样的斜率均小于大孔隙试样，说明疏松试样的动弹性模量随动应变的增加衰减更快，这与不同密实程度卵石的动模量衰减规律一致^[30]。

图10(Fig. 10)

图10 颗粒破碎和孔隙率对斜率b的影响 Fig. 10 Influence of particle breakage and porosity on the value of b

不同孔隙率试样的最大动弹性模量与平均有效固结主应力之间的关系如图11（a）所示。小孔隙率试样的最大动弹性模量比大孔隙率试样大。平均有效固结主应力相同时，孔隙率为0.52的试样最大动弹性模量约为孔隙率0.55试样的1.2倍。采用式（5）拟合，k值对应于平均有效围压为P_a时的最大动弹性模量。孔隙率为0.55的试样对应的k值为（5 615），明显小于孔隙率为0.52试样的k值（7 028）。n值的差异不显著。

图11(Fig. 11)

图11 孔隙率对最大动弹性模量的影响 Fig. 11 Influence of porosity on the maximum dynamic elastic modulus

Hardin等^[31-32]认为材料的变形特性主要与孔隙比和应力状态相关，提出可以用经验公式（6）来描述变形模量与孔隙比和平均有效应力之间的关系：

式中， $A $ 为材料常数，e为孔隙比，f(e)为孔隙比的函数。不同的研究者针对不同的岩土材料提出不同的表达形式用来消除孔隙率的影响，常用的函数形式为 $f\left( e \right) = {{{{\left( {2.973 - e} \right)}^2}}/ {\left( {1 + e} \right)}}$ 。

采用函数f（e）对最大动弹性模量进行归一化，如图11（b）所示。由图11（b）可知，不同孔隙不同围压下的最大动弹性模量可以采用统一的公式描述： ${E_{{\rm{dmax}}}}= $ $ 4\;085f\left( e \right){P_{\rm{a}}}{\left( {{{{\sigma _{\rm{m}}}}/ {{P_{\rm{a}}}}}} \right)^{{\rm{0}}{\rm{.34}}}}$ 。

整体而言，颗粒破碎主要影响最大动弹性模量的衰减速率。孔隙率小的试样受力性能更好，同样应力条件下颗粒不易发生破碎，模量随动应变的衰减速率也相应较慢。这对土石坝的抗震具有较好的借鉴意义，施工中可通过提高填筑密度来改善坝体的抗震性能。

采用PFC3D离散元程序对糯扎渡Ⅱ区粗堆石料循环荷载作用下的变形特性进行模拟。研究了小应变条件下，颗粒破碎和孔隙率对动弹性模量的影响。主要结论如下：

1）不同围压下，数值模拟的骨架曲线与室内堆石料的试验结果吻合较好。围压越大，相同动应变条件下试样的动弹性模量越大。

2）循环振动和颗粒破碎均会降低集合体内部的有效配位数。

3）孔隙率小的试样动力特性要优于孔隙率大的试样。相同围压和动应力条件下，孔隙率小的试样有效配位数大，颗粒破碎量小。

4）颗粒破碎和孔隙率均影响动弹性模量的衰减速率。最大动弹性模量主要受平均有效应力和孔隙率的影响，可以用Hardin等提出的经验公式描述孔隙率对最大动弹性模量的影响。