基于内在特征和能量约束的曲面网格投影新算法

基于内在特征和能量约束的曲面网格投影新算法

有限单元法是在当今技术科学发展和工程分析中应用最广泛的数值方法。在有限单元法之前，网格的三维分类是一项重要的技术。由于工作量大、耗时大、复杂，它已成为有限单元法的主要瓶颈。然而，鉴于三维网格的复杂性，四面体网格的三维模拟具有四面体网格的可能性，因此六面体网格的结构复杂，没有一般算法，但已成为研究对象之一。调查显示,目前60%以上的可六面体划分实体是利用扫掠法生成六面体网格,这充分说明扫掠法是当前运用最为广泛的六面体网格划分算法.扫掠法适用于2.5维实体构成的扫掠体,它由源曲面、目标曲面和一系列的连接曲面组成.对于一个可扫掠的实体,其六面体网格生成算法的一般步骤为:在源曲面上生成四边形网格;在连接曲面上生成结构化的四边形网格;将源曲面网格投影到目标曲面(确保源网格与目标网格拓扑一致);生成内部结点并连接成六面体单元.目前,四边形结构化网格生成算法可以运用映射法和子映射法,非结构化网格生成算法有铺砌法等.文献给出了具体的六面体内节点生成步骤.上述步骤是扫掠法生成高质量六面体网格的关键一步,学者们在此方面做了大量研究.Goodrich将源曲面网格上的每个内节点正交投影到目标曲面上,得到真正位于目标曲面上的节点,该过程被称作“寻根”,但是“寻根”过程计算复杂且耗时.Knupp通过线性变换和光顺的方法实现了源曲面网格投影.Roca等通过源和目标曲面参数域之间的最小二乘近似仿射变换将源曲面上的结点投影到目标曲面.文献的方法借鉴了仿射变换的思想,由于仿射变换保持了模型的“平直性”,即模型间各元素的相对位置不变,这样文献的算法并不适用于源和目标曲面在参数域内形状变化比较大的情况.另外,该方法要求曲面具有参数域,这样就使得源和目标曲面只能是单张曲面.翟建军等利用背景网格上的重心坐标和最小二乘变换将源曲面网格结点投影到各个扫掠层中,该方法仍然不能解决源和目标变形较大的投影问题.Sederberg等提出了多边形的伸缩内在量表示方法,即用边长和相邻边夹角表示多边形,该内在量表示方法具有平移、旋转等刚体运动不变性.Igarashi等将三角形伸缩内在量方法应用到形状针对在变形中,但该方法在形状变化较大且约束过多的情况下会产生网格自交.本研究针对上述算法的不足,结合网格伸缩内在特征方程和网格变形能量,提出了一种保形网格投影算法.1基于连续轴类目标网格的同胚目标网格通过相应的算法,由源网格得到与其拓扑一致的目标网格的过程,称为目标网格的生成,它是扫掠法生成高质量六面体网格的关键一步.文中提出的基于网格内在特征和网格变形能量的保形网格投影算法可以处理多张曲面组成的源和目标面之间的投影,以及变形较大的源面和目标面之间的投影,快速生成保持源网格拓扑一致的目标网格.假设源曲面为S,源网格为SM,目标曲面为T,得到的目标网格为TM,该算法的基本流程如下:(1)将目标曲面T三角化,在此采用Delaunay三角化方法,快速生成保持目标曲面形状的三角网格,称为目标曲面基准网格BM.(2)将目标曲面基准网格BM和源曲面网格SM投影到同一个二维平面上,分别得到平面网格MS和MT,由于本研究是在不固定平面网格边界的情况下得到二维平面网格的,所以本研究采用文献的网格参数化方法.(3)由于MS的边界节点和MT的边界节点是一一对应的,将MS映射到由MT的边界围成的区域,得到网格MS的节点在MT的三角形单元上的重心坐标.(4)利用二维平面网格MS上每个节点在MT三角形上的重心坐标,将MS映射到三维目标曲面基准网格BM上,得到同胚目标网格TM.其中,基于Delaunay的三角网格生成算法以及三维网格参数化算法都是当前研究的热门问题,很多学者也提出了很多成熟的算法,在此不多做阐述.接下来将详细阐述步骤(3)中的网格映射算法.2网格结构的自交性以及网格特征的实现网格之间的映射是目标曲面网格生成过程中的关键一步,很多学者对此做了研究.Floater提出的凸组合算法解决了外边界为凸边的网格之间的映射,但是不能解决外边界为凹并且网格形状变化比较大的情况.本研究结合平面网格的内在特征和有限元变形能量,使映射后的网格保持了原有的拓扑结构,又尽可能保持了之前网格的特征,并避免了网格的自交.网格M可定义为一个二元组M=(V,K),V是顶点集合,V={v0,v1,…,vn},n表示网格顶点的个数;K是单纯复形,它表示M的顶点、边和面的拓扑关系,它的子集称为单纯形,可分为顶点、边、面3类:顶点v={i}ue84aSymbolNCpK,边e={i,j}ue84aSymbolNCpK,面f={i,j,k}ue84aSymbolNCpK.对于顶点v,它的一阶邻域是与v拓扑相连的所有顶点的集合,记作N(v)={j{i,j}ue84aSymbolNCpK},顶点的度数N(v)定义为所有与v相邻的顶点个数.2.1维几何特征方程的构建结合文献中对于任意多边形的均值坐标定义,本研究提出平面网格的伸缩内在量表示方法,该表示方法可以表示三角形和四边形的平面混合网格.对于网格M中顶点v0及其一阶邻域N(v0)构成的多边形,如图1所示.式(1)共有|N(v0)|个方程,其中R(α)是旋转矩阵:R(α)=[cosα-sinαsinαcosα](2)由此可以知道,对于顶点及其一阶邻域点构成的多边形,式(1)是与平移、缩放和旋转无关的方程组,它决定了顶点及其一阶领域点之间的相互几何位置.在二维平面上,顶点v0坐标记为v0=(x,y),则式(1)中的向量方程可以表示为方程组(1)由N(v0)个式(3)组成,空间网格中的每个顶点可得到一个邻域多边形特征方程组,联立可得如下的网格伸缩内在量特征方程组:AX=0(4)式中:X为待求的平面网格顶点坐标向量,X=(x0,y0,x1,y1,…,xn,yn)T;系数矩阵A为2m×2n阶稀疏矩阵,其中m为所有顶点的度数之和,m=∑vi■SymbolΝCpΜN(vi).2.2单元应变矩阵为更合理地表示平面网格的映射,本研究引入了平面网格单元的有限元变形能.为简化计算,本研究将四边形单元分割成两个三角形单元,所以只需计算三节点单元的变形方程.由文献可知,假设典型三节点单元节点整体编码为f,g,h,以逆时针方向编码为正向,节点坐标分别为(xf,yf),(xg,yg),(xh,yh).平面三角形单元应变矩阵为B=[BfBgBh]=12A[bf0bg0bh00cf0cg0chcfbfcgbgchbh](5)式中:A是三角形的面积;bf=yg-yh,cf=xh-xg,bg=yh-yf,cg=xf-xh,bh=yf-yg,ch=xg-xf.平面三角形单元的刚度矩阵为Κe=BΤDBA=[ΚffΚfgΚfhΚgfΚggΚghΚhfΚhgΚhh](6)式中:矩阵D是文献中给出的弹性矩阵,矩阵Ke中的任一分块矩阵可表示为Κrs=BΤrDBsA=E04(1-ν2)A[e1e3e2e4](7)(r,s=f,g,h)式中:E0、ν分别是材料的弹性模量和泊松比;e1=brbs+1-ν2crcs‚e2=νcrbs+1-ν2brcs‚e3=νbrcs+1-ν2⋅crbs,e4=crcs+1-ν2brbs.由单个三角形单元的刚度矩阵Ke组合可得到三角网格的有限元变形能量方程:KΔX=0(8)式中:K为由矩阵Ke组合而成的刚度矩阵,ΔX=(Δx0,Δy0,Δx1,Δy1,…,Δxn,Δyn)T,是待求的平面网格节点的偏移量.2.3约束矩阵及约束矩阵假设MS的边界节点Pu(u=0,1,…,l-1)和MT的边界节点Pt(t=0,1,…,l-1)一一对应,为了使映射后的MS保持原有特征,而又不产生自交等改变网格拓扑的情况.本研究采用罚函数法求解映射后的平面网格.在文中提出的算法中,通过让网格的内在特征变化最小,同时又满足有限元变形能的约束来建立罚方程:G(X)=XΤAΤAX+θi2∑i=1C2i(9)式中:θi表示给定的罚因子,一般取值为105;Ci是约束方程;A由式(4)求得.由于边界节点Pu的坐标等于Pt的坐标,因此可得到2×l个约束条件,可以得到约束矩阵:C1=JAX-V(10)式中:JA为2l×2n阶稀疏矩阵;V为2l×1阶矩阵,由MT的边界节点Pt的坐标组成.同样,由式(8)可以得到网格的有限元变形能约束矩阵:C2=KΔX(11)式中:ΔX=X-X0,X0是网格节点初始坐标向量,给定约束方程的罚因子θ1和θ2,代入式(9)得G(X)=XTATAX+θ1(JAX-V)T(JAX-V)+θ2(X-X0)TKTK(X-X0)=XTAKX+θ1(VTV-2XTJTAV)+θ2(XT0X0-2XTKTKX0).这里,系数矩阵AK=ATA+θ1JTAJA+θ2KTK,令∂G(X)∂X=0,因此,有AKX=θ1JTAV+θ2KTKX0(12)得到X=A-1K(θ1JTAV+θ2KTKX0).因为矩阵A、JA和K都是稀疏矩阵,所以方程组(12)的系数矩阵AK也是稀疏矩阵,该稀疏方程组可以采用TAUCS库快速求解.求得网格MS映射后的节点坐标后,MS的边界节点Pu越来越接近MT的边界节点Pt,计算Pu和Pt每个对应点的距离,如果对应点的距离大于给定的容差ε,则利用映射后的MS继续建立罚方程,求解映射后的网格坐标,直到Pu和Pt之间的距离小于给定的ε.图2和图3所示为平面网格映射实例.其中图2所示实例中网格外边界保持不动,内孔位置变化较大,由图2(b)可看出,映射后的网格不会产生自交,能保持原始平面网格的网格形状.图3的平面网格外部轮廓变化较大,应用文中提出的方法仍然可以得到较好的映射结果.3源网格投影到目标曲面为验证文中提出的方法,本研究使用VisualStudio2008和三维实体几何引擎OpenCascade实现上述算法,实验用的微机配置为Pentium(R)4CPU3.00GHz,2.00GB内存,并已成功应用到具有自主知识产权的三维实体网格划分软件3D-Mesher系统中.图4反映的是源网格投影到目标曲面上的部分中间过程,其中源曲面由两张裁剪曲面组成.首先将源网格投影到平面上,同时目标基准网格也投影到同一平面上;然后将平面源网格映射到平面目标基准网格上;最终根据平面映射结果得到目标网格、以及模型的六面体网格划分结果,图中只标示了3层六面体网格.图5所示为文中提出的算法在复杂扫掠体六面体网格生成中的应用,图中所示实例源曲面和目标曲面形状都很复杂,曲面边界都是凹多边形且变化剧烈,利用文献的方法和文献的方法得到的目标网格都会产生自交.图5(a)所示是模型的六面体网格划分结果,总共有36层网格,图5(b)-(d)分别是模型的第20层、第30层和第36层六面体划分结果.图6为应用文中提出的网格投影算法得到的典型零件六面体网格划分结果,图6(a)是螺母模型的六面体网格划分结果,模型的源曲面和目标曲面都包含有多张曲面.图6(b)是套筒模型的六面体网格划分结果.从上述实例中可以看出,文中提出的算法可以实现任意复杂扫掠体的六面体网格划分,并得到质量较高的六面体网格.4多张源曲面并进行目标网格生成的向量程式本研究针对已有曲面网格投影算法的不足,提出了一种新的保形网格投影算法,并已成功应用于实体网格生成软件中.大量实践表明,该算法稳定可靠;与基于最小二乘近似仿射变换的网格投影算法不同,该算法不但可以实现多张源曲面或多张目标