张亚辉，马永彬

(大连理工大学 工业装备结构分析国家重点实验室，大连 116023)

基于模态思想的有限元方法在进行结构振动分析时，在结构一个振动波长内需要划分6至15甚至更多单元才能准确地模拟结构的振动,而在高频振动下，结构振动波长非常小，应用有限元方法不得不采用大量的自由度来分析结构的振动,因此，高频振动问题需要寻求更为有效的分析方法。统计能量分析(SEA)作为高频振动分析的典型方法[1]，自上世纪60年代提出以来，已经推广到多个领域并得到成功的应用。采用SEA方法进行高频振动分析的计算成本极小，不过只适用于初步验证阶段。因为SEA按振动模式将结构分为若干子系统，分析结果只能给出各子系统能量均方值，随着对结果的需求更加精细化，还需要借助别的方法进行辅助分析。从另一个角度来看，结构的振动可以用波的传播、反射以及传递的形式来表述[2]。这种表述尤其在中、高频域内有优势，因为它只需要很小的计算成本并且有着很高的精确性,因此，近年来也获得结构的弹性波属性近成为了中、高频振动研究领域的热门。针对波导结构的全频域振动响应分析问题，Mace等[3-4]提出了一个新方法：波有限元(WFE)方法。WFE方法结合了有限元方法和周期结构理论，只需选取一小段结构建立有限元模型，便能获得波导结构的波传播属性，进而基于波传播及波散射关系得到问题的解答。相比传统有限元方法，WFE方法在计算成本方面有着很大的优势。目前，这种方法已经成功应用于波色散分析[5]，压电材料分析[6]等领域。然而，由于应用WFE方法时，需要建立一小段结构有限元模型，因而模型的单元属性、网格大小等参数都直接影响到分析结果的精确性，甚至造成严重的误差[7]；于是，一个具体问题的WFE分析需要多次修正有限元模型才能得到满意的结果。因此，如何精确地获得波传播参数是波传播方法首要问题。

钟万勰等[8-10]将Hamilton体系及辛状态空间理论应用到弹性力学，改变了以往以半逆法为主的凑合解法来解决弹性力学问题的现状，从而使得分离变量、辛本征展开等方法能够应用到弹性力学问题。目前已成功应用于梁[11-12]、板[13-14]、壳[15-16]等基本结构原件的弹性力学问题及薄板和中厚板的自由振动分析中[17-19]。将板振动的控制方程导入辛对偶体系，通过求解正则方程，得到的本征值与本征向量恰好分别是波传播方法需要的波传播参数与波形。它们均由完全理性的推导得出，没有引入任何假设。对复杂边界条件，相较于传统方法只能在四边简支的边界条件下才能给出Navier形式的闭合解，辛方法在任意边界条件下均能给出辛解析解，具有显著的优势。

本文基于辛对偶体系和波传播理论，提出了一个分析薄板结构稳态振动响应的新思路。首先，将板振动的控制方程导入辛对偶体系，求解得到本征值(波传播参数)与本征向量(波形)；然后，通过波传播过程中的入射、反射以及传递关系得到各参与波的幅值，进而求得薄板任意位置的响应。与现有文献中基于辛对偶体系求解板的自由振动问题，以及基于传统模态思想求解振动问题不同的是：本文引入波空间的概念，结合波传播理论来研究结构的稳态振动响应。同时，本文以矩阵形式清晰地给出计算公式，在应用时具有更强的通用性。与WFE方法相比，本文方法基于辛本征解，具有更好的精度与数值稳定性；另外，考虑到通用性，相比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(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分别为板横截面单位长度上的剪力、弯矩和扭矩，正方向的规定如图1所示，并且有如下关系式

(3)

图1 薄板及坐标系示意图

定义等效剪力

(4)

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

(5)

式(5)可写为如下矩阵形式

(6)

记为

(7)

式中z={wθFyMy}T为状态向量，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)的解，但在后面的分析中将不予考虑。另外注意到本文研究的是稳态受迫振动问题，激励频率不为零，所以不用考虑静力分析中波传播参数的重根问题。

通过求解超越方程(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]：首先，确定外激励在作用处向两侧无穷波导产生的波的波幅，即直接激励波的波幅；其次，计算波导方向上的边界或者不连续处的波反射系数矩阵；最后，综合直接激励、波反射以及波传播之间的关系，求解可得任意位置正、负向波的波幅c+和c-，进而通过式(28)叠加得到任意位置的位移和内力响应。

2.1 确定直接激励波波幅

外部作用在无穷波导上会向两侧产生直接激励波，波幅用e+和e-表示。正、负向直接激励波应使结构在外载作用位置满足位移连续性以及力平衡条件。在波空间下，根据式(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-，结合各波形的波传播参数，则结构任意位置的位移和内力都可以得到。而波幅的求解可以根据上小节给出的波在边界的反射系数以及本节给出的波传播关系得到。下面通过图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中给出。

采用模态叠加法、波有限元法(WFE)以及本文方法分别求解了激励点输入导纳的幅值以及板的动能、应变能时间均值。

对于四边简支矩形板，可以通过模态叠加得到薄板稳态振动响应的解析解

(45)

事实上，四边简支矩形薄板的辛解析解经过简单的变换，在形式上与式(45)完全一致。然而，由于求解空间不同，本文方法只用较少截断波形便可得到非常精确的结果。在本算例下，本文方法的计算时间4.80 s比WFE方法计算时间242.94 s也少很多。

图3 模态叠加法中不同模态数下输入导纳幅值的相对误差(|Y0|对应1000阶模态)

图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方法除了计算花费时间，还要进行有限元建模及模型的重复修正工作。因此，就分析成本来说，本文方法更有优势。

图6 板的应变能

4 结 论

本文将薄板稳态强迫振动问题引入到辛对偶体系中，结合波传播理论给出了一个求解新思路。和以往文献基于辛对偶体系求解弹性力学以及自由振动问题相比，本文更加强调了波传播的思想，使得辛方法在处理波导结构振动问题时有更为广阔的视野，比如可分析多个波导结构的耦合振动响应。另外，本文以矩阵形式给出了所有计算公式。一方面，使得求解过程更为清楚；另一方面，也使得方法的应用具有更好的通用性。常规有限元方法随着模型的增大计算成本会越来越高，对于高频振动问题有时甚至不能求解。与此相比，本文方法并不受模型大小的影响。和波有限元方法(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] 姚伟岸，孙贞.环扇形薄板弯曲问题的环向辛对偶求解方法[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] 马源.哈密顿体系下弹性圆柱壳的动态屈曲研究[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.