靜止方柱和圓柱繞流的二維數(shù)值分析

畢繼紅 余化軍 任洪鵬

（1.天津大學(xué) 建筑工程學(xué)院，天津 300072；2.濱海土木工程結(jié)構(gòu)與安全教育部重點(diǎn)實(shí)驗(yàn)室（天津大學(xué)），天津 300072）

在一定的雷諾數(shù)范圍內(nèi)，粘性流體流經(jīng)高層建筑物、橋塔、海上石油鉆井平臺(tái)、大跨度的懸索橋或斜拉橋、輸電線等具有鈍體截面形狀的結(jié)構(gòu)物時(shí)，會(huì)發(fā)生邊界層的分離現(xiàn)象，進(jìn)而對(duì)結(jié)構(gòu)產(chǎn)生周期性的氣動(dòng)力作用，這就是所謂的繞流問題.繞流問題一直是眾多學(xué)者研究的熱點(diǎn)，現(xiàn)有的研究成果大都以圓柱為研究對(duì)象.史里希廷［1]依據(jù)實(shí)驗(yàn)結(jié)果繪制了圓柱繞流的斯托羅哈數(shù)St、阻力系數(shù)Cd隨雷諾數(shù)（Re從1到107）的變化曲線，發(fā)現(xiàn)圓柱繞流的阻力系數(shù)值在Re＝5×105左右時(shí)突然由1.2下降到0.3，這種阻力值突然顯著減小的現(xiàn)象被稱為阻力危機(jī)現(xiàn)象.曹豐產(chǎn)、項(xiàng)海帆［2]用基于一般曲線坐標(biāo)系和交錯(cuò)網(wǎng)格的差分法求解原始變量二維不可壓縮粘性流體的N－S方程，計(jì)算了雷諾數(shù)Re從100到1×105范圍內(nèi)圓柱的非定常繞流，得到了渦脫頻率f、阻力系數(shù)Cd，不足之處是在雷諾數(shù)較高時(shí)其阻力系數(shù)與實(shí)驗(yàn)結(jié)果偏差較大.詹昊等［3]利用Fluent數(shù)值模擬的方法對(duì)圓柱在不同雷諾數(shù)范圍內(nèi)繞流的氣動(dòng)力參數(shù)進(jìn)行了系統(tǒng)的分析，研究了渦脫落規(guī)律，驗(yàn)證了圓柱繞流阻力危機(jī)的存在.

方形截面也是一種典型的鈍體截面形式，現(xiàn)有文獻(xiàn)一般是研究雷諾數(shù)較低時(shí)的情形，對(duì)于高雷諾數(shù)時(shí)的方柱繞流則甚少涉及.Davis and Moore［4]研究了不同縱橫比矩形截面柱體的阻力以及尾流，但研究?jī)H限于較低雷諾數(shù)下通過實(shí)驗(yàn)研究了不同截面縱橫比的矩形截面柱體在Re＝70～20 000范圍內(nèi)斯托羅哈數(shù)St與雷諾數(shù)Re的關(guān)系，但沒有涉及另一重要的參數(shù)Cd.Kim D H［6]等采用大渦模擬的方法研究了Re＝3 000時(shí)的方柱，比較了方柱位于固定的狹小通道中和位于無限區(qū)域中的情形，發(fā)現(xiàn)當(dāng)柱體位于狹小通道中時(shí)，由于固壁邊界的存在，阻力值有著明顯的增大.

繞流中周期性的氣動(dòng)力作用有可能引起結(jié)構(gòu)物的破壞，因此采用流體計(jì)算軟件CFX進(jìn)行數(shù)值模擬，系統(tǒng)研究方柱和圓柱的St、Cd隨Re（從層流區(qū)到超臨界區(qū)）的變化規(guī)律，比較方柱和圓柱的繞流參數(shù)，對(duì)于進(jìn)一步認(rèn)識(shí)繞流問題以及工程結(jié)構(gòu)的設(shè)計(jì)具有重要意義.

1 控制方程

在直角坐標(biāo)系下，對(duì)二維不可壓縮粘性流體，其運(yùn)動(dòng)規(guī)律用納維斯－斯托克斯方程（N－S方程）［7]描述：

式中，U、V為流體沿x、y向的分速度；P為流體承受的壓力；Fx、Fy分別為流體所受體積力的x、y向分量.

2 參數(shù)定義

Re、Cd、St是描述繞流問題的幾個(gè)常用的參數(shù)，它們的定義如下：

雷諾數(shù)：表征慣性力與粘性力的一個(gè)無量綱參數(shù)：

式中，ρ為流體密度；U為來流風(fēng)速；D為柱體特征長(zhǎng)度；μ、v分別為流體的動(dòng)力、運(yùn)動(dòng)粘性系數(shù).

阻力系數(shù)：阻力是流體對(duì)柱體的順流方向的作用力，將阻力無量綱化得到阻力系數(shù)Cd

式中，F(xiàn)d為流體作用于單位長(zhǎng)度柱體上沿順流方向的分力.

斯托羅哈數(shù)：描述旋渦脫落的一個(gè)重要的無量綱參數(shù)

式中，f為渦脫頻率.

3 數(shù)值模型

當(dāng)流體流過鈍體的時(shí)候，由于流體本身粘性的作用，使得物體表面上的流體速度為零，而離開物體表面一定距離的流體速度則不受粘性影響，此處的流動(dòng)可以按照無粘來處理.在物面和可以按無粘處理的流體之間的這一部分流體就是邊界層.邊界層通常是很薄的一層，但在這一薄層中各項(xiàng)流體參數(shù)卻發(fā)生著劇烈的變化，存在著較大的速度梯度.因此，在繞流問題中精確地模擬邊界層對(duì)于結(jié)果的重要性至關(guān)重要，這就要求邊界層網(wǎng)格要足夠密，如圖1所示.本文流體計(jì)算域大小取為42D×24D，方柱與圓柱的特征長(zhǎng)度D相等，沿展向的厚度取D／200.

圖1 方柱、圓柱計(jì)算域及網(wǎng)格示意圖

在通常的兩方程渦粘性湍流模型中，k－ε模型能夠較好地模擬遠(yuǎn)離壁面處已充分發(fā)展的湍流流動(dòng)，kω模型則更廣泛地運(yùn)用于各種壓力梯度下的邊界層問題.而綜合了兩種模型各自優(yōu)勢(shì)的SST湍流模型，在近壁面處保留了原始的k－ω模型，而在遠(yuǎn)離壁面處應(yīng)用了k－ε模型，這一特點(diǎn)使得采用SST湍流模型來模擬繞流問題結(jié)果會(huì)更為精確.CFX數(shù)值模擬計(jì)算時(shí)沿柱體厚度方向（展向）取為很薄的D／200，即用很薄的三維模型來模擬二維問題.Re≤500時(shí)采用層流模型，Re＞500時(shí)采用SST湍流模型（中等湍流強(qiáng)度5%；壁面取光滑壁面，無粗糙度）.計(jì)算域邊界條件設(shè)置如下：

1）左側(cè)：速度入口，U＝U0，V＝0；

2）右側(cè)：出口邊界條件，平均靜壓P為零；

3）沿展向兩邊界采用對(duì)稱邊界條件：各變量沿法向的分量為零；

4）其余邊界采用無滑移固壁邊界條件：U＝0，V＝0.

4 模型驗(yàn)證及參數(shù)分析

4.1 方柱繞流模型驗(yàn)證及參數(shù)分析

表1為方柱繞流模擬計(jì)算的各參數(shù)及相應(yīng)的文獻(xiàn)結(jié)果，其中Cd為瞬態(tài)計(jì)算時(shí)對(duì)應(yīng)的阻力系數(shù)的平均值.

表1 方柱繞流數(shù)值計(jì)算各雷諾數(shù)下對(duì)應(yīng)的參數(shù)結(jié)果及相關(guān)文獻(xiàn)結(jié)果

續(xù)表1 方柱繞流數(shù)值計(jì)算各雷諾數(shù)下對(duì)應(yīng)的參數(shù)結(jié)果及相關(guān)文獻(xiàn)結(jié)果

由表1可知方柱仿真計(jì)算結(jié)果與現(xiàn)有文獻(xiàn)結(jié)果吻合較好，最大誤差為8%.從圖2、3中也可以形象地看出本次計(jì)算結(jié)果與文獻(xiàn)結(jié)果的擬合程度，這驗(yàn)證了本文方柱繞流模擬計(jì)算方法的正確性.方柱的St、Cd隨Re的變化規(guī)律如圖2、3所示，其中曲線為數(shù)值結(jié)果的擬合曲線，數(shù)據(jù)點(diǎn)為已有的文獻(xiàn)數(shù)據(jù).

由圖2可知，方柱的St在Re＜50時(shí)為零，這是由于繞流還未產(chǎn)生旋渦的脫落；在50＜Re＜100時(shí)，St由零急劇增大，在Re＝200附近時(shí)達(dá)到最大；在200＜Re＜1 000時(shí)，St隨著Re的增大又迅速減小，而在Re＞1 000后波動(dòng)較小，穩(wěn)定在0.126附近.由圖3可知，對(duì)于方柱繞流，在Re＜100時(shí)Cd隨著Re數(shù)的增大急劇減小.Re＞1 000以后，方柱的平均阻力系數(shù)值變得比較平穩(wěn)，波動(dòng)較小，穩(wěn)定在2.1附近.

4.2 圓柱繞流模型驗(yàn)證及參數(shù)分析

由表2可知圓柱仿真計(jì)算結(jié)果與現(xiàn)有文獻(xiàn)結(jié)果吻合得也比較好，從圖4、5中也可以形象地看出本次計(jì)算結(jié)果與文獻(xiàn)結(jié)果的擬合程度，這驗(yàn)證了本文圓柱繞流模擬計(jì)算方法的正確性.

表2 圓柱繞流數(shù)值計(jì)算各雷諾數(shù)下對(duì)應(yīng)的參數(shù)結(jié)果及相關(guān)文獻(xiàn)結(jié)果

圓柱的St、Cd隨Re的變化規(guī)律如圖4、5所示，其中曲線為數(shù)值結(jié)果的擬合曲線，數(shù)據(jù)點(diǎn)為已有的文獻(xiàn)［1]的數(shù)據(jù).

圖4 圓柱的St－Re曲線

由圖4可知，圓柱的St在Re＜20時(shí)為零，因?yàn)榇藭r(shí)繞流還未產(chǎn)生脫落的旋渦；在20＜Re＜300時(shí)，St由零急劇增大到0.2；在300＜Re＜3×105時(shí)，St在0.2附近輕微波動(dòng)；而后，在Re＝1×106附近渦脫頻率很混亂，沒有穩(wěn)定的St值，可推知此時(shí)已經(jīng)不存在規(guī)則的渦街了；當(dāng)Re＞3.5×106后St值又出現(xiàn)并呈現(xiàn)逐漸增大的趨勢(shì)，可知此時(shí)規(guī)則的渦街又重新建立起來了.

圖5 圓柱的Cd－Re曲線

由圖5可知，對(duì)于圓柱繞流，在Re＜1 000時(shí)Cd隨著Re數(shù)的增大由12.5急劇減小到1.2左右；Re＞1 000以后，圓柱的平均阻力系數(shù)值變得比較平穩(wěn)，波動(dòng)較小，穩(wěn)定在1.2附近；在Re＝1×106附近，可以看到Cd值突然減小.此時(shí)邊界層正經(jīng)歷著從層流分離向湍流分離的過渡階段，旋渦脫落趨于不規(guī)則化，圓柱前方的正壓值不變，后方旋渦脫落區(qū)的負(fù)壓絕對(duì)值減小，導(dǎo)致流體對(duì)柱體的總壓減小，阻力系數(shù)也隨之減小.當(dāng)邊界層分離完全轉(zhuǎn)化為湍流分離后，阻力值又會(huì)增大并趨于穩(wěn)定.

4.3 方柱與圓柱繞流數(shù)值結(jié)果對(duì)比

為了分析方柱和圓柱繞流的不同，本文比較了方柱和圓柱數(shù)值模擬結(jié)果中的St、Cd隨Re的變化規(guī)律，如圖6～7所示.

圖6 方柱、圓柱的St－Re曲線

圖6為方柱和圓柱繞流St隨雷諾數(shù)Re變化的比較，可以分為3階段：第1階段，旋渦開始脫落后，St隨雷諾數(shù)Re的增大而迅速增加，方柱在Re＝200時(shí)達(dá)到極大值點(diǎn)，圓柱要在Re＝300左右才達(dá)到極大值點(diǎn)；第2階段方柱在200＜Re＜1 000范圍內(nèi)的St開始緩慢減小，而圓柱的St在300＜Re＜3×105范圍內(nèi)基本上穩(wěn)定在0.2；第3階段方柱的St隨雷諾數(shù)Re的變化趨于穩(wěn)定的值，而圓柱的St卻又開始增大.此外，第2、3階段方柱的St一直小于圓柱的St.

圖7為方柱和圓柱繞流阻力系數(shù)Cd隨雷諾數(shù)Re變化曲線的比較，可見其整體趨勢(shì)大體一致，但計(jì)算發(fā)現(xiàn)方柱不存在明顯的阻力危機(jī)現(xiàn)象，這與兩者的分離點(diǎn)位置的不同有很大的關(guān)系.圓柱繞流中流體邊界層分離點(diǎn)會(huì)隨著雷諾數(shù)的增大而逐漸后移，而方柱的分離點(diǎn)則固定位于其兩個(gè)前角點(diǎn)處.在Re＜100時(shí)兩者的Cd均隨著Re數(shù)的增大急劇減小.Re＞1 000以后，方柱與圓柱的阻力值基本上變化幅度不大，但方柱的阻力值總是大于相應(yīng)的圓柱的阻力值.

圖7 方柱、圓柱的Cd－Re曲線

5 結(jié) 論

運(yùn)用CFX可以較為精確地模擬方柱和圓柱的繞流問題，通過以上分析，可以得出以下結(jié)論：

1）方柱繞流中St值在較高雷諾數(shù)Re時(shí)是趨于穩(wěn)定的，而圓柱繞流卻不同，在穩(wěn)定變化階段過后，隨著Re的增大St值又開始增大了.并且，在Re＞200后方柱的St值比圓柱的St值小37%以上.

2）對(duì)于方柱和圓柱繞流，在Re＜100時(shí)Cd隨Re變化得很劇烈，而在Re＞1 000后，Cd隨Re的變化基本上趨于穩(wěn)定.在Re＞1 000后方柱的Cd值一直比圓柱的Cd值高60%以上，這是由方柱表面存在著4個(gè)尖銳的角點(diǎn)而圓柱表面比較光滑所致.

3）由于方柱的分離點(diǎn)始終位于其兩個(gè)前角點(diǎn)不變，使得在圓柱繞流中明顯存在的阻力危機(jī)現(xiàn)象，在方柱繞流中并不明顯，亦或要在更大的雷諾數(shù)下才會(huì)發(fā)生.

［1]史里希廷 H.邊界層理論（上冊(cè)）［M].7版.徐燕侯，譯.北京：科學(xué)出版社，1988.

［2]曹豐產(chǎn)，項(xiàng)海帆.圓柱非定常繞流及渦激振動(dòng)的數(shù)值計(jì)算［J].水動(dòng)力學(xué)研究與進(jìn)展（A 輯），2001，16（1）：111－118.

［3]詹 吳，李萬平.不同雷諾數(shù)下圓柱繞流仿真計(jì)算［J].武漢理工大學(xué)學(xué)報(bào)，2008，12（30）：129－132.

［4]Davis R W，Moore E F.A Numerical Study of Vortex Shedding from Rectangles［J].Fluid Mechanics，1982，116：475－506.

［5]Okajima A.Strouhal Numbers of Rectangular Cylinders［J].Fluid Mechanics，1982，123：379－398.

［6]Kim D H，Yang K S，Senda M.Large－eddy Simulation of Turbulent Flow Past a Square Cylinder Confined in a Channel［J].Computers ＆Fluids，2004，33：81－96.

［7]陳政清，項(xiàng)海帆.橋梁風(fēng)工程［M].北京：人民交通出版社，2005.

［8]Singh A P，De A K，Carpenter V K.Flow Past a Transversely Oscillating Square Cylinder in Free Stream at Low Reynolds Numbers［J].International Journal for Numerical Methods in Fluids，2009，61：658－682.

［9]Wan D C，Patnaik B S W，Wei G W.Discrete Singular Convolution－finite Subdomain Method for the Solution of Incompressible Viscous Flows［J].Computational Physics，2002，18：229－255.

［10]Norberg C.Flow Around Rectangular Cylinder：Pressure Forces and Wake Frequencies［J].Journal of Wind Engineering and Industrial Aerodynamics，1993，49：187－196.

［11]Lyn D A，Einav S，Rodi W.A Laser Doppler Velocimetry Study of Ensemble－averaged Characteristic of the Turbulent Flow Near Wake of a Square Cylinder［J].Fluid Mechanics，1995，304：285－319.

［12]Ansumali S，Chikatamarla S S，F(xiàn)rouzakis C M.Entropic Lattice Boltzmann Simulation of the Flow Past Square Cylinder［J].Modern Physics，2004，15（3）：435－445.