基于有限元数值模拟对青东矿10#煤底板破坏深度的研究*

张海涛 赵景宇

(1宿州学院资源与土木工程学院；2宿州学院安徽省煤矿勘探工程技术研究中心安徽宿州 234000)

随着青东矿10#煤开采的不断进行，研究矿区该煤层底板的采动破坏规律以及底板最大破坏深度可以为该矿区的安全开采提供科学的理论依据。有限元数值模拟主要用于分析岩石中的应力、位移、塑性区等特征量[1]，在中国通过建立数值模型研究开挖工作面及矿区情况的案例很多，其中有限差分法的使用最为常见，刘伟涛先生曾用这种方法研究煤层底板破坏因素对破坏深度的敏感性影响。根据半无限体理论以及利用经验公式计算法和塑性理论计算法针对底板采动破坏深度进行求解并以此开展一些研究也很常见[2-3]。在了解青东矿10#煤的水文地质条件以及其他地层信息的基础上获取力学参数，采用数值模拟软件FLAC3D对10#煤底板的破坏深度进行模拟，根据所做的数值模型可以直观的看出煤层塑性区破坏情况以及破坏深度，并判断出达到最大破坏深度的工作面位置也是最容易发生突水事故的地方。

1地质概况

青东矿位于安徽省淮北市濉溪县境内，交通便利，矿井四周分布共约200余条断层，除其中4条为逆断层外其余皆为正断层。矿区西部以F9断层为界，东临大刘家断层，北至太原组顶部K4灰岩露头线，南连F19断层，整个矿区面积约54km2[4]。根据地质资料显示区内主要含煤地层为二叠系，在研究区约883m地层厚度范围内共分布十个煤层，平均煤厚5.7m左右，其中10#煤层为五个开采煤层之一。该煤层平均厚度达到1.7m，定义为薄-中厚煤层，分布面积达到40km2，占整个矿区面积的75%，煤层结构简单，基本不见夹矸，除少部分受岩浆岩侵入影响煤层性质发生改变外整体变化不大，顶板以砂岩为主[5]。

2基本力学模型计算原理

通过FLAC3D软件编辑命令流，使用网格堆砌的方法建立数值分析模型，根据有限差分法进行数值计算。为了更加直观的判断矿井在开挖后塑性区的分布状态和破坏规律，选择摩尔-库仑本构模型作为理论基础[6-7]，该模型在岩土工程领域中相比于其他力学模型运用更为广泛和贴近实际[8]。另外矿体开挖的部分采用空单元(null)模型来充分模拟矿体的开采过程和采空区的充填过程，充填过程直接对空单元(null)模型重新赋值，当材料剪切面上的剪应力与正应力之比达到最大时，认为材料服从摩尔-库仑屈服原则，即材料将发生破坏屈服[9-11]，摩尔-库仑模型破坏准则的基本原理如图1所示。

图1 摩尔-库伦模型破坏原理

3结果与讨论

3.1建立数值模型

该研究范围主要为10#煤，具体建模尺寸需要考虑FLAC3D模型大小对数值模拟结果的影响：模型取得太小容易偏离矿山开采实际真实性，影响计算结果的可靠性；而如果取得太大，则使单元划分太多会影响计算速度，甚至计算机无法计算。由弹塑性力学理论-圣维南原理可知，在承受均匀载荷的无限大弹性体中开挖一圆孔后，孔边的应力状况将发生显著变化，但这种变化的影响范围实际上只限于附近的局部区域；在1倍孔径外的区域，应力只比开孔前的应力稍大、扰动较小；3倍孔径外的应力变化在工程上就可以忽略不计。因此，数值模型的计算范围可以根据开挖范围的尺寸关系进行合理确定。

根据上述原则，根据矿区内6013工作面地质原型、钻孔资料和综合柱状图加以简化建立工程地质物理模型。该模型为沿工作面走向的垂直剖面，边界条件采用左、右边界铰支，底边界固定的方式，上部施加相当于上覆岩土层自重(10.6MPa)的补偿荷载，模型长度450m，纵向最大尺寸60m。为使建立的数值模型适应FLAC3D计算要求，对研究区域进行加密处理，建立的模型范围x方向：0-1000；y方向：0-180；z方向：0-118，模型单总数51.84 万个，节点为 54.2901万个。使用FLAC3D软件编辑命令流，使用网格堆砌的方法建立的原始数值分析模型如图2所示，图中图例表示四种岩层属性，分别为10#煤、砂质泥岩、细砂岩、中砂岩。图3即开挖水平面模型图显示的是沿着开挖水平面所切的剖面，其中图例表示沿工作面分十步开挖的塑性区状态。建模所需的数据参数见表1。

图2 数值分析模型

表1 有限元模拟所需参数表

3.2结果分析

根据前面已建立的数值模型、原岩应力场条件下只考虑自重应力(10.6MPa)的垂直应力如图4所示。计算过程(10步开挖)的最大不平衡力监测曲线如图5所示，从图中可以看出塑性区的破坏由刚开始的不稳定逐渐趋近于平衡。数值计算过程中拉伸或者剪切破坏状态将直接在模型单元的塑性区分布状况中呈现出来，可以清晰的反映出煤层底板的破坏情况。

图4 原岩应力场条件下垂直应力云图

图5 计算过程(10步开挖)的最大不平衡力监测曲线

煤层开挖之前岩体处于平衡状态，原岩应力不发生改变，开挖以后岩体附近平衡状态遭到破坏，原有的应力增大，随着开挖面的向前推进，采空区不断扩大，采空区域逐渐被顶板上垮落下的岩块充填和压实。煤层平行开挖，按每步20m的距离向前推进以及开采，共分10步，为精简篇幅，只列出部分图形。

图6表示第一次开挖后塑性区的纵剖面分布状态，图例中none表示无塑性区，shear-n表示当前剪切破坏塑性区，shear-p表示过去剪切破坏塑性区，tension-n表示当前拉伸破坏塑性区，tension-p表示过去拉伸破坏塑性区。参照图例，可以从图7看出，煤层开挖区域的两侧以及尖角位置容易发生剪切破坏，开挖区域顶底板容易发生拉伸破坏。

图6 第1步开挖后塑性区(纵剖面)

图7表示第一次开挖后塑性区横剖面分布状态，提取的是底板最低层的破坏横剖面单元，从图中可以看出煤层底部破坏深度的影响范围。第一次开挖后，破坏深度约为1m。

随着工作面向前推进(步距20m)，绘制出第3、5、7、9、10步的开挖塑性区分布图。对比分析可知，随着开挖向前推进，沿开挖方向，塑性区逐渐增加，顶底板拉伸破坏和两端及尖角处的剪切破坏单元也逐渐增加。第3步开挖相对第1步而言，底板塑性区破坏单元明显要多。第5步和第7步开挖相对第3和第5步，也有同样的规律。在第7步开挖以后，虽然整体塑性区破坏单元明显还是呈递增的趋势，但底板塑性区破坏深度却不再增加。提取相关数据进行分析后，得出底板破坏的深度随开挖距离的关系如图8所示，底板破坏的深度最大值在12～13m之间。

图7 第1步开挖后塑性区(横剖面Z=74m)

图8 底板破坏深度随开挖距离的关系图

分析开挖后塑性区分布横剖面图，随着开挖的进行，受开采影响的围岩范围逐渐增大，沿煤层走向方向，煤层底板塑性破坏区大致呈一个正面“勺子”形状，随着开挖的进行，底板塑性破坏区逐渐增大，即逐步向勺底发展，之后塑性破坏区不再随着开挖的而继续增大，逐渐向勺柄发展。随着工作面的推进，煤层顶板塑性破坏区在开挖工作面四周分布较为均匀，在两端较为集中。由此可见，煤层开采后，开采的应力对工作面中部底板、开始切割位置和停止的位置附近底板的破坏程度较大。

综合以上分析可知，塑性区主要集中在采空区顶底板中，模拟计算过程正处于塑性流动状态且全部处于tension-n和shear-n塑性区，这表明煤层开挖过程中均出现了拉伸破坏和剪切破坏，但剪切破坏为底板破坏的主要方式。通过分析得出，当工作面推进140m(第7步开挖)时，底板破坏深度最大为12～13m，在此之后破坏深度不继续因为工作面的推进而产生变化，底板采动破坏深度不再纵向加深，但破坏范围仍继续横向扩大。由此可知，由数值模拟计算所获得的青东矿10#煤底板破坏深度最大值约为13m；当工作面推进到140m时，该位置最容易发生突水灾。

4结论

(1)根据建立的数值模型并通过数值模拟计算，得出青东矿10#煤底板的最大破坏深度约为13m；

(2)当工作面推进到140m时，出现青东矿10 #煤底板的最大破坏深度，该位置也最容易发生突水事故；

(3)根据得知的煤层底板最大破坏深度以及最容易发生突水事故的工作面位置，可以在相应位置加强支护，以减小围岩变形与垮落。