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

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

预处理解法

姚 松 田红旗

1中南大学轨道交通安全教育部重点实验室，湖南长沙，410075

dynacn@

摘要：薄板是轨道车辆结构的主要形式，本文基于离散Kirchhoff 假设的DKT 弯曲板单元推导了四边形弯曲板单元DKQ 的构造过程，并进一步阐述了用于一般薄板问题分析的平板单元的构造。提出了一种“多因子不完全分解” 的预处理方法，与共轭梯度迭代法结合能够大大加快薄板问题大型稀疏方程组的收敛速度，经过数值试验，说明该方法是稳定可靠的。该方法避免了常规不完全分解不适用于薄板这样的 “病态”结构的情况。在此基础上，编写了一般薄板问题分析的有限元程序，程序对结构刚度矩阵采用压缩存贮的方法，节约了大量内存空间。本文还对分解算法中的可选参数进行了优化研究。通过一个数值试验，本程序计算结果与商业有限元软件ANSYS5.7的结果完全一致。

关键词：薄板结构，DKQ 单元，预处理，不完全分解，共轭梯度法

1 概 述

有限元单元法已经成为结构分析的重要方法，薄板结构是轨道车辆的主要结构形式，因此薄板结构有限元分析已成为车辆结构分析中的重大课题。早期的弯曲板单元大多基于经典的薄板理论，在以该理论为基础的板单元的能量泛函中，包含位移的二阶偏导数，要求位移为类连续。这给构造板单元带来了困难，由此研究人员将注意力转向了中厚板单元，大多采用中厚板理论，其能量泛函仅包含位移的一阶导数，只要求位移是类连续，但是用厚板理论建立的单元仅对中厚板有效，当板逐渐变薄时，单元刚度矩阵中的剪切项占主导地位，计算出的弯曲变形远小于实际变形；当板非常薄时，求得的位移趋向于零，从而产生了“剪切闭锁”现象。

1C issner Mindlin Re −0C 基于离散的假设，

通过挠度和转角分别独立插值，然后在若干个离散点上强迫挠度与转角满足薄板经典理论中的约束，构造出三角形（DKT ）和四边形（DKQ ）薄板弯曲单元，其泛函的表达式又回复为经典薄板理论的泛函表达式，又自然解决了“剪切闭锁现象”问题。多个文献表明DKT 元与DKQ 元在求解薄板弯曲问题时都显示出良好的性能，具有较高的精度。在对实际车辆结构进行有限元分析时，由于结构受力复杂，在承受板平面内的载荷的同时，也有可能板平面外的载荷，因此在进行分析时所采用的平板单元是平面应力单元与DKT 弯曲单元的组合而成。由于三角形平面应力单元为常应变单元，为了提高分析的精度，在本文中我们讨论由四边形膜单元和DKQ 单元组合而成的平板单元。

Kirchhoff Kirchhoff

1 教育部博士点基金（20020533007）项目资助

1

采用有限元法求解薄板弯曲问题最终归结为求解一组稀疏对称正定的线性方程组，，b u K =⋅K 为整体刚度矩阵，为待求解的位移向量，b 为载荷向量。有限元求解主要分为直接求解器和迭代求解两大类u [1-2]，直接求解是当前应用最为广泛的求解技术，其存贮方案多采用一维变带宽，通过对总体刚度矩阵直接进行或分解，然后再回代求解，采用这类技术经过长期使用比较成熟，但是该方法的劣势在于：分解后不再是稀疏矩阵，在分解的过程中会产生大量的“填入”元，因此对于大型稀疏矩阵的分解不仅耗费时间，而且占用内存。在求解大型结构问题时速度比较缓慢，而且所需存贮空间和计算量随结构规模增大而急剧增加，以致于限制了求解规模。

T LL T

LDL 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<

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

j i ii ij A A i d j d ij ij A W

1>=t t

d 为可选择参数 对： n i ,,2L =∑−==112

2i 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 =σ

j

j m jm m im ij ij L L W L σσ⎟⎟⎠⎞⎜⎜⎝⎛⋅⋅−=∑−=11 m i m im ii i L W σσ⋅−=∑−=1

12

得到上述预条件矩阵后，将薄板结构有限元方程转化为：，其中，，，通过这样的变换，由于误差矩阵的忽略，实际上T y C =⋅()T

L 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