DOI: 10.11817/j.issn.1672-7207.2019.06.011

刚度连续、单调且光滑变化的粗糙界面法向弹塑性接触模型

肖会芳¹，孙韵韵²，徐金梧²

(1. 北京科技大学 机械工程学院，北京，100083；

2. 北京科技大学 钢铁共性技术协同创新中心，北京，100083)

摘要：提出一种新的基于低阶椭圆曲线方程的微凸体法向弹塑性接触刚度模型。进一步基于粗糙表面形貌的Greenwood-Williamson统计模型和微凸体高度分布概率密度函数，推导整个粗糙界面接触刚度和接触载荷表达式，建立粗糙界面接触的总刚度模型，并将模型计算结果与实验测量结果进行对比分析。研究结果表明：该模型考虑粗糙界面的弹性变形、弹塑性变形和完全塑性变形接触状态，且同时满足不同接触状态之间微凸体的接触刚度、接触载荷和接触面积均连续、单调且光滑变化的条件，克服了以往模型存在的接触刚度非连续、非单调的缺点；本文建立的粗糙界面接触刚度模型有效。

关键词：粗糙界面；接触刚度；弹塑性变形；统计学模型

中图分类号：TH11              文献标志码：A         文章编号：1672-7207(2019)06-1343-08

An elastoplastic contact model for rough surface contact with continuous, monotonous and smooth stiffness

XIAO Huifang¹, SUN Yunyun², XU Jinwu²

(1. School of Mechanical Engineering, University of Science and Technology Beijing, Beijing 100083, China;

2. Collaborative Innovation Center of Steel Technology,University of Science and Technology Beijing, Beijing 100083, China)

Abstract: A new elastoplastic contact stiffness model based on the low order elliptic curve was proposed for single asperity contact. The Greenwood-Williamson model was used for surface topography description and the whole contact stiffness for the rough surface was further derived based on the height distribution density function for asperities, and the whole contact stiffness model was established for the rough interface contact. The model predictions were compared with the experimental results. The results show that the proposed model can be used to describe the different contact status of elastic, elstoplastic and fully plastic contact and make the contact parameters of stiffness, force and area all continuous, monotonous and smooth between different contact status, which overcomes the shortcomings of non-continuous and non-monotonous in previous models, the proposed model has validity.

Key words: rough surface; contact stiffness; elastoplastic deformation; statistical model

粗糙接触界面广泛存在于机械结构与机器设备中，例如齿轮的轮齿啮合界面、轴承的滚动体-滚道界面、机床螺栓联接界面、高速列车的车轮-轨道形成的轮轨界面等^[1-4]。粗糙界面接触刚度是机械结构整体刚度的重要组成部分，也是描述界面特性最重要的参数之一，其变化直接影响界面以及机械装备系统的静态特性和动力学特性，包括接触压力分布、振动响应特性、磨损特性以及工作稳定性等^[5-7]。因此，粗糙界面法向接触刚度的准确计算对机械结构与系统的性能分析与预测至关重要。在微观尺度下，粗糙表面可以看成是具有一定曲率半径的单个微凸体按照不同的分布特征形成的整体。随着外部载荷增加，微凸体的变形量增大，单个微凸体的变形状态经历弹性变形、弹塑性变形和完全塑性变形。其中，弹性变形阶段的接触特征通常采用Hertizan接触模型表征^[8]；完全塑性变形阶段的接触特征通常采用线性硬度模型表征^[9]。而处于中间过渡区的弹塑性变形阶段接触模型，需要保证微凸体的接触状态变量，包括接触压力、接触载荷和接触刚度等，从弹性阶段单调、连续且光滑过渡到完全塑性阶段。为了满足不同接触状态之间单调、连续且光滑的条件，研究人员提出了不同的弹塑性模型来表征界面的法向接触特性。ZHAO等^[10]基于CHANG等^[11]的弹-塑性模型，利用插值方法建立了弹塑性变形状态下实际接触面积、平均接触压力与法向接近量之间的关系，提出了一种粗糙界面的弹塑性接触模型(ZMC模型)。但是，平均接触压力和接触载荷的变化在临界点处不光滑，导致接触刚度在临界点处不连续，且刚度变化不是单调的。赵永武等^[12]进一步对ZMC模型进行改进，基于三次样板函数插值的方法，提出了一种在临界点满足连续性和光滑性要求的法向接触模型(MZMC模型)。虽然接触刚度随法向变形量的变化关系是连续的，但是不光滑和非单调。BRAKE^[13]采用Hermit多项式插值的方法，提出了一种满足不同接触状态之间光滑和连续过渡的粗糙表面弹塑性接触模型。但是，由于Hermit多项式的阶次较高，插值获得的弹塑性阶段平均接触压力与法向变形量之间的关系曲线出现振荡，使得接触压力和接触刚度的变化过程并非单调。徐超等^[14]提出一种利用低阶椭圆曲线插值进行解析建模的方法，建立了粗糙界面弹塑性阶段接触模型，该模型能够描述微凸体的接触压力和接触面积随法向接近量单调、连续和光滑的变化过程，但是，界面的法向接触刚度仍然具有非单调的缺点。针对上述弹塑性接触模型的不足，本文作者提出一种新的粗糙界面法向接触弹塑性模型，该模型能够描述微凸体的接触状态变量，包括接触载荷、接触面积和接触刚度，在不同接触状态之间连续、光滑且单调变化。在此基础上，基于粗糙界面微凸体高度分布概率密度函数，推导获得整个粗糙界面的法向接触刚度表达式。

1  单个微凸体接触刚度模型

单个微凸体与刚性平面的接触示意图如图1所示。表面粗糙体的形状为球体且各粗糙体的曲率半径均为β，单个微凸体的高度z定义为微凸体顶点与微凸体高度均线间的距离。微凸体均线与粗糙面均线之间的距离为d_d，d_d=1.15σ_s^[15]，σ_s为粗糙微凸体高度分布的标准偏差。

图1  单个微凸体与刚性平面接触示意图

Fig. 1  Geometry of single asperity contact

1.1  弹性接触(w＜w₁)

根据Hertzian弹性接触理论，其弹性接触力可以表示为^[8]

          (1)

式中：E为接触界面的等效弹性模量；w为弹性变形量，可以表示为

              (2)

接触面积可以表示为

             (3)

平均接触压力可以表示为

             (4)

单个微凸体的法向接触刚度可以表示为

         (5)

由式(5)可知，弹性变形阶段，单个微凸体的法向接触刚度与变形量呈非线性变化关系。微凸体从弹性变形过渡到弹塑性变形的临界变形量w₁^[16]为

             (6)

式中：H为材料硬度；K为硬度系数，与材料的泊松比相关，K=0.454+0.41v。

1.2  完全塑性接触(w＞w₂)

当微凸体的法向变形量w＞w₂时，微凸体发生完全塑性变形，其临界变形量w₂^[17]为

               (7)

此时，界面的平均接触压力p_p为材料的硬度H，即

                (8)

接触面积为

             (9)

界面的接触力可以表示为

           (10)

微凸体的法向接触刚度k_p可以表示为

               (11)

由式(11)可知，完全塑性变形阶段，单个微凸体的法向接触刚度为常数，与变形量无关。

1.3  弹塑性接触(w₁＜w＜w₂)

当微凸体的变形量w₁＜w＜w₂时，界面处于弹塑性变形阶段。随着微凸体变形量增大，微凸体从弹性变形过渡到弹塑性变形，直至完全塑性变形阶段的整个过程，接触状态参量应该满足连续、光滑且单调的条件。为了满足不同接触状态之间单调、连续且光滑的条件，研究人员提出了不同的弹塑性模型，以表征界面的法向接触特性。

1.3.1  ZMC模型

根据ZMC弹塑性接触模型，弹塑性变形阶段的平均接触压力为^[10]

    (12)

接触面积为^[10]

 (13)

联立式(12)和式(13)，接触载荷可以表示为

      (14)

单个微凸体的法向接触刚度可以表示为

         (15)

1.3.2  改进的ZMC模型(MZMC模型)

在改进的ZMC模型中，弹塑性变形阶段的平均接触压力为^[12]

 (16)

其中：

接触面积为^[12]

        (17)

联立式(16)和式(17)，接触载荷可以表示为

       (18)

单个微凸体的法向接触刚度可以表示为

   (19)

1.3.3  Brake模型

根据Brake模型，弹塑性变形阶段的平均接触压力为^[13]

  (20)

其中：

接触面积为^[13]

 (21)

则接触刚度可以表示为

(22)

1.3.4  Xu模型

根据Xu模型，弹塑性变形阶段平均接触压力^[14]为

   (23)

其中：

接触面积^[14]为

 (24)

则接触刚度可以表示为

           (25)

1.3.5  本文模型

为了满足单个微凸体接触刚度在弹性、弹塑性和完全塑性的不同接触之间变化的连续、光滑和单调性要求，本文提出采用低阶椭圆曲线来表征微凸体弹塑性阶段接触刚度，并通过对刚度表达式积分，获得微凸体弹塑性阶段载荷关系式的弹塑性接触模型。

单个微凸体接触刚度的表达式为

        (26)

其中：C₀和(k_p-C₁)分别为椭圆函数的长半轴和短半轴，k_p=2πβH，C₀和C₁为待定系数，由弹性变形到弹塑性变形临界变形量w₁处刚度的连续性和光滑性条件确定。

      (27)

对式(26)求导，可得

      (28)

将连续性和光滑性边界条件式(27)代入式(26)和(28)，则有

，

        (29)

单个微凸体的接触刚度可以表示为

   (30)

因此，弹塑性变形阶段的载荷函数为

  (31)

其中：C₂为待定系数，由弹性变形到弹塑性变形临界变形量w₁处载荷的连续性条件确定，即

     (32)

将式(32)代入式(31)，有

  (33)

单个微凸体的接触面积表达式采用三次样板函数描述，表达式为

  (34)

其中：C₃和C₄为待定系数，由临界变形量w₁和w₂处面积的连续性条件确定，

，

         (35)

将式(35)代入式(34)，有

    (36)

1.4  不同模型计算结果分析

采用不同弹塑性模型计算获得的单个微凸体的接触面积、接触载荷和接触刚度随法向变形量的变化关系曲线，分别如图2～4所示，计算参数为E=52 GPa，v=0.25，H=6.2 GPa，β=0.182 mm。由图2可知：对5种不同的弹塑性接触模型，在整个弹性、弹塑性和完全塑性变形状态，接触面积随法向变形量的变化关系曲线均为单调、连续且光滑的。

图2  不同模型计算获得的单个微凸体接触面积

Fig. 2  Contact area of single asperity contact for different models

接触载荷的变化关系曲线却并非均为单调、连续且光滑，如图3所示。对ZMC模型，接触载荷随法向变形量的变化曲线单调且连续，但是在临界变形量w₂处不光滑，由此导致接触刚度在临界变形量w₂处不连续，如图4所示。对MZMC模型和Brake模型，由于接触载荷表达式多项式的阶次较高，使得接触载荷随法向变形量的变化曲线非单调，导致接触刚度随法向变形量的变化出现震荡，甚至出现负刚度。对Xu模型，接触载荷随法向变形量的变化曲线单调、连续且光滑，但是其接触刚度随法向变形量的变化曲线却非单调。对本文提出的模型，接触载荷和接触刚度随法向变形量的变化曲线均为单调、连续且光滑变化。

图3  不同模型计算获得的单个微凸体接触载荷对比

Fig. 3  Contact force of single asperity contact for different models

图4  不同模型计算获得的单个微凸体接触刚度

Fig. 4  Contact stiffness of single asperity contact

2  粗糙界面接触刚度模型

采用GW统计学模型计算粗糙界面上接触微凸体数^[18]，并将单个微凸体接触刚度和载荷对微凸体数进行积分，获得粗糙界面接触的总刚度和总载荷。粗糙界面接触模型示意图如图5所示。

图5  粗糙表面与刚性平面的统计接触模型

Fig. 5  Statistical micro-contact model of rough surface

引入如下量纲一的量：

当微凸体高度zn大于粗糙面均线与刚性平面之间的距离h_n-d_n时，微凸体与刚性平面接触，其接触概率为^[18]

      (37)

其中：，是高度为标准高斯分布的概率密度函数。则整个粗糙界面的接触刚度和总载荷可以分别表示为

     (38)

     (39)

式中：N为粗糙表面微凸体总数，N=nA_n；n为微凸体分布密度；A_n为名义接触面积。

将不同变形阶段单个微凸体的刚度表达式和载荷表达式分别代入式(38)和(39)，则粗糙界面的接触刚度和载荷可以分别表示为

      (40)

  (41)

其中：。

3  粗糙界面接触刚度模型验证与结果分析

3.1  模型验证

文献[19]采用超声反射系数法测量了不同接触面积钛合金Ti-6Al-4V粗糙界面的接触刚度，钛合金的材料参数为E=115 GPa, v=0.31，界面的微凸体高度标准偏差σ_s=2.3 μm。采用本文接触刚度模型式(40)和式(41)计算获得的法向固体接触刚度-载荷变化关系曲线与文献[19]实验测量获得的结果对比图，如图6所示。由图6可知，模型的计算结果与实验测试结果基本一致，说明本文的界面固体接触刚度的模型是有效的、可靠的。

图6  粗糙界面固体接触刚度的模型计算结果与实验结果^[19]对比

Fig. 6  Comparison of solid contact stiffness versus pressure between model results and experimental results^[19] for different nominal area

3.2  不同粗糙度界面刚度

光滑界面、中等粗糙界面和粗糙界面的表面形貌参数^[20]如表1所示。对不同粗糙度界面，采用本文提出的弹性-弹塑性-完全塑性刚度模型式(40)和式(41)计算获得的界面接触总刚度随接触载荷的变化关系曲线如图7所示，并与纯弹性接触GW模型计算获得的刚度-载荷关系曲线进行对比。

表1  不同粗糙度表面形貌参数

Table 1  Parameters of surfaces with different roughness

图7  不同粗糙度界面的法向接触刚度随载荷的变化关系曲线

Fig. 7  Variation of normal contact stiffness with normal load for different roughness

由图7可知：考虑弹性-弹塑性-完全塑性不同接触状态刚度模型计算获得的刚度大于纯弹性GW接触模型的刚度计算值，且随着表面粗糙度增大，2种不同模型获得的刚度差值逐渐增大。这是由于随着表面粗糙度增大，界面从弹性变形过渡到弹塑性变形的临界变形量w₁逐渐减小，弹塑性变形阶段对界面接触状态的影响逐渐增大，导致2种模型获得的刚度差值逐渐增大。

由图7还可知：粗糙界面接触刚度随法向载荷呈非线性递增，且增大的幅度随载荷的增大逐渐减小(曲线的斜率逐渐减小)；随着表面粗糙度增大，法向接触刚度逐渐减小。

4  结论

1) 提出一种新的基于低阶椭圆曲线的微凸体法向弹塑性接触刚度模型，该模型考虑粗糙界面的弹性变形、弹塑性变形和完全塑性变形接触状态，且保证不同接触状态之间微凸体的接触刚度、接触载荷和接触面积均连续、单调且光滑变化，克服了以往模型的接触刚度非连续、非单调的缺点。

2) 基于微凸体高度分布概率密度函数，推导获得了整个粗糙界面接触刚度和接触载荷表达式，建立了整个粗糙界面接触刚度模型。

3) 模型计算结果与实验测量结果基本一致，验证了本文建立的粗糙界面接触刚度模型的有效性。

4) 考虑弹性-弹塑性-完全塑性不同接触状态刚度模型计算获得的刚度大于纯弹性接触GW模型刚度计算值，且随着表面粗糙度增大，2种模型获得的刚度差值递增。

5) 粗糙界面接触刚度随法向载荷呈非线性递增，随着表面粗糙度的增大而减小。

参考文献：

[1] SAXENA A, PAREY A, CHOUKSEY M. Time varying mesh stiffness calculation of spur gear pair considering sliding friction and spalling defects[J]. Engineering Failure Analysis, 2016, 70: 200-211.

[2] LIU Jing, SHAO Yimin. Dynamic modeling for rigid rotor bearing systems with a localized defect considering additional deformations at the sharp edges[J]. Journal of Sound and Vibration, 2017, 398: 84-102.

[3] XIAO Huifang, SHAO Yimin, XU Jinwu. Investigation into the energy dissipation of a lap joint using the one-dimensional microslip friction model[J]. European Journal of Mechanics A: Solids, 2014, 43: 1-8.

[4] LUNDBERG O E, NORDBORG A, LOPEI ARTEAGA I. The influence of surface roughness on the contact stiffness and the contact filter effect in nonlinear wheel-track interaction[J]. Journal of Sound and Vibration, 2016, 366: 429-446.

[5] CARRELLA A, BRENNAN M J, WATERS T P, et al. Force and displacement transmissibility of a nonlinear isolator with high-static-low-dynamic-stiffness[J].International Journal of Mechanical Sciences, 2012, 55(1): 22-29.

[6] ZOU H T, WANG B L. Investigation of the contact stiffness variation of linear rolling guides due to the effects of friction and wear during operation[J]. Tribology International, 2015, 92: 472-484.

[7] 李小彭, 潘五九, 高建卓, 等. 结合面形貌特性对模态耦合不稳定系统的影响[J]. 机械工程学报, 2017, 53 (5): 116-127.

LI Xiaopeng, PAN Wujiu, GAO Jianzhuo, et al. Influence of surface topography characteristics on mode coupling instability system[J]. Journal of Mechanical Engineering, 2017, 53(5): 116-127.

[8] JOHNSON K L. Contact mechanics[M]. Cambridge: Cambridge University Press, 1985.

[9] LIU Peng, ZHAO Han, HUANG Kang, et al. Research on normal contact stiffness of rough surface considering friction based on fractal theory[J]. Applied Surface Science, 2015, 349: 43-48.

[10] ZHAO Y W, MAIETTA D M, CHANG L. An asperity microcontact model incorporating the transition from elastic deformation to fully plastic flow[J]. Journal of Tribology, 2000, 122(1): 86-93.

[11] 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.

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

ZHAO Yongwu, LU Yanming, JIANG Jianzhong. New elastic-plastic model for the contact of rough surfaces[J]. Journal of Mechanical Engineering, 2007, 43(3): 95-101.

[13] BRAKE M R. An analytical elastic-perfectly plastic contact model[J]. International Journal of Solids and Structures, 2012, 49(22): 3129-3141.

[14] 徐超, 王东. 一种改进的粗糙表面法向弹塑性接触解析模型[J]. 西安交通大学学报, 2014, 48(11): 115-121.

XU Chao, WANG Dong. An improved analytical model for normal elastic-plastic contact of rough surfaces[J]. Journal of Xi’an JiaoTong University, 2014, 48(11): 115-121.

[15] SOJOUDI H, KHONSARI M M. On the modeling of quasi-steady and unsteady dynamic friction in sliding lubricated line contact[J]. Journal of Tribology, 2010, 132(1): 012101.

[16] LIOU J L, LIN J F. A modified fractal microcontact model developed for asperity heights with variable morphology parameters [J]. Wear, 2010, 268(1/2): 133-144.

[17] 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.

[18] GREENWOOD J A, WILLIAMSON J B P. Contact of nominally flat surfaces[J]. Proceedings of the Royal Society of London, Series A, 1966, 295: 300-319.

[19] MULVIHILL D M, BRUNSKILL H, KARTAL M E, et al. A comparison of contact stiffness measurements obtained by the digital image correlation and ultrasound techniques[J]. Experimental Mechanics, 2013, 53(7): 1245-1263.

[20] BEHESHTI A, KHONSARI M M. Asperity micro-contact models as applied to the deformation of rough line contact[J]. Tribology International, 2012, 52: 61-74.

(编辑  赵俊)

收稿日期：2018-06-07；修回日期：2018-08-07

基金项目(Foundation item)：国家自然科学基金资助项目(51775037)(Project(51775037) supported by the National Natural Science Foundation of China)

通信作者：肖会芳，博士，副教授，从事机械结构界面动力学研究；E-mail：huifangxiao@ustb.edu.cn