有限元网格结点编号_欧阳兴

收稿日期:2000-10-22
基金项目:国家高技术863计划资助项目(863-511-820-020) 作者简介:欧阳兴(1968-),男,江西都昌人,博士生,100083,北京.
有限元网格结点编号
欧阳兴 陈中奎 施法中
(北京航空航天大学机械工程及自动化学院)
摘 要:在有限元分析中,求解高阶线性代数方程组时整体刚度矩阵所需存储与由网格结点编号决定的顺序有关.在基于等带宽存储的求解法与基于变带宽存储的求解法的基础上推导出它们的关系.据此,提出了有限元网格结点编号的前沿法与矩形法,并给出了这两种编号法的内存消耗与结点数量的关系.理论分析和实例表明这两种编号法能有效地减少计算机内存消耗.
关 键 词:有限法;刚度矩阵;线性方程;结点编号;排序中图分类号:TB 115
文献标识码:A 文章编号:1001-5965(2002)03-0339-04
1 问题的提出
在有限元分析中,由于求解如下形式的高阶线性代数方程组,需要耗费计算机大量的内存和计算时间:
K x =f
(1)
其中 K 是整体刚度矩阵;x 是未知向量;f 是已知向量.K 具有大型、对称、带状、稀疏、正定和主
元占优等特点.有限元分析的求解效率很大程度上取决于方程组(1)的求解方法,包括K 的存储方法.因为(1)的求解时间在有限元分析的求解时间中占了很大的比重.当单元增多、结点增多、网格加密和未知数大量增加时,尤为如此.为了提高有限元分析效率,必须结合K 的特点,研究出方程组(1)的高效求解方法.
方程组(1)的求解方法主要有两大类:直接求解法和迭代求解法.直接求解法以高斯消元法为基础,求解效率高.迭代求解法有赛德尔迭代法和超松弛迭代法.高斯消元法是目前求解线性方程组最普遍的方法,有一般高斯消元法和基于带状存储的高斯消元法.当K 的阶数非常高时,为了减小内存消耗,考虑到K 的对称性和带状性,带状存储法只存储以主对角线为中心的斜带形区域的半边.带状存储法有等带宽存储法和变带宽存储法.基于变带宽存储的求解法与基于等带宽存储的求解法相比算法更复杂,但内存消耗更少,因
而有限元分析中越来越多地采用基于变带宽存储
的求解法,如活动列求解器(active column solv -ers)
[1]
和前沿求解法(frontal solution methods)
[2]
.
带状存储法中,(1)所需要的存储空间不仅与(1)中x 的个数有关,而且与x 的顺序有关.而
(1)中x 的分量次序与有限元网格结点的编号有关,不同的有限元网格结点编号很可能使(1)所需要的存储空间差别很大.文献[3],[4]提出了有限元网格结点的编号算法.这些算法是通过减小K 的带宽来减少方程组(1)所需要的存储空间,因而只适合于等带宽存储法.本文提出了一种新的有限元网格结点编号算法,这种算法仅与网格模型中结点空间位置及结点的连接关系有关,与网格模型中结点初始编号无关,既适合于等带宽存储法,又适合于变带宽存储法,并且使(1)需要的内存空间非常少.
2 结点编号问题的数学描述
设一个有限元网格模型总共n 个结点和e 个单元,结点编号和单元编号都从0开始,对于i =0,1,,,n -1,具有编号为i 的结点记为v i ,对于i =0,1,,,e -1,具有编号为i 的单元记为u i .设A =(a ij )n @n 为表示有限元网格模型的结点相连矩阵,a ij 它的第i 行第j 列元素,因而:
1)a ij =1,如果i =j 或者v i ,v j 在同一单元内;
2002年6月第28卷第3期北京航空航天大学学报
Journal of Beijing University of Aeronautics and Astronautics June 2002Vol.28 No 13
2)a ij=0,其它.
A是n阶对称方阵.如果每个结点都仅有一个自由度,则不仅K与A阶数相同(都是n),而且它们的零元素与非零元素对应位置完全相同,即它们具有相同的稀疏结构.设每个结点的自由度为f,在A=(a ij)n@n中,将每个1换成f@f个1组成的方阵,将每个0换成f@f个0组成的方阵,这样得到f@n阶方阵记B,则B与K具有相同的稀疏结构,因此有限元网格结点编号算法只需考虑A.
以下元素都是针对A,行标和列标都从0开始.
对于i=0,1,,,n-1,记第i行第1个元素1的列标为s i,显然
s i=min{j|v j与v i在同一单元内,
j=0,1,,,n-1}(2)又记
l i=i-s i+1(3) s i称为结点v i的最小相关结点号,l i称为第i行半行带宽或称为结点v i的半行带宽.
对于j=0,1,,,n-1,记第j列第1个元素1的行标为t j,称c j=j-t j+1为第j列半列带宽.
A的半带宽B可以从以下两式得到:
B=max{l i|i=0,1,,,n-1}(4)
B=max{c j|j=0,1,,,n-1}(5)对应于等带宽存储法,A的半带形区域(包括主对角线)元素个数为
A0=nB-B(B-1)P2(6)此时,K的半带形区域(包括主对角线)元素个数为
K0=[nB-B(B-1)P2]f2(7)对应于变带宽存储法,A的半带形区域(包括主对角线)元素个数为
A1=E n-1i=0l i=E n-1i=0(i-s i+1)(8)此时,K的半带形区域(包括主对角线)元素个数为
K1=f2E n-1i=0(i-s i+1)(9)显然A0[A1,K0[K1,因此,变带宽存储法比等带宽存储法更节省内存.
由n个结点组成的有限元网格有n!个不同的结点编号,对于这n!个不同的结点编号中任意一个结点编号(记为h,h=0,1,,,n!-1),由有限元网格模型中单元与结点的相互关系,可以立即计算其对应的B,A0,A1,K0和K1.可将B, A0,A1,K0和K1看成以h为自变量的函数,记为B(h),A0(h),A1(h),K0(h)和K1(h).因此,对应于等带宽存储法和变带宽存储法,有限元网格模型的结点最佳编号是分别求解g,使
B(g)=min{B(h),h=0,1,,,n!-1}
(10)和
A1(g)=min{A1(h),h=0,1,,,n!-1}
(11)
虽然n!是有限数,但对于一般的有限元网格模型来说,n!非常大,所以不可能通过直接求出n!个B(h)及A1(h),再比较n!个B(h)或A1(h)的大小,来求出h.本文根据结点与单元的相互关系导出以下算法.
3算法原理
算法用到的数据结构包括:整个网格模型的数据包括结点数组和单元数组,每个结点的数据包括结点坐标、结点编号和结点所属的单元数组,每个单元的数据包括单元编号和依序构成单元的结点数组.
由(2)、(3)、(4)、(6)、(8)、(10)、(11)等式可归纳出两种结点编号方法.
311前沿法
先寻找一个边角结点.因边角结点所在的单元数少,可以通过比较结点所属的单元数组的大小,使结点所属的单元数组最小的结点编号为0,依次给0,1,,,n-2号结点所在单元的未编号结点进行编号.图1采用4@5个四边形单元组成的网格模型的结点编号作为前沿法结点编号的示意.
242322212029
151413121928
876111827
325101726
01491625
图1前沿法结点编号
312矩形法
对矩形区域来说,采用前沿法的图1不是最佳结点编号,图2才是最佳结点编号.一般地,若整个网格模型是一矩形区域,由(a-1)(b-1)个
340北京航空航天大学学报2002年
四边形单元组成,共有ab个结点,排列成a行b 列.若按图2所示进行结点编号,可得如下结论.
4914192429
3813182328
2712172227
1611162126
0510152025
图2矩形区域的最佳结点编号
由(4)式得:
B=a+2(12)由(6)式得:
A0=(a+2)ab-(a+2)(a+1)P2(13)
由(8)式得:
A1=1+2+,+2+(a+1)+(a+2)+,+ (a+2)+,+(a+1)+(a+2)+,+
(a+2)=ba2+2ba-a2-b(14)当b=a时,由(14)式得:
A1=a3+a2-a(15)在(14)式中,交换a与b可得:
A1=ab2+2ab-b2-a(16)因为当a<b时,ba2+2ba-a2-b<ab2+ 2ab-b2-a.因此,图3中先按行对结点进行编号比图2中先按列对结点进行编号得到的A1和B要大.
在图2中,将从左到右作为x轴正向,从下到上作为y轴正向,定义关系/<0为p(x0,y0)< q(x1,y1),当且仅当x0<x1或x0=x1且y0<y1.则图2所示的结点编号可看作结点从小到大的排序.
对于一般的网格模型,虽然单元和结点不如图1~图3所示的矩形区域的单元和结点那样排列得非常规则,但由于同一单元的所有结点的空间位置接近,因而可以采用图2所示的方法进行结点编号.若一般的网格模型中结点、单元的连接关系不与图2所示的矩形区域的单元和结点连接关系相似,则需计算多个排序方法所得编号对应的B或A1,再比较求出最好的结点编号.以下叙述求一般的网格模型的最佳的结点编号方法.
242526272829
181920212223
121314151617
67891011
012345
图3矩形区域的行优先编号
如果待编号网格模型是平面网格模型,即所有结点在同一平面内,则可认为网格模型中所有结点有一个坐标分量是同一常数.否则,经平移、旋转坐标系变换可将所有结点的一个坐标分量变换成同一常数.不妨设所有结点z坐标相同,因而z坐标与结点编号无关,下面只有两个坐标的点可认为省写了z坐标.分别定义关系/<0为
1)p(x0,y0)<q(x1,y1)当且仅当x0<x1或x0=x1且y0<y1.
2)p(x0,y0)<q(x1,y1)当且仅当x0<x1或x0=x1且y0>y1.
3)p(x0,y0)<q(x1,y1)当且仅当y0<y1或y0=y1且x0<x1.
4)p(x0,y0)<q(x1,y1)当且仅当y0<y1或y0=y1且x0>x1.
按这4种结点/<0关系的定义,分别排序得到对应的4种结点编号及对应的B或A1,再比较求出最好的结点编号.
如果待编号网格模型不是平面网格模型,则必须把结点按x,y,z3个坐标排序编号.和上述方法一样,因为x,y,z3个坐标分量的排列次序可以任意,而且每个坐标分量可以从大到小也可以从小到大.因而有48种3坐标字典排序方法来定义结点/<0关系.同样分别排序得到对应的48种结点编号及对应的B或A1,再比较求出最好的结点编号.
313前沿与矩形结合法
大部分情况下矩形法比前沿法更节省内存,只有在一些特殊情况下前沿法比矩形法更节省内存.无论如何,都可以将两种方法结合起来,即通过计算和比较两种方法对应的B或A1,求出最好的结点编号.
4实例与总结
本文处理的一个实例中,整个网格模型有6684个结点和7321个单元,其中三角形单元1452个,四边形单元5869个.
采用变带宽存储法并且实数用双精度浮点数(double)存贮.每个结点自由度为5个.未经本文结点编号而采用初始结点编号需要5801511926MB内存;应用本文前沿法,结果只需要1271346289MB内存;应用本文矩形法,结果只需要861866638MB内存.
341
第3期欧阳兴等:有限元网格结点编号
由(15)式可知,当网格模型基本呈正方形时应用矩形法,K所需要存贮元素个数为O(n3P2).同样可推出,当网格模型基本呈正方形时,应用本文前沿法,K所需要存贮元素个数为O(n3P2).由(14)式可知,当网格模型为长条矩形,K所需要存贮元素个数比网格模型基本呈正方形时还要少.应用矩形法,当网格模型形状非常复杂时与网格模型形状为正方形时基本接近.而大多数结点编号算法K所需要存贮元素个数为O(n2).因而本文算法明显优于这些算法.
一般地,设每个实数需要r个字节内存(采用双精度浮点实数r=8,采用单精度浮点实数r= 4),每个结点的自由度为f,采用矩形法K所需要内存大约为r f2n3P2,采用前沿法K所需要内存大约为4P3r f2n3P2.
参考文献
[1]Bathe K J.Fi nite elenent procedures in engineeri ng analysis[M].
Prentice-Hall,Inc,Enge wood Cli ffs,NJ,1982.441~449.
[2]Irons B M.A frontal solution program for finite element analysis[J].
Int J Numer M eth Engng,1986,23:239~256.
[3]Cuthill E,M ckee J.Reducing the bandwidth of s parse symme tric
ma xtrices[A].Proc24th Nat Conf of the ACM[C].1969.157~ 172.
[4]Gibbs N E,Poole W G,Stockmeyer P K.An al gori thm for reduce
the bandwidth and profile of a sparse matrix[J].SIMA J Num Anal, 1976,13(2):236~250.
Numbering of Fin ite Element Mesh Nodes
OUYANG Xing C HEN Zhong-kui SHI Fa-zhong
(Beijing University of Aeronautics and As tronautics,School of Mechanical Engineering and Automati on) Abstract:In finite element analysis,the storage needed by a total stiffness matrix for solving a large-scale sys-tem of linear equations is related to the sequenace determined by numbering of mesh nodes.Based on the solutions for both constant and varible bandwidth storages,their relationship is derived.On the above basis,the frontal meth-od and rectangle method are proposed,the relationships between memory spending for both methods and a mount of nodes are given.Theoretical analyzing and practical examples have proved that the two methods can efficiently de-crease the memory spending of computer.
Key words:finite element methods;stiffness matrix;linear equations;numbering of nodes;sorting 342北京航空航天大学学报2002年。