雙轉(zhuǎn)子燃氣輪機穩(wěn)態(tài)數(shù)值計算方法研究

黃心剛，鐘易成

(南京航空航天大學 能源與動力學院，江蘇 南京 210016)

?

雙轉(zhuǎn)子燃氣輪機穩(wěn)態(tài)數(shù)值計算方法研究

黃心剛，鐘易成

(南京航空航天大學 能源與動力學院，江蘇 南京 210016)

摘要：以雙轉(zhuǎn)子燃機為研究對象，用部件法對燃機進行氣動熱力性能研究；介紹了燃機特性計算中運用牛頓法尋找燃機的共同工作點，根據(jù)Newton-Raphson方法的迭代思想，研究基于高斯消去法的燃機穩(wěn)態(tài)共同工作方程組數(shù)值求解方法。實際計算發(fā)現(xiàn)：試給若干參數(shù)對發(fā)動機各部件逐個進行特性計算，并迭代尋找共同工作點，會出現(xiàn)多種問題(迭代不收斂、高斯方程組系數(shù)矩陣奇異等)導致計算失敗。針對上述問題，采用方程組降階法、試給參數(shù)選擇范圍彈性限制法改進程序中的算法。結(jié)果表明，計算穩(wěn)定性提高，迭代易收斂。

關(guān)鍵詞：雙轉(zhuǎn)子燃氣輪機；部件法；數(shù)值模擬

0引言

燃氣渦輪發(fā)動機的特性可以用實驗方法和計算方法獲得。但實驗的方法需要研制復雜的設(shè)備、投入巨額的資金和消耗巨大的能源，因此實驗的方法不可能經(jīng)常采用。隨著計算機運算能力的不斷提高、發(fā)動機數(shù)學模型研究的不斷深入，計算機仿真精度也在不斷提高，一定程度上彌補了實驗方法的不足，尤其是在發(fā)動機型號研制過程中，燃氣渦輪發(fā)動機計算機仿真技術(shù)發(fā)揮了不可替代的作用。以部件法為燃機特性計算的基本方法，對某型雙轉(zhuǎn)子燃氣渦輪發(fā)動機進行了特性仿真研究。該型燃機是由進氣道、低壓壓氣機、高壓壓氣機、燃燒室、高壓渦輪、低壓渦輪、自由渦輪和噴管組成的，用計算機對這些部件的性能分別進行準確的氣動熱力性能模擬，則可以準確地模擬整臺燃機的氣動熱力性能，這種建立在準確模擬發(fā)動機各部件性能的基礎(chǔ)上的發(fā)動機性能計算方法，即為部件法。該方法是建立在發(fā)動機各部件特性已知的基礎(chǔ)上的，因此是計算精度較高的一種方法[1]。

在VS2005平臺上用C++語言開發(fā)出用于發(fā)動機部件以及整機性能通用計算程序和MFC界面，以研究發(fā)動機的穩(wěn)態(tài)特性，獲取不同工作狀態(tài)下的各部件共同工作點，為之后的CFD仿真提供參考。在運用Newton- Raphson法求解穩(wěn)態(tài)特性方程組時，經(jīng)常出現(xiàn)高斯方程組系數(shù)矩陣奇異、迭代不收斂等問題，文獻[2]中針對上述問題提出降階法、試給參數(shù)選擇范圍彈性限制法是解決這類問題的有效方法。

本文結(jié)合現(xiàn)有燃氣輪機的實際情況，并運用以上方法對程序進行了改進，獲得了適用于該型燃機穩(wěn)態(tài)特性數(shù)值計算的工具。

1發(fā)動機數(shù)學模型和部件特性

1.1　各部件數(shù)學模型的建立

根據(jù)本算例中的雙轉(zhuǎn)子渦軸燃氣輪機，流路簡圖如圖1所示，運用變比熱算法計算整機理論特性，并建立各部件熱力循環(huán)數(shù)學模型，包括進氣道模型、壓氣機模型、主燃燒室模型、燃氣渦輪模型、動力渦輪模型和尾噴管模型。

圖1　雙轉(zhuǎn)子渦軸發(fā)動機流路簡圖

以下列出上述各部件的簡化數(shù)學模型：

1) 進氣道模型

類名:PsInlet

進氣道進口總參數(shù)：

(1)

(2)

進氣道出口總參數(shù)：

T2=T1

(3)

P2=P1σin

(4)

其中：σin為進氣道總壓恢復系數(shù)。

2) 壓氣機模型

類名:PsCompressor

低壓壓氣機模型和高壓壓氣機模型的氣動熱力計算過程基本相同。已知壓氣機特性方程：

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

(5)

C=fC2(ncor.πC)

(6)

壓氣機出口總參數(shù)：

PC-out=PC-inπC

(7)

TC-out=ft(HC-out)

(8)

其中：nC-cor為換算轉(zhuǎn)速，WC-cor為換算流量，πC為壓比，HC-out為壓氣機出口氣體的總焓。

3) 燃燒室模型

類名:PsCombustor

燃燒室出口總參數(shù)：

P8=P7σB

(9)

T8=ft-B(H8,fa8)

(10)

燃燒室出口流量：

W8=W7+Wf

(11)

其中：σB為燃燒室總壓恢復系數(shù)，H8為燃燒室出口氣體總焓，fa8為油氣比，Wf為供油量。

4) 燃氣渦輪模型及動力渦輪模型

類名:PsTurb&PsFreeTurb

低壓渦輪模型、高壓渦輪模型以及動力渦輪模型的氣動熱力計算過程基本相同。已知渦輪特性方程：

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

(12)

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

(13)

渦輪出口總參數(shù)：

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

燃機的排氣噴管采用的是亞聲速擴壓型排氣管，渦輪后的燃氣通過排氣管膨脹至環(huán)境大氣壓P0。

噴管出口總參數(shù)：

P16=P15σNZ

(17)

T16=T15

(18)

噴管排氣面積：

(19)

排氣噴管剩余推力：

(20)

其中：σNZ為總壓恢復系數(shù)，ρ16、V16分別為噴管出口氣流密度及速度，φNZ為噴管速度分布系數(shù)。

針對本文主要研究的渦軸發(fā)動機而言，是由與燃氣渦輪沒有機械聯(lián)系的動力渦輪輸出功率，故在計算發(fā)動機總體性能時需要計算動力渦輪功率，即:

(21)

則整臺發(fā)動機的軸輸出功率為:

Psh=PDTg

(22)

耗油率

(23)

1.2　發(fā)動機各部件特性的錄入

部件法之所以是計算精度比較高的一種模擬發(fā)動機性能的方法，是因為該方法是建立在發(fā)動機各部件特性已知的基礎(chǔ)上的。發(fā)動機各部件特性是由若干條連續(xù)的曲線組成的，計算機中只能 存儲離散的點信息，部件特性計算要解決的首要問題是將部件特性曲線進行離散化，同時充分保留原始曲線的變化范圍、變化趨勢等重要信息[2]。

這里以本算例中的低壓壓氣機為例，介紹其特性錄入的方法。從燃機幾何模型中，取出低壓壓氣機部分先用TurboGrid軟件劃分用網(wǎng)格后，再用CFX軟件進行CFD計算，得到低壓壓氣機特性曲線圖，如圖2、圖3所示。

圖2　低壓壓氣機流量-壓比曲線

圖3　低壓壓氣機流量-效率曲線

圖2及圖3中不同的曲線代表不同物理轉(zhuǎn)速下低壓壓氣機的流量-壓比和流量-效率曲線，物理轉(zhuǎn)速范圍從22000~28500轉(zhuǎn)(共計8個)，每個物理轉(zhuǎn)速下又分別計算了8個不同的工作點。在低壓壓氣機對應的PsCom- pressor類創(chuàng)建對象LCompressor，曲線上每一個點對應的各個參數(shù)存放在一個二維數(shù)組中，為了適用于各種不同的大氣條件，將物理轉(zhuǎn)速、流量用換算轉(zhuǎn)速、換算流量來代替，所以此二維數(shù)組中的參數(shù)需要按照如下算法進行換算：

(24)

(25)

由圖2、圖3可知，低轉(zhuǎn)速下，流量越小時，曲線越垂直于y軸；高轉(zhuǎn)速下，流量越大時，曲線越垂直于x軸，都將導致插值運算的不穩(wěn)定，從而引起較大的誤差，參考文獻[4]中的方法，將每條曲線上離散各點用切比雪夫多項式擬合，用一組變換平緩的直線截取曲線上各點，將二維數(shù)組內(nèi)的數(shù)據(jù)重新整合，并形成一組新的曲線——β輔助線，如圖4所示。

圖4　壓氣機特性輔助線示意圖

算例中燃機渦輪部件的參數(shù)導入方法與所介紹的壓氣機部件導入方法基本類似，導入完畢后，根據(jù)給定的轉(zhuǎn)速，用插值模塊(類名:PsInterpolate)中的二元不規(guī)則插值法計算出該轉(zhuǎn)速下，β值從0變至100(此兩點對應CFD計算給出的同一轉(zhuǎn)速下的參數(shù)始末端點)所對應的各個參數(shù)值(Wcor、π、)。

2算法研究和程序設(shè)計

2.1　穩(wěn)態(tài)算法研究

燃機各部件的特性，大多是由實驗或是CFD計算得到的非線性曲線，要想寫出需求參數(shù)與已知參數(shù)之間的具體函數(shù)關(guān)系式十分困難。因此只能根據(jù)已知參數(shù)、部件特性和各部件共同工作條件，對燃機各部件按氣流方向逐個進行計算。在各部件的特性圖上，每一個坐標對應一個工作點，由于發(fā)動機各部件氣動熱力方程組制約關(guān)系，無法首先確定工作點在圖上對應的坐標，導致有些參數(shù)不能求出，所以需要先試給一些未知參數(shù)，進行試湊迭代計算，迭代過程工作量較大，但結(jié)果比較精確[2]。

對于本算例中的雙轉(zhuǎn)子燃氣輪機，給定環(huán)境條件(高度H、馬赫數(shù)Ma)、高壓轉(zhuǎn)子轉(zhuǎn)速(nh)、動力渦輪輸出功率(PDT)、噴管出口面積(Ae)，如此確定穩(wěn)定的工作狀態(tài)。但還缺少6個參數(shù)(設(shè)為y1~y6)使得發(fā)動機氣動熱力方程組有唯一解，所以現(xiàn)在試給這六個參數(shù)：

1) 進氣道進口流量(Win)；

2) 低壓轉(zhuǎn)子轉(zhuǎn)速(n1)；

3) 高壓渦輪前燃氣溫度(T8)；

4) 高壓渦輪落壓比(πTh)；

5) 低壓渦輪落壓比(πTl)；

6) 動力渦輪落壓比(πTd)。

給出上述參數(shù)后，通過檢查試給值，若試給值滿足則發(fā)動機各部件氣動熱力方程組有唯一解，若不滿足則需要重新給定試給參數(shù)進行迭代求解，可由以下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個殘量有一個不滿足要求，則要重新給定試給參數(shù)再進行計算，從而求出新的一組殘量。顯然，殘量都為試給值的函數(shù)，寫成：

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為向量函數(shù):

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

(30)

要使6個試給值滿足共同工作條件，就要求解下列使Z1~Z6符合條件的方程組：

z=f(y)=0

(31)

求解上述非線性方程組，采用線性化方法能較容易實現(xiàn)，即將非線性方程近似成一組線性方程，構(gòu)造一種迭代格式，來逼近所求的解。本算例中采用Newton-Raph son方法進行求解，根據(jù)此非線性方程組，用該方法的迭代公式：

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

(32)

其中，yi為迭代至當前步數(shù)的試給值向量，y(i+1)為下一次迭代完畢后得到的修正值，J為雅克比矩陣[2]。

2.2　程序設(shè)計

針對該燃機算例的程序設(shè)計，主要圍繞各部件模塊本身以及部件連接截面的氣動熱力計算、以轉(zhuǎn)軸相連的部件之間的機械參數(shù)計算、特性圖插值計算和最后的整機工作點迭代計算。程序基本設(shè)計思路流程如圖5[5]。

圖5　程序計算流程

為了形成一款通用渦軸發(fā)動機穩(wěn)態(tài)計算工具，進行了相關(guān)程序的MFC界面開發(fā)，部分界面如圖6所示。

圖6　計算軟件界面

3計算結(jié)果分析

針對算例燃機，通過所開發(fā)軟件與CFD計算結(jié)果的對比，驗證軟件計算的準確度，同時針對軟件設(shè)計的不足進行了改進。

現(xiàn)給出給定的設(shè)計點參數(shù)以及為了使算例燃機各部件氣動熱力計算方程組成立的6個試給參數(shù)。

具體參數(shù)見表1。

表1　算例燃機穩(wěn)態(tài)計算給定值和試給值

將上述參數(shù)輸入軟件主界面中，并依次導入各部件特性曲線圖以及本燃機所采用的燃氣、燃油的熱值等參數(shù)。軟件計算和CFD計算得到的部分參數(shù)對比如表2所示。

通過算例驗證了Newton-Raphson法的可行性和準確性，該算法能夠準確獲得發(fā)動機給定工作狀態(tài)下的穩(wěn)態(tài)性能參數(shù)。但是通過多次使用發(fā)現(xiàn)，在軟件計算試給值偏離真實值較大的情況下，常常出現(xiàn)迭代過程長時間不收斂，甚至出現(xiàn)高斯方程組系數(shù)矩陣奇異導致無法求解的情況，查看有關(guān)資料可知，其原因可能是該高斯方程組的階數(shù)過高；試給參數(shù)時未考慮實際工況。根據(jù)文獻提供的方法：1) 降階法；2) 試給參數(shù)取值范圍的彈性限制法[2]，并結(jié)合算例燃機本身的實際情況，現(xiàn)對軟件的設(shè)計及迭代算法進行改進。

表2　軟件與CFD設(shè)計點計算結(jié)果對比(壓氣機、渦輪部分出口截面參數(shù))

圖7　采用降階法程序計算流程

在檢查調(diào)整試給值時，需要規(guī)定試給值的取值范圍，即使用彈性限制法對參數(shù)進行約束。在每一次迭代過程中，對試給參數(shù)的調(diào)整并不是單調(diào)遞增或是單調(diào)遞減，參加下一輪迭代的參數(shù)往往是按照迭代格式所給出的因子呈現(xiàn)出波動的情形。所以極有可能出現(xiàn)2輪迭代過程中出現(xiàn)相同的取值(或者參數(shù)差值的絕對值很小)，在建立線性方程時出現(xiàn)等價方程，從而導致高斯方程組系數(shù)矩陣奇異，彈性限制法能使參數(shù)的兩次取值不同來規(guī)避這一情況。

本文給出穩(wěn)態(tài)計算程序改進前和改進后出口截面面積參數(shù)的迭代殘差圖，如圖8所示。

圖8　噴管出口面積計算殘差圖

圖上曲線分別代表噴管出口面積在算法改進前和改進后計算值與真實值的百分比誤差，從圖8可以看出采用降階法和參數(shù)取值范圍的彈性限制法后，迭代次數(shù)較之前減少4步，而在實際的軟件使用中由于對試給值有效的限制，出現(xiàn)使計算中止的壞值的頻率明顯下降。

在軟件計算穩(wěn)定性有所提高的基礎(chǔ)上，針對給定噴管出口面積的情況，計算出壓氣機特性圖上在Ae=0.098m2時的共同工作線如圖9所示。

圖9　Ae=0.098m2時，低壓壓氣機特性圖上共同工作線

改進后的軟件同樣計算了算例燃機在設(shè)計點的情況，給出了改進前所需求出的6個試給值，其值如表3。

表3　改進前的試給值

4結(jié)論

研究了雙轉(zhuǎn)子燃氣輪機穩(wěn)態(tài)數(shù)值計算方法，通過算法研究以及相應的軟件設(shè)計，可以得出以下結(jié)論：

1) 在燃氣輪機特性計算中，國內(nèi)外普遍采用的對多元非線性方程組進行線性化迭代的求解方法，通過算例驗證了其能較準確的獲得各部件在大部分工作狀態(tài)下的參數(shù)，但在計算過程中，易出現(xiàn)迭代計算不穩(wěn)定，長時間不收斂等問題。

2) 針對所選用的算例燃機，參考相應文獻對計算程

序進行改進，采用了降階法以及試給參數(shù)選擇范圍彈性限制法，發(fā)現(xiàn)迭代收斂速度較改進前計算成功的情況有所加快，迭代過程中因出現(xiàn)高斯方程組系數(shù)矩陣奇異而無法進行計算的情況明顯減少。

3) 以Newton-Raphson方法為主要核心算法，進行相關(guān)的軟件開發(fā)，并以實際存在的某臺雙轉(zhuǎn)子燃氣輪機為驗證算例，軟件能夠準確地獲得不同工作狀態(tài)下發(fā)動機各部件的氣動熱力參數(shù)，并為整臺發(fā)動機的性能分析、設(shè)計優(yōu)化提供幫助。

參考文獻:

[1] 駱廣琦，桑增產(chǎn).航空燃氣渦輪發(fā)動機數(shù)值仿真[M]. 北京：國防工業(yè)出版社，2007.4.

[2] 朱行健，王雪瑜.燃氣輪機工作原理及性能 [M]. 北京：科學出版社，1992.12.

[3] 航空發(fā)動機設(shè)計手冊.第6冊[M]. 北京:航空工業(yè)出版社，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-)，男，江蘇無錫人，碩士研究生，從事發(fā)動機總體性能研究。

收稿日期：2014-10-30