薄板振動分析的辛空間波傳播方法

張亞輝，馬永彬

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

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

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

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

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

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

1.1 薄板彎曲振動問題導入辛對偶體系

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

(1)

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

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

(2)

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

(3)

圖1 薄板及坐標系示意圖

定義等效剪力

(4)

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

(5)

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

(6)

記為

(7)

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

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

(8)

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

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

(9)

1.2 辛本征問題的求解

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

η(x)=φeμxx

(10)

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

(11)

式中，kb=(ρhω2/D)1/4為板的自由彎曲波波數。求解式(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)為待定系數，記

s={s1s2s3s4}T

(18)

為基本系數向量。

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

1.3 求解波傳播參數與基本系數向量

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

對邊簡支邊界條件為

(19)

代入式(17)可得到

Ks=0

(20)

其中系數矩陣K為

(21)

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

sin(k1a)sin(k2a)=0

(22)

因此波傳播參數為

(23)

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

(24)

對邊固支邊界條件為

(25)

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

k1k2cos(k1a)cos(k2a)=0

(26)

以及基本系數向量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，意味著板沒有變形，也沒有內力，是一組平凡解。雖然是超越方程(22)或(26)的解，但在后面的分析中將不予考慮。另外注意到本文研究的是穩(wěn)態(tài)受迫振動問題，激勵頻率不為零，所以不用考慮靜力分析中波傳播參數的重根問題。

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

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

(28)

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

(29)

2 受迫振動的波傳播分析

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

2.1 確定直接激勵波波幅

外部作用在無窮波導上會向兩側產生直接激勵波，波幅用e+和e-表示。正、負向直接激勵波應使結構在外載作用位置滿足位移連續(xù)性以及力平衡條件。在波空間下，根據式(28)，上述條件可以表達為

(30)

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

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

(31)

2.2 邊界反射系數矩陣

正向波v+入射到邊界上會產生負向波v-。v-與v+之間滿足的關系需要根據邊界條件給出。以簡支邊界條件為例，假設在y=b處簡支。相應邊界條件的變分式為

(32)

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

k=1,2,m

(33)

式(33)寫成矩陣形式為

(34)

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

(35)

由式(34)可以得到

(36)

式中，R為邊界反射系數矩陣。

注意，這里的波幅右上標的正負號分別代表的是入射波和反射波。對于對稱的另一條側邊，v-為入射波，v+為反射波。于是，右側邊界的波反射系數為

(37)

如果左側y=0處也是簡支邊界條件，則對應的邊界波反射系數為

(38)

2.3 建立波傳播關系并求解

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

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

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

(39)

由波的傳播與反射關系可知

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)，可以推導得到

圖2 波導中的波傳播示意圖

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處的波幅為

c+=T(y-ye)a+

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

(44)

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

3 算 例

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

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

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

(45)

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

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

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

圖6 板的應變能

4 結 論

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

[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] 姚偉岸，鐘萬勰.辛彈性力學[M].北京:高等教育出版社，2002.

[9] 鐘萬勰.應用力學對偶體系[M].北京: 科學出版社，2002.

[10] 鐘萬勰.應用力學的辛數學方法[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] 馬國軍,徐新生,郭杏林.旋轉運動中彈性梁耦合振動的辛方法[J].計算力學學報，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] 談梅蘭,吳光,王鑫偉.矩形薄板面內非線性分布載荷下的辛彈性力學解[J].工程力學，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].力學學報，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].大連:大連理工大學，2009.

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

[17] 鮑四元,鄧子辰.哈密頓體系下矩形薄板自由振動的一般解[J].動力學與控制學報，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].大連:大連理工大學，2012.

[19] 鐘陽,李銳,田斌.矩形中厚板自由振動問題的哈密頓體系與辛幾何解法[J].動力學與控制學報，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.