基于傳遞矩陣法的軌道結(jié)構(gòu)動力分析方法

朱志輝，王盈瑩，龔 威，張 磊，蔣麗忠

(1.中南大學(xué) 土木工程學(xué)院，湖南 長沙 410075；2.中南大學(xué) 高速鐵路建造技術(shù)國家工程實驗室，湖南 長沙 410075)

軌道結(jié)構(gòu)作為鐵路系統(tǒng)中最重要的結(jié)構(gòu)之一，具有引導(dǎo)列車運行方向與傳遞列車荷載至下部結(jié)構(gòu)的作用[1]。作為列車的直接支撐結(jié)構(gòu)，軌道結(jié)構(gòu)振動響應(yīng)直接影響輪軌相互作用力、列車行車安全性與乘坐舒適度的評估[2]。因此考慮軌道結(jié)構(gòu)的柔性支撐作用，建立準(zhǔn)確的軌道結(jié)構(gòu)動力學(xué)模型一直是列車-軌道-橋梁耦合動力學(xué)理論研究重要部分。

由于軌道結(jié)構(gòu)中存在局部與高頻振動[3-4]，為避免模態(tài)疊加法在有效振型選擇上的困難，部分學(xué)者基于精細(xì)有限元模型，采用有限元直接剛度法(Direct Stiffness Method，DSM)建立其動力方程[5-6]。這種方法雖然可以準(zhǔn)確計算軌道結(jié)構(gòu)的局部與高頻振動[7]，但涉及超長軌道結(jié)構(gòu)動力分析時，DSM不可避免會產(chǎn)生大量結(jié)構(gòu)自由度，使得動力方程系數(shù)矩陣規(guī)模龐大，導(dǎo)致計算困難[8]。針對該問題，移動軌道技術(shù)[9-10]被用于降低軌道結(jié)構(gòu)自由度數(shù)目，該方法通過保留荷載作用處與附近一段長度范圍內(nèi)的軌道結(jié)構(gòu)，刪除其余非荷載作用處軌道結(jié)構(gòu)以達(dá)到降低自由度數(shù)目的目的。但該方法在引入新增軌道結(jié)構(gòu)時，由于新增軌道和原有軌道結(jié)構(gòu)在連接處存在位移不協(xié)調(diào)問題，從而導(dǎo)致新增單元處動力響應(yīng)突變的問題[11]。

由于軌道結(jié)構(gòu)具有明顯的周期特征，很多學(xué)者針對其頻域動力分析提出了周期Floquet變換法[12]、周期Fourier模態(tài)法[13]等求解方法，有效解決了軌道結(jié)構(gòu)頻響分析問題。但該類頻域方法無法考慮軌道構(gòu)成的差異性及非線性特征。因此，本文利用傳遞矩陣方法[14]，建立周期性軌道結(jié)構(gòu)模型，并在時域范圍內(nèi)求解其動力響應(yīng)，以更好地解決前述問題。

傳遞矩陣法(Transfer Matrix Method，TMM)的基本思想是針對周期性結(jié)構(gòu)，通過將其重復(fù)性部分劃分為若干元胞，建立元胞內(nèi)部及相鄰元胞間的狀態(tài)矢量傳遞關(guān)系，實現(xiàn)整個周期性結(jié)構(gòu)的求解。因此TMM并不適用于非周期性結(jié)構(gòu)，但對于周期性結(jié)構(gòu)，由于TMM無需建立結(jié)構(gòu)的整體動力方程，避免了大型線性方程組的求解[15]；同時，由于元胞本身自由度非常小，因此該方法可有效減少結(jié)構(gòu)的自由度數(shù)目，降低計算過程中的內(nèi)存占用，并使得模型的建立更為便捷[16]，能較好地解決由于軌道結(jié)構(gòu)自由度數(shù)目過多導(dǎo)致的計算困難和計算效率低下的問題。文獻(xiàn)[17-20]采用該方法，分別針對周期性的彈性支承梁、桁架結(jié)構(gòu)、框架結(jié)構(gòu)和圓柱殼體進(jìn)行了研究，取得了良好的效果。

上述研究中，系統(tǒng)傳遞關(guān)系通過傳遞矩陣的連乘方式建立，廣泛用于結(jié)構(gòu)靜力及動力特性問題的分析求解。但對于結(jié)構(gòu)瞬態(tài)動力分析，由于傳遞矩陣中考慮了等效剛度和等效荷載，且相較于剛度矩陣和荷載向量二者元素更大，因此連乘算法中的累積誤差被加劇，導(dǎo)致時域數(shù)值積分過程中存在計算不穩(wěn)定現(xiàn)象[21]。針對該問題，Horner等[22]提出了適用于鏈?zhǔn)浇Y(jié)構(gòu)的Riccati傳遞矩陣法，通過傳遞矩陣的遞推建立系統(tǒng)傳遞關(guān)系，有效提高了TMM的數(shù)值穩(wěn)定性。Xue[14]在相鄰元胞交界處引入了剛度方程假定，克服了上述方法元胞邊界自由度相等的限制，將其擴(kuò)展至一般結(jié)構(gòu)。上述關(guān)于周期性結(jié)構(gòu)的研究，很好地推動了TMM的推廣應(yīng)用，但其在軌道結(jié)構(gòu)中多用于波傳播特性分析[23]，涉及時域分析的研究還較少。

本文針對具有顯著周期性特征、且通常需要考慮較長長度的軌道結(jié)構(gòu)，將軌道結(jié)構(gòu)的周期性重復(fù)部分劃分為若干類元胞結(jié)構(gòu)，基于鏈?zhǔn)浇Y(jié)構(gòu)的Riccati傳遞矩陣法，確定元胞內(nèi)部及相鄰元胞間的傳遞關(guān)系，從而建立周期性軌道結(jié)構(gòu)傳遞矩陣模型。以有砟軌道和CRTSⅡ型板式無砟軌道為例，基于本文方法分別提出了對應(yīng)的軌道結(jié)構(gòu)元胞劃分方案，通過TMM計算所得移動荷載引起的有砟軌道結(jié)構(gòu)動力響應(yīng)結(jié)果，對比文獻(xiàn)數(shù)據(jù)驗證了本文方法的準(zhǔn)確性。以移動荷載列通過無砟軌道結(jié)構(gòu)為例，對比TMM和DSM計算結(jié)果驗證了本文方法的高效性。

1 時域傳遞矩陣法

1.1 元胞單元時域傳遞關(guān)系

以周期支承的平面梁為例，采用傳遞矩陣法時，可將其劃分為A、B、C三類元胞結(jié)構(gòu)，其中，第1個元胞為A類，第2～n-1個元胞為B類，第n個元胞為C類。對于元胞B，由于相鄰兩元胞包含同一支承結(jié)構(gòu)，其支承剛度為原始剛度的1/2；對于元胞A，由于其處于整個結(jié)構(gòu)的輸入端，故元胞輸入端支承剛度保持不變，輸出端支承剛度為原始剛度的1/2；同理，元胞C輸入端支承剛度為原始剛度的1/2，輸出端支承剛度與原始取值相同，見圖1。

圖1 周期性結(jié)構(gòu)的元胞結(jié)構(gòu)劃分示意

任意取其中第i個元胞，則t+Δt時刻其運動方程為

(1)

采用Newmark-β法對上述方程進(jìn)行求解，則有

(2)

圖2 元胞i輸入端I與輸出端O向量傳遞關(guān)系示意

則對于元胞i有

(3)

將其代入式(2)可得

(4)

(5)

展開式(5)，則有

(6)

(7)

考慮位移和內(nèi)力的連續(xù)性[25]，則剖面兩側(cè)的位移相等，內(nèi)力相反，故相鄰兩元胞間狀態(tài)矢量的傳遞關(guān)系可表示為

Ui,I=Ui-1,ONi,I=-Ni-1,O

(8)

類似狀態(tài)矢量的廣義Riccita變換[22]，假設(shè)聯(lián)系元胞i輸入端內(nèi)力與位移向量的廣義剛度方程為

Ni,I=SiUi,I+Eii>2

(9)

式中：Si和Ei分別為第i個剖面的剛度方程系數(shù)矩陣和等效外荷載向量。

將式(8)和式(9)代入式(6)，可得元胞i+1輸入端內(nèi)力與位移向量的廣義剛度方程為

(10)

將式(10)代入式(7)，可得元胞i+1輸入端內(nèi)力與位移向量的廣義剛度方程為

Ni+1,I=Si+1Ui+1,I+Ei+1

(11)

(12)

(13)

式(12)和式(13)即為相鄰剖面間，剛度方程系數(shù)矩陣和等效外荷載向量的傳遞關(guān)系。

若S2和E2已知，則通過上述兩式可依次解得各剖面的Si和Ei，進(jìn)而建立整體結(jié)構(gòu)的狀態(tài)矢量傳遞模型(以下簡稱為“傳遞模型”)。結(jié)合式(9)可知，第n+1個剖面處(結(jié)構(gòu)輸出端)內(nèi)力與位移的廣義剛度方程可表示為

Nn+1,I=Sn+1Un+1,I+En+1

(14)

式中：n為系統(tǒng)元胞個數(shù)。

考慮結(jié)構(gòu)輸出端邊界條件，通過式(14)可得該處位移響應(yīng)。將式(14)代入式(10)，由輸出端至輸入端依次求解，可得各元胞的位移響應(yīng)，并通過Newmark-β方法解得其速度和加速度響應(yīng)。

1.2 傳遞系數(shù)求解

式(11)中，令i=1，可得

N2,I=S2U2,I+E2

(15)

將式(8)代入式(15)可得

N1,O=-S2U1,O-E2

(16)

對于元胞1，擴(kuò)展式(5)可得

(17)

(18)

考慮結(jié)構(gòu)輸入端邊界條件，結(jié)合式(17)、式(18)，可解得U1,O和N1,O間的關(guān)系，進(jìn)而對照式(15)可解得對應(yīng)邊界條件下的S2和E2。

對于位移邊界條件，U1,I是已知的，因此由式(18)可得元胞1的輸出端內(nèi)力為

(19)

對照式(16)可知

(20)

(21)

對于荷載邊界條件，N1,I是已知的，因此由式(17)可得元胞1的輸入端位移為

(22)

將式(22)代入式(18)，可得元胞1的輸出端內(nèi)力為

(23)

對照式(16)可知

(24)

(25)

(26)

類似地，展開式(26)，求解N1,O和U1,O間的關(guān)系，對照式(16)可得

(27)

(28)

1.3 計算流程

根據(jù)上述理論，基于Matlab平臺編制了外荷載作用下的周期結(jié)構(gòu)瞬態(tài)響應(yīng)的TMM求解程序，主要計算流程見圖3。

圖3 傳遞矩陣法計算流程

2 軌道結(jié)構(gòu)的傳遞矩陣模型

目前，我國鐵路中鋪設(shè)有多種類型的軌道結(jié)構(gòu)，其結(jié)構(gòu)形式和力學(xué)特征均有所不同。本節(jié)分別針對有砟軌道和CRTSⅡ型板式無砟軌道，提出了對應(yīng)的元胞劃分方案，并基于TMM方法建立了元胞間的動力傳遞關(guān)系。為簡便起見，采用垂向軌道模型對方法進(jìn)行介紹，只考慮軌道結(jié)構(gòu)在豎向平面內(nèi)的動力響應(yīng)。

2.1 有砟軌道元胞模型及傳遞關(guān)系

有砟軌道由鋼軌、扣件、軌枕和道砟組成，其垂向整體模型見圖4(a)，單跨模型(相鄰兩軌枕間軌道結(jié)構(gòu))及其節(jié)點自由度如圖4(b)。其中，鋼軌視作連續(xù)彈性離散點支撐的垂向Timoshenko梁模型，考慮豎向位移Zr及轉(zhuǎn)動角θr，長度為L；軌枕(及道砟質(zhì)量)等效為集中質(zhì)量ms[30]，考慮豎向位移Zs，相鄰兩軌枕間距為l；扣件支撐和道砟支撐等效為彈簧-阻尼器，考慮豎向剛度系數(shù)kf、kb及阻尼系數(shù)cf、cb。

圖4 有砟軌道垂向模型示意

有砟軌道以軌枕間距l(xiāng)為周期，其構(gòu)成具有重復(fù)性，因此應(yīng)用傳遞矩陣法時，僅需對單跨軌道進(jìn)行劃分。為對鋼軌進(jìn)行較為精細(xì)的劃分，如圖5所示，將單跨軌道劃分為3類元胞結(jié)構(gòu)：元胞A為左側(cè)支撐結(jié)構(gòu)，由長度為l1的鋼軌、軌枕、扣件和道砟組成；元胞B為中間連接結(jié)構(gòu)，由長度為l2的鋼軌組成；元胞C為右側(cè)支撐結(jié)構(gòu)，除鋼軌長度為l3外，其余構(gòu)成與元胞A相同，且二者參數(shù)均為整體軌道結(jié)構(gòu)的1/2(軌道輸入及輸出端處理方式與1.1節(jié)相同，由于整個系統(tǒng)僅包含2個端部元胞，數(shù)量極少，不再單獨列出，下同)。由圖可知，元胞A、B和C的輸入/輸出端節(jié)點自由度數(shù)目分別為3/2、2/2、2/3。

圖5 有砟軌道結(jié)構(gòu)劃分示意

根據(jù)相鄰元胞結(jié)構(gòu)間的傳遞關(guān)系，進(jìn)而建立整體軌道結(jié)構(gòu)的傳遞模型。有砟軌道傳遞模型中存在A?B、B?B、B?C和C?A共4種傳遞關(guān)系，其中，求解傳遞系數(shù)時方向由輸入端至輸出端(A→B，B→B，B→C，C→A)，求解響應(yīng)時方向由輸出端至輸入端(A→C，C→B，B→B，B→A)。以A?B為例，假設(shè)第i和i+1個元胞的類型分別為A和B，則有

(29)

對于由梁、質(zhì)量塊、線性彈簧阻尼器等簡單構(gòu)件組成的元胞結(jié)構(gòu)，可利用有限元理論直接推導(dǎo)其系數(shù)矩陣[15]，對于由薄板、實體等復(fù)雜構(gòu)件組成的元胞結(jié)構(gòu)，可在通用軟件(如Ansys)中建立有限元模型，導(dǎo)出其系數(shù)矩陣[26]。

元胞i輸入和輸出端位移和內(nèi)力可表示為

(30)

式中：Z和θ分別為位移和轉(zhuǎn)角；Q和M分別為剪力和彎矩；下標(biāo)r和s分別為鋼軌和軌枕。

將式(29)、式(30)代入式(5)，可得元胞i的動力方程為

(31)

將式(31)代入式(12)和式(13)，可實現(xiàn)建模過程中元胞i→i+1(A→B)等效剛度及荷載的傳遞，即

(32)

(33)

軌道結(jié)構(gòu)兩端無內(nèi)力作用，適用于荷載邊界條件，式(32)和式(33)中Si和Ei可在由式(24)和式(25)解得傳遞系數(shù)S2和E2后，循環(huán)應(yīng)用式(12)和式(13)計算得到。

將式(32)和式(33)代入式(10)，可得求解過程中元胞i+1→i(B→A)位移響應(yīng)的傳遞，即

(34)

各元胞位移響應(yīng)可在由式(14)解得邊界處的響應(yīng)后，循環(huán)應(yīng)用式(34)計算得到。至此，已完成元胞A?B傳遞關(guān)系的建立，其他3種傳遞關(guān)系也可通過上述方式導(dǎo)得。結(jié)合軌道結(jié)構(gòu)中各類元胞結(jié)構(gòu)的排列順序，應(yīng)用對應(yīng)的傳遞關(guān)系，可實現(xiàn)相鄰元胞間狀態(tài)矢量的傳遞，從而完成有砟軌道的建模和求解。

由上述推導(dǎo)過程可知，采用傳遞矩陣法對有砟軌道結(jié)構(gòu)進(jìn)行動力分析時，無需建立其整體模型。僅需對重復(fù)構(gòu)成的單跨軌道進(jìn)行拆分，并依據(jù)各類元胞結(jié)構(gòu)間的傳遞關(guān)系和其在整體結(jié)構(gòu)中的排列順序，建立整體軌道結(jié)構(gòu)的傳遞模型即可。對于垂向有砟軌道模型，無論軌道長度L如何變化，計算過程中所涉及的最大矩陣階數(shù)始終為3。

2.2 無砟軌道元胞模型及傳遞關(guān)系

與傳統(tǒng)的有砟軌道相比，無砟軌道具有結(jié)構(gòu)穩(wěn)定性高、剛度均勻性好等突出優(yōu)點[27-28]。CRTSⅡ型板式無砟軌道是我國廣泛應(yīng)用的無砟軌道結(jié)構(gòu)形式之一，京津城際、京滬高鐵等均采用了此種軌道結(jié)構(gòu)[29]。本文以CRTSⅡ型縱連板式無砟軌道垂向模型為例，介紹基于傳遞矩陣法的無砟軌道動力分析方法。

CRTSⅡ型板式無砟軌道垂向模型由鋼軌、扣件、預(yù)制混凝土軌道板、水泥乳化瀝青調(diào)整層(CA砂漿)和混凝土底座板等部分組成，見圖6(a)。其中，鋼軌、軌道板和底座板視作連續(xù)彈性離散點支撐的Timoshenko梁模型，長度為Lr，考慮豎向位移Zr、Zt、Zc及轉(zhuǎn)動角θr、θt、θc；扣件、CA砂漿和基礎(chǔ)以離散分布的黏滯阻尼和線性彈簧表示，間距為l，考慮豎向剛度系數(shù)kf、kc、kb和阻尼系數(shù)cf、cc、cb，見圖6(b)。

圖6 CRTSⅡ型板式無砟軌道垂向模型示意

CRTSⅡ型板式無砟軌道可認(rèn)為由圖6(b)所示的單跨軌道重復(fù)構(gòu)成，因此，應(yīng)用傳遞矩陣法時，僅需要對單跨軌道進(jìn)行劃分。如圖7所示，類似有砟軌道，將其劃分為A、B和C三類元胞結(jié)構(gòu)。其中，元胞A和C均由鋼軌、軌道板、底座板、扣件、CA砂漿和地基組成，長度分別為l1和l2，連接件參數(shù)取值相同，均為整體模型的1/2；元胞B由鋼軌、軌道板和底座板組成，長度為l3。由圖6可知，元胞A、B和C的輸入/出端的節(jié)點自由度數(shù)目均為6/6。

圖7 CRTSⅡ型板式無砟軌道劃分示意

類似式(34)，元胞i和i+1間位移響應(yīng)的傳遞關(guān)系可表示為

(35)

由上述推導(dǎo)過程可知，采用傳遞矩陣法對無砟軌道進(jìn)行動力分析時，同樣無需建立整體模型。僅需對單跨軌道進(jìn)行劃分，并對各類元胞結(jié)構(gòu)進(jìn)行計算即可。對于垂向無砟軌道模型，無論軌道長度L如何變化，計算過程中所涉及的最大矩陣階數(shù)始終為6。

3 移動荷載作用下軌道結(jié)構(gòu)動力響應(yīng)分析

為驗證上述關(guān)于軌道結(jié)構(gòu)傳遞矩陣算法的正確性，本節(jié)以移動荷載作用下的有砟軌道結(jié)構(gòu)為例，對比了TMM、DSM和文獻(xiàn)中的鋼軌動力響應(yīng)計算結(jié)果。并以自由度數(shù)目較多的空間無砟軌道結(jié)構(gòu)為例，對比了TMM和DSM計算所得的移動荷載作用下軌道結(jié)構(gòu)動力響應(yīng)結(jié)果，進(jìn)一步驗證了本文方法的高效性。

3.1 有砟軌道動力響應(yīng)分析

文獻(xiàn)[30]給出了移動荷載F=1 kN以速度60 km/h移動時的跨中(坐標(biāo)x=60.3 m)鋼軌位移時程曲線，本文通過TMM和DSM計算相同的算例，將得出的結(jié)果與文獻(xiàn)結(jié)果進(jìn)行對照，以驗本文方法的正確性。

TMM和DSM計算的鋼軌跨中豎向位移的結(jié)果和文獻(xiàn)中的鋼軌豎向位移結(jié)果見圖8。由圖8可知，兩種方法計算所得鋼軌位移時程曲線幾乎完全一致，且與文獻(xiàn)[30]的結(jié)果在幅值和趨勢上吻合良好，表明本文方法具有較高的計算精度，以及其用于軌道結(jié)構(gòu)的動力分析是可行的。

圖8 有砟軌道跨中鋼軌豎向位移

3.2 無砟軌道動力響應(yīng)分析

為進(jìn)一步對比TMM和DSM的計算效率，本節(jié)以移動荷載作用下的空間無砟軌道為例，通過不同方法計算耗時的對比，對TMM的計算效率進(jìn)行說明。

CRTSⅡ型板式無砟軌道長度為30 m，扣件縱向間距為0.6 m，軌道板長6 m。對于空間的無砟軌道，其元胞劃分方案與平面無砟軌道結(jié)構(gòu)類似，元胞A、B、C的長度均為0.1 m。相對于平面軌道結(jié)構(gòu)，空間無砟軌道的構(gòu)成較為復(fù)雜，其組成構(gòu)件中包含軌道板和底座板，因此難以通過理論方式準(zhǔn)確地描述其動力特性。本文采用有限元法建立元胞結(jié)構(gòu)的模型，并從中導(dǎo)出其質(zhì)量、剛度和阻尼矩陣進(jìn)行計算。無砟軌道結(jié)構(gòu)詳細(xì)參數(shù)見表1。

表1 CRTSⅡ型板式無砟軌道參數(shù)表

根據(jù)單節(jié)CRH2動車的軸荷載分布及軸重，構(gòu)造如圖9所示移動荷載。假定荷載以速度300 km/h沿鋼軌從左至右移動，數(shù)值積分步長取1/10 000 s，分別采用TMM和DSM分析了無砟軌道在移動荷載作用下的動力響應(yīng)。

圖9 CRH2動車輪重荷載配置(單位：m)

TMM和DSM計算的軌道跨中豎向動力響應(yīng)時程曲線見圖10，由圖10可知，兩種方法計算所得的軌道、軌道板和底座板動力響應(yīng)時程曲線均吻合良好。其中，跨中處鋼軌豎向位移、鋼軌豎向加速度、軌道板豎向位移和底座板豎向位移的最大值分別為0.681、79.561、0.221、0.131 mm。

圖10 CRTSⅡ型板式無砟軌道跨中豎向響應(yīng)

為說明本文方法具有更高的計算效率，對比了TMM和DSM兩種方法的計算耗時，見表2。由表2可知，采用TMM時，由于無需建立整體模型，僅需對元胞結(jié)構(gòu)進(jìn)行計算，極大地減少了軌道自由度數(shù)目，降低了計算過程中涉及矩陣的階數(shù)，因此計算內(nèi)存占用小，計算耗時較短，相對于DSM計算耗時降低了83.1%。

表2 兩種方法計算耗時

4 結(jié)論

本文針對軌道結(jié)構(gòu)的動力分析問題，基于傳遞矩陣法建立周期性軌道結(jié)構(gòu)模型，并在時域范圍內(nèi)求解其動力響應(yīng)。分別針對有砟和無砟軌道結(jié)構(gòu)，提出了對應(yīng)的元胞劃分方案，建立了二者的傳遞模型，并分別以移動荷載通過有砟軌道和無砟軌道為例，驗證了該方法的正確性和高效性。主要得到了以下結(jié)論：

(1) 基于軌道結(jié)構(gòu)的周期特性，本文方法僅需建立有限類元胞結(jié)構(gòu)模型，即可實現(xiàn)軌道結(jié)構(gòu)的整體求解，建模便捷且有效降低了結(jié)構(gòu)的自由度數(shù)目，及計算過程中的內(nèi)存占用。

(2) 本文方法僅需建立元胞結(jié)構(gòu)動力方程，而無需建立軌道整體動力方程，有效降低了計算過程中涉及矩陣的階數(shù)，避免了大型方程組的求解。

(3) 針對有砟和無砟軌道，分別提出了對應(yīng)的元胞劃分方案；對比直接剛度法計算結(jié)果，表明該方法具有較高的計算精度和計算效率，并且可用于多種軌道結(jié)構(gòu)的動力響應(yīng)分析。