, 马凯1,2, 赵永1,2 1. 东北大学 资源与土木工程学院, 辽宁 沈阳 110819;
2. 东北大学 深部金属矿山安全开采教育部重点实验室，辽宁 沈阳 110819
收稿日期：2022-05-13
基金项目：中央高校基本科研业务费专项资金资助项目(N2201003); 国家自然科学基金资助项目(U1710253, U1903216, 52004052)。
作者简介：侯俊旭(1996-), 男, 辽宁本溪人, 东北大学硕士研究生;
杨天鸿(1968-), 男, 辽宁抚顺人, 东北大学教授, 博士生导师。
摘要：为解决数值模拟方法在大型岩体工程稳定性分析中精度低、误差大的问题，开发了一种“亿级”自由度数值模拟技术.以大红山铁矿露天与地下联采地表岩移工程问题为例，详细阐述了“亿级”自由度数值模拟技术流程.将模拟计算结果与实测地表裂隙范围相结合，提出了一种模糊沉降值反演地表裂隙范围的方法.将该方法与“亿级”自由度数值模拟方法结合使用，可对地表裂隙发育范围进行准确估算，位置绝对误差在1.59~4.71 m，总长度误差率在1.141 %，验证了方法的可靠性，较常用的百万级自由度数值模拟方法有了提升.同时，对5年后露坑联采地表裂隙范围进行了估算.随着5年后地下300~720 m标高矿体采出，虽然塌陷坑废石充填对围岩变形具有抑制作用，但地表沉降及开裂范围仍会波及露天境界，预计裂隙区域至少将影响露天矿990 m标高台阶以上范围.
关键词：“亿级”自由度数值模拟精细化建模露坑联采地表沉降
More Than 100 Million DOF Numerical Simulation Technique and Its Engineering Application
HOU Jun-xu1,2, YANG Tian-hong1,2
, MA Kai1,2, ZHAO Yong1,2 1. School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China;
2. Key Laboratory of Ministry of Education on Safe Mining of Deep Metal Mines, Northeastern University, Shenyang 110819, China
Corresponding author: YANG Tian-hong, E-mail: yangtianhong@mail.neu.edu.cn.
Abstract: In order to solve the problems of low accuracy and abnormal error in numerical simulation method for stability analysis of large-scale rock mass engineering, a more than 100 million degree-of-freedom(DOF)numerical simulation technique was developed. Taking the surface rock movement engineering of open-pit and underground combined mining at Dahongshan iron mine as an example, this technique is elaborated in detail. Combining the simulated results with the fissure of ground surface range, a method of the fissure of ground surface range inversed by fuzzy surface subsidence value is proposed. Combined with the more than 100 million DOF numerical simulation technique, this method accurately estimates the development of the fissure of ground surface range. The absolute error of position is 1.59~4.71 m with a span error rate of 1.141 %, verifing the reliability of the method and improving accuracy compared to the commonly numerical simulation method of millions DOF. At the same time, the fissure of ground surface range is estimated after 5 years. With the mining of underground 300~720 m level ore body five years later, although the waste rock filling of collapse pit can inhibit the deformation of surrounding rock, the surface movement and cracking range will still affect the open-pit, and it is estimated that the crack area will affect the range above 990 m steps at least in the open-pit.
Key words: more than 100 million DOFnumerical simulationrefined modelingopen-pit and underground combined miningsurface subsidence
金属矿山开采与地质条件复杂，很多工程灾害的诱发与多种因素有关.例如金属矿地下崩落法开采造成的岩体冒落与地表沉降问题相比煤矿地下开采就复杂得多，这是因为矿体赋存、结构面、开采方法、岩层结构、地表地形、构造应力等多种因素的极大差异与相互影响造成的[1-4].目前，研究与防治矿山灾害的方法有很多，主要有工程类比法(经验法)、物理模拟试验、多手段现场监测、数值模拟等[5-6].工程类比法是参照类似工程的经验，并依据规范给出相应参考意见的方法，本质上是一种概念上的、定性的方法，没有严格的推理体系，不同工程类比结果存在较大差异，考虑因素也相对较少[7].物理模拟试验是基于相似理论，通过不同配比的材料模拟实际岩层与矿体，分步开挖模拟开采过程，埋设传感器与观测系统捕捉采动围岩物理场响应的方法[8-9].物理试验可直观反映采动效应，但仍无法完全映射真实的开采情况，这是由于尺度效应与复杂地质工程的关系[10].现场监测是通过现场布置一些传感器，如位移计、应力计、微震等，实时获取监测数据并动态分析，从而实现对地质灾害的跟踪监测[11-12].该方法源于实际，真实可靠，但监测数据为历史数据，且数据分析工作繁琐，不同长度时间线上的数据集成分析很难得到有指导意义的统一规律，甚至会出现难以解释的自矛盾现象，大幅降低灾害预测预警的可靠性.
随着计算机软、硬件技术的飞速发展，数值模拟方法在矿山开采中的应用越来越广泛，涵盖了岩体破坏、渗流突水、岩层移动等方面[13-15].数值模拟对复杂工程地质条件具有较高的适用性，可通过精细化建模考虑断层、弱面的影响，选用不同的本构方程可充分反映岩体材料的物理性质，采用有限元、离散元、边界元等合适的模拟方法可获得更加接近实际的分析结果.朱万成等提出了岩石细观单元的拉伸与剪切损伤本构方程，模拟与试验结果具有较高的一致性[16]；Nguyen等运用改进的有限元Monte-Carlo法模拟分析了相关矿山生产安全性[17]；其他****结合数值模拟方法对矿山安全等进行了详细论述[18-19].数值模拟方法虽然经济可靠，但对于一些大型的复杂工程问题，采用常规的数值模拟方法却很难达到符合相应精度标准的解答.首先，大型工程的尺度与范围较大，当要求的计算精度较高时，常规数值模拟方法可用于计算的网格与单元数目上限很难达到相应的精度要求，该条件下往往需要“亿级”自由度数值计算；其次，大范围工程内的地质条件复杂，要想得到精确的解答，就必须对影响计算精度的多尺度地质体进行精细化建模，普遍简化的常规数值模拟方法很难达到该要求；最重要的是，高精度解对计算能力的要求，常规数值模拟中与计算单元上限相匹配的算法与算力很难实现，需要使用高性能计算机并开发相应的并行计算方法.因此，开展“亿级”自由度数值模拟技术研究对分析大型工程稳定性具有重要意义.
本文提出了一种“亿级”自由度数值模拟技术，并以大红山铁矿露坑联采地表岩移工程问题为实际案例，对“亿级”自由度数值模拟技术流程进行详细阐述.同时，将模拟计算结果与实测地表裂隙范围相结合，提出了一种模糊沉降值反演地表裂隙范围的方法，验证了模拟计算及方法的可靠性.研究可为大红山铁矿露坑联采安全提供有益指导，亦可为“亿级”自由度数值模拟技术在大型岩土工程应用提供参考.
1 技术应用背景玉溪大红山铁矿位于云南省玉溪市新平彝族傣族自治县戛洒镇，现隶属于宝武钢铁集团，已探明的铁矿石储量为4.58亿t.矿区主要赋存有浅部与深部两条矿脉，分别采用露天开采与地下无底柱分段崩落法同时开采.地下矿体呈一轴向近乎东西向延伸的断块向斜产出，总体东高西低、中部厚边部薄、南北翘起、似船型，矿体东西长1 969 m，南北宽440~640 m，面积1.02 km2，埋深362.48~988.31 m，赋存标高为25.72~945.00 m.露采矿体赋存标高为680~1 012 m，走向北西、倾向南西，矿体平均厚度11~18 m，采用“陡帮剥离，缓帮采矿”工艺进行开采，露天采场上口尺寸长1 130 m，宽724 m，底部尺寸长523 m，宽172 m.优化设计后的露天开采生产能力为550万t/年，地下主采区生产能力为400万t/年，按标高划分4个区段同时开采，分别为头部、中部、下部与深部，除头部矿体采用空场法开采外，其余区段均采用无底柱分段崩落法开采，如图 1所示.
图 1(Fig. 1)
 | 图 1 大红山铁矿矿体分布与地下开采区段示意图Fig.1 Schematic diagram of ore-body distribution and underground mining section in Dahongshan iron mine |
目前，影响大红山铁矿安全高效生产的主要难题是露坑联采的相互影响问题.设计的露天最终境界距地下主采区的平面最短距离约250 m，随着地下主采区(下部与深部矿体)持续向深部扩展，以及靠露天边坡西北部矿体与矿柱的回收，覆岩冒落与地表变形范围的扩大将对露采边坡的稳定性产生重要影响.此外，深部矿体品位是浅部矿体的2倍，保障产能的前提下地下开采计划很难为露天开采做出让步，而浅部矿体现仍存留约3 000万t的矿体，露天停产将造成巨大的资源浪费与经济损失.可见大红山铁矿高效开采与安全生产的矛盾突出.
为解决大红山铁矿露坑联采问题，确保开采安全与产能保障，最大程度地回收矿产资源，就必须对未来5年间地下崩落法开采造成的岩移影响范围进行较为准确的预计.然而，露坑联采范围大，且矿区受大小断层交错切割，地质条件较为复杂，影响函数法、工程类比法、物理模拟等常规的岩移估算方法很难有效指导生产设计.鉴于数值模拟方法对复杂工程地质条件的独特优势，对地表岩移范围进行数值模拟分析不失为一种经济可靠的方法.然而，常规数值模拟在计算大范围大型工程时无法达到相应的精度要求，较大的计算偏差将会直接影响后续开采进度计划的设计，产生一定的致灾风险.因此，开展大红山铁矿露坑联采“亿级”自由度数值计算分析，更加精确地获取后续开采岩移发展范围与趋势，对指导后续安全生产设计，保障露坑联采平稳过渡具有重要的实际意义.本文以大红山铁矿露坑联采为实际案例，详述“亿级”自由度数值模拟分析的技术流程，为类似大型工程的精确计算提供有益借鉴.
2 “亿级”自由度数值模拟“亿级”自由度数值模拟技术主要包括：1)建立精细化的几何模型；2)高精度网格划分；3)单元加密；4)并行有限元数值计算.通过4个主要步骤，在高性能计算机上即可实现“亿级”自由度并行有限元数值计算，具体流程如图 2所示.
图 2(Fig. 2)
 | 图 2 “亿级”自由度数值模拟计算实现流程图Fig.2 More than 100 million DOF numerical simulation technique chart |
2.1 精细化几何模型建立几何模型的建立是最为关键的一步，几何模型的质量会直接影响后续高精度网格的划分.考虑技术软件间的兼容性与格式转换，几何模型建立选用功能强大的Rhino7.0软件，3DMine作为辅助处理工具.为减弱和消除露坑联采大范围模拟计算的边界效应，在露天与地下开采的投影平面上外推100 m作为几何模型的边界，采用1∶ 1的比例进行建模.为保障模型精度，绝对公差设置0.01，角度公差设置为1°.大量研究表明，地形地貌对地表沉降有较大影响[20-21].因此，计算模型必须能准确反映地表的起伏状态.矿山采用无人机倾斜摄影测量与激光扫描技术获取了高分辨率优于0.05 m的地表点云数据及模型，将其转化为dxf格式文件即可导入Rhino软件进行处理，采用布帘(网格等间距拟合曲面指令)等曲面工具处理后的高精度地表几何模型如图 3所示.
图 3(Fig. 3)
 | 图 3 大红山铁矿高精度地表几何模型Fig.3 High precision surface geometry model of Dahongshan iron mine |
提取图 3地表模型的边界线，可由边界线竖直向下正交延伸形成几何模型的侧面，再对底面进行加盖(即对底面空洞进行封闭以形成闭合的体模型)处理后，建成3 450 m×2 210 m×2 020 m的基础几何模型(即用来进行后续操作的包含地表的封闭几何体)如图 4所示.
图 4(Fig. 4)
 | 图 4 大红山铁矿高精度基础几何模型Fig.4 High precision basic geometric model of Dahongshan iron mine |
大红山铁矿主要岩性有辉长辉绿岩、变钠质熔岩、白云石大理岩、铜矿带、铁矿带，影响露坑联采的主要构造为F1与F2断层.F1断层走向北西西或近东西，倾向南或南西西，倾角60°~85°，断距大于500 m.F2断层早期为逆断层，晚期为正断层，走向近东西，延长1.1 km以上，倾向南，倾角80°左右.两断层带均由辉长辉绿岩侵入充填，F2断层是浅部与深部矿体的自然分解，F1断层是深部矿体的南界，如图 5所示.
图 5(Fig. 5)
 | 图 5 大红山铁矿岩性与地质构造分布示意图Fig.5 Schematic diagram of lithology and geological structure distribution in Dahongshan iron mine |
建立不同岩性与断层的几何模型可通过剖面与块体切割的方法实现.选取含有多个剖面线的特定标高平面作为放置不同剖面的参考面(此处选取180 m标高平面)，将带有岩性分界线的各剖面依次放置于参考面上(图 6)，参与本次建模的勘探剖面为A28~A45，共18个，相邻剖面水平间距50 m.参与建模的剖面越多，建立的几何模型越接近地质实际.
图 6(Fig. 6)
 | 图 6 参考面选取与剖面放置示意图Fig.6 Schematic diagram of reference surface selection and section placement |
选取每个剖面中相同位置同种岩性的轮廓线，使用放样功能由多条轮廓线形成对应的岩性外轮廓面.这里依据岩体与矿体的距离，以及不同岩体的尺寸差异分别选择不同的放样指标.距离开采矿体较远或空间形态简单的岩体选择“松弛(曲面拟合过程中根据控制点数量定义拟合逼近程度)”，例如矿区分布广泛的大理岩与变钠质熔岩.而对重点关注及空间形态复杂的区域则选择“紧绷(曲面拟合过程中根据控制点间距定义拟合逼近程度)”，例如F1与F2断层(辉长辉绿岩)、矿体，并通过控制逼近公差提高精确度.由此在保障可算性的同时不失模型精度，具体参数选择见表 1.
表 1(Table 1)
 表 1 不同尺寸与岩性建模放样标准参数Table 1 Lofting standard parameters for modeling different sizes and lithologies
| 表 1 不同尺寸与岩性建模放样标准参数 Table 1 Lofting standard parameters for modeling different sizes and lithologies |
使用已经建立的各岩体外轮廓面对建成的基础几何模型(图 4)进行切割，即可获得不同岩性岩体的几何模型.整个几何模型建立过程中，严格按照矿山精确勘查获取的地质剖面完成，实现了对大红山铁矿各个岩性、各个岩体的精准描述.
露天采场的建立主要依据露天最终境界平面图(CAD图)，将平面图导入矿业软件3DMine，对不同标高的台阶线(坡底线与坡顶线)赋值高程，形成露天采场边坡的大致轮廓.而后将具有高程的边坡轮廓线导入Rhino软件，借助“放样”、“扫掠”功能将轮廓线相连接，通过设置“放样”与“扫掠”控制点的个数以提高精度，参照模型尺寸大小(图 4)选取控制点个数为1 000，由此可生成整个露天坑的轮廓面.由轮廓面对已建成的含不同岩性的几何模型(图 7)进行切割，即可生成精细化的最终露天境界几何模型，如图 8所示.软件之间的兼容性保障了模型建立过程中坐标与尺寸的统一.
图 7(Fig. 7)
 | 图 7 精细化建模过程示意图Fig.7 Sketch diagram of fine modeling process |
图 8(Fig. 8)
 | 图 8 露天矿最终境界精细化几何模型图Fig.8 Fine geometric model drawing of the open-pit final boundary |
目前，大红山铁矿地下主采区造成的岩体冒落已发育至地表，且为限制岩移持续发育威胁露天边坡，利用露采剥离的废石对塌陷坑进行了回填.为考虑废石充填的影响，按照上述方法在几何模型中对冒落充填体进行单独建模，计算时便于赋值参数，如图 9所示.冒落区域的几何形态是根据矿山地表钻孔与巷道观测结果以及微震信息推测获取的，具体可参阅文献[22].
图 9(Fig. 9)
 | 图 9 冒落区精细化几何模型图Fig.9 Fine geometric model drawing of the caving zone |
2.2 高精度网格划分建立精细化几何模型后，需将其划分为三维体网格才能用于数值模拟计算.高精度网格划分在Hypermesh软件中进行.通过将几何模型导入，对模型进行实体几何修复、岩性共面连续性处理、边缘平滑处理，随后可进行二维面网格的划分.期间可通过调节网格大小与密度，提高网格质量.此外，Hypermesh软件提供了网格质量检查、畸变网格修补、自由边处理等功能，保障了网格的可算性与计算精度，网格处理前后对比如图 10所示.最后，统一面网格法向量、进行三维体网格划分(四面体网格)，在Hypermesh软件中对不同岩体进行分组，可用于数值模拟软件中不同岩体的赋值参数.至此，可完成高精度网格的划分(图 11)，依据选用的数值模拟软件类型可对应导出不同类型的网格信息文件.
图 10(Fig. 10)
 | 图 10 处理前后的高精度网格Fig.10 High precision mesh before and after processing (a)—处理前；(b)—处理后. |
图 11(Fig. 11)
 | 图 11 高精度三维体网格(四面体网格)图Fig.11 High precision 3D volume mesh(tetrahedral mesh)diagram |
2.3 单元加密Hypermesh软件完成的高精度网格数量无法达到“亿级”自由度标准，需进行单元加密处理.单元加密在有限元数值模拟软件COMSOL Multiphysics中完成.Hypermesh软件中网格导出选择Nastran格式可与COMSOL Multiphysics兼容.通过COMSOL模型开发器中的网格细化工具进行多次细化，即可实现单元加密.最终密化后的单元数量达2亿，自由度可达到“亿级”水平.由模型尺寸估算四面体单元尺寸为5.49 m，重点区域单元尺寸仅为2.46 m.加密处理后的网格模型贴合大红山铁矿露坑联采实际，可大幅提高数值计算精度，获得更加准确的地表影响范围，如图 12所示.
图 12(Fig. 12)
 | 图 12 密化处理后的网格模型图Fig.12 The grid model diagram after densification treatment |
2.4 并行有限元计算“亿级”自由度数值计算对软件及硬件条件要求极高，除配置性能计算机外，需采用并行计算方法.COMSOL Multiphysics经过几代升级，已集成两种并行计算方法，分别为单机多核与多机多核计算.单机多核是将求解域的刚度矩阵划分为若干个子块，让同一CPU的多个核通过数据分布存储方式进行并行代数运算，通过选择专门的优化求解器进行迭代求解计算.单机多核并行计算模式对计算机性能要求很高，尤其是存储容量，因为初始输入的数值模型仅靠一个CPU读取，读取效率较慢，占用内存大.这种计算机价格昂贵，而且主板必须满足高内存兼容性和稳定性.多机多核并行计算采用多核多节点的方法求解子域.并行计算系统结构复杂，从网格生成、区域分解到并行求解都需要一个完整的集成方案，且求解速度通常与CPU节点之间的通信速度有关.
本次模拟基于单机多核的并行计算，计算机为高性能的华硕服务器，其配置为酷睿至强W3245M CPU(16核32线程)、RTX2080Super显卡、512 GB内存；系统环境为64位Win10专业版；并行环境选用COMSOL Multiphysics自带的并行库.经初步试算，可满足“亿级”自由度有限元数值模拟分析.
3 模拟计算分析3.1 模拟方案模拟计算在有限元数值软件COMSOL Multiphysics中进行.模型底面采用固定约束，侧面法向位移为0，整个模型全局在重力作用下进行弹性求解.通过现场取样与室内试验，并由Hoek-Brown方法获取相应的岩体力学参数[23]，如表 2所示.
表 2(Table 2)
表 2 岩体力学参数Table 2 Mechanical parameters of rock mass
| 表 2 岩体力学参数 Table 2 Mechanical parameters of rock mass |
基于大红山铁矿露坑联采实际工况，本次模拟分析可分为以下几个工况进行，如表 3所示.
表 3(Table 3)
表 3 工况设置与说明Table 3 Working condition setting and description
| 表 3 工况设置与说明 Table 3 Working condition setting and description |
为便于分析露坑联采对地表变形产生的影响，以实际勘探剖面A36对应的地表投影线为监测线，对比分析不同工况条件下的地表变形情况.图 13为A36剖面岩性情况，图 14为数值模型中由A36剖面确定的地表监测线.
图 13(Fig. 13)
 | 图 13 A36剖面岩性分布图Fig.13 Lithologic distribution map at section A36 |
图 14(Fig. 14)
 | 图 14 监测线位置示意图Fig.14 Schematic diagram of monitoring line position |
A36监测线总长度2 687.7 m，本次模拟计算使用每间隔1 m设置一个监测点的方式对监测线进行点拟合.自监测线起点位置，依次为露天坑(838.86~1 510.73 m)，F2断层(1 704.59~1 794.23 m)，地下主矿体上方岩层(1 794.23~2 235.87 m)，F1断层(2 235.87~2 479.59 m).
3.2 结果分析3.2.1 露坑联采地表移动分析以工况2为基础，计算工况3，4的z向位移差，分别见图 15，图 16.可以看出，地下420~720 m主矿体开采过后，其产生的地表沉降范围基本波及至露天最终境界，露天境界附近处的沉降值达到了0.1~0.2 m，最大沉降值出现在塌陷区附近，在0.5 m以上.地表沉降形态与地下矿体开采区较为相似，随着400 m分段持续采出，虽然对地表塌陷坑进行了废石充填，对围岩变形会产生一定的抑制作用[8]，但地表沉降仍在扩展，露天境界附近的沉降值已达到0.2~0.3 m.
图 15(Fig. 15)
 | 图 15 冒落区域未充填模拟结果图(工况3)Fig.15 Simulation result diagram of caving zone unfilled(working condition 3) (a)—冒落区域未充填模拟俯视图；(b)—冒落区域未充填模拟立体图. |
图 16(Fig. 16)
 | 图 16 冒落区域废石充填模拟结果图(工况4)Fig.16 Simulation result diagram of waste rock filling in caving zone(working condition 4) (a)—冒落区域废石充填模拟俯视图；(b)—冒落区域废石充填模拟立体图. |
由现场实际勘察GPS定位结果，2017年4月(未充填前8个月)与2018年4月(充填后4个月)矿山地表裂隙范围如图 17所示.工况3(未充填)与工况4(充填)条件下，数值模拟A36监测线上的沉降值分布如图 18所示.由图 17可以看出，充填后，随着地下持续向深部开采，地表裂隙向露天坑方向仍有小幅扩展，这与数值模拟获取的地表沉降结果一致.在图 18中的地表沉降曲线可以看出，虽然F1与F2断层岩性相同，但由于F1断层厚度明显大于F2，沉降曲线的平均斜率也明显较小.同时可以看出，在露天范围内，沉降曲线的曲率明显减小.
图 17(Fig. 17)
 | 图 17 冒落区域废石充填前后地表裂隙发育实测结果图Fig.17 Measured results of surface fissure development before and after waste rock filling in caving zone |
图 18(Fig. 18)
 | 图 18 冒落区域废石充填前后A36监测线地表沉降分布曲线图Fig.18 Surface subsidence distribution curves of monitoring line A36 before and after waste rock filling in caving zone |
现场监测仅统计了地表裂隙的发育范围，未获取系统的地表位移数据，但已有研究表明，地表裂隙发育与沉降位移具有密切的伴生关系[26].为此，在高精度数值模拟计算下，可借助z向位移值对地表裂隙范围进行反分析.由于在横向与竖向上，岩层厚度、围岩物理参数等的差异性，不能用统一的固定z向位移值确定裂缝产生的范围.这里使用非固定的z向位移值来确定地表裂缝范围(即模糊的z向位移值).选取2017年地表裂缝范围在A36测线两交点处对应的z向位移值作为标准，分别为-0.219 59 m和-0.246 54 m，由此反演获取的工况4地表裂隙范围在A36测线上的位置见图 19，总长度X=572 m.
图 19(Fig. 19)
 | 图 19 由沉降值非固定标准反演废石充填后地表裂隙范围图Fig.19 The surface fissure range map after waste rock filling inverted from the surface subsidence floating standard |
结合实际测得的2018年4月(充填后4个月)GPS监测裂隙范围在A36测线上的位置，以及总长度L=578.3 m，由式(1)进行误差分析.分析结果显示：位置绝对误差在1.59~4.71 m，总长度绝对误差为6.3 m，且误差率仅为1.141 %.这不仅对沉降与裂隙发育的伴生关系进行了验证，同时也表明“亿级”自由度数值计算与z向位移值标准反演地表裂隙范围方法的可靠性.
 | (1) |
同样的分析过程在使用常用的百万级自由度(320万)进行分析计算的过程中，得到位置绝对误差在24.32~35.21 m，总长度绝对误差为59.53 m，误差率为10.312 %.较之“亿级”自由度数值模拟计算的结果放大了一个量级，这也证明该技术较十万级自由度有了较大的提升.在实际工程应用中，可有效提高预测范围的准确度，为安全开采提供可靠保障.
3.2.2 露坑联采地表沉降预测数值模拟得到工况3，4，5的A36测线地表沉降分布汇总如图 20所示.可以看出5年后300~720 m标高矿体采出后，露天坑范围内的沉降分布曲线的曲率有所增大，曲率增大易在地表产生张拉裂缝.
图 20(Fig. 20)
 | 图 20 数值模拟得到工况3，4，5的A36测线地表沉降分布曲线图Fig.20 The surface subsidence distribution curves of monitoring line in A36 under conditions 3, 4 and 5 obtained by numerical simulation |
使用的实测数据越多，预测得到的地表裂隙范围越准确.为此，分别基于2017年与2018年的实测地表裂隙范围，采用z向位移值非固定指标对工况5的地表裂隙范围进行估算，估算的地表开裂位置及总长度分别如图 21与图 22所示.综合两组模糊沉降值的数值模拟预测结果，预测5年后地表裂隙范围在A36测线上的位置为：露天坑一侧起裂点在距离测线起点1 364~1 380 m处，另一侧起裂点在距离测线起点2 583~2 585 m处.裂隙总长度为1 203~1 221 m.据此，预计裂隙区域至少将影响露天矿990 m标高台阶及以上范围.
图 21(Fig. 21)
 | 图 21 基于2017年地表裂隙范围估算工况5(5年后)地表裂隙范围图Fig.21 The surface fracture range map of working condition 5(5 years later) estimated based on the surface fissure range in 2017 |
图 22(Fig. 22)
 | 图 22 基于2018年地表裂隙范围估算工况5(5年后)地表裂隙范围图Fig.22 The surface fracture range map of working condition 5(5 years later) estimated based on the surface fissure range in 2018 |
4 结论1) 针对数值模拟方法在大型岩体工程稳定性分析中精度低、误差大的问题，开发了一种“亿级”自由度数值模拟技术，并以大红山铁矿露坑联采地表岩移工程问题为例，详细阐述了“亿级”自由度数值模拟技术流程.
2) 将模拟结果与实测地表裂隙范围相结合，提出了一种模糊沉降值反演地表裂隙范围的方法.将该方法与“亿级”自由度数值模拟方法结合使用，可对地表裂隙发育范围进行准确估算，位置绝对误差在1.59~4.71 m，总长度误差率在1.141 %，验证了模拟结果与方法的可靠性.
3) 对大红山铁矿5年后露坑联采地表裂隙范围进行了估算，地表沉降及开裂范围会波及露天境界，预计裂隙区域至少将影响露天矿990 m标高台阶以上范围，可为大红山铁矿露坑联采安全提供有益指导.
参考文献
| [1] | Woo K S, Eberhardt E, Elmo D, et al. Empirical investigation and characterization of surface subsidence related to block cave mining[J]. International Journal of Rock Mechanics and Mining Sciences, 2013, 61: 31-42. DOI:10.1016/j.ijrmms.2013.01.015 |
| [2] | Yao L, Xie G Q, Mao J W, et al. The relationship between granitic rock and skarn Fe deposit: a case study from the Chengchao Fe deposit, Edong ore district, middle-lower Yangtze River Valley Metallogenic Belt, Eastern China[J]. Acta Geologica Sinica (English Edition), 2014, 88(sup2): 56-57. |
| [3] | Tan Y Y, Guo M C, Hao Y M, et al. Structural parameter optimization for large spacing sublevel caving in Chengchao Iron Mine[J]. Metals, 2021, 11(10): 1619. DOI:10.3390/met11101619 |
| [4] | Wang Z W, Song G F, Ding K. Study on the ground movement in an open-pit mine in the case of combined surface and underground mining[J]. Advances in Materials Science and Engineering, 2020, 2020: 8728653. |
| [5] | Sheorey P R, Loui J P, Singh K B. Ground subsidence observations and a modified influence function method for complete subsidence prediction[J]. International Journal of Rock Mechanics and Mining Sciences, 2000, 37(5): 801-818. DOI:10.1016/S1365-1609(00)00023-X |
| [6] | Christina D L, Erling N, Tomas V. Hanging wall surface subsidence at the Kiirunavaara Mine, Sweden[J]. Engineering Geology, 2011, 121(1/2): 18-27. |
| [7] | Jin S L, Zheng T J, Lei B, et al. Design of underground engineering analogy based on fuzzy comprehensive evaluation model[J]. Applied Mechanics and Materials, 2012, 1800: 1726-1730. |
| [8] | Ren F Y, Zhou Y J, He R X, et al. Similarity model test on the spatiotemporal evolution law of deformation and failure of surrounding rock-induced caving in multi-mined-out areas[J]. Advances in Civil Engineering, 2021, 2021: 1224658. |
| [9] | Ghabraie B, Ren G, Zhang X Y, et al. Physical modelling of subsidence from sequential extraction of partially overlapping longwall panels and study of substrata movement characteristics[J]. International Journal of Coal Geology, 2015, 140: 71-83. DOI:10.1016/j.coal.2015.01.004 |
| [10] | 赵阳升. 岩体力学发展的一些回顾与若干未解之百年问题[J]. 岩石力学与工程学报, 2021, 40(7): 1297-1336. (Zhao Yang-sheng. Retrospection on the development of rock mass mechanics and the summary of some unsolved centennial problems[J]. Chinese Journal of Rock Mechanics and Engineering, 2021, 40(7): 1297-1336.) |
| [11] | Stevens A W, Moritz H R, Elias E P L, et al. Monitoring and modeling dispersal of a submerged nearshore berm at the mouth of the Columbia River, USA[J]. Coastal Engineering, 2023, 181: 104285. DOI:10.1016/j.coastaleng.2023.104285 |
| [12] | 杨天鸿, 王赫, 董鑫, 等. 露天矿边坡稳定性智能评价研究现状、存在问题及对策[J]. 煤炭学报, 2020, 45(6): 2277-2295. (Yang Tian-hong, Wang He, Dong Xin, et al. Current situation, problems and countermeasures of intelligent evaluation of slope stability in open pit[J]. Journal of China Coal Society, 2020, 45(6): 2277-2295.) |
| [13] | 张超, 宋卫东, 付建新, 等. 深部复杂金矿体充填开采对地表建筑物稳定性影响[J]. 东北大学学报(自然科学版), 2021, 42(8): 1143-1151, 1158. (Zhang Chao, Song Wei-dong, Fu Jian-xin, et al. Impact on the stability of surface construction by filling method of deep complex gold ore body[J]. Journal of Northeastern University (Natural Science), 2021, 42(8): 1143-1151, 1158.) |
| [14] | Yang W M, Geng Y, Zhou Z Q, et al. True triaxial hydraulic fracturing test and numerical simulation of limestone[J]. Journal of Central South University, 2020, 27(10): 3026-3039. |
| [15] | Xu J M, Zhu W B, Xu J L, et al. High-intensity longwall mining-induced ground subsidence in Shendong coalfield, China[J]. International Journal of Rock Mechanics & Mining Sciences, 2021, 141: 104730. |
| [16] | 朱万成, 唐春安, 杨天鸿, 等. 岩石破裂过程分析用RFPA-2D系统的细观单元本构关系及验证[J]. 岩石力学与工程学报, 2003, 22(1): 24-29. (Zhu Wan-cheng, Tang Chun-an, Yang Tian-hong, et al. Constitutive relationship of mesocopic elements used in RFPA-2D and its validations[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(1): 24-29.) |
| [17] | Nguyen V U, Chowdhury R N. Simulation for risk analysis with correlated variables[J]. Géotechnique, 1985, 35(1): 47-58. DOI:10.1680/geot.1985.35.1.47 |
| [18] | Cheng G W, Yang T H, Liu H Y, et al. Characteristics of stratum movement induced by downward longwall mining activities in middle-distance multi-seam[J]. International Journal of Rock Mechanics and Mining Sciences, 2020, 136: 104517. |
| [19] | 王连庆, 高谦, 王建国, 等. 自然崩落采矿法的颗粒流数值模拟[J]. 北京科技大学学报, 2007, 29(6): 557-561. (Wang Lian-qing, Gao Qian, Wang Jian-guo, et al. Numerical simulation of natural caving method based on particle flow code in two dimensions[J]. Journal of University of Science and Technology Beijing, 2007, 29(6): 557-561.) |
| [20] | Benko B. Numerical modelling of complex slope deformations[D]. Saskatoon: University of Saskatchewan, 1997. |
| [21] | Zhou Y, Feng S W, Li J W. Study on the failure mechanism of rock mass around a mined-out area above a highway tunnel—similarity model test and numerical analysis[J]. Tunnelling and Underground Space Technology Incorporating Trenchless Technology Research, 2021, 118: 104182. |
| [22] | 李庶林, 胡静云, 周爱民, 等. 崩落法上覆厚大岩层崩落及破裂特性综合研究[J]. 岩石力学与工程学报, 2016, 35(9): 1729-1739. (Li Shu-lin, Hu Jing-yun, Zhou Ai-min, et al. Comprehensive research on character of collapse and fracture of thick and large overburden rock in cave mining[J]. Chinese Journal of Rock Mechanics and Engineering, 2016, 35(9): 1729-1739.) |
| [23] | Kang K, Huang S Y, Liu W H, et al. Sandstone slope stability analysis under wetting-drying cycles based on generalized Hoek-Brown failure criterion[J]. Frontiers in Earth Science, 2022, 2022(10): 838862. |
| [24] | Bai L Y, Liao C L, Wang C X, et al. Study on height prediction of water flowing fractured zone in deep mines based on Weka platform[J]. Sustainability, 2022, 15(1): 737. |
| [25] | Kuznetsov G N, Bud'ko M N, Vasil'ev Y N, et al. Models of rock pressure phenomena[J]. Soviet Mining, 1970, 6(2): 239-240. |
| [26] | 孙月敏, 杨天亮, 卢全中, 等. 基于SBAS-InSAR的西安市鱼化寨地区地面沉降与地裂缝时空演变特征研究[J]. 工程地质学报, 2022, 30(2): 553-564. (Sun Yue-min, Yang Tian-liang, Lu Quan-zhong, et al. Spatial-temporal evolution characteristics of land subsidence and ground fissure in Yuhua Village of Xi'an City using SBAS-InSAR technique[J]. Journal of Engineering Geology, 2022, 30(2): 553-564.) |
