袁 猛,張新玉,柳貢民,張文平
(哈爾濱工程大學(xué)動(dòng)力與能源工程學(xué)院,哈爾濱150001)
計(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è)效果。
獲取動(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ì)角線上的元素。
建立雙圓柱模型并使用流體動(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
在不穩(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
對(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
奇異值的選擇個(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)還原效果較差。
本文利用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è)。