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

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

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

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

格點(diǎn)型有限體積法(Cell-Vertex Finite Volume Method，下面用CVFVM表示)是CFD領(lǐng)域中常用的方法，該算法既具有有限體積法強(qiáng)守恒性的特點(diǎn)，又具有FEM對網(wǎng)格類型靈活性的特點(diǎn)，近些年來在傳熱和+結(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)行了對比和驗(yà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為源項(xiàng)。

1.2 液膜厚度方程

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

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)和時(shí)間的導(dǎo)數(shù)，可用如下形式表示

(8)

因此，源項(xiàng)的另一種表達(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é)點(diǎn)構(gòu)成，節(jié)點(diǎn)儲(chǔ)存壓力P，其他變量存儲(chǔ)在單元中心，如圖1所示。

圖1 CVFVM 控制體示意圖

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

(11)

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

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

(12)

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

(13)

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

本算法目前采用兩種單元形式：常應(yīng)變?nèi)切螁卧碗p線性四邊形單元，下面分別給出各個(gè)單元形函數(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)系下三個(gè)節(jié)點(diǎn)的坐標(biāo)。

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

(16)

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

形函數(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é)點(diǎn)1,2,3,4在全局坐標(biāo)系下的坐標(biāo)。

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

(21a)

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

(21b)

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

(21c)

形函數(shù)圍繞節(jié)點(diǎn)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 算例驗(yàn)證

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

3.1 壓力分布的對比

本節(jié)分別選用0.01°和0.03°工況下的壓力分布結(jié)果進(jìn)行對比，計(jì)算模型的偏位角均為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°時(shí)，壓力在軸承中截面處最大，并沿兩側(cè)方向之間減??；傾斜角為0.03°時(shí)，壓力在軸承兩個(gè)端面處最大，并沿中間位置方向逐漸減小。

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

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

表2 峰值壓力的對比

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

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

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

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

(a) Quad_1

(c) Quad_3

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

4 內(nèi)燃機(jī)主軸承壓力特性求解

主軸承是內(nèi)燃機(jī)結(jié)構(gòu)中的關(guān)鍵部件，其壓力特性的求 解是整機(jī)潤滑性能分析的關(guān)鍵組成部分，本節(jié)選用的主軸承參數(shù)如下表4所示，主要用于分析偏心率和軸頸傾角對壓力特性的影響，其中偏位角為0°。

表4 內(nèi)燃機(jī)主軸承參數(shù)

由圖8可知，隨著偏心率的增大，最小油膜厚度逐漸減小，同時(shí)最大油膜壓力逐漸增大，在偏心率在0.7～0.9之間變化時(shí)，液膜壓力的增大更為明顯。

(a) 油膜厚度的變化

(b) 油膜壓力的變化

由圖9可知，隨著軸頸傾斜角的增大，軸承端部附近的油膜厚度逐漸較小，油膜壓力的峰值由中間位置逐漸向軸承端部移動(dòng)，同時(shí)液膜壓力的幅值逐漸增大。

(a) 傾斜角0.01°

(c) 傾斜角0.03°

5 結(jié) 論

本文基于格點(diǎn)型有限體積法，并采用矢量化方式離散和求解了雷諾方程，通過對不同工況下滑動(dòng)軸承潤滑算例的對比和分析，得出以下結(jié)論：

(1) CVFVM算法用于求解雷諾方程是合理可行的，同時(shí)計(jì)算精度滿足要求。

(2) 與文獻(xiàn)[3]方法相比，CVFVM算法可用于處理非結(jié)構(gòu)化網(wǎng)格劃分的求解域，因此可用于計(jì)算幾何形狀復(fù)雜的潤滑區(qū)域。

(3) 由于有限體積法對局部區(qū)域強(qiáng)守恒性的特點(diǎn)，可以通過對潤滑區(qū)域網(wǎng)格密度的合理分配來提高計(jì)算效率。

(4) 主軸承的偏心距的增大會(huì)引起軸承油膜壓力峰值的增大，軸頸傾斜角的增大會(huì)導(dǎo)致峰值壓力向端部移動(dòng)，會(huì)加快軸承端部的磨損。