并列雙圓柱繞流的動力學模態(tài)分解

袁 猛，張新玉，柳貢民，張文平

（哈爾濱工程大學動力與能源工程學院，哈爾濱150001）

0 引 言

計算機性能的提高和快速算法的研發(fā)使得流場的研究和分析變得更加簡便快捷，將復雜的流場結構分解為模態(tài)結構的分析方法越來越多地用在流場分析和穩(wěn)定性判斷等問題中。在獲取流場模態(tài)信息方面，基于流場數(shù)據(jù)的算法比基于控制方程的算法具有更高的效率[1]。傳統(tǒng)的本征正交分解方法可以用于獲取流動結構，但是對于求解可壓縮流動的模態(tài)信息有一定困難。在2009年美國物理學會會議上Schmid[2]首次提出了DMD方法，并以氦氣射流為例進行了驗證，在之后的研究中又進一步說明了該方法在處理仿真和試驗數(shù)據(jù)中的作用。DMD方法在一定程度上彌補了本征正交分解方法的不足。

DMD 方法自從被提出后的十年時間里，已經(jīng)被廣泛地應用于多個領域的分析研究中，例如金融[3-4]、能源[5-6]、人工智能[7]、醫(yī)療[8]等。在流體力學的研究方面，Zhu等[9]利用DMD 方法提取了渦輪增壓器蝸殼在設計工況和溫和喘振條件下的主導的流動結構和對應的頻率，動態(tài)展示了沖擊條件下蝸殼內(nèi)的流動發(fā)展過程。Sarkar 等[10]研究了納米流體通過方柱時混合對流的穩(wěn)定性，發(fā)現(xiàn)銅-水納米流體的高頻模態(tài)比氧化鋁-水納米流體具有更加小的尺度結構。Muld 等[11]利用動模態(tài)分解提取了高速列車周圍的流結構。DMD方法具有數(shù)據(jù)量大、噪聲和誤差影響精度的問題。在優(yōu)化和改善DMD方法的相關研究方面，Brunton 等[12]對DMD 方法進行了改進來消除潛在的動力學和驅(qū)動之間的影響，以便得到更加準確的輸入和輸出數(shù)據(jù)之間的關系。Duke等[13]進行了DMD 方法的誤差分析，提出了減小誤差的方法。國內(nèi)也有一些學者對DMD進行了應用研究[14-15]，但總體來講，關于動力學模態(tài)分解研究的文獻資料還比較少。

本文首先對并列雙圓柱繞流卡門渦街的形成和發(fā)展過程進行非穩(wěn)態(tài)數(shù)值仿真；然后使用DMD 方法對仿真所得的渦量數(shù)據(jù)進行穩(wěn)定性分析，得到各階模態(tài)以及對應模態(tài)的頻率和振幅；再利用模態(tài)數(shù)據(jù)還原了渦量場并預測快照之外時刻的流場，驗證動態(tài)模態(tài)分解方法在預測流場方面的作用；最后對比不同奇異值截斷階數(shù)下的預測效果。

1 動力學模態(tài)分解方法

獲取動力學降階模型可以采取兩種方法，分別是伴隨矩陣法和近似矩陣法，而后者比前者具有更好的魯棒性[16]。關于DMD 的原理推導，可以參見文獻[17]。本文的分析采用近似矩陣法，在此簡要敘述該方法的計算過程。

設在初始時刻系統(tǒng)的觀測值或者實驗值為一個列向量x1，其為初始時刻的快照，列向量的元素中包含了系統(tǒng)不同空間位置的同一個物理量的值，之后每隔相等的時間步長Δt進行一次數(shù)據(jù)采樣，直到第k × Δt時刻停止，得到共計k組數(shù)據(jù)（即k組快照）構成數(shù)據(jù)集合，記為

假設當k增大到一定值時，上述數(shù)據(jù)集合構成的矩陣的秩不再變化，即第k + 1個觀測值可以由前面的k個觀測值線性表示，則可以用矩陣A將兩個相鄰的快照聯(lián)系起來：

對于試驗或仿真過程中的所有k組采樣數(shù)據(jù)，有如下關系：

式中，W 為A'的特征向量構成的矩陣，Λ 為對角矩陣，其對角線元素由A'的特征值構成。使用A'的特征值和特征向量可以得到DMD模態(tài)：

式（8）中Φ 的列向量φi即為DMD 模態(tài)，模態(tài)對應的特征值就是Λ 的對角線上的元素λi，ln(λi)/Δt的實部是模態(tài)的衰減率（放大率），虛部表示模態(tài)出現(xiàn)的頻率。利用DMD 模態(tài)，由式（9）重構或預測t + Δt時刻的數(shù)據(jù)快照[17]：

式中，bi為第i 個模態(tài)的振幅，b 為由元素bi構成的僅包含一列元素的矩陣，b 可以根據(jù)Φb = x1得到。也就是對于第k個數(shù)據(jù)快照，用n階模態(tài)可以表示為

式中，φi為由式（8）得到的矩陣Φ的第i個列向量，λi為由式（7）計算得到的Λ的對角線上的元素。

2 雙圓柱繞流流場數(shù)值仿真

建立雙圓柱模型并使用流體動力學計算軟件Fluent進行仿真。計算域模型網(wǎng)格如圖1所示，模型總高為20 m，總寬度為50 m。兩圓柱直徑均為1 m，圓心縱向間距為2.2 m，圓柱模型的圓心距離左側入口邊界為20 m。

模型入口邊界條件為速度入口，流速為4 m/s。出口設置為壓力出口，工作流體為空氣，對應的物性參數(shù)均按照0 ℃、一個大氣壓力下選取。監(jiān)測量為出口壓力和圓柱面上的升力和阻力。非穩(wěn)態(tài)計算時間步長為0.01 s，計算5 000 個時間步后的升力和阻力隨計算時間的變化，如圖2 所示，達到穩(wěn)定時的渦量分布如圖3 所示。從圖3 可以看出，在當前圓柱間距和參數(shù)設置下圓柱下游產(chǎn)生的渦街方向一致，屬于同步同相模式，與文獻[18]中所述吻合，說明了仿真結果的正確性。

圖1 數(shù)值仿真二維模型圖Fig.1 Numerical simulation of two-dimensional model diagram

圖2 圓柱受到的升力阻力隨迭代次數(shù)的變化（紅：上方圓柱；藍：下方圓柱）Fig.2 The change of lift and drag force of the two cylinders with the number of iteration (Red:upper cylinder;Blue:lower one)

從圖2 中可以看出：兩圓柱的升力、阻力變化周期保持一致，約為1.25 s；曲線波動規(guī)律較為相似，均在19 s 左右突然產(chǎn)生了明顯波動；在流動的初始階段，兩個圓柱受到的升力作用相差半個周期，而在19 s 之后，兩個圓柱受到升力作用的時刻保持同步，阻力作用的變化規(guī)律與之相反。根據(jù)升阻力曲線的變化規(guī)律，將雙圓柱繞流的流動過程分為不穩(wěn)定周期階段和穩(wěn)定周期階段（圖2（b）中繪制了階段分界線），分階段描述有助于更加清楚地分析圓柱繞流的演化。

圖3 50 s時渦量分布圖Fig.3 Vorticity distribution map in 50 s

3 DMD計算分析

3.1 動態(tài)模態(tài)分析

在不穩(wěn)定的周期階段取初始的0～5 s 時間跨度作為DMD 分析的快照數(shù)據(jù)，共計500 個快照，如圖2（a）中快照空間所示。圖4（a）為特征值的實部和虛部在單位圓的分布，位于特征圓外的點共有17個（綠色標識），分別對應8對共軛模態(tài)和一個主導模態(tài)，其均為不穩(wěn)定模態(tài)。其余特征值均位于單位元內(nèi)部或圓上，屬于穩(wěn)定模態(tài)和周期性模態(tài)。值得注意的是，由于計算機舍入誤差的影響，嚴格意義上，到圓心距離完全等于半徑的點并不存在，此問題在穩(wěn)定的周期階段的分析中也會遇到，因此特征圓的方法分析僅為定性分析判斷模態(tài)穩(wěn)定性的判據(jù)。圖4（b）中是在穩(wěn)定的周期階段時間跨度20～25 s 下同樣取500 個快照作為分析數(shù)據(jù)得到的特征值圓，從圖中可以看出幾乎所有的特征值都落在圓內(nèi)和圓上，而與原點的距離小于1.001的特征值個數(shù)為0，即除去舍入誤差的影響以后沒有特征值落在單位圓外部，此時幾乎所有模態(tài)都屬于穩(wěn)定模態(tài)或周期性模態(tài)。圖4（b）與圖4（a）相比，穩(wěn)定性模態(tài)較少，而周期性模態(tài)較多，幾乎所有特征值都落在單位圓上，屬于穩(wěn)定的極限環(huán)階段。

圖4 兩個流動階段的特征值圓Fig.4 The distribution of real and imaginary parts of eigenvalues in the unit element

根據(jù)Xk-11矩陣奇異值的變化，使用硬閾值方法作為判據(jù)截斷奇異值[17]。在不穩(wěn)定的周期階段的計算中丟棄21階之后的奇異值，按照振幅從大到小對前五階模態(tài)從A-E進行編號，模態(tài)參數(shù)如表1所示，振幅和放大率的關系如圖5所示。

表1 不穩(wěn)定周期階段按振幅排列的前5個模態(tài)參數(shù)Tab.1 Mode parameters of the first five modes arranged by amplitude in unstable period

圖5 不穩(wěn)定的周期階段放大率與振幅的關系Fig.5 Magnification in relation to amplitude in unstable period

由圖5 可以看到，A 點的振幅遠大于其他模態(tài)的振幅，即A 點對應的流場能量相對于其他模態(tài)較高，但是其對應的模態(tài)衰減率也最大，因此該模態(tài)是逐漸衰減的模態(tài)，模態(tài)能量也會逐漸降低，其頻率為0，說明其還屬于發(fā)展階段的不穩(wěn)定模態(tài)，該模態(tài)對各個不同時刻的流場結構貢獻相同。B 點對應的模態(tài)能量強度較大，且放大率為-1.27，其頻率為0，其能量強度呈現(xiàn)逐漸減小的趨勢，也屬于發(fā)展中的模態(tài)。從頻率上看，C 和D 兩個點對應的模態(tài)對各個流場的作用是有差異的，D 點對應的模態(tài)對不同時刻的流場影響大于C 點對應的模態(tài)。另外，就模態(tài)本身而言，兩個模態(tài)能量相差不大，但是D 點對應的模態(tài)變化要快于C 點。E 點的振幅相對于其他模態(tài)較大，其衰減率最小且頻率為0，因此可以判斷E 點對應的是平均流場狀態(tài)，即靜態(tài)模態(tài)，如圖6 所示。對穩(wěn)定的周期階段進行DMD 分析，前五階模態(tài)參數(shù)如表2所示。

圖6 E點對應的模態(tài)云圖Fig.6 The modal cloud diagram corresponding to Point E

表2 穩(wěn)定周期階段按振幅排列的前5個模態(tài)參數(shù)Tab.2 Mode parameter of the first five modes arranged by amplitude in stable period

振幅和放大率的關系如圖7 所示。從圖中可以看出，A 點對應的模態(tài)系數(shù)幅值最大，且其放大率為0，頻率極小，模態(tài)系數(shù)幅值基本不變且對各時刻流場的作用基本相同，其對應的是流場的平均狀態(tài)，對應的模態(tài)云圖如圖8所示。B點和C點對應的模態(tài)放大率為正值，幅值較大，且其模態(tài)系數(shù)幅值還有繼續(xù)增大的趨勢，屬于發(fā)散的模態(tài)。D點和E點對應的模態(tài)放大率不大，模態(tài)幅值較大，其反應了流場隨時間的變化特征。對比表1和表2處于兩個階段的模態(tài)放大率，可以看出穩(wěn)定的周期階段放大率比不穩(wěn)定的周期階段小一個數(shù)量級以上，說明不穩(wěn)定的周期階段流場穩(wěn)定性較差。下文將利用本節(jié)提取的模態(tài)信息對原始的渦量場進行重構和預測，以驗證DMD 方法提取流場特征的效果，研究所提取的模態(tài)用于預測渦量場發(fā)展的作用。

圖7 穩(wěn)定的周期階段放大率與振幅的關系Fig.7 Magnification in relation to amplitude in stable period

圖8 A點對應的模態(tài)云圖Fig.8 The modal cloud diagram corresponding to Point A

3.2 流場重構與預測

對圖2（b）中三角標示對應的流場時刻，即4s、8s 和16 s 處進行渦量場的還原及預測，其中8 s 和16 s 時刻的渦量數(shù)據(jù)沒有包含在快照中，故屬于渦量場的預測過程。渦量場云圖還原及預測結果如圖9所示。由圖中可以看出，對于快照空間時刻內(nèi)的渦量場，使用動模態(tài)分解方法已經(jīng)能夠捕捉關鍵的流場特征并進行流場還原，原始流場與還原流場高度吻合。對于“未來時刻”的流場，從圖9（e）中可以看出，用模態(tài)信息進行預測得到的渦量場還是按照原來的發(fā)展規(guī)律向后方不斷變化，但是對于實際渦量已由原來的渦街對稱的形式變化為同步形式，這說明對于流型轉(zhuǎn)換時的流場特征，前面的500個快照顯然沒有捕捉到關鍵的流場變化特征，這是由數(shù)據(jù)快照提供的流場特征不足導致的。圖2（b）中矩形標示的時刻對應的流場處于穩(wěn)定的周期階段，對該三個時刻的渦量場進行還原及預測，結果如圖10所示。從圖中可以看出，盡管只選取了7個奇異值對原始數(shù)據(jù)矩陣進行近似，但是其對快照內(nèi)流場的還原效果以及對流場發(fā)展的預測效果良好，丟失的流場特征較少。

圖9 不穩(wěn)定的周期階段模態(tài)預測渦量場與真實渦量場對比Fig.9 Comparison between the prediction vorticity field and the real field in unstable period

圖10 穩(wěn)定的周期階段模態(tài)預測渦量場與真實渦量場對比Fig.10 Comparison between the prediction vorticity field and the real vorticity field in stable period

3.3 奇異值截斷階數(shù)對流場還原效果的影響

奇異值的選擇個數(shù)會影響流場還原和預測的精度，對不穩(wěn)定的周期階段，奇異值的截斷個數(shù)分別為7、21、50、130 時的渦量值變化進行對比，圖11 為位于上方的圓柱下游方向，距離圓柱圓心1.5 m 處的采樣點的渦量值變化的真實值與還原值對比。

圖11 采樣點在不同奇異值截斷階數(shù)的渦量還原值與真實值對比（紅線對應真實值，藍線對應還原值）Fig.11 Prediction vorticity and the real vorticity at Point M under different truncation orders (Red line: real values,Blue line:prediction data)

從圖11可以看出，截斷個數(shù)為7，即取前7個主要的奇異值時，還原的值與原始值差距較大，DMD方法的還原結果與采樣點渦量真實值的變化周期一致，但是同步性較差。截斷個數(shù)取21 階時，雖然在數(shù)值準確性上稍差，但是兩曲線的變化完全同步，此時通過DMD 方法的還原結果可以預測快照外的近似值以及變化規(guī)律。當截斷個數(shù)為50時，雖然奇異值取的較多，用于還原流場的模態(tài)個數(shù)增加，但是可能是引入了不必要的干擾模態(tài)，導致預測值周期性雖然與原始流場一致，但預測值的變化幅值整體呈發(fā)散趨勢。當取130個奇異值時，預測值與原始值在1 600個快照時刻以內(nèi)都高度吻合。對比圖11 的四個圖可以看出，在第1 600 個快照時刻，渦量還原值與真實值兩條曲線的周期始終不能一致，其與奇異值截斷大小無關。真實的渦量值在第1 600 個快照時刻發(fā)生了劇烈變化，此時為并列雙圓柱繞流流型的轉(zhuǎn)換階段，由于該階段沒有在快照空間內(nèi)，因此單純增加奇異值截斷參數(shù)也未能捕捉到此時刻的流場特征，導致對該時刻的渦量場還原效果較差。

4 結 論

本文利用DMD 方法對仿真得到的雙圓柱繞流流場的渦量數(shù)據(jù)進行了分析，研究了不同流動階段的繞流流場的穩(wěn)定性和奇異值截斷階數(shù)對于DMD預測效果的影響，得出以下主要結論：

（1）DMD方法可以準確提取雙圓柱繞流流動過程中的流場特征和各主導模態(tài)。

（2）對于圓柱繞流穩(wěn)定的周期階段，流場分布較為規(guī)律，流場還原所需模態(tài)較少；對于不穩(wěn)定的周期階段，判斷穩(wěn)定性需要更多的模態(tài)，即對于發(fā)展中的流場，重構和預測流場所需模態(tài)更多，而且效果較差。

（3）增加奇異值的截斷階數(shù)并不能始終對預測效果起到積極作用，需要合理判斷來避免使用對流場特征的貢獻呈現(xiàn)消極作用的模態(tài)進行數(shù)據(jù)還原和預測。