复杂开挖过程FLAC3D力学仿真代码生成系统研究

马长年，徐国元,，江文武，刘晓明

（1. 中南大学 资源与安全工程学院，长沙 430083；2. 华南理工大学 土木与交通学院，广州 510641）

1 引 言

伴随着人类对自然资源、能源需求的不断增加，各类岩土工程、采矿工程规模越来越大，岩土工程结构和开挖过程越来越复杂，对所涉及结构和过程进行更为贴近实际的三维力学分析、优化结构体力学结构和开挖过程、保证岩土工程的稳定性和安全性有着重要意义。FLAC3D是美国Itasca Consulting Group Inc.开发的三维有限差分力学计算程序，是一种基于拉格朗日算法的数值计算方法。FLAC3D提供了7种材料本构模型，适用于模拟材料的大变形和扭曲转动，已被广泛应用于模拟计算岩土类工程地质材料的力学行为[1－5]。

FLAC3D提供了一组命令集和内嵌式解释语言Fish，用于编写模拟各类连续体力学过程程序，但对于复杂形态工程几何体、多属性、几十万至几百万单元、大量结构体和复杂开挖过程，由用户手工编写程序代码对上述过程进行三维力学分析存在诸多困难，主要表现在：（1）模型的建立只能用数据文件来实现，不能直观显示，建模过程可视化程度低；（2）对于复杂的大型工程地质体建模，需要控制各个边界节点的坐标数据，容易出错，检查也不方便[3－5]；（3）对于大量单元、结构体对其特定属性进行检索、查询、分类和赋值等操作困难；（4）手工编制上述大规模FLAC3D程序效率低、耗时多、出错概率大、查错纠错困难，常导致程序运行不收敛、结果失真等异常。

本文所述 FLAC3D力学仿真代码生成系统是对解决此类问题进行的研究与尝试。

2 复杂开挖过程剖分

2.1 剖分方法分类

任何复杂开挖过程都是由一些具体的、特定的、可划分成模块或阶段的过程构成的，通过每一个模块或阶段，开挖过程模拟可实现整个开挖过程的力学仿真。以下按4种不同方法对复杂开挖过程进行剖分。

（1）按照坐标（X、Y、Z及其组合）、平面（或曲面）相对位置（上下/前后/左右）、边界（内/外）对模型进行剖分；

（2）按一定的开挖顺序（串行或并行）对一组或多组单元进行剖分；

（3）按不同的结构功能可对开挖过程进行多种方式的剖分，如隧洞的开挖过程可按结构功能剖分为：上拱、侧墙、底拱、核部、锚杆和衬砌等结构；

（4）按不同的岩性、开挖状态（开挖/未开挖）、充填状态、矿体/岩石、原岩/破碎带等单元属性进行剖分。

2.2 复杂开挖过程剖分实例

（1）按矿体/岩石属性剖分模型

[4－6]的方法，针对某采用下向进路胶结充填采矿法矿山，本文利用Surpac软件为该矿矿体和地表地形按照模型整体要求的范围分别建立3 D和DTM几何模型[6-7]，所建模型如图1所示。

图1 矿山三维几何模型Fig.1 3D model of mine

在此基础上建立模型的Surpac块体模型，扩展其属性得到结构如表1所示的包含单元质心坐标数据的块体单元数据库。

表1 块体单元数据库基本结构Table1 The basic construction of zone table

可通过矿体块体单元约束条件向块体单元数据库岩性属性赋值，矿体内赋1，矿体外赋0。矿体块体单元约束中的块体尺寸必须与FLAC3D模型的单元尺寸保持一致。由Surpac块体模型质心坐标数据可生成六面体单元的FLAC3D模型[4]。

（2）按开挖顺序剖分模型

参照矿山实际生产情况，按矿体高程以 20 m划分分段；在分段内沿矿体走向按整数勘探线每100 m划分为一个盘区；进路规格与分段高保持一致，进路布置有沿走向和垂直走向2种方式。

（3）按结构体剖分模型

进路中布设的钢筋（吊筋和水平筋）可由FLAC3D所提供的锚索结构体模拟，系统默认所有回采进路都布设钢筋，所有矿体内单元和部分围岩中单元将与一条或多条锚索结构体关联。这样模型单元被划分为布设锚索和不布设锚索2类。

3 FLAC3D力学仿真代码生成系统开发

3.1 模型单元操作

（1）对模型单元直接操作的意义

大型 FLAC3D力学模型通常是由几万至几百万甚至更多的单元组成的，由于地质工程体的原始应力状态、材料属性、模型几何形态、开挖体几何形态、开挖顺序和步骤、工程结构体等的多种多样，每一个单元都有可能具有自身独特的属性。对任一单元进行直接操作，主要涉及对特定单元进行定位、捕捉和重组的问题，如此才可能对各类复杂开挖过程进行有效的力学模拟与仿真。

（2）FLAC3D操作模型单元的方法[2]

FLAC3D对模型单元定位的方法归纳起来有 2种：RANGE命令和单元操作函数。

RANGE命令有一组子关键字，可按环、柱、组、方向、平面相对位置、模型类型、坐标、单元号等多种方式确定单元位置或单元范围。其中按坐标和单元号操作单元是最为基础的。

单元操作函数有：

（3）基于单元数据库的模型单元操作方法

将模型所有单元的质心坐标、岩性、开挖标志、单元上施加的结构体等属性（还可以增加其他属性）通过Surpac形成块体单元数据库，其单元数据表基本结构如表1所示。

通过数据库查询语言可对单元一个或多个属性进行组合检索，将单元质心坐标传递给FLAC3D的单元捕捉函数，如gp_near(x,y,z),z_near(x,y,z)等，再利用函数z_id(p_z)得到单元在模型中的单元号。这样可利用单元号对需要相同操作的单元通过重组进行成批操作和处理。单元的重组和批处理操作可利用FLAC3D中数组、函数和循环语句完成。

3.2 应用程序变量与FLAC3D关键字、变量、函数的映射和参数传递

3.2.1 Fish语言和FLAC3D命令特点[2]

（1）Fish语言是内嵌于FLAC3D内部的语言，它可以让服务用户自己定义变量和函数；

（2）Fish共有4种数据类型：整型、浮点、字符串和指针；

（3）Fish函数由一对关键字Define与End定义；

（4）在Fish语言中所有的变量和函数都是全局的（与Basic相似），函数没有参数；

（5）一条完整的Fish语句占用一行，无法续行。如果公式太长，需要用临时变量将其分开；

（6）Fish包含条件分支、条件选择、循环语句等程序控制语句；

（7）Fish预定义了一些函数和变量，使用此类变量和相关函数，可获得所有节点、单元或全局空间坐标点在内存中的地址，并对其进行操作；

（8）函数可重定义，重定义后，原先的函数将不能再被调用，但原函数所形成的内存空间仍存在，其中的内存变量不变，数组不可重定义。

以上这些Fish语言和FLAC3D命令特点和功能是开发代码生成系统所必须了解和注意的。

3.2.2 字典变量与FLAC3D关键字映射

为实现由 FLAC3D代码生成系统自动生成FLAC3D程序，首先要建立应用程序中变量与FLAC3D语句关键字的映射关系，即翻译字典。应用程序由 VC++开发，在此将专门用于处理这一映射关系的字符串变量称为字典变量，数据类型为CString，字典变量映射格式如表2所示。这样在应用程序中建立起字典变量与所有FLAC3D关键字之间一一对应的关系，这是代码生成系统自动生成FLAC3D应用程序代码的基本前提。

表2 字典变量与FLAC3D关键字映射表Table2 The mapping table between dictionary variables and FLAC3Dkeywords

3.2.3 传递程序变量值到FLAC3D变量、数组和函数

下面以布设2条锚索的程序代码为例来说明如何将程序变量值传递到FLAC3D变量、数组和函数，并形成2条锚索布设的FLAC3D程序代码。

VC++源程序代码片段

至此，我们建立了代码生成系统与FLAC3D关键字、变量、数组、函数和控制语句之间的对应关系，这是应用程序生成FLAC3D程序代码的基础。

3.3 仿真代码生成系统基本原理

（1）仿真代码生成系统的输入和输出

仿真代码生成系统的输入是由 Surpac块体模型导出的包含单元质心数据、单元属性及单元约束状态的Access数据库。仿真代码生成系统的输出是按用户的要求，由系统生成的FLAC3D程序代码。

（2）仿真代码生成系统的工作原理

仿真代码生成系统程序通过 Ado访问 Surpac块体模型导出的块体单元数据库，利用数据库所提供的函数可以方便地对其中的数据进行多种类型的检索和操作。同时程序提供一组交互式表单与用户进行交互，由用户对一类特定岩体回采开挖过程所需步骤和参数进行选择和确认。系统接收这些步骤和参数后生成*.txt格式的FLAC3D仿真程序代码，同时修改块体单元数据库中相应单元属性值。

3.4 系统结构与功能设计

（1）系统结构

①单元数据由数据库结构组织；②应用程序由单文档应用程序框架来实现；③属性表单向导结构：应用程序将采用向导方式按生产过程划分的步骤，一步一步引导用户完成满足要求的FLAC3D代码。每一步采用属性表单的方式，为用户提供多个选项，以确定 FLAC3D力学仿真程序的各关键属性和参数。

（3）功能设计

本文所述 FLAC3D仿真代码生成系统由 VC++开发，系统具体功能是为模拟某采用下向充填法矿山的实际生产过程力学仿真而设计的，主要功能如图2所示。

图2 仿真代码生成系统功能设计Fig.2 The function design of system of generating code

各页面主要完成以下功能：

①主页面为一单文档应用程序，主要功能包含在主菜单“仿真代码生成”菜单中，如图 2(a)所示。打开“仿真代码生成菜单”，程序将引导用户逐步完成单元数据库调用、分段水平划分、盘区划分及进路回采的 FLAC3D力学仿真代码，并以 FLAC3D可识别的*.txt文件格式输出到用户指定位置。

②点击“矿体分段划分”按钮，程序将要求用户调入矿体单元数据库。通过访问块体单元数据库，程序可按分段高度从数据库中提取矿山开采分段的高程数据，本步骤是在垂直方向对矿体进行划分，用户可选择将开采分段，界面功能设计如图2(b)，该步骤同时向程序变量提供单元的Z坐标。

③本步骤是对上一步所选择的开采分段进行盘区划分，通常只对矿体部分进行划分。盘区宽度默认为100 m，盘区编号从西到东。盘区是按生产实际情况划分的，盘区划分有利于实现不同开采顺序的力学模拟。界面功能设计如图2(c)，该步骤向程序变量提供盘区内单元X、Y坐标。

④由用户确定盘区内进路布置的方向，选择待回采进路，确定充填进路内是否布设钢筋等参数（在本实例中钢筋的力学参数以现场实际使用值为默认值）。本步骤主要用来实现生产实际中上下分段（分层）进路交错布置的仿真，在一个盘区内可模拟隔一采一、隔二采一等回采方式。界面功能设计如图2(d)，该步骤向程序变量传递待采进路所有单元质心坐标及钢筋（锚索）的起点、终点单元的质心坐标，同时将单元数据库中相应单元的开挖标志和充填标志置为1。

⑤用户可指定仿真代码的存储位置，每个盘区内的进路可按一定回采顺序生成一个*.txt文件。完成后程序将进入下一个选单，点击回采新盘区可对新的盘区进行回采，点击回采新分段可对新的分段进行回采，直到按用户的方案回采完毕，点击完成按钮退出代码生成器，其界面功能设计如图 2(e)、(f)。

3.5 类结构

（1）基本类

单文档 MFC应用程序都具有的 4个类：CCodeGenApp、CCodeGenDoc、CCodeGenView、CMainFrame。

（2）表单类

该类为程序功能的实现提供与用户交互的界面（见表 3），用户通过这些界面的选项将自己的意图一步步传递给程序，程序从用户获得所需的参数，最终由程序自动地形成FLAC3D力学仿真代码。

（3）功能类

CDBP类：该类是用来访问数据库和进行数据处理的类，CDBP通过接收 CProPageLevel、CProPageSection、CProPageMining、CProPageReturn和CProSheet等类传递来的消息和参数，对Surpac块体单元数据库进行访问和检索。

表3 类的种类与功能Table3 The type and function of class

3.6 仿真程序运算结果

应用本系统成功地生成了下向进路式充填采矿法矿山生产过程的FLAC3D力学仿真代码，可实现垂直方向、水平方向的开采顺序，不同进路布置方式，充填体内钢筋布置等主要生产过程的力学仿真。程序中所采用的有关矿山原岩应力场、矿岩力学参数、开采工艺[8－11]基本与矿山实际保持一致。模型共有139626个单元，锚索结构体1256个，是1个较为复杂的三维力学模型，所生成程序代码调入FLAC3D运行正常。矿体上部5个分段回采充填后，矿体附近第1主应力分布和钢筋（锚索）受力状态仿真结果见图3，图3（b）中锚索浅色表示受压，深色表示受拉，粗细表示受力大小。

图3 部分仿真结果图Fig.3 The part of simulating results

4 结 论

（1）FLAC3D代码生成系统从分析 FLAC3D命令集和Fish语言特点入手，着重解决了模型单元质心定位、单元捕捉和重组，应用程序变量到FLAC3D变量、函数、语句、数据结构之间映射等问题，为编制FLAC3D程序提供了一种新的方法，能够大幅提高编制复杂开挖过程FLAC3D程序代码的效率，降低代码编制出错率。

（2）GCS系统可对大量模型单元进行直接和有效的操作，具有良好的适应性和可扩展性。

（3）在建立块体单元数据库过程中利用了矿用地质软件Surpac的许多功能，架起了地质三维几何造型软件与三维力学计算软件之间的桥梁。

参考文献

[1]谢和平,陈忠辉. 岩石力学[M]. 北京: 科学出版社,2004.

[2]Itasca Consulting Group. Fast Lagrangian analysis of continua in three dimensions[M]. America: Itasca Consulting Group.,1997.

[3]廖秋林,曾钱帮. 基于 ANSYS平台复杂地质体FLAC3D模型的自动生成[J]. 岩石力学与工程学报,2005,24(6): 1010－1013.LIAO Qiu-lin,ZENG Qian-bang. Automatic model generation of complex geologic body with FLAC3Dbased on ANSYS platform[J]. Chinese Journal of Rock Mechanics and Engineering,2005,24(6): 1010－1013.

[4]罗周全,吴亚斌,刘晓明,等. 基于Surpac的复杂地质体 FLAC3D模型生成技术[J]. 岩土力学,2008,29(5):1334－1338.LUO Zhou-quan,WU Ya-bin,LIU Xiao-ming,et al.FLAC3Dmodeling for complex geologic body based on Surpac[J]. Rock and Soil Mechanics,2008,29(5): 1334－1338.

[5]何忠明,彭振斌,曹平,等. 双层空区开挖顶板稳定性的 FLAC3D数值分析[J]. 中南大学学报(自然科学版),2009,40(4): 1066－1071.HE Zhong-ming,PENG Zhen-bin,CAO Ping,et a1.Numerical analysis for roof stability of double gob area after excavation by FLAC3D[J]. Journal of Central South University (Science and Technology),2009,40(4): 1066－1071.

[6]Surpac Software International国际软件公司.SURPAC Vision软件用户使用手册(第四版)[M]. [s. l.]: SURPAC Software International国际软件公司,2000.

[7]刘晓明,罗周全,杨 彪,等. 复杂矿区三维地质可视化及数值模型构建[J]. 岩土力学,2010,31(12): 4006－4011.LIU Xiao-ming,LUO Zhou-quan,YANG Biao,et al.Numerical modeling and geological body visualization for complex mine [J]. Rock and Soil Mechanics,2010,31(12): 4006－4011.

[8]路德维格,宋恕夏,田永妥.中国-瑞典关于金川二矿区岩石力学研究报告[R]. 北京: 北京科技大学,1998.

[9]陈俊彦,N G W,库克,等. 金川镍矿二矿区采矿方法的岩石力学研究报告[R]. 长沙: 中南矿冶学院,1986.

[10]廖椿庭,施兆贤. 金川矿区原岩应力实测及在矿山设计中的应用[J]. 岩石力学与工程学报,1983,2(1): 103－112.LIAO Chun-ting,SHI Zhao-xian. In-situ stress measurements and their application to engineering design in the Jinchuan mine[J]. Chinese Journal of Rock Mechanics and Engineering,1983,2(1): 103－112.

[11]吴满路,马宇,瘳椿庭,等. 金川二矿深部1000 m中段地应力测量及应力状态研究[J]. 岩石力学与工程学报,2008,27(增刊2): 3785－3790.WU Man-lu,MA Yu,LIAO Chun-ting,et al. Study of recent state of stress in depth 1000 m of Jinchuan mine[J].Chinese Journal of Rock Mechanics and Engineering,2008,27(Supp.2): 3785－3790.