0 引言

混合土作为一种特殊岩土材料，常出现在道路与铁道工程的路堑边坡中。混合土可能由几个成分组成，如黏土、粉土、砂土和砾石等。为了简化，一些研究人员将混合土视作二元混合体。粉土-砂混合料、黏土-砂混合料、砂-砾石混合料^[1]和砂-卵石混合料^[2]等均属于混合土。因为混合土由多种岩土材料组成，所以它的力学特性比单种岩土材料的更加复杂。混合土的剪切特性直接影响到路堑边坡工程的安全运营，故混合土的剪切特性需要仔细研究。

研究人员通过试验方法研究了混合土的剪切特性，他们发现混合土的剪切特性受初始孔隙比^[3]、围压^[4]、粗细颗粒粒径比^[5]、粗颗粒含量^[6-8]和粗颗粒形状^[9]等因素影响。然而，通过试验方法只能获得混合土的宏观剪切特性，很难获得细观剪切特性。离散元法的出现为获得混合土的细观剪切特性提供了一种新思路^[10-11]。研究人员^[12-15]使用离散元法研究了混合土的剪切特性，研究的影响因素包括粗(或细)颗粒含量、粗颗粒朝向、粗颗粒形状和粗细颗粒粒径比。需要注意，这些研究只关注了混合土的排水剪切特性。另有研究人员^[16-17]使用离散元法研究了细颗粒含量和粗(或细)颗粒的滚动阻力系数对混合土的不排水剪切特性的影响。然而，上两个研究中的颗粒形状为球形。综上所述，目前几乎没有离散元研究关注含真实形状粗颗粒的混合土的不排水剪切特性。因此，本研究使用离散元法探究了含真实角砾形状粗颗粒的混合土的宏细观不排水剪切特性。

1 离散元模拟概述 1.1 接触模型

颗粒间、颗粒与墙体间的接触模型均为线性模型。表 1给出了线性模型参数。Rowe^[18]发现饱和石英颗粒的摩擦系数为0.40~0.60。本研究颗粒间的摩擦系数设置为0.5。有效模量设置为1.0×10⁸ Pa，和Gong等^[19]的取值相同，近似于Guo和Zhao^[20]的取值。Goldenberg和Goldhirsch^[21]认为真实颗粒材料的法-切向刚度比的取值范围应为(1.0, 1.5)。本研究的法-切向刚度比设置为4/3。Jiang等^[22]指出阻尼系数对准静态测试无显著的影响，并通过试错法确定其取值为0.7。本研究的阻尼系数取值为0.7。

1.2 试样制备

使用离散元软件PFC^3D进行模拟。本研究中混合土的粗、细颗粒可分别视作砾石和砂。使用三维扫描仪扫描真实砾石获得其形状。在PFC^3D中，砾石轮廓内填充重叠的子球来获得模拟中的粗颗粒。所有粗颗粒的子球数量介于168~362之间。自然界中，有些砂的形状偏球形，例如Ottawa砂^[23]。此外，使用球形细颗粒能显著提高计算效率。因此，本研究使用球形细颗粒。如Zhu等^[14]推荐，本研究使用粗细颗粒粒径比4.45。粗颗粒体积和直径15 mm球体的体积相等，细颗粒直径为3.37 mm。细颗粒含量FC为0%，10%，20%，30%，40%，50%，70%，90%和100%。在生成试样前，先给定粗颗粒数量。于是，细颗粒数量可通过计算粗、细颗粒体积获得。

首先，在6面刚性墙围成的立方体空间中随机生成无重叠的粗、细颗粒。为获得各向同性的密实试样，重力加速度、颗粒间和墙-颗粒间摩擦系数暂设置为0。然后，6面刚性墙和所有颗粒向内移动。为阻止颗粒溢出该立方体空间，每10个分析步对所有颗粒的速度和角速度清零，直到运行完20 000分析步。随后，将6面墙上的目标应力设置为200 kPa，并使用伺服机制来各向同性压缩试样。最后，当不平衡力之比小于1.0×10^-4且墙上应力和目标应力的容差小于0.5%时，试样便制备完成，如图 1(a)所示。由于使用无摩擦颗粒，所有试样的相对密度均为100%，和Lopez等^[24]一致。对所有试样而言，最小试样尺寸和最大颗粒直径的比值不小于8.7。根据Cao等^[25]，本研究的尺寸效应非常有限。表 2给出了试样的基本信息。

1.3 不排水三轴压缩

耦合流体的离散元法和等体积法是研究颗粒状土不排水剪切特性的两种常用方法。对耦合流体的离散元法而言，直接模拟流体，不排水三轴压缩试验模拟中的土颗粒和流体的运动是耦合的。因此，它会耗费大量计算时间。对等体积法而言，不直接模拟流体，在不排水三轴压缩试验模拟中，试样体积保持恒定。Zhang等^[26]认为对密实颗粒状土而言，通过耦合流体的离散元法和等体积法获得的剪切特性相似，故用等体积法研究颗粒状土的不排水剪切特性是合理的。

本研究使用等体积法研究混合土的不排水剪切特性。当试样制备完成，将颗粒间摩擦系数设置为0.5。在不排水三轴压缩试验模拟中，上墙以恒定速率竖直向下压缩试样，下墙保持不动。通过水平向外移动四面侧墙使试样体积保持不变，如图 1(b)所示。先前研究^[27-28]采用和本论文相似的方法对颗粒状土的不排水剪切特性进行了数值模拟研究。为保证测试是准静态的，加载速率应足够小，其惯性参数I的取值应小于2.5×10^-3[29]。惯性参数的表达式如下：

(1)

式中，为轴向加载速率；d为所有颗粒的平均直径；ρ为颗粒密度。因上墙加载速率为0.05 m/s，故I小于1.0×10^-3。当轴向应变达到40%，测试终止。

2 结果和讨论 2.1 不排水剪切强度

使用试样的细观量将应力张量表示为^[30]：

(2)

式中，V为试样体积；N_c为接触数量；和分别为第c个接触所对应接触力的第i个分量和枝向量的第j个分量。平均应力和偏应力定义为：

(3)

(4)

式中，为σ_ij的偏张量。内摩擦角可按下式计算：

(5)

轴向应变和体应变可按下式计算：

(6)

(7)

式中，L_{1, 0}和L₁分别为压缩前和压缩中试样高度；V₀和V分别为压缩前和压缩中试样体积。

图 2绘出了应力比随轴向应变的变化趋势。随着ε₁增大，所有曲线均呈现出先增大后减小至基本稳定值的趋势。此外，细颗粒含量FC影响q/p。具体而言，随着FC增大，q/p单调减小。这说明了细颗粒含量影响混合土的不排水剪切强度。图 3绘出了峰值和临界状态的内摩擦角随细颗粒含量的变化趋势。在峰值状态，当FC≤10%时，ϕ_p基本不变；当10%≤FC≤50%时，ϕ_p随着FC增大而逐渐减小；当FC≥50%时，ϕ_p在很小的范围内变化。在临界状态，当FC≤10%时，ϕ_c基本不变；当10%≤FC≤40%时，ϕ_c随着FC增大而逐渐减小；当FC≥40%时，ϕ_c基本不变。需要注意，Cao等^[25]的研究结果也纳入了图 3中。Cao等^[25]研究了含真实角砾形状粗颗粒的混合土排水剪切特性。由图可知，加载模式(排水和不排水)对峰值状态内摩擦角影响很小，对临界状态内摩擦角有一定的影响。

考虑到计算效率，本研究将细颗粒理想化为球形颗粒，会导致细颗粒含量较高的混合土比常见真实混合土的内摩擦角偏低。为便于比较，图 4给出了先前研究和本研究的峰值内摩擦角差值对比，ϕ_p100表示纯细颗粒试样的峰值内摩擦角。由图可知，本研究的数值计算结果和Lindquist^[31]、李维树等^[32]、Xu等^[33]的试验研究结果的变化趋势大体上一致，说明本研究结果是合理的。

2.2 接触类别的作用

混合土试样中的接触可分为粗颗粒-粗颗粒(CC)、粗颗粒-细颗粒(CF)和细颗粒-细颗粒(FF)接触。应力比可写成3个应力比之和的形式：

(8)

式中，q_cc，q_cf和q_ff分别为CC，CF和FF接触的偏应力，可通过对式(2)和式(4)分别采用CC，CF和FF接触获得。根据Lopez等^[24]，考虑CC，CF和FF接触的百分比贡献，使用细颗粒含量界限值(即图 5中的FC₁，FC₂和FC₃)将混合土分为几种类型。具体而言，当FC < FC₁时，混合土为“未填满”结构。当FC₁ < FC < FC₂时，混合土为“相互作用的未填满”结构。当FC₂ < FC < FC₃时，混合土为“相互作用的过度填满”结构。当FC > FC₃时，混合土为“过度填满”结构。

图 5绘出了峰值和临界状态的每种接触对剪切强度的百分比贡献随细颗粒含量变化的趋势。需要注意，Cao等^[25]的研究结果也纳入了图 5中。Cao等^[25]研究了含真实角砾形状粗颗粒的混合土排水剪切特性。在峰值状态，加载模式(排水和不排水)几乎不影响细颗粒含量界限值。在临界状态，不排水条件下的FC₁明显大于排水条件下的FC₁，两者的差异约为10%。不排水条件下的FC₃明显小于排水条件下的FC₃，两者的差异约为5%。FC₂几乎不受加载模式的影响。上述结果表明，加载模式影响临界状态混合土的工程分类，几乎不影响峰值状态混合土的工程分类。

2.3 力-组构异性

先前研究表明试样的应力比来源于力(法向和切向接触力)和组构(接触、法向和切向枝向量)的各向异性。根据文献[20]，应力-力-组构关系表示为：

(9)

式中，a_c为接触的各向异性；a_n和a_t分别为法向和切向接触力的各向异性；a_dn和a_dt分别为法向和切向枝向量的各向异性。

图 6绘出了峰值和临界状态的整体接触网络各向异性随细颗粒含量变化的趋势。图 6(a)中，当FC≤50%时，a_c和a_t随着FC增大而逐渐减小；随着FC进一步增大，a_c和a_t基本不变。随着FC增大，a_n粗略增大。随着FC增大，a_dn先增大后减小，在FC=20%处达到最大值；a_dt逐渐减小。图 6(a)阐明了图 3中峰值状态内摩擦角依赖于细颗粒含量的细观机制。当FC≤10%时，随着FC增大，a_c, a_t和a_dt的减小程度之和补偿了a_n和a_dn的增大程度之和，造成ϕ_p基本不变。当10%≤FC≤20%时，随着FC增大，a_c, a_t和a_dt的减小程度之和超过了a_n和a_dn的增大程度之和，造成ϕ_p减小。当20%≤FC≤50%时，随着FC增大，a_c，a_t，a_dn和a_dt的减小程度之和超过了a_n的增大程度，造成ϕ_p减小。当FC≥50%时，随着FC增大，a_dn和a_dt的减小程度之和补偿了a_n的粗略增大程度，造成ϕ_p粗略不变。

图 6(b)中，随着FC增大，a_c逐渐减小。当FC≤40%时，a_t和a_dt随着FC增大而逐渐减小；当FC≥40%时，a_t和a_dt随着FC增大而基本不变。随着FC增大，a_n逐渐增大。随着FC增大，a_dn先增大后减小，在FC=40%处达到最大值。图 6(b)阐明了图 3中临界状态内摩擦角依赖于细颗粒含量的细观机制。当FC≤10%时，随着FC增大，a_c，a_t和a_dt的减小程度之和补偿了a_n和a_dn的增大程度之和，导致ϕ_c基本不变。当10%≤FC≤40%时，随着FC增大，a_c，a_t和a_dt的减小程度之和超过了a_n和a_dn的增大程度，导致ϕ_c减小。当FC≥40%时，随着FC增大，a_c和a_dn的减小程度之和补偿了a_n的增大程度，导致ϕ_c基本不变。

混合土试样中的整体接触网络可分为m个子接触网络，应力-力-组构关系也可表示为^[34]:

(10)

式中，ξ_m为第m个子接触网络的权重参数，表达式为：

(11)

式中，和N_c分别为第m个子接触网络和整体接触网络的接触数量。和分别为第m个子接触网络和整体接触网络的平均法向接触力。当m=3，忽略很小的项和，式(10)变为：

(12)

式中，ξ_cc，，和分别为CC子接触网络的权重参数、接触、法向和切向接触力的各向异性，ξ_cf，，和分别为CF子接触网络的权重参数、接触、法向和切向接触力的各向异性，ξ_ff，，和分别为FF子接触网络的权重参数、接触、法向和切向接触力的各向异性。

图 7绘出了峰值和临界状态的每种接触权重参数随细颗粒含量变化的趋势。随着FC增大，ξ_cc从1减小至0。此外，当FC≥50%时，ξ_cc非常接近0。这一结果可归因于当FC≥50%时，CC子接触网络中的接触数量和法向接触力非常小。随着FC增大，ξ_ff从0增大至1。此外，当FC≤10%时，ξ_ff非常接近0。随着FC增大，ξ_cf先增大后减小，在FC=20%~30%处达到最大值，约为0.5。

图 8绘出了峰值和临界状态的子接触网络的各向异性随细颗粒含量变化的趋势。在峰值和临界状态，CC子接触网络的各向异性中，最大。在FF子接触网络的各向异性中，最大。该发现归因于粗颗粒由重叠子球构成而细颗粒是球形颗粒。随着FC增大，，和减小；，和基本增大；，和先增大后减小。子接触网络的各向异性分析表明了图 6中整体接触网络的各向异性随细颗粒含量变化的原因。随着FC增大，a_c单调减小主要归因于的减小程度超过了的增大程度。随着FC增大，a_n单调增大主要归因于的增大程度超过了的减小程度。随着FC增大，a_t单调减小主要归因于的减小程度超过了的增大程度。

3 结论

本研究通过离散元法探究了真实角砾形状粗颗粒和球形细颗粒构成的密实混合土的不排水剪切特性。主要结论如下：

(1) 混合土的不排水剪切强度受细颗粒含量影响。在峰值状态，当FC≤10%时，内摩擦角基本不变；当10%≤FC≤50%时，内摩擦角随着FC增大而逐渐减小；当FC≥50%时，内摩擦角在小范围内变化。在临界状态，当FC≤10%时，内摩擦角基本不变；当10%≤FC≤40%时，内摩擦角随FC增大而逐渐减小；当FC≥40%时，内摩擦角基本不变。

(2) 发现了加载模式(排水或不排水)影响临界状态混合土的工程分类，但几乎不影响峰值状态的混合土工程分类。

(3) 整体接触网络的力和组构的各向异性分析表明了内摩擦角依赖于细颗粒含量的细观机制。内摩擦角随细颗粒含量的变化源自接触、法向和切向接触力的各向异性、法向和切向枝向量的各向异性的综合作用。子接触网络的各向异性分析表明了整体接触网络的各向异性随细颗粒含量变化的原因。