有限流道內(nèi)低雷諾數(shù)二維圓柱繞流數(shù)值模擬

桑文慧，孫志強(qiáng)，周孑民

(中南大學(xué) 能源科學(xué)與工程學(xué)院，湖南 長沙，410083)

鈍體繞流，即流體橫向流過非流線型物體，廣泛存在于自然界和工業(yè)生產(chǎn)之中。鈍體繞流涉及流動(dòng)分離、旋渦脫落及其相互作用等諸多問題，尾跡流場(chǎng)極為復(fù)雜，不少流動(dòng)機(jī)理尚未完全認(rèn)識(shí)清楚[1?2]，因而對(duì)鈍體繞流的研究具有十分重要的學(xué)術(shù)價(jià)值。同時(shí)，對(duì)于存在繞流的工業(yè)設(shè)備，流體還可能誘發(fā)鈍體產(chǎn)生振動(dòng)，使鈍體疲勞損傷導(dǎo)致設(shè)備故障[3]，因此，研究鈍體繞流對(duì)提高現(xiàn)代工業(yè)設(shè)備長期運(yùn)行的可靠性與安全性也具有特別重要的現(xiàn)實(shí)意義。鈍體繞流是強(qiáng)烈的非線性時(shí)變湍流，難以用解析的方法求得其流場(chǎng)分布，所以，至今人們對(duì)鈍體繞流的認(rèn)識(shí)幾乎全部依賴于經(jīng)驗(yàn)和實(shí)驗(yàn)。計(jì)算機(jī)數(shù)值模擬技術(shù)特別是計(jì)算流體力學(xué)的飛速發(fā)展為鈍體繞流研究提供了一種有效的新方法。在模擬對(duì)象方面，Kahawita等[4]模擬了梯形鈍體繞流，討論了鈍體形狀對(duì)旋渦脫落的影響；Loc等[5]模擬了圓柱繞流中的二次渦結(jié)構(gòu)；Blackburn等[6]模擬了二維振蕩圓柱繞流的尾跡結(jié)構(gòu)；Yu等[7]模擬了多孔方柱繞流，討論了雷諾數(shù)和達(dá)西數(shù)對(duì)旋渦尺寸和位置的影響；孫志強(qiáng)等[8]通過對(duì)渦街流量計(jì)內(nèi)鈍體繞流的模擬，揭示了在鈍體后部一定范圍內(nèi)壁面壓力分布規(guī)律。在模擬方法方面，Breuer[9]對(duì)雷諾數(shù)為140 000的圓柱繞流做了大渦模擬；Yang等[10]采用動(dòng)網(wǎng)格和Galerkin有限元法模擬了有限流道內(nèi)振蕩矩形柱繞流的尾跡結(jié)構(gòu)；黃鈺期等[11]使用分塊耦合求解方法模擬了圓柱繞流；王亞玲等[12]采用SIMPLE算法和QUICK格式對(duì)圓柱繞流進(jìn)行了計(jì)算，發(fā)現(xiàn)三維模擬的升力和阻力系數(shù)均小于二維模擬；Cao等[13]對(duì)剪切流的圓柱繞流進(jìn)行了直接數(shù)值模擬和大渦模擬，分析了雷諾數(shù)的影響；Rocchi等[14]比較了圓柱繞流的實(shí)驗(yàn)與數(shù)值模擬結(jié)果。上述研究豐富了人們對(duì)鈍體繞流的認(rèn)識(shí)。工業(yè)上鈍體繞流常發(fā)生在有限流道內(nèi)，并且大多數(shù)鈍體的長度遠(yuǎn)大于其寬度，因此，研究有限流道內(nèi)二維鈍體繞流具有重要的實(shí)際意義。然而，目前針對(duì)有限流道內(nèi)鈍體繞流的研究仍比較少。本文作者以計(jì)算流體力學(xué)軟件FLUENT為平臺(tái)，對(duì)有限流道內(nèi)低雷諾數(shù)二維圓柱繞流開展數(shù)值模擬，著重討論壓力?速度耦合算法、壓力離散格式和動(dòng)量方程離散格式對(duì)模擬結(jié)果的影響，尋找較佳的數(shù)值格式組合，并根據(jù)模擬結(jié)果對(duì)圓柱繞流流場(chǎng)進(jìn)行分析。

1 計(jì)算模型

當(dāng)雷諾數(shù)較小時(shí)，流體緊貼圓柱體表面流動(dòng)，尾跡中無旋渦；隨著雷諾數(shù)增大，尾跡中逐漸形成周期的旋渦脫落現(xiàn)象。圓柱繞流過程受到流體流動(dòng)狀況和鈍體結(jié)構(gòu)參數(shù)等因素的影響。本文研究的二維圓柱繞流計(jì)算模型如圖1所示。流體沿x方向流經(jīng)流道內(nèi)的圓柱鈍體，圓柱直徑d=14 mm，流道寬度D=50 mm，流道長 20D，圓柱中心距流道入口和出口分別為 3D和17D。由于流道壁面壓力便于實(shí)驗(yàn)測(cè)量與驗(yàn)證，故設(shè)置P1(160,1)和P2(160,49)為監(jiān)視點(diǎn)。流動(dòng)介質(zhì)為水，密度為998.2 kg/m3，動(dòng)力黏度為1.003 g/(m·s)。

圖1 計(jì)算模型Fig.1 Calculation model

考慮到圓柱繞流流場(chǎng)的特點(diǎn)，采用多塊結(jié)構(gòu)化網(wǎng)格，將計(jì)算區(qū)域劃分為圓柱周圍區(qū)①、圓柱前區(qū)②、近壁面區(qū)③和尾跡區(qū)④等4個(gè)子區(qū)域，如圖2所示。由于旋渦的產(chǎn)生和脫落發(fā)生在圓柱表面，故邊長為0.8D的正方形圓柱周圍區(qū)①采用正交性較好的輻射型網(wǎng)格。圓柱前和靠近壁面處的流場(chǎng)較平穩(wěn)，為節(jié)省計(jì)算時(shí)間，圓柱前區(qū)②采用較稀疏的四邊形網(wǎng)格；近壁面區(qū)③也采用四邊形網(wǎng)格，但較圓柱前區(qū)②稍密集。在圓柱尾跡中旋渦強(qiáng)度沿流動(dòng)方向由強(qiáng)變?nèi)?，故尾跡區(qū)④采用由密漸疏型網(wǎng)格。流道內(nèi)壁和圓柱表面網(wǎng)格設(shè)置為邊界層，并滿足無滑移條件。計(jì)算區(qū)域內(nèi)的網(wǎng)格總數(shù)為44 620，最小和最大網(wǎng)格的邊長分別為0.1 mm和18.4 mm。各求解變量收斂標(biāo)準(zhǔn)設(shè)置為殘差值小于10?8。為了使每個(gè)旋渦脫落周期內(nèi)的采樣點(diǎn)不小于50，時(shí)間步長的設(shè)置隨入口速度不同而相應(yīng)改變。流道入口設(shè)置為速度入口，出口設(shè)置為壓力出口。模擬的雷諾數(shù)為 200～600，在該范圍內(nèi)流動(dòng)為層流，故采用二維非定常不可壓縮N-S方程組。

圖2 多塊輻射型網(wǎng)格Fig.2 Multi-block radial gird

2 數(shù)值格式

FLUENT將控制方程轉(zhuǎn)換為代數(shù)方程來求解，通過對(duì)每一個(gè)控制體積分，得到基于控制體的守恒離散方程，故合理選擇離散格式至關(guān)重要。此外，由于壓力是圓柱繞流求解的重要參數(shù)，數(shù)值計(jì)算中需要根據(jù)動(dòng)量方程和連續(xù)性方程推導(dǎo)出壓力方程，因此，需合理選擇壓力?速度耦合算法。

壓力離散格式主要有 Standard、二階離散、PRESTO、線性和體積力加權(quán)格式。由于有限流道內(nèi)低雷諾數(shù)圓柱繞流控制體中心的壓力梯度變化不大，且計(jì)算時(shí)采用了四邊形網(wǎng)格，因此，Standard、二階離散和PRESTO格式等較為合適。動(dòng)量方程離散格式主要有一階迎風(fēng)、二階迎風(fēng)、冪率和 QUICK格式，根據(jù)圓柱繞流的特點(diǎn)，一階迎風(fēng)、二階迎風(fēng)和QUICK格式較為適合。壓力?速度耦合算法主要有SIMPLE，SIMPLEC和PISO算法。由于進(jìn)出流道的流體質(zhì)量守恒，因此，3種方法均比較合適。但由于控制方程速度校正項(xiàng)的不同，使得三者在計(jì)算速度上存在一定的差異。

綜上所述，適于有限流道內(nèi)低雷諾數(shù)圓柱繞流的壓力?速度耦合算法、壓力離散格式和動(dòng)量方程離散格式各有3種，通過組合可以有27種求解組合。本文首先確定較優(yōu)的組合方法，實(shí)現(xiàn)對(duì)有限流道內(nèi)低雷諾數(shù)二維圓柱繞流計(jì)算精度和速度俱佳的模擬。為了便于后續(xù)描述，不同數(shù)值算法采用如表1所示的代碼表示。

3 斯特勞哈爾數(shù)

斯特勞哈爾數(shù)St是用于表征旋渦脫落特性的最重要的相似準(zhǔn)則數(shù)之一，可視為旋渦脫落的無量綱頻率。因此，本文以斯特勞哈爾數(shù)來判斷數(shù)值模擬結(jié)果的準(zhǔn)確性。對(duì)于有限流道內(nèi)圓柱繞流，St的計(jì)算式為：

式中：m為圓柱兩側(cè)弓形流通面積與流道橫截面積之比；d為圓柱直徑；f為旋渦脫落頻率；U為來流速度。本文中m= 0.648。

人們對(duì)圓柱繞流斯特勞哈爾數(shù)St與雷諾數(shù)Re的關(guān)系做了大量研究。目前比較公認(rèn)的St與Re經(jīng)驗(yàn)關(guān)系式為[15]：

該式在47＜Re＜1 000時(shí)誤差為0.000 5。

本文采用式(2)作為斯特勞哈爾數(shù)參考值的計(jì)算式，通過比較斯特勞哈爾數(shù)的仿真值與式(2)計(jì)算所得的參考值之間的相對(duì)誤差來檢驗(yàn)?zāi)M結(jié)果的準(zhǔn)確性。

表1 數(shù)值算法代碼Table 1 Code of numerical algorithms

4 計(jì)算結(jié)果與討論

4.1 算法選擇

本文模擬了表1所列的不同壓力?速度耦合算法、壓力離散格式和動(dòng)量方程離散格式時(shí)有限流道內(nèi)二維圓柱繞流流場(chǎng)。通過FLUENT求解控制方程，得到監(jiān)視點(diǎn)P1(160,1)和P2(160,49)的靜壓，再對(duì)靜壓進(jìn)行傅里葉分析，求得旋渦脫落頻率，最后由式(1)和(2)計(jì)算得到斯特勞哈爾數(shù)的仿真值和參考值。表2給出了不同數(shù)值格式組合下斯特勞哈爾數(shù)仿真值與參考值的平均相對(duì)誤差。

表2 不同數(shù)值格式組合下斯特勞哈爾數(shù)的平均相對(duì)誤差Table 2 Average relative errors of Strouhal numbers under different numerical scheme combinations

從表2可以看出：A1，A2和A3在各種精度的結(jié)果中出現(xiàn)，三者對(duì)計(jì)算結(jié)果無影響，表明模擬時(shí)可在SIMPLE，SIMPLEC和PISO 3種壓力?速度耦合算法中任選其一。B1和 B3在Ⅰ，Ⅲ和Ⅴ類中出現(xiàn)，表明Standard和二階壓力離散也可任意選擇。分別比較Ⅰ與Ⅱ，Ⅲ與Ⅳ，Ⅴ與Ⅵ可以看出：在A和C相同時(shí)，Standard(B1)和二階壓力離散(B3)的計(jì)算精度高于PRESTO(B2)，這是由于 PRESTO格式不太適用于低速流動(dòng)，不同壓力離散格式引入的計(jì)算誤差相差不大。在A和B相同時(shí)，由于計(jì)算的是渦流且采用了四邊形網(wǎng)格，因此QUICK(C3)的計(jì)算精度最高。從Ⅴ和Ⅵ可知，較大計(jì)算誤差主要是由一階迎風(fēng)(C1)所引入。動(dòng)量方程離散格式對(duì)于數(shù)值模擬誤差的形成起主要作用，壓力離散格式的影響次之，這是由于二維圓柱繞流模型的求解過程主要是二維N-S動(dòng)量方程的求解。

根據(jù)以上分析，模擬精度最高的數(shù)值格式組合有6種。以下通過比較旋渦脫落穩(wěn)定時(shí)間，從計(jì)算速度方面進(jìn)一步確定較佳的數(shù)值格式組合。當(dāng)監(jiān)視點(diǎn)靜壓的某一峰值與前后兩峰值相差小于0.1%時(shí)，認(rèn)為旋渦脫落達(dá)到穩(wěn)定，故將從開始計(jì)算到此峰值的前一峰值的時(shí)間間隔作為旋渦脫落穩(wěn)定時(shí)間。表3給出了6種數(shù)值格式組合在不同雷諾數(shù)下的旋渦脫落穩(wěn)定時(shí)間。比較各個(gè)雷諾數(shù)下的結(jié)果可以看出：壓力?速度耦合算法中 SIMPLEC(A2)穩(wěn)定時(shí)間最短，這是由于SIMPLEC與SIMPLE(A1)速度校正項(xiàng)不同導(dǎo)致的計(jì)算速度不同，而 PISO(A3)在每個(gè)迭代中要花費(fèi)較多的CPU時(shí)間，雖然能減少到達(dá)收斂所需要的迭代步數(shù)，但對(duì)于有限流道內(nèi)圓柱繞流模擬，這一優(yōu)勢(shì)并不明顯。比較同一雷諾數(shù)下不同數(shù)值格式組合的穩(wěn)定時(shí)間可見，壓力離散格式中二階(B3)基本上都優(yōu)于Standard(B1)，這是因?yàn)槎A離散格式適合于流過曲面邊界的計(jì)算。

表3 A1B1C3，A2B1C3，A3B1C3，A1B3C3，A2B3C3和 A3B3C3的穩(wěn)定時(shí)間Table 3 Convergence times of A1B1C3,A2B1C3,A3B1C3,A1B3C3,A2B3C3 and A3B3C3

綜合考慮模擬精度和計(jì)算速度，對(duì)于有限流道內(nèi)低雷諾數(shù)二維圓柱繞流的計(jì)算，建議采用 SIMPLEC壓力?速度耦合，二階壓力離散和QUICK動(dòng)量方程離散格式的組合。

4.2 流場(chǎng)分析

采用SIMPLEC壓力?速度耦合、二階壓力離散和QUICK動(dòng)量方程離散格式，計(jì)算不同雷諾數(shù)下圓柱繞流的渦量，其等值線分布如圖3所示。各個(gè)雷諾數(shù)下圓柱尾跡中都出現(xiàn)了規(guī)則的旋渦脫落，在向下游運(yùn)動(dòng)的同時(shí)，旋渦強(qiáng)度不斷減弱。隨著雷諾數(shù)的增大，相同位置處的旋渦強(qiáng)度逐漸增大。

圖4所示為流道內(nèi)y=0.02D直線上不同雷諾數(shù)下的軸向靜壓分布，圓柱中心位于x/D=3。不同雷諾數(shù)下的軸向靜壓變化趨勢(shì)基本相同。當(dāng)x/D＜1.6時(shí)，流動(dòng)幾乎不受圓柱影響，靜壓保持為常數(shù)；當(dāng)1.6＜x/D＜2.6時(shí)，靜壓逐漸增大；當(dāng)x/D＞2.6以后，隨著圓柱邊界層分離，靜壓出現(xiàn)波動(dòng)，并且隨流速增加，靜壓變化幅度增大，在x/D=3.7附近壓降達(dá)到最大。流速越大，軸向壓降最大點(diǎn)越靠近圓柱。在2.6＜x/D＜12范圍內(nèi)靜壓波動(dòng)較顯著，隨著旋渦向下游運(yùn)動(dòng)，靜壓波動(dòng)性逐漸減弱，直至穩(wěn)定。

圖3 不同雷諾數(shù)下的渦量等值線圖Fig.3 Contours of vorticity at different Reynolds numbers

圖4 軸向靜壓分布Fig.4 Axial static pressure distribution

圖5所示為流道內(nèi)x=3.6D直線上不同雷諾數(shù)下的徑向靜壓分布，圓柱中心位于y/D=0.5。不同雷諾數(shù)下的徑向靜壓變化趨勢(shì)也基本相同。流速越大，徑向靜壓幅度變化越大，這是因?yàn)榱魉僭龃髸r(shí)圓柱體后生成的旋渦強(qiáng)度越大，因而引起的尾跡波動(dòng)也越顯著。對(duì)于雷諾數(shù)為200，300和400的情形，靜壓幅度變化較小，壓降最大值出現(xiàn)在圓柱后部上側(cè)。對(duì)于雷諾數(shù)為500和600的情形，靜壓幅度變化較大，最大壓降位于圓柱后部下側(cè)y/D=0.47附近。

圖5 徑向靜壓分布Fig.5 Radial static pressure distribution

5 結(jié)論

(1)采用多塊輻射型網(wǎng)格可較好地模擬圓柱繞流尾跡旋渦的產(chǎn)生與脫落。

(2)動(dòng)量方程離散格式對(duì)數(shù)值模擬誤差的形成起主要作用，壓力離散格式的影響次之。模擬精度和計(jì)算速度俱佳的有限流道內(nèi)低雷諾數(shù)二維圓柱繞流數(shù)值模擬算法為 SIMPLEC算法、二階壓力離散格式和QUICK動(dòng)量方程離散格式的組合。

(3)流速越大，軸向壓降最大點(diǎn)越靠近圓柱，徑向靜壓幅度變化越大。在2.6＜x/D＜12范圍內(nèi)靜壓波動(dòng)較顯著。

[1]Cai J S,Chng T L.On vortex shedding from bluff bodies with base cavities[J].Physics of Fluids,2009,21(3): 034109.

[2]Williamson C H K.Vortex dynamics in the cylinder wake[J].Annual Review of Fluid Mechanics,1996,28(1): 477?539.

[3]Williamson C H K,Govardhan R.Vortex-induced vibrations[J].Annual Review of Fluid Mechanics,2004,36(1): 413?455.

[4]Kahawita R,Wang P.Numerical simulation of the wake flow behind trapezoidal bluff bodies[J].Computers and Fluids,2002,31(1): 99?112.

[5]Loc T P,Bouard R.Numerical solution of the early stage of unsteady viscous flow around a circular cylinder: A comparison with experimental visualization and measurements[J].Journal of Fluid Mechanics,1985,160: 93?117.

[6]Blackburn H M,Henderson R D.A study of two-dimensional flow past an oscillating cylinder[J].Journal of Fluid Mechanics,1999,385: 255?286.

[7]Yu P,Zeng Y,Lee T S,et al.Wake structure for flow past and through a porous square cylinder[J].International Journal of Heat and Fluid Flow,2010,31(2): 141?153.

[8]孫志強(qiáng),張宏建,周孑民.渦街流量計(jì)內(nèi)壁面壓力分布的數(shù)值模擬[J].傳感器與微系統(tǒng),2008,27(10): 70?75.SUN Zhi-qiang,ZHANG Hong-jian,ZHOU Jie-min.Numerical simulation of wall pressure distribution in vortex flowmeter[J].Transducer and Microsystem Technologies,2008,27(10):70?75.

[9]Breuer M.A challenging test case for large eddy simulation:high Reynolds number circular cylinder flow[J].International Journal of Heat and Fluid Flow,2000,21(5): 648?654.

[10]Yang S J,Chang T R,Fu W S.Numerical simulation of flow structures around an oscillating rectangular cylinder in a channel flow[J].Computational Mechanics,2005,35(5): 342?351.

[11]黃鈺期,鄧見,任安祿.黏性非定常圓柱繞流的升阻力研究[J].浙江大學(xué)學(xué)報(bào): 工學(xué)版,2003,37(5): 596?601.HUANG Yu-qi,DENG Jian,REN An-lu.Research on lift and drag in unsteady viscous flow around circular cylinders[J].Journal of Zhejiang University: Engineering Science,2003,37(5): 596?601.

[12]王亞玲,劉應(yīng)中,繆國平.圓柱繞流的三維數(shù)值模擬[J].上海交通大學(xué)學(xué)報(bào),2001,35(10): 118?122.WANG Ya-ling,LIU Ying-zhong,MIU Guo-ping.Three-dimensional numerical simulation of viscous flow around circular cylinder[J].Journal of Shanghai Jiaotong University,2001,35(10): 118?122.

[13]Cao S Y,Ozono S,Tamura Y,et al.Numerical simulation of Reynolds number effects on velocity shear flow around a circular cylinder[J].Journal of Fluids and Structures,2010,26(5):685?702.

[14]Rocchi D,Zasso A.Vortex shedding from a circular cylinder in a smooth and wired configuration: Comparison between 3D LES simulation and experimental analysis[J].Journal of Wind Engineering and Industrial Aerodynamics,2002,90(4/5):475?489.

[15]Henderson R.Nonlinear dynamics and pattern formation in turbulent wake transition[J].Journal of Fluid Mechanics,1997,352(1): 65?112.