基于本征正交分解和動態模態分解的尾渦激振現象瞬態過程的模態分析

謝海潤， 吳亞東， 歐陽華， 王安正

(上海交通大學 機械與動力工程學院， 上海 200240)

尾渦激振是指處于流場上游物體的尾跡區域內的下游物體受到激勵產生強迫振動的一種氣動彈性現象，會導致被激勵結構的強迫振動.當上游的渦脫落頻率接近被激勵結構的自然頻率時，會產生頻率鎖定現象，振動的幅值會急劇增大，從而危及結構的完整性與疲勞壽命.頻率鎖定現象主要是由于下游物體的振動會在一定程度上改變尾跡區域的流動狀態，從而使激勵力與被激物體之間產生相互作用，使單一的強迫振動現象變成了一種流固耦合的復雜機制.Brika等[1]采用實驗方法研究固定的上游圓柱與彈性的下游圓柱之間的流動耦合現象，指出雙圓柱的同步振動區域比單圓柱的同步振動區域大，同時指出頻率鎖定時同步振動發生在較高的折合速度下.Assi[2]分析了兩個處于上下游的圓柱之間的尾渦激振現象，對比不同雷諾數下4種直徑的圓柱與圓柱相對位置，指出尾渦激振隨著下游圓柱偏離中心線而減弱.劉丹青[3]對振蕩來流下的單圓柱以及串聯雙圓柱進行了數值模擬研究，指出當串聯雙圓柱間距增大時，雙圓柱之間的相互作同更加明顯，特征是在圓柱間隙形成橫流向渦街.

尾渦激振現象中通常包含多個流動特征頻率，具有相對復雜的流動特征.而流場模態分解方法是一種能夠提取流場主要特征、化繁為簡的分析方法.典型的流場模態分析方法有本征正交分解(POD)和動態模態分解(DMD).POD是一種較早提出的模態分解方法，Lumley[4]將其引入到流場分析中，以提取湍流場中的相關性流動結構.POD方法廣泛應用于翼型上的湍流分離流動分析[5]、可壓縮開放腔體流動分析[6]、基于粒子圖像測速法的繞流尾跡動態特性分析[7-8]等.但是POD方法也存在一些不足，為了保證空間模態的正交性，空間模態所對應的時間系數序列往往包含多個頻率成分.此外，POD方法按照能量大小進行排序而沒有按照動態特性的影響進行排序.針對POD方法的弱點，研究者提出了一些改進的算法，比如采用平衡截斷的平衡本征正交方解(BPOD)[9]和結合譜分析方法的譜本征正交分解(SPOD)[10].Schmid[11]最早提出通過動態特性的特征值進行流場模態分析的方法，即DMD方法.Rowley等[12]將Koopman算子與DMD算法相結合，揭示了DMD算法對于非線性系統的意義.DMD方法也存在一些缺點，比如缺少一個單一方法來對特征值的進行排序，即無法確定哪些模態的物理相關性最強，并且DMD模態并非正交模態，導致各個模態之間的內積非零，對于建立降階模型增加了額外的復雜度.DMD方法廣泛應用于橫向射流[12]、離心壓氣機擴壓器的流場分析[13]、波包的穩定性分析[14]等.針對DMD方法的弱點，研究者也提出了一些改進的算法，比如稀疏改進DMD，考慮非正交基截斷誤差的遞歸RDMD方法[15].兩種分解方法的對比研究能夠更好地揭示流動物理現象.Noack等[15]針對圓柱繞流的尾跡，比較了這兩種模態分解方法以及RDMD方法的優缺點.寇家慶等[16]將兩種方法應用于跨聲速抖振的模態分析并進行了比較.

本文針對圓柱與葉片的尾渦激振模型，采用Fluent軟件對翼型從靜止到振蕩過程進行了瞬態數值模擬，并通過POD與DMD方法對振蕩翼型附近的壓力場進行了分解與重構，得到了尾渦激勵現象的主要流動特征，并進一步比較了兩種模態分析方法的特點.

1 數值計算對象與方法

本文數值分析的對象為圓柱與葉片的尾渦激勵模型.由于流動區域中圓柱渦脫落為主要的流動現象，所以采用了文獻[17]中的Realizablek-ε模型作為湍流模型.

計算域尺寸與網格細節如圖1 所示.計算域左側為速度入口，右側為壓力出口，上下壁面、圓柱以及葉型為無滑移壁面.圓柱直徑為50 mm，翼型為弦長36 mm 的NACA6510翼型，圓柱與翼型之間的距離為112 mm.圓柱中心與翼型型心之間的間距約為圓柱直徑的3倍.

流場計算域分為兩個子區域，采用結構化網格劃分.圓柱貼體網格尺寸按照y+=0.173Re0.9Δy/L進行估算.其中：Re為雷諾數；Δy為第1層網格高度；L為邊界層參考尺寸.根據y+=1得到第1層網格高度，建立流場網格.對流場網格的無關性驗證分別針對節點數為 38 656、56 153、68 244及116 842 的4套網格進行.圓柱升力系數的計算結果表明，與網格數為 38 656 的網格相比，其余各套網格的計算結果差別均小于1%，因此滿足網格無關性要求.但是為了更好地體現翼型附近的流場特征，在計算條件允許的情況下，選擇了網格數為 116 842 的網格進行計算.其中計算域2的網格節點數為 45 627，代表了后續模態分解時的空間尺度.

圖1 計算域尺寸和網格細節Fig.1 Computational domain size and mesh details

為了實現對翼型振蕩的模擬，流體計算中采用了動網格策略.計算域2包含翼型，為實現翼型振蕩，將計算域2通過用戶自定義函數整體設置為剛體運動，以前期實驗測量得到的葉片一階自然頻率(26 Hz)進行簡諧振動.該過程中計算域2內部網格相對位置不變， 通過計算域1內網格的變形來實現包含翼型的計算域2的運動.

圖2 瞬態計算結果的監控點與其對應的頻譜Fig.2 Monitoring points of transient calculation and the corresponding spectrum

2 POD與DMD方法

POD和DMD方法具有很多相關性，均可以通過流場各個時間的快照提取主要的流動結構.在流場x∈RN(RN表示實數域)中，對應于時間t的某個變量(速度、壓力及渦量等)的向量場u(x,t)可以用M個以相同采樣頻率提取的流動快照來表示，其中第m個時間點的向量場記作um=u(x,tm)，m=1，2，…，M.

2.1 POD方法

POD方法將流場分解為若干空間正交模態，按照各個模態的能量(特征值)大小進行排序，選擇盡量少的基函數或者模態來捕捉流場中的盡量多的能量，有助于建立流場的降階模型.

(1)

POD方法中流場的脈動值可以表示為所有正交模態與其對應的時間系數的線性疊加：

(2)

式中：aj(t)為第j階模態的時間系數；φj(x)為其對應特征向量.將快照POD的數據以矩陣形式表示：

X=

(3)

POD分析的目的是找出最優的基向量來分解流場數據.特征向量φj(x)可以通過最少的模態數量來表示原始的流場.可以通過快照方法求解φj和其對應的特征值λj：

XTXφj=λjφj,φj∈RM,M

(4)

快照方式的POD特征向量可進一步轉換為原始POD特征向量：

(5)

式中：ψj為POD模態所對應的流場動能.針對流場瞬態計算的特點，即通常網格尺寸遠大于時間步的數量，采用快照方式POD比傳統方式POD求解效率更高.快照方式POD所求解的特征矩陣比傳統方式POD的特征矩陣具有更小的尺度，因此能夠顯著減少計算量.

按照特征值的大小進行排序，假設前r階模態的能量接近流場總體動能：

(6)

則通過前r階模態重構流場：

(7)

其中時間系數反映了POD模態隨時間變化的趨勢，能夠體現POD模態對應的流動結構在重構流場中的能量占比，時間系數可以表示為

aj(t)=〈v(x,t),ψj(x)〉

(8)

2.2 DMD方法

DMD方法基于一個最佳擬合流場動態特性的線性算子，將流場分解到若干具有單一特征頻率和增長/衰減率的模態上.因此DMD可以得到不同頻率的流動結構對于流場的貢獻.

DMD方法需要兩組快照序列作為輸入，快照之間有著恒定的時間差.兩組快照X、X#通過原始時間序列數據以下列方式構建：

(9)

X∈RN×(M-1)

(10)

X#∈RN×(M-1)

假設X和X#之間存在線性關系，即

X#=AX

(11)

而線性算子A∈RN×N可以通過A=X#X+得到，其中X+是矩陣X的偽逆矩陣.DMD分解的模態和特征值定義為矩陣A的特征向量和特征值.通過以下步驟進行計算.

首先，對矩陣X進行奇異值分解：

X=UΣVT

(12)

得到3個矩陣：U∈CN×N，V∈C(M-1)×(M-1)(C表示虛數域)，Σ∈RN×(M-1)，令

(13)

λj=lgμj/Δt

(14)

λj的實部和虛部分別代表了對應DMD模態的增長/衰減率和頻率值.而DMD模態定義為[18]

(15)

用r階DMD模態對流場進行重構，記作

Ψdiag(exp(λjt))b

(16)

式中：bj(0)為各個模態的初始幅值；Ψ為DMD模態的矩陣形式；b為代表模態幅值的向量.將初始時間步的快照v(x,t1)帶入，可以求得各個DMD模態所對應的幅值：

b=Ψ+v(x,t1)

(17)

3 流場的模態分解與結果分析

計算域2共有4 600個快照文件，快照之間的時間步為0.001 s.圖3所示為計算域壓力場某一時刻的瞬態結果以及的計算域2內壓力場的時均結果，圖中p為壓力.可以看出，上游圓柱后已經形成了成對出現的渦街，在計算所處的亞臨界區域內，渦脫落的形態與文獻[19]中相似工況下的結果一致.圓柱尾跡產生的渦團會沿著圓柱中線附近向下游發展，上下交替產生的渦團會持續撞擊在翼型上，從而給翼型施加了一個持續的脈動力.翼型進一步將渦團打散，使其失去主要渦街的一些特征.而時均壓力場則顯示在時間平均的條件下，翼型的上方與下方為一個低壓區域，且吸力面與壓力面的壓力存在差異，體現了翼型對于壓力分布的非對稱性的影響.

圖3 流場的快照與計算域2的時均結果Fig.3 Snapshot of flowfield and time-averaged result of Domain 2

圖4 前50階模態的殘差Fig.4 Residuals of the first 50 order modes

3.1 流場的POD分解

POD分解針對計算域2進行.分解得到的模態根據各階模態的特征值大小進行排序.定義前r階模態所重構的流場的殘差為

(18)

圖4所示為前50階模態的殘差.可以看出，前10階模態的殘差已經低至10-4，因此可以認為前10階模態可以反映流場的主要能量.

圖5所示為以壓力云圖顯示的前10階POD模態.可以看出，模態1和3以翼型為分界，上下呈現相反的壓力.模態2和4則呈現出被翼型一分為二的渦團.模態5、6和模態7、8與低階模態呈現高度的相似性，渦團更加密集.模態9和10 則較為奇異，體現了較復雜的流動形態.

圖6所示為前10階POD模態所對應的時間系數的頻譜，圖中δ為幅值.在針對標準的圓柱繞流渦脫落的流場POD模態分解中，POD模態通常會呈現出以相同幅值和頻率成對出現[15].成對出現的POD模態表示同一個頻率的渦脫落的現象，本質上反映了同一個渦脫落的模式在空間上的發展.由圖6可以看出，模態1和3具有主要的頻率成分，均包含有30 Hz的渦脫落頻率和26 Hz的翼型震蕩頻率，其中翼型震蕩頻率幅值較低.模態2、4，模態5、6和模態7、8可以理解為模態1、3的倍頻.模態9、10的頻譜與前8階模態不同，包含了更多的頻率成分.POD以較少的模態數反映流場的流動特征，對于流場的重構與降解模型的建立具有很大的作用.POD以能量(特征值)進行排序的方式，能夠在一定程度上反映各個模態的流動結構對于流場的貢獻.但該排序方式會忽視各個模態所對應的流動現象，各個模態均為具有多個頻率成分，對于脈動流場的物理意義不夠明確.

圖5 以壓力云圖顯示的前10階POD模態Fig.5 The first 10 POD modes displayed by pressure contour

圖6 前10階POD模態時間系數的頻譜Fig.6 Spectrum of time coefficient of the first 10 POD modes

3.2 流場的DMD分解

DMD分解的零階零頻模態代表了平均流場.圖7所示為DMD模態的特征值分布以及前10階模態所對應的頻率與幅值，圖中下標real，imag分別表示實部及虛部.可以看出，幾乎所有的特征值都分布于單位圓上，說明流場中主要流動結構均處于穩定的狀態，與流場設置的邊界條件符合.其中較大的點是按照2.2節中模態幅值定義進行排序的前10階模態.將前10階模態所對應的頻率與歸一化的幅值以柱狀圖表示，可見DMD第1階模態的幅值與平均流場的幅值較為接近.第1階模態對應渦脫落頻率30 Hz，而模態2，3，4，7對應1階模態的倍頻，但模態幅值隨模態數增加而迅速下降.第9階模態對應翼型振蕩頻率26 Hz，其模態幅值相對較低.

DMD模態以共軛形式成對出現，每一對共軛的特征值可以視為一個DMD模態.圖8所示為前10階的DMD模態.可以看出，DMD的第1階模態和POD的第1階模態非常相似，而其對應的頻率為渦脫落頻率30 Hz.而模態2、3、4以及7則可視為模態1的倍頻.代表翼型振蕩頻率26 Hz的第9階模態與第1階模態在壓力分布的形態上具有高度相似性.

圖7 DMD模態的特征值分布以及前10階模態所對應的頻率與幅值Fig.7 Distribution of eigenvalues and corresponding frequencies and amplitudes of the first 10 DMD modes

圖8 以壓力云圖顯示的前10階DMD模態Fig.8 The first 10 DMD modes displayed by pressure contour

4 POD與DMD對比

為了將POD和DMD的結果進行對比，將反映流場主要特征頻率的模態進行比較.POD是將流場按照空間和時間進行分解，得到POD模態與其對應的時間系數.各個POD模態的幅值反應了瞬時流場的結構，時間系數則表示對應模態隨著時間變化的情況.而DMD是將流場按照空間和頻率進行分解，得到的DMD模態都有著對應的增長/衰減率和頻率，進一步可以建立每個模態所對應的時間系數.圖9所示為第1、2和9階POD模態以及DMD模態的時間系數.參考圖6中POD模態的時間系數頻譜，第1階模態中主要包含30 Hz和26 Hz的頻率成分，其分別對應于圓柱的渦脫落頻率和翼型的振蕩頻率，而渦脫落頻率幅值相對較高.第2階模態中主要頻率成分為第1階模態的倍頻.第9階模態的頻譜中具有較多的頻率尖峰，反映其包含相對復雜的頻率成分，但是其中30 Hz處無尖峰，說明該模態與渦脫落現象的相關性較低.從時間系數可以看出：第1、2階POD模態的時間系數在t<1 s時保持30 Hz穩定振蕩；t>1 s時，隨著26 Hz的翼型振蕩頻率的進入，形成了較為明顯的節拍.第9階POD模態的時間系數更明顯地反映新的頻率成分進入穩定流場的過程，1 s前后的振蕩形態有著明顯的差異.DMD模態對應著單一的頻率成分，第1階DMD模態對應30 Hz的渦脫落頻率，第3階模態對應60 Hz的渦脫落頻率的二倍頻，第9階模態對應了26 Hz的翼型震蕩頻率.可以從時間系數中看出，第1階和第2階模態的幅值在整個計算周期中基本保持穩定，而第9階模態的幅值則逐漸增大，表明翼型振蕩對流場的影響是從無到有，從小到大的.對比POD模態和DMD模態可以看出，POD模態的時間系數明顯反映出翼型振蕩前后流場所發生的變化；而DMD模態的時間系數中則沒有明顯的分界線.但DMD模態反映出各個具有物理含義的頻率成分的變化規律，有利于針對性地研究流場中的特定流動結構.

基于前文得到的前10階POD模態與DMD模態對流場進行重構，基本能較好地反映流場特征.但進一步的誤差分析顯示，POD模態和DMD模態重構的流場具有差異性.圖10所示為用前10階POD和DMD模態重構流場壓力的方均根誤差(RMSE)的瞬時結果.可以看出，POD模態重構的流場誤差基本處于一個穩定水平，但是圖中t<1 s時的脈動幅值小于t>1 s的脈動.而DMD模態重構的流場誤差起初較小，然后隨著時間增加逐漸增大，t>1 s時流場的誤差增加到一個較大的值并以較大的振幅振蕩.比較兩者的差異，可以看出，t<0.5 s時，DMD誤差小于POD誤差，但是t繼續增大時，POD誤差小于DMD誤差.

圖10 用前10階POD和DMD模態重構流場的方均根誤差的瞬時結果Fig.10 RMSE of reconstructed flow field with the first 10 POD and DMD modes

圖11 前10階POD和DMD模態重構流場方均根誤差的時均云圖Fig.11 Time-averaged RMSE contour of reconstructed flowfield with the first 10 POD and DMD modes

圖11所示為用前10階POD和DMD模態重構流場的RMSE歸一化后的時均云圖.兩者采用相同的標尺表示.可以看出，時均結果中POD的誤差最大值顯著小于DMD誤差.而兩者誤差較大的區域均處于翼型前緣上下的部分.這是由于該區域流動現象較為復雜，存在著多種其他頻率成分的流動結構.DMD模態分解的誤差主要分布區域與DMD模態1、5、9誤差區域較為重合，說明誤差可能主要來源于相關的頻率成分.

通過對比可以看出，POD模態分解在處理截斷誤差時比DMD模態的整體表現更為優秀.主要原因為POD模態分解時以正交形式構建的基向量，因此較少的模態數即可較好地表示整體流場的脈動.盡管DMD的最大誤差大于POD分解的誤差，但在流場穩定脈動的過程中，DMD分解誤差較小.

5 結論

本文通過Fluent軟件針對翼型由靜止到振蕩過程的瞬態數值模擬，提取了振蕩翼型周圍的流場，通過POD和DMD兩種模態分解方法，對流場進行了分解與重構，捕捉到尾渦激振現象的主要流動特征，得到上游的渦脫落頻率與翼型的振蕩頻率以及其對應的模態特征.基于尾渦脫落現象與流場瞬態變化的特點，對兩種模態分解方法進行了對比，主要結論有：

(1) POD模態分解是將流場按照空間與時間的方式進行分解，其中各個空間模態是正交的.POD模態按照能量(特征值)進行排序，較少的POD模態數即可將流場的殘差快速降低.POD模態的時間系數可能包含多個頻率成分，不利于研究流場中的主要脈動結構，但可以較為清晰地在時間上顯示變化的關鍵節點.

(2) DMD模態分解是將流場按照空間和頻率的方式進行分解，其中各個空間模態并非正交.DMD模態的排序不具有單一的正確方式，因此在物理上客觀地確定模態的最相關性比較困難.DMD模態對應著單一的頻率成分，可以提取特定頻率的動態流動結構.模態的增長/衰減率可以反映頻率成分的變化情況，但難以表示關鍵的時間節點.

(3) 在本文的算例中，POD模態分解的最大誤差小于DMD模態分解的誤差，但DMD分解在穩定脈動的流場中可以得到更小的誤差.

POD和DMD兩種模態分解方法具有各自的優點和缺點.在流場穩定脈動時可以采取DMD方法，而伴隨其他隨時間變化的流動現象時，需要結合POD方法或其他方法.