圓弧翼型跨聲速流動(dòng)的動(dòng)態(tài)模態(tài)分析

胡萬(wàn)林， 于劍， 劉宏康， 閻超

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

在跨聲速機(jī)翼表面存在著明顯的附面層分離和激波與附面層干擾等復(fù)雜的流動(dòng)現(xiàn)象?？缏曀贍顟B(tài)下激波附面層干擾引起的跨聲速抖振問(wèn)題具有復(fù)雜的非定常、非線性特征，極大地制約了跨聲速飛行器的飛行包線[1]。抖振機(jī)理的研究有利于揭示激波自激振蕩的機(jī)制、激波附面層干擾的機(jī)理以及壓力脈動(dòng)的產(chǎn)生等，從而實(shí)現(xiàn)對(duì)激波抖振的有效控制。1976年，McDevitt等[2]通過(guò)實(shí)驗(yàn)系統(tǒng)研究了激波抖振現(xiàn)象，發(fā)現(xiàn)對(duì)于厚度18%的雙圓弧翼型，在一定的馬赫數(shù)下，激波會(huì)沿著翼型上下表面周期性往復(fù)運(yùn)動(dòng)。Tijdeman和Seebass[3]總結(jié)了以往的跨聲速激波抖振的研究，將其分為A、B和C三類，其中繞厚度18%對(duì)稱雙圓弧翼型的流動(dòng)即為C型，激波在翼型上下面?zhèn)鞑?，向上游運(yùn)動(dòng)的過(guò)程中逐漸減弱，并最終離開(kāi)翼型前緣。2001年，Lee[4]提出了一種自激反饋模型，對(duì)激波抖振的機(jī)理進(jìn)行解釋，認(rèn)為激波的周期運(yùn)動(dòng)誘導(dǎo)的壓力波在附面層的分離區(qū)中向下游傳播，到達(dá)尾緣時(shí)產(chǎn)生的聲波又向上游傳播，為激波的振蕩提供能量，維持系統(tǒng)的穩(wěn)定。Lee的模型得到了許多的研究者[5-10]的驗(yàn)證，對(duì)以后的研究起著重要的指導(dǎo)作用，但是其對(duì)后緣處聲波的產(chǎn)生未給出確切解釋，仍需進(jìn)一步研究。

模態(tài)分析方法能夠?qū)?fù)雜高維流場(chǎng)進(jìn)行低維分解，捕捉對(duì)于非定常流動(dòng)發(fā)展起重要作用的模態(tài)，并在時(shí)間和空間尺度上研究流動(dòng)的發(fā)展。近年來(lái)，基于流場(chǎng)特征提取的模態(tài)分析方法，本征正交分解(POD)[11]和動(dòng)態(tài)模態(tài)分解(DMD)[12]，常被應(yīng)用于復(fù)雜非定常流場(chǎng)的分析中。潘翀等[13]運(yùn)用DMD方法從時(shí)間和空間維度上分析了帶襟翼翼型尾流的流動(dòng)特征，捕獲了尾流場(chǎng)基波和高階諧波的頻率、波長(zhǎng)和傳播速度等豐富的信息。劉宏康等[14]詳細(xì)研究了開(kāi)口空腔的非定常流動(dòng)，結(jié)合POD和DMD方法分析了流場(chǎng)的壓力脈動(dòng)特征，揭示了空腔后緣聲波輻射的2種不同路徑。寇家慶等[15]基于POD和DMD方法，研究了A型跨聲速激波抖振的壓力場(chǎng)，對(duì)比了2種模態(tài)分解方法的特點(diǎn)。陳立為等[10]則研究了雙圓弧翼型激波的周期性運(yùn)動(dòng)，提出了針對(duì)C型激波抖振的自激反饋模型，并用POD分析了主要模態(tài)與激波運(yùn)動(dòng)以及尾緣壓力脈動(dòng)的關(guān)系。但基于系統(tǒng)能量排序的POD方法并不能準(zhǔn)確捕捉與抖振頻率相關(guān)的主要模態(tài)，而基于頻率排序的DMD方法在保留POD方法優(yōu)勢(shì)的同時(shí)，更包含了系統(tǒng)的動(dòng)力學(xué)特性，能提取流場(chǎng)的各階特征頻率。本文針對(duì)厚度18%的對(duì)稱雙圓弧翼型，在時(shí)間維度上對(duì)于跨聲速瞬時(shí)流場(chǎng)運(yùn)用DMD方法，提取包含流場(chǎng)各階主頻的模態(tài)，分析了流場(chǎng)的動(dòng)力學(xué)特性。

1 數(shù)值方法

本文計(jì)算使用基于SA(Spalart-Allmaras)一方程湍流模型[16]的延遲分離渦模擬(Delayed Detached-Eddy Simulation,DDES)[17]求解繞厚度18%對(duì)稱雙圓弧翼型的非定常流動(dòng)，重構(gòu)格式為三階MUSCL(Monotonic Up-stream-centered Scheme for Conservation Laws)[18]，并使用隱式雙時(shí)間步推進(jìn)[19-20]。

2 DMD方法

DMD基于流場(chǎng)的動(dòng)力學(xué)特性，提取的各階模態(tài)在時(shí)間和空間演化特性上相互正交。本文截取流場(chǎng)進(jìn)入穩(wěn)定發(fā)展的非定常狀態(tài)后的500個(gè)樣本進(jìn)行模態(tài)分解，采樣頻率為5 555.6 Hz，后續(xù)討論樣本數(shù)量的關(guān)系，也基于此樣本。

(1)

式中：ui為第i個(gè)瞬時(shí)流場(chǎng)數(shù)據(jù)列向量；N為所取的瞬時(shí)流場(chǎng)樣本總數(shù)目。相鄰瞬時(shí)流場(chǎng)之間的時(shí)間間隔為Δt，假定Δt很小，則兩相鄰瞬時(shí)流場(chǎng)滿足線性映射關(guān)系：

(2)

(3)

式中：W和V均為酉矩陣；Σ為對(duì)角矩陣。再求近似矩陣F，使得A=WFWH，近似矩陣F可通過(guò)求解式(4)得到：

(4)

(5)

(6)

3 計(jì)算網(wǎng)格與算例驗(yàn)證

圖1給出了厚度18%對(duì)稱雙圓弧翼型物理模型及網(wǎng)格，翼型弦長(zhǎng)c=203 mm，展向長(zhǎng)度為0.2c，遠(yuǎn)場(chǎng)為40c。本文使用多塊對(duì)接結(jié)構(gòu)網(wǎng)格，壁面設(shè)置無(wú)滑移絕熱壁，展向?yàn)橹芷谶吔?，遠(yuǎn)場(chǎng)為壓力遠(yuǎn)場(chǎng)邊界。使用兩套網(wǎng)格計(jì)算，grid1、grid2網(wǎng)格分辨率周向、法向和展向分別為465×171×41、625×201×41，第1層網(wǎng)格高度使得y+<1。

本文計(jì)算來(lái)流馬赫數(shù)Ma=0.76，基于弦長(zhǎng)的雷諾數(shù)Rec=1.1×107，迎角為0°，無(wú)側(cè)滑，物理時(shí)間步長(zhǎng)為0.002c/a∞，a∞為來(lái)流聲速。圖2給出了翼型上表面時(shí)均壓力系數(shù)Cp沿流向分布，grid2的計(jì)算結(jié)果與文獻(xiàn)[2]和文獻(xiàn)[10]結(jié)果基本吻合。圖3給出了升力系數(shù)的功率譜密度(PSD)曲線，使用最大熵功率譜估計(jì)方法，在幅值上本文計(jì)算結(jié)果稍大，但變化趨勢(shì)與文獻(xiàn)一致，主頻相近，本文計(jì)算得到的斯特勞哈爾數(shù)St為0.150(f=StU∞/c,f為脈動(dòng)頻率，U∞為來(lái)流速度)，文獻(xiàn)[10]結(jié)果為St=0.148，該主頻對(duì)應(yīng)激波的自激振蕩頻率。

圖1 對(duì)稱雙圓弧翼型壁面及對(duì)稱面網(wǎng)格Fig.1 Symmetric circular-arc airfoil wall and symmetry plane mesh

圖2 壁面壓力系數(shù)分布Fig.2 Wall pressure coefficient distribution

圖3 升力系數(shù)功率譜密度Fig.3 Lift coefficient power spectral density

4 計(jì)算結(jié)果與討論

4.1 DMD分析

對(duì)所取樣本進(jìn)行DMD分析，得到流場(chǎng)的各階主模態(tài)。圖4給出了樣本數(shù)N=100、200、500下的DMD譜，對(duì)應(yīng)的采樣頻率分別為1 388.9、2 777.8、5 555.6 Hz。隨著采樣頻率升高，捕捉流場(chǎng)特征模態(tài)的頻率也越高，且樣本數(shù)200和500的結(jié)果趨近，本文取樣本數(shù)500的結(jié)果進(jìn)行分析。圖4中空心方框?yàn)楸疚乃〉那?階模態(tài)，其中頻率為0的模態(tài)為靜模態(tài)，表征流場(chǎng)的平均特性，其余模態(tài)均成對(duì)出現(xiàn)，其特征值為共軛復(fù)數(shù)。所取模態(tài)增長(zhǎng)率/衰減率均在0附近，說(shuō)明對(duì)流場(chǎng)的非定常發(fā)展過(guò)程具有持續(xù)作用。

圖5給出了DMD的Ritz值以及各階模態(tài)能量與頻率關(guān)系。各階模態(tài)的特征值均處在單位圓附近，個(gè)別模態(tài)位于單位圓內(nèi)，說(shuō)明所取樣本處于準(zhǔn)中性穩(wěn)定狀態(tài)。以往研究者通常根據(jù)模態(tài)幅值或者模態(tài)能量進(jìn)行排列和選取特定的模態(tài)。本文根據(jù)后者選取圖5中前5階模態(tài)進(jìn)行分析，各階模態(tài)主頻分別為激波自激振蕩頻率的倍數(shù)，對(duì)激波的自激振蕩起主要作用。

圖4 頻率與增長(zhǎng)率/衰減率的關(guān)系Fig.4 Relationship of frequency with growth rate/decay rate

圖5 DMD的Ritz 值和模態(tài)能量與頻率關(guān)系Fig.5 Relationship of Ritz value and mode energy with frequency of DMD

圖6給出了DMD前4階模態(tài)(Mode 1～Mode 4)的模態(tài)系數(shù)a(t)隨時(shí)間t的變化及其功率譜密度曲線。結(jié)果表明，前4階模態(tài)系數(shù)具有簡(jiǎn)諧振動(dòng)特征，Mode 1幅值隨時(shí)間衰減，Mode 2～Mode 4幅值則隨時(shí)間增長(zhǎng)，與圖4中的模態(tài)的增長(zhǎng)率/衰減率一致，同時(shí)模態(tài)系數(shù)的幅值隨階數(shù)的升高而減小。對(duì)各階模態(tài)系數(shù)進(jìn)行傅里葉分析，各階模態(tài)包含單一主頻，其對(duì)應(yīng)的St分別為0.141、0.295、0.436、0.577，對(duì)應(yīng)激波自激振蕩St的約1～4倍，且幅值逐漸減小，其中Mode 1在激波抖振過(guò)程中占主導(dǎo)作用。

圖7(a)～(d)分別給出了Mode0～Mode3的壓力脈動(dòng)(實(shí)部)的空間分布，p∞為來(lái)流壓力，文獻(xiàn)[22]研究表明實(shí)部和虛部在流動(dòng)特征上區(qū)別不大,圖中無(wú)量綱壓力值的正負(fù)代表著正負(fù)壓力脈動(dòng)。DMD各階主模態(tài)階數(shù)越高，空間維數(shù)也越高，在翼型上下面，正負(fù)壓力脈動(dòng)相間分布，奇數(shù)階主模態(tài)壓力脈動(dòng)反對(duì)稱分布(上下面壓力脈動(dòng)符號(hào)相反幅值相等)，偶數(shù)階主模態(tài)壓力脈動(dòng)對(duì)稱分布(上下面壓力脈動(dòng)符號(hào)相同幅值相等)。壓力脈動(dòng)的區(qū)間主要集中于翼型的0.5c～0.85c處，即激波抖振的區(qū)間，與文獻(xiàn)[10]結(jié)果相吻合。

圖6 DMD前4階模態(tài)系數(shù)隨時(shí)間的變化及其功率譜密度曲線Fig.6 Variation of coefficient of the first four modes of DMD with time and its power spectral density curves

圖7 DMD模態(tài)的空間分布Fig.7 Spatial distribution of DMD modes

圓弧翼型的激波自激振蕩屬于C型激波抖振，激波在翼型上下面周期性往復(fù)運(yùn)動(dòng)，激發(fā)壓力波的傳播。取t0、t0+1/6T、t0+2/6T、t0+3/6T4個(gè)時(shí)刻觀察激波抖振的過(guò)程，其中t0為激波自激振蕩某個(gè)周期的初始時(shí)刻，T為振蕩周期。圖8給出了激波運(yùn)動(dòng)半周期內(nèi)對(duì)稱面數(shù)值紋影圖和DMD第1階模態(tài)隨時(shí)間的演化過(guò)程，Δρ為密度梯度。觀察圖8紋影圖激波的運(yùn)動(dòng)，對(duì)稱翼型，激波在翼型上下面運(yùn)動(dòng)的過(guò)程是一致的，單獨(dú)分析其在翼型上面的運(yùn)動(dòng),半個(gè)周期內(nèi)，激波從靠近翼型尾緣開(kāi)始向上游運(yùn)動(dòng)至翼型中部，運(yùn)動(dòng)區(qū)間大致為0.5c～0.85c，如圖8紋影圖所示。初始時(shí)刻激波分叉較多，激波較強(qiáng)，如圖8(a)紋影圖所示；向上游運(yùn)動(dòng)過(guò)程中，激波逐漸變?nèi)?，激波后分離區(qū)逐漸增大，至翼型中部時(shí)，激波壓縮成一道，分離區(qū)擴(kuò)大至最大，如圖8(d)紋影圖所示，這一結(jié)果與文獻(xiàn)中描述相似[10]。觀察圖8空間分布，w1、w2分別為負(fù)正壓力脈動(dòng)，隨著激波的抖振，壓力脈動(dòng)也隨之呈現(xiàn)周期性運(yùn)動(dòng)，如圖8紋影圖所示，t0至t0+2/6T時(shí)刻，負(fù)的壓力脈動(dòng)w1向上游移動(dòng)，脈動(dòng)幅值及區(qū)間增大，到t0+3/6T時(shí)刻，w1運(yùn)動(dòng)至翼型中部，脈動(dòng)幅值及區(qū)間有所縮減，同時(shí)壓力脈動(dòng)也將在翼型下面進(jìn)行同樣的周期性運(yùn)動(dòng)。

圖8 翼型對(duì)稱面數(shù)值紋影圖與DMD第3階模態(tài)不同時(shí)刻空間分布Fig.8 Numerical schlieren of symmetry plane of airfoil and spatial distribution of third-order mode of DMD at different moments

4.2 流場(chǎng)重構(gòu)

DMD模態(tài)包含單一的流場(chǎng)特征頻率，少數(shù)的幾階模態(tài)就能刻畫流場(chǎng)的重要信息。圖9給出了損失函數(shù)隨重構(gòu)模態(tài)數(shù)目的變化曲線，Eloss為損失函數(shù)值，n為重構(gòu)模態(tài)數(shù)目，損失函數(shù)定義由文獻(xiàn)[23]給出，可以看出使用前4階模態(tài)的7個(gè)模態(tài)(1個(gè)靜模態(tài)加3對(duì)共軛模態(tài))重構(gòu)就能使損失函數(shù)降到4%以內(nèi)，繼續(xù)增加重構(gòu)的模態(tài)數(shù)目損失函數(shù)變化較小。采用7個(gè)模態(tài)重構(gòu)原始流場(chǎng)，與數(shù)值計(jì)算結(jié)果進(jìn)行對(duì)比，流場(chǎng)重構(gòu)的均方根誤差(相對(duì)于CFD計(jì)算結(jié)果)分布見(jiàn)圖10，RMSE表示均方根誤差值。在激波運(yùn)動(dòng)的區(qū)域以及尾跡區(qū)，重構(gòu)的流場(chǎng)誤差較大，取圖10中A、B、C三處作為測(cè)點(diǎn)，考查重構(gòu)流場(chǎng)對(duì)壓力場(chǎng)的捕捉能力，壓力隨時(shí)間的脈動(dòng)如圖11所示。A點(diǎn)靠近翼型前緣，重構(gòu)誤差較小，重構(gòu)的壓力值與計(jì)算值基本吻合；B、C兩點(diǎn)分別靠近翼型中部和尾緣，處于激波運(yùn)動(dòng)的區(qū)域內(nèi)，重構(gòu)流場(chǎng)整體上刻畫了壓力隨時(shí)間的脈動(dòng)，但在壓力脈動(dòng)的峰值上重構(gòu)誤差較大。對(duì)壓力場(chǎng)的重構(gòu)，表明DMD對(duì)于激波間斷處的捕捉誤差較大。

圖9 損失函數(shù)隨模態(tài)數(shù)目的變化Fig.9 Variation of loss function with number of modes

圖10 流場(chǎng)重構(gòu)的均方根誤差Fig.10 Root mean square errors of flow reconstruction

圖11 測(cè)點(diǎn)壓力隨時(shí)間的變化Fig.11 Variation of observation point pressure with time

5 結(jié) 論

1) DMD能準(zhǔn)確捕捉包含流場(chǎng)特征頻率的各階主模態(tài)。第1階模態(tài)的主頻與激波自激振蕩的頻率相同，在激波的自激振蕩過(guò)程占主導(dǎo)作用，且其正負(fù)壓力脈動(dòng)隨激波的抖振在翼型上下面周期性運(yùn)動(dòng)。

2) DMD提取的前5階模態(tài)的增長(zhǎng)率/衰減率在0附近，對(duì)流場(chǎng)非定常流的發(fā)展起著持續(xù)作用，且其模態(tài)系數(shù)隨時(shí)間變化呈現(xiàn)簡(jiǎn)諧振動(dòng)特征。

3) DMD對(duì)于激波間斷處的捕捉能力較弱。采用前7階模態(tài)重構(gòu)流場(chǎng)，能使損失函數(shù)降到4%以內(nèi)，在激波運(yùn)動(dòng)區(qū)域重構(gòu)誤差較大。