硬岩的复合黏弹塑性流变模型

熊良宵^{1, 2}，杨林德²，张尧²

(1. 宁波大学 建筑工程与环境学院，浙江 宁波，315211；

2. 同济大学 岩土及地下工程教育部重点实验室，上海，200092)

摘  要：基于应变软化模型得到塑性元件，将该塑性元件与六元件黏弹性流变模型组合得到适合于硬岩的复合黏弹塑性流变模型。采用FLAC^3D所提供的二次开发接口程序，开发复合黏弹塑性流变模型的数值程序，对锦屏二级水电站的绿片岩单轴压缩蠕变试验曲线进行拟合，采用拟合得到的参数进行计算以验证接口程序的正确性，并对锦屏二级水电站辅助交通洞进行工程流变计算分析。研究结果表明：复合黏弹塑性流变模型具有程序可置入性，采用该模型进行计算时，计算值与理论值完全一致，从而验证了接口程序的正确性；辅助交通洞开挖后，围岩体的最大变形出现在边墙位置，随着时间的增加，边墙位置处围岩体内的位移和最大主应力呈增加的趋势。

关键词：六元件黏弹性流变模型；复合黏弹塑性流变模型；FLAC^3D；绿片岩；单轴压缩蠕变试验

中图分类号：TU452          文献标志码：A         文章编号：1672-7207(2010)04-1540-09

Composite viscoelasto-plastic rheological model for hard rock

XIONG Liang-xiao^{1, 2}, YANG Lin-de², ZHANG Yao²

(1. Faculty of Architectural, Civil Engineering and Environment, Ningbo University, Ningbo 315211, China;

2. Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University,
Shanghai 200092, China)

Abstract: A plastic component was built based on strain softening model, and composite viscoelasto-plastic model for hard rock was obtained as the combination of six-component viscoelastic rheological model with this plastic component. By adopting the secondary development routine interface of FLAC^3D, the numerical analysis codes of composite viscoelasto-plastic model were developed, and calculation was done by using the identified parameters of uniaxial compression creep test curve of green schist specimen from Jinping Second Stage Hydropower Station. In the end, the numerical analysis was conducted for auxiliary tunnel at Jinping Second Stage Hydropower Station. The research results show that numerical codes are easy to be embedded into composite viscoelasto-plastic model, and the interface codes for this model are verified by the agreement of calculation result with theoretical result. The distortion of surrounding rock near the side wall is the largest after auxiliary tunnel is excavated. The displacement and maximum principal stress of surrounding rock near the side wall increase with time.

Key words: six-component viscoelastic rheological model; composite viscoelasto-plastic rheological model; FLAC^3D; greens chist; uniaxial compression creep test

岩石流变力学本构模型一直是岩石力学领域的研究重点，迄今为止，有关岩石流变力学本构模型的研究成果较多，包含经验流变模型、元件流变模型、损伤断裂流变模型和弹/黏塑性模型等^[1-4]。从近几年的相关研究成果看，研究者普遍注重建立可以反映岩石加速蠕变的非线性流变元件模型^[5-10]，且对Bingham体中的线性黏滞体进行改进，使其成为能反映加速蠕变的非线性黏滞体，再将改进的Bingham体与其他线性黏弹性流变元件模型组合得到非线性黏弹塑性流变模型。曹树刚等^[5-10]提出的这类模型概念比较清楚，能通过对蠕变全程试验曲线进行拟合得到可以明确反映加速蠕变的相关参数。这类模型将塑性屈服元件与可反映加速蠕变的非线性黏滞体进行并联，在程序编制时，塑性屈服准则一般采用Mohr-Coulomb或者Drucker-Prager准则，这就意味着如果岩石受到的外荷载过大，将导致反映加速蠕变的非线性黏滞体起作用，岩石的蠕变特性出现加速蠕变阶段，因此，能否出现加速蠕变完全取决于岩石所受到的外荷载。然而，有些硬脆岩的蠕变试验曲线往往无明显加速阶段，而是由衰减蠕变阶段或者稳态蠕变阶段直接进入破坏阶段，同时，即使岩石受到的外荷载小于其长期强度时，也有可能存在塑性变形，即往往黏弹性和黏塑性变形特性共存。因此，许多研究者采用FLAC^3D的Cvisc模型来分析岩石黏弹性与黏塑性共存的流变特性，该模型是基于Mohr-Coulomb(M-C)模型建立塑性元件，并与Burgers模型组合得到^[11]。陈炳瑞等^[12]则提出将广义Burgers模型与M-C塑性元件组合得到黏弹塑性普适模型。但对深部硬脆岩，采用M-C模型并不合适，因此，本文作者基于M-C应变软化模型得到塑性元件，将该塑性元件与六元件线性黏弹性模型组合得到复合黏弹塑性流变模型，该模型可考虑强度参数随塑性应变增加而弱化的特性。同时，按照UDM接口计算程序的要求，将复合黏弹塑性流变模型置入FLAC^3D中，对锦屏二级水电站辅助交通洞绿片岩进行了单轴压缩蠕变特性试验，采用拟合试验曲线得到的参数进行计算以验证接口程序的正确性，并以锦屏二级水电站辅助交通洞为例进行工程流变数值分析。

1  复合黏弹塑性流变模型

1.1 六元件线性黏弹性流变模型

当硬岩在低应力水平下表现为黏弹性蠕变特性，且出现了衰减和稳态蠕变阶段时，若采用Burgers模型进行分析会存在一定的误差，原因是Burgers主要用于软岩的黏弹性蠕变特性分析^[13]。相对而言，广义Burgers模型的适用范围比Burgers模型的适用范围要广^[14]，因此，本文作者采用广义Burgers模型中的六元件黏弹性流变模型来分析硬岩的衰减和稳态蠕变阶段。六元件黏弹性流变模型为 Kelvin模型和Burgers模型的组合(简称为KB模型)，如图1所示。其中：为受到的应力；E₁，E₂和E₃分别为Maxwell体、第1个Kelvin体和第2个Kelvin体的材料弹性参数；，和分别为Maxwell体、第1个Kelvin体和第2个Kelvin体的材料黏性参数。

图1  六元件黏弹性模型

Fig.1  Six-component viscoelastic rheological model

六元件流变模型的一维蠕变方程为：

             (1)

式中：为应变。

三维应力状态下的偏应变为：

             (2)

且

,           (3)

式中：为偏应变；为偏应力；G₁，G₂和G₃为三维剪切模量；K为体积模量；为泊松比。

式(2)两边同时加上球应变张量后，可以得到岩石的蠕变应变。在三维应力状态下，六元件黏弹性流变模型的轴向蠕变方程为：

        (4)

式中：，和分别为第1主应力、第2主应力和第3主应力；为第1主应变。

1.2 复合黏弹塑性流变模型

参考FLAC^3D中的Cvisc模型^[11]，将基于Mohr-Coulomb模型得到的塑性元件与六元件黏弹性流变模型组合，得到适合于硬岩的复合黏弹塑性流变模型，模型如图2所示，该模型是陈炳瑞等^[12]所提出的黏弹塑性普适模型中的一种。

图2  复合黏弹塑性流变模型

Fig.2  Composite viscoelasto-plastic rheological model

 (1) 当＜时，该模型则退化为六元件黏弹性流变模型。

(2) 当≥时，复合黏弹塑性流变模型总的偏应变速率可表示为：

           (5)

式中：为偏应变速率；上标K1，K2，M和P分别代表第1个Kelvin体、第2个Kelvin体、Maxwell体和Mohr-Coulomb体对应力应变的贡献，其中2个Kelvin体的应变速率关系为：

            (6)

            (7)

式中：和为三维应力状态下第1个Kelvin体的剪切模量和黏滞系数；和为三维应力状态下第2个Kelvin体的剪切模量和黏滞系数。

Maxwell体和Mohr-Coulomb塑性体塑性应变速率，以及Mohr-Coulomb的屈服函数见文献[11]。

由于图2所示模型可以看成由Burgers体、Kelvin体和M-C塑性体共同组成，故可命名为KBMC模型。

1.3 修正复合黏弹塑性流变模型

对于硬脆岩体，尤其是深部硬岩体，采用Mohr-Coulomb应变软化模型(Strain softening model)比Mohr-Coulomb模型更适合。Mohr-Coulomb应变软化模型是Mohr-Coulomb模型的一种特殊形式。在应变软化模型中，由预先定义的硬化参数，根据分段线性原则，在塑性应变产生后，部分或所有单元的强度参数都可发生变化^[11]，如图3所示。其中：为塑性剪切应变；为黏聚力；为摩擦角。

图3  黏聚力、摩擦角与塑性应变关系曲线

Fig.3  Cohesion-plastic strain curve and friction-plastic strain curve

在实际岩石的蠕变试验过程中，c和等强度参数随着塑性应变的增加不断衰减。因此，可基于应变软化模型得到新的塑性元件(SS元件)，并替换图2中的塑性屈服元件，得到修正复合黏弹塑性流变模型，如图4所示。由于该模型可以看成由Burgers体、Kelvin体和SS塑性体共同组成，故可命名为KBSS模型。

图4  修正复合黏弹塑性流变模型

Fig.4  Modified composite viscoelasto-plastic rheological model

应变软化模型(Strain-Softening model)中的破坏准则、屈服函数、势函数、塑性流动准则与应力校正等都与相应的Mohr-Coulomb模型一致。

2  复合黏弹塑性流变模型程序开发

六元件黏弹性流变模型的偏应变增量表达式为：

          (8)

式中：，和分别为Maxwell体、第1个Kelvin体和第2个Kelvin体的偏应变增量。

据文献[11]可得到六元件黏弹性流变模型新偏应力的表达式为：

         (9)

其中：

；      (10)

；     (11)

    (12a)

    (12b)

和分别为新偏应力和旧偏应力；和分别为三维应力状态下第1个Kelvin体的剪切模量和黏滞系数；和分别为三维应力状态下第2个Kelvin体的剪切模量和黏滞系数；和分别为三维应力状态下Maxwell体的剪切模量和黏滞系数；和分别表示第1个和第2个Kelvin体的旧偏应变；为时间增量。

复合黏弹性塑性流变模型新偏应力的表达式为：

  (13)

式中：为塑性偏应变增量。

可根据上述六元件黏弹性流变模型和复合黏弹塑性流变模型的应力增量表达式，采用VC6.0编写其程序代码，通过UDM接口程序嵌入FLAC^3D中。

3  程序验证

3.1  单轴压缩蠕变试验分析

单轴压缩蠕变试验所采用的绿片岩取自锦屏二级水电站辅助交通洞B洞西端6号横通道，具片状构造，常具有灰白色大理岩条带及透镜体，岩块也较硬，在干燥状态下平均单轴抗压强度为60~70 MPa，长期强度为40~50 MPa，岩样尺寸(长×宽×高)为100 mm×100 mm×100 mm，蠕变试验曲线如图5所示。

由图5可知：在第4级应力(48 MPa)作用下，蠕变试验曲线由稳态蠕变阶段直接进入破坏阶段，这属于蠕变韧-脆破坏。

图5  绿片岩的单轴压缩蠕变试验曲线

Fig.5  Uniaxial compression creep test curve of green schist

采用六元件黏弹性流变模型的三维蠕变方程对试验曲线进行拟合时，将泊松比取为0.25，同时在确定每级的瞬时弹性变形时，扣除前面分级加载累计的蠕变变形量。采用六元件黏弹性流变模型进行拟合，结果如图6所示，拟合得到的参数见表1。

图6  六元件流变模型拟合结果

Fig.6  Fitting results based on six-component rheological model

表1  六元件流变模型拟合参数值

Table 1  Rheological parameters on base of six-component rheological model

由图6和表1可知：采用六元件流变模型对试验曲线进行拟合时，拟合曲线和试验曲线较吻合，且拟合的相关系数均在0.990以上，这表明采用六元件黏弹性流变模型分析硬脆岩的黏弹性蠕变特性是合适的。

3.2  接口程序验证

采用FLAC^3D计算时，模型尺寸(长×宽×高)也为100 mm×100 mm×100 mm，模型如图7所示。在计算时观察上表面中心位置A点的竖向位移与时间的关系曲线。

以第1级试验曲线为例，采用六元件流变模型进行计算，将计算得到的A点轴向位移除以试件的高度0.1 m得到该点的轴向应变，并与试验曲线的拟合值进行对比，结果如图8所示。

由图8可知：采用六元件黏弹性流变模型计算时，在初始阶段均有突变，但很快达到稳定，计算得到的轴向应变完全与拟合值一致，验证了KB模型接口程序的正确性。

当采用修正复合黏弹塑性流变模型计算时，将黏聚力和抗拉强度取为1×10²⁰ Pa，不考虑其塑性变形时，得到的结果与采用六元件黏弹性流变模型计算得到的结果完全一致。这说明在没有塑性变形的情况下，采用修正复合黏弹塑性流变模型计算时可以退化为六元件黏弹性流变模型，从而验证了修正复合黏弹塑性流变模型接口程序的正确性。

图7  FLAC^3D计算模型

Fig.7  Computing model of FLAC^3D

图8  基于六元件流变模型的计算结果

Fig.8  Computing result on base of six-component rheological model

4  工程应用

4.1  锦屏二级水电站辅助洞工程简介

锦屏二级水电站位于四川省凉山洲境内的雅砻江干流上锦屏大河弯处，共有4条引水隧洞，洞线平均长度约16.67 km，上覆岩体一般埋深1.5~2.0 km，最大埋深约为2.525 km。辅助洞位于引水隧洞线的南侧，全长约17.5 km，由2条平行的单车道隧道即Ⅰ洞和Ⅱ洞组成，两洞中心距为35 m。引水隧洞和辅助洞均具有埋深大、洞线长和洞径大的特点。

4.2  计算模型和边界条件

现以岩样采集点Ⅱ洞为例进行工程流变数值计算。在辅助洞实际施工过程中，Ⅱ洞开挖通常超前Ⅰ洞，且两洞中心间距达35 m，Ⅰ洞的开挖过程对Ⅱ洞的影响几乎可以不予考虑，因此，本文只是按照Ⅱ洞的断面尺寸来建立计算模型。Ⅱ洞的断面型式为马蹄形，断面计算宽度取6 m，高度取7 m，计算模型如图9所示，其中：AB，AC，AD，AE，BC和DE均为监测线。

图9  Ⅱ洞计算模型

Fig.9  Computational model of tunnel Ⅱ

整个模型水平方向(x轴方向)和竖直方向(z轴方向)各取80 m，沿隧洞轴线(y轴方向)取单位厚度。所建模型共1 504个单元，3 046个节点。计算埋深取为1.5 km，其自重应力约为42 MPa，根据地应力测试报告，西端辅助洞第4横通道的为42.33 MPa，为26.06 MPa，为24.48 MPa。

4.3  参数设置

采用修正复合黏弹塑性流变模型计算时，岩石的瞬时弹性模量E取18 GPa，泊松比取0.25，密度为    2 800 kg·m^-3。体积模量K取为12 GPa，剪切模量G₁为7.2 GPa。根据绿片岩的三轴强度试验可确定计算时初始黏聚力c取8 MPa，初始内摩擦角φ取50°，抗拉强度取2 MPa，黏聚力c和内摩擦角φ随塑性应变增加而衰减的关系如表2所示。

黏弹性参数取5×10⁴ GPa·h，G₂取400 GPa，取20 GPa·h，G₂取400 GPa，取800 GPa·h。现场隧洞基本未进行支护，因此，计算时不考虑设置支护，隧洞全断面一次性开挖后进行流变计算，计算时间为30 d。

表2  c,φ与塑性应变ε^P的关系

Table 2  Relationship among c,φ and ε^P

4.4  计算结果分析

计算得到各监测线收敛位移及变形速率如图10  所示。由图10可知：30 d后，监测线DE和BC的收敛位移最大，监测线AB和AC的收敛位移最小，这说明边墙是最容易发生变形的部位；在10 d时基本达到稳态蠕变阶段，各监测线的收敛变形速率为0.04~   0.08 mm/d，围岩体的变形已趋于稳定。

Fig.10  Displacement and deformation rate convergence curves

不同计算时间得到的位移矢量场如图11所示。由图11可知：隧洞刚开挖完时，最大位移发生在左、右边墙围岩体中，之后，随着时间的增加，最大位移呈增加的趋势。

不同计算时间得到的最大、中间和最小主应力场分别如图12~14所示。

由图12~14可知：围岩主应力也具有时效特性，如隧洞刚开挖完时，主要在边墙两侧围岩体中出现最大主应力集中现象，而随着时间的增加，出现最大主应力集中的位置逐渐向岩体深部位置转移，且应力集中区的最大主应力随着时间的增加也呈增加的趋势。

不同计算时间得到的塑性屈服状态如图15所示。图15中：none表示该区域一直处于弹性变形状态；shear-n和tension-n表示该区域当前发生剪切和拉伸屈服状态；shear-p和tension-p表示该区域之前发生了剪切和拉伸屈服状态，后来由于应力重分布，这一区域退出了塑性屈服状态。

由图15可知：当隧洞刚开挖完时，左、右边墙位置围岩体出现了塑性屈服，且两侧的塑性屈服范围并不完全对称，这是初始应力不对称造成的；随着时间的增加和应力的调整，左、右边墙位置围岩体中塑性屈服的范围和深度逐渐增大，且两侧的塑性屈服范围逐渐趋于对称，但超过一定时间后，塑性屈服的范围和深度增加的幅度很小，变形基本达到稳定。在整个计算过程中，只有仰拱位置围岩体中的个别单元出现了拉伸屈服，其余则均为剪切屈服。

图11  位移矢量场

Fig.11  Displacement vector field

图12  最大主应力场

Fig.12  Major principal stress

图13  中间主应力场

Fig.13  Middle principal stress

图14  最小主应力场

Fig.14  Minor principal stress

图15  塑性屈服状态

Fig.15  Plastic yield state

5  结论

(1) 基于Mohr-Coulomb应变软化模型得到塑性元件，将该塑性元件与六元件黏弹性流变模型组合得到适合于硬岩的KBSS复合黏弹塑性流变模型，该模型可以反映强度参数随塑性应变的增加而衰减的特性。

(2) 六元件黏弹性流变模型和KBSS复合黏弹塑性流变模型均具有程序可置入性。对锦屏二级水电站辅助交通洞绿片岩进行了单轴压缩蠕变特性试验，采用拟合试验曲线得到的参数进行计算，计算值与理论值一致，验证了接口程序的正确性。

(3) 锦屏二级水电站辅助交通洞开挖后，边墙位置水平测线之间的收敛位移最大，且随着时间的增加，边墙位置围岩体的位移、应力集中区的最大主应力和塑性屈服范围也呈增加的趋势，均具有明显的时效特性，但超过10 d时，围岩体的变形开始趋于稳定。

参考文献：

[1] 孙钧. 岩石流变力学及其工程应用研究的若干进展[J]. 岩石力学与工程学报, 2007, 26(6): 1081-1106.
SUN Jun. Rock rheological mechanics and its advance in engineering applications[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(6): 1081-1106.

[2] 徐卫亚, 杨圣奇, 褚卫江. 岩石非线性黏弹塑性流变模型(河海模型)及其应用[J]. 岩石力学与工程学报, 2006, 25(3): 433-447.
XU Wei-ya, YANG Sheng-qi, CHU Wei-jiang. Nonlinear viscoelasto-plastic rheological model (Hohai model) of rock and its engineering application[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(3): 433-447.

[3] 范庆忠, 高延法, 崔希海, 等. 软岩非线性蠕变模型研究[J]. 岩土工程学报, 2007, 29(4): 505-509.
FAN Qing-zhong, GAO Yan-fa, CUI Xi-hai, et al. Study on nonlinear creep model of soft rock[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(4): 505-509.

[4] Cristescu N D. A general constitutive equation for transient and stationary creep of rock salt[J]. International Journal of Rock Mechanics Sciences and Geomechanics Abstracts, 1993, 30(2): 125-140.

[5] 曹树刚, 边金, 李鹏. 岩石蠕变本构关系及其改进的西原正夫模型[J]. 岩石力学与工程学报, 2002, 21(5): 632-634.
CAO Shu-gang, BIAN Jin, LI Peng. Rheological constitutive relationship of rocks and a modificial model[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(5): 632-634.

[6] 陈沅江, 潘长良, 曹平, 等. 一种软岩流变模型[J]. 中南工业大学学报: 自然科学版, 2003, 34(1): 16-20.
CHEN Yuan-jiang, PAN Chang-liang, CAO Ping, et al. A new rheological model for soft rocks[J]. Journal of Central South University of Technology: Natural Science, 2003, 34(1): 16-20.

[7] 宋飞, 赵法锁, 卢全中. 石膏角砾岩流变特性及流变模型研究[J]. 岩石力学与工程学报, 2005, 24(15): 2569-2664.
SONG Fei, ZHAO Fa-suo, LU Quan-zhong. Study on rheological properties and model for gypsum breccias[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(15): 2569-2664.

[8] 徐卫亚, 杨圣奇, 谢守益, 等. 绿片岩三轴流变力学特性的研究(Ⅱ): 模型分析[J]. 岩土力学, 2005, 26(5): 693-698.
XU Wei-ya, YANG Sheng-qi, XIE Shou-yi, et al. Investigation on triaxial rheological mechanical properties of greenschist specimen (Ⅱ): Model analysis[J]. Rock and Soil Mechanics, 2005, 26(5): 693-698.

[9] 赵延林, 曹平, 文有道, 等. 岩石弹黏塑性流变试验和非线性流变模型研究[J]. 岩石力学与工程学报, 2008, 27(3): 477-486.
ZHAO Yan-lin, CAO Ping, WEN You-dao, et al. Viscoelasto-plastic rheological experimental and nonlinear rheological model of rocks[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(3): 477-486.

[10] 蒋昱州, 张明鸣, 李良权. 岩石非线性黏弹塑性蠕变模型研究及其参数识别[J]. 岩石力学与工程学报, 2008, 27(4): 832-839.
JIANG Yu-zhou, ZHANG Ming-ming, LI Liang-quan. Study on nonlinear viscoelasto-plastic creep model of rock and its parameter identification[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(4): 832-839.

[11] Itasca Consulting Group Inc. FLAC-3D (Fast Lagrangian Analysis of Continua in 3 Dimensions), Version3.00, Users Manual[R]. USA: Itasca Consulting Group Inc, 2005.

[12] 陈炳瑞, 冯夏庭. 黏弹塑性普适组合模型及其工程应用[J]. 岩石力学与工程学报, 2008, 27(5): 1028-1035.
CHEN Bing-rui, FENG Xia-ting. Universal viscoelastic combination model and its engineering applications[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(5): 1028-1035.

[13] 陶波, 伍法权, 郭改梅, 等. 西原模型对岩石流变特性的适应性及其参数确定[J]. 岩石力学与工程学报, 2005, 24(17): 3165-3171.
TAO Bo, WU Fa-quan, GUO Gai-mei, et al. Flexibility of visco-elastoplastic model to rheological characteristic of rock and solution of rheological parameter[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(17): 3165-3171.

[14] 周德培. 岩石流变性态研究[D]. 成都: 西南交通大学土木工程学院, 1986: 41-46.
ZHOU De-pei. The study of rock rheology[D]. Chengdu: Southwest Jiaotong University，College of Civil Engineering, 1986: 41-46.

收稿日期：2009-08-03；修回日期：2009-11-10

基金项目：国家自然科学基金、雅砻江水电开发联合研究基金重点资助项目(50639090)

通信作者：熊良宵(1982-)，男，江西九江人，博士，讲师，从事地下工程和岩石力学方面的研究；电话：0574-87600316；E-mail: xiongliangxiao@163.com

(编辑 刘华森)