基于POD和DMD的60°交叉管繞流分析

張 嶔，楊青青，唐也婷，王天源

（中國海洋大學(xué)工程學(xué)院，山東青島 266100）

0 引 言

圓柱繞流作為流體力學(xué)中的經(jīng)典問題，在工程應(yīng)用領(lǐng)域十分廣泛，如熱力管道、樁基礎(chǔ)及風(fēng)或水流中的桁架結(jié)構(gòu)等。當(dāng)流體流經(jīng)圓柱表面并從邊界層分離，隨之產(chǎn)生的渦脫落現(xiàn)象是導(dǎo)致尾流場復(fù)雜的主因。而雷諾數(shù)Re是影響圓柱尾流形態(tài)的關(guān)鍵因素。Zdravkovich[1]通過實(shí)驗(yàn)發(fā)現(xiàn)：單圓柱尾渦形態(tài)隨雷諾數(shù)Re改變明顯，圓柱尾渦在40

傳統(tǒng)的渦識(shí)別方法（如Q 準(zhǔn)則、λ2方法[5]等）呈現(xiàn)的流場具有時(shí)空耦合性，尤其對于嚴(yán)重失穩(wěn)的流場，在小尺度渦結(jié)構(gòu)的干擾下難以捕獲主導(dǎo)流場演變的關(guān)鍵流態(tài)。針對傳統(tǒng)渦識(shí)別方法的局限性，諸如本征正交分解（POD）[6-7]和動(dòng)力學(xué)模態(tài)分解（DMD）[8]等降階模型（ROMs），可分別從空間正交性和頻率獨(dú)立性角度降維分解復(fù)雜流場并提取關(guān)鍵模態(tài)信息。

王智慧等[9]利用DMD 提取了橢圓柱尾流中的相干結(jié)構(gòu)，結(jié)果表明大振幅模態(tài)包含了流場的主要信息。孫婉容等[10]利用DMD重建單、雙圓柱流場時(shí)發(fā)現(xiàn)，模態(tài)分解結(jié)果和重建精度取決于數(shù)據(jù)量。袁猛等[11]利用DMD 方法分析并列雙圓柱尾流場的渦量數(shù)據(jù)，發(fā)現(xiàn)增加奇異值的截?cái)嚯A數(shù)并不能對流場預(yù)測效果起到積極作用。Sakai 等[12]利用POD 和DMD 分析了單、串聯(lián)及并排圓柱尾流，發(fā)現(xiàn)間隙比G影響尾流模態(tài)的對稱性。Bai等[13]通過POD對比了波形和光滑圓柱的三維POD模態(tài)特征，發(fā)現(xiàn)僅憑少量低階模態(tài)便可近似重建尾流場。Sirisup 等[14]提取了交錯(cuò)圓柱三維尾流場的主導(dǎo)POD 模態(tài)，并探究了上、下游圓柱尾流形態(tài)對溫度場分布的影響。現(xiàn)階段，POD 和DMD 常用于分析二維圓柱尾流，尚未涉及三維交叉管算例。基于張嶔等[15]通過直接數(shù)值模擬（DNS）獲取的數(shù)值結(jié)果，本文利用POD 和DMD對比分析60°交叉管在Re=200、G=4下的渦量數(shù)據(jù)以探究其尾渦演變規(guī)律。

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

POD[6，7]和DMD[8]可通過降維分解方式識(shí)別時(shí)空耦合尾流場無法呈現(xiàn)的流動(dòng)信息，本章簡要回顧了POD和DMD的基本理論，具體理論可見相應(yīng)引用文獻(xiàn)。

1.1 POD方法

POD 本質(zhì)上是對快照矩陣X的空間協(xié)方差矩陣R=XXT執(zhí)行特征值分解（EVD）以尋求可最佳表示給定流場數(shù)據(jù)的基函數(shù)，等同于在最小二乘意義上對快照矩陣X執(zhí)行奇異值分解（SVD）：

1.2 DMD 方法

Schmid[8]提出的DMD理論假設(shè)相鄰快照之間存在一個(gè)時(shí)間不變線性算子A，在離散時(shí)間的動(dòng)力學(xué)系統(tǒng)中表示為

式中，A∈?n×n，?表示Moore-Penrose偽逆，相當(dāng)于在最小二乘意義下尋求矩陣A的最佳擬合解。

與POD 不同，DMD 假設(shè)流場的模態(tài)和動(dòng)力學(xué)信息分別包含在算子A的特征向量和特征值中。由于矩陣A一般過大，直接對矩陣A執(zhí)行EVD 需大量的內(nèi)存資源，因此通常采用Tu 等[16]提出的DMD 算法以降低計(jì)算成本，最終可得各DMD模態(tài)及其振幅（能量）αi、增長/衰減率gi及頻率fi等信息。

2 算例概述

2.1 模型建立

本文選用三維（3D）Navier-Stokes（N-S）方程，流體運(yùn)動(dòng)粘性系數(shù)ν=0.005 m2/s，雷諾數(shù)Re=UxD/ν=200，流場處于過渡階段[17]。

如圖1 所示，所選研究對象為直徑D=1 m，長度為40D的交叉布置雙管，其中上游圓柱（UC）水平布置，下游圓柱（DC）按一定間隙比G=L/D（L為兩圓柱外徑間隙長度）傾斜60°布置。計(jì)算域長、寬、高分別為（32+G）D、40D、40D。

圖1 60°交叉雙圓柱計(jì)算域示意圖[15]Fig.1 Schematic diagram of two cylinders in 60°crossing arrangement and computational domain[15]

2.2 頻譜及瞬時(shí)流場

在確定計(jì)算模型及網(wǎng)格（圖2）后，基于OpenFOAM 的有限體積法，對交叉雙管（60°，G=4）在雷諾數(shù)Re=200下進(jìn)行直接數(shù)值模擬DNS。計(jì)算時(shí)間步長為dt=0.01 s，在流場穩(wěn)定的第400 s后開始輸出渦量數(shù)據(jù)以作后處理。圖3展示了420 s瞬時(shí)時(shí)刻的渦量等值面及其在z=0平面上的云圖。對20 s（400～420 s）內(nèi)上游UC 及下游DC 圓柱的升力系數(shù)進(jìn)行快速傅里葉變換（FFT）以計(jì)算其頻譜，結(jié)果如圖4 所示。UC 及DC 的頻譜峰值主要集中于0.19 Hz 附近，即上下游圓柱的渦脫落周期約為5.26 s。由于UC 渦脫落對DC 造成影響，導(dǎo)致兩者頻譜峰值略有差異。UC 頻譜相對光滑，但DC 在來自UC 的渦脫落影響下于0.09 Hz、0.28 Hz、0.38 Hz 和0.49 Hz 附近出現(xiàn)頻譜峰值。

圖2 交叉管（60°，G=4）在z=0平面上的網(wǎng)格劃分Fig.2 Computational mesh around two crossing cylinders in 60°,G=4 on the z=0 plane

圖3 基于合渦量展示的60°交叉管附近的流動(dòng)Fig.3 Flow around two crossing cylinders in 60°arrangement based on vorticity magnitude

圖4 交叉圓柱的升力頻譜圖Fig.4 Spectra of lift coefficients on upstream and downstream cylinders

2.3 模態(tài)收斂性分析

數(shù)據(jù)集收斂性涉及快照采樣頻率（采樣時(shí)間間隔的倒數(shù)）和數(shù)量兩個(gè)參數(shù)，完整的數(shù)據(jù)集既要確保低頻（大周期）流動(dòng)現(xiàn)象完整呈現(xiàn)，又要不缺失高頻流動(dòng)現(xiàn)象。據(jù)Nyquist-Shannon采樣準(zhǔn)則，POD和DMD所能捕獲的流體最大頻率fmax為快照采樣頻率fsamp的一半，即

鑒于快照采樣頻率為fsamp=100 Hz，最大捕獲頻率fmax=50 Hz，設(shè)置不同時(shí)間間隔Δt=[0.01，0.02，0.03，0.04，0.05]s（可依次捕獲50 Hz、25 Hz、16.67 Hz、12.5 Hz 和10 Hz 頻率范圍內(nèi)的模態(tài)）和快照數(shù)量m=[100:100:2000]以研究其對模態(tài)收斂性的影響。

由圖4可知，上下游圓柱的升力頻譜峰值約為0.19 Hz，故本文針對不同時(shí)間間隔和快照數(shù)量統(tǒng)計(jì)0.19 Hz 附近DMD 模態(tài)的增長率，頻率誤差保持在10%以內(nèi)則予以考慮。統(tǒng)計(jì)結(jié)果（圖5）表明，Δt=0.01 s、0.02 s和0.03 s分別在m=100～300，100～200和100未捕獲到0.19 Hz附近的模態(tài)，說明數(shù)據(jù)集至少需覆蓋一個(gè)渦脫落周期（5.26 s）。在完整覆蓋一個(gè)渦脫落周期后，模態(tài)增長率趨于收斂并不再受時(shí)間間隔影響。最終，在計(jì)算資源允許范圍內(nèi)，本文選擇時(shí)間間隔Δt=0.01 s 和快照數(shù)量m=2000的數(shù)據(jù)集，該數(shù)據(jù)集可捕獲0～50 Hz范圍內(nèi)及約3.8個(gè)渦脫落周期內(nèi)的2000個(gè)流動(dòng)模態(tài)。

圖5 模態(tài)收斂性分析Fig.5 Analysis of modal convergence

3 模態(tài)分解結(jié)果

3.1 POD統(tǒng)計(jì)分析

圖6展示了前25個(gè)POD 模態(tài)的能量占比E（i式（3））和累積能量占比Eicum（式（4））。其中第1階模態(tài)能量占比最高，占流場總能量的68.57%；第2、3 階和第4、5 階模態(tài)能量成對相同，呈“階梯”狀下降特征；隨著模態(tài)階數(shù)的增加，相鄰模態(tài)能量不再成對相同，且能量占比均低于1%。此外，累積能量占比分布表明，僅需15個(gè)模態(tài)便可捕獲90%的流場總能量。

圖6 POD模態(tài)能量占比及累積能量占比Fig.6 Energy proportion and cumulative energy proportion of POD modes

圖7展示了前100階POD模態(tài)的頻譜分布，其中頻譜峰值作歸一化處理并由黑色像素點(diǎn)代表。統(tǒng)計(jì)結(jié)果表明，模態(tài)頻譜峰值為模態(tài)階數(shù)的線性函數(shù)。綜合圖6和圖7可知，低階POD模態(tài)對應(yīng)高能、低頻流動(dòng)現(xiàn)象；反之，高階POD模態(tài)則對應(yīng)低能、高頻現(xiàn)象。

圖7 POD模態(tài)頻譜分布Fig.7 Spectrum of POD modes

3.2 DMD統(tǒng)計(jì)分析

為避免大振幅但高衰減率的DMD 模態(tài)影響分析，選擇振幅絕對值除以增長/衰減率絕對值|αi|/ |gi|以衡量DMD 模態(tài)對尾流動(dòng)力學(xué)的貢獻(xiàn)程度。圖8（a）展示了排除零頻f=0 Hz 模態(tài)后 |αi|/ |gi|關(guān)于頻率fi的歸一化統(tǒng)計(jì)結(jié)果，圖8（b）進(jìn)一步將 |αi|/ |gi|作為模態(tài)的排序標(biāo)準(zhǔn)對每一DMD 模態(tài)染色以統(tǒng)計(jì)模態(tài)增長率gi關(guān)于頻率fi的統(tǒng)計(jì)結(jié)果，排序越靠前的模態(tài)顏色越深。

圖8 DMD模態(tài)的頻譜Fig.8 Spectra of DMD modes

由圖8（a）可知，f=0.1861 Hz 模態(tài)對流場的貢獻(xiàn)程度最高，該頻率接近上、下游圓柱渦脫落頻率0.19 Hz，其次，f=0.0985 Hz、0.1722 Hz、0.3522 Hz、0.3738 Hz 等模態(tài)也相當(dāng)重要。由圖8（b）可知，模態(tài)的重要程度（顏色）隨頻率增加而減弱（變淺）。對比圖7可知，POD模態(tài)具有多頻性，而DMD模態(tài)則具有單頻特性。

3.3 模態(tài)分析

由上、下游圓柱的升力頻譜（圖4）、POD 頻譜（圖7）和DMD 頻譜（圖8）可知，主導(dǎo)（第2、3 階）POD模態(tài)及主導(dǎo)（f=0.1861 Hz）DMD模態(tài)均強(qiáng)調(diào)f≈0.19 Hz附近的模態(tài)，圖9與圖10展示了第2階POD模態(tài)及f=0.1861 Hz時(shí)DMD 模態(tài)的實(shí)部。由三維等值面（圖9）可知，兩模態(tài)的空間形態(tài)相似，主要表征上、下游圓柱的P形態(tài)渦脫落現(xiàn)象（Zhao[4]及張嶔[15]）。

圖9 POD和DMD模態(tài)的等值面圖Fig.9 Iso-surfaces of POD and DMD modes

圖10 POD和DMD模態(tài)在z=0平面上的云圖Fig.10 Contours of POD and DMD modes on the z=0 plane

圖11（a）進(jìn)一步對比了第2、3階POD模態(tài)和f=0.1861 Hz時(shí)DMD模態(tài)實(shí)部及虛部的時(shí)間系數(shù)。可以看出，第2階POD模態(tài)與DMD模態(tài)實(shí)部、第3階POD模態(tài)與DMD模態(tài)虛部的振蕩頻率相同，兩者間存在一定相位差。圖11（b）分別利用第2、3階POD模態(tài)和f=0.1861 Hz處DMD模態(tài)實(shí)部及虛部的時(shí)間系數(shù)繪制了相位軌跡圖，表現(xiàn)為不同半徑的同心圓，說明時(shí)間系數(shù)為存在90°相位差的正弦曲線，而半徑差異與時(shí)間系數(shù)的不同振幅有關(guān)。POD 軌跡為非規(guī)則圓，主要由時(shí)間系數(shù)的振蕩引起，DMD 軌跡為螺旋圓，螺距表征衰減率gi。

圖11 POD和DMD模態(tài)關(guān)于時(shí)間系數(shù)的對比Fig.11 Comparisons of time coefficients between POD and DMD modes

POD 從空間正交性角度分解模態(tài)，而DMD 則基于頻率獨(dú)立性，因此POD 模態(tài)多以多頻耦合形式出現(xiàn)，而DMD 模態(tài)始終具有單頻特性。POD 可將周期性流動(dòng)現(xiàn)象分解為空間形態(tài)相似且能量相同的一組模態(tài)對，對應(yīng)圖6 中第2、3 階模態(tài)的階梯狀能量分布。該模態(tài)對具有單頻性，且二者時(shí)間系數(shù)間存在90°相位移（Dietmar 和Fasel[18]）。若流動(dòng)現(xiàn)象由于尾流失穩(wěn)而喪失周期性或周期性減弱，則相鄰POD 模態(tài)的空間形態(tài)與能量便不再成對相同（如第6、7階模態(tài)），且其時(shí)間系數(shù)存在波動(dòng)，導(dǎo)致POD 模態(tài)具有多頻性（圖7）。此外，Schmid 等[19]的研究表明，對于周期性流動(dòng)現(xiàn)象，復(fù)DMD 模態(tài)的實(shí)部和虛部可由一組同頻POD模態(tài)ui+iui+1表示，但對于非周期性流動(dòng)現(xiàn)象沒有同等表示。

鑒于POD 模態(tài)和DMD 模態(tài)的空間特征類似，且DMD 模態(tài)的單頻特性便于流體機(jī)理的分析，故下面將展示其他具有典型特征的DMD模態(tài)。

DMD 頻譜（圖8）中存在衰減率及頻率均為0的模態(tài)，由圖12（a）及圖13（a）可知，該模態(tài)主要表征尾流場中的平均渦量分布，對尾流動(dòng)力學(xué)沒有貢獻(xiàn)。隨著頻率的增加，f=0.0985 Hz 和0.1722 Hz 模態(tài)（圖12（b）～（c）和圖13（b）～（c））與下游圓柱的大尺度流向渦有關(guān)，而f=0.3522 Hz和0.3738 Hz模態(tài)（圖12（d）～（e）和圖13（d）～（e））與全局跨向渦脫落有關(guān)。經(jīng)對比，模態(tài)的空間尺度隨頻率的增加而減小。圖14進(jìn)一步展示了上述模態(tài)的時(shí)間系數(shù)，其包絡(luò)線斜率與模態(tài)衰減率（圖14（b））有關(guān)。

圖12 DMD模態(tài)的等值面圖Fig.12 Iso-surfaces of DMD modes

圖13 DMD模態(tài)在z=0平面上的云圖Fig.13 Contours of DMD modes on z=0 plane

圖14 DMD模態(tài)的時(shí)間系數(shù)演變Fig.14 Time coefficients evolution of DMD modes

3.4 降階模型

基于上述分析，利用f=[0，0.0985，0.1722，0.1861，0.3522，0.3738]Hz 等模態(tài)對尾流場進(jìn)行低階重建，結(jié)果如圖15 所示。與f=0.1861 Hz 模態(tài)（圖9（b）和圖10（b））相比，f=[0.0985，0.1722，0.3522，0.3738]Hz 模態(tài)豐富了交叉點(diǎn)后下游圓柱的流向渦演變及遠(yuǎn)離交叉點(diǎn)處連接相鄰平行P 形態(tài)跨向渦的流向二次渦。

圖15 基于主導(dǎo)DMD模態(tài)的重建尾流場Fig.15 Reconstructed wake field based on dominant DMD modes

4 結(jié) 論

本文基于張嶔等[15]通過DNS模擬的60°交叉管在間隙比G=4和雷諾數(shù)Re=200下的渦量數(shù)據(jù)，利用POD和DMD兩種模態(tài)分解方式分析了尾渦的演變規(guī)律，得到以下結(jié)論：

（1）POD 和DMD 模態(tài)的重要性和空間尺度隨頻率的增加而減小，其中大尺度尾渦可由極少數(shù)低頻DMD模態(tài)重建，而高頻模態(tài)則表征小尺度流動(dòng)形態(tài)；

（2）上、下游圓柱的平行P 形態(tài)渦以0.19 Hz 附近的頻率從上、下游圓柱脫落，并按照相同的頻率向下游演變，直至破碎；

（3）下游圓柱上的多個(gè)高頻升力頻譜峰值由來自上游的脫落渦和下游圓柱相互作用引起，并最終導(dǎo)致下游圓柱出現(xiàn)渦激振動(dòng)現(xiàn)象。