第21卷 第3期应用力学学报Vol.21 No.3 2004年9月CHINESE JOURNAL OF APPLIED MECHANICS Sep.2004文章编号:1000 4939(2004)03 0125 04任意平面区域的高质量有限元网格自动剖分李梅娥 周进雄(西安交通大学 西安 710049)摘要:实现了任意平面区域的高质量三角形有限元网格自动剖分。
初始化时利用广度优先搜索查找孔及凹槽中需删除的三角形,不需定义有向边界,可十分方便地定出实际剖分区域;给出了一种方便、快捷的点定位法,可大大提高程序运行速度;利用Delaunay网格优化算法实现了对三角形形状的控制,可使三角形最小角达到35 ;分割坏三角形时,按最小角由小到大的顺序处理,使三角形网格在满足质量要求的同时网格数量较少。
关键词:高质量三角化;Delaunay优化;任意平面区域中图分类号:TB115 文献标识码: A1 引 言在有限元分析过程中,三角形有限元网格应用最广,如何快速、高效地生成高质量三角形网格是人们研究的重要课题之一。
衡量三角形网格质量的常用指标有三角形最小角、最大角、边长比等。
在有限元计算中,三角形网格最小角过小,会增大离散误差,降低求解精度,甚至导致病态代数方程组,因此,三角形形状应尽量接近正三角形,避免尖角。
最早可以控制三角形形状的算法是由Baker等人[1]提出的,该算法可保证生成的三角形网格没有钝角,即所有角度都不超过90 ,且最小角不小于30 。
随后又出现了一些基于四叉树的三角形形状优化算法[2]。
上述方法都使用了背景网格,且生成的网格数量较大。
1995年,Ruppert提出了一种完全不同的网格形状优化方法,称为Delaunay优化(Delaunay refinement)[3],这是因为该算法在优化过程中仍保持Delaunay三角化性质,选取插入点位置时也遵照Delaunay三角化准则。
该方法具有以下特点:不需要背景网格,算法简单、易实现;在满足相同网格质量要求的前提下,生成的网格三角形数量少;生成的网格没有方向性。
本文利用Delaunay优化算法实现了对三角形形状的控制,并在确定计算区域、点定位等几个方面进行了改进,实现了任意复杂平面区域的高质量有限元网格剖分。
2 初始三角形网格形成平面直线图PSLG(Planar Straight Line Graph)是任意平面区域的归一化描述,因为任何二维曲线均可以使用小直线段进行拟合。
PSLG由一个点集和线段集组成,且线段的端点必须包含在点集中,这些线段称为约束边。
本文的有限元网格自动剖分就以PSLG描述的任意平面区域作为研究对象。
初始三角化就是将PSLG点集内的点作为三角形的顶点,在约束边界给定的区域内形成一组很粗的三角形。
首先采用Lawson逐点插入法(Incremental In sertion Algorithm)将PSLG点集内的点联接成Delaunay三角形网,即无约束Delaunay三角化。
其次,插入PSLG中的约束边线段。
较常用的方法是通过对角线交换去掉所有与该约束边相交的基金项目:国家自然科学基金项目(项目编号10202018) 来稿日期:2003 01 10 修回日期:2003 04 11第一作者简介:李梅娥,女,1972年生,西安交通大学材料科学与工程学院讲师,博士;研究方向:铸造凝固过程宏、微组织数值模拟三角形棱边,最终使得该约束边成为三角形网中的一条或几条边,并恢复约束Delaunay三角化[4]。
上述过程都是在PSLG给定点集的凸包中进行,但实际剖分区域可能并不是凸多边形,而且多边形内部可能有孔,因此,需删除给定剖分区域以外的多余三角形,包括孔及凹槽中的三角形。
本文采用广度优先搜索查找孔及凹槽中的三角形,即对任一多余三角形,分别判断它相邻的三个三角形是否也为多余三角形,再分别对这三个三角形的相邻三角形都进行判断,重复这一过程,直到遇到约束边。
因此,对孔而言,只需给出孔中某一点的坐标,找到该点所在的三角形,即可找出孔中所有三角形;对凹槽的处理也是如此,对凸包边界上某一三角形而言,若它在凸包边界上的边不是约束边,表明该三角形需删除,从该三角形开始,可找出该处凹槽中的所有三角形。
采用上述方法不需定义有向边界,可十分方便地定出实际剖分区域。
3 点的定位在三角化过程中,判断某一点位于网格的哪一个三角形中是程序中占用时间比例较大的操作之一,因此其运行效率对程序整体速度十分重要,如果每定位一个点时对所有三角形都进行判断,势必浪费大量时间,为此本文提出了一种比较好的点定位算法。
点定位时,从某一已知三角形开始查找,进入一个三角形后,可以从另两条边离开这个三角形,由哪一条边离开采用如下方法判断:如图1,若由BC边进入ABC,先判断P与边AB、AC的关系,若点P 位于有向线段AB的右侧,则由AB边离开,如P1点;若点P位于有向线段C A的右侧,则由AC边离开,如P2点;若上述两条件均满足(如P3点),则经过A点作一条垂直于BC的直线,根据点在该垂线的哪一侧来确定由哪一条边离开。
由于P3位于垂线右侧,因此,从AC边离开。
图2为查找点P所在三角形的过程示意图。
由于先后插入的点多数情况下相距较近,为提高查找效率,将最后找到的三角形记录下来,作为下次查找的起始三角形。
4 网格优化方法初始三角形网格形成以后,需要对网格进一步优化,使其满足形状和尺寸要求。
图1 点定位示意图图2 查找P点所在的三角形4 1 Delaunay优化算法在任一三角形网中,若某结点位于以某约束边界线段为直径的圆内或圆上,而且该点对该约束边可见,就称该点侵入(encroach)了该约束边,该约束边称为被侵入边(encroached segment)[3]。
Delaunay优化算法是通过不断插入结点,使所有三角形都满足用户指定的形状或尺寸要求。
在插入结点时根据以下两点进行:1) 如果一条约束边被侵入,就在其中点处插入一个结点,将该约束边分割为两条较小的约束边,并通过对角线交换恢复Delaunay三角化(如图3,其中粗实线代表约束边);2) 如果某三角形形状或尺寸不满足要求,称该三角形为坏三角形,在其外接圆圆心处插入结点,并恢复Delaunay三角化,如图4。
由Delaunay三角形空外接圆性质,该操作必然使得该三角形被删除;但如果插入点侵入某约束边,就不插入该点,而是分割所有被该点侵入的约束边。
图3 分割被侵入约束边图4 去除坏三角形可以证明[5],如果T是所有边界由约束边组成的某平面区域的Delaunay三角剖分,且T没有被侵入约束边,则T中任意三角形的外接圆圆心必然位于T 内。
因此,分割被侵入约束边实质是清除了外接圆圆心在剖分域外的三角形,使得随后处理坏三角形126应用力学学报第20卷时,插入的结点都在剖分域内。
所以,在该算法中,被侵入约束边比坏三角形优先处理。
4 2 网格优化的实现本文根据Delaunay优化算法的思想,提出了一套网格优化方法,具体实施步骤如下所述。
4 2 1 第一步 分割被侵入约束边1) 初始化被侵入约束边encroached segment链表;2) 检查所有约束边线段,判断其是否被侵入,若是,加入encroached segment链表;3) 从表头取出一条被侵入边进行分割,对分割后的两个小线段也要判断其是否被侵入,若是,也加入encroached segment链表尾部;4) 重复执行3),直至链表为空;5) 重复2)~4),直到不再有被侵入边。
4 2 2 第二步 分割坏三角形1) 初始化坏三角形bad triangles链表;2) 检查所有三角形,判断其最小角或面积是否符合要求,若不符合,根据其角度大小插入到链表中合适位置,链表根据角度由小到大顺序排列,即角度小的优先处理;3) 若bad triangles链表不为空,从表头取出一个三角形t;4) 判断t是否仍在三角形网T中,若不在,返回3),否则,执行5);5) 判断t的外接圆圆心是否侵入某一条或几条约束边,!若没有约束边被侵入,就在t的外接圆圆心处插入一个结点,然后依次检查该点相邻的三角形是否符合质量要求,若不符合,将该三角形插入bad triang les链表;∀若有,将t插入bad triangles 链表中,等待重新处理,并将这些被侵入约束边放入encroached segment链表;6) 若encroached segment链表不为空,就执行以下操作:!从表头取出一条被侵入约束边,在其中点插入一个结点,分割该约束边;∀依次检查该新插入点相邻的三角形是否符合质量要求,若不符合,将该三角形插入bad triangles链表,同时依次检查该点所对的边是否为约束边,若为约束边,且该点又侵入该边,将该边加入encroached segment链表;#检查该线段分割后的两个小线段是否被侵入,如果被侵入,也加入encroached seg ment链表尾部;∃重复!、∀、#,直到encroached segment链表为空;7) 重复3)、4)、5)、6),直至bad triangles表为空。
处理坏三角形时采用排序的链表结构,按最小角由小到大的顺序进行,相比不排序的随机处理,或采用先进先出的队列,同样网格质量情况下,最后所得的网格数目小得多,只是排序需花费一定机时。
5 应用实例利用前面所述的算法编制了网格自动剖分程序,以下为该程序的一些剖分实例。
图5为某汽缸体截面的剖分结果,最小角为30 时三角形数量为217个,最小角为35 时三角形数量增加为1044个,约束条件超过35 3 ,程序就不能终止,表明利用本程序该平面区域的三角剖分最小角可达到35 3 ;图6为一只梅花鹿的剖分结果,最小角为34 时三角形数量为1279个,35 时激增为3169个,三角形最小角可达到35 1 。
图5 某汽缸体截面的三角剖分结果图6 梅花鹿的三角剖分结果6 结 论本文提供了一套高质量有限元网格自动生成方法,可对任意复杂平面区域进行剖分,并可对三角形形状和尺寸进行控制,可使三角形最小角达到35 ;127第3期 李梅娥,等:任意平面区域的高质量有限元网格自动剖分采用广度优先搜索查找孔及凹槽中的三角形,不需定义有向边界,可十分方便地定出实际剖分区域;同时给出了一种方便、快捷的点定位法;分割坏三角形时,按角度由小到大的顺序处理,使三角形网格在满足质量要求的同时网格数量较少。
参考文献1 Baker B,Grosse E and Rafferty C S.Nonobtuse triangulation of polygons[J],Discrete Comput Geom,1988,3:147~1682 Bern M,Eppstein D and Gilbert J R.Provably good mesh generation[A],In:Proceedings of th e31st Annual Symposium on Foundationsof Computer S cience[C],IEE E,1990:231~2413 Ruppert J.A Delaunay refinement algorithm for quality2-dimensional mesh generation[J],J Algorithms,1995,18(3):548~5854 胡恩球,陈贤珍,周克定等.有限元网格全自动生成中的初始三角化新方法[J],华中理工大学学报,1996,24(5):19~225 Shew chuk J R.Delaunay refinement algorithms for triangular meshgeneration[J],Comput Geom,2002,22:21~74128应用力学学报第20卷Exact Integration of the linear Element of Laplace EquationIn Boundary Element MethodYuan Zhengqiang 1 H uang Jian 1 Zhu Jialin 2(College of Civil Engineering,Chongqing University,400045,China)1(College of M athematical and Physical S cience,Chongqing Univers i ty,400045,China)2Abstract:The boundary integral in Boundary Element Method affects the precision and the speed of the method.The nonsingular integrals are popularly calculated by the Gauss numerical integ ral,and they are low in precision w hen the source po ints approach the element.The solution by the interpolation of unknow n function is continous,w hich makes the calculation more difficult.This paper presents an alternative way to transform the double integral in Laplace problem on 3 d into the linear integ rals on the boundary of each subdomain,so that all the singular integ rals and nonsingular integrals are calculated by analytical method.It m akes the precision and the speed of BEM im proved.Keywords:boundary element method ,lap lace equations ,ex act integral.Computation of Natural Convection in the Annulus between HorizontalConcentric Cylinders Using Spectral Element methodsChen X uej iang Qin Guoliang Wu L ong X u Zhong(Institute of Fluid M achinery,Xi %an Jiaotong University,Xi %an 710049,C hi na)Abstract:The computation of natural convection in the annulus betw een horizontal concentric cylinders using spectral element method in polar coordinates is presented in this paper.The incompressible Navier Stokes equa tions expressed in terms of primitive variables velocity and pressure and the energ y equation are computed.T he time splitting method for temporal discretization combining w ith the spectral element method for spatial dis cretization,are used.It is show n that the obtained results are in ex cellent ag reement w ith the accepted bench mark solution.Keywords:sp ectral element method,p olar coordinate,natur al convection,time sp litting method.Quality Mesh Generation of Arbitrary Planar Domain for FEML i M ei %e 1 Zhou Jinx iong 2(S chool of M aterial Science and Engi n eering,Xi %an J i aotong University,Xi %an 710049)1(S chool of Civil Engineeri ng and M echanics Xi %an J i aotong University,Xi %an 710049)2Abstract:Quality triangulation is desirable for FEM analysis.In this paper,the Delaunay refinement algorithm is utilized to control the shape of triangles.The data structure of linked list is used to store the bad triangles in dexed by their smallest angle,so that w e can alw ays split the w orst existing triangle firstly.The triangles in holes and concav ities to be removed are found by w idth first search,w hich elim inates the need of oriented11No.3 CH INESE JOURNAL OF APPLIED M ECH ANICSboundaries.An efficient point location algorithm is presented to improve the speed of triangulation.T he pro g ram is used to triangulate several PSLGs(planar straight line graph).The results indicate that the program is capable of triangulating arbitrary complex planar domain and ultimately halts for an angle constraint of up to 350degrees.Keywords:quality tr iangulation,Delaunay r ef inement algor ithm ,ar bitrary p lanar dom ain.Study on Effective Elastic Constants of TubesheetsHomogenization Based TheoryXie Guilan 1,2 Zhang Ping 2 Gong Shuguang 1,2 Chen Yanp ing 2(CAE Institute,School of M echanical Engineering,XiangTan University,411105)1(Institute for computational and applied M athematics,XiangT an U niversity 411105)2(Insti tute for Basic M echanics and M aterial Engineeri ng,XiangTan University 411105)3Abstract :The finite elem ent model based on multiple scale homogenization theory is built for studying on effec tive elastic constant of tube sheet.The numerical result obtained is compared w ith effective elastic constants adopted by the ASME code,and relative error is less than 0.2%.T he expanded tubes reinforcing effect is re searched,the effective elastic constants curves are obtained w hile the rate of pitch of the holes pattern is 0.4.T he numerical result implies that homogenization method is feasible to calculate the effective elastic constants of tubesheets.Keywords:design of tubesheet,eff ective elastic constant ,homogeniz ation theory ,f inite element method.Hydraulic Characteristics of Water Jets by Air Collision in Lateral Direction and Estimation on Limiting Riverbed ScourS un Jian 1 L i Yuz hu 2(Xi %an University of T echnology,Xi %an 710048)1 (T singhua University Beijing 100048,China)2Abstract:By applying the theory of the turbulent jet and the momentum equation to the air collision in lateral direction of two water jets,ex pressions of the velocity vector of the joint jet and the colliding efficiency are pro posed.Both kinds of collision in lateral direction and in up dow n one of w ater jets are analyzed and compared to g ive their differences in colliding plane angle,the colliding efficiency and the air diffusion of the joint jet,and to set up a w ay of estimation on the limiting riverbed scour.The verification through some hydraulic model show s that the forecasted basically tallies with the measured.T hus,it is show n that such an analysis is reasonable and the estimation is feasible.Research results can be used as a reference in desig ning a w ater cushion pool acted by the w ater head from high dam.Keywords:high dam water j ets later al collision p lane angle limiting riverbed scour.12 CH INESE JOURNAL OF APPLIED M ECH ANICS Vol.21。