第四章 結(jié)構(gòu)固有振動特征值問題的數(shù)值解

1、第四章 結(jié)構(gòu)固有振動特征值問題的數(shù)值解4.1 概述 根據(jù)結(jié)構(gòu)振動的數(shù)學(xué)模型，即振動微分方程所形成的矩陣特征值問題，求解結(jié)構(gòu)的固有振動特性固有頻率與固有振型，是結(jié)構(gòu)振動分析的一個主要任務(wù)。結(jié)構(gòu)的固有振動特性是結(jié)構(gòu)振動的內(nèi)因。固有振動特性也是進行結(jié)構(gòu)振動響應(yīng)分析和結(jié)構(gòu)動力學(xué)設(shè)計的基礎(chǔ)。對于簡單的結(jié)構(gòu)，如均勻直梁、均勻直桿等，可以用解析的方法解得其固有振動特性。對于一般結(jié)構(gòu)，如果只需獲得結(jié)構(gòu)有限階的固有振動特性，也可以采用試驗測試（模態(tài)識別）的方法來獲得。但是對于大型復(fù)雜結(jié)構(gòu)，不可能用解析分析的方法得到其固有振動特性，而采用試驗測試的方法不僅花費高，而且周期長，對于處于設(shè)計狀態(tài)的結(jié)構(gòu)，顯然也無法進行

2、試驗。所以對復(fù)雜的工程結(jié)構(gòu)，常用的方法是建立結(jié)構(gòu)的數(shù)學(xué)模型，用數(shù)值求解的方法獲得結(jié)構(gòu)的固有振動特性。隨著計算技術(shù)飛速發(fā)展和特征值計算方法的研究進展，通過矩陣特征值問題的求解來獲得結(jié)構(gòu)固有振動特性，是已經(jīng)被振動工程界普遍接受的一個有效和可靠的途徑。從數(shù)學(xué)理論上也可以證明，許多特征值計算方法具有相當(dāng)好的精度，并且獲得了實踐和實驗的證明。由于結(jié)構(gòu)固有振動特性求解與矩陣特征值求解問題的密切關(guān)系，在結(jié)構(gòu)振動分析中，矩陣特征值問題已經(jīng)成為結(jié)構(gòu)固有振動特性分析的一個代名詞。所以在本章中，只要不作說明，一般講的矩陣特征值問題就是指結(jié)構(gòu)的固有振動特性求解問題。所謂系統(tǒng)的特征值就代指結(jié)構(gòu)的固有頻率，特征向量代指結(jié)

3、構(gòu)的固有振型（固有模態(tài)）矩陣特征值問題的數(shù)值求解方法可以分為三類：矩陣分解法、迭代法和矩陣變換法。由于矩陣（代數(shù)）特征值問題本身就是一個完整的系統(tǒng)，本章只能根據(jù)結(jié)構(gòu)固有振動分析問題的需要，介紹一些常用的求解方法。詳盡的矩陣特征值問題的數(shù)值求解方法可以參考威爾金森的名著代數(shù)特征值問題。本章的論述是建立在已經(jīng)用有限元素法建立了結(jié)構(gòu)振動運動數(shù)學(xué)模型的基礎(chǔ)上。4.2 結(jié)構(gòu)振動特征值問題的性質(zhì)根據(jù)結(jié)構(gòu)振動方程，可以得到結(jié)構(gòu)固有振動的代數(shù)特征值問題： （41）或 （） （42）振動特征值問題除了第二章所述的性質(zhì)外，在特征值問題的數(shù)值求解中，還要用到如下一些性質(zhì)：1 移軸特性對特征值問題 （43）若為一已知

4、實數(shù)，則有： （44）新的特征值問題可寫為： （45） （） （46）顯然，上面兩個特征值問題具有相同的特征向量，而特征值間的關(guān)系為： （47）稱為移軸量。在結(jié)構(gòu)振動分析中，移軸特性常用來消除剛度矩陣的奇異性，也可用來加速迭代求解的收斂速度。2 特征值對合同變換的不變性合同變換：若為一個非奇異矩陣，則變換稱為對矩陣的合同變換。對矩陣特征值問題： （48）中的質(zhì)量陣和剛度陣，用矩陣作如下合同變換： （49） （410）則得到一個新的特征值問題 (411)從而有： (412)由于，故有 （413）顯然變換后的特征值不變：（414）且可以容易地證明，變換前后兩個特征值問題的特征向量之間具有關(guān)系： （

5、415）即合同變換不改變矩陣的特征值，特征向量具有轉(zhuǎn)換關(guān)系。3 特征值對相似變換的不變性相似變換：對非奇異陣，變換稱為相似變換。對矩陣特征值問題： （416）中的質(zhì)量陣和剛度陣，用矩陣作如下相似變換： （417） （418）則得到一個新的特征值問題： （419）從而有： （420）由于，故有 （421）顯然變換后的特征值不變： （422）且可以證明，變換前后的特征向量間具有關(guān)系： （423）即相似變換不改變矩陣的特征值，特征向量具有轉(zhuǎn)換關(guān)系。4 特征值對正交變換的不變性正交變換：在相似變換中，若對于非奇異矩陣，有，即，則用進行的相似變換稱為正交變換。正交變換不僅不改變矩陣的特征值，而且正交變換

6、后矩陣的對稱性不變。5 特征值對旋轉(zhuǎn)變換的不變性旋轉(zhuǎn)變換：取變換矩陣為： （424）可以證明為正交矩陣，用矩陣進行的相似變換，稱為旋轉(zhuǎn)變換。通過一系列的旋轉(zhuǎn)變換，可以使對稱陣變?yōu)閷顷?，使一般矩陣變?yōu)槿顷嚒＠蒙鲜鎏卣髦祮栴}的性質(zhì)，可以得到不同的特征值求解方法。結(jié)構(gòu)振動特征值的數(shù)值求解中常用的方法主要有：【多項式迭代法】多項式迭代法是利用特征值使特征多項式等于零的性質(zhì)，即 （425）【利用特征多項式的Sturm序列性質(zhì)的方法】通過Sturm序列的性質(zhì)，分離出每個特征值所在的區(qū)間，通過二分法和其它加速查找方法，逐步縮小特征值所在的分隔區(qū)間，最后得到所需的特征值?！臼噶康ā渴噶康ɡ玫幕?/p>

7、本關(guān)系式為： （426）具體的迭代方法有：正迭代法、逆迭代法、子空間迭代（同時迭代）法?！揪仃囎儞Q法】矩陣變換法利用的是“矩陣特征值對合同變換的不變性”這一性質(zhì)。利用具有某種性質(zhì)的變換矩陣，經(jīng)過一系列變換，得到一個形式簡單的、容易進行特征值求解的矩陣特征值問題。常用的有雅可比方法：它是經(jīng)過一系列旋轉(zhuǎn)變換，最后使一個對稱陣成為一個對角陣，然后很容易地求解該對角陣的特征值。QR法：QR法是經(jīng)過一系列的相似變換，將一個實數(shù)矩陣轉(zhuǎn)化為一個三角陣。然后求解其特征值。豪斯霍爾德方法：豪斯霍爾德方法 是通過正交相似變換，將一個實對稱矩陣化為三對角陣。然后用QR法求其特征值。還有其它一些方法。本章重點介紹特征

8、多項式迭代法、基于Sturm序列性質(zhì)的方法、矢量迭代法、子空間迭代法，對于矩陣變換方法僅對變換格式作簡略介紹。在求解結(jié)構(gòu)固有振動特性時，有時各種方法也可以聯(lián)合使用，來加快求解的速度。4.2 多項式迭代方法使特征多項式（427）等于零的根就是特征值問題（428）的特征值解。當(dāng)階數(shù)時，無法用公式求根，只能用數(shù)值計算方法求解上方程。而求解方程 （429）的方法通常有顯式和隱式兩類。并且一般都采用“頻率掃描”方法求解。【顯式迭代】將展開成的多項式，然后進行迭代求解，因而，如果陣階數(shù)很大時，將多項式展開就有一個缺點方程的系數(shù)計算誤差引起根的穩(wěn)定性很差，因而在實際結(jié)構(gòu)振動分析中很少采用。【隱式迭代】先給定

9、值，對矩陣（）進行高斯三角分解，再計算出的值。即做分解（430）和分別是單位下三角陣和上三角陣。從而： （431）實際計算時，是從的初值開始，依次計算， 的值，根據(jù)它們的正負號，確定出根所在的區(qū)間（隔根區(qū)間），然后對各隔根區(qū)間用二分法（或0.618優(yōu)選法）進行細分直到求出符合給定精度的解，也可以用更快的割線法來加速尋根：（432）但要注意，使用割線法時，迭代的初值很重要，一般與二分法聯(lián)合使用，以保證迭代的迅速收斂?！纠坑嬎憔仃囂卣髦祮栴}的第一個特征值。，【解】由于矩陣都是正定陣，故取迭代初值，。，迭代法的優(yōu)缺點：理論簡明，方法簡單，可以對任一個特征根進行迭代搜索，但計算量大，且容易漏根。4.

10、3 基于Sturm序列性質(zhì)的方法【Sturm序列】對廣義特征值問題 （433）其對應(yīng)的特征多項式為： （434）若將，中的最后行和最后列去掉，得到相應(yīng)的一個新廣義特征值問題： （435）稱為原系統(tǒng)的“第階伴隨約束系統(tǒng)的特征值問題”，這相當(dāng)于對原系統(tǒng)增加了個約束。定義第r階伴隨約束系統(tǒng)的特征多項式：取，得到個伴隨特征多項式： （436）這個特征多項式組成的序列稱為Sturm序列，根據(jù)瑞利約束原理，第階伴隨系統(tǒng)與第階伴隨系統(tǒng)特征值和之間有關(guān)系： （437）這一性質(zhì)又叫特征值隔離定理，在第一章已經(jīng)介紹過。00000【Sturm序列的性質(zhì)】上面所說的各階系統(tǒng)的特征值，就是Sturm序列中各個多項式的根

11、，把上述各個特征多項式的根所在區(qū)間的情況用曲線畫出來，如圖所示。用直線段把曲線上所有的根連接起來，得到一個折線簇。現(xiàn)考慮任意一個非特征值，并設(shè)，在曲線上作一條的直線，則這一直線就確定了在各伴隨系統(tǒng)的特征值關(guān)系中的位置。對應(yīng)于的直線必與折線簇相交，且根據(jù)特征值隔離定理，它一定不會與折線簇相交。根據(jù)(4-31)式，由高斯消元過程和伴隨特征多項式的定義，知： （438）所以，直線與折線簇相交一次，相應(yīng)的就在的對角線上出現(xiàn)一個負元素，故的對角線上正好有個負元素。因此，Sturm序列的性質(zhì)可以表述為：對角線上負元素的個數(shù)，正好等于小于的特征值的個數(shù)。利用這個性質(zhì)，我們可以來求解系統(tǒng)的特征值。【基于Stu

12、rm序列性質(zhì)的特征值求解】（1）取兩個任意的初值，求和，實際做法是用高斯消元法，直接進行消元，得到上三角陣。（2）由中對角線上負元素的個數(shù)，得到隔根區(qū)間中特征值的個數(shù)（3）若，取，求，根據(jù)對角元負數(shù)的個數(shù)，確定出區(qū)間和中特征值的個數(shù)為：和（4）按照與第三步相同的方法，對各個隔根區(qū)間進行二分，直到每個隔根區(qū)間內(nèi)的特征值個數(shù)為止。（5）對每個只含一個特征值的隔根區(qū)間，用和作為初值進行多項式迭代法（如割線迭代法），求出具體的特征值。從上述步驟可見，實際上只是利用Sturm序列的性質(zhì)，分離出各個特征值所在的區(qū)間，真正求解特征值還得用多項式迭代法?！纠坷肧turm序列性質(zhì)分離特征值問題的特征值。并用

13、割線迭代公式求出。其中 【解】：設(shè)，計算對角線上沒有負元素，所以設(shè)，計算可見，對角線上負元素有3個，所以，設(shè)，計算可見，對角線上負元素有2個，所以設(shè)，計算可見，對角線上負元素有1個，所以從而知道，取初值，進行割線迭代，代入公式（432）得故 4.4矢量迭代法矢量迭代法有正迭代和逆迭代兩種格式，矢量迭代法是逐個迭代出系統(tǒng)的 各個特征值。正迭代法最先得到的是系統(tǒng)的最大特征值，而逆迭代最先得到的是系統(tǒng)的最低階特征值。由于結(jié)構(gòu)振動分析中，關(guān)心的是結(jié)構(gòu)的低階固有頻率，并且由結(jié)構(gòu)有限元模型計算結(jié)構(gòu)高階特征值的精度一般較差，所以，通常是用逆迭代法求解系統(tǒng)的低階特征值。矢量逆迭代有幾種不同的具體格式，但基本思

14、想和迭代關(guān)系是相同的?！臼噶磕娴?基本公式 （439） （440）稱為動力矩陣2迭代具體步驟（i）選取迭代的初始矢量 （441）（ii）將歸一化，得到（442）為中某個特定的元素，為了避免出現(xiàn)零元素，通常取或。（當(dāng)然也可以用其他歸一化方法） （iii）求 （443） （444）根據(jù)展開定理： （445）按照（441）的格式重復(fù)次后，有： （446）當(dāng)時，有 （447）或 （448）（449）可以利用位移展開定理，證明上述結(jié)論。用矢量迭代法求系統(tǒng)第二階以上的特征值時，要先進行動力矩陣的“掃?！薄癝weep Mode”。所謂“掃?！保词侵笍倪x擇的初始向量（模態(tài)）中，抽?。ā皰摺保┏銮懊娴碗A模

15、態(tài)的成分。選初始向量，根據(jù)展開定理: （450）前乘動力矩陣，得到： （451）由固有振型的正交性，知 （452） （453）從中減去的分量，（454） （455）稱為掃模動力矩陣，可以證明，用陣對任意選定的初始向量進行逆迭代，得到的都是系統(tǒng)的第二階特征值和第二階特征向量。與上述推導(dǎo)過程相同，如果已經(jīng)求得了前階特征對，那么，求系統(tǒng)第階特征對的掃模動力矩陣為：（456）或?qū)懗蛇f推形式： （457）【矢量正迭代】矢量正迭代法與矢量逆迭代過程完全一樣，只是迭代開始時的基本方程為：選取任意初始向量進行迭代，可得到系統(tǒng)的最高階特征對?！娟P(guān)于矢量迭代法的討論】1 對于半正定系統(tǒng)，由于剛度矩陣奇異，不能直接

16、用逆迭代方法，如果要用矢量逆迭代法進行特征對的迭代求解，可以采用移軸技術(shù)，使剛度矩陣正定。2 由于在迭代過程中的計算數(shù)值誤差積累，迭代得到的高階特征對精度比低階特征對的精度差，因此，矩陣迭代法通常只用來求解系統(tǒng)的低階特征對。3 對重頻情況，迭代法仍然有效，只是選用不同的初始向量，得到相互線性無關(guān)的特征向量，因此，在實際迭代求解時，如果用不同的初始向量進行迭代，得到相同的特征值和線性無關(guān)的特征向量，則可以判定其為重頻，并通過正交化得到相互正交的特征矢量。4 由于數(shù)值計算誤差，進行掃模動力矩陣計算時，實際上是不能把前面的低階模態(tài)完全掃除干凈的。所以用掃模動力矩陣進行迭代得到的高階特征矢量與前面的低

17、階特征矢量的正交性會受到一定的影響，可以在求出所需的特征對后，用特征矢量的正交化方法（如許密特正交化方法）進行進一步的正交化處理。5矢量迭代法的優(yōu)點是原理清楚，方法簡單，易于程序?qū)崿F(xiàn)，對初始矢量的選擇沒有限制，但缺點是計算量大，且必須從第一階特征對開始逐對進行迭代。計算效率低，而且高階特征對的結(jié)果精度較差。6為了加速迭代的收斂速度，可以利用特征值的移軸定理，即在迭代的基本方程兩邊同時減去一個矢量： （458）因此，相當(dāng)于將動力矩陣的對角線元素都減去一個常數(shù)，得到的新特征值問題的特征值為，而特征矢量不變。記，按照（446）式，有：（459）顯然，特征值迭代的收斂速度取決于比值： （460）因為（

18、461）所以，對于滿足條件的移軸量，必有： （462）考慮到所有高階項的影響，一般取 （463）可以取得很好的加速迭代效果。即在同樣的迭代次數(shù)下用（459）式進行迭代，比用（446）式迭代的收斂速度快。4.5 子空間迭代法（矢量同時迭代法）由矢量迭代法的基本原理可知，如果開始選用一組彼此線性無關(guān)的初始向量進行迭代，只要它們與第一階特征矢量不正交，而且在迭代中不進行正交化處理，那么迭代的結(jié)果是它們都會收斂于第一階特征矢量。反之，如果在迭代過程中，不斷的進行正交化處理，則迭代最后可以得到若干低階特征對。這就是子空間迭代法的基本思想。由于它是同時對幾個特征向量進行迭代求解，因此，又叫做同時迭代法。在

19、子空間迭代中，一般采用的正交化方法是李茲法。因此可以認(rèn)為矢量逆迭代法和李茲法的結(jié)合就是子空間迭代法。對系統(tǒng)的若干階特征對，它們滿足的特征方程寫成矩陣形式為： （464）或 （465）其中， （466） （467）引入動力矩陣 （468） （469）子空間迭代的步驟如下：（1） 選擇個任意矢量組成初始矢量矩陣，要求組成該矩陣的各矢量與所要求解的前階特征矢量不正交，而且相互獨立。記為。通常取為： （470）為了保證前階特征對的精度，推薦的選擇按如下規(guī)則確定 （471）（2） 進行迭代運算 （472）在實際計算求解時，為了減少數(shù)值誤差，不是用剛度矩陣求逆在乘以質(zhì)量陣，來得到動力矩陣，然后計算上式。而

20、是用高斯消元法求解下面形式的矩陣代數(shù)方程： （473）得到 （474）（3） 進行李茲計算以為李茲基底（假設(shè)模態(tài)），進行李茲坐標(biāo)變換，即取 （475）得到一個新的、降階后的特征值問題： （476）（477）（478）（4） 對解得的進行歸一化處理，通常采用模態(tài)質(zhì)量為1的方法：（479）（5） 形成下一次迭代的矢量矩陣（480）重復(fù)（2）（5）步驟，直到求出的系統(tǒng)前個特征值都滿足精度要求為止?？梢宰C明，經(jīng)過足夠多次數(shù)的迭代后，有：（481）迭代過程中得到的對于是加權(quán)正交的，因而對M和K滿足加權(quán)正交。由于在迭代過程中不斷進行正交化處理，故最后得到的特征矢量是相互正交的。4.6 Lanczos 方法

21、Lanczos方法是目前求解大型特征值問題的最有效的方法，其計算量比子空間迭代法要少數(shù)倍。Lanczos方法產(chǎn)生于1950年代，其基本思想是用遞推公式產(chǎn)生一個正交的矢量矩陣Lanczos矢量矩陣，然后通過Lanczos矢量矩陣將原來對稱矩陣的特征值問題變換成一個三對角矩陣的特征值問題。用Lanczos方法求解特征值問題，包括兩個步驟：Lanczos矢量矩陣的形成和三對角矩陣的特征值求解。Lanczos方法求解特征值問題的優(yōu)點在于，僅僅通過矩陣的相乘運算即可獲得一個對于系統(tǒng)真實模態(tài)空間的、品質(zhì)優(yōu)良的假設(shè)模態(tài)矩陣（稱為截斷的Lanczos矢量矩陣），它所張成的模態(tài)空間，可以很好的覆蓋結(jié)構(gòu)的離散化模

22、型的低維真實模態(tài)空間。因此，利用截斷的Lanczos矢量矩陣，不僅可以使系統(tǒng)模型降維，而且可以用來求解系統(tǒng)低階模態(tài)特性，還可以用來求解系統(tǒng)的強迫振動響應(yīng)。【Lanczos矢量矩陣】分別記結(jié)構(gòu)的質(zhì)量矩陣和剛度矩陣為和，首先選擇一個初始矢量，并進行歸一化，使： （481）然后按如下遞推公式，求出第個Lanczos矢量： （482）其中 （483）若迭代到，則算法精確收斂，停止迭代計算?？紤]到編程的方便，在實際計算時，采用如下的迭代格式：選擇適當(dāng)?shù)某踔凳噶?，并使計算，令，對作如下計算：?） （484）（2） （485）（3） （486）（4） （487）（5） （488）在時，就完成了第（1）步，

23、即求出，然后就停止迭代。由于計算誤差，Lanczos矢量間的正交性將會隨著運算次數(shù)的增加而變壞，需要隨時進行檢查，并進行正交化處理。【Gram-Schmidt正交化方法】定義一個衡量正交性的量：（489）令為一個設(shè)定的小量，若，則需對進行正交化處理。令第個Lanczos矢量正交化后為（488）由正交性條件來確定，即（488）【截斷的Lanczos矢量矩陣】經(jīng)過上述的正交化處理后，就可得到截斷的Lanczos矢量矩陣： （489）按照假設(shè)模態(tài)法，得到（489）式的假設(shè)模態(tài)陣后，假設(shè)模態(tài)空間中的廣義質(zhì)量陣和廣義剛度陣為： （490）注意：Lanczos方法并不是直接去解方程（490），而是用Lan

24、czos矢量矩陣對系統(tǒng)的物理坐標(biāo)進行坐標(biāo)變換：（491）代入自由振動方程進行變換，化簡得到：（492）（493）用左乘上式：（494）注意到在形成時Lanczos矢量已經(jīng)對質(zhì)量陣歸一化了，所以上式寫成：（495）其中 （496）可以證明，為三對角矩陣，且： （497）從而，系統(tǒng)的廣義特征值問題，就化為三對角陣的標(biāo)準(zhǔn)特征值問題：（498）解得上述特征值問題的全部特征對，就得到原系統(tǒng)的階低階特征對的近似值。具體求解可利用結(jié)合Sturm序列性質(zhì)的二分法或QR法求解?！緦嵗龑Ρ取?對一個600階的特征值問題，在Siemens7760計算機上計算，用子空間迭代法和Lanczos方法進行求解得到前21階特征值，兩種方法的 CPU時間分別為535.11秒（子空間迭代法）和116.56秒（Lanczos方法）。4.7 矩陣變換法（1） 將廣義特征值問題變?yōu)闃?biāo)準(zhǔn)特征值問題對廣義特征值問題（499）如果采用如下變換過程：（4100）就得到一個標(biāo)準(zhǔn)特征