考慮頻變阻尼的黏彈性阻尼層板結構模態分析方法

周 航， 校金友， 徐 超

(1.中國航天科技集團公司第四研究院第41研究所，西安 710025；2.西北工業大學 航天學院，西安 710072)

黏彈性阻尼減振技術已被廣泛應用于航空、航天、汽車等工程領域，是解決振動控制問題的一項經濟而有效的途徑[1-2]。黏彈性阻尼結構主要通過黏彈性阻尼層的剪切變形耗散能量，以達到減振降噪的目的，因此在結構有限元建模和分析時必須能夠準確反映出這一特性[3-5]。目前黏彈性阻尼結構的有限元建模主要有兩種思路：一是對黏彈性阻尼層采用實體單元進行建模，另一種是基于離散多層板殼理論的Layerwise板殼單元建模法。采用后者可以大大簡化黏彈性阻尼結構的建模過程，同時也便于參數優化[6]。

黏彈性阻尼結構有限元分析中，黏彈性阻尼材料的頻變特性將導致結構模態分析是求解一個非線性特征值問題。目前工程實際中的普遍做法是在某一頻段內假設材料的剪切模量和阻尼參數為常值，不考慮頻變特性，從而得到一個常規二次特征值問題，利用現有Lanczos或子空間迭代方法求解，或者考慮頻變特性，采用迭代法求解非線性特征值問題。不考慮頻變特性必然帶來分析結果的不準確性，而現有非線性特征值問題的迭代解法則可能存在結果不收斂或者大規模問題計算耗費過大的問題[7]。

圍道積分方法是過去十幾年中發展的一種求解特征值問題的新方法[8-11]，相比其他方法，它不僅可以同時求解得到特征值搜索區域的所有特征對，而且適合于大規模并行計算。但是圍道積分法求解特征值時要計算形如T(z)-1U的矩陣的高階多項式矩，并基于此形成塊Hankel矩陣。當矩的階數K較高時，所形成的Hankel矩陣形態變差，甚至接近奇異，限制了圍道積分法對大規模問題的求解精度[12]。為克服此問題，Yokota等提出了投影型圍道積分法SSRR。它利用圍道積分法構造特征向量搜索空間，而后采用Rayleigh-Ritz投影法求解非線性特征值問題。對于矩陣的高階T(z)-1U矩條件數較差的情況，SSRR算法一定程度上改善了特征值解法的穩定性和精度，但是對于高階矩接近奇異的情況，SSRR仍無能為力，甚至會出現漏根[13]。

本文基于上述研究進展，首先提出了一種新的特征空間構造方法，稱為預解矩陣采樣法(Resolvent Sampling Method)。它完全避免了計算T(z)-1U矩陣的矩，從而有效解決了現有圍道積分法和SSRR方法中由于高階矩所帶來的缺陷，顯著改善了求解精度和穩定性。其次，將預解矩陣采樣法和Rayleigh-Ritz投影法相結合，建立了一種新的非線性特征值數值解法RSRR，并將其應用于黏彈性阻尼結構的大規模模態分析，考慮了黏彈性阻尼材料的頻變特性，形成了一種通用的黏彈性阻尼結構模態分析方法。通過工程實際算例證明了RSRR方法的正確性和求解效率。

1 黏彈性阻尼結構Layerwise板單元建模方法

一般而言，典型的黏彈性阻尼結構為三明治夾層結構，即黏彈性阻尼層嵌入在兩層理想彈性層之間，如圖 1所示。由于黏彈性阻尼結構主要通過阻尼層的剪切變形耗散能量，從而達到減振降噪的目的，基于經典層合板理論和一階剪切變形的有限元板殼單元不能反映出黏彈性阻尼結構的此類特性，因此采用商業有限元軟件分析黏彈性阻尼結構動態特性時，阻尼層的建模必須采用實體單元。實體單元建模較板殼單元建模復雜，不利于參數優化，而且阻尼層厚度較小，導致阻尼層單元的長厚比較大，影響計算結果的準確性。

圖1 典型黏彈性阻尼結構Fig.1 A typical viscoelastic damping structure

針對阻尼層采用實體單元建模方法存在的問題，徐超等人推導了一種基于Layerwise離散層理論的四節點Layerwise板單元用于黏彈性阻尼結構的有限元建模，通過采用統一的位移場，確保了各層之間位移的連續性，考慮了每一層的剪切變形，具有建模簡單、計算耗費小和便于參數優化的特點。

在厚度方向上，Layerwise離散層理論將多層結構分為層，每一層可以是一個單層也可以包含多個單層。這層中的每一層均視為考慮面內彎曲和橫向剪切變形的Mindlin板，不考慮厚度方向上的變形，各層之間滿足位移連續條件，并且各層材料均勻，滿足線彈性和線黏彈性理論。

根據Layerwise板單元理論可得多層結構的自由振動微分方程：

(1)

在黏彈性阻尼結構中，剛度矩陣K為復數矩陣，且同時是頻率ω的函數，即K=K(ω)。一般情況下，K(ω)是由線彈性剛度矩陣Kε和黏彈性剛度矩陣KV兩部分組成的，KE和EV均為實常數對稱矩陣。

假設黏彈性阻尼材料為3M ISD112[14]，其剪切模量隨頻率的變化關系為：

(2)

在27℃下，式中的各個參數取值見表1。3M ISD112材料的泊松比和密度分別為0.5和1 600 kg/m3，G0=0.5 MPa。此時，黏彈性阻尼結構的剛度矩陣K(ω)可表示為：

K(ω)=KE+G′(ω)KV

(3)

本文僅考慮黏彈性阻尼材料的阻尼而忽略其他形式的阻尼。由自由振動方程可得到黏彈性阻尼結構的自由振動方程：

[KE+G′(λ)KV-λ2M]V=0

(4)

式中：λ和V分別為特征值和對應的特征向量。

表1 27℃下3M ISD112材料性能參數

特征值問題的系數矩陣是頻率ω的非線性函數，數學上稱為非線性特征值問題，具有如下一般形式：

T(λ)V=0

(5)

當系數矩陣T(ω)為常值矩陣或ω的線性函數時，上式即為有限元分析中常見的標準特征值問題和廣義特征值問題。而在上述3MISD112黏彈性阻尼材料結構中，T(ω)=KE+G′(ω)KV-ω2M。

下面介紹一種新型穩定的、求解非線性特征值問題的數值方法RSRR。

2 基于預解矩陣采樣法的Rayleigh-Ritz方法

本文的非線性特征值問題解法RSRR本質上是一種Rayleigh-Ritz投影法，主要包括兩步：首先，生成一個特征向量搜索子空間，此空間應包含所有待求的特征向量，令Q∈Ck×n為該子空間的一組正交基；然后，將原特征值問題投影到Q空間上，得到一個縮減后的小規模非線性特征值問題：

TQ(λ)g=0

(6)

式中：TQ(Z)=QHT(Z)Q∈Ck×k。采用現有方法求解小型特征值問題，所得的特征值和原問題相同，特征向量滿足關系V=Qg。

RSRR算法的創新性主要在于提出了一個新的而且有效的生成特征空間的方法，即預解矩陣采樣法。

2.1 預解矩陣采樣法

在RSRR方法中，特征空間的構造是通過預解矩陣采樣法來實現的，即通過對系統矩陣逆矩陣的采樣來構造可靠的特征空間。構造特征空間之前，首先需要確定特征值的搜索區域，二維區域可以由橢圓或者矩形圍道來定義。常用的橢圓形圍道可以很容易的由以下參數化方程定義：

φ(α)=γ+[αcos(α)+ibsin(α)],α∈[0,2π)

(7)

預解矩陣采樣法主要包括以下兩個步驟：

(1)在圍道C上或者其內部選取N個采樣點zi,i=1,2,…,N，并生成一個隨機矩陣U∈Cn×L；

(2)計算T(zi)-1U，并組裝矩陣S∈Cn∈N·L，矩陣S的表達式為：

S=[T(z1)-1U,T(z2)-1U,…,T(zN)-1U]

(8)

span(S)即為本文所構造的特征空間。

下面說明span(S)是一個有效的特征空間。系統矩陣T(z)為z的矩陣函數，假設T(z)在所有特征值附近是正則的，大多數實際問題均滿足此假設，那么系統矩陣的逆矩陣T(z)-1可表示為：

T(z)-1=

(9)

式(9)表明：矩陣T(z)-1的特征值即為其逆矩陣T(z)-1的極點；T(z)-1的奇異部分可以由矩陣T(z)和其共軛轉置矩陣T(z)H的廣義特征向量表示。

式(9)的右端第一項可以寫成矩陣形式：

(10)

于是有：

(11)

式中：Φ(z)是矩陣T(z)特征值的Jordan塊組成的矩陣，VC和WC分別是由矩陣T(z)和T(z)H的廣義特征向量組成的矩陣，本文構造特征空間的主要目標就是提取VC的列空間。

為了從T(z)-1中準確提取VC的列空間信息，必須消除RC(z)的影響。由于RC(z)是解析函數，對于一組確定的基函數gj(z)，可將其表示為：

(12)

式中：Rj為展開系數矩陣。當特征值搜索區域為一個區間時，通?？扇j(z)為該區間上的Chebyshev多項式，此時展開式(12)收斂很快。對于其他形式的區域，也可以通過差值等途徑來構造gj(z)。在本文中并不要求顯式地計算gj(z)。

使用隨機矩陣U的目的是從T(z)-1中探測特征向量的信息，并壓縮RC(z)的列空間。U包含L個線性獨立的列向量。為了準確的提取所有特征值，包括重根，L至少要等于矩陣T在圍道C所包圍區域內特征值的最大代數重數[10-11]，即：

(13)

根據式(12)可得：

(14)

(15)

矩陣Z(z)的表達式為：

(16)

式中：Dj(z)=diag(gj(z),gj(z),…,gj(z))∈CL×L。

(17)

通過選擇采樣點zi,i=1,2,…,N，由式(8)計算矩陣S，矩陣S可以寫成如下的形式：

(18)

(19)

因此，矩陣S的列張成了一個有效的特征空間。實際Rayleigh-Ritz投影中，需要該列空間的正交基Q，可以通過對矩陣S的截斷奇異值分解獲得，本文設奇異值截斷誤差為δ。由于S的列數NL通常是一個比較小的數值，且和問題的規模n無關，奇異值分解的計算耗費隨n線性變化。令k為S的數值秩，即Q的列數。為了準確提取圍道C所包圍區域內的所有特征值，必須滿足k≥rank(VC)。這里應注意，為提高方法的穩定性，在對S的奇異值分解之前，對其中的列向量進行歸一化處理很重要。

2.2 基于預解矩陣采樣法的Rayleigh-Ritz方法

表2給出了本文非線性特征值解法RSRR算法的基本步驟。RSRR屬于標準的Rayleigh-Ritz投影方法，其中的特征空間由預解矩陣采樣法構造。

RSRR可以很容易的和其它軟件聯合進行計算，如有限元軟件ANSYS等。實際上，只有步驟(2)和步驟(4)包含對矩陣T的操作，需要在軟件中執行，其余的步驟都可以看作是前處理和后處理過程。

RSRR算法的主要計算量在于第2步求解N個線性方程組T(zi)-1U。這N個方程組以及每個方程組的個右端向量都是相互獨立的，可以在若干個計算核上并行求解，這樣有望大幅提高本方法的計算效率，這也是RSRR算法后續研究的一個方向。

表2 基于預解矩陣采樣法的Rayleigh-Ritz方法

3 算 例

本文黏彈性阻尼層板結構模態分析方法的求解過程分為三步：

(1) 采用Layerwise板單元建模程序建立黏彈性阻尼結構的有限元模型；

(2) 生成質量矩陣M、線彈性剛度矩陣KE和黏彈性剛度矩陣KV，并組裝系統矩陣T；

(3) 應用RSRR算法求解非線性特征值問題，獲得結構模態頻率和振型。

本文的黏彈性阻尼結構Layerwise板單元建模程序和RSRR算法均采用MATLAB程序設計語言實現。所有計算在一臺四核Intel Core i5-3350P CPU(3.10 GHz)的計算機上進行。為了考察算法的求解精度，定義特征對的相對誤差為‖T(λ)v‖2/‖v‖2。對于非線性特征值求解方法誤差的考量，標準并不唯一，‖T(λ)v‖2/‖v‖2是Asakura和Yokota等在文獻中應用較多的一種。因此，本文采用‖T(λ)v‖2/‖v‖2考察算法的求解精度。

3.1 阻尼復合材料慣導設備支架

如圖2所示的阻尼復合材料慣導設備支架，在面板的中部以及底板的中部嵌有0.2 mm厚的黏彈性阻尼材料3M ISD112。慣導設備支架通過外部螺栓與飛機艙體連接，通過內部螺栓與慣導設備連接。

圖2 阻尼復合材料慣導設備支架CAD模型Fig.2 CAD model of the damping composite inertial navigation equipment support

圖3 阻尼復合材料慣導設備支架Layerwise板單元模型Fig.3 Layerwsie plate element model of the damping composite inertial navigation equipment support

圖 3為慣導設備支架的Layerwise板單元模型，在建立模型的過程中沒有考慮螺栓孔，位移邊界條件為外部螺栓孔附件的節點固定約束，整個模型的自由度為41 620。欲求解[0,600]Hz頻率范圍之內的特征值，在式(7)定義的橢圓圍道上取100個采樣點，隨機矩陣列數設置為L=2，橢圓圍道的幾何參數設置為：γ=300，a=300，b=0.4a。

利用RSRR算法在[0,600]Hz范圍之內求得了10個特征值，表 3給出了這10階模態的參數以及相應的特征對的相對誤差?？梢钥闯?，慣導設備支架的前4階模態頻率比較接近，第6階和第7階模態頻率比較接近，第8、9、10階的模態頻率也很接近，慣導設備支架的模態較為密集，而特征對的最大相對誤差僅為1.12×10-6，由此證明本文的RSRR算法能夠用于模態密集結構的非線性特征值求解，精度比較高。

同時，通過本算例比較RSRR和文獻[7]中SSRR算法的精度和穩定性。SSRR算法計算中T(z)-1U矩陣矩的階數取為K=100，此時SSRR所構造的特征空間最大，精度也達到上限。其余參數均和RSRR算法相同。表 3最后兩列為RSRR和SSRR所得到的慣導設備阻尼復合材料支架特征對的誤差，顯然，RSRR算法的求解精度高于SSRR算法。

表3 慣導設備支架前10階模態參數及特征值相對誤差

RSRR較高的求解精度源于其特征空間矩陣較好的性態。圖5顯示了SSRR的特征空間矩陣和RSRR特征空間矩陣奇異值的衰減情況。為方便比較，圖中對所有奇異值進行歸一化，使兩個矩陣的第一個奇異值均為1?？梢?，SSRR奇異值的衰減速度明顯快于RSRR，這說明SSRR的特征空間矩陣的性態隨著特征值數目的增加迅速惡化，由此導致所構造的特征空間誤差較大。相比之下，RSRR的特征空間矩陣則性態較好，因此特征空間穩定性好，求解精度也較高。

3.2 阻尼復合材料紅外設備支架

圖5為一阻尼復合材料紅外設備支架的CAD模型，在支架的面板、側板以及底板中部分別嵌入厚度為0.2 mm的黏彈性阻尼材料3M ISD112，支架通過外部螺栓孔與載機艙體連接，內部螺栓孔與紅外設備連接。

圖4 算例3.1中RSRR和SSRR算法特征空間矩陣奇異值比較Fig.4 The singular value behavior of the eigenspace matrices in the RSRR and SSRR algorithms in example 3.1

圖5 阻尼復合材料紅外設備支架CAD模型Fig.5 CAD model of the damping composite infrared equipment support

圖6 阻尼復合材料紅外設備支架Layerwise板單元模型Fig.6 Layerwise plate element model of the damping composite infrared equipment support

圖 6為支架的Layerwise板單元模型，在建模過程中不考慮螺栓孔，位移邊界條件為外部螺栓孔附件的節點固定約束，整個模型的自由度約為24 220。欲求解結構在[0,1 000]Hz頻率范圍之內的特征值，直接在橢圓圍道上取100個采樣點，隨機矩陣列數設置為L=2，橢圓圍道的幾何參數設置為：γ=500，a=300，b=0.4a。

RSRR算法求得了11個特征值，相應的模態參數和特征對相對誤差見表 4?？梢?，該結構存在比較接近的特征值，設置隨機矩陣列數為2是合理的；特征對的最大相對誤差僅為2.75×10-6，證明了本文非線性特征值求解算法的精度是比較高的。

該算例中也進行了RSRR與SSRR方法計算精度和穩定性方面的對比。SSRR方法中K=100，其余參數與RSRR相同。表 4最后兩列給出了RSRR和SSRR方法獲得的特征對誤差，RSRR的求解精度仍明顯高于SSRR。圖 7為兩種方法特征空間矩陣奇異值的對比情況，可見該算例中SSRR的奇異值衰減速度仍明顯快于RSRR，因此它的特征空間穩定性和求解精度不及RSRR算法。

表4 阻尼復合材料紅外設備支架前11階模態參數

圖7 算例3.2中RSRR和SSRR算法特征空間矩陣奇異值比較Fig.7 The singular value behavior of the eigenspace matrices in the RSRR and SSRR algorithms in example 3.2

4 結 論

本文采用基于Layerwise離散層理論的Layerwise板單元建立黏彈性阻尼結構的有限元模型，比采用混合單元建模方法簡單方便，同時研究了一種求解非線性特征值問題的RSRR方法，結合Layerwise板單元建模方法，將RSRR拓展應用到了黏彈性阻尼結構的模態分析當中。通過兩個工程實際算例，證明了本文的RSRR在實際應用中具有較高的求解精度和好的穩定性。本文工作為考慮阻尼頻變特性的黏彈性阻尼結構模態分析，提供了有效的數值方法。