第四章 结构固有振动特征值问题的数值解

1、第四章 结构固有振动特征值问题的数值解4.1 概述 根据结构振动的数学模型，即振动微分方程所形成的矩阵特征值问题，求解结构的固有振动特性固有频率与固有振型，是结构振动分析的一个主要任务。结构的固有振动特性是结构振动的内因。固有振动特性也是进行结构振动响应分析和结构动力学设计的基础。对于简单的结构，如均匀直梁、均匀直杆等，可以用解析的方法解得其固有振动特性。对于一般结构，如果只需获得结构有限阶的固有振动特性，也可以采用试验测试（模态识别）的方法来获得。但是对于大型复杂结构，不可能用解析分析的方法得到其固有振动特性，而采用试验测试的方法不仅花费高，而且周期长，对于处于设计状态的结构，显然也无法进行

2、试验。所以对复杂的工程结构，常用的方法是建立结构的数学模型，用数值求解的方法获得结构的固有振动特性。随着计算技术飞速发展和特征值计算方法的研究进展，通过矩阵特征值问题的求解来获得结构固有振动特性，是已经被振动工程界普遍接受的一个有效和可靠的途径。从数学理论上也可以证明，许多特征值计算方法具有相当好的精度，并且获得了实践和实验的证明。由于结构固有振动特性求解与矩阵特征值求解问题的密切关系，在结构振动分析中，矩阵特征值问题已经成为结构固有振动特性分析的一个代名词。所以在本章中，只要不作说明，一般讲的矩阵特征值问题就是指结构的固有振动特性求解问题。所谓系统的特征值就代指结构的固有频率，特征向量代指结

3、构的固有振型（固有模态）矩阵特征值问题的数值求解方法可以分为三类：矩阵分解法、迭代法和矩阵变换法。由于矩阵（代数）特征值问题本身就是一个完整的系统，本章只能根据结构固有振动分析问题的需要，介绍一些常用的求解方法。详尽的矩阵特征值问题的数值求解方法可以参考威尔金森的名著代数特征值问题。本章的论述是建立在已经用有限元素法建立了结构振动运动数学模型的基础上。4.2 结构振动特征值问题的性质根据结构振动方程，可以得到结构固有振动的代数特征值问题： （41）或 （） （42）振动特征值问题除了第二章所述的性质外，在特征值问题的数值求解中，还要用到如下一些性质：1 移轴特性对特征值问题 （43）若为一已知

4、实数，则有： （44）新的特征值问题可写为： （45） （） （46）显然，上面两个特征值问题具有相同的特征向量，而特征值间的关系为： （47）称为移轴量。在结构振动分析中，移轴特性常用来消除刚度矩阵的奇异性，也可用来加速迭代求解的收敛速度。2 特征值对合同变换的不变性合同变换：若为一个非奇异矩阵，则变换称为对矩阵的合同变换。对矩阵特征值问题： （48）中的质量阵和刚度阵，用矩阵作如下合同变换： （49） （410）则得到一个新的特征值问题 (411)从而有： (412)由于，故有 （413）显然变换后的特征值不变：（414）且可以容易地证明，变换前后两个特征值问题的特征向量之间具有关系： （

5、415）即合同变换不改变矩阵的特征值，特征向量具有转换关系。3 特征值对相似变换的不变性相似变换：对非奇异阵，变换称为相似变换。对矩阵特征值问题： （416）中的质量阵和刚度阵，用矩阵作如下相似变换： （417） （418）则得到一个新的特征值问题： （419）从而有： （420）由于，故有 （421）显然变换后的特征值不变： （422）且可以证明，变换前后的特征向量间具有关系： （423）即相似变换不改变矩阵的特征值，特征向量具有转换关系。4 特征值对正交变换的不变性正交变换：在相似变换中，若对于非奇异矩阵，有，即，则用进行的相似变换称为正交变换。正交变换不仅不改变矩阵的特征值，而且正交变换

6、后矩阵的对称性不变。5 特征值对旋转变换的不变性旋转变换：取变换矩阵为： （424）可以证明为正交矩阵，用矩阵进行的相似变换，称为旋转变换。通过一系列的旋转变换，可以使对称阵变为对角阵，使一般矩阵变为三角阵。利用上述特征值问题的性质，可以得到不同的特征值求解方法。结构振动特征值的数值求解中常用的方法主要有：【多项式迭代法】多项式迭代法是利用特征值使特征多项式等于零的性质，即 （425）【利用特征多项式的Sturm序列性质的方法】通过Sturm序列的性质，分离出每个特征值所在的区间，通过二分法和其它加速查找方法，逐步缩小特征值所在的分隔区间，最后得到所需的特征值。【矢量迭代法】矢量迭代法利用的基

7、本关系式为： （426）具体的迭代方法有：正迭代法、逆迭代法、子空间迭代（同时迭代）法。【矩阵变换法】矩阵变换法利用的是“矩阵特征值对合同变换的不变性”这一性质。利用具有某种性质的变换矩阵，经过一系列变换，得到一个形式简单的、容易进行特征值求解的矩阵特征值问题。常用的有雅可比方法：它是经过一系列旋转变换，最后使一个对称阵成为一个对角阵，然后很容易地求解该对角阵的特征值。QR法：QR法是经过一系列的相似变换，将一个实数矩阵转化为一个三角阵。然后求解其特征值。豪斯霍尔德方法：豪斯霍尔德方法 是通过正交相似变换，将一个实对称矩阵化为三对角阵。然后用QR法求其特征值。还有其它一些方法。本章重点介绍特征

8、多项式迭代法、基于Sturm序列性质的方法、矢量迭代法、子空间迭代法，对于矩阵变换方法仅对变换格式作简略介绍。在求解结构固有振动特性时，有时各种方法也可以联合使用，来加快求解的速度。4.2 多项式迭代方法使特征多项式（427）等于零的根就是特征值问题（428）的特征值解。当阶数时，无法用公式求根，只能用数值计算方法求解上方程。而求解方程 （429）的方法通常有显式和隐式两类。并且一般都采用“频率扫描”方法求解。【显式迭代】将展开成的多项式，然后进行迭代求解，因而，如果阵阶数很大时，将多项式展开就有一个缺点方程的系数计算误差引起根的稳定性很差，因而在实际结构振动分析中很少采用。【隐式迭代】先给定

9、值，对矩阵（）进行高斯三角分解，再计算出的值。即做分解（430）和分别是单位下三角阵和上三角阵。从而： （431）实际计算时，是从的初值开始，依次计算， 的值，根据它们的正负号，确定出根所在的区间（隔根区间），然后对各隔根区间用二分法（或0.618优选法）进行细分直到求出符合给定精度的解，也可以用更快的割线法来加速寻根：（432）但要注意，使用割线法时，迭代的初值很重要，一般与二分法联合使用，以保证迭代的迅速收敛。【例】计算矩阵特征值问题的第一个特征值。，【解】由于矩阵都是正定阵，故取迭代初值，。，迭代法的优缺点：理论简明，方法简单，可以对任一个特征根进行迭代搜索，但计算量大，且容易漏根。4.

10、3 基于Sturm序列性质的方法【Sturm序列】对广义特征值问题 （433）其对应的特征多项式为： （434）若将，中的最后行和最后列去掉，得到相应的一个新广义特征值问题： （435）称为原系统的“第阶伴随约束系统的特征值问题”，这相当于对原系统增加了个约束。定义第r阶伴随约束系统的特征多项式：取，得到个伴随特征多项式： （436）这个特征多项式组成的序列称为Sturm序列，根据瑞利约束原理，第阶伴随系统与第阶伴随系统特征值和之间有关系： （437）这一性质又叫特征值隔离定理，在第一章已经介绍过。00000【Sturm序列的性质】上面所说的各阶系统的特征值，就是Sturm序列中各个多项式的根

11、，把上述各个特征多项式的根所在区间的情况用曲线画出来，如图所示。用直线段把曲线上所有的根连接起来，得到一个折线簇。现考虑任意一个非特征值，并设，在曲线上作一条的直线，则这一直线就确定了在各伴随系统的特征值关系中的位置。对应于的直线必与折线簇相交，且根据特征值隔离定理，它一定不会与折线簇相交。根据(4-31)式，由高斯消元过程和伴随特征多项式的定义，知： （438）所以，直线与折线簇相交一次，相应的就在的对角线上出现一个负元素，故的对角线上正好有个负元素。因此，Sturm序列的性质可以表述为：对角线上负元素的个数，正好等于小于的特征值的个数。利用这个性质，我们可以来求解系统的特征值。【基于Stu

12、rm序列性质的特征值求解】（1）取两个任意的初值，求和，实际做法是用高斯消元法，直接进行消元，得到上三角阵。（2）由中对角线上负元素的个数，得到隔根区间中特征值的个数（3）若，取，求，根据对角元负数的个数，确定出区间和中特征值的个数为：和（4）按照与第三步相同的方法，对各个隔根区间进行二分，直到每个隔根区间内的特征值个数为止。（5）对每个只含一个特征值的隔根区间，用和作为初值进行多项式迭代法（如割线迭代法），求出具体的特征值。从上述步骤可见，实际上只是利用Sturm序列的性质，分离出各个特征值所在的区间，真正求解特征值还得用多项式迭代法。【例】利用Sturm序列性质分离特征值问题的特征值。并用

13、割线迭代公式求出。其中 【解】：设，计算对角线上没有负元素，所以设，计算可见，对角线上负元素有3个，所以，设，计算可见，对角线上负元素有2个，所以设，计算可见，对角线上负元素有1个，所以从而知道，取初值，进行割线迭代，代入公式（432）得故 4.4矢量迭代法矢量迭代法有正迭代和逆迭代两种格式，矢量迭代法是逐个迭代出系统的 各个特征值。正迭代法最先得到的是系统的最大特征值，而逆迭代最先得到的是系统的最低阶特征值。由于结构振动分析中，关心的是结构的低阶固有频率，并且由结构有限元模型计算结构高阶特征值的精度一般较差，所以，通常是用逆迭代法求解系统的低阶特征值。矢量逆迭代有几种不同的具体格式，但基本思

14、想和迭代关系是相同的。【矢量逆迭代】1基本公式 （439） （440）称为动力矩阵2迭代具体步骤（i）选取迭代的初始矢量 （441）（ii）将归一化，得到（442）为中某个特定的元素，为了避免出现零元素，通常取或。（当然也可以用其他归一化方法） （iii）求 （443） （444）根据展开定理： （445）按照（441）的格式重复次后，有： （446）当时，有 （447）或 （448）（449）可以利用位移展开定理，证明上述结论。用矢量迭代法求系统第二阶以上的特征值时，要先进行动力矩阵的“扫模”“Sweep Mode”。所谓“扫模”，即是指从选择的初始向量（模态）中，抽取（“扫”）出前面低阶模

15、态的成分。选初始向量，根据展开定理: （450）前乘动力矩阵，得到： （451）由固有振型的正交性，知 （452） （453）从中减去的分量，（454） （455）称为扫模动力矩阵，可以证明，用阵对任意选定的初始向量进行逆迭代，得到的都是系统的第二阶特征值和第二阶特征向量。与上述推导过程相同，如果已经求得了前阶特征对，那么，求系统第阶特征对的扫模动力矩阵为：（456）或写成递推形式： （457）【矢量正迭代】矢量正迭代法与矢量逆迭代过程完全一样，只是迭代开始时的基本方程为：选取任意初始向量进行迭代，可得到系统的最高阶特征对。【关于矢量迭代法的讨论】1 对于半正定系统，由于刚度矩阵奇异，不能直接

16、用逆迭代方法，如果要用矢量逆迭代法进行特征对的迭代求解，可以采用移轴技术，使刚度矩阵正定。2 由于在迭代过程中的计算数值误差积累，迭代得到的高阶特征对精度比低阶特征对的精度差，因此，矩阵迭代法通常只用来求解系统的低阶特征对。3 对重频情况，迭代法仍然有效，只是选用不同的初始向量，得到相互线性无关的特征向量，因此，在实际迭代求解时，如果用不同的初始向量进行迭代，得到相同的特征值和线性无关的特征向量，则可以判定其为重频，并通过正交化得到相互正交的特征矢量。4 由于数值计算误差，进行扫模动力矩阵计算时，实际上是不能把前面的低阶模态完全扫除干净的。所以用扫模动力矩阵进行迭代得到的高阶特征矢量与前面的低

17、阶特征矢量的正交性会受到一定的影响，可以在求出所需的特征对后，用特征矢量的正交化方法（如许密特正交化方法）进行进一步的正交化处理。5矢量迭代法的优点是原理清楚，方法简单，易于程序实现，对初始矢量的选择没有限制，但缺点是计算量大，且必须从第一阶特征对开始逐对进行迭代。计算效率低，而且高阶特征对的结果精度较差。6为了加速迭代的收敛速度，可以利用特征值的移轴定理，即在迭代的基本方程两边同时减去一个矢量： （458）因此，相当于将动力矩阵的对角线元素都减去一个常数，得到的新特征值问题的特征值为，而特征矢量不变。记，按照（446）式，有：（459）显然，特征值迭代的收敛速度取决于比值： （460）因为（

18、461）所以，对于满足条件的移轴量，必有： （462）考虑到所有高阶项的影响，一般取 （463）可以取得很好的加速迭代效果。即在同样的迭代次数下用（459）式进行迭代，比用（446）式迭代的收敛速度快。4.5 子空间迭代法（矢量同时迭代法）由矢量迭代法的基本原理可知，如果开始选用一组彼此线性无关的初始向量进行迭代，只要它们与第一阶特征矢量不正交，而且在迭代中不进行正交化处理，那么迭代的结果是它们都会收敛于第一阶特征矢量。反之，如果在迭代过程中，不断的进行正交化处理，则迭代最后可以得到若干低阶特征对。这就是子空间迭代法的基本思想。由于它是同时对几个特征向量进行迭代求解，因此，又叫做同时迭代法。在

19、子空间迭代中，一般采用的正交化方法是李兹法。因此可以认为矢量逆迭代法和李兹法的结合就是子空间迭代法。对系统的若干阶特征对，它们满足的特征方程写成矩阵形式为： （464）或 （465）其中， （466） （467）引入动力矩阵 （468） （469）子空间迭代的步骤如下：（1） 选择个任意矢量组成初始矢量矩阵，要求组成该矩阵的各矢量与所要求解的前阶特征矢量不正交，而且相互独立。记为。通常取为： （470）为了保证前阶特征对的精度，推荐的选择按如下规则确定 （471）（2） 进行迭代运算 （472）在实际计算求解时，为了减少数值误差，不是用刚度矩阵求逆在乘以质量阵，来得到动力矩阵，然后计算上式。而

20、是用高斯消元法求解下面形式的矩阵代数方程： （473）得到 （474）（3） 进行李兹计算以为李兹基底（假设模态），进行李兹坐标变换，即取 （475）得到一个新的、降阶后的特征值问题： （476）（477）（478）（4） 对解得的进行归一化处理，通常采用模态质量为1的方法：（479）（5） 形成下一次迭代的矢量矩阵（480）重复（2）（5）步骤，直到求出的系统前个特征值都满足精度要求为止。可以证明，经过足够多次数的迭代后，有：（481）迭代过程中得到的对于是加权正交的，因而对M和K满足加权正交。由于在迭代过程中不断进行正交化处理，故最后得到的特征矢量是相互正交的。4.6 Lanczos 方法

21、Lanczos方法是目前求解大型特征值问题的最有效的方法，其计算量比子空间迭代法要少数倍。Lanczos方法产生于1950年代，其基本思想是用递推公式产生一个正交的矢量矩阵Lanczos矢量矩阵，然后通过Lanczos矢量矩阵将原来对称矩阵的特征值问题变换成一个三对角矩阵的特征值问题。用Lanczos方法求解特征值问题，包括两个步骤：Lanczos矢量矩阵的形成和三对角矩阵的特征值求解。Lanczos方法求解特征值问题的优点在于，仅仅通过矩阵的相乘运算即可获得一个对于系统真实模态空间的、品质优良的假设模态矩阵（称为截断的Lanczos矢量矩阵），它所张成的模态空间，可以很好的覆盖结构的离散化模

22、型的低维真实模态空间。因此，利用截断的Lanczos矢量矩阵，不仅可以使系统模型降维，而且可以用来求解系统低阶模态特性，还可以用来求解系统的强迫振动响应。【Lanczos矢量矩阵】分别记结构的质量矩阵和刚度矩阵为和，首先选择一个初始矢量，并进行归一化，使： （481）然后按如下递推公式，求出第个Lanczos矢量： （482）其中 （483）若迭代到，则算法精确收敛，停止迭代计算。考虑到编程的方便，在实际计算时，采用如下的迭代格式：选择适当的初值矢量，并使计算，令，对作如下计算：（1） （484）（2） （485）（3） （486）（4） （487）（5） （488）在时，就完成了第（1）步，

23、即求出，然后就停止迭代。由于计算误差，Lanczos矢量间的正交性将会随着运算次数的增加而变坏，需要随时进行检查，并进行正交化处理。【Gram-Schmidt正交化方法】定义一个衡量正交性的量：（489）令为一个设定的小量，若，则需对进行正交化处理。令第个Lanczos矢量正交化后为（488）由正交性条件来确定，即（488）【截断的Lanczos矢量矩阵】经过上述的正交化处理后，就可得到截断的Lanczos矢量矩阵： （489）按照假设模态法，得到（489）式的假设模态阵后，假设模态空间中的广义质量阵和广义刚度阵为： （490）注意：Lanczos方法并不是直接去解方程（490），而是用Lan

24、czos矢量矩阵对系统的物理坐标进行坐标变换：（491）代入自由振动方程进行变换，化简得到：（492）（493）用左乘上式：（494）注意到在形成时Lanczos矢量已经对质量阵归一化了，所以上式写成：（495）其中 （496）可以证明，为三对角矩阵，且： （497）从而，系统的广义特征值问题，就化为三对角阵的标准特征值问题：（498）解得上述特征值问题的全部特征对，就得到原系统的阶低阶特征对的近似值。具体求解可利用结合Sturm序列性质的二分法或QR法求解。【实例对比】 对一个600阶的特征值问题，在Siemens7760计算机上计算，用子空间迭代法和Lanczos方法进行求解得到前21阶特征值，两种方法的 CPU时间分别为535.11秒（子空间迭代法）和116.56秒（Lanczos方法）。4.7 矩阵变换法（1） 将广义特征值问题变为标准特征值问题对广义特征值问题（499）如果采用如下变换过程：（4100）就得到一个标准特征