经过装配的部件间的接触面称之为结合部.结合部大量存在于机械系统中, 且系统的60%刚度和90%的阻尼皆产生于结合部^[1-2].正常工作的机器, 其结合部由于外载传递将不可避免地产生摩擦现象, 反过来摩擦又对系统动态稳定性产生影响^[3-5].因为工业界约有30%~50%的能量会通过热传导、振动、噪声等形式消耗在结合部摩擦上^[6]，因此, 预测和研究结合部的静摩擦特性是必要的.摩擦因数的概念首次由达芬奇给出并引入理论研究范畴, 至此很多学者进行了结合部静摩擦因数特性的研究.Chang等^[7]考虑结合部上粘着力的影响, 并基于统计的Greenwood-Williamson(G-W)模型, 给出了结合部静摩擦因数预测模型.田红亮等^[6]改进了分形理论公式推导, 得出了最大静摩擦力及摩擦因数等与微凸体实际接触面积间的函数预测模型; 兰国生等^[8]根据传统Majumdar-Bhushan (M-B)模型, 给出了静摩擦因数与法向载荷之间的函数变化趋势.

经典的G-W模型受限于测量仪器分辨率和采样长度, 所以推导出的相关模型^[7]具有局限性.M-B分形模型弥补了这一局限, 因其对表面形貌进行分形表征具有尺度独立性, 将不再受仪器分辨率和采样长度限制.但是以上基于分形理论建立的静摩擦因数模型^[6-8]均是建立在二维分形曲线(1 < D < 2) 基础上, 用二维分形模型计算预测三维形貌(2 < D < 3) 的接触特性, 本身就存在理论上的不合理性^[9-10].因此本文旨在修正的M-B分形接触理论和赫兹接触理论基础上, 考虑三维分形接触特性, 建立了新的结合部静摩擦因数非线性数学模型.

1 结合部建模

由Hertz理论可知, 结合部上两微凸体在外载下相互接触可等效为一刚性光滑平面和一球形等效微凸体相接触, 如图 1所示.

图 1中三角形obc根据勾股定理易得

(1)

式(1) 变形可得

(2)

由于微凸体变形量ω远小于R, 故设, 可得

(3)

图 1中没有变形的等效微凸体和刚性光滑平面接触的横截面积为

(4)

由经典赫兹弹性接触理论知, 微凸体处于弹性变形阶段时所承受的法向载荷为^[6]

(5)

当微凸体接触时, 实际变形的微凸体与刚性平面所形成的接触面的半径为^[11]

(6)

(7)

式中:E为微凸体的等效弹性模量; E₁和E₂分别为两相互接触的微凸体的弹性模量; ν₁和ν₂则分别为它们的泊松比.

将式(5) 代入式(6) 可得

(8)

将式(3) 代入式(8) 可得实际变形的微凸体与刚性平面所形成的接触面半径r与等效刚性平面截得等效微凸体后所形成的横截面半径r′之间的关系:

(9)

由式(4) 和式(9) 可知, 图 1中等效微凸体和刚性光滑平面实际接触面积为

(10)

1.1 各变形阶段微凸体的临界接触面积

Yan等^[11]改进单变量的W-M函数, 得到了能模拟三维表面形貌的曲面函数.修改后的双变量W-M函数为

(11)

式中:L为形貌取样长度; D为结合部形貌分形维数(2 < D < 3);G为分形尺度系数; γ(γ > 1, 一般取1.5) 为决定频率密度的参数; M为组成粗糙表面形貌时峰脊叠加的数目; n为频率指数, n_max=Int[lg(L/L_s)/lgγ]; L_s为最低截断长度; x, y为表面微凸体的直角坐标系; φ_{m, n}为乱数相位.

由式(11) 知, 该余弦函数的波峰与波谷幅值差可用来定义微凸体受载后的变形量, 即

(12)

由式(4) 和式(10) 将式(12) 变形为面积a的函数:

(13)

将式(13) 代入式(10) 中, 得等效微凸体曲率半径:

(14)

弹性变形阶段微凸体的临界变形量为^[12]

(15)

式中:为材料特征参数; σ_y为较软材料的屈服强度; k_μ为摩擦力的修正参数, μ为动摩擦系数.当0≤μ≤0.3时, k_μ=1－0.228μ; 当0.3 < μ < 1时, k_μ=0.932 e^{－1.58(μ－0.3)}.

由式(13)~式(15) 可推导出微凸体的临界弹性变形面积:

(16)

1.2 各变形阶段微凸体的法向接触载荷

由式(5)、式(13) 与式(14) 可将微凸体弹性变形阶段法向载荷表示成接触面积a的函数:

(17)

微凸体塑性变形阶段法向载荷为

(18)

式中:λ=H/σ_y, 为定义的系数, H为较软材料的硬度.

1.3 结合部的实际接触面积

结合部上所有微凸体的面积分布函数n(a)^[11]与其中最大的微凸体接触面积a_l间的关系为

(19)

对面积分布函数n(a)在微凸体各变形区间进行连续积分, 可得整个结合面上的实际接触面积A_r.它应包含结合部上微凸体经历从弹性变形到塑性变形阶段的接触面积总和:

(20)

2 法向接触载荷三维分形模型

对单个微凸体在其弹性和塑性变形阶段与面积分布函数n(a)进行连续积分, 由式(17)~式(19) 可得总法向载荷.

1) 当2 < D < 3, D≠2.5时,

(21)

2) 当2 < D < 3, D=2.5时,

(22)

对式(21) 无量纲化可得

(23)

式中:P^*=P/(EA_a); a₁^*=a₁/A_a; a_c^*=a_c/A_a; ; A_a为结合面名义接触面积.

对式(22) 进行无量纲化可得

(24)

3 切向接触载荷三维分形模型

根据汉米尔顿假设^[13], 若微凸体的屈服发生在接触点边缘, 接触边缘的3个主应力分别为

(25)

(26)

(27)

式中:Q, p分别为单个微凸体受到的切向与法向载荷.

由Tresca屈服条件知:

(28)

将式(25)、式(26) 代入式(28), 得

(29)

由式(29) 和式(19) 可知所有处于弹性变形阶段微凸体的切向载荷.

1) 当2 < D < 3, D≠2.5时,

(30)

2) 当2 < D < 3, D=2.5时,

(31)

对于已经发生完全塑性变形的微凸体, 没有承受切向载荷的能力, 故在计算结合部切向载荷时不涉及微凸体完全塑性变形阶段.那么结合部总切向接触载荷为

(32)

对式(30) 和式(31) 分别进行无量纲化可得

(33)

式中, T^*=T/EA_a.

(34)

4 静摩擦因数非线性分形模型

随着切向力的增大, 结合部将处于临界滑动状态, 此时结合部上微凸体所受的总切向载荷T就是最大静摩擦力^[6-8].

因此, 由式(21)、式(22) 及式(30)~式(32) 可得考虑真实三维形貌静摩擦因数的非线性分形模型:

(35)

根据式(23)、式(24)、式(33)、式(34) 可得其无量纲化模型:

(36)

5 仿真与讨论

分别分析静摩擦因数在不同D和不同G下的变化规律; 对比分析了二维分形和三维分形下静摩擦因数模型间的差异, 以及理论计算值与试验值间的误差.

5.1 D对三维静摩擦因数的影响

令φ=1, ν=0.3, G^*=10^－5, λ=2.8, 由式(35) 可得D对静摩擦因数的影响, 如图 2所示.可知:静摩擦因数f均随无量纲法向载荷P^*的增大, 单调增大; 当D小于2.5时, 静摩擦因数随分形维数的增大而增大; 当D大于2.5时, 静摩擦因数随分形维数的增大而减小.

5.2 G^*对三维静摩擦因数的影响

令D分别为2.3, 2.5, 2.7时, 分形尺度系数G^*对f的影响如图 3所示.可知:f随着G^*的减小而增大; f的增长率随G^*的增大而减小.

5.3 本文模型的验证

为了验证文中所提模型的有效性，将其与现有试验数据进行对比.田红亮等^[14]以YK31320型滚齿机床的静压导轨结合面为研究对象, 给出了结合面上静摩擦因数随不同法向载荷的变化值.结合面相关参数如表 1所示.结合面的分形参数D, G分别为2.503 3, 2.006 6×10^-9 m.当结合面法向载荷P=100 kN时, 由表 1中数据及式(15) 和式(20) 可得此时微凸体的最大接触面积a₁=2.100 676×10^－4 m².再由式(29)、式(31) 和式(34) 可计算出此时的f=0.523 0.如图 4所示, 当P=100 kN时的试验数据为f=0.466 93.从图 4a可知, 静摩擦因数与试验数据的变化趋势一致, 均随总法向载荷呈增长趋势.在图 4b和图 4c中, 分别给出了理论计算值与试验值的绝对和相对误差.绝对误差和相对误差的变化范围分别为-0.139 6~0.056 07和-38.78%~12.01%.另外, 试验值的最大与最小差值:f_max－f_min=0.466 93－0.36=0.106 93, 整体波动范围较小, 因此分析两者之间的绝对误差较为合理.基于以上分析, 本文模型对于工程结合面上静摩擦因数预测有一定的合理性.

6 结论

1) 静摩擦因数与分形维数D、分形尺度系数G有着较复杂的非线性关系.

2) 静摩擦因数随D的增大先增大后减小并以D=2.5为界; 静摩擦因数随G的减小而增大.

3) 对比本文理论模型计算值和试验值, 绝对误差范围为-0.139 6~0.056 07.