張宇飛, 陳海昕, 方曉明

(清華大學(xué) 航天航空學(xué)院, 北京 100084)

0 引 言

現(xiàn)代飛機(jī)的設(shè)計(jì)和制造日趨精細(xì)化，研制過(guò)程強(qiáng)烈依賴于多家制造商的分工協(xié)作。飛機(jī)機(jī)體與發(fā)動(dòng)機(jī)這兩大部件一般都由不同公司完成研制。在真實(shí)飛行中，飛機(jī)和發(fā)動(dòng)機(jī)安裝在一起，一方面受到各自的氣動(dòng)阻力，另一方面會(huì)產(chǎn)生新的干擾阻力。但是飛行試驗(yàn)僅能測(cè)得推力與阻力的合力，必須采用輔助手

段標(biāo)定推力和阻力。在平衡飛行狀態(tài)下，飛機(jī)的推力與阻力大小相等，方向相反，所以發(fā)動(dòng)機(jī)推力大小的確定會(huì)影響真實(shí)飛行升阻比的判斷，從而影響著民用客機(jī)的市場(chǎng)競(jìng)爭(zhēng)力。在極端情況下，飛行包線邊界處還可能出現(xiàn)推力不能克服阻力的情況，此時(shí)就更需要定量分析推力與阻力。飛行推力與阻力的確定長(zhǎng)期以來(lái)一直受到國(guó)外飛機(jī)和發(fā)動(dòng)機(jī)公司的重視，美國(guó)Society of Automotive Engineers(SAE)針對(duì)渦扇發(fā)動(dòng)機(jī)、螺旋槳、涵道風(fēng)扇、推力矢量發(fā)動(dòng)機(jī)等，均出臺(tái)了一系列標(biāo)準(zhǔn)可供參考，比如AIR1419、AIR1678、AIR1703、AIR4065、AIR5450、AIR5771、AIR6007等。并且這些標(biāo)準(zhǔn)一直隨著技術(shù)發(fā)展在持續(xù)更新，以AIR1703[2]為例，在1985年出了第一版之后，又在2012年重新制定了修訂版。

推阻分解方法(Thrust/Drag Bookkeeping)是SAE AIR1703標(biāo)準(zhǔn)中確定飛行中推力和阻力的一種通用方法[1-2]，該方法能夠定量計(jì)算作用在飛行系統(tǒng)內(nèi)的氣動(dòng)阻力和發(fā)動(dòng)機(jī)推力，提出一種飛機(jī)和發(fā)動(dòng)機(jī)制造商都能夠認(rèn)可的推力/阻力分解方案，并獲得一套不依賴于發(fā)動(dòng)機(jī)油門狀態(tài)的飛機(jī)配平升力、阻力曲線[2]。此方法已經(jīng)得到波音、通用電氣、普拉特·惠特尼、羅爾斯·羅伊斯等國(guó)外飛機(jī)和發(fā)動(dòng)機(jī)公司的認(rèn)可，并在波音777飛機(jī)研制中應(yīng)用[2]。然而在國(guó)內(nèi)，由于民用飛機(jī)的發(fā)展相對(duì)落后，推阻分解方法的研究和應(yīng)用還相對(duì)較少，尚缺乏一種行業(yè)通用的飛機(jī)和發(fā)動(dòng)機(jī)公司均能接受的推力和阻力計(jì)算方法，同時(shí)也缺乏可靠的標(biāo)準(zhǔn)算例數(shù)據(jù)供試驗(yàn)和數(shù)值計(jì)算作為驗(yàn)證參考。

推阻分解方法的主要原理是利用流體力學(xué)的控制體理論，針對(duì)發(fā)動(dòng)機(jī)短艙的內(nèi)外流道分別進(jìn)行控制體分析，辨別推力與阻力產(chǎn)生的區(qū)域[3-4]。在一般的推力/阻力劃分體系下，發(fā)動(dòng)機(jī)短艙外罩的阻力算做飛機(jī)阻力；進(jìn)氣道和內(nèi)外涵噴管的阻力算做發(fā)動(dòng)機(jī)的推力損失[5]，其損失大小決定于發(fā)動(dòng)機(jī)的設(shè)計(jì)水平；機(jī)體與發(fā)動(dòng)機(jī)安裝后的干擾造成的額外推力損失(如噴流與掛架的干擾、外流速度對(duì)噴流的抑制作用等)稱為干擾阻力，計(jì)入飛機(jī)阻力范圍[2]。Dusa等人[6]給出了各項(xiàng)損失的大小范圍，其中風(fēng)扇噴管的推力損失約占總推力的0.5%～0.8%，核心機(jī)噴管的損失約0.2%～0.3%，噴流作用在掛架上的沖刷阻力約0.5%～1.0%，氣流膨脹、噴流/掛架/機(jī)翼之間的干擾等損失約0.3%～0.9%。進(jìn)氣道阻力和噴管/后體的阻力均與發(fā)動(dòng)機(jī)狀態(tài)相關(guān)，被稱為油門依賴力，進(jìn)氣道溢流阻力和噴管/后體的噴流效應(yīng)阻力均屬于此類。為了確定飛行狀態(tài)的氣動(dòng)力，一般需要選取一個(gè)典型狀態(tài)作為參考狀態(tài)，獲得該狀態(tài)的詳細(xì)氣動(dòng)力參數(shù)，其他飛行狀態(tài)即在此狀態(tài)下附加油門依賴的氣動(dòng)力。為了保持飛機(jī)設(shè)計(jì)過(guò)程的一致性，推阻分解需要遵循的準(zhǔn)則有[1]：

1) 所有部件受力不能遺漏也不能重復(fù)；

2) 分解方法應(yīng)盡量反映部件性能；

3) 選定的參考狀態(tài)應(yīng)接近飛行狀態(tài)，且具有可行的修正方法；

4) 從初步設(shè)計(jì)直到飛行試驗(yàn)中推阻分解方法應(yīng)保持統(tǒng)一，以便追溯研制過(guò)程的性能變化。

推阻分解的核心是準(zhǔn)確計(jì)算發(fā)動(dòng)機(jī)的推力損失。在發(fā)動(dòng)機(jī)設(shè)計(jì)過(guò)程中，發(fā)動(dòng)機(jī)推力及損失可以由地面試車臺(tái)或者高空模擬試車臺(tái)確定[2]。兩種方式下發(fā)動(dòng)機(jī)進(jìn)氣道均需要采用鐘形罩進(jìn)行引流，并保持短艙外罩流速接近于0。推力測(cè)量裝置與鐘形罩不連接，以排除進(jìn)氣阻力的影響。采用這樣的測(cè)試設(shè)備，可以盡量排除短艙外罩阻力對(duì)推力的干擾，獲得發(fā)動(dòng)機(jī)的凈推力，并確定進(jìn)氣道和噴管的推力損失，最終獲得噴管的兩個(gè)重要的性能系數(shù)：流量系數(shù)和速度系數(shù)(也稱作推力系數(shù))。在飛行狀態(tài)下，可以利用地面測(cè)試所獲得的流量系數(shù)和速度系數(shù)，以及飛行溫度、壓力、速度等參數(shù)，計(jì)算飛行狀態(tài)的發(fā)動(dòng)機(jī)推力。為了獲得更準(zhǔn)確的推力數(shù)據(jù)，地面測(cè)試獲得的性能系數(shù)需要經(jīng)過(guò)一系列的修正才能用于飛行推力的確定[2]。但是，修正也僅考慮了飛行的外流速度對(duì)于噴管流量的抑制效應(yīng)，而忽略了外流速度對(duì)于速度系數(shù)的影響[7]。在通行的推阻分解框架下，外流速度對(duì)于速度系數(shù)的影響所帶來(lái)的損失歸于飛機(jī)阻力部分，并得到了飛機(jī)公司的認(rèn)可[2]。

隨著計(jì)算流體力學(xué)(CFD)技術(shù)和計(jì)算機(jī)硬件技術(shù)的發(fā)展，飛行器設(shè)計(jì)越來(lái)越依賴數(shù)值仿真。采用CFD方法進(jìn)行推阻分解，可以很大程度上避免試驗(yàn)測(cè)試的設(shè)備限制、尺度效應(yīng)以及誤差修正問(wèn)題，能夠盡可能接近飛行狀態(tài)，為機(jī)體/動(dòng)力集成帶來(lái)極大便利。近年來(lái)，基于CFD的推阻分解計(jì)算精度研究在國(guó)外已經(jīng)逐漸成為了一個(gè)研究熱點(diǎn)，美國(guó)航空航天學(xué)會(huì)(AIAA)自2011年至今已經(jīng)開展了三屆推進(jìn)空氣動(dòng)力研討會(huì)(Propulsion Aerodynamics Workshop，PAW)，其主要目的就是探討CFD方法對(duì)于噴管推力損失的計(jì)算精度，著重關(guān)注流量系數(shù)和速度系數(shù)的計(jì)算誤差。第一屆PAW的算例是一個(gè)簡(jiǎn)單的軸對(duì)稱收縮噴管[8-10]，第二、三屆算例是民機(jī)雙涵道發(fā)動(dòng)機(jī)噴管[11-12]。

然而在國(guó)內(nèi)，對(duì)這方面研究的關(guān)注還相對(duì)較少，相關(guān)文獻(xiàn)并不多見。本論文主要開展基于CFD方法的推阻分解研究，首先介紹推阻分解的控制體分析方法及原理，然后采用三個(gè)典型算例進(jìn)行推阻分解方法進(jìn)行應(yīng)用驗(yàn)證，分析和驗(yàn)證CFD方法對(duì)于推阻分解的計(jì)算效果，為飛機(jī)設(shè)計(jì)提供參考。

1 發(fā)動(dòng)機(jī)短艙控制體受力分析

針對(duì)典型的民用航空發(fā)動(dòng)機(jī)，推阻分解的控制體分析時(shí)一般不考慮核心機(jī)的旋轉(zhuǎn)部件。現(xiàn)代民機(jī)為了提高推進(jìn)效率，一般采用高涵道比發(fā)動(dòng)機(jī)，噴管分為內(nèi)涵和外涵兩個(gè)流道，圖1顯示了一個(gè)典型的民用發(fā)動(dòng)機(jī)的控制圖示意圖。以唇口駐點(diǎn)為分界線，進(jìn)入發(fā)動(dòng)機(jī)內(nèi)部的流體產(chǎn)生推力，從外部流過(guò)的流體產(chǎn)生阻力。如圖1所示，控制體可以分為三個(gè)區(qū)域：進(jìn)氣流管(從E0到E1截面)；發(fā)動(dòng)機(jī)機(jī)體區(qū)域(從E1到E8和E18截面)；噴流流管(E8到E9截面，E18到E19截面)。

在發(fā)動(dòng)機(jī)具有外流速度的狀態(tài)下，根據(jù)受力關(guān)系，可知發(fā)動(dòng)機(jī)的凈推力FNPF為：

(1)

根據(jù)控制體表面的受力平衡關(guān)系，尾跡區(qū)從E8/E18到E9/E19截面的受力平衡可表示為式(2)，入口從E0到E1截面的受力平衡可以表示為式(3)。其中u9和u19分別表示E8和E18截面流體自由膨脹到環(huán)境大氣后的速度，u0表示自由來(lái)流氣體的速度；ΦPost和ΦPre分別為尾跡和入口流管的側(cè)面受力。

(2)

(3)

整理上述式(1)～式(3)，可以得到發(fā)動(dòng)機(jī)凈推力如式(4)，這也是選取E0到E19截面間的區(qū)域作為一個(gè)控制體的受力形式。

(ΦPre+ΦExt+ΦPost)

=FN-D

(4)

在試車臺(tái)試驗(yàn)中，發(fā)動(dòng)機(jī)性能標(biāo)定過(guò)程均盡量保證短艙外部的流速接近0，并且采用鐘形罩將進(jìn)氣流管的作用力分離，因此可近似認(rèn)為ΦPre=ΦExt=0，定義此時(shí)的凈推力為Fc。

=Fc

(5)

(6)

(7)

由于外流速度的不同，試車臺(tái)狀態(tài)和飛行狀態(tài)的流動(dòng)具有較大的差異。但是如引言所述，在通行的推阻分解方法中隱含了一個(gè)約定，即僅計(jì)入外流速度對(duì)于噴管流量的抑制效應(yīng)，而忽略了外流對(duì)于速度系數(shù)的影響，這種方式的誤差所帶來(lái)的推力損失由飛機(jī)公司承擔(dān)[2]。在這個(gè)約定下，ΦPost被忽略，就可以采用試驗(yàn)測(cè)試所獲得的速度系數(shù)用于求解飛行狀態(tài)下的推力損失。此時(shí)試驗(yàn)測(cè)試的推力系數(shù)可表達(dá)為：

(8)

所以在飛行狀態(tài)時(shí)，發(fā)動(dòng)機(jī)的總凈推力FN為式(9)，而氣動(dòng)阻力相應(yīng)為式(10)。其中Cv由地面測(cè)試而來(lái)，而FNPF采用全機(jī)合力測(cè)量獲得。

(9)

D=FN-FNPF

(10)

從速度系數(shù)的表達(dá)式可知，發(fā)動(dòng)機(jī)噴管的推力損失已經(jīng)在速度系數(shù)中反映出來(lái)了。式(6)給出的是內(nèi)外涵道混合的速度系數(shù)，而在真實(shí)發(fā)動(dòng)機(jī)設(shè)計(jì)過(guò)程中，一般還需要針對(duì)內(nèi)、外涵道分別計(jì)算各自的速度系數(shù)，用于判斷各自的推力損失。此外，進(jìn)氣道和噴管的流量系數(shù)也是發(fā)動(dòng)機(jī)設(shè)計(jì)過(guò)程的重要參數(shù)，它可以表示不同來(lái)流速度下外流對(duì)于噴管的流量抑制效應(yīng)，對(duì)于飛行試驗(yàn)的流量修正具有重要作用。但對(duì)于CFD計(jì)算而言，實(shí)際流量可以通過(guò)計(jì)算網(wǎng)格進(jìn)行精確積分，因此對(duì)于推阻分解的精度影響相對(duì)較小。

綜合上述分析，基于CFD方法的推阻分解過(guò)程分為三個(gè)步驟：首先進(jìn)行地面狀態(tài)分析，獲得噴管性能參數(shù)(流量系數(shù)和速度系數(shù))；然后進(jìn)行飛行狀態(tài)的全機(jī)推阻合力分析；最后將地面狀態(tài)的性能參數(shù)用于推阻分解計(jì)算，最終獲得氣動(dòng)阻力。

2 算例及應(yīng)用

推阻分解方法可以應(yīng)用于進(jìn)/排氣系統(tǒng)的阻力分析，如進(jìn)氣道溢流阻力和噴管噴流推力損失分析，以及機(jī)翼與動(dòng)力短艙之間的干擾評(píng)估等。本節(jié)采用三個(gè)算例對(duì)本文方法進(jìn)行驗(yàn)證和測(cè)試。第一個(gè)是NACA 1-85-75進(jìn)氣道算例，驗(yàn)證方法對(duì)于單獨(dú)進(jìn)氣道溢流阻力的計(jì)算效果；第二個(gè)算例是PAW2雙涵道噴管，分析本文計(jì)算方法對(duì)于噴管性能預(yù)測(cè)的效果；第三個(gè)算例是NAL-AERO-02-01動(dòng)力短艙算例，包含了進(jìn)/排氣性能，對(duì)基于CFD計(jì)算的推阻分解方法進(jìn)行驗(yàn)證和分析。

2.1 NACA 1-85-75進(jìn)氣道算例

進(jìn)氣道溢流阻力的分析是機(jī)體/動(dòng)力一體化集成的重要方面，主要關(guān)心進(jìn)氣道在不同流量系數(shù)下的進(jìn)氣道阻力性能，在動(dòng)力集成中用于不同發(fā)動(dòng)機(jī)推力狀態(tài)之間的性能修正。張美紅等[13]和張兆等[14]采用CFD方法進(jìn)行溢流阻力計(jì)算分析，可作為溢流阻力計(jì)算方法的參考。本文為了測(cè)試推阻分解方法在溢流阻力計(jì)算上的效果，采用NACA 1-85-75進(jìn)氣道進(jìn)行驗(yàn)證。NACA 1系列進(jìn)氣道是NACA早期設(shè)計(jì)的高臨界速度的進(jìn)氣道[15]，擁有較多的試驗(yàn)數(shù)據(jù)[16-17]可供CFD計(jì)算參考，其編號(hào)中的1表示系列，85表示喉道直徑為外罩最大直徑的85%，75表示從唇口到最大直徑處的長(zhǎng)度為最大直徑的75%。

由于NACA 1-85-75進(jìn)氣道算例是軸對(duì)稱外形，且不考慮攻角變化，因此可以采用軸對(duì)稱方法進(jìn)行計(jì)算。計(jì)算模型的最大直徑為45.72 cm，長(zhǎng)度為34.29 cm，喉道直徑為38.86 cm，唇口與喉道的面積之比為1.009，具體模型外形和安裝方式可見文獻(xiàn)[16]。計(jì)算網(wǎng)格如圖2所示，二維截面網(wǎng)格總數(shù)約76 000點(diǎn)。為了減小邊界條件對(duì)計(jì)算結(jié)果的影響，上游遠(yuǎn)場(chǎng)邊界與唇口的距離為500 cm，下游進(jìn)氣道的出口采用等直段延伸100 cm。壁面第一層網(wǎng)格的法向距離為0.0002 cm，以保證第一層的y+小于1。由于模型具有軸對(duì)稱的特點(diǎn)，圓周方向僅計(jì)算3°角的扇形區(qū)域，網(wǎng)格采用二維網(wǎng)格旋轉(zhuǎn)而成。

圖4顯示了幾個(gè)不同流量系數(shù)下的馬赫數(shù)分布云圖和流線圖?？梢钥闯觯?dāng)φ=0.2時(shí)，短艙外壁面已經(jīng)出現(xiàn)明顯的流動(dòng)分離；隨著流量系數(shù)的增加，φ=0.4時(shí)流動(dòng)分離消失，但由于溢流效應(yīng)的影響，表面出現(xiàn)明顯的激波；流量系數(shù)繼續(xù)增加，激波削弱，直至φ=0.8時(shí)激波和流動(dòng)分離均消失，流場(chǎng)特性較好。

為了定量分析不同來(lái)流速度下的溢流阻力特性，本文分別計(jì)算了Ma=0.60、0.85、0.94和1.20四個(gè)速度下的不同流量系數(shù)，基于進(jìn)氣道外罩最大直徑的雷諾數(shù)分別為4.5×106、5.4×106、5.5×106、5.7×106。

圖5給出了不同馬赫數(shù)下φ=0.60的短艙外罩表面壓力分布與試驗(yàn)結(jié)果的對(duì)比。從圖上可以看出，不同馬赫數(shù)下壓力分布均與試驗(yàn)結(jié)果吻合良好。在Ma=0.60時(shí)，由于溢流造成的表面加速，使進(jìn)氣道前緣出現(xiàn)了明顯的吸力峰值；Ma=0.85和0.94時(shí)外壁面都存在明顯的激波。

溢流阻力是進(jìn)氣道設(shè)計(jì)關(guān)注的重要性能，其定義是按照流量系數(shù)為1時(shí)溢流阻力為0，其他流量系數(shù)下的阻力增量即為溢流阻力。這就表示流量系數(shù)為1時(shí)的阻力系數(shù)是一個(gè)重要的參考量，其預(yù)測(cè)精度對(duì)進(jìn)氣道設(shè)計(jì)也有重要意義。文獻(xiàn)[10,13-14]給出了溢流阻力的計(jì)算方法，但對(duì)于阻力系數(shù)絕對(duì)值的計(jì)算效果并未驗(yàn)證。在本文的計(jì)算中采用推阻分解方法，

利用式(4)可以獲得進(jìn)氣道阻力，結(jié)果如圖6所示，圖上給出了不同馬赫數(shù)下的進(jìn)氣道阻力系數(shù)與試驗(yàn)的對(duì)比，其參考面積為進(jìn)氣道的外罩最大截面積。從計(jì)算結(jié)果可以看出，采用本文的推阻分解計(jì)算方法，不僅能夠準(zhǔn)確預(yù)測(cè)進(jìn)氣道隨流量系數(shù)變化的溢流阻力，其進(jìn)氣道阻力系數(shù)的絕對(duì)值也與試驗(yàn)符合良好。

2.2 DSFR噴管驗(yàn)證算例

噴管性能是機(jī)體/推進(jìn)系統(tǒng)集成所關(guān)注的另一重要方面，著重關(guān)注噴管在地面狀態(tài)下不同總壓的推力系數(shù)，獲得推力系統(tǒng)隨總壓的變化曲線，在飛行狀態(tài)下需要利用該性能曲線進(jìn)行推力與阻力的分解。因此噴管性能曲線的計(jì)算精度對(duì)于飛行推力的預(yù)測(cè)精度有直接的影響。AIAA自2011年起組織了三屆PAW研討會(huì)，旨在探討CFD方法對(duì)于推阻分解的計(jì)算精度。第二屆和第三屆的噴管算例是典型的高涵道比民用航空發(fā)動(dòng)機(jī)DSFR噴管(Dual Separate Flow Reference Nozzle)[11-12]。本文為了驗(yàn)證和演示推阻分解方法在噴管性能的計(jì)算過(guò)程，采用軸對(duì)稱的高涵道比噴管進(jìn)行CFD驗(yàn)證。

圖7顯示了噴管算例的模型、邊界條件與計(jì)算網(wǎng)格，其中坐標(biāo)的單位為英寸，與原始算例一致。噴管入口采用總溫、總壓邊界條件，遠(yuǎn)場(chǎng)上游采用速度進(jìn)口邊界條件，下游采用無(wú)反射邊界條件，壁面采用絕熱固壁條件。為了更真實(shí)地模擬發(fā)動(dòng)機(jī)試車臺(tái)的靜止外流條件，但同時(shí)也避免CFD的計(jì)算穩(wěn)定性問(wèn)題，外流需要采用一個(gè)極低速度的來(lái)流條件。該算例中外涵道總壓始終為內(nèi)涵道總壓的1.2倍，內(nèi)、外涵道總溫均為環(huán)境溫度的1.019倍。計(jì)算采用三套不同的計(jì)算網(wǎng)格測(cè)試網(wǎng)格敏感性，研究網(wǎng)格密度對(duì)于計(jì)算結(jié)果的影響，網(wǎng)格點(diǎn)數(shù)分別為87 500、172 000、345 000，不同網(wǎng)格之間的加密方式為x和y方向各加密1.4倍，這種加密方式與AIAA所舉辦的阻力預(yù)測(cè)研討會(huì)(DPW)所推薦的網(wǎng)格加密方式[18]相似。

圖7顯示了計(jì)算所采用的稀網(wǎng)格，網(wǎng)格在噴管入口進(jìn)行了適當(dāng)加密。采用軸對(duì)稱方式計(jì)算，圓周方向采用1個(gè)網(wǎng)格單元。為了保持x、y、z三個(gè)方向網(wǎng)格的密度匹配，從稀網(wǎng)格、中網(wǎng)格到密網(wǎng)格的軸對(duì)稱計(jì)算扇形角度分別為4°、3°、2°。空間離散采用Roe迎風(fēng)格式，湍流模式采用SST模式。

圖8顯示了不同壓比下的流線分布和馬赫數(shù)云圖，這是外流馬赫數(shù)為0.01時(shí)的密網(wǎng)格計(jì)算結(jié)果。從流線看，射流對(duì)于外流有明顯的引射作用，外場(chǎng)流線向射流核心區(qū)集中。當(dāng)外涵道總壓低于1.89時(shí)，整個(gè)流道速度低于聲速；當(dāng)壓比高于1.89時(shí)，流場(chǎng)出現(xiàn)超聲速區(qū)域，并出現(xiàn)激波。

前文已經(jīng)分析過(guò)，推阻分解中噴管性能系數(shù)是評(píng)價(jià)推力損失的重要參數(shù)，需要利用地面狀態(tài)的性能參數(shù)評(píng)估飛行狀態(tài)的推力損失，因此其計(jì)算精度對(duì)結(jié)果影響十分重要。波音公司的Wright[19]的誤差估計(jì)顯示，對(duì)于速度系數(shù)為0.99的噴管，其速度系數(shù)0.001的誤差就將引起飛機(jī)全機(jī)1 count(即阻力系數(shù)0.0001)的阻力系數(shù)誤差。

(11)

(12)

稀網(wǎng)格的風(fēng)扇噴管流量系數(shù)稍有差別，其余結(jié)果具有很好的網(wǎng)格收斂性。計(jì)算網(wǎng)格引起的速度系數(shù)的誤差小于0.0005。從計(jì)算結(jié)果看，中等網(wǎng)格計(jì)算結(jié)果和密網(wǎng)格計(jì)算結(jié)果幾乎完全一致，由于密網(wǎng)格計(jì)算資源消耗較大，所以選擇中等網(wǎng)格密度作為噴管性能分析的參考。

在等熵流動(dòng)假設(shè)下，壓比1.89時(shí)喉道就達(dá)到聲速，出現(xiàn)堵塞，但由于邊界層的作用，真實(shí)噴管達(dá)到堵塞的壓比更高。圖9顯示當(dāng)風(fēng)扇噴管的壓比達(dá)到約2.2時(shí)才趨于堵塞，此時(shí)流量系數(shù)和速度系數(shù)均較高，所以是典型民用航空發(fā)動(dòng)機(jī)的典型設(shè)計(jì)狀態(tài)。

在發(fā)動(dòng)機(jī)地面試驗(yàn)中，通過(guò)測(cè)試設(shè)備將短艙外罩速度近似給定為0，但是在可壓縮CFD分析中，大范圍速度為0的計(jì)算域帶來(lái)計(jì)算魯棒性問(wèn)題，所以需要在外場(chǎng)區(qū)域設(shè)定一個(gè)較小的流動(dòng)速度。但外流流動(dòng)會(huì)給短艙外罩帶來(lái)阻力，從而使推力損失增加，速度系數(shù)降低。下面測(cè)試外場(chǎng)速度對(duì)計(jì)算精度和數(shù)值穩(wěn)定性的影響。

圖10給出了外場(chǎng)馬赫數(shù)為Ma∞=0.01、0.02和0.05時(shí)中等網(wǎng)格的流量系數(shù)和速度系數(shù)與試驗(yàn)結(jié)果的對(duì)比?？梢钥闯?，在這個(gè)馬赫數(shù)范圍內(nèi)，核心噴管和風(fēng)扇噴管的流量系數(shù)變化很小，尤其是核心噴管，流量系數(shù)幾乎不受外流速度的影響。從流量系數(shù)對(duì)比結(jié)果看，核心噴管的誤差最大達(dá)到1%，但由于核心噴管流量很小，其誤差權(quán)重也相對(duì)較??；風(fēng)扇噴管的流量系數(shù)誤差約為0.2%，與試驗(yàn)吻合較好。不同外流馬赫數(shù)結(jié)果顯示，速度系數(shù)明顯受到外流速度的影響，這是由于外部流動(dòng)在短艙外罩上造成的損失所致。結(jié)果顯示，計(jì)算和試驗(yàn)的速度系數(shù)的變化趨勢(shì)非

圖9不同網(wǎng)格的流量系數(shù)和速度系數(shù)與試驗(yàn)對(duì)比
Fig.9Comparisonsofdischargecoefficientsandvelocitycoefficient

圖10不同外流速度下的流量系數(shù)和速度系數(shù)與試驗(yàn)結(jié)果對(duì)比
Fig.10Comparisonsofdischargecoefficientsandvelocitycoefficientofdifferentexternalflowspeeds

常一致，這表明流動(dòng)現(xiàn)象被準(zhǔn)確捕捉。Ma∞=0.01、0.02和0.05的計(jì)算結(jié)果表明，外流速度越小，速度系數(shù)越大，當(dāng)速度趨于0時(shí)趨于一個(gè)收斂的值。在真實(shí)發(fā)動(dòng)機(jī)分析時(shí)往往需要根據(jù)不同外流速度下的速度系數(shù)外插出速度為0時(shí)的值。在本文計(jì)算結(jié)果中，不同速度下的Cv差別大約為0.002，考慮試驗(yàn)的誤差帶，計(jì)算與試驗(yàn)值的誤差大約0.003～0.004左右，并且顯示為整體的系統(tǒng)誤差，便于進(jìn)行修正，所以整體來(lái)看計(jì)算與試驗(yàn)有較好的吻合。

2.3 NAL-AERO-02-01 TPS驗(yàn)證算例

針對(duì)機(jī)體/動(dòng)力系統(tǒng)的集成問(wèn)題，推力與阻力的分解既要分析進(jìn)氣道性能，又要準(zhǔn)確模擬噴管特性。推阻分解過(guò)程需要有三個(gè)步驟：首先計(jì)算發(fā)動(dòng)機(jī)短艙在地面狀態(tài)的性能系數(shù)，然后再計(jì)算飛行狀態(tài)下的氣動(dòng)性能，最后利用推阻分解方法對(duì)推力與阻力進(jìn)行后處理。針對(duì)機(jī)體與動(dòng)力短艙的干擾問(wèn)題，需要先分析獨(dú)立短艙在地面狀態(tài)的性能特性，然后分析短艙/吊掛/機(jī)翼構(gòu)型的性能，獲得各構(gòu)型在地面狀態(tài)下由于干擾量所帶來(lái)的推力損失，最后利用推阻分解方法獲得推力和阻力?？偟膩?lái)說(shuō)，推阻分解計(jì)算結(jié)果的精確程度受到兩方面影響，一是標(biāo)定狀態(tài)的速度系數(shù)精度，二是飛行狀態(tài)的全機(jī)阻力預(yù)測(cè)精度。

本節(jié)以一個(gè)軸對(duì)稱TPS (Turbine Powered Simulator)短艙模型NAL-AERO-02-01模型[20]為例，演示利用推阻分解方法計(jì)算帶動(dòng)力短艙的氣動(dòng)力過(guò)程。為了與真實(shí)民用航空發(fā)動(dòng)機(jī)的短艙更為接近，本文將原始模型進(jìn)行了放大，使其基于短艙外罩的雷諾數(shù)達(dá)到1.5×107(原始模型為2.0×106)，與典型窄體客機(jī)的發(fā)動(dòng)機(jī)的雷諾數(shù)相當(dāng)，其余計(jì)算參數(shù)如馬赫數(shù)、噴流總溫總壓等與文獻(xiàn)[20]相同。

圖11顯示了TPS短艙的計(jì)算網(wǎng)格。進(jìn)氣道唇口處采用C型網(wǎng)格準(zhǔn)確捕捉進(jìn)氣道的溢流流動(dòng)，噴管網(wǎng)格進(jìn)行加密以分辨高壓流體的膨脹過(guò)程，壁面法向第一層網(wǎng)格的y+約為1，采用軸對(duì)稱方法計(jì)算，二維截面網(wǎng)格總數(shù)為7.0×104。采用ROE格式和SA湍流模式進(jìn)行CFD分析。

以跨聲速M(fèi)a∞=0.801的流動(dòng)狀態(tài)為例分析推阻分解的過(guò)程。以遠(yuǎn)場(chǎng)溫度、壓力作為參考量，風(fēng)扇噴管的無(wú)量綱總溫、總壓分別為1.251和2.049，核心機(jī)噴管的總溫、總壓分別為0.690和1.406，進(jìn)氣道的流量系數(shù)為0.497。計(jì)算參考面積選取短艙外罩的最大截面積。為了獲得地面試車狀態(tài)下的噴管性能參數(shù)，首先以來(lái)流馬赫數(shù)0.005、0.01、0.02三個(gè)狀態(tài)進(jìn)行短艙凈推力分析，按照式(8)進(jìn)行速度系數(shù)計(jì)算，三個(gè)狀態(tài)的速度系數(shù)分別為0.9887，0.9886和0.9884，其中來(lái)流馬赫數(shù)0.005最為接近地面靜止?fàn)顟B(tài)，采用外插方法確定此壓比下的噴管速度系數(shù)為0.9888。然后計(jì)算Ma=0.801時(shí)的流場(chǎng)，采用流量自動(dòng)調(diào)節(jié)的方式獲得指定的流量系數(shù)。圖12顯示了來(lái)流馬赫數(shù)為0.005和0.801時(shí)兩個(gè)流場(chǎng)的對(duì)比。其中在模擬地面測(cè)試狀態(tài)時(shí)為了排除進(jìn)氣道對(duì)參數(shù)標(biāo)定的影響，在風(fēng)扇面采用了壁面條件。從圖12可以看出，兩個(gè)狀態(tài)的噴管流場(chǎng)非常相似，但由于外流速度不同，噴流對(duì)外部流體的卷吸現(xiàn)象差別較大，導(dǎo)致外場(chǎng)流線差別較大。

圖13給出了Ma=0.801狀態(tài)下的壓力分布計(jì)算與文獻(xiàn)[20]試驗(yàn)值的對(duì)比，其中橫坐標(biāo)為相對(duì)坐標(biāo)，C為短艙外罩的長(zhǎng)度。雖然本文的計(jì)算雷諾數(shù)更大，但壓力分布與試驗(yàn)也符合較好。

表1給出了根據(jù)CFD計(jì)算和推阻分解方法獲得的計(jì)算結(jié)果。表1中進(jìn)氣道和噴管的流量以及各控制面的受力均基于CFD直接積分獲得，推力和阻力采用推阻分解的方法獲得。其中推力采用式(9)計(jì)算，當(dāng)Cv=1時(shí)為理想推力，當(dāng)Cv=0.9888即為扣除內(nèi)阻損失后的實(shí)際推力。短艙所受合力與理想推力的差量即為短艙所受的總阻力。

表1 Ma=0.801時(shí)推阻積分各受力分量分解(無(wú)量綱)Table 1 Thrust and drag components of Ma=0.801 case (non-dimensional)

當(dāng)需要進(jìn)行不同狀態(tài)下的推阻分解分析時(shí)，若噴管的總溫、總壓沒(méi)有發(fā)生變化，則可以用同樣的速度系數(shù)計(jì)算推力，倘若噴管狀態(tài)發(fā)生變化，就需要重新標(biāo)定該壓比下的地面試車狀態(tài)的性能。表2給出了幾個(gè)不同狀態(tài)下的推阻分解結(jié)果，其中前5個(gè)狀態(tài)是固定噴管流動(dòng)狀態(tài)、改變來(lái)流速度，最后1個(gè)狀態(tài)改變了噴管流動(dòng)狀態(tài)。

表2 不同計(jì)算狀態(tài)下的短艙推阻分解結(jié)果(無(wú)量綱)Table 2 Thrust/drag bookkeeping results of different flow conditions (non-dimensional)

從變馬赫數(shù)的結(jié)果可以看出，由于以來(lái)流速度作為參考值，因此無(wú)量綱的發(fā)動(dòng)機(jī)理想推力隨來(lái)流速度減小越來(lái)越大，相應(yīng)的短艙內(nèi)阻也越來(lái)越大。而進(jìn)行推阻分解后的短艙外阻隨馬赫數(shù)增加先減小后增加。從圖14的不同馬赫數(shù)下的壓力分布對(duì)比可以看出，在馬赫數(shù)較小時(shí)沒(méi)有激波，并且前緣吸力峰值較低；當(dāng)馬赫數(shù)逐漸提高，短艙前緣的吸力峰值提高，因此能夠產(chǎn)生一定的推力分量，使短艙外阻減?。划?dāng)馬赫數(shù)繼續(xù)增加，短艙外罩出現(xiàn)激波，并且激波逐漸增強(qiáng)，所以短艙外阻增加。

表2的最后一個(gè)狀態(tài)是改變噴管流動(dòng)狀態(tài)的結(jié)果，此時(shí)需要重新標(biāo)定噴管性能參數(shù)用于推阻分解分析。其中狀態(tài)6與狀態(tài)4相比發(fā)動(dòng)機(jī)推力增大，噴流與外流之間的干擾增大，所以短艙外阻有所增大。

3 總 結(jié)

本文介紹了一種基于CFD的推阻分解方法，可用于帶動(dòng)力短艙計(jì)算時(shí)的推阻分解。推阻分解方法廣泛應(yīng)用于飛機(jī)公司與發(fā)動(dòng)機(jī)公司，其主要原理是采用流體力學(xué)的控制體理論，針對(duì)推力和阻力產(chǎn)生的各自區(qū)域進(jìn)行受力分解，然后利用風(fēng)洞測(cè)試或者CFD分析確定各項(xiàng)阻力分量。

本文首先分析了推阻分解方法的控制體理論，明確各項(xiàng)推力與阻力分量的計(jì)算方法。然后采用三個(gè)算例對(duì)計(jì)算方法進(jìn)行詳細(xì)分析。其中進(jìn)氣道算例用于驗(yàn)證推阻分解方法對(duì)于溢流阻力的計(jì)算效果；噴管算例用于驗(yàn)證和分析速度系數(shù)和流量系數(shù)的計(jì)算精度；最后是一個(gè)完整的TPS短艙分析，用于演示推阻分解方法對(duì)于動(dòng)力短艙阻力分量的計(jì)算過(guò)程。三個(gè)不同算例的結(jié)果表明，本文CFD計(jì)算結(jié)果與試驗(yàn)有較好的吻合，并且推阻分解過(guò)程能夠準(zhǔn)確分辨推力與阻力分量。

參 考 文 獻(xiàn)：

[1]Rooney E C. Thrust and drag: Its prediction and verification, Chapter 2, Thrust-Drag Accounting Methodology[M]//Covert E E, edited. Progress in Astronautics and Aeronautics, vol. 98. American Institute of Aeronautics and Astronautics, 1985.

[2]SAE Technical Standards. In-flight thrust determination, appendix c thrust and drag accounting, and appendix e case studies and applications[R]. SAE Aerospace Information Report AIR1703 rev. A, 2012.

[3]MIDAP Study Group. Guide to in-flight thrust measurement of turbojets and fan engines, Chapter 2 Propulsion system thrust and drag bookkeeping[R].Advisory Group for Aerospace Research and Development. AGARD-AG-237, 1979: 27-62.

[4]VonGeyr H F, Rossow C C. A correct thrust determination method for turbine powered simulators in-wind tunnel testing[R]. AIAA 2005-3707.

[5]Hunt D N. Experimental techniques used to evaluate propulsion system interference effects on the cruise configuration of the boeing C-14[R]. AIAA 1979-0335, 1979.

[6]Dusa D, Lahti D J, Berry D. Investigation of subsonic nacelle performance improvement concept[R]. AIAA 1982-1042, 1982.

[7]Bousquet J M. Survey of engine integration testing in ONERA wind tunnels[R]. AIAA 2005-3705, 2005.

[8]Winkler C M, Dorgan A J. BCFD analysis for the 1st AIAA propulsion workshop: Nozzle results[R]. AIAA 2013-3731.

[9]Zhang Y F, Chen H X, Zhang M, et al. Performance prediction of conical nozzle using Navier-Stokes computation[J]. Journal of Propulsion and Power, 2015; 31(1): 192-203.

[10]Spotts N, Guzik S, Gao X F. A CFD analysis of compressible flow through convergent-conical nozzles[R]. AIAA 2013-3734.

[11]Mikkelsen K L, Myren D J, Dahl D G, et al. Initial subscale performance measurements of the AIAA dual separate flow reference (DSFR) Nozzle[R]. AIAA 2015-3883.

[12]Li Z, Chen H X, Zhang Y F. NSAWET results of the dual separate flow reference nozzle from AIAA PAW02[R]. AIAA 2015-3779.

[13]張美紅, 張冬云, 王美黎, 等. 基于CFD和推阻分解技術(shù)的全機(jī)溢流阻力預(yù)測(cè)與分析[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2016, 34(5): 625-630.

[14]張兆, 陶洋, 黃國(guó)川. 發(fā)動(dòng)機(jī)短艙溢流阻力的數(shù)值模擬[J]. 航空學(xué)報(bào), 2013, 34(3): 547-553.

[15]Baals D D, Smith N F, Wright J B. The development and application of high-critical-speed nose inlets[R]. Langley Aeronautical Laboratory National Advisory Committee for Aeronautics, NACA TR-920, 1948.

[16]Richard J R. An investigation of several NACA 1-series axisymmetric inlets at Mach numbers from 0.4 to 1.29[R]. NASA TM X-2917, 1974.

[17]Richard J R, William K A. A wind tunnel investigation of three NACA 1-series inlets at Mach numbers up to 0.92[R]. NASA TM-110300, 1996.

[18]Mavriplis D J, Vassberg J C, Tinoco E N, et al. Grid quality and resolution issues from the drag prediction workshop series[R]. AIAA 2008-930.

[19]Wright F L. Comparison of least squares curve fit and individual sample statistical analysis results of calibration data for the velocity coefficient of a flow nacelle[R]. AIAA 1994-2587, 1994.

[20]Hirose N, Asai K, Ikawa K. Transonic 3-D Euler analysis of flows around fanjet engine and TPS (Turbine Powered Simulator), comparison with wind tunnel experiment, evaluation of TPS testing method and 3-D flow[R]. National Aerospace Laboratory, TR-1045, 1989.