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

1．西北工業(yè)大學(xué) 動力與能源學(xué)院，西安 710072

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

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

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

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

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

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

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

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

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

1 數(shù)值方法

1.1 流動求解器

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

在以恒定角速度ω轉(zhuǎn)動的固定在轉(zhuǎn)子上的相對坐標(biāo)系下，積分形式的流動控制方程為

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

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

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

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

1.2 特征值求解器

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

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

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

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

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

2 計算結(jié)果

2.1 算例介紹

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

表1 環(huán)形壓氣機葉柵關(guān)鍵參數(shù)Table 1 Key parameters of annular compressor cascade

圖1 NASA Rotor 67 與環(huán)形葉柵的幾何示意圖和流道內(nèi)近峰值效率點壓力云圖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 穩(wěn)態(tài)流場分析

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

圖2 3 套計算網(wǎng)格Fig.2 3 sets of computational meshes

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

圖3 3 套網(wǎng)格上計算的壓氣機特性線（帶有實心紅色標(biāo)記的點A、D、E 具有幾乎相同的質(zhì)量流量，標(biāo)有A、B、C 和D 的4 個數(shù)據(jù)點被用來進行特征值分析）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 不同網(wǎng)格下沿弦長的表面壓力系數(shù)分布Fig.5 Surface pressure coefficient distributions along chord for different meshes

2.3 特征值分析

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

若將雅可比矩陣根據(jù)該角頻率進行縮放，則所有與旋轉(zhuǎn)失速模態(tài)相關(guān)的特征值都被縮放為模長量級為1 的復(fù)數(shù)。

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

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

圖6 工況A、B 和C 的特征譜（在所有特征值中，工況B的11 個特征值標(biāo)記為λ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

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

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

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

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

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

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

圖9 工況B 下絕對參考系中的歸一化轉(zhuǎn)速與模態(tài)波的節(jié)徑之間的相關(guān)性及其與計算［7］和試驗［36］結(jié)果的比較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 時域非定常模擬

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

圖10 數(shù)值探針監(jiān)控點位置示意圖Fig.10 Position schematic diagram of numerical sensor

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

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

則重構(gòu)非定常流場可以表示為

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

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

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

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

圖12 數(shù)值探針位置處非定常壓力擾動與采用模態(tài)v3 與 重構(gòu)得到的壓力擾動的比較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 下定常流場計算和整個特性線計算各自的計算成本。對于特征值計算，由于對于一個新的算例，在初期需要一定的試錯，以確保關(guān)鍵特征值都被找到，具體的試錯成本與算例高度相關(guān)。因此特征值的計算成本僅考慮了計算最不穩(wěn)定的10 個特征值和特征向量的計算成本。

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

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

3 結(jié) 論

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

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