（工程力学专业论文）三维四面体网格优化及其程序实现.pdf

大连理工大学学位论文 摘要 现有的全自动有限元网格生成的方法，如a f t ( a d v m c h l gf r o n tt c 曲n i q u c ) 、 跳u n a y 、八叉树法等，在生成网格时都不可避免地会产生形状恶劣的单元，导致有限 元数值计算精度得不到保证。三维四面体有限元网格剖分较之二维三角形网格剖分，出 现形状恶劣单元的几率要大得多，原因是四面体单元扭曲变形的形式比三角形单元要 多。相对于三角形网格，四面体网格质最差的单元所占的比例更大、严重变形的单元更 多。因此，三维四面体有限元网格优化是一个特别值得关注的问题。 本文的方法重点解决了三维四面体网格剖分结果中的薄元问题，很好的提高了最差 单元的质量。本文首先概要地介绍了国内外广泛采用的一些有限元网格优化方法，然后 着重研究了三维四面体网格的一种优化算法i 坤】a c i a n 光顺和广义薄元分解的结合算 法，最后给出了这种优化算法的几个算例。 薄元问题在三维四面体网格剖分过程中普遍存在，它的出现不依赖于具体的网格生 成方法。只是在不同的网格生成方法中，广义薄元生成过程的形式略有差别而已。由于 广义薄元的形态恶劣，难以满足有限元法对单元质量的要求，所以必须对其加以处理。 本文提出了一种有效的三维四面体网格质量优化方法。改进了薄元分解方法以更全 面地考虑各种劣质单元类型，能够对三维实体网格剖分中产生的各种类别的孤立劣质单 元进行有效的分解：将改进的薄元分解法与脚l a c i 跚光顺优化方法相结合，修正各种 网格剖分算法所产生的非孤立劣质单元，提高最差单元的质量，以满足有限元分析对网 格质量的要求。 在提出上述优化算法的基础上，作者用c 语言编制了相应的程序。经过实例检验， 本文提出的四面体单元网格优化算法健壮有效、易于实现，优化后的网格质量高。 关键词：簿元分解；l a p l 8 c i 髓光颓；优化：四面体 三维四面体网格优化及其程序实现 t h r e e - d i m e n s i o nt e t r a h e d r a lm e s h o p t i m i z a t i o nm e t h o d a n d p r o g r a m a c h i e v e d a b s t r a c t a u t o m a t i cm e s h g e n e r a t i o n m c t h o d ，s u c h a s a f r ( a d v a n c i n gf r o n tt e c h n i q u e ) d e h u n a yt r i a n g u l a t i o na l g o r i t h m , a n d o c t t e ea l g o r i t h mm a y g e n e r a t ep o o r l ys h a p e d o rd i s t o r t e d e l e m e n t s ，w h i c hc a t l s en u m e r i c a ld i f f i c u l t i e sd u r i n gt h es o l u t i o np r o c e s s a n df o rt h em e s h e s w i t h h i g h l yd i s t o r t e de l e m e n t s ，t h es o l u t i o ni sb o t hl e s sa c c u r a t ea n d m o r ed i f f i c u l tt oc o m p u t e 1 i sp r o b l e mi sm o r es e v e r ei nt h r e ed i m e n s i o n st h a ni nt w od i m e n s i o n s ，b e c a u s et e _ l r a h e d r a l c a nb ed i s t o r t e dt op o o r q u a l i t yi nm o r ew a y s t h a nt r i a n g l e s c o m p a r e dw i t ht r i a n g u l a rm e s h e s ， t e t r a h e d r a lm e s h e st e n dt oh a v eal a r g e r p r o p o r t i o no fp o o rq u a l i t ye l e m e n t sa n dt oh a v e e l e m e n t st h a ta l em o r e s e v e r e l y d i s t o r t e d an u m b e ro fa l g o r i t h m so ff i n i t ee l e m e n tm e s ho p t i m i z a t i o nu s e d e x t e n s i v e l y a r e i n t r o d u c e ds u m m a r i l ya tf i r s t t h e nt h e a l g o r i t h mc o m b i n i n g t h ee x t e n d e ds l i v e r d e c o m p o s i t i o n m e t h o da n dl a p l a c i a ns m o o t h i n gm e t h o di nt h r e e - d i m e n s i o n a lp r o n e ma r ci n v e s t i g a t e dw i t h m u c h e m p h a s i s b y t h ee n d ，v e r a lm e s h g e n e r a t i o ne x a m p l e su s i n g t h e a l g o r i t h m 1 h ep r o b l e mo fs l i v e re l e m e n to c c l r 苫c o m m o n l yd u r i n gt h et e t r a h e d r a lm e s h g e n e r a t i o n p r o c e s s f o rt h ed i f f e r e n tm e s hg e n e r a t i o nt e c h n i q u e ，t h eg e n e r a t i o nf o r mo f s l i v e re l e m e n t si s d i f f e r e n tb o t , a n s es l i v e re l e m e mi ss h a p e do rd i s t o r t e de l e m e n t , w h i c hc a u s en u m e r i c a l d i f f i c u l t i e sd u r i n gt h es o l u t i o n p r o c e s s ，s o i tm u s tb e r e p a i r e d t om e e tt h en e e do ff i n i t ee l e m e n ta n a l y s i s ，a l le f f e c t i v ea l g o r i t h mf o rt h r e e d i m e n s i o n a l t e t r a h e d r a lm e s ho p t i m i z a t i o ni sp r o p o s e di nt h i sp a p e r 日1 es l i v e rd e c o m p o s i t i o nm e t h o di s e x t e n d e dt od e a lw i t ha 1 1k i n d so f p o o r q u a l i t ye l e m e n t a l lt h ep o o rq u a l i t ye l e m e n t sb e l o w t h e s p e c i f i e dq u a l i t yf a c t o ri nf i n i t ee l e m e n tm e s h c a nb ee f f e c t i v e l yr e m o v e d a n dt h e n , an e w m e s h o p t i m i z a t i o n m e t h o d c o m b i n i n g t h ee x t e n d e ds l i v e rd e c o m p o s i t i o nm e t h o da n d l a p l a c i a n s m o o t h i n gm e t h o di sp r o p o s e dt oo v e r c o m et h ep r o b l e mo f u n - i s o l a t e dp o o rq u a l i t ye l e m e n t s f r e q u e n t l yo c c u r r e d i nf i n i t ee l e m e n tm e s h g e n e r a t i o n a c c o r d i n gt ot h ea b o v e - m e n t i o n e da l g o r i t h m , ap r o g r a mi si m p l e m e n t e di ncl a n g u a g e c o m p u m t i o n a le x p e r i m e n t ss h o w t h a tt h ea l g o r i t h m p m p o s e d i sr o b u s t , e f f i c i e n ta n d e a s y t ob e i m p l e m e n t e d i np r a c t i c a l a p p l i c a t i o n s k e y w o r d s ：s l i v e r d e c o m p o s i t i o n ：l a p l a c i a ns m o o t h i n g ；o p t i m i z a t i o n ：t e t r a h a d r a l 独创性说明 作者郑重声明：本硕士学位论文是我个人在导师指导下进行的研究 工作及取得研究成果。尽我所知，除了文中特别加以标注和致谢的地方 外，论文中不包含其他人已经发表或撰写的研究成果，也不包含为获得 大连理工大学或其他单位的学位或证书所使用过的材料。与我一同工作 的同志对本研究所做的贡献均已在论文中做了明确的说明并表示了谢 意。 作者签名：日期： 大连理工大学学位论文 1 绪论 1 1 有限元法 在近代工程科学技术的发展过程中，由于飞行器、船舶、车辆、机械、水坝、桥 梁、房屋等工程设计上的需要，固体力学始终受到人们的重视。在2 0 世纪4 0 年代以 前，人们只能采用变分法、差分法、松弛法等计算方法，这些方法用来分析形状简单的 结构尚可采用，对于实际工程中很复杂的结构很难进行比较精确的分析。到了4 0 年代 中期，电子计算机出现以后，人们想到采用计算机求解杆件结构力学中力法和位移法的 基本方程，形成了矩阵力法和矩阵位移法。在类似的思想的指导下，到了5 0 年代中期 人们提出了有限元法，将连续介质离散成一组有限个单元，使无限自由度问题转化成有 限自由度问题，再利用计算机求解。 有限单元法的出现，是数值分析方法研究领域内重大的突破性的进展。早先发展起 来的数值方法，如有限差分法、配点法、最小二乘法、力矩法等都在不同的领域内取得 了成功的应用，但是它们对于几何形状复杂的阿繇，求解精度很低，甚至难以使用。而 有限单元法则是将连续的求解区域离散为一组有限个、且按一定方式相互连接在一起的 单元的组合体。由于单元能够按照不同的连接方式进行组合，且单元本身又可以有不同 的形状，这样就可以模型化几何形状复杂的求解区域，并且随着单元数目的增加，单元 自由度的增加以及插值函数次数的提高，解的近似程度将不断提高。因此，有限单元法 可以用来分析十分复杂的结构，并且可以在计算机中模拟各种复杂的材料本构关系、荷 载和条件。 三十多年来，有限单元法获得了广泛的应用，显示出了极大的功效和旺盛的生命 力。有限元方法应用已经由弹性力学平面问题扩展到空间问题，板壳问题，由静力平衡 闯题扩展到稳定问题，动力问题和波动问题。分析的对象从弹性材料扩展到塑性、粘弹 性、粘塑性和复合材料等，从固体力学发展到流体力学、传热学、电磁学等连续介质力 学领域。在工程分析中的应用已经从分析和校核发展到优化设计并和计算机辅助设计相 结合。随着计算机水平的发展和计算机辅助前处理技术和计算机辅助后处理技术的发 展，以及有限单元法与计算机的结合，有限单元法必将会获得更加广泛的应用。 1 2 计算机的发展 计算机产生于2 0 世纪4 0 年代，在计算机应用的初期，由于体积庞大，造价昂贵， 操作繁琐，应用仅仅局限于政府、军事以及商务领域，- - d 部分应用在科学计算领域。 三维四面体网格优化及其程序实现 随着科学技术的发展，计算机的体积不断减小，价格不断下降，自身性能却得到大 幅度的提高，由原来的每秒钟运算5 0 0 0 次到现在的每秒钟运算千兆次。随着计算机硬 件技术的发展，软件行业也发生了翻天覆地的变化，操作系统日趋成熟，新型的操作界 面对用户更加友好，同时许多应用软件也如雨后春笋一般纷纷出现，它们可以满足用户 的大部分需求。这样，计算机的价格和可应用范围已经到了广大用户可以接受的程度， 对计算机的操作和使用不再仅仅局限于受过特殊训练的专家和操作员，普通的用户稍加 培训即可对计算机进行操作，完成自己所要完成的任务，因此，计算机得到了广泛的应 用，包括工程计算、航空航天、信息管理、财务管理、生物医学、工业控制等等。其 中，在工程领域里，计算机的应用前景尤为广阔，如计算机辅助设计、计算机辅助制 造、计算机辅助工程等，而与有限单元法关系最为密切的前处理部分即可与计算机辅助 设计结合起来，完成手工难以完成的复杂结构的实体造型、网格生成和分析计算。 1 3 几种最常用的网格剖分算法 ( 1 ) 映射法 映射法既是结构化网格生成方法。又是非结构化网格生成方法。它的基本步骤是； 通过适当的映射函数将待剖分物理域映射到参数空间中形成规则参数域；对规则参数域 进行网格剖分；将参数域的网格反向映射回物理空间，从而得到物理域的有限元网格。 映射法可分为保角映射法( c o n f o r m a lm a p p i n gm e t h o d ) 、基于偏微分方程法( p d e - - b a s e d m e t h o d ) 以及代数插值法( a i g e b r a i ci n t e r p o l a t i o nm e t h o d ) _ 三大类。保角映射法能够处理多于 4 条边界的单连域问题，但难于控制单元形状和单元密度，且不能直接应用于三维问 题。t h o m p s o n 等首次提出的基于偏微分方程法，通过数值求解偏微分方程而得到参数 空间与物理空间之间的映射关系。基于偏微分方程法可细分为椭圆型、抛物线型和双曲 线型三大类，其中最常用的是椭圆型p o i s s o n 方程和l a p l a c e 方程网格生成方法。基于 偏微分方程法主要应用于c f d ( c o m p u t a t i o n a lf l u i dd y n a m i c s ) 领域的贴体坐( b o u n d a r y f i t t e dc o o r d i n a t e ) 网格生成问题，目前仍在深入研究和发展中。代数插值法是一种应用 相当广泛的网格生成方法，它通过代数插值来描述参数空间与物理空间之间的映射关 系，其中最重要的一个类别是超限插值( 1 r a n s f i n i t ei n t e r p o l a t i o n ) ，z i e n k i e w i c z 提出的等 参映射就是超限插值的一种特殊形式。超限插值既可适应特殊的区域边界形状，又可控 制所生成单元的形状和密度，成为代数插值网格生成的一个标准方法。映射法的优点是 算法简单、速度快、单元质量好、密度可控制，它既可生成结构化网格又可生成非结构 化网格；既可生成四边形单元网格又可生成六面体单元网格，可用于曲面网格生成，可 与形状优化算法集成等。因此，映射法在众多的商业有限元分析软件中占有重要的地 2 大连理工大学学位论文 位。映射法一般可直接处理单连通域问题，但对于复杂多连通域问题，需要首先用手工 或自动方法将待剖分域分解成几何形状规则的可映射子区域，然后在每个子区域内应用 映射法。 ( 2 ) d e l a u n a y 三角化方法 d e l a u n a y 三角剖分( 简称d t ) 是目前最流行的通用的全自动网格生成方法之一。 d t 有两个重要特性；最大一最小角特性和空外接圆特性。d t 的最大一最小角特性使 它在二维情况下自动地避免了生成小内角的长薄单元，因此特别适用于有限元网格生 成。所谓空外接圆特性，就是d t 中的每个三角形单元或四面体单元的外按圆( 二维) 或外接球( - - 维) 都不包含其它节点，b o w y e rw a t s o n 算法正是利用了这一特性。三维 b o w y e rw a t s o n 算法的基本步骤：首先定义一个包含所有节点的初始网格，最简单的情 形是单个四面体；向网格中插入一个节点，找出其外接球包含此节点的所有四面体单 元，删除这些单元形成一个包含插入节点的空腔 晒劝：将该插入节点与空腔的每个表 面相连，形成新的四面体网格；重复进行上述的节点插入过程，直到全部节点插入完 毕。 d e l a u n a y 三角剖分算法的计算效率与具体的实现相关。大体上可将d t 算法分为 三大类：分治算法，逐点插入法和三角网生长法。b o w y e rw a t s o n 算法是一种典型的逐 点插入法，其时间复杂度为o ( u 3 ”) ，采用四( 八) 叉树数据结构的b o w y e rw a t s o n 算 法可达到o ( n l o g n ) ，分治算法的时间复杂度为o ( n l o g n ) ，三角网生长法的时间复杂 度为o ( n 3 p 2 ) ，其中n 为生成单元总数。经典d t 技术已经相当成熟，近年来的研究 重点是约束d t 的边界恢复算法，以及如何克服b o w y e rw a t s o n 算法退化现象所产生的 薄元( s l i v e rd e m e n t ) 问题。 ( 3 ) 推进波前法 经过近年来的发展，推进波前法( a d v a n c i n gf r o n tt e c h n i q u e ，a f t ) 已成为目前最流 行的通用的全自动网格生成方法之一，该方法的提出应归功于l o h n e r 和l o 。a f t 方 法没有d e l a a n a y 三角剖分算法那样成熟的理论依据，在很多情形下a f t 方法是靠经验 解决问题，但是这并不妨碍它的成功应用。a f t 算法的时间复杂度为o ( n 1 0 9 r 0 ，与 d e l a u n a y 三角化算法、有限四( 八) 叉树法相当，但其生成单元的质量是三者中最好的。 a f t 方法的基本流程是：首先离散待剖分域的边界，二维待割分域的边乔离散后 是首尾相连的线段的集合，三维待剖分域的边界离散后是拓扑相容的三角形面片的集 合，这种离散后的域边界称为前沿；然后从前沿开始，依次插入一个节点，并连接生成 3 三维四面体网格优化及其程序实现 一个新的单元；更新前沿，这样前沿即可向待剖分域的内部推进。这种插入节点、生成 新单元、更新前沿的过程循环进行，当前沿为空时表孵整个域剖分结束。 1 4 生成网格的质量 网格质量好坏的直接度量标准就是有限元分析模块利用此网格进行求解是否精确， 以及求解速度的快慢。质量好的网格可以为分析阶段提供比较好的单元剐度矩阵，因此 也将得到比较准确的解。 网格质量的评价标准有很多，例如可以根据生成单元的几何形态来度量，常用的参 数有单元最小边长与最大边长的比值，最小角度与最大角度的比值以及单元平面的翘曲 程度等，这些参数在网格生成后评价网格质量方面是非常重要的。 网格的局部单元质量同整体单元质量一样都是不容忽视的，即使整体网格平均质量 比较高，如果有一些局部单元的质量很坏( 如扁平单元) ，或出现拓扑结构错误的情 况，那么必将导致分析阶段的矩阵形态恶劣乃至错误，难以计算出正确的解。 另一个有关网格质量的度量方式就是检查边界处网格的质量和着重分析处嘲格的质 量，这些地方的网格质量也将在很大程度上影响分析结果。 1 5 课题的研究意义 有限元网格生成是工程科学与计算科学相交叉的一个重要研究领域，在经历了3 0 多年发展后的今天依然十分活跃。一方面，有限元法已成为一种能够有效地求解各类工 程和科学计算问题的通用数值分析方法；另一方面，计算机硬件运算能力的不断提高也 容许人们对工程和科学计算的规模、复杂度、效率、精度等方面提出更高的要求。 为了保证有限元分析结果的精度，人们总是希望获得高质量的网格。对几何形状较 为简单的区域，可以采用半自动方式生成网格，在剖分的同时能够直接控制网格的质 量。由于所要分析的问题越来越复杂，工程人员逐渐倾向于使用网格自动生成工具来剖 分网格。但自动生成技术得到的初始网格( 三维情况) 往往质量很差，后续的网格质量 优化就十分必要了，因而网格优化已成为目前网格生成领域中的一个重要课题。而三维 四面体有限元网格剖分较之二维三角形网格剖分，出现形状恶劣单元的几率要大得多， 原因是四面体单元扭曲变形的形式比三角形单元要多。因此，三维四面体有限元网格优 化是一个特别值得关注的问题。 1 6 本文的主要工作 本文的主要工作有以下三个方面： 1 针对传统l a p l a c i a n 光顺算法存在的不足提出了改进的l a p l a c i a n 光顺算法。 4 大连理工大学学位论文 2 定义了广义薄元( 本文中广义薄元的质量系数小于0 0 2 ) ，将广义薄元分成四种不 同的类型，并提出了四种相应的广义薄元分解算法。 3 将改进的l a p l a c i a n 光顺算法和广义薄元分解算法结合使用对四面体网格进行优 化。 四亟佳围揸篮匙凶蔓翅塞皇塞瑰属计算工程与科学计算可视化研究领域，课题 源于： ( 1 ) 国家自然科学基金创新研究群体科学基金( 1 0 4 2 1 0 0 2 ) 计算力学与工程科学计算， 2 0 0 5 2 0 0 7 。 ( 2 ) 国家杰出青年科学基：金( 1 0 2 2 5 2 1 2 ) ：多相多场藕合作用下多孔材料的非线性与失效分 析，2 0 0 3 - 2 0 0 6 。 ( 3 ) 大连市科学技术基金：面向机械零组件精细化设计的三维参数化形状优化方法， 2 0 0 4 2 0 0 5 。 - 5 - 一三壅婴垩堡翌巡垦薹堡壁塞里 ，_ _ ，_ _ - - _ - _ ，_ _ _ 一 一 现今，提高非结构化网格质量的最基本方法【l 】有三种： ( 1 ) 插删除点； ( 2 ) 拓扑优化； ( 3 )网格光顺。 2 1 插入删除点 网格单元中，往往要衡量单元某一边与边中点处期望长度的比例，二维中这个比例 的理想值是1 5 3 ，三维网格中的理想值是2 0 。如果这个比例太大就要将边分解，插入 点可以用来改变这个比例，并且通过仔细选择点插入的位置可以控制单元的几何形状， 达到单元优化的效果。 删除点的操作主要是作为种重定义技术。许多实例当中，通过删除一个区域的网 格完成删除点的操作，然后重新选择长度比例对这个空腔进行网格重新生成。点删除技 术的不同应用也可以用来改善单元的几何形状达到网格优化的目的。 2 2 嗣格光顺 网格光顺是通过改变节点坐标的方式来提高有限元网格整体质量的一项技术。 对于某一个节点来讲，有两方面的属性这里比较关心。 一个属性是这个节点所处的位置。一个节点在最后生成的网格中总的来讲有三个位 置，一个是节点在模型的内部，称之为内部节点；一个是节点在模型的表面，称之为表 面节点；最后一个是节点在模型的外表面的交接线上，称之为边界节点。 网格光顺所能够移动的只能是内部节点和表面的节点，其中表面节点也不是随便可 以移动的，只有处在一个平面上的那些表面节点才可以移动。 另一个属性是所有与这个节点相连接的节点组成一个多边形( 二维) 或多面体( 三 维) 。而最常用的l a p l a c i a n 光顺是将某一节点移动到与之相联系的多边形或多面体的 形心处。这种移动方式对于二维问题获得了极大的成功，但是对于三维问题却未必如 此。在三维问题的处理中，节点的移动可能导致新近生成的单元出现病态，即新近生成 的单元出现畸形。因此，在三维问题中应用变换必须逐步进行，不可能一步到位。那 么，如何应用变换逐步进行单元网格质量的提高呢? 假设k 点进行移动之前笛卡尔坐标为( x ，y ，z ) ，应用变换得到的笛卡尔坐标 为( x ，y ，z ) ，两者的距离为d ，将k 点分别移动至1 3 4 、d 2 、3 d 4 、d 处，分 ，6 - 大连理工大学学位论文 别检查与k 相联系的单元的质量，保留新生成的单元的质量得到提高的最多的移动后 节点坐标。用这种方法在三维问题中应用变换可以得到比较令人满意的结果。 2 2 1 传统的l a p l a c i a n 光顺 传统l a p l a c i a n 光顺做法是移动目标节点到相邻节点的中心。这种方法，算法简 单，但却存在着很大不足：计算量大，单元质量还不一定会得到提高，最致命的是容易 产生不合法单元( 单元反向或者单元体积为负) 。 2 2 2 智能l a p l a c i a a 光顺 智能l a p l a c i a n 光顺不同于传统l a p l a c i a n 光顺的是 当节点移动到相邻节点的几何 中心时，若网格质量得不到提高，则移动操作不予执行。遗憾的是，智能l a p l a c i a n 光 顺一样会产生不合法单元。 2 2 3 基于最优的光顺( o p t i m i z a t i o n - b a s e ds m o o t l l i n g ) 在l a p l a c i a n 光顺之后出现了一种新的网格光顺方法受到了越来越多的关注，它就 是o p t i m i z a t i o n - b a s e d 光顺。节点移动不再采用l a p l a c i a n 光顺中的些做法，而是采用 对目标函数求最优解的方法确定节点的新位置。 o p t i m i z a t i o n - b a s e d 光顺算法的目标函数是( 3 1 ： ( 工) 2 器鲁，( x ) ( 2 1 ) z ( 算) 是单元质量函数，庐( x ) 是日标函数( 节点相邻单元集中最小单元质量溺 数) 。妒( x ) 的优化解就是节点移动后的位置坐标。 对于光顺函数，使用一个类似最速下降法的方法来解决非光顺优化问题。通过解一 个二次规划问题，根据x 处的所有可能梯度的凸线性组合给出最速下降方向，并且计算 每一步的搜索方向。通过预测激励方程将会改变的点和基于一阶t a y l o r 级数的取) 的近 似值，可以解沿着最速下降方向s 的线性搜索予问题。设鬻初始步长口为从当前位置到 激励方程被预测会改变的点的距离。标准步长和迭代终止准则保证了算法的成功应用。 下面给出了o p t i m i z a t i o n - b r e e d 的光顺算子口】。 i - 0 计算f ( x 。) 和a ( x 0 ) w h i l e “一x ) a n d ( 岱 m ns t e p ) ) a n d ( i m i n _ s t e p ) 计算x = x + l 计算厂( t + ，) 和a ( t + 1 ) 尝试可接受步长口= d 2 e n d w h i l e i 鼍+ l e n d w h i l e 这种方法同a 瑚曲t a 即等推出的线性规划法是一样的，因此凸规划准则也可以用来 判断是否有唯一的解x 。a m e n t a 等通过对几个单元质量的描述，发现它们中的许多都 是凸域，有唯一的解。有很多人( 如：b a n k 、s m i t h 、s h e 曲a r d 和g e o r g e s ( 1 9 9 1 ) 等) 提出了相似的局部o p t i m i z a t i o n - b a s e d 光顺算法，但这些算法在最优求解过程和单元质 量评价标准方面都是不同的。 在f r e i t a g t l 3 1 等的文章中，实例检验都表明了o p t i m i z a t i o n - b a s e d 光顺相比于 l a p l a c i a n 光顺算法，无论是在二维还是三维的简单网格单元中的应用都是更有效的。 o p t i m i z a t i o n - b a s e d 光顺在从网格中删除有极限角的单元的操作中是非常有效的，而 l a p l a e i a n 光顺在处理变形恶劣的单元时却显得无能为力。 o p t i m i z a t i o n - b a s e d 光顺计算代价比l a p l a c i a n 光顺要大，但是这种方法的优化结果 更好，特别是对于有凹域的情况。但是，大多数人还是喜欢用l a p l a c i a n 光顺，只是在 很必要的时候才使用o p t i m i z a t i o n - b a s e d 光顺。 c o u g r l 1 研究了最初的o p t i m i z a t i o n b a s e d 光顺算法。他的方法能够提高八叉树法 生成的四面体网格结果中边界附近的扭曲变形的四面体单元的质量。单元的变形程度由 单元的体积和单元表面面积的比例来衡量，针对这种比例他提出了四种属性，来证明这 种判断准则非常适合o p t i m i z a t i o n - b a s e d 光顺。然后又找到了使最大变形尺寸最小的节 点的位置，并且在当前位置和最优节点位置实旌了线性查找。 p a n m y 【4 】也研究了一种针对三角形和四面体网格优化的o p t i m i z a t i o n - b a s e d 光 顺算法，用来修复四叉树和八叉树方法生成的网格单元。对非线性约束的全局优化问题 进行求解，目标函数是单元面比例的最小值，用不等式约束来确保单元的表面面积( 三 维是体积) 都不小于某一( 初始) 值。用改进的可行方向法来进行优化问题求解。 8 一 大连理工大学学位论文 p a r t h a s a r a t h y 发现：求解这个全局优化问题与在网格生成阶段所花费的时间一样多。但 这种光顾算法所得到的优化结果是令人振奋的。 c a n a n n 4 研究了一种全局优化方法，这种方法用的是o d d s 变形准则。虽然这种方 法主要适用于六面体网格，但也很容易扩展到其它的网格单元类型。像p a r t h a s a r a t h y 研 究的方法一样，这种方法也不适用于大网格，因为这两种方法用的都是全局优化算法。 递归的局部优化算法被证明是更可行的。 f r e i 斛研究了一种用最大最小角为准则来优化三角形和四面体单元的方法。优化 的目标函数是一个非连续、非光顺的，用类似最速下降法求解的函数。这种方法被实例 证明是可行的，而且效果好、合理有效。f r e i t a g 也在实际应用中将此方法和l a p l a c i a n 光顺方法结合使用，收到了不错的优化结果。 2 2 4b o u n d a r yl a y e r 光顺 a n s y s 提供了生成2 d 自由两格的功能，而在几何边附近用户往往设定一定层数 的尺寸非常小的单元，在向内生成网格时单元尺寸快速增加到很大。这类单元主要 用于c f d 边界层或者e m a g 表皮效应的研究。为了保证边界层单元光滑过渡到其余 的单元，必须要对光顺过程做一定的改进。在光顾之前，标记离几何边最远的边界 层节点。在光顺过程中，这些节点始终保持离几何边有一个相同的常量边界层距 离。在光顺过程中，使用l a p l a c i a n 光j 暖算法只是对外层节点进行光顺，把节点投影 到几何边上，然后垂直从几何边开始移动节点一定的距离，这个距离等于边界层的 厚度。 2 2 5p h y s i c s - b a s e d 光顺 质量越好的网格应该越接近自然界中的实物，多名权威人士为了解决简单的物 理问题已经研究出了多种光顺算法。l o h n e r 等研究了一种光顺算法，算法中将网格 看成两节点问弹簧连接的系统，这样节点间的力是实际的网格尺寸和理想尺寸的比 值。这种技术会产生基于密度梯度的拉伸单元。 s h i m a d a 和b o s s e n 研究了一种网格划分技术，包括网格生成、网格光顺和质量差 单元的清除，这些操作都是p h y s i c s - b a s e d 方法【4 】的部分。s h i m a d a 将节点视为腔体 的中心，每一个腔体和其它的腔体相交。b o s s s n 研究发展t p l i a n t 方法，这种方法和 s h i m a d a 的方法类似，是用节点间的引力和排斥力来决定节点的位置。这两种方法都 不能算是真正意义上的光顺技术，因为它们在对质量差单元进行优化时都要进行节 点插入和删除。每一种p h y s i c b a s e d 技术都能够保证在各向同性网格中有拉伸单元。 - 9 * 三维四面体网格优化及其程序实现 2 3 拓扑优化( t o p o l o g i c a lo p t i m i z a t i o n ) 拓扑优化算法是一种通过改变网格的拓扑结构来提高网格质量的算法，无论是在二 维三角形网格中还是在三维四面体网格中都有过成功的应用。 2 3 1 边交换( e d g es w a p p i n g ) ( 1 ) 二维三角形网格模型中边交换的应用 二维三角形网格中，一条内边最多被两个三角形单元共享。因此在三角形网格中， 边交换通常用来对共事一条内边的两个三角形进行局部重构，重构后的单元个数不会发 生改变，仍然是两个。 三角形单元可以用三角形的内切圆半径的3 倍和外接圆半径的比值作为单元的质量 评价标准。如果共享一条内边的两个三角形都是质量差单元，那么通过使用边交换可能 会提高网格单元的质最。图2 1 是一个三角形网格中边交换的图示： a b b 图2 1 三角形网格中边交换示意图 f i g 2 1s k e t c hm a po f e d g es w a p p i n g i n t r i a n g u l a rm e s h c 三角形单元a b c 、a c d 共享内边a c 。删除共享边a c ，连接b d ，初始的两个三 角形单元a b c 、a c d 被新单元a b d 和b c d 取代。 ( 2 ) 三维四面体网格模型中边交换的应用 文献 2 】，边交换通常可以用来对共享一条边的质量差的四面体进行重构。如果 有k 个四面体都包含同一个内边e ，那么边e 被删除并且初始的k 个四面体被2 ( k - - 2 ) 个新的四面体替代，替代要采取如下的方式：k 个四面体中除了e 以外的其 它边构成一个k 多边形。这个k 多边形可以分解成k - 2 个三角形。 可以用二面角和固体角作为四面体网格的质量评价标准。如果上述的几个四面 体都是质量差单元，那么通过使用边交换可能会提高网格单元的质量。图2 2 是一 个三维四面体网格中边交换的图示： 一l o 大连理工大学学位论文 3 图2 2 三维网格中边交换示意图 r i g 2 2s k e t c hm a po fe d g es w a p p i n gi nt h r e e d i m e n s i o nm e s h 最初，五个四面体0 1 b t ，1 2 b t ，2 3 b t ，3 4 b t , 和4 0 b t ，共享边t b 。删除共享 边b t ，连接0 2 和2 4 ，初始的四面体0 1 b t ，1 2 b t ，2 3 b t ，3 4 b t ，和4 0 b t 被六个新 的单元0 1 2 t ，0 2 1 b ，0 2 4 t ，0 4 2 b ，2 3 4 t ，和2 3 4 b 取代。 通常，如果一条边被k 个四面体共享，那么这k 个四面体会被2 k - 4 个四面体取 代。然而，当k 很大时，这种重新构造的局部单元数量会比初始单元数量多很多， 因此当k 值大于7 时不会使用这种边交换方法。 边交换方法有凡个不足。每一个初始的结构，特别是当k 很大时，要对比多个 结构来决定哪一个是最优的，这是项计算量很大的过程。由于几何存在对称和冗 余，使得对比算法节省了不少时间。然而，这样的对比算法却会占用一定的内存。 因此，边交换算法在计算时间和内存占用之间也存在一个平衡关系。 2 3 ，2 面交换( f a c es w a p p h a g ) 面交换也是通过改变单元的局部连接关系来提高网格质量的，主要是用在四面 体网格结果中。五个点，两个四面体共享一个内部面，这是一个典型的使用面交换 提高网格质量的结构。可能许多非重叠的四面体结构是由这五个点构成的，但仅仅 有两个四面体单元可以被合理的重新连接或者交换。如图2 3 所示，左图，两个或 者三个四面体单元能够分割这个五个点构成的凸域。要把这个由两个四面体组成的 凸域转换成三个四面体构成的结构。需要在这个凸域中额外加入一条内边。右图， 是一个两个四面体的结构，能够转换成不同的两个四面体构成的结构。图中的阴影 部分是共面的两个三角形，交换算法交换了这个四边形的对角线。这两个共面的三 角形所在的面如果不是边界面或者是被另外一对四面体共享的面，就可以使用交换 算法进行这种2 2 转换。否则，这种交换方法所产生的新边会使得网格不相容。 三维四面体网格优化及其程序实现 图2 3 面交换示意图 f i g 2 3s k e t c hm a po f f a c es w a p p i n g 因为每一个例子都仅有两个合理的形式，但是对比发现却都只有一种的质量是 更好的。如果在网格结果中存在的不是质量更好的这种形式，就要进行重新连接得 到更好的那种结构形式。如果两种形式的质量是相等的，我们选择两个四面体构成 的形式，如果再进行2 屯转换是无意义的。 决定是否要对一个四面体单元进行局部重新连接，可以用两种不同的质量判断 准则：最大最小二面角准则和球形准则。最大最小二面角准则适用于有五个节点， 且仅有的两个四面体共享一个面的结构。这种准则使得两个四面体中的目标四面体 的最大二面角最小。球形准则适用于组成四面体的四个节点的外接球中不包含其它 节点的四面体。感觉上这种准则有点类似于d e t a u n a y 三角化网格划分方法中的插入 新节点的准则。但是无论对于那一种准则，通过面交换实现优化的算法都是在局部 实施而不是全局。最近j o e 又提出来一种更高级的通过局部变换来提高网格质量的 技术。然而这种方法是否能够通过一系列的局部变换达到全局优化还不得而知。 2 4 一些其它的优化算法 2 4 1t d 法( t r a p e z i u md r a w i n g ) 文献 2 8 】，这是一种新型的简单、快速的网格优化算法。这种算法易于实现，且是 一种迭代算法，主要应用于三角形表面和平面网格。 p 。是初始网格的一个节点，p = p ， ( 净1 ，k ) 是和节点p 0 相邻节点的集合。对于节 点集合中的每一个节点p ，节点p 。都可以根据如下过程获得新的位置： 一1 2 大连理工大学学位论文 首先，找到节点p j p ，使得向量p o a 和p 。p j 包含一个最大角( 如图2 4 ) 。节 点p 。相对于p ，可以找到一个新位置p p 。，是这个四边形的一个顶点，使得三角形 p ，p 。，p ，是一个等腰三角形a 图2 4 查找使向量p 。p ，和 p o p j 夹角最大的节点 f i g 2 4s e a r c hf o rt h en o d es u c h a sv e c t o r sp o p ，a n d p o p j c o m p o s e am a x i m u ma n g l e 圈2 5 查找节点p o 相对于 节点p 。新位置p 。， f i g 2 5s e a r c hf o rt h en e wl o c a t i o n o ft h en o d ep fw i t hr e g a r dt o t h en o d e p o 。 这里会重点强调算法的两个方面。对于p o p , h p o p j ，新的节点p 。离初始的近似连 续光滑表面很近，而三角形p ，p 。p 是等腰三角形，这样就使相应的网格单元质量得到 了提高。 用p ，表示线段p fp j 的中点。新节点p o ，的坐标可以通过如下的公式求得： 图2 6 查找节点p 。，的坐标 脚2 6s e a r c hf o r t h ec o o r d i n a t e so f p ”矧。尚嘞+ c o 网。尚确+ 专铲。尚 - 薹 三维四面体网格优化及其程序实现 ( 如图2 6 ) 当所有节点p 。，的坐标值都已经求到