有限元网格剖分原理
来源：机械自动化前沿，检索发现文章最早发布于笨笨鸟新浪博客。
一、引言
有限元法是求解复杂工程问题的一种近似数值解法，现已广泛应用到力学、热学、电磁学等各个学科，主要分析工作环境下物体的线性和非线性静动态特性等性能。
有限元法求解问题的基本过程主要包括：分析对象的离散化，有限元求解和计算结果的处理三部分。
曾经有人做过统计：三个阶段所用的时间分别占总时间的40%~50%、5%及50%~55%。也就是说，当利用有限元分析对象时，主要时间是用于对象的离散及结果的处理。如果采用人工方法离散对象和处理计算结果，势必费力、费时且极易出错，尤其当分析模型复杂时，采用人工方法甚至很难进行，这将严重影响高级有限元分析程序的推广和使用。因此，开展自动离散对象及结果的计算机可视化显示的研究是一项重要而紧迫的任务。
可喜的是，随着计算机及计算技术的飞速发展，出现了开发对象的自动离散及有限元分析结果的计算机可视化显示的热潮，使有限元分析的“瓶颈”现象得以逐步解决，对象的离散从手工到半自动到全自动，从简单对象的单维单一网格到复杂对象的多维多种网格单元，从单材料到多种材料，从单纯的离散到自适应离散，从对象的性能校核到自动自适应动态设计、分析，这些重大发展使有限元分析摆脱了仅为性能校核工具的原始阶段，计算结果的计算机可视化显示从简单的应力、位移和温度等场的静动态显示、彩色调色显示，一跃成为对受载对象可能出现缺陷（裂纹等）的位置、形状、大小及其可能波及区域的显示等，这种从抽象数据到计算机形象化显示的飞跃是现在甚至将来计算机集成设计、分析的重要组成部分。
二、有限元分析对网格剖分的要求
有限元网格生成就是将工作环境下的物体离散成简单单元的过程，常用的简单单元包括：一维杆元及集中质量元、二维三角形、四边形元和三维四面体元、五面体元和六面体元。他们的边界形状主要有直线型、曲线型和曲面型。对于边界为曲线（面）型的单元，有限元分析要求各边或面上有若干点，这样，既可保证单元的形状，又可提高求解精度、准确性及加快收敛速度。不同维数的同一物体可以剖分为由多种单元混合而成的网格。网格剖分应满足以下要求：
合法性。一个单元的结点不能落入其他单元内部，在单元边界上的结点均应作为单元的结点，不可丢弃。
相容性。单元必须落在待分区域内部，不能落入外部，且单元并集等于待分区域。
逼近精确性。待分区域的顶点（包括特殊点）必须是单元的结点，待分区域的边界（包括特殊边及面）被单元边界所逼近。
良好的单元形状。单元最佳形状是正多边形或正多面体。
良好的剖分过渡性。单元之间过渡应相对平稳，否则，将影响计算结果的准确性甚至使有限元计算无法计算下去。
网格剖分的自适应性。在几何尖角处、应力温度等变化大处网格应密，其他部位应较稀疏，这样可保证计算解精确可靠。
三、现有有限元网格剖分方法
K. Ho-Le 对网格生成算法进行了系统分类，该分类方法可沿用至今，它们是拓扑分解法、结点连元法、网格模板法、映射法和几何分解法五种。目前，主要是上述方法的混合使用及现代技术的综合应用。
1. 映射法
映射法是一种半自动网格生成方法，根据映射函数的不同，主要可分为超限映射和等参映射。因前一种映射在几何逼近精度上比后一种高，故被广泛采用。
映射法的基本思想是：在简单区域内采用某种映射函数构造简单区域的边界点和内点，并按某种规则连接结点构成网格单元。这种方法可以很方便地生成四边形和六面体单元，若需要，也很容易转换成三角形和四面体单元。
该法的主要缺点：首先必须将待分区域子划分为所要求的简单区域，这是一个十分复杂且很难实现自动化的过程。对复杂域采用手工方法划分甚至不可能，通常各简单区域边界采用等份划分。另外，该法在控制单元形状及网格密度方面是困难的。鉴于简单区域自动划分的困难性，Blacker试图采用知识系统和联合体素方法解决，但在复杂多孔域上仍难以处理，主要是体素数量和形状有限，将待分区域全自动划分为有限的预定体素并集是很难完全实现的。
2. 拓扑分解法
在不考虑网格单元大小和形状情况下，Wordenweber提出了使用三种算子连接多边形各顶点形成粗三角形的二维拓扑分解法，然后细化粗单元至预定规定的网格密度为止。三种算子使用顺序为：opj、opl、op0 。同时，Wordenweber也提出了在三维域使用opj(i=0,1,2,3) 和opp五种算子剖分实体。Woo、Thomasma和Saxena等扩充了该法，并将其有效地应用到多面体实体有限元自动网格生成中。Saxena称该法为EE法，并已与RSD法混合使用构成RSD/EE法。单一的拓扑分解法因只依赖于几何体的拓扑结构使网格剖分不理想，有时甚至很差。
3. 几何分解法
凡在产生结点的同时也确定结点间连接关系的方法均称为几何分解法，常用的有两种：递归法和迭代法。
递归法：Tracy、左建政和Chae等先离散二维物体边界，然后沿离散边界向物体内挖掉一个、两个或三个三角形，重复此操作直到区域挖空为止。Lindhom、Blacker和B.P. Johnston等使用的迭代法不同于前者，首先从物体中挖掉边界层而不是单元，然后三角化边界层。上述为二维迭代法，Chae在此基础上发展了三维迭代几何分解法，主要分两步：采用二维迭代几何分解法生成表面三角形，然后采用三种算子挖切凸体为四面体。在挖切时，突出的特点在于采用新方法生成关键点。关键点的生成分两步考虑：一是考虑新点对周围面和边的影响，二是通过调整比例因子来确定新点位置。Chae也将所提出的算法成功地应用于自适应网格生成中，但由于被剖分物体形状必须是单连通凸域，因此，不能实现全自动网格生成。
迭代法：Bykat采用该法。他首先将物体划分为凸体（手工或自动），随后根据网格密度分布，在每个凸体边界上插入结点，然后将物体中间“最长轴”一分为二，在该轴上插入结点，继续对两部分做递归分割直到最后子域均为三角形为止。商业网格生成软件Triquamesh仍采用该法，只是分割线的选取与Bykat不同。几何分解法的最大优点是在离散物体时考虑网格单元的形状和大小，因此，所生成的网格单元形状和分布均较好。最大缺点是自动化程度低，不利于复杂件网格生成。
4. 网格模板法（RSD法）
Shephard、Perucchio、Saxena、Sapidis和Yerry等，是这种方法成功运用的主要代表。网格模板法生成有限元网格主要分两步（以介绍三维实体为主）：
第一，将待分实体用适当大小的立方体箱（树根）完全包容，按“一化八”原则递归离散，然后对每个八分块按如下方法进行分类：Procedure ModClassCell(Cell,S)=('IN','OUT','NIO') If (八分块中至少有一个顶点为'OUT'且至少有一个顶点为'IN') then 'NIO' Else if (Cell (*S=() then 'OUT' Else if (Cell (*S=Cell) the 'IN' Else 'NIO' End; {procedure} 对于IN的八分块继续递归离散直到预定水平级为止，OUT的八分块不再划分，NIO的八分块进一步子划分，且分类直到预定水平级为止。称终了IN和NIO八分块的并集为RSD模型。
第二，对已经形成的RSD模型，目前已有多种生成网格的处理方法，主要有RSD/GDT法、RSD/EE法和RSD/DDT法。它们主要有以下特点：
RSD/EE法不能处理曲面实体、非流形体和不连通实体。与此相反，RSD/DDT法却能处理有孔的任意曲面实体、非流形体和不连通实体，而且所形成四面体形状质量良好。
RSD/DDT法根据需要以满足条件为准则插入新点，因此所插入的新点数量少，而RSD/GDT法则会插入许多冗余点。
RSD/GDT法使用点/实体分类，使时间复杂性至少大一个数量级，而RSD/DDT法不使用点/实体分类，因此，RSD/DDT法平均时间复杂性为O(N2)，N为实体S的总表面数。RSD/EE法具有不确定的时间复杂性。
RSD/DDT法完全建立网格图素拓扑一一对应，因此拓扑是健全的，与此相反，RSD/GDT法是拓扑不健全的。各种RSD法的优点是网格生成完全自动，网格剖分速度快，非常适用于自适应网格生成。主要缺点是边界单元形状难于完全保证。另外，RSD法对物体的方向特别敏感。
5. 结点连元法
结点连元法是先生成结点，然后连接结点构成单元。最常用的是DT法和AFM法。
DT法的基本原理：任意给定N个平面点Pi(i=1,2,…,N) 构成的点集为S，称满足下列条件的点集Vi 为Voronoi多边形。其中，Vi 满足下列条件：Vi={ X:|X- Pi|(|X- Pj|,X(R2,i(j,j=1,2,…,N }Vi为凸多边形，称{Vi}mi=1为Dirichlet Tesselation 图或对偶的Voronoi 图。连接相邻Voronoi 多边形的内核点可构成三角形Tk，称集 合{Tk} 为Delaunay 三角剖分。DT法的最大优点是遵循“最小角最大”和“空球”准则。因此，在各种二维三角剖分中，只有Delaunay三角剖分才同时满足全局和局部最优。“最小角最大”准则是在不出现奇异性的情况下，Delaunay三角剖分最小角之和均大于任何非Delaunay剖分所形成三角形最小角之和。“空球”准则是Delaunay三角剖分中任意三角形的外接圆（四面体为外接球）内不包括其他结点。实现Delaunay三角剖分有多钟方法，Lee和Schachter操作很有效，但很难实现；而Watson、Cline和Renka、Sloan因操作容易、时间效率较好等优点而被广泛采用。为了进一步提高效率，Sloan研究其算法操作，提出了时间复杂性为O(N)（N为结点总数）的操作方法，从而为快速Delaunay 三角剖分提供了有效途径。虽然DT法既适用于二维域也适用于三维域，但直接的Delaunay 三角剖分只适用于凸域，不适用于非凸域，因此发展了多种非凸域的Delaunay剖分。
AFM法的基本原理：设区域的有向离散外边界集和边界前沿点集已经确定，按某种条件沿区域边界向区域内部扣除三角形（四面体）直到区域为空集。
AFM法的关键技术有两个：一是区域的边界离散和和内点的合理生成，二是扣除三角形条件。目前，扣除三角形的条件有多种。
最短距离条件。选取到该区域边界前沿垂直距离最短的点或到边界前沿端点距离平方和最小的点构成三角形(四面体)。
最大角条件。在平面区域选与有向边界前沿BC边构成角BAC最大的A点）。实体区域选与有向边界前沿三角形ABC构成的四面体ABCD在D点实体角最大的点。这种选取点构成三角形（四面体）方法可靠，但由于对应于同一边界前沿可能存在多个最大角情况，我们称这种现象为奇异现象。
最大形状质量条件。选取与边界前沿所构成的三角形（四面体）形状质量最大的点构成有效剖分。
最小外接圆（球）条件，即空球条件。选取与边界前沿构成的三角形（四面体）中外接圆（球）半径最小的点构成有效剖分。即在所形成的三角形（四面体）中不包含任何其他边界前沿点集。单一使用上述四种四面体扣除条件均会出现奇异情况。使用后两者扣除单元都将可能引起剖分不可靠，如不可剖分及单元相交等。
AFM法最大优点是，不仅在区域内部而且在区域边界所生成的网格单元形状均优良，网格生成全自动，可剖分任意实体。如果将板/壳、实体和梁采用统一的数据结构，则可采用该原理实现不同维数和多种材料等混合工况结构件的网格自动剖分。若配合误差估计，则这种方法在自适应网格再生技术中使用效果甚佳。目前的发展趋势是采用AFM/DT混合法。在平面域已得到了成功地实现，但三维实体区域仍存在多种问题，例如：可能出现剖分不可靠和奇异等现象。
四、自动自适应网格剖分
有限元的自适应性就是在现有网格基础上，根据有限元计算结果估计计算误差、重新剖分网格和再计算的一个闭路循环过程。当误差达到预规定值时，自适应过程结束。因此，有效的误差估计和良好的自适应网格生成是自适应有限元分析两大关键技术。就目前国外研究来看，自动自适应网格生成从大的方面可分为：网格增加技术和网格再生技术。
1. 网格增加技术
该法主要依靠增加自由度总数来提高有限元分析的精度。目前主要采用三种类型方法提高有限元分析精度：h-型、p-型和h-p-型。
h-型，采用有选择地进一步子划分网格单元来细化网格以提高自由度。该法使用特别广泛，RSD模型的网格改进正是利用该法。
p-型，在保持网格划分不变的情况下，通过提高插值函数的阶数获得高的求解精度。
h-p-型，是将h-型和p-型两种结合的一种方法。该法虽然实现不容易，但它却可使收敛速率明显加快。
实践表明，在获得同一精度时，上述三种类型收敛速率是按h-型、p-型、h-p-型顺序增加的。
2. 网格再生技术
根据现有网格并配合误差估计确定新的结点密度分布，然后重新划分网格，再计算并重复上述过程直到求解精度达到预定目标为止。目前，网格再生技术在平面区域已得到了较好地实现，现概括如下：
确定结点密度的最大和最小值。
按等比数列计算等值线分段总数N和每一条等值线结点密度值。
用N 1个结点密度将区域划分为相应的等值线。
光顺每条等值线。消除等值线与区域边界及等值线之间的尖角，并沿等值线生成结点。
形成封闭环。从内部等值线和区域边界中获得边界链，并按正确顺序连接边界链和内部等值线链构成封闭环。
形成子区域。子区域的结点密度等于区域边界结点密度的平均值。
网格生成。采用二维平面网格生成方法构成三角形。
从理论上讲，该原理可扩充到三维实体域，但由于三维实体域难以完全自动用等结点密度曲面来分割任意实体，因此在三维域的扩充至今仍未实现。实践表明，网格再生技术比网格增加技术具有更大的优点，主要表现在前者收敛速度快、网格单元形状稳定。
另外，还有两种自适应网格生成类型，它们是r-型和h-r-型。
r-型，是保持网格划分和插值函数阶数不变情况下，通过调整结点位置以改善求解精度。该法收敛速度低，因此，目前很少直接使用这种类型。
h-r-型是上述r-型和h-型两种方法的综合。
事实上，在h-型网格改进中所使用的光顺技术就是一种r-型，所以可以把h-r-型看成是h-型。
