空区贯通对大型地下采场开挖变形稳定性的影响

林杭，曹平，余健，余巍伟

(1. 中南大学 资源与安全工程学院，湖南 长沙，410083；2. 南昌大学 科学技术学院，江西 南昌，330029)

采矿工程变形稳定性是十分复杂的岩体工程动态问题，目前岩石力学研究中已有的理论解，只能解决圆形或椭圆形坑洞等简单形状的问题[1]，而对于矿床开采这样复杂的空区显得无能为力。某大型铁矿属于缓倾斜条带状厚大矿床，矿体埋藏深且矿岩坚硬，其资源丰富、开采条件良好。由于该铁矿属深部开采(开采深度在850 m左右)，其突出的问题是应力高，变形大，地压控制难。为了更好地指导生产，确保井下作业安全，需对开挖稳定性进行相应研究。近年来，随着计算机技术的不断发展，数值计算方法已成为岩石力学研究和工程计算的重要手段[2-6]，其中快速拉格朗日元法(FLAC3D)适用于绝大多数的工程力学问题[7-9]，尤其适用于材料的弹塑性、大变形分析和施工过程的数值模拟。因此，本文作者运用 FLAC3D建立数值分析模型，从宏观角度揭示出采场开挖后，不同开挖方案对采场变形稳定性的影响。

1 工程概况

矿段工程地质类型属于以坚硬半坚硬岩层为主，工程地质条件主要特点为：(1) 矿段位于复杂的区域构造复合地区，岩层由沉积岩，火成岩和变质岩构成，存在多种结构面，处于地震活动区，有迹象表明有的构造破碎带目前依然处于地应力活动状态。(2) 矿体及其顶板为火山岩和变质岩，大部分属于坚硬岩石，少部分属于半坚硬岩石，比较完整，因而比较稳定。岩层风化带深，但是风化程度弱。存在软弱结构面F3，局部存在泥化夹层，风化夹层和松软夹层；(3) 铁矿石单轴抗压强度属半坚硬岩类；岩芯岩石质量指标(RQD)较高，岩(矿)石裂隙虽发育，多被后期脉石矿物充填或再生胶结。已开拓的大量坑道很少需要支护，岩体局部稳定。

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

2 计算模型

由于FLAC3D在建模以及单元网格划分等前处理问题上存在一定困难[10-12]，因此，本文首先利用ANSYS建立数值模型，然后，通过自编的 ANSYSFLAC3D的转换程序，导入FLAC3D进行计算(由于SURPAC单元细分，导致某些区域网格不对齐，但细分单元和原始单元成整数倍的关系，因此，可通过FLAC3D中的attach face命令实现网格不对齐区域的信息传递)。图1所示为三维计算模型，图2所示为采场三维模型。其中，对于软弱结构面的模拟采用实体单元建立模型，然后弱化单元参数，以达到模拟效果。模型共33 700个节点，19 780个单元，整体边界尺寸(长×宽×高)为：1 320 m×1 080 m×620 m；采场总体尺寸为：676 m×342 m×305 m。模型包括以下几个部分：(1) 岩体；(2) 矿体；(3) 采场。模型上表面边界采用法向应力约束，底部采用固定位移约束，4个侧面采用法向位移约束。在初始应力场的取值上采用平均构造应力场和重力场叠加，并将所有初始应力投影到模型坐标系后，形成计算中的初始应力场。为了较真实地模拟开挖变形过程，模拟分2步加载。第1步，仅考虑岩体自重情况；第2步，清除第1步产生的岩体位移，以模拟开挖过程中的应力变形状态；计算收敛准则为不平衡力比率(节点平均内力与最大不平衡力的比值)满足10-5的求解要求。

表1 岩体计算参数Table 1 Calculation parameters for rock mass

图1 三维计算模型Fig.1 Three dimensional calculation model

图2 采场三维模型Fig.2 Three dimensional model for stope

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

图3 沿采场进路方向监测点布置Fig.3 Monitoring location along stope route direction

3 计算结果与讨论

为便于分析和文字说明，分析主要针对2个剖面，即过开挖部分中央的横和纵剖面，其中，横剖面为每步开挖后垂直空区监测方向剖面。开挖中部2采区520分段空区和中部2采区500分段空区使主采场与2采场空区相互贯通，为了研究贯通与非贯通情况的不同，对比分析开挖(方案1)与不开挖这 2个采场情况下(方案2)空区周围岩体应力、变形和破坏区的响应。

3.1 应力变化规律对比

图4 横剖面上的最大主应力分布Fig.4 Maximum principle stress distribution on cross section plane

图5 纵剖面上的最大主应力分布Fig.5 Maximum principle stress distribution on longitudinal section plane

图4和5所示分别为开挖完毕后横、纵剖面上最大主应力的分布。从图4和5可知：纵剖面上，两者应力分布形式基本相同，只在顶板480 m和460 m开挖部分的拐角处拉应力分布略有不同，方案1的拉应力区略大于方案2，且其数值也较方案2大，其中方案1的拉应力为4.21 MPa，方案2的拉应力为4.18 MPa。对于横剖面上的应力云图：由于方案1主采区和2采区相互贯通，方案1顶板的拉应力区更靠近拐角处，而方案2拉应力区主要位于顶板中央，且方案2高压力区域面积大于方案1高压应力区域面积，说明空区间存在隔离层时，更有利于应力传递，空区整体更加稳定。

3.2 位移变化规律对比

图6所示为开挖完毕后横剖面上竖直位移的分布。从图6可见：在空区以内2种方案得到的位移等值线云图的分布形式基本相同，但其数值有较大差别，约为25 cm。在主采区和2采区贯通处，2种方案的位移形态有较大差别，方案1的位移趋势指向空区的拐角处，而方案2的位移趋势指向空区内部，并且位于拐角处的位移云图显示，方案1的位移值明显大于方案2的位移值。

图6 横剖面上的竖直位移分布Fig.6 Vertical displacement distribution on cross section plane

图7所示为开挖完毕后各监测线的竖直位移曲线。从图7可见：位移曲线的分布形式基本相同，即先增大再减小的趋势，最大值位于空区进路方向位置200 m左右。方案1的位移大于方案2的位移，两者之间的差别约为60 cm，另外从计算结果可知：2个方案之间的水平位移差别不大。说明是否开挖中部2采区520分段空区和中部2采区500分段空区对顶板沉降的影响较大。

3.3 塑性区分布规律对比

图7 监测点竖直位移Fig.7 Vertical displacement of monitoring point

图8所示为开挖过程中空区周围岩体的破坏状态。从图8可以看出：主采场空区顶板主要发生拉剪破坏，空区正上部顶板发生拉伸破坏，再往上部分主要发生剪切破坏，并且随着开挖的进行，顶板破坏区域面积明显增大，而边帮破坏区域面积基本不变。开挖过程中，空区顶板首先拉裂，出现冒落现象，其上部岩体随着顶板冒落也不断往下变形，从而出现一定高度的破坏岩体，此部分岩体以上部分由于应力重分布及其以下岩体冒落完毕而处于稳定状态；2种方案的塑性区分布基本相同，即沿着空区中央以弧状形式不断向外发展，但方案1空区顶板的破坏区域大于方案2空区顶板的破坏区域。方案1的破坏高度大约为160 m，方案2的顶板破坏高度为120 m左右，可见不开挖中部2采区520分段空区和中部2采区500分段空区能够大大提高空区的稳定性。

图8 横剖面上的塑性区分布Fig.8 Plastic zone distribution on cross section plane

4 结论

(1) 纵剖面上，两者应力分布形式基本相同，且数值也大体相同。横剖面上，方案1顶板的拉应力区更靠近拐角处，方案2拉应力区主要位于顶板中央，且方案2高压应力区域面积大于方案1高压应力区域面积，即空区间存在隔离层时，更有利于应力传递，空区整体更加稳定。

(2) 2种方案的位移等值线云图分布形式基本相同，但其位移有较大差别，约为25 cm。

(3) 2种方案的塑性区分布基本相同，即沿着空区中央以弧状形式不断向外发展，但方案1空区顶板的破坏区域大于方案2空区顶板的破坏区域。方案1的破坏高度约为160 m，方案2的顶板破坏高度为120 m左右，空区贯通对采场的稳定性有较大影响。

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

[2]Won J,You K,Jeong S,et al. Coupled effects in stability analysis of pile-slope systems[J]. Computers and Geotechnics,2005,32(4): 304-315.

[3]Cai F,Ugai K. Reinforcing mechanism of anchors in slopes: A numerical comparison of results of LEM and FEM[J]. Int J Numer Anal Meth Geomech,2003,27(7): 549-564.

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

[5]Kasper T,Meschke G. A numerical study of the effect of soil and grout material properties and cover depth in shield tunnelling[J].Computers and Geotechnics,2006,33(4/5): 234-247.

[6]Griffiths D V,Lane P A. Slope stability analysis by finite elements[J]. Geotechnique,1999,49(3): 387-403.

[7]Grasselli G. 3D Behaviour of bolted rock joints: experimental and numerical study[J]. International Journal of Rock Mechanics& Mining Sciences,2005,42(1): 13-24.

[8]LIN Hang,CAO Ping,GONG Feng-qiang,et al. The directly searching method for slip plane and its influential factors based on the critical state of slope[J]. Journal of Central South University,2009,16(1): 131-135.

[9]范文,俞茂宏,李同录,等. 层状岩体边坡变形破坏模式及滑坡稳定性数值分析[J]. 岩石力学与工程学报,2000,19(s1):983-986.FAN Wen,YU Mao-hong,LI Tong-lu. Failure pattern and numerical simulation of landslide stability of stratified rock[J].Chinese Journal of Rock Mechanics and Engineering,2000,19(s1): 983-986.

[10]林杭,曹平,李江腾,等. 基于SURPAC的FLAC3D三维模型自动构建[J]. 中国矿业大学学报,2008,37(3): 339-342.LIN Hang,CAO Ping,LI Jiang-teng,et al. Automatic generation of FLAC3D model based on SURPAC[J]. Journal of China University of Mining and Technology,2008,37(3): 339-342.

[11]罗周全,吴亚斌,刘晓明,等. 基于 SURPAC的复杂地质体FLAC3D模型生成技术[J]. 岩土力学,2008,29(5): 1334-1338.LUO Zhou-quan,WU Ya-bin,LIU Xiao-ming,et al. FLAC3D modeling for complex geologic body based on SURPAC[J].Rock and Soil Mechanics,2008,29(5): 1334-1338.

[12]廖秋林,曾钱帮,刘彤,等. 基于 ANSYS平台复杂地质体FLAC3D模型的自动生成[J]. 岩石力学与工程学报,2005,24(6): 1010-1013.LIAO Qiu-lin,ZENG Qian-bang,LIU Tong,et al. Automatic model generation of complex geologic body with FLAC3D based on ANSYS platform[J]. Chinese Journal of Rock Mechanics and Engineering,2005,24(6): 1010-1013.