微動接觸中分形粗糙表面的接觸應力研究

阮曉光，麻詩韻，李 玲，蔡安江

（西安建筑科技大學機電工程學院，陜西 西安 710055）

1 引言

微動廣泛地發(fā)生于緊固配合的零部件之間，會引起部件的磨損和疲勞失效。影響微動的因素有很多，主要包括載荷形式、材料屬性、潤滑條件、環(huán)境溫度以及接觸面表面粗糙度[1]。這些因素主要是通過影響接觸面的接觸壓力和剪切摩擦力分布進而引發(fā)結合面的微動損傷。

許多學者已經對各種微動接觸下結合面的微動行為進行了研究，得到了很多有價值的結論[2-3]。但大多數(shù)文獻所采用的模型都是光滑接觸表面，忽略了表面粗糙度對微動行為的影響。部分學者實驗研究了表面粗糙度對摩擦行為的影響[4-5]。文獻[4]通過平面接觸的扭轉微動磨損裝置在不同的表面粗糙度和紋理方向下檢測摩擦力矩，結果表明摩擦力矩和累積耗散能量隨著表面粗糙度的增加而增加?？紤]結合面的表面粗糙度時，結合面接觸區(qū)域離散化，采用實驗的方法很難對接觸區(qū)域內每一點的接觸應力進行測量。隨著計算機運算速度的提高，有限元仿真技術取得了長足的發(fā)展。一部分學者采用有限元方法分析表面粗糙度對摩擦行為的影響[6-7]。文獻[7]提出了一種多尺度方法來研究微動磨損中的粗糙度效應。結果表明在完全滑移條件下粗糙度對微動磨損的影響較小。

建立包含粗糙表面的接觸模型的關鍵在于如何表征粗糙表面。通常用于描述表面粗糙度的方法大都依賴于粗糙表面的統(tǒng)計參數(shù)，而這些統(tǒng)計參數(shù)受樣本大小的影響，不能夠準確地反映粗糙表面的全部特征。而利用分形理論模擬的分形表面具有連續(xù)性，不可分性，尺度不變性和自親和性的特征[8-9]，克服了表面粗糙度表征對樣本大小的依賴，對粗糙表面的描述更穩(wěn)定、更精確。因此，一些學者利用分形表面模擬粗糙表面來研究結合面的接觸特性[10-13]。文獻[11]基于分形理論建立了剛性光滑平面與粗糙平面三維接觸模型，獲得外加載荷與結合面變形之間的函數(shù)表達式。文獻[12]以分形理論為基礎，推導了考慮表面微凸體相互作用影響的結合面接觸剛度分形模型。

這里的目的在于研究微動接觸中接觸面的表面粗糙度對接觸壓力和剪切摩擦力分布的影響，同時討論了材料的彈性、彈-塑性和載荷幅值對剪切摩擦力的影響。為了達到這個目的，編寫Python 腳本將Matlab 中利用Weierstrass-Mandelbrot 函數(shù)（W-M函數(shù)）構造的分形曲線的輪廓坐標數(shù)據(jù)導入ABAQUS，并使用樣條曲線擬合輪廓坐標，從而構建包含粗糙表面的二維柱面/平面接觸模型，進行有限元仿真。

2 二維粗糙表面的表征

分形理論模擬的粗糙表面的二維（2-D）表面輪廓高度由WM 函數(shù)給出[14]。

式中：z（x）—分形表面輪廓高度；x—輪廓的位置坐標；G—特征尺度系數(shù)；D—分形維數(shù)；γn—輪廓的空間頻率，γ—常數(shù)，通常取γ=1.5 較適于表面輪廓特征；輪廓結構的最低頻率與粗糙度樣本長度關系為γn1≈1/L1；n1—W-M 函數(shù)的初始項；L1—粗糙度樣本取樣長度。

根據(jù)式（1）的函數(shù)，當G=0.01nm，L=1mm，γ=1.5，n1=5，n2=70時，取不同的分形維數(shù)D，在Matlab 中分別繪制二維表面輪廓，如圖1 所示。圖1 中的二維粗糙表面分形曲線，能夠反映粗糙表面的微觀形貌。圖1 表明，隨著分形維數(shù)D 的增大，曲線波動趨于一致，表現(xiàn)出自相似性，處處連續(xù)但不可導；表面輪廓的峰值和谷值變小，表面愈加光滑，粗糙度減小。根據(jù)式（1）的函數(shù)，當D=1.35，L=0.3mm，γ=1.5，n1=5，n2=70 時，分別取不同的特征尺度系數(shù)G，仿真生成二維表面形貌，如圖2 所示。L 取0.3mm 是為了清晰地觀察圖中分形曲線輪廓特征。圖2 表明，隨著特征尺度系數(shù)G的不斷增大，圖像波動形式一致，但表面峰值和谷值逐漸增大，表面輪廓愈加粗糙，粗糙度增大。

圖1 不同分形維數(shù)下的分形曲線Fig.1 Fractal Curves Under Different Fractal Dimensions

圖2 不同特征尺度系數(shù)下的分形曲線Fig.2 Fractal Curves Under Different Characteristic Scale Factors

3 仿真分析

3.1 有限元模型

基于微動實驗模型建立二維柱面/平面接觸有限元模型，如圖3 所示。在圓柱試件的頂面施加均布的法向載荷P，并通過多點約束選項（MPC）在頂面中心點施加切向位移載荷S。下試件兩側約束水平方向位移，底部約束豎直方向位移。

圖3 有限元模型Fig.3 Finite Element Model

對有限元模型的接觸面進行接觸屬性設置，法向定義為硬接觸；切向定義為具有各向同性摩擦的庫侖摩擦，摩擦系數(shù)為0.6。為了獲得更加準確的接觸參數(shù)，采用有限滑移算法，并將切向約束定義為拉格朗日乘數(shù)法。在有限元模型中，采用主-從面的設置定義兩個接觸面，主面為圓柱試件的下表面，從面為下試件的上表面。這種設置便于提取接觸區(qū)域的接觸變量。圓柱試件的半徑為5mm，下試件的長寬分別為10mm 和5mm。為了更加深入地研究表面粗糙度對結合面接觸應力的影響，分別建立了光滑表面接觸模型和考慮粗糙表面的接觸模型。當建立考慮粗糙表面的接觸模型時，創(chuàng)建Python 腳本將Matlab 中利用W-M 函數(shù)構造的分形曲線的輪廓坐標數(shù)據(jù)導入ABAQUS，并使用樣條曲線擬合輪廓坐標，從而在下試件頂面中心區(qū)域引入長度為0.6mm 的分形曲線來模擬粗糙表面。有限元模型的網格類型為四節(jié)點平面應變單元（CPE4），并對模型接觸區(qū)域附近的網格進行細化處理。光滑接觸模型網格細化后的單元尺寸為（5×5）μm，粗糙表面接觸模型細化后的網格單元為（1×1）μm。一次仿真分析具體的加載歷程，如圖4 所示。其包括兩個分析步，法向載荷P 在分析步1 施加，并在分析步2 中保持恒定施加，切向位移載荷S 在分析步2 施加，并在分析步結尾達到最大值。

圖4 加載歷程Fig.4 Loading History

3.2 材料屬性

為了研究不同接觸材料和材料屬性對粗糙表面接觸模型接觸應力分布的影響，考慮了steel 和Al 兩種材料，并分別模擬了完全彈性和彈-塑性兩種情況，上下試件采用同種材料。兩種材料的力學屬性，如表1 所示。兩種材料的真實應力和塑性應變的關系曲線，如圖5 所示。

表1 材料屬性Tab.1 Material Properties

圖5 真實應力和塑性應變的關系曲線Fig.5 The Relationship Curve of Between True Stress and Plastic Strain

4 結果與討論

4.1 不同條件下的增量步數(shù)

當模型的材料屬性和載荷條件不同時，求解模型所需的增量步數(shù)也會發(fā)生改變。表2 給出了各種不同條件下，求解一次仿真所需的增量步數(shù)。表中1 代表求解分析步1 所需的增量步數(shù)，2代表分析步2 所需增量步數(shù)。從表2 可以看出，考慮材料彈-塑性（E-P）時，分析步1 所需增量步數(shù)明顯增大。當材料屬性為完全彈性（E）時，分析步2 所需增量步數(shù)不隨切向位移幅值發(fā)生改變，而考慮材料彈-塑性時，分析步2 所需增量步數(shù)隨著位移載荷的增大而增大。這是因為考慮材料彈-塑性時，由于塑性變形的存在，試件發(fā)生變形的區(qū)域會增大，計算的收斂速率減小，導致增量步數(shù)增加。

表2 不同條件下的增量步數(shù)Tab.2 Number of Incremental Steps Under Different Conditions

此外，還觀察到當材料為完全彈性時，求解一次仿真兩種材料所需的增量步數(shù)基本相同，而考慮材料彈-塑性時，steel 模型所需的增量步數(shù)明顯多于Al 模型所需的增量步數(shù)。這是因為相同載荷條件下，steel 模型的接觸壓力值較大，steel 模型率先產生塑性變形，收斂速率減小，因而所需的增量步數(shù)更多。

4.2 表面粗糙度對接觸壓力的影響

當法向載荷相同時，分別討論光滑接觸模型和粗糙接觸模型的接觸壓力分布。光滑接觸模型，外部法向載荷為20MPa 時，模型的接觸壓力分布，如圖6（a）所示。當下試件接觸區(qū)域輪廓為D=1.35 的分形曲線，法向載荷為20MPa 時，模型的接觸壓力分布，如圖6（b）所示。從圖中可以看出當接觸表面光滑時，接觸區(qū)域均勻接觸，當接觸區(qū)域為粗糙表面時，接觸發(fā)生在個別凸起的表面峰上。注意到，由于粗糙表面接觸中法向載荷被個別凸起的表面峰承擔，導致局部應力集中。因此，圖6（b）中的最大應力值遠大于圖6（a）中的應力最大值。為了進一步研究表面粗糙度對接觸表面接觸壓力分布的影響，對不同粗糙度接觸表面的接觸壓力分布進行了研究，如圖7 所示。根據(jù)圖1 所述，通過改變分形維數(shù)來獲得不同粗糙度的接觸表面。下試件接觸表面較為光滑時的接觸壓力分布，如圖7（a）所示。下試件接觸區(qū)相對粗糙時的接觸壓力分布，如圖7（b）所示。圖7（a）和圖7（b）所對應接觸模型的外部載荷和約束條件均相同。從圖7 可以看出，粗糙表面的存在導致接觸壓力分布為非光滑曲線，局部應力集中程度高。當表面粗糙度較大時，結合面上接觸壓力的分布是離散的。這是因為當接觸表面粗糙度較小時整個接觸表面完全接觸，但由于接觸表面凹凸不平引起接觸表面各個接觸點的變形量各不相同，從而導致整個接觸面上接觸壓力分布連續(xù)而不光滑以及局部應力集中。當表面粗糙度較大時，相同的法向載荷作用下，整個接觸區(qū)域只有個別凸起的表面峰發(fā)生接觸，從而引起接觸壓力分布離散化，并且局部應力集中嚴重，具體表現(xiàn)為個別接觸點的應力值非常大。

圖6 接觸區(qū)域的壓力分布（P=20MPa）Fig.6 Pressure Distribution in Contact Area（P=20MPa）

圖7 不同表面粗糙度下的接觸壓力分布Fig.7 Contact Pressure Distribution Under Different Surface Roughness

4.3 不同條件下的剪切摩擦力分布

接觸面的剪切摩擦力分布直接影響接觸面的磨損形貌和微動裂紋的初始位置，同時其影響因素眾多。為了盡可能全面地研究粗糙表面接觸模型的剪切摩擦力分布形式，仿真分析了不同接觸副材料、不同材料屬性、切向載荷幅值以及法向載荷幅值對切向剪切摩擦力分布的影響規(guī)律。

當法向載荷為P=20MPa，切向位移載荷為S=1μm，接觸副材料為steel 時，分別考慮完全彈性和彈-塑性時的剪切摩擦力分布，如圖8 所示。接觸副材料為Al 時分別考慮完全彈性和彈-塑性時的剪切摩擦力分布，如圖9 所示。圖中虛線和實線分別代表粗糙接觸模型和光滑接觸模型。從圖8、圖9 可以看出，當接觸副材料為完全彈性時，相比于光滑接觸模型，粗糙接觸接觸模型的接觸區(qū)域邊緣存在很大的應力集中，而考慮材料彈-塑性時，接觸區(qū)域邊緣的應力集中要小的多。這是因為考慮材料塑性時，相同法向載荷下接觸區(qū)域寬度明顯增加，而相應的接觸壓力值比材料完全彈性時小且波動幅度減小。

從圖8 和圖9 的對比不難發(fā)現(xiàn)，在相同的載荷條件下，彈性模量較小的Al 分別考慮完全彈性和彈-塑性時，接觸區(qū)域的剪切摩擦力波動很小，而彈性模量較大的steel 考慮材料彈-塑性時，接觸區(qū)域剪切摩擦力劇烈波動。法向載荷一定，接觸副材料為Al時，剪切摩擦力隨切向位移幅值的變化曲線，如圖10 所示。圖10表明材料為完全彈性或考慮彈-塑性時，隨著切向位移幅值的增大，剪切摩擦力都逐漸增大。此外，可以看出，材料完全彈性時接觸邊緣的剪切摩擦力大于考慮材料彈-塑性時接觸區(qū)域邊緣的剪切摩擦力。

圖8 steel 不同材料屬性時的剪切摩擦力分布Fig.8 Shear Force Distribution Under Different Material Properties of Steel

圖9 Al 不同材料屬性時的剪切摩擦力分布Fig.9 Shear Force Distribution Under Different Material Properties of Al

圖10 不同切向位移載荷下的剪切摩擦力分布（法向載荷P=20MPa）Fig.10 Shear Force Distribution under Different Tangential Displacement Loads（Normal Load P=20MPa）

切向位移載荷一定，接觸副材料為Al 時，剪切摩擦力隨法向載荷幅值的變化曲線，如圖11 所示。從圖11 中可以看出，當材料為完全彈性或考慮彈-塑性時，隨著法向載荷的增大，接觸區(qū)域寬度明顯增加，曲線分布形式變化不大。此外，從圖11 還可以看出，相同法向載荷下，考慮材料彈-塑性形時，剪切摩擦力峰值較小，而接觸區(qū)域較大。這是因為考慮材料塑性時接觸壓力峰值小，接觸區(qū)域寬度較大，而接觸壓力是剪切摩擦力的決定因素之一。

圖11 不同法向載荷下的剪切摩擦力分布（位移載荷S=1μm）Fig.11 Shear Force Distribution Under Different Normal Loads（Displacement Load S=1μm）

5 結論

利用有限元方法研究了微動接觸中分形粗糙表面的接觸應力分布規(guī)律，可以得出以下幾點：

（1）考慮材料彈-塑性時，接觸模型的接觸區(qū)域明顯增大，仿真一次所需的增量步數(shù)也明顯增大；當接觸副的材料不同時，其仿真一次所需的增量步數(shù)差異明顯。（2）接觸模型中粗糙表面的引入，導致接觸面上的接觸壓力分布為非光滑曲線，局部應力集中程度高；當表面粗糙度較大時，接觸面上接觸壓力的分布會發(fā)生離散化。（3）在粗糙接觸模型中，考慮材料彈-塑性時，不同接觸副材料的剪切摩擦力沿粗糙表面的分布差異明顯。