特征值理論在穩定性預測中的應用研究進展綜述

孫曉峰，董旭，張光宇，王卓，孫大坤，*

1. 北京航空航天大學 能源與動力工程學院，北京 100191 2. 北京航空航天大學 航空發動機研究院，北京 100191

長期以來，中國在航空發動機領域的發展始終徘徊在舉步維艱的探索階段，在發動機設計、制造、試驗以及材料工藝方面同歐美國家存在著一定的差距，這導致了目前以殲20、運20和C919為代表的軍用民用先進飛機遲遲無法裝備“中國心”。在“兩機”專項的重大部署和國家科技發展前瞻性布局的支持引導下，與航空發動機相關的諸多工程和基礎研究領域的問題已經得到了解決，而進一步深入發掘工程問題背后的基礎理論與方法，積極儲備面向未來的先進技術，對于中國現有型號的研制和未來型號預研都具有重要意義。

現代先進航空發動機的發展方向不僅是追求更大的推力和更高的效率，大涵道比民用航空發動機也在同時追求更低的排放和噪聲，而這些追求的背后也蘊藏了諸多工程和基礎研究領域的關鍵科學問題，不僅包括如何利用氣動熱力學基本原理實現高性能設計、如何通過模擬燃燒與化學反應動力學過程提高燃燒效率和減少污染物排放、以及如何利用氣動聲學相關理論降低風扇和噴流噪聲等，更是引發了各種亟待解決的復雜流固熱聲耦合問題，尤其是高推重比與高效率氣動設計帶來了葉輪機械負荷的增加與流動復雜性的增強，這勢必會造成壓氣機穩定性裕度不足、風扇噪聲增強、葉片流致/渦致振動和燃燒室熱聲振蕩等一系列危及發動機穩定工作和飛機飛行安全的工程問題。

作為航空發動機最主要的3大部件，風扇/壓氣機、燃燒室和渦輪，其各自均對發動機整機性能和可靠性起著決定性的影響。對于壓縮部件，高壓比、高性能的氣動設計的工程應用主要受到穩定性問題的制約。通常，風扇/壓氣機的工作圖譜存在著左側失穩邊界，而發動機穩定工作在共同工作線上，為了避免由于過渡態以及進氣畸變等造成工作點偏移越過失穩邊界遭遇失速和喘振，工作線與失穩邊界需要保持一定的距離，也就是所謂的穩定性裕度。半個多世紀的研究與實踐均著眼于如何有效地將失穩邊界左移，并在航空發動機壓縮系統所需的工作范圍內避免失穩現象(失速和喘振)的發生，研究不僅涵蓋了其非定常演化的特點，同時也發展了一系列的主/被動穩定性控制方法。早期的關于葉輪機械內部流動失穩問題的各類解析模型也大都關注于失速喘振的起始階段，從而建立各種各樣的特征值問題。然而，為了得到流場擾動的解析解，這些模型中需要做大量的簡化和假設。模型的預測精度取決于風扇/壓氣機損失、落后角和壓比等特性曲線的準確性，而這些特性曲線常常又不具有足夠的精度，特別是在新的壓氣機型號設計階段。大量對經驗關系的需求限制了這些模型的應用價值，特別是在實際工程中的應用。因此，發展能夠考慮復雜流場及幾何條件，不嚴重依賴經驗關系，適用于葉輪機械設計階段，并且對流動失穩邊界可以較為快速準確地給出清晰判據的穩定性模型勢在必行，非常富有理論和工程指導意義。

不過即便在設計階段保留了足夠的氣動穩定裕度并采取了穩定性控制措施，真實工況下運轉的壓氣機和風扇在部分中低轉速時會在近失速區發生葉片顫振，英國帝國理工大學Vahdati教授團隊與Rolls-Royce開展了長期的合作，通過復雜的CFD計算探究了這種氣彈失穩現象的特性，提出了失速邊界在部分中低轉速氣彈失穩影響下存在“Flutter Bite”現象[1]，二者共同決定了風扇的穩定邊界。氣動彈性穩定性問題涉及到復雜的氣體和固體相互作用，早期學者的研究大多基于解耦假設，因此將氣動彈性穩定性問題重點放在求解葉片表面氣動力上，進而利用特征值法及能量法對葉片運動穩定性進行判定。由于受到計算能力的限制，先后發展了線性氣動力模型、全勢流模型、跨聲速小擾動模型以及升力面模型等用于葉片表面氣動力預測，從而利用特征值方法或者能量法對葉片氣動彈性穩定性進行判斷。直到20世紀80年代開始利用CFD方法進行流固耦合求解，先后求解Euler方程、RANS方程，但大多數均是弱耦合方法，在此階段為了降低計算難度，Hall等[2-3]發展了處理葉輪機氣動彈性穩定性降階模型 (Reduced Order Models)，使 CFD 在氣動彈性穩定性領域開始占據極其重要的地位。

隨著商用航空發動機追求更低污染物排放，軍用航空發動機朝向更高推重比目標發展，民用航空發動機環形燃燒室和軍用航空發動機加力燃燒室均面臨著愈發復雜凸顯的燃燒不穩定性問題。燃燒不穩定性問題往往發生在高溫高壓環境下，涉及燃燒化學反應、湍流流動及復雜環境聲學等多時空尺度問題，通過實驗和高精度數值方法能夠清晰地捕捉局部火焰反應、火焰和流動干涉、火焰對聲波響應等局部關鍵機理信息，但是高成本的實驗、高難度的測試、高成本的計算資源決定了不能完全依賴這些方法用于重復性計算的參數優化設計來支撐燃燒室的研制。因此，通過模型方法對燃燒不穩定性發生頻率、模態進行準確預測，并發展可靠的控制手段成為了重要研究方向。

諸多物理和工程問題在數學上都歸結為求解矩陣的特征值和特征向量的問題，這主要是由于所研究的對象往往都是以動力學微分方程組作為控制方程來描述的，并且由于非線性模型通常可以被近似為線性模型，因此，求解代數微分方程的解析解和數值解就成為了自然科學乃至社會科學問題模型研究的必要手段，相應的研究對象也逐步由代數方程轉移到矩陣本身，依靠矩陣相似變換來求矩陣的特征值與特征向量，并以此來反映線性系統的某些物理特征。特征值理論在求解多種穩定性問題上得到了應用，人們在求解諸如Poiseuille流動和Couette流動的穩定性和轉捩問題時采用的方法是線性全局穩定性分析方法[4]，假定我們可以精確已知或者估計得到Navier-Stokes方程的定常流動解，如果向這個流動解中引入小擾動，那么線化小擾動方程以及相關的邊界條件就是齊次的且與exp(σt)相關，當擾動方程的解存在且Re(σ)>0時，流場就是線化不穩定的。以此為基礎，通過小擾動假設和特征值理論，可以研究多種葉輪機領域常見的穩定性問題，因此，本文將簡要綜述特征值理論在航空發動機流動、氣彈和燃燒穩定性預測方面的應用研究現狀，并給出筆者所在課題組在相關理論建模的分析與效果驗證。

1 基于特征值理論的壓氣機流動穩定性預測

1.1 線化小擾動穩定性模型

自然界中存在多種流動不穩定性問題，比如與重力和加熱相關的穩定性問題、分層流穩定性問題、與離心慣性力有關的穩定性問題、與表面張力有關的穩定性問題、湍流轉捩問題。各類穩定性問題由不同的物理機制所支配。研究流動穩定性問題的方法主要分為線性穩定性理論和非線性穩定性理論兩種。線性穩定性理論是從無限小擾動的假設出發，流場中的各種變量被分解為平均量和擾動量兩部分，并且忽略了小擾動的高階項，即不考慮擾動之間的相互作用，從而使方程線性化。以線性化方程為基礎建立的理論即為線性穩定性理論。線性穩定性理論主要關注的是流體失穩的起始階段，而不能用來描述失穩后轉捩的全過程，因為這個過程主要是非線性效應起決定性作用的階段。但線性理論可以說明哪種速度剖面在擾動作用下會產生失穩，哪些頻率的振動增長最快等，這些對我們預測失穩的發生有著重要意義。而非線性穩定性理論是研究有限振幅擾動波的情況，這時擾動速度與主流平均速度比較起來就不再是很小的量，因此在線化過程中省略掉的擾動量的二階項不能再被忽略，小擾動方法所計算的結果逐漸偏離實際擾動發展的過程，非線性效應開始顯現。按照線性穩定性理論，一旦失穩開始后，擾動應該按照指數關系迅速增長，而實際情況決不可能如此。

實際上，對于湍流轉捩問題而言，流動失穩往往是轉捩過程的開始，而并非其全部過程，穩定性理論也不能描述其全過程。然而，目前沒有任何一種完整的理論分析能解析地描述轉捩的全過程。雖然線性穩定性理論不能用于湍流轉捩過程的非線性階段，但卻可以說明哪種速度剖面最不穩定，哪些頻率波動增長最快。這和葉輪機械內部的流動失穩問題十分類似，又略有不同。

葉輪機械內部流動的不穩定現象屬于流動不穩定性問題中的一種。同時，作為旋轉機械，其變工況特征、復雜的幾何形狀和固壁邊界使這一問題更加棘手。迄今為止，葉輪機流動穩定性理論預測研究方法大致可以分為兩類：數值計算方法[7-15]和半解析半數值的理論模型方法[16-27]。以Emmons模型[16]和Stenning模型[17]為代表的線化小擾動理論模型、以Greitzer模型[21-22]和Moore-Greitzer (M-G模型[23-24])為代表的系統穩定性模型均是典型的半解析半數值理論模型，其優點在于方法判據清晰、求解快速，適合處理多級壓縮系統失速起始判定問題，線化小擾動穩定性模型基于的就是二維不可壓小擾動方法，這里不妨展開簡要介紹：1955年Emmons[16]采用二維不可壓小擾動方法，給出了葉片排出口氣流堵塞面積隨攻角的變化，建立了壓氣機失速預測模型，并給出了旋轉失速物理現象的經典描述。Nenni和Ludwig[19]建立了能夠考慮單排和雙排葉柵影響的二維不可壓旋轉失速穩定性模型，同樣給出了穩定性判定的解析表達式，并指出總壓損失系數及其導數是影響壓氣機穩定性的最重要參數。

而系統穩定性模型基于彈簧、質量、阻尼系統類比，Greitzer建立了最初的一維系統穩定性模型，后來Moore和Greitzer將該模型(M-G模型[23-24])推廣到二維不可壓流，能夠進行喘振、失速的穩定性分析。M-G模型采用了大量簡化處理，將整個壓縮系統作為激盤，無法包含壓氣機復雜幾何造型等因素對穩定性的影響效果，難以適用于具有明顯三維特征的旋轉失速問題，并且其假設壓氣機都在靜壓升系數斜率為零的條件下失穩，用B參數作為判斷失穩模式的關鍵參數，這不是本文重點討論的內容，因此不再贅述。

采用小擾動方法分析葉輪機流動穩定性問題在很大程度上完全符合研究對象的物理特征，圖2給出某單級風扇失速瞬態壁面靜壓信號空間傅里葉變換幅值曲線[28]，可以看出：擾動波在初始階段的幅值很小，隨后開始逐漸增長直至某一時刻突然非線性放大進入失穩狀態。因此，若需要研究失速演化的過程，或者對于某些大尺度擾動(如進口流場畸變引起的流動不穩定性)對流場的影響，可以采用非線性穩定性理論或者數值模擬的方式。然而，如果我們只關注失速起始階段的系統穩定特性，完全可以采用線性小擾動理論來預測失速的產生。

壓氣機的流動失穩過程尤其是在失穩起始階段，通常被認為是擾動演化的開始，此時的擾動幅值較小滿足線性假設，通過線性穩定性分析的方式對小擾動的演化趨勢進行分析判別可以在一定程度上判斷流動的穩定趨勢。在穩定性模型發展的早期階段，經典的線化模型在揭示旋轉失速發生機理層面具有重要的意義，但是大量的諸如不可壓、無黏、無旋和均勻背景流動等簡化假設決定了這些將葉片簡化為激盤或者半激盤的集中參數模型無法包含流動非均勻性和葉片力的影響。

為了克服早期經典穩定性模型的這一缺陷，人們開始利用數值線性穩定性分析方法開展失穩預測的嘗試，美國麻省理工學院Gordon提出了一個三維不可壓數值穩定性模型[25]，在基本流軸對稱的假設基礎上，用體積力模型替代葉片對流體的作用，包含了子午流面內流動非均勻性。不過，該模型采用的葉片力模型較為粗糙，并且不能包含壓縮性的影響，其預測精度和使用范圍受到了限制。然而，針對現代高負荷高速壓氣機而言，發展三維可壓縮穩定性模型勢在必行。1996年，孫曉峰[26]從三維可壓縮線化非定常Euler方程出發，通過與研究流動穩定性所用的三維O-S方程的相應工作類比，采用無窮維系統的截斷技術以及有限空間波傳播界面間的模態匹配方法，發展了可以考慮任意階徑向擾動的三維可壓縮旋轉失速穩定性模型。在此工作基礎上，本課題組[27-28]結合等價分布源法和傳遞單元法，將激波匹配關系成功引入半激盤模型中，發展了跨聲速風扇/壓氣機穩定性預測與擴穩理論模型。該模型不關心擾動量的具體數值，因為擾動量的大小是隨時間變化的。而失穩是因為擾動量隨時間放大所導致，因此該模型只根據擾動量隨時間放大還是衰減來判斷系統穩定與否。為了更精確地反應流場內的物理畫面，使之包含更多的流動細節(尤其是小尺度的擾動)，葉片造型和黏性的影響就要在計算模型中得到考慮，因此求解非定常雷諾平均Navier-Stokes方程的方法開始被引入到壓氣機穩定性模型中來。

一個動力學系統的穩定性取決于其自身對于系統內部或外部的任何小擾動的響應，而并非針對某一特定的擾動。對于大部分流動穩定性問題而言，擾動總是經歷由小到大的發展過程。如果我們更關注流動失穩的開始階段，那么采用線化處理是可行且合理的。這也就是所謂的理論上任何流動失穩開始問題都可以描述為一個基于Navier-Stokes方程的特征值問題。

該預測方法的有效性在Stage37和NASA兩級風扇上得到了驗證[28]，Stage37在設計轉速下葉尖馬赫數為1.4，遠遠超出了亞聲速模型的假定范圍，因此，亞模型不適合該轉速。即使90%的設計轉速下，葉尖部分的大部分流動也是跨聲速，只有在70%的設計轉速下，葉尖的進口相對馬赫數才小于1，因此可以試用文獻[27]所發展的模型來預測其失速起始點。

模型預測的結果見圖3，圖中有兩簇曲線，一簇對應延遲時間為0的情況，一簇對應取正常延遲時間的情況；延遲時間在高亞聲速壓氣機中的作用相比低速壓氣機要明顯的多，從圖3中可以看出，加上延遲時間后，中性穩定點要略向高流量點偏移，高階模態下的偏移要大一些；而延遲時間對相對周向傳播速度的影響似乎要更明顯一些，由此可見，該亞聲速穩定性模型在高亞聲速條件下也能夠做出合乎實際的預測。

NASA兩級風扇在第1級轉子采用小展弦比設計，目的是消除原有大展弦比風扇葉片的弊端，改善通道流場，提高級壓比，減少葉片數目，并用于研究小展弦比葉片的優越性。此外其相對充分的試驗數據，也可以用來評估模型對多級壓氣機穩定性失速先兆點的預測能力。其中，NASA兩級風扇失速實驗通過在不同的展向位置同時采集徑向11處不同點的數據來研究設計轉速從50%～100%的失速邊界。在完全設計轉速情況下，質量流量為33.248 kg/s，總壓力比為2.399，絕熱效率為0.849(圖4)。

圖4的結果顯示了模型對NASA兩級風扇100%設計轉速下的穩定性預測。當衰減因子趨于0，壓縮機系統的穩定性由穩定狀態變為不穩定狀態，從而得到失速點，通過分析可以看出理論預測和實驗結果之間的誤差約5%。

1.2 穩定性預測通用理論

考慮到解析模型和數值模型各自的特點，近年來，一種新的建立穩定性模型的方式得到了發展[29-30]，一方面像解析模型那樣，能夠為判斷失速邊界或者失穩點提供清晰的判斷依據或準則，同時像數值模型那樣，能夠包含葉片造型對流場結構的影響，并且不像解析模型那樣嚴重依賴經驗關系。這樣一種方式，可以提高模型預測的準確性，同時在設計階段為設計人員提供一種更好的理論預測工具，考慮葉片彎掠、流道結構等復雜造型參數對于穩定性的定量影響，這是以往任何模型所難以比擬的?？梢栽O想，通過應用這樣一種在數值解法上可以快速實現、物理層面上可以包含更多因素、失穩點判斷清晰明確的模型方法，使其作為穩定性設計而成為壓氣機設計階段中的關鍵一環，而這樣一環是當前葉輪機械設計體系中所欠缺的。為了更好地解決風扇/壓氣機流動穩定性帶來的一系列問題，有必要將穩定性設計納入到現代風扇/壓氣機設計體系中(圖5[29-30])。葉輪機內流通用穩定性模型就是針對這一體系所做的探索實踐。

類比層流到滯流的轉捩過程問題的處理方法[31-32]，將壓氣機流場的定常數值計算結果作為某一流動狀態的平均流場，并在此基礎上引入一個包含所有物理量的小擾動，通過計算將壓氣機流動作為動力學系統的特征頻率來判斷其流動穩定性。對于任意的離心或是軸流壓氣機系統，從最基本的Navier-Stokes方程出發，應用線化小擾動理論，結合浸入式邊界方法和體積力方程，發展通用的流動穩定性模型，該模型可以用于任意進出口條件和邊界情況下的復雜流動失穩起始點的預測。

該方法的具體操作流程為：從控制方程出發——CFD主流流場——體積力模型——建立特征值問題——求解特征值問題。

根據浸入式邊界理論，葉輪機葉片用分布的力源項來代替，從而，葉輪機內部流動由帶力源項的Navier-Stokes方程來控制，即：

質量守恒方程：

(1)

式中：ρ為密度；u為速度矢量。

動量守恒方程：

(2)

式中：F為單位質量流量的葉片力矢量；Π為二階應力張量。

能量守恒方程：

(3)

式中：λ為導熱系數；WF表示葉片力做功。

如1.1節所述，發展葉輪機流動通用穩定性理論的核心目的是進行失速起始點的預測，因此，基于小擾動假設，把原始三維非定常流動分解成定?；玖骱头嵌ǔＰ_動之和：

(4)

根據多元函數泰勒展開，葉片力可以線化為

(5)

將式(4)、式(5)代入式(1)～式(3)，減掉基本流所滿足的零階方程并且忽略高階無窮小量，則可得到線化形式的Navier-Stokes方程。將線化Navier-Stokes方程在圓柱坐標系下展開，可得：

(6)

式中：Φ=[ρ′u′v′w′p′]T，而A、B、C、E、G、H、M、N、Q、R以及S是由基本流所確定的系數矩陣。

一般地，葉輪機內部基本流在徑向、周向以及軸向都是非均勻的，假定基本流是定常的，因此，只能對小擾動量作關于時間的正則展開，即：

(7)

將式(7)代入到式(6)，可以得到以復頻率ω為特征值的特征值方程：

(8)

根據齊次線性偏微分方程理論，式(8)要有非零解，則必須滿足：

det[L(ω)]=0

(9)

通過求解式(9)即可得到復頻率ω=ωr+iωi，復頻率的實部ωr表征失速先兆波的頻率，而復頻率的虛部ωi則表征小擾動的演化趨勢：ωi為正，則擾動隨時間放大，系統不穩定；ωi為負，則擾動隨時間衰減，系統穩定。

考慮到葉片力建模以及求解特征方程的復雜性，根據實際工程的要求以及當前的計算能力、設計體系，需要對穩定性理論進行不同層次的簡化，得到不同層次的簡化模型。當馬赫數較低情況下，壓縮性對于流動本身影響不大，可以將流動主控方程簡化為三維不可壓縮Navier-Stokes方程組；對于近似軸對稱流動情況，可以重點關注子午平均流面；當徑向方向流動可以忽略時，研究二維軸對稱流動便能抓住問題的本質和關鍵；為了進一步簡化計算量，準二維流線也可以作為研究主體。這些簡化方式要根據實際問題的需要而適當選擇，所發展的葉輪機流動穩定性理論模型可被用于軸流壓氣機和離心壓氣機的流動穩定性預測，并分析探討氣流壓縮性、葉尖間隙及葉片造型對壓氣機穩定性的影響。

該方法在中國科學院工程熱物理研究所的低速軸流壓氣機實驗臺和Rotor37穩定性預測上得到了應用驗證。圖6針對前者低速壓氣機的計算結果顯示：隨著壓氣機截流過程，流量逐漸減小，所有9個模態的衰減因子不斷增長。周向傳播速度較慢的前5個模態比后4個模態更快地趨于衰減因子臨界值。模態1的衰減因子最先由負變正，代表了壓氣機系統由穩定變為不穩定，因而模態1為最不穩定模態。臨界點的質量流量為2.324 kg/s，與實驗測得的失穩開始點質量流量的相對誤差為0.65%。

Rotor37采用跨聲速葉型設計，據了解，目前國際上公開發表的文獻中鮮有研究一般壓氣機結構下跨聲速流動失速起始的模型預測。因為在跨聲速流場中，葉柵中的激波對流動有重要的影響，使得很難在壓氣機系統中對其模化。該方法采用周向平均的定常流場分布，葉柵通道內的激波所導致的損失和壓升都能包含進入平均流場中，并且經過平均處理后，通道內也不存在流場的跳躍式突變。因此，本模型有可能克服激波這一突躍條件帶來的挑戰。圖7為Rotor37的穩定性預測結果，并在此基礎上發展了可以考慮葉片彎掠造型對穩定性影響的預測分析方法。圖中：NTD代表無量綱延遲時間。

最近幾年，本課題組選擇Rotor37為基準葉型，通過沿弦向移動基元葉型來獲得一系列的掠葉片，進而采用子午面模型計算預測穩定裕度，并結合定常流場的計算結果，分析討論氣動掠設計對壓氣機穩定性的影響[33]。

如圖8所示[33]，前掠并沒有對轉子的穩定性造成顯著的改變。然而，后掠轉子相比Rotor37穩定裕度有明顯降低。進一步而言，后掠程度越大，穩定裕度隨之單調遞減。這一結果與過去的研究定性上是吻合的[34-36]。圖中：BS代表Backward Sweep;FS代表Forward Sweep。

以北京航空航天大學單級低速壓氣機實驗臺TA36為算例，通過在流線曲率通流設計中調整葉片環量分布設計得到不同負荷分布的轉子葉片。采用子午面模型預測其失穩邊界，結合定常和非定常流場結果，分析討論轉子葉片軸向負荷分布對壓氣機穩定性的影響[37]。

對3種負荷分布轉子的壓氣機進行失穩邊界預測的結果展示在圖9中[37]?？梢钥吹?，隨著節流過程進行，轉子前加載壓氣機最先失穩，其次為均勻加載壓氣機，轉子后加載壓氣機最后失穩。這表明，該算例中的轉子軸向負荷后移對流動穩定性而言是有利的。

2 基于特征值理論的壓氣機氣彈穩定性預測

航空領域的氣動彈性(簡稱氣彈)穩定性問題涉及到復雜的流動和結構相互作用，預測這一問題需要分別對流動以及結構建立相應的模型，然后將二者結合來對氣動彈性穩定性進行預測。Theodorsen[38]首次給出了采用特征值理論對氣動彈性穩定性進行判斷的方法，該方法采用勢流理論結合庫塔條件描述葉片所受到的流動作用力，同時采用剛體運動方程來描述翼型的振動，聯立二者可建立葉片的動力學方程，通過對方程進行特征值分析來判斷翼型在特定工況下是否會出現氣動彈性穩定性問題。在今天看來，雖然Theodorsen的工作采用了極其簡化的模型來描述流動以及結構的動力學特性，但是其采用特征值方法對于氣動彈性穩定性進行判斷的思想卻對隨后的研究工作具有深刻的啟示意義。

Dowell[39]系統總結了對于單個孤立翼型的特征值分析方法，除了對氣動力進行理論建模以外，還可以通過實驗數據對氣動力進行擬合。Whitehead[40]曾對振蕩葉柵的升力系數進行測量以此作為葉片動力學方程中的輸入來對其氣彈穩定性進行特征值分析。Bendiksen和Friedmann[41]采用特征值方法對一個二維葉柵彎扭復合顫振進行了穩定性分析，Kielb和Kaza[42]通過采用二維梁以及有限元模型來對葉片構建結構模型，結合特征值方法對帶掠風扇葉片的顫振穩定性進行了分析，二者對于不可壓縮流動的描述均采用Whitehead[40]所給出振蕩葉柵氣動力模型。Kileb 和Ramsey[43]則采用基于拉普拉斯變換的線化氣動力模型對葉柵氣動載荷進行建模，通過特征值方法給出對超聲速流動下葉柵的氣彈穩定性變化趨勢。

在早期計算能力不發達的年代，進行特征值分析的前提是需要將氣動作用力解析表達成頻域形式，而對于壓氣機內的相關問題，葉片氣動載荷的建模多取自于外流孤立翼型中相關的研究結果，這也導致理論結果與真實現象存在較大的出入[44]。壓氣機葉片由合金制造，質量比較大，此時顫振通常表現為單一的模態，而外流機翼中多為復合模態；另一方面，壓氣機內部的流動受到葉片和流道復雜幾何、固壁以及上下游結構的影響，其復雜程度遠高于外流問題，并且表現出明顯的流固耦合特性，對葉片氣動作用力進行準確的理論建模難度較大，因此適用于外流氣彈分析的理論解耦模型通常難以準確描述壓氣機內部的真實流動情況。雖然部分基于解析方法的預測模型被用于設計階段的快速迭代，但是通常認為這一類方法預測的結果與實際問題有著較大差別。

Carta[45]提出了一種能量法的判據來預測壓氣機葉片是否會出現顫振問題，這一方法所采用的氣動力及葉片結構模型與Theodorsen的工作是類似的，但是區別在于Carta并沒有對控制方程的特征值進行求解，而是只提取了結構方程部分的特征頻率，在此基礎上假設葉片在特征頻率處做簡諧運動，這一做法與Whitehead[46-47]的工作中對流體域與固體域做解耦處理是一致的，通過計算一個周期內氣動作用力對葉片所做的功來判斷葉片與流動之間的能量傳遞情況，即氣動阻尼功，如果氣動力對葉片輸入負功，則認為流動對葉片的振動起到阻尼作用，此時葉片處于穩定狀態，否則即可能發生顫振。這一基于能量法的顫振判斷準則在壓氣機的氣動彈性問題預測中一直沿用至今。

2.1 基于能量法的快速氣彈預測算法

能量法通過指定葉片振動規律來對氣彈問題進行解耦，其核心問題在于如何對流動的非定常響應進行預測。在壓氣機設計階段，對于預測方法的主要要求便是能夠實現高效迭代以完成快速設計，此時流場對于葉片非定常輸入的響應多采用簡化的解析方法來建模。Sun等[48]采用三維升力面理論來求解由葉片振動誘導產生的非定常載荷響應：在葉片表面建立無穿透邊界條件，要求來流擾動速度場與葉片振動速度在葉片表面法向和速度為零，通過將速度脈動轉化為壓力脈動的積分方程，如圖10所示，結合廣義聲類比積分公式，可將轉子散射的壓力脈動場寫為

(10)

假設壓力脈動存在時空周期性，則葉片上下表面壓差可展開成傅里葉級數的形式：

(11)

進一步結合格林函數，散射的壓力脈動場可寫為積分方程的形式：

(12)

近年來，部分研究結果表明轉子進口段的長度以及聲處理可能會對轉子的顫振穩定性產生較大的影響[48-49]，為了考慮進口段長度以及聲襯邊界條件的影響，Sun等[48]結合傳遞單元方法，通過等價分布源方法來描述聲學邊界對于葉片顫振的非定常響應，同時采用邊界元方法來刻畫管道邊界的反射效應，成功建立了有限/無限長管道內的流固聲耦合的預測模型。

圖11給出了不同聲學邊界條件下轉子葉片氣動功Waero的變化結果，葉片振動呈現為一階彎扭復合模態，當振動節徑ND=2時，氣動功的剛壁值小于0，即轉子在剛壁管道中是穩定的，而當阻抗為Zs=(0.05,-3.9)的聲襯被加入管道后，氣動功進一步下降為一個絕對值比原剛壁值大兩倍的負數，這也意味著氣流對轉子的阻尼效應增強。而在同一條件下，若聲襯阻抗為Zs=(0.05,-3.2)，則聲襯反而會對振蕩轉子的氣動彈性響應產生負面影響，觸發葉片產生顫振。這一對比說明，壁面聲處理對于轉子顫振穩定性產生重要的影響。

與此同時，轉子進口管道的長度也可能對振蕩轉子的氣動阻尼產生重要的影響，圖12給出了不同管道長度下，進口距葉片前緣面距離長度對于振蕩轉子氣動功的影響結果。圖中：FLD和LID分別表示有限長和無限長；Lduct為管道長度；Ca為葉片軸向弦長。從圖中可以看到，基于有限長和無限長管道模型所得到的非定常氣動功存在明顯差異。當管道總長度為5倍葉片軸向弦長時，如圖12(a)所示，非定常氣動功會隨著進口段長度增加而由負值逐漸變為正，這也意味著進口段過長可能導致轉子出現顫振。而在圖12(b)所給出的管道長度為15倍葉片軸向弦長的情形下，當轉子位于管道中間位置時，非定常氣動功的數值最小，此時轉子葉片最穩定，而當轉子逐漸靠近進口或者出口位置時，非定常氣動功的數值呈現出逐漸增大的趨勢直至變為正，這也意味著對于圖中所考慮的振動模態，轉子過于靠近管道進口或者出口可能導致顫振的發生。

2.2 氣流非定常載荷預測方法

近年來，隨著計算方法以及軟硬件的發展，基于CFD方法求解更為復雜的非線性流動模型在氣流非定常載荷的預測中扮演著越來越重要的角色。在20世紀80～90年代，受制于計算能力的限制，Hall和Crawley[50]通過對非線性歐拉方程進行線化的方式來提升計算效率。線化后的歐拉方程可以在特征頻率處做周期性假設，對于氣彈問題，這一頻率通常是葉片的振動頻率，進而可以求得該特征頻率所對應的空間特征向量，特征向量的物理意義是周期性的流動在空間各處脈動的幅值。計算效率是評估壓氣機氣彈預測方法的核心指標之一，而基于特征頻率的處理手段可以避免對時間項進行迭代求解，因而可以大幅提升計算效率。當流場中不存在強非線性結構時，比如流動分離以及激波，基于線化方程的計算方法能夠對由葉片振動引起的非定?，F象進行合理的預測，而當流場中非線性作用占據主導時，頻域線化方法則不再適用。

為了進一步考慮包含非線性的影響同時兼顧頻域方法在計算效率上的優勢，He和Ning[51]提出了一種非線性諧波方法，這一方法的核心思想是在對Navier-Stokes方程進行線化的過程中保留擾動交叉項，以特征頻率對擾動進行周期性假設，在對擾動求解的過程中同時對包含非線性因素的時均流場進行同步修正直至結果收斂。相較于非線性時間推進，非線性諧波方法的計算效率仍然有1～2個數量級的提升?；陬愃频乃枷?，Hall等[52]提出了一種求解具有時空周期性非線性流場的諧波平衡方法，這一方法通過迭代求解特征頻率及其各階諧波所對應的空間特征分布來對流場進行時空重構，能夠實現較高的非定常計算效率。非線性諧波法以及諧波平衡方法是目前研究壓氣機葉片氣彈問題的主要模擬工具。Clark[53]采用諧波平衡方法對壓氣機葉片的流動進行快速模擬，在此基礎上采用POD(Proper Orthogonal Decomposition)方法提取流場中的主要模態來構建降階模型并嘗試將其應用于轉子非同步振動問題的預測。Wang和Huang[54]提出了一種針對諧波平衡方法的加速收斂方法并將其應用在跨聲速轉子的氣彈穩定性預測。Huang等[55]采用諧波平衡方法結合影響系數法求解了葉片振動模態處的氣動作用力，在此基礎上結合僅包含指定模態的結構動力學方程對葉片的顫振實現基于小擾動的特征值分析?；谥C波平衡的思想，Du和Ning[56]提出一種采用多項式來簡化時間項的方式，同樣能夠準確地對具有時空周期性的流動進行預測。Sandberg和Michelassi[57]與Casoni和Benini[58]針對軸流葉輪機械領域典型的數值模型以及頻域和其他各類降階模型給出了全面的總結。

2.3 耦合計算的應用場景

采用時空周期性假設的求解方法一般要利用特征頻率對時間項做傅里葉展開，而當流場中同時包含多個特征頻率時，如多葉片排干涉上下游葉片通過頻率不同，以及轉靜干涉的同時還存在葉片振動，這一類方法通常只能對不同頻率做線性疊加而無法考慮彼此之間的非線性耦合效應。當不同頻率線性疊加，多個頻率互相靠近產生“拍”現象時，計算收斂時間也會大幅延長[59]。這一類方法實際上是解耦的，即沒有考慮流場演化對于葉片振動頻率、振幅以及葉片間相角的影響。由于目前復合材料的大量應用，葉片質量相較于早期已經大幅降低，此時葉片振動頻率及模態與不考慮流固耦合效應的結果會存在差異，與此同時，葉片振動響應可能包含多個模態，采用頻率解耦方法不符合真實的物理演化過程。Chahine等[60]的計算結果表明，隨著質量比以及結構剛性的降低，耦合效應的作用越來越顯著，由此導致耦合方法與解耦的能量法所給出的結果差異越來越大。在流固耦合效應起主導作用時，基于全非線性方程的耦合計算則顯得十分必要[61-62]。

另一類必須采用耦合模型計算的振動問題是轉子非同步振動問題。近年來國外研究機構在不同的實驗平臺上都曾測量到了葉片非同步振動的信號[63-65]，這一現象一般呈現出高周向模態數，并且伴隨著沿周向傳播的分離流動結構，同時壁面壓力信號存在多個非軸頻整數倍的異常成分。Sun[48]和Vahdati[49]等的工作都表明，進口管道的反射會對轉子風扇的非同步振動穩定性產生明顯的影響。對于這一現象必須采用耦合處理的原因在于，從實驗結果來看，受到葉片振動與流動之間干涉效應的影響，流動的周向模態特性在振動前后會出現明顯的改變[64-65]，即流固耦合效應在這一現象中同樣扮演著重要的角色。目前針對非同步振動問題已有較多的研究結果，在判斷葉片非同步振動的穩定性時，絕大多數工作[66-69]仍然采用解耦的處理方式，即給定頻率和振幅并且使用氣動阻尼這一參數。氣動阻尼這一概念源于線性范疇，需要參照給定的葉片振幅做無量綱從而保證氣動阻尼與振幅無關而非同步振動往往伴隨著沿周向傳播的分離流動，表現出強非線性特征，在此情形下，振幅對于氣動阻尼計算結果的影響目前尚不明確。雖然葉片的振動頻率可能仍然是固有頻率，但是非同步振動的周向模態特性與流動密切相關，同樣無法在計算前獲得，并且葉片的振動會顯著改變流動的周向模態特性，這也意味著無法準確地給出周期性邊界條件，因此大多數情況下需要開展全環計算。較于上述基于頻域假設的快速非定常求解算法，采用時間推進方式求解完全非線性的Navier-Stokes方程能夠包含更多復雜流動因素，因此也被廣泛應用于壓氣機氣彈問題預測，同時其對于計算量的要求也進一步提升。

3 基于特征值理論的燃燒不穩定性預測

燃燒不穩定性的發生是燃燒系統內火焰非定常熱釋放和聲波的充分耦合導致的[70]。因此，燃燒不穩定性的發生關鍵在于火焰非定常熱釋放和聲波的耦合，將復雜問題降階簡化，保留影響燃燒穩定性發生的關鍵因素，即火焰對聲波的響應，以及影響聲波反射的燃燒室的阻抗壁面刻畫，在其他部分忽略次要因素，比如整個燃燒室內復雜的湍流流動、燃燒反應等，而只考慮聲波傳播[71]。通過模型實驗和高精度數值模擬可以得到火焰對聲波的響應，在線性頻率內得到火焰傳遞函數(Flame Transfer Function, FTF)[72]，在非線性頻域范圍內，Dowling等發現燃燒不穩定性非線性的主要因素在于火焰對不同幅值聲波響應的不同，因此可以通過建立包含聲波幅值影響的火焰函數，也稱火焰描述函數(Flame Describing Function)[73-74]用于非線性燃燒不穩性分析。至于燃燒室其他部分則通過建立聲學模型，結合火焰傳遞函數在不同界面上進行守恒方程匹配。

燃燒不穩定性分析本質上是基于小擾動特征值理論的穩定性分析，即分析燃燒系統有燃油當量比脈動、旋渦脫落及燃燒噪聲等小擾動存在的情況下，系統的響應特征，如果在施加起始小擾動后的一段時間內，系統的脈動幅值不斷放大，則表明燃燒系統是線性不穩定的，反之，則是線性穩定的。基于這樣的認識，我們可以對系統控制方程進行線化處理，得到關于小擾動量的控制方程，并求解其復特征頻率，其中復特征頻率的實部表示系統的熱聲共振頻率，虛部則代表系統不穩定性的增長率，通過這種理論模型方法，可以快速進行燃燒不穩定性的預測或研究控制手段[75-76]。

國際上通用的聲網絡低階模型方法在認識燃燒不穩定性發生機理及對燃燒室進行穩定性預測方面發揮了重要作用，能夠為燃燒不穩定性控制提供準確的目標，不穩定模態及頻率。在考慮控制手段設計時，以上模型方法能夠直接在特征頻率分析中考慮具有集總參數特點的Helmholtz共振器[77-78]。但是，當需要對燃燒室中常常使用的兼具冷卻壁面和提供聲學耗散作用的含冷卻偏流穿孔板聲襯[79]等具有分布參數特點的控制壁面設計時，需要配合額外聲學計算[80]或實驗[81]。但是這種根據頻率、模態信息再通過聲學計算或實驗解耦地設計燃燒不穩定性抑制手段往往會帶來誤差。這是因為，首先，穿孔板聲襯等聲學軟壁面在燃燒室內會改變系統的熱聲共振頻率和模態分布，因此，在硬壁面條件下建立的預測模型不能準確預測帶有聲襯的燃燒室共振頻率和模態分布。另外，通過聲學計算或實驗設計聲襯過程中，聲襯面臨的聲學環境也很難完全實現其在燃燒室內的環境。因此，為了研究聲襯控制熱聲不穩定性效果，并準確地對聲襯參數進行優化設計，建立一種能夠將聲襯耦合考慮的燃燒不穩定性模型就顯得十分重要而且必要了。

3.1 耦合考慮聲襯壁面的三維燃燒不穩定性預測模型

將聲襯耦合考慮進燃燒不穩定性模型中，難點在于如何描述燃燒室內聲波和壁面的相互干涉作用。在軟壁面環境下，壁面阻抗、系統頻率、管道特征值形成復雜的耦合迭代關系[82]。為了解決這一難題，Sun等提出一種基于等價分布源思想[83]的傳遞單元方法[84]，將聲襯穿孔板阻抗邊界視為等效連續分布的單極子聲源，這樣能夠保留管道特征值的正交特性，避免軟壁面管道徑向特征值的復雜迭代求解。傳遞單元方法，可以用以考慮有限長聲襯的聲傳播問題，也可以用以將聲襯作為聲學傳遞單元包含到燃燒不穩定性模型。

該方法在處理聲襯段聲學單元時，為了可以將之耦合到燃燒不穩定性分析模型中，需要得到其傳遞矩陣，即聲襯段內聲波關于界面入射波的顯示表達式。具體地，將聲襯軟壁面視為等價連續分布的單極子聲源，聲襯段的內壓力擾動和速度擾動表達式可以寫為界面入射聲波和散射聲波之和：

p′=p′I+p′d

(13)

其中：

(14)

聲襯背腔管道內由聲襯單極子聲源產生的散射聲場為

(15)

p′c-p′=Zv′n

(16)

式中：Z為穿孔板聲阻抗；v′n為垂直于聲襯表面的聲質點速度。通過阻抗方程以及聲襯前后端界面匹配條件即可求解得到等階聲質點速度的表達式，再將之代入方程(13)和方程(14)即可以得到聲襯段內聲場關于界面入射波的顯示表達式，即完成了聲襯段聲傳遞單元的建立詳細推導過程和表達式參考文獻[84]接著可以和燃燒室內其他硬壁面條件下聲網絡單元建立起燃燒不穩定性分析模型。關于模型建立的細節可以參考文獻[85]。應用該特征值分析模型，可以研究不同形式的聲軟壁面控制對燃燒不穩定性的影響，并進行快速的參數優化設計。

從Rayleigh準則出發可知，控制燃燒不穩定性的發生可以從兩個方面著手，一是削弱熱聲耦合強度[86]；二是增加系統的聲學耗散[82,85]。下文將從增加系統聲學耗散角度研究穿孔板聲襯在環形燃燒室中控制熱聲不穩定性來闡述本三維模型方法的應用。

3.2 穿孔板聲襯控制環形燃燒室熱聲不穩定性

環形燃燒室由于其能產生更均勻的燃燒，具有結構更緊湊的優點，目前是航空發動機和地面燃氣輪機領域最為常用的燃燒室類型。為了滿足日益嚴苛的排放標準，環形燃燒室多采用貧油預混燃燒，具有較低的NOx(氮氧化物)排放，但更易出現燃燒不穩定性。此外，由于燃燒溫度較高，燃燒室火焰筒壁面往往開有小孔，從高壓壓氣機后引入部分冷氣從小孔內通過，形成小孔冷卻偏流，自然形成了具有聲學耗散的穿孔板聲襯結構。穿孔板聲襯在滿足冷卻需求的情況下，還能提供聲學耗散在線性范圍內抑制燃燒不穩定性的發生，或盡可能的降低振蕩幅值。目前，由于巨大的計算資源和時間成本，還難以通過數值模擬如大渦模擬(LES)研究不同聲襯結構的控制效果。而采用本文所介紹的包含軟壁面聲襯的燃燒不穩定性模型則可以進行聲襯控制效果的研究以及參數的優化。

通過與有限元分析軟件COMSOL對比無熱釋放條件下系統的特征頻率及模態分布情況，驗證本模型的有效性，頻率對比如圖14所示。圖中模態名稱中，L表示軸向模態，A表示周向模態，字母前數字表示模態階數，如2L0A代表第2階軸向模態混合0階周向模態。從圖中計算結果還可以發現，聲襯的引入使得必須考慮徑向耦合模態，才能準確預測頻率以及模態分布。也就是說，當我們考慮包含聲軟壁面的燃燒不穩性預測時，有必要發展三維模型，才能夠準確預測系統的頻率及增長率。

為了研究聲襯控制效果背后機理，給出系統內聲壓分布，對于模態2L0A聲襯控制效果最好的聲襯位置大小L2/Lc=0，L3/Lc=0.47，系統內聲壓分布如圖15(a)和圖15(b)所示。圖中：

以上研究結果表明，對于包含聲襯的燃燒室系統，有必要建立三維燃燒不穩定性模型通過特征值分析理論，預測系統的穩定性并研究聲襯的控制效果。聲襯的大小和位置對于其控制燃燒不穩定性效果有重要影響，可以通過改變聲襯背腔結構提高聲襯的控制效果。另外，值得一提的是，在實際系統中，由于火焰的非線性，還可以通過本模型方法結合非線性火焰函數模型[88]，以研究不同聲襯結構對降低熱聲不穩定性極限環幅值的影響，以滿足工程實際需求。

4 總結與展望

隨著航空發動機設計理論和三維數值模擬技術的不斷發展，針對發動機內流、氣彈和燃燒領域面臨的工程技術問題，各種理論模型和實驗數據庫也更加完善，設計人員普遍習慣于通過所積累的設計經驗來預估和判斷風扇/壓氣機以及主燃/加力燃燒室的穩定工作范圍，但基于經驗的評估方法很難被納入到一體化設計工具中并對設計方向給予合理的引導，因此，快速準確的穩定性評估工具對于航空發動機的設計至關重要。本文簡要闡述了基于小擾動方法的特征值理論在航空發動機氣動、氣彈和燃燒穩定性預測方面的應用，對從經典的線化小擾動模型到發展較為成熟的半解析模型方法的研究發展進行了綜述，并重點以本課題組近年來的研究結果為例，給出了該方法在流動失穩預測、葉輪機設計和燃燒不穩定性預測方面的應用案例。實踐表明，該方法可以作為發動機一體化設計分析工具成為發動機穩定性設計過程中的一環。

此外，鑒于壁面邊界條件可以影響小擾動在動力學系統內的演化行為，在本文所介紹的理論模型基礎上，相繼有通過壁面阻抗調控方法控制流動失穩和熱聲不穩定性的技術手段被發展，為失速喘振和燃燒不穩定性的控制提供了新的實現手段。面向未來更先進航空發動機研制的需求，基于特征值理論的穩定性模型方法需要考慮更加復雜的壁面邊界條件、幾何影響因素、更多非定常擾動干涉機制，以及非線性因素等，以快速低成本的方法優勢支撐行業發展可靠、高效的穩定性評估工具。