弧齿锥齿轮有限元建模与接触分析

侯祥颖，方宗德，邓效忠，蔡香伟，宁程丰

（1.西北工业大学机电学院，陕西西安710072；2.河南科技大学机电工程学院，河南洛阳471023）

弧齿锥齿轮有限元建模与接触分析

侯祥颖1，方宗德1，邓效忠2，蔡香伟1，宁程丰1

（1.西北工业大学机电学院，陕西西安710072；2.河南科技大学机电工程学院，河南洛阳471023）

为了使有限元软件能够对齿面微米量级变化的齿轮进行分析，提出了基于高精度建模的有限元分析方法。将齿面的理论计算公式转化成数值计算算法，按照有限元理论和有限元软件的编程语法将齿轮实体进行节点坐标求解和网格单元划分，可以求得一对装配好的有限元齿轮模型。高精度建模的高重合度弧齿锥齿轮算例中，将有限元计算结果与实验结果对比，发现两者的一致性：弯曲应力曲线形状一致，变化规律相近；接触印痕区域一致，变化趋势相同。高精度建模有限元分析能够清晰的反映齿面微米量级改变导致的印痕和应力的变化。

弧齿锥齿轮；有限元；接触分析；弯曲应力；接触印痕

弧齿锥齿轮是一种不完全共轭局部点接触的齿轮副，齿面形状和加工机床的调整都比较复杂。目前国外已经较多地将弧齿锥齿轮应用于船舰传动系统，国内在船舰上的应用相对较少，这方面的研究正在进行。弧齿锥齿轮传动平稳，是动力传输的关键部件，在强度、动态性能与可靠性方面等方面都有很高的要求。邓效忠等［1］曾提出高重合度弧齿锥齿轮的设计方法，采用LTCA理论，通过有限元应力影响矩阵法对其相关性能进行了详细讨论［2］。高重合度弧齿锥齿轮是通过齿面修形改善其啮合性能的典型，不同的弧齿锥齿轮齿面精度会使弧齿锥齿轮有不同的传动性能。本文以一对高重合度弧齿锥齿轮为例，通过高精度齿面建模和装配，导入有限元分析软件分析，并与实验结果进行了对比验证分析。

1 高重合度弧齿锥齿轮

文献［3⁃5］等利用NURBS曲面函数对齿面进行重构来建立弧齿锥齿轮模型。但是此方法计算求得模型的齿面与理论齿面仍有一定误差；在对一对齿轮进行有限元接触分析时，采用手动装配，可能存在干涉和装配误差等现象，产生计算误差或错误。对于齿面修形处于微米量级的齿轮来说，上述方法建模精度远远达不到要求。采用曲面拟合方法的误差分布如图1所示，误差可能达到十几甚至几十微米，且有很大随机性，分布不均匀；而且此方法需要进行三维建模，然后再将模型导入有限元软件进行实体分割、网格划分等操作，工作量大，可重复性差。

在保证齿轮基本参数（如表1所示）不变的条件下，采用局部综合法使接触路径倾斜和改变齿面的失配量，设计了3种重合度不同的弧齿锥齿轮。（详见文献［2］）高重合度设计方法不改变齿轮轮坯基本参数，只对小轮齿面加工参数进行调整，而保持大轮的加工参数不变，通过小轮齿面的改变使接触路径倾斜，接触路径延长，以达到增大重合度的目的。小轮齿面的修形量往往只有十几微米甚至几微米，故要求有限元模型要有极高的精度才能反映出齿面变化的影响。

图1 NURBS曲面重构误差分析Fig.1 Error analysis of NURBS surface reconstruction

表1 弧齿锥齿轮参数Table 1 Basic parameters of spiral bevel gear

本文算例采取的倾斜角分别为20°、50°、80°接触路径，通过TCA分析得到传动误差如图2所示，图2（a）理论设计重合度只有1.25左右，图2（b）理论设计重合度增加至1.75左右，图2（c）理论设计重合度增加至2.0以上，传动误差曲线明显可以看出存在三齿啮合区。

图2 3种不同齿轮的传动误差Fig.2 Transmission error of three different designs

2 弧齿锥齿轮有限元模型建立

2.1 通过数值计算求解有限元模型

弧齿锥齿轮齿面较为为复杂，其加工原理和方法繁琐，加工参数较多，一般很难得到齿面方程的解析表达式，多采用数值计算的方法，求得齿面上离散的点，即可用离散点构建有限元模型［4⁃5］。

将刀具切削面方程从刀具坐标系变换到轮坯坐标系下，联立齿面的位矢方程和啮合方程如下（文献［6］）：

式中：M为刀具坐标系到轮坯坐标系的变换矩阵，r为刀盘方程，μ、θ分别为切削刃产生圆锥曲面坐标，φ为转角，n、v分别为齿面点法矢和相对速度矢。由式（1）通过编程进行数值计算，可得出齿面上的离散点在轮坯坐标系下的坐标。

按照有限元软件ABAQUS对单元网格以及节点的书写语法要求，首先计算理论齿面上的点，将齿面点坐标绕轴线旋转，对两齿面之间采用点阵进行填充，然后对所有节点按照一定顺序进行自然数编号；之后按顺序编写8节点单元，单元由8个节点自然序号构成，故点阵的密度决定了网格单元的密度。

装配时，选取大轮接触迹线与齿顶线的交点为啮合参考点，即恰好进入啮合的位置。利用轮齿接触分析求得啮合点的坐标以及小、大轮对应的转角φ1、φ2（式（2）），采用坐标转换矩阵，对大小轮分别旋转φ1、φ2，即能通过理论计算完成一对弧齿锥齿轮的装配。

式中：下标1、2分别表示大小齿轮，下标p、g分别表示大小轮的齿面参数，上标h表示在固定坐标系Sh下，分别表示齿面离散点位矢在旋转投影面上的X、Y坐标，ha2表示大轮齿顶高，AFO表示外锥距，γ2为大轮根锥角。

通过编写inp文件将装配坐标系下大小轮的节点和单元数据导入ABAQUS软件中，可直接得到有限元装配模型（如图3所示）。

图3 有限元装配模型Fig.3 Assembly model of finite element

2.2 有限元模型前处理

进行有限元齿轮接触分析，单元一般采用六面体一次减缩积分单元（C3D8R），单元属性在编写inp文件生成有限元模型时定义（文献［7］）；将已建好的弧齿锥齿轮的五齿有限元模型导入ABAQUS软件中，需要进行以下设置（文献［5，8］）：

1）定义材料：将大小轮赋予相同的材料属性，其弹性模量为2.06×105MPa，泊松比为0.3，密度为7.8×10-9t／mm3。

2）定义接触对：分别将大小轮的齿面、过渡曲面、齿顶等形成的连续面选中，并将此两对连续面设置为接触对。在接触属性中，采用无摩擦硬接触，计算方法采用Kinematic接触方法。

3）设置参考点：取轴线上任一点为参考点，并在参考点和小轮内圈和剖面间建立耦合刚体约束。

4）分析类型的定义：采用显示动力学分析算法，设置分析时间并根据网格数量选择适当的质量放大系数。显式算法基于动力学方程，无需迭代，有较好的稳定性。

5）场变量和历史变量设置：设置接触压强、接触力、应力、位移等为场变量输出量，动能、内能为历史变量输出量。

6）边界条件和载荷：大小轮分别只释放绕自身轴线旋转的自由度，施加小轮的绕轴线的转速为0.2 rad／s，大轮施加阻力矩为7.5×105N·mm，小轮转速和大轮转矩均采用平缓加载方式；从而模拟出以小轮驱动大轮，同时小轮受到大轮的阻力矩的作用，在齿面上产生压力的工况。

同样将以上设置步骤编写入inp文件中，这样inp文件就包括了所有模型信息和条件设置信息。直接将inp文件导入ABAQUS中，提交计算即可，无需再人工操作。

3 结果分析

3.1 有限元仿真结果

前处理结束后，用ABAQUS软件分析可得到弧齿锥齿轮的有限元分析结果。图4为某一时刻的小轮凹面和大轮凸面（工作面）的Mises应力云图，图中可以看到3对齿轮同时啮合的状态。

图4 有限元计算Misse应力云图Fig.4 Mises stress nephogram of finite element calculation

改变运动时间，当给小轮施加逆时针转速时，小轮的凹面和大轮凸面发生接触，小轮驱动大轮顺时针旋转，同时大轮的阻力矩模拟负载，使两齿啮合点处产生应力和变形。小轮凹面接触椭圆从大端齿根逐渐向小端齿顶移动，大轮凸面接触区域接触椭圆从大端齿顶逐渐向小端齿根移动；最大主应力位于接触椭圆区域的中心，接触椭圆在齿面中间面积最大，表示在齿面中间轮齿承受载荷最大。

3.2 高重合度对齿面接触力的影响

有限元计算结束生成结果文件，对运行编译好的脚本文件，从计算结果模型中提取所需要的节点坐标，接触压强和Mises应力等，并通过数据处理软件将应力曲线绘出，如图5所示。

图5 3种接触路径的接触应力对比Fig.5 Contrast of contact stress with three different designs

通过对比容易发现，在计算载荷一定的情况下（大轮扭矩7.5×105N·mm），接触路径倾斜度越大，即倾斜角越小，则重合度越大，接触应力越小。倾斜角80°设计中接触应力最大值为1 111.84 MPa，而倾斜角20°设计中的最大值为1 164.98 MPa；接触路径倾斜前后相比，接触应力下降4.56%，并且，随着接触路径的倾斜，重合度的增大，齿面接触应力变化趋于平缓，曲线上升下降都明显变化缓慢，即加载明显趋于平稳。

4 弯曲应力实验对比

本实验由调频电机（30 kW）提供实验驱动力，经减速器减速，利用转矩转速传感器将动力输给试验齿轮小轮；试验齿轮由磁粉制动器施加负荷，装配在试验箱内。通过测量齿根应变获得齿根弯曲应力，通过转矩转速测试仪测得试验齿轮的转矩和转速。电阻应变计贴在大轮凸面的齿根中点处，位于非啮合区，在远离啮合齿面处的齿坯上贴2个温度补偿片，分别与凸面和凹面的应变计组成半桥连接。引线由集流环引出，经应变仪放大，通过多功能信号测试仪纪录［9］。通过实验测得齿根应力测试曲线如图6所示。图7为对应齿轮的有限元仿真计算所得到的大轮齿根中点弯曲应力过程，详细数值比较如表2。通过图6与图7以及表2对比，大轮齿根中点应力变化过程中，实验测得波形倾斜角20°设计最大电压信号达1.053 V，而倾斜角80°设计最大电压信号0.597 V，二者相比峰值下降43.3%；有限元分析结果中，倾斜角20°设计最大弯曲应力为126.08 MPa，而倾斜角80°设计最大弯曲应力值仅为83.28 MPa；接触路径倾斜前后相比，弯曲应力下降34%，且当接触路径越倾斜时，即重合度增大，弯曲应力数值明显下降。随着接触路径倾斜角加大，接触路径延伸增长，齿根弯曲应力逐渐降低，变化趋于平缓，仿真结果与实验结果曲线变化规律一致，形状近似。实验数据与有限元分析结果有一定的误差，误差的主要来源有：实验仪器条件，实验系统误差，应变计所贴位置，温度补偿片的选用；有限元所选用的单元类型、网格密度、计算方法、质量放大系数等。

图6 三对齿轮齿根应力测定曲线Fig.6 Bending stress curves of experiment results

图7 有限元大轮齿根中点弯曲应力过程Fig.7 Bending stress changes of midpoint on gear root

表2 大轮齿根应力实验数据与有限元分析比较Table 2 Bending stress comparison between experi⁃ment results and finite element analysis

5 接触印痕实验对比

有限元分析结束后，编写读取abaqus结果文件（.odb文件）的脚本文件，求得某齿一个啮合周期内每一瞬时的接触区域内的点坐标，利用数据处理软件编程，将所有点叠加到同一坐标系下，并将其投影到旋转投影面，即得到有限元分析的接触印痕区域。有限元计算和实验得到接触印痕区域分别如图8和图9所示。

图8 有限元印痕提取Fig.8 Contact pattern of finite element analysis

图9 实验印痕Fig.9 Contact pattern of experiment

通过对比图8和图9，随着倾斜度逐渐变大，接触区域倾斜，接触路径增长，接触椭圆长轴变短，有限元计算结果与实验结果印痕区域一致，变化规律相同，由此证明此高精度建模方法及有限元分析方法的正确性。

6 结论

1）本文有限元建模采用数值计算的方法，得到节点坐标，并通过编号导入inp文件中直接建立有限元模型，省去了三维建模以及划分网格的过程，通过与NURBS曲面重构结果对比，证明此方法有限元齿面精度及装配精度极高，可重复性好，适用于高精度的齿轮接触分析。

2）通过有限元计算与实验结果对比，随着接触路径倾斜度逐渐变大，有限元计算印痕接触路径逐渐增长，接触椭圆长轴逐渐变短，印痕区域和变化规律与实验一致；有限元计算齿根弯曲应力逐渐降低，变化趋于平缓，与实验弯曲应力变化趋势一致，形状近似，验证了本文方法的计算精度。

［1］邓效忠，方宗德，张金良等.弧齿锥齿轮的高重合度设计［J］.中国机械工程，2002，14（19）：791⁃795.DENG Xiaozhong，FANG Zongde，ZHANG Jinliang，et al.Design of spiral bevel gears with high contact ratio［J］.Chi⁃na Mechanical Engineering，2002，14（19）：791⁃795.

［2］邓效忠，方宗德，魏冰阳.高重合度弧齿锥齿轮的接触应力与弯曲应力分析［C］／／中国工程机械学会2003年年会论文集，2003：107⁃111.DENG Xiaozhong，FANG Zongde，WEI Bingyang.Analysis of contact stress and bending stress on spiral bevel gears of high contact ratio［C］／／China Mechanical Engineering Soci⁃ety 2003 Annual Meeting Proceedings，2003：107⁃111.

［3］唐进元，蒲太平.弧齿锥齿轮动态啮合有限元数值分析［J］.机械科学与技术，2009，8（8）：981⁃985.TANG Jinyuan，PU Taiping.FEM numerical analysis of the dynamic engagement of spiral bevel gears［J］.Mechanical Science and Technology for Aerospace Engineering，2009，8（8）：981⁃985.

［4］ASTOUL J，GENEIX J，MERMOZ E，et al.A simple and robust method for spiral bevel gear generation and tooth con⁃tact analysis［J］.International Journal on Interactive Design and Manufacturing，2013，7（1）：37⁃49.

［5］唐进元，蒲太平.基于有限元法的螺旋锥齿轮啮合刚度计算［J］.机械工程学报，2011，47（11）：23⁃29.TANG Jinyuan，PU Taiping.Spiral bevel gear meshing stiff⁃ness calculations based on the finiteelement method［J］.Journal of Mechanical Engineering，2011，47（11）：23⁃29.

［6］LITVIN F L.Gear geometry and applied theory［M］.Cam⁃bridge University Press，2004：613⁃635.

［7］GONZALEZ⁃PEREZ I，ISERTE J L，FUENTES A.Imple⁃mentation 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.

［8］SHEVELEVA G I，VOLKOV A E，MEDYEDEV V I.Algo⁃rithms for analysis of meshing and contact of spiral bevel gear［J］.Mechanism and Machine Theory，2007，42（2）：198⁃215.

［9］邓效忠.高重合度弧齿锥齿轮的设计理论及实验研究［D］.西安：西北工业大学，2002：49⁃77.DENG Xiaozhong.Research of design theory and experi⁃ments on spiral bevel gears with high contact ratio［D］.Xi'an：Northwestern Polytechnical University，2002：49⁃77.

［10］WANG P Y，FAN S C，HUANG Z G.Spiral bevel gear dy⁃namic contact and tooth impact analysis［J］.Journal of Me⁃chanical Design，2011，133：084501.

［11］GONZALEZ⁃PEREZ I，RODA⁃CASANOVA V，FUENTES A，et al.A finite element model for consideration of the tor⁃sional effect on the bearing contact of gear drives［J］.Jour⁃nal of Mechanical Design，2012，134：071007.

Contact analysis of spiral bevel gears based on finite element model

HOU Xiangying1，FANG Zongde1，DENG Xiaozhong2，CAI Xiangwei1，NING Chengfeng1

（1.Department of Mechanical Engineering，Northwestern Polytechnical University，Xi'an 710072，China；2.Department of Mechani⁃cal Engineering，Henan University of Science and Technology，Luoyang 471023，China）

To make analysis on gears whose tooth surfaces have micron level changes，the finite element method based on high⁃precision finite element model is proposed.The theoretical formulas of tooth surfaces are converted into numerical algorithm，the coordinates of nodes and finite mesh elements of the gears are obtained based on the finite element theory and program syntax of finite element software，and the high⁃precision finite element model of a pair of gears is built.Through analysis of spiral bevel gears with high contact ratio，the simulation results are in good agreement with the experimental ones.The bending stress curves show the same shape and accordant changing trend，but the contact patterns show the same location and changing trend.The finite element analysis with high⁃precision modeling could reflect the changes of stresses and contact pattern caused by modification to the surfaces of gears in micron level.

spiral bevel gear；finite element method；contact analysis；bending stress；contact pattern

10.3969／j.issn.1006⁃7043.201312055

TH132.41

：A

：1006⁃7043（2015）06⁃0826⁃06

http：／／www.cnki.net／kcms／detail／23.1390.u.20150428.0910.012.html

2013⁃12⁃27.网络出版时间：2015⁃04⁃28.

国家自然科学基金资助项目（51175423，51375384，51205310）.

侯祥颖（1990⁃），男，博士生；方宗德（1948⁃），男，教授，博士生导师.

侯祥颖，E⁃mail：houxiangying＠126.com.