CN107423462B - 工件考虑三维粗糙表面形貌的疲劳寿命预测方法及系统 - Google Patents

Abstract

本发明涉及工件疲劳寿命预测领域，公开一种工件考虑三维粗糙表面形貌的疲劳寿命预测方法及系统，以提高预测结果的准确性。本发明方法包括：获取工件危险区域的第一表面形貌高度序列矩阵；对第一表面形貌高度序列矩阵进行插值运算得到第二高度序列矩阵，并以三点共圆法计算插值后占主导地位波谷谷底的曲率半径，进而得到相应的危险区域的特征参数，然后根据该特征参数和基于第二高度序列矩阵构建的三维粗糙表面形貌有限元模型进行有限元计算求解得到表面形貌应力云图结果，对应力云图进行后处理得到应力集中系数和疲劳缺口系数，最后基于局部应力应变法对工件考虑三维粗糙表面形貌的疲劳寿命进行预测；提高了预测结果的准确性和精度。

Description

工件考虑三维粗糙表面形貌的疲劳寿命预测方法及系统

技术领域

本发明涉及工件疲劳寿命预测领域，尤其涉及一种工件考虑三维粗糙表面形貌的疲劳寿命预测方法及系统。

背景技术

工件的疲劳性能取决于表面完整性，直接引入表面系数考虑表面的影响是不精确的。表面形貌是表面完整性中的一个重要参数；表面形貌产生的应力集中效应对工件的疲劳寿命产生了重要影响，不同表面形貌对工件疲劳性能的影响差异非常大。现在大部分对疲劳寿命的预测是通过采取大量样本进行实验的方法，这花费了大量的时间和成本，而且因为样本与实际工件以及工作条件存在差异，导致预测与实验结果存在很大的不确定性。

长期以来，国内外学者对工件考虑表面形貌的疲劳寿命进行了研究。基于统计参数将模型简化为缺口模型，研究了表面形貌对应力集中效应和疲劳寿命的影响。统计参数只能描述表面形貌的平均特征，即使针对极值特征的参数Ry和Rz也会丢失表面形貌的许多信息。基于统计参数的简化模型难以反映实际表面形貌特征，这极大影响了预测结果的准确性。

发明内容

本发明目的在于公开一种工件考虑三维粗糙表面形貌的疲劳寿命预测方法及系统，以提高预测结果的准确性。

为实现上述目的，本发明公开了一种工件考虑三维粗糙表面形貌的疲劳寿命预测方法，包括：

步骤S1、通过宏观有限元分析确定工件的危险区域；

步骤S2、获取所述危险区域的第一表面形貌高度序列矩阵；

步骤S3、对所述第一表面形貌高度序列矩阵进行插值运算得到第二高度序列矩阵；

步骤S4、对所述第二高度序列矩阵内部各区域的波谷点按深度进行排序，将波谷视为球形凹坑，并以含波谷点z(i,j)的矩阵z(i-1:i+1,j-1:j+1)九个点所围成区域所对应对角线方向及中线方向分别构建出过波谷点的四个三角形，计算该四个三角形外接圆半径的平均值作为波谷的曲率半径；并以深度最深的波谷点所计算得出的曲率半径计算所述危险区域的特征参数，或以综合两个或两个以上深度最深的波谷点相对应曲率半径进行折中处理的结果计算所述危险区域的特征参数；

步骤S5、根据所述第二高度序列矩阵构建三维粗糙表面形貌有限元模型；

步骤S6、根据所述危险区域的特征参数和所述三维粗糙表面形貌有限元模型，进行有限元计算求解得到表面形貌应力云图结果，对应力云图进行后处理得到应力集中系数和疲劳缺口系数，然后基于局部应力应变法对工件考虑三维粗糙表面形貌的疲劳寿命进行预测；

其中，i、j为水平面X、Y轴向的分布点坐标，z为高度方向的坐标轴。

为实现上述目的，本发明还公开一种工件考虑三维粗糙表面形貌的疲劳寿命预测系统，包括：

第一处理模块、用于通过宏观有限元分析确定工件的危险区域；

第二处理模块、用于获取所述危险区域的第一表面形貌高度序列矩阵；

第三处理模块、用于对所述第一表面形貌高度序列矩阵进行插值运算得到第二高度序列矩阵；

第四处理模块、用于对所述第二高度序列矩阵内部各区域的波谷点按深度进行排序，将波谷视为球形凹坑，并以含波谷点z(i,j)的矩阵z(i-1:i+1,j-1:j+1)九个点所围成区域所对应对角线方向及中线方向分别构建出过波谷点的四个三角形，计算该四个三角形外接圆半径的平均值作为波谷的曲率半径；并以深度最深的波谷点所计算得出的曲率半径计算所述危险区域的特征参数，或以综合两个或两个以上深度最深的波谷点相对应曲率半径进行折中处理的结果计算所述危险区域的特征参数；

第五处理模块、用于根据所述第二高度序列矩阵构建三维粗糙表面形貌有限元模型；

第六处理模块、用于根据所述危险区域的特征参数和所述三维粗糙表面形貌有限元模型，进行有限元计算求解得到表面形貌应力云图结果，对应力云图进行后处理得到应力集中系数和疲劳缺口系数，然后基于局部应力应变法对工件考虑三维粗糙表面形貌的疲劳寿命进行预测；

其中，i、j为水平面X、Y轴向的分布点坐标，z为高度方向的坐标轴。

本发明具有以下有益效果：

对第一表面形貌高度序列矩阵进行插值运算得到第二高度序列矩阵，并以三点共圆的方法计算插值后占主导地位波谷的曲率半径，进而得到相应的危险区域的特征参数，然后根据该特征参数和基于第二高度序列矩阵构建的三维粗糙表面形貌有限元模型进行有限元计算求解得到表面形貌应力云图结果，对应力云图进行后处理得到应力集中系数和疲劳缺口系数，最后基于局部应力应变法对工件考虑三维粗糙表面形貌的疲劳寿命进行预测；其中，第二高度序列矩阵及以此生成的表面形貌有限元模型能精确反应工件微观层面的真实表面形貌，三点共圆的方法求解波谷曲率及整体方案都具有计算简便、依据合理科学等特点，从而提高了预测结果的准确性和精度。

下面将参照附图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是本发明优选实施例公开的工件考虑三维粗糙表面形貌的疲劳寿命预测方法流程图；

图2是本发明优选实施例公开的计算波谷点曲率时所构建三角形的分布示意图。

具体实施方式

以下结合附图对本发明的实施例进行详细说明，但是本发明可以由权利要求限定和覆盖的多种不同方式实施。

实施例1

本实施例公开一种工件考虑三维粗糙表面形貌的疲劳寿命预测方法，如图1所示，包括：

步骤S1、通过宏观有限元分析确定工件的危险区域。

步骤S2、获取危险区域的第一表面形貌高度序列矩阵。

该步骤可以采用白光干涉仪获取第一表面形貌高度序列矩阵，并将所获取的第一表面形貌高度序列矩阵记为：z(m,n)，其中，m、n分别为该水平面X、Y轴向的分布点数量。

步骤S3、对第一表面形貌高度序列矩阵进行插值运算得到第二高度序列矩阵。

可选的，该步骤采用贝塞尔曲线(Bezier Curve)插值。第一表面形貌高度序列矩阵z(m,n)差值后可得到第二高度序列矩阵为z((t+1)m-t,(t+1)n-t)，具体包括如下步骤：

步骤S31、将高度序列矩阵z(m,n)划分为(m-1)×(n-1)个区域,各区域分别对应水平面上相邻四个点所组成的四边形，并确定切矢k:

式中p为曲线上的n个离散数据点坐标矩阵，gradient为相邻点间的中心差分，norm为向量范数。

步骤S32、在曲线上相邻两个离散点之间确定两个控制点e,f：

式中α＝μ·norm(gradient(p))为控制系数；μ＝1/2。

步骤S33、将每个小区域的四个点两两连线得到共六条曲线，计算得到每条曲线上两个控制点，共得到十二个控制点，结合原来的四个测量点，构成控制矩阵V_4×4，以提供插值点的方向信息；并在每个小区域插值4t+t²个点，其中t∈{0,1,2,…}，利用控制矩阵V_4×4和Bezier曲面的矩阵形式计算每个小区域插值后的高度序列矩阵z(t+2,t+2)；计算公式如下：

其中，B(x)为Bernstein基函数，z(i,j)为小区域插值后第i行第j列元素，记录点的位置和高度信息。

步骤S4、对第二高度序列矩阵内部各区域的波谷点按深度进行排序，将波谷视为球形凹坑，并以波谷点为中心的矩阵九个点所围成区域所对应对角线方向及中线方向分别构建出过波谷点的四个三角形，计算该四个三角形外接圆半径的平均值作为波谷的曲率半径；并以深度最深的波谷点所计算得出的曲率半径计算危险区域的特征参数，或以综合两个或两个以上深度最深的波谷点相对应曲率半径进行折中处理的结果计算危险区域的特征参数(即：计算出占主导地位波谷谷底曲率半径的平均值)。优选的，一般取前五个最深的波谷作为占主导地位的波谷，并以计算占主导地位波谷得出的曲率半径的平均值计算所述危险区域的特征参数；通常所依据的波谷点越多，则越能反应粗糙表面的整体状况。

在该步骤中，若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最小值，则该点为波谷点。另一方面，以z(i,j)为波谷，则上述的几个点在矩阵中的分布为z(i-1:i+1,j-1:j+1)。

上述三角形的构建可以参照附图2，分别为过波谷的Δ1、Δ2、Δ3和Δ4四个三角形。在所构建的各个三角形中，因为各点的高度和间距已知，而在所构建三角形的平面内，通过勾股定理可求出三角形的三边长a,b,c，使用余弦定理和外接圆半径公式可求出半径R，然后计算四个三角形外接圆半径的平均值作为该波谷的曲率半径(为描述方便，下述简称该方法为“三点共圆法”)。在所构建的三角形ABC中，以z(i,j-1)所在的点为A，波谷点z(i,j)为点C，z(i,j+1)所在的点为B，则相关的计算公式如下：

在该步骤中，危险区域相关的特征参数包括敏感系数q，计算公式如下：

q＝1/(1+w/ρ)；

其中，w为材料系数，ρ为上述三点共圆法所确定的曲率半径。

步骤S5、根据第二高度序列矩阵构建三维粗糙表面形貌有限元模型。

可选的，该步骤的具体实现如下：

1)、根据第二高度序列矩阵z((t+1)m-t,(t+1)n-t)，生成包含3列的点文件，保存为CATIA(Computer Aided Three-dimensional Interactive Application，交互式CAD/CAE/CAM系统)能识别的asc格式。其中每一行的数据代表一个点的坐标，第1,2,3列分别表示对应点的xyz坐标。

2)、将点文件导入CATIA生成点云模型，通过曲面连接、曲面快速生成和拉伸功能构造粗糙表面形貌三维实体模型，保存为stp格式。

3)、将粗糙表面形貌三维实体模型stp格式文件导入ABAQUS(一套功能强大的工程模拟的有限元软件)中，对模型进行四面体网格划分，通过取决于表面曲率的单元大小来控制网格密度，指定网格单元类型与材料属性；在表面形貌的两侧施加一个微小拉伸应力σ，底部添加固定约束，构建出三维粗糙表面形貌有限元模型。

步骤S6、根据危险区域的特征参数和三维粗糙表面形貌有限元模型，进行有限元计算求解得到表面形貌应力云图结果，对应力云图进行后处理得到应力集中系数和疲劳缺口系数，然后基于局部应力应变法对工件考虑三维粗糙表面形貌的疲劳寿命进行预测。

可选的，该步骤具体包括：

步骤S61、根据有限元计算求解得到表面形貌应力云图结果，提取占主导地位波谷的应力大小，计算其平均值，作为最大应力σ_max。

步骤S62、将施加在步骤S5中构建的三维粗糙表面形貌有限元模型两侧的应力σ作为名义应力σ₀，采用Kt＝σ_max/σ₀计算应力集中系数。

步骤S63、根据应力集中系数与疲劳缺口系数之间的关系，K_f＝1+q(K_t-1)，计算疲劳缺口系数K_f，q为敏感系数。

步骤S64、采用线性理论ε₂＝K_fε₁计算考虑粗糙表面形貌引起应力集中时的真实应变ε₂，其中ε₁为光滑情况下应变。

步骤S65、根据材料疲劳特性参数，基于局部应力应变法将参数ε₂和ε₁分别代入曼森—科芬方程中的ε进行求解，分别得到考虑粗糙表面形貌和光滑状态下的疲劳寿命N_f；计算公式如下：

其中：ε为总应变，ε_e为弹性应变，ε_p为塑性应变，σ'_f为疲劳强度系数，b为疲劳强度指数，N_f以循环计数的疲劳寿命，ε'_f为疲劳延性系数，c为疲劳延性指数，E为弹性模量。

【应用例一】

以齿轮工件宏观有限元分析为例，17CrNiMo6样本三维粗糙表面形貌为例，对三维粗糙表面形貌测量和参数计算进行说明，结合CATIA和ABAQUS软件，依托上述方法进行相关数据处理。

A：以齿轮工件为对象，齿轮进行六面体网格划分，受力分析找到齿根最大弯曲应力位置，该位置为危险区域，容易出现疲劳破坏。

B：采用白光干涉仪测量17CrNiMo6样本的表面形貌，采样间隔为1μm，得到100×100点数的粗糙表面形貌，生成高度序列矩阵z(100,100)。17CrNiMo6参数如表1所示。

表1：

C：将第一表面形貌高度序列矩阵z(100,100)划分为99×99个区域,X/Y方向相邻四点为一个区域，每个小区域插值12个点；构成第二高度序列矩阵z(298,298)，并计算出占主导地位波谷谷底曲率半径的平均值ρ为3.0451μm，采用Peterson公式计算出敏感系数q为0.0468。

D：在根据第二高度序列矩阵构建三维粗糙表面形貌有限元模型的过程中，在表面形貌的两侧施加一个微小拉伸应力σ＝100MPA。

E：根据有限元计算求解得到表面形貌应力云图结果，提取占主导地位波谷的应力大小，计算出平均值为218MPa，作为最大应力σ_max。

F：将施加在模型两侧的应力σ＝100MPa作为名义应力σ₀，采用公式计算应力集中系数为2.18。

G：根据应力集中系数与疲劳缺口系数之间的关系计算出疲劳缺口系数为：

K_f＝1+q(K_t-1)＝1.0552。

H：采用线性理论ε₂＝1.0552ε₁计算考虑粗糙表面形貌引起应力集中时的真实应变ε₂，其中ε₁为光滑情况下应变；然后根据材料疲劳特性参数，基于局部应力应变法将参数ε₂和ε₁分别代入曼森—科芬方程中的ε进行求解，分别得到考虑粗糙表面形貌和光滑状态下的疲劳寿命N_f如表2所示。进一步的，计算几组不同应变下的寿命，可以作出相应的疲劳寿命曲线。

表2：

综上，本实施例公开的工件考虑三维粗糙表面形貌的疲劳寿命预测方法，对第一表面形貌高度序列矩阵进行插值运算得到第二高度序列矩阵，并以三点共圆的方法计算插值后占主导地位波谷的曲率半径，进而得到相应的危险区域的特征参数，然后根据该特征参数和基于第二高度序列矩阵构建的三维粗糙表面形貌有限元模型进行有限元计算求解得到表面形貌应力云图结果，对应力云图进行后处理得到应力集中系数和疲劳缺口系数，最后基于局部应力应变法对工件考虑三维粗糙表面形貌的疲劳寿命进行预测；其中，第二高度序列矩阵及以此生成的表面形貌有限元模型能精确反应工件微观层面的真实表面形貌，三点共圆的方法求解波谷曲率及整体方案都具有计算简便、依据合理科学等特点，从而提高了预测结果的准确性和精度。

实施例2

与上述方法实施例相对应的，本实施例公开一种工件考虑三维粗糙表面形貌的疲劳寿命预测系统，包括：

第一处理模块、用于通过宏观有限元分析确定工件的危险区域；

第二处理模块、用于获取危险区域的第一表面形貌高度序列矩阵；

第三处理模块、用于对第一表面形貌高度序列矩阵进行插值运算得到第二高度序列矩阵；

第四处理模块、用于对第二高度序列矩阵内部各区域的波谷点按深度进行排序，将波谷视为球形凹坑，并以含波谷点z(i,j)的矩阵z(i-1:i+1,j-1:j+1)九个点所围成区域所对应对角线方向及中线方向分别构建出过波谷点的四个三角形，计算该四个三角形外接圆半径的平均值作为波谷的曲率半径；并以深度最深的波谷点所计算得出的曲率半径计算危险区域的特征参数，或以综合两个或两个以上深度最深的波谷点相对应曲率半径进行折中处理的结果计算危险区域的特征参数；

第五处理模块、用于根据第二高度序列矩阵构建三维粗糙表面形貌有限元模型；

第六处理模块、用于根据危险区域的特征参数和三维粗糙表面形貌有限元模型，进行有限元计算求解得到表面形貌应力云图结果，对应力云图进行后处理得到应力集中系数和疲劳缺口系数，然后基于局部应力应变法对工件考虑三维粗糙表面形貌的疲劳寿命进行预测。

可选的，上述第二处理模块通过白光干涉仪获取危险区域的第一表面形貌高度序列矩阵。以第一表面形貌高度序列矩阵为z(m,n)，第二高度序列矩阵为z((t+1)m-t,(t+1)n-t)为例进行说明；可选的，第三处理模块包括：

子模块一、用于将高度序列矩阵z(m,n)划分为(m-1)×(n-1)个区域,各区域分别对应水平面上相邻四个点所组成的四边形，m、n分别为该水平面X、Y轴向的分布点数量；

确定切矢k:

式中p为曲线上的n个离散数据点坐标矩阵，gradient为相邻点间的中心差分，norm为向量范数；

子模块二、用于在曲线上相邻两个离散点之间确定两个控制点e,f：

式中α＝μ·norm(gradient(p))为控制系数；μ＝1/2；

子模块三、用于将每个小区域的四个点两两连线得到共六条曲线，计算得到每条曲线上两个控制点，共得到十二个控制点，结合原来的四个测量点，构成控制矩阵V_4×4，以提供插值点的方向信息；并在每个小区域插值4t+t²个点，其中t∈{0,1,2,…}，利用控制矩阵V_4×4和Bezier曲面的矩阵形式计算每个小区域插值后的高度序列矩阵z(t+2,t+2)；计算公式如下：

其中，B(x)为Bernstein基函数，z(i,j)为小区域插值后第i行第j列元素，记录点的位置和高度信息。

可选的，上述第四处理模块所计算的特征参数包括敏感系数q，计算公式如下：

q＝1/(1+w/ρ)；

其中，w为材料系数，ρ为步骤S4所确定的曲率半径。

可选的，上述第六处理模块包括：

第一单元、用于根据有限元计算求解得到表面形貌应力云图结果，提取占主导地位波谷的应力大小，计算其平均值，作为最大应力σ_max；

第二单元、用于将施加在第五处理模块构建的三维粗糙表面形貌有限元模型两侧的应力σ作为名义应力σ₀，采用Kt＝σ_max/σ₀计算应力集中系数；

第三单元、用于根据应力集中系数与疲劳缺口系数之间的关系，K_f＝1+q(K_t-1)，计算疲劳缺口系数K_f，q为敏感系数；

第四单元、用于采用线性理论ε₂＝K_fε₁计算考虑粗糙表面形貌引起应力集中时的真实应变ε₂，其中ε₁为光滑情况下应变；

第五单元、用于根据材料疲劳特性参数，基于局部应力应变法将参数ε₂和ε₁分别代入曼森—科芬方程中的ε进行求解，分别得到考虑粗糙表面形貌和光滑状态下的疲劳寿命N_f；计算公式如下：

其中：ε为总应变，ε_e为弹性应变，ε_p为塑性应变，σ'_f为疲劳强度系数，b为疲劳强度指数，N_f以循环计数的疲劳寿命，ε'_f为疲劳延性系数，c为疲劳延性指数，E为弹性模量。

同理，本实施例公开的工件考虑三维粗糙表面形貌的疲劳寿命预测系统，对第一表面形貌高度序列矩阵进行插值运算得到第二高度序列矩阵，并以三点共圆的方法计算插值后占主导地位波谷的曲率半径，进而得到相应的危险区域的特征参数，然后根据该特征参数和基于第二高度序列矩阵构建的三维粗糙表面形貌有限元模型进行有限元计算求解得到表面形貌应力云图结果，对应力云图进行后处理得到应力集中系数和疲劳缺口系数，最后基于局部应力应变法对工件考虑三维粗糙表面形貌的疲劳寿命进行预测；其中，第二高度序列矩阵及以此生成的表面形貌有限元模型能精确反应工件微观层面的真实表面形貌，三点共圆的方法求解波谷曲率及整体方案都具有计算简便、依据合理科学等特点，从而提高了预测结果的准确性和精度。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种工件考虑三维粗糙表面形貌的疲劳寿命预测方法，其特征在于，包括：

步骤S1、通过宏观有限元分析确定工件的危险区域；

步骤S2、获取所述危险区域的第一表面形貌高度序列矩阵；

步骤S3、对所述第一表面形貌高度序列矩阵进行插值运算得到第二高度序列矩阵；

步骤S4、对所述第二高度序列矩阵内部各区域的波谷点按深度进行排序，将波谷视为球形凹坑，并以含波谷点z(i,j)的矩阵z(i-1:i+1,j-1:j+1)九个点所围成区域所对应对角线方向及中线方向分别构建出过波谷点的四个三角形，计算该四个三角形外接圆半径的平均值作为波谷的曲率半径；并以深度最深的波谷点所计算得出的曲率半径计算所述危险区域的特征参数，或以综合两个或两个以上深度最深的波谷点相对应曲率半径进行折中处理的结果计算所述危险区域的特征参数；

步骤S5、根据所述第二高度序列矩阵构建三维粗糙表面形貌有限元模型；

步骤S6、根据所述危险区域的特征参数和所述三维粗糙表面形貌有限元模型，进行有限元计算求解得到表面形貌应力云图结果，对应力云图进行后处理得到应力集中系数和疲劳缺口系数，然后基于局部应力应变法对工件考虑三维粗糙表面形貌的疲劳寿命进行预测；

其中，i、j为水平面X、Y轴向的分布点坐标，z为高度方向的坐标轴。

2.根据权利要求1所述的工件考虑三维粗糙表面形貌的疲劳寿命预测方法，其特征在于，所述第一表面形貌高度序列矩阵为z(m,n)，所述第二高度序列矩阵为z((t+1)m-t,(t+1)n-t)；所述步骤S3包括：

步骤S31、将高度序列矩阵z(m,n)划分为(m-1)×(n-1)个区域,各区域分别对应水平面上相邻四个点所组成的四边形，m、n分别为该水平面X、Y轴向的分布点数量；

确定切矢k:

式中p为曲线上的n个离散数据点坐标矩阵，gradient为相邻点间的中心差分，norm为向量范数；

步骤S32、在曲线上相邻两个离散点之间确定两个控制点e,f：

式中α＝μ·norm(gradient(p))为控制系数；μ＝1/2；

步骤S33、将每个小区域的四个点两两连线得到共六条曲线，计算得到每条曲线上两个控制点，共得到十二个控制点，结合原来的四个测量点，构成控制矩阵V_4×4，以提供插值点的方向信息；并在每个小区域插值4t+t²个点，其中t∈{0,1,2,…}，利用控制矩阵V_4×4和Bezier曲面的矩阵形式计算每个小区域插值后的高度序列矩阵z(t+2,t+2)；计算公式如下：

其中，B(x)为Bernstein基函数，z(i,j)为小区域插值后第i行第j列元素，记录点的位置和高度信息。

3.根据权利要求1或2所述的工件考虑三维粗糙表面形貌的疲劳寿命预测方法，其特征在于，所述步骤S4中的特征参数包括敏感系数q，计算公式如下：

q＝1/(1+w/ρ)；

其中，w为材料系数，ρ为步骤S4所确定的曲率半径。

4.根据权利要求3所述的工件考虑三维粗糙表面形貌的疲劳寿命预测方法，其特征在于，所述步骤S6包括：

步骤S61、根据有限元计算求解得到表面形貌应力云图结果，提取占主导地位波谷的应力大小，计算其平均值，作为最大应力σ_max；

步骤S62、将施加在步骤S5中构建的三维粗糙表面形貌有限元模型两侧的应力σ作为名义应力σ₀，采用Kt＝σ_max/σ₀计算应力集中系数；

步骤S63、根据应力集中系数与疲劳缺口系数之间的关系，K_f＝1+q(K_t-1)，计算疲劳缺口系数K_f，q为敏感系数；

步骤S64、采用线性理论ε₂＝K_fε₁计算考虑粗糙表面形貌引起应力集中时的真实应变ε₂，其中ε₁为光滑情况下应变；

步骤S65、根据材料疲劳特性参数，基于局部应力应变法将参数ε₂和ε₁分别代入曼森—科芬方程中的ε进行求解，分别得到考虑粗糙表面形貌和光滑状态下的疲劳寿命N_f；计算公式如下：

其中：ε为总应变，ε_e为弹性应变，ε_p为塑性应变，σ'_f为疲劳强度系数，b为疲劳强度指数，N_f以循环计数的疲劳寿命，ε'_f为疲劳延性系数，c为疲劳延性指数，E为弹性模量。

5.根据权利要求4所述的工件考虑三维粗糙表面形貌的疲劳寿命预测方法，其特征在于，通过白光干涉仪获取所述危险区域的第一表面形貌高度序列矩阵。

6.一种工件考虑三维粗糙表面形貌的疲劳寿命预测系统，其特征在于，包括：

第一处理模块、用于通过宏观有限元分析确定工件的危险区域；

第二处理模块、用于获取所述危险区域的第一表面形貌高度序列矩阵；

第三处理模块、用于对所述第一表面形貌高度序列矩阵进行插值运算得到第二高度序列矩阵；

第四处理模块、用于对所述第二高度序列矩阵内部各区域的波谷点按深度进行排序，将波谷视为球形凹坑，并以含波谷点z(i,j)的矩阵z(i-1:i+1,j-1:j+1)九个点所围成区域所对应对角线方向及中线方向分别构建出过波谷点的四个三角形，计算该四个三角形外接圆半径的平均值作为波谷的曲率半径；并以深度最深的波谷点所计算得出的曲率半径计算所述危险区域的特征参数，或以综合两个或两个以上深度最深的波谷点相对应曲率半径进行折中处理的结果计算所述危险区域的特征参数；

第五处理模块、用于根据所述第二高度序列矩阵构建三维粗糙表面形貌有限元模型；

第六处理模块、用于根据所述危险区域的特征参数和所述三维粗糙表面形貌有限元模型，进行有限元计算求解得到表面形貌应力云图结果，对应力云图进行后处理得到应力集中系数和疲劳缺口系数，然后基于局部应力应变法对工件考虑三维粗糙表面形貌的疲劳寿命进行预测；

其中，i、j为水平面X、Y轴向的分布点坐标，z为高度方向的坐标轴。

7.根据权利要求6所述的工件考虑三维粗糙表面形貌的疲劳寿命预测系统，其特征在于，所述第一表面形貌高度序列矩阵为z(m,n)，所述第二高度序列矩阵为z((t+1)m-t,(t+1)n-t)；所述第三处理模块包括：

子模块一、用于将高度序列矩阵z(m,n)划分为(m-1)×(n-1)个区域,各区域分别对应水平面上相邻四个点所组成的四边形，m、n分别为该水平面X、Y轴向的分布点数量；

确定切矢k:

式中p为曲线上的n个离散数据点坐标矩阵，gradient为相邻点间的中心差分，norm为向量范数；

子模块二、用于在曲线上相邻两个离散点之间确定两个控制点e,f：

式中α＝μ·norm(gradient(p))为控制系数；μ＝1/2；

子模块三、用于将每个小区域的四个点两两连线得到共六条曲线，计算得到每条曲线上两个控制点，共得到十二个控制点，结合原来的四个测量点，构成控制矩阵V_4×4，以提供插值点的方向信息；并在每个小区域插值4t+t²个点，其中t∈{0,1,2,…}，利用控制矩阵V_4×4和Bezier曲面的矩阵形式计算每个小区域插值后的高度序列矩阵z(t+2,t+2)；计算公式如下：

其中，B(x)为Bernstein基函数，z(i,j)为小区域插值后第i行第j列元素，记录点的位置和高度信息。

8.根据权利要求6或7所述的工件考虑三维粗糙表面形貌的疲劳寿命预测系统，其特征在于，所述第四处理模块所计算的特征参数包括敏感系数q，计算公式如下：

q＝1/(1+w/ρ)；

其中，w为材料系数，ρ为步骤S4所确定的曲率半径。

9.根据权利要求8所述的工件考虑三维粗糙表面形貌的疲劳寿命预测系统，其特征在于，所述第六处理模块包括：

第一单元、用于根据有限元计算求解得到表面形貌应力云图结果，提取占主导地位波谷的应力大小，计算其平均值，作为最大应力σ_max；

第二单元、用于将施加在第五处理模块构建的三维粗糙表面形貌有限元模型两侧的应力σ作为名义应力σ₀，采用Kt＝σ_max/σ₀计算应力集中系数；

第三单元、用于根据应力集中系数与疲劳缺口系数之间的关系，K_f＝1+q(K_t-1)，计算疲劳缺口系数K_f，q为敏感系数；

第四单元、用于采用线性理论ε₂＝K_fε₁计算考虑粗糙表面形貌引起应力集中时的真实应变ε₂，其中ε₁为光滑情况下应变；

第五单元、用于根据材料疲劳特性参数，基于局部应力应变法将参数ε₂和ε₁分别代入曼森—科芬方程中的ε进行求解，分别得到考虑粗糙表面形貌和光滑状态下的疲劳寿命N_f；计算公式如下：

其中：ε为总应变，ε_e为弹性应变，ε_p为塑性应变，σ'_f为疲劳强度系数，b为疲劳强度指数，N_f以循环计数的疲劳寿命，ε'_f为疲劳延性系数，c为疲劳延性指数，E为弹性模量。

10.根据权利要求9所述的工件考虑三维粗糙表面形貌的疲劳寿命预测系统，其特征在于，所述第二处理模块通过白光干涉仪获取所述危险区域的第一表面形貌高度序列矩阵。