沈立龍，劉明維，吳林鍵，李鵬飛，舒 丹

（重慶交通大學(xué)國家內(nèi)河航道整治工程技術(shù)研究中心，水利水運工程教育部重點實驗室，重慶 400074；重慶交通大學(xué)河海學(xué)院，重慶 400074）

亞臨界雷諾數(shù)下圓柱和方柱繞流數(shù)值模擬

沈立龍，劉明維，吳林鍵，李鵬飛，舒 丹

（重慶交通大學(xué)國家內(nèi)河航道整治工程技術(shù)研究中心，水利水運工程教育部重點實驗室，重慶 400074；重慶交通大學(xué)河海學(xué)院，重慶 400074）

基于RNG k?ε模型，采用有限體積法對亞臨界雷諾數(shù)條件下（Re=3×102～3×105）的二維單圓柱和單方柱繞流進行數(shù)值模擬與仿真。得到了圓柱和方柱繞流阻力系數(shù)Cd與Strouhal數(shù)隨雷諾數(shù)的變化規(guī)律。同時對雷諾數(shù)Re=22 000下的圓柱、方柱繞流相關(guān)特性進行詳細對比分析。計算結(jié)果表明：同一雷諾數(shù)下，單圓柱繞流阻力系數(shù)Cd較單方柱低，但圓柱的Strouhal數(shù)較單方柱則要高。雖然二者邊界層分離點不同，但流場的演變與漩渦的脫落具有一定相似特性。

亞臨界雷諾數(shù)；圓柱和方柱；繞流；數(shù)值模擬；RNGk?ε模型

Biography：SHEN Li?long（1990-）,male,master student.

鈍體繞流現(xiàn)象普遍存在于實際生活中，如橋墩繞流，深海輸油立管繞流，飛機機翼繞流，化學(xué)反應(yīng)塔繞流等。繞過鈍體的流體流動極不穩(wěn)定，周圍的流場也變得，當渦街發(fā)生時，交替產(chǎn)生并脫落的漩渦誘發(fā)交變力，使得被繞流的物體產(chǎn)生振動及噪聲，嚴重時產(chǎn)生共振及聲振，對建筑物正常使用安全造成了極大的威脅。由此，鈍體繞流問題現(xiàn)已成為流體力學(xué)領(lǐng)域所研究熱點問題之一。國內(nèi)外學(xué)者對此進行了大量的模擬研究，Standsby等人［1］采用隨機渦法對放置方法不同的串列、并排雙圓柱進行了數(shù)值分析，得到結(jié)果與實驗基本相符；蘇明德等人［2］采用二階精度有限體積法和Smagorinsky渦粘性模式對雷諾數(shù)Re為200和20 000下的圓柱繞流進行了大渦模擬，同樣得到了與實驗結(jié)果相近的升、阻力系數(shù)等動力學(xué)參數(shù)；Kawai Hiromasa等人［3］采用差分法對雷諾數(shù)Re為200的串列方柱的繞流做了數(shù)值模擬；王遠成等人［4］采用湍流RNGk?ε模型對方柱繞流流場的相關(guān)特性進行了數(shù)值模擬，很好的反映出了漩渦脫落的尾跡結(jié)構(gòu)。

樁柱的繞流受到邊界條件、湍流強度、柱面粗糙度等因素的影響［5］，但起決定作用的是雷諾數(shù)。隨著雷諾數(shù)的增加，柱后的流動相應(yīng)由定常逐漸失穩(wěn)，慢慢變得越發(fā)復(fù)雜，最后到達超臨界區(qū)的湍流分離。對于300≤Re≤3×105的亞臨界區(qū)，柱面邊界層為層流，而尾流卻已變?yōu)橥牧鳒u街［6］。天然河道等實際情況中，雷諾數(shù)相對較大，通過實驗研究發(fā)現(xiàn)，亞臨界區(qū)樁柱的繞流情況已基本能反映出較高雷諾數(shù)下天然河道的繞流特性。

本文運用CFD軟件、基于RNGk?ε模型對亞臨界雷諾數(shù)下單圓柱和單方柱繞流進行了數(shù)值模擬，計算得到了這兩種情況下的阻力系數(shù)、Strouhal數(shù)等重要動力學(xué)參數(shù)，同時模擬出了樁柱繞流過程中的漩渦脫落、流場演變等情況，為樁柱繞流問題的實際物理機制進行了合理的仿真分析，在實際工程中具有良好的借鑒意義。

1 控制方程

不可壓縮粘性流體流體的運動可用Navier?Stokes方程［7］來描述，連續(xù)性方程與動量方程分別為

式中：ui為速度分量；p為壓力；ρ為流體的密度；ν為流體的動力黏性系數(shù)。

對于湍流情況，本文采用RNGk?ε模型，RNGk?ε模型是k?ε模型的改進方案，由Yakhot和Orzag［8］提出，在RNGk?ε模型中，通過在大尺度運動和修正后的粘度項體現(xiàn)小尺度的影響，而使這些小尺度運動有系統(tǒng)地從控制方程中去除。所得到的k方程和ε方程，與標準k?ε模型非常相似［9］，其表達式如下

式中：Gk為由于平均速度梯度引起的湍動能k的產(chǎn)生項，μeff=μ+μt，μt=ρCμk2/ε；經(jīng)驗常數(shù)Cμ=0.084 5，αk=αε=1.39，C2ε=1.68。

相對于標準k?ε模型，RNG k?ε模型通過修正湍動粘度，考慮了平均流動中的旋轉(zhuǎn)及旋轉(zhuǎn)流動情況，同時在ε方程中加入了一項Rε，從而反映了主流的時均應(yīng)變率Eij，這樣該模型中產(chǎn)生項不僅與流動情況有關(guān)，而且在同一問題中也還是空間坐標的函數(shù)。從而，RNGk?ε模型可以更好的處理高應(yīng)變率及流線彎曲程度較大的流動［10］。

2 數(shù)值算法

2.1 有限元建模

建模過程中，運用大型通用流體力學(xué)計算軟件Fluent輔助計算。計算區(qū)域采用分塊結(jié)構(gòu)化網(wǎng)格，柱體表面網(wǎng)格做加密處理，邊界區(qū)網(wǎng)格相對稀疏，如圖1所示，計算區(qū)域為33D×15D（D為圓柱直徑和方柱邊長），模型上游來流區(qū)域為7D，下游尾流區(qū)域為25D，上下邊界取7D。流體自左向右流動。

圖1 單圓柱、單方柱計算區(qū)域及柱面網(wǎng)格Fig.1 Computational domain and surface grid of a circular cylinder and a square cylinder

2.2 邊界條件

入流邊界條件：u=V，v=0。

出流邊界條件：自由出流。

壁面條件：上下邊界及柱體表面為無滑移固壁邊界。

2.3 計算模型參數(shù)

計算模型采用RNGk?ε模型，壓力速度耦合方式為PISO算法，動量、湍動能和湍動耗散率均采用二階迎風(fēng)格式。其中，雷諾數(shù)、來流速度等參數(shù)見表1所示。

表1 計算模型參數(shù)Tab.1 Computational model parameters

2.4 模型驗證

由邊界層理論可知，液體對所繞流物體的阻力包括固體表面摩擦阻力和由壓強分布不均造成的壓強阻力。同樣被繞流物體阻力系數(shù)Cd也來源于壁面粘性摩擦與表面壓力系數(shù)，Strouhal數(shù)是一個反應(yīng)漩渦脫落的頻率相對于來流速度快慢的重要無量綱參數(shù)，二者定義為

式中：Fd為繞流阻力，D為圓柱直徑或方柱邊長；u為來流速度；f為渦街脫落頻率；ρ為流體密度。

為驗證本文中數(shù)值模擬的正確性，計算得到了亞臨界雷諾數(shù)下單圓柱及單方柱在不同雷諾數(shù)下的繞流阻力系數(shù)Cd及Strouhal值，并將其與文獻中的計算結(jié)果相對比，結(jié)果如表2、表3所示。其中，Cd代表柱體阻力系數(shù)平均值。

通過表2、表3中數(shù)據(jù)發(fā)現(xiàn)，對亞臨界雷諾數(shù)下圓柱、方柱繞流模擬計算得出的結(jié)果能很好的與參考文獻中數(shù)據(jù)相吻合，即驗證了本文計算結(jié)果的正確性。

同時，根據(jù)表2、表3中的數(shù)據(jù)繪制亞臨界雷諾數(shù)下圓柱和方柱的繞流阻力系數(shù)Cd及Strouhal值隨雷諾數(shù)的變化曲線，如圖2所示。

表2 單圓柱繞流模擬計算結(jié)果Tab.2 Simulation results of flow around a circular cylinder

表3 單方柱繞流模擬計算結(jié)果Tab.3 Simulation results of flow around a square cylinder

圖2 阻力系數(shù)與Strouhal數(shù)隨雷諾數(shù)走勢Fig.2 Trends of drag coefficient and Strouhal number with the change of Reynolds numbers

從圖2中可以看出亞臨界雷諾數(shù)下圓柱和方柱繞流阻力系數(shù)Cd及Strouhal數(shù)值曲線都相對平穩(wěn)，其值變化不大，圓柱和方柱的Strouhal數(shù)分別在0.2和0.13上下微小波動。同一雷諾數(shù)下圓柱繞流阻力系數(shù)較方柱要小，但圓柱的Strouhal數(shù)卻略大于方柱，說明圓柱繞流阻力弱于方柱，漩渦脫落頻率相對于來流速度卻更快。

3 結(jié)果分析

由于文章篇幅有限，結(jié)果數(shù)據(jù)較多。在此，結(jié)合大量文獻資料及實際計算數(shù)據(jù)，只選擇當雷諾數(shù)Re=22 000時的典型情況進行著重分析討論。

3.1 繞流阻力系數(shù)及升力系數(shù)

當雷諾數(shù)Re=22 000時，壓強阻力顯著大于摩擦阻力，占主導(dǎo)地位，此時的阻力系數(shù)主要來源于柱體表面的壓力系數(shù)。繞流升力是由垂直于來流方向柱體兩側(cè)壓力不均造成。由圖3可以看出，圓柱、方柱繞流阻力系數(shù)Cd、升力系數(shù)Cl都呈周期性單一頻率振動，但方柱的Cd、Cl迭代周期明顯大于圓柱。方柱升力系數(shù)Cl幅值較圓柱更大，說明脈動較強。

3.2 流場演變

如圖4、圖5所示，為單圓柱與單方柱在雷諾數(shù)Re=22 000時一個漩渦脫落周期內(nèi)5個典型時刻的流線圖及渦量圖。

圖3 單圓柱與單方柱繞流阻力、升力系數(shù)Fig.3 Lift and drag coefficients of a circular cylinder and a square cylinder

圖4 單圓柱與單方柱繞流漩渦脫落周期內(nèi)不同時刻流線圖Fig.4 Streamlines of the flow around a circular cylinder and a square cylinder at different time in a vortex shedding cycle

圖5 單圓柱與單方柱繞流漩渦脫落周期內(nèi)不同時刻渦量等值線圖Fig.5 Vorticity contours of the flow around a circular cylinder and a square cylinder at different time in a vortex shedding cycle

圖4、圖5為一個周期內(nèi)5個典型時刻的流線圖和渦量圖，由圖可以看出0和T時刻的流線圖、渦量圖基本相同，且二者對稱于T/2時刻圖像，T/4時刻圖像同樣對稱于3T/4時刻圖像。此5個典型時刻圖像生動的表明出了單圓柱、單方柱繞流流場運動及漩渦脫落呈周期性變化的特性。在0時刻，柱后偏上位置產(chǎn)生一個大漩渦，致使柱的繞流升力達到最大，新生的漩渦隨著水流向下方遷移，至T/4時刻，到達柱的正后方，此時升力為零。由于液體的粘滯性，從柱上方遷移下來的漩渦在運動過程中能量逐漸衰減，最后消失，漩渦完成脫落。T/2時刻柱后偏下位置產(chǎn)生漩渦，此時升力同樣最大，漩渦的脫落等同于柱后上方漩渦的脫落特性。一個周期完成，下個周期開始，由此在柱后形成兩列交錯且周期性擺動的漩渦，即Karman渦街，圓柱、方柱都是如此，說明二者繞流流場運動及漩渦脫落具有相似特性，只是兩者柱體表面旋渦脫落周期大小不同，圓柱較方柱偏小，即同一雷諾數(shù)下圓柱繞流旋渦脫落頻率大于方柱繞流旋渦脫落頻率。

3.3 邊界層分離

亞臨界雷諾數(shù)下圓柱、方柱繞流邊界層分離同為層流分離，且邊界層內(nèi)液體運動伴有能量損失，但通過模擬發(fā)現(xiàn)二者分離點有很大區(qū)別，對于圓柱，隨著雷諾數(shù)的增大分離點逐漸向后推移。液流與柱面分離原理（圖6）：在壓強最大的P點至柱上下A、B兩點的邊界層內(nèi)液流處于加速減壓狀態(tài)，壓能除部分轉(zhuǎn)化為動能外主要補償能量損失，在A、B兩點時壓強最小，流速最大，兩點之后液流處于減速增壓狀態(tài)，當動能減小至零而無法提供給壓能時，完成分離。

而對于方柱，由于柱體有凸出的棱角，繞流情況就變的不如圓柱繞流那樣復(fù)雜，分離點固定在方柱的前兩個棱角位置，如圖6中C、D兩點。

圖6 圓柱和方柱繞流示意圖Fig.6 Sketch of flow around circular and square cylinder

4 結(jié)論

運用CFD軟件，通過對亞臨界雷諾數(shù)下的單圓柱和單方柱繞流進行數(shù)值模擬與對比，本文得出以下結(jié)論：

（1）柱體繞流模擬中，網(wǎng)格的精度與劃分方式對計算結(jié)果有很大影響，為了得到更精確的結(jié)果同時又節(jié)省計算時間，要選取數(shù)量適當?shù)木W(wǎng)格及合適的網(wǎng)格劃分方式進行計算，同時對柱體表面及柱后方軸線區(qū)域網(wǎng)格做加密處理。

（2）亞臨界雷諾數(shù)下，圓柱和方柱繞流阻力系數(shù)都變化不大，Strouhal數(shù)更是趨向于定值，阻力系數(shù)及升力系數(shù)都呈周期性單一頻率振動。同一雷諾數(shù)下圓柱繞流阻力系數(shù)較方柱繞流阻力系數(shù)要小，但圓柱的Strouhal數(shù)較方柱的要高。

（3）同一雷諾數(shù)下，單圓柱和單方柱繞流生成的旋渦脫落周期不同，方柱的偏大，但二者漩渦脫落特征與流場運動具有相似特性。

（4）亞臨界雷諾數(shù)下，圓柱和方柱繞流邊界層分離點不同，圓柱的分離點隨著雷諾數(shù)的增大而逐漸向柱后方推移，方柱的分離點則固定在柱前兩個棱角位置。

［1］Slaoutiand A,Stansby P K.Flow around two circular cylinders by the random?vortex method［J］.Journal of Fluids and Structures,1992,6：641-670.

［2］蘇銘德,康欽軍.亞臨界雷諾數(shù)下圓柱繞流的大渦模擬［J］.力學(xué)學(xué)報,1999,31（1）：100-105.

SU M D,KANG Q J.Large eddy simulation of the turbulent flow around acircular cylinder at subcritical Reynolds numbers［J］.Acta Mechanica Sinica,1999,31（1）：100-105.

［3］Kawai,Hiromasa,Fujinamic,et al.Numerical simulation of flow around square prisms in tandem arrangements［C］//Gulf.9th IC?ME.Houston Tx：Gulf Publishing,1995：185-186.

［4］王遠成,吳文權(quán).方柱繞流流場的RNG方法模擬研究［J］.水動力學(xué)研究與進展,2004（19）：916-920.

WANG Y C,WU W Q.Numerical simulation of flow around square cylinder using RNGk?εturbulence model［J］.Journal Hydro?dynamics,2004（19）：916-920.

［5］顧志福,孫天風(fēng),林榮生.高雷諾數(shù)時串列雙圓柱平均壓力的實驗研究［J］.空氣動力學(xué)學(xué)報,1997,15（3）：393-399.

GU Z F,SUN T F,LIN R S.Experimental research of average pressure of two circular cylinder in tandem arrangement at high Reynolds number［J］.Acta Aerodynamica Sinica,1992,15（3）：393-399.

［6］楊爍,吳寶山.二維圓柱繞流數(shù)值模擬［J］.中國造船,2007（48）：533-540.

YANG S,WU B S.The numerical simulation of two dimensional circular flow［J］.China Shipbuilding,2007（48）：534-540.

［7］張遠君.流體力學(xué)大全［M］.北京：北京航空航天大學(xué)出版社,1991.

［8］Yakhot V,Orzag S A.Renormalization grop analysis of turbulence［J］.Basic theory.J Scient Comput,1986,1：3-11.

［9］Versteeg H K,Malalasekera W.An Introduction to Computational Fluid Dynamics［M］.Wiley：The Finite Volume Method,1995.

［10］王福軍.計算流體動力學(xué)分析［M］.北京：清華大學(xué)出版社,2004.

［11］史里希廷H.邊界層理論［M］.北京：科學(xué)出版社,1988.

［12］Lyn D A,Einav S,Rodi W.A Laser Doppler Velocimetry Study of Ensemble?averaged Chsrscteristic of the Turbulent Flow Near of a Square Cylinder［J］.Fluid Mechanics,1995,304：285-319.

［13］Lyn D A,Rodi W.The flapping shear layer formed by flow separation from the forward corner of a square cylinder［J］.Fluid Mesh,1994,267：353-376.

Numerical simulation of the flow around circular and square cylinder at sub?critical Reynolds numbers

SHEN Li?long，LIU Ming?wei，WU Lin?jian，LI Peng?fei，SHU Dan
（National Inland Waterway Regulation Engineering Research Center,Key Laboratory of Hydraulic&Waterway Engineering of the Ministry of Education,Chongqing Jiaotong University,Chongqing400074,China；School of River&Ocean Engineering，Chongqing Jiaotong University，Chongqing400074,China）

With theRNG k?εmodel and the finite volume method，the flow around a two?dimensional circular cylinder and a two?dimensional square cylinder at sub?critical Reynolds numbers was simulated and evaluated,and the trends of drag coefficients and Strouhal numbers with the change of Reynolds numbers from 3×102to 3×105were obtained.Meanwhile,characteristics of the flow around a circular cylinder and a square cylinder atRe=22 000 have been compared and analyzed in detail.Calculation results show that drag coefficient of flow around a circular cylinder is lower than the drag coefficient of flow around a square cylinder with the same Reynold number.Howev?er,the Strouhal number of flow around a circular cylinder is slightly higher than the Strouhal number of flow around a square cylinder.Though their boundary layer separation points are different,the flow field evolution and the vor?tex shedding of two cylinders have some similarities.

sub?critical Reynolds numbers；circular and square cylinder；flow around body；numerical simula?tion；RNG k?εmodel

TV 143+.2；O 351

A

1005-8443（2014）03-0227-07

2013-11-28；

2014-01-06

國家科技支撐計劃項目（2012BAB05B04）

沈立龍（1990-），男，山東省菏澤市人，碩士研究生，主要從事港口海岸工程結(jié)構(gòu)與基礎(chǔ)方面的研究。