金 微,郭 峰,荊兆剛

(青島理工大學(xué) 機(jī)械與汽車工程學(xué)院,山東 青島 266520)

橢圓接觸彈流潤滑廣泛存在于機(jī)械零部件中，如各類滾動軸承和齒輪.在上世紀(jì)70年代，Hamrock和Dowson發(fā)表了關(guān)于橢圓接觸等溫彈流潤滑的研究成果[1-3]，其中包括著名的Hamrock-Dowson膜厚公式(簡稱H-D公式).80年代，朱東和溫詩鑄[4]研究了橢圓度及乏油對潤滑的影響.Chittenden和Dowson等[5-6]研究了卷吸速度與接觸橢圓的長軸重合以及卷吸速度與接觸橢圓的短軸成一定夾角的點(diǎn)接觸問題.在此之后，國內(nèi)外眾多研究人員研究了卷吸速度[7-9]、熱效應(yīng)[10-13]和乏油[14-15]等對彈流潤滑的影響.橢圓比是橢圓接觸彈流潤滑中的1個重要的幾何參數(shù)，在工程實(shí)際中有重要意義，如發(fā)動機(jī)凸輪-挺桿和鼓形軸承滾子的潤滑等.多數(shù)彈流文獻(xiàn)研究的重點(diǎn)大都放在載荷、速度和材料參數(shù)對潤滑狀態(tài)的影響上，對橢圓比的研究相對較少，最具影響的是H-D公式關(guān)于點(diǎn)接觸彈流的膜厚與橢圓比關(guān)系的表達(dá)及Chittenden等的后續(xù)工作.

在潤滑工程設(shè)計中，最常用文獻(xiàn)[2]中的H-D公式估算彈流潤滑膜厚.文中結(jié)果表明無量綱載荷、速度和材料參數(shù)不變時，橢圓比大于6之后，中心膜厚不變.將H-D公式中無量綱載荷用接觸區(qū)中心赫茲壓力表示時，發(fā)現(xiàn)當(dāng)維持接觸區(qū)中心赫茲壓力為常數(shù)，增加橢圓比ke超過一定值時，得到的中心膜厚不升反降，這不符合物理規(guī)律，引起H-D中心膜厚公式的預(yù)測誤差會增大.本文作者應(yīng)用有限單元法(FEM)，借助Comsol Multiphysics軟件求解等溫橢圓接觸彈流問題.研究了橢圓比對中心膜厚的影響，并與H-D中心膜厚公式進(jìn)行對比.構(gòu)造出表征橢圓接觸的綜合幾何參數(shù)，對H-D中心膜厚公式進(jìn)行了修正，使得膜厚的計算更加合理.

1 橢圓接觸彈流潤滑控制方程

穩(wěn)態(tài)等溫橢圓接觸彈流潤滑模型如圖1所示，固體1和2在接觸區(qū)表面沿x方向的速度分別為u1和u2，卷吸速度為ue.令Rx≤ Ry，即令卷吸速度方向與接觸橢圓短軸重合.

Fig.1 Geometrical description of elliptical contact EHL 圖1 橢圓接觸彈流潤滑的幾何描述

無量綱雷諾方程如下，

上述未列參數(shù)見符號表.

1.1 膜厚方程

膜厚方程有三部分組成，包括剛體位移H0，幾何間隙G(X,Y)和接觸面的彈性變形W(X,Y).

對橢圓形滾子，無量綱的幾何間隙為

兩固體間的彈性變形W(X,Y)應(yīng)用三維彈性方程求解.為了簡化模型，應(yīng)用等效的彈性模型計算1個接觸面的彈性變形代替兩個接觸面彈性變形的總和.這樣，固體材料的特性可簡化為

因此，三維彈性方程無量綱形式改寫成式（5）.

其中，λs和μs均為拉梅參數(shù)：

通過上述方程可求解出彈性變形W(X,Y)，更多公式推導(dǎo)細(xì)節(jié)可參考Habchi等[16].

確定方程合適的求解區(qū)域非常關(guān)鍵.橢圓接觸彈流潤滑問題具有對稱性，因此選取半域計算，可減少計算量和縮短計算時間.如圖2所示，計算域?yàn)?Ω，固體立方結(jié)構(gòu)的尺寸至少要大于60倍的赫茲接觸半徑才能滿足半無限空間求解的彈性形變；二維雷諾方程的計算域 Ωc,取-4.5≤X≤1.5和 -3≤Y≤3能夠滿足邊界處壓力為零；?Ωs為對稱面，在XZ平面內(nèi).

因此，上述方程需要滿足的邊界條件為

Fig.2 Computation domain of the elliptical contact EHL圖2 EHL橢圓接觸問題的計算域

1.2 載荷平衡方程

計算域內(nèi)流體動壓產(chǎn)生的壓力總和用于平衡外載荷，考慮到對稱性建立的半域模型，在接觸區(qū)Ωc內(nèi)壓力積分為 π/3，而不是2π/3.

2 數(shù)值計算

同時求解雷諾方程(1)、膜厚方程(2)、線彈性方程(5)和力平衡方程式(6)，可得到接觸區(qū)內(nèi)壓力分布和油膜外形.采用有限元方法求解這些非線性微分方程.在高載荷條件下，雷諾方程是不穩(wěn)定的，會出現(xiàn)振蕩.采用伽遼金最小二乘有限元法(GLS)和各向同性擴(kuò)散技術(shù)對雷諾方程(4)進(jìn)行了修正.因此，修正后的離散的穩(wěn)態(tài)雷諾方程弱表達(dá)式為

其中：前兩項(xiàng)是雷諾方程的伽遼金形式，第三項(xiàng)是懲罰項(xiàng)，使接觸區(qū)內(nèi)負(fù)壓為零，第四項(xiàng)是GLS項(xiàng)，最后1項(xiàng)是ID項(xiàng)，兩者用于消除偽震蕩現(xiàn)象，并且無殘差.Wp是壓力的權(quán)函數(shù)， Ωce是計算域中Ωc單元的集合.方程式(7)中的參數(shù)表達(dá)式如下：

其中：ρID是各向同性擴(kuò)散系數(shù)，he和Pe分別是特征長度和局部Pelect數(shù).

2.1 全耦合模擬

采用通用的有限元軟件進(jìn)行有限元分析.把所有方程聯(lián)立成非線性方程組，作為1個整體，將所有的因變量P、U、V、W和H0整理到1個未知向量中，線性化后用New-Raphson迭代法求解.流體動壓部分采用拉格朗日五次插值函數(shù)，彈性形變部分采用拉格朗日二次插值函數(shù).在接觸區(qū)Ωc內(nèi)，中心赫茲接觸部分網(wǎng)格的最大尺寸不超過0.1，并在靠近出口區(qū)處加密，因?yàn)樵趶椓鳚櫥瑫r此處會出現(xiàn)二次壓力峰.其他部分的網(wǎng)格尺寸最大尺寸不超過0.5.在整個接觸域Ω內(nèi)選擇比較粗糙的網(wǎng)格就能滿足計算要求.自由度大約為60 419，相對容差在10?3～10?4之間.使用Inter(R) Xeon(R) CPU E5-2630處理器，迭代4次，經(jīng)過20 s即可收斂.如果選擇合適的壓力初值，能減少迭代步數(shù)，縮短計算時間，因此合理地選擇壓力初值尤為重要.

2.2 驗(yàn)證

本文作者應(yīng)用有限單元法計算彈流潤滑的壓力分布和油膜厚度，與多重網(wǎng)格法(MG)的計算結(jié)果相互比較，如圖3所示.計算參數(shù)為ke=1、α=2.19×10?8Pa?1、Rx=0.02 m、η0=0.08 Pa·s、UR=5.0×10?11、W=1.0×10?6和G=5 000.圖中對比了Y=0處的壓力分布和油膜厚度，兩種方法的計算結(jié)果吻合，證明了本文方法的正確性.

Fig.3 Comparison of results between finite element method with Multigrid method圖3 有限單元法和多重網(wǎng)格法計算結(jié)果的比較

3 結(jié)果與討論

3.1 載荷不變時中心膜厚隨橢圓比的變化

圖4所示ke= 2、4、6、10、15時，橢圓接觸油膜厚度和壓力分布等值云圖.在油膜厚度等值云圖中，橢圓比對馬蹄形油膜的形成具有重要影響.隨著橢圓比ke的增加，油膜的馬蹄形特征越不明顯，兩個側(cè)緣越來越小，之間的距離越來越遠(yuǎn)；油膜的平行部分縮短.在壓力分布等值云圖中，隨著橢圓比ke的增加，二次壓力峰先升高到最大值之后開始下降，并向接觸區(qū)中心移動.這表明當(dāng)外載荷一定時，增加ke，承載面積增加，潤滑狀態(tài)由彈流潤滑變化到動壓潤滑.

Fig.4 Pressure and film thickness vs elliptic ratios under constant load (UR =5.0×10?11，W=1.0×10?6)圖4 載荷不變時油膜厚度和壓力隨橢圓比的變化 (UR =5.0×10?11，W=1.0×10?6)

圖5所示為本文數(shù)值計算和H-D公式得到的中心膜厚Hcen隨橢圓比的變化.ke的取值范圍從1(點(diǎn)接觸)到16(中心接觸區(qū)為線接觸狀態(tài))，可以看出，兩者趨勢相差很大.在ke=1～8之間，H-D公式和數(shù)值計算表現(xiàn)出相似的變化趨勢，Hcen隨ke的變化表現(xiàn)為先快后慢，數(shù)值的差別應(yīng)來自于計算方法和網(wǎng)格數(shù)的不同.ke>8以后，H-D公式中的Hcen基本不變；當(dāng)UR=5.0×10?11，W=1.0×10?6時，在本文數(shù)值計算中，Hcen隨ke的增加而逐漸下降.隨橢圓比ke增加承載區(qū)域增加，當(dāng)外載荷不變時，油膜壓力水平降低，導(dǎo)致油膜厚度增加，潤滑狀態(tài)由彈流潤滑轉(zhuǎn)變到動壓潤滑，同時油膜壓力的降低會使固體彈性變形變小，對油膜厚度的貢獻(xiàn)降低，又會使得油膜降低.最終的油膜厚度是以上兩種效應(yīng)共同作用的結(jié)果.需要指出的是，Hamrock和Dowson在H-D中心膜厚公式回歸過程中，使用的橢圓比最大為8.橢圓比ke= 8～16之間的油膜厚度公式?jīng)]有數(shù)值計算點(diǎn)的約束，只是回歸公式的外插值.

3.2 接觸中心Hertz壓力不變時中心膜厚隨橢圓比的變化

為了進(jìn)一步檢驗(yàn)H-D公式中Hcen與ke的關(guān)系，保持接觸區(qū)中心赫茲接觸壓力不變，就ke= 1～16時對應(yīng)的Hcen值進(jìn)行了計算，結(jié)果見圖6，其中表示由H-D公式計算出的中心膜厚，表示有限單元法數(shù)值計算的結(jié)果.以ph=0.4 GPa，UR= 5.0×10?12為例，由圖6可知，隨著ke的增加，緩慢增加，在ke>10以后趨近于定值.在ke=1～5時隨ke的增加而增加，而ke>5以后呈現(xiàn)下降的趨勢，ke=16時的計算值(1.549)與ke=5時的計算值(1.654)相比，下降了6.3%，這與物理現(xiàn)象是相悖的.當(dāng)接觸中心壓力不變時，隨著ke增加，端泄對中心膜厚影響減小，中心膜厚不應(yīng)出現(xiàn)減小的趨勢.

定義幾何參數(shù)θ

最大接觸Hertz壓力可改寫為

Fig.5 Central film thickness vs ellipticity ratio under constant load by the method of H-D formula and FEM圖5 由H-D公式(a)和FEM法(b)計算出恒定載荷時中心膜厚與橢圓比的關(guān)系

Fig.6 Central film thickness vs.ellipticity ratio under constant Hertz pressure calculated by H-D formula and FEM圖6 由H-D公式(a)和FEM法(b)計算出恒定赫茲壓力時中心膜厚與橢圓比的關(guān)系

引用參數(shù)θ對H-D公式等效變換，θ僅與橢圓比有關(guān)，而與Rx、Ry無關(guān).由式(10)得

H-D中心膜厚公式:

因E′與E*相差2倍，所以H-D公式中的為

代入H-D中心膜厚公式(12)中，得

其中：esl稱為中心膜厚端泄因子

從理論上講，esl應(yīng)為ke的增函數(shù)，而實(shí)際上，由表達(dá)式(15)求出當(dāng)ke≈4.6時，esl取得最大值.圖7為按式(15)繪出的esl與ke的關(guān)系圖.可以看出在ke≈4.6之后，esl為ke的減函數(shù).因此，可經(jīng)過大量的數(shù)值計算，構(gòu)造出以θ為基本參數(shù)的端泄因子表達(dá)式，對現(xiàn)有的端泄因子進(jìn)行修正.

Fig.7 Variation of leakage factor esl as a function of ellipticity ratio ke圖7 端泄因子esl與橢圓比ke的關(guān)系

Fig.8 Variation of modified leakage factor eslm as a function of ellipticity ratio ke by numerical fitting圖8 數(shù)值擬合出修正端泄因子eslm與橢圓比ke的關(guān)系曲線

圖8分別給出了不同的載荷和速度時，根據(jù)數(shù)值計算得到的修正端泄因子eslm與ke的關(guān)系.其中散點(diǎn)是計算值，實(shí)線是擬合曲線.根據(jù)離散點(diǎn)擬合出如下修正的中心膜厚端泄因子:

同時在圖8上方繪制出了與ke對應(yīng)的θ值.可以看出，θ與esl成反比的關(guān)系.但是，由H-D中心膜厚公式推導(dǎo)出的端泄因子式(15)中，θ與esl是正比的關(guān)系.由此可見，修正后的端泄因子式(16)更加合理.

3.3 橢圓接觸彈流潤滑中心膜厚修正公式

根據(jù)以上數(shù)值分析的數(shù)據(jù)，得到以下橢圓接觸彈流潤滑中心膜厚的修正公式.

與經(jīng)典的Hamrock-Dowson公式相比，式(17)的中心厚度對端泄因子進(jìn)行了修正.對修正后的公式精度進(jìn)行了檢驗(yàn)，結(jié)果列于表1中.

由表1可見，通過上述修正公式得到結(jié)果與數(shù)值計算結(jié)果具有較好吻合性.

表1 中心膜厚數(shù)值計算結(jié)果及修正公式計算結(jié)果比較Table 1 Comparison between the central film thickness of the modified formula and numerical results

4 結(jié)論

在等溫穩(wěn)態(tài)橢圓接觸彈流潤滑條件下，研究了油膜壓力和油膜厚度隨橢圓比的變化規(guī)律.對比中心膜厚的數(shù)值計算值和Hamrock-Dowson公式計算值，得到如下結(jié)論：

a.在無量綱載荷一定時，隨著橢圓比的增加，中心膜厚先增加后減小，潤滑狀態(tài)從彈流潤滑轉(zhuǎn)變到動壓潤滑.

b.中心最大Hertz壓力一定時，隨著橢圓比的增加，Hamrock-Dowson公式計算出的中心膜厚出現(xiàn)了與物理相悖的現(xiàn)象，體現(xiàn)在ke≈4.6之后，端泄因子esl為ke的減函數(shù).

c.依據(jù)數(shù)值計算結(jié)果擬合出關(guān)于構(gòu)造參數(shù)θ的修正端泄因子和修正Hamrock-Dowson公式，并檢驗(yàn)了公式精度，兩者吻合度良好.

符號說明