2.2 航空发动机非线性部件级模型建立
本节以涡轴发动机为例,简要介绍航空发动机部件级建模方法。
涡轴发动机部件级数学模型必须能够反映出真实发动机的工作特性和工作原理[44],[55-57]。具体表现为:
(1)进气道气体的压力和温度;
(2)从进气道进来的气体经过压气机后的温度变化以及增压比;
(3)燃烧室内所进行的化学反应引起的燃烧气体焓值、气体压力、气体流量、油气比以及温度值变化;
(4)涡轮前燃气温度,经过燃气涡轮后燃气的焓值、温度、压力以及油气比变化;
(5)燃气经过动力涡轮后气体的焓值、温度、压力以及油气比变化;
(6)尾喷管内气体静压变化以及扩压损失。
部件级模型的缺点:
(1)部件级模型依赖于部件特性,因此部件特性的准确度会影响模型精度;
(2)在建模的过程中会对一些问题做适当简化或者假设,简化的程度也会在一定程度上影响部件级模型的精度;
(3)采用Newton-Raphson方法来求解非线性程组的时候,解的初猜值的选取也会影响模型的收敛。
由于航空发动机构造非常复杂而且发动机在运转的过程中受到很多内部、外部因素的影响,所以在建模之前,需要对发动机工作过程中的一些问题在不影响发动机的基本性能的前提下做一些适当的简化,从而减少建模的工作量,把工作重心放在基本模型的建立上。
本书的假设如下:
(1)忽略热惯性和部件通道容积动力学以及摩擦因素的影响;
(2)航空发动机中的气体流动使用准一元流模型描述;
(3)转子的转动惯量为定常数;
(4)在发动机中气流速度比声速小得多的部件间,气流的绝热指数和密度按滞止气流计算。
考虑如下涡轴发动机,其结构如图2-1所示。
图2-1 涡轴发动机结构[44]
截面:0—远方未扰动截面;1—进气道进口截面;2—进气道出口截面(压气机进口截面);3—压气机出口截面(燃烧室进口截面);4—燃烧室出口截面(燃气涡轮进口截面);
5—燃气涡轮出口截面(动力涡轮进口截面);6—动力涡轮出口截面(尾喷管进口截面);7—尾喷管出口截面。
部件:Ⅰ—进气道;Ⅱ—压气机;Ⅲ—燃烧室;Ⅳ—燃气涡轮;Ⅴ—动力涡轮;Ⅵ—尾喷管。
本节按照气体在涡轴发动机内部流动状态建立发动机的各部件气动热力模型。涡轴发动机的气体流动顺序:空气依次经过进气道、压气机、燃烧室、燃气涡轮、动力涡轮,最后经尾喷管排出。气体流动示意图见图2-2。
图2-2 涡轴发动机内部气体流动情况
本节采用S-Function的C语言模式,根据S-Function语法规则建立了涡轴发动机各个部件的模型。对压气机模型、燃烧室模型、燃气涡轮模型、动力涡轮模型、尾喷管模型、旋翼模型分别进行函数级封装。另外建立了初猜值模块、特性数据插值模块、系数修正模块、燃气热力学计算模块、Newton-Raphson方法求解这些非线性方程组模块等通用模块。
2.2.1 涡轴发动机各个部件模型
2.2.1.1 进气道模型
进气道是喷气发动机所需空气的进口和通道,进气道在负责给发动机提供一定流量的空气的同时,还需要保证进气道出口也就是压气机进口的气体的流场能够保证压气机以及燃烧室等部件的正常工作。其中在超声速飞行器中进气道还需要对高速气体进行加速增压,其增压作用在超声速下甚至超过压气机。
已知海平面附近的大气温度Ts=288.15K,海平面标准大气压强 ps=1.01225×105Pa,11km高度处的大气压强为p11=2.26×104Pa,设飞机的飞行高度为H,则:
(1)当高度H<11km时,当地的静温T0、大气静压p0是:
(2)当高度H>11km时,当地的静温T0、大气静压p0是:
涡轴发动机主要配备在直升机上,设直升机飞行的马赫数为Ma,空气绝热指数k=1.4,则可以得到发动机进气道进口截面的气体总温T1和总压p1的计算方法:
直升机一般飞行高度为0~6km,飞行马赫数为0~0.4,因而在整个飞行包线内,压气机进口所需的气体的流速均大于飞机的前飞速度,压气机一般为收敛型,收敛型的压气机会引起流经气体的总压损失。一般在收敛型进气道的设计中,气体的总压恢复系数σ0是衡量进气道性能的最为重要的特性参数之一,由于直升机的飞行马赫数比较小,所以在飞行包线内总压恢复系数的摄动也不大,为了减少建模复杂度,我们一般认为总压恢复系数σ0=0.998,为一个恒定值。在上节的假设中我们假设进气道模型中气流在进气道内是个绝热过程,因此在压气机出口气体总温T2与压气机入口气体总温T1相同,压气机出口气体总压p2只需要考虑总压损失系数。所以压气机出口气体总温T2和总压p2参数为
T2=T1
p2=p1σ0
2.2.1.2 压气机模型
压气机在航空发动机中利用高速旋转压气机叶片做功,提高压气机出口的气体总压,使得燃烧室内的混合燃气的总压达到有利于燃烧的状态。压气机进口气流流场状态对压气机的性能有着很大的影响,为了减少建模的复杂度,我们假设压气机进口流场是稳定、均匀的。压气机的设计点一般给出压气机运转的转速、空气流量、气体增压比。但是压气机在实际的工作过程中,一般都是偏离设计点的,同时它的工作状态是不断变化的,压气机在各种工作状态下的性能数据称为压气机特性。在发动机的各部件特性中,压气机特性对发动机的性能影响最大。压气机特性曲线即压比πC、效率ηC、进口相似流量WC2、相对转速NCcor以及压比系数ZC之间的换算关系曲线,根据压气机特性图上的数据进行插值,可得到压气机相关的压比、效率和流量。
(1)根据压气机百分比转速NC和压气机进口总温T2的初始值,可计算压气机相对转速:
其中,NCds为压气机设计点百分比转速,T2ds为压气机设计点进口温度。
本节采用的涡轴发动机模型的压气机设计点百分比转速为NCd s=100,压气机设计点进口温度T2ds=288.15K。
(2)根据实验测得的压气机特性数据,由压气机设计点百分比转速NCds和初猜压比系数ZC通过对压气机特性数据进行插值计算得到压气机压比πC、压气机效率ηC和压气机空气流量WC2。
(3)当发动机状态发生变化时,发动机的工作状态会偏离设计点,因此需要对非设计点的特性参数进行修正,以保证部件级模型能够模拟真实发动机的工作状态。本节采用的涡轴发动机模型的设计点为:NC=100, NCcords=1.0, πCds=10.4, ηCds=0.798, WC2ds=4.44kg/s(其中NCcords为压气机设计点相对转速,πCds为压气机设计点压比,ηCds为压气机设计点效率,WC2ds为压气机设计点进口相似流量)。相应的压比、效率、流量等插值数据的修正系数计算方法分别为
其中,Cfπ、Cfη、CfW分别为压比、效率、进口相似流量的修正系数。
(4)使用以上三式中的3个修正系数Cfπ、Cfη、CfW修正压气机在偏离设计点时的特性数据,可得出发动机实际工作中的压比、效率和进口相似流量的计算方法,如下:
πaC=(πC-1)Cfπ+1
ηaC=ηCCfη
WaC2=WC2CfW
其中,πaC、ηaC、WaC2分别为压气机实际压比、效率和进口相似流量。
(5)根据前述的简化和假设,可以认为气体在压气机内部的流动状态是一个等熵过程,所以压气机出口气体热力学参数计算方法如下:
压气机进口气体熵S2、理想出口气体熵S3I:
S2=f(T2)
S3I=S2+lnπC
其中,f表示部件级模型中的热力学求解函数。
压气机理想出口气体焓H3I、实际出口焓H3:
H3I=f(S3I)
H3=H2+(H3I-H2)/ηC
其中,H2为压气机进口焓,f表示部件级模型中的热力学求解函数。
压气机实际出口总温、总压、进出口流量:
T3=f(H3)
p3=p2πC
其中,T3、p3、WaC3、WaC2分别为压气机实际出口总温、出口总压、出口相似流量、进口相似流量。f表示部件级模型中的热力学求解函数。
压气机消耗功HPC:
HPC=WaC2(H3-H2)
2.2.1.3 燃烧室模型
燃烧室是燃料在其中燃烧生成高温高压燃气的部件,其主要作用是把压气机加压后的空气与燃油混合后进行点火燃烧,并保证燃烧的稳定性。它把燃油的化学能转化为燃气的内能,燃气膨胀通过高低压涡轮做功来带动压气机和输出轴功率。燃烧室模块主要参数为:燃烧室进口气体流量WaC3、供油量WFB、进口气体总温T3、总压 p3等,燃烧室设计点WFB=0.10898kg/s、燃油热值Hu=42900000J/kg、燃烧室燃烧效率ηb=0.98、总压恢复系数σb=0.965。
燃烧室出口流量Wg4、出口气体焓H4、总温T4、总压p4:
Wg4=WFB+WaC3
H4=(WaC3H3+WFB×Huηb)/Wg4
T4=f(H4)
p4=p3σb
2.2.1.4 燃气涡轮模型
从燃烧室出来的高温高压燃气需要在燃气涡轮中膨胀,推动燃气涡轮叶片做功,进而带动压气机以及其他部件旋转。燃气涡轮是将燃气的内能转化为机械能的装置,经过燃气涡轮后燃气的总温、总压等都将降低。根据假设从燃烧室出口到燃气涡轮进口,不计摩擦,我们认为燃气的总压、总温等没有损失。
(1)燃气涡轮进口气体参数
燃气涡轮进口总压p41:
p41=p4
来自压气机的第一路引气Wc321、Wc322和燃气Wg4混合后,构成燃气涡轮进口流量Wg41:
Wg41=Wg4+Wc321+Wc322
燃气涡轮进口焓值H41、总温T41、熵值S41:
H41=[Wg4 H4+(Wc321+Wc322)H3 ]/Wg41
T41=f(H41)
S41=f(H41, T41)
其中,f表示部件级模型中的热力学求解函数。
(2)燃气涡轮特性数据插值计算
燃气涡轮特性数据是指燃气涡轮的效率ηg、流经涡轮的燃气流量Wg 41、涡轮前后燃气总压的落压比系数πg以及燃气涡轮转子相对换算转速Ngcor之间的换算关系。与压气机类似,对于燃气涡轮特性需要对偏离设计点的工作状态进行参数修正。已知燃气涡轮的设计点参数为设计点转子转速百分比 Ngds=100,设计点转子相对转速Ngcords=1.0,设计点燃气落压比πgds=4.69,设计点燃气流量Wg4ds=4.2077kg/s,设计点燃气涡轮效率ηgds=0.85211,设计点燃气总温Tg4ds=1533K。将这些参数通过燃气涡轮特性数据进行插值后得到工作点落压比πg、效率ηg、燃气流量Wg 41。
根据初猜百分比转速Ng 和燃气涡轮进口总温T41计算燃气涡轮相对转速Ngcor方法为
根据实验测得的燃气涡轮特性数据,由设计点转子转速百分比Ngds和初猜压比系数Zg通过对燃气涡轮特性数据进行插值计算得到燃气涡轮落压比πg、压气机效率ηg 和压气机空气流量Wg41c。
(3)燃气涡轮非设计点数据修正
燃气涡轮落压比、效率、燃气流量的修正系数Cfgπ、Cfgη、CfgW的计算方法如下:
使用以上三式中的修正系数Cfgπ、Cfgη、CfgW来修正非设计的特性参数,可计算出真实发动机工作状态下的压比、效率和相似流量等特性,计算公式如下:
πag=Cfgπ(πg-1)+1
ηag=Cfgηηg
Wag 41c=CfgWWg41c
其中,πag、ηag、Wag 41c分别为燃气涡轮的实际压比、效率和相似流量。
(4)燃气涡轮出口热力学参数计算方法
燃气涡轮出口理想熵值S43I、理想焓值H43I:
S43I=S41-lnπg
H43I=f(H41)
其中,f表示部件级模型中的热力学求解函数。
考虑引气作用:
Wg43=Wg41+Wc241+Wc242
其中,Wc241、Wc242分别为燃气涡轮的第一路引气和第二路引气。
燃气涡轮出口实际焓值H43:
油气比FAR43、燃气涡轮出口截面总温T43、总压p43:
FAR43=WFB/(Wg43-WFB)
T43=f(H43, FAR43)
p43=p41/πg
燃气涡轮输出功HPg:
HPg=Wg41(H41-H43)-Pext
其中,Pext为抽取功。
2.2.1.5 动力涡轮模型
涡轴发动机动力涡轮主要用来输出轴功率,轴功率经涡轮动力轴输出至直升机旋翼、尾翼等所有动力负荷,动力涡轮轴与压气机轴在物理上没有关系,能够独立旋转,在本节中假设动力涡轮轴与压气机轴之间没有相互关系。从燃气涡轮出来的燃气能够在动力涡轴中继续膨胀,对动力涡轴做功,这也是动力涡轮能够输出轴功率的根本原因。燃气在动力涡轮中继续膨胀做功,将燃气剩余的大部分内能转化为动力涡轮可以输出的机械能。
(1)动力涡轮进口相关热力学参数计算方法
燃气涡轮出口到动力涡轮进口具有一定的压力损失,总压恢复系数σp=0.945,动力涡轮进口总压p44、来自压气机的引气Wc221和燃气Wg43混合后的动力涡轮进口流量Wg44、动力涡轮进口燃气焓值H44、动力涡轮进口总温T44、动力涡轮进口熵值S44分别为
p44=p43σp
Wg44=Wg43+Wc221
H44=(Wg43 H43+Wc221 H2)/Wg44
T44=f(H44)
S44=f(H44, T44)
(2)动力涡轮特性数据插值计算
动力涡轮的特性数据是指动力涡轮的效率ηp、流经动力涡轮的燃气流量Wp4、动力涡轮前后燃气的落压比πp,以及动力涡轮相对换算转速Npcor之间的换算关系。同样的,对于动力涡轮非设计点工作状态的特性数据也需要进行修正。已知动力涡轮的设计点数据为转子转速百分比Npd s=100,设计点相对转速Npcords=1.0,设计点落压比πpds=3.7,设计点燃气流量Wp4ds=4.7 kg/s,设计点燃气涡轮效率ηp d s=0.865。将这些参数在燃气涡轮特性数据上进行插值得到πp、ηp、Wp 44。
根据动力涡轮百分比转速Np 和动力涡轮进口总温T43计算动力涡轮相对转速Npcor的方法为
注:燃气涡轮出口截面总温与动力涡轮进口总温相等,因此统一用T43表示。
根据实验测得的动力涡轮特性数据,由Npds和初猜压比系数Zp并通过对动力涡轮特性数据进行插值计算得到动力涡轮落压比πp、涡轮效率ηp和涡轮空气流量Wg44c。
(3)动力涡轮非设计点数据修正
动力涡轮落压比修正系数Cfpπ、效率修正系数Cfpη、燃气流量修正系数CfpW分别为
使用以上三式中的修正参数Cfpπ、Cfpη、CfpW可以计算出发动机实际工作状态下的压比πap、效率ηap和相似流量Wap 44c等特性参数,计算公式为
πap=Cfpπ(πp-1)+1
ηap=Cfpηηp
Wap 44c=CfpWWp44c
(4)动力涡轮出口热力学参数计算方法
燃气涡轮出口理想熵值S45I、理想焓值H45I:
S45I=S44-lnπp
H45I=f(H44)
其中,f表示部件级模型中的热力学求解函数。
考虑引气作用:
Wg45=Wg44+Wc222
其中,Wc222为动力涡轮的引气。
动力涡轮出口燃气混合后焓值H5:
动力涡轮出口燃气油气比FAR5、燃气总温T5、燃气总压p5:
FAR5=WFB/(Wg5-WFB)
T5=f(H5, FAR5)
p5=p44/πp
动力涡轮的输出轴功HPp:
HPp=Wg44(H44-H45)
2.2.1.6 尾喷管模型
从动力涡轮出口出来的燃气的静压仍然低于当地大气气压,所以尾喷管的主要作用就是使燃气能够继续膨胀做功,提高燃气静压,降低燃气流速,将燃气的可用功尽可能地转化为机械能推动发动机做功。同时高温燃气经过尾喷管后燃气的温度大大下降,排出体外能增加其隐身性能。涡轴发动机的尾喷管通常行程比较长,部分设计中尾喷管能够向下弯曲,在直升机起飞或飞行过程中能够提供部分推力。
(1)根据假设,燃气流过尾喷管过程可以看作是一个等熵绝热的过程。尾喷管进口截面与动力涡轮出口截面为同一个截面,因此可以计算出尾喷管进口燃气的热力学参数。
尾喷管进口总压p6、总温T6、燃气流量Wg6、油气比FAR6分别为
p6=p5
T6=T5
Wg6=Wg5
FAR6=FAR5
(2)压力损失系数σN是衡量尾喷管性能的重要指标,通过σN能够计算尾喷管出口总压p7:
p7=p6σN
2.2.1.7 旋翼模型
直升机的旋翼是直升机最主要的动力负载,也是直升机相对于别的飞行器的一个明显的特质,旋翼赋予了直升机垂直起降、空中悬停、离地数米的超低空超低速飞行等相对于固定翼飞行器特异的性能。旋翼的动力来自涡轴发动机输出的轴功率,直升机的飞行姿态大多数通过调整旋翼和尾翼来实现。直升机的升力靠调整旋翼的总距杆来实现,前后向、左右向等水平动力靠调整旋翼的周期变距杆来实现。
通过数值拟合的方法(简化处理)可以计算出旋翼所需的扭矩M:
M=C1[C2-pMain(C3-C4pMain)]×N2p
式中,pMain为旋翼负载,C1~C4为拟合系数。
旋翼负载功率HPairscrew的计算方法:
2.2.2 涡轴发动机稳态模型
2.2.2.1 稳态控制方程
航空发动机稳定平衡条件如下。
(1)压气机与涡轮功率平衡:
(2)负载与动力涡轮功率平衡:
(3)燃气涡轮进口流量连续:
(4)动力涡轮进口流量连续:
(5)尾喷口流量出口压力连续:
其中,PL 为负载功率,Wg41X为燃气涡轮进口流量猜值,Wg44X为动力涡轮进口流量猜值。
发动机达到稳态平衡要保证以上五个方程Eq0~Eq4都为0。
2.2.2.2 稳态模型求解
本节非线性方程组求解采用经典Newton-Raphson方法。上一节中所描述的稳态平衡控制方程是发动机达到稳态的必要条件,方程求解成功即表示模型的各性能指标参数达到稳态要求,因此需要先试取五个性能指标参数作为方程要求解的参数,这些试取的未知参数称为猜值参数,猜值参数选择是否恰当将关系到方程能否求解成功,或者模型达到稳态时的性能优劣。发动机各截面流量W、各转动部件功率PW皆是压气机相对转速NCrcs、动力涡轮相对转速NPrcs、燃气涡轮压比系数CGπ、动力涡轮压比系数CPπ以及压气机压比系数CCπ的非线性函数,具体的非线性关系由发动机各部件数学模型和特性决定。因此在这里选择NCrcs、NPrcs、CGπ、CPπ、CCπ为稳态平衡控制方程的猜值参数,式(2-1)至式(2-5)可记作如下形式:
通过稳态共同工作方程的约束和联系,在供油量给定时对发动机稳态工作点的确定就转化为解一组以NCrcs、NPrcs、CGπ、CPπ、CCπ为独立变量的非线性方程组的问题。本节采用经典的Newton-Raphson方法,以迭代方式求解非线性方程组的解,设置一定的精度,当方程都满足精度要求时,迭代完毕,表示方程求解成功。稳态计算流程如图2-3所示。
图2-3 部件级模型稳态计算流程
用Newton-Raphson方法修正一组初猜参数NCrcs、NPrcs、CGπ、CPπ、CCπ的值,使非线性方程组(2-6)在误差绝对值εmin<10-6意义下成立,即
fi(NCrcs, NPrcs, CGπ, CPπ, CCπ)=εi≤εmin i=1,2,3,4,5
用Newton-Raphson方法按偏导数方向修正猜值参数,设第k+1步初猜参数的值为
ni|k+1=nik+Δnii=1,2,3,4,5
n1=CCπ, n2=CGπ, n3=CPπ, n4=NCrcs, n5=NPrcs
其中:
A称为雅克比矩阵,具体表达式为
式中偏导数的计算按中心差分法求,即
2.2.3 涡轴发动机动态模型
2.2.3.1 动态控制方程
在发动机动态过程中,由于发动机工作处在非平衡状态,压气机与燃气涡轮功率以及旋翼负载与动力涡轮功率不平衡,但同时各截面流过的流量还满足连续条件,在动态平衡过程中,动力涡轮转子要保持恒定或者在微小范围波动以维持旋翼的转速恒定,可得动态过程中的三个平衡控制方程:
(1)燃气涡轮进口流量连续,同式(2-3)。
(2)动力涡轮进口流量连续,同式(2-4)。
(3)尾喷管出口压力平衡,同式(2-5)。
2.2.3.2 动态模型求解
发动机模型在动态过程中有三个平衡控制方程,因为动态过程中满足流量连续以及压力平衡,发动机各部件流量与转子功率皆可为压气机压比系数CCπ、燃气涡轮压比系数CGπ、动力涡轮压比系数CPπ所表示的非线性函数,因此选择CCπ、CGπ、CPπ为动态平衡控制方程的猜值参数。当发动机模型从稳态开始动态计算时,发动机模型根据马赫数Ma、飞行高度H、总距θ进行各部件模型计算,这里同样采用Newton-Raphson法求解非线性方程组,动态计算流程如图2-4所示。
图2-4 部件级模型动态计算流程
平衡条件可表示为
fi(CCπ, CGπ, CPπ)=0 i=1,2,3
非线性方程组在误差绝对值εmin<10-6意义下成立,即
fi(CCπ, CGπ, CPπ)=εi≤εmin i=1,2,3
满足上式条件可认为模型得到了在动态点的解。用Newton-Raphson法计算发动机模型得到第k+1步猜值参数为
ni|k+1=ni|k+Δnii=1,2,3
其中,n1=CCπ, n2=CGπ, n3=CPπ。式中Δni:
雅克比矩阵A为
偏导数同样采用中心差分法求取:
基于本节方法设计的涡轴发动机模型Matlab/Simulink混合编程模型模块图如图2-5所示,中间的矩形为采用Matlab/C++混合编程编写的发动机模型。
图2-5 Matlab/Simulink混合编程模型模块