余華峰， 劉宏康， 陳樹生， 閻超

(北京航空航天大學(xué)航空科學(xué)與工程學(xué)院, 北京 100083)

一直以來高超聲速飛行器因其重要的軍事意義而受到世界各軍事強國的廣泛關(guān)注。進入21世紀以來，以超燃沖壓發(fā)動機為動力的吸氣式高超聲速飛行器逐漸成為了高超聲速飛行器中最為熱門的研究領(lǐng)域之一。美國相繼進行了X-43A[1]和X-51A[2]2款飛行器的飛行演示實驗并取得成功，預(yù)示著高超聲速飛行器技術(shù)在工程實用化方面取得了重大進展[3]。

吸氣式高超聲速飛行器由機體和推進系統(tǒng)高度一體化組成，其構(gòu)型往往較為復(fù)雜。如美國的X-51系列飛行器包含了乘波體構(gòu)型機身，吸氣式?jīng)_壓發(fā)動機和控制舵面等部分。而在飛行器從亞聲速到高超聲速飛行的過程中，存在較大迎角的爬升狀態(tài)，此時飛行器前體壓縮面?zhèn)染?、局部機身突起及控制舵面會誘導(dǎo)產(chǎn)生大范圍的流動分離，形成復(fù)雜的分離渦結(jié)構(gòu)。分離渦結(jié)構(gòu)具有強烈的非定常特性，會對飛行器氣動性能產(chǎn)生明顯影響：一方面分離渦干擾破壞了舵面隨迎角線性變化的氣動特性，嚴重影響飛行器的操控性能；另一方面，分離渦脫落的低頻脈動有可能與飛行器材料的結(jié)構(gòu)固有頻率一致引發(fā)飛行器共振，嚴重時甚至?xí)斐娠w行器的結(jié)構(gòu)損壞。

飛行器大迎角分離流動對于數(shù)值計算方法提出了很高的要求。目前工程應(yīng)用最為廣泛的雷諾平均(Reynolds-Averaged Navier-Stokes,RANS)方程建立于附體邊界層類型流動，對于分離流動預(yù)測精度低，難以滿足精細模擬的需求[4]。直接數(shù)值模擬(Direct Numerical Simulation,DNS)雖然可以精細模擬分離流動，但目前還無法應(yīng)用于高雷諾數(shù)的實際工程領(lǐng)域。近年來，隨著計算能力的不斷提升，大渦模擬(Large-Eddy Simulation,LES)以及RANS/LES混合方法被廣泛應(yīng)用于分離流動的數(shù)值模擬。其中分離渦模擬(Detached-Eddy Simulation,DES)類方法是目前工程領(lǐng)域應(yīng)用最為成熟的RANS/LES混合方法之一。其在近壁區(qū)采用RANS模型，在分離區(qū)則類似于LES方法的Smagorinsky模型[5]，已經(jīng)在諸多復(fù)雜的工程應(yīng)用中得到檢驗[6-7]。

國內(nèi)外針對吸氣式高超聲速飛行器開展了一系列數(shù)值模擬研究。X-51A飛行器設(shè)計過程中曾采用CART3D和OVERFLOW軟件進行了大量數(shù)值計算[2]；鄧艷丹等[8]建立了一個類X-51A模型,對其氣動特性進行了初步研究；羅文莉等[9]研究了一種吸氣式高超聲速飛行器模型在大迎角下的氣動特性。然而，這些研究內(nèi)容大多采用RANS方法，無法獲得非定常特性以及精細湍流。而近年來航空航天領(lǐng)域?qū)Υ擞雨P(guān)注，對復(fù)雜工程外形的湍流數(shù)值模擬提出了更高的要求。因此有必要開展高超聲速飛行器的精細化湍流模擬研究。本文選取類X-51A外形飛行器作為研究對象，針對其在超聲速大迎角狀態(tài)下存在的大范圍非定常分離流動現(xiàn)象，采用延遲DES(DDES)方法開展了精細化湍流數(shù)值模擬研究。分析其在大迎角下的分離流動特點、非定常氣動特性以及壁面壓力脈動規(guī)律。

1 數(shù)值方法與驗證

1.1 控制方程

守恒形式的三維可壓縮Navier-Stokes方程表示為

(1)

式中:Q為守恒變量；F(Fv)、G(Gv)和H(Hv)分別為3個坐標方向的無黏(黏性)通量。一般來說，Navier-Stokes求解時分別獨立離散時間和空間。本文黏性通量(Fv,Gv,Hv)采用中心格式，而無黏通量(F,G,H)則采用考慮物理傳播特性的Roe格式[10]。界面左右變量采用5階WENO插值[11]得到。時間推進則采用LU-SGS格式[12]，為了保證時間精度，選取雙時間迭代進行計算。

1.2 DES類方法

Spalart等[13]在1997年最早提出基于SA(Spalart-Allmaras)方程[14]的DES方法，其主要思想是在湍動能方程的耗散項中引入DES長度尺度LDES=min(dw,CDESΔ)，其中dw為當前網(wǎng)格單 元中心到最近壁面的法向距離，網(wǎng)格尺度Δ=max(Δx,Δy,Δz)，系數(shù)CDES通過算例標定，取為0.65。在壁面附近，流向與展向網(wǎng)格尺度一般遠大于到壁面的距離，所以dw

1.3 算例驗證

為了驗證數(shù)值方法的可靠性，本文選取斜坡空腔算例驗證方法對于超聲速分離再附流動的模擬能力。表1給出了算例的計算狀態(tài)。關(guān)于該算例，前人做了詳細的實驗研究，獲取了豐富可靠的實驗數(shù)據(jù)。Settles等[16]給出了剪切層速度型和壁面壓力系數(shù)分布；Hayakawa等[17]分析了湍流相關(guān)參數(shù)；Shen等[18]則研究了壁面壓力脈動規(guī)律。計算網(wǎng)格如圖1所示，上游和下游網(wǎng)格點數(shù)分別為40×55和560×155，展向長度50.8 mm，均布60網(wǎng)格點，總網(wǎng)格量為534萬。通過調(diào)節(jié)上游長度確?？涨贿吔鐥l件與實驗一致，兩側(cè)采用周期邊條。數(shù)值格式按照之前所述設(shè)置，物理時間步長為0.5 μs，總統(tǒng)計時長為20 ms。

首先給出某一時刻的密度梯度如圖2所示，表明本文采用的模擬方法能夠精細地模擬流動分離再附的結(jié)構(gòu)。圖3給出4個不同站位的剪切層速度型，CFD計算值和文獻[16]實驗值接近，較好模擬了剪切層的發(fā)展過程，其中x和y分別為來流流向距離和垂直來流向上的法向距離，空腔高度D=2.54 cm，U/U∞為當?shù)亓飨蛩俣扰c來流速度的比值。圖4給出的是沿斜坡的壁面壓力時均值以及脈動均方根，P、P∞、Prms分別為當?shù)貢r均壓力、來流壓力和當?shù)貕毫γ}動均方根。對比發(fā)現(xiàn)無論是時均量還是脈動量都與實驗吻合，表明本文計算所得的壁面壓力信息具有一定可信度。因此，本文方法對于超聲速分離流動模擬具有較高的精細度和準確度。

表1 計算狀態(tài)Table 1 Computational condition

圖1 計算網(wǎng)格Fig.1 Computational grid

圖2 瞬時密度梯度云圖Fig.2 Instantaneous density gradient contour

圖3 剪切層速度型Fig.3 Velocity profiles in shear layer

圖4 沿斜坡的時均壓力和壓力脈動均方根Fig.4 Mean and root mean square of pressure fluctuation along ramp

2 計算結(jié)果和分析

2.1 模型與網(wǎng)格

本文首先建立了一種類似于X-51A外形的飛行器模型，如圖5所示。超燃沖壓發(fā)動機安置于乘波體構(gòu)型機體的下側(cè)，飛行器前體壓縮面同時也是發(fā)動機進氣道的組成部分。該模型尺寸與文獻[2]中給出的數(shù)據(jù)接近，能夠反映X-51A這類吸氣式高超聲速飛行器的基本氣動特征。本文主要研究飛行器縱向氣動特性，模型左右對稱且不考慮側(cè)滑，所以采取半模計算。計算采用對接結(jié)構(gòu)網(wǎng)格，對稱面網(wǎng)格如圖6所示。在網(wǎng)格生成時盡量滿足DES類方法的要求，圖7展示了部分壁面網(wǎng)格分布情況。由圖可知在飛行器進氣道入口附近以及飛行器背部等可能存在的分離流動區(qū)域保證網(wǎng)格具有良好的正交性、各向同性及足夠分辨率，以滿足相應(yīng)的物理分辨率。而網(wǎng)格分辨率的選取主要參照斜坡空腔驗證算例。在驗證算例中，關(guān)鍵分離再附區(qū)的單位網(wǎng)格長度取L=0.2 mm，以空腔高度D=25.4 mm作為參考長度，則相對網(wǎng)格尺度Δ=L/H=0.007 87。類比至X-51A飛行器，以尾舵高度H=250 mm作為參考長度，因此單位網(wǎng)格長度保證在1～5 mm的量級。最終，計算網(wǎng)格的總網(wǎng)格量約為3 500萬。

高超聲速飛行器大迎角飛行狀態(tài)主要存在于初始爬升階段以及最終降落狀態(tài)。在爬升階段飛行器馬赫較低，迎角相對較大，此狀態(tài)下流動結(jié)構(gòu)更加復(fù)雜。本文選取飛行器爬升階段的一組狀態(tài)參數(shù)作為計算狀態(tài)，馬赫數(shù)Ma=2.5，飛行高度h=18.5 km，迎角α=10°，單位雷諾數(shù)為5.84× 106/m。計算時的物理時間步長取為2 μs,計算總時長50 ms，確保流動充分發(fā)展且有足夠長時間進行統(tǒng)計平均。

圖5 計算模型Fig.5 Computational model

圖6 飛行器計算網(wǎng)格Fig.6 Computational grids of vehicle

圖7 壁面計算網(wǎng)格Fig.7 Wall surface computational grid

2.2 瞬時流場

圖8首先給出瞬時流場結(jié)構(gòu)圖(Q=200/s2等值面)?？梢杂^測到在10°迎角下，流場中主要的漩渦結(jié)構(gòu)來自于飛行器側(cè)緣的繞流。前體壓縮面附近的氣流繞過側(cè)緣卷起強烈的脫體渦。該脫體渦的漩渦破裂點相當靠前，破裂后仍保持足夠強度一直向后發(fā)展，同時伴隨著眾多小的渦環(huán)結(jié)構(gòu)的生成和耗散。此外還可以看出“V”型尾舵兩側(cè)受到干擾情況明顯不同，相比之下，水平尾舵則基本沒有受到干擾。這將會導(dǎo)致2個舵面的氣動特性明顯不同，具體在后面進一步分析。由于進氣道入口位置的擾動，飛行器底部側(cè)面的流動已經(jīng)表現(xiàn)出強烈的湍流特性，存在大量的壁面湍流結(jié)構(gòu)。

圖9展示了流向3個不同站位的密度梯度云圖，觀察側(cè)緣附近渦結(jié)構(gòu)沿著流向發(fā)展的過程。x=1.5 m截面處于壓縮面根部位置，此時在壓縮面?zhèn)染壐浇嬖谝粋€明顯的分離區(qū)，但是其影響范圍被限制，飛行器背部流動仍然附體；當發(fā)展到飛行器中段x=2.7 m位置時，分離區(qū)逐漸擴大，在正迎角的影響下向飛行器背部發(fā)展，有逐漸演化為脫體渦結(jié)構(gòu)的趨勢；在x=3.9 m位置，渦結(jié)構(gòu)與 “V”型尾舵發(fā)生干擾，進而對舵面的氣動性能產(chǎn)生影響。

圖8 瞬態(tài)流場結(jié)構(gòu)(Q=200/s2)Fig.8 Structure of transient flow field (Q=200/s2)

圖9 瞬時密度梯度截面Fig.9 Instantaneous density gradient slices

2.3 氣動特性

由于分離渦的干擾，導(dǎo)致舵面氣動性能發(fā)生了改變，比如舵面氣動力系數(shù)的線性關(guān)系被破壞。為了驗證這一現(xiàn)象，本文計算了相同狀態(tài)下迎角由1°增至13°過程中飛行器舵面的氣動特性。

圖10給出了2個舵面的升力系數(shù)(CL)、側(cè)向力系數(shù)(Cs)，可以明顯看出當迎角增大至10°以后，“V”型尾舵力系數(shù)存在明顯的非線性特征，而水平尾舵則不存在這一情況。其他力系數(shù)也表現(xiàn)出相同的特性，限于篇幅在此沒有一一給出。

側(cè)緣卷起的脫體渦還具有明顯的非定常特性，其與舵面的相互干擾會導(dǎo)致舵面氣動力強烈振蕩，可能會帶來機械共振及疲勞破壞等問題。為了評估舵面力系數(shù)的脈動強度，對其進行歸一化處理，顯示脈動量占平均量的百分比。給出歸一化后的整機和尾部舵面的脈動升力系數(shù)如圖11所示，CLavg為各部件升力系數(shù)的統(tǒng)計平均值。受到分離渦干擾的“V”型尾舵升力系數(shù)脈動幅值達到了15%左右，明顯大于整機和水平尾舵力系數(shù)的脈動幅值。

圖10 尾舵升、側(cè)力系數(shù)隨迎角變化Fig.10 Variation of lift and side force coefficient of tail rudder with angle of attack

圖11 升力系數(shù)脈動曲線Fig.11 Lift coefficient fluctuation curves

2.4 壓力脈動

吳子牛等[19]指出高超聲速飛行器在20～40 km高度范圍內(nèi)的壓力脈動現(xiàn)象會導(dǎo)致飛行器表面出現(xiàn)局部大載荷，誘導(dǎo)抖振響應(yīng)導(dǎo)致結(jié)構(gòu)破壞，縮短飛行器使用壽命；同時脈動壓力會造成嚴重的氣動噪聲。因此，有必要對于飛行器壁面壓力脈動情況進行研究。

飛行器整體壁面壓力脈動均方根分布如圖12所示，Cprms為整機壓力系數(shù)均方根值。飛行器前體流動附體，因此壓力脈動強度不高。接近進氣道入口位置，壓力脈動幅值發(fā)生突越，由第2.2節(jié)瞬時流場的分析可知在該位置繞過壓縮面?zhèn)染壍牧鲃影l(fā)生了分離。分離區(qū)內(nèi)存在明顯的非定常流動，且分離點和再附點不穩(wěn)定，這些均會導(dǎo)致脈動壓力增大。同時分離反作用于激波導(dǎo)致自激振蕩，造成強烈低頻高幅值脈動壓力。相比之下，飛行器后半部分側(cè)面的壓力脈動相對減弱，與該區(qū)域主導(dǎo)的渦強度弱于壓縮面附近的結(jié)論吻合。再觀察飛行器背部，沿著側(cè)緣分離渦的發(fā)展路徑有較強的壓力脈動，而靠近對稱面附近區(qū)域已經(jīng)遠離分離渦干擾，壓力脈動也相對小許多。

尾部2個舵面的非定常干擾情況相對更加復(fù)雜，圖13更加清楚地展示2個舵面迎風(fēng)和背風(fēng)面的壓力系數(shù)脈動分布情況，其中左側(cè)為背風(fēng)面而右側(cè)為迎風(fēng)面。具體有以下幾個特點：①2個舵面的高壓力脈動位置均處于舵前緣根部位置，主要原因是舵前緣產(chǎn)生脫體激波，與機身壁面附近的分離渦相互干擾，在高馬赫數(shù)飛行狀態(tài)下，激波較為貼體，甚至?xí)c邊界層發(fā)生干擾。激波/邊界層干擾會誘導(dǎo)高強度壓力脈動。②迎風(fēng)側(cè)壓力脈動強于背風(fēng)側(cè)，主要是因為背風(fēng)側(cè)激波與壁面距離較遠，不容易發(fā)生激波干擾。③“V”型尾舵的壓力脈動高于水平尾舵，原因是“V”型尾舵受分離渦干擾更加強烈；可以看到“V”型尾舵背風(fēng)側(cè)高脈動區(qū)沿著舵面向上發(fā)展，對應(yīng)分離渦在該位置的發(fā)展路徑逐漸遠離機身。

飛行器的典型蒙皮結(jié)構(gòu)共振頻段為100～ 500 Hz[20]。當壓力脈動主頻接近這個頻段，其危害將會十分嚴重。高超聲速飛行器飛行過程中由于湍流脈動、分離渦脈動和激波/邊界層干擾等原因，存在不同頻段的壓力脈動。因此需要通過功率譜分析壓力脈動的頻率分布情況。圖14給出了本文在計算模型“V”型尾舵上布置的監(jiān)測點位置，其中括號內(nèi)帶上角標表示處于相同位置背風(fēng)面。另外，為了方便觀察“V”型尾舵附近流動結(jié)構(gòu)，給出監(jiān)測點所在3個截面的密度梯度云圖如圖15所示。由圖15可知，不同截面處脫體渦與 舵前緣激波的干擾位置不同：Slice 1截面脫體渦與舵前緣激波正面干擾，導(dǎo)致舵根部迎、背風(fēng)面的激波不穩(wěn)定；在Slice 2截面迎風(fēng)面激波受到的擾動相對減?。籗lice 3截面前緣激波受渦干擾明顯減小。

圖12 整機壓力系數(shù)脈動均方根云圖Fig.12 Contours of root mean square of pressure coefficient fluctuation for complete vehicle

圖13 尾舵壓力系數(shù)脈動均方根云圖Fig.13 Contours of root mean square of pressure coefficient fluctuation for tail rudders

圖14 “V”型尾舵監(jiān)測點Fig.14 Monitoring points on “V” shape tail rudder

圖15 “V”型尾舵不同高度瞬時密度梯度截面Fig.15 Instantaneous density gradient slices of “V” shape tail rudder at different height

接下來通過圖16給出了2 組監(jiān)測點的壓力功率譜密度(PSD)對比曲線。監(jiān)測點數(shù)據(jù)采樣頻率為250 kHz，總采樣點數(shù)為15 000，采用Burg算法進行功率譜估計。圖16(a)對比了舵前緣處6個監(jiān)測點a(a′)、d(d′)、g(g′)的功率譜，由圖可知,a和d點迎風(fēng)側(cè)脈動能量明顯強于背風(fēng)側(cè)，主要原因是迎風(fēng)面激波較背風(fēng)面更貼近壁面，而激波不穩(wěn)定導(dǎo)致的壓力脈動強于背風(fēng)面渦結(jié)構(gòu)干擾的壓力脈動。相比之下處于舵前緣根部的g和g′，兩點則具有接近的壓力脈動強度，主要是因為舵前緣根部兩側(cè)激波受干擾強度接近。此外值得注意的是，g和g′，兩點壓力脈動存在明顯的主頻Fz=289 Hz，該頻率與蒙皮材料的響應(yīng)頻率接近，可能對結(jié)構(gòu)造成影響。圖16(b)進一步展示了 “V”型尾舵根部沿流向g(g′)、h(h′)、i(i′)6個監(jiān)測點的頻譜。隨著位置后移，壓力脈動強度降低，相應(yīng)的低頻區(qū)域的主頻在200～300 Hz的范圍內(nèi)，強度逐漸減弱。

圖16 不同位置功率譜密度對比Fig.16 Comparison of power spectrum density at different locations

從監(jiān)測點功率譜分析可知，尾舵前緣根部位置存在低頻高強度脈動，針對該現(xiàn)象有必要進一步分析。圖17展示了25.0～28.6 ms之間4個不同時刻的瞬時流場渦結(jié)構(gòu)圖。如圖中標志所示，該段時間內(nèi)，從壓縮面開始產(chǎn)生一個展向渦結(jié)構(gòu)，該結(jié)構(gòu)與沿流向發(fā)展，在“V”型尾舵前緣根部位置大部分被耗散。這樣的發(fā)展周期一直貫穿整個模擬過程，且頻率約為208 Hz，這與舵前緣根部監(jiān)測點的壓力脈動主頻較為吻合，可以認為該流動結(jié)構(gòu)是低頻脈動的主要來源。顯然該周期干擾現(xiàn)象的源頭在前體壓縮面附近。圖18展示了對應(yīng)時刻的壓縮面附近對稱面馬赫數(shù)分布云圖。在該計算狀態(tài)下，壓縮面存在分離區(qū)，而分離區(qū)前緣邊界不穩(wěn)定，導(dǎo)致經(jīng)過壓縮面的流動周期性改變。實際上，如果飛行馬赫數(shù)較低，沒有達到超燃沖壓發(fā)動機起動的要求，進氣道入口往往會存在明顯的分離泡，嚴重時會導(dǎo)致內(nèi)流道產(chǎn)生明顯振蕩。顯然，這樣的非定常流動現(xiàn)象不僅會造成內(nèi)流道結(jié)構(gòu)的破壞，可能也會影響到飛行器外部結(jié)構(gòu)。

圖17 四個不同時刻Q等值面(Q=400/s2)Fig.17 Iso-surface of Q-criterion at four different moments(Q=400/s2)

圖18 四個不同時刻壓縮面附近對稱面馬赫數(shù)云圖Fig.18 Ma contours of symmetry plane near compressed surface at four different moments

3 結(jié) 論

本文采用DDES方法對類X-51A外形高超聲速飛行器大迎角爬升狀態(tài)進行了非定常湍流精細數(shù)值模擬。

1) 在爬升階段，較大的迎角會誘導(dǎo)飛行器側(cè)緣發(fā)生分離以及脫體渦的形成，并且對于尾部舵面等部件產(chǎn)生干擾。

2) 受到分離渦的干擾，尾部“V”型尾舵表現(xiàn)出明顯的非線性氣動特性；同時舵面氣動力脈動幅值明顯增強。

3) 分離渦的存在導(dǎo)致飛行器壁面壓力脈動顯著增強；在舵前緣根部位置存在200～300 Hz的高幅值低頻脈動，該現(xiàn)象由前體壓縮面附近分離區(qū)不穩(wěn)定導(dǎo)致，可能會造成結(jié)構(gòu)破壞。