王世军, 何花兰, 郭璞, 赵金娟

(1.西安理工大学 机械与精密仪器工程学院，陕西 西安 710048；2.西安理工大学 印刷包装工程学院，陕西 西安 710048)

零件之间相互接触的表面称为机械结合面。宏观上光滑的机械结合面，在微观上仍然是粗糙的，其宏观的载荷-变形关系与基体部分明显不同，例如机床60%～80%的刚度、90%以上的阻尼都来自结合面。要建立准确的机械整机结构模型，就必须了解结合面的静态和动态特性。结合面静态和动态特性的研究一直是机械整机性能分析的难点之一，目前仍没有一套比较成熟的理论与方法，大多数情况下仍然是靠经验估计。建立机械结合面模型，准确辨识结合面参数便成为机械结构性能分析中的一项关键技术。机械结合面研究要达到文献[1]提出的“结合面特性参数应反映结合面的基本特性，具有通用性，类似零件的特性常数如弹性模量、损耗因子，能方便地供设计解析使用”这一目标，道路依然还很漫长[2]。

粗糙表面的接触面积、接触压力分布等接触特性会直接影响结合面的接触刚度、接触热阻、配合精度、密封等结合面的性能。目前在微观粗糙表面的接触研究中提出了许多理论解析模型，这些接触研究模型主要有统计模型、分形模型和有限元模型三类。统计模型的基础是Greenwood 和Willimson建立的GW模型[3]。GW模型将粗糙峰简化为弹性半球体，并且假定粗糙峰的高度符合正态分布，利用赫兹接触理论分析接触表面的载荷-变形关系。GW模型被许多研究者改进，其中比较重要的是Chang等人在GW模型的基础上计入了塑性变形的影响后提出的CEB模型[4]。Majumdar和Bhushan认为粗糙表面具有分形特征，提出了以分形几何为基础的MB接触模型[5]。该模型采用Weierstrass-Mandelbrot分形函数描述粗糙表面的几何形貌特征，采用赫兹方法分析接触表面的载荷-变形关系。统计模型和分形模型均依赖于一定的假设条件,根据这些假设条件可以对表面轮廓特征进行简化，如微凸体都具有相同的峰顶曲率半径，其变形相互独立，微凸体之间的接触是峰对峰的正接触。这些假设与简化限制了理论解析模型的应用范围，同时影响了粗糙表面接触分析的准确性[2-6]。

文献[7]采用有限元法研究了具有几个圆柱形粗糙峰的刚性表面与弹性半无限体相互接触的情况。文献[8]采用具有半径和高度相同的一定数目圆形粗糙峰的刚性表面与弹塑性半无限体相互接触的模型来模拟多粗糙峰接触，并用有限元法对该模型进行了弹塑性分析。文献[9-10]利用由分形几何模拟的三维粗糙表面将粗糙表面的接触等效成粗糙表面与刚性平面的接触并进行了相关的有限元分析。文献[11]则认为对于具有许多任意形状微凸体的三维粗糙表面的接触问题，采用有限元法需要大量的单元，计算量过大。虽然随着计算机技术的发展,有限元计算的能力已经大幅度提高，采用有限元法求解具有强烈非线性的三维粗糙表面的接触问题仍然存在名义接触面积过小、计算规模过大的问题，例如文献[6]中三维有限元模型的名义接触面积最大为0.25×0.25 mm2，单元数高达52万，文献[12]中的名义接触面积仅为2.5×2.5 μm2，并且两个接触面只有一个是粗糙表面，另一个表面为刚性平面，单元数为23 534。一般认为，过小的名义接触面积难以反映机械加工表面的纹理特征。文献[13]直接采用表面形貌测量仪获得的粗糙表面二维数据建立具有两个粗糙接触表面的有限元弹塑性模型，模型的名义接触长度为5.6 mm，32 236个单元。与三维模型相比，二维模型接触面的名义长度可以较长，单元数量较少，而且结合面法向的载荷-变形计算结果与实验结果也能够很好地吻合。

接触问题研究中，一般认为三维有限元接触模型比二维模型更能反映真实的结合面接触状况，但是因为存在前面所述的问题，目前的研究中主要还是采用二维模型。因为二维模型的局限性，对实际接触性质的反映可能是不全面的，甚至是不正确的。一个二维模型在多大程度上反映真实的接触性质，需要从接触性质的各个方面具体考察。已有的研究表明[13]，笔者提出的这种二维模型的法向变形-载荷关系与实验结果是一致的。在本文中，笔者进一步考察了微观接触面积与载荷的关系，并与已有的电接触法[14]获得的实验结果进行了比较，更深入、广泛地考察这一模型的可行性。

1 粗糙表面的轮廓

图1是通过TR300型表面轮廓测量仪获得的磨削试样的表面轮廓曲线。表面轮廓的测量是沿试样表面垂直于磨削纹理方向扫描，测头直径2 μm,数据采样间距1 μm,采样长度5.6 mm。测量获得的数据经过最小二乘处理后,将数据点用直线连接起来形成图1所示的表面轮廓。

2 有限元模型

图2是利用图1的表面轮廓数据建立的二维有限元平面应变模型。图2(a)是两个接触试样的整体网格，图2(b)是接触表面附近的局部网格。为了尽可能减少单元数量，减少计算量，本文沿垂直方向分层逐步加密网格，以便网格疏密过渡合适。为了保证接触计算的收敛性和分析精度，在接触表面附近采用均匀的四边形网格，同时，为了保证网格水平方向的均匀性，将试样沿水平方向也进行了均匀分块，在接触表面附近获得了比较高的网格质量和合适的网格密度。整个模型的单元数为148 661，结点数为148 650。

图1 磨削表面的轮廓

图2 接触试样的有限元模型

图2(a)中下试样底部的法向位移约束为0，同时左下角的水平位移也为0，两个试样在接触表面间通过拉格朗日乘子法建立接触关系[16]。模型的法向载荷施加在上试样的顶面，最大载荷为1 MPa，按120个载荷步等步长加载。

试样材料为45钢，由单轴拉伸实验获得的真应力-应变关系如图3所示，本文将其简化为图中的多线性弹塑性应变硬化材料模型，采用Von Mises判据作为材料的塑性判据。材料弹性模量E=207 GPa，材料屈服极限σs= 320 MPa，泊松比μ=0.3。

3 应力分析的结果

图4显示了法向载荷为1 MPa时模型中的Mises应力分布，图4(a)是整体结构中的应力分布，图4(b)是图4(a)中左侧高应力区的局部应力分布。由图4(a)可以看到，接触层存在两个明显的高应力接触区，试样中的应力分布不均匀，上下试样中的应力分布相似。图5显示了上试样在不同载荷时接触区各结点的Mises应力。图中的接触结点序号是将接触区内的结点按应力值排序后的序号。由图中可以看出，随着载荷的增大，进入接触的结点数在增多，应力水平升高，同时各结点的应力水平也趋于一致，而且即便是在载荷比较低的情况下，也总是存在接触点的应力超过屈服极限的情况。从图5可以看到，随着载荷增大，弹性接触点数量减少，塑性接触点增加，即使在载荷为0.25 MPa时，也有结点的最高应力超过了材料的屈服极限，最高应力大约是试样名义应力的1 400倍。下试样的分析结果也有类似的情况。

图3 材料的应力-应变曲线

图4 1 MPa时的应力分布

图5 不同载荷时的结点应力

图6是图4(b)中上试样最高应力点由接触表面向上至试样内部路径上的Mises应力分布。由图6可以看出，高应力区只出现在很小范围内，远离接触点的试样内部的应力远低于接触点的应力，在距离接触点50 μm处，Mises应力已经降低到接触点应力的10%以下，在100 μm处则降低到5%以下。结合图4的结果，说明试样内部的应力并不均匀，接触层的应力状态实际由这些真实的接触点的应力状态决定。此外，在法向载荷低于1 MPa时，接触点附近的Mises应力分布趋势与1 MPa时相同，只是应力水平较低，高应力区的范围更小一些。

图6 接触点附近的应力分布

4 承载力的分析结果

图7是图4中上试样接触表面各接触结点在不同载荷作用下的接触力沿垂直方向分量占法向总载荷的比重，反映了接触表面上的接触结点在不同载荷下的法向承载能力，即各结点对法向载荷的贡献大小。图7中横轴是接触中的结点的排序号，按结点法向力比重从大到小排列。纵轴是接触结点垂直分力占试样法向载荷的比重。由图中可以看到，在法向载荷较低的情况下，接触点数很少，各结点承载力差别很大。随着载荷增大，进入接触的结点数增多，承载力趋于一致，这与图5的Mises应力变化趋势一致，但与图5不同的是，不同载荷条件下总是存在少数承载力很低的结点，在法向载荷升高时，这类结点所占的比重在减小。结合图5的结果可以发现，在载荷为1 MPa时，接触点的承载能力趋于一致的同时，接触点都处在塑性应力状态。

图7 不同载荷时结点的法向承载力的变化

由于轮廓测量是等间距采样，进入接触的结点数反映了法向接触面积的大小，即真实承载面积的大小。图7中法向载荷为1 MPa时，接触结点只有7个，占上试样表面结点数的0.13%，真实接触面积只占名义接触面积很小的一部分，使得图5中接触结点的应力水平远高于试样的名义应力，真实接触区的应力水平和应力状态与远离接触点的试样内部有很大不同，接触层的性质与试样内部明显不同。

5 法向变形与法向载荷的关系

图8给出了计算获得的法向位移与载荷关系曲线，法向的位移取自图4(a)中上试样的顶面中心点。图中的试验曲线采用文献[13]中的试验方法获得，试验中的试样采用环形试样，环形接触面的外径为45 mm，内径为25 mm，试样高度为50 mm。计算获得的位移与载荷关系大体呈线性，与图中的试验曲线比较，计算结果与试验结果能较好地吻合，表明模型的宏观计算结果是可靠的。

图8 法向位移与载荷的关系

图9是载荷为1 MPa时图4(b)中接触点附近上试样沿垂直方向路径上的法向变形，是路径上各点相对于路径起点的相对位移，反映了接触点附近垂直方向的变形程度。由图9可以看出，接触点附近的变形主要集中在距离接触点0～10 μm的范围内，这个结果与图6所示的路径上的应力分布规律一致。

图9 接触点附近的法向变形

6 接触面积与载荷的关系

图10给出了有限元计算获得的接触面积比与法向载荷的关系：随着载荷的升高，真实接触面积大体按线性规律增大，载荷为1 MPa时，真实接触区只占名义接触面积的0.13%。作为对照，图中也给出了文献[14]中采用接触电阻法获得的真实接触面积的测量结果。由于接触问题的特殊性，直接观察和测量真实接触区的接触状态和接触面积比较困难。接触电阻法是通过测量接触表面的接触电阻估算真实接触面积的一种方法，该方法认为真实接触面积与接触电阻成正比，通过测量试样间的接触电阻即可确定真实的接触面积。文献[14]中给出的测量结果是采用磨削过的、面积为21 cm2钢平面的测量结果。由图中可见，结果测量与计算结果一致，说明真实接触面积的有限元计算结果是可靠的。

图10 有限元法与电阻法的面积比较

接触区数目的增多和接触区面积的增加都可以导致接触面上真实接触面积增加。结合图4、5的进一步分析表明，在研究的载荷范围内，真实接触面积的增大并不是由于接触区数目的增多，而是微观接触区面积的增大。

7 真实接触区的高度分布

图11为图4中上试样接触表面轮廓高度的概率密度分布图，图12为对应的累积概率图。由图11可以看出，磨削表面轮廓大体符合正态分布，但是也存在一定出入。由图12可以看出，在轮廓的平均位置附近，轮廓高度分布比较接近理想的正态分布，但是在曲线的两端，即轮廓高度的两端与理想的正态分布还是有一定差距。

图11 轮廓高度的概率密度分布

图12 轮廓高度的累积概率图

图12中标出了上试样中接触点的高度在累积概率图中的位置。由图中可以看出，接触点都处于平均高度以上的位置而不是散布在整个轮廓高度范围内，轮廓平均高度以下的位置不产生接触，下试样以及5组同一批次的试样都有同样的特点。由于真实接触点出现在轮廓高度的上部位置，而上部位置的轮廓高度与正态分布的相对偏差通常比中部大，这使得基于轮廓特征符合正态分布假定的GW模型的计算结果与实际情况相比具有比较大的误差。此外，由于真实接触点的数量很少并且不出现在轮廓高度的平均值附近以及平均值以下的位置，这些位置上轮廓数据的分布规律不会影响接触的状态和性质，图11中轮廓高度分布的左侧和中间的形状即使不满足正态分布，接触层也可以有相同的接触性质，如法向变形和刚度等。这一观点也通过下面的计算得到了验证：通过修改图2中表面轮廓高度中平均位置以下的数据点的高度值以后建立新的有限元模型，其接触区的应力分布与宏观的变形和修改之前的图2中的模型并没有不同。当然，这一结论的正确性和适用范围尚需更多、更广泛的计算和试验结果的验证以及理论分析。

8 结 论

采用粗糙表面轮廓数据建立的二维平面应变接触模型不仅能够准确反映接触层宏观的法向位移-载荷关系，由其计算获得的真实接触面积与载荷的关系也能够跟已有的接触电阻法的测量结果吻合，说明该模型的微观计算结果也具有一定的可靠性，可以用作结合面微观接触状态的观察和接触性质研究的较为直观的工具。

相对于接触问题的三维有限元模型，本文的二维模型在规模和效率上都更有优势。图2中的二维模型有44.6万自由度，如果采用矩形的三维模型计算，接触面尺寸为5.6×5.6 mm2，采用同样的网格密度，自由度数将超过22.5亿。目前有限元软件的最大处理能力只有2亿自由度，常用的有限元软件如ANSYS，其商业版本目前的处理能力只有2千万自由度。因此，现阶段采用三维有限元模型求解较大尺寸的接触问题仍然存在困难，相比之下，二维模型还有较大的自由度空间。

现有的计算结果表明，真实接触面积很少，载荷为1 MPa时仅占名义接触面积的0.13%，并且真实接触点只出现在轮廓高度平均值以上的位置，该部分轮廓高度的分布与正态分布有一定差异，轮廓高度平均位置附近及其以下的高度分布不影响接触性质，表面轮廓符合正态分布应该不是GW模型的必要条件。

当法向载荷在1 MPa以下时，粗糙表面接触面积随载荷的增加而线性增大，接触区的法向承载力随载荷增大趋于一致。在很小的法向载荷下，接触区就出现塑性，可以认为接触状态下，接触区一直存在塑性。随着载荷的增加，塑性区占的比重逐渐增大。

参考文献：

[1]黄玉美,付卫平,董立新,等.结合面法向动态特性参数研究[J].机械工程学报,1993,29(3): 74-78.

Huang Yumei, Fu Weiping, Dong Lixin, et al.Research on the dynamic normal characteristic parameters of joint surface[J].Journal of Mechanical Engineering, 1993, 29(3): 74-78.

[2]田红亮,朱大林,秦红玲.机械结合面接触模型的研究进展及存在的问题[J].机械工程与技术, 2013, (2): 1-10.

Tian Hongliang,Zhu Dalin,Qin Hongling.Advances and existing issues in contact model of mechanical joint surface[J].Mechanical Engineering and Technology, 2013,(2):1-10.

[3]Greenwood J A, Williams J B.Contact of nominally flat surfaces[J].P.Roy.Soc.Lond.: A Mat., 1966, 295(1442): 300-319.

[4]Chang W R,Etsion I, Bogy D B.An elastic-plastic model for the contact of rough surfaces[J].ASME J.Tribol, 1987, 109(2):257-263.

[5]Majumdar A, Bhushan B.Fractal model of elastic-plastic contact between rough surfaces[J].ASME J.Tribol, 1991, 113(1):1-11.

[6]杨国庆,熊美华,洪军,等.3D粗糙表面的数字化表征与接触特性分析[J].西安交通大学学报, 2012,46(11): 58-63.

Yang Guoqing, Xiong Meihua, Hong Jun,et al.Numerical charecterization and contact performances for 3D rough surfaces[J].Journal of Xi’an Jiaotong University, 2012,46(11):58-63.

[7]Komvopoulos K，Choi D H.Elastic finite element analysis of multi-asperity contact[J].ASME Journal of Tribology,1992，114:823-831.

[8]杨楠，陈大融，孔宪梅.多粗糙峰弹塑性接触的有限元分析[J].摩擦学学报, 2000, 20(3): 202-206.

Yang Nan, Chen Darong, Kong Xianmei.Elastic-plastic finite element analysis of multi-asperity contacts[J].Tribology, 2000, 20(3): 202-206.

[9]Sahoo P, Ghosh N.Finite element contact analysis of fractal surfaces[J].J.Phys: D Appl.Phys., 2007, 40(14): 4245-4252.

[10]Pei L, Hyun S, Molinari J F, et al.Finite element modeling of elasto-plastic contact between rough surfaces[J].J.Mech.Phys.Solids,2005,53(11):2385-2409.

[11]Bhushan B.Contact mechanics of rough surfaces in tribology：Multiple asperity contact[J].Tribology Letters，1998，4:1-35.

[12]王霄,梁春,刘会霞,等.车削真实粗糙表面的弹塑性接触有限元分析[J].润滑与密封,2008,33(12):72-74.

Wang Xiao, Liang Chun, Liu Huixia, et al.Finite element analysis based on elastoplastic contact of real turning rough surface[J].Lubrication Engineering, 2008,33(12):72-74.

[13]王世军,赵金娟,张慧军,等.一种结合部法向刚度的预估方法[J].机械工程学报, 2011, 47(21):111-115.

Wang Shijun, Zhao Jinjuan, Zhang Huijun, et al.A method of estimating normal stiffness of joint[J].Journal of Mechanical Engineering, 2011, 47(21): 111-115.

[14]鲍登F P,泰伯D.固体的摩擦与润滑[M].陈绍澧,袁汉昌,丁雪加,译.北京:机械工业出版社, 1982:35-37.

[15]Poulios K, Klit P.Implementation and applications of a finite-element model for the contact between rough surfaces[J].Wear, 2013, 303:1-8.

[16]Bathe K J, Bouzinov P A.On the constraint function method for contact problems[J].Computers & Structures, 1997,64(6)：1069-1085.