基于FLAC3D的煤岩体单轴压缩数值模拟研究

褚佳琪

（安徽理工大学 矿业工程学院，安徽 淮南 232001）

0 引言

近年来，随着我国煤炭资源开采深度和开采强度的增加，发生冲击地压的矿井数量快速增多，事故频度和强度也明显提高。由于我国矿山工程地质灾害问题日渐突出，社会的生产活动以及生态环境均受到严重制约，而矿山工程地质灾害受到人为活动条件、气候植被条件、地形地貌条件、地质条件及历史灾害发生的规模、频次、密度等许多因素的影响，近年来发生了不同规模大小的地质灾害，其中煤岩动力灾害问题尤为突出。煤岩动力灾害的主要表现为地震、火山喷发、山崩滑坡、煤与瓦斯突出、冲击地压等，其中冲击地压是一种与地应力有关的地质灾害。据悉，矿山进入1 000 m 以下进行深部开采时，由于高地应力、高地温、高渗透压以及强烈开采扰动的影响，煤岩体往往发生冲击地压。

冲击地压的发生，与煤岩体的力学特性、组合形式以及受力环境等各种因素密切相关，众多学者针对煤岩组合体的特性进行了研究探讨，付斌等利用数值模拟软件，对不同围压和不同组合倾角条件下的煤岩组合体进行了力学特性和声发射特征研究；胡顺银等对煤岩模型破坏时的应力及变形情况进行了模拟分析；张凯文等模拟分析了真三轴加载过程中岩石的损伤与破裂演化过程及特征；蒋中明等对不同围压下对岩石进行了不同变形参数及强度参数的数值模拟研究；周天白等采用数值模拟反演了煤岩材料的参数，对煤岩材料变形破坏的主控参数进行了研究；贾善坡等研究了层状岩体单轴压缩情况下的变形破坏特征及其影响其破损的规律；周元超等针对不同高度比的煤岩组合体，模拟研究了不同组合方式下煤岩组合体的力学特性和声发射特征。

本文基于FLAC3D 的数值模拟研究，采用摩尔-库伦力学模型，建立煤岩体力学数值计算模型，进行单轴加载数值模拟，对比分析不同加载时步下（step200、step400、step600）的煤岩体应力及位移演化规律，同时对煤岩组合体的应力—应变曲线进行分析，从数值模拟的角度对煤岩单轴加载下的力学特性进行认识。

1 煤岩力学机理

矿山在进行深部开采时，由于高地应力的存在，使得煤岩体在一定的环境条件下往往发生冲击地压。煤岩体的变形是衡量深部矿山稳定性的重要指标，故探讨煤岩体的力学机理有着重要的意义，若煤岩体的变形破坏较大，对深部矿山的危害也会相对较严重。在深部矿山中，煤岩体内任一点的应力状态可由正应力σ、σ、σ和剪应力τ、τ、τ、τ、τ、τ组成的应力分量表示，且可以用应力张量σ表示：

式中，，=，，，且σ可用力学模型表示，如图1所示。

图1 煤岩体力学模型示意图

由图1 可知，煤岩体在深部矿山中处于三向受力的状态，对于煤岩体所受到来自、、方向的正应力与剪应力，符合广义的胡克定律，即煤岩体在受力之后，煤岩体的应力与应变之间成线性相关关系：同时，煤岩体在深部矿山中的应力及应变也可用张量分量表示：

式中，、为拉梅常数；煤岩体的力学本构关系可通过拉梅常数表示：

2 模拟过程及分析

2.1 模型建立

本次数值模拟采用由美国ITASCA 公司开发的仿真计算软件—FLAC3D 软件，该软件采用三维连续体的快速拉格朗日分析，同时其包含的11 种材料力学模型可以更准确地对要研究的对象进行数值模拟及力学分析。本次采用的力学数值模型是基于在进行试验时的实际试样，运用FLAC3D 数值模拟软件的导入命令模式，将编写的模型命令导入FLAC3D。基于煤岩体的受力变形特征，建立一个高度为100 mm（Y 方向），直径为50 mm 的标准岩-煤-岩组合体，其中煤体受到岩体的夹持。在FLAC3D 中，首先利用generate zone 的命令将该煤岩体模型总共划分了1 000 个单元，并生成标准圆柱形的模型，同时保证各个网格节点的连接，使得在后续的模拟过程中保证其受力均匀。该模型与试验使用的标准试样一致，能够满足模型的模拟。在完成了初步的试验模型以后，本文单轴模拟试验采用摩尔-库伦塑性模型，该模型可应用于煤岩体的实际承载能力和失效荷载的计算，以及其他以煤岩体破坏为关键因素的数值计算，其中group1 及group3 为岩体，group2 为煤体。初步模型如图2所示。

图2 初步模型图

对模型底部施加约束条件=0，其余方向为自由边界，对该模型进行单轴压缩数值模拟，并利用FISH 语言自动记录应力—应变的数值大小，煤岩体中各物理力学参数如表1所示。

表1 煤岩物理力学参数

在完成煤岩体数值模型建立和参数赋值以后进行模拟计算，数值模拟加载方式采用沿着轴线方向施加0.5 MPa/s的轴向应力，在其余边界上不施加任何应力，即在轴向上σ=，而在环向及其他方向上==0；这样的加载方式有利于保证与实际的实验加载方式一致，以确保实验的准确性，同时确保数值模拟计算能够顺利进行。

2.2 煤岩组合体应力及位移分析

以200 时步递增的顺序对煤岩体进行单轴压缩试验，并以穿过煤岩体试件轴线的平面对模拟结果进行切片，分别记录step200、step400 以及step600 时的应力与位移演化云图，具体模拟过程如图3、图4所示，其中图3 为应力分布切片，图4为位移分布切片。总体来看，煤岩体试件应力及应变变化呈现以圆柱试件中心轴为对称轴的高度对称变化，煤岩体试件在Y轴方向上的位移上部呈下沉式形变，底部呈现上升式形变，且顶底面形变量最大，向煤岩体试件中部位置依次衰减，随着加载时间的增加，高位移区逐渐向试件中心部分延展。

图3 Y 方向应力演化云图

图4 Y 方向位移演化云图

煤岩体试件应力变化情况大致可分为Ⅰ、Ⅱ、Ⅲ三个阶段：

（1）当加载运行至200 时步时，此时为煤岩组合体单轴加载的初期，即变化的第Ⅰ阶段。受上下部荷载影响，岩体试件在Y 轴方向的应力方向为Y 轴负方向，且岩体中心上部中心部位应力最大，煤体内部由于受岩体应力的作用，呈现与岩体应力方向相反的应力，并在煤体中心区域形成较大的应力集中区，此时煤体内部的应力略大于岩体内部的应力，而煤体内部暂无位移变化，岩体沿着Y 方向开始向煤体进行位移，且上下偏移量的数值相同。

（2）当加载运行至400 时步时，此时为煤岩体组合单轴加载的中期，即变化的第Ⅱ阶段。随着加载时间的增加，岩体试件的高应力区范围逐渐增加，且影响范围逐渐向煤岩组合体的中心处延展。煤体中心处应力转变为与岩体相同的Y 轴负方向，且煤体两侧的应力达到20.47 MPa，煤体中心处应力最大，最大可达到62.3 MPa。

（3）当加载运行至600 时步时，此时为煤岩体组合体单轴压缩的后期，即变化的第Ⅲ阶段，此时煤体上下部岩体的应力集中逐渐消失，加载对岩体产生的影响能够较快速度的进行传递，煤体两侧的应力达到20.32 MPa，相比较于第Ⅱ阶段，煤体两侧的应力逐渐向煤体中间靠拢，煤岩组合体整体中间部分应力最大，此时最大可达到72.9 MPa，岩体的位移变化也逐渐向内部扩展。

煤岩组合体在三个阶段都具有明显的变化特征，即Ⅰ阶段煤体与岩体试件具有应力方向相反的应力集中区域，Ⅱ阶段煤体应力与Ⅰ阶段媒体应力相反，Ⅲ阶段岩体应力集中区域消失，同时，应力呈现出由岩体传递向煤体，由煤体两侧传递向煤体中心的趋势，通过图4 同样可以看出，随着加载时步的增加，煤岩体Y 方向的位移由岩体逐渐向煤体转移，大小呈现递增的趋势。

图4 中国独立动画理论研究突显词视图（左）与控制面板参数设置（右）

2.3 煤岩组合体应力与应变变化关系分析

强度作为衡量煤岩体破坏状态的重要力学性质之一，指煤岩体受到外荷载作用时抵抗破坏时的最大能力，为研究模拟过程中煤岩体所承受的最大破坏能力以及煤岩体加载过程中应力与应变之间的关系，以模拟过程中煤岩体的应变为自变量，煤岩体的应力作为因变量，得出煤岩体单轴加载过程中的应力—应变曲线如图5所示。

由图5 可以看出煤岩体应力随应变大致呈现先快速增加后逐渐减缓的趋势。随着加载时间的增加，固定荷载对煤岩体施加轴向应力使得煤体产生弹性形变，煤岩体内部应力随着应变快速增加，但随着应变的逐渐增大，应力的增长速度逐渐减缓，在应变量达到5.0×10之前，整个煤岩体处于一个弹性变形阶段，当应变达到1.0×10时，煤体应力随应变增加具有先小幅下降后小幅上升的趋势，此时煤岩体处于屈服阶段，当应变增加至1.5×10时，应力基本维持在一个较高的水平，其数值大致为32.5 MPa，即该煤岩体所承受的最大破坏能力约为32.5 MPa，若超过该值，则煤岩体发生延性破坏。

图5 应力—应变曲线图

3 结论

通过FlAC3D 的数值模拟结果可以看出，随着加载时间的增加，煤岩组合体内部应力及位移变化具有较强的规律性：

（1）煤岩组合体试件在Y 轴方向上的位移上部呈下沉式形变，底部呈现上升式形变，且顶底面形变量最大，向试件中部位置依次衰减。

（2）煤岩组合体在加载过程中应力变化呈现特征较为明显的Ⅰ、Ⅱ、Ⅲ三个阶段的变化，且Ⅰ阶段煤体与岩体试件具有应力方向相反的应力集中区域，Ⅱ阶段煤体应力与Ⅰ阶段媒体应力相反，Ⅲ阶段岩体应力集中区域消失，同时应力呈现出由岩体传递向煤体，由煤体两侧传递向煤体中心的趋势。

（3）煤岩组合在加载的过程中，煤岩体应力随应变大致呈现先快速增加后逐渐减缓的趋势，且随着加载时间的增加，煤岩体经历了弹性变形、屈服阶段、延性破坏三个变形阶段。