湿陷性黄土的流变本构模型研究*

欧湘萍 闫志濠 李 塘 朱云升 吕乃芝 郭慧峰 鲁 聪

(武汉理工大学交通学院1) 武汉 430063) (浙江省交通规划设计研究院有限责任公司2) 杭州 310000)(内蒙古高等级公路建设开发有限责任公司3) 呼和浩特 010060)(湖北省交通规划设计院股份有限公司4) 武汉 430051)

0 引 言

近年来，随着西部地区高速公路的快速发展，大量黄土隧道建成并投入使用，由湿陷性黄土流变特性所引发的隧道等构造物的长期稳定性问题越来越受到人们的重视[1-2].由于黄土分布的地域性，不同区域的黄土的力学性质与结构特性会有所差异，因此，无论是应力、应变分析还是本构模型[3]的建立都不具有普适性，必须结合当地黄土的实际情况并通过试验数据分析得出[4-5].文中以兰州定远隧道的湿陷性黄土为研究对象，在黄土三轴蠕变试验的基础上，对该湿陷性黄土的本构模型进行了研究，以反映湿陷性黄土的流变特性.研究可为解决黄土地区因湿陷性黄土流变而造成的隧道稳定性问题提供借鉴[6].

1 经验本构模型

经验模型是能够反映试样流变曲线形状的数学表达式，它直接建立在对试验数据的分析上，具有简单明了、参数较少的优点，在工程实践中应用起来方便快捷.但经验模型也有其不足之处，譬如缺乏严谨理论依据的支持，受试验条件等因素影响较大，在土体类型和应力状态改变时，经验模型的适用性仍有待提高等等.

通过对黄土试样进行含水率试验，得到定远隧道湿陷性黄土的含水率范围为1.6%～11.7%，由此选取1.6%,5%,11.7%三个含水率下的湿陷性黄土进行经验本构模型拟合分析.

拟合优度记为R2，常作为反映本构模型拟合效果的表征指标，可以良好的表征模型的适用性.R2取值范围为0～1，拟合优度越大(越靠近1)，则表示所建模型对试验曲线的拟合程度越高，越能描述其土样的实际流变特性.拟合优度的计算公式为

(1)

式中：SSE为残差平方和;SST为总离差平方和.SSE，SST的计算公式如下.

(2)

1.1 Singh-Mitchell模型

在现有的经验模型中，由辛格和米切尔在1968年建立的Singh-Mitchell模型应用较为广泛，模型数学表达式为[7-8]

(3)

对式(3)进行换元，令B=At1/(1-m),λ=1-m，可将模型公式转化为

(4)

对式(4) 两端取自然对数，得

(5)

ln (ε-ε0)=λlnt+β

(6)

以1.6%，5%含水率的湿陷性黄土试样为例，绘出ln (ε-ε0)-lnt散点图并拟合,见图1.

图1 Singh-Mitchell模型的ln (ε-ε0)-ln t图

由图1可知，其应变差-时间双对数曲线用一次函数拟合时，拟合优度R2均在0.988以上，具有良好的线性关系，由此求得Singh-Mitchell应变-时间本构模型中的参数λ，β取值见表1.

表1 模型参数λ，β取值

获取参数λ，β后，即可根据以上由式(6)转化而来的应变-时间式(7)，求出各时间点对应的模型拟合应变，其与实际应变的对比见图2.

ε=ε0+eβtλ

(7)

图2 Singh-Mitchell模型拟合图

计算拟合应变与实际应变间拟合优度R2的值，见表2.

表2 拟合优度R2统计表

由上述相关图表可知，采用的Singh-Mitchell模型对湿陷性黄土试样曲线拟合优度值都在0.98以上，线形重合度良好，仅在最末端拟合值稍高于实际应变量.

当时间t与参考时间t1相等时，式(4)可化简为：

(8)

式(8)两端取自然对数，求得应力应变关系为

(9)

ln (ε-ε0)=aσ3+lnB

ε-ε0=Beaσ3

(10)

取参考时间t1为1 h，绘出1.6%，5%两种含水率下湿陷性黄土试样的ln (ε-ε0)-σ3散点图，拟合情况见图3.

图3 ln(ε-ε0)-σ3关系图

由图3中的拟合情况及拟合优度R2取值中可知，当试验时间t=t1=1 h=3 600 s时，本试验所用土样的应变差对数值与试验所用围压σ3有较好的对应关系，参数a，B值计算见表3.

表3 参数a，B计算值表

因此，1.6%，5%两种含水率下的湿陷性黄土试样在t=3 600 s时所对应的Singh-Mitchell拟合方程为

1.6%含水率土样：ε-ε0=0.010 8e-0.001 4σ3

5%含水率土样：ε-ε0=0.004 6e0.008 1σ3

(11)

上述所讨论的Singh-Mitchell模型拟合精度相对较高，且模型参数可与围压建立关联，具有一定的物理意义，但也存在着适用范围不广的问题.经过进一步计算发现，对本试验中11.7%含水率试样的试验曲线用Singh-Mitchell模型拟合时，其拟合优度较低.Singh-Mitchell经验本构模型无法完整描述本试验中1.6%，5%，11.7%三种含水率黄土试样的应变-时间关系，原因在于所取黄土在不同含水率时的应变-时间曲线会表现出较大差别：较低含水率(1.6%与5%)下，试样蠕变曲线的应变率趋向于均匀减小；而对11.7%的较高含水率试样，由于其结构强度偏低导致蠕变总量本就普遍高于较低含水率试样，因此可以看出Singh-Mitchell经验模型适用于拟合后期仍有较可观应变的末端倾斜曲线.

1.2 对数经验本构模型

为了建立能够同时反应本试验中较低含水率(1.6%，5%)与较高含水率(11.7%)的湿陷性黄土试样流变特性的经验本构模型，在分析得知模型拟合效果与流变后期曲线应变率有关后，决定以最常用的对数模型ε=A+Blnt为蓝本，增加一个可自由控制曲线末端斜率且与时间相关的Ct因子，所得模型为

ε=A+Blnt+Ct

(12)

本对数模型中有A，B，C三个参数，选择采用OrignPro 8.0分析软件进行拟合，得到的对数经验本构模型拟合结果见图4.

图4 对数经验模型拟合图

文中采用的对数经验模型对本试验中1.6%，5%，11.7%三种含水率试样的应变-时间曲线都能达到较高的拟合程度，且拟合优度处于0.94～0.97范围内，拟合效果比较稳定.由于加入了与时间相关的Ct单项式，末端的计算应变值依旧普遍略大于试验应变且仍保持有稳定的增加趋势，这与实际流变试验中试样应变最终趋于稳定的现象有所出入.

对数经验模型参数与拟合优度R2见表4.

表4 对数经验模型参数及拟合优度R2表

由表4可知，参数A取值与试样初始应变ε0大小相近，变化规律也与之相同：较低围压(100 kPa)下A随含水率的增大而变小，较高围压(300 kPa)下则相反，可知本模型中A值代表的含义等同于初始应变.参数B作为与对数时间相乘的因子，在100，200和300 kPa三种围压下都随含水率的增大而增大；而表征曲线末端第二流变阶段常应变速率的参数C则与围压(偏应力)的变化有明显关联，对同一含水率试样，流变趋向于稳定前的应变率随偏应力的增大而增大，见图5.

图5 参数B，C变化规律

2 Burgers元件本构模型

元件模型法是将Hooke弹性体、Newton黏壶及St. Venant塑性体等分别反映弹性、黏滞性、塑性等不同力学特性的元件按照土体实际性质组合成力学模型，然后根据试验数据求出组合模型各项参数的方法.

首先结合卸载回弹阶段的应变情况，可知所需元件模型应至少包含造成初始应变的串联弹性体与一个可产生弹性后效的开尔文体，此时即构成三单元的广义开尔文模型，但该模型的曲线末端趋于平缓，仅可反映第一阶段蠕变趋于稳定的情况.从建立对数模型的经验中可知，文中所得试验曲线的末端普遍较为倾斜，应处于具有常应变速率的第二蠕变阶段，与增加可自由控制曲线末端斜率且与时间相关的Ct项类似，因此，决定在广义开尔文模型的基础上，增加一个反映常应变速率的串联粘壶元件，最终形成了反映土样粘弹性的四参数Burgers(伯格斯)模型，其结构见图6.

图6 Burgers模型结构

由图6可知，所选用的Burgers元件模型由四个单元构成：虎克弹簧G1与牛顿黏壶η1串联形成Maxwell(麦克斯韦)模型，弹性体G2和黏壶η2并联构成Kelvin(开尔文)模型，最后再将Maxwell体与Kelvin体串联，最终形成Burgers模型.

本文中的湿陷性黄土流变试验是在三向应力状态下进行，对Burgers模型中的弹性体变形情况应采用剪切模量G与体积模量K来表示，模型变形计算公式为

(13)

式中：σⅠ，σⅡ与流变的三向应力大小有关，其公式为

σⅠ=σ1-σ3

σⅡ=σ1+2σ3

(14)

式中：剪切模量G体现的是三向加压状态下由于偏应力引起虎克弹性体形状改变的性质，体积模量K则体现了在三向等压状态下的体变特性，其方程为

(15)

(16)

式中：E为弹性模量；μ为土体泊松比，本次试验所取湿陷性黄土泊松比参考值μ为0.44，再通过式(15)～(16)可将Burgers模型中需拟合的未知量转化为E1，E2，η1，η2四个参数，它们分别代表着本构关系中两个弹性体的弹性模量与两个粘壶的粘滞系数.

2.1 模型参数计算

对Maxwell体中串联弹性体参数E1，令试验时间t=0，综合式(13)～(16)可得

式中：t=0时的初始应变ε0可从试验中获得，其余参数σ1，σ3及μ值均已知，则可求出E1.

对串联粘壶系数η1，从伯格斯本构模型式(13)可知，σⅠ/3η1即为流变曲线斜率.对所得应变-时间曲线，当t较大时，将曲线末端视为直线并通过线性回归求出斜率值k，则可通过式η1=σⅠ/3k计算出η1.

求取Kelvin体中的模型参数E2，η2时，需在Burgers本构方程(13)中令t趋于无穷大，则有

(18)

用式(18)减去式(13)可得q值：

(19)

q值为蠕变末端渐近线与实际应变之差，对式(19)两端取对数，可得

(20)

绘出lnq-t关系散点图，按照线性回归对其添加趋势线，求得拟合直线的截距和斜率，再结合式(14)～(15)分别求出E2,η2的数值.

2.2 模型拟合结果分析

按照2.1的三个步骤可计算出不同含水率(1.6%,5%,11.7%)和不同围压(100,200和300 kPa)下湿陷性黄土试样的伯格斯参数值，将其代入伯格斯模型式(13)中，可得各试验条件下的应变-时间函数，从而求出各时间点的拟合应变与实际应变值，见图7.

图7 Burgers元件模型拟合图

由图7可知，由于Burgers模型的力学特性决定了拟合曲线末端是具有常应变速率的的第二流变阶段，而实际试验中应变时间曲线后期则是趋于平缓，应变速率逐渐减小,因此,拟合曲线的尾部均有稍高于实际应变量的现象.但从整体拟合情况上看，伯格斯模型拟合应变与实际应变间的拟合精度R2值都大于0.96，见表5.可以较好地描述本试验中所有湿陷性黄土试样的流变性质.

3 结 论

1) 在两种经验本构模型中，Singh-Mitchell模型对两种含水率下(1.6%,5%)的湿陷性黄土流变试验曲线的拟合优度R2值均在0.98以上，在两种经验本构模型中拟合优度较高，但适用范围有限，无法反映试验中较高含水率下(11.7%)试样的流变特性；经过修改的对数经验模型，其拟合优度R2在两种经验模型中处于中等水平，但可适用于本试验中三种含水率下(1.6%,5%,11.7%)的湿陷性黄土试样全部应变曲线，是拟合效果更好的经验本构模型.

2) Burgers元件模型对湿陷性黄土流变试验曲线的拟合优度R2高于对数经验本构模型，也可适用于不同含水率(1.6%,5%,11.7%)与不同围压(100,200和300 kPa)下的湿陷性黄土试样的流变特性分析，且各元件参数意义较为明确.因此，在对定远隧道进行稳定性分析与变形预测中，Burgers元件模型比对数经验本构模型更为合适.

表5 Burgers元件模型参数表