車輛多體系統振動方程建立探討

全順喜，王 平，趙才友

(1.中鐵第四勘察設計院集團有限公司，武漢 430063;2.西南交通大學 高速鐵路線路工程教育部重點實驗室，成都 610031)

進行車輛-軌道耦合動力學研究時，需先建立車輛-軌道耦合系統振動方程。通常將車輛模擬成由車體、轉向架、輪對組成的多剛體子系統，將軌道模擬成由梁、質量塊等組成的子系統，分別建立振動方程，并通過復雜輪軌接觸實現兩子系統動態耦合［1］。建立車輛多體系統振動方程常用方法有三種:①利用達朗貝爾原理直接平衡法［2］，該法對簡單問題較直接、方便，但對自由度多的車輛系統須先對車輛各部件受力詳細分析后依次寫出各剛體運動方程。②矩陣組裝法［3-4］，該法雖可利用計算機自動形成車輛振動方程，但須先手工推導出各類變換矩陣。③能量法，如哈密爾頓原理、拉格朗日方程、勢能駐值原理等［5-8］，該法從功、能角度建立系統振動方程，較適用于難直接寫出力平衡方程的復雜結構，而哈密爾頓原理因只含純粹標量，較方便應用。

用哈密爾頓與勢能駐值原理時，需用“對號入座”法則建立系統振動方程［8-10］。但目前該法則尚未用有限元法中計算機編碼法，不能將單元矩陣自動組建成總體矩陣。為此，文獻［11］采用有限單元法建立車輛多體系統振動方程。而本文對“對號入座”法則進行改進，使之與有限元分析中計算機編碼法統一，快速利用計算機將單元矩陣自動建成系統總體矩陣，方便了應用哈密爾頓原理建立車輛多體系統振動方程過程。

1 哈密爾頓原理與改進的“對號入座”法則

1.1 哈密爾頓原理簡介［8］

據哈密爾頓原理提供的準則，利用動能、勢能等標量可將實際發生的運動從可能發生的(約束所許可的)運動中挑選出來。哈密爾頓原理表達式為:

式中:T為系統總動能;U為系統總勢能;δ為變分或虛位移符號;δW為保守力及非保守力所做虛功總和;t1，t2為積分起止時間。式(1)表示在任何時間段內，系統總動能與總勢能之差的變分加保守力與非保守力所做功的變分等于零。應用此原理，可方便地據系統總動能、總勢能及相應的虛功直接導出運動方程。

1.2 應用哈密爾頓原理建立振動方程

以圖1振動系統為例。其中m1～m3分別為三參振體質量;k1～k3分別為三彈簧剛度值;c1～c3分別為三阻尼器阻尼值;z1～z3分別為三參振動體振動位移，向下為正;P(t)為作用于系統隨時間變化外荷載。三自由度振動系統總動能一階變分表達式為:

圖1 振動系統Fig.1 Vibration system

由式(1)關于系統總動能積分形式，利用逐步積分得:

據“對號入座”法則，將與加速度對應項放入總質量矩陣中，其對應系統質量矩陣［M］可由式(3)加速度與變分角標獲得，即:

系統總勢能一階變分式為:

為利用“對號入座”法則，將式(5)改寫成位移與變分相乘形式，即:

將與位移對應項放入總剛度矩陣中，其對應系統剛度矩陣［K］可由式(6)中位移與變分角標獲得，即:

同理可得系統阻尼力所做虛功，利用“對號入座”法則，將之改寫成速度與變分相乘形式，即:

將與速度對應項放入總阻尼矩陣，其對應系統阻尼矩陣［C］可由式(8)中速度與變分角標獲得，即:

外荷載所做虛功δW2=δz1P(t)對應系統的外荷載列陣［Q］，即:

至此，系統總動能、總勢能變分及阻尼力、外力所做虛功均已求出，令 {q}=［z1z2z3］T、{δq}=［δz1δz2δz3］，利用哈密爾頓原理得:

由于{δq}的任意性，須式(11)括號內表達式為零時方可滿足，由此得系統振動方程為:

1.3 “對號入座”法則改進

由以上推導知，在利用“對號入座”法則獲取系統質量、剛度、阻尼及荷載矩陣時，需先將系統總動能、總勢能、阻尼及外力所做虛功寫成加速度與變分、位移與變分、速度與變分、外力與變分相乘形式，再據其角標將各項分別放入對應的質量、剛度、阻尼、荷載矩陣中，而未利用有限元法中用計算機將單元矩陣自動組成總體矩陣優勢。系統自由度較多時，此過程仍較繁瑣。為此，本文對其進行改進，將式(3)、式(5)、式(8)及外力所做虛功寫成矩陣形式，再將各參振質量塊、彈簧、阻尼及外力對應的矩陣分別合并獲取系統總質量、總剛度、總阻尼、總荷載矩陣。以形成系統總剛度矩陣為例，將式(5)改寫成矩陣形式為:

式中:K1，K2，K3分別為彈簧k1，k2，k3對應的剛度矩陣，由式(5)直接得出為:k1［-1，1］T［-1，1］、k2［-1，1］T［-1，1］、k3［0，1］T［0，1］。將彈簧k1，k2，k3對應的剛度矩陣合并即得系統總剛度矩陣，即:

系統總質量、總阻尼、總荷載矩陣均可按該方法獲取。由此，① 該方法中各參振質量塊、彈簧、阻尼、外力對應的子矩陣較易寫出，且具有類似性，方便計算機輸入;② 系統自由度較多時，通過改進后方法能借鑒有限元分析思想，可據編好的各參振質量塊、彈簧、阻尼、外力對應的子矩陣在總質量、總剛度、總阻尼矩陣中序號，利用計算機將子矩陣自動組成系統總體矩陣。極大方便了自由度多的車輛多體系統振動方程建立。

2 建立全車-軌道垂向耦合振動方程

2.1 全車-軌道垂向振動模型

為建立全車-軌道垂向耦合振動方程，先建立該系統振動模型，見圖2。該模型中車輛子系統由一個車體、兩個轉向架、四個輪對共七個剛體及一、二系懸掛組成。對車體及轉向架考慮沉浮、側滾、點頭3個自由度，對輪對考慮沉浮、側滾2個自由度，整個車輛子系統共17個自由度。軌道子系統包括以歐拉梁模擬的鋼軌、軌枕、以彈簧阻尼裝置模擬的扣件、道床，且鋼軌采用彈性點支承模型，軌枕采用連續支承模型，并用有限單元法進行離散。

圖2 全車-軌道垂向振動模型Fig.2 Vertical model for vehicle-track coupling system

2.2 全車-軌道垂向耦合振動方程建立

由圖2容易獲得車輛子系統總動能，考慮到哈密爾頓原理中關于系統總動能積分形式，將其逐步積分得:

車輛子系統總勢能為:

令二系懸掛連接的車體、轉向架位移{qcti}=［zc，φc，βc，zti，φti，βti］T、變分{δqcti}=［δzc，δφc，δβc，δzti，δφti，δβti］;一系懸掛連接的轉向架、輪對位移{qtwim}=［zti，φti，βti，zwim，φwim］、變分{δqtwim}= ［δzti，δφti，δβti，δzwim，δφwim］。利用改進后的“對號入座”法則，式(16)一階變分寫為:

式(15)～式(17)中:i表示車體下有兩個轉向架;m表示轉向架下有兩個輪對;j=1表示左側懸掛，j=2表示右側懸掛;mc，Icx，Icy為車體質量、側滾、點頭轉動慣量;mt，Itx，Ity為轉向架質量、側滾、點頭轉動慣量;mw，Iwx為輪對質量、側滾轉動慣量;zc，φc，βc表示車體沉浮、側股、點頭;zti，φti，βti表示轉向架沉浮、側股、點頭;zwim，φwim表示輪對沉浮、側股;Lc，Lt，d1，d2表示一、二系懸掛縱向距離及橫向距離之半;k1z，k2z，{K1zimj}，{K2zij}分別為一、二系懸掛垂向剛度及其剛度子矩陣，由式(16)可得{K1zimj}=k1z［1，(-1)jd1，(-1)mLt， -1，(-1)j+1d1］T［1，(-1)jd1，(-1)mLt，-1，(-1)j+1d1］，{K2zij}=k2z［1，(-1)jd2，(-1)iLc，-1，(-1)j+1d2，0］T［1，(-1)jd2，(-1)iLc，-1，(-1)j+1d2，0］。

車輛子系統中阻尼力所做虛功為:

為利用改進后“對號入座”，將式(18)改寫為:

式中:c1z，c2z，{C1zimj}，{C2zij}分別為一、二系懸掛的垂向阻尼及其阻尼子矩陣，由式(18)可得{C1zimj}={K1zimj}c1z/k1z，{C2zij}={K2zij}c2z/k2z。

輪軌垂向力對車輛子系統所做虛功為:

式中:Pimj為輪對左右側輪軌垂向力，向下為正。同理，也可將式(20)寫成荷載子列陣形式。

由式(15)、式(17)、式(19)、式(20)，據哈密爾頓原理，用以上組建總質量、總剛度、總阻尼、總荷載矩陣方法，利用計算機自動組建車輛系統振動方程，即:

式中:［M］c，［C］c，［K］c，{Pc}為車輛子系統總質量、總阻尼、總剛度、總荷載矩陣;{uc}為車輛系統位移矩陣。對軌道子系統，可用有限單元法建立振動方程［8］。

2.3 振動方程求解

車輛、軌道子系統是相互作用、密不可分的兩個部分，主要通過輪軌間作用力相互聯系。對求輪軌間作用力而言，文獻［1］采用新型積分方法［12-13］求解系統振動方程，該方法中下一時刻位移列陣只與本時刻位移、速度、加速度列陣及上一時刻加速度列陣有關，下一時刻輪軌作用力可據本時刻及上一時刻振動響應求出，勿需迭代求解;文獻［8-14］采用Park法求解系統振動方程組，而據Park法，下一時刻位移列陣不但與本時刻及上一時刻位移、速度、加速度列陣有關，且與下一時刻荷載列向量即輪軌作用力有關，因此，下一時刻輪軌作用力需反復迭代，對輪軌間垂向力，可采用牛頓迭代法求解。

由于本文用有限單元法離散軌道子系統，其質量矩陣非對角陣，即不能用新型積分方法求解振動方程。故采用Park法求解全車-軌道垂向耦合振動方程?？紤]到用牛頓迭代法求解輪軌垂向力時要寫出迭代雅克比陣，計算較復雜，且在車輛-軌道空間耦合振動分析中，輪軌間法向力不但與輪對垂向位移、輪對側滾角、鋼軌垂向位移、垂向不平順有關，且與輪對橫向位移、輪對搖頭角、鋼軌橫向位移、車輪踏面、輪下鋼軌輪廓等有關，該迭代雅克比陣無法寫出。因此先將非線性赫茲彈簧簡化成線性彈簧，采用直接迭代法求解輪軌垂向力。

設車輛靜輪重P0，輪軌間接觸等效線性剛度為［15］:

式中:G為輪軌接觸常數，踏面為錐形時，G=4.57R-0.149×10-8;踏面為磨耗型時，G=4.57R-0.149×10-8(R為車輪半徑)。數值試驗表明，積分步長小于0.2 ms時，迭代過程收斂，約迭代5次即可滿足1%精度要求。

據全車-軌道垂向振動模型及求解方法，用MATLAB編制的計算程序，可用于各種條件下全車-軌道垂向耦合振動分析。

3 焊接凹接頭激勵下車輛-軌道垂向振動特性

3.1 計算參數

為驗證建立車輛多體系統振動方程方法的正確性，對焊接凹接頭激勵下車輛-軌道垂向振動特性進行分析。采用CRH3型動車組、CHN60鋼軌，彈性模量E=2.1 ×105MPa，熱膨脹系數 α =1.18 ×10-5，泊松比0.3，扣件間距0.625 m。軌道子系統相關參數見表1。

鋼軌焊接凹接頭可用兩波長與幅值不同的余弦波疊加表示，取兩余弦波波長分別為1 m、0.1 m，幅值分別為0.2 mm，0.1 mm，如圖3 所示［16］。

圖3 鋼軌焊接凹接頭不平順Fig.3 Track irregularities of welded concave joints

3.2 計算結果

車輛以250 km/h通過焊接凹接頭時，車輛-軌道垂向振動響應如圖4所示。

由圖4看出，焊接凹接頭對車輛-軌道垂向振動特性影響較大。車輛進入焊接凹接頭區域時，凹接頭側輪軌垂向力迅速減載至28.5 kN，隨后又增大到117.3 kN，約為靜輪載的1.7倍，而非凹接頭側輪軌垂向力波動較小;車體垂向加速度存在兩個明顯波動區，由車輛前后轉向架通過焊接凹接頭時所致，但由于車輛懸掛系統隔振效果，波動幅值較小，最大值不超過0.04 m/s2;車輛通過時，焊接凹接頭附近鋼軌、軌枕加速度、位移較大，加速度最大值分別為1 496 m/s2、127 m/s2，位移最大值分別為 0.77 mm、0.23 mm，比無不平順時增加許多，會影響鋼軌使用壽命。本文計算結果在波形上或量值上均與文獻［16］結果接近，說明本文所提建立車輛多體系統振動方程方法正確、可行。

圖4 焊接凹接頭不平順下系統振動響應Fig.4 Vibration responses of the vehicle-track coupling system under the welded concave joints irregularity

4 結論

為方便應用哈密爾頓原理建立車輛多體系統振動方程，本文工作為:

(1)以三自由度振動系統為例，對哈密爾頓原理“對號入座”進行改進，使之與有限元分析中計算機編碼法統一，快速利用計算機將單元矩陣自動組成系統總體矩陣。

(2)建立全車-軌道垂向振動模型，應用改進后“對號入座”法則建立該振動系統方程，并應用MATLAB編制相應計算程序。

(3)用所編程序對焊接凹接頭激勵下車輛-軌道垂向振動特性進行分析，結果表明，本文計算結果與之前研究結論較一致，表明本文所提方法的正確性。

［1］翟婉明.車輛-軌道耦合動力學(第三版)［M］.北京:科學出版社，2007.

［2］翟婉明.車輛-軌道垂向系統的統一模型及其耦合動力學原理［J］.鐵道學報，1992，14(3):10-21.

ZHAI Wan-ming.The vertical model of vehicle-track system and its coupling dynamics［J］.Journal of The China Railway Society，1992，14(3):10-21.

［3］張定賢.機車車輛軌道系統動力學［M］.北京:中國鐵道出版社，1996.

［4］倪純雙，洪嘉振，賀啟庸.車輛系統動力學的符號化多體方法［J］.中國鐵道科學，1996，17(4):12-17.

NI Chun-shuang，HONG Jia-zhen，HE Qi-yong.Symbolic multibody approach of vehicle system dynamics［J］.China Railway Science，1996，17(4):12-17.

［5］曾慶元.彈性系統動力學總勢能不變值原理［J］.華中理工大學學報，2000，28(1):1-3.

ZENG Qing-yuan.The principle of total potential energy with stationary value in elastic system dynamics［J］.Journal of Huazhong University of Science and Technology，2000，28(1):1-3.

［6］曾慶元.彈性系統動力學總勢能不變值原理與列車橋梁時變系統振動分析［J］.鐵道建筑技術，2001(1):1-6.

ZENG Qing-yuan.The principle of total potential energy with stationary value in elastic system dynamics and the vibration analysis of hain-bridge time-varying system［J］.Railway Construction Technology，2001(1):1-6.

［7］張亞輝，林家浩.結構動力學基礎［M］.大連:大連理工大學出版社，2007.

［8］劉學毅，王 平.車輛-軌道-路基系統動力學［M］.成都:西南交通大學出版社，2010.

［9］曾慶元，楊 平.形成矩陣的“對號入座”法則與桁梁空間分析的桁段有限元法［J］.鐵道學報，1986，8(2):48-59.

ZENG Qing-yuan，YANG Ping.The"set-in-right-position"rule for forming structural matrices and the finite trusselement method for space analysis of truss bridges［J］.Journal of The China Railway Society，1986，8(2):48-59.

［10］向 俊，曾慶元.關于機車車輛-軌道系統運動方程的建立［J］.長沙鐵道學院學報，2000，18(4):1-5.

XIANG Jun，ZENG Qing-yuan.On building equations of motion of locomotive vehicle-track system［J］.Journal of Changsha Railway University，2000，18(4):1-5.

［11］李東平，曾慶元，婁 平.車輛多體系統動力學方程的有限元法［J］.中國鐵道科學，2004，25(5):33-38.

LI Dong-ping，ZENG Qing-yuan，LOU Ping.The finite element method forestablishingdynamicalequationsof vehicle multi-body system［J］.China Railway Science，2004，25(5):33-38.

［12］翟婉明.非線性結構動力分析的Newmark預測-校正積分模式［J］.計算結構力學及其應用，1990，7(2):51-58.

ZHAI Wan-ming.The predictor-corrector scheme based on the newmarks integration algorithm for nonlinear structural dynamic analyses［J］.Chinese Journal of Computational Mechanics，1990，7(2):51-58.

［13］Zhai W M.Two simple fast integration methods for large-scale dynamic problems in engineering［J］.International Journal for Numerical Methods in Engineering，1996，39(24):4199-4214.

［14］王 平.道岔區輪軌系統動力學的研究［D］.成都:西南交通大學，1998.

［15］李成輝，萬復光.軌道結構隨機振動［J］.鐵道學報，1998(3):97-101.

LI Cheng-hui，WAN Fu-guang.Random vibration of track structure［J］. Journal of The China Railway Society，1998(3):97-101.

［16］王偉華.土路基上雙塊式無砟軌道垂向動力特性分析［D］.成都:西南交通大學，2009.