文章快速检索  
  高级检索
基于通流方法的某涡喷发动机整机数值仿真
刘晓恒1, 周成华1, 宋满祥1, 金东海1,2, 桂幸民1,2     
1. 北京航空航天大学 能源与动力工程学院, 北京 100083;
2. 北京航空航天大学 先进航空发动机协同创新中心, 北京 100083
摘要: 采用轴对称方法对带有黏性力的三维Euler方程组进行降维,利用时间推进方法求解,得到适用于航空发动机整机计算的准三维数值仿真软件,并对某涡喷发动机整机进行设计点和非设计点数值模拟。首先,对地面静止状态节流特性进行研究,将计算结果与试验数据对比可知:推力的最大相对误差为-5.1%,单位燃油消耗率的最大相对误差为+4.8%,相对转速为95%时,单位燃油消耗率最低;其次,获取了飞行马赫数为0.7工况下的高度特性以及飞行高度为3 km工况下的速度特性,将计算结果与设计参数对比可知:对于高度特性,推力的最大相对误差为-4.61%,单位燃油消耗率的最大相对误差为+5%,对于速度特性,推力的最大相对误差为-5.83%,单位燃油消耗率的最大相对误差为+5.92%;再次,分别对压气机/涡轮进行部件模拟,预测了发动机的共同工作线;最后,对发动机设计工况下的流场以及气动参数的展向分布进行了分析。
关键词: 叶轮机械    整机仿真    涡喷发动机    通流方法    地面节流特性    高度/速度特性    
Overall simulation of a turbojet engine based on throughflow method
LIU Xiaoheng1, ZHOU Chenghua1, SONG Manxiang1, JIN Donghai1,2, GUI Xingmin1,2     
1. School of Energy and Power Engineering, Beihang University, Beijing 100083, China;
2. Collaborative Innovation Center for Advanced Aero-Engine, Beihang University, Beijing 100083, China
Abstract: Based on the axisymmetric method, three-dimensional Euler equations with viscous force are turned into two-dimensional problems. The problems are solved using time-marching method, and the obtained software can be applied to the overall simulation of aircraft engine. By means of this tool, a turbojet engine on design point as well as off-design point is simulated. First, the throttle characteristics on the ground are firstly studied and are compared with the experimental data. The results reveal that the maximum error for thrust is -5.1% and that of specific fuel consumption is +4.8%. The specific fuel consumption is the smallest at 95% rotational speed. Second, altitude characteristic at flighting Mach number 0.7 and velocity characteristic on the height of 3 km are obtained through this method. The comparison between the designed value and the computation results show that, for altitude characteristic, the maximum errors of thrust and specific fuel consumption are -4.61% and +5%. For velocity characteristic, the maximum errors of both parameters are -5.83% and +5.92%. Third, Co-operating line of the engine is acquired through simulations of compressor and turbine individually. At last, flow field and spanwise distribution of aerodynamic parameters of the engine on design point are analyzed.
Keywords: turbomachinery    overall simulation    turbojet engine    throughflow method    ground throttle characteristic    altitude/velocity characteristic    

随着计算流体力学的发展以及计算机计算能力的提升,全三维Navier-Stokes数值计算越来越多地应用到叶轮机械的设计过程中。但是,这需要性能极高的计算设备以及大量的计算时间,以至于无法将发动机整机全三维数值仿真计算方法应用到发动机整机日常设计过程中[1-4]。另一方面,零维计算也经常应用在航空发动机整机性能计算中,它依赖于各部件的性能特性图[5]。但是,零维计算无法提供气动参数的展向分布结果。因此,通流方法在现代燃气轮机和航空发动机设计过程中仍然是非常重要的一部分,尤其是在初始设计阶段[6]。这种方法不仅可以应用在给定气动参数的设计阶段,还可以应用在给定几何结构分析其性能参数的分析阶段。

子午面通流方法的概念应当追溯到吴仲华的研究[7-8],借助准三维流面概念,推导出叶片间流动以及子午面流动的统一模型。通过研究两组相对流面,三维流动的正问题与反问题均可以得到求解。第1套使用Euler方程的通流程序是由Spurr[9]提出的,他利用一种通流模型以及叶片间流动的求解器建立一种准三维程序,这种方法利用时间推进求解Euler方程。将这种程序应用到跨声导向器案例中,计算结果与三维求解Euler方程计算结果相一致。Yao和Hirsch[10]提出一种方法,将三维计算程序转换为Euler通流方法,他们研究的目的是将二维设计方法与三维分析方法相结合,并且应用了损失系数和落后角的经验关系式。Dawes[11]提出一种方法,对于多级部件计算,可以将通流方法与全三维数值模拟相结合,将所研究部件利用三维黏性计算进行模拟,同时其他部件利用轴对称方法进行建模。通过这种方法,单级叶片排可以在多级的环境下进行设计和分析。在Damle等[12]的研究中,描述了利用Euler模型设计跨声/超声叶轮机械的方法,这种Euler通流方法成功应用在一个超声风扇与一级跨声涡轮部件的设计过程中。

通流方法的另一个应用领域,是在发动机整机的分析过程和初始设计阶段。由Ivanov等[13-15]提出的通流模型是应用于发动机分析问题中的,他们阐述了在通流计算中描述发动机工作过程的现代数学模型。该模型使用单调保持的二阶Godunov格式,这种格式可以准确求解Riemann问题。在Petrovic等[16-17]的研究中,提出一种基于分别求解压气机与涡轮的求解程序进而求解整机流场的方法。利用简单的燃烧室模型以及二次流动模型来连接压气机与涡轮流道,它可以将发动机整机进行自动匹配。基于以上模型对一台双轴燃气轮机进行整机数值模拟,获得了发动机流场。

国内相关单位对发动机整机计算进行了初步的研究。施发树和刘兴洲[18]、黄家骅等[19]对小型双涵道涡扇发动机进行二维稳态数值模拟,采用Godunov格式求解带黏性力项的非定常Euler方程组。冯国泰等[20]讨论了适用于发动机叶轮部件多场耦合数值仿真计算的统一数学物理模型,分析了包括发动机三维多功能数值仿真数学模型、精度可靠性与并行算法、发动机仿真软件平台的框架等在内的6个关于建立航空发动机数值仿真试验台的关键技术问题。曹志鹏等[21]借助二维仿真软件对某小型涡喷发动机设计点进行数值仿真,对其性能参数及流场进行分析。万科等[22]借助周向平均降维方法对某高通流风扇/增压级进行了性能分析,并与三维数值模拟结果进行对比分析,为航空发动机整机通流计算奠定了基础。

借助通流方法不仅可以获得叶片内部的流场信息,同时可以获得叶片间的匹配信息,这在发动机初始设计阶段是非常重要的。本文利用通流计算软件,对某涡喷发动机设计点与非设计点进行性能与流场研究。首先,对发动机地面静止状态节流特性进行计算,将计算结果与试验数据和设计参数进行对比分析,以此来验证该通流软件计算的可靠性;其次,为了探索发动机的高度/速度特性,进行了相关工况的数值仿真,将计算结果与设计值进行对比;再次,对压气机/涡轮进行部件仿真,获取了发动机共同工作线;最后,对发动机设计工况下的流场参数以及展向分布进行了分析。

1 通流计算方法 1.1 控制方程组

对于柱面坐标系(z, r, φ)下带有黏性力的非定常Euler方程组,利用贴体坐标系(ξ, η, ζ)进行坐标转换,可以得到如下形式的控制方程组:

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {\frac{{r\mathit{\boldsymbol{\bar U}}}}{J}} \right) + \frac{\partial }{{\partial \xi }}\left[ {\frac{r}{J}\left( {\mathit{\boldsymbol{\bar F}}{\xi _z} + \mathit{\boldsymbol{\bar G}}{\xi _r} + \mathit{\boldsymbol{\bar H}}\frac{1}{r}{\xi _\phi }} \right)} \right] + \\ \;\;\;\;\;\;\;\;\;\;\frac{\partial }{{\partial \eta }}\left[ {\frac{r}{J}\left( {\mathit{\boldsymbol{\bar F}}{\eta _z} + \mathit{\boldsymbol{\bar G}}{\eta _r} + \mathit{\boldsymbol{\bar H}}\frac{1}{r}{\eta _\phi }} \right)} \right] + \\ \;\;\;\;\;\;\;\;\;\;\frac{\partial }{{\partial \zeta }}\left[ {\frac{r}{J}\left( {\mathit{\boldsymbol{\bar F}}{\zeta _z} + \mathit{\boldsymbol{\bar G}}{\zeta _r} + \mathit{\boldsymbol{\bar H}}\frac{1}{r}{\zeta _\phi }} \right)} \right] = \frac{{\mathit{\boldsymbol{\bar h}}}}{{J}} \end{array} $ (1)
 

式中:

$ \mathit{\boldsymbol{\bar U}} = \left[ {\begin{array}{*{20}{c}} \rho \\ {\rho u}\\ {\rho v}\\ {\rho w}\\ e \end{array}} \right], \;\;\;\;\mathit{\boldsymbol{\bar F}} = \mathit{\boldsymbol{\bar F}}\left( {\mathit{\boldsymbol{\bar U}}} \right) = \left[ {\begin{array}{*{20}{c}} {\rho u}\\ {\rho {u^2} + p}\\ {\rho uv}\\ {\rho uw}\\ {\left( {e + p} \right)u} \end{array}} \right] $
$ \mathit{\boldsymbol{\bar G}} = \mathit{\boldsymbol{\bar G}}\left( {\mathit{\boldsymbol{\bar U}}} \right) = \left[ {\begin{array}{*{20}{c}} {\rho u}\\ {\rho uv}\\ {\rho {v^2} + p}\\ {\rho vw}\\ {\left( {e + p} \right)v} \end{array}} \right] $
$ \mathit{\boldsymbol{\bar H}} = \mathit{\boldsymbol{\bar H}}\left( {\mathit{\boldsymbol{\bar U}}} \right) = \left[ {\begin{array}{*{20}{c}} {\rho w}\\ {\rho uw}\\ {\rho vw}\\ {\rho {w^2} + p}\\ {\left( {e + p} \right)w} \end{array}} \right] $
$ \mathit{\boldsymbol{\bar h}} = \left[ {\begin{array}{*{20}{c}} {r\dot m}\\ {r\dot m{V_z} + {f_z}}\\ {r\dot m{V_r} + {f_r} + p + \rho {{\left( {w + \omega r} \right)}^2}}\\ {r\dot m{V_\varphi } + {f_\varphi } - \rho v\left( {w + 2\omega r} \right) - {r^2}\rho \frac{{{\rm{d}}\omega }}{{{\rm{d}}t}}}\\ {r\dot mH' + {f_h} + {\omega ^2}{r^2}\rho v - {r^2}\rho w\frac{{{\rm{d}}\omega }}{{{\rm{d}}t}}} \end{array}} \right] $

其中:fzfrfφ为黏性力分量;fh为黏性损失;ρ为密度;uvw分别为轴向、径向、周向速度;p为静压;e为比总内能;ω为在参考相对坐标系内的周向角速度;$\dot m$为抽吸气单位体积流量;VzVrVφ为抽吸气速度分量;H′为抽吸气总焓;ξz代表偏导数$\partial \xi /\partial z$,其余类似符号含义与之类似;J为雅克比行列式,其表达式为

$ J = \left| {\begin{array}{*{20}{c}} {{\xi _z}}&{{\xi _r}}&{{\xi _\varphi }}\\ {{\eta _z}}&{{\eta _r}}&{{\eta _\varphi }}\\ {{\zeta _z}}&{{\zeta _r}}&{{\zeta _\varphi }} \end{array}} \right| $

利用状态方程将式(1)封闭,其表达式为

$ p = \rho RT $ (2)
 
$ e = \rho \left[ {\varepsilon + \frac{1}{2}\left( {{u^2} + {v^2} + {w^2}} \right)} \right] $ (3)
 
$ \varepsilon = \varepsilon \left( T \right) = \int\limits_{{T_0}}^T {{c_v}\left( \tau \right){\rm{d}}\tau + {\rm{Const}}} $ (4)
 

式中:R为气体常数;T为温度;ε为比内能;T0为参考温度;cv(τ)为定容比热容;Const表示常数。

对于叶轮机械,可以将计算域分为两类:包含叶片区域和不含叶片区域。

对于不含叶片区域(如叶片排间隙和燃烧室区域),可认为流动是轴对称的,因此存在如下假设:

$ \zeta \equiv \varphi , \;\;\;\;{\xi _\varphi } = {\eta _\varphi } = 0, \;\;\;\;\frac{\partial }{{\partial \zeta }} \equiv 0 $ (5)
 

借此,可将控制方程组式(1)推导为

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {\frac{{r\mathit{\boldsymbol{\bar U}}}}{J}} \right) + \frac{\partial }{{\partial \xi }}\left[ {\frac{r}{J}\left( {\mathit{\boldsymbol{\bar F}}{\xi _z} + \mathit{\boldsymbol{\bar G}}{\xi _r}} \right)} \right] + \\ \;\;\;\;\;\;\;\;\frac{\partial }{{\partial \eta }}\left[ {\frac{r}{J}\left( {\mathit{\boldsymbol{\bar F}}{\eta _z} + \mathit{\boldsymbol{\bar G}}{\eta _r}} \right)} \right] = \frac{{\mathit{\boldsymbol{\bar h}}}}{J} \end{array} $ (6)
 

对于存在叶片的计算区域,可以认为ζ=Const定义的曲面为相邻叶片间的流面,这样存在如下关系:

$ u{\zeta _z} + v{\zeta _r} + w\frac{1}{r}{\zeta _\varphi } \equiv 0 $ (7)
 

基于此,可将坐标(ξ, η)选取为满足如下关系式的参数:

$ \left\{ \begin{array}{l} {\xi _z}{\zeta _z} + {\xi _r}{\zeta _r} + \frac{1}{{{r^2}}}{\xi _\varphi }{\zeta _\varphi } = 0\\ {\eta _z}{\zeta _z} + {\eta _r}{\zeta _r} + \frac{1}{{{r^2}}}{\eta _\varphi }{\zeta _\varphi } = 0 \end{array} \right. $ (8)
 

最终,应用于包含叶片区域的控制方程推导为

$ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {\frac{{r\mathit{\boldsymbol{\bar U}}}}{J}} \right) + \frac{\partial }{{\partial \xi }}\left[ {\frac{r}{J}\left( {\mathit{\boldsymbol{\bar F}}{\xi _z} + \mathit{\boldsymbol{\bar G}}{\xi _r} + \mathit{\boldsymbol{\bar H}}\frac{1}{r}{\xi _\varphi }} \right)} \right] + \\ \frac{\partial }{{\partial \eta }}\left[ {\frac{r}{J}\left( {\mathit{\boldsymbol{\bar F}}{\eta _z} + \mathit{\boldsymbol{\bar G}}{\eta _r} + \mathit{\boldsymbol{\bar H}}\frac{1}{r}{\eta _\varphi }} \right)} \right] = \frac{{\mathit{\boldsymbol{\bar h}}}}{J} + {{\mathit{\boldsymbol{\bar h}}}_1} \end{array} $ (9)
 

式中:

$ {{\mathit{\boldsymbol{\bar h}}}_1} = \left[ {\begin{array}{*{20}{c}} 0\\ { - \frac{{r{\zeta _z}}}{J}\left( {\frac{{\partial p}}{{\partial \zeta }}} \right) - p\frac{\partial }{{\partial \zeta }}\left( {\frac{r}{J}{\zeta _z}} \right)}\\ { - \frac{{r{\zeta _r}}}{J}\left( {\frac{{\partial p}}{{\partial \zeta }}} \right) - p\frac{\partial }{{\partial \zeta }}\left( {\frac{r}{J}{\zeta _r}} \right)}\\ { - \frac{{r{\zeta _\varphi }}}{J}\left( {\frac{{\partial p}}{{\partial \zeta }}} \right) - p\frac{\partial }{{\partial \zeta }}\left( {\frac{{{\zeta _\varphi }}}{J}} \right)}\\ 0 \end{array}} \right] $

由于航空发动机各部件空间几何结构的复杂性,为了保证计算域边界处的计算精度,本程序采用贴体坐标系对控制方程进行转换求解[23]。关于曲线坐标系以及轴对称降维方法的详细介绍可以参考文献[15]。

在数值模拟过程中,需要给定以下边界条件:转速,进口总温/总压,大气温度/压力,出口背压,燃油流量,燃油燃烧热值以及对冷却空气提取和排放位置的设定。在计算过程中,对流动过程中真实的三维流动细节进行建模模拟,比如:黏性损失,间隙溢流,流道中通过封严装置的泄漏流,冷却空气的提取以及排放。

1.2 通流模型

整机仿真程序求解的是二维欧拉方程, 因此损失及落后角模型是决定计算精度的关键因素。

包含在式(1)中的(fz, fr, fφ)是表征方程黏性力的3个参数:

$ \left( {{f_z}, {f_r}, {f_\varphi }} \right) = - \frac{{r\mathit{\Phi }}}{{{u^2} + {v^2} + {w^2}}}\left( {u, v, w} \right) $ (10)
 

式中:Φ用以模拟黏性损失。

利用σ描述由于黏性损失导致的熵增,则存在如下关系:

$ \mathit{\Phi } = \rho T\frac{{{\rm{d}}\sigma }}{{{\rm{d}}t}} $ (11)
 

函数σ随空间变化,并且与熵变化是同步的,损失值的大小由试验数据或者经验关系式计算得到。

根据部件结构不同,叶轮机械损失模型可以分为轴流压气机、离心压气机、轴流涡轮及径向涡轮等类型。在仿真程序中,损失模型可以分为设计点损失模型和非设计点损失模型。在设计点,总损失由叶型损失、激波损失以及其他损失构成,前面二者可以准确计算得到,其他损失需要利用损失系数进行修正得到。在非设计点, 除了设计总损失之外, 还需要确定由攻角、端壁、尾缘等引起的附加损失修正系数, 与设计损失一起构成非设计点损失。最终,估算得到的损失会转化为熵增,计入黏性叶片力中。对于轴流压气机,其使用模型包括:Koch-Smith叶型损失模型[24]、Creveling非设计点损失模型[25]、Howell端壁损失模型[26]、Aungier叶尖间隙泄漏损失模型[27];对于离心压气机,使用Galvas损失模型[28];对于轴流涡轮,使用AMDCKO损失模型[29]

程序中落后角模型用于计算设计工况下落后角,在非设计工况需要利用系数对落后角进行修正。对于设计工况和非设计工况,分别采用Carter落后角模型[30]和Creveling落后角模型[25]

2 发动机结构以及试验系统

本文所研究发动机模型在子午平面内的简图如图 1所示,其结构包括:一级轴流压气机、一级离心压气机、折流燃烧室以及一级轴流涡轮。

图 1 涡喷发动机的几何结构 Fig. 1 Geometry of turbojet engine

测试发动机台架如图 2所示,台架采用支撑式摇臂结构,它具有2个底座支撑点以及1个具有半圆形结构的悬挂装置,运动的试验台架通过2根轴和4个摇臂与固定的试验台架相连。当发动机产生推力时,移动的台架向前移动,然后与固定的台架相撞,固定台架上的压力传感器即可测量得到发动机产生的推力。

图 2 测试发动机与测试台架 Fig. 2 Test engine and bench

借助该试验系统,可以测量得到如下数据:推力、燃油流量、进口总压、离心压气机出口总压以及涡轮后的总温。利用这些试验数据,可以对地面工况数值仿真结果进行对比分析进而考核准三维数值仿真软件的计算准确性。

3 计算结果分析与讨论 3.1 地面静止状态节流特性分析与校核

首先,对航空发动机的地面节流特性进行计算,并将其与试验数据以及初始设计值进行对比,以考核软件计算的精度。需要说明的是,该软件计算结果是振荡收敛的,对于收敛结果采用算术平均处理方法,将参数标准差作为误差带进行作图。图 3(a)~图 3(d)分别为发动机推力、单位燃油消耗率、组合压气机压比以及涡轮后总温随转速的变化规律(图中给出了试验数据和计算结果的方差)。其中,试验过程中发动机选取4个相对转速作为稳定状态,分别是80%、90%、95%(额定转速)、98%,因此对这4个工况进行考核分析。

图 3 地面节流特性 Fig. 3 Throttle characteristics on ground

图 3(a)所示,当转速增加时,发动机的推力随之增加。通过计算结果与试验数据和设计值对比发现,3组数据的变化趋势是一致的。与试验数据对比,最大相对误差是-5.1%;与设计值进行对比,最大相对误差是-3.6%。如图 3(b)所示,当转速增加时,发动机的单位燃油消耗率随之降低。与试验数据对比,最大相对误差是+4.8%;与设计参数进行对比,最大相对误差是+4.7%。另外,在95%转速工况下,发动机的单位燃油消耗率最低,这也说明了飞机处于巡航状态时为何选择95%转速作为巡航速度。

图 3(c)图 3(d)展示了压气机压比和涡轮后总温随转速的变化曲线,随着发动机转速的增加,压气机的压比以及涡轮后总温都随之增加。同时可以看到,对于压气机压比,计算结果相较于试验数据较高,这导致在计算结果中涡轮的输出功率较高,因此,对于涡轮后总温,计算结果低于试验数据。

通过上述分析,验证了该计算方法对发动机模型计算的准确性,进而可以对其他工况进行数值模拟。

3.2 高度/速度特性

航空发动机的飞行试验耗费巨大,为了研究发动机的高度/速度特性,利用该通流方法计算得到了发动机在不同高度/速度的特性参数,并将计算结果与初始设计参数进行了对比。

3.2.1 高度特性

进行发动机高度特性数值模拟时,工况设置为:相对转速为100%,飞行马赫数为0.7,飞行高度分别为3、6、9、12、15、18 km。图 4(a)为发动机的相对推力与高度的关系曲线,图 4(b)为发动机的单位燃油消耗率与高度的关系曲线(计算结果均利用100%转速工况下设计参数进行了无量纲化)。通过两组数据可知:计算结果与设计参数的趋势是一致的。对比计算结果与设计值,推力的最大相对误差位于高度3 km处,数值为-4.61%;单位燃油消耗率的最大相对误差出现在同一高度,数值为+5.0%。

图 4 高度特性 Fig. 4 Altitude characteristics

对于曲线变化趋势的解释如下:当飞行高度低于11 km时,随着飞行高度的增加,温度、压力以及密度均降低;同时,相较于温度,压力、密度降低得更快。因为飞行马赫数是给定的,所以当高度增加时,进口流量降低。发动机的推力同时取决于流量和单位推力,流量起主要作用。因此,发动机推力随飞行高度的增加而降低。

当讨论发动机的单位燃油消耗率时,将发动机的单位推力作为一个考虑因素。其中,单位燃油消耗率的定义为

$ {\rm{Sfc}} = \frac{{3\;600f}}{{{F_{\rm{s}}}}} $ (12)
 

式中:f为油气比;Fs为发动机的单位推力。随着飞行高度增加,当高度低于11 km时,大气温度随之降低,当飞行高度高于11 km时,大气温度维持不变。因为发动机物理转速维持不变,在飞行高度低于11 km时,折合转速随高度增加而增加,当高度大于11 km时,折合转速不再发生变化。结合压气机特性可知,随着高度的增加,压气机压比先增大然后维持不变。假设涡轮导向器以及尾喷管均处于堵塞状态,则涡轮的膨胀比维持不变。因此,在飞行高度增加的过程中,发动机出口的静压先升高后维持不变,产生的单位推力先增大后不变,单位燃油消耗率先降低然后维持不变。

3.2.2 速度特性

为了研究发动机的速度特性,在不同飞行马赫数工况下,对发动机进行数值模拟。计算工况设置为:物理转速保持不变,飞行高度为3 km,飞行马赫数的变化范围为0.2~0.8,间隔设置为0.1。图 5(a)为发动机推力与飞行马赫数的关系曲线,图 5(b)为单位燃油消耗率与飞行马赫数的关系曲线。由图可知计算结果与设计参数两组数据的变化趋势是一致的;对于计算精度,推力和单位燃油消耗率的最大相对误差分别是-5.83%和+5.92%,并且二者最大值均出现在马赫数为0.8时。

图 5 速度特性 Fig. 5 Velocity characteristics

呈现该变化趋势的原因如下:一方面,当飞行高度保持不变时,随着飞行马赫数的增加,进气道空气流量随之增加;另一方面,进口总温的增加使发动机的折合转速降低,压气机的工作点在共同工作线上向低压比方向移动,进而发动机出口的压力降低,单位推力降低。在流量与单位推力这两个因素的共同作用下,发动机的推力随着飞行马赫数的增大首先降低然后升高。根据式(12),随着单位推力的减少,发动机的单位燃油消耗率增加。

3.3 共同工作线预测

借助该通流方法可以首先对发动机压气机/涡轮进行部件计算,进而获得发动机的共同工作线,发动机转速的变化范围为75%~100%。在进行压气机和涡轮匹配时,利用流量和功率两个参数。首先,在各个转速,计算得到压气机的特性曲线;然后,计算得到涡轮的特性曲线。相同转速下压气机与涡轮特性曲线的交点即为二者具有相同流量和功率的工作点,也就是发动机的共同工作点。需要说明的是,在进行功率匹配时,涡轮的输出功率有一部分用于克服涡轮轴的损失,该部分占涡轮输出总功率的1%。得到发动机共同工作点的流量和功率后,对应得到无量纲密流和组合压气机压比,即可在压气机特性线上绘制发动机共同工作线,如图 6所示。

图 6 发动机的共同工作线 Fig. 6 Workline of engine

图 6可知:对于无量纲密流,计算结果与设计参数相对误差最大值出现在相对转速为75%位置,相对误差为+5.4%;对于组合压气机压比,计算结果与设计参数相对误差最大值在相对转速为80%位置,相对误差为+9.4%。考虑到这种计算方法的便捷性,相对误差是可以接受的。另一方面,在75%转速处密流和组合压气机压比的相对误差相对于其他转速较高的原因是:在远离100%转速的工况,计算使用的损失和落后角模型相较于100%转速处不够准确。

3.4 整机流场分析

相较于零维/一维计算,通流计算方法不仅可以获得发动机的性能参数,同时可以得到发动机主流道内各位置气动参数,进而获取气动参数的展向分布,为发动机部件匹配提供更详细的信息。本部分对发动机整机地面设计点工况进行计算,并对流场参数进行分析。计算机具有8个处理器,处理器型号为Intel(R) Core(TM) i7 CPU 870@2.93 GHz。整机计算网格数量为8 195,完成200 000步迭代计算耗时为2 h。图 7给出了发动机整机计算流量和推力迭代史,从图中可知,计算结果呈现振荡收敛,因此,需对气动参数进行算术平均处理。

图 7 整机计算参数迭代史 Fig. 7 Iteration history of parameters in simulation ofengine

借助100%转速设计工况下涡轮功率,对计算所得压气机和涡轮功率进行无量纲化,并对振荡结果进行算术平均可得:组合压气机功率为0.960 2,涡轮功率为0.975 3,二者差值占涡轮输出功率的1.55%,因此认为压气机与涡轮功率平衡。表 1给出了设计工况下性能参数的相对误差,其中相对误差的定义为:(计算结果-设计参数)/设计参数。通过表中数据可知,各部件流量的相对误差均小于5%,各部件压比的相对误差均小于8%,基于这些数据可以认为计算结果是合理的。

表 1 性能参数对比 Table 1 Comparison of performance parameters
部件相对误差/%
流量压比
轴流压气机-1.136.42
离心压气机-1.07-7.34
轴流涡轮-3.966.52

图 8为气动参数的流场等值线图。其中,图 8(a)为相对马赫数的等值线图,观察可知,在轴流压气机转子叶片尖部区域,存在马赫数大于1的区域,可以证实该轴流压气机转子叶片是跨声转子。在涡轮进口导流叶片中,同样存在马赫数大于1的区域,这与涡轮叶片设计理念是一致的:超声区域存在于涡轮进口导叶出口根部。在轴向扩压器出口,存在低马赫数区域。图 8(c)为总温的等值线图,可以看到经过风斗后气体的总温降低,这与在折流燃烧室中借助风斗引入掺混冷气起到冷却高温气体的作用是一致的。

图 8 气动参数等值线图 Fig. 8 Contours of aerodynamic parameters

在获得了发动机全流场的气动参数后,即可获得发动机在某一计算站气动参数的展向分布。例如:借助转子进出口气动参数得到转子叶片的效率;借助静子叶片进出口的压力可以获得静子叶片的总压恢复系数。对于整级叶轮机械而言,可以获得其级参数(总压升/降、热力反力度等)。这些级参数的定义为

转子叶片排的效率:

$ {\eta ^ * } = \frac{{h\left( {T_{2{\rm{is}}}^*} \right) - h\left( {T_1^*} \right)}}{{h\left( {T_2^*} \right)h\left( {T_1^*} \right)}} $ (13)
 

式中:h为焓值;T1*为转子叶片进口总温;T2*为转子叶片出口总温;T2is*s(p1*, T1*) = s(p2*, T2is*)计算得到,其中:s为熵;p1*p2*分别为转子叶片进、出口总压。

热力反力度:

$ {\mathit{\Omega }_{\rm{K}}} = \frac{{h\left( {{T_{2{\rm{is}}}}} \right) - h\left( {{T_1}} \right)}}{{h\left( {{T_{3{\rm{is}}}}} \right) - h\left( {{T_1}} \right)}} $ (14)
 

式中:下标1表示级进口位置;下标2表示级叶片排轴向间隙位置;下标3表示级出口位置;下标is表示对应的温度为绝热温度。另外,式(14)中所有的温度均为静温。

以轴流压气机级为例,计算得到气动参数的展向分布如图 9所示。由图 9可知:①转子叶片效率及静子叶片总压恢复系数呈现出叶尖较低、叶中较高的现象,并且二者在叶片根部10%展高范围内均存在降低的趋势;②对于级压比及热力反力度,它们呈现出随叶高增加而升高的趋势。借助这些参数,设计人员即可了解在不同位置叶片的设计结果,同时借此得到相邻叶片排的匹配效果,在发动机的初始设计阶段这些参数具有重要意义。

图 9 气动参数的展向分布 Fig. 9 Spanwise distribution of aerodynamic parameters
4 结论

基于轴对称通流方法,对某涡喷发动机设计点及非设计点进行了整机数值仿真,对其性能及流场参数进行了对比分析,得到如下结论:

1) 对发动机地面静止状态节流特性进行数值仿真,发动机转速设定为80%、90%、95%及98%,获得了推力、单位燃油消耗率、压气机压比以及涡轮后总温随转速的变化规律,并将计算结果与发动机整机试验数据进行对比。发动机推力的最大相对误差为-5.1%,单位燃油消耗率的最大相对误差为+4.8%。另外,计算得到在95%转速时发动机单位燃油消耗率最低,这也说明了飞机在巡航状态下为何选择95%转速作为巡航速度。

2) 对发动机在飞行马赫数0.7处进行高度特性数值模拟,将计算结果与设计参数对比,推力的最大相对误差为-4.61%,单位燃油消耗率的最大相对误差为+5.0%,二者均位于3 km处。对发动机在3 km位置进行速度特性模拟,推力和单位燃油消耗率的最大相对误差分别是-5.83%和+5.92%,二者均位于马赫数为0.8时。

3) 分别进行压气机/涡轮部件模拟,利用流量和功率参数进行整机匹配,获取发动机的共同工作线,发动机转速范围为75%~100%。将计算结果与设计值进行对比,各转速流量相对误差小于6%,组合压气机压比相对误差小于10%。

4) 较于零维/一维计算,通流计算除了获取性能参数外,还可以得到发动机流场。对轴流压气机各气动参数的展向分布特征进行分析,可为部件匹配提供技术支持。

参考文献
[1] TURNER M G, NORRIS A, VERES J. High fidelity 3D simulation of the GE90: AIAA-2003-3996[R]. Reston, VA: AIAA, 2003.
[2] ADAMCZYK J J, MULAC R A, CELESTINA M L. A model for closing the inviscid form of the average-passage equation system[J]. Journal of Turbomachinery, 1986, 108(2): 180-186.
Click to display the text
[3] ADAMCZYK J J. Aerodynamic analysis of multistage turbomachinery flows in support of aerodynamic design[J]. Journal of Turbomachinery, 2000, 122(2): 189-217.
Click to display the text
[4] LIU N S. On the comprehensive modeling and simulation of combustion systems: AIAA-2001-0805[R]. Reston, VA: AIAA, 2001.
[5] SIMON J F. Contribution to throughflow modelling for axial flow turbomachines[D]. Liege: University of Liege, 2007: 7-18.
[6] HORLOCK J H, DENTON J D. A review of some early design practice using computational fluid dynamics and a current perspective[J]. Journal of Turbomachinery, 2005, 127(1): 5-13.
Click to display the text
[7] WU C H. A general through-flow theory of fluid flow with subsonic or supersonic velocity in turbomachines of arbitrary hub and casing shapes: NACA-TN-2302[R].Washington, D.C.: NACA, 1951.
[8] WU C H. A general theory of three-dimensional flow in subsonic and supersonic turbomachines of axial-, radial-, and mixed-flow type: NACA-TN-2604[R]. Washington, D.C.: NACA, 1952.
[9] SPURR A. The prediction of 3D transonnic flow in turbomachinery using a combined throughflow and blade-to-blade time marching method[J]. International Journal of Heat Fluid Flow, 1980, 2(4): 189-199.
Click to display the text
[10] YAO Z, HIRSCH C H. Throughflow model using 3D Euler or Navier-Stokes solver[C]//1st European Conference on Turbomachinery Fluid Dynamics and Thermodynamics, 1995: 51-61.
[11] DAWES W N. Toward improved throughflow capability:The use of three dimensional viscous flow solvers in a multistage environment[J]. Journal of Turbomachinery, 1992, 114(1): 8-17.
Click to display the text
[12] DAMLE S V, DANG T Q, REDDY D R. Throughflow method for turbomachines applicable for all regimes[J]. Journal of Turbomachinery, 1997, 119(2): 256-262.
Click to display the text
[13] IVANOV M J, MAMAEV B I. United modeling of working process in aircraft gas turbine engines: GT2008-50185[R]. New York: ASME, 2008.
[14] IVANOV M J, NIGMATULLIN R Z. Quasi-3D numerical model of a flow passage of the aviation gas turbine engines[C]//10th International Symposium on Air Breathing Engines (ISABE), 1991: 299-305.
[15] NIGMATULLIN R Z, IVANOV M J. The mathematical models of flow passage for gas turbine engines and their components: AGARD-LS-198[R]. Paris: AGARD, 1994.
[16] PETROVIC M V, RAHMAN A A, WIEDERMANN A. A quick method for full flange-to-flange industrial gas turbine analysis based on through-flow modelling[J]. International Journal of Gas Turbine, Propulsion and Power Systems, 2016, 8(1): 9-18.
[17] PETROVIC M V, WIEDERMANN A. Fully coupled through-flow method for industrial gas turbine analysis: GT2015-42111[R]. New York: ASME, 2015.
[18] 施发树, 刘兴洲. 多部件模型在全尺寸小型双函道涡扇发动机气流数值模拟中的应用[J]. 推进技术, 1998, 19(4): 22-26.
SHI F S, LIU X Z. Multicomponent models in application to numerical simulation of a small full-sized by-pass turbofan engine[J]. Journal of Propulsion Technology, 1998, 19(4): 22-26. (in Chinese)
Cited By in Cnki (14) | Click to display the text
[19] 黄家骅, 于廷臣, 冯国泰. 某小型涡扇发动机全流道准三维数值解法[J]. 航空发动机, 2005, 31(2): 42-45.
HUANG J H, YU T C, FENG G T. Quasi-3D numerical method for a small turbofan engine flow passage[J]. Aeroengine, 2005, 31(2): 42-45. (in Chinese)
Cited By in Cnki (6) | Click to display the text
[20] 冯国泰, 黄家骅, 王松涛. 航空发动机数值仿真试验台建立中几个关键技术问题的讨论[J]. 航空动力学报, 2002, 17(4): 483-488.
FENG G T, HUANG J H, WANG S T. Discussion on key technical problems in developing numerical simulation test-bed for aeroengine[J]. Journal of Aerospace Power, 2002, 17(4): 483-488. (in Chinese)
Cited By in Cnki (25) | Click to display the text
[21] 曹志鹏, 刘大响, 桂幸民, 等. 某小型涡喷发动机二维数值仿真[J]. 航空动力学报, 2009, 24(2): 439-444.
CAO Z P, LIU D X, GUI X M, et al. Two dimensional numerical simulation of small turbojet engine[J]. Journal of Aerospace Power, 2009, 24(2): 439-444. (in Chinese)
Cited By in Cnki (2) | Click to display the text
[22] 万科, 朱芳, 金东海, 等. 周向平均方法在某风扇/增压级分析中的应用[J]. 航空学报, 2014, 35(1): 132-140.
WAN K, ZHU F, JIN D H, et al. Application of circumferentially averaged method in fan/booster[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(1): 132-140. (in Chinese)
Cited By in Cnki (2) | Click to display the text
[23] BLAZEK J. Computational fluid dynamics:Principles and applications [M]. 2nd ed. Amsterdam: Elsevier, 2005: 29-31.
[24] KOCH C C, SMITH L H. Loss sources and magnitudes in axial-flow compressors[J]. Journal of Engineering for Gas Turbines and Power, 1976, 98(3): 411-424.
Click to display the text
[25] CREVELING H, CARMODY R. Axial flow compressor computer program for calculating off design performance: NASA-CR-72427[R]. Washington, D.C.: NASA, 1968.
[26] HOWELL A R. Fluid dynamics of axial compressors[J]. Proceedings of the Institution of Mechanical Engineers, 1945, 153(1): 441-452.
Click to display the text
[27] AUNGIER R, FAROKHI S. Axial-flow compressors:A strategy for aerodynamic design and analysis[J]. Applied Mechanics Reviews, 2004, 57(4).
Click to display the text
[28] GALVAS M R. Analytical correlation of centrifugal compressor design geometry for maximum efficiency with specific speed: NASA-TN-D-6729[R]. Washington, D.C.: NASA, 1972.
[29] KACKER S C, OKAPUU U. A mean line prediction method for axial flow turbine efficiency[J]. Journal of Engineering for Gas Turbines and Power, 1982, 104(1): 111-119.
Click to display the text
[30] CARTER A, HUGHES H. A theoretical investigation into the effect of profile shape on the performance of aerofoils in cascade: No. 2384[R]. London: British Aeronautical Research Council, 1946.
http://dx.doi.org/10.7527/S1000-6893.2019.23199
中国航空学会和北京航空航天大学主办。
0

文章信息

刘晓恒, 周成华, 宋满祥, 金东海, 桂幸民
LIU Xiaoheng, ZHOU Chenghua, SONG Manxiang, JIN Donghai, GUI Xingmin
基于通流方法的某涡喷发动机整机数值仿真
Overall simulation of a turbojet engine based on throughflow method
航空学报, 2020, 41(1): 123199.
Acta Aeronautica et Astronautica Sinica, 2020, 41(1): 123199.
http://dx.doi.org/10.7527/S1000-6893.2019.23199

文章历史

收稿日期: 2019-06-06
退修日期: 2019-07-03
录用日期: 2019-09-26
网络出版时间: 2019-10-12 11:21

相关文章

工作空间