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

(中國飛行試驗研究院發(fā)動機所， 西安 710089)

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

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

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

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

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

1.1 計算方法

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

1.2 模型及網(wǎng)格

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

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

圖1 簡化后的發(fā)動機模型Fig.1 Engine model after simplify

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

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

BODY為體圖3 整體計算域Fig.3 Intergral computational domain

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

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

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

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

(1)

(2)

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

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

2 渦扇發(fā)動機安裝阻力計算

2.1 風(fēng)洞試驗中安裝阻力計算步驟

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

圖5 飛機機體受力分析Fig.5 Force analysis of plane

整個試驗中側(cè)滑角為0，推力與機體參考軸線的夾角也為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)系中天平測取的升力方向的力；FG0為阻力；DW為風(fēng)軸坐標(biāo)系下的阻力；ΦD-wind為風(fēng)軸坐標(biāo)系中天平測取的推力方向的力；ΔFE-wind-L為與動力裝置狀態(tài)變化相關(guān)的增量力在風(fēng)軸升力方向上的分量；ΔFE-wind-D為與動力裝置狀態(tài)變化相關(guān)的增量力在風(fēng)軸阻力方向上的分量；FG9為發(fā)動機推力。

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

(5)

(6)

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

(7)

(8)

(9)

(10)

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

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

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

(11)

(12)

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

2.2 CFD仿真中安裝阻力計算

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

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

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

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

(13)

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

3 計算結(jié)果

3.1 半模飛機的升力-阻力特性

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

3.2 流場分布

圖7 飛機表面流場分布Fig.7 Flow field distributed on surface of plane

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

3.3 風(fēng)洞試驗值與修正后CFD計算值對比

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

(14)

(15)

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

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

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

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

4 結(jié)論

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

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

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

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