李 強(qiáng)， 唐心昊， 張 碩， 王玉君， 許偉偉， 王振波

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

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

1 控制方程

1.1 流體控制方程

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

質(zhì)量守恒方程為

(1)

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

動(dòng)量守恒方程為

(2)

其中

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

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

(3)

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

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

(4)

(5)

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

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

(6)

(7)

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

1.2 熱傳遞控制方程

能量守恒方程為

(8)

式中，v為運(yùn)動(dòng)黏度，m2/s；keff為有效導(dǎo)熱率，W/(m·K)；hk為第k相焓，kJ/kg；Shf為體積熱源項(xiàng)，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 ℃時(shí)潤滑油黏度計(jì)算得到，分別為11.88和4.679。

1.3 邊界條件及離散方法

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

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

Ren=ρvl/μ.

(11)

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

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

2 流域模型及計(jì)算方法

2.1 流域模型

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

圖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)計(jì)算過程

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

假設(shè)轉(zhuǎn)子僅受到由轉(zhuǎn)子質(zhì)量產(chǎn)生的重力及軸承產(chǎn)生的油膜力作用。轉(zhuǎn)子所受合力為零時(shí)，可認(rèn)為其處于靜平衡位置。其動(dòng)力學(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。

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

(14)

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

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

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

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

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

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

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

3 結(jié)果與討論

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

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

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

3.2 軸承壓力、溫度分布

圖6給出了軸承在供油溫度27 ℃、未發(fā)生磨損工況下瓦塊上表面壓力、溫度分布云圖。由圖6(a)可知，在未發(fā)生磨損的工況下，由于瓦塊的擺動(dòng)，每個(gè)瓦塊均存在高壓區(qū)，產(chǎn)生動(dòng)壓效應(yīng)，從而提高了可傾瓦軸承的穩(wěn)定性。瓦塊上表面高壓區(qū)沿轉(zhuǎn)子旋轉(zhuǎn)方向(逆時(shí)針方向)略有偏移，加載瓦塊(3、4)壓力明顯高于其他瓦塊，起主要承載作用，非加載瓦塊1壓力最低。由圖6(b)可知，在周向方向瓦塊上表面溫度沿轉(zhuǎn)子旋轉(zhuǎn)方向逐漸升高，且在瓦塊后緣處達(dá)到最高，瓦塊間間隙供油，潤滑油在間隙內(nèi)混合使得在進(jìn)入下一油膜間隙時(shí)溫度得到了有效降低；在軸向方向，由于供油位置在瓦塊間隙的中間，瓦塊上表面溫度出現(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 磨損深度對軸承承載力影響

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

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

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

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

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

圖10給出了不同磨損深度下軸心運(yùn)動(dòng)軌跡。

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

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

4 結(jié) 論

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

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

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

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