一种弧齿锥齿轮时变啮合刚度和传动误差半解析计算方法

0 引言

在齿轮传动系统中，齿轮副的啮合刚度和传动误差是引起传动系统振动和噪声的主要原因[1-2]。精确计算弧齿锥齿轮时变啮合刚度与传动误差是其动态性能研究的基础。

弧齿锥齿轮啮合刚度和传动误差计算方面，国内外学者做了大量研究。Gosselin C 等[3]采用有限条法计算了弧齿锥齿轮轮齿刚度，但是只计算了一定点在载荷作用下弧齿锥齿轮轮齿的变形量，且没有研究啮合周期内啮合刚度时变特性。Mennem[4]使用有限元法求解轮齿接触柔度，计算得到了不同载荷下的时变啮合刚度。CHENG Y 等[5]采用准静态分析方法对准双曲面齿轮进行了轮齿加载接触分析（LTCA），并进一步得到了齿轮副啮合刚度和传动误差。曹雪梅等[6]基于Westinghouse 悬臂梁计算公式计算了单齿对啮合刚度，研究了接触轨迹方向和几何传动误差对承载传动误差的影响。蔡守宇等[7]基于线性规划法与有限元法提出了弧齿锥齿轮扭转啮合刚度计算方法。王延忠等[8]针对高速重载螺旋锥齿轮，在轮齿接触分析（TCA）和LTCA 的基础上提出了时变啮合刚度计算数值模型，研究了不同设计加工参数情况下时变啮合刚度和传动误差的变化规律。唐进元等[9]提出了基于有限元方法的螺旋锥齿轮时变啮合刚度计算模型，运用有限元静态分析提取螺旋锥齿轮法向接触合力和综合弹性变形量，求得了单齿啮合刚度，研究了不同负载转矩对综合啮合刚度的影响。汪中厚等[10]提出了螺旋锥齿轮动传动误差分析方法，对比了动力学分析与静力学分析计算的传动误差不同，分析了载荷对传动误差的影响。周驰等[11]采用Abaqus 隐式算法模拟准双曲面齿轮的准静态啮合过程，计算求得了齿轮副时变啮合刚度和传动误差等时变啮合参数，并研究了负载转矩大小对齿轮副时变啮合参数的影响。

关于弧齿锥齿轮副时变啮合刚度和传动误差研究，目前大多学者采用有限元静态分析方法，1个啮合周期内各个啮合转角都要进行调整三维模型、有限元网格划分及求解分析等过程，计算多次，且计算时间周期长；此外，大多学者采用单齿面法向接触合力与节点弹性变形量的平均值之比计算单齿啮合刚度，但椭圆接触区域内各个节点弹性变形量并不相同，使用节点弹性变形量平均值代替各个节点的平均值必定会产生计算误差，影响时变啮合刚度计算精度。本文中改进了弧齿锥齿轮副时变啮合刚度计算方法，采用单个节点啮合刚度叠加求得单齿啮合刚度，且基于显式动态分析方法对齿轮副动态啮合传动过程进行仿真，只需进行1次有限元分析计算，即可求得1个啮合周期内的时变啮合刚度和传动误差。并进一步研究了弧齿锥齿轮接触椭圆长轴长度、接触轨迹方向两个接触参数及不同负载转矩工况对啮合刚度和传动误差的影响规律，为弧齿锥齿轮系统动力学研究提供依据。

1 啮合刚度/传动误差计算模型

1.1 时变啮合刚度

目前研究中，在计算单齿啮合刚度kn时，大多采用齿面法向接触合力与弹性变形量平均值之比计算，即

式中，i 为工作齿面任意节点；n 为接触椭圆区域共节点数量；fi 为作用于节点i 的法向接触力；ui 为节点i的综合弹性变形。

由于弧齿锥齿轮接触椭圆区域各个节点的变形均不相同，式（1）中采用弹性变形量平均值计算的单齿啮合刚度存在较大误差，如图1所示。本文中引入单个节点啮合刚度，任意节点i啮合刚度ki表达式为

图1 接触椭圆区域法向接触力和弹性变形量

有学者[12-13]将齿面接触区域视为一条接触线，在计算齿轮副单齿啮合刚度时只考虑接触线上节点的变形；但是在实际载荷作用下，齿轮副接触区域将变成椭圆区域，只考虑接触线上节点弹性变形也会造成计算的单齿啮合刚度存在较大误差。为了准确计算单齿啮合刚度，考虑整个工作齿面各个节点的弹性变形，将单个节点啮合刚度叠加得到单齿啮合刚度，则单齿啮合刚度表达式为

式中，kp、kg分别为主、从动齿轮的单齿啮合刚度；m为整个工作齿面的节点数量。

齿轮副1 个啮合周期内某一啮合转角下，存在t对齿接触，第j(1≤j≤t)对轮齿接触对综合啮合刚度表达式为

齿轮副的多齿综合啮合刚度表达式为

1.2 传动误差

本文中基于显式动态分析方法提取弧齿锥齿轮在实际载荷作用下小、大轮实际啮合转角，传动误差定义为从动轮的实际转角与按照理论传动比计算的理论转角差值。弧齿锥齿轮副传动误差表达式为[14]

式中，φ1、φ2 分别为小、大轮啮合转角； 分别为小、大轮初始啮合转角；z1、z2分别为小、大轮的齿数。

2 显式动态分析计算时变啮合刚度和传动误差

2.1 模型参数及网格划分

以一对弧齿锥齿轮为例，相关几何参数如表1所示，齿轮副几何三维模型如图2所示。

图2 弧齿锥齿轮三维几何模型

表1 齿轮几何参数

为了减少运算量，将全齿模型简化为5 齿模型，然后使用有限元网格划分软件Hypermesh 对5 齿模型进行网格划分，图3所示即为所建立的5齿对有限元模型。

图3 5齿对有限元模型

2.2 显式动态分析参数设置

前处理参数设置是影响分析结果正确性的关键因素之一。为了准确仿真弧齿锥齿轮小、大轮动态啮合过程，并确保分析结果的正确性和运算的收敛性，显式动态分析参数设置如下：

（1）材料属性定义：材料属性可根据齿轮材料查取，动态分析需要设置3 个材料属性。本文中定义弹性模量2.06×105 MPa，泊松比0.3，材料密度7 850 kg/m3。

（2）分析类型定义：采用显式动力学分析算法，分析总时间为0.001 765 s。

（3）输出变量定义：定义输出小、大轮平动位移、角位移、角速度、接触力、接触应力等输出变量。设置对总时间均分为给定的间隔数，间隔数为250，即将整个分析均分为250 个时间步，每个时间步输出有限元结果1次，可以更精确和完整地研究弧齿锥齿轮啮合传动过程中的接触啮合特性。

（4）接触对定义：使用的弧齿锥齿轮小轮为主动轮，小轮凹面与大轮凸面啮合接触，采用面-面接触类型，分别5对接触，摩擦因数设为0.1。

（5）边界条件和载荷定义：依次在小、大轮旋转轴中心建立参考点，然后建立参考点与齿轮内表面之间的耦合约束，再分别将小、大轮参考点沿轴线转动自由度释放，约束其余5 个自由度。定义小轮绕自身轴线的转动角速度628.318 rad/s，定义大轮绕自身轴线转动自由度上负载方向负载转矩为1 913.277 N∙m，由于小轮转速和大轮负载都很大，故对小轮转速和大轮负载转矩都采用线性加载，降低载荷施加造成的初始冲击。

（6）单元类型定义：弧齿锥齿轮网格划分时不可避免会出现网格扭曲，而线性缩减积分单元（C3D8R）具有网格扭曲对分析精度影响不大的优点，故采用C3D8R单元。

2.3 齿面网格节点接触力与弹性变形量

采用显式动态分析过程中，各个时间步提取的节点位移包含了轮齿啮合传动中的刚体位移，因此，不能作为节点的弹性变形量。解决方法是分别建立一个与小、大轮随动的坐标系，随动坐标系相对于齿轮是静止的，即将分析结果变换到随动坐标系下可以消除节点位移包含的刚体位移，在随动坐标系下提取的节点位移就是节点弹性变形量。图4所示为某一时间步小轮凹面节点接触力和弹性变形云图；图5所示为同一时间步对应的大轮凸面节点接触力和弹性变形云图。

图4 小轮凹面节点接触力和综合弹性变形

图5 大轮凸面节点接触力和综合弹性变形

2.4 时变啮合刚度与传动误差

图6 1个啮合周期内时变啮合刚度

图6 所示为在所选定高速重载工况（输入转速为6 000 r/min，负载为1 913.277 N·m）下1 个啮合周期内综合啮合刚度随小轮转角变化曲线，下一周期的齿轮副时变啮合刚度继续重复当前周期的时变啮合刚度曲线。从图6中可以看出，弧齿锥齿轮副从2齿啮合区进入3 齿啮合区时，综合啮合刚度会逐渐增加，而在3齿啮合区内时变啮合刚度存在小的波动，其原因是在3齿啮合区内，小轮转角的不同取值、工作齿面啮合位置和轮齿间载荷分配会发生变化，导致工作齿面接触椭圆位置及区域大小也会变化，相应的齿轮副综合啮合刚度在同一齿数啮合区内也会出现波动。分别提取各个时间步小、大轮实际转角值，按照式（7）计算传动误差，图7 所示即为1 个啮合周期内的传动误差，单位为″（角秒），传动误差值为负值，表示大轮转角存在延迟。

图7 1个啮合周期内的传动误差

3 不同负载转矩下时变啮合刚度和传动误差的变化规律研究

保持小轮输入转速恒定不变，分别给大轮施加负载转矩956.64 N·m、1 913.28 N·m、2 868.92 N·m，计算得到不同负载转矩工况下的时变啮合刚度和传动误差，如图8、图10所示。图9、图11所示分别为不同负载转矩下时变啮合刚度和传动误差的峰-峰值。随着负载转矩变大，综合啮合刚度和传动误差值逐渐增大，而时变啮合刚度峰-峰值和传动误差峰-峰值（PPTE）波动逐渐变小。峰-峰值降低的主要原因是随着负载转矩增大，轮齿所受载荷增大，变形量增大，变形量的补偿作用使时变啮合刚度和传动误差波动减小。

图8 不同负载转矩下时变啮合刚度

图9 负载转矩对时变啮合刚度峰-峰值的影响

图10 不同负载转矩下传动误差

图11 负载转矩对传动误差峰-峰值的影响

4 预置接触参数对啮合刚度和传动误差的影响

由于弧齿锥齿轮齿面的空间复杂曲面特性，需要结合齿轮几何参数和加工参数才能生成大、小轮精确的齿面几何三维模型，而在计算弧齿锥齿轮加工参数时，为了确保齿轮副传动过程中具有良好的接触性能，需要预先选定一些如接触椭圆长轴长度、接触轨迹方向及传动比1阶导数等接触参数。不同接触参数会生成不同弧齿锥齿轮的齿面几何模型，齿轮副啮合刚度和传动误差也会相应地改变。目前，研究弧齿锥齿轮接触参数对啮合刚度和传动误差影响的相关文献还很少。本文中分别改变接触椭圆长轴长度和接触轨迹方向两个接触参数，在同一载荷工况下计算齿轮副的时变啮合刚度和传动误差。其中，接触椭圆长轴长度为瞬时接触区域椭圆的长轴长度，而接触轨迹方向为大轮齿面接触轨迹线与大轮根锥的夹角。图12 所示为大轮工作齿面接触椭圆长轴长度和接触轨迹方向示意图。

图12 大轮工作齿面接触椭圆长轴长度和接触轨迹方向

4.1 接触椭圆长轴长度对时变啮合刚度和传动误差的影响

图13、图14 所示分别为不同接触椭圆长轴长度系数对应的弧齿锥齿轮副时变啮合刚度和传动误差，其中，接触椭圆长轴长度系数为齿宽与接触椭圆长轴长度之比，即接触椭圆长轴长度系数越大，接触椭圆长轴长度越长。从图13、图14 中可以看出，随着接触椭圆长轴长度增加，综合啮合刚度增大，而传动误差幅值减小。综合啮合刚度增加的原因是接触椭圆长轴长度变长，瞬时接触椭圆区域变大，椭圆接触区域抵抗弹性变形的能力增大，啮合刚度增大，轮齿自身变形引起的角位移滞后量减小，传动误差减小。随着接触椭圆长轴长度变大，时变啮合刚度峰-峰值先增大再接近于稳定，而传动误差峰-峰值逐渐趋于稳定，如图15所示。

图13 不同接触椭圆长轴长度下时变啮合刚度

4.2 接触轨迹方向对时变啮合刚度和传动误差的影响

改变弧齿锥齿轮接触参数接触轨迹方向角，得到不同接触轨迹方向角下的时变啮合刚度、传动误差及实际承载重合度。其中，实际承载重合度计算表达式为

式中，T 为1 个啮合周期时间；n1、n2分别为1 个啮合周期内交替啮合的齿数；t1、t2分别为基于显式动态分析结果获得的n1、n2对齿啮合时间。

图14 不同接触椭圆长轴长度下传动误差

图15 接触椭圆长轴长度对时变啮合刚度和传动误差峰-峰值的影响

图16 不同接触轨迹方向下时变啮合刚度

如图16 所示，当接触轨迹方向角为15°、25°和35°时，综合啮合刚度值没有明显变化，而当方向角为45°、55°时，综合啮合刚度值存在突增现象。存在这种现象的原因是在接触轨迹方向角小于35°时，齿轮副重合度随接触轨迹方向角增加而增加的斜率很小，齿轮副重合度基本没有变化，在方向角大于35°时，重合度增大的斜率明显增大，如图18所示。所以，接触轨迹方向角为45°的重合度比方向角为35°的重合度增大了很多，综合啮合刚度明显增大。由图17 可知，接触轨迹方向对传动误差的幅值影响很小，随着方向角增大，传动误差变化很小。随着接触轨迹方向角增大，时变啮合刚度峰-峰值先增大后减小，而传动误差峰-峰值先减小后增大，如图19所示。接触轨迹方向对传动误差幅值影响较小，但对传动误差峰-峰值影响较大，传动误差峰-峰值在接触轨迹方向角25°很大的原因是出现了边缘接触，传动误差波动加剧。

图17 不同接触轨迹方向下传动误差

图18 接触轨迹方向对重合度的影响

图19 接触轨迹方向对时变啮合刚度峰-峰值和传动误差峰-峰值的影响

5 结论

通过改进弧齿锥齿轮时变啮合刚度计算方法，由单个节点啮合刚度叠加计算单齿啮合刚度；基于显式动态分析方法计算弧齿锥齿轮时变啮合刚度和传动误差，解决了计算多次、计算量大等问题。文中分析了负载转矩、接触椭圆长轴长度及接触轨迹方向对弧齿锥齿轮时变啮合刚度和传动误差的影响，结论如下：

（1）弧齿锥齿轮综合啮合刚度和传动误差随负载转矩增大呈增大趋势，而时变啮合刚度和传动误差的峰-峰值随负载转矩增大呈减小趋势，轮齿变形量增大的补偿作用减小了时变啮合刚度和传动误差波动。

（2）接触椭圆长轴长度对时变啮合刚度峰-峰值和传动误差峰-峰值的影响呈相反趋势，随接触椭圆长轴长度增大，时变啮合刚度峰-峰值逐渐增大再接近于稳定，而传动误差峰-峰值逐渐减小到趋于稳定。此外，随接触椭圆长轴长度增大，综合啮合刚度增大，而传动误差幅值减小，设计时应根据实际情况选取接触椭圆长轴长度。

（3）在不同接触轨迹方向范围内，时变啮合刚度存在突增现象，原因是随接触轨迹方向角增大过程中，齿轮副重合度增大斜率不一致，重合度的快速增大使综合啮合刚度突增。不同接触轨迹方向下传动误差幅值变化很小，但接触轨迹方向对时变啮合刚度峰-峰值和传动误差峰-峰值的影响也呈相反趋势，随着接触轨迹方向角增大，时变啮合刚度峰-峰值先增大后减小，传动误差峰-峰值反之。