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

鄭罡, 唐宇, 蔡汶秀, 曾廣榕

(重慶交通大學(xué)省部共建山區(qū)橋梁及隧道工程國家重點實驗室, 重慶 400074)

Bernoulli-Euler梁理論僅考慮了梁的轉(zhuǎn)動和截面慣性力,形式簡單易用,對于細長梁,采用該模型即可以滿足精度要求[1]。 因此,Bernoulli-Euler梁被廣泛用于工程結(jié)構(gòu)中。文獻[2-3]根據(jù)Bernoulli-Euler梁理論以及非局部理論,研究了充流碳納米管的自由振動和波動問題。 張濤等[4]將隧道簡化為Bernoulli-Euler梁模型,提出了一種臨近堆載對隧道變形響應(yīng)的計算方法。徐梅玲等[5]研究了忽略剪切變形對Bernoulli-Euler梁振型的影響。付艷艷等[6]基于回傳射線矩陣法給出了黏彈性 Pasternak 地基上6種經(jīng)典邊界下的 Bernoulli-Euler 梁的自振頻率方程和模態(tài)函數(shù)表達式,并分析了不同邊界條件對自振頻率和模態(tài)的影響。文獻[7-9]基于Bernoulli-Euler梁理論,研究了聲子晶體的振動特性。張靖華等[10]在Hamilton體系下探討了Bernoulli-Euler梁的動力問題,展現(xiàn)了Hamilton體系在解決彈性力學(xué)問題獨特的優(yōu)越性。

應(yīng)用力學(xué)從Lagrange系統(tǒng)轉(zhuǎn)換為Hamilton系統(tǒng),從而實現(xiàn)了從傳統(tǒng)的歐幾里得幾何空間到辛幾何空間的轉(zhuǎn)變[11],避免了傳統(tǒng)方法求解高階偏微分方程的困難[12]。辛體系下的Hamilton矩陣的本征向量間滿足共軛辛正交[13],從而本征函數(shù)的正交性和完備性得到了保證,極大地拓寬了Sturm-Liouville問題,即分離變量法的適用范圍[14]。

能帶分析是實際工程中用于分析結(jié)構(gòu)振動特性的常用方法,如聲子晶體[7-9]、桁架結(jié)構(gòu)[15]、鐵路軌道[16]等,因此研究結(jié)構(gòu)的能帶是重要的基礎(chǔ)工作。振動與波是緊密結(jié)合在一起的[11],因此在研究能帶的基礎(chǔ)上,進一步可以考慮波的傳播問題,進而引出波的散射問題。Zhong等[17-18]采用辛本征展開、辛正交歸一等方法,對波傳播問題中的能帶結(jié)構(gòu)和散射做了詳細研究。Zhang等[19]基于辛理論開發(fā)了一種用于數(shù)值模擬碳納米管的色散關(guān)系的全新計算方法和結(jié)構(gòu)模型,其與傳統(tǒng)模型相比更加有效、效率更高,說明了辛數(shù)學(xué)理論在波的界帶分析的研究中具有潛力。鐘萬勰[11]首先給出了彈性力學(xué)中的辛數(shù)學(xué)方法,并建立了Hamilton系統(tǒng)在辛體系下的一般求解方法,推導(dǎo)了Timoshenko梁的Hamilton對偶方程。吳鋒等[20]則在此基礎(chǔ)上深入研究了帶有端部散射體的Timoshenko梁的能帶結(jié)構(gòu)與波散射問題。李淵等[21]又基于此將載流碳納米管等效為Timoshenko梁,將辛彈性理論用于研究碳納米管波的傳播問題。關(guān)于辛體系下的Timoshenko梁波的散射的研究并不少見[20-21],但是對于理論形式更簡單的Bernoulli-Euler梁在辛體系下研究其波的散射問題尚鮮見報道。

鑒于此,以Bernoulli-Euler梁理論為基礎(chǔ),基于應(yīng)用力學(xué)的辛數(shù)學(xué)方法[11]及Lagrange乘子法[22],推導(dǎo)Bernoulli-Euler梁的Hamilton對偶微分方程,研究Bernoulli-Euler梁的能帶結(jié)構(gòu),分析波的散射問題,將禁帶和通帶解耦為兩個部分,降低計算復(fù)雜性,為辛體系下研究波在梁中的傳播問題提供了一種思路。

1 Hamilton對偶方程

圖1 Bernoulli-Euler梁和剛體的連接模型Fig.1 The connection model of Bernoulli-Euler beam and rigid body

為了建立Bernoulli-Euler梁的辛對偶求解體系,需要得到梁的Hamilton對偶方程。首先,構(gòu)造其Lagrange密度函數(shù)為

L=T-V

(1)

式(1)中:梁的動能密度函數(shù)T和應(yīng)變能密度函數(shù)V分別表示為

(2)

式(2)中:t和x分別為時間、空間自變量;ρ為梁的密度;A為截面面積;E為彈性模量;I為截面慣性矩。

在分析梁的波動問題時,經(jīng)常采用將時域轉(zhuǎn)化為頻域的方法,令

(3)

式(3)中:ω為圓頻率;u(x,ω)和φ(x,ω)分別為頻域下的線位移和角位移;i為虛數(shù)單位。

(4)

Lagrange乘子法是構(gòu)造泛函時用于解除約束的常用方法,利用Bernoulli-Euler梁的假設(shè)條件,令

u′=φ,φ′=s

(5)

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

在Lagrange密度函數(shù)L上添加輔助項p1(u′-φ)及p2(φ′-s),其中p1、p2為Lagrange乘子,于是構(gòu)造廣義Lagrange密度函數(shù)L*為

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

(6)

根據(jù)勢能變分原理得

(7)

式(7)中:δ表示對函數(shù)變分。

于是,可以推導(dǎo)出以下約束條件。

(8)

為了方便分析,引入位移向量q=(u,φ)T和應(yīng)力向量p=(p1,p2)T,則根據(jù)Hamilton變分原理有

(9)

Hamilton函數(shù)可表示為

(10)

因此,對式(10)求偏導(dǎo)便可以得到對偶變量為

(11)

把約束條件及式(10)代入式(11),可以得到Hamilton對偶方程為

(12)

式(12)中:AT為A的轉(zhuǎn)置矩陣。

引入狀態(tài)向量v及矩陣H,可表示為

(13)

則Hamilton對偶方程[式(12)]可以寫為

v′=Hv

(14)

易知,矩陣H滿足JH=(JH)T,其中J為辛單位矩陣,因此矩陣H為Hamilton矩陣。

(15)

式(15)中:I2為2×2的單位矩陣。

2 動力方程

如圖1所示的結(jié)構(gòu)模型,其右端的Bernoulli-Euler梁的內(nèi)力可分別表示為

(16)

梁的動力方程為

(17)

結(jié)合式(16)和式(17),通過消除內(nèi)力和應(yīng)變,可以得到以位移表示的Bernoulli-Euler梁動力方程為

(18)

假設(shè)Bernoulli-Euler梁的左端在x=0處有

(19)

左端剛體的動力方程同樣也可寫為

(20)

(21)

對于剛體,其線位移為線性變化,角位移為常數(shù),因此根據(jù)剛體在x=0處的線位移和角位移可得

(22)

將式(22)代入式(20)可以得到剛體的動力方程為

(23)

把式(3)代入式(18),于是Bernoulli-Euler梁在頻域下的動力方程可以寫為

(24)

結(jié)合約束條件和式(16)應(yīng)力向量可表示為

(25)

為了方便分析,把內(nèi)力轉(zhuǎn)化到頻域下

(26)

于是,頻域下x=0處的對偶變量有

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

(27)

因此,左端剛體的動力方程可表示為

(28)

結(jié)合式(27)和式(28)可以把剛體的動力方程寫為

p(0)=Krq(0)

(29)

于是,根據(jù)式(29)便可以得到剛體的動力剛度矩陣Kr,可表示為

(30)

3 Bernoulli-Euler梁波的散射

通過本征向量展開法求解Hamilton對偶方程[式(14)],給定一個圓頻率ω,即可求得相應(yīng)的本征值與本征向量。令行列式det(H-μI)=0,展開得

(31)

由式(31)可知,μ2的解為一個正根和一個負根,則矩陣H的本征值μ可以分為兩類,在虛軸上的2個本征值對應(yīng)的是通帶(Passband)本征值,這一類本征值表示波可以在梁中傳播;在實軸上的2個本征值對應(yīng)的是禁帶(Stopband)本征值,表現(xiàn)為波的局部振動,此時波無法在梁體中傳播,自然也不能發(fā)生波的散射。因此,在分析Bernoulli-Euler梁波的散射時,首先需要分離禁帶。

對于圖1所示的結(jié)構(gòu),假設(shè)有頻率為ω的波在右端半無窮長Bernoulli-Euler梁中傳播,則其在x=0處的狀態(tài)向量v0可由辛體系下的Hamilton矩陣的本征向量構(gòu)成,因此可以采用本征向量展開法對其展開求解。

令Bernoulli-Euler梁的Hamilton矩陣的兩個通帶本征值記為iμp和-iμp,對應(yīng)的本征向量分別記為φap和φbp,下標(biāo)“a”和“b”分別用于區(qū)分入射波和反射波,且要求這兩個本征向量滿足辛歸一性,即

(32)

若不滿足,則需要先進行辛歸一化[20],令

(33)

通過式(33)對本征向量進行轉(zhuǎn)化,便能實現(xiàn)本征向量的辛歸一。其中,φar和φai分別為本征向量φap的實部和虛部,向量的上標(biāo)“-”則表示該向量的共軛向量。

由式(33)可知, Bernoulli-Euler梁的狀態(tài)向量是由1對通帶本征向量和1對禁帶本征向量所組成的,而因為禁帶本征向量不能傳遞波,只會發(fā)生波的局部振動,因此,在分析Bernoulli-Euler梁的散射問題時,狀態(tài)向量可以只考慮通帶的1個入射波和1個反射波。

v=aφap+bφbp

(34)

式(34)中:a和b分別為反射波向量和入射波向量,它們之間存在關(guān)系a=Scab。

因此,為了研究波的散射問題,需要解出a和b,得到散射矩陣Sca。

由于是基于本征向量展開法求解Bernoulli-Euler梁的Hamilton對偶方程,所以應(yīng)該認為在x方向上梁體是均勻的,波的傳播是彼此獨立的,只有在梁體與剛體不均勻的連接部分才會發(fā)生波的散射,因此用狀態(tài)向量法研究波的散射,要求梁和剛體在連接處的狀態(tài)向量一致。

Hamilton矩陣的2個通帶本征向量組成了一個二維通帶子空間,通帶子空間的基為通帶本征向量φap和φbp,同時為了滿足辛歸一性,根據(jù)式(33)可知,通帶子空間的基也可認為是由通帶本征向量φap的實部φar和虛部φai所組成的;2個禁帶本征向量組成了一個二維禁帶子空間。為了簡化計算,在計算出通帶子空間的基后,可以通過任意構(gòu)造一組禁帶子空間的基底向量,再采用辛Gram-Schmidt方法,就可以得到禁帶子空間的基。為了方便分析,將通帶子空間的基和禁帶子空間的基組合為矩陣形式,可表示為

C(ω2)=(VapVasVbpVbs)

(35)

式(35)中:

(36)

(37)

令

(38)

式(38)中:a1、a2、b1和b2為待定系數(shù),可通過辛正交求得

(39)

于是便可求出待定系數(shù)a1、a2、b1和b2,將其代入式(39)即可得到禁帶基向量Vas和Vbs。經(jīng)過辛Gram-Schmidt正交歸一化后得到的矩陣C(ω2)是辛矩陣,把通帶基Vbp和禁帶基Vas重新排列得

Cex(ω2)=(VapVbpVasVbs)

(40)

將Cex(ω2)用于轉(zhuǎn)化Hamilton矩陣H得

(41)

式(41)結(jié)果表明,Bernoulli-Euler梁的Hamilton矩陣經(jīng)過相似變換后得到的矩陣Hex為對角矩陣,且其對角線上的子矩陣為被拆分為兩個對立部分的通帶子空間矩陣Hp和禁帶子空間矩陣Hs。其中,只有通帶子空間矩陣Hp能傳遞波,而禁帶子空間矩陣Hs則會把波限制在剛體附近,發(fā)生局部振動。

把式(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有關(guān),而與無窮遠處的狀態(tài)向量無關(guān),qsb,0與psb,0之間有psb,0=-Q∞qsb,0的關(guān)系,其中Q∞為半無窮長梁左端禁帶動力剛度矩陣,可以通過精細積分法進行計算Q∞,該算法在文獻[11]中有詳細介紹,此處不再贅述。對于剛體,因為其右端與Bernoulli-Euler梁相連,因此它的動力方程[式(29)]也需要進行相同的變換。

首先,把C(ω2)用分塊矩陣表示為

(45)

式(45)中:C(ω2)的子矩陣Qa、Qb、Pa和Pb均為2×2的實數(shù)矩陣,分別表示Cex(ω2)為C(ω2)的一個排列變換,結(jié)合式(42)和式(44)可得

(46)

式(46)中:ac和bc分別為廣義位移向量和廣義力向量,結(jié)合式(29)和式(46)可得

(47)

式(47)中:Krc為廣義動力剛度矩陣,可表示為

(48)

式(48)中:廣義動力剛度矩陣Krc的子矩陣Krcpp和Krcsp為通帶子空間的廣義剛度矩陣;Krcps和Krcss為禁帶子空間的廣義剛度矩陣。

注意到禁帶子空間的廣義位移psb和廣義力qsb之間有psb=-Q∞qsb的關(guān)系,由于禁帶不能傳遞波,因此可以利用禁帶動力剛度矩陣Q∞通過凝聚作用將qsb引起的位移psb凝聚掉

(49)

此時,就可以得到不涉及禁帶部分的剛體動力方程為

Krpqpb=ppb

(50)

式(50)中:Krp為通帶部分剛體動力剛度矩陣。

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

(51)

至此,便得到了只包含通帶部分的Bernoulli-Euler梁的動力方程和左端剛體的動力方程。

為了分析Bernoulli-Euler梁波的散射,則需要聯(lián)立式(42)和式(50)。假設(shè)通帶子空間矩陣Hp的兩個本征值分別為iμpb和-iμpb,其對應(yīng)的兩個本征向量分別為φ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之間有關(guān)系a=Scab,因此可以得到散射矩陣

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

(55)

4 算例分析

算例采用2根半無窮長Bernoulli-Euler工字型截面梁,其中1#梁[20]和2#梁[5]的材料參數(shù)和截面積分別如表1所示。

表1 工字型截面梁的參數(shù)

梁的左端與一個長度l為0.1 m的剛體連接,其截面慣性矩Ir、密度ρr和截面積Ar均與梁的材料參數(shù)和截面積相同。利用不同的本征值μ代入式(31)計算該梁的能帶,能帶曲線如圖2所示。

圖2 Bernoulli-Euler梁能帶結(jié)構(gòu)Fig.2 The energy band of Bernoulli-Euler beam

表2 散射矩陣

文獻[20]分析了Timoshenko梁的能帶結(jié)構(gòu),其剪切模量G和剪切系數(shù)K分別取為G=4×104Pa和K=5/6,其余參數(shù)均和1#梁相同,能帶圖如圖3所示。

圖3 Timoshenko梁能帶結(jié)構(gòu)Fig.3 The energy band of Timoshenko beam

通過對比圖2(a)和圖3可以發(fā)現(xiàn),在低頻(0～7 907 rad/s)范圍內(nèi),Timoshenko梁理論增加的轉(zhuǎn)動慣量項和剪切變形項所帶來的影響較小,其能帶結(jié)構(gòu)(圖3)和Bernoulli-Euler梁的能帶結(jié)構(gòu)圖類似,即狀態(tài)向量均是由一對通帶本征向量和一對禁帶本征向量構(gòu)成。同時也可以發(fā)現(xiàn)隨著μ的增加Bernoulli-Euler梁的頻率明顯高于Timoshenko梁的頻率,這說明Bernoulli-Euler梁由于忽略了轉(zhuǎn)動慣量和剪切變形的影響高估了梁的頻率[5],因此在實際工程中需要慎重選擇梁模型來滿足工程精度要求。

觀察圖3可以發(fā)現(xiàn),當(dāng)Timoshenko梁的頻率超過7 907 rad/s時,其能帶結(jié)構(gòu)圖增加了一條剪切波的能帶曲線(Curve2),此時其狀態(tài)向量是由兩對通帶本征向量構(gòu)成,這是Timoshenko梁理論中多了一個關(guān)于時間的四階導(dǎo)數(shù)項所導(dǎo)致的。可以得到結(jié)論:Bernoulli-Euler梁和Timoshenko梁的能帶結(jié)構(gòu)形式在低頻(0 ～7 907 rad/s)范圍無明顯差異,但是在高頻(7 907 rad/s以上)范圍,二者有顯著不同。

5 結(jié)論

基于辛彈性力學(xué)理論和Bernoulli-Euler梁理論,研究Bernoulli-Euler梁波的傳播問題。通過引入Lagrange乘子,導(dǎo)出對偶變量,推導(dǎo)Bernoulli-Euler梁在Hamilton辛對偶體系下的動力學(xué)方程。采用辛數(shù)學(xué)方法得到Bernoulli-Euler梁的能帶結(jié)構(gòu)表達式,計算Bernoulli-Euler梁波的散射矩陣,揭示了狀態(tài)向量的物理意義,體現(xiàn)出辛體系的優(yōu)越性。通過數(shù)值算例體現(xiàn)出Bernoulli-Euler梁和Timoshenko梁在不同頻率范圍內(nèi)能帶結(jié)構(gòu)的差異性。本文方法可以為車輛軌道、聲子晶體及大跨度空間結(jié)構(gòu)等具有周期性特征結(jié)構(gòu)的能帶研究和波散射分析提供理論依據(jù)。