含硬涂层结构的表面微观接触新模型

李 玲, 李治强， 张 健, 蔡安江

(西安建筑科技大学 机电工程学院, 西安 710055)

制造业的不断发展，对机械加工精度和加工效率提出了更高的要求，机床的动态性能在实际生产中就显得非常重要。机床结合面的动态参数对机床整体动态特性起着决定性影响，它们决定机床的加工精度和稳定性。因此，开展结合面动态特性研究对揭示机床整体动态特性具有重要的理论意义和实际应用价值。

结合面的研究主要是粗糙表面间接触行为的研究，而关于粗糙表面间接触行为的研究以微观表面微凸体的研究为主[1]。国内外大量学者对粗糙表面微凸体采用不同的研究方法，提出了不同的粗糙表面接触模型。其中应用最为广泛的表面接触模型是Greenwood等[2]提出的统计接触模型(Greenwood-Williamson, GW模型)。GW模型是一个根据Hertz理论提出的粗糙表面统计接触模型，模型假设：① 微凸体的峰顶是球体；② 具有相同的曲率半径；③ 微凸体高度服从高斯分布；④ 表面粗糙度各向同性；⑤ 微凸体只进行弹性接触，且各个微凸体的变形相互独立。大多学者基于GW模型，对微观接触模型进行了拓展，Kogut等[3]提出描述粗糙表面的弹性、弹塑性和全塑性形变的接触模型；Chang等[4]根据微凸体塑性变形体积守恒原理，提出了塑性接触模型；赵永武等[5-7]基于接触力学理论提出一种新的粗糙表面弹塑性微观接触模型(简称ZMC模型)；温淑花等[8-9]基于球体与平面的接触理论和粗糙表面的接触分形理论，提出接触刚度和接触阻尼模型；尤晋闽等[10-11]利用分形接触理论，建立了结合面法向接触分形模型；李玲等[12]根据在变形状态转变的临界点处微凸体真实接触面积与接触载荷转化均满足连续条件，建立了结合面真实接触面积、接触载荷与接触刚度模型；肖会芳等[13]采用GW模型描述粗糙表面形貌，建立了不同润滑状态下的滑动粗糙界面模型，研究了界面载荷分配、油膜厚度、摩擦因数和临界速度随界面形貌参数、法向载荷、润滑油属性参数的变化规律；杨红平等[14]基于分形几何理论和接触力学理论，用分形理论表征粗糙表面微凸体参数，建立各变形阶段微凸体的接触刚度模型。

上述学者都通过单微凸体接触来描述粗糙表面的变形，都假设粗糙表面的微观形态是许多微凸体附着在一个基体上，且只考虑了接触过程中微凸体的变形情况，并没有考虑微凸体所附着的基体的变形。当表面微凸体的材料和基体材料相同时，GW接触模型和它的扩展改进模型通常能有效地分析接触行为。然而，对于一些有硬涂层的粗糙表面，其接触特性主要由软基体上的硬涂层的材料决定，较软基体的变形对粗糙表面变形也有相当大的影响[15]。传统的GW接触模型它的扩展改进模型不能完全描述它的接触行为。Shi等[16]测量了一种磁记录材料的接触刚度，测得刚度值大约是106N/m，而通过拓展GW模型计算得出的接触刚度值为107～108N/m。测量得到的接触刚度测量值比拓展GW模型得到的计算值低，这也证明了基体变形对粗糙表面接触特性的影响。

为了准确地预测有硬涂层粗糙表面的接触行为，接触模型中应包含较软的基体变形和较硬的微凸体变形。在各位学者研究的基础上，通过分析较软基体上有较硬涂层的粗糙表面接触行为，提出了一个在GW接触模型基础上改进的，并同时考虑基体变形和微凸体变形的新弹性接触模型。通过对新模型进行有限元模拟，并将新模型与有限元分析结果、Hertz接触模型进行对比，发现新模型可以有效地描述有硬涂层粗糙表面的真实接触情况。

1 接触过程中的接触变形和应力分布

如图1所示，当两个球体发生轴对称接触时，R1和R2为上下两个半球体的曲率半径；E1和E2为上下两个半球体材料的弹性模量；υ1和υ2为上下两个半球体材料的泊松比。假设没有切向滑移，最大的剪切应力将出现在较软材料球体的接触表面上。在弹性变形范围内，随着载荷的增加，最大剪应力的位置会离接触表面越来越远。当接触球体有较大的曲率半径和较高的屈服强度或硬度时，最大剪应力值将更大，其位置也将距接触表面更远。如果材料表面涂有硬涂层，接触时应力分布将会不均匀，简单的Hertz接触模型无法准确地描述这种接触。

图1 球体-球体接触图

为证明这个观点，将硬涂层结构简化为一个球面微凸体和一个平整基体结合，并与一个刚性平面进行接触，具体形式如图2所示，采用ABAQUS有限元软件对新接触模型进行分析模拟。图2中微凸体的半径都取R=15 μm，微凸体的高度都取ha=100 nm，不同的是，一个系统中微凸体和基体的材料是相同的，它们的弹性模量Ea=Eb=100 GPa；另一个系统中微凸体材料较硬于基体材料，即Ea=200 GPa，Eb=100 GPa。为了简化对比，两个系统中微凸体和基体材料的泊松比取υa=υb=0.3。

图2 单微凸体/基体系统

假设微凸体和基体均是弹性体，将微凸体和基体模型都采用轴对称八节点等参数单元方法进行有限元分析。为了加快分析速度，网格划分时，接近接触表面的区域使用较密的面积为1 nm×1 nm的网格，远离测试区域的网格可以相对粗化一些。单微凸体/基体系统和它接触的边界条件如图2所示。微凸体绑定在基体之上，微凸体和基体之间无切向滑移，微凸体和基体的x轴方向和基体的y轴方向是固定的。微凸体上方的刚性平面上下移动使微凸体和基体发生形变。为避免由边界条件引起的误差，将基体的体积设置得比微凸体体积大一些，这样就形成了一个无限弹性半空间体。

图3是有限元模拟时，当实际位移量δ=4.5 nm时，微凸体系统所产生的接触变形和应力分布图。其中，图3(a)为微凸体和基体是同种材料的情况，从图可以看出它的应力场是连续且均匀的散开在微凸体和基体之中。这种情况与传统的微凸体包含在基体中的Hertz接触模型一致，从而证明传统GW模型的球面接触使用Hertz接触解法的正确性。图3(b)为微凸体的材料较硬于基体的材料，从图中可以看出，微凸体系统的应力场分布不均匀，且最大应力出现在微凸体和基体的交界处。由于微凸体和基体材料特性的差异，Hertz接触模型无法解决这种复杂情况。因此，当接触表面微凸体材料比基体材料有较高的强度和硬度时，应该考虑较软基体对接触分析的影响。

(a) Ea=Eb=100 GPa

(b) Ea=200 GPa， Eb=100 GPa

从图3(b)可看出，有限元模拟时，由于假设微凸体和基体之间没有发生相对滑移，且微凸体和基体的材料不同，所以在加载条件下，微凸体和基体之间的应力分布不连续。但是，本文研究对象由于所施加的负载或位移相对较小，接触变形均是在弹性范围内，所以将忽略微凸体和基体之间的相对滑移。

2 单微凸体接触分析建模

第1节有限元分析结果表明，微凸体和基体的变形对粗糙表面微观接触特性的影响不可忽略。考虑到基体材料和基体变形对整体接触的影响，可以分开计算微凸体和基体的变形量，分别建立微凸体和基体的刚度模型，然后将这两个刚度模型耦合形成系统总的刚度模型，基体变形量与微凸体/基体系统总变形量之间关系，微凸体刚度与基体刚度之间关系如图4所示。这种微凸体基体耦合模型得到的接触结果并不是理论上的精确解，因为它不包括材料的切向接触，接触压力主要取决于微凸体和基体的几何和物理特性。当微凸体材料比基体材料硬时，这种把微凸体和基体刚度分开计算的方法是体现微凸体和基体变形的一种简单方法。

图4 接触模型原理图

根据文献[17-18]，假设粗糙表面的微凸体为球体，其半径和高度可从粗糙表面的形貌测量得出。微凸体的刚度可以通过Hertz接触理论获得其表达式[19]

(1)

需要注意的是，微凸体的法向位移量δa和系统总的位移量δ并不是同一个量，系统总的位移量是微凸体位移量δa和基体位移量δb之和。当上部的微凸体受到载荷时，微凸体和基体的接触区是一个半径为rb的圆形区域，rb的值可通过微凸体的曲率半径R和微凸体的高度ha求得，可表示为

(2)

基体圆形接触区域内各点受到的压力各不相同，接触区域中心位置受到的压力最大，最大压力P0为

(3)

距离接触区域中心的距离r不同，基体表面发生的法向位移不同，当基体受到的最大压力P0在不断变化时，基体表面的法向位移为

(4)

根据平衡条件，圆形区域实际接触力F值和作用在微凸体上的法向接触力相同。所以，基体与微凸体所接触圆形触区域中心点的位移为

(5)

通过式(3)和式(5)，可以得到基体的刚度计算公式为

(6)

大部分的微凸体的球形区域被认为是包含在基体中，所以式(1)中包括了基体接触变形的影响。基体的影响分别出现在式(1)和式(6)中，因此当微凸体和基体应用相同材料的时候，提出的接触模型计算出接触变形量将高于实际接触变形量。为了得到一个确切的接触刚度解，需对式(1)进行修改，排除其中基体对微凸体刚度的影响。然而，当微凸体材料比基体的材料硬时，较软基体的接触变形在整个接触过程中占主导地位，因此式(1)使用Hertz解的影响将会很小，最终导致的误差可以忽略。为了简化和方便应用，该微凸体的刚度模型使用Hertz解。

当微凸体和基体的刚度值是常数时，等效刚度k可以通过简单的公式1/k=1/ka+1/kb确定。然而，如式(1)所示，由于ka不是一个定值，而是一个关于微凸体形变量δa的一个方程，所以在求等效刚度之前，应该事先确定获得微凸体形变量δa的值。

根据力的相互作用原理，可得

k×δ=ka(δa)×δa=kb×δb

(7)

δ=δa+δb

(8)

将式(1)和式(6)代入式(7)式(8)，可以得到δa关于δ的公式为

(9)

通过化简和转化，可以将求解δa的函数为

δa=f(δ)

(10)

(11)

将式(11)中的δa替换为与δ、κ、ξ有关的函数g(δ)，式(11)为

(12)

利用不动点迭代法，可以得到式(12)的近似解，式(12)又为

(13)

理论上来说，起始点g0可以是任何值，但是为了避免符号运算造成的复杂性，假设初始点时微凸体和基体拥有相同的几何形状和不同材料特性，即ξ=1，从而得出g0为

(14)

因此迭代序列采用以下方式不断迭代

(15)

式中：n=0，1，…，∞，迭代的次数越多得到的结果越精确。本文只考虑弹性变形的前提下，经过计算发现，只需迭代两次得到的近似解就非常接近于计算得到的数值解，且误差不超过5%。因此微凸体的形变量可以近似地表示为一个关于实际位移量δ的函数，其为

(16)

图5描述了采用数值解法的式(10)和近似解法的式(16)求得的微凸体形变量在各种几何和物理特性之下的对比。微凸体材料比基体材料软时κ的值小于1，微凸体材料比基体材料硬时κ的值大于1，如图5(a)所示，实际位移量不变的情况下，随着微凸体材料不断变硬，微凸体的形变量变得越来越小，基体的形变量变得越来越大，这也证明当微凸体材料硬于基体材料时，微凸体/基体系统中基体的变形量占主导地位。由于一个平滑的表面拥有较大的微凸体半径R和较小的微凸体高度ha，通过式(11)可以得出，表面越光滑，表面粗糙参数ξ的值就越大。如图5(b)所示，实际位移量不变的情况下，ξ值越大，微凸体的形变量变得越来越小，基体的形变量变得越来越大。从这些结果中，可以清楚地看到，表面涂层较硬时，基体的变形对接触行为的影响很大，这和第1节中有限元模拟的分析结果一致。从图5还可以清楚地看出，由式(16)得出的微凸体形变量近似解(虚线)与式(10)得出的微凸体形变量数值解(实线)是基本相同的，所以我们在建模时使用其简单的近似解。

(a) 微凸体物理特性对形变量的影响

(b) 微凸体几何特性对形变量的影响

Fig.5 Comparisons of single asperity deformation with applied displacement between the numerical solution and approximate solution

综合式(1)、式(6)和式(16)，单微凸体系统的等效接触刚度可以用一个与实际位移量有关的方程表示

(17)

接触力可以通过式(18)确定

F(δ)=k(δ)×δ

(18)

3 模型验证与分析

为了验证第2节所提出接触模型的正确性，利用有限元软件模拟一个刚性平面同微凸体和基体的接触。对于单微凸体模型，微凸体半径为R，微凸体高度设为刚性面实际位移的三倍，即ha=3δ。其中基体上圆形接触区域半径rb可以通过式(3)得到。

表1 四种不同几何和物理特性的微凸体

如表1所示，4组不同的微凸体高度ha和微凸体平均曲率半径R代表现实中的不同粗糙表面。R值取5～20 μm，微凸体高度ha的范围为3～750 nm，微凸体和基体弹性模量的比值κ取1～2，为了简化计算，材料的泊松比取υa=υb=0.3，其余物理量rb，ξ可以通过公式计算得到。ξ的四个不同的取值代表了四种不同等级的接触表面。

由于所有节点和单元的接触的信息可以直接从有限元得到分析结果，所以微凸体和基体的接触参数可以很容易地确定。将从有限元分析模拟的结果同提出的新的分析模型和Hertz接触模型进行比较，以确定使用Hertz解决方案的局限性和所提出的模型的有效性。

当ξ=19.37×103m-1/2，κ=1，δ/rb=6.82×10-3时，微凸体/基体系统的应力分布图如图6(a)所示，由图可知应力主要集中在基体内，形变量主要发生在微凸体上。图6(b)既包含了Hertz模型解(实线)，新模型解(虚线)和有限元模拟的解(°点)，还包含了在不同的实际位移下，微凸体/基体系统应力应变。正如预期的那样，在这种情况下，微凸体和基体的材料相同时(κ=1)，有限元分析结果与Hertz解一致，而新提出模型预测的接触力与Hertz解不一致的原因是在第2节中讨论过的位移的双重计数的原因。

(a) δ/rb=6.82×10-3时的应力分布

(b) Hertz解、新模型解和有限元模拟解的比较

Fig.6 Stress distribution diagram and the comparison of the model results atξ=19.37×103m-1/2，κ=1

图7(a)表明，当ξ=8.66×103m-1/2，κ=1.25，δ/rb=0.011时，微凸体/基体系统的应力分布图，从图中可以看出大部分的应力依旧集中在基体内，但是还是有一部分应力出现在微凸体内。从图7(b)无量纲的位移和无量纲的接触力的关系中可以看出，当δ/rb≤0.003 5时，代表有限元分析解，Hertz解，新模型解的三条曲线基本吻合；在δ/rb≥0.003 5时，有限元分析的解开始与Hertz解开始有差别，并慢慢趋近于新提出的分析模型。很显然，在δ/rb的值较高时，新提出的模型优于Hertz模型。

(a) δ/rb=0.011时的应力分布

(b) Hertz解，新模型解，有限元模拟解的比较

Fig.7 Stress distribution diagram and the comparison of the model results atξ=8.66×103m-1/2，κ=1.25

(a) δ/rb=0.023时的应力分布

(b) Hertz解，新模型解，有限元模拟解的比较

Fig.8 Stress distribution diagram and the comparison of the model results atξ=3.36×103m-1/2，κ=1.5

图8(a)表明，当ξ=3.66×103m-1/2，κ=1.5，δ/rb=0.023时，微凸体/基体系统的应力分布图，从图中可以看出，虽然最大的应力依旧集中在基体内，但是有相当大的应力出现在微凸体内。图8(b)中的无量纲位移和无量纲接触力关系图也很好地显示出了有限元分析解的趋势，在δ/rb≤0.003时，有限元分析解，新模型解同Hertz解是基本相同的。但是δ/rb≥0.003之后，有限元分析解和新模型解开始同时偏离Hertz解，并且有限元分析解和新模型解的前进趋势是几乎相同。

图9(a)表明，当ξ=1.24×103m-1/2，κ=2，δ/rb=0.055 3时，微凸体/基体系统的应力分布图，从图中可以看出，当实际位移量最大时，最大应力同时出现在微凸体和基体中，但是从图9(b)中三张代表变形过程的图可以看出，最大应力首先出现在微凸体中，随着位移量的不断增大，最大应力慢慢向下传递到了基体中。这说明将微凸体和基体的变形分开研究对于粗糙表面接触研究非常重要。从图9(b)无量纲位移和无量纲接触力关系可知，粗糙度R=15 μm与R=5 μm、10 μm和15 μm 3种情况类似，在δ/rb值很小的时候，Hertz解和新模型解都很好的与有限元分析结果相吻合，随着δ/rb值的变大之后，有限元分析解和新模型解开始同时偏离Hertz解。

(a) δ/rb=0.055 3时的应力分布

(b) Hertz解，新模型解，有限元模拟解的比较

Fig.9 Stress distribution diagram and the comparison of the model results atξ=1.24×103m-1/2，κ=2

通过4种情况分析对比，可以发现当实际位移δ的值很小时，Hertz解和新提出的模型解都和有限元分析解吻合，所以可以推测当基体的弹性模量或者微凸体和基体的复合弹性模量很小时，Hertz模型解法可以代替新提出的模型。当实际位移δ的值不断增大时，Hertz解法就无法准确的表示真实地接触情况。新提出的模型可以有效地描述表层比基体硬的接触面的真实接触情况。

4 结 论

(1) 当粗糙表面有硬涂层或硬薄膜时，基体材料和微凸体材料的不同将会影响粗糙表面的接触特性。研究表明，基体的接触变形在微凸体接触研究中不可忽略。

(2) 基于GW模型和Hertz接触理论，分别计算微凸体和基体的变形量，从而建立微凸体和基体的刚度模型，通过耦合的刚度模型提出了考虑基体变形的微凸体弹性接触模型。

(3) 利用Hertz接触理论求得了微凸体变形量与微凸体/基体系统总变形量之间的函数关系，并应用不动点迭代法验证得到误差不超过5%的简单函数关系。实际位移量不变的情况下，随着微凸体材料不断变硬，微凸体的形变量越来越小；随着微凸体不断变小，微凸体的形变量越来越小。

(4) 利用Hertz模型和有限元分析结果验证了新模型的正确性。在变形量很小的时候，Hertz模型解和新模型解都很好的与有限元分析结果相吻合，随着变形量的变大，有限元分析解和新模型解开始同时偏离Hertz模型解，且有限元分析解和新模型解的前进趋势相同。证明新提出的接触模型适用于有硬涂层的接触表面。

(5) 对有限元模型仿真结果进行分析，表明当基体材料和微凸体材料不同时，微凸体/基体系统的应力分布会不均匀，微凸体表面的接触力比材料相同时的接触力小，最大应力比材料相同时的最大应力大。

[1] BARBER J R，CIAVARELLA M. Contact mechanics[J]. International Journal of Solids and Structures，2000，37：29-43.

[2] GREENWOOD J A，WILLIAMSON J B P. Contact of nominally flat surfaces[J]. Mathematical and Physical Sciences，1966，295(1442)：300-319.

[3] KOGUT L，ETSION I. Elastic-plastic contact analysis of a sphere and a rigid flat[J]. Journal of Applied Mechanics, 2002，69(5)：657-662.

[4] 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(2)：257-263.

[5] 赵永武，吕彦明，蒋建忠. 新的粗糙表面弹塑接触力学模型[J]. 机械工程学报，2007，43(3):95-101.

ZHAO Yongwu，LÜ Yanming，JIANG Jianzhong. New elastic-plastic model for the contact of rough surfaces[J]. Chinese Journal of Mechanical Engineering，2007，43(3)：95-101.

[6] JENG Y R，PENG S R. Elastic-plastic contact behavior considering asperity interactions for surfaces with various height distributions[J]. ASME Journal of Tribology，2006，128(2)：665-697.

[7] ALMQVIST A，SAHLIN F，LARSSON R，et al. On the dry elastic-plastic contact of nominally flat surfaces[J]. ASME Journal of Tribology，2007，40(10)：574-579.

[8] 温淑花，张学良，武美先，等. 结合面法向接触刚度分形模型建立与仿真[J]. 农业机械学报, 2009，40(11)： 197-202.

WEN Shuhua，ZHANG Xueliang，WU Meixian，et al. Fractal model and simulation of normal contact stiffness of joint interfaces and its simulation[J]. Transactions of the Chinese Society for Agricultural Machinery，2009，40(11)：197-202.

[9] 张学良，温淑花，徐格宁，等. 结合部切向接触刚度分形模型研究[J]. 应用力学学报，2003，20(1)：70-72.

ZHANG Xueliang，WEN Shuhua，XU Gening，et al. Fractal model of the tangential contact stiffness of machined surfaces in contact[J]. Chinese Journal of Applied Mechanics，2003，20(1)：70-72.

[10] 尤晋闽，陈天宁. 结合面法向动态参数的分形模型[J].西安交通大学学报，2009，43(9)：91-94.

YOU Jinmin，CHEN Tianning. Fractal model for normal dynamic parameters of joint surfaces[J]. Journal of Xi’an Jiao tong University，2009，43(9)：91-94.

[11] 尤晋闽，陈天宁. 基于分形接触理论的结合面法向接触参数预估[J]. 上海交通大学学报，2011，45(9)： 1275-1280.

YOU Jinmin，CHEN Tianning. Estimation for normal parameters of joint surfaces based on fractal theory[J]. Journal of Shanghai Jiaotong University，2011，45(9)：1275-1280.

[12] 李玲，蔡安江，蔡力钢，等. 螺栓结合面微观接触模型[J]．机械工程学报，2016，52(7): 205-212．

LI Ling，CAI Anjiang，CAI Ligang，et al．Micro-contact model of bolted-joint interface[J]. Journal of Mechanical Engineering，2016，52(7): 205-212．

[13] 肖会芳，杨荃，邵毅敏，等. 润滑状态下线接触滑动粗糙界面的动摩擦特性研究[J]. 振动与冲击，2016，35(1)：188-194.

XIAO Huifang，YANG Quan，SHAO Yimin，et al. Dynamic friction characteristics of sliding rough interfaces in line contact under lubrication[J]. Journal of Vibration and Shock，2016，35 (1): 188-194.

[14] 杨红平，傅卫平，王雯，等. 基于分形几何与接触力学理论的结合面法向接触刚度计算模型[J]. 机械工程学报，2013，49(1)：102-107.

YANG Hongping，FU Weiping，WANG Wen，et al. Calculation model of the normal contact stiffness of joints based on the fractal geometry and contact theory[J].Journal of Mechanical Engineering， 2013， 49(1)：102-107.

[15] IIDA K，ONO K. Design consideration of contact/near-contact sliders based on a rough surface contact model[J]. ASME Journal of Tribology， 2003，125(3)：1518-1525.

[16] SHI X，ANDREAS A. Investigation of contact stiffness and contact damping for magnetic storage head-disk interfaces[J]. ASME Journal of Tribology，2008，130(2)：843-854.

[17] ZHAO Y W，DAVID D M，CHANG L. An asperity micro contact model incorporating the transition from elastic deformation to fully plastic flow[J]. ASME Journal of Tribology，2000，122(2)：86-93.

[18] 田红亮，钟先友，秦红玲，等. 依据各向异性分形几何理论的固定结合部法向接触力学模型[J]. 机械工程学报，2013，49(21)：108-122.

TIAN Hongliang，ZHONG Xianyou，QIN Hongling，et al. Normal contact mechanics model of fixed joint interface adopting anisotropic fractal geometrical theory[J]. Journal of Mechanical Engineering，2013，49(21)：108-122.

[19] JOHNSON K L. Contact mechanics[M]. Cambridge: Cambridge University Press，1985.