基于映射法的六面体网格生成算法

基于映射法的六面体网格生成算法王东风,翟建军,陈文亮(南京航空航天大学机电学院,江苏南京　210016)摘要:六面体网格划分技术是三维有限元仿真软件处理的关键环节之一,等参映射法既可适应特殊的区域边界形状,又可控制所生成单元的形状和密度。

对基于等参映射法的六面体网格划分原理进行了深入研究,并在此研究基础上对等参映射法的计算过程进行了细致的分析,利用VC++开发了该算法的相应程序,最后给出了2个等参映射法具体的应用实例,计算结果表明该程序的计算精度已经达到了工程要求。

关键词:等参映射法;六面体网格;有限元中图分类号:TP391 文献标识码:A 文章编号:1672-1616(2009)05-0025-03 在有限元仿真过程中,单元类型的选择对整个有限元仿真的计算效率、自动化程度、计算精度都将产生重要影响。

六面体单元由于变形特性好、计算精度高等优点在三维有限元仿真领域中得到了广泛应用[1]。

映射法是三维网格划分中最早使用的方法,和扫略法、基于栅格法等其他方法相比,该方法生成网格速度快、生成的网格单元质量好、网格密度可控制[2～4]。

映射法对复杂实体生成三维有限元网格有两大难点,一是子区域划分问题,二是子区域之间网格相容性问题。

Price与Armstrong等提出中面法,将三维复杂区域分解成可映射子区域[5～7],但是该算法存在一些问题,特别是几何适应能力问题。

李华和程耿东提出了三维组合式模板,一定条件下解决了子区域之间的网格相容性问题[2]。

还有学者提出了Embedded Voronoi Graph[8]和BLOBs[9],对复杂实体利用映射法划分六面体网格。

映射法在众多有限元分析软件中占有重要地位,美国Altair公司Hyper-Mesh软件中的Solid Mesh Panel就是利用映射法生成六面体网格。

本文对基于等参映射法的六面体网格划分技术进行了详细研究。

通过形函数映射技术将物理域映射到参数空间域,对规则参数域进行网格剖分,将参数域的网格反向映射回物理空间,从而得到物理空间六面体网格。

利用VC++实现了映射过程,在输入边界信息和划分信息后,即得到了六面体网格的节点信息和单元信息。

1　映射法生成四边形网格和六面体网格本文主要讨论的是怎样在一个子区域中划分六面体网格。

这里的子区域指的是具有6个面12条边,每条边的特征点已知的区域。

求子区域六面体网格节点的步骤:(1)利用积累弦长参数化法对每条边进行参数化。

(2)利用拉格朗日插值公式求边界函数。

(3)利用边界函数,由双线性混合孔斯曲面片公式求曲面节点坐标。

(4)由孔斯线性混合插值公式求子区域节点坐标。

在计算的过程中要用到2种坐标系,即笛卡尔坐标系和自然坐标系。

笛卡尔坐标系用x,y,z表示,自然坐标系用一组不超过1的无量纲参数r,s, t表示,边界点分别对应自然坐标等于1或0的点。

如图1所示。

图1　自然坐标和笛卡尔坐标之间的变换收稿日期:2008-08-08作者简介:王东风(1979-),男,河南商丘人,南京航空航天大学硕士研究生,主要研究方向为CAD/CAM/CAE。

1.1　四边形网格在边1,8,4,5构成的曲面中,s =0表示边界1,s =1表示边界4,同理r =0表示边界5,r =1表示边界8(如图1所示)。

定义f 1(r ),g 1(r ),h 1(r )作为边界函数,它表示在第一条边根据不同r 值分别得到不同的x ,y ,z 值。

因此,边界1(s =0):x =f 1(r ),y =g 1(r ),z =h 1(r );边界4(s =1):x =f 4(r ),y =g 4(r ),z =h 4(r );边界5(r =0):x =f 5(s ),y =g 5(s ),z =h 5(s );边界8(r =1):x =f 8(s ),y =g 8(s ),z =h 8(s )。

1.1.1　求边界函数每条边的特征点坐标已知,利用积累弦长参数化法对每条边进行参数化[10],即:r ij min =∑i m =imin g+1[(Δxm )2+(Δy m )2+(Δz m )2]1/2∑imaxm =imin g+1[(Δx m )2+(Δy m )2+(Δz m )2]1/2(1)式中:i ,j 为两个索引,i 随着r 的增大而增大,j 随着s 的增大而增大。

同理可以求得r ij ,s ij 。

定义f 1(r ij min )=x ij min ,g 1(r ij min )=y ij min ,利用拉格朗日插值求边界函数。

f (r )=∑i =ni =0xi・(r -r 0)(r -r 1)…(r -r i -1)(r -r i +1)…(r -r n )(r i -r 0)(r i -r 1)…(r i -r i -1)(r i -r i +1)…(r i -r n )(2)1.1.2　利用边界函数求曲面节点坐标由双线性混合孔斯曲面片计算公式求曲面节点坐标:x (r ,s )=(1-s )f 1(r )+sf 2(r )+(1-r )f 3(s )+rf 4(s )-x (0,0)(1-r )(1-s )-x (1,0)r (1-s )-x (0,1)(1-r )s -x (1,1)rs(3)类似方法可以求得y (r ,s ),z (r ,s )。

将r ,s 按一定规律分割可以得到所有节点坐标,由于坐标变换不改变几何拓扑关系,按一定顺序连接节点生成四边形网格。

1.2　六面体网格的生成如图1所示子区域的边界函数为:1　(s =0,t =0)　x =f 1(r ),y =g 1(r ),z =h 1(r )2　(s =0,t =1)　x =f 2(r ),y =g 2(r ),z =h 2(r )3　(s =1,t =1)　x =f 3(r ),y =g 3(r ),z =h 3(r )4　(s =1,t =0)　x =f 4(r ),y =g 4(r ),z =h 4(r )5　(r =0,t =0)　x =f 5(s ),y =g 5(s ),z =h 5(s )6　(r =0,t =1)　x =f 6(s ),y =g 6(s ),z =h 6(s )7　(r =1,t =1)　x =f 7(s ),y =g 7(s ),z =h 7(s )8　(r =1,t =0)　x =f 8(s ),y =g 8(s ),z =h 8(s )9　(r =0,s =0)　x =f 9(t ),y =g 9(t ),z =h 9(t )10　(r =0,s =1)　x =f 10(t ),y =g 10(t ),z =h 10(t )11　(r =1,s =1)　x =f 11(t ),y =g 11(t ),z =h 11(t )12　(r =1,s =0)　x =f 12(t ),y =g 12(t ),z =h 12(t ) 利用孔斯线性混合插值公式,可求子区域笛卡尔坐标值[10]:x (r ,s ,t )=(1-s )(1-t )f 1(r )+(1-s )tf 2(r )+stf 3(r )+s (1-t )f 4(r )+(1-r )(1-t )f 5(s )+(1-r )tf 6(s )+rtf 7(s )+r (1-t )f 8(s )+(1-r )(1-s )f 9(t )+(1-r )sf 10(t )+rsf 11(t )+r (1-s )f 12(r )+c (r ,s ,t )(4)其中c (r ,s ,t )=-2[(1-r )(1-s )(1-t )x (0,0,0)+(1-r )(1-s )t x (0,0,1)+(1-r )s (1-t )x (0,1,0)+(1-r )st x (0,1,1)+r (1-s )(1-t )x (1,0,0)+(r (1-s )t x (1,0,1)+rs (1-t )x (1,1,0)+rst x (1,1,1)]。

同理可求y (r ,s ,t ),z (r ,s ,t )。

如图1所示,几何体有6个面:1　(s =0)　x =f 1(r ,t ),y =g 1(r ,t ),z =h 1(r ,t )2　(s =1)　x =f 2(r ,t ),y =g 2(r ,t ),z =h 2(r ,t )3　(r =0)　x =f 3(s ,t ),y =g 3(s ,t ),z =h 3(s ,t )4　(r=1)　x=f4(s,t),y=g4(s,t),z= h4(s,t)5　(t=0)　x=f5(r,s),y=g5(r,s),z= h5(r,s)6　(t=1)　x=f6(r,s),y=g6(r,s),z= h6(r,s)由公式(3)可知:f1(r,t)=(1-t)f1(r)+tf2(r)+(1-r)・f9(t)+rf12(r)-(1-r)(1-t)x(0,0,0)-(1-r)t x(0,0,1)-r(1-t)x(1,0,0)-rt x(1,0,1)(5)同理求出其他面函数。

公式(4)改写为:x(r,s,t)=12(1-s)[(1-t)f1(r)+tf2(r)+(1-r)f9(t)+rf12(t)]+12s[tf3(r)+(1-t)f4(r)+(1-r)f10(t)+rf11(t)]+12(1-r)[(1-t)f5(s)+tf6(s)+(1-s)f9(t)+sf10(t)]+12r[tf7(s)+(1-t)f8(s)+sf11(t)+(1-s)f12(t)]+12(1-t)[(1-s)f1(r)+sf4(r)+(1-r)f5(s)+rf8(s)]+12t[(1-s)f2(r)+sf3(r)+(1-r)f6(s)+rf7(s)]+ c(r,s,t)(6)将公式(5)等其他面函数公式代入公式(6)得:x(r,s,t)=12[(1-s)f1(r,t)+sf2(r,t)+(1-r)f3(s,t)+rf4(s,t)+(1-t)f5(r,s)+tf6(r,s)]+14c(r,s,t)(7)同理可求y(r,s,t),z(r,s,t)。

通过计算,几何体的边界节点和内部节点都已知,连接节点生成六面体网格。

2　应用实例分别选择开始区域的4条边,结束区域的4条边,连接开始区域和结束区域的4条边,输入x,y, z方向网格数量为10,生成六面体网格(如图2所示)。

通过计算单元最小尺寸、单元纵横比、扭曲度、正交性、单元绝对体积,发现所得到的网格绝大部分满足工程要求。

图2　三维有限元网格生成实例3　结束语本文对基于等参映射法的六面体网格划分技术进行了详细研究。

由应用实例可知,等参映射法所得到的六面体网格质量高、密度可控制,但对于曲率较大部位要进行局部网格疏密调整,增加单元数量,提高网格质量。