司家营研山铁矿采场三维边坡稳定性数值模拟及协同监测分析

杨意德 杨天鸿 刘 洋 邓文学 李 华 田益琳 叶会师

(1.东北大学资源土木工程学院,辽宁 沈阳 110819;2.河北钢铁集团司家营研山铁矿有限公司,河北 唐山 063701)

露天矿边坡稳定性分析是边坡安全设计的关键环节[1],就目前边坡稳定性分析来看,多数专家、学者分析二维模型较多,其优点为可了解露天边坡某个剖面在水位、爆破振动及人工开挖等影响下的稳定性[2-4],分析多个剖面便可了解实际边坡一定范围内的稳定性状况,给研究者提供不少便利。然而,边坡岩层的走向和倾向影响边坡空间应力分布,尤其是在层状边坡中,边坡应力应变分析简化为平面应变问题分析有明显的不足之处,边坡破坏产生的滑坡绝大多数是复杂的三维问题。

在三维边坡稳定性分析方面,Duncan[5]曾从20篇文献中归纳了三维边坡稳定性分析的特点和局限性,曹兰柱[6]、宋子岭[7]等对不同的三维边坡分析方法作对比分析,前人研究表明三维有限元和三维极限平衡方法的应用均有诸多前提条件[8-9],如此一来,也就限制了三维稳定性分析理论的应用范围,从而使三维边坡稳定性分析不能有效发展和应用。相对极限平衡法,有限元强度折减法在分析边坡稳定性时具有的优势在其可自动搜索滑面,同时满足了力的平衡条件并且考虑了岩土体应力—应变关系、岩土体和加固结构的耦合作用,更为严谨。有限元强度折减法在二维稳定性分析中已经取得大量研究成果并且已经付诸于工程实践当中[10-14],其在三维边坡稳定性分析中也有很多应用。在有限元分析中,网格划分至关重要,网格质量的好坏直接关系到分析是否能够顺利、快速地完成,也关系到是否能够得到高精度的分析结果[15-17]。FLAC3D快速拉格朗日有限差分法方法是将差分方程近似表示微分方法的一种数值方法,由于其计算收敛性较好,在岩土边坡工程中应用广泛,许多学者通过工程实践进行了研究[18-21],但是当需要建立较复杂的三维模型时,前处理往往为主要工作[22],但由于FLAC3D并无配套的三维前处理建模及网格划分软件,因此,需借助其他的工具建模并划分网格并通过接口程序转换为FLAC3D可使用的数据格式。

以司家营研山铁矿为研究对象,基于矿山地质剖面建立三维地质模型和地表DTM地表模型,并借助Rhino软件及Griddle插件对复杂地表以及岩性界面进行多轨扫略放样、布帘简化、几何布尔运算以及有限元网格划分等工作,将现状边坡分步开挖至最终境界,进行FLAC3D三维数值模拟计算,从宏观角度分析边坡的潜在塑性区分布及其演化规律后,再借助协同监测系统对局部区域进行了验证,并得到潜在滑区的进一步认识,为边坡优化设计及矿山安全生产提供了依据。

1 工程地质概况及边坡滑坡现状

1.1 工程地质

司家营研山铁矿位于河北唐山滦州市,是国内大型露天铁矿之一。矿区被第四系地层大面积覆盖,基岩露头除在矿区东部和尚山～铁石山一带有较连续的分布外,其他均为零星出露。矿区采场内中元古界长城系大红峪组(Pt2chd)石英砂岩作为沉积盖层呈角度不整合覆盖在新太古界单塔子岩群白庙子组黑云变粒岩(Ar3gnt)、混合岩化黑云变粒岩(Ar3mgnt)、钾长石化云母片岩(Ar3sch)等变质岩系及太古代沉积变质铁矿(BIF)之上。除水平覆盖于岩层上部的第四系之外,其他地层均为倾斜岩层,总体倾向西,其中大红峪组石英砂岩倾角15°±5°,变粒岩系倾角40°±5°,构造线近南北,形成了东帮顺倾,西帮反倾的岩层结构。根据目前开采现状,东边坡由变粒岩系和第四系冲积层组成,西边坡由石英砂岩和第四系冲积层组成,而南北两边坡则由其东半部的变粒岩系和西半部的石英砂岩构成。

1.2 滑坡概况

近几年司家营研山铁矿发生过多次规模不一的边坡滑坡,滑坡的空间分布及部分滑坡现状照片如图1,其中规模较大的滑坡统计如表1所示。在经过现有资料统计和现场勘查后发现,司家营研山铁矿东帮主要由第四系土层、黑云变粒岩及白云母片岩组成,黑云变粒岩按照其风化程度的不同从上到下又可以划分成强风化黑云变粒岩,中风化黑云变粒岩及微风化黑云变粒岩,且东帮临近新河,白云母片岩遇水容易泥化,不同风化程度的黑云变粒岩分界面倾向与边坡倾向一致,极易发生顺层滑坡破坏。西帮边坡稳定性优于东帮,主要滑坡分布在第四系表土层。

图1 司家营研山铁矿历史滑坡分布Fig.1 Historical landslide distribution map of Sijiaying Yanshan Iron Mine

表1 司家营研山铁矿露天边坡历史滑坡区域统计Table 1 Statistics on the historical landslide area of the open slope of Sijiaying Yanshan Iron Mine

2 三维地质模型建立及网格划分

2.1 三维地质模型建立

基于矿山开采现状平面图、勘探线剖面图建立露天采场三维地质模型。详细的方法流程如图2(a),具体步骤为:①基于3Dmine中的坐标转换将所有剖面分布于各对应勘探线,形成各剖面在三维空间上的分布状态,提取各岩层分界线;②根据各岩层分界线,利用Rhino建模软件中NURBS曲面建模功能,通过多轨扫略放样生成各不同岩层间分界面;③根据现状、采剥计划和最终境界平面图,生成三角网DTM面,基于Rhino布帘工具建立简化的现状面、2022年采剥计划面以及最终境界设计平面(如图2(b)、图2(c)所示);④建立相应尺度的立方体,基于Rhino布尔运算的几何切割功能进行岩层及现状境界等分界面的切割,最后形成司家营铁研山铁矿三维地质模型,如图2(d)所示。

图2 司家营研山铁矿地质模型Fig.2 Geological model of Sijiaying Yanshan Iron Mine

2.2 数值网格划分

复杂有限元网格划分需保证界面之间网格节点的重合,即分界面两侧网格在该面上完全重合,此为三维复杂几何网格建立的关键。由于Rhino的几何分割存在一定的精度误差,相邻几何分界面实际上无法完全重合,本研究首先划分所有界面的表面网格(如图3(a)),在确保表面网格连续、分界面上网格重合后(基于Rhino二次开发的Griddle插件进行网格划分,可自动对一定范围内的两NURBS曲面上的网格进行合并),进一步划分三维网格。司家营研山铁矿模型的整体尺寸为长4 693m,宽为3 065m,高为1 066 m,采用四面体网格,共划分单元数为16 065 000个,节点数为2 830 020个,导入FLAC3D中如图3(b)所示,共计完成19个不同的复杂几何网格划分。

图3 有限元网格划分Fig.3 Finite element mesh generation process

3 基于FLAC3D的三维边坡稳定性计算

3.1 力学参数确定及边界条件

模拟所需岩体力学参数主要来源于工程勘察报告以及相关的文献中[23],本文不详细赘述。主要岩体的物理力学参数指标取值见表2所示。计算采用边界条件为:对底部边界采用位移约束,对四周采用法向位移约束。共分三步进行平衡应力计算,第一步采用弹性本构模型计算现状平衡应力场,第二步采用弹塑性本构模型,在上一步基础上开挖至2022年的排产境界,计算其平衡应力场并判断塑性区,第三步采用弹塑性本构模型,并在第二步基础上开挖至最终的设计境界,进行进一步的计算分析。

表2 岩体物理力学参数指标Table 2 Physical and mechanical parameters of rock mass

3.2 平衡应力计算

边坡稳定性计算主要取决于最大不平衡力收敛以及不平衡比率情况,FLAC3D默认当不平衡比率达到10-5计算达到平衡并停止计算。图4即为对现状边坡模型计算过程中的最大不平衡力与计算步的关系曲线,该曲线变化趋于稳定表明模型收敛较好,计算到一定的步数最大不平衡力便趋于零。图5(a)和图5(b)为现状模型典型剖面的Z方向应力云图,边坡垂直应力在高程上呈现渐变趋势,底部应力比顶部应力大,呈现等值线分布,单斜构造的岩层倾斜分布导致两侧应力分布的不对称,较符合现场实际状况。

图4 最大不平衡力与计算步的关系曲线Fig.4 Relation curve between maximun unbalance force and calculating time step

图5 现状Z方向应力云图剖面Fig.5 Stress in direction Z of mining status

3.3 塑性区分析

塑性区分布是判断边坡破坏的主要判据。根据开挖至2022年排产境界以及最终境界的三维塑性区分布图(如图6(a)、图6(b))可以看出,潜在的塑性破坏区域主要分布在N14勘探线西帮附近和N26勘探线东帮附近,与滑坡现状较为吻合,尤其是东帮N26线的塑性变形区域在2018年发生了较大的多台阶的顺层滑坡,危及矿山的安全生产。对比东西帮塑性区分布(图6(c)～图6(f))表明:①东帮塑性区的分布范围以及深度都远大于西帮,说明东帮的顺倾边坡稳定性小于西帮;②随着边坡开挖深部的不断更加,东帮的塑性区集中区域的分布范围和深度增加,贯通分布在整个最终境界边帮;西帮边坡塑性区域分布演化则由浅部第四系土层转移至深部的基岩上,但分布区域范围有限,因此西帮的整体边坡安全系数高于东帮。

图6 2022年排产境界及最终境界塑性区分布Fig.6 Distribution of plastic zone of 2022 production schedule pit limit and the ultimate pit limit

3.4 位移云图分析

X向(即东西向)的位移云图可以直观地反映在开挖过程中东西帮边坡向临空面方向的发展趋势。图7(a)～图7(f)表明:①随着开挖进行,东西帮的边坡向临空向位移的位移量逐渐增加,最大X向位移绝对值分别由0.392 m和0.321 m增加至1.267 m和1.313 m;②东帮向临空面的X向最大位移值大于西帮,变形区域遍布整个东帮,西帮边坡的变形区域则主要集中在边坡下部,这主要是由于东帮顺倾和西帮反倾的岩体结构,导致东西帮变形具有各向异性特征,较大的变形亦不利于东帮边坡的整体稳定性。因此,为确保矿山的安全生产,应对东帮变形及破坏区采取适当加固防治措施。

4 协同监测分析

4.1 协同监测介绍

通过数值模拟结果及现场情况来看,研山铁矿东帮有较大的滑坡风险,因此对研山铁矿东帮N26线进行应力计、微震协同监测,以进一步评估东帮边坡的稳定性,监测布置如图8所示。微震监测系统为6通道,在-157、-187 m各布置3个,锚杆应力计在-157、-187 m平台各布置1个孔,每个孔分别布置4个应力计,每2个应力计之间间隔为4 m。

4.2 数据分析

对研山铁矿东帮顺层边坡进行了2021-04-29—2121-08-27的协同监测,微震监测主要针对的是开挖工作对围岩的扰动情况,每次的爆破开挖都会对边坡围岩造成一定程度的损伤,而每次扰动都会有对应的微震信号,如图9(a)所示,不同灰度代表能量的大小,在-127 m以下有大量的微震监测信号,而浅部的微震信号多于深部,说明浅部的损伤大于深部,且在坡体内6～7m深度处出现较为集中的微震信号贯通,说明爆破振动在边坡表层形成了一定深度的损伤区。

图9 协同监测数据分析Fig.9 Analysis of collaborative monitoring data

如图9(c)所示,1～4号(-187 m)应力计中整体应力变化分别为 3.8、2.5、0.1、2.8 kN,降雨(7月13日)时最大应力变化分别为 4.7、1.9、1.6、1.3 kN;如图9(d)所示,5～8号(-157 m)应力计中整体应力变化分别为 5.6、6.9、1.1、0.8 kN,降雨时最大应力变化 1.3、1.9、0.6、0.6 kN。 东帮布置的8个应力计均有应力变化,尤其是在降雨期间传感器应力有明显变化,在雨后随着坡内静水压力下降后,应力又逐渐恢复原水平,且上下2组应力计均体现边坡浅部应力大于深部。综上所述,微震与应力计监测到的结果相对应,浅部的损伤与应力变化大于深部,爆破振动对边坡表层形成损伤区,且降雨对于边坡内部应力变化作用明显。

5 结 论

在收集司家营研山铁矿的大量地质剖面资料基础上,结合矿山采剥计划和最终境界设计方案,基于3Dmine和Rhino软件,完成了矿山三维地质模型的建立,并划分出FLAC3D可直接使用的有限元网格以作东西帮整体采场的三维稳定性分析,分两步开挖进行数值模拟计算,并对东帮局部进行协同监测验证及进一步分析,得到以下结论:

(1)利用 Rhino的 NURBS曲面建模功能及其Griddle插件的表面网格重合处理功能,可较好地实现露天矿复杂地质体的三维建模及其有限元网格划分,网格质量可保证计算收敛性。

(2)潜在的塑性及变形区域,在东西帮均有分布,这主要是由单斜构造影响下形成了边坡顺层与反倾2种截然不同的结构,应力、变形及潜在破坏范围均具有各向异性特征,主要分布于东帮边坡,而西帮边坡的整体稳定性高于东帮。

(3)通过微震、应力计协同监测分析可知,研山东帮N26线有较明显的损伤与应力变化,与上述数值模拟结果吻合,且降雨对于边坡内部应力变化作用明显,浅部的应力与损伤变化大于深部,应采取相应治理措施,尤其是浅部,以确保东帮边坡的长期稳定性。

(4)本研究存在不足之处,未充分考虑层状岩体原生各向异性的变形及强度差异,当前同一类型岩体力学参数按同一参数取值,因此,对采场整体三维边坡稳定性分析方法仍具有进一步改进的空间。