顧乾磊， 張萬福， 張 堯， 陳璐琪， 李 春， 楊建剛

(1. 上海理工大學(xué) 能源與動力工程學(xué)院,上海 200093；2. 東南大學(xué) 火電機組振動國家工程研究中心,南京 210096)

密封是汽輪機、燃?xì)廨啓C、壓縮機等大型旋轉(zhuǎn)機械中抑制流體泄漏的關(guān)鍵部件，直接影響機組的高效運行[1]。然而，隨著機組容量及工質(zhì)參數(shù)的不斷提高，密封引起的氣流激振問題越來越突出，對高參數(shù)大型透平機械安全與穩(wěn)定運行形成極大挑戰(zhàn)[2-5]。近年，對于密封的研究重點也因此由傳統(tǒng)的靜力特性(如泄漏性能)逐漸轉(zhuǎn)向動力特性。與滑動軸承類似，通常用8個或12個動力特性系數(shù)(剛度、阻尼及質(zhì)量慣性項)來描述密封的動力特性，但密封間隙通常是軸承間隙的2～10倍，密封內(nèi)部流動為更加復(fù)雜的三維可壓縮湍流。因此，密封動力特性系數(shù)的高精度識別及氣流激振的準(zhǔn)確定位是目前國內(nèi)外研究的熱點與難點。

目前對于密封動力特性系數(shù)的識別方法主要有理論、試驗及計算流體力學(xué)(Computational Fluid Dynamics, CFD)。密封理論計算在早期的Bulk-Flow模型[6-7]之后，主要以控制體模型為主，學(xué)者先后提出單控制體模型[8-9]、雙控制體模型[10-13]、三控制體模型[14-15]，取得了一定的預(yù)測效果，然而控制體模型在處理偏心、腔室周向速度、能量損耗系數(shù)、邊界條件等方面做了大量假設(shè)，特別是隨先進密封幾何結(jié)構(gòu)復(fù)雜性與工質(zhì)參數(shù)的提高，以及阻塞流動帶來的較差斂散性，限制了其應(yīng)用的廣度和深度。試驗是密封動靜特性研究的重要手段[16-17]，主要方法有動態(tài)壓力測試法和機械阻抗法，后者由于具有較好的信噪比，在實際中被廣泛應(yīng)用。Texas A & M大學(xué)的Childs教授和Vance教授等在實驗方面做了大量研究工作，從20世紀(jì)80年代開始至今一直從事密封(如梳齒、蜂窩、孔型、刷式、袋式等)動靜特性的試驗識別。1986年第一次給出了密封直接阻尼系數(shù)(最高壓力8.25 bar、最高轉(zhuǎn)速8 000 r/min的試驗值[18]；1988年對試驗臺進行了改進(最高轉(zhuǎn)速16 000 r/min)，并與控制體模型進行了對比，發(fā)現(xiàn)雙控制體模型的預(yù)測結(jié)果要優(yōu)于單控制體模型[19]；2002年Dawson等[20-21]在已有軸承試驗臺基礎(chǔ)上對其改進，密封試驗最高壓力達17.2 bar，2003年試驗參數(shù)進一步提高至最高壓力84 bar、最高轉(zhuǎn)速29 000 r/min。Tiwari[22]對密封試驗的研究進展做了很好的綜述。國內(nèi)關(guān)于密封動力特性的試驗研究相對較少，楊建剛等[23-24]提出密封動力特性系數(shù)識別的影響系數(shù)法和基于不平衡同頻激勵的密封動力特性系數(shù)識別方法。然而，開展試驗研究，特別是不同類型、多種尺寸與工況下的密封性能測試需要大量時間和人力物力投入，代價較高，而且試驗在流場細(xì)節(jié)等方面的測試較為困難。

隨著計算機硬件的高速發(fā)展及其計算能力的不斷提高，計算流體力學(xué)方法正在受到越來越多的應(yīng)用。從早期基于有限差分、有限體積等方法的小規(guī)模編程計算[25-26]，到目前以大型商業(yè)計算軟件(如FLUENT、CFX、STAR CD等)為主的CFD三維湍流模擬。密封動力特性系數(shù)的CFD識別方法主要有如下三種：

(1)有限小擾動法。該方法與滑動軸承動力特性系數(shù)的識別類似，通過在平衡點附近給轉(zhuǎn)子正交方向上較小的位移擾動和速度擾動，進而得到剛度和阻尼系數(shù)，如式(1)所示

(1)

(2)旋轉(zhuǎn)坐標(biāo)系法[27-28]。如圖1所示，該方法應(yīng)用旋轉(zhuǎn)坐標(biāo)系將轉(zhuǎn)子實際的非穩(wěn)態(tài)渦動問題轉(zhuǎn)變?yōu)闇?zhǔn)穩(wěn)態(tài)模型，計算速度較快。

瞬態(tài)渦動法。該方法基于CFD非穩(wěn)態(tài)求解及動網(wǎng)格方法，對密封內(nèi)流場進行計算。西安交通大學(xué)李軍教授團隊做了大量的密封數(shù)值模擬工作，提出了多頻渦動瞬態(tài)計算方法[29-31]，并針對蜂窩密封、孔型密封、袋式密封等開展了系統(tǒng)的數(shù)值仿真研究。

圖1 準(zhǔn)穩(wěn)態(tài)模型示意圖Fig.1 Schematic diagram of quasi-steady state model

有限小擾動法和旋轉(zhuǎn)坐標(biāo)系方法使用較為方便，然而在實際應(yīng)用過程中局限性較多，如：① 轉(zhuǎn)子實際渦動軌跡較為復(fù)雜，如存在偏心、橢圓渦動等；② 密封動力特性系數(shù)與偏心率、渦動頻率、轉(zhuǎn)速等都存在依賴關(guān)系；③ 復(fù)雜密封壁面(蜂窩、孔型、袋式等)的模擬較難實現(xiàn)。瞬態(tài)渦動法避免了上述問題，通過在不同方向施加激振計算密封氣流力及密封動力特性系數(shù)，然而針對轉(zhuǎn)子實際任意軌跡的計算還需要更深入地研究。

本文提出基于微元軌跡的密封動力特性系數(shù)理論識別方法，建立任一頻率激勵下流體動力學(xué)識別模型，應(yīng)用CFD瞬態(tài)求解及動網(wǎng)格方法得到密封氣流力特性及動力特性系數(shù)，并與試驗值進行對比，驗證模型與計算方法的準(zhǔn)確性。

1 密封動力特性系數(shù)識別方法

圖2 氣缸-密封系統(tǒng)模型Fig.2 Cylinder-seal system model

圖2中：1為轉(zhuǎn)子渦動軌跡；2為渦動轉(zhuǎn)子；3為靜子(氣缸)；4為靜平衡位置的轉(zhuǎn)子；Fη,Fξ為轉(zhuǎn)子在η和ξ方向所受氣流力；Fe,Fα為轉(zhuǎn)子在e和α方向所受氣流力。

轉(zhuǎn)子做橢圓軌跡渦動，渦動轉(zhuǎn)速為Ωi(i=1, 2, 3, …,n)，(X,Y)坐標(biāo)系下的渦動軌跡為

(2)

式中:θ為橢圓軌跡傾斜角θ∈(0°,360°)；a和b分別為橢圓軌跡的長、短半軸長度；X0和Y0分別為渦動軌跡中心在(X,Y)坐標(biāo)系下的橫縱坐標(biāo)。

在(e,α)坐標(biāo)系中，橢圓軌跡方程為

(3)

相應(yīng)的轉(zhuǎn)子渦動速度為

(4)

以O(shè)1點為原點建立坐標(biāo)系，忽略氣體質(zhì)量慣性力影響，轉(zhuǎn)子小渦動軌跡下密封動力學(xué)模型可線性簡化為

(5)

將式(3)與式(4)代入式(5)得

(6)

應(yīng)用瞬態(tài)法得到任意渦動轉(zhuǎn)速Ωi(i=1, 2, 3, 4, …,n)轉(zhuǎn)子所受氣流力，在t=0和t=T/4時刻，轉(zhuǎn)子在e方向和α方向氣流力為Fe(t=0, Ω=Ωi),Fα(t=0,Ω=Ωi),Fe(t=T/4,Ω=Ωi),Fα(t=T/4,Ω=Ωi)。

將t=0，t=T/4代入式(6)可得

ΔFe(t=0,Ω=Ωi)=-a·Kee-b·Ceα·Ωi

(7)

ΔFα(t=0,Ω=Ωi)=-a·Kαe-b·Cαα·Ωi

(8)

ΔFe(t=T/4,Ω=Ωi)=-b·Keα+a·Cee·Ωi

(9)

ΔFα(t=T/4,Ω=Ωi)=-b·Kαα+a·Cαe·Ωi

(10)

以渦動頻率為橫坐標(biāo)，ΔFe(t=0,Ω=Ωi),ΔFα(t=0,Ω=Ωi),ΔFe(t=T/4,Ω=Ωi)和ΔFα(t=T/4,Ω=Ωi)為縱坐標(biāo)，如圖3所示繪制氣流力之差隨渦動頻率變化圖，分別記作曲線Ω-ΔFe(t=0,Ω=Ωi),Ω-ΔFα(t=0,Ω=Ωi),Ω-ΔFe(t=T/4,Ω=Ωi),Ω-ΔFα(t=T/4,Ω=Ωi)。

對于任意渦動頻率Ωi，取微元軌跡(Ωi,Ωi+ΔΩ)，

圖3 氣流力之差隨渦動頻率變化Fig.3 Variation of seal force difference with whirling frequency

當(dāng)ΔΩ→0時，微元段內(nèi)動力特性系數(shù)不隨渦動轉(zhuǎn)速變化而發(fā)生改變，即當(dāng)渦動轉(zhuǎn)速為Ωi和Ωi+ΔΩ時密封動力特性系數(shù)Kee,Keα,Kαe,Kαα,Cee,Ceα,Cαe,Cαα不發(fā)生改變。將渦動轉(zhuǎn)速為Ωi時對應(yīng)的作用在密封段轉(zhuǎn)子上的力ΔFe(t=0,Ω=Ωi)和渦動轉(zhuǎn)速為Ωi+ΔΩ時對應(yīng)的作用在密封段轉(zhuǎn)子上的力ΔFe(t=0,Ω=Ωi+ΔΩ)代入式(7)得

ΔFe(t=0,Ω=Ωi)=-a·Kee-b·Ceα·Ωi

(11)

ΔFe(t=0,Ω=Ωi+ΔΩ)=-a·Kee-b·Ceα·(Ωi+ΔΩ)

(12)

對微元段ΔΩ取極值得

(13)

由式(13)得擬合曲線在橫坐標(biāo)為Ωi處切線的斜率。因此過曲線Ω-ΔFe(t=0,Ω=Ωi)上任意點作切線，切線斜率為式(7)中的系數(shù)-bCeα，解得Ceα即切點橫坐標(biāo)對應(yīng)的渦動轉(zhuǎn)速下的交叉阻尼系數(shù)，切線在縱坐標(biāo)上的截距即為式(7)中的常數(shù)項-aKee，解得Kee即切點橫坐標(biāo)對應(yīng)的渦動轉(zhuǎn)速下的直接剛度系數(shù)。同理可得不同渦動轉(zhuǎn)速下的動力特性系數(shù)Kαe,Cαα,Keα,Cee,Kαα,Cαe。

由于(e,α)坐標(biāo)比(η,ξ)坐標(biāo)超前一個θ角，兩坐標(biāo)系關(guān)系為

(14)

(15)

將式(15)用矩陣形式表示

(16)

由圖2可知，用矩陣的形式表達ΔFe,ΔFα,ΔFη,ΔFξ之間的關(guān)系

(17)

易得

(18)

將式(18)用矩陣形式表示

(19)

將式(16)與式(19)合并得剛度系數(shù)為

(20)

同理，阻尼系數(shù)為

(21)

2 數(shù)值計算方法

2.1 幾何模型

本文數(shù)值模擬的密封幾何尺寸為Ertas試驗數(shù)據(jù)。表1給出了密封的具體尺寸，密封齒處細(xì)節(jié)尺寸如圖4所示。

表1 密封幾何尺寸

(a)整體模型

(b)局部尺寸圖4 密封幾何模型Fig.4 Seal geometry

2.2 數(shù)值模型

網(wǎng)格劃分是CFD計算中的重要環(huán)節(jié)，網(wǎng)格質(zhì)量和布置的合理性對計算精確性和收斂速度具有很大影響。為了提高計算精度，采用了結(jié)構(gòu)化網(wǎng)格，并對流動變化劇烈的齒頂處進行適當(dāng)加密。本文采用壁面函數(shù)法將壁面的物理量與湍流核心區(qū)聯(lián)系，并將y+值控制在30～150區(qū)間內(nèi)，密封流場網(wǎng)格如圖5所示，網(wǎng)格節(jié)點數(shù)共計114萬個。

圖5 密封流場網(wǎng)格劃分Fig.5 Grid distribution

2.3 計算工況

表2為本文數(shù)值計算的參數(shù)細(xì)節(jié)。工質(zhì)為理想空氣，采用k-ε湍流模型，湍流強度5%，轉(zhuǎn)子和靜子壁面設(shè)置為絕熱光滑無滑移壁面。在密封進口設(shè)置總溫、總壓，出口設(shè)置平均靜壓出口。瞬態(tài)計算的時間步長依據(jù)渦動頻率進行調(diào)整，本文將轉(zhuǎn)子渦動軌跡簡化為圓軌跡，渦動中心為密封幾何中心。在計算模擬時采用了動網(wǎng)格方法，運用CEL設(shè)定網(wǎng)格運動軌跡，網(wǎng)格的運動由轉(zhuǎn)子的渦動方程決定。當(dāng)轉(zhuǎn)子上偏心和垂直于偏心方向受力呈現(xiàn)周期性變化、不同周期內(nèi)對應(yīng)點受力值相差小于0.01 N、進出口流量相差小于0.1%時，認(rèn)為計算收斂。

表2 工況參數(shù)

3 計算結(jié)果與討論

將數(shù)值模擬得到的直接剛度系數(shù)Kavg、直接阻尼系數(shù)Cavg、交叉剛度系數(shù)Kxy和有效阻尼系數(shù)Ceff與Ertas等的試驗結(jié)果進行對比分析。其中Cavg,Kavg,Ceff的定義為

(22)

式中:Cxx和Cyy分別為x與y方向上的直接阻尼系數(shù)；Kxx和Kyy分別為x與y方向上的直接剛度系數(shù)。

3.1 密封動力特性系數(shù)

圖6給出了直接剛度系數(shù)隨渦動頻率的變化趨勢。模擬所得直接剛度系數(shù)Kavg在數(shù)值上小于試驗值，特別是在低頻段，符號呈相反趨勢，即密封表現(xiàn)為負(fù)剛度，且符號不隨渦動頻率變化而改變，該結(jié)論與文獻[32]一致。隨渦動頻率的增加，二者逐漸吻合，剛度系數(shù)絕對值越來越大。

圖6 直接剛度Kavg隨渦動頻率變化Fig.6 Variation of direct stiffness with swirl frequency

圖7給出了直接阻尼系數(shù)Cavg隨渦動頻率的變化趨勢?？梢钥闯?，本文計算結(jié)果與試驗基本吻合，約為600 N·s/m。

圖7 直接阻尼Cavg隨渦動頻率變化Fig.7 Variation of direct damping with swirl frequency

圖8給出了交叉剛度Kxy隨渦動頻率的變化趨勢，計算與試驗值隨渦動頻率增加，基本保持不變，約為50 kN/m。

圖8 交叉剛度Kxy隨渦動頻率變化Fig.8 Variation of cross stiffness with swirl frequency

有效阻尼系數(shù)是評價密封性能的關(guān)鍵參數(shù)，圖9給出了有效阻尼隨渦動頻率的變化趨勢?？梢钥闯觯嬎憬Y(jié)果略大于試驗值，當(dāng)渦動頻率大于100 Hz時，有效阻尼系數(shù)隨渦動頻率增加基本保持不變。

圖9 有效阻尼Ceff隨渦動頻率變化Fig.9 Variation of effective damping with swirl frequency

3.2 密封靜態(tài)穩(wěn)定性分析

上述理論與試驗結(jié)果都表明，該密封具有較好的動態(tài)穩(wěn)定性，即具有正的有效阻尼系數(shù)。然而，本文與Ertas等的研究結(jié)果均表明直接剛度系數(shù)呈現(xiàn)負(fù)值，如圖6所示。密封靜態(tài)力是轉(zhuǎn)子零轉(zhuǎn)速或沒有渦動情況下受到的靜態(tài)激振力，該力對轉(zhuǎn)子的影響主要表現(xiàn)為改變轉(zhuǎn)子系統(tǒng)剛度，引起轉(zhuǎn)子臨界轉(zhuǎn)速的變化，進而導(dǎo)致轉(zhuǎn)子不穩(wěn)定振動區(qū)域的改變；另一方面對試驗氣缸有負(fù)面影響，即會導(dǎo)致氣缸與轉(zhuǎn)子的碰摩，Picardo等[33-34]為了避免這種負(fù)面影響，專門在氣缸上額外增加支撐結(jié)構(gòu)(stiffener)。為此下文對該密封內(nèi)靜態(tài)力的形成原因開展進一步分析。

表3給出了靜態(tài)時轉(zhuǎn)子受力與泄漏情況。靜態(tài)工況下，其在偏心狀態(tài)下所受氣流激振力均為正值，即直接剛度為負(fù)值，與圖6識別結(jié)果相吻合。轉(zhuǎn)子在偏心自轉(zhuǎn)時會產(chǎn)生垂直于偏心方向的切向力，即產(chǎn)生正的交叉剛度系數(shù)。密封泄漏量在轉(zhuǎn)子轉(zhuǎn)動情況下會略微降低。

表3 密封靜態(tài)特性

圖10給出了靜態(tài)下密封內(nèi)小間隙與大間隙壓差沿泄漏方向變化情況。除第一個密封腔，其余部位的小間隙與大間隙的壓力差值小于0，即產(chǎn)生負(fù)剛度，加劇轉(zhuǎn)子偏離中心程度。

圖10 密封間隙壓差沿泄漏路徑分布Fig.10 Pressure difference distribution in the seal clearance along leakage path

圖11給出了密封最大和最小間隙內(nèi)氣流速度沿泄漏方向分布情況。小間隙處速度均大于大間隙，慣性力占主導(dǎo)作用，小間隙壓力小于大間隙，導(dǎo)致出現(xiàn)負(fù)剛度，即靜態(tài)不穩(wěn)定。

圖11 密封內(nèi)大小間隙處速度分布Fig.11 Velocity distribution in the seal clearance

4 結(jié) 論

(1)應(yīng)用瞬態(tài)動網(wǎng)格方法，提出基于微元軌跡的密封動力特性系數(shù)理論識別方法，可滿足任意轉(zhuǎn)子橢圓渦動軌跡與任意渦動頻率的密封動力特性系數(shù)求解。

(2)研究表明，本文理論計算與試驗結(jié)果具有很高的吻合度，尤其對衡量密封系統(tǒng)穩(wěn)定性的有效阻尼系數(shù)，具有較高的識別精確度。渦動頻率對試驗密封直接阻尼系數(shù)、有效阻尼系數(shù)和交叉剛度系數(shù)影響較小。

(3)試驗密封具有靜態(tài)不穩(wěn)定性。計算表明密封最大間隙氣流速度小于小間隙，慣性力占主導(dǎo)作用，大間隙處平均壓力大于小間隙處，壓力差加劇了轉(zhuǎn)子偏心，造成密封系統(tǒng)靜態(tài)不穩(wěn)定。