胡萬林， 于劍， 劉宏康， 閻超

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

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

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

1 數(shù)值方法

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

2 DMD方法

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

(1)

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

(2)

(3)

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

(4)

(5)

(6)

3 計算網(wǎng)格與算例驗證

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

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

圖1 對稱雙圓弧翼型壁面及對稱面網(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 計算結(jié)果與討論

4.1 DMD分析

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

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

圖4 頻率與增長率/衰減率的關(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)隨時間t的變化及其功率譜密度曲線。結(jié)果表明，前4階模態(tài)系數(shù)具有簡諧振動特征，Mode 1幅值隨時間衰減，Mode 2～Mode 4幅值則隨時間增長，與圖4中的模態(tài)的增長率/衰減率一致，同時模態(tài)系數(shù)的幅值隨階數(shù)的升高而減小。對各階模態(tài)系數(shù)進行傅里葉分析，各階模態(tài)包含單一主頻，其對應的St分別為0.141、0.295、0.436、0.577，對應激波自激振蕩St的約1～4倍，且幅值逐漸減小，其中Mode 1在激波抖振過程中占主導作用。

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

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

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

4.2 流場重構(gòu)

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

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

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

圖11 測點壓力隨時間的變化Fig.11 Variation of observation point pressure with time

5 結(jié) 論

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

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

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