瓦斯对冲击倾向性能量指标影响的数值模拟

李少刚，张广辉，邓志刚

（1.煤炭科学技术研究院有限公司 安全分院，北京 100013；2.煤炭资源高效开采与洁净利用国家重点实验室（煤炭科学研究总院），北京 100013）

冲击倾向性是煤岩突然破坏的固有属性，是冲击地压发生的内在因素[1]。针对煤层冲击倾向性，国内外学者从冲击倾向性指标的关联性[2]、组合煤岩样[3]、煤样含水率[4]、试件介质成分[5]、应变速率[6]和破坏机制[7]等方面开展了大量的研究工作，但均没有将瓦斯作为冲击倾向性的影响因素开展研究，同时伴随开采深度的增加煤层瓦斯压力大、瓦斯含量高，若深部含瓦斯煤层发生冲击地压，必将引发更为严重的后果[8]。《煤层冲击倾向性分类及指数的测定方法》指出冲击倾向性指标需在单轴压缩下测定[9]，而现阶段瓦斯对煤力学性质的影响多是在假三轴下进行，在围压作用下煤样的物理力学参数势必有所改变。因此，为探究瓦斯对煤冲击倾向性能量指标的影响，建立含瓦斯煤样气-固耦合控制方程，并利用COMSOL多场耦合数值模拟软件，通过研究瓦斯对煤样中塑性变形区的作用，开展瓦斯对冲击能量指数和弹性能量指数影响研究。

1 含瓦斯煤气-固耦合控制方程

含瓦斯煤样的变形破坏是外部施加的应力与内部瓦斯共同作用的结果，其涉及到固体力学、表面物理化学和流体力学等众多学科，在建模过程中，做出如下假设[10]：①含瓦斯煤是各向同性弹塑性材料；②煤层瓦斯含量符合朗格缪尔方程；③瓦斯在煤样中的流动负荷Darcy定律；④含瓦斯煤发生的弹塑性变形为小变形且具有连续性；⑤孔隙系统内游离态气体和吸附气体温度相同且恒定。

在以上假设的基础上建立方程。

1.1 应力平衡方程

建立煤样模型中假设含瓦斯煤为理想线弹性体，同时考虑温度场的改变会致使多孔介质体积应变的变化，瓦斯在煤层中的吸附与解吸将导致其骨架的膨胀与收缩，作用在煤样上的外部载荷会引发其正应变，孔隙压力作用与多孔介质上也会引起其骨架的膨胀。考虑外部应力、孔隙压力、瓦斯吸附及温度场的广义胡克定律可表示为[11]：

式中：εij为应变张量；G 为剪切模量，GPa；σij、σkk为应力张量；K为体积模量，GPa；α为孔隙压裂系数；δij为Kronecker符号；εs为瓦斯吸附引起的体积应变；αs为煤的体积热膨胀系数；T为温度变化量，℃；p为瓦斯压力。

在模拟中为提升计算速度，简化运算流程，假设外部温度场恒定且不考虑瓦斯吸附和解吸对煤样骨架的影响，式（1）可改写为[11]：

式中：v为煤的泊松比；εkk为应变张量。

式（2）中将瓦斯压力置于应力应变方程中，实现了孔隙压力对煤样骨架变形破坏的影响。

1.2 孔隙率和渗透率动态演化方程

在模拟单轴压缩过程时，煤样受到外部变化载荷和孔隙压力的共同作用，其变形破坏处于动态变化，进而导致孔隙率和渗透率变化。J.P.Davies和D.K.Davies[12]根据实验室内实验得出了孔隙率与平均有效应力的关系，即：

式中：φ为煤样中动态孔隙率；φ0为煤样初始孔隙率；φr为煤样残余孔隙率；αφ为渗透敏感系数；σv为平均有效应力，MPa。

根据文献[12]，多孔介质中渗透率与孔隙率满足以下关系：

式中：k∞为煤样中动态的渗透率，m2；k∞0为煤样初始渗透率，m2。

1.3 瓦斯流动和含量方程

假设瓦斯气体在煤样中的流动符合达西定律，其渗流速度可表示为[11]：

式中：g为重力加速度m/s2；kg为考虑 Klinkenberg 效应的煤层渗透率，m2；ρg为气体密度，kg/m3；qg为达西渗流速度，m/s；μg为气体动力黏度，Pa·s；z为高差。

由于瓦斯气体密度与高差的乘积相对于瓦斯压力梯度较小，因此瓦斯的重力效应可忽略。

含瓦斯煤是一种多孔介质，煤样中的瓦斯由游离瓦斯和吸附瓦斯共同组成，游离瓦斯含量与煤样中的孔隙率相关，而吸附瓦斯含量可由Langmuir方程计算，单位体积煤内瓦斯含量可表示为[11]：

式中：m为气体含量，kg/m3；Cf为游离瓦斯质量，kg/m3；Ca为吸附瓦斯质量，kg/m3；VL为朗格缪尔体积参数，m3/kg；pL为朗格缪尔压力参数，MPa-1；ρs为煤的密度，kg/m3。

在等温条件下多孔介质中气体流动满足质量守恒定律[11]：

将式（1）～式（6）代入到式（7）中可得到含瓦斯煤样在变形破坏时本构方程[11]：

式中：εv为体积应变；ks为煤骨架的体积模量；pa为标准大气压力；t为时间。

式（8）中煤样体积应变εv和瓦斯压力p相互影响，实现了气-固双向耦合。

2 含瓦斯煤单轴压缩模型参数及边界条件

在COMSOL软件中选用固体力学模块实现对煤样的应力应变求解，利用Darcy模块完成煤样内瓦斯的流动，并在2个模块间设置变量实现多场耦合。由于煤样为轴对称的圆柱体，在建立模型过程中，选择二维轴对称，在该模式下可提升模型迭代平衡的速度；因此将软件模型宽度设置为真实宽度的1/2，即25 mm，以边长为2 mm的尺寸划分单元格（图1）。在模型下端部使用固定约束，右侧边界无约束，且右侧边界无气体流动，瞬态条件下在上端部施加应力载荷或位移载荷，模型基本参数见表1。

图1 单轴试件尺寸、边界条件简图

表1 模型基本参数

为充分模拟煤样中微孔隙微裂隙非均匀分布的赋存状态，还原煤样各向异性的力学性质，因此在COMSOL软件中定义二元正态函数rnl（x，y），用该函数修饰煤样的黏聚力分布（图2），由随机分布的黏聚力可充分反映煤样内部缺陷的随机分布及煤样的各向异性，在低应力作用下煤样即可出现部分区域的塑形变形。

图2 煤样黏聚力随机分布图

3 模拟结果及讨论

3.1 瓦斯对冲击能量指数影响分析

冲击能量指数KE是指用伺服实验机以0.5×10-5mm/s应变速率对煤样进行加载，冲击能量指数计算示意图如图3。

图3 冲击能量指数计算示意图

煤样峰前积聚总应变能量W1与峰后破坏消耗变形能W2比值即为冲击能量指数，冲击能量指数表征了煤样试件完全破坏后盈余的能量；煤样变形破坏消耗能量越少，KE值越大，完全破坏后剩余能量越多且将以动能和辐射能形式释放。

基于现有冲击倾向性鉴定标准，不考虑围压作用下瓦斯对煤冲击倾向性影响的物理实验难以实现，且实验中由于试件本身的唯一性和制样过程中的机械损伤，必将导致数据结果具有一定的离散性，因此在数值模拟过程中采用煤样塑性应变区域及数值反映煤样能量耗散情况和储存弹性应变能的能力。在4种瓦斯压力下煤样内部塑性应变分布图如图4。

从图4中可以看出，在4种瓦斯压力下煤样的塑形应变均呈现出由右上端角和右下端角向煤样中心处发展的趋势，由于模拟中采用轴对称模型，所以应变云图中塑形应变的走势对应“X”型破坏的右半部分，与实验室内煤样实际破坏形式和国内外其他学者的实验结果也相符。其次伴随瓦斯压力的增加煤样塑性变形区域逐渐扩大，塑性变形区域内塑性变形程度也呈现加剧的趋势。模拟结果表征伴随瓦斯压力的增加煤样内部由损伤和塑性变形引起的耗散能增加，内部储存的弹性应变能减少，煤样完全破坏后盈余能量减少，进而导致煤样冲击能量指数降低，单项冲倾指标弱化。

图4 4种工况下煤样塑性应变分布图

3.2 应力加载时瓦斯压力对塑形变形的影响

弹性能量指数WET是多级循环加载下，极限强度前可释放的弹性能与损耗能之比，其表征了煤样在破坏前的储蓄弹性能的能力。冲击能量指数计算值为：

式中：Ue为储存于煤样中可释放的弹性应变能；Ud为煤样永久变形引起的损耗能。

多级循环加载中每级加卸载曲线负荷-变形示意图如图5。

在本次模拟中通过在COMSOL中定义分段函数实现了对煤样单次循环加载，并提取了初次卸载后煤样内部的塑形应变云图（图6）。

由图6可知，WET模拟中煤样塑性应变图的特征与冲击能量指数模拟中应变图特征类似，伴随瓦斯压力的增加，相同煤样中塑性变形区域有所扩展，且塑性应变有所加剧，在循环加载中孔隙瓦斯压力促进了煤样内部孔隙和裂隙的发育，消耗了部分弹性应变能，导致弹性能量指数降低，而在实验室测定中则表现为负荷-变形曲线的滞回环越来越大，煤样弹性能量指数减小。由于外部载荷并未达到煤样破坏所需应力，因此试样内部塑性变形范围和数值相比冲击能量指数时较小。

图5 煤样加卸载曲线负荷-变形示意图

图6 WET模拟中煤样塑性应变分布图

由3.1节和3.2节可看出伴随瓦斯压力的增加煤样内部塑性应变区域逐渐扩展，且塑性应变值不断增长，其原因为瓦斯弱化了煤的强度[13-14]，在受到相同外部载荷时，由损伤和塑性变形引起的能量耗散增多，煤样本身储存弹性应变能的能力下降，煤样破坏形式由脆性向脆塑性转变，在煤样变形破坏过程中消耗能量增多，完全破坏时盈余能量减少，致使煤样弹性能量指数和冲击能量指数降低，导致冲击倾向性由强向弱或由弱向无转变，瓦斯弱化了煤的冲击特性，而在深部高瓦斯煤层开展冲击倾向性鉴定或冲击危险性评价过程中应考虑瓦斯的影响。

4 结 论

1）伴随瓦斯压力的增加煤样内部塑性应变区域扩展，且塑性应变值有一定程度的增加，煤样储存能量减少，导致煤样冲击能量指数和弹性能量指数降低。

2）瓦斯的存在致使煤的强度降低，弱化了煤的冲击倾向性，在含瓦斯煤层冲击倾向性鉴定或冲击危险性评价过程中应充分考虑瓦斯对煤层冲击特性的影响。