基于0D/2D耦合的渦扇發(fā)動機(jī)總體性能計算研究

李 燕，葛 寧

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

1 引言

為提高航空發(fā)動機(jī)的性能、效率及可靠性，同時縮短研制周期、降低研制成本，全球多個航空發(fā)動機(jī)研制機(jī)構(gòu)在開發(fā)并完善發(fā)動機(jī)總體性能仿真系統(tǒng)。傳統(tǒng)的性能仿真軟件，如GasTurb、TURBOMATCH、PYTHIA、T-MATS等，都是以零維(0D)模型為基礎(chǔ)，通過熱力學(xué)平衡方程計算發(fā)動機(jī)各部件進(jìn)出口的循環(huán)參數(shù)，從而得到推力、耗油率等總體性能參數(shù)。但這些軟件有兩個缺陷：①未考慮發(fā)動機(jī)的幾何構(gòu)型及部件內(nèi)部氣體的流動狀況；②非設(shè)計點性能預(yù)估的精度很大程度上依賴于部件特性圖的準(zhǔn)確性，而部件特性要么通過通用特性圖縮放得到，要么通過大量實驗獲取。前者太過于依賴經(jīng)驗，計算精度低；后者一旦部件有所改動，就需重新實驗來獲得部件特性，耗資巨大。

鑒于此，美國在NPSS(推進(jìn)系統(tǒng)數(shù)值仿真)計劃中提出了數(shù)值縮放技術(shù)[1-2]，將發(fā)動機(jī)內(nèi)某一部件的高精度模型(三維(3D))和其他低精度模型(0D)耦合在一起進(jìn)行數(shù)值模擬。目前這項技術(shù)在國內(nèi)外均有應(yīng)用和發(fā)展：Reed等[3]將風(fēng)扇的3D模型集成在傳統(tǒng)的0D仿真平臺上，實現(xiàn)了0D/3D耦合的發(fā)動機(jī)整機(jī)仿真；陳玉春等[4]將3D尾噴管的CFD計算程序和0D發(fā)動機(jī)性能計算程序結(jié)合在一起，實現(xiàn)了數(shù)值縮放技術(shù)在發(fā)動機(jī)混合器和尾噴管中的應(yīng)用；Mund等[5]借鑒數(shù)值縮放技術(shù)中構(gòu)建混合精度模型的思路，從流道設(shè)計出發(fā)，建立了發(fā)動機(jī)部分部件(進(jìn)氣道、風(fēng)扇、外涵)的二維(2D)幾何模型，發(fā)展了0D/2D耦合的發(fā)動機(jī)總體性能計算方法。

本文基于S1流面對某型雙軸混合排氣渦扇發(fā)動機(jī)低壓壓氣機(jī)和低壓渦輪進(jìn)行了初步設(shè)計，并將基于周向平均N-S方程的通流模型與T-MATS平臺上的0D模型結(jié)合起來，改數(shù)值縮放技術(shù)中的3D高精度模型為2D模型，解決了3D數(shù)值模擬計算難度大且耗時長的問題，實現(xiàn)了0D/2D耦合方法在低壓壓氣機(jī)和低壓渦輪設(shè)計中的應(yīng)用。

2 基于0D/2D耦合的渦扇發(fā)動機(jī)總體性能計算方法

2.1 0D的T-MATS仿真平臺

T-MATS(Toolbox for the Modeling and Analysis of Thermodynamic Systems)是NASA格林研究中心在2014年公開發(fā)表的一套用于熱力系統(tǒng)建模、仿真及控制的工具箱[6]。作為MATLAB/Simulink平臺上的一個子庫，其模塊化界面清晰明了，使用靈活；內(nèi)部程序的開源性不僅方便用戶創(chuàng)建或優(yōu)化模塊，而且還便于其與用戶開發(fā)的程序集成。

2.2 基于周向平均N-S方程的2D通流計算

采用課題組開發(fā)的通流計算程序[7]，作為發(fā)動機(jī)部件2D數(shù)值仿真和分析的工具。程序以周向平均N-S方程作為控制方程；附加項中，采用Simon模型[8]模擬無粘葉片力，采用分布式損失模型[9-11]模擬粘性葉片力，非設(shè)計點的落后角和損失模型采用Miller[12]提出的方法?？刂品匠糖蠼馐褂糜邢摅w積法，時間離散采用多步Runge-Kutta顯示格式[13]。該程序計算速度快，能反映真實子午流場。

2.3 基于0D/2D耦合的渦扇發(fā)動機(jī)設(shè)計點總體性能計算方法

圖1為0D/2D耦合計算模型示意圖。圖中，灰色區(qū)域是基于T-MATS平臺搭建的0D模型；白色區(qū)域是2D通流模型。

設(shè)計點求解時，首先由通流計算得出低壓壓氣機(jī)S2流面的流場，出口總溫、總壓經(jīng)過質(zhì)量加權(quán)平均后代入T-MATS中作為核心機(jī)進(jìn)口邊界條件，再將高壓渦輪出口參數(shù)作為低壓渦輪進(jìn)口條件代入通流計算程序。部件匹配必須保證流量平衡和功率平衡兩個守恒條件。其中T-MATS本身通過流量平衡方程和功率平衡方程求解各個部件出口參數(shù)，而低壓壓氣機(jī)和低壓渦輪采用通流程序分別計算，需要驗證是否滿足平衡條件。具體計算流程(圖2)為：

(1)按設(shè)計點要求給定低壓壓氣機(jī)出口背壓p2和低壓轉(zhuǎn)速NL，由通流計算得到低壓壓氣機(jī)流量w1、功耗Wcl及出口總溫、總壓。

(2)T-MATS平臺上，初設(shè)定一個內(nèi)涵流量和燃油量wf(在0D計算結(jié)果附近取值)，計算得到高壓渦輪出口總溫、總壓。

(3)以高壓渦輪出口狀態(tài)作為低壓渦輪進(jìn)口條件，低壓轉(zhuǎn)速仍設(shè)為NL，按設(shè)計點要求設(shè)置低壓渦輪出口背壓，通過通流計算得到低壓渦輪流量做功量Wtl及出口總溫、總壓。

(4)檢驗是否滿足流量平衡和功率平衡條件，如果不滿足則分別重新修改和wf，重復(fù)步驟(2)～(4)，直到同時滿足兩個平衡條件。

(5) 通過T-MATS計算外涵、混合器、尾噴管，得到總體性能參數(shù)。

2.4 基于0D/2D耦合的渦扇發(fā)動機(jī)非設(shè)計點總體性能計算方法

與設(shè)計點計算不同，非設(shè)計點計算無設(shè)計要求限制，所以除了保證流量平衡和功率平衡外，還需滿足混合器進(jìn)口靜壓相等，即保證外涵總壓與內(nèi)涵總壓之比在0.95～1.05[14]，以減小混合器損失。因此，非設(shè)計點計算流程中，選定某一低壓轉(zhuǎn)速后，首先初設(shè)定一個 p2，執(zhí)行設(shè)計點步驟中的(2)～(4)；然后調(diào)整 p2，直到滿足壓力平衡要求(如圖2虛線內(nèi)指示)，最后執(zhí)行步驟(5)。

3 低壓壓氣機(jī)和低壓渦輪的初步設(shè)計

利用吳仲華[15]的兩類流面理論，假設(shè)壓氣機(jī)和渦輪在葉根、葉中、葉尖處存在S1流面，在這3個流面上構(gòu)造2D葉型，通過重心徑向線積疊得到完整葉片，再將得到的初步葉型代入通流計算程序。

3.1 基于S1流面的低壓壓氣機(jī)初步設(shè)計

低壓壓氣機(jī)的初步設(shè)計及性能驗證參考了Denton等[16]的設(shè)計思路，具體方法是：①根據(jù)設(shè)計要求(表1)和0D算出的低壓壓氣機(jī)出口參數(shù)，確定壓氣機(jī)的流道和中徑處的氣流角；②給定等環(huán)量分布規(guī)律，計算出從葉根到葉尖不同半徑處的氣流角；③根據(jù)氣流角，采用雙圓弧葉型構(gòu)造方法設(shè)計二維葉型，再通過徑向積疊構(gòu)建三維葉型。

表1 低壓壓氣機(jī)初步設(shè)計要求Table 1 Preliminary design requirements of the low pressure compressor

初步設(shè)計時不考慮葉片的三維特性，不涉及葉片的彎掠，只在葉根、葉中、葉尖3個流面設(shè)計二維葉型。通過調(diào)整影響壓氣機(jī)性能的輪轂比、壓氣機(jī)進(jìn)出口馬赫數(shù)、葉片展弦比、反力度及雙圓弧葉型特征參數(shù)[17]修改葉型。采用商業(yè)軟件Numeca驗算葉型的氣動性能，根據(jù)葉片通道內(nèi)的激波、分離等流場細(xì)節(jié)不斷優(yōu)化葉型，直到3個S1流面上的氣動性能均滿足設(shè)計要求。表2和圖3為S1流面流場計算結(jié)果。由圖表可知，葉片3個流面的壓比和效率與設(shè)計要求的誤差均在1%以內(nèi)，各葉排的葉根、葉中、葉尖沒有明顯的流動分離，滿足初步設(shè)計要求。因此，該低壓壓氣機(jī)葉型可作為0D/2D耦合計算中的2D模型。

表2 低壓壓氣機(jī)S1流面計算結(jié)果Table 2 Results of the low pressure compressor on S1 stream surface

3.2 基于S1流面的低壓渦輪初步設(shè)計

與壓氣機(jī)設(shè)計類似，低壓渦輪設(shè)計也需要確定渦輪的流道、葉片數(shù)及葉型。具體途徑是：①根據(jù)設(shè)計要求(表3)及0D計算的低壓渦輪進(jìn)出口參數(shù)，得到低壓渦輪流道和中徑處的氣流角；②給定從葉根到葉中為等氣流角的分布規(guī)律；③根據(jù)氣流角，采用基于三階貝塞爾曲線[18]的渦輪葉型參數(shù)化設(shè)計方法設(shè)計二維葉型，再通過徑向積疊構(gòu)建三維葉型。

表3 低壓渦輪初步設(shè)計要求Table 3 Preliminary design requirements of the low pressure turbine

調(diào)整影響渦輪氣動性能的輪轂比、渦輪進(jìn)出口馬赫數(shù)、葉片展弦比、反力度及貝塞爾曲線5個特征點的參數(shù)，不斷改進(jìn)渦輪葉型，直到性能滿足設(shè)計要求。表4和圖4為低壓渦輪S1流面流場計算結(jié)果。由圖表可知，葉片3個流面的落壓比和效率與設(shè)計要求的誤差均在5%以內(nèi)，且各葉排的葉根、葉中、葉尖沒有流動分離現(xiàn)象，滿足初步設(shè)計要求。因此，該低壓渦輪葉型可作為0D/2D耦合計算中的2D模型。

表4 低壓渦輪S1流面計算結(jié)果Table 4 Results of the low pressure turbine on S1 stream surface

4 0D/2D耦合的渦扇發(fā)動機(jī)總體性能計算結(jié)果

4.1 設(shè)計點

根據(jù)該型渦扇發(fā)動機(jī)的要求，選取地面起飛狀態(tài)為設(shè)計點。設(shè)計點低壓轉(zhuǎn)速為11 300 r/min；低壓壓氣機(jī)進(jìn)口壓力為101 325 Pa，溫度為288.13 K，出口背壓為153 000 Pa；低壓渦輪進(jìn)口壓力為362 100 Pa，溫度為1 030 K，出口背壓為150 500 Pa。

圖5為低壓壓氣機(jī)S2流面相對馬赫數(shù)分布?？煽闯龅谝患壓偷诙夀D(zhuǎn)子有一道較強(qiáng)激波，整個流場無明顯分離。表5則顯示周向平均后的通流計算結(jié)果與設(shè)計要求誤差均在±5%以內(nèi)。

表5 低壓壓氣機(jī)S2流面計算結(jié)果Table 5 Results of the low pressure compressor on S2 stream surface

圖6為低壓渦輪S2流面相對馬赫數(shù)分布。由圖可知，整個渦輪都處在亞聲速狀態(tài)，無明顯分離。結(jié)合表6可見周向平均后的通流計算結(jié)果與設(shè)計要求誤差在±5%以內(nèi)。

表7為發(fā)動機(jī)設(shè)計點總體性能計算結(jié)果與實驗結(jié)果的對比?？梢姡?D/2D耦合計算得到的耗油率、推力等性能參數(shù)與實驗結(jié)果之間的誤差均在5%以內(nèi)。這說明該耦合模型精度滿足工程要求，可用于非設(shè)計點性能計算。

表6 低壓渦輪S2流面計算結(jié)果Table 6 Results of the low pressure turbine on S2 stream surface

表7 渦扇發(fā)動機(jī)設(shè)計點總體性能計算結(jié)果Table 7 Simulation results of turbofan engine overall performance on the design point

4.2 非設(shè)計點

非設(shè)計點計算時，假設(shè)發(fā)動機(jī)幾何不可調(diào)，通過調(diào)節(jié)燃油量控制渦輪前溫度，從而求解發(fā)動機(jī)共同工作線以及非設(shè)計點總體性能。0D/2D耦合計算采用上述方法，確定了除1.0轉(zhuǎn)速外的其他5個轉(zhuǎn)速下的共同工作點，每個工作點都至少需20次迭代。0D計算則是利用T-MATS縮放通用特性圖得到性能結(jié)果。

圖7為采用不同方法得到的低壓壓氣機(jī)共同工作線。通過對比可知：0D計算以及0D/2D耦合計算得到的喘振裕度均與實驗值很接近。圖8為非設(shè)計點總體性能對比。由圖可知：同一渦輪前溫度下，推力、燃油量、涵道比均有差異，說明三者對應(yīng)的熱力學(xué)平衡狀態(tài)不同。但相對于0D計算，0D/2D耦合的計算結(jié)果更接近實驗值。其中推力的預(yù)估精度平均提高了3.96%，涵道比的預(yù)估精度平均提高了1.20%，燃油量的預(yù)估精度平均提高了2.74%，這表明0D/2D耦合方法可提高渦扇發(fā)動機(jī)總體性能預(yù)估精度。但0D/2D耦合計算也存在誤差，這主要源于初步設(shè)計的部件和真實部件的差異。由于受壓氣機(jī)設(shè)計方法的限制，本文設(shè)計的低壓壓氣機(jī)特性線與原部件特性線有偏差(圖7)，并且該低壓壓氣機(jī)非設(shè)計點的效率均低于實驗結(jié)果，所以葉型仍有待優(yōu)化。

5 結(jié)論

(1)基于S1流面初步設(shè)計的渦扇發(fā)動機(jī)低壓壓氣機(jī)和低壓渦輪的葉型滿足要求，可用于建立0D/2D耦合模型。

(2)在T-MATS平臺上，結(jié)合基于周向平均N-S方程的通流計算程序，實現(xiàn)了0D/2D耦合方法在渦扇發(fā)動機(jī)總體性能計算中的應(yīng)用。該方法考慮了流動過程中的粘性效應(yīng)和流動損失，真實反映了低壓壓氣機(jī)和低壓渦輪內(nèi)部的流動特性，使得總體性能計算精度比0D的更高。

(3)作為數(shù)值縮放技術(shù)的簡化應(yīng)用，0D/2D耦合計算方法不僅能較為準(zhǔn)確地預(yù)估發(fā)動機(jī)的總體性能，還便于部件設(shè)計人員快速評估葉型變化對總體性能的影響，具有工程應(yīng)用價值。

參考文獻(xiàn)：

[1]Follen G，Aubuchon M.Numerical zooming between a NPSS engine system simulation and a one dimensional high compressor analysis code[R].NASA/TM-2000-209913，2000.

[2]John K L，Gregory F，Cynthia N，et al.Numerical propulsion system simulation review[R].NASA/TM-2002-209795，2002.

[3]Reed J A，Afjeh A A.Distributed and parallel program?ming in support of zooming in numerical propulsion sys?tem simulation[C]//.Proceeding of Symposium on Applica?tions of Parallel and Distributed Computing.USA：1994.

[4]陳玉春，黃 興，高本兵，等.發(fā)動機(jī)總體與尾噴管三維并行設(shè)計研究[J].航空動力學(xué)報，2007，22(10)：1695—1699.

[5]Mund F C，Doulgeris G，Pilidis P.Enhanced gas turbine performance simulation using CFD modules in a 2D repre?sentation of the low-pressure system for a high-bypass tur?bofan[J].Journal of Engineering for Gas Turbines and Pow?er，2007，129(3)：761—768.

[6]Chapman J W，Lavelle T M，May R D.Toolbox for the modeling and analysis of thermodynamic systems user’s guide[R].NASA/TM-2014-216628，2014.

[7]楊曉鋒.通流模型在組合壓氣機(jī)分析與設(shè)計中的應(yīng)用[D].南京：南京航空航天大學(xué)，2016.

[8]Simon J F.Contribution to throughflow modelling for axial flow turbomachines[D].Belgium：University of Liege，2007.

[9]Horlock J H.On entropy production adiabatic flow in tur?bomachines[J].Journal of Basic Engineering，1971，93：587—593.

[10]Bosman C，Marsh H.An improved method for calculating the flow turbo-machines,including a consistent loss model[J].Journal of Mechanical Engineering Sciences，1971，16：25—31.

[11]Hirsch C.Numerical computation of internal and external flows[M].New York：John Wiley&Sons，1988.

[12]Miller D C，Wasdell D L.Off-design prediction of com?pressor blade losses[C]//.Proceedings of the Institution of Mechanical Engineers International Conference:Turboma?chinery-Efficiency Prediction and Improvement.1987：249—260.

[13]Jameson A，Schmidt W，Turkel E.Numerical solutions of the Euler equation by finite volume method using Run?ge-Kutta time stepping schemes[R].AIAA 81-1259，1981.

[14]《航空發(fā)動機(jī)設(shè)計手冊》總編委會.航空發(fā)動機(jī)設(shè)計手冊：第7冊——渦噴及渦扇發(fā)動機(jī)總體[K].北京：航空工業(yè)出版社，2001.

[15]Wu C H.A general theory of three-dimensional flow in subsonic and supersonic turbomachines of axial-,radial-,and mixed-flow types[R].NACA TN 2604，1952.

[16]Denton J D.Multall-An open source,CFD based,turboma?chinery design system[R].ASME GT2017-63993，2017.

[17]《航空發(fā)動機(jī)設(shè)計手冊》總編委會.航空發(fā)動機(jī)設(shè)計手冊：第8冊——壓氣機(jī)[K].北京：航空工業(yè)出版社，2002.

[18]Ramanamurthy S V，Kumar S K.Development and valida?tion of a bezier curve based profile generation method for axial flow turbines[J].International Journal of Scientific Research，2014，3(12)：187.