多年冻土区路基沉降变形是冻土工程研究中的核心，也是铁路正常运营之安全隐患，更是制约寒区工程建设和资源开发的关键因素.冻土因其冰的存在，即使在很小的荷载下也会出现强度和变形随时间而变化的蠕变现象^[1-2].青藏铁路通车近十年的现场实测资料表明，部分多年冻土路段路基不断补碴、抬道使道床厚度甚至达到1.7 m以上^[3]，路基稳定性、行车安全问题日渐凸显，冻土蠕变机理与沉降控制措施研究迫切需要深入开展.

冻土的蠕变特性受土的类型、密度、含水率、温度、应力水平等多种因素影响^[4-5]，在变形过程中还同时存在强化和弱化作用^[6].蠕变曲线分为衰减蠕变和非衰减蠕变两种类型，建立冻土蠕变模型常用3种途径^[7]：1)基于试验结果通过数学公式推导、拟合建立的经验模型.2)应力-应变-时间模型.Bray^[8]基于松弛试验得到了冻土的蠕变参数，并研究了温度和应力对冻土和多晶冰蠕变特性的影响; Arenson等^[9]通过冻土蠕变试验，得到温度接近于0 ℃富冰冻土的蠕变模型; 周志伟等^[10]进行不同温度和围压条件下的三轴压缩和蠕变试验，基于超塑性理论和过应力原理建立了冻结黄土的率相关本构模型; 朱元林等^[11]通过冻土动三轴蠕变试验，回归分析得出了冻土动蠕变模型.3)利用基本的虎克弹性体、牛顿黏性体、圣-维南塑性体等力学模型元件，构建的流变模型.经典的流变模型有Kelvin模型、Burgers模型、Bingham模型以及Nishihara模型等，其模型由一系列力学模型元件组成，概念直观被广泛应用于描述岩土材料的蠕变特性.但这些模型采用的都是线性元件，无法很好地描述蠕变的全过程，尤其是无法描述加速蠕变^[12].故此，学者们提出采用非线性流变元件代替常规的线性流变元件，或者引入损伤断裂力学等新的理论建立冻土蠕变模型.康永刚等^[13]用幂律流体元件替代牛顿黏壶，提出一种经验函数描述加速蠕变阶段; 孙凯等^[14]将黏弹性元件中黏滞系数修正为关于时间的函数来描述非线性蠕变过程，并将一维蠕变模型推广到三维中; 赵晓东等^[15]通过GFC(不经历K₀固结，不均匀温度冻结)的方法，对不同温度梯度和蠕变应力下的冻结饱和黏土进行了一系列单轴蠕变试验，并建立了考虑热梯度的修正Kelvin模型; 赵延林等^[16]在Burgers模型基础上，串联一个能描述蠕变过程中强度参数非线性衰减特性的M-C塑性元件，并利用FLAC^3D实现了新建蠕变损伤模型的二次开发.李栋伟等^[17]引入服从D-P屈服准则的黏塑损伤变量，结合西原模型推导出冻土黏弹塑性损伤耦合本构方程.各学者针对建立的蠕变模型参数的影响因素也进行了深入探讨，获得大量有益认识.

冻土的蠕变是典型的非线性流变，应力和时间都能引起非线性，且对非线性的影响是耦合的.要想合理描述冻土非线性蠕变特征，修正的非定常黏滞系数应是流变时间和应力水平的函数^[18]，然而大部分改进的非线性黏弹性元件很少有考虑应力水平对黏滞系数的影响.鉴于此，以Nishihara模型为基础，将黏弹性元件中的黏滞系数修正为应力水平和流变时间的函数，并在黏塑性元件中引入损伤变量D，进而建立简单而适用的冻结砂土的蠕变本构关系.不同温度和干密度条件下的冻结砂土三轴蠕变试验结果对比表明，该模型能较好地描述不同温度、应力水平条件下冻结砂土的衰减、稳态和加速蠕变特性.本模型可为冻结砂土蠕变预测提供一种新的选择.

1 冻土三轴蠕变试验 1.1 试验条件

试验用土为青藏铁路沿线砂土，初始含水率15.6%，颗粒级配曲线见图 1.采用标准方法批量制备试件^[19]，制成61.8 mm×125 mm的圆柱形试件.试验在MTS Landmark370.10材料试验机上进行(见图 2).仪器控温精度为±0.1 ℃，最大轴向加载100 kN.在0.1 MPa围压下开展5个温度(-0.5, -0.9, -1.2, -1.5和-2.0 ℃)、两种干密度(1.75和1.92 g/cm³)的三轴蠕变试验.当土体发生非衰减蠕变时，轴向应变达24%试验终止; 当土体发生衰减蠕变时，轴向应变速率≤0.000 5/h试验终止.

1.2 试验结果和分析

冻土在不同温度、干密度和剪切应力条件下的蠕变曲线见图 3，ε₁为轴向应变，q为广义剪应力，t为蠕变时间.可以看出，当剪应力水平较低时，硬化作用占优势，表现为衰减蠕变; 随着剪应力增大，弱化作用占优势，表现为非衰减蠕变.当冻土干密度为1.75 g/cm³时，非衰减蠕变过程主要表现为非稳定和稳定蠕变两个阶段，加速阶段并不明显; 而当冻土干密度为1.92 g/cm³时，加速蠕变阶段明显.不同温度和相同应力条件下冻土的三轴蠕变曲线见图 4，θ为试验温度.可以看出，温度越高，未冻水含量增多，冻土内部矿物颗粒和冰胶结作用削弱，蠕变量较大.

2 冻土一维蠕变模型

Nishihara模型由一个广义Kelvin体和一个黏塑性体串联而成(见图 5)，可用来描述岩土材料的蠕变特性^{[14, 18]}，式(1)为蠕变方程:

$ \varepsilon = \left\{ {\begin{array}{*{20}{l}} {\frac{\sigma }{{{E_0}}} + \frac{\sigma }{{{E_1}}}\left( {1 - {{\rm{e}}^{ - \frac{{{E_1}}}{{{\eta _1}}}t}}} \right),\sigma < {\sigma _{\rm{s}}};}\\ {\frac{\sigma }{{{E_0}}} + \frac{\sigma }{{{E_1}}}\left( {1 - {{\rm{e}}^{ - \frac{{{E_1}}}{{{\eta _1}}}t}}} \right) + \frac{{\sigma - {\sigma _{\rm{s}}}}}{{{\eta _2}}}t,\sigma \ge {\sigma _{\rm{s}}}.} \end{array}} \right. $

(1)

式中：E₀为虎克弹性体的弹性模量，E₁和η₁分别为Kelvin体中的弹性模量和黏滞系数，η₂为黏塑性体中的黏滞系数，σ为蠕变应力，σ_s为屈服应力.

由式(1)可知，该模型能够反映材料的弹性、黏弹和黏塑等蠕变特性，适用范围较广^[18].但由于其组成元件均为线性的，无法描述加速蠕变过程.故采用非线性元件替换线性元件，对Nishihara模型进行改进，见图 6.

2.1 黏弹性部分

黏弹性部分由一个弹簧和一个非线性黏性元件并联组成.冻土蠕变具有明显的非线性特征，其非线性程度与荷载大小及荷载持续时间有关.为反映冻土蠕变的非线性特性，参考文献[14, 20-21]将黏滞系数修正为应力水平和流变时间的函数，即

$ {\eta _1}\left( {t,\sigma } \right) = {\eta _{10}}\left( {\frac{{2{\sigma _0} - \sigma }}{\sigma }} \right)\frac{{{{\left( {t/{t_0}} \right)}^{1/3}}}}{{\left[ {{{\left( {t/{t_0}} \right)}^{1/3}} + 1} \right]}}. $

(2)

式中：σ₀为相同温度和围压条件下的常规三轴抗压强度，σ为蠕变应力大小，η₁₀为黏弹性体的初始黏滞系数，t₀为流变时间参考值(h), 本试验t₀=1.

式(2)对时间求导得

$ {{\dot \eta }_1}\left( {t,\sigma } \right) = {\eta _{10}}\left( {\frac{{2{\sigma _0} - \sigma }}{\sigma }} \right)\frac{{\frac{1}{3}{{\left( {t/{t_0}} \right)}^{ - 2/3}}}}{{{{\left[ {{{\left( {t/{t_0}} \right)}^{1/3}} + 1} \right]}^2}}} \ge 0. $

(3)

由式(2)，(3)可知，当t=0时，η₁(t, σ)=0;当t→∞，σ不变时，$\eta_{1}(t, \sigma) \rightarrow \eta_{10}\left(2 \sigma_{0}-\sigma\right) / \sigma $，说明当t从0→∞时，η₁(t, σ)单调递增至η₁₀(2σ₀-σ)/σ，黏滞系数增大会导致土的蠕变速率逐渐减小，能够描述土体在低应力状态下的衰减蠕变过程^[22].而当应力水平增大时，黏滞系数减小，能够反映冻土蠕变速率随应力增大而增大.改进的Kelvin体本构方程为

$ \sigma = {E_1}\varepsilon \left( t \right) + {\eta _1}\left( {t,\sigma } \right)\frac{{{\rm{d}}\varepsilon \left( t \right)}}{{{\rm{d}}t}}. $

(4)

在常应力条件下对式(4)分离变量求定积分(t=0，ε(t)=0)，得到非定常Kelvin体的蠕变方程为

$ \varepsilon \left( t \right) = \frac{{{\sigma _0}}}{{{E_1}}}\left[ {1 - \exp \left( { - \alpha t - \frac{3}{2}\alpha {{\left( {t/{t_0}} \right)}^{2/3}}} \right)} \right]. $

(5)

其中

$ \alpha = \frac{{{E_1}}}{{{\eta _{10}}}}\left( {\frac{\sigma }{{2{\sigma _0} - \sigma }}} \right). $

(6)

修正后的具有非线性元件的Kelvin体可用于描述冻结砂土的衰减蠕变特性.

2.2 黏塑性部分

大量研究表明^{[4, 23]}，冻土中存在损伤效应，当剪应力大于某一个临界值时，损伤效应就会出现，当损伤效应占优势时，会发生非衰减蠕变^[24].通过引入损伤变量D定量描述非衰减蠕变过程中的稳定和加速蠕变阶段，基于损伤机理定义损伤变量^{[4, 19]}

$ D\left( {\sigma ,t} \right) = 1 - \frac{{E\left( {\sigma ,t} \right)}}{{{E_0}}}. $

(7)

$ E(\sigma ,t) = \left\{ {\begin{array}{*{20}{l}} {{E_0},\sigma - {\sigma _\infty } < 0;}\\ {{E_0}{{\rm{e}}^{ - c{t^R}}},\sigma - {\sigma _\infty } \ge 0.} \end{array}} \right. $

(8)

$ D(\sigma ,t) = \left\{ {\begin{array}{*{20}{l}} {0,\sigma - {\sigma _\infty } < 0;}\\ {1 - {{\rm{e}}^{ - c{t^R}}},\sigma - {\sigma _\infty } \ge 0.} \end{array}} \right. $

(9)

式中：D(σ, t)为损伤变量，弹性模量E(σ, t)为随时间变化的函数，E₀为初始弹性模量，c，R为反映损伤程度的材料参数，σ_∞为长期强度，是指稳定蠕变阶段的轴向蠕变速率为0时对应的最大剪应力.可以基于蠕变稳定阶段的速率法来确定长期强度^[19]，即

$ q = {\sigma _\infty } + k\ln \left( {1 + \frac{\nu }{{{\nu _0}}}} \right). $

(10)

式中：q为剪切应力，k为试验参数(MPa)，v为稳定阶段蠕变速率，v₀为蠕变速率参考值，本试验v₀=1.0.

根据有效应力定义^[25]：

$ {\sigma ^\prime } = \frac{\sigma }{{1 - D}}. $

(11)

$ {\sigma ^\prime } = \left\{ {\begin{array}{*{20}{l}} {\sigma ,\sigma - {\sigma _\infty } < 0;}\\ {\sigma {{\rm{e}}^{c{t^R}}},\sigma - {\sigma _\infty } \ge 0.} \end{array}} \right. $

(12)

式中:σ′为有效应力，σ为蠕变应力.

采用由非定常黏性体与塑性体并联的黏塑性元件来描述冻土加速蠕变阶段的非线性特征.当σ＜σ_∞，塑性体不发生变形，此时，黏性体不起作用; 当σ≥σ_∞时，塑性体屈服.因此，非定常黏性元件的黏滞系数可以表示

$ {\eta _2}(t,\sigma ) = {\eta _{20}}\left[ {1 + \frac{{2{\sigma _\infty }}}{\sigma }\left( {1 - {{\rm{e}}^{ - t/{t_0}}}} \right)} \right]. $

(13)

η₂₀为黏塑性元件的初始黏滞系数.

根据串联时合应变为各部分应变之和得到蠕变函数：

$ \begin{array}{*{20}{c}} {\varepsilon = \frac{\sigma }{{{E_0}}} + \frac{\sigma }{{{E_1}}}\left[ {1 - \exp \left( { - \alpha t - \frac{3}{2}\alpha {{\left( {t/{t_0}} \right)}^{2/3}}} \right)} \right],}\\ {\sigma < {\sigma _\infty },} \end{array} $

(14)

$ \begin{array}{l} \varepsilon = \frac{\sigma }{{{E_0}}} + \frac{\sigma }{{{E_1}}}\left[ {1 - \exp \left( { - \alpha t - \frac{3}{2}\alpha {{\left( {t/{t_0}} \right)}^{2/3}}} \right)} \right] + \\ \;\;\;\;\;\frac{{{\sigma ^\prime } - {\sigma _\infty }}}{{{\eta _2}\left( {t,\sigma } \right)}}t,\sigma \ge {\sigma _\infty }. \end{array} $

(15)

3 冻土三维蠕变模型

三维应力条件下，总应变可以表示为

$ {\varepsilon _{ij}} = \varepsilon _{ij}^{\rm{e}} + \varepsilon _{ij}^{{\rm{ve}}} + \varepsilon _{ij}^{{\rm{vp}}}. $

(16)

根据广义虎克定律，弹性体三维本构关系为

$ {\sigma _{kk}} = 3K{\varepsilon _{kk}},{S_{ij}} = 2G{e_{ij}}, $

(17)

式中：S_ij、e_ij分别为应力偏张量和应变偏张量，σ_kk、ε_kk分别为应力张量和应变张量第一不变量，G和K分别为剪切模量和体积模量.

结合式(17)，弹性体、黏弹性体和黏塑性体的应变分别为

$ \varepsilon _{ij}^{\rm{e}} = \frac{1}{{3K}}{\sigma _{\rm{m}}}{\delta _{ij}} + \frac{{{S_{ij}}}}{{2{G_0}}}, $

(18)

$ \varepsilon _{ij}^{{\rm{ve}}} = \frac{{{S_{ij}}}}{{2{G_1}}}\left[ {1 - \exp \left( { - \beta t - \frac{3}{2}\beta {{\left( {t/{t_0}} \right)}^{2/3}}} \right)} \right], $

(19)

$ \varepsilon _{ij}^{{\rm{vp}}} = \frac{1}{{{\eta _2}\left( {t,q} \right)}}\left[ {\varphi \left( {\frac{F}{{{F_0}}}} \right)\frac{{\partial Q}}{{\partial \sigma _{ij}^\prime }}} \right]t. $

(20)

将式(18)、(19)和(20)代入式(16)，三维蠕变本构模型可统一表示为

$ \begin{array}{l} {\varepsilon _{ij}} = \frac{1}{{3K}}{\sigma _{\rm{m}}}{\delta _{ij}} + \frac{{{S_{ij}}}}{{2{G_0}}} + \\ \frac{{{S_{ij}}}}{{2{G_1}}}\left[ {1 - \exp \left( { - \beta t - \frac{3}{2}\beta {{\left( {t/{t_0}} \right)}^{2/3}}} \right)} \right],q < {\sigma _\infty }; \end{array} $

(21)

$ \begin{array}{l} {\varepsilon _{ij}} = \frac{1}{{3K}}{\sigma _{\rm{m}}}{\delta _{ij}} + \frac{{{S_{ij}}}}{{2{G_0}}} + \\ \;\;\;\;\;\;\;\frac{{{S_{ij}}}}{{2{G_1}}}\left[ {1 - \exp \left( { - \beta t - \frac{3}{2}\beta {{\left( {t/{t_0}} \right)}^{2/3}}} \right)} \right] + \\ \;\;\;\;\;\;\;\frac{1}{{{\eta _2}\left( {t,q} \right)}}\left[ {\varphi \left( {\frac{F}{{{F_0}}}} \right)\frac{{\partial Q}}{{\partial \sigma _{ij}^\prime }}} \right]t,q \ge {\sigma _\infty }. \end{array} $

(22)

式中：σ_mδ_ij为球应力张量，G₀为弹性元件的剪切模量，G₀=q/3ε₀，ε₀为初始应变值，q为剪应力，G₁为黏弹性元件的剪切模量，F为屈服函数，F₀为屈服函数初始参考值，本试验F₀=1，Q为塑性势函数，假设试验材料满足相关联流动法则，F=Q，φ(·)为幂函数形式，其指数符号为1，σ′_ij为柯西有效应力张量.

本次试验假定应力偏量在冻土蠕变中起主要作用，因此，屈服函数表达式为^[26-27]

$ F = \sqrt {{J_2}} - \frac{{{\sigma _{\rm{s}}}}}{{\sqrt 3 }}. $

(23)

其中:J₂为第二偏应力不变量.在常规三轴试验中，σ₂=σ₃，因此，

$ {\sigma _{\rm{m}}} = \frac{1}{3}\left( {{\sigma _1} + 2{\sigma _3}} \right), $

(24)

$ {S_{11}} = {\sigma _1} - {\sigma _{\rm{m}}} = \frac{2}{3}\left( {{\sigma _1} - {\sigma _3}} \right), $

(25)

$ \sqrt {{J_2}} = \frac{1}{{\sqrt 3 }}\left( {{\sigma _1} - {\sigma _3}} \right). $

(26)

最终三维蠕变模型为

$ \begin{array}{*{20}{c}} {{\varepsilon _1} = \frac{{{\sigma _1} + 2{\sigma _3}}}{{9K}} + \frac{{{\sigma _1} - {\sigma _3}}}{{3{G_0}}} + }\\ {\frac{{{\sigma _1} - {\sigma _3}}}{{3{G_1}}}\left[ {1 - \exp \left( { - \beta t - \frac{3}{2}\beta {{\left( {t/{t_0}} \right)}^{2/3}}} \right)} \right],}\\ {q < {\sigma _\infty };} \end{array} $

(27)

$ \begin{array}{*{20}{c}} {{\varepsilon _1} = \frac{{{\sigma _1} + 2{\sigma _3}}}{{9K}} + \frac{{{\sigma _1} - {\sigma _3}}}{{3{G_0}}} + }\\ {\frac{{{\sigma _1} - {\sigma _3}}}{{3{G_1}}}\left[ {1 - \exp \left( { - \beta t - \frac{3}{2}\beta {{\left( {t/{t_0}} \right)}^{2/3}}} \right)} \right] + }\\ {\frac{{{{\left( {{\sigma _1} - {\sigma _3}} \right)}^\prime } - {\sigma _\infty }}}{{3{\eta _2}\left( {t,q} \right)}}t,q \ge {\sigma _\infty }.} \end{array} $

(28)

其中

$ \beta = \frac{{{G_1}q}}{{{\eta _{10}}\left( {3{\sigma _0} - q} \right)}}, $

$ {\eta _2}\left( {t,q} \right) = {\eta _{20}}\left[ {1 + \frac{{3{\sigma _\infty }}}{q}\left( {1 - {{\rm{e}}^{ - t/{t_0}}}} \right)} \right], $

$ {\left( {{\sigma _1} - {\sigma _3}} \right)^\prime } = \frac{{\left( {{\sigma _1} - {\sigma _3}} \right)}}{{1 - D}} = \left( {{\sigma _1} - {\sigma _3}} \right){{\rm{e}}^{c{t^R}}}. $

4 模型参数确定及模型验证

改进的模型有K、G₀、G₁、c、R、η₁₀、η₂₀共7个参数，能通过蠕变试验曲线直接测得或者回归分析获得.基于试验数据，模型参数见表 1.ρ_d为土样干密度，拟合曲线与试验曲线对比见图 7(仅列出部分温度条件下的参数值和拟合结果).为了验证模型的适用性，将温度为-2.0 ℃、干密度为1.92 g/cm³的冻结砂土试验结果，与Nishihara模型、孙凯等^[14]建立的模型和改进模型的模拟结果绘在同一图幅内，结果见图 8.可以看出，3种模型均能够较好地描述低应力条件下冻结砂土的衰减蠕变特性，但Nishihara模型不能很好地反映非衰减蠕变过程; 孙凯等^[14]建立的模型，考虑了时间对模型元件的影响，能够描述加速蠕变过程，与孙凯等^[14]建立的模型相比，改进模型考虑了时间和应力对模型元件的耦合影响，不仅能描述衰减和稳态蠕变过程，还能反映加速蠕变特征，拟合精度也要高于Nishihara模型和孙凯等^[14]建立的模型.由图 7, 8可知，改进模型拟合精度较高，数据吻合良好，表明该模型能较好地描述不同温度和剪应力条件下冻结砂土的蠕变过程.

温度和剪应力对模型参数的影响见图 9，可以看出，改进模型中黏弹性元件的黏滞系数η₁₀、剪切模量G₁和黏塑性元件的黏滞系数η₂₀均随着温度降低而增大，随应力水平增大而减小，损伤变量参数c随温度降低而减小，随应力水平增大而增大.

为厘清温度、剪应力与模型参数的关系以及为不同温度、剪应力下模型参数的确定提供参考，回归分析得到模型参数与剪应力q、温度θ的函数关系：

$ \begin{array}{*{20}{c}} {\ln {G_1} = {A_1} + {A_2}q,}\\ {\ln {\eta _{10}} = {B_1} + {B_2}q.} \end{array} $

(29)

其中

$ {A_1} = 0.84 - 1.035\theta ,R = 0.9878, $

$ {A_2} = - 0.715 - 0.306{\theta ^{ - 1}},R = 0.9378, $

$ {B_1} = - 0.242 - 1.369\theta ,R = 0.9786, $

$ {B_2} = {\left( {79.775 + 67.83\theta } \right)^{ - 1}},R = 0.8040. $

$ c = {C_1} + {C_2}q + {C_3}{q^3}, $

$ R = {D_1} + {D_2}q + {D_3}{q^3}. $

(30)

其中

$ {C_1} = - 1.805 - 2.944\theta ,R = 0.9486, $

$ {C_2} = 2.345 + 3.207\theta ,R = 0.9052, $

$ {C_3} = 3.039 + 9.563\theta + 9.338{\theta ^2} + 2.598{\theta ^3},R = 0.9945, $

$ {D_1} = \frac{{ - 0.077\theta }}{{1 + 0.694\theta }},R = 0.9626, $

$ {D_2} = \frac{{0.193 + 0.211\theta }}{{1 + 0.698\theta }},R = 0.9384, $

$ {D_3} = {\left( {14.828 + 16.755{\theta ^3}} \right)^{ - 1}},R = 0.8964. $

$ \ln {\eta _{20}} = {N_1} + {N_2}q, $

$ \ln K = {W_1} + {W_2}q. $

(E31)

其中

$ {N_1} = \frac{{2.663 - 2.792{\theta ^2}}}{{1 - 0.81{\theta ^2}}},R = 0.9224, $

$ {N_2} = 0.023{\left( {1.195 + \theta } \right)^{ - 1}},R = 0.9761, $

$ {W_1} = 4.247 + 0.613{\theta ^{ - 1}},R = 0.9623, $

$ {W_2} = \frac{{ - 1.186 - 0.935\theta }}{{1 + 0.802\theta }},R = 0.9049. $

5 结论

1) 试验表明温度对冻结砂土的蠕变特性影响显著，尤其当温度高于-1.5 ℃时，即使在很小的荷载作用下也能产生明显的蠕变变形.

2) 以Nishihara模型为基础，将黏弹性元件中定常的黏滞系数修正为时间和剪切应力的函数，在黏塑性元件中引入损伤变量D，来描述非衰减蠕变的稳定和加速蠕变阶段，建立了冻结砂土弹黏塑性蠕变模型.该模型能较好地描述不同温度、应力水平下冻结砂土的衰减、稳态和加速蠕变特性.拟合结果与试验数据吻合较好.

3) 改进模型中黏弹性元件的黏滞系数、剪切模量和黏塑性元件的黏滞系数均随着应力水平增大而减小，随温度降低而增大.损伤变量随着应力水平增大而增大，随温度降低而减小.本研究可为冻土蠕变理论研究积累资料.