基于辛理論的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梁在不同頻率范圍內能帶結構的差異性。本文方法可以為車輛軌道、聲子晶體及大跨度空間結構等具有周期性特征結構的能帶研究和波散射分析提供理論依據。