微動(dòng)接觸中分形粗糙表面的接觸應(yīng)力研究

阮曉光，麻詩(shī)韻，李 玲，蔡安江

（西安建筑科技大學(xué)機(jī)電工程學(xué)院，陜西 西安 710055）

1 引言

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

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

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

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

2 二維粗糙表面的表征

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

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

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

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

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

3 仿真分析

3.1 有限元模型

基于微動(dòng)實(shí)驗(yàn)?zāi)Ｐ徒⒍S柱面/平面接觸有限元模型，如圖3 所示。在圓柱試件的頂面施加均布的法向載荷P，并通過(guò)多點(diǎn)約束選項(xiàng)（MPC）在頂面中心點(diǎn)施加切向位移載荷S。下試件兩側(cè)約束水平方向位移，底部約束豎直方向位移。

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

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

圖4 加載歷程Fig.4 Loading History

3.2 材料屬性

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

表1 材料屬性Tab.1 Material Properties

圖5 真實(shí)應(yīng)力和塑性應(yīng)變的關(guān)系曲線Fig.5 The Relationship Curve of Between True Stress and Plastic Strain

4 結(jié)果與討論

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

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

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

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

4.2 表面粗糙度對(duì)接觸壓力的影響

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

圖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 不同條件下的剪切摩擦力分布

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

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

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

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

圖9 Al 不同材料屬性時(shí)的剪切摩擦力分布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 時(shí)，剪切摩擦力隨法向載荷幅值的變化曲線，如圖11 所示。從圖11 中可以看出，當(dāng)材料為完全彈性或考慮彈-塑性時(shí)，隨著法向載荷的增大，接觸區(qū)域?qū)挾让黠@增加，曲線分布形式變化不大。此外，從圖11 還可以看出，相同法向載荷下，考慮材料彈-塑性形時(shí)，剪切摩擦力峰值較小，而接觸區(qū)域較大。這是因?yàn)榭紤]材料塑性時(shí)接觸壓力峰值小，接觸區(qū)域?qū)挾容^大，而接觸壓力是剪切摩擦力的決定因素之一。

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

5 結(jié)論

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

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