基于通流方法的某渦噴發動機整機數值仿真

劉曉恒，周成華，宋滿祥，金東海,2,*，桂幸民,2

1. 北京航空航天大學 能源與動力工程學院，北京 100083 2. 北京航空航天大學 先進航空發動機協同創新中心，北京 100083

隨著計算流體力學的發展以及計算機計算能力的提升，全三維Navier-Stokes數值計算越來越多地應用到葉輪機械的設計過程中。但是，這需要性能極高的計算設備以及大量的計算時間，以至于無法將發動機整機全三維數值仿真計算方法應用到發動機整機日常設計過程中[1-4]。另一方面，零維計算也經常應用在航空發動機整機性能計算中，它依賴于各部件的性能特性圖[5]。但是，零維計算無法提供氣動參數的展向分布結果。因此，通流方法在現代燃氣輪機和航空發動機設計過程中仍然是非常重要的一部分，尤其是在初始設計階段[6]。這種方法不僅可以應用在給定氣動參數的設計階段，還可以應用在給定幾何結構分析其性能參數的分析階段。

子午面通流方法的概念應當追溯到吳仲華的研究[7-8]，借助準三維流面概念，推導出葉片間流動以及子午面流動的統一模型。通過研究兩組相對流面，三維流動的正問題與反問題均可以得到求解。第1套使用Euler方程的通流程序是由Spurr[9]提出的，他利用一種通流模型以及葉片間流動的求解器建立一種準三維程序，這種方法利用時間推進求解Euler方程。將這種程序應用到跨聲導向器案例中，計算結果與三維求解Euler方程計算結果相一致。Yao和Hirsch[10]提出一種方法，將三維計算程序轉換為Euler通流方法，他們研究的目的是將二維設計方法與三維分析方法相結合，并且應用了損失系數和落后角的經驗關系式。Dawes[11]提出一種方法，對于多級部件計算，可以將通流方法與全三維數值模擬相結合，將所研究部件利用三維黏性計算進行模擬，同時其他部件利用軸對稱方法進行建模。通過這種方法，單級葉片排可以在多級的環境下進行設計和分析。在Damle等[12]的研究中，描述了利用Euler模型設計跨聲/超聲葉輪機械的方法，這種Euler通流方法成功應用在一個超聲風扇與一級跨聲渦輪部件的設計過程中。

通流方法的另一個應用領域，是在發動機整機的分析過程和初始設計階段。由Ivanov等[13-15]提出的通流模型是應用于發動機分析問題中的，他們闡述了在通流計算中描述發動機工作過程的現代數學模型。該模型使用單調保持的二階Godunov格式，這種格式可以準確求解Riemann問題。在Petrovic等[16-17]的研究中，提出一種基于分別求解壓氣機與渦輪的求解程序進而求解整機流場的方法。利用簡單的燃燒室模型以及二次流動模型來連接壓氣機與渦輪流道，它可以將發動機整機進行自動匹配?；谝陨夏Ｐ蛯σ慌_雙軸燃氣輪機進行整機數值模擬，獲得了發動機流場。

國內相關單位對發動機整機計算進行了初步的研究。施發樹和劉興洲[18]、黃家驊等[19]對小型雙涵道渦扇發動機進行二維穩態數值模擬，采用Godunov格式求解帶黏性力項的非定常Euler方程組。馮國泰等[20]討論了適用于發動機葉輪部件多場耦合數值仿真計算的統一數學物理模型，分析了包括發動機三維多功能數值仿真數學模型、精度可靠性與并行算法、發動機仿真軟件平臺的框架等在內的6個關于建立航空發動機數值仿真試驗臺的關鍵技術問題。曹志鵬等[21]借助二維仿真軟件對某小型渦噴發動機設計點進行數值仿真，對其性能參數及流場進行分析。萬科等[22]借助周向平均降維方法對某高通流風扇/增壓級進行了性能分析，并與三維數值模擬結果進行對比分析，為航空發動機整機通流計算奠定了基礎。

借助通流方法不僅可以獲得葉片內部的流場信息，同時可以獲得葉片間的匹配信息，這在發動機初始設計階段是非常重要的。本文利用通流計算軟件，對某渦噴發動機設計點與非設計點進行性能與流場研究。首先，對發動機地面靜止狀態節流特性進行計算，將計算結果與試驗數據和設計參數進行對比分析，以此來驗證該通流軟件計算的可靠性；其次，為了探索發動機的高度/速度特性，進行了相關工況的數值仿真，將計算結果與設計值進行對比；再次，對壓氣機/渦輪進行部件仿真，獲取了發動機共同工作線；最后，對發動機設計工況下的流場參數以及展向分布進行了分析。

1 通流計算方法

1.1 控制方程組

對于柱面坐標系(z,r,φ)下帶有黏性力的非定常Euler方程組，利用貼體坐標系(ξ,η,ζ)進行坐標轉換，可以得到如下形式的控制方程組：

(1)

式中：

利用狀態方程將式(1)封閉，其表達式為

p=ρRT

(2)

(3)

(4)

式中：R為氣體常數；T為溫度；ε為比內能；T0為參考溫度；cv(τ)為定容比熱容；Const表示常數。

對于葉輪機械，可以將計算域分為兩類：包含葉片區域和不含葉片區域。

對于不含葉片區域(如葉片排間隙和燃燒室區域)，可認為流動是軸對稱的，因此存在如下假設：

(5)

借此，可將控制方程組式(1)推導為

(6)

對于存在葉片的計算區域，可以認為ζ=Const定義的曲面為相鄰葉片間的流面，這樣存在如下關系：

(7)

基于此，可將坐標(ξ,η)選取為滿足如下關系式的參數：

(8)

最終，應用于包含葉片區域的控制方程推導為

(9)

式中：

由于航空發動機各部件空間幾何結構的復雜性，為了保證計算域邊界處的計算精度，本程序采用貼體坐標系對控制方程進行轉換求解[23]。關于曲線坐標系以及軸對稱降維方法的詳細介紹可以參考文獻[15]。

在數值模擬過程中，需要給定以下邊界條件：轉速，進口總溫/總壓，大氣溫度/壓力，出口背壓，燃油流量，燃油燃燒熱值以及對冷卻空氣提取和排放位置的設定。在計算過程中，對流動過程中真實的三維流動細節進行建模模擬，比如：黏性損失，間隙溢流，流道中通過封嚴裝置的泄漏流，冷卻空氣的提取以及排放。

1.2 通流模型

整機仿真程序求解的是二維歐拉方程, 因此損失及落后角模型是決定計算精度的關鍵因素。

包含在式(1)中的(fz,fr,fφ)是表征方程黏性力的3個參數：

(10)

式中：Φ用以模擬黏性損失。

利用σ描述由于黏性損失導致的熵增，則存在如下關系：

(11)

函數σ隨空間變化，并且與熵變化是同步的，損失值的大小由試驗數據或者經驗關系式計算得到。

根據部件結構不同，葉輪機械損失模型可以分為軸流壓氣機、離心壓氣機、軸流渦輪及徑向渦輪等類型。在仿真程序中，損失模型可以分為設計點損失模型和非設計點損失模型。在設計點，總損失由葉型損失、激波損失以及其他損失構成，前面二者可以準確計算得到，其他損失需要利用損失系數進行修正得到。在非設計點,除了設計總損失之外,還需要確定由攻角、端壁、尾緣等引起的附加損失修正系數,與設計損失一起構成非設計點損失。最終，估算得到的損失會轉化為熵增，計入黏性葉片力中。對于軸流壓氣機，其使用模型包括：Koch-Smith葉型損失模型[24]、Creveling非設計點損失模型[25]、Howell端壁損失模型[26]、Aungier葉尖間隙泄漏損失模型[27]；對于離心壓氣機，使用Galvas損失模型[28]；對于軸流渦輪，使用AMDCKO損失模型[29]。

程序中落后角模型用于計算設計工況下落后角，在非設計工況需要利用系數對落后角進行修正。對于設計工況和非設計工況，分別采用Carter落后角模型[30]和Creveling落后角模型[25]。

2 發動機結構以及試驗系統

本文所研究發動機模型在子午平面內的簡圖如圖1所示，其結構包括：一級軸流壓氣機、一級離心壓氣機、折流燃燒室以及一級軸流渦輪。

測試發動機臺架如圖2所示，臺架采用支撐式搖臂結構，它具有2個底座支撐點以及1個具有半圓形結構的懸掛裝置，運動的試驗臺架通過2根軸和4個搖臂與固定的試驗臺架相連。當發動機產生推力時，移動的臺架向前移動，然后與固定的臺架相撞，固定臺架上的壓力傳感器即可測量得到發動機產生的推力。

借助該試驗系統，可以測量得到如下數據：推力、燃油流量、進口總壓、離心壓氣機出口總壓以及渦輪后的總溫。利用這些試驗數據，可以對地面工況數值仿真結果進行對比分析進而考核準三維數值仿真軟件的計算準確性。

圖1 渦噴發動機的幾何結構Fig.1 Geometry of turbojet engine

圖2 測試發動機與測試臺架Fig.2 Test engine and bench

3 計算結果分析與討論

3.1 地面靜止狀態節流特性分析與校核

首先，對航空發動機的地面節流特性進行計算，并將其與試驗數據以及初始設計值進行對比，以考核軟件計算的精度。需要說明的是，該軟件計算結果是振蕩收斂的，對于收斂結果采用算術平均處理方法，將參數標準差作為誤差帶進行作圖。圖3(a)～圖3(d)分別為發動機推力、單位燃油消耗率、組合壓氣機壓比以及渦輪后總溫隨轉速的變化規律(圖中給出了試驗數據和計算結果的方差)。其中，試驗過程中發動機選取4個相對轉速作為穩定狀態，分別是80%、90%、95%(額定轉速)、98%，因此對這4個工況進行考核分析。

如圖3(a)所示，當轉速增加時，發動機的推力隨之增加。通過計算結果與試驗數據和設計值對比發現，3組數據的變化趨勢是一致的。與試驗數據對比，最大相對誤差是-5.1%；與設計值進行對比，最大相對誤差是-3.6%。如圖3(b)所示，當轉速增加時，發動機的單位燃油消耗率隨之降低。與試驗數據對比，最大相對誤差是+4.8%；與設計參數進行對比，最大相對誤差是+4.7%。另外，在95%轉速工況下，發動機的單位燃油消耗率最低，這也說明了飛機處于巡航狀態時為何選擇95%轉速作為巡航速度。

圖3 地面節流特性Fig.3 Throttle characteristics on ground

圖3(c)和圖3(d)展示了壓氣機壓比和渦輪后總溫隨轉速的變化曲線，隨著發動機轉速的增加，壓氣機的壓比以及渦輪后總溫都隨之增加。同時可以看到，對于壓氣機壓比，計算結果相較于試驗數據較高，這導致在計算結果中渦輪的輸出功率較高，因此，對于渦輪后總溫，計算結果低于試驗數據。

通過上述分析，驗證了該計算方法對發動機模型計算的準確性，進而可以對其他工況進行數值模擬。

3.2 高度/速度特性

航空發動機的飛行試驗耗費巨大，為了研究發動機的高度/速度特性，利用該通流方法計算得到了發動機在不同高度/速度的特性參數，并將計算結果與初始設計參數進行了對比。

3.2.1 高度特性

進行發動機高度特性數值模擬時，工況設置為：相對轉速為100%，飛行馬赫數為0.7，飛行高度分別為3、6、9、12、15、18 km。圖4(a)為發動機的相對推力與高度的關系曲線，圖4(b)為發動機的單位燃油消耗率與高度的關系曲線(計算結果均利用100%轉速工況下設計參數進行了無量綱化)。通過兩組數據可知：計算結果與設計參數的趨勢是一致的。對比計算結果與設計值，推力的最大相對誤差位于高度3 km處，數值為-4.61%；單位燃油消耗率的最大相對誤差出現在同一高度，數值為+5.0%。

對于曲線變化趨勢的解釋如下：當飛行高度低于11 km時，隨著飛行高度的增加，溫度、壓力以及密度均降低；同時，相較于溫度，壓力、密度降低得更快。因為飛行馬赫數是給定的，所以當高度增加時，進口流量降低。發動機的推力同時取決于流量和單位推力，流量起主要作用。因此，發動機推力隨飛行高度的增加而降低。

當討論發動機的單位燃油消耗率時，將發動機的單位推力作為一個考慮因素。其中，單位燃油消耗率的定義為

(12)

式中：f為油氣比；Fs為發動機的單位推力。隨著飛行高度增加，當高度低于11 km時，大氣溫度隨之降低，當飛行高度高于11 km時，大氣溫度維持不變。因為發動機物理轉速維持不變，在飛行高度低于11 km時，折合轉速隨高度增加而增加，當高度大于11 km時，折合轉速不再發生變化。結合壓氣機特性可知，隨著高度的增加，壓氣機壓比先增大然后維持不變。假設渦輪導向器以及尾噴管均處于堵塞狀態，則渦輪的膨脹比維持不變。因此，在飛行高度增加的過程中，發動機出口的靜壓先升高后維持不變，產生的單位推力先增大后不變，單位燃油消耗率先降低然后維持不變。

圖4 高度特性Fig.4 Altitude characteristics

3.2.2 速度特性

為了研究發動機的速度特性，在不同飛行馬赫數工況下，對發動機進行數值模擬。計算工況設置為：物理轉速保持不變，飛行高度為3 km，飛行馬赫數的變化范圍為0.2～0.8，間隔設置為0.1。圖5(a)為發動機推力與飛行馬赫數的關系曲線，圖5(b)為單位燃油消耗率與飛行馬赫數的關系曲線。由圖可知計算結果與設計參數兩組數據的變化趨勢是一致的；對于計算精度，推力和單位燃油消耗率的最大相對誤差分別是-5.83%和+5.92%，并且二者最大值均出現在馬赫數為0.8時。

呈現該變化趨勢的原因如下：一方面，當飛行高度保持不變時，隨著飛行馬赫數的增加，進氣道空氣流量隨之增加；另一方面，進口總溫的增加使發動機的折合轉速降低，壓氣機的工作點在共同工作線上向低壓比方向移動，進而發動機出口的壓力降低，單位推力降低。在流量與單位推力這兩個因素的共同作用下，發動機的推力隨著飛行馬赫數的增大首先降低然后升高。根據式(12)，隨著單位推力的減少，發動機的單位燃油消耗率增加。

圖5 速度特性Fig.5 Velocity characteristics

3.3 共同工作線預測

借助該通流方法可以首先對發動機壓氣機/渦輪進行部件計算，進而獲得發動機的共同工作線，發動機轉速的變化范圍為75%～100%。在進行壓氣機和渦輪匹配時，利用流量和功率兩個參數。首先，在各個轉速，計算得到壓氣機的特性曲線；然后，計算得到渦輪的特性曲線。相同轉速下壓氣機與渦輪特性曲線的交點即為二者具有相同流量和功率的工作點，也就是發動機的共同工作點。需要說明的是，在進行功率匹配時，渦輪的輸出功率有一部分用于克服渦輪軸的損失，該部分占渦輪輸出總功率的1%。得到發動機共同工作點的流量和功率后，對應得到無量綱密流和組合壓氣機壓比，即可在壓氣機特性線上繪制發動機共同工作線，如圖6所示。

由圖6可知：對于無量綱密流，計算結果與設計參數相對誤差最大值出現在相對轉速為75%位置，相對誤差為+5.4%；對于組合壓氣機壓比，計算結果與設計參數相對誤差最大值在相對轉速為80%位置，相對誤差為+9.4%?？紤]到這種計算方法的便捷性，相對誤差是可以接受的。另一方面，在75%轉速處密流和組合壓氣機壓比的相對誤差相對于其他轉速較高的原因是：在遠離100%轉速的工況，計算使用的損失和落后角模型相較于100%轉速處不夠準確。

圖6 發動機的共同工作線Fig.6 Workline of engine

3.4 整機流場分析

相較于零維/一維計算，通流計算方法不僅可以獲得發動機的性能參數，同時可以得到發動機主流道內各位置氣動參數，進而獲取氣動參數的展向分布，為發動機部件匹配提供更詳細的信息。本部分對發動機整機地面設計點工況進行計算，并對流場參數進行分析。計算機具有8個處理器，處理器型號為Intel(R) Core(TM) i7 CPU 870@2.93 GHz。整機計算網格數量為8 195，完成200 000步迭代計算耗時為2 h。圖7給出了發動機整機計算流量和推力迭代史，從圖中可知，計算結果呈現振蕩收斂，因此，需對氣動參數進行算術平均處理。

圖7 整機計算參數迭代史Fig.7 Iteration history of parameters in simulation of engine

借助100%轉速設計工況下渦輪功率，對計算所得壓氣機和渦輪功率進行無量綱化，并對振蕩結果進行算術平均可得：組合壓氣機功率為0.960 2，渦輪功率為0.975 3，二者差值占渦輪輸出功率的1.55%，因此認為壓氣機與渦輪功率平衡。表1給出了設計工況下性能參數的相對誤差，其中相對誤差的定義為：(計算結果-設計參數)/設計參數。通過表中數據可知，各部件流量的相對誤差均小于5%，各部件壓比的相對誤差均小于8%，基于這些數據可以認為計算結果是合理的。

表1 性能參數對比Table 1 Comparison of performance parameters

圖8為氣動參數的流場等值線圖。其中，圖8(a)為相對馬赫數的等值線圖，觀察可知，在軸流壓氣機轉子葉片尖部區域，存在馬赫數大于1的區域，可以證實該軸流壓氣機轉子葉片是跨聲轉子。在渦輪進口導流葉片中，同樣存在馬赫數大于1的區域，這與渦輪葉片設計理念是一致的：超聲區域存在于渦輪進口導葉出口根部。在軸向擴壓器出口，存在低馬赫數區域。圖8(c)為總溫的等值線圖，可以看到經過風斗后氣體的總溫降低，這與在折流燃燒室中借助風斗引入摻混冷氣起到冷卻高溫氣體的作用是一致的。

在獲得了發動機全流場的氣動參數后，即可獲得發動機在某一計算站氣動參數的展向分布。例如：借助轉子進出口氣動參數得到轉子葉片的效率；借助靜子葉片進出口的壓力可以獲得靜子葉片的總壓恢復系數。對于整級葉輪機械而言，可以獲得其級參數(總壓升/降、熱力反力度等)。這些級參數的定義為

轉子葉片排的效率：

(13)

熱力反力度：

(14)

式中：下標1表示級進口位置；下標2表示級葉片排軸向間隙位置；下標3表示級出口位置；下標is表示對應的溫度為絕熱溫度。另外，式(14)中所有的溫度均為靜溫。

以軸流壓氣機級為例，計算得到氣動參數的展向分布如圖9所示。由圖9可知：① 轉子葉片效率及靜子葉片總壓恢復系數呈現出葉尖較低、葉中較高的現象，并且二者在葉片根部10%展高范圍內均存在降低的趨勢；② 對于級壓比及熱力反力度，它們呈現出隨葉高增加而升高的趨勢。借助這些參數，設計人員即可了解在不同位置葉片的設計結果，同時借此得到相鄰葉片排的匹配效果，在發動機的初始設計階段這些參數具有重要意義。

圖8 氣動參數等值線圖Fig.8 Contours of aerodynamic parameters

圖9 氣動參數的展向分布Fig.9 Spanwise distribution of aerodynamic parameters

4 結 論

基于軸對稱通流方法，對某渦噴發動機設計點及非設計點進行了整機數值仿真，對其性能及流場參數進行了對比分析，得到如下結論：

1) 對發動機地面靜止狀態節流特性進行數值仿真，發動機轉速設定為80%、90%、95%及98%，獲得了推力、單位燃油消耗率、壓氣機壓比以及渦輪后總溫隨轉速的變化規律，并將計算結果與發動機整機試驗數據進行對比。發動機推力的最大相對誤差為-5.1%，單位燃油消耗率的最大相對誤差為+4.8%。另外，計算得到在95%轉速時發動機單位燃油消耗率最低，這也說明了飛機在巡航狀態下為何選擇95%轉速作為巡航速度。

2) 對發動機在飛行馬赫數0.7處進行高度特性數值模擬，將計算結果與設計參數對比，推力的最大相對誤差為-4.61%，單位燃油消耗率的最大相對誤差為+5.0%，二者均位于3 km處。對發動機在3 km位置進行速度特性模擬，推力和單位燃油消耗率的最大相對誤差分別是-5.83%和+5.92%，二者均位于馬赫數為0.8時。

3) 分別進行壓氣機/渦輪部件模擬，利用流量和功率參數進行整機匹配，獲取發動機的共同工作線，發動機轉速范圍為75%～100%。將計算結果與設計值進行對比，各轉速流量相對誤差小于6%，組合壓氣機壓比相對誤差小于10%。

4) 較于零維/一維計算，通流計算除了獲取性能參數外，還可以得到發動機流場。對軸流壓氣機各氣動參數的展向分布特征進行分析，可為部件匹配提供技術支持。