粗糙度各向異性對潤滑接觸面摩擦的影響*

張園園 李元慶，2 黃 培，2 付紹云，2

（1. 重慶大學(xué)航空航天學(xué)院 重慶 400044； 2. 輸配電裝備及系統(tǒng)安全與新技術(shù)國家重點實驗室 重慶 400044）

機械零部件之間的摩擦行為往往伴隨著能量的耗散， 使得一些關(guān)鍵零部件的傳動效率以及疲勞性能降低[1］。 齒輪、 軸承、 凸輪等摩擦副是重要的機械基礎(chǔ)件， 其潤滑接觸特性成為決定部件及裝備性能與可靠性的關(guān)鍵因素[2］。 摩擦副之間適當(dāng)?shù)臐櫥瑺顟B(tài)一方面阻止兩接觸面的直接接觸， 一方面起到冷卻效果， 從而改善兩接觸體的摩擦行為。 深入理解繼而控制、 優(yōu)化摩擦副之間的摩擦行為， 能夠有效控制機械零部件運行過程中的能量損失[3］， 降低接觸面之間的磨損，延長摩擦副的服役壽命。

摩擦副之間的接觸面并不是絕對光滑的， 接觸面的表面粗糙度通常與潤滑油膜的厚度具有相同數(shù)量級[4］， 對潤滑接觸面的摩擦行為有重要影響[5-10］。Stribeck 曲線（摩擦因數(shù)隨“λ比” 變化的曲線） 可有效描述粗糙潤滑接觸面的摩擦行為[11］， 其中“λ比” 被定義為潤滑油膜厚度與粗糙度初始幅值的比值。 早期的Stribeck 曲線主要是在輕載情況下通過試驗的方式測得[12-14］， 其中粗糙度對摩擦因數(shù)的影響只與粗糙度初始幅值相關(guān)， 并未涉及粗糙度的形貌變化。 目前大多數(shù)關(guān)鍵機械零部件經(jīng)常在高負(fù)載情況下運行， 單純考慮粗糙度初始幅值對厘清粗糙度與摩擦因數(shù)的定量關(guān)系明顯不足。

隨著流體動力潤滑理論的發(fā)展， 使得重載情況下的潤滑粗糙接觸行為數(shù)值分析成為可能[15］。 研究發(fā)現(xiàn)潤滑條件下摩擦副的摩擦行為與多種因素有關(guān)。 例如圓錐滾子軸承擋邊面與滾子端面橢圓接觸的數(shù)值模擬研究發(fā)現(xiàn)， 接觸面的彈性變形對摩擦因數(shù)的影響不可忽略[16］。 同時， 潤滑劑的本構(gòu)對接觸面之間的摩擦行為也具有顯著的影響[17-18］。 除了上述因素， 摩擦副接觸面的粗糙度對于接觸面間的摩擦行為也有重要影響， 接觸面摩擦行為主要由油膜厚度與粗糙度參數(shù)決定。 油膜厚度作為主要影響因素之一， 與粗糙度幅值、 波長以及各向異性也息息相關(guān)[19-22］。 輕載（協(xié)調(diào)性接觸） 情況下的潤滑接觸研究表明： 當(dāng)“λ比” 等于3 時， 全膜潤滑開始向混合潤滑過渡[23］。最近基于數(shù)值模擬研究發(fā)現(xiàn)， 在重載（非協(xié)調(diào)性接觸） 情況下， 接觸面粗糙度會發(fā)生彈性形變， 其形變量依賴于工況條件以及粗糙度形貌[24-25］（該結(jié)果也被光彈實驗驗證[26］）。 這就意味著傳統(tǒng)的參數(shù)“λ比” 已經(jīng)不再適用于對重載情況下潤滑接觸面的摩擦因數(shù)變化情況進(jìn)行描述， 亟需找到一個新的參數(shù)對重載下的潤滑接觸面的摩擦因數(shù)變化情況進(jìn)行描述。

復(fù)雜的表面粗糙度可以分解成一系列正弦粗糙度分量疊加的形式， 本文作者基于瞬態(tài)彈性流體動力潤滑理論， 構(gòu)建正弦表面粗糙度， 并采用多重網(wǎng)格算法求解， 系統(tǒng)研究表面粗糙度各向異性對潤滑接觸面摩擦因數(shù)的影響， 并進(jìn)行影響量化分析。 文中為解決表面粗糙度的引入造成控制方程的不連續(xù)性， 在多重網(wǎng)格算法中使用ALCOUFFE 等[27］構(gòu)造連續(xù)插值算子的方式提高算法的穩(wěn)健性。

1 數(shù)值模型

1.1 量綱一化方程

潤滑接觸問題首先需要推導(dǎo)待求解的控制方程組[28］， 考慮時變效應(yīng)的量綱一化雷諾方程為

式中：P為量綱一壓力；H為接觸面間的量綱一油膜厚度為潤滑油量綱一密度；T為量綱一時間；X為量綱一坐標(biāo)， 其方向為接觸面的相對運動方向；Y為垂直于相對運動方向的量綱一坐標(biāo)。

量綱一化的具體過程參見文獻(xiàn)[28］。 為簡化雷諾方程， 使用了參數(shù)為潤滑油的量綱一動力黏度； 變量其中為接觸面1 和2 的平均線速度，η0為環(huán)境黏度，Rx為接觸點的當(dāng)量曲率半徑，ah與ph分別為Hertz 接觸半寬以及最大接觸壓力[28］，ah＝其中F是法向載荷，E′是等效彈性模量， 其定義為（1- v22）／E2。 同時雷諾方程的邊界條件為

式中：Xa、Xb、Ya、Yb分別代表接觸區(qū)域的邊界。

在該模型中， 采用Dowson 和Higginson 的密度-壓力方程[29］來描述潤滑油的密度變化情況； 使用Roelands 的黏度-壓力方程[30］來描述潤滑油的黏度與壓力的關(guān)系。

考慮表面粗糙度的量綱一油膜厚度方程為

式中： 指數(shù)項κ ＝10-10｛max[0，（X-X＾）／（λx／ah）］2｝用于加速數(shù)值計算程序的收斂， 其中為未變形粗糙度量綱一幅值，λx／ah與λy／ah分別為X和Y坐標(biāo)方向的粗糙度量綱一波長。

與此同時， 采用參數(shù)r來表示粗糙度各向異性程度， 即：

當(dāng)力平衡方程式（6） 滿足時， 公式（3） 中的法向剛體位移H0即可確定。

在全膜潤滑區(qū)域， 泊肅項與楔形項均對潤滑剪切力有貢獻(xiàn)， 對于牛頓流體， 當(dāng)滑動速度適中時摩擦力主要由楔形項決定[31］， 即：

式中：S為滑滾比， 即滑動速度與卷吸速度之比。基于公式（7）， 相對摩擦因數(shù)可以定義為

式中： 下標(biāo)rr 與ss 分別代表粗糙與光滑接觸兩種情況。

對于點接觸問題， BOSMA 和MOES[32］為了進(jìn)一步減少變量數(shù)目， 分別定義了2 個量綱一化參數(shù)來表征載荷、 速度與材料參數(shù)的影響， 即Moes 量綱一化參數(shù)M和L[32］：

式中：α為黏度-壓力系數(shù)。

1.2 求解算法

對上述控制方程組進(jìn)行離散。 量綱一化方程將在計算區(qū)域內(nèi)的均勻網(wǎng)格上采用二階近似的方式進(jìn)行離散， 即： 泊肅項使用二階中心差分格式， 楔形項與瞬態(tài)項使用二階向前差分格式。 同時使用多重網(wǎng)格法來加快迭代收斂速度， 使用多重積分法（MLMI） 來計算彈性形變。 多重網(wǎng)格算法的實施可參考文獻(xiàn)[28］。油膜厚度公式（3） 中的粗糙度項Ra波動劇烈， 使得傳統(tǒng)的EHL 多重網(wǎng)格求解器不足以應(yīng)對粗糙接觸問題， 因此文中提出了基于Alcouffe 算法的粗網(wǎng)格構(gòu)造方法來提高程序的穩(wěn)健性。 Alcouffe 粗網(wǎng)格構(gòu)造方法的具體實施過程見文獻(xiàn)[33］。

接觸模型量綱一化計算域選取為X∈[-2.5，1.5］ 和Y∈[-2.0， 2.0］ 的矩形區(qū)域， 同時最細(xì)網(wǎng)格上的節(jié)點數(shù)確定為513×513。 計算模型中的離散時間步長與空間域步長相同， 即0.007 812 5。 粗糙度設(shè)置的初始位置為X＝-2.5， 同時保證計算時間足夠長使得計算結(jié)果展現(xiàn)出周期性。為了避免摩擦因數(shù)比的分母出現(xiàn)0 的情況， 使用一個較小的滑滾比S＝0.02。 該模型數(shù)值求解的過程如圖1 所示。

圖1 TEHL 模型數(shù)值求解流程Fig.1 Flow of the numerical simulation process for the TEHL model

2 計算結(jié)果與討論

為討論表面粗糙度各向異性對相對摩擦因數(shù)的影響規(guī)律， 分別研究縱向粗糙度（λx ＞λy）， 橫向粗糙度（λy ＞λx） 以及純縱向粗糙度（λx ＝∞） 對摩擦副接觸面摩擦行為的影響。 所需數(shù)值模擬參數(shù)如表1 所示。

表1 數(shù)值模型參數(shù)Table 1 Parameters for the numerical model

2.1 數(shù)值求解

表2 所示為載荷在M ＝1 000 與L ＝10 的情況下，表面粗糙度波長分別為λx／ah＝1，λy／ah＝0.5， 中心油膜厚度與粗糙度初始幅值之比為3， 即Hc／Ai＝3時， 經(jīng)過30 個“V” 循環(huán)后， 不同網(wǎng)格節(jié)點數(shù)下，經(jīng)典多重網(wǎng)格算法[28］與Alcouffe 粗網(wǎng)格算子實施后的算法， 其計算誤差的對比。 可以看出， 改進(jìn)算法的計算精度在各個網(wǎng)格節(jié)點情況下都要比傳統(tǒng)算法高， 即雷諾方程與力平衡方程的殘差均比傳統(tǒng)算法小3 ～4個數(shù)量級。 尤其當(dāng)網(wǎng)格節(jié)點數(shù)為1 025×1 025 時， 經(jīng)典多重網(wǎng)格算法不收斂， 改進(jìn)后的算法依舊給出了較高的計算精度。 因此， 使用Alcouffe 方法改進(jìn)后的多重網(wǎng)格算法在求解粗糙表面接觸問題時具有更高的穩(wěn)健性。

表2 經(jīng)典多重網(wǎng)格算法與實施Alcouffe 粗網(wǎng)格構(gòu)造法的多重網(wǎng)格算法計算誤差對比Table 2 Comparison of calculation errors between the classical multilevel method and multilevel method implementing Alcouffe’s coarse grid construction

網(wǎng)格密度的選取與數(shù)值解的精度以及計算時間息息相關(guān)， 網(wǎng)格密度過大， 計算精度不夠， 將導(dǎo)致計算結(jié)果失真； 網(wǎng)格密度過小， 計算時間長， 計算代價太大。 表3 所示為相對摩擦因數(shù)與網(wǎng)格節(jié)點數(shù)的關(guān)系，結(jié)果表明當(dāng)計算區(qū)域內(nèi)的節(jié)點數(shù)設(shè)置為513×513 時，具有優(yōu)于1%的計算精度以及較小的計算代價， 因此文中的網(wǎng)格節(jié)點數(shù)選擇為513×513。

表3 網(wǎng)格節(jié)點數(shù)對相對摩擦因數(shù)的影響Table 3 Effects of the mesh point on the relative friction coefficient

將文中數(shù)值模型計算結(jié)果與文獻(xiàn)[24］擬合的計算結(jié)果進(jìn)行對比來驗證數(shù)值模型的正確性。 圖2 所示為文中求解的粗糙度相對幅值變化與文獻(xiàn)[24］中擬合出的粗糙度幅值變形公式的對比， 結(jié)果顯示文中數(shù)值模型計算的結(jié)果全部落在了文獻(xiàn)擬合的主曲線上，說明提出的數(shù)值模型能夠有效用于潤滑粗糙接觸面的摩擦行為分析。

圖2 粗糙度幅值變化與f（r）的關(guān)系：當(dāng)前模型與文獻(xiàn)[24］結(jié)果對比Fig.2 The relative deformed amplitude as a function of f（r）： results comparison of the current model with those of Ref. [24］

2.2 瞬時解

圖3 所示為摩擦副接觸面表面粗糙度分別表現(xiàn)出各向同性 （λx／ah＝λy／ah＝0.5）、 縱向 （λx／ah＝1，λy／ah＝0.5）、 橫向（λx／ah＝0.5，λy／ah＝1） 時， 相同載荷（M ＝1 000，L ＝10 和Hc／Ai＝2） 下的壓力、油膜厚度分布與相對摩擦因數(shù)隨時間的變化情況。 圖3 左圖顯示由于接觸面形貌發(fā)生了改變， 對應(yīng)的油膜厚度分布發(fā)生變化， 此時壓力分布也發(fā)生變化。 根據(jù)公式（7）， 潤滑狀態(tài)下的摩擦力與油膜厚度和黏度相關(guān)， 因此接觸面形貌變化會引起接觸面摩擦行為的改變。 同時從油膜厚度的干涉圖中可以看出， 粗糙度形貌各向異性也會影響油膜厚度的分布， 因此可以推斷接觸面摩擦與粗糙度形貌各向異性也相關(guān)。 由于接觸面粗糙度為周期性函數(shù)， 相對摩擦因數(shù)也具有周期性， 如公式（8） 顯示相對摩擦因數(shù)是時間依賴參數(shù)。 因此在以下的研究中， 使用相對摩擦因數(shù)在一個完整周期內(nèi)的平均值作為其在指定工況下的數(shù)值， 即從圖3 中可以看出， 盡管載荷、 速度、初始粗糙度幅值以及波長都一致， 但是對于不同的r， 相對摩擦因數(shù)變化， 即： 當(dāng)r＝1 時μR＝1.790； 當(dāng)r＝2 時μR＝1.566； 當(dāng)r＝0.5 時μR＝1.484。 該結(jié)果直接表明粗糙度表面各向異性特性影響潤滑情況下接觸面的摩擦情況。 下文將會進(jìn)一步探討粗糙度表面的各向異性程度對接觸面摩擦因數(shù)的影響情況。

圖3 3 種粗糙度各向異性情況下的粗糙度形貌、壓力分布、 油膜厚度分布（左圖） 與相對摩擦因數(shù)隨時間的變化情況（右圖）Fig.3 Surface roughness topographies， pressure distributions， film thickness distributions （left） and relative friction coefficient as a function of the time （right） for three different roughness anisotropy conditions： （a） isotropic case；（b） longitudinal case； （c） transverse case

2.3 縱向和橫向表面粗糙度的影響

圖3 表明相同工況下， 不同的表面粗糙度各向異性參數(shù)r造成不同的相對摩擦因數(shù)， 因此， 文中研究了表面粗糙度各向異性對接觸面摩擦因數(shù)的影響， 并將該影響量化。 為了簡便起見， 以下的研究采用相同的載荷、 速度、 粗糙度y方向的波長等工況參數(shù)， 即M＝1 000，L＝10 以及λy／ah＝0.5。 圖4 所示為不同縱向表面粗糙度（參數(shù)r ＞1） 下相對摩擦因數(shù)隨“λ比”， 即Hc／Ai的變化情況。

從圖4 可以看出， 各向異性表面粗糙度情況下的摩擦因數(shù)曲線與各向同性表面粗糙度情況下的曲線具有相同的變化趨勢， 即隨著Hc／Ai的增大， 摩擦因數(shù)逐漸減小。 這是因為隨著Hc／Ai的逐漸增大， 接觸面粗糙度波動逐漸減小， 接觸面逐漸趨于“光滑”， 因此相對摩擦因數(shù)趨于1。 對于每個粗糙度各向異性參數(shù)r， 對應(yīng)著一條光滑的摩擦因數(shù)曲線， 所有的曲線均分布在各向同性r ＝1 曲線的左側(cè)。 從圖4 放大的右圖中可以看出，r＝16 與r＝32 的曲線基本重合，這就說明曲線的左移是有限的。 這些情況與AR（Amplitude Reduction） 理論[24］里的情況相似， 因此如果找到一個和AR 理論相似的函數(shù)將所有的接觸面表面粗糙度各向異性對摩擦因數(shù)曲線的影響統(tǒng)一， 即可以得到一條統(tǒng)一的相對摩擦因數(shù)曲線。 利用縱向表面粗糙度情況下摩擦因數(shù)曲線的右移， 得到了公式（10）， 從而使得摩擦曲線統(tǒng)一為如圖5 所示。

圖4 不同粗糙度各向異性情況（1≤r≤32） 下相對摩擦因數(shù)μR 隨Hc／Ai 的變化情況（左圖）， Hc／Ai∈[2.3， 2.7］ 時的相對摩擦因數(shù)放大情況（右圖）Fig.4 Relative friction coefficient as a function of Hc／Ai for different r （1≤r≤32） values （left），zoom from 2.3 to 2.7 （right）

圖5 縱向表面粗糙度情況下的統(tǒng)一相對摩擦因數(shù)曲線：相對摩擦因數(shù)μR 隨ff（r）·（Hc／Ai） 的變化情況Fig.5 The unified friction curve for longitudinal surface roughness cases： μR as a function of ff（r）·（Hc／Ai）

采用相同的載荷、 速度、 表面粗糙度x方向的波長等工況參數(shù)， 即M＝1 000，L＝10 以及λx／ah＝0.5來研究橫向表面粗糙度對摩擦因數(shù)的影響。 圖6 所示為不同橫向表面粗糙度（參數(shù)r ＜1） 下， 相對摩擦因數(shù)隨“λ比”， 即Hc／Ai的變化情況。 與縱向表面粗糙度情況相比， 橫向表面粗糙度情況較復(fù)雜。 從圖6 中可以看出， 摩擦因數(shù)曲線分布在各向同性r ＝1 曲線的兩側(cè)， 當(dāng)參數(shù)r從0 變化到0.33 時， 相對摩擦因數(shù)逐漸減?。ㄌ摼€）； 當(dāng)參數(shù)r從0.4 變化到1 時，相對摩擦因數(shù)逐漸增大（實線）。

圖6 不同粗糙度各向異性情況（0≤r≤1） 下相對摩擦因數(shù)μR 隨Hc／Ai 的變化情況（上圖）， Hc／Ai∈[2.4， 2.6］時的相對摩擦因數(shù)放大情況（下圖）Fig.6 μR as a function of Hc／Ai for different r （0≤r≤1）values （upper）， zoom from 2.4 to 2.6 （lower）

與縱向表面粗糙度情況相似， 提出一個描述表面粗糙度各向異性程度的放縮函數(shù)， 可以將所有橫向表面粗糙度情況下相對摩擦因數(shù)曲線統(tǒng)一成單一曲線（如圖7 所示）， 該函數(shù)為

圖7 橫向表面粗糙度情況下的統(tǒng)一相對摩擦因數(shù)曲線：相對摩擦因數(shù)μR 隨ff（r）·（Hc／Ai） 的變化情況Fig.7 The unified friction curve for transverse surface roughness cases： μR as a function of ff（r）·（Hc／Ai）

前期研究工作[31］表明， 相對摩擦因數(shù)與工況參數(shù)相關(guān)， 針對各向同性接觸面表面粗糙度情況， 得出了一個放縮系數(shù)θ2， 在較大的工況范圍內(nèi)（500≤M≤2 000， 5≤L≤15 以及0.25≤λ／ah≤1.0） 將所有的摩擦因數(shù)曲線組織成了一條單一曲線。 進(jìn)一步將各向同性與各向異性粗糙度情況結(jié)合， 可得到如公式（12） 所示的統(tǒng)一方程， 其中粗糙度各向異性特性可使用分段函數(shù)（見式（13） ） （如圖8 （a） 所示）來描述。 最后在較大工況范圍內(nèi)將所有的相對摩擦因數(shù)統(tǒng)一成一條單一曲線， 得到非協(xié)調(diào)性表面摩擦副潤滑粗糙接觸情況下的Stribeck 曲線， 如圖8 （b）所示。

圖8 描述粗糙度各向異性特性的函數(shù)ff（r） 隨r 的變化情況（a）， 高負(fù)載非協(xié)調(diào)性接觸下的Stribeck 曲線（b）Fig.8 ff（r） as a function of r（a）， the Stribeck curve for the highly loaded non-conformal contacts （b）

2.4 純縱向表面粗糙度的影響

研究發(fā)現(xiàn)[34］， 瞬態(tài)純縱向表面粗糙度情況下的潤滑接觸情況與穩(wěn)態(tài)純縱向表面粗糙度情況下的結(jié)果相同。 當(dāng)工況為M＝1 000，L＝10，Hc／Ai＝2，λy／ah＝0.5 以及λx＝∞時瞬態(tài)純縱向表面粗糙度情況下， 相對摩擦因數(shù)與穩(wěn)態(tài)結(jié)果的對比如圖9 所示。 正如文獻(xiàn)[34］所述， 除去磨合階段， 瞬態(tài)的結(jié)果和穩(wěn)態(tài)的結(jié)果確實相同。 上述研究發(fā)現(xiàn)， 摩擦因數(shù)的變化主要是由粗糙度（公式（4） ） 的變化引起的， 當(dāng)λx＝∞時，公式（4） 中橫向粗糙度項消失， 表面粗糙度項將變成非時間依賴項， 即此時雷諾方程中的瞬態(tài)項關(guān)于油膜厚度的一階偏導(dǎo)也消失， 即瞬態(tài)雷諾方程退化成穩(wěn)態(tài)雷諾方程。

圖9 瞬態(tài)純縱向表面粗糙度情況下相對摩擦因數(shù)與穩(wěn)態(tài)結(jié)果的對比情況Fig.9 Comparison between the transient relative friction coefficient and that of the stationary case

表4 所示為縱向各向異性粗糙度與純縱向各向異性粗糙度對相對摩擦因數(shù)的影響。 數(shù)據(jù)顯示， 隨著表面粗糙度各向異性參數(shù)r從1 增長到32， 相對摩擦因數(shù)逐漸減小， 然而純縱向各向異性粗糙度情況具有最高的摩擦因數(shù)， 情況與2.3 小節(jié)的結(jié)果不同。 因此有必要將該情況單獨分析。

表4 相對摩擦因數(shù)隨縱向表面粗糙度各向異性程度的變化Table 4 Relative friction coefficient versus different surface roughness orientations

圖10 所示為純縱向表面粗糙度下， 相對摩擦因數(shù)的變化情況（工況范圍：M∈[500， 2 000］，L∈[5， 15］ 與λy／ah∈[0.25， 0.8］ ）。 從圖10 （a）中可以看出， 對于每組工況， 都存在一條光滑的相對摩擦因數(shù)曲線。 由于文獻(xiàn)[31］提出的參數(shù)θ2已經(jīng)不能統(tǒng)一純縱向表面粗糙度情況下的摩擦曲線， 提出了一個新的單獨參數(shù)（表達(dá)式見公式（14） ） 將所有純縱向表面粗糙度情況下的摩擦因數(shù)曲線簇整合成一條單一曲線， 如圖10 （b） 所示。

圖10 純縱向表面粗糙度情況下的相對摩擦因數(shù)Fig.10 Relative friction coefficient for the purely longitudinal wavy case： （a） μR as a function of the original parameter Hc／Ai； （b） μR as a function of the new parameter

3 結(jié)論

（1） 在非協(xié)調(diào)性表面接觸情況下， 潤滑接觸面粗糙度的各向異性對粗糙度形變量影響顯著， 使得接觸面壓力分布、 油膜厚度以及潤滑油黏度發(fā)生變化，進(jìn)而導(dǎo)致摩擦因數(shù)發(fā)生變化。

（2） 接觸面粗糙度各向異性的程度影響潤滑狀態(tài)的摩擦因數(shù)， 該影響可以通過一個函數(shù)ff（r） 量化， 表明全膜潤滑到混合潤滑的過渡不僅與載荷、 速度等工況參數(shù)相關(guān)， 還與粗糙度各向異性相關(guān)。

（3） 純縱向粗糙度情況下， 瞬態(tài)效應(yīng)消失使得潤滑模型轉(zhuǎn)變?yōu)榉€(wěn)態(tài)模型， 參數(shù)θ2不能描繪純縱向表面粗糙度情況， 擬合得到一個單獨的放縮系數(shù)實現(xiàn)純縱向粗糙度情況下摩擦曲線簇的統(tǒng)一。