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

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

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

0 引言

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

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

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

1 諧響應分析方法

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

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

(1)

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

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

(2)

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

1.1 模態疊加法

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

假定結構特征值和模態分別是Λ和Φ:

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

(3)

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

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

(4)

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

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

(5)

式中ξ為瑞利阻尼比。

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

(6)

可以看出,模態疊加法通過模態的線性組合來近似結構的位移,避免了動剛度陣的分解,提高了計算效率。實用中,模態疊加法一般是通過模態截斷選取有限數目的模態近似表示結構的位移:

(7)

式中l為選取的模態數,l<

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

1.2 模態加速度法

對于無阻尼結構,式(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)

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

(10)

可以看出,模態加速度法通過引入擬靜力響應來補償截斷模態的貢獻,從而達到用較少模態獲得高精度響應的效果。

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 阻尼涂層板的諧響應分析

本節采用阻尼涂層板驗證自適應降階策略有效性,阻尼涂層板結構長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。基體結構右端固支,左端施加簡諧激勵,激勵頻率的范圍為[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頻段內的預測精度較差。其余降階模型具有類似規律。

提取表 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時,降階模型的計算精度趨于穩定,這意味著只增加正交階數而不增加展開頻點數量不能進一步提高降階模型的計算精度;而在增加展開頻點后,降階模型的計算精度迅速提升,并且隨著展開頻點數目的增加,降階模型的計算精度越來越高。當展開頻點數目增加到5個時,降階模型的最大動柔順度誤差最大值為1×10-7。因此,僅在單個展開頻點生成多個標準正交基不能滿足對降階模型計算精度的要求,本文提出的自適應增加展開頻點和標準正交基的策略可建立滿足精度要求的降階模型。

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

本節利用網格加筋筒結構進一步驗證自適應模型降階方法的有效性。結構模型如圖5所示,幾何尺寸如下:圓柱筒直徑D=5000 mm,長度L=16 000 mm,蒙皮厚度ts=5 mm,環筋數目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階模態時模態疊加法(Modal Superposition Method,MSM)和模態加速度法(Modal Acceleration Method,MAM)的計算結果。由圖6可以看出,本文建立的降階模型計算結果與直接頻響法十分吻合,最大相對誤差為1×10-7;模態疊加法和模態加速度法在低頻段內計算精度能夠滿足要求,但在高頻段內計算誤差較大,模態疊加法的在高頻段內預測誤差最大值為421.97%,模態加速度法的預測誤差最大值為46.63%。這是由于本文中的網格加筋筒殼結構是由兩種材料組成的非比例阻尼結構,其阻尼矩陣不再滿足模態正交化條件,導致模態疊加法和模態加速度法的求解精度降低。

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

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

4 結論

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

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

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