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

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

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

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

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

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

矢量形式的雷諾方程為：

(1)

軸承的幾何關(guān)系如圖1所示。

圖1 軸承幾何關(guān)系Fig.1 Geometric position of journal bearing

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

h=c+ecos(θ-φ)=

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

c+xcosθ+ysinθ

經(jīng)過調(diào)研分析，對于系統(tǒng)功能的需求主要分為對基礎(chǔ)地理信息的查詢和操作功能、配方施肥決策功能和用戶管理功能3個方面?；A(chǔ)GIS功能主要包括地圖縮放及漫游、查詢、圖層顯示及控制、測量功能；配方施肥決策功能主要包括養(yǎng)分含量查詢功能、施肥配方?jīng)Q策功能、土壤肥力評價功能等；用戶管理功能主要包括用戶注冊、認證、權(quán)限管理等功能。系統(tǒng)功能結(jié)構(gòu)見圖2。

(2)

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

平均誤差計算公式為：

(3)

式中：m表示潤滑區(qū)域的節(jié)點總數(shù)；pi表示各個算例中的節(jié)點壓強；pi,ref表示參考節(jié)點壓強。

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

(4)

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

(5)

式中：nc表示節(jié)點周圍的單元數(shù)；ncni表示第i個單元的節(jié)點總數(shù)；N代表形函數(shù)；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)可轉(zhuǎn)化為：

(7)

2 形函數(shù)的處理

2.1 常應(yīng)變?nèi)切螁卧?/h3>

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

圖4 三節(jié)點三角形單元的坐標(biāo)轉(zhuǎn)換Fig.4 Coordinate transformation of 3-node triangular element

在節(jié)點周圍的第i個單元內(nèi)，形函數(shù)導(dǎo)數(shù)沿控制體邊界的積分為：

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

(8)

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

利用雅克比矩陣，經(jīng)過推導(dǎo)可得型函數(shù)對全局坐標(biāo)的導(dǎo)數(shù)為：

(9)

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

2.2 六節(jié)點三角形單元

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

在節(jié)點周圍的第i個單元內(nèi)，形函數(shù)導(dǎo)數(shù)沿控制體邊界的積分為(其余推導(dǎo)同2.1):

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

(10)

圖5 六節(jié)點三角形單元的坐標(biāo)轉(zhuǎn)換Fig.5 Coordinate transformation of 6-node triangular element

3 軸承實例分析

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

3.1 算法準確性驗證

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

表1 網(wǎng)格劃分結(jié)果Table 1 Grid meshing results

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

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

3.2 高階單元的結(jié)果對比

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

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

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

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

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

3.3 實例計算

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

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

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

4 結(jié)論

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

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

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

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