基于扫掠法的六面体网格生成算法

基于扫掠法的六面体网格生成算法

几何分解法随着计算机技术的发展，区块技术被广泛应用于各种工程的实际问题。作为有限元前处理中非常重要的网格生成也成为了一个十分重要的研究方向。目前在实体网格生成算法方面,四面体网格生成算法已经十分成熟,但由于六面体网格的单元个数少、可精确反映几何模型的结构和更接近分析的要求等原因,人们更倾向于生成六面体网格。扫掠法(Sweeping)是一种通过把曲面网格拉伸为实体网格从而将2.5维实体生成六面体网格的方法。研究表明,在加入了几何分解(Geometrydecomposition)的前提下,有50%以上的三维六面体网格都用扫掠算法生成。过去10年里扫掠算法受到人们广泛的研究:如Blacker的“Thecoopertool”采用了最小二乘平均的技术来解决变扫掠截面和扫掠方向的问题;Staten等人的“BMSweep”方法利用了背景网格上的重心坐标来将源结点投影到各个扫掠层中;Roca等通过源和目标曲面参数域之间的最小二乘近似仿射变换来将源曲面上的结点投影到目标曲面,并进一步将此算法扩展到三维空间来得到内结点;Knupp则分别通过线性变换和光顺(Smoothing)的方法实现了扫掠。本文基于Roca与Staten的方法实现了一种新的算法,该算法利用BMSweep算法的内结点生成解决了Roca算法中的重复解线性方程组而造成的低效率,又利用了Roca算法中的网格映射方法来避免源网格映射到目标网格时的“寻根(Rootfinding)”问题。1基于扫描法的体格设计算法1.1扫掠体的六面体网格扫掠算法适用于2.5维实体(沿扫掠轴的剖面拓扑结构为常量的实体)构成的“一到一”(One-to-one)的扫掠体的六面体网格生成。具体约束如下:(1)扫掠体由一个源曲面、一个目标曲面和一系列的连接曲面组成,如图1所示;(2)源和目标曲面必须拓扑结构相同(有相同数目的孔洞和逻辑边);(3)源和目标曲面的参数化存在;(4)连接曲面必须由4个逻辑边定义(可以用映射法划分结构化网格);(5)扫掠体有且只有一条扫掠轴。1.2建立结构化网格,确保封扫掠法生成六面体网格算法的一般步骤如下:(1)在源曲面上生成四边形网格,详见2.1节;(2)在连接曲面上生成结构化的四边形网格,详见2.2节;(3)将源曲面网格投影到目标曲面(确保源网格与目标网格具有一致的连接性),详见2.3节;(4)生成内部结点并连接成单元,详见第3节。2网格结构2.1在其他网格生成器中的应用在本算法中,源曲面网格是应用已有的基于Paving算法的四边形网格生成器得到,该网格生成器已经应用于ETA公司的板料成型分析软件DYNAFORM中。2.2边形网格生成算法因为连接曲面由4个逻辑边定义(每个逻辑边可以由多个边组成),所以可以很方便地利用四边形网格生成算法生成结构化的网格,这里使用的是超限映射法(Transfinitemapping)。已有基于此方法的网格生成器,该算法不在本文讨论范围内。需要注意的是,连接曲面网格与源曲面网格在边界处必须保持一致(Conformed)。2.3源曲面网格ms的建立目标曲面网格要满足与源曲面网格具有相同的拓扑,一般都使用映射的方法将源曲面的网格映射到目标网格,问题的关键就是找到这个映射˜φ。一般的源与目标曲面间的网格映射都是在三维空间内直接展开的,即有S⊂R3˜φ→Τ⊂R3(1)式中:S为源曲面;T为目标曲面;˜φ为S到T的映射。这种算法必须进行“寻根”运算来找到映射后的结点在曲面上的投影,因而增加了计算量,降低了效率。为了避免这一问题,本文采用参数域之间的映射φ来间接地实现˜φ,即有S⊂R3˜φ→Τ⊂R3φS|DS↑↑φΤ|DΤDS⊂R2φ→DΤ⊂R2(2)式中:DS,DT分别为S与T的参数域,φS|DS,φT|DT分别为由DS到S的映射,DT到T的映射。由式(2)可知,˜φ可以通过式(3)得到˜φ=φΤ|DΤ⋅φ⋅φS|DS-1(3)其中:φS|DS-1为DS到S的逆映射,但却不需要用解析的方法求得,只要在生成源曲面网格MS时既存储网格的结点坐标,又存储其参数信息即可。而映射φT|DT为由曲面上点的参数信息得到该点坐标,每一个曲面造型系统都会提供此功能。故只要能解决从DS到DT之间的映射φ,˜φ也就能得到。令n(n≥3)为源/目标曲面边界上的结点个数,令US={uiS}i=1,…,n⊂R2,UT={uiΤ}i=1,…,n⊂R2分别为源和目标网格边界结点的参数坐标,只要能找到一个函数φ使得φ(uiS)=uiΤ,i=1,⋯,n(4)可以用一个仿射变换来近似得到φuΤ=φ(uS)≈Α(uS-uarbS)+b(5)式中:A为一个原点在uarbS的线性变换,而b为一个平移向量。但是给定源和目标的边界结点,一般找不到相应的仿射变换,所以通过最小二乘拟合来寻找一个线性变换A和一个平移向量b。这样的A与b应该满足使式(6)取最小值。F(A,b)=n∑i=1∥uiΤ-(A(uiS-uarbS)+b)∥2(6)如果取uarbS=ucS=1nn∑i=1uiS,那么b=ucΤ=1nn∑i=1uiΤ(7)因此作者定义下述坐标系ˉuS=uS-ucS,ˉuΤ=uΤ-ucΤ(8)在此坐标系下,式(5)可以写成ˉuΤ=ˉφ(ˉuS)≈AˉuS(9)式中:ˉφ为φ在新坐标系下的形式。此时式(6)变为F(A)=n∑i=1∥ˉuiΤ-A(ˉuiS)∥2(10)解此最小二乘问题得A0,使得下述仿射变换成立。φ0(uS)=A0(uS-ucS)+ucΤ(11)到此已经得到了一个由DS到DT之间的映射φ0,根据式(3),定义˜φ0(p)=φΤ|DΤ(φ0(φS|DS-1(p))),p∈S(12)式中p为源曲面上任意一点。因此由源曲面网格MS到目标曲面网格MT的函数为ΜΤ=˜φ0(ΜS)(13)3曲面网格的连接经过上述3个步骤,已经得到了同拓扑的源/目标曲面网格和连接面上结构化的四边形网格。四边形网格上沿扫掠方向的结点连接成众多的“肋线”,如图2所示。有了这些条件,就可以计算出内部结点。3.1角网格形成背景网格提供了用于计算内层结点位置的框架。背景网格是通过将源网格的边界结点生成三角网格而形成的。由于不要求背景网格的质量,但要求更快速、更健壮的生成算法,所以采用二维的Bowyer/WatsonDelaunay三角化算法在源曲面参数域生成平面网格,再映射到三维空间得到要求的背景网格。3.2计算插值信息每个源网格内结点的插值信息包括:包含该结点的背景网格三角形;该结点在该三角形中的重心坐标;该结点相对于源/目标网格的偏移距离。3.2.1可建立三角形通过每个结点NiS的参数值在源曲面的参数域背景网格中可以很方便地找到该结点所在的三角形Ti。这就建立了结点与包含它的三角形间的映射关系。显然,一个三角形可能包含多个结点,所以该映射为一对多的映射。3.2.2重心坐标的计算目标曲面内结点NiS在背景网格三角形中的重心坐标用于计算体网格内结点的位置。对于每个NiS需要计算出两组重心坐标值:(1){aiS,biS,ciS}:NiS相对于三角形Ti在源曲面参数域中的重心坐标。(2){aiΤ,biΤ,ciΤ}:NiΤ相对于三角形Ti在目标曲面参数域中的重心坐标。3.2.3格内终点nis的位置偏移距离是一个源/目标网格内结点到包含它的背景网格三角形Ti的空间距离。如果源和目标曲面不平行,那么就需要用偏移距离来调整内结点的位置。令PiS为源曲面网格内结点NiS的位置,(PiS1,PiS2,PiS3)为该结点的背景网格三角形Ti在源曲面网格上的3个顶点的位置;PiΤ为目标曲面网格内结点NiΤ的位置,(PiΤ1,PiΤ2,PiΤ3)为该结点的背景网格三角形Ti在目标曲面网格上的3个顶点的位置。则diS=∥(aiS⋅ΡiS1+biS⋅ΡiS2+ciS⋅ΡiS3)-ΡiS∥(14)diΤ=∥(aiΤ⋅ΡiΤ1+biΤ⋅ΡiΤ2+ciΤ⋅ΡiΤ3)-ΡiΤ∥(15)式中:diS,diΤ分别为源和目标网格的偏移距离。3.3构建背景网格三角形ti一旦上述的插值信息都计算完毕,即可得到第L层的内结点NiL的位置PiL。计算过程分为4步:(1)由源和目标曲面的信息线性插值出重心坐标{aiL,biL,ciL}和偏移值diL,插值参数可以由弦长参数化方法确定。(2)在当前层构建背景网格三角形Ti。可以通过取得源曲面背景网格三角形顶点所在的肋线,在当前层L的结点构成三角形(PiL1,PiL2,PiL3)。(3)计算偏移方向ViL。偏移方向可用三角形(PL1i,PL2i,PiL3)顶点处的扫掠方向单位向量{VL1i,VL2i,VL3i}线性插值得到VLi=aLi⋅VL1i+bLi⋅VL2i+cLi⋅VL3i(16)(4)计算内结点NLi的位置PLi。ΡLi=(aLi⋅ΡL1i+bLi⋅ΡL2i+cLi⋅ΡL3i)+dLi⋅VLi(17)式(17)就是最终的结点位置,不需要优化。3.4生成单元当所有的内结点都已经生成,并且组织成肋线后,就可以生成单元。源网格中的每个四边形都与下一层相同肋线上的4个结点构成一个六面体单元。重复此过程直到所有单元都生成为止。4面体网格划分图3为图1～2中几何体的网格,注意该几何体的源和目标曲面具有不同的