粗糙球形表面的分形接觸力學模型

原園,張利華,徐穎強

(1.西安理工大學機械與精密儀器工程學院,710048,西安;2.西北工業大學機電學院,710072,西安)

接觸現象廣泛存在于工程領域中。掌握物體之間的接觸力學性能,對研究潤滑、摩擦、磨損及熱傳導等具體的工程實際問題具有十分重要的作用。經典接觸力學中普遍認為物體的接觸表面是光滑連續的,當兩物體接觸時,其間的實際接觸面積與名義接觸面積是相等的。然而,實驗觀測表明,物體的接觸表面是由眾多幾何尺寸不同的微凸體構成,即接觸表面是粗糙表面。由于兩個物體的接觸實際是兩接觸面上微凸體的相互作用,使得實際接觸面積遠遠小于名義接觸面積,真實接觸面積承受過大的接觸載荷,導致接觸面的壓潰、磨損,引起接觸零部件的失效。因此,研究粗糙表面在接觸過程中真實接觸面積與接觸載荷的關系,對提高接觸零部件的承載能力和使用壽命具有十分重要的理論指導作用。

目前,粗糙表面接觸問題的研究對象以平面接觸為主,研究方法主要從實驗測試[1-2]和理論研究兩方面進行。實驗測試可以直觀地得到粗糙表面接觸過程中的各項力學性能規律,驗證接觸零部件是否滿足設計要求,但實驗測試結果對儀器的測量精度和實驗環境較為敏感,測量結果只在一定范圍(材料屬性、表面形貌、載荷歷程等)內有效,不具有普遍性。理論研究所建立的力學模型和計算方法具有一定的普遍適用性,因此得到長足發展。在理論研究中,粗糙表面接觸力學模型主要分為統計學模型、分形模型、多尺度模型[3-4]、確定性模型[5]及有限元模型[6],這些模型各具優勢。多尺度、確定性和有限元模型考慮到微凸體相互疊加的因素,更準確地反映了粗糙表面的真實面貌。統計學模型和分形模型將微凸體簡化為半球體,獲得了粗糙表面接觸過程中真實接觸面積和接觸載荷的關系。由于這兩種模型求解簡便,易于擴展到不同的工程領域,應用較為廣泛。

Greenwood和Williamson于1966年首次建立了粗糙表面接觸的統計學模型,即GW模型[7]。在該模型中,所有微凸體的曲率半徑被認為是相同的,其高度服從高斯分布,并且引入塑性指數作為評判微凸體發生彈性或塑性變形的依據,進而推導出粗糙表面真實接觸面積和接觸載荷的關系。隨后,GW模型被不斷擴展來滿足不同的工況要求,例如材料的彈塑性變形[8]、材料的黏彈性變形[9]、彈塑性變形的卸載[10]、微凸體的相互作用[11]等,這些研究成果使得GW模型日趨完善。Mandelbrot于1982年觀測到物體表面形貌輪廓呈現出連續、不可微和自相似的特征,滿足分形的幾何特性[12]。隨后,Majumdar等采用Weierstrass-Mandelbrot函數模擬出二維粗糙表面輪廓[13],Majumdar和Bhushan根據二維分形表面輪廓的特征,認為微凸體的分布形式服從海洋島嶼的分布規律,建立了二維粗糙表面接觸力學模型,即MB模型,獲得了真實接觸面積與接觸載荷的解析表達式[14]。由于MB模型易于直接得到真實接觸面積與接觸載荷的解析解,很快被應用到工程中的眾多具體問題,例如:三維粗糙表面接觸[15]、焊接材料界面的接觸率[16]、接觸面的黏著磨損[17]、結合部的法向剛度[18]等。然而,在MB模型中,粗糙表面的微凸體在接觸過程中首先發生完全塑性變形,隨著接觸下壓量的增大,微凸體的變形逐漸轉變為彈性,該現象與經典的赫茲接觸模型矛盾。對此,Morag等修正了分形模型中單個微凸體的彈塑性接觸模型,使微凸體的變形過程與赫茲接觸模型保持一致[19]。Yuan等人修正了傳統的MB模型[20]。在修正模型中,微凸體的變形過程符合赫茲接觸理論,微凸體在接觸過程中承擔的接觸面積和接觸載荷與其曲率半徑成正比。該模型也被應用到三維接觸[21]、線接觸[22]及接觸卸載[23]等具體問題。

基于上述研究成果,根據分形理論,本文建立了球形分形粗糙表面的接觸力學模型。球形分形粗糙表面是一個曲面,將粗糙表面的接觸區域分成N份,每份接觸區域都可以等效成一個平面接觸問題。推導每份接觸區域上各個頻率指數的微凸體的截斷面積密度分布函數,獲得真實接觸面積與總接觸載荷的解析表達式,即可得到接觸半寬上的接觸壓力分布規律。本文可為機械工程領域的點接觸問題的理論研究、零部件設計提供一定的理論依據。

1 球形分形粗糙表面的表征

在球坐標系(r,φ,φ)中,各向同性的球形分形粗糙表面的方程[24]為

W(r,φ,φ)=r+r(G/r)D-2(lnγ/MS)1/2·

(1)

式中:W(r,φ,φ)表示粗糙表面上任意一點到圓心的距離,r為球形分形粗糙表面的基圓半徑,φ為從z軸正方向看r在xoy平面上的投影與x軸正方向的夾角,φ為r與z軸正方向的夾角;ψm,s,n為隨機相位,m和s分別為φ方向和φ方向上重疊隆起部的數量,M和S分別為其最大值,n為微凸體的頻率指數;γ為隨機輪廓的空間頻率,通常情況下取γ=1.5;D為分形維數,對于三維表面輪廓,D的取值范圍為2

nm,su=r[sin(φα)sin(φ)cos(φ-φα)+

cos(φα)cos(φ)]

(2)

式中:φα=πs/S;φα=2πm/M。

若不考慮微凸體之間的相互重疊,則粗糙球形表面上對應φ方向與φ方向上的重疊隆起部的個數為1。令M=S=1,得到φα=π,φα=2π,帶入式(2)可得nm,su=rcosφ。基于上述簡化過程,式(1)可變換為

W(r,φ,φ)=r+r(G/r)D-2(lnγ)1/2·

(3)

式中:ψnψ1,1,n;粗糙球形表面上微凸體頻率指數范圍從nmin變化到nmax,二者之間的關系為nmax=nmin+int[log(L/La)/logγ],int[]表示取括號內的值的整數部分,L為粗糙表面的取樣長度,對于粗糙球形表面,L=r[24],La為截止長度。采用式(3)生成D=2.3、G=2×10-5m時的粗糙球形表面,并在直角坐標系xyz下進行表示,結果如圖1所示。

圖1 三維粗糙球形表面輪廓

2 單個微凸體的接觸力學特性

粗糙表面上的微凸體的力學性能直接影響粗糙表面的力學性能,因此必須先確定單個微凸體在接觸過程中接觸載荷和接觸面積的關系。假設粗糙球形表面上的微凸體之間不發生相互作用,微凸體的幾何外形為半球形。在第i份接觸區域中(接觸區域的劃分在第3節中詳細描述),根據式(3),直角坐標下頻率指數為n的微凸體未變形前的輪廓方程為

(4)

式中:Zi,n表示單個微凸體輪廓曲線上任意一點的高度;ψi,n為(0,2π)范圍內隨機相位。

圖2表示了該微凸體與一剛性平面接觸的情況,虛線部分為微凸體的原始輪廓,實線部分為微凸體的變形輪廓。

圖2 單個微凸體接觸模型

圖2中:li,n為微凸體的基底長度,與取樣長度L的關系為li,n=L/γn[24];ωi,n為剛性平面的下壓量;δi,n為微凸體的高度,表達式為

δi,n=GD-2(lnγ)1/2(L/γn)3-D

(5)

Ri,n為對應的頂端曲率半徑,表達式為

(6)

(7)

(8)

當n

(9)

(10)

(11)

從式(10)(11)可以看出,ai,nec≤ai,nep≤ai,nepc,ai,nepc為臨界彈塑性接觸面積,ai,nepc=153.6ai,nec。當n

(12)

(13)

(14)

3 真實接觸面積和總接觸載荷

3.1 接觸模型的建立

一個粗糙球形表面與一個粗糙平面接觸可以轉化為一個粗糙球形表面與一個剛性平面接觸,接觸過程中基圓不發生變形,如圖3所示。

圖3 粗糙球形表面與剛性平面接觸示意圖

圖3中:d為剛性平面與基圓底部之間的距離;ω為剛性平面的下壓量;δnmin表示頻率指數為nmin的微凸體高度,δnmin與d之間的關系為ω+d=δnmin。粗糙球形表面與剛性平面之間的接觸區域是一個半徑為Rc的圓。在接觸區域與非接觸區域的邊界上,假設剛性平面剛好與頻率指數為nmin的一個微凸體頂端相切,則Rc可表示為

(15)

與接觸半寬Rc對應的圓心角為θ。若將圓心角θ平均分成N份,則接觸區域也被分為N份。s1表示接觸區域的第一份,s2表示接觸區域的第二份,si表示接觸區域的第i份,其中接觸區域的第一份為圓形,其余各份為環形。至此,粗糙球形表面接觸問題轉化成N個平面接觸問題。當剛性平面的下壓量為ω時,每份對應的下壓量也不同,第i份的下壓量可表示為

ωi=δnmin-(δnmin-ω+diz)/cosθi

(16)

(17)

其中i∈[1,N]

3.2 微凸體的截斷面積密度分布函數

整個接觸區域包含頻率指數從nmin到nmax的微凸體,將整個接觸區域分成N份,每一小份接觸區域中同樣包含頻率指數從nmin到nmax的微凸體。若用a′表示微凸體的接觸截斷面積,則整個接觸區域上微凸體截斷面積密度分布函數n′(a′)的表達式[15]為

(18)

(19)

(20)

式中ξi(N)為第i份的微凸體截斷面積密度分布函數修正系數,與第1份的修正系數ξ1(N)的關系為

ξi(N)=Δi(N)ξ1(N)

(21)

其中Δi(N)=si(N)/s1(N)=2i-1,s1(N)表示第一份接觸區域的名義面積,si(N)表示第i份接觸區域的名義面積。

(22)

(23)

(24)

(25)

為簡化計算,令Λi,min=Λi,min+1=…=Λi,max-1=Λi,max=Λi,根據文獻[20],可得

(26)

(27)

3.3 真實接觸面積和總接觸載荷

球形分形粗糙表面與剛性平面接觸時,整個接觸區域上的所有微凸體可能發生彈性變形、彈塑性變形或完全塑性變形。整個接觸區域的真實接觸面積Ar可以表示為

Ar=Are+Arep+Arp

(28)

式中Are、Arep和Arp分別是由彈性變形、彈塑性變形和完全塑性變形引起的真實接觸面積,分別表示為

(29)

(30)

(31)

與真實接觸面積對應的總接觸載荷可表示為

Fr=Fre+Frep+Frp

(32)

式中Fre、Frep和Frp分別是由彈性變形、彈塑性變形和完全塑性變形引起的總接觸載荷,這三個部分可以表示為

(33)

(34)

(35)

第i份接觸區域上的接觸壓力pi為

pi=Fi,r/Ai,r, 1≤i≤N

(36)

式中Ai,r和Fi,r分別是第i份接觸區域上的真實接觸面積和總接觸載荷,可以表示為

(37)

(38)

4 計算結果與討論

根據推導出的解析結果,研究粗糙球形表面接觸過程中真實接觸面積和總接觸載荷之間的關系。粗糙球形表面的各項計算參數[15]如表1所示。

表1 計算參數

接觸區域被分成N份,每份都包含不同頻率指數的微凸體,分別獨立給出每個頻率指數微凸體的最大截斷面積會導致結果數據冗余且無規律可循。為了使計算結果符合實際情況,規定在第i份接觸區域內,各個頻率指數的微凸體中具有最大截斷面積的微凸體,其接觸下壓量與自身高度的比值都相同。根據表1的數據,可以得到微凸體的臨界彈性頻率指數nec=33,臨界彈塑性頻率指數nepc=46。

(a)ω=0.5δnmin

(b)ω=0.4δnmin和ω=0.8δnmin圖4 不同頻率指數范圍微凸體產生的真實接觸面積與總真實接觸面積的比值變化

圖5表示當微凸體頻率指數范圍為17～65、ω=0.95δnmin時,粗糙球形表面接觸過程中總真實接觸面積Ar與總接觸載荷Fr的關系,可以看出:粗糙球形表面在整個接觸過程中近似表現為彈性性質,總接觸載荷與真實接觸面積的3/2次方成正比,即Fr～Ar3/2。這是因為臨界彈性頻率指數nec=33,所以頻率指數處于17～33的微凸體只能發生彈性變形。又由于在分形粗糙表面接觸過程中,前6個頻率指數的微凸體對整個接觸過程起主導作用,即使頻率指數處于34～65的微凸體發生彈塑性或完全塑性變形,其比例在總接觸載荷和總真實接觸面積中小于0.01。

圖5 真實接觸面積與總接觸載荷的關系

圖6與圖7表示微凸體頻率指數范圍為38～86時,粗糙球形表面接觸過程中的變形狀況。頻率指數處于38～46的微凸體在接觸過程中可以發生彈性變形和彈塑性變形,頻率指數處于47～86的微凸體可以發生彈性變形、彈塑性變形及完全塑性變形。

圖6 真實接觸面積與總接觸載荷的關系

圖7為接觸區域的變形狀況。圖7a是當ω=0.35δnmin時的接觸區域變形情況,可以看出:接觸區域被分成N1份(N1=100);第1份接觸區域發生彈塑性變形,其余各份接觸區域都呈現近似彈性變形。因為彈塑性變形與總變形的比值十分小,此時接觸變形狀態仍呈現為彈性性質。圖7b和圖7c分別是當ω=0.52δnmin和ω=0.68δnmin時的接觸區域變形情況,可以看出:隨著接觸下壓量的增加(接觸區域增大,接觸區域所分的份數也逐漸增加,以確保不同下壓量下第i份接觸區域的面積近似保持一致),彈塑性變形區域逐漸增大,彈塑性變形占總變形的比例也逐漸增大。當這個比例大于60%時,粗糙球形表面逐漸顯現彈塑性變形性質。

(a)ω=0.35δnmin

(b)ω=0.52δnmin

(c)ω=0.68δnmin圖7 接觸區域的變形狀況

圖8 真實接觸面積與總接觸載荷的關系

圖9表示接觸區域上接觸壓力的分布規律。圖9a表示當微凸體頻率指數范圍為17～65時的情況,可以分析得出:粗糙球形表面的接觸半寬與接觸壓力成正比;粗糙球形表面在接觸過程中呈現彈性性質,壓力分布曲線近似拋物線,與赫茲接觸模型的壓力分布較為相似。圖9b表示當微凸體頻率指數范圍為38～86時的情況,可以分析得出:當ω=0.3δnmin時,粗糙球形表面呈現彈性變形性質,接觸壓力分布近似為拋物線;當ω=0.6δnmin和ω=0.9δnmin時,粗糙球形表面呈現彈塑性變形性質,接觸壓力分布曲線的頂端變得較為平緩。圖9c表示當微凸體頻率指數范圍為53～101時的情況,可以分析得出:當ω=0.3δnmin時,粗糙球形表面呈現彈塑性變形性質,接觸壓力分布曲線的頂端較為平緩;當ω=0.6δnmin和ω=0.9δnmin時,粗糙球形表面呈現彎曲塑性變形性質,接觸壓力分布曲線的頂端趨于直線,其輪廓近似為一個矩形,這表明接觸壓力在接觸半寬上近似均勻分布。

(a)頻率指數范圍為17～65

(b)頻率指數范圍為38～86

(c)頻率指數范圍為53～101圖9 接觸區域上接觸壓力的分布規律

(a)接觸半寬上壓力分布的對比

(b)接觸壓力峰值的對比圖10 本文模型與赫茲接觸模型的接觸壓力對比

圖11 本文模型與赫茲接觸模型的壓力峰值相等時各分形參數之間的關系

為了獲得粗糙球形表面接觸模型與赫茲接觸模型之間的內在聯系,只考慮粗糙球形表面接觸的彈性變形階段。為尋求本文模型與赫茲接觸模型在接觸下壓量相等時,接觸壓力峰值也同樣相等的存在條件,構建兩模型中各分形參數的關系,如圖11所示,圖中Rmin表示最小頻率指數的微凸體頂端曲率半徑,可以看出,G在1.36×10-9～1.36×10-12m之間變化,對D與Rmin/r的關系影響很小。忽略G的影響,對圖11中的曲線進行擬合可得

Rmin/r=f(D)=-0.063 9D2+0.197D-0.028 7

(39)

對于給定的分形維數D,最小頻率指數的微凸體頂端曲率半徑與基圓半徑比值Rmin/r滿足式(39)。對于相同的接觸下壓量,粗糙球形表面接觸模型的接觸壓力峰值始終等于赫茲接觸模型的接觸壓力峰值。如果Rmin/r>f(D),則粗糙球形表面接觸模型的壓力峰值小于赫茲接觸模型的壓力峰值;否則,粗糙球形表面接觸模型的壓力峰值大于赫茲接觸模型的壓力峰值。從工程應用考慮,為了確保接觸零件的安全使用,在額定下壓量的情況下,應盡量增大表面微凸體頂端曲率半徑,降低表面粗糙度。

5 結 論

(1)微凸體的頻率指數范圍直接影響粗糙球形表面的接觸力學性質,前6個頻率指數的微凸體在整個接觸過程中起主導作用。最小頻率指數nmin與臨界彈性頻率指數nec滿足nmin+5≤nec,粗糙球形表面在整個接觸過程表現出彈性變形行為。最小頻率指數nmin處于臨界彈性頻率指數nec與臨界彈塑性頻率指數nepc之間,粗糙表面在初始接觸時表現為彈性變形行為,隨著接觸載荷的增大,逐漸呈現為彈塑性變形行為。當nmin>nepc時,粗糙球形表面在整個接觸過程中都呈現非彈性變形行為,這是因為彈性變量形遠遠小于總變形量。

(2)粗糙球形表面的接觸半寬主要取決于基圓半徑,微凸體的幾何尺寸對其影響很小。當粗糙球形表面處于彈性或彈塑性變形階段時,接觸壓力分布曲線近似為拋物線,接觸壓力在接觸區域中心達到最大,向接觸區域邊緣方向遞減。在完全塑性變形階段,接觸壓力分布曲線的輪廓近似為矩形,接觸壓力在整個接觸區域上均勻分布。