全順喜,王 平,趙才友
(1.中鐵第四勘察設(shè)計(jì)院集團(tuán)有限公司,武漢 430063;2.西南交通大學(xué) 高速鐵路線路工程教育部重點(diǎn)實(shí)驗(yàn)室,成都 610031)
進(jìn)行車(chē)輛-軌道耦合動(dòng)力學(xué)研究時(shí),需先建立車(chē)輛-軌道耦合系統(tǒng)振動(dòng)方程。通常將車(chē)輛模擬成由車(chē)體、轉(zhuǎn)向架、輪對(duì)組成的多剛體子系統(tǒng),將軌道模擬成由梁、質(zhì)量塊等組成的子系統(tǒng),分別建立振動(dòng)方程,并通過(guò)復(fù)雜輪軌接觸實(shí)現(xiàn)兩子系統(tǒng)動(dòng)態(tài)耦合[1]。建立車(chē)輛多體系統(tǒng)振動(dòng)方程常用方法有三種:①利用達(dá)朗貝爾原理直接平衡法[2],該法對(duì)簡(jiǎn)單問(wèn)題較直接、方便,但對(duì)自由度多的車(chē)輛系統(tǒng)須先對(duì)車(chē)輛各部件受力詳細(xì)分析后依次寫(xiě)出各剛體運(yùn)動(dòng)方程。②矩陣組裝法[3-4],該法雖可利用計(jì)算機(jī)自動(dòng)形成車(chē)輛振動(dòng)方程,但須先手工推導(dǎo)出各類(lèi)變換矩陣。③能量法,如哈密爾頓原理、拉格朗日方程、勢(shì)能駐值原理等[5-8],該法從功、能角度建立系統(tǒng)振動(dòng)方程,較適用于難直接寫(xiě)出力平衡方程的復(fù)雜結(jié)構(gòu),而哈密爾頓原理因只含純粹標(biāo)量,較方便應(yīng)用。
用哈密爾頓與勢(shì)能駐值原理時(shí),需用“對(duì)號(hào)入座”法則建立系統(tǒng)振動(dòng)方程[8-10]。但目前該法則尚未用有限元法中計(jì)算機(jī)編碼法,不能將單元矩陣自動(dòng)組建成總體矩陣。為此,文獻(xiàn)[11]采用有限單元法建立車(chē)輛多體系統(tǒng)振動(dòng)方程。而本文對(duì)“對(duì)號(hào)入座”法則進(jìn)行改進(jìn),使之與有限元分析中計(jì)算機(jī)編碼法統(tǒng)一,快速利用計(jì)算機(jī)將單元矩陣自動(dòng)建成系統(tǒng)總體矩陣,方便了應(yīng)用哈密爾頓原理建立車(chē)輛多體系統(tǒng)振動(dòng)方程過(guò)程。
據(jù)哈密爾頓原理提供的準(zhǔn)則,利用動(dòng)能、勢(shì)能等標(biāo)量可將實(shí)際發(fā)生的運(yùn)動(dòng)從可能發(fā)生的(約束所許可的)運(yùn)動(dòng)中挑選出來(lái)。哈密爾頓原理表達(dá)式為:
式中:T為系統(tǒng)總動(dòng)能;U為系統(tǒng)總勢(shì)能;δ為變分或虛位移符號(hào);δW為保守力及非保守力所做虛功總和;t1,t2為積分起止時(shí)間。式(1)表示在任何時(shí)間段內(nèi),系統(tǒng)總動(dòng)能與總勢(shì)能之差的變分加保守力與非保守力所做功的變分等于零。應(yīng)用此原理,可方便地?fù)?jù)系統(tǒng)總動(dòng)能、總勢(shì)能及相應(yīng)的虛功直接導(dǎo)出運(yùn)動(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ù)“對(duì)號(hào)入座”法則,將與加速度對(duì)應(yīng)項(xiàng)放入總質(zhì)量矩陣中,其對(duì)應(yīng)系統(tǒng)質(zhì)量矩陣[M]可由式(3)加速度與變分角標(biāo)獲得,即:
系統(tǒng)總勢(shì)能一階變分式為:
為利用“對(duì)號(hào)入座”法則,將式(5)改寫(xiě)成位移與變分相乘形式,即:
將與位移對(duì)應(yīng)項(xiàng)放入總剛度矩陣中,其對(duì)應(yīng)系統(tǒng)剛度矩陣[K]可由式(6)中位移與變分角標(biāo)獲得,即:
同理可得系統(tǒng)阻尼力所做虛功,利用“對(duì)號(hào)入座”法則,將之改寫(xiě)成速度與變分相乘形式,即:
將與速度對(duì)應(yīng)項(xiàng)放入總阻尼矩陣,其對(duì)應(yīng)系統(tǒng)阻尼矩陣[C]可由式(8)中速度與變分角標(biāo)獲得,即:
外荷載所做虛功δW2=δz1P(t)對(duì)應(yīng)系統(tǒng)的外荷載列陣[Q],即:
至此,系統(tǒng)總動(dòng)能、總勢(shì)能變分及阻尼力、外力所做虛功均已求出,令 {q}=[z1z2z3]T、{δq}=[δz1δz2δz3],利用哈密爾頓原理得:
由于{δq}的任意性,須式(11)括號(hào)內(nèi)表達(dá)式為零時(shí)方可滿(mǎn)足,由此得系統(tǒng)振動(dòng)方程為:
由以上推導(dǎo)知,在利用“對(duì)號(hào)入座”法則獲取系統(tǒng)質(zhì)量、剛度、阻尼及荷載矩陣時(shí),需先將系統(tǒng)總動(dòng)能、總勢(shì)能、阻尼及外力所做虛功寫(xiě)成加速度與變分、位移與變分、速度與變分、外力與變分相乘形式,再據(jù)其角標(biāo)將各項(xiàng)分別放入對(duì)應(yīng)的質(zhì)量、剛度、阻尼、荷載矩陣中,而未利用有限元法中用計(jì)算機(jī)將單元矩陣自動(dòng)組成總體矩陣優(yōu)勢(shì)。系統(tǒng)自由度較多時(shí),此過(guò)程仍較繁瑣。為此,本文對(duì)其進(jìn)行改進(jìn),將式(3)、式(5)、式(8)及外力所做虛功寫(xiě)成矩陣形式,再將各參振質(zhì)量塊、彈簧、阻尼及外力對(duì)應(yīng)的矩陣分別合并獲取系統(tǒng)總質(zhì)量、總剛度、總阻尼、總荷載矩陣。以形成系統(tǒng)總剛度矩陣為例,將式(5)改寫(xiě)成矩陣形式為:
式中:K1,K2,K3分別為彈簧k1,k2,k3對(duì)應(yīng)的剛度矩陣,由式(5)直接得出為:k1[-1,1]T[-1,1]、k2[-1,1]T[-1,1]、k3[0,1]T[0,1]。將彈簧k1,k2,k3對(duì)應(yīng)的剛度矩陣合并即得系統(tǒng)總剛度矩陣,即:
系統(tǒng)總質(zhì)量、總阻尼、總荷載矩陣均可按該方法獲取。由此,① 該方法中各參振質(zhì)量塊、彈簧、阻尼、外力對(duì)應(yīng)的子矩陣較易寫(xiě)出,且具有類(lèi)似性,方便計(jì)算機(jī)輸入;② 系統(tǒng)自由度較多時(shí),通過(guò)改進(jìn)后方法能借鑒有限元分析思想,可據(jù)編好的各參振質(zhì)量塊、彈簧、阻尼、外力對(duì)應(yīng)的子矩陣在總質(zhì)量、總剛度、總阻尼矩陣中序號(hào),利用計(jì)算機(jī)將子矩陣自動(dòng)組成系統(tǒng)總體矩陣。極大方便了自由度多的車(chē)輛多體系統(tǒng)振動(dòng)方程建立。
為建立全車(chē)-軌道垂向耦合振動(dòng)方程,先建立該系統(tǒng)振動(dòng)模型,見(jiàn)圖2。該模型中車(chē)輛子系統(tǒng)由一個(gè)車(chē)體、兩個(gè)轉(zhuǎn)向架、四個(gè)輪對(duì)共七個(gè)剛體及一、二系懸掛組成。對(duì)車(chē)體及轉(zhuǎn)向架考慮沉浮、側(cè)滾、點(diǎn)頭3個(gè)自由度,對(duì)輪對(duì)考慮沉浮、側(cè)滾2個(gè)自由度,整個(gè)車(chē)輛子系統(tǒng)共17個(gè)自由度。軌道子系統(tǒng)包括以歐拉梁模擬的鋼軌、軌枕、以彈簧阻尼裝置模擬的扣件、道床,且鋼軌采用彈性點(diǎn)支承模型,軌枕采用連續(xù)支承模型,并用有限單元法進(jìn)行離散。
圖2 全車(chē)-軌道垂向振動(dòng)模型Fig.2 Vertical model for vehicle-track coupling system
由圖2容易獲得車(chē)輛子系統(tǒng)總動(dòng)能,考慮到哈密爾頓原理中關(guān)于系統(tǒng)總動(dòng)能積分形式,將其逐步積分得:
車(chē)輛子系統(tǒng)總勢(shì)能為:
令二系懸掛連接的車(chē)體、轉(zhuǎn)向架位移{qcti}=[zc,φc,βc,zti,φti,βti]T、變分{δqcti}=[δzc,δφc,δβc,δzti,δφti,δβti];一系懸掛連接的轉(zhuǎn)向架、輪對(duì)位移{qtwim}=[zti,φti,βti,zwim,φwim]、變分{δqtwim}= [δzti,δφti,δβti,δzwim,δφwim]。利用改進(jìn)后的“對(duì)號(hào)入座”法則,式(16)一階變分寫(xiě)為:
式(15)~式(17)中:i表示車(chē)體下有兩個(gè)轉(zhuǎn)向架;m表示轉(zhuǎn)向架下有兩個(gè)輪對(duì);j=1表示左側(cè)懸掛,j=2表示右側(cè)懸掛;mc,Icx,Icy為車(chē)體質(zhì)量、側(cè)滾、點(diǎn)頭轉(zhuǎn)動(dòng)慣量;mt,Itx,Ity為轉(zhuǎn)向架質(zhì)量、側(cè)滾、點(diǎn)頭轉(zhuǎn)動(dòng)慣量;mw,Iwx為輪對(duì)質(zhì)量、側(cè)滾轉(zhuǎn)動(dòng)慣量;zc,φc,βc表示車(chē)體沉浮、側(cè)股、點(diǎn)頭;zti,φti,βti表示轉(zhuǎn)向架沉浮、側(cè)股、點(diǎn)頭;zwim,φwim表示輪對(duì)沉浮、側(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]。
車(chē)輛子系統(tǒng)中阻尼力所做虛功為:
為利用改進(jìn)后“對(duì)號(hào)入座”,將式(18)改寫(xiě)為:
式中:c1z,c2z,{C1zimj},{C2zij}分別為一、二系懸掛的垂向阻尼及其阻尼子矩陣,由式(18)可得{C1zimj}={K1zimj}c1z/k1z,{C2zij}={K2zij}c2z/k2z。
輪軌垂向力對(duì)車(chē)輛子系統(tǒng)所做虛功為:
式中:Pimj為輪對(duì)左右側(cè)輪軌垂向力,向下為正。同理,也可將式(20)寫(xiě)成荷載子列陣形式。
由式(15)、式(17)、式(19)、式(20),據(jù)哈密爾頓原理,用以上組建總質(zhì)量、總剛度、總阻尼、總荷載矩陣方法,利用計(jì)算機(jī)自動(dòng)組建車(chē)輛系統(tǒng)振動(dòng)方程,即:
式中:[M]c,[C]c,[K]c,{Pc}為車(chē)輛子系統(tǒng)總質(zhì)量、總阻尼、總剛度、總荷載矩陣;{uc}為車(chē)輛系統(tǒng)位移矩陣。對(duì)軌道子系統(tǒng),可用有限單元法建立振動(dòng)方程[8]。
車(chē)輛、軌道子系統(tǒng)是相互作用、密不可分的兩個(gè)部分,主要通過(guò)輪軌間作用力相互聯(lián)系。對(duì)求輪軌間作用力而言,文獻(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ù)迭代,對(duì)輪軌間垂向力,可采用牛頓迭代法求解。
由于本文用有限單元法離散軌道子系統(tǒng),其質(zhì)量矩陣非對(duì)角陣,即不能用新型積分方法求解振動(dòng)方程。故采用Park法求解全車(chē)-軌道垂向耦合振動(dòng)方程??紤]到用牛頓迭代法求解輪軌垂向力時(shí)要寫(xiě)出迭代雅克比陣,計(jì)算較復(fù)雜,且在車(chē)輛-軌道空間耦合振動(dòng)分析中,輪軌間法向力不但與輪對(duì)垂向位移、輪對(duì)側(cè)滾角、鋼軌垂向位移、垂向不平順有關(guān),且與輪對(duì)橫向位移、輪對(duì)搖頭角、鋼軌橫向位移、車(chē)輪踏面、輪下鋼軌輪廓等有關(guān),該迭代雅克比陣無(wú)法寫(xiě)出。因此先將非線性赫茲彈簧簡(jiǎn)化成線性彈簧,采用直接迭代法求解輪軌垂向力。
設(shè)車(chē)輛靜輪重P0,輪軌間接觸等效線性剛度為[15]:
式中:G為輪軌接觸常數(shù),踏面為錐形時(shí),G=4.57R-0.149×10-8;踏面為磨耗型時(shí),G=4.57R-0.149×10-8(R為車(chē)輪半徑)。數(shù)值試驗(yàn)表明,積分步長(zhǎng)小于0.2 ms時(shí),迭代過(guò)程收斂,約迭代5次即可滿(mǎn)足1%精度要求。
據(jù)全車(chē)-軌道垂向振動(dòng)模型及求解方法,用MATLAB編制的計(jì)算程序,可用于各種條件下全車(chē)-軌道垂向耦合振動(dòng)分析。
為驗(yàn)證建立車(chē)輛多體系統(tǒng)振動(dòng)方程方法的正確性,對(duì)焊接凹接頭激勵(lì)下車(chē)輛-軌道垂向振動(dòng)特性進(jìn)行分析。采用CRH3型動(dòng)車(chē)組、CHN60鋼軌,彈性模量E=2.1 ×105MPa,熱膨脹系數(shù) α =1.18 ×10-5,泊松比0.3,扣件間距0.625 m。軌道子系統(tǒng)相關(guān)參數(shù)見(jiàn)表1。
鋼軌焊接凹接頭可用兩波長(zhǎng)與幅值不同的余弦波疊加表示,取兩余弦波波長(zhǎng)分別為1 m、0.1 m,幅值分別為0.2 mm,0.1 mm,如圖3 所示[16]。
圖3 鋼軌焊接凹接頭不平順Fig.3 Track irregularities of welded concave joints
車(chē)輛以250 km/h通過(guò)焊接凹接頭時(shí),車(chē)輛-軌道垂向振動(dòng)響應(yīng)如圖4所示。
由圖4看出,焊接凹接頭對(duì)車(chē)輛-軌道垂向振動(dòng)特性影響較大。車(chē)輛進(jìn)入焊接凹接頭區(qū)域時(shí),凹接頭側(cè)輪軌垂向力迅速減載至28.5 kN,隨后又增大到117.3 kN,約為靜輪載的1.7倍,而非凹接頭側(cè)輪軌垂向力波動(dòng)較小;車(chē)體垂向加速度存在兩個(gè)明顯波動(dòng)區(qū),由車(chē)輛前后轉(zhuǎn)向架通過(guò)焊接凹接頭時(shí)所致,但由于車(chē)輛懸掛系統(tǒng)隔振效果,波動(dòng)幅值較小,最大值不超過(guò)0.04 m/s2;車(chē)輛通過(guò)時(shí),焊接凹接頭附近鋼軌、軌枕加速度、位移較大,加速度最大值分別為1 496 m/s2、127 m/s2,位移最大值分別為 0.77 mm、0.23 mm,比無(wú)不平順時(shí)增加許多,會(huì)影響鋼軌使用壽命。本文計(jì)算結(jié)果在波形上或量值上均與文獻(xiàn)[16]結(jié)果接近,說(shuō)明本文所提建立車(chē)輛多體系統(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
為方便應(yīng)用哈密爾頓原理建立車(chē)輛多體系統(tǒng)振動(dòng)方程,本文工作為:
(1)以三自由度振動(dòng)系統(tǒng)為例,對(duì)哈密爾頓原理“對(duì)號(hào)入座”進(jìn)行改進(jìn),使之與有限元分析中計(jì)算機(jī)編碼法統(tǒng)一,快速利用計(jì)算機(jī)將單元矩陣自動(dòng)組成系統(tǒng)總體矩陣。
(2)建立全車(chē)-軌道垂向振動(dòng)模型,應(yīng)用改進(jìn)后“對(duì)號(hào)入座”法則建立該振動(dòng)系統(tǒng)方程,并應(yīng)用MATLAB編制相應(yīng)計(jì)算程序。
(3)用所編程序?qū)附影冀宇^激勵(lì)下車(chē)輛-軌道垂向振動(dòng)特性進(jìn)行分析,結(jié)果表明,本文計(jì)算結(jié)果與之前研究結(jié)論較一致,表明本文所提方法的正確性。
[1]翟婉明.車(chē)輛-軌道耦合動(dòng)力學(xué)(第三版)[M].北京:科學(xué)出版社,2007.
[2]翟婉明.車(chē)輛-軌道垂向系統(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ī)車(chē)車(chē)輛軌道系統(tǒng)動(dòng)力學(xué)[M].北京:中國(guó)鐵道出版社,1996.
[4]倪純雙,洪嘉振,賀啟庸.車(chē)輛系統(tǒng)動(dòng)力學(xué)的符號(hào)化多體方法[J].中國(guó)鐵道科學(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é)總勢(shì)能不變值原理[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ì)能不變值原理與列車(chē)橋梁時(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é)毅,王 平.車(chē)輛-軌道-路基系統(tǒng)動(dòng)力學(xué)[M].成都:西南交通大學(xué)出版社,2010.
[9]曾慶元,楊 平.形成矩陣的“對(duì)號(hào)入座”法則與桁梁空間分析的桁段有限元法[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ī)車(chē)車(chē)輛-軌道系統(tǒng)運(yùn)動(dòng)方程的建立[J].長(zhǎng)沙鐵道學(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]李東平,曾慶元,婁 平.車(chē)輛多體系統(tǒng)動(dòng)力學(xué)方程的有限元法[J].中國(guó)鐵道科學(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ù)測(cè)-校正積分模式[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]李成輝,萬(wàn)復(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]王偉華.土路基上雙塊式無(wú)砟軌道垂向動(dòng)力特性分析[D].成都:西南交通大學(xué),2009.