2020年国内开通地铁的城市将达到32个，其中大部分城市在滨海滨河地区，地下覆盖深厚淤泥质软黏土. 现有监测数据显示^[1-2]，淤泥质软黏土由于水的质量分数大、压缩性大，在地铁列车长期运行作用下容易产生不均匀、难稳定的沉降. 目前，针对地铁列车荷载长期作用下的饱和软黏土动力特性研究，通常采用室内试验方法^[3-5]. 在模拟列车荷载方面，传统室内研究通常沿用地震荷载的振动方式，采用连续的基础波形循环往复振动上万次^[6-10]. 实际上，地铁列车车速约为40~80 km/h，8节车厢列车组的长度约为180 m，可以估算出列车全速经过某点所需时间约为8~16 s，每趟列车之间有几十秒到几分钟不等的间隔. 因此，如采用连续不间断的振动方式模拟列车荷载长期作用，试样的受力情况与实际工程有较大区别.

目前，国内外学者针对连续振动长期作用下软黏土应变的研究已取得较为丰富的成果^[11-13]. Parr^[11]分析不排水条件下Kaolin重塑黏土试验结果，发现各振次的塑性应变增长率 ${\dot \varepsilon _N}$对第1次塑性应变 ${ \varepsilon _{\rm{1}}}$的比值与循环应力加载周期的关系在双对数坐标系下呈线性相关.基于Parr^[11]的研究，黄茂松等^[12]对上海淤泥质软黏土进行不同先期固结围压p_c下的不排水循环三轴试验，发现归一化塑性应变率 ${\dot \varepsilon _N}$/ ${\dot \varepsilon _{\rm{1}}}$和加载次数N的双对数曲线斜率b一致分布在−0.644左右，与前期固结围压、静偏应力水平、动偏应力水平无关；第1次加载塑性应变 ${\varepsilon _{\rm{1}}}$反映连续振动循环作用的围压、动应力大小和应力历史等其他因素的全部影响. 该发现肯定了归一化应变增长速率 ${\dot \varepsilon _{\rm{N}}}$/ ${\dot \varepsilon _{\rm{1}}}$的研究价值，表示不同动应力、围压、静偏应力下塑性应变率比值 ${\dot \varepsilon _N}$/ ${\dot \varepsilon _{\rm{1}}}$和加载次数N在双对数坐标系下的关系是统一的. 黄茂松等^[13]基于上海淤泥质软黏土在等向、偏压固结循环三轴加载下的试验结果，提出修正后的循环应变显式模型. 该模型基于塑性应变率的唯一性建立，可以同时考虑不同动偏应力水平、前期固结方式、围压的影响，经验证效果较好.

以上研究仅针对连续振动下的塑性应变发展，未考虑间歇对软土循环应变的影响. 已有非连续性振动加载^[14-16]仅模拟振动1 h停振1 h的长期振动-停振交替情况，不能模拟与地铁列车荷载类似的短期振动-停振交替情况，也尚未提出适用于非连续振动长期作用下的应变模型. 何绍衡等^[17-19]通过试验研究发现考虑间歇效应的塑性应变发展规律与连续振动条件下的相差较大，间歇会导致后续加载稳定时的最大塑性应变明显降低，使应变响应趋于稳定所需的振次减少，且停振时间与单位振动时长的比值越大，间歇效应越明显，但未建立考虑可停振比的应变模型，无法对间歇性循环加载下的塑性应变发展进行精确的经验预测.

本研究基于黄茂松^[12-13]的研究结论，针对连续振动下塑性应变增长率进行分析，验证塑性应变增长率的唯一性；在连续振动结论的基础上，研究考虑间歇的循环加载作用下塑性应变的经验模型.

试验土样为第4相海相沉积原状淤泥质软黏土，取自杭州拟建地铁2号线沿线，基本性质如表1所示. 表中，γ为重度，w为水的质量分数，G_s为比重，w_p为塑限，w_L为液限，I_p为塑性指数，I_L为液性指数. 取土深度均为4~6 m，取土区域的地下水位深度为1~2 m，因此土样处于饱和状态. 为了保证土样的原状性，采用I级取土器TB3型活塞式薄壁取土器，按照相关技术标准进行钻孔、贯入、取样、运输、密封、储存，所取土样符合《建筑工程地质勘探与取样技术规程》^[20]关于I级黏土的规定. 在制样时为了保证土样原状性，先把薄壁取土器中的土挤出，再将顶部和靠近薄壁边缘部分的土去掉，选取中心部分的土制成直径为38 mm、高为76 mm的标准圆柱体试件.

试验采用英国GDS公司的电伺服动三轴仪，动三轴试验步骤如下.

1）在制成标准圆柱形试样后，在试样上、下各贴一片滤纸，再各铺一块透水石，将试样放在饱和器中，进行一次饱和. 具体操作为：参照《土工试验规程》^[21]，将装有试样的饱和器放在真空桶中，用真空泵吸出空气维持负压保持3 h以上，再注入无气水浸泡12 h以上.

2）将试样从饱和器取出装入动三轴仪中进行二次饱和，原理为利用反压控制器将无气水强制注入试样，消融剩余空气. 当系统监测到无气水不能继续注入试样时，用sketpmon检测法检测土样饱和度b，当b=0.95时视为饱和完成. 本研究中试样均一次性达到饱和要求.

3）采用有效围压为90 kPa的应力路径对试样进行固结. 应力路径与K₀线重合，K₀线即固结阶段σ₃′/σ₁′=K₀的应力路径，其中，σ₃′为固结阶段施加的有效径向应力，σ₁′为固结阶段施加的有效轴向应力，K₀为土体的静止侧压力系数，经前期试验测得K₀=0.52.

4）进行间歇性循环加载. 如图1（a）所示为典型的地震荷载，如图1（b）所示为文献[22]给出的列车动荷载形式. 图中，L_N为归一化荷载，t为加载时长. 与地震荷载这类动荷载相比，列车运行引起的动荷载为单向脉冲式而非双向正弦模式^[23-25]. 根据国内列车运营普遍规律可知，列车全速通过某一点的时间约为6~18 s，相邻列车之间的发车间隙大于1 min（发车间隙通常是加载时长的几十倍）. 因此，地铁列车运营期荷载是由短暂的低幅振动和相对较长的停歇交替形成的周期性非连续荷载，并非连续振动. 用于间歇性循环加载来模拟地铁列车荷载的基础波形如图1（c）所示. 图中，T为单位振停循环的振动时长，ΔT为单位振停循环的停歇时长. 对于连续振动，选用基础的半正弦等幅连续波形. 考虑到列车荷载的实际情况^[23]，振动采用轴力控制式，振动频率为1 Hz，加载幅值为19~27 kPa，振动次数为10 000次.

如图2所示为循环三轴试验的有效应力路径. 图中，p′为平均有效主应力，p_s′为饱和阶段平均有效主应力，p₀′为振动阶段起始平均有效主应力，q为偏应力，q_s为振动阶段初始静偏应力，q_d为动偏应力，q_f为达到临界状态的偏应力（又称静破坏强度），M为临界状态线（critical state line，CSL）的斜率. 在三轴试验中，有

(1) ${p{'}} = {{\left( {\sigma _1{'} + 2\sigma _3{'}} \right)} /3},$

(2) $q = {\sigma _1} - {\sigma _3}.$

式中：σ₁、σ₃分别为轴向应力和径向应力. 如图2所示，饱和阶段的有效应力路径是一条与p′轴重合的直线. 为了保证试样合轴，设置q=0~1.0 kPa. 为了保证试样径向不受拉应力，对试样施加有效围压σ_3s′=10 kPa. 因此p′由0增至p_s′，由式（1）可知 ${p_{\rm{s}}}{'} \approx \sigma _{\rm{3s}}{'}$.

固结步的有效应力路径（即K₀线）是斜率小于临界状态线CSL斜率的直线. 在固结过程中p′、q成固定比例增加，p_s′增至p₀′，q从0增至q_s. 根据式（1）、（2）及K₀的定义，有

(3) ${q_{\rm{s}}} = \sigma _{\rm{1c}}{'} - \sigma _{\rm{3c}}{'} = \left( {{1 /{{K_0}}} - 1} \right)\sigma _{\rm{3c}}{'},$

(4) $p_0{'} = {{\left( {\sigma _{\rm{1c}}{'} + 2\sigma _{\rm{3c}}{'}} \right)} /3} + p_{\rm{s}}{'} = \left( {{1 /{{K_0}}} + 2} \right){{\sigma _{\rm{3c}}{'}} /3} + p_{\rm{s}}{'}.$

式中： $\sigma _{\rm{1c}}{'}$、 $\sigma _{\rm{3c}}{'}$分别为固结阶段施加的有效轴压和有效围压. 利用式（3）可以得到，振动初始静偏压q_s=83 kPa.

循环加载步的有效应力路径是众多斜率为3的短斜线圈，其纵向行程即振幅为q_d. 循环加载累积产生的残余孔隙水压力为u，u增大使平均有效主应力p′减小，导致短斜线逐渐向左移动，直至达到稳定或逼近临界状态线而破坏.

对有效应力路径发展进行分析可知，振动起始点到CSL的竖向距离为q_f−q_s，即振动幅值q_d的容许范围. 黄茂松^[12]基于此提出描述动应力水平的参量，相对偏应力水平：

(5) $D = {q_{\rm{d}}}/{\left( {q_{\rm{f}}} - {q_{\rm{s}}} \right)}.$

经前期试验测得q_f=146 kPa.

选用相对偏应力水平作为动应力水平的描述参量，可以同时考虑静偏应力、动偏应力、静破坏偏应力的耦合作用，具有较强的理论意义.

为了研究振动方式对塑性应变增长速率的影响，分别设置连续振动组和间歇振动组，试验方案如表2所示. 对于连续振动组仅考虑相对偏应力水平的影响，对于间歇振动组则同时考虑相对偏应力水平与停振比的影响.

如图3所示，在每次振动作用下土体产生的实时应变ε_i包括可恢复的弹性应变ε_e和不可恢复的塑性应变ε_p. 塑性应变曲线以振动阶段起始值为零点，由单位循环内应变的最低点连接可得. 由于本研究仅讨论塑性应变，将塑性应变ε_p直接记作ε.

如图4所示为在正常坐标系中，当相对偏应力D不同时，长期连续振动作用下的土体塑性应变发展曲线. 可以看出，在连续振动作用过程中塑性应变随振动发展单调递增，曲线斜率随振动发展逐渐减小；相对偏应力水平对塑性应变累积影响显著，在其余条件相同时，D越大，塑性应变累积值越大. 以上规律与已有研究成果一致^[3-13].

如图5所示为双对数坐标系下连续振动组的塑性应变发展曲线. 可以看出，在双对数坐标系下连续振动塑性应变累积曲线近似斜直线，振动次数越大，斜率越稳定，与郭林^[26]研究温州软黏土连续振动所得到的塑性应变发展规律趋势相似. 对连续振动下的等效初次应变和归一化应变增长率分别进行分析，结果如下.

1）等效初次塑性应变. 在GDS动三轴试验系统中，对某次循环下不可恢复应变的测量有不可避免的误差，原因如下：在振动过程中试样刚性底座上下移动，对于试样顶部的位移传感器读取轴向变形数据有一定的影响；在间歇加载下，对振动过程的采样精度不高. 塑性应变随振动次数增加而累积，测量偏差对累积塑性应变的影响减小，因此塑性应变累积曲线具有较高的可信度. 基于此，为了得到塑性应变在各振次的等效增量dε_N，采取数学微分方法^[27]：对应变累积曲线进行拟合，得到累积塑性应变与振次的关系函数ε（N）；求导得到ε^′（N），即各振次应变等效增量dε_N与振次N的关系函数；取N=1，可以得到等效初次循环塑性应变增长率dε₁，与等效初次塑性应变ε₁在数值上相等. 这种方法求得的等效初次应变在一定程度上可以表征真实的初次塑性应变.为了保证其准确度须尽量降低应变累积曲线的拟合误差，须选取相关系数高、趋势吻合度高的拟合函数. 分别采用幂函数y=ax^b、双曲函数y=x/（ax+b）这2种常见函数^[28]对试验结果进行拟合，效果如图6所示.可以看出，试验结果与双曲函数较吻合. 类似地，对各条塑性应变累积曲线进行双曲函数拟合，得到参数a、b及拟合情况，如表3所示. 表中，R为相关系数.

绘制连续振动作用下等效初次塑性应变ε₁与相对偏应力水平D的关系，如图7所示. 可以看出，等效初次塑性应变受相对偏应力水平的影响较显著，随着相对偏应力水平增加，等效初次塑性应变呈指数型增长，与黄茂松等^[12-13]所得的初次应变与相对偏应力水平的回归关系相似. 通过幂函数拟合上述关系，R²=0.99，拟合效果如图7所示.

2）归一化应变增长率. 对得到的拟合公式求导，即可得到等效塑性应变增长率dε_N与振次N的函数关系. 将各振次下的等效塑性应变增长率dε_N对初次增长率dε₁归一化得到归一化增长率dε_N/dε₁，物理意义为各振次等效应变增长率对初次增长率的衰减程度. 绘制归一化增长率dε_N/dε₁与振次N的关系，如图8所示. 可以看出，dε_N/dε₁的发展可以分为2个阶段：在N=1~1 000时处于快速衰减期，曲线近似为线性；在N=1 000~10 000时处于缓慢衰减期，曲线坡度不断减小. 在连续振动条件下，最终归一化等效塑性应变增长率约为12.8%，说明此时的塑性应变势仍然较为显著. 如果以dε_N/dε₁<10%作为应变增长进入稳定期的标志，则在连续振动下应变在振动过程中尚未进入稳定期，因为该淤泥质软土黏粒质量分数较高，土体结构单元亲水性较高，土体具有低渗透性、高压缩性的特点，在长期连续振动作用下塑性应变持续增长、不易稳定. 不同相对偏应力水平D的连续振动作用下得到的归一化增长率曲线重合度较高，说明D对连续振动归一化应变增长率几乎没有影响. 以双曲函数的一次导函数y=n/（mx+n）²拟合归一化增长率曲线，其中，m，n为曲线的形状参数. 拟合情况如表4所示. 可以看出，不同相对偏应力水平下的拟合结果较接近，但是与D没有相关性. 由此可知，当其余条件相同，仅D（等效于动应力幅值）不同时，连续振动下塑性应变的归一化增长率曲线具有唯一性，与黄茂松等^[12-13]对塑性应变的研究结论一致. 连续振动下曲线的形状参数m、n可以取平均值1.798×10⁻⁴、0.999 64.

在正常坐标系中，间歇性循环加载作用下的塑性应变累积曲线如图9所示. 可以看出，曲线有2处明显交叉：ΔT/T=10.0在N=2 000、7 000处，分别被ΔT/T=1.0、2.0反超；ΔT/T=5.0在N=1 800处被ΔT/T=1.0反超. 该现象表明不同停振比对应的应变增长率衰减速度不一样，无法平行增长，会出现交叉现象.

如图10所示为双对数坐标下的间歇振动塑性应变累积曲线. 可以看出，间歇振动作用下的结果与连续振动下的发展趋势在前期相似，均为近似线性规律，在后期出现明显差异. 在振动后期，间歇振动的塑性应变曲线逐渐偏离原有线性趋势、偏向横坐标轴，曲线坡度减小，其中ΔT/T=5.0、10.0时的偏离现象最为明显，在加载后期曲线近似水平，表征应变发展进入稳定期. 由此可知，在停振比较大时，应变稳定所需的振次较少，与K₀固结下间歇性循环加载结果的发展趋势一致^[17]. 下面对间歇振动下的等效初次应变和归一化应变增长率分别进行分析.

1）等效初次塑性应变. 间歇加载在双对数坐标系上出现明显弯曲，可知采用幂函数的拟合效果较差，因此选用y=x/（a+bx）对各曲线进行拟合，参数a、b及dε₁如表5所示.如图11所示为等效塑性应变ε₁与相对偏应力的关系. 可以看出，在间歇性循环加载下，与连续加载情况相比，ε₁随D增加而增大的趋势没有改变. 汇总间歇振动和连续振动试验结果，发现不同振动方式下的（ε₁，D）散点与同一条幂函数曲线较吻合，说明等效初次塑性应变ε₁与振动方式无明显相关性. 理论上，间歇不影响初次振动下的土体响应，该结果侧面反映了采用微分法拟合求解等效初次塑性应变的合理性.

2）归一化塑性应变率. 如图12所示为不同停振比ΔT/T对应的归一化应变增长率曲线. 如图12（a）所示为双对数坐标系中的试验结果. 可以看出，间歇加载下的归一化应变增长率普遍小于连续加载下的，且ΔT/T越大，差值越大. 当振动结束时，连续振动试样的残余归一化应变增长率约为13.1%，停振比ΔT/T=0.5、1.0、2.0、4.0、5.0、10.0时的残余归一化增长率分别为11.8%、8.2%、4.2%、1.8%、1.4%、0.6%. 可以看出，相同振动次数下残余归一化应变增长率随ΔT/T的增大单调减少. 在ΔT/T较大的情况下，dε_N/dε₁<10%，表明在间歇加载下应变增长比连续振动提早进入稳定期. 如图12（b）所示为正常坐标系中的试验结果. 可以看出，随停振比增大，塑性应变增长率的衰减趋势加剧，仅需较少的振次塑性应变增长率即可以进入稳定期，孔塑性应变累积曲线偏离原有增长趋势（见图9、10）表现为“屈服”形态.

由上述分析可知，连续振动下的塑性应变率仅随振动次数改变，间歇振动下的累积塑性应变率与振动次数、停振比均有关，为本研究模型的建立提供了重要依据.

相对偏应力水平仅影响等效初次塑性应变，而振动次数、停振比仅影响归一化塑性应变增长速率，因此可以建立间歇加载塑性应变模型. 由图10可知，当D、ΔT/T为常数时，等效初次塑性应变 ${\varepsilon _1}(D)$为常数，在后续振动过程中各振次下的塑性应变增量 ${\rm{d}}{\varepsilon _N}(D)$与归一化塑性应变增长率 ${{{\rm{d}}{\varepsilon _N}(D)}} / {{{\rm{d}}{\varepsilon _1}(D)}}$有如下关系：

(6) ${\rm{d}}{\varepsilon _N}(D) = \frac{{{\rm{d}}{\varepsilon _N}(D)}}{{{\rm{d}}{\varepsilon _1}(D)}}{\rm{d}}{\varepsilon _1}(D).$

${\rm{d}}{\varepsilon _1}(D)$为常数且与 ${\varepsilon _1}(D)$在数值上相等，因此式（6）可以写作：

(7) ${\rm{d}}{\varepsilon _N}(D) = \frac{{{\rm{d}}{\varepsilon _N}(D)}}{{{\rm{d}}{\varepsilon _1}(D)}}{\varepsilon _1}(D).$

由归一化塑性应变增长率的计算过程可知， ${{{\rm{d}}{\varepsilon _N}(D)}}/{{{\rm{d}}{\varepsilon _1}(D)}}$与振次N的关系与双曲函数导函数较吻合，因此以导函数 $y = {n}/{{{{(mx + n)}^2}}}$拟合各归一化塑性应变增长率曲线，拟合结果如表6所示. 据前文可知， m、n仅受ΔT/T影响，当D、ΔT/T为常数时，m、n是常数. 因此式（7）可以进一步写作：

(8) ${\rm{d}}{\varepsilon _N}(D) = \frac{n}{{{{(mN + n)}^2}}}{\varepsilon _1}(D).$

ε₁的表达式为

(9) ${\varepsilon _1} = 4.387 \times {10^{ - 6}}\times{{\rm{e}}^{12.141D}}.$

结合式（8）、（9），可以得到

(10) $\begin{split}{\rm{d}}{\varepsilon _N}(D) = & \frac{n}{{{{(mN + n)}^2}}}\times \\ & 4.387 \times {10^{ - 6}}\times{{\rm{e}}^{12.141D}}.\end{split}$

对式（10）中振动次数N进行积分可以得到累积塑性应变，表达式为

(11) $\begin{split} \varepsilon = & \int\limits_0^N {{\rm{d}}{\varepsilon _N}(D){\rm{d}}} N = \int\limits_0^N {\frac{n}{{(mN + n)}^2}}\times \\ & 4.387 \times {{10}^{ - 6}}\times{{\rm{e}}^{12.141D}} {\rm{d}}N. \end{split}$

m、n、D均为常数，并且当N=0时，ε=0，可以得到式（11）的计算结果：

(12) $\varepsilon = 4.387 \times {10^{ - 6}}\times{{\rm{e}}^{12.141D}}\frac{N}{{mN + n}}.$

当T固定不变时，参数m、n与ΔT/T的关系如图13所示. 可以看出，m、n与ΔT/T有较好的线性回归关系，以线性函数y=kx+c拟合，结果如下：

(13) $m = 1.06 \times {10^{ - 4}}\left( {{{\Delta T}/ T}} \right) + 2.055 \times {10^{ - 4}},$

(14) $n = - 2.03 \times {10^{ - 4}}\left( {{{\Delta T} / T}} \right) + 0.999\;59.$

结合表6和表4可知，拟合参数n的变化范围较小（0.997 6~0.999 6）. 为了保证模型的实用性和简便性，认为参数n可以设定为表4、6中n的平均值，即n=0.999 2.

根据式（12）和上述分析，可以计算不同停振比、相对偏应力水平下的间歇性循环加载长期应变发展，预测结果与试验结果对比如图14所示. 可以看出，预测结果与试验结果较吻合. 说明该模型同时考虑相对偏应力水平和停振比的影响，能较好地模拟不同停振比引起的塑性应变发展曲线交叉现象，可以解释间歇对循环应变发展的影响机制. 该模型由两部分组成，第1部分表示等效初次塑性应变，第2部分表示振动次数对塑性应变增长率的影响，相比传统的应变预测显式模型，具有物理含义清晰、操作简便、参数明确的优点. 鉴于本研究中主要采用双曲函数进行式（12）的推导，而在不同相对偏应力水平下塑性应变发展曲线的形态差异较大^[12-13]，结合本研究试验情况，该模型的建议适用范围为0.2<D<0.4.

（1）在双对数坐标系下，连续振动下的累积应变曲线近似为一条直线；在间歇振动作用下，累积应变曲线在前期近似为直线，在后期逐渐向横坐标偏移，坡度逐渐减小.

（2）在连续振动下归一化塑性应变增长率随振动次数增加而单调减小，与振次的对应关系具有唯一性，与相对偏应力水平无关；在间歇振动下归一化塑性应变增长率与振次的关系受停振比影响显著，随停振比增大，增长率衰减加剧，应变增长提前进入稳定期.

（3）等效初次塑性应变随相对偏应力水平的增大呈指数型增长，不同振动方式对等效初次塑性应变的影响不大，与理论上间歇不影响初次循环应变的结论相符，表明以微分方法取得的等效初次应变具有合理性.

（4）基于上述结论，可以将塑性应变模型分解为：等效初次塑性应变和应变归一化增长率两部分，并建立间歇加载塑性应变预测显式模型，该模型能同时考虑停振比、相对偏应力水平的影响，经验证预测结果与试验结果较吻合，对于在地铁列车荷载长期作用下饱和软黏土的塑性应变预测具有指导意义.

（5）结合本研究结论可知，忽略地铁荷载间歇与考虑间歇的试验结果有一定差距，因此在地铁沉降预测及稳定性评估方面，未来研究有必要考虑地铁荷载间歇及其他符合实际荷载的特点，以观测到更加符合现实工况的室内试验结果，从而得到更具实际应用意义的结论。