粗糙表面接觸面積和承載規(guī)律的研究

王世軍, 何花蘭, 郭璞, 趙金娟

(1.西安理工大學 機械與精密儀器工程學院，陜西 西安 710048；2.西安理工大學 印刷包裝工程學院，陜西 西安 710048)

零件之間相互接觸的表面稱為機械結合面。宏觀上光滑的機械結合面，在微觀上仍然是粗糙的，其宏觀的載荷-變形關系與基體部分明顯不同，例如機床60%～80%的剛度、90%以上的阻尼都來自結合面。要建立準確的機械整機結構模型，就必須了解結合面的靜態(tài)和動態(tài)特性。結合面靜態(tài)和動態(tài)特性的研究一直是機械整機性能分析的難點之一，目前仍沒有一套比較成熟的理論與方法，大多數情況下仍然是靠經驗估計。建立機械結合面模型，準確辨識結合面參數便成為機械結構性能分析中的一項關鍵技術。機械結合面研究要達到文獻[1]提出的“結合面特性參數應反映結合面的基本特性，具有通用性，類似零件的特性常數如彈性模量、損耗因子，能方便地供設計解析使用”這一目標，道路依然還很漫長[2]。

粗糙表面的接觸面積、接觸壓力分布等接觸特性會直接影響結合面的接觸剛度、接觸熱阻、配合精度、密封等結合面的性能。目前在微觀粗糙表面的接觸研究中提出了許多理論解析模型，這些接觸研究模型主要有統(tǒng)計模型、分形模型和有限元模型三類。統(tǒng)計模型的基礎是Greenwood 和Willimson建立的GW模型[3]。GW模型將粗糙峰簡化為彈性半球體，并且假定粗糙峰的高度符合正態(tài)分布，利用赫茲接觸理論分析接觸表面的載荷-變形關系。GW模型被許多研究者改進，其中比較重要的是Chang等人在GW模型的基礎上計入了塑性變形的影響后提出的CEB模型[4]。Majumdar和Bhushan認為粗糙表面具有分形特征，提出了以分形幾何為基礎的MB接觸模型[5]。該模型采用Weierstrass-Mandelbrot分形函數描述粗糙表面的幾何形貌特征，采用赫茲方法分析接觸表面的載荷-變形關系。統(tǒng)計模型和分形模型均依賴于一定的假設條件,根據這些假設條件可以對表面輪廓特征進行簡化，如微凸體都具有相同的峰頂曲率半徑，其變形相互獨立，微凸體之間的接觸是峰對峰的正接觸。這些假設與簡化限制了理論解析模型的應用范圍，同時影響了粗糙表面接觸分析的準確性[2-6]。

文獻[7]采用有限元法研究了具有幾個圓柱形粗糙峰的剛性表面與彈性半無限體相互接觸的情況。文獻[8]采用具有半徑和高度相同的一定數目圓形粗糙峰的剛性表面與彈塑性半無限體相互接觸的模型來模擬多粗糙峰接觸，并用有限元法對該模型進行了彈塑性分析。文獻[9-10]利用由分形幾何模擬的三維粗糙表面將粗糙表面的接觸等效成粗糙表面與剛性平面的接觸并進行了相關的有限元分析。文獻[11]則認為對于具有許多任意形狀微凸體的三維粗糙表面的接觸問題，采用有限元法需要大量的單元，計算量過大。雖然隨著計算機技術的發(fā)展,有限元計算的能力已經大幅度提高，采用有限元法求解具有強烈非線性的三維粗糙表面的接觸問題仍然存在名義接觸面積過小、計算規(guī)模過大的問題，例如文獻[6]中三維有限元模型的名義接觸面積最大為0.25×0.25 mm2，單元數高達52萬，文獻[12]中的名義接觸面積僅為2.5×2.5 μm2，并且兩個接觸面只有一個是粗糙表面，另一個表面為剛性平面，單元數為23 534。一般認為，過小的名義接觸面積難以反映機械加工表面的紋理特征。文獻[13]直接采用表面形貌測量儀獲得的粗糙表面二維數據建立具有兩個粗糙接觸表面的有限元彈塑性模型，模型的名義接觸長度為5.6 mm，32 236個單元。與三維模型相比，二維模型接觸面的名義長度可以較長，單元數量較少，而且結合面法向的載荷-變形計算結果與實驗結果也能夠很好地吻合。

接觸問題研究中，一般認為三維有限元接觸模型比二維模型更能反映真實的結合面接觸狀況，但是因為存在前面所述的問題，目前的研究中主要還是采用二維模型。因為二維模型的局限性，對實際接觸性質的反映可能是不全面的，甚至是不正確的。一個二維模型在多大程度上反映真實的接觸性質，需要從接觸性質的各個方面具體考察。已有的研究表明[13]，筆者提出的這種二維模型的法向變形-載荷關系與實驗結果是一致的。在本文中，筆者進一步考察了微觀接觸面積與載荷的關系，并與已有的電接觸法[14]獲得的實驗結果進行了比較，更深入、廣泛地考察這一模型的可行性。

1 粗糙表面的輪廓

圖1是通過TR300型表面輪廓測量儀獲得的磨削試樣的表面輪廓曲線。表面輪廓的測量是沿試樣表面垂直于磨削紋理方向掃描，測頭直徑2 μm,數據采樣間距1 μm,采樣長度5.6 mm。測量獲得的數據經過最小二乘處理后,將數據點用直線連接起來形成圖1所示的表面輪廓。

2 有限元模型

圖2是利用圖1的表面輪廓數據建立的二維有限元平面應變模型。圖2(a)是兩個接觸試樣的整體網格，圖2(b)是接觸表面附近的局部網格。為了盡可能減少單元數量，減少計算量，本文沿垂直方向分層逐步加密網格，以便網格疏密過渡合適。為了保證接觸計算的收斂性和分析精度，在接觸表面附近采用均勻的四邊形網格，同時，為了保證網格水平方向的均勻性，將試樣沿水平方向也進行了均勻分塊，在接觸表面附近獲得了比較高的網格質量和合適的網格密度。整個模型的單元數為148 661，結點數為148 650。

圖1 磨削表面的輪廓

圖2 接觸試樣的有限元模型

圖2(a)中下試樣底部的法向位移約束為0，同時左下角的水平位移也為0，兩個試樣在接觸表面間通過拉格朗日乘子法建立接觸關系[16]。模型的法向載荷施加在上試樣的頂面，最大載荷為1 MPa，按120個載荷步等步長加載。

試樣材料為45鋼，由單軸拉伸實驗獲得的真應力-應變關系如圖3所示，本文將其簡化為圖中的多線性彈塑性應變硬化材料模型，采用Von Mises判據作為材料的塑性判據。材料彈性模量E=207 GPa，材料屈服極限σs= 320 MPa，泊松比μ=0.3。

3 應力分析的結果

圖4顯示了法向載荷為1 MPa時模型中的Mises應力分布，圖4(a)是整體結構中的應力分布，圖4(b)是圖4(a)中左側高應力區(qū)的局部應力分布。由圖4(a)可以看到，接觸層存在兩個明顯的高應力接觸區(qū)，試樣中的應力分布不均勻，上下試樣中的應力分布相似。圖5顯示了上試樣在不同載荷時接觸區(qū)各結點的Mises應力。圖中的接觸結點序號是將接觸區(qū)內的結點按應力值排序后的序號。由圖中可以看出，隨著載荷的增大，進入接觸的結點數在增多，應力水平升高，同時各結點的應力水平也趨于一致，而且即便是在載荷比較低的情況下，也總是存在接觸點的應力超過屈服極限的情況。從圖5可以看到，隨著載荷增大，彈性接觸點數量減少，塑性接觸點增加，即使在載荷為0.25 MPa時，也有結點的最高應力超過了材料的屈服極限，最高應力大約是試樣名義應力的1 400倍。下試樣的分析結果也有類似的情況。

圖3 材料的應力-應變曲線

圖4 1 MPa時的應力分布

圖5 不同載荷時的結點應力

圖6是圖4(b)中上試樣最高應力點由接觸表面向上至試樣內部路徑上的Mises應力分布。由圖6可以看出，高應力區(qū)只出現在很小范圍內，遠離接觸點的試樣內部的應力遠低于接觸點的應力，在距離接觸點50 μm處，Mises應力已經降低到接觸點應力的10%以下，在100 μm處則降低到5%以下。結合圖4的結果，說明試樣內部的應力并不均勻，接觸層的應力狀態(tài)實際由這些真實的接觸點的應力狀態(tài)決定。此外，在法向載荷低于1 MPa時，接觸點附近的Mises應力分布趨勢與1 MPa時相同，只是應力水平較低，高應力區(qū)的范圍更小一些。

圖6 接觸點附近的應力分布

4 承載力的分析結果

圖7是圖4中上試樣接觸表面各接觸結點在不同載荷作用下的接觸力沿垂直方向分量占法向總載荷的比重，反映了接觸表面上的接觸結點在不同載荷下的法向承載能力，即各結點對法向載荷的貢獻大小。圖7中橫軸是接觸中的結點的排序號，按結點法向力比重從大到小排列。縱軸是接觸結點垂直分力占試樣法向載荷的比重。由圖中可以看到，在法向載荷較低的情況下，接觸點數很少，各結點承載力差別很大。隨著載荷增大，進入接觸的結點數增多，承載力趨于一致，這與圖5的Mises應力變化趨勢一致，但與圖5不同的是，不同載荷條件下總是存在少數承載力很低的結點，在法向載荷升高時，這類結點所占的比重在減小。結合圖5的結果可以發(fā)現，在載荷為1 MPa時，接觸點的承載能力趨于一致的同時，接觸點都處在塑性應力狀態(tài)。

圖7 不同載荷時結點的法向承載力的變化

由于輪廓測量是等間距采樣，進入接觸的結點數反映了法向接觸面積的大小，即真實承載面積的大小。圖7中法向載荷為1 MPa時，接觸結點只有7個，占上試樣表面結點數的0.13%，真實接觸面積只占名義接觸面積很小的一部分，使得圖5中接觸結點的應力水平遠高于試樣的名義應力，真實接觸區(qū)的應力水平和應力狀態(tài)與遠離接觸點的試樣內部有很大不同，接觸層的性質與試樣內部明顯不同。

5 法向變形與法向載荷的關系

圖8給出了計算獲得的法向位移與載荷關系曲線，法向的位移取自圖4(a)中上試樣的頂面中心點。圖中的試驗曲線采用文獻[13]中的試驗方法獲得，試驗中的試樣采用環(huán)形試樣，環(huán)形接觸面的外徑為45 mm，內徑為25 mm，試樣高度為50 mm。計算獲得的位移與載荷關系大體呈線性，與圖中的試驗曲線比較，計算結果與試驗結果能較好地吻合，表明模型的宏觀計算結果是可靠的。

圖8 法向位移與載荷的關系

圖9是載荷為1 MPa時圖4(b)中接觸點附近上試樣沿垂直方向路徑上的法向變形，是路徑上各點相對于路徑起點的相對位移，反映了接觸點附近垂直方向的變形程度。由圖9可以看出，接觸點附近的變形主要集中在距離接觸點0～10 μm的范圍內，這個結果與圖6所示的路徑上的應力分布規(guī)律一致。

圖9 接觸點附近的法向變形

6 接觸面積與載荷的關系

圖10給出了有限元計算獲得的接觸面積比與法向載荷的關系：隨著載荷的升高，真實接觸面積大體按線性規(guī)律增大，載荷為1 MPa時，真實接觸區(qū)只占名義接觸面積的0.13%。作為對照，圖中也給出了文獻[14]中采用接觸電阻法獲得的真實接觸面積的測量結果。由于接觸問題的特殊性，直接觀察和測量真實接觸區(qū)的接觸狀態(tài)和接觸面積比較困難。接觸電阻法是通過測量接觸表面的接觸電阻估算真實接觸面積的一種方法，該方法認為真實接觸面積與接觸電阻成正比，通過測量試樣間的接觸電阻即可確定真實的接觸面積。文獻[14]中給出的測量結果是采用磨削過的、面積為21 cm2鋼平面的測量結果。由圖中可見，結果測量與計算結果一致，說明真實接觸面積的有限元計算結果是可靠的。

圖10 有限元法與電阻法的面積比較

接觸區(qū)數目的增多和接觸區(qū)面積的增加都可以導致接觸面上真實接觸面積增加。結合圖4、5的進一步分析表明，在研究的載荷范圍內，真實接觸面積的增大并不是由于接觸區(qū)數目的增多，而是微觀接觸區(qū)面積的增大。

7 真實接觸區(qū)的高度分布

圖11為圖4中上試樣接觸表面輪廓高度的概率密度分布圖，圖12為對應的累積概率圖。由圖11可以看出，磨削表面輪廓大體符合正態(tài)分布，但是也存在一定出入。由圖12可以看出，在輪廓的平均位置附近，輪廓高度分布比較接近理想的正態(tài)分布，但是在曲線的兩端，即輪廓高度的兩端與理想的正態(tài)分布還是有一定差距。

圖11 輪廓高度的概率密度分布

圖12 輪廓高度的累積概率圖

圖12中標出了上試樣中接觸點的高度在累積概率圖中的位置。由圖中可以看出，接觸點都處于平均高度以上的位置而不是散布在整個輪廓高度范圍內，輪廓平均高度以下的位置不產生接觸，下試樣以及5組同一批次的試樣都有同樣的特點。由于真實接觸點出現在輪廓高度的上部位置，而上部位置的輪廓高度與正態(tài)分布的相對偏差通常比中部大，這使得基于輪廓特征符合正態(tài)分布假定的GW模型的計算結果與實際情況相比具有比較大的誤差。此外，由于真實接觸點的數量很少并且不出現在輪廓高度的平均值附近以及平均值以下的位置，這些位置上輪廓數據的分布規(guī)律不會影響接觸的狀態(tài)和性質，圖11中輪廓高度分布的左側和中間的形狀即使不滿足正態(tài)分布，接觸層也可以有相同的接觸性質，如法向變形和剛度等。這一觀點也通過下面的計算得到了驗證：通過修改圖2中表面輪廓高度中平均位置以下的數據點的高度值以后建立新的有限元模型，其接觸區(qū)的應力分布與宏觀的變形和修改之前的圖2中的模型并沒有不同。當然，這一結論的正確性和適用范圍尚需更多、更廣泛的計算和試驗結果的驗證以及理論分析。

8 結 論

采用粗糙表面輪廓數據建立的二維平面應變接觸模型不僅能夠準確反映接觸層宏觀的法向位移-載荷關系，由其計算獲得的真實接觸面積與載荷的關系也能夠跟已有的接觸電阻法的測量結果吻合，說明該模型的微觀計算結果也具有一定的可靠性，可以用作結合面微觀接觸狀態(tài)的觀察和接觸性質研究的較為直觀的工具。

相對于接觸問題的三維有限元模型，本文的二維模型在規(guī)模和效率上都更有優(yōu)勢。圖2中的二維模型有44.6萬自由度，如果采用矩形的三維模型計算，接觸面尺寸為5.6×5.6 mm2，采用同樣的網格密度，自由度數將超過22.5億。目前有限元軟件的最大處理能力只有2億自由度，常用的有限元軟件如ANSYS，其商業(yè)版本目前的處理能力只有2千萬自由度。因此，現階段采用三維有限元模型求解較大尺寸的接觸問題仍然存在困難，相比之下，二維模型還有較大的自由度空間。

現有的計算結果表明，真實接觸面積很少，載荷為1 MPa時僅占名義接觸面積的0.13%，并且真實接觸點只出現在輪廓高度平均值以上的位置，該部分輪廓高度的分布與正態(tài)分布有一定差異，輪廓高度平均位置附近及其以下的高度分布不影響接觸性質，表面輪廓符合正態(tài)分布應該不是GW模型的必要條件。

當法向載荷在1 MPa以下時，粗糙表面接觸面積隨載荷的增加而線性增大，接觸區(qū)的法向承載力隨載荷增大趨于一致。在很小的法向載荷下，接觸區(qū)就出現塑性，可以認為接觸狀態(tài)下，接觸區(qū)一直存在塑性。隨著載荷的增加，塑性區(qū)占的比重逐漸增大。

參考文獻：

[1]黃玉美,付衛(wèi)平,董立新,等.結合面法向動態(tài)特性參數研究[J].機械工程學報,1993,29(3): 74-78.

Huang Yumei, Fu Weiping, Dong Lixin, et al.Research on the dynamic normal characteristic parameters of joint surface[J].Journal of Mechanical Engineering, 1993, 29(3): 74-78.

[2]田紅亮,朱大林,秦紅玲.機械結合面接觸模型的研究進展及存在的問題[J].機械工程與技術, 2013, (2): 1-10.

Tian Hongliang,Zhu Dalin,Qin Hongling.Advances and existing issues in contact model of mechanical joint surface[J].Mechanical Engineering and Technology, 2013,(2):1-10.

[3]Greenwood J A, Williams J B.Contact of nominally flat surfaces[J].P.Roy.Soc.Lond.: A Mat., 1966, 295(1442): 300-319.

[4]Chang W R,Etsion I, Bogy D B.An elastic-plastic model for the contact of rough surfaces[J].ASME J.Tribol, 1987, 109(2):257-263.

[5]Majumdar A, Bhushan B.Fractal model of elastic-plastic contact between rough surfaces[J].ASME J.Tribol, 1991, 113(1):1-11.

[6]楊國慶,熊美華,洪軍,等.3D粗糙表面的數字化表征與接觸特性分析[J].西安交通大學學報, 2012,46(11): 58-63.

Yang Guoqing, Xiong Meihua, Hong Jun,et al.Numerical charecterization and contact performances for 3D rough surfaces[J].Journal of Xi’an Jiaotong University, 2012,46(11):58-63.

[7]Komvopoulos K，Choi D H.Elastic finite element analysis of multi-asperity contact[J].ASME Journal of Tribology,1992，114:823-831.

[8]楊楠，陳大融，孔憲梅.多粗糙峰彈塑性接觸的有限元分析[J].摩擦學學報, 2000, 20(3): 202-206.

Yang Nan, Chen Darong, Kong Xianmei.Elastic-plastic finite element analysis of multi-asperity contacts[J].Tribology, 2000, 20(3): 202-206.

[9]Sahoo P, Ghosh N.Finite element contact analysis of fractal surfaces[J].J.Phys: D Appl.Phys., 2007, 40(14): 4245-4252.

[10]Pei L, Hyun S, Molinari J F, et al.Finite element modeling of elasto-plastic contact between rough surfaces[J].J.Mech.Phys.Solids,2005,53(11):2385-2409.

[11]Bhushan B.Contact mechanics of rough surfaces in tribology：Multiple asperity contact[J].Tribology Letters，1998，4:1-35.

[12]王霄,梁春,劉會霞,等.車削真實粗糙表面的彈塑性接觸有限元分析[J].潤滑與密封,2008,33(12):72-74.

Wang Xiao, Liang Chun, Liu Huixia, et al.Finite element analysis based on elastoplastic contact of real turning rough surface[J].Lubrication Engineering, 2008,33(12):72-74.

[13]王世軍,趙金娟,張慧軍,等.一種結合部法向剛度的預估方法[J].機械工程學報, 2011, 47(21):111-115.

Wang Shijun, Zhao Jinjuan, Zhang Huijun, et al.A method of estimating normal stiffness of joint[J].Journal of Mechanical Engineering, 2011, 47(21): 111-115.

[14]鮑登F P,泰伯D.固體的摩擦與潤滑[M].陳紹澧,袁漢昌,丁雪加,譯.北京:機械工業(yè)出版社, 1982:35-37.

[15]Poulios K, Klit P.Implementation and applications of a finite-element model for the contact between rough surfaces[J].Wear, 2013, 303:1-8.

[16]Bathe K J, Bouzinov P A.On the constraint function method for contact problems[J].Computers & Structures, 1997,64(6)：1069-1085.