基于齿面网格节点变形量的齿廓修形计算方法

李剑敏 吴跃成 宣海枫 陈文华 潘晓东 王 威 周新腾

1.浙江理工大学浙江省机电产品可靠性重点实验室，杭州，3100182.杭州前进齿轮箱集团股份有限公司，杭州，310027

基于齿面网格节点变形量的齿廓修形计算方法

李剑敏1吴跃成1宣海枫1陈文华1潘晓东2王 威1周新腾1

1.浙江理工大学浙江省机电产品可靠性重点实验室，杭州，3100182.杭州前进齿轮箱集团股份有限公司，杭州，310027

针对重载齿轮修形计算问题，提出了按照齿顶、齿面在载荷作用下的变形量进行修形的方法，以有限元动力学计算结果为依据，提取轮齿相关节点的位移，处理得到其弹性变形，作为齿轮最大修形量进行修形；编制了计算程序，实现了齿轮修形轮廓线的自动计算；某风电齿轮箱三级齿轮经该方法修形后，齿轮最大应力减小约20%，延长了齿轮寿命，提高了齿轮可靠性。

齿轮修形；有限元分析；节点变形；应力

0 引言

齿轮是机械传动的主要方式，在风电、航空、舰船、车辆等领域都有着极为广泛的应用。经典齿轮齿面轮廓一般按渐开线设计，但在啮入/啮出阶段的冲击增大了轮齿表面应力；啮合过程中,轮齿受载荷作用而变形，使得齿面曲线发生变化，也造成啮合应力有所增大。以上原因造成齿轮啮合应力比预期的值要大，从而使齿轮的寿命缩短，可靠性降低。为延长齿轮寿命，提高其可靠性，通过齿轮修形来减小和均衡齿面接触应力成为目前主要的技术手段[1-3]。

目前在齿轮设计研究过程中，企业通常按照产品技术、工艺要求设计齿轮，然后根据齿轮的形制、载荷、材料等技术参数，进行轮齿修形[4-6]，但在修形计算中，没有对轮齿的表面应力、变形状况进行分析，而修形与轮齿的啮合应力、变形等直接相关，因此，该修形方式不能得到较好的修形结果。近年来，基于有限元分析技术可进行齿轮啮合过程的数值仿真[7-11]，文献[12-14]提出了基于有限元数值仿真结果的齿轮修形计算方法，但在上述修形计算中，齿轮最大修形量的取值及其物理含义尚不够明确，大量的计算数据不能从有限元结果中自动提取，数据交换比较复杂，需要大量的人工参与，在齿轮设计生产企业进行推广有一定困难。

本文提出了一种基于有限元结果的齿轮修形量计算方法，能够在有限元分析基础上实现修形量、齿面轮廓的自动计算。该修形方法使用ANSYS有限元程序的APDL语言生成齿轮三维模型，分析得到轮齿的应力、应变、位移等数值结果；根据上述结果计算轮齿修形量，并根据修形量移动模型节点坐标形成新的轮齿曲面。对某风电齿轮箱第三级齿轮副进行修形前后对比计算，结果表明齿轮啮合最大应力减小约20%。

1 齿轮三维建模

由机械原理可知，渐开线在极坐标系和直角坐标系下的方程分别如下[15]：

(1)

(2)

式中，rb为基圆半径；rk为渐开线在任意点k的向径；αk为渐开线在k点的压力角；θ、θk为展角。

齿轮模型除了渐开线建模之外，根据实际加工情况，需建立齿根过渡曲线。值得注意的是，过渡曲线对于齿轮弯曲强度的求解有重要意义。根据普通齿条刀具切割过程，在直角坐标系下，可得到渐开线过渡曲线的方程[16-18]。

2 齿轮修形设计

齿轮运行中的危险应力一般发生在齿面刚进入啮合时，主要原因是啮合过程中会因碰撞(冲击)、挤压而产生弹塑性变形。根据齿轮啮合工况设计，当一个齿面刚啮入时，啮合处两齿面相切是最佳啮合设计工况，也就是渐开线齿轮的基本要求。但当齿轮受力变形时，其齿面轮廓发生了一定的变化，从补偿变形角度来看，本文齿轮的最大修形量ep直接取齿顶啮入变形量。但实际分析中，齿顶变形较难直接得到，因为在齿轮运行过程中，轮齿整体的转动与受力引起的弹性变形耦合在一起，而修形量仅与其弹性变形相关，因此，需要从有限元分析所提取的位移中剔除轮齿转动效应，从而得到轮齿的“真实”变形。通常齿轮受载荷引起的轮齿变形相对于齿轮整体转动引起的位置变化是很小的，因此，在分析中，不考虑转动与弹性变形之间的相互耦合，仅认为在提取的位移中，既包含了转动位移，也包含了弹性变形位移。转动位移视为刚性位移，可以按照刚体转动理论进行计算，并从总位移中扣除，得到所需要的相对弹性位移。

如图1所示，设A、C两点为齿轮同一横截面上两个点，其中，A位于齿根处，C为齿面上的点。修形量计算需要得到C点相对于A点的弹性变形。设A、C点的坐标分别为(xA,yA,zA)和(xC,yC,zC)。当齿轮运转之后，设齿轮转动了角θ，节点A的坐标变为A4(xA+ΔxA,yA+ΔyA,zA+ΔzA)。A4在oxy面上的投影为A2(xA+ΔxA,yA+ΔyA,0)，其中，θ=∠AoA2=∠CoC1，其计算公式如下：

(3)

式中，a为中心距。

图1 向量在oxy平面内的位移计算Fig.1 Displacement of the vector on oxy plane

点A、C在oxy面内的投影，仅旋转θ角之后所得节点为A1和C1。点A1、C1坐标分量为

(4)

(5)

通过平移向量A1C1，即可获得轮齿上C点在齿面上的相对位移。由于在齿轮啮合动力学分析中，各节点沿z轴方向的位移都很小(1μm以内)，故C点在z轴方向的理论位移直接取ΔzA，至此即可获得C点刚性转动角时的理论位置C3(xA+ΔxA+xC1-xA1,yA+ΔyA+yC1-yA1,0)，如图2所示。而节点C可通过有限元动力学分析，得到包含弹性变形与刚性转动的综合位移，即其坐标值为C4(xC+ΔxC,yC+ΔyC,zC+ΔzC)，节点C4与C3的距离即为点C的弹性变形量。在修形量计算中，通过有限元软件ANSYS的APDL语言获得齿面上相应节点的弹性变形，作为齿廓最

图2 节点C理论未变形旋转位置Fig.2 Location place on rotation about node C

大修形量ep和齿向最大修形量e1、e2的计算参考数值。

齿廓修形长度主要分为长修形和短修形。一般可根据基节长度和重合的关系获得，即

λ=pb(εα-1)

(6)

式中，pb为齿轮基节长度；εα为齿轮端面重合度。

齿廓修形的几何形态如图3所示，其齿廓沿渐开线发生线方向上的修形量en可按下式确定[13,19]：

(7)

其中，ep可以取齿顶最大相对变形量；Le是最大修形量所在位置的渐开线发生线长度；Ls是修形起始的渐开线发生线的长度；n为齿廓修形函数指数(1

图3 齿廓修形Fig.3 Tooth profile modification

对于齿向修形，主要考虑修鼓半径、齿端最大修形量。本文考虑到修鼓量沿渐开线方向的变化情况，将修形量递减曲线和径向修形长度两种因素考虑进来。其中，修形量递减曲线为齿端最大修形量随渐开线长度变化而变化的曲线；径向修形长度为齿向修形时在渐开线方向上的修形长度。

图4所示为渐开线齿廓发生面内的鼓形修形，齿面上各节点可通过下式确定[13]：

(8)

(9)

(10)

式中，e1、e2为齿端最大修形量；b1、b2为齿廓表面鼓形起始位置到齿端的距离(图4)；zi为各节点所在z坐标；ei为各节点齿向修形量。

图4 齿廓发生面内鼓形修形Fig4 Tooth drum modification in plane

本文以上述公式为基础，编制了齿轮和修形的计算程序，实现了齿轮修形的自动计算，其程序流程如图5所示。在ANSYS有限元软件里进行齿轮建模，使用APDL提取齿面各节点编号及其三维坐标；根据齿轮运行工况进行瞬态动力学计算，得到齿面啮合应力、变形的分布规律；根据有限元的计算结果，对齿面应力较大的啮合处、齿顶等计算修形量；由计算得到的修形量，转换坐标，进行节点移位，并形成修形后的轮齿轮廓线；对新的齿轮进行有限元计算，验证修形前后啮合应力的变化，完成齿轮修形分析过程。

图5 齿轮建模分析流程图Fig.5 Flow diagram of gear modeling

3 齿轮修形分析实例

本文以杭州前进齿轮箱股份有限公司的某型风电齿轮箱第三级齿轮对为研究对象，其主要参数见表1和表2。

表1 齿轮基本参数Tab.1 Basic parameters of the gear

表2 齿轮材料性质参数[20]

将未修形的模型导入ANSYS中进行齿轮瞬态动力学分析，齿轮材料为20CrMnTi，主动轮转速为275.97 r/min、齿轮输出力矩为9160 N·m，其装配模型如图6所示。

图6 齿轮模型以及约束定义Fig.6 Gears and constraint

根据计算结果，获得某齿面在设计载荷工况下的啮合应力、应变分布情况。当仿真时间为0.1739 s时，齿顶进入啮合状态，提取此刻最大应力节点编号，以此节点作为起始节点分别沿轴向和齿面渐开线方向设计路径，获得啮合路径上的应力、应变分布情况，如图7和图8所示。从图7可知，轴向路径上最大应力所在节点为z=0处节点，即齿顶外缘的齿端处，其应力值为443.8 MPa；最小应力在z=46 mm处，即齿顶中心位置处，其应力值为351.717 MPa。从图8可知，在渐开线方向上，应力最大处即齿顶处，其应力值为443.8 MPa；最小应力值为63.45 MPa，由于弯曲应力的作用，沿齿根方向应力逐渐增大。通过分析即可得到齿顶啮合时危险区域最大应力所在位置及其应力值。

图7 齿轮轴向啮合路径上的应力分布(0.1739 s时)Fig.7 Stress on path of engage axial direction of gears when time is 0.1739 s

图8 啮合齿轮沿渐开线方向应力分布(0.1739 s时)Fig.8 Stress on path of evolvent direction of gears when time is 0.1739 s

由式(3)～式(5)所示节点弹性变形量的计算方法，可以直接计算得到随时间变化的齿顶两端面危险节点和齿顶中间节点的弹性变形量。当齿顶中心节点处于啮合状态时，其变形量λ1=6 μm；齿顶两端节点处于啮合状态时，其变形量λ2≈8μm和λ2≈7.2μm；其中以上位置获得的节点变形量与各最大修形量ep、e1、e2之间的关系为

(11)则可得：齿廓修形最大修形量ep取6μm，齿向最大修形量e1、e2分别取2μm和1.2μm。通过上述确定的各最大修形量和修形曲线，即可获得各节点因综合修形而减少的渐开线发生线长度，并将其转化为该节点的坐标变化，直接生成新的节点坐标和齿的轮廓线。由于之前已提取了各齿面节点坐标和编号，故在修形中，仅仅通过APDL就可以实现对节点的坐标变化，无需对几何结构进行重新设计和有限元网格重新划分，从而大大提高了修形分析的计算效率。由于修形量都很小，故节点移动所造成的单元形状改变也很小，不会因有限元网格的畸变而导致计算不收敛或计算精度下降。

对于本算例，修改各齿面的节点三维坐标，获得修形后的齿轮模型。按原有外载和各约束参数重新进行瞬态分析，结果如图9、图10所示。

图9 修形之后齿顶啮合时轴向应力分布Fig.9 Stress on path of axial direction of gears after profile modification

图10 修形之后齿顶啮合时沿渐开线方向的应力分布Fig.10 Stress on path of evolvent direction of gears after profile modification

提取原来齿面在原来路径上处于啮合状态时的应力和应变变化情况，发现沿渐开线方向上的最大应力从443.8MPa减小到363MPa，最小应力从63.455MPa减小到51.031MPa，沿渐开线方向齿根处应力也从158.55MPa减小到129.024MPa，相应节点的应力值对比结果见表3。通过对比，发现齿轮综合修形能明显地改善齿面上的应力大小和应力分布情况。

表3 齿顶啮合时修形后轴向啮合应力分布比较

疲劳应力S与寿命N的函数关系为

N=CS-M

(12)

考虑到风力发电机通常处于沙漠、海岛等，维修成本极高，对可靠性要求也较高。不同损伤概率下的齿轮材料20GrMnTi钢的参数C、M值见表4，在表4中取损伤概率为1%进行寿命计算。由于齿轮接触应力是脉冲式的应力分布，按照脉动循环应力计算得到齿轮疲劳寿命，并进行Goodman修正以消除非0均值应力的影响。其中两个关键点(齿顶边缘、齿顶中间)修形前后疲劳寿命的变化情况见表5。可以看到，经过齿轮的综合修形，齿面最大应力减小约20%，其危险节点寿命增加了2个数量级，普通齿顶节点寿命增加了4倍左右，从而延长了齿轮寿命，提高了可靠性，达到了齿轮修形的目的。

表4 20GrMnTi钢P-S-N曲线参数估计

表5 两节点修形前后疲劳寿命比较

4 结论

(1)本文提出了根据有限元齿轮瞬态动力学结果计算齿面关键点的弹性变形的计算方法，可以方便地得到齿轮的最大修形量。

(2)本文建立了基于有限元应力结果的齿轮修形程序，通过有限元模型的建立、动力学分析、关键节点的变形提取、修形计算、模型修改等步骤实现了程序化自动计算，从而降低了修形分析的难度，有利于本方法在企业的推广。

(3)本文以某风电三级齿轮为研究对象，进行了综合修形，通过结果对比分析可知，齿面综合修形能显著降低齿面啮合过程中的应力，增加齿轮疲劳寿命，提高齿轮可靠性。提出了一种瞬态分析齿轮变形量的计算方法。

[1] CHEN P C, HUANG A C.Adaptive Multiple-surface Sliding Control of Hydraulic Active Suspension Systems Based on the Function Approximation Technique[J].Journal of Vibration and Control,2005,11(5):685-706.

[2] MARZBANRAD J.Optimal Active Control of Vehicle Suspension System Including Time Delay and Preview for Rough Roads[J].Journal of Vibration and Control,2002,8(7)：967-991.

[3] 汪强,周锦进.齿轮修形技术及工艺的发展[J].电加工与模具,2001(5): 21-24. WANG Qiang, ZHOU Jinjing. The Development of the Gear Tooth Modification and Technology[J]. Electromachining & Mould,2001,(5):21-24.

[4] 袁哲,孙志礼,郭瑜.直齿圆柱齿轮齿廓修形曲线优化设计[J].机械传动,2010,34(5): 1-4. YUAN Zhe, SUN Zhili, GUO Yu. Optimal Design of Profile Modification Curves for Spur Gears[J]. Journal of Mechanical Transmission,2013,34(5):1-4.

[5] 陶燕光.高速齿轮热变形修形的试验研究[J].齿轮,1988,12(2):25-28. TAO Yanguang. The Test Research of Heat-deformation Profile Modification about High-speed Gears[J].Gears, 1988,12(2):25-28.

[6] 杨廷力,叶新.渐开线高速齿轮的齿高修形[J].齿轮，1982,6(3):14-24. YANG Tingli, YE Xin. The Tooth High Profile Modification about Evolvent High-speed Gears[J]. Gears,1982,6(3):14-24.

[7] 周志峰.渐开线圆柱直齿轮修形的分析与研究[D].北京:北京交通大学,2014:32-50. ZHOU Zhifeng. Analysis and Research on Modification Involute Spur Gear[D].Beijing:Beijing Jiaotong University,2014:32-50.

[8] 刘旦.2.5 MW风电齿轮箱齿轮修形研究[D].银川：宁夏大学,2013:10-32. LIU Dan.Research on Tooth Modification of 2.5 MW Wind Power Gearbox[D].Yinchuan：Ningxia University, 2013:10-32.

[9] 王亮.渐开线圆柱齿轮齿廓修形设计及动态接触分析[D]. 太原:太原理工大学,2012:18-51. WANG Liang. Modification and Design of Tooth Profile and Analysis of Dynamic Contact for a Cylinder Gear of Involute[D]. Taiyuan:Taiyuan University of Technology, 2012:18-51.

[10] WANG Jiande, Howard L. Finite Element Analysis of High Contact Ratio Spur Gears in Mesh[J]. Transactions of the ASME,2005,127(3):11-15.

[11] 颜海燕.直齿轮轮齿变形计算的数值积分法[J].机械传动，2005,9(2):7-9. YAN Haiyan. Numerical Integral Method for Gear Deformation Calculation[J]. Journal of Mechanical Transmission，2005,9(2):7-9.

[12] 尚振国.风力发电机增速器齿轮修形技术研究[D].大连：大连理工大学,2010. SHANG Zhenguo. Research on Gear Modification of High-power Speed Increasing Gearboxes for Wind Turbines[D]. Dalian:Dalian University of Technology, 2010.

[13] 郝东升.齿轮啮合数值分析建模方法及其应用研究[D].大连：大连理工大学,2012:85-121. HAO Dongsheng. Research on Numerical Analysis Modeling Method for Gear Meshing and Its Applications[D]. Dalian:Dalian University of Technology,2012:85-121.

[14] 李润方.齿轮传动的刚度分析和修形方法[M].重庆:重庆大学出版社,1998. LI Runfang. The Stiffness Analysis and Profile Modification Method about Gear Transmission[M]. Chongqing: Chongqing University Press,1998.

[15] 孙恒.机械原理[M].5版.北京:高等教育出版社,2006:256. SUN Heng. Theory of Machines and Mechanisms[M].5th ed.Beijing:Higer Education Press,2006:256.

[16] 王文奎.机械原理[M].北京:电子工业出版社, 2007. WANG Wenkui. Theory of Machines and Mechanisms[M]. Beijing:Publishing House of Electronics Industry,2007.

[17] SIGG N. Tooth Profile Modification of High Speed Duty Gear[C]//Proceedings of International Conference on Gearing. New York, 1958:313-316.

[18] 博弈创作室.APDL参数化有限元分析技术及其应用实例 [M].北京:中国水利水电出版社,2004. BOYI Booker. FEM Analysis and Application Example Using APDL[M]. Beijing:China Water & Power Press,2004.

[19] 田雪.大型风电机组齿轮箱齿轮的疲劳寿命分析与预测 [D].北京:华北电力大学,2014:23-30. TIAN Xue. Analysis and Prediction of Gears Fatigue in Large Wind Turbine Gearbox[D]. Beijing:North China Electric Power University, 2014:23-30.

[20] 周志刚. 随机风作用下风力发电机齿轮传动系统动力学及动态可靠性研究[D].重庆:重庆大学,2012:85-95. ZHOU Zhigang. Study on Dynamics and Time-dependent Reliability of Gear Transmission System of Wind Turbine under Random Wind Conditions[D]. Chongqing:Chongqing University,2012:85-95.

(编辑 陈 勇)

Calculation Method of Tooth Profile Modification Based on Mesh Node Deformations

LI Jianmin1WU Yuecheng1XUAN Haifeng1CHEN Wenhua1PAN Xiaodong2WANG Wei1ZHOU Xinteng1

1.Key Laboratory of Reliability Technology for Mechanical and Electronic Product, Zhejiang Sci-Tech University,Hangzhou ,310018 2.Hangzhou Advance Gearbox Group Co.,Ltd.,Hangzhou, 310027

In order to deal with the problems of gear modification, a method to calculate the displacements of the gear tooth related nodes was proposed, which was based on the finite element dynamics calculation. The elastic deformations of the teeth were obtained by extracting the displacements of the gear teeth, which were used as the gear’s maximum profile.The calculation program was used to calculate the gear tooth profiles. Taking a wind power gear box for example, the maximum stress of the gear is reduced by about 20%, and the life of the gear is prolonged, and the reliability of the gear box is increased.

gear modification; finite element analysis; node deformation; stress

2015-10-08

国家国际科技合作专项(2015DFA71400)；浙江省自然科学基金资助项目(LY13E050025)

TH132.413;TP319

10.3969/j.issn.1004-132X.2017.03.007

李剑敏，男，1962年生。浙江理工大学机械与自动控制学院教授。主要研究方向为机械结构强度与可靠性分析。获省级科技进步三等奖1项。发表论文30余篇。E-mail:ljmzrz@163.com。吴跃成，男，1968年生。浙江理工大学机械与自动控制学院副教授。宣海枫，男，1989年生。浙江理工大学机械与自动控制学院硕士研究生。陈文华，男，1963年生。浙江理工大学机械与自动控制学院教授、博士研究生导师。潘晓东，男，1969年生。杭州前进齿轮箱集团股份有限公司教授级高级工程师。王 威，男，1990年生。浙江理工大学机械与自动控制学院博士研究生。周新腾，男，1992年生。浙江理工大学机械与自动控制学院硕士研究生。