压气机作为航空发动机等机械装置的核心部件, 对整机性能有着至关重要的影响。S1、S2流面理论的提出与计算流体力学(computational fluid dynamics, CFD)的发展为叶轮机械的仿真模拟提供了有力的理论依据和技术支撑, 使压气机的研制降低了对试验的依赖, 大幅缩短了研制周期, 节约了研究成本。
在多级轴流压气机的设计体系中，一维设计和特性分析是压气机设计的关键步骤, 其快速高效的特点, 使其广泛用于参数筛选和初步设计方案的确定。现行的大多数一维计算, 无论用于设计问题[1-2]还是分析及优化问题[3], 均是在叶片排进、出口的平均半径处设计算站, 进行参数迭代使其满足流量、能量守恒。这种一维方案的缺陷在于:①将模型简化过多, 无法反映叶片内部的流动情况; ②完全依赖于经验模型, 无法自动捕获激波及自动预测堵塞工况。为解决传统一维计算方法面临的问题, 研究者们开始将先进CFD方法引入平均半径处的计算中, 并开发了基于时间推进的准二维计算模型。比利时列日大学的Léonard等[4]将准二维模型应用于亚声速压气机上, 对落后角和损失模型进行了自动校准, 得到了与试验数据吻合较好的结果。杜文海[5]在Léonard等工作的基础上, 将准二维模型应用于非定常计算并对失速工况进行了研究, 但目前不稳定边界很难采用准二维模型自动预测, 失速工况仍需要经验模型进行判定。除此之外, 国内外的相关研究并不多见。
本文发展了一种能够模拟叶片通道内部流动的CFD准二维计算模型(以下简称准二维模型), 该模型基于Euler方程, 采用有限体积方法进行空间离散并采用时间推进方法进行求解。同时发展了一种基于时间微分方程的无黏叶片力模型, 用以模拟叶片对气流的偏转作用。以3.5级亚声速轴流压气机PW3S1与1.5级高速轴流压气机作为验证算例, 并将计算结果与试验数据或全三维计算结果进行比较。结果表明, 准二维模型能够较为准确地预测压气机特性并自动捕获激波以及预测堵塞工况, 为多级轴流压气机的快速性能分析提供了新思路。
1 准二维计算基本模型 1.1 控制方程准二维控制方程可由S2流面的准三维Euler方程[6]沿径向平均得到
|
(1) |
式中
|
式中:(z, r, ϕ)为圆柱坐标系坐标分量; m为子午坐标分量; (ξ, ζ, η)为所选取贴体坐标系坐标分量; rξ, mξ, ξm为坐标系变换相关度量; J为坐标系变换的雅各比矩阵。ρ为密度, p为静压, ω为角速度, ε为流体内能。Vm, Vϕ分别为流速在m, ϕ方向的分量。fb, m, fb, ϕ为无黏叶片力在m, ϕ方向的分量。fv, m, fv, ϕ为黏性力在m, ϕ方向的分量。
上述(1)式为叶片区的控制方程, 在叶片区, 选取平均半径处m-ϕ平面内的叶型中弧线为计算流线, 由于压气机叶型较薄, 叶型中弧线可采用叶型吸、压力面坐标平均进行近似处理。在非叶片区, 控制方程也为(1)式, 只是源项中无黏叶片力项等于零。
控制方程采用基于网格中心的有限体积法进行离散, 使用Roe高阶迎风格式进行对流通量计算, 通过5级显式Runge-Kutta进行时间推进, 并进一步嵌入当地时间步长和隐式残差光顺方法加速收敛。
1.2 无黏叶片力模型准二维模型中模拟真实叶片对气流的偏转是通过无黏叶片力实现的, 即无黏叶片力应使叶片区气流流动方向与计算流线相切
|
(2) |
式中:W为气流的相对速度; n为叶片区流线的当地法向量。
由于叶片对气流的偏转作用主要通过叶片表面压力作用产生, 因此通常认为无黏叶片力方向垂直于当地流线。对于无黏叶片力的大小, Léonard等[4]在其准二维模型中通过解速度三角形得到出口速度, 由叶片排做功得到叶片力, 再将叶片力线性分布进叶片区。基于时间推进方法的模型天然具有自动捕获激波的优势, 但基于设计问题的叶片力模型只有在子午马赫数大于1时才能捕捉到激波[7], 传统的压气机中极少出现子午马赫数超声的情况, 因此通常认为该种叶片力计算方法无法自动捕捉激波。Baralon[8]在通流计算中发展了一种基于时间微分的正问题叶片力模型, 无需提前求解叶片排出口速度, 通过规定平均流线, 并采用理论推导的方式阐明了这种正问题叶片力能够捕获通道正激波, 即当叶片区相对马赫数大于1即具备捕获激波的能力。本文基于该种叶片力模型发展了一种适用于准二维计算的叶片力
|
(3) |
式中, τ为虚拟时间步长。这种叶片力的含义为:如果不满足流线相切条件, 即等式右端不为零, 就会在当前时间步自动调整叶片力大小, 调整量用以抵消其不相切的程度, 直到满足相切条件。C为松弛因子, 用来控制叶片力的变化量以保证求解的稳定性, C的量纲为s-2, 即C~Δτ-2, Δτ为当地时间步长, 当地时间步长与当地网格尺寸和速度有关, 因此在不同算例中, 需要对松弛因子进行试取, 一般数量级为0.01或0.1。(3)式可以整合到控制方程(1)中一起进行求解。
1.3 黏性力模型本文的准二维模型基于Euler方程, 因此需要在模型中加入黏性力源项以模拟黏性损失。由于黏性力的作用效果是阻碍气流流动, 因此其方向应与当地流动速度方向相反。借鉴通流计算中的一种广泛适用的黏性力模型[9]
|
(4) |
式中，fv为黏性力的模。黏性力模的大小与熵增有关[14]
|
(5) |
式中，T为当地温度。通过损失模型得到各叶片排的损失系数, 进而得到各叶片排的熵增, 再将熵增转化为相应的熵增分布, 加入黏性力模型, 即可模拟黏性损失效应。
1.4 叶片堵塞由于进行了周向平均, 准二维模型需要使用堵塞因子模拟叶片厚度造成的流道堵塞效应。通过对三维控制方程进行周向平均可得到如下周向堵塞因子[7]
|
(6) |
式中:ϕp, ϕs为叶片压力面、吸力面的ϕ坐标; Nb为叶片数。
Baralon等[8]认为法向堵塞因子比周向堵塞因子更能反映真实叶片厚度对流动的影响, 并给出了法向堵塞因子的计算表达式
|
(7) |
式中:S为栅距; t1, t2为对应点的叶片厚度; γ为流线角度。如图 1所示, BC段的堵塞因子可以用(7)式进行计算, AB、CD段可采用线性分布。
|
| 图 1 法向堵塞因子计算示意图 |
准二维模型能够描述超声速流动并具备捕捉激波的能力, 因而可以自动预测堵塞工况。堵塞因子描述了叶片的厚度分布, 给出了叶片通道最小截面的位置, 因此堵塞因子分布的合理性对堵塞工况预测的准确性有重要影响。
1.5 叶片前、尾缘流线修正处理在准二维模型中, 气流在非叶片区时不受叶片力约束, 而一旦进入叶片区即要在叶片力作用下满足流线相切条件, 这样周向速度会在前缘(LE)发生突变, 产生速度间断并导致非物理熵增。本文采取修改前缘流线使其与来流速度相切的方法, 以避免速度间断的产生。对于尾缘(TE), 由于是从叶片力约束区进入非约束区, 不会产生这样的速度间断, 但是需要加入落后角的影响, 因此也要对尾缘流线进行修正。图 2所示为流线修正示意图, 采用三次样条插值进行修正以使前、尾缘流线满足如下条件
|
(8) |
|
| 图 2 前、尾缘流线修正示意图 |
图 3为无黏流动情况下, 1.5级压气机算例中叶片区流线修正前、后的熵增示意图, 可以看出修正基本消除了叶片前缘的非物理熵增。
|
| 图 3 流线修正前、后熵增 |
Baralon[8]在通流计算中指出, 影响性能计算的关键在于气流本身的偏转而非叶型中弧线的偏转角度, 因此计算的准确性主要取决于落后角模拟的准确性, 流线的修改可能会影响叶片载荷分布、马赫数分布等流动细节参数, 但并不会对总体性能参数计算造成太大影响。
1.6 损失及落后角模型损失和落后角模拟的准确性是准二维计算的关键, 本文所用经验模型如表 1所示。其中, 对于超声速流动, 准二维模型可以自动计算激波造成的熵增, 无需添加额外的激波损失模型。
为使模拟结果更接近真实流动, 对落后角模型考虑径向不均匀性的影响。本文对Robert等试验得出的落后角修正曲线做出线性简化, 如图 4所示, 图中关键节点的数值见文献[17]。对轴向速度沿径向的分布做出如图 5所示的假设[4], 并对落后角进行流量平均, 即
|
(9) |
|
| 图 4 转子、静子中落后角沿径向分布 |
|
| 图 5 轴向速度沿径向分布 |
PW3S1为某八级轴流压气机的后3级[18], 带有进口导叶, 转子和静子叶型为非标准叶型, 中弧线为圆弧, 厚度分布采用NACA 65系列的厚度分布。PW3S1设计点流量为4.3 kg/s, 设计压比为1.357。PW3S1的计算域网格示意图如图 6所示, 为保证流线修正的效果, 在前后缘处对网格进行了加密处理。图 7所示为一个状态点的收敛史。本文对该压气机的85%, 100%和105%转速的整体气动性能进行了计算, 并与试验结果做了对比, 如图 8至9所示。为了将试验值与计算结果进行量化分析, 提出平均误差的计算方法
|
(10) |
|
| 图 6 PW3S1计算网格示意图 |
|
| 图 7 PW3S1某状态点的收敛史 |
|
| 图 8 PW3S1压比特性对比 |
|
| 图 9 PW3S1效率特性对比 |
式中, ei=|xc, i-xe, i|/xe, i×100%为特性图中, 流量相同时对应的准二维计算结果xc和试验结果xe的误差, 每条特性线上n≥10。
从压比特性图可以看出:在85%, 100%和105%转速下的平均误差分别为0.20%, 0.43%和0.94%。从效率特性图可以看出:3个转速下的平均误差分别为0.10%, 0.74%和1.06%。计算所得曲线与试验曲线的变化趋势基本一致, 85%转速下效率预测整体偏低, 而在100%转速与105%转速下效率均偏高。计算结果相比试验值, 在高转速、大流量工况下计算流量偏大, 压比偏高, 效率偏高。在计算过程中发现预测的落后角相比试验值偏小, 这是造成了压比预测偏大、流量预测偏大的直接原因。此外, 作为高压压气机的后三级, 附面层对流动影响较大, 准二维模型将整个径向的流动简化到平均半径处, 即用流动状态较好的中径流动代替流动状态较恶劣的轮毂、机匣处流动, 一定程度上会造成损失估计偏小, 流量、压比估计偏大的情况。
2.2 1.5级高速轴流压气机该算例为德国亚琛工业大学设计的3级跨声速轴流压气机第一级加进口导叶[19], 设计转速为17 000 r/min, 流量13.4 kg/s, 压比1.3, 在中高转速下的大流量工况, 静子中会出现超声速流动以及激波现象。其计算域网格示意图与一个状态点的收敛状况如图 10至11所示。图 12为基于设计问题的叶片力与本文采用的正问题叶片力计算的绝对马赫数对比, 可以看出采用正问题叶片力可以捕获到静子叶片中的激波, 而采用设计问题的叶片力时, 虽然流动超声速, 但静子叶片中马赫数分布光滑连续, 并未捕获到激波。本文对该压气机的68%, 88%和100%转速进行了特性计算, 由于缺少试验数据, 在与仅有的压比-流量试验值对比外, 还与NUMECA的全三维模拟结果进行了对比, 如图 13至14所示。
|
| 图 10 1.5级压气机计算网格示意图 |
|
| 图 11 1.5级压气机某状态点收敛史 |
|
| 图 12 2种叶片力计算的绝对马赫数对比 |
|
| 图 13 1.5级压气机压比特性对比 |
|
| 图 14 1.5级压气机效率特性对比 |
从压比特性图可以看出:与试验值相比, 68%与100%转速下压比的平均误差分别为1.19%和1.06%。与NUMECA计算值相比, 3个转速的平均误差分别为0.72%, 1.04%和3.43%。在68%和100%转速下, NUMECA计算值相对试验值偏小, 可能是由于多级环境下非定常效应明显, 采用的掺混面模型、湍流模型等并不能完全反映实际流动现象, 这也是单通道定常三维计算的局限性。从效率特性图可以看出:与NUMECA计算值相比, 3个转速的平均误差分别为2.68%, 3.01%和3.04%。准二维模型预测的最大效率偏高。此外, 从上述特性图可以看出, 准二维模型可以自动预测堵塞工况, 预测的堵塞流量与全三维结果较为接近, 88%与100%转速下误差分别为0.79%和1.04%。
3 结论本文结论如下：
1) 本文开发的准二维性能分析程序求解快速、稳定性好, 计算残差能在10 s内下降10个数量级, 能为轴流压气机的初始设计分析提供快速的方案筛选工具。
2) 在PW3S1压气机算例中, 准二维模型的总特性预测结果与试验值较为吻合, 压比与效率特性的平均误差均小于1%。
3) 在1.5级压气机算例中, 准二维模型自动捕捉到了激波, 自动预测的堵塞流量与全三维结果相比误差在1%左右。计算的压比特性与试验值的误差在1%左右, 压比和效率特性与全三维计算结果的误差均在3%左右。
| [1] | NICLAS Falck. Axial Flow Compressor Mean Line Design[D]. Sweden: Lund University, 2008 |
| [2] |
张军.多级轴流压气机方案设计与特性计算研究[D].北京: 北京理工大学, 2016 ZHANG Jun. Preliminary Design and Performance Calculation for Multistage Axial Flow Compressor[D]. Beijing: Beijing Institute of Technology, 2016(in Chinese) |
| [3] |
钟勇健.多级轴流压气机一维性能计算与优化研究[D].上海: 上海交通大学, 2014 ZHONG Yongjian. Mean Line Performance Calculation and Optimization of Multi-Stage Axial Compressor[D]. Shanghai: Shanghai Jiaotong University, 2014(in Chinese) http://d.wanfangdata.com.cn/Thesis/D576872 |
| [4] | LÉONARD O, ADAM O. A Quasi-One-Dimensional CFD Model for Multistage Turbomachines[J]. Journal of Thermal Science, 2008, 17(1): 7-20. DOI:10.1007/s11630-008-0007-z |
| [5] | DU W H, LÉONARD O. A Quasi-One-Dimensional CFD Model for Multistage Compressors[C]//Proceedings of ASME Turbo Expo: Turbine Technical Conference & Exposition. Vancouver, 2011: 101-112 |
| [6] | YANG C, WU H, YANG J, et al. Time-Marching Throughflow Analysis of Multistage Axial Compressors Based on a Novel Inviscid Blade Force Model[J]. Journal of Aerospace Engineering, 2019, 233(14): 5239-5252. |
| [7] |
杨金广, 王春雪, 王大磊, 等. 基于时间推进的通流计算方法:现状及展望[J]. 航空学报, 2017, 38(9): 58-70.
YANG Jinguang, WANG Chunxue, WANG Dalei, et al. Time Marching Based Throughflow Method:Current Status and Future Development[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(9): 58-70. (in Chinese) |
| [8] | BARALON S, ERIKSSON L, HALL U. Validation of a Throughflow Time-Marching Finite-Volume Solver for Transonic Compressors[C]//International Gas Turbine & Aeroengine Congress & Exhibition, Stockholm, 1998: 98 |
| [9] | HIRSCH C, WARZEE G. A Finite-element Method for throughflow Calculations in Turbomachines[J]. Journal of Fluid Engineering, 1976, 98(3): 414-419. DOI:10.1115/1.3448342 |
| [10] | CETIN M, UCER A, HIRSCH C, et al. Application of Modified Loss and Deviation Correlations to Transonic Axial Compressors[R]. AGARD R-745, 1987 |
| [11] | LIEBLEIN S. Incidence and Deviation-Angle Correlations for Compressor Cascades[J]. Journal of Fluids Engineering, 1960, 82(3): 584-585. |
| [12] | CREVELING H, CARMODY R. Axial Flow Compressor Computer Program for Calculating Off-Design Performance[R]. NASA, CR-72427, 1968 |
| [13] | KOCH C, SMITH L. Loss Sources and Magnitudes in Axial Flow Compressors[J]. ASME J Eng for Power, 1976, 98(3): 411-424. DOI:10.1115/1.3446202 |
| [14] | HOWELL A. Fluid Dynamics of Axial Compressors[J]. Preceedings of the Institution of Mechanical Engineers, 1945, 153: 441-452. DOI:10.1243/PIME_PROC_1945_153_049_02 |
| [15] | AUNGIER R, FAROKHI S. Axial-Flow Compressors:a Strategy for Aerodynamic Design and Analysis[J]. Applied Mechanics Reviews, 2004, 57(4): 22. |
| [16] | SARAYANAMUTTOO H, ROGERS G, COHEN H. Gas Turbine Theory 6th Edition[M]. Harlow: Pearson Education, 2009. |
| [17] | ROBERTS W B, SEROVY G K, SANDERCOCK D M. Modeling the 3D Flow Effects on Deviation Angle for Axial Compressor Middle Stages[J]. ASME Journal of Engineering for Gas Turbines and Power, 1986, 108: 131-137. DOI:10.1115/1.3239859 |
| [18] | BURDSALL E A, CANAL E J, LYONS K A. Core Compressor Exit Stage Study-I. Aerodynamic and Mechanical Design[R]. NASA Contract Report CR-159714, 1979 |
| [19] | NIEHUIS R, BOHNE A, HOYNACKI A. Experimental Investigation of Unsteady Flow Phenomena in a Three Stage Axial Compressor[J]. Proc IMechE, Part A:J Power and Energy, 2003, 217: 341-348. DOI:10.1243/095765003322315397 |
