基于大渦模擬的圓柱繞流剪切層不穩(wěn)定性

郭志遠，虞培祥,b，歐陽華,b

(上海交通大學 a.機械與動力工程學院；b.燃氣輪機與民用航空發(fā)動機教育部工程研究中心,上海 200240)

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

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

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

1 研究方法

1.1 數值方法與驗證

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

(1)

(2)

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

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

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

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

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

圖2 Re=3 900圓柱繞流的流場軸線處ux沿流向分布的時均圖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的尾跡不同位置的時均流向速度剖面圖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 數據后處理DMD方法與驗證

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

根據特征值能量的大小提取DMD前5階的模態(tài)與文獻[24]的Re=100的DMD前5階模態(tài)作對比，結果如圖8所示.圖8左側由上到下依次是文獻[24]前5階模態(tài)(M1、M2、M3、M4、M5)，圖8右側由上到下依次是本文得出的前5階模態(tài).由圖8容易看出，兩組DMD結果均捕捉到了圓柱下游尾跡處復雜的渦結構，每一階模態(tài)的流場特征都是對稱分布的，同時在圓柱下游的遠處也都捕捉到了渦系的分布特征.兩組DMD的結果基本一致，也證明了本文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 結果與分析

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

2.1.1傳統(tǒng)的頻譜分析結果 剪切層的速度梯度較大，是一個位置狹小的窄帶，fSL對監(jiān)控點的分布位置比較敏感，因此本文布置大量的監(jiān)測點來捕捉fSL.通過比較不同位置處的速度的能量譜密度α，當監(jiān)控點分布在如圖9所示圓柱下游紅色區(qū)域內的剪切層中時會得到比較好的結果.圖10和圖11分別是Re=3 900流場中的監(jiān)控點A(0.62,0.7,0.5π)和B(0.8,0.9,0.5π)的y方向速度uy值的頻譜圖，圖12和圖13分別是Re為2 000和5 000的頻譜圖.根據圖10和圖11可以清楚地觀察到，剪切層內的點A(0.62,0.7,0.5π)清晰地捕捉到了fSL的寬頻信號，而位于剪切層外的點B(0.8,0.9,0.5π)雖然離剪切層距離較近，但僅僅捕捉到fVS的頻率.根據對剪切層內的流場頻譜分析，本文得出Re=2 000,3 900,5 000的fSL分別為1.051,1.463,1.584，相對應的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流場中點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流場中點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流場中點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流場中點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的結果分析 將大渦模擬計算得到的等時間間隔的速度模值|U|和y方向的速度uy的“快照”進行DMD分析.以Re=3 900為例，從各階模態(tài)特征值的實部和虛部的分布圖14可以看出，絕大部分的特征值位于單位圓上，因此這些模態(tài)都是周期性模態(tài)，不會隨著流動的發(fā)展而增長或者消耗.從圖15看出，頻率f=0是流場的主模態(tài)，其對應的能量值也最大，說明其主導了流場的基本形態(tài).Re=2 000,3 900,5 000圓柱繞流的St≈0.21，從圖15～17中可以清晰地看到在f=0.21附近有一個孤立的頻率，這個頻率的值與卡門渦街fVS相對應.由圖15～17，在f=1.003,1.458,1.578附近也可以看到由模態(tài)特征頻率組成的一個“寬峰”，這與頻譜圖中觀測到的寬頻信號fSL范圍一致.Re=2 000,3 900,5 000相對應的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的后處理方法都比較準確地捕捉到fSL，并與已有的DNS結果相一致.DMD方法與頻譜分析的方法相比，無需在流場中精心布置監(jiān)測點即可得到fSL的寬頻信號，同時DMD特征頻率是離散的點，峰值fSL的選取更加方便和準確，沒有位置選取這樣的人為因素干擾.DMD頻率離散化的特點為fSL的判斷提供了一個更好的手段.

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

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

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

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

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

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

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

3 結論

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

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

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