全順喜，王 平，趙才友

(1.中鐵第四勘察設(shè)計(jì)院集團(tuán)有限公司，武漢 430063;2.西南交通大學(xué) 高速鐵路線路工程教育部重點(diǎn)實(shí)驗(yàn)室，成都 610031)

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

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

1 哈密爾頓原理與改進(jìn)的“對號入座”法則

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

據(jù)哈密爾頓原理提供的準(zhǔn)則，利用動(dòng)能、勢能等標(biāo)量可將實(shí)際發(fā)生的運(yùn)動(dòng)從可能發(fā)生的(約束所許可的)運(yùn)動(dòng)中挑選出來。哈密爾頓原理表達(dá)式為:

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

1.2 應(yīng)用哈密爾頓原理建立振動(dòng)方程

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

圖1 振動(dòng)系統(tǒng)Fig.1 Vibration system

由式(1)關(guān)于系統(tǒng)總動(dòng)能積分形式，利用逐步積分得:

據(jù)“對號入座”法則，將與加速度對應(yīng)項(xiàng)放入總質(zhì)量矩陣中，其對應(yīng)系統(tǒng)質(zhì)量矩陣［M］可由式(3)加速度與變分角標(biāo)獲得，即:

系統(tǒng)總勢能一階變分式為:

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

將與位移對應(yīng)項(xiàng)放入總剛度矩陣中，其對應(yīng)系統(tǒng)剛度矩陣［K］可由式(6)中位移與變分角標(biāo)獲得，即:

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

將與速度對應(yīng)項(xiàng)放入總阻尼矩陣，其對應(yīng)系統(tǒng)阻尼矩陣［C］可由式(8)中速度與變分角標(biāo)獲得，即:

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

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

由于{δq}的任意性，須式(11)括號內(nèi)表達(dá)式為零時(shí)方可滿足，由此得系統(tǒng)振動(dòng)方程為:

1.3 “對號入座”法則改進(jìn)

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

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

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

2 建立全車-軌道垂向耦合振動(dòng)方程

2.1 全車-軌道垂向振動(dòng)模型

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

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

2.2 全車-軌道垂向耦合振動(dòng)方程建立

由圖2容易獲得車輛子系統(tǒng)總動(dòng)能，考慮到哈密爾頓原理中關(guān)于系統(tǒng)總動(dòng)能積分形式，將其逐步積分得:

車輛子系統(tǒng)總勢能為:

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

式(15)～式(17)中:i表示車體下有兩個(gè)轉(zhuǎn)向架;m表示轉(zhuǎn)向架下有兩個(gè)輪對;j=1表示左側(cè)懸掛，j=2表示右側(cè)懸掛;mc，Icx，Icy為車體質(zhì)量、側(cè)滾、點(diǎn)頭轉(zhuǎn)動(dòng)慣量;mt，Itx，Ity為轉(zhuǎn)向架質(zhì)量、側(cè)滾、點(diǎn)頭轉(zhuǎn)動(dòng)慣量;mw，Iwx為輪對質(zhì)量、側(cè)滾轉(zhuǎn)動(dòng)慣量;zc，φc，βc表示車體沉浮、側(cè)股、點(diǎn)頭;zti，φti，βti表示轉(zhuǎn)向架沉浮、側(cè)股、點(diǎn)頭;zwim，φwim表示輪對沉浮、側(cè)股;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］。

車輛子系統(tǒng)中阻尼力所做虛功為:

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

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

輪軌垂向力對車輛子系統(tǒng)所做虛功為:

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

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

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

2.3 振動(dòng)方程求解

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

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

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

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

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

3 焊接凹接頭激勵(lì)下車輛-軌道垂向振動(dòng)特性

3.1 計(jì)算參數(shù)

為驗(yàn)證建立車輛多體系統(tǒng)振動(dòng)方程方法的正確性，對焊接凹接頭激勵(lì)下車輛-軌道垂向振動(dòng)特性進(jìn)行分析。采用CRH3型動(dòng)車組、CHN60鋼軌，彈性模量E=2.1 ×105MPa，熱膨脹系數(shù) α =1.18 ×10-5，泊松比0.3，扣件間距0.625 m。軌道子系統(tǒng)相關(guān)參數(shù)見表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 計(jì)算結(jié)果

車輛以250 km/h通過焊接凹接頭時(shí)，車輛-軌道垂向振動(dòng)響應(yīng)如圖4所示。

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

圖4 焊接凹接頭不平順下系統(tǒng)振動(dòng)響應(yīng)Fig.4 Vibration responses of the vehicle-track coupling system under the welded concave joints irregularity

4 結(jié)論

為方便應(yīng)用哈密爾頓原理建立車輛多體系統(tǒng)振動(dòng)方程，本文工作為:

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

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

(3)用所編程序?qū)附影冀宇^激勵(lì)下車輛-軌道垂向振動(dòng)特性進(jìn)行分析，結(jié)果表明，本文計(jì)算結(jié)果與之前研究結(jié)論較一致，表明本文所提方法的正確性。

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

［2］翟婉明.車輛-軌道垂向系統(tǒng)的統(tǒng)一模型及其耦合動(dòng)力學(xué)原理［J］.鐵道學(xué)報(bào)，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］張定賢.機(jī)車車輛軌道系統(tǒng)動(dòng)力學(xué)［M］.北京:中國鐵道出版社，1996.

［4］倪純雙，洪嘉振，賀啟庸.車輛系統(tǒng)動(dòng)力學(xué)的符號化多體方法［J］.中國鐵道科學(xué)，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］曾慶元.彈性系統(tǒng)動(dòng)力學(xué)總勢能不變值原理［J］.華中理工大學(xué)學(xué)報(bào)，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］曾慶元.彈性系統(tǒng)動(dòng)力學(xué)總勢能不變值原理與列車橋梁時(shí)變系統(tǒng)振動(dòng)分析［J］.鐵道建筑技術(shù)，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］張亞輝，林家浩.結(jié)構(gòu)動(dòng)力學(xué)基礎(chǔ)［M］.大連:大連理工大學(xué)出版社，2007.

［8］劉學(xué)毅，王 平.車輛-軌道-路基系統(tǒng)動(dòng)力學(xué)［M］.成都:西南交通大學(xué)出版社，2010.

［9］曾慶元，楊 平.形成矩陣的“對號入座”法則與桁梁空間分析的桁段有限元法［J］.鐵道學(xué)報(bào)，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］向 俊，曾慶元.關(guān)于機(jī)車車輛-軌道系統(tǒng)運(yùn)動(dòng)方程的建立［J］.長沙鐵道學(xué)院學(xué)報(bào)，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］李東平，曾慶元，婁 平.車輛多體系統(tǒng)動(dòng)力學(xué)方程的有限元法［J］.中國鐵道科學(xué)，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］翟婉明.非線性結(jié)構(gòu)動(dòng)力分析的Newmark預(yù)測-校正積分模式［J］.計(jì)算結(jié)構(gòu)力學(xué)及其應(yīng)用，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］王 平.道岔區(qū)輪軌系統(tǒng)動(dòng)力學(xué)的研究［D］.成都:西南交通大學(xué)，1998.

［15］李成輝，萬復(fù)光.軌道結(jié)構(gòu)隨機(jī)振動(dòng)［J］.鐵道學(xué)報(bào)，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òng)力特性分析［D］.成都:西南交通大學(xué)，2009.