陶善聰， 郝克理， 周 毅

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

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

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

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

1 數(shù)值方法

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

(1)

(2)

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

流場(chǎng)布置如圖1所示，X、Y、Z方向?qū)?yīng)的域長(zhǎng)分別為L(zhǎng)X= 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ù)見(jiàn)表1。入口和出口邊界條件采用Inflow/Outflow，法向和展向采用周期性邊界條件。時(shí)間遞進(jìn)采用三階Runge-Kutta方法求解，對(duì)流項(xiàng)采用四階中心差分求解，采用六階中心差分緊致格式(X方向)和譜方法(Y、Z方向)求解黏性項(xiàng)，壓力項(xiàng)通過(guò)泊松方程求解。方柱模型采用浸入式邊界法。為精確模擬邊界層附近的小尺度運(yùn)動(dòng)，方柱附近沿X方向采用函數(shù)加密，網(wǎng)格布置及其加密示意圖如圖2(a)、圖2(b)所示。

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

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

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

圖3 瞬時(shí)渦量云圖

圖7、圖8分別給出了時(shí)均流場(chǎng)流向速度與表面壓強(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é)果較為接近。演變趨勢(shì)與Bearman等[16]和Mizota等[17]相同。綜上所述，數(shù)值模擬結(jié)果與前期實(shí)驗(yàn)相吻合。

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

圖5 速度脈動(dòng)能譜

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

圖7 中心線上平均流速

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

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

3.1 快速本征正交分解

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

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

(3)

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

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

(4)

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

(5)

(6)

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

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

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

圖10 原始流場(chǎng)與重構(gòu)流場(chǎng)對(duì)比

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

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

3.2 速度梯度張量不變量

速度梯度可定義為

(7)

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

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

(8)

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

P=-Si

(9)

(10)

(11)

圖13 第四、第五模態(tài)系數(shù)隨時(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的單方柱繞流特性，通過(guò)將斯特勞哈爾數(shù)、平均流速和表面壓強(qiáng)系數(shù)與文獻(xiàn)中結(jié)果[11-17]對(duì)比，驗(yàn)證數(shù)值方法的正確性。進(jìn)一步采用本征正交分解方法提取湍流場(chǎng)中不同尺度結(jié)構(gòu)，分析相干結(jié)構(gòu)對(duì)湍流流動(dòng)的貢獻(xiàn)。最后借助(R，Q)和(QW，-QS)聯(lián)合概率密度函數(shù)探討流場(chǎng)不同流向位置處的湍流流動(dòng)機(jī)理，結(jié)論如下。

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

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