非饱和土蠕变力学特性试验及经验模型研究

魏建柄，刘卫斌

(陕西铁道工程勘察有限公司，陕西 西安 710043)

蠕变特性是土力学中的核心内容之一，与基坑工程、边坡工程等长期稳定性紧密关联[1-3]。非饱和土是一种包含固、液、气的三相土，相比饱和土而言，非饱和土蠕变特性更为复杂。在库岸边坡中，由于降雨入渗及库水位的变化，土体在饱和及非饱和状态之间转化，土体具备非饱和特性，其蠕变变形逐渐累积，对库岸边坡的长期稳定性造成潜在威胁[4-5]。现对于非饱和土蠕变特性的研究已有一定的进展，李滨锷[6]针对非饱和堆积土，进行压缩蠕变状态下的土-水特征曲线试验，以此建立非饱和土固-液-气三相耦合模型；李晓宁等[7]研究地基土的非饱和蠕变特性，发现含水率的增加会引起回弹变形和压缩蠕变变形，前者变形大于后者；王新刚等[8]进行了基质吸力控制条件下的非饱和黄土蠕变试验，研究表明基质吸力越小，蠕变曲线达到稳定所需的时间越长，且当应力水平较低时，黄土的变形较小，具有一定程度的可恢复性；李冬等[9]研究了应变率对非饱和土变形和强度影响，将屈服应力表示为应变率的函数，从而推导一维条件下非饱和土的增量方程，以此描述非饱和土随时间变化的长期变形行为。

目前关于土体蠕变模型的研究已取得一定成果，现有蠕变模型主要分为经验模型和元件模型，其中经验模型因其针对性、灵活性强而得到广泛应用。比较典型的经验模型有Singh-Mitchell模型和Mesri模型[10-11]。Singh-Mitchell模型的应力-应变关系采用指数函数来描述，而Mesri模型采用双曲线函数来描述，两者的应变-时间关系都采用幂函数来描述[12]。邹良超等[13]利用Morgan Mercer Florin函数描述应变-时间关系，黄海峰等[14]基于Log-Modified函数对Mesri模型进行改进，提出Log-Modified经验模型。但是大多数蠕变经验模型都是基于饱和土的蠕变变形特征而建，无法反映非饱和土体中基质吸力对蠕变特性的影响。鉴于此，本文以某库岸边坡中非饱和粉质黏土为研究对象，进行考虑基质吸力控制条件下的三轴压缩固结排水蠕变试验，在Mesri模型及Log-Modified模型的基础上，通过初始切线模量建立与吸力的关系，从而构建能反映基质吸力变化的蠕变经验模型，并对改进后的两个模型进行对比分析，证明改进后的Mesri模型能较好地描述非饱和粉质黏土的蠕变特性。

1 试验材料和方案

1.1 试验材料及试样制备

粉质黏土试样取自某库岸边坡，基本物理力学参数如表1所列。土样经过自然风干后碾散，过2 mm筛，为了便于试样成形，使用蒸馏水将其配制成含水率为20%的土样，用保鲜膜包裹土样后静置1 d以使水分扩散均匀，再利用削土器制备成为Φ60 mm×120 mm的重塑圆柱样。

表1 土样基本物理力学参数

1.2 试验步骤及方案

试验设备采用FSR-6型非饱和土三轴蠕变仪，进行基质吸力控制条件下的非饱和三轴压缩蠕变试验，在蠕变试验之前先开展三轴排水剪切试验确定排水剪切强度τf。围压σ3设置为100 kPa，基质吸力s分别为100，200，300和400 kPa，试验采用分级加载方式，每1级加载历时200 h以上。应力水平D从0.55开始，对应的偏应力差为0.55τf，每一级递增0.05直至破坏，破坏偏应力为(σ1-σ3)f，三轴排水剪切试验结果及蠕变试验加载方案见表2。

表2 蠕变试验加载方案

2 试验成果分析

图1为s=400 kPa时的分级加载蠕变曲线，利用Boltzmann线性叠加原理[15]处理后如图2所示。选取图2中1 ，26，51，76，101，126，151，176，201 h共9个时间节点的偏应力和应变绘制成等时应力-应变曲线，如图3所示。图4为不同基质吸力下的等时应力-应变曲线，由于该曲线较为密集，为了方便对比，只保留时间节点为1 h和201 h的曲线。

图1 s=400 kPa时的分级加载蠕变曲线

图2 Boltzmann线性叠加原理处理后的s=400 kPa分别加载蠕变曲线

图3 等时应力-应变曲线

图4 不同基质吸力的等时应力-应变曲线

根据图1、图2分析如下：土样在轴向加载瞬间表现有一定量的弹性瞬时应变，随着加载时间的累积，蠕变变形不断地增长，应力水平越高，蠕变变形越大。由图3、图4可看出，当加载时间超过51 h以后，等时应力-应变曲线具有一定非线性特征，曲线簇表现出逐渐偏向应变轴(横轴)的趋势。实际上该非饱和粉质黏土在不同应力水平和基质吸力下，其蠕变曲线具有较强的相似性，因此可采用相同的应力-应变-时间关系函数来描述其蠕变力学特性。

3 经验蠕变模型

目前广泛使用的经验蠕变模型主要有Singh-Mitchell模型和Mesri模型，文献[14]提出一种基于Mesri模型的Log-Modified经验模型，而Singh-Mitchell模型已被证明无法预测低应力水平下的蠕变变形，由于Mesri模型和Log-Modified模型未能考虑基质吸力对蠕变特性的影响，故本文分别改进Mesri模型和Log-Modified模型，并对改进后的两个模型对非饱和土蠕变特性的辨识能力进行分析。

3.1 改进的Mesri模型

3.1.1应力-应变-时间关系

总应变由瞬时应变和排水蠕应变两部分组成，若不考虑触变效应、固结比和老化等因素影响[13]，土体蠕变的应力-应变-时间关系可以为：

ε=f1(k)f2(t)

(1)

式中：ε——应变；

f1(k)，f2(t)——应力-应变关系函数和应变-时间关系函数。

3.1.2应力-应变关系

从图3、图4可看出，曲线形态近似双曲线，故引入Kondner提出的双曲线型应力-应变方程[13]，该等轴双曲线可写为：

(2)

式中：σ1，σ3——最大、最小主应力；

a，b——双曲线方程参数。

对式(2)进行微分，可得初始切线模量Eu为：

(3)

当ε→∞时，式(2)的极限值即为最终偏应力差(σ1-σ3)ult：

(4)

实际上应变不可能趋于无穷，一旦达到排水剪切强度τf便破坏。为了使该双曲线经过[εf，(σ1-σ3)f]，引入破坏比Rf[10-11]：

(5)

将式(3)、(5)同时代入式(2)可得：

(6)

式中：D——应力水平，D=(σ1-σ3)/(σ1-σ3)f。

为了方便求取参数，现对式(6)做如下变换：

(7)

(8)

(9)

3.1.3应变-时间关系

应变-时间函数可采用多种函数形式，包括双曲线函数、幂次函数和对数函数等[12-13]。Mesri模型采用幂次函数：

(10)

式中：t1——初始蠕变的参考时间，由于1h时土样蠕变不明显，为了便于计算，本文将t1取为1h；

ε1——t=t1时的初始蠕应变；

m——lnε-lnt曲线的斜率。

3.1.4改进的Mesri模型建立

在这一案例中，学生最初猜测“钉子板上围占2格的图形，可以围多少个”的时候，多数学生认为可以围2～3种，猜可以围4种或4种以上的学生不足20%。学生最初围出的图形大致如图1所示。

为了使非饱和土经验模型可以反映定量化含水率作用，将基质吸力作为独立变量体现到模型中，需建立应力-基质吸力-应变-时间关系模型。Janbu[16]研究发现初始切线模型和围压σ3在双对数坐标中线性相关，即初始切线模量Eu是σ3的幂函数。故此，Eu与基质吸力s之间可采用同样函数：

(11)

式中：pa——大气压(101.33 kPa)；

F，n——材料常数；

s——基质吸力。

将式(6)取t=t1代入式(10)可得：

(12)

式中：D1——当t=t1时的D值。

破坏偏应力(σ1-σ3)f、初始切线模量Eu和破坏比Rf的取值与t无关。

再将式(11)代入式(12)可得：

(13)

3.2 改进的Log-Modified模型

文献[14]中的Log-Modified模型，其应力-应变关系函数与Mesri模型一致，不同之处在于应变-时间关系采用一种三参数的新型幂函数来描述，其表达式为：

(14)

式中：p，q，r——模型参数。

将式(6)取t=t1代入式(14)可得：

(15)

再将式(11)代入式(15)可得：

(16)

式(16)即为本文考虑基质吸力的改进Log-Modified模型。

4 模型参数求解及验证

4.1 应力-应变关系参数求解

由于Mesri模型和Log-Modified模型应力-应变函数相同，故其应力-应变关系参数也相同。由式(9)可看出，ε/D与ε线性相关，Rf为斜率，(σ1-σ3)f/Eu为截距。选取1～201 h共9个时间节点的蠕变数据，以s=400 kPa为例，绘制不同时刻下ε/D-ε的关系曲线(图5)，为了便于观察，省略部分节点。从图5中可看出ε/D与ε有明显的线性相关性，并且(σ1-σ3)f/Eu和Rf与应力和时间无关，限于篇幅，仅给出在s=400 kPa下的参数取值，如表3所列。

图5 不同时刻下的ε/D-ε曲线

表3 应力-应变关系参数

由表3可知，Rf随时间的增长呈减小趋势，而(σ1-σ3)f/Eu随其增长都呈增大趋势。由于Mesri模型和Log-Modified模型假定模型参数与应力和时间无关，所以不同基质吸力下的(σ1-σ3)f/Eu和Rf取9个时间节点的平均值，如表4所示。

表4 应力-应变关系参数平均值

由表4可看出，(σ1-σ3)f/Eu平均值和Eu随着基质吸力s的增大而递增，Rf平均值随着基质吸力s的增大的递减。随着s的减小，库岸边坡内非饱和土初始切线模量不断减小，土体变软，说明随着库岸边坡内含水率的增加，土体蠕变变形更加显著。

4.2 参数F，n的确定

改进Mesri模型及Log-Modified模型皆是基于式(11)建立与吸力的关系，式(11)为Eu与参数F，n之间的幂次关系。因此改进后的2个模型的参数F，n也一致。通过表4中的Eu对s/pa进行拟合，如图6所示。由图6可知，曲线幂次拟合效果较好，R2达到0.9602，参数F取477.47，n取0.1882。

图6 基质吸力s与初始切线模量Eu关系

4.3 应变-时间关系参数求解

4.3.1改进的Mesri模型

式(10)中，lnε与lnt线性相关，m为其斜率。lnε-lnt曲线(以s=400 kPa为例)如图7所示。Mesri模型应变-时间关系参数如表5所列。

图7 lnε与lnt关系曲线

表5 不同应力水平下的lnε-lnt曲线斜率

由图7可看出，lnε-lnt线性相关较好，R2平均达到0.9742。由式(10)可看出，参数m与应力和时间无关，完全由lnε-lnt曲线的斜率决定，所以取不同应力水平下m的平均值作为该基质吸力条件下的Mesri模型参数，如表6所列。

将表2中的破坏偏应力(σ1-σ3)f、表4中的Rf平均值、表6中的参数m平均值以及参数F，n代入式(13)得：

表6 参数m平均值

(17)

式(17)即为本文非饱和粉质黏土在σ3=100 kPa、s=400 kPa下的改进Mesri蠕变经验模型。限于篇幅，仅给出s=400 kPa下的模型。

4.3.2改进的Log-Modified模型

将表2中的破坏偏应力(σ1-σ3)f、表4中的Rf平均值以及参数F，n代入式(16)，利用数学优化软件1stOpt基于BFGS算法和通用全局优化法求取每一级应力水平下的参数p，q和r，取其平均值，改进Log-Modified模型参数p，q和r如表7所列。

表7 参数p，q和r

将表4和表7中的模型参数以及参数F,n代入式(16)可得到如下经验型蠕变方程：

(18)

式(18)即为本文非饱和粉质黏土在σ3=100 kPa、s=400 kPa下的改进Log-Modified蠕变经验模型。限于篇幅，仅给出s=400 kPa下的模型。

4.4 模型验证

以s=300 kPa和400 kPa的试验数据为例，分别对改进的2个模型进行对比验证，如图8所示，为了便于观察，图中分别将改进Mesri模型及Log-Modified模型简写为M模型和L-M模型。

图8 试验值与理论值对比曲线

图8中改进Mesri模型辨识试验数据的R2平均值达到0.9912，改进Mesri模型的预测效果优良，能较好地反映非饱和粉质黏土的蠕变特性。而改进Log-Modified模型辨识试验数据的R2平均值为0.9634，对于衰减蠕变阶段的辨识能力较差，该段理论曲线总是低于试验曲线，且稳定蠕变阶段快结束的部分，其理论曲线高于试验曲线。

改进Mesri模型中应变-时间关系仅有1个未知参数m，且m可通过线性拟合求取。而改进Log-Modified模型中应变-时间关系包含3个未知参数，这3个参数通过数学软件1stOpt基于某种算法求取，相较前者模型具有一定的不确定性。

综合分析，改进Mesri模型的预测效果较好，能够较为准确地描述非饱和粉质黏土的蠕变特性。

5 结论

(1)本文进行了非饱和粉质黏土固结排水三轴压缩蠕变试验，该土样蠕变特性显著，在不同基质吸力和应力水平下，蠕变曲线较为相似，由此分别建立考虑基质吸力的两种经验模型来描述不同应力环境下的蠕变特性。

(2)在同一应力水平下，随着基质吸力的减小，非饱和土蠕变变形逐渐增大。库岸边坡工程中含水率变化对于非饱和土时效变形影响显著，含水率的增加加剧蠕变变形。

(3)本文模型通过初始切线模量与基质吸力之间的关系，分别建立了改进后的Mesri模型及Log-Modified模型。基于试验数据进行对比验证，改进Mesri模型能较好地描述非饱和土的蠕变特性。

由于本文研究成果基于非饱和粉质黏土蠕变特性，所以本文所建蠕变经验模型对于其它土体的适用性还待进一步研究。