基于格點型有限體積法的滑動軸承油膜壓力特性分析

2019-09-17 06:46:26楊國棟曹貽鵬明平劍張文平李燎原

振動與沖擊 2019年17期

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

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

滑動軸承是最常見的流體潤滑軸承，具有形式簡單、接觸面積和承載力大等優(yōu)點，在工業(yè)領(lǐng)域中廣泛應(yīng)用[1]。當(dāng)軸頸轉(zhuǎn)動時，軸與軸瓦之間的楔形間隙內(nèi)會形成一層動壓液膜，起到平衡外載荷和減小摩擦損耗的作用。

雷諾方程是求解軸承油膜壓力的核心方程，其求解方法主要有解析法、有限差分法(FDM)、有限元法(FEM)等[2]。解析法只適用于無限短或無限長軸承，此時雷諾方程退化為一維方程，無法準(zhǔn)確描述軸承內(nèi)的壓力分布；FDM是最常用的數(shù)值解法，形式簡單，方程中的偏導(dǎo)數(shù)項可以通過差商形式表示，潤滑區(qū)域的離散采用正交化的四邊形網(wǎng)格，Sun等[3-10]基于該算法研究了滑動軸承的潤滑特性，但是對于結(jié)構(gòu)復(fù)雜的潤滑區(qū)域，無法采用該算法求解；FEM是從彈性力學(xué)發(fā)展而來的算法，之后應(yīng)用于潤滑計算領(lǐng)域[11-15]，潤滑區(qū)域的劃分可以采用非結(jié)構(gòu)化網(wǎng)格，單元大小和節(jié)點數(shù)可以任意選取，對潤滑區(qū)域形狀的要求較少，但是該算法需要構(gòu)建一個總剛矩陣，計算量較大。

格點型有限體積法(Cell-Vertex Finite Volume Method，下面用CVFVM表示)是CFD領(lǐng)域中常用的方法，該算法既具有有限體積法強(qiáng)守恒性的特點，又具有FEM對網(wǎng)格類型靈活性的特點，近些年來在傳熱和+結(jié)構(gòu)等領(lǐng)域也有較多的應(yīng)用[16-17]，但是在潤滑領(lǐng)域應(yīng)用較少。

本文基于CVFVM離散雷諾方程，并利用矢量化方式推導(dǎo)得出了適用于非結(jié)構(gòu)化網(wǎng)格的方程形式，并對文獻(xiàn)[3]中的算例進(jìn)行了對比和驗證，在此基礎(chǔ)上，研究了偏心率和軸頸傾斜角度對內(nèi)燃機(jī)主軸承壓力特性的影響。

1 理論方程

1.1 雷諾方程

在潤滑現(xiàn)象的研究中，雷諾方程的求解是關(guān)鍵步驟，不可壓縮流體潤的瞬態(tài)雷諾方程如下

(1)

式中:θ為周向角度坐標(biāo)，degree;z為軸承寬度坐標(biāo)，m;R為軸承半徑，m;p為油膜壓力，Pa;h為油膜厚度，m;U為軸頸速度，m/s;η為潤滑油黏度，Pa·s。

根據(jù)下式進(jìn)行無量綱化

(2)

式中:B為軸承寬度，m;c為軸承間隙，m；φ偏位角，degree;ε為偏心率。

可得

(3)

為了使式(3)適用于非結(jié)構(gòu)化網(wǎng)格，雷諾方程轉(zhuǎn)化為矢量形式如下

▽·(Γ▽p)=Φs

(4)

其中

式中:Γ為壓力系數(shù)；Φs為源項。

1.2 液膜厚度方程

設(shè)軸心坐標(biāo)為(x,y)，徑向滑動軸承的液膜厚度可用下式表示

h=c+ecos(θ-φ)=

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

c+xcosθ+ysinθ

(5)

式中:e為偏心距。

無量綱形式為

H=1+Xcosθ+Ysinθ

(6)

X=x/c;Y=y/c

(7)

油膜厚度位置對坐標(biāo)和時間的導(dǎo)數(shù)，可用如下形式表示

(8)

因此，源項的另一種表達(dá)如下

Φs=6(-Xsinθ+Ycosθ)

(9)

2 雷諾方程的離散

2.1 方程的積分形式

在控制體內(nèi)，對式(4)進(jìn)行積分可得

(10)

式中:V為控制體體積。

2.2 方程的離散

CVFVM算法將有限體積法(FVM)的守恒性和有限單元法(FEM)處理非結(jié)構(gòu)單元的靈活性結(jié)合在了一起。其控制體是圍繞節(jié)點構(gòu)成，節(jié)點儲存壓力P，其他變量存儲在單元中心，如圖1所示。

圖1 CVFVM 控制體示意圖

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

(11)

式中:nc為節(jié)點周圍的單元數(shù)；ncni為第i個單元的節(jié)點總數(shù)；j為單元內(nèi)的節(jié)點編號；N為形函數(shù)；α為x和y兩個坐標(biāo)方向；S為圍繞控制體的積分線。

假設(shè)源項在單元內(nèi)均勻分布

(12)

最終，式(10)可轉(zhuǎn)化為

(13)

2.3 形函數(shù)積分的求解

本算法目前采用兩種單元形式：常應(yīng)變?nèi)切螁卧碗p線性四邊形單元，下面分別給出各個單元形函數(shù)導(dǎo)數(shù)積分的推導(dǎo)過程：

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

三角形單元的示意圖如圖2所示。

三角形單元形函數(shù)

N1=1-ξ-η;N2=ξ;N3=η

(14)

式中:ξ,η為局部坐標(biāo)系中的坐標(biāo)。

整體坐標(biāo)系下的坐標(biāo)可用下式表示

圖2 任意三角形單元的坐標(biāo)轉(zhuǎn)換

y=y1N1+y2N2+y3N3

(15)

其中，(x1,y1)、(x2,y2)、(x3,y3)為整體坐標(biāo)系下三個節(jié)點的坐標(biāo)。

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

(16)

式中:j為單元的節(jié)點編號；Lijα為第i個單元圍繞第j個節(jié)點的積分線在α方向上的投影長度。

形函數(shù)在全局坐標(biāo)下的導(dǎo)數(shù)如下

(17)

各條積分線的長度可用式(18)求解

(18)

2.3.2 雙線性四邊形單元

四邊形單元的示意圖如圖3所示。

圖3 任意直邊四邊形單元的坐標(biāo)轉(zhuǎn)換

四邊形單元形函數(shù)

(19)

整體坐標(biāo)可用下式表示

x=x1N1+x2N2+x3N3+x4N4

y=y1N1+y2N2+y3N3+y4N4

(20)

式中：(x1,y1)、(x2,y2)、(x3,y3)、(x4,y4) 分別為四邊形節(jié)點1,2,3,4在全局坐標(biāo)系下的坐標(biāo)。

形函數(shù)圍繞節(jié)點1的線積分

(21a)

形函數(shù)圍繞節(jié)點2的線積分

(21b)

形函數(shù)圍繞節(jié)點3的線積分

(21c)

形函數(shù)圍繞節(jié)點4的線積分

(21d)

線積分的投影長度可以用下式表示

Lx|AO=(x3+x4-x1-x2)/4

Ly|AO=(y3+y4-y1-y2)/4

Lx|BO=(x1+x4-x2-x3)/4

Ly|BO=(y1+y4-y2-y3)/4

Lx|CO=(x1+x2-x3-x4)/4

Ly|CO=(y1+y2-y3-y4)/4

Lx|DO=(x2+x3-x1-x4)/4

Ly|DO=(y2+y3-y1-y4)/4

(22)

3 算例驗證

本文的計算程序是在哈爾濱工程大學(xué)動力裝置工程技術(shù)研究所自主開發(fā)的通用輸運方程求解器GTEA軟件的基礎(chǔ)上開發(fā)的，采用Fortran90語言編程，在參數(shù)為4 GB RAM、Intel Core i5-2400、CPU3.1 GHz的計算機(jī)中運行程序，采用Gambit等網(wǎng)格劃分軟件進(jìn)行網(wǎng)格劃分，通過GTEA的前處理模塊讀入網(wǎng)格信息。為了驗證該算法的適用性和精度等特點，本文采用文獻(xiàn)[3]中的徑向滑動軸承算例，分別計算了不同工況下的壓力結(jié)果，軸承參數(shù)如表1所示。

3.1 壓力分布的對比

本節(jié)分別選用0.01°和0.03°工況下的壓力分布結(jié)果進(jìn)行對比，計算模型的偏位角均為90°。具體結(jié)果如圖4～圖6所示。

表1 軸承基本參數(shù)

(a) 參考值

(b) 當(dāng)前值

(a) 參考值

(b) 當(dāng)前值

(a) 傾斜角=0.01°

(b) 傾斜角=0.03°

由圖4～圖5可知，不同傾角情況下，主壓力區(qū)均出現(xiàn)在90°～270°之間，與偏位角值相對應(yīng)；傾斜角為0.01°時，壓力在軸承中截面處最大，并沿兩側(cè)方向之間減小；傾斜角為0.03°時，壓力在軸承兩個端面處最大，并沿中間位置方向逐漸減小。

由于當(dāng)傾斜角為0.01°時，最大壓力值出現(xiàn)在軸承中截面附近，因此取33 mm處的壓力分布進(jìn)行對比，當(dāng)傾斜角為0.03°時，最大壓力出現(xiàn)在軸端端部附近，因此取63.8 mm處的壓力值進(jìn)行對比。如圖6所示，在兩個傾斜角情況下，CVFVM算法的結(jié)果與參考值基本一致。

由表2可知，在不同傾角下，CVFVM算法與參考結(jié)果下的做大壓力值基本一致，最大誤差0.67%。

表2 峰值壓力的對比

經(jīng)過上面的對比分析，說明CVFVM用于雷諾方程的求解是適用的，同時也能保證較好的計算精度。

3.2 網(wǎng)格適應(yīng)性的驗證

由于雷諾方程的離散過程是采用矢量化方式推導(dǎo)的，因此該算法對于非結(jié)構(gòu)化網(wǎng)格有很好的適應(yīng)性，同時，基于有限體積法的強(qiáng)守恒特性，也可通過合理分配求解域的網(wǎng)格密度來提高計算效率。下面分別采用不同密度的四邊形網(wǎng)格和三角形網(wǎng)格劃分計算域，進(jìn)而對上述兩個特點進(jìn)行驗證。表3給出了不同模型的網(wǎng)格參數(shù)及計算耗時，其中Quad表示四邊形網(wǎng)格，Tri代表三角形網(wǎng)格，計算結(jié)果如圖7所示。

表3 不同算例的網(wǎng)格參數(shù)及耗時

(a) Quad_1

圖7中，方框區(qū)域內(nèi)油膜壓力值及梯度值較大，在油膜壓力分析中為重點研究的區(qū)域，因此，保持該區(qū)域網(wǎng)格密度不變，分別逐步減小該區(qū)域兩側(cè)的網(wǎng)格密度。由圖7可知，無論采用四邊形網(wǎng)格還是三角形網(wǎng)格劃分區(qū)域，油膜壓力分布基本一致；同時由表3可知，隨著兩側(cè)網(wǎng)格密度的減小，節(jié)點數(shù)逐漸下降，相應(yīng)模型的計算時間逐漸減小。由此，CVFVM算法對非結(jié)構(gòu)化網(wǎng)格的適應(yīng)性及對計算效率的提高得以驗證。