陶善聰， 郝克理， 周 毅

(南京理工大學(xué)能源與動力工程學(xué)院，南京 210094)

繞流運(yùn)動是指流體繞物體外部的流動。在日常生活和工程應(yīng)用中常見此類問題，如氣流吹過建筑群、海水流過海洋平臺等。由于邊界層的分離，固體表面會有漩渦交替生成、脫落、發(fā)展和消亡。當(dāng)漩渦脫落頻率與固體障礙物自振頻率一致時，易產(chǎn)生共振而對物體自身造成損毀，因此研究繞流問題意義重大。由于方柱繞流的分離點(diǎn)固定，且制造工藝和網(wǎng)格生成簡便，因此方柱被廣泛用于研究繞流問題。

近年來，眾多學(xué)者對方柱繞流問題開展了廣泛研究。外國開展此類問題研究起步較早。Lyn等[1]通過激光多普勒測速儀定量分析單方柱繞流(雷諾數(shù)Re=21 400)的湍流統(tǒng)計(jì)特性，精確的實(shí)驗(yàn)結(jié)果成為后續(xù)學(xué)者驗(yàn)證實(shí)驗(yàn)和數(shù)值計(jì)算結(jié)果的依據(jù)。實(shí)際流動是非穩(wěn)態(tài)的，Kurtulus等[2]通過粒子圖像測速儀研究單方柱在Re=4 900下的非穩(wěn)態(tài)氣動力，指出無黏流下的非穩(wěn)態(tài)項(xiàng)可以忽略且二維流動下的非穩(wěn)態(tài)氣動力可以通過時間分辨粒子圖像測速(time-resolved particle image velocimetry，TR-PIV)精確測量。最近許多研究結(jié)果表明Kolmogorov的論述并不完全適用，-5/3能譜冪次規(guī)律同樣可以在非均勻和各向異性的剪切流或邊界層中出現(xiàn)。Portela等[3]利用KHMH(Kármán-Howarth-Monin-Hill)方程研究繞流單方柱下游的能量傳遞行為，發(fā)現(xiàn)非均勻和各項(xiàng)異性湍流區(qū)出現(xiàn)的-5/3能譜伴隨著能量反向串級行為的發(fā)生，證實(shí)在剪切流或邊界層中能量不僅由大尺度湍渦向小尺度湍渦傳遞且存在由小尺度湍渦向大尺度湍渦傳遞。鈍體繞流問題涉及到湍流流動的不穩(wěn)定性，對其進(jìn)行深入研究能夠增強(qiáng)關(guān)于流動本質(zhì)的認(rèn)識。例如，單方柱繞流的遠(yuǎn)場湍動能預(yù)算[4]、串列并列雙方柱[5-6]和三方柱繞流的本質(zhì)特征等[7]研究均較好地反映出流場特性。工程實(shí)踐當(dāng)中，人們關(guān)注的是繞流柱體的氣動力特性。郜陽等[8]采用大渦模擬方法比較不同阻塞率下二維方柱繞流的整體氣動力和表面風(fēng)壓，結(jié)果表明側(cè)面及背風(fēng)面負(fù)風(fēng)壓系數(shù)和柱體兩側(cè)流場特性對阻塞比變化響應(yīng)明顯。杜曉慶等[9]著重探討了圓角化方形截面柱體的流場作用機(jī)理，發(fā)現(xiàn)方柱氣動力減小和斯特勞哈爾數(shù)(St)增大是由尾流寬度減小和旋渦脫落強(qiáng)度削弱造成的。

現(xiàn)擬從單方柱繞流出發(fā)，通過直接數(shù)值模擬(direct numerical simulation, DNS)方法研究湍流場物理特性。基于快速本征正交分解(snap-proper orthogonal decomposition,Snap-pod)將瞬時流場分解為不同含能尺度結(jié)構(gòu)的模態(tài)，定性地分析多尺度渦結(jié)構(gòu)對流場的貢獻(xiàn)；借助速度梯度張量不變量(R,Q)和(QW, -QS)聯(lián)合概率密度函數(shù)(joint probability density function, JPDF)探討流場不同流向位置處的湍流場特性。

1 數(shù)值方法

直接數(shù)值模擬方法的控制方程為連續(xù)性方程和不可壓縮Navier-Stokes方程組，張量表達(dá)式如下：

(1)

(2)

式中：u為速度；ρ為密度；P為壓力；υ為運(yùn)動黏度；t為時間；X、Y、Z分別為流向、法向和展向方向。

流場布置如圖1所示，X、Y、Z方向?qū)?yīng)的域長分別為LX= 18T0，LY= 16T0，LZ= 16T0，T0為方柱流向和法向尺寸。方柱中心布置在距離入口Xbar= 8T0處，入口雷諾數(shù)為Rein=UinT0/v=1 200，其中Uin為入口流速，為3.75 m·s-1。阻塞率σ=T0/LY×100%= 6.25%，各方向網(wǎng)格量分別為NXNYNZ=421×337×337。具體計(jì)算參數(shù)見表1。入口和出口邊界條件采用Inflow/Outflow，法向和展向采用周期性邊界條件。時間遞進(jìn)采用三階Runge-Kutta方法求解，對流項(xiàng)采用四階中心差分求解，采用六階中心差分緊致格式(X方向)和譜方法(Y、Z方向)求解黏性項(xiàng)，壓力項(xiàng)通過泊松方程求解。方柱模型采用浸入式邊界法。為精確模擬邊界層附近的小尺度運(yùn)動，方柱附近沿X方向采用函數(shù)加密，網(wǎng)格布置及其加密示意圖如圖2(a)、圖2(b)所示。

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

方柱繞流某一瞬時渦量場如圖3所示。從圖3中可以清晰地看到旋渦分別從方柱上下表面脫落而向下游發(fā)展，在下游形成交替的渦街。圖4給出了流場網(wǎng)格的空間分辨率。可以看出，流場中最大空間分辨率為(ΔXΔYΔZ)1/3/η=3.3，表現(xiàn)在X/T0=2.1左右，方柱下游出口段空間分辨率最小，其值為(ΔXΔYΔZ)1/3/η=2.3。其中η為Kolmogorov長度尺度，是湍流場中的最小尺度。Laizet等[10]指出，當(dāng)流場空間分辨率小于5時，利用DNS求解一階及二階動量的結(jié)果與實(shí)驗(yàn)相比誤差不超過5%，因此計(jì)算精度滿足空間分辨率要求。圖5(a)、圖5(b)分別展示了中心線上不同位置處(X/T0=2, 4, 6, 8)流向和法向速度脈動能譜在頻域的變化?？梢钥吹?，流向和法向速度脈動能譜在頻域演化趨勢一致，但峰值點(diǎn)對應(yīng)的fT0/Uin不同(f為取樣率)。法向脈動能譜極值點(diǎn)為fT0/Uin=0.14，等于單方柱繞流的漩渦脫落頻率斯特勞哈爾數(shù)(St)。流向速度脈動能譜極值點(diǎn)fT0/Uin=0.28為漩渦脫落頻率St的兩倍。圖6為St與Davis等[11]的實(shí)驗(yàn)結(jié)果對比，結(jié)果基本吻合。

圖2 網(wǎng)格及加密示意圖

圖3 瞬時渦量云圖

圖7、圖8分別給出了時均流場流向速度與表面壓強(qiáng)系數(shù)分布，其中黑色空心方框表示方柱模型。

由于各實(shí)驗(yàn)布置條件有所差別，最大回流速度略小于Hu等[12]、Lyn 等[1]和Durao 等[13]的計(jì)算結(jié)果，恢復(fù)速度與低雷諾數(shù)下實(shí)驗(yàn)Lee等[14]結(jié)果非常吻合，而方柱表面壓強(qiáng)分布與Chen等[15]實(shí)驗(yàn)結(jié)果較為接近。演變趨勢與Bearman等[16]和Mizota等[17]相同。綜上所述，數(shù)值模擬結(jié)果與前期實(shí)驗(yàn)相吻合。

圖4 空間分辨率的流向演變示意圖

圖5 速度脈動能譜

圖6 斯特勞哈爾數(shù)(St)與實(shí)驗(yàn)對比

圖7 中心線上平均流速

圖8 方柱表面壓強(qiáng)系數(shù)

3 結(jié)果分析與討論

3.1 快速本征正交分解

現(xiàn)通過快速本征正交分解方法將某一時刻的流場分解為具有不同漩渦特征的r個模態(tài)以對比分析方柱繞流場中的相干結(jié)構(gòu)。其基本思想是，通過包含N(1, 2, …,N)個時間步流場的矩陣U求解相關(guān)矩陣C，再找出前r個模態(tài)φ，這r個模態(tài)累積的能量基本等效于全空間的總能量。具體求法如下。

通過相對能量ε(r)的大小確定模態(tài)個數(shù)r的取值：

(3)

式(3)中：λi是時間相關(guān)矩陣C的特征值從大到小的排列。

再計(jì)算每一模態(tài)：

(4)

最后通過每個模態(tài)的系數(shù)α重構(gòu)原始流場。

(5)

(6)

相關(guān)矩陣特征值λ與相對能量ε隨模態(tài)r的累積如圖9所示。由圖9可知，前50個模態(tài)累積的能量達(dá)到總能量的99.99%。故提取前50個模態(tài)重構(gòu)任意時刻流向速度場，該時刻重構(gòu)流場與原始流場對比如圖10所示。不難發(fā)現(xiàn)，前50個模態(tài)基本已提取到原始流場全部的流向速度細(xì)節(jié)，也側(cè)面驗(yàn)證了所用方法的準(zhǔn)確及精確性。由于利用本征正交分解(proper orthogonal decomposition，POD)分解速度場，因而第一階模態(tài)代表平均速度場的特征，剩余模態(tài)對應(yīng)分解脈動速度場的特征。前五模態(tài)的速度等值線圖如圖11所示。第一階模態(tài)沿Y/T0=0對稱分布，柱體上下表面附近產(chǎn)生加速區(qū)而回流區(qū)在柱體下游生成，這表明流動隨時間呈現(xiàn)周期性變化。第二和第三階模態(tài)沿Y/T0=0表現(xiàn)出反對稱分布，而沿流向生成正負(fù)交替的渦街。值得注意的是，第三階模態(tài)從位置和形狀來看類似于第二階模態(tài)渦街的“成長”。第四和第五模態(tài)速度等值線圖與前三個模態(tài)相比更加紊亂，小尺度旋渦充滿整個流場，且與文獻(xiàn)[18]的結(jié)果存在較大出入。這可能由于文獻(xiàn)[18]模擬的是Re= 100的二維層流工況，無法提取湍流場中的多尺度旋渦結(jié)構(gòu)。

第二、第三、第四、第五階模態(tài)的系數(shù)α隨時間演化分別如圖12和圖13所示。圖12中第二、三階模態(tài)系數(shù)隨時間呈現(xiàn)出正弦函數(shù)演變，且脫落周期為St1=1/(t2-t1)=1/(t4-t3)≈0.14，恰好等于圖5所求的St，這更加證明第二、三階模態(tài)提取的是流場大尺度旋渦脫落特征。圖13中第四、五階模態(tài)的旋渦脫落周期顯然小于第二、三階模態(tài)流場中的旋渦脫落周期，結(jié)合速度等值線圖可知這兩個模態(tài)提取的是流場中的高頻小尺度旋渦。圖14給出了(第二、第三)和(第四、第五)模態(tài)系數(shù)的利薩如圖形。模態(tài)系數(shù)的分散程度實(shí)際上代表大尺度相干結(jié)構(gòu)的周期性強(qiáng)弱，顯而易見第二、第三模態(tài)周期性更強(qiáng)。且利薩如圖圓形分布半徑的平方與前兩個模態(tài)對總能量的貢獻(xiàn)成正比，表明小尺度運(yùn)動(第四、五模態(tài))蘊(yùn)含的能量較弱。

圖9 相關(guān)矩陣特征值λ與相對能量ε隨模態(tài)r的累積

圖10 原始流場與重構(gòu)流場對比

圖11 各模態(tài)流向速度等值線圖

圖12 第二、第三模態(tài)系數(shù)隨時間演變

3.2 速度梯度張量不變量

速度梯度可定義為

(7)

式(7)中：V0為初始速度；V為終止速度。

速度梯度張量的特征方程為

(8)

式(8)中：λi為速度梯度張量的特征值；第一不變量P、第二不變量Q和第三不變量R可表達(dá)為

P=-Si

(9)

(10)

(11)

圖13 第四、第五模態(tài)系數(shù)隨時間演變

圖14 模態(tài)系數(shù)利薩如圖

圖15 (Q , R),(-QS,QW) 拓?fù)鋱D[19]

圖16 各流向位置(Q, R), (-QS,QW) JPDF云圖

4 結(jié)論

采用直接數(shù)值模擬方法研究Re=1 200的單方柱繞流特性，通過將斯特勞哈爾數(shù)、平均流速和表面壓強(qiáng)系數(shù)與文獻(xiàn)中結(jié)果[11-17]對比，驗(yàn)證數(shù)值方法的正確性。進(jìn)一步采用本征正交分解方法提取湍流場中不同尺度結(jié)構(gòu)，分析相干結(jié)構(gòu)對湍流流動的貢獻(xiàn)。最后借助(R，Q)和(QW，-QS)聯(lián)合概率密度函數(shù)探討流場不同流向位置處的湍流流動機(jī)理，結(jié)論如下。

(1)對于速度場的模態(tài)分解，第一模態(tài)代表平均速度場的特征；第二、第三階模態(tài)提取的是流場中的低頻大尺度旋渦特征；第四、第五模態(tài)提取的是流場中的高頻小尺度旋渦特征。第二、第三模態(tài)均有強(qiáng)周期性且蘊(yùn)含了大部分流場能量。

(2)通過對速度梯度張量的第二不變量Q和第三不變量R的分析，柱體下游大致可分為兩個流動階段：發(fā)展階段(X/T0=0.51-1.9)，流體以渦流層結(jié)構(gòu)和耗散作用為主，渦流管結(jié)構(gòu)逐漸生成；成熟階段(X/T0=1.9-9.0)，流場中湍流結(jié)構(gòu)高渦量擬能和高能量耗散率并存。