劉　丹　王　琥

湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國家重點(diǎn)實(shí)驗(yàn)室，長沙，410082

基于Lanczos法的模態(tài)重分析法在拓?fù)鋬?yōu)化中的應(yīng)用

劉丹王琥

湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國家重點(diǎn)實(shí)驗(yàn)室，長沙，410082

在拓?fù)鋬?yōu)化中,經(jīng)常要求對(duì)結(jié)構(gòu)進(jìn)行修改,快速準(zhǔn)確地計(jì)算修改后結(jié)構(gòu)的低階特征值對(duì)于提高整個(gè)結(jié)構(gòu)優(yōu)化的效率非常重要。將基于Lanczos算法的模態(tài)重分析法應(yīng)用于拓?fù)鋬?yōu)化過程中,利用初始結(jié)構(gòu)模態(tài)分析結(jié)果,結(jié)合Lanczos算法和投影技術(shù),采用縮減基方法求解修改結(jié)構(gòu)的固有頻率和振型, 則該方法同時(shí)具備了Lanczos向量快速收斂的優(yōu)點(diǎn)和基于全局近似的縮減基向量的高精度。剛架算例驗(yàn)證了該重分析法的高精度。固支方形板和車架結(jié)構(gòu)優(yōu)化結(jié)果表明,該方法在保證求解精度的同時(shí)能夠在一定程度上提高優(yōu)化迭代速度。

Lanczos算法;模態(tài)重分析;拓?fù)鋬?yōu)化;頻率優(yōu)化

0　引言

在拓?fù)鋬?yōu)化中,隨著結(jié)構(gòu)的不斷修改,需要求解一系列形如KiX=P的大規(guī)模線性方程組,其中Ki是在結(jié)構(gòu)拓?fù)鋬?yōu)化中第i步結(jié)構(gòu)剛度陣,X是待求結(jié)構(gòu)的響應(yīng),P是載荷向量。在結(jié)構(gòu)拓?fù)鋬?yōu)化中求解一系列大規(guī)模線性方程組已成為計(jì)算速度的瓶頸。結(jié)構(gòu)重分析技術(shù)就是利用初始結(jié)構(gòu)響應(yīng),設(shè)計(jì)高效重分析算法來求解修改后結(jié)構(gòu)的新響應(yīng)而不對(duì)新結(jié)構(gòu)進(jìn)行完整分析的技術(shù)，可以顯著降低計(jì)算成本,從而加速優(yōu)化過程。而修改結(jié)構(gòu)的模態(tài)分析比靜力分析更加耗時(shí),所以模態(tài)重分析就顯得更加重要。

近年來,很多學(xué)者研究了結(jié)合直接法和近似法形成改善的動(dòng)力學(xué)重分析方法的途徑,如黃海等[1]將攝動(dòng)和Padé逼近法、迭代和組合逼近相結(jié)合形成動(dòng)力學(xué)重分析方法。這些方法可用于結(jié)構(gòu)形狀修改和拓?fù)湫薷?包括自由度增加)的動(dòng)力學(xué)重分析。但到目前為止,還沒有真正形成一些十分簡便、快速、高效和通用的拓?fù)湫薷膭?dòng)力學(xué)重分析方法。Kirsch[2]將靜力學(xué)組合近似法引入結(jié)構(gòu)振動(dòng)分析過程,將模態(tài)方程和靜力平衡方程進(jìn)行了等效處理,在Krylov子空間內(nèi)構(gòu)造縮減基,求解降階的模態(tài)方程,并給出了誤差評(píng)估[3]。同樣,組合近似法對(duì)于結(jié)構(gòu)小修的計(jì)算效果較好，但是結(jié)構(gòu)修改量增大或需要計(jì)算高階模態(tài)時(shí)，常常會(huì)出現(xiàn)較大的誤差。為了提高組合近似法在結(jié)構(gòu)大修改時(shí)的求解精度,文獻(xiàn)[4-6]提出了擴(kuò)展的組合近似法,利用瑞利商提高了組合近似法的求解精度,但卻增加了計(jì)算量。Huang等[7]提出了基于擴(kuò)展基向量和瑞利-里茲分析的模態(tài)重分析方法。

結(jié)構(gòu)的頻率在許多工程領(lǐng)域都是關(guān)注的重點(diǎn)。工程結(jié)構(gòu)在動(dòng)載荷作用下的動(dòng)響應(yīng)在很大程度上依賴于結(jié)構(gòu)前幾階固有頻率。當(dāng)動(dòng)載荷頻率接近于被作用結(jié)構(gòu)的固有頻率之一時(shí)，該結(jié)構(gòu)會(huì)出現(xiàn)過度振動(dòng)。為了避免嚴(yán)重振動(dòng)，常常需要調(diào)整結(jié)構(gòu)基頻或前幾階頻率，使之遠(yuǎn)離動(dòng)載荷頻率范圍。因此，頻率優(yōu)化是結(jié)構(gòu)拓?fù)鋬?yōu)化中很重要的一部分。本文采用雙向漸進(jìn)優(yōu)化方法(evolutionary structural optimization,BESO)解決頻率優(yōu)化問題,該方法是在一定體積約束條件下，最大化第一階固有頻率。逐漸刪除材料的雙向漸進(jìn)結(jié)構(gòu)優(yōu)化概念在頻率優(yōu)化中已得到推廣和有效應(yīng)用。優(yōu)化過程中求解特征值問題時(shí)應(yīng)用基于Lanczos算法的模態(tài)重分析法可以在保證精度的同時(shí)，在一定程度上提高求解效率。

1　模態(tài)重分析

考慮初始無阻尼系統(tǒng)特征值問題：

K0Φ0=Λ0M0Φ0

(1)

其中,K0和M0分別為初始結(jié)構(gòu)剛度矩陣和質(zhì)量矩陣,Λ0和Φ0分別為初始結(jié)構(gòu)的特征值矩陣和相應(yīng)的模態(tài)向量組,結(jié)構(gòu)發(fā)生修改時(shí),修改后系統(tǒng)的特征值問題可以表示為

KmΦm=ΛmMmΦm

(2)

(3)

(4)

下標(biāo)n表示初始結(jié)構(gòu)的自由度,m表示修改結(jié)構(gòu)新增加的自由度。Λm和Φm為修改后結(jié)構(gòu)的特征值矩陣和相應(yīng)的模態(tài)向量組。直接求解式(2)可以得到精確解,但是為提高效率,我們可以利用初始分析的結(jié)果近似求解結(jié)構(gòu)修改之后的方程。文獻(xiàn)[7]提出基于擴(kuò)展基向量和瑞利-里茲分析的模態(tài)重分析方法，主要步驟如下。

(1)假設(shè)Φm和Λm的近似表達(dá)為

(5)

Λm≈Λ0

(6)

(2)把式(3)～式(6)代入式(2),展開可得

(ΔKmn-λ0iΔMmn)φ0i=-(ΔKmm-λ0iΔMmm)Δφi

(7)

所以有

Δφi=-(ΔKmm-λ0iΔMmm)-1·

(ΔKmn-λ0iΔMmn)φ0i

(8)

(3)三角分解剛度矩陣Km為

Km=LDLT

(9)

(4)求解靜態(tài)方程,獲得改善的里茲基向量:

(10)

(5)構(gòu)造縮減方程:

(11)

(12)

(6)最終的近似特征值和特征向量分別為

(13)

(14)

在此基礎(chǔ)上，改進(jìn)的單步攝動(dòng)瑞利商逆迭代法[8]將單步攝動(dòng)法和瑞利商逆迭代法結(jié)合起來,獲得了更好的效果。

2　基于Lanczos算法的模態(tài)重分析法

本文提供了一種將Lanczos算法、投影技術(shù)和縮減基法相結(jié)合的模態(tài)重分析方法。同主流重分析方法相比，本方法并沒有直接采用初始特征向量構(gòu)造縮減基向量組,而是基于Lanczos向量良好的正交性和快速收斂性[9],利用Lanczos法迭代出前幾階Lanczos向量做基向量。主要步驟如下。

(1)計(jì)算修改結(jié)構(gòu)的剛度矩陣K和質(zhì)量矩陣M。

(2)對(duì)剛度矩陣K進(jìn)行三角分解。

(3)選取初始分析所得特征向量Φ0為初始基向量構(gòu)造Lanczos向量：

(15)

q0是初始結(jié)構(gòu)的特征向量,正交化得

(16)

(17)

正則化得

(18)

qi即為迭代出的Lanczos向量。這就將廣義特征值問題KΦ=ΛMΦ轉(zhuǎn)換為Lanczos向量空間內(nèi)對(duì)角矩陣Tm的標(biāo)準(zhǔn)特征值問題，即求解

TmZ=Z λ

(19)

λ=diag(λi)Λ=λ-1

(4)構(gòu)造投影向量矩陣T。

(5)計(jì)算縮減矩陣KR和MR：

(20)

(6)求解特征值方程KRy=λMRy,得到前幾階固有頻率,再用線性組合法計(jì)算出相應(yīng)的主振型。

2.1截?cái)嗟腖anczos向量

選取已求得初始結(jié)構(gòu)的特征向量作為初始基向量，利用Lanczos算法迭代出前幾階向量。求解過程中，只需要采用PARDISO程序包一次三角分解稀疏矩陣K,之后就是前后回代計(jì)算。設(shè)v0是要求解的修改結(jié)構(gòu)的特征向量,r0是初始結(jié)構(gòu)求得的特征向量,令

(21)

則

(22)

迭代出基向量：

(23)

將式(22)代入式(23)得

(24)

正交化得

(25)

將式(20)和(23)代入式(24),正交化結(jié)果為

(26)

最后可簡化得

r1=H·Δv+ΔvTpv0

(27)

因此

|r1|≤|H||Δv|+|p||v0||Δv|=

(|H|+|p||v0|)|Δv|=c|Δv|

其中,c為一定常數(shù)。

2.2投影方法

(28)

(29)

設(shè)n是Lanczos向量截?cái)嚯A次,投影基矩陣的大小為nsdof×(n+1)j,如此縮減矩陣之后，需要求解的特征值矩陣規(guī)模就減小為(n+1)j×(n+1)j,大大縮短了CPU的運(yùn)算時(shí)間。

3　頻率優(yōu)化

動(dòng)力問題設(shè)計(jì)中,通常需要將結(jié)構(gòu)的基頻或低頻控制在一定范圍內(nèi)以避免共振,BESO在處理頻率優(yōu)化時(shí)簡單有效。與剛度優(yōu)化類似，頻率優(yōu)化也要根據(jù)靈敏度來確定需要?jiǎng)h除的單元。

3.1拓?fù)鋬?yōu)化問題

為了辨識(shí)結(jié)構(gòu)修改的最好位置，經(jīng)常需要進(jìn)行靈敏度分析。在有限元分析中，結(jié)構(gòu)動(dòng)特性一般特征值問題為

(30)

式中，ωi為第i階固有頻率；φi為對(duì)應(yīng)于ωi的特征矢量。

固有頻率及相應(yīng)的特征矢量通過瑞利商相互關(guān)聯(lián):

(31)

其中,模態(tài)剛度ki和模態(tài)質(zhì)量mi分別定義為

這里考慮拓?fù)鋬?yōu)化，最大化一階固有頻率，優(yōu)化問題表示如下：

(32)

式中,Vi為單個(gè)單元的體積;V*為預(yù)先設(shè)定的結(jié)構(gòu)體積;N為結(jié)構(gòu)中總的單元數(shù)目。

3.2材料插值方法

為了獲得設(shè)計(jì)變量的梯度信息,我們將材料在xmin和1之間進(jìn)行插值。假設(shè)原始的材料彈性模量和密度分別為E0和ρ0,則有：

(33)

式中,p為懲罰系數(shù)。

為了避免當(dāng)xi=xmin時(shí),由于對(duì)剛度和質(zhì)量的懲罰比值太大而引起局部奇異模態(tài)[10],因此,我們將材料插值表示為

(34)

根據(jù)式(31),可以求得目標(biāo)函數(shù)ω1的靈敏度為

(35)

根據(jù)式(31)以及特征向量關(guān)于質(zhì)量矩陣正交的性質(zhì),可以將式(35)簡化為

(36)

3.3拓?fù)鋬?yōu)化的進(jìn)化過程

采用雙向漸進(jìn)結(jié)構(gòu)優(yōu)化方法，逐漸從結(jié)構(gòu)中刪除材料,使結(jié)構(gòu)頻率能移向期望值。在優(yōu)化迭代過程的每一步，已經(jīng)算得當(dāng)前結(jié)構(gòu)的固有頻率與主振型，并確定出可刪除單元集合?？蓜h除單元是當(dāng)前結(jié)構(gòu)的一部分，利用已知的主振型即可對(duì)刪除某個(gè)單元所引起的特征值改變量進(jìn)行估計(jì)。本文的模態(tài)重分析方法可以加快優(yōu)化迭代速度,主要步驟如下：

(1)用有限元精細(xì)網(wǎng)格離散結(jié)構(gòu)；

(2)確定BESO方法的參數(shù)，如目標(biāo)體積V*、進(jìn)化率er和懲罰值p；

(3)采用有限元Lanczos方法求解特征值問題；

(4)用式(36)計(jì)算每個(gè)單元的靈敏度,通過平均單元過濾器內(nèi)的單元靈敏度和平均單元?dú)v史靈敏度來更新單元靈敏度，并刪除一定數(shù)量具有最小靈敏度的單元；

(5)采用本文基于Lanczos算法的模態(tài)重分析方法近似求解修改結(jié)構(gòu)的特征值；

(6)在結(jié)構(gòu)體積約束條件限制下,重復(fù)步驟(3)～步驟(5)直至達(dá)到最佳結(jié)構(gòu)需求。

4　數(shù)值算例

為了驗(yàn)證該方法的有效性，本文給出三個(gè)數(shù)值算例，第一個(gè)算例驗(yàn)證該重分析方法的精確性，后兩個(gè)算例驗(yàn)證該方法在結(jié)構(gòu)拓?fù)鋬?yōu)化過程中的應(yīng)用可行性。

4.1剛架結(jié)構(gòu)

如圖1a所示的平面剛架結(jié)構(gòu)，剛架的材料參數(shù)為：E=2.1×1011Pa,ρ=7800 kg/m3,梁單元的橫截面半徑r=0.02 m,節(jié)點(diǎn)1到2之間的長度是0.5 m,節(jié)點(diǎn)1到3之間的長度也是0.5 m,節(jié)點(diǎn)1和節(jié)點(diǎn)2固定。假設(shè)結(jié)構(gòu)拓?fù)湫薷闹笤黾恿?個(gè)節(jié)點(diǎn)和3個(gè)梁單元，如圖1b所示。同時(shí)剛架的材料參數(shù)變化為E=1.8×1011Pa,質(zhì)量密度變?yōu)棣?6800 kg/m3,節(jié)點(diǎn)1到2和節(jié)點(diǎn)1到3的長度都變?yōu)?.625 m。

(a)剛架初始結(jié)構(gòu)(b)剛架修改后結(jié)構(gòu)圖1　剛架結(jié)構(gòu)拓?fù)湫薷?/p>

以λi表示拓?fù)湫薷暮笥邢拊?jì)算參考解,λi0表示近似方法求得的解,下標(biāo)i表示第i階特征值。相對(duì)誤差定義為

從表1可以看出，本文方法求得前6階特征值近似值與有限元計(jì)算參考解是一致的，相對(duì)于其他方法，精度有一定提高。

表1　剛架結(jié)構(gòu)拓?fù)湫薷那昂笄?階特征值比較

4.2固支方形板

邊長為20 m的正方形板，厚度為1 m，四邊固支。在ANSYS中細(xì)分模型為40×40的四邊形殼單元網(wǎng)格。材料參數(shù)如下：彈性模量E=100 GPa，泊松比υ=0.3,密度ρ=7800 kg/m3,在板中心施加集中非結(jié)構(gòu)質(zhì)量m=3.12×105kg。BESO參數(shù)為：進(jìn)化率er=2%,最大增加率RA,max=2%，單元過濾半徑rmin=1.5 m,懲罰因子p=3.0,xmin=10-6,體積約束為初始體積的50%。

最終的拓?fù)鋬?yōu)化結(jié)果如圖2所示。

(a)固支方板初始設(shè)計(jì)域(b)固支方板拓?fù)鋬?yōu)化后模型圖2　固支方板

圖3顯示了基頻和體積分?jǐn)?shù)的迭代歷史進(jìn)程，最初基頻隨著體積的減小呈緩慢下降趨勢，在達(dá)到體積約束后，基頻逐漸收斂到一定值。其中ω1表示全分析求得的一階角頻率，ω2表示本文重分析方法求得的一階角頻率，可以看出兩者結(jié)果十分相近。

圖3　一階角頻率隨迭代次數(shù)變化曲線

優(yōu)化過程中，每近似分析十步進(jìn)行一次全分析，保證了解的正確性。

4.3車架結(jié)構(gòu)

輕型載貨汽車車架設(shè)計(jì)的初始模型如圖4所示，本算例是優(yōu)化車架橫梁的位置，車架的基本參數(shù)為：E=200 GPa,ρ=7800 kg/m3,μ=0.3。

圖4　車架設(shè)計(jì)初始模型

BESO參數(shù)為：er=2%,RA,max=2%，Rmin=75 mm,p=3.0,xmin=10-2，優(yōu)化區(qū)域?yàn)檐嚰軝M梁,設(shè)置縱梁為非優(yōu)化區(qū)域。優(yōu)化目標(biāo)為maxωi;約束條件為V*≥30%Vo。

優(yōu)化迭代歷史次數(shù)與一階角頻率和體積分?jǐn)?shù)的關(guān)系曲線如圖5所示。

圖5　一階角頻率參考解和近似解隨迭代次數(shù)變化曲線

最后拓?fù)鋬?yōu)化后的車架結(jié)構(gòu)如圖6所示，符合實(shí)體貨車車架的基本構(gòu)造，橫梁一般都布置在車架前部和后部以承收貨車大部分載荷。

圖6　車架拓?fù)鋬?yōu)化后模型

為保證特征解的收斂，優(yōu)化過程中，每近似分析五步就進(jìn)行一次全分析。由于每次最多只刪除總體積的2%，經(jīng)過60步迭代后，體積減小到原優(yōu)化區(qū)域體積的30%，在之后的幾步迭代中目標(biāo)值逐漸趨于穩(wěn)定。迭代停止的條件為

其中，k是當(dāng)前迭代步數(shù)；τ是允許的收斂誤差；N是一整數(shù)，本算例中取5，表示最后的10步迭代，特征值的改變量已經(jīng)足夠小。

從圖5中可以看出,采用本文方法優(yōu)化固有頻率所得近似解和參考解十分接近，一階固有角頻率從82.96 rad/s增加到110.36 rad/s,增加33%。同時(shí)，求解時(shí)間從1946 s縮短到1494 s,相較于全分析縮短了23%,保證求解精度的同時(shí)提高了求解效率。

5　結(jié)論

本文將Lanczos算法和縮減基法相結(jié)合,提出了一種高精度的模態(tài)重分析方法。該方法在結(jié)構(gòu)拓?fù)鋬?yōu)化中得到很好應(yīng)用，可以在保證求解精度的同時(shí)有效節(jié)省求解時(shí)間。

[1]黃海,陳塑寰,孟光.攝動(dòng)法結(jié)合Padé逼近在結(jié)構(gòu)拓?fù)渲胤治鲋械膽?yīng)用[J]. 應(yīng)用力學(xué)學(xué)報(bào),2005，22(2)： 155-158.

Huang Hai,Chen Suhuan,Meng Guang.The Application of Perturbation Method Combined with Padé Approximation in the Structural Topology Reanalysis[J].Chinese Journal of Applied Mechanics,2005,22(2):155-158.

[2]Kirsch U.Approximate Vibration Reanalysis of Structures[J].AIAA Journal,2003,41(3)：504-511.

[3]Kirsch U,Bogomolni M.Error Evaluation in Approximate Reanalysis of Structures[J].Structural and Multidisciplinary Optimization,2004,28(2/3):77-86.

[4]呂振華.結(jié)構(gòu)動(dòng)力學(xué)修改重分析方法的發(fā)展[J].計(jì)算結(jié)構(gòu)力學(xué)及其應(yīng)用,1994,11(1):85-91.

Lü Zhenhua.Development of Reanalysis Method for Structural Dynamics Modifications[J].Computional Structural Mechanics and Application,1994,11(1):85-91. [5]Chen Suhuan, Yang Xiaowei. Extended Kirsch Combined Method for Eigenvalue Reanalysis[J].AIAA Journal, 2000, 38(5): 927-930.[6]Feng Rong, Chen Suhuan, Chen Yudong. Structural Modal Reanalysis for Topological Modifications with Extended Kirsch Method[J]. Computer Methods in Applied Mechanics and Engineering，2003,192(5):697-707.

[7]Huang Cheng，Chen Suhuan，Liu Zhongsheng. Structural Modal Reanalysis for Topological Modifications of Finite Element Systems[J]. Engineering Structures，2000，22：304-310.

[8]何建軍. 結(jié)構(gòu)拓?fù)湫薷膭?dòng)力學(xué)重分析方法的研究[D]. 西安：西北工業(yè)大學(xué)，2006.

[9]王欣欣. 求解對(duì)稱矩陣特征值問題的Lanczos算法的改進(jìn)及分析[D]. 哈爾濱：哈爾濱工業(yè)大學(xué)，2008.

[10]Huang Xiaodong,Zuo Zhihao,Xie Yimin.Evolutionary Topological Optimization of Vibrating Continuum Structures for Natural Frequencies[J]. Computers and Structures, 2010，88(5):357-364.

(編輯袁興玲)

Application of Lanczos-based Modal Reanalysis Algorithm in Topological Optimization

Liu DanWang Hu

State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body,Hunan University,Changsha,410082

In topological optimization,the structures were often required to make some modifications,and it was very important for improving the effciency of optimization that the low-order eigenvalues could be calculated accurately and quickly.The Lanczos-based modal reanalysis method was applied to topological modification.The presented method which combined Lanczos algorithm and projection techniques,solved the natural frequency and mode shapes for the modified structures utilizing the already obtained eigenvectors with reduced-basis algorithm.The major advantages laid that the fast convergence of Lanczos vector and high-accuracy of reduced-basis method based on global approximations.A stiff frame example was demonstrated the high accuracy of the reanalysis method.The optimization results of the clamped square plate and vehicle frame indicate that the method can also ensure the calculation accuracy and optimization iteration speed to a certain extent.

Lanczos algorithm;modal reanalysis;topological optimization;frequency optimization

2014-01-21

國家自然科學(xué)基金資助項(xiàng)目(11172097)；教育部新世紀(jì)優(yōu)秀人才支持計(jì)劃資助項(xiàng)目(NCET-11-0131);湖南省自然科學(xué)基金資助項(xiàng)目(11JJA001)

TH11;O32DOI：10.3969/j.issn.1004-132X.2015.11.015

劉丹,女,1987年生。湖南大學(xué)機(jī)械與運(yùn)載工程學(xué)院碩士研究生。研究方向?yàn)槟B(tài)重分析及結(jié)構(gòu)優(yōu)化。王琥，男，1975年生。湖南大學(xué)機(jī)械與運(yùn)載工程學(xué)院副教授。