多級軸流壓氣機二維性能預測方法

楊曉飛，蔣永松，潘若癡，孟德君

中國航發沈陽發動機研究所，沈陽 110015

壓氣機特性計算主要有一維計算、二維計算、三維數值模擬3種方法。一維計算是在壓氣機平均流面上進行，當壓氣機非設計點的流動參數徑向發生扭曲[1]，或流路收縮嚴重時，其計算精度降低。三維數值模擬在20世紀90年代商業化后迅速普及，并在壓氣機氣動性能分析中扮演重要角色，尤其是三維定常雷諾平均數值模擬已是葉輪機氣動領域不可或缺的計算方法。然而，湍流模型、摻混面模型等所進行的簡化會造成計算誤差，邊界條件的不同、幾何細節差異也會對計算結果產生影響[2-3]。隨著壓氣機級數增多、負荷升高，三維數值模擬的精度下降，性能分析完全依賴三維計算存在巨大的風險；除此之外，三維計算花費時間長、數據量大，在壓氣機方案設計階段完全采用三維計算使得設計周期延長。在三維數值模擬“浪潮”漸漸回歸平靜后，S2流面正問題計算仍是壓氣機特性評估及流場分析的重要手段。

吳仲華先生在1952年提出了兩類流面迭代計算葉輪機械內部流動的理論，基于此理論發展的流線曲率法或通流矩陣法仍是S2流面求解的主要方法[4-6]。兩種方法都采用定常、絕熱、無黏假設，黏性的影響通過引進經驗或半經驗的模型來彌補。這些模型必須經過大量試驗數據和設計經驗不斷改進、校核和驗證。隨著壓氣機設計技術的進步及對內部流動機理認識的深入，Koch和Smith[7]、Konig等[8-9]研究了葉型損失計算方法，Bloch[10]、Boyer和Brien[11]對雙激波模型計算方法做了重要論述。近幾年，Milan[12]、Kim[13]、Azamar Aguirre[14]等持續對經驗模型進行了研究；同時，為提升計算精度，適應不同理念葉型的損失、落后角預估，S1流面計算在模型發展中得到應用[15-17]。中國對于經驗模型的研究亦在不斷推進，劉波[18-19]、季路成[20]、滕金芳[21]等開展了有關工作。對于工程設計，中國的壓氣機技術體系是較為完備的，但更新完善卻進展緩慢[22]。隨著中國多級高壓壓氣機的發展，穩定、高精度的S2正問題計算程序需求日益迫切。

無論采用何種流線曲率方程，其求解過程都是相似的，計算的精度主要受經驗、半經驗模型的影響，這些模型也正是本研究的重點。本文研究面向工程需要，為了提高計算精度，滿足方案評估要求，在一個流線曲率計算程序的基礎上，通過對比分析，優化了攻角、落后角及損失計算方法，尤其改進了端壁損失計算方法；同時，為克服基于傳統平面葉柵試驗數據的攻角、落后角模型與新設計理念葉型之間的不匹配問題，采用S1流面計算修正S2的方法，最終發展出了一個高精度的S2正問題計算程序；利用3個不同負荷水平的、經試驗驗證的多級壓氣機進行了檢驗計算，證實了方法的計算精度和穩定性。目前，基于該方法的流線曲率程序已在工程設計及分析中開始應用。

1 流線曲率程序

利用流線曲率法求解完全徑向平衡方程是S2流面計算的主要方法之一，在工程中應用也最廣泛。一個常見的完全徑向平衡方程為

(1)

式中：ρ為密度；p為壓力；V、W分別為絕對和相對速度；r為半徑；rm為子午流線的曲率半徑；φ為流線與軸向的夾角；“m”表示子午方向，“θ”表示圓周方向。式(1)中忽略了葉片力。

本文研究基于的平臺為中國航發沈陽發動機研究所的流線曲率程序，該程序采用基于軸對稱假設的流線曲率法求解S2流場，具有包括靜子角度調節在內的全轉速特性計算功能。輸入參數主要包括流路、進出口葉型角、稠度、最大相對厚度及位置、最大撓度位置、葉尖間隙以及轉速等。因不具備摻混計算，且模型未曾改進，所以程序對當前多級壓氣機的計算精度較低。該程序較為完整的框架為本文的經驗、半經驗模型研究提供了良好基礎。

2 經驗模型

2.1 最小損失攻角、落后角模型

目前，最小損失攻角的主流計算方法為NASA SP-36[23]中的方法，本文仍沿用這一方法。本程序的最小損失攻角在來流馬赫數Ma≤1.0時，利用式(2)計算；在Ma≥1.2時，使用葉背相切法計算，見式(3)；在1.0

Ma≤1.0

(2)

i*=im2Ma≥1.2

(3)

i*=im2+(im1-im2)f(Ma) 1.0

(4)

落后角計算模型為

(5)

(6)

(7)

2.2 設計點損失模型

與攻角、落后角模型的匱乏相反，自流線曲率法在工程中應用以來，損失模型不斷被提出，葉型損失、激波損失、二次流損失、端壁損失計算模型相對豐富，但每種模型都有使用范圍，不同模型之間的良好匹配是二維正問題程序計算的基礎。本文發展的設計點損失計算方法為

ωt=ωp+ωs+ωsec+ωe

(8)

式中：ωt表示總損失；等號右側4項依次為葉型損失、激波損失、二次流損失和端壁損失。

2.2.1 葉型損失

葉型損失由葉片表面附面層和尾跡摻混形成。Lieblein和Roudebush根據葉柵試驗建立了計算葉型損失的基本方法，即根據擴散程度與尾跡動量厚度之間的關系計算葉型損失[24-25]，后續研究者大都沿用此方法，只是模型的擴散程度與尾跡動量厚度之間的曲線略有不同。本文采用Wright和Miller[26]發展的葉型損失模型，該模型利用葉型損失參數ωpara與當量擴散因子Deq之間的關系曲線計算損失，并考慮了馬赫數帶來的影響，如圖1所示。當量擴散因子和葉型損失參數計算方法分別為

(9)

(10)

圖1 葉型損失參數關系曲線Fig.1 Correlation curves for profile loss parameter

2.2.2 激波損失

高負荷多級高壓壓氣機出于效率和可靠性的要求，一級轉子葉尖進口相對馬赫數一般不高于1.4，因此，本文的激波損失采用經典的Miller模型計算[27]。

Miller模型的關鍵在于如何計算激波前吸力面馬赫數，文獻中大都采用普朗特-邁耶函數計算，但這不可避免地會遇到如何求解吸力面氣流轉角的問題。根據計算校核，本文中激波前吸力面馬赫數MaS的計算表達式為[28]

(11)

式中:Δβ為氣流流經葉片的轉角；Ma1為葉柵進口馬赫數。

2.2.3 二次流損失

葉片表面附面層的徑向流動及葉片之間的通道渦會造成氣流的二次流損失，這種二次流動存在于由葉根到葉尖的整個葉片通道。本文利用Griepentrog[29]方法計算二次流損失，二次流損失與總損失之間的關系為

(12)

式中:fs為二次流損失因子，與葉片彎角θ′及展弦比h/c有關，如圖2所示。

圖2 二次流損失因子曲線Fig.2 Curve for secondary loss factor

2.2.4 端壁損失

端壁區域是壓氣機內部流動最復雜、流動機理認識最淺薄的區域。葉尖間隙流動、端壁附面層及其與葉片表面附面層，乃至激波的相互作用，使得此部分區域的損失很難預測。一種常用的端壁損失模型認為根、尖流面端壁損失為葉型損失的倍數，用端壁損失系數表示，由端壁附面層引起的附加損失沿葉高呈三次曲線分布[11]。此模型的不足之處在于端壁損失系數的給定只憑經驗，而端壁區域損失在總損失的占比較大且在二維正問題計算中熵積累會對計算收斂性產生影響。

為了確定端壁損失系數的值，本文利用2個已有的模型，提出了新的端壁損失計算方法。

(13)

(14)

對于轉子，式(13)和式(14)中的速度采用相對速度。

圖3 端壁損失參數隨擴散因子和間隙的變化Fig.3 Variation of endwall loss parameter with diffusion factor and clearance

已知根、尖流面的端壁損失后，可計算出葉根和葉尖的端壁損失系數KH和KT：

Kj=ωe/ωpj=H,T

(15)

在求出KH和KT的值以后，由端壁附面層引起的附加損失仍根據Hearsey方法計算，但附面層影響范圍控制在距葉根、葉尖30%葉高區域，即端壁損失計算方法為

(16)

(17)

2.3 非設計點損失模型

非設計點損失計算采用Aungier[1]模型，即

ω=ωt[5-4(ξ+2)]ξ<-2

(18)

ω=ωt(1+ξ2) -2≤ξ≤1

(19)

ω=2ωtξξ>1

(20)

式中:ξ=(i-i*)/WIDTH，WIDTH為可用攻角的范圍。

可用攻角范圍的計算參考文獻[18]中的方法。同時，考慮到負攻角狀態葉柵通道阻塞的影響，對可用攻角范圍進行修正。當葉柵喉部阻塞時，由進口到喉部應用連續方程，經變換得到：

(21)

式中：o為葉柵喉道寬度；s為柵距；λ為速度系數；k為比熱比；βcr為阻塞狀態的相對氣流角。

葉柵負攻角狀態的可用攻角范圍應滿足以下要求：

(22)

根據計算經驗，由葉柵阻塞引起的可用攻角范圍修正量不應超過2.5°。一個轉子葉片修正前后的各流面負攻角狀態可用攻角范圍及與S1計算結果的對比如圖4所示，可見修正后與S1計算結果吻合得更好。

圖4 可用攻角范圍對比Fig.4 Comparison of available angle of incidence range

2.4 徑向摻混

為了降低端壁區域的低能量氣流堆積，提高程序的計算穩定性，在S2正問題程序中有必要進行摻混計算，這也符合壓氣機內部的流動規律。本文采用Aungier[1]的方法模擬計算內部摻混流動。

3 S1流面計算的應用

隨著壓氣機設計技術的發展，不斷涌現出新的葉型，比如多圓弧、控制擴散、定制葉型等，但公開文獻中針對這些新設計理念葉型的最小損失攻角與落后角模型尚未涉及，不同的造型方法也給解析形式的經驗模型發展帶來困難。圖5和圖6分別為NACA-65葉型和定制葉型利用經驗關系式計算的最小損失攻角和落后角與S1流面計算結果的對比，其中，最小損失攻角采用式(2)計算，落后角采用式(6)計算，并利用式(7)計算密流比的影響。對于NACA-65葉型，模型計算的最小損失攻角與S1流面計算結果相當，如圖5(a)所示，落后角比S1結果偏高0.5°，相對偏高8.6%，如圖5(b)所示；而對于定制葉型，模型計算的最小損失攻角偏正，如圖6(a)，落后角偏高1.8°，相對偏高37.6%左右，如圖6(b)所示。最小損失攻角、落后角計算方法對NACA-65葉型計算精度相對較高，對新理念葉型的計算精度相對較低。

圖5 NACA-65葉型計算結果Fig.5 Calculation results of NACA-65 blade profile

圖6 定制葉型計算結果Fig.6 Calculation results of tailored blade profile

為了解決豐富的葉型類型與相對匱乏的最小損失攻角、落后角模型之間的不匹配問題，引入了S1流面計算。事實上，利用數值計算已成為當前改進經驗模型的手段之一，例如，Schnoes和Nicke[17]利用葉柵數值模擬重新校驗了傳統的落后角模型、損失模型的系數。本文提出的利用S1流面計算修正S2的計算方法僅限于對最小損失攻角和對應狀態落后角的修正。修正方法為

(23)

(24)

(25)

式中:δ*為最小損失狀態的落后角；下標“M”“H”“T”分別代表S2計算中的葉中、葉根和葉尖流面；“S1”表示S1流面計算結果；C1、C2為S1流面計算與實際流動之間的修正量，即流動的三維效應對落后角的影響。根據與三維計算、試驗特性的比較，C1、C2取值在[0.5°，2.5°]。

圖7為利用S1修正S2計算的一個多級壓氣機的跨聲速轉子、亞聲速轉子和一排靜子葉片的氣流角分布及其與三維計算的對比，該壓氣機采用定制葉型。對比的狀態點均在計算的最高效率或總壓恢復系數附近。對比表明，在10%～90%相對葉高之間的主流區，S2計算的出口氣流角與三維計算結果吻合較好，尤其是亞聲速轉子和靜子葉片。跨聲速轉子和靜子的進口氣流角與三維計算的展向分布不同，但狀態點在最小損失點附近，因此，由攻角偏離最小損失攻角引起的落后角變化量較小。在近端壁區域，三維計算考慮了黏性的影響，壁面采用無滑移邊界；而S2計算出于收斂性的考慮，會避免熵增在端壁區域的堆積，近壁面的子午速度相對主流區不會迅速降低，所以，S2的計算結果與三維計算在端壁區有所不同。

圖7 跨聲速轉子、亞聲速轉子和靜子的氣流角對比Fig.7 Comparison of airflow angles of transonic rotor, subsonic rotor and stator

在目前攻角模型、落后角模型匱乏的狀況下，該方法解決了傳統的攻角、落后角模型與種類繁多、不同設計理念葉型之間的適用性問題，既不受葉型類型的約束，又保留了S2正問題程序中馬赫數、葉尖間隙等因素對最小損失攻角的影響以及攻角變化對落后角的影響。該方法及不同模型的聯合實施有效提升了二維正問題程序的精度和可靠性。

4 失速模型

本研究采用最大總壓升準則來預估多級壓氣機的穩定工作邊界。對任意一條等轉速線，當壓氣機流量減小后，計算總壓比不再升高，則認為壓氣機發生失速。

為了驗證失速模型，統計了4臺多級壓氣機在設計轉速和0.8相對換算轉速的試驗流量-壓比特性。4臺壓氣機依次用C1、C2、C3、C4表示，其平均級壓比如圖8所示，利用堵點流量和失速點壓比作歸一化處理后的特性如圖9所示。由圖9可知，設計轉速特性線的最大壓比點即為失速點；在0.8相對換算轉速下，C1、C2、C3壓氣機的等轉速線均在最高壓比點發生失速，而C4壓氣機的等轉速線在經過最高壓比點(A點)后，試驗仍錄取一個狀態點；經對比，由A點計算的喘振裕度比失速點計算的喘振裕度僅偏低0.5%。對于多級壓氣機的喘振裕度而言，數值計算產生0.5%的偏差是可接受的。

本研究以計算的最高壓比點作為失速點，下文的計算證實此種方法是有效的，預測的多級壓氣機中高轉速失速邊界與試驗結果吻合良好。

圖8 4臺壓氣機的平均級壓比Fig.8 Average stage pressure ratio of four compressors

圖9 歸一化試驗特性Fig.9 Normalized test characteristics

5 校驗計算

采用3個不同負荷水平的、經試驗驗證的壓氣機對發展的方法進行校核計算。壓氣機的負荷水平用葉尖負荷系數表示，計算表達式為

(26)

式中：H為負荷系數；Lu為壓縮功；UT為葉片葉尖旋轉速度。3個壓氣機(依次用A、B、C表示)的平均負荷系數分布如圖10所示。

圖10 壓氣機負荷系數Fig.10 Work coefficient of compressors

為了驗證方法的精度及穩定性，首先利用改進前的程序計算了A壓氣機特性，并與改進后程序的計算結果進行了對比。對于B、C兩臺壓氣機，因為總壓比、負荷提高，改進前的程序未能計算收斂。此外，針對C壓氣機開展了級間參數對比，以驗證級間匹配狀況。

5.1 A壓氣機

A壓氣機為面向燃氣輪機的低壓壓氣機，壓氣機前兩級轉子存在局部超聲速區，一級轉子葉尖進口相對馬赫數為1.16。壓氣機轉、靜子葉片均采用修正圓弧中線、NACA-65葉厚分布。為適應不同葉片的進口馬赫數及流通能力，中弧線采用了不同程度的加彎或減彎。

改進前后的計算特性及其與試驗的對比如圖11所示，其中試驗相對換算轉速nr為1.0和0.9時未錄取到喘振邊界。對比表明，在1.0和0.9轉速程序改進前后計算的流量-壓比、流量-效率特性基本重合；轉速降低后，改進前程序計算的效率降低，計算的流量-壓比線在改進后計算特性線的左下方；值得注意的是，改進前的程序計算的穩定邊界水平較低，本研究發展的計算方法顯著提升了程序的計算穩定性。

圖11 A壓氣機特性Fig.11 Overall performance of Compressor A

改進后程序計算的堵點流量比試驗偏低4%左右，最高壓比水平與試驗相當；各轉速計算的峰值效率與試驗相當，從0.8和0.75轉速來看，計算穩定邊界與試驗較為接近。

5.2 B壓氣機

B壓氣機為E3十級高壓壓氣機第2輪修改后的壓氣機。E3十級高壓壓氣機在1.0換算轉速的設計流量為54.4 kg/s，總壓比為25.0，效率為0.847[30]，一級轉子葉尖進口相對馬赫數為1.353。

壓氣機前4級轉子為跨聲速葉排，采用定制流面葉型的方法設計。5～10級轉子采用修正圓弧中線，5、6級轉子采用多圓弧厚度分布，7～10級采用NACA-65厚度分布。1～4級靜子具有修正的圓弧中線和多圓弧厚度分布，5～9級靜子采用NACA-65葉型。最后一級靜子負責將氣流轉回軸向，彎角較大，采用了特殊設計的葉型。為保證喘振裕度，前6級靜子可調。

計算的壓氣機0.9相對換算轉速以上特性及與試驗的對比如圖12所示。對比可得，計算的流量-壓比特性與試驗吻合良好，計算效率比試驗略高，計算穩定邊界與試驗基本一致。1.0與1.025相對換算轉速的計算效率偏高，這一方面是由于轉子葉尖進口相對馬赫數較高，激波模型計算精度降低，另一方面，文獻[30]認為進口測量儀表及多次高速失速引起的性能退化使測試效率偏低。

圖12 B壓氣機(E3十級壓氣機)特性Fig.12 Overall performance of Compressor B (E3 ten-stage compressor)

0.9相對轉速計算與試驗偏差較大，這除了計算給定的間隙較小外，還與進口導流葉片的性能計算有關，轉速降低后進口導流葉片關角度值很大，出口氣流角及總壓恢復系數計算精度降低。

5.3 C壓氣機

C壓氣機為高負荷多級壓氣機，前3級轉子葉片存在超聲速區域，一級轉子葉尖進口相對馬赫數為1.22。各排葉片均采用多圓弧葉型。在特性計算前，首先計算了各排葉片根、中、尖3個流面的葉型攻角特性，對最小損失攻角和落后角進行了修正，隨后開展了全轉速特性預估。計算特性如圖13所示，計算與試驗吻合良好。

試驗對前5級轉子后的總壓和總溫進行了測量。設計轉速由堵點到失速點的級間壓比、溫升比分布的計算結果與試驗對比如圖14～圖16所示，圖中的壓比、溫升比利用試驗出口總壓比、溫升比進行了無量綱處理。對比表明，堵點、工作點的計算壓比和溫升比與試驗基本吻合。在失速點，計算的溫升比更低，這與圖13中的總特性相應，在設計轉速近失速點，當計算與試驗壓比相當時，計算的效率更高。

圖13 C壓氣機特性Fig.13 Overall performance of Compressor C

圖14 近堵塞點的壓力和溫度對比Fig.14 Comparison of pressure and temperature near choked operating point

圖15 近設計點的壓力和溫度對比Fig.15 Comparison of pressure and temperature near design operating point

圖16 近失速點的壓力和溫度對比Fig.16 Comparison of pressure and temperature near stall operating point

6 結 論

二維計算仍是多級壓氣機設計與分析的核心，經驗、半經驗模型更需隨壓氣機設計技術的進步而不斷更新發展。本文通過對比分析，優化了攻角、落后角及損失計算模塊，提升了經驗、半經驗計算模塊的精度，并提出了利用S1流面修正S2正問題計算的方法，解決了攻角、落后角模型與新理念葉型之間不匹配的問題。最后，利用3個經試驗驗證的多級壓氣機對發展的計算方法進行校驗。結果表明：

1) 除低負荷壓氣機的計算流量偏低外，采用本研究發展的方法計算的壓氣機性能與試驗結果吻合良好，展現出相對較高的計算精度。

2) 本文發展的方法合理，較大幅度提高了程序的計算穩定性，可用于多級軸流壓氣機性能分析，為壓氣機多程序校核提供了新手段。