張亞輝，馬永彬

(大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，大連 116023)

基于模態(tài)思想的有限元方法在進(jìn)行結(jié)構(gòu)振動分析時，在結(jié)構(gòu)一個振動波長內(nèi)需要劃分6至15甚至更多單元才能準(zhǔn)確地模擬結(jié)構(gòu)的振動,而在高頻振動下，結(jié)構(gòu)振動波長非常小，應(yīng)用有限元方法不得不采用大量的自由度來分析結(jié)構(gòu)的振動,因此，高頻振動問題需要尋求更為有效的分析方法。統(tǒng)計能量分析(SEA)作為高頻振動分析的典型方法[1]，自上世紀(jì)60年代提出以來，已經(jīng)推廣到多個領(lǐng)域并得到成功的應(yīng)用。采用SEA方法進(jìn)行高頻振動分析的計算成本極小，不過只適用于初步驗證階段。因為SEA按振動模式將結(jié)構(gòu)分為若干子系統(tǒng)，分析結(jié)果只能給出各子系統(tǒng)能量均方值，隨著對結(jié)果的需求更加精細(xì)化，還需要借助別的方法進(jìn)行輔助分析。從另一個角度來看，結(jié)構(gòu)的振動可以用波的傳播、反射以及傳遞的形式來表述[2]。這種表述尤其在中、高頻域內(nèi)有優(yōu)勢，因為它只需要很小的計算成本并且有著很高的精確性,因此，近年來也獲得結(jié)構(gòu)的彈性波屬性近成為了中、高頻振動研究領(lǐng)域的熱門。針對波導(dǎo)結(jié)構(gòu)的全頻域振動響應(yīng)分析問題，Mace等[3-4]提出了一個新方法：波有限元(WFE)方法。WFE方法結(jié)合了有限元方法和周期結(jié)構(gòu)理論，只需選取一小段結(jié)構(gòu)建立有限元模型，便能獲得波導(dǎo)結(jié)構(gòu)的波傳播屬性，進(jìn)而基于波傳播及波散射關(guān)系得到問題的解答。相比傳統(tǒng)有限元方法，WFE方法在計算成本方面有著很大的優(yōu)勢。目前，這種方法已經(jīng)成功應(yīng)用于波色散分析[5]，壓電材料分析[6]等領(lǐng)域。然而，由于應(yīng)用WFE方法時，需要建立一小段結(jié)構(gòu)有限元模型，因而模型的單元屬性、網(wǎng)格大小等參數(shù)都直接影響到分析結(jié)果的精確性，甚至造成嚴(yán)重的誤差[7]；于是，一個具體問題的WFE分析需要多次修正有限元模型才能得到滿意的結(jié)果。因此，如何精確地獲得波傳播參數(shù)是波傳播方法首要問題。

鐘萬勰等[8-10]將Hamilton體系及辛狀態(tài)空間理論應(yīng)用到彈性力學(xué)，改變了以往以半逆法為主的湊合解法來解決彈性力學(xué)問題的現(xiàn)狀，從而使得分離變量、辛本征展開等方法能夠應(yīng)用到彈性力學(xué)問題。目前已成功應(yīng)用于梁[11-12]、板[13-14]、殼[15-16]等基本結(jié)構(gòu)原件的彈性力學(xué)問題及薄板和中厚板的自由振動分析中[17-19]。將板振動的控制方程導(dǎo)入辛對偶體系，通過求解正則方程，得到的本征值與本征向量恰好分別是波傳播方法需要的波傳播參數(shù)與波形。它們均由完全理性的推導(dǎo)得出，沒有引入任何假設(shè)。對復(fù)雜邊界條件，相較于傳統(tǒng)方法只能在四邊簡支的邊界條件下才能給出Navier形式的閉合解，辛方法在任意邊界條件下均能給出辛解析解，具有顯著的優(yōu)勢。

本文基于辛對偶體系和波傳播理論，提出了一個分析薄板結(jié)構(gòu)穩(wěn)態(tài)振動響應(yīng)的新思路。首先，將板振動的控制方程導(dǎo)入辛對偶體系，求解得到本征值(波傳播參數(shù))與本征向量(波形)；然后，通過波傳播過程中的入射、反射以及傳遞關(guān)系得到各參與波的幅值，進(jìn)而求得薄板任意位置的響應(yīng)。與現(xiàn)有文獻(xiàn)中基于辛對偶體系求解板的自由振動問題，以及基于傳統(tǒng)模態(tài)思想求解振動問題不同的是：本文引入波空間的概念，結(jié)合波傳播理論來研究結(jié)構(gòu)的穩(wěn)態(tài)振動響應(yīng)。同時，本文以矩陣形式清晰地給出計算公式，在應(yīng)用時具有更強(qiáng)的通用性。與WFE方法相比，本文方法基于辛本征解，具有更好的精度與數(shù)值穩(wěn)定性；另外，考慮到通用性，相比WFE方法需要根據(jù)具體問題多次劃分網(wǎng)格并調(diào)整程序，應(yīng)用本文方法使得分析過程更為方便快捷。最后還需要注意到，本文公式的推導(dǎo)是基于矩形薄板結(jié)構(gòu)，如果激勵頻率過高則薄板假設(shè)不再適用，此時需要考慮板的剪切變形的影響[2]。

1 薄板波傳播的辛空間描述

首先，將彈性薄板振動的控制方程導(dǎo)入辛對偶體系。然后，求解辛本征值問題得到本征值與本征向量。本征值即為薄板的波傳播參數(shù)，本征向量即為各階波形。獲得波傳播參數(shù)與波形后就可以將物理空間的受迫振動問題映射到波空間進(jìn)行求解。

1.1 薄板彎曲振動問題導(dǎo)入辛對偶體系

考慮圖1所示矩形薄板，其彎曲自由振動方程為[20]

(1)

式中W(x,y,t)為板的撓度，D=E(1+ηi)h3/12(1-v2)為板的彎曲剛度,E、v為彈性模量和泊松比，ρ、h為板的密度和厚度,η,i分別為阻尼損耗因子和虛數(shù)單位。

對于穩(wěn)態(tài)振動問題，板的撓度可以寫為W(x,y,t)=Re{w(x,y)eiωt},w(x,y)為W(x,y,t)的幅值，因此振動方程可通過w(x,y)在頻域內(nèi)求解。根據(jù)力(矩)平衡可得下述方程

(2)

式中，Qx，Qy，Mx，My，Mxy分別為板橫截面單位長度上的剪力、彎矩和扭矩，正方向的規(guī)定如圖1所示，并且有如下關(guān)系式

(3)

圖1 薄板及坐標(biāo)系示意圖

定義等效剪力

(4)

令θ=?w/?y，并由式(2)、(3)和(4)可得

(5)

式(5)可寫為如下矩陣形式

(6)

記為

(7)

式中z={wθFyMy}T為狀態(tài)向量，H為哈密頓算子矩陣，(·)表示對y的導(dǎo)數(shù)。因此方程(7)的解為

z(x,y)=η(x)eμyy

(8)

式中η(x)為僅與x有關(guān)的向量，μy為y方向的波傳播參數(shù)。將式(8)代入方程(7)得到

Hη(x)=μyη(x)

(9)

1.2 辛本征問題的求解

求解方程(9)，考慮η(x)中各變量在x方向上有相同變化形式，即

η(x)=φeμxx

(10)

式中φ為與x無關(guān)的常向量，μx為x方向的波傳播參數(shù)。將其代入方程(9)，得到

(11)

式中，kb=(ρhω2/D)1/4為板的自由彎曲波波數(shù)。求解式(11)得到如下的特征方程

(12)

解得

μx1=-ik1,μx2=ik1

μx3=-ik2,μx4=ik2

(13)

φ1,2={1μyχ1χ2}T

φ3,4={1μyχ3χ4}T

(14)

式中

(15)

(16)

于是，方程(9)的通解可以寫為

(17)

式中si(i=1,2,3,4)為待定系數(shù)，記

s={s1s2s3s4}T

(18)

為基本系數(shù)向量。

至此，基于辛對偶體系給出了薄板的波形向量的一般表達(dá)式(17)。下面將通過矩形薄板的對邊邊界條件得到波傳播參數(shù)μy與基本系數(shù)向量s，進(jìn)而得到波形矩陣。

1.3 求解波傳播參數(shù)與基本系數(shù)向量

篇幅所限，這里僅給出對邊簡支與對邊固支兩種邊界條件下傳播參數(shù)μy和基本系數(shù)向量s，其他邊界條件可按相同過程推導(dǎo)得出。

對邊簡支邊界條件為

(19)

代入式(17)可得到

Ks=0

(20)

其中系數(shù)矩陣K為

(21)

基本系數(shù)向量s存在非平凡解的條件是系數(shù)矩陣行列式為零。于是，可得到波傳播參數(shù)μy滿足的超越方程

sin(k1a)sin(k2a)=0

(22)

因此波傳播參數(shù)為

(23)

由方程(20)得到基本系數(shù)向量s的一組非平凡解

(24)

對邊固支邊界條件為

(25)

得到波傳播參數(shù)μy滿足的超越方程

k1k2cos(k1a)cos(k2a)=0

(26)

以及基本系數(shù)向量s的一組非平凡解

s1=k2eik1a+ik1sin(k2a)-k2cos(k2a)

s2=-k2e-ik1a-ik1sin(k2a)+k2cos(k2a)

(27)

s3=-k1[cos(k1a)-eik2a]-ik2sin(k1a)

s4=k1[cos(k1a)-cos(k2a)+

isin(k2a)]-ik2sin(k1a)

這里需要注意兩個特殊的情況：①k1=0；②k2=0。由式(13)、(17)、(24)和(27)知這兩種情況下η(x)=0，意味著板沒有變形，也沒有內(nèi)力，是一組平凡解。雖然是超越方程(22)或(26)的解，但在后面的分析中將不予考慮。另外注意到本文研究的是穩(wěn)態(tài)受迫振動問題，激勵頻率不為零，所以不用考慮靜力分析中波傳播參數(shù)的重根問題。

通過求解超越方程(22)和(26)便可得到各階波的傳播參數(shù)μy，進(jìn)而求得基本系數(shù)向量s，再由式(17)便可得到各階波形η(x)。對邊簡支邊界可以得到解析解形式的波傳播參數(shù)，而其它邊界條件下的超越方程不能給出顯式解，此時可通過圍線積分等方法求解。從μy的超越方程可以看出，μy與-μy是成對出現(xiàn)的，意味著正向波和負(fù)向波是成對出現(xiàn)的，波形也相應(yīng)的分為正向波形和負(fù)向波形。在后面的具體應(yīng)用中，對傳播參數(shù)按照共軛辛正交性質(zhì)[8]排序，即μy1,μy2,…,μym,-μy1,-μy2,…,-μym,其中m是正向波的個數(shù)。正向波的波傳播參數(shù)μy,i滿足Re{μyi}<0,或者Re{μyi}=0,lm{μyi}<0且按虛部從小到大排序?；诟麟A波形之間的共軛辛正交關(guān)系，可對波形進(jìn)行歸一化處理。得到各階波形后，便可將物理空間下的受迫振動問題，轉(zhuǎn)換到以正、負(fù)波的波幅c+、c-所表征的波空間，即

z=A+c++A-c-=Ac

(28)

式中，波形矩陣A=[η1(x)η2(x) …η2m(x)]，波形矩陣的共軛辛正交性質(zhì)為

(29)

2 受迫振動的波傳播分析

波導(dǎo)結(jié)構(gòu)受外部作用的振動響應(yīng)分析在波空間下可分三個步驟來進(jìn)行[21]：首先，確定外激勵在作用處向兩側(cè)無窮波導(dǎo)產(chǎn)生的波的波幅，即直接激勵波的波幅；其次，計算波導(dǎo)方向上的邊界或者不連續(xù)處的波反射系數(shù)矩陣；最后，綜合直接激勵、波反射以及波傳播之間的關(guān)系，求解可得任意位置正、負(fù)向波的波幅c+和c-，進(jìn)而通過式(28)疊加得到任意位置的位移和內(nèi)力響應(yīng)。

2.1 確定直接激勵波波幅

外部作用在無窮波導(dǎo)上會向兩側(cè)產(chǎn)生直接激勵波，波幅用e+和e-表示。正、負(fù)向直接激勵波應(yīng)使結(jié)構(gòu)在外載作用位置滿足位移連續(xù)性以及力平衡條件。在波空間下，根據(jù)式(28)，上述條件可以表達(dá)為

(30)

式中，fext為外部作用向量。

利用波形矩陣的共軛辛正交性質(zhì)(29)，式(30)兩邊分別左乘ATJ2，并關(guān)于板x向積分。得到直接激勵波波幅的計算表達(dá)式

(31)

2.2 邊界反射系數(shù)矩陣

正向波v+入射到邊界上會產(chǎn)生負(fù)向波v-。v-與v+之間滿足的關(guān)系需要根據(jù)邊界條件給出。以簡支邊界條件為例，假設(shè)在y=b處簡支。相應(yīng)邊界條件的變分式為

(32)

將式(17)代入式(32)后得到

k=1,2,m

(33)

式(33)寫成矩陣形式為

(34)

式中矩陣U1和U2中的元素分別為

(35)

由式(34)可以得到

(36)

式中，R為邊界反射系數(shù)矩陣。

注意，這里的波幅右上標(biāo)的正負(fù)號分別代表的是入射波和反射波。對于對稱的另一條側(cè)邊，v-為入射波，v+為反射波。于是，右側(cè)邊界的波反射系數(shù)為

(37)

如果左側(cè)y=0處也是簡支邊界條件，則對應(yīng)的邊界波反射系數(shù)為

(38)

2.3 建立波傳播關(guān)系并求解

由式(28)知，一旦確定某個位置的波幅c+和c-，結(jié)合各波形的波傳播參數(shù)，則結(jié)構(gòu)任意位置的位移和內(nèi)力都可以得到。而波幅的求解可以根據(jù)上小節(jié)給出的波在邊界的反射系數(shù)以及本節(jié)給出的波傳播關(guān)系得到。下面通過圖2所示的一個簡單的波導(dǎo)結(jié)構(gòu)來說明波幅的求解過程。

波導(dǎo)結(jié)構(gòu)在ye位置受外部作用，此處的波幅a+和g-可以用直接激勵波e+和e-以及入射波g+和a-給出[7]，即

a+=e++g+,g-=e-+a-

(39)

由波的傳播與反射關(guān)系可知

b+=T(b-ye)a+,a-=T(b-ye)b-

b-=RRb+,d+=RLd-

(40)

式中T(y)=diag[eμy1yeμy2y… eμymy]是波傳播矩陣。由式(39)和(40)，可以推導(dǎo)得到

圖2 波導(dǎo)中的波傳播示意圖

a+=e++g+=e++T(ye)d+=e++T(ye)RLd-=

e++T(ye)RLT(ye)g-=

e++T(ye)RLT(ye)(e-+a-)

(41)

而

a-=T(b-ye)b-=T(b-ye)RRb+=

T(b-ye)RRT(b-ye)a+

(42)

將式(42)代入式(41)則有

a+=[I-T(ye)RLT(b)RRT(b-ye)]-1

[e++T(ye)RLT(ye)e-]

(43)

因此，響應(yīng)點y處的波幅為

c+=T(y-ye)a+

c-=T(b-y)RRT(b-y)c+

(44)

得到了各波的波幅，由式(28)便可得到響應(yīng)位置的位移和內(nèi)力。

3 算 例

考慮如圖1所示的四邊簡支矩形薄板[7]，板的材料屬性為E=2.0×1011，ρ=7 800，v=0.3，幾何參數(shù)為a=0.18，b=0.6，厚度h=0.001 8，單位均為國際標(biāo)準(zhǔn)單位。板的損耗因子η=0.03。橫向力F的作用位置(xe,ye)也在圖1中給出。

采用模態(tài)疊加法、波有限元法(WFE)以及本文方法分別求解了激勵點輸入導(dǎo)納的幅值以及板的動能、應(yīng)變能時間均值。

對于四邊簡支矩形板，可以通過模態(tài)疊加得到薄板穩(wěn)態(tài)振動響應(yīng)的解析解

(45)

事實上，四邊簡支矩形薄板的辛解析解經(jīng)過簡單的變換，在形式上與式(45)完全一致。然而，由于求解空間不同，本文方法只用較少截斷波形便可得到非常精確的結(jié)果。在本算例下，本文方法的計算時間4.80 s比WFE方法計算時間242.94 s也少很多。

圖3 模態(tài)疊加法中不同模態(tài)數(shù)下輸入導(dǎo)納幅值的相對誤差(|Y0|對應(yīng)1000階模態(tài))

圖8給出了上述薄板在對邊(x=0,a)固支，另一對邊(y=0,b)簡支邊界下，由本文方法、波有限元法以及ABAQUS三種方法給出的輸入點導(dǎo)納響應(yīng)幅值曲線。其中ABAQUS有限元模型的單元類型為S4R，單元個數(shù)為33400，選取2000階振型；本文方法選取56對波，WFE方法選取105對波。此時，波傳播參數(shù)超越方程的解沒有顯式表達(dá)式，可采用圍線積分方法求解?？梢钥闯霰疚姆椒ㄅcABAQUS給出結(jié)果吻合很好，而WFE在700 Hz后有些許誤差。計算效率方面，本文方法與WFE的計算時間分別為3611.23 s和853.43 s。本文方法花費了更多的計算時間，且主要花費在求解波傳播參數(shù)的超越方程。不過應(yīng)該注意到，本文方法是完全解析推導(dǎo)計算的，幾乎不需要任何前期準(zhǔn)備；而WFE方法除了計算花費時間，還要進(jìn)行有限元建模及模型的重復(fù)修正工作。因此，就分析成本來說，本文方法更有優(yōu)勢。

圖6 板的應(yīng)變能

4 結(jié) 論

本文將薄板穩(wěn)態(tài)強(qiáng)迫振動問題引入到辛對偶體系中，結(jié)合波傳播理論給出了一個求解新思路。和以往文獻(xiàn)基于辛對偶體系求解彈性力學(xué)以及自由振動問題相比，本文更加強(qiáng)調(diào)了波傳播的思想，使得辛方法在處理波導(dǎo)結(jié)構(gòu)振動問題時有更為廣闊的視野，比如可分析多個波導(dǎo)結(jié)構(gòu)的耦合振動響應(yīng)。另外，本文以矩陣形式給出了所有計算公式。一方面，使得求解過程更為清楚；另一方面，也使得方法的應(yīng)用具有更好的通用性。常規(guī)有限元方法隨著模型的增大計算成本會越來越高，對于高頻振動問題有時甚至不能求解。與此相比，本文方法并不受模型大小的影響。和波有限元方法(WFE)相比，本文對波傳播參數(shù)及波形的求解完全是理性推導(dǎo)得出，并沒有引入任何試函數(shù)。僅有的一處近似是在求解波傳播參數(shù)時作出的，因為對于復(fù)雜邊界條件不能得到波傳播參數(shù)的解析表達(dá)式。不過采用圍線積分法可以很容易得到方程的精確數(shù)值解。因此得到的波傳播參數(shù)與波形向量是辛解析解的，具有高度的精確性。總結(jié)來看，本文方法精確度高，分析成本非常低，適合推廣于其他波導(dǎo)問題中。

[1] Lyon R H，DeJong R G.Theory and application of statistical energy analysis [M].Butterworth-Heinemann，1995.

[2] Cremer L，Heckl M，Petersson B A T.Structure-borne sound: structural vibrations and sound radiation at audio frequencies [M].Springer，2005.

[3] Mace B R，Duhamel D，Brennan M J，et al.Finite element prediction of wave motion in structural waveguides [J].Journal of the Acoustical Society of America，2005，117(5): 2835-2843.

[4] Mencik J M，Ichchou M N.Multi-mode propagation and diffusion in structures through finite elements [J].European Journal of Mechanics-A/Solids，2005，24(5): 877-898.

[5] Waki Y，Mace B R，Brennan M J.Free and forced vibrations of a tyre using a wave/finite element approach [J].Journal of Sound and Vibration，2009，323(3-5): 737-756.

[6] Bareille O，Kharrat M，Zhou W，et al.Distributed piezoelectric guided-T-wave generator，design and analysis [J].Mechatronics，2012，22(5): 544-551.

[7] Waki Y，Mace B R，Brennan M J.Numerical issues concerning the wave and finite element method for free and forced vibrations of waveguides [J].Journal of Sound and Vibration，2009，327(1): 92-108.

[8] 姚偉岸，鐘萬勰.辛彈性力學(xué)[M].北京:高等教育出版社，2002.

[9] 鐘萬勰.應(yīng)用力學(xué)對偶體系[M].北京: 科學(xué)出版社，2002.

[10] 鐘萬勰.應(yīng)用力學(xué)的辛數(shù)學(xué)方法[M].北京: 高等教育出版社，2006.

[11] Leung A Y T，Mao S G.Symplectic integration of an accurate beam finite element in non-linear vibration [J].Computers & Structures，1995，54(6): 1135-1147.

[12] 馬國軍,徐新生,郭杏林.旋轉(zhuǎn)運動中彈性梁耦合振動的辛方法[J].計算力學(xué)學(xué)報，2004，21(6): 671-677.

MA Guo-jun，XU Xin-sheng，GUO Xing-lin.A symplectic method for the coupling vibration of elastic beams in the revolution system [J].Chinese Journal of Computational Mechanics，2004，21(6): 671-677.

[13] 談梅蘭,吳光,王鑫偉.矩形薄板面內(nèi)非線性分布載荷下的辛彈性力學(xué)解[J].工程力學(xué)，2008，25(10): 50-53.

TAN Mei-lan，WU Guang，WANG Xin-wei.Symplectic elasticity solutions for thin rectangular plates subjected to non-linear distributed in-plane loadings [J].Engineering Mechanics，2008，25(10): 50-53.

[14] 姚偉岸，孫貞.環(huán)扇形薄板彎曲問題的環(huán)向辛對偶求解方法[J].力學(xué)學(xué)報，2008，40(4): 557-563.

YAO Wei-an，SUN Zhen.Symplectic solutions for the bending of annular sector plate in circumferential direction [J].Chinese Journal of Theoretical and Applied Mechanics，2008，40(4): 557-563.

[15] 褚洪杰.彈性圓柱殼動力和熱屈曲中的辛方法[D].大連:大連理工大學(xué)，2009.

[16] 馬源.哈密頓體系下彈性圓柱殼的動態(tài)屈曲研究[D].大連理工大學(xué),2007.

[17] 鮑四元,鄧子辰.哈密頓體系下矩形薄板自由振動的一般解[J].動力學(xué)與控制學(xué)報，2005，3(2): 12-18.

BAO Si-yuan，DENG Zi-chen.A general solution of free vibration for rectangular thin plates in Hamilton systems [J].Journal of Dynamics and Control，2005，3(2): 12-18.

[18] 李銳.矩形板問題的Hamilton求解方法[D].大連:大連理工大學(xué)，2012.

[19] 鐘陽,李銳,田斌.矩形中厚板自由振動問題的哈密頓體系與辛幾何解法[J].動力學(xué)與控制學(xué)報，2009，7(4): 302-307.

ZHONG Yang，LI Rui，TIAN Bin.On Hamilton system and new symplectic approach for free vibration of moderately thick rectangular plates[J].Journal of Dynamics and Control，2009，7(4): 302-307.

[20] Leissa A W.The free vibration of rectangular plates [J].Journal of Sound and Vibration，1973，31(3): 257-293.

[21] Renno J M，Mace B R.Calculation of reflection and transmission coefficients of joints using a hybrid finite element/wave and finite element approach [J].Journal of Sound and Vibration，2013，332(9): 2149-2164.