基于高效特征值分析方法的旋轉失速先兆預測

徐慎忍，何晨，孫大坤，袁蔡嘉，曹東明，趙家資，王丁喜

1．西北工業大學 動力與能源學院，西安 710072

2．中國空氣動力研究與發展中心，綿陽 621000

3．北京航空航天大學 能源與動力工程學院，北京 100191

旋轉失速起始點決定了壓氣機能夠穩定運行的最小流量。失速發生時，壓氣機中出現自激流動失穩現象，可導致氣動性能惡化并誘發流致振動，造成結構破壞。因此，過去的約80 年見證了失速機理、失速預測和失速控制等方面的大量試驗、理論和數值研究［1］。早期的一些試驗研究中發現旋轉失速發生前可在葉片前緣上游監測到幅值逐漸增加并最終導致旋轉失速甚至喘振的長波擾動，這一現象后來被稱為模態型旋轉失速［2-4］。麻省理工學院的Greitzer 提出了可用線性失穩理論來解釋旋轉失速先兆擾動［5］。該理論認為，在旋轉失速發生初期，壓氣機內流動以Hopf 分岔形式發生失穩。這一線性失穩理論定性解釋了很多早期試驗研究中觀察到的模態型失速現象。然而在1990 年前后針對負荷更高的現代壓氣機的諸多試驗測量中發現了一種與模態型失速截然不同的失速模式，此類失速發生時流場擾動的空間尺度遠小于模態型失速且發生的時間尺度更?。?］。具體而言，在流場中產生振幅與主流相當的劇烈波動前并未監測到模態型長波擾動。由于這類失速發生時所監測到的壓力脈動信號往往存在突尖波形，因此這種失速現象被稱為突尖型失速。這對經典的線性失穩理論提出了挑戰，因為該理論似乎無法解釋突尖型失速的觸發機理。

雖然學術界和工業界普遍認為突尖型失速不是由模態波引起的，然而這2 種現象并非毫無關聯。劍橋大學Day 在試驗中發現只要減小葉頂間隙即可實現從模態型到突尖型失速的切換［6］。由此可以猜測：突尖型失速是否亦是由模態波觸發的？雖然這一猜測由于在突尖型失速試驗中未觀察到模態波而很早就在文獻中被否定了，然而作者團隊認為這一判斷并無嚴謹的理論根據。對于在試驗中未在突尖信號前發現模態波信號，并不能排除測量儀器精度有限的原因。可以想象，模態型失速中，模態波幅值增長到某臨界值才會觸發大規模流動失穩并導致旋轉失速。該臨界值與壓氣機的負荷水平顯然具有很強的相關性。對于高負荷壓氣機，如果觸發流場大規模失穩所需的擾動閾值小于儀器測量精度，自然就無法在所謂的突尖失速發生前監測到模態波信號［6］。與風洞試驗相比，數值模擬中可通過監測某些特定位置的壓力脈動來實現數值傳感器的效果，其精度比物理傳感器高得多。例如，對于雙精度非定常RANS 求解器，10-10Pa 量級的壓力脈動可以很容易被監測到，而這在試驗測量中幾乎不可能。因此，數值模擬方法為研究突尖型失速是否亦由模態波擾動觸發提供了可能。

隨著數值算法和高性能計算（High Performance Computing， HPC）硬件的快速發展，基于RANS（Reynolds Averaged Navier-Stokes）方程的計算流體力學（Computational Fluid Dynamics，CFD）方法目前已被廣泛應用于壓氣機氣動設計和分析中。在旋轉失速研究中，傳統測量方法時空分辨率低，而基于RANS 方程的數值模擬方法很好地克服了這一弱點。牛津大學的He 在某亞聲速壓氣機環形葉柵［7］與He 和Ismael 在某跨聲速軸流壓氣機［8］上的非定常模擬結果顯示，非定常RANS 的CFD 手段可定性地模擬旋轉失速現象。為了進一步提高分析效率，劍橋大學的Brandvik 和Pullan 開發了可在圖形處理單元（Graphics Processing Unit， GPU）集群上運行的高效大規模非定常流場模擬程序TurboStream［9］，并基于此程序對NASA E3壓氣機中某轉子和劍橋大學某低速壓氣機進行了非定常模擬［10］。該研究一方面如實地再現了突尖型旋轉失速，提出突尖型旋轉失速是由從葉尖前緣吸力面延伸到機匣上的龍卷風狀渦流結構導致的；另一方面也指出突尖失速發生時伴隨著葉尖區域局部沖角過大的情況，從而引發前緣流動分離。并基于此現象，進一步提出了局部沖角過大而非葉頂間隙才是發生突尖型失速的必要條件。然而，盡管使用了高效的GPU 求解器，此類高精度非定常分析的計算成本仍然極高?；诋斍暗母咝阅苡嬎阌布Y源，只能對某些特定算例的個別工況進行分析，若要將其應用于參數化研究從而確定旋轉失速關鍵影響因素，依然顯得力不從心。除高昂的計算成本外，大多數旋轉失速非定常模擬都采用了增大個別葉片的安裝角以人為制造葉片失諧或施加外部擾動加速觸發失速的做法。這樣做的目的之一是將失速發展過程控制在一段相對較短的時間內，從而降低計算成本。否則，從非定常模擬開始到突尖型失速信號出現可能需要耗費極長的時間。然而，若要徹底弄清突尖失速是否由緩慢增長的模態波增長到一定程度觸發的，必須開展從穩定工況開始直到突尖失速發生之間整個物理過程的非定常模擬，并且此過程中不應引入葉片失諧或施加外界擾動，可想而知，相應的非定常模擬成本會更高。

本工作的主要目的是建立一種高效的穩定性分析方法，使其能夠以遠低于非定常模擬計算成本的代價研究無葉片失諧、無外部擾動時流場小擾動自發增長的過程，并基于此工具準確診斷所謂的突尖型失速是否亦由模態波觸發。由于本文工作的遠期目標是為了驗證模態型和突尖型失速是否都由線性失穩誘發，因此可引入小擾動假設，將非定常小擾動分析轉換為特征值問題，并發展適用于大規模流動問題的高效特征值分析方法?；谔卣髦档娜址€定性分析方法使得能夠以極低的計算成本對失速前的流動失穩細節進行精細化定量研究。

特征值方法已成功應用于外流穩定性研究，尤其是針對翼型、機翼甚至整個飛行器的激波抖振現象的研究［11-12］。外流中激波邊界層干擾所導致的抖振現象的根源近幾十年來一直懸而未決。直到最近幾年，針對其失穩機理的觀點才逐漸統一，即抖振本質上是Hopf 分岔線性失穩現象。而這個觀點也在時域非定常模擬［13-14］和特征值分析研究［12，15］中得以印證。與時域非定常模擬手段相比，特征值分析得到的特征值和特征向量可用于重構流場小擾動的時空演化，但其計算速度卻比前者快幾個數量級，并且在流場擾動增長至開始體現非線性之前能夠揭示與成本極高的非定常模擬同樣豐富的流動信息。此外，特征值分析結果還包含了激波抖振模態的空間結構，將其與非定常模擬結果或試驗結果進行對比，也可驗證計算所得特征值和特征向量的準確性。

基于特征值計算的穩定性分析方法在壓氣機流動穩定性研究中也已有應用?；诓豢蓧嚎s、無黏和軸對稱假設，通過引入忽略葉片復雜幾何外形但仍保留其對流動擾動作用的葉片力模型，流動失穩問題即可轉化為特征值問題［16］。類似的，北京航空航天大學孫曉峰團隊通過引入浸沒邊界法，并分別將葉片對流動轉折和黏滯的作用力視為源項代入到線化RANS 方程中，建立了基于特征值方法的可壓縮流動穩定性通用理論［17］。該方法已成功應用于軸流壓氣機［18］和離心壓氣機［19］失速研究中。此外，壓氣機葉片前掠和后掠對于失速起始點的影響也可以通過此基于特征值的穩定性分析方法進行研究［20］。然而，需要注意的是，由于在上述研究中引入了背景流場周向均勻的假設，因此不能解析流動在周向的空間分布。事實上，只有完全基于三維流動方程的特征向量才能夠為流動失穩研究提供與三維非定常流動模擬同樣豐富的信息。在充分考慮葉輪機械中幾何結構和流動循環對稱特征的基礎上，帝國理工大學的Schmid 等［21］提出了一種高效的基于三維流場的特征值分析方法。由于全周計算域具有循環對稱這一空間特性，基于全周計算域的特征值分析可以在考慮鄰近通道的影響作用下僅根據基于單通道計算域進行的特征值分析結果計算得到。然而，文獻［21］中使用了稀疏矩陣直接求解器進行大型稀疏陣的直接求逆［22］，計算的時間成本和內存消耗極高。如利物浦大學的Timme 所指出的，其極高的內存消耗可能成為研究計算量巨大的真實三維壓氣機流動失穩問題的瓶頸［12］。

為了克服試驗測量和非定常數值模擬各自的缺點，受到北京航空航天大學孫曉峰團隊，帝國理工Schmid 團隊和利物浦大學Timme 團隊等相關研究工作［12，17，21，23］的啟發，本文提出了一種基于三維RANS 方程的能夠精確解析真實壓氣機中三維幾何和流動細節的高效特征值分析方法。由于這是對壓氣機流動進行完全基于三維RANS 方程的特征值分析的首次嘗試，不失一般性，本文僅將該方法應用于某跨聲速壓氣機環形葉柵算例中。實際上，該方法可以直接擴展到全三維壓氣機流動失穩問題的研究中。本文的總體結構如下：第1節對特征值分析方法的具體算法進行了介紹，包括定常流計算以及特征值分析算法2 部分。第2節對定常流動計算和特征值分析結果進行了討論。最后，本研究的結論在第3 節中進行總結。

1 數值方法

1.1 流動求解器

本研究中所開展的定常、非定常模擬和特征值分析研究皆基于自研有限體積RANS 的CFD求解器NutsCFD。NutsCFD 能夠在旋轉坐標系下任意非結構化網格上實現定常、非定常流場求解。其算法和實現細節已在文獻［24］中詳細說明，此處不再詳細論述，僅作簡要介紹。

在以恒定角速度ω轉動的固定在轉子上的相對坐標系下，積分形式的流動控制方程為

式中：W為守恒變量；向量和分別為固定于轉子參考系中的對流和黏性通量；Fω為由于旋轉而產生的附加源項。它們的定義分別為

式中：ρ為密度；p為壓力；T為溫度；u為絕對速度；E為總能量；H為總焓；n為通量面單位法向量；τ為應力張量；κ為導熱系數。除了平均流場的控制方程外，還包括由Spalart-Allmaras 湍流模型產生的一個額外方程［25-26］。

平均流場和湍流方程的對流通量均采用Roe格式進行離散。原始變量的梯度由格林高斯公式計算獲得，用于從節點到通量面進行流動變量的外插值，從而使平均流場具有空間二階精度。湍流方程的對流通量使用一階迎風格式進行離散，以獲得更好的收斂魯棒性。時間離散基于全局化牛頓方法，采用偽瞬態延拓技術獲得更為穩定的殘差收斂性質。將所有空間離散通量和源項集中在一個向量R中，解向量W（包括平均流場和湍流變量，每個網格點6 個變量）基于全局化牛頓方法進行如式（3）所示的全隱式時間推進，直到找到穩態解。

在定常流場求解器的基礎上，對時間離散項稍作修改即可得到利用二階向后差分（Backward Difference Formula 2， BDF2）實現的雙時間步時域非定常求解程序。對于一個給定的物理時間步長，內迭代僅需求解附加了時間源項的定常流場問題，同樣可通過由ILU（0）預處理的GMRES 進行求解。

1.2 特征值求解器

空間離散后，半離散形式的流動方程可以表達為

式（4）可看作一個高維動力學系統。需要注意的是，為了與用于研究動力系統的符號使用習慣保持一致，右端殘差項的符號與流動控制方程中的殘差符號相反。此外，有限體積法產生的體積加權以及邊界條件都已包含在殘差函數R中。因此，右端項對于自變量W線化所產生的雅可比矩陣可準確描述非定常小擾動流場的時空變化信息。式（4）的特征值分解為

式中：J為由非線性殘差精確線化得到的雅可比矩陣；v為對應于特征值λ的右特征向量。穩態流場計算完全收斂后，在最后一次非線性迭代期間計算的流場雅可比矩陣按相應節點的控制體積將每一行進行縮放并乘以負號后即可得到J。特征值分析是基于系統矩陣J進行的。

為計算矩陣J的特征值和特征向量，本研究中使用了隱式重啟Arnoldi 方法（Implicitly Restarted Arnoldi Method， IRAM）［30］。所使用的廣義特征值問題數學庫為ARPACK［31］。為了更高效地使用內存并實現高效并行，ARPACK 與流場求解器通過子程序實現耦合，直接通過內存實現數據傳遞。與其他同樣調用ARPACK 數學庫的通用科學計算軟件（如MATLAB［32］或SciPy［33］）相比，將ARPACK 程序包與流場求解器緊密耦合的優勢在于更有效的內存使用和更容易的并行算法實現，這2 點優勢對于大規模流動特征值計算都至關重要。IRAM 通過構建Krylov 子空間{b，Jb，J2b，…，Jm b}，并在子空間中找到最優近似解作為特征向量。如果需要內特征值（即離原點最近的特征值）和特征向量，則將Arnoldi 過程應用于系統矩陣的逆，即構建Krylov 子空間{b，J-1b，J-2b，…，J-m b}。此時，矩陣J的內特征值成為矩陣J-1的外特征值，而IRAM 應用于外特征值問題時更易收斂。由于高雷諾數流動問題空間離散產生的雅可比矩陣直接求逆的成本極高［22］，在IRAM 方法中，并不直接計算大型稀疏陣J的逆陣。事實上，每個Arnoldi 步驟本質上需要的僅僅是矩陣向量積J-1x，而非J-1。因此，可以使用非線性流場求解器中已有的大型稀疏線性系統求解器來執行每一步Arnoldi 過程，無需對矩陣J求逆。在這項工作中，流場求解器中的ILU（0）預處理的GMRES（程序需稍作改動以求解復值線性方程組）用于求解大規模稀疏線性方程組Jy=x以獲得J-1x，用于在ARPACK 中構建IRAM 方法所需要的Krylov 子空間。當使用IRAM 方法進行特征值分析時，另一個需要提供的輸入量是種子向量b。本研究中，種子向量b設為隨機向量，原因是希望可以盡量保證所有特征向量在種子向量內都有相對平均的權重，從而不遺漏個別潛在的重要特征向量。此外，由于在特征值分析中并不計算所有特征值和特征向量，而是重點關注靠近虛軸的特征值，即接近失穩邊界的模態。然而這些特征值不一定是最靠近原點的，因此有必要在原矩陣基礎上施加一個復值的偏移量ζ，使得需要被求解的特征值在偏移后成為矩陣的內特征值從而被快速計算。實際操作中，首先根據工程經驗預估關鍵特征值所在范圍，選擇某個相近的復數偏移量ζ并將其施加于原始系統矩陣J。矩陣J中最接近于ζ的特征值即轉化為(J-ζI)的內特征值，其中I為單位矩陣。而矩陣(J-ζI)的內特征值可通過建立Krylov 子空間得到。同時，對矩陣施加一個復值偏移僅導致原始矩陣J的特征值偏移ζ，并不改變其特征向量。一旦獲得了相關的特征值和特征向量，就可以在背景流場附近重構非定常流場。如果相關的特征值以共軛復數成對出現，即λ1，λ2，…，λn和，，…，，則重構的非定常小擾動場可表示為

式中：復系數ci和由擾動初始值W(t)|t=0確定。

2 計算結果

2.1 算例介紹

本文研究的壓氣機算例是通過將NASA Rotor 67［34］的三維計算域與一圓柱面相交而生成的環形壓氣機葉柵。所選用圓柱面半徑為0.18 m，接近壓氣機50%葉高，此葉高處環形葉柵內流動為跨聲速。三維轉子和環形葉柵的幾何及z-rθ坐標系下的展開視圖如圖1 所示。表1 列出了環形壓氣機葉柵的一些關鍵參數。此環形壓氣機葉柵算例流場中來流相對馬赫數范圍約為1.05～1.2，代表了現代航空發動機跨聲速風扇葉尖附近截面的典型流場，因此本文中針對此算例開展研究。

表1 環形壓氣機葉柵關鍵參數Table 1 Key parameters of annular compressor cascade

圖1 NASA Rotor 67 與環形葉柵的幾何示意圖和流道內近峰值效率點壓力云圖Fig.1 Geometric diagram of NASA Rotor 67 and synthetic annular compressor cascade and pressure contours of flow path at near peak efficiency condition

2.2 穩態流場分析

為了計算穩態流場，使用NUMECA Autogrid 軟件（版本13.2）生成了3 套逐漸加密的全周22 通道計算網格（如圖2 所示），分別包含27 萬、49 萬和110 萬個網格點，所有網格近壁單元厚度為3×10-6m，滿足y+＜1。計算域進口邊界位于葉片前緣上游1.5 倍弦長處，進口條件為總壓101 325 Pa，總溫288.15 K，且規定進氣方向為軸向。計算域出口位于葉片尾緣下游2 倍弦長處，出口邊界使用常數背壓。所有計算中假定流動為完全湍流。為了獲得壓氣機特性曲線，通過增加背壓逐個計算穩態流場。當背壓增量為50 Pa 就導致流場求解器不收斂時，即認為找到了失速邊界。

圖2 3 套計算網格Fig.2 3 sets of computational meshes

使用3 套網格計算的特性曲線如圖3 所示。可見中網格與密網格特性計算誤差不大且失速邊界也很接近，因此可認為中網格已實現網格收斂。由于近失速點的流動特性是本文的主要關注點，因此詳細檢查了圖3 中用實心紅色符號標記的3 套網格上流量接近的工況（A， D， E）的流場細節。圖4 比較了該流量工況不同網格上的壓力云圖，未觀察到顯著差異。圖5 所示為沿葉片表面的壓力分布對比，其中ptin為進口總壓。可以看出，中、細網格的壓力分布，包括前緣附近的壓力尖峰，非常一致。而粗網格的結果與中、細網格結果有較大差異。這進一步表明，中網格足以捕捉足夠的流動物理細節，可用于后續特征值分析，因此本文的分析皆基于粗網格和中網格開展。

圖3 3 套網格上計算的壓氣機特性線（帶有實心紅色標記的點A、D、E 具有幾乎相同的質量流量，標有A、B、C 和D 的4 個數據點被用來進行特征值分析）Fig.3 Compressor characteristic curves computed on three meshes （Points A， D， E with solid red markers have nearly identical mass flow rates，and conditions labeled with A， B， C， and D are those for which eigenvalue analyses are performed）

圖4 工況A、D、E 的壓力云圖Fig.4 Pressure contours for flow solutions under Conditions A， D， E

圖5 不同網格下沿弦長的表面壓力系數分布Fig.5 Surface pressure coefficient distributions along chord for different meshes

2.3 特征值分析

針對粗網格上3 個近失速工況（A、B、C）以及中網格上質量流量與粗網格失速邊界接近的一個工況（D）進行特征值分析。為了尋找與旋轉失速現象相關的特征值和特征向量，需依賴一定的工程經驗，即模態波應具有與軸頻率接近的特征頻率。由于軸頻率為

若將雅可比矩陣根據該角頻率進行縮放，則所有與旋轉失速模態相關的特征值都被縮放為模長量級為1 的復數。

首先對粗網格上近失速工況（性能曲線上標記為A）的流場進行特征值計算。由于沒有關于特征譜的先驗信息，因此首次特征值搜索具有一定的隨意性，須通過試錯來獲得關鍵的特征值，即那些接近穩定邊界（虛軸）的特征值。具體來說，將復值偏移量ζ設為(0.5n，0.5m)，其中整數n、m取值為n=-6，-5，…，2 和m=0，1，…，6。對每個偏移值計算10 個特征值/特征向量。全部特征值計算完成后，進行檢查以確保找到復平面中矩形域( -3，1)×(0，3)內的所有特征值。所計算得到的特征值如圖6 所示。可以看出，對于近失速工況A，存在6 個不穩定模態，這意味著通過強隱式求解器得到的全環穩態流場在物理上是不穩定的［35］。隨后，按照與工況A 特征值分析相同的步驟，在流量較大的工況B 和工況C 進行特征值分析。工況B 和工況C 的特征值也畫在圖6中以進行比較。可以看出，隨著工況從失速邊界向更高的流量系數移動，特征值整體逐漸向左移動，意味著流動趨于穩定。表2 比較了3 個工況下最不穩定特征值的實部值?？墒褂梅侄尉€性插值來確定物理失速點流量為0.225 02 kg/s。根據相同的計算流程，還分析了中網格上工況D 的特征值，如圖7 所示?？梢钥闯?，盡管A、D 這2 種工況流量幾乎相同，但是中網格上工況D 不存在不穩定的特征值，比粗網格上工況A 更穩定。這與特性線計算時中網格上流場計算可以在更低的質量流量收斂到穩態解一致。

表2 工況A、B、C 下的最不穩定特征值Table 2 Least stable eigenvalues for Conditions A， B，and C

圖6 工況A、B 和C 的特征譜（在所有特征值中，工況B的11 個特征值標記為λ1～λ11）Fig.6 Eigenvalue spectra under Conditions A， B， and C （among all eigenvalues shown， eleven for Condition B are labeled with λ1-λ11）

圖7 工況A 與工況D 的特征譜Fig.7 Eigenvalue spectra under Conditions A and D

為了分析計算得到的特征值和特征模態，首先詳細檢查圖6 中標記的11 個特征值及其對應模態。所選特征值分別記為λ1，λ2，…，λ11。編號的理由將在后續段落中解釋。為了便于討論，首先基于模態2～模態6（即5 個最不穩定的模態）重構時間相關的擾動場。具體來說，對于任一特征值λ及其特征向量v，可構造以下非定常流動擾動場為

式中：角頻率ω為特征值λ的虛部。由于特征值實部僅決定模態振幅增長率，因此在重構中不予考慮。

基于式（8），可重構非定常流場擾動場?(t)在一個周期的任意時刻的取值。在此工作中，僅在t=0，T4，T2，3T4 這4 個時刻重構流場，根據奈奎斯特采樣定理，足以解析非定常擾動場的時空演化。其中T=2πω是對應模態的周期。需要注意的是，T對于每個模態都是不同的。

圖8 展示了使用第2～第6 個模態重構的一個周期內的流場擾動場。首先可以看出，失穩模態在激波與激波后邊界層及尾跡處幅值遠高于流場中其他位置，體現出與激波邊界層干涉極強的關聯性。同時還可以發現，所有重構的非定常擾動場都是從壓氣機葉片的壓力面傳播到吸力面的行波，波數分別為2、3、4、5 和6，可以認為此即模態型失速試驗中觀測到的模態波。換言之，根據Schmid 等的理論推導［21］，它們是循環對稱系統中某主模態在不同節徑下的表示。盡管這里只顯示了11 個被標記的模態中的5 個，但其余模態都遵循相同的規律，即具有不同節徑，而它們的下標即每個模態各自的節徑。表3 列出了在粗網格工況B 下不同模態波的節徑和歸一化轉速。如圖9 所示，與其他研究人員進行的典型計算和試驗研究［36］中所發現的模態波轉速與節徑的關系相比，盡管壓氣機幾何、運行條件和流動特性存在巨大差異，但仍能觀察到模態波轉速與模態節徑數（即周向波數）的高度正相關性。這是首次使用數值方法系統地揭示了這種相關性，其與其他研究者所發現規律的高度一致性表明本文所提出的高保真特征值分析方法確實捕捉到了模態擾動的固有特性。需要強調的是，非定常模擬很難明確地獲得多個失穩模態和接近失穩的穩定模態，特征值方法在揭示高維動力學系統的內蘊結構方面具有巨大的優勢。

表3 工況B 下粗網格上計算的不同模態波的節徑和基于轉軸速度歸一化后的轉速Table 3 Nodal diameter and normalized rotating speed of modal waves which is normalized by shaft speed on coarse mesh under Condition B

圖8 使用特征模態2～6 在一個周期內重構的軸向動量擾動場的時間演化Fig.8 Time evolution of axial momentum perturbation field reconstructed using Eigenmodes 2-6 over one period

圖9 工況B 下絕對參考系中的歸一化轉速與模態波的節徑之間的相關性及其與計算［7］和試驗［36］結果的比較Fig.9 Correlation between normalized rotating speed in absolute frame of reference and nodal diameter of modal waves for Condition B compared to computational［7］ and experimental［36］ data

2.4 時域非定常模擬

為驗證本文所提出的高效特征值分析方法的計算高效性和結果準確性，進一步針對粗網格近失速工況B進行了非定常模擬。非定常模擬由此工況完全收斂的定常流場進行初始化，通過雙時間步方法進行時間推進，每一個物理時間步內的迭代采用Newton-Krylov 方法實現殘差收斂。通過步長無關性研究后確定了物理時間步為0.36 μs，對應特征模態v3一個時間周期的1/200。非定常模擬總物理時間為3.3 s。在非定常模擬的基礎上，監測圖10 中葉片前緣附近綠色圓點所在位置靜壓值隨時間的演化，壓力擾動監測結果如圖11 所示。首先可以發現，在最初的1.5 s 物理時間內，壓力脈動值始終在10-10Pa 附近隨機振蕩。可以想象，若未先進行特征值穩定性分析，只通過非定常模擬進行研究，這一工況很可能被誤判為穩定工況，這一現象體現了特征值方法進行穩定性分析的優勢，即可以明確地判定流場的穩定性。從1.5 s 開始，監測點靜壓值開始指數增長。在約3.3 s 時，非定常流場中由于在計算域進口處出現回流，與所施加的總壓、總溫、軸向進氣等進口邊界條件沖突，導致計算迭代發散。由于本研究主要關注流場小擾動分析，而非充分發展的旋轉失速現象，因此本工作中未再作調整計算域或邊界條件以獲得穩定的旋轉失速現象的嘗試。

圖10 數值探針監控點位置示意圖Fig.10 Position schematic diagram of numerical sensor

圖11 數值探針位置處靜壓擾動變化曲線Fig.11 History of pressure perturbation at numerical probe position

為驗證計算所得特征值和特征向量的準確性，將非定常流場投影到特征值分析獲得的特征向量上進行非定常流場重構。為減小計算量，僅考慮圖6 中的11 個特征模態及其共軛模態，即若令

則重構非定常流場可以表示為

為了獲得η(t)，需要將式（10）兩邊分別左乘VH，即

系數向量η(t)的計算式為

式中：VHV為22×22 的方陣，只需對其求逆一次。對于VH?(t)項，在每一物理時刻，將非定常擾動向模態矩陣V投影即可得到，計算量相比于流場計算也可忽略不計。在獲得所有物理時刻的系數向量η(t)之后，即可通過獲得基于最不穩定模態的重構流場。

圖12 中顯示了通過非定常模擬得到的流場(t)和重構流場中在監測點的靜壓擾動值各自隨時間的演化歷史。可見其在3.3 s 前吻合很好。3.3 s 后由于非定常流場本身出現非線性，重構信號逐漸偏離時域模擬信號。通過此比對可以確認特征分析的準確性。

圖12 數值探針位置處非定常壓力擾動與采用模態v3 與 重構得到的壓力擾動的比較Fig.12 Comparison between unsteady pressure perturbation and reconstructed pressure perturbation by v3 and at numerical probe location

2.5 計算成本比較

本文工作中定常計算、特征值計算、非定常模擬都在24 核工作站上開展，所使用的中央處理器（Central Processing Unit， CPU）型號為Intel Xeon E52678W v4。為了定量分析計算效率，表4 列出了各類分析的實際計算時間。由于定常計算中各個工況的收斂速度依賴于初始化流場，且近堵點、最高效率點、近失速點的流場特性和收斂難度皆有所不同，因此對于定常分析單獨列出了工況B 下定常流場計算和整個特性線計算各自的計算成本。對于特征值計算，由于對于一個新的算例，在初期需要一定的試錯，以確保關鍵特征值都被找到，具體的試錯成本與算例高度相關。因此特征值的計算成本僅考慮了計算最不穩定的10 個特征值和特征向量的計算成本。

表4 粗網格上3 類分析的計算成本Table 4 Computational cost of three types of analysis on coarse grid

從表4 中計算成本的分析可以看出，特征值分析與非定常模擬相比加速約155 倍，卻獲得了遠比時域非定常模擬更為豐富的非定常小擾動流場信息。同時，特征值分析計算成本為單次定常流場計算的3.8 倍，且比特性線計算的成本要低很多，大約為計算整條特性線所需時間的28%。

3 結 論

1） 本文提出了一種適用于壓氣機三維流場流動穩定性研究的高效特征值分析方法，并用于分析代表了典型現代跨聲速壓氣機葉片截面的環形壓氣機葉柵內流場的特征值和特征向量。將隱式重啟Arnoldi 方法（IRAM）應用于內特征值分析，以識別與旋轉失速現象關聯最緊密的特征模態。內特征值分析使用了偏移求逆技巧。由于大型稀疏陣求解成本極大，因此利用大型稀疏線性方程組求解器，即ILU（0）預處理的GMRES 求解器的復值版本，等效地實現了求逆的效果。流場求解器中所采用的Newton-Krylov方法本就會計算雅可比矩陣，且此非線性流場求解方法已得到廣泛驗證，其計算效率可與最先進的CFD 求解器相媲美，經過簡單修改就能夠處理復值線性系統求解問題，從而進行大規模流動問題的特征值分析。與此同時，將較為成熟的特征值分析程序包ARPACK 與非線性流場求解器中的雅可比矩陣計算和線性求解器等現有程序相耦合，可實現工具庫間的緊密數據傳遞，從而實現更高效的內存管理和更高的并行效率。本文所發展的高效特征值計算方法為基于三維流動控制方程的特征值分析在壓氣機流動穩定性問題上的首次嘗試，由于其計算效率遠高于時域非定常模擬計算方法（對于本文所研究算例，小擾動穩定性分析比非定常模擬加速約155 倍，僅為特性線計算的28%），使得基于小擾動分析的穩定性研究在考慮三維效應的壓氣機流動穩定性上成為可能。

2） 對于所研究的壓氣機算例，首先在失速點附近計算了若干深度收斂的穩態流場，然后使用上述特征值分析方法對部分近失速工況流場進行了特征值分析。由于特征值分析方法完全基于三維流動方程的精確雅可比矩陣，因此可獲得極其豐富且與非定常小擾動流場分析完全一致的信息。同時，由于沒有對壓氣機及其流動進行任何幾何和流場簡化，因此特征值和特征向量與非線性流動求解器所解析的流動物理完全一致。特征值分析結果表明，此壓氣機算例不穩定模態波對于流場的擾動主要體現在激波與激波后邊界層，意味著旋轉失速的發生與激波邊界層干擾有關，為解釋旋轉失速失穩機理提供了一種重要思路。同時，所有特征值方法可以一次獲得大量失穩或接近失穩的特征模態，發現在失速邊界附近存在一系列具有連續變化的節徑的特征模態。不同模態根據各自的節徑和轉速以行波形式朝葉片旋轉相反方向傳播。對于略超出失速邊界的工況，數個特征模態會同時失穩，預示了失穩過程和失速后流場的高度復雜性。此外，本研究還探究了模態波轉速與模態節徑之間的關系，解釋了試驗中大量出現的模態波節徑越大絕對坐標系下轉速越高這一現象。通過本文研究可以看出，由于所發展穩定性分析方法的準確性和高效性，相比于非定常計算具有巨大的優勢，為針對真實三維壓氣機模型的流動穩定性分析和失穩機理研究提供了重要的工具。