李 強， 唐心昊， 張 碩， 王玉君， 許偉偉， 王振波

(中國石油大學(xué)(華東) 新能源學(xué)院，山東青島 266580)

隨著石油機械行業(yè)的發(fā)展，壓縮機、離心泵等旋轉(zhuǎn)機械逐漸大型化、高速化，對滑動軸承性能要求不斷提高，可傾瓦軸承因其高穩(wěn)定性得到了廣泛應(yīng)用?；瑒虞S承性能的研究一直受到學(xué)者們的高度重視。Kayar等[1]和Li等[2]基于Reynolds方程的求解研究了偏心率、轉(zhuǎn)速和潤滑介質(zhì)特性對軸承性能的影響。在求解Reynolds方程時耦合求解能量方程進而研究考慮熱效應(yīng)下固定瓦軸承[3-4]和可傾瓦軸承[5-6]性能。計算機技術(shù)的發(fā)展使求解N-S方程成為可能，Keogh等[7]開發(fā)了CFD程序，研究了軸向開槽的軸承溫度分布；李強等[8-11]采用CFD方法，結(jié)合動網(wǎng)格技術(shù)研究了軸承動靜特性以及軸承內(nèi)不凝結(jié)氣體；Ding等[12]基于CFD方法研究了考慮空氣夾帶下可傾瓦軸承的摩擦損失；Armentrout等[13]分別采用Reynolds方程和CFD方法分析了湍流及對流慣性項對可傾瓦軸承的影響；Parovay等[14]采用CFD方法研究了可傾瓦軸承性能與載荷、轉(zhuǎn)速和徑向間隙的關(guān)系?；赗eynolds方程和CFD方法，逐漸形成了滑動軸承熱流體潤滑(THD)模型。在啟?；蜣D(zhuǎn)子失穩(wěn)的工況下，軸瓦發(fā)生磨損引起其形狀的變化，進而改變軸承性能。磨損工況下軸承性能的研究逐漸展開[15-16]。但目前軸承磨損的研究大多為基于Reynolds方程的固定瓦軸承[17-18]，采用CFD方法的可傾瓦軸承的研究較為少見。筆者采用CFD方法、結(jié)合自編動網(wǎng)格程序建立可傾瓦軸承三維熱流體動壓潤滑模型(thermo-hydrodynamic, THD)，實現(xiàn)潤滑流場、轉(zhuǎn)子運動以及瓦塊轉(zhuǎn)動三者之間的耦合。在此基礎(chǔ)上研究不同磨損深度對可傾瓦軸承潤滑性能的影響。

1 控制方程

1.1 流體控制方程

基于CFD方法對可傾瓦軸承三維流場進行求解，同時考慮潤滑介質(zhì)黏溫效應(yīng)及空化的影響。

質(zhì)量守恒方程為

(1)

式中，ρm為混合相密度，kg/m3；vm為混合相速度矢量，m/s。

動量守恒方程為

(2)

其中

式中，p為壓力，Pa；μm為混合相黏度，Pa·s；g為自由落體加速度矢量，m/s2；F為外體積力，N；ρv為氣相密度，kg/m3。

轉(zhuǎn)軸的偏心以及各瓦塊的轉(zhuǎn)動，導(dǎo)致軸頸與軸承之間油膜間隙分布不均勻，進而產(chǎn)生收斂楔和發(fā)散楔。潤滑油進入軸承發(fā)散楔，油膜壓力降低。當(dāng)油膜壓力低于潤滑油的飽和蒸氣壓時，則會發(fā)生空化，產(chǎn)生氣態(tài)潤滑油。為了更好地預(yù)測軸承性能，采用組分輸運方程描述潤滑介質(zhì)氣、液兩相組分之間的輸運，表示為

(3)

式中，fk為第k相質(zhì)量分?jǐn)?shù)；Re和Rc分別為蒸發(fā)和冷凝速率。

質(zhì)量平均速度和混合密度可描述為

(4)

(5)

式中，ρk為第k相密度，kg/m3；vk為第k相速度矢量，m/s。

組分輸運方程中源項由Singhal等[19]提出的“全空化模型”計算。該模型基于多相流框架，考慮了所有一階效應(yīng)，可以處理與流體相變相關(guān)的大密度變化，而無需預(yù)先確定空化的位置、程度或類型。空化模型描述為

(6)

(7)

式中，Ce和Cc為經(jīng)驗常數(shù)；fg為空氣質(zhì)量分?jǐn)?shù)；pv為空化壓力，Pa。

1.2 熱傳遞控制方程

能量守恒方程為

(8)

式中，v為運動黏度，m2/s；keff為有效導(dǎo)熱率，W/(m·K)；hk為第k相焓，kJ/kg；Shf為體積熱源項，W/m3。

流體與邊界之間能量的傳遞采用對流換熱邊界條件。其熱通量可表示為

q=hf(Tw-Tf)+qrad=hext(Text-Tw).

(9)

式中，q為熱通量，W/m2；hf為流體側(cè)局部傳熱系數(shù)，W/(m2·K)；Tw為邊界表面溫度，K；Tf為局部流體溫度，K；qrad為輻射熱通量，W/m2；hext為外部傳熱系數(shù)，W/(m2·K)；Text為外部散熱溫度，K。

潤滑油黏度對溫度變化十分敏感，黏度隨溫度升高而迅速降低。黏溫效應(yīng)可以用Walther[20]方程表征，其形式為

loglog(v+c)=a-blog(T).

(10)

式中，T表示絕對溫度，K；c通常取0.7；a、b由40、100 ℃時潤滑油黏度計算得到，分別為11.88和4.679。

1.3 邊界條件及離散方法

軸承進出口邊界條件為壓力邊界條件，進口壓力為0.2 MPa，出口壓力為大氣壓，瓦塊上表面導(dǎo)熱系數(shù)為350 W/(m2·℃)，其他為絕熱，固壁為無滑移邊界。基于有限體積法離散控制方程，選擇基于壓力求解器，壓力速度耦合采用SIMPLE算法，壓力差分格式采用PRESTO!。

在本文所有的計算工況下，可傾瓦軸承內(nèi)流場雷諾數(shù)均小于400，因此流動狀態(tài)可認(rèn)為層流。為進行區(qū)分，雷諾數(shù)Ren表示為

Ren=ρvl/μ.

(11)

式中，ρ為潤滑介質(zhì)密度，kg/m3；v為特征速度，m/s；l為特征長度，m；μ為動力黏度，Pa·s。

由式(11)計算最大雷諾數(shù)為238.31。

2 流域模型及計算方法

2.1 流域模型

本文的研究對象為工程上應(yīng)用最廣泛的五瓦塊可傾瓦軸承，結(jié)構(gòu)參考文獻[18]，軸承及潤滑介質(zhì)主要參數(shù)如下：軸承半徑和長度分別為30和35 mm，軸瓦數(shù)為5，支點偏心為0.5，半徑間隙為0.09 mm, 瓦塊包角為57°, 瓦塊厚度為20 mm, 瓦塊轉(zhuǎn)動慣量為0.000 12 kg·m2, 轉(zhuǎn)子質(zhì)量為55.7 kg, 轉(zhuǎn)速為10 200 r/min。三維流域計算模型如圖1所示，包括進油口、瓦塊間間隙、油膜間隙以及瓦塊背

圖1 流域模型Fig.1 Watershed model

部間隙等。軸承加載方式為瓦間加載，供油方式為瓦間徑向供油，潤滑油從軸承兩側(cè)間隙流出。32#潤滑油主要參數(shù)為密度、定壓比熱容、導(dǎo)熱系數(shù)、氣態(tài)密度和空化壓力,其值分別為860 kg/m3、2 025 J/(kg·℃)、0.129 W/(m·℃)、1.2 kg/m3和29 185 Pa。

加載瓦塊油膜厚度較小，容易發(fā)生磨損。因此假設(shè)加載瓦塊(3、4)發(fā)生磨損，磨損位置為瓦塊上表面中心，磨損面積為瓦塊上表面總面積的20%。不同的磨損程度通過磨損深度表征，磨損深度分別為10、20、30、40、50 μm。為滿足網(wǎng)格要求，磨損邊緣采用直邊段過渡。瓦塊磨損示意圖及磨損處網(wǎng)格如圖2所示。

圖2 瓦塊磨損示意圖及網(wǎng)格劃分Fig.2 Schematic diagram and mesh of pad wear

2.2 瞬態(tài)計算過程

在轉(zhuǎn)子由初始位置運動至各工況下靜平衡位置的瞬態(tài)過程中，既包括轉(zhuǎn)子的轉(zhuǎn)動和擾動，也包括各瓦塊的自由擺動。為了準(zhǔn)確模擬上述瞬態(tài)過程，通過自編動網(wǎng)格程序?qū)崿F(xiàn)潤滑流場、轉(zhuǎn)子運動以及各瓦塊轉(zhuǎn)動三者之間的耦合。

假設(shè)轉(zhuǎn)子僅受到由轉(zhuǎn)子質(zhì)量產(chǎn)生的重力及軸承產(chǎn)生的油膜力作用。轉(zhuǎn)子所受合力為零時，可認(rèn)為其處于靜平衡位置。其動力學(xué)方程為

(12)

式中，M為轉(zhuǎn)子質(zhì)量，kg；Fx和Fy分別為油膜力x、y方向的分力，N。

由轉(zhuǎn)子表面油膜壓力積分得

(13)

式中，R為轉(zhuǎn)子半徑，m；L為軸承長度，m；p為油膜壓力，Pa；θ為壓力與x軸正方向夾角，rad。

各瓦塊繞支點的轉(zhuǎn)動角速度由剛體繞定軸轉(zhuǎn)動微分方程得

(14)

式中，Jop為瓦塊對支點的轉(zhuǎn)動慣量，kg·m2；ωp為瓦塊轉(zhuǎn)動角速度，rad/s；Mop(Fp)為瓦塊受到的合力矩，N·m。

瞬態(tài)計算過程如圖3所示，通過式(13)積分得轉(zhuǎn)子表面的瞬態(tài)油膜力，求解式(12)得到轉(zhuǎn)子x、y方向的加速度、速度及位移。由瓦塊上下表面單位面積受到的油膜力及與瓦塊支點的距離得到瓦塊受到的力矩，代入式(14)得到瓦塊的角加速度、角速度及轉(zhuǎn)動角度。通過自編動網(wǎng)格程序(UDF)實現(xiàn)可傾瓦軸承網(wǎng)格更新及瓦塊轉(zhuǎn)動。瓦塊轉(zhuǎn)動后網(wǎng)格如圖4所示，為更好顯示網(wǎng)格更新效果，將軸承間隙及瓦塊背部間隙放大并加大了瓦塊轉(zhuǎn)動角度。在轉(zhuǎn)子移動和各個瓦塊轉(zhuǎn)動后，網(wǎng)格節(jié)點實現(xiàn)了精確移動，即使較大的瓦塊轉(zhuǎn)動角度也能保證較高的網(wǎng)格質(zhì)量。

圖3 瞬態(tài)計算過程Fig.3 Transient calculation process

圖4 網(wǎng)格更新效果Fig.4 Mesh update effect

2.3 網(wǎng)格無關(guān)性驗證

采用的動網(wǎng)格程序需要實時讀取網(wǎng)格節(jié)點位置并進行更新。因此需要對計算模型進行網(wǎng)格無關(guān)性驗證，以在保證計算準(zhǔn)確性的同時提高計算效率。參考文獻[21]，表1給出了偏心率為0.4工況下軸承最大油膜壓力隨油膜間隙網(wǎng)格數(shù)變化情況，在5層徑向網(wǎng)格時軸承最大壓力變化低于1%，因此選擇徑向網(wǎng)格層數(shù)為5層，網(wǎng)格數(shù)為336 000。

表1 網(wǎng)格無關(guān)性

3 結(jié)果與討論

3.1 模型準(zhǔn)確性驗證

將供油溫度37 ℃工況下軸承瓦塊(1、3、4)最高溫度數(shù)值計算結(jié)果與試驗數(shù)據(jù)[18]進行對比，如圖5所示。由圖5可知，不同磨損深度下瓦塊溫度模擬結(jié)果與試驗數(shù)據(jù)吻合較好，加載瓦塊(3、4)溫度隨磨損深度增加逐漸升高，非加載瓦塊(1)反之，最大溫度誤差出現(xiàn)在瓦塊3，數(shù)值為2.62 ℃，小于文獻[18]中瓦塊1最大溫度誤差5 ℃，證明了采用CFD方法建立的三維THD模型可以更好地預(yù)測可傾瓦軸承性能。

圖5 溫度模擬與試驗結(jié)果對比Fig.5 Comparison of temperature simulation and experimental results

3.2 軸承壓力、溫度分布

圖6給出了軸承在供油溫度27 ℃、未發(fā)生磨損工況下瓦塊上表面壓力、溫度分布云圖。由圖6(a)可知，在未發(fā)生磨損的工況下，由于瓦塊的擺動，每個瓦塊均存在高壓區(qū)，產(chǎn)生動壓效應(yīng)，從而提高了可傾瓦軸承的穩(wěn)定性。瓦塊上表面高壓區(qū)沿轉(zhuǎn)子旋轉(zhuǎn)方向(逆時針方向)略有偏移，加載瓦塊(3、4)壓力明顯高于其他瓦塊，起主要承載作用，非加載瓦塊1壓力最低。由圖6(b)可知，在周向方向瓦塊上表面溫度沿轉(zhuǎn)子旋轉(zhuǎn)方向逐漸升高，且在瓦塊后緣處達到最高，瓦塊間間隙供油，潤滑油在間隙內(nèi)混合使得在進入下一油膜間隙時溫度得到了有效降低；在軸向方向，由于供油位置在瓦塊間隙的中間，瓦塊上表面溫度出現(xiàn)了中間低兩邊高的現(xiàn)象；轉(zhuǎn)子偏心導(dǎo)致加載瓦塊(3、4)油膜厚度減小，黏性耗散更加嚴(yán)重，因此其溫升明顯高于非加載瓦塊。

圖6 瓦塊上表面壓力、溫度分布云圖Fig.6 Pressure and temperature contours of pad

3.3 磨損深度對壓力、溫度影響

圖7給出了供油溫度27 ℃、不同磨損深度下各瓦塊上表面壓力及溫度變化，數(shù)據(jù)取自軸承中面與瓦塊上表面交線，其中橫坐標(biāo)表示軸承圓周方向瓦塊位置。

圖7 不同磨損深度下各瓦塊表面壓力及溫度Fig.7 Surface pressure and temperature of pad for various wear depth

由圖7可知，在未發(fā)生磨損工況下，瓦塊表面壓力、溫度分布連續(xù)性變化；在磨損工況下，加載瓦塊(3、4)上表面壓力及溫度的連續(xù)性遭到破壞，且磨損深度越大連續(xù)性破壞越嚴(yán)重。在瓦塊磨損處壓力及溫度產(chǎn)生局部降低，但與未磨損工況相比，加載瓦塊最高壓力和最高溫度明顯升高，非加載瓦塊反之。隨著磨損深度增加，靜平衡位置逐漸下降，加載瓦塊油膜厚度減小，非加載瓦塊油膜厚度增大，因此加載瓦塊的最高壓力和溫度逐漸升高，非加載瓦塊最高壓力和溫度逐漸降低。

3.4 磨損深度對軸承承載力影響

為了研究磨損深度對軸承承載力影響，通過自編動網(wǎng)格程序?qū)⑥D(zhuǎn)子移動至偏心率0.5處，保證潤滑流場、瓦塊轉(zhuǎn)動與轉(zhuǎn)子運動之間的耦合，最終形成穩(wěn)定的油膜間隙和瓦塊轉(zhuǎn)動角，積分轉(zhuǎn)子表面油膜力即為軸承承載力。圖8給出了不同磨損深度下軸承承載力，由圖8可知，隨著磨損深度增加，軸承的承載力不斷降低，但軸承承載力的變化趨勢逐漸減緩。磨損導(dǎo)致瓦塊表面壓力連續(xù)性遭到破壞并引起局部油膜厚度增大，使得加載瓦塊局部壓力降低，因此軸承承載力不斷降低。由于磨損面積僅占瓦塊上表面面積的20%，磨損深度的變化對壓力的影響有限，因此承載力變化趨勢逐漸減緩。在一定的磨損面積(20%)下，軸承承載力對磨損深度的敏感性逐漸降低。

圖8 不同磨損深度下承載力Fig.8 Loading capacity under different wear depths

3.5 磨損深度對偏心率和靜平衡位置影響

圖9給出了供油溫度27 ℃時不同磨損深度下軸心偏心率，由圖9可知:隨著磨損深度增加，軸心的偏心率不斷增大，但變化趨勢逐漸減緩，磨損深度在0～30 μm時，軸心偏心率變化約為10%；而磨損深度在30～50 μm時，軸心偏心率變化率約為2%。隨著磨損深度增加，軸承承載力降低，因此需要更大的偏心以獲得足夠的油膜力支撐轉(zhuǎn)子。在磨損深度為0～30 μm時，一方面，軸承承載力降低且變化較大，另一方面，在較小的偏心率下，由偏心變化引起壓力升高較小，則需要更大的偏心變化獲得足夠的壓力升高，因此在這兩方面因素的共同作用下，偏心率變化較大。磨損深度為30～50 μm時，承載力變化較小，且在高偏心率下，偏心變化引起的壓力升高較大，因此軸心偏心率變化較小。

圖9 不同磨損深度下軸心偏心率Fig.9 Eccentricity under different wear depths

圖10給出了不同磨損深度下軸心運動軌跡。

圖10 不同磨損深度下軸心靜平衡位置及瓦塊轉(zhuǎn)動角Fig.10 Static balance position of shaft and rotation angle of pad under different wear depths

由圖10可知，軸心由起始位置經(jīng)過擾動后最終到達靜平衡位置。隨著磨損深度增加，軸心運動軌跡不斷擴大，靜平衡位置逐漸降低且在x軸方向先增大后減小。為進一步分析靜平衡位置在x軸方向的變化，圖10(b)給出了不同靜平衡位置下各瓦塊的轉(zhuǎn)動角度。由圖10可知，隨著磨損深度增加，加載方向左側(cè)的瓦塊(2、3)由順時針轉(zhuǎn)動逐漸變?yōu)槟鏁r針轉(zhuǎn)動，加載方向右側(cè)的瓦塊(4、5)始終逆時針轉(zhuǎn)動且轉(zhuǎn)動角度不斷增大，由于不同磨損深度下各瓦塊產(chǎn)生不同的轉(zhuǎn)動，改變了轉(zhuǎn)子表面受力，使得軸心靜平衡位置在x軸方向呈現(xiàn)先增大后減小的變化趨勢。

4 結(jié) 論

(1)在所建立的磨損工況下可傾瓦軸承的三維熱流體動壓潤滑模型中，通過自編動網(wǎng)格程序?qū)崿F(xiàn)了潤滑流場、瓦塊轉(zhuǎn)動、轉(zhuǎn)子運動的固體-流體-固體的多向耦合計算，計算結(jié)果與試驗結(jié)果吻合較好，驗證了該計算模型的有效性。

(2)可傾瓦軸承每個瓦塊均存在高壓區(qū)，產(chǎn)生動壓效應(yīng)，瓦塊上表面溫度沿轉(zhuǎn)子旋轉(zhuǎn)方向逐漸升高，且出現(xiàn)中間低兩邊高的現(xiàn)象。

(3)瓦塊的磨損破壞了瓦塊表面壓力和溫度的連續(xù)性，出現(xiàn)了局部的壓力、溫度降低，磨損深度越大，連續(xù)性破壞越嚴(yán)重。

(4)隨著磨損深度增加，軸承承載力逐漸降低，偏心率不斷增大，但變化趨勢逐漸減緩，軸心靜平衡位置逐漸降低且在x軸方向先增大后減小。