郭志遠(yuǎn)，虞培祥,b，歐陽華,b

(上海交通大學(xué) a.機(jī)械與動力工程學(xué)院；b.燃?xì)廨啓C(jī)與民用航空發(fā)動機(jī)教育部工程研究中心,上海 200240)

開爾文-亥姆霍茲不穩(wěn)定性(Kelvin-Helmholtz Instability)是指在一條速度不連續(xù)的切變線上產(chǎn)生渦度集中而導(dǎo)致的流動不穩(wěn)定性，簡稱K-H不穩(wěn)定性.在自然界中存在諸多K-H不穩(wěn)定性，比如高空的波狀云、風(fēng)引起的水面波和云團(tuán)鋒面的不穩(wěn)定現(xiàn)象等.在實(shí)際工程中，有學(xué)者發(fā)現(xiàn)后向臺階流[1]、翼型表面的分離流動[2]和軸流壓氣機(jī)旋轉(zhuǎn)不穩(wěn)定性[3]等現(xiàn)象都與K-H不穩(wěn)定性存在一定的關(guān)聯(lián).對于經(jīng)典的圓柱繞流問題,當(dāng)雷諾數(shù)Re<49時(shí)，流動是定常的，在圓柱下游存在一對穩(wěn)定的回流區(qū)；當(dāng)49194時(shí)，流動的三維效應(yīng)開始出現(xiàn).文獻(xiàn)[4]指出，當(dāng)Re達(dá)到 1 200 以上時(shí)，圓柱下游尾跡的剪切層中就會出現(xiàn)K-H不穩(wěn)定性現(xiàn)象.

本征正交分解(POD)是最早提出的一種模態(tài)分解方法，POD能將高階的、非線性的流場分解為若干個(gè)空間正交模態(tài)，按照模態(tài)的能量(即特征值)大小進(jìn)行排序，進(jìn)而選出流動的主要模態(tài).動態(tài)模態(tài)分解(DMD)則是近年來新興的一種數(shù)據(jù)處理方法，按照頻率對系統(tǒng)進(jìn)行排序，提取系統(tǒng)的特征頻率，從而得出不同頻率的模態(tài)對流場的貢獻(xiàn)，進(jìn)而分析復(fù)雜流場的非定常脈動的主要特征或者建立低階的流場動力模型.POD將系統(tǒng)分解為低階的空間正交模態(tài)，但這些模態(tài)會混疊多種頻率成分，因此不利于流場的物理機(jī)理分析，而DMD方法可以得到模態(tài)的單一頻率和擴(kuò)散率(或耗散率)，則有利于流場的結(jié)構(gòu)分析和動力學(xué)分析[9-10].文獻(xiàn)[11-16]使用DMD方法對50≤Re≤13 000圓柱繞流的數(shù)值與實(shí)驗(yàn)結(jié)果進(jìn)行了模態(tài)分解，進(jìn)而得到流場不同頻率的模態(tài)特征，達(dá)到了分解流場和捕捉流場渦結(jié)構(gòu)特征的研究目的.這些文獻(xiàn)均準(zhǔn)確地捕捉到卡門渦街的頻率特征和模態(tài)，部分學(xué)者還對流場中的高頻和低頻特征開展了研究.因此，DMD方法是一種有效的數(shù)據(jù)處理方法，可以用來分析流場中的頻率特征和模態(tài)特征.不過需要指出的是，這些學(xué)者并未關(guān)注剪切層K-H不穩(wěn)定性.

本文將對Re=2 000,3 900,5 000的圓柱繞流開展大渦模擬，利用大渦模擬得到的速度場，通過傳統(tǒng)的選擇合適監(jiān)測點(diǎn)的頻譜分析和DMD分析的方法來探究圓柱繞流剪切層K-H不穩(wěn)定性的特征，并對比分析兩種數(shù)據(jù)處理方法的特點(diǎn).

1 研究方法

1.1 數(shù)值方法與驗(yàn)證

本文考慮不可壓縮流動，為實(shí)現(xiàn)高雷諾數(shù)的流動模擬，采用大渦模擬方法進(jìn)行數(shù)值計(jì)算，因此，以濾波后的不可壓縮Navier-Stokes方程作為控制方程.采用張量和愛因斯坦求和標(biāo)記法，控制方程可表述為

(1)

(2)

本文數(shù)值方法基于有限體積法，采用PISO(Pressure Implicit Split Operator)算法，時(shí)間和空間都采用二階精度格式.

圖1是圓柱繞流x-y平面的網(wǎng)格示意圖.本文除角度以外，所有變量均采用無量綱單位(包括長度、速度等).其中，圓柱直徑定義為D=1，流體域展向長度為πD，x方向長度為40D，y方向在圓柱兩側(cè)各有15D的距離，網(wǎng)格節(jié)點(diǎn)總數(shù)目是275萬.算例的進(jìn)出口邊界條件分別是速度入口和自由出口，來流速度設(shè)為U0=1，流體域兩側(cè)和展向的邊界條件均是對稱邊界條件，圓柱壁面是無滑移壁面.

圖1 圓柱繞流x-y平面的網(wǎng)格示意圖Fig.1 Computational meshes for flow around a cylinder

本文先對Re=3 900的圓柱繞流問題開展數(shù)值模擬.對于該問題，目前不僅有大量LES的結(jié)果，而且也存在DNS的“標(biāo)準(zhǔn)解”，可用來檢驗(yàn)不同LES方法的準(zhǔn)確性.表1是斯特勞哈爾數(shù)(St=fD/U)、阻力系數(shù)(Cd)、升力系數(shù)均方根(Cl,rms)和圓柱表面背壓系數(shù)(Cp,b)的比較結(jié)果.圖2是圓柱下游軸線上(坐標(biāo)為y=0，z=0.5πD)x方向的速度ux的時(shí)均分布.圖3是圓柱壁面不同角度(θ)處的壓力系數(shù)(Cp)分布圖.圖4和5分別是速度ux在不同位置處的時(shí)均值和脈動量均方根分布圖.大渦模擬直接計(jì)算大尺度的渦旋而對小尺度渦旋使用亞格子應(yīng)力的方式處理，因此在速度脈動值的統(tǒng)計(jì)上與DNS的結(jié)果略有差異，但根據(jù)表1和圖2～5可以看出，本文大渦模擬的各個(gè)特征參數(shù)與已有研究結(jié)果非常接近，說明本文的LES結(jié)果是準(zhǔn)確可信的.

表1 Re=3 900的圓柱繞流特征結(jié)果對比Tab.1 Typical flow results at Re=3 900

圖2 Re=3 900圓柱繞流的流場軸線處ux沿流向分布的時(shí)均圖Fig.2 Distribution of mean streamwise velocity ux along wake centerline at Re=3 900

圖3 Re=3 900圓柱繞流的Cp沿圓柱表面周向分布圖Fig.3 Pressure coefficient on cylinder surface at different angles at Re=3 900

圖4 Re=3 900的尾跡不同位置的時(shí)均流向速度剖面圖Fig.4 Mean streamwise velocity at three locations in the wake of circular cylinder at Re=3 900

圖5 Re=3 900的不同位置的流向速度脈動值均方根分布圖Fig.5 Variance of streamwise velocity fluctuations at three locations in the wake of circular cylinder at Re=3 900

1.2 數(shù)據(jù)后處理DMD方法與驗(yàn)證

S=[s(t1)s(t2) …s(tm)]∈Rn×m

(3)

S#=[s(t2)s(t3)…s(tm+1)]∈Rn×m

(4)

圖6 Re=100圓柱繞流特征值分布圖Fig.6 Distribution of eigenvalues of DMD at Re=100

本文選擇Re=100的二維圓柱繞流的算例來驗(yàn)證DMD算法.Re=100圓柱繞流的St=0.165，選取2個(gè)卡門渦街周期的流場速度快照作為DMD的輸入，每個(gè)卡門渦街周期內(nèi)包含60個(gè)等時(shí)間間隔的取樣.圖6是DMD各階模態(tài)特征值的實(shí)部和虛部所構(gòu)成的點(diǎn)在單位圓附近的分布，其中橫軸為特征值的實(shí)部，縱軸是特征值的虛部.在DMD特征值的分析中，當(dāng)特征值對應(yīng)的點(diǎn)位于單位圓內(nèi)時(shí)，表示該模態(tài)是穩(wěn)定模態(tài)；當(dāng)特征值對應(yīng)的點(diǎn)位于單位圓上時(shí)，表示該模態(tài)是周期性模態(tài)；當(dāng)特征值對應(yīng)的點(diǎn)在單位圓外，則表示該模態(tài)是不穩(wěn)定(發(fā)散)模態(tài).由圖6看出，各階模態(tài)特征值都位于單位圓上或者位于單位圓內(nèi)同時(shí)又非?？拷鼏挝粓A，因此所有的模態(tài)都是周期性模態(tài)，此處的圓柱繞流處于穩(wěn)定極限環(huán)狀態(tài).按能量‖Φ‖的大小排序，如圖7所示.圖中前5個(gè)頻率(f)依次對應(yīng)0,St0,2St0,3St0,4St0，其中St0=0.165.由此看出DMD方法可準(zhǔn)確捕捉到圓柱繞流卡門渦街的主頻及其各階倍頻.

根據(jù)特征值能量的大小提取DMD前5階的模態(tài)與文獻(xiàn)[24]的Re=100的DMD前5階模態(tài)作對比，結(jié)果如圖8所示.圖8左側(cè)由上到下依次是文獻(xiàn)[24]前5階模態(tài)(M1、M2、M3、M4、M5)，圖8右側(cè)由上到下依次是本文得出的前5階模態(tài).由圖8容易看出，兩組DMD結(jié)果均捕捉到了圓柱下游尾跡處復(fù)雜的渦結(jié)構(gòu)，每一階模態(tài)的流場特征都是對稱分布的，同時(shí)在圓柱下游的遠(yuǎn)處也都捕捉到了渦系的分布特征.兩組DMD的結(jié)果基本一致，也證明了本文DMD算法的可靠性.

圖7 Re=100圓柱繞流DMD能量譜Fig.7 DMD energy spectra at Re=100

圖8 Re=100圓柱繞流的前5階DMD模態(tài)對比圖Fig.8 Comparison of first five dominant DMD modes for flow around a cylinder at Re=100

2 結(jié)果與分析

2.1 K-H不穩(wěn)定性的特征頻率

2.1.1傳統(tǒng)的頻譜分析結(jié)果 剪切層的速度梯度較大，是一個(gè)位置狹小的窄帶，fSL對監(jiān)控點(diǎn)的分布位置比較敏感，因此本文布置大量的監(jiān)測點(diǎn)來捕捉fSL.通過比較不同位置處的速度的能量譜密度α，當(dāng)監(jiān)控點(diǎn)分布在如圖9所示圓柱下游紅色區(qū)域內(nèi)的剪切層中時(shí)會得到比較好的結(jié)果.圖10和圖11分別是Re=3 900流場中的監(jiān)控點(diǎn)A(0.62,0.7,0.5π)和B(0.8,0.9,0.5π)的y方向速度uy值的頻譜圖，圖12和圖13分別是Re為2 000和5 000的頻譜圖.根據(jù)圖10和圖11可以清楚地觀察到，剪切層內(nèi)的點(diǎn)A(0.62,0.7,0.5π)清晰地捕捉到了fSL的寬頻信號，而位于剪切層外的點(diǎn)B(0.8,0.9,0.5π)雖然離剪切層距離較近，但僅僅捕捉到fVS的頻率.根據(jù)對剪切層內(nèi)的流場頻譜分析，本文得出Re=2 000,3 900,5 000的fSL分別為1.051,1.463,1.584，相對應(yīng)的fSL/fVS分別為4.98,7.05,7.75.

圖9 Re=3 900的圓柱繞流剪切層示意圖Fig.9 Schematic diagram of shear layer downstream a cylinder at Re=3 900

圖10 Re=3 900流場中點(diǎn)A(0.62,0.7,0.5π)的頻譜圖Fig.10 Point A(0.62,0.7,0.5π)velocity power spectra at Re=3 900

圖11 Re=3 900流場中點(diǎn)B(0.8,0.9,0.5π)的頻譜圖Fig.11 Velocity power spectra of Point B(0.8,0.9,0.5π)at Re=3 900

圖12 Re=2 000流場中點(diǎn)A(0.62,0.7,0.5π)的頻譜圖Fig.12 Velocity power spectra of Point A(0.62,0.7,0.5π)at Re=2 000

圖13 Re=5 000流場中點(diǎn)A(0.62,0.7,0.5π)的頻譜圖Fig.13 Velocity power spectra of Point A(0.62,0.7,0.5π)at Re=5 000

2.1.2DMD的結(jié)果分析 將大渦模擬計(jì)算得到的等時(shí)間間隔的速度模值|U|和y方向的速度uy的“快照”進(jìn)行DMD分析.以Re=3 900為例，從各階模態(tài)特征值的實(shí)部和虛部的分布圖14可以看出，絕大部分的特征值位于單位圓上，因此這些模態(tài)都是周期性模態(tài)，不會隨著流動的發(fā)展而增長或者消耗.從圖15看出，頻率f=0是流場的主模態(tài)，其對應(yīng)的能量值也最大，說明其主導(dǎo)了流場的基本形態(tài).Re=2 000,3 900,5 000圓柱繞流的St≈0.21，從圖15～17中可以清晰地看到在f=0.21附近有一個(gè)孤立的頻率，這個(gè)頻率的值與卡門渦街fVS相對應(yīng).由圖15～17，在f=1.003,1.458,1.578附近也可以看到由模態(tài)特征頻率組成的一個(gè)“寬峰”，這與頻譜圖中觀測到的寬頻信號fSL范圍一致.Re=2 000,3 900,5 000相對應(yīng)的fSL/fVS分別為4.63,7.03,7.52.

圖14 Re=3 900的DMD特征值分布圖Fig.14 Distribution of eigenvalues of DMD at Re=3 900

圖15 Re=2 000 的DMD能量譜Fig.15 DMD energy spectra at Re=2 000

圖16 Re=3 900 的DMD能量譜Fig.16 DMD energy spectra at Re=3 900

圖17 Re=5 000 的DMD能量譜Fig.17 DMD energy spectra at Re=5 000

2.1.3兩種后處理方法的比較 由圖18看出，頻譜分析和DMD的后處理方法都比較準(zhǔn)確地捕捉到fSL，并與已有的DNS結(jié)果相一致.DMD方法與頻譜分析的方法相比，無需在流場中精心布置監(jiān)測點(diǎn)即可得到fSL的寬頻信號，同時(shí)DMD特征頻率是離散的點(diǎn)，峰值fSL的選取更加方便和準(zhǔn)確，沒有位置選取這樣的人為因素干擾.DMD頻率離散化的特點(diǎn)為fSL的判斷提供了一個(gè)更好的手段.

圖18 fSL/fVS與Re的關(guān)系圖Fig.18 Ratio of shear layer frequency/vortex shedding frequency versus Re

2.2 K-H不穩(wěn)定性的特征流場

選取f=0,fVS,fSL分別對應(yīng)的模態(tài)M0,MVS,MSL作圖，如圖19所示.從圖19中看出，M0是基本對稱分布的流場形態(tài)，表征了因圓柱的存在阻礙流體流動而造成的整個(gè)流場波動的特征.卡門渦街對應(yīng)的模態(tài)MVS也呈現(xiàn)出對稱的特性，并沿著來流方向往下游傳播，與Re=100的卡門渦街頻率fVS所對應(yīng)的模態(tài)有許多相似之處.剪切層K-H不穩(wěn)定性MSL的模態(tài)圖在圓柱下游的剪切層區(qū)域內(nèi)捕捉到許多細(xì)小的渦結(jié)構(gòu)，這些渦團(tuán)的位置關(guān)于圓柱也顯示出一定的對稱性.剪切層兩側(cè)的速度差誘發(fā)了K-H不穩(wěn)定性，導(dǎo)致大量細(xì)微的渦旋產(chǎn)生于此，其在頻率圖上的表現(xiàn)則為寬頻信號的凸起.

比較不同雷諾數(shù)下fSL對應(yīng)的模態(tài)圖MSL，可以發(fā)現(xiàn)隨著雷諾數(shù)的增加，剪切層不穩(wěn)定性伊始的位置向著來流方向移動，更加地靠近圓柱，并且較高雷諾數(shù)下的剪切層不穩(wěn)定性消逝的位置比低雷諾數(shù)的位置也更加靠近圓柱.這個(gè)現(xiàn)象說明，雷諾數(shù)越大，剪切層K-H不穩(wěn)定性發(fā)展得越快，消逝得越早，意味著剪切層的在x方向上的跨度越小.

在剪切層區(qū)域附近，uy的波動比模值 |U|更加顯著，因此可以選擇uy作為DMD的輸入來進(jìn)一步研究剪切層模態(tài)的特性.從圖20看出不同雷諾數(shù)下的MSL模態(tài)圖表現(xiàn)出很大的相似性：在圓柱下游的兩側(cè)均有序地排列著若干個(gè)渦旋并明顯地表現(xiàn)出了一定的對稱性，渦旋所在的位置也對應(yīng)于以模值 |U|作DMD分解得到的MSL模態(tài)中的剪切層區(qū)域.不論是以模值 |U|還是uy作為DMD的輸入，模態(tài)分解得到的剪切層不穩(wěn)定性的區(qū)域是統(tǒng)一的.

剪切層區(qū)域內(nèi)各個(gè)渦旋的直徑大小差異并不大，與文獻(xiàn)[6]中某個(gè)瞬時(shí)速度矢量場中剪切層渦旋產(chǎn)生的位置和大小相吻合.速度梯度的剪切作用在剪切層處產(chǎn)生了大小不同的渦旋，由此推測剪切層K-H不穩(wěn)定性寬頻信號里的不同頻率對應(yīng)著直徑大小不同的渦旋.

圖20 不同Re下uy的剪切層模態(tài)圖Fig.20 Shear layer modes of uy at different Re values

3 結(jié)論

本文使用大渦模擬的數(shù)值計(jì)算方法，對Re=2 000,3 900,5 000的圓柱繞流問題開展了研究，通過傅里葉變換頻譜分析和DMD分析兩種不同的數(shù)據(jù)處理方法，得到了比較一致的剪切層K-H不穩(wěn)定性頻率.本文同時(shí)用DMD方法研究了不同雷諾數(shù)對剪切層模態(tài)特征的影響.通過對比分析，發(fā)現(xiàn)在剪切層附近的位置生成了若干個(gè)平行于來流方向的渦旋并向下游傳播，渦旋的位置關(guān)于y=0的軸線對稱分布.隨著雷諾數(shù)的增加，剪切層的不穩(wěn)定性發(fā)生和消逝的位置都更加地靠近上游，并且剪切層的區(qū)域變得更小.

與傅里葉變換的方法相比，DMD方法無需在流場中精心布置監(jiān)測點(diǎn)，僅通過流場快照即可得到剪切層K-H不穩(wěn)定性的發(fā)生位置、寬頻信號和模態(tài)特征等豐富的信息.當(dāng)流場中剪切層的位置未知時(shí)(即非圓柱繞流流場時(shí))，DMD的這些優(yōu)勢將更為明顯.此外DMD頻率離散化的特點(diǎn)為fSL提供了一個(gè)更好的判斷手段.

我們注意到，在目前的尾跡剪切層中，存在著由K-H不穩(wěn)定性導(dǎo)致的渦旋尺度相近的DMD模態(tài)現(xiàn)象，為了更細(xì)致地分析該現(xiàn)象，尤其是小尺度流動的模態(tài)細(xì)節(jié)，今后有必要采用直接數(shù)值模擬的方法對該問題開展更深入的研究.