汽輪機(jī)大尺寸動(dòng)壓軸承特性數(shù)值模擬研究*

(華中科技大學(xué)能源與動(dòng)力工程學(xué)院 湖北武漢 430074)

2016年，我國核電累計(jì)發(fā)電量為2 105.19億kW·h，約占全國累計(jì)發(fā)電量的3.56%，比2015年同期上升了25.07%[1]。按照我國電力行業(yè)發(fā)展規(guī)劃，我國核電裝機(jī)容量還將繼續(xù)增加。隨著核電在我國的加速布局，核電汽輪機(jī)也在不斷向更高容量、更大尺寸和重載荷的方向發(fā)展。如東方電機(jī)股份有限公司某核電機(jī)組軸承直徑已經(jīng)達(dá)到800 mm，比壓達(dá)到2.63 MPa[2]。

雖然用于軸承潤滑性能研究的流體動(dòng)力潤滑理論及軸承計(jì)算的傳統(tǒng)方法已經(jīng)趨于成熟，但是隨著工業(yè)實(shí)踐對(duì)動(dòng)壓軸承結(jié)構(gòu)的更高要求，按照傳統(tǒng)動(dòng)壓理論假設(shè)設(shè)計(jì)的大尺寸動(dòng)壓軸承在運(yùn)行中出現(xiàn)了一些問題，比如碰摩等不穩(wěn)定現(xiàn)象時(shí)有發(fā)生，表明傳統(tǒng)動(dòng)壓理論對(duì)于大尺寸動(dòng)壓軸承的設(shè)計(jì)可能存在某些不適用的地方。

傳統(tǒng)上，國內(nèi)外學(xué)者一般利用求解基于層流假設(shè)的雷諾方程的方法研究油膜軸承動(dòng)靜特性[3-5]。這種方法具有計(jì)算時(shí)間短的突出優(yōu)點(diǎn)，但忽略了慣性項(xiàng)、油膜曲率等因素的影響。而且基于長軸承或者短軸承理論對(duì)滑動(dòng)的求解也存在一定的局限性，比如忽略了油膜軸承的具體局部結(jié)構(gòu)，例如油槽等。隨著計(jì)算流體動(dòng)力學(xué)技術(shù)的發(fā)展和計(jì)算機(jī)性能的提高，通過求解N-S方程的方法研究軸承動(dòng)特性[6-8]已經(jīng)越來越普遍了。采用CFD三維建模，可以考慮任意結(jié)構(gòu)形式軸承及油膜內(nèi)部空化效應(yīng)的影響[9]。通過CFD軟件提供的自定義函數(shù)實(shí)現(xiàn)軸頸的速度和位移擾動(dòng)，便于分析動(dòng)壓軸承剛度阻尼特性[10]。存在的缺點(diǎn)主要是油膜圓周尺寸相對(duì)油膜間隙的大比例增加了高質(zhì)量網(wǎng)格劃分的難度，以及針對(duì)流動(dòng)、空化、兩相流不合理不匹配的模型所導(dǎo)致收斂性差的問題。

目前，針對(duì)核電汽輪機(jī)等行業(yè)的大尺寸動(dòng)壓軸承的研究分析還并不充分，文獻(xiàn)[2]基于經(jīng)典潤滑理論的廣義雷諾方程，結(jié)合流固耦合方法考慮軸瓦的熱彈性變形，分析了直徑為800 mm、1 750 MW級(jí)核電半速發(fā)電機(jī)試驗(yàn)用滑動(dòng)軸承的潤滑特性，軸瓦溫度和油膜壓力的仿真值與實(shí)驗(yàn)值偏差較小，最小油膜厚度和仿真值與測試值偏差15%～22%。文獻(xiàn)[11]指出，在大尺寸油膜軸承的實(shí)際運(yùn)行過程中，軸承發(fā)散區(qū)油楔會(huì)產(chǎn)生負(fù)壓，更易導(dǎo)致潤滑油氣泡析出，發(fā)生空化效應(yīng)并影響壓力場的分布。文獻(xiàn)[12]指出，油膜的空化效應(yīng)將降低承載力并導(dǎo)致軸承結(jié)構(gòu)損傷。

本文作者基于有限體積法，以某核電軸承結(jié)構(gòu)為對(duì)象，考慮空化效應(yīng)和溫黏效應(yīng)，對(duì)典型大尺寸動(dòng)壓軸承流動(dòng)問題進(jìn)行模擬研究。

1 基本方程和計(jì)算方法

1.1 空化理論

動(dòng)壓軸承中的油膜壓力是不均勻的，存在高壓區(qū)和低壓區(qū)。當(dāng)油膜壓力低于一定值時(shí)，潤滑油中溶解的空氣會(huì)析出，該現(xiàn)象稱為空化效應(yīng)。因此在對(duì)流場進(jìn)行求解之前，引入多相流模型。從計(jì)算效率和模型適用范圍的角度考慮，文中采用Mixture模型來模擬潤滑油和空氣之間的兩相耦合作用。

Singhal et al.模型，又稱為全空化模型，是由SINGHAL等[13]于2001年提出。該模型考慮了相變、氣泡動(dòng)力學(xué)、湍流壓力波動(dòng)和不凝結(jié)氣體的影響。兩相連續(xù)方程為

液相：

(1)

氣相：

(2)

混合：

(3)

式中：l表示液相；v表示氣相;ρ表示混合密度，

ρ=αρv+(1-α)ρl

(4)

綜合式(1)—(3)，混合密度(ρ)和氣相體積分?jǐn)?shù)(α)表示為

(5)

在流體和氣泡之間具有零速度滑移的流動(dòng)液體中，氣泡動(dòng)力學(xué)方程源于廣義Rayleigh-Plesset方程：

(6)

式中：RB為氣泡半徑；σ為流體表面張力系數(shù);ρl為流體密度；pB為氣泡表面壓力；p為當(dāng)?shù)剡h(yuǎn)場壓力。

忽略二階項(xiàng)和表面張力，方程可簡化為

(7)

該方程提供了將氣泡動(dòng)力學(xué)的效應(yīng)引入空化模型的物理方法，也是混合物密度的方程。蒸汽體積分?jǐn)?shù)可以用氣泡數(shù)密度(n)和氣泡半徑(RB)表示為

(8)

由式(7)和(8)可得

(9)

由式(7)和(9)，得相變速率R

(10)

根據(jù)氣泡半徑重建相變速率表達(dá)式

(11)

式(11)精確表達(dá)了空化中從液相到氣相(氣泡的生長或氣化)的質(zhì)量傳遞，單位體積的質(zhì)量交換率與氣相密度(ρv)、液相密度(ρl)和混合密度(ρ)相關(guān)。該式通常作為氣泡破裂建模的第一近似值。方程右邊采用絕對(duì)壓力差來修正并被視為源匯相。在不考慮非凝結(jié)氣體、質(zhì)量傳輸和黏性阻尼時(shí)，氣泡壓力(pB)等于飽和蒸汽壓力(pv)。SINGHAL等[13]提出的模型計(jì)入了紊流和不凝結(jié)氣體的影響?？刂品匠虨?/p>

(12)

式中：fv為氣相質(zhì)量分?jǐn)?shù)；fg為不凝結(jié)氣體質(zhì)量分?jǐn)?shù)；Γ為擴(kuò)散系數(shù)。

當(dāng)p≤pv，

(13)

當(dāng)p>pv，

(14)

通過對(duì)湍流壓力波動(dòng)的局部值的估計(jì)來校正飽和壓力

(15)

式中：Fvap=0.02，F(xiàn)cond=0.01。

1.2 剛度阻尼計(jì)算方法

油膜的剛度阻尼是動(dòng)壓軸承的重要特性參數(shù)。轉(zhuǎn)子在平衡位置受外載荷作用，對(duì)油膜產(chǎn)生位移擾動(dòng)和速度擾動(dòng)，軸承油膜力將發(fā)生變化。由攝動(dòng)理論，將油膜力線性化處理后，油膜力的變化與轉(zhuǎn)子擾動(dòng)位移和擾動(dòng)速度之間的關(guān)系如下：

(16)

(17)

采用動(dòng)網(wǎng)格技術(shù)，令轉(zhuǎn)子分別在水平方向和垂直方向移動(dòng)微小位移ΔX、ΔY，并計(jì)算出有無擾動(dòng)位移前后油膜在水平和垂直方向的油膜力的變化ΔFX、ΔFY，通過式(16)、(17)可以求得KXX、KYX、KXY和KYY。

(18)

(19)

2 計(jì)算域網(wǎng)格與邊界條件

文中分析的油膜軸承是核電汽輪機(jī)動(dòng)壓軸承，相比抽油機(jī)、小型油膜軸承，具有更大的尺寸和相對(duì)更小的半徑間隙。軸承直徑D=800 mm，寬度B=680 mm，相對(duì)頂隙比ψ=0.15%，載荷為F=1 374.65 kN。入口直徑d=160 mm，方形槽寬度為軸承寬度的0.7倍。潤滑油為20#汽輪機(jī)油，分析中計(jì)入溫黏效應(yīng)。入口油壓為0.2 MPa，出口油壓為0.1 MPa?？栈Ｐ椭锌栈瘔毫Σ捎梦墨I(xiàn)[14]的推薦值29 158 Pa。

油膜流體域網(wǎng)格在ICEM CFD中采用六面體結(jié)構(gòu)化網(wǎng)格劃分，進(jìn)口處采用O-BLOCK結(jié)構(gòu)劃分，采用Determinant 2×2×2標(biāo)準(zhǔn)檢查全局網(wǎng)格質(zhì)量大于0.6。采取50萬、100萬、200萬、400萬4種網(wǎng)格數(shù)量方案，如圖1和圖2所示為最大壓力和承載力隨網(wǎng)格數(shù)量變化曲線，可知當(dāng)網(wǎng)格數(shù)量大于200萬后，油膜最大壓力以及油膜承載力變化較小，已經(jīng)滿足計(jì)算精度的要求。兼顧計(jì)算精度、計(jì)算時(shí)間，后續(xù)分析時(shí)網(wǎng)格數(shù)量統(tǒng)一取為200萬。

圖1 最大壓-網(wǎng)格數(shù)量曲線

圖2 承載力-網(wǎng)格數(shù)量曲線

3 結(jié)果與分析

3.1 流動(dòng)模型和空化模型的比較討論

表1 各轉(zhuǎn)速工況下雷諾數(shù)

表2示出了軸承在轉(zhuǎn)速為1 500 r/min、偏心率為0.7時(shí)，分別在層流和紊流realizable k-epsilon模型下的靜態(tài)參數(shù)。可見，層流模型計(jì)算的數(shù)值與經(jīng)驗(yàn)值偏離較遠(yuǎn)，紊流模型的結(jié)果與經(jīng)驗(yàn)數(shù)值更為接近。同時(shí)，采用紊流模型計(jì)算的最大溫度值、壓力值與文獻(xiàn)[2]同直徑的大尺寸軸承試驗(yàn)得到的最大溫度值82.8 ℃、壓力值6.65 MPa也較為接近。因此，采用湍流模型比較適合所研究的軸承的油膜流動(dòng)狀態(tài)。

表2 不同流動(dòng)模型分析結(jié)果

空化現(xiàn)象常存在于軸承油膜中，常用的數(shù)值空化模型有3種：Singhal et al.模型，Zwart-Gerber-Belamri模型和Schnerr and Sauer模型。文獻(xiàn)[17]論述了各個(gè)空化模型的特點(diǎn)，對(duì)于旋轉(zhuǎn)機(jī)械的空化現(xiàn)象模擬，該文獻(xiàn)建議使用后2種空化模型。而在文中的計(jì)算分

析中，通過調(diào)整初始化方式和松弛因子，Zwart-Gerber-Belamri模型和Schnerr and Sauer模型收斂異常困難，epsilon和氣相的殘差很難降低至10-2量級(jí)，而Singhal et al.模型收斂精度高，各殘差均在10-5量級(jí)以下。因此文中采用Singhal et al.模型模擬油膜中空化效應(yīng)。

3.2 基于力平衡的油膜參數(shù)求解

為了得到額定轉(zhuǎn)速1 500 r/min及給定外載荷1 374.65 kN下的軸承最小油膜厚度及相關(guān)參數(shù)，建立了不同偏心量的油膜軸承模型并進(jìn)行試算，如表3和圖3所示。當(dāng)某偏心率在該轉(zhuǎn)速下剛好滿足外荷載等于承載力時(shí)，此時(shí)油膜偏心率即為所求偏心率。其他靜態(tài)特性及相關(guān)參數(shù)在此偏心基礎(chǔ)上進(jìn)行計(jì)算求取。

表3 不同偏心下承載力計(jì)算

圖3 承載力-偏心率曲線

試算結(jié)果表明，當(dāng)偏心率為0.69時(shí)軸承的承載力為1 365.00 kN，與給定外載荷僅偏差約0.73%。此時(shí)計(jì)算的最小油膜厚度約為186 μm。圖4所示為此時(shí)軸頸壓力分布，最大壓力約為6.35 MPa。圖5所示為流體域中速度場的分布，最大線速度為62.74 m/s，入口方形槽中速度較低但存在旋流引起的渦。

圖4 軸頸側(cè)油膜壓力分布

圖5 速度場分布

圖6所示為軸瓦溫度分布云圖，最大溫度為81.52 ℃，在接近入口油槽處。圖7所示為軸頸和軸瓦側(cè)油膜的氣穴體積分?jǐn)?shù)云圖，其中軸頸表面的氣穴體積分?jǐn)?shù)更高，空化效應(yīng)較軸瓦側(cè)更加明顯。氣穴體積分?jǐn)?shù)最大處在接近入口油槽處，最大約為0.81。

圖6 軸瓦側(cè)溫度分布云圖

圖7 軸瓦側(cè)(左)和軸頸側(cè)(右)氣穴體積分?jǐn)?shù)

3.3 轉(zhuǎn)速波動(dòng)下的軸承靜特性

工程應(yīng)用中，除了考察額定轉(zhuǎn)速下的靜特性外，也關(guān)注其他各重要轉(zhuǎn)速如1 000、1 200、1 750 r/min時(shí)軸承的最小油膜厚度、溫度、壓力偏位角等靜態(tài)特性參數(shù)。同樣采用基于承載力和外載荷平衡的原則進(jìn)行多轉(zhuǎn)速工況下，多偏心尺寸的油膜軸承試算。表4給出了各轉(zhuǎn)速下滿足外載荷要求的最小油膜厚度值。由圖8可知，在該轉(zhuǎn)速范圍內(nèi)，轉(zhuǎn)子轉(zhuǎn)速越高軸承最小油膜厚度越大。轉(zhuǎn)速過低時(shí)，如1 000 r/min時(shí)，最小油膜厚度僅為132 μm。

表4 定載荷各轉(zhuǎn)速下的最小油膜厚度

圖8 不同轉(zhuǎn)速下的最小油膜厚度曲線

Fig 8 Minimum oil film thickness curve at different rotation speeds

圖9示出了1 000～1 750 r/min轉(zhuǎn)速范圍內(nèi)偏位角的變化趨勢，轉(zhuǎn)子轉(zhuǎn)速越高軸承的偏位角越小，偏位角范圍為40°～55°。圖10示出了各轉(zhuǎn)速下油膜中心周向的壓力分布，轉(zhuǎn)子轉(zhuǎn)速越高，壓力峰值越低，壓力峰值的位置出現(xiàn)得越早。在各轉(zhuǎn)速下壓力峰值的范圍為6～8 MPa。圖11示出了各轉(zhuǎn)速下油膜中心周向的溫度分布，溫度趨勢均隨角度增大而增大，在接近入口方形槽溫度急速降低，最大溫度在接近入口油槽處。

圖9 不同轉(zhuǎn)速下的偏位角曲線

圖10 軸頸中心壓力曲線

圖11 軸瓦中心溫度分布曲線

3.4 油膜動(dòng)態(tài)特性求解

圖12示出了各轉(zhuǎn)速下油膜的剛度系數(shù)，圖13示出了各轉(zhuǎn)速下的阻尼系數(shù)。因計(jì)算剛度和阻尼采用相同的擾動(dòng)速度和擾動(dòng)位移，因此剛度和阻尼在隨轉(zhuǎn)速的變化趨勢上相同。在1 000～1 750 r/min轉(zhuǎn)速范圍內(nèi)，剛度系數(shù)的量級(jí)為109～1010，阻尼系數(shù)的量級(jí)為107～108。

圖12 各轉(zhuǎn)速下油膜剛度系數(shù)

圖13 各轉(zhuǎn)速下油膜阻尼系數(shù)

4 結(jié)論

(1)針對(duì)大尺寸動(dòng)壓軸承油膜流場，可根據(jù)臨界雷諾數(shù)準(zhǔn)則判斷流動(dòng)模態(tài)，而不應(yīng)直接采用動(dòng)壓軸承設(shè)計(jì)手冊(cè)的層流假設(shè)。針對(duì)文中研究對(duì)象在1 000、1 250、1 500、1 750 r/min時(shí)采用realizable k-epsilon模型計(jì)算更為準(zhǔn)確。

(2)針對(duì)文中研究對(duì)象，在低于額定轉(zhuǎn)速1 500 r/min以下工況時(shí)，潤滑油溫度由常溫到較高溫度轉(zhuǎn)變時(shí)，流動(dòng)模態(tài)有可能發(fā)生轉(zhuǎn)換，此時(shí)軸承油膜的承載力和其他參數(shù)均有可能發(fā)生改變，易于導(dǎo)致軸承運(yùn)行不穩(wěn)，值得引起關(guān)注。

(3)在對(duì)動(dòng)壓軸承進(jìn)行數(shù)值空化模擬收斂困難時(shí)，空化模擬中可采用Singhal et al.模型，該模型擁有更高的收斂精度和更快的收斂速度。

(4)計(jì)算得到了典型尺寸軸承在不同轉(zhuǎn)速工況條件下油膜厚度規(guī)律、壓力分布、溫度分布、滑油流量等靜態(tài)特性和軸承剛度、阻尼動(dòng)態(tài)特性等數(shù)據(jù)，可為類似大型動(dòng)壓油膜軸承的設(shè)計(jì)和軸承的載荷分配提供參考。