雙轉子燃氣輪機穩態數值計算方法研究

黃心剛，鐘易成

(南京航空航天大學 能源與動力學院，江蘇 南京 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