0.   引言

随着西部开发的拓展，我国的工程活动重心偏向西北、西南等地，而这些地区的地形地貌复杂，严峻的工程地质条件给工程建设带来了极大困难，同时也导致了很多灾难，如滑坡和泥石流等。这些区域具有代表性的工程介质就是碎石土，其物质组成的复杂性和结构分布的不规则给工程和科研带来极大挑战。有很多学者通过依托工程或现场及室内物理实验对其力学特性进行了相应的研究，也获得了一定的研究成果^[1~8]。李翀等^[9]和阳云华等^[10]分别对粗粒料和膨胀土进行了现场和室内试验，发现存在尺寸效应的影响。李翀等^[9]分析了试样直径以及最大粒径的影响，得知主应力差、内摩擦角和初始切线模量随试样直径的减小而增大。但由于只做了4种工况的试验，考虑的因素不全，不能充分反应力学性质。阳云华等^[10]提出室内试验抗剪强度取值的折减系数，具有一定参考价值，但不能通用。杨圣奇等^[11]、吕兆兴等^[12]发现不同长度的岩石，其长度对峰值应力前的变形特性没有显著影响，但明显改变峰后的变形特性，长度越大的岩样，岩石峰后越脆；并指出了岩石材料强度长度效应是由于岩样端部摩擦效应所致，而并非根源于材料的非均质性。吕兆兴等^[12]则通过非均质系数衡量尺寸效应。物理实验虽然直观，但实验过程繁琐复杂、耗时长且受尺寸限制，工作量大。为此，一部分学者借助于软件从数值模拟角度进行研究，但对于碎石土这种高度离散物质要实现模拟也存在一定困难，尤其是碎石土中块体的随机生成难以实现。肖昭然等^[13]利用PFC2D将球形单元粘结在一起模拟粗集料，虽然能够实现随机生成和多边形化，但是从二维的角度模拟三维实体存在一定差距。李耀旭^[14]则通过数字图像直接转换成PFC2D数据格式，生成土石混合体模型，这种方法在一定程度上保留了试样土体的真实性，但也局限于二维模拟。王新^[15]运用PFC3D颗粒流软件模拟土石混合体的三轴试验，从二维突破到三维空间，但是其对块体的模拟直接用不同粒径大单个球体代替，模拟过于粗糙，难以真实反映块体的特性，对模拟结果有一定的影响。

针对碎石土三轴试验中碎石的三维仿真随机生成难点，本文基于离散单元理论与PFC程序平台，考虑碎石的粒度、形状对物理实验的影响，试图建立几种典型碎石外形的生成法则，并以室内三轴试验标定细观参数为基础验证模型的可靠性。此外，从碎石土试样的直径和高度两个方面重点研究碎石土应力-应变曲线特性、峰值应力变化、峰值体应变、粘聚力及内摩擦角的尺寸效应。研究成果有助于提高三维离散仿真精度，更深入理解试件尺寸效应，便于获得可靠的残坡积碎石土物理力学指标参数，为相关设计和施工提供科学依据。

1.   材料与方法

1.1   实验背景

实验材料取自湖南昭山某具有级配代表性的堆积碎石土斜坡，天然状态含石为83%，含黏粒量17%，含水率11.4%；母岩为粉质砂岩，黏粒为红黏土；最大碎石粒径84 mm。碎石土及配额曲线如图 1。为了获取碎石土真实的物理力学参数(弹性模量、内摩擦角和粘聚力)，将取回的土样在不同含石量(20%，40%，60%)和不同含水率(7%，9%，11%，13%)下进行常规三轴试验，试验装置和结果见图 2。

图  1  碎石土颗粒组成

Figure  1.  Granular composition of gravel soil

下载: 全尺寸图片 幻灯片

图  2  碎石土试验结果

Figure  2.  Test results of gravel soil

下载: 全尺寸图片 幻灯片

1.2   仿真模型生成

基于三轴试验结果，以60%含石量、7%含水率一组为研究对象进行数值模拟。为了更加接近真实物理实验并减小模拟误差，数值模型尺寸与粒径大小应尽可能接近真实尺寸。因此，模型尺寸设置为直径101 mm、高度200 mm的圆柱体，与实际三轴试验试样尺寸一致。模拟过程中颗粒分为两种(土和碎石)，模型中土颗粒以及各粒径碎石的含量严格遵守三轴试样配置的百分比。在程序运行效率与颗粒生成总数之间折中考虑，模拟实验中小于2 mm粒径的颗粒统一用2.5 mm代替，2~5 mm的则用2.5~5 mm代替，这样既克服了计算机容量和速度的限制，又满足了计算精度的要求。数值模拟模型的具体级配组成见表 1。

表  1  模拟级配与试样级配对比

Table  1.  Comparison of simulated grading and sample grading

三轴试验试样尺寸和级配				数值模拟模型尺寸和级配
试样尺寸/mm	粒径/mm	含量/%		模型尺寸/mm	粒径/mm	含量/%
101×200	20~10	21.72		101×200	20~10	21.72
	10~5	20.19			10~5	20.19
	5~2	18.09			5.0~2.5	18.09
	＜2	40			2.5	40

下载: 导出CSV 

| 显示表格

碎石块的形状对模拟结果起着至关重要的作用。而实际碎石土中碎石块的形状大小极不规则，对碎石块进行模拟有一定困难，因此很多研究者在碎石土或土石混合体数值模拟中直接用大球代替碎石块，虽然能够获得应力-应变关系，但与实际情况有一定差距，仿真效果也不佳。本文在原有的Fish函数的基础上进行二次开发，根据母岩特性和三轴剪切后块体的状态特征，利用clump命令进行碎石块模拟。观测碎石土中块体的组成大致可分成3种(见图 3a)，因此分3种情况对碎石块进行模拟。3种块体的生成法则：① 按碎石百分含量在墙内(wall)随机生成碎石块体信息球，记录每一个信息球的三维坐标及半径；② 删除信息球，通过遍历记录信息球确定组成块体小球的坐标及半径，按一定生成规律(见图 3b)重新生成组成块体的小球，利用clump命令将小球合成块体，标识id；③ 计算块体体积，当生成的块体体积达到粒径组所占的百分比时，结束块的生成。生成效果见图 3c。

图  3  碎石块生成法则

Figure  3.  Law of crushed stone formation

下载: 全尺寸图片 幻灯片

块体生成后，根据试验实际孔隙率进行土颗粒的生成。利用fish函数，以试验测得的孔隙率为基准，生成含量为40%的土颗粒。为了使土颗粒生成后达到密实状态，先将所有土颗粒的尺寸缩小到设定尺寸的一半，然后通过孔隙率计算出半径扩大系数，最后进行半径扩张。仿真模型生成结果如图 4。

图  4  仿真模型生成

Figure  4.  Simulation model generation

下载: 全尺寸图片 幻灯片

1.3   宏微观参数的确定

在颗粒流模拟中，选择合适的微观参数是相当重要的，直接影响模拟结果的正确性。本文以三轴试验的应力-应变曲线为基础，对仿真模型进行微观参数的设置。通过反复调整模型中各参数，使得模拟应力-应变曲线最大程度接近试验应力-应变曲线。采用200 kPa围压的应力-应变曲线进行标定，再以100 kPa和400 kPa时的应力-应变关系进行验证(见图 5)。反复调整后最终确定的标定参数见表 2。

图  5  参数标定及验证

Figure  5.  Parameter calibration and verification

下载: 全尺寸图片 幻灯片

表  2  参数标定值

Table  2.  The calibration parameters

参数含义	程序代码	参数取值
杨氏模量	md_Ec	4.5e6
平行粘结模量	pb_Ec	1e7
密度	md_dens	2600
土颗粒法向刚度	kn	7.2e4.
块体法向刚度	kn	4e5
法向刚度与切向刚度比	kn/ks	0.8
接触粘结法向力平均值	cb_sn_mean	1.35e5
接触粘结切向力平均值	cb_ss_mean	1.35e5
平行粘结法向力平均值	pb_sn_mean	7e4
平行粘结切向力平均值	pb_ss_mean	7e4

下载: 导出CSV 

| 显示表格

图 5显示，围压200 kPa时模拟曲线与试验曲线吻合很好；100 kPa时应力-应变曲线稍有偏差，但整个曲线走势基本一致，最大值偏差为6.95%，400 kPa时虽然曲线后段吻合不是很好，但是两条曲线整体基本一致，最大值的误差为8.9%。初步分析出现这种现象的原因主要是颗粒的粒径以及模型试样内部的块体形状与分布的影响。从图 5可知，该模型的微观参数设置基本正确。

为了进一步验证模型参数的正确性以及分析模拟实验与三轴试验的误差，计算出模拟实验和三轴试验的粘聚力c与内摩擦角φ值(见表 3)，从表中可知模拟的c、φ值与试验的c、φ值非常接近，误差分别为14.3%和7.45%，每个围压下对应峰值应力偏差也不大(最大为9.18%)。据此认为模型参数较为合理，块体生成方法能较好地实现对碎石土的数值模拟。

表  3  模拟实验与三轴试验误差分析

Table  3.  Error analysis in simulation experiment and triaxial test

实验类型	峰值偏应力/kPa			粘聚力c/kPa	内摩擦角φ/(°)
	围压100 kPa	围压200 kPa	围压400 kPa
三轴试验	501.80	677.90	1119.20	79.89	30.67
模拟实验	466.91	708.41	1016.43	91.33	28.38
误差/%	6.95	4.50	9.18	14.30	7.45

下载: 导出CSV 

| 显示表格

2.   结果与讨论

碎石土中碎石块的粒径很不均匀，最大粒径与最小粒径相差极大，且受物理实验条件限制。为了进一步探究尺寸效应对碎石土物理力学特性的影响，以室内直径101 mm、高度200 mm圆柱体试样的物理实验为基础，借助离散元软件平台，建立不同高度和直径的数值模型对碎石土进行尺寸效应的对比研究，具体工况见表 4。

表  4  数值模拟工况组合

Table  4.  Numerical simulation case combinations

编号	试样尺寸/mm
	直径	高度
试样1	50	200
试样2	70	200
试样3	101	200
试样4	150	200
试样5	200	200
试样6	250	200
试样7	101	100
试样8	101	150
试样9	101	200
试样10	101	250
试样11	101	300
试样12	101	350

下载: 导出CSV 

| 显示表格

2.1   应力-应变关系

不同直径、不同高度应力-应变曲线(见图 6，图 7)显示，不同尺寸碎石土的应力-应变曲线在前一部分(2%应变)范围内的形状基本一致，主要区别在于曲线后段。随着试样直径和高度的变大，碎石土的峰值应力明显增大(在直径50 mm时，由于试样中最大粒径的尺寸与试样直径相差不大，剪切过程中直接形成碎石块骨架效应，该规律不显著)。峰值应力之后，应力-应变曲线出现明显的变化，直径越小、围压越低，应变软化越明显。不同试样高度时也有同样的规律，区别在于随着试样高度的增大，破坏形式由劈裂破坏转化成端部剪切破坏。当试样的直径和高度增大到一定值(直径250 mm，高度350 mm)时，应力-应变曲线几乎不再出现应变软化现象，残余应力保持恒定。

图  6  不同直径应力-应变曲线

Figure  6.  Stress strain curve of different diameters

下载: 全尺寸图片 幻灯片

图  7  不同高度应力-应变曲线

Figure  7.  Stress strain curve of different heights

下载: 全尺寸图片 幻灯片

不同直径试样峰值应力曲线(见图 8a)显示，碎石土峰值应力随着试样直径的增长呈线性增长趋势(直径50 mm除外)，在不同围压下曲线几乎平行，说明围压对尺寸效应的影响不明显；相同围压下峰值应力稳步增长，增长率为25%左右。不同高度试样峰值应力曲线(见图 8b)有所不同，曲线在前一段(试样高度小于200 mm)成直线增长，试样高度大于200 mm后，曲线明显上扬，斜率变大，应力峰值的增长随试样高度增长变得更加剧烈。从100 mm高到200 mm高的峰值应力增长率为11.9%，200 mm后的增长率为28.9%。分析其原因，是由于试样高度的增大，以及端部摩擦力的影响，使得试验在剪切变化过程中试样两端的变化大于中部，端部的碎石块先发生滚动和咬合，形成骨架效应，最终峰值应力明显升高。对比两图，试样高度在200 mm以下时对峰值应力的影响小于直径的影响，大于200 mm后高度的变化对碎石土峰值应力的影响更大。

图  8  偏应力峰值

Figure  8.  Deviatoric stress peak

下载: 全尺寸图片 幻灯片

2.2   体应变尺寸效应

对模拟过程进行体应变监测，将不同试样直径和高度时峰值应力对应的体应变列于表 5和表 6中。分析两表数据发现，围压越大，体应变越小。随着试样直径和高度的增大，破坏时体应变值都呈先增大后减小的趋势，可知碎石土试样存在某一高度和直径的阈值使其容易发生劈裂破坏。由表 5可知，体应变在试样直径为101 mm左右达到最大值，之后下降，且围压越大，下降速率越大，400 kPa时的下降速率是100 kPa时下降速率的约1.85倍；试样直径250 mm时体应变为负值，表明体应变由剪胀变成剪缩。表 6显示，不同围压下试样高度在200~250 mm时峰值体应变较大，在该高度区间，峰值体应变排在前三的分别为0.0912、0.0587、0.0568。围压在400 kPa时后期波动较大，说明试样在剪切过程中体应变不均匀，试样高度越大，端部影响更明显。

表  5  不同试样直径下峰值体应变

Table  5.  Peak body strain under different diameter

试样直径/ mm	峰值体应变
	100 kPa	200 kPa	400 kPa
50	0.0264	0.0481	0.0174
70	0.0326	0.0477	0.0452
101	0.0568	0.0587	0.0372
150	0.0286	0.02840	0.03640
200	0.0113	0.00983	0.00282
250	-0.0021	-0.01360	-0.03350

下载: 导出CSV 

| 显示表格

表  6  不同试样高度下峰值体应变

Table  6.  Peak body strain under different height

试样高度/ mm	峰值体应变
	100 kPa	200 kPa	400 kPa
100	0.0279	0.0121	0.0074
150	0.0379	0.0309	0.0151
200	0.0568	0.0587	0.0372
250	0.0557	0.0448	0.0912
300	0.0563	0.0462	0.0396
350	0.0133	0.0137	0.0888

下载: 导出CSV 

| 显示表格

2.3   粘聚力与内摩擦角影响

计算每组工况三种不同围压下的粘聚力和内摩擦角，并绘制成图(见图 9)。从图 9可以看出，随着试样直径的增大，碎石土的粘聚力几乎呈直线增长，趋势稳定；内摩擦角虽然有个别突出点，但整体上呈线性递减趋势。这是因为直径越小，碎石土中最大粒径碎石块的影响越大，碎石比表面积越大，试样主要体现碎石的性质，表现为高内摩擦角、低粘聚力；随着试样直径的增大，碎石比表面积减小，试样体现土的力学性质，表现为高粘聚力、低内摩擦角。不同高度的试样，其粘聚力和内摩擦角体现出与不同直径试样相反的特征。

图  9  粘聚力及内摩擦角

Figure  9.  Cohesive force and inner friction angle

下载: 全尺寸图片 幻灯片

从研究结果可知，碎石土的尺寸效应与粗粒土或岩石及堆石料等材料明显不同。这主要是因为碎石土组成复杂，土体内的碎石含量及分布不均匀，碎石的母岩成分不一致(软岩和硬质岩)，因此表现出区别于其他材料的不稳定的力学特性。通过上文的研究，可初步了解尺寸效应对碎石土物理力学特性的影响，为今后的三维离散仿真研究和相关工程的设计及防灾减灾提供一定的参考依据。

3.   结论

本文提出的三种块体生成法则能较好地模拟碎石土中的碎石块，通过微观参数的标定，PFC3D离散元软件能很好地实现三轴试验的数值模拟。在不同的试样直径和高度时，应力应变曲线前段(2%应变)几乎一致，曲线后段出现明显的变化。试样直径和高度越小、围压越低，应变软化越明显，当高度和直径达到一定值(直径250 mm，高度350 mm)时，曲线残余应力保持恒定。碎石土的峰值应力随试样直径和高度的增大而增大。峰值应力随试样直径增大以25%的增长率呈线性增长，随高度的增大则呈非线性增长，试样高度小于200 mm时增长率为11.9%，大于200 mm后为28.9%。围压对碎石土力学性质的尺寸效应不明显。随着试样直径与高度的增大，破坏时体应变值均表现为先增大后减小的趋势。体应变随直径的增大而下降，且围压越大，下降速率越快(400 MPa时的下降速率是100 MPa时的约1.85倍)；不同围压下试样高度为200~250 mm时峰值体应变保持最大，之后呈下降趋势。随着直径的增大，碎石土的粘聚力几乎呈直线增长，内摩擦角整体上成线性递减趋势，而随高度的增加则表现出相反的特性。