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

（哈爾濱工程大學(xué)動(dòng)力與能源工程學(xué)院，哈爾濱150001）

0 引 言

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

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

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

1 動(dòng)力學(xué)模態(tài)分解方法

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

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

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

對(duì)于試驗(yàn)或仿真過(guò)程中的所有k組采樣數(shù)據(jù)，有如下關(guān)系：

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

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

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

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

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

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

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

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

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

圖3 50 s時(shí)渦量分布圖Fig.3 Vorticity distribution map in 50 s

3 DMD計(jì)算分析

3.1 動(dòng)態(tài)模態(tài)分析

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

圖4 兩個(gè)流動(dòng)階段的特征值圓Fig.4 The distribution of real and imaginary parts of eigenvalues in the unit element

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

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

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

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

圖6 E點(diǎn)對(duì)應(yīng)的模態(tài)云圖Fig.6 The modal cloud diagram corresponding to Point E

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

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

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

圖8 A點(diǎn)對(duì)應(yīng)的模態(tài)云圖Fig.8 The modal cloud diagram corresponding to Point A

3.2 流場(chǎng)重構(gòu)與預(yù)測(cè)

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

圖9 不穩(wěn)定的周期階段模態(tài)預(yù)測(cè)渦量場(chǎng)與真實(shí)渦量場(chǎng)對(duì)比Fig.9 Comparison between the prediction vorticity field and the real field in unstable period

圖10 穩(wěn)定的周期階段模態(tài)預(yù)測(cè)渦量場(chǎng)與真實(shí)渦量場(chǎng)對(duì)比Fig.10 Comparison between the prediction vorticity field and the real vorticity field in stable period

3.3 奇異值截?cái)嚯A數(shù)對(duì)流場(chǎng)還原效果的影響

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

圖11 采樣點(diǎn)在不同奇異值截?cái)嚯A數(shù)的渦量還原值與真實(shí)值對(duì)比（紅線對(duì)應(yīng)真實(shí)值，藍(lán)線對(duì)應(yīng)還原值）Fig.11 Prediction vorticity and the real vorticity at Point M under different truncation orders (Red line: real values,Blue line:prediction data)

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

4 結(jié) 論

本文利用DMD 方法對(duì)仿真得到的雙圓柱繞流流場(chǎng)的渦量數(shù)據(jù)進(jìn)行了分析，研究了不同流動(dòng)階段的繞流流場(chǎng)的穩(wěn)定性和奇異值截?cái)嚯A數(shù)對(duì)于DMD預(yù)測(cè)效果的影響，得出以下主要結(jié)論：

（1）DMD方法可以準(zhǔn)確提取雙圓柱繞流流動(dòng)過(guò)程中的流場(chǎng)特征和各主導(dǎo)模態(tài)。

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

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