基于ANSYS 随机粗糙表面接触分析

周朝朋，杨照

（200093 上海市 上海理工大学 机械工程学院）

0 引言

随着国民经济不断发展，高新技术快速更新迭代，涌现了以光刻机、数控设备、机器人、航空航天设备等为代表的现代装备，为了满足人民的需求，设备也越加复杂化，所有的装备都离不开接触面，而接触面的特性对于机械设备的正常运行有着极大的影响。1966 年，Williamson 和Greenwood[1]提出了著名的GW 接触模型，模型假设粗糙表面的微凸体的分布符合高斯分布，并将结合部等效为光滑的刚性面和粗糙的表面；在国内，张学良[2]采用螺栓固定的悬臂梁模型，通过敲击法获得了模型的固有频率，利用MATLAB 软件识别随着面压载荷的变化法向接触刚度的变化规律，经过计算获得不同面压载荷下的单位结合面上的法向接触刚度的数值。但由于加工工艺的多样性和复杂性，到目前为止还没有对接触面有明确的理论基础，本文以MATLAB 生成模拟随机粗糙表面关键点数据并导入仿真软件中进行有限元分析，就接触力、接触面积随着位移的关系进行研究。研究表明，接触力、接触面积随位移呈三次函数关系，并对之后研究接触面特性有积极的参考价值。

1 随机粗糙表面生成

1.1 生成粗糙表面数据函数

为了模拟复杂的接触表面，经研究得知，加工表面凸起峰值符合高斯分布，现阶段生成表面数据方式有2 种：第1 种利用Rand（）函数，y=sqrt（0.09）×rand（1 000，1）即可生成1 000个点方差为0.09 的随机数据；第2 种利用wgn（）函数，y=wgn（1000，1，10×log10（0.09））即生成1 000 个点方差为0.09 的符合高斯分布的随机数据。本文采用第2 种方式。

1.2 生成随机表面数据

使用MATLAB 软件，写M 文件并运行生成随机表面数据[3]。图1、图2 所示分别为粗糙度为Ra=3.2、Ra=6.3 的轮廓数据图，表1 为局部数据表。从图中可以看出随机数据的复杂性。为了研究生成数据的分布规律，绘制数据分布规律直方图，如图3、图4 所示。由统计学来看，数据分布具有规律性，完全可以模拟符合高斯分布的复杂曲面，分析结果也具有代表性。

图1 粗糙度为3.2 轮廓数据图Fig.1 Contour data diagram with roughness of 3.2

图2 粗糙度为6.3 轮廓数据图Fig.2 Contour data with roughness of 6.3

表1 MATLAB 生成的粗糙表面数据点Tab.1 Rough surface data points generated by MATLAB

图3 第1 组白噪声分布直方图Fig.3 Histogram of white noise distribution

图4 第1 组白噪声分布箱线图Fig.4 Boxed diagram of white noise distribution

2 Ansys 接触特性分析

2.1 模型的建立

将MATLAB 生成的为粗糙表面的关键数据点直接导入ANSYS 中，在软件中利用Spline 函数生成光滑的曲线。图5 所示为粗糙表面局部图。光滑的曲线可以避免在之后的网格划分中畸形网格对分析结果产生影响。在建好光滑曲线之后，通过K，L 命令建立模型其他的轮廓线，由于直线的建立较为简单，此处不再赘述。如图6 为建立好的关键线轮廓。在粗糙接触面曲线和关键轮廓线建好之后，利用AL 命令，建立分析所需要的面，如图7 所示。

图5 ANSYS 生成的光滑曲线Fig.5 Smooth curve generated by ANSYS

图6 关键线轮廓线图Fig.6 Contour diagram of key lines

图7 模型面图Fig.7 Model plane diagram

2.2 定义模型属性及网格划分

本次分析采用45#，其弹性模量：E=207 GPa，屈服强度：σs=0.355 GPa，泊松比：υ=0.3，密度：ρ=7 850 kg/mm3。分析类型为平面应变，使用平面182 单元。由于此次分析为接触分析，分析的关键在接触面[4]，为了节约计算空间且对计算结果不产生影响，所以在Y 方向采用分层并网格逐步加密的方法。图8 所示为生成的网格图，网格数：313 446；网格节点：306 011。

图8 网格图Fig.8 Grid diagram

2.3 接触与约束

分析的塑性选项采用双线性等向强化（BISO），此屈服模式下的 Von Mises 屈服准则也是适用的。随机粗糙表面与刚性平面的接触采用拉格朗日法则建立接触关系。在粗糙模型Y 方向定义位移约束，为了清晰地得到接触力、接触面积随位移的变化趋势，采用300 个等步长加载，并关闭自动判断选项，增加步长的同时也使求解结果更易收敛，且结果采用Von Mises 判断依据。

2.4 结果分析

以粗糙度为Ra=3.2 为例，图9、图10 为接触分析结果位移云图和局部接触应力图。

图9 位移云图Fig.9 Displacement cloud map

图10 局部接触应力图Fig.10 Local contact stress diagram

图11、图12 为应力应变曲线。从曲线可以看出，当应力达到0.3 GPa 时，零件由弹性阶段转变为弹塑性阶段，即随机粗糙表面接触前期处于纯弹性阶段，达到屈服应力时，此时零件既有弹性变形，也有塑性变形[5]。

图11 应力应变曲线Fig.11 Stress-strain curves

图12 应力应变曲线局部放大图Fig.12 Local enlargement of stress-strain curves

以压入深度为横坐标，分别以接触力、接触面积为纵坐标，做压入深度-接触力关系图和压入深度-接触面积关系图，如图13、图14 所示。

图13 压入深度和接触力曲线Fig.13 Pressing depth and contact force curve

图14 压入深度和接触面积曲线Fig.14 Pressing depth and contact area curve

以式（1）为拟合函数拟合[6]曲线，得到压入深度接触力拟合参数，见表2。

表2 压入深度接触力拟合参数Tab.2 Fitting parameters of pressing depth and contact force

由表2 可知，拟合精度完全符合要求。

表2 中数据分别代入式（1），可得：

3 结论

通过ANSYS 仿真分析，得到粗糙度分别为Ra=3.2，Ra=6.3，Ra=12.5 的随机粗糙表面下的接触应力、接触面积随着位移的变化关系，并根据变化曲线以三次函数拟合得到拟合曲线、拟合方程。由拟合误差来看，拟合符合要求，可为今后研究接触特性提供数据支持和参考。