分形接触干摩擦模型及在围带叶片受迫响应分析中的应用

邱恒斌 徐自力

（1.北京强度环境研究所，北京100076）（2.西安交通大学机械结构强度与振动国家重点实验室，西安710049）

引言

在透平机械中，叶片系统经常采用围带结构增加叶片间的干摩擦阻尼来达到减振的目的.为了研究围带干摩擦阻尼叶片的振动特性，建立一个高效准确的干摩擦模型是十分必要的.

在已有干摩擦阻尼系统的研究中，微滑移模型由于能较好的模拟实际摩擦接触行为而被广泛采用.Iwan［1］采用线性弹簧和库仑摩擦元组合，在接触面分布多个接触点对，提出串联和并联微滑移模型.Griffin［2］基于库仑摩擦定律，采用恒定摩擦系数，并采用线性弹簧表征切向接触刚度，研究了恒定正压力下的一维摩擦相对运动.Yang［3］提出了二维的干摩擦接触模型，采用两个无质量的弹簧表征切向和法向的接触刚度，研究了变正压力下的二维摩擦相对运动.Chen，Yang和Menq［4］在二维模型基础上提出了三维的摩擦接触模型 .Petrov［5，6］采用三维摩擦模型研究了多种形式的干摩擦阻尼结构的减振特性.在已有的微滑移模型中，通常采用线性弹簧模拟接触面发生相对滑移之前的弹性变形，用库仑摩擦元模拟接触面发生滑移时的摩擦力.

接触刚度和摩擦系数是计算摩擦力的两个重要参数，在先前大部分摩擦接触模型中，摩擦系数和接触刚度大都采用经验值.实际上，摩擦系数和接触刚度均会随着接触正压力，接触面形貌，弹性模量，泊松比等参数的变化而改变.在计算接触刚度和摩擦系数过程中，接触正压力，弹性模量，泊松比等参数较容易获得，但是准确的模拟接触面形貌是较为困难的.

在以往的接触研究中，大多采用GW（Greenwood-Williamson模型）和分形理论模拟接触面形貌.Chang等［7］基于GW接触模型，提出金属表面静摩擦系数的计算模型，得出摩擦系数与正压力的关系.尤晋闽等［8］研究了分形接触表面的摩擦模型，得出了静摩擦系数与接触面材料参数、接触面分形参数和正压力的关系.Jones［9］利用GW模型研究了法向载荷、切向载荷、表面形貌等对切向刚度的影响.Jiang等［10］采用分形几何模拟接触面形貌，研究了机械结合面的接触刚度，理论计算结果与实验结果相接近.GW模型获得的表面特征参数受限于分辨率及采样长度等，不能准确反映出粗糙表面的全部信息［11］.粗糙表面轮廓具有自相似性，且与尺度无关，因此可以由分形几何来表征，利用表面分形特性建立接触模型，可使表面接触的分析结果具有确定性和唯一性［11，12］.

本文基于三维干摩擦微滑移模型，发展了一种分形接触干摩擦微滑移模型，用于计算考虑接触面形貌的摩擦力.采用分形几何模拟接触表面的形貌，基于分形理论和赫兹接触理论，建立接触面粗糙度、正压力、弹性模量、泊松比等参数与接触刚度和摩擦系数的关系.由接触刚度、摩擦系数和接触点对的相对位移计算摩擦力.采用该发展的模型研究真实围带叶片的受迫振动响应.

1 分形接触

1.1 分形接触表面

分形具有自相似性和尺寸独立性，使其能够描述具有多尺度特性的粗糙表面.具有分形特征的三维粗糙表面轮廓曲线可由W-M（Weierstrass-Mandelbrot）函数来模拟，该函数具有处处连续，处处不可导，具有无限细节特征［13］.该数学表达式为：

式中：z(x)为粗糙表面轮廓的高度；x，y为粗糙表面轮廓的位移坐标；G为反映表面轮廓z的特征尺度系数，G越大越粗糙；D为表面轮廓的分形维数，表征表面轮廓的不规则性和复杂程度，且1

图1 三维分形接触表面模拟Fig.1 Simulated 3D fractal surfaces

1.2 接触刚度

围带接触面的接触刚度源于接触表面的一系列微凸体，研究表明［13］可将两个粗糙表面的接触等效为一个粗糙表面和一个光滑的刚性平面接触的如图2.则该系统的等效弹性模量为：

图2 粗糙表面微凸体与刚性平面的接触Fig.2 Contact between an asperity of composite rough surface and a rigid plane

式中v1、v2和E1、E2分别为两个接触微凸体的泊松比和弹性模量.微凸体变形量由式（1）W-M分形函数［13］决定：

Majumdar和 Bhushan［12］研究表明接触点的面积分布密度函数与海洋面上的岛屿面积分布密度函数相似，得到切割面积a'的统计学分布函数为：

式中，a'=π(r')2是微接触的切割面积，al'是最大微接触切割面积.微凸体的曲率半径为R2=(R-δ)2+(r')2，忽略二阶小量可写为(r')2=2Rδ，代入到式（3）可得：

分形维数D、分形参数G与接触面粗糙度的关系可近似为［14］：D=1.515∕Ra0.088，G=10(-8.259∕Ra0.088).Kogut等［15］将微凸体的变形可分为四个阶段：完全弹性、弹塑性和完全塑性变形.即当δ<δc时微凸体处于完全弹性状态；当δc≤δ≤6δc时微凸体处于第一弹塑性状态；当6δc<δ≤110δc时微凸体处于第二弹塑性状态；当δ>110δc时微凸体处于完全塑性状态.微凸体处于各个接触状态的接触载荷如表1所示［15］.

表1 接触变形及对应的接触载荷Table 1 Contact deformation and corresponding contact load

微凸体接触面由完全弹性到弹塑性变形的临界变形量为：

式中，H是材料的硬度，K是材料的硬度系数，K=0.454+0.41v，v是软材料的泊松比.将式（3）和式（5）代入式（6），微凸体由完全弹性到弹塑性变形的临界切割面积可表示为：

整个接触面的接触刚度可认为是处于完全弹性、第一和第二弹塑性区的单个微凸体接触刚度的总合.因此整个接触面的法向接触刚度可表达为：

同样的，整个接触面的切向接触刚度可表达为：

切向和法向接触刚度随正压力和粗糙度的变化如图3和图4所示.由图可知，接触刚度随着正压力的增加而增加，随着粗糙度的增加而减小.接触表面越粗糙意味着接触表面越尖锐，接触表面微凸体在受到接触载荷时越容易发生塑性变形，因此接触刚度会减小.

图3 接触刚度随正压力的变化Fig.3 Contact stiffness vs normal pressure

图4 接触刚度随粗糙度的变化Fig.4 Contact stiffness vs roughness

1.3 摩擦系数

当微凸体同时受到法向载荷F和切向载荷T共同作用时，微凸体表面的应力场目前还没有完整的解析解，现有文献大多基于一定的假设条件对表达式进行简化.本文采用文献［8］中提出的假设，认为微凸体的屈服发生在微凸体表面接触点的边缘，接触点边缘的应力为：

式中σs为较软材料的屈服强度.在法向载荷作用下，接触界面上处于塑性变形的微凸体由于受到局部接触载荷作用而发生塑性流动，将不能继续承受切向载荷［17］.因此在计算静摩擦力只包括前两个状态的微凸体，当接触面处于临界滑动状态时，接触面将达到完全屈服，此时整个接触面所有微凸体所承受的切向载荷即为最大静摩擦力Tt：

接触面摩擦系数随正压力和粗糙度的变化如图5和图6所示.由图可知，摩擦系数随着正压力的增加而增加，随着粗糙度的增加而先增大后减小.

图6 摩擦系数随粗糙度的变化Fig.6 Contact stiffness vs roughness

2 干摩擦模型

叶片振动时，围带接触面会发生三维相对运动，产生摩擦约束力.为了模拟实际摩擦接触过程中的微滑移现象，将围带接触区域被离散成一系列接触单元，每个接触单元由一组接触点对来表征该接触单元在摩擦接触运动过程中粘滞、滑移和分离的摩擦运动状态.

接触面上单个接触点对的模型如图7所示.接触点对的接触刚度由式（10）和式（11）确定，摩擦系数由式（15）确定.A，B两点分别为互相接触的两个平面上对应的接触点，用两个无质量的线性弹簧表征接触刚度.弹簧Kt表征接触面的切向接触刚度，弹簧Kn表征接触面间的法向接触刚度.当接触面做三维相对运动时，假设点A与面1始终保持粘滞，而点B会沿着面2产生粘滞或滑移的运动状态.在面2的接触面上建立局部坐标系，在局部坐标系中，点A和点B的位移分别定义为和，点A和点B之间的相对位移定义为

图7 单个摩擦接触点对模型Fig.7 The model of a friction pair

接触面间的法向接触正压力由接触面初始正压力n0、法向接触刚度kn和法向相对位移zk确定，可由下式来计算：

若接触面初始正压力n0为负，表示接触面存在初始间隙；等于零，则表示接触面在此时刻发生分离，接触点摩擦力为零.接触面摩擦约束力就可由点A和点B间的相对位移计算获得.切向摩擦力fkt定义为：

定义点A与点B之间的切向拉力Tk为：

如果Tk小于临界摩擦力，则点B保持不动，若Tk大于，那么点B将向点A滑动直至拉力Tk等于临界摩擦力.则切向摩擦力可以表示为：

3 围带叶片响应

某真实围带叶片的有限元模型如图8所示，采用该发展的分形接触三维干摩擦模型和高阶谐波平衡法［18］预测该围带叶片的受迫振动响应.对叶片轮盘施加固定约束，在叶片轴向上施加p(t)=50sinωt的激振力，接触面刚度及摩擦系数采用分形接触模型计算，选取前5阶谐波和前5阶模态进行响应计算.为了研究接触面粗糙度对围带叶片减振效果的影响，本文选取了表面粗糙度为Ra=1.6μm和Ra=2.2μm的两种情况，计算两种表面粗糙度在不同初始正压力下叶片系统的振动响应，两种粗糙度下叶片系统的幅频曲线如图9和图10所示.

图8 围带叶片的有限元模型Fig.8 The finite element model of a shroud blade

图9 Ra=2.2μm时的幅频曲线Fig.9 Amplitude-frequency curves at Ra=2.2μm

由图9和图10可知，当接触面粗糙度不同时，叶片的幅频曲线、共振幅值和共振频率均发生变化，说明叶片的振动响应对接触面粗糙度较为敏感.就本文选取参数计算所得结果，在相同的初始正压力下，表面粗糙度1.6μm计算的共振幅值较表面粗糙度2.2μm的小.同时发现该真实围带叶片的响应没有出现摩擦锁死的动力学现象，这是因为实际机组为了在振动过程中达到减振的目的，会选取合适的初始正压力，防止出现摩擦锁死现象，这与机组实际运行的情况相符.

图10 Ra=1.6μm时的幅频曲线Fig.10 Amplitude-frequency curves at Ra=1.6μm

4 结论

本文基于三维干摩擦微滑移模型，发展了一种分形接触干摩擦微滑移模型，用于计算考虑接触面形貌的摩擦力.在该模型中，摩擦接触面被离散成一系列接触单元，每个接触单元由一组接触点对来表征该接触单元的粘滞、滑移、分离的摩擦运动状态.采用分形几何模拟接触面形貌，基于分形理论和赫兹接触理论，建立接触面粗糙度、正压力、弹性模量、泊松比等参数与接触刚度和摩擦系数的关系.由接触刚度、摩擦系数和接触点对的相对位移计算摩擦力.

结果表明：该分形接触三维干摩擦模型能较好地模拟接触面形貌，有效地考虑接触面形貌、正压力等参数对摩擦力的影响.其中接触刚度随着正压力的增加而增加，随着粗糙度的增加而减小；摩擦系数随着正压力增加而增加，随着粗糙度的增加，先增大后减小；存在最优的初始正压力使叶片共振幅值最小，共振频率随着初始正压力的增加而增加；接触面形貌对围带叶片受迫振动响应有显著影响.