车辆薄板有限元分析中的多因子不完全分解预处理解法

车辆薄板有限元分析中的多因子不完全分解预处理解法姚 松 田红旗1中南大学轨道交通安全教育部重点实验室，湖南长沙，410075dynacn@摘要：薄板是轨道车辆结构的主要形式，本文基于离散Kirchhoff 假设的DKT 弯曲板单元推导了四边形弯曲板单元DKQ 的构造过程，并进一步阐述了用于一般薄板问题分析的平板单元的构造。

提出了一种“多因子不完全分解” 的预处理方法，与共轭梯度迭代法结合能够大大加快薄板问题大型稀疏方程组的收敛速度，经过数值试验，说明该方法是稳定可靠的。

该方法避免了常规不完全分解不适用于薄板这样的 “病态”结构的情况。

在此基础上，编写了一般薄板问题分析的有限元程序，程序对结构刚度矩阵采用压缩存贮的方法，节约了大量内存空间。

本文还对分解算法中的可选参数进行了优化研究。

通过一个数值试验，本程序计算结果与商业有限元软件ANSYS5.7的结果完全一致。

关键词：薄板结构，DKQ 单元，预处理，不完全分解，共轭梯度法1 概 述有限元单元法已经成为结构分析的重要方法，薄板结构是轨道车辆的主要结构形式，因此薄板结构有限元分析已成为车辆结构分析中的重大课题。

早期的弯曲板单元大多基于经典的薄板理论，在以该理论为基础的板单元的能量泛函中，包含位移的二阶偏导数，要求位移为类连续。

这给构造板单元带来了困难，由此研究人员将注意力转向了中厚板单元，大多采用中厚板理论，其能量泛函仅包含位移的一阶导数，只要求位移是类连续，但是用厚板理论建立的单元仅对中厚板有效，当板逐渐变薄时，单元刚度矩阵中的剪切项占主导地位，计算出的弯曲变形远小于实际变形；当板非常薄时，求得的位移趋向于零，从而产生了“剪切闭锁”现象。

1C issner Mindlin Re −0C 基于离散的假设，通过挠度和转角分别独立插值，然后在若干个离散点上强迫挠度与转角满足薄板经典理论中的约束，构造出三角形（DKT ）和四边形（DKQ ）薄板弯曲单元，其泛函的表达式又回复为经典薄板理论的泛函表达式，又自然解决了“剪切闭锁现象”问题。

多个文献表明DKT 元与DKQ 元在求解薄板弯曲问题时都显示出良好的性能，具有较高的精度。

在对实际车辆结构进行有限元分析时，由于结构受力复杂，在承受板平面内的载荷的同时，也有可能板平面外的载荷，因此在进行分析时所采用的平板单元是平面应力单元与DKT 弯曲单元的组合而成。

由于三角形平面应力单元为常应变单元，为了提高分析的精度，在本文中我们讨论由四边形膜单元和DKQ 单元组合而成的平板单元。

Kirchhoff Kirchhoff1 教育部博士点基金（20020533007）项目资助1采用有限元法求解薄板弯曲问题最终归结为求解一组稀疏对称正定的线性方程组，，b u K =⋅K 为整体刚度矩阵，为待求解的位移向量，b 为载荷向量。

有限元求解主要分为直接求解器和迭代求解两大类u [1-2]，直接求解是当前应用最为广泛的求解技术，其存贮方案多采用一维变带宽，通过对总体刚度矩阵直接进行或分解，然后再回代求解，采用这类技术经过长期使用比较成熟，但是该方法的劣势在于：分解后不再是稀疏矩阵，在分解的过程中会产生大量的“填入”元，因此对于大型稀疏矩阵的分解不仅耗费时间，而且占用内存。

在求解大型结构问题时速度比较缓慢，而且所需存贮空间和计算量随结构规模增大而急剧增加，以致于限制了求解规模。

T LL TLDL L Cholesky 对于象整车结构分析问题，有限元离散方程组的阶数可以达到几十万阶，采用迭代算法可以仅仅保存刚度矩阵中的非零元素，由给定初值通过若干迭代步骤获得满足一定精度的近似解。

传统的迭代法包括：，Jacobi seidel Gauss −，SOR 等等，这些方法收敛速度过慢而且没有保证，实际运用不多。

以预条件共轭梯度法[4-11]（preconditioned Conjugate Gradient Method ，简称为PCG ）为代表的迭代法是近十几年来逐渐兴起并开始得到应用的一类迭代方法，PCG 法的基础是共轭梯度法（CG ）， CG 方法的收敛速率取决于条件数，当矩阵K 的条件数接近1时，CG 法的迭代速度很快，而当矩阵K 的条件数()210>K Con 时，CG 法的迭代就非常慢。

为了提高收敛速度，必须通过使用预条件技术把原先的方程组转换成一个等价的，但是系数矩阵条件数更小、更易于收敛的方程组。

即选择对称正定矩阵M ，考虑等价方程，若b M Ku M 11−−=K M 1−的条件数比K 要好的话，再运用CG 法求解收敛速度就会很快，问题的关键是如何选取M ，使得谱条件数能够得到较大改善。

在本文中针对薄板这样的病态问题提出了一种预条件器。

2 多因子不完全分解预处理算法Meijerink 和于1977年基于“不完全” 分解提出了预条件算子Vorst der van Cholesky M ，T L D L M ⋅⋅=，其中矩阵为单位下三角阵，，L ()()⎪⎩⎪⎨⎧==≠≠==00001ij ij ij K if K if j i L ()0≠=ij ij ij K if K M ，D 为对角矩阵，该方法对于对称正定且对角占优的矩阵比较有效，能够大大提高CG 法的收敛速度，计算又非常简单。

但是对于薄板这样的“病态”结构，刚度矩阵并不是严格的对角占优矩阵，在分解过程中矩阵的对角线元素可能会出现负值，从而不能保证预条件矩阵D M 为正定矩阵。

在本文中针对薄板问题提出了如下的“多因子不完全分解法”[3]。

将K 矩阵分解为对角矩阵和非对角矩阵的叠加，B D K −=，其中矩阵为D K 矩阵的对角矩阵，构造出下式：()()FBF B FBF D K −−−=，其中：，()，其中为可选择的参数乘子。

记矩阵i d n d d d d diag F L L 321,,=10<<i d 2()FBF D −为W ，可以看到：()ij j i ij n k n m mj km ik ij ij ij ij B d d D F B F D FBF D W −=−=−=∑∑==11由上式可以得到：，可以看出W 与有相同的稀疏结构，其对角元素与矩阵相同，在非对角元素上由于有和两个小乘子，因此有。

因此矩阵也是一个对称正定矩阵。

的计算公式如下：⎩⎨⎧≠==j i K d d j i D W ijj i ii ij A A i d j d ij ij A W <W i d 1,11>=t td 为可选择参数 对： n i ,,2L =∑−==1122i m mm im m a K d s ，若2t K s ii >，则有s K t d ii i 1=，若2t K s ii ≤，则t d i 1= 当K 对称正定，在的情况下，W 的不完全分解存在，记为，1>t E L L W T−⋅⋅=σE 为误差矩阵，且对角元素i σ均大于零。

那么有()()FBF B E L L FBF B W K T −+−⋅⋅=−−=σ，将的不完全分解也可作为总体刚度矩阵W T L L ⋅⋅σK 的不完全分解，而且由于对角元素i σ均大于零，该分解可以作为预条件算子，其算法如下：111W =σjj m jm m im ij ij L L W L σσ⎟⎟⎠⎞⎜⎜⎝⎛⋅⋅−=∑−=11 m i m im ii i L W σσ⋅−=∑−=112得到上述预条件矩阵后，将薄板结构有限元方程转化为：，其中，，，通过这样的变换，由于误差矩阵的忽略，实际上T y C =⋅()TL K L C 111−−−⋅⋅⋅=σx L y T ⋅=b L T ⋅⋅=−−11σ(FBF B E −+)I C ≠，但是其条件数与K 相比有明显改善。

以上述变换公式代入CG 方法，即可以求出结构位移向量，流程如下：00x K b r ⋅−=， 为第一次迭代值的余差010r M p ⋅=−0r 0x 那么第次迭代为：k ()()()k T k k T k k p K p p r ⋅⋅=αk k k k p x x ⋅+=+α1 k k k k p K r r ⋅⋅−=+α1()(()())k T k k T k k r r r r ⋅⋅=++11β k k k k r r p ⋅+=++β11迭代过程直至和二者的差别“足够”小时，计算过程就收敛了。

1+k x k x 3从上面的流程可见，整个计算流程只要求计算刚度矩阵K 与一个列向量的乘积，而刚度矩阵K 对称而且各行的零元素对于乘积是没有任何贡献的，因此可以按照压缩存贮的格式将整体刚度稀疏矩阵的下三角矩阵“按行或按列”的方式压缩存贮，在本文中采用了实数型数组存放非零元素的值，整型数组存放非零元素所在的列，整型数组存放“每一行或每一列”起始非零元素的位置，其中为结构刚度矩阵下三角矩阵中非零元素的总数，ND 为结构总自由度数。

采用这样的压缩存贮方式能够最大限度地节约内存空间。

)(NZ stiffvalue )(NZ newcolumn )1(_+ND lower total NZ 3 DKQ 弯曲板单元及平板单元DKQ 四边形弯曲板单元基于四个DKT 三角形弯曲板单元[12]，如下图1所示：图1 DKQ 单元 图2 DKT 单元根据虚功原理，四边形单元（1234）的虚功可以由上图所划分的①、②、③、④四个三角形的虚功相叠加之后除2。

因此DKQ 单元的构造是通过组成四边形的四个三角形DKT 单元的构造组集以后得到。

DKT 单元采用和w x θ，y θ的独立插值，w 和x θ，y θ之间的约束方程通过在三角形的三个角节点和三个中节点位置强迫实现，因此其泛函表达式恢复为经典薄板理论的泛函表达式：∫∫∫∫ΩΩ⋅⋅−⋅⋅⋅=Πdxdy w q dxdy D T κκ21， κ为板弯曲的广义应变 每个角节点有参数,i w xi θ，yi θ()3,2,1=i ,一共9个自由度构成单元的节点位移向量，边中节点有参数xi θ，yi θ()6,5,4=i ，直法线的假设通过下述方式引入：Kirchhoff 在角节点：(3,2,100=⎪⎪⎩⎪⎪⎨⎧=−⎟⎟⎠⎞⎜⎜⎝⎛∂∂=−⎟⎠⎞⎜⎝⎛∂∂i y w x w yi i xi i θθ) 在中间节点：()()6,5,4210=⎪⎪⎩⎪⎪⎨⎧+==−⎟⎠⎞⎜⎝⎛∂∂k s w nj ni nk sk k θθθθ 通过约束条件将各边中间节点的转角参数“凝聚”掉，不但提高了位移场的精度，而且保证了相邻单元之间的协调性。

将弯曲板单元和平面膜单元组合起来就得到了平板单元的单元刚度矩阵，在局部坐标系中，节点位移参数本不包含zi θ，但是为了进行总刚集成时避免共面单元组集后出现zi θ方向刚度为零的情况，在组合过程中给zi θ方向赋予单位刚度，组集过程如下图3：4图3 板单元刚度组集示意图图4 薄板例题4 数值试验针对薄板有限元问题，本文编制了一个有限元程序，该程序采用压缩存贮的方法，将整体刚度矩阵中的非零元素存贮在三个一维数组中。