早期，众多研究人员利用有限元法( Finite Element Method，FEM)来解决裂纹扩展问题。LEWICKI基于FEM开展大量仿真计算，研究了轮缘厚度对齿轮裂纹扩展寿命和路径的影响，发现轮缘厚度增加易导致齿轮轮缘断裂，从而缩短齿轮寿命。除轮缘厚度外，初始裂纹位置与初始裂纹方向也是影响裂纹扩展路径和寿命的重要因素。与初始裂纹方向相比，初始裂纹位置影响更加显著。由于裂纹更新过程中需要不断重划分网格，因此FEM求解效率较为低下。为解决这一问题，扩展有限元法(Extended Finite Element Method，XFEM)被广泛应用，以实现不连续场描述与网格边界的分离，从而避免了网格重划分，显著提升了计算效率。后续研究普遍采用XFEM解决裂纹扩展问题。ＲAGHAV等建立了齿根裂纹扩展有限元模型，通过与实测裂纹扩展路径和扩展寿命进行对比，验证了计算方法的有效性。 

然而，上述研究仅局限于齿轮几何、结构参数等因素，没有考虑粗糙度的影响。实际上，粗糙度也是影响齿轮寿命的重要因素。粗糙度的变化影响着工件的疲劳极限，且粗糙度在构件中的不同位置会导致材料疲劳性能形成较大差异，进而影响裂纹扩展速率和方向。LI等通过试验对比，证明了粗糙度对裂纹扩展路径有较大影响。贺迪等的研究表明，齿轮弯曲疲劳寿命随表面粗糙度Ｒa的增大而降低。然而，这些研究主要集中在粗糙度对应力集中的影响，以及应力集中如何导致表面裂纹的萌生。关于粗糙度与裂纹扩展之间的关联规律，相关研究报道少见。 

为从理论上预测表面粗糙度对齿根裂纹扩展的影响，本文提出一种考虑表面粗糙度的齿根裂纹扩展计算方法。对粗糙表面重构所得粗糙度高度序列，先以最小拟合误差为原则进行拟合，再通过空间坐标变换对齿根轮廓节点进行坐标调整，构建齿根粗糙齿面几何模型。为兼顾计算精度与效率，对网格尺寸和粗糙度采样长度的选取进行分析。基于XFEM，借助Fortran编制裂纹扩展计算程序。为展示所提方法，应用所编制程序，对不同粗糙度下齿根裂纹扩展路径和扩展寿命进行研究。

在传统FEM中，位移场uFEM通过式(1)所示的形函数与单元节点位移来表达。 

式中:D为所有节点集合；i为节点序号；Ni(x)为形函数；ui为节点位移。XFEM通过在单元内引入具有不连续特性的扩充函数H(x)和F(x)，用于描述位移场的不连续性和应力奇异性，其位移场uXFEM可表示为

式中:Dh和Dc分别为由H(x)与F(x)扩充的节点集合；j和k为扩充节点的序号；cj和dk分别为裂纹贯穿单元和裂尖单元的节点附加自由度；Nj(x)和Nk(x)为扩充单元形函数。扩展有限元法相关原理可参看文献。

二维齿根裂纹可归为Ⅰ、Ⅱ复合型裂纹扩展问题。依据最大周向拉应力准则，裂纹沿最大周向拉应力方向扩展至失稳断裂，扩展角度∂和等效应力强度因子Keq分别为

式中:KⅠ、KⅡ分别为Ⅰ型、Ⅱ型应力强度因子，可由裂尖应力场经相互作用积分确定，计算过程参看文献。 应用当前广泛使用的Paris公式，得到裂纹扩展速率与应力强度因子幅度的关系，如式(5)所示。

式中:a为裂纹长度；FN为裂纹扩展寿命；C、m为试验确定的常数；ΔKeq为应力强度因子幅度，对于齿轮脉动循环加载，ΔKeq等于最大应力强度因子Kmax。

考虑到有限元软件难以直接输出应力强度因子、裂纹扩展方向角等参数，导致疲劳裂纹扩展路径和扩 展寿命的预测存在困难。为此，课题组基于XFEM，借助Fortran编制裂纹扩展计算程序，为后续裂纹扩展分析提供了计算手段。 

实测表面形貌易受测量方向、采样长度和表面缺陷影响，导致数据误差较大。为解决这一问题，数值重构方法展现出明显优势。它能突破测量范围限制，避免随机误差，并可高效批量生成结果，显著节省成本和时间。因此，本文基于随机过程理论，利用采样高度和自相关分布，开展粗糙表面重构。相关重构算法可参看文献。通过给定不同统计分布参数，模拟得到不同统计分布下的粗糙度高度序列Z'。为构建齿根粗糙齿面几何模型，将Z'与光滑齿面叠加。由于Z'起伏较大，为保证网格划分质量，参照文献的处理方法，采用抛物线方程对Z'进行拟合，得到较为平滑的粗糙度轮廓Z。拟合前后粗糙度高度序列的相对误差为

式中:ε为相对误差；N为Z'的总点数；xi为第i个点的位置；z'i为拟合前第i个点的高度；f(xi) 为利用抛物线方程拟合后所得高度值。

对拟合所得高度序列Z，通过空间坐标变换对齿根轮廓节点坐标进行调整，调整方法如图1所示。图1中，P0为齿廓渐开线上任意一点，Zi为该点高度序列值，坐标调整后变为点P1。为避免齿根网格坐标调整幅度过大而影响网格质量，构建多层网格节点，各层节点坐标按公式(7)计算，详见参考文献。

图1 齿根轮廓节点坐标调整示意图

记P0=(x0，y0)，P1=(x1，y1)，法向量n=(exi，eyj)，则网格节点坐标调整计算公式如下:

式中:x0、y0、x1、y1分别为点P0和P1的横、纵坐标；exi、eyj为方向矢量。齿轮网格节点坐标调整示意图如图2所示。

图2 齿轮网格节点坐标调整示意图

建立齿根粗糙齿面几何模型后，在齿根插入初始裂纹，得到考虑表面粗糙度的齿根弯曲疲劳裂纹扩展几何模型。其中，齿根区域网格细化尺寸和高度序列的采样间隔均为1μm。依据前述裂纹扩展计算原理，基于XFEM，借助Fortran编制裂纹扩展计算程序。限于试验条件，为间接验证编制程序的有效性，对光滑齿面与较小粗糙度(轮廓均方根偏差Ｒq=0.05)下的齿根裂纹扩展进行对比分析。计算所得应力强度因子与裂纹扩展寿命如图3所示。不难发现，两者应力强度因子与裂纹扩展寿命较为接近。一方面，粗糙度的存在导致应力强度因子增大，裂纹扩展寿命减小，这与预期结果相吻合。另一方面，由于粗糙度较小，其影响较弱，使得两者计算偏差有限。这间接证明了本文所提方法是可信的。

图3 考虑粗糙度的裂纹扩展模型验证

为兼顾计算精度与效率，尚需合理确定裂纹扩展区网格尺寸。为此，将裂纹扩展区网格尺寸分别设置为0.05、0.10、0.20、0.50mm，分别由Ms1、Ms2、Ms3、Ms4表示。利用所编制程序，计算得到不同网格尺寸下的应力强度因子如图4所示。为减小随机误差的影响，同一网格尺寸下随机重构5个粗糙度样本。以最小网格尺寸Ms1为标准，分析其他网格尺寸下应力强度因子的计算误差。相较而言，Ms4网格尺寸最大，其整体误差也较大，误差为1.68%～11.53%。Ms3与Ms2网格尺寸较小，计算误差也相对较小。其中，Ms3误差为0.15%～1.47%，Ms2误差为0%～0.35%。为获得较高的计算精度与效率，本文裂纹扩展区网格尺寸取为0.2mm。

粗糙度的引入使得该区域网格尺寸较小，导致计算效率大大降低。依据圣维南原理，在远离裂纹扩展的区域，粗糙度对裂纹扩展影响较小。因此，合理限制粗糙度的采样长度，保留更多的光滑齿面，既能保证足够的计算精度，也能获得较高的计算效率。为合理选取粗糙度的采样长度，将其分别设置为2.048、1.024、0.512mm，对应的采样点数分别为2048、1024、512。以裂纹为中心，在裂纹左右两侧添加粗糙度构建粗糙齿面。计算得到不同采样长度下Ⅰ型应力强度因子随裂纹长度的变化规律如图5所示。为减小随机误差的影响，同一采样长度下也随机重构5个粗糙度样本。以最大采样长度2.048mm为标准，分析其他采样长度下应力强度因子的计算误差。其中，采样点数1024对应的误差为0.2%～1.5%，采样点数512对应的误差为0.722%～3.368%。综合计算精度与效率，本文选取采样长度为1.024mm。

图4 网格质量验证

图5 粗糙度采样长度验证

为展示本文所提方法在齿根裂纹扩展计算方面的应用，对一薄轮缘齿轮，利用快速傅里叶变换(Fast Fourier Transform，FFT)重构方法随机重构粗糙度样本，分析轮廓均方根偏差Rq对其弯曲疲劳裂纹扩展的影响。模型与相关参数取值参看文献。在单齿啮合最高点施加载荷，大小取为5500N。为分析粗糙度对Mises应力的影响，分别对Rq=0.7和光滑齿面2种情况进行计算，齿根Mises应力对比结果见表1。表中列出了齿根区域节点中Mises应力值最大的前5个节点。由于齿根处粗糙度导致应力集中，Mises应力增大，最大应力从约228MPa(光滑齿面)增大到约365MPa(Rq=0.7)，增幅为60.09%。此外，相比于光滑齿面，Rq=0.7时Mises应力较大的节点空间分布更为分散，且主要出现在齿面粗糙度的波谷处，这些位置类似于缺口，产生显著的应力集中效应。

表1  齿根 Mises  应力值  单位:MPa

为进一步分析Rq对裂纹扩展路径与扩展寿命的影响，分别取Rq为0、0.1、0.3、0.5和0.7。其中，Rq=0代表光滑齿面。为减小随机误差的影响，同一Rq下随机重构5个粗糙度样本。计算得到裂纹扩展路径如图6所示。同一Rq下，不同样本所得路径分布在一定范围内。为便于清晰展示并从整体上分析Rq对裂纹扩展路径的影响，图6中只给出裂纹扩展路径的两边界值，分别用max(下边界)和min(上边界)表示，代表了裂纹扩展路径的分布范围。为讨论裂纹扩展路径的变化趋势，分析裂尖纵向位置随Rq的变化规律。以光滑齿面裂尖纵向位置为参考，Rq为0.1、0.3、0.5、0.7时，max裂尖位置分别降低了22.39、67.16、134.32、173.09μm，min裂尖位置分别降低了13.953、24.035、77.199、118.370μm。这表明，随着Rq增大，齿根裂纹扩展转向齿轮径向的趋势逐渐增大，导致轮缘断裂的概率随之增加。

图6 Rq对裂纹扩展路径影响

裂纹扩展寿命随Rq的变化如图7(a)所示。由图7(a)可见，裂纹扩展寿命随着Rq的增大而减小。相较于光滑齿面，轮廓均方根偏差Rq为0.1、0.3、0.5、0.7时，裂纹扩展寿命分别降低了1.02%、2.33%、6.67%、11.03%。由于齿根裂纹扩展类型以Ⅰ型为主，裂纹扩展寿命主要受KⅠ影响。为探究Rq对裂纹扩展寿命的影响机制，计算得到不同Rq下的Ⅰ型应力强度因子如图7(b)所示。计算表明，Rq为0.1、0.3、0.5、0.7时，KⅠ相较于光滑齿面结果分别增大了2.10%、4.49%、6.68%、8.03%。由此可见，Rq增大使得Ⅰ型应力强度因子增大，导致裂纹扩展寿命减小。

图7 Rq对裂纹扩展寿命与Ⅰ型应力强度因子影响

1)基于粗糙表面重构所得粗糙度高度序列，以最小拟合误差为原则进行拟合，通过空间坐标变换对齿根轮廓节点进行坐标调整，构建了齿根粗糙齿面几何模型，并基于扩展有限元法提出了考虑表面粗糙度的齿根裂纹扩展计算方法。通过对光滑齿面与较小粗糙度下的齿根裂纹扩展进行对比分析，间接验证了所提方法的有效性。

2)研究了网格尺寸和粗糙度采样长度对裂纹扩展计算的影响。为兼顾计算精度与效率，裂纹扩展区网格尺寸取0.2mm，此时计算误差为0.15%～1.47%；粗糙度采样长度取1.024mm，相应计算误差为0.2%～1.5%。

3)随着轮廓均方根偏差Rq的增大，齿根裂纹扩展转向齿轮径向的趋势逐渐增大，致使轮缘断裂的概率随之增加。相较于光滑齿面，Rq为0.1、0.3、0.5和0.7时，Ⅰ型应力强度因子KⅠ分别增大了2.10%、4.49%、6.68%、8.03%，导致裂纹扩展寿命分别降低了1.02%、2.33%、6.67%、11.03%。