硬脆性岩石多尺度损伤蠕变模型及长期强度研究

赵伦洋，赖远明，牛富俊，李鹏飞，朱其志

(1.华南理工大学亚热带建筑科学国家重点实验室，广东广州，510641；2.华南理工大学华南岩土工程研究院，广东广州，511442；3.中国科学院西北生态环境资源研究院冻土工程国家重点实验室，甘肃兰州，730000；4.巴黎东部大学多尺度模拟与仿真实验室，法国巴黎，77454；5.河海大学岩土力学与堤坝工程教育部重点实验室，江苏南京，210098)

随着我国对地下空间开发利用和对地下资源需求的日益增大，大量的深部重大地下工程已陆续启动实施，如地热开发[1]、高放废物地质处置库[2]、川藏铁路[3]等。大量研究表明[4−7]，深部高地应力硬脆性工程围岩在承受小于其抗压强度的荷载作用下，有可能产生较大的蠕变变形或诱发围岩迟滞性岩爆灾害，对工程的稳定性造成显著的不利影响。由于硬脆性岩石材料的复杂性(非均质性、多相性、工程时效性、体积膨胀现象等)及工程环境的不确定性，建立硬脆性岩石蠕变本构模型可以为深部岩体工程的长期稳定性分析提供重要的理论依据。现有的岩石蠕变本构模型主要有元件模型和经验模型两大类。元件模型采用弹簧、滑块和黏壶分别描述岩石的弹性、塑性和黏性行为，并通过它们之间的串并联来模拟岩石的应力、应变与时间的关系，其中，西原模型、Bingham模型和Burgers 模型得到了较为广泛的应用。这类模型虽然概念明确，表达直观，但由于其参数确定和模型识别的复杂性以及对于加速蠕变阶段描述的局限性，限制了其在工程上的应用。经验模型是基于不同试验条件下的试验结果和现场观察来确定应力、应变及时间函数关系的一种唯象模型，包括老化理论、时空理论、流动理论和绩效理论等。由于对试验的依赖性较大，经验模型无法适应环境因素的复杂变化。此外，元件模型和经验模型均是基于岩石受力变形的宏观表征建立，属于宏观流变学范畴。目前，从硬脆性岩石蠕变细观机理出发，建立能够描述硬脆性岩石蠕变过程中非线性特征的蠕变模型成为岩石流变力学研究的重点和难点。

长期强度作为岩石流变力学范畴中又一重要的科学问题一直备受关注。岩石长期强度的定义为其强度随时间而持续降低，并逐渐趋近于一个稳定收敛的低限定值[8]。目前确定岩石长期强度的方法主要有3 种：1)通过绘制应力−应变等时曲线簇，由曲线簇上的拐点确定长期强度[9]；2)根据稳态蠕变速率和应力水平的关系曲线，将曲线过渡到直线段的拐点确定为长期强度[10]；3)将体积应变开始反向时对应的裂纹损伤应力作为岩石的长期强度[11−12]。对于硬脆性岩石而言，裂纹损伤应力确定长期强度方法受到国内外学者的认同，MARTIN[11−12]通过对加拿大Lac du Bonnet花岗岩短期和长期强度进行系统研究，提出了以试样开始扩容时裂纹损伤应力作为该试样的长期强度；SZCZEPANIK 等[13]根据该理论得出花岗岩的单轴长期强度为常规单轴压缩试验峰值强度的70%～80%；SCHMIDTKE 等[14]的研究也验证了这一方法的可行性；李良权等[15]基于向家坝砂岩的常规三轴压缩试验和三轴压缩流变试验成果，指出岩石的裂纹损伤应力可反映长期强度所在应力水平。尽管国内外学者对于长期强度进行了大量的试验研究和理论分析，但是通过力学模型直接获得硬脆性岩石长期强度解析表达的研究甚少。

近年来，细观试验研究表明：微裂纹是硬脆性岩石非线性力学特性、蠕变变形和损伤的主要载体；微裂纹的扩展与当前的应力水平以及时间有着密切的关系[16−18]。在理论模型方面，SHAO等[19]指出硬脆性岩石内部微裂纹扩展方式包含：应力诱导的瞬时扩展和应力侵蚀导致的亚临界(subcritical)扩展，并指出后一种方式引起了材料的时效变形行为，进而建立了与裂纹扩展有关的时效损伤本构模型，从细观损伤角度为研究硬脆性岩石的时效特性提供了新的思路。在硬脆性岩石多尺度本构模型研究方面，朱其志等[20]基于Mori-Tanaka 均匀化方法和不可逆热动力学理论，构建了多尺度力学损伤摩擦耦合模型；ZHAO等[21−22]从岩石细观损伤累积和物理力学性质随时间劣化的角度，建立了硬脆性岩石统一多尺度损伤模型，但是该模型对加速蠕变阶段描述尚不充分，此外，模型参数的标定以及硬脆性岩石的长期强度研究均有待深入。本文作者在前期研究工作[19−23]的基础上，从细观角度出发，分析微裂纹扩展引起的宏观变形响应，建立可充分描述硬脆性岩石蠕变三阶段的多尺度损伤蠕变模型，并通过对摩擦损伤耦合和细观时效损伤演化规律的研究，构建模型参数的跨尺度标定方法，推导硬脆性岩石长期强度的解析表达式，为深部岩体工程问题的长期稳定性和耐久性分析提供帮助。

1 多尺度损伤蠕变模型

微裂纹成核和扩展引起的损伤是导致硬脆性岩石变形和破坏的主要力学机制。研究表明[24−25]这类损伤包括在应力变化下的微裂纹成核和扩展产生的瞬时损伤以及常应力下微裂纹发生亚临界扩展引起的时效损伤。据此，考虑到时间效应，对损伤d进行如下分解：

式中：di为瞬时损伤，dt为时效损伤。继而基于均匀化方法、热动力学理论和亚临界扩展理论建立硬脆性岩石多尺度损伤蠕变模型。

1.1 Helmholtz自由能与状态方程

硬脆性岩石宏细观结构[18]示意图如图1所示。从硬脆性岩石的细观结构出发(图1(b))，选取如图1(c)所示的代表性体积单元体(RVE)为研究对象，RVE的边界为∂Ω，整个RVE由线弹性固体基质和大量币形微裂纹组成。根据ZHU等[26]的研究成果，含微裂纹RVE的Helmholtz自由能Π为基质弹性自由能和与微裂纹相关的塑性功之和：

图1 硬脆性岩石宏细观结构示意图Fig.1 Schematic diagram of macroscopic and mesoscopic structure of hard brittle rock

式中：E为总应变；Ec为非弹性应变；Cm为各项同性基质弹性张量，Cm=2μmK+3kmJ，km和μm分别为基质体积压缩模量和剪切模量，通过引入二阶单位张量δ，四阶张量Kijkl=(δikδjl+δilδjk)/2-δijδkl/3，Jijkl=δijδkl/3；Cb为四阶塑性刚度张量，基于Mori-Tanaka均匀化方法可得其表达式[20]为

根据热动力学基本理论，应变E、非弹性应变Ec、损伤d的共轭力分别为

式中：Σ为宏观应力；Σc为作用在微裂纹表面的局部应力，可分解为偏应力Sc=K:Σc和球应力Pc=trΣc3；Dd为损伤驱动力。

1.2 摩擦准则和和瞬时损伤演化准则

在压应力作用下，硬脆性岩石材料的破坏以压剪类型为主，故可选取库仑摩擦准则来描述非弹性变形的演化。假设裂纹面内局部应力分布是均匀的，库仑摩擦滑动准则F可表示为[26]

式中：α为粗糙微裂纹面的摩擦因数。

出于简化考虑，选取相关联的流动法则，塑性势函数G=F，根据正交法则计算非弹性应变增量

式中：λc为塑性乘子；V=Sc/‖Sc‖。

在不可逆热动力学理论框架下，对于损伤行为，基于应变能释放率的瞬时损伤演化准则可表述为[26]

式中：R(d)为损伤发展抗力，根据损伤抗力函数的基本特征[22]，其表达式如下：

式中：ξ=d/dc；dc为损伤临界值，对应于材料的峰值应力处的损伤，当且仅当d=dc时，R(d)=R(dc)。

运用正交化准则计算瞬时损伤增量：

式中：λd为损伤乘子。结合式(4)和式(7)可以看出，损伤准则只与摩擦产生的非弹性应变Ec和当前损伤d有关，即损伤由微裂纹面上的摩擦滑移驱动。反之，损伤演化通过Cb对控制摩擦滑移的局部应力Σc产生影响。因此，摩擦和损伤的演化是强耦合的，可以通过联立准则(5)和(7)的一致性条件方程求解塑性乘子λc和损伤乘子λd。

1.3 时效损伤演化准则

时效损伤变量dt是关于时间t和损伤d的函数，即

当t→∞时，dt→。为t时刻岩石细观结构达到平衡状态时的时效损伤。当dt<时，系统尚未达到平衡状态，材料细观结构将不断向自平衡状态发展。文献[16]指出：细观结构演化在动力学方面可以解释为系统偏离平衡状态的距离，即(-dt)，因此，时效损伤演化准则可用线性形式来描述[16−17]：

式中：γ为与材料性质有关的常数；∊[0,∞)；

对于硬脆性岩石，当偏应力水平低于岩石长期强度时，岩石仅会发生衰减蠕变和稳态蠕变，并且稳态蠕变的速率很小，变形最终会趋于稳定；当偏应力水平高于岩石长期强度时，将会发生加速蠕变，并导致岩石迅速破坏。在微裂纹各向同性分布假定下，非弹性应变Ec可分解为

式中：标量β为非弹性体积应变，β=trEc；二阶张量Γ=K:Ec为非弹性偏应变。

在应力空间，式(5)可表示为

式中：S=K:Σ；P= trΣ3。

由式(13)可知：当岩石受到的应力水平小于岩石初始屈服应力时(F<0)，损伤未发生；当应力水平高于初始屈服应力时(F=0)，在蠕变试验过程中应力保持恒定，则和均为常数，即宏观非弹性应变Ec与细观损伤d的比值保持恒定。忽略弹性后效作用，这一理论结果与MOGI[27]和OHNAKA[28]通过试验观测到的结果相一致。因此，为反映硬脆性岩石在各蠕变阶段的变形规律，对加速蠕变阶段进行更深入的研究，在文献[16]的基础上引入细观参数的表达式：

式中：B和n为模型参数，B>0，控制瞬态蠕变和稳态蠕变阶段时效损伤的演化速度，n>0，控制加速蠕变阶段时效损伤的演化速度。显然，当损伤达到损伤临界值dc时，加速蠕变启动。

1.4 增量蠕变本构关系

在建立蠕变本构模型过程中，假定弹性应变全部是由应力变化产生，非弹性应变包括与时间无关的非弹性变形和与时间相关的非弹性变形，其中，与时间相关的非弹性变形由时效损伤所引起。

考虑时间效应后，当应力保持不变时，库仑型摩擦准则式(5)依然适用。与时间相关的非弹性应变增量也可根据流动法则求得

式中：λct为与时间相关的塑性乘子。

根据摩擦准则的一致性条件：

求解可得：

综上，蠕变本构关系的增量形式可写为

式中：Sm=(Cm)-1为基质柔度张量。

2 长期强度研究

MARTIN等[12,29]研究发现：硬脆性岩石材料在压缩荷载作用下的应力−应变关系可分为裂纹闭合阶段、线弹性阶段、稳定裂纹扩展阶段和不稳定裂纹扩展阶段[12]，分别对应岩石裂纹闭合应力Σcc、裂纹初始应力Σci、裂纹损伤应力Σcd和峰值强度Σf。裂纹损伤应力Σcd对应于岩石加载过程中体积应变−轴向应变曲线上体积变形的拐点。MARTIN等[12,14]研究表明：裂纹损伤应力Σcd可作为硬脆性岩石的长期强度。

由弹塑性理论可知，体积应变Ev包含弹性体积应变Ev,e和非弹性体积应变Ev,p

对于常规三轴压缩试验，忽略围压所产生的体积应变和裂纹闭合阶段产生的非弹性应变，在轴向加载过程中

式中：P0为围压；Λc为偏压引起的累计塑性乘子，表达式为[20,30]

ZHU[30]通过摩擦损伤耦合分析推导了在常规三轴加载情况下，摩擦准则(式(5))在主应力空间的表达形式为

将式(21)和式(22)代入式(20)并化简可得

对损伤变量d进行求导，将=0(体积压缩/膨胀转化点)化简，并结合边界条件得

通过求解式(24)，得到满足条件的解ξ∞(见图2)，从而求得长期强度Σ∞对应的瞬时损伤为

图2 损伤抗力和轴向应力演化示意图Fig.2 Diagram of damage resistance and axial stress evolution

将式(25)代入式(22)，即可得到不同围压下硬脆性岩石长期强度的解析表达式为

3 模型验证

为验证所提出的多尺度损伤蠕变模型的准确性和有效性，本文采用该模型模拟向家坝砂岩在不同围压下的三轴压缩蠕变试验。李良权等[15]在围压分别为3，5 和7 MPa 下对向家坝砂岩进行了常规三轴压缩试验和三轴压缩蠕变试验，并对向家坝砂岩的长期强度进行了分析，其常规三轴压缩试验和三轴压缩蠕变试验的围压和偏压加载速率均为7.5 MPa/min，其中，蠕变试验采用多级加载方法，详细的试验方案和试验曲线见文献[15,31]。

3.1 模型参数跨尺度标定

本模型涉及的5个力学参数和3个时效参数均可通过室内试验和理论推导确定。

1)基质的弹性模量Em和泊松比νm可通过常规三轴压缩试验应力−应变曲线线性段确定。

2)由式(22)得到岩石在受压状态下的强度为

显然，式(27)可以与摩尔−库仑强度准则直接关联，继而建立细观参数(α和R(dc))与材料宏观强度参数(内摩擦角ϕ和黏聚力c)的跨尺度关系

3)损伤临界值dc对应于应力峰值处的损伤，可通过声发射试验确定，即在常规三轴压缩试验时记录不同围压下峰值应力处声发射信号的累积计数。在本构模型中，dc在应力−应变关系上的宏观表现为控制峰前和峰后非弹性变形，不影响岩石的峰值强度。LOCKNER[32]指出，dc与围压大致呈线性关系。针对向家坝砂岩，由于其围压差别不大，本文将dc取为常数。运用上述方法最终确定向家坝砂岩的5个力学参数如表1所示。

表1 基本力学参数Table 1 Basic mechanical parameters

4)当岩样受到长期强度荷载时，瞬时损伤和时效损伤满足如下关系

结合式(14)和式(28)，可以得到模型参数

将表1中的力学参数代入式(24)和式(25)，计算得到长期强度对应的瞬时损伤d∞=0.62。将dc和d∞代入式(30)，得到B=0.9。

5)蠕变参数γ可通过三轴压缩蠕变试验中瞬态蠕变的速率来确定[16]，其取值与材料的性质有关，一般，软岩的γ取值在10-5～10-6量级，硬岩的γ取值在10-4量级，针对向家坝砂岩取γ=10-4。蠕变参数n可通过三轴压缩蠕变试验加速蠕变持续时长来确定，为简便取n=1。

3.2 单级加载下长期强度预测

通过式(26)可以得到在单级加载情况下围压为3，5和7 MPa砂岩的长期强度。表2所示为多尺度模型确定的长期强度与李良权等[15]运用应力−应变等时曲线簇法和稳态蠕变速率−应力水平的关系曲线法确定的长期强度对比。从表2可以看出：除围压为3 MPa下多尺度模型对长期强度的预测结果较常用方法确定的长期强度偏小，围压为5 MPa 和7 MPa下多尺度模型预测的长期强度与上述2种方法得到的结果相一致。

表2 向家坝砂岩在不同围压下的长期强度Table 2 Long-term strengths of Xiangjiaba sandstone for different confining pressure MPa

图3所示为不同围压、长期强度荷载下砂岩的蠕变演化模拟曲线。从图3可以看到典型蠕变的三阶段特征，即瞬态蠕变、稳态蠕变和加速蠕变，加速蠕变发生时间在200～250 h之间。

图3 不同围压下长期强度荷载下砂岩的蠕变演化模拟曲线Fig.3 Simulation curves of creep evolution of sandstone under long-term strength with different confining pressures

图4所示为围压5 MPa时，长期荷载作用下砂岩的轴向蠕变、轴向蠕变速率以及细观损伤演化的模拟曲线。从图4可见：细观损伤演化曲线与蠕变曲线具有明显的一一对应关系，尤其是根据临界损伤可以很直观地在三阶段蠕变曲线中找到加速蠕变启动点。此外，还可以发现蠕变应变的3个阶段与应变速率紧密联系。

图4 围压5 MPa时长期强度荷载下砂岩的轴向蠕变、轴向蠕变速率以及细观损伤演化的模拟曲线Fig.4 Simulation curves of creep,creep rate and microscopic damage of sandstone under long-term strength with confining pressure of 5 MPa

图5所示为围压5 MPa时不同轴向荷载作用下砂岩的蠕变曲线。从图5可以看到，当轴压为155 MPa时，不会发生加速蠕变，即不会产生蠕变失效；当轴压为157.1 MPa时，发生蠕变失效，且失效时间约为240 h；当轴压为160 MPa 时，蠕变破坏时间缩短约40 h。说明当荷载低于长期强度时，岩石不会发生蠕变失效，而当荷载高于长期强度时，随着轴压增大，蠕变失效时间在逐渐减小，这与众多试验观察结果相一致。

图5 围压5 MPa时不同轴向荷载下砂岩的蠕变模拟曲线Fig.5 Simulation curves of creep evolution of sandstone under different axial load with confining pressure of 5 MPa

当围压为5 MPa，长期强度荷载下n=0.5，1.0，1.5 和2.0 时多尺度蠕变模型获得的砂岩轴向蠕变和细观损伤演化曲线如图6所示。从图6可以看出，参数n仅对加速蠕变段产生影响，随着n增大，加速蠕变阶段减少。换言之，参数n对蠕变失效时间有很大的影响。

图6 不同参数n下砂岩的轴向蠕变和细观损伤模拟曲线Fig.6 Simulation curves of creep and microscopic damage of sandstone with different parameter n

3.3 三轴压缩蠕变试验模拟

为进一步对多尺度损伤蠕变模型进行验证，采用3.1 节的模型参数对围压为5 MPa 和7 MPa 时向家坝砂岩三轴压缩蠕变试验应变−时间全过程曲线进行数值模拟(图中蓝色曲线)，并与试验数据进行对比，结果如图7所示。从图7可以看出，模型可以较为完整地描述不同围压下硬脆性岩石分级加载蠕变的力学行为，包括在低应力水平下衰减蠕变和稳态蠕变的特征、在高应力水平下加速蠕变破坏的特征。值得指出的是，用分级加载的方式，上一级加载的应力水平会对岩样造成不同程度的损伤，而随着加载级数的增多，损伤逐渐增加。在图7(a)的模拟中，轴向应力为155 MPa等级下仅经历约50 h后加速蠕变发生。而从图5可以看出，在单级加载情况下，轴向应力为155 MPa等级下不会发生加速蠕变破坏。说明从岩石材料随时间损伤累积和物理力学性质劣化的角度出发建立的多尺度蠕变损伤模型很好地考虑了多级加载下损伤的累积效应。此外，根据本文预测的长期强度可以对分级加载蠕变试验的荷载分级提供一个可靠的参考依据。例如，在围压为5 MPa的分级加载蠕变试验中，将轴向应力为155 MPa的荷载等级持续时间延长至约150 h，若岩样未破坏，再加荷载至165 MPa；在围压为7 MPa的分级加载蠕变试验中，在荷载等级为177 MPa与187 MPa之间建议增加一个持续时长为75 h的182 MPa荷载等级。此外，结合图7的分析结果，通过调整参数n可以提高模型预测结果的准确性，如图7中红色曲线所示。

图7 向家坝砂岩三轴压缩蠕变试验结果[15]与数值模拟对比Fig.7 Comparison between numerical simulation and experimental data[15]of triaxial compression creep test of Xiangjiaba Sandstone

4 结论

1)基于均匀化方法建立多尺度损伤蠕变模型是可行的，和通常的宏观损伤力学蠕变模型相比，本模型所体现的力学机制更为明确，而且所含参数少，物理意义明确，并均可通过室内试验和理论推导予以确定。

2)通过摩擦损伤耦合分析和细观时效损伤演化规律研究，推导了硬脆性岩石长期强度的解析表达式。通过对向家坝砂岩长期强度预测和对比，验证了运用多尺度模型确定长期强度的可行性，为多级蠕变试验的荷载分级提供理论依据。

3)模型可以较为完整地描述硬脆性岩石在不同围压下分级加载和单级加载蠕变的力学行为，包括在低应力水平下衰减蠕变和稳态蠕变的特征、高应力水平下发生加速蠕变破坏的特征。通过多级加载和单级加载蠕变试验的模拟对比表明模型可以较为准确地考虑到应变历史对损伤演化的影响，说明基于损伤累积和物理性能随时间劣化的角度出发建立蠕变模型的合理性。

4)需要指出的是本文所提出的是一种基于细观力学的各向同性简化模型，而在实际工程中，岩体的应力状态往往较为复杂，各向异性损伤以及与微裂纹张开闭合有关的单边效应都需要予以考虑。