高黏度油潤(rùn)滑彈塑性接觸問題數(shù)值仿真*

張生光 劉 剛 胡文穎

(1.中國(guó)航空發(fā)動(dòng)機(jī)研究院 北京 101304；2.中國(guó)黃金集團(tuán)(上海)貿(mào)易有限公司 上海 200126)

隨著高速機(jī)床、航空航天、鐵路運(yùn)輸?shù)裙こ碳夹g(shù)的迅速發(fā)展，對(duì)高速、重載下的機(jī)械裝備提出了更加苛刻的高精度、長(zhǎng)壽命等要求。滾動(dòng)軸承作為提供回轉(zhuǎn)支撐的關(guān)鍵零部件，是機(jī)械系統(tǒng)中不可或缺的組成部分。在軸承制造和使用過程中，由于重載或振動(dòng)等工況，滾動(dòng)體與滾道接觸時(shí)可能發(fā)生塑性變形，從而造成主軸精度下降，喪失精度保持性。而軸承需要工作在潤(rùn)滑條件下，其發(fā)生塑性變形時(shí)是典型的彈塑性潤(rùn)滑工況，因此，有必要對(duì)彈塑性潤(rùn)滑問題開展研究。

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

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

1 彈塑性潤(rùn)滑建模與求解

潤(rùn)滑條件下，接觸壓力通過雷諾方程求解：

(1)

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

求得接觸壓力后，可進(jìn)一步求得次表層應(yīng)力分布。對(duì)于彈塑性接觸，表層下應(yīng)力由接觸壓力產(chǎn)生的施加應(yīng)力和塑性應(yīng)變產(chǎn)生的殘余應(yīng)力兩部分組成，即：

σ=σ0+σp

(2)

式中：σp為塑性應(yīng)變?chǔ)舙引起的殘余應(yīng)力，可通過特征應(yīng)變的方法進(jìn)行計(jì)算[12-13]；σ0為接觸壓力產(chǎn)生的施加應(yīng)力：

(3)

文中研究，彈塑性模型本構(gòu)方程遵循小變形假設(shè)，應(yīng)力-應(yīng)變本構(gòu)關(guān)系采用圖1所示的線性強(qiáng)化準(zhǔn)則；而對(duì)于理想彈塑體，切向模量ET=0?；诟飨蛲詮?qiáng)化準(zhǔn)則下的J2流動(dòng)理論，采用徑向返回算法[15]求解塑性應(yīng)變?chǔ)舙。

圖1 應(yīng)力-應(yīng)變本構(gòu)關(guān)系

對(duì)于點(diǎn)接觸，兩接觸體間隙方程為

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

(4)

式中：h0為兩接觸體未考慮變形時(shí)的相互接近位移；fg為兩接觸體間未考慮變形時(shí)的幾何間隙；δ為表面粗糙度；ve為表面彈性變形；vp表示塑性應(yīng)變?chǔ)舙引起的塑性變形。

當(dāng)材料為純彈性體時(shí)，塑性變形vp=0，潤(rùn)滑模型由彈塑性流體潤(rùn)滑變?yōu)閺椥粤黧w動(dòng)壓潤(rùn)滑。

潤(rùn)滑油黏度按Barus黏-壓關(guān)系計(jì)算：

η=η0exp(αp)

(5)

式中：η0表示環(huán)境黏度；α為黏-壓系數(shù)。

密度通過Dowson-Higginson密-壓關(guān)系計(jì)算：

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

(6)

式中：ρ0為初始?jí)毫统跏紲囟认碌拿芏?；Ca和Cb為系數(shù)，Ca取6×10-10m2/N，Cb取1.7×10-9m2/N。

載荷平衡方程為

w=?Ωp(x,y)dxdy

(7)

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

2 高黏度潤(rùn)滑油工況下計(jì)算方法分析

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

圖2 不同黏度下彈塑性潤(rùn)滑接觸壓力

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

圖3 彈塑性潤(rùn)滑計(jì)算流程

為驗(yàn)證模型的準(zhǔn)確性，與文獻(xiàn)[11]中的彈塑性潤(rùn)滑計(jì)算結(jié)果進(jìn)行對(duì)比，對(duì)比結(jié)果如圖4所示。通過對(duì)接觸壓力對(duì)比，驗(yàn)證了模型的準(zhǔn)確性。

圖4 模型驗(yàn)證

3 彈塑性潤(rùn)滑計(jì)算結(jié)果及分析

彈塑性潤(rùn)滑模型采用圖2相同的參數(shù)，假設(shè)半無限大體為理想彈塑體，接觸載荷取150 N，對(duì)應(yīng)的赫茲接觸半徑與最大赫茲接觸壓力分別為a=216 μm，pH=1 532 MPa；計(jì)算域?yàn)閧-2

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

(3)對(duì)于較高的速度和黏度工況，塑性應(yīng)變區(qū)域會(huì)向表面及接觸區(qū)出口處移動(dòng)。