基于Rhino-Griddle-Flac3D的矿体建模与地表稳定性分析

2020-09-19石满生张美道

金属矿山 2020年8期

石满生许威张美道

（1.江西天億矿业有限公司，江西贵溪 335400；2.江西理工大学资源与环境工程学院，江西赣州 341000）

矿体地下开采后，周围岩层的地应力会重新分布。随着采空区暴露面积的增大，若不采取有效措施，开采扰动将会对自身开采范围内（地表移动带）的地表建/构筑物以及相临矿山的安全生产产生一定影响［1］，准确预测和有效预防的技术手段就是开采前进行三维模拟。在运用三维数值模拟计算预测方面，国内外学者进行了较多研究。郑文棠等［2］基于AutoCAD平台，借助AutoLisp语言，采用滑动最小二乘法插值拟合，构建三维可视化模型，揭示边坡的宏观破坏模式。杨美宏等［3］采用3Dmine软件对地表及井下工程建立几何模型，运用Flac3D模拟分析采空区对地表稳定性的影响。李占金等［4］基于FLAC3D对马城铁矿深部大采场开采及回填过程中，围岩与充填体的稳定性进行了数值模拟计算。陈庆发等［5］以3DMine数字化模型为基础，提出了3DMine-FLAC3D耦合建模方法和3DMine-Surfer-Rhino-ANSYS-FLAC3D多软件耦合建模方法，构建了锌多金属矿三维数值模型，分析了矿体上行开采地表沉陷规律。余璨等［6］基于DIMINE软件，创建铜厂矿体沿走向、倾向、厚度3个方向上的矿化数学模型，得出Cu品位空间分布模型。程海勇等［7］应用FLAC3D数值分析软件，探讨了龙首矿进路式采矿及充填对地表下沉的影响规律。程立年等［8］采用FLAC3D数值模拟方法对该矿“三下”开采移动范围内典型剖面的地表沉降进行了分析，得出该矿“三下”开采移动带内地表建（构）筑物安全稳定。Nengxiong Xu等［9］通过有限差分法（FDM）进行了开采引起的地面沉降预测，以判断煤层的开采是否会对大坝产生负面影响。Lü Xiangfeng等［10］使用3D有限元代码对长壁板的地层行为进行了3D建模，并使用CFD代码对采空区瓦斯气流进行了模拟，该研究对采矿引起的地层应力变化，裂缝和气体流型之间的复杂动力相互作用产生了许多新见解。

为确定天億矿业开采对地表稳定性以及银海矿业的影响，采用“Rhino-Griddle-Flac3D”联合方法，建立矿体与围岩的三维数值计算模型，通过计算4种工况和3个变形参数予以分析。

1 工程概况

天億矿业是一座基建期地下矿山，以层状斑岩型矿床为主，西南方向与银海矿业、鲍家矿业毗邻（如图1）。其中，鲍家矿体处于浅部，南矿带位于标高+50 m以上、北矿带位于标高+80 m以上，且通地表。银海矿业部分矿权与鲍家南矿带在垂直方向交接，标高+50 m以上为鲍家南矿带、+50～-70 m为隔离岩体、-70～-420 m为银海矿权。银海矿业为生产矿山，当前设计排采计划至2041年；鲍家矿业也处于生产阶段，但资源即将枯竭，故不予考虑。

根据矿山提供的勘探线剖面图、地表地形图，地质勘查报告等资料，应用“Rhino-Griddle-Flac3D”三维数值建模与计算方法，建立天億矿业、银海矿业、鲍家矿业的矿体与围岩计算模型，探究天億矿业开采对自身开采范围地表稳定性以及相临银海矿业的影响。

2 建立模型

根据圣维南原理，将计算范围设为（X、Y、Z）从（39 518 500，3 087 800，-850）至（39 521 000，3 090 600，地表），且设（39 518 500，3 087 800，-850）为模型计算相对原点（0，0，0），X轴朝东为正、模型长度2 500 m，Y轴朝北为正、模型宽度2 800 m，Z轴铅直朝上为正、模型高度从-850 m标高到地表。

2.1 建立矿体模型

将银海矿业（144～120#勘探线）、鲍家矿业（130～100#勘探线）、天億矿业（100～187#勘探线）勘探线剖面图导入Rhino，再通过地表XY轴坐标系、标高等作为参考，使用绕轴旋转、平移等命令调整好各个勘探线剖面图的位置。在Rhino中，使用“控制点曲线”命令将矿体轮廓圈定出来再使用放样、曲线、曲面、组合等命令建立矿体相邻勘探线剖面之间的同一矿体模型，确保每一个矿体模型都为“封闭的实体—多重曲面”。建立的矿体模型如图2所示。

2.2 建立围岩模型

将地质地形图导入Rhino，提取等高线。使用“曲线分段”命令将等高线每隔2m划分成高程点（图3（a）所示），选中所有高程点后使用“嵌面”命令生成地表模型。再用通过计算范围的矩形线框对地表模型进行布尔分割，得到计算范围内的地表—围岩模型。

计算范围内的围岩主要为上盘花岗斑岩和下盘晶屑凝灰岩。将相邻的勘探线剖面图中的地质体分层线使用“线—面”命令生成地质界面。通过地质界面，使用布尔分割命令将计算范围内的地表—围岩模型划分成2种岩性（图3（b）所示）。

2.3 基于Griddle插件的网格导出

通过Griddle插件将Rhino中网格导入到Flac3D。该插件与陈金龙等［11］应用KUBRIX插件构建大型复杂地质体相比具有以下3个优势：①Griddle插件中的GInt功能能够自动生成交互式网格，清理不正确相交的曲面网格，减少因为网格自交而导致导出时的报错问题（图4（a）所示）；②如果生成的网格质量差或者未闭合等情况，Griddle插件能够显示出出错的位置，能够很快对模型进行调整，极大地提高了检查模型的效率（图4（b）所示）；③Griddle插件中的Gsurf命令用于将所选的表面网格重新划分为所需的单元大小和类型，重新划分后的曲面网格作为输入数据到Gvol功能中进行体积剖分，使用表面网格来确定实体单元的大小和类型（图4（c）所示，矿体与围岩网格长度不同），Gvol命令能填充四面体或六面单元（六面体、棱镜），生成用于如Flac3D或3DEC的数值模拟程序的文件格式，该功能满足于本次研究中所需要的“矿体—围岩网格不同单元大小”情况。通过“Rhino-Griddle-Flac3D”数值建模方法，将导出的GRID文件导入到Flac3D中（图4d所示）。

注：银海为黄色模型，鲍家为粉色模型，天億为绿色模型

3 数值模拟计算与结果分析

3.1 边界条件

假设模型中矿岩体均为理想弹塑性连续介质，且采用摩尔—库仑（Mohr-Coulomb）屈服准则，其物理力学模型采用弹塑性力学模型，仅考虑自重应力。数值模型初始边界条件为：前后左右面及底面采用位移约束，即前后左右面限制水平方向位移，底面限制垂直方向位移［12］。

（1）X轴边界：约束X方向的移动，其边界位移为0。

（2）Y轴边界：约束Y方向的移动，其边界位移为0。

（3）Z轴边界：约束Z方向下部边界位移为0，上部设为自由边界。

选用折减后的工程岩体物理力学参数如表1。

3.2 工况设计

为探究天億矿业开采对计算范围内地表稳定性及银海矿业的影响，根据天億矿业和银海矿业排采计划，选择以下4种代表性工况。工况1（最不利情况）：天億矿业和银海矿业全部矿体开采，不充填；工况2（最有利情况）：天億矿业和银海矿业全部矿体开采，嗣后全尾砂胶结充填；工况3（天億矿业开采至+10 m，银海矿业开采至-240 m）：天億矿业和银海矿业开采到2035年，嗣后全尾砂胶结充填；工况4（天億矿业开采至-40 m，银海矿业开采结束）：天億矿业和银海矿业开采到2041年，嗣后全尾砂胶结充填。由于鲍家矿业即将闭坑且矿体通地表，数值模拟中只开采不充填。

在地表40个重点位置设置变形监测点（图5），具体为：村庄CZ1～CZ3、冷水河HL1～HL4、公路GL1～GL5、银海工业场地和隔离矿柱YH1～YH17、银珠山主副井和南北风井等工业场地GC1～GC9、银珠山斜坡道平硐PD1～PD2。

3.3 计算结果分析

在Flac3D中，根据所设计4种工况进行数值模拟计算。地表沉降情况见图6，地表变形量汇总见表2。

由表2可知，工况1～工况4的地表40个监测点的最大倾斜i、最大曲率k及最大水平变形ε均小于《有色金属采矿设计规范》I级保护物、《煤矿测量规程》（2013版）一般砖石结构建筑物规定i=3 mm/m、k=0.2 mm/m/m、ε=2 mm/m的临界变形极限值。说明天億矿业在工况1～工况4条件下开采，对地表无显著不利影响。且对比工况1和工况2发现，充填后天億矿业最大沉降减小了316.68 mm、银海矿业最大沉降减小了65.11 mm，监测点最大倾斜减小了1.138 9 mm/m、最大曲率减小了0.015 8 mm/m/m、最大水平变形减小了0.217 4 mm/m，说明采空区嗣后充填能显著降低地表变形，提高地表稳定性。