基于彈流耦合算法的船舶艉軸承潤(rùn)滑特性分析*

張 斌 朱漢華

(武漢理工大學(xué)能源與動(dòng)力工程學(xué)院 湖北武漢 430063)

大型船舶艉軸承經(jīng)常處于重載工況，交變載荷產(chǎn)生的摩擦熱使?jié)櫥陀湍ず穸葴p小，甚至使軸瓦承載區(qū)無(wú)法形成具有一定承載力的油膜[1]。此時(shí)，彈性變形和油膜厚度非常接近，相對(duì)誤差如此之小的彈性變形勢(shì)必會(huì)對(duì)軸承潤(rùn)滑性能產(chǎn)生重要影響。彈流潤(rùn)滑理論的日趨完善，使廣大學(xué)者漸漸看到了彈性變形的本質(zhì)和影響。

CONWAY和LEE[2]利用彈性位移方程來(lái)研究剛性軸承座中特殊形狀軸瓦表面彈性位移分布的變化，并簡(jiǎn)化了數(shù)值計(jì)算求解彈性變形的過(guò)程。FANTINO和FRENE[3]對(duì)表征內(nèi)燃機(jī)連桿軸承的潤(rùn)滑性能參數(shù)進(jìn)行了迭代求解，對(duì)比分析了是否考慮彈性變形潤(rùn)滑性能產(chǎn)生的差異，結(jié)果表明，計(jì)入彈性變形后油膜壓力呈現(xiàn)雙峰甚至多峰分布，同時(shí)最小油膜厚度近似線性減小。何芝仙和桂長(zhǎng)林[4]優(yōu)化了彈性變形的計(jì)算方法，雖然結(jié)果存在一定誤差，卻大大提高了計(jì)算速度。鄒爭(zhēng)等人[5]通過(guò)數(shù)值計(jì)算的方法得到艉軸承變形的解析解，分析了安裝參數(shù)和工況參數(shù)對(duì)軸承變形的影響，并通過(guò)設(shè)計(jì)正交表來(lái)評(píng)估不同影響因素對(duì)彈性變形的影響程度。馬艷艷和崔子偉[6]發(fā)現(xiàn)應(yīng)力偶軸承更適合研究彈性變形，并在雷諾方程中加入應(yīng)力偶參數(shù)，推導(dǎo)出適合分析應(yīng)力偶軸承變形的雷諾方程。結(jié)果表明：相比于普通軸承，應(yīng)力偶軸承的彈性變形更明顯，且最大彈性變形量與應(yīng)力偶參數(shù)成正比。近年來(lái)，隨著對(duì)彈性變形研究的逐步深入，單一變量求解的方法已經(jīng)不再滿足精度要求，研究人員提出了油膜壓力、厚度、彈性變形的耦合求解算法，使得彈流潤(rùn)滑理論日趨完善。

本文作者以重載工況下的船舶艉軸承為研究對(duì)象，通過(guò)有限單元法結(jié)合耦合算法求解油膜壓力、油膜厚度、彈性變形，進(jìn)一步分析不同影響因素對(duì)艉軸承彈流潤(rùn)滑特性的動(dòng)態(tài)影響趨勢(shì)，得到一些有利于軸承設(shè)計(jì)和安裝的結(jié)論。

1 彈流潤(rùn)滑數(shù)學(xué)模型

1.1 Reynolds方程

對(duì)Navier-Stokes方程(N-S方程)進(jìn)行簡(jiǎn)化得到Reynolds方程[7]：

(1)

式中：h為油膜厚度；p為油膜壓力；η為潤(rùn)滑油動(dòng)力黏度；U為軸頸轉(zhuǎn)速；x為周向坐標(biāo)；y為軸向坐標(biāo)。

式(1)為研究流體潤(rùn)滑提供了數(shù)值計(jì)算的依據(jù)，其意義在于將油膜壓力p、油膜厚度h、潤(rùn)滑油動(dòng)力黏度η、軸頸轉(zhuǎn)速U等看似沒(méi)有關(guān)聯(lián)的參數(shù)有效聯(lián)系起來(lái)，實(shí)現(xiàn)了已知某一參數(shù)來(lái)求解其他未知參數(shù)的目標(biāo)。文中將基于Reynolds方程進(jìn)行潤(rùn)滑分析。

1.2 有限單元法求解軸瓦彈性變形

文中建立了考慮軸瓦厚度的三維有限單元模型，將軸瓦整個(gè)表面的彈性變形分解到每個(gè)節(jié)點(diǎn)來(lái)求解。對(duì)軸瓦模型建立三維坐標(biāo)系，x、y、z分別表示周向、軸向和徑向。所建模型及坐標(biāo)系如圖1所示。

(1)離散和選擇一般單元

在對(duì)多自由度的計(jì)算模型進(jìn)行離散時(shí)通常選用8節(jié)點(diǎn)六面體單元，文中軸瓦模型單元用u、v、w分別表示每個(gè)節(jié)點(diǎn)的3個(gè)自由度，下標(biāo)1～8分別表示六面體的8個(gè)節(jié)點(diǎn)。如圖2所示。

(2)選擇位移函數(shù)

假設(shè)每個(gè)節(jié)點(diǎn)的坐標(biāo)為(xi，yi,zi)，i=1,2,…,8

選擇位移函數(shù)[8]：

u=a1+a2x+a3y+a4z+a5xy+a6xz+a7yz+a8xyz

(2)

各節(jié)點(diǎn)的位移分別表示為u1,u2,…,u8，位移矩陣如式(3)所示。

(3)

定義形狀函數(shù)為N，節(jié)點(diǎn)位移量為S，則：

(4)

將3個(gè)自由度的位移場(chǎng)u、v、w分別表示為向量形式：

(5)

(3)定義應(yīng)變位移和應(yīng)力應(yīng)變的關(guān)系[9]

三維實(shí)體模型的應(yīng)變向量如式(6)所示。

(6)

定義矩陣B(6×24)為

(7)

則應(yīng)力應(yīng)變的關(guān)系可以表示為

ε=BS

(8)

(4)推導(dǎo)單元?jiǎng)偠染仃?/p>

單元?jiǎng)偠染仃囉删仃嘊、D通過(guò)轉(zhuǎn)置變換和積分推導(dǎo)求得[10]：

(9)

式中：Ke為單元?jiǎng)偠染仃?；D為彈性矩陣。

(5)組裝形成整體剛度矩陣

如果沒(méi)有引入適當(dāng)?shù)倪吔鐥l件對(duì)單元?jiǎng)偠染仃囘M(jìn)行自由度限制而直接求解，很難消除多自由度帶來(lái)的不確定性，進(jìn)而引發(fā)位移結(jié)果的多變性。因此，合適的邊界條件對(duì)于剛度矩陣的求解至關(guān)重要[11]。文中引入如下的邊界條件：

①三維有限單元模型表面的等效壓力載荷來(lái)自于軸瓦內(nèi)表面受到的油膜壓力，且油膜壓力連續(xù)分布在整個(gè)軸瓦內(nèi)表面；

②在組裝整體剛度矩陣的過(guò)程中，軸瓦外表面受到的約束為剛性約束，模型各節(jié)點(diǎn)均無(wú)任何自由度的位移[12]。

組裝形成的整體剛度矩陣為

(10)

式中：K為整體剛度矩陣；G為單元節(jié)點(diǎn)轉(zhuǎn)化矩陣。

(6)求解位移

求解位移，需建立有限元方程進(jìn)行求解，進(jìn)而得到整體剛度矩陣中各節(jié)點(diǎn)的位移式[13]。文中建立的有限元方程如下：

KSz=Fz

(11)

式中：Sz為整體結(jié)構(gòu)中各節(jié)點(diǎn)的位移；Fz為整體結(jié)構(gòu)中各節(jié)點(diǎn)的等效載荷。

1.3 彈流潤(rùn)滑耦合算法

彈流潤(rùn)滑耦合算法具體方法為：首先通過(guò)有限差分法求解Reynolds方程得到油膜厚度和壓力，作為彈流潤(rùn)滑油膜厚度和壓力的初值；在初始?jí)毫ψ饔孟螺S瓦會(huì)產(chǎn)生彈性變形，計(jì)入彈性變形后油膜厚度會(huì)改變[14]，此時(shí)需要在初始油膜厚度上加入彈性變形，得到新的油膜厚度；油膜厚度的改變又會(huì)影響油膜壓力，再代入Reynolds方程求解得到新的油膜壓力。以上步驟作為第一次循環(huán)，通過(guò)MATLAB編程進(jìn)行多次迭代，事先設(shè)定收斂精度，每一次循環(huán)結(jié)束后都校核軸承給定外載荷和油膜承載力的誤差[15]，檢驗(yàn)誤差是否滿足收斂精度，這樣才能真正實(shí)現(xiàn)對(duì)潤(rùn)滑性能參數(shù)的耦合求解。耦合算法的求解流程如圖3所示。

圖3 彈流潤(rùn)滑耦合算法流程Fig.3 Flow of elastohydrodynamic lubrication coupling algorithm

2 模型驗(yàn)證

為驗(yàn)證所建彈流潤(rùn)滑數(shù)學(xué)模型的正確性，選取文獻(xiàn)[16]中的全圓軸承參數(shù)進(jìn)行驗(yàn)證，基本參數(shù)如表1所示。

表1 軸承及潤(rùn)滑油基本參數(shù)Table 1 The basic parameters of bearing and lubricating oil

對(duì)比圖4和圖5可知，雖然文中得到的壓力分布圖中從峰值壓力降為0的過(guò)程要比較緩慢，可能是因?yàn)樵谇蠼廨S頸平衡位置的過(guò)程中文獻(xiàn)[16]是采用割線法，而文中采用非線性方程求根的方法，造成一定誤差，但油膜壓力沿周向和軸向的變化趨勢(shì)基本一致，可以證明文中數(shù)學(xué)模型的正確性。

圖4 文中油膜壓力分布Fig.4 Oil film pressure distribution in this paper

圖5 文獻(xiàn)[16]油膜壓力分布Fig.5 Oil film pressure distribution in reference[16]

3 動(dòng)態(tài)參數(shù)對(duì)艉軸承彈流潤(rùn)滑特性的影響

基于建立的船舶艉軸承彈流潤(rùn)滑數(shù)學(xué)模型，將彈性模量、軸承間隙、長(zhǎng)徑比3種參數(shù)作為影響因素，分析了3種動(dòng)態(tài)參數(shù)對(duì)艉軸承彈流潤(rùn)滑特性的影響。所選用的艉軸承相關(guān)參數(shù)如表2所示。

表2 艉軸承及潤(rùn)滑油相關(guān)參數(shù)Table 2 Parameters of stern bearing and lubricationg oil

3.1 動(dòng)態(tài)彈性模量對(duì)艉軸承彈流潤(rùn)滑特性的影響

彈性模量在彈性力學(xué)中定義為物體所受的應(yīng)力和產(chǎn)生應(yīng)變的比例系數(shù)，作為彈性體在彈性變形階段的重要物理量，可以用來(lái)衡量材料產(chǎn)生變形的難易程度，也可以看作材料剛度的體現(xiàn)。圖6—13示出了動(dòng)態(tài)彈性模量對(duì)艉軸承彈流潤(rùn)滑特性的影響，彈性模量的范圍是100～150 GPa。

圖6 動(dòng)態(tài)彈性模量影響下彈性變形三維分布Fig.6 Three dimensional distribution of elastic deformationunder the influence of dynamic elastic modulus：(a)100 GPa;(b)110 GPa;(c)120 GPa;(d)130 GPa;(e)140 GPa;(f)150 GPa

圖7 彈性模量對(duì)最大彈性變形量的影響Fig.7 Influence of elastic modulus on maximum elastic deformation

圖8 動(dòng)態(tài)彈性模量影響下油膜壓力三維分布Fig.8 Three dimensional distribution of oil film pressureunder the influence of dynamic elastic modulus：(a)100 GPa;(b)110 GPa;(c)120 GPa;(d)130 GPa;(e)140 GPa;(f)150 GPa

圖9 彈性模量對(duì)最大油膜壓力的影響Fig.9 Effect of elastic modulus onmaximum oil film pressure

圖10 動(dòng)態(tài)彈性模量影響下油膜厚度三維分布Fig.10 Three dimensional distribution of oil film thicknessaffected by dynamic elastic modulus：(a)100 GPa;(b)110 GPa;(c)120 GPa;(d)130 GPa;(e)140 GPa;(f)150 GPa

圖11 彈性模量對(duì)最小油膜厚度的影響Fig.11 Influence of elastic moduluson minimum oil film thickness

圖12 彈性模量對(duì)摩擦力的影響Fig.12 Effect of elastic modulus on friction

圖13 彈性模量對(duì)端泄流量的影響Fig.13 Effect of elastic modulus on end discharge

對(duì)比圖6和圖8可以發(fā)現(xiàn),彈性變形和油膜壓力的分布趨勢(shì)相似，幾乎都是在周向中點(diǎn)附近取得最大值，沿周向和軸向均呈現(xiàn)先增后減的趨勢(shì)。由圖10可知，油膜厚度在周向中點(diǎn)附近取得最小值。軸承在周向180°附近油膜壓力和油膜厚度會(huì)同時(shí)出現(xiàn)最大值和最小值，因此該區(qū)域潤(rùn)滑環(huán)境復(fù)雜，需要引起關(guān)注。通過(guò)觀察圖7、圖9和圖11可知，彈性模量增加時(shí)，油膜峰值壓力增加，最大彈性變形量和最小油膜厚度均減小。原因在于彈性模量增加，軸瓦表面剛性變強(qiáng)，軸瓦不易發(fā)生形變，相同外部條件下最大彈性變形量減小；軸承承載能力增加，意味著油膜峰值壓力增加；增大油膜壓力會(huì)導(dǎo)致偏心率增大，最小油膜厚度在周向0～180°的區(qū)間內(nèi)出現(xiàn)，由油膜厚度的公式可知，在此區(qū)間內(nèi)油膜厚度與偏心率成反比，因此最小油膜厚度會(huì)減小。

由于黏度的作用，潤(rùn)滑油在軸承間隙流動(dòng)時(shí)會(huì)產(chǎn)生摩擦力。如果摩擦力太大，會(huì)產(chǎn)生過(guò)量的摩擦熱和摩擦功耗，影響軸承工作效率。軸承在工作中會(huì)從軸承兩端泄漏部分油量，稱為端泄流量。為保持一種動(dòng)態(tài)平衡，需要不斷補(bǔ)充和端泄流量相等的潤(rùn)滑油，所以計(jì)算端泄流量的意義在于了解補(bǔ)償流量。由圖12和圖13可知，增大彈性模量，材料的摩擦因數(shù)增大，因此摩擦力增加；端泄流量與油膜壓力和厚度都有關(guān)，雖然膜厚減小，但周向兩端壓力梯度增大的作用更明顯，導(dǎo)致端泄流量小幅度增加。

綜合彈性模量對(duì)艉軸承彈流潤(rùn)滑性能的影響可以看出，在一定區(qū)間內(nèi)增大彈性模量能有效減小軸瓦產(chǎn)生的彈性變形，油膜峰值壓力、摩擦力和端泄流量均增加，但最小油膜厚度會(huì)近似線性減小。

3.2 動(dòng)態(tài)軸承間隙對(duì)艉軸承彈流潤(rùn)滑特性的影響

彈性變形和軸承間隙之間有著不可分割的關(guān)系，而間隙改變會(huì)直接影響到油膜厚度。此外，間隙太小可能使軸頸和軸瓦直接接觸，發(fā)生摩擦磨損，導(dǎo)致潤(rùn)滑失效。圖14—21示出了艉軸承潤(rùn)滑特性隨軸承間隙的動(dòng)態(tài)變化趨勢(shì)，軸承間隙的范圍是0.03～0.08 mm。

圖14 動(dòng)態(tài)軸承間隙影響下彈性變形三維分布Fig.14 Three dimensional distribution of elastic deformationunder the influence of dynamic bearing clearance：(a)0.03 mm;(b)0.04 mm;(c)0.05 mm;(d)0.06 mm;(e)0.07 mm;(f)0.08 mm

圖15 軸承間隙對(duì)最大彈性變形量的影響Fig.15 Influence of bearing clearance onmaximum elastic deformation

圖16 動(dòng)態(tài)軸承間隙影響下油膜壓力三維分布Fig.16 Three dimensional distribution of oil film pressureunder the influence of dynamic bearing clearance：(a)0.03 mm;(b)0.04 mm;(c)0.05 mm;(d)0.06 mm;(e)0.07 mm;(f)0.08 mm

圖17 軸承間隙對(duì)最大油膜壓力的影響Fig.17 Influence of bearing clearanceon maximum oil film pressure

圖19 軸承間隙對(duì)最小油膜厚度的影響Fig.19 Influence of bearing clearance on minimum oil film thickness

圖20 軸承間隙對(duì)摩擦力的影響Fig.20 Influence of bearing clearance on friction

圖21 軸承間隙對(duì)端泄流量的影響Fig.21 Influence of bearing clearance on end discharge

由圖14—19可知，彈性變形、油膜壓力沿周向和軸向均近似呈拋物線分布，且基本分布趨勢(shì)不受軸承間隙變化的影響。軸承間隙改變會(huì)直接影響彈性變形量，進(jìn)而改變油膜壓力、厚度，可以看出，當(dāng)軸承間隙增加時(shí)，油膜峰值壓力增加，最大彈性變形量和最小油膜厚度則線性減小，與軸承間隙呈現(xiàn)出較強(qiáng)的負(fù)相關(guān)關(guān)系。從圖20、21可以看出，摩擦力、端泄流量隨軸承間隙增大而增大，但增加幅度不大，均在10%以內(nèi)。這主要是因?yàn)樵黾虞S承間隙和增大材料的泊松比效果相似，會(huì)導(dǎo)致軸瓦的豎直變形量增加，因此最大彈性變形量增加。在定偏心距時(shí)，由偏心率和軸承間隙的關(guān)系可知，間隙增大會(huì)引起偏心率減小，而偏心率是相關(guān)潤(rùn)滑參數(shù)的重要影響因素，其值減小會(huì)直接引起油膜峰值壓力、摩擦力、端泄流量增加，最小油膜厚度減小，尤其是大偏心率情況下，潤(rùn)滑參數(shù)的變化會(huì)更明顯。

綜上所述，軸承間隙增加可以有效減小最大彈性變形，提高油膜峰值壓力，但最小油膜厚度會(huì)減小，而摩擦力和端泄流量受軸承間隙的影響不大。因此在軸承安裝時(shí)，需要將軸承間隙控制在合理的區(qū)間內(nèi)，確保軸承處于優(yōu)良的潤(rùn)滑環(huán)境。

3.3 動(dòng)態(tài)長(zhǎng)徑比對(duì)艉軸承彈流潤(rùn)滑特性的影響

長(zhǎng)徑比作為材料重要的設(shè)計(jì)參數(shù)之一，通常用來(lái)描述材料的形狀，長(zhǎng)徑比大的材料為細(xì)長(zhǎng)狀，長(zhǎng)徑比小的材料為短粗狀。圖22—29示出了艉軸承彈流潤(rùn)滑特性隨長(zhǎng)徑比的動(dòng)態(tài)變化趨勢(shì)，長(zhǎng)徑比的范圍是0.5～5。

圖22 動(dòng)態(tài)長(zhǎng)徑比影響下彈性變形三維分布Fig.22 Three dimensional distribution of elastic deformationunder the influence of dynamic aspect ratio：(a)L/d=0.5;(b)L/d=1;(c)L/d=2;(d)L/d=3;(e)L/d=4;(f)L/d=5

圖23 長(zhǎng)徑比對(duì)最大彈性變形量的影響Fig.23 Influence of length diameter ratioon maximum elastic deformation

圖27 長(zhǎng)徑比對(duì)最小油膜厚度的影響Fig.27 Influence of aspect ratio on minimum oil film thickness

圖28 長(zhǎng)徑比對(duì)摩擦力的影響Fig.28 Effect of length diameter ratio on friction

圖29 長(zhǎng)徑比對(duì)端泄流量的影響Fig.29 Effect of length diameter ratio on end discharge

從圖22、24可看出，彈性變形與油膜壓力出現(xiàn)最大值的位置仍在周向中點(diǎn)附近。從圖23、25、27可知，當(dāng)長(zhǎng)徑比從0.5線性增加至5時(shí)，最大彈性變形量近似線性增加，油膜峰值壓力則剛好相反，而最小油膜厚度僅增加了0.001 3 mm，說(shuō)明長(zhǎng)徑比對(duì)最小油膜厚度的影響較小。從圖28、29可知，摩擦力和端泄流量都與長(zhǎng)徑比呈現(xiàn)出一種負(fù)相關(guān)關(guān)系?？傮w來(lái)看，增大長(zhǎng)徑比，會(huì)使軸瓦內(nèi)表面產(chǎn)生較大的彈性變形，軸承間隙減小，進(jìn)而降低油膜峰值壓力、摩擦力和端泄流量。產(chǎn)生上述變化的原因在于增大長(zhǎng)徑比，軸瓦受力的比面積減小，相同外載條件下會(huì)產(chǎn)生更大的形變，因此最大彈性變形量增加；改變長(zhǎng)徑比的本質(zhì)在于改變了軸承形狀，大長(zhǎng)徑比的軸承會(huì)使油膜壓力分布更均勻，比壓不變時(shí)，油膜峰值壓力減小，最小油膜厚度略有增加。此外，長(zhǎng)徑比也是偏心率的重要影響因素之一，增大長(zhǎng)徑比致使偏心率增加，而摩擦力和端泄流量都會(huì)隨偏心率的增加同步減小。因此，在設(shè)計(jì)艉軸承長(zhǎng)徑比時(shí)，應(yīng)綜合考慮艉軸承在重載工況下可能產(chǎn)生的彈性變形以及彈性變形對(duì)潤(rùn)滑特性的影響。

4 結(jié)論

基于彈流潤(rùn)滑耦合算法分析了艉軸承的彈流潤(rùn)滑特性，探討彈性模量、軸承間隙、長(zhǎng)徑比3種因素對(duì)艉軸承彈流潤(rùn)滑特性的影響，得到以下結(jié)論：

(1)彈性變形和油膜壓力的三維分布趨勢(shì)相似，沿周向和軸向都近似呈拋物線分布，最大值出現(xiàn)在周向中點(diǎn)附近。

(2)隨彈性模量增加，油膜峰值壓力增加，最大彈性變形量和最小厚度均減小，摩擦力和端泄流量同時(shí)增加；隨軸承間隙增加，油膜峰值壓力穩(wěn)步提高，最大彈性變形量和最小油膜厚度卻剛好相反，近似線性減小，而摩擦力、端泄流量受軸承間隙的影響較?。浑S長(zhǎng)徑比增加，彈性變形與油膜壓力的分布更為平緩，最大彈性變形量近似線性增加，油膜峰值壓力減小，最小油膜厚度變化不明顯，而摩擦力和端泄流量都與長(zhǎng)徑比呈現(xiàn)出一種負(fù)相關(guān)關(guān)系。