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

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

王东风,翟建军,陈文亮

(南京航空航天大学机电学院,江苏南京　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 =i

min g+1[(Δx

m )

2

+(Δy m )2+(Δz m )2]

1/2

∑

i

max

m =i

min 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 =n

i =0

x

i

・

(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 )