大孤山露天矿边坡亿级自由度建模与基于RFPA3D的数值模拟

李斓堃 朱万成 代 风 刘溪鸽 邓文学

（东北大学资源与土木工程学院，辽宁 沈阳 110819）

近年来，数值模拟被广泛应用于采矿工程中的岩体稳定性分析。大型岩土工程计算区域常达千米级，并且地质构造复杂、岩石类型多样以及涉及到多场耦合等问题。Park D等［1］建立了长达1 000 m的隧道模型，并得到隧道沿轴线方向的变形响应特征；尚振华等［2］采用FLAC3D软件对尺寸为1 000 m×1 600 m×660 m的某钨矿模型进行模拟，将采空区稳定性问题量化；刘春等［3］建立了尺寸为2 800 m×1 500 m×1 600 m的茂县新磨村滑坡模型，模拟了新磨村滑坡启动、高速下滑和堆积的全过程；张成良等［4］采用三维有限元模拟了3种开采方案对地表公路的影响，得到当开采高度超过1 300 m水平时，公路上部围岩将处于不稳定状态；漆祖芳等［5］建立630 m×1 850 m×900 m的有限元模型，模拟了大岗山水电站坝肩边坡开挖支护过程；马克等［6］用RFPA3D对大岗山水电站坝址区右岸实际微震监测范围（690 m×400 m×1 600 m）等比例的三维模型进行模拟，对抗剪洞加固前、后的边坡稳定性进行分析。

考虑到计算成本及技术方面的问题，前人对计算模型进行了大量简化，且受到计算模型网格数量受限，致使模型网格尺寸较大，导致了数值模型难以反映实际工程的现场真实情况［7］。因此，对于大规模、复杂岩石工程的数值模拟，大多数情况下会将关键位置网格和非关键位置网格差异化对待，即关键位置网格剖分需要细化，而非关键位置可采用过渡单元，将单元尺寸从核心研究区域向边界区域逐渐粗化［8-9］。网格单元尺寸的差异化将客观地对岩石工程内部应力的传播、损伤演化及裂纹扩展等路径造成不可忽视的影响。因此，在不考虑计算机性能和计算成本的前提下，数值模拟计算的期望是将计算单元或网格合理地细化，以期更为准确地表征围岩特征。

在数值模拟计算中，数值结果依赖于网格尺寸，如果模型中相互接触的网格尺寸比值过大，则容易造成数值模拟结果的失真。因此，诸多研究和实践表明，数值计算模型中单元体积的变化率应该控制在30%以内［10］。此外，网格的尺度还应该小于材料变形的最小特征尺寸，以满足细致刻画岩体变形和损伤演变过程的要求。基于上述两点，对于含有多种岩体的大型复杂岩石工程而言，既要满足网格尺寸小于岩石变形的最小特征尺寸，又要满足单元体积变化率控制在30%的要求，这就需要建立尺寸足够小、数量足够多的网格和单元，其自由度往往达到千万级甚至亿级。

随着大型工程对精细数值分析要求的不断增加，超过千万自由度的数值分析不断增多［11］，然而，受到计算机性能，如CPU浮点计算能力、GPU图像处理能力以及内存开销等硬件因素，采用常规数值模拟软件计算一些大的模型需要几天甚至几个月的时间［12］，这样就极大提高了工程分析的时间成本，也使得分析人员不得不减小计算规模。为此，一些研究人员利用高性能并行计算机对复杂岩土工程问题进行模拟。范宣华等［13］构建了大规模模态分析并行计算体系，并验证3种算法构建的并行求解体系均可在1 h内求解千万级自由度量的大规模模态分析问题。张友良等［14］采用网格加密方法生成5 000万自由度的有限元模型，利用对偶原始有限元撕裂内联法（FETI-DP）进行求解。K.Garatani等［15］采用GeoFEM系统对简单形状进行线弹性分析，该计算自由度高达一亿。刘耀儒等［16-17］考虑到耦合分析时步多、计算量大的问题，采用基于element-by-element策略的有限元并行计算方法进行数值模拟。针对岩石破裂及大型岩体失稳的求解问题，张永彬等［18］结合现代有限元方法和数值计算方法，在消息传递并行环境下，利用区域分解和主从编程模式，建立了岩石破裂过程分析RFPA3D-Parallel并行分析系统。

本研究提出了新的结构化六面体网格化方法，使用东北大学AMAX PSC-HC2X GPU服务器，建立了大孤山露天矿亿级自由度三维模型，并基于RFPA3D-Parallel软件平台对该模型进行了数值分析，以期为边坡稳定性评价及灾害防控提供参考。

1 RFPA3D-Parallel精细化数值模拟方法

1.1 RFPA3D-Parallel简介

RFPA3D-Parallel并行分析系统采用分布存储稀疏线性迭代并行求解方法，在Linux机群上实现应力分析模块中有限元计算的并行处理，通过Windows和Linux协调处理策略，有效地把原有的前后处理功能和机群系统强大的计算能力结合起来，建立岩石破裂过程分析平台［19］。经过后期实验验证，RFPA3DParallel在确保高效、准确、强大的计算能力的同时，实现了岩石介质破坏过程的准确还原［20-21］。

为了在数值模拟中表达天然岩石的非均质特征，假设岩石由许多细观单元组成，其力学性质服从Weibull分布［22］。利用Weibull分布确定岩石参数的方法已被大量应用于数值模拟，并取得了很好的结果［23-26］。

岩石的损伤与岩石的应力状态密切相关。RFPA3D中采用简单的弹性损伤本构模型，在达到破坏准则之前，单元保持线弹性的力学本质。当岩石的应力状态满足最大拉应力准则或摩尔—库仑准则时，岩石就开始发生破坏［27］，其力学模型为

式（1）为最大拉应力准则；式（2）为经典的摩尔库伦准则；式（3）中ω为损伤变量，E、E0分别为损伤单元和未损伤单元的弹性模量。

1.2 网格划分与模型构建

为快速、自动化实现复杂岩体六面体网格划分，并提供RFPA3D-Parallel可直接读取的网格文件，本项目开发了一种复杂岩体亿级自由度结构化六面体网格参数化建模方法。基于亿级自由度的结构化六面体网格参数化建模方法包含了背景网格划分、四面体网格数量估算、六面体网格更新、空属性六面体网格检索和剔除、网格文件导入等过程。基于RFPA3DParallel数值模拟软件，针对大规模复杂岩体精细化数值计算的技术路线如图1所示。

（1）非结构化背景网格划分。首先对大型不规则结构的三维复杂地质体模型进行万级自由度的背景网格划分，并赋予编号IDt，t从1开始。该网格划分过程采用非结构化网格，它针对复杂地质体的不同位置采用自适应过渡的四面体网格自动划分方法。

（2）结构化六面体网格建立。根据非结构化的四面体背景网格信息，确定x、y、z3个方向上的坐标最大值和最小值xmax、ymax、zmax、xmin、ymin、zmin。从而确定包含整个复杂岩体最小长方体区域的长宽高X、Y和Z，根据目标自由度的数量m，采用立方体网格的情况下，假设z方向有n个网格，则

由式（4）确定n近似值，取整后可得到x、y、z方向上网格数量，分别为Xnum、Ynum、Znum。故可对该长方体区域进行参数化六面体网分割，获得Xnum×Ynum×Znum个结构化六面体网格，并赋予编号IDs，其网格节点编号s形式为（x、y、z），按照x、y、z坐标从小至大排序。

（3）空属性六面体网格检索。每个六面体网格质心坐标与背景网格中对应的四面体（如图2），通过其在背景网格中的网格编号IDt，检索其所从属的复杂岩体网格编号IDs。判断空间中一点是否位于四面体内，可计算空间点与四面体任意3个顶点组成的4个小四面体体积之和，并比较其与该四面体体积的大小以作判断：当两者相等时，空间点位于四面体内或者四面体边界上；当两者不等时，空间点位于四面体外。将位于所有四面体外的六面体定义为空属性六面体网格，对于上亿级自由度的网格，该过程可采用并行计算加速六面体网格的复杂岩体属性检索过程。

（4）空属性六面体网格剔除。由于背景网格在岩体界面上的非连续性，以及岩体开挖运算的存在，极易导致在岩体界面上出现零星的空属性六面体，甚至在岩体的开挖区出现成团聚集的空属性六面体集合。此时需要人为地将上述空属性六面体网格剔除或赋值，并对全部网格进行重新排序，更新六面体单元的拓扑关系。

（5）亿级自由度网格文件导入RFPA3D-Parallel。编写RFPA3D-Parallel的输出接口，导出节点坐标、边界约束、单元节点编号、节点力4个文本文件，将4个结果文件拷贝到RFPA文件所在目录，新建工程并保存后即可导入。

2 在大孤山露天矿的应用

2.1 工程地质概况及数值建模

大孤山铁矿自1916年开始开采以来逐渐成为亚洲典型的深凹露天生产基地。矿场封闭圈标高+90 m，沿走向长1 700 m，宽1 500 m，现已开采到-330 m标高，开采深度已达420 m，设计开采深度为528 m。

大孤山露天矿地层主要由混合岩岩体、太古代花岗岩岩体、花岗斑岩岩体、绿泥石英片岩岩体、千枚岩岩体、玢岩岩体、小孤山矿体、大孤山矿体和低品位矿石条带矿体等组成，各地层有明显的岩性分界面。

为了构建三维数值模型，进行大孤山矿区三维数值模拟分析，利用现有获取的地质资料对矿区进行了详细三维地质模型构建。采用Rhino软件构建断层、岩性分界面和开采境界的三维NUBERS曲面，建立模型区域三维长方体实体模型（2 000 m×2 000 m×890 m），基于布尔运算分割三维实体模型，最终形成图3所示的露天矿复杂地质体三维模型，共计19种复杂岩体结构，27个露天境界设计台阶，围岩中包含2条断层，多个矿体以及包括断层破碎带在内的9种不同岩性。

2.2 数值模拟

对大孤山三维地质体模型进行万级自由度的四面体网格划分，因其不直接参与有限元计算，无需保证各地质体交界面上网格的严格连续性。因此，全过程可自动化实现，且较少的网格即可准确地表达出复杂地质体的几何形态，如图4所示。经统计，网格数量为119 077个，节点数量为32 515个。

根据1亿自由度的网格目标，采用式（4）可确定x、y和z方向上各有420、420和188个网格，共计33 498 549个节点，此为达到1亿自由度所需的最少节点，建立结构化六面体网格后，采用并行计算加速六面体网格的复杂岩体属性检索过程，对于空属性的六面体集合，根据有限元的计算需要剔除或者赋值开挖属性，空属性网格的空间分布如图5所示。

将剔除空属性的六面体群的结构化网格后，对节点重新编号，并生成网格计算文件，导入RFPA3DParallel软件中，如图6所示，单元数量为33 163 200个，节点33 498 549个，自由度达100 495 647。

大孤山铁矿矿体及其围岩涉及的全部岩石材料的物理力学性质见表1。

2.3 数值模拟结果分析

通过数值求解，得到大孤山铁矿损伤场和弹性模量场如图7所示。首先进行重力平衡计算，随后采用离心加载法［28］逐步增加重力加速度，离心加载系数取值为0.4。计算结果显示，F14断层带与F15断层带之间的楔形体矿石条带出现了滑移现象，如图7所示，数值模拟结果与大孤山铁矿观测到的小范围滑动破坏的位置、规模基本相符。

图7表明大孤山铁矿的西北帮为滑坡灾害发生的高风险区域，并显示了该区域楔形体滑移形成过程中损伤演化和裂纹扩展等过程，据此可将整个楔形体滑移的孕育和发生过程分为多个阶段。

Step1～Step6，模型在自重应力下，整体逐渐发生沉降和压实，损伤场和弹性模量场几乎没有任何变化。

Step7～Step13，因模型各部分岩体强度存在差异，以及各部分岩体之间结构面的存在。模型各部分岩体在自重应力下相互挤压、缓慢变形和错动，这造成了坡体断层面处首先出现应力集中，继而出现纵向剪切损伤和裂纹，但由于其坡脚下部台阶较厚大，抑制了边坡的整体变形。

Step14～Step20，楔形体上部已有裂隙在岩体自重应力作用下持续发育、扩展和贯通，部分裂隙贯穿岩体，坡体表层出现滑塌后，坡体失稳加快。矿石条带沿着F14断层带与F15断层带组成的楔形滑面发生整体错动，局部岩块滚落，并且在滚落过程中迅速解体，滑至坡脚时已碎裂成块。

大孤山铁矿边坡的西北帮的变形及局部滑塌主要是由F14断层带与F15断层带组成的楔形造成的。岩体揭露后，坡体暂时没有出现明显的损伤和裂纹扩展行为；在自重应力的持续作用下，矿石条带与F14断层、F15断层相互挤压、错动，导致应力集中，继而在坡体上部表面出现损伤区和裂纹；随着损伤区和裂纹的继续发育、扩展和贯通，楔形体在自重应力作用下发生了整体滑移。上述现象表明，大孤山露天矿西北帮F14断层与F15断层之间的楔形矿石条带存在较大滑坡风险。

3 结论

针对矿山三维复杂岩体大尺度建模与求解计算难点，本研究在进行大孤山露天矿亿级自由度建模的基础上，利用RFPA3D-Parallel软件平台研究了大孤山露天矿西北帮滑坡灾害形成过程，得到如下主要结论：

（1）实现了露天矿工程尺度复杂岩体的结构化六面体网格精细构建，且基于并行计算平台实现了大孤山露天矿边坡稳定性的亿级自由度的计算分析。

（2）数值模拟结果表明，大孤山露天矿西北帮的F14断层带与F15断层带之间楔形矿石条带出现局部损伤区，该结果与现场实际观测情况完全吻合，即该楔形体有可能引起滑坡灾害，应采取加固、削坡等方式进行治理。

（3）大孤山露天矿边坡楔形体的损伤具有自上而下的发展演化特征，即损伤首先发生在上部台阶，随着裂隙向下逐渐扩展、贯通，将最终导致楔形体沿断层发生错动、滑移，这一认识对现场滑坡监测预警及防控治理具有重要的指导意义。