邵衛(wèi)東,高賢智

(1.中國(guó)航發(fā)商用航空發(fā)動(dòng)機(jī)有限責(zé)任公司,200241,上海;2.上海市航空發(fā)動(dòng)機(jī)數(shù)字孿生重點(diǎn)實(shí)驗(yàn)室,200241,上海)

為降低民用航空發(fā)動(dòng)機(jī)污染物排放,當(dāng)前的燃燒室設(shè)計(jì)一般采用貧油預(yù)混分級(jí)組織噴霧燃燒技術(shù)[1],但其易產(chǎn)生燃燒不穩(wěn)定性。該現(xiàn)象的物理本質(zhì)是火焰動(dòng)力學(xué)與聲學(xué)的雙向耦合[2-7],當(dāng)火焰產(chǎn)生的擾動(dòng)與燃燒室內(nèi)聲波在火焰處的相位差小于某一定值時(shí),擾動(dòng)幅值在初期線性增長(zhǎng),如果非線性激勵(lì)大于非線性阻尼,擾動(dòng)幅值呈指數(shù)增長(zhǎng),從而誘發(fā)低頻高幅振蕩,破壞燃燒室的結(jié)構(gòu)完整性。

對(duì)燃燒室構(gòu)型與氣動(dòng)布局進(jìn)行優(yōu)化,可有效規(guī)避熱聲振蕩的發(fā)生。優(yōu)化的物理模型即為熱聲耦合控制方程,當(dāng)前確定物理解的數(shù)值方法有大渦模擬(LES)[8-11]、混合計(jì)算流體力學(xué)-計(jì)算氣動(dòng)聲學(xué)方法[12-14]、低階模型[15]3大類:

(1)LES方法能夠準(zhǔn)確考慮渦波、熵波和聲波之間的能量轉(zhuǎn)換與非線性過程,但其對(duì)計(jì)算資源的需求仍不適用于實(shí)際燃燒室優(yōu)化設(shè)計(jì);

(2)混合方法將流動(dòng)變量分解為平均值與脈動(dòng)值,并對(duì)脈動(dòng)值采用小擾動(dòng)假設(shè)將控制方程線性化,該方法能極大降低計(jì)算量,但燃燒釋熱源項(xiàng)需要采用LES方法或者統(tǒng)計(jì)模型獲得;

(3)低階模型將實(shí)際燃燒室離散為若干單元,在單元內(nèi)考慮聲傳播的集總特征,將所有單元通過邊界匹配條件形成傳遞矩陣,雖然該方法不能考慮模態(tài)間能量轉(zhuǎn)換,但其具有直觀的物理過程,計(jì)算速度快,關(guān)聯(lián)了燃燒室?guī)缀螀?shù),有利于優(yōu)化設(shè)計(jì),故本文采用低階模型方法。

在預(yù)測(cè)熱聲振蕩過程中,在燃燒室出口通常給定隨頻率變化的阻抗邊界條件,阻抗邊界對(duì)于準(zhǔn)確評(píng)估熱聲振蕩頻率至關(guān)重要。該邊界條件可采用試驗(yàn)方法或跨部件LES方法獲得,對(duì)于不同的燃燒室工況和渦輪構(gòu)型該阻抗邊界不同,在設(shè)計(jì)階段通過反復(fù)獲得阻抗邊界的代價(jià)很大。燃燒室出口阻抗反映了聲波和熵波在燃燒室與渦輪中傳播的綜合過程,該邊界條件是跨部件物理過程的中間產(chǎn)物。

基于此,本文將燃燒室出口邊界作為傳遞變量,采用跨部件耦合方法來(lái)進(jìn)行熱聲振蕩的研究,下面將從熱聲網(wǎng)絡(luò)、激勵(lì)盤模型、耦合策略、模型驗(yàn)證和應(yīng)用等方面展開。

1 數(shù)值方法

1.1 熱聲網(wǎng)絡(luò)

熱聲網(wǎng)絡(luò)是預(yù)測(cè)熱聲振蕩過程低階模型之一,考慮到一般燃燒室頭部等結(jié)構(gòu)尺寸與聲波波長(zhǎng)的比值在0.025左右,可利用低頻近似將燃燒室分為包含直圓/環(huán)管單元、聲學(xué)緊湊單元、火焰模型和邊界條件等聲學(xué)子單元的集合,如圖1所示。由于環(huán)管單元的周向尺寸遠(yuǎn)大于其軸向和徑向尺寸,易激發(fā)起周向模態(tài),頭部等結(jié)構(gòu)的周向尺寸較小,周向模態(tài)的截止頻率相對(duì)較高,不易激發(fā)起周向模態(tài)而以軸向模態(tài)為主。通過傳遞矩陣將單元內(nèi)與單元間的聲場(chǎng)特征聯(lián)結(jié),形成熱聲系統(tǒng)方程并轉(zhuǎn)化為特征值問題,即可得到振蕩頻率及幅值分布。

圖1 環(huán)形燃燒室熱聲網(wǎng)絡(luò)示意圖

基于理想氣體、軸向均勻平均流、無(wú)黏、無(wú)擴(kuò)散、小擾動(dòng)和薄環(huán)形假設(shè),將Navier-Stokes方程在柱坐標(biāo)系下進(jìn)行線性化處理得

(1)

(2)

(3)

(4)

(5)

式中:(x,θ,R)為柱坐標(biāo)系;(ux,uθ,uR)為軸向、周向和徑向速度;ρ、p和γ分別為密度、壓力和比熱容比;變量上方橫線表示平均值,波浪線表示脈動(dòng)值。

對(duì)于不含釋熱源的聲學(xué)子單元,如擴(kuò)壓器、集氣室等,聯(lián)合方程(1)至方程(5)可化簡(jiǎn)得

(6)

式中:c和Mx分別音速和軸向平均流馬赫數(shù)。采用分離變量法,方程(6)的解表示為

(7)

(8)

式中:A為待定常數(shù);上標(biāo)±表示波傳播的方向。將式(7)代入方程(2)和方程(3),即可解析軸向與周向脈動(dòng)速度,利用脈動(dòng)壓力與軸向脈動(dòng)速度解,軸向相鄰聲學(xué)子單元的聲學(xué)特征(單模態(tài))可表示為

(9)

(10)

對(duì)于周向相鄰聲學(xué)子單元,可采用相似的方法進(jìn)行聲學(xué)特征關(guān)聯(lián),則

(11)

(12)

對(duì)于包含釋熱源的聲學(xué)子單元如火焰筒,在薄火焰面處對(duì)脈動(dòng)壓力和軸向脈動(dòng)速度利用質(zhì)量流量和壓力連續(xù)性條件可得

(13)

式中:xF+和xF-分別表示火焰面位置的兩側(cè);S為聲學(xué)單元截面的面積;Ω為釋熱源項(xiàng),一般采用火焰?zhèn)鬟f函數(shù)[16-18]對(duì)其進(jìn)行?；?則

(14)

其中n表示釋熱與聲波干涉因子,τ表示延遲時(shí)間,上標(biāo)ref表示參考位置。

將所有聲學(xué)子單元通過傳遞矩陣連接,即可形成燃燒室熱聲系統(tǒng)的控制方程

Ax=a

(15)

式中:A為系統(tǒng)特征矩陣;x為系統(tǒng)的脈動(dòng)變量;a為系統(tǒng)的源項(xiàng)。求解方程(15)對(duì)應(yīng)的特征值問題,可得到振蕩頻率和振蕩模態(tài)分布。

1.2 激勵(lì)盤模型

由于渦輪葉片軸向弦長(zhǎng)與聲波波長(zhǎng)的比值一般低于0.1,葉柵稠度較高且動(dòng)葉通過頻率及其諧波對(duì)應(yīng)的波長(zhǎng)遠(yuǎn)小于發(fā)生熱聲振蕩時(shí)的聲波波長(zhǎng),激勵(lì)盤模型可利用長(zhǎng)波長(zhǎng)假設(shè)將葉排在軸向簡(jiǎn)化為無(wú)限薄盤片,如圖2所示。入射的熵波、渦波和聲波進(jìn)入渦輪級(jí),熵波與渦波全部透射而一部分聲波被葉片散射,另一部分聲波直接透過葉片通道;熵波與渦波伴隨著當(dāng)?shù)貧饬魉俣认蛳掠蝹鞑?聲波同時(shí)向上游與下游傳播;通過守恒定律確立盤片上下游脈動(dòng)量之間的關(guān)系。雖然渦輪葉片的軸向弦長(zhǎng)較短,但葉排的軸向間距相對(duì)于波長(zhǎng)不一定可以忽略,故必須考慮熵波、渦波和聲波在葉排軸向間隙的傳播過程。

圖2 兩級(jí)渦輪激勵(lì)盤模型示意圖

基于理想氣體、均勻平均流、無(wú)黏、無(wú)熱導(dǎo)、長(zhǎng)波長(zhǎng)、高稠度和忽略葉排內(nèi)聲傳播徑向變化的假設(shè),將Navier-Stokes方程在二維坐標(biāo)系下進(jìn)行線性化處理得

(16)

(17)

(18)

(19)

(20)

式中:(ux,uy)為x和y方向的速度;s為熵。

將熵波、渦波和聲波用傅里葉形式表示為

(21)

(22)

(23)

式中:σ、ω、ks、vs分別為熵波幅值、圓頻率、熵波數(shù)、波數(shù)與x軸的夾角,ω=2πf,f為頻率;ξ、vξ分別為渦波幅值、渦波數(shù)與x軸的夾角;kx、ky分別為聲波數(shù)k在x與y方向的分量;變量上方尖角表示傅里葉分量。將方程(21)～方程(23)代入方程(16)～方程(19)并整理,可得色散關(guān)系為

(24)

(25)

(26)

式中:ca、M、θ分別為聲速、馬赫數(shù)、平均流速度方向與x軸的夾角。

由于熵波不產(chǎn)生伴隨的脈動(dòng)速度與脈動(dòng)壓力,其引起的密度變化為

(27)

渦波不產(chǎn)生伴隨的脈動(dòng)密度和脈動(dòng)壓力,渦脈動(dòng)速度在平均流方向上的投影和相對(duì)于平均流的角度分別為

(28)

(29)

聲波引起的密度變化為

(30)

聲脈動(dòng)速度在平均流方向上的投影和相對(duì)于平均流的角度分別為

(31)

i(kxx+kyy)]

(32)

由于葉排內(nèi)脈動(dòng)具有周期性,熵波、渦波與聲波在y方向具有一致波數(shù)

ky=kssinvs=kξsinvξ

(33)

對(duì)熵守恒、質(zhì)量守恒和總焓守恒關(guān)系進(jìn)行線性化處理可得

(34)

(35)

(36)

式中:下標(biāo)1與2分別表示葉排的上游與下游。當(dāng)葉排出口為亞聲速時(shí),滿足庫(kù)塔條件有

(37)

當(dāng)葉排出口為超聲速時(shí),葉片通道內(nèi)部存在喉口截面,質(zhì)量守恒關(guān)系的線性化形式為

(38)

在單葉排內(nèi),可將守恒關(guān)系表示成

A1q1=A2q2

(39)

(40)

對(duì)于葉片出口為亞聲速平均流情形,矩陣A1和A2分別為

(41)

(42)

式中:ζ=[1+M2(γ-1)/2]-1。對(duì)于葉片出口為超聲速平均流情形,矩陣A1和A2分別為

(43)

(44)

脈動(dòng)幅值與復(fù)合變量q的關(guān)系可直接通過式(21)～式(23)進(jìn)行關(guān)聯(lián)。

對(duì)于渦輪級(jí)間的聲傳播過程,上下游葉排間的脈動(dòng)幅值滿足

zu=Tdzd

(45)

(46)

Td=

(47)

式中:|δx|為上下游葉片間距;下標(biāo)中+、-分別表示聲波傳播的方向。

將所有葉排通過傳遞矩陣連接,即可形成渦輪內(nèi)熵波、渦波和聲波傳播的控制方程

Bx=b

(48)

式中:B為等價(jià)的傳播矩陣;x為系統(tǒng)的脈動(dòng)變量;b為等價(jià)的入射波源項(xiàng)。

1.3 耦合策略

在燃燒室與渦輪系統(tǒng)內(nèi)分別求解方程(15)和方程(48)可得到各自的特征解,在非耦合情形下,方程(15)中矩陣A不能考慮渦輪向燃燒室反射聲波的影響,同時(shí)方程(48)中矩陣B和向量b不能考慮燃燒室向渦輪入射熵波、渦波和聲波的影響。

為求解跨部件系統(tǒng)特征頻率,本文引入介于燃燒室與渦輪界面處聲阻抗作為過渡變量,通過其在燃燒室與渦輪間進(jìn)行信息傳遞,如圖3所示。

圖3 跨部件耦合計(jì)算流程示意圖

渦輪進(jìn)口聲阻抗是振蕩頻率、周向波數(shù)、渦輪幾何與氣動(dòng)參數(shù)的非線性函數(shù),對(duì)于給定渦輪集總參數(shù)構(gòu)型,聲阻抗可看作只是頻率的非線性函數(shù),故渦輪對(duì)燃燒不穩(wěn)定性的影響通過過渡變量間接進(jìn)入熱聲系統(tǒng)控制方程(15)中,火焰?zhèn)鬟f函數(shù)的非線性作用直接體現(xiàn)在矩陣A中。

具體計(jì)算步驟為:先給定渦輪進(jìn)口邊界條件和假定頻率,將這些初邊值條件輸入到激勵(lì)盤模型,結(jié)合渦輪葉排結(jié)構(gòu)與氣動(dòng)參數(shù)計(jì)算并裝配矩陣B和向量b,得到渦輪進(jìn)口的前傳聲波幅值與相位,從而進(jìn)一步計(jì)算得到聲阻抗;將聲阻抗輸入到熱聲網(wǎng)絡(luò)模型,結(jié)合燃燒室結(jié)構(gòu)與氣動(dòng)參數(shù)計(jì)算并裝配矩陣A和向量a,通過內(nèi)迭代求解非線性特征值問題得到振蕩頻率及燃燒室出口相應(yīng)的熵波、渦波和聲波幅值與相位;將得到的振蕩頻率作為耦合系統(tǒng)外迭代的初始值重復(fù)進(jìn)行直至收斂;判斷收斂的依據(jù)是外迭代過程相鄰兩次頻率之差不大于10-5,根據(jù)測(cè)試算例,外迭代過程一般不大于10次即可收斂;耗時(shí)模塊仍然是求解特征值問題的內(nèi)迭代過程。

雖然耦合計(jì)算模型的總體計(jì)算量約為單獨(dú)熱聲網(wǎng)絡(luò)模型計(jì)算量的4～7倍,但由于單獨(dú)熱聲網(wǎng)絡(luò)模型計(jì)算過程的耗時(shí)較短,耦合計(jì)算過程的耗時(shí)仍是工程實(shí)際計(jì)算過程可以接受的。

2 結(jié)果與討論

2.1 模型驗(yàn)證

由于跨部件耦合計(jì)算相對(duì)單獨(dú)燃燒室熱聲振蕩計(jì)算,增加了渦輪激勵(lì)盤模型,故耦合策略的驗(yàn)證分兩步:首先驗(yàn)證激勵(lì)盤模型的準(zhǔn)確性,再驗(yàn)證跨部件耦合計(jì)算的準(zhǔn)確性,后者主要在下面的應(yīng)用中進(jìn)行陳述。

對(duì)于激勵(lì)盤模型的驗(yàn)證,選取Cumpsty的單葉排算例[19],其氣動(dòng)參數(shù)如表1所示。

表1 單葉排氣動(dòng)參數(shù)

當(dāng)葉排上游僅有熵波入射時(shí),上游反射聲波與下游透射聲波隨歸一化頻率的變化關(guān)系如圖4所示,可見本文的計(jì)算結(jié)果與Cumpsty的結(jié)果吻合較好。由于葉排上游平均流為軸向,y方向上的反射聲波模態(tài)峰值一致,均出現(xiàn)在截?cái)囝l率處。由圖4b可見,葉片通道內(nèi)氣流加速與折轉(zhuǎn)引起的周向旋流顯著改變了y方向上的透射聲波。

(a)葉排上游反射聲波與歸一化頻率的關(guān)系

(b)葉排下游透射聲波與歸一化頻率的關(guān)系圖4 葉排上游僅熵波入射工況下聲波與歸一化頻率的關(guān)系

當(dāng)葉排上游僅有聲波入射時(shí),上游反射聲波與下游透射聲波隨歸一化頻率的變化關(guān)系如圖5所示,本文的計(jì)算結(jié)果與Cumpsty的結(jié)果再次吻合較好,驗(yàn)證了激勵(lì)盤模型的準(zhǔn)確性。對(duì)于多葉排或多級(jí)渦輪的情形,測(cè)試結(jié)果表明當(dāng)葉排的個(gè)數(shù)大于等于3或者渦輪級(jí)數(shù)大于等于2時(shí),激勵(lì)盤模型計(jì)算得到的渦輪前等效聲阻抗已不隨葉排的個(gè)數(shù)而變化。對(duì)于多葉排或多級(jí)渦輪各分量的獨(dú)立驗(yàn)證,將在后續(xù)的分層仿真規(guī)劃中考慮。

(a)葉排上游反射聲波與歸一化頻率的關(guān)系

(b)葉排下游透射聲波與歸一化頻率的關(guān)系圖5 葉排上游僅聲波入射工況下聲波與歸一化頻率的關(guān)系

2.2 應(yīng)用

選取某航空發(fā)動(dòng)機(jī)高壓渦輪的前兩級(jí)作為應(yīng)用算例,其氣動(dòng)參數(shù)(氣流角在靜葉中為絕對(duì)值而在動(dòng)葉中為相對(duì)值)如表2所示。

對(duì)于渦輪第一級(jí)靜葉進(jìn)口和第二級(jí)動(dòng)葉出口,給定邊界條件為

(49)

圖6給出了振蕩頻率為400 Hz下等效聲阻抗Z(以渦輪進(jìn)口密度和聲速進(jìn)行歸一化)隨模態(tài)數(shù)(以kyRm表示,Rm為渦輪第一級(jí)中徑)的變化關(guān)系,當(dāng)kyRm=0,±1時(shí),聲阻抗變化劇烈,當(dāng)模態(tài)數(shù)較大時(shí)聲阻變化平緩而聲抗保持不變。這是因?yàn)榻o定頻率下,渦輪內(nèi)只能激發(fā)出3種模態(tài),而大部分高階模態(tài)被截止。

表2 多葉排氣動(dòng)參數(shù)

圖6 渦輪進(jìn)口等效聲阻抗隨模態(tài)數(shù)的變化

圖7給出了給定模態(tài)數(shù)下渦輪進(jìn)口等效聲阻抗隨頻率的變化。由圖可以看出:不同模態(tài)數(shù)下聲阻抗具有相似的變化規(guī)律,聲阻隨著頻率增大而增大,聲抗隨著頻率增大呈現(xiàn)先遞減再緩緩增加然后減小的趨勢(shì),在250～600 Hz范圍內(nèi)的變化梯度較小,而該頻率范圍包含熱聲振蕩發(fā)生的頻率區(qū)間。

圖7 渦輪進(jìn)口等效聲阻抗隨頻率的變化

對(duì)于真實(shí)的燃燒室,化學(xué)反應(yīng)放熱使得燃燒室內(nèi)溫度梯度較大,對(duì)比聲學(xué)網(wǎng)絡(luò)模型計(jì)算得到的燃燒室出口聲波幅值和熵波幅值,發(fā)現(xiàn)燃燒室出口熵波幅值一般為聲波幅值的1～30倍。

(a)單獨(dú)熱聲網(wǎng)絡(luò)模型結(jié)合無(wú)反射邊界條件 計(jì)算所得增益因子與頻率分布

(b)單獨(dú)熱聲網(wǎng)絡(luò)模型結(jié)合全反射邊界條件 計(jì)算所得增益因子與頻率分布

(c)單獨(dú)熱聲網(wǎng)絡(luò)模型結(jié)合一維平面波方程導(dǎo)出 阻抗邊界條件計(jì)算所得增益因子與頻率分布

(d)耦合熱聲網(wǎng)絡(luò)模型與激勵(lì)盤模型計(jì)算 所得頻譜分布圖8 燃燒室發(fā)生熱聲振蕩下的頻譜分布

在較小熵波幅值下,傳統(tǒng)的熱聲網(wǎng)絡(luò)模型預(yù)測(cè)熱聲振蕩需要在燃燒室出口給定邊界條件,通常有無(wú)反射邊界條件、全反射邊界條件和基于一維平面波方程得到的阻抗邊界條件3種方式。不同轉(zhuǎn)速比(核心機(jī)轉(zhuǎn)速與額定轉(zhuǎn)速之比)下單獨(dú)的熱聲網(wǎng)絡(luò)模型結(jié)合3類邊界條件的預(yù)測(cè)結(jié)果如圖8a～8c所示,其中無(wú)反射邊界條件無(wú)法判定熱聲振蕩過程的發(fā)生,全反射邊界條件判定所有工況均發(fā)生熱聲振蕩,第3類邊界條件發(fā)生熱聲振蕩的頻率對(duì)應(yīng)的亥姆霍茲數(shù)He(He=ωRm/ca)為1.02,而試驗(yàn)測(cè)得熱聲振蕩的亥姆霍茲數(shù)為1.36,誤差為25%。

耦合計(jì)算獲得的亥姆霍茲數(shù)如圖8d所示,其中發(fā)生熱聲振蕩的頻率對(duì)應(yīng)的He為1.36,該值與試驗(yàn)結(jié)果相同,二者的一致性證明了耦合計(jì)算模型的可靠性。圖9給出了特征頻率所對(duì)應(yīng)的振蕩模態(tài),可見該模態(tài)近似為周向負(fù)一階模態(tài),但該模態(tài)波節(jié)和波節(jié)之間的相位差略小于180°。

圖9 特征頻率所對(duì)應(yīng)的振蕩模態(tài)

圖10和圖11分別給出了振蕩頻率下歸一化聲波和熵波幅值A(chǔ)+、Ae隨周向模態(tài)數(shù)的變化。當(dāng)周向模態(tài)數(shù)為-1時(shí),聲波和熵波的幅值達(dá)到最大,而其他模態(tài)數(shù)下幅值均可忽略,周向模態(tài)數(shù)為-1時(shí)的熵波幅值大約為聲波幅值的20倍。

圖10 燃燒室出口歸一化聲波幅值隨周向模態(tài)數(shù)的變化

圖11 燃燒室出口歸一化熵波幅值隨周向模態(tài)數(shù)的變化

圖12給出了振蕩頻率下的等效聲阻抗隨周向模態(tài)數(shù)的變化,與圖6中各模態(tài)下聲波和熵波均勻入射情形不同,真實(shí)燃燒室出口聲波和熵波入射條件下等效聲阻抗隨周向模態(tài)數(shù)呈非對(duì)稱特征,當(dāng)周向模態(tài)數(shù)為2時(shí)聲阻最大而聲抗最小,當(dāng)周向模態(tài)數(shù)為-1時(shí)聲阻次大而聲抗次小。這說(shuō)明燃燒室熱聲振蕩主要受到周向正二階模態(tài)和周向負(fù)一階模態(tài)的影響,也是圖9中模態(tài)不完全對(duì)稱的原因。

圖12 燃燒室出口等效聲阻抗隨周向模態(tài)數(shù)的變化

值得一提的是,當(dāng)激勵(lì)盤入射熵波幅值遠(yuǎn)大于聲波幅值時(shí),耦合求解無(wú)法獲得振蕩頻率?？赡艿脑蚴侨肷潇夭ㄔ谧矒舻綔u輪S1葉片前緣時(shí)發(fā)生散射并在后續(xù)葉排間傳播時(shí)發(fā)生耗散,從而熵波與聲波間發(fā)生能量轉(zhuǎn)換。本研究未考慮入射熵波的散射和耗散效應(yīng),在下一步的工作中將考慮該影響。

3 結(jié) 論

本文針對(duì)航空發(fā)動(dòng)機(jī)燃燒室熱聲振蕩預(yù)測(cè)中實(shí)際出口阻抗邊界不定問題,提出了跨部件耦合計(jì)算策略,引入燃燒室與渦輪界面處聲阻抗作為傳遞變量,在求解過程中增加了系統(tǒng)間的外迭代。通過經(jīng)典算例和工程算例的驗(yàn)證,相對(duì)于單獨(dú)燃燒室熱聲計(jì)算,跨部件耦合計(jì)算方法能夠給出準(zhǔn)確的振蕩頻率和模態(tài)分布,該方法可進(jìn)一步推廣至高精度熱聲預(yù)測(cè)模型。