基于辛理论的Bernoulli-Euler梁波散射分析

郑罡, 唐宇, 蔡汶秀, 曾广榕

(重庆交通大学省部共建山区桥梁及隧道工程国家重点实验室, 重庆 400074)

Bernoulli-Euler梁理论仅考虑了梁的转动和截面惯性力,形式简单易用,对于细长梁,采用该模型即可以满足精度要求[1]。 因此,Bernoulli-Euler梁被广泛用于工程结构中。文献[2-3]根据Bernoulli-Euler梁理论以及非局部理论,研究了充流碳纳米管的自由振动和波动问题。 张涛等[4]将隧道简化为Bernoulli-Euler梁模型,提出了一种临近堆载对隧道变形响应的计算方法。徐梅玲等[5]研究了忽略剪切变形对Bernoulli-Euler梁振型的影响。付艳艳等[6]基于回传射线矩阵法给出了黏弹性 Pasternak 地基上6种经典边界下的 Bernoulli-Euler 梁的自振频率方程和模态函数表达式,并分析了不同边界条件对自振频率和模态的影响。文献[7-9]基于Bernoulli-Euler梁理论,研究了声子晶体的振动特性。张靖华等[10]在Hamilton体系下探讨了Bernoulli-Euler梁的动力问题,展现了Hamilton体系在解决弹性力学问题独特的优越性。

应用力学从Lagrange系统转换为Hamilton系统,从而实现了从传统的欧几里得几何空间到辛几何空间的转变[11],避免了传统方法求解高阶偏微分方程的困难[12]。辛体系下的Hamilton矩阵的本征向量间满足共轭辛正交[13],从而本征函数的正交性和完备性得到了保证,极大地拓宽了Sturm-Liouville问题,即分离变量法的适用范围[14]。

能带分析是实际工程中用于分析结构振动特性的常用方法,如声子晶体[7-9]、桁架结构[15]、铁路轨道[16]等,因此研究结构的能带是重要的基础工作。振动与波是紧密结合在一起的[11],因此在研究能带的基础上,进一步可以考虑波的传播问题,进而引出波的散射问题。Zhong等[17-18]采用辛本征展开、辛正交归一等方法,对波传播问题中的能带结构和散射做了详细研究。Zhang等[19]基于辛理论开发了一种用于数值模拟碳纳米管的色散关系的全新计算方法和结构模型,其与传统模型相比更加有效、效率更高,说明了辛数学理论在波的界带分析的研究中具有潜力。钟万勰[11]首先给出了弹性力学中的辛数学方法,并建立了Hamilton系统在辛体系下的一般求解方法,推导了Timoshenko梁的Hamilton对偶方程。吴锋等[20]则在此基础上深入研究了带有端部散射体的Timoshenko梁的能带结构与波散射问题。李渊等[21]又基于此将载流碳纳米管等效为Timoshenko梁,将辛弹性理论用于研究碳纳米管波的传播问题。关于辛体系下的Timoshenko梁波的散射的研究并不少见[20-21],但是对于理论形式更简单的Bernoulli-Euler梁在辛体系下研究其波的散射问题尚鲜见报道。

鉴于此,以Bernoulli-Euler梁理论为基础,基于应用力学的辛数学方法[11]及Lagrange乘子法[22],推导Bernoulli-Euler梁的Hamilton对偶微分方程,研究Bernoulli-Euler梁的能带结构,分析波的散射问题,将禁带和通带解耦为两个部分,降低计算复杂性,为辛体系下研究波在梁中的传播问题提供了一种思路。

1 Hamilton对偶方程

图1 Bernoulli-Euler梁和刚体的连接模型Fig.1 The connection model of Bernoulli-Euler beam and rigid body

为了建立Bernoulli-Euler梁的辛对偶求解体系,需要得到梁的Hamilton对偶方程。首先,构造其Lagrange密度函数为

L=T-V

(1)

式(1)中:梁的动能密度函数T和应变能密度函数V分别表示为

(2)

式(2)中:t和x分别为时间、空间自变量;ρ为梁的密度;A为截面面积;E为弹性模量;I为截面惯性矩。

在分析梁的波动问题时,经常采用将时域转化为频域的方法,令

(3)

式(3)中:ω为圆频率;u(x,ω)和φ(x,ω)分别为频域下的线位移和角位移;i为虚数单位。

(4)

Lagrange乘子法是构造泛函时用于解除约束的常用方法,利用Bernoulli-Euler梁的假设条件,令

u′=φ,φ′=s

(5)

式(5)中:s为角位移φ的斜率。

在Lagrange密度函数L上添加辅助项p1(u′-φ)及p2(φ′-s),其中p1、p2为Lagrange乘子,于是构造广义Lagrange密度函数L*为

L*=L+p1(u′-φ)+p2(φ′-s)

(6)

根据势能变分原理得

(7)

式(7)中:δ表示对函数变分。

于是,可以推导出以下约束条件。

(8)

为了方便分析,引入位移向量q=(u,φ)T和应力向量p=(p1,p2)T,则根据Hamilton变分原理有

(9)

Hamilton函数可表示为

(10)

因此,对式(10)求偏导便可以得到对偶变量为

(11)

把约束条件及式(10)代入式(11),可以得到Hamilton对偶方程为

(12)

式(12)中:AT为A的转置矩阵。

引入状态向量v及矩阵H,可表示为

(13)

则Hamilton对偶方程[式(12)]可以写为

v′=Hv

(14)

易知,矩阵H满足JH=(JH)T,其中J为辛单位矩阵,因此矩阵H为Hamilton矩阵。

(15)

式(15)中:I2为2×2的单位矩阵。

2 动力方程

如图1所示的结构模型,其右端的Bernoulli-Euler梁的内力可分别表示为

(16)

梁的动力方程为

(17)

结合式(16)和式(17),通过消除内力和应变,可以得到以位移表示的Bernoulli-Euler梁动力方程为

(18)

假设Bernoulli-Euler梁的左端在x=0处有

(19)

左端刚体的动力方程同样也可写为

(20)

(21)

对于刚体,其线位移为线性变化,角位移为常数,因此根据刚体在x=0处的线位移和角位移可得

(22)

将式(22)代入式(20)可以得到刚体的动力方程为

(23)

把式(3)代入式(18),于是Bernoulli-Euler梁在频域下的动力方程可以写为

(24)

结合约束条件和式(16)应力向量可表示为

(25)

为了方便分析,把内力转化到频域下

(26)

于是,频域下x=0处的对偶变量有

p(0)=(-Q0,-M0)T,q(0)=(u0,φ0)T

(27)

因此,左端刚体的动力方程可表示为

(28)

结合式(27)和式(28)可以把刚体的动力方程写为

p(0)=Krq(0)

(29)

于是,根据式(29)便可以得到刚体的动力刚度矩阵Kr,可表示为

(30)

3 Bernoulli-Euler梁波的散射

通过本征向量展开法求解Hamilton对偶方程[式(14)],给定一个圆频率ω,即可求得相应的本征值与本征向量。令行列式det(H-μI)=0,展开得

(31)

由式(31)可知,μ2的解为一个正根和一个负根,则矩阵H的本征值μ可以分为两类,在虚轴上的2个本征值对应的是通带(Passband)本征值,这一类本征值表示波可以在梁中传播;在实轴上的2个本征值对应的是禁带(Stopband)本征值,表现为波的局部振动,此时波无法在梁体中传播,自然也不能发生波的散射。因此,在分析Bernoulli-Euler梁波的散射时,首先需要分离禁带。

对于图1所示的结构,假设有频率为ω的波在右端半无穷长Bernoulli-Euler梁中传播,则其在x=0处的状态向量v0可由辛体系下的Hamilton矩阵的本征向量构成,因此可以采用本征向量展开法对其展开求解。

令Bernoulli-Euler梁的Hamilton矩阵的两个通带本征值记为iμp和-iμp,对应的本征向量分别记为φap和φbp,下标“a”和“b”分别用于区分入射波和反射波,且要求这两个本征向量满足辛归一性,即

(32)

若不满足,则需要先进行辛归一化[20],令

(33)

通过式(33)对本征向量进行转化,便能实现本征向量的辛归一。其中,φar和φai分别为本征向量φap的实部和虚部,向量的上标“-”则表示该向量的共轭向量。

由式(33)可知, Bernoulli-Euler梁的状态向量是由1对通带本征向量和1对禁带本征向量所组成的,而因为禁带本征向量不能传递波,只会发生波的局部振动,因此,在分析Bernoulli-Euler梁的散射问题时,状态向量可以只考虑通带的1个入射波和1个反射波。

v=aφap+bφbp

(34)

式(34)中:a和b分别为反射波向量和入射波向量,它们之间存在关系a=Scab。

因此,为了研究波的散射问题,需要解出a和b,得到散射矩阵Sca。

由于是基于本征向量展开法求解Bernoulli-Euler梁的Hamilton对偶方程,所以应该认为在x方向上梁体是均匀的,波的传播是彼此独立的,只有在梁体与刚体不均匀的连接部分才会发生波的散射,因此用状态向量法研究波的散射,要求梁和刚体在连接处的状态向量一致。

Hamilton矩阵的2个通带本征向量组成了一个二维通带子空间,通带子空间的基为通带本征向量φap和φbp,同时为了满足辛归一性,根据式(33)可知,通带子空间的基也可认为是由通带本征向量φap的实部φar和虚部φai所组成的;2个禁带本征向量组成了一个二维禁带子空间。为了简化计算,在计算出通带子空间的基后,可以通过任意构造一组禁带子空间的基底向量,再采用辛Gram-Schmidt方法,就可以得到禁带子空间的基。为了方便分析,将通带子空间的基和禁带子空间的基组合为矩阵形式,可表示为

C(ω2)=(VapVasVbpVbs)

(35)

式(35)中:

(36)

(37)

令

(38)

式(38)中:a1、a2、b1和b2为待定系数,可通过辛正交求得

(39)

于是便可求出待定系数a1、a2、b1和b2,将其代入式(39)即可得到禁带基向量Vas和Vbs。经过辛Gram-Schmidt正交归一化后得到的矩阵C(ω2)是辛矩阵,把通带基Vbp和禁带基Vas重新排列得

Cex(ω2)=(VapVbpVasVbs)

(40)

将Cex(ω2)用于转化Hamilton矩阵H得

(41)

式(41)结果表明,Bernoulli-Euler梁的Hamilton矩阵经过相似变换后得到的矩阵Hex为对角矩阵,且其对角线上的子矩阵为被拆分为两个对立部分的通带子空间矩阵Hp和禁带子空间矩阵Hs。其中,只有通带子空间矩阵Hp能传递波,而禁带子空间矩阵Hs则会把波限制在刚体附近,发生局部振动。

把式(41)代入式(14),则Bernoulli-Euler梁的通带Hamilton对偶方程和禁带Hamilton对偶方程可以分别表示为

v′pex=Hpvpex

(42)

v′sex=Hpvsex

(43)

式中:

(44)

式(44)中:qpb和ppb分别为通带子空间的位移向量和力向量;qsb和psb分别为禁带子空间的位移向量和力向量,这两个部分彼此独立。

禁带Hamilton对偶方程的解表示梁的局部振动,所以对于半无穷长梁来说,在x=0处的力向量qsb,0只与位移向量psb,0有关,而与无穷远处的状态向量无关,qsb,0与psb,0之间有psb,0=-Q∞qsb,0的关系,其中Q∞为半无穷长梁左端禁带动力刚度矩阵,可以通过精细积分法进行计算Q∞,该算法在文献[11]中有详细介绍,此处不再赘述。对于刚体,因为其右端与Bernoulli-Euler梁相连,因此它的动力方程[式(29)]也需要进行相同的变换。

首先,把C(ω2)用分块矩阵表示为

(45)

式(45)中:C(ω2)的子矩阵Qa、Qb、Pa和Pb均为2×2的实数矩阵,分别表示Cex(ω2)为C(ω2)的一个排列变换,结合式(42)和式(44)可得

(46)

式(46)中:ac和bc分别为广义位移向量和广义力向量,结合式(29)和式(46)可得

(47)

式(47)中:Krc为广义动力刚度矩阵,可表示为

(48)

式(48)中:广义动力刚度矩阵Krc的子矩阵Krcpp和Krcsp为通带子空间的广义刚度矩阵;Krcps和Krcss为禁带子空间的广义刚度矩阵。

注意到禁带子空间的广义位移psb和广义力qsb之间有psb=-Q∞qsb的关系,由于禁带不能传递波,因此可以利用禁带动力刚度矩阵Q∞通过凝聚作用将qsb引起的位移psb凝聚掉

(49)

此时,就可以得到不涉及禁带部分的刚体动力方程为

Krpqpb=ppb

(50)

式(50)中:Krp为通带部分刚体动力刚度矩阵。

Krp=Krcpp-Krcps(Krcss+Q∞)-1Krcsp

(51)

至此,便得到了只包含通带部分的Bernoulli-Euler梁的动力方程和左端刚体的动力方程。

为了分析Bernoulli-Euler梁波的散射,则需要联立式(42)和式(50)。假设通带子空间矩阵Hp的两个本征值分别为iμpb和-iμpb,其对应的两个本征向量分别为φapb和φbpb,且φapb和φbpb需要满足辛归一原则φapbJφbpb=1,若不满足可以按照式(33)的做法进行处理。

为了方便处理,将通带子空间的Hamilton矩阵的本征向量φ用矩阵表示为

(52)

把式(52)代入式(34)可以得到用矩阵表示的位移向量和力向量,分别表示为

qpb=Qapa+Qbpb,ppb=Papa+Pbpb

(53)

为了求解待定反射波向量a,给定入射波向量b,把式(53)代入式(50)得

a=-(Pap-KrpQap)-1(Pbp-KrpQbp)b

(54)

由于a和b之间有关系a=Scab,因此可以得到散射矩阵

Sca=-(Pap-KrpQap)-1(Pbp-KrpQbp)

(55)

4 算例分析

算例采用2根半无穷长Bernoulli-Euler工字型截面梁,其中1#梁[20]和2#梁[5]的材料参数和截面积分别如表1所示。

表1 工字型截面梁的参数

梁的左端与一个长度l为0.1 m的刚体连接,其截面惯性矩Ir、密度ρr和截面积Ar均与梁的材料参数和截面积相同。利用不同的本征值μ代入式(31)计算该梁的能带,能带曲线如图2所示。

图2 Bernoulli-Euler梁能带结构Fig.2 The energy band of Bernoulli-Euler beam

表2 散射矩阵

文献[20]分析了Timoshenko梁的能带结构,其剪切模量G和剪切系数K分别取为G=4×104Pa和K=5/6,其余参数均和1#梁相同,能带图如图3所示。

图3 Timoshenko梁能带结构Fig.3 The energy band of Timoshenko beam

通过对比图2(a)和图3可以发现,在低频(0～7 907 rad/s)范围内,Timoshenko梁理论增加的转动惯量项和剪切变形项所带来的影响较小,其能带结构(图3)和Bernoulli-Euler梁的能带结构图类似,即状态向量均是由一对通带本征向量和一对禁带本征向量构成。同时也可以发现随着μ的增加Bernoulli-Euler梁的频率明显高于Timoshenko梁的频率,这说明Bernoulli-Euler梁由于忽略了转动惯量和剪切变形的影响高估了梁的频率[5],因此在实际工程中需要慎重选择梁模型来满足工程精度要求。

观察图3可以发现,当Timoshenko梁的频率超过7 907 rad/s时,其能带结构图增加了一条剪切波的能带曲线(Curve2),此时其状态向量是由两对通带本征向量构成,这是Timoshenko梁理论中多了一个关于时间的四阶导数项所导致的。可以得到结论:Bernoulli-Euler梁和Timoshenko梁的能带结构形式在低频(0 ～7 907 rad/s)范围无明显差异,但是在高频(7 907 rad/s以上)范围,二者有显著不同。

5 结论

基于辛弹性力学理论和Bernoulli-Euler梁理论,研究Bernoulli-Euler梁波的传播问题。通过引入Lagrange乘子,导出对偶变量,推导Bernoulli-Euler梁在Hamilton辛对偶体系下的动力学方程。采用辛数学方法得到Bernoulli-Euler梁的能带结构表达式,计算Bernoulli-Euler梁波的散射矩阵,揭示了状态向量的物理意义,体现出辛体系的优越性。通过数值算例体现出Bernoulli-Euler梁和Timoshenko梁在不同频率范围内能带结构的差异性。本文方法可以为车辆轨道、声子晶体及大跨度空间结构等具有周期性特征结构的能带研究和波散射分析提供理论依据。