胡颖鹏, 任凤玉, 丁航行, 付煜
东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2017-03-28
基金项目：国家自然科学基金重点资助项目(51534003)；“十三五”国家重点研究发展计划项目(2016YFC0801600)。
作者简介：胡颖鹏(1988-)，男，河北邯郸人，东北大学博士研究生;
任凤玉(1956-)，男，内蒙古赤峰人，东北大学教授，博士生导师。

摘要：RFPA^3D软件前处理功能薄弱, 无法有效便捷地建立形态复杂的岩体模型, 而采矿专业软件SURPAC具有强大的复杂地质精确建模能力, 且两款软件模型单元均为六面体单元, 以此作为研究基础, 提出利用SURPAC软件快速准确地进行RFPA^3D模型生成的方法.根据SURPAC块体单元和RFPA^3D模型单元之间的数据关系, 在MATLAB中编制以上两种软件的模型转换程序STOR.m; 并在某露天铁矿挂帮矿开采边坡失稳分析中进行了应用.计算过程及结果表明, 复杂岩体模型的转换高速便捷、准确可靠, 同时拓展了RFPA^3D和SURPAC的应用空间.
关键词：SURPACRFPA^3DMATLAB三维模型转换程序
RFPA^3D Modeling and Application for Complex Rock Masses Based on SURPAC
HU Ying-peng, REN Feng-yu, DING Hang-xing, FU Yu
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: DING Hang-xing, E-mail: hangxingding@mail.neu.edu.cn
Abstract: Pre-processing function for RFPA^3D is not so strong that rock mass model with complex geomorphology could not be effectively and conveniently built. While SURPAC, a professional mining software, has a strong ability in modeling complex geological body. Besides, the requirement of unit types for both two softwares is hexahedral. On this basis, a method to build the model in RFPA^3D from SURPAC was proposed. According to the data relationship between the two softwares, a conversion program named STOR.m was developed in MATLAB to achieve the model conversion from SURPAC to RFPA^3D. The program was used in a slope instability analysis model in an open pit iron mine. Results show that the model conversion method can be applied conveniently and extends the application of SURPAC and RFPA^3D.
Key Words: SURPACRFPA^3DMATLABthree-dimensional modelconversion program
随着计算机技术的高速发展, 数值模拟技术成为分析岩土工程问题的重要手段^[1].但在很多岩土工程数值模拟中, 建模等前处理工作复杂耗时, 特别是面对复杂不规则的地质体, 常规的前处理方法往往无能为力, 这影响了数值模拟技术在岩土工程领域的应用.
RFPA^3D是在RFPA^2D的基础上扩展开发的一种基于有限元应力分析和统计损伤理论的材料破裂过程分析软件, 能够模拟岩土工程中岩质材料渐进破裂直至失稳的过程^[2-4].文献[5-7]分别用RFPA^3D软件对云阳县云阳大道边坡和锦屏水电站左岸岩质边坡稳定性进行模拟, 成功得到两种边坡失稳破坏过程.但RFPA^3D软件前处理功能薄弱, 自身无法建立复杂三维模型, 同时其单元必须为六面体单元, 这又限制了很多前处理软件在RFPA^3D模型建立中的应用.
SURPAC是一款先进的矿业软件, 具有强大的复杂地质体精确建模能力, 可以准确表现出不同地质结构体及其在空间中的分布, 还可建立六面体单元的块体模型^[8].将SURPAC和RFPA^3D联合应用, 可以各取所长, 实现对复杂地质体的RFPA^3D建模.类似的SURPAC与岩土分析软件的联合应用已经得到实现, 林杭等^[9]采用Fortran语言编制SURPAC-FLAC^3D接口程序, 实现模型转换; Li等^[10]通过MATLAB软件编制接口程序, 将SURPAC模型成功导入FLAC^3D中.
本文在SURPAC中建立块体模型, 输出块体模型的单元数据, 利用MATLAB软件的m文件编制SURPAC-RFPA^3D的模型转换接口STOR.m, 将SURPAC块体模型转换为RFPA^3D计算模型.利用该转换方式, 进行某露天铁矿东边坡挂帮矿开采的模拟, 模拟结果与现场观测结果基本吻合, 证明了SURPAC与RFPA^3D耦合计算的有效性.
1 SURPAC模型构建SURPAC对于地质体模型有两种表现形式:一种是实体模型, 另一种是块体单元.实体单元可以形象地展示出地质结构体的形态和分布, 例如岩层、空区和断层等, 但是这种实体模型本质仅为一层蒙皮, 其内部没有划分单元, 即内部无有效的单元数据.基于实体模型可以建立一种更准确的三维块体模型, 它包括大量离散的规则六面体单元, 这些块体模型不仅可以准确地建立地质体特性与空间位置的关系, 而且其本质上是一个数据库, 包含顶点坐标、岩性种类等数据, 恰好可以对应RFPA^3D建模所要求的六面体单元.
SURPAC构建块体模型的步骤如图 1所示, 利用地质平面图、剖面图和地表等高线图等二维图纸, 建立三维实体模型, 根据实体模型和不同岩性的约束情况建立块体模型, 并为块体赋属性.
图 1(Fig. 1)

图 1 块体模型构建步骤Fig.1 Steps of constructing block model

将SURPAC块体模型的单元信息导出为EXCEL格式文件, 该文件包含的信息为块体模型每一单元的质心三维坐标、单元三边长度和该单元的属性.
2 SURPAC与RFPA^3D模型转换SURPAC中, 块体单元模型的形状为长方体, 这满足RFPA^3D所要求的单元必须为六面体的原则.RFPA^3D通过4个文本文件接受外部导入的模型, 分别为coor0.txt, elem0.txt, id0.txt, disp0.txt.
要完成模型的转换, 需要对所提取的SURPAC块体模型数据进行处理, 变为RFPA^3D可接受的4个文本文件.
coor0.txt文件为RFPA^3D模型所有节点的三维空间坐标信息文件.RFPA^3D的单元为六面体, 其8个顶点(P₁~P₈)的编号顺序见图 2.SURPAC块体单元输出的信息为单元质心的坐标(X₀, Y₀, Z₀)、三维尺寸L₁, L₂, L₃和块体单元属性的编号^[9].图 2展现了2种模型数据的关系, 根据几何关系得到两者坐标的转换方程:
图 2(Fig. 2)

图 2 RFPA3D与SURPAC单元数据关系Fig.2 Relationship of element data between RFPA3D and SURPAC (a)—顶点编号图；(b)—质心坐标与单元尺寸.

式中，X_P1, Y_P1, Z_P1分别表示点P₁的X, Y, Z坐标.
同理可得P₂~P₈各点的三维坐标.
将所得到的单元顶点数据进行筛选, 删掉重复的顶点数据, 将筛选后的顶点数据按照G_n, X_n, Y_n, Z_n的格式写入文本coor.txt, n为筛选后的顶点数目, G_n为顶点编号, X_n, Y_n, Z_n为相应顶点的坐标.
elem0.txt为记录模型所有单元数据的文件, 包括RFPA^3D模型中每个单元的顶点编号和该单元的属性编号.在筛选后的顶点中, 查询每个单元包含的顶点编号, 按照Z_m, P₁, P₂, P₃, P4, P5, P₆, P₇, P₈, V的格式将顶点编号写入elem0.txt文件, m为单元数目, Z_m为单元编号, P₁~P₈为单元顶点编号(如图 2所示), V为单元属性编号.
id0.txt为模型边界条件的信息文件, 包含每个顶点X, Y, Z三个方向的约束情况.RFPA^3D规定, 节点位移约束为-1, 应力约束为1.搜索边界的顶点、区分边界处顶点和模型内部顶点, 按照G_n, X_n, Y_n, Z_n的格式, 将顶点信息写入id0.txt文件, X_n, Y_n, Z_n为相应顶点的三个坐标方向的约束条件, 为1或-1.
disp0.txt为模型节点的初始位移信息文件, 包含每个节点X, Y, Z三个方向的位移情况, 一般默认为0.将筛选后的顶点数据按照G_n, X_n, Y_n, Z_n(X_n, Y_n, Z_n一般默认为0)的格式写入disp0.txt.
根据上述的SURPAC到RFPA^3D模型转换的方法, 在MATLAB中利用m语言编制模型转换程序STOR.m.该程序将SURPAC导出的块体模型单元数据转换为RFPA^3D可接受的4个数据文件, 再导入RFPA^3D中，分别给不同组岩体赋予力学参数, 即可进行计算.具体转换流程见图 3.
图 3(Fig. 3)

图 3 SURPAC与RFPA3D转换流程图Fig.3 Flow chart of the conversion from SURPAC to RFPA3D

3 工程应用某大型露天铁矿露天生产进入末期, 为延续产能, 利用无底柱分段崩落法开采东部边坡下部的挂帮矿, 由此引起东部边坡失稳破坏, 边坡及地表产生大规模岩移.本文应用上述方法, 在SURPAC中对该露天矿东部边坡进行建模, 将模型导入RFPA^3D中模拟计算边坡破坏失稳的情况.
选取该矿东边坡作为模拟区域, 所选区域南北方向长700 m, 东西方向长600 m, 由地表至-200 m水平; 开采区域位于边坡下部矿体中, 垂直范围为-51~-123 m, 水平范围为南北走向100~140 m, 东西走向约130 m.利用该矿的平面图和剖面图, 在SURPAC中建立东边坡的实体模型, 结合约束条件将实体模型转化为包含3组属性的块体模型.如图 4所示, 模型北部为千枚岩, 南部为混合岩, 中部为矿体; 边坡下挂帮矿开采区域见图 5.所建模型单元为尺寸统一的立方体单元, 单元南北走向边长7 m, 东西走向边长6 m, 垂直走向边长5 m, 模型单元共计542 158个.
图 4(Fig. 4)

图 4 某露天铁矿东边坡SURPAC块体模型Fig.4 Block model of east slope in a some open pit mine in SURPAC

图 5(Fig. 5)

图 5 东边坡下开采区域Fig.5 Excavating zone under the east slope

将块体模型的质心数据文件导出, 采用本文中提出的转换方法, 在MATLAB中通过SURPAC- RFPA^3D接口程序STOR.m, 导入块体质心文件, 导出coor0.txt, elem0.txt, id0.txt和disp0.txt 4个文件, 将这4个文件导入RFPA^3D中, 赋予各岩性单元相应的岩石力学参数, 如表 1所示.生成的RFPA^3D计算模型如图 6所示.
表 1(Table 1)

表 1 岩体力学参数岩体力学参数Table 1 Mechanical parameters of rock masses 岩性 弹性模量 泊松比 密度 内聚力 摩擦角 GPa g·cm-3 MPa (°) 千枚岩 2.5 0.30 2.77 0.5 38 矿体 4.1 0.26 2.65 1.0 41 混合岩 3.5 0.26 2.60 0.8 38

表 1 岩体力学参数岩体力学参数 Table 1 Mechanical parameters of rock masses

图 6(Fig. 6)

图 6 某露天铁矿东边坡RFPA3D计算模型Fig.6 Numerical model of east slope in RFPA3D

在RFPA^3D中, 对东边坡下的挂帮矿进行开挖(开挖区域见图 5).由图 7所示的计算结果分析可知, 空区距离边坡最近的一侧(面向露天坑一侧)会首先发生破坏, 破坏沿空区边缘向南北两边发展, 产生一道连续的南北走向的破坏带后, 空区上方边坡岩体失去下部的拱脚支撑后, 发生整体滑移, 向采空区塌陷, 地表形成崩落区.
图 7(Fig. 7)

图 7 东边坡破坏过程岩体损伤变化图Fig.7 Damage evolution of rock mass during failure processes in east slope (a)—RFPA3D计算步骤第45步；(b)—RFPA3D计算步骤第63步.

为验证模型转换方法的可靠性, 建立与RFPA^3D模型尺寸和单元划分形式相同的Flac^3D模型, 采用表 1的岩石力学参数, 进行挂帮矿开挖计算, 得到如图 8所示的计算结果, 选取位移≥2 m的区域为崩落区域^[11], 得到挂帮矿开采引起的崩落范围.
图 8(Fig. 8)

图 8 东边坡Flac3D计算位移云图Fig.8 Displacement contour of east slope obtained from Flac3D

将RFPA^3D和Flac^3D模拟计算得到的边坡崩落范围与现场实际破坏情况进行对比分析, 如图 9所示, 模拟结果与现场实际观察的破坏情况基本吻合, 很好地解释了破坏发展的过程, 证明了SURPAC-RFPA^3D接口的有效性.
图 9(Fig. 9)

图 9 东边坡实地崩落范围与计算范围Fig.9 Comparison of scopes obtained from field-based monitor and simulation

4 结论1) SURPAC块体单元和RFPA^3D计算单元的形状均为六面体, 以此为基础, 将SURPAC地质模型转化为RFPA^3D计算模型, 利用SURPAC强大的建模能力来弥补RFPA^3D在复杂模型建模中的不足.
2) 根据所建立的SURPAC和RFPA^3D模型单元数据关系, 在MATLAB中使用m语言编制了转换程序STOR.m.
3) 使用SURPAC和RFPA^3D模型转换程序STOR.m进行某大型露天矿东边坡挂帮矿开采的模拟计算, 结果表明:所编制的转换程序能够快速准确地进行模型转换, 并与Flac^3D计算结果吻合, 转换程序有效、可靠.
参考文献

[1]	侯恩科, 吴立新, 李建民, 等. 三维地学模拟与数值模拟的耦合方法研究[J].煤炭学报, 2002, 27(4): 388–392. ( Hou En-ke, Wu Li-xin, Li Jian-min, et al. Study on the coupling of 3D geoscience modeling with numerical simulation[J].Journal of China Coal Society, 2002, 27(4): 388–392.)
[2]	Tang C A, Tang A B. Applications of rock failure process analysis (RFPA) method[J].Journal of Rock Mechanics and Geotechnical Engineering, 2011, 3(4): 352–372.DOI:10.3724/SP.J.1235.2011.00352
[3]	Liang Z Z, Tang C A, Li H X, et al. Numerical simulation of 3D failure process in heterogeneous rocks[J].International Journal of Rock Mechanics and Mining Sciences, 2004, 41(1): 323–328.
[4]	Dai F, Wei M D, Xu N W, et al. Numerical assessment of the progressive rock fracture mechanism of cracked chevron notched Brazilian disc specimens[J].Rock Mechanics and Rock Engineering, 2015, 48(2): 463–479.DOI:10.1007/s00603-014-0587-8
[5]	Li L C, Tang C A, Zhu W C, et al. Numerical analysis of slope stability based on the gravity increase method[J].Computers and Geotechniques, 2009, 36(7): 1246–1258.DOI:10.1016/j.compgeo.2009.06.004
[6]	徐奴文, 唐春安, 周钟, 等. 基于三维数值模拟和微震监测的水工岩质边坡稳定性分析[J].岩石力学与工程学报, 2013, 32(7): 1371–1381. ( Xu Nu-wen, Tang Chun-an, Zhou Zhong, et al. Stability analysis of hydraulic rock slope based on three-dimensional numerical simulation and microseismic monitoring[J].Chinese Journal of Rock Mechanics and Engineering, 2013, 32(7): 1371–1381.)
[7]	Xu N W, Tang C A, Li L C, et al. Microseismic monitoring and stability analysis of the left bank slope in Jinping first stage hydropower station in southwestern China[J].International Journal of Rock Mechanics and Mining Sciences, 2011, 48(6): 950–963.DOI:10.1016/j.ijrmms.2011.06.009
[8]	Luo Z Q, Xie C Y, Zhou J M, et al. Numerical analysis of stability for mined-out area in multi-field coupling[J].Journal of Central South University, 2015, 22(2): 669–675.DOI:10.1007/s11771-015-2569-8
[9]	林杭, 曹平, 李江腾, 等. 基于SURPAC的FLAC3D三维模型自动构建[J].中国矿业大学学报, 2008, 37(3): 339–342. ( Lin Hang, Cao Ping, Li Jiang-teng, et al. Automatic generation of FLAC3D model based on SURPAC[J].Journal of China University of Mining & Technology, 2008, 37(3): 339–342.)
[10]	Li X B, Li D Y, Liu Z X, et al. Determination of the minimum thickness of crown pillar for safe exploitation of a subsea gold mine based on numerical modeling[J].International Journal of Rock Mechanics & Mining Sciences, 2013, 57(1): 42–56.
[11]	Sainsbury D P, Sainsbury B L, Lorig L J. Investigation of caving induced subsidence at the abandoned Grace Mine[J].Mining Technology, 2013, 119(3): 151–161.