考虑转动阻抗作用的砂土模拟直剪试验

王 壮，段志波，廖新超，刘一鸣

(湖北工业大学 土木建筑与环境学院，湖北 武汉 430068)

砂土是人类进行工程建设的主要材料之一且自然界中分布广泛，工程实践活动都与砂土的力学特性息息相关：如建筑物地基承载力、路桥地基稳定性和土质边坡稳定性等问题[1]。直剪试验是研究土体抗剪强度及力学性能的重要方法，而室内直剪试验通常只能研究砂土宏观的力学特性，无法衡量细观层面砂土颗粒之间相互作用对宏观力学变形特性的影响。Cundall等[2]基于牛顿第二定律和力位移方程提出了针对砂土等理想散体材料为研究对象的离散单元法。随着国内外学者对离散元深化研究及计算机算力的巨大提升，离散单元法在研究土体细观机理方面得到广泛运用。基于离散单元法，学者们做了大量模拟试验，并取得了许多成果：朱遥等[3]探究了颗粒形态对砂的抗剪强度的影响，发现不规则形态砂颗粒具有较大抗剪强度；潘远阳等[4]研究了剪切带形成的细观机理；孔亮等[5]分析了颗粒大小和形状对砂土力学特性的影响，并构建了颗粒外轮廓特征的性状系数。文献[6]表明当颗粒间的摩擦较大时，决定颗粒材料力学机制及微观变形机制的主要因素为颗粒之间相对转动作用。传统的离散元模拟将颗粒设为球形或圆形，忽略了颗粒转动的影响。Iwashita等[7-8]均对离散单元法本构模型进行了改良，考虑了颗粒转动的影响，并引入了转动阻抗的概念。Jiang等[9]提出了抗转动线性接触模型，基于三维接触模型，设定颗粒间的接触为面接触，并考虑了滚动阻力和扭转阻力的影响。Zhu等[10]基于离散单元法，通过模拟三轴试验分析了剪切强度随转动阻抗和细粒含量变化的机理。Zhao等[11]综合比较了滚动阻力和颗粒形状对颗粒材料内部剪切诱发的组构和各向异性的影响，认为滚动阻力模型可有效再现砂土中剪切诱导各向异性主要特征。而对转动阻抗的选取则需综合考虑颗粒的实际形状和颗粒间的相对间距，总体而言对于使用纯圆颗粒建模的试验，实际颗粒形状越接近球形，相对间距越小，需选取的抗转动系数越小[12-13]。

本文利用离散元软件PFC3D对砂样直剪试验进行一系列标准砂直剪数值模拟试验，考虑不同围压下转动阻抗对砂样直剪过程中的宏观受力及变形特性影响，细观层面则通过接触力链和配位数分析颗粒接触状态和剪切带发育的演化规律；结合玫瑰图、偏组构及组构张量分析各向异性演化规律。通过对砂土宏细观力学结合分析，表明转动阻抗作用是一个不容忽视的因素，对揭示砂土剪切的本质机理提供理论参考。

1 离散元建模

1.1 抗转动线性接触模型

抗转动线性接触模型在考虑颗粒转动影响时已经取得良好运用[10-11]。由于砂颗粒材料自身各向异性不明显及不考虑砂颗粒与墙体之间的边界效应，综合选取砂颗粒-墙体接触类型为线性接触模型，砂颗粒-砂颗粒接触类型为抗转动线性接触模型（图1）。线性模型的接触点可以想象为一对弹性弹簧，线性力由具有恒定法向刚度和切向刚度的线性弹簧产生，主要包括法向接触力和切向接触力，其力学原理为：

图1 接触模型Fig. 1 Contact model

在线性模型的基础上，抗转动线性接触模型[7]增加了转动阻抗力矩，其定义为：

式中：kr为 颗粒滚动刚度；R¯为 有效接触半径；Ij为j颗粒惯性矩；rj为j颗粒半径为与j颗粒接触q颗粒的切向力为与j颗粒接触q颗粒的扭矩；R1和R2分别为颗粒q、j的半径。

式中：μr为 抗转动系数；θb为接触颗粒间的相对转动为抗转转动线性接触模型中的法向接触力。

1.2 数值模型及细观参数

离散元建模过程中颗粒数量会影响模拟的结果和计算速率，在保证结果可靠前提下提高计算效率，通常采用扩大颗粒粒径R或者减小试样尺寸L的方法，文献[14]中数值平面应变试验结果表明，当L/R=30～40 时才能保证尺寸效应的影响最小。刘海涛等[15]通过一系列三轴试验的颗粒流模拟，认为当L/R＞40时，尺寸效应对粗粒土抗剪强度没有影响。尹小涛等[16]通过模拟不同内尺度比岩土试样单轴压缩试验，发现试样中的最小粒径dmin与数值试样的长度之比小于0.01时，数值试验的结果具有较好稳定性。综上所述，砂土离散元试样采用尺寸与室内大型直剪仪尺寸一致，尺寸为 305 mm×305 mm×200 mm，砂颗粒粒径相对实际值放大10倍并删除部分小粒径，最终取最小粒径dmin=1 mm，最大粒径dmax=12 mm，这可同时满足dmin/L≤0.01，L/R≥40。按初始孔隙率0.39随机生成颗粒43 775个。离散元数值模型分别采用5块刚性墙体模拟上下剪切盒，并生成2道水平刚性挡板防止剪切过程中的颗粒溢出（图2）。

图2 数值模型（单位：mm）Fig. 2 Numerical model of sand （unit：mm）

通过对比研究前人对砂土建模细观参数的选取，初步确定砂细观参数取值范围，经过反复调试最终确定细观参数见表1。

表1 砂细观参数Tab. 1 Sand parameters

砂土直剪数值模拟主要包括制样、固结和剪切3个阶段，本研究首先采用孔隙率制样随机生成颗粒并消除不平衡力影响，然后对试样施加不同围压固结，最终剪切阶段保持法向应力不变，对下盒施加一个向右的剪切速度模拟砂土剪切。分析法向应力分别为50、100、150和200 kPa时不同抗转动系数条件下砂土直剪的宏细观力学特性。

2 宏观力学特性

2.1 应力-应变响应

由图3中剪应力-剪切位移曲线可知，曲线初始斜率基本一致，随着抗转动系数的增大各曲线峰值随之增大，而增幅逐渐降低；竖向位移-剪切位移曲线均表明试样先微弱剪缩后剪胀的变形特性。相比于抗转动系数为0的试样，其他试样均表现出应变软化特征，抗转动系数越大，软化越明显，且在剪切后期仍有明显剪胀发生。

图3 剪切应力-剪切位移及竖向位移-剪切位移曲线Fig. 3 Curves of shear stress-shear displacement and volume change-shear displacement

以上规律表明：转动阻抗对砂土试样的初始弹性模量无明显影响；对试样的峰值强度有显著增强作用，但随着μr的增大，转动阻抗的影响逐渐趋于饱和；在转动阻抗的作用下，砂体更易表现出应变软化现象，且随着μr的增大，软化逐渐加剧。分析主要原因为：由于转动阻抗的作用，颗粒间发生相对转动需要克服的阻力变大，造成宏观抗剪强度的增大；而在转动阻抗作用下试样表现出的应变软化现象符合中密砂的性状。

2.2 抗剪强度指标

将模拟结果抗剪强度包线与王军等[20-23]所做室内标准砂直剪试验拟合包线进行对比。由图4（a）可知，当抗转动系数为0时，抗剪强度包线与前人室内试验结果基本吻合且为直线，符合库伦强度准则的一般规律，验证了模拟试验的可行性与准确性。抗剪强度包线的斜率随抗转动系数的增大而增大，即转动阻抗的作用可有效增大峰值内摩擦角。对于纯圆颗粒，可根据实际颗粒的不规则程度及实际试样的宏观力学特性选取抗转动系数。

图4 抗剪强度包线对比Fig. 4 Comparison of shear strength envelope

2.3 剪切带

为了探究剪切过程中试样内不同位置的变形，在试样中均匀设置了9个长方体区域，并将其中砂颗粒设置为黑色，监测直剪过程中各组砂颗粒移动情况。对比图5（b）（c）可知，相比于初始阶段，剪应力达到峰值时，靠近两侧的砂颗粒出现相对位移，表明试样内砂颗粒的相对位移由左右两侧开始，变形由外向内发展，最终变形集中于一条带状区域内，即剪切带（如图5（d）），其厚度为79.2 mm，为平均粒径的13～15倍，与潘远阳等[4]研究结果相吻合。

图5 试样变形演化Fig. 5 Evolution of specimen’s deformation

3 细观结果分析

3.1 接触力链

力链网络是土骨架体系上粒间接触力的宏观形式，可体现颗粒体系受外荷载作用时受力响应机制的敏感性。图6为初始状态、峰值状态和最终状态试样内力链的分布情况，其中力链直径与法向接触力的大小成正比。初始状态，由于加载尚未进行，伺服后试样内接触力链分布均匀，各向异性不明显。达到剪应力峰值时，强力链向剪切带区域靠拢，且沿直剪仪的斜对角分布。达到最终状态时，强力链更加集中，且主要分布于剪切带附近。

图6 接触力链网络演化Fig. 6 Contact force chains in specimen

为进一步定量分析试样内部接触力演化规律，依据文献[24]提出的强弱接触力系统，分别对不同条件下试样的法向和切向的强接触力比例进行量化统计。图7和8分别为法向应力为50 kPa时法向强接触力和切向强接触力的演化规律。由图7（a）、8（a）可知，法向与切向强接触力比例均随着剪切进行先降低，当达到极值后逐渐趋于稳定；抗转动系数越大，强接触力比例越小。这是由于转动阻抗作用越大，颗粒间抗转动能力越强，颗粒间传递力矩越大，从而使强接触力更大，而强接触力主要抵抗外部荷载的作用，因此当外部荷载相同时所需的强接触力数越少。由图7（b）、8（b）可明显看出，随着围压的增大，强接触力比例极值呈微弱增大趋势。

图7 法向接触力变化规律Fig. 7 Variation of normal contact force

图8 切向接触力变化规律Fig. 8 Variation of tangential contact force

3.2 配位数

配位数是指试样中某颗粒与周围颗粒的接触数目，描述颗粒材料内部颗粒的接触程度。配位数又称平均接触数目，本文采用力学配位数[25]，其定义如下：

式中：Nc为试样内部总接触数：Ns为试样内部总颗粒数；Ns0为某颗粒周围无接触的颗粒数；Ns1为某颗粒周围接触数目为1的颗粒数。

图9 为200 kPa法向应力条件下的平均力学配位数随剪切位移的变化曲线，可看出配位数均随剪切进行而降低，后趋于稳定。随着抗转动系数的增大，配位数的稳定值随之减小，这是由于转动阻抗的作用使颗粒不易发生相对转动，进而使颗粒间形成较大的接触力，抗转动系数越大，传递能力越强，平均有效接触就越少。

图9 配位数随剪切位移变化曲线Fig. 9 Variation curve of coordination number with shear displacement

3.3 各向异性

接触组构[26]是反映颗粒散料微细观各向异性的常用指标。剪切过程中对试样内颗粒间的法向接触力和切向接触力的数量及大小进行统计分析，图10和11展示了200 kPa法向应力峰值强度时颗粒间法向接触力和切向接触力的各向异性分布玫瑰图，图中每10°为一个区间，每个区间内接触力的平均值与扇形面积相对应。可明显看出，法向接触力分布玫瑰图呈花生状，切向接触力分布玫瑰图呈花瓣状，与有关文献[27]研究结果相吻合。

图10 法向接触力分布Fig. 10 Distribution of normal contact force

图11 切向接触力分布Fig. 11 Distribution of tangential contact force

对比不同抗转动系数条件下的玫瑰图可知，随着抗转动系数的增大，“花生”变得更加细长，“花瓣”亦表现出类似变化，这表明转动阻抗的作用使法向接触力和切向接触力更加集中分布于部分区间，增大了这部分区间内的接触力平均值，进而增大了各向异性程度。此外，随着抗转动系数的增大，法向接触力与切向接触力的主方向均发生了顺时针偏转，其中法向接触力偏转相对明显。

为进一步定量分析颗粒材料的各向异性，参考组构张量的概念[28]

其中：Fi j为组构张量在不同方向上的分量表示组构。在此基础上，张铎[29]根据广义偏应力概念采用组构张量定义“组构偏量”，即：

图12 （a）为法向应力是50 kPa时，偏组构Fq随剪切位移的变化曲线，可以看出随着抗转动系数的增大，偏组构峰值随之增大，即试样整体的各向异性随之增大，但随着剪切进一步进行，抗转动系数不为0的试样的各向异性均逐渐降低，类似于应变软化趋势，且抗转动系数越大，降低幅度越大。依据文献[24]提出的强弱接触力系统，将颗粒间接触分为强弱接触。图12（b）给出了强接触偏组构随剪切位移的变化曲线，其规律与整体偏组构类似，但在整个剪切过程中，其大小均大于整体偏组构，这表明转动阻抗对试样整体各向异性的增强作用主要是增加强接触的各向异性。分析主要原因：转动阻抗增大了强接触传递力的能力和抵抗力链破坏的能力，使强接触处颗粒可传递更大的力，同时减小了接触数量，即使接触力的传递更加集中，各向异性程度更大。图12（c）为强接触偏组构峰值随法向应力的变化曲线，可见随着法向应力的增大，不同抗转动系数试样的各向异性均呈下降趋势。

Oda[30]提出采用组构张量描述颗粒材料内的接触特性，通过分布密度函数E(Ω)定义组构张量，并采用组构张量的第二不变量定义其各向异性系数ac。

式中：φi j为 组构张量；ni、nj为接触法向在笛卡尔坐标系中的分量为描述全局接触法向分布密度的参量；ac为的第二不变量。

与组构张量各向异性系数的描述类似，法向接触力张量也可定义其各向异性系数an。图13和14为各向异性系数随轴向应变的变化曲线。随着剪切进行，试样发生膨胀变形，ac和an的演化规律十分相似，均表现为先急剧增加，相比于无转动阻抗作用的试样，其他试样达到峰值后随轴向应变显示出软化趋势，且抗转动摩擦系数越大，软化趋势越明显，与剪应力-剪切位移曲线相对应；对比图13和14可知，整个剪切过程中an大 于ac，这是由法向接触力大小的非均匀分布导致的，同时表明法向接触力各向异性对试样抗剪强度影响较大；随着抗转动系数的增大，ac与an均随之增大，即在砂试样中产生更显著的剪切诱导各向异性，从而提高了砂的稳定性及抗剪强度。

图13 接触法向各向异性系数αc随剪切位移变化曲线Fig. 13 Curve of contact normal anisotropy coefficient αc with axial strain

图14 法向接触力各向异性系数αn随剪切位移变化曲线Fig. 14 Curve of normal contact force anisotropy coefficient αn with axial strain

4 结 语

（1）数值试验宏观应力-应变特性、强度指标及剪切带变形规律与前人研究结果相类似，表明砂土直剪数值模拟试验存在良好的可行性；相比于无转动阻抗作用的试样，转动阻抗的作用使砂表现出明显的应力软化及剪胀特性，且抗转动系数越大，剪应力峰值越大。

（2）随着剪切进行，力链分布区域化，达到峰值后集中于剪切带区域，且发生了主方向的偏转。由于转动阻抗增大颗粒间传递接触力的能力，外部荷载相同条件下，抗转动系数越大，所需接触数越少，即表现为平均力学配位数随抗转动系数的增大而减小，且强弱接触力比例随之更加悬殊。

（3）玫瑰图表明转动阻抗的作用增强了法向接触力和切向接触力的各向异性，偏组构的变化规律说明转动阻抗主要增大了强接触力的各向异性，而张量的第二不变量演化规律则表明转动阻抗对法向接触力各向异性的增强作用对提高宏观抗剪强度具有显著作用。