CN101404005A - 基于可变预处理迭代获得非比例阻尼力学系统模态频响的方法 - Google Patents

Abstract

基于可变预处理迭代获得非比例阻尼力学系统模态频响的方法，它涉及一种基于Jacobi迭代的非比例阻尼力学系统模态频率响应计算方法，以解决在非比例阻尼系统中计算频率响应的方法存在的计算代价较高、收敛速度和计算精度较差的问题。根据Jacobi迭代的预处理矩阵计算频率响应的初始解，并计算残差，令k＝0；计算第k+1步的预处理矩阵，根据第k+1步预处理矩阵计算第k+1步的频率响应向量，并计算第k+1步的残差；判断第k+1步的残差是否满足残差范数的相对值小于给定的收敛指标，判断结果为是，则迭代计算结束，判断结果为否，则令k＝k+1并重新计算第k+1步的残差。本发明能够加快迭代的收敛速度，保证计算精度，还具有简单易于实施的特点。

Description

基于可变预处理迭代获得非比例阻尼力学系统模态频响的方法

技术领域

本发明涉及一种基于Jacobi迭代的非比例阻尼力学系统模态频率响应计算方法，属于结构动力学特性分析领域。

背景技术

频率响应分析是结构动力学特性分析的主要方法之一。频率响应的计算实际上是对复数线性系统进行求解，对于大型复杂结构，直接对有限元自由度进行求解计算代价非常昂贵，因而通常采用模态频响法。对于非比例阻尼系统，其模态阻尼矩阵通常是满阵，其频率响应在计算求解上同样具有一定的挑战。

目前已有多种求解线性系统的方法，这些方法主要可以分为两类：直接法和迭代法，这些方法都可以用来进行频率响应计算，直接法通过将系统矩阵分解为等价的三角方程组进行求解。这类方法直接、精确但由于矩阵分解需要O(n³)的浮点操作，因而其计算代价非常昂贵；迭代法通过生成一系列收敛于精确解的近似解对系统进行求解，可以用来替代直接法，这类方法具有计算代价小，计算效率高的优点，但收敛速度及精度难以达到要求。

发明内容

本发明为解决在非比例阻尼系统中计算频率响应的方法存在的计算代价较高、收敛速度和计算精度较差的问题，提供一种基于可变预处理迭代获得非比例阻尼力学系统模态频响的方法。本发明包括以下步骤：

步骤一、根据Jacobi迭代的预处理矩阵D⁽⁰⁾＝D，计算频率响应的初始解X⁽⁰⁾，并计算残差r⁽⁰⁾，令k＝0；

步骤二、计算第k+1步的预处理矩阵：

D^(k+1)＝diag(AD^(k)-1r^(k)./D^(k)-1r^(k))

根据第k+1步预处理矩阵计算第k+1步的频率响应向量：

X^(k+1)＝X^(k)+D^(k+1)-1r^(k)

并计算第k+1步的残差：

r^(k+1)＝F-AX^(k+1)

步骤三、判断第k+1步的残差是否满足残差范数的相对值小于给定的收敛指标ε，即 $\frac{\left|\right|r\left|\right|}{\left|\right|F\left|\right|}<\mathrm{\&epsiv;},$ 判断结果为是，则迭代计算结束，此时X^(k+1)即为得到的模态频响；判断结果为否，则令k＝k+1并返回步骤二继续计算，所述k是整数，且k≥0。

有益效果：本发明通过在迭代计算频率响应的过程中选择不同的预处理矩阵，建立了一种求解非比例系统频率响应的可变预处理迭代方法，该方法能够加快Jacobi迭代法的收敛速度，并保证计算精度，同时还具有简单易于实施的特点。

具体实施方式

具体实施方式一：本实施方式由以下步骤组成：

步骤一、根据Jacobi迭代的预处理矩阵D⁽⁰⁾＝D，计算频率响应的初始解X⁽⁰⁾，并计算残差r⁽⁰⁾，令k＝0；

步骤二、计算第k+1步的预处理矩阵：

D^(k+1)＝diag(AD^(k)-1r^(k)./D^(k)-1r^(k))

根据第k+1步预处理矩阵计算第k+1步的频率响应向量：

X^(k+1)＝X^(k)+D^(k+1)-1r^(k)

并计算第k+1步的残差：

r^(k+1)＝F-AX^(k+1)

步骤三、判断第k+1步的残差是否满足残差范数的相对值小于给定的收敛指标ε，即 $\frac{\left|\right|r\left|\right|}{\left|\right|F\left|\right|}<\mathrm{\&epsiv;},$ 判断结果为是，则迭代计算结束，此时X^(k+1)即为得到的模态频响；判断结果为否，则令k＝k+1，并返回步骤二继续计算，所述k是整数，且k≥0。

迭代法的收敛速度取决于系数矩阵的谱特征，因此通常应用预处理技术来改善原系统的谱特征来加速迭代法的收敛速度。本实施方式基于的原理既是在计算频率响应的过程中对残差进行可变预处理迭代计算，直到计算出的收敛指标达到参数，即达到了计算目的。

结构动力学系统频域内的模态运动方程可以表达为：

A(ω)X(ω)＝F(ω)                    (1)

其中A(ω)＝-ω²I+iωB+(1+iη)Λ+iK_s，B和K_s分别是模态粘性阻尼矩阵和模态结构阻尼矩阵，对非比例阻尼系统，二者均为满阵；Λ为系统的特征值对角阵，F(ω)为模态力向量，X(ω)为待求解的模态频率响应向量。

方程(1)可以重新写为：

DX＝CX+F                            (2)

其中矩阵D为A的负对角阵，而C＝D-A，由方程(2)可以建立方程(1)的Jacobi迭代：

DX^(k+1)＝CX^(k)+F                    (3)

或

X^(k+1)＝D^-1CX^(k)+D^-1F               (4)

其中，上标(k)表示第k步迭代，方程(2)减去方程(3)并左乘D^-1得到第(k+1)步的误差为：

X-X^(k+1)＝E^(k+1)

＝D^-1CE^(k)                          (5)

＝(I-D^-1A)E^(k)

由式(5)可知，选择D＝D^*满足下面的条件时，第(k+1)步将得到精确解：

(I-D^*-1A)E^(k)＝0                    (6)

而迭代过程中第(k)步的残差定义为：

r^(k)＝F-AX^(k)                       (7)

将方程(1)代入到方程(7)得到：

r^(k)＝AE^(k)                         (8)

由式(8)可以得到E^(k)的近似计算为：

E^(k)＝A^-1r^(k)≈D^-1r^(k)              (9)

将式(9)代入方程(6)得到：

[I-D^*-1A]D^-1r^(k)＝0                 (10)

由方程(10)可以求得：

D^*＝diag(AD^-1r^(k)./D^-1r^(k))                (11)

其中符号./表示矩阵的元素除法，diag表示将向量扩展为对角矩阵，则根据式(4)和式(11)建立可变预处理迭代过程为：

X^(k+1)＝D^(k+1)-1C^(k+1)X^(k)+D^(k+1)-1F

＝X^(k)+D^(k+1)-1r^(k)                        (12)

D^(k+1)＝diag(AD^(k)-1r^(k)./D^(k)-1r^(k))

具体实施方式二：本实施方式在具体实施方式一的基础上进一步限定了步骤一和步骤二中所述的预处理矩阵通过对迭代过程的误差分析计算得到，其计算表达式为：

D^(k+1)＝diag(AD^(k)-1r^(k)./D^(k)-1r^(k))

式中A为系统的系数矩阵，r^(k)为第k步迭代的残差，D^(k)为第k步迭代的预处理矩阵，D^(k+1)为计算得到的第k+1步迭代的预处理矩阵，./表示矩阵的元素除法。

Claims (2)

1、基于可变预处理迭代获得非比例阻尼力学系统模态频响的方法，其特征在于它包括以下步骤：

步骤一、根据Jacobi迭代的预处理矩阵D⁽⁰⁾＝D，计算频率响应的初始解X⁽⁰⁾，并计算残差r⁽⁰⁾，令k＝0；

步骤二、计算第k+1步的预处理矩阵：

D^(k+1)＝diag(AD^(k)-1r^(k)./D^(k)-1r^(k))

根据第k+1步预处理矩阵计算第k+1步的频率响应向量：

X^(k+1)＝X^(k)+D^(k+1)-1r^(k)

并计算第k+1步的残差：

r^(k+1)＝F-AX^(k+1)

步骤三、判断第k+1步的残差是否满足残差范数的相对值小于给定的收敛指标ε，即 $\frac{\left|\right|r\left|\right|}{\left|\right|F\left|\right|}<\mathrm{\&epsiv;},$ 判断结果为是，则迭代计算结束，此时X^(k+1)即为得到的模态频响；判断结果为否，则令k＝k+1并返回步骤二继续计算，所述k是整数，且k≥0。

2、根据权利要求1所述的基于可变预处理迭代获得非比例阻尼力学系统模态频响的方法，其特征在于步骤一和步骤二中所述的预处理矩阵通过对迭代过程的误差分析计算得到，其计算表达式为：

D^(k+1)＝diag(AD^(k)-1r^(k)./D^(k)-1r^(k))

式中A为系统的系数矩阵，r^(k)为第k步迭代的残差，D^(k)为第k步迭代的预处理矩阵，D^(k+1)为计算得到的第k+1步迭代的预处理矩阵，./表示矩阵的元素除法。