基于能量泛函變分原理的無砟軌道垂向振動帶隙分析

馮青松，廖寶亮，郭文杰，楊 舟，付景文，陸建飛

(1.華東交通大學(xué) 鐵路環(huán)境振動與噪聲教育部工程研究中心，江西 南昌 330013；2.江蘇大學(xué) 土木工程與力學(xué)學(xué)院，江蘇 鎮(zhèn)江 212013)

由基本單元沿線路方向周期性排列而成是鐵路軌道結(jié)構(gòu)的基本構(gòu)成形式。當(dāng)振動以彈性波形式在軌道結(jié)構(gòu)中傳播時，滿足Bloch理論和周期性條件，彈性波會形成相應(yīng)的波帶隙(帶隙中的波會迅速衰減以致無法傳播)。實(shí)際應(yīng)用中，帶隙特性對于結(jié)構(gòu)減振或保障結(jié)構(gòu)的可靠性有很大影響。此外，對軌道結(jié)構(gòu)有關(guān)參數(shù)進(jìn)行優(yōu)化設(shè)計，可達(dá)到優(yōu)化其帶隙頻率的位置與寬度的目的，提高其減振效果，從而為軌道結(jié)構(gòu)減振設(shè)計提供新思路。因此，研究周期性無砟軌道結(jié)構(gòu)中彈性波傳播特性對其振動與噪聲控制技術(shù)發(fā)展具有重要意義。

用于求解周期性結(jié)構(gòu)的帶隙特性的方法主要有有限元法、傳遞矩陣法、平面波展開法等。其中，有限元法適用性最強(qiáng)，一般是通過有限元軟件仿真實(shí)現(xiàn)帶隙計算，但計算過程不同于一般的特征值求解，僅有Comsol[1-3]等少數(shù)有限元軟件可通過直接設(shè)置Floquet周期性邊界并掃描波矢實(shí)現(xiàn)帶隙計算。然而，這類商業(yè)軟件計算模塊尚未完善，以Comsol為例，對于無砟軌道結(jié)構(gòu)(梁-板組合結(jié)構(gòu))需要采用實(shí)體單元建模，計算量較大。也有學(xué)者通過自編有限元程序計算周期支撐下鋼軌頻散特性，但尚未考慮軌道板的影響[4]。

傳遞矩陣法[5]對于求解簡單一維結(jié)構(gòu)帶隙特征是十分便捷有效的，如計算考慮周期性彈簧支撐的鋼軌結(jié)構(gòu)帶隙[6-7]。本文建立的無砟軌道模型同時考慮了鋼軌和軌道板，即梁-板組合周期結(jié)構(gòu)，而傳遞矩陣法難于對軌道板建模，因此該方法并不適用。

平面波展開法[8-11]是用于求解結(jié)構(gòu)帶隙的通用方法，但是受限于直接求解微分方程組的分析思路，在處理梁、板連接部位及軌道板邊界條件時有一定困難。此外，也有學(xué)者提出了一些新的周期結(jié)構(gòu)帶隙計算方法，如譜動剛度法[12]、小波法[13]、漸進(jìn)匹配展開法[14]等，這些方法雖各有特點(diǎn)，但處理梁-板組合模型并不適用。

此外，從振動響應(yīng)頻譜特征判斷周期軌道結(jié)構(gòu)中振動波的傳播與衰減[15-21]是目前廣大研究學(xué)者采用的分析其振動特性的思路與方法。但這類研究思路無法從機(jī)理上揭示帶隙的形成過程，且容易受激勵點(diǎn)和測量點(diǎn)位置的影響。

鑒于此，本文以我國CRTSⅢ型無砟軌道結(jié)構(gòu)為研究對象，基于周期結(jié)構(gòu)帶隙原理，提出一種計算周期性軌道結(jié)構(gòu)帶隙特性的新方法，即基于能量泛函變分原理得到各組分的能量泛函，并將求解微分方程邊值問題轉(zhuǎn)化為泛函極值問題，易于解決組合結(jié)構(gòu)的耦合問題。具體來講，首先，利用平面波級數(shù)和布洛赫定理構(gòu)造鋼軌位移場函數(shù)，使其能夠滿足波傳播條件和周期邊界條件，從而得到鋼軌應(yīng)變能和動能；然后，結(jié)合切比雪夫級數(shù)處理軌道板位移場，得到其動能和勢能，接著通過扣件彈簧建立梁與板位移場聯(lián)系，并得到扣件彈性勢能；最后，將所有各分部能量疊加得到總能量泛函，對所有位移場中未知系數(shù)求導(dǎo)后即可得到特征方程，從而求出特征頻率及振型函數(shù)。本文方法解決了傳統(tǒng)解析方法囿于直接求解微分方程組以及難于對軌道板建模的問題，為軌道結(jié)構(gòu)帶隙計算提供新的思路。

1 分析與計算方法

1.1 分析模型

建立雙層彈性支承軌道結(jié)構(gòu)動力模型，軌道結(jié)構(gòu)簡化為由鋼軌、扣件、軌道板和CA砂漿層組成的無限周期結(jié)構(gòu)，見圖1。扣件與CA砂漿層簡化為彈簧單元，為同時考慮鋼軌和軌道板的彎曲和剪切性質(zhì)，鋼軌簡化為Timoshenko梁單元，軌道板簡化為Mindlin板單元。值得一提的是，阻尼對軌道結(jié)構(gòu)振動響應(yīng)有較為顯著的影響，尤其對振動響應(yīng)頻譜曲線中的峰值有明顯削弱。然而，阻尼對結(jié)構(gòu)振動特征頻率影響相對較小，即對帶隙特性影響較弱，因而很多關(guān)于軌道帶隙特性的理論研究中通常選擇忽略阻尼的影響[6-7]。因此，本文研究選擇忽略阻尼。

圖1 無限周期CRTSⅢ型無砟軌道結(jié)構(gòu)

1.2 計算原理

1.2.1 鋼軌和軌道板計算模型

本文基于簡化的CRTSⅢ型無砟軌道模型，取一個周期的結(jié)構(gòu)進(jìn)行建模，模型見圖2。

圖2 鋼軌和軌道板模型

圖2分別針對鋼軌和軌道板建立兩個笛卡爾坐標(biāo)系roB和xoy。l為單個軌道板的長度；s0為軌道板間縫隙大小，l+s0為鋼軌周期，扣件周期為l0，且l0=(l+s0)/9；c為軌道板寬度的一半；y0為鋼軌距軌道板外側(cè)邊緣的寬度。

本方法考慮系統(tǒng)總能量的組成部分有鋼軌應(yīng)變能與動能、軌道板應(yīng)變能與動能、板下CA砂漿層勢能和梁板連接處扣件彈性勢能。

1.2.2 鋼軌能量泛函

取一根鋼軌進(jìn)行分析，為全面分析鋼軌對帶隙特性的影響，將鋼軌考慮為Timoshenko梁。根據(jù)Bloch定理及平面波級數(shù)展開，Timoshenko梁位移場可展開為

(1)

式中：v(r)為垂向位移；β(r)為梁轉(zhuǎn)角，展開項(xiàng)均為(2η+1)項(xiàng)；k為波數(shù)；r為梁的橫坐標(biāo)；i為虛數(shù)單位；ξ(r)為關(guān)于r的試函數(shù)行向量；A1,A2為未知系數(shù)列向量，且

A1=[v-η,v-η+1,…,v0,…,vη]T

A2=[β-η,β-η+1,…,β0,…,βη]T

則(1)式可寫為

(2)

雙軌因形變產(chǎn)生的應(yīng)變能可通過積分得到

(3)

式中：EB為梁的楊氏模量；IB為梁的截面慣性矩；kB為梁的剪切系數(shù)；GB為梁的剪切模量；AB為梁的截面面積。

將(1)式代入(3)式，可得

(4)

式中：A0H為A0的共軛轉(zhuǎn)置向量，A0H={A1H，A2H}；KB為梁的剛度矩陣。

雙軌動能為

(5)

式中：ρB為梁的線密度；ω為系統(tǒng)圓頻率。

將(1)式代入(5)式，可得

(6)

式中：MB為梁的質(zhì)量矩陣。

1.2.3 軌道板能量泛函

考慮為Mindlin板的軌道板位移場可展開為

(7)

式中：w為板的垂向位移；θx為板沿x向轉(zhuǎn)角；θy為板沿y向轉(zhuǎn)角；Amn、Bmn、Cmn為未知系數(shù)；φm(x)、ψn(y)為板的位移試函數(shù)，可選擇具有任意性的切比雪夫級數(shù)[22-23]，其三階微分可導(dǎo)、四階微分連續(xù)的性質(zhì)可以滿足任意邊界條件，且

φ1(x)=1

(8)

(9)

(10)

ψ1(y)=1

(11)

(12)

(13)

式(7)也可寫為

(14)

式中：Amn、Bmn、Cmn均為M×N維列向量，且Amn=[A11,A12,…,Amn]T，Bmn=[B11,B12,…,Bmn]T,Cmn=[C11,C12,…,Cmn]T;λ=[φ1(x)ψ1(y),φ1(x)ψ2(y),…,φM(x)ψ1(y),…,φM(x)ψN(y)]=[Q1(x,y),Q2(x,y),…,QMN(x,y)]，對其中任意元素，有φm(x)·ψn(y)=Qa(x,y)，且a=(m-1)N+n，m∈[1,M]，n∈[1,N]。

軌道板因形變而產(chǎn)生彈性體運(yùn)動，進(jìn)而產(chǎn)生的應(yīng)變能可由二重積分得到，可表示為

(15)

軌道板因振動產(chǎn)生的剛體運(yùn)動而具有動能，其表達(dá)式為

(16)

式中：TP為軌道板的動能。

1.2.4 軌道板下CA砂漿層勢能

軌道板由其板下CA砂漿層均布支承，同時結(jié)合軌道板垂向位移可以通過二重積分得到支承層勢能，其表達(dá)式為

(17)

式中：U1為支承層勢能；ks為板下CA砂漿層平均剛度。

1.2.5 扣件彈性勢能

由扣件垂向剛度和扣件處鋼軌與軌道板的垂向位移差值可求出扣件彈性勢能。模型包含對稱的兩根鋼軌，因此，扣件的總彈性勢能U2可表示為

(18)

軌-板位移差值為

(19)

將式(19)代入式(18)，可得

(20)

式中：AH為A的共軛轉(zhuǎn)置向量，AH={A1H,A2H,AmnT,BmnT,CmnT}為未知系數(shù)列向量；j為扣件序列編號；j0為單軌扣件個數(shù)；kz為扣件垂向剛度；01為 (2η+1)×(2η+1)維零矩陣；02為(2η+1)×MN維零矩陣；03為MN×(2η+1)維零矩陣；04為MN×MN維零矩陣。

1.2.6 總能量泛函

系統(tǒng)總能量泛函Π為

Π=UB+UP+U1+U2-TB-TP

(21)

對未知系數(shù)求極值，即

(22)

于是，結(jié)構(gòu)振動問題轉(zhuǎn)化為求解特征值問題，可表示為

(K-ω2M)A=0

(23)

式中：K為系統(tǒng)總剛度矩陣；M為系統(tǒng)總質(zhì)量矩陣；A為未知的系數(shù)列向量；ω為圓頻率。

2 算例分析

2.1 收斂性分析

本文將位移進(jìn)行級數(shù)展開，為得到頻率結(jié)果，計算時需將位移級數(shù)進(jìn)行截斷。截斷項(xiàng)個數(shù)的取值不同，將影響計算精度和效率，因此有必要對不同截斷項(xiàng)數(shù)下的頻率進(jìn)行收斂性分析。

本節(jié)采用控制變量法研究文中涉及的3個截斷項(xiàng)數(shù)M、N、η的取值對帶隙頻率的影響。我國CRTSⅢ型軌道板長為5 600 mm，板寬為2 500 mm，板厚為210 mm，板縫寬為70 mm，鋼軌為60 kg/m鋼軌，其余軌道結(jié)構(gòu)參數(shù)取值見表1。其中，CA砂漿剛度為換算后的平均支承剛度。從頻散曲線中隨機(jī)取3條曲線波數(shù)位于π/l處的頻率作為研究對象，其結(jié)果見圖3。從圖3可以看出，當(dāng)M=20，N=η=8時，帶隙頻率基本收斂。因此，將其作為后面的算例分析中截斷項(xiàng)數(shù)取值。

表1 無砟軌道結(jié)構(gòu)參數(shù)

圖3 帶隙頻率收斂性分析

2.2 將軌道板考慮為質(zhì)量塊的垂向振動頻散特性分析

國內(nèi)外現(xiàn)有涉及無砟軌道振動特性的研究中，或忽略軌道板的影響[7]，或?qū)④壍腊蹇紤]為梁[24]，鮮有研究將其作為板進(jìn)行建模。因此，在運(yùn)用本文方法之前，有必要在原方法的基礎(chǔ)上，不考慮軌道板形變將其簡化為具有無限剛度的有限長質(zhì)量塊，對模型垂向振動頻散特性進(jìn)行分析。將結(jié)果作為對照組，以探索本文將軌道板考慮為Mindlin板的分析方法的優(yōu)越性，其帶隙結(jié)果見圖4。為便于分析計算，引入?yún)?shù)m，用于代替波矢k，且滿足

(24)

式中：k表示被限制在第一不可約Brillioun區(qū)內(nèi)的波矢。

軌道結(jié)構(gòu)可看成沿軌道縱向的一維周期結(jié)構(gòu)，對于一維周期結(jié)構(gòu)，其第一不可約Brillioun范圍即k的取值范圍為[-π/l,π/l]，由式(24)可知，參數(shù)m的取值為[-1,1]。

從圖4中可以看出，將軌道板考慮為質(zhì)量塊時，在0～1 200 Hz范圍內(nèi)，結(jié)構(gòu)共存在三階垂向振動帶隙，頻率范圍分別為0～89.7、92.7～217、1 002～1 046 Hz。

圖4 軌道板考慮為質(zhì)量塊時結(jié)構(gòu)帶隙結(jié)果

2.3 將軌道板考慮為板的垂向振動頻散特性分析

將軌道板考慮為Mindlin板時，運(yùn)用理論法計算得到0～1 200 Hz范圍內(nèi)無砟軌道結(jié)構(gòu)垂向振動頻散曲線見圖5。

由圖5可見，無砟軌道結(jié)構(gòu)垂向振動會產(chǎn)生四階帶隙，前三階在低頻產(chǎn)生(0～89.2、104.8.5～115.5、119.5～195.4 Hz)，第四階在高頻產(chǎn)生(1 001.7～1 032.0 Hz)，其中，第一階、第三階帶隙帶寬較大，分別為89.2、75.9 Hz。

圖5 垂向振動頻散曲線

從結(jié)果可以看出，將軌道板考慮為Mindlin板和質(zhì)量塊所得到的帶隙結(jié)果基本吻合，區(qū)別在于考慮為質(zhì)量塊時得到的二階帶隙頻率為92.7～217 Hz，而考慮為Mindlin板時在這一頻率范圍內(nèi)生成了兩條帶隙，分別為104.8～115.5、119.5～195.4 Hz。這是因?yàn)閷④壍腊蹇紤]為Mindlin板后，軌道板不僅存在剛體運(yùn)動，還存在因形變而產(chǎn)生的彈性體運(yùn)動，多種運(yùn)動耦合后在該頻段范圍內(nèi)產(chǎn)生通帶，從而將一條帶隙分割為兩條，因此該頻段帶隙的產(chǎn)生受軌道板影響較大。此外，這一結(jié)果也表明，更為細(xì)致地考慮軌道板在結(jié)構(gòu)振動中的表現(xiàn)，可更全面地展現(xiàn)軌道結(jié)構(gòu)的帶隙特性。

為進(jìn)一步驗(yàn)證前述理論分析的正確性，建立單個周期無砟軌道結(jié)構(gòu)有限元模型，見圖6。利用對稱性建模，將鋼軌與軌道板均考慮為實(shí)體單元，扣件與道砟簡化為連接彈簧，其中，鋼軌與軌道板通過扣件彈簧連接，軌道板與基礎(chǔ)通過道砟彈簧連接。將有限元模型劃分為有限個通過節(jié)點(diǎn)連接的單元，單元類型為Langrage-quadratic，為滿足計算精度需求，將其劃分為521 655個域單元、138 250個邊界元以及24 065個邊單元，求解總自由度為2447 210個。選用多物理場仿真軟件Comsol Multiphysics中的固體力學(xué)模塊進(jìn)行求解，在圖6中板的內(nèi)側(cè)邊界添加對稱邊界，其余為自由邊界，同時鋼軌右端截面定義為源邊界，左端截面定義為目標(biāo)邊界，在鋼軌源邊界添加floquet周期性邊界條件，計算機(jī)配置為2.7 GHz CPU，192 GB內(nèi)存，計算時長約14 h。理論法采用Matlab進(jìn)行編程計算，耗時約32 s，表明本文方法計算效率較高。圖7～圖9給出了解析法與有限元法得到的垂向振動頻散曲線局部對比分析圖，不同方法計算得到的帶隙頻率對比見表2。

圖6 雙層無砟軌道結(jié)構(gòu)有限元模型

圖7 軌道板考慮為板時第一、二階帶隙對比分析

圖8 軌道板考慮為板時第三階帶隙對比分析

圖9 軌道板考慮為板時第四階帶隙對比分析

表2 不同方法計算得到的帶隙頻率對比

由表2可見，比較四階帶隙頻率處頻散曲線可以發(fā)現(xiàn)，理論解和有限元法計算得到的帶隙頻率位置基本一致，整體頻散曲線也基本吻合，說明本文方法不僅計算效率高，能保證計算結(jié)果的可靠性，還表明0～1 200 Hz范圍內(nèi)，采用Timoshenko梁與Mindlin板耦合模型能夠較為準(zhǔn)確地表征雙層周期性無砟軌道結(jié)構(gòu)的帶隙特性。

3 帶隙形成機(jī)理分析

本文2.3節(jié)從無砟軌道結(jié)構(gòu)的帶隙計算結(jié)果中選出了位于0～1 200 Hz范圍內(nèi)的四條帶隙，其中前三階帶隙為局域共振帶隙，第四階帶隙為Bragg帶隙。

鋼軌和周期性排列的局域共振單元是局域共振帶隙形成的主要條件。局域共振帶隙的起始與截止頻率對應(yīng)振動模態(tài)可以用“質(zhì)量-彈簧”模型描述[25-27]。本文建立了與雙層無砟軌道結(jié)構(gòu)相對應(yīng)的“質(zhì)量-彈簧”模型，以進(jìn)一步分析帶隙中局域共振帶隙的產(chǎn)生機(jī)理，見圖10。其中，k1為CA砂漿支承總剛度；k2為扣件垂向總剛度；m1和m2分別為軌道板和鋼軌質(zhì)量，且滿足關(guān)系式：k1=2kslc，m1=2ρPlch，k2=18kz，m2=2ρBAB·(l+s0)。

圖10 雙層無砟軌道結(jié)構(gòu)局域共振帶隙分析模型

由圖10(a)可得局域共振帶隙的起始頻率fs為

fs1=0

(25)

(26)

由圖10(b)可得其特征頻率ω滿足關(guān)系式為

(27)

從而可以得到局域共振帶隙截止頻率fc為

(28)

(29)

(30)

(31)

將各項(xiàng)軌道參數(shù)代入式(26)和式(28)中，得到結(jié)構(gòu)局域共振帶隙起始頻率和截止頻率，從而可估算結(jié)構(gòu)局域共振帶隙頻率范圍為：0～89.3 Hz，111.4～193.7 Hz。估算結(jié)果與表2基本一致。需要強(qiáng)調(diào)的是，考慮為Mindlin板后，軌道板不僅存在剛體運(yùn)動，還存在因形變而產(chǎn)生的彈性體運(yùn)動，多種運(yùn)動耦合后在該頻段范圍內(nèi)產(chǎn)生通帶，從而將其中一條低頻帶隙分割為兩條，證明了2.3中兩種方法得到的前三階帶隙為局域共振帶隙，再次證明本文方法求解軌道結(jié)構(gòu)帶隙特性的準(zhǔn)確性。

根據(jù)Bragg散射機(jī)理[28]，當(dāng)彎曲波以波長λ在周期性軌道結(jié)構(gòu)中傳播時，此時晶格常數(shù)為l0(即扣件間距)，則Bragg條件可表示為：2l0=nλ(n=1,2,…)。當(dāng)n取最小值1時，對應(yīng)的彎曲波頻率為第一階“pinned-pinned”頻率的起始頻率，此時位于扣件處的鋼軌垂向位移為零，而跨中處鋼軌垂向位移達(dá)到最大值[29]，振型見圖11，且可設(shè)鋼軌的位移場滿足關(guān)系：

圖11 第一階“pinned-pinned”頻率振型

v=H1sin(nπx/l0)

(32)

β=H2cos(nπx/l0)

(33)

式中：H1和H2為未知系數(shù)，且不為零。

鋼軌的振動微分方程為[30]

將式(32)和式(33)代入方程組(34)，可得

(35)

要使式(35)成立，則有

(36)

將式(36)行列式展開，代入n=1及各項(xiàng)軌道結(jié)構(gòu)參數(shù)可解得：f1=1 001.7 Hz，f2=6 598.3 Hz。其中，f1為以彎曲位移為主對應(yīng)的“pinned-pinned”頻率，f2為以剪切位移為主對應(yīng)的“pinned-pinned”頻率。因此，第一階“pinned-pinned”頻率的起始頻率為1 001.7 Hz，與表2中各方法所求最后一階帶隙頻率起始頻率相吻合，證明第四階帶隙為Bragg帶隙。從式(36)可以看出，Bragg帶隙的起始頻率僅由扣件間距和鋼軌物理性質(zhì)所決定，這也可以解釋將軌道板由考慮為質(zhì)量塊轉(zhuǎn)變?yōu)榭紤]成Mindlin板后，第四階帶隙起始頻率基本不變的現(xiàn)象。此外，Bragg帶隙的帶寬主要由扣件剛度決定[31]。

4 結(jié)論

本文以我國CRTSⅢ型無砟軌道結(jié)構(gòu)為研究對象，提出了一種基于能量泛函變分原理和平面波級數(shù)展開的混合方法，用于計算周期性無砟軌道結(jié)構(gòu)振動帶隙。理論計算結(jié)果與有限元仿真結(jié)果和簡化模型結(jié)果進(jìn)行了對比分析，得到如下結(jié)論：

(1)本文提出的基于能量泛函變分原理的新方法在保證計算準(zhǔn)確性的前提下，極大地提高了計算效率(同模型下，計算速度約為有限元法的上千倍)。本文方法克服了傳統(tǒng)解析方法需要直接求解微分方程組以及難以對完整軌道板進(jìn)行建模的問題，可更為全面地分析軌道板模態(tài)對結(jié)構(gòu)振動特性的影響。

(2) 0～1200 Hz頻率范圍內(nèi)，用Timoshenko梁和Mindlin板能較準(zhǔn)確地表征周期無砟軌道結(jié)構(gòu)帶隙特性，且垂向振動產(chǎn)生0～89.2、104.8.5～115.5、119.5～195.4、1 001.7～1 032.0 Hz四階帶隙。

(3) 與將軌道板考慮為質(zhì)量塊相比，考慮為Mindlin板后，因軌道板存在多種模態(tài)，使得其不僅存在剛體運(yùn)動，還存在因形變而產(chǎn)生的彈性體運(yùn)動，多種運(yùn)動耦合后產(chǎn)生通帶，從而將在該頻段范圍內(nèi)一條帶隙分割為多條，可更為細(xì)致地考慮軌道板在結(jié)構(gòu)振動中的表現(xiàn)，更全面地展現(xiàn)軌道結(jié)構(gòu)的帶隙特性。

(4) 文中求得帶隙中前三階為局域共振帶隙，第四階為Bragg帶隙。這兩類帶隙頻率均可通過估算得到其大致范圍，其中，局域共振帶隙頻率范圍可由“質(zhì)量-彈簧”模型估算得到，Bragg帶隙頻率范圍可由不同的“pinned-pinned”頻率振型分析計算得到。此外，Bragg帶隙起始頻率僅由扣件間距和鋼軌物理性質(zhì)決定。