陣列四柱繞流的數(shù)值模擬

劉為民, 谷家揚, 盧燕祥

(江蘇科技大學(xué) 船舶與海洋工程學(xué)院, 江蘇 鎮(zhèn)江 212003)

多鈍體繞流問題廣泛存在于各個工程領(lǐng)域中,例如,橋梁工程中的橋墩、土木工程中的建筑群、航天工程中的支撐機翼的雙柱、海洋工程中的柔性立管以及海洋鉆井平臺的立柱等．方柱群繞流是一種典型的鈍體繞流,一直是流體領(lǐng)域熱點研究對象之一,流體會對方柱產(chǎn)生力的作用,繞過流體后的流體也會相互作用．各立柱周圍的流場結(jié)構(gòu)非常復(fù)雜,各立柱動力學(xué)特征和單獨置于流體中的單個立柱動力學(xué)特征完全不同．

風(fēng)洞和水池的實驗研究涉及大量的人力、物力,而數(shù)值模擬具有投資小、周期短、效率高等特點,因此,成為一般研究人員的主要研究手段．多柱體繞流無論是串列、并列布置方式,相鄰圓柱間的渦脫對最終流體力和旋渦脫落會產(chǎn)生很大的影響,與單圓柱相比會有很大區(qū)別．當間距不同時會在柱體后方形成不同尾流．當前對多柱體靜止繞流的數(shù)值及試驗研究,多是基于等直徑樁柱群進行的[1-2]．

文獻[3]采用有限體積法對流向分別為0°及4°的方柱在空氣中的繞流進行了數(shù)值模擬并進行了比較,研究發(fā)現(xiàn)特征長度對柱體尾部泄渦結(jié)構(gòu)有很大影響．文獻[4]采用大渦模擬與二階全展開ETG有限元離散格式相結(jié)合的方法對并列方柱繞流進行了數(shù)值模擬,模型中,柱間距比為1.5,對兩柱體的水動力系數(shù)及對稱布置點上的流向速度進行了時域和頻域分析,研究結(jié)果表明在來流對稱條件下,兩方柱繞流時域不對稱而頻域?qū)ΨQ．文獻[5]運用有限體積法對間距比在1.5～5.0之間每間隔0.5共計8個典型間距比的串列雙圓柱的二維繞流進行了數(shù)值計算,重點考察間距比對圓柱相互作用和流場的影響,并發(fā)現(xiàn)間距大小很大程度上決定著流動特性．文獻[6]采用有限元法對相似串列方柱構(gòu)建物在雷諾數(shù)等于100時的繞流進行了數(shù)值研究,重點分析了間距比對流場及水動力系數(shù)的影響,研究發(fā)現(xiàn),方形圓柱串列時,流體阻力有效降低,并獲得臨界間距比為4.50～4.75,此時流體各項力學(xué)性能發(fā)生跳躍．文獻[7]采用分區(qū)覆蓋網(wǎng)格算法,對附屬小圓柱置于主圓柱適當位置時的繞流特性進行了數(shù)值模擬,探討了附屬小圓柱擺放位置、大小、個數(shù)及雷諾數(shù)對主圓柱繞流的影響及流場演變規(guī)律．文獻[8]采用表面渦法對高雷諾數(shù)下不同排列方式雙圓柱繞流的流動狀態(tài)進行了計算,考慮了并列、串列及級列下流動結(jié)構(gòu)、尾渦結(jié)構(gòu)變化及圓柱受力特性．文獻[9]對空氣中兩串列方柱的干擾問題進行了數(shù)值模擬,采用改進的標志網(wǎng)格法,重點考察了間距比的影響．文獻[10]在北京大學(xué)理學(xué)與工程科學(xué)中型大氣邊界層風(fēng)洞中對成等邊三角形排列的3個同直徑圓柱在不同風(fēng)向角下的表面壓力進行了實驗研究,研究發(fā)現(xiàn)在小間距比時,風(fēng)向角變化會產(chǎn)生有剪切層、鄰近和尾流3種不同機理力引起的干擾,且處于下游位置圓柱的內(nèi)側(cè)或外側(cè)會出現(xiàn)大面積的負壓分布從而產(chǎn)生很大的橫向力．文獻[11]基于有限體積法對低雷諾數(shù)下正三角形排列等直徑三圓柱的繞流和渦激運動進行了數(shù)值模擬,間距比取為1.5～6.0,重點對各圓柱的氣動力、響應(yīng)及相關(guān)頻率特性、尾流結(jié)構(gòu)模式隨間距比的變化進行了研究,研究表明：三圓柱振蕩氣動力遠大于固定繞流,間距比為3.5～4.0時能產(chǎn)生最強的啟動干擾．文獻[12]通過實驗對2個串列方柱在不同間距比下的流場特性進行了數(shù)值研究,其采用的雷諾數(shù)為5 300和16 000．文獻[13]數(shù)值分析了雷諾數(shù)100時2個串列方柱的繞流特性,計算中考慮了5 個間距比,但未給出臨界間距比．

1 數(shù)值計算模型

1.1 流體力學(xué)控制方程

不可壓縮粘性流體的控制方程為質(zhì)量和動量守恒方程:

(1)

(2)

湍流模型采用SSTk-ω模型,湍動能k與比耗散率ω的輸運方程為：

(3)

(4)

1.2 計算模型與邊界條件設(shè)置

文中以四邊形張力腿平臺為例．此平臺由4個方形浮筒和4個方形浮箱組成,浮筒直徑17.40 m,浮筒中心間距51.40 m,浮箱高度8.70 m,浮箱寬度11.60 m,浮箱長度34.00 m,排水量45 536 t,設(shè)計吃水25.35 m．針對張力腿平臺4個方形立柱,分別建立了來流方向為0°，22.5°，45°模型．網(wǎng)格模型采用縮尺比為1∶40模型．模型均為50D×20D(D為立柱邊長)大小的矩形區(qū)域,其中原點位于四立柱兩斜對角線交點處,原點距離上游入口處為15D,距離下游出口處為35D,距離上下兩邊界為15D．流體域采用了混合型網(wǎng)格,對立柱表面采用邊界層,近場區(qū)域采用結(jié)構(gòu)化網(wǎng)格,遠場區(qū)域采用非結(jié)構(gòu)化網(wǎng)格;同時對立柱表面和尾流區(qū)域等參數(shù)梯度變化較大的區(qū)域進行了局部網(wǎng)格加密．來流方向為0°，22.5°及45°的網(wǎng)格單元總數(shù)分別為90 972,90 464,89 932．網(wǎng)格模型具體如圖1,為了對柱間干擾進行分析,這里定義上游上方立柱為column1,上游下方立柱為column2,下游上方立柱為column3,下游下方立柱為column4．

a) 全局計算網(wǎng)格(0°)

b) 近場網(wǎng)格(0°)

c) 全局計算網(wǎng)格(22.5°)

d) 近場網(wǎng)格(22.5°)

e) 全局計算網(wǎng)格(45°)

f) 近場網(wǎng)格(45°)

圖1網(wǎng)格模型

Fig.1Meshmodels

流動方向自左向右,不同來流方向模型的邊界條件設(shè)置為:左側(cè)為速度入口邊界,右側(cè)為完全發(fā)展出流邊界,上下兩側(cè)為自由滑移邊界,立柱表面為無滑移邊界．由于采用縮尺比模型,因此在數(shù)值模擬中需要對來流速度進行轉(zhuǎn)換,實際來流速度與模型來流速度轉(zhuǎn)換關(guān)系為:

(5)

式中:Um為模型來流速度;Us為實際來流速度;λ為縮尺比,λ=40．

2010年馬來西亞科研機構(gòu)對中國南海Borneo附近垂向流速隨水深的變化進行了實測[14],根據(jù)上述文獻結(jié)合實際,選取6個典型流速,見表1．

表1 計算的典型流速

2 計算結(jié)果與分析

2.1 單方柱繞流特性

作為與陣列方柱繞流問題對比的依據(jù)和對照,對上述6個流速下單方柱繞流問題首先進行計算,并對升力系數(shù)曲線做傅里葉變換得到頻譜圖,其峰值即為斯特哈爾數(shù)(St)．研究中發(fā)現(xiàn),同一來流方向下,固定方柱的渦泄頻率隨來流速度的增加而不斷增大,這和實驗中觀測到的柱體后方渦脫頻率隨流速成正比相吻合．將文中方柱在0°浪向6個不同流速的斯特哈爾數(shù)進行平均,與相關(guān)實驗和數(shù)值對比(表2),采用的SSTk-ω計算的結(jié)果跟實驗結(jié)果已經(jīng)十分接近,可見該方法在預(yù)報方柱繞流相關(guān)特性方面具有很高的可信度．圖2給出了不同來流方向時斯特哈爾數(shù)隨模型速度的變化,從圖中可以看出,同一來流St隨速度有輕微振蕩,但是變化幅度不大,因為斯特哈爾數(shù)是一個無因次物理量,反映了渦脫相對于來流速度快慢．對不同來流方向:0°時,St=0.138；22.5°時，St=0.161；而45°時St最大,達到了0.201.斯特哈爾數(shù)隨這3個流向角度增大而不斷增加．圖3給出了22.5°流向，Um=0.158 m·s-1時升力系數(shù)及拖曳力系數(shù)時歷曲線．

圖2 單立柱不同來流及流速下的斯特哈爾數(shù)

圖3 升力系數(shù)及拖曳力系數(shù)時歷曲線

模型斯特哈爾數(shù)(St)模型斯特哈爾數(shù)(St) 文獻[15]試驗結(jié)果0.132±0.004LES模型[18]0.134 文獻[16]試驗結(jié)果0.139顯示代數(shù)應(yīng)力模型[19]0.150 k-ε模型[17]0.124LES模型[18]0.130 雷諾應(yīng)力方程(壁面函數(shù))[17]0.136SSTk-ω模型[20]0.129 雷諾應(yīng)力方程(2層)[17]0.159文中SSTk-ω模型0.138

圖4給出了Um=0.158 m·s-1,對應(yīng)于Us=1.0 m·s-1,柱體后方渦脫結(jié)構(gòu)圖．從圖中可以看出,旋渦在柱體后方交替脫落．由于在柱體前方壓力逐漸減小,而柱體后方壓力不斷增大,在流動區(qū)域產(chǎn)生分離,在柱體后方產(chǎn)生了波動度很高的速度從而形成渦,被稱為卡門渦街．由于45°流向時,柱體繞流的特征長度為0°流向的1.414倍,所以流向45°時,渦的影響要比0°大的多,這也證明了沿流向特征長度對柱體繞流的影響．不同流向下渦泄結(jié)構(gòu)也有不同,在0°流向時,其尾渦結(jié)構(gòu)為經(jīng)典的卡門渦街,尾渦主要由立柱的上方泄放出,雖然下方也釋放一個單獨的渦,但是逐漸減弱,并被上方釋放出的尾渦吞噬.22.5°時尾渦從立柱下方泄放,開始時為2P模式,但是隨著后方的渦強逐漸減弱,逐漸演化為單渦結(jié)構(gòu).45°時尾渦泄放點為立柱上方,開始時尾渦結(jié)構(gòu)為P+S模式,但是走勢很快與0°，22.5°下游尾渦結(jié)構(gòu)相同,顯示為單渦．

a) 流向0°

b) 流向22.5°

c) 流向45°

圖4Um=0.158m·s-1不同流向的尾渦結(jié)構(gòu)

Fig.4VortexsheddingmodeatcurrentvelocityUm=0.158m·s-1underdifferentdirections

2.2 四柱繞流特性

在海洋工程實際問題中,研究的熱點和重點是平臺所受到的流體力以及相應(yīng)的動力響應(yīng)．然而,在深入研究張力腿平臺渦激運動這一復(fù)雜的流固耦合現(xiàn)象時,對流場自身也需加以特別地關(guān)注,平臺立柱周圍的尾流軌跡,流動分離及渦泄結(jié)構(gòu)是研究流場的重要手段．

對單方柱進行繞流模擬及驗證基礎(chǔ)之上,對張力腿平臺四柱繞流進行深一步的研究和分析．張力腿平臺4個立柱排列與單個立柱的流場特性有著本質(zhì)區(qū)別,當前研究重點主要為串列或并列柱體的繞流,而張力腿四立柱從形式上來說既是并列又是串列,故其流場特性,例如邊界層分離、柱體表面壓力系數(shù)的分布、柱間相互干擾、尾渦脫落模式等都有待進行詳細探討．

a)Um=0.118 6 m·s-1

b)Um=0.197 6 m·s-1

c)Um=0.276 7 m·s-1

圖50°來流時升力系數(shù)時歷曲線

Fig.5Timehistoryofliftcoefficientatcurrentdirection0°

由于文章篇幅有限,僅列出了0°來流條件下,3個典型速度下4個立柱的升力系數(shù)和拖曳力系數(shù)曲線(圖5，6)．無論升力系數(shù)還是拖曳力系數(shù)隨著來流速度的增大,其周期都不斷減?。魉佥^小時,0°來流時各立柱的升力系數(shù)時歷曲線不是規(guī)則的曲線,見圖5a)；隨著流速的增加,各柱的升力系數(shù)逐漸規(guī)則,周期性逐漸增強,如圖5b)；當流速繼續(xù)增加到0.276 7 m·s-1時,由于由于column1,column2位于上游,來流對其影響較小,下游立柱受到上游立柱后方來流的影響,流動形態(tài)波動性增大．縱觀上述3個典型流速下,column1，column2的升力系數(shù)振幅較小,而column3，column4振幅明顯加大．再看拖曳力系數(shù)隨時歷程曲線,在某些流速下,如6b)由于上游立柱的影響,下游立柱的拖曳力系數(shù)曲線具有鋸齒狀,很好地證明了上下游立柱流場之間的耦合作用,此流速下也觀測到了“同相”和“反相”現(xiàn)象,此時column1與column2, column3與column4拖曳力系數(shù)基本吻合在一起,相位也基本相同,同時出現(xiàn)波峰和波谷,而前后兩柱column1與column3,column2與column4卻出現(xiàn)了反相現(xiàn)象,例如column1出現(xiàn)峰值時,column3卻為谷值．

圖7為45°來流3個典型流速時張力腿平臺四立柱所受到的總體升力系數(shù)及其對應(yīng)的頻譜變化,盡管四立柱繞流之間存在耦合作用,但是從頻譜轉(zhuǎn)化上來看,其每個流速都只對應(yīng)著一個頻率峰值．

a)Um=0.118 6 m·s-1

b)Um=0.197 6 m·s-1

c)Um=0.276 7 m·s-1

圖60°來流時拖曳力系數(shù)時歷曲線

Fig.6Timehistoryofdragcoefficientatcurrentdirection0°

a)Um=0.118 6 m·s-1

b)Um=0.197 6 m·s-1

c)Um=0.276 7 m·s-1

圖745°來流時四立柱升力系數(shù)時歷曲線

Fig.7Timehistoryofliftcoefficientoffourcolumnsatcurrentdirection45°

渦泄頻率跟單個方柱的繞流渦泄頻率變化趨勢相同,隨著流速增大而增加,但是斯特哈爾數(shù)45°流向時最大,6個速度平均值達到了0.299 m·s-1,而0°流向次之,為0.218 m·s-1,22.5°流向最小僅為0.197 m·s-1．

a)Um=0.079 m·s-1

b)Um=0.119 m·s-1

c)Um=0.158 m·s-1

d)Um=0.198 m·s-1

e)Um=0.237 m·s-1

f)Um=0.227 m·s-1

圖845°流向尾渦結(jié)構(gòu)

Fig.8Vortexsheddingmodeatcurrentdirection45°

文獻[20]于2000年在對圓柱渦激運動流場研究時,根據(jù)尾流旋渦結(jié)構(gòu)的生成和演化方式,大致將圓柱后方的尾流流場分為三組不同的形態(tài),分別對應(yīng)著2S模式,2P模式以及P+S模式．文中僅列出了45°流向6個典型速度下的尾渦結(jié)構(gòu)圖,從圖8中可以發(fā)現(xiàn),雖然來流速度存在很大差別,旋渦脫離立柱的方式基本一致,開始脫離時,都是從立柱上、下方各泄放出一個旋渦,但隨著遠離流場,逐漸形成上下兩行獨立的旋渦,尾流結(jié)構(gòu)形式為典型的2S模式．

3 結(jié)論

文中在不考慮浮箱對流場影響的情況下,引入雷諾平均法求解NS方程，結(jié)合SSTk-ω湍流模型對張力腿平臺四立柱繞流特性進行了數(shù)值模擬,為簡化研究,僅考慮了二維立柱不同來流和不同流速下的繞流特性．在對單個方柱繞流進行實驗數(shù)據(jù)對比的基礎(chǔ)上,采用相同方法對四立柱繞流進行了計算,得出了如下結(jié)論:

1) 無論單立柱還是張力腿平臺的四立柱,來流速度對升力系數(shù)和阻力系數(shù)的時歷曲線影響很大,在來流方向上,由于立柱流場之間的耦合效應(yīng),上游立柱的升力系數(shù)和拖曳力系數(shù)的時歷曲線較后方立柱較規(guī)則；

2) 陣列四立柱升力系數(shù)頻譜變換后得到的渦泄頻率在0°及45°流向時均大于單方柱的渦泄頻率,而在22.5°流向時,四立柱渦泄頻率基本等于單立柱的渦泄頻率；

3) 張力腿平臺四柱繞流時,立柱表面渦強最大值集中在上游立柱迎流面的拐點處,此處為流場分離點,但隨流向迅速減?。煌瑏砹鞣较蚣皝砹魉俣认?立柱后方的尾流結(jié)構(gòu)形式均為2S模式．

[1] 朱仁慶,侯玲.LNG船液艙晃蕩數(shù)值模擬[J].江蘇科技大學(xué)學(xué)報:自然科學(xué)版,2010,24(1):1-6.

Zhu Renqing, Hou Ling. Numerical simulation of liguid sloshing in the tanks of LNG carrier[J].JournalofJiangsuUniversityofScienceandTechnology:NaturalScienceEdition,2010,24(1):1-6.(in Chinese)

[2] 朱仁慶,李晏丞,倪永燕,等.氣泡在水中上升運動的數(shù)值模擬[J].江蘇科技大學(xué)學(xué)報:自然科學(xué)版,2010,24(5):417-422.

Zhu Renqing, Li Yancheng, Ni Yongyan, et al.Numerical simulation of bubble rising in the water[J].JournalofJiangsuUniversityofScienceandTechnology:NaturalScienceEdition,2010,24(5):417-422. (in Chinese)

[3] 黃娟,張志國, Mudaliar A V,等.2D非穩(wěn)態(tài)方柱繞流的數(shù)值模擬與分析[J].武漢理工大學(xué)學(xué)報:交通科學(xué)與工程版,2004,33(1):192-195.

Huang Juan,Zhang Zhiguo, Mudaliar A V,et al. 2-D unsteady computational analysis of square cylinders in cross flow[J].JournalofWuhanUniversityofTechnology:TransportationScience&Engineering, 2004,33(1):192-195. (in Chinese)

[4] 魏英杰,朱蒙生,何鐘怡.并列雙方柱繞流的大渦模擬及頻譜分析[J].應(yīng)用數(shù)學(xué)和力學(xué),2009,25(8):824-830.

Wei Yingjie,Zhu Mengsheng,He Zhongyi.Large eddy simulation and spectrum analysis of flow around two square cylinders arranged side by side[J].AppliedMathematicsandMechanics,2009,25(8):824-830. (in Chinese)

[5] 吳七二,周岱,付功義.串列方形鈍體構(gòu)筑物繞流的數(shù)值模擬[J].上海交通大學(xué)學(xué)報,2012,46(1):130-135.

Wu Qier,Zhou Dai,Fu Gongyi.Numerical analysis of flow past two bluff bodies with square-section in tandem arrangement[J].JournalofShanghaiJiaotongUniversity,2012,46(1):130-135. (in Chinese)

[6] 劉松,符松,串列雙圓柱擾流問題的數(shù)值模擬[J].計算力學(xué)學(xué)報,2000,17(3):260-266.

Liu Song,Fu Song. Numerical simulation of flow past two cylinders in tandem arrangement[J].ChineseJournalofComputationalMechanics,2000, 17(3):260-266. (in Chinese)

[7] 李偆萱,彭少波,吳子牛.附屬小圓柱對主圓柱繞流影響的數(shù)值模擬[J].北京航空航天大學(xué)學(xué)報,2003,29(11):951-958.

Li Chunxuan,Peng Shaobo,Wu Ziniu.Numerical study of flow around a main cylinder by controlled satellite cylinders[J].JournalofBeijingUniversityofAreonauticsandAstronautics,2003,29(11):951-958. (in Chinese)

[8] 廖俊,景思睿.高雷諾數(shù)下圓柱繞流的數(shù)值模擬[J].水動力學(xué)研究與進展,2001,16(1):101-110.

Liao Jun,Jing Sirui.Numerical simulation of flow around two cylinders at high Reynolds number[J].JournalofHydrodynamics,2001,16(1):101-110. (in Chinese)

[9] 陳素琴,黃自萍,沈劍華,等.兩串列方柱繞流的干擾數(shù)值研究[J].同濟大學(xué)學(xué)報,2001,29(3):320-325.

Chen Suqin,Huang Ziping,Shen Jianhua,et al.Numerical computation of flow around two square cylinders in tandem arrangement[J].JournalofTongjiUniversity,2001,29(3):320-325. (in Chinese)

[10] 顧志福,孫天風(fēng).三圓柱繞流的實驗研究[J].空氣動力學(xué)學(xué)報,2000,18(4):442-447.

Gu Zhifu,Sun Tianfeng.Experimental research on flow around three circular cylinders[J].ActaAerodynamicaSinica,2000,18(4):442-447. (in Chinese)

[11] 徐楓,歐進萍.正三角形排列三圓柱繞流與渦激振動數(shù)值模擬[J].空氣動力學(xué)學(xué)報,2010,28(5):582-590.

Xu Feng,Ou Jinping. Numerical simulation of vortex-induced vibration of three cylinders subjected to a cross flow in equilateral arrangement[J].ActaAerodynamicaSinica,2010,28(5):582-590. (in Chinese)

[12] Kim M K,Kim D K,Yoon S H,et al.Measur ements of the flow fields around tw o square cylinders in a tandem arrangement[J].JournalofMechanicalScienceandTechnology,2008,22(2):397-407.

[13] Lankadasu A,Veng adesan S.Interference effect of two equal sized square cylinders in tandem arrang ement:with planar shear flow[J].InternationalJournalofNumericalMethodsinFluids,2008,57(8):1005-1021.

[14] Rizwan Sheikh,Alison Brown．Extreme vertical deepwater current profiles in the south Chinese,offshore Borneo[C]∥ProceedingsoftheASME2010 29thInternationalConferenceonOcean,OffshoreandArcticEngineering.[s.l.]:OMAE,2010-49151.

[15] Lyn D A,Einav S,Rodi W.Laser-doppler velocimetry study of an semble-averaged characteristics of the turbulent near wake of a square cylinder[J].JournalofFluidMechanics,1995,304:285-292.

[16] Durao D F G,Heitor M V,Pereira J C F.Measurements of turbulent and periodic flows around a square cross-secti on cylinder[J].ExperimentsinFluids,1988,6(5):298-309.

[17] Franke R, Rodi W. Calculati on of vortex shedding past a square cylinder with various turbulence models[C]∥Proceedingsofthe8thInternationalSymposiumonTurbulentShearFlows.New York:Springer-Verlag Press,1991:189-204.

[18] Bouris D,Bergeles G.2D LES of vortex shedding from a square cylinder[J].JournalofWindEngineeringandIndustrialAerodynamics,1999,80:31-42.

[19] Lbcke H.Comparison of LES and RANS in bluff-body flows[J].JournalofWindEngineeringandIndustrialAerodynamics, 2001,89:1471-1483.

[20] Rodi W, Ferziger J,Breuer M.Status of large eddy simulation: results of a workshop[J].JournalofFluidsEngineering, 1997,119(2):248-261.

[21] 張偉,葛耀君.方柱繞流粒子圖像測速試驗與數(shù)值模擬[J].同濟大學(xué)學(xué)報:自然科學(xué)版,2009,37(7):857-892.

Zhang Wei,Ge Yaojun.Partical image velocity study and numerical simulation of turbulent near wake of square cylinder[J].JournalofTongjiUniversity:NaturalScience,2009,37(7):857-892. (in Chinese)

[22] Govardban R,Williamson C H K.Modes of vortex formation and frequency response of a freely vibrating cylinder[J].JournalofFluidMechanics,2000,420:85-130.