桑文慧，孫志強，周孑民

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

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

1 計算模型

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

圖1 計算模型Fig.1 Calculation model

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

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

2 數(shù)值格式

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

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

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

3 斯特勞哈爾數(shù)

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

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

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

該式在47＜Re＜1 000時誤差為0.000 5。

本文采用式(2)作為斯特勞哈爾數(shù)參考值的計算式，通過比較斯特勞哈爾數(shù)的仿真值與式(2)計算所得的參考值之間的相對誤差來檢驗模擬結(jié)果的準確性。

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

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

4.1 算法選擇

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

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

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

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

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

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

4.2 流場分析

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

圖4所示為流道內(nèi)y=0.02D直線上不同雷諾數(shù)下的軸向靜壓分布，圓柱中心位于x/D=3。不同雷諾數(shù)下的軸向靜壓變化趨勢基本相同。當x/D＜1.6時，流動幾乎不受圓柱影響，靜壓保持為常數(shù)；當1.6＜x/D＜2.6時，靜壓逐漸增大；當x/D＞2.6以后，隨著圓柱邊界層分離，靜壓出現(xiàn)波動，并且隨流速增加，靜壓變化幅度增大，在x/D=3.7附近壓降達到最大。流速越大，軸向壓降最大點越靠近圓柱。在2.6＜x/D＜12范圍內(nèi)靜壓波動較顯著，隨著旋渦向下游運動，靜壓波動性逐漸減弱，直至穩(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ù)為200，300和400的情形，靜壓幅度變化較小，壓降最大值出現(xiàn)在圓柱后部上側(cè)。對于雷諾數(shù)為500和600的情形，靜壓幅度變化較大，最大壓降位于圓柱后部下側(cè)y/D=0.47附近。

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

5 結(jié)論

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

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

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

[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]孫志強,張宏建,周孑民.渦街流量計內(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].浙江大學學報: 工學版,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]王亞玲,劉應中,繆國平.圓柱繞流的三維數(shù)值模擬[J].上海交通大學學報,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.