基于分形理论的两粗糙表面接触的黏滑摩擦模型

周 华，龙新华，孟 光

（上海交通大学机械系统与振动国家重点实验室，上海 200240）

引 言

机械结构中存在大量的连接界面，例如螺栓连接。这些连接界面受到法向紧固载荷后，能够承受切向载荷。从微观角度看，在切向载荷作用下，连接表面大量的微凸体接触会发生黏着和滑移等行为，使连接界面呈现出复杂的非线性行为。建立两粗糙表面接触的切向黏滑摩擦模型一直是具有挑战性的问题。

两粗糙表面的接触模型可以分为统计模型和分形模型。Greenwood 等［1］假设微凸体高度满足高斯分布，并基于统计方法建立了法向接触模型（GW 模型）。不少学者在此模型上进行了改进，但是这种方法受到测量仪器的分辨率和采样长度的影响。因此，Majumda 等［2］用两个分形参数表征了表面的形貌，并提出了分形接触模型（MB 模型），克服了统计学方法的不足。Wang 等［3］、Yan 等［4］和刘文威［5］对MB 模型进行了改进，使分形接触模型更加合理，然而这些模型并未考虑微凸体的相互作用。Zhao 等［6］考虑了微凸体相互作用造成的基底面下降，并应用在GW 统计模型中。此外，张学良等［7］、温淑花等［8］和田红亮等［9］基于分形理论建立了结合面的法向接触刚度模型。Zhou 等［10］考虑微凸体相互作用造成的基底面下降，建立了两粗糙表面的法向接触模型。

对于两粗糙表面接触的黏滑摩擦研究，Mind⁃lin［11］推导了在给定恒定法向载荷时，球接触的切向载荷与切向位移之间的非线性关系。文中认为圆形的接触区域由两部分组成，接触的中心区域为Stick黏着区，接触的外层为Slip 滑移区。当切向载荷增大时，滑移区会不断增大，直至发生宏观滑动。Mindlin 等［12］给出了处于弹性接触的球体在循环加载时的响应。Ödfalk 等［13］把弹性接触扩展成弹性⁃塑性接触情况。Björklund［14］基于单个微凸体的Mindlin 模型和GW 统计模型分析了接触粗糙表面的黏滑现象，但仅考虑了微凸体的弹性变形。Er⁃iten 等［15⁃17］考虑了微凸体的弹塑性变形，并且根据不同的摩擦系数模型建立了黏滑摩擦模型。Wang等［18⁃19］考虑了微凸体的弹性和完全塑性变形，但也是根据统计模型分析了粗糙表面的切向接触特性。Jamshidi 等［20］考虑了微凸体的斜接触，即切向和法向变形耦合，提出了改进的两粗糙表面接触模型。该模型直接分析两个粗糙表面接触，无需等效成一个光滑表面和一个等效粗糙表面，但也是基于GW统计模型。Song 等［21］尝试基于分形模型建立界面黏滑摩擦模型，但并未考虑微凸体的相互作用。

本文采用文献［10］中的法向接触模型，并确定在给定法向预紧载荷下微接触截面积的概率密度函数。同时，本文基于Mindlin 模型和Masing 准则获得了单个微凸体的切向模型，然后根据分形理论对微接触截面积进行积分，建立了整个粗糙表面接触的黏滑摩擦模型，并研究了分形参数D和G以及法向预紧力对接触表面能量耗散的影响。

1 两粗糙表面法向接触模型

1.1 考虑微凸体相互作用的分形接触模型

两个粗糙表面接触可以等效成一个光滑表面和一个等效的粗糙表面［22］，等效表面的等效弹性模量为：

式中E1和E2分别为两表面材料的弹性模量，ν1和ν2分别是两表面材料的泊松比。

图1为一个光滑刚性表面和一个等效粗糙表面接触的示意图。δn表示微凸体的总变形量。作用在其他微凸体上的压力会对给定微凸体产生额外的变形，用微凸体平均高度的下降距离u表示微凸体相互作用产生的影响，此时单个微凸体真实的干涉量ω=δn-u。

图1 一个刚性光滑表面与粗糙表面的接触Fig.1 Contact of a rough surface with a rigid smooth flat

微凸体的总变形量δn可以由粗糙表面的形貌获得。Yan 等［4］使用改进的两变量W⁃M 函数来表征粗糙表面的形貌。他们假设接触表面是各向同性表面，并且推导了等效的单变量等式。

式中Ds（2

式中 int［…］表示取中括号数值的最大整数。若微凸体的截半径为r′，则微凸体波形的最长波长为2r′。可以合理地假设微接触力主要是由于基波长的微凸体尖端的小变形造成的，并且相应的频率指数为n0=ln(L/2r′)/lnγ［4］。代入到式（2），可以得到基波轮廓曲线函数z0（x），表示为：

δn定义为关于x轴的余弦函数的峰值高度，写为：

式中c1=22-Dπ(D-2)/2GD-1(lnγ)1/2，a′为微接触截面积，即a′=πr′2。

由于微凸体相互作用造成的微凸体平均高度下降的位移u可以表示为［6］：

式中Fn为总法向载荷，Aa为名义接触面积，Fas为一个微凸体上的力。

考虑微凸体相互作用后，单个微凸体的真实干涉量可以表示为：

需要指出的是，当微凸体处于不同变形阶段时，上式中的Fas是不同的。

研究表明，微接触截面积的分布满足地球岛屿面积分布，微接触截面积的尺寸分布函数n(a′)，或称为概率密度函数，可以表示为：

式中a′l为最大的微凸体截面积，ψ为域扩展因子，与分形维数相关，表示为：

对等式（9）求解并拟合数据，可以直观表示为：

微凸体的曲率半径R可以从等式（4）获得：

式中c2=[22-Dπ2+D/2GD-1(lnγ)1/2]-1。

微凸体有三个变形阶段：弹性变形、弹塑性变形和完全塑性变形。ωyc定义为微凸体从弹性变形阶段转为弹塑性变形阶段的临界屈服变形值，ωpc定义为微凸体从弹塑性变形阶段转为完全塑性变形阶段的临界变形值，可以表示为［22⁃23］：

式中H表示较软材料的硬度，H≈2.8σb；K为硬度系数［24］，且K= 0.4645 + 0.3141ν+ 0.1943ν2表示初始屈服点的平均压力，=1.1σy；表示完全塑性流动转变点处的均压［23］，=2.8σy。

把等式（11）代入等式（12），临界屈服变形量为：

（i）弹性变形阶段

微凸体在此阶段的微接触面积ae和接触载荷Fe可以表示为［2］：

此时，真实干涉量ω可以表示为：

根据等式（17）中a′和ω的关系，临界微接触屈服截面积a′yc可以由等式（12）中的ωyc求得；类似的，临界微接触截面积a′pc可以由等式（13）中的ωpc求得。

因为c4与法向载荷Fn相关，所以这两个临界接触面积并不是常值，而是与法向载荷相关。

（ii）弹塑性阶段

微凸体在此阶段的微接触面积aep和接触载荷Fep可以表示为［5］：

式中b和s为常数，可以根据微接触面积和载荷在弹性、弹塑性和完全塑性变形的连续性条件计算出这两个参数。

此时，真实干涉量ω可以表示为：

（iii）完全塑性阶段

微凸体在此阶段的微接触面积ap和接触载荷Fp可以表示为［22］：

此时，真实干涉量ω可以表示为：

1.2 法向总载荷

单个微凸体在弹性、弹塑性和完全塑性阶段的法向载荷分别用Fe，Fep和Fp表示。根据尺寸分布函数n(a′)，这三个阶段总的接触载荷分别用Fne，Fnep和Fnp表示，由下面等式给出：

需要指出的是，这些积分的被积函数是隐式表达式，三个阶段的ω和a′的关系分别由等式（17），（23）和（26）给出，并基于Maltab 软件的fzero 函数求解这些积分。

连接界面总的法向载荷是三个变形阶段法向载荷的总和，可以表示为：

2 两粗糙表面的黏滑摩擦建模

在法向载荷的作用下，两接触表面被压紧，受到较小的往复切向载荷或切向位移时，会产生往复的切向位移，出现滞回现象。从单个微凸体角度看，有的微凸体接触处于黏着状态，而有的微凸体接触处于滑移状态。当切向载荷不断增大时，短的微凸体受到的法向载荷比高的微凸体小，所以更容易先滑动。同理，微凸体接触截面积越小，越容易滑动。在首次加载时，接触截面积大的黏着微凸体构成了切向刚度，而接触截面积小的滑移微凸体导致了能量耗散。卸载时，切向载荷减小，而预紧力不变，即极限切向载荷不变，所以一些在首次加载时滑移的微凸体会转变为黏着状态。

图2和3 展示了微凸体在不同接触状态时微接触截面积概率密度函数的区间。首次加载时，接触的微凸体状态可以分成两组。而卸载时，接触的微凸体状态分成三组［17］。

（i）微凸体在卸载和加载时都滑移。

（ii）微凸体在加载时滑移，但卸载时处于黏着状态。

（iii）微凸体在卸载和加载时都处于黏着状态。

图2和3 中，a′s是最小的微凸体接触截面积，可以近似为0。a′0，a′1和a′2都是临界接触截面积，用于区别首次加载与卸载时微凸体的接触状态，下文给出表达式。

图2 首次加载时，微凸体在不同接触状态时微接触截面积概率密度函数的区间Fig.2 Categories of asperities when initial loading based on the truncated microcontact area distribution

图3 卸载时，微凸体在不同接触状态时微接触截面积概率密度函数的区间。Fig.3 Categories of asperities when unloading based on the truncated microcontact area distribution

Mindlin［12］给出了切向载荷与切向相对位移之间的关系：

式中δ为相对切向位移，T为切向载荷，r为接触半径，Fas为单个微凸体承受的法向力，μ为摩擦系数，G′为等效剪切模量。

根据等式（31），切向载荷与位移之间的关系为：

无量纲切向相对位移用弹性阶段的接触面积与法向载荷表示为：

当微凸体处于弹塑性阶段时，等式（34）中的比值r/Fas会不同，但是至多导致10% 的误差［17］。为了便于计算，本文都采用了弹性接触时的表达式。 摩擦系数μ采用CFC 恒定值模型，并取μ= 0.2。

当切向载荷不足以引起整体滑动时，处于滑移状态的微凸体，其切向载荷满足Coulomb 摩擦定律；处于黏着状态的微凸体，其切向载荷与相对位移满足Mindlin 理型。此时，可以得到：

（i）首次加载

式中Tas为单个微凸体的切向载荷，a′0为临界微接触截面积，对应的临界变形量ω0=4G′δ/(μE)，由于ω与a′存在隐式关系，无法直接给出显式表达式，但可以由等式（17），（23），（26）确定。

根据微凸体截面积的尺寸分布函数n(a′)，面⁃面接触的总切向载荷为：

需要指出的是，等式（36）需要进行分段积分。通过判断a′0，a′pc和a′yc的大小关系，确定积分区间。例如当a′pc

（ii）卸载

式中a′1和a′2为临界接触截面积，对应的临界变形量为ω1=2G′(δmax-δ) /(μE) 和ω2=4G′δmax/(μE)，ω与a′的关系由等式（17），（23），（26）确定。

同理，计算卸载时面⁃面接触的总切向载荷时，需要判断a′1，a′2，a′pc和a′yc的大小关系。例如当a′1≤a′pc

（iii）再次加载

微凸体的卸载方程和再次加载方程满足Mas⁃ing 映射准则，可以得到再次加载时面⁃面接触的总切向载荷：

（iv）能量耗散

在周期载荷作用下，面⁃面接触的力⁃位移的曲线形成了封闭区域，封闭区域的面积即为一个周期内的能量耗散，可以表示为：

3 数值仿真

为了进一步研究所提出的基于分形理论的切向黏滑摩擦模型的特性，开展了数值仿真研究。切向载荷与相对位移关系曲线的计算流程如图4所示。无量纲化公式为：

图4 切向总载荷的计算流程图Fig.4 Calculation of total tangential load

仿真参数为：E1=E2= 210 GPa，ν1=ν2=0.3，H= 2.744 GPa，σy= 785 MPa，Aa= π ×（10/1000）2m2，μ= 0.2。其他参数分成3 组，用于对比研究不同参数对切向黏滑摩擦模型的影响。

（1）G= 1×10-7m，Fn= 1000 N，D分别取1.60，1.62，1.64。

（2）D= 1.60，Fn= 1000 N，G分别取0.8×10-7，1×10-7，1.2×10-7m。

（3）D= 1.60，G= 1×10-7m，Fn分别取1000，1050，1100 N。

图5展示了不同的分形维数D时，连接界面的迟滞曲线。从图中看出，在相同切向载荷下，D增大时，迟滞曲线所围面积增大，即一个周期的能量耗散增大。D= 1.60，1.62 和1.64 的能量耗散分别为1.202×10-4，1.408×10-4和1.591×10-4J。分形维数D反映表面轮廓的复杂性和自相似性，一般来说，D越大，表面越光滑，实际接触面积增加，滑移微凸体的个数增加，所以能量耗散也增大，该趋势与文献［25⁃26］中的结果一致。

图5 不同D 时，G = 1×10-7 m，Fn = 1000 N 时的切向迟滞曲线Fig.5 Hysteresis loop for different fractal dimension D with G = 1×10-7 m，Fn = 1000 N

不同分形粗糙度参数G时，连接界面的迟滞曲线如图6所示。从图中可以看出，在相同切向载荷下，G增大时，迟滞曲线所围面积减小，即能量耗散减小。G= 0.8×10-7，1×10-7和1.2×10-7的能量耗散分别为1.370×10-4，1.202×10-4和1.061×10-4J。分形粗糙度参数G反映粗糙表面的垂直方向的光滑程度，G越大，则表面越粗糙。所以图6说明表面越粗糙，能量耗散越小，与前文结果一致。

图6 不同G 时，D = 1.60，Fn= 1000 N 时的切向迟滞曲线Fig.6 Hysteresis loop for different fractal roughness G with D = 1.60，Fn = 1000 N

图7展示了不同法向紧固载荷时，连接界面的迟滞曲线。从图中可以看出，在相同切向载荷下，Fn增大时，迟滞曲线所围面积减小，即能量耗散减小。Fn= 1000，1050，1100 N 时的能量耗散分别为1.202×10-4，8.447×10-5和6.287×10-5J。这是因为法向载荷的增大会导致微凸体滑移的最大切向载荷增大，在切向外载荷不变时，滑移微凸体的个数会减少，所以能量耗散减小，该趋势与文献［17，26］的结果一致。

图7 不同Fn时，D=1.60，G=1×10-7 m 时的切向迟滞曲线Fig.7 Hysteresis loop for different normal load Fn with D =1.60，G = 1×10-7 m

4 结 论

本文考虑了微凸体相互作用对基底面的下降，并基于分形理论建立了两粗糙表面的法向接触模型，从而确定在给定法向预紧载荷下微接触截面积的概率密度函数。然后根据Mindlin 模型、Masing准则和分形理论，建立了两粗糙表面接触的切向黏滑摩擦模型。最后开展了数值仿真研究，分析了不同参数对系统能量耗散的影响。结论如下：

（1）分形维数D越大，能量耗散越大；分形粗糙度参数G越大，能量耗散越小，即表面越粗糙，能量耗散越小。

（2）法向预紧载荷增大时，微凸体滑移的最大切向载荷会增大，在切向外载荷不变时，滑移微凸体的个数会减少，所以能量耗散越小。