弧齿锥齿轮具有噪声小，传动效率高等优点，广泛用于航空航天、军用飞机，船舶等大型工业设备中。齿轮的 5 种失效类型当中，轮齿断裂发生率最大，而裂纹轮齿发生断裂的关键因素，它影响整个齿轮的可靠性和使用寿命。
杜国君等采用 ANSYS 有限元软件对 20CrMnTi 直齿轮进行不同渗碳浓度淬火数值分析，表明淬火后齿轮残余应力表现为“外压内拉”；高启林等研究直齿锥齿轮在不同的淬火介质下的温度场于应力场的变化规律;罗家元等研究淬火残余应力对铝合金厚板裂纹扩展的影响，表明淬火残余压应力对裂纹扩展起抑制作用;王延忠等采用 DEFORM 有限元软件研究不同淬火温度对轮齿残余应力、热变形的影响规律;侯兴隆采用 ABAQUS 有限元软件对 9310 钢直齿锥齿轮进行淬火数值模拟，采用模具压力加载方式等方法有效的控制齿轮整体变形量;刘双等基于 ABAQUS 软件研究不同初始裂纹长度下直齿轮齿根裂纹扩展特性;李强等通过 FRANC3D 软件研究不同因素下弧齿锥齿轮裂纹扩展规律，并求解应力强度因子，表明弧齿锥齿轮裂纹扩展以张开型为主，且当齿顶载荷增大时，不会改变原有裂纹扩展趋势;裴未迟等采用 ABAQUS 探究直齿圆柱齿轮的裂纹扩展情况，并求解应力强度因子;余洋等基于 XFEM 研究离心力对直齿轮裂纹扩展路径的影响规律;KATO 等研究残余应力对裂纹萌生和扩展的影响;LEWICKI 等预测各种轮齿和轮辋的裂纹扩展路径，分析初始裂纹位置、裂纹长度的影响规律。
在上述论述中，学者们主要聚焦于二维裂纹扩展和圆柱齿轮齿根裂纹扩展研究，对弧齿锥齿轮等复杂模型三维裂纹扩展研究较少，且未考虑残余应力对齿轮裂纹扩展的影响。因此，本文针对已有研究局限，以弧齿锥齿轮为研究对象，探究淬火残余应力对齿根裂纹扩展的影响规律，揭示残余应力作用下弧齿锥齿轮裂纹扩展机理。该研究成果为复杂模型三维裂纹扩展研究提供一种新思路，为进一步理解弧齿锥齿轮断裂失效提供参考与依据。
一、轮齿断裂有限元模型的建立
本文以弧齿锥齿轮为研究对象，由于 Abaqus 建立弧齿锥齿轮有限元模型和划分网格有一定的难度，因此，使用 SolidWorks 建立弧齿锥齿轮三维实体有限元模型并将其保存为stp.格式，齿轮几何参数如表 1 所示，将其导入 Hypermesh 中划分网格，轮齿接触部位的网格需细化，主动轮网格较为密集，减少沙漏现象，即 ALLAE /ALLIE<1% ，以确保分析结果的精准。
应力分析及危险位置确定：网格划分完成之后，将 inp.文件导入 Abaqus 中施加如下边界条件：
(1) 设置材料属性：密度为 7850 kg /m3，弹性模量为 207 GPa，泊松比为 0.3，单元类型为 C3D8R。
(2) 在大小齿轮轴孔中心建立参考点，并与之内表面相耦合。
(3) 分别约束大、小齿轮除绕轴旋转以外的所有自由度。
(4) 大齿轮施加 2220 rmp，小齿轮施加 200 N/m 的负载力矩。
通过有限元求得小齿轮应力云图如图 1 所示，在大齿轮转动过程中，啮合时小齿轮中间轮齿的齿根弯曲应力与接触力最大，由此可得，弧齿锥齿轮啮合时，中间轮齿发生断裂失效的概率高于其余轮齿，因此，以中间轮齿为研究对象，研究弧齿锥齿轮三维裂纹的扩展机理。
裂纹扩展分析：根据线弹性断裂力学理论，裂纹扩展有 2 种观点，即能量守恒观点与应力强度因子观点;能量守恒观点为裂纹扩展断裂的能量与裂纹新生面所需的能量相互平衡;应力强度因子观点为裂纹尖端应力理论值为无穷大，故以应力强度因子来表征裂纹尖端的弹性应力场。首先，基于扩展有限元法( XFEM) 探讨裂纹扩展路径。XFEM 以中心分解法为基础，在传统扩展有限元的基础上引入富集函数，其位移表达式中引入了跳跃函数与渐近位移场函数，并以此来处理裂纹扩展时单元界面不断演化而引起的复杂的非连续位移问题，其位移近似函数表达式描述如下：
式中，I 为所有单元节点集合；nc 为裂纹完全贯穿时所有单元节点集合；np 为裂纹尖端附近的所有单元节点集合；Ni (x) 为标准有限元单元形函数；ui 为标准单元节点位移矢量；ai 为裂纹完全贯穿时的位移矢量；bji 为裂纹尖端的位移矢量;H(x) 为跳跃函数，表达裂纹尖端的间断特性；φjnp为裂纹尖端渐近位移场分支函数，其表达式如下：
式中，
式中，r、θ 为裂纹尖端坐标原点的极坐标值。
裂纹扩展可分为 3 个阶段，分别为裂纹萌生阶段，裂纹稳定扩展阶段，裂纹失稳扩展阶段；基于上述应力场的分析结果，以啮合时小齿轮中间轮齿为研究对象，加载方式为单齿对啮合时的最高点加载，如图 2f 所示。本文设置如下 4 种工况研究裂纹扩展：
通过加载上述工况进行求解，图 2a～图 2d 为工况 4 的 4 个裂纹扩展过程;可以看出，当施加载荷时，裂纹先从齿根中部萌生，与应力场的分析结果相符合，且裂纹先沿齿面方向扩展，再沿深度方向扩展，裂纹沿齿面扩展速率大于裂纹沿深度扩展速率。图 2f～图 2h 为前 3 种工况下的裂纹扩展最终路径;当施加的载荷逐渐增大时，虽然裂纹扩展轨迹变化较大，但裂纹萌生部位均为齿根中部，最终贯穿整个齿根，扩展方向与工况 4 相同，区别在于，工况 1 的裂纹扩展路径较为平缓，从轮齿受拉一侧向受压一侧缓慢扩展，直至断裂失效；工况 2、3 的裂纹扩展路径较工况 1 相比，均先沿轮齿受压一侧扩展，再逐渐沿轮缘方向扩展，且工况 3 的扩展幅度大于工况 2，工况 4 的扩展幅度最大，最后断裂失效。
二、残余应力对齿根裂纹扩展的影响
弧齿锥齿轮淬火数值模拟：弧齿锥齿轮淬火过程属于复杂的非线性过程，包括多种耦合关系；淬火过程中，不同材料的齿轮与不同的淬火介质实现热交换的过程是非线性的，因此，只考虑温度对齿轮的影响。本文所采用齿轮材料为 20CrMnTi，淬火方式为分级淬火，即将初始温度为 850 ℃的齿轮放入淬火油中冷却至 80 ℃，然后在空气中冷却至 25 ℃。20CrMnTi 齿轮的弹性模量、比热容、淬火介质等塑性属性参数取自文献。20CrMnTi 平面应变断裂韧性为 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，齿芯应力表现为“先压后拉”。
应力强度因子的求解：应力强度因子(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 所示。
完成上述模型的建立，求解 2 种情形下的应力强度因子：
(1) 初始裂纹长度为0.2 mm 时，施加上述4 种载荷。
(2) 求解裂纹长度为 0.2 mm、0.4 mm、0.6 mm、0.8 mm，施加工况 4 为载荷，求解应力强度因子。
1.不同载荷对 SIFs 的影响
当初始裂纹长度为 0.2 mm 时，施加工况 3 作为载荷得到 K1、K2 有无残余应力的变化规律，如图 6 所示。
从图 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 曲线变化平稳。这表明，随着载荷增大，裂纹先沿齿面朝轮齿两端扩展，轮齿中间节点扩展速率最大，且裂纹在靠近大端沿齿面的扩展速率大于小端的扩展速率。
2.不同裂纹长度对 K1 的影响
为研究不同裂纹长度下残余应力对裂纹扩展的影响规律，施加工况 4 的载荷，研究裂纹长度为 0.2 mm、 0.4 mm、0.6 mm、0.8 mm 时 K1 的变化规律。
图 8 为不同裂纹长度下有无残余应力时 K1 的变化曲线。随着裂纹长度的增加，曲线仍保持中间大两头小的变化规律。考虑残余应力时，曲线峰值增幅较大，表明裂纹沿着深度的扩展速率随着裂纹长度的增大而增大。
图 9a ～ 图 9c 为轮齿大端节点、中间节点、小端节点随裂纹长度的变化曲线。
可以看到，中间节点与大端节点曲线呈线性变化;由于轮齿大、小端齿面载荷分布不均匀，大端承受载荷较大，随着裂纹长度的增加，小端节点先增大，后下降，考虑残余应力时小端应力强度因子的最大值是未考虑残余应力时的 1.92 倍，这表明，残余拉应力对轮齿基体裂纹扩展起促进作用，裂纹沿深度方向稳定扩展;弧齿锥齿轮大端齿厚大于小端齿厚，且裂纹在靠近轮齿大端沿深度的扩展速率大于小端的扩展速率，这有利于延长轮齿完全断裂失效时间，提高齿轮疲劳寿命。
三、结论
本文在已有研究成果的基础上，基于 XFEM 模拟了弧齿锥齿轮裂纹扩展路径，同时，探讨不同参数下淬火残余应力对应力强度因子的影响规律，得出以下结论：
(1) 裂纹萌生在齿根中部，先沿齿面扩展，在沿深度扩展，随着载荷的增大，逐渐朝轮缘扩展，且裂纹沿齿面的扩展速率大于沿深度的扩展速率。
(2) 考虑残余应力时，K1 值随着载荷的增大而增大，但不改变裂纹原有的扩展趋势;裂纹沿深度的扩展速率随着裂纹长度的增大而增大;裂纹在靠近轮齿大端的扩展速率总是大于小端的扩展速率。
参考文献略.