姜 鹏1，潘鹏志2,3，赵善坤4，吴振华2,3，陈 刚1
(1.武汉理工大学 理学院，湖北 武汉 430071; 2.中国科学院 武汉岩土力学研究所 岩土力学与工程国家重点实验室，湖北 武汉 430071; 3.中国科学院大学，北京 100049; 4.煤炭科学技术研究院有限公司 安全分院，北京 100013)
摘 要:在Perzyna黏弹塑性理论的基础上，引入基于应变能理论的岩石细观单元强度损伤模型，同时考虑岩石蠕变过程中蠕变速率随时间变化的特性，构建了能描述岩石从初始蠕变、稳定蠕变到非线性加速蠕变整个蠕变过程的细观黏弹塑性损伤耦合蠕变本构模型，并将该模型嵌入到三维弹塑性细胞自动机模拟系统(EPCA3D)中，通过实验数据验证该模型的正确性。应用该模型进行不同轴向应力、不同围压和不同均质度系数条件下的单、三轴压缩蠕变过程数值实验，结果表明:① 轴向应力提高增加了蠕变速率，因此减少了蠕变失效时间;② 围压和均质度系数的增加降低了蠕变速率并增加了蠕变失效时间，较好的再现了典型的实验现象。将该模型用于解释煤矿中“蠕变型”冲击地压的力学机理，较好的反映了煤矿巷道开挖后弹性应变能累积，围岩经历稳定蠕变和加速蠕变的损伤破坏过程，可为岩体工程的长期稳定性研究提供分析工具。
关键词:岩石力学;应变能;黏弹塑性损伤耦合;蠕变;冲击地压
中图分类号:TD313;TD324
文献标志码:A
文章编号:0253-9993(2018)11-2967-13
JIANG Peng1,PAN Pengzhi2,3,ZHAO Shankun4,WU Zhenhua2,3,CHEN Gang1
(1.School of Science,Wuhan University of Technology,Wuhan 430071,China; 2.State Key Laboratory of Geotechnical Engineering,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,Wuhan 430071,China; 3.University of Chinese Academy of Sciences，Beijing 100049，China; 4.Mine Safety Technology Branch,China Coal Research Institute,Beijing 100013,China)
Abstract:On the basis of elasto-viscoplastic theory of Perzyna,a mesoscopic coupled elasto-viscoplastic damage model,describing the whole creep process of rock from initial creep,stable creep to nonlinear acceleration creep,is constructed by introducing a mesoscopic unit strength damage model based on the strain energy theory.The characteristics of the creep rate change with time in the creep process of rock can be taken into consideration,and can be integrated into the EPCA3D.The model was verified by the experimental data and then applied to the numerical experiments of single and triaxial compression creep process under different axial stresses,confining the pressures and homogeneity coefficients.The results showed that,on one hand,the increased axial stress can speed up the creep rate,which reduces the creep failure time.On the other hand,the increased confining pressure and homogeneity coefficient can reduce the creep rate and lead to the increase of creep failure time.This study can better reproduce the typical experimental phenomena.It can be used to explain the mechanical mechanism of creep-induced rock burst in coal mines,and also better reflect the accumulation of elastic strain energy after the excavation of coal mine roadway and the damage process of surrounding rock undergoing stable creep and accelerating creep.This study can provide an analytical tool for the long-term stability of rock mass engineering.
Key words:rock mechanics;strain energy;coupling of elasto-viscoplastic damage;creep;rock bursts
移动阅读
姜鹏，潘鹏志，赵善坤,等.基于应变能的岩石黏弹塑性损伤耦合蠕变本构模型及应用[J].煤炭学报,2018,43(11):2967-2979.doi:10.13225/j.cnki.jccs.2018.8009
JIANG Peng,PAN Pengzhi,ZHAO Shankun,et al.A coupled elasto-viscoplastic damage model based on strain energy theory of rock and application[J].Journal of China Coal Society,2018,43(11):2967-2979.doi:10.13225/j.cnki.jccs.2018.8009
收稿日期:2018-08-07
修回日期:2018-11-15
责任编辑:郭晓炜
基金项目:国家重点研发计划资助项目(2017YFC0804203);中国科学院前沿科学重点研究资助项目(QYZDB-SSW-DQC029);中国科学院国际合作重点资助项目(115242KYSB20160024)
作者简介:姜 鹏(1993—)，男，湖北黄梅人，硕士研究生。E-mail:jiangpeng-lwx@whut.edu.cn
通讯作者:潘鹏志(1976—)，男，福建永春人，博士，研究员。E-mail:pzpan@whrsm.ac.cn
岩石的蠕变性质是岩石的重要力学特性之一，也是影响岩体工程长期稳定性的主要因素，岩石力学研究工作者一直致力于寻找合适的岩石蠕变本构模型来描述岩石的变形破坏与时间之间的演化规律，针对岩石的蠕变行为，国内外学者开展了大量的研究工作。
在岩石蠕变本构模型研究方面，ZHOU Hongwei等[1]在经典西原模型中用分数阶导数Abel dashpot代替牛顿缓冲器，提出了一种新的基于时间分数导数的蠕变本构模型，模拟结果比西原模型更好。XU Tao等[2]引入介观重整化的概念建立了单轴加载条件下非均质脆性岩石蠕变和应力松弛数值模型，模拟表明岩石复杂的蠕变行为可以通过元素和材料退化的小规模相互作用来解释。LI Xiang等[3]引入损伤力学的基本理论，提出了岩石统计损伤本构模型，可以反映应变软化现象。CHEN Liang等[4]基于花岗岩在不同温度下蠕变行为的实验研究，在西元模型的基础上，提出了一种基于损伤机理的蠕变模型。PING Cao等[5]基于岩石的非线性损伤蠕变特性和损伤变量，利用改进的Burgers模型，Hooke模型和St Venant模型，建立了高应力软岩新的非线性损伤蠕变本构模型，较为合理地解释了软岩蠕变变形。WANG Rui等[6]用分形缓冲器代替牛顿缓冲器，考虑损伤效应，提出了一种新的基于时间分形导数的本构模型来描述花岗岩的全蠕变过程。TANG Hao等[7]提出了一种新的基于变阶分数导数和连续损伤力学的四元素蠕变模型，新提出的蠕变模型模拟结果与常山岩盐的实验数据很吻合，认为分段处理是模拟岩石蠕变行为的有效方法。JIA Shanpo等[8]基于修正的Mohr-Coulomb准则，提出了一种新的非线性弹黏塑性损伤模型，用于研究高放废物处置库建造期间黏质岩中的蠕变和渗流。YANG Shengqi等[9]从应变能的角度建立了一个蠕变损伤模型，并与Burgers模型相结合，较好的描述了岩石蠕变的全过程。ZHAO Yanlin等[10]通过连接虎克体，虎克体和圣维南体的并联组合以及开尔文体和宾汉体的并联组合，提出了非线性弹塑性(EVP)蠕变本构模型，所提出的EVP蠕变模型可以精确描述加载和卸载蠕变行为。蔡美峰等[11]提出了岩石细观单元蠕变本构模型，合理的描述了岩石蠕变过程中力学性质的衰减，得到了广泛的应用。康永刚等[12]在经典的岩石蠕变模型基础上，提出一种改进的岩石蠕变模型，较好的描述了岩石蠕变全过程。沈才华等[13]从内能角度分析发生蠕变的机制，基于应变能理论，采用Perzyna黏塑性理论与西原正夫元件模型相结合，建立了一种能描述衰减蠕变、稳定蠕变和加速蠕变3个阶段全过程的蠕变统一本构模型。王其虎等[14]提出裂隙岩石塑性变形体元件，引入初始损伤影响因子，构建模拟岩石加速蠕变的蠕变损伤体元件。蒋昱州等[15]划分出了岩石由变形发展至破裂需经历的不同黏弹塑性状态，用Burgers模型来描述瞬时变形、黏弹变形及黏性变形，用Perzyna黏塑理论来描述黏塑性变形，并考虑后继屈服硬化作用，引入统计损伤因子，构建了能描述非线性加速蠕变特性的黏弹塑性损伤演化模型。
综合上述分析可以看出，现有的研究大多集中于将经典元件模型与损伤模型相结合来描述岩石的蠕变行为，取得了一定的进展。另外，对于蠕变加速阶段如何判别和模拟研究还存在争议，目前大部分得到的都是经验性公式，理论研究不够。SIH GEORGE C[16]采用应变能理论对裂纹扩展特性进行了研究，很好的解释了裂纹扩展引起的宏观应力应变过程，成为分析含裂纹体应变场变化的有效方法，而岩石蠕变过程与裂纹的扩展(即损伤)密不可分，蠕变可能引起损伤，而损伤又会加速蠕变的发展，两者相互影响、相互制约，在研究岩石蠕变行为时，需要考虑两者的耦合效应，而前人在该方面的工作较少涉及。基于应变能理论，考虑黏弹塑性和损伤模型耦合是需要进一步深入研究的问题。
另外，工程实测和模拟试验[17]表明冲击地压现象具有时间效应，即冲击地压在巷道开挖完成后，经过一段时间才发生。当应力水平较高时，煤岩体内局部区域出现裂纹集中现象，裂纹扩展搭接，发生等速蠕变。该区域煤岩进入应变软化变形阶段，抵抗变形的能力随着变形增加而减小，使得变形集中程度加剧，如此相互作用，煤岩变形迅速增加，可能发展为加速蠕变，最终导致煤岩破坏，如果破坏前积聚的应变能较高，有可能发生冲击地压。例如2014-04-15，耿村煤矿 13 采区西翼的 13180工作面下平巷发生了一次冲击事故，冲击位置距本工作面回采位置约 650 m，超出了回采扰动的影响范围，且其东翼的 13120 工作面已停采2个月，附近无其他外力扰动，属于蠕变型冲击地压。徐思朋[18]把这种冲击地压归类为煤岩体流变引起的冲击地压。从能量角度来看，冲击地压是岩体弹性应变能累积和释放的结果，可见针对基于应变能的蠕变模型研究对揭示蠕变型冲击地压发生机理有重要的作用。
因此，笔者在Perzyna黏弹塑性理论的基础上，从应变能反映岩石材料损伤程度的角度出发，引入基于应变能理论的岩石细观单元强度损伤模型，同时考虑岩石蠕变过程中蠕变速率随时间变化的特性，提出跟时间相关的黏塑性流动参数表达式，构建了岩石细观黏弹塑性损伤耦合蠕变本构模型，并将该模型嵌入到三维弹塑性细胞自动机模拟系统(EPCA3D)中，通过一系列数值实验验证该模型的正确性和合理性。并基于EPCA3D，模拟煤矿巷道开挖后的长期力学行为，解释蠕变型冲击地压发生的力学机理。
对于黏弹塑性问题，可以将总的应变ε分解为弹性应变εe和黏塑性应变εvp两部分，因此，总的应变率
可以表达为
(1)
则依赖于弹性应变率的总应力率
为
(2)
式中,
为弹性应变率;
为黏塑性应变率;D为弹性矩阵。
由下列形式的屈服准则判断物体内某一点是否进入黏塑性状态:
F(σ,εvp)-F0=0
(3)
其中，F0为单向屈服应力，是硬化参量κ的函数，与材料性质有关，在具体的计算中，F0与所选择的屈服准则有关。假定当且仅当F>F0的情况下，才产生黏塑性流动，这与弹塑性理论中当F>F0时产生塑性流动是一致的。
下面需要选择特定的规则来定义黏塑性应变，最简单的是认为黏塑性应变率依赖于当前的应力大小，因此存在:
(4)
该关系可以推广到包含应变硬化和温度的依从关系，也可以考虑与状态有关变量的影响。
式(4)的显式格式可以由下式的黏塑性流动法则给出:
(5)
其中，
(6)
其中，Q=Q(σ,εvp,κ)为塑性势函数;γ为黏塑性流动参数，用于控制黏塑性流动率。当x>0时，φ(x)为一个正的单调递增函数。
式(5)与传统非关联塑性流动法则类似。如果采用关联塑性流动法则，有F≡Q，式(5)变为
(7)
式中，a为流动矢量，且存在
(8)
对于函数φ(x)，可以有不同的形式，下面2种形式最常用:
(9)
(10)
式中，M,N为任意指定的常数;F为等效应力。
由式(7)我们可以通过隐式时间步长方案来定义在时间间隔Δtn=tn+1-tn内出现的黏塑性应变增量:
(11)
式中，Θ为时间积分因子，当Θ=0时为显式法，Θ=1时为全隐式法，
时为隐式梯形法或半隐式法，其中后2种方法具有无条件稳定的特点。
为了定义
可利用有限泰勒级数展开式展开为
(12)
其中，
(13)
式中，Δσn为在时间间隔Δtn=tn+1-tn内出现的应力增量。则式(11)变为
(14)
其中，
Cn=ΘΔtnHn
(15)
式中，矩阵H取决于应力水平，其具体形式见文献[19]。
式(2)的增量形式为
(16)
或者用位移增量来表示总应变增量:
Δεn=BnΔdn
(17)
将式(14)代入式(16)得
(18)
其中，
(19)
特别的，通过显式法解决线弹性问题时
此时式(18)变为
(20)
在任一瞬时tn都要满足下列平衡方程:
(21)
式中，fn为由于施加表面力、体力、温度载荷以及水压力等引起的等效节点力矢量;B为应变矩阵;Ω为定义域(积分区域)。这里认为孔隙水压力和温度应力以及外界应力共同作用产生的有效应力仅对瞬时变形产生影响，而对流变变形的影响可以忽略不计，但是水对岩石的软化作用、以及温度对岩石力学性能的影响对岩石的流变变形有重要的影响，例如，由于水对岩石的软化作用，岩石在偏应力作用下的蠕变速率增大。
在一个增量时间内，需要满足由式(21)的增量形式所给出的平衡方程，即
(22)
式中，Δfn为在时间增量Δtn内的载荷变化量，工程中碰到的大多数问题，其载荷增量是按不连续的时间间隔施加的，除了第一次时间间隔有增量变化外，对于其他时间步都有Δfn=0。
天然岩石具有不同尺度的宏观和微观缺陷，它们的变形和破坏不可避免地受到这些固有缺陷的影响,因此，考虑岩石细观单元强度蠕变损伤演化规律的岩石蠕变本构模型可以更合理地描述岩石的蠕变行为;下面的部分将介绍蠕变阶段的损伤过程，并建立岩石细观单元强度蠕变损伤模型[9]。
现有的研究表明[20]，岩石材料的破坏过程本质上是能量积累和耗散的过程，材料损伤不仅受随机分布的内部缺陷控制，而且受内部应力和应变状态的影响;岩石的声发射监测表明，在加速蠕变阶段开始后能量大幅增加，伴随着应变能的快速释放;鉴于上述讨论，应变能可以直接反映岩石的应力和应变状态;因此，用应变能作为内部变量来描述损伤演化规律更为合理。
应变能可以表示为
U=
σijdεij
(23)
式中，σij和εij分别为应力和应变，i,j=1,2,3。
因此提出了一种基于应变能的岩石细观单元强度损伤演化方程，其表达式为
σ=σ0(1-D)
(24)
其中，
(25)
其中，σ为单元强度;σ0为单元初始强度;D为损伤变量，当D=0时，在蠕变过程中不会发生损伤，D=1对应于蠕变失效的材料，因此，0<D<1对应于在蠕变过程中具有不同程度损伤的材料(图1);U为应变能，应变能可以用式(23)计算;U0为与初始蠕变损伤相对应的临界应变能;α为与材料蠕变有关的参数，并且总是大于0;U0和α的量值均取决于岩体的物理力学性质及岩体所处的温度、应力、渗流等环境因素，可由室内岩石实验或现场原位测试实验数据反演确定。
图1 岩石损伤演化规律示意
Fig.1 Rock damage evolution trend
与Perzyna黏弹塑性模型不同，引入了岩石细观单元强度损伤模型可以反映整个蠕变阶段，而Perzyna黏弹塑性模型只能反映初始蠕变阶段和稳定蠕变阶段，如图2所示。
图2 黏弹塑性损伤耦合模型和Perzyna黏弹塑性模型
Fig.2 A coupled elasto-viscoplastic damage model and Perzyna elasto-viscoplastic Model
此外，大量的现场监测和室内试验都表明，岩石材料弹性模量与岩石强度具有相似的变化规律，因此，假定每个单元的弹性模量也遵循类似式(24)的损伤演化方程。
由黏塑性流动系数的定义可知，黏塑性流动参数用于控制黏塑性流动率，也即岩石的蠕变速率，由蠕变试验现象可知，初始蠕变阶段对应着较大水平的蠕变速率，随着试验的进行，蠕变速率减少，进入稳定蠕变阶段后，其蠕变速率稳定在一个较低的水平。
在岩石的黏塑性流动系数方面，孙钧[21]进行了深入的研究，认为在非线性蠕变的情况下，黏塑性系数与岩土的黏滞性、物体所承受的应力和时间有关。潘鹏志等[22]在经典弹黏塑性理论的基础上，提出黏塑性流动系数张量表达式，建立岩石各向异性弹黏塑性蠕变模型。谭万鹏等[23]在引入复变量求导法的前提下，依据弹黏塑性有限元基本原理，采用牛顿-拉普拉斯变换迭代优化法来反算黏塑性流动系数等计算参数。关于岩石的黏塑性流动系数的确定方法较少见诸报道，在实际计算中，关于岩石黏塑性流动系数的确定比较困难，常黏塑性流动系数不能很好的描述蠕变速率的变化规律，鉴于此，基于对岩石蠕变机理的理解，本文提出与时间相关的黏塑性流动系数表达式，认为黏塑性流动系数与时间成指数关系，如图3所示。该表达式可以根据实验数据较快的反演分析确定。
图3 黏塑性流动系数随时间变化示意
Fig.3 Viscoplastic fluidity parameter changes with time
用指数型经验公式表示为
γt=γ∞+(γ0-γ∞)e-bt
(26)
其中，γt为t时刻的黏塑性流动系数;γ0为初始黏塑性流动系数;γ∞为t接近无穷大时的黏塑性流动系数;b为黏塑性流动系数衰减的系数，表征黏塑性流动系数衰减速度的快慢程度;t为时间;k和b的量值均取决于岩体的物理力学性质及岩体所处的温度、应力、渗流等环境因素，可由室内岩石实验或现场原位测试实验数据反演确定。
令
则式(26)变为
γt=γ0[k+(1-k)e-bt]
(27)
图4为不同k和b下黏塑性流动系数随时间变化的关系，由图4可知，通过选择适当的k和b值，可以较好地模拟实际岩石蠕变速率变化特性，从而模拟实际的蠕变变形曲线。
图4 不同k和b下黏塑性流动系数随时间变化的关系
Fig.4 Relationship between viscoplastic fluidity parameter and time under different values of k and b
为了简化计算，式(11)的计算采用显式法，时间积分因子Θ=0，则式(11)变为
(28)
将式(7)代入式(28)得
Δεvp=γ〈φ(F)〉aΔt
(29)
本文采用式(9)的流动方程，
且M=1，则式(29)变为
(30)
引入上述提出的损伤变量(式(25))与黏塑性流动参数表达式(式(27))，得到黏弹塑性损伤耦合模型增量形式的黏塑性应变为
(31)
为了表征岩石材料的非均匀特性，对于不含有肉眼可见裂隙的岩石，采用随机方法来进行表征，认为岩石的力学参数，包括变形参数、强度参数以及水力学参数服从一定的随机分布，通常采用Weibull随机分布[24]来表征岩石的非均匀性。Weibull统计概率密度函数定义为
(32)
式中，u为随机变量，如强度和杨氏模量;u0为比例参数，与平均单元参数有关;m为形状参数，反映了材料均匀性的程度，并被定义为均质度系数。
图5给出了在不同均质度系数(m=1.1，1.5，2，4，6)下，给定比例参数u0=100 MPa的单轴抗压强度统计密度分布曲线。根据Weibull随机分布和均质度系数的定义，由图5可知，较大的m意味着更多的单元将具有类似于给定的比例参数值的量值，从而获得更均匀的材料。图6给出了5个由5 000(50×100)个单元组成数值岩样，给定相同的单轴抗压强度参数，根据Weibull随机分布描述了不同均质度系数m下数值岩样的初始强度图像。在图6中，不同的颜色对应于不同的单元初始强度值，随着均质度系数m的升高，单元初始强度值分布在围绕平均值较小的范围内，从而获得更均匀的数值岩样，反之亦然。因此，均质度系数m是控制数值岩样宏观响应的重要参数。
图5 在不同均质度系数下，给定比例参数u0=100 MPa的单轴抗压强度统计密度分布曲线
Fig.5 Single compressive strength statistical density distribution curves with different homogeneity indices at a given scale parameter of u0=100 MPa
图6 不同均质度系数下的岩石试样强度分布
Fig.6 Rock specimen strength distribution with different homogeneity indices
在岩石破裂过程分析的三维弹塑性细胞自动机模拟系统(该系统的详细介绍见文献[25])的基础上，参考文献[19]中有关黏弹塑性问题的子程序，基于VC++开发工具，将上述三维细观黏弹塑性损伤耦合本构模型嵌入到三维弹塑性细胞自动机系统(EPCA3D)中，其计算流程图如图7所示。
数值模型采用二维平面应变模型，试样尺寸100 mm×50 mm，模型划分为50×25(1 250)个单元(图8)，本文中所有数值模拟的模型数值几何尺寸保持不变。模拟过程中在垂直方向对试样施加均布恒定载荷，在水平方向施加围压载荷(单轴蠕变施加围压为0 MPa)，然后固定保持轴压不变(图8)，试样在垂直方向采用位移约束。
图7 黏弹塑性损伤耦合蠕变本构模型计算流程
Fig.7 Flow chart of the coupled elasto-viscoplastic damage model
图8 数值模型
Fig.8 Numerical model
为了检验模型的适用性，对文献[26]中针对天然页岩样进行2种轴向应力状态下单轴蠕变试验的结果进行分析，所使用的输入参数在模拟中通过反演确定(表1)。值得注意的是，表1中列出的输入参数代表了宏观尺度上天然页岩样本的统计力学性质。图9，10分别为采用表1数据模拟得到的轴压55 MPa和65 MPa下蠕变试验与模拟结果。可以看出，Perzyna黏弹塑性模型的模拟曲线无法描述加速蠕变阶段;常黏塑性流动参数模型不能很好的描述稳定蠕变阶段的蠕变速率，稳定阶段的蠕变速率偏大，从而导致加速蠕变阶段时间提前;而本文所建模型模拟结果与试验结果在趋势和数值大小上均表现较好的一致，验证了该模型的正确性和合理性。
图11为模拟得到的岩石蠕变全过程曲线，模拟结果很好的再现了岩石蠕变破坏的3个典型阶段，即初始蠕变阶段、稳定蠕变阶段和加速蠕变阶段。
同一岩石试样在不同轴向应力下会表现出不同的蠕变特性，因此我们进行了一系列在不同轴向应力条件下的单轴压缩蠕变过程数值实验，施加的恒定应力大小分别为50，52，54，56，58和60 MPa，力学参数见表1，研究不同轴向应力对岩石蠕变过程的影响。图12为不同轴向应力下蠕变曲线，图13为不同轴向应力下蠕变速率随时间的演变曲线。图12和图13清楚地显示了在实验中看到的岩石蠕变的减速-加速现象，模拟结果显示轴向应力增加对岩石蠕变具有强烈的影响。
表1 数值模型参数
Table 1 Parameters of numerical model
图9 轴压55 MPa下蠕变试验与模拟结果
Fig.9 Creep test and simulation results under axial pressure 55 MPa
图10 轴压65 MPa下蠕变试验与模拟结果
Fig.10 Creep test and simulation results under axial pressure 65 MPa
图11 VEPCA3D模拟的岩石蠕变全过程曲线
Fig.11 Simulated rock creep curve by VEPCA3D
图12 同一试样不同轴向应力下蠕变曲线
Fig.12 Creep curves of the same specimens under different axial stress
从整体看来，当应力水平较高时，蠕变速率较高，轴向应力从50 MPa变化到60 MPa将最小蠕变速率提高了一个数量级(图13)。同时，由于较高轴向应力下的较大蠕变速率，蠕变失效破坏时间减小，例如轴向应力54 MPa下的蠕变失效时间为38.8 d，当轴向应力提高到60 MPa，其蠕变失效时间减小至2.4 d(图14)。
图13 同一试样不同轴向应力下蠕变速率随时间的演变曲线
Fig.13 Creep rate curves over time of the same specimens under different axial stress
图14 同一试样不同轴向应力下蠕变失效时间的拟合曲线
Fig.14 Creep failure time fitting curve of the same specimens under different axial stress
工程岩体开挖后，对岩体施加一定的支护，会有效控制围岩的时效变形，可以用不同围压条件下的三轴压缩蠕变过程数值实验来说明支护对蠕变的抑制效应，这里施加围压大小分别为0，1，2和3 MPa，并保持轴向应力恒定为65 MPa，力学参数依旧采用表1的参数，模拟结果显示围压对岩石蠕变具有强烈的影响。类似于之前的模拟，不同围压下蠕变曲线(图15)清楚地显示了在实验中看到的岩石蠕变的减速-加速现象，首先，当围压较高时，蠕变速率较低(图16):围压为0 MPa时的最小轴向和侧向蠕变速率均是围压为3 MPa时的5倍左右;随着围压增加，蠕变速率降低，蠕变失效时间增加。例如，单轴条件下的蠕变失效时间约为1.1 d，但当围压增加到3 MPa时，蠕变失效时间增加到约6.9 d。正如实验观察到的现象一样，围压的适当增加会显著改变蠕变速率的大小和蠕变失效时间。
众所周知，岩石是一种非均匀材料。为了研究均质度系数对岩石蠕变的影响，采用m=2,3,4和5的不同均质度系数进行单轴压缩蠕变过程数值实验，并保持轴向应力恒定为60 MPa，在模拟中，其它力学参数仍然采用表1的参数，只简单地改变均质度系数的值。如前所述，更大的均质度系数意味着更多的单元将具有类似于给定的力学参数值的量值(式(32)和图5)。因此，均值度系数较大的岩石试样将包含较少的低强度单元，从而变得更强且更脆。模拟的不同均质度系数下蠕变曲线(图17)以及蠕变速率随时间的演变曲线(图18)表明，均值度系数的增加导致蠕变速率的降低和相应的蠕变失效时间的增加。例如，均值度系数为2时的最小轴向和侧向蠕变速率大约是均值度系数为5时的4和4.5倍，并且随着均值度系数从2增加到5(图17和图18)，蠕变失效时间从约2.4 d增加到12.2 d。反之亦然，均值度系数的降低具有相反的效果，分别导致蠕变速率和蠕变失效时间的增加和减少。
图15 同一试样轴压65 MPa，不同围压下的蠕变曲线
Fig.15 Creep curves with different confining pressure of the same specimens under axial pressure 65 MPa
表2 巷道数值模型参数
Table 2 Parameters of Tunnel numerical model
图16 同一试样轴压65 MPa，不同围压下蠕变速率随时间的演变曲线
Fig.16 Creep rate curves over time with different confining pressure of the same specimens under axial pressure 65 MPa
图17 同一试样轴压60 MPa，不同均质度系数下蠕变曲线
Fig.17 Creep curves with different homogeneity indices of the same specimens under axial pressure 60 MPa
图18 同一试样轴压60 MPa，不同均质度系数下蠕变速率随时间的演变曲线
Fig.18 Creep rate curve over time with different homogeneity indices of the same specimens under axial pressure 60 MPa
近年来，随着煤矿开采深度的不断增加，冲击地压事故越来越频繁，深部一些新的冲击现象不断出现，蠕变型冲击便是其中之一。多起蠕变型冲击事故分析表明:该类冲击属于典型的隐蔽性灾害，具有自发性和时滞性的特点，发生机理的不明导致其防治极具挑战，给深部特厚煤层的安全开采造成了严重威胁[27]。蠕变型冲击与不稳定蠕变密切相关，开展煤岩巷道围岩的蠕变型冲击地压破坏过程模拟并解释其发生的机理，具有重要实践意义，因此本节以某煤岩巷道为背景，运用本文建立的模型，对巷道开挖后的围岩蠕变型冲击破坏过程进行模拟，研究巷道围岩力学行为随时间演化过程，力图明确巷道围岩变形状况随时间发展的变化规律，从而为巷道特别是煤岩巷道的长期稳定与安全提供合理的评价与建议。
数值模型的尺寸及其边界条件如图19所示，不考虑开挖进程对巷道变形的影响，即考虑为一次性开挖，并且假设巷道挖后不进行支护。模型采用非均质性模型和常黏塑性流动系数，高50 m，宽50 m;巷道采用直墙拱形巷道，宽7.5 m，高4.5 m。巷道布置在煤岩中，煤岩力学参数列于表2中。根据现场测定的地应力结果，施加大小为21.83 MPa的垂直应力和13.24 MPa的水平应力后维持不变。
图19 巷道数值模型及其边界条件
Fig.19 Tunnel numerical model and boundary conditions
图20为模拟得到的煤岩巷道开挖后左右两帮的水平收敛位移曲线。从图20可以看出，巷道开挖形成后，随着围岩应力状态的调整其变形经历了初始蠕变、稳定蠕变和加速蠕变3个阶段，在初始阶段巷帮围岩均有一个瞬时的位移增量，是煤岩体由三向受力变为二向受力后的回弹扩容，之后随着深部静载自重应力作用于巷道两侧，煤体内部受力变形，原生裂纹扩展的同时又促使次生裂纹生成，巷道变形进入初始蠕变阶段，蠕变速率较大。随着时间的延长，煤岩巷道围岩的内部应力重新调整完成，进入稳定蠕变阶段，蠕变速率基本保持稳定，围岩内部损伤变形依旧发展。当发展到一定时间时(31 d)，在巷道外部没有主控应力抑制围岩蠕变的前提下，随着内部累计弹性能的释放而进入加速蠕变阶段，巷道变形速度激增。为进一步说明巷道内部能量释放与蠕变过程的对应关系，对应分析了不同时刻围岩弹性应变能释放(即局部能量释放率LERR[28]:在岩体中开挖地下硐室后，围岩局部集聚的应变能超过岩体的极限储存能时，单位体积的岩体突然释放的能量)分布云图(图21)。从图21可以看出，0.1 d时巷道左右两帮有较少弹性应变能释放，此为开挖瞬时产生的围岩损伤破坏。之后，在稳定蠕变阶段，煤岩体内部弹性应变能释放有所增加。随着时间的推移，到31 d左右，应变能突然大量释放，巷道右帮形成较大的破坏区，具有冲击危险发生的可能。
因此，对于蠕变型冲击地压，在初期冲击地压危险评价时一定要重视静载荷的长时效应和危险评估，注重巷道围岩内部应力状态监测数据在一段时期内的对比分析(以往单纯评价单位时间内应力梯度变化对于蠕变型适用性有待进一步研究)，监控稳定蠕变期间的巷道应力变化，加强巷道围岩初期支护强度，提高锚杆(索)预紧力。同时，巷道卸压解危要以“降低静载强度、控制蠕变速度、转移应力深度”为原则，重视初始蠕变阶段和稳定蠕变阶段的应力控制。
图20 巷道左壁-右壁水平收敛位移曲线
Fig.20 Wall-to-wall horizontal convergence displacement curve of Tunnel
图21 不同时刻围岩弹性应变能释放分布云图
Fig.21 Elastic strain energy release cloud figure of surrounding rock in different time step
(1)本文在Perzyna黏弹塑性理论的基础上，从应变能反映岩石材料损伤程度的角度出发，引入基于应变能理论的岩石细观单元强度损伤模型，同时考虑岩石蠕变过程中蠕变速率随时间变化的特性，提出跟时间相关的黏塑性流动参数表达式，构建了能描述岩石从初始蠕变、稳定蠕变到非线性加速蠕变整个蠕变过程的细观黏弹塑性损伤耦合蠕变本构模型。
(2)通过实验数据验证了该模型的正确性，模拟结果很好的再现了岩石蠕变破坏的3个典型阶段，即初始蠕变阶段、稳定蠕变阶段和加速蠕变阶段。
(3)应用该模型进行不同轴向应力、不同围压和不同均质度系数条件下的单、三轴压缩蠕变过程数值实验，结果表明:① 轴向应力提高增加了蠕变速率，因此减少了蠕变失效时间;② 围压和均质度系数的增加降低了蠕变速率并增加了蠕变失效时间。模拟结果同室内试验所观察到的试验现象十分吻合。
(4)应用本文建立的基于应变能的蠕变模型，以煤岩巷道为工程背景，对煤岩巷道蠕变型冲击地压破坏过程进行数值模拟，模拟显示蠕变性冲击地压的发生与蠕变全过程相对应，从煤岩体弹性应变能的累积与释放角度，较好的解释了蠕变性冲击地压的发生机理。
参考文献:
[1] ZHOU Hongwei,WANG Chunping,HAN B B,et al.A creep constitutive model for salt rock based on fractional derivatives[J].International Journal of Rock Mechanics and Mining Sciences,2011,48(1):116-121.
[2] XU Tao,TANG Chunan,ZHAO Jian,et al.Modelling the time-dependent rheological behaviour of heterogeneous brittle rocks[J].Geophysical Journal International,2012,189(3):1781-1796.
[3] LI Xiang,CAO Wengui,SU Yonghua.A statistical damage constitutive model for softening behavior of rocks[J].Engineering Geology,2012,143:1-17.
[4] CHEN Liang,WANG Chunping,LIU Jianfeng,et al.A damage-mechanism-based creep model considering temperature effect ingranite[J].Mechanics Research Communications,2014,56:76-82.
[5] PING Cao,WEN Youdao,WANG Yixian,et al.Study on nonlinear damage creep constitutive model for high-stress soft rock[J].Environmental Earth Sciences,2016,75(10):900.
[6] WANG Rui,ZHUO Zhuang,ZHOU Hongwei,et al.A fractal derivative constitutive model for three stages in granite creep[J].Results in Physics,2017,7:2632-2638.
[7] TANG Hao,WANG Dongpo,HUANG Runqiu,et al.A new rock creep model based on variable-order fractional derivatives and continuum damage mechanics[J].Bulletin of Engineering Geology and the Environment,2018,77(1):375-383.
[8] JIA Shanpo,ZHANG Liwei,WU Bisheng,et al.A coupled hydromechanical creep damage model for clayey rock and its app-lication to nuclear waste repository[J].Tunnelling and Underground Space Technology,2018,74:230-246.
[9] YANG Shengqi,HU Bo,RANJITH Pathegama G,et al.Multi- step loading creep behavior of red sandstone after thermal treatments and a creep damage model[J].Energies,2018,11(1):212.
[10] ZHAO Yanlin,ZHANG Lianyang,WANG Weijun,et al.Separ-ation of elastoviscoplastic strains of rock and a nonlinear creep model[J].International Journal of Geomechanics,2017,18(1):04017129.
[11] 蔡美峰.岩石力学与工程[M].北京:科学出版社,2002.
[12] 康永刚,张秀娥.一种改进的岩石蠕变本构模型[J].岩土力学,2014,35(4):1049-1055.
KANG Yonggang,ZHANG Xiue.An improved constitutive model for rock creep[J].Rock and Soil Mechanics,2014,35(4):1049-1055.
[13] 沈才华,张兵,王文武.一种基于应变能理论的黏弹塑性蠕变本构模型[J].岩土力学,2014,35(12):3430-3436.
SHEN Caihua,ZHANG Bing,WANG Wenwu.A new viscoela-stoplastic creep constitutive model based on strain energy theory[J].Rock and Soil Mechanics,2014,35(12):3430-3436.
[14] 王其虎,叶义成,刘艳章,等.考虑初始损伤和蠕变损伤的岩石蠕变全过程本构模型[J].岩土力学,2016,37(S1):57-62.
WANG Qihu,YE Yicheng,LIU Yanzhang,et al.A creep constitutive model of rock considering initial damage and creep damage[J].Rock and Soil Mechanics,2016,37(S1):57-62.
[15] 蒋昱州,王奔,王瑞红,等.基于应变屈服临界的岩石黏弹塑性蠕变模型研究[J].长江科学院院报,2017,34(11):89-95.
JIANG Yuzhou,WANG Ben,WANG Ruihong,et al.Viscoelastic-plastic creep model of rock based on strain yield critical criteria[J].Journal of Yangtze River Scientific Research Institute,2017,34(11):89-95.
[16] SIH George C.From monoscale to multiscale modeling of fatigue crack growth:Stress and energy density factor[J].Science China Physics Mechanics and Astronomy,2014,57(1):39-50.
[17] 王洪英,张凤翔.大安山矿房山采区冲击地压模拟试验研究[J].辽宁工程技术大学学报:自然科学版,1999,18(5):491-493.
WANG Hongying,ZHANG Fengxiang.Simulating experiment investigation in Fangshan mining zone of Da Anshan Colliery[J].Journal of Liaoning Technical University:Natural Science,1999,18(5):491-493.
[18] 徐思朋,缪协兴.冲击地压的时间效应研究现状[J].矿业安全与环保,2001,28(2):27-29.
XU Sipeng,MIU Xiexing.Research status of time effect of rock burst[J].Mining Safety and Environmental Protection,2001,28(2):27-29.
[19] OWEN D R J,HINTON E.Finite elements in plasticity[M].Swansea:Pineridge press,1980.
[20] 陈惠发,A F萨里普,萨里普,等.弹性与塑性力学[M].北京:中国建筑工业出版社,2005.
[21] 孙钧.岩土材料流变及其工程应用[M].北京:中国建筑工业出版社,1999.
[22] 潘鹏志,冯夏庭,申林方,等.裂隙花岗岩各向异性蠕变特性研究[J].岩石力学与工程学报,2011,30(1):36-44.
PAN Pengzhi,FENG Xiating,SHEN Linfang,et al.Study of anisotropic creep behavior of fractured granite[J].Chinese Journal of Rock Mechanics and Engineering,2011,30(1):36-44.
[23] 谭万鹏,郑颖人.岩质边坡弹黏塑性计算参数位移反分析研究[J].岩石力学与工程学报,2010,29(S1):2988-2993.
TAN Wanpeng,ZHENG Yingren.Study of displacement back analysis of elasto-viscoplastic parameters of rock slope[J].Chinese Journal of Rock Mechanics and Engineering,2010,29(S1):2988-2993.
[24] WEIBULL W.A statistical distribution function of wide applicability[J].Journal of Applied Mechanics,1951,18(3):293-297.
[25] 潘鹏志.岩石破裂过程及其渗流-应力耦合特性研究的弹塑性细胞自动机模型[D].武汉:中国科学院武汉岩土力学研究所,2006.
PAN Pengzhi.Research on rock fracturing process and coupled hydro-mechanical effect using an elasto-plastic celluar automaton[D].Wuhan:Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,2006.
[26] 沈才华,王浩越,王媛,等.基于COD理论的岩石蠕变加速判别法探究[J].岩石力学与工程学报,2017,36(S2):3752-3759.
SHEN Caihua,WANG Haoyue,WANG Yuan,et al.Study on creep acceleration discrimination method based on COD theory[J].Chinese Journal of Rock Mechanics and Engineering,2017，36(S2):3752-3759.
[27] 姜福兴,冯宇,KOUAME,等.高地应力特厚煤层“蠕变型”冲击机理研究[J].岩土工程学报,2015,37(10):1762-1768.
JIANG Fuxing,FENG Yu,KOUAME,et al.Mechanism of creep-induced rock burst in extrathick coal seam under high ground stress[J].Chinese Journal of Geotechnical Engineering,2015,37(10):1762-1768.
[28] 苏国韶.高地应力下大型地下洞室群稳定性分析与智能优化研究[D].武汉:中国科学院武汉岩土力学研究所,2006.
SU Guoshao.Study on stability analysis and intelligent optimiza-tion for large underground caverns under high geostress condition[D].Wuhan:Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,2006.