多级轴流压气机二维性能预测方法

杨晓飞，蒋永松，潘若痴，孟德君

中国航发沈阳发动机研究所，沈阳 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) 本文发展的方法合理，较大幅度提高了程序的计算稳定性，可用于多级轴流压气机性能分析，为压气机多程序校核提供了新手段。