超聲速全機帶動力短艙對近場壓力信號和地面聲爆的影響

林榕婷，譚廉華，霍 滿，吳宇昂，杜 璽，林大楷

(1.中國商用飛機有限責(zé)任公司北京民用飛機技術(shù)研究中心，北京 102211;2.民用飛機設(shè)計數(shù)字仿真技術(shù)北京市重點實驗室，北京 102211)

引 言

20世紀(jì)60年代以來，圖144[1]、協(xié)和號飛機[2]的出現(xiàn)開啟了超聲速商業(yè)飛行的時代。協(xié)和號因出色的氣動性能造就了世界航空史的里程碑[3]，但最終卻由于超聲速聲爆、油耗等一系列問題，于2003年退役。航空界對超聲速商用飛機的研究從未停止，特別是對環(huán)保低聲爆的研究。近年來主要航空強國采取由政府科研機構(gòu)主導(dǎo)、一些創(chuàng)業(yè)公司積極參與等方式，掀起超聲速商用飛機研發(fā)熱潮，旨在提高運營經(jīng)濟性，實現(xiàn)洲際、跨洋航線上低聲爆超聲速巡航[4-6]。日本宇航研究開發(fā)機構(gòu)JAXA的S4靜音超聲速客機研究項目對低聲爆關(guān)鍵技術(shù)展開了一系列試驗和理論研究[7]。NASA制定的“N+3”計劃[8]列出了超聲速客機的聲爆技術(shù)指標(biāo)。近年來，NASA緊鑼密鼓推進(jìn) X-59“靜音超聲速技術(shù)”(X-59 QueSST)[9]。Aerion Supersonic公司的AS2超聲速公務(wù)機已完成相關(guān)風(fēng)洞試驗[10]，Boom Supersonic公司的超聲速飛機XB-1，已于2020年10月完成整機安裝[11]。俄羅斯茹科夫斯基研究所曾計劃在2021年前完成超聲速客機概念方案設(shè)計工作[12]。我國目前正迎來從制造大國走向制造強國、從航空大國走向航空強國的新時代，發(fā)展國產(chǎn)超聲速商用飛機意義重大。低聲爆氣動設(shè)計是超聲速飛機設(shè)計的關(guān)鍵技術(shù)[13]，帶動力短艙是聲爆的重要影響因素之一，直接關(guān)系平面布局的優(yōu)化等諸多方面。

美國航空航天學(xué)會(AIAA)組織的聲爆預(yù)測研討會迄今為止舉辦了3屆。會議提供了眾多驗證標(biāo)模，NASA，Boeing，Airbus，Ansys，西北工業(yè)大學(xué)等眾多國內(nèi)外知名航空研究單位和制造商參與并針對標(biāo)模計算展開相關(guān)的技術(shù)交流，研究了諸如黏性、網(wǎng)格密度、空間離散格式、湍流模型、噴流等對超聲速聲爆的影響。超聲速噴口Biconvex[14-15]、帶動力短艙的NASA C25D[16]等標(biāo)模為參會者研究噴流對聲爆的影響提供了重要參考。近年來，國內(nèi)研究單位包括中航工業(yè)氣動院[17]、中國航空研究院等針對低聲爆設(shè)計也開展了眾多試驗和數(shù)值研究。中航工業(yè)氣動院于FL-60風(fēng)洞采用Seeb-ALR低聲爆標(biāo)模和自行設(shè)計的帶噴流的旋成體模型開展并成功完成近場壓力特征測量的驗證試驗[18]。西北工業(yè)大學(xué)的王剛等[19]研究了不同空間離散格式、模型尖點精度、黏性對近場壓力信號和遠(yuǎn)場聲爆預(yù)測結(jié)果的影響。喬建領(lǐng)等[20]開發(fā)了可考慮大氣風(fēng)效應(yīng)的遠(yuǎn)場高精度聲爆預(yù)測程序，研究了大氣風(fēng)對聲爆的影響等。目前，國內(nèi)針對黏性、網(wǎng)格密度、聲爆預(yù)測算法對超聲速聲爆預(yù)測的研究較多，但針對動力短艙對近場壓力信號和地面聲爆的影響研究較少，尤其是在超聲速商用飛機研究方向。

中國商飛北研中心自2016年開始超聲速客機的相關(guān)研究，開展了針對超聲速流場的CFD工具驗證、全機復(fù)雜流場CFD數(shù)值模擬[21-22]、地面聲爆信號測試、超聲速自然層流機翼氣動設(shè)計等工作。對標(biāo)模NASA C25D的帶動力短艙對超聲速飛機近場壓力信號和地面聲爆的影響開展了一系列計算研究和規(guī)律探索，為超聲速商用飛機低聲爆設(shè)計做了鋪墊。本文共分為3部分:模型和計算方法介紹、計算分析、結(jié)論。

1 模型及方法

1.1 模型及網(wǎng)格

采用第2屆聲爆預(yù)測研討會的低聲爆驗證標(biāo)模NASA C25D，分別對帶通氣短艙(以下簡稱C25F)和帶動力短艙(以下簡稱C25P)兩種構(gòu)型開展數(shù)值計算研究。

聲爆預(yù)測研討會提供的幾何模型圖如圖1，2所示。模型的參考面積半模為37.16 m2，參考長度即機身長度為32.92 m。

圖1 NASA C25D概念圖[16]Fig.1 NASA C25D concept graph

(a) Diagram of C25F

采用混合網(wǎng)格，模型自帶3.375°攻角。內(nèi)域為圓柱體區(qū)域的非結(jié)構(gòu)化網(wǎng)格，外域為結(jié)構(gòu)化網(wǎng)格，以Mach錐角為傾斜角度，5倍參考長度為縱向高度，外域網(wǎng)格的增長率為1.025。針對通氣短艙和動力短艙，分別繪制無黏的Euler網(wǎng)格和有黏RANS計算網(wǎng)格，每種類型制作3種不同密度的網(wǎng)格，密度化主要是針對內(nèi)域非結(jié)構(gòu)網(wǎng)格。內(nèi)域網(wǎng)格，對于黏性邊界層，首層網(wǎng)格為8.0×10-5m，增長率為1.15，50層。總網(wǎng)格單元數(shù)如表1所示。為驗證網(wǎng)格無關(guān)性并確保計算的準(zhǔn)確程度，提取了計算結(jié)果中的重要氣動參數(shù)升力系數(shù)CL、阻力系數(shù)CD作為指標(biāo)，數(shù)據(jù)如表2所示?？梢钥吹剑珻L的差距在1.0×10-4量級，CD的差距在1.0×10-5量級，因此可認(rèn)為各自3套網(wǎng)格的計算結(jié)果在可忽略的精度范圍內(nèi)，均符合網(wǎng)格無關(guān)性要求。整體和局部截圖如圖3～5所示。

表1 NASA C25D網(wǎng)格數(shù)目Table 1 Grid numer of NASA C25D

表2 力系數(shù)的精度對比Table 2 Accuracy comparison of force coefficient

圖3 網(wǎng)格整體分布Fig.3 Global diagram of grids

(a) C25F invicid(grid No.3)

(a) C25P invicid(grid No.3)

1.2 計算方法

表3是C25F的計算工況條件，C25P的計算條件除了表3的基本參數(shù)之外，還包含發(fā)動機進(jìn)口和出口的壓比和溫比等參數(shù)，如表4所示。表中，P*為總壓，T*為總溫，a為聲速，V為速度，下標(biāo)inf表示來流。

表3 C25F工況信息(海拔15.76 km)Table 3 C25F case information (altitude 15.76 km)

表4 C25P工況信息(海拔15.76 km)Table 4 C25P case information (altitude 15.76 km)

流場計算方面，對C25F和C25P均分別進(jìn)行基于Euler的定常計算和基于RANS的定常計算。時間推進(jìn)為隱式格式，空間離散采用TVD格式，使用 Minmod限制器。湍流模型采用一方程SA模型。邊界條件的設(shè)置:計算域的遠(yuǎn)場入口采用超聲速入口條件，給定靜壓、密度和速度;遠(yuǎn)場出口采用超聲速出口，無需給定條件，參數(shù)由內(nèi)場直接外推得到;動力短艙的發(fā)動機入口給定背壓，出口給定發(fā)動機出口的總溫總壓。遠(yuǎn)場壓力信號的計算是基于線化理論和幾何聲學(xué)波形參數(shù)法的自主開發(fā)工具，得到地面壓力信號，進(jìn)而求解以分貝(dB)為單位的地面可感知聲音響度PL值，與參會者的計算結(jié)果進(jìn)行對比。

2 計算分析

2.1 流場分布

2.1.1 結(jié)果展示

流場計算直接影響近場壓力信號的提取和地面聲爆的計算。圖6～10分別是C25F和C25P兩種構(gòu)型，在Euler計算和RANS計算下，由各自的密網(wǎng)格，即第3套網(wǎng)格(網(wǎng)格密度接近，網(wǎng)格單元總數(shù)均在半模5×107左右)計算得到的流場參數(shù)云圖對比圖，涵蓋了:Mach數(shù)、溫度、壓力、壓力信號、表面壓力系數(shù)等。壓力信號定義為

dp/pinf=(p-pinf)/pinf

(a) C25F invicid(grid No.3)

(a) C25F invicid(grid No.3)

(a) C25F invicid(grid No.3)

(a) C25F invicid(grid No.3)

(a) Main wing surface at 50% span of flat tail

式中，pinf為來流靜壓，p為當(dāng)?shù)仂o壓。

如圖9所示，選取4種組合各自的第3套網(wǎng)格，分別在主翼約60%展長處和平尾約50%展長處提取了機翼表面的壓力系數(shù)，這里的壓力系數(shù)的定義是

式中，uinf為來流速度。

圖10是主翼和平尾在這兩處的上下表面壓力系數(shù)曲線。圖中BATRI代表中國商飛北研中心。

2.1.2 規(guī)律分析

(1)黏性對流場的影響

有黏計算由于壁面采用無滑移邊界條件，壁面處速度為零，相應(yīng)地，壁面處溫度和壓力均升高。黏性計算的機翼上同位置的機體溫度大約是遠(yuǎn)場來流的1.5倍，而無黏計算由于沒有黏性耗散，參數(shù)在壁面的呈現(xiàn)與周邊過渡均勻。從圖6～8可以看到，無黏計算的激波系相比有黏計算更明顯和清晰，尤其是短艙內(nèi)部的激波反射現(xiàn)象得到了較清晰的呈現(xiàn)和保持。黏性計算由于短艙內(nèi)壁邊界層耗散，出口處整體速度低于無黏的出口速度，壓力、溫度均比較高。黏性的影響持續(xù)到了短艙出口外的一段距離，可從尾跡的流動明顯看到，無黏的噴流較能保持均勻傳播，而有黏的出口尾跡逐漸放大，這是黏性耗散使得發(fā)動機噴流更快地減速增溫增壓。從圖9，10的比對，同位置，有黏計算比同構(gòu)型的無黏計算的激波位置更加靠前，且在主翼上表面后半部分，壓力系數(shù)絕對值更小。

(2)動力短艙對流場的影響

C25F與C25P相比，譬如圖8(b)，(d)的對比，C25P在發(fā)動機進(jìn)氣的唇口處出現(xiàn)了明顯的溢流現(xiàn)象，唇口上部出現(xiàn)了激波，這是由動力短艙發(fā)動機入口的流量限制而導(dǎo)致的。C25F的通氣短艙由于溢流較少，唇口無明顯激波出現(xiàn)。C25P在發(fā)動機出口和平尾附近的波系比C25F更為復(fù)雜和多樣。結(jié)合圖9(b)，(d) 和圖10(a)，C25P的溢流明顯多于C25F，因此在主翼上表面的激波更提前。在所提取的平尾50%處，對應(yīng)的主翼上表面后緣處出現(xiàn)了明顯的激波區(qū)域，這是短艙入口唇口的溢流導(dǎo)致對前方高速來流的擠壓而出現(xiàn)的。該展長處，C25P平尾的下表面壓力系數(shù)變化率比C25F更劇烈，這跟平尾下表面的速度變化趨勢有關(guān)。C25P由于唇口溢流明顯，在平尾下表面前緣的氣流速度受到溢流影響而比C25F同位置的速度更低;在平尾下表面后緣處，受到噴流速度的影響，氣流的黏性帶動平尾下表面后緣周邊氣流加速流動，壓力比C25F同位置更小。因此，沿著弦向，C25P的平尾下表面加速度更大，壓力系數(shù)曲線的斜率更大，這也直觀地反映在圖10(b)中。從圖10(c)可以看到，主翼上表面的壓力系數(shù)分布差異較大。C25P的動力短艙溢流造成激波提前，反映在云圖中，這道激波的波后高壓區(qū)域面積比C25F更大，因此在曲線圖中，C25P的上表面后緣處能看到明顯的激波，而C25F則基本沒有。

2.2 近場壓力信號

2.2.1 結(jié)果展示

圖11是本小節(jié)近場壓力信號的提取位置和角度示意圖。圖11(a)為計算的提取位置。一般地，我們以參考長度L進(jìn)行丈量，參考長度即為機身長度，選取1倍機身長度和3倍機身長度為半徑R,進(jìn)行等值面繪制。如圖11(b)所示，在不同半徑的等值面上，提取不同角度(這里稱為“滾轉(zhuǎn)角”，用φ表示)的壓力信號，展開分析研究。選取了0°，10°，20°，30°，40°，50°這6種角度，半徑R=H/L=1和3，進(jìn)行壓力信號提取的比對分析。圖12～15分別為計算結(jié)果提取的不同高度、不同滾轉(zhuǎn)角的近場壓力信號與部分參會者計算結(jié)果的比對。圖中，Xn=0為該高度下激波系的起始點，根據(jù)網(wǎng)格起始點和Mach錐角就可以換算得到。從圖12～15可看到，計算結(jié)果與參會者的趨勢是一致的，由于采用不同網(wǎng)格密度，網(wǎng)格越細(xì)密，數(shù)值耗散越小，對近場壓力信號的捕捉也就越精細(xì)，與參會者的壓力信號峰值峰谷等數(shù)據(jù)更接近，更有一部分的計算結(jié)果比參會者捕捉到的近場壓力信號更精細(xì)。下一節(jié)將分析黏性、有無動力、提取位置差異對近場壓力信號的影響。

(a)Extraction height distribution

(a)R=1，φ=0°[24]

(a)R=1，φ=0°[26]

(a)R=1，φ=0°

(a)R=1，φ=0°

2.2.2 規(guī)律分析

(1)黏性對近場信號的影響

流場差異直接導(dǎo)致近場壓力信號的差異，同時波系受機身下部幾何壁面曲率的影響，產(chǎn)生膨脹波和激波。譬如，對比圖12(a)和13(a)，圖12(b)和13(d)，圖14(b)和15(b)可看到，對同種構(gòu)型，黏性計算與無黏計算同位置的近場壓力信號幅值峰值和峰谷基本無差別。在前半部分機身(對應(yīng)于Xn約20 m之前)的距離，由于激波分布基本一致，因此，近場壓力信號的分布趨勢也基本相同。但在機身變化較劇烈的后半部分，比如在Xn的30～32 m處，由于機身存在較大角度的擴張，黏性計算受邊界層厚度影響使得角度加大，從曲線圖中可看到這部分激波壓縮比無黏更強。

(2)動力短艙對近場信號的影響

對比圖12(a)和14(a)，圖12(b)和14(b)，圖13(a)和15(a)，圖13(d)和15(b)，可以看到，在同種計算條件和同個近場提取位置下，C25P近場壓力信號幅值的峰值要比C25F大，對比流場分布，峰值處位于模型的短艙入口之前，原因在于動力短艙入口唇口的溢流對機身該處來流的減速導(dǎo)致壓力升高，使得該位置處的壓力信號比通氣短艙高。在短艙入口之后的后半部分，近場壓力信號由于C25F與C25P在該區(qū)域激波系的差異而不同。從圖9，10的壓力系數(shù)分布差異可見，由于發(fā)動機噴流的影響，C25P后部的壓力系數(shù)絕對值和壓力系數(shù)變化梯度均比同條件下的C25F大，激波更強烈，因此近場壓力信號在這一部分變化幅值也更大。從圖8對比也可以看到，動力噴流、短艙影響導(dǎo)致下表面的激波更強，因此近場提取的壓力信號也更強。

(3)提取距離和角度對近場信號的影響

對于同種構(gòu)型同種計算條件下，當(dāng)提取的等值面半徑相同時，不同滾轉(zhuǎn)角處的壓力信號是不同的，這是因為滾轉(zhuǎn)角不同，受機體幾何外形的影響而造成激波系形狀也不盡相同，可以從圖11 (a)的示意圖中大致看到差異。從圖12～15可以看到，當(dāng)提取的滾轉(zhuǎn)角相同時，離機體越遠(yuǎn)，捕捉到的近場壓力信號越小，這是黏性耗散和網(wǎng)格數(shù)值耗散導(dǎo)致，符合物理規(guī)律。

2.3 地面聲爆信號

2.3.1 結(jié)果展示

基于所提取的近場壓力信號，本節(jié)選取1倍機身長度和3倍機身長度、滾轉(zhuǎn)角0°處的近場壓力信號作為地面聲爆信號波形參數(shù)法內(nèi)部計算程序的輸入，計算傳播至地面的壓力信號的壓力-時間分布曲線，與參會者的計算范圍和若干參會代表計算結(jié)果進(jìn)行比對分析，如圖16所示。圖例中Outline指代的是官網(wǎng)提供的所有參會者計算值的范圍邊界。表5是根據(jù)計算的地面聲爆信號曲線，轉(zhuǎn)換為地面聲爆的可感知響度PL值，與摘自研討會FTP網(wǎng)站[27]的參會者數(shù)據(jù)的上下范圍比對。為驗證可感知響度轉(zhuǎn)換程序的準(zhǔn)確性，隨機選取官網(wǎng)上部分參會者[27]提供的地面壓力信號值，輸入轉(zhuǎn)換程序，得到的響度值與參會者計算所得值的對比數(shù)據(jù)如表6所示。

(a)C25F invicid，R=1，φ=0°[27-28]

表5 地面聲爆可感知響度PL值

表6 地面聲爆可感知響度程序驗證Table 6 Program verification of perceptible loudness of ground sonic boom

表6中inv和vis后綴1，2，3，4分別指代官網(wǎng)[27]的4份參會者數(shù)據(jù)，數(shù)據(jù)源文件名分別為:“c25d-flowthru-ground-Boeing_Magee-inv-mixed-100” “c25d-flowthru-ground-Park-c25d-flo-visc-tet-200”“c25d-powered-ground-INRIA-c25d-pow-inv-tet-v6-200.b8”“c25d-powered-ground-MorgensternMarconi-vis_2.00Scale_Grid”。

2.3.2 規(guī)律分析

(1)黏性對地面聲爆信號的影響

計算得到的地面聲爆信號分布與參會者的趨勢符合較好，幅值也較為接近。對比圖16(a)和16(b)，圖16(c)和16(d)，圖16(e)和16(f)，圖16(g)和16( h)，可見同構(gòu)型和同種位置壓力信號輸入時，黏性計算的地面聲爆信號幅值較無黏的略大，這與近場壓力信號輸入數(shù)據(jù)的差異有關(guān)，但整體來看，兩者變化趨勢一致。

(2)動力短艙對地面聲爆信號的影響

對比圖16 (a)和16 (e)，圖16 (b)和16 (f)，圖16 (c)和16 (g)，圖16 (d)和16 (h)，可以看到C25F的結(jié)果比C25P得到的幅值峰值較小，受近場壓力信號輸入的影響，地面聲爆信號與近場壓力信號的變化趨勢一致。從表5可以看到，整體上，PL值的差異與地面聲爆信號的差異是一致的，同工況和提取源時，C25P的PL值略大于C25F。

(3)提取距離對地面聲爆信號的影響

對比圖16 (a)和16 (c)，圖16 (b)和16 (d)，圖16 (e)和16 (g)，圖16 (f)和16 (h)，可以看到同個流場，當(dāng)滾轉(zhuǎn)角相同時，不同提取距離的近場壓力信號的輸入對地面聲爆信號的幅值大小影響不大，只造成了曲線后半段的趨勢和形狀不同。

(4)傳播算法的影響

從圖16的對比可見，計算得到的地面聲爆信號曲線的形狀帶有較尖銳的過渡，參會者的曲線過渡則較為圓滑。從表5可以看到，相比參會者的計算結(jié)果，計算得到的地面聲爆可感知響度值偏大。由表6的數(shù)據(jù)對比，經(jīng)過驗證，得到的響度值與參會者計算所得基本一致，差距在1 dB量級，因此排除可感知響度值轉(zhuǎn)換程序引入的誤差，差異原因鎖定為地面壓力信號輸入的差異。由于采用的地面聲爆信號計算程序是基于波形參數(shù)法的內(nèi)部工具[29]，與參會者采用的方法不同，譬如NASA等參會者用的聲爆傳播程序SBoom[30]是求解廣義Burgers方程。本文遠(yuǎn)場至地面聲爆的幅值峰值偏大，因此轉(zhuǎn)換得到的可感知聲爆響度值偏大。后續(xù)可進(jìn)行不同傳播算法的研究。

3 小結(jié)

本文對AIAA聲爆預(yù)測會議的全機帶動力短艙標(biāo)模NASA C25D進(jìn)行了計算與分析，小結(jié)如下:

(1) 由于低聲爆構(gòu)型壓力信號幅值微小(本文計算的NASA C25D近場壓力信號幅值為10-3到10-2量級)，對計算區(qū)域，尤其是內(nèi)域網(wǎng)格加密，有助于減少網(wǎng)格耗散，以捕捉到更精細(xì)的壓力信號。

(2) 黏性計算的發(fā)動機出口尾跡區(qū)域比無黏時更大;短艙入口前的近場壓力信號基本無影響，原本受機身下部分曲率影響而存在的膨脹波或激波波系，在黏性作用下邊界層的存在加大了實際轉(zhuǎn)折角，進(jìn)而影響了對應(yīng)位置的波強度，使黏性計算的近場壓力信號和聲爆信號在對應(yīng)位置更強。

(3) 帶動力短艙構(gòu)型的唇口溢流現(xiàn)象比通氣短艙構(gòu)型更明顯，因此同工況同提取位置時，動力短艙構(gòu)型入口區(qū)域的近場壓力信號的幅值峰值大于通氣短艙構(gòu)型，整體分布趨勢在短艙入口之前基本相似;動力短艙出口區(qū)域的波系比通氣短艙更復(fù)雜多樣，出口之后的區(qū)域受波系影響而呈現(xiàn)較大差異;動力噴流導(dǎo)致的出口區(qū)域近場壓力信號變化幅度比通氣構(gòu)型更大;地面聲爆信號分布趨勢與所輸入的近場壓力信號分布趨勢基本一致。

(4) 同距離下，不同滾轉(zhuǎn)角提取的近場壓力信號和地面聲爆信號因空間波系影響差異較大;同滾轉(zhuǎn)角下，提取位置離機身越遠(yuǎn)，捕捉到的信號越弱。

(5) 不同聲爆傳播計算方法得到的聲爆可感知響度值不同。波形參數(shù)法得到的聲爆曲線變化幅度較求解廣義Burgers方程大。同工況和同近場信號為輸入源時，動力短艙構(gòu)型PL值略大于通氣短艙構(gòu)型。

致謝感謝中國商用飛機有限責(zé)任公司北京民用飛機技術(shù)研究中心、民用飛機設(shè)計數(shù)字仿真技術(shù)北京市重點實驗室提供的平臺及資源的支持。