基于可破碎离散元法的堆石料应力变形及剪胀特性缩尺效应研究

徐 琨，杨启贵，周 伟，马 刚，黄泉水

（1.长江勘测规划设计研究院，武汉 430010；2.水利部土石坝破坏机理与防控技术重点实验室，南京 210029；3.武汉大学水资源与水电工程科学国家重点实验室，武汉 430072）

0 引 言

堆石料作为混凝土面板堆石坝的主要筑坝材料，具有较宽的颗粒级配，最小的颗粒粒径可小于0.1 mm，最大的颗粒粒径可达1 000 mm 以上，国内一些面板堆石坝工程中采用的粒径最大可达800～1 600 mm［1］。室内试验是当前获取堆石料力学变形参数的重要手段，《土工试验规程》（SL237-1999）中规定试样直径应不小于试验土料最大粒径的5～6 倍，由于筑坝堆石料最大粒径过大，受限于当前试验设备和技术水平，只能对缩尺后的堆石料开展力学试验。缩尺后堆石料与原级配堆石料所呈现出的力学变形特性差异，称之为堆石料的缩尺效应。

当前众多学者针对堆石料缩尺效应开展了诸多有益的研究。在力学特性的缩尺效应方面，文献［2-8］通过研究发现不同尺寸试样的抗剪强度一般呈现出随最大粒径的增大而降低，大尺寸试样破碎率高于小尺寸试样破碎率的规律；但也有文献［9-11］研究发现堆石料力学特性呈现了趋势相反的缩尺规律。

在变形特性缩尺效应方面，王继庄［12］和郦能惠等［9］通过研究均发现体积弹性模量随最大粒径的增大而增加，而李翀等［3］研究发现体积弹性模量随最大粒径的增大而减小。褚福永等［13］和武利强等［14］研究表明，堆石料初始弹模和体积切线模量随最大粒径的增加而增大，初始泊松比随最大粒径的增大逐渐减小。Wei等［15］研究认为粗粒料试样的压缩模量随最大粒径的增大而增加。花俊杰等［16］研究认为大尺寸试样的切线模量更高。王永明等［17］研究表明按相同相对密度制样时，最大粒径越大的试样其体积模量和初始弹模越大；按相同干密度制样时，试样的变形参数随最大粒径的增加呈现出先减小后增大的非单调关系。武利强等［7］研究表明，同一干密度制样时，粗粒料的变形模量随最大粒径的增大而减小。孔宪京等［18］基于超大型三轴仪研究了堆石料相变点处的体变缩尺规律，结果表明试样尺寸越大，相变点处的体变也越大，也即剪缩性越强。

以上研究表明，缩尺效应不仅对堆石料力学特性有影响，同时对其变形特性也有显著影响。当前，在堆石料缩尺规律方面取得了丰富的研究成果，在缩尺机理上有一定程度的认知，但由于堆石料缩尺规律的影响因素较多［4，7，15，19，20］，在多个影响因素的耦合作用下，堆石料力学变形特性的缩尺效应仍有待深入研究，缩尺引起的堆石料应力变形差异产生机理还有待明晰。在当前试验条件和技术水平下，室内试验一方面在试样尺寸上有很大的限制，另一方面难以监测到试样细观层面的实时演化情况，无法非常有效的揭示堆石料缩尺效应产生的机理。

离散元方法作为近些年来逐渐兴起的一种重要数值模拟研究手段，为从细观尺度研究堆石料力学响应及其机理提供了有效的途径。为进一步揭示堆石料缩尺效应产生机理，本文基于离散元方法开展堆石料缩尺效应研究。考虑剪胀特性是土体所具有的特殊性质，也是描述堆石料变形特性的关键指标之一，已有研究初步揭示了堆石料剪胀特性的缩尺规律，但缩尺产生的机理尚未有深入研究，仍不清晰。鉴于此，本文通过开展不同尺寸试样的数值三轴剪切试验，探讨堆石料力学变形缩尺规律，并从宏、细观层面深入分析试样剪胀特性的缩尺规律，揭示其产生的细观机理。

1 数值试验准备

1.1 颗粒破碎的模拟

研究采用PFC3D［21］模拟常规三轴剪切试验，模拟中使用“圆球”作为堆石料颗粒，从而消除颗粒形状对试验结果的影响。颗粒破碎是影响缩尺效应的主要原因之一［22，5］，研究中采用碎片替代法（Fragment Replacement Method，FRM）［23］模拟颗粒破碎。采用FRM 模拟颗粒破碎需要考虑两大问题：①颗粒破碎准则；②碎片替换模式。本研究中颗粒破碎准则采用八面体剪应力准则［24］，碎片替换模式采用Apollonian 堆积态碎片替换模式［25］，可以较为真实的反映颗粒材料的力学响应特性，具体细节可参见文献［8］，本文不再赘述。

1.2 数值试样级配曲线

如图1所示，本研究采用相似级配法对原级配缩尺，原级配堆石料最大颗粒粒径（dmax）为600 mm，分别缩尺得到dmax为60、90和120 mm 的三条级配曲线。图1中所示级配曲线由Tyler公式［26］［pi= 100(di/dmax)3-D］计算得到，其中，pi为小于粒径di的所有颗粒的累计质量所占所有颗粒总质量的百分比，D为描述级配曲线分形特征的维数。考虑模拟效率，试样中颗粒的数目不能过多，同时试验所选取的颗粒材料的级配宽度也不能太宽；经反复模拟测试，最终取D为2.0，颗粒最小粒径（dmin）为10 mm，对于粒径10 mm 以下的颗粒以等质量方式替换为粒径10 mm的颗粒。

1.3 数值试样制备

研究采用的试样径径比为5，在圆柱形空间内采用图1 中所示的3 种缩尺后级配，通过粒径膨胀法生成三轴试样。制得试样的尺寸为300×600、450×900 和600×1 200（直径×高度），单位均为mm，各尺寸试样中dmax分别为60、90和120 mm。

图1 原级配曲线及其缩尺后的级配曲线Fig.1 Prototype gradation curve and the scaled gradation curves

由于3 种尺寸试样的级配不同，为使不同试样初始密度状态达到基本相同的水平，制样方法和参数参考文献［27］。最终制得100 kPa 围压下3 种尺寸的初始试样（图2），初始试样的孔隙率分别为30.65%（dmax=60 mm）、28.23%（dmax=90 mm）和27.03%（dmax=120 mm）。结果表明，试样的最大粒径越大，采用相同初始密度状态控制制样得到的试样密度越大，与文献［19，28］研究结果一致。

图2 不同缩尺级配的初始试样示意图Fig.2 Illustration of the initial samples with different particle size distribution

研究中采用的接触模型及离散元参数参考文献［29，30，8］选取，如表1所示。

表1 离散元模拟参数Tab.1 Input parameters for DEM simulation

研究选取0.4，0.8 和1.2 MPa 三种试验围压，对各试样以最大颗粒粒径和试验围压为标识编号，如D60CP400 表示最大颗粒粒径为60 mm，试验围压为0.4 MPa（400 kPa）的试样，其中D代表“Diameter”，CP 代表“Confining Pressure”，其他试样编号均遵循该规则。

2 数值试验结果与分析

2.1 力学特性的缩尺规律

三轴数值剪切试验的应力-应变演化曲线如图3 所示。由图3 可知，不同尺寸试样间的应力-应变曲线存在显著差异，主要表现为尺寸越大试样的峰值强度越小，且围压越高时差异越明显。提取各试样峰值摩擦角φpeak［sinφpeak=(σ1-σ3)/(σ1+σ3)，其中σ1为大主应力，本文研究中也为轴向应力，σ3为小主应力，本文研究中也即围压］表征试样抗剪强度，得到峰值摩擦角与dmax的关系曲线如图4所示。可知，试样峰值摩擦角随试样最大粒径的增大而减小，也即试样尺寸越大，试样抗剪强度越低，且围压越大时越明显。以上规律与室内试验研究结果［31，32］一致，表明研究中开展的三轴数值剪切试验是合理有效的。

图3 应力-应变演化曲线Fig.3 Stress-strain curves

图4 峰值摩擦角与dmax的关系曲线Fig.4 Relationship between peak friction angle and dmax

根据De Mello［33］对粗粒料提出的幂函数非线性公式［τ=A(σn)b，A和b为拟合参数，τ为粗粒料的剪切强度，σn为正应力］来拟合得到研究中各组试样的剪切强度包络线，如图5 所示。可知，拟合得到的b值均小于1，试样剪切强度包络线为一条斜率逐渐变缓的曲线，即随正应力的增加试样剪切强度的增加幅度有减小的趋势，且最大粒径越大的试样，其增幅减小越明显。

图5 三轴试样剪切强度包络线Fig.5 Shear strength envelope of triaxial tests

2.2 变形特性的缩尺规律

三轴数值剪切试验的体变-应变演化曲线如图6所示，本研究中取体变剪缩为负，剪胀为正。由于生成的初始试样是较为密实的试样，加载过程中试样呈现出的剪缩量很小。由图6 可知，不同尺寸试样间的体变演化存在差异，尺寸越大的试样在演化过程中呈现出更弱的剪胀特性，且围压越高时这一差异越明显。

图6 不同围压下试样体变-应变曲线Fig.6 Volumetric strain response under different confining pressures

为更直观对比不同试验围压下各尺寸试样剪胀剪缩的差异，提取相变点处和10%剪应变处试样的体变量的绝对值如图7所示。其中，相变点［34］为剪胀性材料所特有的性质，即体积由体缩变为体胀的特征点，相变点处体变量也即试样的最大剪缩量。由图7（a）可知，10%剪应变处，围压越高的试样，剪胀量越小，且试样剪胀量随着最大粒径的增大而降低；围压越高时，随着最大粒径的增大，剪胀量降低程度越大，如：dmax由60 mm增大至120 mm，围压为0.4 MPa 时，剪胀量降低约15.8%，围压为1.2 MPa 时，剪胀量降低约66.2%。由图7（b）可知，相变点处，试验围压越高，试样的剪缩量越大；试样剪缩量随着最大粒径的增大而增加。以上结果表明，试样尺寸越大，剪缩特性越强。

图7 试样剪胀特性随最大粒径的变化关系Fig.7 Relationship of the shear dilatancy characteristics with dmax

2.3 堆石料剪胀缩尺效应产生的原因分析

现有研究指出颗粒破碎是导致堆石料力学响应缩尺效应的关键因素，试样尺寸越大，颗粒破碎率越高，抗剪强度越弱［5］，而导致这一现象的细观力学机理则是由于剪切过程中，大尺寸试样较小尺寸试样产生了更多的力学不稳定颗粒导致的［8］。

对于堆石料剪胀特性，颗粒破碎同样也是其缩尺规律的关键因素。颗粒材料的剪胀主要来源于颗粒错动、转动等重新排列后空隙的增加，而颗粒破碎会产生更为细小颗粒碎片，细小颗粒碎片可填充于大颗粒间的空隙中，从而减弱了试样的剪胀效应［35］；颗粒破碎率越高，则剪胀效应减弱地越多，从而造成了尺寸越大试样的剪缩特性越强。但受限于与室内试验监测手段的限制，这一分析目前尚缺少试验数据的支撑，细观机理未得到深入分析，为此借助数值试验的优势开展深入研究。

（1）有效配位数演化规律。颗粒配位数是散粒体结构特征的重要表征参数之一，表示与目标颗粒发生接触的颗粒数目［36，37］。Thornton［38］在研究中将配位数小于2 的颗粒为称之为悬浮颗粒，并将试样中颗粒配位数小于2 的悬浮颗粒除去后得到的试样平均配位数定义为有效配位数CNeff，可表示为：

式中：C为试样中各颗粒间接触的总数；Ntot、N0和N1分别为试样中颗粒的总数、配位数为0 的颗粒数目和配位数为1 的颗粒数目。

缩尺试样有效配位数演化情况如图8 所示。由图8 可知，不同围压下，各尺寸试样的有效配位数均呈现出相似的演化规律，表现为：有效配位数随着试样剪切先呈现出陡降，之后缓慢减小，随着剪切继续，有效配位数逐步达到基本稳定状态，这与Duran 等［39］和Gu 等［40］的研究结果相似。同时，也发现各尺寸试样有效配位数在数值上存在较为明显的差异，表现为：同一剪切位移时，有效配位数随试样最大颗粒粒径的增大而降低；结果表明，剪切试验中，试样尺寸越大，其承力结构越稀疏。

图8 试样有效配位数演化曲线Fig.8 Evolution curves of effective coordination number

（2）各粒径组有效配位数。有效配位数的演化结果展示了不同尺寸试样间承力结构的整体性差异，但不足以解释堆石料试样剪胀特性差异产生的原因。为深入分析，将试样中颗粒按粒径分组，然后求出各粒径组内颗粒的有效配位数，得到各粒径组有效配位数分布。分组规则为，粒径小于等于10 mm 的为一组，10 到20 mm 的为第二组，依次类推；所有分组均为左开右闭区间，颗粒最大粒径为60、90 和120 mm 的试样分别分出6、9和12组。

初始试样的各粒径组有效配位数直方图如图9 所示。可知，不同尺寸试样中各粒径组有效配位数均呈现随颗粒粒径的增大而升高的趋势，但最大粒径组的有效配位数处于同一量级，相差不大。对于同一粒径组，dmax越大的试样，其粒径组有效配位数越低，如（50，60］粒径组时，颗粒最大粒径为60、90 和120 mm 的试样其粒径组有效配位数分别为25.09、11.37和6.48。

将粒径组和有效配位数归一化处理，得到初始试样各粒径组有效配位数在半对数坐标下的分布，如图10 所示。由图可知，虽然标准化后的各粒径组有效配位数分布有一些差异，但大体呈现出较为相似的分布，表明采用相似级配法缩尺得到的各尺寸试样在一定程度上保持初始承力结构的相似性。

图10 归一化的初始试样各粒径组有效配位数分布Fig.10 Normalized effective coordination number distribution of each particle size group of the initial sample

图11所示为剪切前后各粒径组有效配位数分布对比情况，各试验围压下分布规律相似，仅列出0.4 和1.2MPa 围压下的结果。由图11可知，不同尺寸试样的各粒径组有效配位数分布在剪切后（ε1= 10%）呈现出不同的演化趋势。剪切后，小尺寸试样（dmax= 60 mm）各粒径组有效配位数均有所降低，整体降低明显；随试样尺寸的增大，其各粒径组有效配位数降低程度减弱，并有部分粒径组下的有效配位数高于剪切前的值。此外，不同试验围压下，试样剪切前后的各粒径组有效配位数变化情况有所差异，表现出围压越大且试样尺寸越大时，剪切后的各粒径组有效配位数较剪切前的上抬越明显。

图11 试样各粒径组有效配位数分布Fig.11 Effective coordination number distribution of each particle size group

另开展一组不考虑颗粒破碎的三轴剪切试样做对比研究，得到试样剪切前后各粒径组有效配位数分布如图12所示，仅展示1.2 MPa 试验围压。对比图12 与图11（b）可知，颗粒破碎对试样剪切后各粒径组有效配位数分布产生了显著影响，颗粒破碎使剪切后的各粒径组有效配位分布发生上抬，且试样尺寸越大时上抬程度越大。采用Marsal［2］对试样颗粒破碎率的定义，提取试样颗粒破碎率如图13 所示。由图可知，试样尺寸越大，颗粒破碎率越大，且试验围压越大时，颗粒破碎率越大。综上，可以证明各粒径组有效配位数分布的上抬与颗粒破碎率是成正相关的，即试样颗粒破碎率越大，各粒径组有效配位数分布上抬越多。

图12 试样各粒径组有效配位数分布（颗粒不破碎）Fig.12 Effective coordination number distribution of each particle size group（uncrushable particle）

图13 颗粒破碎率BmFig.13 Particle breakage factor Bm

有效配位数在一定程度上表征试样的密实程度，配位数越高表明颗粒接触越多，试样越密实。当试样颗粒发生破碎后，碎片填充于颗粒间的空隙中，随着剪切的进行逐渐充满密实，试样的承力结构变得密实，也即有效配位数增大；颗粒破碎越多时，则试样填充越密实，有效配位数增大越多，由此揭示了不同尺寸试样各粒径组有效配位数分布剪切前后演化差异产生的原因。

此外，有效配位数的减小在一定程度上代表试样的剪胀，反之代表试样的剪缩。为更清晰展现剪切前后试样各粒径组有效配位数的变化情况，定义各粒径组有效配位数分布差ΔCNeff=CNeff|ε1=10%-CNeff|ε1=0%。其中，CNeff|ε1=0%、CNeff|ε1=10%分别代表剪切前、后的各粒径组有效配位数分布，ΔCNeff为正值代表剪切后发生了剪缩，负值代表剪切后发生了剪胀。对比颗粒不破碎（图12）和颗粒破碎［图11（b）］两种情况，得到各粒径组有效配位数分布差如图14 所示。可知，不同尺寸试样间ΔCNeff的有明显差异，同时考虑颗粒破碎的试样较不考虑颗粒破碎的试样的ΔCNeff均在较大粒径区间有明显增大，试样尺寸越大则ΔCNeff增大越明显，也即试样发生了更多程度的剪缩。同时，统计得到各试验围压下不考虑颗粒破碎试样剪切后（ε1= 10%）的体变均在+5.00±0.1%范围内，表明不考虑颗粒破碎时各尺寸试样的剪胀量基本相同。综上可知，一方面，颗粒破碎促使ΔCNeff向增大的方向演化；另一方面，尺寸越大试样的大粒径区间的ΔCNeff增多越明显，降低了试样的剪胀量，使得试样尺寸越大，剪胀性越弱；由此揭示了试样剪胀缩尺规律发生的原因。

图14 各粒径组有效配位数分布差Fig.14 Difference of effective coordination number distribution of each particle size group

3 结 论

本文基于可破碎离散元法，从宏、细观层面开展了堆石料应力变形及剪胀特性缩尺效应研究，得到了如下主要结论：

（1）数值模拟试验较好反映了堆石料应力变形的缩尺规律。试样尺寸越大，抗剪强度越低，剪缩特性越强，最终的剪缩量越大；试验围压越高，缩尺效应对堆石料抗剪强度、剪胀特性的影响越大，缩尺现象越显著。

（2）采用相似级配法缩尺得到的不同尺寸试样一定程度上保持着初始承力结构的相似性，但尺寸越大的试样，其承力结构越稀疏。试样尺寸越大，各粒径组有效配位数分布的降低程度减弱，并有部分粒径组下的有效配位数高于剪切前的值，围压越大且试样尺寸越大时，这一现象越明显。

（4）试样的各粒径组有效配位数分布的演化与试样颗粒破碎率成正相关，试样颗粒破碎率越大，各粒径组有效配位数分布上抬越多；通过对比研究，揭示了不同尺寸试样间各粒径组有效配位数分布差（ΔCNeff）演化的差异是试样剪胀特性缩尺规律产生的细观机理之一。 □