蠕变不仅对煤与瓦斯突出、冲击地压、顶板塌陷等矿井煤岩动力灾害具有十分重要的作用[1]，而且对巷道围岩稳定性也具有显著影响[2]。因为蠕变与工程实践联系紧密，所以蠕变成为流变领域中最重要的研究内容。
国内外学者对蠕变做出了大量研究，20世纪50年代Sally(1956)和Rotherham(1951)研究岩石蠕变，发现典型的蠕变3阶段性态，即瞬间蠕变、稳态蠕变、加速蠕变[3]。曹树刚和鲜学福[4]通过单轴和三轴加载在充气和不充气情况下的实验分析了煤岩蠕变损伤特性，认为蠕变都会出现衰减蠕变、稳定蠕变、加速蠕变3阶段。MARANINI E和BRIGNOLI M[5]通过软岩单轴、三轴分级加载实验也出现相同规律。为了得到蠕变的函数表达式，CRUDEN[6]较早就提出了幂函数模型，MA L等[7]和王登科等[8]曾在实验研究中采用幂函数拟合蠕变曲线。为了研究岩石的力学模型，JCMAXWELL提出了弹簧和黏壶串联的麦克斯韦模型;沃格特(VOIGT)提出过弹簧和黏壶并联的沃格特模型;由于这些模型过于简单不能体现蠕变的特征，学者们先后又提出了MAXWELL-WEICHERT模型、VOIGET-KELVIN模型[9]、西原模型和BURGERS模型等。阎岩等[10]以西原模型为基础通过最小二乘法求出实验中的各个蠕变参数。邓荣贵等[11]基于非线性黏滞元件提出了新的流变模型。徐鹏和杨圣奇[12]通过循环加载实验提出了黏弹塑性蠕变模型。但是，各种模型的普适性尚存在一定的缺陷，却没有相关文献对这些蠕变模型进行归纳总结，以得出一种可靠性较好的蠕变模型。而且，学者们研究岩石蠕变时仅仅分析了蠕变速率、蠕变加速率，基本没有更进一步的研究。
因此，本文在前人研究基础上，进行低瓦斯压力下分级加载轴压时煤岩蠕变试验，一方面观察低瓦斯压力时煤岩蠕变规律;另一方面根据实验数据检验各种蠕变模型的适用情况及优缺点，并对相应模型进行改进，以期得到最实用的煤岩蠕变模型。本文在蠕变理论上提出蠕变拐点的观点，结合蠕变模型可为岩石失稳破坏及煤岩动力灾害等工程实践提供参考。
1 含瓦斯煤蠕变试验
1.1 煤样制备
本实验研究煤样取自潞安集团余吾煤业，在取芯机钻取煤样时由于煤质较软导致成品率太低。考虑到原煤内部孔隙结构不均，各向异性，为了减小外界因素引起的系统误差，排除各向异性对实验造成干扰，本文拟采用型煤代替原煤实验，而相关文献发现型煤代替原煤进行力学实验是可行的[13]。制作型煤时为了更加符合原煤实际情况，将煤样加工成煤粉，再将煤粉混合均匀，加入一定的水和煤焦油制成直径为50 mm、高为100 mm的圆柱形煤样，然后以100 ℃高温干燥12 h排除水分干扰。
1.2 实验系统
为了使得实验稳定、数据准确，课题组完善了三轴应力状态下的煤岩应变测试系统。该系统主要包括:大型压力机(天辰实验机)、三轴夹持器、围压加载装置、ES-I型温控仪、气源(浓度为99.99%的甲烷)、JC-4A静态应变仪等组件。具体实物如图1所示。
图1 实验系统实物
Fig.1 Object picture of test system
1.3 实验方案
岩石蠕变试验的加载方法主要有3种:单级加载法、逐级加载法、逐级增量加卸载法[14]。鉴于分级加载法的可操作性强、采用率高，本实验采用分级加载轴压的方法，研究不同瓦斯压力条件下型煤蠕变规律。由于该煤层瓦斯压力较低(0.1～0.4 MPa)，为了真实模拟现场实际，本文设置了4组对照实验(瓦斯压力分别为0.13,0.20,0.30,0.40 MPa)。实验过程中保持温度恒为25 ℃、围压恒为2 MPa。具体实验方案为:先施加较小轴压，然后以0.1 MPa/s的速率加载围压至2 MPa，再打开气瓶进气阀从三轴夹持器上端进气口通入甲烷气体，使进气口瓦斯压力维持在实验预设值，并保持一定时间待煤体吸附平衡形成含瓦斯煤;然后以0.1 kN/s的速率加载轴压至预设值，每级轴压维持足够时间待轴向应变为定值后，再加载下一级轴压，加载到最后一级使得煤样破坏后停止;更换煤样重复上述步骤再次实验，实验过程中记录轴向应变数据。
2 蠕变曲线及实验结果
含瓦斯煤蠕变通常分为衰减蠕变和非衰减蠕变两大类，具体如图2所示。
图2 煤蠕变曲线
Fig.2 Coal creep curves
当轴向应力较小未达到煤体屈服应力时，煤体轴向应变(ε)随时间延长先增大后趋于稳定，蠕变曲线呈上凸状。
当轴向应力较大时，煤体轴向应变随时间延长过程中，一般表现出典型的3阶段:
I阶段:减速蠕变阶段[15-16]，应变增大，蠕变速率逐渐减小，即
蠕变加速度标量
此阶段蠕变曲线呈上凸状。
II阶段:稳定蠕变阶段[17]，随时间延长，应变缓慢增大，蠕变速率
减小。
III阶段:加速蠕变阶段，煤岩应变大于屈服应变时，蠕变速率迅速增大，应变急剧增大。蠕变曲线凸凹性改变，呈下凹状。
经上述分析发现，发生非衰减蠕变时蠕变曲线的凸凹性会发生变化，曲线上至少存在1个拐点，暂且称之为蠕变拐点。未达到蠕变拐点时煤岩只会发生衰减蠕变，超过蠕变拐点时煤岩就会发生加速蠕变。假设蠕变曲线是连续函数，那么，拐点处二阶导可能为0或者不存在，说明蠕变拐点处蠕变加速度不一定是0，但煤岩发生非衰减蠕变时函数曲线凸凹性会改变，蠕变曲线上一定存在蠕变拐点，只不过蠕变拐点处的蠕变速率不一定为0。由此可知，蠕变曲线上至少存在1个蠕变拐点是发生非衰减蠕变的充要条件。
此外，当轴向应力一开始就过大时，煤岩骨架变形急剧增大，直接造成煤岩破坏，短时间内就出现了加速蠕变，轴向应变急剧增大。
本文研究了A1～A4四个型煤煤样在瓦斯压力分别为0.13,0.20,0.30,0.40 MPa时，分级加载轴压直至煤样发生破坏过程中，轴向应变随时间的变化情况。利用MATLAB将实验数据绘制成图3。
图3 分级加载轴压时的蠕变试验
Fig.3 Creep test under axial loading at grade loading
由图3可知，不同瓦斯压力条件下，含瓦斯煤蠕变情况存在着异同点:
(1)轴压和围压相同时，瓦斯压力越大，轴向应变越大。例如，轴压为6 MPa时，A2煤样轴向应变是A1煤样轴向应变的1.4～1.7倍。
(2)同一瓦斯压力下，分级加载轴压过程中煤样轴向应变呈梯度增大，分析原因是上1级轴压加载结束时，煤样发生残余变形。A1和A2煤样分级加载轴压经历36 h后受轴压15 MPa作用时，直接表现出加速蠕变;A3煤样经过60 h后受18 MPa轴压作用时表现出加速蠕变;A4煤样经过48 h后在轴压为15 MPa时表现出典型的非衰减蠕变特征。
(3)A1,A2,A3煤样受到最后1级轴压作用时，立即发生加速蠕变，轴向应变急剧增大，蠕变曲线呈下凹状;A4煤样受到最后1级轴压作用时，轴向应变在前0.2 h内缓慢增大，0.2 h后急剧增大，蠕变曲线上存在蠕变拐点，呈典型的非衰减蠕变特征。由此可见，轴压、围压和瓦斯压力一定时，煤体发生非衰减蠕变时曲线上存在蠕变拐点;若外界应力过大，煤体直接发生加速蠕变时曲线上不存在蠕变拐点。4个煤样受到最后1级轴压作用时的蠕变曲线如图4所示。由于瓦斯压力为0.40 MPa时，A4煤样发生了典型的非衰减蠕变，蠕变曲线上存在蠕变拐点，所以，接下来研究非衰减蠕变的蠕变模型时将重点围绕A4煤样进行论证。
图4 最后1级蠕变
Fig.4 Final creep
3 蠕变模型及对比分析
3.1 函数逼近
函数逼近是根据常用的蠕变经验方程，利用MATLAB的曲线拟合工具(CFTOOL)对实验结果进行逼近。
(1)幂逼近
CRUDEN曾论证单轴压力下岩石初始蠕变阶段符合幂函数形式，LMA采用幂函数拟合衰减蠕变，考虑到煤岩初始应变，本文对幂函数右边项加入初始应变，写成代数形式:
ε=atb+c
(1)
式中，a，b，c为系数。
将4组实验中的蠕变数据进行幂逼近，结果见表1。
表1 幂逼近
Table 1 Power approximation
由表1知，幂函数的拟合精度较高，基本在0.94以上。
对式(1)求导得到蠕变速率:
(2)
由表1拟合系数可知，发生衰减蠕变时，ab>0且b-1<0，蠕变速率
随时间延长蠕变速率减小，符合衰减蠕变的变化特征;发生加速蠕变时，ab>0且b-1>0，蠕变速率且随时间延长而增大，亦符合加速蠕变的变化特征。
若只考虑拟合精度和蠕变速率这两个因素，容易证明幂函数模型可靠，但进一步研究则会发现幂函数拟合存在着一定的不足。
当系数a，b都为正数时，对式(1)求极限:
(3)
可见，随时间延长应变不断增大，这与衰减蠕变应变先增大后趋于稳定的结论相矛盾，因此幂函数拟合衰减蠕变时稍有缺陷。而且，对幂函数求导发现
恒成立，应变单调递增，无法体现非衰减蠕变的变化特征。例如，瓦斯压力P=0.40 MPa、轴压σ1=15 MPa时的幂逼近如图5所示。
图5 非衰减蠕变的幂逼近(P=0.40 MPa,σ1=15 MPa)
Fig.5 Power approximation of non decaying creep(P=0.40 MPa,σ1=15 MPa)
(2)指数逼近
指数模型[18]是在VOIGT(沃格特)和KELVIN(开尔文)模型的基础上进行改进得到的，KELVIN模型由一个黏性元件(Newton体)和弹性元件(Hooke体)并联而成，其力学关系式为
(4)
式中，σ，η，E分别为应力、黏性系数、弹性系数。
求微分方程得:
(5)
当t=0时，ε=0。事实上，当含瓦斯煤受到三轴应力的作用时，存在初始应变，且非0。因此，将上述等式的右边添一常数项，写成代数形式:
ε=d+e[1-exp(-ft)]
(6)
式中，d，e，f均为系数。
将4组实验中的蠕变数据进行指数逼近，结果见表2。
由表2可知，拟合系数均为正数，对式(6)求极限有
(7)
说明随着时间延长，变形量逐渐增大至一定值。
对式(6)求导，得蠕变速率:
(8)
求极限:
(9)
可以发现，指数函数蠕变速率最终趋于0，能体现衰减蠕变的变化特征，但无法体现非衰减蠕变的变化特征。
(3)多项式逼近
由于幂函数和指数函数各有缺陷，故有学者提出三次多项式拟合[19]。本文对多项式的拟合次数提出质疑，认为三次多项式并非一定最优，因为多项式拟合次数不同，得到的拟合结果将会不同，下面进行不同次数的多项式拟合。
表2 指数逼近
Table 2 Exponential approximation
在利用MATLAB的CFTOOL命令进行多项式拟合时，调用LAR格式对实验数据进行降噪处理，以减小极值点或异常波动点的影响。拟合过程中结合实验散点图和拟合度选出最佳拟合次数(n)，拟合多项式通式如下:
(10)
式中，ai，a0均为系数。
对4组实验数据采用多项式逼近得到的结果见表3。
由表3可知，多项式拟合时精度均为0.99左右，整体较幂函数和指数函数好。其中，多项式拟合非衰减蠕变时，拟合曲线上存在蠕变拐点，能较好地体现非衰减蠕变的变化特征。例如，瓦斯压力0.40 MPa、轴压15 MPa时煤样发生了典型的非衰减蠕变，采用5次多项式拟合，精度高达0.999 4，具体拟合情况如图6所示。
多项式拟合精度高，但拟合次数难以确定，次数决定了拟合曲线的形态。而且，多项式拟合衰减蠕变时存在失真问题，例如，瓦斯压力0.40 MPa、轴压6 MPa时，采用5次多项式拟合，结果发现在12 h后的拟合值骤增，如图7所示。
表3 多项式逼近
Table 3 Polynomial approximation
图6 非衰减蠕变的多项式逼近(P=0.40 MPa,σ1=15 MPa)
Fig.6 Polynomial pproximation of non decaying creep(P=0.40 MPa,σ1=15 MPa)
图7 衰减蠕变的多项式逼近(P=0.40 MPa,σ1=6 MPa)
Fig.7 Polynomial approximation of decaying creep(P=0.40 MPa,σ1=6 MPa)
综上所述，多项式逼近法得到的拟合度最优，拟合非衰减蠕变时曲线上存在蠕变拐点，较能反映非衰减蠕变的变化规律，但拟合次数难以确定，且不适合用于衰减蠕变的拟合;而指数逼近和幂逼近均存在缺陷，不太适合用于衰减蠕变或非衰减蠕变的拟合。所以，接下来探讨力学模型。
3.2 力学模型
(1)西原模型及西原加速模型
西原(Nishihara)模型[20-21]是由弹性元件(H)、黏弹性元件(N/H)和理想黏塑性元件(N/St.V)串联而成。
将西原模型化成带参数的形式如下:
当σ<σs时:
ε=A+B(1-exp(-Ct))
(11)
当σ>σs时:
ε=A+B[1-exp(-Ct)]+Dt
(12)
式中，σs为屈服应力;A，B，C，D均为参数。
由于式(11)和式(6)形式一致，二者的拟合结果相同，所以西原模型拟合衰减蠕变时的结果直接参见表2。当煤样受最后1级轴压加载发生加速蠕变或者非衰减蠕变时，利用式(12)进行拟合，拟合精度介于0.88～0.99。
其中，利用西原模型拟合含加速蠕变的非衰减蠕变时，拟合曲线单调递增，图像呈下凹状，不存在蠕变拐点，无法体现加速蠕变的变化规律。例如，瓦斯压力0.40 MPa、轴压15 MPa时用西原模型拟合的结果如图8所示。
图8 西原模型拟合非衰减蠕变(P=0.40 MPa,σ1=15 MPa)
Fig.8 Nishihara model fitting non decaying creep (P=0.40 MPa,σ1=15 MPa)
西原加速模型[22]是在西原模型的基础上串联一个惯性元件，具体模型如图9所示，图中4个子系统对应的应变εe，εve，εvp，ερ分别表示弹性应变、黏弹性应变、黏塑性应变、惯性元件应变。
图9 西原加速模型
Fig.9 Extended Nishihara model
设材料进入加速蠕变阶段时的应变阈值为εa，该惯性元件满足:
(13)
由于惯性元件对应变较为敏感，当煤样的应变ε大于εa时，煤样由稳定蠕变阶段进入到加速蠕变阶段。当ε≥εa时，解微分方程得
(14)
式中，ρ为惯性质量，kg/m;C1和C2是任意常数。
所以，西原加速模型为
当σ<σs时，发生衰减蠕变，蠕变方程式和式(11)相同;
当σ>σs时，发生非衰减蠕变，蠕变方程式为
ε=εe+εve+εvp+ερ
(15)
即:
(16)
将其化成一般式如下:
ε(t)=E+F[1-exp(-Gt)]+Ht+It2
(17)
式中，E，F，G，H，I均为系数。
西原加速模型拟合衰减蠕变时的结果见表2，拟合4组煤样最后1级轴压加载发生加速蠕变或非衰减蠕变时采用式(17)，得到的结果见表4。
表4 西原加速模型拟合最后1级蠕变
Table 4 Accelerated Nishihara model fitting final creep
由表4可知，西原加速模型拟合典型的非衰减蠕变时，拟合曲线存在滞后性，效果并不理想，说明该模型的普适性存在一定缺陷。文献[22]用该模型拟合砂岩非衰减蠕变时，效果较好，但是砂岩与煤岩存在一定的差别。煤岩吸附或解吸瓦斯能力强，且煤岩内部含有复杂的孔隙结构，这些因素直接影响蠕变试验过程中煤样的应变。另外，该模型在求解微分方程式时忽略了系数项的影响，说明西原加速模型有待改进。
(2)Burgers模型和改进的Burgers模型
Burgers模型[23-24]由马可斯伟尔体和开尔文体串联而成，其中，马可斯伟尔体由Hooke体和Newton体串联而成。Burgers模型描述衰减蠕变和非衰减蠕变的方程式相同，且与式(12)形式一致。所以，采用Burgers模型拟合4组煤样最后1级轴压加载时的拟合结果与西原模型一致。通过MATLAB拟合4组煤样衰减蠕变时，拟合结果见表5。
表5 Burgers模型拟合衰减蠕变
Table 5 Burgers model fitting decaying creep
由表5可知，Burgers模型拟合衰减蠕变时精度基本在0.92以上。进一步分析，对式(12)求导:
(18)
由表5可知，系数B，C，D均为正数，说明轴向应变随时间延长不断增大，不符合衰减蠕变的变化特征，可见Burgers模型拟合衰减蠕变时依然存在缺陷。而且由图8可知，Burgers模型拟合非衰减蠕变时亦不适合，拟合曲线不存在蠕变拐点，无法体现煤样由稳定蠕变阶段过渡到加速蠕变阶段。
由于Burgers模型存在一定的缺陷，王新刚等[25]提出改进的Burgers模型:当煤体发生衰减蠕变时，用式(12)拟合;当煤体发生非衰减蠕变时，用下式拟合:
ε=J+K[1-exp(-Lt)]+Mt+Nexp(tp-1)
(19)
式中，J，K，L，M，N，p均为拟合系数。
利用改进的Burgers模型拟合4组煤样最后1级轴压加载时的加速蠕变或非衰减蠕变时，得到的结果见表6。从表6可以看到，改进的Burgers模型拟合典型的非衰减蠕变时，精度0.875，拟合效果不显著。由此可知，Burgers模型拟合及其改进模型在用于煤岩衰减蠕变和非衰减蠕变的模拟时均存在一定的不足之处。
表6 改进的Burgers模型拟合最后1级蠕变
Table 6 Improved Burgers model fitting final creep
4 西原加速模型的改进
综上所述，常用的经验方程(幂函数、指数函数、多项式形式等)都存在一定的不足，拟合情形较单一，只能拟合衰减蠕变或非衰减蠕变中的一种，而不能同时用于衰减蠕变和非衰减蠕变的拟合。
4种力学模型中，Burgers模型及其改进模型模拟衰减蠕变和非衰减蠕变时均存在一定缺陷，西原模型无法拟合非衰减蠕变，西原加速模型拟合非衰减蠕变时存在滞后性。
由于西原加速模型中惯性元件的惯性质量ρ影响着拟合曲线中轴向应变的大小，从而影响拟合曲线的形态，由此可知，滞后性问题主要由式(14)所引起。再者，结合蠕变拐点理论，煤样发生非衰减蠕变时曲线上存在蠕变拐点。因此，建议保留含ρ的二次项去掉一次项，则式(16)变为
ε=E+F[1-exp(-Gt)]+It2
(20)
为了求得蠕变拐点，两端求二阶微分得
(21)
再令
得到
说明，蠕变曲线上存在ta使得:
当t<ta时，
拟合曲线呈上凸状，蠕变速率减小，应变先增大后逐渐趋于平缓;
当t>ta时，
拟合曲线呈下凹状，蠕变速率增大，轴向应变急剧增大。
其中，蠕变拐点为(ta,εa)。于是，就得到了改进的西原加速模型:
当σ<σs时，发生衰减蠕变，蠕变方程式为
ε(t)=A+B[1-exp(-Ct)]
(22)
当σ>σs时，发生非衰减蠕变，蠕变方程式为
ε(t)=E+F[1-exp(-Gt)]+It2
(23)
为了检验改进的西原加速模型的可靠性，利用MATLAB模拟4组煤样最后1级轴压加载的加速蠕变或非衰减蠕变，得到的参数见表7。
表7 改进的西原加速模型拟合最后1级蠕变
Table 7 Improved Nishihara acceleration model fitting final creep
由表7可知，改进的西原加速模型拟合效果显著。4组煤样受最后1级轴压作用时只有A4煤样发生了典型的非衰减蠕变，通过A4煤样实验结果对比西原加速模型和改进的西原加速模型，结果如图10所示。
图10 非衰减蠕变拟合(P=0.40 MPa,σ1=15 MPa)
Fig.10 Non decaying creep fitting(P=0.40 MPa,σ1=15 MPa)
由图10可以清楚看到，西原加速模型的拟合曲线在前0.05 h内，曲线上凸、应变快速增大;在0.05～0.25 h内，应变缓慢增大;0.25 h后，曲线下凹、应变呈增大的趋势，而实际蠕变曲线在0.2 h点后应变急剧增大，其中，蠕变拐点为(0.2,0.069)。西原加速模型拟合非衰减蠕变时具有明显的滞后性，而改进的西原加速模型成功克服了滞后性问题，能真实地模拟出煤样发生非衰减蠕变时的变化特征。
5 讨 论
文献[26]曾提到关于瓦斯压力对煤岩蠕变会造成影响的猜想，本实验结合余吾矿实际，研究了4组低瓦斯压力下煤样受三轴应力作用时的蠕变规律。结果表明煤岩蠕变与瓦斯压力呈正相关，煤岩内部瓦斯的存在使得煤岩应变增大。但是，由于现场实际，本文没有进行高、中、低瓦斯压力下的蠕变实验，在后续研究中可以开展高瓦斯压力的实验，进一步研究瓦斯压力对煤岩蠕变的实验，并在蠕变模型中加入瓦斯压力参数。
当煤岩所受载荷大于长期强度时，会发生失稳破坏，但至今任没有一种较好的方法来确定[27]。文献[15]也曾指出，煤岩发生失稳破坏无明显征兆，加速蠕变阶段与稳定蠕变阶段不明显。从本文实验结果可以发现，A1～A3煤岩在最后1级轴压之前发生稳定蠕变，最后1级轴压加载时立即出现加速蠕变，其失稳破坏前无明显前期征兆。这些特点不仅危及煤层气开采过程中的钻进安全，而且还有可能出现煤与瓦斯突出等煤岩动力灾害。在现场实际过程中，煤岩长期强度难以确定，若不能预测煤岩是否会发生失稳破坏或者发生失稳破坏的时间，无疑对安全作业构成极大威胁。设想若先找到煤岩蠕变的数学模型，再求算蠕变拐点(ta,εa)，若不存在蠕变拐点则不会发生失稳破坏，若存在蠕变拐点则会发生非衰减蠕变，进而求出拐点出现的时间以及对应的轴向应变，通过应变的变化趋势就能预判煤岩是否会发生失稳破坏。因为，煤岩发生非衰减蠕变与蠕变曲线上至少存在1个蠕变拐点是互为充要条件，这在上面已经证实过了。所以，蠕变拐点的提出不仅加强了蠕变的数学理论，还能为煤岩动力灾害、巷道围岩稳定性及水电站坝基岩石稳定性等工程实践提供一定的指导。
6 结 论
(1)岩石发生非衰减蠕变是蠕变曲线上至少存在1个蠕变拐点的充要条件。蠕变拐点理论不仅完善了蠕变加速度的理论，而且为煤岩动力灾害预警等工程实践提供了一条新思路。由此展望，若能预先判断出蠕变拐点(ta,εa)出现的条件，并采取卸压等防护措施就有可能避免煤岩失稳破坏。
(2)煤岩蠕变与瓦斯压力呈正相关，瓦斯压力从0.13 MPa增大至0.40 MPa时，煤样轴向应变随之呈小幅度增大。说明，余吾矿采取预抽采瓦斯可预防煤与瓦斯突出等动力灾害的发生。
(3)各种常见的蠕变模型中，经验函数(多项式、幂函数、指数函数等)均存在些许缺陷，只适用于衰减蠕变或非衰减蠕变，不能同时用于模拟衰减蠕变和非衰减蠕变的拟合。西原模型不适合用于模拟非衰减蠕变，Burgers模型及其改进模型在拟合煤岩衰减蠕变和非衰减蠕变时均存在不足，改进的西原加速模型拟合煤岩非衰减蠕变时存在滞后性。
(4)本文提出改进的西原加速模型克服了滞后性问题，可以用于煤岩的衰减蠕变和非衰减蠕变的模拟。结合蠕变拐点理论，先判断出煤岩是否会发生失稳破坏，进而求出蠕变拐点出现的时间以及对应的应变值，能为预防煤岩失稳破坏提供一定的指导。
参考文献(References):
[1] 何学秋.煤岩流变电磁动力学[M].北京:科学出版社,2003:1-7.
[2] 陈占清.采动岩体蠕变与渗流耦合动力学[M].北京:科学出版社,2010:1.
[3] 耶格.岩石力学基础[M].北京:科学出版社,1981:382-385.
JAEGER.Fundamentals of rock mechanics[M].Science Press,1981:382-385.
[4] 曹树刚,鲜学福.煤岩蠕变损伤特性的实验研究[J].岩石力学与工程学报,2001,20(6):817-821.
CAO Shugang,XIAN Xuefu.Experimental study of rock creep damage characteristics of coal[J].Chinese Journal of Rock Mechanics and Engineering,2001,20(6):817-821.
[5] MARANINI E,BRIGNOLI M.Creep behaviour of a weak rock.Experimental characterization[J].International Journal of Rock Mechanics & Mining Sciences,1999,36(1):127-138.
[6] CRUDEN D M.The form of the creep law for rock under uniaxial compression[J].International Journal of Rock Mechanics & Mining Science & Geomechanics Abstracts,1971,8(2):105-126.
[7] MA L,DAEMEN J J K.An experimental study on creep of welded tuff[J].International Journal of Rock Mechanics & Mining Sciences,2006,43(2):282-291.
[8] WANG D,WEI J,YIN G,et al.Triaxial creep behavior of coal containing gas in laboratory[J].Procedia Engineering,2011,26(4):1001-1010.
[9] 阿克洛尼斯,李仲伯,李怡宁.聚合物粘弹性引论[M].北京:宇航出版社,1984:165-185.
[10] 阎岩,王思敬,王恩志.基于西原模型的变参数蠕变方程[J].岩土力学,2010,31(10):3025-3035.
YAN Yan,WANG Sijing,WANG Enzhi.Creep equation of variable parameters based on Nishihara model[J].Rock and Soil Mechanics,2010,31(10):3025-3035.
[11] 邓荣贵,周德培,张倬元,等.一种新的岩石流变模型[J].岩石力学与工程学报,2001,20(6):780-780.
DENG Ronggui,ZHOU Depei,ZHANG Zhuoyuan,et al.A new rheological model for rocks[J].Chinese Journal of Rock Mechanics and Engineering,2001,20(6):780-780.
[12] 徐鹏,杨圣奇.循环加卸载下煤的黏弹塑性蠕变本构关系研究[J].岩石力学与工程学报,2015,34(3):537-545.
XU Peng,YANG Shengqi.Study on the constitutive relation of viscoelastic plastic creep under cyclic loading and unloading[J].Chinese Journal of Rock Mechanics and Engineering,2015,34(3):537-545.
[13] 尹光志,王登科,张东明,等.两种含瓦斯煤样变形特性与抗压强度的实验分析[J].岩石力学与工程学报,2009,28(2):410-417.
YIN Guangzhi,WANG Dengke,ZHANG Dongming,et al.Test analysis of deformation characteristics and compressive strengths of two types of coal specimens containing gas[J].Chinese Journal of Rock Mechanics and Engineering,2009,28(2):410-417.
[14] 阎岩.渗流作用下岩石蠕变试验与变参数蠕变方程的研究[D].北京:清华大学,2009:36-37.
YAN Yan.Research on rock creep tests under seepage flow and variable parameters creep equation[D].Beijing:Tsinghua University,2009:36-37.
[15] 范翔宇,张千贵,艾巍,等.煤岩储气层岩石蠕变特性与本构模型研究[J].岩石力学与工程学报,2013(S2):3732-3739.
FAN Xiangyu,ZHANG Qiangui,AI Wei,et al.Study on creep characteristics and constitutive model of coal and gas reservoir rock[J].Chinese Journal of Rock Mechanics and Engineering,2013(S2):3732-3739.
[16] MARANINI E,BRIGNOLI M.Creep behaviour of a weak rock experimental characterization[J].International Journal of Rock Mechanics & Mining Sciences,1999,36(1):127-138.
[17] PING C,WEN Y,WANG Y,et al.Study on nonlinear damage creep constitutive model for high-stress soft rock[J].Environmental Earth Sciences,2016,75(10):900.
[18] FENG J,ZHANG C,GANG W,et al.Creep modeling in excavation analysis of a high rock slope[J].Journal of Geotechnical & Geoenvironmental Engineering,2003,129(9):849-857.
[19] 陈绍杰,郭惟嘉,杨永杰.煤岩蠕变模型与破坏特征试验研究[J].岩土力学,2009,30(9):2595-2598.
CHEN Shaojie,GUO Weijia,YANG Yongjie.Experimental study on creep model and failure characteristics of coal rock[J].Rock and Soil Mechanics,2009,30(9):2595-2598.
[20] 沈才华,张兵,王文武.一种基于应变能理论的黏弹塑性蠕变本构模型[J].岩土力学,2014,35(12):3430-3436.
SHEN Caihua,ZHANG Bing,WANG Wenwu.A visco elastic plastic creep constitutive model based on strain energy theory[J].Rock and Soil Mechanics,2014,35(12):3430-3436.
[21] 杨逾,李盈,周小科.基于西原加速模型的煤岩蠕变试验研究[J].煤炭学报,2014,39(11):2190-2194.
YANG Yu,LI Ying,ZHOU Xiaoke.Study on creep test of coal and rock based on the original acceleration model[J].Journal of China Coal Society,2014,39(11):2190-2194.
[22] JIANG Q,QI Y,WANG Z,et al.An extended Nishihara model for the description of three stages of sandstone creep[J].Geophysical Journal International,2013,193(2):841-854.
[23] 杨振伟,金爱兵,周喻,等.伯格斯模型参数调试与岩石蠕变特性颗粒流分析[J].岩土力学,2015,36(1):240-248.
YANG Zhenwei,JIN Aibing,ZHOU Yu,et al.Analysis of Burgers model parameters and particle flow analysis of rock creep[J].Rock and Soil Mechanics,2015,36(1):240-248.
[24] 熊良宵,杨林德,张尧.岩石的非定常Burgers模型[J].中南大学学报(自然科学版),2010,41(2):679-684.
XIONG Liangxiao,YANG Linde,ZHANG Yao.Unsteady Burgers model of rock[J].Journal of Central South University(Science and Technology),2010,41(2):679-684.
[25] WANG X,HU B,TANG H,et al.A constitutive model of granite shear creep under moisture[J].Journal of Earth Science,2016,27(4):677-685.
[26] 尹光志,王登科,张东明,等.含瓦斯煤岩三维蠕变特性及蠕变模型研究[J].岩石力学与工程学报,2008,27(S1):2631-2636.
YIN Guangzhi,WANG Dengke,ZHANG Dongming,et al.Research on triaxial creep properties and creep model of coal containing gas[J].Chinese Journal of Rock Mechanics and Engineering,2008,27(S1):2631-2636.
[27] 杨圣奇.裂隙岩石力学特性研究及时间效应分析[M].北京:科学出版社,2011:267.
10