王蘊馨 馬金奎 陳淑江 路長厚 李 佳 劉志穎

(山東大學機械工程學院,高效潔凈機械制造教育部重點實驗室 山東濟南 250061)

目前滑動軸承動力特性研究常用的計算方法有差分法和偏導數(shù)法。為避免兩次差分近似計算互相疊加的計算誤差[3]，本文作者以圓軸承為對象，采用有限差分法求解雷諾方程，在計算油膜壓力過程中，討論了雷諾邊界條件的實現(xiàn)方法；同時根據(jù)偏導數(shù)法計算擾動壓力并給出了其分布，從而得到滑動軸承的動特性系數(shù)。

在研究滑動軸承流體潤滑時，油膜厚度和動特性系數(shù)的表示通常有2種坐標系，即極坐標系[4]和直角坐標系[5]，本文作者分別建立了2種坐標系下的數(shù)值分析模型以及二者之間的轉(zhuǎn)換矩陣，經(jīng)驗證可獲得相近的結(jié)果。文中的研究完善了經(jīng)典軸承理論的坐標轉(zhuǎn)換方法，對滑動軸承的研究有一定的指導作用。

1 油膜壓力分布

1.1 油膜厚度與雷諾方程的量綱一化形式

圖1為油膜厚度計算簡圖，Ob為軸承中心，Oj為軸的中心；點A為沿y軸正方向逆時針轉(zhuǎn)動θ后的計算點;φ為偏位角，φ=θ-φ；e為偏心距；ω為軸的轉(zhuǎn)速。設c為半徑間隙，可得計算點A的膜厚表達式為

(1)

圖1 油膜厚度計算簡圖Fig 1 A sketch map of oil film thickness calculation

根據(jù)滑動軸承的流體潤滑理論，將雷諾方程量綱一化，得

(2)

1.2 雷諾方程的求解及邊界條件

采用有限差分法求解雷諾方程[6-7]，將油膜劃分為m×n個網(wǎng)格，沿φ方向劃分m等分并用i編號，沿λ方向劃分n等分并用j編號，節(jié)點位置用(i,j)二維編號表示。利用差商取代雷諾方程中的一階偏導數(shù)，整理可得

(3)

其中：

Di,j=Ai,j+Bi,j+Ci,j+Ci,j

Ei,j=3(Hi+1,j-Hi,j)Δφ

計算油膜壓力時，邊界條件的合理采用是影響滑動軸承潤滑性能計算精度的重要因素。文中引入雷諾邊界條件[8-10]：在油膜破裂邊界上，油膜的壓力和壓力的一階導數(shù)均為0。即：

P=0，?P/?φ=0，φ∈[0,2π]

如圖2所示，實現(xiàn)雷諾邊界條件的具體方法是：計算油膜壓力時，由起始點向終止點方向逐點迭代。先假設內(nèi)部各點壓力為0，運用有限差分法依次計算各點的壓力值，如果算出某點壓力值為負數(shù)，則判定此點出現(xiàn)油膜破裂并取為0，這樣就得到了油膜壓力的第1次近似值。利用該近似值再次從第一個點重新計算壓力值，如此反復迭代直至滿足收斂誤差?？梢钥闯觯湍ら_始破裂的位置隨著迭代次數(shù)的增加逐漸向下游推移，則整個壓力曲線隨之增高，于是破裂邊界和壓力分布會逐漸逼近實際破裂邊界和雷諾邊界條件的壓力分布。

圖2 漸變的壓力曲線下游Fig 2 Downstream of gradual pressure curves

從文獻[11-12]中取一組軸承參數(shù)作示例，如表1所示。所得量綱一油膜壓力分布如圖3所示。

表1 軸承參數(shù)Table 1 Parameters of a bearing

圖3 量綱一油膜壓力分布Fig 3 Distribution of dimensionless oil film pressure

可以看出，滑動軸承的量綱一油膜壓力的三維分布近似拋物面分布。在0≤φ≤π的區(qū)域，量綱一油膜壓力在某一段逐漸增大到最大壓力值，之后急劇下降，在φ>π的某一區(qū)域降為0。這與經(jīng)典滑動軸承潤滑理論一致[13-14]。

2 偏導數(shù)法求解滑動軸承動力特性系數(shù)

油膜具有非線性力學特性，其剛度和阻尼會影響系統(tǒng)的運動穩(wěn)定性[15]。油膜剛度和阻尼系數(shù)的2種坐標系下的求解模型如圖4所示。

圖4 滑動軸承動特性系數(shù)計算模型Fig 4 A calculation model of dynamic characteristic coefficients of sliding bearing

圖4中，Ob為軸承中心，O0為靜平衡位置,Oj為軸的瞬時中心，極坐標系為εObφ;直角坐標系為xOby。軸心做微小運動時的量綱一油膜力為

(4)

(5)

2.1 極坐標系求解模型

(6)

定義油膜剛度系數(shù)為單位位移所引起的油膜力增量，即

定義油膜阻尼系數(shù)為單位速度所引起的油膜力增量，即

其中，下標0表示在靜平衡位置處求導；各系數(shù)的第1個下標代表力的方向，第2個下標代表位移或速度的方向(下同)。

(7)

將擾動壓力在完整油膜區(qū)內(nèi)積分，可得到極坐標系中的量綱一剛度系數(shù)和阻尼系數(shù)，即

(8)

(9)

計算時，擾動壓力的邊界條件是：在完整油膜區(qū)的全部周邊上均等于0。根據(jù)表1數(shù)據(jù)及上述分析，可得到滑動軸承量綱一擾動壓力周向分布，如圖5所示。

圖5 量綱一擾動壓力周向分布Fig 5 Circumferential distribution of dimensionless disturbed pressure (a)circumferential distribution of Pε;(b)circumferential distribution of Pφ;(c)circumferential distribution of

2.2 直角坐標系求解模型

同理，將式(5)在靜平衡位置展開為Taylor級數(shù)，僅保留一階微量，可得到

(10)

同樣地，可以定義

(11)

擾動壓力的邊界條件不變，積分后可得到直角坐標系中的量綱一剛度系數(shù)和阻尼系數(shù)，即

(12)

(13)

3 坐標轉(zhuǎn)換矩陣

3.1 坐標轉(zhuǎn)換矩陣推導

(14)

(15)

(16)

已知油膜力可表示為

因此

根據(jù)系數(shù)相等可得2種數(shù)值分析模型的轉(zhuǎn)換方法為

(17)

或

(18)

經(jīng)典軸承理論[1]中的轉(zhuǎn)換矩陣如式(19)所示，由于結(jié)論有所差異，接下來給出驗證分析。

(19)

3.2 驗證分析

首先確定邊界條件，求解壓力分布，然后根據(jù)偏導數(shù)法得到的擾動壓力雷諾方程的量綱一化形式計算各擾動壓力，再經(jīng)積分求出2種計算模型的量綱一動特性系數(shù)，最后通過坐標轉(zhuǎn)換矩陣換算后的結(jié)果，如表2所示。

表2 量綱一動特性系數(shù)(l/d=0.5,φ=30°)Table 2 Dimensionless dynamic characteristic coefficients(l/d=0.5,φ=30°)

2種計算模型得到的量綱一動特性系數(shù)近似相等，最大誤差為6.05%，其中，交叉阻尼系數(shù)近似相等，這有力地論證了上述坐標轉(zhuǎn)換方法的正確性。

4 結(jié)論

(1)分別建立極坐標系和直角坐標系下滑動軸承動特性系數(shù)的計算模型，2種模型均可求解出滑動軸承油膜壓力和擾動壓力分布及其動特性系數(shù)，同時給出了二者之間的坐標轉(zhuǎn)換矩陣，并進行了驗證，說明了方法的正確性。

(2)通過對壓力曲線下游邊界變化情況的分析，結(jié)合迭代算法，研究了雷諾邊界條件的實現(xiàn)方法：先假設內(nèi)部各點壓力為0，運用有限差分法和負壓置零的條件逐點計算出各點的第1次壓力近似值，利用該近似值再次從第一個點重新計算，如此反復迭代直至滿足收斂誤差，從而得到滿足條件的破裂邊界和壓力分布。