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

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

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

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

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

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

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

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

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

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

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

1 控制方程

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

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

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

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

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

2 數(shù)值算法

2.1 有限元建模

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

圖1 單圓柱、單方柱計(jì)算區(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 計(jì)算模型參數(shù)

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

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

2.4 模型驗(yàn)證

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

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

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

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

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

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

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

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

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

3 結(jié)果分析

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

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

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

3.2 流場(chǎng)演變

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

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

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

3.3 邊界層分離

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

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

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

4 結(jié)論

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

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

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

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

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

［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é)報(bào),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］王遠(yuǎn)成,吳文權(quán).方柱繞流流場(chǎng)的RNG方法模擬研究［J］.水動(dòng)力學(xué)研究與進(jìn)展,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ù)時(shí)串列雙圓柱平均壓力的實(shí)驗(yàn)研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào),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］.中國(guó)造船,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］張遠(yuǎn)君.流體力學(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］王福軍.計(jì)算流體動(dòng)力學(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

國(guó)家科技支撐計(jì)劃項(xiàng)目（2012BAB05B04）

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