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

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

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

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

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

矢量形式的雷諾方程為：

(1)

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

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

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

h=c+ecos(θ-φ)=

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

c+xcosθ+ysinθ

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

(2)

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

平均誤差計(jì)算公式為：

(3)

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

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

(4)

擴(kuò)散項(xiàng)可按下式進(jìn)行離散：

(5)

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

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

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

源項(xiàng)離散為：

(6)

則式(4)可轉(zhuǎn)化為：

(7)

2 形函數(shù)的處理

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

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

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

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

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

(8)

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

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

(9)

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

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

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

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

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

(10)

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

3 軸承實(shí)例分析

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

3.1 算法準(zhǔn)確性驗(yàn)證

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

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

表2 油膜壓強(qiáng)最大值的對(duì)比

圖6 Case0壓強(qiáng)分布準(zhǔn)確性驗(yàn)證Fig.6 Verification of the accuracy of pressure distribution of Case0

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

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

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

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

圖7 不同算例壓強(qiáng)分布對(duì)比Fig.7 Comparison of pressure distribution of different cases

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

3.3 實(shí)例計(jì)算

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

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

圖8 傾角對(duì)軸承壓強(qiáng)分布的影響Fig.8 Effects of tilt angle on oil film pressure distribution

4 結(jié)論

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

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

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

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