Lanczos-QR方法在大型非比例阻尼結(jié)構(gòu)復(fù)模態(tài)計(jì)算中的應(yīng)用

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

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

Lanczos-QR方法在大型非比例阻尼結(jié)構(gòu)復(fù)模態(tài)計(jì)算中的應(yīng)用

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

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

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

特征值問(wèn)題;非比例阻尼;Lanczos算法;QR算法

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

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

1 非比例阻尼結(jié)構(gòu)的特征值問(wèn)題

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

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

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

對(duì)于式(3)所示的非對(duì)稱廣義特征值問(wèn)題，若矩陣B非奇異，則可化為如下的標(biāo)準(zhǔn)特征值問(wèn)題的形式:

2 非對(duì)稱矩陣特征問(wèn)題的Lanczos算法

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

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

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

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

3 具有原點(diǎn)位移的雙重步QR算法

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

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

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

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

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

4 復(fù)特征向量的變換

5 標(biāo)準(zhǔn)特征值轉(zhuǎn)化的分塊Cholesky分解方法

對(duì)于標(biāo)準(zhǔn)特征值問(wèn)題:

其一，若剛度矩陣K奇異，則無(wú)法進(jìn)行求逆;

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

其三，對(duì)于大型矩陣，直接求逆運(yùn)算的計(jì)算效率較低。

為了改善求解過(guò)程中的數(shù)值計(jì)算精度和效率，可以采用分塊Cholesky分解來(lái)求出標(biāo)準(zhǔn)特征值問(wèn)題的系數(shù)矩陣。

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

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

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

Lanczos-QR算法的主要思想是利用求解絕對(duì)值最小的特征值的逆冪法思想來(lái)改進(jìn)截?cái)喔袷降腖anczos方法降低矩陣階數(shù)，得到一個(gè)包含原矩陣的低階特征值的三對(duì)角矩陣，再使用帶移位的雙重步QR方法對(duì)降階后的三對(duì)角矩陣進(jìn)行全部特征值求解，以避免計(jì)算過(guò)程中的復(fù)數(shù)運(yùn)算。

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

使用改進(jìn)截?cái)喔袷降腖anczos算法求解大型復(fù)特征值問(wèn)題軟件實(shí)現(xiàn)流程如下。

圖1 程序設(shè)計(jì)流程圖Fig.1 Programming flowchart

7 數(shù)值算例

7.1 約束阻尼梁

為了驗(yàn)證本文提出的算法在復(fù)特征值求解時(shí)的有效性和準(zhǔn)確性，采用含粘彈性阻尼的約束阻尼結(jié)構(gòu)來(lái)進(jìn)行計(jì)算。如圖2所示的約束阻尼懸臂梁，阻尼材料選用ZN－1型粘彈性材料。具體參數(shù)如下:

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

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

表1 頻率和阻尼比計(jì)算結(jié)果Tab.1 Frequency and damping ratio results

7.2 大型衛(wèi)星結(jié)構(gòu)

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

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

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

表2 頻率計(jì)算結(jié)果Tab.2 Frequency results

8 結(jié)論

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

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

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

［2］榮見(jiàn)華，鄭健龍，徐飛鴻.結(jié)構(gòu)動(dòng)力修改及優(yōu)化設(shè)計(jì)［M］.北京:人民交通出版社，2002.

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

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

［5］王勖成.有限單元法［M］.清華大學(xué)出版社，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.范洛恩.矩陣計(jì)算［M］.科學(xué)出版社，2001.

［8］A·R·高爾臘伊G·A·瓦特桑.矩陣特征值問(wèn)題的計(jì)算方法［M］.上海科學(xué)技術(shù)出版社，1980，103－104.

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

［10］陳 前，朱德懋，周蒂蓮.粘彈性阻尼結(jié)構(gòu)的動(dòng)特性分析［J］.航空學(xué)報(bào).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

國(guó)防科技工業(yè)民用專項(xiàng)科研技術(shù)研究項(xiàng)目(C4220061310)

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

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