基于蠕滑理論的蛇形運(yùn)動分岔問題研究

0蛇形運(yùn)動穩(wěn)定性研究當(dāng)鐵路車輛沿軌跡移動時(shí)，具有一定形狀的轉(zhuǎn)向車輪對角（水平移動），旋轉(zhuǎn)至銅軸的棕色軸（旋轉(zhuǎn)）。這兩種運(yùn)動的結(jié)合使車輛沿著正確的軌道前進(jìn)。這種現(xiàn)象通常被稱為蛇形。蛇形是一種非常普遍的自激振橫向擺動不穩(wěn)定現(xiàn)象,主要是由輪軌之間的非線性接觸力及懸掛裝置作用而引起。高速車輛蛇形失穩(wěn)不僅惡化車輛運(yùn)行性能,降低乘坐舒適性,而且還會加劇部件磨損,對車輛及線路造成損傷,甚至?xí)l(fā)脫軌的重大事故,因此,車輛系統(tǒng)蛇形運(yùn)動穩(wěn)定性的研究具有重要意義。國內(nèi)外學(xué)者對車輛系統(tǒng)橫向穩(wěn)定性及分岔問題進(jìn)行了大量研究:Wickens推導(dǎo)了單個(gè)輪對的運(yùn)動方程,并從動力學(xué)角度詳細(xì)分析了輪對的失穩(wěn)問題;Kaas-Petersen分析了Cooperrider構(gòu)架在有輪緣和無輪緣情況下的對稱擺振、不對稱擺振以及混沌行為;True在Kaas-Petersen研究的基礎(chǔ)對穩(wěn)定性及分岔問題提出了很多好的方法和建議;曾京使用QR算法與黃金分割法確定車輛系統(tǒng)的Hopf分岔點(diǎn),用打靶法進(jìn)行極限環(huán)的數(shù)值求解,較好地解決了相關(guān)的數(shù)值計(jì)算問題;Ahmadian等用漸近法研究了單個(gè)轉(zhuǎn)向架的穩(wěn)定性,并對諸多影響因素進(jìn)行了分析;司道林等應(yīng)用多剛體動力學(xué)軟件SIMPACK和數(shù)值積分技術(shù),得到了牽引力對車輛臨界速度的影響規(guī)律。在這些研究中,對簡單模型研究比較深入,對復(fù)雜鐵道車輛系統(tǒng)的研究往往側(cè)重Hopf分岔的參數(shù)條件,對Hopf分岔后周期解的變化過程研究較少?；诖?本文對一六軸機(jī)車系統(tǒng)進(jìn)行研究,運(yùn)用延續(xù)算法對Hopf分岔及分岔后極限環(huán)的解分支進(jìn)行連續(xù)跟蹤,討論Hopf分岔的類型和分岔過程中出現(xiàn)的多種非線性動力學(xué)現(xiàn)象,并闡述其中的工程意義。1力學(xué)模型1.1橫向運(yùn)動方程圖1中的車輛模型包含車體(Mc,Icx,Icz)、2個(gè)三軸轉(zhuǎn)向架(Mt,Itx,Itz)、6個(gè)輪對(Mw,Iwz)和一系懸掛(Kpx,Kpy,Kpz,Cpx,Cpy,Cpz)和二系懸掛(Ksx,Ksy,Ksz,Csx,Csy,Csz)的六軸機(jī)車系統(tǒng),暫不考慮牽引電機(jī)的振動,一系懸掛、二系懸掛均視為線性彈簧和阻尼的組合。車輛系統(tǒng)不考慮軌道不平順帶來的影響,假設(shè)輪子在足夠光滑、水平的直線鋼軌上滾動。由于機(jī)車系統(tǒng)前后、左右對稱,可以近似地把橫向運(yùn)動方程和垂向運(yùn)動方程解耦,得到獨(dú)立的橫向運(yùn)動方程組。整個(gè)車輛系統(tǒng)共有21個(gè)自由度:車體考慮三個(gè)自由度,即橫移yc、側(cè)滾?c和搖頭ψc運(yùn)動;每個(gè)轉(zhuǎn)向架考慮構(gòu)架的橫移ytj(j=1,2)、側(cè)滾?tj及搖頭ψtj運(yùn)動;對于輪對則只考慮橫移ywi(i=1,…,6)和搖頭ψwi運(yùn)動。應(yīng)用牛頓定律可導(dǎo)出整個(gè)系統(tǒng)的運(yùn)動微分方程,用如下形式的一階非線性常微分方程組來表示dzdt=f(z,V)(1)dzdt=f(z,V)(1)式中:z∈R42為包含系統(tǒng)各剛體位移與速度的狀態(tài)向量;V∈R+為車輛運(yùn)行速度;t為時(shí)間。1.2輪軌接觸力與輪緣力輪軌的相互作用關(guān)系主要包括輪軌接觸幾何關(guān)系和非線性輪軌接觸力。傳統(tǒng)的機(jī)車車輪踏面都設(shè)計(jì)成錐面,以減小通過曲線時(shí)的磨耗和有助于自導(dǎo)向,輪軌接觸幾何參數(shù),如左、右輪滾動圓半徑(rL,rR),輪軌接觸角(δL,δR)以及輪對側(cè)滾角位移(φw)都可近似認(rèn)為是輪對橫移量yw的非線性函數(shù){rL=r0+κywrR=r0-κywδL=δR=κφw=κyw/a(2)???????????rL=r0+κywrR=r0?κywδL=δR=κφw=κyw/a(2)式中:r0為車輪名義滾動圓半徑;a為輪軌接觸點(diǎn)橫向距離之半;κ為接觸角。非線性輪軌接觸力主要包括輪軌滾動接觸蠕滑力和輪緣力。本文采用Vermeulen-Johnson理論計(jì)算輪軌滾動接觸蠕滑力,并充分考慮蠕滑力的飽和效應(yīng)(但不考慮自旋蠕滑力)。輪對的縱向和橫向蠕滑率分別為{ξx=a˙ψw/V+κyw/r0ξy=˙yw/V-ψw(3){ξx=aψ˙w/V+κyw/r0ξy=y˙w/V?ψw(3)式中:˙ywy˙w為輪對橫移速度;ψw、˙ψwψ˙w分別為輪對搖頭角位移與速度。合成蠕滑率為ξR=√(ξx/φ)2+(ξy/ψ)2(4)ξR=(ξx/φ)2+(ξy/ψ)2???????????????√(4)式中:φ、ψ為利用Hertz接觸理論從Johnson公式中計(jì)算出的權(quán)重系數(shù)?？v向和橫向蠕滑力分別為{Fx=ξxFR/(ξRφ)Fy=ξyFR/(ξRψ)(5)FR=μΝ?{u-u2/3+u3/27u<31u≥3(6)u=ξRGπaebe/(μΝ)(7){Fx=ξxFR/(ξRφ)Fy=ξyFR/(ξRψ)(5)FR=μN(yùn)?{u?u2/3+u3/271u<3u≥3(6)u=ξRGπaebe/(μN(yùn))(7)式中:N為輪軌接觸斑的法向力;G為輪軌合成剪切模量;ae、be分別為輪軌接觸橢圓的長、短半軸長度;μ為輪軌間的庫侖摩擦因數(shù);FR為修正后的合成蠕滑力;u為經(jīng)歸一化修正后的合成蠕滑率。輪緣力一般用有死區(qū)的剛性彈簧來模擬,即可表達(dá)為一分段線性函數(shù)形式FΤ={k0(yw-η)yw>η0|yw|≤ηk0(yw+η)yw<-η(8)FT=?????k0(yw?η)0k0(yw+η)yw>η|yw|≤ηyw<?η(8)式中:k0為彈性系數(shù);η為輪緣間隙,不同輪對的輪緣間隙大小可能不一樣。2基于生物活性的新能源更新算法一般來說,在一定速度下,由式(1)確定的車輛系統(tǒng)存在穩(wěn)態(tài)擺動和周期擺動,對應(yīng)于微分方程組的定常解和周期解,因此,下面討論如何求解式(1)的定常解和周期解。定常解與周期解分支曲線都屬于解的分岔問題。在這里,不失一般性,假設(shè)x為n維狀態(tài)向量,λ為一般參數(shù),定義映射F:Rn×R→Rn,力求找到所有的(x,λ),使其滿足式為F(x,λ)=0x,F∈Rn,λ∈R(9)F(x,λ)=0x,F∈Rn,λ∈R(9)如果要確定式(1)的定常解,可令F(x,λ)=f(z,V)(10)如果要尋找式(1)的周期解,只須讓F(x,λ)=Ρ(z,V)-z(11)式中:P為Poincare映射的函數(shù)。換句話說,現(xiàn)在只需關(guān)心式(9)的求解問題了。延續(xù)算法是從系統(tǒng)的一個(gè)已知解出發(fā),通過預(yù)測-校正的機(jī)制可逐步求解出整個(gè)系統(tǒng)的解分支曲線。假如(x0,λ0)是解分支曲線上的已知點(diǎn),為了找到λ=λ0+Δλ(Δλ是足夠小步長)處解分支上的點(diǎn),對F(x,λ)=0關(guān)于λ進(jìn)行微分、變換,可得(dxdλ)0=-[■SymbolΡropBΤvBpF(x0,λ0)■SymbolΡropBΤvBpx]-1■SymbolΡropBΤvBpF(x0,λ0)■SymbolΡropBΤvBpλ(12)故λ0+Δλ處的x預(yù)測值為xp=x0+Δλ(dxdλ)0(13)以這個(gè)值作為初始值,利用Newton-Raphson迭代法求解式(9),可得到λ0+Δλ處x的計(jì)算值。如果ue851SymbolPropBTvBpF/ue851SymbolPropBTvBpx奇異,或者遇到解曲線的分岔點(diǎn)時(shí),上面的算法可能會失效,這時(shí)可考慮如下形式的擴(kuò)展系統(tǒng){G(y,s)=[F(x,λ)n(x,λ,s)]=0y=(x,λ)(14)式中:y與s分別為擴(kuò)展系統(tǒng)的狀態(tài)向量與參數(shù);n(x,λ,s)為函數(shù)。這樣處理后即便ue851SymbolPropBTvBpF/ue851SymbolPropBTvBpx奇異,通過選擇合適的n(x,λ,s)可使G(y0,s0)非奇異,因此y=(x,λ)在(y0,s0)附近仍是s的光滑函數(shù),預(yù)測校正算法仍可以應(yīng)用于以(y0,s0)為初始點(diǎn)G(y,s)=0的系統(tǒng),只不過此時(shí)預(yù)測值表達(dá)式變?yōu)閧xp=x0+Δs(dxds)0λp=λ0+Δs(dλds)0(15)s為通常的弧長參數(shù),該方法因此也稱為弧長延續(xù)算法。當(dāng)然,n(x,λ,s)的選擇并不唯一,Doedel使用的n(x,λ,s)形式為n(x,λ,s)=(dxds)Τ0(x-x0)+(dλds)0(λ-λ0)-(s-s0)=0(16)以(x0,λ0)為初始點(diǎn),利用式(15)進(jìn)行預(yù)測,然后求解式(14),可得到解分支上的下一點(diǎn),不斷重復(fù)此過程,可以得到系統(tǒng)的解分支曲線。對復(fù)雜的多自由度多強(qiáng)非線性的鐵道車輛系統(tǒng)而言,用理論方法確定式(11)的Poincare映射函數(shù)是非常困難的,因此,大多數(shù)情況下需要用數(shù)值方法得到系統(tǒng)的Poincare映射,本文取前轉(zhuǎn)向架1位輪對與鋼軌碰撞后橫向運(yùn)動速度為零且橫向位移為正的面為Poincare截面σ={(z,V)∈R42×R+,˙yw1=0,yw1>0}(17)通過數(shù)值積分來建立系統(tǒng)的Poincare映射函數(shù)。3車輛運(yùn)行速度的數(shù)值取值本文以戚墅堰機(jī)車車輛廠1990年生產(chǎn)的東風(fēng)9型(DF9)內(nèi)燃機(jī)車為分析對象,機(jī)車系統(tǒng)各剛體的質(zhì)量、剛度系數(shù)、阻尼系數(shù)等參數(shù)的取值見文獻(xiàn),其中抗蛇形減震器卸荷速度取0.043m·s-1,減震器飽和阻力為13kN。輪軌計(jì)算參數(shù)在文獻(xiàn)中均可查到,車輛運(yùn)行速度變化范圍為0<V≤150m·s-1。圖2是應(yīng)用延續(xù)算法得到的以車輛運(yùn)行速度作為控制參數(shù)的1位輪對相對軌道的最大橫向位移分岔圖。點(diǎn)線代表不穩(wěn)定的運(yùn)動,具體分析如下。(1)系統(tǒng)jacabi矩陣特征值實(shí)部點(diǎn)A(VA=53.700m·s-1,α1,2=-1.751×10-10±5.311i,αi代表系統(tǒng)特征值,下同)是系統(tǒng)的第1個(gè)Hopf分岔點(diǎn),此時(shí)有一對復(fù)共軛特征值正向穿越虛軸,定常解因此而失去其原有的穩(wěn)定性;隨著速度的增大,系統(tǒng)Jacobi矩陣特征值實(shí)部基本是增加的;當(dāng)速度達(dá)到VD(VD=114.657m·s-1,α1,2=1.585±8.903i,α3,4=3.520×10-11±9.266i)時(shí),又有一對復(fù)共軛特征值正向穿越虛軸而使定常解變得“更加不穩(wěn)定”;隨著速度的繼續(xù)增大直到速度終止值,再沒有發(fā)現(xiàn)其他的Hopf分岔點(diǎn)。(2)系統(tǒng)的周期運(yùn)動狀態(tài)分別以Hopf分岔點(diǎn)A和點(diǎn)D為初始點(diǎn)進(jìn)行周期解的延拓計(jì)算,發(fā)現(xiàn)系統(tǒng)從點(diǎn)A處以亞臨界的方式分岔出一不穩(wěn)定的周期解AB(有1個(gè)floquet乘子位于單位圓外,對應(yīng)于車輛系統(tǒng)不穩(wěn)定的蛇形運(yùn)動),分岔出的不穩(wěn)定周期解在鞍結(jié)分岔點(diǎn)B(VB=50.855m·s-1,max|yB|=1.500mm,TB=1.2345s,yi為點(diǎn)i的橫移,Ti代表系統(tǒng)運(yùn)動周期,下同)恢復(fù)了穩(wěn)定,此處由于輪對的大幅擺振出現(xiàn)了輕微的輪緣接觸;隨著車輛運(yùn)行速度的增加,穩(wěn)定的周期解(對應(yīng)于車輛系統(tǒng)穩(wěn)定的蛇形運(yùn)動)幅值也繼續(xù)增加,直到點(diǎn)C(VC=106.619m·s-1,max|yC|=1.732mm,TC=0.6400s),此后系統(tǒng)進(jìn)入擬周期運(yùn)動狀態(tài)(后面進(jìn)行詳細(xì)說明);在另一Hopf分岔點(diǎn)D處,系統(tǒng)仍以亞臨界的方式分岔出一不穩(wěn)定的周期解DE(有1個(gè)實(shí)、1對復(fù)共軛共3個(gè)floquet乘子位于單位圓處),該周期解的幅值隨著速度的減小而逐漸增加直到點(diǎn)E(VE=111.426m·s-1,max|yE|=0.734mm,TE=0.6941s),此時(shí)一個(gè)實(shí)floquet乘子進(jìn)入單位圓內(nèi)。(3)穩(wěn)定極限環(huán)c當(dāng)車輛運(yùn)行速度從小逐漸增加,速度低于VA時(shí),系統(tǒng)最終會穩(wěn)定在相應(yīng)速度對應(yīng)的平衡點(diǎn)上,振幅始終保持為0;當(dāng)速度大于VA時(shí),平衡點(diǎn)失穩(wěn),這時(shí)系統(tǒng)會躍變到穩(wěn)定的極限環(huán)BA′C上,從圖2中的A點(diǎn)躍變到A′點(diǎn),造成系統(tǒng)的“突跳”;當(dāng)速度逐漸降低小于VA時(shí),系統(tǒng)仍會運(yùn)行在穩(wěn)定的極限環(huán)解分支BA′C上,直到速度小于VB時(shí),極限環(huán)消失,系統(tǒng)才會落到平衡點(diǎn)解分支上,從B點(diǎn)突降到B′點(diǎn)。顯然這個(gè)振幅突變BB′比升速時(shí)的突跳AA′延遲了,形成“遲滯”現(xiàn)象。在車輛系統(tǒng)中,這種雙穩(wěn)態(tài)引起的突跳使系統(tǒng)在接近失穩(wěn)點(diǎn)時(shí)由于受到擾動而使振幅急劇變化,會帶來意想不到的問題,因此,應(yīng)該盡量避免。(4)車輛運(yùn)行速度的極限環(huán)從圖2中可以看出,當(dāng)車輛運(yùn)行速度VB<V<VA時(shí),存在穩(wěn)定的定常擺動、不穩(wěn)定的周期擺動與穩(wěn)定的周期擺動共存的現(xiàn)象;當(dāng)車輛運(yùn)行速度VE<V<VD時(shí),則同時(shí)存在3種形式的運(yùn)動:不穩(wěn)定的定常運(yùn)動、不穩(wěn)定的周期運(yùn)動和穩(wěn)定的擬周期運(yùn)動。圖3給出了車輛運(yùn)行速度V為51.5m·s-1時(shí)1位輪對橫向位移-橫向速度的2種極限環(huán),其中不穩(wěn)定的極限環(huán)只能用延續(xù)算法算出,普通的數(shù)值積分方法是無能為力的。因車輛運(yùn)行速度滿足關(guān)系式VB<V<VA,故機(jī)車系統(tǒng)在該速度下存在2種穩(wěn)定的運(yùn)動:一種是穩(wěn)態(tài)運(yùn)動(對應(yīng)于圖3中的定常吸引子點(diǎn)),另一種則是穩(wěn)定的極限環(huán)運(yùn)動(圖3中的實(shí)線吸引極限環(huán)),它們中間被一個(gè)不穩(wěn)定的極限環(huán)(圖3中的點(diǎn)線排斥極限環(huán))隔開。車輛系統(tǒng)最終是趨于穩(wěn)態(tài)運(yùn)動還是趨于穩(wěn)定的極限環(huán)運(yùn)動,主要與系統(tǒng)的初始擾動有關(guān)。如果系統(tǒng)各剛體的擾動均位于其不穩(wěn)定的極限環(huán)內(nèi),則運(yùn)動趨于穩(wěn)定的平衡位置;反之,如果系統(tǒng)的擾動均位于點(diǎn)線以外,則運(yùn)動趨于穩(wěn)定的極限環(huán)運(yùn)動。圖4給出了車輛運(yùn)行速度為107.0m·s-1時(shí)1位輪對橫向位移功率譜和yw1-yw4Poincare截面,以此來說明圖2中速度大于VC后車輛的振動形式。從功率譜可看出譜線是分立的和離散的,不但包含了基頻1.53Hz,而且在這個(gè)頻率周圍還衍生出了另外的頻率1.44Hz和1.63Hz