不同工況下有限長(zhǎng)圓柱繞流實(shí)驗(yàn)測(cè)量及分析

梁圣召，周軍偉，梅 蕾

(哈爾濱工業(yè)大學(xué)（威海）船舶與海洋工程學(xué)院，山東 威海 264209)

0 引 言

不同長(zhǎng)徑比、倒立入水的有限長(zhǎng)圓柱流體流動(dòng)特性對(duì)很多工程應(yīng)用，包括浮式海上結(jié)構(gòu)和橋梁渦振都有重要的研究?jī)r(jià)值。因?yàn)橛邢揲L(zhǎng)圓柱的尾流和自由端的流動(dòng)相當(dāng)復(fù)雜，僅用實(shí)驗(yàn)方法不能滿(mǎn)足相似參數(shù)和相似定律的要求，而且各項(xiàng)數(shù)據(jù)的測(cè)量以及尾流和自由端流動(dòng)狀態(tài)的觀察比較困難。實(shí)驗(yàn)的這些不足可以由數(shù)值模擬來(lái)補(bǔ)充，但是截至目前，研究不同長(zhǎng)徑比倒立入水的有限長(zhǎng)圓柱繞流流動(dòng)特性的文獻(xiàn)還不夠充分。

早期的研究以實(shí)驗(yàn)為主，Okamoto 和Yagita[1]通過(guò)測(cè)量長(zhǎng)徑比AR=1-12.5 圓柱表面的壓力分布得到了圓柱平均阻力系數(shù)，并分析了自由端的分離現(xiàn)象。H.Akilli和D.Rockwell[2]利用粒子圖像測(cè)速法（PIV）對(duì)倒立入水的長(zhǎng)徑比為AR=1 和AR=2 的有限長(zhǎng)圓柱進(jìn)行了研究，研究結(jié)果形象解釋了大尺度漩渦的形成。J.R.Chaplin等[3]在雷諾數(shù)Re與傅汝德數(shù)Fr的比值恒等于2.79×105的條件下對(duì)有自由液面的倒立入水圓柱在拖曳水池中進(jìn)行了阻力實(shí)驗(yàn)研究。張榮譽(yù)[4]實(shí)驗(yàn)研究了具有自由液面的圓柱繞流現(xiàn)象，測(cè)量了不同F(xiàn)r數(shù)工況下尾流場(chǎng)的變化規(guī)律。

隨著CFD 技術(shù)的發(fā)展，Lee 等[5]采用CFD 方法（大渦模擬）對(duì)2 種長(zhǎng)徑比AR=2.5 和10 的壁面固定的有限長(zhǎng)圓柱進(jìn)行了數(shù)值模擬研究。Afgan 等[6]也采用大渦模擬在雷諾數(shù)Re=2×104和長(zhǎng)徑比AR=6 和10 的情況下對(duì)壁面固定的有限長(zhǎng)圓柱進(jìn)行了研究。Guilherme F.Rosetti 等[7]在雷諾數(shù)Re=4.3×104和傅汝德數(shù)Fr=0.31 的條件下，對(duì)長(zhǎng)徑比ui=2倒立入水的有限長(zhǎng)圓柱從實(shí)驗(yàn)和數(shù)值模擬兩方面進(jìn)行了研究。發(fā)現(xiàn)自由表面的存在會(huì)影響流動(dòng)，但是占主導(dǎo)作用的還是自由端的影響。M.A.Benitz 等[8]在雷諾數(shù)Re=2 900 和傅汝德數(shù)Fr=0.65 的條件下，對(duì)長(zhǎng)徑比從1～19 倒立入水的有限長(zhǎng)圓柱進(jìn)行了實(shí)驗(yàn)研究，并用大渦模擬進(jìn)行了對(duì)比驗(yàn)證。

目前學(xué)者們已經(jīng)探討了部分雷諾數(shù)，傅汝德數(shù)以及長(zhǎng)徑比下圓柱繞流的流動(dòng)特性，同時(shí)開(kāi)展了以大渦模擬為主的數(shù)值方法驗(yàn)證。但是Re和Fr的覆蓋范圍還不夠全面，探討還不夠充分。本文首先采用實(shí)驗(yàn)測(cè)量和數(shù)值模擬相結(jié)合的方法對(duì)不同長(zhǎng)徑比圓柱的阻力系數(shù)進(jìn)行對(duì)比分析，而后對(duì)有限長(zhǎng)圓柱的尾流場(chǎng)特征進(jìn)行討論，為海工結(jié)構(gòu)設(shè)計(jì)提供一定的數(shù)據(jù)參考。

1 數(shù)值模擬方法

1.1 控制方程

對(duì)于三維非定常不可壓縮粘性流動(dòng)，其連續(xù)性方程和動(dòng)量方程為：

式中：ui，uj為速度分量；p為壓力；v為總的流體運(yùn)動(dòng)黏性系數(shù)，表示為層流黏性系數(shù)與湍流黏性系數(shù)之和，其中湍流黏性系數(shù)由湍流模型獲得；fi為外力分量。

本文采用Wilcox 的k-ω湍流模型，具體可參照文獻(xiàn)[9]。基于CFX 求解器實(shí)現(xiàn)方程的求解。

1.2 計(jì)算域模型

幾何模型如圖1 所示，長(zhǎng)×寬×高為4 m×2 m×1.2 m，相當(dāng)于80D×40D×24D，其中特征長(zhǎng)度D=0.05 m。坐標(biāo)原點(diǎn)在圓柱底面的中心處。計(jì)算域入口距離圓柱中心為20D，出口距離圓柱中心60D，兩側(cè)面距離圓柱中心20D。

圖1 有限長(zhǎng)圓柱繞流幾何模型Fig.1 Model of computational domain for finite cylinder

流體域進(jìn)口設(shè)定為速度進(jìn)口（inlet）；出口設(shè)定為壓力出口邊界（outlet）；流體域上表面以及2 個(gè)側(cè)面均設(shè)定為自由滑移壁面（free slip wall）,下表面設(shè)置為無(wú)滑移壁面（no slip wall）；圓柱表面包括自由端面設(shè)定為無(wú)滑移壁面（no slip wall）。

1.3 數(shù)值方法驗(yàn)證

1）分別進(jìn)行網(wǎng)格數(shù)量和湍流模型對(duì)計(jì)算結(jié)果和計(jì)算時(shí)間的對(duì)比分析。計(jì)算條件為：流場(chǎng)無(wú)窮遠(yuǎn)處的均勻來(lái)流速度為U=0.6 m/s，流體為水，密度ρ=1 000 kg/m3，運(yùn)動(dòng)粘性系數(shù)ν=10-6m2/s，計(jì)算得雷諾數(shù)Re=3×104，時(shí)間步長(zhǎng)0.001 s。網(wǎng)格在圓柱周向均布100 個(gè)節(jié)點(diǎn)，法向擠出20 層，線(xiàn)性增長(zhǎng)率1.3。近壁面處保證y+＜2。3 種網(wǎng)格密度采用同樣的加密區(qū)域7D×3D，采用四面體、棱柱型與六面體網(wǎng)格相結(jié)合的方式，以改善網(wǎng)格質(zhì)量，計(jì)算域網(wǎng)格劃分如圖2 所示。3 套網(wǎng)格分別命名為Coarse，Medium 和Fine，網(wǎng)格數(shù)分別為49 萬(wàn)、119 萬(wàn)和200 萬(wàn)。表1 為k-ω湍流模型不同網(wǎng)格密度條件下，以及在Medium 網(wǎng)格下分別采用DES，k-ω和k-ε三種湍流模型所得到的升阻力系數(shù)以及計(jì)算所耗時(shí)間?？梢钥闯霎?dāng)采用k-ω模型時(shí)，使用Medium 和Fine 網(wǎng)格，Cd計(jì)算結(jié)果相差2.7%，Cl計(jì)算結(jié)果相差9.93%，但是計(jì)算時(shí)間相差1 倍，所以Medium 網(wǎng)格數(shù)計(jì)算精度足夠。而在Medium 網(wǎng)格數(shù)下，k-ε的計(jì)算結(jié)果與DES 和k-ω的結(jié)果相差較大，DES 和k-ω的計(jì)算精度最大誤差在15%以?xún)?nèi)，但是計(jì)算時(shí)間相差近1 倍。綜合考慮計(jì)算精度與計(jì)算時(shí)間成本，本文選取Medium 網(wǎng)格數(shù)119 萬(wàn)，湍流模型選擇k-ω。

圖2 不同密度的網(wǎng)格劃分Fig.2 Grid with different densities

表1 不同網(wǎng)格密度和湍流模型的比較Tab.1 Numerical results of different mesh refinement and different turbulence model

2）鑒于文獻(xiàn)在研究相關(guān)問(wèn)題時(shí)多采用大渦模擬（LES）湍流模型，為了進(jìn)一步驗(yàn)證本文計(jì)算模型的可用性，在同等初始條件下，對(duì)不同工況及長(zhǎng)徑比的有限長(zhǎng)圓柱進(jìn)行計(jì)算，并與前人所做相近工況的實(shí)驗(yàn)和數(shù)值模擬結(jié)果進(jìn)行對(duì)比。首先，本文計(jì)算和實(shí)驗(yàn)結(jié)果與前人所做相近工況的實(shí)驗(yàn)和數(shù)值模擬結(jié)果進(jìn)行對(duì)比，如表2 所示?？梢钥闯觯涸陂L(zhǎng)徑比AR=6，雷諾數(shù)Re=3×104工況下，本文計(jì)算結(jié)果比H.Sakamoto[10]的實(shí)驗(yàn)值大6.59%，實(shí)驗(yàn)結(jié)果比H.Sakamoto[10]的實(shí)驗(yàn)值小12.3%；在長(zhǎng)徑比AR=10，雷諾數(shù)Re=3×104工況下本文的計(jì)算結(jié)果比LEE[5]的計(jì)算值大5.18%，實(shí)驗(yàn)值比其小2.56%。與相似工況下的文獻(xiàn)結(jié)果的誤差均在可以接受范圍內(nèi)。同時(shí)，在長(zhǎng)徑比AR=10，雷諾數(shù)Re=3×104工況下，本文的實(shí)驗(yàn)結(jié)果與計(jì)算結(jié)果也進(jìn)行了對(duì)比，如圖3 所示。可知實(shí)驗(yàn)值與計(jì)算值的升阻力系數(shù)時(shí)程圖比較接近，而橫流向的振動(dòng)頻率約是順流向振動(dòng)頻率的2 倍，這與王曉聰[11]的研究結(jié)論相符，也側(cè)面佐證了本文研究的準(zhǔn)確性。最終選取k-ω湍流模型進(jìn)行數(shù)值分析。

圖3 實(shí)驗(yàn)與計(jì)算升阻力系數(shù)時(shí)程圖Fig.3 Distribution of the time-dependent lift and drag coefficient for EXP and CFD

表2 不同作者相似工況計(jì)算結(jié)果Tab.2 Calculation results of similar conditions for different authors

2 實(shí)驗(yàn)測(cè)試

實(shí)驗(yàn)在哈爾濱工業(yè)大學(xué)（威海）海洋工程學(xué)院的循環(huán)水槽實(shí)驗(yàn)室進(jìn)行，三維模型如圖4 所示。實(shí)驗(yàn)所用圓柱材質(zhì)為均質(zhì)6061 鋁合金，直徑50 mm，圓柱表面經(jīng)過(guò)高速車(chē)削處理，以達(dá)到較小的粗糙度。實(shí)驗(yàn)中，測(cè)力單元處于水面之上，而底板略低于水面，以盡量減小自由液面對(duì)測(cè)量結(jié)果的影響。

圖4 三維模型Fig.4 3D model

實(shí)驗(yàn)來(lái)流速度U選擇為0.2，0.4 和0.6 m/s，對(duì)應(yīng)的雷諾數(shù)分別是1×104，2×104和3×104。測(cè)試3 個(gè)長(zhǎng)徑比的圓柱，其高度分別為H=150，300 和500 mm，對(duì)應(yīng)長(zhǎng)徑比為AR=3，6 和10。實(shí)驗(yàn)中對(duì)測(cè)力計(jì)進(jìn)行細(xì)致的標(biāo)定，以保證測(cè)量結(jié)果的準(zhǔn)確度。3 個(gè)不同的來(lái)流速度以及3 個(gè)長(zhǎng)徑比組合共9 組實(shí)驗(yàn)，每組實(shí)驗(yàn)進(jìn)行4 次測(cè)量，取平均值。

3 結(jié)果與分析

3.1 阻力系數(shù)分析

將不同長(zhǎng)徑比下，數(shù)值模擬與實(shí)驗(yàn)研究所得阻力系數(shù)的均方根值進(jìn)行對(duì)比，其與雷諾數(shù)的關(guān)系如圖5所示。研究發(fā)現(xiàn)，當(dāng)長(zhǎng)徑比AR=3，6 和10 時(shí)，阻力系數(shù)均隨流速的增加有1 個(gè)逐漸減小的趨勢(shì)，這與R.T.Gon?alves[13]的實(shí)驗(yàn)結(jié)果相一致，原因是流速越大，自由端對(duì)流動(dòng)的影響越大。但阻力系數(shù)總體變化不大，表明在長(zhǎng)徑比大于3 的工況下，有限長(zhǎng)圓柱繞流的阻力系數(shù)趨于平穩(wěn)，其值在0.61 與0.86 之間。但同時(shí)也發(fā)現(xiàn)，本文的實(shí)驗(yàn)結(jié)果與數(shù)值模擬結(jié)果存在5%～25%的誤差，其產(chǎn)生的原因有2 個(gè)，一是在實(shí)驗(yàn)過(guò)程中循環(huán)水槽存在槽體共振現(xiàn)象。而實(shí)驗(yàn)中阻力的最大值只有4N 左右，所以輕微的槽體振動(dòng)也會(huì)對(duì)實(shí)驗(yàn)結(jié)果產(chǎn)生一定的影響。第2 個(gè)原因是本文數(shù)值模擬采用的是壁面固定式圓柱，而實(shí)驗(yàn)則是采用的“倒立入水”式。為了盡量減小自由液面對(duì)實(shí)驗(yàn)結(jié)果影響，實(shí)驗(yàn)中在與自由液面相切的位置增加了一塊底板，防止液面波動(dòng)對(duì)圓柱受力產(chǎn)生影響，但仍不能完全避免自由液面的影響。

圖5 雷諾數(shù)與阻力系數(shù)的關(guān)系Fig.5 Relationship between Re and drag coefficient

3.2 流場(chǎng)分析

流速U=0.6 m/s 工況下，長(zhǎng)徑比AR=3，6 和10 的三維流線(xiàn)圖如圖6 所示。長(zhǎng)徑比AR=3 的圓柱在展向幾乎沒(méi)有渦流，也沒(méi)有漩渦脫落。而長(zhǎng)徑比AR=6 和10 的圓柱在展向則有明顯的渦流以及漩渦脫落現(xiàn)象。同時(shí)發(fā)現(xiàn)，當(dāng)長(zhǎng)徑比AR=3 時(shí)，圓柱后方有極其混亂的湍流，其影響范圍從圓柱自由端一直到固定壁面。分析原因應(yīng)該是當(dāng)流體流經(jīng)圓柱自由端時(shí)產(chǎn)生“下洗”現(xiàn)象，當(dāng)長(zhǎng)徑比較小時(shí)，“下洗”運(yùn)動(dòng)直接沖擊到圓柱的根部，同時(shí)也破壞了圓柱后方有規(guī)律的漩渦脫落，產(chǎn)生了極其混亂的湍流；當(dāng)長(zhǎng)徑比AR為6 和10 時(shí)，圓柱后方均有一個(gè)較小的回流區(qū)，并且也會(huì)產(chǎn)生“下洗”現(xiàn)象，但因?yàn)殚L(zhǎng)徑比較大，“下洗”運(yùn)動(dòng)只能影響圓柱的自由端，破壞自由端附近的漩渦脫落，而對(duì)圓柱中間以及根部沒(méi)有影響。而且相比固定壁面處，圓柱中間部分的漩渦脫落更為明顯，說(shuō)明固定壁面的存在一定程度上抑制了漩渦脫落的產(chǎn)生。即圓柱的長(zhǎng)徑比越長(zhǎng)，“下洗”運(yùn)動(dòng)影響的長(zhǎng)度占展向長(zhǎng)度的比例就越小，所以隨著長(zhǎng)徑比的增加，阻力系數(shù)會(huì)緩慢增加。而有限長(zhǎng)圓柱繞流的阻力系數(shù)比無(wú)限長(zhǎng)圓柱低的主要原因同樣可能是有限長(zhǎng)圓柱自由端的存在使得流體可以繞過(guò)自由端流動(dòng)，而在無(wú)限長(zhǎng)的情況下流體只能從兩側(cè)繞著圓柱流動(dòng)。

圖6 不同長(zhǎng)徑比三維流線(xiàn)圖Fig.6 3D streamline with different AR

為了能清晰地得到圓柱后方漩渦脫落的流動(dòng)，采用Q準(zhǔn)則（Q-criterion）對(duì)尾渦進(jìn)行可視化。Q的定義為：

不同長(zhǎng)徑比下三維渦對(duì)比如圖7 所示?？梢钥闯觯?dāng)長(zhǎng)徑比AR=3 時(shí)，漩渦結(jié)構(gòu)沒(méi)有震蕩，因?yàn)槠溟L(zhǎng)徑比較小，所以自由端產(chǎn)生的高渦量控制了整個(gè)圓柱體的流動(dòng)行為；當(dāng)長(zhǎng)徑比AR=6 和10 時(shí)，出現(xiàn)了不同的漩渦脫落區(qū)域；漩渦在圓柱的中下部周期性震蕩，并且?guī)缀跖c圓柱平行；在圓柱的中上部，漩渦的周期性震蕩仍然存在，但是不再與圓柱平行，而是呈現(xiàn)出一定的角度。這種現(xiàn)象可能是由圓柱中上部分不同展向位置處的振幅不同引起的。長(zhǎng)徑比AR=10 時(shí)不同流速下三維渦對(duì)比如圖8 所示。隨著流速的增加，渦黏度最大值也在增加。當(dāng)從0.2 m/s 增加到0.4 m/s時(shí)，渦黏度最大值增加了92.4%，從0.4 m/s 增加到0.6 m/s時(shí)，渦黏度最大值增加了12.4%。說(shuō)明隨著流速的增加，渦黏度最大值的增速在逐漸降低。從圖中也可以看出，隨著流速的增加，圓柱后方不同位置處渦的相互作用也增強(qiáng)了。

圖7 不同長(zhǎng)徑比下三維渦粘度圖（Q=7）Fig.7 3D Eddy viscosity with different AR（Q=7）

圖8 長(zhǎng)徑比AR=10 時(shí)不同流速下三維渦粘度圖（Q=7）Fig.8 3D Eddy viscosity at different velocity when AR=10（Q=7）

4 結(jié) 語(yǔ)

本文對(duì)雷諾數(shù)Re=1×104-3×104和長(zhǎng)徑比AR=3，6 和10 的三維有限長(zhǎng)圓柱繞流進(jìn)行了數(shù)值模擬和實(shí)驗(yàn)研究，得到如下結(jié)論：

1）在所涉及的工況范圍內(nèi)，采用k-ω模型對(duì)三維有限長(zhǎng)圓柱繞流數(shù)值分析是可行的，計(jì)算穩(wěn)定，耗時(shí)少，計(jì)算誤差在可以接受范圍內(nèi)。

2）當(dāng)長(zhǎng)徑比AR=3，6 和10 時(shí)，數(shù)值模擬和實(shí)驗(yàn)的結(jié)果均表明，阻力系數(shù)隨流速的增加有一個(gè)逐漸減小的趨勢(shì)，數(shù)值模擬結(jié)果穩(wěn)定在0.8 左右，實(shí)驗(yàn)分析得到的阻力系數(shù)落在0.61～0.93 區(qū)間內(nèi)。表明在長(zhǎng)徑比大于3 的工況下，有限長(zhǎng)圓柱繞流的阻力系數(shù)趨于平穩(wěn)。

3）圓柱自由端的存在使流動(dòng)具有更明顯的三維特性，自由端的“下洗”作用會(huì)繞過(guò)自由端沖擊圓柱后方的漩渦脫落，導(dǎo)致阻力系數(shù)小于無(wú)限長(zhǎng)圓柱。在長(zhǎng)徑比AR=6 和10 時(shí)，出現(xiàn)了不同的漩渦脫落區(qū)域。漩渦在圓柱的中下部周期性震蕩，并且?guī)缀跖c圓柱平行。在圓柱的中上部，漩渦的周期性震蕩仍然存在，但是不再與圓柱平行，而是呈現(xiàn)出一定的角度。

4）渦黏度最大值隨著流速的增加也在增加。當(dāng)從0.2 m/s 增加到0.4 m/s 時(shí)，渦黏度最大值增加了92.4%，從0.4 m/s 增加到0.6 m/s 時(shí)，渦黏度最大值增加了12.4%。