基于人工彈簧模型的周期結(jié)構(gòu)帶隙計(jì)算方法研究1)

馮青松 楊 舟 郭文杰,2) 陸建飛 梁玉雄

*(華東交通大學(xué)鐵路環(huán)境振動與噪聲教育部工程研究中心,南昌 330013)

?(江蘇大學(xué)土木工程與力學(xué)學(xué)院,江蘇鎮(zhèn)江 212013)

引言

周期結(jié)構(gòu)因其具有良好的彈性波衰減特性[1-2](帶隙特性),在結(jié)構(gòu)振動與噪聲控制領(lǐng)域體現(xiàn)出了巨大的潛力.近年來,周期性結(jié)構(gòu)中彈性波的波動特性研究發(fā)展迅速,從單一的周期性梁[3]、板結(jié)構(gòu)[4-5]拓展到復(fù)雜的三維結(jié)構(gòu)[6-7].與此同時(shí),作為周期結(jié)構(gòu)最為重要的特性,帶隙頻率的放大和調(diào)諧技術(shù)已成為當(dāng)前研究的一個(gè)熱點(diǎn),大量的學(xué)者立足于此,以梁板結(jié)構(gòu)為基礎(chǔ),將智能材料(如壓電材料[8-9]、磁流變材料[10]、功能梯度材料[11]等) 和彈簧-振子系統(tǒng)[12-15]引入周期結(jié)構(gòu)中,以實(shí)現(xiàn)帶隙的主動調(diào)控和帶寬的增大,并取得了很好的效果.

在過去的幾十年里,學(xué)者們提出了大量的周期結(jié)構(gòu)帶隙計(jì)算的數(shù)值和解析方法.有限元法[16-19]作為一種純數(shù)值方法,其在計(jì)算復(fù)雜周期結(jié)構(gòu)帶隙方面表現(xiàn)出強(qiáng)大的能力.但是,有限元方法的精度依賴于離散單元的網(wǎng)格尺寸,高精度往往意味著很大的計(jì)算成本.且有限元法中構(gòu)件耦合和端部邊界條件的處理過程比較繁瑣.后續(xù)也有部分學(xué)者針對復(fù)雜周期結(jié)構(gòu)的帶隙計(jì)算,為提高傳統(tǒng)有限元法的計(jì)算效率,提出了一些新的數(shù)值方法,如波動有限元法[20-21],并將其迅速推廣到超材料桿[22]、多孔板[23]、周期性加筋殼[24]等結(jié)構(gòu)的帶隙計(jì)算之中.但是無論從檢驗(yàn)數(shù)值法的角度,還是從揭示復(fù)雜系統(tǒng)振動機(jī)理的角度,解析或半解析方法不可或缺.目前帶隙計(jì)算的常用解析方法主要有傳遞矩陣法[25-26]、平面波展開法[27-28]等.傳遞矩陣法對研究一維周期結(jié)構(gòu)的帶隙特性非常方便高效,例如帶有周期彈簧支撐的軌道結(jié)構(gòu)[29-30].對于一維結(jié)構(gòu),傳遞矩陣法是一種通用的方法,但用傳遞矩陣法難以處理多維以上的結(jié)構(gòu).平面波展開法是計(jì)算一維、二維和三維周期結(jié)構(gòu)帶隙的傳統(tǒng)方法,但是由于該方法基于對微分方程組的直接求解,因此,在處理一些具有復(fù)雜邊界條件的周期性組合結(jié)構(gòu)時(shí)存在一定的困難.此外,也有學(xué)者提出了一些新的周期結(jié)構(gòu)帶隙計(jì)算方法,諸如譜動剛度法[31]、小波法[32]、漸進(jìn)匹配展開法[33]等,這些方法雖各有特點(diǎn),但處理一些復(fù)雜周期結(jié)構(gòu)模型(如周期性組合結(jié)構(gòu))并不適用.因此,有必要發(fā)展分析復(fù)雜周期結(jié)構(gòu)振動帶隙特性的解析或半解析方法.

能量法具有將求解微分方程邊值問題轉(zhuǎn)化為泛函極值問題的優(yōu)點(diǎn),這有利于復(fù)雜結(jié)構(gòu)系統(tǒng)耦合問題的求解,因而傳統(tǒng)能量法(如Rayleigh-Ritz 法[34-35])被廣泛應(yīng)用于分析一些復(fù)雜組合結(jié)構(gòu)的自由振動問題[36-37],近年來也被引入周期結(jié)構(gòu)帶隙計(jì)算中[38-39].但在利用傳統(tǒng)能量法計(jì)算周期結(jié)構(gòu)帶隙時(shí),存在如下問題:(1) 能量法需要基于Bloch 定理構(gòu)造滿足周期性邊界的位移場形函數(shù)[40-41],從數(shù)學(xué)角度而言,形函數(shù)的周期性重構(gòu)難度較大,而且不同形函數(shù)的構(gòu)造方式不一定相同;(2)重構(gòu)以后的位移場形函數(shù)包含波數(shù),這會導(dǎo)致涉及形函數(shù)的結(jié)構(gòu)的質(zhì)量和剛度矩陣中含有波數(shù),在計(jì)算帶隙時(shí),結(jié)構(gòu)的質(zhì)量和剛度矩陣需隨著波數(shù)的變化進(jìn)行反復(fù)計(jì)算,隨著結(jié)構(gòu)質(zhì)量和剛度矩陣維度的增大,或者掃描波數(shù)點(diǎn)數(shù)的增多,計(jì)算成本也會隨之增大.因此,對于能量法而言,以較低的計(jì)算成本對復(fù)雜周期結(jié)構(gòu)的帶隙進(jìn)行精確計(jì)算是一項(xiàng)具有挑戰(zhàn)性的工作.

鑒于此,本文對傳統(tǒng)能量法進(jìn)行了改進(jìn),提出了一種利用人工彈簧[42-44]來模擬周期邊界條件的方法,將形函數(shù)與周期邊界條件分離,無需構(gòu)造滿足周期邊界的形函數(shù),這樣可解決傳統(tǒng)能量法中構(gòu)造滿足周期邊界條件的形函數(shù)的困難問題;同時(shí),利用人工彈簧模擬周期邊界條件,將周期邊界約束轉(zhuǎn)化為人工彈簧的彈性勢能,這樣只有包含波數(shù)的周期邊界彈性勢能矩陣需參與隨波數(shù)掃描的循環(huán)計(jì)算,其他矩陣只需計(jì)算一次,因而能提高傳統(tǒng)能量法計(jì)算帶隙的效率.此外,人工彈簧模型使用靈活、便捷,可進(jìn)一步拓展到更為復(fù)雜的周期性組合結(jié)構(gòu)的帶隙分析中.

1 人工彈簧模擬周期邊界的基本原理

1.1 問題陳述

以圖1 所示周期性薄板為例,對引言中提到的傳統(tǒng)能量法計(jì)算周期結(jié)構(gòu)帶隙時(shí)存在的問題進(jìn)行具體闡述.

圖1 周期性薄板自由振動分析模型Fig.1 Vibration analysis model of a periodically thin plate with free boundary

對于圖1 所示的周期薄板結(jié)構(gòu),取出一個(gè)周期單元進(jìn)行分析.根據(jù)Kirchhoff-Love 薄板理論,平面內(nèi)位移u和v可以通過一階泰勒級數(shù)圍繞中面層展開,表示為平面外位移w,位移向量可以表示為

將式(1)中的橫向位移w表示為基函數(shù)ξi(x,y)和一個(gè)未知的權(quán)重系數(shù)di(t)的組合,即

式中,符號?表示克羅內(nèi)克積.

根據(jù)Bloch 定理,周期單元在x,y方向上的位移、轉(zhuǎn)角、彎矩和剪力需要滿足Bloch 周期性邊界條件,即

式中,kx,ky表示沿著x,y方向的波數(shù).方程(3)～(6)要求式(2)中的形函數(shù)ξi(x,y)滿足Bloch-Floquet 周期性條件.從數(shù)學(xué)角度而言,形函數(shù)的重構(gòu)過程具有相當(dāng)大的難度,且不同的形函數(shù)構(gòu)造方式不一定相同,這給利用傳統(tǒng)能量法求解周期結(jié)構(gòu)振動帶隙問題帶來了困難.

其中,ρ,h,ν 分別表示板的材料密度、板的厚度以及板材料泊松比;D=E*h3/[12(1-ν2)]表示板的抗彎剛度.由式(7)、式(8)可以得到周期單元的整體拉格朗日量

進(jìn)一步地,定義未知的與時(shí)間相關(guān)向量d(t)=結(jié)合Eular-Lagrange 方程

就能導(dǎo)出周期薄板結(jié)構(gòu)的運(yùn)動方程

式(10) 是一個(gè)標(biāo)準(zhǔn)的特征值方程問題,在第一布里淵區(qū)kx×ky=[-π/a,π/a]×[-π/a,π/a]掃描波數(shù),即能得到周期薄板結(jié)構(gòu)彎曲振動的頻散曲線.值得注意的是,由于重構(gòu)以后的形函數(shù)必然包含波數(shù),這就導(dǎo)致利用位移形函數(shù)表征的結(jié)構(gòu)動能表達(dá)式(7)和應(yīng)變能表達(dá)式(8)中也含有波數(shù),進(jìn)而使得薄板結(jié)構(gòu)運(yùn)動方程(10)中的質(zhì)量矩陣剛度矩陣同樣包含波數(shù).因此,在掃描波數(shù)進(jìn)行帶隙計(jì)算時(shí),結(jié)構(gòu)的質(zhì)量和剛度矩陣需要隨著波數(shù)的變化進(jìn)行循環(huán)計(jì)算,計(jì)算時(shí)間會隨著結(jié)構(gòu)的質(zhì)量和剛度矩陣維度以及掃描波數(shù)點(diǎn)數(shù)的增大而呈倍數(shù)增大,計(jì)算成本也會急劇增大.因此,如上所述,對于傳統(tǒng)能量法而言,以較低的計(jì)算成本對復(fù)雜周期結(jié)構(gòu)的帶隙進(jìn)行精確計(jì)算是一項(xiàng)具有挑戰(zhàn)性的工作.

1.2 人工彈簧模擬周期性邊界條件

嘗試將特殊的周期性邊界條件與位移形函數(shù)完全分離,即無需對形函數(shù)進(jìn)行重新構(gòu)造,用人工彈簧的形式來模擬周期性邊界條件.如圖2 所示,在周期單元的4 條邊界處定義線彈簧和轉(zhuǎn)角彈簧兩類彈簧,它的物理意義是在x和y兩個(gè)方向通過人工彈簧將上一周期單元的尾端與下一周期單元的首端相連.文獻(xiàn)[3]證明了,只需要位移和轉(zhuǎn)角滿足周期性條件,解即收斂于真實(shí)解.

圖2 人工彈簧定義周期性邊界的分析示意圖Fig.2 Analysis schematic diagram of periodic boundary defined by artificial spring

以x方向?yàn)槔?此時(shí)邊界勢能可以表示為

式中w(-a,y)-w(a,y)e-ikx2a,?w(-a,y)/(?x)-[?w(a,y)/(?x)]e-ik2a分別表示第一個(gè)周期首端與第二個(gè)周期首端的位移、轉(zhuǎn)角差值.同時(shí),由式(3)、式(4)可知位移和轉(zhuǎn)角的差值趨于0,即

結(jié)構(gòu)的剪力和彎矩都是一個(gè)有限值,因此式(12) 若要成立,則需要滿足limktx=∞,limktx=∞,同理考慮y方向的周期性邊界,式(11)改寫為

在式(13) 中,周期性邊界條件表示為包含波數(shù)在內(nèi)的彈性勢能的形式,與位移形函數(shù)分離.此時(shí),形函數(shù)ξi(x,y) 無需進(jìn)行構(gòu)造以滿足周期性邊界,周期單元的動能和應(yīng)變能表示為

進(jìn)一步地,周期單元的整體拉格朗日量表示為

結(jié)合前述定義的與時(shí)間相關(guān)未知向量d(t)=以及Eular-Lagrange 方程

式(10)中周期薄板結(jié)構(gòu)的運(yùn)動方程改寫為

通過比較方程(14) 和(15) 與方程(7) 和(8) 可以發(fā)現(xiàn),利用人工彈簧模擬周期性邊界,將位移形函數(shù)與周期性邊界條件分離,選取的形函數(shù)無需滿足周期性邊界,進(jìn)而用形函數(shù)表征的結(jié)構(gòu)動能和應(yīng)變能表達(dá)式中沒有包含波數(shù),同時(shí)通過對比最后的特征值方程(17) 和(10) 可以看出,僅有周期邊界勢能剛度矩陣中包含波數(shù),質(zhì)量與應(yīng)變能剛度矩陣中都不包含波數(shù),在掃描波數(shù)進(jìn)行求解時(shí),除了周期邊界彈性勢能矩陣參與波數(shù)掃描循環(huán)計(jì)算,其他的矩陣僅需計(jì)算一次,這將大幅度提高計(jì)算效率,進(jìn)而降低計(jì)算成本.后續(xù)的算例驗(yàn)證中,將進(jìn)一步證明本文方法的準(zhǔn)確性和效率性.

2 算例驗(yàn)證:周期開口板彎曲振動帶隙計(jì)算

開口板廣泛應(yīng)用于土木、船舶等結(jié)構(gòu)中.開孔板的振動在能量傳遞中起著重要作用,但也與輻射噪聲對周圍環(huán)境的影響有直接關(guān)系.因此,本節(jié)以周期開口板為例,驗(yàn)證人工彈簧模擬周期性邊界的準(zhǔn)確性,同時(shí)說明該方法與傳統(tǒng)能量法相比,在計(jì)算效率上的優(yōu)勢性.

2.1 位移形函數(shù)的選取及特征值方程求解

由于開口板出現(xiàn)了局部厚度突變,若選用改進(jìn)傅里葉級數(shù)、切比雪夫級數(shù)等(適用于描述連續(xù)性較好的結(jié)構(gòu)) 作為位移形函數(shù)可能需要多項(xiàng)級數(shù)來進(jìn)行擬合以達(dá)到精度要求,進(jìn)而降低計(jì)算效率.因此,本文選取具有局域化特性的高斯小波函數(shù)[45-46]為位移形函數(shù),以確保能夠準(zhǔn)確捕捉到厚度變化區(qū)域的波動特征,即對于式(2)定義

式中,p,q和k,r分別表示伸縮因子和平移因子.伸縮因子p,q控制著解的精度,采用越大的尺度因子使得形函數(shù)具有更高的分辨率,進(jìn)而能更精確地模擬更高頻的波場,但同時(shí)要求更多的計(jì)算資源.平移因子k,r由伸縮因子p,q決定,其數(shù)量控制著式(18)中i的取值范圍,也就是形函數(shù)的個(gè)數(shù)和矩陣的維度.伸縮因子與平移因子的取值滿足

式中,ceil(*)表示最接近并大于*的整數(shù),floor(*)表示最接近并小于*的整數(shù).關(guān)于高斯小波函數(shù)的詳細(xì)介紹可以參考文獻(xiàn)[45-46].定義一個(gè)厚度函數(shù)h(x,y)滿足

其中,x0和y0分別表示圖3 所示開口板坐標(biāo)系內(nèi)任意一點(diǎn)的橫坐標(biāo)和縱坐標(biāo);hu表示板的標(biāo)準(zhǔn)厚度;r表示圓形開口半徑.利用式(21) 中的h(x,y) 替換式(14)、式(15) 中的h,掃描波數(shù)求解特征值方程(17)即能得到周期開口板彎曲振動的頻散曲線.

圖3 周期開口板分析模型Fig.3 Analysis model of periodic plate with opening

2.2 準(zhǔn)確性分析

2.2.1 人工彈簧剛度取值收斂性分析

由前述分析可知,當(dāng)用人工彈簧模擬周期邊界條件時(shí),人工彈簧剛度需要取為無窮,但在實(shí)際計(jì)算過程中需要用一個(gè)較大的值來代替,因此位移約束彈簧和轉(zhuǎn)角約束彈簧剛度的取值在很大程度上影響了計(jì)算結(jié)果的準(zhǔn)確性,需要對人工彈簧剛度取值進(jìn)行收斂性分析,找到解收斂時(shí)的彈簧剛度取值,從而能夠準(zhǔn)確替代人工彈簧剛度取無窮這一情況.板的長度和寬度都為0.12 m,板的均質(zhì)厚度取0.001 m,圓形開口半徑為0.036 m,選取鋼材料為板的材料,密度為7800 kg/m3,彈性模量為210 GPa,泊松比為0.3.從頻散曲線中隨機(jī)選取3 條曲線波數(shù)位于kx=ky=π/a處的頻率為研究對象,其結(jié)果如圖4 所示.

圖4 人工彈簧剛度取值收斂性分析Fig.4 Convergence analysis of stiffness of artificial spring

由圖4 可以看出,當(dāng)彈簧剛度系數(shù)數(shù)值大于1010時(shí),得到的波模態(tài)頻率值已經(jīng)趨于穩(wěn)定,基本不再變化,認(rèn)為結(jié)果已收斂于周期性邊界.表明當(dāng)利用人工彈簧周期性邊界條件時(shí),實(shí)際計(jì)算中只要將剛度系數(shù)的數(shù)值取為ktx=kty=1010N/m2,krx=kry=1010N/rad 時(shí)即可得到收斂的解.

2.2.2 頻散特性對比分析

圖5 給出了本文方法計(jì)算得到的周期開口板彎曲振動頻散曲線,為了對比,同時(shí)給出了利用文獻(xiàn)[41] 中重構(gòu)位移形函數(shù)的方法計(jì)算得到的結(jié)果.其中兩種方法采用的形函數(shù)均為高斯函數(shù),板的幾何和材料參數(shù)一致,伸縮因子、平移因子取值以及掃描波數(shù)均保持一致.

圖5 周期開口板彎曲振動頻散曲線Fig.5 Dispersion curve of bending vibration of periodic plate with opening

由圖5 可以看出,兩種方法計(jì)算得到的結(jié)果基本一致,即0～1500 Hz 范圍內(nèi),圖3 所示的這種周期開口板結(jié)構(gòu)沒有產(chǎn)生完全帶隙和方向帶隙,進(jìn)一步說明了利用人工彈簧模擬周期邊界的準(zhǔn)確性.

2.3 效率對比分析

2.3.1 掃描波數(shù)的影響

掃描波數(shù)決定了頻散曲線上點(diǎn)的個(gè)數(shù),越大的掃描波數(shù)能夠得到更為平滑的頻散曲線.本節(jié)分析了不同掃描波數(shù)下,兩種方法計(jì)算頻散曲線的時(shí)效性,將對比結(jié)果列于表1.分析中保持伸縮因子為6不變.

由表1 可以看出,當(dāng)掃描波數(shù)為30 時(shí),本文方法求解僅需要4.5 s,而傳統(tǒng)能量法需要27.5 s,約為本文方法的6.1 倍;隨著掃描波數(shù)的增大,時(shí)間倍數(shù)越來越大,本文方法效率性體現(xiàn)得更加明顯,當(dāng)掃描波數(shù)為300 時(shí),傳統(tǒng)能量法需要268.2 s,約為本文方法的21.8 倍.這是因?yàn)閭鹘y(tǒng)能量法需要對位移場形函數(shù)進(jìn)行重構(gòu)以滿足周期性邊界條件,這就導(dǎo)致結(jié)構(gòu)的動能和應(yīng)變能矩陣中都含有波數(shù),每掃描一次波數(shù),結(jié)構(gòu)的動能和應(yīng)變能矩陣都要計(jì)算一次,隨著掃描波數(shù)的增大,動能和應(yīng)變能矩陣的計(jì)算次數(shù)也隨之增大,這極大地延長了求解時(shí)間.而本文提出的利用人工彈簧模擬周期性邊界條件的方法,通過將周期性約束條件轉(zhuǎn)變?yōu)閺椥詣菽?僅在有邊界彈性勢能矩陣中包含波數(shù),結(jié)構(gòu)的動能和應(yīng)變能矩陣只需要計(jì)算一次,只有彈性勢能矩陣參與掃描波數(shù)的循環(huán)計(jì)算,這就使得求解時(shí)間大大縮短.

表1 不同掃描波數(shù)下計(jì)算時(shí)效對比Table 1 Comparison of calculation time under different scanning wave number

2.3.2 矩陣維度的影響

前文中提到,伸縮因子控制著解的精度,采用越大的尺度因子使形函數(shù)具有更高的分辨率,但同時(shí)要求更多的計(jì)算資源.平移因子由伸縮因子決定,其數(shù)量控制著形函數(shù)的個(gè)數(shù)和矩陣的維度.本節(jié)分析了不同伸縮因子(矩陣維度) 條件下,兩種方法計(jì)算頻散曲線的時(shí)效性,將對比結(jié)果列于表2.分析中保持掃描波數(shù)為100 不變.

表2 不同矩陣維度下計(jì)算時(shí)效對比Table 2 Comparison of computational effectiveness under different matrix dimensions

傳統(tǒng)能量法需要對位移場形函數(shù)進(jìn)行重構(gòu)以滿足周期性邊界條件,這就導(dǎo)致結(jié)構(gòu)的動能和應(yīng)變能矩陣中都含有波數(shù),結(jié)構(gòu)的動能和應(yīng)變能矩陣需要隨著波數(shù)變化進(jìn)行反復(fù)計(jì)算,隨著伸縮因子取值的增大,平移因子的取值范圍也隨之增大,動能和應(yīng)變能矩陣的維度也隨之增大,這就導(dǎo)致每次循環(huán)計(jì)算動能和應(yīng)變能矩陣的時(shí)間增大.

本文提出的利用人工彈簧定義周期性邊界的方法,僅在邊界彈性勢能矩陣中包含波數(shù),結(jié)構(gòu)的動能和應(yīng)變能矩陣不需要參與波數(shù)循環(huán),只需要計(jì)算一次,僅有邊界彈性勢能矩陣參與隨波數(shù)的循環(huán)計(jì)算.雖然伸縮因子的增大也導(dǎo)致了邊界彈性勢能矩陣維度的增大,但其計(jì)算量遠(yuǎn)小于增大結(jié)構(gòu)的動能和應(yīng)變能矩陣維度的計(jì)算量,這就使得隨著伸縮因子的增大,本文提出的人工彈簧定義周期邊界的方法求解時(shí)間增長量遠(yuǎn)小于傳統(tǒng)能量法:由表2 可以看出,當(dāng)伸縮因子由4 增大到5 時(shí),本文方法求解時(shí)間由3.3 s 延長到3.7 s,求解時(shí)間延長0.4 s;而傳統(tǒng)能量法求解時(shí)間由35.7 s 延長到47.5 s,求解時(shí)間延長11.8 s,約為本文方法的29.5 倍.同時(shí),在編程計(jì)算中,隨著矩陣維度的增大,不僅結(jié)構(gòu)的動能和應(yīng)變能計(jì)算時(shí)間增大,最后特征值求解的時(shí)間也隨之增大,因此在表3 中,當(dāng)伸縮因子由7 增大到8 時(shí),本文方法求解時(shí)間從55.8 s 延長到735.4 s,延長了679.6 s.但相比于傳統(tǒng)能量法,優(yōu)勢還是很明顯,對傳統(tǒng)能量法,伸縮因子由7 增大到8 時(shí),求解時(shí)間從608.5 s 延長到5748.6 s,延長了5140.1 s.因此,相較于傳統(tǒng)能量法,隨著模型復(fù)雜程度的增加以及分析精度要求的提高,結(jié)構(gòu)的質(zhì)量剛度矩陣維度以及掃描波數(shù)點(diǎn)數(shù)隨之增大,本文提出的用人工彈簧模擬周期性邊界的方法在計(jì)算效率上的優(yōu)勢也會越來越明顯.

3 周期組合結(jié)構(gòu)振動帶隙計(jì)算

在鐵路工程、土木工程、機(jī)械工程等領(lǐng)域,組合結(jié)構(gòu)是一種常見的結(jié)構(gòu)類型.這類結(jié)構(gòu)通常也具有周期性,研究這類組合結(jié)構(gòu)的帶隙特性,能夠更全面地理解彈性波在結(jié)構(gòu)中的傳播特性,從而更好地指導(dǎo)其減振降噪設(shè)計(jì).在簡化建模過程中,由于組合結(jié)構(gòu)的復(fù)雜性,可能是由梁-板耦合、多塊板耦合甚至是板-圓柱殼耦合等,這就導(dǎo)致在利用位移形函數(shù)進(jìn)行結(jié)構(gòu)動能和勢能的表征時(shí),形成的質(zhì)量和剛度矩陣維度很大.針對這種大維度矩陣問題,若通過重構(gòu)滿足周期性邊界的位移形函數(shù)的方法來進(jìn)行帶隙計(jì)算,將會面臨相當(dāng)大的計(jì)算成本.因此,本節(jié)以我國CRSTIII 型無砟軌道結(jié)構(gòu)為例,采用人工彈簧定義周期性邊界的方法,計(jì)算了周期性組合結(jié)構(gòu)彎曲振動的帶隙,并與有限元結(jié)果進(jìn)行了比較.

3.1 模型簡化

建立圖6 所示的由鋼軌、扣件、軌道板和CA 砂漿層組成的無限周期結(jié)構(gòu).扣件與CA 砂漿層簡化為彈簧單元;為同時(shí)考慮鋼軌和軌道板的彎曲和剪切性質(zhì),鋼軌簡化為Timoshenko 梁單元,軌道板簡化為Mindlin 板單元,鋼軌選用60 kg/m 鋼軌.本文基于鋼軌和軌道板的周期性與對稱性(假設(shè)結(jié)構(gòu)對稱及載荷形式對稱),可以取一個(gè)周期下結(jié)構(gòu)的一半進(jìn)行建模,模型如圖7 所示.l為單個(gè)軌道板的長度,相鄰軌道板之間是不連續(xù)的,s0為軌道板間縫隙的大小.l+s0為鋼軌周期,其小周期為l0,且l0=(l+s0)/9.c為軌道板寬度的一半,y*為鋼軌距軌道板外側(cè)邊緣的寬度.圖中分別針對鋼軌和軌道板建立了兩個(gè)笛卡爾坐標(biāo)系x′O′y′和xOy.

圖6 CRTSIII 無砟軌道結(jié)構(gòu)縱斷面示意圖Fig.6 Schematic diagram of the longitudinal section of the CRTSIII track structure

圖7 CRTSIII 型無砟軌道結(jié)構(gòu)典型單元示意圖Fig.7 Schematic diagram of CRTSIII track structure of a cell

3.2 運(yùn)動方程的建立

鋼軌和軌道板的位移場可以分別寫為與式(2)類似的形式,即

其中,wb(x′,t) 和θb(x′,t) 分別表示鋼軌的垂向位移和截面轉(zhuǎn)角函數(shù);wp(x,y,t),θpx(x,y,t),θpy(x,y,t) 分別表示軌道板彎曲振動位移、沿x方向轉(zhuǎn)角以及沿y方向轉(zhuǎn)角函數(shù).鋼軌與軌道板位移場形函數(shù)都采用第2 節(jié)提到的高斯函數(shù).

3.2.1 鋼軌的動能和應(yīng)變能

分析中鋼軌考慮為Timoshenko 梁,因此其動能包含了平動動能和轉(zhuǎn)動動能,應(yīng)變能包含了彎曲應(yīng)變能和剪切應(yīng)變能,即

3.2.2 軌道板的動能和應(yīng)變能

3.2.3 扣件彈性勢能與軌道板下彈性地基勢能

模型中9 個(gè)扣件的總彈性勢能可表示為

3.2.4 對稱邊界及周期性邊界條件的定義

圖7 中xOy坐標(biāo)系y=c處邊界為正對稱邊界,可認(rèn)為板邊界上轉(zhuǎn)角為0,位移自由,利用人工彈簧模擬正對稱邊界條件,線彈簧和轉(zhuǎn)角彈簧剛度分別表示為kts和krs,剛度取值參考文獻(xiàn)[43].在x′O′y坐標(biāo)系上x′=0 和x′=l+s0處為Floquet 周期性邊界,由于CRSTIII 型軌道板結(jié)構(gòu)是不連續(xù)的,相鄰軌道板之間存在板縫,因此軌道板沒有Floquet 周期性邊界.整體結(jié)構(gòu)的邊界勢能可表示為

3.2.5 單個(gè)周期的拉格朗日能量泛函

單個(gè)周期結(jié)構(gòu)的總動能和總勢能可表示為

單個(gè)周期結(jié)構(gòu)整體的拉格朗日量可寫為

結(jié)合前文中提到的歐拉-拉格朗日方程,可以將式(34)轉(zhuǎn)化為標(biāo)準(zhǔn)特征值方程問題,即

在第一布里淵區(qū)掃描波矢k=[-π/(l+s0),π/(l+s0)],求解特征值方程(35)就能得到周期梁-板組合結(jié)構(gòu)彎曲振動的頻散曲線.

3.3 算例分析

為了進(jìn)一步驗(yàn)證本文方法的正確性,建立了單個(gè)周期CRSTIII 型無砟軌道結(jié)構(gòu)有限元模型,如圖8所示.鋼軌與軌道板均考慮為實(shí)體單元,扣件與CA砂漿層分別簡化為點(diǎn)彈簧和均布彈簧,其中,鋼軌與軌道板通過扣件彈簧連接,軌道板與基礎(chǔ)通過CA 砂漿均布彈簧連接.軌道板的內(nèi)邊界添加正對稱邊界條件,在鋼軌兩端添加Floquet 周期性邊界條件,其余邊界為自由邊界.鋼軌、軌道板的幾何和材料參數(shù)以及扣件剛度、CA 砂漿均布剛度取值與上述半解析模型一致,鋼軌密度ρb=7850 kg/m3,截面積Ab=77.45 cm2,截面慣性比Ib=3217 cm4,彈性模量Eb=210 GPa,泊松比νb=0.3;軌道比密度ρp=2300 kg/m3,半寬度C=1.25 m,長度l=5.6 m,厚度h=0.21 m,彈性模量Ep=25 GPa,泊松比νb=0.2;軌道板間縫系S0=0.07 m,軌道與FSB 間距y*=0.492 2 m,垂直連接鋼度kf=7.2×107N/m,CA砂漿層的均勻支撐鋼度kca=1.723×108N/m3.選用多物理場仿真軟件Comsol Multiphysics 中的固體力學(xué)模塊進(jìn)行求解,通過設(shè)置掃描波數(shù)求解得到了周期性無砟軌道結(jié)構(gòu)的振動模態(tài),從中挑選出垂向振動模態(tài)進(jìn)行繪制即能得到垂向振動頻散曲線.0～1500 Hz范圍內(nèi),CRSTIII 型無砟軌道結(jié)構(gòu)垂向振動產(chǎn)生了四階帶隙,圖9 給出了兩種方法計(jì)算結(jié)果對比.

圖8 CRTSIII 無砟軌道結(jié)構(gòu)典型單元有限元模型Fig.8 FEM model of typical cell of CRTSIII track structure

圖9 CRTSIII 型無砟軌道結(jié)構(gòu)垂向振動頻散曲線Fig.9 Vertical vibration dispersion curve of CRTSIII track structure

圖9 CRTSIII 型無砟軌道結(jié)構(gòu)垂向振動頻散曲線(續(xù))Fig.9 Vertical vibration dispersion curve of CRTSIII track structure(continued)

由圖9 可以看出,利用人工彈簧定義周期性邊界這種方法進(jìn)行帶隙計(jì)算與有限元計(jì)算結(jié)果吻合度很好,再次證明了該方法的準(zhǔn)確性.目前,針對此類復(fù)雜的周期性組合結(jié)構(gòu)帶隙計(jì)算問題,主要是通過商業(yè)有限元軟件來實(shí)現(xiàn).其中僅有Comsol 等少數(shù)有限元軟件可通過直接設(shè)置Floquet 周期性邊界并掃描波矢來實(shí)現(xiàn)帶隙計(jì)算.然而由于這類商業(yè)軟件目前計(jì)算模塊尚未完善,以Comsol 為例,對于圖6 所示的無砟軌道結(jié)構(gòu)(梁-板組合結(jié)構(gòu))需要采用實(shí)體單元建模,計(jì)算量較大.利用人工彈簧模擬周期性邊界,除周期邊界處人工彈簧彈性勢能矩陣需隨波數(shù)變化而重復(fù)計(jì)算,其他矩陣均只需計(jì)算一次,可以解決計(jì)算量大的問題.

4 結(jié)論

本文基于能量法的基本框架,對傳統(tǒng)能量法進(jìn)行了改進(jìn),提出了利用人工彈簧模擬周期性邊界條件的方法,實(shí)現(xiàn)了周期性邊界條件與位移形函數(shù)的分離.首先以周期開口板為例,驗(yàn)證了本文方法的準(zhǔn)確性,并與傳統(tǒng)能量法比較了計(jì)算效率;然后,針對周期性組合結(jié)構(gòu),以我國CRTSIII 型無砟軌道結(jié)構(gòu)為例,分析了其垂向振動帶隙特性,并與有限元法進(jìn)行了比較.通過本文的研究,可得如下結(jié)論:

(1)與傳統(tǒng)能量法以及有限元法的對比結(jié)果表明了本文方法具有很好的準(zhǔn)確性.同時(shí),相較于傳統(tǒng)能量法,利用人工彈簧模擬周期性邊界條件,無需構(gòu)造滿足周期邊界的位移形函數(shù),進(jìn)一步實(shí)現(xiàn)了周期邊界與位移形函數(shù)的分離,成功解決了傳統(tǒng)能量法中構(gòu)造滿足周期性邊界條件形函數(shù)的困難問題,提高了能量法的適用性.

(2)與傳統(tǒng)能量法相比,利用人工彈簧模擬周期邊界,除周期邊界處彈性勢能矩陣需隨波數(shù)變化而重復(fù)計(jì)算,其他矩陣均只需計(jì)算一次,這極大地提高了計(jì)算效率,實(shí)現(xiàn)了以較低的計(jì)算成本對復(fù)雜周期結(jié)構(gòu)的帶隙進(jìn)行精確計(jì)算.隨著模型復(fù)雜程度的增加以及分析精度要求的提高,結(jié)構(gòu)質(zhì)量、剛度矩陣的維度以及掃描波數(shù)點(diǎn)數(shù)也隨之增多,本文方法在計(jì)算效率上的優(yōu)勢會更加明顯.