高學(xué)軍，李映輝，樂(lè) 源

(1.成都理工大學(xué) 環(huán)境與土木工程學(xué)院，成都 610059;2.西南交通大學(xué) 力學(xué)與工程學(xué)院，成都 610031)

隨著高速列車(chē)的普及和各種動(dòng)車(chē)組的逐步開(kāi)行，車(chē)輛動(dòng)力學(xué)的研究也變得相當(dāng)重要和迫切，其中橫向運(yùn)動(dòng)穩(wěn)定性研究[1－3]是車(chē)輛系統(tǒng)動(dòng)力學(xué)性能非常重要的一個(gè)方面，因?yàn)樗婕暗杰?chē)輛臨界速度的確定，直接影響到列車(chē)允許的最高運(yùn)行速度。Isaksen等[4－5]基于Cooperrider轉(zhuǎn)向架建立了自己的研究模型并對(duì)其分岔行為分析后發(fā)現(xiàn)，亞臨界Hopf分岔現(xiàn)象在車(chē)輛系統(tǒng)中很常見(jiàn)，這一結(jié)論對(duì)傳統(tǒng)的由線性穩(wěn)定性方法確定臨界速度的方式提出了疑問(wèn)。Ahmadian等[6]用漸近法研究了單個(gè)轉(zhuǎn)向架的穩(wěn)定性，并對(duì)諸多影響因素進(jìn)行了分析。楊紹普等[7]在其著作中則研究了具有滯后非線性懸掛的轉(zhuǎn)向架和機(jī)車(chē)車(chē)輛的Hopf分岔行為和運(yùn)動(dòng)穩(wěn)定性。Gao等[8]則研究了四軸客車(chē)系統(tǒng)橫向運(yùn)動(dòng)的分岔行為，發(fā)現(xiàn)系統(tǒng)存在大量的跳躍和遲滯現(xiàn)象，進(jìn)一步的研究結(jié)果表明系統(tǒng)經(jīng)過(guò)多次的音叉分岔反復(fù)的經(jīng)歷對(duì)稱性破壞和對(duì)稱性恢復(fù)的過(guò)程，最終進(jìn)入不對(duì)稱的運(yùn)動(dòng)狀態(tài)。這些研究中均未考慮軌道結(jié)構(gòu)的振動(dòng)影響，計(jì)算模型與實(shí)際存在一定差別，因其自由度較少，可用來(lái)分析車(chē)輛系統(tǒng)的各種分岔行為，為臨界速度的確定提供指導(dǎo)和依據(jù)。實(shí)際上，在車(chē)輛運(yùn)行中，鋼軌所受的橫向載荷將使其產(chǎn)生橫向變形，Garg等[9]注意到了這一問(wèn)題并將鋼軌簡(jiǎn)化為橫向剛度和阻尼連接的元件計(jì)算鋼軌在橫向發(fā)生的撓曲，該模型可在考慮軌道整體橫向彈性和阻尼的情況下，研究車(chē)輛系統(tǒng)的運(yùn)行穩(wěn)定性和曲線通過(guò)性能。

基于此，本文在傳統(tǒng)車(chē)輛動(dòng)力學(xué)模型基礎(chǔ)上，將鋼軌考慮成具有一定橫向剛度和阻尼，且隨各個(gè)輪對(duì)一起向前運(yùn)動(dòng)的離散剛體。應(yīng)用延續(xù)算法對(duì)包含有簡(jiǎn)單軌道的客車(chē)系統(tǒng)橫向運(yùn)動(dòng)的分岔問(wèn)題進(jìn)行研究，探討該軌道模型對(duì)車(chē)輛系統(tǒng)分岔行為的影響規(guī)律，為車(chē)輛系統(tǒng)穩(wěn)定性研究提供部分參考和依據(jù)。

1 車(chē)輛系統(tǒng)數(shù)學(xué)模型

1.1 簡(jiǎn)單軌道模型

圖1 考慮鋼軌彈性和阻尼的簡(jiǎn)單軌道模型Fig.1 The track model considering the elasticity and damping of rails

如圖1所示將鋼軌視為具有一定橫向剛度和阻尼，且隨各個(gè)輪對(duì)一起向前運(yùn)動(dòng)的離散剛體[9]，每個(gè)輪對(duì)下左、右各存在這樣一個(gè)離散鋼軌，其實(shí)際作用相當(dāng)于車(chē)輛懸掛的最后幾個(gè)剛體。由于鋼軌的彈性力與阻尼力項(xiàng)遠(yuǎn)大于慣性力項(xiàng)，因此忽略鋼軌的等效質(zhì)量，故左右鋼軌(i=1－4，下同)的橫向運(yùn)動(dòng)微分方程可表示為:

式中:Cr、Kr為鋼軌的橫向阻尼和橫向剛度;yrli，yrri分別為第 i個(gè)輪對(duì)下左、右鋼軌的橫向位移;Fcyli，F(xiàn)nyli，F(xiàn)tli分別為輪對(duì)對(duì)左軌的蠕滑力、法向力、輪緣力在橫向的分量;Fcyri，F(xiàn)nyri，F(xiàn)tri則為輪對(duì)對(duì)右軌的蠕滑力、法向力、輪緣力在橫向的分量。

1.2 輪軌接觸幾何關(guān)系

輪軌接觸幾何參數(shù)主要包括左/右輪滾動(dòng)圓半徑(rl，rr)，左/右輪輪軌接觸角(δl，δr)以及輪對(duì)側(cè)滾角(φw)等，它們是輪對(duì)橫移量 yw的函數(shù)[9]:

式中:r0是車(chē)輪名義滾動(dòng)圓半徑，λ是車(chē)輪踏面等效錐度，δ0是輪對(duì)居中時(shí)的左/右輪接觸角，ε0則是輪對(duì)橫移引起的接觸角變化參數(shù)，a是輪軌接觸點(diǎn)橫向距離之半，σ為輪對(duì)側(cè)滾角參數(shù)。

1.3 非線性輪軌接觸力

輪軌之間的非線性接觸力主要包括輪軌滾動(dòng)接觸點(diǎn)的蠕滑力和輪緣力。

輪軌滾動(dòng)接觸蠕滑力和力矩首先根據(jù)Kalker線性蠕滑理論計(jì)算，然后采用沈氏蠕滑理論進(jìn)行非線性修正。輪軌蠕滑力與蠕滑率有關(guān)，縱向、橫向、自旋蠕滑率對(duì)左輪(l，?取上面的符號(hào))，右輪(r，?取下面的符號(hào))可依次表示為[10]:

式中:V是車(chē)輛前行速度，ψw是輪對(duì)搖頭角位移。

根據(jù)Kalker線性蠕滑理論[11]，輪軌之間的接觸斑蠕滑力在線性范圍內(nèi)可表達(dá)為:

式中:Fx、Fy、Mz分別為接觸斑縱向、橫向蠕滑力和旋轉(zhuǎn)蠕滑力矩，而 f11、f22、f23、f33分別為與輪軌接觸幾何參數(shù)有關(guān)的縱向、橫向、橫向/旋轉(zhuǎn)、自旋蠕滑系數(shù)。

考慮到輪軌接觸面可能存在大蠕滑的情形，輪軌蠕滑力和蠕滑力矩可采用沈氏蠕滑理論[12]進(jìn)行非線性修正，將修正后的接觸斑蠕滑力通過(guò)坐標(biāo)變換轉(zhuǎn)換到軌道坐標(biāo)系內(nèi)[13]，即可用于運(yùn)動(dòng)微分方程的建立。

左、右輪輪緣力Ft(l，r)一般用有死區(qū)的剛性彈簧來(lái)模擬[14]，表達(dá)為:

式中:k0是彈性系數(shù)，η代表輪緣間隙大小。Yl=ywyrl，Yr=yw－yrr分別為輪對(duì)相對(duì)于左、右鋼軌的有效橫移量。

1.4 車(chē)輛系統(tǒng)

考慮一四軸鐵道客車(chē)系統(tǒng)如圖2所示，主要包括車(chē)體(Mc，Icx，Icz)、兩個(gè)構(gòu)架(Mt，Itx，Itz)、4 個(gè)輪對(duì)(Mw，Iwx，Iwz)等剛體和一系(Kpx，Kpy，Kpz，Cpx，Cpy，Cpz)及二系(Ksx，Ksy，Ksz，Csx，Csy，Csz)懸掛的多剛體動(dòng)力學(xué)系統(tǒng)，除懸掛系統(tǒng)具有線彈性特性外，其它部分都是剛性的。為簡(jiǎn)化分析計(jì)算，假設(shè)垂向運(yùn)動(dòng)和橫向運(yùn)動(dòng)是解耦的，研究中主要考慮整個(gè)系統(tǒng)的橫向運(yùn)動(dòng)。

圖2 四軸客車(chē)橫向穩(wěn)定性模型Fig.2 Model of four-axle railway passenger car system

整個(gè)系統(tǒng)的運(yùn)動(dòng)主要包括:4個(gè)輪對(duì)的橫移ywi(i=1～4)和搖頭Ψwi運(yùn)動(dòng);兩個(gè)構(gòu)架的橫移ytj(j=1～2)、側(cè)滾 φtj及搖頭 Ψtj運(yùn)動(dòng);車(chē)體的橫移 yc，側(cè)滾 φc及搖頭Ψc運(yùn)動(dòng);各個(gè)輪對(duì)下左右鋼軌的橫移運(yùn)動(dòng)yrli，yrri。至此則可以導(dǎo)出整個(gè)系統(tǒng)的運(yùn)動(dòng)微分方程為:

式中:x∈R42為包含系統(tǒng)各剛體位移與速度的狀態(tài)向量，V∈R+代表車(chē)輛運(yùn)行速度，t表示時(shí)間。

2 分析方法

在一定速度下，由式(6)確定的車(chē)輛系統(tǒng)存在定常運(yùn)動(dòng)和周期運(yùn)動(dòng)，甚至還有可能出現(xiàn)混沌運(yùn)動(dòng)。在這里，通過(guò)討論如下含有參數(shù)λ，更具普遍意義的常微分方程組:

解的計(jì)算問(wèn)題，以此來(lái)說(shuō)明車(chē)輛系統(tǒng)定常解和周期解的延續(xù)計(jì)算問(wèn)題。

2.1 定常解的求解

方程組(7)的定常解由下式確定:

一般來(lái)說(shuō)，除去一些特別的點(diǎn)，定常解是參數(shù)λ的函數(shù)，這些解形成了系統(tǒng)的解分支曲線。如果f(x，λ)是非線性的，對(duì)確定的參數(shù)值λ，系統(tǒng)可能不止有一個(gè)零點(diǎn)，假設(shè)初始參數(shù)值λ0是確定的并有 x0=x(λ0)，則式f(x0，λ0)=0成立。令:

式中:fx實(shí)際是系統(tǒng)的Jacobi矩陣。

為了找到λ =λ0+Δλ(Δλ是足夠小步長(zhǎng))處解分支上的點(diǎn)，使用預(yù)測(cè)校正算法，對(duì)f[x(λ)，λ]=0關(guān)于λ進(jìn)行微分并變換得:

故x(λ0+Δλ)處的預(yù)測(cè)值為:

以該值作為初始值利用Newton-Raphson法求解式(8)即可得到x(λ0+Δλ)處的計(jì)算值。

如果fx奇異，或遇到解曲線的分岔點(diǎn)時(shí)，上面的算法可能會(huì)失效。這時(shí)可以考慮如下形式的擴(kuò)展系統(tǒng)[15]:

處理后即便fx奇異，通過(guò)選擇合適的n(x，λ，s)仍可保證Fy(y0，s0)非奇異，因此預(yù)測(cè)校正算法仍可以應(yīng)用于以(y0，s0)為初始點(diǎn)F(y，s)=0的系統(tǒng)，此時(shí)預(yù)測(cè)值表達(dá)式變?yōu)?

式中s為弧長(zhǎng)參數(shù)，該方法因此也稱為“偽弧長(zhǎng)延續(xù)算法”

n(x，λ，s)的選擇并不唯一，Doedel[16]使用的 n(x，λ，s)形式為:

以(x0，λ0)為初始點(diǎn)，利用式(13)進(jìn)行預(yù)測(cè)，后求解式(12)即可得到解分支上的下一點(diǎn)，不斷重復(fù)此過(guò)程則可以得到系統(tǒng)的定常解分支曲線，定常解的穩(wěn)定性則可以通過(guò)系統(tǒng)Jacobi矩陣特征值實(shí)部的正負(fù)來(lái)判定。

2.2 周期解計(jì)算

在一定條件下，方程組(7)存在周期解且滿足邊值問(wèn)題常微分方程組:

式中:T為系統(tǒng)周期，無(wú)量綱時(shí)間τ=t/T。

顯然上式是因缺少方程而無(wú)法直接求解，需要增加一個(gè)積分相位條件，即:

考慮到偽弧長(zhǎng)延續(xù)算法的優(yōu)勢(shì)，在上述方程組的基礎(chǔ)上增加一偽弧長(zhǎng)延續(xù)方程[17]:

即可應(yīng)用偽弧長(zhǎng)延續(xù)算法來(lái)求解系統(tǒng)的周期解，周期解的穩(wěn)定性可通過(guò)系統(tǒng)的特征乘子是否位于單位圓內(nèi)來(lái)說(shuō)明。

3 數(shù)值結(jié)果

以四方車(chē)輛廠生產(chǎn)的某準(zhǔn)高速車(chē)為分析對(duì)象，車(chē)輛系統(tǒng)中各剛體的質(zhì)量值、剛度系數(shù)、阻尼系數(shù)等參數(shù)的取值見(jiàn)文獻(xiàn)[13]。輪軌接觸幾何參數(shù)為:車(chē)輪踏面等效錐度λ=0.05;輪對(duì)居中時(shí)的左/右輪輪軌接觸角δ0=0.05，輪對(duì)橫移引起的接觸角變化參數(shù) ε0=0，輪對(duì)側(cè)滾角參數(shù)σ=0.05。輪軌接觸力計(jì)算參數(shù)為:輪緣力剛性彈簧系數(shù)k0=14.6×107N/m，輪緣間隙η=9.1 mm;蠕滑系數(shù) f11=1.023 ×107N，f22=9.43 ×106N，f23=1 200 N·m，f33=1 000 N·m2;輪軌粘著系數(shù)μ=0.15。鋼軌橫向彈性系數(shù) Kr=14.6×107N/m，阻尼系數(shù)Cr=14.6×104N·s/m。車(chē)輛運(yùn)行速度變化范圍0＜V＜160.0 m/s。

圖3是應(yīng)用延續(xù)算法計(jì)算得到的車(chē)輛運(yùn)行速度作為控制參數(shù)與前轉(zhuǎn)向架前導(dǎo)輪對(duì)最大橫向幅值分岔圖，其中實(shí)線代表穩(wěn)定的運(yùn)動(dòng)，而點(diǎn)線則代表不穩(wěn)定的運(yùn)動(dòng)。具體可從以下兩個(gè)方面進(jìn)行分析:

(1)首先對(duì)車(chē)輛系統(tǒng)的定常解進(jìn)行說(shuō)明:由于分析的是直線軌道，因此零解必定是其定常解，這從圖中零附近的定常解分支曲線OQ也可反映出來(lái)。點(diǎn)A(VA=119.385 m/s，α1，2= － 3.794 × 10－9± 22.499i，αi代表系統(tǒng)Jacobi矩陣實(shí)部最大的特征值，下同)是車(chē)輛系統(tǒng)第1個(gè)Hopf分岔點(diǎn)，此時(shí)第一對(duì)復(fù)共軛特征值正向穿越虛軸，穩(wěn)定的定常運(yùn)動(dòng)(OA段實(shí)線)因此變得不穩(wěn)定(AC段點(diǎn)線);隨著速度的繼續(xù)增大，系統(tǒng)Jacobi矩陣特征值實(shí)部基本是增加的，當(dāng)車(chē)輛運(yùn)行速度達(dá)到VC(VC=121.430 m/s，α1，2=0.167 ± 22.751i，α3，4=2.547 ×10－9±22.985i)時(shí)，又有一對(duì)復(fù)共軛特征值正向穿越虛軸而使定常解(CQ段點(diǎn)線)變得“更加不穩(wěn)定”。之后隨著速度的繼續(xù)增大直到速度終止值，再?zèng)]有發(fā)現(xiàn)其它的Hopf分岔點(diǎn)。

圖3 前轉(zhuǎn)向架前導(dǎo)輪對(duì)最大橫向幅值分岔圖Fig.3 Bifurcation diagram showing the maximum lateral amplitude of the leading wheel set of the front bogie)

(2)其次對(duì)車(chē)輛系統(tǒng)的周期解進(jìn)行分析:分別以Hopf分岔點(diǎn)A和點(diǎn)C為初始點(diǎn)進(jìn)行周期解的延拓計(jì)算發(fā)現(xiàn)，系統(tǒng)從點(diǎn)A處以亞臨界方式分岔出一不穩(wěn)定的周期解曲線(AB段點(diǎn)線)，該周期解的幅值隨著速度的減小而逐漸增加并在點(diǎn) B(VB=88.209 m/s，s，T表示系統(tǒng)振動(dòng)周期，下同)處前導(dǎo)輪對(duì)與鋼軌出現(xiàn)了輪緣接觸。與此同時(shí)，系統(tǒng)從點(diǎn)C處仍以亞臨界方式分岔出另一不穩(wěn)定的周期解曲線(CD段點(diǎn)線)，在分岔點(diǎn)0.273 4 s，TD=0.395 1 s)處不穩(wěn)定的周期運(yùn)動(dòng)運(yùn)動(dòng)逐漸恢復(fù)了穩(wěn)定。之后隨著速度的增加，穩(wěn)定的周期運(yùn)動(dòng)(DE段實(shí)線)幅值逐漸增加直到到達(dá)非周期點(diǎn)因車(chē)輛臨界速度的確定實(shí)際上是在運(yùn)行速度－狀態(tài)向量空間中尋找定常解是唯一漸近穩(wěn)定解時(shí)車(chē)輛運(yùn)行速度最高值的過(guò)程[5]，故包含有簡(jiǎn)單軌道的該客車(chē)系統(tǒng)的臨界速度為圖3中點(diǎn)D處的速度Vcr=84.344 m/s。

圖4(a)(b)分別給出了圖3中AB解分支和CD解分支上不穩(wěn)定的極限環(huán)圖，其中橫坐標(biāo)均表示前轉(zhuǎn)向架前導(dǎo)輪對(duì)橫向位移，而縱坐標(biāo)均表示該輪對(duì)的橫移速度。圖4(a)中極限環(huán)由外向內(nèi)對(duì)應(yīng)的車(chē)輛運(yùn)行速度分別為88.209，90，95，100，105，110，115 m/s，該分支的極限環(huán)最后收縮到中心的定常吸引子點(diǎn)VA=119.385 m/s上;而圖4(b)中極限環(huán)由外向內(nèi)對(duì)應(yīng)的車(chē)輛運(yùn)行速度分別為90，84.344，95，100，105，110，115，120 m/s，該分支的極限環(huán)最后收縮到中心的定常吸引子點(diǎn)VC=121.430 m/s上。需要說(shuō)明的是，這些不穩(wěn)定的解集只能應(yīng)用延續(xù)算法算出，普通的數(shù)值積分方法是無(wú)法得到不穩(wěn)定解集的。

通過(guò)對(duì)不考慮鋼軌彈性和阻尼、但同等懸掛條件下的傳統(tǒng)車(chē)輛動(dòng)力學(xué)模型研究發(fā)現(xiàn)其分岔曲線與包含簡(jiǎn)單軌道的車(chē)輛動(dòng)力學(xué)模型相似，只是分岔點(diǎn)的速度值有所區(qū)別。表1給出了兩種模型分岔點(diǎn)速度計(jì)算結(jié)果比較，從中可以看出，不管是第1，2個(gè)Hopf分岔點(diǎn)還是車(chē)輛臨界速度值，兩種模型計(jì)算結(jié)果相差甚微，且簡(jiǎn)單軌道模型各個(gè)速度值比傳統(tǒng)車(chē)輛模型相對(duì)應(yīng)速度值要稍高些，這與文獻(xiàn)[2]的結(jié)論也一致。與文獻(xiàn)[2]不同之處是，文獻(xiàn)[2]得到的往往是第一個(gè)Hopf分岔點(diǎn)值和臨界速度值，并不能確定系統(tǒng)中不穩(wěn)定的周期解，而本文應(yīng)用延續(xù)算法不但確定了系統(tǒng)兩個(gè)Hopf分岔點(diǎn)值，而且還得到了從這兩個(gè)點(diǎn)分岔出的不穩(wěn)定周期解曲線，最后根據(jù)解的分岔情況很容易確定了車(chē)輛系統(tǒng)的臨界速度。

圖4 兩個(gè)解分支上的不穩(wěn)定極限環(huán)Fig.4 Unstable limit cycles of the two branches

表1 兩種模型分岔點(diǎn)速度比較Tab.1 Comparison of the bifurcation speeds in two models

圖5給出了鋼軌橫向阻尼Cr=6.5×104N·s/m的情況下，鋼軌橫向剛度對(duì)分岔速度的影響曲線。其中HP1表示第1個(gè)Hopf分岔點(diǎn)值，HP2表示第2個(gè)Hopf分岔點(diǎn)值，CR則表示車(chē)輛臨界速度值。從中可以看出，這三種分岔速度值都隨著鋼軌橫向剛度的增大而減小，當(dāng)鋼軌橫向剛度值超過(guò)50×106N/m后，該剛度值的大小對(duì)分岔速度值影響不再顯著。

圖5 鋼軌橫向剛度對(duì)分岔速度的影響Fig.5 Effect of lateral stiffness of rails on the bifurcation speed

圖6 則給出了鋼軌橫向剛度Kr=29.4×106N/m的情況下，鋼軌橫向阻尼對(duì)分岔速度的影響曲線。由圖可見(jiàn)，增大或減小橫向阻尼對(duì)三種分岔點(diǎn)速度值的影響都很有限，即三種分岔點(diǎn)速度對(duì)鋼軌橫向阻尼值的變化不敏感。

圖6 鋼軌橫向阻尼對(duì)分岔速度的影響Fig.6 Effect of lateral damping of rails on the bifurcation speed

4 結(jié)論

本文討論了車(chē)輛系統(tǒng)運(yùn)動(dòng)微分方程組定常解和周期解的延續(xù)計(jì)算方法，并應(yīng)用該方法研究了包含有簡(jiǎn)單軌道的客車(chē)系統(tǒng)的橫向運(yùn)動(dòng)分岔問(wèn)題。分析中將軌道考慮成具有整體橫向剛度和阻尼，且隨各個(gè)輪對(duì)一起向前運(yùn)動(dòng)的離散鋼軌。在車(chē)輛運(yùn)行速度變化范圍內(nèi)，借助延續(xù)算法方便快捷的得到了車(chē)輛系統(tǒng)的兩個(gè)Hopf分岔點(diǎn)以及從這兩個(gè)點(diǎn)分岔出來(lái)的周期解曲線，根據(jù)分岔曲線的走向確定了車(chē)輛系統(tǒng)的臨界速度。同時(shí)，也將相關(guān)結(jié)果與相同懸掛條件下，不考慮鋼軌的傳統(tǒng)車(chē)輛動(dòng)力學(xué)模型計(jì)算結(jié)果進(jìn)行了對(duì)比。結(jié)果表明簡(jiǎn)單軌道模型對(duì)車(chē)輛系統(tǒng)的分岔曲線變化趨勢(shì)及臨界速度值影響不大，且考慮簡(jiǎn)單軌道的車(chē)輛系統(tǒng)兩個(gè)Hopf分岔點(diǎn)值和臨界速度值都要比同等懸掛條件下傳統(tǒng)車(chē)輛動(dòng)力學(xué)模型的計(jì)算結(jié)果稍高些。而高出的程度主要與鋼軌的橫向剛度大小有關(guān)，鋼軌的橫向阻尼對(duì)其影響則很小。研究結(jié)果可為高速客車(chē)的設(shè)計(jì)、穩(wěn)定性分析中車(chē)輛系統(tǒng)的建模以及分岔的研究提供部分理論依據(jù)。

[1]Knothe K，Bohm F.History of stability of railway and road vehicles[J].Vehicle System Dynamics，1999，31(5):283－323.

[2]劉宏友.高速列車(chē)中的關(guān)鍵動(dòng)力學(xué)問(wèn)題研究[D].成都:西南交通大學(xué)，2003.

[3]Zeng J，Wu P B.Stability analysis of high speed railway vehicles[J].JSME International Journal，Series C，2004，47(2):464－470.

[4]Isaksen P，True H.On the ultimate transition to chaos in the dynamics of Cooperrider's bogie[J].Chaos，Solitons ＆Fractals，1997，8(4):559 －581.

[5]True H.On the critical speed of high-speed railway vehicles[J].Noise and Vibration on High-Speed Railways.FEUP Porto Portugal，2008:149 －166.

[6]Ahmadian M，Yang S P.Effect of system nonlinearities on locomotive bogie hunting stability[J].Vehicle System Dynamics，1998，29(6):365－384.

[7]楊紹普，申永軍.滯后非線性系統(tǒng)的分岔與奇異性[M].北京:科學(xué)出版社，2003.

[8]Gao X J，Li Y H，Gao Q.Lateral bifurcation behavior of a four-axle railway passenger car[J].Transaction of the ASME，Journal of Applied Mechanics，2010，77(6):061001:1－8.

[9]Garg V K，Dukkipati R V.Dynamics of railway vehicle systems[M].New York:Academic Press，1984.

[10]Petersen D E，Hoffmann M.Curving dynamics of railway vehicles[R].Informatics and Mathematical Modellng，The Technical University of Demnark，Lyngby，2002.

[11]Kalker J J.A fast algorithm for the simplified theory of rolling contact[J].Vehicle System Dynamics，1982，11(1):1 －13.

[12]Shen Z Y，Hedrick J K，Elkins J A.A comparison of alternative creep force models for rail vehicle dynamic analysis[C].Proceeding of 8th IAVSD Symposium on Vehicle System Dynamics，Dynamics of Vehicles on Roads and Tracks.MIT，Cambridge:Swets and Zeitlinger，1984:591－605.

[13]翟婉明.車(chē)輛－軌道耦合動(dòng)力學(xué)[M].北京:中國(guó)鐵道出版社，1997.

[14]Kaas-Petersen C.Chaos in a railway bogie[J].Acta Mechanica，1986，61(1－4):89－107.

[15]Brindley J，Kaas-Petersen C，Spence A.Path-following methods in bifurcation problems[J].Physica D:Nonlinear Phenomena，1989，34(3):456－461.

[16]Doedel E J.Auto-07P:continuation and bifurcation software for ordinary differential equations[R].California Institute of Technology，2008.

[17]Krauskopf B，Osinga H M，Galan-Vioque J.Numerical continuation methods for dynamical systems[M].Path following and Boundary Value Problems:Springer，2007.