基于多重多級動力子結構的Lanczos算法

張 盛，方 杰， 張洪武， 陳飆松

(1.大連理工大學 工業裝備結構分析國家重點實驗室，大連 116024;2.大連理工大學 運載工程與力學學部工程力學系，大連 116024)

隨著工程技術的飛速發展，結構系統越來越龐大，越來越復雜，如飛機、大型輪船、高層建筑、大型機械和各種航天器。在分析計算大型復雜結構的動力特性和動力響應時，有限元離散化所得到的系統自由度是成千上萬階，有時甚至高達幾十萬階、幾百萬階。對于這種龐大的多自由度系統，用傳統的求解特征值方法求解是十分困難的，甚至是不可能的。子結構方法是計算大型復雜結構動態特性十分有效的方法。

對于廣義特征值問題［K］{φ}=λ［M］{φ}的求解，經典的子結構方法如模態綜合法［1］，計算精度的好壞直接由選取的位移表達式即假設模態決定。界面位移凝聚法［2］是通過對子結構剛度陣和質量陣凝聚處理，達到降階目的。胡海昌［3］利用解析的模態分析方法構造了約束界面模態綜合法。邱吉寶等［4］采用半解析法提出了三類精確動態子結構方法。而對于一般廣義特征值問題的求解，里茲向量法、子空間迭代法［5］、Lanczos算法［6］都是很實用的近似解法。喻永聲等［7］提出利用動態子結構周游技術實現子空間迭代，求解大型結構的廣義特征值方程，并取得較好的計算精度，但其需要迭代收斂判斷;文獻［8］提出利用子結構凝聚實現Lanczos算法的反迭代過程，實際上其求解仍然是在整體結構中進行。

Lanczos算法被認為是求解大型矩陣特征值問題的一種最有效的算法，由截斷 Lanczos過程產生的Lanczos矢量空間能有效地逼近結構離散化模型的低維狀態空間［9］，從而利用Lanczos矢量矩陣對數學模型進行降階，求解結構低階特征值。本文提出一種將多重多級子結構技術和Lanczos算法結合的方法。在子結構Lanczos迭代過程中，分別對每個子結構求解正交化系數和歸一化系數，然后累加形成總體正交化系數和歸一化系數，形成最終的三對角矩陣。數值算例結果表明該方法具有很高的計算效率和精度。

本文提出的基于多重多級子結構動力特征值分析的Lanczos算法已在具有自主版權的CAE軟件JIGFEX中實現。

1 多重多級子結構的Lanczos方法

對于廣義特征值方程:

其中［K］和［M］分別為結構的剛度矩陣和質量矩陣。文獻［7］提出了一種基于多重子結構的子空間迭代法，利用子結構周游樹技術實現子空間迭代。本文利用文獻中提到的子結構方法，實現基于多重子結構的Lanczos方法，計算效率較之于前者取得很大的進步。

Lanczos算法本質與子空間迭代法類似，都屬于向量反迭代法和Rayleigh-Ritz法相結合的一種方法，但它結合得巧妙，使計算過程大大簡化。其基本思想主要包括三步:反迭代、正交化和模歸一化處理。本文主要是利用多重多級子結構周游樹技術實現這三部分，并求解方程(1)。

考慮圖1所示的子結構模式劃分，整體結構為一平板，被離散成為一系列子結構模式，其中1、2、3、4、5為基本子結構模式，通過組裝調用形成6、7、8、9、10及11子結構模式，其中Sub11為頂層子結構模式。子結構周游樹總節點數為18，子結構周游樹中子結構調用關系見圖2所示。

方程(1)右端可看成廣義外力向量{F}，對于k≥1，Lanczos反迭代可依下式進行:

其中，

利用多重子結構靜力分析方法可以很方便求解式(2)中的位移向量}k+1。設子結構成員數為comp，s為子結構成員號。首先對子結構成員序列中的每一個子結構成員，利用式(3)形成廣義外力向量{F}k+1，并將方程(2)寫成如下形式:

其中，

求解頂層子結構靜力方程:

通過子結構前序周游回代求解下層每個子結構的內部位移向量:

圖2 中的子結構前序周游順序是 18，13，1，2，14，3，4，15，5，6，16，7，8，17，9，10，11，12。于是所有子結構位移自由度即可求出，完成式(2)中的反迭代。

Lancozs算法正交化由下式進行:

其中αk、βk是正交化系數，由下式確定:

式(10)中的位移向量都為整體結構下的位移向量，也是全體子結構位移向量的組合，顯然式(10)可以轉化為對每個子結構的正交化:

［Ms］為質量陣，且對于集中和協調質量陣同樣有效，于是可以得到:

計算中，為避免式(10)中的正交化不徹底，往往需要加入重正交化步驟，用上述的步驟同樣可以完成重正交化:

Lanczos向量的歸一化系數也可以利用類似的步驟完成。先考慮下式:

其中:

類似于式(14)、(15)，分別計算每個子結構位移向量的歸一化系數，并累加得:

分別對每個子結構歸一化:

對每個子結構進行模歸一化即可完成對整體結構的歸一化處理。事實上，上述的正交化和歸一化中，所有子結構位移自由度都參與了計算，這與整體結構參加計算完全相同，計算精度并不受子結構劃分的影響，計算量卻大大降低。

設求解特征值個數為n，特征子空間階數為m(一般取m=n×2)，k是迭代次數。程序流程如圖3所示，多重子結構Lanczos方法計算可按如下步驟進行:

(1)迭代次數k取0，對每個子結構s，隨機選取初始矢量{xs}0。

(2)迭代次數k進1。

(3)對每個子結構形成式(3)中的廣義外力向量{Fs}k。

(6)利用式(14)、(15)求解每個子結構的正交化系數并累加得到全局正交化系數。

(7)利用式(13)對每個子結構進行正交化。

(8)利用式(16)、(17)對每個子結構進行重正交化。

圖3 本文方法主要程序流程Fig.3 The main program flowchart of the proposed method

(9)利用式(20)求得每個子結構成員的歸一化系數，并累加得到全局歸一化系數。

(10)利用式(21)對每個子結構成員進行歸一化計算，并得到第k個Lanczos向量。

(11)返回步驟(2)，進行下一步迭代，直至k=m+1。

(12)形成由正交化系數組成的三對角陣。

對m×m階三對角陣的特征值方程:

求解得:

則方程(1)的解為:

選擇將隨機向量{xs}0進行一次反迭代和模歸一化后作為第一個Lanczos向量{xs}1。

本文算法與全結構算法計算效率的對比與靜力分析情況類似。算法的關鍵公式為平衡方程式(2)求解、式(10)正交化及式(18)歸一化，因此處未引入近似，故算法迭代次數以及式(10)、式(18)的計算復雜度與全結構Lanczos算法一致。在計算效率方面，主要考察式(2)，設全結構共有n個自由度，則全結構模型求解式(2)的計算復雜約為O(n3);若利用本文算法將結構等分為兩個子結構，子結構自由度為n/2，子結構的出口節點自由度為w，則按多重多級子結構求解式(2)的計算復雜度約為O((n/2)3+(n/2)3+w3)=O(n3/4+w3);一般w遠小于n，因此采用本文算法計算效率有優勢。對于復雜結構若子結構劃分得當，優勢將更顯著。

本文工作是在基于我國自主開發的CAE軟件系統JIGFEX 中實現的［10-11］。

2 計算實例

算例一 考慮圖1所示的矩形板，幾何尺寸10 m×20 m×0.05 m，邊界條件為:平面內位移固定，即dx=dy=θz=0;沿x=0 和x=10，dz=θy=0;沿x=0 和x=10，dz=θx=0;彈性模量 200 GPa，泊松比 0.3，密度8.0 ×103kg/m3。

采用圖1中的子結構模式劃分，利用本文方法求解前10階模態，同時對整體結構用MSC.Nastran求解，兩者都采用大小相同的三角形殼單元，總單元個數相等，前十階自振頻率計算結果如表1所示。計算結果表明本文子結構求解的計算精度與整體結構求解相當，顯示了本文方法具有很好的計算精度。圖4為平板結構的前4階振型。

表1 平板前10階頻率Tab.1 The free vibration frequencies of the flat

圖4 平板前四階振型Fig.4 The first four vibration model of plate

算例二 圖5所示為某型號運載火箭整體結構。箭體主體為圓筒加錐段組成，并附有夾筋梁，分為一子級、二子級、整流罩等幾個部分。模型的幾何尺寸為:一子級直徑為3.2 m，長度為16 m;二子級直徑3.2 m，長度8 m;三子級直徑為2.8 m，長度為6 m;整流罩直徑為4 m，長度為4 m。選用的材料為鋁合金:彈性模量 69 GPa，泊松比 0.3，密度 2.7 ×103kg/m3，梁截面尺寸為0.02 m×0.02 m。圖5(a)中整體結構有限元網格劃分共有節點7 037個，梁單元4 184個，三角形板單元13 840個。

圖5 三級火箭體子結構模式Fig.5 The substructure models of the three-leval rocket

按子結構方式劃分火箭，將整體火箭分為一子級、二子級、三子級、前錐段、后椎段，整流罩筒體段和整流罩錐段7段，根據結構圓周對稱的特點再將這7段沿圓周分為4部分，對這7個1/4火箭段進行建模，利用子結構的組裝調用最后形成圖5(b)中的整體結構。整個結構分為15個子結構模式，其中7個基本子結構模式，共有36個超級單元。圖5(b)中為火箭體的子結構模式序列。利用多重多級子結構Lanczos方法求解子結構模型前10階模態，并在MSC.Nastran中計算整體結構前10階自振頻率。

表2 火箭前10階自振頻率Tab.2 The free vibration frequencies of the rocket

表2為分別利用子結構模型與整體結構模型計算的前10階自振頻率，二者相對差值均不到0.5%，表現了多重多級子結構Lanczos方法良好的計算精度。

3 結論

本文方法在求解大型結構動力特性時具有較高的計算精度和效率。數值算例表明利用本文方法求解與整體結構計算結果相對差值皆在1%以下;較之多重子結構子空間迭代［7］，本文方法不需要迭代收斂判斷，每次迭代僅一個向量參與計算，效率更高。

［1］Hurty W C.Dynamic analysis of structural systems using component modes［J］.AIAA Journal，1965，3(4):678-685.

［2］ Guyen R J.Reduction of stiffness and mass matrice［J］.AIAA Journal，1965，3(2):380-380.

［3］胡海昌.多自由度結構固有振動理論［M］.北京:科學出版社，1987.

［4］邱吉寶，向樹紅，張正平.計算結構動力學［M］.合肥:中國科學技術大學出版社，2009.

［5］Bathe K J，Wilson E L.Solution methods for eigenva lue problems in structural mechanics［J］.International Journal for Numerical Methods in Engineering， 1973， 6(2):213-226.

［6］徐稼軒，鄭鐵生.結構動力分析的數值方法［M］.西安:西安交通大學出版社，1993.

［7］喻永聲，林家浩.超大型結構特征值問題求解的多重子結構子空間迭代［J］.工程力學，2003，20(6):149-156.

［8］張汝清，胡 寧.大型結構特征值問題的Lanczos子結構并行算法［J］.計算結構力學及其應用，1991，8(4):359-364.

［9］劉 豫，孫 秦.大型結構動力特性的Lanczos數值計算方法及程序設計［J］.科學技術與工程，2008，8(4):2010-1015.

［10］鐘萬勰，李錫夔.JIGFEX系統及其在大型結構分析中的應用［J］.土木工程學報.1982，15(3):20-28.

［11］ 張洪武，喻永聲，紀 崢，等.結構分析有限元程序JIGFEX和它的新發展［J］.計算結構力學及其應用，1995，12(3):292-297.