基于0D/2D耦合的渦扇發動機總體性能計算研究

李 燕，葛 寧

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

1 引言

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

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

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

2 基于0D/2D耦合的渦扇發動機總體性能計算方法

2.1 0D的T-MATS仿真平臺

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

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

采用課題組開發的通流計算程序[7]，作為發動機部件2D數值仿真和分析的工具。程序以周向平均N-S方程作為控制方程；附加項中，采用Simon模型[8]模擬無粘葉片力，采用分布式損失模型[9-11]模擬粘性葉片力，非設計點的落后角和損失模型采用Miller[12]提出的方法。控制方程求解使用有限體積法，時間離散采用多步Runge-Kutta顯示格式[13]。該程序計算速度快，能反映真實子午流場。

2.3 基于0D/2D耦合的渦扇發動機設計點總體性能計算方法

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

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

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

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

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

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

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

2.4 基于0D/2D耦合的渦扇發動機非設計點總體性能計算方法

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

3 低壓壓氣機和低壓渦輪的初步設計

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

3.1 基于S1流面的低壓壓氣機初步設計

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

表1 低壓壓氣機初步設計要求Table 1 Preliminary design requirements of the low pressure compressor

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

表2 低壓壓氣機S1流面計算結果Table 2 Results of the low pressure compressor on S1 stream surface

3.2 基于S1流面的低壓渦輪初步設計

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

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

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

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

4 0D/2D耦合的渦扇發動機總體性能計算結果

4.1 設計點

根據該型渦扇發動機的要求，選取地面起飛狀態為設計點。設計點低壓轉速為11 300 r/min；低壓壓氣機進口壓力為101 325 Pa，溫度為288.13 K，出口背壓為153 000 Pa；低壓渦輪進口壓力為362 100 Pa，溫度為1 030 K，出口背壓為150 500 Pa。

圖5為低壓壓氣機S2流面相對馬赫數分布。可看出第一級和第二級轉子有一道較強激波，整個流場無明顯分離。表5則顯示周向平均后的通流計算結果與設計要求誤差均在±5%以內。

表5 低壓壓氣機S2流面計算結果Table 5 Results of the low pressure compressor on S2 stream surface

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

表7為發動機設計點總體性能計算結果與實驗結果的對比。可見，0D/2D耦合計算得到的耗油率、推力等性能參數與實驗結果之間的誤差均在5%以內。這說明該耦合模型精度滿足工程要求，可用于非設計點性能計算。

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

表7 渦扇發動機設計點總體性能計算結果Table 7 Simulation results of turbofan engine overall performance on the design point

4.2 非設計點

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

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

5 結論

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

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

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

參考文獻：

[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]陳玉春，黃 興，高本兵，等.發動機總體與尾噴管三維并行設計研究[J].航空動力學報，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]楊曉鋒.通流模型在組合壓氣機分析與設計中的應用[D].南京：南京航空航天大學，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]《航空發動機設計手冊》總編委會.航空發動機設計手冊：第7冊——渦噴及渦扇發動機總體[K].北京：航空工業出版社，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]《航空發動機設計手冊》總編委會.航空發動機設計手冊：第8冊——壓氣機[K].北京：航空工業出版社，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.