王成彥，巫凱旋，陳沛楷，黃 技

（廣東海洋大學(xué)船舶與海運(yùn)學(xué)院，廣東湛江 524005）

0 引 言

圓柱繞流是流體力學(xué)研究領(lǐng)域的一個(gè)經(jīng)典問題，在海洋工程和土木工程等領(lǐng)域中有著廣泛應(yīng)用［1］。近年來(lái)，隨著陸地上的石油和天然氣資源日益枯竭，石油開采不斷向深海發(fā)展，多柱體系統(tǒng)在海洋平臺(tái)及海底管線中得到了廣泛應(yīng)用［2-3］。方形排列的四圓柱是海洋工程中的一種典型結(jié)構(gòu)，廣泛存在于各種海洋工程結(jié)構(gòu)物中，如延伸式張力腿平臺(tái)。柱體繞流產(chǎn)生的疲勞破壞可能會(huì)威脅人員的安全，造成巨大的財(cái)產(chǎn)損失，危害海洋環(huán)境，因此研究四圓柱的擾流現(xiàn)象具有重要意義。

近年來(lái)，國(guó)內(nèi)外學(xué)者針對(duì)圓柱繞流問題開展了大量研究。隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，數(shù)值模擬手段成本低、實(shí)施方便的優(yōu)點(diǎn)逐漸顯現(xiàn)，各種數(shù)值方法得到了較快發(fā)展。韓兆龍等［4］采用譜單元法對(duì)多圓柱繞流進(jìn)行了研究；龔志軍等［5］采用格子Boltzmann法對(duì)低雷諾數(shù)下四圓柱的繞流模式及受力特性進(jìn)行了模擬；賈曉荷等［6］采用大渦模擬的方法對(duì)采用不同排列方式的雙圓柱的繞流進(jìn)行了仿真和計(jì)算；謝啟迪等［7］在OpenFOAM環(huán)境下采用RANS方法對(duì)不同雷諾數(shù)下的單圓柱繞流進(jìn)行了分析；WANG等［8］從試驗(yàn)的角度出發(fā)，對(duì)雷諾數(shù)為8 000 時(shí)的方形排列四圓柱繞流進(jìn)行了物理實(shí)驗(yàn)；殷長(zhǎng)山［9］采用粒子圖像測(cè)速（Particle Image Velocimetry，PIV）技術(shù)對(duì)不同迎流角下的四圓柱繞流進(jìn)行了試驗(yàn)研究；LAM等［10］對(duì)不同間距比（1.6 ～5.0）下雷諾數(shù)為200 時(shí)的四圓柱繞流現(xiàn)象進(jìn)行了研究，得到遮蔽模式轉(zhuǎn)化為再附著模式的臨界間距比為2.0，再附著模式轉(zhuǎn)化為漩渦拍擊模式的臨界間距比為3.5。

雖然已有學(xué)者［9，11］對(duì)不同迎流角下的方形排列四圓柱進(jìn)行研究，但缺少對(duì)不同迎流角下的方形排列四圓柱更全面、系統(tǒng)的研究。同時(shí)，LAM等［10］和于定勇等［12］在研究間距比對(duì)方形排列四圓柱的繞流影響時(shí)均指出尾流模式轉(zhuǎn)換的臨界間距比為2.0 和3.5，因此對(duì)臨界間距比和迎流角進(jìn)行研究非常有意義。

真實(shí)的流場(chǎng)大多是在超高雷諾數(shù)下出現(xiàn)的，采用大渦模擬（Large Eddy Simulation，LES）或直接數(shù)值模擬（Direct Numerical Simulation，DNS）方法會(huì)消耗大量計(jì)算資源，對(duì)計(jì)算機(jī)的計(jì)算能力有很高的要求。本文的研究是在雷諾數(shù)為3 900 的情況下開展的，屬于亞臨界區(qū)，盡管該雷諾數(shù)與實(shí)際雷諾數(shù)相差較大，但超高雷諾數(shù)下湍流場(chǎng)的特征往往也能在較低的雷諾數(shù)下得到體現(xiàn)，因此該研究可供實(shí)際工程及以后更高雷諾數(shù)的研究參考。本文采用計(jì)算流體動(dòng)力學(xué)（Computational Fluid Dynamics，CFD）軟件Fluent，重點(diǎn)分析來(lái)流角度為0°、15°、30°和45°，間距比分別為2.0 和3.5 時(shí)四圓柱的水動(dòng)力特性。

1 理論基礎(chǔ)

1.1 流體控制方程

根據(jù)流體力學(xué)理論，本文所研究的二維模型的連續(xù)性方程為

式（1）中：ui，uj，uk分別為x，y，z三個(gè)方向上流體的速度分量，ρ為流體的密度。

假設(shè)流場(chǎng)內(nèi)的溫度不變，則只考慮動(dòng)量方程

式（2）中：u，p，υ分別為速度矢量、壓力和運(yùn)動(dòng)粘度；?為哈密頓算子

1.2 湍流模型

本文研究的工況為亞臨界區(qū)（雷諾數(shù)Re=3 900），湍流模型選用RNGκ-ε湍流模型。該模型是采用基于重整化群理論的統(tǒng)計(jì)方法推導(dǎo)得出的，其形式類似于standard κ-ε模型，但在ε方程中增加了1 項(xiàng)，提高了高速流動(dòng)的準(zhǔn)確性，考慮了渦流對(duì)湍流的影響，提升了漩渦流動(dòng)的精度，同時(shí)為Prandtl 數(shù)提供了一個(gè)解析公式，而standard κ-ε模型只能由用戶指定常值。

1）湍流動(dòng)能κ輸運(yùn)方程為

2）湍流耗散率ε輸運(yùn)方程為

式（3）～式（5）中：μeff=μ +μt；η =Sκ/ε；S 為應(yīng)變率，S =（2Eij·Eij）1/2；0.012；Gk為由層流速度梯度而產(chǎn)生的湍流能,αε分別為κ方程和ε方程的湍流Prandtl數(shù)。

2 計(jì)算模型及算例驗(yàn)證

2.1 計(jì)算模型

本文共進(jìn)行8 個(gè)組次的數(shù)值模擬，各組次工況設(shè)置見表1。所設(shè)置的流體域大小為24D×24D（D為圓柱直徑），2 個(gè)相鄰圓柱中心之間的距離為L(zhǎng)，Re=3 900，4 個(gè)圓柱呈方形排列。具體計(jì)算區(qū)域見圖1。

表1 各組次工況設(shè)置

圖1 計(jì)算區(qū)域

本文的壓力與速度耦合采用SIMPLE算法實(shí)現(xiàn)，壓力計(jì)算格式采用二階格式，動(dòng)量計(jì)算格式、湍動(dòng)能和湍流耗散率均采用二階迎風(fēng)格式。瞬態(tài)求解采用二階隱式格式實(shí)現(xiàn)。整個(gè)流體域采用結(jié)構(gòu)化網(wǎng)格，并在圓柱表面設(shè)置邊界層網(wǎng)格，以提高計(jì)算精度。對(duì)圓柱和尾流區(qū)域進(jìn)行局部加密，以得到較精確的尾流形態(tài)。整體網(wǎng)格劃分和圓柱邊界層網(wǎng)格示意見圖2 和圖3。

圖3 圓柱邊界層網(wǎng)格示意

2.2 邊界條件與初始條件

1）邊界條件：入口邊界條件為速度入口；出口邊界條件為壓力出口；上下邊界設(shè)為對(duì)稱邊界；圓柱表面為無(wú)滑移壁面。

2）初始條件：為節(jié)省計(jì)算時(shí)間和提高計(jì)算精度，首先進(jìn)行500 步穩(wěn)態(tài)模擬，并將其作為初始條件，隨后進(jìn)行瞬態(tài)計(jì)算。

3 數(shù)值計(jì)算方法驗(yàn)證

為驗(yàn)證湍流模型、網(wǎng)格劃分和參數(shù)設(shè)置的準(zhǔn)確性，首先在Re=3 900 的情況下對(duì)單圓柱繞流進(jìn)行模擬，并將所得結(jié)果與其他研究結(jié)果相對(duì)比，結(jié)果見表2。隨后對(duì)Re=3 900、間距比為3.5 的四圓柱繞流過(guò)程進(jìn)行模擬，并將仿真結(jié)果與其他研究結(jié)果相對(duì)比，具體見表3。

表2 Re=3 900 情況下的單圓柱結(jié)果驗(yàn)證

表3 四圓柱結(jié)果驗(yàn)證

由表2 和表3 可知，本文得到的單圓柱斯特魯哈數(shù)與已有研究結(jié)果較為接近，而得到的阻力系數(shù)與已有研究結(jié)果相比較小，這是由于雷諾平均方程只能提供湍流的平均信息，忽略了湍流脈動(dòng)量。通過(guò)對(duì)四圓柱的驗(yàn)證結(jié)果進(jìn)行分析可知，1 號(hào)圓柱的阻力系數(shù)與文獻(xiàn)［12］的研究結(jié)果較為接近，3 號(hào)圓柱的斯特魯哈數(shù)與其他文獻(xiàn)的研究成果均吻合得較好。因此，可驗(yàn)證本文的網(wǎng)格劃分、湍流模型選擇和參數(shù)設(shè)置是較為準(zhǔn)確的，可進(jìn)行下一步的研究。

4 計(jì)算結(jié)果及分析

4.1 流場(chǎng)速度發(fā)展云圖

當(dāng)Re=3 900，間距比為2.0 和3.5 時(shí)，方形排列四圓柱在不同來(lái)流角度下的流場(chǎng)速度發(fā)展云圖見圖4。由于來(lái)流角度不同，流場(chǎng)速度云圖存在明顯差異。當(dāng)間距比為2.0，來(lái)流角度不為0°時(shí)，偏斜尾流現(xiàn)象非常明顯，上游圓柱沒有形成明顯的尾渦，下游圓柱脫落的尾渦被延長(zhǎng)，四圓柱的尾跡相比來(lái)流角度為0°時(shí)更明顯地被混合在一起。

圖4 方形排列四圓柱在不同來(lái)流角度下的流場(chǎng)速度發(fā)展云圖

當(dāng)間距比為3.5 時(shí)，隨著間距比的增加，上游圓柱脫落的尾渦“拍擊”在下游圓柱上。當(dāng)來(lái)流角度為15°和30°時(shí)，尾流區(qū)域流場(chǎng)更紊亂，沒有對(duì)稱性，尤其是當(dāng)來(lái)流角度為30°時(shí)表現(xiàn)得更強(qiáng)烈。當(dāng)來(lái)流角度為45°時(shí)，2 號(hào)圓柱與4 號(hào)圓柱的尾渦基本對(duì)稱。由此可見，來(lái)流角度不同會(huì)極大地影響流場(chǎng)和尾渦的形態(tài)。

4.2 升力系數(shù)和阻力系數(shù)

在圓柱繞流過(guò)程中，周期性脫落的漩渦會(huì)對(duì)圓柱體產(chǎn)生繞流阻力Fd和繞流升力Fl，為方便研究，常進(jìn)行無(wú)量綱化處理，得到阻力系數(shù)Cd和升力系數(shù)Cl，分別表征圓柱橫向和縱向的受力變化。本文通過(guò)數(shù)值模擬得到臨界間距比（L/D=2，L/D=3.5）下，不同來(lái)流角度（α=0°，α=15°，α=30°，α=45°）對(duì)應(yīng)的各圓柱阻力系數(shù)和升力系數(shù)的變化情況。

1）當(dāng)間距比為2.0 和3.5 時(shí)，各圓柱的平均阻力系數(shù)隨來(lái)流角度的變化情況分別見圖5a和圖5b。

由圖5a可知，當(dāng)間距比為2.0 時(shí)，上游圓柱對(duì)下游圓柱的遮蔽作用較強(qiáng)，因此上游圓柱的阻力系數(shù)相比下游圓柱更大。但是，隨著來(lái)流角度的增大，2 號(hào)圓柱的阻力系數(shù)增大，超過(guò)了1 號(hào)圓柱，這是由于2 號(hào)圓柱沒有受到1 號(hào)圓柱的遮蔽作用，受到了其尾流的影響，且2 號(hào)圓柱的阻力系數(shù)在來(lái)流角度為15°時(shí)達(dá)到最大之后隨來(lái)流角度的增大逐漸減小，但仍保持在較高水平，這說(shuō)明隨著來(lái)流角度的增大，2 號(hào)圓柱實(shí)際受到的1號(hào)圓柱的尾流影響在減弱。3 號(hào)圓柱仍處于1 號(hào)圓柱的尾流遮蔽作用下，因此其阻力系數(shù)仍較小，在來(lái)流角度為30°時(shí)也較小。但是，當(dāng)來(lái)流角度增大到45°時(shí)，3 號(hào)圓柱的阻力系數(shù)顯著增大并超過(guò)1 號(hào)圓柱，與2 號(hào)圓柱的阻力系數(shù)接近。這說(shuō)明，當(dāng)下游圓柱擺脫上游圓柱的遮蔽作用之后，仍會(huì)受到上游圓柱尾流的影響，導(dǎo)致其阻力系數(shù)增大。4 號(hào)圓柱受到的遮蔽作用最強(qiáng)烈，因此其阻力系數(shù)最小，且隨著來(lái)流角度的改變，其阻力系數(shù)一直保持在較低的水平。由圖5b可知，當(dāng)間距比為3.5，來(lái)流角度為0°時(shí)，下游圓柱的阻力系數(shù)相比間距比為2.0 時(shí)的大，這說(shuō)明上游圓柱對(duì)下游圓柱仍有遮蔽作用，但是不如間距比為2.0 時(shí)強(qiáng)烈。當(dāng)來(lái)流角度為15°時(shí)，2 號(hào)圓柱和3 號(hào)圓柱的阻力系數(shù)與1 號(hào)圓柱接近，且當(dāng)來(lái)流角度繼續(xù)增大時(shí)，2 號(hào)圓柱和3 號(hào)圓柱的阻力系數(shù)增大至超過(guò)1 號(hào)圓柱。4 號(hào)圓柱的阻力系數(shù)隨來(lái)流角度的增大先增大后減小，由此可觀察到下游圓柱交替受到上游圓柱的遮蔽作用和尾流的影響。

圖5 各圓柱的平均阻力系數(shù)隨來(lái)流角度的變化情況

2）當(dāng)間距比為2.0 和3.5 時(shí)，各圓柱的平均升力系數(shù)隨來(lái)流角度的變化情況分別見圖6a和圖6b。

圖6 各圓柱的平均升力系數(shù)隨來(lái)流角度的變化情況

由圖6a可知：1 號(hào)圓柱的升力系數(shù)隨來(lái)流角度的增大不斷增大，2 號(hào)、3 號(hào)和4 號(hào)圓柱的升力系數(shù)均是在來(lái)流角度為15°時(shí)達(dá)到最大，隨后不斷減小；當(dāng)來(lái)流角度為45°時(shí)，1 號(hào)圓柱與4 號(hào)圓柱的升力系數(shù)十分接近，與單圓柱的升力系數(shù)也十分接近；2 號(hào)圓柱與3 號(hào)圓柱的升力系數(shù)更接近一種相反數(shù)的關(guān)系，在圖像上表現(xiàn)為關(guān)于1 號(hào)圓柱和4 號(hào)圓柱對(duì)稱。由此可知，2 號(hào)圓柱與3 號(hào)圓柱是交替向后泄放漩渦的，有一定的對(duì)稱關(guān)系。根據(jù)圖6b，1 號(hào)圓柱的升力系數(shù)保持逐漸增大的趨勢(shì)。2 號(hào)圓柱的升力系數(shù)隨來(lái)流角度的增大逐漸增大，并在來(lái)流角度為30°時(shí)達(dá)到最大，隨后逐漸減小。3 號(hào)圓柱和4 號(hào)圓柱的升力系數(shù)隨來(lái)流角度的增大逐漸減小。當(dāng)來(lái)流角度為45°時(shí)，各圓柱的升力系數(shù)與單圓柱的升力系數(shù)趨于相等，這表現(xiàn)出了其在間距比為3.5 時(shí)所具有的特點(diǎn)。

4.3 斯特魯哈數(shù)St

St實(shí)際上表現(xiàn)了圓柱無(wú)量綱化的旋渦脫落頻率［4］。通過(guò)對(duì)所得升力系數(shù)作快速傅里葉變換即可得到旋渦脫落頻率fx，再進(jìn)行無(wú)量綱化處理即可得到St。St與fx的關(guān)系為

式（6）中：D為圓柱直徑；U為來(lái)流速度；fx為旋渦脫落頻率。

本文采用Fluent軟件自帶的后處理功能，通過(guò)導(dǎo)入計(jì)算時(shí)所得升力系數(shù)即可獲得St。圖7a 和圖7b 分別為間距比為2.0 和3.5 時(shí)各圓柱的St隨來(lái)流角度的變化情況。

圖7 各圓柱的St隨來(lái)流角度的變化情況

由圖7a可知，當(dāng)間距比為2.0 時(shí)，隨著來(lái)流角度的增大，1 號(hào)、2 號(hào)和3 號(hào)圓柱的St均表現(xiàn)為先上升后下降，而4 號(hào)圓柱的St表現(xiàn)為逐漸下降。造成這種現(xiàn)象的原因是在來(lái)流角度增大過(guò)程中，上游圓柱對(duì)4 號(hào)圓柱的影響逐漸加大。由圖7b可知，雖然隨著來(lái)流角度的變化，St 也有波動(dòng)，但沒有發(fā)生較大的變化，且各圓柱的St均保持相同。這是由于方形排列的四圓柱的漩渦脫落由再附著模式變?yōu)殇鰷u拍擊模式的臨界間距比剛好為3.5［9］。因此，各圓柱的St是相等的，這與已有研究觀察到的現(xiàn)象相符。

5 結(jié) 語(yǔ)

本文采用有限體積法，利用Fluent軟件中的RNG κ-ε湍流模型，對(duì)臨界間距比、不同來(lái)流角度下四圓柱的繞流進(jìn)行數(shù)值模擬，通過(guò)分析各圓柱的阻力系數(shù)和升力系數(shù)，主要得到以下結(jié)論：

1）隨著來(lái)流角度的增大，2 號(hào)圓柱和3 號(hào)圓柱會(huì)逐漸擺脫上游圓柱的遮蔽作用，而受到上游圓柱引起的流場(chǎng)變化的作用，這2 個(gè)圓柱的阻力系數(shù)增大。4 號(hào)圓柱仍受到遮蔽作用的影響，其阻力系數(shù)偏小。

2）在來(lái)流角度增大過(guò)程中，升力系數(shù)呈現(xiàn)出由離散逐漸收斂的趨勢(shì)。當(dāng)間距比為2.0，來(lái)流角度為45°時(shí)，2 號(hào)圓柱和3 號(hào)圓柱的升力系數(shù)的幅值在圖像上表現(xiàn)為關(guān)于1 號(hào)圓柱和4 號(hào)圓柱對(duì)稱。當(dāng)間距比為3.5，來(lái)流角度為45°時(shí)，4 個(gè)圓柱的升力系數(shù)趨于相等，且與單圓柱的升力系數(shù)相近。

3）當(dāng)間距比為2.0 時(shí)，1 號(hào)、2 號(hào)和3 號(hào)圓柱的St均表現(xiàn)為先上升后下降，4 號(hào)圓柱的St一直下降。當(dāng)間距比為3.5 時(shí)，各圓柱的St相等，只是隨來(lái)流角度的變化而略有不同。