残余应力作用下弧齿锥齿轮裂纹扩展仿真分析*

2022-09-27罗家元李明学邓蒲实

组合机床与自动化加工技术 2022年9期

罗家元，李明学，邓蒲实

(重庆交通大学机电与车辆工程学院，重庆 400074)

0 引言

弧齿锥齿轮具有噪声小，传动效率高等优点，广泛用于航空航天、军用飞机，船舶等大型工业设备中。齿轮的5种失效类型当中，轮齿断裂发生率最大，而裂纹是轮齿发生断裂的关键因素，它影响整个齿轮的可靠性和使用寿命。

杜国君等[1]采用ANSYS有限元软件对20CrMnTi直齿轮进行不同渗碳浓度淬火数值分析，表明淬火后齿轮残余应力表现为“外压内拉”；高启林等[2]研究直齿锥齿轮在不同的淬火介质下的温度场于应力场的变化规律；罗家元等[3]研究淬火残余应力对铝合金厚板裂纹扩展的影响，表明淬火残余压应力对裂纹扩展起抑制作用；王延忠等[4]采用DEFORM有限元软件研究不同淬火温度对轮齿残余应力、热变形的影响规律；侯兴隆[5]采用ABAQUS有限元软件对9310钢直齿锥齿轮进行淬火数值模拟，采用模具压力加载方式等方法有效的控制齿轮整体变形量；刘双等[6]基于ABAQUS软件研究不同初始裂纹长度下直齿轮齿根裂纹扩展特性；李强等[7]通过FRANC3D软件研究不同因素下弧齿锥齿轮裂纹扩展规律，并求解应力强度因子，表明弧齿锥齿轮裂纹扩展以张开型为主，且当齿顶载荷增大时，不会改变原有裂纹扩展趋势；裴未迟等[8]采用ABAQUS探究直齿圆柱齿轮的裂纹扩展情况，并求解应力强度因子；余洋等[9]基于XFEM研究离心力对直齿轮裂纹扩展路径的影响规律；KATO等[10]研究残余应力对裂纹萌生和扩展的影响；LEWICKI等[11]预测各种轮齿和轮辋的裂纹扩展路径，分析初始裂纹位置、裂纹长度的影响规律。

在上述论述中，学者们主要聚焦于二维裂纹扩展和圆柱齿轮齿根裂纹扩展研究，对弧齿锥齿轮等复杂模型三维裂纹扩展研究较少，且未考虑残余应力对齿轮裂纹扩展的影响。因此，本文针对已有研究局限，以弧齿锥齿轮为研究对象，探究淬火残余应力对齿根裂纹扩展的影响规律，揭示残余应力作用下弧齿锥齿轮裂纹扩展机理。该研究成果为复杂模型三维裂纹扩展研究提供一种新思路，为进一步理解弧齿锥齿轮断裂失效提供参考与依据。

1 轮齿断裂有限元模型的建立

本文以弧齿锥齿轮为研究对象，由于Abaqus建立弧齿锥齿轮有限元模型和划分网格有一定的难度，因此，使用SolidWorks建立弧齿锥齿轮三维实体有限元模型并将其保存为stp.格式，齿轮几何参数如表1所示，将其导入Hypermesh中划分网格,轮齿接触部位的网格需细化，主动轮网格较为密集，减少沙漏现象，即ALLAE/ALLIE<1%，以确保分析结果的精准。

表1 齿轮几何参数

1.1 应力分析及危险位置确定

网格划分完成之后，将inp.文件导入Abaqus中施加如下边界条件：

(1)设置材料属性：密度为7850 kg/m3，弹性模量为207 GPa，泊松比为0.3，单元类型为C3D8R。

(2)在大小齿轮轴孔中心建立参考点，并与之内表面相耦合。

(3)分别约束大、小齿轮除绕轴旋转以外的所有自由度。

(4)大齿轮施加2220 rmp，小齿轮施加200 N/m的负载力矩。

通过有限元求得小齿轮应力云图如图1所示，在大齿轮转动过程中，啮合时小齿轮中间轮齿的齿根弯曲应力与接触力最大，由此可得，弧齿锥齿轮啮合时，中间轮齿发生断裂失效的概率高于其余轮齿，因此，以中间轮齿为研究对象，研究弧齿锥齿轮三维裂纹的扩展机理。

图1 小齿轮应力云图

1.2 裂纹扩展分析

根据线弹性断裂力学理论，裂纹扩展有2种观点，即能量守恒观点与应力强度因子观点；能量守恒观点为裂纹扩展断裂的能量与裂纹新生面所需的能量相互平衡；应力强度因子观点为裂纹尖端应力理论值为无穷大，故以应力强度因子来表征裂纹尖端的弹性应力场。首先，基于扩展有限元法(XFEM)探讨裂纹扩展路径。XFEM以中心分解法为基础，在传统扩展有限元的基础上引入富集函数,其位移表达式中引入了跳跃函数与渐近位移场函数，并以此来处理裂纹扩展时单元界面不断演化而引起的复杂的非连续位移问题，其位移近似函数表达式[13]描述如下：

(1)

(2)

式中，

式中，r、θ为裂纹尖端坐标原点的极坐标值。

唐飞霄咬了咬牙，这八只节足中，最灵活的便是两只前足，此刻两足先后而断，蜘蛛机甲的攻击力无疑大大受损。他望着天葬师，瞳中闪过一抹戾色，而后身形下压，伴着一阵“嘎啦嘎啦”的响声，蜘蛛的六只节足，开始出现了变化。那些节足本是三段结构，中间以关节相连，此刻，每一段都开始拉伸，各自又分化为三小段，这令整只节足，转眼变成了九段结构，段与段间的连接，都有一扎长的利刃裸露在外。

裂纹扩展可分为3个阶段，分别为裂纹萌生阶段，裂纹稳定扩展阶段，裂纹失稳扩展阶段；基于上述应力场的分析结果，以啮合时小齿轮中间轮齿为研究对象，加载方式为单齿对啮合时的最高点加载，如图2f所示。本文设置如下4种工况研究裂纹扩展：

工况1：Ft=4000 N，Fr=1000 N
工况2：Ft=5000 N，Fr=2000 N
工况3：Ft=6000 N，Fr=2000 N
工况4：Ft=7000 N，Fr=3000 N

通过加载上述工况进行求解，图2a～图2d为工况4的4个裂纹扩展过程；可以看出，当施加载荷时，裂纹先从齿根中部萌生，与应力场的分析结果相符合，且裂纹先沿齿面方向扩展，再沿深度方向扩展，裂纹沿齿面扩展速率大于裂纹沿深度扩展速率。图2f～图2h为前3种工况下的裂纹扩展最终路径；当施加的载荷逐渐增大时，虽然裂纹扩展轨迹变化较大，但裂纹萌生部位均为齿根中部，最终贯穿整个齿根，扩展方向与工况4相同，区别在于，工况1的裂纹扩展路径较为平缓，从轮齿受拉一侧向受压一侧缓慢扩展，直至断裂失效；工况2、3的裂纹扩展路径较工况1相比，均先沿轮齿受压一侧扩展，再逐渐沿轮缘方向扩展，且工况3的扩展幅度大于工况2，工况4的扩展幅度最大，最后断裂失效。

(a) 裂纹萌生阶段 (b) 裂纹沿齿面稳定扩展

2 残余应力对齿根裂纹扩展的影响

2.1 弧齿锥齿轮淬火数值模拟

弧齿锥齿轮淬火过程属于复杂的非线性过程，包括多种耦合关系；淬火过程中，不同材料的齿轮与不同的淬火介质实现热交换的过程是非线性的，因此，只考虑温度对齿轮的影响。本文所采用齿轮材料为20CrMnTi，淬火方式为分级淬火，即将初始温度为850 ℃的齿轮放入淬火油中冷却至80 ℃，然后在空气中冷却至25 ℃。20CrMnTi齿轮的弹性模量、比热容、淬火介质等塑性属性参数取自文献[1]。20CrMnTi平面应变断裂韧性[14]为128 MPa·m1/2。在ABAQUS中采用顺序热力耦合进行模拟，齿轮2个截面采用对称约束边界条件，其余面为对流面，将对流换热系数作为热边界条件施加在对流面，分析步时间为500 s，油淬时间200 s，空气中冷却300 s，齿轮淬火模型及边界条件的施加如图3a所示，图中A、B、C分别代表齿顶、齿根以及齿芯部位。

轮齿在淬火介质中的各点冷却曲线如图3b所示，从中可以看出，3个特征点的冷却曲线各不相同，齿根与齿顶冷却速率大于内部的冷却速率；油淬阶段冷却至200 s时结束，最低温度位于齿顶部位，且曲线变化较陡，空冷阶段冷却速度较为平缓，在300 s时冷却至30 ℃，最高温度位于齿心；在淬火冷却初期，轮齿各个位置急速冷却，且齿顶部位的冷却速度大于齿根和齿心部位，在4 s时，齿顶温度达到310 ℃，而齿根与齿心较为接近，达到365 ℃，伴随着冷却的继续进行，三者冷却速度更为接近，逐渐到达80 ℃。

图3c为特征点B沿x、y、z方向的应力变化曲线。从中可知，淬火冷却初期，因轮齿齿根及齿面的冷却速率高于内部，齿面向内收缩，受到拉应力作用，随着表面冷却速率的下降，芯部与表面相反；轮齿齿根表面产生拉应力时，单向拉应力最高达到335 MPa，而后下降转变为压应力，最大单向压应力达到413 MPa，随着冷却曲线上升且在100 s左右达到平稳，此时，齿根表面残余压应力最大值为239 MPa，最小值为78 MPa，齿根表面应力表现为“先拉后压”。

图3d为特征点C沿x、y、z方向的应力变化曲线。齿心应力变化与齿根相反，淬火起始时，内部因阻碍表层组织向内收缩而产生压应力，之后，轮齿表面冷却速率下降，而芯部的冷却速率高于表面，此时，内部收到拉应力。齿芯产生压应力时，单向压应力最高达到408 MPa，而后曲线上升转变为拉应力时，拉应力最大为430 MPa，随着冷却曲线下降，在100 s左右时达到稳定值，单向拉应力最大值为269 MPa，最小值为192 MPa，齿芯应力表现为“先压后拉”。

(a) 齿轮淬火模型 (b) 各点冷却曲线

2.2 应力强度因子的求解

应力强度因子(SIF)是反映弹性体裂纹尖端的应力场强度的物理量，是评价裂纹扩展的重要指标，其有3种类型：张开型K1、滑开型K2、撕开型K3。在上文中利用XFEM研究了弧齿锥齿轮的裂纹扩展，得到了裂纹扩展路径，可知裂纹由齿根中部萌生逐渐贯穿整个齿根，直至断裂，因此，为更进一步探究裂纹尖端的应力场变化以及残余应力对齿根裂纹扩展的影响，建立含初始裂纹模型，模型网格如图3a所示，模型处理如下。

奇异单元的构建：利用ABAQUS计算裂纹尖端应力强度因子时，通常会将裂纹尖端附近区域采用奇异单元划分。奇异单元是在ABAQUS的GUI中构建，因此,尽管裂纹尖端采用三角形单元划分，ABAQUS会通过自动调整单元节点来完成奇异单元的构建，但裂纹尖端的网格如果是从外部(Hypermesh等)导入的孤立网格，ABAQUS会将其视为三角形网格，而非奇异等参单元，而围线积分不能采用三角形单元计算应力强度因子，因此，首先将需要求解的网格转换为奇异单元。转换方法为将hypermesh画好的网格导出为inp格式，通过查找裂纹尖端的单元节点，将六面体网格被折叠的边线上的节点复制即可，在ABAQUS中将单元类型设置C3D20，在网格编辑功能区选择调整中间节点，将其中间节点调整为靠近裂纹尖端的1/4处，获得r1/2的奇异性，如图4所示。

载荷、边界条件与裂纹参数：载荷施加上述4种工况，载荷施加位置为单齿对啮合时的最高点处，轮齿2个截面施加对称边界约束条件，约束内表面6个自由度，q为裂纹扩展方向，a为裂纹长度，本文设定初始裂纹长度为0.2 mm，采用围线积分法定义三维裂纹如图5所示。

图4 奇异等参单元图5 三维裂纹模型

完成上述模型的建立，求解2种情形下的应力强度因子：

(1)初始裂纹长度为0.2 mm时，施加上述4种载荷。

(2)求解裂纹长度为0.2 mm、0.4 mm、0.6 mm、0.8 mm，施加工况4为载荷，求解应力强度因子。

2.2.1 不同载荷对SIFs的影响

当初始裂纹长度为0.2 mm时，施加工况3作为载荷得到K1、K2有无残余应力的变化规律，如图6所示。

(a) K1的分布 (b) K2的分布

从图6a中可以看出，轮齿齿根中部的应力强度因子最大，考虑残余应力时最大值为1679 MPa·mm1/2，最小值为179.2 MPa·mm1/2,未考虑残余应力时最大值为1598 MPa·mm1/2，最小值为91.96 MPa·mm1/2，二者最大值相差81 MPa·mm1/2，最小值相差87.24 MPa·mm1/2，可以看出二者差值相差较大，且最大值均位于轮齿中部，轮齿两端的应力强度因子较小，且小端比大端小，即中间拉应力大于轮齿两端拉应力，从而曲线为拱形，与上述裂纹扩展结论相符合；图6b中未施加残余应力的曲线逐渐增大，呈线性分布，最大值为148.7 MPa·mm1/2，最小值为-180.8 MPa·mm1/2，考虑残余应力时，最大值为133 MPa·mm1/2，最小值为-153.4 MPa·mm1/2。从两图中可以得到，当初始裂纹长度相同时，齿面齿根部位残余应力主要表现为“外压内拉”，而齿芯部位主要表现为内部拉应力，裂纹受到两种不同情形的应力，且两种应力分布并非绝对，故残余应力对于齿根裂纹扩展有两种影响，即裂纹前缘节点受残余压应力时，抑制裂纹扩展，裂纹前缘节点受残余拉应力时，促进齿根裂纹扩展。由图可知，K1远大于K2，且K2的最大值为K1最大值的8%，所以，对裂纹扩展起决定性作用的是K1，故弧齿锥齿轮裂纹扩展以张开型为主。

为得到不同载荷作用下K1的变化规律，施加其余3种工况，得到图7a中K1的分布，可知，无论载荷如何变化，K1的最大值总位于齿根中部，且裂纹前缘应力强度因子随着载荷的增大而增大，但曲线变化趋势相同，考虑残余应力时，应力强度因子也随之增大，曲线保持中间大，两头小的变化规律。

图7b～图7d为轮齿大端节点、中间节点、小端节点随载荷的变化规律。考虑残余应力时，大端节点应力强度因子最大值为768.4 MPa·mm1/2，小端最大值为179.2 MPa·mm1/2，大端为小端的4.29倍；3个节点的K1曲线变化平稳。这表明，随着载荷增大，裂纹先沿齿面朝轮齿两端扩展，轮齿中间节点扩展速率最大，且裂纹在靠近大端沿齿面的扩展速率大于小端的扩展速率。

(a) 其余3种工况K1的分布 (b) 大端节点

2.2.2 不同裂纹长度对K1的影响

为研究不同裂纹长度下残余应力对裂纹扩展的影响规律，施加工况4的载荷，研究裂纹长度为0.2 mm、0.4 mm、0.6 mm、0.8 mm时K1的变化规律。

图8为不同裂纹长度下有无残余应力时K1的变化曲线。随着裂纹长度的增加，曲线仍保持中间大两头小的变化规律。考虑残余应力时，曲线峰值增幅较大，表明裂纹沿着深度的扩展速率随着裂纹长度的增大而增大。

图8 不同裂纹长度下K1的变化规律

图9a～图9c为轮齿大端节点、中间节点、小端节点随裂纹长度的变化曲线。

(a) 大端节点 (b) 中间节点

可以看到，中间节点与大端节点曲线呈线性变化；由于轮齿大、小端齿面载荷分布不均匀，大端承受载荷较大，随着裂纹长度的增加，小端节点先增大，后下降，考虑残余应力时小端应力强度因子的最大值是未考虑残余应力时的1.92倍，这表明，残余拉应力对轮齿基体裂纹扩展起促进作用，裂纹沿深度方向稳定扩展；弧齿锥齿轮大端齿厚大于小端齿厚，且裂纹在靠近轮齿大端沿深度的扩展速率大于小端的扩展速率，这有利于延长轮齿完全断裂失效时间，提高齿轮疲劳寿命。