基于高階單元的滑動軸承潤滑特性計算

楊國棟， 張文平， 曹貽鵬， 明平劍， 李燎原

(1.哈爾濱工程大學 動力與能源工程學院，黑龍江 哈爾濱 150001；2.中國艦船研究設計中心，湖北 武漢 430064)

滑動軸承在工業生產和生活領域中應用較廣，具有承載力大、穩定性好等特點[1]。正常工作時，軸和軸瓦之間會形成楔形的幾何空間，潤滑油在界面運動的帶動下進入該間隙，進而形成動壓強，起到承受外載荷和減小摩擦力的作用。雷諾方程是求解滑動軸承潤滑特性的基本方程，一維雷諾方程具有解析解，但是這類方程左側缺少了1個方向的泊肅葉流項，會導致計算得到的壓強值偏大，不符合實際情況。二維雷諾方程可以很好地描述有限寬軸承的壓強分布，但是沒有解析解，通常采用數值法求解。有限差分法[2-7]是較常用的數值解法，該算法用差商來代替導數，形式簡單，但只能處理正交化的四邊形網格，不適用于復雜的計算域；有限元法[8-10]從彈性力學發展而來，后來應用于潤滑分析中，求解域的網格劃分更為靈活，節點壓強梯度可以采用單元形函數導數的方式處理，但該算法需要構建較大的剛度矩陣，計算消耗較大；非結構有限體積法兼具有限體積法守恒性的特點和對復雜區域良好適用性的特點，在傳熱學[11]和結構領域[12-13]應用較多，在潤滑領域應用較少[14-15]，并且尚未涉及高階單元，而隨著對計算精度和效率的要求不斷提高，有必要從高階單元的角度出發開展相關研究。

本文使用六節點三角形單元劃分求解域，采用多重網格法求解離散的代數方程組，進而研究軸承的潤滑特性。

1 潤滑理論方程及雷諾方程的離散

矢量形式的雷諾方程為：

(1)

軸承的幾何關系如圖1所示。

圖1 軸承幾何關系Fig.1 Geometric position of journal bearing

設軸心坐標為(x,y)，徑向滑動軸承的液膜厚度為：

h=c+ecos(θ-φ)=

c+ecosφcosθ+esinφsinθ=

c+xcosθ+ysinθ

經過調研分析，對于系統功能的需求主要分為對基礎地理信息的查詢和操作功能、配方施肥決策功能和用戶管理功能3個方面。基礎GIS功能主要包括地圖縮放及漫游、查詢、圖層顯示及控制、測量功能；配方施肥決策功能主要包括養分含量查詢功能、施肥配方決策功能、土壤肥力評價功能等；用戶管理功能主要包括用戶注冊、認證、權限管理等功能。系統功能結構見圖2。

(2)

式中：c為軸承間隙；e為偏心距；φ為偏位角。油膜反力方程、端泄方程和摩擦系數方程參考文獻[1-2]。

平均誤差計算公式為：

(3)

式中：m表示潤滑區域的節點總數；pi表示各個算例中的節點壓強；pi,ref表示參考節點壓強。

圖2為本文采用的2類三角形單元及其節點的示意圖，A～F表示節點，內部空心點為積分點。格點型有限體積法的控制體是圍繞節點形成的，如圖3所示，在控制體內，對式(1)積分可得：

(4)

擴散項可按下式進行離散：

(5)

式中：nc表示節點周圍的單元數；ncni表示第i個單元的節點總數；N代表形函數；nα(α=x,y)為面單位外法線失量n沿α方向的分量；ψ表示圍繞控制體的積分線。

圖2 2種三角形單元示意Fig.2 Two kinds of triangular elements

圖3 控制體和積分點示意Fig.3 The control volume and its integration points

源項離散為：

(6)

則式(4)可轉化為：

(7)

2 形函數的處理

2.1 常應變三角形單元

如圖4所示，1～3為三角形單元的3個節點，A、B、C分別為邊L12、L23、L31的中點，O為單元中心，型函數等信息可參考文獻[17]，積分點坐標可由幾何關系獲得。

圖4 三節點三角形單元的坐標轉換Fig.4 Coordinate transformation of 3-node triangular element

在節點周圍的第i個單元內，形函數導數沿控制體邊界的積分為：

(j=1,2,3;α=x,y)

(8)

式中Ljα為圍繞第j個節點的積分線在α方向上的投影長度。

利用雅克比矩陣，經過推導可得型函數對全局坐標的導數為：

(9)

式中：(x1,y1)、(x2,y2)、(x3,y3)分別為全局坐標下三角形節點1、2、3的節點坐標；S123為全局坐標下三角形單元的面積。

2.2 六節點三角形單元

如圖5所示，1～6為高階三角形單元的6個節點編號，S1～S9為高斯積分點的編號，六節點三角形單元的形函數和積分點坐標見文獻[17]。

在節點周圍的第i個單元內，形函數導數沿控制體邊界的積分為(其余推導同2.1):

(j=1,2,3,4,5,6;α=x,y)

(10)

圖5 六節點三角形單元的坐標轉換Fig.5 Coordinate transformation of 6-node triangular element

3 軸承實例分析

本文的計算程序是在哈爾濱工程大學動力裝置工程技術研究所自主開發的通用輸運方程求解器GTEA軟件的基礎上開發的，采用Fortran90語言編程，在參數為4 GB RAM、Intel Core i5-2400、CPU 3.1 GHz的計算機中運行程序。

3.1 算法準確性驗證

本節采用文獻[2]中滑動軸承的計算結果來驗證算法的準確性，軸承半徑為30 mm，長度為66 mm，間隙為0.03 mm，粘度為0.009 Pa·s，轉速為3 000 r/min。網格劃分如表1所示，其中Case0為驗證算例所用的模型，采用四邊形網格劃分，其結果也將用于下一節中的分析。由圖6可得，Case0的壓強分布與文獻[2]中的結果基本一致；由表2可得Case0的最大壓強、端泄流量、摩擦系數與文獻均較接近，因此，非結構格點型有限體積法的正確性得以驗證，該算法可用于下一步的分析。

表1 網格劃分結果Table 1 Grid meshing results

表2 油膜壓強最大值的對比

圖6 Case0壓強分布準確性驗證Fig.6 Verification of the accuracy of pressure distribution of Case0

3.2 高階單元的結果對比

本節采用表1中的算例研究高階單元的影響，其中，Case1和Case2采用三節點三角形單元，Case3和Case4采用六節點三角形單元，由表1可知，Case1和Case3的節點數一致，Case2和Case4的節點數一致。由于算例Case0結果的正確性已經得到驗證，因此Case0的結果可作為參照。

由圖7可知，Case1～Case4的壓強分布與Case0基本一致，說明采用高階單元時，雖然單元數目較少，但是可以得到較為準確的壓強分布。

由表2可知，各算例的最大壓強均與文獻[2]較為接近，說明在網格數較少時，高階單元也可以得到較準確的峰值壓強；同時，Case3的反力大于Case1，Case4的反力大于Case2，且與參考值更接近；Case3的端泄流量和摩擦系數均小于Case1，與參考值較為接近，且Case4和Case2亦如此，說明節點數一樣時，高階單元計算得到的潤滑結果更準確；由表3可知，Case3的內存和計算耗時均大于Case1，Case4的內存和計算耗時均大于Case2，但是Case3和Case4的計算誤差更小，說明在節點數一致時，雖然高階單元的單元數較少，但是精度更高。

圖7 不同算例壓強分布對比Fig.7 Comparison of pressure distribution of different cases

表3 計算時間、內存及誤差的對比Table 3 Comparison of calculation time, memory and error

3.3 實例計算

以某船用滑動軸承為例，采用六節點三角形單元和非結構格點型有限體積法研究軸承壓強分布隨傾角的變化規律，設定偏位角為0°，粘度為0.01 Pa·s，軸承半徑為0.05 m，長度為0.2 m，間隙為0.04 mm，偏心率為0.6，轉速為300 r/min。

由圖8可知，隨著傾斜角度的增加，軸承壓強峰值逐漸向軸承的端部移動，同時軸承壓強值逐漸增大，這是由于軸頸傾斜時，端部位置的油膜厚度急劇減小造成的。

圖8 傾角對軸承壓強分布的影響Fig.8 Effects of tilt angle on oil film pressure distribution

4 結論

1) 高階單元對軸承壓強分布和峰值壓強的計算影響較小，但是會使油膜反力值更接近參考值。

2) 節點數一致時，高階單元得出的端泄流量和摩擦系數更接近參考值。

3) 采用高階單元時，滑動軸承液膜壓強計算的精度會相應提高，但是內存消耗會相應增大，計算效率下降。

4)高階單元模型可以適用于傾斜狀態下的船用滑動軸承。