汽輪機大尺寸動壓軸承特性數值模擬研究*

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

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

雖然用于軸承潤滑性能研究的流體動力潤滑理論及軸承計算的傳統方法已經趨于成熟，但是隨著工業實踐對動壓軸承結構的更高要求，按照傳統動壓理論假設設計的大尺寸動壓軸承在運行中出現了一些問題，比如碰摩等不穩定現象時有發生，表明傳統動壓理論對于大尺寸動壓軸承的設計可能存在某些不適用的地方。

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

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

本文作者基于有限體積法，以某核電軸承結構為對象，考慮空化效應和溫黏效應，對典型大尺寸動壓軸承流動問題進行模擬研究。

1 基本方程和計算方法

1.1 空化理論

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

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

液相：

(1)

氣相：

(2)

混合：

(3)

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

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

(4)

綜合式(1)—(3)，混合密度(ρ)和氣相體積分數(α)表示為

(5)

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

(6)

式中：RB為氣泡半徑；σ為流體表面張力系數;ρl為流體密度；pB為氣泡表面壓力；p為當地遠場壓力。

忽略二階項和表面張力，方程可簡化為

(7)

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

(8)

由式(7)和(8)可得

(9)

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

(10)

根據氣泡半徑重建相變速率表達式

(11)

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

(12)

式中：fv為氣相質量分數；fg為不凝結氣體質量分數；Γ為擴散系數。

當p≤pv，

(13)

當p>pv，

(14)

通過對湍流壓力波動的局部值的估計來校正飽和壓力

(15)

式中：Fvap=0.02，Fcond=0.01。

1.2 剛度阻尼計算方法

油膜的剛度阻尼是動壓軸承的重要特性參數。轉子在平衡位置受外載荷作用，對油膜產生位移擾動和速度擾動，軸承油膜力將發生變化。由攝動理論，將油膜力線性化處理后，油膜力的變化與轉子擾動位移和擾動速度之間的關系如下：

(16)

(17)

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

(18)

(19)

2 計算域網格與邊界條件

文中分析的油膜軸承是核電汽輪機動壓軸承，相比抽油機、小型油膜軸承，具有更大的尺寸和相對更小的半徑間隙。軸承直徑D=800 mm，寬度B=680 mm，相對頂隙比ψ=0.15%，載荷為F=1 374.65 kN。入口直徑d=160 mm，方形槽寬度為軸承寬度的0.7倍。潤滑油為20#汽輪機油，分析中計入溫黏效應。入口油壓為0.2 MPa，出口油壓為0.1 MPa。空化模型中空化壓力采用文獻[14]的推薦值29 158 Pa。

油膜流體域網格在ICEM CFD中采用六面體結構化網格劃分，進口處采用O-BLOCK結構劃分，采用Determinant 2×2×2標準檢查全局網格質量大于0.6。采取50萬、100萬、200萬、400萬4種網格數量方案，如圖1和圖2所示為最大壓力和承載力隨網格數量變化曲線，可知當網格數量大于200萬后，油膜最大壓力以及油膜承載力變化較小，已經滿足計算精度的要求。兼顧計算精度、計算時間，后續分析時網格數量統一取為200萬。

圖1 最大壓-網格數量曲線

圖2 承載力-網格數量曲線

3 結果與分析

3.1 流動模型和空化模型的比較討論

表1 各轉速工況下雷諾數

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

表2 不同流動模型分析結果

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

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

3.2 基于力平衡的油膜參數求解

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

表3 不同偏心下承載力計算

圖3 承載力-偏心率曲線

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

圖4 軸頸側油膜壓力分布

圖5 速度場分布

圖6所示為軸瓦溫度分布云圖，最大溫度為81.52 ℃，在接近入口油槽處。圖7所示為軸頸和軸瓦側油膜的氣穴體積分數云圖，其中軸頸表面的氣穴體積分數更高，空化效應較軸瓦側更加明顯。氣穴體積分數最大處在接近入口油槽處，最大約為0.81。

圖6 軸瓦側溫度分布云圖

圖7 軸瓦側(左)和軸頸側(右)氣穴體積分數

3.3 轉速波動下的軸承靜特性

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

表4 定載荷各轉速下的最小油膜厚度

圖8 不同轉速下的最小油膜厚度曲線

Fig 8 Minimum oil film thickness curve at different rotation speeds

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

圖9 不同轉速下的偏位角曲線

圖10 軸頸中心壓力曲線

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

3.4 油膜動態特性求解

圖12示出了各轉速下油膜的剛度系數，圖13示出了各轉速下的阻尼系數。因計算剛度和阻尼采用相同的擾動速度和擾動位移，因此剛度和阻尼在隨轉速的變化趨勢上相同。在1 000～1 750 r/min轉速范圍內，剛度系數的量級為109～1010，阻尼系數的量級為107～108。

圖12 各轉速下油膜剛度系數

圖13 各轉速下油膜阻尼系數

4 結論

(1)針對大尺寸動壓軸承油膜流場，可根據臨界雷諾數準則判斷流動模態，而不應直接采用動壓軸承設計手冊的層流假設。針對文中研究對象在1 000、1 250、1 500、1 750 r/min時采用realizable k-epsilon模型計算更為準確。

(2)針對文中研究對象，在低于額定轉速1 500 r/min以下工況時，潤滑油溫度由常溫到較高溫度轉變時，流動模態有可能發生轉換，此時軸承油膜的承載力和其他參數均有可能發生改變，易于導致軸承運行不穩，值得引起關注。

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

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