齿轮作为传动系统的核心部件，被广泛运用于航空航天、汽车坦克和高精密仪器制造等领域，喷丸强化是高性能齿轮制造的关键工序，通过齿根喷丸强化引入残余压应力可以有效提高齿轮使用寿命。喷丸与其他机械加工方法一样，最终零件表面是凹凸不平的微观形貌，这些微观形貌的谷底区域在工件工作时会产生应力集中，表面粗糙度越大，产生的应力集中越严重，工件的疲劳寿命下降越快。因此，精确计算齿根喷丸粗糙表面的应力与应力集中系数对评估齿轮弯曲疲劳寿命极为关键。
对于粗糙表面与应力集中系数的关联性研究，NEUBER将粗糙表面轮廓近似为连续相邻的缺口, 并由缺口深度和宽度是影响缺口应力集中的主要因素，推广到基于粗糙表面参数 Rz 和谷底曲率半径 ρ 的应力集中系数计算公式。AROLA 等依据 NEUBER 提出的经验公式建立了 Arola-Ramulu 模型，将实际的表面形貌近似为连续相邻的正弦状缺口展开研究，并通过实验验证了该模型的准确性。CHAN对电子束加工得到的板块表面进行测量，研究了具有凹痕特征的形貌的应力集中，提出了加工制造过程中对表面形貌的要求;KANTZOS借助神经网络，将应力集中和表面形貌进行关联分析，实现了裂纹萌生的预测。在数值计算方面，学者们通过对具有特定规则的形貌进行简化研究，采用一阶边界摄动法，基于 Hilbert 变换和能量守恒原理，预测了拉伸作用下二维表面形貌上的应力分布并进行实验验证。然而，目前可见文献对于粗糙表面与应力集中系数的关联规律大多基于二维表面形貌展开研究，对于三维实测表面形貌引起的应力集中研究较少，MEREUTA分别将 Neuber 模型和 Arola-Ramulu 模型的公式中二维粗糙表面参数转换为对应的三维粗糙表面参数进行计算，并经过试验验证了修正公式的可靠性。LI基于表面形貌重构技术对重构表面引起的应力集中进行分析，利用少数几个重构平面模型分析了粗糙表面参数 Sa 疲劳寿命之间的关系。XU研究了腐蚀对 Q235 钢板疲劳寿命的影响，借助 ANSYS 软件生成了三维粗糙形貌并分析了三种腐蚀坑与疲劳寿命的关系，最后与疲劳实验进行对比验证。在 GB/T 3480—1997 中通过大量的疲劳实验提出了相对齿根表面状况系数经验公式,然而该经验公式没有给出粗糙度表面应力集中系数的计算表达式，而是通过相对齿根应力集中系数的比值拟合的公式间接获得应力集中对齿轮弯曲疲劳寿命的影响，且该经验公式仅用 1 个二维粗糙度参数 Rz 进行表征，显然无法精确评估三维实际表面产生的应力集中。
随着对三维表面的认识，用于表征粗糙表面形貌特征的三维粗糙度参数多达几十个，不同的粗糙度参数表征表面的不同性能，而对于表征应力集中的三维粗糙度参数，目前尚未得到共识。且无论是二维还是三维粗糙表面问题，现阶段所有研究对象均是以平板或者简化为平面模型进行研究，解决齿根粗糙表面下实际应力精确计算，得到表征粗糙表面应力集中合适的形貌参数问题是高端齿轮、高性能齿轮设计制造的核心技术之一，具有较大的理论与工程应用价值。
论文以磨削后喷丸的直齿轮为研究对象，提取齿轮齿根三维粗糙表面形貌信息，依据空间坐标变换原理，研究齿根过渡曲面添加三维粗糙面的建模方法，从而突破现有齿轮强度分析只能创建光滑表面模型的限制，建立基于实测三维表面形貌的齿根区域有限元模型，对有限元计算得到的数据进行非线性曲线拟合，得到齿根三维粗糙表面参数与应力集中系数 Kt 的关联公式，为后续的齿轮弯曲疲劳寿命精准预测和齿轮表面完整性参数设计提供相应的参考。
一、齿根表面形貌测量
为研究齿根喷丸表面危险区域的表面粗糙度与应力集中系数的关系，首先确定齿根危险区域的具体位置，基于磨削喷丸加工的齿轮(表 1)，借助 Matlab 编程生成直齿轮齿廓，将齿廓导入 CATIA，通过阵列、拉伸等操作生成直齿轮实体模型，将齿轮实体模型导入 ABAQUS 中进行准静态仿真分析，查看仿真结果应力云图，确定了齿根危险区域(应力较大区域)。
其次，测量表面形貌。本文借助 ABAQUS 研究齿根危险区域表面粗糙度对应力集中的影响，通过调整网格节点坐标值实现粗糙度的添加，而齿根危险区域的网格数量取决于形貌测量的采样点数，因此，选择较为合适的采样间隔减小网格数量、提高仿真效率是很有必要的。文中采用白光干涉仪 Wyko NT9100 对粗糙形貌进行测量。对于表面形貌的测量，采样间隔越小，表面形貌还原越完整。考虑文中对有限元网格数量的要求，在保留形貌基本特征前提下，适当选择较大的采样间隔，舍弃部分形貌信息。一般地，采样频率大于形貌频率 4 倍时，采样离散数据点可以基本复原形貌的大致轮廓特征。
文中以喷丸后的粗糙表面为研究对象，主要研究喷丸处理后产生的丸坑形貌对应力集中的影响。一般地，表面形貌的应力集中最有可能产生在粗糙形貌最深的几个谷底区域之一，已知齿轮喷丸直径约为 120 μm，喷丸后产生的几个最深丸坑深度平均值绝大部分在5 μm以上(详细数据见表2中的 S10z )，则丸坑的表面直径大于 40 μm。因此，在满足采样频率大于 4 倍的形貌频率(丸坑表面直径)的基础上，论文测量时选择采样间隔为 10 μm，采样面积为1 mm ×1 mm。得到磨削后喷丸表面形貌如图 1 所示。
二、考虑三维粗糙度的齿根应力有限元计算
齿轮网格划分及细化
齿根应力有限元计算需要精确的齿根形貌与齿面形貌，渐开线齿面的计算按渐开线方程计算，齿根过渡曲面的构建按课题组前期提出的建模方法构建，详细齿根过渡曲面构建方法可参考文献。
本文通过调整网格节点坐标实现粗糙度的添加工作，在齿根区域的网格节点数目应与采样离散点的数目一致。考虑直齿轮实体模型是由齿廓方程生成的，沿齿廓方向各点的受力状态基本相同。因此在建立齿轮有限元模型时将表 1 直齿轮参数的齿宽设定为与采样面积宽度一致的 1 mm，大大降低重复的网格数量。此外，为了验证齿宽为 1 mm 时齿轮有限元模型的可靠性，考虑模型边缘效应等问题，分别对齿宽为 1 mm、5 mm、10 mm 的有限元模型进行仿真，保证其他边界、材料等条件相同，载荷扭矩分别为T、5T 、10T 进行仿真，并对仿真结果进行对比分析，仿真结果表明三种齿宽条件下齿根弯曲应力的最大值相差在 5%以内，证实 1 mm 齿宽条件下齿轮模型的可靠性。其次将齿轮切分成两部分，一部分包含齿根危险区域，需要进行网格细化，另一部分不需要调整网格节点，因此，网格数量不做要求。切分操作在 CATIA 中完成，具体切分过程如图 2 所示。
将切分后的齿轮两部分分别导入 Hypermesh 中进行网格划分。对齿根危险区域网格进行局部细化，齿根细化区域网格尺寸为10 μm ×10 μm，细化层数为六层。同时，为保证网格质量，在齿轮啮合区域网格与齿根危险区域网格的尺寸相差不宜过大，未包含齿根危险区域部分齿轮对仿真过程关注的弯曲应力无影响，网格尺寸可尽量调大，最后划分得到有限元网格数约为 25 万，直齿轮的网格划分模型如图 3 所示。
齿根危险区域表面粗糙形貌的生成
将 Hypermesh 划分好的齿轮模型保存为 inp 文件导入 ABAQUS 中，借助 Python 对 ABAQUS 进行二次开发实现齿根粗糙形貌的添加，具体操作步骤如下所示。
(1) 删除直齿轮有限元模型中非齿根细化区域网格，保留齿根网格细化区域。
(2) 对齿根网格细化区域表层创建节点集，创建节点集后删除表层网格选择第二层网格表面创建节点集，依次创建不同层数的节点集，此处论文创建六层节点集。
(3) 利用 Python 对 ABAQUS 进行二次开发，提取步骤(2)中创建的节点集的节点坐标信息;
(4) 由步骤(3)获得的节点坐标信息是杂乱无序的，对节点坐标依次按照节点沿齿廓方向以及节点高度值大小进行排序，获得规则的网格节点编号以及对应节点坐标;
(5) 基于空间坐标变换原理，将采样离散点高度信息以垂直于齿根危险区域曲面的方式分别对应添加到 ABAQUS 提取的节点集中，即改变相应节点的坐标值。具体变换过程如图 4 所示。
如图 4 所示，直齿轮由齿廓曲线沿齿宽方向拉伸生成，因此用平面齿廓曲线代替直齿轮模型。其中，P0 为齿廓上齿根危险区域的一点，该点坐标由步骤(4)获得，P0 点对应曲线上的法向量
由齿廓方程得到，ΔZi 为采样离散点高度值，P1 为 P0 沿其对应法向量 n 方向上添加形貌高度值 ΔZi 后的坐标。其中，第一层网格节点垂直曲面高度增加值为 ΔZi 的 6/6 ，第二层网格节点垂直曲面高度增加值为 ΔZi 的 5/6 ，第三层增加值为 4/6 ，依次类推，直至离散点高度信息全部添加完成。
(6) 将修改完成的网格节点信息替换前面的齿轮模型 inp 文件对应的网格节点，至此，齿根危险区域表面粗糙形貌生成。
生成的齿轮形貌示意如图 5 所示。
直齿轮的齿根应力仿真计算
为了获得粗糙表面参数与应力集中系数的关联规律，根据 ISO 25178 标准，选择其中几个常用的空间粗糙度表征参数进行计算分析：算术平均高度 Sa 、最大谷深 Sv 和十点区域高度 S10z 。各参数的定义如下
式中，z 为采集形貌离散点的高度，i和 j 分别表示喷丸形貌的凸峰和凹谷。三维粗糙表面参数 Sa 、Sv 、S10z 与之对应的二维粗糙表面参数为 Ra、Rv、Rz，二维粗糙表面参数是在三维粗糙表面上截取一条线段进行计算，因此，二维粗糙表面参数会丢失部分形貌信息，为了更完整地研究粗糙表面产生的应力集中，论文选择三维粗糙表面参数进行形貌表征。
完成齿根危险区域添加表面粗糙度的工作后，论文选择三组齿轮副的啮合模型作为研究对象，其中小齿轮的中间齿轮为添加粗糙度的齿轮，两边齿轮则未添加。建立耦合和接触对，将切分的包含粗糙度的中间齿轮用绑定约束进行连接，添加载荷并设置约束条件，采用准静态分析方法进行仿真分析，齿轮对的有限元模型如图 6 所示。
为了研究齿根危险区域表面粗糙度与应力集中系数的关联规律，本文采用的应力集中系数定义如下
式中，σmax 表示带有缺口模型的缺口处应力最大值，σn 表示光滑试样无缺口的名义应力。此处取σmax 为添加粗糙度齿轮模型齿根弯曲应力最大值，σn 为无粗糙度时齿根弯曲应力最大值。由该定义可知，通过两者的比值可消除 1 mm 齿宽与实际齿宽不同带来的影响。
将添加粗糙度的齿轮仿真得到的齿根弯曲应力最大值与未添加粗糙度的齿轮弯曲应力最大值相比，即得到该粗糙表面参数下对应的应力集中系数 Kt 。有无粗糙度的齿轮弯曲应力分布对比如图 7 所示。为了更清晰地看出齿根粗糙度表面的应力分布，只选取齿轮切分包含齿根危险区域部分的进行展示。
如图 7 所示，带有表面粗糙形貌的齿根应力分布与无粗糙度时的齿根应力分布相差较大，粗糙表面的齿根应力集中主要分布在表面形貌的波谷处，符合实际情况。取 8 个形貌样本进行仿真分析，得到各粗糙表面参数与应力集中系数之间的仿真结果如表 2 所示。
三、粗糙表面参数与齿根应力集中系数关联规律
不同于现有对平面粗糙形貌应力集中研究的计算方法，论文提出了一种考虑齿根过渡曲线宏观几何结构上的微观形貌产生的应力集中求解方法。为了验证论文所研究的求解方法，将论文的仿真结果与 GB/T 3480—1997 中的相对齿根表面状况系数 YRrelT (通过试验获得)进行比较。
在 GB/T 3480—1997 中，齿根表面状况系数的定义是“考虑齿轮根部的表面状况，主要是齿根圆角(即齿根危险区域)处的表面粗糙度对齿根弯曲强度的影响”，齿根表面状况系数即为齿根粗糙表面产生的应力集中系数。而相对齿根表面状况系数定义是“所计算齿轮的齿根表面状况系数与试验齿轮的齿根表面状况系数的比值”，给出的计算经验公式是
式中，Rz 为二维粗糙表面参数中的微观不平度十点高度，Rz 对应论文中的十点区域高度 S10z，基于 S10z 的相对齿根表面状况系数公式为
论文仿真计算的表面粗糙度产生的应力集中系数(齿根表面状况系数)，与相对齿根表面状况系数的定义不同，因此，不能直接将粗糙表面参数代入求解，需要进行转换计算，转换计算的原理方法如下：将国标中定义的试验齿轮(粗糙表面参数 S10z =11.27 μm )对应的应力集中系数除以所需要计算齿轮的应力集中系数即为该计算齿轮的相对齿根表面状况系数。试验齿轮的粗糙表面参数 S10z 与样本 8 基本相同。因此，应力集中系数按样本 8 数据进行计算，经转换计算，表 2 中各粗糙度对应的仿真计算获得的相对齿根表面状况系数为样本 8 的应力集中系数除以对应样本的应力集中系数。将各个样本对应的十点区域高度 S10z 代入式(6)中，即得到国标经验公式所计算的YRrelT 值。对比结果如图 8 所示。
如图 8 所示，仿真计算得到的相对齿根表面状况系数与国标经验公式计算得到的系数整体趋势相同，且两者计算得到的系数值也十分接近。此外，在国标中注释了经过强化处理(如喷丸)的齿轮，其相对齿根表面状况系数要稍微大于计算经验公式所确定的数值，从图 8 整体趋势不难看出仿真值稍微大于经验公式计算的数值，与国标中的注释说法一致。图 8 从侧面论证了论文仿真求解方法的可信度。此外，国标经验公式间接表示了齿轮齿根处形貌引起的应力集中，而本文则直接求解对应形貌的应力集中，表现更加直观、具体。此外，与实验获得应力集中时耗时耗力相比，通过论文可以快速、低成本计算出应力集中系数，为后续精确计算齿轮弯曲疲劳寿命提供技术支持。
为了得到喷丸表面三维粗糙表面参数与应力集中系数的关联规律，借助 Matlab 软件的 cftool 工具包对表 2 中的数据进行关系式拟合，得到拟合关系如图 9 所示。
图 9a、9b、9c 分别表示粗糙表面参数 Sa 、Sv 、S10z 与应力集中系数 Kt 的拟合曲线图，对应的拟合相关系数 R 分别为 0.799，0.784，0.914。在对 Kt 的拟合关系中，十点区域高度 S10z 的拟合效果最好，拟合公式如式(7)所示
式(7)的拟合公式可以较好表征粗糙表面形貌与应力集中系数的定量关系。需要说明的是，加工零部件图纸上标注的粗糙表面参数 Ra 是粗糙表面的二维表征参数，表征三维粗糙表面具有一定的局限性，选择三维粗糙表面参数更为合适。从论文的分析结果看，为更加科学的表征三维粗糙表面对应力集中的影响，使用十点区域高度参数 S10z 更加科学，工程应用中，可优先使用 S10z 对应力集中系数进行计算。
四、结论
(1) 本文提出了一种基于齿根三维粗糙表面形貌的有限元建模与分析方法，解决了三维粗糙齿根表面应力集中系数计算问题，为齿轮表面完整性与抗疲劳设计研究提供了一种基础方法。
(2) 论文采用有限元计算得到的表面粗糙度引起的应力集中系数与 GB/T 3480—1997 中给出的相对齿根表面状况系数经验公式进行间接对比，两者计算结果基本相同。
(3) 对有限元计算得到的数据进行拟合分析，结果表明粗糙表面参数 Sa 、Sv 、S10z 拟合应力集中系数的拟合公式的相关系数分别为 0.799，0.784， 0.914，十点区域高度参数 S10z 能较好地表征齿根表面的应力集中，工程应用中可用十点区域高度 S10z 对喷丸粗糙表面产生的应力集中进行计算与表征。
参考文献略.