周 航， 校金友， 徐 超

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

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

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

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

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

1 黏彈性阻尼結(jié)構(gòu)Layerwise板單元建模方法

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

圖1 典型黏彈性阻尼結(jié)構(gòu)Fig.1 A typical viscoelastic damping structure

針對(duì)阻尼層采用實(shí)體單元建模方法存在的問題，徐超等人推導(dǎo)了一種基于Layerwise離散層理論的四節(jié)點(diǎn)Layerwise板單元用于黏彈性阻尼結(jié)構(gòu)的有限元建模，通過采用統(tǒng)一的位移場(chǎng)，確保了各層之間位移的連續(xù)性，考慮了每一層的剪切變形，具有建模簡(jiǎn)單、計(jì)算耗費(fèi)小和便于參數(shù)優(yōu)化的特點(diǎn)。

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

根據(jù)Layerwise板單元理論可得多層結(jié)構(gòu)的自由振動(dòng)微分方程：

(1)

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

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

(2)

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

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

(3)

本文僅考慮黏彈性阻尼材料的阻尼而忽略其他形式的阻尼。由自由振動(dòng)方程可得到黏彈性阻尼結(jié)構(gòu)的自由振動(dòng)方程：

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

(4)

式中：λ和V分別為特征值和對(duì)應(yīng)的特征向量。

表1 27℃下3M ISD112材料性能參數(shù)

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

T(λ)V=0

(5)

當(dāng)系數(shù)矩陣T(ω)為常值矩陣或ω的線性函數(shù)時(shí)，上式即為有限元分析中常見的標(biāo)準(zhǔn)特征值問題和廣義特征值問題。而在上述3MISD112黏彈性阻尼材料結(jié)構(gòu)中，T(ω)=KE+G′(ω)KV-ω2M。

下面介紹一種新型穩(wěn)定的、求解非線性特征值問題的數(shù)值方法RSRR。

2 基于預(yù)解矩陣采樣法的Rayleigh-Ritz方法

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

TQ(λ)g=0

(6)

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

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

2.1 預(yù)解矩陣采樣法

在RSRR方法中，特征空間的構(gòu)造是通過預(yù)解矩陣采樣法來實(shí)現(xiàn)的，即通過對(duì)系統(tǒng)矩陣逆矩陣的采樣來構(gòu)造可靠的特征空間。構(gòu)造特征空間之前，首先需要確定特征值的搜索區(qū)域，二維區(qū)域可以由橢圓或者矩形圍道來定義。常用的橢圓形圍道可以很容易的由以下參數(shù)化方程定義：

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

(7)

預(yù)解矩陣采樣法主要包括以下兩個(gè)步驟：

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

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

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

(8)

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

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

T(z)-1=

(9)

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

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

(10)

于是有：

(11)

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

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

(12)

式中：Rj為展開系數(shù)矩陣。當(dāng)特征值搜索區(qū)域?yàn)橐粋€(gè)區(qū)間時(shí)，通常可取gj(z)為該區(qū)間上的Chebyshev多項(xiàng)式，此時(shí)展開式(12)收斂很快。對(duì)于其他形式的區(qū)域，也可以通過差值等途徑來構(gòu)造gj(z)。在本文中并不要求顯式地計(jì)算gj(z)。

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

(13)

根據(jù)式(12)可得：

(14)

(15)

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

(16)

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

(17)

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

(18)

(19)

因此，矩陣S的列張成了一個(gè)有效的特征空間。實(shí)際Rayleigh-Ritz投影中，需要該列空間的正交基Q，可以通過對(duì)矩陣S的截?cái)嗥娈愔捣纸猥@得，本文設(shè)奇異值截?cái)嗾`差為δ。由于S的列數(shù)NL通常是一個(gè)比較小的數(shù)值，且和問題的規(guī)模n無關(guān)，奇異值分解的計(jì)算耗費(fèi)隨n線性變化。令k為S的數(shù)值秩，即Q的列數(shù)。為了準(zhǔn)確提取圍道C所包圍區(qū)域內(nèi)的所有特征值，必須滿足k≥rank(VC)。這里應(yīng)注意，為提高方法的穩(wěn)定性，在對(duì)S的奇異值分解之前，對(duì)其中的列向量進(jìn)行歸一化處理很重要。

2.2 基于預(yù)解矩陣采樣法的Rayleigh-Ritz方法

表2給出了本文非線性特征值解法RSRR算法的基本步驟。RSRR屬于標(biāo)準(zhǔn)的Rayleigh-Ritz投影方法，其中的特征空間由預(yù)解矩陣采樣法構(gòu)造。

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

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

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

3 算 例

本文黏彈性阻尼層板結(jié)構(gòu)模態(tài)分析方法的求解過程分為三步：

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

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

(3) 應(yīng)用RSRR算法求解非線性特征值問題，獲得結(jié)構(gòu)模態(tài)頻率和振型。

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

3.1 阻尼復(fù)合材料慣導(dǎo)設(shè)備支架

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

圖2 阻尼復(fù)合材料慣導(dǎo)設(shè)備支架CAD模型Fig.2 CAD model of the damping composite inertial navigation equipment support

圖3 阻尼復(fù)合材料慣導(dǎo)設(shè)備支架Layerwise板單元模型Fig.3 Layerwsie plate element model of the damping composite inertial navigation equipment support

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

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

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

表3 慣導(dǎo)設(shè)備支架前10階模態(tài)參數(shù)及特征值相對(duì)誤差

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

3.2 阻尼復(fù)合材料紅外設(shè)備支架

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

圖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 阻尼復(fù)合材料紅外設(shè)備支架CAD模型Fig.5 CAD model of the damping composite infrared equipment support

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

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

RSRR算法求得了11個(gè)特征值，相應(yīng)的模態(tài)參數(shù)和特征對(duì)相對(duì)誤差見表 4?？梢?，該結(jié)構(gòu)存在比較接近的特征值，設(shè)置隨機(jī)矩陣列數(shù)為2是合理的；特征對(duì)的最大相對(duì)誤差僅為2.75×10-6，證明了本文非線性特征值求解算法的精度是比較高的。

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

表4 阻尼復(fù)合材料紅外設(shè)備支架前11階模態(tài)參數(shù)

圖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 結(jié) 論

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