張 靜，劉明輝，鄭鋼鐵，2

(1.哈爾濱工業(yè)大學 航天學院，哈爾濱 150001;2.清華大學 宇航學院，北京 100084)

Lanczos-QR方法在大型非比例阻尼結構復模態(tài)計算中的應用

張 靜1，劉明輝1，鄭鋼鐵1，2

(1.哈爾濱工業(yè)大學 航天學院，哈爾濱 150001;2.清華大學 宇航學院，北京 100084)

根據(jù)大型非比例阻尼結構的動力學模型特點，發(fā)展了適用于大型復特征值問題求解的Lanczos－QR方法。針對大規(guī)模非對稱實矩陣的標準特征值問題，此方法首先采用Lanczos迭代方法對矩陣進行降階;然后采用帶移位的雙重步位移QR方法求解降階后的實三對角矩陣的全部特征值和特征向量，最后經(jīng)過相應的變換得到原矩陣的特征解。另外，提出了分塊矩陣三角分解方法將廣義特征值進行標準化，該方既提高了計算效率，又避免了范數(shù)誤差。還對奇異問題的處理進行了討論。最后在FORTRAN語言中進行了程序實現(xiàn)，完成初步軟件化，并且通過算例對算法及其程序的精度和收斂速度進行了驗證。

特征值問題;非比例阻尼;Lanczos算法;QR算法

對于大型結構，除了建立精確的有限元模型［1］外，動力學計算方法的精確性也是保證結構動力分析精度的關鍵。結構模態(tài)分析又是各種動力分析的基礎，其計算結果直接影響結構動力分析的精度［2］。常用求解特征值的方法有很多種:子空間迭代法［3］、乘冪(逆冪)法、Lanczos方法［4－5］、QR 方法，迭代Ritz向量法、Jacobi方法等，然而，對于大型非比例阻尼系統(tǒng)，上述方法均存在局限性。本文討論系數(shù)矩陣具有如下特點的大型非比例阻尼結構的特征值問題:大型、稀疏、非對稱、正定或半正定。這類大規(guī)模非對稱實矩陣的標準特征值問題，解為復特征值和復特征向量，且只需要求解前若干階較低的特征值和模態(tài)向量。

針對以上要求以及各種算法的特點，本文對求解大型非比例阻尼結構復模態(tài)問題的Lanczos-QR迭代法進行了發(fā)展。此方法結合了Lanczos算法快速高效和QR方法穩(wěn)定可靠的優(yōu)點，首先采用Lanczos迭代方法對矩陣進行降階;再采用帶移位的雙重步位移QR方法求解降階后的特征問題。本文還提出使用分塊三角分解方法得到標準特征值問題，該方法在廣義二次特征問題的標準化過程中具有明顯的優(yōu)越性，不僅減小了矩陣求逆的規(guī)模，并且避免了由于矩陣元素差異大引起的范數(shù)誤差。另外本文還對剛度矩陣奇異時出現(xiàn)的問題進行了討論和完善，從而得到標準的計算流程，并設置了與通用有限元軟件NASTRAN的接口。最后通過算例驗證了本文中的算法及程序的可實現(xiàn)性，結果表明本算法具有較強的收斂性，且穩(wěn)定、可靠。

1 非比例阻尼結構的特征值問題

含有非比例阻尼系統(tǒng)的動力學方程為:

對于一般結構，矩陣M對稱正定，矩陣K對稱非負定。令y={x}T，將其化成等價的狀態(tài)空間形式:

經(jīng)過變換，得到廣義特征值問題:

對于式(3)所示的非對稱廣義特征值問題，若矩陣B非奇異，則可化為如下的標準特征值問題的形式:

2 非對稱矩陣特征問題的Lanczos算法

n階非對稱矩陣A的特征值問題(如式(4)所示)的Lanczos算法［6］的主要思想是利用兩組雙正交向量將矩陣A逐次三對角化，得到一個ns(ns≤n)階的三對角矩陣T，將矩陣T的特征值作為原矩陣A的ns個極大特征值的近似，將矩陣T的特征向量通過變換作為原矩陣A相應特征向量的近似。該方法最大的優(yōu)勢是大大減少了迭代次數(shù)，降低了存儲代價。

為了求出矩陣A的前m個極小特征值，可做如下變換:的特征問題轉化成求解如下的低階特征問題:

經(jīng)過以上過程，矩陣

首先設定一小量ε，計算Gram－Schmidt系數(shù):

3 具有原點位移的雙重步QR算法

具有原點位移的QR算法通過選取近似特征值作為位因子sk來提高QR算法的收斂速度。對于含復特征值的系統(tǒng)，采用雙重步QR算法可以避免復數(shù)運算［7］，雙重步QR算法的移位因子sk，sk+1可以通過求解系數(shù)矩陣T的右下角2×2階子陣的特征值來獲得，即:

由QR算法的性質，對于三對角矩陣T，迭代過程中產(chǎn)生的所有T(k)均為三對角矩陣［8］。迭代結束判定準則為出現(xiàn)次對角線元素 tl，l－1＜ ε(ε(tl，l+tl－1，l－1)(ε為設置的精度控制參數(shù))。

(1)當l=m(m為T的階數(shù))時，T的一個高階特征值可以由其最后一個主對角線元素獲得;

(2)當l=m－1時，T的兩個高階特征值可以由其2階后主子獲得。

當求出若干個特征值后，矩陣收縮相應的階數(shù)，迭代重新開始。通常迭代10次后還未確定出一個特征值時，可通過改變位移量來調整求解，一種可行的調整策略為:最后求得降階后的矩陣T的所有特征值，i=1，2，…，ns。代入式(15)，即可得到矩陣T的所有特征向量 φt。

4 復特征向量的變換

5 標準特征值轉化的分塊Cholesky分解方法

對于標準特征值問題:

其一，若剛度矩陣K奇異，則無法進行求逆;

其二，直接求逆運算得到的矩陣嚴重不對稱，會給數(shù)值計算中帶來很大的舍入誤差;

其三，對于大型矩陣，直接求逆運算的計算效率較低。

為了改善求解過程中的數(shù)值計算精度和效率，可以采用分塊Cholesky分解來求出標準特征值問題的系數(shù)矩陣。

首先，為了避免矩陣K奇異，可引入移位因子α，得到非奇異矩陣～B:

該方法只需要求兩次Cholesky分解，一次下三角矩陣求逆，三次矩陣乘法即可得到矩陣A。對于奇異問題也只需要多進行兩次矩陣減法。另一方面，該方法形成的系數(shù)矩陣近似反對稱，便于進行壓縮存儲，也不會在計算過程中引入很大的范數(shù)誤差，不需要再對系數(shù)矩陣進行平衡化處理，相當于減少了一個相似變換的過程。

6 Lanczos-QR算法實現(xiàn)

Lanczos-QR算法的主要思想是利用求解絕對值最小的特征值的逆冪法思想來改進截斷格式的Lanczos方法降低矩陣階數(shù)，得到一個包含原矩陣的低階特征值的三對角矩陣，再使用帶移位的雙重步QR方法對降階后的三對角矩陣進行全部特征值求解，以避免計算過程中的復數(shù)運算。

本文對算法進行了編程實現(xiàn)，程序采用科學計算程序FORTRAN語言進行實現(xiàn)。在程序設計中，遵循了以下的設計思想:① 設計中用子程序來分別實現(xiàn)各個功能，用主程序控制流程，這樣使得算法思路很清晰，易于進一步的改進、優(yōu)化;② 在程序中，設置了與現(xiàn)有有限元軟件的數(shù)據(jù)接口，方便對其進行讀取。

使用改進截斷格式的Lanczos算法求解大型復特征值問題軟件實現(xiàn)流程如下。

圖1 程序設計流程圖Fig.1 Programming flowchart

7 數(shù)值算例

7.1 約束阻尼梁

為了驗證本文提出的算法在復特征值求解時的有效性和準確性，采用含粘彈性阻尼的約束阻尼結構來進行計算。如圖2所示的約束阻尼懸臂梁，阻尼材料選用ZN－1型粘彈性材料。具體參數(shù)如下:

首先在MATLAB軟件中建立該問題的有限元模型，單元類型及離散化方程可參見文獻［10］。粘彈性阻尼等效使用文獻［11］中所提到的Biot模型進行建模，經(jīng)過等效之后得到復特征值問題。分別用MATLAB和本文算法求解，計算結果如表1所示。由表1中數(shù)據(jù)可見，與MATLAB結果相比，本文程序的計算結果具有很高的精度。

圖2 約束阻尼梁Fig.2 A Sandwich Beam with Constrained Layer Damping

表1 頻率和阻尼比計算結果Tab.1 Frequency and damping ratio results

7.2 大型衛(wèi)星結構

該算例的目的是驗證本文算法在求解大型結構時的精度。計算用衛(wèi)星整體有限元模型(見圖3)在大型有限元軟件MSC/NASTRAN中建立。計算模型包括:節(jié)點:7162個;單元:9 964個;自由度:36 289個。具體材料參數(shù)與單元特性參數(shù)均由工程單位提供。

對底部進行約束，分別使用MSC/NASTRAN軟件和本文程序計算得到整星前10階固有頻率，并與工程單位提供的衛(wèi)星整星模型前四階模態(tài)測試結果數(shù)據(jù)對比，結果如表2所示。結果表明本文程序在計算大型結構特征值問題中也是有效的。

圖3 衛(wèi)星有限元模型Fig.3 Satellite FEM model

表2 頻率計算結果Tab.2 Frequency results

8 結論

針對含有非比例阻尼的大型結構動力學分析中出現(xiàn)的大規(guī)模復特征值的求解問題，本文發(fā)展了適用于大規(guī)模非對稱實矩陣的標準特征值問題的求解方法，該方法集合了Lanczos算法的高效率和QR算法的穩(wěn)定性，可以有效解決大規(guī)模復特征值問題的求解。本文算法具有如下優(yōu)點:① 采用分塊矩陣的三角分解來對廣義特征值問題進行標準化，可以避免大型矩陣求逆運算和非對稱矩陣引起的舍入誤差，從而可大大降低算法的計算量，并且明顯改善了求解精度;② 可以通過人為設置控制求解精度，從而在效率和精度的矛盾中取得平衡。

本文還給出了使用此算法計算時的求解流程，并使用科學計算軟件FORTRAN對該算法進行編程實現(xiàn)，對其模塊化和軟件化提供了算法支持，為開發(fā)有自主知識產(chǎn)權的有限元求解器奠定基礎。實例證明了該算法對于大型非對稱矩陣特征值問題是有效的和精確的。

［1］淡丹輝，孫利民.結構動力有限元分析的阻尼建模與評價［J］.振動與沖擊，2007，26(2):121－125.

［2］榮見華，鄭健龍，徐飛鴻.結構動力修改及優(yōu)化設計［M］.北京:人民交通出版社，2002.

［3］宮玉才，周洪偉，陳 璞，等.快速子空間迭代法、迭代Ritz向量法與迭代 Lanczos法的比較［J］.振動工程學報.2005，18(2):228－232.

［4］張 琪.利用有限元和Lanczos法的細長彈體模態(tài)分析［J］.彈箭與制導學報，2007，27(4):61－63.

［5］王勖成.有限單元法［M］.清華大學出版社，2003.

［6］De Samblanx G，Bultheel A.Nested Lanczos:implicitly restarting an unsymmetric Lanczos algorithm［J］.Numerical Algorithms，1998，18:31 －50.

［7］ G.H.戈盧布，C.F.范洛恩.矩陣計算［M］.科學出版社，2001.

［8］A·R·高爾臘伊G·A·瓦特桑.矩陣特征值問題的計算方法［M］.上?？茖W技術出版社，1980，103－104.

［9］閆慶友，魏小鵬.求解大規(guī)模非對稱矩陣特征問題的精化雙正交Lanczos算法［J］.高等學校計算數(shù)學學報，2004，26(2):110－124.

［10］陳 前，朱德懋，周蒂蓮.粘彈性阻尼結構的動特性分析［J］.航空學報.1986，7(1):37－45.

［11］ Zhang J，Zheng G T.The Biot Model and Its Application in Viscoelastic Composite Structures［J］.Journal of Vibration and Acoustics，2007，129(5):533 －540.

Application of lanczos-QR algorithm in large non-classical damping structures'modal analysis

ZHANG Jing1，LIU Ming-hui1，ZHENG Gang-tie1，2

(1.Harbin Institute of Technology，Harbin 150001，China;2.Tsinghua University，Beijing 100084，China)

Here，Lanczos-QR algorithm was introduced into the computation of a large complex eigenvalue problem.Combining efficiency of Lanczos algorithm and stability of QR algorithm，this method firstly reduced matrix scale with Lanczos algorithm，and then solved complex eigenvalues and complex eigenvectors of the reduced matrix with QR algorithm.In order to improve the algorithm's computing stability and precision，the partitioned matrix triangle decomposition was introduced to get a standard eigenvalue problem.A step-by-step procedure of this method was summarized and programmed with FORTRAN language.Numerical examples showed that the improved Lanczos-QR algorithm is precise and stable.

complex eigenvalue problem;non-symmetrical damping;Lanczos algorithm;QR algorithm

TB123;V414

A

國防科技工業(yè)民用專項科研技術研究項目(C4220061310)

2009－12－14 修改稿收到日期:2010－03－30

張 靜 女，博士，1981年4月生