多软件联合的煤矿三维数值建模方法及其应用*

朱权洁，刘金海，王胜开，高林生，梁　娟

(1. 华北科技学院 安全工程学院，北京 东燕郊 101601；2. 驰宏锌锗股份有限公司，云南 曲靖 655011；3. 防灾科技学院，北京 东燕郊 101601)

0　引言

矿井地应力场测量与反演是确定工程岩体力学属性、进行围岩稳定性分析以及地下工程开挖设计科学化的必要前提[1-2]，也是有效防治冲击地压等灾害、保证矿井安全有序生产的关键。目前，地应力测试都是在硬岩层中进行的，煤层应力值除可根据已知点对应岩层地应力实测结果进行估算外，目前尚无其他较好办法，但这种方法也具有局限性和不稳定性。要全盘掌握矿区内不同深度、不同位置、同一煤层不同深度上的应力分布，需进行区域地应力场的反演分析，这其中地质模型的数值计算是关键步骤。

数值计算是当前解决科研难题、实际工程问题的有效手段[3-5]。在岩土工程领域，有限元法、离散元法以及有限差分法等是数值分析计算的重要方法之一。Itasca公司研发的FLAC3D数值模拟软件是其中的优秀代表之一，该软件对于模拟不产生离散、连续变形模型具有良好效果，因此，在岩土工程、采矿工程等行业和领域应用广泛。

在具有强大数值分析能力的同时，FLAC3D的建模尤其是复杂模型的构建能力是其最大短板。对于规则模型而言，利用fish语言编写命令流可实现快速建模，但不规则、复杂的地质模型依靠命令流实现是非常困难的，这也是FLAC3D面临的一大挑战[6]。现有研究成果中多选择“先外部构建几何模型、后导入FLAC3D”的方式开展[7-9]。如林杭等编制了Surpac到FLAC3D的接口程序，直接实现Surpac三维模型到FLAC3D的转换；纪洪广等[[10]提出一种利用VC++开发的复杂矿山地质形体网格精细划分方法，解决了FLAC3D构建复杂地质体费时、繁琐等缺陷；赵文等[11]提出利用MATLAB和Surfer划分网格，并生成FLAC3D命令流语句，提高了FLAC3D建模速度。这些方法中前者在网格划分方面存在缺陷，后者对使用者的编程能力要求较高。

此外，袁海平等[12]提出联合AutoCAD，ANSYS和FLAC3D联合构建三维地质模型，该方法在CAD中制作面元素，仍然受限于构建大型复杂三维模型；李翔等[13]基于钻孔柱状图资料，提出利用SURPAC，CAD和ARCGIS建立复杂三维地址模型。该方法一定程度上解决上了上述问题，但在利用ANSYS的APDL语言生成模型时涉及大量命令流操作，影响模型构建的速度，如果遇到命令错误会导致模型无法完整生成。由此可见，如何有效、快速构建复杂三维地质模型进而开展数值计算是地应力反演的前提和重要研究课题[14-15]。

鉴于上述原因，本文提出了一种复杂三维地质模型的构建方法，该方法联合3DMINE，GOCAD，SURFER，RHIRO，ANSYS以及FLAC3D以及MATLAB共7种软件共同建模，减少了编程工作量、兼顾了ANSYS强大的前处理功能，弥补了大量命令流编写、无法合理划分网格等缺陷，实现了复杂三维地质体的快速构建。

1　建模思路

以山东郓城煤矿为例，对整个建模思路进行介绍。为了客观反演地应力的分布规律，综合考虑地形、地质构造以及地层介质等因素，结合前期开展的地质勘探工作，最后确立了目标研究区域的范围：以坐标点A(7 344.8，3 958.1)，B(9 324.8，3 958.1)，C(9 324.8，9 998.1)和D(7 344.8，9 998.1)圈定的长方形区域为计算区域，如图1中所示。

图1　几何模型构建区域Fig.1　Construction area of geometric model

为了建立该的三维地质模型，使该模型更贴进实际的地形地貌，减少因数值模拟不精确带来的地应力反演误差，开展了矿区复杂三维地质模型精细化构建研究_建模原理基于“点-线-面-体”逐步构建的思路，在构建过程中使用了GOCAD，ANSYS，FLAC3D等多种软件，具体过程及思路如图2所示。

图2　三维地质模型构建过程及思路Fig.2　Process of 3D numerical modeling

2　点阵数据和面的形成

常规数值模拟一般采用规则建模方法构建模型，这类建模方式由于软件自身和建模方式不灵活从而无法完整表现地质体的起伏、断层构造等信息。为了精细地描述复杂的地层结构，提出利用钻孔信息生成岩层界面等高线，并“以点构线、以线构面、以面构体”的三级建模法构建模型，以形成更真实的复杂的地质体模型，提高数值模拟的精确度和可靠性。

2.1　等高线的生成

勘探钻孔是构建岩层界面等值线的数据来源。钻孔含有大量信息，通过钻孔柱状图可以了解矿体、岩层界线以及断层构造等隐藏于地下的信息。根据这些信息，可以生成相应的岩层界面等高线。利用CAD建立三维等高线操作较为麻烦，因此，选择借助SURPAC或3DMINE等数字矿山三维软件直接生成等高线，其操作步骤如下。

1)准备钻孔原始数据。从柱状图中读取相应的岩性信息表、钻孔空间坐标，并导入到三维软件之中。

2)三维钻孔形态的生成。根据数据库中存储的钻孔信息，生成三维形态的钻孔立体图，如图3(a)所示。图中以开孔位置为起点、终孔位置为终点连接成线，并按照钻孔岩性表中的岩性特点分段标记。

3)三维等高线的生成。根据这些不同地层的不同标记，提取相应地层的点坐标，这些点坐标经过分步取样、插值即可生成等高线图，如图3(b)所示的煤层顶板等高线。

图3　利用钻孔信息生成等值线Fig.3　Contour map generation using borehole data

钻孔数量的多少影响数值建模的精细程度，因此，可根据工程实际要求设计不同的钻孔数量需求。如存在大量外部数据(如等高线)，可直接进行加载，导入GOCAD进行离散处理。

2.2　离散点阵的生成

由于前期钻孔信息有限，所生成的等高线上仅存有少量离散的拟合点坐标数据，这些数据如果直接导入ANSYS或FLAC3D进行建模，一方面将造成模型偏差较大(严重失真)，另一方面利用少量点构建体操作性不强。因此，需要对等高线进行离散化和细分处理，生成相应的等高线拟合点数据。

经过3DMINE处理后生成dxf文件，并导入到GOCAD中。GOCAD读取dxf数据，提取出等高线中拟合点的数据。如图4所示，为煤层顶板等高线的网格点阵图，其中包括等高线和阵列点的空间坐标信息。

图4　利用GOCAD生成网格点阵坐标(煤层顶板)Fig.4　Grid lattice coordinates using GOCAD

2.3　数据网格化处理

前期得到的GOCAD导出的等高线坐标点数据有限且较为离散，直接导入ANSYS构建实体模型较为复杂，因此，借助SURFER软件对离散数据进行插值处理，使数据点均匀化，即网格化操作。

利用SURFER软件“数据”选项读取上一步获得的离散等高线拟合点数据，选择“克里金插值算法”对离散数据进行插值和均匀化处理，即得到等高线的点阵数据集，生成线框图如图5所示。

图5　点阵数据构建的线框图Fig.5　Wireframe generation using grid lattice coordinates

需要注意的是，由于地质模型在导入ANSYS中时会考虑模型大小问题，因此，需要对钻孔的绝对坐标进行偏移，以保证模型的顺利生成。坐标系偏移为：X偏移-20479 680，Y偏移-3911 100(这是后文进行反演时的重要参数，以此对比CAD图纸上任一点与几何模型坐标的关联性)。

2.4　不规则面的生成

SURFER提供的点坐标数据无法直接导入ANSYS模型，借助命令流过程复杂。因此，为减少编程带来的工作量，将SURFER生成的dat文件导入RHINO三维软件，并利用RHINO提供的srfptgrid命令将数据点连线成面，如图6所示，最后将该面文件转换为ANSYS可读的iegs格式文件(该格式可供ANSYS直接读取)。如图6所示，利用RHINO软件生成的三维面模型，分别展示了该面模型的俯视图、前视图、右视图和透视图。

图6　模型的生成与格式转换Fig.6　Face model generation and data format conversion

3　模型的建立与数值计算

FLAC3D是采矿、岩土工程等领域的专业数值模拟软件，在数值计算和分析方面有先天优势，但其建模功能受制于命令流式的手动建模方式，对于构建复杂模型过程复杂。ANSYS有限元程序具有强大的前处理功能，因此，借助ANSYS软件建立三维实体模型、划分网格单元，然后导入FLAC3D进行数值计算，将大大简化建模过程，提高数值模拟的效率和精度。

3.1　几何模型构建

ANSYS提供了2种建模方式，一是创建或导入实体模型，然后进行网格划分，并生成有限元模型；另一种是利用单元和节点编号，依靠ANSYS提供的命令流进行生成操作。本文选择前者，其建模过程可简述如下。

1)将前面处理好的关键岩层分界面iges文件导入WORKBENCH的Geometry之中，形成如图7(a)线框图，(b)面域图所示的面文件。

2)将相邻两面选取，点击Skin命令蒙皮成体，“面面构体”生成体单元，如图7(c)所示，完成后组成新的Part部分。

3)转入到“Model”模式，选择合理参数，对模型进行网格划分。

4)将Model中划分网格后的模型导入到Mechanical APDL模式。

5)然后利用河海大学郑文棠博士编写的ANSYS_TO_FLAC3D接口程序将ANSYS模型信息导出，另存为节点(Node)和单元(Element)2个文件，并经过转换文件转换为*.flac3d文件。FLAC3D通过Import Grid即可直接读取该文件，得到如图8所示的三维地质模型。

图7　ANSYS中由面构体形成的三维地质模型Fig.7　3D model generation using the face data

3.2　网格划分

矿山工程岩体应力具有三维特征，然而三维数值计算及成图复杂，对计算软件要求高，再加上反演区域范围大，划分单元多，给计算带来难度。为此，模型大小及网格划分尤为关键。

ANSYS提供了强大的网格划分功能，相比FLAC3D更为智能和便利。综合考虑有限元数值模拟的速度和构建效果，网格划分时选定的网格大小为80 m，最小网格为0.5 m，采用六面体网格划分法将计算模型进行划分，共划分得到114 240个单元，488 179个节点，网格划分完毕后的效果如图8所示。

图8　网格划分结果Fig.8　Mesh partition

3.3　FLAC3D中模型建立与数值计算

在ANSYS WORKBENCH中建立好三维地质模型后，对实体模型进行单元划分，将模型导入到ANSYS APDL中，然后利用河海大学郑文棠博士编写的ANSYS_TO_FLAC3D接口程序将ANSYS模型信息导出，另存为节点(Node)和单元(Element)两个文件，并经过转换文件转换为*.flac3d文件。FLAC3D通过Import Grid即可直接读取该文件，得到如图9所示的三维地质模型。

图9　有限元数值模型的建立与导入Fig.9　Establishment and loading of finite element model

根据实际情况对不同岩层分别赋予材料属性，并进行分组处理；然后依据前期分析确立相应的边界条件和计算参数，最后导入模型进行数值计算。

4　地应力反演结果展示与分析

以山东郓城煤矿为例，根据该矿前期测试的4个测点的地应力数据，通过建立三维数值计算模型，利用先进的数值模拟软件(FLAC3D)和应力反分析技术，反演和重构郓城煤矿一采区初始地应力场。最终获得的郓城煤矿初始地应力场回归模型，回归模型的相关性系数为88.26%，表明拟合效果较好。

图10所示为利用回归系数重构的初始地应力场SZZ的分布情况，取其中3组岩层上覆岩层、煤层和底部岩层展示。

图10　煤岩层SZZ应力分布规律Fig.10　SZZ stress distribution of coal seam and rock strata

图11　-800 m最大主应力SIG1分布云图Fig.11　Nephogram of maximum principal stress on -800 m level

为了更清晰展示郓城煤矿一采区区域地应力概况，将反演重构的煤层地应力场投影到CAD图上，得到如图11所示的最大主应力SIG1分布云图。可以看出，最大主应力与埋深和构造密切相关。图中深色区域为主应力较大区域，这些区域正好位于褶曲构造位置，由此可以推断，褶曲是引起区域内岩体应力高度集中的主因；同时，区域内的水平应力大于垂直应力。

利用上述建模方法，实现了对复杂矿山地质模型的快速构建，这为准确计算和展示地应力场分布规律提供了依据，为区域地应力场反演奠定了基础，同时也为矿井深部开采和采掘规划提供了依据。

5　结论

1)利用3DMINE，SUFER，ANSYS以及FLAC3D等软件联合建立大型精细三维有限元模型的方法，通过“点-线-面-体”的建模思路，实现了复杂三维地质模型的构建。

2)以山东郓城煤矿为例，通过构建了该矿的三维地质模型，结合现场实测测点地应力值，反演出该矿的地应力场分布，并以应力云图的方式进行直观展示，效果良好。

3)实际工程应用表明，该方法可实现大尺度、复杂三维地质模型的构建，有效解决了复杂模型构建效率低、网格划分难度大等问题，通过数据的逐步对接，减少了中间过程的人为干扰和误差，一定程度缩短了建模时间、提高了数值计算精度。

[1]康红普，林健，张晓. 深部矿井地应力测量方法研究与应用[J]. 岩石力学与工程学报，2007，26(5)：929-933.

KANG Hongpu, LIN Jian, ZHANG Xiao. Research and application of in-situ stress measurement in deep mines[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(5): 929-933.

[2]杨德林. 岩土工程问题的反演理论与工程实践[M]. 北京：科学出版社，1996.

[3]龚晓南．对岩土工程数值分析的几点思考[J]．岩土力学，2011，32(2)：321-325.

GONG Xiaonan. Reflections on numerical analysis of geotechnical engineering[J]. Rock and Soil Mechanics, 2011, 32(2): 321-325.

[4]朱志洁，张宏伟，韩军，等. 不同应力条件下巷道稳定性研究[J]. 中国安全生产科学技术，2015，11(11)：11-16.

ZHU Zhijie, ZHANG Hongwei, HAN Jun, et al. Research on stability of roadway under different conditions of stress field[J]. Journal of Safety Science and Technology, 2015, 11(11): 11-16.

[5]刘宁，朱维申，辛小丽. 双江口水电站初始地应力场反演回归分析[J]. 山东大学学报(工学版)，2008，38(6)：121-126.

LIU Ning, ZHU Weishen, XIN Xiaoli. Back regression analysis on the initial geostress field of the Shuangjiangkou hydropower station[J]. Journal of Shandong university: Engineering Science, 2008, 38(6): 121-126.

[6]郑文棠，徐卫亚，童富果，等．复杂边坡三维地质可视化和数值模型构建[J]．岩石力学与工程学报，2007，26(6)：1633-1644．

ZHENG Wentang, XU Weiya, TONG Fuguo, et al. 3D geological visualization and numerical modeling of complicated slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(6): 1633-1644.

[7]裴启涛，李海波，刘亚群，等. 复杂地质条件下坝区初始地应力场二次反演分析[J]. 岩石力学与工程学报，2014，33(S1)：2779-2785.

PEI Qitao, LI Haibo, LIU Yaqun, et al. Two-stage back analysis of initial geostress field of dam areas under complex geological conditions[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(S1): 2779-2785.

[8]崔芳鹏，胡瑞林，刘照连，等．基于Surfer平台的FLAC3D复杂三维地质建模研究[J]．工程地质学报，2008(5)：699-702.

CUI Fangpeng, HU Ruilin, LIU Zhaolian, et al. Surfer software platform based on complex three dimensional geological digital models for pre-processing of FLAC3D[J]. Journal of Engineering Geology, 2008(5): 699-702.

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

LIAO Qiulin, ZENG Qianbang, LIU Tong, et al. Automatic model generation of complex geologic body with FLAC3D based on ansys platform[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(6): 1010-1013.

[10]向鹏，纪洪广，邹静，等. 复杂矿山地质形体的三维数值网格划分自动化平台[J]. 岩土力学，2015，36(4)：1211-1216.

XIANG Peng, JI Hongguang, ZOU Jing, et al. An automation platform for generating the 3D numerical grids of mining fields with complex geomorphology[J]. Rock and Soil Mechanics, 2015, 36(4): 1211-1216.

[11]赵文，段魁，王彦东，等. 基于MATLAB、SURFER平台的FLAC3D复杂三维地质建模研究[J]. 工程地质学报，2016，24(s1)：208-213.

ZHAO Wen, DUAN Kui, WANG Yandong, et al. Study on complex three-dimensional geological modeling of Flac3D based on Matlab and Surfer platform[J]. Journal of Engineering Geology, 2016, 24(s1): 208-213.

[12]袁海平，王金安，蔡美峰. 复杂地应力场几何跨尺度反演与重构方法研究[J]. 采矿与安全工程学报，2011，28(4)：589-595.

YUAN Haiping, WANG Jinan, CAI Meifeng. Geometric cross-scale back-analysis and reconstruction menthod for complicated geo-stress field[J]. Journal of Mining & Safety Engineering, 2011, 28(4): 589-595.

[13]李翔，王金安，张少杰. 复杂地质体三维数值建模方法研究[J]. 西安科技大学学报，2012，32(6)：676-681.

LI Xiang, WANG Jinan, ZHANG Shaojie. 3D modeling method of complicated geological body[J]. Journal of Xi’an University of Science and Technology, 2012, 32(6): 676-681.

[14]王金安，李飞. 复杂地应力场反演优化算法及研究新进展[J]. 中国矿业大学学报，2015，44(2)：189-205.

WANG Jinan, LI Fei. Review of inverse optimal algorithm of in-situ stress and new achievement[J]. Journal of China University of Mining & Technology, 2015, 44(2): 189-205.

[15]赵自强，王学潮，随裕红. 南水北调西线工程深埋隧洞岩爆与地应力研究[J]. 华北水利水电学院学报，2002，23(1)：48-50.

ZHAO Ziqiang, WANG Xuechao, SUI Yuhong. The deep tunnel rock blasting and crustal stress research of west route project of water transfer from south to north[J]. Journal of North China Institute of Water Conservancy and Hyadroelectric Power, 2002, 23(1): 48-50.