朱潤野，刁明軍，劉 偉

(1.四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610064;2.四川水利干部學(xué)校，四川 成都 610072)

低雷諾數(shù)下五圓柱繞流的數(shù)值模擬研究

朱潤野1，刁明軍1，劉 偉2

(1.四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610064;2.四川水利干部學(xué)校，四川 成都 610072)

基于四叉樹自適應(yīng)網(wǎng)格，采用有限體積法求解二維不可壓粘性N-S方程，對X型排列和十字型排列的五圓柱繞流進(jìn)行了模擬，計(jì)算均在Re=100的條件下進(jìn)行.在文中給出了升力系數(shù)及阻力系數(shù)隨圓柱間距比和排列方式變化的規(guī)律，并結(jié)合不同間距比下的流動(dòng)特征進(jìn)行了分析.計(jì)算結(jié)果表明，在不同的間距比下，五圓柱繞流存在三種不同的流動(dòng)形式;在同一間距比下，十字型布置和X型布置的流動(dòng)形式則既有類似之處又有各自的特點(diǎn).在不同排列方式下，升阻力系數(shù)隨圓柱間距比的變化規(guī)律有較大不同.

五圓柱繞流;四叉樹網(wǎng)格;數(shù)值模擬;升力系數(shù);阻力系數(shù)

在工程應(yīng)用中，多圓柱繞流現(xiàn)象十分常見，如海洋石油平臺樁柱、斜拉橋拉索、海底輸油管線、電廠冷卻塔群、桅桿、天線等.同時(shí)，因?yàn)楫?dāng)流體流過多圓柱時(shí)，由于各個(gè)圓柱之間存在相互影響，多圓柱繞流往往會(huì)出現(xiàn)復(fù)雜且多樣的流動(dòng)形態(tài)，與單圓柱繞流十分不同，故針對多圓柱繞流的研究有著廣泛的實(shí)際意義.

前人對多圓柱繞流的研究往往集中在雙圓柱、三圓柱和四圓柱繞流.對雙圓柱繞流問題，Willianmson［1］通過實(shí)驗(yàn)，研究了并列雙圓柱后的尾流演變;鄭庭輝等［1］采用有限體積法對，對雷諾數(shù)為100～300的串列雙圓柱的尾流進(jìn)行了研究.顧志福和孫天風(fēng)［3］采用實(shí)驗(yàn)的方法，對雷諾數(shù)為5.5×104和1.4×104的三圓柱繞流問題進(jìn)行了研究.Yan Bao等［4］采用二階CBS有限元算法對呈正三角形排列的圓柱繞流問題進(jìn)行了研究，并進(jìn)行了有益的討論.Lam［5］通過實(shí)驗(yàn)，研究了亞臨界雷諾數(shù)下，四圓柱繞流的流動(dòng)情況，同時(shí)采用有限體積法模擬了三維情況下的四圓柱繞流問題［6－12］.

國內(nèi)外對五圓柱繞流問題的研究非常少見，考慮到五圓柱繞流問題在實(shí)際工程中的潛在價(jià)值，對其的研究具有一定的意義.本文利用數(shù)值方法，對低雷諾數(shù)下，兩種不同排列方式的五圓柱繞流問題進(jìn)行了研究，給出了在不同間距比下一些重要的流動(dòng)參數(shù)，并對流態(tài)進(jìn)行了分析.

1 數(shù)值模型

1.1 控制方程

二維不可壓粘性非定常流動(dòng)運(yùn)動(dòng)方程為

1.2 數(shù)值方法

采用分步預(yù)估法離散時(shí)間項(xiàng)，控制方程可寫作:

采用二階迎風(fēng)格式離散對流項(xiàng).該方法在時(shí)間及空間上均具有二階精度.具體細(xì)節(jié)參考文獻(xiàn)［7］［8］.

1.3 計(jì)算區(qū)域及邊界條件

本文討論的五圓柱的布置情況分為X型分布和十字型分布，如圖1(a)和圖1(b)0～4分別代表五個(gè)柱的編號，下文以C0～C4指代，其中C0為位于中間的圓柱.L代表兩個(gè)圓柱之間的距離.五個(gè)圓柱具有相同的直徑，用D表示圓柱的直徑.

圖1(a)X型布置Fig.1(a)Five cylinders in X arrangement

圖1(b) 十字型布置Fig1.(b)Five cylinders in cruciform arrangement

計(jì)算區(qū)域及邊界條件如圖1(c)所示.為了盡可能減少邊界的影響，在y方向上寬度為50D，在x方向上長度為70D.在圓柱表面采用無滑移邊界條件.

圖1(c) 計(jì)算區(qū)域及邊界條件Fig.1(c)Computational domain and boundary condition

本文采用四叉樹自適應(yīng)網(wǎng)格進(jìn)行計(jì)算，網(wǎng)格在計(jì)算過程中將根據(jù)流速梯度和壓力變化自動(dòng)加密.圖2(a)為計(jì)算完成后單個(gè)圓柱周圍的網(wǎng)格;圖2(b)為計(jì)算完成后計(jì)算域的網(wǎng)格.

圖2(a) 單個(gè)圓柱周圍網(wǎng)格Fig.2(a)Mesh around single cylinder

圖2(b) 計(jì)算域網(wǎng)格Fig.2(b)Mesh of computational domain

1.4 升、阻力系數(shù)的計(jì)算

用Cd和Ci表示圓柱的升、阻力.每個(gè)圓柱的升阻力以下標(biāo)表示.根據(jù)文獻(xiàn)［5］隨時(shí)間變化的升阻力定義為:

阻力系數(shù)和升力系數(shù)的時(shí)均值和定義為:

阻力系數(shù)和升力系數(shù)的均方根值和定義為:

升力系數(shù)的峰值定義為:

2 數(shù)值有效性驗(yàn)證

為了驗(yàn)證計(jì)算的有效性，同時(shí)為了與五圓柱繞流進(jìn)行比較，本文計(jì)算了在Re=100時(shí)單圓柱的繞流問題，計(jì)算結(jié)果見表1.從表1可知，本文計(jì)算與前人結(jié)果的最大誤差為3.7%，誤差較小.

3 流動(dòng)形式

本文計(jì)算了 L/D=1.0，1.5，2.0，2.5，3.0，4.0，5.0等幾種間距比情況下的繞流問題，得到了各間距比下的流動(dòng)形式，以下將根據(jù)兩種布置分別討論流動(dòng)形式隨間距比的變化并比較它們之間的異同.

表1 單圓柱繞流計(jì)算結(jié)果對比Tab.1 Comparisons of mean drag，amplitude of lift force and Strouhal number at Re=100 for single cylinder

3.1 X型布置

在小間距比下，X型布置的流動(dòng)形式類似于單圓柱繞流.由于圓柱之間間隙流的影響，圓柱C2的尾渦偏向上方，圓柱C3的尾渦偏向下方.

圖3 L/D=1時(shí)渦量等值線圖Fig.3 Snapshot of instantaneous vorticity contours for L/D=1

圖4和圖5為五圓柱處于中等間距比時(shí)，渦量等值線圖和流線圖.從圖4可以看到，C1和C4后的渦呈對稱反相脫落.脫落的渦相對于來流方向?qū)ΨQ分布.從圖5可以看到，由于圓柱C1和C4的影響，中間圓柱C0的尾流受到擠壓，并在C0后形成了穩(wěn)定且對稱的回流區(qū).根據(jù)文獻(xiàn)［13］，當(dāng)L/D=2時(shí)，上游圓柱不會(huì)產(chǎn)生渦脫落，其剪切層將附著在下游圓柱上，但在本文的情況中，由于中間圓柱與C2、C3之間的間隙流，C2與C3后的尾流將會(huì)偏向柱群外側(cè)，只有C2下邊緣的剪切層和C3上邊緣的剪切層會(huì)附著在下游圓柱上.在遠(yuǎn)離柱群的下游，上下兩側(cè)的渦會(huì)相互靠近融合形成單個(gè)的渦.

圖4 L/D=2時(shí)渦量等值線圖Fig.4 Snapshot of instantaneous vorticity contours for L/D=2

圖5L/D=2時(shí)流線圖Fig.5 Snapshot of streamlines for L/D=2

圖6為L/D=4時(shí)的渦量等值線圖.當(dāng)間距比較大時(shí)，圓柱之間的相互影響減小，在C2、C3和C0后形成了充分發(fā)展的渦脫落.圓柱C2和C3的渦呈反相脫落，并會(huì)撞擊下游圓柱.下游圓柱C1和C4的渦脫落依然會(huì)受到中間圓柱C0的影響而變得不規(guī)律.C1下緣和C4上緣脫落的渦會(huì)偏向中間.C1和C4的脫落的渦在柱群后會(huì)形成兩條渦帶.

圖6 L/D=4時(shí)渦量等值線圖Fig.6 Snapshot of instantaneous vorticity contours for L/D=4

圖7為L/D=1.5時(shí)的流線圖.C1和C4后形成了一大一小兩個(gè)渦，是因?yàn)镃1和C4兩圓柱間的間隙流偏向C1，從而使C1后的渦受到影響，其尾渦的長度明顯短于圓柱C4后的渦.圖7(a)和圖7(b)為無量綱時(shí)間分別為200和225時(shí)的流線圖，可以看到當(dāng)兩圓柱間的間隙流偏向某一圓柱后，間隙流的偏向就不會(huì)再發(fā)生改變，形成所謂的單穩(wěn)態(tài)流動(dòng).

圖7(a)Ut/D=200時(shí)流線圖Fig.7(a)Snapshot of streamlineswhen Ut/D=200

圖7(b)Ut/D=225時(shí)流線圖Fig.7(b)Snapshot of streamlineswhen Ut/D=225

3.2 十字型布置

當(dāng)間距比較小時(shí)，與X型布置類似，十字型布置的流動(dòng)形態(tài)類似于單個(gè)圓柱.渦會(huì)從C1、C3和C4后產(chǎn)生，并在柱群后一定距離內(nèi)融合形成一個(gè)較大的渦.從圖8可以看到，由于圓柱間的距離較近，C2和C0后的渦并沒有脫落，而是粘附在圓柱C4上;圓柱C4由于間隙流和周圍圓柱渦的影響，分離點(diǎn)前移.

隨著間距比的增大，圓柱C1和C3的渦開始脫落，圖9為L/D=2時(shí)渦量等值線圖.C1和C3的渦脫落會(huì)影響到C4后的流動(dòng)，使得C4后的流動(dòng)較混亂.C4后的渦會(huì)與C1和C3后的渦融合.C1和C3后脫落的渦隨著流動(dòng)，很快就會(huì)破碎，失去完整的形狀.

圖8 L/D=1時(shí)渦量等值線圖Fig.8 Snapshot of instantaneous vorticity contours for L/D=1

圖9 L/D=2時(shí)渦量等值線圖Fig.9 Snapshot of instantaneous vorticity contours for L/D=2

當(dāng)間距比較大時(shí)，受其他圓柱影響較小的圓柱會(huì)形成充分發(fā)展的渦脫落.圖10為L/D=4時(shí)的渦量等值線圖.可以看到，圓柱C1和C3的渦脫落受中間圓柱的影響減小，并呈反相脫落.圓柱C2的渦開始脫落，并影響到中間圓柱C0.中間圓柱C0后的渦撞擊到圓柱C4上，會(huì)與圓柱C4的渦融合并使C4的渦脫落更不規(guī)則.

4 間距比對升、阻力的影響

圖7為X型布置下升、阻力系數(shù)隨間距比的變化規(guī)律.圖8為十字型布置下升、阻力系數(shù)隨間距比的變化規(guī)律.根據(jù)文獻(xiàn)［13－15］，多圓柱繞流問題中，圓柱間的相互影響可分為三類:鄰近影響、剪切層影響和尾流影響.以下將基于這三種影響因素，對兩種布置下五圓柱繞流問題中，升、阻力系數(shù)的變化規(guī)律進(jìn)行分析.

4.1 X型布置

從圖11(a)中可以看到，C1和C4的時(shí)均阻力系數(shù)值小于C2和C3的時(shí)均阻力系數(shù)值，這是因?yàn)楫?dāng)間距比較小時(shí)，C1和C4位于上游圓柱剪切層形成的低壓區(qū)中，圓柱前后緣的壓差較小，此時(shí)和較小，隨著間距比增大，剪切層的影響變小，阻力系數(shù)便開始上升.同時(shí)來流主要作用在C2和C3上，故C2和C3的時(shí)均阻力系數(shù)大于C1和C4.圓柱C0由于在小間距比的情況下，部分受到圓柱C2和C3剪切層的影響，故其時(shí)均阻力系數(shù)較小但大于C1和C4的時(shí)均阻力系數(shù)和在L/D＞3后有一定下降，這與雙圓柱繞流的變化規(guī)律一致［16］.

從圖11(b)可以看到由于中間圓柱的存在，除去C0外的所有圓柱的時(shí)均升力系數(shù)不再等于0.而由于流態(tài)的對稱性，則幾乎為0.在間距比2＜L/D＜4時(shí)，作用于C2和C3的升力大于作用于C1和C4上的升力.

從圖11(c)可以看到，當(dāng)間距比大于1.5時(shí)，和隨著間距比增大而增大，但當(dāng)間距比大于4后則開始下降.這是因?yàn)殡S著間距比的增大，上游圓柱C2和C3后的渦開始脫落，同時(shí)C0的渦也開始脫落，脫落的渦撞擊到C1和C4上，造成和較大.圓柱C0也會(huì)受到C2和C3的影響，故變化規(guī)律與C1和C4類似，但大小小于和

從圖11(d)可以看到，當(dāng)間距比1.5≤L/D≤4時(shí)，由于C1和C4后的渦脫落和上游圓柱脫落的渦的影響和隨著間距比增大而增大和也會(huì)在L/D大于某個(gè)值后開始增加，并逐漸接近單個(gè)圓柱的升力系數(shù)均方根值.

圖11 X型布置下升、阻力系數(shù)變化規(guī)律Fig.11 variation of force coefficients with gap spacing at X arrangement

4.2 十字型布置

圖12(a)為時(shí)均阻力系數(shù)隨間距比的變化規(guī)律.隨著間距比的增加而增加.當(dāng)L/D小于1.5時(shí)，由于C1和C3鄰近影響，C1、C3和C0之間的間隙流偏向圓柱C0，在C0前間形成了一個(gè)高壓區(qū)，使得大于，但隨著L/D的增大，鄰近影響逐漸減小，使得隨著L/D的增大而減小，當(dāng)L/D≥3時(shí)又開始增大，并接近于單個(gè)圓柱繞流時(shí)的時(shí)均阻力系數(shù)大小.當(dāng)L/D≤3時(shí)接近于0.

圖12(b)為時(shí)均升力系數(shù)隨間距比變化的曲線.和的絕對值隨著L/D的增大而減小.而和則幾乎接近于0.

圖12(c)為阻力系數(shù)的均方根值隨D/L的變化規(guī)律.隨著L/D的增加突然上升，達(dá)到0.174，并且在L/D=2～L/D=4的范圍內(nèi)一直較大，這主要是因?yàn)樵贚/D較小時(shí)，C1和C3后脫落的渦影響了C4周圍的壓力分布，而隨著L/D增大，雖然C1和C3的影響減小，但C0后的渦開始脫落，從而持續(xù)影響C4周圍的壓力分布.在L/D=3的位置，C1和C3的影響減小，而C0后的尾流還未充分發(fā)展，因此有小幅度的下降.

圖12(d)為升力系數(shù)的均方根值隨L/D的變化規(guī)律.和隨這L/D的增加而增加，并趨向于單個(gè)圓柱的升力系數(shù)均方根值大小.當(dāng)L/D＞3，和突然開始增大，且的值大于，因?yàn)镃4同時(shí)受到C2和C0脫落的渦的影響，而C0只受到C2脫落的渦的影響.

圖12 十字型布置下升、阻力系數(shù)變化規(guī)律Fig.12 variation of force coefficients with gap spacing at cruciform arrangement

5 結(jié)論

本文采用在時(shí)間和空間上皆具有二階精度的分步預(yù)估法，基于四叉樹自適應(yīng)網(wǎng)格，對Re=100條件下X型布置和十字型布置的五圓柱繞流問題進(jìn)行了計(jì)算，得到了如下結(jié)論:

1)低雷諾數(shù)下，布置相同的五圓柱，根據(jù)L/D的大小可以將流動(dòng)形式分為三種類型.當(dāng)L/D較小時(shí)，流動(dòng)形式類似于單個(gè)圓柱的繞流，在五圓柱后形成單個(gè)脫落的渦.在中等間距下，自柱群上下緣產(chǎn)生渦脫落，而處于中間位置的圓柱由于受到其他圓柱的影響尾流無法充分發(fā)展成脫落渦.在大間距比下，柱群中不受上下游影響的圓柱的尾流會(huì)充分發(fā)展，形成渦脫落.

2)在同一間距比下，兩種布置形式的流動(dòng)形態(tài)會(huì)有所不同.小間距比下，十字型布置后渦開始脫落的距離要長于X型布置.在中等間距比和大間距比下，兩種布置形式后渦的分布有較大區(qū)別.對于X型布置，在特定間距比下則會(huì)出現(xiàn)單穩(wěn)態(tài)流動(dòng).

3)當(dāng)L/D較小時(shí)，升、阻力系數(shù)主要受剪切層和鄰近圓柱的影響，而當(dāng)L/D增大，圓柱后的渦開始脫落，脫落渦的影響將逐漸增大，受到上游圓柱脫落渦影響的圓柱，Crmsd和Crms1會(huì)增大.當(dāng)五圓柱呈十字型布置時(shí)，存在一臨界L/D值，當(dāng)L/D大于該值時(shí)，升、阻力系數(shù)的變化規(guī)律會(huì)與L/D小于該值時(shí)的規(guī)律有明顯不同.該值在本計(jì)算中大約在3附近.

4)采用四叉樹自適應(yīng)網(wǎng)格可以減少計(jì)算時(shí)間，提高計(jì)算效率，且對多圓柱繞流問題的計(jì)算結(jié)果令人滿意.

［1］WILLIAMSON C H K.Evolution of a single wake behind a pair of bluff bodies［J］.Journal of Fluid Mechanics，1985，159:1-18.

［2］鄭庭輝，費(fèi)寶玲.串列雙圓柱尾跡流的數(shù)值分析［J］.西南交通大學(xué)學(xué)報(bào)，2008，43(6):747-750.

ZHENG TING HUI，F(xiàn)EI BAO LING.Numerical analysis of flow around two tandem circular cylinders［J］.Journal of Southwest University for Nationalities 2008，43(6):747-750.

［3］GU Z，SUN T.Classifications of flow pattern on three circular cylinders in equilateral-triangular arrangements［J］.Journal of Wind Engineering＆ Industrial Aerodynamics，2001，89(6):553-568.

［4］BAO Y，ZHOU D，HUANG C.Numerical simulation of flow over three circular cylinders in equilateral arrangements at low Reynolds number by a second-order characteristic-based split finite element method［J］.Computers ＆ Fluids，2010，39(5):882 –899.

［5］LAM K，LI J Y，SO R M C.Force coefficients and Strouhal numbers of four cylinders in cross flow［J］.Journal of Fluids ＆ Structures，2003，18:305–324.

［6］BARANYIL.Lift And Drag Evaluation In Translating And Rotating Non-Inertial Systems［J］.Journal of Fluids＆ Structures，2005，20(1):25–34.

［7］BROWN D L，CORTEZ R，MINION M L.Accurate Projection Methods for the Incompressible Navier-Stokes Equations［J］.Journal of Compu tational Physics，2001，168(2):464 – 499.

［8］POPINET S.Gerris:a tree-based adaptive solver for the incompressible Euler equations in complex geometries［J］.Journal of Computational Physics，2003，190(2):572 – 600.

［9］LINNICK M N，F(xiàn)ASEL H F.A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains［J］.Journal of Computational Physics，2005，204(1):157 – 192.

［10］FUJISAWA N，KAWAJI Y，IKEMOTO K.Feedback Control Of Vortex Shedding From A Circular Cylinder By Rotational Oscillations［J］.Journal of Fluids ＆ Structures，2001，15(1):23-37.

［11］HE J-，GLOWINSKI R，METCALFE R，et al.Active Control and Drag Optimization for Flow Past a Circular Cylinder.I.Oscillatory Cylinder Rotation［J］.Journal of Computational Physics，2000，163(1):83-117(35).

［12］HENDERSON R D.Details of the drag curve near the onset of vortex shedding［J］.Physics of Fluids，1995，7(9):2102-2104.

［13］MENEGHINI J，SALTARA F.Numerical Simulation Of Flow Interference Between Two Circular Cylinders In Tandem And Side-By-Side Arrangements［J］.Journal of Fluids ＆ Structures，2001，15(2):327-350.

［14］LAM K，CHEUNG W C.Phenomena of vortex shedding and flow interference of three cylinders in different equilateral arrangements［J］.Journal of Fluid Mechanics，1988，196(196):1-26.

［15］顧志福，孫天風(fēng).三圓柱繞流的實(shí)驗(yàn)實(shí)驗(yàn)研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2000，18(4)，441-447.

GU ZHI FU，SUN TIAN FENG.Experimental Research on Flow A-round Three Circular Cylinders［J］.Acta Aerodynamica Sinica，200，18(4):441-447.

［16］SHARMAN B，LIEN F S，DAVIDSON L，et al.Numerical predictions of low Reynolds number flows over two tandem circular cylinders［J］.International Journal for Numerical Methods in Fluids，2005，47(5):423–447.

Numerical investigation of flow past five circular cylinders at low Reynolds number

ZHU Run-ye，DIAO Ming-jun
(State Key Lab of Hydraulics and Mountain River Engineering，Sichuan University，Chengdu610065，P.R.C.)

In order to have a deeper understanding of flow around multiple circular cylinders，a parallel code which discretizes the 2D incompressible Navier-Stokes equations is employed in this paper and an adaptive octree mesh is also used.The simulation is operated at Re=100.The influence of cylinder-to-cylinder spacing ratio and the arrangements of five cylinders are investigated.The results of the simulation suggested that there were three different flow patterns connected to various spacing ratio;there were also differences between two arrangements of cylinders when the spacing ration was the same.The spacing ration has different effects on the coefficients of drag and lift at X arrangement and cruciform arrangement.

flow pastfive circular cylinders;octree mesh;numerical simulation;fluid force

TV131.4

A

2095-4271(2015)06-0758-09

10.11920/xnmdzk.2015.06.020

2015-09-23

朱潤野(1991－)，男，浙江溫州人，碩士研究生，研究方向:水工水力學(xué)，Email:273216779@qq.com

刁明軍(1968－)，男，四川簡陽人，教授，博士生導(dǎo)師，研究方向:工程水力學(xué)的科研與教學(xué)工作.Email:diaomingjun@scu.edu.cn

四川省學(xué)術(shù)和技術(shù)帶頭人培養(yǎng)基金(2012DTY020)

(責(zé)任編輯:付強(qiáng)，張陽，李建忠，羅敏;英文編輯:周序林)