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

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

2．中國(guó)空氣動(dòng)力研究與發(fā)展中心，綿陽(yáng) 621000

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

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

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

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

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

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

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

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

1 數(shù)值方法

1.1 流動(dòng)求解器

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

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

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

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

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

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

1.2 特征值求解器

空間離散后，半離散形式的流動(dòng)方程可以表達(dá)為

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

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

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

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

2 計(jì)算結(jié)果

2.1 算例介紹

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

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

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

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

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

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

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

2.3 特征值分析

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

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

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

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

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

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

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

表3 工況B 下粗網(wǎng)格上計(jì)算的不同模態(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 在一個(gè)周期內(nèi)重構(gòu)的軸向動(dòng)量擾動(dòng)場(chǎng)的時(shí)間演化Fig.8 Time evolution of axial momentum perturbation field reconstructed using Eigenmodes 2-6 over one period

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

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

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

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

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

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

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

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

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

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

圖12 數(shù)值探針位置處非定常壓力擾動(dòng)與采用模態(tài)v3 與 重構(gòu)得到的壓力擾動(dòng)的比較Fig.12 Comparison between unsteady pressure perturbation and reconstructed pressure perturbation by v3 and at numerical probe location

2.5 計(jì)算成本比較

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

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

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

3 結(jié) 論

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

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