核電超長末葉氣動穩(wěn)定性的流固耦合數(shù)值研究

庾明達, 徐自力, 范志飛, 曹守洪, 范小平, 方宇

(1.西安交通大學航天航空學院, 710049, 西安; 2.西安交通大學機械結(jié)構(gòu)強度與振動國家重點實驗室, 710049, 西安; 3.東方電氣集團東方汽輪機有限公司, 618000, 四川德陽)

汽輪機中的顫振是一種持續(xù)、不衰減的自激振動現(xiàn)象,是由汽輪機流道內(nèi)的氣流力與葉片發(fā)生相互作用引起,所以研究葉片的氣動穩(wěn)定性需考慮流固耦合的作用。目前,研究顫振穩(wěn)定性應(yīng)用最多的流固耦合方法為能量法,仲繼澤等通過能量法預(yù)測了相位角、流量等參數(shù)對壓氣機葉片的顫振穩(wěn)定性的影響[1-2]。能量法思路清晰、簡便,但由于假設(shè)過多,與實際情況不很相符[3]。另有很多學者采用流固耦合方法來模擬葉片的實際振動響應(yīng)特征,Vahdati等用強耦合CFD方法研究了寬弦風扇的失速顫振現(xiàn)象,但此方法求解難度大,占用資源高[4];Dowell等基于帶副翼操縱面的機翼3DOF模型,運用ROM降階技術(shù)研究了臨界音速條件下馬赫數(shù)對機翼顫振速度和穩(wěn)定性的影響[5],這種氣動力降階模型目前多局限于機翼等具有簡單幾何外形結(jié)構(gòu)的流固耦合計算;文獻[6-7]利用模態(tài)疊加法將葉片響應(yīng)形式分解為單階模態(tài)的線性振動,研究了跨音速壓氣機的顫振邊界。還有不少學者使用二維模型模擬壓氣機的流場特征[8-10],但二維模型難以完全模擬高負載葉輪機械的流動狀況。

目前,許多研究文獻雖利用上述簡化耦合模型研究了流場流動和結(jié)構(gòu)振動特征,然而大多數(shù)文獻集中于研究機翼或壓氣機葉片的顫振穩(wěn)定性[11-14],對核電汽輪機末級葉片,尤其是本文所研究的超長核電末級葉片的氣動特征鮮有涉及。

因此,本文采用三維流固耦合時域響應(yīng)算法,研究了新開發(fā)的核電汽輪機超長末級葉片在設(shè)計背壓3 600 Pa、高背壓7 200 Pa和10 800 Pa工況下的氣動穩(wěn)定性問題。計算結(jié)果表明,在所計算的背壓工況內(nèi)葉片維持氣動穩(wěn)定,未發(fā)現(xiàn)顫振現(xiàn)象,但穩(wěn)定性逐漸降低。

1 計算原理和方法

1.1 流體控制方程

為完全模擬核電葉片不同葉高的流動特征,采用三維Navier-Stokes(N-S)方程作為控制方程,該方程組表示流體運動所應(yīng)滿足的質(zhì)量守恒,動量守恒和能量守恒三大物理學定律,可寫為

(1)

式中:Q為瞬態(tài)矢量;F、G、H分別為3個方向上的通量;S為源性通量。

雖然瞬時N-S方程可以用來描述湍流,但其非線性特征不易求解,工程上常引入雷諾應(yīng)力的方法將N-S方程中瞬態(tài)的變化量以時均化的形式表現(xiàn),引入湍流控制方程把湍流的脈動值與時均值聯(lián)系起來,方程變?yōu)槔字Z平均納維斯托克斯方程(RANS)。本文采用的k-ε湍流控制方程為

(2)

(3)

1.2 固體求解算法

固體域的動力學方程可寫為

(4)

式中:M、C、K分別為質(zhì)量、阻尼和剛度矩陣;q(t)為位移;p(t)為作用在流固耦合交界面上的氣動壓力。采用直接積分法求解葉片的振動方程得到葉片的有限元節(jié)點振動位移。

1.3 流固耦合的方法及計算流程

流固耦合計算的實質(zhì)是流體和固體的交界面上發(fā)生相互作用,進而產(chǎn)生耦合并相互影響。由于流場網(wǎng)格與固體網(wǎng)格并不一一對應(yīng),所以在流固耦合交界面上,需運用樣條方法,將流固交界面節(jié)點位移做轉(zhuǎn)換

Δxa=GΔxs

(5)

式中:xa為流場位移;xs為固體位移;G為樣條矩陣。反之,流場氣動力需要施加到固體交界面上,根據(jù)做功等效原理有

qs=GTqa

(6)

式中:qa、qs為施加在流固交界面上的氣動力。在流固耦合計算過程中,流場完成內(nèi)部迭代計算后,按照式(6),與固體交界面交換壓力等變量,之后固體產(chǎn)生位移,此時流場網(wǎng)格變形,以便進行下一次的迭代計算。在流固耦合交界面處,流場網(wǎng)格節(jié)點直接進行變形,其他較遠處的網(wǎng)格變形遵循以下位移擴散方程

·(Γdispδ)=0

(7)

(8)

圖1 時域響應(yīng)流固耦合算法的流程圖

2 末級汽輪機有限元網(wǎng)格模型

本文研究的對象為新開發(fā)的核電汽輪機,末級葉片長度超過2 m,靜葉數(shù)為66,動葉數(shù)為76,蒸汽從末級葉片的靜葉入口流入,隨后經(jīng)動葉出口流出,幾何示意圖如圖2所示。

圖2 末級葉片通道幾何示意圖

本文采用三維計算模型和循環(huán)對稱邊界的單通道流道作為計算對象,用以計算末級流場,這樣可有效減少計算規(guī)模,提高計算效率。為了選取合適的網(wǎng)格,本文采用由疏到密具有不同單元數(shù)量、不同網(wǎng)格尺寸的4套流場網(wǎng)格,對設(shè)計背壓工況的定常流動進行了計算,得到了流體施加在動葉片中部截面背弧面的壓力分布,如圖3所示。圖3中,X為葉片中部展弦方向任意一點的位置坐標,C為葉片中部的弦長。

圖3 設(shè)計背壓工況下的動葉中部背弧面流體壓力分布

由圖3可知,隨著網(wǎng)格單元數(shù)的不斷增多,計算得到的背弧面靜壓力分布趨于定值,在流場網(wǎng)格單元數(shù)為707 100時,繼續(xù)增加網(wǎng)格單元數(shù),背弧面的靜壓力計算結(jié)果也基本保持不變。因此,對該核電末級葉片,單元數(shù)選取707 100。

圖4給出了核電末級葉片的網(wǎng)格,其中葉片固體網(wǎng)格含有13 041個單元。動靜葉流體計算域的網(wǎng)格均為結(jié)構(gòu)化網(wǎng)格,對葉片周圍的流動區(qū)域設(shè)置了邊界層,內(nèi)部計算域網(wǎng)格環(huán)繞邊界層分布直至充滿整個計算域。動、靜葉網(wǎng)格的平均第一層邊界層網(wǎng)格高度約為3.1×10-5m,對第一層網(wǎng)格高度做無量綱化處理,可得流場網(wǎng)格的y+≈4<5。

(a)末級流道整體 (b)葉片固體

(c)靜葉葉頂流道 (d)動葉葉頂流道圖4 核電末級葉片的計算網(wǎng)格

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

3.1 流場特征分析

采用單通道模型,在流場計算域中,動靜葉流道兩側(cè)均施加循環(huán)對稱邊界,將靜葉與動葉流場域統(tǒng)一視為內(nèi)部區(qū)域,以此模擬動葉與靜葉之間的相對旋轉(zhuǎn)運動,葉片表面、葉頂、葉轂為無滑移壁面,靜葉進口設(shè)置蒸汽的流量和溫度,在動葉出口設(shè)置壓力值。對設(shè)計背壓3 600 Pa和7 200、10 800 Pa兩個高背壓工況下的定常流動進行了計算,得到了3 600、7 200 Pa工況下末級流道截面流線分布,如圖5所示,可知在設(shè)計背壓時,流體從動葉流出后,流動保持平穩(wěn),而在提高背壓后,流動未發(fā)生分離或是方向改變。而當背壓提高至10 800 Pa時,動葉尾部流道內(nèi)發(fā)生了旋渦現(xiàn)象,此時流線分布及壓力分布如圖6所示。由圖6可知,流體是順著壓力梯度向下游流動,其流動趨勢與壓力梯度保持一致,然而在該背壓工況下動葉流道中部的壓力梯度發(fā)生了明顯的波動,出現(xiàn)幾個局部低壓凸口,流體順著壓力梯度梯度流經(jīng)這些區(qū)域時,由于動能已經(jīng)減小,不足以克服壓差導(dǎo)致流向改變,少量流體在動葉流道近出口處形成了旋渦。

(a)背壓為3 600 Pa (b)背壓為7 200 Pa圖5 子午面流線分布

(a)流線分布 (b)壓力分布圖6 背壓為10 800 Pa時子午面流線分布及壓力分布

3.2 葉片氣動穩(wěn)定性分析

(a)3 600 Pa響應(yīng)曲線 (b)3 600 Pa頻譜圖

(c)7 200 Pa響應(yīng)曲線 (d)7 200 Pa頻譜圖

(e)10 800 Pa響應(yīng)曲線 (f)10 800 Pa頻譜圖圖7 各背壓工況下的尾緣響應(yīng)曲線及其頻譜圖

在不同背壓工況下動葉頂尾緣點上的時間歷程響應(yīng)曲線及其頻譜圖如圖7所示。由圖7可知,3個背壓工況下的響應(yīng)曲線全部現(xiàn)為衰減趨勢,響應(yīng)曲線中的主頻成分與動葉的第1階固有動頻ω1的比值接近1,即葉片是以動葉的第1階固有動頻為主進行自由衰減響應(yīng),對葉片進行的模態(tài)分析結(jié)果顯示,葉片第2,3,4階模態(tài)與第1階模態(tài)的固有頻率之比分別為1.61,2.57,3.39。由圖7d可知,頻譜圖中含有葉片2階固有頻率附近的振動成分,表明在背壓為7 200 Pa時,葉片表現(xiàn)為1階、2階模態(tài)的疊加振動。隨著背壓的逐步增大,葉片的幅值均分別呈現(xiàn)出逐漸減小的趨勢,即葉片在這3個背壓工況下都將保持氣動穩(wěn)定,不會發(fā)生顫振。

表1、表2給出了在3種工況下,葉片振動至最大位移和最小位移時流場的軸向力和周向力。由表1、表2可知,在3個背壓條件下,當葉片從最大位移處振動至最小位移處時,流場的軸向力和周向力均發(fā)生了變化,即固體的振動影響了流體氣動力的分布。

表1 不同背壓下葉片振動至最大、最小位移時的流場軸向力

表2 不同背壓下葉片振動至最大、最小位移時的流場周向力

對數(shù)阻尼比的計算公式為

(9)

式中:m為周期間隔;A為振動幅值。計算時,由于該核電葉片的材料阻尼很小,式(7)的阻尼比主要為氣動阻尼比。氣動阻尼比與背壓的關(guān)系如圖8所示,可知隨著背壓的提高,3個工況下氣動阻尼均為正,但提高背壓后阻尼比隨之減小,葉片的氣動穩(wěn)定性隨之減弱。

圖8 氣動阻尼比與背壓的關(guān)系

背壓工況下振動歷程的相平面圖如圖9所示。由圖9可知:在背壓為3 600、7 200 Pa時葉片振動的速度和位移均呈現(xiàn)出雜亂的變化趨勢,相平面圖呈現(xiàn)出多個頻率和幅值的疊加運動,氣動阻尼對結(jié)構(gòu)有顯著影響,相平面呈現(xiàn)不規(guī)則的收斂趨勢,體現(xiàn)了強烈的非線性特征;背壓升至10 800 Pa時,頻率成分變得單一,此時的阻尼比相當小,有限時間段內(nèi)葉片在相平面圖上收斂趨勢變得更加緩慢,隨著時間的進一步推進,葉片最終收斂到靜平衡位置。

(a)背壓為3 600 Pa (b)背壓為7 200 Pa

(c)背壓為10 800 Pa圖9 各背壓下振動響應(yīng)的相平面圖

80%葉高截面流場馬赫數(shù)云圖和截面展弦向葉片壓力分布如圖10、11所示。由圖10、11可知,葉柵流道內(nèi)出現(xiàn)了激波,激波不斷朝葉片前緣移動,在背壓為3 600、7 200 Pa時,流場的背弧面均可觀察到穩(wěn)定的激波;到背壓為10 800 Pa時,葉片背弧面上出現(xiàn)了低能流團,同時邊界層出現(xiàn)了分離,并促使激波結(jié)構(gòu)發(fā)生了變化,此時的激波分為兩支,一支為前緣正激波,一支與流道平行使整個激波呈“人”字形分布。

(a)背壓為3 600 Pa (b)背壓為7 200 Pa

(c)背壓為10 800 Pa圖10 80%葉高截面流場馬赫數(shù)云圖

圖11 80%葉高截面展弦向葉片壓力分布

文獻[15-16]指出,流道內(nèi)激波的位置變化是導(dǎo)致葉片氣動彈性穩(wěn)定性差異的主要誘因,當激波從葉尖前緣脫體后,葉片從氣流中吸收能量導(dǎo)致振動發(fā)散并迅速喪失其氣動穩(wěn)定性。本文研究中,背壓為10 800 Pa工況時動葉流道內(nèi)雖已出現(xiàn)靠近前緣的分岔激波,但并沒有脫離動葉前緣葉柵流道,背壓的提高使動葉柵內(nèi)的激波向前移動,但激波沒有發(fā)生脫體,這在一定程度上解釋了3個工況下核電葉片保持氣動穩(wěn)定的原因;在葉片振動相對劇烈的葉高中上部流道內(nèi),流動保持穩(wěn)定,沒有紊流產(chǎn)生。當背壓為10 800 Pa及葉片開始振動后的最大、最小位移時,提取葉片0.5、0.7葉高截面和靠近葉尖振動觀察點的0.95葉高截面的流動分布情況,如圖12所示。由圖12可知,葉片最大、最小位移時葉片3個截面的動葉葉柵的氣流分別沿著葉片內(nèi)、背弧面向下游流動,未發(fā)生附面層分離,葉片保持氣動穩(wěn)定,同時流場在各個截面處的流動平穩(wěn)。

(a)最大位移時 (b)最小位移時圖12 背壓為10 800 Pa時葉片不同振動時刻的瞬態(tài)流線分布

4 結(jié) 論

本文采用三維時域響應(yīng)流固耦合算法,研究了新開發(fā)核電汽輪機末級葉片從設(shè)計背壓為3 600 Pa逐漸提高到高背壓為7 200、10 800 Pa條件下的氣動穩(wěn)定性問題。計算結(jié)果表明,隨著背壓的升高,在葉柵流道內(nèi)的近出口處引起了流體壓力的擾動,在壓差的影響下逐漸形成了旋渦現(xiàn)象。在從設(shè)計背壓提高至兩個高背壓工況下,核電汽輪機葉片的時域響應(yīng)均以1階動頻為主分別進行自由衰減,表明該核電超長葉片在所計算的背壓范圍內(nèi)將保持氣動穩(wěn)定。背壓提高后,動葉柵和流道內(nèi)的激波隨著背壓的增大不斷向前緣移動,但并未發(fā)生脫體,葉片振動至不同時刻時在葉高各截面流道內(nèi)流動沒有發(fā)生分離,這解釋了葉片保持氣動穩(wěn)定的原因。但是,氣動阻尼比隨背壓的提高逐漸消減,氣動穩(wěn)定性隨背壓的提高而逐漸減弱。