孫 逸，葛 寧，舒 杰

（南京航空航天大學能源與動力學院，南京210016）

0 引言

在航空發(fā)動機總體設(shè)計初期，通流方法可以快速得到發(fā)動機性能以及子午面流動特征[1-2]，仍是現(xiàn)代航空發(fā)動機設(shè)計工作的基石和重要的設(shè)計工具[3-4]。Simon[5]提出基于Navier-Stokes方程的周向平均通流模型，實現(xiàn)適用于軸流壓氣機的定常與非定常通流計算，并對周向平均Navier-Stokes方程中出現(xiàn)的各種附加項做了詳細論證；金海良[6]對某跨聲速單級風扇、某子午加速風機和某跨聲速雙級風扇的研究表明，基于Navier-Stokes方程的周向平均通流模型計算得到的子午流場和徑向參數(shù)分布與3維計算結(jié)果十分吻合，在近端壁區(qū)域具有較高的計算精度；周向平均通流模型的另一重要應(yīng)用是發(fā)動機整機數(shù)值仿真計算，可用于發(fā)動機各部件間的流場匹配，研究發(fā)動機各部件間的相互影響[7]；Friederike[8]針對某大涵道比渦扇發(fā)動機提出1種0D/2D耦合航空發(fā)動機性能預(yù)估方法；Denton等[9]基于40多年研究經(jīng)驗開發(fā)了多層次葉輪機械設(shè)計程序Multall，其中通流計算作為重要環(huán)節(jié)，可快速分析葉片損失和效率。由于現(xiàn)階段受CFD技術(shù)和計算機性能的限制，采用3維數(shù)值模擬方法計算發(fā)動機整機需要大量的計算資源和時間，而采用通流計算方法能夠節(jié)約計算時間，快速發(fā)現(xiàn)各部件匹配過程中的問題，具有一定的工程應(yīng)用價值。本文選取某渦扇發(fā)動機低壓部件為研究對象，采用通流模型耦合計算，研究各部件相互作用以及過渡段幾何參數(shù)對發(fā)動機總體性能的影響。通流模型可用于發(fā)動機總體設(shè)計初期對發(fā)動機性能快速預(yù)估，實現(xiàn)各部件性能匹配。

1 基于Navier-Stokes方程的周向平均通流模型

在相對柱坐標系中，對Navier-Stokes方程進行周向平均[10]，得到通流模型的控制方程為

其中

式中：Q為守恒量；F和H為對流（無黏）通量；Fv和Hv為擴散（黏性）通量；S為Navier-Stokes方程組在相對柱坐標系下導(dǎo)出的源項；SB為無黏葉片力項；Sv為黏性葉片力項；fB和fv分別為?；蟮臒o黏葉片力和黏性葉片力；x、r、θ分別為軸向、徑向和周向坐標方向；ρ為密度；v為相對速度；p為靜壓；e為單位質(zhì)量總能；w為轉(zhuǎn)速；τ為黏性應(yīng)力；q為熱流量；定義葉片堵塞系數(shù)

式中：θs和θp分別為葉片吸力面和壓力面上的角坐標；N為葉片數(shù)。

堵塞系數(shù)b在葉片區(qū)小于1，在非葉片區(qū)等于1?；跓o黏葉片力始終垂直于平均流面的假設(shè)，本文將無黏葉片力的求解分為2步[11-12]：（1）在動量方程中不計入無黏葉片力；（2）在動量方程中僅計入無黏葉片力，從而可以修正第1步求得的動量，使其與平均流面相切。這就避免了顯式求解無黏葉片力。黏性葉片力的?；捎弥姆植紦p失模型[13]，引入1種分布的體積力來計入流動損失的影響。黏性葉片力平行于相對速度并且反向，所作功僅產(chǎn)生熵增，采用總壓損失系數(shù)計算葉排進口到出口熵增。本文采用時間推進的有限體積法求解上述控制方程，用具有較高間斷分辨率的Roe通量差分分裂方法[14]對無黏對流通量進行離散，而黏性通量的離散則是利用網(wǎng)格中心的流場變量將網(wǎng)格面上的黏性通量表示出來，紊流模型采用B-L模型。

2 通流模型驗證

通流耦合計算中首先要保證通流模型的準確性，激波位置和激波強度對壓氣機流量、效率均有較大影響。因此，選取NASA Rotor67和某3級風扇為驗證算例。

2.1 Rotor67算例

Rotor67是1個低展弦比跨聲速轉(zhuǎn)子，Strazisar[15]于1989年做了細致的試驗。設(shè)計轉(zhuǎn)速為16043 r/min，設(shè)計轉(zhuǎn)速下最高效率點流量為34.573 kg/s，壓比為1.642，效率為0.93。最高效率點通流模型計算結(jié)果與試驗值對比見表1，可見通流計算結(jié)果與試驗值基本吻合。根據(jù)Strazisar的描述，由于在轉(zhuǎn)子出口邊界層處測點較少，試驗測量的等熵效率可能偏高。

表1 Rotor67通流計算與試驗結(jié)果對比

Rotor67在最高效率點時，通流模型計算得到的葉排出口截面參數(shù)（分別為氣流角、總溫、靜壓、總壓）徑向分布與試驗值對比如圖1所示。從圖中可見，通流計算得到的靜壓分布與試驗值十分吻合，但由于葉尖存在激波，氣流角和總溫在葉尖處與試驗值有一定偏差?？倝簭较蚍植寂c試驗基本吻合，而且能夠預(yù)測出端壁附面層附近的總壓變化趨勢。

圖1 Rotor67葉排出口徑向參數(shù)分布

Rotor67在最高效率點的馬赫數(shù)等值線如圖2所示，從圖中可見，葉片區(qū)域存在2道激波，第1道為葉尖進口斜激波，第2道為尾緣附近的正激波，波前馬赫數(shù)約為1.4。通流計算得到的Rotor67在100%轉(zhuǎn)速下的特性線如圖3、4所示。從圖中可見，二者結(jié)果基本吻合，驗證了通流模型的準確性和可靠性。

圖2 Rotor67最高效率點馬赫數(shù)等值線

圖3 流量-效率特性線

圖4 流量-壓比特性線

2.2 某軸流風扇算例

選取某帶進口導(dǎo)葉的3級軸流風扇進行驗證，其設(shè)計狀態(tài)下壓比為3.5，流量為81.0 kg/s。風扇子午流場相對馬赫數(shù)如圖5所示。從圖中可見，第1級轉(zhuǎn)子葉尖最高波前馬赫數(shù)約為1.5，較強的激波使得葉尖壓比相對試驗值偏高，同時較強的激波-邊界層相互干擾損失導(dǎo)致轉(zhuǎn)子葉尖效率較低。

圖5 風扇子午流場相對馬赫數(shù)

計算了該軸流風扇在100%轉(zhuǎn)速下的性能特性，并與試驗值進行了對比，流量-效率特性線與流量-壓比特性線分別如圖6、7所示。從圖中可見，在風扇近設(shè)計點計算得到的壓比為3.53，效率為0.84。由于周向平均通流模型捕捉到的激波其物理本質(zhì)與3維計算有所區(qū)別，因此捕捉到的激波比實際稍強，導(dǎo)致計算結(jié)果效率偏低。通流模型計算得到的特性線與試驗值基本吻合，可用于多級壓氣機特性計算。

圖6 風扇流量-效率特性

圖7 風扇流量-壓比特性

3 通流耦合計算

3.1 研究對象

選取某渦扇發(fā)動機除核心機3大部件外的低壓部件為研究對象，包括進口帶導(dǎo)葉的3級風扇、出口帶支板的2級低壓渦輪、內(nèi)外涵道、摻混段以及尾噴管。其流道和葉型根據(jù)該發(fā)動機總體及部件性能參數(shù)試驗值（見表2），采用流道設(shè)計程序計算得到。子午面流道如圖8所示。風扇出口分為內(nèi)外涵道2個計算域，并通過外涵道與渦輪后摻混段相連，摻混段下游為尾噴管。內(nèi)、外涵冷熱流體摻混是1個重要的氣動熱力過程，本文采用流量加權(quán)后定比熱方法計算摻混后流場。耦合計算邊界條件均按試驗值給定，包括風扇進口總溫、總壓，內(nèi)涵出口靜壓，低壓渦輪進口總溫、總壓，噴管出口給定標準大氣壓力。

表2 通流計算結(jié)果與試驗值對比

圖8 幾何模型子午面

3.2 計算結(jié)果及分析

計算結(jié)果表明，各部件單獨計算時性能均能滿足總體設(shè)計要求。但多個部件耦合數(shù)值計算時，發(fā)動機總體性能參數(shù)往往會偏離設(shè)計狀態(tài)。本文通流耦合計算的核心就是研究各部件匹配以及過渡段對發(fā)動機總體性能的影響。

研究發(fā)現(xiàn)，風扇出口內(nèi)、外涵分流板、渦輪出口摻混段整流板是影響部件匹配和發(fā)動機總體性能最重要的2個過渡段幾何參數(shù)。分流板不僅分配內(nèi)、外涵流量，控制涵道比，并且對風扇的工作狀態(tài)有較大影響。摻混段整流板對內(nèi)、外涵摻混有一定影響，而且可以抑制尾椎流動分離，對渦輪和噴管工作狀態(tài)也至關(guān)重要。

通過調(diào)節(jié)分流板和整流板幾何參數(shù)，使得各部件性能結(jié)算結(jié)果接近設(shè)計狀態(tài)。通流計算子午流場相對馬赫數(shù)如圖9所示。從圖中可見，風扇轉(zhuǎn)子葉尖均存在激波，第1級轉(zhuǎn)子中波前馬赫數(shù)約為1.4，較強的激波損失使得風扇效率計算結(jié)果偏低3.57%。

圖9 子午流場相對馬赫數(shù)

從表2列出的風扇、低壓渦輪、內(nèi)外涵以及尾噴管近設(shè)計點性能參數(shù)與試驗值對比可見，風扇計算流量偏大3.14%，增壓比偏小0.56%，效率偏低3.57%。低壓渦輪計算流量偏大4.93%，落壓比偏大0.92%，效率偏低3.41%。通流模型得到的內(nèi)涵流量偏大2.28%，外涵流量偏小0.48%，導(dǎo)致涵道比偏小2.55%。噴管出口計算流量和排氣速度比試驗值偏大，導(dǎo)致發(fā)動機計算推力比試驗值偏大2.67%，耗油率偏低2.52%。從誤差分析來看，本文通流耦合計算在近設(shè)計點能夠滿足發(fā)動機總體性能初步要求，最大誤差不超過5.0%。

4 結(jié)論

本文通過對某渦扇發(fā)動機低壓部件通流耦合計算，得到如下結(jié)論：

（1）采用周向平均通流模型對Rotor67和某3級風扇進行了計算并與試驗值對比，結(jié)果證明通流模型能夠正確捕捉激波強度與位置，預(yù)測的部件特性曲線趨勢與試驗基本吻合。

（2）對某渦扇發(fā)動機風扇、低壓渦輪、內(nèi)外涵道、摻混段和尾噴管等部件進行了通流耦合計算，研究了各部件匹配與過渡段幾何參數(shù)對發(fā)動機總體性能的影響，在設(shè)計狀態(tài)下發(fā)動機沿程各截面計算特征參數(shù)與試驗值最大誤差不超過5.0%。

（3）對比分析發(fā)現(xiàn)，風扇出口分流板控制涵道比并對風扇性能有一定影響。渦輪出口摻混段整流板對風扇和低壓渦輪性能均有較大影響。通過調(diào)節(jié)過渡段幾何參數(shù)，實現(xiàn)各部件性能匹配，達到發(fā)動機總體性能要求。

該通流模型能夠快速預(yù)估發(fā)動機總體性能，發(fā)現(xiàn)各部件匹配中存在的問題，具有一定的工程應(yīng)用價值。內(nèi)、外涵摻混段涉及變比熱問題，本文采用內(nèi)、外涵流量加權(quán)平均比熱比，噴管出口總溫平均值比試驗值偏高3.5%。摻混段計算有待進一步完善，使噴管溫度場與實際情況更加吻合。

[1]Smith L H.The radial-equilibrium equation of turbomachinery[J].ASME Journal of Engineering for Power，1966，88（1）：1.

[2]Jin H L，Jin D H，Li X J.A time-marching throughflow model and its application in transonic axial compressor [J].Journal of Thermal Sinence，2010，19（6）：519-525.

[3]Denton J D，Horlock J H.A review of some early design practice using computational fluid dynamics and a current perspective[J].Journal of Turbomachinery，2005，127（1）：5-13.

[4]Denton J D，Dawes W N.Computational fluid dynamics for turbomachinery design[J].Archive Proceedings of the Institution of Mechanical Engineers，1998，213（2）：107-124.

[5]Simon J F，Leonard O.Modeling of 3D losses and deviations in a throughflow analysis tool[J].Journal of Thermal Science，2007，16（3）：208-214.

[6]金海良.周向平均方法在多級軸流風扇/壓氣機設(shè)計與分析中的應(yīng)用[D].北京：北京航空航天大學，2011.JIN Hailiang.Application of circumferentional average method in multistage axial fan/compressor design and analysis[D].Beijing：Beijing University of Aeronautics and Astronautics，2011.（in Chinese）

[7]Wan K，Jin H，Jin D.Influence of nonaxisymmetric terms on circumferentially averaged method in fan/compressor[J].Journal of Thermal Science，2013，22（1）：13-22.

[8]Mund F C，Doulgeris G，Pilidis P.Enhanced gas turbine performance simulation using CFD modules in a 2D representation of the low-pressure system for a high-bypass turbofan [J].Gas Turbines Power,2006，129（3）,761-768.

[9]John D Denton.Multall-an open source，CFD based，turbomachinery design system[R].ASME 2017-GT-63993.

[10]Baralon S，Eriksson L E，Hall U.Evaluation of higher-order terms in the throughflow approximation using 3D Navier-Stokes computations of a transonic compressor rotor[R].ASME 1999-GT-074.

[11]Gentry R A，Martin R E，Daly B J.An eulerian differencing method for unsteady compressible flow problems[J].Journal of Computational Physics，1966，1（1）：87-118.

[12]寧方飛.考慮真實幾何復(fù)雜性的跨音壓氣機內(nèi)部流動的數(shù)值模擬[D].北京：北京航空航天大學，2002.Ning F F.Numerical investigations of flows in transonic compressors with real geometrical complexities[D].Beijing：Beijing University of Aeronautics and Astonautics，2002.（in Chinese）

[13]Bosman C，Marsh H.An improved method for calculating the flow in turbomachines，including a consistent loss model[J].Archive Journal of Mechanical Engineering Science，1974，31（16）：25-31.

[14]Roe P L.Approximate riemann solvers，parameter vectors，and difference schemes.[J].Journal of Computational Physics，1981，43（2）：357-372.

[15]Strazisar A J，Wood J R，Hathaway M D.Laser anemometer measurements in a transonic axial-flow fan rotor[J].NASA STI/Recon Technical Report N，1981，90（2）：430-437.