轮齿承载能力是评价齿轮传动品质的一项重要标准，齿间载荷分配系数和面接触应力的模型精度是齿轮承载接触力学计算结果精确性的前提与关键。齿轮修形可以有效改善齿轮动载荷变化梯度，减少冲击、振动与噪音，对提高齿轮传动品质有十分重要的作用。然而，已有的齿轮承载接触分析力学解析模型大都只能考虑鼓形的齿向修形，较难计入齿廓修形的影响。
对于齿轮承载过程中产生的齿面接触应力，学者常采用有限元法和解析法。其中，有限元法建模精度高，且计算结果精确，因而得到较为广泛的应用。朱才朝利用基于势能法的三段函数式齿间载荷分配系数模型分析齿轮的弹流润滑特性; AZNAR 等采用 MPC 算法与局部网格细化法优化齿轮的承载接触有限元模型，提高了计算精度与效率; MAPER 等建立了包含齿廓修形的渐开线直齿圆柱齿轮承载接触有限元模型; 唐进元等利用有限元分析了正交面齿轮的重合度、载荷分布、接触应力等承载接触性能参数; 白恩军等建立了斜齿圆柱齿轮承载接触有限元模型。然而，有限元算法计算精度依赖于模型的网格质量和局部网格细化程度，存在计算效率低、收敛难的问题。
相较于有限元法，基于赫兹弹性接触理论的力学解析算法，计算快捷高效受到了越来越多的关注与研究。王会良等基于齿轮 TCA 法和齿面柔度系数，分析齿轮承载接触特性，建立了齿轮 LTCA 理论; 曹雪梅等提出一种简化轮齿接触分析中非线性方程数量的分解算法，提高了计算效率; 王羽达等基于轮齿 TCA 模型和赫兹接触理论，计算了渐开线齿轮的动态接触应力。
综上所述，传统齿间载荷分配系数解析建模未考虑齿廓修形及其引起的单双齿接触区域变化，无法准确计算承载啮合齿对的齿间载荷分配系数与齿面接触应力，同时，LTCA 方法大都仅在 TCA 分析阶段考虑齿廓修形，而在齿面柔度系数或承载接触分析时采用有限元方法。为此，本文以齿廓修形的渐开线直齿轮为研究对象，提出修形齿轮单双齿啮合状态的 s 判定方法与相应的啮合刚度计算方法，建立更为精确的修形齿轮齿间载荷分配系数与齿面接触应力解析算法，对不同修形高度、修形曲线与配对齿数条件下的齿轮算例进行承载接触分析，并与同参数下的有限元结果进行对比，验证解析算法的精度与稳定性。
一、齿廓修形的渐开线直齿轮模型
为建立齿廓修形渐开线直齿轮的齿廓曲线方程，需确定齿廓最大修形量 Δmaxi、修形高度 h、修形曲线，根据文献最大修形量 Δmaxi为:
式中: Δmaxi、ci 分别为最大修形量与修形修正量(i 为 1、2，用以区分参与啮合的主/从动轮) ，相应的主/从动轮修形修正量 c1、c2 分别取 9 与 4 ，Ft 为接触点的圆周力，B 为齿轮宽度。
建立坐标系 XdOYd，如图 1 所示，以 Δx 作为增量等距离散啮合线，并在圆周方向上计算出轮齿 1 上的齿廓对应点，若该点位于修形齿廓区域内，啮合线离散点的齿廓对应点处的修形量 Δj 为:
式中: 修形曲线变化系数 k 由选取的修形曲线确定，Linear、Yoshio、Walker 与抛物线修形曲线变化系数 k 分别为 1.0、1.2、1.5 与 2.0，x 为啮合点 M 在啮合线方向上距离单双齿变化点 A 的距离，L 为轮齿 2 的修形区域在高度方向上对应至啮合线 A 点至 C 点的线段长度( 长修形) ，C 点为轮齿 1 的啮合起始点; 另外，根据文献，正弦修形曲线的啮合点修形量为:
确定齿廓对应点处的修形量 Δj 后，根据渐开线齿廓曲线展成原理，可建立修形齿廓曲线参数方程为:
式中: u 为啮合线离散点对应渐开线齿廓处的参角，u = θ + αt，θ 为啮合点处展角，αt 为压力角。
二、修形齿轮单双齿啮合状态判定与啮合刚度
标准渐开线齿轮副在承载啮合传动过程中齿轮啮合刚度主要包括 3 部分: 轮齿刚度(弯曲刚度 kb、剪切刚度 ks 与径向压缩刚度 ka ) 、齿轮接触刚度 kc 与齿轮基体刚度 kf。单齿啮合刚度 Ki 为上述各刚度的串联形式，双齿啮合刚度 Kall为两对轮齿啮合刚度的并联形式。
式中: i 为 1、2，用于区分参与啮合的不同齿对，下标 g、e 分别表示主动轮与从动轮。
当采用齿廓修形时，原单齿啮合区间的齿轮啮合刚度根据式(5) 进行计算，原双齿啮合区间由于修形齿廓参与啮合，啮合状态会发生变化，需要重新判定啮合状态: 根据齿对啮合力与啮合刚度计算接触点处沿啮合线方向的变形量，然后与接触点处修形量进行对比，基于二者的大小关系判定修形齿轮副的单双齿啮合状态，并计算单双齿的啮合区间与啮合刚度，其中，修形齿轮啮合刚度定义为 KTPM。
如图 2a 所示，当主动轮顺时针转动时，齿对 1 首先进入啮合承担载荷，由于齿对 2 主动轮处存在修形齿廓，当齿对 1 沿啮合线方向的变形量 δ1 大于齿对 2 处主动轮的齿廓修形量 Δ1 时，齿对 2 将进入啮合并参与传动，此时两对轮齿共同承担负载扭矩 Tc，并在啮合线方向上的两处作用点产生啮合力 F1、F2，根据力矩平衡和广义胡克定律可得:
式中: K1、K2 与 δ1、δ2 分别为齿对 1、齿对 2 的单齿啮合刚度与接触点沿啮合线方向的变形量，rb2为从动轮的基圆半径。
与此同时，齿对 1 沿啮合线方向的变形量 δ1 即为整体啮合变形量 δall，则 δall = δ1 > δ2，相应的修形齿轮啮合刚度 KTPM为:
式中: Fn 为齿轮副法向力，根据式(7) 将式(8) 推导为:
同时，齿对 2 在啮合点处啮合变形量 δ2 与修形量 Δj 之和等于齿对 1 处啮合变形量为 δ1。
将式(10) 带入式(9) ，可得到:
将式(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) 计算双齿区间内的齿轮啮合刚度。
如图 3 所示，基于主动轮齿廓，对齿轮副啮合区间进行划分，分析齿对 1、2 的接触啮合状态及相应的啮合刚度，点 a1 为啮入点，点 a3、a4 为单双齿啮合变化点，点 a2、a5 为原双齿接触区域平分点，点 a1 ～ 点 a5 将接触区域划分为 A1、A2、A3、A4、A5。
如图 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 所述，即可得到齿廓修形齿轮在各啮合位置的单双齿啮合状态判定与相应啮合刚度计算方法。
三、齿廓修形渐开线直齿轮承载接触分析
齿间载荷分配系数
未修形齿轮的齿间载荷分配系数：根据式(6) 确定的双齿啮合刚度 Kall与齿轮整体变形量 δall，法向力 Fn 可表示为:
根据式(7) 中确定的载荷、刚度以及变形关系，带入式(13) 等式左侧，可推导出:
由于啮合齿对在啮合线方向上变形量 δ1、δ2 中的较大值与双齿啮合整体变形量 δall相等，此时无论δall = δ1≥δ2 或 δall = δ2≥δ1，根据式(14) 皆可推导出各啮合齿对在啮合线方向上的变形量 δ1、δ2 与双齿整体啮合变形量 δall相等，即:
进一步可根据式(7) 中啮合力、刚度与变形的数学关系推导出:
式中: Fi 为接触齿对的啮合力( i 为 1、2) 。
对式(16) 进行比例互换，可得到未修形下的齿轮齿间载荷分配系数 S 为:
齿廓修形齿轮的齿间载荷分配系数：当齿廓修形后，原双齿啮合区域的啮合状态发生变化，根据标准渐开线齿轮时变啮合刚度式(5) 与齿廓修形齿轮啮合刚度式(12) ，重新确定齿轮的单双齿啮合状态与相应啮合刚度，继而根据式(17) 计算修形齿轮不同齿廓区域参与啮合的齿间载荷分配系数 SLDF。根据图 3，以齿对 1 为例，分析齿对 1 承载啮合下的齿间载荷分配情况，具体为:
式中: Ku 为齿对 1 的啮合刚度，当齿对 1 处于单齿啮合状态时，其啮合刚度依据式(5) 计算，即 Ku = K1 ; 当齿对 1、2 同时参与啮合时，此时计算 Ku 需考虑齿对 2 啮合刚度 K2 的叠加影响，即:
当齿对 1 在齿廓区域 A1 处啮合时，由于齿对 2 接触变形量 δ2 的影响，齿对 1 存在非接触与接触两种状态，若 δ2 > Δ2，则齿对 1 接触，此时为双齿啮合状态，将式(18) 代入式(19) 可得到齿对 1 的齿间载荷分配系数 SLDF。
若 δ2 < Δ2，则齿对 1 为非接触状态，此时齿对 1 不分担载荷，则齿间载荷分配系数为 0，即:
当齿对 1 在齿廓区域 A2 处啮合时，由于齿对 1 接触变形量 δ1 的影响，齿对 2 存在非接触与接触两种啮合状态，若 δ1 > Δ1，则齿对 1 接触，此时为双齿啮合状态，齿对 1 的啮合齿廓为非修形区域，其齿间载荷分配系数根据式(22) 进行计算。
若 δ1 < Δ1，则齿对 2 为非接触，此时齿对 1 承担全部载荷，齿间载荷分配系数为 1，即:
当齿对 1 在齿廓 A3 处啮合时，此时为单齿啮合区域，齿间载荷分配系数为 1，即:
当齿廓进入 A4、A5 段参与啮合时，齿对 1 皆为非修形齿廓参与啮合，此时仅考虑齿对 1 的啮合刚度，即 Ku = K1，相应的齿间载荷分配系数为:
综上所述，在一个啮合周期内，根据齿对 1 在齿廓区域 A1、A2、A3、A4、A5 是否接触参与啮合及其啮合状态判定，计算参与啮合齿对的啮合刚度，根据式(18) ～ 式(25) ，可得到修形渐开线直齿轮齿间载荷分配系数 SLDF计算模型。
齿面接触应力
根据上述建立的修形齿轮齿间载荷分配系数模型，可计算啮合力 Fi，继而根据赫兹接触理论计算齿面接触应力。
首先，需确定不同齿廓赫兹接触类型下接触点处的综合曲率半径 Re，针对未修形的齿廓区域，根据欧拉萨瓦利公式，通过计算沿啮合线方向上基圆至啮合点的距离求得其曲率半径 ρ1、ρ2，对于齿廓修形区域，以主动轮为例，根据微分几何与修形齿廓方程式 (4) ，可推导出修形齿廓曲率 CTPM与曲率半径 ρTPM。
针对齿廓修形齿轮承载啮合时存在的不同赫兹接触类型，可推导出相应的综合曲率半径 Re。
通过式(26) 确定的齿间载荷分配系数 SLDF，可计算出参与啮合齿对的接触力 Fi，根据式(29) 所得综合曲率半径 Re，基于赫兹弹性接触理论，可推导出修形齿轮承载接触时的齿面接触应力 σH 与接触半宽 b。
式中: Ee 为综合弹性模量，具体表达式为:
式中: E1、E2 与 ν1、ν2 分别为主、从动轮的弹性模量与泊松比。
如图 4 所示为该解析算法的流程图。
四、计算结果分析与讨论
基于上述建立的算法，求解不同修形高度、修形曲线与配对齿数 3 种条件下不同算例的齿间载荷分配系数与齿面接触应力，并与同参数下的有限元计算结果进行对比评价，验证本文解析算法计算结果的准确性和稳定性。基于表 1 中抛物线修形齿轮的主要几何参数，每种条件下选取 5 个算例求解以上不同齿廓修形齿轮算例的齿间载荷分配系数(图 5 ～ 图 7) 和齿面接触应力(图 12 ～ 图 14) 。使用有限元法对比评价解析 算法结果的准确性和稳定性，以验证算法的正确性。
分析如图 5 ～ 图 7 所示的齿间载荷分配系数结果可知:
(1) 本文解析算法与有限元两种计算方法的齿间载荷分配系数具有相同变化趋势: 由最小值逐渐增大至最大值 1，此过程为双齿啮合区间; 随后进入单齿啮合区间，即 SLDF = 1 的水平线; 再进入双齿啮合区，下降至最小值。同时啮合线各点所对应的误差数值较小。
(2) 根据图 5 ～ 图 7 数据结果，可得到不同齿廓修形算例在一个啮合区间内的齿间载荷分配系数误差均值，如图 8 所示，每种条件下 5 个算例误差均值变化范围分别为 7.5% ～ 9.2%、7.5% ～ 8.0%、5.6% ～ 7.8%， 表明误差均值较小且不同算例的误差均值相近。算例结果分析表明，对于不同算例情况下的齿轮副，本文解析算法的齿间载荷分配系数计算精度较高，且算法具有较好的误差稳定性。
(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 种条件下不同齿廓修形算例在一个啮合区间内的齿面接触应力结果可知:
(1) 解析算法的齿面接触应力结果与有限元结果在变化趋势具有较好的一致性; 根据图 12 ～ 图 14 的齿面接触应力结果，可得到 3 种条件下齿廓修形算例在一个啮合区间内的齿面接触应力误差均值( 图 15a) 变化范围分别为7.0% ～ 9.0%、8.4% ～ 10.4%、9.0% ～ 10.7%，误差均值较小，同时不同算例间误差均值相近，表明解析算法具有较高的齿面接触应力计算精度与误差稳定性。
(2) 与此同时，在修形齿轮刚开始进入啮合、即将退出啮合以及修形曲线与渐开线过渡区域，有限元结果产生了局部应力突变，以图 12 ～ 图 14 中的 A、B 区域为例，主要由以下两个原因导致: ①刚开始啮入与即将退出啮合时，接触齿对啮合力较小，需要极小网格尺寸( 低于微米 1 ～2 个量级) 才能有效捕捉到该处的接触应力信息，建模困难，且啮入时齿顶部位存在几何尖点; ②修形齿廓过渡区域对几何曲率产生一定程度的啮合畸变。因此，当求解步长位于上述位置附近时会直接影响有限元计算结果，导致有限元局部应力突变、计算误差较大，这是不可避免的。除去由有限元结果突变区域所带来的计算误差后，解析算法的齿面接触应力在一个啮合区间内的误差均值(图 15b) 变化范围分别为 4.8% ～ 6.4% 、4.5% ～ 6.8% 、6.1% ～ 8.6% ，再次表明解析算法具有较高的齿面接触应力计算精度与误差稳定性。
五、结论
(1) 本文针对齿廓修形的渐开线直齿轮，通过对比承载齿轮接触点沿啮合线方向的变形量与修形量，得到了修形齿轮单双齿啮合状态的判定方法与相应啮合刚度的计算方法，通过推导修形齿轮啮合力、刚度与变形的力学关系，建立了修形齿轮各啮合齿对精确的 齿间载荷分配系数计算模型。
(2) 基于齿间载荷分配系数计算了啮合齿对实际承担载荷，应用欧拉萨瓦利公式和赫兹接触理论，建立了齿廓修形齿轮啮合过程中的齿面接触应力模型。
(3) 不同修形高度、修形曲线与配对齿数条件下多种齿廓修形齿轮的算例结果表明: 本文所提出解析算法的齿间载荷分配系数、齿面接触应力与相应单双齿啮合区间结果与有限元方法结果变化规律相一致，相对误差较小，具有较高的计算精度与稳定性，为齿廓修形渐开线直齿轮的承载接触分析提供了一种高效便捷的新方法。
参考文献略.