韋開君，吳 晴，劉立超

(中國(guó)電子科技集團(tuán)公司第三十八研究所 浮空平臺(tái)部，合肥 230088)

0 引言

圓柱擾流作為流體力學(xué)中的經(jīng)典問題之一，在化工、建筑、水利、海洋等眾多工程領(lǐng)域都有非常重大的研究意義。國(guó)內(nèi)外學(xué)者通過實(shí)驗(yàn)和數(shù)值模擬的方法對(duì)圓柱擾流問題進(jìn)行了大量研究，研究結(jié)果顯示，盡管流體在經(jīng)過圓柱后，會(huì)產(chǎn)生流動(dòng)分離、渦脫落以及尾跡相互作用等一系列復(fù)雜的流動(dòng)現(xiàn)象，但其基本特點(diǎn)是存在較為明確的渦脫落頻率[1-2]。相比于常規(guī)基于時(shí)間歷程的非定常流場(chǎng)分析方法，模態(tài)分解方法可以分辨出流場(chǎng)中存在的大尺度有序運(yùn)動(dòng)，更加直觀地觀察不同特征頻率下的流動(dòng)模式及其發(fā)展變化規(guī)律，特別適合用來描述圓柱擾流的非定常流動(dòng)特征[3]。

2010年Schmid[4]首次將動(dòng)力學(xué)模態(tài)分解(Dynamic Mode Decomposition，DMD)方法引入湍流研究，DMD方法已成功應(yīng)用于分析各種流體力學(xué)問題，例如凹腔流動(dòng)[5]、動(dòng)失速尾流[6]、橫向射流[7]、有限平板流動(dòng)分離[8]、噴管流動(dòng)分離[9]等。這表明DMD方法是探索復(fù)雜非穩(wěn)態(tài)時(shí)空流場(chǎng)本質(zhì)的成熟且有前途的工具。

Tissot等[10]對(duì)比了傳統(tǒng)DMD和殘差最小優(yōu)化DMD方法對(duì)圓柱擾流尾跡PIV測(cè)試結(jié)果分析的影響。張賓等[11]和王智慧等[12]分別應(yīng)用DMD對(duì)旋轉(zhuǎn)圓柱擾流和橢圓柱擾流的PIV測(cè)試結(jié)果進(jìn)行了研究。Stankiewicz等[13]采用DNS直接數(shù)值模擬方法對(duì)二維和三維圓柱擾流進(jìn)行了仿真，并應(yīng)用POD和DMD方法對(duì)仿真流場(chǎng)進(jìn)行分析。Noack等[14]應(yīng)用DNS方法對(duì)旋轉(zhuǎn)圓柱擾流進(jìn)行數(shù)值仿真，對(duì)比了傳統(tǒng)DMD和遞歸DMD對(duì)仿真流場(chǎng)分析的影響。在實(shí)際工程應(yīng)用中，PIV測(cè)試成本過高，而DNS數(shù)值模擬目前尚難以用于真正意義上的工程計(jì)算。因此，需要討論采用工程上常用的湍流模型對(duì)圓柱擾流進(jìn)行數(shù)值仿真，應(yīng)用DMD方法進(jìn)行流場(chǎng)分析的可行性。

本文基于SSTk-ω湍流模型，進(jìn)行了圓柱擾流的非定常仿真，在時(shí)頻分析的基礎(chǔ)上，引入DMD方法對(duì)非定常流場(chǎng)進(jìn)行分析，獲取了各階流場(chǎng)模態(tài)的頻率及分布，提取了引起渦脫落的主導(dǎo)模態(tài)，并通過模態(tài)重構(gòu)來分析圓柱擾流非定常流動(dòng)的特征及變化規(guī)律。

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

對(duì)于一定時(shí)間和空間范圍內(nèi)的流場(chǎng)數(shù)據(jù)，在整個(gè)可用于分析的時(shí)間序列中，流場(chǎng)可表征為一系列離散時(shí)間步的“快照”組合(Snapshots) ，如圖1所示。其中，第i個(gè)時(shí)間步的瞬時(shí)流場(chǎng)分布特性可用列向量vi表示：

圖1 流場(chǎng)動(dòng)力學(xué)模態(tài)分解原理示意

vi=[u1,u2,…,um-1,um]T∈Rm(1)

式中：uj為j(j∈[1,m])節(jié)點(diǎn)上表征流場(chǎng)的物理量(壓力、速度、或密度等)；m為流場(chǎng)的網(wǎng)格節(jié)點(diǎn)數(shù)量。

式中，n代表離散時(shí)間步，則矩陣的維數(shù)為m×n。

對(duì)流場(chǎng)做線性時(shí)不變假設(shè)，則可認(rèn)為vi和vi+1存在線性Koopman映射[15]：

vi+1=Avi(3)

式中，矩陣A∈￡m×m，其特征向量反映了流場(chǎng)內(nèi)各階模態(tài)的空間結(jié)構(gòu)，其特征值則反映了各階模態(tài)對(duì)應(yīng)的頻率。因此，基于全局模態(tài)分解方法求解流場(chǎng)相干結(jié)構(gòu)的關(guān)鍵在于求解矩陣A的特征值和特征向量。一般對(duì)于2000階以下的中小矩陣，采用高斯法、QR法等直接方法即可求解。然而，對(duì)于一個(gè)復(fù)雜流場(chǎng)而言，其網(wǎng)格節(jié)點(diǎn)m通常遠(yuǎn)大于2000，而且矩陣A通常還是稀疏矩陣。對(duì)于此類矩陣，通常采用投影類算法，將高維矩陣投影至低維子空間，通過迭代計(jì)算其特征值和特征向量。

對(duì)于矩陣A和列向量v1，其張成的n維Krylov子空間為：

kn(A,v1)=span{v1,Av1,…,An-1v1}

=span{v1,v2,…,vn} (4)

對(duì)于矩陣A和列向量v2，其張成的n維Krylov子空間為：

kn(A,v2)=span{v2,Av2,…,An-1v2}

=span{v2,v3,…,vn+1} (5)

由式(4)和式(5)可知：

當(dāng)n增加到一定數(shù)量時(shí)，可認(rèn)為子空間基底存在線性相關(guān)，則有：

Vn+1=h1v1+h2v2+…+hnvn+r(7)

寫成矩陣形式則有：

式中:r為殘差向量;e為n階單位向量，H∈Cn×n為上Hessenberg矩陣，具體可寫為：

vi≈Uμi(10)

FDMD=U*AU(11)

式中:Σ為r×r(r≤n) 的對(duì)角陣;U∈￡m×r;W∈￡n×r;星號(hào)代表矩陣的復(fù)共軛轉(zhuǎn)置。將式(6)和(12)式帶入式(11)中，則有：

則可直接對(duì)FDMD∈￡r×r進(jìn)行特征值和特征向量求解：

FDMD=YΛY-1(14)

將POD基底空間上的特征向量Y投影回矩陣A，則流場(chǎng)的各階動(dòng)力學(xué)模態(tài)為：

Φ=UY(15)

與式(3)類似，在低維POD投影子空間內(nèi)，有：

μi+1=FDMDμi(16)

將式(14)、式(15)和式(16)帶入式(10)中，則可基于流場(chǎng)的各階動(dòng)力學(xué)模態(tài)對(duì)瞬時(shí)流場(chǎng)進(jìn)行重構(gòu)：

式中:向量φ為流場(chǎng)的各階動(dòng)力學(xué)模態(tài);α表示各階動(dòng)力學(xué)模態(tài)的貢獻(xiàn)量;λ為矩陣FDMD的各階特征值，則整個(gè)時(shí)間序列的流場(chǎng)重構(gòu)可寫為如下形式：

將式(18)帶入式(12)中，可采用GMRES求解得到diag{α}：

上述過程即為流場(chǎng)的動(dòng)力學(xué)模態(tài)分解方法。

2 算例與分析

2.1 數(shù)值方法

基于二維數(shù)值計(jì)算模擬了圓柱繞流流場(chǎng)，計(jì)算幾何模型如圖2所示。圓柱模型的直徑d=0.04 m，幾何的坐標(biāo)原點(diǎn)為圓柱圓心，坐標(biāo)軸X方向?yàn)閬砹鞣较?，進(jìn)口距離軸心5d，出口距離軸心15d。計(jì)算域上下截面關(guān)于圓柱對(duì)稱，距離軸心5d。

圖2 計(jì)算域幾何模型

數(shù)值計(jì)算采用結(jié)構(gòu)化網(wǎng)格劃分，網(wǎng)格總數(shù)為34.7萬。入口邊界條件為沿X方向速度入口，vx=0.15 m/s。出口邊界條件為流出出口，出口壓力為標(biāo)準(zhǔn)大氣壓。上下表面為滑移壁面，滑移速度的大小和方向和進(jìn)口來流速度相同。圓柱邊界為無滑移壁面。計(jì)算的流動(dòng)介質(zhì)為水，密度ρ=998.2 kg/m3，動(dòng)力粘度系數(shù)μ=0.001003 Pa·s。本文計(jì)算的雷諾數(shù)Re=ρvxd/μ=5971，屬于亞臨界區(qū)。數(shù)值計(jì)算采用SSTk-ω湍流模型，邊界層網(wǎng)格劃分保證y+≤0.5。非定常計(jì)算時(shí)間步長(zhǎng)為0.05 s，時(shí)間離散采用二階隱式，空間離散采用二階迎風(fēng)，壓力速度耦合算法采用SIMPLEC算法。計(jì)算域設(shè)置2個(gè)壓力脈動(dòng)監(jiān)測(cè)點(diǎn)，分別位于距離圓柱后5d的正后方和下方2d處。

2.2 壓力脈動(dòng)時(shí)頻分析

圖3所示為某時(shí)刻卡門渦街分布情況。由圖3可知，流體在給定的進(jìn)口速度流向圓柱后，尾跡發(fā)生周期性渦脫落產(chǎn)生了交替的卡門渦街。

圖3 卡門渦街分布

圖4所示為監(jiān)測(cè)點(diǎn)1的壓力脈動(dòng)時(shí)域分布和頻域分布特征?？芍獔A柱繞流后渦脫落頻率為0.7812 Hz。

圖5所示為監(jiān)測(cè)點(diǎn)2的壓力脈動(dòng)時(shí)域分布和頻域分布特征。監(jiān)測(cè)點(diǎn)2位于圓柱正后方，受到卡門渦街兩側(cè)渦脫落影響，其壓力脈動(dòng)頻率為1.5620 Hz，是渦脫落頻率的2倍頻。

斯特勞哈爾數(shù)St是當(dāng)?shù)貞T性力與遷移慣性力的比值，反映非定常流動(dòng)的周期性[16]。針對(duì)圓柱后尾跡周期性渦脫落呈現(xiàn)的“卡門渦街”現(xiàn)象，可用St表示其周期特性。

式(20)中，f為渦脫落頻率，d為圓柱直徑，V0為來流速度。由前文分析可知，仿真計(jì)算得到的渦脫落頻率f=0.7812 Hz，帶入式計(jì)算得St=0.208。實(shí)驗(yàn)和計(jì)算結(jié)果表明，在 1500 ≤Re≤ 7200 范圍內(nèi)，基準(zhǔn)圓柱的St數(shù)穩(wěn)定在0.2～ 0.21之間。本文計(jì)算的雷諾數(shù)Re=5971，斯特勞哈爾數(shù)St=0.208，符合文獻(xiàn)[17]給出的St數(shù)范圍。

2.3 流場(chǎng)動(dòng)力學(xué)模態(tài)分解

采用DMD方法對(duì)圓柱繞流流場(chǎng)進(jìn)行分析，結(jié)果如圖6所示。需要指出的是，DMD分解的第一階模態(tài)為零頻率模態(tài)，表示非定常流場(chǎng)的一個(gè)周期內(nèi)的平均信息，圖6給出的DMD分解結(jié)果為去除零頻率模態(tài)后的結(jié)果。由圖6可知，流場(chǎng)中存在2種非定常流動(dòng)模態(tài)，1階主模態(tài)頻率為0.7863 Hz，2階模態(tài)頻率為1.5726 Hz。1階模態(tài)是流動(dòng)中最不穩(wěn)定的模態(tài)，全局能量最大，其模態(tài)頻率與前文時(shí)頻分析中監(jiān)測(cè)點(diǎn)1的壓力脈動(dòng)頻率0.7812 Hz基本一致，該模態(tài)流動(dòng)特征在圓柱擾流的非定常流動(dòng)中占主導(dǎo)地位。2階模態(tài)頻率約為1階模態(tài)頻率的2倍頻，與前文時(shí)頻分析中監(jiān)測(cè)點(diǎn)2的壓力脈動(dòng)頻率1.562 Hz基本一致。

圖6 流場(chǎng)標(biāo)準(zhǔn)動(dòng)力學(xué)模態(tài)分解結(jié)果

圖7為經(jīng)DMD分解的1階和2階流場(chǎng)壓力分布特征。由圖7可知，0.7863 Hz的1階模態(tài)表現(xiàn)為流體經(jīng)過圓柱后，在兩側(cè)形成周期性對(duì)稱排列的旋轉(zhuǎn)方向相反的雙列渦線，即卡門渦街。1.5726 Hz的2階模態(tài)主要表現(xiàn)為在圓柱正后方，對(duì)稱排列旋向相反的雙列渦線相互作用形成的類駐波特征。

3 結(jié)論

本文采用非定常數(shù)值仿真方法對(duì)圓柱擾流進(jìn)行了研究，對(duì)流場(chǎng)內(nèi)典型監(jiān)測(cè)點(diǎn)進(jìn)行了時(shí)頻分析，并引入DMD方法對(duì)非定常流場(chǎng)進(jìn)行分析，獲取了各階流場(chǎng)模態(tài)的頻率及分布特征。主要結(jié)論如下：

(1)DMD方法可準(zhǔn)確提取圓柱后方流場(chǎng)的振蕩頻率，與圓柱后方典型監(jiān)測(cè)點(diǎn)的時(shí)頻分析得到的頻率吻合較好；

(2)去除零頻率模態(tài)，DMD提取的1階模態(tài)表現(xiàn)為圓柱后方兩側(cè)形成的雙列渦脫落特征，2階模態(tài)表現(xiàn)為兩側(cè)渦團(tuán)在圓柱正后方相互作用形成的類駐波特征。

(3)DMD方法與CFD非定常計(jì)算相結(jié)合可以作為研究圓柱擾流非定常流動(dòng)特征的有效工具，并推廣至其他非定常流動(dòng)問題的分析應(yīng)用。