基于分形理论和蒙特卡罗方法的粗糙表面形貌模拟方法

技术领域

本发明属于机加工表面微观形貌表征与模拟研究领域，具体是一种基于分形理论和蒙特卡罗方法的粗糙表面形貌模拟方法。

背景技术

广泛存在于各种工程领域中的接触现象，是摩擦学研究的重要课题之一。从微观角度来看，接触表面中的粗糙度、加工纹理和波纹度等复杂几何形貌，显著影响接触零件的力学性能和物理性能。然而，表面形貌研究多数基于大量的实验数据，实验成本高且表面形貌难以重复验证。为了建立微观形貌与界面性能的定量映射规律并简化实验成本，有必要寻求一种三维表面形貌的精确描述和模拟构建方法。基于该方法精确模拟机械元件摩擦学特性，可在节省实验研究时间和成本的同时，得到在特定摩擦学条件下获得最佳表面形貌。同时有助于分析零件接触承载能力、降低磨损以及改善摩擦状况。

对于表面微观形貌数值模拟的研究，一方面，部分研究学者从统计学的角度入手，从轮廓幅度、纹理走向、高度分布等多角度评价表面形貌。然而，许多参数并非表面形貌的固有参数，而是依赖于测量仪器的探测精度。为了更好的表征形貌结构，不断引入新的参数成为了唯一的方法。“参数爆炸”的后果是难以定量和客观的确定最相关的形貌参数。另一方面的研究是，基于自相似性和尺度独立性的分形理论表面模拟方法。该方法被运用于接触刚度、摩擦粘滑、电导率和热导率等场合中。然而，部分问题仍然存在。例如，分形参数由于数量有限，只能反映部分特征信息而不能完全表征粗糙表面形貌。同时，尺度依赖性这一基本困难仍未得到彻底解决。

针对现存的表面模拟方法的不足，很有必要设计一种将分形理论和传统统计学参数相结合的表面形貌模拟方法来准确地描述加工表面微观形貌特征。

发明内容

针对现有技术的不足，本发明拟解决的技术问题是，提供一种基于分形理论和蒙特卡罗方法的粗糙表面形貌模拟方法。

本发明解决所述技术问题的技术方案是，提供一种基于分形理论和蒙特卡罗方法的粗糙表面形貌模拟方法，其特征在于，该方法包括以下步骤：

步骤1、测量并提取粗糙表面的原始三维形貌数据；再根据原始三维形貌数据获得XYZ三维点云数据；再根据XYZ三维点云数据获得分形维数D、特征尺度参数G、实际表面最大峰谷差Sz和表面高度分布规律，并将这四个参数作为表面模拟的起始参数；

步骤2、以步骤1获得的起始参数为基础，采用分形理论和蒙特卡罗方法确定粗糙表面形貌结构堆叠状微凸体的尺寸和空间位置分布；

步骤3、遍历步骤2中生成的不同尺度的堆叠状微凸体的高度、基底直径和位置坐标，以此数据为基础得到粗糙表面的模拟三维表面形貌数据。

与现有技术相比，本发明的有益效果在于：

(1)本发明在保留自仿射和随机性特征的同时，采用实际表面最大峰谷差Sz和表面高度分布规律约束模拟表面，利用四个起始参数，解决了传统统计学的参数爆炸的问题，实现了试件粗糙表面形貌特征的模拟；同时，将分形参数与传统统计学参数相结合，解决了分形理论中参数不足导致不能完全表征粗糙表面形貌的问题；本发明不仅为加工零件的表面形貌模拟研究提供了新的思路，而且为后续接触模型的建立和分析提供了理论依据。

(2)本发明提出的表面形貌数值模拟方法便于编程实现，在PC机运算起来非常快，100,000次只需几分钟。

(3)本发明提出的表面形貌数值模拟方法，其独特之处在于使用了随机数来确定粗糙表面上的微凸体大小和位置，产生的表面上的微凸体大小和高度的分布具有随机性，可生成具有任意方向分布的非高斯表面，这比其它许多模型更符合实际情况。

(4)本发明通过对比实验结果，从表面幅值特性和空间分布两个角度分析了试件模拟和实际表面的形貌差异，验证了模拟结果的准确性。

附图说明

图1为本发明的方法流程框图；

图2为本发明实施例1的不同粗糙度的电火花加工方式试件的表面形貌图；

图3为本发明实施例1的不同粗糙度的喷砂加工方式试件的表面形貌图；

图4为本发明实施例1的变差法原理图；

图5为本发明实施例1的利用变差法计算电火花加工试件表面的分形维数D和特征尺度参数G；

图6为本发明实施例1的利用变差法计算喷砂加工试件表面的分形维数D和特征尺度参数G；

图7为本发明实施例1的试件表面以及堆叠状微凸体高度概率密度图；

图8为本发明实施例1的电火花加工表面形貌高度的概率密度图；

图9为本发明实施例1的喷砂加工表面形貌高度的概率密度图；

图10为本发明实施例1的模型模拟的电火花加工表面形貌图；

图11为本发明实施例1的模型模拟的喷砂加工表面形貌图。

具体实施方式

下面给出本发明的具体实施例。具体实施例仅用于进一步详细说明本发明，不限制本发明权利要求的保护范围。

本发明提供了一种基于分形理论和蒙特卡罗方法的粗糙表面形貌模拟方法(简称方法)，其特征在于，该方法包括以下步骤：

步骤1、测量并提取粗糙表面的原始三维形貌数据；再根据原始三维形貌数据获得XYZ三维点云数据；再根据XYZ三维点云数据获得分形维数D、特征尺度参数G、实际表面最大峰谷差Sz和表面高度分布规律，并将这四个参数作为表面模拟的起始参数；

优选地，步骤1中，使用三维形貌测量仪测量并提取粗糙表面的原始三维形貌数据。

优选地，步骤1中，根据原始三维形貌数据获得XYZ三维点云数据的方法是：将原始三维形貌数据经过降噪、精简和滤波处理后，获得XYZ三维点云数据。

降噪用于排除由于环境干扰或信号损失造成的数据存在的噪声；精简用于减小数据的大小和复杂度，避免为后续处理和存储造成负担；滤波用于去除数据中存在的一些异常点、离群值和不规则形状；经过降噪、精简、滤波保留了粗糙表面的原始形貌特征。

优选地，步骤1中，采用变差法计算获得分形维数D和特征尺度参数G的具体步骤如下：

S1.1、将宽度为R的矩形框置于XYZ三维点云数据沿X轴方向的二维轮廓截线上；

S1.2、移动矩形框以遍历整个二维轮廓截线，矩形框的高度由二维轮廓截线位于矩形框内的最高点和最低点决定，进而得到每一个矩形框的面积；再将所有矩形框的面积相加求和，得到矩形框总面积S(R)；

S1.3、改变宽度R的大小，重复S1.1和S1.2；

S1.4、作LnS(R)-LnR曲线，再进行线性回归得其拟合直线，拟合直线的斜率为W，得到分形维数D为：

D＝2-W(1)

进而得到特征尺度参数G为：

式(2)中，R

优选地，步骤1中，采用通过统计学方法计算获得实际表面最大峰谷差Sz和表面高度分布规律的具体方法如下：

其中，实际表面最大峰谷差Sz为：

Sz＝max(Z(x

式(3)中，Z(x

其中，表面高度分布规律是：表面形貌高度Z(x

步骤2、以步骤1获得的起始参数为基础，采用分形理论和蒙特卡罗方法确定粗糙表面形貌结构堆叠状微凸体的尺寸和空间位置分布；

优选地，步骤2中，分形理论和蒙特卡罗方法的具体步骤如下：

S2.1、确定堆叠状微凸体的尺寸；

S2.1.1、余弦波状的微凸体是构成粗糙表面的最基本单元，单频率γ

式(4)中，l为微凸体的基底直径；δ为基底直径l对应的微凸体的高度，δ＝G

S2.1.2、明确微凸体为余弦结构后，确定微凸体的基底直径l

微凸体基底尺寸分布遵循分形标度律，基于海平面岛屿分布规律和MonteCarlo方法，尺度i的微凸体的基底直径l

式(5)中，lmax和lmin分别对应最大和最小微凸体基底直径，lmax和lmin受到实际表面最大峰谷差Sz和表面高度分布规律的约束；R

S2.1.3、基于S2.1.2可获得i＝1～N尺度下的余弦波状的微凸体的基底直径和对应高度，以自仿射的形式，将最大高度的微凸体作为最底层，再将次最大高度的微凸体堆叠在最底层微凸体之上，依此类推，共生成N个不同层数的堆叠状微凸体；h

lmax最终可通过式(7)迭代求解：

S2.1.4、对S2.1.3获得的堆叠状微凸体进行高度分布修正，修正过程遵循表面高度分布规律，具体方式为：将不同尺度的堆叠状微凸体依据表面形貌高度划分若干区间，每个区间选取若干个堆叠状微凸体进行复制，复制倍数依据概率密度。选取复制后的总堆叠状微凸体数目N

S2.2、确定堆叠状微凸体的位置分布；

假设宏观模拟表面为一个边长为L

式(8)中，R

堆叠状微凸体之间在平面分布满足相互独立，若相邻堆叠状微凸体发生重叠，则重新生成位置坐标，直至全平面各个堆叠状微凸体之间不存在交集为止；具体的判断标准可通过比较相邻堆叠状微凸体圆心的距离和基底半径之和来实现：

步骤3、遍历步骤2中生成的不同尺度的堆叠状微凸体的高度、基底直径和位置坐标，以此数据为基础得到粗糙表面的模拟三维表面形貌数据。

为了验证本发明方法的正确性，对步骤3得到的模拟三维表面形貌数据进行模拟精度验证：分别计算步骤1中的XYZ三维点云数据的表面幅值特性和空间分布以及步骤3中的模拟三维表面形貌数据的表面幅值特性和空间分布，通过比较XYZ三维点云数据和模拟三维表面形貌数据的表面幅值特性和空间分布的误差来判断模拟效果。

优选地，表面幅值特性包含两个参数，分别为算数平均高度Sa和尺度限制的表面均方根高度Sq；空间分布包含两个参数，分别为尺度限制的表面均方根斜率Sdq和展开表面积比Sdr；四个参数的具体计算方法如下：

算数平均高度Sa：

式(10)中，m是采样区域A内沿x方向的离散点数；n是采样区域A内沿y方向的离散点数；Z(x

尺度限制的表面均方根高度Sq：

尺度限制的表面均方根斜率Sdq：

将式(12)离散化后，获得离散数学表达式为：

式(13)中，Δx和Δy为x和y方向上的采样间隔；

展开表面积比Sdr：

将式(14)离散化后，获得离散数学表达式为：

式(15)中，A

XYZ三维点云数据对应的是试件实际粗糙表面，因此误差的计算公式为：

实施例1

步骤1、使用三维形貌测量仪测量并提取8个镍合金试件表面的原始三维形貌数据，依据粗糙度不同对试件进行编号，8个镍合金试件中，4个为电火花加工方式，另外4个为喷砂加工方式，试件表面形貌如图2和图3所示；

实验中所用的三维形貌测量仪主要为两种仪器：LCuborAST非接触式工业3D相机，垂直分辨率为1μm，中心视野范围36×27mm，用于试件较大范围的轮廓提取；另一台LeicaDVM6数码显微镜拥有更高的垂直分辨率10nm，用于试件局部区域形貌的观察与分析；

将8个镍合金试件的原始三维形貌数据经过降噪、精简、滤波处理后获得XYZ三维点云数据；将XYZ三维点云数据通过变差法和统计学方法分别计算获得分形维数D、特征尺度参数G、实际表面最大峰谷差Sz和表面高度分布规律，共四个参数作为表面模拟的起始参数；

其中，变差法的计算方法如图4所示，选取不同宽度R的矩形框，通过遍历二维轮廓截线的方式获得LnS(R)-LnR曲线，再进行线性回归得其拟合直线，通过式(1)和式(2)计算获得电火花和喷砂两种加工方式下，镍合金试件表面的分形维数D和特征尺度参数G，如图5和图6所示；

利用式(3)计算XYZ三维点云数据的实际表面最大峰谷差Sz，计算结果如表1所示；

表1

表面形貌高度Z(x

步骤2、将步骤1计算获得的镍合金试件表面的分形维数D和特征尺度参数G作为起始参数，同时给定lmax和N的初值，分别为1×10

利用公式N＝[(l

再以自仿射的形式，将最大高度的微凸体作为最底层，再将次最大高度的微凸体堆叠在最底层微凸体之上，依此类推，共生成N个不同层数的堆叠状微凸体，利用式(6)可以计算出最大堆叠状微凸体高度h

参数h

确定了lmax后，需要对堆叠状微凸体进行高度分布修正，修正过程遵循表面高度分布规律，具体方式为：将不同尺度的堆叠状微凸体依据表面形貌高度Z(x

再将复制后的总堆叠状微凸体数目N

将上述高度分布修正过程运用在电火花和喷砂试件中，电火花和喷砂试件表面高度概率密度和模型表面的高度概率密度计算结果分别如图8和9所示；

随后，确定堆叠状微凸体的位置分布；根据式(8)和式(9)可确定N

步骤3、遍历步骤2的表面形貌特征点(包括堆叠状微凸体的尺寸和空间位置分布)，将特征点打包进三维形貌分析软件Gwyddion中，分别得到电火花和喷砂试件表面的模拟三维表面形貌数据，将数据分别绘制为电火花和喷砂表面形貌图，如图10和11所示；

对模拟数据进行精度验证：利用式(10)-(17)分别计算步骤1中的XYZ三维点云数据和步骤3中的模拟三维表面形貌数据的表面幅值特性和空间分布以及误差，计算结果如表2所示；

表2

由表2可以看出，电火花表面的模拟精度高于喷砂表面，且存在表面越平坦，模拟效果越好的特点；同时，本发明方法表面幅值特性的计算误差普遍低于空间分布，表面幅值特性两个参数Sa和Sq的误差小于8％，而空间分布两个参数Sdq和Sdr在电火花编号3和喷砂编号3和4的误差超过了10％；对于算数平均高度Sa分别为1.611μm、3.721μm、4.589μm和7.581μm的电火花试件，模拟表面和试件表面的表面幅值特性误差和空间分布误差均在12％以内。此外，Sa为2.508μm和3.342μm的喷砂试件，表面幅值特性误差和空间分布误差在5％以内；上述结果说明，采用本发明方法可以用来代替具有相同表面形态特征的实际加工表面，证明了本发明方法在粗糙表面建模的准确性。

本发明未述及之处适用于现有技术。