大型客機渦扇發動機動力特性模擬

譚偉偉，顏洪，聶智軍，馬涂亮，梁益華,*

1. 西北工業大學 航空學院，西安 710072 2. 航空工業西安航空計算技術研究所，西安 710065 3. 中國商飛上海飛機設計研究院，上海 201210

在現代大型客機設計中，機體與發動機之間的相互干擾影響對評估和改善飛機氣動性能至關重要[1-3]。20世紀80年代至今，國內外研究學者針對各種發動機進排氣氣動干擾問題，開展了理論、試驗和數值計算等方面的研究。由于發動機內部進氣、燃燒和噴氣等過程十分復雜，導致試驗研究存在難度高、周期長和費用高等問題，美國NASA(National Aeronautics and Space Administration)蘭利研究中心用試驗的方法[4]，采用渦輪動力模擬器對發動機短艙在機翼下不同安裝位置的情況進行了大量研究，以達到減少動力效應帶來的干擾阻力的目的。

隨著計算機計算能力的逐年提升，計算流體力學(Computational Fluid Dynamics, CFD)技術也逐漸成為模擬發動機進排氣氣動干擾問題的最有效手段之一[5]，Hirose[6]和Deese[7]等通過求解Euler方程，模擬了發動機進排氣效應，得到了進排氣效應會引起發動機唇口激波強度變化的結論；上海飛機設計研究院張美紅和王志棟[8]利用Ansys公司商用軟件CFX，對帶動力民用飛機進行了Navier-Stokes方程數值模擬，并應用到帶動力氣動設計中，譚兆光[9]和郝海兵[10]等建立模型，研究了機體/動力裝置一體化分析中的動力效應。模擬方法上，不論是商業軟件，還是研究機構的in-house代碼，比如FUN3D、USM3D、NSU3D、TAU、FLUENT、CFX等，都是建立發動機動力特性模型來模擬，即在進氣口設置質量流量比，排氣口和外涵道設置總溫比/總壓比來等效模擬發動機動力特性。西北工業大學劉李濤等[11]提出了一種新的渦扇發動機動力特性計算模型，并基于特征變量法構造進排氣邊界條件，最終對外吹式動力吹氣襟翼進行Navier-Stokes方程數值模擬，取得了較好的結果，但是該方法在整個邊界面網格上指定均勻的質量流量，因此會影響計算的收斂性和精度。

本文作者在實際評估大型客機動力特性工作中，分析并實現了文獻中的這些模型，并對這些模型進行了比較研究[12]，發現對大涵道比發動機，當前的這些模型計算所得發動機進排氣口質量流量(絕對值)存在較大偏差，有時相對誤差達到16%之多。因此，本文借鑒遠場特征邊界思想(保證流入、流出流場流體質量守恒)，將外涵道指定為特征邊界，建立了一種新的發動機動力特性模型，并采用標模、某型客機等構型對該模型進行了系統的驗證和確認。

1 研究方法

本文研究工作基礎為自研非結構網格并行計算軟件WoF90，是一款Euler/Navier-Stokes流動計算軟件，基于格點格式有限體積法(所有離散均在對偶網格[13-16]上進行)，可處理多種單元類型(三角形、四邊形、四面體、金字塔、三棱柱、六面體)的非結構/混合網格，支持中心格式、Roe、AUSM(Advection Upstream Splitting Method)系列以及LDFSS(Low Diffusion Flux Splitting Scheme)等多種空間離散格式，時間推進采用顯式、點/線隱式Runge-Kutta方法，同時采用多重網格、當地時間步長和殘值光順等加速收斂技術。

1.1 基本求解器

在笛卡兒坐標系下，可壓縮Navier-Stokes方程的表達式為

(1)

式中：U為流場守恒變量；FI與FV分別為無黏與黏性通量；Q為源項。將式(1)兩端積分，則可得到可壓縮Navier-Stokes方程的積分表達式為

(2)

其中:Ω為任意封閉參考區域；S為邊界面；n為邊界面單位外法向。

由式(2)可得有限體積法的數值離散方程為

(3)

式中：U0為網格頂點v0處的流場守恒變量；m0為與網格頂點v0相鄰網格點的總數；S0k為連接網格頂點v0與相鄰網格點vk(k=1,2,…,m0)的邊所對應控制面的面積；n0k為該控制面的單位法向量；FI0k與FV0k為穿過該控制面的無黏通量和黏性通量；Q0和V0分別為網格頂點v0對應控制體上的源項和控制體體積。

1.2 渦扇發動機動力特性模型

由于發動機內部風扇進氣、燃燒、噴氣等過程十分復雜，工程上評估機體/發動機干擾問題時，通常需要采用簡化模型，即在發動機進排氣口設置等效的邊界條件，建立發動機動力特性模型，來模擬發動機進排氣效應，而不直接去模擬發動機內部的復雜流場。圖1為典型民用渦扇發動機動力邊界示意圖。

在WoF90軟件中，編程實現了文獻[10]中的渦扇發動機動力特性模型A；同時分析、構造了一種新的渦扇發動機動力特性模型B。

如表1所示,分析、構造并實現的兩種渦扇發動機動力特性模型主要涉及3種邊界條件：發動機進氣口邊界條件(Nacelle_in)、發動機噴口邊界條件(Nacelle_out)以及外涵道(特征)邊界條件(Characteristic)。

圖1 渦扇發動機動力邊界示意圖Fig.1 Sketch of boundary of turbofan engine with power

表1 兩種渦扇發動機動力特性模型Table 1 Two models of turbofan engine with power

1.2.1 進氣口邊界條件

圖2為亞聲速發動機進氣道的幾何特性和進氣道流動，圖2(a)為高速飛行時捕捉到的流管；圖2(b)為起飛和低速飛行時捕捉到的流管。進入壓氣機的流管有3個主要橫截面面積(圖2(a))：前方無窮遠處的面積A∞(近似為發動機入口上游處流管剖面面積)，喉道面積A*和風扇入口面積Ainlet。

定義捕獲面積ε為流過進氣道捕獲流管上游剖面質量流量與風扇入口質量流量之比：

(4)

發動機進氣口邊界采用靜壓條件，程序根據事先設定的目標流量(通過捕獲面積ε來表達)，在計算過程中不斷調節邊界面靜壓，使實際流量盡快收斂到目標流量。

圖2 亞聲速發動機進氣道特征和進氣道流動[1]Fig.2 Characteristics of a subsonic engine intake and intake flow[1]

由等熵原理，風扇入口總壓p0,inlet等于自由來流總壓p0,∞。進而，風扇入口處靜壓為

(5)

式中：Mainlet為風扇入口馬赫數；p0,∞的表達式為

(6)

因此，計算風扇入口馬赫數Mainlet后，即可確定風扇入口靜壓。

由連續性方程ρAU=ρ*A*U*(A、U分別表示面積和速度大小，上標“*”表示聲速喉道處的物理量)可知，發動機截面面積與聲速喉道處的面積之比為

(7)

式中:T為溫度;下標“0”表示滯止參數。

根據等熵關系式：

(8)

則有：

(9)

將式(9)代入式(7)，可以直接推導出面積比與馬赫數之間的對應關系：

(10)

進一步引入拉瓦數，其定義為當地速度與臨界聲速之比，即La=U/c*,則自由來流拉瓦數和馬赫數之間的關系為

(11)

式中：R為理想氣體參數。

根據等熵關系式，可以確定自由來流拉瓦數和馬赫數之間的關系為

(12)

將式(12)代入式(10)，最終推導出面積比與拉瓦數的關系為

(13)

為了獲取風扇入口處拉瓦數Lainlet，重新定義了估算拉瓦數：

(14)

(15)

則由式(14)～式(15)可推導出下一步迭代的拉瓦數為

(16)

式中：松弛因子k可取k=0.875。通過對式(16)進行迭代求解收斂后(Δ≤1.0×10-10)，即可確定Lainlet，由拉瓦數-馬赫數關系式(12)，可計算風扇入口處的馬赫數：

(17)

進而可根據式(5)計算風扇入口處靜壓。

值得注意的是：① 進氣口邊界點上除了靜壓，其他流場變量(包括湍流變量)均通過流場內點外插得到；② 進氣口邊界條件實為流場出流條件，而下一節的噴口邊界條件實為流場入流條件。

1.2.2 噴口邊界條件

發動機噴口需給定出口處總溫與自由來流總溫之比TEx,ratio=T0,Ex/T0,∞和出口處總壓與自由來流靜壓之比pEx,ratio=p0,Ex/p∞，給定噴口的總溫比及總壓比，即可求出總溫和總壓：

(18)

利用等熵關系式，可計算噴口處馬赫數為

(19)

進而可計算得到噴口處的靜溫、靜壓和密度為

(20)

1.2.3 特征邊界條件

文獻[10]中的發動機動力特性模型(模型A)，其外涵道一般采用噴口邊界條件。本文構造的模型B，借鑒特征邊界思想，將外涵道指定為特征邊界，得到一新的發動機動力特性模型，該模型基于質量守恒定律，保證了流入、流出流場的質量守恒，因此物理上可行。

L=L(V∞,n)=

(21)

式中：V∞為自由來流原始變量；c為聲速;n=[nxnynz]T。

L-1=L-1(V∞,n)=

(22)

然后根據特征值的符號確定特征變量W1的各分量：

(23)

式中:λi為第i個特征向量對應的特征值，具體表達式為

(24)

其中:u=[uvw]。

2 計算模型驗證

為驗證本文所構造和實現渦扇發動機動力特性模型的正確性和可靠性，分別對軸對稱超高涵道比渦扇發動機模型(Counter Rotating Ultra-high-bypass Fan Simulator, CRUFS)，軸對稱渦輪動力模型(Turbine Powered Simulator, TPS)的黏性流場進行了動力特性模擬，并與試驗數據進行對比分析。在數值模擬中，主控方程對流項采用二階中心格式，黏性項采用中心差分格式離散，時間推進采用線隱Runge-Kutta迭代法，并利用多重網格、當地時間步長及殘值光順等加速收斂方法，湍流模型均采用Spalart-Allmaras(S-A)一方程模型[17]。

若無特殊說明，本文所有算例中對稱平面指定為對稱邊界條件，物面指定為絕熱壁面邊界條件，遠場指定為特征邊界條件。

2.1 CRUFS模型

CRUFS模型[18]為軸對稱超高涵道比渦扇發動機模型。由于其超高涵道比，因此在計算中通常忽略內涵道，認為由進氣道進入的氣體全部由外涵道排出[19]。

2.1.1 計算網格及狀態

CRUFS為對稱模型，可采用半模計算。網格采用ICEM CFD網格生成器制作，附面層第1層網格高度給定為5×10-6m，拉伸比為1.2，網格總點數約為273萬，遠場長度約為50倍短艙外罩長度，支桿一直延伸到遠場。發動機表面、進排氣邊界網格以及網格中心截面如圖3所示。CRUFS模型仿真參數如表2所示,表中Re/L為單位長度雷諾數。

圖3 CRUFS模型計算網格Fig.3 Computational grid of CRUFS model

表2 CRUFS模型仿真參數Table 2 Simulation parameters of CRUFS model

2.1.2 兩種動力特性模型比較

驗證渦扇發動機動力特性模型正確性的一個重要標準是：排除計算數值誤差影響后，發動機進排氣口質量流量是否一致。

首先考察收斂性，圖4給出了風扇入口處質量流量收斂歷史比較，兩種動力特性模型雖然風扇入口邊界條件一樣，但不同的外涵邊界條件導致模擬初期風扇入口質量流量有些許不同，迭代約1 800步后，兩者逐步穩定并趨于一致，很快收斂到目標值2.001 kg/s。圖中升力系數CL收斂歷史特征類似風扇入口質量流量。

圖4 CRUFS模型進氣口流量收斂速度比較Fig.4 Comparison of convergence rates of mass flow through engine fan of CRUFS model

表3 外涵道邊界類型對質量流量的影響(CRUFS)

圖5 CRUFS模型計算得到的壓力系數分布比較Fig.5 Comparison of pressure coefficient distributions calculated by CRUFS model

2.1.3 流場分析

兩種動力特性模型模擬所得流場比較接近，圖6是動力特性模型B計算得到的流場馬赫數等值線分布圖。可以看到，由于風扇進氣質量流量較大，進氣道內喉道處加速區已經超過聲速，最大馬赫數約為1.2，風扇入口中心馬赫數也達到了0.65。在短艙外表面有一加速區，但由于受風扇較大進氣質量流量的影響，加速區最大馬赫數未能超過聲速。

此外，風扇出口處中心馬赫數也在0.6附近，由于外涵道內氣流量較大且為收縮型噴管，使得馬赫數沿流向逐漸增大，且在外涵道出口附近達到最大值，約為0.85，隨著流動向下游擴散，馬赫數又逐漸降低到接近來流馬赫數0.8附近。同時，在外涵道出口處，噴管排出的高速氣流與外部氣流混合，形成了較強的剪切流動。

圖6 CRUFS模型計算得到的馬赫數等值線Fig.6 Mach number contours calculated by CRUFS model

2.2 TPS模型

TPS模型[10,20]是一軸對稱渦扇發動機模型，是日本航空宇宙技術研究所設計的風洞試驗模型。模型包括風扇入口、內涵(噴口)及外涵3個部分。

2.2.1 計算網格及狀態

TPS為對稱模型，可采用半模計算。網格采用ICEM CFD網格生成器制作，附面層第1層網格高度給定為5.8×10-6m，拉伸比為1.2，網格總點數約為273萬，遠場長度約為50倍短艙外罩長度。發動機表面、進排氣邊界網格以及網格中心截面如圖7所示，TPS模型仿真參數如表4所示。

圖7 TPS模型計算網格Fig.7 Computational grid of TPS model

表4 TPS模型仿真參數Table 4 Simulation parameters of TPS model

2.2.2 兩種動力特性模型比較

圖8給出了TPS模型風扇入口處質量流量收斂歷史比較，類似于CRUFS模型，兩種動力特性模型雖然風扇入口邊界條件一樣，但不同的外涵邊界條件導致模擬初期風扇入口質量流量有些許不同，迭代約2 000步后，兩者逐步穩定并趨于一致，很快收斂到目標質量流量8.789 kg/s。圖中升力系數收斂歷史特征類似風扇入口質量流量，但是升力系數在1 000步以后就趨于穩定。

圖8 TPS模型進氣口流量收斂速度比較Fig.8 Comparison of convergence rates of mass flow through engine fan of TPS model

正如2.1.2節中所述，驗證渦扇發動機動力特性模型正確性的一個重要標準是：排除數值誤差影響后，發動機進排氣口質量流量是否一致。表5給出了兩種動力特性模型計算得到的進排氣口質量流量(絕對值)比較。圖9給出了短艙中心截面壓力系數分布與試驗值的比較，結合表5可知，兩種動力特性模型計算所得進排氣口質量流量、短艙表面壓力系數分布一致性較好，與風洞試驗也吻合較好，說明兩種動力特性模型精度相當，可靠性好。

表5 外涵道邊界類型對質量流量的影響(TPS)

圖9 TPS模型計算得到的壓力系數分布比較Fig.9 Comparison of pressure coefficient distributions calculated by TPS model

2.2.3 流場分析

圖10給出了TPS模型中心截面處馬赫數等值線分布(動力特性模型B模擬結果)，可以看出計算得到的中心截面馬赫數分布合理，在外涵道和內涵道出口附近，短艙外氣流與外涵道氣流、外涵道氣流與內涵道氣流之間形成較強的剪切流動并且隨著流動向下游擴散，3股氣流逐漸混合，且馬赫數在下游逐步接近來流馬赫數。

圖10 TPS模型計算得到的馬赫數等值線Fig.10 Mach number contours calculated by TPS model

3 某型客機動力特性評估

為驗證本文構造的兩類動力特性模型的工程適用性，對某型客機帶翼吊渦扇發動機模型巡航狀態下有無動力、機體/發動機干擾問題進行驗證、評估，計算方法和參數與CRUFS、TPS模型一致。

3.1 計算網格及狀態

計算網格采用ICEM CFD網格生成器制作，先生成結構化網格，然后轉換為非結構網格進行計算。帶動力模型需要模擬發動機噴流，因此對發動機附近網格進行了局部修改。除發動機外，其余部件網格拓撲結構與無動力模型網格基本一致，只是對發動機部分拓撲進行了重新構造，其他位置的網格保持與無動力模型基本相同。

無動力模型為C4，網格總點數為2 589.6萬；動力模型為C4PN，網格總點數為2 577.5萬。遠場取100倍平均氣動弦長，保證遠場的氣流不干擾飛機表面的氣流，以減少對阻力的影響。機翼前緣網格間距約為0.1%平均氣動弦長。后緣分布20個網格點，機翼展向在翼根附近的網格間距約為0.8%半展長，在翼稍附近的網格間距約為0.5%半展長，機身前端和后端為1.5%平均氣動弦長。附面層第1層網格高度約為1×10-5m，拉伸比≤1.25。局部計算網格如圖11所示，某型客機仿真參數如表6所示。

3.2 兩種動力特性模型比較

與CRUFS、TPS模型一樣，首先驗證發動機進排氣口質量流量是否一致，表7給出了兩種動力特性模型計算得到的進排氣口質量流量(絕對值)的比較。從表7可知，兩種動力特性模型計算所得進排氣口質量流量一致性較好。

3.3 有無動力計算比較

工程上，評估飛行器機體/發動機干擾問題，一般先計算通氣模型無動力流場，進而模擬發動機動力特性，評估發動機動力特性對機體氣動特性的影響。因此，本節對比了有無動力對客機氣動性能的影響，圖12是有無動力模型計算所得機翼2個展向剖面上的壓力系數分布曲線(圖中b為展長)。

圖13是有無動力計算得到的展向剖面(Y=6.18 m)、法線剖面(Z=-2.12 m)馬赫數等值線分布比較圖，結合圖12可知：① 有無動力模型風扇入口處馬赫數均約為0.55；② 發動機噴流對機翼上下表面流場影響較大，使得上翼面激波位置前移、強度稍有減弱，發動機唇口激波強度亦減弱；③ 下翼面短艙外氣流與外涵道氣流、外涵道氣流與內涵道氣流之間形成較強的剪切流動并且隨著流動向下游擴散，3股氣流逐漸混合，形成了較強的剪切流動，并一直向下游延伸；④ 外涵道內氣流逐步加速，馬赫數沿流向逐漸增大，在外涵道出口附近達到最大值，約為0.92，3股氣流混合后，馬赫數又逐漸降低到接近來流馬赫數0.785附近。

圖11 某型客機計算網格Fig.11 Computational grid of civil aircraft

表6 某型客機模型仿真參數Table 6 Simulation parameters of civil aircraft model

表7 外涵道邊界類型對質量流量的影響(C4)

圖12 機翼展向剖面有無動力計算得到的壓力系數分布Fig.12 Surface pressure coefficient distribution on wing span profiles obtained by models with or without power

圖14為有無動力模型計算得到的空間渦黏性系數分布，圖中物面采用壓力系數分布渲染。翼稍小翼、后緣襟翼滑軌，機身尾部、發動機短艙均不同程度產生了尾渦，尤其是發動機噴流尾渦強度較大且一直向下游延伸，但是對平尾幾乎沒有干擾。

圖13 有無動力模型計算得到的馬赫數等值線比較Fig.13 Comparison of Mach number contours calculated by models with and without power

圖14 有無動力模型計算得到的空間渦黏性系數分布Fig.14 Eddy viscosity distribution calculated by models with or without power

4 結 論

本文以等熵流動原理為出發點，分析、構造并實現了3類發動機動力邊界條件，據此實現了文獻中的發動機動力特性模型，同時，基于特征邊界思想，建立了一種新的渦扇發動機動力特性模型。采用軸對稱超高涵道比渦扇發動機模型CRUFS、軸對稱渦輪動力模型TPS對這兩種動力特性模型進行了驗證，并采用某型客機帶翼吊渦扇發動機模型整機進行了工程評估，結果表明：

1) 兩種發動機動力特性模型風扇入口質量流量能很快收斂到目標質量流量，且均能很好地模擬渦扇發動機動力效應，計算結果與試驗數據吻合較好，驗證了本文構造的渦扇發動機動力特性模型的正確性、可靠性以及工程適用性。

2) 兩種動力特性模型計算得到的發動機進、排氣口質量流量(絕對值)相當，短艙附近機翼表面壓力系數分布一致性較好。新構造的模型B更適合于模擬大涵道比渦扇發動機動力特性，且該模型無需事先給定外涵道總溫比、總壓比等參數，在實際工程中應用更廣泛。

3) 翼吊發動機動力特性(噴流)使機翼上表面激波位置前移、強度稍有減弱，發動機唇口激波強度亦減弱，對下表面遠離短艙位置(翼根、翼稍)的影響較小。