通过对式(16)求关于r的一阶导数，有

式中：H1与Y1分别为1阶斯特鲁夫函数与第二类1阶的贝塞尔函数。

由于式(18)为奇异积分，难以获得它的解析解，故本文采用高斯-切比雪夫求积公式对单个微凸体模型进行数值求解，详细的求解步骤如下，图3给出了数值求解流程步骤为：

步骤1：先给定一个法向载荷F的大小。

步骤2：合理假设载荷F作用下接触半径a的最大取值amax与最小值amin，并采用二分法取其初始值a0=(amax+amin)/2。

步骤3：使用高斯-切比雪夫求积公式将式(18)及其约束方程式(13)化简成积分区间为[-1,1]的高斯切比雪夫型求积方程，然后将a0代入方程中，并使用MATLAB进行数值迭代求解。

步骤4：通过“步骤3”的求解即可获得接触区域的压力分布p(t)，并判断接触边缘处的压力是否为有限的正值，同时需要保证接触区域内的压力p(t)是光滑连续变化的，否则继续执行“步骤2”与“步骤3”，直至获得合理的压力分布为止。

步骤5：类似地，采用高斯-切比雪夫求积公式对式(17)进行化简，然后将上述获得的真实接触半径a以及压力分布p(t)代入方程，即可求解出弹性半空间的压入深度ω。

实际工程中不同粗糙表面对应的微凸体曲率半径各不相同，在上述的数值求解过程中，以微凸体曲率半径R=10 mm为例来研究表面张力对接触特性的影响，其他参数G=1 MPa, μ=0.4, β=0.1 J/m2。获得了无量纲接触载荷F/(4/3 E* R1/2 ω3/2)和无量纲接触面积A/πRω与无量纲参数(Rω)1/2/s的数值关系，通过对数值解的拟合，得到了考虑表面张力时载荷F和真实接触面积A与压入深度ω的关系式分别如式(19)和式(20)所示，且不同的R值只会影响式(19)和式(20)中的拟合系数。

图4为无量纲接触载荷F/(4/3 E* R1/2 ω3/2)与无量纲参数(Rω)1/2/s的关系。由图4可知，式(19)与数值结果吻合，说明了考虑表面张力时的接触载荷与压入深度关系模型的正确性。同时发现，随着(Rω)1/2/s的减小，F/(4/3 E* R1/2 ω3/2)逐渐增大。

图5为无量纲真实接触面积A/πRω与无量纲参数(Rω)1/2/s的关系。类似地，数值结果与式(20)吻合，说明了考虑表面张力时的真实接触面积与压入深度关系模型的正确性。此外，随着(Rω)1/2/s的减小，A/πRω也逐渐减小。

同时，从图4和图5可知，当不计表面张力的影响时，即s趋于无穷小，参数(Rω)1/2/s趋于无穷大时，本模型逐渐趋近于赫兹接触模型；相反，当考虑表面张力的影响时，接触载荷与真实接触面积明显不同于赫兹接触模型的分析结果；当压入深度相同时，与赫兹求解接触模型的分析结果相比，由于表面张力的存在，接触载荷较大，真实接触面积较小。

根据刚度的定义，可得单个微凸体与刚性平面接触时的法向接触刚度为

图6显示了无量纲接触刚度与无量纲参数(Rω)1/2/s的关系曲线。由图6可知，数值结果与式(21)吻合。当忽略表面张力的影响时，即参数(Rω)1/2/s趋于无穷大，接触刚度逐渐接近于赫兹接触分析结果；相反，当表面张力的影响不可忽略时，此时接触刚度明显偏离于赫兹结果。

3 结合面新模型

受外界载荷作用，微凸体与刚性平面产生接触，单个微凸体的变形量ω=z-d。采用统计学方法并结合微凸体高度与曲率分布的概率密度函数式(10)，将单个微凸体的计算模型扩展到整个结合面上，获得了结合面的接触载荷Ft、真实接触面积At、接触刚度Kt分别为

将式(19)、式(20)和式(21)分别代入式(22)、式(23)和式(24)，并进行无量纲处理有

F_t* = (√m₂)/(π^(3/2)) √(α^(5/2)/(6(α-1))) ∫∫ (λ*-u)^(3/2) ρ^(5/2) [1+0.8102(√((λ*-u)/ρ)/s*)^(-0.946)] * exp[ (3/2)ρ² - (αλ*²)/(2(α-1)) ] * erfc( 3ηρ - (η√α)/(α-1) λ* ) dλ* dρ

(26) 无量纲化后的总真实接触面积表达式：

A_t* = √(3α²/(32π(α-1))) ∫∫ (λ*-u) ρ² [ (1+0.2181(√((λ*-u)/ρ)/s*)^(-0.8505)) / (1+1.355(√((λ*-u)/ρ)/s*)^(-0.7641)) ] * exp[ (3/2)ρ² - (αλ*²)/(2(α-1)) ] * erfc( 3ηρ - (η√α)/(α-1) λ* ) dλ* dρ

(27) 无量纲化后的总接触刚度表达式：

K_t* = 1/(√m₂ π^(3/2)) √(3α^(3/2)/(8(α-1))) ∫∫ (λ*-u)^(1/2) ρ^(5/2) [1+0.5548(√((λ*-u)/ρ)/s*)^(-0.946)] * exp[ (3/2)ρ² - (αλ*²)/(2(α-1)) ] * erfc( 3ηρ - (η√α)/(α-1) λ* ) dλ* dρ

(式中：s* = s/α^(1/4) * √(m₄/m₂) 为无量纲的表面张力参数；u = d/√m₀ 为无量纲的两表面间平均距离。)