通過對式(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₀ 為無量綱的兩表麵間平均距離。)