基于3DMine-Phase2分析复杂空区稳定性

杨 宁 ，尹贤刚 ，林 尧

（1.长沙矿山研究院有限责任公司，湖南 长沙 410012；2.金属矿山安全技术国家重点实验室，湖南 长沙 410012；3.国家金属采矿工程技术研究中心，湖南长沙410012；4.洛阳坤宇矿业有限公司，河南 洛阳 471700）

0 引言

采空区稳定性分析可通过现场观测、模型试验以及数值模拟等手段进行[1-5]。第一种方法即在现场建立观测点，施工时工作人员进行追踪，察看其应力-应变随施工而变化的情况，该法可记录现场最真实的资料，总结变化规律得出的结论可指导后续工作，同时能够作为反分析的基础数据源，因此该法在有些情况下显得尤为重要。第二种方法需靠大量的资源以及复杂的结构条件去支撑其实现与真实情况一致的简易模型，通过对施工步骤进行模拟得知其发生改变的动态过程。该法取得的效果视觉冲击力大，但在试验过程中所受约束因素太多，不易实现。随着科学技术的进步，渐渐被方案三代替。数值模拟法是通过电脑中输入的各项数据模拟实际情况开展的一项试验。适应能力强，可实现多种复杂模型的模拟工作，可随时改动，变通性强，成本也比较低，因此受到广大学者和生产企业的青睐[6-8]。

研究结合某矿地下采空区的实际情况进行三维数值模拟，主要用于对采空区稳定性及其地压活动规律进行模拟研究，从而对该矿在开采扰动下的矿柱、上下盘围岩的应力、位移状态等进行模拟分析。

1 工程概况

该矿矿体形态为层状、似层状与透镜状，倾角为25°～88°，厚度为1～108 m。矿体赋存于下盘大理岩与上盘绿泥石英片岩之间或赋存于大理岩、绿泥石英片岩之中。矿山采用平硐+竖井+斜坡道联合开拓，留矿法和分段空场法开采，地表允许塌陷。目前8中段及以上已回采完，三柱已及时回收，并崩落围岩形成了缓冲垫层。9～11中段正在进行回采出矿。根据矿山规划，10中段以下逐步转为充填法开采。

矿山经多年开采，已在地下多中段形成了大小、形态各异，多中段联通或单独的采空区，采用3DMine软件主要依据矿山提供的实测图件以及各中段的采空区分层平面图，从图中提取了空区实测线，并将其导入转换到真三维空间中对采空区进行了三维数值模拟。经与矿方技术人员讨论后决定：（1）在各采空区分层（中段）平面图上，当同一矿体存在多个相邻采空区时，空区之间的间柱已回采，采空区按照贯通处理，按各分层（中段）内采空区的并集参与建模；（2）在高度方向上，当同一矿体上下有相邻采空区时，根据矿山采矿方法及三柱回采情况，采空区之间的顶底柱已回采，相邻的采空区按照上下对应采空区平面图直接相连处理。由模型统计，该矿山开采30多年来，因采矿共形成的贯通性采空区共277个，体积总量约782万m3，但是由于采矿时间长，大量采空区已经冒落，剩余有空区和不确定是否冒落的空区总量127.3万m3，是该矿采空区治理的重点。空区三维模型如图1所示，采空区稳定性分析也是基于该三维数值模拟的结果。

图1 采空区立体模型及地表复合图Fig.1 Three-dimensional model and surface composite diagram of the mine goaf

2 有限元模拟分析模型

2.1 岩体力学参数

根据某矿工程地质特征，主要岩组为大理岩及片岩，由岩体结构特征及岩体的分级指标，经工程处理后综合选取的岩体力学参数如表1所示。

2.2 力学模型与边界条件

据工程地质、采矿方法特点及已存的采空区相关基础资料，为模拟该开采过程、围岩受采动影响的岩层移动过程、采空区的稳定性以及开采过程中采场的地压活动规律，反映多水平间采空区相互间的影响程度[9-13]，以采空区调查的三维数值模拟结果为依据，选取采空区分布较为密集的8线、31线、39线及48线地质剖面图，建立二维弹塑性有限元分析模型，各剖面的分析图见图2。从图中看出，因剖面的位置不同所切到采空区分布情况各有不同，不仅是标高不同，数量等也各有不同，相邻采空区的间距也不尽相同。为各采空区编号，即8线包括1～7号、31线包括8～20号、39线包括21～32号及48线包括33～47号。

因为前期并没有通过一些手段掌握矿区的原始地应力，则按照矿岩石自身的重应力场作为有限元模拟分析时的依据。为计算所创建的模型规模大于采空区开采规模的5倍，各计算模型范围为采空区开采范围的5倍以上，不同剖面在横向上的领域也各不一样，纵向上从地表至2 400 m之间。横向边界面上的所有节点在X方向上为确定的位置，即滑动铰支座；纵向2 400 m边界平面上的所有节点在X、Y方向上为确定的位置，即固定铰支座，地表是自由面，通过四节点等参单元展开分析。

2.3 模拟分析步骤

为简化起见，各模型的模拟开采步骤相同，第1步计算在自重应力场作用下的初始应力场，第2步计算在各中段已采空条件下的应力、位移、破坏区等分布特征。各剖面2个模拟计算步骤所划分单元和网格节点数分别为：8线在1934～1727间和1001～983间；31线在2555～2280间和1324～1286间；39线在 12255～9181间和 6254～5152间；48线在2403～2086间和1225～1196间。

表1 地下开采岩体力学参数Tab.1 Mechanics parameters for underground mining rock

图2 各勘探线剖面的有限元模型图Fig.2 Finite element model of each section

3 分析计算结果

3.1 应力及屈服区域

第一模拟步骤计算初始原岩应力场，从典型的31线剖面模型的主应力等值线（约定拉负压正）分布如图3所示，σ1、σ3分别代表最大、最小主应力。在模拟开采区域内矿岩的主应力等值线表示岩体内的应力大小随深度的增加而增大，呈水平状，在地表区域的主应力大小与方向受山形地势的影响，为沿山坡的顺坡方向。

第二模拟步骤计算后，在已采空部分的围岩处产生了应力重分布或拉、剪的破坏区，其中各剖面模型主应力等值线的分布状态如图4～图7所示。

图3 31线剖面模型的初始主应力等值线Fig.3 Initial principal stress contour map of line 31 profile model

图4 8线剖面主应力等值线及屈服区域Fig.4 Contour lines and yield areas of main stress in line 8

图5 31线剖面主应力等值线及屈服区域Fig.5 Contour lines and yield areas of main stress in line 31

图6 39线剖面主应力等值线及屈服区域Fig.6 Contour lines and yield areas of main stress in line 39

图7 48线剖面主应力等值线及屈服区域Fig.7 Contour lines and yield areas of main stress in line 48

3.1.1 8线剖面模型

在2 702 m以上存在7个单独或多中段联通的采空区。从应力分布及破坏区分布的情况来看，在该剖面最大的2号采空区上盘3 100 m以上形成了范围较大的拉破坏区，在其下与3号采空区相交处，有应力集中，且产生了拉、剪破坏区，但从整个剖面来看，应力集中最大部位为4号采空区与其下7号采空相交的部位（约在6中段、7中段），在7号采空区与5号采空区的顶板有拉破坏区，5号采空区的剪破坏也较明显。6号采空区目前的影响区仅在该采空区的周边，诱生应力一般为3.5～25.8 MPa。

除2 920 m的采空区周边存在压应力集中外，一般都在采空区周边形成应力松弛释放区，有较小拉应力值存在，最大值为0.31 MPa。上部2号大采空区会塌陷与地表贯通，其他采空区一般影响范围为15～25 m。在2 920 m附近采场产生的应力集中，最大压应力值为56.93 MPa，小于岩体的抗压强度，有压裂与其上部采空区联通的危险。

该剖面采空区稳定性最差的部位为：上部的4号采空区，其次为5号、7号采空区。

3.1.2 31线剖面模型

在2 702 m以上存在13个采空区，各采空区较为独立，相互间的影响较小。从应力分布及破坏区分布的情况整体影响来看：

3 102 m以上采空区对围岩稳定性的影响，没有其下部采空区对围岩稳定性的影响大。

一般是在采空区上下盘形成范围在5～28 m间的应力松弛区；在标高为3 032～3 100 m、2 890～2 945 m、2 948～2 964 m存在一个压剪区，以应力集中为主。

该剖面靠地表的9号采空区存在垮塌危险，会透地表垮塌。危险性最大的采空区为：2 762 m上的11号采空区、其上部16号采空区以及18号采空区，这3个采空区主要是联通的中段与采场较多，高度也较大。

在采空区周边诱生的压应力一般为2.5～18.7MPa，形成应力松弛释放区，有较小拉应力值存在，最大值为0.31 MPa。

3.1.3 39线剖面模型

在2 702 m以上存在12个采空区，各采空区较独立，26号采空区最大。从破坏区分布的整体来看：

从26号采空区到临地表的32号，其破坏区基本上是连在一起的，说明在该部位的采空区是垮塌在一起的。

大致从7中段以上，采空区的破坏是以张破坏为主，7中段以下在采空区的周边还存在剪切破坏区，且在采空区周边有应力集中现象，提示对该矿深部地压活动规律的认识需考虑到深部地应力存在。

在采空区周边诱生的压应力一般为2.5～20MPa，形成应力松弛释放区，有较小拉应力值存在，最大值为0.35 MPa。

3.1.4 48线剖面模型

在2 702 m以上存在15个体积较大的采空区，本剖面各采空区体积较大且相互距离较近一些，相互间的影响明显。从应力分布及破坏区分布的情况整体影响来看，在3 182 m的45号采空区周边破坏较大。

在3 108～2 905 m间的采空区（33号、35号、38号、41号～43号）等地段的破坏较集中与明显，该部位的采空区稳定性较差，应力的松弛释放也较明显。

在42号采空区右下角及37号采空区顶板、底板存在应力集中。

在采空区周边诱生的压应力一般为4.5～22.4 MPa，形成应力松弛释放区，有较小拉应力值存在，最大值为0.45 MPa。

3.2 位 移

第二模拟步骤计算后，各剖面位移分布受采空区位置与大小的不同，而呈现不同的分布特征，且位移量值也有不同，各剖面位移等值线与位移方向的分布状态如图8所示。

8线剖面位移最大的为上部2号采空区，表现为顶板下沉，最大下沉量为9.5 m，且由于上部采空区离地表较近，该部位的采空区会垮塌到地表，其他采空区位移较小。

图8 各剖面模型第二模拟步骤的总位移等值线与位移方向Fig.8 Contour lines and displacement directions of the total displacement in simulation step 2 of each profile model

31线剖面位移，受到下部11号采空区开采影响，这两个采空区位移量较大，最大值为1.3 m，位移方向为采空区上下盘沿空区方向。再从位移移动方向来看，本剖面的位移会发展到地表，从而引起地表开裂变形等情况。

39线剖面位移，主要受26～28号采空区影响，这些采空区体积较大，位移方向为沿空区的上下盘方向，指向空区，且会发展到地表形成塌陷。

48线剖面最大位移部位，与应力松弛区对应，在3 108～2 905 m间的采空区（33号、35号、38号、41号～43号）等地段，位移方向为沿空区的上下盘方向，指向空区，最大位移为0.57m。再从位移移动方向来看，该剖面的位移在地表也较小，地表有小的变形。

4 结论

通过3DMine软件联合二维弹塑性有限元方法对8线、31线、39线和48线地质剖面采空区的稳定性进行了应力、位移及屈服区的分析，结果表明：

（1）某矿采空区群的复杂、不规则地存在，使区域内应力分布异常复杂，从采空区周边的应力集中或应力松弛情况即可反映出来，有时是采空区两侧为应力松弛区，有时为采空区顶底板，空区间的顶底板位置存在应力集中等现象。一般在7中段以上的采空区以拉破坏为主，7中段以下的采空区的破坏还有剪破坏，且在采空区周边有应力集中现象，提示对该矿深部地压活动规律的认识需考虑到深部地应力的存在。

（2）从4个模型的2个模拟步骤的分析来看，本矿不稳定区域主要是根据采空区位置及采空区是否为位于同一矿带上，且上、下中段采场采空区直接或间接连成一片的情况而定。总体来看，采空区相对独立，相互的影响较小，若多中段采空区直接或间接连成一片，其稳定性较差，位移分布情况也类似。

（3）从4个剖面的稳定性分析来看，8线剖面会塌陷到地表，31线剖面会在地表形成开裂甚至塌陷等，39线剖面会塌陷到地表，48线在地表的位移较小，重点关注区逐渐从上部向下移动。