劉　丹　王　琥

湖南大學汽車車身先進設計制造國家重點實驗室，長沙，410082

基于Lanczos法的模態(tài)重分析法在拓撲優(yōu)化中的應用

劉丹王琥

湖南大學汽車車身先進設計制造國家重點實驗室，長沙，410082

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

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

0　引言

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

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

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

1　模態(tài)重分析

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

K0Φ0=Λ0M0Φ0

(1)

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

KmΦm=ΛmMmΦm

(2)

(3)

(4)

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

(1)假設Φm和Λm的近似表達為

(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)構造縮減方程:

(11)

(12)

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

(13)

(14)

在此基礎上，改進的單步攝動瑞利商逆迭代法[8]將單步攝動法和瑞利商逆迭代法結合起來,獲得了更好的效果。

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

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

(1)計算修改結構的剛度矩陣K和質(zhì)量矩陣M。

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

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

(15)

q0是初始結構的特征向量,正交化得

(16)

(17)

正則化得

(18)

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

TmZ=Z λ

(19)

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

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

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

(20)

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

2.1截斷的Lanczos向量

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

(21)

則

(22)

迭代出基向量：

(23)

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

(24)

正交化得

(25)

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

(26)

最后可簡化得

r1=H·Δv+ΔvTpv0

(27)

因此

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

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

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

2.2投影方法

(28)

(29)

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

3　頻率優(yōu)化

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

3.1拓撲優(yōu)化問題

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

(30)

式中，ωi為第i階固有頻率；φi為對應于ωi的特征矢量。

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

(31)

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

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

(32)

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

3.2材料插值方法

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

(33)

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

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

(34)

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

(35)

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

(36)

3.3拓撲優(yōu)化的進化過程

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

(1)用有限元精細網(wǎng)格離散結構；

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

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

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

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

(6)在結構體積約束條件限制下,重復步驟(3)～步驟(5)直至達到最佳結構需求。

4　數(shù)值算例

為了驗證該方法的有效性，本文給出三個數(shù)值算例，第一個算例驗證該重分析方法的精確性，后兩個算例驗證該方法在結構拓撲優(yōu)化過程中的應用可行性。

4.1剛架結構

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

(a)剛架初始結構(b)剛架修改后結構圖1　剛架結構拓撲修改

以λi表示拓撲修改后有限元計算參考解,λi0表示近似方法求得的解,下標i表示第i階特征值。相對誤差定義為

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

表1　剛架結構拓撲修改前后前6階特征值比較

4.2固支方形板

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

最終的拓撲優(yōu)化結果如圖2所示。

(a)固支方板初始設計域(b)固支方板拓撲優(yōu)化后模型圖2　固支方板

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

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

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

4.3車架結構

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

圖4　車架設計初始模型

BESO參數(shù)為：er=2%,RA,max=2%，Rmin=75 mm,p=3.0,xmin=10-2，優(yōu)化區(qū)域為車架橫梁,設置縱梁為非優(yōu)化區(qū)域。優(yōu)化目標為maxωi;約束條件為V*≥30%Vo。

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

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

最后拓撲優(yōu)化后的車架結構如圖6所示，符合實體貨車車架的基本構造，橫梁一般都布置在車架前部和后部以承收貨車大部分載荷。

圖6　車架拓撲優(yōu)化后模型

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

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

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

5　結論

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

[1]黃海,陳塑寰,孟光.攝動法結合Padé逼近在結構拓撲重分析中的應用[J]. 應用力學學報,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]呂振華.結構動力學修改重分析方法的發(fā)展[J].計算結構力學及其應用,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]何建軍. 結構拓撲修改動力學重分析方法的研究[D]. 西安：西北工業(yè)大學，2006.

[9]王欣欣. 求解對稱矩陣特征值問題的Lanczos算法的改進及分析[D]. 哈爾濱：哈爾濱工業(yè)大學，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

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

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

劉丹,女,1987年生。湖南大學機械與運載工程學院碩士研究生。研究方向為模態(tài)重分析及結構優(yōu)化。王琥，男，1975年生。湖南大學機械與運載工程學院副教授。