0 引 言
岩石蠕变是指岩石矿物组构(骨架)随时间增长而不断调整重组，导致其应力、应变状态亦随时间而持续地增长变化[1-2]。随着矿山开采工程的扩大以及朝深部发展，矿山巷道在高应力状态下面临的开挖、支护、变形与维护问题愈加剧烈，在工程的设计、施工与运营维护过程中都需要研究围岩的蠕变力学性质[3-4]。
许多学者对传统的弹塑黏性元件进行串并联得到了许多线性元件组合模型，然而模型并不能描述加速蠕变阶段。为此，许多学者进行了非线性元件的研究[5]，并与传统元件进行串并联建立了非线性模型，用以描述岩石蠕变的全过程，为研究岩石的蠕变特性提供了新的方法与思路。然而非线性组合模型似乎只是为了使理论曲线与试验曲线相吻合而缺乏理论依据，且非线性元件的参数缺乏具体的物理意义[6-8]。岩石的蠕变破坏本质是岩石在长期荷载作用下，内部裂纹不断产生与发展，岩石材料力学性质逐渐劣化损伤的过程。在蠕变模型中考虑材料力学参数随应力状态以及时间的损伤变化，建立岩石的蠕变损伤模型更能反映蠕变的真实状态。张树光等[9]将强度参数引入到三维状态下的西原体模型中，并进一步通过损伤理论将模型整体蠕变参数进行非定常化，从而分析花岗岩在外力作用下蠕变损伤特性，最终应用反演法得出建立模型与实际状况基本相符，为后期围岩支护提供一定的理论依据。王明旭等[10]通过对预制裂缝花岗岩进行三轴压缩蠕变试验，分析在不同荷载作用、不同预制裂缝下岩石的蠕变规律，并用改进的开尔文模型对蠕变曲线进行拟合，得出在高压条件下花岗岩会出现松弛现象。
上述研究中蠕变模型并不能较好地描述加速蠕变阶段的特性，或者对于岩石加载过程中损伤演化特性未进行较深的研究，笔者拟通过研究花岗岩三轴蠕变特性，基于Burgers模型，引入改进的Kachanov损伤模型，建立描述加速蠕变阶段的蠕变损伤模型，并验证所建模型的正确性与合理性。
1 三轴蠕变试验分析
高德煤矿位于辽宁省阜新市，煤田所在区域大部分为第四系，地势呈波状起伏，地形、地貌简单，地质构造较复杂，断层构造较发育，煤层产状变化较小、厚度变化较大，连续性差，但是没有含煤地层出露，煤层深大约350～1 400 m，煤层构造结构单斜，倾斜方向为由北至南，倾斜角度为5°～15°，试样采自埋深为856.4 m的花岗岩层，质地均匀，细粒结构，无明显的层理，完整性与均匀性均较好。为确保试验的真实有效性，试样均采自同一块岩石，将试样加工制作成直径50 mm、高100 mm的圆柱体。测得围压15 MPa情况下平均抗压强度为292.6 MPa，泊松比为0.244，弹性模量为45.72 GPa。蠕变试验采用TAW-2000电液伺服岩石三轴试验仪，采用单体分级加载的方式，围压15 MPa，起始荷载125 MPa，分级荷载依次增加25 MPa，每级加载速率设定为0.5 MPa/s,各级荷载持续时间为24 h，采样记录间隔为1 min，得到的试样蠕变曲线如图1所示。
图1 花岗岩试样蠕变曲线
Fig.1 Creep curves of granite sample
由图1a可得，在各级荷载作用下，轴向、径向与体积均出现了瞬时变形，其后表现为不同的蠕变规律。存在长期强度值，当应力水平低于长期强度时，仅出现蠕变衰减阶段，蠕变速率很快衰减为0，当应力水平高于长期强度时，才呈现稳定蠕变与加速蠕变阶段，且轴向长期强度大于径向长期强度。在最后1级荷载作用下(持续时间6.7 h)，可以完整地观察到轴、径向衰减蠕变、稳定蠕变与加速蠕变阶段。通过观察图1b轴向蠕变曲线，当轴向应力水平低于175 MPa时(前3级荷载作用下)轴向蠕变曲线最终趋于稳定值，而高于200 MPa时(第4级荷载作用下)轴向蠕变曲线出现了稳定蠕变阶段。通过观察图1c径向蠕变曲线，前3级荷载作用下的径向蠕变曲线最终趋于稳定值，而第4级荷载作用下的径向蠕变曲线均出现稳定蠕变阶段。由图1d可得，试样在第1级荷载作用下的体积应变为正值，且逐渐增大，说明岩石在低应力水平条件下不会出现扩容现象，低应力仅能使岩石更加致密。试样轴向受压处于压缩状态，径向受拉处于拉伸状态，在低应力水平条件下轴向受压占据主要地位。在第2级荷载作用时，体积应变基本不随时间发生变化，体积应变后期蠕变速率基本为0。当试样处于第3级荷载作用下时，试样出现明显的扩容现象，体积应变逐渐减小。当试样处于第4级荷载作用下时，试样的体积应变由正值过渡到负值。因此在高应力状态下径向受拉占据主要地位。且随着应力水平的增加，可以发现体积应变的速率逐渐减小，蠕变曲线的斜率逐渐由正值过渡到负值。与轴径向蠕变对比分析，可以发现当应力水平处于径向长期强度(150～175 MPa)时，岩石材料开始扩容，当应力水平处于轴向长期强度(175～200 MPa)时，岩石材料体积应变开始由正值过渡到负值。各级荷载作用下的蠕变试验结果见表1。
表1 花岗岩三轴蠕变试验结果
Table 1 Experimental results of granite triaxial creep
轴向应力/MPa瞬时应变/10-2轴向应变ε1径向应变ε2体积应变εV最终应变/10-2轴向应变ε1径向应变ε2体积应变εV1250.357 8 -0.080 2 0.197 4 0.381 2 -0.085 00.211 11500.413 7 -0.107 3 0.199 2 0.427 3 -0.116 70.194 31750.446 0 -0.152 9 0.140 3 0.457 4 -0.166 90.123 72000.482 3 -0.209 1 0.064 2 0.506 0 -0.307 2-0.108 42250.578 0 -0.392 1 -0.206 3 0.609 6 -0.501 3-0.392 92500.658 4 -0.644 0 -0.629 6 0.720 8 -0.966 8-1.212 8
2 岩石蠕变损伤模型的建立
Burgers蠕变体是一种弹黏性体，由Maxwell体与Kelvin体串联而成，其力学模型如图2所示。
图2 Burgers体蠕变模型
Fig.2 Burgers creep mode
由Burgers模型可以推导出三维状态下的蠕变方程为
(1)
式中：εij为蠕变应变；σm为球应力，MPa；Sij为偏应力张量，MPa；K为体积模量，GPa；G1为Maxwell体剪切模量，GPa；G2为Kelvin体剪切模量，GPa；η1为Maxwell体黏滞性系数，GPa·h；η2为Kelvin体黏滞性系数，GPa·h；t为时间，h。
当采用Burgers模型时，能很好地描述岩石的瞬时弹性变形、衰减蠕变与稳定蠕变阶段，且模型简单实用，但是Burgers模型由于使用传统的线性元件，难以描述岩石的加速蠕变阶段。在衰减蠕变阶段，岩石呈现出明显的非线性硬化现象；在稳定蠕变阶段，岩石变形呈现线性变化现象；在加速蠕变阶段，岩石表现出非线性加速增大的现象，这表明岩石内部产生了新的损伤，且随时间的增长而加大直至岩石破坏，因此可以认为蠕变的损伤主要发生在加速蠕变阶段。Kachanov 在岩石损伤方面做了大量研究[11-12]，提出Kachanov损伤模型如下：
(2)
式中：λ和r为材料损伤参数；D为损伤变量；σ为应力，MPa。
当应力σ为常量时，得到蠕变损伤变量与时间以及蠕变寿命tF的关系式为
(3)
式中:tF为蠕变寿命，蠕变试验中直接获得。
为了仅反映岩石在加速蠕变阶段的损伤特性，将上式改为
(4)
式中：a为常数，可以表示为a=1/(r+1)；t1为加速蠕变起始时间，可根据蠕变速率曲线取得，蠕变速率呈现非线性增加的时刻即为t1。当t=t1时，岩石非线性损伤开始发生，当t=tF，损伤量为1，表明岩石结构已经完全损伤破坏。
一般情况下，岩石处于复杂的三向应力状态下，为了使得模型更好地描述这种应力状态，建立相关的三维状态蠕变模型很有必要[13-14]。因此，通过将式(1)进行进一步分解得到基于Burgers模型的蠕变损伤方程为
(5)
式中：σ1为轴向应力，MPa；σ3为围压，MPa。
当t<t1，模型退变为普通Burgers模型，用以描述衰减蠕变与稳定蠕变阶段，当t≥t1，模型为蠕变损伤模型，用以描述加速蠕变阶段。
3 模型参数的确定及其验证
根据实验曲线，当t→0时，εij=σm/3K+Sij/2G1，其为瞬时弹性应变。由于试验为室内假三轴试验，则球应力σm和偏应力张量Sij为
(6)
通过表1中的数据，结合式(6)和式(7)可以确定体积模量K与Maxwell体剪切模量G1为
(7)
然后将所得到的K与G1代入蠕变损伤方程，并通过观察轴、径向试验曲线是否进入加速蠕变阶段，分别选择Burgers模型与蠕变损伤方程进行拟合[15-17]。最后针对各级应力水平下的蠕变曲线，运用Origin软件基于最小二乘法原理，选择Levenberg Marquardt算法对模型参数进行辨识，得到围压为15 MPa作用下模型的蠕变参数见表2。通过拟合结果可以看出拟合参数决定系数均在0.95以上，说明拟合效果比较理想。
表2 Burgers模型的蠕变参数
Table 2 Creep parameters of Burgers model
轴向应力/MPa轴向蠕变参数K/GPaG1/GPaη1/GPaG2/GPaη2/GPat1/htF/haR212529.165 13.989 4 031.066 2 026.667 659.081 5.734 6.875 0.244 0 0.9815033.560 14.436 45 039.159 389.587 2 002.446 6.672 7.174 0.248 7 0.9817553.933 14.792 35 574.399 916.076 2 179.683 7.396 7.837 0.250 2 0.95200132.310 14.816 7 657.737 1 930.969 564.802 7.408 8.768 0.150 2 0.97225-45.633 11.987 9 998.836 597.948 1 529.682 5.993 9.126 0.140 6 0.99250-14.848 8.950 1 477.625 33 621.775 12.276 5.544 6.507 0.127 00.99
将辨识得到的蠕变参数代入蠕变方程得到损伤模型理论曲线，并与蠕变试验曲线对比，如图3所示。从对比曲线可以看出二者吻合效果相当理想。表明提出的损伤模型能够很好地反映花岗岩体的衰减蠕变、稳定蠕变与加速蠕变阶段的蠕变特性，验证了所建立蠕变损伤模型的合理性。
根据所建立的花岗岩蠕变损伤方程[18]，将最后1级荷载作用下蠕变参数，代入Kachanov损伤演化方程式(4)中。并通过数值计算绘制蠕变损伤变量与时间关系的曲线，如图4所示。
由图4可知，损伤变量D随着时间t的推移而增大，在初始衰减蠕变阶段产生较小的蠕变损伤，其曲线基本与时间轴重合，待到稳定蠕变阶段时，岩石内部裂隙进一步发育导致损伤程度加大，损伤速率呈现快速增大趋势，最后在加速蠕变阶段时损伤趋势随着内部缺陷发育和积累而变大，同时由于裂纹在岩石内部软弱区域逐渐聚集，最终在外荷载作用下裂纹贯通形成宏观裂缝，导致损伤变量急剧上升，呈现几乎竖直状态。因此，上述损伤变量D与时间t的变化规律较好地反映了花岗岩在荷载作用下的损伤特性，反映了Kachanov损伤模型适用于高德煤矿围岩的花岗岩蠕变损伤变形规律的描述。
图3 花岗岩蠕变试验曲线与模型拟合曲线对比
Fig.3 Comparison between fitting curves and creep test curves of granite
图4 损伤变量与时间关系
Fig.4 Relationship between damage variable and time
4 蠕变损伤模型的应用
将得到的花岗岩蠕变损伤参数代入到本构模型中，对实际围岩进行数值模型分析，其岩石峰值强度84.3 MPa，弹性模量19.83 GPa，泊松比0.165 2，黏聚力11.2 MPa，内摩擦角53.1°。
计算在34 m×30 m×2 000 m范围内巷道围岩，模型的地表取为自由边界，其余各侧面施加法向约束，其断面为半圆拱形，断面宽度为5 m，直墙高为2 m，拱高为2.5 m。
为便于分析说明蠕变损伤模型对高德煤矿围岩控制的意义，准确的预测出考虑损伤和不考虑损伤围岩的稳定性，选取巷道拱顶为监测点，检测其在不同条件下的蠕变情况如图5所示。因此，将所建立的花岗岩损伤本构模型代入到FLAC3D中进行数值计算，其中模型各个参数选取见表2。
计算得到的考虑和不考虑损伤的控制点的蠕变曲线如图6所示。由图6可知，在模拟二维巷道围岩蠕变过程中，建立的蠕变模型考虑损伤比不考虑损伤的位移要大许多，考虑损伤最大变形比不考虑损伤最大变形增大约0.055 27 mm，说明了蠕变损伤模型可以更好地反映巷道围岩的劣化效应，对于巷道的后期支护与安全具有重大指导意义。
图5 巷道控制点
Fig.5 Lane control point
图6 考虑损伤与否控制点蠕变曲线
Fig.6 Damage curves of monitoring points for consider damage or not
5 结 论
1)由径向蠕变确定的长期强度低于由轴向蠕变确定的长期强度，且在高应力状态下径向蠕变较轴向蠕变更加敏感，蠕变速率更大；在低应力水平条件下，花岗岩体积蠕变仅表现为压密；当应力水平达到径向长期强度时，开始出现扩容；当应力水平达到轴向长期强度时，体积应变开始由正值变为负值。
2)基于Burgers模型，引入改进的Kachanov损伤模型，建立了蠕变损伤模型，可以较好地描述岩石的衰减蠕变、稳定蠕变与加速蠕变阶段。且基于试验结果，采用Levenberg Marquardt算法对损伤模型参数进行辨识，识别效果理想。试验曲线与拟合曲线吻合度高，验证了所建立蠕变损伤模型的合理性。
3)在模拟二维巷道围岩蠕变过程中，建立的蠕变模型考虑损伤最大变形比不考虑损伤要大约0.055 27 mm，说明了蠕变损伤模型可以更好地反映巷道围岩的劣化效应，对于巷道的后期支护与安全具有重大指导意义。
[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].金属矿山,2009(11):119-123.
SONG Jinhu，HE Ruixia，BAI Zhe,et al.Deformation law simulation of deep tunnel in stress concentration zone[J].Metal Mine，2009,(11):119-123.
[3] 陈晓祥,勾攀峰,范增哲,等.深井高应力高突区域回采巷道变形特征及控制[J].采矿与安全工程学报,2013,30(3):363-368.
CHEN Xiaoxiang，GOU Panfeng，FAN Zengzhe,et al.Deformation characteristics and controlling of tailgate in high stress high gas and outburst region of deep mine[J].Journal of Mining & Safety Engineering，2013,30(3):363-368.
[4] 杨晓杰,庞杰文,张保童,等.回风石门软岩巷道变形破坏机理及其支护对策[J].煤炭学报,2014,39(6) :1000-1008.
YANG Xiaojie，PANG Jiewen，ZHANG Baotong,et al.Deformation and failure mechanism and support measures of the soft rock roadway in the air return laneway[J].Journal of China Coal Society,2014,39(6):1000-1008.
[5] 范庆忠,高延法,崔希海,等.软岩非线性蠕变模型研究[J].岩土工程学报,2007,3(4):505-509.
FAN Qingzhong,GAO Yanfa,CUI Xihai,et al.Study on nonlinear creep model of soft rock[J].Chinese Journal of Geotechnical Engineering,2007,3(4):505-509.
[6] 张耀平,曹 平,赵延林.软岩黏弹塑性蠕变特性及非线性蠕变模型[J].中国矿业大学学报,2009,38(1):34-40.
ZHANG Yaoping，CAO ping，ZHAO Yanlin.Visco-plastic rheological properties and a nonlinear creep model of soft rocks [J].Journal of China University of Mining & Technology,2009,38(1):34-40.
[7] 蒋昱州,张明鸣,李良权.岩石非线性黏弹塑性蠕变模型研究及其参数识别[J].岩石力学与工程学报,2008,27(4):832-839.
JIANG Yuzhou，ZHANG Mingming，LI Liangquan.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.
[8] 黄 明,刘新荣.三峡库区泥质粉砂岩非线性蠕变模型辨识[J].岩石力学与工程学报,2011,9(S2):3804-3810.
HUANG Ming，LIU Xinrong.Research on nonlinear creep model identification of siltite in Three Gorges reservoir area[J].Chinese Journal of Rock Mechanics and Engineering，2011,9(S2):3804-3810.
[9] 张树光,刘文博,王有涛.花岗岩非定常蠕变及其参数确定[J].防灾减灾工程学报,2017,37(3):435-441.
ZHANG Shuguang，LIU Wenbo，WANG Youtao.Study on non-stationary creep and parameter determination of granite[J].Journal of Disaster Preventtion and Mitigation Engineering，2017,37(3):435-441.
[10] 王明旭,许梦国,杜宇翔.三轴加载作用下矿化花岗岩的蠕变试验研究[J].化工矿物与加工,2017,46(10):49-54.
WANG Mingxu，XU Mengguo，DU Yuxiang.Study on creep test of grannite under triaxial loading[J].Journal of China Coal Society,2017,46(10):49-54.
[11] Kachanov M，Montagut E .Interaction of a crack with certain microcrack arrays[J].Engineering Fracture Mechanics，1986，25(5/6):625-636.
[12] Kachanov M.Introduction to continuum damage mechanics[M].The Netherlands:Martinus Nijhoff Publishers,1986.
[13] 朱元广,刘泉声,康永水,等.考虑温度效应的花岗岩蠕变损伤本构关系研究[J].岩石力学与工程学报,2011,30(9):1882-1888.
ZHU Yuanguang，LIU Quansheng，KANG Yongshui,et al.Study on creep damage constitutive relation of granite considering thermal effect[J].Chinese Journal of Rock Mechanics and Engineering,2011,30(9):1882-1888.
[14] 王来贵,赵 娜,何 峰,等.岩石蠕变损伤模型及其稳定性分析[J].煤炭学报,2009,34(1):64-68.
WANG Laigui，ZHAO Na，HE Feng,et al.Rock creep damage model and its stability analysis[J].Journal of China Coal Society,2009,34(1):64-68.
[15] 万志军，赵阳升，董付科，等.高温及三轴应力下花岗岩体力学特性的实验研究[J].岩石力学与工程学报,2008，17(1):72-77.
WAN Zhijun，ZHAO Yangsheng，DONG Fuke,et al.Experimental study on mechanical characteristics of granite under high temperatures and triaxial stresses[J].Chinese Journal of Rock Mechanics and Engineering,2008，17(1):72-77.
[16] 张 宁，赵阳升，万志军，等.高温作用下花岗岩三轴蠕变特征的实验研究[J].岩土工程学报200920(8):1309-1313.
ZHANG Ning，ZHAO Yangsheng，WAN Zhijun,et al.Experimental study on triaxial creep properties of granite under high temperature[J].Chinese Journal of Geotechnical Engineering，2009，20(8):1309-1313.
[17] 张志沛,王芝银,彭 惠.陕南泥岩三轴压缩蠕变试验及其数值模拟研究[J].水文地质工程地质,2011,38(1) :53-58.
ZHANG Zhipei,WANG Zhiyin,PENG Hui.Study on triaxial compression creep test and numerical simulation about mud rock in Southern Shaanxi[J].Hydrology & Geology,2011,38(1) :53-58.
[18] SHAO J F，CHAU KT，FENG X T.Modeling of anisotropic damage and creep deformation in brittle rocks[J].International Journal of Rock Mechanics and Mining Sciences，2006，43(4):582-592.
