工程实践表明，机械设备结合面接触特性对装备振动等指标有非常重要的影响^[1]。准确计算高端精密装备结合面接触特性是提高装备精度和质量的重要手段^[2]。Bowden等^[3]在1939年提出了微凸体接触理论并被广泛接受，许多研究者基于该理论建立结合面接触模型^[4-8]。至今，已有许多结合面模型被提出，如GW模型^[4]、MB模型^[7]及相关的改进型模型WA(Whitehouse&Archard)模型^[5]、CEB(Chang, Etsion & Bogy)模型^[9]、ZMC(Zhao, Maietta & Chang)模型^[10]等。然而，结合面模型仍存在理论上的争议，有的甚至存在完全相反的物理解释，如GW模型认为微凸体先弹性变形后塑性变形，而MB模型认为微凸体先塑性变形后弹性变形。

导致争论长期存在的一个重要原因是以往研究人员大多通过电导率、接触刚度^{[2, 4, 7]}等间接方法对模型进行实验验证，而直接观测结合面实际接触面变化的物理规律是验证接触模型的更准确、直接的方法。近年来，随着技术的进步，研究人员尝试采用更先进的技术直接观测结合面^[11-14]。Rubinstein等^[11]在2004年提出了一种基于光学全反射(TR)原理的方法来观测真实接触面，发现在结合面上大于临界角的入射光会穿透接触的微凸体，而在没有接触的地方进行全反射，透射的光强正比于接触的面积。文[11]研究了静摩擦到动摩擦转化过程中的前驱波的演化规律，并对法向力与真实接触面积之间的关系进行了研究，得到真实接触面积与法向力成线性关系的结论。2010年，Rubinsstein等^[14]又利用TR方法研究了摩擦黏滑运动中接触面积的演化规律。基于TR方法的一些诸多突破性发现也证明了该方法的有效性^[11-14]。2017年，Song等^[15]利用TR方法研究了PMMA材料的真实接触面积与法向力的关系，发现在轻载下，真实接触面积随法向力线性增加，而重载下增量逐步减少趋于饱和^[15]，即真实接触面积与法向力的线性关系只发生在轻载下，重载下会表现出非线性特征。然而，目前基于TR方法的真实接触面积观测实验集中于整体趋势规律的分析，很少进行接触斑点的量化分析及基于接触斑点变化的模型分析。

本文基于TR方法，构建了真实接触面积观测实验台，并采用图像处理方法分析了光学接触斑点的变化规律；基于光学图像分析结果，比较了GW和MB模型，对争议性问题进行了探讨。

1 真实接触面观测装置

当2个透明物体相互接触时，表面微凸体在压力作用下相互接触，形成连续的界面，而未能产生接触的部分被空气填充。当光线以大于临界角从PMMA射入界面时，将穿透相互接触的微凸体，而在充满空气的部分被全反射回PMMA内^[15]，理论上，透射的光强将与实际接触的微凸体面积成正比，而透射的图案将反映接触的实际状态。TR方法原理如图 1所示。

采用TR方法需要射入角界面的光线入射角大于临界角，例如，对于PMMA材料，入射角需要大于42°。进行全反射实验需要经常对入射角进行测量，给实验造成不便，为克服以往实验中需要进行全反射角测量的问题，设计采用光路如图 2所示，入射光通过一块金属镀膜反光镜射入PMMA。

$ \frac{\sin \alpha}{\sin \gamma}=\frac{n_2}{n_1}, $

(1)

$ \sin c=\frac{1}{n_2}. $

(2)

其中：α为入射光相对于PMMA的入射角，等于入射光相对于发射镜的入射角；γ为PMMA的折射角；n₁为空气折射率，约为1；n₂为PMMA折射率，约为1.49；c为全反射角，由式(2)得到约为42°。

若要形成全反射，需要γ < 90°-c，于是有

$ \sin \alpha=\frac{n_2}{n_1} \sin \gamma \leqslant \frac{n_2}{n_1} \sin \left(90^{\circ}-c\right) \approx 1.1. $

(3)

这表明，只要光线可以通过反射镜到达PMMA，就可满足全反射条件，从而避免了需要进行精确角度调整的问题。实验光路效果图如图 3所示，红色激光射在金属镀膜的反射镜上，反射光穿透PMMA投影在屏幕上，一个量角器贴在挡块上位于光路旁，实现对入射角度的测量。

基于上述光路原理，设计实验装置如图 4所示。

实验装置包含如下模块：

1) 法向加载模块：包含S形力传感器、X形法向平台和高精度步进电机。S形力传感器量程500 N；X形法向平台具有良好的直线度，可以保证力在加载范围内的直线加载；高精度步进电机可以通过输入频率来调整加载速度，提供最小2 μm的步进距离，并实现最大500 N的法向加载力。

2) 测试组件：包含上块和下块，均是PMMA。为提供足够的接触面压强，接触面尺寸选择为15.35 mm×5.00 mm，使得面压最大可到6.5 MPa。入射光经反射镜从下块入，上块出。试样表面形貌峰值高度分布如图 5所示，形貌数据由ZYGO白光干涉仪扫描得到，可见，标准差在10 ns左右，表面光滑，能有效保证透光度。

3) 真实接触面积观测模块：包含半导体激光器、相机、投影屏。半导体激光器为650 nm一字线激光器，线宽可调；相机为NIKON D610，2 400万像素；投影屏采用金属白色涂层专用光学白屏。

实验装置原理为：激光器发射一字线激光到金属反射镜上，反射光射入PMMA试样下块，入射光穿透上下块的接触界面，由上块折射出来投影在白屏上，相机连接计算机对投影屏上的图案进行采集，同时计算机对法向力进行采集。

2 图像处理方法

由于采用红色激光器，仅对红色部分灰度图进行分析。图 6是采集的到接触斑点图，展现出了较清楚的接触斑点。截取整张图像中的约1 mm²面积的2个区域进行局部分析，区域1为图 6中右上角矩形蓝色框，区域2为图中间橙色框。区域1在不同法向力作用下的接触斑点如图 7所示。可见，接触斑点明显；随着法向力的增大，接触斑点逐渐增多，变亮。仔细观察，可以看到随着法向力的增大，部分接触斑点存在融合的特征，而大部分接触斑点保持了比较好的独立性，但是边缘普遍比较模糊。设计采用计算机自适应分割算法对接触斑点图像进行分析。

文[16-19]中，普遍采用了阈值分割的方法对接触斑点图像进行处理。为此本文延续前人研究，采用基于阈值的图像分割方法。考虑到激光器本身的光强不均衡，选用局部阈值自适应分割方法。局部阈值分割方法主要有分块Otsu方法^[20]、Niblack法^[21]、Sauvola法^[22]。分块Otsu法是将整个图片分成多个小块，在这些小块中光照基本均衡，对小块进行Otsu运算，然后将各个小块拼接成原图像大小。Niblack法对每一个像素都单独计算二值化阈值，二值化阈值的确定方法由式(4)确定。该方法在比较干净的文本图像处理上效果较好^[23]，但由于Niblack法对每一个像素的都进行计算，将可能导致背景块内的像素被过度处理，产生伪影^[23]。

$ T_{x y}=\sigma_{x y}+k m_{x y}. $

(4)

其中: T_xy为阈值，σ_xy和m_xy分别为像素点(x, y)在邻域S内的标准差和均值，k∈(0, 1)为修正参数。

为改进Niblack法，Sauvola法采用了式(5)来确定阈值。该方法具有更好的鲁棒性和抗噪声性，成为了目前文本二值化领域最有效的算法之一^[24]。但是接触斑点的无规则性可能给该类算法带来重大挑战。

$ T_{x y}=m_{x y}\left[1+k\left(\frac{\sigma_{x y}}{R}-1\right)\right] . $

(5)

其中R是标准差的动态范围，当图像为8位灰度时为128。

将3种算法分别应用于区域1的图像处理，结果如图 8所示。分块Otsu处理时，采用50×150像素的网格划分；Niblack法处理时，k为0.3，模板窗大小设置为61×61像素；Sauvola法处理时，参数与Niblack法相同。

从图 8可以看出，3种方法都较好地分割出了接触斑点。其中，分块Otsu法由于分块过细，在高亮点处导致了接触斑点的分裂；Niblack法在局部亮点的连通性处理上更接近原始图片，且没有分块Otsu方法存在的边缘切割现象；Sauvola法由于平滑了背景块内的像素，导致同样参数下，接触斑点过度融合，效果相对最差。综合分析，分块Otsu法、Niblack法更适合本实验场景。

作为对比验证，将全局Otsu法、分块Otsu法、Niblack法分割得到的微凸体数目及图像归一化整体灰度(归一化到4 300)绘制于同一图上，如图 9所示。可以看到，接触斑点数目基本与整体灰度趋势一致；全局Otsu法由于不能适应局部接触区域的灰度不均匀，导致接触斑点数目偏少；分块Otsu法由于接触斑点的分裂，导致接触斑点数目比Niblack法多，但是与Niblack法的曲线趋势基本一致。总体上看，可以将微凸体变化区域分为3个区域：线性增长区Ⅰ、数目下降区Ⅱ、缓变区Ⅲ。开始加载时，接触微凸体数目增加，亮度增加，此时为线性增长区；当法向力较大时，靠近的接触斑点开始融合，形成更大的斑点，接触数目降低，但是亮度依然增加；当靠近的微凸体融合一部分后，较远的微凸体越来越难以融合，导致接触微凸体数目保持恒定，进入稳定区，同时接触面积也呈现饱和状态，增长非常微弱。为验证上述分析，本文对局部区域1的接触斑点数目进行分析。

采用MATLAB中的regionprops函数对全局Otsu法二值化后的区域1接触斑点、接触面积等进行计算并显示，如图 10所示。结合图 9分析，可以看到在法向力小于50 N时，随着法向力的增加，接触斑点显著增多，此阶段对应线性增长区Ⅰ；在法向力为50~100 N，斑点增加逐渐减少，存在明显的斑点融合现象(见红色实线框)，这导致了接触斑点数目减少，但是整体亮度提高，增速变缓，此阶段对应数目下降区Ⅱ；法向力大于100 N以后，依然会有少量的微凸体融合，但是融合速度非常慢，融合区域很小，对整体亮度贡献很小，接触微凸体数目在法向力增大到300 N的范围内，依然保持基本稳定，接触面积增长缓慢，此阶段对应缓变区Ⅲ。对区域1和2的接触斑点数目进行统计，如图 11所示，可以看到2个局部区域均能看到3个阶段的变化，但是由于局部应力分布不均，导致阶段出现的时间点并不一致。在接触理论中，认为轻载状态下主要以弹性变形为主，微凸体尚未融合，由对图 9和10的分析，本文定义法向力小于50 N时，即性增长区Ⅰ为轻载阶段；法向力不小于50 N时，即数目下降区Ⅱ、缓变区Ⅲ为重载阶段，接触面积呈现非线性增加特征。

综合上述分析，本文认为：轻载下，接触面积的线性快速增加主要是因为接触斑点数目的增多。重载下接触斑点会产生融合，接触数目减少，接触面积增加变缓，呈现非线性；接触斑点的融合会逐渐变慢变少，导致接触斑点数目保持稳定、接触面积增加缓慢。

3 接触模型对比分析

GW模型认为微凸体接触数目总量一定，正压力的增长导致接触的微凸体数目增加，引起面积线性增大^[4]，这一点和实验观测一致，GW模型的物理解释具有实际意义，但是重载段接触面积预测效果较差。MB模型认为微凸体首先会发生塑性融合，随着融合面积的增大，平均压强降低，转变为弹性变形^[7]。可见，MB模型对物理过程的解释恰好与GW模型相反。MB模型认为接触斑点数目N满足式(7)。可以看到，只要面积比率a₁/a恒定，接触面积A_a大于a接触微凸体数目N就是恒定的。

$ N\left(A_{\mathrm{a}}>a\right)=\left(\frac{a_1}{a}\right)^{D / 2} . $

(7)

其中：a_l为最大接触斑点面积，A_a是任一微凸体接触面积，a是指定的微凸体接触面积，D是分形维数。

为此，本文统计每个图像中最大接触斑点面积0.3倍以上的接触斑点数目。结果如图 12所示，可以看到：接触斑点的数目并非恒定，轻载阶段接触斑点数目急剧增加，然后经过一个短暂的下降后，接触斑点数目趋于稳定。虽然3种图像处理算法得到的接触斑点数目不一样，但是总体趋势一致，并且接触斑点数目的变化趋势与图 9整体区域接触斑点数目变化趋势基本保持一致。为此，本结果至少可以作为一个针对MB模型的特例，或者说明在某些尺度下，微凸体融合不明显，彼此独立性较好，使得MB模型并不适用于本实验。Morag等^[25]指出，MB模型错误地将微凸体的整个高度当成了真实的接触变形，从而导致得出了先塑性变形后弹性变形的错误结论。在MB模型中，微凸体轮廓表示为

$ \begin{aligned} z(x)= & G^{(D-1)} \cdot l^{(2-D)} \cdot \cos \left(\frac{\mathsf{π} x}{l}\right), \\ & -\frac{1}{2} < x < -\frac{l}{2} . \end{aligned} $

(8)

其中：l是微凸体的直径，由A_a=l²得到；x是以微凸体中心为原点的横坐标变量，G是特征长度。

MB模型错误认为微凸体的变形为微凸体中心点的高度，将微凸体的顶点高度当成了微凸体的变形量δ，得到

$ \delta=G^{(D-1)} A_{\mathrm{a}}\left(\frac{2-D}{2}\right) . $

(9)

MB模型同样沿用了GW模型中关于弹塑性变形的临界值，并得到微凸体临界变形量

$ \delta_{\mathrm{c}}=\left(\frac{H}{E}\right)^2 \cdot \frac{A_{\mathrm{a}}\left(\frac{2-D}{2}\right)}{\mathsf{π}^2 G^{(D-1)}} . $

(10)

其中：H是PMMA硬度，E是PMMA弹性模量。

当δ_c/δ=1时，化简得临界接触面积

$ a_{\mathrm{c}}=\left(\frac{G^2}{(H / 2 E)^{2 /(D-1)}}\right) \text {. } $

(11)

文[25]指出，由此MB模型导出了临界接触面积独立于尺度的结论，进而得到了先塑性变形后弹性变形的错误结论。

Morag等^[25]进一步指出，微凸体变形尺度关系如图 13所示，其中，δ是整体高度，ω是微凸体实际变形量, δ_c是临界变形量，R是微凸体半径，l_t是微凸体接触截面直径。由此，MB分形模型中的临界接触面积实际应该为

$ a_{\mathrm{c}}=\mathsf{π} R \delta_{\mathrm{c}}=\frac{1}{\mathsf{π}}\left(\frac{K H}{2 E}\right)^2\left(\frac{l^D}{G^{D-1}}\right)^2 . $

(12)

其中K为由硬度H与PMMA屈服强度σ_y的比例关系H=Kσ_y得到的比例因子。而ω可表示为

$ \omega=G^{(D-1)} l^{(2-D)}\left[1-\cos \left(\mathsf{π} l_{\mathrm{t}} / 2 l\right)\right] . $

(13)

由式(12)、(13)可见，a_c和ω是与l和l_t有关的量，也即与尺度有关的量，同时，a_c与δ_c同步变化，也即当a>a_c时，接触为塑性变形，弹塑性变形关系与GW模型保持一致。

综上所述，结合实验分析结果，本文认为GW模型对线性增长区的物理过程解释更接近真实的微凸体接触状态，而MB模型目前仍然存在理论上的争议，需要进一步研究。

4 重载下接触面积非线性问题分析

GW模型给出了轻载下较为吻合的物理解释，所预测的接触面积随法向力线性增加，但是不能有效解释重载下接触面积的非线性缓变趋向极限值的现象。按照分子力的理论，全接触将导致极大的Van der Waals力，而使得两个物体几乎连接在一起。而实际中，表面微凸体的存在，使得真实接触面积会远小于名义接触面积，全接触在实际中不可能实现。Persson等^[26]认为，微凸体上还有更小的微凸体，随着观测分辨的增大，更小的微凸体参与接触，真实接触面积会变小，这些无限小的物理的微凸体接触将导致真实的接触面积远远小于名义接触面积，而GW模型解决的只是某一个尺度的问题。但是, 即便在一定尺度上，接触面积也不可能随法向力持续增加达到全接触，也就是重载下接触面积并不会像GW模型预测的那样无限制线性增加，而应该是非线性趋向一个小于名义面积的极限值。实验观测结果也证明了该问题，重载下接触面积增加会变得非常缓慢并趋向于一个极限状态。在这个阶段，接触的微凸体数目几乎停止增长，微凸体的接触面积增长缓慢。Popov等^[27]认为存在一个饱和接触刚度，由于基体变形的影响，真实的界面接触刚度应该不会超过基体本身的刚度。其他的一些研究者也考虑了基体变形的影响，认为实际接触刚度应该考虑界面和基体两者刚度的作用^[28-29]。采用GW模型对试样结合面刚度进行计算。考虑篇幅及GW模型的广泛应用，此处仅给出预测结果，如图 14所示，结合面平均法向刚度在10⁹量级。而PMMA试样的体积刚度为EA/L=15.35×10⁶ N/m，其中E等于3 GPa，A为上下PMMA试样接触面面积76.75 mm²，L为下块高度15 mm。可见，体积刚度远小于结合面刚度。为此，一个合理的假设是，由于表面粗糙度比较小，结合面刚度非常大，远大于PMMA软性材料本身强度，在重载状态下，叠加上克服表面波纹度、微凸体融合、排除空气等需要的能量，导致界面的实际刚度快速增大，于是重载下PMMA块的变形主要以基体变形为主，基体下降的程度大于微凸体的变形，导致微凸体接触数目几乎保持不变，接触面积也很难增加。

5 结论

本文基于TR方法构建了PMMA材料真实接触面积观测实验台，采用图像处理方法，应用全局Otsu法、分块Otsu法及Niblack法3种方法对真实接触面图像中接触斑点的变化进行了分析，基于图像分析结果对主要的GW模型与MB接触模型进行了对比分析。研究结果表明：

1) 接触斑点在法向力作用下变化规律呈现3个区域：线性增长区Ⅰ、数目下降区Ⅱ、缓变区Ⅲ。线性增长区内，接触斑点数目、接触面积与法向力同步线性增加；数目下降区Ⅱ内接触斑点数目轻微下降，接触面积继续增加，增速变缓，呈现非线性；缓变区Ⅲ内接触斑点数目几乎停止增加，接触面积增加缓慢。

2) 在对接触斑点变化规律的预测上，GW模型轻载下与实验结果较为一致，但是无法有效预示重载下接触面积增加几乎停滞的现象；MB模型对于接触斑点的变化规律的解释、预测与实验结果相差较大，对MB模型仍然需要进一步的理论研究。

3) 基体变形会造成参考平面整体下降，造成参考平面分离量增大，对于软性材料，在重载下基体变形影响可能占主要成分，导致微凸体变形停滞，若能有效模拟基体变形造成的影响，将可获得更准确的接触模型。

本文仅分析了特定PMMA结合面材料，尚未进行大量实验；对于基体变形的假设仅进行了定性分析，需要进一步研究有效的数学模型。