賈洪印，鄧有奇，馬明生，張耀冰

（中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所，四川 綿陽(yáng) 621000）

0 引言

飛機(jī)／發(fā)動(dòng)機(jī)一體化是現(xiàn)代飛機(jī)設(shè)計(jì)中一個(gè)重要方面，準(zhǔn)確模擬和預(yù)測(cè)機(jī)體與動(dòng)力裝置之間的相互干擾影響，對(duì)于評(píng)估和改善飛機(jī)氣動(dòng)性能是十分必要的。在以往的研究中，人們多是依賴試驗(yàn)的方法來(lái)模擬進(jìn)排氣條件下飛機(jī)的氣動(dòng)性能。如今，隨著計(jì)算機(jī)技術(shù)的發(fā)展，CFD 技術(shù)已經(jīng)廣泛應(yīng)用于飛機(jī)的設(shè)計(jì)過(guò)程中，這也使得對(duì)推進(jìn)系統(tǒng)和機(jī)體一體化進(jìn)行數(shù)值模擬，建立民機(jī)動(dòng)力影響分析系統(tǒng)成為可能。

所謂發(fā)動(dòng)機(jī)進(jìn)排氣動(dòng)力影響，是指對(duì)于航空發(fā)動(dòng)機(jī)，一般其前部都要配置進(jìn)氣道，而后部配置尾噴管，這樣進(jìn)氣道前面的進(jìn)氣流和尾噴管后面的尾噴流，都會(huì)對(duì)飛行器的外部流動(dòng)產(chǎn)生干擾影響，從而改變飛行器的氣動(dòng)特性。從20世紀(jì)80年代開(kāi)始，國(guó)外就針對(duì)各種發(fā)動(dòng)機(jī)進(jìn)排氣效應(yīng)進(jìn)行了研究，NASA 的Langley研究中心用試驗(yàn)的方法，采用渦輪動(dòng)力模擬器，對(duì)發(fā)動(dòng)機(jī)短艙裝在機(jī)翼下的的布局進(jìn)行了大量的研究，以減少動(dòng)力效應(yīng)帶來(lái)的干擾阻力［1］；在數(shù)值模擬方面，Hirose N［2］、Deese J E［3］等人通過(guò)數(shù)值求解Euler方程，模擬了發(fā)動(dòng)機(jī)的進(jìn)排氣效應(yīng)，得到了進(jìn)排氣效應(yīng)引起唇口激波強(qiáng)度變化的結(jié)論。

本文采用非結(jié)構(gòu)混合網(wǎng)格，通過(guò)數(shù)值求解NS方程，對(duì)發(fā)動(dòng)機(jī)進(jìn)、排氣效應(yīng)進(jìn)行了模擬。首先通過(guò)對(duì)DLR－F4翼身組合體計(jì)算，驗(yàn)證了程序的可靠性。然后采用單獨(dú)TPS風(fēng)洞試驗(yàn)?zāi)Ｐ?，考察了不同進(jìn)排氣條件對(duì)發(fā)動(dòng)機(jī)表面壓力分布的影響，證明本文采用的進(jìn)排氣模擬技術(shù)是可行的。在此基礎(chǔ)上，通過(guò)對(duì)某民機(jī)帶動(dòng)力狀態(tài)進(jìn)行模擬研究，分析了進(jìn)排氣效應(yīng)對(duì)民用大飛機(jī)流場(chǎng)的干擾影響。

1 數(shù)值方法

本文采用CARDC 自主研制的亞跨超聲速流場(chǎng)解算器MFlow 進(jìn)行計(jì)算。MFlow 解算器是基于格心的非結(jié)構(gòu)混合網(wǎng)格和雷諾平均NS方程的大規(guī)模并行流場(chǎng)解算器。它可以使用任意形狀的網(wǎng)格單元，具有較大的靈活性。采用有限體積法對(duì)空間進(jìn)行離散，未知變量位于網(wǎng)格單元的體心。離散方程組的求解采用隱式LU－SGS方法或顯式Runge－Kutta 方法，采用FAS融合多重網(wǎng)格方法加速收斂。MFlow解算器有各種不同的選項(xiàng)可以使用，例如各種空間對(duì)流項(xiàng)和擴(kuò)散項(xiàng)離散格式、各種時(shí)間迭代方法、不同的湍流模型等。

在本文的研究中，采用四面體和三棱柱單元混合的非結(jié)構(gòu)網(wǎng)格，主控方程對(duì)流項(xiàng)采用二階迎風(fēng)Roe通量差分裂格式進(jìn)行離散，采用隱式LU－SGS 時(shí)間項(xiàng)離散方法求解。湍流模型采用SA 一方程湍流模型。

1.1 控制方程

守恒形式的非定?？蓧嚎sNS方程［4］：

其中，Ω表示控制體的體積，?Ω表示控制體封閉面的面積，W為守恒變量，F(xiàn)c為無(wú)粘通量，F(xiàn)v為粘性通量。

1.2 邊界條件

邊界條件的給定及其離散處理方式是數(shù)值求解Euler／NS 方程的重要問(wèn)題之一，不合適的邊界條件會(huì)引起對(duì)真實(shí)系統(tǒng)的不正確模擬，而且對(duì)解的收斂速度和穩(wěn)定性也有很大影響。本文中主要采用的邊界條件有：

（1）遠(yuǎn)場(chǎng)邊界條件：采用基于局部一維Riemann不變量的無(wú)反射邊界條件；

（2）無(wú)滑移物面邊界條件：無(wú)滑移、絕熱條件；

（3）對(duì)稱面條件：無(wú)穿透條件；

（4）發(fā)動(dòng)機(jī)進(jìn)排氣條件：由于發(fā)動(dòng)機(jī)內(nèi)部燃燒和工作過(guò)程相當(dāng)復(fù)雜，我們?cè)跀?shù)值模擬中可以通過(guò)設(shè)定合適的邊界條件，使發(fā)動(dòng)機(jī)進(jìn)排氣效應(yīng)與實(shí)際情況一致，而不去詳細(xì)模擬發(fā)動(dòng)機(jī)內(nèi)部的工況。本文采用的進(jìn)排氣條件如圖1所示。

圖1 發(fā)動(dòng)機(jī)進(jìn)排氣邊界條件示意圖Fig.1 Boundary condition of turbo－fan engine

對(duì)于發(fā)動(dòng)機(jī)進(jìn)氣口，此時(shí)相對(duì)計(jì)算流場(chǎng)來(lái)說(shuō)為出流條件。通過(guò)給定邊界面上的目標(biāo)流量mtarget，在每個(gè)時(shí)間迭代步中，計(jì)算出邊界面的實(shí)際流量mreal，然后根據(jù)mtarget／mreal調(diào)節(jié)邊界面的速度，使其滿足目標(biāo)流量。具體推導(dǎo)如下：

同時(shí)，根據(jù)等熵關(guān)系，有：

其中下標(biāo)c1表示邊界面內(nèi)側(cè)單元的體心值。

所以，有：

邊界面上壓力和密度可表示為：

這種方法的好處是在流場(chǎng)計(jì)算收斂后，發(fā)動(dòng)機(jī)的進(jìn)氣量與實(shí)際情況相一致。

對(duì)于發(fā)動(dòng)機(jī)出流邊界，此時(shí)相對(duì)計(jì)算流場(chǎng)來(lái)說(shuō)為入流條件，我們指定出口邊界面上的總溫、總壓和出口速度方向與x、y、z軸的夾角，壓力p采用外插，其他變量可根據(jù)等熵關(guān)系求得：

邊界面上溫度T可表示為：

當(dāng)?shù)芈曀俸退俣瓤杀硎緸椋?/p>

所以邊界面上的其他變量值可求得：

2 程序驗(yàn)證

2.1 DLR－F4標(biāo)模計(jì)算

為了對(duì)本文采用的數(shù)值計(jì)算方法進(jìn)行考核和驗(yàn)證，我們首先對(duì)DLR－F4標(biāo)模［5］進(jìn)行了計(jì)算，并將試驗(yàn)值和第一屆阻力會(huì)議提供的不同程序計(jì)算結(jié)果進(jìn)行了比較［6］。

計(jì)算條件：M＝0.75；攻角α＝－3°，－2°，－1°，0°，1°，2°；溫度T＝283.15K；雷諾數(shù)Re＝3.0×106（基于平均氣動(dòng)弦長(zhǎng)Cref＝0.1412m）。

圖2 是F4 的計(jì)算網(wǎng)格示意圖。網(wǎng)格單元總數(shù)為2164萬(wàn)個(gè)。其中四面體單元1368萬(wàn)個(gè)，三棱柱單元795萬(wàn)個(gè)。物面單元數(shù)為29.5萬(wàn)個(gè)，物面法向三棱柱網(wǎng)格數(shù)為27個(gè)，物面法向第一層間距約為1.0×10－6m。

圖3給出了本文計(jì)算得到的極曲線和試驗(yàn)結(jié)果以及Tau、NSU3D、USM3Dns等不同軟件計(jì)算得到的結(jié)果比較?？梢钥闯觯疚挠?jì)算結(jié)果落在其它幾個(gè)程序計(jì)算結(jié)果之間，與試驗(yàn)結(jié)果吻合較好。

圖2 DLR－F4網(wǎng)格Fig.2 Grid of DLR－F4

圖3 DLR－F4極曲線Fig.3 Polor curve of DLR－F4

圖4 給出了F4的機(jī)翼典型剖面壓力分布比較?？梢钥闯?，不同程序的計(jì)算結(jié)果都非常接近，與試驗(yàn)結(jié)果符合得也比較好，幾個(gè)程序計(jì)算結(jié)果的前緣吸力峰值都要比試驗(yàn)值低，激波位置靠前，波后壓力系數(shù)偏低，本文的計(jì)算結(jié)果落在其他幾個(gè)程序計(jì)算結(jié)果之間。

圖4 DLR－F4壓力分布Fig.4 Pressure distribution of DLR－F4

2.2 TPS風(fēng)洞試驗(yàn)?zāi)Ｐ陀?jì)算

為了考核本文采取的進(jìn)排氣邊界條件，我們對(duì)日本宇航技術(shù)研究所“NAL－AERO－02－01”TPS（Turbine Powered Simulator）風(fēng)洞試驗(yàn)?zāi)Ｐ停?］進(jìn)行了計(jì)算，并與試驗(yàn)值進(jìn)行比較。計(jì)算模型網(wǎng)格分布如圖5所示。

圖5 TPS風(fēng)洞試驗(yàn)?zāi)Ｐ途W(wǎng)格分布Fig.5 The TPS model and surface grid

為了考察不同流量條件對(duì)結(jié)果的影響，我們選取了兩個(gè)不同流量條件狀態(tài)進(jìn)行了計(jì)算。計(jì)算馬赫數(shù)均為M＝0.8，攻角為0°，具體計(jì)算條件如表1所示。

表1 TPS模型計(jì)算狀態(tài)Table 1 The condition of TPS model

圖6給出了兩種不同工作狀態(tài)下計(jì)算值與試驗(yàn)值表面壓力的結(jié)果對(duì)比，可以看出計(jì)算值和試驗(yàn)值吻合的較好。同時(shí)，從壓力分布的峰值可以看出，隨著發(fā)動(dòng)機(jī)進(jìn)氣流量的增大，發(fā)動(dòng)機(jī)唇口處的激波強(qiáng)度逐漸減弱。

圖7給出了兩種狀態(tài)下子午面馬赫數(shù)分布，可以看出計(jì)算得到的子午面馬赫數(shù)分布合理，在發(fā)動(dòng)機(jī)的出口，由于發(fā)動(dòng)機(jī)出口噴流的影響，形成了較強(qiáng)的剪切，尤其是對(duì)于狀態(tài)1，發(fā)動(dòng)機(jī)出口壓力較高的情形，在局部甚至出現(xiàn)了超聲速。

圖6 兩種狀態(tài)下計(jì)算與試驗(yàn)值表面壓力對(duì)比Fig.6 Comparison of surface pressure coefficients between CFD and experiment

圖7 兩種狀態(tài)下子午面馬赫數(shù)分布Fig.7 The Mach number contours on meridian configuration

通過(guò)以上的計(jì)算對(duì)比分析，可以看出，本文采用的計(jì)算方法可以較好地模擬流場(chǎng)的結(jié)構(gòu)，得到的壓力分布與試驗(yàn)值吻合較好，說(shuō)明本文采用的發(fā)動(dòng)機(jī)進(jìn)排氣模擬技術(shù)是可行的，驗(yàn)證了程序的準(zhǔn)確性和可靠性。

3 某民機(jī)動(dòng)力影響計(jì)算分析

3.1 計(jì)算模型及網(wǎng)格

為了研究發(fā)動(dòng)機(jī)進(jìn)排氣效應(yīng)對(duì)民機(jī)流場(chǎng)的影響，我們選取某典型翼吊式民機(jī)外形，對(duì)有無(wú)動(dòng)力情況進(jìn)行了計(jì)算分析，計(jì)算外形及網(wǎng)格分布如圖8所示。在空間生成四面體單元，在附面層內(nèi)生成三棱柱單元，中間通過(guò)金字塔單元過(guò)渡。半機(jī)網(wǎng)格量為1174萬(wàn)，其中四面體670 萬(wàn)，三棱柱478 萬(wàn)，物面單元為20萬(wàn)，物面法向三棱柱網(wǎng)格數(shù)為30個(gè)，物面法向第一層間距約為1.0×10－5m。

圖8 某民機(jī)外形及網(wǎng)格分布Fig.8 The civil aircraft model and surface gird

計(jì)算狀態(tài)為馬赫數(shù)M＝0.74，攻角α＝8°，發(fā)動(dòng)機(jī)進(jìn)口流量為445.5kg／s，內(nèi)外涵道總壓分別為54.4kPa、55.5kPa，總溫分別為763.3K、287K。

3.2 計(jì)算結(jié)果及討論

圖9給出了有無(wú)動(dòng)力情況下短艙吊架內(nèi)、外兩側(cè)機(jī)翼上下表面的壓力分布?？梢钥闯?，與無(wú)動(dòng)力狀態(tài)相比，帶動(dòng)力情形由于發(fā)動(dòng)機(jī)出口噴流的引射作用，使得機(jī)翼上表面激波位置發(fā)生后移，而對(duì)激波前的壓力分布影響不大。對(duì)于機(jī)翼下表面，由于發(fā)動(dòng)機(jī)出口噴射出來(lái)的氣流壓力較高，導(dǎo)致下表面壓力略有增大，而且越靠近發(fā)動(dòng)機(jī)出口位置影響相對(duì)越明顯。發(fā)動(dòng)機(jī)進(jìn)排氣效應(yīng)對(duì)發(fā)動(dòng)機(jī)外側(cè)機(jī)翼的影響要比內(nèi)側(cè)明顯，且這種影響量隨著離開(kāi)發(fā)動(dòng)機(jī)展向距離的增加呈遞減趨勢(shì)。圖10給出了有無(wú)動(dòng)力情況下機(jī)翼上表面壓力云圖。

圖9 有無(wú)動(dòng)力情況下不同剖面壓力分布對(duì)比Fig.9 Comparison of surface pressure coefficients at different sections between power on and power off

圖10 有無(wú)動(dòng)力情況下機(jī)翼表面壓力云圖Fig.10 Comparison of surface pressure contours on the wing between power on and power off

4 結(jié)論

本文通過(guò)采用非結(jié)構(gòu)混合網(wǎng)格，數(shù)值求解NS方程的方法，對(duì)發(fā)動(dòng)機(jī)進(jìn)排氣效應(yīng)進(jìn)行了數(shù)值模擬分析，得到以下結(jié)論：

（1）本文采用的數(shù)值計(jì)算方法，可以較好地模擬發(fā)動(dòng)機(jī)進(jìn)排氣效應(yīng)下的動(dòng)力影響。

（2）對(duì)于發(fā)動(dòng)機(jī)入口，隨著進(jìn)氣流量的增大，發(fā)動(dòng)機(jī)唇口處的激波強(qiáng)度逐漸減弱。發(fā)動(dòng)機(jī)出口由于噴流的影響，會(huì)形成較強(qiáng)的剪切，局部可能達(dá)到超聲速。

（3）對(duì)于翼吊式民機(jī)外形，在本文的計(jì)算條件下，由于發(fā)動(dòng)機(jī)出口噴流的引射作用，機(jī)翼上表面激波位置發(fā)生后移。發(fā)動(dòng)機(jī)進(jìn)排氣效應(yīng)對(duì)發(fā)動(dòng)機(jī)外側(cè)機(jī)翼的影響比內(nèi)側(cè)明顯。

［1］HENDERSON W P，PATTERSON J C.Propulsion installation characteristics for turbofan transports［R］.AIAA 83－0087，1983.

［2］HIROSE N，ASAI K.3－D Euler flow analysis of fanjet engine and turbine powered simulator with experimental comparison in transonic speed［R］.AIAA 89－1835，1989.

［3］DEESE J E，AGARWAL R K.Calculation of axisymmetric inlet flowfield using the Euler equations［R］.AIAA 83－1853，1983.

［4］朱自強(qiáng).應(yīng)用計(jì)算流體力學(xué)［M］.北京航空航天大學(xué)出版社，1998.

［5］REDEKER G.DLR－F4wing body configuration，in chapter B of A selection of experimental test cases for the validation of CFD codes［R］.AGARD AR－303 21994，1994.

［6］LEVY D W，ZICKUHR T，VASSBERG J.Summary of data from the first AIAA CFD drag prediction workshop［R］.AIAA 2002－0841，2002.

［7］HIROSE N，ASAI K，IKAWA K.Transonic 3－D Euler analysis of flows around fan－jet engine and TPS（turbine powered simulator）［R］.NAL－TR－1045，1989.

［8］LI J，LI F，CHEN H.3－D Flow simulations for generalpowered engine nacelles using Euler equations［R］.AIAA 98－0929，1998.

［9］CHEN H C，YU N L，RUBBERT P E.Flow simulations for general nacelle configurations using Euler equations［R］.AIAA83－0539，1983.

［10］WAI J C，SUN C C，YOSHIHARA H.Transonic perturbation analysis of wing－fuselage－nacelle－pylon configurations with powered jet exhausts［R］.NASA－CR－165852，1982.

［11］BOPPE C W，STERN M A.Simulated transonic flows for aircraft with nacelles，pylons and winglets［R］.AIAA 80－0130，1980.