基于蠕滑理论的蛇形运动分岔问题研究

基于蠕滑理论的蛇形运动分岔问题研究

0蛇形运动稳定性研究当铁路车辆沿轨迹移动时，具有一定形状的转向车轮对角（水平移动），旋转至铜轴的棕色轴（旋转）。这两种运动的结合使车辆沿着正确的轨道前进。这种现象通常被称为蛇形。蛇形是一种非常普遍的自激振横向摆动不稳定现象,主要是由轮轨之间的非线性接触力及悬挂装置作用而引起。高速车辆蛇形失稳不仅恶化车辆运行性能,降低乘坐舒适性,而且还会加剧部件磨损,对车辆及线路造成损伤,甚至会引发脱轨的重大事故,因此,车辆系统蛇形运动稳定性的研究具有重要意义。国内外学者对车辆系统横向稳定性及分岔问题进行了大量研究:Wickens推导了单个轮对的运动方程,并从动力学角度详细分析了轮对的失稳问题;Kaas-Petersen分析了Cooperrider构架在有轮缘和无轮缘情况下的对称摆振、不对称摆振以及混沌行为;True在Kaas-Petersen研究的基础对稳定性及分岔问题提出了很多好的方法和建议;曾京使用QR算法与黄金分割法确定车辆系统的Hopf分岔点,用打靶法进行极限环的数值求解,较好地解决了相关的数值计算问题;Ahmadian等用渐近法研究了单个转向架的稳定性,并对诸多影响因素进行了分析;司道林等应用多刚体动力学软件SIMPACK和数值积分技术,得到了牵引力对车辆临界速度的影响规律。在这些研究中,对简单模型研究比较深入,对复杂铁道车辆系统的研究往往侧重Hopf分岔的参数条件,对Hopf分岔后周期解的变化过程研究较少。基于此,本文对一六轴机车系统进行研究,运用延续算法对Hopf分岔及分岔后极限环的解分支进行连续跟踪,讨论Hopf分岔的类型和分岔过程中出现的多种非线性动力学现象,并阐述其中的工程意义。1力学模型1.1横向运动方程图1中的车辆模型包含车体(Mc,Icx,Icz)、2个三轴转向架(Mt,Itx,Itz)、6个轮对(Mw,Iwz)和一系悬挂(Kpx,Kpy,Kpz,Cpx,Cpy,Cpz)和二系悬挂(Ksx,Ksy,Ksz,Csx,Csy,Csz)的六轴机车系统,暂不考虑牵引电机的振动,一系悬挂、二系悬挂均视为线性弹簧和阻尼的组合。车辆系统不考虑轨道不平顺带来的影响,假设轮子在足够光滑、水平的直线钢轨上滚动。由于机车系统前后、左右对称,可以近似地把横向运动方程和垂向运动方程解耦,得到独立的横向运动方程组。整个车辆系统共有21个自由度:车体考虑三个自由度,即横移yc、侧滚ϕc和摇头ψc运动;每个转向架考虑构架的横移ytj(j=1,2)、侧滚ϕtj及摇头ψtj运动;对于轮对则只考虑横移ywi(i=1,…,6)和摇头ψwi运动。应用牛顿定律可导出整个系统的运动微分方程,用如下形式的一阶非线性常微分方程组来表示dzdt=f(z,V)(1)dzdt=f(z,V)(1)式中:z∈R42为包含系统各刚体位移与速度的状态向量;V∈R+为车辆运行速度;t为时间。1.2轮轨接触力与轮缘力轮轨的相互作用关系主要包括轮轨接触几何关系和非线性轮轨接触力。传统的机车车轮踏面都设计成锥面,以减小通过曲线时的磨耗和有助于自导向,轮轨接触几何参数,如左、右轮滚动圆半径(rL,rR),轮轨接触角(δL,δR)以及轮对侧滚角位移(φw)都可近似认为是轮对横移量yw的非线性函数{rL=r0+κywrR=r0-κywδL=δR=κφw=κyw/a(2)⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪rL=r0+κywrR=r0−κywδL=δR=κφw=κyw/a(2)式中:r0为车轮名义滚动圆半径;a为轮轨接触点横向距离之半;κ为接触角。非线性轮轨接触力主要包括轮轨滚动接触蠕滑力和轮缘力。本文采用Vermeulen-Johnson理论计算轮轨滚动接触蠕滑力,并充分考虑蠕滑力的饱和效应(但不考虑自旋蠕滑力)。轮对的纵向和横向蠕滑率分别为{ξ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公式中计算出的权重系数。纵向和横向蠕滑力分别为{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⋅{u−u2/3+u3/271u<3u≥3(6)u=ξRGπaebe/(μN)(7)式中:N为轮轨接触斑的法向力;G为轮轨合成剪切模量;ae、be分别为轮轨接触椭圆的长、短半轴长度;μ为轮轨间的库仑摩擦因数;FR为修正后的合成蠕滑力;u为经归一化修正后的合成蠕滑率。轮缘力一般用有死区的刚性弹簧来模拟,即可表达为一分段线性函数形式FΤ={k0(yw-η)yw>η0|yw|≤ηk0(yw+η)yw<-η(8)FT=⎧⎩⎨⎪⎪k0(yw−η)0k0(yw+η)yw>η|yw|≤ηyw<−η(8)式中:k0为弹性系数;η为轮缘间隙,不同轮对的轮缘间隙大小可能不一样。2基于生物活性的新能源更新算法一般来说,在一定速度下,由式(1)确定的车辆系统存在稳态摆动和周期摆动,对应于微分方程组的定常解和周期解,因此,下面讨论如何求解式(1)的定常解和周期解。定常解与周期解分支曲线都属于解的分岔问题。在这里,不失一般性,假设x为n维状态向量,λ为一般参数,定义映射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映射的函数。换句话说,现在只需关心式(9)的求解问题了。延续算法是从系统的一个已知解出发,通过预测-校正的机制可逐步求解出整个系统的解分支曲线。假如(x0,λ0)是解分支曲线上的已知点,为了找到λ=λ0+Δλ(Δλ是足够小步长)处解分支上的点,对F(x,λ)=0关于λ进行微分、变换,可得(dxdλ)0=-[■SymbolΡropBΤvBpF(x0,λ0)■SymbolΡropBΤvBpx]-1■SymbolΡropBΤvBpF(x0,λ0)■SymbolΡropBΤvBpλ(12)故λ0+Δλ处的x预测值为xp=x0+Δλ(dxdλ)0(13)以这个值作为初始值,利用Newton-Raphson迭代法求解式(9),可得到λ0+Δλ处x的计算值。如果ue851SymbolPropBTvBpF/ue851SymbolPropBTvBpx奇异,或者遇到解曲线的分岔点时,上面的算法可能会失效,这时可考虑如下形式的扩展系统{G(y,s)=[F(x,λ)n(x,λ,s)]=0y=(x,λ)(14)式中:y与s分别为扩展系统的状态向量与参数;n(x,λ,s)为函数。这样处理后即便ue851SymbolPropBTvBpF/ue851SymbolPropBTvBpx奇异,通过选择合适的n(x,λ,s)可使G(y0,s0)非奇异,因此y=(x,λ)在(y0,s0)附近仍是s的光滑函数,预测校正算法仍可以应用于以(y0,s0)为初始点G(y,s)=0的系统,只不过此时预测值表达式变为{xp=x0+Δs(dxds)0λp=λ0+Δs(dλds)0(15)s为通常的弧长参数,该方法因此也称为弧长延续算法。当然,n(x,λ,s)的选择并不唯一,Doedel使用的n(x,λ,s)形式为n(x,λ,s)=(dxds)Τ0(x-x0)+(dλds)0(λ-λ0)-(s-s0)=0(16)以(x0,λ0)为初始点,利用式(15)进行预测,然后求解式(14),可得到解分支上的下一点,不断重复此过程,可以得到系统的解分支曲线。对复杂的多自由度多强非线性的铁道车辆系统而言,用理论方法确定式(11)的Poincare映射函数是非常困难的,因此,大多数情况下需要用数值方法得到系统的Poincare映射,本文取前转向架1位轮对与钢轨碰撞后横向运动速度为零且横向位移为正的面为Poincare截面σ={(z,V)∈R42×R+,˙yw1=0,yw1>0}(17)通过数值积分来建立系统的Poincare映射函数。3车辆运行速度的数值取值本文以戚墅堰机车车辆厂1990年生产的东风9型(DF9)内燃机车为分析对象,机车系统各刚体的质量、刚度系数、阻尼系数等参数的取值见文献,其中抗蛇形减震器卸荷速度取0.043m·s-1,减震器饱和阻力为13kN。轮轨计算参数在文献中均可查到,车辆运行速度变化范围为0<V≤150m·s-1。图2是应用延续算法得到的以车辆运行速度作为控制参数的1位轮对相对轨道的最大横向位移分岔图。点线代表不稳定的运动,具体分析如下。(1)系统jacabi矩阵特征值实部点A(VA=53.700m·s-1,α1,2=-1.751×10-10±5.311i,αi代表系统特征值,下同)是系统的第1个Hopf分岔点,此时有一对复共轭特征值正向穿越虚轴,定常解因此而失去其原有的稳定性;随着速度的增大,系统Jacobi矩阵特征值实部基本是增加的;当速度达到VD(VD=114.657m·s-1,α1,2=1.585±8.903i,α3,4=3.520×10-11±9.266i)时,又有一对复共轭特征值正向穿越虚轴而使定常解变得“更加不稳定”;随着速度的继续增大直到速度终止值,再没有发现其他的Hopf分岔点。(2)系统的周期运动状态分别以Hopf分岔点A和点D为初始点进行周期解的延拓计算,发现系统从点A处以亚临界的方式分岔出一不稳定的周期解AB(有1个floquet乘子位于单位圆外,对应于车辆系统不稳定的蛇形运动),分岔出的不稳定周期解在鞍结分岔点B(VB=50.855m·s-1,max|yB|=1.500mm,TB=1.2345s,yi为点i的横移,Ti代表系统运动周期,下同)恢复了稳定,此处由于轮对的大幅摆振出现了轻微的轮缘接触;随着车辆运行速度的增加,稳定的周期解(对应于车辆系统稳定的蛇形运动)幅值也继续增加,直到点C(VC=106.619m·s-1,max|yC|=1.732mm,TC=0.6400s),此后系统进入拟周期运动状态(后面进行详细说明);在另一Hopf分岔点D处,系统仍以亚临界的方式分岔出一不稳定的周期解DE(有1个实、1对复共轭共3个floquet乘子位于单位圆处),该周期解的幅值随着速度的减小而逐渐增加直到点E(VE=111.426m·s-1,max|yE|=0.734mm,TE=0.6941s),此时一个实floquet乘子进入单位圆内。(3)稳定极限环c当车辆运行速度从小逐渐增加,速度低于VA时,系统最终会稳定在相应速度对应的平衡点上,振幅始终保持为0;当速度大于VA时,平衡点失稳,这时系统会跃变到稳定的极限环BA′C上,从图2中的A点跃变到A′点,造成系统的“突跳”;当速度逐渐降低小于VA时,系统仍会运行在稳定的极限环解分支BA′C上,直到速度小于VB时,极限环消失,系统才会落到平衡点解分支上,从B点突降到B′点。显然这个振幅突变BB′比升速时的突跳AA′延迟了,形成“迟滞”现象。在车辆系统中,这种双稳态引起的突跳使系统在接近失稳点时由于受到扰动而使振幅急剧变化,会带来意想不到的问题,因此,应该尽量避免。(4)车辆运行速度的极限环从图2中可以看出,当车辆运行速度VB<V<VA时,存在稳定的定常摆动、不稳定的周期摆动与稳定的周期摆动共存的现象;当车辆运行速度VE<V<VD时,则同时存在3种形式的运动:不稳定的定常运动、不稳定的周期运动和稳定的拟周期运动。图3给出了车辆运行速度V为51.5m·s-1时1位轮对横向位移-横向速度的2种极限环,其中不稳定的极限环只能用延续算法算出,普通的数值积分方法是无能为力的。因车辆运行速度满足关系式VB<V<VA,故机车系统在该速度下存在2种稳定的运动:一种是稳态运动(对应于图3中的定常吸引子点),另一种则是稳定的极限环运动(图3中的实线吸引极限环),它们中间被一个不稳定的极限环(图3中的点线排斥极限环)隔开。车辆系统最终是趋于稳态运动还是趋于稳定的极限环运动,主要与系统的初始扰动有关。如果系统各刚体的扰动均位于其不稳定的极限环内,则运动趋于稳定的平衡位置;反之,如果系统的扰动均位于点线以外,则运动趋于稳定的极限环运动。图4给出了车辆运行速度为107.0m·s-1时1位轮对横向位移功率谱和yw1-yw4Poincare截面,以此来说明图2中速度大于VC后车辆的振动形式。从功率谱可看出谱线是分立的和离散的,不但包含了基频1.53Hz,而且在这个频率周围还衍生出了另外的频率1.44Hz和1.63Hz