通過對(duì)式(16)求關(guān)于r的一階導(dǎo)數(shù)，有

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

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

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

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

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

步驟4：通過“步驟3”的求解即可獲得接觸區(qū)域的壓力分布p(t)，并判斷接觸邊緣處的壓力是否為有限的正值，同時(shí)需要保證接觸區(qū)域內(nèi)的壓力p(t)是光滑連續(xù)變化的，否則繼續(xù)執(zhí)行“步驟2”與“步驟3”，直至獲得合理的壓力分布為止。

步驟5：類似地，采用高斯-切比雪夫求積公式對(duì)式(17)進(jìn)行化簡(jiǎn)，然后將上述獲得的真實(shí)接觸半徑a以及壓力分布p(t)代入方程，即可求解出彈性半空間的壓入深度ω。

實(shí)際工程中不同粗糙表面對(duì)應(yīng)的微凸體曲率半徑各不相同，在上述的數(shù)值求解過程中，以微凸體曲率半徑R=10 mm為例來研究表面張力對(duì)接觸特性的影響，其他參數(shù)G=1 MPa, μ=0.4, β=0.1 J/m2。獲得了無量綱接觸載荷F/(4/3 E* R1/2 ω3/2)和無量綱接觸面積A/πRω與無量綱參數(shù)(Rω)1/2/s的數(shù)值關(guān)系，通過對(duì)數(shù)值解的擬合，得到了考慮表面張力時(shí)載荷F和真實(shí)接觸面積A與壓入深度ω的關(guān)系式分別如式(19)和式(20)所示，且不同的R值只會(huì)影響式(19)和式(20)中的擬合系數(shù)。

圖4為無量綱接觸載荷F/(4/3 E* R1/2 ω3/2)與無量綱參數(shù)(Rω)1/2/s的關(guān)系。由圖4可知，式(19)與數(shù)值結(jié)果吻合，說明了考慮表面張力時(shí)的接觸載荷與壓入深度關(guān)系模型的正確性。同時(shí)發(fā)現(xiàn)，隨著(Rω)1/2/s的減小，F(xiàn)/(4/3 E* R1/2 ω3/2)逐漸增大。

圖5為無量綱真實(shí)接觸面積A/πRω與無量綱參數(shù)(Rω)1/2/s的關(guān)系。類似地，數(shù)值結(jié)果與式(20)吻合，說明了考慮表面張力時(shí)的真實(shí)接觸面積與壓入深度關(guān)系模型的正確性。此外，隨著(Rω)1/2/s的減小，A/πRω也逐漸減小。

同時(shí)，從圖4和圖5可知，當(dāng)不計(jì)表面張力的影響時(shí)，即s趨于無窮小，參數(shù)(Rω)1/2/s趨于無窮大時(shí)，本模型逐漸趨近于赫茲接觸模型；相反，當(dāng)考慮表面張力的影響時(shí)，接觸載荷與真實(shí)接觸面積明顯不同于赫茲接觸模型的分析結(jié)果；當(dāng)壓入深度相同時(shí)，與赫茲求解接觸模型的分析結(jié)果相比，由于表面張力的存在，接觸載荷較大，真實(shí)接觸面積較小。

根據(jù)剛度的定義，可得單個(gè)微凸體與剛性平面接觸時(shí)的法向接觸剛度為

圖6顯示了無量綱接觸剛度與無量綱參數(shù)(Rω)1/2/s的關(guān)系曲線。由圖6可知，數(shù)值結(jié)果與式(21)吻合。當(dāng)忽略表面張力的影響時(shí)，即參數(shù)(Rω)1/2/s趨于無窮大，接觸剛度逐漸接近于赫茲接觸分析結(jié)果；相反，當(dāng)表面張力的影響不可忽略時(shí)，此時(shí)接觸剛度明顯偏離于赫茲結(jié)果。

3 結(jié)合面新模型

受外界載荷作用，微凸體與剛性平面產(chǎn)生接觸，單個(gè)微凸體的變形量ω=z-d。采用統(tǒng)計(jì)學(xué)方法并結(jié)合微凸體高度與曲率分布的概率密度函數(shù)式(10)，將單個(gè)微凸體的計(jì)算模型擴(kuò)展到整個(gè)結(jié)合面上，獲得了結(jié)合面的接觸載荷Ft、真實(shí)接觸面積At、接觸剛度Kt分別為

將式(19)、式(20)和式(21)分別代入式(22)、式(23)和式(24)，并進(jìn)行無量綱處理有

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) 無量綱化后的總真實(shí)接觸面積表達(dá)式：

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) 無量綱化后的總接觸剛度表達(dá)式：

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?) 為無量綱的表面張力參數(shù)；u = d/√m? 為無量綱的兩表面間平均距離。)