2. 西北工业大学 机电学院, 陕西 西安 710072
粗糙表面间接触特性的研究对分析其摩擦、磨损、导电、导热等性能具有重要影响。早期的研究主要是基于统计学分析的接触模型,采用的统计参数与采样长度和仪器分辨率相关,进而导致对确定粗糙表面的表征和分析结果不具有唯一性[1, 2]。分形几何理论提出后,迅速应用到粗糙表面的接触问题,利用分形理论建立的粗糙表面接触模型可以提供多尺度的接触行为预测分析。Majumdar等[3]提出以分形几何为基础的分形接触模型(MB模型),但该模型中未考虑微凸体的弹塑性变形,认为微凸体的临界弹性接触面积与尺度无关,得到接触过程中微凸体先发生塑性变形,后发生弹性变形,这一结果与传统的接触力学结果相反;Kogut等[4]用有限元法分析了粗糙表面上单个球状微凸体与刚性平面的接触情况(KE模型);Morag等[5]基于分形模型,应用Hertz理论证明了微凸体的临界接触面积与微凸体的尺寸相关,推导出了接触变形过程中微凸体先发生弹性变形,再发生非弹性变形;然而上述2种模型都只研究了粗糙表面上单个微凸体的变形机制,并没有考虑整个粗糙表面上的接触载荷与真实接触面积之间的关系。Liou等[6]研究了柱状粗糙表面与刚性光滑平面的接触,同样得到了随着变形量的不断增加,微凸体先发生弹性变形,再发生弹塑性变形或是完全塑性变形;缪小梅和丁雪兴等[7, 8]依据分形理论也推导出了单个微凸体的弹性临界接触面积与微凸体的尺度相关。
基于上述研究现状,依据分形理论,考虑微凸体的等级,确定了弹性临界等级、第一弹塑性临界等级和第二弹塑性临界等级的表达式,研究单个微凸体的力学特性,确定了不同等级微凸体的弹性临界接触面积、第一弹塑性临界接触面积和第二弹塑性临界接触面积。结合海洋岛屿面积分布规律建立了粗糙表面弹塑性接触力学模型,获得了不同等级微凸体下整个粗糙表面上接触总载荷与真实接触面积之间的关系。
1 粗糙表面几何模型的建立粗糙表面具有自仿射与多尺度的分形特性,Majumdar等[3]用Weierstrass-Mandelbrot函数(W-M函数)来表征粗糙表面二维轮廓曲线,表达式如下
式中,z(x)表示粗糙表面轮廓曲线的高度;x为轮廓的位移坐标;D为表面轮廓分形维数;G为轮廓特征尺度系数;γ为大于1的常数,对于服从正态分布的表面,取γ=1.5[3, 5];γn表示轮廓曲线的空间频率,nmin为与轮廓结构最低截止频率对应的序数。
轮廓曲线由D、G和nmin3个参数决定,由于表面轮廓为非平稳的随机过程,最低截止频率跟取样长度有关,D和G可由二维W-M函数的功率谱获得。
2 两粗糙表面接触两粗糙表面间的接触可以等效为一个粗糙表面与一个刚性光滑平面的接触,且等效粗糙表面满足各向同性的分形特性;忽略接触过程中微凸体之间的相互作用;假设变形时只有微凸体发生变形且不考虑接触过程中接触硬化和硬度随深度的变化;不考虑摩擦。
2.1 单个微凸体力学模型的建立粗糙表面是由一系列不同尺寸的余弦波状微凸体叠加而成的,对于任何一个微凸体l=1/γn,变形前的轮廓曲线可定义为
图 1为单个接触微凸体变形示意图,其中l为微凸体基底的尺寸,δ为微凸体的幅值,ω为实际变形量,取值范围为0≤ω≤δ,具体大小与所受载荷有关,l′为用一刚性平面截微凸体所得的截断长度,lr为微凸体变形量为ω时的实际接触长度。
 |
| 图 1 接触微凸体变形示意图 |
由(2)式可得微凸体顶端曲率半径R为
由(2)式可得微凸体的幅值δ为
当截断长度为l′时,微凸体的实际变形量ω为 2.2 微凸体的变形机制
两粗糙表面接触过程中,随着变形量的不断增加,表面上的微凸体可能存在3种变形状态,分别为弹性变形、弹塑性变形和完全塑性变形。下面分别讨论单个微凸体处于3种不同变形状态时的接触面积和接触载荷。
2.2.1 弹性变形状态存在条件根据Hertz理论[9],将开始塑性变形的临界变形量为
式中，H为较软材料的硬度,K为硬度系数,与材料的泊松比v相关,满足K=0.454+0.41v,E为等效弹性模量,ϕ=H/E为材料属性。
当ω≤ωec时,微凸体发生弹性变形,由Hertz理论可得,单个微凸体的接触面积a和接触载荷F分别为
当微凸体变形量为ωec时,可认为微凸体处于弹性变形范围内,由(7)式可知,微凸体的弹性临界接触面积为
对比(7)式和(9)式可得:当ω≤ωec时,有a≤aec,即单个微凸体的接触面积小于或等于弹性临界接触面积时发生弹性变形。
联立(7)式和(8)式得单个微凸体发生弹性变形时,接触载荷F与接触面积a之间的关系为
由(10)式可知,单个微凸体发生弹性变形时,接触载荷与接触面积的3/2次方成正比。
2.2.2 弹塑性变形状态存在条件Kogut等[4]研究了粗糙表面单峰接触的力学特性,结果表明:当ωec＜ω≤110ωec时,接触微凸体发生弹塑性变形,并将微凸体的弹塑性变形分为2个区域:
第一弹塑性变形区(ωec＜ω≤6ωec)
第二弹塑性变形区(6ωec＜ω≤110ωec)
易得第一弹塑性和第二弹塑性的临界变形量分别为
将(6)式、(8)式和(9)式代入(11)式和(12)式中,可得单个微凸体发生弹塑性变形时,接触载荷F与接触面积a之间的关系为
式中，aepc和apc分别为第一弹塑性临界接触面积和第二弹塑性临界接触面积,其大小分别为aepc=7.119 7aec、apc=205.382 7aec由(15)式和(16)式可知,单个微凸体发生弹塑性变形时,接触载荷与接触面积近似成线性关系。 2.2.3 完全塑性变形存在条件
当ω>ωpc时,微凸体发生完全塑性变形,接触面积a和接触载荷F的表达式如下
综上所述,单个微凸体的临界接触面积(弹性临界接触面积、第一弹塑性临界接触面积和第二弹塑性临界接触面积)是尺度相关的,符合经典赫兹接触理论[7],与传统的分形接触模型不同。对于确定的材料属性和分形参数,微凸体临界弹性接触面积与微凸体的等级成反比。任一微凸体受力发生变形,接触面积a满足:a≤aec时为弹性变形;aec＜a≤aepc时为第一弹塑性变形;aepc＜a≤apc时为第二弹塑性;a>apc时为完全塑性变形。
2.3 微凸体等级n的计算当δ≤ωec,微凸体发生弹性变形,即:
同理可以得到第一弹塑性临界等级nepc为
第二弹塑性临界等级npc为
综上所述,当nmin≤n≤nec时,微凸体只发生弹性变形;当nec＜n≤nepc时,微凸体发生弹性变形和第一弹塑性变形;当nepc＜n≤npc时,微凸体发生弹性变形、第一弹塑性变形和第二弹塑性变形,当n>npc时,微凸体发生弹性变形、第一弹塑性变形、第二弹塑性变形和完全塑性变形。
2.4 微凸体的面积分布密度函数Wang和Komvopoulos等[10]在修正的M-B分形接触模型中指出,分形粗糙表面与刚性光滑平面接触时,微凸体的接触面积分布密度函数为
真实接触面积Ar为 式中，al为最大微凸体的接触面积,φ为分形区域扩展系数,与轮廓分形维数D之间的函数表达式[8]为 2.5 真实接触面积和接触载荷
两粗糙表面接触时的真实接触面积Ar和总的接触载荷Fr是所有接触微凸体的接触面积和接触载荷的总和,计算如下
真实接触面积为
式中，Are为弹性变形部分的接触面积,Arep1为第一弹塑性变形部分的接触面积,Arep2为第二弹塑性变形部分的接触面积,Arp为完全塑性变形部分的接触面积。下面分别计算这四部分的真实接触面积。
总的接触载荷为
同理,Fre为弹性变形部分接触载荷,Frep1为第一弹塑性变形部分接触载荷,Frep2为第二弹塑性变形部分接触载荷,Frp为完全塑性变形部分接触载荷。
将(31)~(34)式代入(30)式并进行无量纲化处理即两边同时除以EAa,可得总的接触载荷与真实接触面积之间的无量纲关系
式中,
,为无量纲参数;lmax为最大微凸体尺寸,Aa为名义接触面积; 
,为与表面分形维数相关的常数,因此,对于特定的材料属性ϕ和分形参数D、G,无量纲接触载荷Fr*是微凸体的无量纲弹性临界接触面积aec*=aec/Aa和无量纲真实接触面积Ar*的函数,避免对粗糙表面进行弹塑性及完全塑性的重复计算。
由于单个微凸体临界接触面积的尺度相关性,在求解真实接触面积和接触载荷过程中,临界接触面积选取影响结果的准确性,在本文中我们选择最大微凸体的临界接触面积(对应等级n最小)来计算整个粗糙表面的接触载荷和真实接触面积,以下对其误差进行分析。计算参数选取D=1.5,G=2.5×10-9m,ϕ=0.001。
1) 当所有微凸体都处于同种变形状态时,同时处于弹性变形状态、第一弹塑性变形状态、第二弹塑性变形状态或完全塑性变形状态。粗糙表面的真实接触面积和接触载荷可以用最大微凸体的临界接触接触面积来计算[7]。
当an≤anec,所有微凸都体处于弹性变形状态;当anec＜an≤anepc所有微凸体都处于第一弹塑性变形状态;当anepc＜an≤anpc,所以微凸体都处于第二弹塑性变形状态;当an≥anpc,所有微凸体都处于完全塑性变形状态。
2)、当微凸体处于不同变形状态时,分为以下3种情况:
(1) 当最大微凸体处于弹性变形状态时,其余微凸体可能处于第一弹塑性变形状态。以本文计算结果,真实接触面积为
(2) 当最大微凸体处于弹性变形状态时,其余微凸体可能处于第二弹塑性变形阶段。同理可得真实接触面积的比值为:
(3) 当最大微凸体处于弹性变形状态时,其余微凸体可能处于完全塑形变形阶段。同理可得真实接触面积的比值为
为了说明我们选择最大微凸体来计算真实接触面积和载荷的合理性。我们又选择了小一级和小两级的微凸体,计算了其真实接触面积的与最大微凸体真实接触面积的比值。
经过计算,选择最大微凸体相邻的微凸体得到的总的真实接触面积与最大微凸体所得的真实接触面积的比值都小于50%,为了获得准确的计算结果,可以用最大微凸体的临界接触面积来计算整个粗糙表面的真实接触面积和载荷。
经过误差分析,在微凸体下压量与微凸体高度比值相同的情况下,最大微凸体的接触面积是最大的,因而粗糙表面中最大微凸体的力学性能直接决定整个粗糙表面的力学性能。当最大微凸体发生弹性变形,即使较大等级的微凸体发生弹塑性变形,对其整个粗糙表面的真实接触面积的影响很小,从工程应用和保守计算角度分析,选取最大微凸体的临界接触面积计算整个粗糙表面的真实接触面积和接触载荷是合理的。
3 数值仿真结果与分析基于前文所推导出的计算结果,分析粗糙表面接触过程中的真实接触面积与微凸体总接触载荷之间的关系。粗糙表面相关参数为:轮廓分形维数D分别取1.5、1.7、1.9,轮廓特征尺度系数取G=2.5×10-16～2.5×10-7,泊松比ν=0.3,材料的弹性模量E=2.3×1011N/m2,硬度为H=6.58×108N/m2。
图 2为微凸体的临界接触面积(弹性临界接触面积、第一弹塑性临界接触面积和第二弹塑性临界接触面积)与微凸体等级之间的关系。
 |
| 图 2 微凸体的临界接触面积与等级之间的关系 |
粗糙表面的分形参数为D=1.5,G=2.5×10-9。其中:Ⅰ为弹性变形区;Ⅱ为第一弹塑性变形区;Ⅲ为第二弹塑性变形区;Ⅳ为完全塑性变形区。在Ⅰ区域中,这些等级的微凸体只能发生弹性变形;在区域Ⅱ中,这些等级的微凸体可以发生弹性变形和第一弹塑性变形;在区域Ⅲ中,这些等级的微凸体可以发生弹性变形、第一弹塑性变形和第二弹塑性变形;在区域Ⅳ中,这些等级的微凸体可以发生弹性变形、弹塑性变形和完全塑性变形。从图 2中可以看出,当表面的分形参数一定时,单个微凸体的临界接触面积aec,aepc,apc是变化的,是和微凸体的尺度相关的,随着微凸体的等级n的增大而减小。对于确定等级的微凸体,随着接触面积的增大,微凸体先发生弹性变形,再发生弹塑性变形,最后发生完全塑性变形。
图 3为在轮廓的分形维数D=1.5的情况下,微凸体的临界等级(弹性临界等级、第一弹塑性临界等级和第二弹塑性临界等级)与轮廓的特征尺度参数之间的关系。由图 3可知,3种临界等级都随着轮廓尺度参数的增大而减小。在相同的特征尺度参数下,第二弹塑性临界等级最大,弹性临界等级最小,而第一弹塑性临界等级介于这2种临界等级之间。当轮廓尺度参数G增大到一定值时,临界等级小于零,但事实上,最大微凸体的等级大于零的,因此当分形维数为D=1.5时,轮廓的特征尺度参数不能大于2.5×10-7m。
 |
| 图 3 微凸体的临界等级与轮廓特征尺度参数之间的关系 |
图 4a)~图 4c)分别表示最大微凸体的等级处于不同的临界等级范围内时无量纲接触载荷Fr*与无量纲真实接触面积Ar*之间的关系。如图 4a)所示:当最大微凸体的等级小于弹性临界等级时,即n=10＜nec=13,微凸体只能发生弹性变形,因此在整个粗糙表面接触过程中,Fr*与Ar*的3/2次方成正比关系,粗糙表面表现为弹性性质;如图 4b)所示,当最大微凸体的等级大于弹性临界等级而小于第一弹塑性临界等级时,即n=16＜nec=18,随着下压量的不断增大,微凸体先发生弹性变形,再发生弹塑性变形,因此在整个粗糙表面接触过程中,在载荷较小的情况下,Fr*与Ar*的3/2次方成正比关系,随着载荷的增加,Fr*与Ar*的1.254 4次方成正比,因而粗糙表面先表现为弹性性质,随着载荷的增大表现为弹塑性性质。如图 4c)所示,当最大微凸体的等级大于第二弹塑性临界等级时,即n=34>npc=25,微凸体的变形顺序为弹性变形,弹塑性变形和完全塑性变形,因此在整个粗糙表面接触过程中,粗糙表面的接触性质也是这一顺序。图 4c)中Fr*与Ar*之间成线性正比的关系,这是因为接触过程中弹性变形和弹塑性变形所占的比例很小。由图 4可知,随着轮廓分形维数D的增大,Fr*与Ar*的直线斜率逐渐增大,在相同接触面积上接触载荷减小,这是因为:表面分形维数D越大,粗糙表面越精细,单位面积上的微凸体数量增多。
 |
| 图 4 无量纲接触载荷与无量纲真实接触面积之间的关系 |
图 5为本文模型、GW模型[1]、MB分形模型[3]与Bhushan所测的实验数据[11]的对比,所有分形模型的计算参数均取[3]:D=1.49,G*=10-10,ϕ=0.05。本文的模型与其他模型具有一致的变化趋势,即无量纲载荷增加的情况下,无量纲真实接触面积也增大。从图 5a)可得,与MB模型选取相同基底长度的微凸体,从而可以确定该微凸体所对应的等级n=34。在较小载荷的情况下,本文模型的计算结果与Bhushan实验数据比较接近,GW模型的计算结果偏差较大;从图 5b)可得,当最大微凸体的等
 |
| 图 5 不同接触模型与实验对比 |
级为n=8时,在载荷较小时,本文模型的计算结果误差比较大:当无量纲接触载荷满足2×10-4＜Fr*＜5×10-4时,本文模型与Bhushan实验数据最为吻合,MB和GW模型的计算结果与实验数据有一定的偏差。从图 5可以看出,不管最大微凸体的等级在什么范围内,在高的载荷的情形下,G-W模型的计算结果与实验数据较为接近,本文模型和MB分形模型的结果与实验结果偏差较大。比较结果说明,接触载荷处于中低载荷范围时,分形模型的计算结果优于统计结果,接触载荷处于高载范围,统计模型计算结果较为准确。
5 结 论1) 分形粗糙表面中单个微凸体的临界接触面积aec,aepc,apc是和微凸体的尺度相关的,临界接触面积是由材料属性和微凸体等级n共同决定。在确定的材料属性条件下,单个微凸体临界弹性接触面积与微凸体等级成反比。
2) 材料属性和分形参数一定时,粗糙表面上的接触载荷与最大微凸体接触面积和其对应的临界弹性接触面积相关,实际计算时,只要得到分形粗糙表面上最大微凸体接触面积和其对应的弹性临界接触面积就可获得整个粗糙表面上的接触载荷与真实接触面积,避免对粗糙表面进行弹塑性及完全塑性的重复计算,简化计算步骤。
3) 粗糙表面接触过程中,最大微凸体所处的等级范围不同,粗糙表面表现出来的力学性能也不相同。在整个粗糙表面接触过程中,粗糙表面变形过程与单个微凸体的变形过程一致。接触面积一定时,轮廓分形维数增大,接触载荷减小。
| [1] | Greenwood J A, Williamson J B P. Contact of Nominally Flat Surfaces[J]. Mathematical and Physical Sciences, 1966, 295(1442): 300-319 |
| Click to display the text | |
| [2] | Chang W R, Etsion I, Bogy D B. An Elastic-Plastic Model for the Contact of Rough Surfaces[J]. ASME Journal of Tribology, 1987, 109: 257-263 |
| Click to display the text | |
| [3] | Majumdar A, Bhushan B. Fractal Model of Elastic-Plastic Contact between Rough Surfaces[J]. ASME Journal of Tribology, 1991, 113: 1-11 |
| Click to display the text | |
| [4] | Kogut L, Etsion I. Elastic-Plastic Contact Analysis of a Sphere and a Rigid Flat[J]. ASME Journal of Applied Mechanics, 2002, 69(5): 657-662 |
| Click to display the text | |
| [5] | Morag Y, Etsion I. Resolving the Contradiction of Asperities Plastic to Elastic Mode Transition in Current Contact Models of Fractal Rough Surfaces[J]. Wear, 2007, 262(5/6): 624-629 |
| Click to display the text | |
| [6] | Jeng Luen Liou, Chi Ming Tsai, Lin Jenfin. A Microcontact Model Developed for Sphere-and Cylinder-Based Fractal Bodies in Contact with a Rigid Flat Surface[J]. Wear, 2010, 268: 431-442 |
| Click to display the text | |
| [7] | Miao Xiaomei, Huang Xiaodiao. A Complete Contact Model of a Fractal Rough Surface[J]. Wear, 2014, 309: 146-151 |
| Click to display the text | |
| [8] |
丁雪兴, 严如奇, 贾永磊. 基于基底长度的粗糙表面分形接触模型的构建与分析[J]. 摩擦学学报, 2014, 34(4): 341-347 Ding Xuexing, Yan Ruqi, Jia Yonglei. Construction and Analysis of Fractal Contact Mechanics Model for Rough Surface Based on Base Length[J]. Tribology, 2014, 34(4): 341-347 (in Chinese) |
| Cited By in Cnki (1) | Click to display the text | |
| [9] | Johnson K L. Contact Mechanics[M]. London: Cambridge University Press, 1985: 79-128 |
| [10] | Wang S, Komvopoulos K. A Fractal Theory of the Interfacial Temperature Distribution in the Slow Sliding Regime: PartⅡ——Multiple Domains, Elastoplastic Contacts and Applications[J]. ASME Journal of Tribology, 1994, 116(4): 824-832 |
| Click to display the text | |
| [11] | Bhushan B. The Real Area of Contact in Polymeric Magnetic Media-PartⅡ: Experimental Data Analysis[J]. ASLE Transactions, 1985, 28: 181-197 |
| Click to display the text |
2. School of Mechanical Engineering, Northwestern Polytechnical University, Xi'an 710072, China
