煤与瓦斯突出多物理场分布特征的数值模拟研究

曹 偈，赵旭生，刘延保

(1.瓦斯灾害监控与应急技术国家重点实验室，重庆 400037； 2.中煤科工集团重庆研究院有限公司，重庆 400037)

2019年我国煤炭消费量占能源消费总量的57.7%[1]，煤炭在未来相当长一段时间内还将作为主要能源存在。我国是世界上煤与瓦斯突出灾害最为严重的国家，近年来在大量关停突出煤矿,以及煤与瓦斯突出防治技术不断创新的背景下，煤与瓦斯突出事故呈现下降趋势。但由于煤矿井下条件的复杂性及煤与瓦斯突出防治技术的局限性，突出事故仍是影响煤矿安全生产的主要灾害之一。

目前，学者们普遍认可煤与瓦斯突出是地应力、瓦斯、煤岩物理力学性质综合作用的结果，并利用自主研制的大型煤与瓦斯突出模拟装置开展物理相似模拟试验[2-8]，对突出机理进行进一步的量化分析，为灾害的预防控制提供了基础数据。但由于开展大型突出试验的难度较大且测试参数有限，在试验基础上利用数值分析方法对突出过程煤与瓦斯等参数的变化规律进行研究具有重要意义。徐涛[9]、王路军[10]、颜爱华[11]、刘超[12]等科技工作者利用RFPA数值分析软件研究了突出过程应力场、渗流场等的变化规律；刘雪琴等[13]构建了突出过程含瓦斯煤压力场控制方程，以阳煤五矿赵家分区煤矿为例，利用Fluent软件分析了掘进过程中瓦斯压力、渗透率演变规律；崔聪[14]利用Abaqus软件构建三维有限元模型，分析了不同水平应力条件下突出洞穴形体的孕育特征，并应用平煤八矿突出现场洞穴特征与数值模拟结果进行分析比较，证实了该方法的准确性；刘永立[15]、王锐[16]等利用颗粒流理论对煤与瓦斯突出的临界状态、突出颗粒运动轨迹等进行了模拟分析；蓝航等[17]建立了煤岩动力灾害统一能量方程，并利用FLAC模拟软件计算分析了不同初始瓦斯压力下围岩能量及瓦斯涌出量变化，认为瓦斯的不同参与程度导致了不同的动力灾害显现特征；AN F H等[18]建立了突出发生时考虑煤损伤的瓦斯膨胀能计算方法，利用COMSOL软件模拟分析了不同气压、渗透率及加载条件下的突出瓦斯膨胀能。

综上所述，目前学者们针对其研究重点，采用不同数值分析软件对煤与瓦斯突出过程的各物理场分布规律进行了初步探索。笔者在前期开展的煤与瓦斯突出物理模拟试验基础上，利用COMSOL模拟软件对含瓦斯煤的气固耦合方程进行解算，分析突出后煤体应力—变形—瓦斯场的分布规律，为突出机理的量化分析提供基础数据。

1 煤与瓦斯突出数值模拟

1.1 含瓦斯煤气固耦合方程

煤与瓦斯突出是一个复杂的力学作用过程[1]，在开展数值分析过程中提出基本假设及简化：①煤体为各向同性的弹塑性多孔介质；②不考虑煤岩变形破坏及气体解吸渗流过程中的温度因素；③煤与瓦斯突出后突出孔洞壁后方煤体变形为小变形；④气体渗流服从广义幂定律。

1)应力场控制方程

由前人研究可知，考虑气体压力及吸附膨胀变形影响的煤体平衡微分方程[19]如下：

(1)

式中:G为切变模量；ui为煤体位移；ν为泊松比；α为有效应力系数；p为气体压力；δij为Kronecker参数；E为煤骨架弹性模量；εs为气体吸附引起的骨架变形；fi为体积力。

本研究中采用Drucker-Prager匹配Mohr-Coulomb准则来表征煤体的塑性失稳，其关系式为：

(2)

式中：I1、I2为应力张量第一、第二不变量；J2为偏应力张量第二不变量；C为黏聚力；φ为内摩擦角；σi为主应力分量，i=1，2，3。

2)气体流动场控制方程

煤中气体流动场控制方程由气体状态方程、气体含量方程、气体运动方程及连续性方程构成。

①气体状态方程

(3)

式中：ρg为气体密度；ρn为标准状态下的气体密度；pn为标准状态下的气体压力；z、zn为气体压缩因子。

②气体含量方程

煤中甲烷主要以游离态与吸附态存在，其中吸附瓦斯含量占90%以上。瓦斯吸附解吸满足朗格缪尔方程，则煤中气体含量计算公式为：

(4)

式中：M为煤体中的气体含量；Ma为吸附气体含量；Mb为游离气体含量；a、b为吸附常数；η为煤质校正系数；q为孔隙率。

③气体运动方程

煤中由于裂隙发育的离散性导致气体流动形态不同，可根据当地雷诺数决定，其运动可用广义幂定律描述：

(5)

式中:v为气体渗流速度；k为渗透率；μ为气体动力黏度；m为气体状态指数，m=1～2，当m=1时，式(5)为达西渗流模型。

由于含瓦斯煤岩在应力环境中，应力的改变导致孔隙率及渗透率呈现动态变化。孔隙率计算公式如下：

(6)

(7)

式中k0为初始渗透率。

④气体运动的连续性方程

(8)

将公式(3)～(8)组合变换成偏微分形式，利用Comsol Multiphysics中的固体力学和PDE自定义模块，对方程进行数值求解，模拟突出煤层损伤与气体渗流的耦合规律，并与前期开展的物理模拟试验结果相互验证。

1.2 几何模型及模型参数

研究团队前期利用自主研制的大型煤与瓦斯突出模拟试验装置开展了垂直应力16 MPa、气体压力0.5 MPa条件下的突出模拟试验[20]。试验模型尺寸为1 500 mm×800 mm×800 mm，上覆200 mm厚的刚性压板。突出煤粉15.86 kg，形成深10 cm、宽30 cm左右的突出孔洞。基于此次突出物理模拟试验构建了二维数值分析模型，模型煤层长1 500 mm、高800 mm，上覆刚性压板长1 490 mm、高200 mm，将突出口孔洞简化为半径分别为150、100 mm的半椭圆形。突出几何模型及网络划分如图1所示。

图1 突出几何模型及网格划分

根据物理试验条件，几何模型上部施加16 MPa垂直应力，左右两边及底部采用辊支承边界，突出后突出孔洞内(图1红线所示)为自由边界。突出前，煤层内部初始气压p0=0.5 MPa，煤层四周边界均为零通量不透气边界(dp/dxi=0)；突出后，突出孔洞边界为大气压。数值模拟以突出物理模拟试验模型中取样测试得到的煤层物性参数[20]为基础，如表1所示。模拟分为2个步骤，首先利用固体力学模块计算得到突出前煤层应力状态；然后利用固体力学与PDE模块计算突出后煤层应力、变形，以及气体流动规律。

表1 模型物性参数值

2 模拟结果分析

2.1 数值模拟结果验证

前期开展的物理模拟试验为了采集突出过程轴向应力变化，在煤层垂直高度600 mm的水平面中心线上距突出口100、350、550、750 mm位置布置有应力传感器。由试验结果可知，在距突出口100 mm处的煤层应力由于下方煤体破碎并抛出，形成突出孔洞，其应力呈下降趋势。而孔洞壁后方煤体由于应力重新分布，垂直应力有所升高，且越靠近孔洞壁处的应力变化越大。为验证数值分析结果的可靠性，在模型相同位置设置观测线。数值模拟与物理试验的突出前后煤层应力分布对比如图2所示。

图2 突出前后煤层应力分布的对比图

由图2可以看出，由于数值模型的简化，数值模拟结果与试验结果有一定差异。但数值模拟结果与物理试验结果具有一致的变化规律，突出发生后在距突出口100～350 mm区域应力呈现降低趋势，距突出口350～1 150 mm区域应力呈增大趋势。证实了数值模型可用于对突出前后煤层物理场分布规律进行分析。

2.2 突出后各物理因素的分布规律

1)应力变形规律

突出后煤层应力、体积应变、塑性区分布云图如图3所示。突出发生前由于轴向应力施加在刚性压板中心700 mm范围内，形成了煤体上表面实际施加了非均匀分布力的特点，模拟了实际工作面煤体应力集中现象。

(a)垂直应力

(b)体积应变

(c)塑性区

由图3可以得出，突出后煤体应力达到新的平衡状态。孔洞壁上下面由于突出过程煤体被抛出而卸压，应力降低；孔洞壁后方形成应力增高区，此处煤体体积应变增大，产生塑性变形，破坏较为严重。而到煤体深部受突出影响较小，煤体体积应变并不显著。

2)气体流动规律

突出后煤体渗透率及气体压力分布云图如图4所示。

(a)渗透率

(b)气体压力

由图4(a)可知，由于突出孔洞壁后方煤体破坏严重，其渗透率在此处达到最高，从孔洞壁面到煤体深部呈减小趋势，在体积应变最大值处达到最小，随后恢复至原始渗透率。

由图4(b)可知，煤体内气体压力从孔洞壁面到煤体深部急剧上升，在应力集中处达到最大值，且高于煤体原始气压，随后恢复至原始气压，如图5所示。这是由于应力集中使煤体被进一步压实，渗透性降低，气体难以向外流动，使瓦斯在短时间内迅速积聚，在突出孔洞壁面附近产生较高的压力梯度。若此时突出壁面附近煤体破碎严重，且产生的气体压力梯度具有较高能量将煤体抛出，突出则会持续发生。

图5 随着时间推移突出口中轴线上气体压力分布曲线

由图5还可以看出，随着时间推移，突出孔洞壁面后方气体压力逐渐降低，气体压力梯度减小，5 s时间内气压峰值消失，说明突出的持续发生具有时效性。同一位置处，气体压力随时间的衰减速率也呈现降低趋势。

不同垂直应力条件下突出口中轴线上煤体渗透率的分布曲线如图6所示。

图6 不同垂直应力条件下的煤体渗透率分布曲线

由图6可以看出，不同应力条件下煤体渗透率受突出影响程度不同。当应力较小时煤体塑性变形不显著，渗透率变化较小；随着应力增大，突出孔洞壁处煤体渗透率变化明显，当应力大于 16 MPa 时，比煤体原始渗透率增加2～3个数量级。由此，对于高应力环境突出，较大的渗透率变化导致煤体气体流动及压力衰减较快，突出发生后在极短的时间内气体压力梯度减小，不具备突出持续发生的条件，突出煤量相对较小。

3 结论

利用Comsol Multiphysics中的固体力学和PDE自定义模块构建了煤与瓦斯突出的数值分析模型，获得了突出后煤体应力场、变形场、气体流动场分布规律，并与前期开展的煤与瓦斯突出物理模拟试验结果进行对比验证，证实了数值分析方法的可行性。分析得到突出后孔洞壁附近煤体发生塑性变形，导致此处渗透率增大，且随着应力增大，突出孔洞壁处煤体渗透率可比原始渗透率增加2～3个数量级。而气体压力从突出孔洞壁面到煤体深部呈现先急剧上升后降低至原始气压的变化趋势，在孔洞壁处形成较高的压力梯度，这是突出持续发生的必要条件。