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

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

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

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

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

本文提出的基于多重多級子結(jié)構(gòu)動力特征值分析的Lanczos算法已在具有自主版權(quán)的CAE軟件JIGFEX中實現(xiàn)。

1 多重多級子結(jié)構(gòu)的Lanczos方法

對于廣義特征值方程:

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

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

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

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

其中，

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

其中，

求解頂層子結(jié)構(gòu)靜力方程:

通過子結(jié)構(gòu)前序周游回代求解下層每個子結(jié)構(gòu)的內(nèi)部位移向量:

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

Lancozs算法正交化由下式進(jìn)行:

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

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

［Ms］為質(zhì)量陣，且對于集中和協(xié)調(diào)質(zhì)量陣同樣有效，于是可以得到:

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

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

其中:

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

分別對每個子結(jié)構(gòu)歸一化:

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

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

(1)迭代次數(shù)k取0，對每個子結(jié)構(gòu)s，隨機(jī)選取初始矢量{xs}0。

(2)迭代次數(shù)k進(jìn)1。

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

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

(7)利用式(13)對每個子結(jié)構(gòu)進(jìn)行正交化。

(8)利用式(16)、(17)對每個子結(jié)構(gòu)進(jìn)行重正交化。

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

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

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

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

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

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

求解得:

則方程(1)的解為:

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

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

本文工作是在基于我國自主開發(fā)的CAE軟件系統(tǒng)JIGFEX 中實現(xiàn)的［10-11］。

2 計算實例

算例一 考慮圖1所示的矩形板，幾何尺寸10 m×20 m×0.05 m，邊界條件為:平面內(nèi)位移固定，即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中的子結(jié)構(gòu)模式劃分，利用本文方法求解前10階模態(tài)，同時對整體結(jié)構(gòu)用MSC.Nastran求解，兩者都采用大小相同的三角形殼單元，總單元個數(shù)相等，前十階自振頻率計算結(jié)果如表1所示。計算結(jié)果表明本文子結(jié)構(gòu)求解的計算精度與整體結(jié)構(gòu)求解相當(dāng)，顯示了本文方法具有很好的計算精度。圖4為平板結(jié)構(gòu)的前4階振型。

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

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

算例二 圖5所示為某型號運載火箭整體結(jié)構(gòu)。箭體主體為圓筒加錐段組成，并附有夾筋梁，分為一子級、二子級、整流罩等幾個部分。模型的幾何尺寸為:一子級直徑為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)中整體結(jié)構(gòu)有限元網(wǎng)格劃分共有節(jié)點7 037個，梁單元4 184個，三角形板單元13 840個。

圖5 三級火箭體子結(jié)構(gòu)模式Fig.5 The substructure models of the three-leval rocket

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

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

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

3 結(jié)論

本文方法在求解大型結(jié)構(gòu)動力特性時具有較高的計算精度和效率。數(shù)值算例表明利用本文方法求解與整體結(jié)構(gòu)計算結(jié)果相對差值皆在1%以下;較之多重子結(jié)構(gòu)子空間迭代［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］胡海昌.多自由度結(jié)構(gòu)固有振動理論［M］.北京:科學(xué)出版社，1987.

［4］邱吉寶，向樹紅，張正平.計算結(jié)構(gòu)動力學(xué)［M］.合肥:中國科學(xué)技術(shù)大學(xué)出版社，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］徐稼軒，鄭鐵生.結(jié)構(gòu)動力分析的數(shù)值方法［M］.西安:西安交通大學(xué)出版社，1993.

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

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

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

［10］鐘萬勰，李錫夔.JIGFEX系統(tǒng)及其在大型結(jié)構(gòu)分析中的應(yīng)用［J］.土木工程學(xué)報.1982，15(3):20-28.

［11］ 張洪武，喻永聲，紀(jì) 崢，等.結(jié)構(gòu)分析有限元程序JIGFEX和它的新發(fā)展［J］.計算結(jié)構(gòu)力學(xué)及其應(yīng)用，1995，12(3):292-297.