李修磊，李起伟,2，李 倩

(1.重庆交通大学交通运输学院，重庆 400074；2.中交第一公路勘察设计研究院有限公司，陕西 西安 710075；3.四川大学锦江学院，四川 眉山 620860)

随着深部资源的开采(如油井、气井、矿井、地热等)，深埋硐室开挖以及高放射核废料深埋处置等深部地下工程已逐渐向地下延伸数千米[1-2]，这些工程往往处于高温、高应力环境，尽管初始阶段呈稳定状态，但由于长时间岩石的蠕变作用可能导致围岩变形发生破坏。因此，研究岩石在长时间温度-应力耦合作用下的蠕变力学行为和变形规律，对工程的稳定性具有重要的意义[3-4]。

关于岩石温度效应的蠕变性能研究，已有相关试验开展。Chopra等[5]采用高分辨率气体介质的试验装置，研究纯橄榄石在围压为300 MPa、2种温度分别为1 100℃和1 300℃下的蠕变特性，试验发现高温、高应力下橄榄石的蠕变行为可用Burgers模型进行描述。Kinoshitan等[6]通过试验研究了花岗岩在20～100℃下单轴蠕变特性，结果表明温度的升高能够加速花岗岩的蠕变破坏。Chen等[7]单轴和三轴蠕变试验结果表明，花岗岩的时效破坏与其断裂过程所引起的损伤演化密切相关；相同加载应力下，围压的增加能够降低应变速率，延缓破坏时间；相同围压条件下，温度和应力水平的升高不会对花岗岩破坏时裂纹扩展形态有显著影响，但会加快裂纹扩展的速度，进而缩短蠕变破坏时间。刘泉声等[8]针对三峡花岗岩开展了20～300℃范围内的单轴和三轴蠕变试验，给出了花岗岩蠕变特性随温度的变化规律。张强勇等[9]研究了不同温度下片麻状花岗岩的三轴蠕变特性，试验结果表明，片麻状花岗岩存在蠕变应力阈值，且受温度和围压的双重影响；中、低应力条下，岩石只发生减速和等速蠕变；高应力条件下岩石会出现加速蠕变过程；温度升高会导致岩石的蠕变速率增大，应力阈值越低，蠕变破坏时间越短。张宁等[10]和陈亮等[11]分别对鲁灰花岗岩和北山花岗岩的蠕变特性进行了试验研究，温度和应力水平都会影响岩石的蠕变特性。上述试验成果为建立岩石的蠕变本构关系奠定了良好基础。

目前，国内外已有学者考虑温度效应的岩石蠕变本构关系理论研究。如，Chen[12]等采用西原模型对花岗岩蠕变特性进行了描述，分析了加速蠕变阶段损伤演化规律，建立了各模型参数与温度之间函数关系；唐皓等[13]和曹丽丽等[14]分别基于函数阶微积分理论建立了适用于盐岩和泥岩的函数阶蠕变本构模型。胡其志等[15]基于广义Bingham模型在蠕变衰减和稳态蠕变阶段引入非线性函数，对加速蠕变过程进行损伤劣化分析，建立了考虑温度的盐岩损伤蠕变本构模型；王春萍等[16]建立了温度-应力耦合作用下的西原模型，并用试验结果进行了验证；张强勇等[17]基于不同温度、不同应力水平下片麻花岗岩的三轴蠕变试验结果，建立了热力耦合作用的热黏弹塑性损伤蠕变模型，提出了高、中、低应力状态下温度对片麻花岗岩蠕变特性的影响规律；梁玉雷等[18]考虑温度变化周期对岩石蠕变的影响，建立了适合描述大理岩的热力耦合的改进Burgers模型；Xu等[19]基于热力学原理，建立了温度-应力耦合作用下脆性岩石的蠕变损伤模型，并以最大拉应力和摩尔库伦准则作为单元破坏准则，在有限元软件Comsol的基础上对模型进行二次开发，较为准确地模拟了不同温度条件下花岗岩典型的三阶段蠕变变形全过程。

由以上分析可知，现有考虑温度作用的岩石蠕变模型大多是在经典蠕变模型的基础上，根据不同温度的试验结果建立模型参数与温度的函数关系，进而提出了考虑温度效应的岩石蠕变本构关系。根据文献[11]，岩石的蠕变破坏与裂纹起裂、扩展引起的损伤变化有关，而现有蠕变模型中并没有考虑外力作用下岩石内部裂纹扩展引起的损伤演化。为此，本文将基于岩石断裂理论，考虑裂纹扩展引起的损伤，推导热力耦合作用下的一维损伤蠕变方程，并将其推广到三维应力状态，采用非线性最小二乘法确定模型参数，进而利用不同温度下花岗岩的蠕变试验结果验证本文模型的适用性和合理性。

1 应力与裂纹扩展之间的关系

图1给出了脆性岩石单元体微裂纹扩展的细观力学模型[20]，假定岩石为各向同性弹性材料，单个初始裂纹的长度为a，初始裂纹与主应力σ1方向的夹角为β，岩石所受围压σ2=σ3。

图1 压力作用下裂纹扩展力学模型

利用莫尔圆(图2)，得到应力作用下裂纹表面的正应力σn和剪应力τ分别为：

(1)

(2)

当外部应力大于临界应力σ1C(也就是说，裂纹表面的剪应力τ大于材料固有的抗剪切强度τ0)时，翼型裂纹开始产生；当外部应力小于临界应力σ1C，没有翼型裂纹的产生。

图2 摩尔-库伦模型

材料固有的抗剪切强度τ0为：

(3a)

对于单轴压缩试验，τ0为：

(3b)

式中：φ——内摩擦角；

C0——无侧限抗压强度。

在常规三轴试验中(图1)，由应力分量坐标变换，远场的应力状态为：

(4)

(5)

σs=σ1C-σ3

(6a)

单轴应力状态下的屈服应力为：

(6b)

式中：σs——屈服应力或称为临界损伤应力。

2 岩石热黏弹塑性蠕变损伤模型

根据不同温度下岩石的蠕变试验结果[17]，发现岩石的蠕变过程具有时效性(图3，其中分级加载应力差σ1-σ3分别为160，180，200，220和240 MPa)，具体特征如下：加载过程中岩石会发生瞬时应变增量；低应力水平下岩石发生减速和等速蠕变，变形随时间趋于稳定；高应力水平下(应力超过临界值)时，岩石除了有减速和等速蠕变，还会发生加速蠕变，为典型的非线性蠕变行为。岩石的蠕变特性以及相应的力学参数与其受到的温度和外部应力密切相关。

图3 不同温度下片麻花岗岩的轴向蠕变曲线(σ3=30 MPa)[17]

岩石经典的蠕变模型(如Kelvin模型、Bingham模型、Burgers模型和Nishihara(西原)模型等)仅能反映蠕变变形的瞬时和稳定阶段[21]，无法有效考虑温度和应力耦合作用下岩石蠕变变形的全过程。目前，多数试验结果[7,11,22]证实，加速蠕变阶段的发生主要与岩石内部微裂纹扩展的损伤演化有关。基于此，本文将结合常用的Burgers模型和西原模型引入岩石的断裂损伤，并考虑温度和应力的耦合(图4)。

图4 模型元件的组成

改进后的本文模型可看作是在Burgers模型的基础增加了1个考虑损伤的黏塑性元件，也可看作是在西原模型的基础上增加了1个串联的牛顿黏壶元件。

2.1 热-黏塑性损伤元件

近些年，已有试验研究成果表征了岩石的损伤演化过程[23-24]。结果表明，蠕变试验过程中岩石的损伤演化可以用负指数函数来表征[12,17]：

D=1-exp(-αt)

(7)

式中：D——损伤变量，由0逐渐趋近于1，分别对应着初始和完全损伤状态；

α——与损伤演化过程相关的参数。

本文研究中将采用式(7)来模拟随时间变化的岩石内部裂纹扩展引起的损伤演化过程。

Li等[25]的试验数据证实只有当施加的应力达到一定程度和超过某一阈值时，岩石才会发生损伤演化。此外，Chen等[26]通过电子显微镜对岩石蠕变破坏机理的研究同样证实，只有当施加的应力超过屈服极限时，微裂纹扩展才会加速。在此基础上，本文给出一种新的热耦合元件(即热黏塑性体)，来表征岩石在不同温度下的损伤演化对其蠕变变形的影响(图4c)。试验结果表明，由于微裂纹扩展，岩石的蠕变速率呈不断增加的趋势，尤其是高温环境下这种现象更加明显。该热黏塑性体中的黏壶元件为受温度和应力影响的牛顿体，其本构关系如下：

(8)

式中：σv——作用在黏壶的应力；

η3(T,D)——与温度T和损伤变量D有关的黏滞系数；

εvp——黏塑性应变。

为了考虑温度T和损伤演化对蠕变变形的影响，提出黏度系数η3的表达式：

η3(T,D)=η3(T)(1-D)

(9)

热黏塑性损伤耦合元件的总应力σ为作用在牛顿黏壶上的应力σv与塑性元件上的应力σp之和。因此，塑性元件上的应力σp可表示为：

(10)

式中：σs——临界损伤应力。

σs可用式(6)进行描述。联立公式(7)～(10)可得：

(11)

令折减函数p(t)=exp(-αt)，α取不同值时，p(t)随时间t的变化规律如图5所示。可以看出，折减函数随时间逐渐减小，黏滞系数η3(T,D)因岩石损伤蠕变而逐渐衰减。固定应力水平下，初始条件为时间t= 0时的黏塑性应变εvp= 0。因而，可获得热黏塑性耦合元件的本构关系如下：

(12)

图5 α不同时，p(t)随时间的变化曲线

大量的蠕变试验[27-29]结果表明岩石蠕变变形的稳定阶段具有明显的非线性特征(即不同应力水平下稳定蠕变阶段的应变率与应力的比值不是常数)。宋飞等[27]通过大量的岩石蠕变试验，分析了稳定阶段蠕变速率与应力的关系，提出了非线性的黏滞分量，其中黏滞系数与应力呈指数变化关系。基于此，本文提出一种新的反映临界应力并与温度有关的非线性黏滞分量，本构关系如下：

(13)

式中：η1(T)——热黏弹性体的黏滞系数；

λ——与温度有关的黏性参数；

εv——黏塑性应变。

固定应力水平下的蠕变方程为：

(14)

2.2 热-黏-弹-塑性损伤蠕变本构模型

图6给出了本文模型蠕变变形随时间曲线的示意图。可以看出，相比传统的西原模型，当σ<σs时，稳态阶段的蠕变速率不为0，非线性Newton体的黏滞系数与温度和施加应力的大小有关，而在Burgers模型中Newton体的黏滞系数为常数；当σ≥σs时，考虑了岩石脆性断裂损伤引发的加速蠕变，明显不同于传统的西原模型和Burgers模型。

图6 改进后本文模型的应力-应变关系

图4c中模型的瞬时热弹性模量为E1(T)，热黏弹性模量为E2(T)，热黏滞系数为η1(T)，热黏弹性体的黏滞系数为η2(T)，热黏塑性体的黏滞系数为η3(T,D)，元件中的临界应力为σs。图4c中四部分所受到的应力分别为σe，σv，σve和σvp，对应的应变分别为εe，εv，εve和εvp。本文模型满足以下应力-应变关系：

(15)

整理公式(15)，得到一维的本构关系为：

当σ<σs时，

(16)

当σ>σs时，

(17)

当应力为常数时，可以根据叠加原理求得一维情况下的热黏弹塑性损伤蠕变本构方程：

当σ<σs时，

(18)

当σ>σs时，

(19)

式中各符号的意义均与前述相同。

2.3 热-黏-弹-塑性损伤蠕变本构方程的三维表达式

虽然可视化的物理元件能够方便描述单轴应力状态下的一维蠕变模型，但对于围压、轴压同时存在的三轴应力状态，一维蠕变模型并不适用。为此，通过引入弹、塑性理论，将应力分解为球应力张量和偏应力张量。球应力张量只产生弹性体积应变，而流变部分由偏应力张量产生[21]。

(20)

(1)热弹性体

应力张量σij可分解为球应力张量δijσm和偏应力张量sij，即：

σij=sij+δijσm

(21)

式中δij为Kronecker符号；球应力张量δijσm只改变物体体积不改变物体形状；偏应力张量sij改变物体形状对体积不产生影响；平均正应力σm为：

(22)

应变张量也可分解偏应变张量eij和球应变张量δijεm：

(23)

根据广义Hooker定律，热弹性体的三维本构关系可描述为：

(24)

式中：G，K——剪切模量和体积模量。

根据弹性力学知识，剪切模量G、体积模量K、弹性模量E之间的关系为：

(25)

(26)

(2)热黏性体

对于热黏性体，考虑塑性流动法则，三维应力条件下的热黏弹性应变为：

(27)

式中：F——屈服函数；

F0——初始状态屈服函数值；

Q——塑性势函数。

采用相关流动性法则，F=Q。根据Tresca屈服准则，函数F的表达式为：

(28)

(3)热黏弹性体

(29)

(4)热黏塑性体

考虑塑性流动法则，三维应力状态下热黏塑性应变为：

(30)

式中：F——屈服函数；

F0——初始状态屈服函数值；

Q——塑性势函数；

< >——开关函数。

(31)

根据相关流动性法则，F≥0时，F=Q，对式(30)进行积分，可得：

(32)

联立式(20)、(26)、(27)、(29)和(31)，可得三维应力条件下岩石热黏弹塑性损伤蠕变方程：

(33)

令式(33)中i=j=1，三向应力分别为σ1，σ2，σ3，传统三轴压缩试验中σ2=σ3，有：

(34)

令初始状态屈服函数值F0=1，则有：

(35)

将式(34)、(35)代入式(33)中，整理得到三轴应力状态下的热-力耦合损伤蠕变本构方程：

当σ1-σ3<σs时，

(36)

当σ1-σ3≥σs时，

(37)

3 模型参数的确定

岩石蠕变模型参数确定的方法已有很多，通常采用两种类型：图形拟合法[30]和优化分析法[17,31]。前者是根据蠕变曲线几何形态与蠕变参数物理意义之间的对应关系，通过数据拟合来确定相应的模型参数；后者通常采用回归分析和最小二乘法来确定模型参数。与第一类方法相比，优化分析法能够得到更高的拟合精度，并且由于简单适用性被广泛应用于蠕变试验参数的确定。为此，本文研究中将采用Levenberg-Marquardt非线性最小二乘法，对式(19)和式(37)中的热黏弹塑性损伤蠕变模型的参数进行确定。

为了简化模型参数确定的复杂性，分别对一维应力条件下的蠕变模型公式(19)和三维应力条件下的蠕变模型公式(37)进行简化，形式如下：

+〈E′〉[eα(T)t-1]

(38)

式中：A′，B′，C′，D′，E′——拟合参数。

一维应力条件：

(39)

三维应力条件：

(40)

4 模型验证及参数分析

4.1 模型验证

为了验证本文蠕变模型的合理性，以下将采用2组试验数据(图3中σ3=30 MPa、T=70℃的三轴蠕变试验结果，以及文献[7]中不同温度下的单轴蠕变试验结果)与本文模型的计算结果对比分析，分别如图7和图8所示。模型的有效性取决于对试验数据的拟合程度，本文将采用非线性最小二乘法对模型参数进行确定，通过反演得到的热黏弹塑性损伤蠕变模型的力学参数分别见表1和表2。

图7 本文模型与花岗岩三轴蠕变试验值[17]的对比

由图7和图8可以看出，无论是单轴蠕变试验还是三轴蠕变试验，本文提出的热黏弹塑性损伤蠕变模型均可准确地反映岩石在不同温度、不同应力状态下的蠕变变形的全过程，尤其是对加速蠕变阶段的模拟效果，与试验结果相吻合较好。说明本文建立的热黏弹塑性损伤蠕变模型及通过反演得到的损伤蠕变力学参数是合理可靠的。

表1 围压为30 mPa、温度为70℃时的损伤蠕变参数

表2 不同温度下的损伤蠕变参数(σ3=0 mPa)

图8 本文模型与花岗岩单轴蠕变试验值[7]的对比

根据表2中不同温度下对试验数据进行反演得到的蠕变模型力学参数，从而得到模型参数与温度之间的函数关系：

(41)

由上式可知，弹性模量E2、黏滞系数η1、η2和η3随温度的增加逐渐减小；λ为常数；损伤参数α随温度的增加而增大，表明在温度越高同等外部应力作用下岩石的损伤程度越明显。

4.2 模型参数讨论

在实例验证的基础上，进一步对本文蠕变模型中的主要参数进行分析讨论，包括黏滞系数η1，η2，η3和损伤变量系数α。为了更好理解这些参数对岩石蠕变特性的影响，取表2中T=23℃时相应的参数，在分析某一参数的影响时，其他参数均保持不变。由式(6)可知，岩石发生损伤破坏时的临界应力σs与岩石的固有剪切强度τ0、内摩擦角φ和断裂韧度KIC有关，而且这3个参数有内在关联，因而为了便于分析模型参数敏感性，取临界应力σs为定值(即取σs=104 MPa)。

在其他参数保持不变的情况下，图9给出了不同黏滞系数η1，η2和η3分别对应的应变随时间变化关系曲线。由图9可知，随着黏滞系数η1的增加，蠕变变形稳定阶段的蠕变速率逐渐增大，而对初始段和加速破坏段蠕变速率的影响相对较小，η1< 500GPa·h时岩石蠕变变形的差异性很小。初始段的瞬态蠕变速率随着黏滞系数η2的增加而减小，且过渡到稳态阶段的时间逐渐增长，而η2的变化并不会对加速变形阶段的蠕变速率及最终的变形量产生影响。黏滞系数η3主要是影响加速变形阶段的蠕变速率和变形量，而不会对初始瞬态阶段和稳定阶段的变形情况产生作用。

图9 模型参数η1, η2, η3对蠕变特性的影响

不同α值下岩石的蠕变随时间关系曲线如图10所示。可知，参数α越大，岩石的损伤演化过程越明显，加速蠕变阶段的启动越早；当α较小时(如α<0.41h-1)，加速蠕变阶段并不显著。说明在蠕变模型中引入损伤变量可以很好地反映岩石蠕变变形的全过程。

图10 模型参数α对蠕变特性的影响

保持表2中T=23℃时的模型参数不变，不同临界损伤应力σs对应的蠕变变形曲线如图11所示，可以看出，当σ>σs时，恒定外力作用下岩石才具有明显的加速蠕变阶段，且σs越大加速蠕变过程越显著，加速蠕变阶段的启动越早；当σ<σs时，恒定外力作用下岩石不存在加速蠕变过程，只包括初始瞬态阶段和稳定变化阶段，且相同时间下岩石蠕变变形随σs增加而增大，很好地反映了岩石物理力学指标(如固有剪切强度τ0、内摩擦角φ和断裂韧度KIC)对其蠕变特性的影响，说明本文蠕变模型在稳态蠕变阶段考虑临界损伤应力和外部荷载的非线性影响是合理的。

由上述试验结果比对和模型参数分析可知，本文所建立的热-力耦合作用岩石损伤蠕变模型能够有效合理地反映岩石蠕变3个阶段各自的力学形态特征以及岩石蠕变变形的全过程特征。本文损伤蠕变模型只考虑了温度和外部应力作用的情况，并未涉及岩石内部存有裂隙以及裂隙内部存在高孔隙水压的情况。因而，本文蠕变模型适用于分析处于高温、高应力环境下较完整岩石的长期蠕变变形。

5 结论

(1)基于断裂力学推导了岩石临界损伤应力σs的表达式，提出了反映岩石稳态蠕变阶段与σs相关的非线性黏性分量，并将该黏性分量、σs和指数形式的损伤变量引入到岩石的流变本构关系和蠕变方程中，建立了新的热-力耦合的岩石损伤蠕变模型；当外部应力σ>σs时，岩石蠕变变形全过程具有明显三阶段特征：初始瞬态阶段、稳态阶段和加速蠕变阶段；当σ<σs时，岩石蠕变只包含前两个阶段。

(2)利用非线性最小二乘法得到了蠕变模型的拟合表达式，并通过反演获得了岩石损伤蠕变模型的力学参数，通过计算验证了本文蠕变模型和反演所得蠕变力学参数的可靠性。

(3)比对不同温度、不同应力条件下花岗岩三轴蠕变试验数据的结果表明，相比传统西原模型和Burgers模型，本文蠕变模型的计算结果能够更好地反映岩石在初始瞬态、稳态和加速蠕变阶段全过程的变形规律。

(4)模型参数的分析情况如下：黏滞系数η3、损伤变量系数α、和临界损伤应力σs的大小显著影响岩石的蠕变特性；η3和α越小、σs越大，岩石的损伤演化过程越明显，加速蠕变阶段启动越早；黏滞系数η2越小由初始瞬态过渡到稳态阶段的时间越早，但不会影响岩石后期的蠕变速率和最终的蠕变值；黏滞系数η1越大，对应的蠕变速率和蠕变变形越小，当η1增大到一定值时，其影响可以忽略。

本文建立的岩石损伤蠕变模型适用于分析高温、高应力环境下较完整岩石的长期蠕变变形，并未涉及岩石内部存有裂隙以及裂隙内部含有高孔隙水压的情况。