六面体单元生成方法及相关技术

1、3.1 引言有限元法是求解偏微分方程描述的连续体问题的一种近似工程方法。为克服实际连续体问题难于处理，它将分析区域离散化，将偏微分方程转化为线性方，采用数值计算方法求出连续体问题的近似解。用有限元法进行工程分析的主要过程包括三个阶段：（1）有限元模型的建立和数据准备；（2）用分析计算；（3）分析结果的判断和评定。迅速而合理地划分有限元网格是完成有限元分析的前提和保证。有限元网格生成技术发展到现在，已经出现了大量的、不同的实现方法 22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39。网格生成方法一般可以分为两类：结构化网格生成方法和非结构网

2、格生成方法。所谓结构化网格，严格的讲，是指网格节点具有相同的邻接单元的网格，即网格节点具有相同的拓扑结构。结构化网格生成通常使用复合的迭代平滑技术，用边界或物理区域来排列单元，划分区域的边界不能太复杂，以便能将区域分解为拓扑结构的块。所谓非结构化网格，放松了对节点的连接要求，允许一个节点属于任意数量的单元。虽然四边形和六面体单元也可以是非结构网格，但通常非结构网格是指三角形和四面体网格。目前二维区域网格自动划分发展较为成熟，许多商业都提供了成工具。而三维区域的网格自动划分，除个别商业提供了四面体网格的自动划分外，还不成熟，特别是形状较规范的六面体网格的自动生成。一般来说，由程序自动划分的网格总

3、有一部分的单元的形状不是很好，需要对网格的质量进行优化。节点序号的的标注方法直接影响刚度矩阵容量的大小从而影响计算的效率，因而需要对网格的节点重新，以得到具有较小带宽的网格系统。本章主要了三维六面体网格生成的方法及其相关技术。3.2 结构化网格生成结构化网格生成方法主要分为两类：代数法和 PDE（Partial Differential Equation，偏微分方法和保角法；PDE 方法包括椭圆系统方程、Poisson程）方法。代数法包括超限插值法、等参系统方程、抛物线系统方程和双曲线系统方程方法。应用最多的是等参法。结构化网格的优点是易于实现，在每个子区域内网格可以得到很好的控制，生成规则的

4、结构化网格，而且能生成曲面网格。缺点是对不规则的形体有时生成的网格质量很差，需要事先根据所要产生的网格类型将目标区域分割为一系列可的子区域，这一步通常需要手工完成，因而自动化程度低，不适合网格全自动生成。提供结构化网格生成的有553DMAGGS（NASA，Langley & Lockhead，四边形和六面体单元）、CAGI（ERC，Miss. Se，四边形和六面体单元）、CSCMDO（CSC/NASA LaRC，四边形单元）、EMC2（INRIA，四边形单元）、GENIE+（ERC，Miss. Se，四边形和六19面体单元）、ICEM CFD（ICEM CFD Engineering，四边形和

5、六面体单元）、QulkGrid（Edge，四边形单元）、TrueGrid（XYZ Scientific Applications，Inc.，四边形和六面体单元）、VGM（NASA，Langley & Lockhead，四边形和六面体单元）等。3.2.1 代数法网格生成代数法网格生成技术40是通过插值函数将理想的直角坐标系表示的计算区域变换为实际物理区域来实现的（图3-1）。它的具体做法是将计算区域的直角坐标用均匀的间隔划分为计算区域的网格，通过变换，计算区域上均匀分布的坐标被为物理区域上的坐标（图3-2）。网格点数目的多寡并不影响变换的特性。从计算区域到物理区域的变换是通过一个向量函数完成的物

6、理区域计算区域zyx图3-1计算区域和物理区域Fig. 3-1 Compuionaland physical计算网格物理网格x = x(, , )y = y(, , )z = z(, , )zyx图3-2计算区域网格和物理区域网格Fig. 3-2 Gridshe compu ionalandhe physical20 x , , F , , y , , 0 1, 0 1, 0 1(3-1) z , , 多变量的坐标变换式(3-1)使用无限插值(Transfiniteolation)来完成。在三维情况下，使用混合函数和与之相联系的参数（特定点的位置及偏导数）来显式地决定式(3-1)，然后通过对每

7、个单变量的循环完成无限插值。一般来说，这些混合函数和参数都选取为区域边界处的函数和参数。文献40给出了几种无限插值函数及其参数的选取。代数法无限插值网格生成法的特点是计算简单，可以采用中间变量的方法方便地控制网格的密度，对边界简单的区域，可以生成质量较高的网格，但缺点是不适应复杂边界的划分，边界不规则时生成的网格的质量很差，并可能产生奇异性。可以通过将划分区域分解为子区域的方法，在子区域上应用，可以在一定程度上克服这些缺点，但不易实现自动划分。3.2.2 偏微分方程法类似于代数法，偏微分方程法也是通过将曲的边界为网格坐标线，只不过这里使用的函数满足某一类微分方程。以二维椭圆偏微分方程为例41，

8、满足最大条件的椭圆型微分方可以为, Q , (3-2)xxyy右端的强迫函数为 MN 12P , a ecm2m2iim i m1i 1mm(3-3) d MN12Q , a m m1 e cm m2b2mmeiii i i 1mm在变换空间中，这些函数写成为 Qx (3-4)Py Qy y 2y y J2 yy ，而这里 J 是变换的 Jacobi 矩阵行列式的值，J = x y , x x y y , x 2 y 222 和代数法相比，偏微分方程法的计算较为复杂，需要求解偏微分方程，右端强迫函数的选取也不方便，可以通过强迫函数的选取得到希望的网格。213.2.3 超单元法等参超单元网格生成

9、法是由 Zienkiewicz 和31,42，其基本Phillips 最早是将划分区域分为自然坐标系19成更简单的子块（超单元），子块被5中的直角形（正方形），根据每个方向给定的分级权系数，将自然坐标系中的直角形（正方形）离散并反向变换回超单元。在三维情况使用超单元生成法时，要划分的实体用 20 节点六面体单元分成数个子区域（图3-3）。在超单元内一点的坐标 x，y 和 z 和自然坐标之间的关系由下式确定163202020 x Ni xi , y Ni yi , z Ni zi1(3-5)92i 1i1i 1式中，Ni(, , )是与按曲线坐标系， 和定义的与每个节点相联系的形函数， 和的取值

10、范围为1 到 1。 20 节点六面体单元形函数在各角点定义为图3-3 20 节点六面体单元Fig. 3-3 20-node hexahedral element, i 1, 2, , 8(3-6)在各棱边的中间节点定义为N ( , , ) 1 1 2 1 1 , i 9,11, 17, 19iii41 1N ( , , ) 1 1 2, i 10,12, 18, 20(3-7)iii4N ( , , ) 1 1 1 1 2 , i 13,14, 15,16iii4如节点的坐标已知，则与特定坐标， 和相对应的点的坐标可由式(3-5)计算出。沿每个局部自然坐标，给定划分的数目和划分的间距，这些子区

11、域（超单元）被划分为更小的单元。以为例，划分时自然坐标按下式增加 i (3-8)0 = 1。法计算简单，通过对自然坐标的不同划分，可以得到密度变化的网格，但超单元这里W超单元描述的边界和划分区域的实际边界存在误差，不能适应复杂的边界。223.3 非结构网格四边形/六面体网格生成法由于算法的特点，自动非结构网格生成算法主要是三角形或四面体网格划分，有许多文献讨论了三角形或四面体网格划分方法和相关的技术43,44,45,46,47,48,49,50,51,52,53,54，许多也提供了三角形或四面体网格划分功能55。非结构网格生成的特点是能适应复杂的边界，需要的人工干预少，计算费时。Steven.

12、 J. Owen55总结了当前学术界和工业界使用的网格划分方法和可以利用的。非结构六面体网格生成方法主要有法、转换法、铺层算法和须段编织算法等。不象三角形网格划分算法向四面体网格划分算法扩展那样直接，从二维四边形网格划分算法向三维六面体网格算法的扩展并不是直接的。3.3.1法当划分区域的几何特性合适时，四边形或六面体网格生成法可以得到相当好的网格56。虽然法被认为是结构化网格生成法，但许多非结构网格划分都提供这种方法。当应用映射法时，要求划分区域的对边具有相同的子划分。在三维情况下，将要划分区域的拓扑结构的立方体的每个相对的面必须有相同的表面网格。这种情况在几何形状复杂是不可能达到，或者需要人

13、工干预将几何区域划分成可的子区域并分配适当的间隔。等57,58使用的方法是将划分区域分成子区域，每个子区域对应于不同的模块，每个模块中的网格已预先划分，模块中的网格为子区域中的实际网格。这里模块中的网格可以是结构化的，也可以是非结构化的，但要求相邻子区域的交接面上的表面网格具有相同的拓扑结构。该方法的特点是可以针对子区域的几何特性和物理特性选择相应的模块，生成的网格质量有 ANSYS（ANSYS Inc.）、CUBIT（Sandia高，缺点是通用性差。提供National Laboratories）、goe Inc.）等。法生成网格的Wrap（Raindrop GeoInc.）、Hyperme

14、sh（Altair Computing扫描线目标源图3-4 扫描法生成六面体单元Fig. 3-4 Hex elements generated by s）网格生成法是另一类六面体网格生成方法59,60,61，有时被认为是二扫描法（S维半网格生成方法。该方法是将四边形网格沿空间一条曲线扫描，给定一定的间隔，使用该四边形网格的拓扑结构生成给定层的六面体网格。这种方法归纳为给定源表面和目标表面划分某些类型的体积区域的技术。假定源表面和目标表面具有相同的拓扑结构并被可划分表面相连接，23源表面上的四边形网格扫过体积区域生成六面体网格（图3-4）。提供扫描法生成网格的有ANSYS（ANSYS Inc.）

15、、CUBIT（Sandia National Laboratories）、Hypermesh（Altair Computing Inc.）等。3.3.2 非结构四边形网格划分非结构四边形网格划分算法一般可分为两类：直接方法和间接方法。间接方法是先将划分区域划分为三角形，然后采用各种算法将三角形转换为四边形。直接法是在划分区域上直接生成四边形，不需要先生成三角形。图3-5 三角形分为四边形生成的网格Fig. 3-5 Quad mesh generated by splitting each triangleo three quads图3-6 三角形组合为四边形生成的网格Fig. 3-6 Quad

16、-dominant mesh generated by combining triangles1三角形分解/合并成四边形（间接法）最简单的间接四边形网格生成是将每个三角形划分为三个四边形，如图3-5所示，这种方法保证可以生成全四边形的网格，但网格中大量存在的不规则节点导致低的网格质量。另一种算法是组合相邻的两个三角形形成一个四边形，如图3-6所示，这种方法可以增加网格的质量，但同时会o62给出的算法留下大量的三角形单元。在三角形组合次序上采取某些措施，可以改进该方是通过分析可能的合并，按照合并后四边形的优劣结果决定合并的次序，得到一个包含最少三角等63对Lo 的算法进行了改进。Johnson

17、等64提出采用增加形单元的四边形占优的网格。闵对局部单元的划分和交换的办法来提高四边形单元的数量和质量。Lee 等65通过引入局部三角形的分解扩展了 Lo 的算法62，并使用了波前的概念，如果区域边界上初始的单元边数是偶数，这种方法可以组全四边形网格。Owen 等66Q-Morph 算法用于从三维的三角形网格中生24成四边形网格，该算法采用决定三角形转换的序列，使用三角形网格中已存在的边或额外的节点或对局部三角形进行变换来形成四边形，任意数量的三角形可以合并生成一个四边形。Q-Morph 算法可以生成沿区域边界排列的很好的四边形，但在区域还是存在数目有限的不规则节点。Q-Morph 算法在有限

18、元ANSYS 上得到了实现。由于三角形分解合并为四边形的所有操作都限于局部，故间接算法的速度非常快，不需要象直接算法那样进行相交检查。其缺点是网格中会留下许多不规则的节点。对某些应用而言，即使留下的不规则节点是少数，也不能保证沿区域边界的单元能以希望的形式排列。可以通过拓扑结构的调整去掉这些不规则的节点从而提高网格的质量。2栅格法Baehmann 等67采用的四叉树分解方法是将目标区域用尽量小的正方形圈定，然后把这个正方形分解为四个子区域，对每个子区域测试其是否完全落在目标区域之外，或是否满足密度控制要求，若满足则停止对区域的细分，否则细分下去，递归执行直到到达预定的要求，留下合适的四边形单元

19、，调节节点位置保持与边界的一致。栅格法（Grid Method）68,69是将一包含目标区域的栅格放置在目标区域上面，除去落在区域外的栅格单元，并将与区域边界相交的栅格单元进行调整或裁剪，然后在栅格边界单元和区域边界之间生成四边形单元（图3-7）。区域的网格最终的网格图3-7 基于栅格的四边形网格生成Fig. 3-7 Quadrilateral mesh generated by rigid-based methodTalbert 等70使用的分解方法是将区域迭代划分为简单的多边形子区域，这些子多边形的形状与有限数量的模板相匹配，然后在这些多边形中四边形单元。Tam 和 Armstrong71

20、利用区域的中间轴（Medial Axis）将区域分解为简单凸多边形区域，然后使用一套模板将四边形网格区域。平面对象的中间轴是对象中可放入最大圆的圆心轨迹，即沿对象滚动的最大直径的圆的圆心轨迹（图3-8）。中间轴由数个直线或曲线段组成，采用一定的规则用这些线段将区域分解为凸的子区域，再采用界拓扑相容的网格。的方法在子区域上生成与子区域边有限四叉树分解法和基于网栅的分解法在区域可以生成质量很高的单元，但在区域边界上要满足边界一致性要求，边界单元的质量较差，不适合塑性成形变形大、网格容易发生畸变的应用领域。分解为子区域的划分方法，通过控制所划分的子区域的形状和增加可模板的数量，25可以提高子区域和模

21、板的相似性，从而使网格质量得到提高，然而对于一般几何形状的区域，难于做到通用性。中间轴图3-8 二维对象的中间轴Fig. 3-8 Medial axis of a 2D object图3-9 铺路算法示意Fig. 3-9 Demo of paving algorithm3四边形网格生成的Zhu 等72最早提出了使用的四边形网格的生成算法。该算法是先在边界上生成初始节点，然后向区域投射边段，先使用传统的三角形网格产生两个三角形单元，再将两个三角形组合为一个四边形。Blacker 和 Stephenson铺路算法73（Paving Algorithm）是重复地在区域边界放一等74将铺路层或者说铺一

22、排单元，算法从区域边界开始向形成四边形单元（图3-9）。Cass算法推广到了三维表面的划分，White 和 Kinney 重新设计了铺路算法75，单元不是生成一排，而是一个接一个的生成。方法。CUBIT（Sandia National Laboratories）等提供了铺路算法的网格划分铺路算法能生成边界上具有较高质量的网格，是一种较好的非结构网格生成方法。3.3.3 非结构六面体网格生成类似于非结构四边形网格生成方法，非结构六面体网格生成也可以分为直接法和间接法。图3-10 四面体分解为六面体ition of a tetrahedrono four hexahedral elementsFi

23、g. 3-10261四面体分解为六面体76（间接法）如果体积区域可以划分为四面体网格，则每个四面体单元可以划分为四个六面体单元，如图77采用 103-10所示。一般情况下这种方法生成的单元质量较差，在有限元分析时很少采用。节点的四面体单元分解为四个六面体单元的方法，可以提高边界处单元与划分区域边界的一致性，但对单元质量并没有实质性的提高。还没有文献对四面体单元组。2基于栅格的方法六面体单元的方法进行基于栅格的六面体网格生成方法包括（有限）八叉树分解法78和正则网格覆盖法79。这里先采用八叉树分解或网格覆盖的方法，在划分区域生成合适的三维六面体单元网格，然后在边界处加入或调整六面体网格，将规则网

24、格与边界之间的空白区域充满，并满足边界一致性要求。基于栅格的方法在所划分的体积区域的边界处生成的单元质量不高，通常单元沿边界的排列也不好，同时生成的网格取决于成的单元除边界处外大小相同，具有明显单元尺寸过渡的网格。于栅格的六面体网格生成方法。六面体网格的取向，即与坐标轴的方向有关。正则网格法生等80和Schneider81对八叉树算法进行了修改，可以生成Men （MARCysis Research Corpoion）等提供了基3中间面法（Medial Surface Method）中间面法82,83,84是划分四边形网格中间轴法在三维的直接扩展。这里划分区域被分解为由中间面分开的子区域的集合，

25、中间面可以认为是在区域滚过的最大球的球心点组成的面。由中间面分开所形成的子区域形成可划分网格的区域，这些区域对应于一系列的具有所期望的拓扑结构的模板，整个区域用后的六面体网格充满，使用线性编程来保证从一个区域到另一个区域单元划分之间的匹配。这种方法适合于某些具有特定几何形状的区域的划分，但对一般几何形状的区域可靠性差，生成中间面及规定中间面定义的所有区域还是一个难于解决。中间面划分算法在 CADfix（FEGS.）、TurboMesh（SolidPo）等中得到实现。4铺层算法（Plastering Algorithm）铺层算法85是铺路算法在三的拓展，这种算法是先将单元放在边界处，然后逐层向划

26、分区域中心推进（图3-11），定义一套探索程序决定单元形成的次序。和其它推进类似，所有的四边形当前的波前，每个的四边形向体积区域投射形成六面体。此外，算法必须判断与表面的相交、何时和如何连接原来存在的节点来缝合留下的空间。随着算法的继续，在区域可能留下复杂的空隙，在某些情况下这些区域不能用全六面体单元填充。为了方便地向内放入六面体，铺层算法有时必须修改已生成的单元。铺层算法在大规模问题上应用的可靠性还没有被证明。虽然可以成功地在边界上放上数层六面体单图3-11 铺层法生成单元的过程Fig. 3-11 Plastering Pros forming elements at the boundar

27、y元，但相交和封闭算法的可靠性不高。等提供了铺层算法划分网格的方法。CUBIT275须段编织算法Tautges 等86须段编织算法（Whisker Weaving Algorithm）是基于空间缠连集（SpatialTwist Continuum，STC）的一种拓扑划分方法。STC 是建立在全六面体网格对等体上的一个整体构造集，在 STC 中，用相交面的排列来对等地表示网格的实体，这些面在各个方向上平分六面体单元并形成相交的曲线和顶点，STC 封装了全六面体网格的约束信息。平分每个单元的面有三个，这些面称为须段面（Whisker Sheet），可以认为它是一个一般的的表面，这个表面在网格中延续

28、，在相对的两个面的中间平分每个六面体。图3-12是由两个单元组成的网格，它有 4 个须段面，编号 1 至 4，这些面的排列就称为 STC。一个须段面代表一层相邻的单元之间具有共享面的单元。两个须段面相交形成的曲线称为须段弦（Whisker Chord），图中有 5 个须段弦，用 A 至 E 标识，须段弦表示网格中相邻单元之间具有共享面的一列单元，如须段线 E 定义了这个网格中两个单元组成的单元列。每个六面体单元对应于一个对等的称为 STC 顶点（STC Vertex）的实体，它是三个须段面相交形成的交点（同样是三条须段线的交点）。须段面与几何表面的交线称为环（Loop），它是一条封闭曲线，表示

29、一圈表面四边形，每个相邻四边形具有共享的边。BD12aAbC34图3-12 两个单元的 STCFig. 3-12 The STC for a solid comed of two hexahedra须段面是空间中和相交的二维表面，它们之间的相交的几何位置是隐含的，在须段编织中，须段面用平面面图（Sheet Diagram）表示（图3-13）。当须段面变平时，它的环形成一个多边形，它的须段弦是多边形内的线段。须段面上的须段弦在面图上用一个边列表示，为与须段弦相区别称为面弦（Sheetchord）。须段编织算法从定义一个封闭三维体积实体的四边形表面网格开始的，因此必须先生成划分区域表面的四边形网格

30、。须段编织算法分为四步构造初始的环、须段面、须段弦和面弦。从给定的表面网格取一未遍历过的边，从这条边开始，向四边形的对边遍历，直到从另一边回到开始边就形成了一个环。当所有的四边形都被两个环穿过时所有的环都形成了。环生成后，构造每个须段面的须段弦和对应的面弦。找出三条两两相邻须段弦。28E在须段编织算法中，三条须段弦相交形成 STC 顶点组成六面体，在面图上表现为三张须段面上三对相交的面弦。使三条须段弦相交形成一个新的须段六面体。检查连接的有效性，连接须段弦。重复步骤（2）（4）直到完成编织。完成的标志是没有悬挂的须段留下。完成的编织表示全六面体网格的整体的连接性，而须段编织并没有产生六面体单元

31、几何位置的任何信息，环须段弦也没有决定体积区域实际单元或节点的位置，这些位置信息可以通过使用初始构造算法（Primal Construction Algorithm）得到。须段编织算法可以生成质量较高的网格，但对范围广泛的问题算法的稳定性（Robustness）和可靠性还没有得到证实。6六面体单元占优网格图3-13 须段面图Fig. 3-13 A whisker sheet diagram由于大多数六面体网格生成方法的稳定性较差，有些研究者建议使用六面体四面体混合网格。在许多情况下，六面体占优的网格划分的方法是令人满意的。Owen 给出的一个简单方法是手动划分区域成可划分的子区域和几何上较复杂

32、的子区域，在复杂的子区域上定义四面体单元，在四面体单元和六面体单元的的交界处使用过渡单元。在 ANSYS 网格生成这种算法划分网格。中可以选择Tuchinsky 等87,88使用的方法是组合铺层算法和四面体网格划分算法，使用铺层算法尽可能地将六面体单元推进入划分体积区域，在区域内留下的空隙用四面体单元填充，也可选择五面体单元作为六面体单元和四面体单元的过渡。 CUBIT 可以选择这种方法生成六面体占优的网格。3.4 表面网格生成当前许多网格生成与任意表面单元的生成有关。典型地，这些表面是商业 CAD包生成的 NURBS 曲面。表面单元可以直接作为壳单元使用，也可作为体积网格的输入数据。将二维网

33、格生成算法推广到三维表面网格的划分需要作某些修改。表面网格生成算法可分为参数空间法和直接划分法两类。3.4.1 参数空间法所有的 NURBS 曲面都是基于二维参数空间(u, v)表示的，而二中的网格划分通常效率较高。参数空间算法（Parametric Space Algorithm）是先在划分表面的二维参数空间(u, v)中生成网格，然后将参数空间的坐标(u, v)回实际坐标(x, y, z)。这种方法的缺点是在参数空间中生成的单元回划分表面时，并不总是形成三中的好形状的单元。可以采取两种措施来解决这个问题：修改或重新参数化曲面的参数表示，以得到从参数空间到实际空间的合适的；修改网格生成算法，

34、使在二维区域划分的扁的或不等向的单元能元。为三维好形状、等向的单29为了得到一个好的参数化表示，要求在定义域内表面偏导数和(u, v)变化不是太大。文献 89中有精确弧长的重新参数化定义，但花费过大。通过有选择计算区域内的表面偏导数，调节局部的(u, v)的值保持(u, v)的大小大致为常数，可以定义近似的弧长参数化或称翘曲参数空间（Warped Parametric Space）。在许多情况下，翘曲参数空间能生成合理的网格，但还有许多问题是重新参数化不能适当解决的，由于这种原因，多数有关表面网格划分的文献集中了第二种方法：在二中形成不等向的网格，然后成三的等向网格。实践中使用的一个通用的方法

35、是利用表面偏导数、u 和v 这些可以方便地从 NURBS 曲面中计算出的量。Gee 和 Borouchaki55使用了一个由表面第一基本形式导出的表面度量，该度量是一个 22 矩阵，用于变换参数空间的向量和距离。Cuillire90使用节点密度优先的波前算法完成了参数表面的三角化网格划分，其波前是在划分表面的边界上形成，然后反回参数空间，在参数空间中使用通常的波前算法将网格划分，划分时同样采用了表面度量来决定新点的位置。Lee 和 Joun 给出了一种基于的四边形表面网格生成方法91，这种算法是将复杂的、不利于直接的表面裁剪成三角形，三角形边界上的一个点为参数区域中的一条线段，为避免出现此种不

36、利情形，将区域划分四边形网格，然后反MARC/Men、ANSYS 等参数区域再次，成为三角形的再区域，在再回初始区域的实际网格。中提供了参数空间法生成网格的算法。3.4.2 直接划分法直接法三维表面网格生成方法是直接在划分的几何区域上形成网格而不管几何区域的参数表示。在没有参数化表示可用或表面参数化非常差时，可以使用直接法。Lau 和 Lo92使用波前算法在任意空间曲面生成了三角形网格。他们的方法是给定一个波前的线段，检查与边界上或区域内的的一个节点生成新三角形的可能性，根据曲面曲率、单元质量和给定的单元密度分布函数来优化生成的单元，为保持新生成的节点位于曲面之上，需要将生成的节点投射到曲面上

37、。Cass 等74通过引入节点投射、决定波前网格 、支持周期性表面和单元的平滑等技术，将铺路算法扩展到了三维曲面的四边形网格生成。在CUBIT 中可以利用该方法生成表面四边形网格。3.5 网格后处理一般来说，无论用何种方法生成网格，总会有一些形状不好的单元产生，如内角太大的单元、长条单元、薄单元等。这些单元的存在本身会引起有限元解答的精度降低，甚至会导致解的不收敛。在金属塑性成形有限元模拟过程中，这些形状不好的单元会很快发生，导致频繁的网格再划分。为了减少模拟时间，提高模拟精度，必须采用某种或几种方法，将这些性状不好的单元改性以改进网格的整体质量。对六面体单元，可用于判断单元质量的数据有单元的

38、长宽比，即单元的最大与最小边长的比值。每个节点球面度的误差，球面度的定义参见 6.4 节。可以取每个节点的球面度与标准值差值的百分比绝对值作为判断标准。体积测量。六面体每个节点连接三条边，以三条边为起始于该节点三个向量 r1，r2 和 r3，30对每个节点计算 V = (r1r2)r3，采用 C = |V Ve| / Ve 作为体积测量的判断标准，Ve 为单元的体积。4单元面的翘曲度。六面体单元每个面由四个节点组成，有四条边，每三个节点确定一个平面，平面与另一个节点的距离为 d，单元面的翘曲度定义为最大的 d 与最短边长的比值。5单元的度。六面体单元由三对相对的面，单元的度定义为相对面夹角的锐

39、角。改进网格质量的方法主要有两类：节点平滑和网格整理。节点平滑是在不改变网格拓扑结构的条件下，通过调整节点的位置来提高网格质量。网格整理通常是指通过合并单元、划分单元、单元节点的重新分派等改变网格拓扑结构的方法改善网格质量的方法，由于网格结构的特点，六面体网格很少采用网格整理的方法来提高网格的质量。节点平滑处理技术大多数节点平滑程序是采用迭代过程将节点重定位来改善局部单元质量的。有多种平滑技术已被提出，通常这些方法可分为平均方法、基于优化的方法、基于物理的方法、中间节点定位等四种方法。平均方法在所有平滑算法中，Laplace 平滑算法是应用最广的一种，这种算法是把网格的一个节点放在与该节点用边

40、相连的节点的平均位置上。如图3-14a 所示，每个节点的新位置由与该节点相连节点的平均坐标决定1nP P(3-9)ini1这里 P 是要平滑的节点的位置向量，Pi 直接连接到节点 P 的邻近的节点的位置向量，n 为连接到该节点的节点数。Laplace 平滑是共享边的平滑，而等参平滑(Isoparametric Smoothing)是共享面的平滑，如图 3-14b，使用所有与节点 PS 相连的面，直接与 PS 相连的节点坐标取为正，而在面上对角的节点坐标取为减。图中括号内的加减号数目表示同一节点的加减次数。经过小量的修改，这些方法可应用于任何形状的单元。绝大多数这样的平滑程序需要在网格中所有节点

41、进行多次迭代，直到每个节点移动的距离小于给定的偏差。虽然这些方法存在某些问题，但程序实现简单，因而得到了广泛应用。类似于 Laplace 节点平滑方法，基于周围环绕。Canann 等93对常用节点优节点和单元几何特性的化方法进行了总结。平均，产生了多种迭代节点重单独使用平均方式平滑节点有时会降低局部单元的质量，为避免降低单元质量的现象出现，平均方法通常也要对节点的移动施加某种形式的约束，对节点移动前后局部单元的质量进行比较，只有单元质量有所改善时才移动节点的位置。这些方法通常称为约束Laplace 平滑（ConstrainedLaplan Smoothing）。31P3()P4()P6P2()

42、P3P4PP5()SP1()P13()PP10()P6()P ()P9()8PP1P52P7()P ()P12()11（b）（a）图3-14 Laplace 平滑和等参平滑Fig. 3-14 Laplan smoothing and isoparametric smoothing基于优化的方法基于优化的平滑技术测量环绕一个节点单元的质量，通过计算单元质量关于节点位置的局部梯度来优化单元的质量，节点沿梯度增加的方向移动，直到达到优化的目标94,77。单元质量好坏由单元长度方向的相对关系和单元各个面（边）之间的相对关系来衡量，反映在单元变形上与拉压变形和剪切变形有关，而单元的 Jacobi 矩阵中

43、的各项反映了这些变形，故常用 Jocobi 矩阵分量来衡量单元的质量。有关的算法如下将单元的 Jacobi 矩阵正则化，去掉单元尺寸的影响J ijJ (3-10)ij1Jn这里 n 是问题的维数。2引入能描述单元畸变和应变之间模拟关系的 Green 应变，它的一项可表示为 Jacobi 矩阵的函数Cij Jki J kj(3-11)3单元畸变的度量D C C 1 C2(3-12)ij ijkkn它是 Jacobi 矩阵的函数，也即是节点坐标的函数。对整个区域得到一个畸变的度量值，以此作为目标函数进行优化。目标函数变量较多，可采用共轭梯度法(Conjugate Gradient Method)等

44、进行优化，得到优化的节点坐标值。除此之外，还有多种关于单元质量度量的分析表达可用于单元质量的优化，这些量包括单元32每条边连接的两个面的夹角、每个节点处的球面度（对二维单元为平面角）、单元的扁率（ElementAspect Ratio ）、单元面积（周长）与体积（面积）的比值等。所有这些分析量可以向上面所述那样单独使用，但由于单元和网格的质量不是任意一个分析量可以完全描述的，故一般的做法是将几个分析量采用组合的方法，组成一个目标函数，对这个目标函数进行优化。一般来说目标函数的自变量是网格节点的坐标，变量较多，可采用适合多变量优化的方法进行优化，得到优化的节点坐标。基于优化的方法可以得到质量优良

45、的网格，但由于变量多、约束条件多，导致计算量过大且不易收敛。有的研究者95建议将 Laplace 节点平滑技术和基于优化的方法组合使用来优化网格，通常大部分时间采用 Laplace 平滑，仅3.5.1.3 基于物理的方法部单元形状度量低于一定的限值才转到优化的方法。这种改进网格质量的方法是通过模拟物理上任一节点和它周围的节点之间存在着与距离有关的或排斥力的方法来重定位节点的。网格中或排斥力，重定位的结果是该点处于使它所受的力处于平衡的位置。Lohner 等96用一个相互作用的弹簧系统来模拟相邻节点之间的关系，Shimada 等97和sen 等98把节点当作气泡的中心，重新定位这些气泡的位置使它

46、们处于平衡，改变这些粒子（气泡）之间力的大小和方向可以改变单元等向性和尺寸的大小。3.5.1.4 中间节点定位法多数平滑方法都是针对单元的角点进行的。对存在中间节点的单元，Salem 等99提出了重定位中间节点来提高单元质量的方法，这种方法计算环绕中间节点的一个称为中间节点容许空间（Mid-node Admissible Space，MAS）的区域，在这个区域内可以安全地移动中间节点，保持或改进单元的质量。3.5.2 单元细分单元细分定义为在网格上完成的任何能有效减小单元尺寸的操作。减小单元尺寸可能是为了能得到更详细的答数据，也可能仅仅是为了改善局部单元的质量。有些细分方法本身可以认为是网格生

47、成算法，从一个粗分的网格开始，用一个网格细分过程细分直到达到希望的节点密度。细分算法通常被作为自适应分析的一部分，由先前的解答结果提供用于网格细分的准则。图3-15 局部四边形网格细分Fig. 3-15 Exle of local quad refinement由于四边形六面体网格本身的特性，通常不能应用点和边等分方法来细分网格。四边33形和六面体网格细分的主要方法是基于一套预定义的模板来分解单元。文献100描述了一系列的单元分解模板和算法，图3-15显示了一个四边形局部细分的例子。为了维持网格的一致性，某些四边形和六面体细分方法场需要引入三角形单元或四面体、五面体等单元。3.6 带宽优化有限

48、元分析等许多科学和工程的数值计算问题都需要求解一个线性代数方Kx f(3-13)这里已知向量 f 和方阵 K，求解得到未知的向量 x。在有限元分析中，K 是一个稀疏矩阵（矩阵中的零元素远多于非零元素），矩阵 K 中非零元素的位置取决于未知量排列的次序，在有限元分析中未知量排列的次序取决于节点的次序。有限元分析中导出的刚度矩阵等的阶数很高，一般是对称矩阵，完全需要很大的空间，一般采用带宽技术来刚度矩阵。带宽技术的量取决于矩阵中主对角元素与该行非零量的最大距离（按行为 1 时，对节点进行重新），而带宽与单元节点的最大差值直接相关。当每个节点的度等价于对方程(3-13)进行重新排序；当每个节点的度为

49、 m 时，则每个节点对应于 m 个方程，对节点重新小带宽的结果是一样的。只是一个比对度重标号小 O(m)，其减减小带宽的算法主要有三种：直接减小算法（Direct Reduction Algorithms）、图理论算法（Graph Theory Algorithms）和混合算法（Hybrid Algorithms）。直接减小算法是通过直接交换矩阵中的行和列来减小带宽的，这种方法可以得到好的结果，但计算花费的时间很大，这种算法的例子可参考文献101。图理论算法和图理论算法与直接减小算法的混合方法计算效率高、结果好，以下主要对后两种方法进行说明102,103,104,105,106,107,108

50、,109。3.6.1 有关的定义最合适的方程序号重标算法取决于求解方程所使用的求解方法。有限元分析中根据刚度矩阵等所使用的方式，可分为二维等带宽方法、一维变带宽方法和。为评估算法的有效性，需要定义和矩阵有关的参数。设矩阵 K 为 n 阶对称方阵，定义第 i 行的行带宽i(K)为该行对角线上的元素到该行最后一个不为零的元素的列数（包括对角元素） K maxj i 1, j i | A 0(3-14)iij这里假定 Aii 0，0 i n。矩阵 K 的带宽(K)（Bandwidth）和廓 p(K)（Profile）定义为 K max(3-15)npK i K i1(3-16)34在矩阵 K 第 i

51、 行中，如果第 j 列满足 j i，同时在任一行 k i 中存在 Kkj 0，则说 i 行中的第j 列是活动的。用 ci(K)表示行 i 中的活动列数，则矩阵 K 的最大波前定义为wmaxK maxci K , 0 i n方阵 K 是对称矩阵，故有(3-17)nnpK i K ci K (3-18)i 1i 1平均波前 wavg(K)和均波前 wrms(K)定义为pK 1nwK c K i 1(3-19)avginnn1nwK c K i2(3-20)rmsi 1从以上的定义中，有wavg wrms wmax n(3-21)定义节点 i 的度 di 为连接到节点 i 的节点数目，如果不计节点的

52、度，di 是矩阵 K 中的 i行中除对角元素以外的非零项数。这里暗含着在四边形或六面体单元中处于对角的节点是相互连接的。最大的节点度为m maxdi , 1 i n独立的边数 e 定义为对角线以上非零元素的个数(3-22)1ne d(3-23)i2i1这样矩阵 K 中总的非零元素的个数为 2en，矩阵 K 的密度 2e n(3-24)n 2需要注意的是在i 和 ci 的定义中包含矩阵 K 的对角元素。根据这些量的定义，可以方便的将各种参数从一种表示（包括对角元素）转换为另一种表示（不包括对角元素）。同样地，这里所说的矩阵 K 的阶 n 有时就是指节点的数目。在一般的有限元分析中，每个节点有数个

53、度，例如如果每个节点有 6 个度，则在没有约束的情况下，wavg，wrms，wmax是仅考虑网格节点时对应的数值的 6 倍。353.6.2 来源于图理论的有关定义定义一个图 G = (V, E)，V = v1, v2, , v n是顶点(Vertex)的有限非空集，vi, vjE，i j 定义了称为边(Edge)的有限集，V = V(G)，E = E(G)。这里 G 是一个无向图（Undirected Graph）。如果vi, vj，则称顶点 vi 和 vj 是相邻接的。对于有限元分析来说，V 相当于节点集合，E 相当于节点之间的联系。图 G(V, E)的子图（Subgraph）S(V, E)

54、定义为 SG，VV，EE。一个顶点 v 的度（Degree）g(v)是与 v 相连的边数。一个长度为 t 的路径（Path）是边的序列v1,v2，v2,v3， vt-1,vt。两个顶点 u 和 v 之间存在一路径，则说 u 和 v 是连通的（Connected）。u，v 之间的距离（Distance） (u, v)是连通 u，v 两顶点的最短路径。如果图 G 中的任一顶点对是连通的，则图 G 是连通的，否则是不连通的。如果图 G 有 n 个顶点 V = v1, v2, , vn，一对一的f 将v1, v2, , vn为1，2，n，f 称为图 G 的标号（Numbering）。对每一个标号 f，

55、定义相对于 f 的图 G 的带宽为 f G max f v1 f v2 , v1, v2 E(3-25)一个图 G 的层次结构（Level Structure）L(G)是将图 G 按如下方式分为子集 L1，L2， Ld与层 L1 中顶点相邻接的每个顶点处于层 L1 或 L2 中；与层 Lj（1 j 0，做 f-i 步。选择将被标号的顶点。搜索优先队列，根据 pt = maxpj，jQ定位顶点 t，m 为顶点 t的索引，qm = t。更新队列和优先级。通过设置 qm = qn 将顶点 t 从优先队列中删除，nn 1。如果顶点t 不是预活动的，则进行下一步；否则检查与顶点 t 邻接的每个顶点 x

56、并按照 px = px + W2增加它的优先级（对应于将顶点 x 的度减小一个）。如果顶点 x 是不活动的，通过设置 nn + 1，qn = x 以预活动状态将 x优先队列中。h.标号下一个顶点。置顶点计数 ii + 1，用 i 标号顶点 t；置 rt = i，R 为新顶点标号表。置顶点 t 的状态为过活动的。i. 更新优先级和队列。检查与顶点 t 邻接的每个顶点 x，如果 x 不是预活动的，置为不活动；否则，置顶点 x 为活动状态，px = px + W2，检查与顶点 x 相邻接的每个顶点 y。如果顶点 y 不是过活动的，增加它的优先级 py = py + W2；如果 y 是不活动的，通过设

57、置 nn +1，qn = y 以预活动状态将 y优先队列中。j. 结束。新的顶点标号保存在 R 中，rt 是顶点 t 标号。上面的算法中权 W1 和 W2 反映了当前顶点与端部顶点的距离和当前顶点度的相对重要性，根据经验一般取 W1= 1 和 W2 = 2，并且 W2 W1。以上的重标号算法适用于使用带宽求解的求解方案。使用求解时，由于方程的集合和删除是以单元依次进行的，故必须决定一个高效的单元次序。按已优化的节点，以单元中最小节点的升序分类和标号单元，可以导出一个好的单元标号。这样做保证了删除方程是以类似于优化的节点的次序进行的，实践证明也是有效的。Sloan 算法的效果较好，且适应带宽解法

58、和波前解法，编程也较容易。3.6.3.4 混合方法Armstrong104混合算法先使用 GPS 算法将节点重标号，然后使用行列交换的的方法来改进标号的次序。Armstrong 观察到基于图理论的算法能快速给出一个粗略的优化节点标号，但不能象行列交换算法那样得到更好的、详细的节点。Armstrong 的混合算法可以得到好于 GPS算法的结果。3.7 小结本章对有限元网格生成技术，特别是三维非结构六面体网格生成技术进行了全面的总结，并对有关的技术进行了较为详细的描述。对有限元网格后处理所使用的技术，如网格节点位置平滑、带宽优化等技术进行了总结，并详细的描述了有关的算法。40K. Ho-Le, F

59、inite Element Mesh Generation Methods: A Review and Classification, Computer Aided Design, 1988, 20:2738.M. S. Shepherd, K. R. Grice, J. O. Lo and W. J. Schroeder, Trends in Automatic ThreeDimenal Mesh Generation, Computers & Structures, 30(1-2):42142924 D. A. Field, The Legacy of Automatic Mesh Gen

60、eration from Solid Ming, Computer Aided Geometric Design, 1995, 12:65167325，有限元网格生成方法发展综述，计算机辅助设计和图形学学报，1997，9(4):37838326，(2):2831，有限元网格自动生成典型方法及发展方向，计算机辅助设计和制造，1996，27 闵，有限元网格划分技术，计算机研究与发展，1995，32(7):3742，28 C. Gay, P. Montmitonnet, N. Soyris and J. L. Chenot, Validation of a 3-D Finite Ele ment P