多介质复杂区域四边形网格自动剖分算法及应用

付成华，周洪波

(1.西华大学能源与环境学院，成都 610039;2.二滩水电开发有限责任公司，成都 610051)

1 概述

有限元计算程序不断涌现出来并成功应用于实际工程中，围绕有限元计算程序集成开发的软件也越来越多。在有限元分析中，网格生成技术一直是制约科学和工程计算精度和效率的一个重要因素，所以根据计算区域的几何结构信息，利用计算机进行自动剖分，产生有限元计算所必须的网格数据，进一步提高有限元数值计算的效率，加大有限单元法及其集成软件的推广应用是非常必要的。

本文提出一种四边形网格自动剖分的改进行波算法。无需设置背景网格，将多介质复杂区域分成多个封闭子域，先生成平面三角形网格，优化后通过合并、分解生成四边形单元，再次对四边形单元优化处理即可获得过渡光滑且质量良好的四边形网格，从而实现多介质复杂区域的四边形网格自动剖分。最后通过应用于某铁路路堤断面和瀑布沟电站厂房2#机组断面计算域的网格剖分，验证了该方法的可行性及可靠性。

2 算法步骤

行波法最初由Lo提出并用于平面三角形网格的自动生成，后来通过引入背景网格来控制内部节点的生成，并推广到任意形状三维区域的四面体网格生成。行波法主要从区域边界开始，定义一个“前沿”，在前沿上满足一定条件的地方生成一系列网格，同时更新“前沿”，如此不断地重复循环，直到前沿缩为一点，结束网格的生成。行波法剖分网格的自动化程度高，对复杂几何形状与边界的二维(或三维)区域中的三角形(或四面体)网格生成具有很大的灵活性及可靠性［1－2］。

本文改进行波法是基于计算域AutoCAD模型和程序语言实现的，四边形网格自动剖分的具体步骤如下:

(1)将多介质复杂区域看成是多个封闭子区域的组合，利用编制的AutoCAD应用程序ptline.lsp从计算域的AutoCAD图中提取各封闭子区域边界的点线信息。

(2)离散区域边界，以离散的各子域边界作为“前沿”，利用改进行波法在各个子域内生成三角形网格，并进行三角形网格质量评价和优化。

(3)分子域将2个三角形合并为1个四边形，合并后的四边形分解为2个或4个小四边形，未合并三角形分解为3个四边形。

(4)对四边形网格进行质量评价及优化，实现整个计算域四边形网格的自动剖分。

3 关键技术改进

3.1 子域内部节点生成

将多介质复杂区域看作多个封闭子域进行剖分，每个子域在剖分过程中又被分为已剖分域和未剖分域，封闭域内部节点的生成步骤(图1)如下:

(1)根据网格密度控制参数离散各子域边界线得到各子域的初始生成波。

(2)从子域生成波中选取内角最小的点，作为新三角形的基点。

若基点A内角小于等于90°，则与基点相邻的2点B，C直接作为新三角形的顶点，生成三角形ABC，记录子域未剖分域，更新子域生成波。

若基点B内角＞90°，在∠DBC平分线上取点P，LBP=(LBD+LBC)/2，且P点在子域内部。若P点不在子域内，则取∠DBC平分线与子域边界的交点为P。判别∠FDP和∠PCE，若∠FDP、∠PCE均≥30°，则P点作为新生成的内部节点，生成三角形 BDP 和BPC;若∠FDP＜30°、∠PCE≥30°，则直接生成三角形BDF和BFC;若∠FDP≥30°、∠PCE＜30°，则直接生成三角形BDE和BEC。

(3)记录子域未剖分域，更新子域生成波。重复(2)，直至生成波为空。

图1 内部节点生成示意图Fig.1 Schematic diagram of the generation of internal nodes

3.2 网格密度控制

网格密度指单位面积内单元的平均个数，这里网格密度控制采用总体密度控制和局部过度控制相结合的方式实现。

(1)总体密度控制。由用户根据自身需要设定待剖分区域网格总体密度的控制参数，参数在(0，1)之间取值，数值越小说明网格越密。

(2)局部过渡控制。各子域按面积由小到大排序，依次离散各子域的边界线。如图2，离散子域边界线AB时，以2端点A，B为垂足分别向子域内做垂线AI和BJ，与子域边界的交点为I和J。A端离散长度l取总体控制长度、相邻边界线FA的长度 LFA和LAI三者中的最小者。

图2 边界离散控制示意图Fig.2 Schematic diagram of border discrete control

3.3 网格质量评价及优化

网格自动剖分过程中，单元质量控制非常重要。这里通过建立网格质量评价体系和多重优化措施来保证最终获得良好的四边形网格。

3.3.1 网格质量评价

(1)三角形网格质量评价。三角形单元质量主要通过单元内角来判断。三角形单元内角介于45°和135°之间，质量好;内角都介于30°和150°之间，质量一般;若最小内角＜30°或最大内角＞150°，质量较差［1］。

(2)四边形网格质量评价。按逆时针走向的四边形ABCD被对角线AC分成2个三角形ABC和ACD，由公式(1)计算其 α值分别为和;同时四边形ABCD又可被对角线BD分成2个三角形ABD和BCD，其α值分别为和。将4个三角形的α 值，，按降序重新排列成为α1，α2，α3，α4。则四边形 ABCD的质量判断采用文献［3－4］准则:

矩形的β值为1，β值越小，四边形质量越差。

3.3.2 网格优化

(1)三角形网格优化。本算法中三角形单元优化采用文献［1］中的优化方式:①删除节点;②对角线调整;③拉普拉斯优化。

(2)四边形网格优化。对网格质量较差的四边形采用拉普拉斯光顺和内部调整进行优化［1］。

4 应用实例

(1)某铁路路堤断面计算域如图3所示。包括＜0＞、＜4＞、＜6＞、＜13－1W4＞、＜13－W4＞、＜13－W2＞ 共6个不规则地层，P1，P2和P3，P4之间分别为2轨道宽度。加载应用程序ptline.lsp从AutoCAD模型提取6个封闭子域的点线信息，运行自动剖分程序即得计算域的网格自动剖分结果如图4至图6所示，粗网格四边形单元数215，节点数249;细网格四边形单元数756，节点数811。优化前后单元质量分布情况见表1。

(2)瀑布沟电站地下厂房2#机组剖面计算域如图7至图8所示:水平方向沿主厂房中心线向上游取150m;向下游取250m;铅直方向由500m高程处取至地表。计算模型中考虑的断层、挤压破碎带主要包括:f13，f29，f19，g(1)，f(7)，f(8)，f(1)，f28。网格剖分结果如图9至图10所示。优化后三角形网格单元数4 744，节点数2 492，其中质量好的单元4 382个，质量较差的单元66个。优化后四边形网格单元数9 678，节点数9 920，其中β≥0.5的单元6 235个，网格整体剖分效果较好。

图3 某铁路路堤断面计算域Fig.3 Computational domain of the railway embankment section

图4 粗网格Fig.4 Coarse grids

图5 细网格Fig.5 Fine grids

图6 优化后的局部四角形网格Fig.6 Optimized local quadrilateral grids

表1 优化前后单元质量分布Table 1 Distribution of quadrilateral element quality before and after the optimization

图7 瀑布沟电站厂房2#机组剖面Fig.7 Computatonal domain of the profile of No.2 unit of Pubugou power station underground plant

图8 四边形网格Fig.8 Quadrilateral grids

图9 局部三角形网格Fig.9 Local triangular grids

图10 局部四边形网格Fig.10 Local quadrilateral grids

5 结语

本文提出的四边形网格自动剖分的改进行波法，无需设置背景网格，直接由AutoCAD模型出发，加强了网格剖分的过程控制，能迅速而有力地生成高质量的三角形网格和四边形网格，实现多介质复杂区域网格的自动剖分，简单易行，适应性强，既可满足有限元计算的需要，又可用于大型有限元计算系统的集成开发，具有很好的应用价值。

［1］冯道雨.边坡稳定分析的自适应有限元法研究［D］.武汉:武汉大学，2004.(FENG Dao-yu.Self-adaptive Finite Element for Slope Stability Analysis［D］.Wuhan:Wuhan University，2004.(in Chinese))

［2］冯道雨，陈尚法，陈胜宏.复杂区域四边形网格生成的一种改进方法［J］.岩土力学，2004，25(6):917－921.(FENG Dao-yu，CHEN Shang-fa，CHEN Sheng-hong.A New Quadrilateral Mesh Generation Method Based on Advancing Front Technique［J］.Rock and Soil Mechanics，2004，25(6):917－921.(in Chinese))

［3］刘春太，杨晓东，陈静波，等，任意平面域自渐变三角形网格的自动划分［J］.计算机辅助设计与图形学学报，1999，11(4):293－295.(LIU Chun-tai，YANG Xiao-dong，CHEN Jing-bo，etal.Automatic Generation of Self-changing Triangle Grid for Arbitrary Planar Domains［J］.Journal of Computer-Aided Design ＆ Computer Graphics，1999，11(4):293－295.(in Chinese))

［4］杨晓东，申长雨，陈静波，等，任意平面区域的变密度四边形网格生成方法［J］.工程力学，2001，18(2):135－139.(YANG Xiao-dong，SHEN Chang-yu，CHEN Jing-bo，etal.Generation of Graded Quadrilateral Meshes for Arbitrary Planar Domains［J］.Engineering Mechanics，2001，18(2):135－139.(in Chinese))

［5］CHEN Jian-jun，ZHENG Yao.A Robust and Quality Guaranteed Pattern Module Scheme for Multi-subdomain Methods in Mesh Generation［J］.Journal of Computer-Aided Design and Computer Graphics，2005，17(10):2286－2292.