王定奇， 李密， 高翔， 李秋鋒

(中國(guó)飛行試驗(yàn)研究院發(fā)動(dòng)機(jī)所， 西安 710089)

軍用運(yùn)輸機(jī)最重要的用途之一是將一定的載荷(包括人員、設(shè)備等)運(yùn)抵至規(guī)定地點(diǎn)，有效載荷、航程是衡量大型運(yùn)輸機(jī)性能的重要指標(biāo)。因此軍用運(yùn)輸機(jī)在研發(fā)階段以及交付使用前，需對(duì)有效載荷、航程等關(guān)鍵指標(biāo)給出鑒定結(jié)論，對(duì)飛/發(fā)一體化設(shè)計(jì)水平和動(dòng)力裝置性能特性進(jìn)行評(píng)估，而運(yùn)輸機(jī)平臺(tái)性能的優(yōu)劣與飛/發(fā)一體化設(shè)計(jì)水平高低的直觀體現(xiàn)是能否準(zhǔn)確獲得飛機(jī)升阻極曲線。在飛機(jī)研制的設(shè)計(jì)階段通常借助于比例模型風(fēng)洞試驗(yàn)獲得飛機(jī)系統(tǒng)和短艙通流模型的升阻極曲線[1-3]，然而，短艙通流模型風(fēng)洞試驗(yàn)中沒有考慮到發(fā)動(dòng)機(jī)狀態(tài)變化引起的飛機(jī)升力、阻力及俯仰力矩的變化，對(duì)于飛機(jī)升阻極曲線帶來誤差。

當(dāng)短艙為通流模型情況下，風(fēng)洞天平中獲取的風(fēng)軸下的阻力相當(dāng)于機(jī)體阻力；而短艙帶動(dòng)力情況下，風(fēng)洞中天平所獲取的力相當(dāng)于額外推力，此時(shí)機(jī)體系統(tǒng)的阻力需要通過安裝凈推力和風(fēng)洞天平力來確定。在帶動(dòng)力進(jìn)排氣氣動(dòng)特性的影響試驗(yàn)中，國(guó)外主要采用進(jìn)排氣模擬裝置(turbine powered simulator, TPS)和引射噴管裝置兩種形式進(jìn)行了研究[4-5]。TPS和真實(shí)發(fā)動(dòng)機(jī)相比，除不能模擬高溫噴流以及進(jìn)排氣流量比小于1.0外，其余特征與真實(shí)發(fā)動(dòng)機(jī)非常接近，但針對(duì)本次試驗(yàn)條件，TPS尺寸較小，無法在與TPS配套的短艙外表面布置足夠的靜壓測(cè)點(diǎn)?？紤]到本項(xiàng)目的成本和試驗(yàn)周期，選取引射器短艙作為風(fēng)洞動(dòng)力模擬試驗(yàn)的發(fā)動(dòng)機(jī)模擬器。

本文所研究的大涵道比渦扇發(fā)動(dòng)機(jī)，短艙的軸向尺寸較短，進(jìn)排氣流動(dòng)耦合程度較高，試驗(yàn)中無法單獨(dú)將進(jìn)氣溢流阻力與排氣干擾阻力進(jìn)行計(jì)算，因此必須借助全三維仿真(computational fluid dynamics, CFD)方法和縮比模型試驗(yàn)確定各項(xiàng)力的組成。國(guó)外較早開展了針對(duì)帶動(dòng)力飛機(jī)的數(shù)值模擬研究[6-10]，分析了發(fā)動(dòng)機(jī)安裝位置，發(fā)動(dòng)機(jī)狀態(tài)及飛行狀態(tài)等對(duì)安裝性能的影響，并開展了風(fēng)洞試驗(yàn)和數(shù)值仿真相關(guān)性研究，形成了相關(guān)推阻力劃分體系[11-12]。國(guó)內(nèi)這方面起步較晚，譚兆光等[13]就動(dòng)力效應(yīng)做了初步探索，主要為了確定(Navier-Stokes equations, N-S)方程的可行性；但帶動(dòng)力與組合體之間的相互干擾，沒有深入研究。李強(qiáng)等[14]利用CFM56發(fā)動(dòng)機(jī)噴流試驗(yàn)結(jié)果對(duì)比，由于沒有公開的數(shù)據(jù)對(duì)比，存在一定的誤差。Zhang等[15]通過數(shù)值模擬，對(duì)民機(jī)中短艙通流模型和帶動(dòng)力模型的阻力特性進(jìn)行了對(duì)比分析。高翔等[16]利用風(fēng)洞試驗(yàn)數(shù)據(jù)對(duì)數(shù)值計(jì)算方法進(jìn)行驗(yàn)證,總結(jié)了溢流阻力及排氣干擾阻力的變化規(guī)律。

通過上述研究可以看出，目前主要針對(duì)通流模型或帶動(dòng)力模型的流場(chǎng)干擾問題開展了一定研究，而對(duì)于對(duì)帶動(dòng)力運(yùn)輸機(jī)升阻特性CFD與風(fēng)洞試驗(yàn)相關(guān)性研究的較少?，F(xiàn)針對(duì)帶動(dòng)力運(yùn)輸機(jī)CFD仿真，對(duì)不同馬赫數(shù)、不同攻角和不同發(fā)動(dòng)機(jī)狀態(tài)進(jìn)行計(jì)算分析，提取出各項(xiàng)阻力，并通過風(fēng)洞試驗(yàn)進(jìn)行對(duì)比，驗(yàn)證CFD方法獲取飛機(jī)升阻特性正確性。隨后以風(fēng)扇壓比FNPR=1.61時(shí)為基準(zhǔn)，獲取該發(fā)動(dòng)機(jī)狀態(tài)下不同攻角時(shí)，帶動(dòng)力運(yùn)輸機(jī)CFD計(jì)算結(jié)果的修正因子；將該修正因子運(yùn)用到其他發(fā)動(dòng)機(jī)狀態(tài)，使各飛行狀態(tài)及發(fā)動(dòng)機(jī)狀態(tài)下CFD計(jì)算得到升阻特性達(dá)到最優(yōu)解，進(jìn)一步驗(yàn)證帶動(dòng)力運(yùn)輸機(jī)升阻特性的相關(guān)性方法。

1 數(shù)值計(jì)算方法和網(wǎng)格生成

1.1 計(jì)算方法

對(duì)于黏性起主導(dǎo)作用的機(jī)翼擾流問題，其流場(chǎng)伴隨著尾跡混合、流動(dòng)分離機(jī)附面層的干擾等復(fù)雜流動(dòng)特性[17-21]。本文研究采用有限體積法求解該方程，空間離散格式為二階迎風(fēng)Roe格式，時(shí)間推進(jìn)格式為L(zhǎng)U-SGS格式。綜合考慮計(jì)算效率和計(jì)算精度，流場(chǎng)模擬采用加強(qiáng)型S-A湍流模型[22]。

1.2 模型及網(wǎng)格

試驗(yàn)采用的模型為M3飛機(jī)的1∶10動(dòng)力半模模型，參數(shù)如下：參考面積為0.455 m2，平均氣動(dòng)弦長(zhǎng)：0.341 5 m，機(jī)翼展長(zhǎng)：1.444 m，翼展2.5 m，機(jī)翼面積1.2 m2，展弦比6.7。數(shù)值計(jì)算中需要對(duì)幾何模型進(jìn)行相應(yīng)的簡(jiǎn)化：①刪減試驗(yàn)件幾何模型中的引氣管路、測(cè)耙管路等附屬結(jié)構(gòu)，僅保氣動(dòng)影響的幾何型面；②刪減尾噴管進(jìn)口至進(jìn)氣道出口中間的引射器幾何構(gòu)型；③忽略試驗(yàn)件表面縫隙、凸臺(tái)等微小因素的影響。簡(jiǎn)化后的模型如圖1所示。

動(dòng)力系統(tǒng)的內(nèi)外涵采用非結(jié)構(gòu)四面體網(wǎng)格，涵道壁面第1層網(wǎng)格0.4 mm，增長(zhǎng)率1.1，共11層；內(nèi)外涵網(wǎng)格量為300萬，發(fā)動(dòng)機(jī)內(nèi)外涵表面網(wǎng)格分布如圖2所示。

圖1 簡(jiǎn)化后的發(fā)動(dòng)機(jī)模型Fig.1 Engine model after simplify

圖2 發(fā)動(dòng)機(jī)內(nèi)涵表面網(wǎng)格分布Fig.2 Distribution of grid on surface of engine

根據(jù)該飛機(jī)的結(jié)構(gòu)尺寸，半模計(jì)算區(qū)域選擇長(zhǎng)50 m、寬15 m和高15 m的長(zhǎng)方體作為計(jì)算控制域，結(jié)構(gòu)體形狀及飛機(jī)模型在計(jì)算域中的位置如圖3所示。對(duì)整機(jī)采用非結(jié)構(gòu)化網(wǎng)格進(jìn)行劃分，靠近機(jī)翼/短艙/吊掛區(qū)域進(jìn)行網(wǎng)格加密，遠(yuǎn)離機(jī)身區(qū)域計(jì)算網(wǎng)格逐漸稀疏，固壁面均采用7層網(wǎng)格進(jìn)行加密，第1層網(wǎng)格高度為1 mm，增長(zhǎng)率1.2，飛機(jī)半模網(wǎng)格量為550萬，飛機(jī)表面流場(chǎng)分布如圖4所示。

BODY為體圖3 整體計(jì)算域Fig.3 Intergral computational domain

圖4 飛機(jī)表面網(wǎng)格分布Fig.4 Distribution of grid on surface of plane

1.3 數(shù)據(jù)選取

試驗(yàn)中總溫、總壓測(cè)量耙的參數(shù)是以穩(wěn)定狀態(tài)點(diǎn)的直接測(cè)量參數(shù)作為輸入條件的，但試驗(yàn)中的穩(wěn)態(tài)數(shù)據(jù)采集，會(huì)因各種不確定的原因出現(xiàn)一些異常值，如瞬時(shí)或者間斷性的測(cè)量系統(tǒng)故障或參數(shù)真實(shí)的波動(dòng)等，從而無法得到發(fā)動(dòng)機(jī)關(guān)鍵參數(shù)的真實(shí)值，必須剔除所謂的異常值。對(duì)于風(fēng)洞中常規(guī)測(cè)力試驗(yàn)，數(shù)據(jù)采集測(cè)量時(shí)，待速壓穩(wěn)定后，等待2 s，測(cè)試系統(tǒng)采集1 s，每點(diǎn)采集1 000次，所以需對(duì)風(fēng)洞驗(yàn)證試驗(yàn)的原始數(shù)據(jù)進(jìn)行異常值剔除和數(shù)據(jù)平均等預(yù)處理。

根據(jù)阿諾德工程研究中心提出的改進(jìn)的萊茵達(dá)準(zhǔn)則判斷異常值[23]，其基本思想與萊因達(dá)法則相同，以C倍測(cè)量值樣本均值的標(biāo)準(zhǔn)偏差作為置信區(qū)間，超過此區(qū)間的予以剔除。具體表達(dá)式為

(1)

(2)

式(2)中：N為樣本數(shù)；若N≥65，則C等于3，與萊茵達(dá)準(zhǔn)則相同。將測(cè)試系統(tǒng)采集的穩(wěn)態(tài)試驗(yàn)數(shù)據(jù)按照改進(jìn)的萊因達(dá)準(zhǔn)則進(jìn)行剔點(diǎn)和算術(shù)平均，得到對(duì)應(yīng)風(fēng)洞試驗(yàn)點(diǎn)的參數(shù)平均值。以來流Ma=0.2、α=0°的工況下，不同發(fā)動(dòng)機(jī)狀態(tài)下各個(gè)截面參數(shù)，如表1所示。

表1 計(jì)算邊界條件設(shè)置Table 1 Boundary condition of calculation

2 渦扇發(fā)動(dòng)機(jī)安裝阻力計(jì)算

2.1 風(fēng)洞試驗(yàn)中安裝阻力計(jì)算步驟

風(fēng)洞試驗(yàn)采用半模模型，受力如圖5所示。

圖5 飛機(jī)機(jī)體受力分析Fig.5 Force analysis of plane

整個(gè)試驗(yàn)中側(cè)滑角為0，推力與機(jī)體參考軸線的夾角也為0。

LW+ΔFE-wind-L+FG9sinα=ΦL-wind

(3)

FG9cosα-FG0-ΔFE-wind-D-DW=ΦD-wind

(4)

式中：LW為風(fēng)軸坐標(biāo)系中的升力；ΦL-wind為風(fēng)軸坐標(biāo)系中天平測(cè)取的升力方向的力；FG0為阻力；DW為風(fēng)軸坐標(biāo)系下的阻力；ΦD-wind為風(fēng)軸坐標(biāo)系中天平測(cè)取的推力方向的力；ΔFE-wind-L為與動(dòng)力裝置狀態(tài)變化相關(guān)的增量力在風(fēng)軸升力方向上的分量；ΔFE-wind-D為與動(dòng)力裝置狀態(tài)變化相關(guān)的增量力在風(fēng)軸阻力方向上的分量；FG9為發(fā)動(dòng)機(jī)推力。

當(dāng)動(dòng)力裝置工作參考狀態(tài)，在風(fēng)軸坐標(biāo)系下ΔFE-wind-L=0、ΔFE-wind-D=0，可得

(5)

(6)

當(dāng)動(dòng)力裝置偏離工作參考狀態(tài)時(shí)，在風(fēng)軸坐標(biāo)系下，可得

(7)

(8)

(9)

(10)

式中：上標(biāo)0表示工作參考狀態(tài)；上標(biāo)X表示實(shí)際工作狀態(tài)。

圖6 半模飛機(jī)升阻特性曲線Fig.6 Curve of lift and drag for semi-plane

根據(jù)風(fēng)洞試驗(yàn)數(shù)據(jù)，可得

(11)

(12)

式中：Ma為馬赫數(shù)；FNPR為風(fēng)扇壓比；A0/A1為進(jìn)氣道捕獲面積比；f為映射函數(shù)。

2.2 CFD仿真中安裝阻力計(jì)算

外部阻力特性CFD仿真數(shù)據(jù)處理中的氣動(dòng)參數(shù)均可在算例結(jié)果中提取，包括內(nèi)外涵總溫、總壓，環(huán)境溫度、壓力，內(nèi)外涵尾噴管出口流量等，并利用數(shù)值積分獲取半模飛機(jī)所有固壁面的壓差阻力和摩擦阻力。風(fēng)洞中獲取的阻力和升力為附加前體力Φpre與半模飛機(jī)所有固壁面外部阻力Φa/c分別在升力和阻力方向的分量。全三維仿真計(jì)算結(jié)果中提取的升力和阻力應(yīng)與上式保持一致，數(shù)據(jù)處理分為兩部分。

(1)半模飛機(jī)所有固壁面外部阻力Φa/c計(jì)算。體軸坐標(biāo)系下，在CFD計(jì)算結(jié)果中直接提取半模飛機(jī)各個(gè)部分固壁面所受外力在升力和阻力方向的分力，進(jìn)行坐標(biāo)轉(zhuǎn)換，獲得半模飛機(jī)所有固壁面外部阻力在風(fēng)軸坐標(biāo)系下升力和阻力方向的分力。

(2)附加前體力計(jì)算。即CFD獲取的阻力為

D=FN-ΦD-wind=ΔFE+DW=Φpre+Φa/c

(13)

通過上述分析，可知風(fēng)洞試驗(yàn)中，通過六分力天平獲得飛機(jī)在風(fēng)軸下的額外推力；前期校準(zhǔn)箱試驗(yàn)獲得的發(fā)動(dòng)機(jī)噴管特性結(jié)合試驗(yàn)狀態(tài)點(diǎn)可得到發(fā)動(dòng)機(jī)的標(biāo)準(zhǔn)凈推力。當(dāng)選定發(fā)動(dòng)機(jī)某一狀態(tài)為參考狀態(tài)時(shí)，可以通過標(biāo)準(zhǔn)凈推力、附加前提力及天平力，計(jì)算出風(fēng)軸下的機(jī)體阻力。

3 計(jì)算結(jié)果

3.1 半模飛機(jī)的升力-阻力特性

分別對(duì)比了9個(gè)不同攻角、5個(gè)發(fā)動(dòng)機(jī)狀態(tài)和3個(gè)馬赫數(shù)條件下的飛機(jī)升阻特性曲線(圖6)，試驗(yàn)結(jié)果表明：①升力系數(shù)隨著攻角的增大而變大，當(dāng)攻角大于12°后，升力系數(shù)增幅逐漸減小，最大為1.2；阻力系數(shù)隨著攻角的增大而變大，當(dāng)攻角大于12°時(shí)，阻力系數(shù)迅速增大；主要是由于攻角大于12°后，氣流通過機(jī)翼后分離較大，出現(xiàn)了一定的失速；②相同馬赫數(shù)下，不同發(fā)動(dòng)機(jī)工作狀態(tài)下，升力系數(shù)隨著攻角的變化曲線基本重合，而阻力系數(shù)隨著發(fā)動(dòng)機(jī)工作狀態(tài)的增大而增大，表明發(fā)動(dòng)機(jī)工作狀態(tài)對(duì)升力系數(shù)影響較小，而阻力曲線影響較大。

3.2 流場(chǎng)分布

圖7 飛機(jī)表面流場(chǎng)分布Fig.7 Flow field distributed on surface of plane

通過飛機(jī)表面壓力云圖(圖7)可以看出，相同攻角和發(fā)動(dòng)機(jī)狀態(tài)下，隨著來流馬赫數(shù)的增大，在短艙前緣、機(jī)翼前緣壓力升高，來流馬赫數(shù)越大機(jī)翼上表面壓力越低，飛機(jī)的升力隨之增大。相同馬赫數(shù)、相同發(fā)動(dòng)機(jī)狀態(tài)下，隨著飛行攻角的增大，機(jī)翼下表面壓力隨之增大，且起飛構(gòu)型下副翼和襟翼的壓力增大程度明顯，通過增大Ma和攻角α可以快速增大飛機(jī)的升力。

3.3 風(fēng)洞試驗(yàn)值與修正后CFD計(jì)算值對(duì)比

由于本次試驗(yàn)?zāi)Ｐ蜁r(shí)不僅存在飛機(jī)外流場(chǎng)的流動(dòng)，還存在發(fā)動(dòng)機(jī)噴管的內(nèi)流與外流還有耦合情況，且整體計(jì)算模型尺寸較小，當(dāng)發(fā)動(dòng)機(jī)狀態(tài)變化時(shí)，升力、阻力系數(shù)變化的絕對(duì)量有限，同時(shí)CFD仿真過程不可避免的與風(fēng)洞中實(shí)際流動(dòng)情況存在著差異，因此需對(duì)CFD計(jì)算結(jié)果進(jìn)行修正。以FNPR=1.61時(shí)，試驗(yàn)獲得的升力系數(shù)(CL)與阻力系數(shù)(CD)作為基準(zhǔn)，在不同攻角下，獲得CFD計(jì)算結(jié)果的修正因子，即

(14)

(15)

式中：下標(biāo)FNPR表示風(fēng)扇壓比；下標(biāo)EXP表示試驗(yàn)數(shù)據(jù)；下標(biāo)CFD表示仿真計(jì)算數(shù)據(jù)。

根據(jù)以上修正因子，選取發(fā)動(dòng)機(jī)狀態(tài)點(diǎn)(FNPR=1.53)的CFD計(jì)算結(jié)果進(jìn)行相關(guān)性驗(yàn)證，結(jié)果如圖8所示。

通過對(duì)比可以看出：①CFD仿真與風(fēng)洞試驗(yàn)獲得的升力系數(shù)趨勢(shì)一致，Ma=0.15與Ma=0.2情況下，升力系數(shù)最大誤差僅為1.8%；Ma=0.1時(shí)，在大攻角條件下(α=15°)誤差較大，達(dá)到7.1%，其余攻角條件下最大誤差為2.5%；②CFD仿真與風(fēng)洞試驗(yàn)獲得的阻力系數(shù)趨勢(shì)一致，Ma=0.15與Ma=0.2情況下，阻力系數(shù)與試驗(yàn)值基本貼合；Ma=0.1時(shí)，在小攻角(α<5°)曲線誤差較大。

修正參考狀態(tài)為FNPR=1.6圖8 FNPR=1.53狀態(tài)下升阻特性對(duì)比Fig.8 Comparison of lift and drag characteristic at FNPR=1.53

4 結(jié)論

以M3半模飛機(jī)配裝某大涵道比渦扇發(fā)動(dòng)機(jī)三維模型為研究對(duì)象，基于分區(qū)拼接網(wǎng)格對(duì)發(fā)動(dòng)機(jī)內(nèi)流和飛機(jī)外流場(chǎng)網(wǎng)格進(jìn)行拼接，并設(shè)置與風(fēng)洞試驗(yàn)相同的工況點(diǎn)進(jìn)行流場(chǎng)數(shù)值仿真。分別對(duì)比了9個(gè)不同攻角、5種發(fā)動(dòng)機(jī)狀態(tài)和3個(gè)馬赫數(shù)條件下的飛機(jī)升阻特性曲線進(jìn)行對(duì)比，得到如下結(jié)論。

(1)升力系數(shù)隨攻角增大而增大，當(dāng)攻角大于12°后由于氣流分離嚴(yán)重升力系數(shù)逐漸降低，最大升力系數(shù)為1.2；阻力系數(shù)隨攻角增大而增大，攻角大于12°后，阻力系數(shù)迅速增大；

(2)相同Ma，隨著發(fā)動(dòng)機(jī)狀態(tài)的增大飛機(jī)升力系數(shù)變化不明顯，表明翼吊形式的發(fā)動(dòng)機(jī)排氣系統(tǒng)對(duì)飛機(jī)升力影響?。浑S著發(fā)動(dòng)機(jī)狀態(tài)增大阻力系數(shù)增大，表明發(fā)動(dòng)機(jī)工作狀態(tài)對(duì)阻力曲線影響較大；

(3)通過引入升阻修正系數(shù)，將CFD與風(fēng)洞試驗(yàn)結(jié)果進(jìn)行對(duì)比，升力系數(shù)相對(duì)比整體誤差在4%以內(nèi)；而阻力系數(shù)量級(jí)較小，與試驗(yàn)計(jì)算的趨勢(shì)基本一致，表明本文的CFD計(jì)算結(jié)果可滿足后續(xù)分析的基本要求。