循环加卸载条件下考虑温度的煤体本构模型

周宏伟 ， 侯 魏 ， 张龙丹 ， 谢森林 ， 贾文豪 ， 史贝贝

(1.中国矿业大学(北京) 能源与矿业学院, 北京 100083；2.中国矿业大学(北京) 煤炭精细勘探与智能开发全国重点实验室, 北京 100083；3.中国矿业大学(北京) 力学与土木工程学院, 北京 100083)

深部煤炭资源赋存条件除高地应力外，同时伴随着高地温的情况[1]。深部煤层赋存环境中所存在的高地温会显著影响煤体强度、变形特征以及围岩的稳定性[2-5]。其次，在对深部煤体的开采过程中会受到长期的采动应力作用影响，其循环加卸载作用特点导致煤体的力学行为相较于单调的单轴或者三轴加载或者卸载条件会更加地复杂多变[6-8]。因此研究处于不同温度的深部煤体在循环加卸载作用下的变形破坏规律及演化特征对深部煤层的稳定性分析和绿色安全开采具有重要的理论意义和工程参考价值。

深部工程作业中，岩石将受到热-力学耦合的复杂作用[2,9-10]。冯子军等[11]认为热力耦合效应促进了无烟煤煤样基质活化，更易于进入塑性区。低于临界温度时，温度对煤样的变形以热膨胀为主，而高于临界温度值后，由于高温而挥发的气体使得煤样产生更大的变形。而周长冰等[12]发现有烟煤在热力耦合作用下的变形机制与无烟煤有所不同。薛东杰等[13-14]通过热力耦合作用下的岩样声发射特征指出了岩样强度随着温度升高呈指数型降低的趋势，另外细观有限元模拟结果显示出应变能及耗散能随着温度的升高而增大。WANG 等[6]对热力耦合下的花岗岩力学行为进行了实验研究，发现随着温度的升高，岩石的弹性模量、裂纹损伤应力和破坏强度等力学特性略有降 低。 XU 等[15]、 YIN 等[16]和WANG 等[17]利 用Weibull 分布和连续损伤理论建立岩石的热-力损伤本构模型，拟合曲线显示出峰值应力对应的岩石损伤值和应变能释放率均随围压和温度的增加呈指数增长，同时XU 等[18]基于Weibull 分布和Lemaitre 应变等效原理建立的花岗岩损伤本构模型表明温度和围压对岩石试件的延性有增强作用。综上研究结果均显示出热损伤随温度的升高而持续增长，温度场的存在促进了岩石裂纹发育。

在循环加卸载条件下煤体的损伤机理研究工作越来越多，包括损伤本构模型和能量演化等方面。彭瑞东等[19]和李涛等[20]分别设置不同的轴向峰值荷载对煤体和冻结砂岩进行重复加卸载，通过变换围压发现弹性模量的变化并不是单一的。彭瑞东等[19]基于此定义了以耗散能为主的损伤变量，发现损伤随着应力单调增大。基于此，王向宇等[21]对饱和深部煤样展开了循环加卸载实验，引入修正定义的损伤变量有效表征了弹塑性煤体的损伤演化程度。杨磊等[22]以煤岩组合体为研究对象，发现了峰前阶段能量以弹性能为主。孟庆彬等[8]指出基于弹性模量恒定的能量计算方法的不足，建立了以应力-应变曲线面积为准的能量表征方法。为了对峰后阶段岩样的本构关系有更深刻的认识，杨小彬等[23]组建了峰后重复加卸载实验，通过对岩样裂纹演化规律的研究，建立了摩擦耗散能以表征结构面之间的错动、滑移而损失的能量。陈结等[24]利用声发射对含有瓦斯的煤样进行了多种应力路径下的加卸载实验，总结了气体吸附下煤样的能量密度演化过程，并指出了弹性能与总输入能之间的线性关系。

由于分数阶导数对描述非线性问题有着很好的优势[25]，近年来众多学者基于分数阶流变理论建立了煤体循环加卸载条件下的岩石本构模型[26-29]。王军保等[26]引入分数阶Burgers 模型用于研究盐岩在低频循环荷载作用下的蠕变特性，并通过实验结果对所建轴向蠕变的Burgers 模型进行了验证。郭建强等[27]提出弹性、黏性和塑性疲劳元件组合建立了循环荷载作用下岩石的一维非线性黏弹塑性疲劳模型。PU 等[28]将循环荷载分解为静荷载和平均应力为零的循环荷载，并根据流变力学理论和黏弹性力学理论分别建立对应的流变本构方程，再将两者进行叠加得到新的岩石循环加卸载动态本构方程。

综上，分数阶在描述煤体在循环加卸载条件下的力学特性具有一定的优势，并且温度对于深部煤体力学特性的影响不可忽略。鉴于此，笔者拟开展不同温度下的循环加卸载试验，通过煤体的全应力应变曲线，结合能量原理分析煤体在试验过程中弹性参数、能量演化及损伤特性。利用连续介质损伤理论引入热力耦合损伤变量，建立了基于分数阶导数的深部煤岩考虑温度效应的黏弹塑性本构模型，并通过试验结果验证了其适用性。

1 试验流程

1.1 试样制备

试验所用煤样取自中国平煤十二矿己16-17-17200工作面，取样深度约为1 050 m。依据国际岩石力学学会(ISRM)建议方法将采集的煤体经过钻孔取心、切割、打磨等工序，加工成ϕ50 mm (直径)×100 mm(高度)的圆柱体标准试件，保证试样的上下端面直径偏差不大于0.3 mm 以满足试验要求。如图1 所示。

图1 煤样制备示意Fig.1 Diagram of coal samples preparation

1.2 试验设备及方案

本次试验采用MTS 815 Flex Text GT 岩石力学试验设备开展不同温度下的循环加卸载试验，围压设定为25 MPa，参考平顶山矿区平均地温梯度值在32.0～46.0 ℃/km[30]，分别设置为30、50、70 ℃，其中30 ℃为地表浅层煤体所处环境温度，50 ℃为试验煤体所处环境温度，70 ℃为该地区煤体极端温度，并可为未来更深处煤炭资源开发利用提供参考。根据煤样常规三轴试验预设30、50、70 ℃条件下的峰值为30、25 和20 MPa。本试验按照以下步骤进行，具体过程如图2 所示，图中σ1为轴向应力，σ2=σ3为环向应力：① 按照20 ℃/h 的升温速率将目标煤体从常温加到预定温度30、50、70 ℃并保持恒定；② 以20 kN/min 的加载速率将围压和轴压逐渐加载至25 MPa，使煤体处于静水压力状态；③ 采用轴向应力控制方式，以20 kN/min 的速率加载至轴向预定值；④ 以0.15 mm/min 的速率卸载至静水压力，完成一次循环试验；⑤ 重复步骤③和④，逐级增加轴向预定值，考虑到煤体的离散性，可适当增加循环次数，直至峰值；⑥峰后按照峰值的80%、60%，继续进行循环加卸载试验，直致其失去承载能力。

图2 循环加卸载示意Fig.2 Schematic diagram of cyclic loading and unloading

2 试验结果与分析

2.1 不同温度下煤体力学特性分析

通过对试验所得数据进行后处理，得到了各个煤体的循环加卸载应力应变曲线如图3 所示。

图3 不同温度下煤体循环加卸载的应力-应变曲线Fig.3 Stress-strain curves of cyclic loading and unloading of coal and rock at different temperatures

为了便于观察分析，将原图中前4 个循环部分进行放大并置于原图右方，如图3 中的(ii)所示。从图3可以看出，轴向加卸载应力-应变曲线外包络线与常规三轴压缩曲线类似，在整个加卸载过程中，再加载曲线仍沿着单调加载曲线上升，即煤体具有记忆效应。根据已有研究成果[8,21]，可将轴向加卸载应力-应变曲线分为线弹性阶段、稳定破裂发展阶段、非稳定破裂发展阶段以及应变软化阶段。在30 ℃与50 ℃条件下，第4 次循环轴向应力峰值为实际峰值强度的43.55%和55.72%，可认为煤体仍处于线弹性阶段；第5 次循环过程中试验煤样内部开始出现塑性区，应力应变曲线下弯，煤体内部的原生微裂纹的端部、微缺陷等引起局部应力集中，使得裂隙面产生剪切运动，从而导致微裂纹稳定扩展。而70 ℃条件下，从第2个循环开始出现塑性区。通过比较不同温度下的应力-应变外包络线可以发现，随着温度的升高，曲线会更快出现下弯，说明温度的升高增强了煤体的塑性。30 ℃条件下煤体在第5 次循环过程中失稳破坏，而50 ℃和70 ℃条件下的煤体的失稳破坏出现在第6次循环过程中。在超过峰值强度后，煤体内部出现破裂，承载能力下降，但由于围压的存在，煤体仍具有一定的强度。不同温度下峰后第1 次循环应力峰值与其峰值强度相比较，可知30 ℃时下降幅度最大，为45.88%，50 ℃次之，为29.27%，70 ℃时下降幅度最小，为6.77%。说明温度的升高降低了煤体的脆性。

在循环加卸载过程中，受载煤体内部原生裂隙以及内部微缺陷的闭合、扩展、贯通乃至破坏等都会消耗能量且不可逆。在三轴循环加卸载偏应力-应变曲线中表现为卸载曲线始终位于加载曲线下方，并且下一个循环对应的再加载曲线与卸载曲线相交形成一个闭合的环形区域，此为滞回环。在试验过程中煤体的损伤程度越大，相应的耗散能也就越多，滞回环的面积也就越大。煤体的损伤程度与滞回环的面积呈现出明显的正相关关系，因此可以通过观察计算应力应变曲线中滞回环的面积来表征煤体的损伤程度。

从图3 可以明显看出，在前4 次循环加卸载中，滞回环面积随着循环次数的增加而逐渐增大。30 ℃条件下在第5 次循环加卸载过程中试件到达峰值强度，滞回环面积到达最大，说明在这个阶段试件达到储能极限，前期积聚的弹性能快速释放并转化为煤样破裂形成宏观裂缝所需的表面能，同时伴随着剧烈的声响，说明部分转化为声发射等辐射能，此阶段耗散能达到最大，随后在残余应变阶段逐渐减小。通过对比不同温度下试件达到的峰值强度，30 ℃对应的试件峰值强度为48.61 MPa、50 ℃对应的试件峰值强度为39.19 MPa、70 ℃对应的试件峰值强度为19.64 MPa，可以看出，随着温度的增加，峰值强度逐渐减小，说明温度的升高会加速煤体孔裂隙的扩展贯通，导致损伤的积累加速，进而造成煤体承载能力的降低。

2.2 不同温度下煤体变形特性分析

通过对比不同温度下煤体的轴向应变可知，随着温度的增加，煤体的轴向应变最大值越来越大，分别为2.70%、3.72%、4.38%，表明随着温度的升高，煤体延性逐渐增强；通过观察图3 中不同温度下的环向应变和体积应变随着轴向应变的变化情况，可知不同温度下前4 个循环加卸载过程中，体积应变与偏应力呈线性关系，环向应变稳定增长。在第5 次循环的加载阶段以及之后的循环过程中，煤体的横向体积膨胀抵消了部分轴向体积压缩，然后达到扩容点，此时体积应变增量近乎为0，煤体开始进入屈服阶段，之后体积膨胀大于体积压缩，煤体内部出现宏观裂隙，体积明显增大，最终导致煤体破坏。

变形模量以及泊松比在循环加卸载条件下的演化能直观地表征煤体的变形特性。由于卸载曲线平均斜率与加载曲线直线段斜率相同[21]，因此试验所得的变形模量采用卸载割线斜率(卸载曲线的斜率)。泊松比为每阶段卸载点到静水压力点对应的环向应变之差与此阶段两点对应的轴向应变之差的除数的相反数。试样每阶段对应的变形模量及泊松比如图4 所示。

图4 不同温度煤体变形模量、泊松比与循环次数关系Fig.4 Relationship between deformation modulus and Poisson’s ratio of coal and rock at different temperatures with the number of cycles

从图4 可以看出，不同温度下的试样的变形模量随着循环次数的增加整体呈现线性减小的趋势。30 ℃时变形模量谷值相对初始弹性模量下降了24%；50、70 ℃时试样相比下降了11%、21%；值得注意的是在50 ℃最后2 个循环对应的变形模量略有上升，此阶段位于煤体破坏形成宏观破裂面之后，此时煤体的残余强度主要来自于接触面之间的滑动摩擦力提供。此阶段变形模量的上升主要是由于形成了具有较稳定接触面所致。

不同温度下煤体的泊松比变化趋势大致相同。在线弹性阶段，煤样泊松比增长较为平缓，说明随着循环次数的增加，试样的环向变形速率大于轴向变形速率；当煤样加载到峰值强度后，试样的泊松比急剧上升，这是由于试样在此阶段已经破环，出现了宏观裂缝，从而体积膨胀，由于环向应力相对轴向应力较小，煤体环向开始膨胀，宏观上表现为煤样侧向变形量急剧增大。峰后阶段煤体泊松比开始下降。总体来看，煤体泊松比在到达峰值应力前，呈现出减速上升趋势，在峰值应力处达到最大值，在峰后阶段略有下降。

2.3 不同温度下煤体能量演化

在循环加卸载试验中，煤样从初始状态加载到围压稳定的静水压力状态，煤体内部产生了一定量的能量变化，但这部分能量变化较小，笔者忽略这部分能量。

考虑煤样单元在外力作用下产生变形，假设在循环加卸载试验过程中处于一个封闭系统，与外界没有发生能量交换，由能量守恒定律可知

式中，U为输入能密度(外力对煤体所做的功)；Ue为弹性能密度；Ud为耗散能密度。

从轴向方向看，围压在加载过程中做负功，在卸载过程中做正功；在横向方向上，围压在加载过程中做正功，在卸载过程做负功。在此次循环加卸载试验中，围压水平不高并且固定不变，使得在轴向和横向所做的功绝对值相差较小，笔者予以忽略。

对于试验中第i个循环来说，其输入能密度、弹性能密度和耗散能密度可根据偏应力-轴向应变曲线进行积分得到，如图5 及式(2) ～ (4)所示。其中A为峰值强度点，B为最大轴向应变点对应的应力点，需要注意的是由于“滞回效应”，当到达峰值A后的AB段，仪器仍然在对煤体做正功，故B亦为卸载点。C为卸载到静水压力状态后的轴向应变，D为卸载点B对应的轴向应变。

图5 循环加卸载下的能量计算示意Fig.5 Schematic diagram of energy calculation under cyclic loading and unloading

式中，σi为第i次循环中应力-应变曲线中任意一点对应的偏应力；εi为第i次循环中σi对应的轴向应变。

通过上述方法计算出不同温度下每个阶段煤体内部的能量状态，结果见表1，将表1 中的数据进行具象化，可直观地展示出变化趋势，如图6 所示。首先观察随着循环次数增加各能量的变化趋势。在达到峰值强度前，随着循环应力的增加，输入能密度稳定增加。在前两次循环中，输入能密度以及耗散能密度增长平缓，弹性能密度高于耗散能密度，说明煤体处于弹性阶段，内部裂隙尚未扩展。进入第2、3 次循环后，弹性能密度与耗散能密度都有所增加，但耗散能密度增长的速度比弹性能密度增长速度更快，说明此阶段煤体已经进入塑性阶段，并且内部裂隙已经开始扩展延伸。到了第5 次循环，荷载达到煤体的峰值强度并导致煤体迅速破坏，输入能密度和耗散能密度都急剧增加，弹性能密度增加较为平缓，此为破坏阶段。此阶段煤体内部裂隙扩展并贯通，输入能密度的持续增加使得煤体达到破坏临界值，然后所储存的弹性能突然释放，导致煤体破坏并出现明显的破裂面，加上结构面之间的相互摩擦消耗了大量的能量，这也是耗散能密度急剧增加的原因。在破坏后的2 次循环中，由于围压的存在，煤体仍然具有一定的承载能力，输入能密度以及耗散能密度急剧降低，峰后阶段的耗散能多为结构面相互摩擦而消耗的能量，即摩擦耗散能。

表1 各能量密度随循环次数和温度的变化情况Table 1 Variation of each energy density with number of cycles and temperature

图6 能量随循环次数与温度的变化曲面Fig.6 Energy variation surface with the number of cycles and temperature

其次观察各能量密度随着温度的增加的变化趋势。总体来看，不同温度下能量的演化趋势相似，但大小却不同。弹性能密度可表征煤体在当前条件下释放能量的程度。在前4 次循环中，弹性能密度随着温度的增加而递减，说明温度对煤体的弹性储存能力具有劣化作用。在达到峰值强度以及峰后阶段，弹性能密度与温度变化未表现出相关性，主要是因为不同温度下煤体产生损伤破裂的临界条件不同，并且破坏过程中产生的裂隙，不同结构面之间的摩擦等导致耗散能的增加具有较强的随机性，此时煤体的结构对弹性能密度影响占主导地位。也正是由于不同煤体结构的差异性，使得表现出温度的变化对输入能密度与耗散能密度的影响并不明显。

从上文分析可知通过不同温度下各能量密度的大小对比并不能准确地判断煤体的损伤破裂程度，因此引入无量纲参数[21]：弹性能比例Rie(弹性能在输入能中的占比)和耗散能比例Rdi(耗散能在输入能中的占比)来分析不同阶段弹性能密度和耗散能密度的大小分配规律。

图7 为不同温度下弹性能及耗散能演化规律。在弹性阶段，弹性能占据输入能量的80%以上，耗散能显著低于弹性能。在第3、4 循环加卸载中，30 ℃对应的弹性能比例高于耗散能比例，但50 ℃以及70 ℃对应的耗散能就已经开始高于弹性能，耗散能比例达到了60%左右，说明随着温度的升高，煤体内部热应力对裂隙的增长有促进作用。进入第5 次循环后，煤体达到峰值强度，煤体内部裂隙急剧扩展形成宏观裂缝并伴随着较大的声响，释放并消耗了大量的能量，在图7 中表现为耗散能比例急剧增加，弹性能比例大幅衰减，在30 和50 ℃条件下煤体耗散能比例达到了峰值，分别为0.87 和0.73。而此时70 ℃条件下的耗散能比例并未到达峰值，由于试验加载过程中采用的是轴向应力控制方式，峰值判定条件为轴向应力降低到预设峰值强度的2%时，即认为达到峰值，随后进行卸载过程。由于煤样的孔裂隙的随机分布，轴向应力的增加导致部分裂隙产生一定的错动，导致轴向应力的小幅降低，并达到了系统判定条件，而此时煤体仍未形成宏观裂缝，仍具有承载能力，使得70 ℃条件下其耗散能密度在第6 次循环后到达峰值，为0.83。在峰后阶段，由于循环峰值应力的减小，输入能密度也开始较少，弹性能密度以及耗散能密度也相对的开始减小，耗散能比例开始下降或趋于稳定，弹性能比例相对增长，但耗散能比例仍然大于弹性能比例，说明在围压的作用下，煤体在破坏后又形成了新的相对稳定结构。耗散能比例在第7 次循环时达到最低，30、50 和70 ℃条件下的耗散能比例分别为0.56、0.58 和0.81，相较峰值时下降了35.6%、20.5%和0.02%。随着温度的升高，耗散能比例减小的速率减缓，相对地弹性能比例增速变缓。这是因为当温度增加时，能量的输入相对减小，同时煤体内的损伤和变形受到的限制减弱，因此使得耗散能比例减速下降。

图7 不同温度下弹性能及耗散能演化Fig.7 Evolution of elastic energy and dissipative energy at different temperatures

3 考虑温度影响的分数阶黏弹塑性模型的建立与验证

3.1 考虑温度影响的分数阶黏弹塑性模型的建立

从以上对煤体循环加卸载过程中的能量演化可知，煤体在变形的全过程中损伤是存在并不断累计的。基于连续介质力学理论，单独使用塑性模型或者损伤模型能表征煤体在单调加载中力学演化过程，但并不能准确地表征卸载过程中的强度演化趋势[29]。因此需要考虑塑性-损伤的耦合作用，将虎克体中弹性模量的损伤和塑性应变相结合。并且试验过程中外部荷载具有一定的速率，与时间相关联，而分数阶微分算子具有长期的历史依赖性或记忆效应[28]，被认为是描述黏弹和黏塑性行为的有力工具，特别是用于建立时间相关本构模型，故引入分数阶Abel 元件，得到黏弹性体的一般本构方程[31]为

由前一节的实验结果可知，高温作用产生的热应力会破坏煤体内部结构，导致岩石力学参数的劣化。因此有必要建立考虑温度影响的本构关系。而其中最关键部分是建立合适的损伤变量去表征温度对煤体力学参数的劣化程度。考虑温度对弹模的影响较为显著，根据宏观损伤力学理论[39]，采用间接描述的方法对热损伤DT(0 ≤DT＜1)进行表征，即

式中，ET为材料在T℃下弹性模量；E0为常温下的无损煤体的弹性模量，可通过实验获取。

热力耦合作用下，循环加卸载中煤体的损伤累积在Abel 元件中主要为黏度系数的劣化[33]，在这个过程中温度与外部应力共同作用影响黏性系数。因此我们引入了考虑热-力耦合作用的损伤变量来定义黏度系数的衰减过程。结合式(9)、(10)，可得黏性系数的本构方程表达形式为

其中，DTM为热-力耦合损伤变量， 0 ＜DTM＜1。将式(11)代入到式(7)中，即可得考虑温度影响的变系数分数阶的本构方程可定义为

3.2 循环加卸载条件下的解析解

在三轴循环加卸载试验中，煤体在加载阶段受恒定应力速率控制，在卸载阶段受恒定应变速率控制，如图2 所示。由于2 者都与时间相关联，可将2 者转化为时间-应力-应变问题。在加载阶段，恒定应力速率c1与施加的应力的关系可表示为

同理，在卸载阶段，恒定应变速率c2与煤岩应变的关系可表示为

3.3 分数阶本构模型验证分析

为了验证上述考虑温度影响的分数阶黏弹塑性本构模型的有效性，选取深部煤体在不同温度下循环加卸载的试验数据对模型进行验证，如图8 所示。在加载阶段，使用式(17)对恒应力速率试验结果进行拟合，而卸载阶段需使用式(18)对恒应变速率试验结果进行拟合。模型中的参数通过拟合试验数据获取，具体数值见表2。

表2 不同温度下的试验曲线拟合结果Table 2 Fitting results of test curves in different temperature

图8 不同温度下的试验数据和理论曲线对比Fig.8 Comparison of experimental data and fitting curves at different temperatures

不同温度下的拟合结果如图8 所示，其中图8(a)～(c)中的(i)为整个循环过程中试验数据与拟合曲线的对比。便于观察分析，将前4 次循环部分放大并置于原图右方，如图8(a)～(c)中的(ii)所示。从图8 可以看出总体上本构模型拟合曲线与试验结果吻合良好，相关度较高，表明上述考虑了温度影响的黏弹塑性本构模型能较好地表征深部煤体在循环加卸载条件下的变形过程，特别是能够反映出处于弹塑性阶段煤体的变形情况。

4 结 论

(1) 煤体变形模量随着循环次数的增加呈线性减小的趋势，而煤体泊松比减速上升，在峰值应力处达到极值，然后降低，表明循环加卸载作用对煤体的力学性能有劣化作用。随着温度的升高，损伤累积变缓，对应的应变逐渐增大，表明温度对煤体的延性有增强作用。相同应变条件下煤体的损伤程度加剧，温度的上升加速了煤体损伤的发展，降低了岩石的平均强度。

(2) 将不同加卸载速率相关的应力-应变问题转化为时间-应力-应变问题，同时引入分数阶微积分理论，建立了考虑热力耦合作用的分数阶本构模型，可以有效表征深部煤体的黏弹塑性。

(3) 针对应力控制和应变控制2 种不同基本加载方式，得到了相应模型的解析解，并通过三轴循环加卸载实验数据验证了所提出的分数阶黏弹塑性模型。结果显示考虑温度影响的分数阶黏弹塑性本构模型能很好地表征深部煤体在不同温度下的应力应变关系。