结合面首次加载过程的法向接触刚度分形模型*

殷东华，张学良，温淑花，兰国生，陈永会

(1.太原科技大学机械工程学院，太原 030024；2.山西工学院智能制造产业学院，朔州 036000)

0 引言

结合面广泛存在于机械结构中，在动载荷作用下，结合面间会产生微小的相对运动，这使得结合面既存储能量又消耗能量，表现出既有弹性又有阻尼，即存在接触刚度和接触阻尼，并严重影响着机械结构的静动态特性[1]。结合面的接触刚度常常是机械结构整体刚度的重要组成部分。如对机床结构而言，机床振动问题中60%以上的问题源于结合面，并且结合面的接触刚度约占机床总刚度的60%～80%。许多学者针对结合面接触刚度进行了一系列的研究。在结合面接触刚度的研究中，最困难的任务之一是建立表征接触面并准确模拟接触行为的理论模型。MJUMDAR等[2]基于Weierstrass-Mandelbrot函数(W-M函数)首先提出粗糙表面接触的二维接触分形模型。YAN等[3]首先提出了粗糙表面接触的三维接触分形模型。张学良等[4]最早建立了结合面法向接触刚度分形模型。尤晋闽等[5-6]建立了结合面法向动态参数的分形模型，揭示了接触刚度和接触阻尼与作用在粗糙表面上的法向载荷、粗糙表面材料属性参数以及分形参数之间的复杂关系。随后，田红亮等[7]再次对结合面法向接触刚度模型进行了改进，但只考虑了弹性接触变形阶段和塑性接触变形阶段。然后，张学良等[8]提出了包含弹性变形、弹塑性变形和塑性变形情况的结合面法向接触刚度分形模型。MIAO等[9]提出一个完整的尺度相关的分形表面接触修正模型，研究了分形粗糙表面的总载荷、面积和刚度。田红亮等[10]基于修正分形几何学理论与Hertz法向接触力学方程，建立了微凸体与平面接触的法向接触刚度方程。张伟等[11]基于经典Greenwood-Willamson(GW)统计模型计算刚度的方法，建立了一种综合考虑微凸体基体变形和相互作用的结合面法向接触刚度模型。综上所述，关于固定结合面加卸载过程中的接触刚度分形理论建模研究工作非常少见。KOGUT等[12]基于有限元法建立了单个球体与刚性平面的弹塑性接触模型，简称KE模型，但该模型获得的接触载荷是不连续的，从而引起其刚度模型的不连续性。JACKSON等[13]采用比KE模型更细的网格，建立了新的单个微凸体与刚性平面的有限元模型，该模型考虑了硬度的变化，但没有给出开始发生塑性变形时的临界点。

因此，本文基于接触分形理论，采用Hermite插值法建立单个微凸体弹塑性变形阶段的接触面积模型，同时引入平均接触压力参数，弥补了单一微凸体接触载荷建模时存在的不连续缺陷，进而建立了接触刚度模型，借此提出一种固定结合面加载过程中的法向接触刚度分形模型，并通过仿真揭示了各相关参数对结合面法向接触刚度的影响，以及法向接触刚度和法向接触载荷的关系。本文研究为今后深入并准确地进行结合面加卸载法向接触刚度的研究提供了理论基础。

1 分形理论

两粗糙表面之间的接触往往可以等效为一个粗糙表面与一个刚性光滑平面的接触，进而等效成研究半个球形微凸体与刚性平面的接触，如图1所示。

图1 两粗糙表面接触示意图

根据等效粗糙表面的各向同性的分形特性，以及表征粗糙表面的二维轮廓曲线的Weierstrass-Mandelbrot函数(W-M函数)[14]，表面上任何一个微凸体变形前的轮廓曲线都可描述为：

(1)

式中，x为轮廓的位移坐标；D为表面轮廓分形维数，1

对于满足式(1)的微凸体，其峰顶曲率半径为：

(2)

微凸体的变形量ω与接触面积a的关系为：

ω=GD-1a(2-D)/2

(3)

根据Mandelbrot对海岸线上的岛屿面积分布性质，固定结合面之间微凸体的接触面积大小分布密度函数可以表达为：

(4)

式中，aL为微凸体的最大接触面积。

那么结合面间真实接触面积为：

(5)

单个微凸体从弹性向弹塑性转变时的临界变形量为[12]：

(6)

式中，H为较软材料的硬度，与材料的屈服强度σy相关，一般取H=2.8σy；K为硬度系数，取K=0.454+0.41υ，其中υ为较软材料的泊松比。

由式(2)和式(6)可得临界接触变形量：

(7)

比较ω和ωc得到：

(8)

当ω=ωc时，所对应的接触面积即为微凸体从弹性向弹塑性转变时的临界接触面积ac，可表达为：

(9)

因此，

(10)

2 单个微凸体加载过程的力学模型

2.1 弹性接触阶段

当ω≤ωc时，微凸体发生弹性变形。根据Hertz理论，刚性平面与峰顶曲率半径为R的微凸体的弹性接触中，单个微凸体的接触面积ae、接触载荷Fe、最大接触压力pm和平均接触压力pea与接触干涉ω之间的关系可以表示为：

ael=πRω

(11)

(12)

pm=(2E/π)(ω/R)1/2

(13)

(14)

由式(12)可得单个微凸体的法向接触刚度为：

(15)

2.2 全塑性接触阶段

当ω≥110ωc时，微凸体发生完全塑性变形，其平均接触压力达到硬度值H并保持恒定。根据文献[15]，单个微凸体的接触面积apl、接触载荷Fpl和平均接触压力ppa可以表示为：

apl=2πRω

(16)

Fpl=Hapl=2πRωH

(17)

ppa=H

(18)

但当去除载荷后，发生完全塑性变形的微凸体没有恢复过程，故忽略此阶段的接触刚度。

2.3 弹塑性接触阶段

当ωc≤ω≤110ωc时，微凸体发生弹塑性变形。根据文献[16]，在微凸体的弹塑性变形范围内，平均接触压力pepa与接触干涉ω之间的关系可以表示为：

(19)

接触面积aep和接触干涉ω之间的关系可以用一个多项式来表示。通过满足以下边界条件：①在ω=ωc处aepl=ael，daepl/dω=dael/dω；②在ω=110ωc处aepl=apl，daepl/dω=dapl/dω，本文采用Hermite插值法对定义域为ωc≤ω≤110ωc的aepl-ω函数进行逼近得到：

(20)

然后利用式(19)和式(20)，在弹塑性变形状态下微凸体的接触载荷可以表示为：

(21)

对KE模型和式(21)进行数值仿真，得到两种模型下接触载荷和变形量的关系，如图2所示。

图2 两种模型下接触载荷和变形量的关系

由图可知，KE模型中的接触载荷在ω/ωc=6处是不连续的，而且众所周知接触载荷在ω/ωc=1和ω/ωc=110处也是不连续的。而本文模型弥补了这种不连续性缺陷。

由式(21)可得加载过程中单个微凸体的法向接触刚度为：

(22)

3 结合面加载过程的分形模型

结合微凸体的接触面积大小分布密度函数，加载过程中结合面的法向接触载荷和法向接触刚度分别为：

(23)

(24)

对上式进行无量纲化得：

(25)

(26)

式中，

4 仿真及结果

(c) φ=2.5

基于考虑微接触加载过程中的接触载荷和接触刚度分析，本文结合粗糙表面的分形参数，分析了整个结合面接触刚度模型对接触载荷的依赖关系。由本文模型得到的加载过程中法向接触刚度与接触载荷的关系如图4所示，结合面无量纲法向接触刚度随着法向接触载荷的增大而增大，但在不同的分形维数和特征尺度系数的范围内其变化速率不同。根据文献[17]中实验例证得到的结论：由于实际接触面积的增加，接触刚度随接触载荷的增加而增加。可见，本文模型与结合面实际情况相符合。

(a) 1.1≤D≤1.3 (b) 1.3≤D≤1.5

(c) 1.5≤D≤1.7 (d) 1.7≤D≤1.9

5 结论

基于接触分形理论，本文采用Hermite插值法建立单个微凸体弹塑性变形阶段的接触面积模型，同时引入平均接触压力参数，弥补了单一微凸体接触载荷建模时存在的不连续缺陷，进而建立了接触刚度模型，借此提出一种固定结合面加载过程中的法向接触刚度分形模型。得到结论如下：

(1)加载过程中固定结合面无量纲法向接触刚度随着分形维数的增大而增大，且分形维数取值范围不同时增长速率不同。

(2)加载过程中固定结合面无量纲法向接触刚度随着特征尺度系数的增大而减小，但是在分形维数取值1.1～1.4时，特征尺度系数的变化对无量纲法向接触刚度的影响不明显。

(3)加载过程中固定结合面无量纲法向接触刚度随着塑性指数的增大而增大，且随法向接触载荷的增大而增大。