直齿行星传动齿轮修形优化设计方法

郝 东 升， 董 惠 敏＊， 毛 范 海， 张 永 帅

（1.大连理工大学 机械工程与材料能源学部，辽宁 大连 116024；2.大连重工·起重集团有限公司 齿轮箱研究所，辽宁 大连 116013）

0 引 言

行星齿轮传动由多个行星轮分担荷载，具有传动功率大、体积小、重量轻等优点，广泛应用于冶金矿山、水泥、火力发电等行业的动力传输设备.由于太阳轮功率密度高，多对齿轮副同时参与啮合，扭转变形使太阳轮齿宽方向荷载分布不均，齿轮受载变形和负载波动也会加剧齿轮啮入、啮出冲击，降低齿轮传动的可靠性.

齿轮修形能够改善齿轮啮合平稳性，减小齿轮受载变形和制造、装配误差引起的啮合冲击，改善齿面的润滑状态，获得较为均匀的齿面荷载分布.杨廷力等[1]定性给出了确定齿廓修形曲线、修形长度和修形量的一般原则；王统等[2]归纳了国内外计算齿廓修形和齿向修形的不少方法；李润方[3]以弹性啮合数值方法为基础，对齿轮修形方法进行了全面系统的论述.目前齿轮修形方法主要针对定轴直齿和斜齿传动，对行星传动齿轮修形研究较少.近年来有限元接触分析方法广泛应用于齿轮修形研究，但需要解决两个关键问题才能实现齿轮修形参数优化：（1）渐开线圆柱齿轮啮合理论上为线接触，接触线邻近区域应力梯度大，需要非常细密的网格划分来反映接触区的真实应力、接触变形和修形齿面的接触状态[4].为了减小齿轮有限元规模，提高非线性接触分析求解效率，通常采用疏密网格逐渐过渡方式划分网格，也可以通过多点约束方程实现疏密网格过渡有限元建模.文献[5，6]采用分级剖分法较好地实现了网格过渡，但对于直齿轮采用多点约束方程网格过渡建模，有限元模型的网格畸变小[7]，计算精度高.本文采用多点约束方程实现疏密网格过渡建模，提出具有齿面节点坐标修正的直齿轮接触区局部网格细化建模方法.（2）应用有限元接触分析优化齿轮修形参数，需要多次变更修形齿轮有限元模型，求解齿面上的应力分布，因此修形齿轮快速精确有限元建模成为齿轮修形参数优化的关键.

本文提出一种直齿轮修形齿面快速精确有限元建模方法：首先在ANSYS有限元分析平台上应用APDL语言建立局部网格细化的齿轮啮合有限元接触分析模型，施加荷载和位移边界条件；然后在此基础上按修形函数控制渐开线发生线长度修改齿面上的节点坐标，实现修形齿面快速精确重构.由于在一个啮合相位下只建立一个接触分析模型，修形参数迭代优化过程中的模型重建仅通过修形齿面重构实现，极大地提高了修形齿轮有限元建模效率.这种方法优势在于有限元模

型的建模效率高，修形区齿廓几何逼近精度高，修形齿面接触状态求解准确，有限元规模小，求解效率高.

对于高精度硬齿面低速重载的行星传动系统，由于齿轮的啮合频率与系统的低阶固有频率相差很远，动态激励不显著[8].在准静态条件下只有少数轮齿参与啮合，考虑到行星架和内齿圈通常具有较大刚度，因此只研究太阳轮和行星轮受载变形对齿面应力分布的影响.假设行星传动系统的各行星轮传递荷载相等，对太阳轮剖分取其中一路传动子系统模拟，退化为一对定轴齿轮传动，简化了行星传动系统有限元模型，但是需要建立太阳轮子域边界位移协调关系，模拟太阳轮受载的轴向扭转变形.

本文以某4 600kW立式磨机主减速机直齿行星传动为例，阐述行星传动系统行星轮修形参数优化过程.

1 行星传动系统支路有限元建模

立式磨机主减速机由三级齿轮传动组成，如图1所示，其中最后一级为5个行星轮均布的高精度硬齿面低速重载直齿行星传动.动力通过鼓形齿联轴器传递给浮动的太阳轮，行星架作为动力输出，径向轴承约束，内齿圈用键和螺栓固定在机架上.太阳轮输入转速105r/min，扭矩418 381 N·m，行星架转速24r/min，啮合频率35.6Hz，喷油强制润滑方式控制齿轮温升.行星传动系统齿轮参数见表1.

太阳轮齿数为27，行星轮数量为5，所以每个传动支路分担的齿数不为整数，无法应用循环对称约束建模求解.取太阳轮5齿局部模型与行星轮啮合，并对太阳轮两侧边界面上的节点对应网格划分，应用位移约束方程建立太阳轮两边界面的位移协调关系，近似模拟太阳轮受载后的轴向扭转变形（见图2）.

图1 主减速机传动系统Fig.1 Driven system of the main reducer

表1 行星传动系统齿轮参数Tab.1 Gear parameters of the planetary gear train

图2 太阳轮与行星轮啮合有限元模型Fig.2 FE model of a sun and planet gear mesh

旋转太阳轮所有节点的节点坐标系与全局圆柱坐标系方位一致，两边界上节点位移满足如下约束方程：

其中Hi和Lj分别为两边界上的节点对编号.

太阳轮内孔表层节点施加径向位移约束，输入端花键轴所在轴段表层节点施加切向荷载，其中边界面一侧棱边上的节点不施加荷载.

节点切向荷载计算：

式中：T为太阳轮输入扭矩；Nc为行星轮个数；d为简化的花键轴段直径；N为施加荷载的节点数.

在行星轮回转中心建立局部圆柱坐标系，旋转行星轮所有节点的节点坐标系与该局部圆柱坐标系方位一致.行星轮与内齿圈啮合的邻近区域节点施加切向位移约束，内环节点施加径向位移约束和轴向位移约束，防止行星轮产生刚体位移.

2 接触区细化单元边界位移耦合建模

齿轮啮合接触线邻近区域应力梯度大，对齿面网格划分敏感，需要高密度网格划分提高齿廓几何逼近精度，真实反映修形齿轮的齿面应力和接触状态.齿轮有限元详细建模过程参见文献[9]，

接触区细化单元边界位移耦合建模过程如下：

（1）在Matlab中模拟齿条型刀具加工齿轮，生成轮齿廓线上的离散点坐标.

（2）在ANSYS中应用APDL程序分段样条拟合齿廓坐标点，生成端面齿廓线.

（3）定义材料属性，选择单元类型.

（4）端面齿廓线区域划分，生成端面网格，拉伸生成六面体网格.

（5）旋转两齿轮到分析啮合相位.

（6）求解装配相位下啮合点坐标.

（7）按式（2）修改齿面上与啮合点最邻近节点的坐标为啮合点坐标（齿顶上节点不调整），按齿面上渐开线弧长坐标调整齿面上与啮合点相邻的节点间距.

（8）根据接触线上啮合点的节点编号获得相关联的单元，轮寻单元上的节点序列，分别以节点坐标生成Keypoints，并以此生成六面体实体；设置体上附着线的分段数和分段比例，映射网格生成细化六面体单元，最后删除原位置上的母单元，合并相同位置的节点编号，实现接触区单元细化，如图3所示.

图3 啮合区细化单元Fig.3 Refined elements near the mesh zone

（9）建立齿轮啮合接触对，定义接触分析实常数.

（10）选择细化单元边界上的节点和相邻大单元，执行CEINTF命令，按照节点在单元中的位置自动建立位移约束方程，按邻域容差移动节点到大单元表面上.

（11）施加荷载和位移边界条件.

由于细化单元按线性插值生成，单元的表层节点并不在设计齿面上，需要按式（2）修改齿面上的节点坐标重建精确齿面.

3 修形齿轮有限元快速精确建模方法

在局部网格细化的齿轮啮合有限元接触分析模型的基础上，按修形函数控制齿面上节点的渐开线发生线长度，修改齿面上的节点坐标，遍历齿面上的节点重构修形齿面，实现修形齿轮有限元快速精确建模.

设分析相位Φ，如图4所示，齿面上的节点识别过程如下：

检索获得轮齿上大于基圆半径的节点坐标.令z＝zi的齿轮轴截面上中心点O坐标为（x0，y0），任一节点p在齿轮局部坐标系下的坐标为（x1，y1），映射与p点具有相同半径rp的p′点到渐开线上.

p′点对应的压力角

基圆上渐开线起点A对应相位

B点对应相位

δ＝±1对应的齿廓如图4所示.

图4 轮齿节点映射Fig.4 Mapping of a node on a tooth profile

对于修形齿廓，渐开线的发生线长度L由修形函数确定.

（1）齿向修形（图5）

p′点坐标

式中：z为齿面上节点的轴向坐标；z0为齿轮端面的轴向坐标；lw为齿向修形长度；ew为齿向修形最大修形量；m为齿向修形函数指数；δw在z坐标正向增长时取1，反向取－1.

图5 齿向修形Fig.5 Lead modification of gear

（2）齿廓修形（图6）

式中：L为齿面上节点的渐开线发生线长度；Le为最大修形量所在位置的渐开线发生线长度；Ls为修形起点的渐开线发生线长度；l为齿廓修形所在区间啮合线上度量长度；ep为齿廓修形法向最大修形量；n为齿廓修形函数指数；δp在渐开线伸长方向取1，反向取－1.

图6 齿廓修形Fig.6 Profile modification of gear

齿廓修形参数工程标注计算如下：

（3）齿向齿廓综合修形

齿向修形和齿廓修形都是在节点的渐开线发生面内修改齿面上的节点坐标实现，因此可以统一齿向和齿廓修形算法，并且可以实现齿向齿廓综合修形建模.

由于接触区细化单元尺寸小，当一个节点上的总修形量较大时，直接沿渐开线发生线方向修改节点坐标会使单元形态变差，影响有限元模型求解精度，可以通过修正渐开线发生线向量角来避免这个问题.

令eb为齿面上一节点的齿向修形量，en为该点的齿廓修形量，在渐开线发生线方向上的总修形量

齿面上一节点K的渐开线发生线如图7所示，发生线长度

对式（7）两边取微分

由于总修形量在0.001～0.010mm量级，近似取微分量与变动量等价.

修形后齿面上节点的发生线向量角

修形后齿面上节点的发生线长度

图7 发生线向量角修正Fig.7 Vector angle correction of a generation line

发生线向量角调整引起的附加齿廓修形变动量可以通过对式（5）取微分获得：

通常齿廓修形函数指数1≤n≤2，mm量级，l＞1，｜den/dL｜ 1，因此发生线长度变动量dL引起发生线向量角变动，进而由齿廓修形函数引起的附加齿廓修形变动量den为高阶小量，式（2）忽略此项仍然具有很高精度.

4 齿轮修形参数优化

高精度硬齿面低速重载的齿轮传动系统通常以齿面应力分布均匀作为齿轮修形的设计目标，齿轮有限元面－面接触分析方法能够求解轮齿上的应力分布.本文的齿形优化方法按照齿面上的应力分布趋势变更修形函数，以齿面重构算法快速建立精确修形的齿轮有限元接触分析模型，迭代求解修形齿轮接触问题，直到获得理想的齿面应力分布，完成齿轮修形参数优化.修形参数优化流程如图8所示.

这种基于有限元结构修改的齿轮修形参数优化设计过程相当于对齿轮箱做负荷试验，每次检测齿面上的接触斑，按接触斑分布重新加工修形齿轮继续做负荷试验直到获得理想的接触斑为止.

4.1 齿轮修形基本参数确定

接触分析非线性求解时间较长，为了减少修形齿轮分析模型的建模和求解次数，提高齿轮修形参数优化效率，齿廓修形基本参数初值按文献[1]推荐值选取.

行星轮既参与外啮合又参与内啮合，为了简化修形齿轮加工工艺，仅对行星轮修形.高精度硬齿面修形齿轮通常以磨齿作为最后加工工序，采用数控成型磨齿机加工修形齿轮，对齿轮修形参数适应能力强.

图8 齿轮修形参数优化流程Fig.8 Flow chart for optimizing gear modification

直齿轮传动重合度较低，齿廓修形推荐为短修形，基本参数如下：双齿啮合区单侧长度Γ＝（ε－1）pb，其中基圆齿距pb＝πmcosα，端面重合度为ε.啮合线上修形长度l＝0.5Γ；齿廓修形函数指数n＝1.5.

（1）齿顶修形参数

（2）齿根修形参数

齿向修形参数与齿轮的结构和荷载密切相关，修形长度按ISO[10]推荐值取lw＝0.1b，b为齿轮的齿宽.

图9 齿轮啮合有限元模型Fig.9 FE model of gear meshing

图10 行星轮齿顶修形Fig.10 Modified profile at the tip of the planet gear

4.2 优化结果分析

应用ANSYS 10.0的增广拉格朗日法求解修形齿轮啮合接触问题，惩罚刚度系数取1，穿透系数取0.1，忽略齿面接触摩擦，采用高斯积分点作为接触探测点，关闭自动接触调整选项以避免修改齿面上的节点坐标.

图12 行星轮齿根修形Fig.12 Modified profile at the root of the planet gear

图13 行星轮齿顶啮合相位下额定荷载不修形和修形齿轮等效应力Fig.13 von Mises stress of the none modified and modified gear at rated load under meshing of the tip of the planet gear

建立行星轮齿顶啮合、单齿啮合和行星轮齿根啮合3个关键啮合相位下的齿轮有限元接触分析模型，如图9所示，求解修形齿轮的应力分布，以齿面应力分布状况评价齿轮修形效果.

齿轮修形参数优化是一个迭代过程，限于篇幅仅列出修形前的齿轮应力分布和修形参数优化后的齿轮应力分布，中间迭代结果不再展示.模型坐标系定义如图2所示，对应以上3个啮合相位的修形齿轮有限元模型局部放大图见图10～12，应力分布分别见图13～15.

图11 单齿啮合相位行星轮修形Fig.11 Modified planet gear at single tooth mesh

图14 单齿啮合相位下修形齿轮等效应力Fig.14 von Mises stress of the modified gear under meshing of single tooth

图15 行星轮齿根啮合下额定荷载不修形和修形齿轮等效应力Fig.15 von Mises stress of the none modified and modified gear at rated load under meshing of the root of the planet gear

齿轮修形量的下限由齿轮的弹性变形量确定，上限由磨齿机的加工精度保证，修形齿轮优化参数如图16所示，参考方向见图1.

考虑到磨机工作时单向运转，所以行星轮两侧齿面具有不同修形参数.内啮合齿向变形比外啮合小，仅做齿端修薄，避免边缘接触，齿廓修形两侧取相同参数.修形前后应力对比（表2）表明：修形后齿面上最大应力值显著降低，齿面上应力分布得到明显改善，如图13～15所示；应力对荷载波动适应性增强，1.2倍额定荷载下应力为331 MPa，应力比为1.05.

图16 齿轮修形参数Fig.16 Parameters of the modified gear

表2 额定荷载下最大等效应力比较Tab.2 Comparison of max von Mises stress at rated load

5 结 论

（1）针对多个行星轮均布的直齿行星传动系统，应用位移约束方程建立太阳轮子域边界位移协调关系，以支路模型等效轮系整体模型，简化了行星传动系统有限元模型.

（2）提出具有齿面节点坐标修正的直齿轮接触区网格局部细化建模方法.通过单元形函数建立细化单元边界上的节点位移与相邻大单元节点位移之间的约束方程实现疏密网格过渡，既提高了接触区网格的齿廓几何逼近精度，又降低了齿轮有限元规模，为准确模拟直齿轮弹性啮合齿面的接触状态提供了新方法.

（3）按齿面重构算法修改齿面上的节点坐标，统一了齿向修形、齿廓修形和综合修形的修形齿面建模方法，实现了修形齿轮有限元接触分析快速精确建模，解决了齿轮修形参数优化过程中修形齿轮有限元模型更新这一关键问题.

[1] 杨廷力，叶 新，王玉璞，等.渐开线高速齿轮的齿高修形（一）[J].齿轮，1982（3）：14－24.YANG Ting－li，YE Xin， WANG Yu－pu，etal.Profile modification of high speed involute gears（1）[J].Gear，1982（3）：14－24.（in Chinese）

[2] 王 统，贾 毅，邱良恒，等.渐开线齿轮修形方法的进展[J].上海交通大学学报，1998，32（5）：133－137.WANG Tong，JIA Yi，QIU Liang－heng，etal.Progress outline of gear tooth profile modification[J].Journal of Shanghai Jiaotong University，1998，32（5）：133－137.（in Chinese）

[3] 李润方.齿轮传动的刚度分析和修形方法[M].重庆：重庆大学出版社，1998.LI Run－fang.Stiffness Analysis and Modification Method on Gear Transmission [M].Chongqing：Chongqing University Press，1998.（in Chinese）

[4] Coy J J，CHAO C H.A method of selecting grid size to account for Hertz deformation in finite element analysis of spur gears [J].Journal of Mechanical Design，1982，104（4）：759－764.

[5] 王文杰.风机增速器斜齿轮啮合载荷及齿向修形的研究[D].大连：大连理工大学，2008.WANG Wen－jie. The engaging load and axial modification of helical gear used in increasing gearbox of wind turbine generator [D].Dalian：Dalian University of Technology，2008.（in Chinese）

[6] Moriwaki I，Fukuda T，Watabe Y，etal.Global local finite element method （GLFEM）in gear tooth stress analysis[J].Transactions of ASEM，Journal of Mechanical Design，1993，115（4）：1008－1012.

[7] Mao K. Gear tooth contact analysis and its application in the reduction of fatigue wear [J].Wear，2007，262（11－12）：1281－1288.

[8] Parker R G，Agashe V，Vijayakar S M.Dynamic response of a planetary gear system using a finite element/contact mechanics model [J].Journal of Mechanical Design，2000，122（3）：304－310.

[9] 郝东升，王德伦.斜齿轮精确接触分析有限元建模方法[J].大连理工大学学报，2011，51（6）：825－831.HAO Dong－sheng，WANG De－lun.An approach to build a FEM of helical gears with accuracy contact analysis [J].Journal of Dalian University of Technology，2011，51（6）：825－831.

[10] ISO.ISO 6336－1－2007Calculation of Load Capacity of Spur and Helical Gears：Basic Principle and Influence Factor[S].Geneva：ISO，1990.