非比例阻尼結構諧響應分析的自適應模型降階方法①

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

(1.中國航發(fā)沈陽發(fā)動機研究所,沈陽 110015;2. 大連理工大學 工程力學系,大連 116024)

0 引言

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

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

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

1 諧響應分析方法

簡諧激勵載荷下結構的頻響方程為

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

(1)

直接頻響法將一系列離散頻點直接代入式(1)得到結構在頻點ω處的響應:

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

(2)

可以看出,直接頻響法需要得到每個頻點ω處動剛度陣的逆矩陣,計算規(guī)模較大,不適用于激勵頻率密集或自由度數目較多的問題。

1.1 模態(tài)疊加法

模態(tài)疊加法是把對應無阻尼結構的模態(tài)為空間基底,經過坐標變換使得原動力學方程解耦,通過求解n個獨立的方程獲得模態(tài)位移,進而通過疊加各階模態(tài)的貢獻求得結構的響應。

假定結構特征值和模態(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)疊加法通過模態(tài)的線性組合來近似結構的位移,避免了動剛度陣的分解,提高了計算效率。實用中,模態(tài)疊加法一般是通過模態(tài)截斷選取有限數目的模態(tài)近似表示結構的位移:

(7)

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

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

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

對于無阻尼結構,式(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)加速度法通過引入擬靜力響應來補償截斷模態(tài)的貢獻,從而達到用較少模態(tài)獲得高精度響應的效果。

2 基于SOAR的自適應模型降階方法

2.1 SOAR算法

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

(11)

其中,

(12)

通過二階Krylov子空間對式(12)進行降階,對應減縮基矩陣如下:

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

(13)

其中,r0(ωa),r1(ωa),…,rk-1(ωa)代表在頻點ωa處展開的k個標準正交基向量:

(14)

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

(15)

其中,

(16)

可以看出,基于二階Krylov子空間建立降階模型與展開頻點數及標準正交基階數有關,展開頻點數和標準正交基階數越多,降階模型的計算精度越高。為平衡降階模型的計算精度和效率,本文提出了自適應確定展開頻點數和標準正交基階數的降階策略。

2.2 自適應降階策略

結構動柔順度[17]誤差常被用于衡量降階模型的計算精度:

(17)

式(17)需要獲得所有頻點的位移來計算結構的動柔順度誤差,這使得建立降階模型過程中需要浪費大量的時間去評估降階模型的計算精度,失去了利用降階模型提高計算效率的意義。本文提出利用交叉驗證的策略來評估降階模型的預測精度:

(18)

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

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

第一步:設置每個頻點的初始標準正交基階數r和最大正交基數目Maxorder,正交基階數增加的步長p;設置結構動柔順度誤差閾值Tol;初始化集合S1={ω1}和集合S2={ω2},其中ω1和ω2分別為目標頻段的上下限。

第二步:利用SOAR算法得到展開頻點ωa處的r個標準正交基,并組裝成減縮基矩陣Tr(ωa):

ωa=max{S1}

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

(19)

ωb=min{S2}

(20)

第四步:通過SOAR算法擴充展開頻點ωa處的標準正交基數目,并組裝成減縮基矩陣Tr+p(ωa):

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

(21)

第五步:利用減縮基矩陣Tr+p(ωa)計算展開頻點 處的動柔順度誤差:

(22)

S1=S1∪{ωc}

(23)

S2=S2∪{ωa}

S1=S1{ωa}

(24)

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

迭代結束后,將集合S2中所有頻點的標準正交基進行正交化得到減縮基矩陣T:

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

(25)

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

3 自適應模型降階方法驗證

3.1 阻尼涂層板的諧響應分析

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

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

選取初始標準正交基階數r=3、步長p=1,最大正交基數目Maxorder=20,根據圖1給出的自適應模型降階流程,在目標頻段內確定了5個展開頻點以及每個頻點處所需標準正交基的數目,如表 1所示。將5個展開頻點生成的標準正交基進行SVD正交化得到降階模型的減縮基矩陣T。

表1 自適應確定的展開頻點數和標準正交基階數Table 1 Expansion frequency points and its related orders of orthonormal basis frequency points

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

(a) Dynamic compliance

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

為了驗證本文提出的交叉驗證策略有效性,圖3 (b)分別給出了由單個展開頻點的正交基建立的降階模型預測誤差。圖中括號內的數字表示展開頻點,下標表示標準正交基的階數,例如,T4(200)表示展開頻點為200 Hz,該展開頻點處的標準正交基階數為4。圖中豎直灰色實線代表展開頻點,水平灰色虛線代表閾值Tol?？梢钥闯?降階模型(ROM_T4(200))在頻點200 Hz附近頻段內的預測精度較高,在遠離200 Hz的頻段內的預測精度較差。同樣,降階模型(ROM_T6(400))在頻點400 Hz附近頻段內的預測精度較高,在遠離400 Hz頻段內的預測精度較差。其余降階模型具有類似規(guī)律。

提取表 1中5個降階模型(ROM_T20(0)、ROM_T3(100)、ROM_T4(200)、ROM_T6(400) 、ROM_T5(800))計算動柔順度誤差的下包絡,如圖3(b)中的紅色實線所示?？梢钥闯?紅色實線完全位于水平灰色虛線的下方,表明展開頻點間的預測誤差亦滿足精度要求,這也證明了本文利用交叉驗證策略來評估降階模型的預測精度是合理有效的。

圖4進一步討論了增加展開頻點數量及正交基階數的必要性。圖4中縱軸表示動柔順度相對誤差的最大值,橫軸表示正交基的階數,紅色矩形框表示增加展開頻點的位置,水平灰色虛線表示誤差閾值。

圖4 動柔順度最大相對誤差與標準正交基數目的關系Fig.4 Maximum relative error of dynamic compliance with respect to the order of orthonormal basis

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

3.2 網格加筋筒殼的諧響應分析

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

圖5 正置正交網格加筋筒結構Fig.5 Orthogonal grid stiffened cylinder

選取初始標準正交基的階數r=3、步長p=1,最大正交基數目Maxorder=20。根據圖1給出的自適應模型降階流程,在目標頻段內確定了3個展開頻點以及每個頻點處所需標準正交基的數目,如表 2所示。將3個展開頻點生成的標準正交基進行SVD正交化得到降階模型的減縮基矩陣T。利用降階模型計算結構的動柔順度及相對誤差如圖6所示。

表2 自適應確定的展開頻點和標準正交基階數Table 2 Expansion frequency points and its related orders of orthonormal basis Frequency points

(a) Dynamic compliance

(b) Relative error圖6 網格加筋筒結構的動柔順度及誤差Fig.6 Dynamic compliance and relative error of the orthogonal grid stiffened cylinder

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

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

表3 不同方法計算頻響函數的時間Table 3 Computational time to calculate frequency response function by different methods

4 結論

(1)針對非比例阻尼結構諧響應分析效率低的問題,本文闡述了基于SOAR方法建立降階模型的主要流程以及研究確定展開頻點數目和標準正交基階數的必要性,提出了基于二分法和交叉驗證的自適應降階策略來提高非比例阻尼結構諧響應計算精度和效率。

(2)通過阻尼涂層板和網格加筋筒殼數值算例驗證了提出的自適應降階策略的有效性。數值結果表明,自適應降階模型預測結構動柔順度相對誤差最大值不超過1×10-7,相比模態(tài)疊加法和模態(tài)加速度法表現出更高的計算精度和效率。

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