0.   引言

土体失稳是岩土工程的研究重点之一,传统观念认为土体在达到破坏线之前不会发生失稳,但Castro^[1]发现在不排水条件下,松砂在其应力状态到达破坏线之前便会发生类似液化的流滑失稳现象,这种失稳现象不伴随剪切带的出现,被称之为分散性失稳。失稳发生时土体未达到强度破坏线,在该状态下施加任一微小的应力增量将导致土体产生很大的变形,强度迅速降低或丧失。关于砂土分散性失稳的早期探究大多集中于不排水试验^[2-7],并且发展了相应的本构模型和数值方法用以模拟和预测这种失稳类型^[8-11]。这些研究为分析砂土的液化流滑失稳提供了重要的理论基础。

砂土在不排水条件下发生的静力液化只是分散性失稳的一个特例。随着研究的深入,人们开始探究在其他排水条件下的分散性失稳问题,常采用等比例应变试验^[12-14]。近些年的研究表明,在完全排水条件下砂土也会发生分散性失稳现象。Eckersley^[15]和Olson等^[16]的试验结果表明,在准静态和完全排水的条件下,土体可以发生类似完全不排水试验中的分散性失稳现象,并且土体在失稳前没有产生超静孔隙水压力,失稳后孔隙水压力才开始增长。其中Olson 等^[16]的研究结果表明,Wachusett水库在破坏前经历的应力路径类似于常偏应力排水剪切应力路径,偏应力保持不变有效平均主应力逐渐减小。实际工程中,地下水位上升和暴雨、雪融、灌溉造成的孔隙水重新分布均会导致有效平均主应力的降低,进而引发土体特别是斜坡土体的失稳,这种破坏都可近似看作在偏应力不变,有效平均主应力不断减小的应力路径下的破坏。在室内可采用常偏应力剪切（constant shear drained,简称CSD）的试验进行研究。

砂土在CSD应力路径下的分散性失稳研究主要存在两个问题：一是失稳点的应力状态与临界状态线相对位置；一是Hill^[17]提出的二阶功准则能否用来判定这种应力路径下的失稳。关于第一点,Sasitharan等^[18]、Skopek等^[19]、Anderson 等^[20]以及Chu等^[21]的研究结果表明在CSD试验中松砂土在应力状态到达破坏线前便发生失稳；而Lade等^[22]、Leong等^[23]、Monkul等^[24]的结果则相反。关于第二点,Sasitharan等^[18]、Skopek等^[19]认为可以采用二阶功准则判定失稳点,而Monkul等^[24]、Chu等^[25]的研究结果表明二阶功的零点与土体失稳点并不十分匹配。Monkul等^[24]认为产生矛盾的原因在于试验过程中剪应力q并未保持恒定不变,控制偏应力的误差是得出正确结论的关键。基于此,Dong等^[26-27]采用较精准的变围压三轴试验系统进行了CSD试验探究,得出了较可靠的结果。

目前CSD试验大部分集中于松砂试样^{[20-21, 28-33]},密砂试样在CSD应力路径下的力学行为探究较少^{[21, 34]}。此外受限于设备和人工操作,室内试验中很难控制加载过程中剪应力的误差,而离散元数值模拟可以实现对试样应力路径的精准控制。基于此,本文采用离散元数值方法分别对松散和密实试样在不同剪应力水平下的CSD试验进行了模拟,基于Hill^[17]二阶功准则分析了在CSD加载条件下的I型和II型两种失稳模式,分析了不同失稳模式下的应力和应变特点,研究了剪应力水平和初始孔隙比等因素对砂土失稳的影响,讨论了失稳线斜率与初始孔隙比之间的关系,并进一步提出了排水与不排水条件下砂土失稳的统一性分析方法。模拟结果有助于进一步理解砂土失稳的一般规律和内在机理。

1.   数值模拟与加载条件控制

1.1   离散元数值模型

采用三维离散元程序PFC^3D建立砂土数值模型。数值试样为圆柱形,高40 mm,直径20 mm,颗粒粒径在0.5～0.9 mm间均匀分布,平均粒径d₅₀为0.7 mm,接触模型采用接触滑动模型,DEM微观参数的选择和校验过程可参考文献[35],模拟参数见表1,生成的试样如图1所示。

表  1  离散元模拟微观力学参数

Table  1.  Micromechanical parameters of numerical samples

参数	颗粒粒径/mm	颗粒密度/(kg·m-3)	颗粒-颗粒摩擦系数	颗粒接触法向刚度/(N·m-1)	颗粒接触刚度比
量值	0.5~0.9	2600	0.75	1×107	1.5

下载: 导出CSV 

| 显示表格

图  1  t数值试样

Figure  1.  Typical numerical sample

下载: 全尺寸图片 幻灯片

试样生成后,通过伺服程序可实现将固结和CSD应力路径下的加载,亦可实现等比应变路径、常规的排水和不排水三轴剪切模拟等。

1.2   常偏应力剪切加载条件控制

CSD试验控制条件为

伺服控制过程如下：①固结后对试样进行常规的完全排水试验,排水加载至一定的剪应力水平,测出此时的偏应力q；②赋予围墙恒定的径向膨胀速率$v_{\text{r}}$；③打开轴向方向的伺服开关,将轴向的目标应力设置为${\sigma }_3+q$；④实时监测上下压盘的当前应力${\sigma }_1$,计算出与目标应力${\sigma }_3+q$之间的应力差$\Delta \sigma$；⑤实时计算颗粒与墙体的接触数和平均接触刚度,得出增益系数G；⑥计算墙体移动速率$u^{\text{w}}$,并赋予上下压盘移动速率$v_z=u^{\text{w}}$。

采用上述加载控制程序,本文对初始孔隙比为e=0.618,e=0.658,e=0.779,e=0.818的4组试样进行了CSD试验模拟。试样固结后进行完全排水剪切至预定的剪应力水平,即30,50,70 kPa,然后保持剪应力不变进行CSD路径加载,直至临界状态。需要指出的是,本文的初始孔隙比是指固结后的孔隙比。

2.   两种失稳类型

Hill认为如果材料是稳定的,那么其二阶功必须严格大于0。二阶功准则在砂土稳定性分析中得到了广泛应用^[26-32]。材料点的二阶功计算公式如下：

对于三轴试验来说,式（2）可以写成如下形式：

式中 ${{\sigma }^{\prime }}_1$为轴向有效应力；${{\sigma }^{\prime }}_3$为水平有效应力；${\epsilon }_1$为轴向应变；${\epsilon }_3$为横（径）向应变；$q={\sigma }_1-{\sigma }_3$为偏应力；${\epsilon }_{\text{v}}={\epsilon }_1+2{\epsilon }_3$为体积应变。

在CSD试验中：

因此,

结合式（3）～（6）可知

因此对于CSD试验来说,其二阶功${\text{d}}^{2}w$小于0或者失稳分为两种情况：

类型I：$\text{d}{p}^{\prime }<0$同时$\text{d}{\epsilon }_{\text{v}}>0$,

类型II：$\text{d}{p}^{\prime }>0$同时$\text{d}{\epsilon }_{\text{v}}<0$。

在CSD试验中,平均主应力是逐渐减小的,因此对于失稳类型I,其失稳点发生在临界状态线CSL之下,为分散性失稳；对于失稳类型II,其失稳点处在临界状态线之上,不属于分散性失稳,如图2所示。

图  2  CSD试验两种失稳类型示意图

Figure  2.  Two types of instability in CSD tests

下载: 全尺寸图片 幻灯片

3.   临界状态与失稳点的确定

3.1   临界状态

上述两种类型的失稳应力状态分别在临界状态应力之前和之后,因此在CSD试验中首先需要确定临界状态线CSL。试样的松密程度用状态参数$\psi =e-e_{\text{c}}$表示,$e_{\text{c}}$为临界状态孔隙比。首先通过常规的排水和不排水试验模拟来确定临界状态线,为此对孔隙比为e=0.618,e=0.658,e=0.715,e=0.748,e=0.754,e=0.767,e=0.779,e=0.818由密到松的8组试样进行了模拟,两种排水条件都加载至临界状态。模拟得到的p-q应力空间和e-p空间的临界状态线如图3所示。临界状态线斜率为0.85,在CSD试验分析中均以此值为临界状态线的斜率,模拟得到的临界孔隙比$e_{\text{c}}$为0.76。

图  3  排水、不排水和CSD试验预测的临界状态线

Figure  3.  Critical state lines (CSL) from simulations of drained triaxial, uandrained triaxial and constant shear drained stress paths

下载: 全尺寸图片 幻灯片

本文试样颗粒级配相同,固结历史相同,初始组构也相同,可以认为临界状态线具有唯一性。图3同时给出了后面进行的CSD试验得到的临界状态点,可以看出,其结果与排水和不排水三轴剪切模拟得到的结果非常吻合。

3.2   失稳点确定

由式（7）可知CSD试验中可根据体积应变随平均主应力的变化曲线来确定出失稳点,而在室内试验中当轴向应变发展速率加大导致试验完全排水条件不可控时便视为失稳。为了对比,这里也给出了轴向应变随平均主应力的变化曲线。

（1）松散试样

图4给出了孔隙比e=0.779,e=0.818两个松散试样（$\psi$分别为0.019和0.058）的体积应变和轴向应变随平均主应力的发展曲线。从图4中可以看出,两个试样在排水剪切阶段体积应变上升$\text{d}{\epsilon }_{\text{v}}>0$,同时有效平均主应力增加$\text{d}{p}^{\prime }>0$,因此二阶功大于零,试样保持稳定。随后在进行CSD加载过程中平均主应力始终减小$\text{d}{p}^{\prime }<0$,此过程中体变发展发生转折,此转折点就是失稳点。此外,从图4中轴向应变发展曲线可以看出,失稳点处曲线斜率虽不是最大,但失稳点之后轴向应变速率发展速率明显加大,这与一般室内试验^[26]结果相同。

图  4  松散试样失稳点的确定

Figure  4.  Instability points of loose samples under CSD stress path

下载: 全尺寸图片 幻灯片

(2) 密实试样

图5为e=0.618和e=0.658两个密实试样（$\psi$分别为-0.142和-0.102）的体积应变和轴向应变随平均主应力的变化曲线。试样在常规三轴排水剪切过程中体积应变少许上升$\text{d}{\epsilon }_{\text{v}}>0$,同时有效平均主应力增加$\text{d}{p}^{\prime }>0$,因此二阶功大于零,试样保持稳定。随后在进行CSD加载过程中体积应变始终减小,但此过程中体变发展发生了转折,即出现了失稳点。从图5中轴向应变发展曲线可以看出,失稳点处曲线斜率接近于无穷大,密实试样的失稳表现与室内松砂不排水试验所观察的失稳物理现象一致^[27]。

图  5  密实试样失稳点的确定

Figure  5.  Instability points of dense samples under CSD stress path

下载: 全尺寸图片 幻灯片

4.   两种失稳类型的模拟结果分析

4.1   应力-应变关系与应力路径

（1）Ⅰ型失稳

图6给出了2个松散试样的应力-应变关系曲线和应力路径,从图中可以看出,偏应力在加载过程中始终保持不变,表明本文离散元模拟实现了对试样应变路径和应力路径的精确控制。从图6中的应力路径中还可以看出,失稳点位于临界状态线以内。因此,2个松散试样在CSD加载条件下均发生了I型失稳,即分散性失稳。将失稳点连接起来即达到所谓的失稳线（instability line,IL）。

图  6  松散试样应力应变关系与应力路径

Figure  6.  Stress-strain curves and stress paths for loose samples

下载: 全尺寸图片 幻灯片

由于CSD试验模拟的失稳现象为有效平均主应力的逐渐降低,在室内试验中的控制参数通常为平均主应力的降低速率,即平均主应力均匀减小。因此平均主应力与时间等效,从CSD剪切开始至失稳点的应力路径曲线长度（平均主应力的减少量$\Delta p$）可代表失稳触发时间的长短。从图6中可以看出,对于松散试样,随着剪应力水平的提高,失稳发生的越来越快。当$q$=70 kPa时,两个松散试样均表现为在CSD剪切的起点便发生失稳；而剪应力水平相同时,e=0.818试样的失稳发生明显早于e=0.779试样,其失稳线斜率分别为0.484和0.575,这表明试样越松散其失稳发生的越快。

（2）Ⅱ型失稳

图7给出了2个密实试样的应力-应变关系曲线和应力路径发展图,从图中可以看出,失稳点后$\text{d}{p}^{\prime }>0$同时$\text{d}{\epsilon }_{\text{v}}<0$,因此2个密实试样在CSD加载条件下发生了Ⅱ型失稳,在应力路径上失稳点均越过了临界状态线,失稳线IL位于临界状态线CSL之上。

图  7  密实试样应力应变关系与应力路径

Figure  7.  Stress-strain curves and stress paths for dense samples

下载: 全尺寸图片 幻灯片

模拟结果显示两个试样的失稳均表现出随着剪应力水平的提高,失稳发生的越快。而对于同一剪应力水平,无论是在应变水平或是应力水平上,e=0.658试样的失稳发生均早于e=0.618试样。从图中还可以看出,e=0.658试样的失稳线斜率M_IL小于e=0.618试样,这与松散试样的模拟结果一致。

4.2   体积应变与孔隙比发展规律

（1）Ⅰ型失稳

图8给出了松散试样的体变与孔隙比随轴向应变的发展曲线。从图8中可以看出,试样的体变始终保持为正值,在CSD剪切初期发生了剪胀,在失稳点处达到极小值,随后发生剪缩,体积应变增加,直至临界状态。不同剪应力水平的松散试样体变曲线形态相似,但体变率变化较大,其原因在于CSD剪切起点试样孔隙比不同,这一孔隙比差别是由同一初始孔隙比试样排水剪切至不同剪应力水平所致,剪应力水平越大,其剪缩程度越大,CSD剪切起点的孔隙比越小。

图  8  松散试样体积应变与孔隙比发展曲线

Figure  8.  Volumetric strain and void ratio curves for loose samples

下载: 全尺寸图片 幻灯片

由孔隙比发展曲线可以确定试样失稳时的孔隙比。初始孔隙比e=0.779的试样在3个剪应力水平下的失稳孔隙比分别为0.774,0.765,0.758；初始孔隙比e=0.818的试样在3个剪应力水平下的失稳孔隙比分别为0.784,0.773,0.764。对于I型失稳,失稳点的孔隙比随剪应力水平的增加有也少许降低。

（2）Ⅱ型失稳

图9给出了两个密实试样的体应变和孔隙比随轴向应变的发展曲线。从图9中看出,对于密实试样排水剪切阶段的剪缩程度均很小,在CSD剪切起点的孔隙比相差不大,在CSD加载过程中两个试样均表现出剪胀,3个剪应力水平下的体变曲线变化趋势一致,剪胀率也很接近。

图  9  密实试样体积应变与孔隙比发展曲线

Figure  9.  Volumetric strain and void ratio curves for dense samples

下载: 全尺寸图片 幻灯片

从图9中还可以看出,e=0.618试样在3种剪应力水平下的失稳点孔隙比几乎相同,为0.675；而e=0.658试样在对应3种剪应力水平下的失稳点孔隙比也比较接近,分别为0.692,0.679,0.671。与I型失稳类似,对于II型失稳,失稳点的孔隙比随剪应力水平的增加同样有少许降低,当试样比较密实时这种降低几乎可以忽略不计。

5.   关于两种类型失稳的进一步讨论

5.1   CSD路径下的初始孔隙比与失稳类型

不同初始孔隙比数值试样在同一剪应力水平下模拟结果表明,在CSD加载条件下初始孔隙比大小（或状态参数$\psi =e-e_{\text{c}}$的正负）决定了砂土的失稳类型。对于密实试样,其在破坏线之前不会发生失稳,但对于松散试样则在破坏线之前发生分散性失稳,同时初始孔隙比还影响了失稳线的斜率,图10是试样的失稳线斜率M_IL随初始孔隙比的变化曲线。显然,试样初始孔隙比越大,失稳线斜率越小,不管对于Ⅰ型失稳还是Ⅱ型失稳均是如此。

图  10  失稳线斜率与初始孔隙比关系

Figure  10.  Variation of M_L of instability line with initial void ratio

下载: 全尺寸图片 幻灯片

初始孔隙比同时影响着从CSD加载开始（同时也是排水应力路径的结束点）至失稳发生时平均有效应力的减小量$\Delta p$的大小,$\Delta p$的大小反应了失稳发展的快慢。图11$\Delta p$随初始孔隙比的变化图。从图中可以看出,随着初始孔隙比的增加,达到失稳时所需的平均主应力减小量（负增量）逐渐降低,不同的剪应力水平的变化规律一致。初始孔隙比增大到一定程度（如70 kPa剪应力水平初始孔隙比为0.818的试样）,$\Delta p$几乎为零,即失稳毫无征兆地在CSD剪切的起点便发生失稳了。

图  11  CSD剪切下的平均主应力减少量与初始孔隙比的关系

Figure  11.  Variation of mean principal stress with initial void ratio under different shearing stresses

下载: 全尺寸图片 幻灯片

5.2   CSD路径中剪应力水平与失稳发展

前面的分析表明,剪应力水平对失稳线的斜率没有影响,但其影响了试样失稳发展的快慢,对于初始密实状态相同的试样,剪应力水平不同,失稳所需的平均主应力减少量和轴向应变发展不同。

图12是平均有效应力减小量$\Delta p$与轴向应变减少量$\Delta {\epsilon }_1$（其定义与$\Delta p$类似,即CSD加载开始至失稳点的轴向应变减少量）随剪应力水平的变化关系。不同孔隙比的试样其变化规律一致,以孔隙比为e=0.818的模拟结果为例,随着剪应力水平的提高,从CSD剪切开始至失稳所需的平均主应力的降低量$\Delta p$越来越少,在失稳线和排水应力路径的交点之上（见图6）,平均有效应力仅有微小的降低,失稳随即发生。对于较密实的试样,剪应力水平的增加也同样使得失稳更容易更快地发生。

图  12  剪应力水平对失稳发展的影响

Figure  12.  Variation of instability with shearing stress for samples with different initial void ratios

下载: 全尺寸图片 幻灯片

轴向应变的发展也有同样的规律,即随着剪应力水平的增加,从CSD剪切开始至失稳发生所需的轴向应变逐渐减小,不管是松散试样还是密实试样均如此。显然不管是从轴向应变还是平均主应力的降低度来看,CSD试验随着剪应力水平的增加,失稳更容易发生,特别是松散砂土在高剪应力水平下,应力路径的微小变化可能导致土体的突然失稳,这在实际工程中具有重要意义。从这个意义上说,CSD剪切之前土体所处的排水应力路径位置对可能发生的失稳有重要影响。

5.3   排水与不排水条件下砂土失稳的统一性分析

不管是传统的三轴试验还是CSD试验都是在完全排水条件下进行的,上述讨论结果都未考虑排水条件的影响或者说仅适用于完全排水条件。完全不排水条件下松散砂土也会发生I型的失稳,即所谓的静力液化,密实砂土在不排水剪切条件下由于剪胀趋势强而不发生失稳。首先讨论不排水剪切与CSD应力路径得到的失稳线之间的关系。

图13（a）是孔隙比e=0.779的松散试样在不排水条件下和其在CSD应力路径下的模拟结果。图中空心圆点是不排水剪切得到的失稳点,实心圆点是前述CSD应力路径下得到的失稳点。从图中可以看出,松散试样在不排水和CSD应力路径下均发生了Ⅰ型失稳,而且两种应力路径下得到的失稳线非常一致,其斜率的大小仅取决于初始孔隙比。换句话说,松散试样CSD应力路径的失稳线可以从其不排水剪切的失稳线得到,排水条件并未影响其失稳线的位置。

图  13  排水与不排水条件下的失稳线

Figure  13.  Instability lines in both drained and undrained conditions

下载: 全尺寸图片 幻灯片

前已述及,密实试样在不排水剪切下不发生失稳但在排水条件下可能发生失稳。图13（b）是孔隙比e=0.618的密实试样在排水剪切和CSD应力路径下的模拟结果,同样空心和实心圆点代表两种应力路径得到的失稳点。从图中可以看出,密实试样在排水剪切和CSD应力路径下均发生了Ⅱ型失稳。而且两种应力路径下得到的失稳线基本重合,其斜率的大小也取决于试样的状态参数。换句话说,密实试样在CSD应力路径的Ⅱ型失稳线可以从其常规排水剪切的失稳线得到,排水条件的不同也未影响到其失稳线的位置。

5.4   与室内试验的比较分析

Daouadji等^[28-29]对松散Hostun砂进行了三轴不排水（CU）和完全排水的CSD试验,讨论了松砂的分散性失稳问题,试验结果如图14（a）所示（图中三角形和菱形数据点）。与图13讨论的一致,两类试验的失稳线IL重合,反映了在发生I型时两类试验满足相同的内在条件。图中矩形空心和实心数据点是本文DEM模拟结果,与Hostun砂的失稳线表现出很好的一致性。

图  14  几种砂土试验结果与DEM模拟结果比较

Figure  14.  Comparison between test data and DEM results

下载: 全尺寸图片 幻灯片

Chu等^[21]对Changi砂进行了一系列三轴排水（CD）、不排水（CU）和常剪应力剪切试验（CSD）试验。图15（b）中菱形数据点是Changi密实试样的CD试验和CSD试验结果,同样与5.3节讨论的一致,CD和CSD试验的Ⅱ型失稳线基本重合。图14（b）中还给出了Feia等^[36]所做的HostunⅢ（HNM）密实砂的CSD试验结果（三角形数据点）以及本文密实试样的模拟结果（正方形数据点）。从图中可以看出,本文的模拟结果与Changi砂和Hoston砂的试验结果规律一致,图中不同砂失稳线斜率的差异反映了试样初始孔隙比的影响。

图  15  初始孔隙比与失稳线：Changi砂与DEM模拟结果

Figure  15.  Relationship between M_IL and e: Changi sand and DEM results

下载: 全尺寸图片 幻灯片

显然,室内的试验结果与本文数值模拟结果均表明：松砂发生Ⅰ型失稳,其失稳点处在临界状态线以内；而密砂发生了Ⅱ型失稳,失稳点越过了临界状态线。不管是Ⅰ型还是Ⅱ型失稳,失稳线位置与排水和不排水应力路径无关。

为了进一步分析初始孔隙比对失稳线斜率M_IL的影响,将本文得到的失稳线斜率与初始孔隙比关系与Chu等^[21]对松、密Changi砂的试验结果进行比较,如图15所示。从图中可以看出,本文DEM模拟结果与Changi砂的试验结果规律一致,大致成反“S”型的曲线关系。但需要注意的是,二者的临界孔隙比不同,模拟的结果和试验数据点并不重合。

6.   结论

本文采用离散元方法模拟了密砂和松砂在不同剪应力水平下的常偏应力剪切（CSD）试验,实现了偏应力的较精准控制。基于Hill的二阶功失稳准则分析了砂土的失稳特征,讨论了初始孔隙比和剪应力水平对CSD加载条件下砂土失稳的影响,并进一步探讨了砂土失稳的排水不排水统一分析方法。最后将模拟结果与室内多种砂的试验结果进行了比较,研究表明：

（1）Hill的二阶功准则可以用来确定CSD试验中的失稳点,根据失稳特征不同可将CSD试验中砂土的失稳分成两种类型：$\text{d}{p}^{\prime }<0$同时$\text{d}{\epsilon }_{\text{v}}>0$的Ⅰ型失稳和$\text{d}{p}^{\prime }>0$同时$\text{d}{\epsilon }_{\text{v}}<0$的II型失稳。

（2）在完全排水的CSD加载条件下,固结后的初始孔隙比决定了砂土的失稳类型：松散试样发生了Ⅰ型失稳,失稳线在临界状态线下；而密实试样发生了Ⅱ型失稳,失稳线在临界状态线上。且初始孔隙比e还影响了失稳线的斜率M_IL,本文的数值模拟结果和Changi砂试验结果一致,在e-M_IL坐标系中二者呈反“S”形曲线关系。

（3）砂土初始孔隙比越大,失稳越容易发生；同一初始孔隙比条件下,剪应力水平越高,失稳发生的越快,当剪应力q足够大时,失稳在CSD剪切起点便可毫无征兆地发生,这在实际工程中具有重要意义。

（4）松砂和密砂的CSD试验失稳线分别与同一初始孔隙比CU和CD试验失稳线重合,表明可用统一的方法分析砂土在不同排水应力路径下的失稳,室内多种砂的试验结果与本文离散元模拟结果一致,均说明了这一点。