螺旋锥齿轮是核心传动元件，广泛应用于车辆、飞机和船舶等装备。在传递动力的过程中，传动系统受载变形会导致齿轮副的啮合状态产生偏移，影响齿轮副的传动性能。通常要结合齿轮副的啮合错位结果开展齿面设计与修形。因此快速准确分析传动系统的啮合错位情况，对评估齿轮副的加载啮合性能和优化齿面设计非常重要。
国内外学者对齿轮传动系统的啮合错位量分析方法进行了大量的研究。王钦等建立了多支撑轴系耦合分析模型，综合考虑轴承刚度、齿轮、轴和壳体的影响求解齿轮啮合错位量，并通过实验进行验证。高洁等综合考虑了机匣、轴系及轴承的支撑变形对齿轮轴线偏移的影响对当量错位量叠加，在此基础上对弧齿锥齿轮齿面修形。覃秋霞等考虑了机匣的加工误差、轴承的过盈配合等因素建立啮合错位量计算模型，基于计算结果对弧齿锥齿轮齿面进行设计，用实验验证啮合错位后的印痕和传动误差。赖长发等借助 Masta 软件建立了驱动主减速器的仿真分析模型计算螺旋锥齿轮的啮合错位量，将仿真分析结果与实验作对比。然而为了快速获得计算结果，上述研究在分析传动轴的变形时，大多将传动轴简化为梁单元模型。但这种简化处理忽略了欧拉梁、铁木辛柯梁理论推导的假设条件，当传动轴结构复杂时，其变形分析结果可能会产生较大的误差，导致齿轮传动系统的啮合错位量计算不准确，从而影响齿面优化设计效果。
针对以上问题，本文提出一种基于实体有限元法的分步螺旋锥齿轮传动系统啮合错位分析方法。该方法考虑轴承刚度和齿轮啮合力的影响，采用实体有限元模型进行传动轴的变形分析保证齿轮副啮合错位量计算的准确性，以此为基础建立包含啮合错位影响的轮齿加载接触分析以大幅度提高计算分析效率。
一、齿轮副啮合错位量计算
啮合错位量的定义
齿轮传动系统包括轴、轴承和齿轮等零件，且结构复杂。传动系统工作时，齿轮副的啮合力传递到轴和轴承上，导致系统产生变形，从而使齿轮的相对安装位置发生偏移。齿轮副在载荷作用下安装位置的偏移量即为齿轮副的啮合错位量，如图 1 所示。
轴交错点变形前为 O、Q，变形后为A、B 。齿轮轴线变形前用细实线表示，理论变形后用点划线表示，ROA、ROB、ROQ 是点O 到 A、B、Q的向量。轴交错点的偏移矢量为
将其分解到初始齿轮坐标系中，则啮合错位量的计算公式为
式中：∆P为齿轮副沿小轮轴向的相对位移量;∆W 为沿大轮轴向的相对位移量;∆E为沿偏置方向的相对位移量;∆Σ为沿轴交角方向的相对变化量;i、j、k为初始坐标系中坐标轴的单位矢量;θ1、θ2 为变形前后的轴交角。
啮合错位量的计算方法
图 2 所示为齿轮副啮合错位量的计算原理图。初始坐标系Σq = {O; x, y,z} x、y、z为X、Y、Z轴的矢量;弯矩坐标系Σqy = {O; x ′, y ′,z ′} x ′、y ′、z为X ′、 Y ′、Z ′轴的矢量。
P、G 是变形前小、大齿轮齿宽中点在齿轮轴线上的投影点。受载变形后点 P、G 变成了点P ′、G ′ ，δP、δW分别是点 P到P ′ 以及G 到 G ′ 的距离。实际变形后轴线变成曲线用点划线表示，在点P ′、G ′处分别作两条曲线的切线，其方向向量为s1 和 s2。则两个 方向向量的法矢 n 为
设s3 是点 P ′ 和 G ′连线的矢量，设向量 s1和 s2之间的距离为lE 和角度为θS ，计算公式为
设小轮等效轴线上点A 与 P ′之间的距离为lP2，则建立包含 lP2 的方程表达式为
式中：ROP′为点O 和 P ′连线的矢量。
构建包含向量 s2 和法向量 n 的平面，其中点 A、B、G ′ 在该平面内，则RG′A 满足
式中：ROG′ 为点O 和G′ 连线的矢量;RG′A 为点G ′ 和 A 连线的矢量。
设大轮等效轴线上点B 与 G ′ 之间的距离为 lW2，计算公式为
将系统变形前后的参数叠加，获得齿轮副的啮合错位量。
式中：RQP 为点 Q和 P连线的矢量;ROG为点O 和 G 连线的矢量。
二、不同计算模型的传动轴变形分析对比
传动轴的受载变形是导致齿轮副啮合错位的主要影响因素，因此准确计算传动轴的变形量至关重要。传动轴大多简化为梁单元模型，其变形分析原理如图 3 所示。
根据传动轴的截面情况将其分段，用两个节点和一个梁单元表示相同截面的轴段。整根传动轴则可表示为由若干个不同截面的梁单元组合，其节点个数比单元个数多 1。每个节点包含 6 个自由度， 因此每个相同截面的梁单元可以用一个 12×12 的刚度矩阵表示。
将整根传动轴中每一段两节点的梁单元刚度矩阵进行叠加，得到整根传动轴的总刚度矩阵如下：
式中：K(n−1)n为 2×2 的刚度矩阵;矩阵 K(n−1)n包含
是节点 (n−1)和 n所在梁单元中在节点 n−1处的刚度;
为节点n处的刚度;A(n−1)n、An(n−1)为耦合刚度，且每个刚度都是 6×6 的矩阵，其余刚度以此类推。
根据力与变形之间的关系，构建梁的力平衡方程如下：
式中：Fn 为节点n处的载荷子矩阵;δn为变形子矩阵，两者都为 6×6 的矩阵。
已知外力和总刚度矩阵，通过平衡方程计算出梁单元在各节点上的变形量。
以某个复杂传动轴为例，首先根据传动轴的轴段数量来确定梁单元的初始节点数目，然后以载荷施加位置、边界条件及轴的变形提取位置进一步增加节点数量，最后通过这些节点去建立完整的梁单元模型如图 4 所示。节点 2 和节点 16 是支撑位置。节点 5～节点 9 是位移变形量的提取位置，节点 7 为载荷施加位置，其余节点是轴段的分界位置。
轴的材料为 20CrMnTi，泊松比为 0.3。在节点 7 处施加 3 个方向上的集中力，分别为 1 700、1 200、2 000 N。在节点 2 和节点 16 处添加边界条件，需要约束除了轴线回转方向之外的其余 5 个自由度。
为了对比不同模型分析轴变形的误差情况。同时按照图 4 的结构去构建轴的实体有限元模型，其材料、受力情况以及边界条件与梁单元模型的设置相同。通过静力学分析得到两种计算模型传动轴的变形如图 5 所示。
为了直观地对比两种模型的传动轴变形情况，提取轴线上节点的变形量，见表 1。
由表 1 可知，采用梁单元模型计算的传动轴变形量在各方向上与实体有限元模型相比均存在一定误差，径向方向(Xb、Yb)误差约 13%～16%，沿轴线方向(Zb)误差更大，约 30%。因此传动轴结构复杂时，不能用梁单元模型代替实体有限元模型去分析其受载变形情况。这种误差可能来源于梁单元模型受载变形后的梁截面与轴线关系的假设。
三、分步式啮合错位分析模型
综上所述，采用梁单元模型计算传动轴受载变形会在各方向上产生较大误差。为了提高变形分析的准确度，有学者研究通过实体有限元构建包含传动轴、齿轮以及轴承的传动系统全有限元模型，通过接触对定义齿面的啮合关系进行系统啮合错位的分析。这种方法能较好地提高计算精度，但是网格数量多，划分困难，非常耗费计算资源。为了保证计算分析的准确性，同时提高计算效率，本文提出一种基于实体有限元法的分步螺旋锥齿轮传动系统啮合错位分析方法：仅将传动轴构建为实体有限元模型，轴承等效为刚度矩阵耦合于传动轴上;将齿轮副传递的载荷等效为啮合力;采用静力学有限元法分析传动轴变形，进而计算齿轮副的啮合错位量;根据啮合错位量构建轮齿的加载接触分析模型分析齿面的接触偏移情况。
传动轴耦合分析模型
首先建立包含轴承影响的传动轴变形耦合分析模型，如图 6 所示。依据传动轴不同轴段去划分网格，为了减少网格数量，全部划分为六面体单元。再根据轴承安装位置、啮合力施加位置以及变形量提取位置对网格进行调整。
图 6 中轴段 1～5 为轴承安装位置，将这些轴段上轴承等效受力点处建立导向点，并将轴承与对应轴段上表面网格进行耦合。考虑轴承的弹性特征可以用弹簧来等效。每个导向点会通过 5 个弹簧约束其除回转方向之外的 5 个运动自由度，每个导向点在各方向弹簧的刚度由轴承的刚度矩阵中主对角元的数值确定。
轴段 6 和轴段 7 为大、小齿轮安装位置，在该轴段需要施加齿轮啮合力，同时需要提取传动轴受载变形量。通常将齿宽中部的节点位置定义为接触区的中心位置。若将啮合力施加在传动轴上时，会在齿宽中点 P 处切割出一个轴截面(P 点的位置如图 6 所示)，在轴心位置建立导向点并与截面耦合，并在该导向点上添加啮合力和弯矩。另外将该轴段等分切割成若干个子轴段，按照从左到右的顺序依次在轴心位置建立若干个导向点并与对应的轴截面耦合。通过以上处理，齿轮在安装位置的变形可由各导向点的位移表示。
等效啮合力计算
为了减少齿面接触分析的计算量，将齿轮接触作用等效为啮合力。根据传动扭矩计算节点处的切向力Fmt，考虑螺旋锥齿轮参数进一步计算出轴向力 Fax和径向力 Frd。以上计算的啮合力是作用在齿面上的，而 3.1 节所述的计算模型其啮合力作用于轴线上。力的作用点改变后需要采用力的平移定理重新计算载荷。其中径向力朝着轴心位置平移不产生弯矩，切向力产生绕轴线回转方向的弯矩不影响变形，轴向力产生其余方向弯矩会影响变形。
小轮坐标系 ΣP={Q; xP, yP, zP}与初始坐标系 Σq 在 z 方向的距离为偏置距，xP 是小轮回转轴，小轮轴向力产生 zP 方向的弯矩：
当齿轮副的偏置距为 0 时，大轮的弯矩计算方式与小轮相同。但偏置距不为 0 时，大轮轴向力产 生弯矩的方向与 z 方向有一定的角度，则在大轮轴心导向点处的弯矩坐标系Σqy 与啮合力坐标系Σq 不同。为了避免在同一个导向点施加啮合力和弯矩需要添加两个不同的坐标系问题，以弯矩坐标系Σqy 为主，将切向力和径向力分解到坐标系Σqy 中，计算公式为
式中：FX′、FZ' 为坐标系Σqy 中 x ′与 z'方向的力;MZ ′ 为 z'方向的弯矩。
在小轮坐标系 ΣP中齿宽中点所在轴截面的导向点处施加啮合力与弯矩;在大轮坐标系 Σqy 中齿宽中点所在轴截面的导向点处施加啮合力与弯矩。
考虑啮合错位的轮齿加载接触分析
将 3.2 节计算的等效啮合力代入 3.1 节所述的传动轴变形耦合分析模型，导出该模型中传动轴轴线上多个导向点的位移变形量，将其代入 1.2 节的公式求出系统的啮合错位量。通过啮合错位量去调整多齿啮合模型中齿轮的位置及姿态。首先将大轮沿着全局坐标系的Y 轴移动∆W，然后以大轮为基准将小轮分别沿着局部坐标系的 X1、Z1轴分别移动 ∆P、∆E，最后将小轮绕轴 Y1旋转 ∆Σ，从而构建出包含啮合错位的轮齿加载接触分析模型，如图 7 所示。
轮齿齿面通过接触对建立啮合作用。在大小齿轮宽中部的轴线中心建立载荷作用参考点，分别与齿轮的轮毂面建立耦合作用以便施加载荷。大轮参考点施加工作扭矩，小轮参考点施加转角，同时约束除齿轮回转自由度外的其余自由度。采用静力分析方法进行计算，则可获得某啮合瞬时齿面的接触状态以及齿轮的实际转角。改变小轮转角进行多步计算，则可得整个啮合过程中齿轮副的加载传动误差及接触印痕。
四、计算实例
啮合错位量计算对比
为了检验本文提出的实体有限元法分步螺旋锥齿轮传动系统啮合错位分析方法的准确性，以某客车驱动桥为例进行齿轮副啮合错位分析。其中螺旋锥齿轮基本参数见表 2。
齿轮副的额定工作扭矩为 15 000 N·m，利用 3.2 节所述的计算公式得到啮合力和弯矩计算结果，见表 3。
参考多支撑轴系中轴承单元刚度计算公式[11] 计算传动系统中轴承刚度矩阵，导出各轴承刚度矩阵中在主对角元上的数值，见表 4。
采用上述参数与传动轴图纸可以建立传动轴耦合分析有限元模型，计算获得传动轴的变形量及齿轮副的啮合错位量。同时建立图 8 所示的传动系统全有限元模型与本文方法做对比分析。由于齿轮与传动轴结合部位六面体网格划分困难，因此轮齿部分采用四面体进行网格划分。该模型中的材料属性、传动轴的设置方式与 3.1 节所述模型相同。齿面啮合作用参考 3.3 节接触对设置。添加边界条件时，需注意在施加工作扭矩和转角的参考点处约束自由度，除了轴线回转自由度外，其余自由度将由轴承的等效弹簧约束。
同时以梁单元构建传动系统模型进行系统的啮合错位分析，分别提取三种模型传动轴轴线上节点处的变形量。其中全有限元模型齿轮啮合位置是时刻变化的，因此提取五齿模型中第三齿在节点啮合时的变形量。最终绘制出传动轴的变形结果，如图 9 所示。
对比分析梁单元模型与全有限元模型的传动轴变形结果，其平均误差输入轴约为 12.4%、输出轴约为 10.2%。比较实体有限元传动轴耦合分析模型与全有限元模型的传动轴变形结果，其平均误差输入轴约为 1.9%、输出轴约为 1.2%。考虑到梁单元计算结果误差较大，将传动轴的变形量代入 1.2 节 的计算公式得到错位量，见表 5。
由表 5 可知，传动轴耦合分析模型与全有限元模型的啮合错位量误差在 2.5% 以内。
加载印痕和传动误差对比
将每一啮合瞬时的齿面接触位置的应力提出叠加到齿面上可获得齿面的加载接触印痕，本文所述方法计算的加载印痕如图 10a 所示，全有限元模型计算的加载印痕如图 10b 所示，两者的传动误差结果如图 11 所示。
由图 11 可知，分步螺旋锥齿轮传动系统啮合错位分析方法得到的加载印痕与全有限元模型计算的加载印痕位置及形状均一致，仅在齿宽方向有细微差异，其加载传动误差曲线形态吻合，幅值偏差在 1.3% 左右。
分步螺旋锥齿轮传动系统啮合错位分析方法分为两部分：传动轴耦合变形分析模型部分的网格数约为 18.5 万;考虑啮合错位的轮齿加载接触分析模型网格数约为 80.4 万(齿面接触部分网格大小为 0.2 mm)。全有限元模型网格数量约为 308 万(齿面接触部分网格大小为 0.4 mm)。统计两种分析模型的计算时间，见表 6。
由表 6 可知，分步螺旋锥齿轮传动系统啮合错位分析方法的分析时间相比全有限元模型大幅减少。若将全有限元模型的齿面网格划分为 0.2 mm，则分析时间还将进一步增加。
五、结语
针对传动轴结构复杂时采用梁单元模型求解轴的变形差异较大，会导致齿轮副啮合错位分析不准确的问题，本文提出一种基于实体有限元法的分步螺旋锥齿轮传动系统啮合错位分析方法。通过理论研究以及仿真分析得出以下结论：
(1)传动轴模型结构复杂时，采用梁单元模型求解轴的受载变形与实体有限元模型计算结果存在较大误差。
(2)提出的基于实体有限元模型的传动轴变形耦合分析模型，综合考虑了齿轮啮合力作用以及传动轴结构柔性与轴承支撑柔性耦合状况，能够更为准确地分析螺旋锥齿轮传动系统传动轴的受载变形情况。
(3)提出的基于实体有限元法的分步螺旋锥齿轮传动系统啮合错位分析方法能够准确分析齿轮副的加载印痕及传动误差，在保证计算精度的同时，计算效率明显提升。
参考文献略.