基于分形理論的兩粗糙表面接觸的黏滑摩擦模型

周 華，龍新華，孟 光

（上海交通大學(xué)機(jī)械系統(tǒng)與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240）

引 言

機(jī)械結(jié)構(gòu)中存在大量的連接界面，例如螺栓連接。這些連接界面受到法向緊固載荷后，能夠承受切向載荷。從微觀角度看，在切向載荷作用下，連接表面大量的微凸體接觸會(huì)發(fā)生黏著和滑移等行為，使連接界面呈現(xiàn)出復(fù)雜的非線性行為。建立兩粗糙表面接觸的切向黏滑摩擦模型一直是具有挑戰(zhàn)性的問(wèn)題。

兩粗糙表面的接觸模型可以分為統(tǒng)計(jì)模型和分形模型。Greenwood 等［1］假設(shè)微凸體高度滿足高斯分布，并基于統(tǒng)計(jì)方法建立了法向接觸模型（GW 模型）。不少學(xué)者在此模型上進(jìn)行了改進(jìn)，但是這種方法受到測(cè)量?jī)x器的分辨率和采樣長(zhǎng)度的影響。因此，Majumda 等［2］用兩個(gè)分形參數(shù)表征了表面的形貌，并提出了分形接觸模型（MB 模型），克服了統(tǒng)計(jì)學(xué)方法的不足。Wang 等［3］、Yan 等［4］和劉文威［5］對(duì)MB 模型進(jìn)行了改進(jìn)，使分形接觸模型更加合理，然而這些模型并未考慮微凸體的相互作用。Zhao 等［6］考慮了微凸體相互作用造成的基底面下降，并應(yīng)用在GW 統(tǒng)計(jì)模型中。此外，張學(xué)良等［7］、溫淑花等［8］和田紅亮等［9］基于分形理論建立了結(jié)合面的法向接觸剛度模型。Zhou 等［10］考慮微凸體相互作用造成的基底面下降，建立了兩粗糙表面的法向接觸模型。

對(duì)于兩粗糙表面接觸的黏滑摩擦研究，Mind?lin［11］推導(dǎo)了在給定恒定法向載荷時(shí)，球接觸的切向載荷與切向位移之間的非線性關(guān)系。文中認(rèn)為圓形的接觸區(qū)域由兩部分組成，接觸的中心區(qū)域?yàn)镾tick黏著區(qū)，接觸的外層為Slip 滑移區(qū)。當(dāng)切向載荷增大時(shí)，滑移區(qū)會(huì)不斷增大，直至發(fā)生宏觀滑動(dòng)。Mindlin 等［12］給出了處于彈性接觸的球體在循環(huán)加載時(shí)的響應(yīng)。?dfalk 等［13］把彈性接觸擴(kuò)展成彈性?塑性接觸情況。Bj?rklund［14］基于單個(gè)微凸體的Mindlin 模型和GW 統(tǒng)計(jì)模型分析了接觸粗糙表面的黏滑現(xiàn)象，但僅考慮了微凸體的彈性變形。Er?iten 等［15?17］考慮了微凸體的彈塑性變形，并且根據(jù)不同的摩擦系數(shù)模型建立了黏滑摩擦模型。Wang等［18?19］考慮了微凸體的彈性和完全塑性變形，但也是根據(jù)統(tǒng)計(jì)模型分析了粗糙表面的切向接觸特性。Jamshidi 等［20］考慮了微凸體的斜接觸，即切向和法向變形耦合，提出了改進(jìn)的兩粗糙表面接觸模型。該模型直接分析兩個(gè)粗糙表面接觸，無(wú)需等效成一個(gè)光滑表面和一個(gè)等效粗糙表面，但也是基于GW統(tǒng)計(jì)模型。Song 等［21］嘗試基于分形模型建立界面黏滑摩擦模型，但并未考慮微凸體的相互作用。

本文采用文獻(xiàn)［10］中的法向接觸模型，并確定在給定法向預(yù)緊載荷下微接觸截面積的概率密度函數(shù)。同時(shí)，本文基于Mindlin 模型和Masing 準(zhǔn)則獲得了單個(gè)微凸體的切向模型，然后根據(jù)分形理論對(duì)微接觸截面積進(jìn)行積分，建立了整個(gè)粗糙表面接觸的黏滑摩擦模型，并研究了分形參數(shù)D和G以及法向預(yù)緊力對(duì)接觸表面能量耗散的影響。

1 兩粗糙表面法向接觸模型

1.1 考慮微凸體相互作用的分形接觸模型

兩個(gè)粗糙表面接觸可以等效成一個(gè)光滑表面和一個(gè)等效的粗糙表面［22］，等效表面的等效彈性模量為：

式中E1和E2分別為兩表面材料的彈性模量，ν1和ν2分別是兩表面材料的泊松比。

圖1為一個(gè)光滑剛性表面和一個(gè)等效粗糙表面接觸的示意圖。δn表示微凸體的總變形量。作用在其他微凸體上的壓力會(huì)對(duì)給定微凸體產(chǎn)生額外的變形，用微凸體平均高度的下降距離u表示微凸體相互作用產(chǎn)生的影響，此時(shí)單個(gè)微凸體真實(shí)的干涉量ω=δn-u。

圖1 一個(gè)剛性光滑表面與粗糙表面的接觸Fig.1 Contact of a rough surface with a rigid smooth flat

微凸體的總變形量δn可以由粗糙表面的形貌獲得。Yan 等［4］使用改進(jìn)的兩變量W?M 函數(shù)來(lái)表征粗糙表面的形貌。他們假設(shè)接觸表面是各向同性表面，并且推導(dǎo)了等效的單變量等式。

式中Ds（2

式中 int［…］表示取中括號(hào)數(shù)值的最大整數(shù)。若微凸體的截半徑為r′，則微凸體波形的最長(zhǎng)波長(zhǎng)為2r′?？梢院侠淼丶僭O(shè)微接觸力主要是由于基波長(zhǎng)的微凸體尖端的小變形造成的，并且相應(yīng)的頻率指數(shù)為n0=ln(L/2r′)/lnγ［4］。代入到式（2），可以得到基波輪廓曲線函數(shù)z0（x），表示為：

δn定義為關(guān)于x軸的余弦函數(shù)的峰值高度，寫(xiě)為：

式中c1=22-Dπ(D-2)/2GD-1(lnγ)1/2，a′為微接觸截面積，即a′=πr′2。

由于微凸體相互作用造成的微凸體平均高度下降的位移u可以表示為［6］：

式中Fn為總法向載荷，Aa為名義接觸面積，F(xiàn)as為一個(gè)微凸體上的力。

考慮微凸體相互作用后，單個(gè)微凸體的真實(shí)干涉量可以表示為：

需要指出的是，當(dāng)微凸體處于不同變形階段時(shí)，上式中的Fas是不同的。

研究表明，微接觸截面積的分布滿足地球島嶼面積分布，微接觸截面積的尺寸分布函數(shù)n(a′)，或稱為概率密度函數(shù)，可以表示為：

式中a′l為最大的微凸體截面積，ψ為域擴(kuò)展因子，與分形維數(shù)相關(guān)，表示為：

對(duì)等式（9）求解并擬合數(shù)據(jù)，可以直觀表示為：

微凸體的曲率半徑R可以從等式（4）獲得：

式中c2=[22-Dπ2+D/2GD-1(lnγ)1/2]-1。

微凸體有三個(gè)變形階段：彈性變形、彈塑性變形和完全塑性變形。ωyc定義為微凸體從彈性變形階段轉(zhuǎn)為彈塑性變形階段的臨界屈服變形值，ωpc定義為微凸體從彈塑性變形階段轉(zhuǎn)為完全塑性變形階段的臨界變形值，可以表示為［22?23］：

式中H表示較軟材料的硬度，H≈2.8σb；K為硬度系數(shù)［24］，且K= 0.4645 + 0.3141ν+ 0.1943ν2表示初始屈服點(diǎn)的平均壓力，=1.1σy；表示完全塑性流動(dòng)轉(zhuǎn)變點(diǎn)處的均壓［23］，=2.8σy。

把等式（11）代入等式（12），臨界屈服變形量為：

（i）彈性變形階段

微凸體在此階段的微接觸面積ae和接觸載荷Fe可以表示為［2］：

此時(shí)，真實(shí)干涉量ω可以表示為：

根據(jù)等式（17）中a′和ω的關(guān)系，臨界微接觸屈服截面積a′yc可以由等式（12）中的ωyc求得；類似的，臨界微接觸截面積a′pc可以由等式（13）中的ωpc求得。

因?yàn)閏4與法向載荷Fn相關(guān)，所以這兩個(gè)臨界接觸面積并不是常值，而是與法向載荷相關(guān)。

（ii）彈塑性階段

微凸體在此階段的微接觸面積aep和接觸載荷Fep可以表示為［5］：

式中b和s為常數(shù)，可以根據(jù)微接觸面積和載荷在彈性、彈塑性和完全塑性變形的連續(xù)性條件計(jì)算出這兩個(gè)參數(shù)。

此時(shí)，真實(shí)干涉量ω可以表示為：

（iii）完全塑性階段

微凸體在此階段的微接觸面積ap和接觸載荷Fp可以表示為［22］：

此時(shí)，真實(shí)干涉量ω可以表示為：

1.2 法向總載荷

單個(gè)微凸體在彈性、彈塑性和完全塑性階段的法向載荷分別用Fe，F(xiàn)ep和Fp表示。根據(jù)尺寸分布函數(shù)n(a′)，這三個(gè)階段總的接觸載荷分別用Fne，F(xiàn)nep和Fnp表示，由下面等式給出：

需要指出的是，這些積分的被積函數(shù)是隱式表達(dá)式，三個(gè)階段的ω和a′的關(guān)系分別由等式（17），（23）和（26）給出，并基于Maltab 軟件的fzero 函數(shù)求解這些積分。

連接界面總的法向載荷是三個(gè)變形階段法向載荷的總和，可以表示為：

2 兩粗糙表面的黏滑摩擦建模

在法向載荷的作用下，兩接觸表面被壓緊，受到較小的往復(fù)切向載荷或切向位移時(shí)，會(huì)產(chǎn)生往復(fù)的切向位移，出現(xiàn)滯回現(xiàn)象。從單個(gè)微凸體角度看，有的微凸體接觸處于黏著狀態(tài)，而有的微凸體接觸處于滑移狀態(tài)。當(dāng)切向載荷不斷增大時(shí)，短的微凸體受到的法向載荷比高的微凸體小，所以更容易先滑動(dòng)。同理，微凸體接觸截面積越小，越容易滑動(dòng)。在首次加載時(shí)，接觸截面積大的黏著微凸體構(gòu)成了切向剛度，而接觸截面積小的滑移微凸體導(dǎo)致了能量耗散。卸載時(shí)，切向載荷減小，而預(yù)緊力不變，即極限切向載荷不變，所以一些在首次加載時(shí)滑移的微凸體會(huì)轉(zhuǎn)變?yōu)轲ぶ鵂顟B(tài)。

圖2和3 展示了微凸體在不同接觸狀態(tài)時(shí)微接觸截面積概率密度函數(shù)的區(qū)間。首次加載時(shí)，接觸的微凸體狀態(tài)可以分成兩組。而卸載時(shí)，接觸的微凸體狀態(tài)分成三組［17］。

（i）微凸體在卸載和加載時(shí)都滑移。

（ii）微凸體在加載時(shí)滑移，但卸載時(shí)處于黏著狀態(tài)。

（iii）微凸體在卸載和加載時(shí)都處于黏著狀態(tài)。

圖2和3 中，a′s是最小的微凸體接觸截面積，可以近似為0。a′0，a′1和a′2都是臨界接觸截面積，用于區(qū)別首次加載與卸載時(shí)微凸體的接觸狀態(tài)，下文給出表達(dá)式。

圖2 首次加載時(shí)，微凸體在不同接觸狀態(tài)時(shí)微接觸截面積概率密度函數(shù)的區(qū)間Fig.2 Categories of asperities when initial loading based on the truncated microcontact area distribution

圖3 卸載時(shí)，微凸體在不同接觸狀態(tài)時(shí)微接觸截面積概率密度函數(shù)的區(qū)間。Fig.3 Categories of asperities when unloading based on the truncated microcontact area distribution

Mindlin［12］給出了切向載荷與切向相對(duì)位移之間的關(guān)系：

式中δ為相對(duì)切向位移，T為切向載荷，r為接觸半徑，F(xiàn)as為單個(gè)微凸體承受的法向力，μ為摩擦系數(shù)，G′為等效剪切模量。

根據(jù)等式（31），切向載荷與位移之間的關(guān)系為：

無(wú)量綱切向相對(duì)位移用彈性階段的接觸面積與法向載荷表示為：

當(dāng)微凸體處于彈塑性階段時(shí)，等式（34）中的比值r/Fas會(huì)不同，但是至多導(dǎo)致10% 的誤差［17］。為了便于計(jì)算，本文都采用了彈性接觸時(shí)的表達(dá)式。 摩擦系數(shù)μ采用CFC 恒定值模型，并取μ= 0.2。

當(dāng)切向載荷不足以引起整體滑動(dòng)時(shí)，處于滑移狀態(tài)的微凸體，其切向載荷滿足Coulomb 摩擦定律；處于黏著狀態(tài)的微凸體，其切向載荷與相對(duì)位移滿足Mindlin 理型。此時(shí)，可以得到：

（i）首次加載

式中Tas為單個(gè)微凸體的切向載荷，a′0為臨界微接觸截面積，對(duì)應(yīng)的臨界變形量ω0=4G′δ/(μE)，由于ω與a′存在隱式關(guān)系，無(wú)法直接給出顯式表達(dá)式，但可以由等式（17），（23），（26）確定。

根據(jù)微凸體截面積的尺寸分布函數(shù)n(a′)，面?面接觸的總切向載荷為：

需要指出的是，等式（36）需要進(jìn)行分段積分。通過(guò)判斷a′0，a′pc和a′yc的大小關(guān)系，確定積分區(qū)間。例如當(dāng)a′pc

（ii）卸載

式中a′1和a′2為臨界接觸截面積，對(duì)應(yīng)的臨界變形量為ω1=2G′(δmax-δ) /(μE) 和ω2=4G′δmax/(μE)，ω與a′的關(guān)系由等式（17），（23），（26）確定。

同理，計(jì)算卸載時(shí)面?面接觸的總切向載荷時(shí)，需要判斷a′1，a′2，a′pc和a′yc的大小關(guān)系。例如當(dāng)a′1≤a′pc

（iii）再次加載

微凸體的卸載方程和再次加載方程滿足Mas?ing 映射準(zhǔn)則，可以得到再次加載時(shí)面?面接觸的總切向載荷：

（iv）能量耗散

在周期載荷作用下，面?面接觸的力?位移的曲線形成了封閉區(qū)域，封閉區(qū)域的面積即為一個(gè)周期內(nèi)的能量耗散，可以表示為：

3 數(shù)值仿真

為了進(jìn)一步研究所提出的基于分形理論的切向黏滑摩擦模型的特性，開(kāi)展了數(shù)值仿真研究。切向載荷與相對(duì)位移關(guān)系曲線的計(jì)算流程如圖4所示。無(wú)量綱化公式為：

圖4 切向總載荷的計(jì)算流程圖Fig.4 Calculation of total tangential load

仿真參數(shù)為：E1=E2= 210 GPa，ν1=ν2=0.3，H= 2.744 GPa，σy= 785 MPa，Aa= π ×（10/1000）2m2，μ= 0.2。其他參數(shù)分成3 組，用于對(duì)比研究不同參數(shù)對(duì)切向黏滑摩擦模型的影響。

（1）G= 1×10-7m，F(xiàn)n= 1000 N，D分別取1.60，1.62，1.64。

（2）D= 1.60，F(xiàn)n= 1000 N，G分別取0.8×10-7，1×10-7，1.2×10-7m。

（3）D= 1.60，G= 1×10-7m，F(xiàn)n分別取1000，1050，1100 N。

圖5展示了不同的分形維數(shù)D時(shí)，連接界面的遲滯曲線。從圖中看出，在相同切向載荷下，D增大時(shí)，遲滯曲線所圍面積增大，即一個(gè)周期的能量耗散增大。D= 1.60，1.62 和1.64 的能量耗散分別為1.202×10-4，1.408×10-4和1.591×10-4J。分形維數(shù)D反映表面輪廓的復(fù)雜性和自相似性，一般來(lái)說(shuō)，D越大，表面越光滑，實(shí)際接觸面積增加，滑移微凸體的個(gè)數(shù)增加，所以能量耗散也增大，該趨勢(shì)與文獻(xiàn)［25?26］中的結(jié)果一致。

圖5 不同D 時(shí)，G = 1×10-7 m，F(xiàn)n = 1000 N 時(shí)的切向遲滯曲線Fig.5 Hysteresis loop for different fractal dimension D with G = 1×10-7 m，F(xiàn)n = 1000 N

不同分形粗糙度參數(shù)G時(shí)，連接界面的遲滯曲線如圖6所示。從圖中可以看出，在相同切向載荷下，G增大時(shí)，遲滯曲線所圍面積減小，即能量耗散減小。G= 0.8×10-7，1×10-7和1.2×10-7的能量耗散分別為1.370×10-4，1.202×10-4和1.061×10-4J。分形粗糙度參數(shù)G反映粗糙表面的垂直方向的光滑程度，G越大，則表面越粗糙。所以圖6說(shuō)明表面越粗糙，能量耗散越小，與前文結(jié)果一致。

圖6 不同G 時(shí)，D = 1.60，F(xiàn)n= 1000 N 時(shí)的切向遲滯曲線Fig.6 Hysteresis loop for different fractal roughness G with D = 1.60，F(xiàn)n = 1000 N

圖7展示了不同法向緊固載荷時(shí)，連接界面的遲滯曲線。從圖中可以看出，在相同切向載荷下，F(xiàn)n增大時(shí)，遲滯曲線所圍面積減小，即能量耗散減小。Fn= 1000，1050，1100 N 時(shí)的能量耗散分別為1.202×10-4，8.447×10-5和6.287×10-5J。這是因?yàn)榉ㄏ蜉d荷的增大會(huì)導(dǎo)致微凸體滑移的最大切向載荷增大，在切向外載荷不變時(shí)，滑移微凸體的個(gè)數(shù)會(huì)減少，所以能量耗散減小，該趨勢(shì)與文獻(xiàn)［17，26］的結(jié)果一致。

圖7 不同F(xiàn)n時(shí)，D=1.60，G=1×10-7 m 時(shí)的切向遲滯曲線Fig.7 Hysteresis loop for different normal load Fn with D =1.60，G = 1×10-7 m

4 結(jié) 論

本文考慮了微凸體相互作用對(duì)基底面的下降，并基于分形理論建立了兩粗糙表面的法向接觸模型，從而確定在給定法向預(yù)緊載荷下微接觸截面積的概率密度函數(shù)。然后根據(jù)Mindlin 模型、Masing準(zhǔn)則和分形理論，建立了兩粗糙表面接觸的切向黏滑摩擦模型。最后開(kāi)展了數(shù)值仿真研究，分析了不同參數(shù)對(duì)系統(tǒng)能量耗散的影響。結(jié)論如下：

（1）分形維數(shù)D越大，能量耗散越大；分形粗糙度參數(shù)G越大，能量耗散越小，即表面越粗糙，能量耗散越小。

（2）法向預(yù)緊載荷增大時(shí)，微凸體滑移的最大切向載荷會(huì)增大，在切向外載荷不變時(shí)，滑移微凸體的個(gè)數(shù)會(huì)減少，所以能量耗散越小。