直齿圆柱齿轮啮合刚度计算方法研究

陈思宇，谭儒龙，郭晓东，阚 磊

（1．重庆理工大学 机械工程学院，重庆 400054；2．陆装驻重庆地区第六军代室，重庆 400042）

刚度激励是齿轮动态激励的主要来源之一，是由齿轮传动过程中随时间周期变化的啮合刚度即时变啮合刚度所产生的［1－2］，确定齿轮的啮合刚度一直都是齿轮系统动力学研究的热点，对此国内外已有大量学者进行研究，并取得了不少成果。在已发表的文献中，齿轮啮合刚度的计算多采用材料力学法和有限元法，其中材料力学法是应用最早且是最快速的求解方法［3］。

在基于材料力学的啮合刚度计算研究中，日本学者通过对轮齿齿形的大幅简化，提出了能够实现啮合刚度快速计算的石川公式，但由于该计算模型未考虑轮齿基体变形对啮合刚度的影响，计算精度不够理想［4］。李亚鹏等［5］在石川公式的基础上对其进行改进，将轮齿基体的变形量考虑在内提出了新的啮合刚度计算方法。Weber［6］根据利用能量积分导出包含轮齿弯曲、剪切与压缩的综合变形计算方法。Cornell［7］在Weber的研究基础上，进一步提出了包含齿根圆角以及轮齿基础弹性变形的数值积分法。YANG等［8］基于势能原理，提出了齿轮啮合刚度计算模型，WU等［9］和TIAN［10］在该模型的基础上进一步研究了包含剪切变形的影响。陈再刚等［11］基于刚度与轮齿误差相互作用的思想，提出了更加符合实际的时变齿侧间隙啮合刚度计算模型。贵新成等［12］采用修正的齿面赫兹接触刚度计算方法，基于势能法建立了单轮齿对啮合综合刚度模型。CHEN Zaigang等［13］通过考虑齿形变形、齿接触变形、圆角－基础变形、齿轮体结构耦合效应以及齿形偏差等各种变形，提出了一种综合的、通用的解析齿轮啮合模型。

与材料力学理论模型相比，有限元模型可以模拟复杂的几何形状、边界条件及工况，能够获得更加接近实际的变形及详细的应力分布情况［11］。WANG J等［14］采用有限元法分析了不同载荷对渐开线直齿圆柱齿轮啮合刚度的影响。唐进元等［15］和雷敦财等［16］利用有限元法中的准静态计算方法，针对齿根裂纹、齿廓修形等对啮合刚度的影响进行了计算分析。陆凤霞等［17］针对时变啮合刚度的计算，在有限元法基础上提出了一种新的基于齿轮啮合特性的有限元网格划分方法。冯正玖等［18］针对齿轮副的啮合刚度，运用三维有限元法获得了直齿轮齿面的柔度矩阵，建立了直齿轮啮合副的线接触和面接触2种轮齿承载接触分析模型。

在齿轮啮合刚度计算中，材料力学方法计算耗时短，但计算准确度较低，以有限元方法为代表的数值计算方法计算准确度高，但计算耗时较长。因此，本文力图探寻一种在精度损失较小情况下，直齿圆柱齿轮啮合刚度的快速计算方法。本文从能量等效思想出发，提出单齿啮合刚度快速计算模型，并结合位移相容原理，形成直齿圆柱齿轮时变啮合刚度计算方法。最终对比本文中计算方法与有限元分析方法的结果表明，本文中提出的方法能够在不显著降低计算准确度的情况下，实现直齿圆柱齿轮时变啮合刚度的快速计算。

1 渐开线直齿圆柱齿轮单齿啮合刚度计算模型

由于齿轮在啮合过程中情况比较复杂，作为对实际啮合状况的一个合理简化，在分析齿轮变形时进行以下设定：载荷沿接触线的均匀分布；无安装制造误差及动载荷；渐开线直齿圆柱齿轮齿廓为二维平面曲线。

将渐开线直齿圆柱齿轮轮齿视为如图1所示的变截面悬臂梁模型。如图1所示，F是沿啮合线方向上的啮合作用力，α1是啮合力与齿厚方向的夹角，h是啮合作用点的齿厚的一半，hx是变截面处齿厚的一半，x是截面到啮合作用点的距离，dx则是截面微段长度，Rint是齿轮内孔半径。

图1 渐开线直齿圆柱齿轮轮齿截面悬臂梁模型示意图

将轮齿的变形等效成啮合力F引起的沿啮合线方向的弹簧变形，弹簧做工的弹性势能表达式为：

刚度K乘上弹簧变形量μ等于这个弹簧的弹簧反力：

结合式（1）～（2），在啮合力F作用下，轮齿发生的弯曲、剪切与沿齿高方向的轴向压缩变形而存储的等效弹性势能为：

式中：Kb、Ks、Ka分别是沿啮合线方向与轮齿弯曲、剪切、轴向压缩变形相对应的等效弹簧刚度；Eb、Es、Ea分别是啮合力沿啮合线方向做的功分别等效为轮齿弯曲变形、剪切变形以及轴向压缩变形引起的等效弹簧弹性势能。

结合式（2）～（5）联立求得弯曲、剪切、轴向压缩变形的等效弹性势能：

根据悬臂梁理论，将啮合力F处的应变能分为弯曲、剪切、轴向压缩应变能：

式中：K是量纲为1的量，与横截面形状和尺寸有关，渐开线直齿轮齿截面为矩形截面，故K取值1 2。E和G是材料的弹性模量和剪切模量，其中：

Fb、Fa、M分别是啮合力F沿轮齿齿厚方向的分力、垂直于齿厚方向的分力以及相对于宽度为dx的微截面力矩，结合图1表示为：

Ix、Ax是距离啮合点x处的截面惯性矩与截面面积，结合图1表示为：

式中的b是轮齿的宽度。其中的h与hx的值可通过计算直齿轮任意弦齿厚的方法得到，表达式为：

结合式（9）～（17）与图示1可得弯曲刚度Kb的计算表达式为：

剪切刚度Ks的计算表达式为：

轴向压缩刚度Ka的计算表达式为：

赫兹接触刚度Kh计算表达式为：

齿轮啮合刚度不仅受到轮齿变形的影响，也会受到轮体变形的影响。轮体变形引起的啮合线上等效刚度Kf为：

式中，δf为齿轮的轮体变形量［11］。

由式（22）～（26）可得渐开线齿轮单齿啮合刚度Ke的表达式：

式中：脚标1、2分别对应主动轮、从动轮。

2 渐开线直齿圆柱齿轮时变啮合刚度计算模型

为保证齿轮正确传动，任何瞬间至少要有一对轮齿啮合。如图2所示，在齿轮正常传动过程中，一对轮齿还未完全退出啮合处于啮合点B1时，下一对轮齿已经于啮合点B2开始进入啮合。

图2 渐开线直齿圆柱齿轮连续啮合传动示意图

利用重合度可得出单齿啮合区与多齿啮合区，重合度计算式为：

式中：αa1和αa2为齿顶圆的压力角；α′为分度圆压力角。

接触点B1和B2到圆心的距离为：

由于齿轮传动任意啮合点都处于如图2所示啮合线N1N2上，通过重合度可计算出当第二对轮齿开始进入啮合时，第一对轮齿所在啮合点B1位置。

基于位移相容原理，在从动轮变形转角为θ的情况下，参与啮合的第一对齿的变形量δ1与同时参与啮合的第二对齿的变形量δ2的表达式：

式中α1和 α2是啮合点B1与B2的压力角，由式（30）～（31）表明变形量是啮合副状态的函数。通过变形量与单齿啮合刚度值可得啮合齿对的啮合力表达式：

式中的脚标i表示参与啮合的第i对齿。

当全部啮合力与外力处于平衡时有：

根据啮合刚度定义，每单位齿宽的齿面法向载荷和每个轮齿齿面法向变形量的总和的比值，称为一对轮齿的啮合刚度Ke，表达式为：

由式（32）～（34）可得多齿啮合刚度表达式为：

结合式（30）～（35）简化可得：

式中的单齿啮合刚度Ki通过第1节公式获得。由该式可见，Ki、Li与αi是齿轮系统运动状态变量的函数，综合啮合刚度K具有非线性。

3 算例分析

通过有限元法的具体实例来验证文中所提出算法的准确性，采用如表1所示的3组齿对参数。

表1 3组齿对参数

3．1 齿轮有限元前处理模型

将齿对的三维实体几何模型导入到有限元前处理软件中，得到其有限元网格模型。由接触力学可知，接触应力在接触区附近高度集中，应力梯度较大，因此接触部分网格进行密化。根据弹性力学圣维南原理，在远离接触区部分，其受力分布情况对接触区应力分布影响较小，误差小于1%。因此，考虑到实例中的齿轮副参数重合度大于1小于2，故任何瞬间最多2对轮齿最少1对轮齿参与啮合，为合理分配计算量，则采用3对轮齿进行分析。对于单元类型，由于六面体在接触问题中的计算代价远小于四面体网格，同时在单元积分形式的选择上，由于网格的扭曲变形对分析精度会有影响，故可选择产生影响相对较小的减缩积分单元，因此，在本文中，综合考虑计算效率和计算精度后，采用六面体一次缩减积分单元（C3D8R）建立齿轮有限元模型。三齿渐开线直齿圆柱齿轮有限元模型如图3所示，并采用文献［19］的有限元方法对其进行静力学分析。

图3 模数1．5齿对有限元模型示意图

3．2 计算结果处理与分析

齿轮单齿啮合刚度的一般表达式为

式中：Fn是作用于齿廓曲面的法向接触力；un是单对轮齿的综合弹性变形。通过有限元后处理得到法向接触力Fn以及综合弹性变形量un之后，由式（37）计算出单齿一个啮合周期在不同啮合位置的啮合刚度，再通过多项式插值得到单齿啮合刚度曲线。将单齿啮合刚度曲线左右移动nΔα角度可得到多齿对的啮合刚度，并通过式（39）对其进行叠加得到综合啮合刚度曲线。

根据表1中第1组齿对参数，利用本文的计算方法完成主动轮和从动轮的轮齿刚度计算，结果如图4所示。从图4中可以看出，根据本文提出的计算方法所得结果与有限元法计算结果基本相符，误差低于6%。

图5为本文方法、有限元法获得的单齿啮合刚度曲线。由图5可知，在整个轮齿啮合过程中，啮入啮出区域单齿刚度较小，而中间部分对应的啮合刚度较强。图6为综合啮合刚度曲线。由图6可知，综合啮合刚度是在一定范围内呈周期性交替变化的。具体啮合刚度计算误差如表2所示。通过3组实例的对比分析，可以认为本文所提出的直齿圆柱齿轮啮合刚度计算方法计算结果较为准确。

图4 轮齿刚度曲线

图5 单齿啮合刚度曲线

图6 综合啮合刚度曲线

表2 啮合刚度计算误差

4 结论

1）基于能量等效原理，可以有效构建渐开线直齿圆柱齿轮单齿啮合刚度计算模型。

2）基于位移相容的条件，从齿轮单齿啮合刚度中可以导出综合啮合刚度计算模型。

3）提出的直齿圆柱齿轮啮合刚度计算方法计算结果较为准确。由于采用的是简单积分计算式，可以实现直齿圆柱齿轮啮合刚度快速计算。