張亞輝，馬永彬

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

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

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

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

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

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

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

(1)

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

對(duì)于穩(wěn)態(tài)振動(dòng)問題，板的撓度可以寫為W(x,y,t)=Re{w(x,y)eiωt},w(x,y)為W(x,y,t)的幅值，因此振動(dòng)方程可通過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為哈密頓算子矩陣，(·)表示對(duì)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ù)向量。

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

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

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

對(duì)邊簡支邊界條件為

(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)

對(duì)邊固支邊界條件為

(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)

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

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

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

(28)

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

(29)

2 受迫振動(dòng)的波傳播分析

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

2.1 確定直接激勵(lì)波波幅

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

(30)

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

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

(31)

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

正向波v+入射到邊界上會(huì)產(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ù)號(hào)分別代表的是入射波和反射波。對(duì)于對(duì)稱的另一條側(cè)邊，v-為入射波，v+為反射波。于是，右側(cè)邊界的波反射系數(shù)為

(37)

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

(38)

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

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

波導(dǎo)結(jié)構(gòu)在ye位置受外部作用，此處的波幅a+和g-可以用直接激勵(lì)波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)點(diǎn)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)以及本文方法分別求解了激勵(lì)點(diǎn)輸入導(dǎo)納的幅值以及板的動(dòng)能、應(yīng)變能時(shí)間均值。

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

(45)

事實(shí)上，四邊簡支矩形薄板的辛解析解經(jīng)過簡單的變換，在形式上與式(45)完全一致。然而，由于求解空間不同，本文方法只用較少截?cái)嗖ㄐ伪憧傻玫椒浅＞_的結(jié)果。在本算例下，本文方法的計(jì)算時(shí)間4.80 s比WFE方法計(jì)算時(shí)間242.94 s也少很多。

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

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

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

4 結(jié) 論

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

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

[17] 鮑四元,鄧子辰.哈密頓體系下矩形薄板自由振動(dòng)的一般解[J].動(dòng)力學(xué)與控制學(xué)報(bào)，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] 鐘陽,李銳,田斌.矩形中厚板自由振動(dòng)問題的哈密頓體系與辛幾何解法[J].動(dòng)力學(xué)與控制學(xué)報(bào)，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.