金属矿山三维应力计算建模技术

梅婷婷，张电吉，周春梅，杨 曼,王家国

(1.武汉工程大学环境与城市建设学院,湖北 武汉 430074；2.湖北宜化集团矿业有限责任公司，湖北 宜昌 443000)

0 引 言

矿山三维地质体建模技术研究，是“数字地球”、“数字矿山”的核心组成部分，是现代矿山信息化研究的热点和难点.近年来，国内外学者吸收科学可视化、计算机图形图像学中的可视化技术在三维地质建模理论与方法方面进行了许多有益的探索[1].侯恩科[2]、吴立新、王占刚[3]、齐安文[4]、侯遵译[5]、张菊明[6]、柴贺军[7]、毛善军[8]等学者对建模理论和方法的研究和探讨，为三维可视化模型的研制和开发奠定了基础.针对金属矿山开挖的建模过程都极为复杂，尤其是建立三维实体模型，如果采取传统的建模方法，不仅需要花费大量的时间，而且还容易出现错误[9].FLAC3D在解决岩土问题上具有很大的优越性，是求解岩土工程问题的最理想工具之一,但是在模型建立以及单元网格划分等前处理方面却一直是软肋.相比之下，ANSYS通用有限元软件对于复杂模型的实体建模，其前处理功能优点尤为突出.笔者以“点→线→面→体”自底而上的实体建模技术路线为指导，以某金矿为应用实例，试图利用ANSYS 软件对复杂工程地质体建立相应数值模型，再通过数据转化实现 FLAC3D有限元模型的自动生成.针对建模过程中所遇到的诸如数据获取和模型简化等问题提出了具体的解决方法，以供其它类似矿山建立应力场计算模型时参考.

1 实体建模理论与方法

目前，三维地质建模使用的方法基本是基于空间曲面插值拟合方法来实现，主要有趋势面分析法、克里格法、形函数法和人工神经网络模拟法等[10].本模型主要在形函数和趋势面分析法的理论基础上，借助ANSYS自底向上的建模技术，建立三维地质体模型.ANSYS各图元之间的关系如下：

(低→高)

1.1 由关键点生成样条曲线

以研究目的为前提，对模型进行适当的简化，结合具体的资料，确定模型数据采样点(关键点)的数量及三维坐标.在ANSYS实体建模中，所有线均用非均匀有理B样条表示[11]，利用公式(1)对数据采样点进行插值计算，形成由关键点生成的样条曲线.

(1)

Pi为数据采样点；Ni,k(u)是k次规范B样条基函数；wi为权因子(权重系数)；若基函数的节点是非均匀的，则称为非均匀有理B样条.

1.2 由样条曲线构造样条曲面

对于复杂金属矿山地表模型的建立，主要运用“蒙皮法”.蒙皮法传统上称为放样，其任务就是拟合一张曲面(即“皮”)，使其通过一组有序的截面线.本文引用文献[12-13]对样条曲线进行蒙皮操作，借助公式(2)生成一张双三次B样条曲面.具有(m+1)(Nmax+2)个控制点的双三次B样曲面的方程为：

(2)

式(2)中：di,j(i=0,1,…,Nmax+1)是曲面的控制点；Ni,3(u)是定义在同一节点矢量上的B样条函数，Nj,3(v)是定义在另一个参数方向的节点矢量V上的B样条函数.

针对矿岩地质界面拟合，主要采用趋势面分析方法.本模型采用白润才[14]提出的趋势面分析法的改进方法，由关键点插值生成曲面.如果平面上分布的某参数用Z表示，那么趋势面分析可简单看成寻求一个方程，按二次多项式拟合空间曲面为：

Z=f(X,Y)

(3)

P(X,Y)=C00+C10X+C20Y+C30X2+

C40XY+C50Y2

(4)

∂Q/∂Crs(r,s=0,1,2)

(5)

可得到6个未知数的6个方程组，求得式(5)的系数Crs.

式(3)～(5)中，(Xi,Yi)为岩体所在范围某平面的数据点坐标；Zi为点(Xi,Yi)上的估测点；Q为所有数据点的拟合值P(Xi,Yi)和估测值Zi的误差的平方和；n为估点数；w(di)=w[(Xi-a)2+(Yi-b)2]为权值，是距离的函数.

1.3 三维实体的生成

由已生成的关键点建立样条曲线，再根据样条曲线构造相应的曲面，最后由一系列连续封闭的曲面定义体.强大的布尔运算工具可以实现实体之间加、减、分类、搭接和分割等复杂运算，大大提高了建立复杂地质体三维模型的效率.对于实体模型的网格划分，ANSYS提供了功能强大的控制工具，还可对实体模型图直接划分网格[15].

2 采用ANSYS实体建模的步骤

对于整个有限元模型的建立，其主要的建模步骤如下：

a.建模方案的确定.在开始进入ANSYS之前，首先收集资料并采集多元数据，选择合适的单元类型，考虑适当的网格密度，并确定单元的实常数、材料属性和单元坐标系.

b.模型简化.在建模前，首先熟悉工程地质特征和明确自己的建模对象，通过经验定性分析影响矿体开挖的主要因素和次要因素，确定是否具备简化条件，即所谓的模型简化.

c.几何模型的建立.进入前处理(PREP 7)模块，利用几何数据，遵循从“点→线→面→体”自外而内、自底而上的建模技术，生成符合几何形状的独立实体模型，并借助布尔操作将各个独立的实体模型适当的连接在一起.

d.网格剖分.考虑研究对象的影响程度，依次设置不同的网格划分精度，对实体模型划分网格，生成节点和单元.

e.接口与转换.运用ANSYS to FLAC程序,将d中的节点和单元导入到FLAC有限差分软件，进行计算.

3 建模过程中存在的难点及解决方法

ANSYS的建模是一个复杂的过程，计算模型的建立是应力计算中最为困难的一步，模型最终将直接影响到计算分析结果的真实性和可靠性[16].

3．1 数据获取

建模的关键取决于数据采样点坐标的准确性，所以首先要解决的问题是如何获取精确数据.由于资料的详尽程度不同，图形上可能有部分数据缺失，此时可依据周围地形地貌的走势，借助“线性内插法”对数据采样点的高程值进行估算.在结合水平平面图以及剖面图读取数据点坐标时，要注意角度转化问题.在地表建模过程中，参照现场实际情况，对描述露天采场特征的数据要进行适当的修正，力求与工程实际开挖情况相符.

3．2 模型简化

复杂矿山的三维地质体建模最大的难点就在模型简化.目前对断层的三维表示与模拟研究还处于研究和探索阶段，要模拟形态各异的金矿，也存在技术上的困难.鉴于模型主要用于模拟矿体开挖，从系统工程的角度出发，坚持“重点研究局部，局部服从整体”的原则，可对矿体模型和断裂带模型进行适当的简化.

a.三维实体的建模技术.该模型主要用于应力场计算，借助基于面的构模方法的思路，将三维实体假想成不规则封闭曲面构成的六面体模型.虽然矿体的形状、产状、品味和大小各异，在建模过程中可其简化为六面体模型.根据矿体大致呈脉状以及矿体品位可圈定拟建矿体模型的外轮廓线，构造封闭曲面，从而建立三维实体模型.值得注意的是：第一，将圈定范围内不同品位的矿体与存在的夹层，均视为品味相同的矿体，尽量使矿体形状由不规则变为规则.第二，为了便于建模和对模型进行修改，一般选取实体轮廓线的特征点为数据采样点(关键点)，宜尽量保持采样点数量的统一、采样点编号顺序的统一.

b.三维曲面的建模技术.断层一般由断层面来表示其产状、断层带的宽度与延伸尺寸由控制点处得几何特征来定义.由地质构造判断断裂破碎带的稳定性，确定其对数值计算的影响程度.如果该断裂带较为稳定，属于控矿断裂带且矿体伴随断裂带产生，可假定该断裂带靠近矿体的一侧与矿体大致重合，避免夹层对网格剖分产生不利影响.考虑到断裂带的控制长度、宽度变化以及影响程度，取适当的数值作为断裂带的平均宽度.以假定重合的曲面为基准面，利用ANSYS的VDRAG命令将基准面向外平移，平移的距离即为平均宽度，最终生成断裂带模型.

除此之外，ANSYS软件的建模方法较为灵活，可利用菜单创建有限元模型或利用命令流的方式建模.在进行布尔运算之前，要判断模型是否满足运算的要求，特别是在边界处.而且尽量减少布尔运算的次数和充分利用模型的对称性，使模型的几何形式尽量简单.

4 应用实例

某金矿产于伊犁-中天山板块北部的古生代中天山北缘活动大陆边缘的吐拉素断陷盆地中[17],区域地层及构造线方向均呈北西西向展布，是古生代的低硫型浅成低温热液金矿.该金矿中心点地理坐标为：东经81°36′30″，北纬44°13′45″，矿区面积0.348 km2.如图1所示，北露天采矿场为山坡地形，整体呈一个长轴近于南北的椭圆形凹陷露天采坑.矿体在地表呈-向西突出的带状分布，呈脉状产出，倾向东-北东，倾角60°～80°.在露天开采范围内，与采矿工程相关的地质构造有2组，即F2断裂破碎带和不整合接触面.在矿区的北部F2产状倾向为100°，向南部倾向变为50°，倾角为55°～85°，局部倾角为90°，大体与矿体延伸方向一致.结合上述工程地质特性，笔者主要以地表建模为例进行详细介绍.

图1 北露天采矿场构造地质图

4.1 地表建模

为了建模方便，地表模型采用相对坐标系XOY，取某金矿地形地质及井上下对照图(1∶2000)，将图上的大地坐标(4899000,27548200)定为XOY平面的起点，起点坐标为(200.00,0.00)，坐标精度保留两位小数，以米为单位，地表模型范围为1 200 m×1 600 m.地表建模的思路：在地形图上构造一系列大小相等、规则的M×N栅格(GRID)，然后根据已知等高线和控制点的高程值，通过线性内插的方法求得每个栅格点的高程值，从而得到各栅格点的三维坐标，按照点成线、线成面的步骤，建立地表光滑曲面.

其中，几何模型的建模过程如下：

a.在地形图上建立20.00 m×20.00 m的规则栅格，将每个栅格四周的拐角点作为数据采集点，本图采集点共4 941个.

b.将数据采集点相邻的等高线值作为已知值，通过数据采集点做两条等高线的公垂线，测量公垂线以及采样点与其中一条等高线的最短垂直距离，利用线性内插的方法，估算各采样点的高程值.最后结合所在地形图的XOY平面，确定各点的X坐标和Y坐标，并按照ANSYS软件指定格式编辑命令流.

c.依据文中提到的非均匀有理B样条插值方法，由关键点生成样条曲线，如图2所示；借助“蒙皮法”将样条曲线生成蒙皮似的光滑曲面，如图3所示；由光滑曲面向下拖拉1 200 m生成体，如图4所示；在800 m高程处建立工作面，利用布尔运算切割各拖拉体后，再将其合并，如图5所示.

d.由于北露天采矿场的地形较为复杂，建模方法是依据北露天采矿场各勘探线剖面图，读取各台阶拐点处的X坐标和Z坐标，再将该点投影至该金矿地形地质及井上下对照图(1∶2000)中各勘探线上的相应位置，读取该点Y坐标.同样利用多次线性内插法，由各台阶拐点的坐标估算北露天采矿场区域内各删格四周的拐角点三维坐标，做法同b所示.

图2 点成线

图3 线成面

图4 面成体

图5 修改后的体

e.在建模过程中考虑不整合界线，并对两组不同的岩性分别进行物理力学性质参数赋值.不整合界限模型主要的建模思路同地表模型一致，即以地层不整合界线的特征点作为关键点，按照点成线和线成面的原理生成边界不整合界面，最后利用VSBA命令，对地层不整合界面与整个地表模型布尔运算，从而形成含边界不整合界面的地表模型，如图6所示.

图6 简化后的边界不整合界面模型

4.2 计算模型的生成

根据已有的地质勘探、矿体资源储量、采矿工艺和流程等资料，依照相同的原理，在地表模型建成的基础上，分别建立断裂构造、矿体、塌落体、上覆岩梁(隔离层)、矿房和矿柱等三维地质体模型.通过定义实体模型的各种属性以及借助ANSYS软件强大的网格划分功能，将该实体模型转化为有限元计算模型.运用ANSYS to FLAC3D程序,将有限元计算模型中的节点和单元导入到FLAC3D有限差分软件，进行矿体开挖的数值模拟计算，如图7所示.

图7 有限元计算模型

5 结 语

以上研究的基于ANSYS的金属矿山三维地质建模方法，不仅克服了建立复杂计算模型的困难，而且操作简单、易行，所要的基本数据仅为几何图形数据点，大大减少了建模所需的时间、精力，具体如下:

a.采用“点→线→面→体”自外而内、自底而上的建模技术，依次建立地表、断裂构造、矿体、塌落体、上覆岩梁(隔离层)、矿房和矿柱等三维地质体模型，再利用布尔运算组装各个独立的实体模型.不仅保证各地质体模型之间的独立性和有限元模型的整体性，而且便于模型的修改和更新，提高了建模效率.

b.ANSYS软件的前处理功能适用于建立庞大复杂的模型,而FLAC3D有限差分软件计算模拟速度较快,充分结合两种软件的优点，解决了复杂工程地质体模型的数值计算难题.该方法已运用于某金矿的实际生产中，实践证明其对矿山企业的安全生产及采矿方法的优化具有指导意义.

c.遵循“重点研究局部，局部服从整体”的简化原则，以建模的目的为出发点，以实际的工程地质情况为前提，对模型进行简化.但是在矿山三维地质建模的过程中，至今还没有形成一套完善的简化规律，该问题仍待进一步探讨和研究.

参考文献：

[1] Simon W Houlding. 3D Geosciences modeling：computer technique for geological characterization[M]．Hong Kong: Springer.1993.

[2] 侯恩科，吴立新．面向地质建模的三维体元拓扑数据模型研究[J]．武汉大学学报：信息科学版，2002,10(5):467-472.

[3] 王占刚，曹代勇．基于改进三棱柱模型的复杂地质体3D建模方法[J]．中国煤田地质，2004,16(1):3-6.

[4] 齐安文，吴立新．基于类三棱柱的三维地质模拟与拓扑研究[J]．矿山测量，2003(3):65-67.

[5] 侯遵译，许振邦．建立在三维空间上的地学计算机图形图像系统[J]．中国数学地质进展，1995，25(7):168-179.

[6] 张菊明．空间地质模型的设计和显示[J]．中国数学地质进展，1995,25(4):189-195.

[7] 柴贺军，黄地龙，黄润秋，等．岩体结构三维可视化及其工程应用研究[J]．岩土工程学报，2001,23(2):217-220.

[8] 毛善军，许友志，张海荣，等．空间地质模型及其可视化系统[J]．中国数学地质进展，2003，25(7):168-179.

[9] 卢宏建，罗涛，甘德清．基于ANSYS的矿山三维动态建模系统研究与开发[J]．金属矿山，2009(1):106-107.

[10] 潘结南，孟召平，甘莉．矿山三维地质建模与可视化研究[J]．煤田地质与勘探，2005,33(1):16-18.

[11] 潘玉松．浅谈使用ANSYS时的建模方法[J]．成都电子机械高等专科学校学报，2000(2):28-30.

[12] Woodward C D． Skinning techniques for interactive B-spline surface interpolation[J].Computer-Aided Design,1988,20(8).

[13] 张秀萍．基于点云数据复杂曲面的三维重构[D]．新疆大学，2004.

[14] 白润才，赵树贤，刘淑让．关于趋势面分析法在矿岩界面拟合应用中的改进[J]．阜新矿业学院学报，1995,14(2):101-103.

[15] 廖秋林，曾钱帮．基于ANSYS平台复杂地质体FLAC3D模型的自动生成[J]．岩土力学与工程学报，2005,24(6):1010-1013.

[16] 李元松，李新平，姜天华．板桁结合梁斜拉桥3D-ANSYS建模技术研究[J]．武汉工程大学学报，2008,30(4):65-68.

[17] 翟伟，孙晓明，苏丽薇，等．新疆阿希金矿：古生代的低硫型浅成低温热液金矿床[J]．地学前缘(中国地质大学)，2010,17(2):266-285.