全国研究生数学建模竞赛参赛作品解读

2023年12月30日发(作者：天津2013高考数学试卷)

参赛密码

（由组委会填写）

第十届华为杯全国研究生数学建模竞赛

学 校

参赛队号

队员姓名

参赛密码

（由组委会填写）

第十届华为杯全国研究生数学建模竞赛

题 目

变循环发动机部件法建模及优化

摘 要：

本文利用附录1、2、3、4给出的特性数据以及计算公式得出了风扇特性表中流量随风扇压比函数值的变化图形规律，求出了给定条件下风扇和CDFS出口的总温、总压和流量，并进一步建立了发动机非线性平衡方程组求解的遗传算法模型，并得到了给定条件下非线性方程组求解结果，进而了解了给定条件下变循环发动机双涵道模式的工作状况；同时还给出了计算发动机性能最优对应的发动机CDFS导叶角度、低压涡轮导叶角度和喷管喉道面积的求法。

针对问题一。第一小问，借用附录3中的压气机压比函数值定义式，求得附录4中风扇特性数据表中各换算转速下增压比对应的压比函数值（见附表一），然后由所得压比函数值及表中流量数据画出了流量随压比函数值变化的图形(见图5-1、5-2)；第二小问，首先借助题中给定的物理转速及压比函数值利用附录4中风扇和CDFS特性数据表插值并利用附录一中公式（2.7）进行修正得到二者对应的增压比prc、效率c和换算流量Wc，由于风扇进口总温、总压=进气道出口总温、总压，CDFS进口总温、总压=风扇出口总温、总压。借助修正的增压比、效率、换算流量及题中所给初始条件并利用附录中的相关公式可依次得到了风扇和CDFS的出口总温、总压及流量分别为378.333、1.288、19.048；431.803、1.774、16.940。

针对问题二。首先利用附录1、2、3、4中的信息得到了题中7个非线性方程的基本参数表达形式。将题中的初始条件代入得到了仅含nH、ZCL、ZCDFS、ZCH、ZTH、ZTL、T4*七个基本参数的具体形式非线性方程组。本文采用了遗传- 1 -

算法对所得非线性方程组进行了求解，依次进行随机生成初始化群体（即解的初始值）、选择、交叉、变异计算、停止判断操作。这些流程通过C语言编码实现。通过计算得到所求非线性方程组解为nH=0.87,ZCL=0.70,ZCDFS=0.79,

ZCH=0.31,ZTH=0.52,ZTL=0.22,T4*=1419.9指明了问题二给定条件下的变循环发动机双涵道模式运行工作状况。

针对问题三。根据本问给定初始条件，可调未知参数共有11个，按照飞机发动机常用的调节规律，选择其中的4个参数赋予初值。实际发动机非线性数学模型中，通常采用数值计算方法，首先给出一组7个独立变量的初值，代入发动机模型，计算各平衡方程的误差，通过一定的非线性方程组迭代解法，不断地修正独立变量的试给值，直至各平衡方程的误差满足设定的迭代精度要求，从而得出各个参数的值，然后由所得参数值，结合附录1中计算发动机性能参数中的推力、单位推力、耗油率公式对变循环发动机的性能进行取优评价，此时对应的发动机CDFS导叶角度、低压涡轮导叶角度和喷管喉道面积即为所求。

关键词：变循环发动机；非线性方程组；遗传算法；单/双涵道模式

- 2 -

一、问题重述

1.1 问题背景

由飞机/发动机设计原理可知，对于持续高马赫数飞行任务，需要高单位推力的涡喷循环，反之，如果任务强调低马赫数和长航程，就需要低耗油率的涡扇循环。双涵道变循环发动机可以同时具备高速时的大推力与低速时的低油耗。变循环发动机的内在性能优势，受到了各航空强国的重视，因此有必要对变循环发动机部件法建模问题进行深入研究，进而推动航空事业的发展。

1.2 现有信息和数据

变循环发动机有两种工作模式，分别为涡喷模式和涡扇模式。发动机在亚音速巡航的低功率工作状态，风扇后的模式转换活门因为副外涵与风扇后的压差打开，使更多空气进入副外涵，同时前混合器面积开大，打开后混合器，增大涵道比，降低油耗，此时为发动机的涡扇模式。发动机在超音速巡航、加速、爬升状态时，前混合器面积关小，副外涵压力增大，选择活门关闭，迫使绝大部分气体进入核心机，产生高的推力，此时为发动机的涡喷模式。

燃气涡轮发动机是由进气道、压气机、主燃烧室、涡轮、喷管等部件组成的。如果计算机能够对这些部件的性能进行准确的模拟，那么也就能准确地模拟整个发动机的性能。这种建立在准确模拟发动机各部件性能基础上的发动机性能计算方法，称为部件法。该方法是建立在发动机各部件特性已知的基础上的，因此是计算精度较高的一种方法。

发动机各部件匹配工作时，受低压轴功率平衡方程、高压轴功率平衡方程、高压涡轮进口截面流量平衡方程、低压涡轮进口截面流量平衡方程、后混合器静压平衡方程、尾喷管面积平衡方程和风扇出口流量平衡方程7个平衡方程的制约，同时该问题给出了发动机部件计算公式、工质热物理性质参数、气动函数和压气机及涡轮特性等数据。

1.3 需要解决的问题

1.3.1 流量随压比函数值变化图形的求解以及风扇和CDFS的出口总温、总压和流量的求解。

（1）、请画出附录4中风扇特性数据表中流量随压比函数值变化的图形。

（2）、设在发动机飞行高度H11km，飞行马赫数Ma0.8的亚音速巡航点，导叶角度均设置为0°，风扇和CDFS的物理转速都为0.95，风扇和CDFS的压比函数值都为0.5，求风扇和CDFS的出口总温、总压和流量。

1.3.2 请运用或设计适当的算法求解由发动机7个平衡方程组成的非线性方程组。要求陈述算法的关键步骤及其解释，尽可能讨论算法的有效性。

设在发动机飞行高度H11km，飞行马赫数Ma0.8的亚音速巡航点，采用双涵道模式，导叶角度均设置为0°，选择活门完全打开，副外涵道面积设为1.8395e+003，后混合器出口总面积设置为2.8518e+004，尾喷管喉道面积A89.5544e+003，nL=0.85。请运用或设计适当的算法求解由发动机7个平衡方- 3 -

程组成的非线性方程组。要求陈述算法的关键步骤及其解释，尽可能讨论算法的有效性。

1.3.3 发动机性能最优的求解以及尾喷管喉道面积随飞行马赫数变化规律的求解。

（1）、设在发动机飞行高度H11km，飞行马赫数Ma1.5的超音速巡航点，发动机采用单涵道模式，将选择活门面积设置为0，风扇导叶角度、高压压气机导叶角度、高压涡轮导叶角度均设置为0，后混合器面积设置为2.8518e+004。请问发动机CDFS导叶角度、低压涡轮导叶角度和喷管喉道面积3个量为多少时，发动机的性能最优?

（2）、试研究发动机飞行高度H11km，飞行马赫数从Ma1.1变化到Ma1.6，发动机特性最优时，CDFS导叶角度、低压涡轮导叶角度，尾喷管喉道面积随飞行马赫数的变化规律。此时发动机采用单涵道模式，将选择活门面积设置为0，风扇导叶角度、高压压气机导叶角度、高压涡轮导叶角度均设置为0，后混合器出口总面积设置为2.8518e+004，后混合器内、外涵道面积可调（即不受附录1后混合器给定的内、外涵道面积值的约束）。

二、问题分析

2.1 问题一

针对本问题。第一小问可借用附录3中的压比函数值定义式，求得附录4中风扇特性数据表中增压比对应的压比函数值，然后由所得压比函数值及流量数据画出流量随压比函数值变化的图形。第二小问采用附录1中发动机进气道、压气机的计算公式，依据假定1按进气道-风扇-CDFS的顺序求解[1]。其中风扇和CDFS均为压气机部件，在计算中特征流量、增压比、效率均为核心转数、压比函数值、导叶角角度的函数，需采用附录4中风扇特性数据表和CDFS特性数据表插值求解。最后即可得到总温、总压和流量。

2.2 问题二

针对本问题。分两步求解，首先依据附录1中的发动机部件计算公式再代入附录2、3中的辅助公式列出题目中给出的7个发动机平衡方程，建立仅含nH、ZCL、ZCDFS、ZCH、ZTH、ZTL、T4*七个变量的待求的非线性方程组。然后针对得到的非线性方程组的特点，选用标准遗传算法，按照算法的步骤编程求解。

2.3 问题三

当给定发动机部件特性时，同时即描述了发动机的工作过程，但由已建立的方程组唯一确定满足发动机共同工作方程的解，还需给出发动机的调节规律- 4 -

[7]。针对本问题。本文采用2组具有合理初值的参数分别计算得到发动机各参数的最优解，从而得到优化发动机CDFS导叶角度、低压涡轮导叶角度和喷管喉道面积3个量的方法，使发动机性能最优。

在马赫数从Ma1.1变化到Ma1.6过程中，若算法计算较慢，则采用几个离散的马赫数求得对于的最优发动机CDFS导叶角度、低压涡轮导叶角度和喷管喉道面积，然后寻找各自的变化规律；若能得到快速收敛的算法，则取较小马赫数增量步长连续得到一系列最优发动机CDFS导叶角度、低压涡轮导叶角度和喷管喉道面积，得到随马赫数变化的规律，分别画出三个参数随马赫数变化的图形。

三、基本假设

通过阅读有关参考资料和对数据进行简单的分析，我们对该变循环发动机模型做如下假设：

3.1 高压压气机后不经主燃烧室的分流气流为冷却气流，在本题中忽略不计。

3.2 如果计算机能够对变循环发动机部件的性能进行准确的模拟，那么也就能够准确地模拟整个发动机的性能。

3.3 在前混合器计算中，假定CDFS涵出口流量，即为通过CDFS特性数据线性插值得到的CDFS计算流量减去高压压气机特性数据线性插值得到的高压压气机计算流量。

3.4 如无特殊要求，发动机中直接相连的下一级主要部件进口总温、总压假定为上一级主要部件出口处的总温、总压（给出总压恢复系数的另外考虑），直接相连的发动机主要部件顺序如下：风扇-CDFS-高压压气机-主燃烧室-高压涡轮-低压涡轮。例如：风扇计算得到风扇出口的总温Tout总压pout1，这一总温、1、总压即为CDFS的进口总温、总压，以此类推。

四、符号说明

进气道参数：进气道进口总温 T0f1(H)进气道进口总压 p0f2(H,Ma)

进气道总压恢复系数

if3(Ma)进气道出口总压 ppi10

- 5 -

风扇参数：风扇进口总温 TinT1风扇进口总压 pinp1*风扇换算转速 ncorf4(nL,Tin)

风扇特性 prcf5(nL,ZCL,L)

cf6(nL,ZCL,L)换算空气流量 Waf7(nL,ZCL,L)**风扇出口总压 pout1pinprc*风扇出口总温 Toutf(T18in，c，prc）空气流量 WCLf9(T，p，Wa)风扇消耗功 lCLf10(Tin，c，pr)风扇消耗功率 NCLf11(Tin，c，pr)WCLin*in

CDFS参数：*CDFS换算转数 ncorf12(nH,Tout1)CDFS特性 prcf13(nH,ZCDFS,CDFS)

cf14(nH,ZCDFS,CDFS)CDFS换算空气流量 Waf15(nH,ZCDFS,CDFS)**CDFS出口总压 pout2pout1prc*CDFS出口总温 Tout2f16(Tout1，c，prc）*CDFS空气流量 WCDFSf17(Tout，p1out1，Wa)CDFS消耗功 lCDFSf18(Tout1，c，pr)CDFS消耗功率 NCDFSf19(Tout1，c，pr)WCDFS

高压压气机参数：*高压压气机换算转数 ncorf20(nH,Tout2)高压压气机特性 prcf21(nH,ZCH,H)

cf22(nH,ZCH,H)高压压气机换算空气流量 Waf23(nH,ZCH,H)**高压压气机出口总压 pout3pout2prc*高压压气机出口总温 Tout3f24(Tout2，c，prc）*高压压气机空气流量 WCHf25(Tout2，pout2，Wa)高压压气机消耗功 lCHf26(Tout2，c，pr)高压压气机消耗功率 NCHf27(Tout2，c，pr)WCH

燃烧室参数：**燃烧室出口总压 pout4pout3b油气比 fbf28(T*out3,T)*4

主供油量 WfWCHfb- 6 -

高压涡轮参数：*高压涡轮流量 WTHf29(T4*,pout4,Wc)涡轮特性

cf30(nH,ZTH,H) prcf31(nH,ZTH,H)

**高压涡轮出口总温 Tout5f32(T4, prc,c)**高压涡轮出口总压 pout5f33(pout4, prc)低压涡轮参数：**低压涡轮燃气流量 WTLf34(Tout5,pout5,Wc)涡轮特性

cf35(nL,ZTL,CH) prcf36(nL,ZTL,CH)

**低压涡轮出口总温 Tout6f37(Tout5, prc,c)**低压涡轮出口总压 pout6f38(pout5, prc)前混合器参数：CDFS涵出口流量 Wg125WCDFSWCH后混合器参数：后混合器外涵出口流量 Wg62W15后混合器出口总温 T6*后混合器出口总压 p*6

后混合器出口流量 Wg6总压恢复系数

duct0.98加力燃烧室（略）

尾喷管参数：*尾喷管进口总温 TinT6*尾喷管进口总压 pp*in*6

尾喷管进口流量 WgWg6尾喷管喉部面积 A8尾喷管喉部流量 Wg8

尾喷管尾部面积 A9尾喷管出口静温 T9尾喷管出口气流速度 c9

- 7 -

五、模型建立和求解

5.1 求解问题一

5.1.1 风扇特性数据表中流量随压比函数值变化的图形

由附录3中的压气机压比函数值的定义可知，压比函数值zz可以如下算得：

设压气机某换算转速所对应的增压比数据（见附录4）的最大值为prmax，最小值为prmin，则定义该换算转速对应的压气机增压比pr的压比函数值zz为：

zz=prprmin （5.1.1-1）

prmaxprmin由附录4中风扇特性数据，利用公式（5.1.1-1）可分别算出换算转速为0.4、0.5、0.6、0.7、0.81、0.9、0.95、1和1.075对应的压比函数值。然后利用得到的压比函数值数据与流量数据可以得风扇流量随压比函数值变化的图形。结果如下所示：

流量随压比函数值变化(一)8010000.10.20.30.40.50.60.70.80.9压比函数值11.1换算转速 0.4换算转速 0.5换算转速 0.6换算转速 0.7换算转速 0.81流量

图5-1 风扇流量随压比函数值变化图形(一)

流量随压比函数值变化（二）120000.10.20.30.40.50.60.70.80.9压比函数值11.1换算转速 0.9换算转速 0.95换算转速 1换算转速 1.075流量

图5-2 风扇流量随压比函数值变化图形(二)

- 8 -

5.1.2 给定条件下风扇及CDFS出口的总温、总压和流量

（1）、计算思路

由附录1中公式（2.1）～（2.4）算出进气道的出口总温总压，由“风扇进口总温总压=进气道出口总温总压”结合风扇的特性数据算出风扇出口总温总压，再利用“CDFS进口总温总压=风扇出口总温总压”结合CDFS特性数据及相关公式算出CDFS的出口总温总压；风扇和CDFS出口的流量则可以直接利用附录4中的特性数据进行插值得到。

（2）、给定条件

题中给定的初始条件有：

飞行高度H=11km；

飞行马赫数为Ma=0.8的亚音速巡航点；

各导叶角度均设置为0°；

风扇和CDFS物理转速均为0.95；

风扇和CDFS的压比函数值均为0.5。

（3）、算法流程

①进气道出口总温总压的计算

利用附录1中的公式（2.1）～（2.4）结合上面给定条件算出进气口的总温、总压分别为：244.38、0.345。

②风扇出口总温总压及流量的计算

A、总温总压计算

首先，计算修正后的增压比prc、效率c和换算流量Wc。

利用给定的物理转速n=0.95及附录1中的公式（2.5）计算出换算转速ncor=1.032；利用ncor及给定的压比函数值zz=0.5从附录4中风扇特性数据表线性插值得到增压比prc,map=2.146、效率c,map=0.782和换算流量Wc,map=104.157；再结合附录1公式（2.7）算出修正后的增压比prc=3.737、效率c=0.836和换算流量Wc=51.558。

**之后，计算出口总温Tout和总压pout。

*由假设有：风扇进口总压pin=进气道出口总压p1，则有：

***pout=pinprc=p1prc

*由假设有：风扇进口总温Tin=进气道出口总温T1，再结合附录1中“2.2.2 计*算过程——>3) 计算压气机出口参数”的公式算法步骤算出风扇出口总温Tout。

B、出口流量Wa2计算

利用附录4中的风扇特性数据线性插值并修正求出来的流量Wc结合附录1中的公式（2.8）求出风扇出口流量Wa2=19.048。

③CDFS出口总温总压及流量的计算

按与风扇同样的方法步骤算出修正后的增压比prc=1.377、效率c=0.875和换算流量Wc=43.813；

**由于CDFS进口总压pin=风扇出口总压，CDFS进口总温Tin=风扇出口总温。再利用与风扇总温总压同样的求解方法可以得出CDFS的出口总温总压；

同样采用附录4中CDFS特性数据线性插值并修正获得的流量Wc结合附录- 9 -

1中的公式（2.8）求出CDFS出口流量Wa24=16.940。

（4）、计算结果

经过计算得到风扇和CDFS的出口总温、总压及流量结果如下表：

表5-1 风扇和CDFS的出口总温、总压及流量

风扇

CDFS

5.2 求解问题二

出口总温

378.333

431.803

出口总压

1.288

1.774

出口流量

19.048

16.940

如无特殊要求，发动机中直接相连的下一级主要部件进口总温、总压假定为上一级主要部件出口处的总温、总压（给出总压恢复系数的另外考虑），直接相连的发动机主要部件顺序如下：风扇-CDFS-高压压气机-主燃烧室-高压涡轮-低压涡轮。例如：风扇计算得到风扇出口的总温Tout1

、总压pout1，这一总温、总压即为CDFS的进口总温、总压，以此类推。

5.2.1 展开发动机平衡方程组

（1）、展开发动机平衡方程组一

NCLNTLmL0

（5.2.1-1）进气道计算：5.2553H

p01.013251 一定高度下大气压力、温度44.308T288.156.5H012TT1Ma002 进气道进口总温、总压

112p0p012MaMH1：i1.0 进气道总压恢复系数1.35MH1：i1.00.075(MH1)T1T0

进气道出口总温、总压

pp10i风扇计算：

TinT1

风扇进口总温、总压pp1in

- 10 -

lchouthin 风扇消耗功、功率

NcWalc**Tinkwpin,dWCLCwW(nL,ZCL,L)(1L)* 风扇流量*100Tinpin,d**lCLh(Tout)h(T1in) 风扇消耗功

NCLWCLlCL 风扇消耗功率**(Tout(Tin),ei)*Tout,ei*1(Tout,ei) 风扇出口理想温度kprRlnCpr(prc(nL,ZCL,L)1)(1L)1 风扇出口理想熵M100**h(T)h(T)out,eiin**h(Tout1)h(Tin) 风扇出口焓k2Cc(nL,ZCL,L)(1L)100*1*Tout1h(Tout1) 风扇出口总温

p*out1kprpCpr(prc(nL,ZCL,L)1)(1L)1 风扇出口总压100*in低压涡轮运算：

lThinhout 低压涡轮功、功率NWlgTmTTinTout5

低压涡轮进口总温、总压

ppout5in**Tinkwpin,dWTLCwW(nL,ZTL,CH)(1CH)* 低压涡轮流量*

100Tinpin,d****lTLh(Tin,fb)h(Tout6,fb)h(Tin,fb)h(Tout6,fb) 低压涡轮功NTLWTLlTLm 低压涡轮功率

kpr**poutpC(pr(n,Z,)1)(1)16inprcLTLCHCH 低压涡轮出口总压

100T*out6prc1kprC(pr(n,Z,)1)(1)1prcLTLCHCH 低压涡轮出口总温

100（2）、展开发动机平衡方程二

NCHNCDFSNTHmH0

（5.2.1-2）- 11 -

CDFS计算：

lchouthin CDFS消耗功、功率NWlacCDFSTinTout1

CDFS进口总温、总压

pinpout1**Tinkwpin,dCwW(nH,ZCDFS,CDFS)(1CDFS)* CDFS流量

*100Tinpin,dWCDFS**lCDFSh(Tout2)h(Tin) CDFS消耗功率NCDFSWCDFSlCDFS CDFS消耗功率**(Tout(Tin),ei)

kprRlnCpr(prc(nH,ZCDFS,CDFS)1)(1CDFS)1 CDFS出口理想熵M100*1Tout(Tout*,ei) CDFS出口理想总温,eih(T

*out2*Tout2k2Cc(nH,ZCDFS,CDFS)(1CDFS)100*h1(Tout2) CDFS出口总温)h(T)*in**h(Tout,ei)h(Tin) CDFS出口焓

kpr**poutpC(pr(n,Z,)1)(1)12inprcHCDFSCDFSCDFS CDFS出口总压100高压压气机计算：

lchouthin 高压压气机消耗功、功率NcWalcTinTout2

高压压气机进口总温、总压

pinpout2WCH**Tinkwpin,dCwW(nH,ZCH,H)(1H)* 高压压气机流量

*100Tinpin,d**lCH（hTouth(Tin) 高压压气机消耗功3）NCHWCHlCH 高压压气机消耗功率**(Tout(Tin),ei)

kprRlnCpr(prc(nH,ZCH,H)1)(1H)1 高压压气机理想出口熵M100- 12 -

*1*

Tout,ei(Tout,ei) 高压压气机出口理想总温

h(T*out3**h(Tout,ei)h(Tin))h(T) 高压压气机出口焓

2kCc(nH,ZCH,H)(1H)100*in*1*Tout3h(Tout3) 高压压气机出口总温

kpr**poutpC(pr(n,Z,)1)(1)13inprcHCHHH 高压压气机出口总压

100主燃烧室运算：

TinTout3

主燃烧室进口总温、总压ppout3inh3Hair(Tin) 主燃烧室进口焓h4Hair(T)*4fb1fbHst(T) 主燃烧室出口焓*4

fbh4h3 油气比

bHuh3WfWa3fb 主燃油流量WoutWfWa3 主燃烧室总流量

Wa3WCH**Tinkwpout,dCwW(nH,ZCH,H)(1H)**2 主燃烧室空气流量100Tout2pin,d

***poutpp44in3b 主燃烧室出口总压Tout4 T 主燃烧室出口总温*4

高压涡轮运算：

lThinhout 高压涡轮功、功率NWlgTmTTinTout4

高压涡轮进口总温、总压

ppout4in- 13 -

WTH**Tinkwpin,dCwW(nH,ZTH,H)(1H)* 高压涡轮流量*100Tinpin,d****lTHh(Tin,fb)h(Tout5,fb)h(Tin,fb)h(Tout5,fb) 高压涡轮做功

NTHWTHlTHm 高压涡轮功率kpr**poutpC(pr(n,Z,)1)(1)15inprcHTHHH 高压涡轮出口总压100kpr*1ToutprC(pr(n,Z,)1)(1)15cprcHTHHH 高压涡轮出口总压

100*Toutprc11*Tin/ｃCp/R 高压涡轮落压比（3）、展开发动机平衡方程三

Wg41Wg410

（5.2.1-3）TinTout4

高压涡轮进口总温、总压

ppout4in**Tinkwpin,dCwW(nH,ZTH,H)(1H)* 高压涡轮总流量

*100Tinpin,dWTHW41W41WfWa3 主燃油流量和经主燃烧室空气总流量的和Wa3WCH

**Tinkwpout,dCwW(nH,ZCH,H)(1H)**2 高压压气机总流量100Tout2pin,dh3Hair(Tout3) 主燃烧室进口焓h4Hair(T)*4fb1fbHst(T) 主燃烧室出口焓*4

fbh4h3 油气比bHuh3

WfWa3fb 主燃油流量（4）、展开发动机平衡方程四

Wg45Wg450

（5.2.1-4）

**Tinkwpout,dWTLCwW(nL,ZTL,CH)(1W45CH)**5 低压涡轮总流量

100Tout5pin,dW45W41

- 14 -

（5）、展开发动机平衡方程五

p61p620

（5.2.1-5）前混合器运算：q(125)*Wg125T125*kmp125A125 前混合器CDFS涵气动函数q

A125= 608.4252 前混合器CDFS涵面积km0.0404 空气流量系数**p125pout2CDFS 涵道出口总压**T125Tout2CDFS 涵道出口总温

Wg125WCDFSWCH CDFS总流量与高压压气机总流量的差**Tinkwpout,dCwW(nH,ZCDFS,CDFS)(1CDFS)**1 CDFS总流量100Tout1pin,dWCDFSWCHCwW(nH,ZCH,H)(111TkwH)100T*in,d*out2

pp*out2*in,d 高压压气机总流量q(125)(12)125（1-1-11） 气动函数q-12125125q1q(125) 速度系数*225

*p125p125(125) CDFS涵出口静压p225p125p（225）副外涵出口静压2-1

（225）（1-2251-11225） 气动函数(225) 速度系数*p225

Wg225kmT*225A225q(225) 前混合器副外涵出口流量

**ToutTin

**poutpinduct 涵道进出口总温、总压、总流量关系

Wa,outWa,in

A225=1.8395e+003 副外涵出口面积km0.0404 空气流量系数

- 15 -

**p225pout1duct 副外涵出口总压T*225T*out1 副外涵出口总温

Wg225Wg125Wg15

WhWhWh

g125125g1515g225225p*f()AP*f()Ap*f()(AA)22522515225125225*p15Wg15kmA15q(15)

*T15前混合器总方程组，通过此方程组可求得前混合器混合后出口的总温、总压、总流量f()/q()

Z()2111Z() 气动函数

111 气动函数11-12-1q()(2)（1-1） 气动函数后混合器运算：**ToutTin

**pp 涵道进口出口总温、总压、总流量关系outinductWa,outWa,in

q(61)*Wg61T61kmpA61*61 气动函数A61= 5.3061e+003 内涵出口面积

km0.0397 燃气流量系数**p61pout6duct 内涵出口总压TT*61*out6 内涵出口总温

Wg61Wa,in**Tinkwpout,dWTLCwW(nL,ZTL,CH)(1CH)**5 内涵出口总流量100Tout5pin,d1111-12-1q(61)(2)61（1-） 气动函数161

61q1q(6)1 速度系数p61p61()6

1内涵出口静压*

- 16 -

**ToutTin

**poutpinduct 涵道进口出口总温、总压、总流量关系

Wa,outWa,in

q(62)*Wg62T62kmpA62*62 气动函数

A62= 2.3212e+004 外涵出口面积km0.0404 空气流量系数**p62p15duct 外涵出口总压**T62T15 外涵出口总温

Wg62W15Wg225Wg125 外涵出口总流量q(62)(12)62（1-11-11） 气动函数-1262162q1q(62) 速度系数

*p62p62(62) 外涵出口静压p62p61Wg61Wg62Wg6

WhWhWh

g6262g66g6161p*f()AP*f()Ap*f()(AA)6616261

*p6Wg6kmA6q(6)

*T6据此方程组求得后混合器出口总温、总压、流量（6）、展开发动机平衡方程六

0

（5.2.1-6）

A8A8加力燃烧室（进口总温、总压、流量采用后混合器出口总温、总压、流量）：P6*outP6*in6**T6outT6in

Wg6outWg6in- 17 -

尾喷管运算：

A89.5544e+003 尾喷管面积

A8Wg8T8*kmpq(8)*8 由尾喷管流量公式计算得尾喷管面积18212Ma2 速度系数

Ma212P8*P6*out**T8T6out 尾喷管进口总温、总压、流量

WWg6outg81.33,R287.31 燃气气体绝热指数、气体常数km0.0397 燃气流量系数（7）、展开发动机平衡方程七

Wa2Wa21Wa130

（5.2.1-7）**Tinkwpout,dCwW(nH,ZCDFS,CDFS)(1CDFS)**1 CDFS流量

100Tout1pin,dWa21WCDFSWa13Wg225km*p225T*225A225q(225) 前混合器副外涵出口流量

详见前混合器运算**Tinkwp1,dWa2WCLCwW(nL,ZCL,L)(1L)* 风扇流量

*100T1pin,d对于第二问所给条件，本方程组7个一共有7个待定变量，他们分别是ZCL、ZCDFS、ZCH、T4、ZTL、ZTH、nH。

5.2.2 遗传算法模型的建立与求解

由于发动机的7个平衡方程组成的非线性方程组，利用传统的非线性方程组求解，无论是从算法的选择还是算法本身的构造都与索要解决的问题的特性有很大的关系，很多情况下算法中算子的构造及有效性会成为我们解决问题的巨大障碍。而遗传算法作为一种灵活的自适应算法无需过多的考虑问题的具体- 18 -

*

形式，具备全局收敛性，突破了传统算法的诸多限制与不足，因而本文采用遗传算法进行发动机非线性方程组的求解。

（1）、非线性方程组求解转化成遗传算法的优化问题

设有限空间内的非线性方程组：

f1(x1,x2,......,xn)0f(x,x,......,x)0212n

........fn(x1,x2,......,xn)0其中：aixibi , i=,n ;

常量ai、bi为相应自变量xi的上下限。

记Xx1,x2,......,xn，上述方程组可表示为：

f1(X)0f(X)02

........fn(X)0构造非线性函数

F(X)f1(X)f2(X)fn(X)

XD, D(a1,b1)(a2,b2)(an,bn)为n维欧式空间的一个有界区域。***在整个自变量的有界区域，当上述非线性方程组有解X*函x1,x1,......,x1时，数满足F(X)0；即，函数的极小值为0

minF(X)F(X*)0

很显然上述非线性方程组的求解问题就等价成了函数F(X)的极小值问题，从而，通过构造的函数F(X)，将非线性方程组的求解问题转换成了遗传算法优化问题[6]函数F(X)取极小值时的X即为所求方程组的解。

（2）、遗传算法基本步骤

第一步：采用二进制编码，由变量xi，要求的精度i，计算变量的编码长度li，并确定遗传算法个体编码长度Lli。

i1n第二步：初始化

①确定群体规模N、杂交概率pc、变异概率pm及终止进化准则;

②随机生成初始群体Y(0)Y1(0),Y2(0),,YN(0), (二进制字符串群);

③计算Yi(0)的适应度fYi(0), (i=1,2,,N);

④置进化代数k0;

第三步：遗传操作

①对个体Yi(k)，依据其适应度fYi(k)，计算复制概率

- 19 -

pi(k)fYi(k)/fYi(k)

i1N②以概率pi(k)从Y(k)中选择个体，并保留最佳个体，形成新群体Y(k1)：

(k1)

Y(k1)Y1(k1),Y2(k1),,YN③在新群体Y(k1)中，依次选取个体Yi(k1)，Yi1(k1)以概率pc对Yi(k1)，Yi1(k1)进行一点交叉，产生两个新个体Yi (k1)，Yi1(k1)。循环执行本过程N/2次，形成新群体Y(k1)：

(k1)

Y(k1)Y1 (k1),Y2(k1),,YN④在新群体Y(k1)中，依次选取个体Yi (k1)，以概率pm对Yi (k1)

进行简单变异，产生新个体Yi (k1)。循环执行本过程N次，形成新一

代群体Y(k1)：

Y(k1)Y1(k1),Y2(k1),,YN(k1)

第四步：如果Y(k1)满足进化停止准则，解码计算输出最优解X(t+1)，否则

k=k+1转第三步。

遗传算法的流程图如下：

图5-3 遗传算法的流程图

- 20 -

（3）、遗传算法关键步骤的解释

①对参数进行二进制编码

遗传算法求解问题时，不是直接在问题的解空间上操作，而是利用解的某种编码讲解空间的数据表示成遗传空间的基因型串结构数据，并进行运算。最常用的编码方式是二进制0/1字符编码。二进制编码将解空间映射到0/1空间上，然后在位串空间上进行遗传操作操作结果再通过解码还原成变量原型，以进行适应性评估。

②初始群体是遗传算法搜索的出发点

对于群体规模N，编码长度L的问题，随机产生N*L位编码字符，组成N个初始个体，构成初始群体。

③适应函数设计

本文采用的适应函数为：

U(X)FmaxF(X)k(FmaxFmin)

其中，k(k0)为当前代群体中最大与最小适应度（Umax 、Umin）之比的控制系数；Fmax、Fmin分别为当前代群体中函数F(X)的最大和最小函数值。

④遗传算子设计

遗传算子包括选择算子、交叉算子、变异算子。

选择算子：优胜劣汰的选择机制使得适应度大的个体有较高的存活率，本文的选择算子采用轮盘式的选择方法，该法首先计算每个个体的相对适应值作为对应个体的复制概率，然后根据复制概率将一个圆盘分成N份，对应扇形中心角的大小与其复制概率成正比，选择个体时，圆盘内指针不动，转动圆盘，指针落入第i个扇形内，则选择第i个个体。从统计意义上来说，个体相对适应度越大，则被选择到的机会越多，选择操作无法产生新个体。

交叉算子：交叉操作仿造生物学中的杂交原理，采取一定方法使个体间一个或多个基因交叉，产生新个体，是遗传算法中产生新个体的主要手段。父本杂交比率由交叉概率决定。

变异算子：选择与交叉无法完成在初始群体组合以外的空间进行搜索，使进化过程容易过早陷入局部解而终止进化过程，通过变异操作可确保群体中基因的多样性，以使搜索空间能够尽可能大。

⑤算法控制参数

遗传算法中需要选择的运行参数主要有个体编码串长度、群体规模、交叉概率、变异概率，这些参数对遗传算法的性能影响较大，须认真选取。

⑥算法停止准则

由于目标函数的极小值一定且为0，因此可根据方程组解的精度要求，明确给出如下进化终止准则（收敛判断准则）

- 21 -

f(Xii1n2*)

***其中，X*x1,x1,......,x1为当前进化代最优个体对应解；

；

为方程组要求的精度（收敛误差）5.2.3 遗传算法的有效性论证

遗传算法具有大范围收敛性质，对初始点没有任何要求与限制，而且只要有解的存在，从理论上证明就可以完全收敛于全局最优解；采用群体搜索，具有良好的并行性；遗传算法是一种求解系统优化问题的通用方法，不依赖于所研究问题的具体领域，运算通过种群更新和迭代来搜索全局最优解，与其他优化算法相比，对复杂工程问题的求解具有较强的鲁棒性；仅利用目标函数取值信息，无需诸如梯度及其他辅助信息，适合于大规模、高度非线性及无解析表达式的目标函数优化问题；采用随机化算子而不是严格的确定性运算，可以直接逼近非线性、多约束、多目标优化问题的求解目标；不要求组成方程的函数连续、可导，可求解常规算法难以求解的超越方程组。这些特点都证明了遗传算法的有效性。

5.2.4 本题的计算

本题共涉及到7个非线性方程组成的方程组，含有7个未知量：

nH、ZCL、ZCDFS、ZCH、ZTH、ZTL、T4*

利用MATLAB7.0编码进行计算上述非线性方程组得到7个未知量的值为：

nH=0.87；ZCL=0.70；ZCDFS=0.79、ZCH=0.31、ZTH=0.52、ZTL=0.22、T4*=1419.9

5.3 求解问题三

考察发动机平衡方程函数关系，在一定的飞行条件下(高度与飞行速度一定，即H，Ma为定值)，当给定发动机部件特性时，即描述了发动机的工作过程，但由上述方程唯一确定满足发动机共同工作方程的解，还需给出发动机的调节规律。

对于本文所研究的双涵道变循环发动机，单涵道模式时，可调参数有低压转数、高压转数、风扇压比函数值、CDFS压比函数值、高压压气机压比函数值、主燃烧室出口温度、高压涡轮压比函数值、低压涡轮压比函数值、风扇可调导叶角度、CDFS可调导叶角度、高压压气机可调导叶角度、高压涡轮可调导叶角度、低压涡轮可调导叶角度、喷管喉部面积 共14个参数，本题中风扇导叶角度、高压压气机导叶角度、高压涡轮导叶角度均设置为 ，故还有11个未知参数。上述描述发动机工作过程的一系列函数的参数，并非相互独立，对于发动机非设计状态，在所有未知变量中，一般选取7个独立的变量作为未知量，另外4个设为合理常数初值。参考已有文献制定发动机4个赋予初值的参数方案如下：

T4*Const nLConst ZCHConst ZTHConst

- 22 -

对应地，共同工作方程可表示为所选变量隐形式的(非显函数)非线性方程组，即：

 f1(x1,x2,x3,x4,x5,x6,x7)0

f(x,x,x,x,x,x,x)021234567f3(x1,x2,x3,x4,x5,x6,x7)0f4(x1,x2,x3,x4,x5,x6,x7)0

f(x,x,x,x,x,x,x)051234567f6(x1,x2,x3,x4,x5,x6,x7)0f7(x1,x2,x3,x4,x5,x6,x7)0发动机共同工作点求解问题相应转化为，在给定调节规律下的上述非线性方程组求解问题。由于上述非线性方程组无明显的数学表达形式，因此，不可能得到发动机共同工作点的解析解，在实际发动机非线性数学模型中，通常采用数值计算方法，首先给出一组独立变量的初值，代入发动机模型，计算各平衡方程的误差，通过一定的非线性方程组迭代解法，不断地修正独立变量的试给值，直至各平衡方程的误差满足设定的迭代精度要求。

六、模型评价和优化

6.1 模型的评价

6.1.1 优点

本文第二问采用了遗传算法模型，遗传算法模型的优点如：

（1）、具有大范围收敛性质，对初始点没有任何要求与限制，而且只要有解的存在，都能求得方程组的解。

（2）、对优化的目标函数无限制，不需要函数连续，也不要求可微，既适用于有数学解析表达的显函数，也适用于只有输入输出映射的隐函数，应用范围广。

（3）、采用群体搜索，从多个初始点可行解出发，且算法本身隐含成倍模式，具有良好的计算并行性。

（4）、遗传算法是一种求解系统优化问题的通用方法，不依赖于所研究问题的具体领域，运算通过种群更新和迭代来搜索全局最优解，与其他优化算法相比，对复杂工程问题的求解具有较强的鲁棒性。

（5）、仅利用目标函数取值信息，无需诸如梯度及其他辅助信息，适合于大规模、高度非线性及无解析表达式的目标函数优化问题。

（6）、采用随机化算子而不是严格的确定性运算，可以直接逼近非线性、多约束、多目标优化问题的求解目标。

6.1.2 缺点

遗传算法模型的计算效率还不太理想、算法参数选择不合理时有陷入局部- 23 -

最优解的可能。

6.2 模型的优化

由于遗传操作的随机性，使得新一代的个体产生没有明确的方向，个体编码的基因组合不能保证向最优个体靠近，这大大增加了迭代优化时间，降低了算法效率。为了加快优化到最优解的速度和精度，可以引入“进化方向”算子，在每一代的进化中，该算子的作用是指导个体向更优的个体编码方向进化。这样引入“进化方向”算子可以大大增加计算效率。

- 24 -

七、参考文献

[1] 苟学中，周文祥，黄金泉．变循环发动机部件级建模技术．航空动力学报，

2013，28（1）：104-111．

[2] 王 元，李秋红，黄向华．变循环发动机建模技术研究．航空动力学报，2013，

28（4）：954-960．

[3] 方昌德．变循环发动机．燃气涡轮试验与研究，2004，17(3)：1-5．

[4] 刘增文，王占学，黄红超，蔡元虎．变循环发动机性能数值模拟．航空动

力学报，2010，25（6）：1310-1315．

[5] 苟学中．变循环发动机建模及控制规律研究．南京航空航天大学硕士学位

论文，2011年12月．

[6] 张宏生，吴拓．遗传算法在非线性方程组求解中的应用．肇庆学院学报，

2002，23(2)：16-19．

[7] 苏三买．遗传算法及其在航空发动机非线性数学模型中的应用研究．西北

工业大学博士学位论文，2002年4月．

[8] 谭浩强．C程序设计（第二版），北京：清华大学出版社，1999年．

- 25 -

八、附录

附表一 风扇特性数据表中流量随压比函数值变化对应值表

换算转速 0.4

流量

zz

0 38.42672

0.111741 37.62523

0.215547 36.81818

0.310931 36.00186

0.39919 35.18553

0.480486 34.3692

0.55498 33.55288

0.622996 32.73655

0.68664 31.93692

0.739919 31.10204

0.787045 30.26716

0.830769 29.45083

0.868988 28.63451

0.901862 27.81818

0.927126 26.9833

0.949798 26.16698

0.970364 25.3692

0.986559 24.57328

0.99498 23.75139

1 22.93506

换算转速 0.7

流量

zz

0 57.19109

0.11908 56.83488

0.230352 56.40445

0.336563 55.94434

0.432651 55.39332

0.520859 54.78664

0.600246 54.11874

0.673126 53.41929

0.736678 52.65863

0.792857 51.86085

0.841371 51.02597

0.882293 50.15399

0.917215 49.26345

0.946208 48.35436

换算转速 0.5

流量

zz

0 43.49722

0.117736 42.97959

0.230252 42.44898

0.338476 41.91837

0.434868 41.30612

0.522561 40.65306

0.60573 40

0.675792 39.26531

0.740865 38.53061

0.792367 37.71429

0.834358 36.85714

0.869157 35.97959

0.89746 35.08905

0.921007 34.19852

0.939914 33.30798

0.956502 32.43599

0.968797 31.56401

0.979237 30.71058

0.99014 29.89425

1 29.10019

换算转速 0.81

流量

zz

0 68.57143

0.107614 68.27273

0.211084 67.93043

0.310506 67.55102

0.404723 67.12059

0.494024 66.6475

0.577928 66.12987

0.655904 65.56215

0.72641 64.92764

0.78906 64.22635

0.843614 63.45826

0.891518 62.64564

0.929976 61.7551

0.958651 60.78664

- 26 -

换算转速 0.6

流量

zz

0 49.33117

0.119572 48.97866

0.230852 48.55195

0.335622 48.08163

0.432011 47.55102

0.525011 47.02041

0.613107 46.4746

0.691841 45.86456

0.758181 45.1628

0.814356 44.39889

0.859296 43.56401

0.891485 42.64564

0.918145 41.718

0.940348 40.79035

0.958181 39.86271

0.971823 38.93443

0.983504 38.02531

0.991262 37.1162

0.997414 36.22635

1 35.33581

换算转速 0.9

流量

zz

0 82.98701

0.111511 82.94991

0.221189 82.87941

0.327264 82.74583

0.427935 82.52319

0.522646 82.2115

0.612969 81.84601

0.694744 81.36178

0.768495 80.77737

0.83478 80.10946

0.895431 79.39147

0.943147 78.51577

0.978975 77.50649

0.99738 76.28942

0.96544 47.38776

0.97925 46.40816

0.989082 45.42857

0.995011 44.44898

0.999277 43.4898

1 42.53061

换算转速 0.95

流量

zz

0 91.63451

0.104765 91.63451

0.208261 91.60668

0.310228 91.54731

0.410562 91.45826

0.508357 91.32468

0.601073 91.10204

0.688321 90.79035

0.767718 90.35436

0.837138 89.76623

0.894976 89.00928

0.939856 88.0705

0.972532 86.97217

0.991423 85.69573

1 84.30427

0.999326 82.81737

0.987873 81.21194

0.973025 79.60651

0.956389 78.0249

0.936151 76.4361

0.979904 59.76994

0.993205 58.7013

0.998169 57.57488

0.999855 56.45269

1 55.35251

0.997446 54.25945

换算转速 1

流量

zz

0 101

0.093383 101

0.186786 101

0.280189 101

0.373572 101

0.466975 101

0.560378 101

0.653761 101

0.743448 100.9202

0.827128 100.7199

0.897091 100.2523

0.957745 99.62523

0.993854 98.55288

1 96.95362

0.987094 95.05395

0.966679 93.07813

0.942508 91.09577

0.915175 89.11563

0.884719 87.1355

0.851279 85.15537

1 74.89425

0.99476 73.4397

0.983724 71.95547

0.969117 70.47124

0.953201 69.0167

0.938235 67.62152

换算转速 1.075

流量

zz

0 108.5

0.065917 108.5

0.131834 108.5

0.197751 108.5

0.263669 108.5

0.329598 108.5

0.395515 108.5

0.461433 108.5

0.52735 108.5

0.593267 108.5

0.656739 108.4267

0.714963 108.2041

0.765964 107.7885

0.81086 107.2245

0.850185 106.5362

0.88596 105.7848

0.917752 104.9629

0.947845 104.1317

0.975033 103.2616

1 102.3711

- 27 -