基于CFD的透平膨脹機軸向推力計算

李 超 張國軍 程金明 孫照嵐

(蘭州理工大學石油化工學院)

基于CFD的透平膨脹機軸向推力計算

李 超*張國軍 程金明 孫照嵐

(蘭州理工大學石油化工學院)

為了準確地預測透平膨脹機的軸向推力，建立其力學分析模型和葉輪間隙、迷宮密封三維模型。通過受力分析，采用低Re數k-ε模型和Simplec算法，考慮機械磨損的影響，對透平膨脹機軸向推力進行了模擬計算，并將模擬結果與軸向推力經驗公式計算結果進行比較。

透平膨脹機 軸向推力 CFD 工作輪 增壓輪

透平膨脹機是利用工質的絕熱等熵膨脹流動來獲得低溫和進行能量回收的主要設備[1]，在空調、空氣分離、低溫環境以及廢棄余能等的綜合利用中有著顯著的優勢[2]。隨著透平膨脹機的大型化和高速化，對其機組設計、制造、運行壽命及可靠性評估等方面提出了更高的要求，尤其是如何提高轉子運行穩定性方面。

葉輪作為透平膨脹機的核心部件，其性能和結構的穩定性對整個透平膨脹機安全可靠的運行起著至關重要的作用。但由于在設計階段對轉子軸向推力的計算不準確，導致葉輪軸向竄動，從而影響透平膨脹機的安全、可靠運轉，因此提高軸向推力計算精度、軸系運轉穩定性和機組可靠性，對透平膨脹機理論設計和延長設備使用壽命具有重要理論意義和工程應用價值。

軸向推力的計算主要有經驗公式法、解析法和數值解法3種方法。高原分析了透平膨脹機軸向推力的組成，介紹了透平膨脹機軸向推力估算的有關近似計算公式，并對增壓透平膨脹機轉子軸向推力進行了計算[3,4]。Kurokawa J等通過分析火箭發動機中透平泵旋轉壁面和靜止壁面之間的間隙流流動特性，進而用解析法計算其軸向推力[5]。

近年來，隨著CFD技術的快速發展，越來越多的研究者采用數值法計算軸向推力。趙問銀等對大型內壓縮流程空分裝置用低溫液體膨脹機在不同工況下的整級流場進行了數值模擬計算，預測了不同工況下轉子軸向推力水平[6,7]。王維民等分析了離心壓縮機軸向推力的組成，并采用數值法對其進行計算[8,9]。Han Z X和Cizmas P G A在考慮離心壓縮機葉輪間隙與葉輪耦合的基礎上計算出離心壓縮機的軸向推力[10]。

以上文獻均對軸向推力進行了計算，但未考慮到透平膨脹機在運轉時存在的機械磨損。機械磨損使密封齒齒高減小、密封間隙增大，導致葉輪空腔側壓力分布不同，進而使軸向推力發生變化。筆者分析了某透平膨脹機的軸向推力，通過建立其工作輪、增壓輪間隙和迷宮密封的三維模型，在考慮機械磨損等因素的基礎上，采用CFD對葉輪兩側間隙內壓力和軸向推力進行了分析計算。

1 透平膨脹機葉輪兩側間隙中流體流動CFD分析

1.1透平膨脹機軸向推力組成

整個透平膨脹機轉子所受的軸向推力是增壓葉輪所受的軸向推力與工作輪所受的軸向推力之和。工作輪與增壓輪結構尺寸和軸向推力分布如圖1所示。

a. 工作輪

b. 增壓輪

工作輪軸向推力的計算式為：

Fe=F1e+F2e-F3e-F4e-F5e

(1)

式中F1e、F3e——工作輪輪背側和輪蓋側間隙內泄漏氣體所產生的軸向推力；

F2e——作用在工作輪輪背側環形端面上的軸向推力；

F4e——作用在工作輪出口端面上的軸向推力；

F5e——氣體排出時由徑向流動變為軸向流動引起動量變化所產生的軸向推力。

增壓輪軸向推力的計算公式為：

Fc=F1c+F2c+F3c-F4c-F5c

(2)

式中F1c、F3c——增壓輪輪背側和輪蓋側間隙內泄漏氣體所產生的軸向推力；

F2c——作用在增壓輪輪背側環形端面上的軸向推力；

F4c——增壓輪進口壓力作用于葉輪蓋面和葉輪內盤面所產生的軸向推力；

F5c——增壓輪進口氣體由軸向流動變為徑向流動引起動量變化所產生的軸向推力。

則透平膨脹機轉子所受的總軸向推力Fa為：

Fa=Fc-Fe

(3)

1.2CFD分析模型

葉輪轉速、間隙結構尺寸、間隙內氣體壓力及葉輪出口處氣體壓力等是影響葉輪兩側間隙內壓力分布的主要因素。分別建立工作輪和增壓輪輪蓋側、輪背側間隙和對應的密封模型，以便精確地確定兩側密封的泄漏量和軸向推力。工作輪和增壓輪的結構參數如下：

工作輪輪軸線至工作輪出口端面的距離d1107mm

工作輪外徑d2160mm

轉軸直徑ds130mm

工作輪軸線至增壓輪背側出口端面的距離ds2130mm

工作輪密封間隙 0.12mm

增壓輪軸線至增壓輪進口端面的距離D1113mm

增壓輪外徑D2174mm

轉軸直徑Ds130mm

增壓輪軸線至增壓輪背側出口端面的距離Ds2130mm

增壓輪密封間隙 0.12mm

工作輪和增壓輪的工質均為氮氣，轉速為31 328r/min，但工作輪的工質進口壓力為1.714MPa，出口壓力為0.628MPa；增壓輪的工質進口壓力為2.592MPa，出口壓力為4.100MPa。

透平膨脹機葉輪兩側間隙的幾何結構相對于葉輪轉軸是對稱的，因此將間隙內的流場看作周期對稱結構。采用六面體網格來對計算域進行劃分(圖2、3)，在近壁面處存在較大的速度梯度，采用較細密的網格劃分，網格密度在遠離壁面的方向逐漸變疏。近壁面處相鄰網格的邊長比率設定為1.1。由于在密封出口處具有強大的渦流，為了提高計算精度，筆者增加了下游區域。

a. 輪蓋側

b. 輪背側

a. 輪蓋側

b. 輪背側

圖3 增壓輪網格截面

1.3控制方程和邊界條件

對于定常流動，王福軍給出了其控制方程的通用形式[11]：

(4)

式中S——廣義源項；

φ——通用變量，可代表1，u、v、w及T等求解變量；

ρ——流體密度；

?！獜V義擴散系數。

由于標準k-ε、RNGk-ε不適合描述Re數較低的縫隙內流體流動，因此筆者采用低Re數k-ε模型描述縫隙內的紊流，其輸運方程為[11]：

(5)

(6)

(7)

f1≈1.0

fμ=exp[-2.5/(1+Ret/50)]

Ret=ρk2/(ηε)

式中Gk——由平均流速梯度引起的紊動能k的產生項；

k——紊動能；

n——壁面法向坐標；

Ret——湍流雷諾數；

u——與壁面平行的流速；

ε——紊動能耗散率；

μt——湍動粘度。

模型常數C1ε、C2ε、Cμ、σk、σε的取值分別為1.44、1.92、0.09、1.0、1.3。

由葉輪輪蓋、輪蓋密封、輪盤、輪盤密封和隔板組成的縫隙入口邊界條件為入口壓力邊界條件，出口邊界條件為壓力出口邊界條件。輪蓋側隔板和密封是靜止壁面，輪蓋為轉動壁面；輪背側輪盤和密封為轉動壁面；輪背側隔板為靜止壁面。壓力速度耦合方式采用Simplec算法，欠松弛因子取默認值，采用二階迎風格式離散差分方程，收斂精度設定為10-4。

2 數值模擬結果

圖4為工作輪和增壓輪輪蓋側、輪背側壓力隨葉輪半徑的變化規律。當氣體由泄漏通道的水平方向轉彎進入徑向方向流動時，葉輪側為降壓流動，固壁側為擴壓流動；轉彎后葉輪側為擴壓流動，固壁側為降壓流動，因此壓力呈現出先降低再升高的變化趨勢。隨著向心流動，由于摩擦損耗，壓力沿流動方向逐漸減小，在密封處由于轉彎和流通面積減小，壓力變化呈現為略有上升再下降的趨勢。圖5所示為泄漏通道入口轉彎處同一直徑處氣體壓力沿軸向的分布曲線，從圖5可知，由于彎道的影響，氣體壓力沿軸向從葉輪側向固壁側逐漸增加。

a. 工作輪輪蓋側

b. 工作輪輪背側

c. 增壓輪輪蓋側

d. 增壓輪輪背側

a. 工作輪

b. 增壓輪

當流體在壓差的驅動下進入密封腔內，由于面積突然的擴大，流體的動能轉化為熱能，形成一個強烈的反轉漩渦。在密封齒尖存在噴射流動，此時流體具有較高的速度。由于采用直通型梳齒密封，部分流體不經過旋流而直接進入下一密封腔，呈現出密封直吹效應。密封齒結構型式不同，氣體流動形態和泄漏量也有所不同。

為了驗證計算結果與網格無關，通過改變網格密度，對工作輪輪背側計算區域采用3種不同網格密度進行計算，計算結果見表1，根據表1可以判斷出該計算結果是收斂且具有足夠精度的[10]。

表1 網格數量對計算精度的影響

由于存在機械磨損，隨著機器的運行，迷宮密封的密封齒高度將逐漸減小，使密封齒與轉子之間的間隙增加，導致泄漏量增加，從而影響透平膨脹機的效率和軸向推力。設定密封間隙分別為0.12、0.24、0.36、0.48mm，計算4組不同間隙下葉輪的泄漏量和軸向推力(圖6，表2)。從圖6可以看出：隨著密封間隙的增大泄漏量也逐漸增大，嚴重影響著透平膨脹機的效率；從表2可以看出：隨著間隙的增大，工作輪輪背側的軸向推力變化較大，而輪蓋側的軸向推力變化相對較小。

圖6 不同密封間隙下葉輪間隙內的泄漏量

間隙/mmF1e/NF3e/NF1c/NF3c/N0.1210996.018478.841976.055400.40.2410839.618399.641835.655314.00.3610623.618356.441515.255242.00.4810312.218345.640996.855177.2

3 軸向推力的計算與比較

將數值模擬結果與采用經驗公式法計算的葉輪氣動軸向推力計算結果進行比較。根據式(1)、(2)和文獻[3]可得經驗公式為：

(8)

(9)

式中c1——氣體出工作輪時的軸向速度；

c2——氣體進增壓輪時的軸向速度；

p1e、p2e——工作輪氣體的入口和出口壓力；

p1c、p2c——增壓輪氣體的入口和出口壓力；

qm1——出工作輪的氣體質量流量；

qm2——進增壓輪的氣體質量流量。

計算結果見表3。

表3 工作輪氣動軸向推力計算結果 N

從表3可以看出：兩種計算方法得出的轉子葉輪側的軸向推力相差較大。氣體由泄漏通道入口的水平方向轉彎進入徑向方向、氣體在密封處轉彎和梳齒密封使流通面積減小，進而導致葉輪側壓力分布發生變化，這是造成這種差異的主要原因。

4 結論

4.1軸向推力的大小受葉輪兩側間隙通道轉彎和梳齒密封的影響而發生變化，隨齒縫間隙的增大，氣體泄漏量增大，軸向推力減小。

4.2泄漏通道入口轉彎處同一直徑的氣體壓力沿軸向由葉輪側向固壁側逐漸增大；梳齒密封轉彎處同一直徑的氣體壓力沿軸向由葉輪側向固壁側逐漸減小。

[1] 計光華.透平膨脹機[M].西安：西安交通大學出版社，2006:1～9.

[2] 李公法,石教華,王瑾輝,等.基于ECS-700的透平膨脹機控制系統[J].化工自動化及儀表，2011，38(5)：559～561.

[3] 高原.氣體分離用透平機械[M].北京：石油工業出版社，1991:53～59.

[4] 高原.增壓透平膨脹機變工況轉子軸向力的估算[J].低溫工程，1998，(4):37～41.

[5] Kurokawa J, Kamijo K, Shimura T. Axial Thrust Behavior in LOX-pump of Rocket Engine[J]. Journal of Propulsion and Power, 1994, 10(2):244～250.

[6] 趙問銀，王科，宋鵬.大型空分裝置用低溫液體膨脹機內流及軸向推力數值計算[J].低溫工程，2012，(2)：16～21.

[7] Wang K，Sun J J, He Z L, et al. Prediction of Axial Thrust Load Acting on a Cryogenic Liquid Turbine Impeller [C].ASME 2011 Turbo Expo: Turbine Technical Conference and Exposition.New York: ASME, 2011, 487～498.

[8] 王維民，高金吉，江志農.高壓離心式壓縮機軸向推力研究[J].風機技術，2006，(6):1～4.

[9] 聞蘇平，朱報楨，苗永淼.高壓離心壓縮機軸向推力計算[J].西安交通大學學報，1998，32(11)：63～67.

[10] Han Z X, Cizmas P G A. A CFD Method for Axial Thrust Load Prediction of Centrifugal Compressor[J].International Journal of Turbo and Jet-Engines, 2003, 20(1):1～16.

[11] 王福軍.計算流體動力學分析——CFD軟件原理與應用[M].北京：清華大學出版社，2004：126～132.

EstimationofAxialThrustinTurboExpanderBasedonCFD

LI Chao, ZHANG Guo-jun, CHENG Jin-ming, SUN Zhao-lan

(CollegeofPetrochemicalIndustry,LanzhouUniversityofTechnology,Lanzhou730050,China)

To predict the turbo-expander’s axial thrust accurately, a turbo-expander was chosen as an example, and the mechanical model and the three-dimensional model of the impeller clearance and labyrinth seal were built. Through analyzing the forces and using low Reynolds numberk-εturbulence model and the SIMPLEC algorithm, as well as considering the impact of leakage, the axial thrust acting on the impeller was simulated, and the simulation results and the results calculated by the empirical method were compared.

turbo-expander, axial thrust, CFD, impeller, supercharger vane

*李 超，男，1958年3月生，教授。甘肅省蘭州市，730050。

TQ051.5

A

0254-6094(2015)01-0110-06

2014-03-06，

2015-01-06)