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

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

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

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

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

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

1 粗糙表面的輪廓

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

2 有限元模型

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

圖1 磨削表面的輪廓

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

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

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

3 應(yīng)力分析的結(jié)果

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

圖3 材料的應(yīng)力-應(yīng)變曲線

圖4 1 MPa時的應(yīng)力分布

圖5 不同載荷時的結(jié)點應(yīng)力

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

圖6 接觸點附近的應(yīng)力分布

4 承載力的分析結(jié)果

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

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

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

5 法向變形與法向載荷的關(guān)系

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

圖8 法向位移與載荷的關(guān)系

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

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

6 接觸面積與載荷的關(guān)系

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

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

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

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

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

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

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

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

8 結(jié) 論

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

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

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

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

參考文獻：

[1]黃玉美,付衛(wèi)平,董立新,等.結(jié)合面法向動態(tài)特性參數(shù)研究[J].機械工程學(xué)報,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]田紅亮,朱大林,秦紅玲.機械結(jié)合面接觸模型的研究進展及存在的問題[J].機械工程與技術(shù), 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粗糙表面的數(shù)字化表征與接觸特性分析[J].西安交通大學(xué)學(xué)報, 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].摩擦學(xué)學(xué)報, 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]王世軍,趙金娟,張慧軍,等.一種結(jié)合部法向剛度的預(yù)估方法[J].機械工程學(xué)報, 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.