基于傳遞矩陣法的軌道結構動力分析方法

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

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

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

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

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

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

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

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

1 時域傳遞矩陣法

1.1 元胞單元時域傳遞關系

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

圖1 周期性結構的元胞結構劃分示意

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

(1)

采用Newmark-β法對上述方程進行求解，則有

(2)

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

則對于元胞i有

(3)

將其代入式(2)可得

(4)

(5)

展開式(5)，則有

(6)

(7)

考慮位移和內力的連續性[25]，則剖面兩側的位移相等，內力相反，故相鄰兩元胞間狀態矢量的傳遞關系可表示為

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

(8)

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

Ni,I=SiUi,I+Eii>2

(9)

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

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

(10)

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

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

(11)

(12)

(13)

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

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

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

(14)

式中：n為系統元胞個數。

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

1.2 傳遞系數求解

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

N2,I=S2U2,I+E2

(15)

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

N1,O=-S2U1,O-E2

(16)

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

(17)

(18)

考慮結構輸入端邊界條件，結合式(17)、式(18)，可解得U1,O和N1,O間的關系，進而對照式(15)可解得對應邊界條件下的S2和E2。

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

(19)

對照式(16)可知

(20)

(21)

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

(22)

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

(23)

對照式(16)可知

(24)

(25)

(26)

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

(27)

(28)

1.3 計算流程

根據上述理論，基于Matlab平臺編制了外荷載作用下的周期結構瞬態響應的TMM求解程序，主要計算流程見圖3。

圖3 傳遞矩陣法計算流程

2 軌道結構的傳遞矩陣模型

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

2.1 有砟軌道元胞模型及傳遞關系

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

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

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

圖5 有砟軌道結構劃分示意

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

(29)

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

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

(30)

式中：Z和θ分別為位移和轉角；Q和M分別為剪力和彎矩；下標r和s分別為鋼軌和軌枕。

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

(31)

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

(32)

(33)

軌道結構兩端無內力作用，適用于荷載邊界條件，式(32)和式(33)中Si和Ei可在由式(24)和式(25)解得傳遞系數S2和E2后，循環應用式(12)和式(13)計算得到。

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

(34)

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

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

2.2 無砟軌道元胞模型及傳遞關系

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

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

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

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

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

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

(35)

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

3 移動荷載作用下軌道結構動力響應分析

為驗證上述關于軌道結構傳遞矩陣算法的正確性，本節以移動荷載作用下的有砟軌道結構為例，對比了TMM、DSM和文獻中的鋼軌動力響應計算結果。并以自由度數目較多的空間無砟軌道結構為例，對比了TMM和DSM計算所得的移動荷載作用下軌道結構動力響應結果，進一步驗證了本文方法的高效性。

3.1 有砟軌道動力響應分析

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

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

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

3.2 無砟軌道動力響應分析

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

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

表1 CRTSⅡ型板式無砟軌道參數表

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

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

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

圖10 CRTSⅡ型板式無砟軌道跨中豎向響應

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

表2 兩種方法計算耗時

4 結論

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

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

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

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