密实散粒体宏微观特性的直剪试验离散元数值分析

蒋明镜，王富周，3，朱合华，胡海军，赵 涛

(1.同济大学地下建筑与工程系，上海 200092;2.同济大学岩土及地下工程教育部重点实验室，上海 200092; 3.上海建工(集团)总公司技术中心，上海 200083)

常规室内直剪试验存在有效剪切面积变化等缺点，许多成果对此进行了校正[1-2]，并研制了一系列适用于特殊条件下的直剪仪[3-4].由于常规室内直剪试验操作简单，测试精度能够满足一般工程要求，仍然广泛应用于科研和工程领域.

在土体的剪切破坏中，应变局部化问题已受到国际力学界和岩土工程界的广泛关注[5].然而，直剪试验中土体应变局部化区域被限定于剪切面附近，使得试样内部应力、应变分布具有明显不均匀性[6-10]，从而很难对土体内部应力、应变状态进行准确描述.Cundall等[11]提出的离散单元法，其基本模型由单元以及连接单元的变形元件构成，常用单元为刚性圆盘或三维球体.借助离散单元法可以很好地模拟直剪试验的剪破面，从而对剪切带内外土体的应力、应变特性进行分析[12-13].另外，许多研究也涉及边界效应[14-16]、颗粒形状[17-19]和颗粒破碎[20]等砂土抗剪强度影响因素.这些研究内容与试验结果吻合较好，而且可重复性较好，这对直剪试验本质规律的认识有重要作用.

对土体应力、应变状态的分析，有利于加深对不同试验条件下土体强度指标的认识.笔者近期对直剪试验剪切带的研究结果[13]表明，土体应变局部化区域出现在剪切面附近，厚度为10d50～15d50，而且初步探讨了试样内部的应力偏转问题.本文对主应力偏转进行定量研究，并以剪切带内和剪切带外的应力、应变路径为核心，分析试样内部不同区域的应力和应变状态;最后分析了颗粒群沿试样纵向的微观运动状态.

1 试样制备及离散元模拟流程

传统的离散元样本一般采用2种或3种粒径，或者某个粒径范围内的随机粒径.这种制样方法很难获得符合目标颗粒级配曲线的离散元样本，故本文采用了文献[21]提出的粒径级配离散元实施方法进行制样.对于颗粒总数为N的离散元样本，如果已知某一粒径ri在目标颗粒级配曲线的质量分数为Pi，则该粒径的颗粒数由式(1)决定:

式中s为维数，二维分析s=2，三维分析s=3.变量P由式(2)计算得到:

式中np为离散元样本中所用颗粒的种类数.

为了制备更加均匀的试样，本文采用了文献[21]提出的分层欠压法.这种方法的基本原则为:将样本在比目标空间大的区域内分层生成，按照欠压准则分层压缩当前颗粒群至准则规定的控制孔隙比，直到生成最后一层时，将样本压缩至目标孔隙比.在分n层压实时，整个压缩过程要求满足e1＞e2＞…＞ei＞…＞en= e，其中ei是按照欠压准则生成第i层颗粒时，当前所有颗粒所需压缩到的控制孔隙比.确定e1，e2，…，en的方法，即为欠压准则，取决于孔隙比的目标值e、总分层数n、每层的颗粒数目，详细过程见文献[21].经过笔者研究，发现分层欠压法能够减小层间压缩能量传递的影响，从而制备更加均匀的试样.离散元模型参数颗粒数目为24000，直径6～9mm，密度2600kg/m3，法向刚度150MN/m，切向刚度100MN/m;颗粒间摩擦系数0.5;孔隙比0.196;墙体法向刚度150MN/m，切向刚度100MN/m;颗粒与墙体间的摩擦系数0.根据分层欠压法制备的离散元试样如图1所示，颗粒级配曲线见图2.

图1 离散元试样Fig.1 DEM sample

图2 颗粒级配曲线Fig.2 Grain size distribution

试样生成后首先在12.5kPa的竖向压力作用下固结稳定，然后分别在不同的竖向压力下压缩稳定，最后进行直剪试验测试.离散元分析内容包括试样强度变化，边界摩擦效应，主应力偏转，应力、应变路径和颗粒分布与微观运动状态.

2 宏观力学特性

2.1 应力、应变特性

试样剪切过程中的体变通过试样高度变化计算得到.如图3所示，剪切刚开始时试样受到的剪应力急剧增长，并且迅速达到峰值，最后趋于残余强度.由图4所示的体变曲线可以看出，密实散粒体在剪切过程中先剪缩后剪胀，其拐弯点发生在剪应力峰值之前.本次模拟密实散粒体的宏观响应和以往经典密砂特性相吻合.

图3 剪应力-剪切位移关系曲线Fig.3 Relations between shear stress and shear displacement

图4 体变-剪切位移关系曲线Fig.4 Relationsbetween volumetric strain and shear displacement

2.2 强度特性

在直剪试验中，剪切盒内壁与土体之间一般都存在摩擦.为研究边界摩擦对试验结果的影响，本文对边界有摩擦和无摩擦2种情况进行了比较分析，2种情况下得到的强度包线如图5所示，强度指标见表1[13].对比结果表明，边界摩擦力的存在使密实散粒体的峰值强度内摩擦角提高了1.57°，残余强度内摩擦角提高了1.42°，说明实际的室内直剪试验高估了密实散粒体的抗剪强度.因为当密实散粒体试样发生剪胀时，剪切面以上的土体相对于剪切盒产生了向上运动的趋势，由于存在边界摩擦，上半盒内的试样受到向下的摩擦力，从而增加了剪切面上受到的竖向压力，密实散粒体强度得到提高.这一结果与已往研究结果相同[14-15]

图5 密实散粒体强度包线Fig.5 Shear strength envelopes of dense granular materials

表1 不同边界条件下密实散粒体强度指标Table 1 Parameters of DEM samples

2.3 剪切过程中的主方向变化

图6为密实散粒体在100kPa竖向压力作用下的主应力场的变化情况.由应力场分布的变化情况来看，未发生剪切时大小主应力分布均匀，大主应力沿竖直方向，小主应力沿水平方向;试样受到剪切后，主应力分布在大小和方向上变得很不均匀，其中试样内部9个观察点在变形后的主方向角如图6(b)所示.从主方向角的分布来看，试样内部不同位置处的应力偏转是不同的，呈现出近似反对称的特征，中心点的主方向角稳定在44°.图7为观察点4，5，6主方向角的变化过程曲线.结果表明，未发生剪切时3个观察点的主方向角在90°左右，主方向角在剪切的初始阶段变化较大，并在剪切位移达到5mm左右趋于稳定.

图6 主应力矢量分布Fig.6 Distribution of principal stress vectors

图7 主应力偏转过程曲线Fig.7 Rotation curves of principal stresses

2.4 应力路径

在文献[13]对剪切带研究的基础上，进一步对剪切带内外的应力和应变路径进行分析.图8为观察点2和观察点5的应力路径曲线，其中的强度包线是根据表1边界无摩擦时密实散粒体的强度指标换算得来的.具体换算公式为

剪切带外的观察点2，球应力和偏应力处于增大的趋势，达到峰值后表现出小幅度加卸载特性，这是因为剪切带上的总体应力软化带动了观察点2的卸荷;剪切带内的观察点5，在加载阶段达到峰值后，表现出大幅度反复加卸载的特性，这可能是因为剪切带上各点不是同时达到峰值/残余强度，存在波动所致.在偏应力的变化幅度上，剪切带内明显比剪切带外大，约为剪切带外偏应力变化幅度的2倍，说明直剪试验中的主要受剪区域位于剪切带内部.

为观察剪切带内外强度发展过程与试样整体强度发展规律之间的关系，图8给出了密实散粒体整体强度包线.可以看出，剪切带外应力路径峰值在残余强度包线之下，而剪切带内应力路径峰值在峰值强度包线之上.

图8 剪切带内外应力路径曲线Fig.8 Stress paths within and outside shear band

2.5 应变路径

根据Wang等[12]对应变局部化研究的结果，大应变主要集中在剪切带内部，其中剪切带内应变可以达到-0.7%～-0.5%，剪切带外为-0.1%～0.1%.文献[13]借助方格划分法[22]对试样剪切变形的分析表明，试样变形分布很不均匀.

图9 剪切带内外应变路径曲线Fig.9 Strain paths within and outside shear band

图9为观察点2和观察点5的体变-球应力变化曲线.由图9可以看出，剪切带外先经过一个短暂的剪缩过程，然后发展为剪胀，而剪切带内剪缩不明显，主要处于剪胀状态.剪切带外这一微弱的初始剪缩过程刚好与图4试样体变的初始变化规律一致，说明密实散粒体呈现出的初始剪缩现象主要是由剪切带外颗粒群体积变化引起.体积变化数据显示，剪切带内颗粒群的体积变化达到-0.9%，而剪切带外仅为-0.03%，两者相差一个数量级，说明直剪试验中密实散粒体的体积膨胀主要发生在剪切带内部.

3 微观力学特性分析

3.1 试样均匀性分析

剪切使得颗粒内部运动状态发生分化，这种运动状态的分布必然会对试样的均匀性产生影响.孔隙比是描述试样均匀性的一个常用参量，图10为密实散粒体在100kPa竖向压力作用下的孔隙比变化曲线.从孔隙比变化的整个过程可以看出，试样整体的孔隙比有小幅的增加，剪切带内孔隙比增大幅度较大，而剪切带外孔隙比未发生明显变化.值得注意的是，剪切带内孔隙比初始有减小的趋势.此变化规律与图9所示的宏观应变路径一致，说明密实散粒体试样初始的体积剪缩主要发生在剪切带外部区域，而试样剪胀发生局部化，主要集中在剪切面附近的条带内，这与文献[13]中观察到的剪切带一致.

3.2 颗粒平动分析

由文献[13]研究可知，剪切带在剪切面上下15倍颗粒平均粒径范围内.为深入分析剪切带内外颗粒不同的运动状态，本文将试样沿纵向进行了条带分割，共27行.

如图4所示，密实散粒体在剪切过程中高度发生变化，说明颗粒群在竖向也存在速度.对各条带的平动速度进行分解，所得的水平和竖向速度分布曲线如图11所示.与文献[13]的研究结果一致，剪切带内确实存在大的速度变化梯度(图11(a)).由图11(b)可知，颗粒群在纵向上并非单向运动，而是同时向上下2个方向运动.试样低端至第10条带范围内的颗粒，以近似相等的速度向下运动;试样顶端至第17条带范围内的颗粒以近似相等的速度向上运动;第11～16条带为过渡带，纵向速度变化较大.由于向上运动的颗粒较多且速度较大，因此试样宏观上表现出剪胀现象(图4).对比水平和纵向速度可见，颗粒水平速度比纵向速度高一个数量级，即颗粒平动以水平运动为主.

图10 孔隙比变化曲线Fig.10 Change curves of void ratio

图11 各条带的平动速度分布曲线Fig.11 Velocity distributions of each band

3.3 颗粒转动分析

文献[23]研究表明，在一个拥有N个接触点的样本中，总的能量耗散速率可以分解为平动能量耗散速率和转动能量耗散速率，而前者刚好是经典连续介质力学中的能量耗散速率，因此实际的能量耗散速率与经典连续介质力学中的能量耗散速率是有区别的.笔者在研究过程中提出了接触点处纯转动率的概念，以此来反映转动能量耗散速率，其表达式为

若考虑某一砂土颗粒群的平均纯转动率(averaged pure rotation rate，APR)，可得[23]

式中:ω——砂土颗粒群的APR，正负号代表不同转动方向;rk——第k个接触点2个颗粒的平均半径.APR可以近似视为离散微观力学和连续介质力学的桥梁[23].

图12(a)显示，转动场中剪切带边缘存在一个等势线密集的过渡带，过渡带以外颗粒基本未发生转动，过渡带内部等势线较稀疏并且转动较高，这个区域对应前面应变局部化区域.图12(b)为APR沿纵向的分布曲线，可以看出转动几乎完全集中在过渡带和剪切带内部，过渡带以外的区域转动几乎为零，说明应变局部化和转动局部化共同导致了砂土剪切带的形成.

图12 APR在试样内部的分布Fig.12 Distribution of APR in samples

需要说明的是，如图11和图12所示，在边界条件完全对称的条件下，颗粒速度和转动沿纵向的分布具有微小的不对称性，这是由于试样垂直向的微小刚体平动引起的.

4 结 语

a.离散元分析密实散粒体宏观力学特性和以往经典性质相一致，密实散粒体试样表现出了软化和剪胀特性.

b.传统的直剪试验剪切盒内壁与试样间存在摩擦力.经分析，在边界完全粗糙的条件下室内试验测得的密实散粒体内摩擦角较实际高1.5°左右.排除边界摩擦，本文分析得到密实散粒体试样峰值强度内摩擦角为22.66°，残余强度内摩擦角为17.10°.

c.剪切试样内部的变形和应力分布不均匀，不同部位产生了不同程度的应力偏转.分析发现主应力偏转发生在剪切的初始阶段，发生5mm剪切位移后趋于稳定.本文研究发现，主应力分布呈现出近似反对称的特征，中心点的主方向角稳定在44°附近.

d.对剪切带内外的应力路径研究过程可知，剪切带外主要处于加载状态并有小幅度的加卸载循环;而剪切带内表现为大幅度的加卸载循环，并且剪切带内的偏应力变化幅度为剪切带外的2倍;密实散粒体的应变路径显示，试样整体的初始剪缩现象主要是剪切带外颗粒群体积变化引起的，而试样整体的剪胀现象主要出现在剪切带内部.对应力路径和强度包线关系的研究结果表明，剪切带外应力路径峰值未超越试样整体残余强度包线，而剪切带内应力路径峰值越过了试样整体的峰值强度包线.

e.对孔隙比和颗粒运动的分析表明，剪切带是应变局部化和转动局部化共同作用的结果.试样的剪胀现象和颗粒的转动主要出现在剪切带内部，剪切带外颗粒主要处于平动状态.

由于试样强度主要来源于中间剪切面附近土体的抗剪能力，且距离剪切面远近不同土体呈现不同的力学响应，同时剪切盒内壁摩擦力也会影响砂土的抗剪强度，因此，直剪试验应尽量保持土体的均匀性，且保持剪切盒内壁的光滑性，从而保证测试结果的准确性.

致谢:孙渝刚博士生、颜海滨硕士生、肖俞硕士生参加了部分工作，特此致谢！

[1]詹金林，雷国辉，施建勇，等.有效剪切面积对直剪试验结果的影响分析[J].河海大学学报:自然科学版，2007，35(2):213-216.(ZHAN Jin-lin，LEI Guo-hui，SHI Jian-yong，et al.Effect of effective shearing area on result of direct shear tests[J].Journal of Hohai University:Natural Sciences，2007，35(2):213-216.(in Chinese))

[2]张敏江，郭尧，张丽萍，等.直剪试验中对土抗剪强度的一种修正方法[J].沈阳建筑大学学报:自然科学版，2005，21(2):78-79.(ZHANG Min-jiang，GUO Yao，ZHANG Li-ping，et al.A Revising method of shear strength of soils in direct shear test[J].Journal of Shenyang Arch and Civ Eng Univ:Natural Sciences，2005，21(2):78-79.(in Chinese))

[3]周国庆，赵光思，别小勇.超高压直残剪试验系统及其初步应用[J].中国矿业大学学报，2001，30(2):118-121.(ZHOU Guoqing，ZHAO Guang-si，BIE Xiao-yong.Super-high pressure direct/residual shear testing system and its primary applications[J].Journalof China University of Mining&Technology，2001，30(20):118-121.(in Chinese))

[4]陈正汉，馗胜霞，孙树国，等.非饱和土固结仪和直剪仪的研制及应用[J].岩土工程学报，2004，26(2):161-166.(CHEN Zheng-han，KUI Sheng-xia，SUN Shu-guo，et al.Development&application of consolidation apparatus and direct shear apparatus for unsaturated soils[J].Chinese Journal of Geotechnical Engineering，2004，26(2):161-166.(in Chinese))

[5]蒋明镜，沈珠江.土体应变局部化的研究现状[C]//施建勇，刘汉龙，蒋明镜，等.第三届全国青年岩土力学与工程会议论文集.南京:南京水利科学研究院，1998:134-149.

[6]SCARPELLI G，WOOD D M.Experimental observations of shear band patterns in direct shear test[C]//Proceedings of the IUTAM Conference on Deformation and Failure of Granular Materials.Delft:[s.n.]，1982，472-484.

[7]DOUNIAS G T，POTTS D M.Numerical analysis of drained direct and simple shear tests[J].Journal of Geotechnical Engineering，1993，119(12):1870-1891.

[8]TEJCHMAN J.Fe-simulation of direct wall shear box test[J].Soils and Foundations，2004，44(4):67-81.

[9]TEJCHMAN J，BAUER E.Fe-simulation of a direct shear and a true simple shear test within a polar hypoplasticity[J].Computers and Geotechnics，2005，32:1-16.

[10]TERZAGHI K，PECK R B.Soil mechanics in engineering practice[M].New York:Wiley，1996.

[11]CUNDALL P A，STRACK O D L.The distinct numerical model for granular assemblies[J].Geotechnique，1979，29:47-65.

[12]WANG J F，DOVE J E，GUTIERREZ M S.Discrete-continuum analysis of shear banding in the direct shear test[J].Geotechnique，2007，57(6):513-526.

[13]蒋明镜，王富周，朱合华.单粒组密砂剪切带的直剪试验离散元数值分析[J].岩土力学，2009，31(1):253-257，298.(JIANG Ming-jing，WANG Fu-zhou，ZHU He-hua.Shear band formation in ideal dense sand in direct shear test by DEM analysis[J].Rock and Soil Mechanics，2009，31(1):253-257，298.(in Chinese))

[14]LIU S H，SUN D A，MATSUOKA H.On the interface friction in direct shear test[J].Computers and Geotechnics，2005，32:317-325.

[15]LIU S H.Simulating a direct shear box test by DEM[J].CanGeotech，2006，43:155-168.

[16]刘斯宏，徐永福.粒状体直剪试验的数值模拟与微观考察[J].岩石力学与工程学报，2001，20(3):288-292.(LIU Si-hong，XUN Yong-fu.Numerical simulation for a direct box shear test on granular material and microscopic consideration[J].Chinese Journal of RockMechanics and Engineering，2001，20(3):288-292.(in Chinese))

[17]GUO P J.Modified direct shear test for anisotropic strength of sand[J].Journal of Geotechnical and Geo-Environmental Engineering，2008，134(9):1311-1318.

[18]严颖，季顺迎.碎石料直剪试验的组合颗粒单元数值模拟[J].应用力学学报，2009，26(1):1-7.(YAN Ying，JI Shun-ying. Numerical simulation of direct shear test for rubbles with clumped particles[J].Chinese Journal of Applied Mechanics，2009，26(1):1-7.(in Chinese))

[19]JIANG M J，YU H S，HARRIS D.A novel discrete model granular material incorporating rolling resistance[J].Computers and Geotechnics，2005，32:340-357.

[20]赵光思，周国庆，朱锋盼，等.颗粒破碎影响砂直剪强度的试验研究[J].中国矿业大学学报，2008，37(3):291-294.(ZHAO Guang-si，ZHOU Guo-qing，ZHU Feng-pan，et al.Experimental research on the influence of particle crushing on direct shear strength of sand[J].Journal of China University of Mining&Technology，2008，37(3):291-294.(in Chinese))

[21]JIANG M J，KONRAD J M，LEROUEIL S.An efficient technique for generating homogeneous specimens for DEM studies[J].Computers and Geotechnics，2003，30:579-597.

[22]JIANG M J，YU H S，HARRIS D.Discrete element modeling of deep penetration in granular soils[J].International Journal for Numerical and Analytical Methods in Geomechanics，2006，30:335-361.

[23]JIANG M J，YU H S，HARRIS D.Kinematic variables bridging discrete and continuum granular mechanics[J].Mechanics Research Communications，2006，33:651-666.