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

周 華，龍新華，孟 光

（上海交通大學機械系統與振動國家重點實驗室，上海 200240）

引 言

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

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

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

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

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

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

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

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

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

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

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

式中Ds（2

式中 int［…］表示取中括號數值的最大整數。若微凸體的截半徑為r′，則微凸體波形的最長波長為2r′。可以合理地假設微接觸力主要是由于基波長的微凸體尖端的小變形造成的，并且相應的頻率指數為n0=ln(L/2r′)/lnγ［4］。代入到式（2），可以得到基波輪廓曲線函數z0（x），表示為：

δn定義為關于x軸的余弦函數的峰值高度，寫為：

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

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

式中Fn為總法向載荷，Aa為名義接觸面積，Fas為一個微凸體上的力。

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

需要指出的是，當微凸體處于不同變形階段時，上式中的Fas是不同的。

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

式中a′l為最大的微凸體截面積，ψ為域擴展因子，與分形維數相關，表示為：

對等式（9）求解并擬合數據，可以直觀表示為：

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

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

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

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

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

（i）彈性變形階段

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

此時，真實干涉量ω可以表示為：

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

因為c4與法向載荷Fn相關，所以這兩個臨界接觸面積并不是常值，而是與法向載荷相關。

（ii）彈塑性階段

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

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

此時，真實干涉量ω可以表示為：

（iii）完全塑性階段

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

此時，真實干涉量ω可以表示為：

1.2 法向總載荷

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

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

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

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

在法向載荷的作用下，兩接觸表面被壓緊，受到較小的往復切向載荷或切向位移時，會產生往復的切向位移，出現滯回現象。從單個微凸體角度看，有的微凸體接觸處于黏著狀態，而有的微凸體接觸處于滑移狀態。當切向載荷不斷增大時，短的微凸體受到的法向載荷比高的微凸體小，所以更容易先滑動。同理，微凸體接觸截面積越小，越容易滑動。在首次加載時，接觸截面積大的黏著微凸體構成了切向剛度，而接觸截面積小的滑移微凸體導致了能量耗散。卸載時，切向載荷減小，而預緊力不變，即極限切向載荷不變，所以一些在首次加載時滑移的微凸體會轉變為黏著狀態。

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

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

（ii）微凸體在加載時滑移，但卸載時處于黏著狀態。

（iii）微凸體在卸載和加載時都處于黏著狀態。

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

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

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

Mindlin［12］給出了切向載荷與切向相對位移之間的關系：

式中δ為相對切向位移，T為切向載荷，r為接觸半徑，Fas為單個微凸體承受的法向力，μ為摩擦系數，G′為等效剪切模量。

根據等式（31），切向載荷與位移之間的關系為：

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

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

當切向載荷不足以引起整體滑動時，處于滑移狀態的微凸體，其切向載荷滿足Coulomb 摩擦定律；處于黏著狀態的微凸體，其切向載荷與相對位移滿足Mindlin 理型。此時，可以得到：

（i）首次加載

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

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

需要指出的是，等式（36）需要進行分段積分。通過判斷a′0，a′pc和a′yc的大小關系，確定積分區間。例如當a′pc

（ii）卸載

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

同理，計算卸載時面?面接觸的總切向載荷時，需要判斷a′1，a′2，a′pc和a′yc的大小關系。例如當a′1≤a′pc

（iii）再次加載

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

（iv）能量耗散

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

3 數值仿真

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

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

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

（1）G= 1×10-7m，Fn= 1000 N，D分別取1.60，1.62，1.64。

（2）D= 1.60，Fn= 1000 N，G分別取0.8×10-7，1×10-7，1.2×10-7m。

（3）D= 1.60，G= 1×10-7m，Fn分別取1000，1050，1100 N。

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

圖5 不同D 時，G = 1×10-7 m，Fn = 1000 N 時的切向遲滯曲線Fig.5 Hysteresis loop for different fractal dimension D with G = 1×10-7 m，Fn = 1000 N

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

圖6 不同G 時，D = 1.60，Fn= 1000 N 時的切向遲滯曲線Fig.6 Hysteresis loop for different fractal roughness G with D = 1.60，Fn = 1000 N

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

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

4 結 論

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

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

（2）法向預緊載荷增大時，微凸體滑移的最大切向載荷會增大，在切向外載荷不變時，滑移微凸體的個數會減少，所以能量耗散越小。