低雷諾數(shù)圓柱繞流特性的數(shù)值模擬

吳言超， 陸曉峰

( 南京工業(yè)大學 機械與動力工程學院，江蘇 南京 210009 )

0 引言

鈍體繞流現(xiàn)象在航空航天、化工、海洋石油工程、核工程等領(lǐng)域普遍存在，如化工生產(chǎn)中廣泛應用的圓柱形鋼制塔器和管殼式換熱器等設備，由鈍體繞流造成的經(jīng)濟損失日益引起人們的重視.圓柱繞流作為鈍體繞流的特例，在工程中表現(xiàn)頻繁，成為人們的研究對象.目前，對圓柱繞流的研究主要集中于遠場尾渦和某一特定雷諾數(shù)近壁面附近流動特性[1-7]，對其近場水動力學特性隨雷諾數(shù)變化規(guī)律未見報道.

采用有限體積法、分區(qū)結(jié)構(gòu)化網(wǎng)格和層流模型(Laminar)求解Navier-Stokes方程，對低雷諾數(shù)(30≤Re≤300)時的圓柱繞流進行數(shù)值模擬，再現(xiàn)圓柱體后部旋渦生成、長大和脫落的發(fā)展過程；分析圓柱繞流近壁面水動力學特性，并對比計算結(jié)果與實驗結(jié)果.

1 數(shù)學模型

1.1 控制方程及離散格式

在直角坐標系的xoy平面內(nèi)建立方程組，控制方程包含質(zhì)量守恒的連續(xù)性方程和動量方程：

式中：u，v分別為流體速度矢量在x、y方向的速度分量；ρ、υ分別為流體的密度及運動黏度.

采用有限體積法離散微分方程，為使計算穩(wěn)定加速收斂，先求解定常流動，并將該值作為非定常流動的初始值.分別采用SIMPLE算法、PISO算法求解定常流動時、非定常流動時的離散方程；時間推進采用二階隱式格式；鑒于精度和收斂速度考慮，空間離散采用二階迎風格式.當雷諾數(shù)大于300時，圓柱體后面的渦街已完全轉(zhuǎn)變?yōu)橥牧?，但圓柱體表面上的邊界層為層流[2]，主要研究圓柱體附近流場特性；因此采用層流(Laminar)模型計算.

1.2 計算區(qū)域和網(wǎng)格劃分

計算區(qū)域以圓柱體直徑D為特征尺度，為了消除邊界對圓柱體周圍流場的影響，圓柱體上游取10D，下游取30D，上、下壁面各取10D(見圖1)；為生成較好的結(jié)構(gòu)化網(wǎng)格對計算區(qū)域分區(qū)以及圓柱體周圍10D×10D范圍進行網(wǎng)格加密，獲得正交性較好的高質(zhì)量網(wǎng)格，加速計算收斂；壁面最小網(wǎng)格間距為3.927×10-4m.

圖1 計算區(qū)域劃分(單位：m)

1.3 邊界條件和初始條件

邊界條件：進口條件為均勻來流u=u0，v=0；出口條件為連續(xù)出流條件(outflow)；圓柱體表面及上下固壁為無滑移壁面邊界.以相應雷諾數(shù)下進口速度邊界作為初始條件.

2 數(shù)值計算結(jié)果

2.1 定常渦

當Re=30時，圓柱體后尾跡區(qū)是對稱的定常渦，尾渦尺寸L與圓柱體直徑D之比為1.50(見圖2)，與文獻[7]結(jié)果吻合.通過實驗測量和數(shù)值模擬證實，定常流動失穩(wěn)的臨界Re為40～50[7-10]，實驗條件不同、柱體后近尾跡不穩(wěn)定性及圓柱體縱橫比差異導致臨界Re存在差異.

圓柱體表面不同時刻的無因次壓力沿壁面角θ分布見圖3.由圖3可見，壓力因數(shù)Cp呈現(xiàn)較好的對稱性，與Dennis S C R等獲得的變化規(guī)律一致[11].不同周期T的壓力因數(shù)分布曲線能夠較好重合，表明Re=30時對稱渦具有穩(wěn)定性；最低壓力點位于±90°，說明迎流面為降壓增速流動，背流面為增壓減速流動.

圖2 定常渦流線圖(Re=30)

圖3 圓柱體表面無因次壓力沿壁面角分布(Re=30)

2.2 非定常渦

隨著Re增加，柱體背后繞流尾跡拉長，當Re超過臨界Re時，圓柱體兩側(cè)形成周期性、交替向下游運動的旋渦流.當Re=60時開始出現(xiàn)周期性渦脫.當Re=90時，典型旋渦生成—長大—脫落過程見圖4，1個周期內(nèi)圓柱體表面的無因次壓力分布曲線見圖5.

圖4 不同周期渦脫過程(Re=90)

圖5 典型渦脫周期內(nèi)的無因次壓力分布(Re=90)

由圖4可見，周期性渦脫是一個連續(xù)過程，一側(cè)旋渦在脫落的同時，另一側(cè)的旋渦已在孕育長大.T=297.34 s時圓柱體下側(cè)旋渦脫落使繞流改善，速度較快，靜壓力較低；圓柱體上側(cè)旋渦形成并長大，繞流較差，速度較慢，靜壓力較高.與T=297.34 s時的無因次壓力分布曲線對應(見圖5)，使周期性交變作用力方向始終指向剛釋放完渦流的一側(cè).由圖5可見，由于圓柱體后旋渦交替脫落作用，背流面壓力因數(shù)出現(xiàn)波動在某一時刻不再具有對稱性，但整個周期還是具有較好的對稱性.

流體繞流圓柱體時，作用在圓柱體上的周期性交變作用力由流體作用在圓柱體表面的摩擦力和壓力差產(chǎn)生，可以分解為沿流體運動方向的阻力和垂直于流體運動方向的升力.圓柱體后產(chǎn)生周期性渦脫后的升力因數(shù)Cl、阻力因數(shù)Cd的時程曲線及由頻譜分析獲得的振幅頻圖見圖6.由圖6可見,升力、阻力因數(shù)周期性變化，阻力因數(shù)的震蕩頻率是升力因數(shù)的2倍，升力因數(shù)振幅遠大于阻力因數(shù)的振幅.

圖6 升力、阻力因數(shù)時程及頻譜分析

2.3 參數(shù)變化規(guī)律

2.3.1 壓力因數(shù)

在不同的Re范圍時，圓柱體的時均壓力因數(shù)分布規(guī)律不一致(見圖7).最小壓力點位置隨Re的增加逆流前移，最小壓力因數(shù)的絕對值隨Re增加而增大.

圓柱體前、后駐點時均壓力因數(shù)變化曲線見圖8.由圖8可見,隨Re的增加，前駐點壓力因數(shù)Cps變??；定常渦時變化速率較快，產(chǎn)生周期性旋渦脫落后變化速率減緩，最終穩(wěn)定在1.0附近，與實驗測量數(shù)據(jù)吻合[8].

圖7 時均壓力因數(shù)分布曲線

后駐點壓力因數(shù)Cpb在臨界雷諾數(shù)之前隨Re增加而變大，當產(chǎn)生周期性渦脫時Cpb隨Re增加而減小.文中計算的Cpb偏小是由圓柱體后尾流的復雜性和不穩(wěn)定性引起的[9]，在30≤Re≤60時有拐點，說明Cpb在臨界雷諾數(shù)附近有極值，與文獻[12]結(jié)果一致.

2.3.2 阻力及升力因數(shù)

分析不同Re時形成周期性渦脫后的阻力、升力因數(shù)曲線頻譜，阻力因數(shù)Cd的基波頻率是升力因數(shù)Cl基波頻率的2倍(見圖9).

圖8 前、后駐點壓力因數(shù)隨Re數(shù)變化曲線

圖9 升力、阻力因數(shù)基頻關(guān)系

圖10 平均阻力因數(shù)變化規(guī)律

2.3.3St及邊界層分離角度

理論上求解旋渦脫落頻率困難，因此工程一般用斯特羅哈數(shù)St確定旋渦脫落的頻率.St為量綱一的參數(shù)，與旋渦發(fā)生體形狀及Re有關(guān)，圓柱體St與Re關(guān)系見圖12，模擬結(jié)果位于Lienhard研究獲得的廣闊Re時St范圍內(nèi)，并與實驗數(shù)據(jù)吻合較好[15-17]，變化規(guī)律與文獻[18]符合.

圖11 脈動力均方根變化規(guī)律

圖12 St隨Re變化規(guī)律

伴隨周期性渦脫的層流流動,St隨Re的增加而變大，當二維流動向三維流動過渡時(Re≈150)[9]，由于三維流動的不穩(wěn)定擾動致使尾渦結(jié)構(gòu)軸向不穩(wěn)定，使St降低；形成三維旋渦脫落以后，St變化趨于平緩，流動處于亞臨界狀態(tài)(Re<2×105)時，計算渦流釋放頻率的經(jīng)驗公式中，St為0.20～0.21.

隨著圓柱體兩側(cè)旋渦不間斷地生成、脫落，圓柱體壁面附近邊界層不斷分離；通過對相應Re時整個渦脫周期圓柱體壁面剪切應力的監(jiān)測，得到不同Re時圓柱體壁面邊界層分離的變化范圍.

由平均分離角度及角度波動隨Re變化曲線(見圖13)可知，邊界層分離點隨Re的增大逆流而上，與最小壓力點移動變化趨勢一致，平均分離角度θs不斷增加，但變化速率不斷減小，與WU M H實驗測量數(shù)據(jù)及Park J(Re<160)計算結(jié)果一致[9].由于阻塞比的差異，Grove測量數(shù)據(jù)與模擬結(jié)果有一定差距；當Re>150時，圓柱體后尾流的三維不穩(wěn)定性導致實驗測量的平均分離角度出現(xiàn)波動.

分離角度波動θf隨Re的增加而變大(見圖13(b))，這主要是由旋渦脫落的振蕩特性以及隨Re的增加旋渦脫落的增強導致的.

圖13 平均分離角度及分離角度波動隨Re變化曲線

3 結(jié)論

(1)對于低雷諾數(shù)圓柱繞流,隨Re的增加，流動狀態(tài)由定常渦轉(zhuǎn)變?yōu)橹芷谛孕郎u脫落狀態(tài)，通過頻譜分析獲得周期性變化的阻力、升力因數(shù)的基波頻率為2倍的線性關(guān)系.

(2)周期性渦脫導致壁面壓力因數(shù)不再對稱，隨著Re的增加，最低壓力點逆流前移，邊界層分離區(qū)域變大，時均分離點隨最低壓力點一起前移.

(3)脈動阻力、升力因數(shù)隨Re的增加而變大，平均分離角度及角度波動隨之增大.