一种包含微观表面形貌的结合面接触性能分析方法

摘要

本发明公开了一种包含微观表面形貌的结合面接触性能分析方法，提高了网格划分质量和求解效率，改善了粗糙表面接触性能预测的准确性和可靠性。通过激光共聚焦显微镜测量实际表面或利用三维形貌数字化模拟方法获得粗糙表面，提取3D粗糙表面形貌的高度矩阵，利用点云处理方法将3D粗糙表面离散成有限元软件容易提取的高度矩阵文件，将其导入有限元软件中生成结合面每一连接表面的微观形貌，利用有限元软件关键点的平移和坐标修改功能，建立考虑微观表面形貌的体模型，通过网格控制的方法，对体模型进行六面体网格划分，构建两3D粗糙表面接触的有限元接触模型，逐步施加位移和力载荷边界条件，对结合面的接触特征进行分析。

说明书

技术领域

本发明涉及一种结合面接触性能仿真分析方法，特别是一种考虑两3D粗糙表面接触 的结合面性能仿真分析方法，可用于分析微观粗糙表面的摩擦、磨损、接触刚度、接触热 阻等接触性能。

背景技术

机械系统中存在大量的结合面，这些结合面的存在对机械系统整体特性有着很大影 响，而结合面的接触实际是其粗糙表面的接触，粗糙表面形貌对结合面摩擦、磨损及其接 触变形都有重要的影响。

长期以来，国内外学者从理论层面上对粗糙表面的接触开展了大量的研究。在理论解 析模型方面，构建了多种经典接触模型：GW模型、CEB模型、MB模型和KE模型等， 但这些理论模型的建立都存在一定的假设条件，如所有微凸体都具有相同的峰顶曲率半 径、微凸体接触时其变形相互独立，忽略微凸体之间的相互作用、不存在大变形、微凸体 之间的接触全部是峰对峰正接触等。正是由于这些假设，使其无法贴合结合面间实际的微 凸体接触状态，这极大影响了利用解析模型预测结合面性能的精度。

为了解决上述问题，近年来，有学者提出利用有限元数值模型模拟结合面间微观接触 的方法，但其主要集中在单对微凸体的接触性能分析和光滑表面—粗糙表面的接触性能分 析方面，尚未实现更贴合实际接触状态的两3D粗糙表面接触性能仿真分析。本发明从点 云处理技术、体模型生成技术和网格划分技术三个方面入手，克服了现有技术瓶颈，提高 了网格划分质量，构建了两3D粗糙表面接触模型，从而在一定程度上改善了结合面性能 预测的准确性。

发明内容

本发明的主要目的是提供一种包含微观表面形貌的结合面接触性能分析方法，克服了 两粗糙表面有限元接触模型构建中点云处理、体模型生成和网格划分技术瓶颈，提高了网 格划分质量和求解效率，改善了粗糙表面接触性能预测的准确性和可靠性。

本发明的技术方案是：

一种包含微观表面形貌的结合面接触性能分析方法，包含下述步骤：

1)通过激光共聚焦显微镜测量实际表面或利用三维形貌模拟方法获得高斯或非高斯 粗糙表面，提取3D粗糙表面形貌的高度矩阵为Z(m，n)，其中m、n分别表示x、 y方向取点的数目；

2)以步骤1)中获得的微观表面形貌为对象，选取适宜的间距，利用点云处理方法 将3D粗糙表面离散成有限元软件容易提取的txt文件；

3)利用步骤2)生成的三维点的txt文件，导入有限元软件中生成点云模型，利用有 限元的双重循环命令，连接X/Y方向相邻四点，生成微观面片群，获得结合面每 一连接表面的微观形貌；

4)以步骤3)生成的两组粗糙表面为基础，通过平面拉伸，并利用有限元软件关键 点的平移和坐标修改功能，构建一组考虑连接表面微观表面形貌的体模型；

5)以步骤4)生成的体模型为基础，通过线网格控制的方法，对体模型进行六面体 网格划分；

6)以步骤5)生成的网格模型为基础，定义结合面接触对，固定下连接体下表面z 方向的位移，通过上连接体上表面位移约束施加载荷边界条件，构建两3D粗糙 表面接触的有限元接触模型，并进行结构静力求解，观测连接结合面的接触载荷 和接触压力分布规律。

所述的利用三维形貌模拟方法获得高斯粗糙表面，具体包含以下步骤：

1)利用白噪声发生器randn(m，n)产生白噪声序列η(m，n)，并计算其傅里叶变换A(m， n)：

$A(k+1,l+1)=\underset{r-0}{\overset{m-1}{\Sigma }}\underset{s-0}{\overset{n-1}{\Sigma }}\eta (r+1,s+1)e^{-2\mathrm{\pi i}(\mathrm{kr}/m+\mathrm{ls}/n)}\left(\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}\right)$

2)利用给定自相关函数生成相关矩阵R(m，n)；假定自相关函数为f(x，y)，其中β_x、β_y分别为X、Y方向的自相关长度；

$f(x,y)=\exp (-2.3\sqrt{{\left(\frac{x}{{\beta }_x}\right)}^2+{\left(\frac{y}{{\beta }_y}\right)}^2})$

指定X、Y方向离散间距分别为Δ_x＝1μm、Δ_y＝1μm，在0≤x≤m/2，0≤y≤n/2区域内 对自相关函数进行离散，获得自相关函数的离散矩阵R(m/2+1，n/2+1)，然后扩展 为-m/2≤x≤m/2，-n/2≤y≤n/2区域内对应的自相关矩阵函数离散矩阵R(m，n)：

R(k+1，l+1)＝f(k，l)                      k＝0，1，2...0.5m，l＝0，1，2...0.5n

R(0.5m+2+k，l+1)＝f(0.5m-k，l+1)          k＝0，1，2...0.5m-2，l＝0，1，2...0.5n

R(k+1，0.5n+l+2)＝f(k+1，0.5n-l)          k＝0，1，2...0.5m，l＝0，1，2...0.5n-2

R(0.5m+k+2，0.5n+l+2)＝f(0.5m-k，0.5n-l)  k＝0，1，2...0.5m-2，l＝0，1，2...0.5n-2

3)对自相关函数矩阵R(m，n)进行傅里叶变换，获得模拟表面的功率谱密度P；

$P(1+k,1+l)=\left|\underset{r=0}{\overset{m-1}{\Sigma }}\underset{s=0}{\overset{n-1}{\Sigma }}R(r+1,s+1)e^{-2\mathrm{\pi i}(\mathrm{kr}/m+\mathrm{ls}/n)}\right|\left(\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}\right)$

4)由白噪声的功率谱密度为常数C，假设C＝1，计算白噪声到模拟高斯表面的传递 函数H(m，n)：

$H(m,n)=\sqrt{P}$

5)利用频域点积(A.*H)取反傅里叶变换的方法获得高斯表面的初始高度矩阵序列 z：

$z(k+1,l+1)=\frac{1}{\mathrm{mn}}\underset{r=0}{\overset{m-1}{\Sigma }}\underset{s=0}{\overset{n-1}{\Sigma }}(A.*H)e^{2\mathrm{\pi i}(\mathrm{kr}/m+\mathrm{ls}/n)}\left(\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}\right)$

6)根据给定高度标准偏差σ与模拟表面z的实际标准偏差std2(z)，求白噪声功率谱 密度常数C＝σ/std2(z)，从而模拟生成符合给定标准偏差等要求的高斯表面高度矩 阵z：z＝z*C.

所述的利用三维形貌模拟方法获得非高斯粗糙表面，具体包含以下步骤：

1)直接利用高斯粗糙表面的数字化模拟所生成的高度序列z，采用Pearson或Johnson 非高斯变换系统进行非高斯变换，生成非高斯序列Z；

2)如果所生成的非高斯序列的偏斜度S_k与峰度K_u不满足精度要求，则采用新的白 噪声序列，重复高斯粗糙表面高度序列的模拟及其非高斯转换，直到满足精度要 求，完成非高斯粗糙表面的数字化模拟。

所述的利用点云处理方法将3D粗糙表面离散成有限元软件容易提取的高度矩阵文 件，具体包含以下步骤：

1)对于实测或模拟表面，可以在水平方向按照一定的规律，生成包含3列的点文件， 其中每一行的三个数据代表一个点坐标，第1、2、3列分别表示对应点的x、y、 z坐标；

2)通过Ultraedit软件将点的格式修改为有限元软件识别的点创建命令，从而方便有 限元软件读取此类型的txt点文件，直接生成点云模型。

所述的利用有限元软件关键点的平移和坐标修改功能，构建一组考虑微观表面相貌的 体模型，具体包含以下步骤：

1)分别将两组微观表面沿Z正方向与反方向拉伸，生成具有一定高度的两组微六面 体群；

2)修改每一组微六面体群的点坐标，使得每一组六面体群中连接表面的对面顶点具 有相同的Z坐标值，形成光滑表面；

3)调节某一连接体的Z向位置，确保两连接表面有局部点刚发生初始接触。

所述的通过线网格控制的方法对体模型进行六面体网格划分，具体包含以下步骤：

1)网格单元数量的控制：在上述构建的体模型中选择厚度方向的所有线条，控制线 条高度方向的网格大小或网格层数，同时通过控制水平方向网格的大小，精确控 制网格的数量；

2)对同一个微六面体中，若四条高度方向的线条，高度相差悬殊，对其高度方向的 网格大小进行适当调整，使高度值大的网格数目与高度值小的线条对应网格数目 相差偶数，从而保证顺利生成高质量的六面体网格；

3)指定网格单元类型与材料属性，生成六面体扫描网格。

本发明提出一种考虑两3D粗糙表面接触的结合面性能仿真分析方法，从点云处理技 术、体模型生成技术和网格划分技术三个方面入手，克服了现有技术瓶颈，提高了网格划 分质量和求解效率，构建了两3D粗糙表面接触模型，从而在一定程度上解决了求解问题 中的难收敛性、改善了结合面性能预测的准确性和可靠性。以此，可有效指导结合面的设 计、加工和装配，使其满足结合面性能要求，从而保障机械装备整体性能。

附图说明

图1各向同性高斯表面及其高度分布(β_x＝β_y＝3，μ＝0，σ＝1.6，点数256×256)

图2各向异性高斯表面及其高度分布(β_x＝2，β_y＝100，μ＝0，σ＝1.6，点数256×256)

图3各向同性非高斯表面及其高度分布(β_x＝β_y＝5，μ＝0，σ＝1.6，S_k＝-1，K_u＝4.5，点数128×128)

图4各向异性非高斯表面及其高度分布(β_x＝2，β_y＝100，μ＝0，σ＝1.6，S_kz＝0.5，K_uz＝3.5，点数 256×256)

图5构成结合面的点云(128×128×2)

图6连接表面的面片群

图7上下连接表面的体拉伸、上下连接体群的顶点坐标修改以及连接体的Z轴平移

图8考虑两粗糙表面接触的六面体网格模型图

图9载荷不断增大对应的上连接体Z向应力变化规律

图10载荷不断增大对应的上连接体等效应力变化规律

图11载荷不断增大对应的结合面接触应力与接触区域变化规律

具体实施方式

以满足一定的自相关函数、标准偏差、偏斜度与峰度等统计特征的3D高斯和非高斯 粗糙表面的数字化模拟为例，对本发明的微观表面形貌获取方法进行说明；以一组高斯表 面形貌接触为例，采用ANSYS有限元软件，对本发明的结合面接触性能模型构建和分析 方法进行说明。

1.微观表面形貌获取：

利用三维形貌模拟方法获得高斯或非高斯粗糙表面，提取3D粗糙表面形貌的高度矩 阵为Z(m，n)，其中m、n分别表示x、y方向取点的数目，具体步骤如下：

(1)利用三维形貌模拟方法获得高斯粗糙表面，步骤如下：

1)利用白噪声发生器randn(m，n)产生白噪声序列η(m，n)，且获得其傅里叶变换

$n):A(k+1,l+1)=\underset{r-0}{\overset{m-1}{\Sigma }}\underset{s-0}{\overset{n-1}{\Sigma }}\eta (r+1,s+1)e^{-2\mathrm{\pi i}(\mathrm{kr}/m+\mathrm{ls}/n)}\left(\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}\right)$

2)利用给定自相关函数生成相关矩阵R(m，n)；假定自相关函数为f，其中β_x、β_y分 别为X、Y方向的自相关长度；

$f(x,y)=\exp (-2.3\sqrt{{\left(\frac{x}{{\beta }_x}\right)}^2+{\left(\frac{y}{{\beta }_y}\right)}^2})$

指定X、Y方向离散间距分别为Δ_x＝1μm、Δ_y＝1μm，在0≤x≤m/2，0≤y≤n/2区域内 对自相关函数进行离散，获得自相关函数的离散矩阵R(m/2+1，n/2+1)，然后扩展 为-m/2≤x≤m/2，-n/2≤y≤n/2区域内对应的自相关矩阵函数离散矩阵R(m，n)：

R(k+1，l+1)＝f(k，l)                      k＝0，1，2...0.5m，l＝0，1，2...0.5n

R(0.5m+2+k，l+1)＝f(0.5m-k，l+1)          k＝0，1，2...0.5m-2，l＝0，1，2...0.5n

R(k+1，0.5n+l+2)＝f(k+1，0.5n-l)          k＝0，1，2...0.5m，l＝0，1，2...0.5n-2

R(0.5m+k+2，0.5n+l+2)＝f(0.5m-k，0.5n-l)  k＝0，1，2...0.5m-2，l＝0，1，2...0.5n-2

3)对自相关函数矩阵R(m，n)进行傅里叶变换，获得模拟表面的功率谱密度P；

$P(1+k,1+l)=\left|\underset{r=0}{\overset{m-1}{\Sigma }}\underset{s=0}{\overset{n-1}{\Sigma }}R(r+1,s+1)e^{-2\mathrm{\pi i}(\mathrm{kr}/m+\mathrm{ls}/n)}\right|\left(\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}\right)$

4)由白噪声的功率谱密度为常数C，假设C＝1，计算白噪声到模拟高斯表面的传递 函数H(m，n)：

$H(m,n)=\sqrt{P}$

5)利用频域点积(A.*H)取反傅里叶变换的方法获得高斯表面的初始高度矩阵序列 z：

$z(k+1,l+1)=\frac{1}{\mathrm{mn}}\underset{r=0}{\overset{m-1}{\Sigma }}\underset{s=0}{\overset{n-1}{\Sigma }}(A.*H)e^{2\mathrm{\pi i}(\mathrm{kr}/m+\mathrm{ls}/n)}\left(\begin{array}{c}k=\mathrm{0,1,2}...m-1\\ l=\mathrm{0,1,2}...n-1\end{array}\right)$

6)根据给定高度标准偏差σ与模拟表面z的实际标准偏差std2(z)，求白噪声功率谱 密度常数C＝σ/std2(z)，从而模拟生成符合给定标准偏差等要求的高斯表面高度矩 阵z：z＝z*C。

利用上述给定的自相关函数，x、y方向的自相关函数均取3μm，模拟表面的x、y方 向取点为256×256、表面的标准偏差为1.6μm，模拟获得的各向同性高斯表面形貌及其 高度分布如图1所示。x、y方向的自相关函数分别取2μm、100μm，模拟获得的高斯表 面形貌表现出各向异性，具有一定的纹理，如图2所示。

(2)利用三维形貌模拟方法获得非高斯粗糙表面，步骤如下：

1)直接利用高斯粗糙表面的数字化模拟所生成的高度序列z，采用Pearson或Johnson 非高斯变换系统进行非高斯变换，生成非高斯序列Z；

2)如果所生成的非高斯序列的偏斜度S_k与峰度K_u不满足精度要求，则采用新的白 噪声序列，重复高斯粗糙表面高度序列的模拟及其非高斯转换，直到满足精度要 求，完成非高斯粗糙表面的数字化模拟。

利用上述给定的自相关函数，x、y方向的自相关函数均取5μm，模拟表面的x、y方 向取点为128×128、表面的均值、标准偏差、偏斜度与峰度分别为0、1.6μm、-1与4.5， 模拟获得的各向同性非高斯表面形貌及其高度分布如图3所示。若x、y方向的自相关函 数分别取2μm、100μm，模拟表面的x、y方向取点为256×256、表面的标准偏差、偏斜 度与峰度分别为1.6μm、0.5与3.5，模拟获得的各向异性非高斯表面形貌呈现出一定的纹 理，如图2所示。

2.微观形貌的离散与信息存储

利用步骤1微观形貌的获取方法，获得两128×128点数的高斯粗糙表面(上下表面 标准偏差分别为0.8μm与1.6μm)，x、y方向间距选取1μm，利用点云处理方法将3D粗 糙表面离散成有限元软件容易提取的txt文件，具体步骤如下：

1)对于实测或模拟表面，可以在水平方向按照一定的规律，生成包含3列的点文件， 其中每一行的三个数据代表一个点坐标，第1、2、3列分别表示对应点的x、y、 z坐标；

2)通过Ultraedit软件将点的格式修改为有限元软件识别的点创建命令“k，，X，Y，Z” (k相当于ANSYS中的点创建命令；“，，”两个逗号间省略了点的代号；X、Y、Z 分别表示点的坐标值，通过逗号隔开)，从而方便ANSYS软件读取这种类型的txt 点文件直接生成点云模型。

3.微观表面形貌的面模型构建

将步骤2生成的三维点的txt文件，导入ANSYS软件中生成点云模型(如图5所示)， 利用ANSYS的双重循环命令(*do，...，*end)，连接X/Y方向相邻四点，生成微观面片群， 获得结合面每一连接表面的微观形貌(如图6所示)。

4.微观表面形貌的体模型构建

以步骤3生成的两组粗糙表面为基础，通过平面拉伸，结合ANSYS关键点的平移和 坐标修改功能，构建一组考虑连接表面微观表面形貌的体模型，上下连接表面的拉伸、上 下连接体群顶点Z坐标的修改以及连接体群的Z轴平移如图7所示，具体步骤如下：

1)分别将上、下两组微观表面沿Z轴的正、反方向拉伸，生成具有一定厚度(厚度 大于两粗糙表面的最大轮廓高度，此处约为12μm)的两组微六面体群；

2)分别修改上、下两组微六面体群的上、下表面的顶点坐标，分别使得上、下两组 六面体群中的上、下表面的顶点具有相同的Z坐标值，形成光滑表面；

3)调节上一连接体的Z向位置，确保两连接表面有局部点刚好发生初始接触。

5.体模型网格划分

以步骤4生成的体模型为基础，通过线网格控制的方法，对体模型进行六面体网格划 分；

1)网格单元数量的控制：在上述构建的体模型中选择厚度方向的所有线条，控制线 条高度方向的网格大小或网格层数，同时通过控制水平方向网格的大小，精确控 制网格的数量，如图8的网格数量为127×127×4×2；

2)对同一个微六面体中，由于四条高度方向的线条，高度相差不大，在高度方向统 一定义4层网格单元；

3)指定网格单元类型为Solid45，上下材料为45钢，生成六面体扫描网格，如图8 所示。

6.结合面接触特性分析

以步骤5生成的网格模型为基础，定义结合面接触对，固定下端体模型下表面z方向 的位移，通过上端体模型上表面位移约束施加载荷边界条件(-0.5，0.8，11，1.4，...， 2.3，2.6μm)，构建两3D粗糙表面接触的有限元接触模型，并对两粗糙表面的接触中微 凸体的变形、Z向应力与等效应力以及结合面的接触压力与接触区域进行了分析。其中图 9所示为载荷增大时对应的上连接体Z向应力变化规律，图10所示为载荷增大时对应的 上连接体等效应力变化规律，图11所示为载荷增大时对应的结合面接触压力与接触区域 的变化规律。