并列雙圓柱繞流的動力學模態分解

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

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

0 引 言

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

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

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

1 動力學模態分解方法

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

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

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

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

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

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

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

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

2 雙圓柱繞流流場數值仿真

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

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

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

圖2 圓柱受到的升力阻力隨迭代次數的變化（紅：上方圓柱；藍：下方圓柱）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；曲線波動規律較為相似，均在19 s 左右突然產生了明顯波動；在流動的初始階段，兩個圓柱受到的升力作用相差半個周期，而在19 s 之后，兩個圓柱受到升力作用的時刻保持同步，阻力作用的變化規律與之相反。根據升阻力曲線的變化規律，將雙圓柱繞流的流動過程分為不穩定周期階段和穩定周期階段（圖2（b）中繪制了階段分界線），分階段描述有助于更加清楚地分析圓柱繞流的演化。

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

3 DMD計算分析

3.1 動態模態分析

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

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

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

表1 不穩定周期階段按振幅排列的前5個模態參數Tab.1 Mode parameters of the first five modes arranged by amplitude in unstable period

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

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

圖6 E點對應的模態云圖Fig.6 The modal cloud diagram corresponding to Point E

表2 穩定周期階段按振幅排列的前5個模態參數Tab.2 Mode parameter of the first five modes arranged by amplitude in stable period

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

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

圖8 A點對應的模態云圖Fig.8 The modal cloud diagram corresponding to Point A

3.2 流場重構與預測

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

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

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

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

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

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

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

4 結 論

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

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

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

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