圓弧翼型跨聲速流動的動態模態分析

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

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

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

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

1 數值方法

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

2 DMD方法

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

(1)

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

(2)

(3)

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

(4)

(5)

(6)

3 計算網格與算例驗證

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

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

圖1 對稱雙圓弧翼型壁面及對稱面網格Fig.1 Symmetric circular-arc airfoil wall and symmetry plane mesh

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

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

4 計算結果與討論

4.1 DMD分析

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

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

圖4 頻率與增長率/衰減率的關系Fig.4 Relationship of frequency with growth rate/decay rate

圖5 DMD的Ritz 值和模態能量與頻率關系Fig.5 Relationship of Ritz value and mode energy with frequency of DMD

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

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

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

圖7 DMD模態的空間分布Fig.7 Spatial distribution of DMD modes

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

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

4.2 流場重構

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

圖9 損失函數隨模態數目的變化Fig.9 Variation of loss function with number of modes

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

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

5 結 論

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

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

3) DMD對于激波間斷處的捕捉能力較弱。采用前7階模態重構流場，能使損失函數降到4%以內，在激波運動區域重構誤差較大。