粗粒土大型直剪试验宏细观研究与离散元模拟

王晓帅，王子寒，景晓昆，肖成志

河北工业大学土木与交通学院，天津 300401

目前中国大型交通、水利等民生工程的建设越来越多，保证安全使用是重中之重. 作为主要填料的粗粒土压实性好、抗剪强度高、承载能力高，且在自然界广泛存在，在工程中利用率很高. 国内外岩土工作者对粗粒土基本性质与颗粒破碎等进行了深入试验研究，取得了丰硕成果[1-9]. 在数值计算方面，TEJCHMAN等[10]通过在二维有限元模型中嵌入破碎机制，模拟平面应变条件下粗粒土固定侧压力下的单向压缩实验，得到剪切带处颗粒变形破碎性质. WEI[11]对堆石坝材料进行大三轴试验，总结了粗粒土的垮塌变形规律，并通过修正南水模型进行模拟. 潘家军等[12]通过真三轴试验与平面应变试验，在不同主应力的条件下分析粗粒土的应力变形规律，为研究本构模型奠定基础. 徐晗等[13]使用ABAQUS对试验进行三维模拟，分析侧壁摩擦系数与初始应力对模拟结果的影响. 以上研究多侧重于大三轴试验与有限元模拟的结合，但粗粒土多棱角，常引起三轴仪橡皮膜破裂而导致试验失败，并且三轴仪不利于观察现象. 同时，粗粒料堆体是由离散颗粒填筑而成，采用连续介质理论的有限元模拟和真实情况有差距. 本研究通过在大直剪仪侧壁上开设透明窗口进行细观摄像，用图像处理程序[14]对细观图像进行分析，不仅可以得到土体的常规力学参数，还能对剪切过程中颗粒的细观组构进行研究，最后利用离散元对直剪试验进行数值模拟，对试样整体的宏观表现和颗粒细观运动规律进行研究，进一步研究粗粒土剪切过程中的细观力学机制.

1 试验概况

1.1 试验仪器

试验所用的主要仪器是河北工业大学生产的大直剪仪，配合透明观测窗口的尺寸，在剪切箱内设置刚性隔板，使得整个剪切试样的侧面都能被观测到. 设置隔板后试样尺寸(长×宽×高)为400 mm×600 mm×210 mm，壁厚25 mm；试验时上盒固定不动，拉动下盒进行剪切；最大垂直压力500 kN，最大水平拉(推)力300 kN；透明观测窗口尺寸(长×高)为400 mm×210 mm，满足试验要求，试验中采用单反相机(Canon 350D)进行宏细观图像摄录.

1.2 试验材料及方法

试验材料选用天然石灰岩散体颗粒．为分析粒径d对试验的影响[15-16]，通过颗粒筛分对3组不同级配的试样进行剪切测试，级配曲线如图1. 实际粗粒土工程中土通常必须压密实，相对密实度约在0.7～0.9. 本研究通过控制试样干密度，使粗粒土的最终压实度为0.8，颗粒密度为2.69 g/cm3. 受剪切盒尺寸限制，级配试样的粒径范围为0.025～5.000 cm. 试样的相关物理特性指标见表1.

图1 粗粒土试样级配曲线Fig.1 (Color online) Grading curves of coarse-grained samples

表1 试样的基本物理特性
Table 1 Physical properties of limestone samples

试样干密度/(g·cm-3)不均匀系数曲率系数粒径2～5 cm1.9051.6520.962粒径1～2 cm1.9751.4540.960级配试样2.0803.2738.182

土样分4层制备，每层压实至预定高度，装满后静置24 h，使内部应力分布均匀. 受壁厚25 mm限制，为防止剪切过程中颗粒掉落，在隔板外侧黏贴20 mm厚木板，使有效壁厚达到45 mm. 试验前将下盒由对齐位置反向推进15 mm，再进行装样．当拉出极限位移60 mm时(应变15%)，保证不漏土，对刚性接触面涂抹凡士林来消除边界摩擦. 试验先施加不同法向应力δn(静力控制)，分别为100、200、300和400 kPa，待其稳定后施加切向拉力(位移控制)，速度为1 mm/min，拉出位移达到30 mm时结束试验. 数据采样间隔为1.0 s，摄像间隔为10.0 s.

2 试验结果及分析

2.1 试验现象

通过透明观测窗口可清晰看到粗粒土剪切破坏的全过程，剪切过程中没有明显观测到颗粒发生较大破碎的现象. 由此可推断，直剪试验中绝大多数破碎以微小破碎为主，细观上表现为接触点侵蚀，造成接触刚度下降，因此粗粒土工程中应首要考虑颗粒接触点侵蚀特性，大块破碎是相对次要因素.

2.2 宏观力学特性分析

图2 剪应力-累计应变关系曲线Fig.2 (Color online) Shear stress versus cumulative strain

图2为剪应力与累计拉出应变曲线.由图2可见，曲线变化规律基本相同，都属于剪切硬化型. 剪应力随着累计应变增加而增大，每条曲线的峰值剪应力随着法向应力的增加而增大. 剪切过程伴随着颗粒的翻滚、破碎和重排列，导致曲线存在跳跃点，并不十分平滑. 级配试样曲线随法向应力的增加，累计应变相同时所需剪应力差值比较均匀，性质稳定.

试验中取峰值剪应力为抗剪强度，将抗剪强度与法向应力的关系进行线性拟合，得出粗粒土的强度参数如表2. 由表2可见，各组试验结果的离散性较小，拟合度较大. 级配试样工程性质最佳，拟合得到的假性黏聚力很小，仅5.92 kPa；单一粒径范围试样拟合产生很大的假性黏聚力，强度非线性特性明显增强，工程性质不佳.

表2 各组试样强度参数Table 2 Strength parameters of samples

2.3 颗粒定向性分析

通过透明观测窗口，利用分析程序Geodip对粗粒土剪切过程进行细观图像分析，总结颗粒长轴定向、平面孔隙率和平均配位数等细观组构的演化规律，揭示粗粒土剪切过程的细观机制.

各组试验图像分析结果显示，颗粒细观参数的演化规律基本一致，为了避免重复，本研究仅列出现象最为明显的、粒径范围在2～5 cm的石灰岩粗粒土试样，在法向应力为400 kPa时颗粒细观组构的演化规律. 其他各组试验结果不再赘述.

颗粒的长轴方向排列直接影响粗粒土的宏观力学特性.长轴在各个方向上均匀分布，则粗粒土的宏观表现为均质的各向同性，若颗粒长轴在某个方向或者几个方向集中分布，则表现为各向异性. 利用Geodip对细观照片进行分析，发现粒径在2～5 cm的石灰岩粗粒土试样，在400 kPa法向应力下不同拉出应变时刻的颗粒长轴定向演化如图3. 由图3可见，剪切开始前，颗粒长轴优选于水平方向(累计拉出应变为0)；随着剪切的进行，颗粒开始逐步旋转，颗粒长轴分布比较均匀，但水平方向仍偏多(累计拉出应变为5%)；进一步剪切，颗粒继续旋转，颗粒长轴在竖直方向的分布明显增多(累计拉出应变为10%)；试验结束时颗粒长轴分布又趋于均匀(累计拉出应变为15%)，说明剪切造成颗粒翻滚后又重新排列，达到新的平衡状态.

图3 颗粒长轴定向玫瑰图Fig.3 (Color online) Rose diagrams of long-axis orientations

2.4 平面孔隙率与平均配位数

剪切过程中，颗粒旋转和翻滚造成剪胀，同时接触点破碎成为细小颗粒，不断地填充孔隙，造成孔隙率与颗粒接触点数(配位数)的改变. 利用Geodip进行分析，粒径为2～5 cm的试样，在400 kPa法向应力下不同拉出应变时的平面孔隙率和平均配位数分别见图4和图5. 其中，平面孔隙率是根据试验中摄录图像的孔隙面积与图像面积的比值确定的. 孔隙率短暂增大后(累计应变为5%)，随着应变的继续增大而减小.平均配位数虽然在累计应变达5%以后，受剪胀作用影响增加趋势有所减缓，但仍随着应变的增大而增加；当累计应变达10%以后，细粒填充重新起主导作用，配位数显著增加. 说明在试样的剪胀过程中一直存在颗粒破碎发生，而剪胀作用仅在试验前期对细观参数有显著影响.

图4 平面孔隙率随累计应变变化曲线Fig.4 Porosity versus cumulative strain

图5 平均配位数随累计应变变化曲线Fig.5 Coordination number versus cumulative strain

2.5 位移场分析

利用图像分析程序GEODOG对大直剪试验进行位移场分析，粒径为2～5 cm试样在400 kPa法向应力下的不同应变位移场云图如图6.

图6 位移场云图Fig.6 (Color online)The nephograms of displacement fields

由图6可见，随着剪切的进行，上、下盒之间存在一条明显的水平主剪切面，在剪切过程中下盒的整体位移总是最大的. 但是剪切破坏区域却并不是始终沿着水平剪切缝发展. 随着剪切的进行，在上盒中下部逐渐形成一条向上凸起的弧形副剪切面，土体发生剪切破坏的剪切带可认为是主、副剪切面之间的拱形区域.

3 离散元数值模拟

大直剪仪的观测区域只能研究土体剪切方向的侧表面，不能对试样内部颗粒运动规律进行研究. 离散元数值模拟(distinct element modeling，DEM)是对大直剪试验的很好补充，可对试样内部颗粒细观运动规律进行全面把握，进一步研究粗粒土剪切过程中的细观力学机制.

3.1 数值模型的建立

建立模型先在目标位置生成相应的墙体，尺寸(长×宽×高)与实际直剪试验完全一致，为400 mm×600 mm×210 mm(图7).粒径范围2～5 cm和1～2 cm试样的数值模拟与直剪试验对应，在粒径范围内随机生成颗粒. 级配试样模拟生成小颗粒的数目过多，超出计算机处理能力，因此剔除粒径1 cm以下颗粒，维持级配试样中2～5 cm与1～2 cm的颗粒比例与直剪试验中相同. 使颗粒填满整个剪切盒，在重力加速度g下计算平衡.

参照文献[17]中采用的非线性接触模型，经由试验测试得出的石灰岩颗粒及刚性墙体接触参数，来模拟剪切过程中颗粒接触点的侵蚀过程，法向和切向非线性接触刚度参数分别为

(1)

(2)

Ks= 8.772×10-2·Fn·(Xs/5)-0.5

(3)

其中，Kn[A]、Kn[B]和Ks分别为石灰岩颗粒法向接触刚度、刚性墙体法向接触刚度和切向接触刚度；Xn和Xs分别为法向和切向接触位移；Fn为法向接触力. 直剪试验中在墙体表面涂抹凡士林，因此数值模拟时不考虑刚性墙的切向刚度(摩擦系数为0). 颗粒与颗粒之间的摩擦系数设为1.0. 将位于模型水平向和竖向中心线附近的颗粒染色，以便于观察颗粒运动(图7(b)).

图7 大直剪试验数值模型Fig.7 (Color online)Numerical model of large shear tests

3.2 数值模拟结果分析

图8 数值计算剪应力-累计应变关系曲线Fig.8 (Color online) Shear stress versus cumulative strain in DEM

图8为3组数值试验得出的剪应力-累计拉出应变关系曲线. 由图8可见，数值模拟与物理试验结果较为吻合，都属于剪切硬化型曲线. 粒径为2～5 cm的试样曲线在累计应变达到5%前，剪应力增长缓慢，笔者认为该组试验粒径较大，竖向边界处颗粒数目较少，边界效应明显，而数值计算中颗粒不能如物理试验一样发生破碎，未能形成整体受力骨架；后期累计应变增加，颗粒绕过关键点，重新形成受力骨架，剪应力又开始增大，与物理试验变化规律重归吻合.

表3对比了数值计算和试验测试结果. 由表3可见，粒径在2～5 cm和1～2 cm的试样数值计算出的抗剪强度和实测值相比，相对误差明显较小，最大误差率不足6%；而级配试样抗剪强度的相对误差在14%左右. 这是因为，数值计算中无法采用与直剪试验完全一致的颗粒，尤其是计算容量受限，细小颗粒无法模拟，造成一定计算误差. 数值模拟基本符合实测结果，级配试样由于缺少细小颗粒填充，计算误差相对较大.

表3 数值计算与实测结果对比Table 3 The comparison of the computed and the recorded results

各组试验剪切过程中颗粒细观运动的演化规律基本一致，为了避免重复，将现象最为明显、粒径范围在2～5 cm的数值试样，在法向应力400 kPa时的颗粒速度场及剪切滑动带示于图9. 由图9可见，当累计应变达到5%时，颗粒大都沿着弧形剪切面运动，形成连续的剪切带；此后颗粒沿着剪切带滑动，逐渐趋于稳定. 数值上则表现为，随着累计应变的增加，剪应力的增长趋势减缓.

图9 颗粒速度场及剪切滑动带Fig.9 (Color online)Particle velocity fields and shear zones

4 结 论

通过大型直剪试验与离散元数值模拟，对粗粒土剪切特性进行宏细观研究，可知：

1)粗粒土直剪试验过程中，剪应力与累计应变关系曲线属于剪切硬化型，剪应力随累计应变的增加而增大；随着法向应力的增加，峰值剪应力稳步增加；通过图像观测推断，颗粒破碎主要以棱角发生微小破碎为主，大颗粒破碎成等大两块的现象相对少见.

2)图像分析发现，剪切过程中颗粒长轴从开始优选于水平方向转变为优选于竖直方向再逐渐趋于均匀分布，表明伴随着颗粒滚动破碎，颗粒重新排列形成新的稳定状态. 细观参数的演变规律表明，颗粒翻滚与细粒填充同时作用，造成平面孔隙率先增加后减小，配位数稳定增加.

3)离散元数值模拟中，计算结果与实测结果基本一致. 单一粒径范围的试样数值计算结果较为理想，抗剪强度计算误差较小. 级配试样的计算误差略大，主要是由于受目前计算机计算容量的限制，无法模拟更加细小的颗粒，造成计算误差，使得离散元应用受限.

4)大直剪试验图像分析和离散元数值模拟中均发现，粗粒土剪切过程中(上盒固定，下盒剪切)剪切带的形成与发展并不是严格按照上下盒之间的剪切缝展开；除了存在水平主剪切面外，还存在一条向上突起的弧形副剪切面. 剪切带是主、副剪切面之间的拱形区域. 但是，如何利用粗粒土剪切过程中形成的拱形剪切带，对大型直剪试验进行强度参数修正，仍值得进一步深入研究.