基于齿根圆角圆心所在位置的时变啮合刚度修正模型

张珂铭, 邵毅敏, 许 晋, 何 融, 李 亮

(1.重庆大学 机械传动国家重点实验室，重庆 400044； 2.北方车辆研究所，北京 100072)

齿轮广泛地应用于旋转机械设备，其振动噪声越来越受到各方关注。众所周知，啮合齿对变化产生周期性的刚度激励是主要的内部激励，对研究时变啮合刚度对于齿轮振动控制具有十分重要的意义。

目前对于啮合刚度的研究，可分为有限元法(FEM)和解析法(AM)。有限元法是借助有限元分析软件来计算齿轮的啮合刚度[1-3]，计算精度较高，但计算量巨大，结构参数改变后需花费大量的时间建模与计算。相较于有限元法(FEM)的复杂性和耗时性，解析法(AM)为计算啮合刚度提供了一个更高效的途径，被广泛应用于啮合刚度计算。

目前，很多学者对解析法计算轮齿啮合刚度进行了相应研究，Yang等[4]考虑了弯曲、轴向压缩、赫兹接触刚度，提出用势能法来计算外直齿圆柱齿轮的啮合刚度。 Tian等[5-6]在势能法中加入了剪切刚度，Saxena等[7]考虑了齿轮轮体变形来完善啮合刚度的算法，计算了齿根有剥落或裂纹时的啮合刚度，Ma等[8]提出了修缘齿轮的啮合刚度改进算法。现有的研究一般都是将轮齿考虑为一个基圆开始的悬臂梁。但是，实际齿轮对的啮合，因齿数、变位系数、齿根圆角等参数的不同，其齿轮齿廓线存在不开始于基圆的情况，如图1所示。当基圆大于齿根圆时，目前的啮合刚度解析模型忽略了齿根圆与基圆之间的轮齿部分的势能，将导致计算啮合刚度相对偏大；当基圆小于齿根圆时，因多考虑了基圆与齿根圆之间的变形能会导致计算所得啮合刚度偏小。陈再刚等[9-11]近年来的针对啮合刚度的进行了大量的研究，其研究主要考虑了修形参数以及空间裂纹啮合刚度的影响，提出了‘切片式’裂纹啮合刚度的算法，其研究都考虑了齿根圆与基圆不重合的情况，但并未涉及到齿根倒角以及齿根倒角圆心位置对啮合刚度影响。

针对存在的问题，Liang等[12]提出了考虑了齿根圆半径小于基圆半径，齿根圆半径大于基圆半径两种情况，采用一条直线来表示基圆与齿根圆之间的齿廓，对现有的解析模型进行了修正，如图1所示。但是，由于大多数情况下，齿轮齿根部分都有切于齿廓线与齿根圆的齿根圆角，简单地将其倒有圆角的齿廓曲线用一条直线替代会对刚度计算造成一定误差。万志国等[13]考虑了齿根圆角半径，对啮合刚度进行了修正，但其考虑的为齿根圆圆心恰好位于过基圆与齿廓线交点做z轴垂线的延长线上，如图1文献[13]中假定倒角圆所示，且认为齿根圆角半径r大于基圆与齿根圆之间齿形在z轴距离x1。然而，实际情况中，齿根圆角半径的圆心位置由齿数，模数，压力角，倒角圆半径和变位系数等决定，且存在倒角圆圆心位于基圆内与基圆外两种情况。同时，齿根圆角的位置也会对轮齿刚度计算造成一定影响，且圆角半径r并不一定大于x1，当r

图1 各文献研究中拟定轮齿几何图形

针对以上问题，本文提出了一种考虑倒角圆心位置的啮合刚度修正算法，考虑齿根圆角在基圆内与基圆外两种情况，对以往刚度模型进行了改进，精确计算了齿根圆与基圆之间的势能，减小了目前解析法计算刚度的存在的误差，适合不同齿轮参数，具有更广泛的适用性。

1 基于齿根圆角圆心位置的啮合刚度计算方法

本文基于势能法将轮齿简化为齿根圆上的变截面悬臂梁，如图2所示。图中，F为垂直于齿面的啮合力，d为啮合力在轮齿上的有效作用长度，h代表啮合力作用点位置处的半齿厚，hx为有效作用长度上的齿廓线上任意位置处的半齿厚，齿轮形变等效为悬臂梁在F作用下的弹簧形变存储在齿轮中，其储存的势能有3部分，分别为弯曲势能Ub，剪切势能Us，径向压缩势能Ua。分别可表示为

(1)

式中：ks、kb、ka分别为轮齿在啮合力F的作用下沿剪切、弯曲、沿齿高方向轴压缩变形方向上的等效刚度。

图2 直齿轮变截面悬臂梁模型

传统算法大都假设轮齿齿廓线起始于基圆，即齿根圆即基圆重合，如图2所示，忽略了齿根圆与基圆之间的势能对啮合刚度的影响，且在大多数情况下齿根圆与基圆并不重合，导致计算啮合刚度的误差较大。针对这个问题，文献[12]虽考虑了齿根圆与基圆之间存在的势能，但仅用近似的齿廓线进行模拟，提高了计算啮合刚度的精度。然而，这种算法在某些参数条件下会产生误差，为了更为精确计算在各种参数条件下轮齿啮合刚度，本文提出了基于齿根圆角圆圆心位置的轮齿时变啮合刚度的算法，根据齿根圆角圆心在基圆内与基圆外两种情况，改进了传统算法。

1.1 倒角圆圆心位于基圆外啮合刚度算法

图3 倒角圆圆心在基圆外轮齿示意图1

倒角圆圆心位于基圆外部，倒角圆与齿廓线切点位于基圆以上，存在齿根圆半径小于基圆半径Rd

第一积分项为啮合点到齿廓线与倒角圆切点之间的轮齿部分，(以往方法为啮合点到基圆与齿廓线交点之间的轮齿部分)，如图4中阴影部分①所示，在z轴上积分区间为[0,d′]，根据轮齿几何关系，积分上限d′可表示为

d′=Rb[cosα1+(α1+α2)sinα1]-

(2)

式中，Rb为基圆半径。

图4 倒角圆圆心在基圆外轮齿示意图2

如图3,4所示，根据齿轮啮合角、齿廓曲线等几何关系

(3)

(4)

(5)

式中：r为齿根圆角半径，Rd为齿根圆半径，Rp为分度圆直径，x0为齿轮的变位系数。

第二积分项为齿廓线与倒角圆切点到齿根圆部分，如图4中阴影部分②所示，在z轴上的积分区间为[0,d1](需要注意的是，当Rd≥Rb，r=0,时，d1=0，即齿根圆半径为0时得到第二积分项的积分区间为0，此时在啮合力作用下，由剪切、弯曲、沿齿高方向径向压缩变形所储存的形变能的第二积分项Ub2,Us2,Ua2均等于0)。

考虑基圆与齿根圆之间轮齿的势能，当倒角圆圆心位于基圆外时，在啮合力作用下，由剪切、弯曲、沿齿高方向径向压缩变形所储存的形变能可表示为

(6)

(7)

(8)

式中：Fb为啮合力沿直齿轮齿厚方向上的分力，Fa为啮合力沿垂直于齿厚方向上的分力，M为Fa作用在微截面上的力矩，E为材料弹性模量，G为材料剪切模量，Ix、Ax分别为距离倒角圆与齿廓线切点x处轮齿截面的惯性矩与截面面积，Ix1为距离倒角圆与齿廓线切点x1处轮齿截面的惯性矩与截面面积。

Ax=2hxL,Ax1=2hx′L

(9)

(10)

x=Rb[cosα+(α+α2)sinα-cosα2]

x1=rsinβ-rsinβ0

(11)

hx=Rb[(α+α2)cosα-sinα]

hx1=he+r[1-cosβ]

(12)

dx=Rb(α+α2)cosαdα, dx1=rcosβdβ

(13)

式中：r为倒角圆半径，β为x1对应点与倒角圆心连线与水平方向夹角，β0为x1=0时，x1对应齿廓线的点与倒角圆心连线与水平方向夹角，βend为x1=d1时，x1对应点跟倒角圆心连线与水平方向夹角，根据齿廓曲线与各参数之间的几何关系可知

(14)

根据式(13)和(14)，第一段积分项的积分单元dx可用dα来替代，积分区间[0,d]可用啮合角区间[α1,α1]来表示。第二段积分项的积分单元dx1可用dβ来替代，积分区间[0,d2]可用x1对应齿廓线上的点跟倒角圆心连线与水平方向夹角表区间[β0,βend]来表示。

考虑式(6)～(14)，啮合点处的弯曲势能、剪切势能以及轴向压缩势能可改写为(设α位于z轴左侧为正，右侧为负)

(15)

(16)

(17)

将式(15)～(17)代入式(1)，即可求得可得当Rd+r-Rb>0时，在啮合力F的作用下沿剪切、弯曲、沿齿高方向轴压缩变形方向上的等效刚度kb、ks、ka

(18)

(19)

(20)

啮合时，应考虑赫兹接触刚度，其表达式为

(21)

式中，E为齿轮弹性模量，L为齿宽，v为泊松比，根据Muskhelishvili理论，齿根圆部分受到恒定的随时间变化的作用力，在该作用力下，齿轮轮体变形量为δf

(22)

若忽略齿轮基体的柔性变形，将会造成所求啮合刚度偏大。其式(22)的各参数定义，请参考文献[10]。

齿轮轮体变形引起的啮合线上等效刚度为

(23)

综合轮齿弯曲变形、剪切变形、轴向压缩变形、赫兹接触变形以及轮体变形，其对应的啮合线上等效刚度，可表示为各个变形所对应刚度的串联形式。单齿啮合刚度可表示为

(24)

采用式(24)，可计算倒角圆圆心位于基圆外时的轮齿单齿啮合刚度。

1.2 倒角圆圆心位于基圆内啮合刚度算法

齿根圆角圆圆心在基圆内，即Rd+rr，且倒角圆半径大于等于零，因此此种条件下只包含了齿根圆半径小于Rd

图5 倒角圆圆心在基圆外轮齿示意图1

为了更为精确计算啮合刚度，与现有算法不同，新算法采用三段积分法来计算时变啮合刚度，第一段积分项为啮合点到基圆与齿廓线交点所代表的区间，如图6中阴影部分①所示。

图6 倒角圆圆心在基圆外轮齿示意图2

在z轴上积分区间为[0,d]，根据齿廓曲线与各个参数之间的几何关系，积分上限d可表示式为

d=Rb[cosα1+(α1+α2)sinα1-cosα2]

(25)

第二积分项为基圆与齿廓线交点到齿根圆角圆圆心与z轴垂线与齿廓线交点，如图6中阴影部分②所示，在z轴上积分区间为[0,d1]

(26)

式中，Ra为齿顶圆半径，dmax为齿廓线与齿顶圆交点处的交点处的啮合点处所对应的d。he为第一积分区间积分终点所对应的啮合点到z轴的距离，当Rd+r

he=Rbsinα2

(27)

第三积分项为齿根圆角圆圆心与z轴垂线与齿廓线交点到倒角圆与齿根圆切点, 如图6中阴影部分②所示，在z轴上积分区间为[0,d2]。

考虑基圆与齿根圆之间轮齿的势能，当倒角圆圆心位于基圆内时，基于如上所述的三段积分法，在啮合力作用下，由剪切、弯曲、沿齿高方向径向压缩变形所储存的形变能可表示为

Ub=

(28)

(29)

(30)

式中，E为齿轮弹性模量，Ix、Ax分别为距离基圆x处轮齿截面的惯性矩与截面面积，Ix1、Ax1为第二积分项中距离基圆x1处轮齿截面的惯性矩与截面面积。Ix2、Ax2为第三积分项中距离倒角圆与齿廓线切点x2处轮齿截面的惯性矩与截面面积。

Ax=2hxL,Ax1=2hx1L,Ax2=2hx2L

(31)

(32)

x=Rb[cosα+(α+α2)sinα]-

x2=rsinβ

(33)

hx=Rb[(α+α2)cosα-sinα],hx1=he

hx2=he+r(1-cosβ)

(34)

dx=Rb(α+α2)cosαdα, dx2=rcosβdβ

(35)

式中，r为倒角圆半径，β为x1对应啮合点与倒角圆心连线与水平方向夹角，β0为x2=0时，x2对应啮合点与倒角圆心连线与水平方向夹角，βend为x2=d2时，x2对应点与倒角圆心连线与水平方向夹角，根据倒角圆与齿廓曲线之间的几何关系可知，

(36)

根据式(35)和(36)，第一段积分项的积分单元dx可用dα来替代，积分区间[0,d]可用啮合角区间[α1,α2]来表示。第三段积分项的积分单元dx2可用dβ来替代，积分区间[0,d2]可用倒角到倒角圆心连线与水平方向夹角表区间[0,βend]来表示。

因此啮合点处的弯曲势能、剪切势能以及轴向压缩势能可改写为(设α位于z轴左侧为正，右侧为负)

(37)

(38)

(39)

将式(37)～(39)代入式(1)，即可求得可得当Rd+r

(40)

(41)

(42)

基于赫兹接触刚度与轮体变形刚度的计算，运用式(21)～(23)即可求得当齿根圆角圆圆心位于基圆外时的轮齿时变啮合刚度，再通过式(24)可计算倒角圆圆心位于基圆外时的轮齿单齿啮合刚度。

当有两对齿轮同时参与啮合时，综合啮合刚度可表示为

(43)

式中，i=1表示第一对齿轮啮合，i=2表示第二对齿轮啮合。

2 基于齿根圆位置的啮合刚度算法的验证

如表1所示，为倒角圆圆心在基圆外(Rd+r>Rb)的两对相互啮合齿轮1,2的基本参数，表2所示为倒角圆圆心在基圆内(Rd+r

表1 倒角圆圆心在基圆外(Rd+r>Rb)齿轮基本参数

表2 倒角圆圆心在基圆内(Rd+r

Tab.2 The basic parameters of the gear for center ofchamfer circle is inside of the base circle

采用表1所示参数，建立有限元法算刚度模型，如图7所示。综合考虑求解精度及计算消耗量，细化接触区域的单元尺寸，主动轮1与从动轮2对共有299 213个单元，300 587个节点，单位类型为shell181,并建立齿轮1和2有限元模型的啮合接触对。约束被驱动齿轮2内圈上所有节点x,y,z方向的平动以及转动自由度，在驱动齿轮1内孔周向各个节点上施加切向力，以此施以顺时驱动力矩，提取内孔所有节点的切向位移并求取平均Uy，根据式(44)可计算得到此啮合位置的啮合刚度

(44)

式中，T为输入转矩，Rb(drive)为驱动齿轮的基圆半径，θ为驱动齿轮基于初始啮合位置的转角位移，接触点绕中心的转角，Uy为内孔所有节点沿切线方向的位移，Rhub(drive)为驱动轮轴孔半径。

依次改变模型啮合位置，根据式(44)可分别求得各个啮合位置处的啮合刚度，即可得到通过有限元法求得的表1所示参数下的轮齿啮合刚度。表2所示齿轮参数亦通过同样算法获得其轮齿啮合刚度。

图7 有限元法计算啮合刚度模型示意图

图8、图9分别为三种算法的轮齿啮合刚度结果与有限元法计算结果的对比。即：表1、表2所示齿轮参数，分别采用文献[7](未考虑齿根圆与基圆之间的轮齿部分，考虑了基体变形)和文献[12](将齿根圆与基圆之间的轮齿部分的齿廓用一条直线替代)及本文新算法的结果，与有限元法的对比。

图8显示，表1所示齿轮参数，倒角圆圆心在基圆外时(Rd+r>Rb)，本文算法的结果与有限元计算结果对比，单齿啮合区啮合刚度最大误差值为0.92×107N/m，双齿啮合区所算得的啮合刚度最大误差值小于0.75×107N/m，单双齿区的结果误差在1%以内。文献[12]将基圆与齿根圆的齿廓曲线近似考虑为一构造直线，与有限元计算结果对比，所得的刚度偏大，刚度最大差值达5×107N/m。文献[7]由于并未考虑基圆与齿根圆部分的因素，轮齿啮合刚度误差明显偏大，最大误差值可达11×107N/m。

图8 表1参数下各算法啮合刚度对比

图9显示，表2所示齿轮参数下，当倒角圆圆心在基圆内时(Rd+r

图9 表2参数下各算法啮合刚度对比

3 动力学响应仿真

为研究不同刚度算法对动力学响应的影响，本文建立了一个考虑齿面摩擦的6自由度集中质量模型，如图10所示，其动力学方程可表示为

(45)

动态啮合力Fp(t)与齿面摩擦力Ff(t)可分别表示为

(46)

Ff(t)=λfFp

(47)

式中，Rp，Rg分别为主被动轮的基圆半径，mp，mg为主被动轮质量，Jp，Jg为主被动轮转动惯量，kmi为齿对i的啮合刚度，cmi为齿对i的啮合阻尼，kgx，kgy，kpx，kpy为两齿轮轴承的支撑刚度，ei为齿对i的误差。λ为齿轮摩擦力方向系数，Ff做功与Tp同向时取‘+1’，反之取为‘-1’(即啮合点过节点之前取‘-1’，位于节点取‘0’，过节点后取‘+1’),H为啮合点到节点的距离。f为等效摩擦因数，由于摩擦因数在啮合线上随着运行条件的变化只有轻微的变化，可取其平均值以达到近似目的，其取值计算按照国家指导标准GB/Z 22559.1—2008进行

(48)

其式(48)的各参数定义，请参考文献照国家指导标准GB/Z 22559.2—2008。

图10 两对齿轮啮合的6自由度非线性动力学模型

在表1,2所示参数下，将各方法所求得的离散的刚度拟合为函数，再在动力学方程中根据主动轮所转过的角度从啮合刚度函数中取出对应啮合位置处的刚度，将带入动力学方程(支撑刚度为2×108N/m，输入力矩为2 000 Nm,输入转速1 000 r/min)，采用Runge-Kutta法求解微分方程，得到表1,2所示参数下，轮齿1扭转方向上振动的动力学响应分别如图11,12所示。

图11,12显示，基于本文啮合刚度算法的动力学响应与有限元法结果基本一致，误差都小于1%。采用文献[12]啮合刚度算法的动力学响应，在单齿啮合区与有限元法较为接近，但在双齿啮合区存在一定误差。采用文献[7]算法，齿轮1扭转方向的振动幅值，在单双齿啮合区较其他三种方法明显偏小。为了定量分析图11,12中三种解析算法与有限元算法所求得的扭转方向的动力学响应的误差，利用差分法进行量化，即把两种参数下三种解析算法所求得的齿轮1扭转方向的动力学响应与有限元法所求得的齿轮1扭转方向的动力学响应作差，再求取其RMS值，结果如图13所示。

图11 参数1下轮齿1扭转方向上振动的动力学响应

图12 参数2下轮齿1扭转方向上振动的动力学响应

图13 各种解析算法与有限元法的动力学响应差分值

图13显示，在表1参数下，基于本文啮合刚度算法所求得的齿轮1扭转方向的动力学响应的差分误差为文献[12]算法求得的差分误差的47.53%, 为文献[7]的17.61%。在表2参数下，基于本文啮合刚度算法所求得的齿轮1扭转方向的动力学响应的差分误差为文献[12]的29.11%, 仅为文献[7]的10.01%，以上结果说明了基于本文算法减小了以往算法在动力学分析中求解的误差。

为了进一步研究各种解析法对动力学响应的影响，基于表1，2参数的各种刚度算法以及输入力矩为2 000 Nm的输入条件下，轮齿1旋转方向振动位移去均值RMS值随输入转速上升变化如图14，15所示。

图14 参数1下轮齿1旋转方向振动位移去均值 RMS值随速度上升变化图

Fig.14 The relation between the mean removal RMS value of torsional vibration of gear 1 and the rising speed (parameter 1)

图15 参数2下轮齿1旋转方向振动位移去均值RMS值随速度上升变化图

Fig.15 The relation between the mean removal RMS value of torsional vibration of gear 1 and the rising speed (parameter 2)

图14,15显示，在本文所求得啮合刚度与基于有限元所求得刚度下，轮齿1旋转方向振动位移去均值RMS值随转速变化趋势基本一致，只在幅值上有微小差异。在同转速，参数1条件下，最大误差小于8%；参数2条件下，最大误差小于6%，但基于文献[7]与[12]刚度所求得的随转速变化的RMS值与有限元法相比在幅值和趋势上都存在着一定的误差，且在各转速下的响应存在不一致现象，未考虑齿根圆与基圆之间势能的文献[7]相对其他算法误差最大。

以上结果表明，忽略或粗略考虑齿根圆与基圆之间的势能，会对动力学响应造成一定影响，为更精确计算啮合齿对的振动特性，应根据不同参数条件下的啮合齿对，精确考虑齿根圆与基圆之间的势能，进行啮合刚度的计算。

4 结 论

(1) 本文基于齿根圆角在基圆内与基圆外两种情况(Rf>Rb和Rf

(2) 本文通过与传统算法及有限元仿真结果对比，表明：本文提出的修正算法能更为精确地计算轮齿的时变啮合刚度与有限元法所求得的啮合刚度最大误差小于2%，减小了以往解析法计算啮合刚度的误差，并适用于不同齿轮参数，具有更广泛的适用性。

(3) 动力学仿真结果显示，不同的刚度算法会对分析结果造成一定影响，基于本文算法的动力学响应的结果误差较此前算法至少降低了50%，为了更精确地分析齿对的啮合振动特性，应基于齿根圆角圆心所在位置，精确考虑齿根圆与基圆之间的势能。