开采沉陷FLAC3D参数化建模方法研究

陈元非，王 磊

(1. 池州学院 地理与规划学院，安徽 池州 247000；2. 安徽理工大学 空间信息与测绘工程学院，安徽 淮南 232001)

0 引言

FLAC3D数值模拟软件是美国明尼苏达Itasca软件公司开发的三维显式有限差分程序[1]，并首先将其应用到岩土工程力学计算中，现已成为岩土力学计算中重要的数值计算方法，在采矿、土木、水利工程等领域有着广泛的应用[1]。FLAC3D软件本构模型丰富，可以模拟各种弹塑性模型和开挖模型；具有强大的后处理功能，可以输出位移、应力、塑性区分布等各种信息，并且提供了内嵌的FISH语言，方便用户进行二次开发。但FLAC3D软件应用上也存在明显的短板，一是采用命令流的方式，二是在模型构建上的困难，一直是困扰科研工作者的难题之一[2]。

对此，国内外也有学者展开了复杂FLAC3D模型建模和预处理的研究。高永刚[2]通过分析AutoCAD二维平面图形与FLAC3D网格单元的关系，提出了基于AutoCAD快速建立FLAC3D网格模型的方法。李根[3]等利用ANSYS软件的强大建模功能，验证了ANSYS模型导入FLAC3D的具体操作可行性。罗周全[4]、林杭[5]提出基于SURPAC软件的直观、快速构建FLAC3D模型的技术。荆永滨[6]等利用TetGen网格生成工具实现了FLAC3D模型的可视化构建和交互式操作。冯艳顺[7]尝试了GIS结合EXCEL进行快速建模的方法，实现了矿区地层快速建模。刘秀军[8]利用GOCAD 建模软件，实现了将复杂三维地质模型导入FLAC3D软件中。可以看出，大量的复杂FLAC3D模型构建的研究基本是借助于第三方软件平台，因此用户需要在建模前对几种软件的使用较为熟悉，掌握建模方法并了解其与FLAC3D建模的区别。这无疑增加了额外的负担和使用FLAC3D软件的难度，打消了用户学习和使用软件的积极性[2]。

开采沉陷问题的FLAC3D建模具有它的特殊性，例如开采工作面一般较为规整，工作面采厚基本不变，且位于煤层中。因此可以尝试利用简单的六面体网格构建模型，而不必要引入第三方软件进行复杂的模型构建。为了满足开采沉陷问题研究中FLAC3D建模应用方面的具体需要，同时尽可能减小FLAC3D建模难度，降低使用门槛，本文提出一种参数化构建FLAC3D模型的方法。该方法利用大量实测钻孔数据，通过直接参数输入的方式，自动生成模型构建的命令流文件，避免了命令流代码编辑及后续建模上的困难，可以较为准确的模拟地下矿山开采实际的三维地层模型，对于倾斜煤层、变化厚度的岩层、山区地表等复杂三维模型的开采沉陷预测和规律的研究提供了方法，也对矿山开采沉陷虚拟仿真教学具有一定的参考价值。

1 FLAC3D建模方法

1.1 FLAC3D均匀格网建模方法

在利用FLAC3D建模进行矿区开采沉陷预计问题中，一般的会根据地表实测大量的钻孔数据，选择一个具有代表性的钻孔进行岩层简化合并，最终确定各个岩层的厚度。因此，在大量钻孔数据被遗弃、不同岩层被简化合并的基础上，构建出的FLAC3D地层模型往往不能准确模拟该地区实际地层情况，计算结果也因此具有较大偏差。

为了尽可能多的利用上每个钻孔的实测数据，结合开采沉陷研究领域的特点，将整个模型从平面上进行单元均匀格网化，即将模型底面进行均匀格网划分，通过插值计算得到各层格网节点的坐标，再进行组合建模的思路。利用单元网格模型采用brick六面体中8个角点定位的方式(如图1所示)，构建不同形状的单元格网模型。再利用单元格网模型之间的位置逻辑关系，进行组合计算和建模，如图2所示。

图2 相邻格网间的连接

在图1中，直接由8个角点确定一个单元brick六面体模型，它表示某一岩层的单元矩形格网模型。P0-P1P2P3方向为X、Y、Z轴方向。单元格网模型的尺寸可以设定，根据设置的整体模型XYZ方向的总尺寸和SIZE参数，确定各单元格网的尺寸。对于图1所示单元格网模型，设置得该岩层单元格网模型SIZE=1,1,3。

多个单元格网模型组合形成整体模型，各单元之间有其紧密连接的条件。如图2所示，单从X方向上看，单元1和单元2有4个点共用的约束条件，分别是[1-P1；2-P0]；[1-P4；2-P2]；[1-P6；2-P3]；[1-P7；2-P5]。因此Y方向上各个相邻单元之间的约束条件分别是[1-P2；3-P0]；[1-P4；3-P1]；[1-P5；3-P3]；[1-P7；3-P6]。从Z方向上，下一层岩层的顶面又是上一层岩层的底面，共用下一岩层的P3、P5、P6、P7四个点。因此，采用均匀格网的方法建模时，各节点XY坐标的计算较为简单，可以直接通过均匀格网的划分得到。Z坐标则可以采用多钻孔数据中的岩层厚度信息经过空间插值的方法计算。这样就可以构建由众多单元格网模型组合构成的FLAC3D整体模型。

1.2 网格节点岩层厚度数据插值

多钻孔数据的建模方法的思想是将整体模型格网化，分解成一个个便于建立的单元模型，再将单元模型组合成一个整体模型的思路。其中最为重要的工作是准确插值计算每个单元模型中各个节点的Z值。对于插值方法国内外已经有过大量研究。其中比较典型的插值方法有：克里金(Kriging)插值法、反距离加权插值法(IDW)、样条插值法(Spline)等。其中克里金插值法在地质探矿、地理空间插值、地学统计与建模中应用十分广泛[9,10]。因此，本研究采用克里金插值作为钻孔岩层数据的插值方法。

多钻孔数据建模工作主要是计算单元格网模型的8个角点的三维坐标，其难点在于插值各岩层上界面4点的Z坐标。通过单元格网的底面Z坐标，利用克里金插值方法插值出该节点处的各个岩层厚度，并从下往上相加即可。底面4个节点则可以根据相邻下一岩层单元模型的顶面推算，以此推算出所有节点的Z值。采用克里金插值方法插值的流程如图3所示。

图3 克里金插值法对任一单元节点P的层厚计算流程

1.3 参数化建模方法

均匀格网的建模方法可以实现基于多钻孔数据的岩层模型建模，但是整个建模过程十分繁琐，需要给每一岩层的每一个单元格网模型进行建模、组合、建立接触面、参数赋值、开挖、求解等过程，并且涉及到网格节点的岩层厚度插值计算等内容。因此仅仅利用手工计算编辑命令流的方式是不可能实现的。为了降低使用难度，结合EXCEL VBA程序编译功能，对程序功能进行封装，用户可以输入建模需要主要参数，从而直接生成FLAC3D建模命令流文件，实现建模计算功能的参数化。在利用程序生成命令流文件时，需要输入的主要参数如表1所示。

表1 程序设计输入参数内容

参数化建模过程直接利用命令流文件中网格坐标数据之间的逻辑联系，通过程序设计手段，生成完整的命令流文件。由于整个模型构建是由大量单元格网模型的组合形成，因此文件中会为每一单元格网模型编写一行由brick六面体构建的相应的命令流。最终完成整体模型命令流文件的生成。以平面5×5格网建模为例，其建模过程如图4所示。

图4 建模流程示意图

(1) 首先建立最底层初始单元格网模型a，单元网格模型8个角点坐标中，平面坐标由模型X、Y方向坐标与X、Y单元格数SIZE划分确定。模型Z坐标由钻孔数据中底层岩层厚度数据进行插值计算得到；

(2) 依次沿X、Y轴方向建立其它底层单元格网模型b～f；

(3) 继续依次沿Z轴方向建立其它岩层单元格网模型g。该层单元格网模型底部Z值与相邻下部岩层模型该点顶部Z值相等；顶部Z值由钻孔数据中各层岩层厚度数据进行插值，并从底部向上累计计算得到；

(4) 最终完成整体模型构建h。

2 工程案例

现模拟矿区计划开采一个走向长2000 m，倾向长300 m的工作面。在该区域范围内，共有岩层钻孔数据123个。该区域岩层从下至上分别为粉砂岩(底板)、煤层、中砂岩 、细砂岩 、表土层共5层，煤层厚度4 m。各岩层的力学参数可以参考以往历史数据。现可以根据需求构建一个平面大小为4000 m×4000 m区域的模型。模型开挖煤层中范围为1000 m

利用Excel vba编译的程序自动生成FLAC3D命令流文件，并用FLAC3D软件调用命令流文件，生成矿区该地区三维地层模型，如图5所示。可以看出，该地区地表有一定的起伏状态，地表点距离模型底部高度在170～500 m距离不等，坡度在1/10左右。为分析数据插值效果，对比了直接利用钻孔数据的Surfer8绘制的地表图，如图6所示。开挖工作面后最终计算的下沉和水平移动云图如图7所示。该模型在Core i7、32G内存的计算机上运行，代码自动生成和建模过程用约3 min，模型解算过程约用1 h。可以看出参数化建模方法建模过程所消耗的时间相对于整体解算过程消耗时间可以忽略，验证了该方法的效率。

图5 矿区FLAC3D三维模型

图6 程序插值地表(a)与Surfer8软件绘制的地表云图(b)对比

图7 模型开挖之后地表下沉(a)及水平移动云图(b、c)效果

从模型建立和开挖之后的变形结果来看，模型地表的三维结果与Surfer8基于钻孔数据的插值结果一致。开采变形运行结果沉降最大值约为3.5 m，最大水平移动值约为0.8 m，下沉系数和水平移动系数分别为0.88和0.22，符合一般矿区的经验，验证了该方法可靠性和有效性。

3 结论

(1) 基于开采沉陷预测领域FLAC3D建模问题进行分析，提出利用均匀格网分割、插值计算、组合为一体的参数化建模方法，实现了多钻孔实测数据的建模应用。可以较为准确的模拟地下矿山开采实际的三维地层模型，对于倾斜煤层、变化厚度的岩层、山区地表等复杂三维模型的开采沉陷预测和规律的研究提供了方法，也对矿山开采沉陷虚拟仿真教学具有一定的参考价值。

(2) 针对网格节点坐标计算和组合问题采用克里金插值方法插值岩层厚度数据，并通过编程实现了克里金插值、参数化建模功能，避免了人工处理大量钻孔数据编辑命令流文件的困难。

(3) 利用矿区实测钻孔数据进行参数化建模应用，通过分析三维模型和变形模拟结果，验证了该方法的可靠性和有效性。