黄 敏,崔年生,危剑林,何建华,孔繁成,黄英华
(1.长沙矿山研究院有限责任公司,湖南 长沙 410012;2.金属矿山安全技术国家重点实验室,湖南 长沙 410012;3.福建省新华都工程有限责任公司,福建 厦门 361000;4.湖南有色新田岭钨业有限公司,湖南 郴州 423000;5.绍兴铜都矿业有限公司,浙江 绍兴 312000)
开展金属矿山地下开采引起的地表岩层移动变形规律研究,对评判地表建(构)筑物安全程度、圈定地表影响范围、缓解矿山与周边居民日益紧张的关系、征用土地、搬迁房屋以及保证矿山安全生产具有十分重要的意义[1-4]。当前,众多科研工作者都对其进行了探索和分析,取得了许多有价值的科研成果,如王艳辉等[5]采用正交试验及数值模拟计算方法对不同采高、不同冒落高度、不同岩体力学参数下的岩层移动规律进行了深入研究,具有一定的借鉴作用;夏开宗等[6]通过对具有陡倾结构面的程潮铁矿西区的地表变形监测资料及宏观破坏特征进行了分析,得出矿区岩层移动分为采空区顶板岩体破坏扩展至地表引起的塌陷阶段及采空区周边围岩向采空区的倾倒破坏阶段,并得出了岩体主要发生水平移动的变形规律,呈现出先缓慢变形,后快速变形的特征,两者之间有明显的转折点;赵海军等[7]在自重应力和构造应力两种工况下,采用数值模拟方法对急倾斜矿体开采的岩移特征与变形机理进行了详细分析,得出了构造应力条件下地表出现双沉降中心的,自重应力条件下只存在单沉降中心的现象,同时,两者在地表移动变形量、移动变形影响区规模及地表宏观变形破坏特征上存在较大差异;李文秀等[8]采用FLAC3D软件对金宝山铂钯矿地下采矿过程引起的地表变形和岩层移动变形规律进行了三维仿真模拟,得出了地下开采引起的地表岩层移动是一个连续变化的过程并预测了地表移动范围;袁仁茂等[9]采用理论力学及数值计算的方法对金川二矿区急倾斜厚大矿体地下开采引起的地表岩层移动传递模式及地表岩层移动规律进行了深入研究,并结合地表监测结果,得出矿体厚度对地表岩移的显现特征具有明显的影响;廉海等[10]运用FLAC软件对石人沟铁矿地下开采过程进行分步开挖模拟,得到了不同开挖阶段的地表动态移动曲线;陆清运等[11]等先利用地表监测数据与三维数值模拟计算反演出材料力学参数,后采用三维数值方法分析地下开采冒通地表与未冒通地表两种工况下的地表变形特征,得到了岩层移动角的变化趋势。
矿山地下开采对地表建(构)筑物的影响主要取决于地质和采矿等因素,如矿床的开采深度、采高、倾角、采矿工艺、地质结构、围岩特性等。为探究地下开采的影响范围,评估地表建构物安全程度,指导矿山安全生产,本文采用FLAC3D有限元差分法软件分析矿山开采后地表岩层移动变形情况,评判矿山地下开采对地表建(构)筑物的影响范围和程度。
哈图金矿属中型地下金矿山,采用地下开采方式,设计生产能力为1 200 t/d。齐Ⅰ金矿区地表和深部分布矿脉有近百条,+934 m中段以上主要有L27、L7、L27-8、L27-14、L27-15等35种矿脉,矿体分散但品位高。矿体走向近东西向,平均厚度2.6 m,局部厚度8~15 m,矿体北倾,倾角50°~70°,属热液充填交代型矿床。矿体上下盘围岩主要为蚀变玄武岩,少量蚀变凝灰岩及石英脉岩。矿山采用浅孔留矿法采矿,电耙出矿结构,阶段高度40 m,沿矿体走向布置采场,采场长30~50 m,顶柱高为5~6 m,底柱高5 m,人行天井宽1.2 m,人行天井两边留2~2.5 m矿柱。
矿区地表有少量的工业厂房、柏油路等建(构)筑物。井下采空区主要分布在934 m、974 m、1 014 m及1 054 m中段。
以各中段采空区影响范围为研究区域,具体建模尺寸需要考虑FLAC3D模型大小对数值模拟结果的影响,并兼顾矿山开采实际及计算机计算速度,根据弹塑性力学理论-圣维南原理合理确定。据此对矿区地表、矿体及各采空区数据进行提取和分析。最终建立的模型范围X方向:3 850~5 550;Y方向:5 100~5 800;Z方向:+750 m至地表。
由于矿体形态复杂,采空区数量众多,FLAC3D软件在复杂模型的构建上存在一定缺陷,因此采用AutoCAD、Dimine及Midas GTS NX等软件联合构建地质模型。Dimine软件仅用于三维实体模型生成,实现模型可视化。为了将三维模型应用到数值模拟仿真计算,选用Midas GTS NX软件作为媒介进行模型传递。Dimine模型经AutoCAD处理后导入Midas GTS NX中后进行面缝合、实体扩展、布尔运算、网格划分等操作。地质模型见图1。
数值分析必须将模型计算区域进行离散化处理,通过转换软件将模型的单元、节点和组信息导入FLAC3D软件中。为了合理规划模型单元数量、提高计算效率、适应FLAC3D软件计算要求,对各中段采空区及矿体进行加密处理,其中,各采空区单元尺寸3 m,围岩36 m、次级单元18 m,最终生成的数值分析模型如图2所示,模型单元总数约131.5万个,节点总数约22.5万个。
采用数值模拟分析矿山开采对地表及建构筑物的影响程度和范围。数值模拟计算的前提是确定材料本构模型、力学参数、边界条件等。
1) 本构模型。选择摩尔-库仑本构模型,另外,开挖部分采用空单元(null)模型。
2) 材料力学参数。确定模型的材料属性是FLAC3D软件计算的关键。根据矿围岩室内试验,考虑岩体节理裂隙等客观因素,并通过RocLab软件折减计算,得到各岩体力学参数,见表1。
3) 边界条件。数值分析中初始应力场的施加十分关键,施加有偏差的应力场会导致数值计算的失败。本次按岩体自重应力施加荷载,在模型底部施加Z向约束,四周分别施加X向约束,Y向约束,地表为自由面。
4) 计算方案。本次数值模拟根据矿区934~1 054 m中段实际开采区域进行模拟计算。整个模拟计算分五步进行,首先初始平衡,之后对1 054 m中段、1 014 m中段、974 m中段、934 m中段分四个步骤进行开挖计算。
图1 地质模型Fig.1 Geological model
图2 数值分析模型Fig.2 Numerical analysis model
表1 模型力学参数Table 1 Model mechanics parameters
岩性容重/(kN/m3)体积模量/GPa剪切模量/GPa抗拉强度/MPa泊松比内聚力/MPa内摩擦角/(°)矿体28.504.504.100.180.152.6148.35蚀变玄武岩28.106.096.670.220.274.2043.74
按照地表移动变形值大小及其对建(构)筑物及地表的影响程度,可将地表移动盆地划分为最外边界、移动边界和裂缝边界三个边界(图3)。对应描述地表移动的参数主要有边界角、移动角及裂缝角。
1) 最外边界:一般指下沉1 mm的点圈定的边界,如图3中的ACBD。
2) 移动边界:是指临界变形值确定的边界,根据《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规程》(简称“三下规程”)中规定:使用破坏等级与地表变形的对应关系最低等级的三个参数临界值来划定地表影响范围,即倾斜i=3.0 mm/m、曲率K=0.2 mm/m2、水平变形ε=2.0 mm/m。最终划定的移动边界如图3中的A’C’B’D’。
3) 裂缝边界:是指移动盆地的最外侧的裂缝圈定的边界,如图3中的A”C”B”D”。
根据“三下规程”,地下开采地表建筑物的损坏等级由地表水平变形、曲率、倾斜三个指标判定。地表变形的大小和建筑物本身抵抗变形的能力决定了地表建(构)筑物的损坏程度。地表倾斜反映地表移动盆地沿某一方向的坡度;曲率反映弯曲程度的变化;水平变形表示单位长度线段的拉伸或压缩。
4.2.1 地表位移分析
矿山地表整体位移等值线图如图4所示。为描述地表整体位移变化趋势,选用0.01 m为趋势边界。随着934~1 054 m中段之间多个采场回采后,岩层位移量逐渐增大,最大位移约为12 cm,主要集中在开挖区域的顶底板中央位置,主要为矿体的上盘位置。同时,随着开挖的进行,位移影响范围也逐渐增大,最终在地表形成位移约为1 cm的影响区域。该区域主要集中在L7脉板房至水塔之间,东西范围约320 m,南北范围约265 m。最终,开挖后地表全位移最大值为0.0108 m。
为进一步分析地下开采对地表的影响范围,沿矿体走向选取剖面进行分析。如Y=5 450剖面,见图5。由图5可以看出,地下矿体开挖后的影响范围将波及地表,主要集中在Y=4 500~4 900之间,位移量约为0.01 m。
图3 地表移动盆地边界示意图Fig.3 The schematic diagram of the boundary of the surface subsidence basin
图4 地表整体位移等值线图Fig.4 Contour map of surface whole displacement
图5 Y=5 450处整体位移等值线云图Fig.5 Contour map of whole displacement in Y=5 450
4.2.2 地表倾斜分析
地表倾斜见图6~7。在东西向倾斜方面,地下开采在地表形成了一正一负两个倾斜区域,以X=4 600~4 700为界,分别位于东西两侧。正倾斜、负倾斜区域的倾斜绝对值最大分别为4e-5m/m、3.5e-5m/m,均未达到3 mm/m,小于最低等级要求。在南北向倾斜方面,地表形成了一正一负两个倾斜区域,以Y=5 450~5 550为界,分别位于南北两侧。正倾斜、负倾斜区域的倾斜绝对值最大分别为4e-5m/m、4.5e-5m/m,均未达到3 mm/m。
图6 地表东西向倾斜Fig.6 Surface east-west tilt
图7 地表南北向倾斜Fig.7 Surface north-south tilt
4.2.3 地表曲率分析
地表曲率见图8~9。在东西向曲率方面,地表形成了一正一负两个曲率区域,正曲率、负曲率区域的曲率绝对值最大分别为4e-6/m、3.5e-6/m,均未达到2e-4/m。在南北向曲率方面,地表形成了一正一负两个曲率区域,正曲率、负曲率区域的曲率绝对值最大分别为1.2e-5/m、1.6e-5/m,均未达到2e-4/m。
图8 地表东西向曲率Fig.8 Surface east-west curvature
图9 地表南北向曲率Fig.9 Surface north-south curvature
4.2.4 地表水平变形分析
地表水平变形见图10~11。在东西向水平变形方面,地表形成了一正一负两个水平变形区域,以正变形区域居中,负变形区域位于其两侧,正变形、负变形区域的变形绝对值最大分别为3e-3mm/m、2.1e-2mm/m,均未达到2 mm/m。在南北向水平变形方面,地表形成了一正一负两个水平变形区域,正变形、负变形区域的变形绝对值最大为2.2e-2mm/m、3.6e-2mm/m,均未达到2 mm/m。
为了界定地下开采对地表的影响范围,根据国家的相关标准及规范对岩层移动形成的地表位移、倾斜、曲率及水平变形进行分析,并通过地表变形值的比较,最终圈定地表移动影响范围。
矿区934~1 054 m中段矿体开挖以后,采空区主要影响区域是矿体上盘,同时由于矿体沿东西走向较长,矿体在X=4 500~4 900之间比较集中,且有L27-8厚矿体(矿脉平均厚度达到4.68 m,最厚处可达16.3 m)存在,因此,在地表的开挖影响区域主要集中在X=4 500~4 900之间,这与地表整体位移及竖向位移分析结果一致,同时也与急倾斜矿体地下开采围岩的变形破坏规律[7,9]相吻合。因此,根据地表整体位移及地表竖向位移的分析,按照地表移动边界中最外边界(一般指下沉10 mm的点圈定的边界)的圈定原则,以下沉10 mm作为临界值圈定地表影响范围(图12),可见圈定的地表影响范围具有一定的可信度。
图10 地表东西向水平变形Fig.10 Surface east-west horizontal deformation
图11 地表南北向水平变形Fig.11 Surface north-south horizontal deformation
图12 地表影响范围Fig.12 The affected area of the surface
1) 借助AutoCAD、Dimine、Midas GTS NX等软件联合构建三维数值计算仿真模型,并通过转换软件生成FLAC3D文件,这一方法提高了建模效率,能够满足复杂地质模型的构建及三维计算分析的要求。
2) 地表倾斜、曲率及水平变形均小于砖混结构建筑物破坏等级要求的临界变形值,表明地表倾斜、曲率及水平变形均不会对地表建(构)筑物造成影响。
3) 地下采空区主要影响区域是矿体上盘,最终开挖影响区域主要集中在矿体分布集中区域对应的上盘位置,这与地表位移分析结果一致,同时也与急倾斜矿体地下开采围岩的变形破坏规律相吻合。
4) 以地表移动边界中最外边界(下沉10 mm为临界值)圈定地表影响范围具有一定的可信度。