压气机进气畸变数值模拟技术研究

安玉戈， 刘火星

北京航空航天大学 能源与动力工程学院 航空发动机气动热力科技重点实验室， 北京 100191

压气机进气畸变数值模拟技术研究

安玉戈， 刘火星*

北京航空航天大学 能源与动力工程学院 航空发动机气动热力科技重点实验室， 北京 100191

发展了一种针对进气畸变条件下的风扇/压气机进行性能预估和稳定性分析的计算方法。首先研究了将叶片作用力简化为体积力源项的建模方法，在此基础上开发出一套基于体积力的三维进气畸变数值模拟程序，使用该程序对NASA Rotor 35在均匀进气、进口存在稳态总压畸变及同时存在总压和总温畸变的流场进行了模拟分析。结果表明，该程序获得的压气机特性及参数分布与雷诺平均Navier-Stokes(RANS)计算吻合得很好，同时正确地模拟出了压气机转子与上游畸变来流的耦合作用及其对压气机性能和稳定工作裕度的影响。

航空航天推进系统； 风扇/压气机； 进气畸变； 体积力模型； 稳定性

在飞机实际飞行过程中，机动飞行、发射武器和进气道不起动等工况都会造成压气机进口流场发生畸变，使其压比、效率和稳定裕度下降，影响发动机的正常工作和飞行安全。未来发动机越来越高的性能指标对压气机的失速裕度提出了更高的要求[1]，要求设计者在设计阶段就考虑兼顾压气机的性能和抗畸变能力。

20世纪70年代有学者基于无黏和不可压缩假设，采用将流动控制方程线性化的方法研究进气畸变对压气机的影响[2-5]，这些研究揭示了上游畸变流场和压气机转子之间有着强烈的耦合作用，这种耦合作用引起的流场参数再分布的规律以及畸变流动的大尺度扰动特性。然而由于做了较多简化，这些理论与工程实际还有相当大的距离，于是随着技术的发展研究者们先后提出了平行压气机模型[6]、激盘模型[7]、周向平均体积力模型、黏性体积力模型以及全环非定常雷诺平均Navier-Stockes(RANS)等方法来模拟畸变流动。平行压气机模型[6]将压气机沿周向分成若干个独立运行通过边界条件相联系的一维子压气机；激盘模型[7]将叶片排对流体的作用集中在一片没有厚度的激盘上。它们都对实际的三维流动做了适当的降维处理，往往无法准确刻画叶片通道中影响畸变特征的复杂流动。近些年来计算机技术的发展使对压气机进行全环非定常数值模拟[8-9]成为了可能，但该方法动辄占用数百颗CPU计算几周时间，显然是一般工程设计所无法接受的。周向平均体积力模型和黏性体积力模型都采用源项来模拟对计算量消耗较大的叶片力, 通过求解三维Euler方程来模拟畸变流场。这两种方法对计算资源的需求很小，并契合了压气机进气畸变大尺度的特征，可以较好地刻画出进气畸变在叶轮机械内部的传递。它们的不同在于周向平均体积力模型将叶片的所有作用力均以源项的形式加入流场，而黏性体积力模型在计算过程中计入叶片的形状，仅仅将黏性作用通过源项的形式加入流场，计算量比周向平均体积力模型高出一个量级。Gong[10]采用周向平均体积力方法将叶片的作用力分解为升力系数和阻力系数作为源项模拟畸变流场；Hale和O’Brien[11]将叶片作用力分解为压力项和惯性项作为源项对进口稳态畸变加以模拟；Chima[12]将叶片作用力分解为与相对速度垂直的非耗散部分和平行于相对速度的耗散部分作为源项建立了一套压气机稳定性计算程序CSTALL。国内郑宁[13]采用周向平均体积力模型模化叶片作用模拟畸变流动，获得了畸变流动的主要特征。Xu[14]采用黏性体积力模型作为研究手段，从单通道RANS计算结果中提取叶片表面的摩擦系数，以源项的方式加入Euler方程求解代替黏性力，成功地模拟出了畸变流场的非定常现象。在以上研究方法中周向平均体积力模型对计算资源和经验简化的需求都较小[14],非常适合工程设计阶段使用。

本文基于周向平均体积力模型的思想发展了一套利用较少的计算资源即可快速预估压气机在进气畸变条件下性能和稳定性的计算程序CDIST。该程序采用RANS计算获得的流场环量和熵增数据来构造叶片力源项，通过求解三维Euler方程来模拟畸变流场。

1 体积力模型

在压气机通道中利用三维无黏非定常Euler方程来求解管道流动，叶片作用以体积力的形式加入到方程源项中，控制方程如式(1)所示：

(1)

Φ=[ΦrΦθΦx]T

(2)

式中：ρ为密度；p为压力；e为内能；r为半径;Vr、Vx和Vθ分别为径向、轴向和周向速度；Φr、Φx、Φθ和Wf分别为径向、轴向和周向的体积力源项以及体积力对气流做的功。在叶片区域这些源项由周向平均体积力模型给出，在叶片排以外的区域，这些源项为0。对控制方程空间离散采用有限体积法，离散精度为二阶，时间离散采用改进的具有三阶精度的Lax-Wendroff格式。

将体积力分解为与相对流动方向平行的力f和垂直于相对流动方向的力F，可以写为

Φ=F+f

(3)

F·V′=0

(4)

(5)

式中：V′为叶片参考系中的相对速度；fθ、fx和fr分别为平行于相对速度的体积力在周向、轴向、径向的分量。Chima[12]指出周向的体积力与流场速度环量密切相关，F无耗散作用，f表征着损耗，可以看做是外界对气流的摩擦力，与熵增存在着对应关系：

(6)

(7)

式中：T为静温；Vm为子午面流动速度；∂m(rVθ)和∂ms分别为速度环量和熵增沿子午面流线的导数。当流场参数和叶片排环量熵增变化特性确定后即可唯一确定叶片作用力。

(8)

(9)

图1 参考点Δ(rVθ)、Δs沿叶展方向的分布Fig.1 Spanwise distribution of Δ(rVθ) and Δs

图2 Δ(rVθ)随mcor和ncor的变化规律Fig.2 Correlation of Δ(rVθ) with mcor and ncor

图3 Δs随mcor和ncor变化规律Fig.3 Correlation of Δs with mcor and ncor

在求解过程中，每当迭代步到达更新体积力源项时，位于叶片排前缘的探测面会计算各个周向位置的折合转速ncor、折合流量mcor，并根据预先导入的RANS数据进行插值，对该周向位置的沿叶展方向的速度环量和熵增分布进行整体缩放，Euler方程求解出某一周向位置的折合流量超过其所参考数据的喘点和堵点范围时，即认为此时压气机已经失稳或发生堵塞，如图4所示。由于这种参数匹配方式考虑了气流切向相对马赫数的影响，可以较Chima[12]的原始模型更准确地刻画气流切向速度对加功量和损失的影响，同时使程序进一步具备了模拟不同叶片通道工作在不同折合转速线上的进口温度畸变的能力。

图4 参数周向匹配示意图Fig.4 Schematic of parameters circumferential coupling

在畸变流场中各个不同的周向位置有着不同的折合转速、折合流量，因而该处的叶片作用力是不同的。周向平均体积力模型通过设立探测面来捕捉这一现象，贴合了畸变流场的特征。另外正如式(6)和式(7)所描述的在此体积力模型中，体积力不仅仅取决于由ncor和mcor控制的环量和熵的分布，同时也和当地流场状况密切相关，故最终收敛的结果也应当是进气畸变扰动充分发展后的结果。同时由于求解器求解的是三维Euler方程，可以利用很少的计算资源刻画出畸变流动强三维、大尺度的特点。

2 算例验证

利用CDIST畸变程序计算了NASA Rotor 35在均匀进气、进口存在稳态总压畸变及同时存在总压和总温畸变条件下的流场以验证其对压气机流场和性能的模拟计算能力。NASA Rotor 35基本参数如表1所示，更详细的信息参阅文献[15]。

表1 NASA Rotor 35基本参数Table 1 Basic parameters of NASA Rotor 35

图5给出了模拟流场的子午面和周向的计算网格，轴向为62个网格，周向为60个网格，径向为30个网格，总共有111 600个网格单元。0～2为3个变量的探测截面，在下文变量的数字编号中0为上游参考截面，1为转子前活动截面，2为转子尾缘截面。网格的生成过程即保证了网格具有严格的轴对称性质，同时网格沿周向和径向也是等距的，保证了网格的均匀性和正交性。这种网格利于使用大的时间步加快收敛速度，对于单转子均匀进气计算利用普通PC十几分钟即可收敛。

图5 NASA Rotor 35计算网格Fig.5 Computational grid of NASA Rotor 35

2.1 均匀进气计算结果

首先使用畸变程序CDIST完成了均匀进气条件下NASA Rotor 35在100%、90%以及70%等转速线的计算。图6对比了畸变程序CDIST、RANS计算以及实验测量获得的NASA Rotor 35流量-总压比和流量-效率特性。采用商用程序CFX进行RANS计算，选用标准的k-ε模型，并通过了网格无关性检验。畸变程序和RANS计算都获得了与实验测量非常接近的流量-总压比特性。对于流量-效率特性畸变程序和RANS计算的结果都普遍比实验值高出1%～2%。可以看出，由于正确地加载了叶片排的加功和损失情况，畸变程序以非常小的计算代价获得了和RANS计算几乎完全一致的压气机特性。

图6 NASA Rotor 35特性曲线Fig.6 Performance curves of NASA Rotor 35

图7将在参考点对采用畸变程序CDIST和RANS计算获得的总温比、总压比及效率沿叶展方向的分布与实验测量值进行了对比，可以看出三者获得的参数分布规律是一致的。对于总温比和总压比计算值和测量值的差别均在2%以内；畸变程序CDIST在叶根处计算精度较高，在叶尖激波较强处抹平了部分气动损失，造成对叶尖效率计算结果偏高。

以上研究结果表明在不同流场位置、不同工况下畸变程序都准确地加载了当地的加功和损失情况，以非常小的计算资源模拟出了和RANS计算相近的结果。

图7 参数沿叶展方向的分布对比Fig.7 Comparison of parameters spanwise distribution

2.2 进口存在稳态总压畸变计算结果

使用畸变程序CDIST模拟计算了NASA Rotor 35在90%转速时流场上游(图5中0处)120°～240°周向范围内存在总压畸变时的压气机流场，进口总压分布如图8所示。图中：pt为总压，下标0为上游参考截面，1为转子前活动截面，2为转子尾缘截面，c代表进口来流无畸变区域，因此pt0为上游平均总压，pt0c为上游无畸变区域总压。使用单台PC进行计算20 min左右即可收敛。

图8 流场进口总压分布Fig.8 Total pressure distribution at inlet

2.2.1 转子对上游畸变流动的影响

在实际工作过程中进气畸变条件下转子对流场的影响向上游传播，能够使上游畸变气流参数在进入风扇叶片排前重新分布。根据平行压气机原理[6]，处于低进口总压区域的叶片通道工作在特性线流量偏小的工况，加功量大、气流通过通道产生的压升也相对较大，因此正如图9给出的截面 1在不同位置时参数周向分布情况，在低进口总压区域叶片排对上游气体的抽吸作用更为强烈，使该处静压下降、气流加速，这样风扇转子和畸变流场的相互作用使风扇进口的空气质量流量趋于均匀。图中：Ux为轴向速度，x为图5中转子上游截面 1的不同轴向位置,c为叶片弦长。

图9 压气机上游静压、质量流量分布 Fig.9 Upstream static pressure and mass flux distribution of compressor

动叶前缘不均匀的静压引起了该处气流的周向流动，使畸变区域两侧的气体涌向中间的畸变区域如图10所示，使在120°周向位置附近的气流攻角变小，在240°周向位置附近的气流攻角变大。图11给出了转子前缘绝对气流角的周向分布。进口绝对气流角的周向不均匀是影响进气总压畸变条件下压气机性能的重要因素，本文模型计算方法较原始的Chima[12]模型计算方法增加了准则参数ncor来捕捉进口绝对气流角不均匀对压气机性能所造成的影响，在下文将进行更为详细的讨论。

图10 转子前缘静压分布Fig.10 Static pressure distribution at the leading edge

图11 转子前缘绝对气流角分布Fig.11 Absolute whirl angle distribution at the leading edge

2.2.2 不同位置叶片通道工况

图12为在进气总压畸变条件下NASA Rotor 35 的流量-总压比特性以及各个周向位置的流量-总压比特性。流量-总压比特性与在90%转速均匀进气的特性偏差不大，这与文献[12]中的结果是一致的。由于流量和攻角状态不一致，不同周向位置的工作特性环绕平均特性点构成了一条包络线，在工作过程中，每个叶片通道在旋转一圈时将依次遍历这些工况。同时还可以看出，点e处附近叶片通道距失速边界较近，这表明进气畸变消耗了失速裕度，有可能带来压气机稳定性问题，CDIST畸变程序正是通过判断不同通道工况与失速边界的距离来评估进气畸变对压气机稳定性带来的影响。

图12 叶片通道工况周向分布Fig.12 Operation condition circumferential distribution

图13～图16给出了转子尾缘总温、总压分布以及转子前缘质量流量和绝对气流角的分布，其中Tt为总温，下标LE表示转子前缘位置。就图12中a、b、c、d和e5个典型通道周向位置工况进行分析，在a点处于畸变区域，进口质量流量小，同时受上游气流涌入畸变区域附加周向速度的影响攻角大，因此加功量最大。在随后的b处虽然并未处于畸变区域，但受前方周向气流的影响，有较大的正攻角，加功量较大，并且由于其前方对应的非畸变区域总压较高，在出口达到了最大总压。到c点质量流量达到最大，前缘绝对气流角持续增大，攻角减小，加功量持续下降。随后由于前缘周向气流对攻角的影响压比持续下降，直至d点压比到达极小值。随后进入畸变区域，随着质量流量的减小与负攻角状态的削弱压比开始回升，到e点质量流量达到最小。以上分析结果和文献[1]中采用128颗CPU进行耗时1～2个月的全环非定常数值模拟得到的结论是完全一致的。正是由于畸变程序采用了mcor和ncor两个准则参数来匹配RANS数据，综合考虑了质量流量和进口气流角对于加功量、损失的影响，因此利用非常小的计算资源较为准确地模拟出了压气机的畸变流动。

图13 转子尾缘总温分布Fig.13 Total temperature distribution at the trailing edge

图14 转子尾缘总压分布Fig. 14 Total pressure distribution at the trailing edge

图15 转子前缘质量流量分布图Fig.15 Mass flux distribution at the leading edge

图16 转子前缘绝对气流角分布Fig. 16 Absolute whirl angle distribution at the leading edge

2.3 进口同时存在总压和总温畸变计算结果

使用CDIST畸变程序模拟计算NASA Rotor 35转子上游90°～210°范围内总温上升25%，150°～270°范围内总压降低9%同时存在着总温畸变和总压畸变的流场情况，如图17所示，以进一步考察程序的性能。

图17 进口流场Fig.17 Inlet boundary condition

图18和图19给出了总压畸变和总温畸变同时存在时转子尾缘的总温和总压分布。在畸变重叠区域，由于进口总压低并处于低的折合转速条件下，出口总压达到全场的最低值，然而相对于总温畸变的其他区域该处处于总压畸变下，折合流量低、加功量大，因此该处的总温达到了全场的最高值；在总压畸变区的逆时针下游区域进口总压高，并工作在高折合转速线上，同时受到正攻角状态的影响，形成了高出口总压区；对于无畸变区域，由于进口总温低、工作在大折合流量工况，因此在中心位置形成了全场出口总温最低的区域。

图18 温度压力组合畸变条件下转子尾缘总温分布Fig.18 Total temperature distribution at the trailing edge with combined pressure temperature distortion

图19 温度压力组合畸变条件下转子尾缘总压分布Fig.19 Total pressure distribution at the trailing edge with combined pressure temperature distortion

3 结 论

1) 在Chima[12]原始模型的基础上进一步考虑进口绝对气流角对转子通道工况的影响，采用叶片前缘质量流量和相对切向速度参数作为匹配准则加载RANS计算获得速度环量和熵增数据构造体积力可以更为准确地模拟出进气畸变条件下压气机流场各地的加功和损失情况。

2) 进气畸变流动是以压气机半径为特征尺度的三维流动，周向平均体积力模型通过模拟不同周向位置叶片通道对不均匀来流的响应来捕捉大尺度畸变扰动的传播，采用求解三维Euler方程的方式描述三维流动，契合了畸变流动的特点，可以得到与全环非定常RANS计算一致的结果。

3) 周向平均体积力模型可以利用非常少的计算资源准确地刻画出畸变流动的主要特征及其对压气机性能和稳定工作裕度的影响，可有效提高工程设计阶段抗畸变能力的分析效率，具有较好的发展前景。

[1] Fidalgo V J, Hall C A, Colin Y. A study of fan-distortion interaction within the NASA Rotor 67 transonic stage. ASME Paper, GT-2010-22914, 2010.

[2] Callahan G M, Stenning A H. Attenuation of inlet flow distortion upstream of axial flow compressors. Journal of Aircraft, 1971, 8(4): 227-233.

[3] Henderson R E, Horlock J H. An approximate analysis of the unsteady lift on airfoils in cascade. ASME Paper， GT-1972-5, 1972.

[4] Greitzer E M, Mazzawy R S. Flow field coupling between compression system components in asymmetric flow. ASME Journal of Engineering for Power, 1978, 100(1): 66-72.

[5] Longley J P, Greitzer E M. Inlet distortion effects in aircraft propulsion system integration. AGARD Lecture Series, 1992.

[6] Mazzawy R S. Multiple segment parallel compressor model for circumferential flow distortion. ASME Journal of Engineering for Power, 1977, 99(2): 288-296.

[7] Hawthorne W R. Non-axisymmetric flow through annular actuator disks: inlet distortion problem. ASME Journal of Engineering for Power, 1978, 100(4): 604-617.

[8] Yao J, Gorrell S D，Wadia A R. A time-accurate CFD analysis of inlet distortion induced swirl in multistage fans. AIAA-2007-5059, 2007.

[9] Gorrell S E, Yao J,Wadia A R. High fidelity URANS analysis of swirl generation and fan response to inlet distortion. AIAA-2008-4985, 2008.

[10] Gong Y F. A computational model for rotating stall and inlet distortions in multistage compressors. Cambrige: Department of Aeronautics and Astronautics, Massachusetts Institute of Technology, 1998.

[11] Hale A, O’Brien W. A three-dimensional turbine engine analysis compressor code (TEACC) for steady-state inlet distortion. Journal of Turbo-machinery, 1997, 120(3): 422-430.

[12] Chima R V. A three dimensional unsteady CFD model of compressor stability. ASME Paper, GT-2006-90040, 2006.

[13] Zheng N. 3D unsteady numerical simulation of fan/compressor with inletdistortion. Journal of Aerospace Power, 2007, 22(1): 60-65. (in Chinese)

郑宁. 风扇进气畸变三维非定常模拟技术研究. 航空动力学报, 2007, 22(1): 60-65.

[14] Xu L. Assessing viscous body forces for unsteady calculations. ASME Paper, GT-2002-30359, 2002.

[15] Reid L, Moore R D. Performance of a single-stage axial-flow compressor with rotor and stator aspect ratios of 1.19 and 1.26, respectively, and with design pressure ratio of 1.82. NASA Technical Paper 1338, 1978.

NumericalSimulationofCompressorwithInletDistortion

ANYuge,LIUHuoxing*

NationalKeyLaboratoryofScienceandTechnologyonAero-Engines,SchoolofJetPropulsion,BeihangUniversity,Beijing100191,China

Anumericaltoolisdevelopedtoevaluatetheperformanceandstabilityoffan/compressorwithinletdistortion.First,themodelingmethodtoreplacethebladerowforceswithdistributedbulkbodyforcesourcetermsisinvestigated.Thenathree-dimensionalcomputationalfluiddynamics(CFD)codeisdevelopedbyaddingthemodelintoanEulersolvertosimulatethefan/compressorflowfieldwithinletdistortion.ANASARotor35flowfieldwithcleaninlet,inletsteadytotalpressuredistortionandcombinedtotalpressuretotaltemperaturedistortionaresimulatedrespectivelywiththecode.ItdemonstratesthattheresultsobtainedbythecodeofNASARotor35withcleaninletagreewellwiththesolutionsofReynoldsaverageNavier-Stokes(RANS).Inthecasewithinletdistortion,thiscodecanobtainthekeyfeatureoftheinteractionoftherotorandupstreamflowfield,anditsinfluenceonrotorperformanceandstallmargin.

aerospacepropulsionsystem;fan/compressor;inletdistortion;bodyforcemodel;stability

2011-10-14；Revised2011-12-09；Accepted2011-12-23；Publishedonline2011-12-281818

URL：www.cnki.net/kcms/detail/11.1929.V.20111228.1818.008.html

NationalNaturalScienceFoundationofChina(51106004)

.Tel.：010-82316418E-mailliuhuoxing@buaa.edu.cn

2011-10-14;退修日期2011-12-09;录用日期2011-12-23; < class="emphasis_bold">网络出版时间

时间：2011-12-281818

www.cnki.net/kcms/detail/11.1929.V.20111228.1818.008.html

国家自然科学基金(51106004)

.Tel.：010-82316418E-mailliuhuoxing@buaa.edu.cn

AnYG,LiuHX.Numericalsimulationofcompressorwithinletdistortion.ActaAeronauticaetAstronauticaSinica,2012,33(9):1624-1632. 安玉戈，刘火星．压气机进气畸变数值模拟技术研究．航空学报,2012,33(9):1624-1632.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

1000-6893(2012)09-1624-09

V231.3; V231.1

A

安玉戈男， 博士研究生。主要研究方向： 叶轮机气体动力学。

Tel: 010-82338139-805

E-mail: yugekl@163.com

刘火星男， 博士， 副教授， 博士生导师。主要研究方向： 叶轮机气体动力学。

Tel: 010-82316418

E-mail: liuhuoxing@buaa.edu.cn