组合Clump颗粒加筋砂土三轴剪切试验离散元模拟分析

王恒通 王家全 唐毅 黄文勤

摘  要： 针对常规离散元数值模拟中采用的球形颗粒难以再现颗粒间嵌锁作用的问题，基于不同形状的Clump颗粒，利用离散元软件PFC3D开展了多组三轴剪切数值模拟试验.通过改变球形度和加筋层数研究了颗粒形状和加筋层数对砂土剪切特性的影响，并初步探讨了筋土界面作用机理.结果表明：随着颗粒球形度的增加，试样峰值抗剪强度与残余强度逐渐减小;试样中部孔隙率增加，上部孔隙率减小，其剪胀现象更加明显;抗剪强度指标c、φ值与球形度n存在二次函数关系;黏聚力随着球形度的增加减小速度变缓;而内摩擦角的减小速度随着球形度的增加变快;三轴应力强度峰值随加筋层数增加而增加，但3层加筋的峰值较2层加筋的峰值较为接近，增幅为5.9%;采用组合Clump颗粒可以有效模拟不同形状砂土颗粒，对以后的数值模拟具有一定的借鉴意义.

关键词：三轴数值模型;离散元;颗粒形状;加筋土

中图分类号：TU441                   DOI：10.16375/j.cnki.cn45-1395/t.2021.03.006

0    引言

土工合成材料由于具有造价低、易施工、易运输等优点在工程中被大量使用，因此，土工合成材料加筋土一直是研究的热点[1].三轴试验作为测定土体抗剪强度特性较为理想的室内试验之一，在科学研究和工程实际中的应用都较为广泛[2]，但随着对加筋土内部机理研究的深入，室内三轴试验显现了一定的局限性，即无法从细观层面研究土体的特性，而数值模拟可以填补室内三轴试验的不足，因此，国内外越来越多的学者采用离散元软件来模拟三轴试验.金磊等[3-4]针对含石量进行了一系列土石混合三轴模擬试验，提出了一种较为简便的柔性边界模型创建方法，此外，还发现含石量的增加能够使土石混合体局部剪切带产生不规则变化，且更加分散.Giang等[5]通过数值模型试验发现，硅砂比钙质砂更适合采用球形颗粒模拟，且颗粒级配对钙质砂的小应变剪切模量影响较大.王家全等[6-8]利用离散元PFC对黏性土进行了数值模拟研究，确定了土体的宏观抗剪强度参数与细观参数之间的数学关系.刘勇等[9]通过数值模拟分析了不同应力路径下的力学响应，发现细观参数、摩擦系数和接触强度对试样峰值强度有较大影响.徐小敏等[10]根据离散元PFC模拟三轴试验的结果，得出了接触模型为线性接触模型时颗粒材料的宏观参数杨氏模量、泊松比与细观参数颗粒法向刚度、刚度比之间的经验公式.周伦伦等[11]等采用可破碎三维离散颗粒进行了多组三轴试验，研究了不同应力状态下砂土的破碎行为.崔博等[12-14]通过设置多个参数将不同形态黏土颗粒量化，进行三轴模拟试验发现颗粒形状对土体抗剪强度有较大影响.Salot等[15]在利用离散元模拟三轴试验时采用2个颗粒连接成的簇来模拟土颗粒，得到了与室内三轴试验吻合较好的结果.

国内外的研究主要集中在颗粒粒径和级配对砂土宏观特性的影响以及细观参数与宏观参数的关系，对颗粒形状的研究较少，并且还有一定的局限性，比如定义多个参数表示颗粒形状，需要电镜扫描土颗粒外形，然后再创建不同形状颗粒，在操作上不够简便等.此外，传统的离散元方法中通常采用球形颗粒模拟砂土颗粒，忽略了颗粒形状导致的咬合摩擦，而颗粒表面凹凸产生的嵌锁作用是提供土体抗剪强度不可忽视的一部分.针对以上不足，本文基于室内三轴试验[16]，采用已开发的不同形状的Clump颗粒[17]，利用颗粒流软件PFC3D建立三轴模拟模型，研究球形度值与加筋层数对砂土的宏、细观力学特性的影响.

1    三轴试验数值模拟

1.1    Clump颗粒生成

本文采用文献[17]中Clump颗粒生成的方法.Clump颗粒的生成遵循体积等效原则和质量等效原则，颗粒间接触设置为线性接触.如图1所示，根据上述2个规则，选择一个初始球颗粒，以保证Clump颗粒体积不变，保持球颗粒A、B、C的大小不变且半径相同，同时调整球颗粒D的半径以生成不同形状的Clump颗粒.表示其不同形状的球形度值n的定义为最大内切圆与最小外接圆的比值（即球颗粒D与球颗粒A的比值），n的最大值为1.00，最小值为0.33，当[n=1.00]时为完整的球形颗粒.本次模拟试验选取的球形度值分别为0.33、0.58、0.77、1.00，具体示意图见图2.室内三轴试验土体为粗砂，平均粒径为0.62 mm，因此，对组成不同形状Clump颗粒的球颗粒半径进行调整，如表1所示.

1.2   筋材参数确定

文献[16]中所用的加筋材料为聚丙烯双向土工格栅，网孔尺寸为40 mm×50 mm，纵向抗拉强度19.2 kN/m，横向抗拉强度20.2 kN/m.为建立与试验相同土工格栅模型及再现土工格栅的抗拉性能，格栅颗粒之间的接触为平行粘结接触模型，不断调试相关细观参数，具体结果见表2，得到模拟格栅抗拉强度峰值为19 kN/m，基本可与室内试验结论保持一致.

1.3    数值模型建立及试验方案

本次模拟采用PFC3D加筋砂土三轴试验离散元模型的尺寸与文献[16]室内三轴试验尺寸一致，均为152.8 mm×200.0 mm（直径×高），如图3所示.根据不同工况需要在模型空间内生成已开发出的不同球形度砂土颗粒和相应尺寸的双向土工格栅，并在模型内布置测量球以记录相关数据，不断循环至模型达到平衡.保持轴压不变，将加载速率设定为        1 mm/min，直至模型试样轴向应变达到15%.

1.4   模型验证及模拟试验方案

图4为室内试样和数值模型在加筋层数N=1、数值模型球形度n=1.00时，围压分别为50 kPa、100 kPa时得到的应力-应变曲线.由图4可以看出，利用PFC3D模拟出的三轴剪切试验应力-应变曲线能够较好地吻合文献[16]中室内三轴剪切试验的应力-应变曲线.在达到应力峰值之前，两者发展趋势线大致重合，且取得了较为接近的应力峰值，围压为50 kPa时室内试验和数值模型试样对应的应力峰值分别为440 kPa、465 kPa，围压为100 kPa时室内试验和数值模型试样对应的应力峰值分别为780 kPa、758 kPa，可以看出相同围压下其应力峰值也较为接近.相比室内试验应力-应变曲线，颗粒流模型试样应力达到峰值有所提前，这是由于在室内加筋三轴试验时仪器的局限性使得加筋材料不能在施加荷载的瞬间立刻投入工作状态，而数值模拟不存在这一问题，故而出现了上述应力峰值提前的现象，但试样达到的应力峰值是相近的.在围压为100 kPa时，数值模拟试验在剪切后期表现出数值曲线随应变的增加软化程度持续增加，这是由于土颗粒的粒径和形状单调，不能实现室内试验中颗粒多样性的要求，但利用PFC3D的方法基本可以实现对室内三轴剪切试验的数值模拟.

为了研究球形度及加筋层数对砂土力学特性的影响，其他条件保持一致，模拟了表3中的7组   .

2   颗粒形状对砂土力学特性的影响

2.1   应力-应变及体积应变特性分析

图5、图6给出了相同围压（100 kPa）和加筋层数（N=1）下，颗粒球形度n分别为0.33、0.58、0.77、1.00时的应力-应变曲线和体应变-轴向应变曲线.由图5、图6可得：1）随着球形度的减小，应力-应变曲线峰值强度增大;球形度对残余强度影响较大，球形度的减小会使残余强度增大，同时应变软化程度逐渐减弱;颗粒球形度n影响试样的变形特性，n越大，试样的体积剪胀特征越显著，加荷初期试样发生短暂剪缩，随着应变的发展，试样剪胀显著. 2）球形度越小，Clump颗粒越不规则，相同体积下Clump的表面积越大，增强了颗粒间摩擦及咬合作用，颗粒发生滑动、旋转时其需要的轴向应力越大，即随着球形度的增大在应力-应变曲线上表现为峰值强度减小，在变形上表现为剪胀特性更加明显.

2.2   孔隙率特性分析

为研究三轴剪切试验的细观机理，分别在数值模型的中部及上部位置设置了编号为1、2的2个测量球，见图7，通过该测量球可获得局部孔隙率.图8为围压100 kPa下单层加筋时不同球形度条件下，模型试样的孔隙率随轴向应变的发展变化曲线.由图8可知，同一位置不同球形度的孔隙率变化趋势基本相同，在试样中心位置（1号测量球）处，孔隙率随轴向应变的增加均呈递增的趋势，说明孔隙率在加荷过程中一直在增大，试样砂土颗粒间的孔隙变大;在模型上部（2号测量球）随着轴向应变的增加，不同球形度下的孔隙率均逐渐减小，说明此处孔隙率在加荷过程中一直在减小，试样砂土颗粒间的孔隙变大.这是因为2号测量球处于试样顶部，在端部约束和轴向应力的共同作用下，此处颗粒密度持续增加，使得该处孔隙率持续减小，而试样中间位置处（1号测量球），受到端部约束力影响较小，在剪应力的作用下砂土颗粒间发生较大错动，使得该处孔隙率逐渐增大.此外，从图8中还可以发现随着球形度n的增加，试样中部及上部位置的孔隙率表现为增大趋势，分析原因为：土颗粒球形度的增加，使得Clump颗粒表观形状逐渐接近纯圆颗粒，使得试样在外力作用过程中的摩擦咬合作用力逐渐减弱，导致当球形度较大时土颗粒的移动和旋转较为显著，故而表现为中心部位及上部位置孔隙率均随试样球形度的增加而增大.

2.3   抗剪强度指标分析

图9、图10为不同围压下，加筋土三轴试验颗粒流模拟中土颗粒球形度分别为0.33、0.58、0.77、1.00时所获得抗剪强度指标（黏聚力c与内摩擦角φ）的变化曲线.随着球形度的增加，黏聚力与内摩擦角逐渐减小，但黏聚力前期减小幅度较大，而内摩擦角的减小幅度随着球形度的增加越来越大.抗剪强度指标（黏聚力c及内摩擦角φ）拟合成关于球形度n的多项式如式（1）、式（2）：

[c=4.294 6n2-11.636 0n+24.358 0]           （1）

[φ=-7.776 7n2+5.127 8n+53.652 0]           （2）

式中：n為球形度值，c为黏聚力，φ为内摩擦角.这两条关系曲线拟合相关系数R2分别为0.999 3和1.000 0，说明拟合效果比较好.

3    加筋层数对砂土力学特性的影响

3.1   应力-应变特性分析

根据数值模拟与室内试验的对比结果，所选参数能较合理地模拟室内单层加筋三轴试验，故采用相同的计算参数，等间距改变加筋三轴试样的加筋层数，加筋层数N分别为1层、2层和3层，进一步分析加筋层数变化对三轴宏观力学特性的影响.图11为球形度n=0.58时不同加筋层数的三轴试验应力-应变曲线.

从图11可知，试样均表现为随着轴向应变的增大，主应力差前期呈现线性增长，随后逐渐趋于平稳，达到峰值.随着加筋层数的增加，应力应变曲线的峰值强度在不断增加，但加筋层数由2层增加到3层时峰值强度的增加不太明显.当双向土工格栅加筋层数分别为1层、2层、3层时峰值强度分别为720 kPa、1 020 kPa、1 080 kPa，加筋层数从1层增加到2层时峰值强度增幅为41.7%，而当加筋层数从2层增加到3层时峰值强度增幅为5.9%.由此可见，通过增加加筋层数来提高峰值强度的效果是有限的，不可单纯地认为等量增加筋材的铺设量就可相应地提高土体的抗剪强度.与图4对比可知，由于在室内试验中受试验条件的影响，格栅在围压的作用下受到压缩，因此，在发挥作用之前格栅本身要经过恢复的过程，在此外部条件的影响下，三轴试验的数值模拟应力-应变曲线相对迟缓，但最终峰值取值仍相近似.

3.2   孔隙率特性分析

图12为4种工况下达到最大轴向应变（15%）时中心剖面的孔隙率等值线云图，颜色由深到浅代表孔隙率由小到大，中间孔隙率最大.综合4幅图的发展规律可以发现，初始孔隙率为0.35，孔隙率数值最大的部位集中在试样中部范围（y = -0.02 ～0.02 m），纯砂、单层加筋、2层加筋以及3层加筋对应的最大孔隙率分别为0.465、0.463、0.460、0.451，由此可见，不同层数的格栅对砂土颗粒侧限变形的约束作用大小不同，宏观表现为纯砂时试样中部的膨胀最为明显，而加筋层数为3层时的膨胀变形最小.

当试样为纯砂时，可以看到出现了以试样中心为圆心不断扩散的从内向外逐渐减小的孔隙率，此时试样中部砂土颗粒发生了较为剧烈的翻转、移动或是越过附近颗粒发生的侧向移动.加入1层双向土工格栅后，孔隙率取值最大的位置不再集中在中心位置，在试样高度方向y值为0附近的位置呈均匀分布的现象，且孔隙率最大值减小，究其原因：当承受较大的轴向荷载时，格栅和砂土颗粒均发生一定的位移，但格栅与砂土间互相摩擦咬合使得它们成为整体，相互抑制对方的变形和移动，故发生了如图12（b）所示较纯砂情况并不凸显的峰值位置点.加筋层数为2层时，格栅的布置位置距离试样中心有一定的距离，但由于无论在何种情况下，在承受轴压时总是试样中部发生较大的体积膨胀.因此，在y值为0附近位置试样的孔隙率较加筋层数为1层时有所增加，但变化略小.分布在y =0.04 m和y = -0.04 m处的格栅均具有一定的影响范围，当 y = 0附近的砂土颗粒即将发生较剧烈的运动时与受格栅影响的砂土颗粒相遇，格栅控制了其附近砂土颗粒的运动，进而此处的砂土颗粒将这种抑制作用传递开来，从而导致了如图12（c）所示的现象.当加筋层数为3层时，如图12（d）所示，此时格栅与砂土颗粒的结合形成了一个相对最为稳固的结合体，无论是格栅的变形还是砂土颗粒的移动都不可避免地受到了最大的调控，虽然中心部位依旧发生鼓胀变形，但加筋体的优势已明显展示出来，孔隙率的峰值为4种情况中的最低，且孔隙率数值表现为分布均匀规律，说明格栅作为“骨架”恰当地分布了试样所受的外部荷载.

3.3   砂土颗粒位移分析

图13为颗粒球形度值n=0.58，侧向围压是 100 kPa时不同加筋层数砂土颗粒的位移云图.从  图13可以看出，轴向应变达到15%时，纯砂试样中位移较小的颗粒数目显著少于加筋层数为1层时的砂土试样中位移较小的颗粒数目，与加筋层数增加时的情况一样.在相同范围内，加筋层数为3层的模型中位移较小的砂土颗粒个数要多于加筋层数为2层的砂土颗粒个数.以上现象说明加筋限制了颗粒位移且3层加筋时限制效果最强.依据等效围压的理论，在三轴试验中筋材的铺设之所以可以提升土体的剪切性能，是因为筋材对土颗粒位移的限制相当于对其周围土体有等效围压的作用，当加筋层数由2层增加到3层时，原本2层格栅控制砂土颗粒的转动和位移增加为由3层格栅加以控制，且控制范围也增大.

当加筋层数为1层时，格栅放置在试样中心位置未产生较大变形，而对比加筋层数分别为2层和3层时，发现在这2种工况下格栅发生挠曲变形，在加筋层数为3层时挠度小于2层，这说明加筋层数为2层时筋材发挥的作用最大.分析格栅产生挠曲变形原因为：该数值模拟试验的轴向加载是通过上下2个顶面加载的，大小相等，方向相反，1层格栅加筋时格栅位于试样中间，上下荷载相互抵消，在纵向的位移很小，而多层加筋时，由于试样发生剪胀，在试样中间颗粒横向位移与纵向位移较大，从而导致上下2层格栅产生较大挠度.

4    结论

1）颗粒形状对砂土的力学特性影响较大，应力-应变曲线峰值强度和残余强度均随着球形度n的减小不断增大，试样剪胀现象越明显;三轴试样中部位置孔隙率随轴向应变增大而增大，试样上部位置孔隙率随轴向应变增大而减小.

2）抗剪强度指标c、φ值与球形度n之间存在二次函数关系，均随着球形度的增加而减小，其中黏聚力随着球形度的增加，减小速度變缓，而内摩擦角的减小速度随着球形度的增加越来越快.

3）三轴应力强度峰值随加筋层数增加而增加，但3层加筋的峰值较2层加筋的峰值较为接近，增幅为5.9%，持续增加筋材层数对增加试样强度的贡献有限.

4）组合Clump颗粒相比传统的离散元模拟中的Ball颗粒可以有效模拟不同形状砂土颗粒之间的咬合摩擦，可通过改变Clump颗粒的球形度n模拟不同形状的砂土颗粒.

参考文献

[1]     王家全，张亮亮，陈亚菁，等.土工格栅加筋砂土三轴试验离散元细观分析[J].水利学报，2017，48（4）：426-434，445.

[2]     畅振超，王家全，周圆兀，等.交通动载下砾性土动三轴试验分析[J].广西科技大学学报，2019，30（2）：13-19.

[3]     金磊，曾亚武，李欢，等.基于不规则颗粒离散元的土石混合体大三轴数值模拟[J].岩土工程学报，2015，37（5）：829-838.

[4]     金磊，曾亚武.基于三维柔性薄膜边界的土石混合体大型三轴试验颗粒离散元模拟[J].岩土工程学报，2018，40（12）：2296-2304.

[5]     GIANG P H H， VAN IMPE P O， VAN IMPE W F， et al. Small-strain shear modulus of calcareous sand and its dependence on particle characteristics and gradation[J]. Soil Dynamics and Earthquake Engineering， 2017，100：371-379.

[6]     王家全，施春虎，王宇帆.加筋砂土三轴试验细观参数对宏观特性影响分析[J].广西科技大学学报，2014，25（1）：1-6，99.

[7]     陈建峰，李辉利，周健.黏性土宏细观参数相关性研究[J].力学季刊，2010，31（2）：304-309.

[8]     李可宇，杨果岳，李良吉，等.基于颗粒流模拟的黏性土宏细观参数相关性分析[J].实验力学，2020，35（6）：1147-1156.

[9]     刘勇，朱俊樸，闫斌.基于离散元理论的粗粒土三轴试验细观模拟[J].铁道科学与工程学报，2014，11（4）：58-62.

[10]   徐小敏，凌道盛，陈云敏，等.基于线性接触模型的颗粒材料细–宏观弹性常数相关关系研究[J].岩土工程学报，2010，32（7）：991-998.

[11]   周伦伦，楚锡华，徐远杰.基于离散元法的真三轴应力状态下砂土破碎行为研究[J].岩土工程学报，2017，39（5）：839-847.

[12]   崔博，邓博麒，刘明辉，等.基于不规则颗粒离散元的砾石土三轴数值模拟[J].水力发电学报，2020，39（4）：73-87.

[13]   赵洲，张鹏，宋晶，等.低围压下颗粒形态对软黏土抗剪强度影响的离散元分析[J].工程地质学报，2019，27（5）：1085-1092.

[14]   ABBIREDDY C ，CLAYTON C.The impact of particle form on the packing and shear behaviour of some granular materials： an experimental study[J].Granular Matter，2015，17（4）：427-438.

[15]   SALOT C， GOTTELAND P， VILLARD P. Influence of relative density on granular materials behavior： DEM simulations of triaxial tests[J].Granular Matter，2009，11（4）：221-236.

[16]   吴景海.土工合成材料与土工合成材料加筋砂土的相关特性[J].岩土力学，2005，26（4）：538-541.

[17]   王家全，陈亚菁，陆梦梁，等.砂土堆积试验的组合Clump颗粒离散元模拟[J].广西大学学报（自然科学版），2016，41（4）：1131-1138.

Simulation analysis of triaxial shear test of multi-sphericity Clump particles discrete element reinforced sand

WANG Hengtong， WANG Jiaquan*， TANG Yi， HUANG Wenqin

（School of Civil Engineering and Architecture， Guangxi University of Science and Technology，

Liuzhou 545006， China）

Abstract： In view of the problem that the spherical particles used in conventional discrete element     numerical simulations are difficult to reproduce the inter-particle interlocking effect， multiple sets of three-axis simulation experiments were carried out using the discrete element software PFC3D， based on the Clump particles of different shapes. The influence of particle shape and number of reinforcement layers on the shear behavior of sand was studied by changing the sphericity and number of                   reinforcement layers， and its interface interaction mechanism is preliminarily discussed. The results show that with the increase of particle sphericity， the peak shear strength and residual strength of the sample gradually decrease， the porosity of the middle part increases， and the porosity of the upper part decreases， and the dilatancy phenomenon becomes more obvious; There is a quadratic function           relationship between the shear strength index c， φ and the sphericity n. The cohesion decreases slowly with the increase of the sphericity， while the internal friction angle decreases quickly with the increase of the sphericity. The peak value of the triaxial stress intensity increases with the increase of the number of layers， but the peak value of the three-layer reinforcement is closer to that of the two-layer               reinforcement， with an increase of 5.9%. Multi-sphericity Clump particles can effectively simulate sand particles of different shapes， which has certain reference significance for future numerical simulation.

Key words： triaxial numerical model; discrete element; particle shape; reinforced soil

（責任编辑：罗小芬、黎   娅）