基于精确偏导数的航空发动机部件级模型构建方法与流程

本发明属于航空宇航推进理论与工程中的系统仿真与控制领域，具体涉及一种基于精确偏导数的航空发动机部件级模型构建方法。

背景技术

在航空发动机的研究领域，数学模型发挥着重大作用。而在各种数学模型中，部件级模型凭借其能够大包线使用、能估计发动机截面参数和推力这类不可测量而始终受到青睐。

然而，部件级模型因为其共同方程组不封闭，需要应用数值方法进行求解。其中最为常用的为牛顿法和拟牛顿法，这两类方法中由于迭代的需要，都涉及到计算有关函数对自变量的雅克比矩阵(即偏导数)。这也就意味着部件级模型为了能够实现仿真计算，就需要在每一迭代步内大量地重复计算偏导数。

而长期以来，在航空发动机模型中都没有考虑直接对模型所需的雅克比矩阵进行求解，这也就造成了部件级模型中变量偏导数直接计算方面的空白。再者，由于航空发动机的热力参数相互影响，其某一变量对另一变量的影响不易厘清，加之模型内部存在着局部迭代、查表等，常规的求导方法难以适用，因而在各个部件的热力学模型建立后，直接利用所建立的部件模型结合差分方法成了近似求解相应偏导数关系的通用方法，也即基于差分方法的部件级模型。

但这种基于差分的部件级模型存在较为明显的缺陷，即该模型每次为了更新发动机状态，都需要对变量进行小扰动再重复调用部件模型多次以完成仿真计算。具体而言，一是应用差分方法时，变量扰动幅值的最优值是难以选择的。因为航空发动机强烈的非线性特征，在不同飞行包线内，变量的扰动量过小可能导致因变量没有明显变化，进而导致其导数为零。或是扰动幅值过大，虽可以使得因变量有较为明显变化，却使发动机模型的当前工作状态也发生明显变化，以致偏离当前实际共同工作点而难以实现近似，进而直接影响了导数的精确性。二是应用差分方法时需要多次调用模型进行计算以获得近似值，因变量越多则需要重复调用模型的次数越多。这也就意味着，在该模型的气动热力计算过程中，计算量因为雅克比矩阵的需求而被显著地扩大，进而就使得实时性受到影响。

为了提高导数精度和模型的收敛性能，在基于差分的部件级模型中通常采用中心差分的方法，使得为了计算某一对因变量和自变量导数时，就需要上下摄动变量然后调用模型(即需要额外调用模型2次)进行相应的气动热力计算。以常见的双轴混排涡扇发动机为例，该发动机模型在计算中通常涉及6个初猜值和6个平衡方程，由于中心差分的使用，就需要额外调用模型12次。换句话说，在模型计算过程中，除了自身必须的热力计算需要调用模型进行1次计算外，由于偏导数的需求而额外增加了12倍的计算量。

技术实现要素：

本发明所要解决的技术问题在于克服现有技术不足，提供一种基于精确偏导数的航空发动机部件级模型构建方法，利用该模型，在计算获得各部件截面参数的同时既可以保证更高精度的导数值，又可以避免模型的重复调用有利于降低模型的总体计算量，从而提高实时性。

本发明基于精确偏导数的航空发动机部件级模型构建方法，包括以下步骤：

步骤a.建立航空发动机各部件的气动热力学模型和各部件共同工作方程；

步骤b.分析共同工作方程所需偏导数的因变量和自变量，结合各部件气动热力学模型和链式法则，建立相应的部件级偏导数模型；具体包括：

步骤b1.分析共同工作方程所需要求解的偏导数对应的因变量和自变量及其所属部件；

步骤b2.记共同工作方程中所需的自变量为v，且v∈r^l,l≥1，因变量为f，且f(u)∈rⁿ,n≥1，作为中间变量的所有热力学参数，记为u∈r^m。

将自变量v的微分向量转化为对角化矩阵

式中，dv为v的微分向量，dw为dv对角化后的矩阵；

步骤b3.根据链式法则，建立中间变量的热力学参数u对自变量v的微分关系和因变量f对热力参数u的微分关系，即：

步骤b4.获得因变量f对自变量v的部件级偏导数模型：

步骤c.联合各部件气动热力学模型和部件级偏导数模型共同计算，结合迭代算法，实现部件级模型的稳态计算或动态计算。

进一步地，所述步骤c具体包括：

步骤c1.初始化模型；

步骤c2.给定k时刻的飞行条件和发动机输入；

步骤c3.开展各部件模型和各部件偏导数模型的联合计算，并输出雅克比矩阵计算需要的精确导数值；

步骤c4.判断是否满足收敛条件，若满足，则执行步骤c6，否则继续执行步骤c5；收敛条件是指平衡方程残差是否小于给定的阈值要求，若小于，则满足，否则不满足；

步骤c5.利用获得的偏导数，构造雅克比矩阵，更新部件级模型的初猜值和各截面热力参数，并判断是否达到最大允许的迭代次数，若达到则执行步骤c6，否则返回步骤c3；

步骤c6.动态计算时则进行转子动力学计算，更新发动机转速，并判断动态计算过程是否结束，并输出各个部件截面的热力参数；稳态计算过程不进行转子动力学计算，直接输出各个部件截面的热力参数；

步骤c7.判断模型计算仿真是否结束，若否则返回步骤c2，若是则计算结束。

本发明可采用现有的各种迭代算法，优选地，步骤c中的迭代算法为牛顿法或者拟牛顿法，例如newton-raphson方法等。

自变量的微分初始化值可以取任意常数，优选地，步骤c1中对部件级偏导数模型进行初始化时，各自变量的微分初始化值均为1。

优选地，在步骤b3中，沿着航空发动机气流流程顺序逐部件建立所述微分关系。这样，对于和多个部件相关的微分计算，可以将这些微分关系逐部件递推计算，构成各部件的偏导数模型。

相比现有技术，本发明技术方案具有以下有益效果：

(1)精度高：本发明能够直接计算得到具有更高精度的雅克比矩阵，能更准确地更新猜值，进而更精确地更新发动机状态，有效避免了差分方法需要扰动变量和扰动幅值选择困难的缺陷，也避免了由于导数的近似求解所产生的不利影响。

(2)实时性好且计算简单：本发明通过部件热力计算和偏导数模型的同步执行直接求解得到精确雅克比矩阵，且均沿着气流流程方向进行各部件计算，计算流程简单，又避免了差分方法需要重复调用部件模型的缺点，从而大大减少模型的计算量，提高实时性。

(3)通用性和可移植性：本发明通过结合发动机工作原理、气动热力学理论和部件级模型计算方法实施，对由多个部件组成并共同工作的各类航空发动机机均适用。

附图说明

图1为典型的双轴混排涡扇发动机结构图；

图2为基于精确偏导数的部件级模型稳态计算流程示意图；

图3为基于精确偏导数的部件级模型动态计算流程示意图；

图4为基于差分的部件级模型稳态计算流程示意图；

图5为基于差分的部件级模型动态计算流程示意图；

图6有效性和收敛性测试稳态仿真收敛图。

具体实施方式

针对现有基于差分方法的部件级模型所存在的缺陷，本发明通过结合相应的发动机工作原理、气动热力学理论、转子动力学和相应的部件级模型计算方法等，建立各部件气动热力学模型，通过将实际的求导转化为求微分，结合链式法则建立相应的偏导数模型，联合部件气动热力学模型和偏导数模型进行共同计算，快速构造雅克比矩阵，进而在线实时求解部件级模型。其适用对象为由各个部件组合并共同工作的航空发动机，包括但不限于涡喷发动机、涡扇发动机、涡轴发动机、涡桨发动机、变循环发动机、涡轮基冲压组合发动机等。

为了便于公众理解，下面以双轴混排涡扇发动机为例来对本发明技术方案进行详细说明。

图1所示为典型的双轴混排涡扇发动机结构图。图中1截面为进气道进口；2截面为进气道出口和风扇进口；25截面为内外涵进口；3截面为压气机出口和燃烧室进口；4截面为燃烧室出口；41截面为高压涡轮进口；42截面为高压涡轮出口；45截面为低压涡轮进口；46截面为低压涡轮出口；16和6截面分别为外涵出口和内涵出口；7截面为加力燃烧室进口；75截面为加力燃烧室出口；8截面为尾喷管喉道截面；9截面为尾喷管出口。

首先根据步骤a，对各个部件和共同工作关系进行建模。

建立各部件共同工作方程，即根据流量平衡、功率平衡、压力平衡来建立6个稳态计算时的共同工作方程：

1)风扇和低压涡轮同轴，应当满足功率平衡：

f1＝wf/(wlt·δlt)-1(4)

式中，wf为风扇功率，wlt为低压级涡轮发出的功率，δlt为低压轴机械效率。

2)压气机和高压涡轮同轴，应当满足功率平衡：

f2＝(wc+wex)/(wht·δht)-1(5)

式中，wc为压气机功率，wex为涡轮抽功量，wht为高压涡轮发出的功率，δht为高压轴机械效率。

3)根据高压涡轮转子进口流量平衡，有平衡方程：

f3＝(mx41c-m41c)/m41c(6)

式中，m41c为高压涡轮进口实际流量的相似流量。

4)根据低压级涡轮转子进口流量平衡，有平衡方程：

f4＝(mx45c-mx45c)/m45c(7)

式中，m45c为低压级涡轮进口实际流量的相似流量。

5)根据气流从内外涵道进入掺混室掺混，有压力平衡方程：

f5＝(ps16-ps6)/ps6(8)

式中，ps16为外涵道出口的静压，ps6为内涵道出口的静压。

6)根据尾喷管喉道处总压平衡原则，有平衡方程：

式中，为尾喷管喉道处总压，为喷口背压。

对于上述6个方程构成的方程组，通过newton-raphson方法进行迭代，当6个残差f1到f6小于某一阈值时，认为模型此时收敛。

在动态过程中，由于功率不再平衡，需要考虑转子动力学，因此在平衡方程组中去掉公式(4)(5)两式，并引入转子动力学方程。此时，同样采取采用newton-raphson方法，平衡方程组为公式(6)-(9)。

根据发动机工作原理，确定稳态计算时初猜值为低压涡轮与风扇轴相对物理转速nf、高压涡轮与高压压气机轴相对物理转速nc、风扇压比系数zf、压气机压比系数zc、高低压涡轮进口相似流量mx41c和mx45c，动态计算时初猜值为风扇压比系数zf、压气机压比系数zc、高低压涡轮进口相似流量mx41c和mx45c。

同时，建立各个部件的气动热力学模型。

考虑到后文为了描述方便，将以低压涡轮部件中的低压转子功率平衡方程f1对高压涡轮中的猜值mx41c的导数为例来说明步骤b的过程。因此在此特别列出高压涡轮热力学模型am1、高低压涡轮部件间涵道热力学模型am2和低压涡轮热力学模型am3。其中高压涡轮气动热力学模型am1的输入为高压涡轮进口流量m41、进口焓h41和进口总压输出为高压涡轮出口流量m42、出口焓h42和出口总压高低压涡轮间涵道气动热力学模型am2的进口参数即为高压涡轮出口参数和冷却气流量mcool、冷却气焓hcool，其出口参数即为低压涡轮进口参数，包括低压涡轮进口流量m45、进口焓h45和进口总压低压涡轮气动热力学模型am3出口参数包括出口流量m46、出口焓h46和出口总压

由于高低压涡轮的气动热力计算过程基本相同，故这里仅给出高低压涡轮间涵道气动热力学模型am2和低压涡轮热力学模型am3。

高低压涡轮之间涵道的热力计算模型am2：

考虑涡轮进口冷却气流掺混：

m45＝m42+mcool(10)

f45＝mfb/(m45-mfb)(11)

h45＝(m42h42+mcoolhcool)/m45(12)

式中，m45表示低压涡轮进口实际流量，m42表示高压涡轮出口燃气流量，mcool表示冷却气实际流量，f45表示油气比，h45、h42和hcool表示低压涡轮进口燃气焓、高压涡轮出口燃气焓和冷却气流焓，表示高压压涡轮出口总压和低压涡轮进口总压。

低压涡轮的气动热力学模型am3：

已知低压涡轮焓降dhlt和效率ηlt与低压涡轮折合转速低压涡轮进口折合流量m45c之间的特性关系，猜值mx45c以及相应的设计点参数。

步骤am3.1：根据焓熵总温和油气比关系求出低压涡轮进口总温和熵s45：

s45＝g2(h45,f45)(15)

式中，g1和g2表示相应的焓、熵、总温和油气比关系。

步骤am3.2：计算低压涡轮进口折合流量m45c和低压涡轮进口折合转速

式中，和分别表示低压涡轮进口设计点总压和总温。

步骤am3.3：根据涡轮特性关系求得焓降dhlt和效率ηlt：

式中，g3和g4表示相应的特性关系。

步骤am3.4：计算低压涡轮出口焓h46和理想焓h46i：

h46i＝h45-(h45-h46)/ηlt(21)

步骤am3.5：由焓熵关系g2和理想焓h46i求出理想熵s46i：

s46i＝g2(h46i,f45)(22)

步骤am3.6：由熵增计算出低压涡轮落压比πlt和低压涡轮出口总压以及出口流量m46和油气比f46：

m46＝m45(25)

f46＝f45(26)

步骤am3.7：求出低压涡轮功wlt：

wlt＝m45(h45-h46)(27)

步骤am3.8：计算低压涡轮功率平衡方程f1：

f1＝wf/(wlt·δlt)-1(28)

下面将以低压涡轮部件中的低压转子功率平衡方程f1对高压涡轮中的猜值mx41c的导数为例来说明步骤b的过程。低压涡轮部件对其他初猜值的偏导数模型以及其他部件的偏导数模型都可以按照参照同样方法建立。值得指出的是f1对mx41c的导数需要结合高压涡轮和低压涡轮部件模型的计算过程获得。

根据步骤b1，mx41c为高压涡轮部件内部参数，故属于高压涡轮部件，f1为低压涡轮部件内部参数，故属于低压涡轮。因此沿气流流程方向逐步建立从高压涡轮部件、低压涡轮部件以及高低压涡轮之间涵道的各部件偏导数模型。

根据步骤b2，设猜值mx41c的微分为dmx41c。特别指出，此处仅仅是以mx41c为例，而实际上还存在其他猜值，例如风扇压比比zf的微分dzf等，并将这些微分构成自变量的微分向量用以建立后续的偏导数模型。

根据步骤b3，此处以dmx41c为例，需要建立mx41c到f1的微分关系。根据链式法则和中间变量的分析，则需要构建mx41c到高压涡轮出口焓h42的微分关系式bm1、h42到低压涡轮进口焓h45的微分关系式bm2以及h45到f1的微分关系式bm3。因为这些关系式与各部件热力计算紧密相关，故可以按照部件顺序逐个构建。该三个微分关系式所要描述的内容如下：

由于高低压涡轮的热力学模型基本相同，因此其二者建立的微分关系式也基本相同。因此，下面仅给出微分关系式bm2和bm3。

h42到低压涡轮进口焓h45的微分关系式bm2：

特别指出，这里仅仅给出了h42到h45的微分关系式为例子。但事实上，还存在其他进口变量到出口变量的微分关系式，例如求解f6到zc的偏导数时，需要以高低压涡轮间涵道的进出口热力参数为中间变量建立联系等。这些与高低压涡轮间涵道热力计算紧密相关的微分关系式共同构成了高低压涡轮间涵道偏导数模型bm4。

低压涡轮进口焓h45到f1的微分关系建立过程bm3：

步骤bm3.1：根据步骤am3.1求微分，由此得到对h45的偏微分：

步骤bm3.2：根据步骤am3.2求微分，得到对的微分

步骤bm3.3：根据步骤am3.3求微分，求得焓降dhlt对的偏微分ddhlt：

步骤bm3.4：根据步骤am3.4求出h46对h45的偏微分dh46：

步骤bm3.5：根据步骤am3.7求得低压涡轮功的微分dwlt：

dwlt＝m45(dh45-dh46)(35)

步骤bm3.6：根据步骤am3.8求微分获得残差对h45的导数：

由此，即建立了f1对h45的微分关系式。特别指出，此处仅以h45为例，实际上还存在其他与低压涡轮热力计算相关的微分运算，例如求解f1对zf偏导数时，需要以低压涡轮进口的热力参数h45、m45来传递zf到f1的影响，再如求解f6到zc的偏导数时，需要以低压涡轮的进出口热力参数为中间变量建立联系等。这些与低压涡轮相关的计算共同构成了低压涡轮部件偏导数模型bm5。

此外，由上述过程可知，式(10)-(28)和式(30)-(36)可以同步进行计算，在计算的过程中式(30)-(36)可以大量复用式(10)-(28)中的热力参数。而且，上述过程中，各式都可以沿着发动机气流流程方向从高压涡轮出口起计算至低压涡轮出口。

根据步骤b4，将微分关系式bm1、bm2和bm3利用链式法则串联起来，即得到f1对mx41c的偏导数模型bm6：

类似地，也可以获得其他平衡方程对其他变量的偏导数模型。优选地，取各个自变量微分初始值为1。此外，由于各平衡方程对猜值的偏导数模型中，都涉及到与各个部件的热力计算相关的部分，因此可以仿照bm4和bm5得到各个部件偏导数模型，进而构成了完整的部件级偏导数模型。

由此可以建立稳态迭代时6个猜值、6个平衡方程建立的完整偏导数模型和动态迭代时4个猜值、4个平衡方程建立的完整偏导数模型。这两个模型均可以沿着航发的气流方向，从发动机进口计算至发动机出口，从而有效地避免了传统求导中首先需要完成发动机模型所有热力参数的计算，再从因变量反向计算回自变量的缺陷，亦不需要厘清各个变量之间的影响关系，易于实现。又避免了差分方法所要求的对变量进行扰动和重复调用模型的缺点。

之后实现相应的稳态和动态计算过程。

对于稳态计算，其过程如下：

执行步骤c1.初始化模型；

执行步骤c2.给定稳态飞行条件和发动机输入；

执行步骤c3.优选地取各自变量微分值为1；同步执行部件模型计算过程和部件级偏导数模型，并输出平衡方程组对初猜值的导数；

执行步骤c4.判断模型是否收敛，如收敛，则转步骤c6；如不收敛，继续执行步骤c5；

执行步骤c5.利用获得的偏导数，构造雅克比矩阵，通过相应的迭代算法更新部件级模型的猜值和各部件截面参数；判断是否达到最大允许迭代次数，若没有达到则转步骤c3；

执行步骤c6.输出各个部件截面的热力参数；

执行步骤c7.判断模型计算仿真是否结束，若否则返回步骤c2，若是则计算结束。

对于动态计算，其过程如下：

执行步骤c1.以上一时刻计算结果为初始状态；

执行步骤c2.给定k时刻的飞行条件和发动机输入；

执行步骤c3.优选地取各自变量微分值为1；同步执行部件模型计算过程和部件级偏导数模型，并输出平衡方程组对猜值的导数；

执行步骤c4.判断模型是否收敛，如收敛，则转步骤c6；如不收敛，继续执行步骤c5；

执行步骤c5.利用获得的偏导数，构造雅克比矩阵，通过相应的迭代算法更新部件级模型的猜值和各部件截面参数；判断是否达到最大允许迭代次数，若没有达到则转步骤c3；

执行步骤c6.进行转子动力学计算，并输出各个部件截面的热力参数；

执行步骤c7.判断模型仿真是否结束，若未结束则返回步骤c2，若是则计算结束。

相应地，稳态和动态的计算也发生了变化。图2和图3给出了本发明的基于精确导数的部件级模型计算过程，同时为了便于比较图4和图5给出了传统基于差分的部件级模型进行稳态和动态计算时的流程图。从图中可以看出，传统的模型计算中由于差分方法的应用，需要额外调用多次部件模型进行计算，而应用本发明后，部件气动热力计算和相应的偏导数求解同步进行，不再需要重复调用模型，有利于提高实时性。

为了验证本发明的有效性和收敛性，进行不同飞行条件不同收敛阈值条件下的稳态和动态计算仿真。仿真的飞行条件如下表1所示。

表1仿真飞行条件

仿真前模型初始化在地面点(表1中编号0)，之后逐渐改变飞行条件来进行仿真稳态计算过程。仿真是从初始点起按照表格中编号顺序依次改变飞行条件，且每次改变飞行条件后需保证模型计算收敛(即达到指定的收敛精度)。采用newton-raphson算法以修正初猜值，迭代步长设置为0.1。测试的收敛精度范围为10^-6到10^-14。仿真结果如图6所示，图中横坐标为所设置的收敛精度的以10为底的对数绝对值，纵坐标为模型的平衡方程收敛到该精度所需要的气动热力计算次数。从图6中可以看出，在不同的飞行条件下，气动热力计算次数随收敛精度的对数线性增长，表明在不同飞行条件下本发明能够稳定收敛到所需精度，且保持了牛顿-拉夫逊算法原有的超线性收敛性，证明了本发明的有效性和收敛性。

为了验证本发明的实时性优势，按照表1中所示的飞行条件进行稳态和动态计算仿真。牛顿-拉夫逊算法的迭代步长设置为0.1，收敛精度设置为10^-6。稳态计算的最大允许迭代次数为1000次，动态计算的最大允许迭代次数为1次，即采用一次通过算法。为了便于对于效果，对采用中心差分方法进行模型求解的部件级模型进行相同的仿真。重复仿真50次，记录相关数据并取平均值，展示在下表2-表5中。本仿真为仿真环境为dellt5810win7旗舰版，cpu为intelxeon(tm)1650v43.6ghz，32gb内存。程序运行平台为vs2010旗舰版release模式。计时函数采用c++中的windows.h标准库提供的计时函数。

稳态仿真的结果如下表2所示，动态仿真结果如表4所示。表中飞行条件编号与表1中相对应，模型类型中模型1指本发明所建立的基于精确偏导数的发动机模型，模型2指采用中心差分方法的建立的发动机模型。计算耗时指每次改变飞行条件后，模型计算至收敛所需要的时间。气动热力计算次数指模型收敛到指定收敛精度所需要的部件气动热力计算次数。为了便于比较，计算相关数据模型1相对于模型2的百分比偏差，列于表3和表5中。

表2稳态仿真结果

表3稳态仿真结果对比

表4动态仿真结果

表5动态仿真结果对比

通过表2-表5可知，通过本发明建立的部件级偏导数模型和部件级模型的联合计算，在不同条件下收敛所需的时间均有显著下降，稳态计算耗时下降至少50％，动态计算耗时下降至少49％。相应地，稳态和动态所需要的气动热力计算次数下降约90％。由此可知。本发明显著地降低了模型计算的计算量，极大地提高了模型的实时性。

此外，可以看出，应用本发明建立的航空发动机部件级模型单步气动热力计算时间小于1ms，满足航空发动机模型仿真的实时性要求。