基于Midas/GTS的FLAC3D边坡建模技术及工程应用

钟志辉, 刘祚秋, 杨光华,，张玉成

（1.华南理工大学 土木与交通学院，广东 广州 510641；2.中山大学 工学院，广东 广州 510275；3.广东省水利水电科学研究院，广东 广州 510610）

0 前 言

FLAC3D是一款专业的岩土工程分析软件，其计算功能十分强大。作为有限差分软件，其采用的“混合离散法”来模拟材料的塑性流动和破坏[1]，比有限元采用的“离散集成法”更为准确和合理。然而FLAC3D的前处理功能较弱，尽管其提供了12种基本网格单元，可以通过连接、组合等命令或者内置的FISH语言建立网格模型，但对于复杂的二维或三维模型的建立仍然十分困难。目前工程技术人员更多的偏向于“（快捷）建立几何模型—（智能）划分模型网格”的建模思路，因为这有助于提高建模的效率和模型网格的质量，对模型的调试也十分方便高效。

针对 FLAC3D的前处理缺点，国内外不少学者研究了如何把模型网格从其他前处理功能强大的软件导入FLAC3D的技术方法[2-6]。本文探讨如何基于有限元软件Midas/GTS来实现FLAC3D的复杂建模，主要着重于二维边坡的建模技术，介绍 Midas/GTS与 FLAC3D网格之间的转换关系和建模的技巧，最后将此建模技术应用于乐昌峡某边坡的稳定性分析。

2 Midas/GTS与FLAC3D的模型转换方法[7-9]

2.1 Midas/GTS的建模技术

Midas/GTS是专门为岩土工程分析而设计的有限元软件，拥有强大的几何建模和网格划分技术，用户能在 Midas/GTS中快捷方便地建立复杂的网格模型，对于建立二维的边坡模型尤其方便。一般Midas/GTS建立边坡模型的过程如下：

(1) 通过点—线—面（—体）的方法建立边坡几何模型。对于二维模型，通常先根据边坡的剖面信息在AutoCAD中生成模型的线元素，因为在AutoCAD中能十分方便地建立和修改边坡的剖面轮廓、地层分布、地下水位线等几何信息，然后保存成dxf文件再导入 Midas/GTS，生成几何平面（组）；对于三维模型，还需要通过面（组）生成实体，然后对实体进行平移、嵌入、切割或布尔运算来生成更为复杂的三维实体（组）。Midas/GTS对几何元素的编辑提供预览、撤销、重做等人性化操作，大大提高建模速度，降低出错概率。

(2) 布置网格种子点。即在生成网格之前利用网格尺寸控制命令事先定义要生成网格的对象线的分割单元大小，这样能够使得不同面（或体）的相邻边（或面）的网格大小一致，以保证相邻边（或面）的节点耦合。

(3) 分网格。对于平面问题，Midas/GTS提供（高次）三角形和（高次）四边形单元，内置循环网格法、格栅网格法和德劳内网格法，三种方法均能划分高质量的单元，由于 FLAC3D对网格质量的要求较高，宜尽量采用规则的brick网格（六面体），因此对于较为复杂的边坡模型，建议采用格栅网格法以“四边形+ 三角形”的方式生成混合网格，这种划分方式使得绝大部分为均匀的四边形网格，只有在局部不规则的区域用三角形过渡。最后通过编制程序转换成FLAC3D中的“brick + wedge”混合网格。对于三维问题，一般采用四面体和六面体划分网格。

(4) 检查网格。Midas/GTS的“检查网格”菜单能检查网格的自由边（面），一般采用此功能来检查模型内部是否存在节点不耦合的情况；“检查网格质量”菜单则可以对单元的纵横比、锥度、扭曲、雅克比比率等进行检查，便于发现质量不好的单元。

2.2 Midas/GTS与FLAC3D的网格数据格式

Midas/GTS与FLAC3D有各自的数据结构，根据表 1中的节点对应关系就能实现从 Midas/GTS到FLAC3D的网格转换。目前边坡工程一般等效为平面应变问题或假三维问题分析，此时首先在Midas/GTS中建立平面网格模型，然后对平面单元进行纵向扩展（一般扩展单位长度）以生成FLAC3D的三维网格。如表1所示，平面三角形单元扩展为wedge网格、平面四边形单元扩展为brick网格。对于三维的边坡问题，Midas/GTS中的四面体、楔形体、六面体能直接转换为FLAC3D的tetrahedron网格、wedge网格、brick网格。

表1 Midas/GTS与FLAC3D网格节点对照表

2.3 Midas/GTS模型导入FLAC3D

Midas/GTS能够提供网格模型的节点和单元信息，通过节点表格和单元表格查询。节点表格的表头格式为“节点号，坐标系，X坐标，Y坐标，Z坐标”，每一行代表一个节点信息，单元表格的表头格式为“单元号，单元类型，属性，属性类型，材料，特性，节点1，节点2，节点3……节点N”，每一行代表一个单元信息，其中可以通过指定不同的属性或材料以实现FLAC3D中的分组（GROUP）功能。表格的内容均可复制到Excel表格或文本上以便进行操作。

FLAC3D的网格数据则是以点（GRIDPOINT）、单元（ZONE）、组（GROUP）的形式保存的，保存的文件以*.flac3d为后缀，完整的数据结构如下：

* GRIDPOINTS 节点表头

G 1 0.0 0.0 0.0 节点1坐标信息

G 2 1.0 0.0 0.0 节点2坐标信息

……

* ZONE 单元表头

Z B8 1 1 2 3 4 5 6 7 8 单元1的信息

随着城市轨道交通系统的迅速发展，地铁盾构隧道在运营期将不可避免地出现一系列病害现象，这将影响到地铁运营的安全性、经济性和耐久性。目前，国内外学者对隧道结构的病害问题已有大量研究，但大多是集中于对山岭、公路、铁路隧道结构病害的研究，对盾构隧道结构的病害成因及治理措施的研究相对较少[1]，因而对盾构隧道常见的病害原因及治理措施展开进一步研究显得尤为重要。

……

* GROUPS 组表头

ZGROUP soil 标记组soil

1 单元1属于组soil

……

ZGROUP rock 标记组rock

10 单元10属于组rock

……

* GRIDPOINTS为节点的表头，其下面的行是节点信息，格式为“G 节点号 X坐标 Y坐标 Z坐标”；* ZONE为单元的表头，其下面的行是单元信息，格式为“Z 网格类型 网格编号 节点 1 节点 2 节点3……”；* GROUPS为组的表头，ZGROUP标记一个组，一个组标记后面的行是单元编号，表示该编号的单元都属于这个组。

可见，根据表1的对应关系，将Midas/GTS的网格数据转换成 FLAC3D的网格数据文件，利用FLAC3D的impgrid命令读入数据文件即可生成网格模型。

3 三角形与四边形单元的精度比较

多数边坡工程按平面应变问题进行分析，因此一般在 Midas/GTS中采用平面三角形或四边形单元划分网格，然后转换成FLAC3D的brick或wedge网格。但采用三角形单元生成的网格的计算精度比四边形单元的要低，现选用澳大利亚计算机应用协会(ACADS)发布的一个边坡稳定分析考题来说明。该边坡模型见图 1，侧压力系数K0为 0.65，土体参数见表2，ACADS给出的安全系数为1.00。模型底部固定，左右侧边按σh=K0γh 施加侧压力。

图1 ACADS考题的计算模型Fig.1 Computational model of ACADS exercise.

表2 ACADS考题的计算参数

在 Midas/GTS中分别采用三角形和四边形进行网格划分，单元尺寸统一取0.5，网格分布分别见图2。导入 FLAC3D后采用内置的强度折减法求解安全系数，结果见表3。

可见对于相同的单元尺寸，即使三角形的单元数约为四边形单元的2倍，采用四边形单元的计算精度仍比三角形单元要高，求出的安全系数等于精确值1.00，而三角形单元的计算结果有 3%的误差。说明FLAC3D的计算精度对网格有很强的依赖性，一般认为均匀的brick网格较为适合有限差分法。因此建议尽量采用四边形单元进行网格划分；而对于较为复杂的边坡模型，建议采用格栅网格法以“四边形 + 三角形”的方式生成混合网格，这种划分方式使得绝大部分为四边形单元，而局部不规则的区域用三角形单元过渡，从而保证计算精度。

图2 三角形和四边形网格划分方式Fig.2 Meshing by quadrilateral and triangular manner.

表3 三角形与四边形单元的计算结果

4 工程应用

为了验证由 Midas/GTS生成模型网格然后导入到 FLAC3D的可行性，现采用这种建模技术对广东省乐昌峡的某临水边坡的一个剖面进行建模并分析其稳定性。该边坡于2005年在长时间大暴雨影响下，边坡体周边出现裂缝和滑动，表明其处于极限稳定状态。该边坡的土层分布较为复杂，如图3所示，地下水位线经过全风化土层，计算时水位线以下的区域按静水压力施加孔隙水压力。各土层的参数见表4，其中中风化土采用线弹性模型，其余土层采用Mohr-Coulomb理想弹塑性模型。

模型首先在Midas/GTS中采用格栅网格法以“四边形 + 三角形”的方式生成网格，共7 480个网格单元，其中四边形7 306个、三角形174个，然后把网格导入FLAC3D，图3给出了FLAC3D中的局部网格情况。

图3 乐昌峡边坡的土层分布及局部网格Fig.3 Distribution of soil layers and local grids on Lechangxia slope.

表4 乐昌峡边坡的土层参数[10]

采用文献[11]提出的局部强度折减法对该边坡进行稳定性分析，即仅对滑带土进行强度折减，直到边坡监测点的位移出现突变，此时的折减系数视为安全系数。局部强度折减法得出的滑动面见图4，基本与滑带重合。另外在坡面设置了位移监控点A、B（图4），其折减系数-水平位移曲线如图 5所示，当折减系数为1.06时监测点的水平位移发生突变，所以该边坡的安全系数为1.06，表明该边坡濒临失稳，与实际的边坡稳定状态十分接近。

图4 边坡滑动面位置Fig.4 Location of sliding surface on the slope.

图5 监测点折减系数-水平位移曲线Fig.5 Relation curves of reduced factor and horizontal displacement at monitoring points.

5 结 语

基于 Midas/GTS的强大建模功能，通过编制Midas/GTS与 FLAC3D的模型转换接口，从Midas/GTS导入网格到FLAC3D中，从而实现复杂边坡模型高效、高质量的建立，解决了 FLAC3D前处理的不足和困难，为 FLAC3D构建复杂的工程模型提供新的思路。

此外，针对二维复杂边坡的问题，本文建议在Midas/GTS中采用“四边形 + 三角形”的方式生成混合网格，然后导入FLAC3D中形成“brick + wedge”混合网格进行分析，以保证计算精度。

[1]Marti Joaquin, Cundall Peter. Mixed discretization procedure for accurate modelling of plastic collapse[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1982, 6(1):129-139.

[2]王树仁，张海清. MIDAS/GTS－FLAC3D耦合建模新方法及其应用[J]. 土木建筑与环境工程，2010,(1)：12-17.

[3]丁秀美，黄润秋，刘光士. FLAC-3D前处理程序开发及其工程应用[J]. 地质灾害与环境保护，2004，15(2)：68-73.

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

[5]罗周全，吴哑斌，刘晓明，等.基于SURPAC的复杂地质体FLAC3D模型生成技术[J].岩土力学，2008，29(5)：1334-1338.

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

[7]Itasca Consulting Group.Fast Lagrange Analysis of continua in 3 dimensions, user’s manual[M]. MN, U.S.A. Minneapolis: Itasca Consulting Group, 2005.

[8]MIDAS Information Technology Co, Ltd. Manual of Midas/GTS[M].South Korea: MIDAS Information Technology Co, Ltd, 2005.

[9]陈育民，徐鼎平. FLAC/FLAC3D基础与工程实例[M]. 北京：中国水利水电出版社，2008：267-269.

[10]张有祥. 库岸边坡稳定性及抗滑桩加固研究[D]. 武汉：武汉大学，2009.

[11]杨光华，钟志辉，张玉成，等.用局部强度折减法进行边坡稳定性分析[J].岩土力学，2010，31(增刊2)：53-58.