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

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

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

0 引 言

圓柱繞流作為流體力學中的經典問題，在工程應用領域十分廣泛，如熱力管道、樁基礎及風或水流中的桁架結構等。當流體流經圓柱表面并從邊界層分離，隨之產生的渦脫落現象是導致尾流場復雜的主因。而雷諾數Re是影響圓柱尾流形態的關鍵因素。Zdravkovich[1]通過實驗發現：單圓柱尾渦形態隨雷諾數Re改變明顯，圓柱尾渦在40

傳統的渦識別方法（如Q 準則、λ2方法[5]等）呈現的流場具有時空耦合性，尤其對于嚴重失穩的流場，在小尺度渦結構的干擾下難以捕獲主導流場演變的關鍵流態。針對傳統渦識別方法的局限性，諸如本征正交分解（POD）[6-7]和動力學模態分解（DMD）[8]等降階模型（ROMs），可分別從空間正交性和頻率獨立性角度降維分解復雜流場并提取關鍵模態信息。

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

1 模態分解方法

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

1.1 POD方法

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

1.2 DMD 方法

Schmid[8]提出的DMD理論假設相鄰快照之間存在一個時間不變線性算子A，在離散時間的動力學系統中表示為

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

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

2 算例概述

2.1 模型建立

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

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

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

2.2 頻譜及瞬時流場

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

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

圖3 基于合渦量展示的60°交叉管附近的流動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 模態收斂性分析

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

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

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

圖5 模態收斂性分析Fig.5 Analysis of modal convergence

3 模態分解結果

3.1 POD統計分析

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

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

圖7展示了前100階POD模態的頻譜分布，其中頻譜峰值作歸一化處理并由黑色像素點代表。統計結果表明，模態頻譜峰值為模態階數的線性函數。綜合圖6和圖7可知，低階POD模態對應高能、低頻流動現象；反之，高階POD模態則對應低能、高頻現象。

圖7 POD模態頻譜分布Fig.7 Spectrum of POD modes

3.2 DMD統計分析

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

圖8 DMD模態的頻譜Fig.8 Spectra of DMD modes

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

3.3 模態分析

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

圖9 POD和DMD模態的等值面圖Fig.9 Iso-surfaces of POD and DMD modes

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

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

圖11 POD和DMD模態關于時間系數的對比Fig.11 Comparisons of time coefficients between POD and DMD modes

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

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

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

圖12 DMD模態的等值面圖Fig.12 Iso-surfaces of DMD modes

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

圖14 DMD模態的時間系數演變Fig.14 Time coefficients evolution of DMD modes

3.4 降階模型

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

圖15 基于主導DMD模態的重建尾流場Fig.15 Reconstructed wake field based on dominant DMD modes

4 結 論

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

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

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

（3）下游圓柱上的多個高頻升力頻譜峰值由來自上游的脫落渦和下游圓柱相互作用引起，并最終導致下游圓柱出現渦激振動現象。