有限元刚度矩阵的压缩存贮组集及快速求解(完整版)实用资料.doc

有限元刚度矩阵的压缩存贮组集及快速求解（完整版）实用资料
(可以直接使用，可编辑完整版实用资料，欢迎下载）

有限元刚度矩阵的压缩存贮组集及快速求解
姚松田红旗
1中南大学轨道交通安全教育部重点实验室,湖南长沙,410075

摘要:基于“细胞元”索引存贮方案提出了一种仅组集有限元刚度矩阵中“非零元素”的方法,该方法最突出的特点是计算所需内存空间与有限元网格节点和单元的编号模式无关,适于进行“自适应网格细化”有限元分析。针对刚度矩阵的“一维压缩存贮”格式,对稀疏矩阵直接解法和预处理共轭梯度法进行了探讨,并编制了相应的计算机程序对某地铁车辆有限元模型进行了分析,计算结果与ANSYS5.7的计算结果一致,相差不超过2%,说明提出的存贮方案和求解方法是正确可靠的。
关键词:细胞元;组集;非零元素;稀疏直接求解法;预处理共轭梯度法
1概述
有限元方法自1956年首次应用于飞机结构分析以来,由于其具有几何适应性强、易于处理非线性等优点,已成为科学和工程计算领域应用最为广泛的数值方法之一[1]。在结构有限元分析中,对求解域进行空间离散以后,经常产生大型或巨型的稀疏线性方程组A•X=F,对称正定的系数矩阵A通常称为结构总体刚度矩阵,其大小为ND×ND,ND为结构的总自由度数,刚度矩阵A所需的存贮空间随结构规模增大而增加极快,以至于限制了求解规模。目前,刚度矩阵的存贮方案[2-4]多采用二维等带宽存贮或Profile/Envelope结构(即一维变带宽。采用二维等带宽存贮时将有可能由于局部带宽过大而使整体刚度系数矩阵的存贮量大大增加;一维变带宽存贮虽然会比等带宽存贮节省内存空间,但是对于Skyline下的零元素还是要加以存贮,且存贮的零元素数量取决于单元和节点的编号模式。
毫无疑问,仅对整体刚度矩阵A中的非零元素进行存贮是最节省存贮空间的方法[5],由于仅对其中的非零元素进行操作,从而提高了计算效率。另外该存贮格式具有如下特点:无论是对存贮还是方程的求解来说,所需要的内存空间与有限元网格划分时节点和单元的编号模式无关。因此该方法高度适合“h-自适应网格细化”有限元分析,由于没有必要对新产生的节点和单元进行重新编号从而降低了计算要求。
本文提出了基于“一维压缩存贮”的整体刚度矩阵组集方法,并针对刚度矩阵的压缩存贮格式,讨论了大型有限元方程组的快速求解方法。
2结构稀疏刚度矩阵的存贮
对于某地铁车辆的车体有限元模型,结构总自由度ND为266862。由于刚度矩阵的对称正定性,仅需存贮整体刚度矩阵的下三角(或上三角矩阵。整体刚度矩阵的最大半带宽为263982,表1列出了采用不同存贮格式时所需要存贮的元素总数以及零元素的个数:
1教育部博士点基金(20533007项目资助
表1采用不同存贮格式时所需存贮元素的数目
存贮格式零元素个数存贮元素总数

一维变带宽635136972640022307
本文方法04885335由上可知,无论是二维等带宽还是一维变带宽格式,都存贮了大量多余的零元素,而本文提出的一维压缩存贮格式将所有非零元素压缩成一维向量存贮,最大可能地节省了存贮空间。同时采用辅助向量来描述非零元素在结构矩阵A中的位置等信息[6]。在本文中采取三个向量来存贮整体刚度矩阵。分别是:
Values:长度为NZ的双精度向量,存贮稀疏矩阵下三角矩阵中的NZ个非零元素,元素按“所在行”的顺序依次存贮,每一行以对角元素结束。
Column,长度为NZ的整数向量,存贮对应Values向量中每个元素所在的列号。RowIndex,长度为(ND+1的整数向量,存贮每一行起始元素的指针,根据RowIndex(I+1和RowIndex(I的差值可以计算出第I行非零元素的个数。
3整体刚度矩阵的组集
记节点自由度为DOF,通常DOF>1;考虑到结构有限元分析中节点具有多个自由度的特征(注意:每个节点的自由度是连续编号的,我们将整体刚度矩阵按照节点(而不是节点自由度形成分块子矩阵[7-10],在本文中称该子矩阵为“细胞元”。根据有限元法的基本原理,如果编号为I,J的两节点不共存于同一个单元中,那么在整体刚度矩阵中对应于(I,J位置的“细胞元”(大小为DOF×DOF的子矩阵必然为零矩阵。它们占据整体刚度矩阵中位置如下:[DOF×(I-1+M,DOF×(J-1+N],其中1≤M,N≤DOF。当且仅当I,J两节点共存于同一单元中时,整体刚度矩阵中对应于(I,J位置的“细胞元”才不为零。
由上述原理,根据节点-节点之间的拓扑连接关系就可大致推断出整体刚度矩阵中非零元素的分布情况;由于仅仅需要存贮对称整体刚度矩阵的下三角矩阵(或上三角矩阵就可以了。那么只要对编号为I的节点,找出所有与I节点共处于同一单元且编号比I小的相关节点即可。图1中列出了地铁有限元模型中相关节点数最大的节点390