非比例阻尼結(jié)構(gòu)諧響應(yīng)分析的自適應(yīng)模型降階方法①

李玉韋,魏 曉,曹 航,柏漢松,王 博,郝 鵬

(1.中國(guó)航發(fā)沈陽(yáng)發(fā)動(dòng)機(jī)研究所,沈陽(yáng) 110015;2. 大連理工大學(xué) 工程力學(xué)系,大連 116024)

0 引言

為了在服役過(guò)程中控制結(jié)構(gòu)的振動(dòng)水平,通常需要獲得結(jié)構(gòu)在簡(jiǎn)諧激勵(lì)下的穩(wěn)態(tài)響應(yīng)。工程中常用的諧響應(yīng)分析方法包括直接頻響法、模態(tài)疊加法和模態(tài)加速度法等。直接頻響法因其計(jì)算精度高,常被應(yīng)用于計(jì)算結(jié)構(gòu)的穩(wěn)態(tài)響應(yīng)[1],但該方法計(jì)算效率較低,不適合處理激勵(lì)頻率較為密集或自由度數(shù)目較大的工程問(wèn)題。模態(tài)疊加法通過(guò)引入模態(tài)正交化條件解耦動(dòng)力學(xué)方程,提高了復(fù)雜結(jié)構(gòu)的諧響應(yīng)分析效率,被廣泛應(yīng)用于復(fù)雜工程結(jié)構(gòu)[2-3]。由于模態(tài)疊加法忽略高階模態(tài)的影響,導(dǎo)致該方法計(jì)算精度降低。為提高模態(tài)疊加法的計(jì)算精度,需要選擇更多模態(tài)參與諧響應(yīng)分析,但對(duì)于如何選擇模態(tài)以及選擇多少模態(tài)參與計(jì)算還缺乏科學(xué)的準(zhǔn)則。為了彌補(bǔ)模態(tài)疊加法忽略高階模態(tài)導(dǎo)致計(jì)算精度降低的問(wèn)題,模態(tài)加速度法引入擬靜力響應(yīng)來(lái)補(bǔ)償截?cái)嗄B(tài)的貢獻(xiàn),從而達(dá)到用較少模態(tài)獲得高精度響應(yīng)的效果[4-5]。模態(tài)疊加法或模態(tài)加速度法計(jì)算結(jié)構(gòu)穩(wěn)態(tài)響應(yīng)的高效性體現(xiàn)在結(jié)構(gòu)阻尼矩陣滿足模態(tài)正交化條件,從而可以對(duì)動(dòng)力學(xué)方程進(jìn)行解耦,提高響應(yīng)分析的效率。但是,對(duì)于由不同材料和構(gòu)件組成的非比例阻尼結(jié)構(gòu),阻尼矩陣不再滿足模態(tài)正交化條件,此時(shí),模態(tài)疊加法或模態(tài)加速度法不再適用,否則會(huì)導(dǎo)致錯(cuò)誤的計(jì)算結(jié)果[6-7]。

Krylov子空間法[8]通過(guò)構(gòu)建一組標(biāo)準(zhǔn)正交基來(lái)實(shí)現(xiàn)模型降階,提高結(jié)構(gòu)響應(yīng)計(jì)算效率。這種模型降階方法數(shù)學(xué)理論完善、算法穩(wěn)定,亦不要求結(jié)構(gòu)阻尼矩陣滿足模態(tài)正交化條件,因此該方法被廣泛應(yīng)用于工程結(jié)構(gòu)的模型降階[9-10]。一階Krylov子空間法應(yīng)用于結(jié)構(gòu)模型降階時(shí),需要將二階動(dòng)力學(xué)方程轉(zhuǎn)換為一階狀態(tài)方程,再對(duì)狀態(tài)方程降階[11]。這種線性化處理方式的優(yōu)點(diǎn)是收斂速度快,但線性化過(guò)程破壞了原結(jié)構(gòu)的矩陣特征,計(jì)算規(guī)模是原動(dòng)力學(xué)方程的兩倍,計(jì)算量顯著增加。為此,基于二階Krylov子空間作為前射子空間,柏兆俊和蘇仰峰提出了二階Arnoldi(Second-order Arnoldi Reduction,SOAR)算法[12-14],該算法不會(huì)破壞結(jié)構(gòu)動(dòng)力學(xué)方程的二階特性和結(jié)構(gòu)矩陣特點(diǎn),在保證分析精度的前提下,可以有效減少數(shù)值分析的復(fù)雜度,提高計(jì)算效率。

針對(duì)SOAR方法需要使用昂貴的矩陣反演生成正交投影矩陣且無(wú)法保持原系統(tǒng)穩(wěn)定性的問(wèn)題,MALIK 等[15]用高斯核加權(quán)的插值點(diǎn)在期望的頻率范圍內(nèi)生成減縮基矩陣。TAMRI等[16]提出基于數(shù)值秩性能系數(shù)的概念來(lái)自動(dòng)選擇最優(yōu)的降階模型。本文提出一種簡(jiǎn)單的自適應(yīng)模型降階方法,利用交叉驗(yàn)證和二分法策略確定展開(kāi)頻點(diǎn)數(shù)和正交基階數(shù),自適應(yīng)地建立在目標(biāo)頻段內(nèi)具有更高計(jì)算精度的降階模型,提高非比例阻尼結(jié)構(gòu)諧響應(yīng)分析效率。

1 諧響應(yīng)分析方法

簡(jiǎn)諧激勵(lì)載荷下結(jié)構(gòu)的頻響方程為

(K+jωC-ω2M)U(ω)=F

(1)

直接頻響法將一系列離散頻點(diǎn)直接代入式(1)得到結(jié)構(gòu)在頻點(diǎn)ω處的響應(yīng):

U(ω)=(K+jωC-ω2M)-1F

(2)

可以看出,直接頻響法需要得到每個(gè)頻點(diǎn)ω處動(dòng)剛度陣的逆矩陣,計(jì)算規(guī)模較大,不適用于激勵(lì)頻率密集或自由度數(shù)目較多的問(wèn)題。

1.1 模態(tài)疊加法

模態(tài)疊加法是把對(duì)應(yīng)無(wú)阻尼結(jié)構(gòu)的模態(tài)為空間基底,經(jīng)過(guò)坐標(biāo)變換使得原動(dòng)力學(xué)方程解耦,通過(guò)求解n個(gè)獨(dú)立的方程獲得模態(tài)位移,進(jìn)而通過(guò)疊加各階模態(tài)的貢獻(xiàn)求得結(jié)構(gòu)的響應(yīng)。

假定結(jié)構(gòu)特征值和模態(tài)分別是Λ和Φ:

Φ=[φ1,φ2,…,φn]

(3)

其中,Λ和Φ滿足特征值方程和正交性條件:

KΦ=ΛMΦ,ΦTKΦ=Λ,ΦTMΦ=I

(4)

若阻尼為瑞利阻尼,則阻尼陣C滿足模態(tài)正交化條件:

ΦTCΦ=2ξdiag(ω1,ω2,…,ωn)

(5)

式中ξ為瑞利阻尼比。

利用式(4)和式(5)將式(1)解耦為

(6)

可以看出,模態(tài)疊加法通過(guò)模態(tài)的線性組合來(lái)近似結(jié)構(gòu)的位移,避免了動(dòng)剛度陣的分解,提高了計(jì)算效率。實(shí)用中,模態(tài)疊加法一般是通過(guò)模態(tài)截?cái)噙x取有限數(shù)目的模態(tài)近似表示結(jié)構(gòu)的位移:

(7)

式中l(wèi)為選取的模態(tài)數(shù),l<

模態(tài)疊加法忽略了高階模態(tài),使得該方法在高頻范圍內(nèi)的求解精度降低。

1.2 模態(tài)加速度法

對(duì)于無(wú)阻尼結(jié)構(gòu),式(1)可以改寫為

U(ω)=K-1F+ω2K-1MU(ω)

=K-1F+ω2K-1MΦ(Λ-ω2I)-1ΦTF

(8)

將式(4)代入式(8)可得

U(ω)=K-1F+ω2K-1MU(ω)

=K-1F+ω2K-1MΦ(Λ-ω2I)-1ΦTF

(9)

忽略高階模態(tài),式(9)可變換為如下形式:

(10)

可以看出,模態(tài)加速度法通過(guò)引入擬靜力響應(yīng)來(lái)補(bǔ)償截?cái)嗄B(tài)的貢獻(xiàn),從而達(dá)到用較少模態(tài)獲得高精度響應(yīng)的效果。

2 基于SOAR的自適應(yīng)模型降階方法

2.1 SOAR算法

在任一頻點(diǎn)ωa處,頻響方程(1)的等價(jià)形式如下:

(11)

其中,

(12)

通過(guò)二階Krylov子空間對(duì)式(12)進(jìn)行降階,對(duì)應(yīng)減縮基矩陣如下:

Tk(ωa)=span{r0(ωa),r1(ωa),…,rk-1(ωa)}

(13)

其中,r0(ωa),r1(ωa),…,rk-1(ωa)代表在頻點(diǎn)ωa處展開(kāi)的k個(gè)標(biāo)準(zhǔn)正交基向量:

(14)

建立減縮基矩陣Tk(ωa)后,式(11)可寫為

(15)

其中,

(16)

可以看出,基于二階Krylov子空間建立降階模型與展開(kāi)頻點(diǎn)數(shù)及標(biāo)準(zhǔn)正交基階數(shù)有關(guān),展開(kāi)頻點(diǎn)數(shù)和標(biāo)準(zhǔn)正交基階數(shù)越多,降階模型的計(jì)算精度越高。為平衡降階模型的計(jì)算精度和效率,本文提出了自適應(yīng)確定展開(kāi)頻點(diǎn)數(shù)和標(biāo)準(zhǔn)正交基階數(shù)的降階策略。

2.2 自適應(yīng)降階策略

結(jié)構(gòu)動(dòng)柔順度[17]誤差常被用于衡量降階模型的計(jì)算精度:

(17)

式(17)需要獲得所有頻點(diǎn)的位移來(lái)計(jì)算結(jié)構(gòu)的動(dòng)柔順度誤差,這使得建立降階模型過(guò)程中需要浪費(fèi)大量的時(shí)間去評(píng)估降階模型的計(jì)算精度,失去了利用降階模型提高計(jì)算效率的意義。本文提出利用交叉驗(yàn)證的策略來(lái)評(píng)估降階模型的預(yù)測(cè)精度:

(18)

交叉驗(yàn)證策略僅需計(jì)算展開(kāi)頻點(diǎn)的位移精確解,大大提高了建立降階模型的效率。為了在保證降階模型計(jì)算精度的前提下,盡可能減少展開(kāi)頻點(diǎn)的數(shù)目,本文建立了基于二分法的自適應(yīng)降階策略,具體流程如圖1所示。

圖1 基于SOAR算法的自適應(yīng)模型降階流程Fig.1 Flowchart of the SOAR-based adaptive model order reduction method

第一步:設(shè)置每個(gè)頻點(diǎn)的初始標(biāo)準(zhǔn)正交基階數(shù)r和最大正交基數(shù)目Maxorder,正交基階數(shù)增加的步長(zhǎng)p;設(shè)置結(jié)構(gòu)動(dòng)柔順度誤差閾值Tol;初始化集合S1={ω1}和集合S2={ω2},其中ω1和ω2分別為目標(biāo)頻段的上下限。

第二步:利用SOAR算法得到展開(kāi)頻點(diǎn)ωa處的r個(gè)標(biāo)準(zhǔn)正交基,并組裝成減縮基矩陣Tr(ωa):

ωa=max{S1}

Tr(ωa)=span{r0(ωa),r1(ωa),…,rr-1(ωa)}

(19)

ωb=min{S2}

(20)

第四步:通過(guò)SOAR算法擴(kuò)充展開(kāi)頻點(diǎn)ωa處的標(biāo)準(zhǔn)正交基數(shù)目,并組裝成減縮基矩陣Tr+p(ωa):

Tr+p(ωa)=span{r0(ωa),r1(ωa),…,rr+p-1(ωa)}

(21)

第五步:利用減縮基矩陣Tr+p(ωa)計(jì)算展開(kāi)頻點(diǎn) 處的動(dòng)柔順度誤差:

(22)

S1=S1∪{ωc}

(23)

S2=S2∪{ωa}

S1=S1{ωa}

(24)

若移除頻點(diǎn)ωa后,集合S1仍為非空集合,則返回至第二步,否則迭代結(jié)束。

迭代結(jié)束后,將集合S2中所有頻點(diǎn)的標(biāo)準(zhǔn)正交基進(jìn)行正交化得到減縮基矩陣T:

T=orth[Tr+p(ω1),Tr+p(ω2),…,Tr+p(ωi)]

(25)

其中,orth表示對(duì)矩陣進(jìn)行正交化。本文采用奇異值分解技術(shù)(Singular Value Decomposition,SVD)進(jìn)行正交化。

3 自適應(yīng)模型降階方法驗(yàn)證

3.1 阻尼涂層板的諧響應(yīng)分析

本節(jié)采用阻尼涂層板驗(yàn)證自適應(yīng)降階策略有效性,阻尼涂層板結(jié)構(gòu)長(zhǎng)300 mm,寬25 mm ,厚3.5 mm,其中阻尼涂層厚2 mm,如圖2所示。基體材料的彈性模量E=71 GPa,密度ρ=2700 kg/m3,泊松比μ=0.3,阻尼比ξ=0.02。阻尼涂層材料的彈性模量E=50 GPa,密度ρ=1400 kg/m3,泊松比μ=0.3,阻尼比ξ=0.05。基體結(jié)構(gòu)右端固支,左端施加簡(jiǎn)諧激勵(lì),激勵(lì)頻率的范圍為[0, 800 Hz], 分析步長(zhǎng)為2 Hz。

圖2 阻尼涂層平板Fig.2 A plate with free-layer damping

選取初始標(biāo)準(zhǔn)正交基階數(shù)r=3、步長(zhǎng)p=1,最大正交基數(shù)目Maxorder=20,根據(jù)圖1給出的自適應(yīng)模型降階流程,在目標(biāo)頻段內(nèi)確定了5個(gè)展開(kāi)頻點(diǎn)以及每個(gè)頻點(diǎn)處所需標(biāo)準(zhǔn)正交基的數(shù)目,如表 1所示。將5個(gè)展開(kāi)頻點(diǎn)生成的標(biāo)準(zhǔn)正交基進(jìn)行SVD正交化得到降階模型的減縮基矩陣T。

表1 自適應(yīng)確定的展開(kāi)頻點(diǎn)數(shù)和標(biāo)準(zhǔn)正交基階數(shù)Table 1 Expansion frequency points and its related orders of orthonormal basis frequency points

圖3 (a)分別給出了基于全階模型(FOM)和降階模型(ROM_T)得到的動(dòng)柔順度曲線。可以看出,降階模型與全階模型的計(jì)算結(jié)果十分吻合,相對(duì)誤差最大值不超過(guò)1×10-7(見(jiàn)圖3(b))中黑色實(shí)線)。

(a) Dynamic compliance

(b) Relative error圖3 阻尼涂層平板的動(dòng)柔順度及相對(duì)誤差Fig.3 Dynamic compliance and relative error of the free-layer damping plate

為了驗(yàn)證本文提出的交叉驗(yàn)證策略有效性,圖3 (b)分別給出了由單個(gè)展開(kāi)頻點(diǎn)的正交基建立的降階模型預(yù)測(cè)誤差。圖中括號(hào)內(nèi)的數(shù)字表示展開(kāi)頻點(diǎn),下標(biāo)表示標(biāo)準(zhǔn)正交基的階數(shù),例如,T4(200)表示展開(kāi)頻點(diǎn)為200 Hz,該展開(kāi)頻點(diǎn)處的標(biāo)準(zhǔn)正交基階數(shù)為4。圖中豎直灰色實(shí)線代表展開(kāi)頻點(diǎn),水平灰色虛線代表閾值Tol。可以看出,降階模型(ROM_T4(200))在頻點(diǎn)200 Hz附近頻段內(nèi)的預(yù)測(cè)精度較高,在遠(yuǎn)離200 Hz的頻段內(nèi)的預(yù)測(cè)精度較差。同樣,降階模型(ROM_T6(400))在頻點(diǎn)400 Hz附近頻段內(nèi)的預(yù)測(cè)精度較高,在遠(yuǎn)離400 Hz頻段內(nèi)的預(yù)測(cè)精度較差。其余降階模型具有類似規(guī)律。

提取表 1中5個(gè)降階模型(ROM_T20(0)、ROM_T3(100)、ROM_T4(200)、ROM_T6(400) 、ROM_T5(800))計(jì)算動(dòng)柔順度誤差的下包絡(luò),如圖3(b)中的紅色實(shí)線所示。可以看出,紅色實(shí)線完全位于水平灰色虛線的下方,表明展開(kāi)頻點(diǎn)間的預(yù)測(cè)誤差亦滿足精度要求,這也證明了本文利用交叉驗(yàn)證策略來(lái)評(píng)估降階模型的預(yù)測(cè)精度是合理有效的。

圖4進(jìn)一步討論了增加展開(kāi)頻點(diǎn)數(shù)量及正交基階數(shù)的必要性。圖4中縱軸表示動(dòng)柔順度相對(duì)誤差的最大值,橫軸表示正交基的階數(shù),紅色矩形框表示增加展開(kāi)頻點(diǎn)的位置,水平灰色虛線表示誤差閾值。

圖4 動(dòng)柔順度最大相對(duì)誤差與標(biāo)準(zhǔn)正交基數(shù)目的關(guān)系Fig.4 Maximum relative error of dynamic compliance with respect to the order of orthonormal basis

由圖4可以看出,當(dāng)僅有1個(gè)展開(kāi)頻點(diǎn)時(shí),隨著正交基階數(shù)的增加,降階模型的計(jì)算精度逐漸提高,但當(dāng)正交基數(shù)超過(guò)6時(shí),降階模型的計(jì)算精度趨于穩(wěn)定,這意味著只增加正交階數(shù)而不增加展開(kāi)頻點(diǎn)數(shù)量不能進(jìn)一步提高降階模型的計(jì)算精度;而在增加展開(kāi)頻點(diǎn)后,降階模型的計(jì)算精度迅速提升,并且隨著展開(kāi)頻點(diǎn)數(shù)目的增加,降階模型的計(jì)算精度越來(lái)越高。當(dāng)展開(kāi)頻點(diǎn)數(shù)目增加到5個(gè)時(shí),降階模型的最大動(dòng)柔順度誤差最大值為1×10-7。因此,僅在單個(gè)展開(kāi)頻點(diǎn)生成多個(gè)標(biāo)準(zhǔn)正交基不能滿足對(duì)降階模型計(jì)算精度的要求,本文提出的自適應(yīng)增加展開(kāi)頻點(diǎn)和標(biāo)準(zhǔn)正交基的策略可建立滿足精度要求的降階模型。

3.2 網(wǎng)格加筋筒殼的諧響應(yīng)分析

本節(jié)利用網(wǎng)格加筋筒結(jié)構(gòu)進(jìn)一步驗(yàn)證自適應(yīng)模型降階方法的有效性。結(jié)構(gòu)模型如圖5所示,幾何尺寸如下:圓柱筒直徑D=5000 mm,長(zhǎng)度L=16 000 mm,蒙皮厚度ts=5 mm,環(huán)筋數(shù)目21,縱筋數(shù)目30,筋條高度100 mm,筋條厚度tc=5 mm。蒙皮材料彈性模量E=71 GPa,筋條材料彈性模量為E=50 GPa,密度ρ=2700 kg/m3,泊松比μ=0.3,阻尼比為ξ=5%。邊界條件為一端固支一端自由,在自由側(cè)施加簡(jiǎn)諧激勵(lì),頻率范圍[0, 100 Hz],分析步長(zhǎng)為1 Hz。

圖5 正置正交網(wǎng)格加筋筒結(jié)構(gòu)Fig.5 Orthogonal grid stiffened cylinder

選取初始標(biāo)準(zhǔn)正交基的階數(shù)r=3、步長(zhǎng)p=1,最大正交基數(shù)目Maxorder=20。根據(jù)圖1給出的自適應(yīng)模型降階流程,在目標(biāo)頻段內(nèi)確定了3個(gè)展開(kāi)頻點(diǎn)以及每個(gè)頻點(diǎn)處所需標(biāo)準(zhǔn)正交基的數(shù)目,如表 2所示。將3個(gè)展開(kāi)頻點(diǎn)生成的標(biāo)準(zhǔn)正交基進(jìn)行SVD正交化得到降階模型的減縮基矩陣T。利用降階模型計(jì)算結(jié)構(gòu)的動(dòng)柔順度及相對(duì)誤差如圖6所示。

表2 自適應(yīng)確定的展開(kāi)頻點(diǎn)和標(biāo)準(zhǔn)正交基階數(shù)Table 2 Expansion frequency points and its related orders of orthonormal basis Frequency points

(a) Dynamic compliance

(b) Relative error圖6 網(wǎng)格加筋筒結(jié)構(gòu)的動(dòng)柔順度及誤差Fig.6 Dynamic compliance and relative error of the orthogonal grid stiffened cylinder

圖6中給出了選取結(jié)構(gòu)前100階模態(tài)時(shí)模態(tài)疊加法(Modal Superposition Method,MSM)和模態(tài)加速度法(Modal Acceleration Method,MAM)的計(jì)算結(jié)果。由圖6可以看出,本文建立的降階模型計(jì)算結(jié)果與直接頻響法十分吻合,最大相對(duì)誤差為1×10-7;模態(tài)疊加法和模態(tài)加速度法在低頻段內(nèi)計(jì)算精度能夠滿足要求,但在高頻段內(nèi)計(jì)算誤差較大,模態(tài)疊加法的在高頻段內(nèi)預(yù)測(cè)誤差最大值為421.97%,模態(tài)加速度法的預(yù)測(cè)誤差最大值為46.63%。這是由于本文中的網(wǎng)格加筋筒殼結(jié)構(gòu)是由兩種材料組成的非比例阻尼結(jié)構(gòu),其阻尼矩陣不再滿足模態(tài)正交化條件,導(dǎo)致模態(tài)疊加法和模態(tài)加速度法的求解精度降低。

表 3統(tǒng)計(jì)了直接頻響法、模態(tài)疊加法法、模態(tài)加速度法和本文提出的自適應(yīng)模型降階法的計(jì)算時(shí)間,可以看出,自適應(yīng)模型降階法的計(jì)算時(shí)間為25.80 s,分別是直接頻響法計(jì)算時(shí)間的0.78%,模態(tài)疊加法計(jì)算時(shí)間的37.52%,模態(tài)加速度法計(jì)算時(shí)間的35.47%。綜上可以證明本文提出的自適應(yīng)降階模型表現(xiàn)出更優(yōu)異的計(jì)算精度和效率。

表3 不同方法計(jì)算頻響函數(shù)的時(shí)間Table 3 Computational time to calculate frequency response function by different methods

4 結(jié)論

(1)針對(duì)非比例阻尼結(jié)構(gòu)諧響應(yīng)分析效率低的問(wèn)題,本文闡述了基于SOAR方法建立降階模型的主要流程以及研究確定展開(kāi)頻點(diǎn)數(shù)目和標(biāo)準(zhǔn)正交基階數(shù)的必要性,提出了基于二分法和交叉驗(yàn)證的自適應(yīng)降階策略來(lái)提高非比例阻尼結(jié)構(gòu)諧響應(yīng)計(jì)算精度和效率。

(2)通過(guò)阻尼涂層板和網(wǎng)格加筋筒殼數(shù)值算例驗(yàn)證了提出的自適應(yīng)降階策略的有效性。數(shù)值結(jié)果表明,自適應(yīng)降階模型預(yù)測(cè)結(jié)構(gòu)動(dòng)柔順度相對(duì)誤差最大值不超過(guò)1×10-7,相比模態(tài)疊加法和模態(tài)加速度法表現(xiàn)出更高的計(jì)算精度和效率。

(3)本文建立自適應(yīng)模型降階方法適用于由不同材料組成的非比例阻尼結(jié)構(gòu)快速諧響應(yīng)分析,后續(xù)研究中,將進(jìn)一步擴(kuò)展所提出方法在其他類型的非比例阻尼結(jié)構(gòu)諧響應(yīng)分析的適用性。