含瓦斯水合物煤体力学性质离散元建模与试验验证

高 霞， 王婷婷， 王维亮

(黑龙江科技大学 建筑工程学院， 哈尔滨 150022)

0 引 言

随着煤矿开采深度逐年增加，应力与瓦斯动力灾害情况变得更加复杂与严峻[1]，煤层结构破坏导致煤层应力结构发生破坏，从而引起煤与瓦斯突出等情况[2-5]。吴强课题组[6]提出煤与瓦斯水合固化防突技术，该技术利用水合物特性改善了煤体的力学性能，但是试验无法观察到试样内部的变化情况，故采用离散元方法展开研究，解决非连续介质材料的数值模拟问题[7]。

离散元单元法用于研究非连续体的力学性质[8]。郭放[9]利用离散元方法研究了煤样内部破坏演化过程。蒋明镜等[10]采用离散单元法实现了纯水合物宏观强度特性的力学试验模拟。Brugada[11]、Jung[12]等利用离散元随机生成土颗粒并在孔隙中生成球形水合物颗粒，研究填充型模型水合物沉积物的力学性质。王璇等[13]模拟了含水合物沉积物双轴剪切试验，从粘结断裂、剪切带变化两方面分析加载过程中含水合物沉积物微观响应。蒋明镜等[14]利用离散单元法研究裹覆型能源土在复杂应力路径下的力学性质。针对水合物赋存的形式，有学者在介质中填充单颗粒水合物[15-16]或者通过胶结模型代替水合物[17-19]进行力学性质研究。综上所述，学者对水合物沉积物力学性质进行了大量离散元数值研究，但是关于含瓦斯水合物煤体力学性质模拟研究较少。

笔者在试验基础上采用数值模拟对含瓦斯水合物煤体进行力学性质研究，分析应力-应变曲线、颗粒位移以及速度场的变化规律，揭示了细观变化过程与宏观变化之间的联系，为含瓦斯水合物煤体三轴压缩试验的离散元模拟提供建模思路。

1 试 验

试验材料的煤来自东保卫煤矿41#煤层，取样深度是800 m。试验所采用的仪器是本课题组自主研制的设备——瓦斯水合固化三轴试验测试一体化系统，试验步骤及数据来源于文献[20]。为了进一步分析含瓦斯水合物煤体在不同影响因素下的破坏情况，文中利用同一煤粉粒径0.850～0.425 mm(20～40 目)在16 MPa下进行3种饱和度(7.72%、9.69%、12.36%)的三轴压缩数值试验。

2 数值模型的建立

本文对含瓦斯水合物煤体进行离散元数值建模，如果按照水合物在含瓦斯水合物煤体中真实赋存形式，则生成的水合物既包括胶结模式[21]，又包括填充模式[22]，难度较大，且不容易实现。因此本文在建模时进行两方面假设：假设模型中只有水合物和煤两种物质；假设水合物均匀包裹煤颗粒。图1为含瓦斯水合物煤体中物质简化过程。

图1 数值模型简化Fig. 1 Simplification of numerical model

2.1 三轴压缩试验模拟流程

试验试样尺寸是φb×h=50 mm×100 mm，若按照试样尺寸进行1∶1建模，则颗粒数量超过几十万，影响计算效率，故离散元试样尺寸一般小于物理实验的试样尺寸。部分学者用φb×h=2 mm×4 mm进行建模，如杨期君等[23]对试样进行了圆柱形φ2 mm×4 mm的数值模拟，采用离散元方法开展了不同水合物饱和度沉积物的研究；徐小敏等[24]对砂土采用圆柱形φ2 mm×4 mm进行三轴压缩试验模拟，研究砂土的宏细观参数相关性。因此，本文数值试样尺寸设定为φb×h=2 mm×4 mm，共生成2 819个颗粒，颗粒模拟如图2所示。

图2 颗粒模拟Fig. 2 Particle simulation

(1)墙体建立：用 wall create 命令创建一个圆柱形的侧面和两个平行的上下底面，圆柱形侧面为试样提供围压，上下底面作为加载刚性墙(图2a)。

(2)试样生成：利用ball distribute在圆柱形内生成颗粒(图2c)。颗粒之间模型采取接触黏结模型，颗粒与墙之间采取线性接触模型。

(3)三轴加载：待试样施加围压并稳定一段时间后，需对试样进行加载，加载方式为上下墙体同时加载，确保加载时试样为准静态模式，加载速率的影响就可以忽略[25]。在本数值模拟中，采用的加载速率为0.1 mm/s。如果按照含瓦斯水合物煤体颗粒实际的目数生成颗粒，生成的颗粒数量过于庞大，且运行速度较慢，故依据模型整体尺寸的1/30大于颗粒平均粒径时[26]，颗粒的尺寸效应可以忽略，因此本文将颗粒粒径设为0.075～0.100 mm，随后进行三轴伺服及加载，从不同细观角度分析试样力学性质变化情况。试样生成流程图如图3所示。

图3 颗粒生成步骤Fig. 3 Step of particle generation

2.2 接触模型选取及细观参数确定

模型选取是根据研究对象材料的力学性质决定的，接触模型示意如图4所示。蒋明镜等[14]在对水合物沉积物进行模拟时通过胶结模型将颗粒进行胶结，该方法可以很好体现水合物性质，因此文中选用接触黏结模型。细观参数标定影响离散元模型的吻合程度，与试验数据相对比，调整模拟数据，利用“试错法”[27]对细观参数进行调整得到细观参数的精确值。细观参数标定过程如图5所示，最终得出本试验细观参数的取值，如表3所示。

图4 颗粒间接触黏结示意Fig. 4 Schematic of intergranular contact bonding

图5 细观参数标定流程Fig. 5 Microscopic parameter calibration process

表3 细观参数取值

3 模拟结果分析

3.1 数值模型验证

图6是试验与模拟的应力-应变曲线对比图，从图中可以明显地观察到弹性阶段、屈服阶段和强化阶段，且曲线均呈应变硬化型。

图6 不同饱和度下试验模拟应力-应变对比Fig. 6 Simulated stress-strain contrast under different saturation

由图6可知，饱和度为7.72%时，试验和模拟的峰值强度分别是23.2和21.4 MPa，误差率为8%；饱和度为9.69%时，试验和模拟的峰值强度分别是25.4和23.7 MPa，误差率为7%；饱和度为12.36%时，试验和模拟的峰值强度分别为27.3和26.3 MPa，误差率为4%。3种饱和度的应力-应变曲线误差率均小于10%[26-28]，标定的细观参数数值合适[29]。从应力-应变曲线可以得出，随着饱和度逐渐增大，含瓦斯水合物煤体的峰值强度逐渐增大，煤体强度得到提高，水合物含量越多，煤体的峰值强度提高越多。数值模拟曲线变化趋势与试验曲线变化趋势相似，随着饱和度增大，误差率逐渐减小，数值模拟的曲线更加贴合试验曲线，峰值强度误差率也越来越小。因此，接触黏结模型能较好地对含瓦斯水合物煤体进行模拟。

3.2 颗粒位移场分析

图7为3种饱和度下颗粒位移的变化情况，箭头所指方向代表颗粒运动的方向，箭头颜色表示颗粒位移量的大小；色条从下到上表示颗粒位移由小到大的变化，蓝色表示最小位移量，红色表示最大位移量。

图7a、b的位移场箭头颜色大多偏蓝，图7c位移场的箭头呈绿色居多，随着饱和度增大，位移场颗粒位移变化量逐渐增大。可以看出，3种饱和度颗粒位移场分布规律相似，且均保持着轴对称的特征；总体来看，颗粒位移矢量随着饱和度增加也相继增加。

3种饱和度位移场在第10 000 时步和第20 000 时步时，颗粒位移分布比较均匀，且颗粒均向外侧运动；从第30 000 时步开始，靠近上、下加载板的颗粒开始向内侧运动，中部颗粒保持不变，依然向外侧运动，颗粒从两端往中间变化的趋势是位移逐渐减小。

在加载过程中，可以看出三种饱和度位移场分布规律相似，均保持着轴对称的特征，且颗粒的相对位置变化也不大；总体来看，饱和度为12.36%的颗粒位移矢量大于饱和度为7.72%和9.69%的颗粒位移矢量，分析认为可能是摩擦系数过小所造成的。在数值试验结束后，试样内部的颗粒仍然维持这种运动形式，且靠近上加载板和下加载板的颗粒位移明显变大。含瓦斯水合物煤体颗粒在加载过程中颗粒位移场的变化情况可以解释煤样发生了压胀破坏。

图7 不同饱和度下不同时步的颗粒位移场变化情况Fig. 7 Variation of particle displacement field at different time steps under different saturation

3.3 速度场分析

图8为不同饱和度(7.72%、9.69%、12.36%)的含瓦斯水合物煤在加载过程中不同时步的颗粒速度变化情况，箭头所指的方向表示颗粒速度方向，箭头不同颜色代表颗粒速度量变化大小，色条从下往上表示颗粒速度量由小到大的变化，深蓝色表示最小的速度，红色表示最大的速度。

图8 不同饱和度颗粒速度变化情况Fig. 8 Particle displacement changes at different saturation levels

由图8可知，试样靠近上部的颗粒向下运动，靠近下部的颗粒向上运动，颗粒从两端往中间变化的趋势是速度逐渐减小，速度场总体保持着轴对称的特征。图8a、b、c在第60 000时步时，对应的轴向应变为12%左右，在此时刻可以明显地看到颗粒溢出墙体的现象，在这一时步，可以看到模型的右上方颗粒向下移动，左下方颗粒向左上方移动，中部颗粒向两侧移动。随着饱和度增加，在饱和度为12.36%时，且在第60 000时步可以明显地看到从右上到左下的剪切带，形成一条由左上到右下的剪切带，印证了煤样的破坏状态，如图9所示。

图9 破坏对比图Fig. 9 Damage comparison diagram

4 结 论

(1)利用数值模拟对含瓦斯水合物煤体的建立与分析，发现应力-应变曲线模拟结果与试验结果相似，峰值强度误差率较小，应力-应变呈现应力硬化型。

(2)含瓦斯水合物煤体在加载各个阶段的颗粒位移场呈轴对称向外分布，说明煤体发生了压胀破坏，与试验煤体的压胀破坏的情况相似。

(3)通过对加载过程颗粒速度矢量变化情况分析，发现在加载过程后期，模型形成一条由左上到右下的剪切带，印证了煤样发生了破坏。