考虑齿轮轴变形的斜齿轮接触分析

白恩军，谢里阳，佟安时，白鑫

(1.东北大学 机械工程与自动化学院 现代设计与分析研究所，辽宁 沈阳110819;2.东北大学 航空动力装备振动及控制教育部重点实验室，辽宁沈阳110819)

0 引言

齿轮系统作为变速器的动力和运动传递中枢，其动力学行为和工作性能对变速器有重要影响［1］。齿轮的工作性能主要受齿轮接触状态的影响，而齿轮轴的变形，直接影响齿轮的接触状态［2］。在齿轮轴工作过程中，齿轮轴要承受斜齿轮的轴向力、径向力和切向力，在重载、工况极其恶劣的条件下，会产生很大的弯矩和扭矩，改变原有齿轮的中心距，本为平行的齿轮轴会形成一定角度，直接影响齿轮传动系统的工作性能。分析齿轮轴的变形，对于传动性能有较高要求的传动系统显得尤为重要。

国内外许多学者采用有限元方法分析了齿面接触状态。Hwang 等［3］建立了多齿有限元模型，通过改变齿轮转动角度，计算单个轮齿啮合周期的最大接触应力，并与标准进行了对比分析。Gonzalez-Perez 等［4］建立了2 对轮齿有限元模型，分析了齿面接触应力和Tresca 应力分布。Patil 等［5］建立了3 对轮齿有限元模型，通过改变轮齿螺旋角和摩擦系数，分析了接触应力的变化趋势。李杰等［6］建立了变速器有限元模型，分析了齿轮轴变形对齿面接触状态的影响。这些研究虽然分析了单齿对啮合、多齿对啮合的接触状态，或者含齿轮轴的静态齿轮接触应力分析，但没有分析齿轮轴变形对载荷沿齿向分布、重合度和齿轮时变啮合刚度等的影响。

本文针对典型变速器建立了含齿轮轴的齿轮接触模型，以发动机最大输出扭矩为输入载荷，1 挡输出齿轮为分析对象，分析了齿轮轴在综合变形情况下的齿向载荷分布、接触应力分布、齿根弯曲应力分布、重合度和时变啮合刚度。

1 变速器齿轮轴的变形

1.1 有限差分法计算轴弯曲变形量

在受弯轴(或杆件)为变截面或载荷比较复杂的情况下，用积分法求解弯曲变形计算比较困难，而有限差分法是一种简单有效的方法［7］。

如图1所示，齿轮轴的挠曲线为υ =f(x)，取横坐标分别为0，1，2，…，各相邻点间距离Δx 可以相等，也可以不等。在任一点i 截面处的弯矩为Mi，任意截面处挠度与弯矩可用如下有限差分方程［8］表示:

图1 轴的弯曲变形简图Fig.1 Deformation of gear shaft

式中:υi为xi截面处轴心挠度;Mi为xi截面处轴心弯矩;E 为弹性模量;Ii为xi处轴截面惯性矩。对于实心轴，Ii=πD4/64，D 为轴圆截面的直径。对于空心轴，Ii=πD4(1 -α)4/64，α =d/D，D 和d 分别为空心圆截面的外径和内径。

1.2 轴的扭转变形

轴的扭转变形是横截面绕轴线的相对转动。对于阶梯轴，需要分段求和计算相对扭转角［8］，即

式中:Ti为扭矩;li为等直径段轴长;φi为相对扭转角;G 为剪切弹性模量;Iρi为轴截面极惯性矩;对于实心轴，Iρi=πD4/32，D 为轴圆截面的直径;对于空心轴，Iρi=πD4(1 -α)4/32，α=d/D.

2 齿轮啮合传动有限元模型

以输入齿轮轴与1 挡输出齿轮组成的1 对斜齿轮副为对象，不考虑齿侧间隙、加工及安装误差的影响，研究其在齿轮轴变形影响下的接触特性。斜齿轮副几何参数如表1所示。

齿轮轮齿有限元建模是齿轮接触分析的关键，对重要区域细化出规整的单元，能够得到更为准确的分析结果。根据文献［9］的方法，并进行一定调整，能够更好地控制单元质量，下面详细阐述轮齿的建模方法。

步骤1 将轮齿进行分割，在轮齿的端面做辅助线1 ～9，将端面分割成10 部分，以控制单元节点走势，如图2(a).

步骤2 在两侧端面划分2D 单元，控制沿齿高方向设置的单元数量，在齿面及剖分表面沿齿向方向预置单元数量划分2D 单元，同时考虑齿高方向单元节点数量，划分结果如图2(b).

表1 斜齿轮传动参数Tab.1 Parameters of helical gears

步骤3 将端面2D 单元沿齿向方向单元节点路径进行拉伸，并与另一端面单元对应，拉伸成的3D 单元轮齿模型，如图2(c). 将轮齿3D 单元进行阵列，生成齿圈模型，如图2(d).

图2 齿轮建模方法Fig.2 Gear modelling method

根据标准渐开线方程和渐开线圆柱齿轮啮合原理［10］以及有限元建模方法［11］，考虑整体模型的单元数量，建立如图3所示斜齿轮副啮合接触有限元分析模型。

图3 齿轮啮合有限元模型Fig.3 Finite element model of gear mesh

在有限元接触分析中，由于二次单元会导致等效节点接触力在角节点和边中节点之间的震荡［12］，模型采用8 节点6 面体减缩积分单元。网格扭曲对结果影响小，并且求解速度快。齿轮材料密度ρ=7 800 kg/m3，弹性模量为207 GPa，泊松比为0.25. 模型中输入端花键定义参考点，建立耦合约束，具有转动自由度。齿轮轴与轴承接触处建立耦合约束，具有旋转自由度。从动齿轮与输出轴绑定约束。啮合齿对间设置面-面接触［13］。其他未工作齿轮与轴进行绑定约束。各参考点的耦合约束建立方法参考示意图如图4，在轴心建立坐标系Nx1y1z1，N 为参考点(根据实际情况释放参考点自由度)，将耦合面耦合到点N，T1(ω1)为驱动载荷(若无驱动则不考虑)。为分析齿轮轴变形对齿轮副接触特性的影响，本文建立了刚性齿轮轴和柔性齿轮轴两种分析模型。

图4 边界条件示意图Fig.4 Boundary conditions

3 仿真结果分析与讨论

3.1 齿轮轴变形有限元结果与有限差分法计算结果对比

为模拟齿轮轴的受载变形，将输入扭矩施加于输入端花键齿轮，限制输出轴差速器小齿轮的转动自由度，以实现齿轮副的接触分析。表2为输出齿轮轴弯曲变形和扭转变形的经典计算和有限元计算结果。

表2 齿轮轴变形计算结果Tab.2 Caculated results of gear shaft deformation

比较表2中数据可见，有限差分法结果大于有限元方法计算结果，但二者变化趋势相同。这是因为经典计算中将齿轮视为刚体，没有考虑齿轮的受载变形。有限元法中考虑了主动齿轮轮齿和从动齿轮的变形，这些因素将使有限元计算结果偏小，这也说明有限元法更能真实反应齿轮轴的变形。

3.2 载荷分布

使用传统的方法很难预测齿面载荷分布和接触应力分布，而有限元法通过正确定义齿轮几何形状、载荷和边界条件，能够精确预测载荷和应力分布［14］。

图5显示了刚性与柔性齿轮轴时载荷在齿宽上的分布。当齿轮轴为刚性时，齿1 所受载荷在齿宽中间位置相对边缘位置小，齿2 所受载荷从齿宽中间向两侧递减，比较均匀分布在整个齿宽上。当考虑齿轮轴变形时，齿1 所受载荷在中间位置有所增加，齿2 所受载荷从齿宽左端向右端有递减趋势，载荷发生偏移。这是因为当考虑齿轮轴变形时，由于轴受载发生弯曲和扭转变形，使得齿轮发生倾斜，改变了原有齿轮的接触状态，引起偏载现象。

图6为刚性与柔性齿轮轴在相同扭矩载荷作用下产生的等效应力。从图6中可以看出，当考虑齿轮轴变形时，齿对啮合线上Mises 应力最大值为1 909 MPa，相对刚性轴时增加448 MPa，并且齿对啮合线上Mises 应力分布相对于刚性轴时产生偏移，全齿啮合区域相对于刚性轴时接触区域在齿宽左端略有增加，从左向右呈现递减。这和沿齿宽方向的接触力分布分析趋势相同，并且从轮齿等效应力分布上亦可以看出接触面积从左到右呈现递减趋势。

图5 齿宽方向载荷分布Fig.5 The distribution of contact force along tooth width

图6 齿轮接触Mises 应力Fig.6 The von Mises stress of gear

图7为刚性与柔性齿轮轴双齿啮合状态时的齿根弯曲应力。图8为刚性与柔性齿轮轴接近单齿啮合状态时的齿根弯曲应力。图8中，由于加载，轮齿和轮体变形，使得齿轮啮合重合度大于2，原本为单齿啮合状态变成了双齿啮合，如图9所示。

从图7中可以看出，当齿轮轴为柔性时，齿根弯曲应力分布相对刚性轴时向左侧发生偏移，并且弯曲应力大小由刚性轴时的最大值717.5 MPa 增加到823.3 MPa. 从图8中接近单齿啮合状态也可以得到图7中相似的结论。由于齿轮轴的变形引起的偏载使得齿轮在靠近左侧区域载荷增加，易引起破坏现象。

图7 齿根弯曲应力Fig.7 Bending stress of tooth root

图8 齿根弯曲应力Fig.8 Bending stress of tooth root

当齿轮受载转动时，齿轮轮齿将经历从进入啮合产生接触力到退出啮合接触力消失的过程，这个过程所消耗的时间便是轮齿实际啮合的时间。通过测出单齿参加啮合的时间ΔT 和相邻两齿进入啮合的时间差Δt，则重合度［12］为

式中:ΔT 可以写成角位移Δθ1与角速度ω 的商;Δt可以写成角位移Δθ2与角速度ω 的商。则(3)式可以进一步转化为

式中:Δθ1为一对齿完成啮合过程的主动(从动)齿轮角位移;Δθ2为相邻一对轮齿进入啮合的主动(从动)齿轮角位移差;ω 为主动(从动)齿轮角速度。

图9为柔性齿轮轴时从动齿轮轮齿的接触力随转角的变化曲线。

图9 轮齿接触力曲线Fig.9 Contact force curves of each tooth pair during meshing

图9中从动齿轮齿3 经历了一个完整的从开始接触到脱离接触的过程，同时计算了齿3 接触时相邻的轮齿接触力。将齿3 经过完整啮合的转角Δθ1与相邻齿4 进入啮合的转角差Δθ2代入(4)式，可得考虑齿轮轴变形时的重合度为2.273. 比较文中齿轮设计重合度为1.983，可见由于齿轮轴和齿轮的变形，其实际啮合重合度比设计重合度增大了14.62%，其结果更接近于实际情况［15］。

3.3 啮合刚度

啮合刚度是指使一对或几对同时啮合的轮齿在1 mm 齿宽上产生1 μm 挠度所需的载荷，其计算公式［16］可表达为

式中:cγα为齿轮切向啮合刚度;c'为单齿啮合刚度均值;εT为端面重合度;CM为理论修正系数;CR为轮坯结构系数;CB为基本齿廓系数;β 为螺旋角;zn1和zn2分别为主、从动齿轮当量齿数;q'、C1、C2、C3、C4、C5、C6、C7、C8、C9各系数值见ISO 6336-1 标准［16］。

根据以上公式计算的齿轮啮合刚度是均值啮合刚度，而齿轮的啮合刚度是时变的，通过有限元方法可以求得齿轮在任意啮合位置的齿轮啮合刚度。Amezketa 等［17］通过将齿轮啮合点延圆周方向的转角位移转换到延啮合线方向的位移，进而很容易地求解齿轮啮合刚度，公式为

式中:kgear为齿轮切向啮合刚度;Ft为切向力;B 为齿宽;rF为圆周力作用点半径;rb为基圆半径;δrF为啮合点沿圆周方向位移。

计算齿轮啮合刚度时，大齿轮转角从0° ～22.7561°为1 个周期，即1 个齿从齿根开始接触，到齿顶脱离啮合结束，具体方法参见文献［3］.

从图10 中可见，当齿轮轴为柔性时，啮合刚度相对刚性时有所减小。这是因为柔性轴时，由于轴的变形，使齿轮中心距增大，啮合点相对远离轴心，轮齿啮合刚度下降。当齿轮轴为柔性时，啮合刚度曲线相对平缓，其刚度均值为15.235 4 N/(mm·μm).刚性轴时，啮合刚度曲线相对波动较大，其刚度均值为18.343 1 N/(mm·μm). 通过(5)式～(7)式计算的均值刚度为17.564 5 N/(mm·μm)，可以看出有限元方法计算的刚度值与ISO 6336-1 标准［6］的结果比较接近。通过有限元方法可以反映出齿轮啮合的时变刚度，清楚地显示齿轮在啮合过程中的啮合刚度较小值的位置，为齿轮的可靠性设计提供更有利的依据，而通过公式计算仅得出啮合刚度均值。

图10 啮合刚度曲线Fig.10 Mesh stiffness curves of gear during meshing

4 结论

1)通过与有限差分法计算齿轮轴变形量的对比分析，证明本文有限元方法以及所建模型的正确性、准确性，并能更真实地反应齿轮轴的变形。

2)分析了柔性齿轮轴和刚性齿轮轴两种有限元模型的齿轮载荷分布差异。柔性齿轮轴相对刚性齿轮轴齿宽接触力，齿面等效应力和齿根弯曲应力均呈现左侧增大，右端减小趋势，并且从等效应力分布图上可以明显看出左端接触面积增大，右端减小。这些分布差异可以为齿轮的设计提供更有利的参考。

3)由于齿轮轴和齿轮的变形，使齿轮重合度增大，时变啮合刚度减小。通过有限元方法计算的齿轮啮合刚度均值与ISO 6336-1 标准［16］中齿轮啮合刚度计算结果比较接近(当齿轮轴为柔性时，其刚度均值为15.235 4 N/(mm·μm). 刚性轴时，其刚度均值为18.343 1 N/(mm·μm). 通过ISO 6336-1 标准［16］计算的均值刚度为17.564 5 N/(mm·μm)，有限元方法能反映出齿轮啮合刚度的时变特性，并能清楚显示啮合刚度较小值的啮合位置，为齿轮的优化设计和可靠性设计提供更有利的设计依据。

References)

［1］张展. 渐开线圆柱齿轮传动［M］. 北京:机械工业出版社，2012.ZHANG Zhan. Involute cylindrical gear transmission［M］. Beijing:China Machine Press，2012. (in Chinese)

［2］唐进元，陈嫦，蒲太平. 考虑误差、轴变形和轴承刚度的弧齿锥齿轮LTCA［J］. 机械传动，2010，34(5):5 -9.TANG Jin-yuan，CHEN Chang，PU Tai-ping. LTCA analyses of spiral bevel gears with machining errors ，assembly errors ，axial deformation and bearing stiffness［J］. Journal of Mechanical Transmission，2010，34(5):5 -9.(in Chinese)

［3］Hwang S C，Lee J H，Lee D H，et al. Contact stress analysis for a pair of mating gears［J］. Mathematical and Computer Modelling，2013，57(1/2):40 -49.

［4］Gonzalez-Perez I，Iserte J L，Fuentes A. Implementation of Hertz theory and validation of a finite element model for stress analysis of gear drives with localized bearing contact［J］. Mechanism and Machine Theory，2011，46(6):765 -783.

［5］Patil S S，Karuppanan S，Atanasovska I，et al. Contact stress analysis of helical gear pairs，including frictional coefficients［J］. International Journal of Mechanical Sciences，2014，85(8):205 -211.

［6］李杰，张磊，赵旗，等. 变速器齿轮轴变形对齿轮接触状态的影响［J］.汽车工程，2012，34(2):138 -143.LI Jie，ZHANG Lei，ZHAO Qi，et al. The effects of gear shaft deformation on gear contact state in transmission［J］. Automotive Engineering，2012，34(2):138 -143. (in Chinese)

［7］殷昭云. 用有限差分法求解静定梁和静不定梁的挠度［J］.机械设计，1999，16(8):34 -35，50.YIN Zhao-yun. Solution of the deflection of statically determinate and statically indeterminate beams by means of finite difference method［J］. Journal of Machine Design，1999，16(8):34 -35，50. (in Chinese)

［8］刘鸿文. 材料力学［M］. 北京:高等教育出版社，1991.LIU Hong-wen. Mechanics of materials［M］. Beijing:Higher Education Press，1991. (in Chinese)

［9］Argyris J，Fuentes A，Litvin F L. Computerized integrated approach for design and stress analysis of spiral bevel gears［J］.Computer Methods in Applied Mechanics and Engineering，2002，191 (11/12):1057 -1095.

［10］于靖军. 机械原理［M］. 北京:机械工业出版社，2013.YU Jing-jun. Mechanisms and machine theory ［M］. Beijing:China Machine Press，2013. (in Chinese)

［11］顾守丰，连小珉，颜磊，等. 斜齿轮轮齿三维有限元网格自动生成及细化［J］. 清华大学学报，1996，36(8):77 -82.GU Shou-feng，LIAN Xiao-mian，YAN Lei，et al. Automatic generation and densification of 3-D finite element mesh fot helical gear tooth［J］. Journal of Tsinghua University，1996，36(8):77 -82. (in Chinese)

［12］唐进元，蒲太平. 基于有限元法的螺旋锥齿轮啮合刚度计算［J］.机械工程学报，2011，47(11):23 -29.TANG Jin-yuan，PU Tai-ping. Spiral bevel gear meshing stiffness calculations based on the finite element method［J］. Journal of Mechanical Engineering，2011，47(11):23 -29. (in Chinese)

［13］黄亚玲，秦大同，罗同云，等. 基于ANSYS 的斜齿轮接触非线性有限元分析［J］. 兵工学报，2006，27(4):31 -39.HUANG Ya-ling，QIN Da-tong，LUO Tong-yun，et al. The contact finite element analysis of helical gear based on ANSYS［J］. Acta Armamentarii，2006，27(4):31 -39. (in Chinese)

［14］Lin T J，Ou H，Li R F. A finite element method for 3D static and dynamic contact impact analysis of gear drives［J］. Computer Methods in Applied Mechanics and Engineering，2007，196(9/10/11/12):1716 -1728.

［15］朱增宝，朱如鹏，李应生，等. 啮合刚度对人字齿行星传动系统动态载荷特性的影响研究［J］. 兵工学报，2013，34(6):732 -738.ZHU Zeng-bao，ZHU Ru-peng，LI Ying-sheng，et al. Impact of mesh stiffness on dynamic load characteristics of herringbone planet train［J］. Acta Armamentarii，2013，34(6):732 -738.(in Chinese)

［16］ISO. ISO 6336-1 Calculation of load capacity of spur and helical gears-part 1:basic principles，introduction and general influence factors［S］. Geneva:ISO，2007.

［17］Amezketa M，Iriarte X，Ros J，et al. Dynamic model of a helical gear pair with backlash and angle-varying mesh stiffness［C］∥ECCOMAS Thematic Conference on Multibody Dynamics. Warsaw，Poland:Public University of Navarra，2009:1 -17.