通過對式(16)求關于r的一階導數，有

式中：H1與Y1分別為1階斯特魯夫函數與第二類1階的貝塞爾函數。

由于式(18)為奇異積分，難以獲得它的解析解，故本文采用高斯-切比雪夫求積公式對單個微凸體模型進行數值求解，詳細的求解步驟如下，圖3給出了數值求解流程步驟為：

步驟1：先給定一個法向載荷F的大小。

步驟2：合理假設載荷F作用下接觸半徑a的最大取值amax與最小值amin，并采用二分法取其初始值a0=(amax+amin)/2。

步驟3：使用高斯-切比雪夫求積公式將式(18)及其約束方程式(13)化簡成積分區(qū)間為[-1,1]的高斯切比雪夫型求積方程，然后將a0代入方程中，并使用MATLAB進行數值迭代求解。

步驟4：通過“步驟3”的求解即可獲得接觸區(qū)域的壓力分布p(t)，并判斷接觸邊緣處的壓力是否為有限的正值，同時需要保證接觸區(qū)域內的壓力p(t)是光滑連續(xù)變化的，否則繼續(xù)執(zhí)行“步驟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)與數值結果吻合，說明了考慮表面張力時的接觸載荷與壓入深度關系模型的正確性。同時發(fā)現，隨著(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。采用統(tǒng)計學方法并結合微凸體高度與曲率分布的概率密度函數式(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 - (αλ*2)/(2(α-1)) ] * erfc( 3ηρ - (η√α)/(α-1) λ* ) dλ* dρ

(26) 無量綱化后的總真實接觸面積表達式：

A_t* = √(3α2/(32π(α-1))) ∫∫ (λ*-u) ρ2 [ (1+0.2181(√((λ*-u)/ρ)/s*)^(-0.8505)) / (1+1.355(√((λ*-u)/ρ)/s*)^(-0.7641)) ] * exp[ (3/2)ρ2 - (αλ*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 - (αλ*2)/(2(α-1)) ] * erfc( 3ηρ - (η√α)/(α-1) λ* ) dλ* dρ

(式中：s* = s/α^(1/4) * √(m?/m?) 為無量綱的表面張力參數；u = d/√m? 為無量綱的兩表面間平均距離。)