梁圣召，周軍偉，梅 蕾

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

0 引 言

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

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

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

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

1 數(shù)值模擬方法

1.1 控制方程

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

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

本文采用Wilcox 的k-ω湍流模型，具體可參照文獻[9]?；贑FX 求解器實現(xiàn)方程的求解。

1.2 計算域模型

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

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

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

1.3 數(shù)值方法驗證

1）分別進行網(wǎng)格數(shù)量和湍流模型對計算結(jié)果和計算時間的對比分析。計算條件為：流場無窮遠處的均勻來流速度為U=0.6 m/s，流體為水，密度ρ=1 000 kg/m3，運動粘性系數(shù)ν=10-6m2/s，計算得雷諾數(shù)Re=3×104，時間步長0.001 s。網(wǎng)格在圓柱周向均布100 個節(jié)點，法向擠出20 層，線性增長率1.3。近壁面處保證y+＜2。3 種網(wǎng)格密度采用同樣的加密區(qū)域7D×3D，采用四面體、棱柱型與六面體網(wǎng)格相結(jié)合的方式，以改善網(wǎng)格質(zhì)量，計算域網(wǎng)格劃分如圖2 所示。3 套網(wǎng)格分別命名為Coarse，Medium 和Fine，網(wǎng)格數(shù)分別為49 萬、119 萬和200 萬。表1 為k-ω湍流模型不同網(wǎng)格密度條件下，以及在Medium 網(wǎng)格下分別采用DES，k-ω和k-ε三種湍流模型所得到的升阻力系數(shù)以及計算所耗時間。可以看出當(dāng)采用k-ω模型時，使用Medium 和Fine 網(wǎng)格，Cd計算結(jié)果相差2.7%，Cl計算結(jié)果相差9.93%，但是計算時間相差1 倍，所以Medium 網(wǎng)格數(shù)計算精度足夠。而在Medium 網(wǎng)格數(shù)下，k-ε的計算結(jié)果與DES 和k-ω的結(jié)果相差較大，DES 和k-ω的計算精度最大誤差在15%以內(nèi)，但是計算時間相差近1 倍。綜合考慮計算精度與計算時間成本，本文選取Medium 網(wǎng)格數(shù)119 萬，湍流模型選擇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）鑒于文獻在研究相關(guān)問題時多采用大渦模擬（LES）湍流模型，為了進一步驗證本文計算模型的可用性，在同等初始條件下，對不同工況及長徑比的有限長圓柱進行計算，并與前人所做相近工況的實驗和數(shù)值模擬結(jié)果進行對比。首先，本文計算和實驗結(jié)果與前人所做相近工況的實驗和數(shù)值模擬結(jié)果進行對比，如表2 所示。可以看出：在長徑比AR=6，雷諾數(shù)Re=3×104工況下，本文計算結(jié)果比H.Sakamoto[10]的實驗值大6.59%，實驗結(jié)果比H.Sakamoto[10]的實驗值小12.3%；在長徑比AR=10，雷諾數(shù)Re=3×104工況下本文的計算結(jié)果比LEE[5]的計算值大5.18%，實驗值比其小2.56%。與相似工況下的文獻結(jié)果的誤差均在可以接受范圍內(nèi)。同時，在長徑比AR=10，雷諾數(shù)Re=3×104工況下，本文的實驗結(jié)果與計算結(jié)果也進行了對比，如圖3 所示?？芍獙嶒炛蹬c計算值的升阻力系數(shù)時程圖比較接近，而橫流向的振動頻率約是順流向振動頻率的2 倍，這與王曉聰[11]的研究結(jié)論相符，也側(cè)面佐證了本文研究的準(zhǔn)確性。最終選取k-ω湍流模型進行數(shù)值分析。

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

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

2 實驗測試

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

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

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

3 結(jié)果與分析

3.1 阻力系數(shù)分析

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

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

3.2 流場分析

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

圖6 不同長徑比三維流線圖Fig.6 3D streamline with different AR

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

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

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

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

4 結(jié) 語

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

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

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

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

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