基于时变啮合刚度的修形齿轮承载接触分析*

王 浩,韩振华,周酝晨,石万凯,单文桃,曹忠亮

(1.江苏理工学院机械工程学院,常州 213001;2.南京华秦光声科技有限责任公司,南京 210034;3.重庆大学机械传动国家重点实验室,重庆 400044)

0 引言

轮齿承载能力是评价齿轮传动品质的一项重要标准,齿间载荷分配系数和面接触应力的模型精度是齿轮承载接触力学计算结果精确性的前提与关键[1]。齿轮修形可以有效改善齿轮动载荷变化梯度,减少冲击、振动与噪音,对提高齿轮传动品质有十分重要的作用。然而,已有的齿轮承载接触分析力学解析模型大都只能考虑鼓形的齿向修形,较难计入齿廓修形的影响[2-4]。

对于齿轮承载过程中产生的齿面接触应力,学者常采用有限元法和解析法。其中,有限元法建模精度高,且计算结果精确,因而得到较为广泛的应用。朱才朝[5]利用基于势能法的三段函数式齿间载荷分配系数模型分析齿轮的弹流润滑特性;AZNAR等[6]采用MPC算法与局部网格细化法优化齿轮的承载接触有限元模型,提高了计算精度与效率;MAPER等[7]建立了包含齿廓修形的渐开线直齿圆柱齿轮承载接触有限元模型;唐进元等[8]利用有限元分析了正交面齿轮的重合度、载荷分布、接触应力等承载接触性能参数;白恩军等[9]建立了斜齿圆柱齿轮承载接触有限元模型。然而,有限元算法计算精度依赖于模型的网格质量和局部网格细化程度,存在计算效率低、收敛难的问题。

相较于有限元法,基于赫兹弹性接触理论的力学解析算法,计算快捷高效受到了越来越多的关注与研究。王会良等[10]基于齿轮TCA法和齿面柔度系数,分析齿轮承载接触特性,建立了齿轮LTCA理论;曹雪梅等[11]提出一种简化轮齿接触分析中非线性方程数量的分解算法,提高了计算效率;王羽达等[12]基于轮齿TCA模型和赫兹接触理论,计算了渐开线齿轮的动态接触应力。

综上所述,传统齿间载荷分配系数解析建模未考虑齿廓修形及其引起的单双齿接触区域变化,无法准确计算承载啮合齿对的齿间载荷分配系数与齿面接触应力,同时,LTCA方法大都仅在TCA分析阶段考虑齿廓修形,而在齿面柔度系数或承载接触分析时采用有限元方法。为此,本文以齿廓修形的渐开线直齿轮为研究对象,提出修形齿轮单双齿啮合状态的s判定方法与相应的啮合刚度计算方法,建立更为精确的修形齿轮齿间载荷分配系数与齿面接触应力解析算法,对不同修形高度、修形曲线与配对齿数条件下的齿轮算例进行承载接触分析,并与同参数下的有限元结果进行对比,验证解析算法的精度与稳定性。

1 齿廓修形的渐开线直齿轮模型

为建立齿廓修形渐开线直齿轮的齿廓曲线方程,需确定齿廓最大修形量Δmaxi、修形高度h、修形曲线,根据文献[13]最大修形量Δmaxi为:

Δmaxi=ci+0.04FtB-1

(1)

式中:Δmaxi、ci分别为最大修形量与修形修正量(i为1、2,用以区分参与啮合的主/从动轮),相应的主/从动轮修形修正量c1、c2分别取9与4[13],Ft为接触点的圆周力,B为齿轮宽度。

建立坐标系XdOYd,如图1所示,以Δx作为增量等距离散啮合线,并在圆周方向上计算出轮齿1上的齿廓对应点,若该点位于修形齿廓区域内,啮合线离散点的齿廓对应点处的修形量Δj为:

图1 基于主动轮的啮合线等距离散示意图

Δj=Δmaxi(x/L)k

(2)

式中:修形曲线变化系数k由选取的修形曲线确定[14],Linear、Yoshio、Walker与抛物线修形曲线变化系数k分别为1.0、1.2、1.5与2.0,x为啮合点M在啮合线方向上距离单双齿变化点A的距离,L为轮齿2的修形区域在高度方向上对应至啮合线A点至C点的线段长度(长修形),C点为轮齿1的啮合起始点;另外,根据文献[15],正弦修形曲线的啮合点修形量为:

Δj=Δmaxisin(πx/L)

(3)

确定齿廓对应点处的修形量Δj后,根据渐开线齿廓曲线展成原理,可建立修形齿廓曲线参数方程为:

(4)

式中:u为啮合线离散点对应渐开线齿廓处的参角,u=θ+αt,θ为啮合点处展角,αt为压力角。

2 修形齿轮单双齿啮合状态判定与啮合刚度

标准渐开线齿轮副在承载啮合传动过程中齿轮啮合刚度主要包括3部分:轮齿刚度(弯曲刚度kb、剪切刚度ks与径向压缩刚度ka)、齿轮接触刚度kc与齿轮基体刚度kf。单齿啮合刚度Ki为上述各刚度的串联形式,双齿啮合刚度Kall为两对轮齿啮合刚度的并联形式[16-17]。

(5)

(6)

式中:i为1、2,用于区分参与啮合的不同齿对,下标g、e分别表示主动轮与从动轮。

当采用齿廓修形时,原单齿啮合区间的齿轮啮合刚度根据式(5)进行计算,原双齿啮合区间由于修形齿廓参与啮合,啮合状态会发生变化,需要重新判定啮合状态:根据齿对啮合力与啮合刚度计算接触点处沿啮合线方向的变形量,然后与接触点处修形量进行对比,基于二者的大小关系判定修形齿轮副的单双齿啮合状态,并计算单双齿的啮合区间与啮合刚度,其中,修形齿轮啮合刚度定义为KTPM[18]。

如图2a所示,当主动轮顺时针转动时,齿对1首先进入啮合承担载荷,由于齿对2主动轮处存在修形齿廓,当齿对1沿啮合线方向的变形量δ1大于齿对2处主动轮的齿廓修形量Δ1时,齿对2将进入啮合并参与传动,此时两对轮齿共同承担负载扭矩Tc,并在啮合线方向上的两处作用点产生啮合力F1、F2,根据力矩平衡和广义胡克定律可得:

(7)

(a) 齿对2受力分析 (b) 齿对1受力分析

式中:K1、K2与δ1、δ2分别为齿对1、齿对2的单齿啮合刚度与接触点沿啮合线方向的变形量,rb2为从动轮的基圆半径。

与此同时,齿对1沿啮合线方向的变形量δ1即为整体啮合变形量δall,则δall=δ1>δ2,相应的修形齿轮啮合刚度KTPM为:

KTPM=Fn/δ1,δall=δ1>δ2

(8)

式中:Fn为齿轮副法向力,根据式(7)将式(8)推导为:

KTPM=(K1δ1+K2δ2)/δ1,δall=δ1>δ2

(9)

同时,齿对2在啮合点处啮合变形量δ2与修形量Δj之和等于齿对1处啮合变形量为δ1。

δ1=Δj+δ2

(10)

将式(10)带入式(9),可得到:

KTPM=K1+K2-K2Δj(δ1)-1,δall=δ1>δ2

(11)

将式(8)转化为关于变形量δ1的表达式,并代入式(11),可推导得出图2a所示齿对承载接触状态下的修形齿轮啮合刚度KTPM,即式(12)中条件为δall=δ1>δ2的表达式;同样地,如图2b所示的齿对承载接触状态,此时齿对2沿啮合线方向的变形量δ2即为整体啮合变形量δall,即δall=δ2>δ1,同理可推导出相应的修形齿轮啮合刚度KTPM,即式(12)中条件为δall=δ2>δ1的表达式。因此,图2所示齿对承载接触状态下的修形齿轮啮合刚度KTPM可写为:

(12)

式(12)即为齿廓修形渐开线直齿轮的啮合刚度分析方法,下文将建立单双齿啮合状态的判定模型,确定实际双齿啮合区间,然后依据式(12)计算双齿区间内的齿轮啮合刚度。

如图3所示,基于主动轮齿廓,对齿轮副啮合区间进行划分,分析齿对1、2的接触啮合状态及相应的啮合刚度,点a1为啮入点,点a3、a4为单双齿啮合变化点,点a2、a5为原双齿接触区域平分点,点a1～点a5将接触区域划分为A1、A2、A3、A4、A5。

(a) 接触区域A1参与啮合 (b) 接触区域A2参与啮合

如图3a所示,当齿对1以a1点即将进入接触区域A1参与啮合时,齿对1中从动轮待接触区域存在修形齿廓,若齿对2沿啮合线方向的变形量δ2大于齿对1中从动轮修形量Δ2,则齿对1接触参与啮合,即为双齿啮合状态,此时的修形齿轮啮合刚度KTPM可依据式(12)计算;若齿对2沿啮合线方向的变形量δ2小于齿对1中主动轮修形量Δ2,则齿对1不参与啮合,此时为单齿啮合状态(齿对2单独承载),齿对2的啮合刚度K2依据式(5)进行计算。

如图3b所示,当齿对1以a2点即将进入接触区域A2参与啮合时,齿对2中主动轮待接触区域存在修形齿廓,若齿对1沿啮合线方向的变形量δ1大于齿对2中主动轮修形量Δ1,则齿对2接触参与啮合,即为双齿啮合状态,此时的修形齿轮啮合刚度KTPM可依据式(17)计算;若齿对1沿啮合线方向的变形量δ1小于齿对2中主动轮修形量Δ1,则齿对2不参与啮合,此时为单齿啮合状态(齿对1单独承载),齿对1的啮合刚度K1依据式(5)进行计算。

如图3c所示,当齿对1以a3点即将进入接触区域A3参与啮合时,此时接触区域A3为单齿啮合状态,齿对1的啮合刚度依据式(5)进行计算。

综合图3所述,即可得到齿廓修形齿轮在各啮合位置的单双齿啮合状态判定与相应啮合刚度计算方法。

3 齿廓修形渐开线直齿轮承载接触分析

3.1 齿间载荷分配系数

3.1.1 未修形齿轮的齿间载荷分配系数

根据式(6)确定的双齿啮合刚度Kall与齿轮整体变形量δall,法向力Fn可表示为:

Fn=Kallδall=(K1+K2)δall

(13)

根据式(7)中确定的载荷、刚度以及变形关系,带入式(13)等式左侧,可推导出:

K1(δall-δ1)=K2(δall-δ2)

(14)

由于啮合齿对在啮合线方向上变形量δ1、δ2中的较大值与双齿啮合整体变形量δall相等,此时无论δall=δ1≥δ2或δall=δ2≥δ1,根据式(14)皆可推导出各啮合齿对在啮合线方向上的变形量δ1、δ2与双齿整体啮合变形量δall相等,即:

δall=δ1=δ2

(15)

进一步可根据式(7)中啮合力、刚度与变形的数学关系推导出:

Fi/Ki=Fn/Kall=δall

(16)

式中:Fi为接触齿对的啮合力(i为1、2)。

对式(16)进行比例互换,可得到未修形下的齿轮齿间载荷分配系数S为:

S=Fi/Fn=Ki/Kall

(17)

3.1.2 齿廓修形齿轮的齿间载荷分配系数

当齿廓修形后,原双齿啮合区域的啮合状态发生变化,根据标准渐开线齿轮时变啮合刚度式(5)与齿廓修形齿轮啮合刚度式(12),重新确定齿轮的单双齿啮合状态与相应啮合刚度,继而根据式(17)计算修形齿轮不同齿廓区域参与啮合的齿间载荷分配系数SLDF。根据图3,以齿对1为例,分析齿对1承载啮合下的齿间载荷分配情况,具体为:

(18)

式中:Ku为齿对1的啮合刚度,当齿对1处于单齿啮合状态时,其啮合刚度依据式(5)计算,即Ku=K1;当齿对1、2同时参与啮合时,此时计算Ku需考虑齿对2啮合刚度K2的叠加影响,即:

Ku=KTPM-K2

(19)

当齿对1在齿廓区域A1处啮合时,由于齿对2接触变形量δ2的影响,齿对1存在非接触与接触两种状态,若δ2>Δ2,则齿对1接触,此时为双齿啮合状态,将式(18)代入式(19)可得到齿对1的齿间载荷分配系数SLDF。

(20)

若δ2<Δ2,则齿对1为非接触状态,此时齿对1不分担载荷,则齿间载荷分配系数为0,即:

SLDF=0

(21)

当齿对1在齿廓区域A2处啮合时,由于齿对1接触变形量δ1的影响,齿对2存在非接触与接触两种啮合状态,若δ1>Δ1,则齿对1接触,此时为双齿啮合状态,齿对1的啮合齿廓为非修形区域,其齿间载荷分配系数根据式(22)进行计算。

(22)

若δ1<Δ1,则齿对2为非接触,此时齿对1承担全部载荷,齿间载荷分配系数为1,即:

SLDF=1

(23)

当齿对1在齿廓A3处啮合时,此时为单齿啮合区域,齿间载荷分配系数为1,即:

SLDF=1

(24)

当齿廓进入A4、A5段参与啮合时,齿对1皆为非修形齿廓参与啮合,此时仅考虑齿对1的啮合刚度,即Ku=K1,相应的齿间载荷分配系数为:

(25)

综上所述,在一个啮合周期内,根据齿对1在齿廓区域A1、A2、A3、A4、A5是否接触参与啮合及其啮合状态判定,计算参与啮合齿对的啮合刚度,根据式(18)～式(25),可得到修形渐开线直齿轮齿间载荷分配系数SLDF计算模型。

(26)

3.2 齿面接触应力

根据上述建立的修形齿轮齿间载荷分配系数模型,可计算啮合力Fi,继而根据赫兹接触理论计算齿面接触应力。

首先,需确定不同齿廓赫兹接触类型下接触点处的综合曲率半径Re,针对未修形的齿廓区域,根据欧拉萨瓦利公式[19],通过计算沿啮合线方向上基圆至啮合点的距离求得其曲率半径ρ1、ρ2,对于齿廓修形区域,以主动轮为例,根据微分几何与修形齿廓方程式(4),可推导出修形齿廓曲率CTPM与曲率半径ρTPM。

(27)

ρTPM=1/CTPM

(28)

β=u+tanu,rb1为主动轮基圆半径。

针对齿廓修形齿轮承载啮合时存在的不同赫兹接触类型,可推导出相应的综合曲率半径Re。

(29)

通过式(26)确定的齿间载荷分配系数SLDF,可计算出参与啮合齿对的接触力Fi,根据式(29)所得综合曲率半径Re,基于赫兹弹性接触理论,可推导出修形齿轮承载接触时的齿面接触应力σH与接触半宽b。

(30)

式中:Ee为综合弹性模量,具体表达式为:

(31)

式中:E1、E2与ν1、ν2分别为主、从动轮的弹性模量与泊松比。

如图4所示为该解析算法的流程图。

图4 齿间载荷分配系数与齿面接触应力计算流程图

4 计算结果分析与讨论

基于上述建立的算法,求解不同修形高度、修形曲线与配对齿数3种条件下不同算例的齿间载荷分配系数与齿面接触应力,并与同参数下的有限元计算结果进行对比评价,验证本文解析算法计算结果的准确性和稳定性。基于表1中抛物线修形齿轮的主要几何参数,每种条件下选取5个算例求解以上不同齿廓修形齿轮算例的齿间载荷分配系数(图5～图7)和齿面接触应力(图12～图14)。使用有限元法对比评价解析算法结果的准确性和稳定性,以验证算法的正确性。

表1 齿轮参数

图5 不同修形高度算例的齿间载荷分配系数

图6 不同修形曲线算例的齿间载荷分配系数

图7 不同配对齿数算例的齿间载荷分配系数

分析如图5～图7所示的齿间载荷分配系数结果可知:

(1)本文解析算法与有限元两种计算方法的齿间载荷分配系数具有相同变化趋势:由最小值逐渐增大至最大值1,此过程为双齿啮合区间;随后进入单齿啮合区间,即SLDF=1的水平线;再进入双齿啮合区,下降至最小值。同时啮合线各点所对应的误差数值较小。

(2)根据图5～图7数据结果,可得到不同齿廓修形算例在一个啮合区间内的齿间载荷分配系数误差均值,如图8所示,每种条件下5个算例误差均值变化范围分别为7.5%～9.2%、7.5%～8.0%、5.6%～7.8%,表明误差均值较小且不同算例的误差均值相近。算例结果分析表明,对于不同算例情况下的齿轮副,本文解析算法的齿间载荷分配系数计算精度较高,且算法具有较好的误差稳定性。

图8 一个啮合区间内的齿间载荷分配系数误差均值

图10 单齿啮合区结束点位置误差

(3)根据图5～图7数据结果,分析齿廓修形后每个算例的单齿啮合区位置误差,以此可进一步评价解析算法的计算精度。其中,齿廓修形后的单齿啮合区起始点位置误差、结束点位置误差与啮合区间位置误差分别如图9～图11所示。不同修形高度、修形曲线与配对齿数算例的单齿啮合区起始误差变化范围分别为5.8%～12.4%、5.1%～14.2%、5.8%～10.1%,单齿啮合区结束点位置误差变化范围分别9.2%～13.1%、6.1%～13.3%、9.7%～12.6%,单齿啮合区间位置误差变化范围分别为3.2%～12.5%、3.2%～13.0%、2.2%～8.8%,表明解析算法所得出的单双齿变换点及其区间大小与有限元结果之间误差较小,可较为准确地计算出单双齿啮合变换点以及啮合区间位置。因此,本文解析算法中依据修形量与变形建立的单双齿啮合区间判定方法具有较高的计算精度,可较为准确的计算单双齿变换点以及啮合区间位置。

(4)另外,由于配对齿数的变化,齿轮副有效啮合区间也会随之变化,分析图7中得到的横坐标实际啮合线长度可得到齿廓修形齿轮的有效啮合区间,相对于有限元结果,5个配对齿数算例的有效啮合区间误差分别为1.3%、1.0%、3.0%、3.1%、1.1%,相对误差较小,表明本文解析算法对不同配对齿数的齿轮副啮合区间计算结果的准确度较高。

分析如图12～图14所示的4种条件下不同齿廓修形算例在一个啮合区间内的齿面接触应力结果可知:

图12 不同修形高度算例的齿面接触应力

图13 不同修形曲线算例的齿面接触应力

图14 不同配对齿数算例的齿面接触应力

(1)解析算法的齿面接触应力结果与有限元结果在变化趋势具有较好的一致性;根据图12～图14的齿面接触应力结果,可得到3种条件下齿廓修形算例在一个啮合区间内的齿面接触应力误差均值(图15a)变化范围分别为7.0%～9.0%、8.4%～10.4%、9.0%～10.7%,误差均值较小,同时不同算例间误差均值相近,表明解析算法具有较高的齿面接触应力计算精度与误差稳定性。

(a) 单啮合区间 (b) 去除突变区域

(2)与此同时,在修形齿轮刚开始进入啮合、即将退出啮合以及修形曲线与渐开线过渡区域,有限元结果产生了局部应力突变,以图12～图14中的A、B区域为例,主要由以下两个原因导致:①刚开始啮入与即将退出啮合时,接触齿对啮合力较小,需要极小网格尺寸(低于微米1～2个量级)才能有效捕捉到该处的接触应力信息,建模困难,且啮入时齿顶部位存在几何尖点;②修形齿廓过渡区域对几何曲率产生一定程度的啮合畸变。因此,当求解步长位于上述位置附近时会直接影响有限元计算结果,导致有限元局部应力突变、计算误差较大,这是不可避免的。除去由有限元结果突变区域所带来的计算误差后,解析算法的齿面接触应力在一个啮合区间内的误差均值(图15b)变化范围分别为4.8%～6.4%、4.5%～6.8%、6.1%～8.6%,再次表明解析算法具有较高的齿面接触应力计算精度与误差稳定性。

5 结论

(1)本文针对齿廓修形的渐开线直齿轮,通过对比承载齿轮接触点沿啮合线方向的变形量与修形量,得到了修形齿轮单双齿啮合状态的判定方法与相应啮合刚度的计算方法,通过推导修形齿轮啮合力、刚度与变形的力学关系,建立了修形齿轮各啮合齿对精确的齿间载荷分配系数计算模型。

(2)基于齿间载荷分配系数计算了啮合齿对实际承担载荷,应用欧拉萨瓦利公式和赫兹接触理论,建立了齿廓修形齿轮啮合过程中的齿面接触应力模型。

(3)不同修形高度、修形曲线与配对齿数条件下多种齿廓修形齿轮的算例结果表明:本文所提出解析算法的齿间载荷分配系数、齿面接触应力与相应单双齿啮合区间结果与有限元方法结果变化规律相一致,相对误差较小,具有较高的计算精度与稳定性,为齿廓修形渐开线直齿轮的承载接触分析提供了一种高效便捷的新方法。