第36卷第1期
计算机仿真2019年1月文章编号：1006-9348(2019)01-0076-06
航空发动机气动热力系统模型求解方法研究
伍谦，李秋红，单睿斌
(南京航空航天大学能源与动力学院，江苏南京210016)
摘要：航空发动机气动热力系统模型通过求解部件之间共同工作方程组，来确定发动机相关状态参数。传统的共同工作方
程组求解方法是采用定步长Newton迭代方法，每次迭代需要多次调用模型计算雅可比矩阵，为了保证模型的收敛性，通常
步长设置比较保守，在一个仿真周期内需要进行多次的迭代计算进行方程组求解，导致航空发动机气动热力系统模型的实
时性较差。提出了一种混合变步长Newton-逆Broyden方法求解共同工作方程组。为保证收敛性，在猜值初始偏差较大时
采用小步长的Newton迭代法求解方程。为保证收敛速度，在猜值接近真实解的时候采用较大步长的逆Broyden递推法求解
方程，并引入松弛因子提髙迭代稳定性。应用上述方法求解F404/YJ101变循环发动机的气动热力系统模型，进行了发动机
稳态和动态工作仿真。结果表明，混合Newton-逆Broyden方法在航空发动机气动热力系统模型计算中的收敛能力、实时性
和计算精度有着明显的优越性。
关键词：航空发动机；非线性方程组；拟牛顿法；松弛因子
中图分类号:V231文献标识码：B
Rearch on Methodfor Solving Aero-Engine Thermodynamic System Model
WU Qian,LI Qiu-hong,SHAN Rui-bin
(College of Energy and Power Engineering,Nanjing University of Aeronautics and Astronautics,Nanji
ng Jiangsu210016,China) ABSTRACT:The aero engine thermodynamic system model is ud to get the engine working parameters by solving
the components co-working equations.The traditional method for solving co-working equations is the fixed step New
ton method,each iteration of which needs to repeatedly call the model to calculate the Jacobian matrix,and the step
is conrvative to ensure the convergence,resulting in poor real-time.A hybrid variable step Newton-inver Broy
den method was propod in the paper.To ensure the convergence,the guess value deviation was greater when using
Newton method with small step,and the guess value deviation was clo to the roots when using inver Broyden
method with large step.To improve the method1s stability,relaxation factor was ud in solving process.The method
怎么包包子
was applied to solve the thermodynamic system model of the F404/YJ101variable cycle engine,and the steady and
dynamic simulation of the engine was realized.Simulation results show the superiority of the propod method on con
vergence,real-time performance and accuracy in the solving process of aero-engine thermodynamic model.
KEYWORDS:Aero-engine；Nonlinear equations；Quasi-newton；Relaxation factor
1引言
随着数值仿真技术的快速发展，航空发动机气动热力系统数值模型在型号研究早期有着越来越重要的辅助作用。同时，高实时性、高精度的航空发动机数值模型也是航空发动机智能控制、健康管理等先进技术的重要基础，因而数值仿真模型一直得到国内外的重视「3。
航空发动机部件级数值模型能够对发动机各个截面的气动热力参数和发动机性能参数进行仿真计算，在发动机性能仿真与控制系统设计中具有重要的地位。部件级数学模型在各个部件模型的基础上，通过求解部件间共同工作方程
收稿日期:2017-10-20修回日期:2017-11-08组，来确定发动机相关性能参数，因而求解共同工作方程组算法的收敛能力、收敛效率和收敛精度，直接关系到模型的仿真性能"7。
因航空发动机的共同工作方程组实质是隐式的受约束的非线性方程组，其只能采取数值解法，同时模型本身对精度和实时性要求较高，这就对求解方法的收敛能力和收敛速度提出了较高的要求。
常用的求解非线性方程组的方法有Newton法和拟New-ton法。其中,Newton法因需要计算雅可比矩阵以及其逆矩阵，增加了计算的复杂度，同时算法迭代收敛与否还依赖于初始解的偏离程度和计算步长的大小。而拟Newton法是将雅可比矩阵作逼近处理，并通过递推获得逼近雅可比矩阵，
降低了计算复杂性，但是算法迭代收敛要求初始解的偏离程度更小，并受计算步长影响，收敛能力弱于Newton法。
在航空发动机部件共同工作方程求解中，目前主要采取Newton法H和Broyden拟Newton法皿山。陈玉春等对航空发动机数值模型计算过程的收敛性进行了分析，并采用独立变量初值限制法、变步长Newton法等改善求解的收敛性问题，取得了一定的效果，但仍需计算雅可比矩阵。王元「3)等提出了计算步长自校正Broyden拟Newton法，通过发散趋势判断实现Broyden拟牛顿法校正矩阵的重置，在迭代过程中有可能出现振荡收敛现象，影响模型的收敛速度。为解决模型收敛性和收敛速度的矛盾，提岀了一种混合变步长Newton-Steffenn-逆Broyden方法。同时为提高逆Broy-den方法的稳定性，
在逆Broyden的迭代计算过程中引入松弛因子〔⑷，从而构成具有局部超线性收敛能力的混合New-ton-Steffenn-逆Broyden(NS-iB)算法。在F404/YJ101变循环发动机数值仿真模型共同工作方程组求解中进行相关数值试验f,5-,6]，结果表明该算法能有效降低计算的复杂度，提高迭代解的精度与模型的实时性。
2变循环发动机气动热力系统模型
本文研究所用航空发动机气动热力系统模型是基于美
国通用航空与NASA改装的F404/YJ101变循环发动机建立的⑴-冏，结构如图1所示。
图1变循环结构示意图
该类型发动机稳态性能参数通过求解由7个共同工作条件所对应的非线性方程组获取，该方程组为
/.=-1
fl=(Q c DFS+-1
ft= Q h PT.<；/Q hp T,R—1
ft= Q lpt.c^Q lpt.k-1(1)
fi=(Hc”s*〃”《:)/(%'H hpt)-1
fb=-1
fl=PSd/pS"-1
式中,<?8为喷管气流流量,0；为加力燃烧室出口气流流量，Qcms为核心机驱动风扇气流流量为通过模式选择活门的气流流量，0皿为风扇气流流量.Q hpt,g为高压涡轮气流流量猜值,Q hpt.k为高压涡轮进口气流流量,Q ut.c为低压涡轮气流流量猜值,<2片,”为低压涡轮进口气流流量,H cms为核心机驱动风扇消耗功率，丹”代为压气机消耗功率，为高压涡轮产生功率,％为高压转子传动效率,//柿为风扇消耗功率，也刃为低压涡轮产生功率,％为低压转子传动效率, ps«3为外涵道出口静压,p%3为低压涡轮出口静压。
与之对应的热力循环中的7个未知参量为：风扇增压比叫”，核心机驱动风扇增压比Fcws，压气机增压比叭必，高压涡轮落压比77”叶，低压涡轮落压比77片以及高压转子换算转速pm和低压转子换算转速P”。以此构成该发动机的稳态数学模型，即通过求解非线性偏差方程组F(X“>)=[f t,
=0，以确定模型热力循环参数以及相关性能参数，式中X=[Fan»FCDFS,HPC>HPT»LPT^P n f»pncV
o而该模型的动态工作由5个共同工作条件构成的残差方程组G(X(“)=[£,£，厶，办，易]「=0，即与两个转子动力学方程所决定
dNg—%.Hb-H”5(2)
~dT=N—J「(”/30)2
°闩”斤-日”皿—H cms(3)
弘=几■(”/30)2
式中,J H和J L分别为高/低压转子转动惯量。与之对应的热力循环中的5个未知参量为X®= [丙臥ms E h pc L h e E L PT~\T。
根据上述内容，通过求解部件共同工作方程组，即可求出发动机相关工作状态参数，用以监测、控制和健康管理等任务。
3航空发动机共同工作方程组求解方法
求解航空发动机共同工作方程组实质是求解非线性方程组，因而广泛使用Newton法进行求解。在求解
精度较高时，因需要反复调用模型计算雅可比矩阵，该方法计算量大，模型实时性较差。为此，提出了一种混合Newton-Steffenn-逆Broyden方法用于共同工作方程组求解。
3.1Newton-Steffenn方法
对于非线性方程组
F(X)=0(4)采用Newton法求解时，计算步长及变量的初始值关系到方程的求解效果。当初始值一定时，计算步长关系到了算法的收敛性和收敛效率。
目前普遍应用于航空发动机数值模型求解的非线性方程组求解方法是前后双向差分的离散Newton型迭代法，其迭代公式为
=X®-[J(X(t),A(4))]-'F(X<t))
/1\
誌应(X®+Z
1-£(肝-硏乜)]>
(5)
式中,e，为第i个坐标为1的坐标向量是依赖于X®的已知向量。酒吧dj打碟
因为前后双向差分离散Newton法需要多次热力模型计
算，采用前向差分的离散Newton - Steffenn 法，有助于减小式中，△为算法选择阈值，一般A > l o 为计算步长,y®
算法的计算量，该算法迭代公式为
为算法选择变量
KX* ,h w )
(6)
II “ II 2
(10)
式中，碑> =h, || F(X® ) ||，这里& # 0为常数。
由于Newton 算法在迭代过程中需要计算雅可比矩阵J , 对于许多问题，这是一件困难且耗时的事情。
而航空发动机 的数学模型要求较高的实时性，为避免或降低对雅可比矩阵 的依赖，基于割线法衍生而来的拟Newton 法既保留了 Newton 法的局部超线性收敛能力，又能降低计算复杂度。因
此在许多实际问题的求解过程中，都有采用拟Newton 法进 行问题求解的实例［6'710
即当y® > A 时，使用式(6)前向差分求取雅可比矩阵
即:= (_/(於*>，肚>))7，否则按式(8)
更新B <ul,o 由式(11)确定计算步长
0.9 t > 0. 9
(11)
0. 1
7 < 0. 1
式中,T 为
II F(XE>) ||II F(X “>) ||
(12)
对于航空发动机数值仿真模型来说，模型精度要求较3.2 逆 Broyden 方法
Broyden 拟牛顿法(简称Broyden 法)是Broyden 于1965
年首次提出，并迅速发展，形成Broyden 族方法加-山。其中 逆Broyden 方法以校正矩阵B 逼近逆雅可比矩阵J"求解方 程组，通过递推进行校正矩阵B 的计算，降低了算法的复杂高，算法终止条件通常设为方程组残差分量均小于某一数 值，如le - 5，但如果算法陷入局部最优，将无法达到精度要 求。由此，为避免这种陷入局部最优而无法达到精度要求的
情况出现，同时结合文献［14］中对Broyden 算法稳定性的研
究，本文将松弛因子8引入到逆Broyden 方法中，得到校正矩性。
Broyden 法秩1迭代方法为
何“+” = X ⑷-=
t (y _ 人 s )(s )
阵的计算为
(13)
(7)
式中丿⑴为雅可比逼近矩阵,s"> = X"" -X®，要求s® # 0,y (t) = F(X"*") -F(X">)o 应用 Sherman - Morrison
定理，在<4可逆时，可得逆Broyden 法
乂“" = X W _ B {k}F(X w )
B<“ = b ⑷ +(代-B®产)(s “>)W"
(8)
I I M a I y (t) I < a
(14)
式中,o ■属于(0,1 ) ,sign(0) = l 0
综上所述，本文所提出的具有变步长和自适应松弛参数
的混合 Newton - Steffenn -逆 Broyden 方法(NS - iB)求解
过程中步骤如下：
式中，校正矩阵 B = A-',3.# 0o
可见逆Broyden 方法中，不需要计算逆矩阵，直接递推获 得校正矩阵B"”，简化了计算，是研究的热点［18-,91…但是
拟Newton 法的收敛范围低于Newton 法。而Newton 法的计算 工作量又明显大于逆Broyden 方法。因此，如何结合两者优 势成为本文的主要研究内容。
3.3 混合 Newton - Steffenn -逆 Broyden 方法
方程组收敛与否与初始解的偏离程度和计算步长的大
步骤1：初始化
给定算法选择阈值A，方程组解的初值雅可比矩阵
计算差分系数向量胪。将初值代入部件模型计算残差 方程组，并计算雅可比矩阵J(X <0) ,/»(0))o
步骤2：逼近矩阵B (l+,)计算
if\ I > A，求取雅可比矩阵,h w )
= ,h w )y',
小关，当方程组残差2范数较大时，意味着解的偏离较大，采 用较强收敛能力的Newton 法求解，而残差2范数较小时，意
味着解的偏离较小，则采用较高收敛效率的逆Broyden 法，并 在算法中采用具有偏离自适应能力的可变计算步长，在偏离 较大时采用较小的计算步长，进一步保障模型的收敛能力，
而偏离较小时，则采用较大的计算步长，提高算法的收敛效 率，从而得到混合Newton - Steffenn -逆Broyden 方法
X (i+,) = X (k) - A (k)B (k)F(X (k))
B “+" = I y w I > d (9)
i el
B (U» = B “> + 护(H -B®严)
(s “>)肿严
松弛因子e，按照式(13)计算。
步骤3：按照式(9)迭代计算Xg" ,k = lo
步骤4：将川小> 代入部件模型计算残差方程组及其2范
数 I I F(X"">) I I?。
步骤5 ：如果满足收敛条件或达到最大允许迭代次数，结
束本工作点计算，否则，转到步骤2。
4实验结果分析
在Visual Studio 2013中，使用C#编写数值模型应用程序
吕⑷+磧e,)、
1"(肝)］
丿”
T
1
]_ sign(y “> )<r
1 _严
开展仿真。在Windows10企业版.E5-1620v3,16GB RAM, 256GB SSD的工作站中运行1,000,000次变循环发动机热力模型，总耗时92,022ms，平均1次热力循环模型计算耗时0. 092ms o
将NS-iB方法应用于稳态模型求解中，与定步长X=0.4的Newton-Steffenn法(0.4-N-S)、定步长X=0.8的New-ton-Steffenn法(0.8-N-S)、变步长Newton-Steffenn法(VS-N-S)、使用Newton-Steffenn初始化的标准逆Broyden 法(S-iBroyden)、使用Newton-Steffenn初始化的松弛逆Broyden法(R-i Bro
yden)以及混合Newton-Steffenn-逆Broyden法(NS-iB)进行比较，均以残差方程组各分量的绝对值均小于le-5为迭代终止条件，为确保模型收敛，最大允许迭代次数设为1000次。
在不同初始偏离情况下，稳态计算结果如表1所示，N 为求解过程中调用模型热力循环计算的次数,Iter为求解迭代次数,NC表示求解发散。
表1不同求解方法计算结果
〃F(X0)〃
0.4-N-S0.8-N-S VS-N-S S-iBroyden R-iBroyden NS-iB
N Iter N Iter
N Iter N Iter N Iter N Iter
0.0124511314415577157157157
0.06309137174967391791792510
0.083137174966581791792510
0.087713717496739221422142611
0.0958137174968110201220122611
0.115145184968110191119112611
0.11936153195778911241624162813
0.121145184968110NC20122712
0.203153195778911NC38302914
0.28827161205779712NC47393015
皮肤病治疗0.33976161205778911NC NC2813
0.447516921NC9712NC NC3015
0.533717722NC9712NC NC3116
0.547917722NC9712NC NC3015
故宫旅游攻略
0.6717722NC10513NC NC3116
调用模型热力循环计算次数N和迭代次数Iter的统计
图如图2所示，也就是采用不同方法更新逆雅可比矩阵时，
迭代次数是相同的，但模型调用次数不同，根据式(6)计算时
需要重新调用模型计算雅可比矩阵，调用次数等于猜值个
数，而根据式(13)递推则不需要额外调用模型。
从表1中的数据对比知,NS-iB方法在发动机模型共同
工作方程求解中具有明显的优势。S-iBroyden方法收敛范
围最小,Newton法的收敛范围广。增大计算步长，可提高算
法的计算效率，但大的计算步长在求解中收敛范围变窄。使
用变步长算法可有效解决收敛速率与收敛范围的问题。本
文研究的NS-iB算法在计算量上，部件模型调用次数约为小
步长0.4-N-S算法的1/5，大步长0.8-N-S算法的1/2，变
步长牛顿法的1/3。通过与不同算法的对比，充分验证了本
文研究的混合算法的收敛范围广和收敛速率快的特点。
为了进一步验证算法的适用性，在发动机动态模型中开
展仿真验证。用NS-iB算法有限次通过(最大允许迭代次数
为10次,10-NS-iB)方法求解动态模型，其收敛条件设为残
差方程组分量的绝对值均小于le-5，并与变步长Newton-Steffenn算法有限次通过(最大允许迭代次数分别为2次、5比次和10次，即2-N-S、5-N-S和10-N-S)，通过仿真展开对图2Iter和N的统计示意图
鸡蛋葱花饼在相同的稳态工作点上给出不同的指令，模拟发动机动(Os至8s)内累计热力模型调用次数N，最大残差分量绝对态运行过程，以20ms为仿真步长，记录在固定仿真时间段值Max。不同的动态过程仿真结果如表2所示。
表2不同动态过程计算结果对比
动态过程
2-N-S5-N-S
Max 10-N-S10-NS-iB
Max N
N N Max N Max 地面pnf100%-95%17768.36E-42334 5.11E-52478 1.0E-510419.91E-6地面pnf100%-83%2178 3.95E-32952 2.29E-434629.96E-614529.97E-6高空pnf98%-93%3018 4.36E-34882 2.87E-44056 1.0E-516439.98E-6高空pnf98%-82%3906 3.48E-25550 3.58E-45802 1.0E-522649.98E-6如表2所示，在变循环发动机动态模型中,2-N-S通过
算法由于对解迭代更新较少，解的偏离程度较大，动态仿真
过程会引起较大偏离。而10-N-S算法能够对解进行多次
修正，使解的精度达到或接近阈值，模型的精度得到提高。
而本文所提出的NS-iB方法，由于混合了变步长及递推的梯
度计算，在有限次通过中，计算量和收敛效率、收敛精度明显
优于有限次通过的变步长Newton-Steffenn算法，说明了该
算法在模型动态仿真过程中同样保持较低的计算量和高效
的收敛效率，且可以获得较高的收敛精度。
图3是表2中地面pnf100%-95%的动态仿真过程方
程求解情况，图4是表2中高空pn©8%-82%的动态仿真过程方程求解情况。
从表2和图3、图4综合分析知，在模型动态仿真过程中,NS-iB不仅有高于2-N-S算法的实时性，计算精度比10 -N-S算法高，进一步说明了算法在模型动态求解过程中的优越性。
5总结
本文针对变循环发动机数学模型共同工作方程的求解方法开展研究工作。对常用的Newton法和逆Broyden法进行了介绍，提岀了一种混合Newton-Steffenn-逆Broyden求解方法，利用Newton算法的强收敛能力和逆Broyden算法的高实时性，在不同偏离下采用不同的求解方法，并引入变步长，进一步提高模型的收敛能力和实时性，同时为提高逆Broyden算法的稳定性，将松弛因子引入逆Broyden的递推计算过程。在发动机稳态模型和动态模型的数值仿真表明，本文所提出的NS-iB算法在保持逆Broyden算法结构简单、局部超线性收敛等优势的同时，保持较高的精度的解。
在F404/YJ101变循环发动机的数值模型中的稳态计算试验中表明,NS-iB法求解过程中调用模型热力循环计算的次数约为小步长牛顿法的1/5，大步长牛顿法的1/2，变步长牛顿法的1/3，并接近以高实时性著称的逆Broyden方法，充分验证了算法的有效性。在模型的动态仿真中，在不同指令下的动态响应过程的分析中,NS-iB法在低于变步长Newton -Steffenn有限次通过法的计算量的同时，计算精度也有提
(a)Max|F(A»>)1|K$ft
(b)筋力模型的次敦N梦见哭泣
图3地面pnf100%~95%的动态变化仿真过程
养颜糖水
高。证明了NS-iB法在降低计算量的同时，提高了解的精度，在具有更多可调参数、工况复杂的新型发动机建模过程中，具有良好的应用前景。
参考文献:
[1]Garg Sanjay.Aircraft Turbine Engine Control Rearch at NASA
Glenn Rearch Center[J].Journal of Aerospace Engineering, 2013,26(2):422-43&
[2]潘阳，李秋红，王元.基于Broyden算法的航空发动机气路故障
诊断[J].推进技术,2017,38(1):191-198.
[3]王元，李秋红，黄向华•基于DMOM算法的航空发动机性能寻
题目自拟的作文
