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

(1.中國航發(fā)商用航空發(fā)動機有限責任公司,200241,上海;2.上海市航空發(fā)動機數(shù)字孿生重點實驗室,200241,上海)

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

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

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

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

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

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

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

1 數(shù)值方法

1.1 熱聲網(wǎng)絡

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

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

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

(1)

(2)

(3)

(4)

(5)

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

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

(6)

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

(7)

(8)

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

(9)

(10)

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

(11)

(12)

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

(13)

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

(14)

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

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

Ax=a

(15)

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

1.2 激勵盤模型

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

圖2 兩級渦輪激勵盤模型示意圖

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

(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)生伴隨的脈動速度與脈動壓力,其引起的密度變化為

(27)

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

(28)

(29)

聲波引起的密度變化為

(30)

聲脈動速度在平均流方向上的投影和相對于平均流的角度分別為

(31)

i(kxx+kyy)]

(32)

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

ky=kssinvs=kξsinvξ

(33)

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

(34)

(35)

(36)

式中:下標1與2分別表示葉排的上游與下游。當葉排出口為亞聲速時,滿足庫塔條件有

(37)

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

(38)

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

A1q1=A2q2

(39)

(40)

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

(41)

(42)

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

(43)

(44)

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

對于渦輪級間的聲傳播過程,上下游葉排間的脈動幅值滿足

zu=Tdzd

(45)

(46)

Td=

(47)

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

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

Bx=b

(48)

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

1.3 耦合策略

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

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

圖3 跨部件耦合計算流程示意圖

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

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

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

2 結(jié)果與討論

2.1 模型驗證

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

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

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

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

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

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

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

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

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

2.2 應用

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

對于渦輪第一級靜葉進口和第二級動葉出口,給定邊界條件為

(49)

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

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

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

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

圖7 渦輪進口等效聲阻抗隨頻率的變化

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

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

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

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

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

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

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

圖9 特征頻率所對應的振蕩模態(tài)

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

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

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

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

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

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

3 結(jié) 論

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