基于边缘接触时变刚度的轮齿表面剥落动力学模型与响应特征

邵毅敏, 王新龙, 刘 静, 陈再刚

(重庆大学 机械传动国家重点实验室，重庆 400030)

啮合接触刚度是齿轮系统中影响系统动态响应的重要因素，目前的研究集中在轮齿故障的动态特性计算方面[1-3]，针对齿面剥落轮齿在啮出时，剥落边界和对应轮齿齿面发生边缘接触作用，以及这种边缘接触作用对啮合刚度及其动态响应的影响和能量传递的影响等研究较少。

本文考虑边缘接触作用，建立了“复合刚度计算模型”，提出了轮齿表面剥落缺陷的时变刚度算法，计算了剥落扩展过程中的不同剥落界面形貌的边缘接触时线接触弹性接触刚度和啮合刚度，通过对时变信号和频谱的分析，研究了存在表面剥落的齿轮系统动态响应特征。获得了表面剥落对于齿轮系统动态响应的影响特性。

1 轮齿表面剥落缺陷的时变刚度算法

表面剥落缺陷轮齿与正常轮齿啮合过程中，在啮出剥落边界时，与其边缘发生接触作用，导致齿面接触力产生突变，直接影响齿对啮合动态特性，因此，分析边缘接触问题，建立考虑边缘接触的啮合刚度计算模型，是研究表面剥落对啮合刚度及其动态响应影响的重要问题。

本文提出的复合刚度计算模型是同时包含基于势能原理的啮合刚度计算模型与边缘接触刚度计算模型的综合啮合刚度计算模型。

1.1 轮齿表面剥落的复合接触刚度模型

图1 表面剥落齿对啮合示意图

表面剥落齿对啮合，如图1所示。剥落区域按照啮合线法向和切向投影得到齿对啮合过程中剥落区域的接触过程示意图，如图2和图3所示。其虚线表示发生边缘接触位置。表面剥落齿对啮合过程接触刚度情况，如图4所示，Kn、Ks分别表示沿齿宽方向的渐开线接触刚度、边缘接触刚度。

(1) 在进入剥落边界之前的正常齿面区域，齿对的啮合通过渐开线齿面线接触啮合，接触刚度包括沿齿面方向的Kn(如图4(a))；

(2) 在剥落且非边界区域(如图2(b), 啮合点Y)，齿面与剥落轮齿两端齿面通过渐开线接触啮合，接触刚度包括两端的Kn(如图4(b))；

(3) 在啮出剥落边界时(如图2(c), 啮合点Z)，齿面与剥落边界边缘发生接触，导致接触力产生突变，同时齿面与非剥落区域产生渐开线接触，接触刚度包括Kn、Ks(如图4(c))；

(4) 齿对啮合离开剥落区域，接触情况和(1)相同(如图4(d))。因此表面剥落齿对啮合的接触过程存在一边-两边-三边-一边的变化过程，所以，啮合刚度K(t)的计算必须要用分段函数表示。

图2 剥落区域齿对啮合法向示意图

图3 剥落区域齿对啮合切向示意图

图4 表面剥落齿对啮合过程接触刚度示意图

齿对啮合刚度[6]可以表示为式(1)，其中Kh、Kb、Ks、Ka、Kf分别表示接触刚度、弯曲刚度、剪切刚度、轴向压缩刚度、基体变形刚度。

K(t)=

(1)

式中:Rs表示剥落边界渐开线曲率半径，Rs在剥落扩展的过程中保持常量特性。利用线接触弹性接触变形的解析算法[4,5]，计算边缘接触刚度，计算有限长线接触相对应的计算参数ε:

(π/2-θ1)/ε=(π/2-arccosε)/ε=1

(2)

与接触区域相关的面积系数的计算：

Δ1=(f0-f1)(π/2-θ1)/2

(3)

Δi=(fi-1-fi+1)(θi-θi+1)/2

(4)

(5)

离散积分：

(6)

光滑弹性体接触理论中的法向接触变形计算式为：

δs=ηbp0Kθ1

(7)

式中：η表示和接触体材料属性有关的特性系数：

(8)

Kθ1表示第一类全椭圆积分, 其计算式为：

(9)

由式(6)和式(9)，可获得Kθ1的简化计算方程：

(10)

p0表示接触区最大接触压力, 计算式为：

p0=3P/(2πLab)

(11)

其中：

(12)

(13)

∑ρ表示主曲率之和：

(14)

渐开线齿对啮合的接触刚度沿整个接触线上可视为近似常数[9]。渐开线齿对啮合接触刚度的计算采用式(15)，由式(15)和式(7)可得复合接触刚度式(16):

(15)

(16)

则齿对在整个啮合过程的复合接触刚度Kh计算式表示为如下分段函数：

(17)

式中:α1、αI、αO、LsH、W分别表示任意啮合压力角、啮入剥落边界压力角、啮出剥落边界压力角、非剥落轮齿长度、轮齿宽度。

将轮齿视为变截面悬臂梁计算直齿轮轮齿啮合过程的Kb、Ks、Ka、Kf，则齿对啮合的弯曲、剪切和轴向压缩刚度计算式可以表示为[2,6]：

(18)

(19)

齿轮基体变形对于时变啮合刚度有重要影响，Sainsot等[7]计算得出基体变形刚度计算式：

(20)

(21)

综合上述，将复合接触刚度Kh及Kb、Ks、Ka、Kf计算式带入啮合刚度公式(17)，并设：

(22)

则整个啮合过程的复合刚度计算式表示为如下分段函数K(t):

(23)

文献[6, 8]分片积分计算法，将剥落模型分成正常、故障两部分分别计算变形刚度，整个轮齿的刚度可以将啮合刚度在整个轮齿齿宽W积分获得：

(24)

1.2 接触刚度计算模型验证

为了验证本文边缘接触刚度计算模型的正确性，建立了相应的有限元分析模型。利用有限元法和本文提出的边缘接触刚度算法分别计算了存在齿面剥落故障的齿轮啮合副边缘接触刚度，其结果如图5所示。图5显示了四种不同剥落宽度的8组边缘接触刚度计算结果。两者计算结果一致，最大误差仅为4.54%，证明了本文提出的边缘接触刚度计算模型的正确性。

图5 复合接触刚度算法结果与有限元计算结果

2 仿真计算条件及仿真结果

2.1 仿真条件

仿真的系统参数，如表1所示。仿真分析轮齿表面剥落缺陷位于小齿轮，表面剥落参数，如表2所示。

2.2 仿真结果

不同剥落宽度(D=45 mm,H=1.0 mm)影响的单双齿对啮合刚度曲线、不同剥落长度(A=3 mm,H=1.0 mm)影响的单双齿对啮合刚度曲线，分别如图6和图7、图8和图9所示。横坐标为啮合压力角，实线方块、虚线方格分别表示考虑边缘接触“与”“否”的计算结果，计算结果显示：随着剥落宽度扩展，啮合刚度曲线出现明显衰减现象，而最大的衰减出现在存在剥落缺陷的轮齿退出啮合的位置，特别是齿对在啮出剥落边界时，由于边缘接触作用的影响，导致啮合刚度发生突变现象。

表1 齿轮系统参数

表2 小轮齿表面剥落缺陷仿真参数

图6 不同剥落宽度影响的单齿对啮合刚度曲线

图7 不同剥落宽度影响的综合啮合刚度曲线

图8 不同剥落长度影响的单齿对啮合刚度曲线

图9 不同剥落长度影响的综合啮合刚度曲线

3 轮齿表面剥落缺陷的动力学响应及特征和实验研究

3.1 齿轮系统的动力学模型

考虑边缘接触的刚度激励，建立了6自由度考虑阻尼系数的直齿轮箱系统[6]，动力学模型示意图如图10所示。Y轴平行于轮齿啮合线，TP/TG表示驱动/制动扭矩，ΩP/ΩG表示小齿轮/大齿轮的工作转速，JP/JG表示小齿轮/大齿轮的转动惯量，KiBj/CIBj表示支撑轴承的刚度和阻尼，i分别表示小齿轮和大齿轮，j分别表示x和y坐标轴方向，mP/mG表示小齿轮和大齿轮质量，K(t)表示时变啮合刚度，Cm表示轮齿啮合阻尼，Ff表示存在于啮合轮齿间的摩擦力。

直齿轮系统的参数，如表1所示。施加到直齿轮系统的负载扭矩为60 Nm，仿真转动频率为10 Hz。采用Coulomp摩擦模型，其动力学摩擦系数u参考文献[6]。

建立系统振动方程式为：

Mx″+Cx′+Kx=F

(25)

(26)

(27)

(28)

(29)

图10 齿轮系统动力学模型

3.2 轮齿表面剥落缺陷对系统动态特性影响分析

采用龙格-库塔法求解时变刚度激励的动力学方程，获得系统动态响应。

仿真计算直齿轮系统的小齿轮y向动态响应，其正常动态响应时域图和频谱图，如图11所示，其时域信号平稳，有啮合频率和谐频成分。

图11 正常小齿轮y向动态响应时域图及频谱图(D=0 mm,A=0 mm)

稳定状态的剥落缺陷影响的时域图和频谱分别如图12 (D=45 mm,A=2 mm)、图13 (D=45 mm,A=3 mm) 和图14(D=60 mm,A=2 mm)所示。红色线、蓝色线分别表示考虑边缘接触和不考虑边缘接触的响应。时域分析显示初期剥落轮齿产生的时域信号特性并不明显。然而，当剥落宽度和长度扩展到一定程度，出现明显的冲击，相邻冲击间的间隔时间T(0.1 s)和小齿轮转动周期(60/600=0.1 s)相同。

包络解调频谱显示剥落缺陷宽度和长度扩展，影响啮合频率及其谐频以及边频的变化，其特点是：边频明显增强，而齿轮啮合频率及其谐频幅值增长相对缓慢，边频特征比齿轮啮合频率及其谐频对于剥落缺陷扩展更加敏感，其结果与文献[1]及文献[3]的有限元计算结论一致，从另一个角度验证了本研究方法的正确性。

3.3 轮齿表面剥落缺陷实验与特性研究

为验证理论模型仿真分析的正确性，利用二级齿轮箱开展轮齿表面剥落缺陷实验研究。实验条件如图15所示：输入轴转速600 r/min(即转动频率为10 Hz)，采样频率为10 kHz；齿面剥落缺陷设置在第二级小齿轮轮齿上，测点布置在第二级小齿轮轴承座上；齿面剥落：小缺陷(D=45 mm，A=2 mm,H=1 mm)和大缺陷(D=60 mm，A=2 mm,H=1 mm)。

齿面剥落小缺陷与大缺陷对应的振动信号时域图和频谱分别如图15、图16所示。结果表明，齿面剥落缺陷较小时，齿轮系统振动冲击信号特征不明显。齿面剥落缺陷越大，齿轮系统振动冲击越明显，如图16所示，相邻振动冲击间的间隔时间T(0.101 7 s)与小齿轮转动周期(1/9.83 s)相同。

图12 小齿轮y向动态响应时域图、频谱图(D=45 mm,A=2 mm)

图15 轮齿表面剥落故障实验布置图

图16 小齿轮y向实验动态响应时域图、频谱图(D=45 mm,A=2 mm,H=1 mm)

图17 小齿轮y向实验动态响应时域图、频谱图(D=60 mm,A=2 mm,H=1 mm)

频谱中f1表示第一级啮合频率(230 Hz)，f2表示第二级(缺陷齿轮)啮合频率(147.5 Hz)，含故障的齿轮轴转频为5.9 Hz；频谱显示，随着剥落扩展，边频(与缺陷齿轮转频倍数吻合)明显增强，而缺陷齿轮啮合频率及其谐频幅值增长相对缓慢，边频特性相比齿轮啮合频率及其谐频对剥落缺陷扩展更加敏感。实验结果和文中3.2节的仿真结果吻合良好，验证了建立的齿轮剥落缺陷动力学理论模型的正确性。

4 结 论

提出了齿面剥落复合接触刚度模型，结合势能原理推导了齿轮副啮合刚度算法，分析了边缘接触对时变刚度激励以及系统动态响应特性的影响，开展了轮齿表面剥落缺陷实验研究。仿真与实验结果表明：缺陷边缘接触导致啮合刚度发生突变，其动态响应信号的啮合频率及其倍频的调制边频成分明显。本文提出的复合刚度模型能够准确表征剥落缺陷啮合的动态响应特征，可为在线监测和齿轮故障诊断提供理论支撑。

[1]Chaari F, Baccar W, Abbes M S, Mohamed Haddar, Effect of spalling or tooth breakage on gearmesh stiffness and dynamic response of a one-stage spur gear transmission, European Journal of Mechanics A/Solids,2008,27:691-705.

[2]Wu S Y, Zuo M J, Parey A, Simulation of spur gear dynamics and estimation of fault growth[J].Journal of Sound and Vibration, 2008,317:608-624.

[3]Jia S X, Howard I, Comparison of localised spalling and crack damage from dynamic modelling of spur gear vibrations[J]. Mechanical Systems and Signal Processing,2006,20:332-349.

[4]Xiao H F, Brennan M J, Shao Y M，On the undamped free vibration of a mass interacting with a Hertzian contact stiffness[J].Mechanics Research Communications，2011,38:560-564.

[5]Ding Chang-an, Zhang Lei, Zhou Fu-zhang, et al. Theoretical formula for calculation of line-contact Elastic Contact deformation[J]. Tribology，2001.

[6]Chen Z G, Shao Y M, Dynamic simulation of spur gear with tooth root crack propagating along tooth width and crack depth[J]. Engineering Failure Analysis, 2011.07.006.

[7]Sainsot P, Velex P, Contribution of gear body to tooth deflections-A new bidimensional analytical formula[J]. ASME J. Mec. Des.,2004,126:748-752.SHOCK,2012,31:83-89.

[8]Chen Z G, Shao Y M, Mesh stiffness calculation of a spur gear pair with tooth profile modification and tooth root crack[J]. Mechanism and Machine Theory, 2013,62:63-74.

[9]Yang D C H, Su Z S, A rotary model for spur gear dynamics[J]. ASME J. Mech., Trans. Aut. Des.,1985,107:529-535.