地面效應(yīng)對反推狀態(tài)下發(fā)動機(jī)進(jìn)口流場影響數(shù)值研究

周凱，王志強(qiáng)

(南京航空航天大學(xué) 能源與動力學(xué)院，江蘇 南京 210016)

0 引言

隨著渦扇發(fā)動機(jī)涵道比的增加，風(fēng)扇流量大幅度增加[2]，引起反推力裝置流通面積增加，增大反推力裝置軸向尺寸，縮短了反推力裝置與地面距離，加劇了反推氣流與地面相互作用。

地面應(yīng)該模擬成一個以自由流速度移動的地板，當(dāng)?shù)匕逡苿铀俣冗_(dá)到自由流速度時可消除固定地面產(chǎn)生的邊界層影響。2005年G. H. Hegen利用帶有移動帶的DNW-LLF風(fēng)洞對TPS模型開展不同滑跑速度下地面效應(yīng)對反推氣流重吸入的影響。有無移動帶地板風(fēng)洞試驗示意圖如圖1所示。國內(nèi)反推力裝置性能試驗大都在帶有固定地板的風(fēng)洞測試，對地面效應(yīng)影響的研究較少[3-5]。

圖1 有無移動帶地板風(fēng)洞試驗示意圖

本文針對某大涵道比渦扇發(fā)動機(jī)的簡化模型，開展著陸滑跑狀態(tài)下反推氣流流場的數(shù)值研究，分析不同滑跑速度下地面效應(yīng)對反推狀態(tài)下大涵道比渦扇發(fā)動機(jī)進(jìn)口流場畸變影響。

1 建模與網(wǎng)格劃分

本文研究對象為大涵道比渦扇發(fā)動機(jī)1∶1全尺寸模型。主要由發(fā)動機(jī)短艙、吊架、進(jìn)氣錐、內(nèi)涵噴管和葉柵式反推裝置組成，如圖2所示。

圖2 計算對象示意圖

在發(fā)動機(jī)外建立了半圓柱外流區(qū)域，如圖3(a)所示。發(fā)動機(jī)中軸線距離地面為H，發(fā)動機(jī)直徑D(此時H/D約為2)，發(fā)動機(jī)軸向長度為L。計算域半徑為14D，計算域長度取為發(fā)動機(jī)軸向長度L的10倍，計算域進(jìn)口距離發(fā)動機(jī)進(jìn)口約為5.5L。

采用ICEM對計算域進(jìn)行網(wǎng)格劃分，外流場流域最大網(wǎng)格尺度為1024mm，在流動復(fù)雜區(qū)域發(fā)動機(jī)周圍和發(fā)動機(jī)進(jìn)氣道附近設(shè)置雙加密區(qū)，如圖3(b)所示，發(fā)動機(jī)進(jìn)氣唇口區(qū)域最大網(wǎng)格尺度為16mm，并在壁面設(shè)置了10層邊界層網(wǎng)格，第一層網(wǎng)格為0.2mm，延展比為1.2。反推力裝置域相對尺寸較小，最大網(wǎng)格尺寸設(shè)為32mm(圖3(c))，其中葉柵區(qū)面網(wǎng)格最大尺寸為2mm(圖3(d))，并在外流域與反推力裝置內(nèi)域交界處反推氣流出口設(shè)置相同網(wǎng)格參數(shù)，最后生成網(wǎng)格單元約為2100萬(elements)，記為Grid_2100萬。通過合理改變上述網(wǎng)格參數(shù)，共生成4套網(wǎng)格分別為Grid_1000萬、Grid_1600萬、Grid_2100萬、Grid_5500萬。

圖3 計算網(wǎng)格示意圖

2 數(shù)值計算方法

2.1 計算方法

本文數(shù)值模擬采用ANASYS CFX軟件。為了驗證本文數(shù)值模擬方法的可靠性和所采用網(wǎng)格密度、湍流模型的合理性，首先對臺架試驗的5個不同狀態(tài)點(diǎn)開展反推氣流流場的數(shù)值模擬研究，通過與試驗數(shù)據(jù)對比，驗證本文采用的數(shù)值模擬方法的準(zhǔn)確性。

計算時的邊界條件給定如圖4所示：外流域主要由遠(yuǎn)場、遠(yuǎn)場進(jìn)口、遠(yuǎn)場出口以及地面邊界組成，在驗證數(shù)值模擬方法準(zhǔn)確性計算中三遠(yuǎn)場為開放邊界，給定試驗記錄的環(huán)境大氣總溫、總壓；外流域的下邊界為無滑移壁面，用于模擬試驗中固定地面；發(fā)動機(jī)進(jìn)口截面設(shè)定為質(zhì)量流量速率出口；發(fā)動機(jī)內(nèi)涵噴管進(jìn)口和反推力裝置進(jìn)口為計算域的進(jìn)口邊界，給定相應(yīng)試驗狀態(tài)下的總溫、總壓。

師：拆開容易，拼裝難呀！同學(xué)們能把拆開的魔方裝回去嗎？（生紛紛響應(yīng)）現(xiàn)在就來比一比誰拼得快，完成的同學(xué)老師給他戴上一個皇冠。

圖4 邊界條件示意圖

2.2 網(wǎng)格密度和湍流模型選定

國內(nèi)外已開展了大量的大涵道比渦扇發(fā)動機(jī)反推氣流流場數(shù)值模擬研究。DE Andrade F O等人[6-10]分別選用k-ω、SST、k-ε湍流模型對大涵道比渦扇發(fā)動機(jī)反推擾流流場展開了數(shù)值模擬研究。

從上述研究可以看出，反推流場數(shù)值模擬方法仍存在以下問題：第一，湍流模型的選定仍是關(guān)鍵；第二，研究對象幾何模型復(fù)雜，部件尺寸差異較大，網(wǎng)格參數(shù)是影響數(shù)值模擬精確性的重要因素。

本文分別采用4套網(wǎng)格對臺架試驗的5個狀態(tài)點(diǎn)開展反推氣流流場的數(shù)值模擬研究。

圖5為不同網(wǎng)格量計算得到的不同狀態(tài)點(diǎn)下推力與試驗值的相對偏差。從上述4套網(wǎng)格計算結(jié)果可以看出，各套網(wǎng)格計算得到的推力值與試驗結(jié)果的偏差趨勢基本一致，總體網(wǎng)格量1000萬計算結(jié)果較其他3套網(wǎng)格推力相對誤差更大，基本都超過5%，網(wǎng)格增加至2100萬后，數(shù)值模擬結(jié)果與試驗推力貼合最好，推力相對誤差集中在±5%左右，故采用Grid_2100萬為本次計算網(wǎng)格。其中，推力相對誤差ε如下式：

式中：TCFD表示數(shù)值模擬計算推力；TTest表示推力試驗值。

圖5 不同網(wǎng)格密度反推力相對誤差圖

本文在選用Grid_2100萬網(wǎng)格的前提下選取了SST模型、標(biāo)準(zhǔn)k-ω模型以及標(biāo)準(zhǔn)k-ε模型，分別開展了反推狀態(tài)下的排氣系統(tǒng)流場計算。

從圖6可以看出，各湍流模型計算得到的發(fā)動機(jī)推力與試驗值的偏差隨工況的變化趨勢是一致的。從圖可以看出，標(biāo)準(zhǔn)k-ω模型計算得到的發(fā)動機(jī)推力最大，而標(biāo)準(zhǔn)k-ε模型的計算結(jié)果與試驗結(jié)果最接近。

圖6 不同湍流模型反推力相對誤差圖

3 地面效應(yīng)數(shù)值模擬結(jié)果分析

本文采用Grid_2100萬選擇標(biāo)準(zhǔn)k-ε模型，分別選擇0.03、0.04、0.06、0.08、0.10、0.15、0.22、0.25等來流馬赫數(shù)[11]為多個穩(wěn)態(tài)瞬間進(jìn)行定常計算。計算中改遠(yuǎn)場進(jìn)口為速度進(jìn)口，給定相應(yīng)的來流馬赫數(shù)，其余邊界設(shè)定保持不變。通過改變地面為固定無滑移壁面和移動無滑移壁面(速度與遠(yuǎn)場進(jìn)口相等)，分析地面效應(yīng)對反推狀態(tài)下大涵道比渦扇發(fā)動機(jī)進(jìn)口流場影響。選用ARP1420[12]所規(guī)定的周向穩(wěn)態(tài)總壓畸變指數(shù)Δσ，計算發(fā)動機(jī)進(jìn)口AIP(aerodynamic interface plane)截面上的周向穩(wěn)態(tài)總壓畸變指數(shù)，來評估由地面效應(yīng)引起的發(fā)動機(jī)進(jìn)口流場的總壓不均勻性。

為了獲取AIP截面上的畸變指數(shù)，在AIP截面上布置若干個數(shù)值測點(diǎn)，周向均勻分布的16個點(diǎn)，沿徑向布置10個測點(diǎn)，這10個測點(diǎn)是AIP截面上10個等環(huán)面的面積中心點(diǎn)，再加上截面中心的1個測點(diǎn)，AIP截面上共布置有161個測點(diǎn)，如圖7所示。

圖7 AIP截面測點(diǎn)分布

3.1 固定地板和移動地板計算結(jié)果分析

圖8為地面固定時不同來流馬赫數(shù)的反推流場。圖中給出了發(fā)動機(jī)進(jìn)口和反推出口流線和地面靜壓云圖，圖中紅色線為流經(jīng)地面高靜壓區(qū)域流線(如圖中紅色標(biāo)記區(qū)域)(因本刊為黑白印刷，有疑問之處可咨詢作者)。隨著來流馬赫數(shù)的增加，發(fā)動機(jī)進(jìn)口捕捉流管直徑逐漸減小，反推氣流區(qū)域也相應(yīng)縮減，慢慢形成較為規(guī)整的反推氣流流管。在來流Ma≥0.10時，發(fā)動機(jī)進(jìn)口捕捉流管與反推流管在反推力裝置前后形成分離的兩部分，基本不再相互影響。但是在Ma數(shù)較小時，反推氣流流管影響區(qū)域直徑較大，來自反推裝置底部葉柵的反推射流在發(fā)動機(jī)進(jìn)氣唇口前部區(qū)域接觸地面形成速度滯止區(qū)，即高靜壓區(qū)域。在高靜壓區(qū)附近，近地面來流與反推氣流相遇，在反推氣流和發(fā)動機(jī)進(jìn)口捕捉流管的作用下，形成向上氣流，來自地面的上升氣流(如紅色流線示)對捕捉流管底部產(chǎn)生影響。在滑跑速度較小時，如Ma=0.03近地面氣流被吸入發(fā)動機(jī)進(jìn)口，隨著自由來流馬赫數(shù)的增加，反推氣流流管與發(fā)動機(jī)進(jìn)口捕捉流管直徑減小，反推氣流影響區(qū)域后移，近地面流線上升趨勢減弱，最終趨于一條平直流線帶。

圖8 不同Ma數(shù)下反推氣流擾流流線分布

為了進(jìn)一步說明地面效應(yīng)對于發(fā)動機(jī)進(jìn)口流場的影響，圖9分別給出Ma數(shù)分別為0.08、0.22時發(fā)動機(jī)進(jìn)口截面總壓恢復(fù)系數(shù)σ分布云圖。由圖可知，在來流Ma較小時通過給定地面移動速度對由固定地面邊界層發(fā)展引起的位擾動有抑制作用， AIP截面底部總壓損失區(qū)縮??；在Ma較大時，兩種邊界條件計算結(jié)果基本一致。由于單發(fā)反推流場計算中無機(jī)身、機(jī)翼等干擾，反推氣流無明顯直接再吸入現(xiàn)象，由固定地面附面層低能氣流被卷吸進(jìn)入發(fā)動機(jī)引發(fā)的進(jìn)氣畸變，程度小，區(qū)域僅在AIP截面底部較為明顯。

圖9 不同馬赫數(shù)發(fā)動機(jī)進(jìn)口截面總壓恢復(fù)系數(shù)σ分布云圖

圖10分別給出了兩種地面邊界條件下AIP截面的周向穩(wěn)態(tài)總壓畸變指數(shù)隨來流Ma數(shù)變化曲線，由圖可知，周向穩(wěn)態(tài)總壓畸變指數(shù)隨相對來流Ma數(shù)減小而慢慢增大，大致呈階梯狀，與早年國外學(xué)者JOHN H P[13]與HEGEN G H等人[14]通過進(jìn)口截面溫度特征參數(shù)表示反推氣流對發(fā)動機(jī)進(jìn)口的擾動影響結(jié)果相似。從圖中可以看出，兩種壁面邊界條件下AIP截面周向穩(wěn)態(tài)總壓畸變指數(shù)隨來流Ma數(shù)變化趨勢雖然基本一致，但圖中線2整體略低于線1，說明移動地板可以達(dá)到消除固定地板邊界層發(fā)展進(jìn)氣的發(fā)動機(jī)進(jìn)口底部不均勻性。

圖10 周向穩(wěn)態(tài)總壓畸變指數(shù)隨相對來流Ma數(shù)變化曲線

3.2 離地高度影響分析

在前文研究的基礎(chǔ)上，選擇相對來流Ma=0.08，地面為固定無滑移邊界，通過改變H/D大小，開展離地高度對反推狀態(tài)下發(fā)動機(jī)進(jìn)口流場影響研究。

圖11給出了發(fā)動機(jī)進(jìn)口AIP截面周向穩(wěn)態(tài)總壓畸變指數(shù)隨H/D變化曲線。從圖11可以看出，畸變指數(shù)隨離地高度的降低而減小，到H/D=1.5后，畸變指數(shù)又有所增大，最后在0.004 5附近浮動。

圖11 周向穩(wěn)態(tài)總壓畸變指數(shù)隨H/D變化曲線

圖12分別給出了H/D=1.2、1.8反推流線和發(fā)動機(jī)進(jìn)口流線，地面為靜壓云圖。對比兩圖可以看出：離地高度的下降使反推氣流更早接觸地面，高靜壓區(qū)向發(fā)動機(jī)后部移動，反推流管影響區(qū)域減小，降低了對相對來流的卷吸作用，故H/D在1.5～20之間，畸變指數(shù)隨離地高度的降低而減?。欢S著離地高度進(jìn)一步減小，發(fā)動機(jī)進(jìn)口捕捉流管慢慢接觸地面，吸入邊界層氣體，可能造成地面渦吸入，引起發(fā)動機(jī)進(jìn)口流場畸變，周向穩(wěn)態(tài)總壓畸變增大。

圖12 H/D=1.2、1.8反推氣流擾流流線分布

4 結(jié)語

本文針對反推狀態(tài)下大涵道比渦扇發(fā)動機(jī)開展反推擾流流場數(shù)值模擬研究，分析地面效應(yīng)對反推狀態(tài)下大涵道比渦扇發(fā)動機(jī)進(jìn)口流場影響，并通過穩(wěn)態(tài)周向總壓畸變隨滑跑速度的變化直觀展現(xiàn)，為以后相關(guān)風(fēng)洞試驗以及大涵道比渦扇發(fā)動機(jī)設(shè)計[15]提供一定的有用信息。