壓氣機穩定邊界的數值預測研究

靳 楠，胡 駿，李 駿

（南京航空航天大學能源與動力學院，南京210016）

壓氣機穩定邊界的數值預測研究

靳 楠，胡 駿，李 駿

（南京航空航天大學能源與動力學院，南京210016）

為研究不同壓氣機穩定性預測方法的準確性，采用最大靜壓升系數法與商用CFD軟件，對某4級低速壓氣機的穩定邊界進行了數值模擬，并與試驗結果對比。計算結果表明：最大靜壓升系數法在預測壓氣機穩定邊界上有一定的準確性。使用商用CFD軟件NUMECA對某4級壓氣機進行模擬計算，考慮了不同的網格劃分方案、使用不同的湍流模型和差分格式等因素，以出現數值發散作為壓氣機失穩的標準，結果表明：主流區的網格對計算結果影響較大。不同差分格式的結果顯示，中心差分格式的結果最好。

軸流壓氣機；穩定邊界；數值模擬；最大靜壓升；航空發動機

0 引言

航空燃氣渦輪發動機的風扇/壓氣機的長度、質量在整臺發動機中占比很大，生產、維護費用高昂，其性能不僅影響推重比、耗油率等性能參數，而且對整機穩定性也影響巨大，其性能優劣直接影響整機的研制。發動機的設計總是要求風扇/壓氣機在保持高效率和足夠穩定裕度的條件下，實現以更少的級數達到更高的增壓比，同時還要滿足強度、可靠性、噪聲、隱身等方面的要求[1]。為了增加發動機推重比、降低耗油率，風扇/壓氣機向高壓比和高級負荷方向發展，必然使風扇/壓氣機端壁和葉片吸力面的附面層承受更大的逆壓力梯度，在這種情況下風扇/壓氣機內部流動更容易分離，其穩定性就難以得到保證，因此對氣動穩定性研究越發重要。航空發動機推重比的不斷提高，需要高性能、高效率和高穩定性的風扇/壓氣機，也需要不斷發展完善風扇/壓氣機氣動設計體系中與之相適應的設計和計算方法。

自航空燃氣渦輪發動機出現以來，壓氣機穩定性一直是壓氣機設計和研究人員關注的問題。隨著研究人員對壓氣機內部流動認識的不斷深入，主要有3類針對壓氣機穩定性的分析方法。第1類是經驗關聯法。在研究壓氣機穩定性的初期階段，研究人員通過對大量的試驗數據進行分析和理論的研究，建立起判斷壓氣機失穩的準則參數，并將準則參數與壓氣機的幾何參數和氣動參數等關聯起來，把準則參數作為壓氣機失穩的判斷標準。在這階段，Emmons[2]、Stenning[3]、Marble[4]、Sears[5]、Krieble[6]、Dunham、Lieblein[7]、Koch[8]等人隨著研究的不斷深入，提出了多種經驗準則來評估壓氣機的穩定邊界，其中代表性的有：Dunham提出了將特性線上壓升系數對流量系數的導數為零點作為判定壓氣機失速的準則；Lieblein[8]提出了將擴散因子、擴散率作為葉柵失速發作的近似準則；Koch提出用壓氣機最大靜壓升系數來評估壓氣機的壓升潛力[9]。這類型的方法依據試驗數據，使用經驗關系式來判斷壓氣機狀態點是否失穩，因此有一定的誤差，但是其優點是使用簡單、便捷。第2類方法以小擾動的線性穩定性理論為主，用于預測風扇/壓氣機或壓縮系統的工作穩定性。該類方法已發展出考慮多葉片排的不可壓縮和可壓縮穩定性模型。自 20世紀 50年代以來，Emmons、Stenning[3]、Marble[9]、Sears、Fabri和 Moore 等人的工作都屬于此類，用于穩定性分析的渦方法也大多屬于此類方法。國內外學者對這一類方法進行了較多的研究與發展，具體成果見參考文獻[10-11]；第3類方法基于壓縮系統的動力學模型，不僅關心風扇/壓氣機或壓縮系統在平衡工作點受小擾動的穩定性，還非常重視有限擾動對穩定性的影響及其發展情況，以及風扇/壓氣機失速后類失速流態—旋轉失速和喘振的性質，這類方法為建立在時域上的時間推進法，以Moore和Greitzer的研究工作為代表[12-14],近年來此類方法得到了較大的發展，例如參考文獻[15-16]。

1 方法介紹

C.C.Koch通過對大量的低速軸流式壓氣機試驗結果和有限的高速風扇/壓氣機試驗數據的分析，將風扇/壓氣機葉柵通道類比于2維擴壓器，提出了最大靜壓升系數的概念來評估壓氣機的壓升潛力[9]。將壓氣機的最大靜壓升能力與葉柵無因次長度比L/g2、工作雷諾數、無因次軸向間隙比△z/s和無因次徑向間隙比ε/g關聯起來，分別如圖1、3、4、5所示。其中

t為柵距，而其他幾何參數的定義如圖2所示。

定義風扇/壓氣機運行過程中的氣動有效靜壓升系數（Ch）ef，用來描述壓氣機某一工作狀態的壓升。當有效靜壓升超過壓氣機的壓升潛力，即表示壓氣機失穩。

式中：V、W和U分別為絕對速度、相對速度和牽連速度；V1，ef為葉柵進口的有效動壓頭速度，定義為

其中：

當（α1+β1）≤90°時，Vmin=V1sin（α1+β1）

當（α1+β1）＞90°時，Vmin=V1

當 α1 ＜0° 時 ，Vmin=U（4）

當 β1＜0°時，Vmin=U

各速度分量的定義如圖2所示。

圖1中的結果是在Re=130000、無因次徑向間隙ε/g=0.055、以及無因次軸向間隙△z/s=0.38的情況下獲得的。以此為參考點，當這些參數變化時，需要對圖1的結果進行修正，引入工作雷諾數、無因此軸向間隙、無因此軸向間隙的影響

（Ch）s=KRe×Kε/g×K△z/s×（Ch）s，ref（5）

式中：雷諾數修正系數KRe、葉尖間隙修正系數Kε/g和葉排軸向間隙修正系數K△z/s的取值分別由圖3、4、5中的曲線確定。

2 基元葉柵法與判穩方法結合計算壓氣機穩定邊界

壓氣機氣動穩定性的基本技術問題就是在均勻進氣的條件下穩定邊界的預測以及穩定邊界與其設計參數相關性的問題。穩定邊界預測一般與非設計點性能預測結合使用。基元葉柵法就是利用已獲得的關于基元葉柵的落后角和總壓損失與其主要影響因素（如馬赫數Ma，雷諾數Re和攻角i等）的關聯式，結合無葉區流場的計算預測壓氣機非設計點的氣動性能。有關基元葉柵法的詳細介紹請參見文獻[17]。通過基元葉柵法得到壓氣機流場信息，如速度、壓力、溫度等參數，根據流場信息，依靠最大靜壓升系數法[9]即可判斷壓氣機是否穩定。

本文將基元葉柵法與最大靜壓升系數法相結合，發展了1種預測壓氣機穩定邊界的數值方法，并以某4級低速壓氣機為算例，驗證了該方法的有效性。

所用的壓氣機為4級低速模擬壓氣機，采用類重復級設計，第3級為模擬級，前2級和第4級分別為模擬級提供恰當的進、出口條件，其設計轉速為900 r/min，外徑為1.5 m，輪轂比為0.88，轉子葉尖速度為70.7 m/s，轉、靜子展弦比分別為0.957與1.5，轉子與靜子中徑處葉柵稠度分別為1.528、1.625。

試驗得到的壓氣機特性線與計算結果如圖6所示。分別對低速壓氣機在設計轉速為900、700、500 r/min下的試驗結果與計算結果做了對比如圖7所示。

不同轉速下邊界點流量和壓比的相對誤差如圖8所示。從圖中可見，基元葉柵法與最大靜壓升系數法結合計算所得的穩定邊界與試驗結果比較接近，邊界點相對誤差如圖8所示。

使用最大壓升系數法判斷壓氣機穩定邊界，得到壓氣機在不同轉速 (500、700、900 r/min)下，靜壓升系數的 4 級分布，分別如圖 9（a）、（b）、（c）所示。

從圖中可見，在不同轉速下，靜壓升系數的4級分布趨勢非常一致。4級的失速靜壓升系數幾乎相同，由此表明4級的壓升潛力差別不大，這是因為所用的4級壓氣機為類重復級壓氣機，各級葉片幾何結構近似相同導致的。在不同流量點下，有效靜壓升系數逐級增大，在失速點均為第4級時有效靜壓升系數超過了最大靜壓升系數，計算結果表明壓氣機的第4級首先出現失穩的現象。在失速點處各級的靜壓升系數差別不大，且比較接近最大靜壓升系數。

3 4級低速壓氣機數值模擬研究

為了校驗數值模擬結果的可靠性，同時找到最優的數值模擬方法，使用商用軟件NUMECA對4級壓氣機進行數值計算，以出現數值發散作為壓氣機失穩的標準，得到壓氣機穩定邊界。分別改變網格劃分方式、湍流模型以及差分格式進行計算，將計算結果與試驗結果進行比較，比較不同的網格、湍流模型、差分格式對CFD計算的穩定邊界的影響，并與實驗結果做出對比。

3.1 不同網格劃分方式數值模擬的結果分析

合適的網格生成技術往往是正確模擬一些復雜流動的關鍵所在，通過改變網格的分布情況及控制總網格量，主要從主流區均勻網格百分比、葉尖間隙內網格節點數以及總網格量（改變葉片區網格量）3方面考慮，得到6種不同計算網格（見表1），分別進行計算后尋找理想的網格劃分方式。6種不同網格劃分方式中Grid1～Grid4的總網格量均接近500萬，故近似認為他們總網格量相同；由于模型壓氣機的設計采用類重復級設計，故所有靜子葉片（包括進口導葉）網格量相同，所有轉子葉片網格量也相同。除網格不同外其余計算條件均相同。

表1 網格設計方案

在設計轉速下壓氣機特性如圖10所示。從圖10（a）中可見，大流量點的效率均偏離試驗數據較大，Grid2整體比試驗數據偏低。從圖10（b）中可見，除Grid2外，其他算例的壓比特性均比試驗值高，Grid1與Grid6的壓比特性幾乎重合。考慮CFD方法計算的壓氣機穩定邊界點，Grid5和Grid2的邊界與試驗結果最為接近。

Grid1、Grid2、Grid3、Grid4 以及試驗數據在邊界點的壓比特性比較如圖11所示。從圖中可見，4個算例總網格量均相同，與Grid1相比，Grid2主流區網格占比減少，Grid3葉尖間隙網格數減少，Grid4葉尖間隙網格數增加。主流區網格占比對計算的而穩定邊界影響較大。

Grid5、Grid1、Grid6 的網格量分別為 300 w、500 w、720 w，葉尖間隙網格量、主流區網格量占比均相同。邊界點流量和壓比相對誤差如圖12所示。從圖中可見，Grid5計算的穩定邊界流量與試驗值最為接近，但是邊界點壓比與試驗值相差最大。而Grid1與Grid6穩定邊界流量與試驗值差別較大，但是壓比與試驗值較為接近。500 w網格與720 w網格的計算結果相差不大。

3.2 不同湍流模型數值模擬結果分析

在利用原始變量求解雷諾平均N-S方程組時需用湍流模型對其進行封閉求解，故湍流模型的選擇是影響壓氣機內部流動數值模擬精度的重要因素。在選定Grid1網格劃分方式的基礎上改變計算中所使用的湍流模型，與試驗進行比較分析，得到較為理想的湍流模型。所采用的湍流模型有S-A（Spalart-Allmaras）模型（一方程模型，魯棒性較好，能較好地解決一些復雜流動）、Yang-Shih k-ε模型（兩方程模型，低雷諾數湍流模型，可以更加精確地求解壁面附近的流動）以及SST k-ω模型（兩方程模型，改進后的k-ω模型，在廣泛的流動領域中有更高的精度和可信度）。在計算中，除湍流模型不同外其余計算條件均相同。

不同湍流模型計算所得的壓比特性和效率特性如圖13（a、b）所示。從圖中可見，壓比特性線均高于試驗值，效率特性線在小流量點上與試驗值接近，在大流量點上與試驗值差別較大。就穩定邊界而言，使用k-ε模型所得到的穩定邊界最接近試驗所得邊界。

使用不同湍流模型所得的邊界點流量均與試驗值有一定差別，相對誤差如圖14。

3.3 不同差分格式數值模擬結果分析

在數值模擬中選擇合理的計算格式進行計算可以提高數值模擬的精度，進而解決一些較為復雜的流動問題。本節是在前兩節的基礎上，選定Grid1的網格劃分方式，S-A湍流模型，對N-S方程中空間項的離散方法進行改變，將計算得到的結果與試驗結果進行比較分析，得到較理想的計算格式。本文中差分格式分別采用中心差分格式及迎風格式（TVD格式和通量差分格式，均為2階精度）。在相同計算條件下，中心差分格式計算得到的結果誤差比迎風格式更小，但當流動中出現如激波這樣的間斷時，中心差分格式的數值解會產生波動，得到的結果誤差會較大；TVD格式是建立在非線性限制器的概念上的，力圖阻止數值波動的產生，使得數值解的變化以非線性的方式得以控制；通量差分格式是在每個網格步求解時對相鄰間斷的控制方程進行描述。本節的計算，除空間差分格式不同之外其余計算條件均相同。

改變差分格式計算所得的壓比和效率特性如圖15（a、b）所示，從圖中可見，壓比特性均高于試驗值，中心差分格式與TVD格式在小流量點的效率特性與試驗值較為符合，大流量點的效率特性與試驗值差別較大，通量差分格式效率特性與試驗值符合度較差。就穩定邊界而言，使用TVD格式計算出的邊界點與試驗值差別較大。總而言之，TVD格式對特性的模擬效果較差。

邊界點流量和壓比相對誤差如圖16（a、b）所示，從圖中可見，中心差分格式邊界點流量和壓比與試驗最接近。TVD格式邊界點流量與壓比偏離試驗值最遠。中心差分的計算結果最好。

4 總結與展望

針對4級低速模擬壓氣機的穩定邊界，以最大靜壓升系數法為主要方法，商用CFD軟件數值模擬為參考，進行了一系列的數值計算，并與試驗結果做出對比。得到的主要結果如下：

（1）使用最大靜壓升系數法判斷壓氣機的穩定邊界有一定的準確性。

（2）分析通過改變網格劃分方案的計算結果，主流區網格占比對壓氣機穩定邊界影響較大。分析不同差分格式的結果，中心差分格式的結果最好。使用NUMECA計算的壓氣機穩定邊界與試驗結果有一定的誤差。下一步的工作將會圍繞使用最大靜壓升系數法修正CFD計算結果來展開。

[1]劉永泉.國外戰斗機發動機的發展與研究[M].北京：航空工業出版社，2016，261-262.LIU Yongquan.Development and investigation of foreign fighter engine[M].Beijing:Aviation Industry Press,2016：261-262.(in Chinese)

[2]Emmons H W,Pearson C E,Grant H P C.Compressor surge and stall propagation [J].Transactions of the ASME,1955，77（4）:455-469.

[3]Stenning A H,Kriebel A R,Montgomery S R.Stall propagation in axial-flow compressors[J].Massachusetts Institute of Technology,1956,115:163-8.

[4]Marble F E.Propagation of stall in a compressor blade row[J].Journal of Aeronautical Sciences,1955,22（8）:541-544.

[5]Sears W R.Rotating stall in axial compressors [J].Zeitschrift für Angewandte Mathematik und Physik,Birkh?user Verlag AG Basel,Switzerland,1955,35（6）:429-455.

[6]Kriebel A R,Schuuind R G.Stall propagation in a cascade of airfoils[R].NACA-TRR-61,1960.

[7]Leiblein S,Roudenbush W H.Theoretical Loss relation for low-speed two dimensional cascade flow[R].NACA-TN-3662,1956.

[8]Koch C C.Stalling pressure rise capability of axial flow compressor stages[J].Journal of Engineering for Power,1981,103:645-656.

[9]Marble F E.Propagation of stall in a compressor blade row[J].Journal of Aeronautical Sciences,1955，22（8）：541-544.

[10]劉超.燃氣渦輪發動機整機穩定裕度數值模擬方法研究[D].南京航空航天大學,2012.LIU Chao.Aerodynamic stability numerical simulation method study of compressor component running in gas turbine engine[D].Nanjing University of Aeronautics and Astronautics,2012.(in Chinese)

[11]侯睿煒,孫曉峰.一種基于數值計算的軸流壓氣機旋轉失速穩定性模型——之一:理論模型及算法分析[J].航空動力學報,2009,24（12）:97-104.HOU Ruiwei,Sun Xiaofeng,Liner stability model of rotating stall based on direct numerical techniques:part I-theoretical model and algorithm[J].Journal of Aerospace Power,2009,24(12):97-104.(in Chinese)

[12]Moore F K.A theory of post-stall transients in axial compression systems:part I—development of equations[J].Journal of Engineering for Gas Turbines&Power,1986,108（1）:68-76.

[13]Greitzer E M,Moore F K.A theory of post-stall transients in axial compression systems:part II—application[J].Journal of Engineering for Gas Turbines&Power,1986,108（2）:231-239.

[14]Moore F K.Stall transients of axial compression systems with inlet distortion[J].Journal of Propulsion&Power,2015,2（6）:552-561.

[15]張環.旋轉總壓畸變對壓氣機穩定性影響的研究[D].南京航空航天大學,2009.ZHANG Huan,An investigation of effects of rotating total pressure distortion on aerodynamic stability of compressors[D].Nanjing U-niversity of Aeronautics and Astronautics,2009(in Chinese)

[16]王聰,文彬鶴,司文杰,等.軸流壓氣機旋轉失速建模與檢測:Ⅰ——基于確定學習理論與高階Moore-Greitzer模型的研究[J].自動化學報,2014,40（7）:1265-1277.WANG Cong,WEN Binhe,SI Wenjie,et al,Modeling and detection of rotating stall in axial flow compressors:part I-investigation on high-order M-G models via deterministic learning[J].Acta Automatica Sinica,2014,40（7）:1265-1277.(in Chinese)

[17]Johnsen I A,Bullock R O.Aerodynamic design of axial-flow compressors[R].NASA SP-36.1965.

Numerical Predicting Study on Compressor Surge Margin

JIN Nan,HU Jun,LI Jun
（College of Energy and Power Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,China）

In order to evaluate the accuracy of different method in predicting surge margin of axial compressor,a fourth-stage lowspeed model compressor was numerically simulated by maximum pressure rise method and commercial CFD software MUNECA,and the result was compared to experiment data.The calculation result shows that the maximum pressure rise method has certain accuracy according to experiment data.The simulation calculation of a fourth stage compressor was performed using NUMECA software,considering the influence of meshing means,different turbulence models and different scheme on surge margin,and taking numerical divergence as the instability standard of compressor.The results show that the factors in consideration all influence the accuracy of the simulation results,especially grid percentage of mid-flow cells,and the result of central difference scheme is the best in three schemes.

axial compressor;surge margin;numerical simulation;maximum pressure rise method;aeroengine

V231.3

A

10.13477/j.cnki.aeroengine.2017.04.009

2016-11-16

靳楠（1991），男，在讀碩士研究生，主要研究方向為內流氣動熱力學；Email:jinnan666@foxmail.com。

靳楠，胡駿，李駿.壓氣機穩定邊界的數值預測研究[J].航空發動機，2017，43（4）：48-55.JINNan，HUJun，LI Jun.Numerical predictingstudy on compressor surge margin[J].Aeroengine，2017，43（4）：48-55.

（編輯：張寶玲）