串列雙方柱繞流的數(shù)值模擬

呂啟兵，楊 斌，楊忠超，李鵬浩

(1.重慶交通大學(xué)河海學(xué)院，重慶400074;2.溫州港集團(tuán)有限公司，浙江溫州325100)

鈍體繞流問(wèn)題廣泛存在于橋梁、海洋港口、航空航天、交通運(yùn)輸?shù)裙こ炭茖W(xué)領(lǐng)域，如置于水中的橋墩、石油鉆井平臺(tái)的墩柱、港口碼頭的樁柱、飛機(jī)機(jī)翼繞流、汽車(chē)的風(fēng)阻等。而方柱作為一種典型的鈍體，一直是流體力學(xué)領(lǐng)域繞流問(wèn)題的熱點(diǎn)之一。關(guān)于方柱繞流問(wèn)題的實(shí)驗(yàn)研究也已經(jīng)進(jìn)行得相當(dāng)廣泛，許多研究者通過(guò)各種方法對(duì)方柱的繞流進(jìn)行了實(shí)驗(yàn)?zāi)M研究，D.A.Lyn，等［1-2］提供了大量的實(shí)驗(yàn)?zāi)M資料，為實(shí)際的工程問(wèn)題提供了一定的參考。

隨著計(jì)算機(jī)計(jì)算能力的提升、數(shù)值算法的改進(jìn)與計(jì)算流體力學(xué)(CFD)的突飛猛進(jìn)，數(shù)值模擬技術(shù)已成為鈍體繞流問(wèn)題研究的一種有效手段，樁柱繞流問(wèn)題的數(shù)值模擬也得到了很大的發(fā)展。早期的Hiromasa Kawai，等［3］對(duì)雷諾數(shù)為 200 的串列方柱的繞流用差分方法做了數(shù)值模擬;王遠(yuǎn)成，等［4］人采用湍流RNG k-ε模型對(duì)鈍體繞流流場(chǎng)的不穩(wěn)定、非定常流動(dòng)進(jìn)行了數(shù)值模擬;S.Muralami，等［5］使用標(biāo)準(zhǔn)k-ε模型對(duì)鈍體繞流進(jìn)行了研究，并且利用了大渦模型(LES)進(jìn)行了數(shù)值計(jì)算;陳素琴，等［6］數(shù)值模擬了雷諾數(shù)為10 000情況下兩串列方柱繞流的干擾問(wèn)題，并分析了不同間距下各方柱的升阻力特征及流場(chǎng)特性。

由于內(nèi)河特殊的水文環(huán)境(大水位差變化與較大的流速)，近年來(lái)內(nèi)河碼頭運(yùn)用全直柱框架結(jié)構(gòu)作為一種新的碼頭結(jié)構(gòu)型式，在考慮施工過(guò)程、工程條件與經(jīng)濟(jì)因素等方面時(shí)，全直樁方柱框架碼頭已十分普遍。筆者利用流體軟件、采用RNG k-ε模型對(duì)串聯(lián)雙方柱繞流進(jìn)行了數(shù)值模擬，分析了每根柱升、阻力特性與不同樁距下的流場(chǎng)演變情況，為實(shí)際工程問(wèn)題提供一定參考。

1 數(shù)值模型與模型驗(yàn)證

1.1 控制方程

對(duì)于不可壓縮黏性流體，在直角坐標(biāo)系下，流體的運(yùn)動(dòng)可用Navier-Stokes方程［7］與動(dòng)量方程描述：

模擬計(jì)算采用RNG k-ε模型，該模型由Yakhot和Orzag提出，模型是從瞬態(tài)N-S方程中推出的，通過(guò)大尺度運(yùn)動(dòng)和修正后的黏度項(xiàng)體現(xiàn)小尺度的影響，而使這些小尺度運(yùn)動(dòng)有系統(tǒng)地從控制方程中去除，所得到的k方程和ε方程與標(biāo)準(zhǔn)k-ε模型非常相似：

式中：Gk是由層流速度梯度產(chǎn)生的湍流動(dòng)能k的產(chǎn)生項(xiàng)C1ε=1.44，η。=4.377，β =0.012，C2ε=1.68;ε 為湍動(dòng)耗散率。

RNG k-ε模型來(lái)源于嚴(yán)格的統(tǒng)計(jì)技術(shù)，是對(duì)標(biāo)準(zhǔn)k-ε模型的重要改進(jìn)。RNG模型在ε方程中加入了一個(gè)Rε條件，從而反映了主流的時(shí)均應(yīng)變率Eij，這樣該模型中產(chǎn)生項(xiàng)不僅與流動(dòng)情況有關(guān)，而且在同一問(wèn)題中也還是空間坐標(biāo)的函數(shù)，從而可以更好的處理高應(yīng)變率及流線彎曲程度較大的流動(dòng)。

1.2 數(shù)值算法

計(jì)算區(qū)域網(wǎng)格采用分塊結(jié)構(gòu)化網(wǎng)格(圖1)，網(wǎng)格總量在4.9萬(wàn)～18.9萬(wàn)之間，利用有限體積法(FVM)對(duì)方程進(jìn)行離散，壓力基求解器與壓力2階離散格式，動(dòng)量、湍動(dòng)能和湍動(dòng)耗散率選項(xiàng)全部為2階迎風(fēng)格式，計(jì)算模型選用RNG k-ε模型，PISO算法。

1.3 邊界條件

進(jìn)口邊界條件：u=1，v=0。

出口邊界條件：自由出流。

上下邊界及方柱表面為壁面邊界條件，即u=0，v=0的無(wú)滑移邊界條件。

雷諾數(shù)：Re=UD/υ=22 000。式中：U為來(lái)流流速;D為方柱邊長(zhǎng);υ為流體的動(dòng)力黏性系數(shù)。

1.4 計(jì)算區(qū)域與網(wǎng)格劃分

計(jì)算區(qū)域如圖1，單柱計(jì)算區(qū)域?yàn)?2D×13D(D為方柱邊長(zhǎng))，距入口邊界距離6D，考慮上下邊界對(duì)方柱周?chē)鲌?chǎng)的影響，上下邊界距方柱都為6D，與出口距離為15D;雙方柱間距為L(zhǎng)，其余設(shè)置與單方柱相同。

圖1 單方柱與雙方柱計(jì)算區(qū)域與網(wǎng)格Fig.1 Calculation area and grid of single square cylinder and double square cylinder

1.5 模型驗(yàn)證

方柱繞流常用于模擬研究流體的荷載，即對(duì)建筑物和構(gòu)筑物的阻力、升力和柱上渦街脫落頻率。為了驗(yàn)證數(shù)值模擬的準(zhǔn)確性，首先計(jì)算了Re=22 000情況下的單方柱繞流問(wèn)題。方柱上阻力系數(shù)、升力系數(shù)與斯特羅哈數(shù)的定義為：

式中：Fd為阻力，F(xiàn)l為升力，u為來(lái)流流速，D為方柱邊長(zhǎng)，f為渦街脫落頻率，即方柱上渦街的脫落頻率，、Cl分別代表阻力系數(shù)均值與升力系數(shù)。

從表1可以看出，模擬計(jì)算結(jié)果與 D.A.Lyn，等［1-2］試驗(yàn)結(jié)果非常接近，模擬計(jì)算誤差在7%以內(nèi)，模擬得到的單方柱情形下的升、阻力系數(shù)如圖2。

表1 單方柱繞流模擬計(jì)算結(jié)果比較Table 1 Comparison results of flow around square cylinder simulation

圖2 Re=22 000時(shí)單方柱升、阻力系數(shù)Fig.2 Lift and drag coefficients of single square cylinder when Re=22 000

2 模擬結(jié)果及分析

2.1 升、阻力分析

筆者主要關(guān)注不同間距情況下各方柱的阻力、升力系數(shù)與變化情況，簡(jiǎn)略分析了流場(chǎng)演變情況。

表2給出了不同間距下上游方柱的阻力、升力與斯特勞哈數(shù)。在L/D=1時(shí)，兩方柱的升力都非常的小，表明兩方柱此時(shí)的脈動(dòng)微弱;在L/D=2時(shí)，上游方柱阻力較其他情況時(shí)有所減小，說(shuō)明間距較小時(shí)，兩柱存在干擾作用。隨著間距比的不斷增大，上游方柱阻力(圖3)與斯特勞哈數(shù)St趨于單柱繞流情況下的情形;在L/D=3時(shí)上下方柱的升力都變化較大，此種間距比下游方柱增加的升力有較大增加，主要是因?yàn)椋荷嫌畏街奈擦鲗?duì)下游方柱的來(lái)流有很大影響，此間距下上游方柱的尾渦已開(kāi)始脫落，擴(kuò)散出去，改變了下游方柱的來(lái)流(圖4)，導(dǎo)致上下游方柱的升力都有大幅度的增加。

表2 上下游方柱阻力計(jì)算結(jié)果Table 2 Calculation results of drags of square cylinders at upstream and downstream

圖3 上下游方柱阻力與升力系數(shù)走勢(shì)Fig.3 Trend of drag and lift coefficient of upstream and downstream square cylinder

圖4 流線圖局部放大Fig.4 Local amplification of flow line

從表2可以看出，在L/D=1與2間距下時(shí)下游方柱阻力為負(fù)值，此時(shí)方柱受力為水流的相反方向，其主要原因是：當(dāng)兩方柱間距較小時(shí)(如L/D=1與2)，下游方柱受上游方柱尾流漩渦影響，柱前形成為負(fù)壓，從而導(dǎo)致這樣的結(jié)果，此結(jié)果與一般的實(shí)驗(yàn)情況基本吻合。此時(shí)兩方柱的升力系數(shù)非常小，斯特勞哈數(shù)St較單柱情況時(shí)變小，說(shuō)明方柱的脈動(dòng)較小，而渦的形成頻率較單柱時(shí)加快。隨后方柱脈動(dòng)開(kāi)始變大，上游方柱阻力有所減小，但減小得并不十分明顯，而后兩方柱的脈動(dòng)放緩。兩方柱的間距比L/D值的變化對(duì)上游方柱所受阻力影響不大，主要為上游方柱對(duì)下游方柱的影響;當(dāng)L/D≥5時(shí)，上游方柱的所受阻力、升力、斯特勞哈數(shù)均與單柱情況幾乎一致。隨著方柱間距的增大，下游方柱所受阻力呈變大的趨勢(shì)，升力較單柱時(shí)大了近一倍左右。其原因?yàn)橄掠畏街挥谏嫌畏街蛛x的脈動(dòng)高速區(qū)附近，上游方柱脫落的渦直接作用于下游方柱，產(chǎn)生更小尺度的渦，在其前角產(chǎn)生額外的角動(dòng)量，從而產(chǎn)生這樣的結(jié)果。

當(dāng)兩方柱靠得較近時(shí)，阻力與升力變化都比較劇烈，其變化的趨勢(shì)各不相同，當(dāng)距離較遠(yuǎn)時(shí)，兩方柱的阻力趨于一致，升力也趨于穩(wěn)定，這在工程設(shè)計(jì)建設(shè)中對(duì)水流荷載的考慮提供一定參考價(jià)值。

2.2 流場(chǎng)演變

圖5給出5種情形下的瞬時(shí)流線圖。間距較小時(shí)，只在方柱背后產(chǎn)生漩渦，兩方柱附近的流場(chǎng)改變不大，尾流區(qū)較遠(yuǎn)處沒(méi)有渦的擺動(dòng);隨著L的增大，下游方柱前后流態(tài)變化明顯，尾流區(qū)形成周期性的渦街，兩柱都有規(guī)律性渦脫落，并且下游方柱渦街?jǐn)[動(dòng)越來(lái)越大，說(shuō)明L的增大增加了兩方柱尾流的不穩(wěn)定性，而上游方柱周?chē)鲌?chǎng)幾乎不受下游方柱影響。從圖中對(duì)比可以看出，串列雙方柱情形下流場(chǎng)只有在L/D=1間距比下才出現(xiàn)了對(duì)稱情況，并伴有對(duì)稱尾流，其余情況主要為偏轉(zhuǎn)尾流。在L/D=2時(shí)兩柱之間有較大回流區(qū)，隨著L的增大，回流區(qū)有一定的減弱。

圖5 流線圖Fig.5 Flow line

圖6給出了 L/D 為1，2，3，6，10情況下的渦量圖。可以發(fā)現(xiàn)，隨著間距的增大，渦量較大的地方由方柱兩側(cè)移向方柱的后方，而后方柱卻剛好相反。

圖6 渦量圖Fig.6 The vorticity

3 結(jié)論

筆者主要研究了兩串列方柱擾流的干擾情況。模擬計(jì)算過(guò)程中發(fā)現(xiàn)，兩方柱間網(wǎng)格精度尤為重要，直接關(guān)系到下游方柱繞流特性模擬所得的數(shù)值。經(jīng)研究得到以下結(jié)論：

1)當(dāng)兩方柱串列排列時(shí)，上游方柱所受阻力幾乎無(wú)變化，與單方柱情形時(shí)大致相同。

2)下游方柱所受阻力在隨著間距的增大而先減小后增大，且均小于上游方柱所受阻力，屏蔽效應(yīng)明顯;由于上游方柱的存在，下游方柱升力明顯增大。

3)間距較大時(shí)，上下游方柱的脈動(dòng)頻率幾乎相等。

4)在串聯(lián)時(shí)兩方柱之間間距L的變化對(duì)流場(chǎng)的改變作用明顯。

［1］Lyn D A，Einav S，Rodi W.Laser doppler velocimetry study of ensemble-averaged characteristic of the turbulent flow near of a square cylinder［J］.Fluid Mechanics，1995，304：285-319.

［2］Lyn D A，Rodi W.The flapping shear layer formed by flow separation from the forward corner of a square cylinder［J］.Fluid Mesh，1994，267：353-376.

［3］Hiromasa Kawai，Kiyoshi Fujinamic.Numerical Simulation of Flow around Square Prisms in Tandem Arrangements［C］//9thICME.Houston Texas：Gulf Publishing，1995：185-186.

［4］王遠(yuǎn)成，吳文權(quán).方柱繞流流場(chǎng)的RNG方法模擬研究［J］.水動(dòng)力學(xué)研究與進(jìn)展：A，2004，19(增刊1)：916-920.Wang Yuancheng，Wu Wenquan.Numerical simulation of flow around square cylinder using RNG k-ε turbulence model［J］.Journal Hydrodynamics：A，2004，19(S1)：916-920.

［5］Muralami S，Mochida A.3-D numerical simulation of airflow around a cubic model by means of the k-ε model［J］.Journal of Wind Engineering and Industrial Aerodynamics，1988，31(1)：327-338.

［6］陳素琴，黃自萍，沈劍華，等.兩串列繞流的干擾數(shù)值模擬研究［J］.同濟(jì)大學(xué)學(xué)報(bào)：自然科學(xué)版，2001，29(3)：320-325.Chen Suqin，Huang Ziping，Shen Jianhua，et al.Numerical computation of the flow around two square cylinders in tandem arrangement［J］.Journal of Tongji University：Natural Science，2001，29(3)：320-325.

［7］張遠(yuǎn)君.流體力學(xué)大全［M］.北京：北京航空航天大學(xué)出版社，1991.Zhang Yuanjun.Fluid Mechanics［M］.Beijing：Beijing University of Aeronautics and Astronautics Press，1991.

［8］Rodi W.Simulation of flow past buildings with statistical turbulence models［J］.Wind Climate in Cities，1995，277：649-668.