考虑应变软化的地下采场开挖变形稳定性分析

何忠明^{1, 2}，曹  平²

(1. 中南大学 矿业工程博士后流动站，湖南 长沙，410083；

2. 中南大学 资源与安全工程学院，湖南 长沙，410083)

摘  要：为了分析某大型地下采场开挖稳定性，利用FLAC^3D同时考虑岩土体的应变软化特性，计算开挖过程中的变形和破坏情况。结果表明：随着开挖的进行，周围岩体的位移不断增大，开挖步1_3，2_3和4引起的位移增量最大；开挖步1_1~1_2引起主采场顶板破坏高度为40 m左右，当开挖步为1_3时，顶板破坏高度达到80 m左右；2_1和2_2开挖步引起的破坏区面积变化不大，而2_3开挖步引起破坏高度迅速增长，达到120 m；开挖到第4步时，由于主采场与1和2采场空区相互贯通，破坏高度同样迅速增大，达到150 m左右；此后，由于影响稳定性的开挖阶段结束，破坏高度增量减小，最终破坏高度为160 m左右，得到的结果能够为工程实践提供指导。

关键词：地下采场；变形稳定性；数值分析

中图分类号：TU457         文献标识码：A         文章编号：1672-7207(2008)04-0641-06

Deformation and stability analysis of underground stope after

 excavation considering strain softening

HE Zhong-ming^{1, 2}, CAO Ping²

(1. Postdoctor Station of Mining Engineering, Central South University, Changsha 410083, China;

2. School of Resources and Safety Engineering, Central South University, Changsha 410083, China)

Abstract: In order to study the three dimensional stability of one large underground stope after excavation, the numerical software FLAC^3D was used to simulate the deformation and failure situation of the surrounding rock mass. The characteristics of strain softening were considered during calculation. The calculation results show that during the procedure of excavation, the displacements of surrounding rock mass increase gradually, the excavation steps of 1_3, 2_3 and 4 lead to the largest increment of displacement; the failure height of stope after excavation steps of 1_1-1_2 is about 40 m, while the excavation of 1_3 leads to the failure height of 80 m. The excavation steps of 2_1 and 2_2 lead to few increments of displacement, while the 2_3 excavation step leads to great increment of displacement, which is 120 m displacement for this excavation step. The main stope and 1, 2 stopes connect with each other after the fourth excavation step, which leads to the great increase of plastic zone areas, and the height of failure is about 150 m. Then after the fourth excavation step, the displacement increment is small for the same step. That is because the excavation step which can influence the stability of gob areas has finished. At last, the failure height reaches about 160 m. The results obtained from calculation can give guidance to the real practice.

Key words: underground stope; deformation stability; numerical analysis

某大型铁矿属于缓倾斜条带状厚大矿床，矿体埋藏深、呈多层产出且矿岩坚硬，其资源丰富，储量可靠，开采条件良好，服务周期长，矿石易选，精矿品质优良。由于该铁矿开采范围大，服务年限长，属深部开采(开采深度在850 m左右)，其突出的问题是应力高，变形大，地压控制难。为了更好地指导生产，确保井下作业安全，充分回收矿产资源，需对矿山基建和生产中矿区开挖稳定性进行相应研究。由于采矿工程变形稳定性是十分复杂的岩体工程动态问题，目前岩石力学研究中已有的理论解，只能解决圆形或椭圆形坑洞等简单形状的问题^[1]，而对于矿床开采这样复杂的空区显得无能为力。因此，许多学者寻求其他分析方法，如郭立等^[2]基于岩石工程系统理论(RES)的基本原理，构造了地下工程岩体失稳模式判别的智能预测模型；高谦^[3]提出了基于块体理论的可靠度分析理论和系统突变失稳预测方法；并建立了围岩突变失稳预测的可靠度分析方法；李晓静等^[4]采用层次分析法分析了地下洞室围岩稳定性的各个影响因素；董志宏等^[5]采用数值流行元方法对地下洞室稳定性进行探讨。这些方法一般从数学的角度对问题进行分析。近年来，随着计算机技术的不断发展，数值计算方法已成为岩石力学研究和工程计算的重要手段，并出现了许多种数值分析方法，如有限元法、快速拉格朗日差分法(FLAC^3D)、离散单元法、DDA法等^[6-13]。其中，快速拉格朗日元法(FLAC^3D)适用于绝大多数的工程力学问题，尤其适用于材料的弹塑性、大变形分析和施工过程的数值模拟^[11-13]。本文作者运用FLAC^3D同时考虑岩土体的应变软化特性建立数值分析模型，从宏观角度揭示出采场开挖后不同区域的破坏模式和变形情况，以便为工程实践提供指导。

1  考虑应变软化的数值方法

建立应变软化模型，并分析其在数值计算过程中的实施情况。该模型是Mohr-Coulomb模型的一种特殊形式，两者之间的不同之处在于：在应变软化模型中，由预先定义的硬化参数，根据分段线性原则，在塑性应变产生后，部分或所有单元的屈服参数，如黏结力、内摩擦角、剪胀角和拉伸强度都可发生变化。在每一个时间步内，总的塑性剪应变和拉应变都会被增量硬化参数校验，然后模型参数会调节到与自定义方程一致。该准则在屈服面上，剪切失效应力点的位置由非关联流动准则决定，拉伸失效应力点的位置由关联流动准则决定。在主应力空间σ₁-σ₃平面内，剪切失效的包络线f ^s=0由Mohr-Coulomb屈服准则可表示为^[14]

拉伸失效包络线f_t=0，由拉伸失效准则可表示为：

应变软化模型中的破坏准则、屈服函数、势函数、塑性流动准则和应力校正等都与相应的Mohr-Coulomb模型一致。在应变软化模型中，对于每个单元都要定义2个硬化参数k^s和k^t，分别作为塑性剪应变和拉应变的增量度量的和。

在FLAC^3D计算的基本四面体单元中，可计算出单元剪切和拉伸硬化增量。对于一个指定的四面体，其剪切硬化增量是由塑性剪应变增量张量的第二不变量来定义的，即：

拉伸硬化增量则由塑性拉伸应变增量表示为：

2  计算模型

2.1  工程概况

矿段工程地质类型属于以坚硬半坚硬岩层为主，工程地质条件属于中等条件的矿床，该矿山的采矿方法为崩落法，其主要特点如下：

a. 矿段位于复杂的区域构造复合地区，岩层由沉积岩，火成岩和变质岩构成，存在沉积间断、火成、构造和次生等多种结构面，处于地震活动区，有迹象表明，有的构造破碎带目前依然处于地应力活动状态。

b. 矿体及其顶板为火山岩和变质岩，大部分属于坚硬岩石，少部分属于半坚硬岩石，比较完整，因而比较稳定。岩层风化带深，但是，风化程度弱。辉长辉绿岩侵入体频繁穿插，接触带裂隙发育；Ⅰ号矿带间接顶板方柱石白云石大理岩，局部地段裂隙发育，岩心破碎；存在软弱结构面F3，局部存在泥化夹层、风化夹层和松软夹层；它们的力学强度有较大降低。不良工程地质因素对矿床开采有影响。

c. 铁矿石单轴抗压强度属半坚硬岩类。岩芯取芯率较高，岩(矿)石裂隙虽发育，多被后期脉石矿物充填或再生胶结。已开拓的大量坑道很少需要支护，岩体属局部稳定地段。

根据室内岩土试验，并由相应的岩体力学参数的工程处理，得到研究范围内相关岩土的物理力学工程特性指标如表1所示。

表1  岩体计算参数

Table 1  Calculation parameters for rock mass

由试验结果可知，发生塑性变形后，岩土体相应的黏结力c_p和内摩擦角与原始黏结力c_i和内摩擦角的关系为：；。其中，w_c和为黏结力和内摩擦角的变化因子。w_c和与塑性应变的关系见表2。

2.2  模型介绍

首先由SURPAC建立三维实体、块体模型，然后，

表2  w_c和w_φ与塑性应变ε^p的关系

Table 2  Relationship between w_c, w_φ and ε^p

利用自编的SURPAC-FLAC^3D接口程序建立FLAC^3D三维地质力学模型(由于SURPAC单元细分，导致某些区域网格不对齐，但细分单元与原始单元成整数倍的关系，因此，可通过FLAC^3D中的attach face命令实现网格不对齐区域的信息传递)，如图1所示。模型共33 700个节点，19 780个单元，整体边界尺寸(长×宽×高)为：1 320 m×1 080 m×620 m，采场总体尺寸为：676 m×342 m×305 m。模型包括以下几个部分。

a. 岩体：ptdm4，ptdm3，ptdh3，ptdh2，ptdh1，λω；

b. 矿体：磁铁富矿、磁铁贫矿、铜矿；

c. 采场(按开挖顺序排序)，1_1开挖步：主采区480分段空区_1_1，主采区460分段空区_1_1，主采区440分段空区_1_1，中部2采区540分段空区，中部2采区560分段空区，中部2采区590分段空区，中部1采区625分段空区，中部1采区645分段空区，中部1采区675放顶空区；1_2开挖步：主采区480分段空区_1_2，主采区460分段空区_1_2，主采区440分段空区_1_2；1_3开挖步：主采区480分段空区_1_3，主采区460分段空区_1_3，主采区440分段空区_1_3；2_1开挖步：主采区460分段空区_2_1，主采区440分段空区_2_1；2_2开挖步：主采区460分段空区_2_2，主采区440分段空区_2_2；2_3开挖步：主采区480分段空区_2_3，主采区460分段空区_2_3，主采区440分段空区_2_3；3开挖步：中部2采区520分段空区；4开挖步：主采区420分段空区，中部2采区480分段空区，中部2采区500分段空区，中部1采区605分段空区；5开挖步：中部1采区580分段空区；6开挖步：中部1采区560分段空区；7开挖步：主采区400分段空区。

图1  三维计算模型

Fig.1  Three dimensional calculation model

2.3  边界条件

模型上表面边界采用法向应力约束，底部采用固定位移约束，4个侧面采用法向位移约束。

2.4  地应力施加方法

在初始应力场的取值上采用平均构造应力场和重力场迭加，并将所有初始应力投影到模型坐标系后，形成计算中的初始应力场。为了较真实地模拟开挖变形过程，模拟分2步加载：第1步，仅考虑岩体自重情况；第2步，清除第1步产生的岩体位移，以模拟开挖过程中的应力变形状态；计算收敛准则为不平衡力比率(节点平均内力与最大不平衡力的比值)满足10^-5的求解要求。

2.5  监测点布置

为便于模拟计算的结果后处理与分析，在模拟矿体开挖过程中，首先垂直于采场进路方向布设监测线，监测线之间的距离为70 m，共布设P01~P07监测线，并且考虑到顶板处的几何形状不规则，监测线的位置高于顶板20 m布置；然后，沿进路方向位置以每隔20 m的剖面在监测线上截取监测点，截取46个剖面，共计322个监控点，如图2所示，以监控不同位置点位移变化规律。

图2  沿采场进路方向监测点布置

Fig.2  Location of monitoring points along excavation direction of stope

3  计算结果与讨论

分别从开挖过程中的位移变化规律、塑性区分布规律进行讨论。另外，为便于分析和文字说明，分析主要针对2个剖面，即过开挖部分中央的横、纵剖面，其中，横剖面为每步开挖后垂直空区进路方向剖面，各个剖面为动态剖面，不同开挖步的位置不同。

3.1  变形分析

3.1.1  开挖过程中监测点竖直位移

选取P02监测线为标准，对比其在开挖过程中的竖直位移变化，得到图3。从图3可以看出，随着开挖的进行各个监测点的位移不断增大，开挖步1_1引起的位移增量梯度不大。当开挖到1_3时，位移曲线出现明显的跳跃，位移增量为50 cm左右。继续开挖，位移增量梯度减缓。当开挖到2_2时，又一次出现位移曲线的跳跃。可见，对同一开挖步，开挖初期即使开挖量相同，其引起的位移增量也是不相同的。当空区拉开一定距离后，很小的开挖量就能引起很大的位移增量，即存在位移突变的临界开挖范围。从图3还可见，开挖步1的临界增量出现在1_3开挖步；另一个临界增量出现在开挖步5。开挖步5~7引起位移持续增大，但其趋势减缓。这是由于此时主采场顶板附近的开挖基本完成。另外，图中还显示，当监测点与起始监测点的距离超过200 m后，又迅速上升，然后再回落—上升，这是由于主采场空区开挖分2步。两部分空区高度不同，都引起一定位移，主采场空区1引起的位移较大。

图3  不同开挖步P02监测线竖直位移

Fig.3  Vertical displacement of P02 monitoring line for different excavation steps

3.1.2  开挖过程中监测点水平位移

不同开挖步P02监测线水平位移如图4所示。可见，在开挖过程中，沿采场进路方向监测点位移先增大后减小。最大水平位移值为0.51 m，较竖直位移的1.852 m小，这是由于空区形态为扁平型，水平方向范围远大于竖直方向范围，因此，位移形式以竖直位移为主。同一监测点随着开挖的进行，水平位移不断增大，但增加的梯度不同。开挖步1_3引起的水平位移增量最大，其次是开挖步4。而开挖步1_1和1_3~3引起水平位移较小，开挖步4~7位移基本不增长。这是由于1_3开挖范围超过空区临界开挖跨度。

图4  不同开挖步P02监测线水平位移

Fig.4  Horizontal displacement of P02 monitoring line for different excavation steps

3.2  稳定性分析

图5所示为开挖过程中，空区周围岩体的破坏状态。从各采场的横、纵剖面图中可以看出：

a. 主采场空区顶板主要发生拉剪破坏，空区正上部顶板发生拉伸破坏，再往上部分主要发生剪切破坏，且开挖初期主采场空区周围岩体破坏区域主要位于两侧边帮和顶板处，随着开挖的进行，顶板破坏区域面积明显增大，而边帮破坏区域面积基本不变。采场空区1和2周围岩体的破坏区域主要集中在两侧边帮，与主采场破坏区域分布形式不同。

b. 开挖过程中，空区顶板首先拉裂，出现冒落现象，其上部岩体随着顶板冒落也不断往下变形，从而出现一定高度的破坏岩体，此部分岩体以上部分由于应力重分布及其以下岩体冒落完毕而处于稳定状态。从图5可见，开挖步1_1~1_2引起主采场顶板破坏高度在40 m左右，当开挖步为1_3时，主采场顶板破坏高度达到80 m左右。可见，开挖步1_3引起塑性区面积的突然增大，其对稳定性的影响最大，与位移趋势的分析结果相同。由于采场空区1和2为长条形，其破坏区域位于边帮，厚度为20~30 m。开挖步2_1引起的主采场破坏高度为60 m左右，开挖步2_2引起的破坏区面积变化不大。而开挖步2_3由于主采场整体开挖，空区水平方向面积达到最大，同样引起了破坏区域面积的迅速增长，此时破坏区高度为120 m。当开挖步为4步时，由于主采场与采场空区1和2相互贯通，因此，塑性区同样迅速增大，此时破坏区高度为150 m左右。此后，由于影响稳定性的开挖阶段已经结束，这些开挖步引起的破坏区域面积的增量较小。最终，开挖完毕破坏高度为160 m左右。

(a) 第1_1步采场横剖面；(b) 第1_3步采场横剖面；

(c) 第2_3步采场纵剖面；(d) 第4步采场横剖面；(e) 第7步采场横剖面

图5  开挖过程中采场塑性区分布

Fig.5  Distribution of plastic zone during excavation

4  结  论

a. 随着开挖的进行各个监测点的位移不断增大，当开挖步为1_3时，位移曲线出现明显的跳跃，当开挖步为2_3时，又一次出现位移曲线的跳跃。在开挖过程中，沿进路方向监测点位移先增大后减小。同一监测点随着开挖的进行，水平位移不断增大，但增大的梯度不同。开挖步1_3引起的水平位移增量最大，其次是开挖步4。

b. 在开挖过程中，空区顶板首先拉裂，出现冒落现象，其上部岩体随着顶板冒落也不断往下变形，从而出现一定高度的破坏岩体，此部分岩体以上部分由于应力重分布及其以下岩体冒落完毕而处于稳定状态。主采场空区顶板主要发生拉剪破坏，空区正上部顶板发生拉伸破坏，再往上部分主要发生剪切破坏，且开挖初期主采场空区周围岩体破坏区域主要位于两侧边帮，随着开挖的进行，顶板破坏区域面积明显增大，而边帮破坏区域面积基本不变。

c. 开挖步1_1和1_2引起主采场顶板破坏高度大约在40 m左右，当开挖步为1_3时，主采场顶板破坏高度达到80 m左右。开挖步2_1引起的主采场破坏区高度为60 m左右，开挖步2_2引起的破坏区域面积变化不大。而开挖步2_3由于主采场整体开挖，空区水平面积达到最大，同样引起了破坏区域面积的迅速增长，此时破坏区高度为120 m。当开挖步为4步时，由于主采场与采场空区1和2相互贯通，因此，塑性区同样迅速增大，此时破坏区高度为150 m左右。此后，由于影响稳定性的开挖阶段已经结束，因此，这些开挖步引起的破坏区域面积的增量较小。

参考文献：

[1] Hoek E. Rock engineering[M]. North Vancouver: Evert Hoek Consulting Engineering Inc, 2000.

[2] 郭 立, 吴爱样, 周科平, 等. 基于RES理论的岩体失稳模式判别及其智能实现[J]. 中国有色金属学报, 2003, 13(3): 749-753.
GUO Li, WU Ai-xiang, ZHOU Ke-ping, et al. Pattern recognition and its intelligent realization of probable rock mass failure based on RES approach[J]. The Chinese Journal of Nonferrous Metals, 2003, 13(3): 749-753.

[3] 高 谦. 地下大跨度采场围岩突变失稳风险预测[J]. 岩土工程学报, 2000, 22(5): 523-527.
GAO Qian. Instability foreca st and risk evaluation of the surrounding rock masses for a large space stope[J]. Chinese Journal of Geotechnical Engineering, 2000, 22(5): 523-527.

[4] 李晓静, 朱维申, 陈卫忠, 等. 层次分析法确定影响地下洞室围岩稳定性各因素的权值[J]. 岩石力学与工程学报, 2004, 23(s2): 4731-4734.
LI Xiao-jing, ZHU Wei-shen, CHEN Wei-zhong, et al. Determining weight of factors in stability analysis of underground caverns by analytic hierarchy process[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(s2): 4731-4734.

[5] 董志宏, 邬爱清, 丁秀丽. 基于数值流形元方法的地下洞室稳定性分析[J]. 岩石力学与工程学报, 2004, 23(s2): 4956-4959.
DONG Zhi-hong, WU Ai-qing, DING Xiu-li. Stability study of underground cavern with numerical manifold method[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(s2): 4956-4959.

[6] Thomas K, Gunther M. A numerical study of the effect of soil and grout material properties and cover depth in shield tunnelling[J]. Computers and Geotechnics, 2006, 33: 234-247.

[7] Jinoh W, Kwangho Y. Coupled effects in stability analysis of pile-slope systems[J]. Computers and Geotechnics, 2005, 32: 304-315.

[8] Laouafa F, Darve F. Modelling of slope failure by a material instability mechanism[J]. Computers and Geotechnics, 2002, 29: 301-325.

[9] 姜清辉, 周创兵, 罗先启, 等. 三维DDA与有限元的耦合分析方法及其应用[J]. 岩土工程学报, 2006, 28(8): 998-1001.
JIANG Qing-hui, ZHOU Chuang-bing, LUO Xian-qi, et al. Coupling method of three-dimensional DDA-FEM and its application[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(8): 998-1001.

[10] 黄晚清, 陆 阳. 散粒体重力堆积的三维离散元模拟[J]. 岩土工程学报, 2006, 28(12): 2139-2143.
HUANG Wan-qing, LU Yang. 3D DEM simulation of random packing of particulates under gravity[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(12): 2139-2143.

[11] 林 杭, 曹 平, 宫凤强. 位移突变判据中监测点的位置和位移方式分析[J]. 岩土工程学报, 2007, 29(9): 1433-1438.
LIN Hang, CAO Ping, GONG Feng-qiang. Analysis of location and displacement mode of monitoring point in displacement mutation criterion[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(9): 1433-1438.

[12] 沈金瑞, 林 杭. 多组节理边坡稳定性FLAC3D数值分析[J]. 中国安全科学学报, 2007, 17(1): 29-33.
SHEN Jin-rui, LIN Hang. The numerical analysis on the stability of slope with several sets of joints by FLAC3D[J]. China Safety Science Journal, 2007, 17(1): 29-33.

[13] 王贵成, 曹 平, 林 杭, 等. 用灰色理论确定边坡最优监测点及安全系数[J]. 中南大学学报: 自然科学版, 2007, 38(3): 574-578.
WANG Gui-cheng, CAO Ping, LIN Hang, et al. Determination of optimum monitoring point of slope and safety factor using grey system theory [J]. Journal of Central South University: Science and Technology, 2007, 38(3): 574-578.

[14] 彭文祥, 赵明华, 袁海平, 等. 基于拉格朗日差分法的全长注浆锚杆支护参数优化[J]. 中南大学学报: 自然科学版, 2006, 37(5): 1002-1007.
PENG Wen-xiang, ZHAO Ming-hua, YUAN Hai-ping, et al. Parameters analysis of grouted bolts by Lagrangian difference method[J]. Journal of Central South University: Science and Technology, 2006, 37(5): 1002-1007.

收稿日期：2007-10-20；修回日期：2007-12-12

基金项目：国家自然科学基金资助项目(50774093)；中南大学博士后科学基金资助项目(200717)

通信作者：何忠明(1980-)，男，湖南永兴人，博士后，从事岩土工程的教学与科研工作；电话：13055178045； E-mail: hezhongming45@126.com