对蠕变模型的研究由来已久，近5年来，赵延林等^[1]提出岩石黏弹塑性应变分离的蠕变试验方法和数据处理技巧.在FLAC3D中实现BNMC蠕变损伤本构模型的二次开发, 提出该模型参数的辩识方法.王军保等^[2]通过将表征应力水平对岩石非线性蠕变特性影响的函数引入到常规分数阶黏滞体的本构关系中, 提出一种改进的分数阶黏滞体，能够描述岩石加速蠕变变形.杨振伟等^[3]基于三维颗粒流程序(PFC3D), 采用Burgers流变模型，分析了伯格斯模型中的参数对瞬时强度特性和流变特性的影响.杨红伟等^[4]以细粒砂岩为研究对象, 对三轴压缩条件下岩石孔隙水压力分级加载蠕变试验进行了蠕变特性及模型研究.王军保等^[5]基于损伤力学理论, 通过引入Kachanov蠕变损伤演化规律, 构建了可反映岩石加载瞬时变形和加速蠕变的考虑时效损伤的弹性体元件, 建立了一个新的2元件非线性蠕变损伤模型及其蠕变方程.康永刚等^[6]提出用非牛顿黏壶代替Kelvin模型中的牛顿黏壶, 串联用一种非定常、非牛顿黏壶和塑性体并联为黏塑性体, 黏塑性体关闭时, 所得蠕变函数可以描述衰减蠕变和稳态蠕变.黏塑性体开启时, 蠕变函数可以描述加速蠕变变形.蒋海飞等^[7]采用幂函数、对数函数混合方程对加速蠕变段进行拟合, 再通过类比的方法提出一个新的非线性黏性元件, 并将其与塑性体并联, 得到一个可以反映岩石加速蠕变特性的非线性黏塑性模型, 将该模型与Burgers模型串联, 构建一个新的六元件非线性黏弹塑性蠕变模型.此外，还提出过一个能反映不同孔隙水压作用下加速蠕变特征的二元件黏塑性蠕变模型^[8], 将该模型与修正的广义Kelvin模型串联组成一个新的改进西原模型.范翔宇等^[9]利用饱依丁-汤姆逊体与宾汉姆体, 建立储气层煤岩蠕变模型.刘浪等^[10]分析深部岩石在饱水情况下的流变规律, 并基于统一流变力学模型, 对深部饱水岩石的流变模型进行辨识, 确定流变模型的相关参数.乔丽苹等^[11]提出采用不可恢复应变作为内变量来描述复杂条件下岩石蠕变性质, 并据此提出了岩石的内变量蠕变模型.宋勇军等^[12]将含分数阶导数的软体元件与虎克体串联, 引入能反映应力水平和时间影响的损伤变量, 提出一种四元件非线性蠕变损伤模型.并借鉴经典元件模型的建模思路^[13], 将岩石的初始屈服强度作为蠕变硬化的应力阈值, 岩石的长期强度作为损伤软化的应力阈值, 引入能反映岩石硬化效应的硬化函数和损伤效应的损伤变量, 建立能够全面反映蠕变机制的岩石非线性蠕变模型.高赛红等^[14]基于蠕变过程中的硬化-损伤机制，对高应力条件下衰减蠕变阶段的黏滞系数(硬化黏滞系数)进行修正，结合定常蠕变阶段的损伤蠕变机制, 提出高应力条件下岩石损伤效应的累积-扩散机理, 并据此引入以累积损伤蠕变量为判定准则的加速蠕变触发模型(非线性蠕变体)描述岩石的不稳定蠕变特征.在此基础上, 引入瞬时塑性元件, 与虎克体和黏弹塑性体串联, 建立了一个能够完整描述岩石蠕变全过程的非线性黏弹塑性蠕变模型.

与本文有直接联系的论文成果有：杨圣奇等^[15]基于在岩石剪切流变仪上得到的泥板岩剪切流变试验结果, 采用库仑强度准则获得了岩石长期强度参数(长期黏聚力和长期内摩擦系数), 与瞬时强度参数进行对比, 结果表明, 岩石黏聚力对时间的敏感程度高于内摩擦系数对时间的敏感程度，提出了一个考虑黏聚力和内摩擦系数随时间变化的塑性元件, 将其与线性黏性元件并联起来, 得到一个新的考虑黏聚力和内摩擦系数的非线性黏塑性体(CF-NVPB), 同时，将CF-NVPB模型与三元件线性黏弹性流变模型串联起来, 建立了一个新的考虑黏聚力和内摩擦系数的岩石非线性黏弹塑性流变模型, 并推导了该非线性流变模型的蠕变方程与松弛方程.佘成学等^[16]在西原模型的基础上, 研究建立了一种新的岩石非线性蠕变模型.首先, 在分析岩石宏、细观蠕变过程的基础上, 提出了蠕变参数变化率概念; 以此为基础, 统一建立了等速蠕变和加速蠕变阶段的岩石蠕变参数变化公式.然后, 利用该非线性蠕变参数以及根据岩石的屈服软化过程, 以西原模型为基础, 建立了岩石的非线性弹-黏弹-黏塑性蠕变模型.最后, 将单轴压缩蠕变试验结果与数值计算结果对比, 验证了模型对于描述岩石全过程蠕变变形破坏的正确性.

通常认为岩石强度伴随着结构破损，即塑性剪切应变的增加而降低，与塑性剪应变或时间有关，相对于岩石而言，特殊土的力学性质更为复杂、强度敏感性更强，强度伴随时间的变化规律更复杂，即使没有塑性剪切应变发生，原状土受到扰动，强度也会发生先减后增的变化，探讨特殊土的关于变参数蠕变模型的具体形式及应用方式的论文成果有：

孔令伟等^[17]基于湛江黏土的蠕变变形在低偏应力下较小, 其蠕变等时曲线没有明显的屈服特征而近似为直线, 表现为线性蠕变性状; 当偏应力达到某一临界值后, 在短时间内发展到加速蠕变阶段并发生破坏.从湛江黏土蠕变各阶段的变形机制出发, 通过引入服从Mohr-Coulomb塑性流动规律的变参数塑性元件, 将该塑性元件与广义Kelvin模型串联, 提出描述蠕变全过程的黏弹塑性模型.同时将该模型在FLAC3D中进行二次开发, 与FLAC3D中自带的Cvisc模型进行对比验证, 结果表明, 基于变参数塑性元件的黏弹塑性模型能较好地描述湛江黏土各个蠕变阶段的变形性状.

为了尽可能准确描述湛江黏土蠕变全过程，并试图从土工材料流变学角度能够给出合理解释^[18]，需要引入塑性元件，塑性元件服从的强度准则有很多种，这里还是采用FLAC3D描述土体常用的Mohr-Coulomb塑性元件.传统塑性元件又称为摩阻元件，其是一个线性元件.当施加的应力大于该元件的摩阻力时，元件将随时间产生持续增加的变形，一旦外力卸去，已发生的变形也不会恢复，成为永久应变.然而，由莫尔-库伦强度理论，土体的塑性发生通常与土体的黏聚力c和内摩擦角φ两个强度参数有关.本文提出变强度参数蠕变模型应用于湛江黏土，探讨强度参数c和φ的变化对土体加速蠕变阶段的影响，试图得到一些基本规律为实际工程建设提供一些建议和参考.

本文的研究对于完善非线性蠕变模型理论，促进与推动关于能描述加速蠕变阶段的蠕变全过程理论模型的构建，具有重要的意义和价值.

1.   Cvisc模型

Cvisc模型是由传统经典Burger模型和一个塑性元件串联而成，见图 1，使得该本构模型能够描述材料的弹性、黏性、黏弹及塑性变形.当施加的应力σ＜σ_s时，只会产生前3部分变形，不包含塑性变形；若σ＞σ_s时，则包含全部形式的变形，此时σ_s称为屈服应力.在FLAC3D中，该塑性元件满足莫尔-库伦强度准则.

图  1  Cvisc模型

Figure  1.  Cvisc model

下载: 全尺寸图片 幻灯片

Cvisc模型一维本构关系为

式中ε_v(tσ)表示塑性体变形量是关于时间和应力的函数，其他参数同流变学文献一致^[19]，这里不再赘述.

2.   变强度参数

2.1   变强度参数的引入

在工程建设中，岩土的蠕变问题是个需要重视的棘手问题，尤其在地基处理、边坡施工和围岩开挖与支护等工程中，岩土长时间在恒定外荷载作用下，产生的持续变形会破坏其上建筑(或构筑物)的稳定性和安全性；尤其是加速蠕变的出现，给正在建设或已经完工的工程带来巨大的破坏和损失.而由摩尔-库伦强度准则，在外荷载不变的情况下，强度参数的大小决定土体内部是否存在达到强度屈服准则的区域，即塑性区，而一旦产生塑形区便意味着该部分区域的变形量将可能不断增加.对于具有蠕变特性的土体而言，即会出现加速蠕变现象.所以强度参数对土体加速蠕变阶段出现与否存在一定的影响.

另外大量的土工试验及研究表明^[20-21]：土体在受荷产生变形直至破坏的过程中，其强度参数并非一成不变，而是随着位移的发生，强度有增有减，另外土体的破坏过程区别于岩石，其具有一个明显的渐变过程，属于塑性破坏，而土体受荷产生的位移量是关于时间的函数.这种受荷呈现出强度随时间变化的特征，表现出土体作为一种三相材料体，即具有颗粒滑移的一般特性又具有一定的自身结构和黏结能力的基本特征.

本文根据对湛江黏土三轴试验获得的瞬时强度指标和长期强度指标来看，蠕变引起的强度衰减主要表现为黏聚力c值的明显降低，而摩擦角φ变化不大^[22].

基于上述2个指标在时间变化上表现出的差异，引入变强度参数模型，考虑采用FLAC3D所自带的FISH语言功能定义一个函数，实现土体强度参数c和φ随蠕变时间发生变化，试图通过该模型来描述160 kPa偏应力下湛江黏土的加速蠕变阶段.

2.2   变强度参数的建立

文中变强度参数模型是在前述Cvisc模型的基础上建立起来，通过FLAC3D内置FISH语言定义函数，实现强度参数c和φ随着蠕变时间不断发生变化，此外，实际的土体在受荷后，即使强度产生衰减导致发生破坏，其依然具有一定的残余强度，也即其强度参数c和φ不会为零，所以在建立强度参数随时间变化的过程中，要么给定2个强度指标降低的一个最低限值；要么确定蠕变时间，通过定义它们的降低速率来控制最低值.

下面是通过FISH语言定义一个函数，在函数内部定义2个强度参数值随蠕变时间不断减小，具体实现代码如下：

def cvs_cs_gb

whilestepping

p_z = zone_head

loop while p_z # null

p_coh = y_coh -j_xl * crtime

p_fri = z_frw -f_js * crtime

z_prop(p_z, 'coh') = p_coh

z_prop(p_z, 'fri') = p_fri

p_z = z_next(p_z)

endloop

end

set y_coh = 38e3;初始黏聚力为38 kPa

set j_xl = 2e4

set z_frw = 25;初始内摩擦角为25°

set f_js = 5

cvs_cs_gb

上述代码2个强度参数随蠕变时间是不断降低的，用户可以根据实际情况添加if..else..条件语句来控制它们降低的最低限值.然后将代码插入到三轴蠕变命令流文本中，通过solveagexx来设置2个不同的求解时间点，分别在软件界面通过点击Plootitems-Add-Block-property来选择查看黏聚力和内摩擦角是否发生变化，上述代码经过验证能够实现2个强度参数随蠕变时间发生改变，见图 2.

图  2  2 000时步强度参数c和φ

Figure  2.  Intensity of the parameters c and φ with 2 000 times

下载: 全尺寸图片 幻灯片

3.   c和φ对加速蠕变的影响

3.1   c、φ变化的9种情形

采用的模型为如图 2所示的圆柱体，高为1个单位，直径为0.5个单位，主要分析强度参数对湛江黏土加速蠕变阶段出现的影响，虽然影响土体加速蠕变阶段出现的因素很多，如荷载水平、湿度、颗粒级配和结构等，但在这里暂不考虑.为了通过数值模拟，准确分析2个强度参数对加速蠕变出现的影响程度，需要将强度参数变化分成以下9种情形进行简要分析：

1) 黏聚力和内摩擦角都不变.

2) 黏聚力减小，内摩擦角不变.

3) 黏聚力增大，内摩擦角不变.

4) 内摩擦角减小，黏聚力不变.

5) 内摩擦角增大，黏聚力不变.

6) 黏聚力先增后减，内摩擦角不变.

7) 黏聚力先减后增，内摩擦角不变.

8) 内摩擦角先增后减，黏聚力不变.

9) 内摩擦角先减后增，黏聚力不变.

显然上述9种情形不能完全包含2个强度参数指标可能出现的所有变化范围.此外，还有黏聚力和摩擦角可能呈抛物线形式同时变化，发生先增后减、先减后增或呈指数函数形式变化等等，这些变化形式应当结合工程实际进行具体分析，现仅讨论以上9种情形所产生的影响，对于其他变化情况可参照同样的分析方法.

3.2   数值计算参数

为了便于分析，将数值结果与试验数据进行比较，选用湛江黏土在160 kPa偏应力下初始蠕变阶段和稳定蠕变阶段试验数据，采用Burger三维本构关系拟合相关参数，拟合方法与前述相同，强度参数已由试验测定得到，初始黏聚力可选取为38 kPa，内摩擦角为25°^[17]，拟合出相关参数见表 1.

表  1  160 kPa下Burger/Cvisc模型相关参数

Table  1.  Relevant parameters of the Burger/Cvisc model under 160 kPa

σΙ/kPa	K/kPa	GM/kPa	ηM/(kPa·h-1)	GK/kPa	ηK/(kPa·h-1)	R
160	39 018	18 008	1 858	12 963	227	0.997 2

下载: 导出CSV 

| 显示表格

3.3   9种情形分析比较

1) 强度参数不发生变化，单纯Cvisc模型进行数值计算，其结果与三轴蠕变试验数据对比见图 3.

图  3  Cvisc模型数值计算结果与三轴蠕变试验数据对比

Figure  3.  Comparison of numerical results of Cvisc model with test data of triaxial creep

下载: 全尺寸图片 幻灯片

从图 3可以看出：当依据在偏应力160 kPa下前2阶段试验数据拟合出的参数，在强度参数不变的情况下，整个蠕变曲线是趋向于稳定，没有发生加速蠕变现象，且2条曲线在初始蠕变阶段和稳定蠕变阶段基本一致.

2) 强度参数按情形2、3变化，采用一次函数y=at+b进行描述，a为黏聚力变化率，kPa/h，b为初始黏聚力38 kPa；t为蠕变时间，蠕变步长设为5×10^-5；虽然蠕变是一个历时很长的过程，但试验状态下的加速蠕变历时很短，湛江黏土的加速蠕变时间只能记录到0.55 h，与此对应，计算模拟的总蠕变时间定为0.5 h，内摩擦角25°.情形2、3分别是黏聚力下降和增加，而为了更直观表明黏聚力变化对加速蠕变阶段出现时间的影响，将上述2种情况依据黏聚力的变化速率快慢又各分为4小类，具体如下：

① 黏聚力随蠕变时间不断降低，摩擦角不变；下降情况见图 4，所得数值计算结果见图 5.

图  4  黏聚力递减曲线

Figure  4.  Decreasing cohesion vs. time curves

下载: 全尺寸图片 幻灯片

图  5  黏聚力递减蠕变曲线

Figure  5.  Creep curves with decreasing cohesion

下载: 全尺寸图片 幻灯片

由图 5可以看出，在变强度参数模型中，随着给定黏聚力降低幅度的增大，试样会逐渐出现加速蠕变现象，且加速蠕变量随着黏聚力的降低而不断增大.另外，如果将稳定蠕变阶段的延长线与加速蠕变阶段切线的交点称之为“突变点”，意即由稳定蠕变阶段向加速蠕变阶段转化的明显转折点，那么从图中可以看出随着黏聚力降低速率进一步增大，该“突变点”不断左移，表明黏聚力降低，保持其他因素不变，土体在相同条件下更容易出现加速蠕变达到破坏.最后，从图中显而易见，在降低速率为70×10³时，该数值曲线与三轴试验曲线基本一致.

② 黏聚力随蠕变时间不断升高，内摩擦角不变；由于c=38 kPa时不会出现加速蠕变，此时继续增加黏聚力，显得毫无意义，所以为了表明黏聚力的增加对加速蠕变的影响，这里将黏聚力初始值设为8 kPa.升高情况如图 6，所得数值计算结果见图 7.

图  6  黏聚力递增曲线

Figure  6.  Increased cohesion vs. time curves

下载: 全尺寸图片 幻灯片

图  7  黏聚力递增蠕变曲线

Figure  7.  Creep curves with increased cohesion

下载: 全尺寸图片 幻灯片

从图 7可知：初始黏滞系数为8 kPa，4条模拟曲线都是“S”型曲线，而不同于“反倒S”型的三轴试验数据的，此时蠕变曲线只包含初始蠕变和稳定蠕变2个阶段；曲线的竖直位置次序表明随着黏聚力增长速率的提高，轴向应变量逐渐减小，稳定蠕变的“稳定点”，也就是曲线由斜向上到水平的拐点逐渐左移.

3) 强度参数按情形4、5变化，同样分为2种情况，采用一次函数y=ct+d；d为初始内摩擦角，c为降低或升高速率，t为蠕变时间.

① 内摩擦角随蠕变时间不断降低，给定初始值为25°，黏聚力不变；变化情况见图 8，所得数值计算结果见图 9.

图  8  内摩擦角递减曲线

Figure  8.  Decreasing internal friction angle vs. time curves

下载: 全尺寸图片 幻灯片

图  9  内摩擦角递减蠕变曲线

Figure  9.  Creep curves with decreasing internal friction angle

下载: 全尺寸图片 幻灯片

由图 9显然可知：内摩擦角对土体的加速蠕变阶段存在影响.影响规律为随着内摩擦角降幅的提高，加速蠕变阶段的出现越发明显，随着内摩擦角降幅的增加，蠕变曲线的“突变点”不断左移，意味着加速蠕变阶段的提早出现.

② 内摩擦角随蠕变时间不断增加，黏聚力不变；同样为了表明内摩擦角增加对加速蠕变阶段的影响，给定内摩擦角初始值为5°.具体变化情况见图 10，所得数值计算结果见图 11.

图  10  内摩擦角递增曲线

Figure  10.  Increased internal friction angle vs. time curves

下载: 全尺寸图片 幻灯片

图  11  内摩擦角递增蠕变曲线

Figure  11.  Creep curves with increasesd internal friction angle

下载: 全尺寸图片 幻灯片

图 11表明：在内摩擦角增幅最低时，c=10，得到最早加速蠕变曲线；当增幅达到c=20时，得到稳定蠕变曲线，未出现加速蠕变阶段，所以随着内摩擦角增幅逐渐增大，蠕变曲线会逐渐变缓，导致对应的变形量明显降低，加速蠕变阶段不再显现.

4) 黏聚力按情形6变化，内摩擦角不变.初始黏聚力设为5 kPa，按照先增后减变化规律，分为4种变化幅度，采用二次函数(抛物线形式)来描述黏聚力的变化情况，具体见图 12，依据变化情况所得数值计算曲线见图 13.

图  12  黏聚力先增后减曲线

Figure  12.  First increaseing and then decreasing cohesion vs. time curves

下载: 全尺寸图片 幻灯片

图  13  黏聚力先增后减蠕变曲线

Figure  13.  Creep curves with first increasing and then decreasing cohesion

下载: 全尺寸图片 幻灯片

由图 13可以看出：在黏聚力先增后减变化幅度最小时，出现最大加速蠕变，随着黏聚力变化幅度的增大，蠕变曲线的“突变点”逐渐右移，加速蠕变阶段出现时间推后且对应的蠕变量明显降低.

5) 黏聚力按情形7变化，内摩擦角不变.根据黏聚力先减后增的变化幅度不同，分为4种情况，同样采用二次函数(抛物线形式)来进行描述，具体见图 14，依据变化情况所得数值计算曲线如图 15.

图  14  黏聚力先减后增曲线

Figure  14.  First decreasing and then increasing cohesion vs. time curves

下载: 全尺寸图片 幻灯片

图  15  黏聚力先减后增蠕变曲线

Figure  15.  Creep curves with first decreasing and then increasing cohesion

下载: 全尺寸图片 幻灯片

在图 15中，曲线为“S”型曲线，黏聚力先减后增变化幅度较大时，轴向应变量较大，随着幅度的降低，“稳定点”逐渐左移.

6) 内摩擦角按情形8变化，黏聚力不变.内摩擦角根据先增后减的变化幅度不同，分为4种情况，也是采用二次函数(抛物线形式)来进行描述，具体见图 16，依据变化情况所得数值计算曲线见图 17.

图  16  内摩擦角先增后减曲线

Figure  16.  First increasing and then decreasing internal friciton angles vs. time curves

下载: 全尺寸图片 幻灯片

图  17  内摩擦角先增后减蠕变曲线

Figure  17.  Creep curves with first decreasing and then increasing internal friction angle

下载: 全尺寸图片 幻灯片

在图 17中，4条模拟曲线都是“反倒S”型，“突变点”的水平位置在前3条上并没有明显变化，轴向应变最大的第4条更多的是受内摩擦角偏小的影响.

7) 内摩擦角按情形9变化，黏聚力不变.内摩擦角根据先减后增的变化幅度不同，分为4种情况，同样也是采用二次函数(抛物线形式)来进行描述，具体见图 18，依据变化情况所得数值计算曲线如图 19：

图  18  内摩擦角先减后增曲线

Figure  18.  First decreasing and then increasing internal friction angle vs. time curves

下载: 全尺寸图片 幻灯片

图  19  内摩擦角先减后增蠕变曲线

Figure  19.  Creep curves with first decreasing and then increasing internal friction angle

下载: 全尺寸图片 幻灯片

4条模拟曲线都是“S”型，与试验曲线的“反倒S”型不同.同样的，随先减后增幅度的增大，轴向应变量增大.

综合对上述5种情形的分析，可以得出如下结论：

1) 黏聚力和内摩擦角对土体的加速蠕变阶段存在显著的影响.

2) 黏聚力或内摩擦角增加，会使得土体难以进入加速蠕变阶段；反之，土体会提前进入加速蠕变阶段.

3) 从模拟曲线与三轴试验曲线的变化趋势的对比来看，黏聚力和内摩擦角按线性减小、黏聚力按抛物线先增后减与内摩擦角先增后减时，模拟曲线与试验曲线变化规律能够保持一致.

4) 模拟曲线与三轴试验曲线的变化趋势最不一致的是：黏聚力和内摩擦角按线性减增加、黏聚力按抛物线先减后增与内摩擦角先减后增时，模拟曲线以“S”型曲线为主，而试验曲线是“反倒S”型.

4.   结论

1) 本文从采用Cvisc模型无法描述湛江黏土加速蠕变阶段，且依据对湛江黏土的长期和瞬时强度指标的差异来看，引入变强度参数模型，讨论强度参数变化对土体加速蠕变阶段的影响.然后将强度参数变化分为9种情形，采用FISH语言定义变强度参数模型进行数值计算，将模拟曲线与三轴试验曲线的变化趋势进行对比，结合试验研究得出的湛江黏土内摩擦角变化不大的结论，可以得出湛江黏土在能导致加速蠕变的高偏应力条件下的强度参数变化规律：以黏聚力逐渐减小或先有所增大后减小为主要变化趋势.

2) 综合而言，本文是通过数值模拟探讨强度参数的变化规律对于高偏应力条件下蠕变变形的影响，所做研究工作对于探索加速蠕变的内部机理，具有较为重要的意义和价值.