在离散元软件内建立二维模型，模型剪切盒尺寸为长800 mm、宽600 mm，剪切面位于剪切盒中部，剪切过程中垂直方向施加法向压力，下剪切盒施加切向剪应力。在剪切过程中，上下剪切盒接触面积不断减小，因此在计算剪应力的过程中考虑了剪切面面积的变化，模型示意图如图1所示。

图1 直剪模型示意图Fig.1 Schematic diagram of direct shear test model

颗粒的生成方法采用了颗粒膨胀法。为使得生成颗粒的重叠量最小，颗粒系统更加均匀化，在制样和颗粒预固结时，颗粒与颗粒或与墙体之间的摩擦系数设置为零。在试样的剪切阶段再将摩擦系数设置为相关值。本文模型的细观参数值是来自于Wang等^[12]中研究选取的颗粒参数值。虽并未经过严格的参数标定，但Wang等利用该参数取得了一定研究价值的结果，证明了这些参数在模拟中的价值。因此，本文选取这些参数对不同砾石含量的砂土开展直剪模拟研究。

表1 模拟中颗粒的细观参数Table 1 Micromechanical parameters of particles in the simulation

根据牛首山东麓商业配套综合项目在勘察资料中砾石的含量及尺寸，本次数值模拟重点研究砾石含量在0%~25%的砂土力学行为变化规律，实际工程获得的样品的最小粒径可达0.001 mm，颗粒级配范围很大，在数值模拟中完全按照实际颗粒级配进行建模将导致颗粒数过多从而无法进行计算。因此本次模拟中参考王磊^[13]的方法开展，最小粒径设定为2 mm，以5 mm以下的颗粒代表砂土颗粒，20 mm以上的颗粒代表砾石颗粒，模拟中的砾石含量分别为0%， 5%， 10%， 15%， 20%， 25%，采用颗粒膨胀法生成了六组不同的级配，各含量级配如表2所示。模型建立过程中，砾石在剪切盒内随机分布。不同砾石含量的模型如图2所示。

表2 模拟中粒组含量表Table 2 Grain size fraction contents in DEM simulation

图2 不同砾石含量砂土模型Fig.2 Numerical models of gravelly sand with different gravel contents

直剪模拟过程中，法向压力分别设置为100、200和300 kPa。为了确保模拟试验是在准静态条件下进行，剪切盒的剪切速率应该足够小，以保证在剪切过程中产生的动能忽略不计。本研究采用Da Cruz等^[14]提出的无量纲惯性参数I来量化这种效应，当I<1.0×10^-3时，即进入准静态临界状态。因此，离散元中的所有的剪切模拟试验惯性数I均保持在I<1.0×10^-4以下，即在准静态条件下进行。采用位移控制的剪切模式，剪切速率设定为2 mm/s，剪切位移达到50 mm时终止。

模拟分别得到了100、200和300 kPa三种法向压力6种砾石含量条件下的剪应力与剪切位移之间的曲线关系，如图3所示。

图3 剪切位移与剪应力关系曲线Fig.3 Shear stress-displacement relationship curve

从图3中可以看出，三种法向压力6种不同砾石含量条件下各组模型剪切位移与剪应力曲线的趋势基本一致，随剪切位移的增大剪应力均呈现先增大后减小并趋于平稳，出现应变软化现象即剪切过程中抗剪强度先达到峰值后逐渐降低的现象，这主要是因为剪切初期因颗粒间咬合紧密或结构强度高，抗剪强度达到峰值，但随着剪切进行，颗粒咬合破坏或结构崩解，强度逐渐降低，最终趋于残余强度。

不同法向压力下，抗剪强度与砾石含量的曲线如图4所示，从图中可以看出抗剪强度均表现出砾石含量为0%时最高，5%~20%之间是抗剪强度变化值较小，当砾石含量达到25%时，其抗剪强度陡然下降。在100 kPa法向压力下，抗剪强度从砾石含量为0时的116 kPa下降至99 kPa，下降幅度为15%；在200 kPa法向压力下，抗剪强度从砾石含量为0时的216.5 kPa下降至176.7 kPa，下降幅度为18.4%；在300 kPa法向压力下，抗剪强度从砾石含量为0时的294.3 kPa下降至246.3 kPa，下降幅度为17.3%。这主要是因为大直径砾石含量的增加，使得细颗粒砂土无法充分填充粗颗粒间的空隙，导致混合体形成“点接触骨架”而非“密实嵌锁结构”。此时，粗颗粒间缺乏细粒土应力的传递作用，颗粒易滑动或滚动，整体摩擦阻力下降，从而导致抗剪强度降低。并且从图4中可以看出， 25%的砾石含量为抗剪强度的阈值点，当砾石含量超过25%时其抗剪强度将显著下降。

图4 不同法向压力下砾石含量与抗剪强度曲线Fig.4 Relationship between gravel content and shear strength under different normal pressures

而残余抗剪强度与抗剪强度相比则表现出不同的演化特征，不同法向压力下残余强度基本为抗剪强度的60%~70%。随着砾石含量的增加，残余抗剪强度呈现出先平缓后下降的趋势，当砾石含量达到20%，残余抗剪强度均出现明显的下降，因此，20%的砾石含量为残余抗剪强度的阈值点。

剪应力比是剪应力/法向应力与剪切位移的关系曲线，它是归一化曲线，消除了不同法向应力的影响，因此在分析其剪切力学特性时更具有代表性，当剪应力达到最大值时，相应的得到该种模型下最大剪应力比。图5为不同砾石含量下最大剪应力比与法向压力的关系曲线。从曲线中可以看不同砾石含量条件下，随着法向应力的增加，最大剪应力比均呈现下降趋势，以砾石含量5%为例，法向压力为100 kPa时，最大剪应力比为1.16，法向压力为200 kPa时，最大剪应力比下降至1.08，法向压力为300 kPa时，最大剪应力比再次下降至0.98，最高下降约16%。这说明含砾砂土在高应力条件下表现出非线性强度特性。

图6为拟合计算得到的不同砾石含量下的砂土内摩擦角，可以看出，砾石含量为0时，内摩擦角最大为41°，砾石含量为25%时，内摩擦角最低下降至36°。这主要是因为砾石含量较高时，粒径分布趋于两极分化，缺乏中间粒径颗粒的过渡，导致颗粒间咬合能力下降，剪切面绕石块边缘发展，形成绕流破坏，砾石颗粒作为刚性体阻碍连续剪切带形成，导致能量消耗在局部颗粒重排而非整体摩擦，表观内摩擦角下降。在实际工程中计算边坡稳定性时，应该考虑随着砾石含量的升高对内摩擦角进行相应的折减。

图5 不同砾石含量下法向压力与剪应力比关系曲线Fig.5 Relationship curves between normal pressure and shear stress ratio with different gravel contents

图6 内摩擦角与砾石含量关系曲线Fig.6 Relationship between internal friction angle and gravel content

力链是剪切过程中，颗粒见接触力高度集中的链状网络，形成承载骨架，形象地表征了样品颗粒之间力的传递途径，力链在剪切过程中不断的断裂与重组。不同砾石含量条件下，样品剪切过程中表现出不同的力链演化特征。图7为法向压力100 kPa时不同砾石含量的砂土在剪切位移达到最大值时力链图像，黑色线条为颗粒与颗粒之间的接触力，线条的粗细则代表颗粒与颗粒之间力的大小，黑色线条越粗，代表颗粒间接触力越大，为颗粒间接触力的主要传递路径。当砾石含量为0时，力链分布较均匀，主力链分布在上剪切盒右侧靠近剪切面处，与剪切面的夹角约为45°，随着砾石含量的增加，可以明显看出力链分布特征发生变化，主力链主要分布在砾石颗粒周围，说明力主要通过大颗粒的砾石传递，在砾石颗粒的周围存在一定的空白区域，表明砾石颗粒周围出现架空现象不能有效对力进行传递，这也从微观角度解释了随着砾石含量的增加试样的抗剪强度下降的原因。

图7 法向压力为100 kPa时不同砾石含量砂土剪切力链分布图Fig.7 Distribution of shear force chains in sand with different gravel contents at 100 kPa normal pressure

在剪切过程中，由于颗粒不断运动导致剪切面附近的孔隙比发生变化，孔隙比反映了试样剪切过程中颗粒的体积变化，并且与试样的抗剪强度有着密切的关系。

图8为法向压力为100 kPa时不同砾石含量的砂土剪切面附近某一局部区域内孔隙比随剪切位移的关系曲线。从图中可以看出随着剪切位移的增大，该区域内的孔隙比基本均呈现升高的趋势，不同砾石含量条件下，孔隙比增量最小约30%，最大增加约60%。这主要源于其独特的颗粒组成与力学响应机制。当砾石形成连续骨架时，剪切位移迫使砾石发生滚动或滑动，导致原本咬合紧密的骨架解体。大颗粒砾石之间的相对运动形成局部架空孔隙，小颗粒无法及时填充，孔隙比显著增大。剪切过程中，小颗粒在剪应力作用下被挤入相邻低应力区，导致剪切带内小颗粒占比降低，粗颗粒间残留孔隙未被填充。局部密实小细粒区域会发生剪胀作用，尤其在砾石颗粒边缘的小颗粒因受约束较弱，颗粒翻越咬合导致孔隙比局部增大。含砾石砂土剪切过程中剪切面孔隙比增加的根源在于骨架解体—细粒迁移—局部剪胀的协同作用，这一机理揭示了“微观孔隙演化—宏观力学劣化”的内在关联。

图8 法向压力为100 kPa时不同砾石含量砂土剪切面处剪切位移与孔隙比关系曲线Fig.8 Shear displacement versus void ratio at shear plane for sand with different gravel contents under 100 kPa normal pressure

砾石颗粒引发的力链网络异化与剪切带孔隙比重构形成正反馈，粗颗粒主导的应力传递路径加剧了局部应变集中，导致剪切带持续扩容与整体结构失稳。