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

(哈爾濱工程大學 動力與能源工程學院，哈爾濱 150001)

滑動軸承由于形式簡單，接觸面積和承載能力大，廣泛應用于工業(yè)生產(chǎn)和生活領(lǐng)域中。當前，滑動軸承動壓潤滑現(xiàn)象的求解方法主要包括有限差分法[1]、有限元法[2]及CFD軟件仿真[3-4]等。有限差分方法(Finite Differential Method，F(xiàn)DM)是計算機數(shù)值模擬最早采用的方法，發(fā)展也較為成熟，但是FDM對求解域的幾何特征要求較高，其結(jié)構(gòu)化網(wǎng)格不易處理復雜的區(qū)域，同時在網(wǎng)格數(shù)較少時，較難滿足積分守恒；有限元法(Finite Element Method， FEM)主要應用于結(jié)構(gòu)力學領(lǐng)域，近些年開始應用于流體力學領(lǐng)域，能夠采用非結(jié)構(gòu)化網(wǎng)格處理幾何特征復雜的區(qū)域，但是FEM需要求解剛度矩陣，當單元數(shù)較多時，會占用較大的計算機內(nèi)存，進而使求解速度下降[5]；一些學者也采用CFD軟件仿真的方法來研究軸承的潤滑性能，通過求解N-S方程來得到軸承的壓力分布，這種方法更為準確，但是求解較復雜，而且油膜厚度是微米級，厚度方向的網(wǎng)格劃分較為困難[6]。

李強等[7]基于FVM(Finite Volume Method)研究了考慮氣穴現(xiàn)象的滑動軸承油膜壓力特性，但是其方程推導是基于結(jié)構(gòu)化網(wǎng)格進行的，這種方法本質(zhì)上有限差分法一致，沒有體現(xiàn)出來有限體積法對復雜網(wǎng)格的適應性。

格心型有限體積法(Cell-Centered Finite Volume Method， CCFVM)是基于CFD技術(shù)中的FVM的離散方法建立起來的[8]，所有變量定義在單元中心，要求因變量的積分守恒對任意一組控制體積都得到滿足，進而對整個計算區(qū)域也滿足守恒；該方法對復雜的求解域有很好的適應性，當網(wǎng)格數(shù)較多時，處理效率也比較高，因此該方法廣泛應用于流場數(shù)值計算和力學相關(guān)領(lǐng)域，而在潤滑領(lǐng)域應用較少。

本文基于CCFVM離散雷諾方程，推導適用于求解非結(jié)構(gòu)化網(wǎng)格的方程離散形式，將計算域的網(wǎng)格要求由結(jié)構(gòu)化網(wǎng)格擴展到非結(jié)構(gòu)化網(wǎng)格，可用于復雜幾何表面之間的潤滑特性研究，并對文獻[9]中的算例進行了驗證，在此基礎(chǔ)上，以水潤滑橡膠軸承為例，分析了長徑比、偏心率對其潤滑特性的影響。

1 理論方程

1.1 雷諾方程

雷諾方程是由簡化的納維-斯托克斯方程和連續(xù)方程推導而得的，常用來求解轉(zhuǎn)子機構(gòu)的動靜壓潤滑問題，適用于二維、穩(wěn)態(tài)、不可壓流體動壓潤滑計算的雷諾方程形式為[10]

(1)

式中：x為圓周方向坐標；y為軸承寬度方向；h為油膜厚度；η為滑油黏度；U為軸頸轉(zhuǎn)速；P為油膜壓力。

1.2 液膜厚度方程

滑動軸承液膜厚度h的計算公式為

h=c+ecosφ

(2)

(3)

式中：c為軸承間隙；e為軸承偏心距；φ為軸承圓周方向角坐標；R為軸承半徑。

2 數(shù)值方法分析

2.1 雷諾方程的矢量表達形式

雷諾方程式(1)可以用如下形式表示

(4)

對式(4)，在空間上對控制體V進行積分，有

(5)

2.2 擴散項的離散[11]

根據(jù)高斯散度定理，將體積分轉(zhuǎn)化為面積分并離散可得

(6)

式中：nf為當前單元周圍的單元個數(shù)；Af為面積矢量Af的模。

圖1 CCFVM控制體示意圖Fig.1 Sketch map of the control volume of CCFVM

圖1中：ef為兩個單元中心的距離矢量；d=ef/ef為單位距離矢量；n為單元面的單位外法線矢量。

由于

(7)

將式(7)代入式(6)中可得

(8)

單元面上的變量梯度由當前單元中心和相鄰單元中心的變量梯度線性插值得到。

2.3 源項的處理

式(5)中的源項在控制體內(nèi)積分可得

(9)

2.4 雷諾方程的離散形式

由式(8)和式(9)可得式(1)的離散形式為

(10)

2.5 邊界條件

對于徑向滑動軸承，其求解域可展開為矩形，采用雷諾邊界條件

(11)

式中：x′為待定的出口邊界；B為軸承的寬度。

邊界上變量直接賦值為Dirichlet邊界條件，邊界上變量梯度直接賦值稱之為Neumann邊界條件。這兩類邊界條件對式(10)中的貢獻為

式中：nD為位于Dirichlet邊界SD上所有單元面的個數(shù)；nN為位于Neumann邊界SN上所有單元面的個數(shù)； 下標Bf表示邊界面。

2.6 收斂條件

CCFVM算法的收斂條件為

(12)

式中：ξ=1×10-6。

3 算例驗證與分析

本文的計算程序是在哈爾濱工程大學動力裝置工程技術(shù)研究所自主開發(fā)的通用輸運方程求解器GTEA軟件的基礎(chǔ)上開發(fā)的，采用Fortran語言編程，網(wǎng)格在Gambit 2.4.6中進行劃分。

3.1 驗證算例

為了驗證CCFVM求解二維穩(wěn)態(tài)潤滑問題的適用性和準確性，本文將通過Sun等研究中的算例來驗證。

表1 徑向滑動軸承基本參數(shù)

采用非結(jié)構(gòu)化網(wǎng)格劃分求解域，共1 556個節(jié)點，2 946個三角形單元；孫軍采用5 400個節(jié)點，5 236個單元。本算例中偏位角為0°，軸頸傾角分別為0°，0.004°，0.007°，0.01°。

計算完畢后，以Sun等論文中的壓力結(jié)果(圖2～圖5中左側(cè)的圖)為參考，進行以下對比：

由圖2～圖5和表2可知，當前算法的液膜壓力分布同參考值趨勢一致，最大壓力值基本相當，驗證了CCFVM算法在求解雷諾方程及滑動軸承潤滑特性方面的適用性和準確性。

圖2 傾斜角為0°Fig.2 Slope angle=0°

圖3 傾斜角為0.004°Fig.3 Slope angle=0.004°

圖4 傾斜角為0.007°Fig.4 Slope angle=0.007°

圖5 傾斜角為0.01°Fig.5 Slope angle=0.01°

參數(shù)數(shù)值傾角/(°)00.0040.0070.01參考值/MPa33.0639.663.58415.35當前值/MPa33.0139.3763.62412.7偏差/%0.150.580.060.72

3.2 實例計算

以水潤滑橡膠軸承為例，研究不同工況下的軸承潤滑特性。采用CCFVM研究軸承長徑比、偏心率的變化對軸承潤滑特性的影響。

表3 艉軸承的基本參數(shù)

為了研究長徑比的變化對軸承壓力分布的影響，保持軸承寬度不變，長徑比分別為1，2，3，4，計算得到的壓力分布如圖6所示。

由圖6可知，沿軸承圓周方向，液膜壓力先增大后減??；隨著長徑比的增大，高壓區(qū)沿軸承寬度方向逐漸擴張，壓力峰值逐漸下降，說明長徑比增大會降低軸承的峰值壓力，同時提高軸承的承載力。

為了研究偏心率對液膜壓力的影響，保持長徑比為1，其他參數(shù)不變，偏心率分別取0.1,0.2,0.3，0.4,05,0.6,0.7,0.8,0.9，結(jié)果如圖7所示。

由圖7可知，隨著偏心率的增大，液膜壓力峰值逐漸增大，偏心率在0.7～0.9變化時，壓力增大的更為劇烈；在圓周方向，壓力的分布遵循收斂楔的幾何變化，先增大后減小。

圖6 長徑比對液膜壓力分布的影響Fig.6 Effects of aspect ratio on film pressure distribution

圖7 偏心率對液膜壓力的影響Fig.7 Effects of eccentricity ratio on film pressure

4 結(jié) 論

(1) 提出了數(shù)值求解雷諾方程的CCFVM算法，該算法物理意義更明確，由于采用了非結(jié)構(gòu)化網(wǎng)格處理技術(shù)，可以適用結(jié)構(gòu)復雜的潤滑區(qū)域，同時不會消耗太多存儲空間，相比FDM(結(jié)構(gòu)化網(wǎng)格)和FEM(消耗較大存儲空間)具有明顯的優(yōu)勢。

(2) 通過與文獻的結(jié)果對比，液膜壓力分布和趨勢基本吻合，最大液膜壓力值接近，驗證了該算法的適用性和正確性較好。

(3) 隨著傾角的增大，壓力峰值逐漸向軸承尾端偏移，同時最大壓力值急劇增大；長徑比越大，液膜高壓區(qū)沿軸承寬度方向擴張的越大，同時壓力峰值會減?。黄穆试酱?，軸承壓力峰值越大，當偏心率大于0.7時，峰值變化更為明顯。