双参变量下单级齿轮传动系统分岔/冲击特性分析

田亚平， 褚衍东， 饶晓波

(1.兰州交通大学 机电工程学院，兰州 730070； 2. 兰州交通大学 甘肃省轨道交通装备系统动力学与可靠性重点实验室，兰州 730070)

因齿侧间隙、时变啮合刚度、加工制造误差、质量偏心等因素导致齿轮系统表现出强非线性的动力学振动特征。Kahraman等[1]考虑齿侧间隙和误差激励因素建立了单自由度单级齿轮系统动力学模型，发现了次谐和混沌响应。Vinayak等[2]研究了齿轮啮合综合传递误差对系统非线性动力学行为的影响。刘晓宁等[3]采用增量谐波平衡法对三自由度齿轮系统周期运动进行了分析并判稳。苟向锋等[4]采用相图、Lyapnuov指数谱、Poincaré截面法找出了三自由度单级齿轮副随参数变化通向混沌运动的途径。王晓笋等[5]研究了含磨损故障的单级齿轮传动系统的分岔行为及通向混沌的途径。

结合Floquet稳定性理论的PNF法(Poincaré-Newton-Floquet)是周期运动求解及判稳同步进行的一种高效数值方法。Luo等[6]运用PNF法研究了含裂纹以及碰磨耦合故障的转子系统的周期运动的稳定性问题，Han[7-8]运用PNF法研究了碰磨转子系统的周期运动稳定性问题，李同杰等[9]运用PNF法研究了行星齿轮传动系统周期运动的稳定性问题。上述研究中均采用单参变量分析系统的分岔特性，无法全面了解参数间耦合下系统的非线性特性转迁规律。双参平面内研究系统分岔/冲击特性能更清晰的了解系统参数之间的耦合关系并确定系统稳定运行的参数范围，为工程实际中的齿轮结构的参数优化提供理论指导。Gou等[10-11]采用胞映射和逃逸算法研究了双参平面内的纯扭转齿轮副的分岔特性。但采用PNF法和延续算法研究齿轮系统双参平面分岔特性的文献还鲜有报道。

本文以三自由度单级齿轮传动系统的动力学模型为研究对象，采用PNF法和延续追踪法求解系统周期运动并判别分岔类型，然后采用双参平面伪彩图探寻系统的运动分岔/冲击特性随参数转迁和参数间的耦合关系。

1 单级齿轮非线性动力学模型

单级齿轮传动系统动力学模型如图1所示。图中，mgi，Igi，rgi，θgi(i=1,2)分别代表主从动齿轮的质量、转动惯量、基圆半径、扭转角位移；bi，cbi，kbi，Fbi，ygi分别表示主从动齿轮支承轴承的间隙、阻尼、支承刚度和作用于轴承的径向预载荷和径向位移。Tg1，Tg2为输入输出扭矩均值。齿轮啮合的时变刚度、齿侧间隙、啮合线性阻尼、齿轮啮合误差分别用kh(t)，bh，ch，e(t)表示。

图1 单级齿轮传动系统非线性模型

由文献[3]知其系统的量纲一化动力学方程为

(1)

(2)

式中:bhi为轴承支承径向间隙(i=1,2)和半齿侧间隙(i=3)；bi为其对应的量纲一化间隙。

y3为量纲一化齿轮啮合的动态传动误差与静态传动误差的差值。

y3(τ)=(rg1θg1(τ)-rg2θg2(τ)+yg1-

yg2-e(τ))/bc

(3)

(4)

2 系统周期运动分岔转迁规律的求解过程

双参平面内，系统(4)的分岔/冲击特性及其转迁规律采用PNF法周期不动点的求解和周期数的统计、周期解延续算法[12]、分岔类型的判断(Floquet乘子法[13])、齿轮啮合冲击类型的判断来实现。其具体求解过程如下。

2.1 PNF法周期数统计及追踪

PNF法是一种Floquet理论和打靶法相结合的同时进行周期运动求解并判稳的数值法，在周期运动求解中统计出其周期数。将系统方程式(4)改写为周期激励的n+1维状态方程

Rn×R,n=2N,f(X,τ)=f(X,τ+T0),τ>0

(5)

式中：T0=2π/Ω为系统激振力周期；N为系统的自由度，本文N=3。系统的m倍周期运动定义为

X(τ)=X(τ+T), ∀τ>0

(6)

式中：T=mT0，m为正整数。

在n+1维状态空间中定义n维Poincaré截面

∑={(X(τ),τ)|mod(τ,T)=0}

(7)

在该截面上定义Poincaré映射

(8)

连续系统(5)的周期解等价于求解如下代数方程的不动点XF

Q(X)=P(X)-X=0

(9)

式(9)的求解由打靶法思想进行：事先给出XF的一个近似解Xk∈Σ，然后按照Newton-Raphson方法构建迭代公式，以确定更加接近XF的向量Xk+1。

Xk+1=[I-DP(Xk)]-1[P(Xk)-

DP(Xk)·Xk]

(10)

式中：I为单位方阵；P(Xk)为Xk处的Poincaré映射；DP(Xk)为Poincaré映射在Xk处的Jcaobi矩阵(离散转移矩阵)。实际计算P(Xk)的初值由下式迭代

(11)

式(11)在T时刻的解便是Xk的Poincaré映射，即P(Xk)=X(T)。DP(Xk)的n×n个初值求解

(12)

式(12)在T时刻的解就是Poincaré映射在Xk处的Jacobi矩阵，即Φ(T)=DP(Xk)，∂fX(Xk,τ)/∂Xk为状态方程(5)在Xk处对应的Jacobi矩阵。联立式(11)和式(12)得

(13)

式(13)以(Xk,I)为初值积分一个周期T获得系统的Xk+1和DP(Xk)。令Xk=Xk+1作为下次迭代的初始值，通过式(13)反复迭代，直到满足精度要求的不动点XF，就是满足等式(6)的m倍周期运动X(τ)的不动点。

找到稳定周期运动后，根据延续分岔算法引入分岔参数λ，采用延续追踪的方式进行周期追踪。构建代数方程

G(λ,X)=P(λ,X)-X=0

(14)

分岔参数λ的延续追踪过程就是式(14)关于λ的解曲线求解过程。式(14)的解曲线问题转化为求解常微分方程

(15)

的问题。设已求得λ=λk(k=0,1,2,…)时，周期T的Poincaré截面的不动点Xk，采用Euler型数值积分方法，可得λ=λk+1对应不动点Xk+1的初始不动点Xk+1,0的预测公式

(16)

式中：Δλ为步长，GX(λk,Xk)=I-DP(λk，Xk),Gλ(λk,Xk)=Pλ(λk,Xk)。DP(λk，Xk)由式(12)求得，Pλ(λk,Xk)由打靶法思想建立的如下微分方程组求得

(17)

以(Xk，λk)为初值积分一个Poincaré映射周期T获(Xk，Pλ(λk,Xk))。

从分岔参数λk开始，通过式(13)迭代获得周期不动点后，再以式(16)进行延续追踪得λk+Δλ处的初始迭代点，重复式(13)和(16)的迭代过程，可求得一系列外参数λ值的周期不动点Xk及其离散转移矩阵DP(Xk)，完成周期m的追踪。周期解的稳定性由离散转移矩阵DP(Xk)的特征值(Floquet乘子模)的最大值|λ|max确定。当判断周期解失稳时判断分岔类型并重新寻求新的周期数。

2.2 周期和冲击类型的表示方法

定义符号“I/P”为系统的齿轮副啮合冲击类型和周期数，I表示齿轮啮合的冲击类型：I=0,1,2分别表示无冲击、单边冲击、双边冲击类型，P表示周期数。设量纲一化齿侧间隙为2b，依据文献[14]，当xmin>b时齿面接触无冲击；当|xmin|

3 系统双参分岔仿真分析

为了分析系统(4)的动力学分岔、冲击、幅值跳跃等特性，取仿真参数：ζ11=ζ22=0.1，ζ13=ζ23=0.012 5，ζ33=0.05，k11=k22=1.25，fm=0.2，fah1=0.05，fb1=fb2=0.1，b1=b2=0，b3=1，仿真初值取为0。在系统特定参数下，双参平面内用不同的颜色代表不同的“I/P”类型，其区域的边界线即为分岔曲线。Ω-α双参平面分岔图如图2所示。齿轮系统常稳定运行在短周期状态下，为简化分岔图对P≥65的长周期、拟周期(Qusi-P)和混沌(chaos)运动均用“I/N”表示(白色区域)。当系统从状态“0/P”向“1/P”状态转迁中，系统从无冲击区域在x=b截面处穿越到脱啮区，在该点发生了擦切(G Bif)，其交线为擦切分岔曲线。当系统从“I/P”向“I/2P”转迁时发生了倍化分岔(PD Bif)，其交线为倍化分岔曲线。当系统从“I/1”向“I/n”转迁时系统从单周期经hopf分岔(PF Bif)为拟周期运动或经激变分岔(CIC)为混沌运动，具体由Floquet乘子判断。

齿轮系统的重合度和转速是影响系统运转稳定性的主要因素，而时变啮合刚度幅值系数α是表征重合度大小的主要指标，量纲一频率Ω是表征转速大小的主要指标。因此，选双参(Ω，α)进行系统分岔特性转迁规律分析，如图2所示。横坐标量纲一化频率Ω∈[1,2]，纵坐标啮合刚度幅值系数α∈[0,0.49]。在双参平面内，系统出现了倍化(PD Bif)、擦切(G Bif)、幅值跳跃(JP Bif)、激变、hopf等分岔现象，周期1、周期2、周期4、周期8、周期16、拟周期、混沌等运动状态和三种冲击状态共存。

图2 Ω-α双参平面分岔图

沿Ω方向系统分岔转迁规律较复杂。当α∈[0,0.05]范围内(小幅值刚度波动大重合度状态下)系统表现出擦切和倍化分岔，系统在“0/1”、“1/1”和“1/2”运动状态转迁出现周期泡现象，并在Ω∈[1.4,1.6]区域内系统处于单边冲击脱啮状态。当α∈(0.05,0.1]范围内系统倍化分岔出现了“1/4”、“1/8”运动状态。当α∈(0.1,0.4]范围内系统通过hopf和鞍结激变分岔，单周期运动进入混沌运动后经逆倍化分岔经周期“1/16”、“1/8”、…、到稳定“0/1”运动状态，单边冲击脱啮的范围逐渐扩大。当α∈(0.4,0.49]范围内系统经历了单周期、激变分岔进入双边冲击齿背接触混沌“2/N”状态、混沌吸引子激变演变为单边冲击脱啮混沌“1/N”运动、经短暂的周期运动后Hopf分岔进入拟周期运动、逆Hopf分岔进入单周期运动状态。

(a) Ω-α-x5

(b) α-Ω-x5

沿α方向系统分岔运动的转迁规律较为简单。当Ω∈[1,1.35]∪[1.58,2]范围内系统处于单周期无冲击稳定运动状态，仅在Ω∈[1.6,1.9]区域内出现了擦切分岔形成了单边冲击周期运动状态。在Ω∈(1.35,1.58)的共振频率附近，系统通过倍周期分岔形成周期2、周期4、周期8、…、混沌运动的演化过程。特别在α∈(0.22,0.3)范围内发生小窗口的多次幅值跳跃双边冲击周期2运动状态，即在该区域内系统的系统动力特性非常复杂，设计中应尽量避开该区域选取参数。

为了验证PNF法和延续追踪法求解的双参平面分岔图(图2)的正确性，沿图2中的α和Ω方向等间距选取截面采用4～5阶Runge-kutta数值法仿真获得系统(4)的Ω-α-x5三维分岔图如图3所示。图中，齿轮啮合相对综合误差x5随Ω和α的变化，系统的动力学分岔特性的转迁规律及分岔点位置与图2相一致，说明PNF法和延续追踪法求解双参平面分岔图的方法是有效的。

取图2中Ω=1.45截面的α-x5分岔图如图4(a)所示。系统通过倍化分岔进入混沌运动，其间出现了幅值跳跃分岔现象。取图2中α=0.45截面的Ω-x5分岔图如图4(b)所示。系统由混沌运动激变退化为周期8运动、再经过鞍结分岔幅值跳跃为周期2运动，经hopf分岔进入拟周期运动，逆hopf分岔退化为周期2运动，激变分岔为混沌运动、逆倍化分岔退化为周期1运动的转迁规律。

取图4(b)中Ω=1.5截面关于(x5,x6)Poincaré映射图如图5(a)所示。图5(a)中出现近似周期2的黑点，将其放大可看见两个不稳定的焦点向外扩散，当Ω增大到1.505时系统焦点分岔为两个封闭的hopf圈(图5(b)所示)，即在Ω=1.5处发生了hopf分岔，考查其Floquet乘子为[0.791 1, -0.889 2±0.459 5i, -0.624 2±0.553 1i, -0.310 4], |λ|max=|-0.889 2±0.459 5i|=1.049，即从一对共轭复数穿出单位圆。在图2中取α=0.079，Ω=1.632、1.638、1.642值的(x5,x6)相图如图5(c)、图5(d)、图5(e)所示，相图展示了系统从擦切前的单边冲击运动向擦切后的无冲击运动的转迁过程，很明显在x5=1截面处发生了擦切运动。

采用Runge-kutta数值法仿真的三维/二维分岔图、Poincaré映射图和相图验证了综合PNF法和延续追踪法的双参平面分岔求解的正确性。因延续追踪法省去了大量的瞬态过程求解时间，故是一种有效的双参平面分岔计算方法。

(a) Ω=1.45

(b) α=0.45

(a) Ω=1.50, α=0.45

(b) Ω=1.505, α=0.45

(c) Ω=1.634, α=0.079

(d) Ω=1.638, α=0.079

(e) Ω=1.642, α=0.079

双参平面分岔图揭示了参数Ω、α范围内系统分岔特性的转迁规律，获得了周期运动、拟周期运动、混沌运动的分岔曲线。工程实际应用中，为了延长齿轮设备的运行寿命降低事故的发生的几率和设备维修费用，常常希望齿轮处于“0/1”运动状态，因此在齿轮装备设计中合理选择齿轮结构参数、润滑油黏度及制造精度使系统量纲一化参数在双参平面 “0/1”区域内为佳。

4 结 论

(1) 在双参分岔平面内，提出了一种PNF法和延续算法相结合的非线性动力系统周期运动求解和啮合冲击类型判别(I/P分岔)的新方法，并用Runge-Kutta数值法的相图、三维/二维分岔图等验证了其有效性。

(2) 以三自由度单级齿轮传动系统为分析模型，在时变啮合刚度幅值系数和量纲一化频率双参平面内分析了系统分岔/冲击特性及转迁规律，划分了双参耦合下系统稳定周期运动和失稳混沌运动参数区域。

(3) 在时变刚度幅值系数和量纲一化频率特定参数耦合状态下，避开共振频率范围和小时变啮合刚度幅值系数(大重合度)范围内选取参数能提高齿轮传动的稳定性降低疲劳磨损。