温度-围压-瓦斯压力作用下煤岩力学性质及有限变形行为

高亚楠，高 峰，谢 晶，滕 腾，闫伟城，高 琳

(1.中国矿业大学 深部岩土力学与地下工程国家重点实验室，江苏 徐州 221116; 2.中国矿业大学 力学与土木工程学院，江苏 徐州 221116; 3.四川大学 水利水电学院，四川 成都 610065; 4.中国矿业大学(北京) 能源与矿业学院，北京 100083)

煤炭是我国的主体能源，2019年在我国一次性能源消费结构中约占58%，预计到2050年仍将占到40%以上[1-2]。巨大的消耗导致浅部煤炭资源逐渐枯竭，深部开采已成为国家的战略需求[3]。目前，我国共有千米以上煤矿矿井约50座，最大开采深度已超过1 500 m，并且仍以10～25 m/a的速度延伸[4-5]。深部开采过程中，矿山压力的控制、煤岩体的控稳、灾害的防控，较浅部相比更加困难，其根本在于深部煤岩体变形破坏更加复杂，且机制不明确，使得工程实践上存在一定的盲目性、低效性和不确定性[6]。因此，必须对开采环境下，煤岩体的力学行为、特征及规律开展深入研究。

开采过程中，影响煤岩力学性质的环境因素主要有地应力、瓦斯压力以及地温，煤岩体力学行为最大的特征为大变形[7-8]。针对这些方面，国内外学者开展了大量的研究。实验室内往往采用围压来模拟。地应力对煤岩力学特性影响的研究，随着围压的增加，煤岩除了弹性模量、泊松比、峰值强度会发生变化之外[9-10]，其变形破坏特征和机制、变形破坏过程中的能量转移转化规律等也会发生相应改变，比如由低围压、低应力作用下的小变形破坏逐渐转化到高围压、高应力作用下的塑性以及大变形破坏[11-14]。瓦斯及瓦斯压力的存在，对于煤岩的微结构、孔隙压力所产生的影响，将会在宏观上体现为力学性质、变形特征以及破坏方式等方面，赵洪宝等[15]通过含瓦斯煤岩的三轴压缩试验分析了瓦斯对煤岩强度和破坏方式的影响。谢广祥等[16]建立了含瓦斯煤采动应力-瓦斯压力耦合作用的数学模型，并阐明了瓦斯压力与采动应力协同演化的力学机理。YANG等[17]以连续介质损伤力学和细观损伤理论为基础，考虑了瓦斯压缩能，建立了含瓦斯煤非线性损伤模型。CHEN等[18]基于分形几何理论研究了含瓦斯煤压出过程中表面裂纹的演化规律。温度的变化会改变煤岩内部结构，进而带来煤体力学性质的改变[19]。许江等[20]借助自行研发的多物理场耦合试验系统，指出20～80 ℃内，温度对含瓦斯煤的变形及力学特性在不同的温度区间内分别有不同的变化规律。目前，现场观测和实验室试验都已证明，温度升高会造成煤体基质热开裂和物质热挥发现象。NAKAGAWA等[21]指出温度变化对煤体微观结构存在重要影响。YU等[22]通过CT扫描手段，研究了温度对煤体中细小孔洞和缺陷的影响。

如前所述，目前采动煤岩体极其显著的力学特征就是非线性变形，经典小变形理论对深部煤岩行为无法准确描述[7]，因此，对采动煤岩力学行为的认识必须建立在非线性大变形的框架下。康红普等[23]从围岩物性劣化、偏应力诱导围岩扩容、软岩结构性流变及采动影响等方面,揭示了采动应力和采场高地应力耦合作用下巷道围岩大变形机理;何满潮等[24-25]针对煤矿开采沉陷问题介绍了有限变形理论处理大变形问题的优越性，给出了非线性大变形理论在软岩巷道工程中的应用实例;GUZEV等[26]基于非线性几何场推导了深部围岩位移场的解析解，并且指出非线性几何场内变量是影响弹性位移协调性的决定因素;PARK[27],VRAKAS[28]等采用有限应变理论，计算了围岩变形特征曲线，并指出处理大变形问题时，小变形在方法与理论方面的误差;局部转动不协调是材料大变形破坏的诱因之一[29]，文献[30-31]基于有限变形理论，以局部转动为出发点，建立了岩石的大变形本构方程，并发展了相应的数值计算方法。

目前，温度、地应力(围压)、瓦斯压力等因素对煤岩力学性质影响的研究已经十分丰富，但是同时也可以发现，综合考虑上述因素的研究并不多见;与此同时，有关煤岩大变形的研究大都聚焦在工程尺度和位移场等方面，对于非线性力学框架下煤岩大变形的理论研究相对较少。基于此，笔者借助室内试验，设计了温度-围压-瓦斯压力作用下的煤岩力学试验，旨在探究温度、围压和瓦斯压力对煤岩力学性质的影响程度和规律，进而借助有限变形理论，描述煤岩的大变形行为，引入平均整旋角作为内变量，建立多物理场作用下的煤岩大变形本构模型，为煤炭资源开采设计和煤岩控稳提供理论基础。

1 温度-围压-瓦斯压力作用下煤岩力学试验

1.1 试样与试验设备

本次试验煤样全部采自平煤矿区，工作面埋深580～705 m，瓦斯压力1.5～2.0 MPa，煤岩变形行为复杂(图1)。对煤样进行XRD衍射分析可知[10]，煤样主要由高岭石、石英以及方解石等矿物组成(图2)。按照煤岩力学测试标准[32]，将煤岩样品在实验室钻取切割成φ50 mm×100 mm的标准试样，之后用磨平机和砂纸对试样两端仔细打磨，使上下表面平行度在±0.5 mm以内(图3)。

图1 现场巷道片帮大变形Fig.1 Large deformation of rib spalling of roadway

图3 样品采集加工Fig.3 Sample collection and preparation

试验使用国产TAW-2000型微机伺服岩石试验系统(图4)。该试验系统主要由轴压加载系统、围压加载系统、伺服控制系统以及温度控制系统等部分组成。轴压加载速率0.01～100 kN/s;围压加载速率0.01～0.50 MPa/s;变形测量装置测量精度:轴向0.1 μm、径向0.1 μm;温度加载采用外加热装置，波动范围±1 ℃。

1.2 试验设计

基于工作面开采环境，结合学者在试验方面的经验，试验考虑了3个因素，即温度、围压、瓦斯压力，每个因素选取3个水平，具体为[30,33-34]:① 温度选取20,40和80 ℃三个水平，以表示常温、井下温度、以及温度异常区的预测[30,33];② 围压选取0，4和8 MPa三个水平，以表示单轴压缩状态，中等围压水平以及高围压水平[34];③ 瓦斯压力选取0，1和3 MPa三个水平，以表示无瓦斯压力状态，低瓦斯压力水平，和高瓦斯压力水平[34]。

由于试验需要考虑到3个因素，每个因素3个水平，试验量较大，因此，笔者采用正交试验法[35]，将试验分为9组，每组3个试样，试验设计见表1。试验组中的名称含义:以试样M80-3-4为例，80为试验温度，℃;3为瓦斯压力，MPa;4为围压，MPa。

图4 试验设备及示意Fig.4 Experiment facility and the illustration

表1 试验设计Table 1 Experiment design

试验开始前，预先将煤样置于80 ℃恒温箱内烘干24 h，之后置于瓦斯压力(1或3 MPa)环境中24 h，最后将试样取出按如下过程进行试验，过程为① 将试样两端紧贴压头，并用热缩管密封试样及试样与压头连接处，防止加载介质(液压油)在试验过程中进入试样，同时防止试样破坏后产生的碎屑污染油缸;② 在试样上安装环向变形传感器和轴向变形传感器，并将试样固定在压力室内，准备加压(图4);③ 通过轴压系统，对试样施加1 kN的预应力，以防止加载围压过程中出现试样位置偏移，加载过程选用力控模式，加载速率为50 N/s;④ 开启围压充液，直至压力室充满加载介质，按位移方式加载至目标围压(4或8 MPa)，加载速率为10 mm/min;⑤ 开启轴压加载，选用变形控制模式以0.02 mm/min的速率加载直至试样破坏;⑥ 试验结束后，依次卸载轴压和围压，之后打开回油阀开始回油，待压力室内的加载介质清除后，取出试样，拆除变形传感器并将破坏试样保存。

这里需要补充的是，在试验过程中还需要施加不同温度，对于20 ℃的试验条件，主要依靠室内恒温调整，保持环境温度为20 ℃，同时将加热装置的工作温度调整为20 ℃，作为辅助温度控制;对于40 ℃和80 ℃的试验条件，主要依靠外加热装置保持试验温度。

2 温度-围压-瓦斯压力作用下煤岩力学试验结果分析

图5为部分试样的应力-应变曲线(其中,M20-0-0-1中，1为试验序号，下同)，由图5可知，应力-应变曲线基本都可以分为压密阶段、弹性阶段、弹塑性阶段、屈服破坏阶段。其中压密阶段在高围压时不明显，这是因为围压加载过程中，对煤岩试样有一定的压实作用。不同试验条件下的试样不仅在压密阶段和弹性阶段的变形特征有着明显的区别，而且峰后的应力-应变曲线也存在一定差异，部分试样在峰后屈服时变形较大，表现出较好的残余承载能力，而有些试样峰后曲线则具有较为明显的跌落形态，表现出较强的脆性，可以看出不同试验条件下的煤岩力学参数存在明显差别。因此，下面将讨论各个因素对煤岩力学性质的影响。

2.1 各因素对煤样力学性质影响的方差分析

由试验结果可知，温度、围压和瓦斯压力对煤样力学性质影响明显。为了进一步厘清各因素对煤样力学性质和力学行为的影响规律，基于正交试验设计思想和统计学方法对温度、瓦斯压力以及围压引起的煤岩各力学参数变化进行方差分析[35]，得出各因素对煤岩峰值应力、弹性模量及泊松比影响的程度和显著性。

图5 煤样应力-应变曲线Fig.5 Stress-strain curves of the coal specimens

由方差分析结果(表2～4)可知，温度、围压以及瓦斯压力对煤岩力学性质有显著影响。由表3～5可知，温度对弹性模量E、峰值应力σp和泊松比ν影响的显著性水平均为95%，围压对这3个力学参数影响的显著性水平也为95%，瓦斯压力对峰值应力和弹性模量影响的显著性水平为95%，对泊松比影响的显著性水平为90%。

表2 各因素对弹性模量影响的方差分析Table 2 Analysis of variance of elastic modulus with all factors

值得注意的是，尽管获得了各因素对煤岩力学参数的影响显著性，但是仍属于定性分析，还不能确定各个因素对力学参数的具体影响规律，因此，2.2节将尝试分析温度、围压以及瓦斯压力对煤岩弹性模量，峰值强度以及泊松比的量化影响。

2.2 各因素对煤岩力学性质的影响规律

温度、围压、瓦斯压力等因素对煤岩物理力学性质影响的定量研究中，大多数学者采用了线性模型，或者是类线性模型[30]，文献[36]研究了在不同围压和不同瓦斯压力下型煤和原煤的变形特征，并建立了强度、弹性模量与瓦斯压力间的线性模型。文献[37]基于含瓦斯煤力学试验，发现了瓦斯压力与弹性模量以及泊松比存在线性关系。文献[38]认为煤岩弹性模量随围压变化呈现非线性变化规律，并用多项式来描述弹性模量随围压的变化规律，但给出数学模型中高次项系数比一次项系数小将近2个数量级。文献[20]指出，可以用线性模型描述温度变化对含瓦斯煤强度和弹性模量的影响，同时给出了泊松比随温度变化的分段线性模型。因此，结合前人经验，利用多元线性拟合的方法分析多因素影响下煤岩力学参数变化规律。

表3 各因素对峰值应力影响的方差分析Table 3 Analysis of variance of peak stress with all factors

表4 各因素对泊松比影响的方差分析Table 4 Analysis of variance of Poisson’s ratio with all factors

表5 弹性模量线性拟合结果Table 5 Linear fitting results of elastic modulus

由表5～7可以看出弹性模量E、峰值应力σp、泊松比μ与各因素的关系可以用线性模型来表示:

E(T,C,G)=0.013T+0.131C+0.434G+1.346

(1)

σp(T,C,G)=0.029T+1.28C+2.436G+5.477

(2)

μ(T,C,G)=0.001T+0.006C+0.009G+0.214

(3)

表6 峰值应力线性拟合结果Table 6 Linear fitting results of peak stress

表7 泊松比线性拟合结果Table 7 Linear fitting results of Poisson’s ratio

由表5和式(1)可知，瓦斯压力、围压和温度对弹性模量影响权重依次为0.434,0.131和0.013，均为正相关。各因素与弹性模量线性相关的显著性水平均为95%。

由表6和式(2)可知，瓦斯压力、围压和温度对峰值应力的影响权重依次为2.436,1.280和0.029，均为正相关。瓦斯压力和围压与峰值应力线性相关的显著性水平为95%，温度为90%。

由表7和式(3)可知，瓦斯压力、围压和温度对泊松比影响权重分别为0.009,0.006和0.001，均为正相关。温度、围压、瓦斯压力与泊松比的线性相关的显著性水平为90%。同时需要注意的是，表4中温度和围压对泊松比的影响显著性水平为95%，这2处显著性水平不同的原因在于线性模型的误差。

就温度而言，对煤岩刚度和强度的影响，现有的研究中可以分为正相关、负相关和阈值现象[39-41]，而在本文的研究范围内，温度所产生的热应力没有超过基质颗粒间的拉伸强度，同时又导致基质颗粒膨胀，使得基质更加致密、结构更加完整，进而提高了承载能力和刚度;就围压而言，目前大多的研究结论是:围压的升高，可以提高煤岩体的刚度及强度[42]，这是由于静水压力的约束，抑制了煤岩内微裂纹和孔隙等天然损伤的发育，从而提高了承载能力和刚度。就瓦斯压力而言，经过瓦斯保压的煤样，其孔隙间的瓦斯，会产生孔隙压力，当煤体受到外载荷作用时，孔隙压力会降低骨架受到的有效应力，从而表现出承载能力和刚度的提高。同时值得注意的是，从本文的各个线性模型中可以看到，温度的权重系数最小，围压的权重系数居中，瓦斯压力的权重系数最大，或许这也意味着与温度和地应力相比，瓦斯压力对煤体基质微结构和有效应力的改变，是影响煤岩变形破坏更为关键的因素，高瓦斯压力矿井，煤岩变形行为更加复杂，煤岩体控稳难度更大。

3 温度-围压-瓦斯压力作用下煤岩变形过程的有限变形

3.1 平均整旋角及在三轴试验中的计算方法

开采条件下，煤岩最典型的变形特征就是非线性，经典的小变形理论并不能准确描述，必须采用大变形理论对其力学行为进行分析，尽管试验室内，常规三轴加载路径较为简单，但是煤岩体是含天然损伤及夹杂的非均匀材料，即便是在对称载荷的作用下，依然会呈现如局部变形等在内的复杂力学行为(图6)。因此，需要寻找或建立能够描述复杂变形过程的内变量或者指标，对煤岩的非线性变形行为进行深入研究。建立在有限变形理论基础上的物理量—平均整旋角伴随着材料由变形到破坏的全过程，已经被证明可以用于描述岩石类材料的局部变形等复杂的变形过程[30，43]。

图6 含瓦斯煤岩的局部变形Fig.6 Localized deformation of the coal specimens

下面将结合有限变形理论和三轴压缩试验条件，给出平均整旋角的计算方法。陈至达[29]在拖带坐标系下建立了S-R分解定理，解决了应用经典有限变形极分解理论中存在的左极分解和右极分解不惟一的问题，根据和分解定理，有

(4)

(5)

试验测量时应变分量是基于小变形理论计算得出的，即

(6)

结合小变形应变分量和有限变形应变分量可得

(7)

根据胡克定律，一维弹性无损的本构方程为

σ=ES

(8)

式中，σ为应力；S为应变。

由式(7),(8)可得

(9)

将式(8)推广至三维情况则可以写为

(10)

(11)

在三轴压缩试验中有如下假设[30]:

(12)

则式(11)可以写为

(13)

由式(12)可得

(14)

3.2 平均整旋角的计算结果与分析

表8列出了根据3.1节中式(14)计算的平均整旋角结果，图7为平均整旋角随应力、应变的变化规律。

表8 平均整旋角计算结果Table 8 Summary of the mean rotation angle

由图7可知，平均整旋角随应力、应变的变化过程大体上可以分为快速增长区、平缓增长区和加速增长区(下降区M20-0-0-2和M40-1-8-2)。快速增长区约为轴向应力加载至峰值应力40%之前的阶段，此阶段内，试样处于压密阶段，原生微孔隙、空洞等结构逐渐闭合，试样产生较多的局部变形，平均整旋角增长较快。平缓增长区对应着轴向应力由40%峰值应力加载至70%峰值应力阶段，此阶段内，试样处于弹性阶段，变形过程较为协调稳定，平均整旋角增长也较为平缓。加速增长区则为轴向应力由70%峰值应力加载至接近破坏阶段，此阶段内裂纹发育贯穿较快，局部变形量增加较快，因而，平均整旋角也加速增长。需要注意的是M20-0-0-2和M40-1-8-2两个试样在试验过程中存在下降区，其原因在于:M20-0-0-2试件在峰前由于局部破坏出现了应力小幅跌落，这个跌落对试件整体有一定“卸载”作用，此时平均整旋角也就相应有所回落，而M40-1-8-2在峰前出现了变形的局部协调，因而平均整旋角在数值上有所减小，这2个现象可以在3.3节中图8(a)和图8(e)中看到。

切线模量的变化是材料在受力过程中刚度和脆塑性变化的标志。通过5%～95%峰值应力(每点之间间隔大小为峰值应力的10%)数据来求得切线模量，并绘制出切线模量ET随平均整旋角变化趋势图(图9)。如图9所示，切线模量随着平均整旋角的上升先增加后减小。这是因为，轴向压力加载至75%之前，试样经历了压密和弹性阶段，此时，试样的刚度持续增加，是弹性模量的硬化阶段[45];切线模量的跌落位置起始于峰值应力的75%～85%的部分，在此阶段，试样逐渐由弹脆性进入塑性状态，同时，煤岩试样内部裂纹扩展和局部变形加速，平均整旋角增长加快，试样刚度也开始下降。

图8 模型结果与试验数据对比Fig.8 Comparison of the model results and the experimental data

图9 切线模量随平均整旋角变化趋势Fig.9 Evolution of tangent modulus with the mean rotation angle

3.3 基于平均整旋角的煤岩本构模型

在岩石本构模型方面的研究，国内外学者采用引入应变、弹性模量、声发射变量等物理量作为损伤变量的方法来进行描述[46-48]。值得商榷的是:以往的模型中损伤变量随着外载荷的增加逐渐增大，也就是说材料在加载过程中处于不断弱化的过程，然而煤岩体在受载过程中并非如此，比如压密阶段，煤岩基质结构不断致密，实际上是承载和抗变形能力的强化过程，由图8可知，轴向应力加载至峰值应力75%之前的阶段，切线模量的增加也可以说明这一点;与此同时，损伤变量如何选取、损伤模型的演化是否可以反映材料变形破坏过程的本质，也尚无定论。

由3.2节的分析可知，平均整旋角伴随着煤岩受载-变形的全过程，记录着煤岩变形的“信息”，因此，这里选取平均整旋角作为内变量，建立煤岩本构模型如下:

(15)

在建立岩石类材料本构模型时，Weibull分布函数是应用最广泛的分布函数之一[49]，然而，Weibull分布是根据最弱环模型建立的，实际上该模型应用到岩石类材料时，将裂纹启裂和扩展等同，有异于实际岩石压缩破坏过程[45]。根据上述分析以及试验数据，通过半经验法，选用高斯模型[30]:

(16)

其中,y0为补偿值;A为幅度参数;θc为位置参数;w为形状参数，根据平均整旋角和应变曲线确定。根据试验数据，将式(16)中的参数拟合后汇总于表9中。

表9 模型拟合参数Table 9 Summary of the parameters of the model

将弹性模量、泊松比随各因素的变化规律以及式(16)代入式(15)中得出温度、围压、瓦斯压力作用下的煤岩本构方程:

(17)

本文提出的非线性本构模型基于有限变形(大变形)理论框架，突破了经典柯西应变的假设。由图9试验数据与模型计算的对比可知，本文建立的模型能较好的反应煤岩压密阶段、弹性阶段和屈服阶段3个阶段。这也是以往的损伤本构模型很难从物理意义上实现的，传统的本构模型建立方式，大多都是基于损伤理论/等效损伤模型而建立的，这类损伤模型的建立中，必须引入损伤因子，且该因子的变化从0到1(0为无损伤状态，1为破坏状态)[50-51]，换言之，在此框架下，煤岩在受载过程中，始终处于刚度、强度弱化过程。然而，通过试验和经验可知，煤岩样品必然经过压密、弹性、再到屈服破坏阶段，实际上从应力应变曲线中可以看到，由压密到弹性阶段的转化，是一个刚度硬化的过程[45]，因此，传统的本构模型建立方式也是与实际试验有所出入的。本文的模型引入非线性变形内变量，可以准确反映煤岩峰前变形的各个阶段和行为特征。

对于峰后破坏的行为，从严格的力学意义上来讲，不能称之为变形。因为，峰后的煤岩经过裂隙扩展、贯通，已经破碎为若干个离散体，而变形的概念建立在连续介质的基础上[52]。目前学术界尚不存在对多连体、非连续体的力学本构模型的概念和定论。本文的模型是基于非线性连续介质力学而建立，因此，也需要遵守连续介质力学的基本概念、假设和研究范围。煤岩在峰后破坏后，由连续体转变为散体结构，超出了连续介质力学的概念和范畴，不满足连续介质力学的假设，因此，在此基础上计算得出的平均整旋角不具有可靠性，本文对于本构模型的研究，也就聚焦在峰前范围内。此外，需要说明的是，虽然煤、岩等脆性材料破坏时发生的变形往往相比较塑性材料量级要小很多，甚至依然是“小变形”范围内，但是运用有限变形理论来研究煤、岩材料的变形破坏，其意义和价值不仅仅在于变形量的大小，还在于对变形的复杂和非线性程度的刻画。另外，小变形也是大变形或者有限变形经过简化的特殊情况，因此，在研究小变形问题时采用有限变形理论并不冲突。与此同时，已有研究证明即便是小变形问题，也是常常遇到非线性的变形问题，运用有限变形理论来研究，会更加准确有效[29，31]。

4 结 论

(1)温度对弹性模量、峰值应力和泊松比影响的显著性水平均为95%，围压对这3个力学参数影响的显著性水平也为95%，而瓦斯压力对峰值应力和弹性模量影响的显著性水平为95%，对泊松比影响的显著性水平为90%。各因素对力学参数的影响规律可以用线性模型来描述。

(2)平均整旋角随应力应变的变化基本可以划分为快速增长区、平缓增长区和加速增长区，这3区的出现也是煤岩压密、弹性、屈服破坏3阶段的力学过程体现;切线模量随平均整旋角可以分为增长阶段和跌落阶段，反映了煤岩从压密、硬化到脆塑性转化的过程。

(3)由于平均整旋角的变化紧密关联着煤岩变形的过程，故将其作为内变量，建立了温度、围压以及瓦斯压力作用下煤岩大变形本构方程，通过和试验数据的对比可知，笔者提出的模型可以描述煤岩变形破坏的整个过程。笔者提出的模型是基于有限变形(连续介质力学)理论，然而峰后的试样往往破碎成多个块体，属于非连续体力学，故本文模型没有涉及。