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

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

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

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

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

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

(1)

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

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

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

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

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

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

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

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

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

(2)

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

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

(3)

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

(4)

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

(5)

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

(6)

應(yīng)用瞬態(tài)法得到任意渦動(dòng)轉(zhuǎn)速Ωi(i=1, 2, 3, 4, …,n)轉(zhuǎn)子所受氣流力，在t=0和t=T/4時(shí)刻，轉(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)

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

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

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

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

對(duì)微元段ΔΩ取極值得

(13)

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

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

(14)

(15)

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

(16)

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

(17)

易得

(18)

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

(19)

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

(20)

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

(21)

2 數(shù)值計(jì)算方法

2.1 幾何模型

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

表1 密封幾何尺寸

(a)整體模型

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

2.2 數(shù)值模型

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

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

2.3 計(jì)算工況

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

表2 工況參數(shù)

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

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

(22)

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

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

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

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

圖7給出了直接阻尼系數(shù)Cavg隨渦動(dòng)頻率的變化趨勢(shì)?？梢钥闯觯疚挠?jì)算結(jié)果與試驗(yàn)基本吻合，約為600 N·s/m。

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

圖8給出了交叉剛度Kxy隨渦動(dòng)頻率的變化趨勢(shì)，計(jì)算與試驗(yàn)值隨渦動(dòng)頻率增加，基本保持不變，約為50 kN/m。

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

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

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

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

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

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

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

圖10給出了靜態(tài)下密封內(nèi)小間隙與大間隙壓差沿泄漏方向變化情況。除第一個(gè)密封腔，其余部位的小間隙與大間隙的壓力差值小于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)動(dòng)網(wǎng)格方法，提出基于微元軌跡的密封動(dòng)力特性系數(shù)理論識(shí)別方法，可滿足任意轉(zhuǎn)子橢圓渦動(dòng)軌跡與任意渦動(dòng)頻率的密封動(dòng)力特性系數(shù)求解。

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

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