土体变形具有明显的流变性质。为了理论和数值模拟土体的流变特性，岩土学者已建立了许多土体流变模型。但这些模型主要侧重于研究饱和土的应变速率效应^[1-6]，有关非饱和土的流变模型则相对较少^[7-8]。Gennaro等^[7]率先开展非饱和土率相关本构模型的研究工作.他们把净应力和吸力作为状态变量^[9]，根据非饱和土等应变速率压缩线提出了屈服净应力与应变速率之间的关系式，建立了非饱和土弹黏塑性模型.李冬等^[8]把土骨架平均应力和有效吸力作为状态变量^[10-11]，建立了非饱和土弹黏塑性模型.他们在非饱和土流变模型上的开创性工作，有力地推动了非饱和土流变力学的发展.但这些非饱和土流变模型在描述速率效应对屈服应力的影响时，认为总体应变速率(等于弹性应变速率和黏塑性应变速率之和)影响非饱和土的屈服应力.而现代连续介质力学认为^{[3, 12-13]}，弹性变形属于可恢复变形，其速率不会对引起不可恢复变形的屈服应力产生影响.因此采用黏塑性应变速率而不是总应变速率来建立屈服应力流变方程，更符合热力学规律.其次，非饱和土具有两种塑性机制.一个是净应力加载引起的塑性体积变形，采用LC屈服面来描述；另一个是吸力加载引起的塑性体积变形，采用SI屈服面来反映其变形特征.但现有的非饱和土流变模型只考虑LC屈服面的变形特性，对非饱和土流变特性的研究还欠全面.因此有关非饱和土流变模型的研究还有待日臻完善.

本文借鉴巴塞罗那非饱和土弹塑性模型^[9]，以净应力和吸力作为建模应力变量，根据等效时间法^[4]，分别获得LC和SI两个屈服面的弹黏塑性方程，建立了非饱和土双屈服面流变模型.“净应力和吸力”与“土骨架平均应力和有效吸力”成一一对应关系，故“净应力和吸力”是与“土骨架平均应力和有效吸力”等效的非饱和土本构模型应力变量，采用“净应力和吸力”双变量来建模同样具有热力学意义上的合理性.

1 非饱和土流变模型

非饱和土弹塑性模型存在两个屈服面，一个是LC屈服面，另一个是SI屈服面.根据弹黏塑性流变模型理论，双屈服面的总应变可以表示为弹性应变和两个屈服面的黏塑性应变之和：

$ {{\dot \varepsilon }_{\rm{v}}} = \dot \varepsilon _{\rm{v}}^{\rm{e}} + \dot \varepsilon _{\rm{v}}^{{\rm{Lvp}}} + \dot \varepsilon _{\rm{v}}^{{\rm{svp}}}, $

(1)

$ {{\dot \zeta }_{ij}} = \dot \zeta _{ij}^{\rm{e}} + \dot \zeta _{ij}^{{\rm{Lvp}}} + \dot \zeta _{ij}^{{\rm{svp}}}. $

(2)

式中：ε_v为体应变，ε_v^e为弹性体应变，ε_v^Lvp和ε_v^svp分别是LC和SI屈服面的黏塑性体应变.式(2)中的ζ_ij为偏应变，ζ_ij^e为弹性偏应变，ζ_ij^Lvp和ζ_ij^svp分别是LC和SI屈服面的黏塑性偏应变.

1.1 弹性方程

非饱和土的弹性本构方程可表示为^[9]：

$ \varepsilon _{\rm{v}}^{\rm{e}} = \frac{\kappa }{{{v_0}}}\frac{{\dot {\bar p}}}{{\bar p}} + \frac{{{\kappa _{\rm{s}}}}}{{{v_0}}}\frac{{\dot s}}{{s + {p_{{\rm{atm}}}}}}, $

(3)

$ \dot \zeta _{ij}^{\rm{e}} = {{\dot S}_{ij}}/\left( {2G} \right). $

(4)

式中：p=p-p_G为净应力，s=p_G-p_L为吸力，p、p_G和p_L分别是总应力，孔隙气体和液体压力，p_atm为标准大气压，κ和κ_s为净应力和吸力回弹模量.υ₀=1+e₀为初始比容，e₀为初始孔隙比.S_ij为偏应力，G=[1.5(1-2μ)/(1+μ)](p′υ₀/κ)为剪切模量，μ为泊松比.

1.2 LC流变屈服方程

根据巴塞罗那弹塑性模型，非饱和土的LC屈服方程为

$ F=q^{2}+M^{2}\left(\bar{p}+\bar{p}_{\mathrm{s}}\right)\left(\bar{p}-\bar{p}_{\mathrm{c}}\right)=0. $

(5)

式中：q为广义剪应力，p_c为屈服净应力，p_s=k_as，k_a为常数^[9].在塑性模型中，屈服净应力的硬化参数是塑性体应变和吸力；而在黏塑性流变模型中，屈服净应力的硬化参数是黏塑性体应变速率、黏塑性体应变和吸力，即屈服净应力的硬化规律可表示为p_c=f₁($\dot \varepsilon $_v^Lvp, ε_v^Lvp, s).假定p_c与$\dot \varepsilon $_v^Lvp成一一对应关系，对f₁求$\dot \varepsilon $_v^Lvp的反函数得

$ \dot \varepsilon _{\rm{v}}^{{\rm{Lvp}}} = f\left( {{{\bar p}_{\rm{c}}},s,\varepsilon _{\rm{v}}^{{\rm{Lvp}}}} \right). $

(6)

式(6)表明黏塑性体应变速率可唯一地表示为屈服净应力，吸力和黏塑性体应变的函数.

非饱和土试验研究表明，当吸力保持不变时，不同等效时间的黏塑性体应变随屈服净应力的压缩线近似为平行线且随等效时间对数增长^[7]，见图 1.

选择t_r时刻的等时间压缩线作为参考时间线，t_Le为常应力和吸力作用下从参考时间起算的变化到现时黏塑性体应变所需的蠕变时间^[4]，则非饱和土的黏塑性应变公式可表示为

$ \varepsilon _{\rm{v}}^{{\rm{Lvp}}} = \varepsilon _{\rm{r}}^{{\rm{Lvp}}}\left( s \right) + \frac{{\lambda \left( s \right) - \kappa }}{{{v_0}}}\ln \frac{{{{\bar p}_{\rm{c}}}}}{{{{\bar p}_{\rm{r}}}}} + \frac{{\psi \left( s \right)}}{{{v_0}}}\ln \frac{{{t_{\rm{r}}} + {t_{{\rm{Le}}}}}}{{{t_{\rm{r}}}}}. $

(7)

式中：λ(s)为压缩指数，ψ(s)为次固结系数，p_r为参考净应力，υ₀为初始比容，ε_r^Lvp(s)为吸力s时在参考时间t_r线上参考净应力p_r所对应的参考黏塑性体应变，本文把ε_r^Lvp(s)称为参考时间黏塑性体应变.对式(7)求时间导数得

$ \dot \varepsilon _{\rm{v}}^{{\rm{Lvp}}} = \psi \left( s \right)/\left[ {{v_0}\left( {{t_{\rm{r}}} + {t_{{\rm{Le}}}}} \right)} \right]. $

(8)

令黏塑性指数比β(s)为

$ \beta \left( s \right) = \psi \left( s \right)/\left[ {\lambda \left( s \right) - \kappa } \right]. $

(9)

式(8)中的t_Le用式(7)代替后利用式(9)得

$ \dot \varepsilon _{\rm{v}}^{{\rm{Lvp}}} = \frac{{\psi \left( s \right)}}{{{v_0}{t_{\rm{r}}}}}\exp \left[ {\frac{{\varepsilon _{\rm{r}}^{{\rm{Lvp}}}\left( s \right) - \varepsilon _{\rm{v}}^{{\rm{Lvp}}}}}{{\psi \left( s \right)/{v_0}}}} \right]{\left( {\frac{{{{\bar p}_{\rm{c}}}}}{{{{\bar p}_{\rm{r}}}}}} \right)^{\frac{1}{{\beta \left( s \right)}}}}. $

(10)

由于黏塑性体应变速率是关于屈服净应力，吸力和黏塑性体应变的唯一函数.故式(10)必须与式(6)相等.同时，式(10)中的黏塑性体应变速率是通过蠕变时间t_Le间接获得的，与文献[4]一样，把t_Le称为等效时间.

从式(8)可知，当黏塑性体应变速率相同时，等效时间t_Le会随着参考时间t_r选择的不同而发生变化，这不利于流变模型参数在不同参考时间之间的转换.故引入绝对等效时间t_La^{[3, 13]}，它的定义为参考时间t_r与等效时间L_Le之和.由式(8)可知，绝对等效时间与黏塑性体应变速率成一一对应关系.由于黏塑性体应变速率是一个客观量，故绝对等效时间也是一个客观量，它不会随着参考时间的不同选择而发生变化，因此工程应用时绝对等效时间比一般等效时间要方便得多.把t_La=t_Le+t_r代入式(7)可得黏塑性体应变-净应力-绝对等效时间的唯一性关系方程为

$ \varepsilon _{\rm{v}}^{{\rm{Lvp}}} = \varepsilon _{\rm{r}}^{{\rm{Lvp}}}\left( s \right) + \frac{{\lambda \left( s \right) - \kappa }}{{{v_0}}}\ln \frac{{{{\bar p}_{\rm{c}}}}}{{{{\bar p}_{\rm{r}}}}} + \frac{{\psi \left( s \right)}}{{{v_0}}}\ln \frac{{{t_{{\rm{La}}}}}}{{{t_{\rm{r}}}}}. $

(11)

假定LC屈服准则符合相关联流动法则，根据Perzyna过应力势函数理论，由式(5)得

$ \left( {\dot \varepsilon _{\rm{v}}^{{\rm{Lvp}}}/\dot \varepsilon _{\rm{s}}^{{\rm{Lvp}}}} \right) = {M^2}\left( {2\bar p + {{\bar p}_s} - {{\bar p}_{\rm{c}}}} \right)/\left( {2q} \right), $

(12)

式中ε_s^Lvp为LC屈服面的广义剪应变.由式(10)、(12)得

$ \begin{array}{l} \dot \varepsilon _{\rm{s}}^{{\rm{Lvp}}} = \frac{{2q\psi \left( s \right)/\left( {{v_0}{t_{\rm{r}}}} \right)}}{{{M^2}\left( {2\bar p + {{\bar p}_s} - {{\bar p}_{\rm{c}}}} \right)}} \times \\ \;\;\;\;\;\;\exp \left[ {\frac{{\varepsilon _{\rm{r}}^{{\rm{Lvp}}}\left( s \right) - \varepsilon _{\rm{v}}^{{\rm{Lvp}}}}}{{\psi \left( s \right)/{v_0}}}} \right]{\left( {\frac{{{{\bar p}_{\rm{c}}}}}{{{{\bar p}_{\rm{r}}}}}} \right)^{\frac{1}{{\beta \left( s \right)}}}}. \end{array} $

(13)

当q=0和净应力p为压时，由式(5)可知p_c=p.把p_c=p代入到式(10)和式(11)可得：

$ \dot \varepsilon _{\rm{v}}^{{\rm{Lvp}}} = \frac{{\psi \left( s \right)}}{{{t_{\rm{r}}}}}\exp \left[ {\frac{{\varepsilon _{\rm{r}}^{{\rm{Lvp}}}\left( s \right) - \varepsilon _{\rm{v}}^{{\rm{Lvp}}}}}{{\psi \left( s \right)}}} \right]{\left( {\frac{{\bar p}}{{{{\bar p}_{\rm{r}}}}}} \right)^{\frac{1}{{\beta \left( s \right)}}}}, $

(14)

$ \varepsilon _{\rm{v}}^{{\rm{Lvp}}} = \varepsilon _{\rm{r}}^{{\rm{Lvp}}}\left( s \right) + \frac{{\lambda \left( s \right) - \kappa }}{{{v_0}}}\ln \frac{{\bar p}}{{{{\bar p}_{\rm{r}}}}} + \frac{{\psi \left( s \right)}}{{{v_0}}}\ln \frac{{{t_{{\rm{La}}}}}}{{{t_{\rm{r}}}}}. $

(15)

式(14)是等向压缩条件下黏塑性体应变速率公式，式(15)是等向压缩条件下吸力s恒定时黏塑性体应变-净应力-绝对等效时间的唯一性关系方程.鉴于一维完全侧限压缩变形和等向压缩变形的力学机理均为体积压缩变形，因此把式(14)~(15)中的p和p_r改为竖向净应力后，可以描述一维完全侧限条件下非饱和土的流变特性.

1.3 LC屈服净应力的硬化规律方程

从式(5)、(6)之间的文字说明可看出，黏塑性体应变速率方程和屈服净应力硬化规律方程互为反函数.因此，本节通过LC流变速率方程建立包含黏塑性体应变速率影响的屈服净应力硬化规律.

1) 首先推导黏塑性体应变-净应力-黏塑性体应变速率之间的唯一性关系方程

选择$\dot \varepsilon $_vr^Lvp作为参考黏塑性体应变速率，ε_vr^Lvp(s)是吸力s时p_r在参考$\dot \varepsilon $_vr^Lvp线上所对应的参考黏塑性体应变，本文把ε_vr^Lvp(s)称为参考速率黏塑性体应变.取$\dot \varepsilon $_v^Lvp=$\dot \varepsilon $_vr^Lvp、ε_v^Lvp=ε_vr^Lvp(s)和p_c=p_r代入式(10)得

$ \varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}} = \varepsilon _{\rm{r}}^{{\rm{Lvp}}}\left( s \right) + \frac{{\psi \left( s \right)}}{{{v_0}}}\ln \frac{{\psi \left( s \right)}}{{{v_0}{t_{\rm{r}}}\varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}}}. $

(16)

利用式(16)把式(10)的t_r用$\dot \varepsilon $_vr^Lvp代入后变换得

$ \varepsilon _{\rm{v}}^{{\rm{Lvp}}} = \varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}\left( s \right) + \frac{{\lambda \left( s \right) - \kappa }}{{{v_0}}}\ln \frac{{{{\bar p}_{\rm{c}}}}}{{{{\bar p}_{\rm{r}}}}} + \frac{{\psi \left( s \right)}}{{{v_0}}}\ln \frac{{\dot \varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}}}{{\dot \varepsilon _{\rm{v}}^{{\rm{Lvp}}}}}. $

(17)

式(17)为黏塑性体应变-净应力-黏塑性体应变速率之间的唯一性关系方程.

2) 其次推导同一屈服面上屈服净应力与吸力和黏塑性体应变速率之间的数学关系式

饱和土的s=0，净应力等于有效应力.设饱和土参考$\dot \varepsilon $_vr^Lvp线上黏塑性体应变ε_v^Lvp所对应的屈服净应力为p_{c, r}(0).LC屈服面以黏塑性体应变为硬化参数，故同一个LC屈服面上的黏塑性体应变相等.根据式(17)中的ε_v^Lvp保持不变可得

$ \frac{{{{\bar p}_{\rm{c}}}}}{{{{\bar p}_{\rm{r}}}}} = \exp \left\{ {\frac{{\varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}\left( 0 \right) - \varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}\left( s \right)}}{{\left[ {\lambda \left( s \right) - \kappa } \right]/{v_0}}}} \right\}{\left( {\frac{{\varepsilon _{\rm{v}}^{{\rm{Lvp}}}}}{{\dot \varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}}}} \right)^{\beta \left( s \right)}}{\left[ {\frac{{{{\bar p}_{{\rm{c}},{\rm{r}}}}\left( 0 \right)}}{{{{\bar p}_{\rm{r}}}}}} \right]^{\frac{{\lambda \left( 0 \right) - \kappa }}{{\lambda \left( s \right) - \kappa }}}}. $

(18)

式(18)给出了同一屈服面上s和$\dot \varepsilon $_v^Lvp时的p_c与饱和土$\dot \varepsilon $_vr^Lvp时的p_{c, r}(0)之间的关系式.若同一屈服面s相同时参考$\dot \varepsilon $_vr^Lvp线上的屈服净应力为p_c=p_{c, r}(s)，代入式(18)后再把式(18)原式中的p_r用p_{c, r}(s)代替得

$ {{\bar p}_{\rm{c}}}/{{\bar p}_{{\rm{c}},{\rm{r}}}}\left( s \right) = {\left( {\dot \varepsilon _{\rm{v}}^{{\rm{Lvp}}}/\dot \varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}} \right)^{\beta \left( s \right)}}. $

(19)

式(19)是同一屈服面吸力相同时p_c随$\dot \varepsilon $_v^Lvp的变化公式.显然，当s=0时式(19)与Yin-Graham获得的饱和土屈服应力随黏塑性体应变速率变化公式相同^[14]，这说明本文获得的非饱和土屈服净应力公式可退化为饱和土屈服净应力公式.

3) 最后推导后继屈服面硬化方程

非饱和土LC后继屈服面反映的是屈服净应力随黏塑性体应变发展关系.设饱和土初始黏塑性体应变ε_v0^Lvp=0，在参考$\dot \varepsilon $_vr^Lvp线上对应的初始屈服净应力为p_c0，把p_c=p_c0和ε_v^Lvp=ε_v0^Lvp=0代入到式(17)后再去减式(17)原式可得

$ \begin{array}{*{20}{c}} {\frac{{{{\bar p}_{\rm{c}}}}}{{{{\bar p}_{\rm{r}}}}} = \exp \left\{ {\frac{{\varepsilon _{\rm{v}}^{{\rm{Lvp}}} + \varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}\left( 0 \right) - \varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}\left( s \right)}}{{\left[ {\lambda \left( s \right) - \kappa } \right]/{v_0}}}} \right\} \times }\\ {{{\left( {\frac{{\dot \varepsilon _{\rm{v}}^{{\rm{Lvp}}}}}{{\dot \varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}}}} \right)}^{\beta \left( s \right)}}{{\left( {\frac{{{{\bar p}_{{\rm{c}}0}}}}{{{{\bar p}_{\rm{r}}}}}} \right)}^{\frac{{\lambda \left( 0 \right) - \kappa }}{{\lambda \left( s \right) - \kappa }}}}.} \end{array} $

(20)

式(20)即是LC后继屈服面的屈服净应力方程.

1.4 等效时间LC屈服流变方程的特色

1) LC屈服面硬化方程的退化

假定任何s的ε_vr^Lvp(s)相同，即：

$ \varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}\left( 0 \right) = \varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}\left( s \right). $

(21)

式(21)表明在ε_v^Lvp~p_c图中所有s的等$\dot \varepsilon $_vr^Lvp压缩线穿过同一个点.当不考虑非饱和土蠕变特性时β(s)=0.把式(21)和β(s)=0代入式(18)得

$ \frac{{{{\bar p}_{\rm{c}}}}}{{{{\bar p}_{\rm{r}}}}} = {\left[ {\frac{{{{\bar p}_{{\rm{c}},{\rm{r}}}}\left( 0 \right)}}{{{{\bar p}_{\rm{r}}}}}} \right]^{\frac{{\lambda \left( 0 \right) - \kappa }}{{\lambda \left( s \right) - \kappa }}}}. $

(22)

式(22)与非饱和土巴塞罗那弹塑性模型^[9]的屈服方程相同，这表明本文流变屈服方程可以退化为巴塞罗那LC屈服净应力方程.

2) 与Gennaro-Pereira非饱和土流变模型的比较

Gennaro-Pereira^[7]通过等速率法建立了非饱和土流变模型，获得的屈服方程为

$ \frac{{{{\bar p}_{\rm{c}}}}}{{{{\bar p}_{\rm{r}}}}} = \exp \left[ {\frac{{{v_0}\varepsilon _{\rm{v}}^{{\rm{Lvp}}}}}{{\lambda \left( s \right) - \kappa }}} \right]{\left[ {\frac{{\dot \varepsilon _{\rm{v}}^{\rm{e}} + \dot \varepsilon _{\rm{v}}^{{\rm{Lvp}}}}}{{{{\left( {\dot \varepsilon _{\rm{v}}^{\rm{e}} + \dot \varepsilon _{\rm{v}}^{{\rm{Lvp}}}} \right)}_{{\rm{vr}}}}}}} \right]^{\frac{{\lambda \left( 0 \right) - \mathit{\boldsymbol{\kappa }}}}{{\lambda \left( s \right)}}\beta \left( s \right)}}{\left( {\frac{{{{\bar p}_{{\rm{c}}0}}}}{{{{\bar p}_{\rm{r}}}}}} \right)^{\frac{{\lambda \left( 0 \right) - \kappa }}{{\lambda \left( s \right) - \kappa }}}}. $

(23)

比较式(20)和式(23)可以发现本文LC屈服面的流变屈服方程和Gennaro-Pereira流变屈服方程之间存在明显差异.(Ⅰ)本文屈服方程式(20)采用黏塑性体应变速率表示，而Gennaro-Pereira屈服方程式(23)采用总体应变速率表示；(Ⅱ)式(20)中黏塑性体应变速率的幂是β(s)，而式(23)中总的体应变速率的幂是β(s){[λ(0)-κ]/λ(s)}.由于饱和土的Yin-Graham屈服净应力只与黏塑性体应变速率有关而与弹性体应变速率无关^[13-14]，故式(23)无法退化到饱和土的Yin-Graham屈服净应力公式.

3) 满足式(21)的$\dot \varepsilon $_vr^Lvp和p_r的选择难易探讨

如果式(20)中ε_vvr^Lvp(s)和ε_vvr^Lvp(0)满足式(21)，则LC后继屈服面的净应力方程式会更简单.现证明：最多只有一组($\dot \varepsilon $_vr^Lvp，p_r)对于所有的s满足式(21).

实际上，假定还存在另一组参考黏塑性体应变速率$\dot \varepsilon $_vr^Lvp*和参考净应力p_r^*，对于所有的s也满足式(21).根据($\dot \varepsilon $_vr^Lvp，p_r)和($\dot \varepsilon $_vr^Lvp*，p_r^*)同时满足式(17)和式(21)后取s=s和s=0，可得

$ \bar p_{\rm{r}}^* = {{\bar p}_{\rm{r}}}{\left( {\dot \varepsilon _{{\rm{vr}}}^{{\rm{Lv}}{{\rm{p}}^*}}/\dot \varepsilon _{{\rm{vr}}}^{{\rm{Lvp}}}} \right)^{\frac{{\psi \left( 0 \right) - \psi \left( s \right)}}{{\lambda \left( 0 \right) - \lambda \left( s \right)}}}}. $

(24)

土工试验表明ψ(s)/λ(s)随s发生变化^[7]，故式(24)等式右边也随s发生变化，但等式左边的p_r^*对所有的s取同一个值，故式(24)不成立.这就证明了满足式(21)的参考黏塑性体应变速率和参考净应力最多只有一组.故根据土工试验确定唯一满足式(21)的参考黏塑性体应变速率和参考净应力是困难的.在实际应用中只能人为地假定参考黏塑性体应变速率和参考净应力满足式(21)，从而不可避免地会带来误差.式(20)不需要满足苛刻条件式(21)，因此采用式(20)作为非饱和土的LC流变屈服方程更合理.

1.5 SI流变屈服方程

非饱和土除LC屈服面外还存在SI屈服面.SI屈服面的屈服条件为s=s_c=constant^[9].假定SI屈服面的压缩变形随蠕变时间指数衰减，其表达式为

$ \varepsilon _{\rm{v}}^{{\rm{svp}}} = \varepsilon _{v0}^{{\rm{svp}}} + \frac{{\lambda _{\rm{s}}^{{\rm{vp}}}}}{{{v_0}}}\ln \frac{{s + {p_{{\rm{atm}}}}}}{{{s_{\rm{c}}} + {P_{{\rm{atm}}}}}}\left[ {1 - \exp \left( { - \frac{{{t_{{\rm{sa}}}}}}{{{\psi _{\rm{s}}}}}} \right)} \right]. $

(25)

式中：ε_v0^svp是s_c时SI屈服面的黏塑性体应变，λ_s^vp是吸力黏塑性指数，ψ_s是吸力蠕变系数，t_sa为吸力恒定条件下的蠕变时间.对式(25)关于时间求导得

$ {\psi _{\rm{s}}}\dot \varepsilon _{\rm{v}}^{{\rm{svp}}} = \frac{{\lambda _{\rm{s}}^{{\rm{vp}}}}}{{{v_0}}}\ln \frac{{s + {p_{{\rm{atm}}}}}}{{{s_{\rm{c}}} + {p_{{\rm{atm}}}}}}\exp \left( { - \frac{{{t_{{\rm{sa}}}}}}{{{\psi _{\rm{s}}}}}} \right). $

(26)

把式(25)中的t_sa代入到式(26)得

$ {\psi _{\rm{s}}}\dot \varepsilon _{\rm{v}}^{{\rm{svp}}} = \varepsilon _{v0}^{{\rm{svp}}} - \varepsilon _{\rm{v}}^{{\rm{svp}}} + \frac{{\lambda _{\rm{s}}^{{\rm{vp}}}}}{{{v_0}}}\ln \frac{{s + {p_{{\rm{atm}}}}}}{{{s_{\rm{c}}} + {P_{{\rm{atm}}}}}}. $

(27)

与LC屈服面的等效时间黏塑性体应变速率性质一样，利用等效时间获得的SI屈服面的黏塑性体应变速率方程式(27)也适用于任何一种应力路径.

Alonso等^[9]认为，SI屈服面服从非相关联流动法则，它的黏塑性势函数与剪应力无关，故根据Perzyna理论可知SI屈服面不会产生黏塑性剪应变，即ζ_ij^svp=0，相应的黏塑性广义剪应变为ε_s^svp=0.

1.6 非饱和土弹黏塑性方程

根据上面分析，非饱和土弹黏塑性方程包括：(Ⅰ)体应变流变方程，它由式(1)、(3)、(10)和(27)组成；(Ⅱ)偏应变流变方程，它由式(2)、(4)、(13)和ζ_ij^svp=0组成.偏应变流变方程可合写为

$ \begin{array}{l} {{\dot \zeta }_{ij}} = \frac{{{{\dot S}_{ij}}}}{{2G}} + \dot \lambda \frac{{\partial F}}{{\partial q}}\frac{{\partial q}}{{\partial {S_{ij}}}} = \frac{{{{\dot S}_{ij}}}}{{2G}} + \varepsilon _{\rm{s}}^{{\rm{Lvp}}}\frac{{\partial q}}{{\partial {S_{ij}}}} = \frac{{{{\dot S}_{ij}}}}{{2G}} + \\ \frac{{2\psi \left( s \right)q/\left( {{\mathit{\boldsymbol{v}}_0}{t_{\rm{r}}}} \right)}}{{{M^2}\left( {2\bar p + {{\bar p}_{\rm{s}}} - {{\bar p}_{\rm{c}}}} \right)}}\exp \left[ {\frac{{\varepsilon _{\rm{r}}^{{\rm{Lvp}}}\left( s \right) - \varepsilon _{\rm{v}}^{{\rm{Lvp}}}}}{{\psi \left( s \right)/{v_0}}}} \right]{\left( {\frac{{{{\bar p}_{\rm{c}}}}}{{{{\bar p}_{\rm{r}}}}}} \right)^{\frac{1}{{\beta \left( s \right)}}}}\frac{{\partial q}}{{\partial {S_{ij}}}}. \end{array} $

(28)

在三轴试验中，净应力和净应力速率的Lode角均为-π/6，由式(28)可得

$ \begin{array}{l} {{\dot \zeta }_{\rm{s}}} = \frac{{\dot q}}{{3G}} + \frac{{2\psi \left( s \right)q/\left( {{v_0}{t_{\rm{r}}}} \right)}}{{{M^2}\left( {2\bar p + {{\bar p}_{\rm{s}}} - {{\bar p}_{\rm{c}}}} \right)}} \times \\ \;\;\;\;\;\;\;\;\exp \left[ {\frac{{\varepsilon _{\rm{r}}^{{\rm{Lvp}}}\left( s \right) - \varepsilon _{\rm{v}}^{{\rm{Lvp}}}}}{{\psi \left( s \right)/{v_0}}}} \right]{\left( {\frac{{{{\bar p}_{\rm{c}}}}}{{{{\bar p}_{\rm{r}}}}}} \right)^{\frac{1}{{\beta \left( s \right)}}}}. \end{array} $

(29)

此时第1、第2和第3主应变ε_i(i=1, 2, 3)有：

$ {{\dot \varepsilon }_1} = \left( {G{{\dot \varepsilon }_{\rm{v}}} + \dot q + 3G{{\dot \zeta }_{\rm{s}}}} \right)/3G, $

(30)

$ {{\dot \varepsilon }_2} = {{\dot \varepsilon }_3} = \left( {2G{{\dot \varepsilon }_{\rm{v}}} - \dot q - 3G{{\dot \zeta }_{\rm{s}}}} \right)/6G. $

(31)

1.7 土-水特性方程

试验表明，一些非饱和土的水力和力学之间存在较明显的相互耦合作用^[15]，表现在土-水特性上，非饱和土的饱和度除与吸力相关外，还与土骨架应变有关.它们之间的本构方程可表示为

$ {{\dot S}_{\rm{r}}} = \left( {\partial {S_{\rm{r}}}/\partial s} \right)\dot s + \left( {\partial {S_{\rm{r}}}/\partial {\varepsilon _{\rm{v}}}} \right){{\dot \varepsilon }_{\rm{v}}}. $

(32)

假定土-水特征曲线满足Genuchten方程^[16]，则主脱湿和主吸湿方程可表示为：

主脱湿方程

$ {S_{{\rm{rd}}}} = {\left[ {1 + {{\left( {{s_{\rm{d}}}/{\alpha _{\rm{d}}}} \right)}^m}} \right]^{ - n}},{\rm{d}}{s_{\rm{d}}} > 0; $

(33)

主吸湿方程

$ {S_{{\rm{rw}}}} = {\left[ {1 + {{\left( {{s_{\rm{w}}}/{\alpha _{\rm{w}}}} \right)}^m}} \right]^{ - n}},{\rm{d}}{s_{\rm{w}}} < 0. $

(34)

式中：s_d和s_w分别是主脱湿和主吸湿过程的吸力，S_rd和S_rw分别是主脱湿和主吸湿过程的饱和度，α_d、α_w、m和n分别是Genuchten方程参数.

设一般脱湿(或吸湿)过程起始时刻的饱和度和吸力为S_r0和s₀，S_r0时在主脱湿和主吸湿曲线上所对应的吸力分别为s_d0和s_w0.假定饱和度S_r相同时毛细管道的孔喉概率分布特征随孔隙半径r呈r^b自相似，b为孔喉结构分形维数.根据Young-Laplace方程可知孔隙半径r与吸力s呈反比，故对于相同饱和度S_r，s与式(33)~(34)中s_d和s_w的关系：

$ \left\{ \begin{array}{l} s_0^{ - b} - {s^{ - b}} = s_{{\rm{d}}0}^{ - b} - s_{\rm{d}}^{ - b},{\rm{d}}s > 0\left( {脱湿过程} \right);\\ {s^{ - b}} - s_0^{ - b} = s_{\rm{w}}^{ - b} - s_{{\rm{w}}0}^{ - b},{\rm{d}}s < 0\left( {吸湿过程} \right). \end{array} \right. $

(35)

由式(33)~(35)可得式(32)中等式右边的第一项为

$ \frac{{\partial {S_{\rm{r}}}}}{{\partial s}} = \left\{ {\begin{array}{*{20}{l}} {\left( {\partial {s_{\rm{d}}}/\partial s} \right)\left( {\partial {S_{{\rm{rd}}}}/\partial {s_{\rm{d}}}} \right),{\rm{d}}s > 0\left( {脱湿过程} \right);}\\ {\left( {\partial {s_{\rm{w}}}/\partial s} \right)\left( {\partial {S_{{\rm{rw}}}}/\partial {s_{\rm{w}}}} \right),{\rm{d}}s < 0\left( {吸湿过程} \right).} \end{array}} \right. $

(36)

式(32)中等式右边第二项的表达式建议取为

$ \partial {S_{\rm{r}}}/\partial {\varepsilon _{\rm{v}}} = {S_{\rm{r}}}{\left( {1 - {S_{\rm{r}}}} \right)^{{\alpha _2}}}/n. $

(37)

式中n为孔隙率，α₂为水力-力学耦合参数.

2 模型算例分析

非饱和土三维流变方程涉及的模型参数包括回弹指数κ/υ₀，压缩指数λ(s)/υ₀，次固结系数ψ(s)/υ₀、吸力回弹指数κ_s/υ₀、吸力黏塑性指数λ_s^vp(s)/υ₀、吸力蠕变系数ψ_s(s)、弹性剪切模量G、参数κ_a和M以及土水特征曲线的α_d、α_w、α₂、b、m和n.上述参数的确定方法见文献[4, 7-9, 16].

冯君^[17]对合肥地区非饱和黏土的蠕变特性进行了系统的试验研究，试验土样取自肥东县长临河铁路旁，场区表层1~2 m为第四系松软土，其下为20~25 m厚黏土，属塑性较强的下蜀组黏土，该层黏土的主要物理力学性质指标见表 1^[17].

冯君^[17]在净压力50~300 kPa和吸力0~200 kPa范围内对黏土进行排水固结蠕变和剪切试验，由此获得参数如下：λ(s)/υ₀=[6.03exp(-0.015 s)+2.88]×10^-3，ε_vr^Lvp(s)=[1.175+6.5exp(-0.018 s)]×10^-3，p_r=50 kPa，ψ(s)/υ₀=[2.12exp(-0.02 s)+3.12]×10^-4，t_r=1 d，κ/υ₀=1.05×10^-3，κ_s/υ₀=0.001 56，μ=0.25，k_a=0.3，M=1.01.当s≤s_c时可不考虑吸力黏塑性压缩影响，故取λ_s^vp=0和ψ_s=0^[9].同时，试验表明土骨架变形对土-水特性的影响较小，为节省篇幅，本算例不对土-水特征曲线做具体分析.利用式(30)~(31)可获得各吸力和蠕变持续时间下的理论预测ε_v~p曲线，与实测值对比见图 2.

从图 2可看出，随着吸力增大，非饱和黏土的弹塑性压缩线斜率变缓，压缩变形减少.同时次固结系数降低，蠕变变形也相应减小.图中预测值与实测值基本吻合，说明等效时间弹黏塑性模型能够反映该非饱和黏土的流变力学特性.图中还给出了20 a后次固结变形预测值，应变变化量约为0.31%~0.52%，当黏土层厚度达到30~40 m时将产生约9.5~12.4 cm(按0.31%计算)和15.6~20.8 cm(按0.52%计算)的次固结沉降，对于高等级公路的关键部位，如路基和桥头接触处，上述沉降量超过了公路规范规定，因此不可忽略次固结变形对黏土地基工后沉降的影响.

图 3给出了LC屈服净应力p_c随吸力s和黏塑性体应变速率$\dot \varepsilon $_v^Lvp的变化规律.

从图 3可看出，p_c随着s和$\dot \varepsilon $_v^Lvp增大而增大.当s较小时，p_c增长速率较快；随着s增大，p_c增速逐渐变缓并趋向最大值.当$\dot \varepsilon $_v^Lvp较小时，p_c随s增长幅值较小，随着$\dot \varepsilon $_v^Lvp增大，p_c随s增长幅值明显增大，呈现出随$\dot \varepsilon $_v^Lvp幂函数增大趋势，说明黏塑性应变速率对屈服净应力硬化规律具有较明显的影响，研究土的屈服特性需要考虑流变速率效应.

根据非饱和土三维流变模型可以分析一般应力作用下土体应变随时间的变化特性.设土样先在围压σ₁=σ₂=σ₃=100 kPa下固结24 h，把此时应变作为理论分析的起始零应变.然后保持σ₂=σ₃=100 kPa不变，在σ₁方向分三种工况施加荷载.q=σ₁-σ₃，图 4给出了广义剪应力q和吸力s在三种工况中的加载应力历史.在图 4三种工况中q和s都分三级加载，第一次加载幅值为Δq=60 kPa，Δs=-100 kPa，第二次加载幅值为Δq=60 kPa，Δs=-50 kPa，第三次加载幅值为Δq=30 kPa，Δs=-20 kPa，三次加载后广义剪应力从零增加到150 kPa，吸力从200 kPa减少为25 kPa.三种工况的不同点为，工况1为突然施加广义剪切力和吸力，然而维持72 h保持不变后再施加下级荷载.工况2为线性加载，加载时间12 h，然而维持60 h不变后再施加下级荷载.工况3为线性加载，但加载时间为24 h，然而维持48 h保持不变后再施加下级荷载.

采用4阶Rung-Kutta方法对体应变和广义剪应变进行了编程计算，计算结果见图 5、6.

图 5给出了体应变随时间变化图.从图 5可看出，在加载过程和维持恒载的初始阶段，三种加载历史对体应变有明显的影响，工况1的体应变＞工况2的体应变＞工况3的体应变.但随着维持恒载的时间增长，这三种工况的体应变逐渐接近.这一性质与饱和黏土的流变特性相类似.图 6给出了广义剪应变随时间变化曲线.从图 6可看出，当广义剪应力较小时，这三种加载方式引起的广义剪应变差别不大，但随着广义剪应力逐渐增大，三种加载方式引起的广义剪应变具有明显差别.它们不但在加载过程存在明显不同的增加速率，而且在恒载维持阶段也具有不同的增长幅值.工况1产生的广义剪应变明显大于工况2和工况3的广义剪应变.造成这一情况的主要原因是广义剪应变的增加速率不光与体积流变速率有关，而且还与剪切应力水平密切相关.从式(12)可看出，达到较高应力水平的时间越快，同样体积流变量产生的剪切流变量越大.工况1达到较高应力水平的时间最快，工况2次之，工况3最迟，故对于流变产生的广义剪应变值，工况1＞工况2＞工况3.这说明在流变模型中，不同加载历史会对剪应变产生较明显的影响.

3 结论

1) 根据非饱和土黏塑性体应变-净应力-等效时间关系式，利用等效时间法获得LC屈服面的黏塑性体应变速率方程，建立了LC屈服面的黏塑性广义剪应变速率方程.根据SI屈服面的常吸力压缩蠕变公式，利用等效时间法建立了SI屈服面的黏塑性体应变速率方程，SI屈服面不产生黏塑性剪应变.

2) 当非饱和土退化为饱和土时，本文非饱和土等效时间流变模型退化为饱和土的Yin-Graham流变模型.当不考虑流变特性时，本文非饱和土流变模型能够退化为巴塞罗那非饱和土弹塑性模型.

3) 数值分析了非饱和土一维压缩蠕变特性，研究了三轴不同加载方式作用下非饱和土体应变和广义剪应变随时间变化规律.结果表明，加载方式只对加载期间和恒荷初始阶段的体应变有影响，对长期体应变值影响不大；无论在加载还是在荷载恒定阶段，加载方式均对剪应变产生较明显的影响.