黄心刚，钟易成

(南京航空航天大学 能源与动力学院，江苏 南京 210016)



双转子燃气轮机稳态数值计算方法研究

黄心刚，钟易成

(南京航空航天大学 能源与动力学院，江苏 南京 210016)

摘要：以双转子燃机为研究对象，用部件法对燃机进行气动热力性能研究；介绍了燃机特性计算中运用牛顿法寻找燃机的共同工作点，根据Newton-Raphson方法的迭代思想，研究基于高斯消去法的燃机稳态共同工作方程组数值求解方法。实际计算发现：试给若干参数对发动机各部件逐个进行特性计算，并迭代寻找共同工作点，会出现多种问题(迭代不收敛、高斯方程组系数矩阵奇异等)导致计算失败。针对上述问题，采用方程组降阶法、试给参数选择范围弹性限制法改进程序中的算法。结果表明，计算稳定性提高，迭代易收敛。

关键词：双转子燃气轮机；部件法；数值模拟

0引言

燃气涡轮发动机的特性可以用实验方法和计算方法获得。但实验的方法需要研制复杂的设备、投入巨额的资金和消耗巨大的能源，因此实验的方法不可能经常采用。随着计算机运算能力的不断提高、发动机数学模型研究的不断深入，计算机仿真精度也在不断提高，一定程度上弥补了实验方法的不足，尤其是在发动机型号研制过程中，燃气涡轮发动机计算机仿真技术发挥了不可替代的作用。以部件法为燃机特性计算的基本方法，对某型双转子燃气涡轮发动机进行了特性仿真研究。该型燃机是由进气道、低压压气机、高压压气机、燃烧室、高压涡轮、低压涡轮、自由涡轮和喷管组成的，用计算机对这些部件的性能分别进行准确的气动热力性能模拟，则可以准确地模拟整台燃机的气动热力性能，这种建立在准确模拟发动机各部件性能的基础上的发动机性能计算方法，即为部件法。该方法是建立在发动机各部件特性已知的基础上的，因此是计算精度较高的一种方法[1]。

在VS2005平台上用C++语言开发出用于发动机部件以及整机性能通用计算程序和MFC界面，以研究发动机的稳态特性，获取不同工作状态下的各部件共同工作点，为之后的CFD仿真提供参考。在运用Newton- Raphson法求解稳态特性方程组时，经常出现高斯方程组系数矩阵奇异、迭代不收敛等问题，文献[2]中针对上述问题提出降阶法、试给参数选择范围弹性限制法是解决这类问题的有效方法。

本文结合现有燃气轮机的实际情况，并运用以上方法对程序进行了改进，获得了适用于该型燃机稳态特性数值计算的工具。

1发动机数学模型和部件特性

1.1　各部件数学模型的建立

根据本算例中的双转子涡轴燃气轮机，流路简图如图1所示，运用变比热算法计算整机理论特性，并建立各部件热力循环数学模型，包括进气道模型、压气机模型、主燃烧室模型、燃气涡轮模型、动力涡轮模型和尾喷管模型。

图1　双转子涡轴发动机流路简图

以下列出上述各部件的简化数学模型：

1) 进气道模型

类名:PsInlet

进气道进口总参数：

(1)

(2)

进气道出口总参数：

T2=T1

(3)

P2=P1σin

(4)

其中：σin为进气道总压恢复系数。

2) 压气机模型

类名:PsCompressor

低压压气机模型和高压压气机模型的气动热力计算过程基本相同。已知压气机特性方程：

WC-cor=fC1(nC-cor.πC)

(5)

C=fC2(ncor.πC)

(6)

压气机出口总参数：

PC-out=PC-inπC

(7)

TC-out=ft(HC-out)

(8)

其中：nC-cor为换算转速，WC-cor为换算流量，πC为压比，HC-out为压气机出口气体的总焓。

3) 燃烧室模型

类名:PsCombustor

燃烧室出口总参数：

P8=P7σB

(9)

T8=ft-B(H8,fa8)

(10)

燃烧室出口流量：

W8=W7+Wf

(11)

其中：σB为燃烧室总压恢复系数，H8为燃烧室出口气体总焓，fa8为油气比，Wf为供油量。

4) 燃气涡轮模型及动力涡轮模型

类名:PsTurb&PsFreeTurb

低压涡轮模型、高压涡轮模型以及动力涡轮模型的气动热力计算过程基本相同。已知涡轮特性方程：

WT-cor=fT1(nT-cor,πT)

(12)

ηT=fT2(nT-cor,πT)

(13)

涡轮出口总参数：

PT-out=PT-in/πT

(14)

TT-out=ft-B(HT-out,fa-Tout)

(15)

涡轮出口流量：

Wout=Win+Wcool

(16)

其中：πT为涡轮落压比，HT-out为涡轮出口气体的总焓，fa-Tout为涡轮出口气体的油气比，Wcool为涡轮冷却气流流量。

5) 尾喷管模型

类名:PsNozzle

燃机的排气喷管采用的是亚声速扩压型排气管，涡轮后的燃气通过排气管膨胀至环境大气压P0。

喷管出口总参数：

P16=P15σNZ

(17)

T16=T15

(18)

喷管排气面积：

(19)

排气喷管剩余推力：

(20)

其中：σNZ为总压恢复系数，ρ16、V16分别为喷管出口气流密度及速度，φNZ为喷管速度分布系数。

针对本文主要研究的涡轴发动机而言，是由与燃气涡轮没有机械联系的动力涡轮输出功率，故在计算发动机总体性能时需要计算动力涡轮功率，即:

(21)

则整台发动机的轴输出功率为:

Psh=PDTg

(22)

耗油率

(23)

1.2　发动机各部件特性的录入

部件法之所以是计算精度比较高的一种模拟发动机性能的方法，是因为该方法是建立在发动机各部件特性已知的基础上的。发动机各部件特性是由若干条连续的曲线组成的，计算机中只能 存储离散的点信息，部件特性计算要解决的首要问题是将部件特性曲线进行离散化，同时充分保留原始曲线的变化范围、变化趋势等重要信息[2]。

这里以本算例中的低压压气机为例，介绍其特性录入的方法。从燃机几何模型中，取出低压压气机部分先用TurboGrid软件划分用网格后，再用CFX软件进行CFD计算，得到低压压气机特性曲线图，如图2、图3所示。

图2　低压压气机流量-压比曲线

图3　低压压气机流量-效率曲线

图2及图3中不同的曲线代表不同物理转速下低压压气机的流量-压比和流量-效率曲线，物理转速范围从22000~28500转(共计8个)，每个物理转速下又分别计算了8个不同的工作点。在低压压气机对应的PsCom- pressor类创建对象LCompressor，曲线上每一个点对应的各个参数存放在一个二维数组中，为了适用于各种不同的大气条件，将物理转速、流量用换算转速、换算流量来代替，所以此二维数组中的参数需要按照如下算法进行换算：

(24)

(25)

由图2、图3可知，低转速下，流量越小时，曲线越垂直于y轴；高转速下，流量越大时，曲线越垂直于x轴，都将导致插值运算的不稳定，从而引起较大的误差，参考文献[4]中的方法，将每条曲线上离散各点用切比雪夫多项式拟合，用一组变换平缓的直线截取曲线上各点，将二维数组内的数据重新整合，并形成一组新的曲线——β辅助线，如图4所示。

图4　压气机特性辅助线示意图

算例中燃机涡轮部件的参数导入方法与所介绍的压气机部件导入方法基本类似，导入完毕后，根据给定的转速，用插值模块(类名:PsInterpolate)中的二元不规则插值法计算出该转速下，β值从0变至100(此两点对应CFD计算给出的同一转速下的参数始末端点)所对应的各个参数值(Wcor、π、)。

2算法研究和程序设计

2.1　稳态算法研究

燃机各部件的特性，大多是由实验或是CFD计算得到的非线性曲线，要想写出需求参数与已知参数之间的具体函数关系式十分困难。因此只能根据已知参数、部件特性和各部件共同工作条件，对燃机各部件按气流方向逐个进行计算。在各部件的特性图上，每一个坐标对应一个工作点，由于发动机各部件气动热力方程组制约关系，无法首先确定工作点在图上对应的坐标，导致有些参数不能求出，所以需要先试给一些未知参数，进行试凑迭代计算，迭代过程工作量较大，但结果比较精确[2]。

对于本算例中的双转子燃气轮机，给定环境条件(高度H、马赫数Ma)、高压转子转速(nh)、动力涡轮输出功率(PDT)、喷管出口面积(Ae)，如此确定稳定的工作状态。但还缺少6个参数(设为y1~y6)使得发动机气动热力方程组有唯一解，所以现在试给这六个参数：

1) 进气道进口流量(Win)；

2) 低压转子转速(n1)；

3) 高压涡轮前燃气温度(T8)；

4) 高压涡轮落压比(πTh)；

5) 低压涡轮落压比(πTl)；

6) 动力涡轮落压比(πTd)。

给出上述参数后，通过检查试给值，若试给值满足则发动机各部件气动热力方程组有唯一解，若不满足则需要重新给定试给参数进行迭代求解，可由以下6个共同工作条件方程来检查试给值：

(26)

其中：Pt,h、Pt,l、Pc,h、Pc,l分别为高、低压涡轮和高、低压压气机的功率；m,h、m,l分别为高、低压涡轮的机械效率；Wt,h,in为计算得出的高压涡轮进口燃气流量；Wt,h,map、Wt,l,map、Wt,d,map分别为从涡轮特性图上插值得出的高、低压涡轮以及动力涡轮的燃气流量；为计算得出的尾喷管气流出口面积；Z1~Z6为各平衡方程式残量，若试给值满足则6个残量都应该为0(实际计算过程中只要6个值小于某一允许误差)。若6个残量有一个不满足要求，则要重新给定试给参数再进行计算，从而求出新的一组残量。显然，残量都为试给值的函数，写成：

Z=f(y)

(27)

z和y均为6维向量，即:

Z=[Z1,Z2,Z3,Z4,Z5,Z6]T

(28)

y=[y1,y2,y3,y4,y5,y6]

(29)

f为向量函数:

f=[f1,f2,f3,f4,f5,f6]T

(30)

要使6个试给值满足共同工作条件，就要求解下列使Z1~Z6符合条件的方程组：

z=f(y)=0

(31)

求解上述非线性方程组，采用线性化方法能较容易实现，即将非线性方程近似成一组线性方程，构造一种迭代格式，来逼近所求的解。本算例中采用Newton-Raph son方法进行求解，根据此非线性方程组，用该方法的迭代公式：

y(i+1)=yi-J-1f(yi)i=0,1,…

(32)

其中，yi为迭代至当前步数的试给值向量，y(i+1)为下一次迭代完毕后得到的修正值，J为雅克比矩阵[2]。

2.2　程序设计

针对该燃机算例的程序设计，主要围绕各部件模块本身以及部件连接截面的气动热力计算、以转轴相连的部件之间的机械参数计算、特性图插值计算和最后的整机工作点迭代计算。程序基本设计思路流程如图5[5]。

图5　程序计算流程

为了形成一款通用涡轴发动机稳态计算工具，进行了相关程序的MFC界面开发，部分界面如图6所示。

图6　计算软件界面

3计算结果分析

针对算例燃机，通过所开发软件与CFD计算结果的对比，验证软件计算的准确度，同时针对软件设计的不足进行了改进。

现给出给定的设计点参数以及为了使算例燃机各部件气动热力计算方程组成立的6个试给参数。

具体参数见表1。

表1　算例燃机稳态计算给定值和试给值

将上述参数输入软件主界面中，并依次导入各部件特性曲线图以及本燃机所采用的燃气、燃油的热值等参数。软件计算和CFD计算得到的部分参数对比如表2所示。

通过算例验证了Newton-Raphson法的可行性和准确性，该算法能够准确获得发动机给定工作状态下的稳态性能参数。但是通过多次使用发现，在软件计算试给值偏离真实值较大的情况下，常常出现迭代过程长时间不收敛，甚至出现高斯方程组系数矩阵奇异导致无法求解的情况，查看有关资料可知，其原因可能是该高斯方程组的阶数过高；试给参数时未考虑实际工况。根据文献提供的方法：1) 降阶法；2) 试给参数取值范围的弹性限制法[2]，并结合算例燃机本身的实际情况，现对软件的设计及迭代算法进行改进。

表2　软件与CFD设计点计算结果对比(压气机、涡轮部分出口截面参数)

图7　采用降阶法程序计算流程

在检查调整试给值时，需要规定试给值的取值范围，即使用弹性限制法对参数进行约束。在每一次迭代过程中，对试给参数的调整并不是单调递增或是单调递减，参加下一轮迭代的参数往往是按照迭代格式所给出的因子呈现出波动的情形。所以极有可能出现2轮迭代过程中出现相同的取值(或者参数差值的绝对值很小)，在建立线性方程时出现等价方程，从而导致高斯方程组系数矩阵奇异，弹性限制法能使参数的两次取值不同来规避这一情况。

本文给出稳态计算程序改进前和改进后出口截面面积参数的迭代残差图，如图8所示。

图8　喷管出口面积计算残差图

图上曲线分别代表喷管出口面积在算法改进前和改进后计算值与真实值的百分比误差，从图8可以看出采用降阶法和参数取值范围的弹性限制法后，迭代次数较之前减少4步，而在实际的软件使用中由于对试给值有效的限制，出现使计算中止的坏值的频率明显下降。

在软件计算稳定性有所提高的基础上，针对给定喷管出口面积的情况，计算出压气机特性图上在Ae=0.098m2时的共同工作线如图9所示。

图9　Ae=0.098m2时，低压压气机特性图上共同工作线

改进后的软件同样计算了算例燃机在设计点的情况，给出了改进前所需求出的6个试给值，其值如表3。

表3　改进前的试给值

4结论

研究了双转子燃气轮机稳态数值计算方法，通过算法研究以及相应的软件设计，可以得出以下结论：

1) 在燃气轮机特性计算中，国内外普遍采用的对多元非线性方程组进行线性化迭代的求解方法，通过算例验证了其能较准确的获得各部件在大部分工作状态下的参数，但在计算过程中，易出现迭代计算不稳定，长时间不收敛等问题。

2) 针对所选用的算例燃机，参考相应文献对计算程

序进行改进，采用了降阶法以及试给参数选择范围弹性限制法，发现迭代收敛速度较改进前计算成功的情况有所加快，迭代过程中因出现高斯方程组系数矩阵奇异而无法进行计算的情况明显减少。

3) 以Newton-Raphson方法为主要核心算法，进行相关的软件开发，并以实际存在的某台双转子燃气轮机为验证算例，软件能够准确地获得不同工作状态下发动机各部件的气动热力参数，并为整台发动机的性能分析、设计优化提供帮助。

参考文献:

[1] 骆广琦，桑增产.航空燃气涡轮发动机数值仿真[M]. 北京：国防工业出版社，2007.4.

[2] 朱行健，王雪瑜.燃气轮机工作原理及性能 [M]. 北京：科学出版社，1992.12.

[3] 航空发动机设计手册.第6册[M]. 北京:航空工业出版社，2001.9.

[4] Claus Riegler, Michael Bauer, Joachim Kurzke. Some Aspects of Modeling Compressor Behavior in Gas Turbine Performance Calculations[R]. ASME, Vol.123, 1970, 4.

[5] Brian P., Curlett and James. Object-Oriented Approach for Gas Turbine Engine Simulation[R]. NASA TM-106970.

Research on Steady-State Numerical Calculation Method of A Dual-Rotor Gas Turbine

HUANG Xin-gang, ZHONG Yi-cheng

(College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)

Abstract:The dual-rotor gas turbine is taken for the objecr of study and the gas turbine component method is used to research on the aero and thermodynamic properties. This paper introduces the common combustion engine operating point which is sought using Newton’s method in the turbine characteristic calculation. According to the method of iterative Newton-Raphson ideas, It also researches on the equations numerical method in the steady-state combustion engine common operation based on the Gaussion elimination. It is found in the actual calculation that the parameters are used to try calculating the characteristic of the engine components one by one and the common operating point is found out through the iteration. A variety of problems (iteration does not converge, Gauss equation coefficient matrix singularity, etc) come to existence, which lead the calculation to fail. For these problems, it tries using equations reduction method and the parameter selection Linit method to improve flexibility in the program algorithm. The results show that the calculated stability is improved and it is easy to converge the iteration.

Keywords:dual-rotor gas turbine; component method; numerical simulation

中图分类号：V231.3

文献标志码：A

文章编号：1671-5276(2015)02-0050-05

作者简介：黄心刚(1989-)，男，江苏无锡人，硕士研究生，从事发动机总体性能研究。

收稿日期：2014-10-30