類X-51A飛行器非定常湍流精細(xì)模擬

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

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

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

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

飛行器大迎角分離流動對于數(shù)值計(jì)算方法提出了很高的要求。目前工程應(yīng)用最為廣泛的雷諾平均(Reynolds-Averaged Navier-Stokes,RANS)方程建立于附體邊界層類型流動，對于分離流動預(yù)測精度低，難以滿足精細(xì)模擬的需求[4]。直接數(shù)值模擬(Direct Numerical Simulation,DNS)雖然可以精細(xì)模擬分離流動，但目前還無法應(yīng)用于高雷諾數(shù)的實(shí)際工程領(lǐng)域。近年來，隨著計(jì)算能力的不斷提升，大渦模擬(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)用中得到檢驗(yàn)[6-7]。

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

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

1.1 控制方程

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

(1)

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

1.2 DES類方法

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

1.3 算例驗(yàn)證

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

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

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

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

圖2 瞬時(shí)密度梯度云圖Fig.2 Instantaneous density gradient contour

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

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

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

2.1 模型與網(wǎng)格

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

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

圖5 計(jì)算模型Fig.5 Computational model

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

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

2.2 瞬時(shí)流場

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

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

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

圖9 瞬時(shí)密度梯度截面Fig.9 Instantaneous density gradient slices

2.3 氣動特性

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

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

側(cè)緣卷起的脫體渦還具有明顯的非定常特性，其與舵面的相互干擾會導(dǎo)致舵面氣動力強(qiáng)烈振蕩，可能會帶來機(jī)械共振及疲勞破壞等問題。為了評估舵面力系數(shù)的脈動強(qiáng)度，對其進(jìn)行歸一化處理，顯示脈動量占平均量的百分比。給出歸一化后的整機(jī)和尾部舵面的脈動升力系數(shù)如圖11所示，CLavg為各部件升力系數(shù)的統(tǒng)計(jì)平均值。受到分離渦干擾的“V”型尾舵升力系數(shù)脈動幅值達(dá)到了15%左右，明顯大于整機(jī)和水平尾舵力系數(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)破壞，縮短飛行器使用壽命；同時(shí)脈動壓力會造成嚴(yán)重的氣動噪聲。因此，有必要對于飛行器壁面壓力脈動情況進(jìn)行研究。

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

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

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

圖12 整機(jī)壓力系數(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)測點(diǎn)Fig.14 Monitoring points on “V” shape tail rudder

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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