考虑微凸体相互作用及硬度变化的加卸载接触分形模型

Fractal Model of Loading-Unloading Contact Considering Asperity Interaction and Hardness Changes

摘要: 基于分形理论，建立了考虑微凸体的相互作用及硬度变化的加卸载接触模型。推导了单个微凸体在加卸载过程中接触面积、接触载荷、表面分离、变形量的关系式。根据接触面积的分布密度函数，推导粗糙表面的实际接触特性关系式，并通过数值计算分析了相关参数对接触特性的影响。结果表明：单个微凸体发生弹性变形时，卸载过程与加载过程完全一致；发生弹塑性变形时，随着微凸体加载结束时的变形量的增大，在相同接触载荷下，卸载时的接触面积越大，且不能完全回复；发生完全塑性变形时，微凸体在卸载过程中不回复。随着分形维度D的增大、尺度系数G的减小，粗糙表面越光滑，在相同接触载荷下，接触面积越大，表面分离减小。

Abstract: Based on fractal theory, a loading-unloading contact model considering the interaction of asperities and the change of hardness is established. The relationship between contact area, contact load, sur-face separation and deformation of a single asperity during loading and unloading is derived. Ac-cording to the distribution density function of the contact area, the actual contact characteristics of rough surfaces are derived, and the influence of relevant parameters on the contact characteristics is analyzed through numerical calculation. The results show that the unloading process is consistent with the loading process when a single asperity undergoes elastic deformation; when elastic-plastic deformation occurs, with the increase of deformation at the end of loading, the contact area at un-loading is larger under the same contact load, and cannot be fully recovered; when complete plastic deformation occurs, the asperities do not recover during unloading. With the increase of fractal di-mension D and fractal roughness G, the smoother the rough surface is, the larger the contact area is, and the smaller the surface separation is under the same contact load.

1. 引言

1882年Hertz [1] 提出的弹性接触理论是经典接触力学的基础。Greenwood和Williamson [2] 基于Hertz理论，将两个粗糙表面接触等效为粗糙表面和刚性光滑表面间的接触，提出了GW模型。Chang等 [3] 基于GW模型和体积守恒定律，考虑到塑性变形，提出了CEB模型。但CEB模型没有考虑弹性变形和塑性变形之间的过渡状态，即弹塑性变形状态，在临界屈服点处发生跳跃式变化。Zhao等 [4] 将微凸体的变形分为弹性变形、弹塑性变形、塑性变形三个阶段，并利用对数函数和四次多项式法提出了ZMC模型。Kought和Estion [5] 利用有限元方法，对半球体和粗糙表面间的弹塑性接触进行分析，提出了KE模型。杨红平等 [6] 基于ZMC模型，采取统计参数、分形参数结合的方法建立了一种新的接触刚度模型。赵永武等 [7] 建立结合面法向微观接触模型，符合压力的平滑性、连续性的需求。田红亮等 [8] 基于Hertz理论和GW模型，建立了适用于较大接触载荷下的法向接触模型。肖会芳等 [9] 利用椭圆曲线方程，建立了结合面法向接触刚度模型。

国内外学者在研究粗糙表面形貌时，发现粗糙表面具有连续性、不可微性和自仿射性等分形特征。因此，分形接触模型也成为了一种研究结合面接触问题的有效方法。Majumdar和Bhushan [10] 基于与尺度无关的粗糙表面分形参数和Weierstrass-Mandelbrot函数，建立了MB模型。模型认为随着载荷增大，微凸体由塑性变形转变为弹性变形。Morag和Etsion [11] 基于MB模型，对弹塑性变形临界高度进行修改，提出了ME模型。对MB模型中与传统的接触研究相反的问题进行解答，并证明微凸体临界接触面积的尺度相关性。Wang和Komvopoulos [12] 考虑了接触面积分布域扩展因子的作用，对MB模型进行了修正。葛世荣等 [13] 研究发现经过不同加工工艺的表面轮廓具有鲜明的分形特征。陈国安和葛世荣 [14] 提出了“特征粗糙度”的概念，用以描述粗糙度表面轮廓。李小彭等 [15] 基于静摩擦因数模型 [16]，建立了域扩展因子是分形维数函数的静摩擦因数分形模型。

本文基于分形理论，建立了综合考虑微凸体间的相互作用及弹塑性变形时硬度随深度变化而变化的粗糙表面接触特性模型。

2. 粗糙表面的加卸载接触特性模型

2.1. 粗糙表面的分形表征

两个粗糙表面的接触问题可以等效为刚性光滑平面和粗糙表面的接触问题。粗糙表面的表面形貌具有无序、随机和多尺度的特性，其二维轮廓可以WM函数进行表征，表达式如下：

$Z\left(x\right)=G^{D-1}{\displaystyle {\sum }_{n=n_{\min }}^{\infty }\frac{\cos \left(2\pi {\gamma }^{n}x\right)}{{\gamma }^{\left(2-D\right)n}}}$ (1)

式中： $Z\left(x\right)$ 为表面轮廓曲线在x处的高度；G为轮廓分形特征尺度系数；D为分形维数； $\gamma$ 为频率系数，轮廓高度服从正态分布时， $\gamma =1.5$ ；n为频率指数。

根据式(1)，微凸体变形前的轮廓曲线可描述为：

$z\left(x\right)=G^{D-1}{l}^{2-D}\cos \left(\frac{2\pi x}{l}\right)$ (2)

式中： $l=\frac{1}{{\gamma }^n}$ ；

微凸体顶端的曲率半径的表达式为：

$R=\frac{l^D}{{\pi }^2{G}^{D-1}}$ (3)

2.2. 微凸体间相互作用的机理

微凸体间的相互作用是指微凸体发生变形时，对周围微凸体变形造成的影响。图1为等效后的接触面在相互作用下的示意图。

Figure 1. Schematic of the contact interface

图1. 接触面示意图

如图1所示，z表示未变形前的微凸体的高度， ${\omega }_a$ 表示系统总变形， $u_m$ 表示微凸体平均平面的位移， $u_g$ 表示微凸体相互作用引起的变形，d表示刚性光滑平面与微凸体平均平面之间的距离。从图1中的几何关系可知：

$\omega =z-d+\left(u_m-u_g\right)$ (4)

研究表明 [17]：

$u_m-u_g=\frac{1.12}{E}\sqrt{f{P}_m}$ (5)

式中：E为等效弹性模量， $E=\left(\frac{1-{\nu }_1^2}{E_1}+\frac{1-{\nu }_2^2}{E_2}\right)$， $E_1$ 、 ${\nu }_1$ 、 $E_2$ 、 ${\nu }_2$ 分别为粗糙表面的弹性模量和泊松比；f为微凸体所承受的载荷； $P_m$ 为两表面间的接触压强。

将式(4)代入式(3)得：

$\omega =z-d+\frac{1.12}{E}\sqrt{f{P}_m}$ (6)

2.3. 单个微凸体的接触特性模型

2.3.1. 加载过程

在微凸体加载过程中，随着接触载荷的增大，微凸体可能发生弹性、弹塑性和塑性变形。根据赫兹接触理论 [1]，当 $\omega \le {\omega }_{ec}$ 时，微凸体发生弹性变形，接触面积和接触载荷的表达式为：

$a_e=\pi R\omega =\frac{l^D}{\pi G^{D-1}}\omega$ (7)

$f_e=\frac{4}{3}E{R}^{1/2}{\omega }^{3/2}=\frac{4E{\pi }^{1/2}{G}^{D-1}}{3l^D}{a}_e^{3/2}$ (8)

式中：K为硬度系数， $K=0.454+0.41\nu$ ；H为较软材料的硬度。

根据文献 [18] 可得：

${\omega }_{ec}={\left(\frac{3\pi KH}{4E}\right)}^{2}R$ (9)

将式(8)代入式(6)中，可得考虑微凸体相互作用的变形量 $\omega$ ：

$\omega =z-d+\frac{1.12}{E}\sqrt{\frac{4}{3}E{R}^{1/2}{\omega }^{3/2}{P}_m}$ (10)

将式(9)代入式(7)中，可得弹性临界面积 $a_{ec}$ ：

$a_e=\pi R{\omega }_{ec}=\frac{1}{\pi }{\left(\frac{3KH{l}^D}{4E{G}^{D-1}}\right)}^2$ (11)

当 ${\omega }_{ec}<\omega \le 6{\omega }_{ec}$ 时，微凸体处于第一弹塑性变形阶段，加载过程中的接触面积和接触载荷 [5] 的表达式为：

$a_{ep1}=0.93{\left(\frac{\omega }{{\omega }_{ec}}\right)}^{1.136}{a}_{ec}$ (12)

$f_{ep1}=KH{\left(\frac{a_{ep1}}{a_{ec}}\right)}^{c_2}{a}_{ep1}=KH{a}_{ec}^{-c_2}{a}_{ep1}^{1+c_2}$ (13)

式中： $c_1=2.8K$ ； $c_2=\frac{\ln \left(1.1218\times {7.1197}^{0.2544}\right)}{\ln 7.1197}$。

将式(13)代入式(6)中，可得考虑微凸体相互作用的变形量 $\omega$ ：

$\omega =z-d+\frac{1.12}{E}\sqrt{{0.93}^{1+c_2}KH{a}_{ec}{\left(\frac{\omega }{{\omega }_{ec}}\right)}^{1.136\left(1+c_2\right)}{P}_m}$ (14)

当 $6{\omega }_{ec}<\omega \le 110{\omega }_{ec}$ 时，微凸体处于第二弹塑性变形阶段，加载过程中的接触面积和接触载荷的表达式为：

$a_{ep2}=0.94{\left(\frac{\omega }{{\omega }_{ec}}\right)}^{1.146}{a}_{ec}$ (15)

$f_{ep2}=c_{3}Y{\left(\frac{a_{ep2}}{a_{ec}}\right)}^{c_4}{a}_{ep2}=c_{3}Y{a}_{ec}^{-c_4}{a}_{ep1}^{1+c_4}$ (16)

式中： $c_3=K\times 3.15896\times {7.1197}^{0.2544-c_4}$ ； $c_2=\frac{\ln 1.1218+0.2544\ln 7.1197-\ln 1.4988-0.1021\ln 205.3827}{\ln 7.1197}$。

将式(16)代入式(6)中，可得考虑微凸体相互作用的变形量 $\omega$ ：

$\omega =z-d+\frac{1.12}{E}\sqrt{{0.94}^{1+c_4}{c}_{3}Y{a}_{ec}{\left(\frac{\omega }{{\omega }_{ec}}\right)}^{1.146\left(1+c_2\right)}{P}_m}$ (17)

当 $110{\omega }_{ec}<\omega$ 时，微凸体处于完全塑性变形阶段，加载过程中的接触面积和接触载荷的表达式为：

$a_p=2\pi R\omega =\frac{2l^D}{\pi G^{D-1}}\omega$ (18)

$f_p=H{a}_p$ (19)

将式(19)代入式(6)中，可得考虑微凸体相互作用的变形量 $\omega$ ：

$\omega =z-d+\frac{1.12}{E}\sqrt{\frac{2H{l}^D}{\pi G^{D-1}}\omega P_m}$ (20)

2.3.2. 卸载过程

图2为微凸体的卸载示意图。图中，R为微凸体发生变形前顶端的曲率半径； $R^u$ 为完全卸载后微凸体顶端的曲率半径； ${\omega }_{\max }$ 为微凸体完全加载时的变形量； ${\omega }_{res}$ 为微凸体完全卸载后的残余变形量。

假设微凸体相互作用引起的变形为弹性变形，在卸载过程中完全回复，则由式(6)可得：

${\omega }^u=z-d+\frac{1.12}{E}\sqrt{f^u{P}_m}$ (21)

式中： $f^u$ 为卸载过程中微凸体所承受的载荷。

当 ${\omega }_{\max }\le {\omega }_{ec}$ 时，微凸体只发生弹性变形，在卸载过程中，微凸体完全回复。完全卸载后，微凸体的高度和顶端曲率半径没有变化，卸载过程中的接触面积和接触载荷的表达式为：

$R^u=R=\frac{l^D}{\pi G^{D-1}}$ (22)

$a_e^u=\pi R{\omega }^u=\frac{l^D}{\pi G^{D-1}}{\omega }^u$ (23)

$f_e^u=\frac{4}{3}E{\left(R^u\right)}^{1/2}{\left({\omega }^u\right)}^{3/2}=\frac{4E{\pi }^{1/2}{G}^{D-1}}{3l^D}{\left(a_e^u\right)}^{3/2}$ (24)

${\omega }_{res}=0$ (25)

Figure 2. Loading and unloading diagram of a single asperity

图2. 单个微凸体加卸载示意图

将式(24)代入式(21)中，可得卸载过程中，考虑微凸体相互作用的变形量 ${\omega }^u$ ：

${\omega }^u=z-d+\frac{1.12}{E}\sqrt{\frac{4}{3}E{\left(R^u\right)}^{1/2}{\left({\omega }^u\right)}^{3/2}{P}_m}$ (26)

当 ${\omega }_{ec}<{\omega }_{\max }\le 110{\omega }_{ec}$ 时，微凸体发生弹塑性变形，由于存在塑性变形，微凸体不能完全回复。Etsion [19] 等研究表明，加载过程中的最大变形量 ${\omega }_{\max }$ 与残余变形量 ${\omega }_{res}$ 的关系表达式为：

$\frac{{\omega }_{res}}{{\omega }_{\max }}=\left(1-\frac{1}{{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.28}}\right)\left(1-\frac{1}{{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.69}}\right)$ (27)

微凸体变形前的曲率半径R与卸载后的曲率半径 $R^u$ 的表达式如下：

$\frac{R^u}{R}=1+1.275{\left(\frac{E}{Y}\right)}^{-0.216}\left(\frac{{\omega }_{\max }}{{\omega }_{ec}}-1\right)$ (28)

微凸体在卸载过程中，接触载荷 $f^u$ 、接触面积 $a^u$ 与变形量 ${\omega }^u$ 的表达式如下：

$f^u=f_{\max }{\left(\frac{{\omega }^u-{\omega }_{res}}{{\omega }_{\max }-{\omega }_{res}}\right)}^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{-0.0331}}$ (29)

$a^u=a_{\max }{\left(\frac{{\omega }^u-{\omega }_{res}}{{\omega }_{\max }-{\omega }_{res}}\right)}^{{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{-0.12}}$ (30)

将式(30)代入式(29)中，可得接触载荷 $f^u$ 与接触面积 $a^u$ 的表达式：

$f^u=f_{\max }{\left(\frac{a^u}{a_{\max }}\right)}^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}$ (31)

将式(29)代入式(21)中，可得卸载过程中，考虑微凸体相互作用的变形量 ${\omega }^u$ ：

${\omega }^u=z-d+\frac{1.12}{E}\sqrt{f_{\max }{\left(\frac{{\omega }^u-{\omega }_{res}}{{\omega }_{\max }-{\omega }_{res}}\right)}^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{-0.0331}}{P}_m}$ (32)

当 ${\omega }_{ec}<{\omega }_{\max }\le 6{\omega }_{ec}$ 时，微凸体发生第一弹塑性变形，变形量 ${\omega }^u$ 与接触载荷 $f^u$ 、接触面积 $a^u$ 的表达式如下：

$a_{ep1}^u=0.93{\left(\frac{{\omega }_{\max }}{{\omega }_{ec}}\right)}^{1.136}{\left(\frac{{\omega }^u-{\omega }_{res}}{{\omega }_{\max }-{\omega }_{res}}\right)}^{{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{-0.12}}{a}_{ec}$ (33)

$f_{ep1}^u=KH{\left[0.93{\left(\frac{{\omega }_{\max }}{{\omega }_{ec}}\right)}^{1.136}\right]}^{1+c_2}{\left(\frac{a_{ep1}^u}{0.93{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{1.136\left(1+c_2\right)}}\right)}^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}{a}_{ec}$ (34)

${\omega }^u=z-d+\frac{1.12}{E}\sqrt{{0.93}^{1+c_2}KH{a}_{ec}{\left(\frac{{\omega }_{\max }}{{\omega }_{ec}}\right)}^{1.136\left(1+c_2\right)}{\left(\frac{{\omega }^u-{\omega }_{res}}{{\omega }_{\max }-{\omega }_{res}}\right)}^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{-0.0331}}{P}_m}$ (35)

当 $6{\omega }_{ec}<{\omega }_{\max }\le 110{\omega }_{ec}$ 时，微凸体发生第二弹塑性变形，变形量 ${\omega }^u$ 与接触载荷 $f^u$ 、接触面积 $a^u$ 的表达式如下：

$a_{ep2}^u=0.94{\left(\frac{{\omega }_{\max }}{{\omega }_{ec}}\right)}^{1.146}{\left(\frac{{\omega }^u-{\omega }_{res}}{{\omega }_{\max }-{\omega }_{res}}\right)}^{{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{-0.12}}{a}_{ec}$ (36)

$f_{ep1}^u=c_{3}Y{\left[0.94{\left(\frac{{\omega }_{\max }}{{\omega }_{ec}}\right)}^{1.146}\right]}^{1+c_4}{\left(\frac{a_{ep2}^u}{0.94{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{1.146\left(1+c_4\right)}}\right)}^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}{a}_{ec}$ (37)

${\omega }^u=z-d+\frac{1.12}{E}\sqrt{{0.94}^{1+c_4}{c}_{3}Y{a}_{ec}{\left(\frac{{\omega }_{\max }}{{\omega }_{ec}}\right)}^{1.146\left(1+c_4\right)}{\left(\frac{{\omega }^u-{\omega }_{res}}{{\omega }_{\max }-{\omega }_{res}}\right)}^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{-0.0331}}{P}_m}$ (38)

当 $110{\omega }_{ec}<{\omega }_{\max }$ 时，微凸体发生完全塑性变形，在卸载过程中不回复，卸载过程中的接触载荷和接触面积为0。

2.4. 粗糙表面的实际接触面积和接触载荷

2.4.1. 加载过程

根据文献 [10]，加载过程中接触面积的分布密度函数 $n\left(a\right)$ 为：

$n\left(a\right)=\frac{D}{2}{a}_l^{D/2}{a}^{-\left(D+2\right)/2}$ (39)

式中： $a_l$ 为加载过程中，所有微凸体的最大接触面积。

实际接触面积 $A_r$ 为：

$A_r={\displaystyle {\int }_0^{a_l}n\left(a\right)a\text{d}a}=\frac{D}{2-D}{a}_l$ (40)

当 $a_l\le a_{ec}$ 时，接触面上所有接触的微凸体处于弹性变形阶段，实际接触面积 $A_r$ 和接触载荷 $F_r$ 为：

$A_r=A_{re}={\displaystyle {\int }_0^{a_l}n\left(a\right)a\text{d}a}=\frac{D}{2-D}{a}_l$ (41)

$F_r=F_{re}={\displaystyle {\int }_0^{a_l}{f}_{e}n\left(a\right)a\text{d}a}=\frac{D}{3-D}\frac{4E{\pi }^{1/2}{G}^{D-1}}{3l^D}{a}_l^{3/2}$ (42)

当 $a_{ec}<a_l\le 7.1197a_{ec}$ 时，接触面上所有接触的微凸体处于弹性或第一弹塑性变形阶段，实际接触面积 $A_r$ 和接触载荷 $F_r$ 为：

$A_r=A_{re}+A_{rep1}=\frac{D}{2-D}{a}_l$ (43)

$\begin{array}{c}{F}_r=F_{re}+F_{rep1}\\ F_{re}={\displaystyle {\int }_0^{a_{ec}}{f}_{e}n\left(a\right)a\text{d}a}=\frac{D}{3-D}\frac{4E{\pi }^{1/2}{G}^{D-1}}{3l^D}{a}_{ec}^{\left(3-D\right)/2}{a}_l^{D/2}\\ F_{rep1}={\displaystyle {\int }_{a_{ec}}^{a_l}{f}_{ep1}n\left(a\right)a\text{d}a}=\frac{KHD}{2c_2-D+2}\left(a_{ec}^{-c_2}{a}_l^{1+c_2}-a_{ec}^{\left(2-D\right)/2}{a}_l^{D/2}\right)\end{array}\}$ (44)

式中： $A_{re}$ 、 $F_{re}$ 为加载过程中发生弹性变形的微凸体的实际接触面积和接触载荷； $A_{rep1}$ 、 $F_{rep1}$ 为加载过程中发生第一弹塑性变形的微凸体的实际接触面积和接触载荷。

当 $7.1197a_{ec}<a_l\le 205.3827a_{ec}$ 时，接触面上所有接触的微凸体处于弹性、第一弹塑性或第二弹塑性变形阶段，实际接触面积 $A_r$ 和接触载荷 $F_r$ 为：

$A_r=A_{re}+A_{rep1}+A_{rep2}=\frac{D}{2-D}{a}_l$ (45)

$\begin{array}{c}{F}_r=F_{re}+F_{rep1}+F_{rep2}\\ F_{re}={\displaystyle {\int }_0^{a_{ec}}{f}_{e}n\left(a\right)a\text{d}a}=\frac{D}{3-D}\frac{4E{\pi }^{1/2}{G}^{D-1}}{3l^D}{a}_{ec}^{\left(3-D\right)/2}{a}_l^{D/2}\\ F_{rep1}={\displaystyle {\int }_{a_{ec}}^{a_{epc}}{f}_{ep1}n\left(a\right)a\text{d}a}=\frac{KHD}{2c_2-D+2}\left({7.1197}^{\left(2+2c_2-D/2\right)}-1\right)a_{ec}^{\left(2-D\right)/2}{a}_l^{D/2}\\ F_{rep2}={\displaystyle {\int }_{a_{epc}}^{a_l}{f}_{ep2}n\left(a\right)a\text{d}a}=\frac{c_{3}YD}{2c_4-D+2}\left(a_{ec}^{-c_4}{a}_l^{1+c_4}-{7.1197}^{\left(2+2c_2-D/2\right)}{a}_{ec}^{\left(2-D\right)/2}{a}_l^{D/2}\right)\end{array}\}$ (46)

当 $205.3827a_{ec}<a_l$ 时，接触面上所有接触的微凸体处于弹性、弹塑性或完全塑性变形阶段，实际接触面积 $A_r$ 和接触载荷 $F_r$ 为：

$A_r=A_{re}+A_{rep1}+A_{rep2}+A_p=\frac{D}{2-D}{a}_l$ (47)

$\begin{array}{c}{F}_r=F_{re}+F_{rep1}+F_{rep2}+F_{rp}\\ F_{re}={\displaystyle {\int }_0^{a_{ec}}{f}_{e}n\left(a\right)a\text{d}a}=\frac{D}{3-D}\frac{4E{\pi }^{1/2}{G}^{D-1}}{3l^D}{a}_{ec}^{\left(3-D\right)/2}{a}_l^{D/2}\\ F_{rep1}={\displaystyle {\int }_{a_{ec}}^{a_{epc}}{f}_{ep1}n\left(a\right)a\text{d}a}=\frac{KHD}{2c_2-D+2}\left({7.1197}^{\left(2+2c_2-D/2\right)}-1\right)a_{ec}^{\left(2-D\right)/2}{a}_l^{D/2}\\ F_{rep2}={\displaystyle {\int }_{a_{epc}}^{a_{pc}}{f}_{ep2}n\left(a\right)a\text{d}a}=\frac{c_{3}YD}{2c_4-D+2}\left({205.3827}^{\left(2+2c_4-D/2\right)}-{7.1197}^{\left(2+2c_2-D/2\right)}\right)a_{ec}^{\left(2-D\right)/2}{a}_l^{D/2}\\ F_{rp}={\displaystyle {\int }_{a_{pc}}^{a_l}{f}_{p}n\left(a\right)a\text{d}a}=\frac{HD}{2-D}\left(a_l^{\left(2-D\right)/2}-{205.3827}^{\left(2+2c_4-D/2\right)}{a}_{ec}^{\left(2-D\right)/2}\right)a_l^{D/2}\end{array}\}$ (48)

2.4.2. 卸载过程

假设在卸载过程中，微凸体的接触面积分布形式仍满足分形特征，分布函数为：

$n^u\left(a^u\right)=H\frac{D}{2}{\left(a_l^u\right)}^{D/2}{\left(a^u\right)}^{-\left(D+2\right)/2}$ (49)

式中： $a_l^u$ 为卸载过程中，所有微凸体的最大接触面积；H为接触面积分布函数的修正系数。

当 $a_l\le a_{ec}$ 时，在接触面开始卸载时，所有微凸体处于弹性变形阶段。接触面的接触面积和接触载荷在加载完成时和卸载开始时相等，即：

$\begin{array}{c}{\displaystyle {\int }_0^{a_l}n\left(a\right)a\text{d}a}={\displaystyle {\int }_0^{a_l}{n}_e^u\left(a^u\right)a^u\text{d}{a}^u}\\ {\displaystyle {\int }_0^{a_l}{f}_{e}n\left(a\right)\text{d}a}={\displaystyle {\int }_0^{a_l}{f}_e^u{n}_e^u\left(a^u\right)\text{d}{a}^u}\end{array}\}$ (50)

可得修正指数 $H_e$ ：

$H_e=1$ (51)

则卸载过程中的实际接触面积 $A_r^u$ 和接触载荷 $F_r^u$ 为：

$\begin{array}{c}{A}_r^u=A_{re}^u={\displaystyle {\int }_0^{a_l}{n}_e^u\left(a^u\right)a^u\text{d}{a}^u}\\ =\frac{D}{2-D}{H}_e{a}_l^u\end{array}$ (52)

$\begin{array}{c}{F}_r^u=F_{re}^u={\displaystyle {\int }_0^{a_l}{f}_e^u{n}_e^u\left(a^u\right)\text{d}{a}^u}\\ =\frac{D}{3-D}\frac{4E{\pi }^{1/2}{G}^{D-1}}{3l^D}{H}_e{\left(a_l^u\right)}^{3/2}\end{array}$ (53)

当 $a_{ec}<a_l\le 7.1197a_{ec}$ 时，在接触面开始卸载时，所有接触的微凸体处于弹性或第一弹塑性变形阶段。接触面的各部分接触面积在加载完成时和卸载开始时相等，即：

$\begin{array}{c}{\displaystyle {\int }_0^{a_{ec}}n\left(a\right)a\text{d}a}={\displaystyle {\int }_0^{a_{ec}}{n}_e^u\left(a^u\right)a^u\text{d}{a}^u}\\ {\displaystyle {\int }_{a_{ec}}^{a_l}n\left(a\right)a\text{d}a}={\displaystyle {\int }_0^{a_l}{n}_{ep1}^u\left(a^u\right)a^u\text{d}{a}^u}\end{array}\}$ (54)

可得修正指数 $H_e$ 、 $H_{ep1}$ ：

$\begin{array}{c}{H}_e=1\\ H_{ep1}=1-{\left(\frac{a_{ec}}{a_l}\right)}^{\left(2-D\right)/2}\end{array}\}$ (55)

可得卸载过程中的实际接触面积 $A_r^u$ 为：

$\begin{array}{c}{A}_r^u=A_{re}^u+A_{rep1}^u\\ A_{re}^u={\displaystyle {\int }_0^{a_{el}}{n}_e^u\left(a^u\right)a^u\text{d}{a}^u}=\frac{D}{2-D}{H}_e{\left(a_{el}^u\right)}^{\left(2-D\right)/2}{\left(a_l^u\right)}^{D/2}\\ A_{rep1}^u={\displaystyle {\int }_0^{a_l}{n}_{ep1}^u\left(a^u\right)a^u\text{d}{a}^u}=\frac{D}{2-D}{H}_{ep1}{a}_l^u\end{array}\}$ (56)

接触面的各部分接触载荷在加载完成时和卸载开始时相等，但由于式(27)、(28)为有限元拟合计算的，存在一定误差。因此，陈建江等 [20] 研究表明，可引入载荷误差修正系数K避免载荷关系式不等的问题。卸载过程中的实际接触面积 $F_r^u$ ：

$\begin{array}{c}{F}_r^u=F_{re}^u+F_{rep1}^u\\ F_{re}^u={\displaystyle {\int }_0^{a_{el}^u}{K}_e{f}_e^u{n}_e^u\left(a^u\right)\text{d}{a}^u}=\frac{D}{3-D}\frac{4E{\pi }^{1/2}{G}^{D-1}}{3l^D}{K}_e{H}_e{\left(a_{el}^u\right)}^{\left(3-D\right)/2}{\left(a_l^u\right)}^{D/2}\\ \begin{array}{l}{F}_{rep1}^u={\displaystyle {\int }_0^{a_l^u}{K}_{ep1}{f}_{ep1}^u{n}_{ep1}^u\left(a^u\right)\text{d}{a}^u}\\ =\frac{KHD}{2}{K}_{ep1}{H}_{ep1}\cdot {0.93}^{1+c_2-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}{\left[1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}-\frac{D}{2}\right]}^{-1}\\ \cdot {\left({\omega }_{\max }/{\omega }_{ec}\right)}^{1.136\left[1+c_2-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}\right]}{\left(a_{ec}\right)}^{1-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}{\left(a_l^u\right)}^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}\end{array}\end{array}\}$ (57)

接触面的各部分接触载荷在加载完成时和卸载开始时相等，即：

$\begin{array}{c}{\displaystyle {\int }_0^{a_{ec}}{f}_{e}n\left(a\right)a\text{d}a}={\displaystyle {\int }_0^{a_{ec}}{K}_e{f}_e^u{n}_e^u\left(a^u\right)\text{d}{a}^u}\\ {\displaystyle {\int }_{a_{ec}}^{a_l}{f}_{ep1}n\left(a\right)a\text{d}a}={\displaystyle {\int }_0^{a_l}{K}_{ep1}{f}_{ep1}^u{n}_{ep1}^u\left(a^u\right)\text{d}{a}^u}\end{array}\}$ (58)

可得载荷误差修正系数 $K_e$ 、 $K_{ep1}$ 为：

$\begin{array}{c}{K}_e=1\\ \begin{array}{l}{K}_{ep1}=\frac{a_{ec}^{-c_2}{a}_l^{1+c_2}-a_{ec}^{\left(2-D\right)/2}{a}_l^{D/2}}{\left(c_2-\frac{D}{2}-1\right)H_{ep1}\cdot {0.93}^{1+c_2-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}{\left[1+c_2-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}-\frac{D}{2}\right]}^{-1}}\\ \cdot \frac{1}{{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{1.136\left[1+c_2-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}\right]}{\left(a_{ec}\right)}^{1-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}{\left(a_l\right)}^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}}\end{array}\end{array}\}$ (59)

当 $7.1197a_{ec}<a_l\le 205.3827a_{ec}$ 时，在接触面开始卸载时，所有接触的微凸体处于弹性、第一弹塑性或第二弹塑性变形阶段，卸载过程中的实际接触面积 $A_r^u$ 为：

$\begin{array}{c}{A}_r^u=A_{re}^u+A_{rep1}^u+A_{rep2}^u\\ A_{re}^u={\displaystyle {\int }_0^{a_{el}}{n}_e^u\left(a^u\right)a^u\text{d}{a}^u}=\frac{D}{2-D}{H}_e{\left(a_{el}^u\right)}^{\left(2-D\right)/2}{\left(a_l^u\right)}^{D/2}\\ A_{rep1}^u={\displaystyle {\int }_0^{a_{ep1l}}{n}_{ep1}^u\left(a^u\right)a^u\text{d}{a}^u}=\frac{D}{2-D}{H}_{ep1}{\left(a_{ep1l}^u\right)}^{\left(2-D\right)/2}{\left(a_l^u\right)}^{D/2}\\ A_{rep2}^u={\displaystyle {\int }_0^{a_l}{n}_{ep2}^u\left(a^u\right)a^u\text{d}{a}^u}=\frac{D}{2-D}{H}_{ep2}{a}_l^u\end{array}\}$ (60)

接触面的各部分接触面积在加载完成时和卸载开始时相等，可得修正指数 $H_e$ 、 $H_{ep1}$ 、 $H_{ep2}$ ：

$\begin{array}{c}{H}_e=1\\ H_{ep1}=1-{\left(\frac{1}{7.1197}\right)}^{\left(2-D\right)/2}\\ H_{ep1}=1-{\left(\frac{7.1197a_{ec}}{a_l}\right)}^{\left(2-D\right)/2}\end{array}\}$ (61)

卸载过程中的实际接触面积 $F_r^u$ 为：

$\begin{array}{c}{F}_r^u=F_{re}^u+F_{rep1}^u+F_{rep2}^u\\ F_{re}^u={\displaystyle {\int }_0^{a_{el}^u}{K}_e{f}_e^u{n}_e^u\left(a^u\right)\text{d}{a}^u}=\frac{D}{3-D}\frac{4E{\pi }^{1/2}{G}^{D-1}}{3l^D}{K}_e{H}_e{\left(a_{el}^u\right)}^{\left(3-D\right)/2}{\left(a_l^u\right)}^{D/2}\\ \begin{array}{l}{F}_{rep1}^u={\displaystyle {\int }_0^{a_{ep1l}^u}{K}_{ep1}{f}_{ep1}^u{n}_{ep1}^u\left(a^u\right)\text{d}{a}^u}\\ =\frac{KHD}{2}{K}_{ep1}{H}_{ep1}\cdot {7.1197}^{c_2-0.7527}{\left[1.7527-\frac{D}{2}\right]}^{-1}{\left(a_{ec}\right)}^{-0.7527}{\left(a_{ep1l}^u\right)}^{1.7527-\frac{D}{2}}{\left(a_l^u\right)}^{D/2}\end{array}\\ \begin{array}{l}{F}_{rep2}^u={\displaystyle {\int }_0^{a_l^u}{K}_{ep2}{f}_{ep2}^u{n}_{ep2}^u\left(a^u\right)\text{d}{a}^u}\\ =\frac{c_{3}YD}{2}{K}_{ep2}{H}_{ep2}\cdot {0.94}^{1+c_4-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}{\left[1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}-\frac{D}{2}\right]}^{-1}\\ \cdot {\left({\omega }_{\max }/{\omega }_{ec}\right)}^{1.146\left[1+c_4-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}\right]}{\left(a_{ec}\right)}^{1-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}{\left(a_l^u\right)}^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}\end{array}\end{array}\}$ (62)

可得载荷误差修正系数 $K_e$ 、 $K_{ep1}$ 、 $K_{ep2}$ 为：

$\begin{array}{c}{K}_e=1\\ K_{ep1}=\frac{{7.1197}^{1+c_2-\frac{D}{2}}-1}{\left(c_2-\frac{D}{2}+1\right)H_{ep1}\cdot {7.1197}^{1+c_2-\frac{D}{2}}{\left[1.7527-\frac{D}{2}\right]}^{-1}}\\ \begin{array}{l}{K}_{ep2}=\frac{a_{ec}^{-c_4}{a}_l^{1+c4}-{7.1197}^{1+c_4-\frac{D}{2}}{a}_{ec}^{\left(2-D\right)/2}{a}_l^{D/2}}{\left(c_4-\frac{D}{2}+1\right)H_{ep2}\cdot {0.93}^{1+c_4-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}{\left[1+c_4-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}-\frac{D}{2}\right]}^{-1}}\\ \cdot \frac{1}{{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{1.146\left[1+c_4-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}\right]}{\left(a_{ec}\right)}^{1-1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}{\left(a_l\right)}^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}}\end{array}\end{array}\}$ (63)

当 $205.3827a_{ec}<a_l$ 时，在接触面加载结束时，所有接触的微凸体处于弹性、弹塑性或塑性变形阶段，由于处于塑性阶段的微凸体在卸载过程中不回复，则 $a_l$ 可视为塑性临界接触面积 $205.3827a_{ec}$，并代入式(60)、(61)、(62)、(63)中，可得对应表达式。

3. 数值分析

通过对粗糙面接触模型进行数值计算，分析接触面在加卸载过程中的接触特性演变规律。为了方便数据处理，对以下变量进行无量纲化处理： $F^{*}=F/A_{a}E$， $A^{*}=A/A_a$， $d^{*}=d/l$。其中，F为加卸载过程的实际接触载荷，A为加卸载过程的实际接触面积， $A_a$ 为名义接触面积。

下面是基于本文推导的接触模型，计算接触面的接触特性。相关参数参考文献 [21]：分形维度 $D=1.4$，尺度系数 $G=1.36\times {10}^{-10}\text{m}$，较软材料的泊松比 $\nu =0.17$，等效弹性模量 $E=72\text{GPa}$，较软材料的硬度 $H=5.5\text{GPa}$。

图3为单个微凸体在加卸载过程中无量纲接触面积和接触载荷的关系。在卸载点1卸载时，微凸体只发生弹性变形，在卸载过程中完全恢复，所以加卸载过程中的接触特性变化曲线完全一致，接触载荷和接触面积的关系为 $f~a^{3/2}$。在卸载点2卸载时，微凸体处于第一弹塑性变形阶段，在加载过程中该阶段的接触载荷和接触面积的关系为 $f~a^{1+c_2}$ ；在卸载过程中，由于发生弹塑性变形，微凸体不能完全回复，二者曲线不重合，接触载荷和接触面积的关系为 $f~a^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}$。在卸载点3卸载时，微凸体处于第二弹塑性变形阶段，在加载过程中该阶段的接触载荷和接触面积的关系为 $f~a^{1+c_4}$ ；在卸载过程中，接触载荷和接触面积的关系为 $f~a^{1.5{\left({\omega }_{\max }/{\omega }_{ec}\right)}^{0.0869}}$。从卸载点2、3处卸载的关系曲线可以看出，在相同接触载荷下，卸载时的接触面积大于加载时的接触面积。这是因为由于发生弹塑性变形，微凸体的曲率半径增大，从而导致在相同载荷下接触面积增大，且微凸体的变化量越大，卸载时的接触面积越大。

Figure 3. Relationship between dimensionless contact area and load of a single asperity

图3. 单个微凸体无量纲接触面积和接触载荷的关系

图4为粗糙面相接触时，在不同分形维度D的条件下，无量纲接触载荷与无量纲接触面积、无量纲两表面间距离的关系。在尺度系数 $G=1.36\times {10}^{-10}\text{m}$ 与相同载荷的条件下，随着分形维度D的增加，微凸体的、接触面积增大、表面间的距离减小。

(a) 接触载荷和表面分离 (b) 接触载荷和接触面积

Figure 4. The relationship under the change of fractal dimension D

图4. 分形维度D的变化下的关系

图5为粗糙面相接触时，在不同轮廓分形特征尺度系数G的条件下，无量纲接触载荷与无量纲接触面积、无量纲两表面间距离的关系。在分形维度 $D=1.4$ 与相同载荷的条件下，随着尺度系数G的增加，微凸体的接触面积减小、表面间的距离增大。说明随着尺度系数G的减小、分形维度D的增大，粗糙表面越光滑，在相同接触载荷下，接触面积越大，表面间的距离减小。

(a) 接触载荷和表面分离 (b) 接触载荷和接触面积

Figure 5. The relationship under the change of fractal characteristic length scale parameter G

图5. 分形特征尺度系数G的变化下的关系

4. 结论

1) 本文基于分形理论，建立了综合考虑微凸体间的相互作用及弹塑性变形时硬度随深度变化而变化的粗糙表面接触特性模型。

2) 单个微凸体发生弹性变形时，卸载过程与加载过程完全一致；发生弹塑性变形时，随着微凸体加载结束时的变形量的增大，在相同接触载荷下，卸载时的接触面积越大，且不能完全回复；发生完全塑性变形时，微凸体在卸载过程中不回复。

3) 随着尺度系数G的减小、分形维度D的增大，粗糙表面越光滑，在相同接触载荷下，接触面积越大，表面间的距离减小。