基于曲率球的参数曲面自适应三角化方法研究

1、 PAGE 9基于曲率球的参数曲面自适应三角化方法研究杨伟民1,2，楼天良1，林志伟2，沈洪垚21.义乌工商职业技术学院机电信息分院，义乌，322000 2.浙江大学现代制造工程研究所，杭州，310027摘要：参数曲面的三角化算法是有限元分析和数控机床路径规划的基础。文中提出了一种新的基于分层递推思想的三角化算法，直接在曲面上生成三角网格。根据曲面微分几何定义相对保守的曲率球概念，将局部曲面近似成曲率球面，在球面上建立三角形边和面片的误差模型和递推关系，计算三角形各个顶点在欧式空间以及参数空间的坐标。然后采用逐层递推的方法，从曲线边界开始，逐层生成三角形网格，直至满足收敛准则。文章最后将此算法

2、用C+程序语言编程并通过NURBS曲面进行实例验证，结果表明该算法能够较好的对参数曲面进行三角网格划分。关键词：参数曲面；三角化；曲率球；误差控制中图分类号：TP391.73Parametric surface triangulation based on curvature sphereWeimin Yang1,2，Tianliang Lou1，Zhiwei Lin2，Hongyao Shen21.Yiwu Industrial & Commercial College，Yiwu，322000 2.Zhejiang University，Hangzhou，310027Abstract：Par

3、ametric surface triangulation algorithm is the basis of both the finite element method and the tool path generation in CNC machining. In this paper, a new method is proposed to directly generate triangular meshes on the surface through a layered recursive method. The concept of curvature sphere is d

4、efined. The local surface is approximated to be the spherical surface of the curvature sphere. On the curvature surface, the error models of the triangular mesh and the original surface are established and the coordinates of the triangular vertexes, both in the parametric domain and the Euclidean sp

5、ace, are calculated. Then the layered recursive method is applied to generate the triangular meshes layer by layer from the boundary of the surface until the convergent condition is met. The proposed method is implemented in C+ programing language and is tested by a NURBS surface example. Results sh

6、ow that the proposed method can deal well with the parametric surface triangulation problem.Keywords：Parametric surface; triangulation; curvature sphere; error control0 引言三角化是指对已知曲面方程的参数曲面进行三角形离散划分。它是有限元分析和数控机床路径规划的基础。三角化算法的优劣直接关系到这些计算的计算量和准确度。首先要求保证三角化离散的精度，即保证自由曲面离散后的三角面片能够逼近原始曲面；其次要求离散具有自适应性，即生成的

7、三角面片的疏密程度能够随着曲面曲率变化而变化，在曲面曲率大的地方生成的三角面片致密一些，而在曲面曲率小的地方生成的三角面片相对稀疏一些。由于路径规划中干涉检测算法需要遍历曲面中所有三角面片，所以在保证精度的前提下，曲面三角面片数量越少越好。对参数曲面而言，最简单、最直观的三角化方法是映射法。这种方法先在参数空间实现平面三角化，如Delaunay平面三角剖分 ADDIN EN.CITE ADDIN EN.CITE.DATA HYPERLINK l _ENREF_1 o Cignoni, 1998 #252 1, HYPERLINK l _ENREF_2 o 余杰, 2010 #9 2，然后将得到

8、的三角形映射到欧式空间。这种方法由于参数空间和欧式空间的映射没有空间形态的必然联系，得到的三角面片没有自适应性，难以评价三角化精度 ADDIN EN.CITE 王青2003131117王青程艳旗林海鲍虎军浙江大学国家重点实验室,浙江大学国家重点实验室,浙江大学国家重点实验室,浙江大学国家重点实验室 浙江杭州310027 ,浙江杭州310027 ,浙江杭州310027 ,浙江杭州310027隐函数曲面的实时采样与三角化技术浙江大学学报(理学版)浙江大学学报(理学版)621-62706隐函数曲面多边形化曲面离散20031008-949733-1246/NCnki HYPERLINK l _ENRE

9、F_3 o 王青, 2003 #1 3。徐松等在传统映射法的基础上，在平面三角化算法中环边统一处理，采用自适应三角网格加密法，有效地处理带有特征约束条件的任意曲面三角剖分问题 ADDIN EN.CITE 徐松2000242217徐松王剑英中国科学技术大学计算机系!合肥230027,中国科学技术大学计算机系!合肥230027曲面的自适应三角网格剖分计算机辅助设计与图形学学报计算机辅助设计与图形学学报267-27104三角剖分特征约束自适应网格加密曲面200011-2925/TP11-2925/TPCnki HYPERLINK l _ENREF_4 o 徐松, 2000 #2 4。关振群等提出一种

10、基于映射法的复杂三维组合曲面的有限元网格全自动生成方法，通过引入虚边界解决了闭合曲面在参数域中边界不完整的问题，并通过调节虚边界提高了复杂组合曲面网格生成的质量 ADDIN EN.CITE 关振群2003353317关振群隋晓峰顾元宪李云鹏大连理工大学工业装备结构分析国家重点实验室,大连理工大学工业装备结构分析国家重点实验室,大连理工大学工业装备结构分析国家重点实验室,大连理工大学工业装备结构分析国家重点实验室 大连116024 ,大连116024 ,大连116024 ,大连116024复杂三维组合曲面的有限元网格生成方法计算力学学报计算力学学报409-41604网格剖分闭合曲面虚边界组合曲面

11、20031007-470821-1373/O3Cnki HYPERLINK l _ENREF_5 o 关振群, 2003 #3 5。吴宝海等采用改进波前法对曲面进行三角化细分，从曲面两条边界向曲面中心同时推进，避免了常规波前法中由于曲面角点的不良形态导致网格规划的失败以及生成低质量的网格 ADDIN EN.CITE 吴宝海2005464417吴宝海王尚锦西安交通大学能源与动力工程学院,西安交通大学能源与动力工程学院 西安710049 ,西安710049参数曲面网格生成的改进波前法计算机辅助设计与图形学学报计算机辅助设计与图形学学报1686-169008参数曲面三角剖分波前法20051003-9

12、77511-2925/TPCnki HYPERLINK l _ENREF_6 o 吴宝海, 2005 #4 6。梁义等结合映射法和前沿推进技术，通过曲率计算自动识别曲率特征并创建密度场，对组合参数曲面生成几何自适应的网格 ADDIN EN.CITE 梁义2010676617梁义,陈建军,陈立岗,郑耀,几何自适应参数曲面网格生成计算机辅助设计与图形学学报计算机辅助设计与图形学学报327-33520102曲面网格生成 自适应 Delaunay三角化 黎曼度量 中轴20101003-977597390XVIP HYPERLINK l _ENREF_7 o 梁义, 2010 #6 7。李李等利用均匀采

13、样点来估算曲面曲率，针对裁剪曲面进行三角剖分，提高剖分效率和精度 ADDIN EN.CITE 李李2006585517李李 王亚平 裁剪曲面自适应三角化剖分计算机应用计算机应用2006年S1期TP391.41针对CAD/CAM领域中裁剪曲面的三角剖分问题,从提高效率和满足剖分精度的角度出发,提出了一种较为实用的自适应离散方法,与传统的自适应方法相比,三角片数量有所减少,速度有明显的提高。2006 HYPERLINK l _ENREF_8 o 李李, 2006 #5 8。王飞等基于B样条曲面的适应性分割，将背齐尔曲面片的三角化处理扩展到整张曲面，消除相邻面间边界处可能出现的裂缝 ADDIN EN

14、.CITE 王飞1997797717王飞,董国辉,基于b样条曲面适应性分割的三角化算法北京邮电大学学报北京邮电大学学报95-98202计算机图形学 B样条曲面 三角化算法19971007-532191520AVIP HYPERLINK l _ENREF_9 o 王飞, 1997 #7 9。Shimada等则提出用BMG（Bubbul Mesh Generation）对曲面进行三角化 ADDIN EN.CITE Shimada19988108817Shimada, KenjiGossard, David C.Automatic triangular mesh generation of trim

15、med parametric surfaces for finite element analysisComputer Aided Geometric DesignComputer Aided Geometric Design199-222153TriangulationDelaunay triangulationVoronoi polygonsSphere packingPhysically based modelBubble meshingDynamic simulationFinite element analysis19980167-8396/science/article/pii/S

16、016783969700037X/10.1016/S0167-8396(97)00037-X HYPERLINK l _ENREF_10 o Shimada, 1998 #8 10。在本文中，我们结合曲面的局部微分几何，提出曲率球概念，在局部范围内，用曲率球近似曲面，在曲率圆上直接近似生成三角面片，并利用分层递推思想，将生成的的三角片面覆盖至整张曲面。1 参数曲面的微分几何对于给定的参数曲面，只要知道参数域内的任意一点坐标值，就可以根据相应的公式计算出曲面在欧式空间内对应的点，以及曲面在该点的一阶、二阶偏导、和单位法向量。曲面在点处的法曲率可以表示为： （1）其中：，为曲面第一基本型；，为曲面

17、第二基本型。根据式(1)可以得到曲面在该点的主曲率： （2）其中：；。图1 曲面的微观几何特性不妨设，根据、在数轴上取值正负不同，可以把曲面的微观几何特性分成5种情况，如图1所示。图1b所示为椭圆抛物面，对应图1a中的第1和第5两种情况，即、同时为正或同时为负。当、同时为正时，椭圆抛物面开口朝上；当、同时为负时，椭圆抛物面开口朝下。图1c所示为抛物面，对应图1a中的第2和第4两种情况，即、两者中其中有一个为零。当为零时，抛物面开口朝下；当为零时，抛物面开口朝下。图1d所示为猴鞍面，对应图1a中的第3中情况，即为正，为负。为了简化讨论和计算，引入曲率球的概念。曲率球的定义如下：设、为曲面上处的两

18、个主曲率，取，把球心在过的法线上，半径为的球称为曲面在点的曲率球。由曲率球的定义可知，曲率球是以曲面上某点的极大曲率为半径建立的球面，因此在将曲面附近的点转移到曲率球上考虑，从误差角度上讲，是偏安全的。2 单个三角面片生成及误差控制上节定义了曲率球的概念，本文涉及到的误差都是将曲面局部近似到曲率球面上分析的。本节详述在给定的精度条件下曲面上单个三角面片的生成。2.1 曲线上弦生成在考虑在曲面上获得一条符合精度要求的三角形边之前，先考虑在曲线上获得一条符合精度要求的弦。对空间参数曲线及其曲线上的一个已知点，在曲线上点附近求点，使得弦和原始曲线段之间的弓高误差为不大于给定的弦容差。近似地，本文在点

19、的曲率圆上考虑问题。如图2所示，设曲线上弦长度为，近似用圆上的弦代替。在圆中，为中点，连接交圆于点，则线段就是弦的弓高，也即原始曲线上弦的弓高。设曲线上点曲率圆半径为，在中，易得：图2 在曲线上得到一条符合精度要求的弦 （3）对式（3）忽略2次小量，得：。由上可知，只要保证，就能保证弦的精度。临界情况下，弦的长度可以表示为： （4）下面根据式（4）中得到的值反向近似求解点参数值。设一质点沿着曲线运动，在某点速率为，方向为曲线切线方向。曲线的参数是时间的变量，写成，随着时间递增，因此。由于： （5）所以： （6）令，它们分别代表，点参数值。对在处一阶泰勒展开并将（6）式代入得： （7）为质点走过

20、的路程，近似认为该路程等于弦长度，即：。将式（4）代入（7）式，得： （8）根据式（8）即可求出点的参数值，在曲线上的位置可完全确定。运用上述方法可以求得任意曲线上一条符合精度要求的弦。2.2 曲面上弦生成下面考虑在曲面上获得一条符合精度要求的弦，即三角形的一条边。设是曲面上一点，其参数坐标为；球为点处的曲率球；是曲面上领域内的任意一点。当趋近，即时，近似认为点落在球的球面上；那么曲面上弦的弓高就是弦的两个端点都落在球上时的最大弓高，即在球最大圆上的弓高。图3 在曲面上得到一条符合精度要求的弦如图3所示，考虑过点的、向曲线，由式（8）可在、向曲线上分别找到符合弓高误差的弦和弦： （9） （10

21、）式（9）、（10）中为曲率球的半径，为弦、容许误差，、分别为曲面对、的一阶偏导。式（9）、（10）中用曲率球半径求出的、值是偏安全的。如图3所示，分析点的误差，点坐标可以表示为。令为，的夹角，记弦长度为，、长度分别为、，则由图3可知： （11）式（11）中， 长度随着变化而变化。显然，在曲率球面上，当最长时，弓高误差最大。因此在最坏的情况下，即时，此时。而由式（4）可知，因此可得曲面上弦的最大弓高误差为： （12）要保证弦的弓高误差小于给定弦容差，即，只需保证： （13）即在计算式（9）、（10）中的、时的容许误差要小于给定弦容差的四分之一。2.3 曲面上三角形生成本文中提出的生成三角形面片

22、的算法是递推的，因此认为的一条边的两个端点在三维空间的坐标以及在参数域内值已经在前一步中计算出来。这里要考虑的问题是如何根据这条已知边在曲面上找到第三个点，使得点分别和，连起来之后得到的尽量能够符合精度要求。 如图4所示，用2.1中的方法，对点进行一次运算，在参数域内得到点坐标；同理对点进行运算，在参数域内得到点坐标，这里用是让尽量靠近，从而提高精度。取中点，其在参数域内坐标为：由此得到的就是所要求的三角形，如图4。下面对进行误差分析。在微观情况下，由于三边是在同一点邻域内产生的弦，在这点邻域内曲面的曲率变化很小，近似认为这三边所在的曲率球的半径相同。由于三边弓高一样，所以的三条边长相同，认为

23、是等边三角形。图4 在曲面上得到一个符合精度要求的三角形的三个顶点都落在曲面上，近似情况下以三角形其中一个顶点的曲率球代替附近曲面。如图5所示，垂直于所在平面交球面于，则所在平面到曲面的弓高就近似为长度。图5三角面片在曲率球上的近似弓高图6 三角形面片弓高和三角形边弓高的关系如图6所示，圆弧为球大圆上的圆弧，球半径为，即为在2.2中弦最大弓高误差。则在直角中： （14）在直角中，： （15）在直角中： （16）令，则： （17）忽略平方项，得： （18）所以要想把的弓高误差限定在，只要在2.2中保证： （19）综合式（13）、（19），若要保证三角面片生成精度，只需在计算式（9）、（10）中的

24、、时的容许误差要小于给定弦容差的四分之一。3 三角网格生成第2节中介绍了曲面上单个三角面片的生成以及误差控制方法，下面阐述曲面三角网格生成。三角网格生成时遵循“点对点，边对边”原则，即每个三角形的顶点必须对应另外三角形的顶点，每个三角形的边对应另外三角形的边；不允许出现某个三角形的顶点对应另外三角形的边，否则生成的三角网格会出现裂缝。本文提出一种分层递推生成三角面片的方法，直接在参数曲面上生成三角网格，生成的三角网格无裂缝而且符合精度要求。对一张参数曲面，其中，如图7所示，其边界线分别是线，线，线，线。按照本文提出的分层递推方法对曲面逐层进行三角网格划分，假设第层划分结束后形成的的图形如图7所

25、示，其中左边有三角网格的是已经进行三角划分过的曲面，右边的是没有进行三角划分的曲面，它们之间的分界折线是由当前一层三角形的边首尾依次连接得到的，记这一层折线为第层折线。下面对第层折线进行轮廓更行，同时在曲面上生成三角形，得到第层三角形面片。图7 第k层三角面片生成图8 折线段轮廓更新分析第层折线段、，首先考虑，夹角。当时，则直接连接得到，同时删除、，成为新的，成为，下面的线段序号依次减1。下一个分析对象为。当时，则对线段和夹角考虑。当时，则直接连接得到，同时删除、，成为新的，下面的线段序号都不变。下一个分析对象为。当，那么对根据两端点，按照2.3节中方法计算出曲面上一点，得到的，同时删除，成为

26、新的，成为新的，成为新的。下一个分析对象是。图8所示的情况由于得到的新三角形面片都在折线右侧，即在曲面未三角化区域，此时不存在问题。但当得到的新三角形面片在折线左侧，也就是在曲面已经三角化的区域，此时新三角形会和原来的三角形重叠起来，这种情况应该避免。事实上，图8中的夹角是根据向量，得到的： （20）式（8）无法计算当线段和的夹角大于180时的情况。而新三角形出现在左边正是因为没有检测到和夹角大于180的情况。为此引入凹角以及凸角，直观的，相对于右边还没有三角化了的区域，当角是图8中的情况时，认为是凹角，它是凹向未三角化区域的；当是图9的情况时，认为是凸角，它是凸向未三角化区域的。图9 凸角定

27、义以上根据观察给出了凹凸角的简单定义，下面给出凹凸角精确的数学定义。点位于曲面上，如图8所示，点法向量可由得到。令向量，理论上当生成的三角面片无穷小的时候，都是位于的切平面上的，此时有和方向完全一致。而实际上因为生成的三角形不可能无穷小，和方向不可能真正完全一致。人为规定一个，向量方向一致的判据。设向量和向量的夹角，人为规定当时，和同向；当时，和反向。而当和同向时，是凸角；当和反向时，是凹角。所以在前述的算法中，加上对的判断，如果是凸角，则跳出直接去判断和的夹角是凹角还是凸角。如果是凹角，则进入判断的情况。3.1 边界处理以上讨论中没有涉及到曲面边界，下面针对边界进行讨论。对第层折线的第一条线

28、段，它的起点所在的曲面边界可能有3种情况。如图10所示，分别是在线上，在线上，在线上。这三种情况可以简单地从点的参数坐标判断出来。分析在线上这种情况，其他也是类似的。在图10中，为线和线在的切线的夹角。当时，按照2.1节中方法在线上得到一段符合精度要求的弦,得到的三角形作为一个新的三角形，同时删除，将变为，其他都保持不变。而当时，则在考虑和的所成角的凹凸情况。如果是凸角则直接根据两端点、按2.3节中方法形成新三角形，同时处理；如果和的所成角是凹角，计算和的角度：如果该角度大于90则按照2.3方法根据、形成新三角形，同时处理，否则连接、以及端点得到新的三角形，同时处理。对层折线段的最后一条线处理

29、和第一条线是类似的。也分为三种情况，的终点分别和线，线，线相交；处理也和类似，都是遵循先和前一条线段求角，然后再考虑分别和曲面界面的切线求角，这里不再赘述。图10 网格生成边界处理 图11 网格生成的收敛条件3.2 网格的起始和收敛递推法对曲面进行三角化处理依赖于第一层折线段的生成。本文中第一层折线段取自曲面边界，用2.1节中方法在边界曲线上连续取出一系列折线段。下面分析网格生成的收敛条件。可以想象当网格生成快要结束的时候，对应的那一层折线段会碰到曲面的边界，如图11所示。每一层折线段的起点和终点都落在曲面的边界上，在每一层折线段对每条线段生成三角形之前，先对这条的终点进行判断，如果的终点在边

30、界上，那么的起点也必然在曲面边界上，跳过和，下一个要分析的是；如果的终点不在边界上，则照常处理。对这样一层处理结束之后对新得到逐个扫描，判断终点是否落在边界上，并记录终点落在边界上的个数，记为。当满足： （21）则可以判断曲面三角网格生成结束。最后再针对生成最后一层三角形，结束网格生成。4 实例演示本文提出的算法用C+语言编程，以下用一张NURBS曲面验证本文算法。该曲面的控制顶点信息如下：节点矢量分别为：权值系数为： 。三角化结果如图12所示。由图12可知，本文算法能够较好地对参数曲面进行三角化，在曲面曲率小的地方，能够生成较疏的三角网格；而在曲率大的地方，生成相对密的三角网格。图12 NURBS曲面三角化实例5 总结三角化算法是有限元分析和数控机床路径规划的基础。本文提出了一种新的针对自由曲面三角化的算法。定义了曲率球的概念，将复杂的曲面误差模型转化到曲率球上考虑，大大简化了误差计算。先后在曲面上获得符合精度要求的弦和三角形面片。在此基础上，利用递推算法，从曲面边界开始，逐层递推生成三角网格，并确定了三角网格生成的收敛条件。本文算法经过实例验证，表明算法确实有效。参考文献 ADDIN EN.REFLIST 1.Cignoni P, Montani C, Scopigno R (1998) DeWall: A fast divide and c