技术领域
[0001] 本
发明属于机械零部件疲劳失效风险的预测方法，具体涉及一种预测以
渗碳淬火为代表的表面硬化齿轮副弹塑性接触疲劳点蚀与剥落失效风险的方法。
背景技术
[0002] 在高速、重载机械中，渗碳硬化齿轮的接触疲劳失效问题已成为限制装备可靠性的重要
瓶颈之一，硬化齿轮的接触疲劳失效形式常见的有点蚀和剥落两种，这两种失效问题会带来极大地安全隐患和经济损失，而硬化齿轮的接触疲劳失效影响因素众多，如工况、材料
力学特征的梯度、残余
应力等，此
外齿轮接触过程中伴随时变多轴应力状态，使得目前对于渗碳重载齿轮的接触疲劳机理的认识与控制不够深入，工程实际中对于失效风险的预测存在很大困难。
发明内容
[0003] 本发明所要解决的技术问题是提供一种预测表面硬化齿轮弹塑性接触疲劳点蚀与剥落失效风险的方法，它能预测齿轮在弹塑性接触状态下疲劳点蚀和疲劳剥落的失效风险，预测结果在工程实际中具有指导作用，避免因齿轮失效而导致的突发事故、因装置设备停机造成生产效益的损失。
[0004] 本发明所要解决的技术问题是通过这样的技术方案实现的，它包括以下步骤：
[0005] 步骤1、根据齿轮材料的硬度曲线和硬度与
屈服强度的线性关系，计算得到齿轮材料局部屈服强度；
[0006] 步骤2、依据齿轮副的几何参数，将
啮合位置的接触状态简化为二维接触模型，同时基于ABAQUS平台建立该模型；
[0007] 步骤3、基于Python编程语言和齿轮材料局部屈服强度，为二维接触模型添加材料属性；
[0008] 步骤4、根据齿轮的材料和Dang Van多轴应力准则，结合齿轮的具体工况，计算齿
轮齿面下的Dang Van等效应力；
[0009] 步骤5、以Dang Van等效应力除以齿轮材料局部屈服强度，得出关键接触区域任一点处的疲劳失效风险值；根据失效风险值在不同区域相对大小和最大值所在的深度位置预测点蚀或剥落失效的风险。
[0010] 本发明的技术效果是：
[0011] 解决了机械行业长期存在的技术难题，能在考虑齿轮硬化后材料力学性能梯度变化的条件下，分析齿轮在弹塑性状态下的点蚀和剥落的风险，为齿轮的设计、制造、使用提供依据。
附图说明
[0012] 本发明的附图说明如下：
[0013] 图1为渗碳硬化齿轮失效风险区分布图；
[0014] 图2为表面硬化齿轮硬度随浅层深度变化的曲线图；
[0015] 图3为
实施例中齿轮材料局部屈服强度沿深度分布曲线；
[0016] 图4为齿轮啮合接触状态的简化示意图；
[0017] 图5为实施例中有限元模型示意图；
[0018] 图6为实施例中某兆瓦级风电齿轮箱的传动图；
[0019] 图7为实施例中Dang Van等效应力沿深度分布的曲线；
[0020] 图8为实施例的齿轮的接触疲劳失效风险值曲线。
具体实施方式
[0021] 下面结合附图和实施例对本发明作进一步说明：
[0022] 如图1所示，根据工程经验，在距离
齿面几百微米处为疲劳点蚀的高发区，而距离齿面距离为1～3mm为疲劳剥落的高发区。图示中的曲线代表疲劳失效风险值沿深度方向的变化曲线，根据不同区域失效风险的大小以及最大失效风险值的深度位置可以用于预测齿轮发生疲劳点蚀和疲劳剥落失效的风险，若最大值发生在点蚀高发区，则齿轮在给定工况下的点蚀风险较高，若最大失效风险值在剥落高发区发生，则齿轮剥落失效风险较高。
[0023] 如图2所示，曲线为使用维氏硬度测试法得到的硬度曲线。纵坐标为维氏硬度值，横坐标为距离齿面的深度值，单位为毫米(mm)，总测量深度为5mm，测量值为一系列数据点，根据测量数据点拟合出一条曲线，拟合函数为：
[0024] HV(y)＝a·y4+b·y3+c·y2+d·y+e (1)
[0025] 式(1)中，a＝-1.733,b＝25.37,c＝-118,d＝144.8,e＝589.5，y为深度值。
[0026] 图2所示的硬度值对于不同的齿轮需要重新测量，因为该值与齿轮的材料和
热处理工艺等有关，而测量硬度值的方法众多，这里只提供一个参考的实施例，其他方法所获得的硬度值也可以用于本
专利提供的方法进行失效风险预测，这里不再赘述。
[0027] 本发明包括以下步骤：
[0028] 步骤1、根据齿轮材料的硬度曲线和硬度与屈服强度的线性关系，计算得到齿轮材料局部屈服强度；
[0029] 根据“Correlation of Yield Strength and Tensile Strength with Hardness for Steels”,E.J.Pavlina and C.J.V.Tyne,Journal of Materials Engineering&Performance,vol.17,pp.888-893,2008.(“
钢的
抗拉强度与屈服强度之间的关系”，E.J.Pavlina and C.J.V.Tyne，材料工程与性能，第888-893页，2008年)记载了维氏硬度与屈服极限的线性关系为：
[0030] σYS(y)＝-90.7+2.876HV(y) (2)
[0031] 式中，y代表齿轮深度，σYS为齿轮材料局部屈服强度，HV为沿深度y分布的维氏硬度值。
[0032] 由公式(2)计算得到的齿轮材料屈服强度沿深度的变化如图3所示。
[0033] 步骤2、依据齿轮副的几何参数，将啮合位置的接触状态简化为二维接触模型，同时基于ABAQUS平台建立该模型；
[0034] 齿轮副啮合接触状态简化为二维接触模型的过程如图4所示，左图代表要分析的啮合齿轮副，根据两齿轮啮合点处的
曲率半径将分析对象转化为中间图所示的两圆接触，该模型又可进一步简化为一刚性半圆与柔性体的接触(如右图)。此时的刚性半圆的综合
曲率半径和柔性体的等效弹性模型的计算方法为：
[0035] R＝R1R2/(R1+R2) (3)
[0036]
[0037] 上式中，R1，R2为两齿轮接触位置的曲率半径，R为综合曲率半径，E1，E2为两齿轮的
弹性模量，E为等效弹性模量，υ1，υ2为两齿轮的泊松比。“AGMA information sheet 908-B89,1989“.Geometry factors for determining the pitting resistance and bending strength of spur,helical and herringbone gear teeth”.”(美国齿轮标准1989年的“判定
直齿轮，
斜齿轮，人字形齿轮轮齿的抗点蚀能力和弯曲强度的几何影响因子”)在第5-7页记载了R1，R2的计算方法。
[0038] 步骤3、基于Python编程语言和步骤1得到的齿轮材料局部屈服强度，为二维接触模型添加材料属性；
[0039] 基于Python语言对ABAQUS有限元
软件进行二次开发，编写对柔性体进行分层的程序。
[0040] 如图5所示，左图表示使用Python语言对柔性体上部分进行分层后的局部放大图，
单层厚度为0.01mm，中间图代表在ABAQUS平台中建立的二维接触模型，右图为该模型添加齿轮材料局部屈服强度后的示意图，右图中的曲线代表从齿面到芯部的屈服强度分布。
[0041] 在该步骤中，可以根据实际工程需要，对采用不同的单层厚度，均属于本专利的范畴，在此不再赘述。
[0042] 步骤4、依据步骤3中建立的模型，计算齿面下的Dang Van等效应力；
[0043] "On a new multiaxial fatigue limit criterion:theory and application,"K.Dang Van,B.Griveau,and O.Message,Biaxial and Multiaxial Fatigue,M.W.Brown and K.J.Miller,Eds.,ed:Mechanical Engineering Publications,London,1989,pp.479-496.(“一种新的多轴疲劳极限准则：理论与应用”，K.Dang Van,B.Griveau,and O.Message，两轴与多轴疲劳，M.W.Brown and K.J.Miller，编辑，机械工程出版，伦敦，1989年，479-496页)中提出了Dang Van等效应力的计算公式：
[0044] τmax(θ,t)+ασh(t)＝σequi,DangVan(y) (5)
[0045]
[0046]
[0047] 上述式子中，σequi,DangVan(y)是Dang Van多轴疲劳准则的等效应力，τmax(θ,t)为最大剪应力幅值，σh(t)为
水静应力，α是材料参数；τ-1是扭转疲劳极限，σ-1是弯曲疲劳极限，σ1、σ2和σ3是第一主应力、第二主应力和第三主应力。
[0048] "Prediction of contact fatigue for the rough surface elastohydrodynamic lubrication line contact problem under rolling and sliding conditions",Q.Hua,PhD Thesis,Cardiff University,UK,2005.(“滑动和滚动条件下的粗糙表面弹流润滑线接触问题的接触疲劳预测”，Q.Hua，博士论文，卡迪夫大学，英国，2005年)中指出最大剪应力幅值的求法为：
[0049] τ(θ,t)＝τxy(cos2θ-sin2θ)+(σy-σx)sinθcosθ (8)
[0050] τmax(θ,t)＝max{τa(θ,t)＝|τ(θ,t)-τm|} (9)
[0051] 式(8)、式(9)中，τ(θ,t)是材料点瞬时剪应力，τxy是x-y平面剪应力，σx是x轴向正应力，σy是y轴向正应力，θ代表剪应力幅值达到最大值的平面与滚动方向之间的夹
角；τmax(θ,t)是最大剪应力幅值，τa(θ,t)是剪应力幅值，τm是剪应力均值。
[0052] 步骤5、以Dang Van等效应力除以齿轮局部材料屈服强度，得出关键接触区域任一点处的疲劳失效风险值。疲劳失效风险值越大，疲劳失效的可能性越大，根据风险值在不同区域的大小，预测点蚀或剥落的失效风险。疲劳失效风险值的计算式如下：
[0053]
[0054] 式中Aff(y)为疲劳失效风险值。
[0055] 实施例
[0056] 如6图所示，样品齿轮来自于某兆瓦级
风力发
电机齿轮箱中间及齿轮副，它在工程实际应用时，失效的概率明显高于其他齿轮。
[0057] 齿轮副的主要参数如下：
[0058]
[0059] 图2为该样本所对应的拟合硬度曲线，根据实施例齿轮参数，以齿轮啮合线上
节点处啮合瞬时状态为例，逐步计算出齿轮副接触的疲劳时效风险值：
[0060] 步骤1、根据齿轮硬度曲线和公式(2)，可得齿轮材料局部屈服强度，如图3所示。
[0061] 步骤2、采用图4所示的简化过程将啮合齿轮简化为二维接触模型。根据公式(3)～(4)，由R1＝236.28mm，R2＝44.18mm，计算得到齿轮副综合曲率半径R＝37.22mm；等效弹性模量E＝1.15×1011Pa；同时在ABAQUS平台建立该二维接触模型。
[0062] 步骤3、使用Python语言和步骤1得到的齿轮材料局部屈服强度，为二维接触模型添加材料参数。Python语言程序示例如下:
[0063] for i in range(1,500):
[0064] s.Line(point1＝(-20.0,10-0.01*i),point2＝(20.0,10-0.01*i))
[0065] p＝mdb.models['Model-1'].parts['deformable']
[0066] f＝p.faces
[0067] pickedFaces＝f.findAt(((6.666667,-6.666667,0.0),))
[0068] e1,d2＝p.edges,p.datums
[0069] p.PartitionFaceBySketch(faces＝pickedFaces,sketch＝s)
[0070] s.unsetPrimaryObject()
[0071] del mdb.models['Model-1'].sketches['__profile__']
[0072] 步骤4、在ABAQUS平台使用公式(5)～(9)计算Dang Van等效应力，得到Dang Van等效应力沿深度分布如图7所示。Dang Van等效应力在表面以下达到最大值，最大值深度位置约为500um，整条曲线先上升后下降。
[0073] 步骤5、根据公式(10)计算出材料的局部疲劳失效风险值，如图8所示，失效风险值在齿面以下达到最大值，最大值深度位置约为500um，最大失效风险值约为0.3，失效风险值较低。整条曲线先上升后下降，最大失效风险值在疲劳点蚀高发区。
[0074] “Tooth Flank Fracture–Basic Principles and Calculation Model for a Sub-Surface-Initiated Fatigue Failure Mode of Case-Hardened Gears”，I.Boiadjiev,J.Witzig,T.Tobie and K.Stahl,Gear Technology,pp.59-64,2015(“齿面断裂-表面硬化齿轮的次表面疲劳失效的基本准则与计算模型”，I.Boiadjiev,J.Witzig,T.Tobie和K.Stahl，齿轮技术，第59-64页，2015年)记载了考虑齿面硬化的失效风险值的实验，实验得到的沿深度分布的失效风险值曲线趋势与本发明得到的趋势基本一致，由此验证了本方法发明的实用性和可靠性。
