高黏度油潤滑彈塑性接觸問題數值仿真*

張生光 劉 剛 胡文穎

(1.中國航空發動機研究院 北京 101304；2.中國黃金集團(上海)貿易有限公司 上海 200126)

隨著高速機床、航空航天、鐵路運輸等工程技術的迅速發展，對高速、重載下的機械裝備提出了更加苛刻的高精度、長壽命等要求。滾動軸承作為提供回轉支撐的關鍵零部件，是機械系統中不可或缺的組成部分。在軸承制造和使用過程中，由于重載或振動等工況，滾動體與滾道接觸時可能發生塑性變形，從而造成主軸精度下降，喪失精度保持性。而軸承需要工作在潤滑條件下，其發生塑性變形時是典型的彈塑性潤滑工況，因此，有必要對彈塑性潤滑問題開展研究。

目前，針對彈流潤滑問題，已經開展大量研究，并建立了有效的仿真計算模型，取得了豐富的成果[1-5]。在彈流潤滑模型基礎上，牛榮軍和黃平[6]考慮材料塑性變形影響，研究了線接觸粗糙表面彈塑性潤滑問題；采用余弦函數表征粗糙表面，分析了粗糙度對彈塑性潤滑的影響。隨后，也有一些研究者對彈塑性潤滑問題進行了研究，他們多采用基于各向同性強化準則下的J2流動理論求解塑性應變，應力-應變本構關系多采用線性強化準則或者Swift強化準則。REN等[7]結合統一雷諾方程和半解析法求解了點接觸彈塑性潤滑問題，通過對比彈塑性潤滑與彈流潤滑結果，發現考慮塑性變形后最大壓力明顯降低，并且壓力范圍明顯增大，同時潤滑油膜厚度也有所降低。REN等[8]還研究了單個微凸峰和凹坑通過接觸區時對彈塑性潤滑的影響，結果表明存在單個微凸峰和凹坑時，表層下應力集中較大，即使載荷較小時也較容易發生塑性變形。HE等[9-10]發展了REN等[7]的模型，研究了不同粗糙表面點接觸情況以及有限長與無限長線接觸情況，同時還考慮了塑性應變引起的材料硬化對彈塑性潤滑的影響。在此基礎上，HE等[11]研究了鋼球單次和多次滾過表面后，塑性變形及其累計對接觸行為的影響，給出了等效塑性應變的分布情況，為理解重載軸承滾道變形機制提供了有力的理論支撐。

上述研究成功建立了彈塑性潤滑問題的數學模型并提出了合適的求解方法，這對于低黏度潤滑油的求解是可行的。但如果采用高黏度潤滑油，上述求解方法尚存在一定問題。由于塑性變形計算采用分步加載的方式完成，如果采用高黏度潤滑油并分步加載，計算結果的接觸區壓力曲線會存在較高的第二壓力峰，進而導致較大塑性應變的產生，與實際載荷產生的塑性應變存在較大偏差。因此，為解決高黏度潤滑油工況下的彈塑性潤滑仿真問題，本文作者采用半解析法建立彈塑性潤滑接觸模型，塑性變形采用特征應變方法進行求解，并且通過逐步增加潤滑油初始黏度的方式來進行仿真計算，最終實現高黏度潤滑油工況下的彈塑性潤滑有效仿真分析。

1 彈塑性潤滑建模與求解

潤滑條件下，接觸壓力通過雷諾方程求解：

(1)

式中：x、y表示坐標方向；h為兩表面間隙；p為接觸壓力；u為卷吸速度；t為時間；η和ρ分別為潤滑油黏度與密度。

求得接觸壓力后，可進一步求得次表層應力分布。對于彈塑性接觸，表層下應力由接觸壓力產生的施加應力和塑性應變產生的殘余應力兩部分組成，即：

σ=σ0+σp

(2)

式中：σp為塑性應變εp引起的殘余應力，可通過特征應變的方法進行計算[12-13]；σ0為接觸壓力產生的施加應力：

(3)

文中研究，彈塑性模型本構方程遵循小變形假設，應力-應變本構關系采用圖1所示的線性強化準則；而對于理想彈塑體，切向模量ET=0。基于各向同性強化準則下的J2流動理論，采用徑向返回算法[15]求解塑性應變εp。

圖1 應力-應變本構關系

對于點接觸，兩接觸體間隙方程為

h(x,y)=h0+fg(x,y)+δ(x,y)+ve(x,y)+vp(x,y)

(4)

式中：h0為兩接觸體未考慮變形時的相互接近位移；fg為兩接觸體間未考慮變形時的幾何間隙；δ為表面粗糙度；ve為表面彈性變形；vp表示塑性應變εp引起的塑性變形。

當材料為純彈性體時，塑性變形vp=0，潤滑模型由彈塑性流體潤滑變為彈性流體動壓潤滑。

潤滑油黏度按Barus黏-壓關系計算：

η=η0exp(αp)

(5)

式中：η0表示環境黏度；α為黏-壓系數。

密度通過Dowson-Higginson密-壓關系計算：

ρ=ρ0[1+Cap/(1+Cbp)]

(6)

式中：ρ0為初始壓力和初始溫度下的密度；Ca和Cb為系數，Ca取6×10-10m2/N，Cb取1.7×10-9m2/N。

載荷平衡方程為

w=?Ωp(x,y)dxdy

(7)

為計算方便，需要將各方程量綱一化處理。將計算域劃分網格并用有限差分法離散量綱一化方程。數值求解時，彈性變形、塑性變形、表層下應力等采用FFT快速算法加速計算，壓力求解采用Gauss-Seidal低松弛迭代逐行掃描法求解，具體差分格式及求解步驟可參見文獻[16]。計算時，塑性變形和加載路徑有關，應分步加載。由于塑性應變產生的塑性變形會影響接觸壓力，接觸壓力產生的應力又會影響塑性應變的求解，因此各參數相互耦合，需迭代求解直至收斂。

2 高黏度潤滑油工況下計算方法分析

當選用的潤滑油黏度較小時，如η0=0.03 Pa·s，載荷選取為100 N，采用半徑20 mm的光滑剛性球與彈性模量206 GPa、泊松比0.3的半無限大光滑理想彈塑體接觸，屈服極限取σs=600 MPa，彈塑性潤滑接觸壓力計算結果如圖2(a)所示。由于塑性變形的存在，最大接觸壓力明顯小于最大赫茲壓力pH，并且接觸區域略大于赫茲接觸區，計算結果比較符合實際情況。但是當潤滑油黏度較大時，如η0=0.25 Pa·s，計算結果如圖2(b)所示。在x>0區域壓力存在波動，由于兩表面均為光滑表面，不存在由于表面粗糙度引起壓力波動的可能，因此，該壓力波動的存在并不合理。而通過觀察圖2(c)所示的各分步加載下的壓力分布可見，波動的位置恰好對應各分步載荷產生的二次壓力峰位置。這是因為潤滑油黏度較大，分步加載時，在各個分步載荷下，二次壓力峰都比較大，會在其正下方產生較大塑性應變以及表面殘余變形。隨著分步加載載荷逐漸增大，接觸面積增加，接觸壓力區域增大，原來二次壓力峰位置的壓力實際有所減小，但是由于上一步加載時該位置存在較大塑性應變，無法消除，因此會在該位置產生較大的表面殘余變形，最終導致接觸壓力的波動。

圖2 不同黏度下彈塑性潤滑接觸壓力

為解決上述問題，文中采用了步進黏度法來計算彈塑性潤滑問題。計算過程如下，分步加載時，令初始黏度η0=η0/n。根據黏度的大小，n可以取不同值，黏度較大時，n取值較大，黏度較小時n取值較小或取1。當載荷加至給定值后，在壓力收斂后逐步增大初始黏度即令η0=η0×i/n，直至η0為所需黏度。求解過程中，為使方程簡化，通常引入量綱一化參數，文中模型位移、壓力以及油膜厚度的量綱一化參數為：X=x/a，Y=y/a，P=p/pH，H=h/a。其中a為赫茲接觸半徑，pH為最大赫茲接觸壓力。彈塑性潤滑問題的具體計算流程如圖3所示。采用步進黏度法計算的壓力分布如圖2(d)所示，接觸區內壓力分布光滑無波動，符合實際情況。

圖3 彈塑性潤滑計算流程

為驗證模型的準確性，與文獻[11]中的彈塑性潤滑計算結果進行對比，對比結果如圖4所示。通過對接觸壓力對比，驗證了模型的準確性。

圖4 模型驗證

3 彈塑性潤滑計算結果及分析

彈塑性潤滑模型采用圖2相同的參數，假設半無限大體為理想彈塑體，接觸載荷取150 N，對應的赫茲接觸半徑與最大赫茲接觸壓力分別為a=216 μm，pH=1 532 MPa；計算域為{-2

當潤滑油黏度較低時，取η0=0.02 Pa·s，不同卷吸速度下沿中心線的接觸壓力、油膜厚度以及y=0平面的Mises應力分布如圖5所示。其中圖5 (a)—(c)示出了不同卷吸速度下的彈流潤滑計算結果，可見卷吸速度較低時，油膜厚度較小，接觸壓力與赫茲接觸相似；隨著卷吸速度的增加，潤滑油膜厚度逐漸增加，彈流特征逐漸明顯；由于彈流潤滑出口處二次壓力峰的影響，對應位置的Mises應力略有增大。圖5(d)—(f)示出了不同卷吸速度下的彈塑性潤滑計算結果，與彈流潤滑計算結果相比，壓力主峰明顯減小，二次壓力峰變化較小；隨著卷吸速度的增加，油膜厚度也逐漸增大，其數值與彈流潤滑的結果接近；從Mises應力結果可見，表層下區域已發生部分屈服，由于材料設置為理想彈塑體，因此其最大Mises應力只能達到屈服極限600 MPa。

圖5 η0=0.02 Pa·s時不同速度下沿中心線壓力、膜厚以及y=0平面Mises應力分布

圖6給出了潤滑油黏度η0=0.25 Pa·s時不同卷吸速度下沿中心線的接觸壓力、油膜厚度以及y=0平面的Mises應力分布。潤滑油黏度較高時，PEHL計算需采用步進黏度法進行處理。圖6(a)—(c)示出了不同卷吸速度下的彈流潤滑計算結果，與低黏度潤滑油相比，油膜厚度明顯增大，接觸壓力的二次壓力峰也更加明顯，表層下Mises應力分布與低黏度潤滑油相似，數值上也相差不大。圖6(d)—(f)示出了不同卷吸速度下的彈塑性潤滑計算結果，可見采用步進黏度法后壓力、膜厚都比較平滑，油膜厚度與低黏度結果相比具有較大提升；壓力主峰相比于彈流潤滑計算結果明顯降低，二次壓力峰相比于彈流潤滑計算結果，不再是隨著速度的增加一直增大，而是增大到一定程度后就逐漸平穩。

圖6 η0=0.25 Pa·s時不同速度下沿中心線壓力、膜厚以及y=0平面Mises應力分布

表1給出了η0=0.02 Pa·s和η0=0.25 Pa·s情況下中心膜厚與最小膜厚的結果。當黏度與速度均較低時，最小膜厚為0，發生接觸；隨著速度升高，中心膜厚與最小膜厚均增加，并且高黏度潤滑油的油膜厚度比低黏度潤滑油高出較多。但是彈流潤滑的中心膜厚與彈塑性潤滑中心膜厚十分接近，幾乎沒有變化；彈流潤滑的最小膜厚比彈塑性潤滑的稍高。

表1 彈流潤滑(EHL)與彈塑性潤滑(PEHL)中心膜厚與最小膜厚對比

不同黏度、卷吸速度下的表面殘余變形及等效塑性應變計算結果如圖7所示。圖中表面殘余變形為量綱一化量，表面殘余變形量S與赫茲接觸半徑a的比值。塑性變形與塑性應變區域主要與載荷大小相關，因此對比圖7(a)、(b)、(c)、(d)可見，在宏觀上塑性應變區域變化較小。隨著卷吸速度的增加，表面殘余變形略有增加，并且由于高速條件下二次壓力峰較明顯，導致其下方對應位置會產生較大的應力集中，從而更容易發生屈服。對于較高的速度和黏度工況，塑性應變區域會向表面及接觸區出口處移動。

圖7 不同黏度、速度下表面殘余變形及等效塑性應變

4 結論

(1)研究高黏度潤滑油彈塑性接觸問題數值仿真方法，利用所提出的步進黏度法計算高黏度潤滑油彈塑性接觸模型，消除了接觸區壓力曲線的波動。

(2)彈塑性潤滑壓力主峰相比于彈流潤滑結果明顯降低，二次壓力峰相對變化較小。彈流潤滑的中心膜厚與彈塑性潤滑中心膜厚接近，最小膜厚稍大。

(3)對于較高的速度和黏度工況，塑性應變區域會向表面及接觸區出口處移動。