刘少杭，聂少武，蒋闯，柴福博

(河南科技大学机电工程学院，河南洛阳 471003)

0 前言

准双曲面齿轮具有重合度大、传动平稳、承载能力高等优点，广泛应用于各种传动装置，在汽车驱动后桥中的应用尤为广泛。近些年来，为了对汽车驱动桥准双曲面齿轮的啮合性能进行评价，国内外学者对准双曲面齿轮的齿面接触分析(Tooth Contact Analysis，TCA)进行了一系列的研究。王星等人[1-2]建立HGT准双曲面齿轮切齿加工数学模型，对齿面接触和边缘接触进行了分析。杜进辅等[3]通过分析克林贝格制摆线齿准双曲面齿轮，建立一种改进的齿面接触分析模型并对其进行分析，避免了TCA结果出现几何上的不准确。杜进辅等[4]通过建立奥利康制准双曲面齿轮的全齿面模型，对其进行了齿面接触分析。SHIH等[5]建立了采用端面铣削和端面滚切加工的准双曲面齿轮齿面接触分析的统一数学模型。王峥嵘[6]对采用端面滚齿法加工准双曲面齿轮提出了两种基于数值迭代的齿面接触分析方法，避免了复杂的主曲率主方向的推导。卓耀彬等[7]建立了适合准静态齿面接触分析的准双曲面齿轮传动系统，并对其进行了准静态齿面接触分析。苏进展等[8]对弧齿锥齿轮进行了数值化齿面接触分析，实现了数值化齿面的几何特征修改。

目前，对准双曲面齿轮副啮合性能的评价通常是采用轮齿接触分析(TCA)方法。传统TCA方法是通过解析法求解齿面接触区，因涉及到主方向、主曲率等复杂曲率计算，推导过程十分繁琐且计算量大，甚至有时出现齿面接触椭圆求解异常的情况。因此，本文作者提出一种基于解析法和有限元法融合的齿面接触分析方法，即利用解析法计算齿面传动误差，通过有限元软件分析加载接触，得到齿面接触区，通过与专业软件进行对比，验证该方法的可行性。

1 HFT数学模型的建立

根据砂轮结构建立刀盘坐标系[9]，如图1所示。其中，坐标系St(Xt,Yt,Zt)与刀盘固连；Gg和Gp是内外刀刀刃上的一点；rg和rp分别是内外刀刀尖半径；αg和αp分别为内外刀齿齿形角；αg取正值；αp取负值；ug和up为刀齿曲面参数。

在HFT法加工过程中，根据刀具与工件的相对运动和位置变化，建立了小轮磨齿加工数学模型，如图2所示。

坐标系Sm(Xm,Ym,Zm)与加工机床固连称为机床坐标系，Om与机床中心重合，Sc(Xc,Yc,Zc)与摇台固连称为摇台坐标系，加工前坐标Sm和Sc重合。j为刀倾角，i为刀转角，q1为角向刀位，Sr1为径向刀位，Em1为垂直轮位，Xb为床位，Xp1为水平轮位，γm为工件安装角度，θ为刀具齿面参数，Sd与工件固连称为工件坐标系，通过展成运动关系可知φ=Rbφ，Rb为切齿滚比。

根据图1刀盘坐标系中的矢量关系，可以推导出刀具曲面方程以及单位法矢。

(1)

nt=[-cosαkcosθk-cosαksinθksinαk]

(2)

式中：k=g、p,分别代表内刀和外刀。

由St坐标系中的刀具曲面方程和单位法矢经过一系列坐标变换和矢量运算得到小轮齿面方程和单位法矢。

r1(uk,θk,φ)=M1hMhsMsmMmcMcjMjtrt(uk,θk)

(3)

n1(uk,θk,φ)=LdhLhsLsmLmcLcjLjtnt(uk,θk)

(4)

式中:Ldh、Lhs、Lsm、Lmc、Lcj、Ljt分别为矩阵Mdh、Mhs、Msm、Mmc、Mcj、Mjt去掉最后一行及最后一列得到的3×3矩阵。

大轮采用成形法加工，成形法加工数学模型是刀倾法加工数学模型的一种特殊情况，因此只需令图2中的j=0、Xb=0、Em1=0、Rb=0、φ=0、φ=0，即可推导出大轮齿面方程和单位法矢。

2 解析法求解传动误差

齿面接触分析是通过数学方法模拟出齿轮副的啮合过程，依据准双曲面齿轮副在啮合过程中的位置和运动关系[10]，建立准双曲面齿轮啮合数学模型，如图3所示。

图3中小轮轴线和大轮轴线夹角为90°，轴线之间距离为ΔE，S2与大轮固连，S1与小轮固连，Sb和Sd为辅助坐标系，啮合过程中S2绕大轮轴线旋转φ2，S1绕小轮轴线旋转φ1。

通过上述建立的刀倾法加工数学模型即可推导出小轮齿面方程r1(u1,θ1)、小轮法矢n1(u1,θ1)和大轮齿面方程r2(u2,θ2)、大轮法矢n2(u2,θ2)。分别将它们进行坐标变换，变换到固定坐标系Sa中，得到如下：

ra1(u1,θ1,φ1)=MadMd1r1(u1,θ1)

(5)

ra2(u2,θ2,φ2)=MabMb2r2(u2,θ2)

(6)

na1(u1,θ1,φ1)=LadLd1n1(u1,θ1)

(7)

na2(u2,θ2,φ2)=LabLb2n2(u2,θ2)

(8)

式中：Lad、Ld1、Lab、Lb2分别为矩阵Mad、Md1、Mab、Mb2去掉最后一行及最后一列得到的3×3矩阵。

在齿轮副啮合过程中，齿面Σ1和齿面Σ2持续相切接触，故在任一时刻两齿面都有公共接触点，在公共接触点处两齿面有相同的径矢ram和相同法矢nam。

ra1(u1,θ1,φ1)=ra2(u2,θ2,φ2)

(9)

na1(u1,θ1,φ1)=na2(u2,θ2,φ2)

(10)

由式(9)(10)可得到5个相互独立的标量方程，5个标量方程中含有6个未知数，通过给定φ1值，求解出其他5个未知数(u1,u2,θ1,θ2,φ2)，得到第一个接触点，将它设为参考点。以此为初值，通过以一定的步长改变φ1值求解出齿面上的所有点。

由传动误差的定义，得到其表达式[11-12]如下:

(11)

3 齿面接触区有限元仿真

通过加工数学模型求解出齿面方程，采用旋转投影原理进行网格划分进而求解出齿面点，将齿面点导入UG中建立准双曲面齿轮模型[13]。为便于有限元软件的仿真分析，将准双曲面齿轮模型分割为单齿模型。为便于六面体网格的划分，建立辅助面将单齿模型进行切分，最后将模型和辅助面导出为(.x_t)文件。

ABAQUS齿面接触有限元仿真流程如图4所示。仿真流程主要分为前处理、求解和后处理3个部分。在前处理过程中，通过导入上述(.x_t)文件，进行装配，完成几何模型的构建。通过设置材料密度、弹性模量、泊松比完成材料参数的定义。通过对大小轮各边布点，划分网格完成网格系统的构建。

在求解过程中，求解器的设定采用静力通用类型。连接关系采用面与面接触。边界条件的设定通过设置边界约束以及施加负载扭矩实现。

后处理过程是当上述所有操作设置完毕后，提交作业，得到场输出数据中的接触应力。

接触椭圆是齿轮副啮合产生的瞬时接触区，不同时刻瞬时接触椭圆的集合称为齿面接触区[14-15]。因此，本文作者通过使用Python语言编写脚本，对分析结果进行提取。Python提取ABAQUS中输出结果的流程如图5所示。

Python提取流程图主要分为读取数据和输出数据两部分。在读取数据阶段，首先打开odb文件读取场输出结果中的齿面接触应力，然后定义2个数组，分别储存单元编号和对应的接触应力值，最后关闭结果文件。在输出数据阶段，首先打开新建的储存文件，然后通过比较每个单元在不同时刻的接触应力，找到最大接触应力，并将每个单元的最大接触应力值写入文件中，最后关闭文件。

ABAQUS软件在分析结束后，关闭odb文件，运行上述的脚本文件，即可得到主从动轮接触面上的每个单元在整个分析过程中不同时刻的最大接触应力。每一时刻显示的接触应力云图形状等同于齿面瞬时接触椭圆。通过提取不同时刻的最大接触应力，得到齿面接触区。

4 算例分析

以一对HFT法加工的准双曲面齿轮为例，分别进行建模以及有限元接触分析。表1所示为准双曲面齿轮几何参数，表2所示为加工参数。

表1 准双曲面齿轮几何参数

表2 加工参数

利用解析法对轮齿接触原理进行研究，通过 MATLAB软件编写传动误差程序。图6所示为由MATLAB程序计算得出的传动误差，图7所示为由格里森CAGE软件得到的传动误差。

由图6—图7可以看出：解析法得到的传动误差幅值(交点对应的纵坐标)为28 μrad，格里森CAGE软件得到的传动误差幅值为28.5 μrad，两者传动误差幅值基本相同，验证了解析法求解传动误差的正确性。

在UG中建立准双曲面齿轮模型，如图8所示。

将建立好的模型导入有限元软件ABAQUS中进行负载接触分析，得到分析结果后再进行后处理，即对分析结果进行提取，得到齿面接触区，如图9—图11所示。

然后将准双曲面齿轮的几何参数以及加工参数输入到软件MASTA中，建立齿轮副模型，通过施加不同载荷，得到齿轮副在不同负载下的齿面接触区，如图12所示。

通过对比ABAQUS仿真齿面接触区与MASTA软件仿真结果，发现接触区都是位于齿面中部，成内对角接触，齿面接触区形状大小趋于一致，验证了文中所提出的齿面接触分析方法的正确性。

为进一步对有限元仿真结果进行验证，对表1中准双曲面齿轮进行磨齿加工和齿面误差测量实验，如图13和图14所示。通过对磨齿后的齿面进行测量反调，保证了实际加工齿面与设计齿面一致，最后在滚检机T50上进行了齿面滚检实验(加载扭矩40 N·m)，得到了轻载下大轮凹面接触区，如图15所示。

通过对比滚检结果图15(b)与有限元仿真结果图11(b)可以看出：仿真接触区的形状、位置、大小与滚检结果趋于一致，进一步证明了所提出的齿面接触分析方法的正确性与可行性。

5 结论

基于解析法和有限元法对准双曲面齿轮齿面接触分析进行了研究，通过解析法求解出了传动误差曲线，通过有限元仿真获得了齿面三维接触区，避免了解析法在求解接触区椭圆时复杂的曲率计算，同时使得齿面接触区显示更加真实直观。通过对比仿真结果与MASTA软件结果，以及齿面加工滚检实验，验证了所提齿面接触分析方法的正确性与可行性。