基于傅里叶变换的加工表面应力集中系数建模*

黄 烁，宋俨轩，殷国富

(四川大学机械工程学院，成都 610065)

0 引言

机械加工会对零件的表面形貌造成损伤，使得零件局部高度应力集中，加速了零件疲劳裂纹的萌生，使得零件疲劳强度下降、寿命减短[1]。为了预测及提高航空零件的疲劳寿命，学者们开始研究零件表面几何形貌与寿命之间的联系[2]。据前苏联进行的大量表面完整性研究显示，零件的表面粗糙度对疲劳强度的影响最大，达到了50%～60%，因此降低零件表面粗糙度是提高疲劳强度的重要措施[3]。

表面粗糙度Ra表征了工件的形貌特征。研究表明在无残余应力的情况下，Ra超过0.1 μm对疲劳寿命有显著影响[4]。但随着各种计算需要，表征零件表面形貌特征的参数越来越多，考虑的因素愈渐复杂。学者们认为单一的Ra没有考虑到疲劳损伤发生的物理诱因，即由于加工表面不光滑引起的应力集中。因此，AROLA、NEUBER等[5-7]提出了应力集中系数Kt的概念，将轮廓谷底曲率半径ρ、微观不平度十点高度Rz等因素考虑其中，把表面几何特征与局部应力状态联系在一起，更为客观地表征了表面形貌对于疲劳寿命的影响。

但在使用Arola模型预测应力集中系数时面临两个问题：

一是测量轮廓谷底曲率半径ρ时一般采用半径规比对[6]，这样的测量手段具有极大的主观测量误差与随机性误差[8]，无法真正应用在加工质量评定中。为了解决第一个问题，LI、XIONG、AS等[9-11]采用了贝塞尔插值函数对轮廓曲线拟合重构,通过计算得到ρ，将整个参数识别过程自动化。

二是在实际的测量过程中，由于外界振动与设备精度问题使得测得的曲线夹杂了噪声，这些高频成分的噪声对于工件的疲劳寿命没有影响[4]，但对Arola模型的计算影响非常大。而解决第二个问题需要对所得的轮廓曲线进行滤波去噪处理。

本文基于Arola模型，采用傅里叶变换的方法对表面轮廓曲线进行重构和滤波去噪两步处理，使得Kt的预测结果更准确。本文通过对构造曲面的数值试验和疲劳实验两种方法来验证该算法的有效性和优越性。在疲劳试验中观察Kt与疲劳寿命的关系验证方法预测Kt值的有效性。在对构造曲面的数值试验中，通过和构造曲面的有限元仿真值Kt的对比，体现出本文方法对Kt值预测的准确性。

1 表面应力集中系数建模方法

1.1 应力集中系数

工件的表面是由无数个波峰波谷所组成的，而在轮廓波谷处会产生应力集中的现象。用式(1)来定义应力集中系数Kt：

(1)

式中，σmax为波谷处所受最大应力；σnom为工件横截面处所受的名义应力。

但Kt的定义式实际应用较为困难。同时为了描述加工表面与应力集中之间的关系，许多学者提出了用部分表面参数计算表征Kt。文献[5]提出了一个半经验公式：

(2)

式中，Ry为轮廓峰谷总高度；n为不同的应力状态，当工件受到剪切应力时，n=1；当工件受到了拉应力或者弯曲应力时，n=2。轮廓谷底曲率半径ρ的计算如式(3)所示。

(3)

将采样曲线段中所有纵坐标值小于左右相邻两处值的点记为“波谷”，ρ取最小值记为轮廓谷底最小曲率半径ρmin。

1.2 表面轮廓曲线重建

在使用Alora等模型进行预测时，测量轮廓谷底曲率半径ρ时一般采用半径规比对，这样的测量手段具有极大的主观测量误差与随机性误差，无法真正应用在加工质量评定中。另一方面在使用轮廓仪对工件表面进行采样时，通常采集到的轮廓曲线为由数个点相连接的折线。在预测Kt时，会导致Kt计算所需的ρ过小，产生奇点，从而导致计算结果误差较大。因此，在计算Kt值之前，需要将轮廓表面地形重建，从而更加准确且客观地测量计算出ρ与Kt值。

重构的方法有很多，文献[9-10]采用了样条插值函数的方法进行重建。采用一组三阶贝塞尔样条插值公式进行表面地形重建，使得轮廓曲线更加平滑。

本文采用傅里叶变换的方法来重建表面轮廓。傅里叶变换法的原理是将表面轮廓曲线视为时域信号函数，通过傅里叶变换，将时域信号函数转换为频域信号函数，根据频域函数图的表征的幅频特性分析，提取出各个频率信号的幅值、相位等，可以将原信号函数表示为多个三角函数的叠加。

(4)

(5)

Pi=atan2(b,a)

(6)

(7)

式中，C为重构后的轮廓曲线函数；Ai为第i阶谐波的幅值；a和b分别为原函数点经过傅里叶变换得到的复数的实部与虚部；N为采样点数；fi为第i阶谐波的频率；fs为采样频率，即时域图中相邻两点的距离(采样时间)的倒数；Pi为第i阶谐波的相位。

1.3 表面轮廓曲线的滤波去噪

(8)

(9)

式中，HV为材料的维氏硬度，查阅文献[14]得到TC11材料洛氏硬度HRS=27，根据硬度换算表查得维氏硬度HV=280；cmin为最小缺陷深度，当轮廓曲线在cmin的范围内波动，对疲劳极限不造成影响。所以可对测得的表面轮廓曲线进行滤波去噪处理，优化过的曲线方程R(x)与原始曲线方程W(x)需要满足下式：

|R(x)-W(x)|=c

(10)

滤波去噪的方法很多，最主要分为时域滤波方法与频域滤波方法。时域滤波方法通常采用五点滑动平均法，是利用最小二乘法原理对离散数据进行线性平滑的直线滑动平均法，可以削弱干扰信号的影响，提高振动曲线光滑度。

本文采用的是频域滤波方法——傅里叶变换法。对采样得到的曲线进行离散傅里叶变换，将轮廓曲线函数转换为频域函数。根据频域函数图，将频幅过小的频率点清零，把剩余的频率点代入式(7)，还原成时域函数。频域方法具有较好的频率选择性和灵活性，且运算的速度要高于时域方法。

2 实验验证

提供了一种考虑测量噪声条件下Kt的计算方法。采用了两种方法来验证该算法的优越性：对构造曲面的数值试验和疲劳实验。由于粗糙表面的Kt值越大，在相同载荷条件下的低周疲劳寿命越短。因此，Kt值是否与疲劳寿命呈反比是实验验证计算方法的一个准则。然而，仅有这种尺度的实验验证难以定量说明该算法在计算Kt方面的优越性，因为由于环境与设备因素的影响，“真实”的表面轮廓曲线是无法获知的。

因此，采取人为设计构造真实曲面并随机加入白噪声的方法来获取两种样本模型：真实曲面a和含噪声曲面b。其中真实曲面a代表了工件表面轮廓的实际情况，含噪声曲面b代表由于测量环境和设备等误差因素影响后得到的所测情况。

在验证轮廓曲线重构方法的时候，将真实曲面a的Kt有限元仿真值、傅里叶变换法和其他优化方法计算的真实曲面a的Kt值进行对比，证明本文算法在自动计算ρ值时的准确性。

在验证轮廓曲线滤波去噪方法的时候，将真实曲面a的Kt有限元仿真值，与原始Arola算法、傅里叶变换算法以及现有的其他优化算法三者计算的含噪声曲面b的Kt值作对比来验证本文算法在计算含测量误差表面的Kt值时的优越性。

2.1 轮廓曲线重构方法验证

在探究表面形貌对于Kt值的影响时，为了得到模型的Kt值，利用ABAQUS软件对目标模型进行有限元分析，将分析结果值作为该模型的Kt标准值与模型预测值作对比。查询文献[14]得到材料TC11的杨氏模量为123 GPa，泊松比为0.33。计算结果如图1所示。

图1 ABAQUS有限元仿真计算结果图

为了定量分析方法的有效性，本文设计了两组标准波形，如图2所示，波形参数如表1所示。

(a) 三角波 (b) 正弦波图2 波形示意图

表1 三角波与正弦波实验参数设计表

取表1中两条标准曲线来模拟工件表面轮廓地形，分别使用傅里叶变换与贝塞尔插值函数进行轮廓曲线重构。表2和表3为原始曲线有限元计算结果，与两种仿真方法对原始轮廓曲线重建后的各项表面形貌特征量的对比。

(11)

表2 三角波有限元仿真、傅里叶变换与插值函数重构表面参数对比

表3 正弦波有限元仿真、傅里叶变换与插值函数重构表面参数对比

可以看出，傅里叶变换方法重建的表面轮廓模型预测出的结果ρ和Kt与有限元仿真值对比，相对误差小于10%，较为准确。θ值越小，表明预测值与仿真值越接近。傅里叶变换法预测值与贝塞尔插值函数预测值相比，更为接近有限元仿真结果。

2.2 轮廓曲线滤波去噪方法验证

为了模拟实际测量曲线的情况，在表1的标准波形中加入高斯白噪声。噪声的大小由信噪比SNR衡量。由于在测量的过程中，噪声来源于仪器设备和测试环境的影响，因此使用接触式轮廓仪对仪器校准用的标准模版(Ra=1.6)进行测量。SNR的测量计算公式可以表达为：

M=2k

(12)

(13)

(14)

(15)

式中，M为量化级数，测量模版得到的曲线波峰最大值为2.17 μm，波谷最小值为-2.36 μm，精度为0.01 μm，因此取量化级数M=2k=29=512。

SNR=3×22k-1=3×217=393 216

(16)

SNRdB=1.76+6.02k=55.94 dB

(17)

因此取信噪比为56 dB的高斯白噪声加入正弦信号波，经由MATLAB计算加入噪声后的信号W(x)与优化过后的信号R(x)之间的最大差值cmax=0.128 nm，小于接触式轮廓仪说明书中的示值误差5 nm，因此可以得出56 dB的噪声属于仪器造成的误差。

傅里叶变换法需要根据实际的要求，将需要滤除的频率部分清零。由于过滤清零的比例无法确定，过滤得太少无法减少高频部分对曲率半径ρ的影响，过滤得太多会导致优化后的曲线与原始曲线偏差过大。为了验证过滤比例多少最为合适，本文取表1中曲线加入白噪声，以原曲线的有限元仿真的Kt为标准值，与占频幅最大值不同比例的过滤预测值进行对比，如表4和表5所示。

表4 傅里叶变换法的正弦波预测值与仿真值对比

续表

表5 傅里叶变换法的三角波预测值与仿真值对比

可以看出，当过滤去噪的比例达到0.5%～0.55%的时候，θ的值最小，模型预测出来的结果最接近原始曲线有限元仿真值，因此过滤去噪比例0.525%为最佳值。

为了横向对比傅里叶变换法滤波去噪的效果，采用了五点滑动平均法做了一组实验。利用五点滑动平均法处理曲线时，曲线的光滑程度与平滑处理的次数有关。与傅里叶变换法相同，需要对处理次数进行探讨。图3和图4为两种方法不同处理次数的预测值与有限元仿真值的对比。

图3 正弦波傅里叶变换法与五点滑动平均法优化后值对比图 图4 三角波傅里叶变换法与五点滑动平均法优化后值对比图

由图可以得出当平滑次数达到50～60次左右时，模型预测值最为接近有限元仿真值。但傅里叶变换法的θ值要小于五点滑动平均法。因此使用傅里叶变换法优于五点滑动平滑法，且过滤比例为0.525%最为准确。

2.3 疲劳实验验证

根据Kt值越大，疲劳寿命越短的原理设计疲劳实验，验证本文方法在实际应用之中的可行性和优越性,实验流程如图5所示。

图5 实验流程图

本文采用的TC11棒材如图6所示，共6根。

图6 TC11试件尺寸

实验采用接触式轮廓仪对疲劳试样表面截面的二维轮廓进行测量。由于在疲劳试验过程中，断裂区大多出现在试样的中间位置。因此需要对该区域的二维剖面进行测量和分析。设定采样间隔为1 μm，取样长度为20 mm，共计采样20 000点，尽可能重建零件表面形貌如图7所示。

每根棒材上随机采样3条轮廓曲线，计算Kt值时取平均值。采用PQ-60四点旋转弯曲疲劳试验机进行疲劳试验，加载应力为550 MPa，转速3000 r/min。将未经处理的Alora公式预测结果与傅里叶变换后的预测结果进行对比，如图8所示。

图7 试件表面轮廓部分测量采集点与重构滤波去噪后的曲线对比 图8 傅里叶变换与原Alora公式预测的Kt值和疲劳寿命之间关系的对比图

可以看出，原Alora公式预测的Kt值随着试件疲劳寿命的增加而波动较大。而经过傅里叶变换滤波去噪后Kt值随着试件疲劳寿命的增大而减小，符合事实规律，该方法具有可行性和优越性。

3 结论

(1)提出了利用傅里叶变换重构表面轮廓曲线的方法，以标准正弦波和三角波为模型，以ABAQUS有限元仿真结果作为标准值，验证了该方法的有效性。通过横向对比了三阶贝塞尔样条插值函数对轮廓表面形貌重构的结果，得出傅里叶变换重构出来的表面形貌更为贴近原始函数，预测的Kt更接近仿真值。

(2)提出了利用傅里叶变换法对表面轮廓曲线滤波去噪的方法。通过测量计算得到，人为加入信噪比为56 dB的高斯白噪声的正弦波最接近真实工况。以加入了噪声的标准波形为研究对象，验证了该方法的有效性。通过与有限元仿真值的比较得出，滤波比例为0.525%的时候最为准确。同时以五点滑动平均法为对照组，对比分析得出傅里叶变换法优于五点滑动平均法。

(3)采集6个相同加工参数的试件的表面轮廓数据，并进行疲劳实验。结果得出随着疲劳寿命增大，经过傅里叶变换预测得到的Kt值越小，而原Arola公式预测的Kt值不稳定。证明了本文方法的有效性。