姜 波，王继强，何小龙，胡忠志

（南京航空航天大学江苏省航空动力系统重点实验室，南京210016）

0 引言

航空发动机由静止逐渐加速至慢车状态这一过程称为起动过程。这个过程发动机气动热力参数变换剧烈，如何保证这个过程成功、快速的完成对延长发动机寿命、保证飞机安全有着重要的意义。对该过程进行数学抽象，建立起动过程的数值计算模型不仅能够为设计人员提供许多有价值的信息，而且能够让设计人员在发动机设计的初期就展开合理的起动过程控制规律优化研究。因而建立能够准确反映发动机起动过程的数值仿真模型对发动机的理论研究及工程应用都有着极大地意义。

在20世纪80年代，国外已经展开了发动机慢车以下状态的气动热力学建模研究，M.Yunis和R.K.A-grawal在1980年基于流量平衡的假设，通过外推方法获得低转速部件特性从而建立起了单转子发动机的起动模型[1]。在90年代NASA的Chappell采用慢车以上建模的方法，建立了包含高空起动、风车起动和地面起动的AVEST-V3模型[2]。此外在2001年和2007年，英国弗吉尼亚理工大学[3]英国克兰菲尔德大学[4]的学者在上述研究的基础上针对部件特性外推公式及零转速附近特性修正展开了深入的研究。国内学者对此研究展开的时间较晚，但也取得了大量的研究成果，如：基于大量发动机试车数据，使用系统辨识的方法建立的实时起动模型[5-8]；基于容积动力学和部件特性外推的方法建立双轴混合排气的涡扇发动机的地面及空中起动模型[9-10]，利用压气机及涡轮特性数据的脊柱特性建立能够准确模拟单轴发动机风车起动的部件级模型[11]等。但这些研究大都基于大量的试验数据并且在起动初始点需要根据经验反复的试凑出十分准确的初始值，工作量大且过程较为复杂。

本文利用NASA公开的PW公司JT9D大涵道比双转子涡扇发动机慢车以上部件级模型，通过部件特性外推、起动附件建模、修正部件总压恢复系数和改进模型求解算法建立了JT9D发动机起动至最大转速整个过程的仿真模型，并通过仿真结果验证了模型的准确性及改进求解算法的有效性。

1 JT9D发动机部件级模型

建立的起动过程模型是基于美国PW公司于20世纪60年代年设计制造的大涵道比双转子涡扇发动机JT9D的慢车以上状态的部件级模型，具体结构如图1所示，具体的建模过程这里不再赘述，参照文献[12]。

图1 JT9D发动机部件结构

2 双轴涡扇发动机起动过程

双轴涡扇发动机在起动时没有空气流经整台发动机，因此燃烧室无法点火。只有当转子达到一定的转速后，燃烧室内的气流才能建立稳定燃烧所需的气流压力和温度。因而发动机地面起动过程必须依靠起动机带转，一般情况下，双轴涡扇发动机的起动过程（如图2所示）分为以下3个阶段。

图2 双轴涡扇发动机起动过程扭矩特性

第1阶段：燃烧室未供油，起动机带动高压转子，起动机扭矩MST大于起动过程阻力矩MC/ηm，高压转子转速逐渐增加到n1。

第2阶段：燃烧室开始供油并点火，涡轮开始产生力矩，高压转子转速逐渐增加到n2。

第3阶段：起动机脱开，涡轮产生足够的力矩并将高压转子转速独立带转至n3。

3 起动建模关键问题及技术

3.1 低转速部件特性获取

发动机的部件特性在高转速时差别较大，但在低转速范围内却有着很大的相似性[1]。因而，采用部件特性指数外推的方法，能够在缺少试验特性的情况下得到较为合理的低转速部件特性。采用的指数外推公式如下

式中：p、q、r为外推指数；cor表示换算状态；ref表示参考状态；Wc为换算流量；N为转速；PR为增压比（涡轮特性为膨胀比ER）；η为效率。以风扇为例，选取换算转速为1的特性线为参考线。利用原有的特性数据通过不断变换外推指数的数值，使得不同换算转速线上最高效率状态点连成的线近似为线性（效率仅在低转速Ncor＜0.8时成线性）进而确定外推指数。不同外推指数下风扇流量、压比、效率、随换算转速变化分别如图3～5所示，不同外推指数下的外推结果对比如图6所示。

图3 不同外推指数下风扇流量随换算转速变化

图4 不同外推指数下风扇压比随换算转速变化

图5 不同外推指数下风扇效率随换算转速变化

在确定了外推指数后，以换算转速为0.6的试验特性线为参考，外推换算转速为0.5的特性线并与原有试验特性对比，结果显示：按照不同换算转速线上的最高效率点连成的线近似为线性，所选取的外推指数能更好的反映低转速的部件特性。据此确定外推指数见表1。

图6 不同外推指数下的外推结果对比

表1 部件特性外推指数

在确定部件特性的外推指数后，以各部件试验的最低状态换算转速线为外推参考线，外推各部件的风扇换算流量-压比特性、风扇换算流量-效率特性、高压涡轮换算流量-膨胀比、高压涡轮换算流量效率特性特性分别如图7～10所示。

图7 风扇换算流量-压比特性

3.2 起动机建模

JT9D发动机使用美国汉密尔顿公司提供的PS-700空气涡轮起动机[13]，主要参数见表2。

空气涡轮起动机在发动机起动过程的输出扭矩与高压转子转速近似为线性关系[14]。起动机扭矩特性如图11所示。

图8 风扇换算流量-效率特性

图9 高压涡轮换算流量-膨胀比特性

表2 PS-700空气涡轮起动机主要参数

根据扭矩特性按照如下公式建立起动机模型

式中：P为起动机输出功率；n2为起动机脱开时高压转子转速；N为高压转子转速；T为起动机输出扭矩。

图11 起动机扭矩特性

3.3 改进牛顿法及其应用

在求解JT9D发动机共同工作非线性方程组时NASA提供了1种传统的Newton-Raphson迭代求解方法。此方法在求解发动机慢车以上状态共同工作非线性方程组时能够在选定较为准确的迭代初始值情况下达到快速收敛的效果[15-17]。但发动机起动过程部件特性的相对值变化剧烈，即使在给出相对准确的迭代初始值情况下也需要多次求解才能达到收敛精度的要求。并且每次迭代求解共同工作非线性方程组时计算机都需要对模型进行气动热力计算，然后再对平衡方程偏导数组成的Jacobian矩阵进行中心差分，庞大的计算量降低了模型的执行效率，使得计算时间变长，模型的实时性降低。

本文所述的改进牛顿法是参考广泛用于国外先进数值仿真程序，如美国航空航天局NASA的NPSS（Numerical Propulsion Simulation System）[18]、荷兰国家航空航天实验室（NLR）的GSP（Gas turbine Simulation Program）[19]等在求解非线性方程组时使用的Broyden法[20]：即在传统Newton-Raphson方法基础上改变Jacobian矩阵的迭代更新方式，进而减少模型的气动热力计算次数，提高模型的仿真求解速度。计算方法如下：

对于发动机的共同工作非线性方程组

式中：x=（x1,x2,…,xn）T，n 为平衡方程数，xi为发动机的某参数；εi为平衡方程的误差量；fi（x）代表不同的平衡方程，且F（x）=（f1（x1）,f2（x2）,…,fn（xn））T，Jacobian矩阵为

其中：偏导数项在模型求解计算时需根据中心差分求得。根据Newton-Raphson求解方法

式中：k为迭代次数，此时构造如下矩阵

根据Broyden法求解原理[20]，Jacobian矩阵迭代更新的修正项为

使用Sherman-Morrison公式解析的求Jacobian矩阵的逆，从而可得Jacobian矩阵在更新时

当 |F（xk）|＜δ（设定误差允许量）时，认为迭代收敛。

这种改进的牛顿法在求解发动机共同工作非线性方程组时第1步仍采用中心差分的方法得到Jacobian矩阵，而后续的迭代更新只需通过计算ΔJk就可完成。这样模型的动态仿真过程就避免了每次通过中心差分方法重新计算Jacobian矩阵而带来的庞大计算量，从而大大提高了需要多次迭代求解发动机共同工作方程的起动过程模型的执行效率。

3.4 冷转动过程模型

发动机在点火前的工作状态即起动过程的第一阶段称为冷转动过程。这一过程涡轮和压气机都是负载，建模方法如下：

首先计算起动机和压气机的扭矩，然后根据转子动力学计算高压转子转速

式中：ΔT为起动机扭矩与负载扭矩差；J为高压转子转动惯量；N为高压转子转速；t为时间。

之后，确定1个合适的点火状态，根据此时的转速 nR，WR和压比 PR利用式（13）、（14）计算对应时刻下的参数

最后，低压转子转速及其他截面状态参数按高压转子转速进行线性插值即可得到。

3.5 部件总压恢复系数修正

NASA公开的JT9D发动机模型各部件的总压恢复系数都取了定值，这在慢车以上状态建模时是适用的。但在发动机的起动过程中由于各部件的流量、压比、转速都很小，真实的总压恢复系数更接近于1，并且随着转速的增大而逐渐减小[9]。因而在建立起动模型时，各部件的总压恢复系数在点火前取近似略小于1的常数，而在点火状态和慢车状态之间时按照转速进行线性插值

式中：a为点火状态；b为慢车状态；σ为总压恢复系数。

4 仿真结果分析

NASA公开的JT9D发动机模型的仿真范围为：70%NHmax～100%NHmax。按照上述建模方法，补充建立包含起动过程的全状态性能模型，并在地面标准状况下给定供油规律（如图12所示），仿真结果分别如图13～15所示。

图12 供油曲线

图13 转速变化曲线

图14 推力变化曲线

图15 燃烧室出口温度变化曲线

由仿真结果可以看出，起动过程的第1阶段燃烧室没有供油。高压转子由起动机带动，转速缓慢上升并带动低压转子转速上升，此时推力也缓慢上升，燃烧室出口温度基本不变。在第16 s时开始供油，转子转速和燃烧室出口温度急剧上升至慢车状态。在第40 s以后增加供油，发动机由慢车状态逐渐运行到最大转速状态，转子转速、推力及燃烧室出口温度都增加至100%。整个模型能很好的反映起动过程的各阶段并且与原有慢车以上状态模型有效地连接，验证了建模方法的有效性。

为了进一步验证建模方法的准确性，将模型的仿真结果与仅有的ICAO（国际民航组织）公布的JT9D-7J发动机数据进行对比，结果见表3。

表3 7%推力状态下（52%NHmax）模型仿真结果验证

各部件总压恢复系数修正前（修正前采用固定值σb）与修正后高压转子转速的仿真结果对比如图16所示，燃烧室点火至慢车状态过程的高压转子转速最大相对误差为1%。

图16 总压恢复系数修正前后仿真结果对比

同时，模型在点火后由各部件特性迭代求解的初始点计算时间最长，最难收敛。因而在此处用改进牛顿法进行验证效果最为明显。

单个步长Jacobian矩阵计算时间见表4。表4计算的硬件环境为：AMD A6-3400M（CPU），主频1.4 G Hz，内存4 GB。由结果可以看出改进牛顿法的求解速度要明显高于传统的Newton-Raphson方法。并且，基于改进牛顿法的原理，当发动机共同工作方程和独立变量个数增加以及迭代初始值偏离真值程度越大时，这种效果将会变得更加明显。

表4 单个步长的Jacobian矩阵计算时间

5 结论

通过对JT9D发动机部件级起动过程建模、使用改进牛顿法替换传统的Newton-Raphson方法仿真求解，得到以下结论：

（1）采用部件特性指数外推法并使用不同换算转速线上最高效率状态点连成的线，近似为线性来确定外推指数，可以得到合理的各部件低转速部件特性；

（2）起动附件建模及转子冷转动过程建模能够有效地帮助建立完整的起动过程模型。

（3）发动机起动过程各部件总压恢复系数按照转速变化修正能够有效的提高模型的精度。

（4）采用改进牛顿法替代传统的Newton-Raphson方法可以改善模型的执行效率，大幅提高模型的仿真求解速度。

[1]Agrawal R K,Yunis M.A generalized mathematical model to estimate gas turbine starting characteristics[J].Physical Review B,1981（1）:5668-5687.

[2]Chappell M A,Mclaughlin P W.Approach of modeling continuous turbine engine operation from startup to shut down[J].Journal of Propulsion&Power,2015,9（3）:466-471.

[3]Sexton W R.A new method to control turbofan engine starting by varying compressor surge valve bleed[D].Virginia:Virginia Polytechnic Institute and State University,2001.

[4]Howard,Jason.Sub-idle modeling of gas turbines:altitude relight and wind milling.[D].Cranfield:Cranfield University,2007.

[5]于达仁,郭钰锋,牛军,等.涡喷发动机风车启动工况的神经网络建模[J].推进技术,2001,22（3）:183-186.YU Daren,GUO Yufeng,NIU Jun,et al.Turbojet modeling in wind milling based on radial basis function networks[J].Journal of Propulsion Technology,2001,22（3）:183-186.（in Chinese）

[6]卓刚.航空发动机智能建模与故障诊断研究[D].南京：南京航空航天大学,2004.ZHUO Gang.The research on aeroengine intelligent modeling and fault diagnosis[D].Nanjing:Nanjing University of Aeronautics and Astronautics.（ in Chinese）

[7]刘建勋,王剑影,李应红,等.某型发动机起动模型的支持向量机辨识及应用[J].推进技术,2004,25（5）:401-404.LIU Jianxun,WANG Jianying,LI Yinghong,et al.An identification model of aero engine starting based on support vector machine and its application[J].Journal of Propulsion Technology，2004，25（5）:401-404.（in Chinese）

[8]李永进,张海波,张天宏.一种考虑非线性余项的机载发动机自适应模型建立及其在寻优控制中的应用 [J].推进技术,2016,37（1）:172-180.LI Yongjin，ZHANG Haibo，ZHANG Tianhong.Establishment and application in performance seeking control of an on-board adaptive aero-engine model considering nonlinear remainders[J].Journal of Propulsion Technology，2016，37（1）:172-180.（in Chinese）

[9]周文祥,黄金泉,窦建平.涡扇发动机部件级起动模型[J].航空动力学报,2006,21（2）:248-253.ZHOU Wenxiang,HUANG Jinquan,DOU Jianping.Development of component-level startup model for a turbofan engine[J].Journal of Aerospace Power.2006,21（2）:248-253.（in Chinese）

[10]周文祥,黄金泉.涡扇发动机高空起动模型研究[J].航空动力学报,2007,22（8）:1384-1390.ZHOU Wenxiang,HUANG Jinquan.Research on the startup model of turbofan engine at high altitude[J].Journal of Aerospace Power,2007,22（8）:1384-1390.（in Chinese）

[11]施洋,屠秋野,严红明,等.一种航空发动机全状态性能模型[J].航空动力学报,2017,32（2）:373-381.SHI Yang,TU Qiuye,YAN Hongming,et al.A full states performance model for aero engine[J].Journal of Aerospace Power,2007,22（8）:1384-1390.（in Chinese）

[12]AIAA.A process for the creation of T-MATS propulsion system models from NPSS data[C]//.AIAA Joint Propulsion Conference.Cleveland,OH.2014.

[13]胡晓煜.世界中小型航空发动机手册[M].北京.航空工业出版社,2006:443-445.HU Xiaoyu.World minitype aeroengine manual[M]Beijing:Aviation Industry Press,2006:443-445.（ in Chinese）

[14]杨帆,樊丁,彭凯,等.基于试车数据的航空发动机起动过程建模[J].空军工程大学学报:自然科学版,2013（6）:1-4.YANG Fan,FAN Ding,PENG Kai,et al.Modeling of an aero-engine in starting process based on the test data[J].Journal of Air Force Engineering University:Natural Science Edition,2013（6）:1-4.（in Chinese）

[15]Sellers J F,Daniele C J.Dyngen:a program for calculating steady-state and transient performance of turbojet and turbofan engines[R].NASA-TN－D-7901,1975.

[16]Daniele C J,Krosel S M,Szuch J R,et al.Digital computer program for generating dynamic turbofan engine models [R].NASA-TM-83446,1983.

[17]李庆扬,莫孜中,祁力群.非线性方程组的数值解法[M].北京：科学出版社,1999:108-113.LI Qingyang,MO Zizhong,QI Liqun.Numerical solutions of nonlinear equations[M].Beijing：Science Press,1999:108-113.（in Chinese）

[18]Claus R W,Evans A L,Lylte J K,et al.Numerical propulsion system simulation [J].Computing Systems in Engineering,1991,2（4）:357-364.

[19]Visser W P J,Broomhead M J.GSP,a generic object-oriented gas turbine simulation environment[C]//ASME Turbo Expo 2000:Power for Land,Sea,and Air.2000:V001T01A002.

[20]Broyden C G.Quasi-Newton methods and their application to function minimization[J].Mathematics of Computation,1967,21（99）:368-381.