張志超，楊漢彬

(西南交通大學(xué)土木工程學(xué)院，四川成都 610036)

鈍體繞流(尤其是方形柱體繞流)是流體力學(xué)的經(jīng)典問(wèn)題之一。當(dāng)流體以一定速度流經(jīng)柱體結(jié)構(gòu)時(shí)，流體會(huì)在結(jié)構(gòu)后方產(chǎn)生規(guī)則的旋渦脫落，即卡門(mén)渦街現(xiàn)象。方形截面是典型的柱體截面形式之一，但工程中也存在方形帶倒角的截面形式，如各種大型橋梁的橋墩、高樁碼頭、高層建筑的立柱等。到目前為止，許多學(xué)者對(duì)方柱繞流問(wèn)題進(jìn)行了研究，Lyn等[1]對(duì)雷諾數(shù)Re=21400的方柱繞流流場(chǎng)特性進(jìn)行了詳細(xì)的研究，給出了流動(dòng)速度及其脈動(dòng)量的分布，為數(shù)值模擬提供了事實(shí)依據(jù)和考核標(biāo)準(zhǔn)。秦浩等[2]對(duì)高雷諾數(shù)下方柱繞流問(wèn)題進(jìn)行了PIV(Particle Image Velocimetry)試驗(yàn)和二維數(shù)值模擬，研究表明通過(guò)PIV試驗(yàn)可準(zhǔn)確觀測(cè)到方柱尾流的整個(gè)脈動(dòng)流場(chǎng)且γ-Reθ模型對(duì)方柱繞流流場(chǎng)的模擬十分準(zhǔn)確。付英男等[3]對(duì)方柱繞流進(jìn)行了三維數(shù)值模擬，發(fā)現(xiàn)Re=100、300時(shí)方柱繞流存在三維特性，導(dǎo)致三維數(shù)值模擬計(jì)算的阻力系數(shù)和升力系數(shù)均小于二維計(jì)算結(jié)果。李雪健等[4]對(duì)方柱繞流進(jìn)行了二維和三維數(shù)值模擬，結(jié)果發(fā)現(xiàn)三維數(shù)值模擬優(yōu)于二維數(shù)值模擬且隨著雷諾數(shù)的增加，Cd的平均值增加而St先增大后減小。王遠(yuǎn)成等[5]采用RNGk-ε湍流模型對(duì)方柱繞流流場(chǎng)進(jìn)行了數(shù)值模擬，結(jié)果表明RNGk-ε湍流模型可以成功模擬繞方柱的不穩(wěn)定，非定常和劇烈分離流動(dòng)。朱夢(mèng)楠等[6]對(duì)單個(gè)方柱及間距為1H，2H，3H的兩并列方柱的流動(dòng)進(jìn)行了數(shù)值模擬，研究發(fā)現(xiàn)當(dāng)兩并列方柱之間的間距為2H時(shí)，方柱壁面處的升阻力會(huì)發(fā)生突變，間距為3H時(shí)，兩并列方柱繞流數(shù)值模擬的結(jié)果與單方柱繞流的結(jié)果類(lèi)似。周強(qiáng)等[7]對(duì)高雷諾數(shù)下方柱繞流進(jìn)行了三維數(shù)值模擬，結(jié)果表明：雷諾數(shù)為22 000下方柱尾流區(qū)回轉(zhuǎn)長(zhǎng)度為1.37倍方柱寬帶，St為0.121，脈動(dòng)升力系數(shù)為1.40，展向長(zhǎng)度取8倍方柱寬帶可更準(zhǔn)確地獲得三維湍流特性；沈立龍等[8]等人對(duì)亞臨界區(qū)圓柱和方柱繞流進(jìn)行了數(shù)值模擬，發(fā)現(xiàn)同一雷諾數(shù)下圓柱繞流阻力系數(shù)較方柱繞流阻力系數(shù)要小，但圓柱的Strouhal數(shù)較方柱的要高，單圓柱和單方柱繞流生成的旋渦脫落周期不同，方柱的偏大，二者的旋渦脫落特性與流場(chǎng)運(yùn)動(dòng)具有相似特性。圓柱的分離點(diǎn)隨著雷諾數(shù)的增大而逐漸向柱后方推移，方柱的分離點(diǎn)則固定在柱前兩個(gè)棱角位置。雖然之前對(duì)方柱繞流問(wèn)題進(jìn)行了很多研究，但對(duì)于帶倒角方柱的研究卻很少，杜明倩等[9]采用數(shù)值模擬的方法對(duì)不同倒角半徑下方柱繞流進(jìn)行了二維數(shù)值模擬，研究發(fā)現(xiàn)增加倒角半徑會(huì)改變方柱的分離點(diǎn)，使尾流區(qū)長(zhǎng)度增加，旋渦尺度減小，柱體截面越接近圓形，斯托羅哈數(shù)會(huì)逐漸增大。Kumar等[10]使用粒子圖像流速測(cè)量技術(shù)(PIV)，研究了不同倒角半徑的流場(chǎng)特性，發(fā)現(xiàn)當(dāng)?shù)菇前霃皆酱髸r(shí)，尾流附近旋渦脫落越均勻。

本文運(yùn)用CFD軟件、基于k-ωSST模型對(duì)亞臨界區(qū)Re=22000的單方柱和倒角半徑為0.15D、0.25D和0.35D的方柱進(jìn)行了三維數(shù)值模擬，計(jì)算得到這幾種情況下的升阻力系數(shù)、斯特羅哈數(shù)等重要?jiǎng)恿W(xué)參數(shù)，并對(duì)圓柱、方柱和帶倒角方柱繞流機(jī)制進(jìn)行了分析。

1 控制方程與建模

1.1 控制方程

不可壓縮粘性流體的運(yùn)動(dòng)可用Navier-Stokes方程來(lái)描述，連續(xù)性方程與動(dòng)量方程分別為：

式中：ui，uj為速度分量；p為壓力；ρ為流體的密度；v為流體的動(dòng)力粘性系數(shù)。

對(duì)于湍流情況，本文采用k-ωSST模型，該模型綜合了k-ω模型在近壁區(qū)計(jì)算的優(yōu)點(diǎn)和k-ε模型在遠(yuǎn)場(chǎng)計(jì)算的優(yōu)點(diǎn)，將k-ω模型和標(biāo)準(zhǔn)k-ε模型都乘以一個(gè)混合函數(shù)后再相加就得到這個(gè)模型。在近壁區(qū)，混合函數(shù)的值等于1，因此在近壁區(qū)等價(jià)于k-ω模型。在遠(yuǎn)離壁面的區(qū)域混合函數(shù)的值則等于0，因此自動(dòng)轉(zhuǎn)換為標(biāo)準(zhǔn)k-ε模型。與標(biāo)準(zhǔn)k-ω模型相比，k-ωSST模型中增加了橫向耗散導(dǎo)數(shù)項(xiàng)，同時(shí)在湍流粘度定義中考慮了湍流剪切應(yīng)力的輸運(yùn)過(guò)程，模型中使用的湍流常數(shù)也有所不同。這些特點(diǎn)使得k-ωSST模型的適用范圍更廣，如可以用于帶逆壓梯度的流動(dòng)計(jì)算、翼型計(jì)算、跨音速激波計(jì)算等。模型方程如下：

(3)

(4)

式中：k為湍流動(dòng)能，ω為湍流比耗散率，σω2=0.856，β*=0.09，混合函數(shù)`定義如下：

(5)

(6)

式中：y表示距離壁面的距離，CDkω代表比耗散率輸運(yùn)方程中交錯(cuò)擴(kuò)散項(xiàng)的正值部分，渦粘系數(shù)定義為

(7)

式中：a1=0.31，Ω為渦量，混合函數(shù)F2定義如下：

(8)

模型中常數(shù)的取值為OpenFOAM中的默認(rèn)值。

1.2 參數(shù)定義

1.3 邊界條件和初始條件設(shè)置

計(jì)算流場(chǎng)的邊界條件和初始條件設(shè)置如下：

初始條件：初始設(shè)定一個(gè)從左向右的速度場(chǎng)，u=U0=0.22m/s，v=0，w=0；其中，u、v和w分別為x、y和z方向的速度，方柱邊長(zhǎng)D=0.1m。

邊界條件：來(lái)流方向?yàn)榫鶆蛄魉賃0，出口為自由出流，方柱采用無(wú)滑移固壁邊界，其余邊界采用對(duì)稱(chēng)邊界(圖1)。

(a)XY方向計(jì)算域示意

(b)Z方向計(jì)算域示意

網(wǎng)格劃分采用六面體結(jié)構(gòu)化網(wǎng)格，方柱壁面處y+取20，具體網(wǎng)格劃分如圖2所示。

(a)方柱

(b)倒角0.15D

(c)倒角0.25D

(d)倒角0.35D

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

2.1 升阻力系數(shù)

圖3為利用OpenFOAM軟件模擬Re為22 000時(shí)方柱和倒角半徑為0.15D，0.25D和0.35D方柱阻力系數(shù)Cd及升力系數(shù)Cl隨時(shí)間t變化。

由圖3(a)可以看出，單方柱的阻力系數(shù)和升力系數(shù)在20s附近開(kāi)始穩(wěn)定，阻力系數(shù)穩(wěn)定在2.12，升力系數(shù)在-1.92到1.92之間波動(dòng)。當(dāng)?shù)菇前霃皆龃蟮?.15D時(shí)，方柱的升力系數(shù)和阻力系數(shù)都開(kāi)始減小且阻力系數(shù)降低約46 %，而后當(dāng)?shù)菇前霃嚼^續(xù)增加時(shí)，方柱的阻力系數(shù)繼續(xù)減小，升力系數(shù)變化不大，當(dāng)?shù)菇前霃綖?.35D時(shí)，阻力系數(shù)相對(duì)于倒角半徑為0.25D時(shí)無(wú)明顯變化，但升力系數(shù)卻開(kāi)始增大。

表1為三維計(jì)算結(jié)果與實(shí)驗(yàn)及以往文獻(xiàn)結(jié)果的對(duì)比，結(jié)果顯示，Re=22000時(shí)的模擬結(jié)果能很好的與參考文獻(xiàn)中數(shù)據(jù)相吻合，即說(shuō)明本文模擬結(jié)果的正確性。

2.2 斯托羅哈數(shù)

斯托羅哈數(shù)又被稱(chēng)為渦脫落的無(wú)量綱頻率。根據(jù)模擬所得結(jié)果繪制了Re=22000時(shí)單方柱和倒角半徑為0.15D、0.25D和0.35D方柱繞流的St隨倒角半徑的變化曲線(xiàn)如圖4所示。

(a)方柱

(b)倒角半徑0.15D

(c)倒角半徑0.25D

(d)倒角半徑0.35D

由上圖可知，單方柱的St為0.124，當(dāng)?shù)菇前霃皆龃蟮?.15D時(shí)，St以較快速度增大到0.212，當(dāng)?shù)菇前霃嚼^續(xù)增大時(shí)，St的增速變得平緩，說(shuō)明方柱增加倒角后，旋渦脫落頻率

表1 Re=22 000時(shí)單方柱繞流模擬計(jì)算結(jié)果

圖4 方柱繞流的斯托羅哈數(shù)隨倒角半徑變化曲線(xiàn)

會(huì)加快，脈動(dòng)更強(qiáng)，這也與文獻(xiàn)[9]中的結(jié)果相一致。

2.3 三維渦量圖

方柱及倒角半徑為0.15D、0.25D和0.35D方柱的三維渦量圖如圖5所示。

(a)方柱

(b)倒角半徑0.15D

(c)倒角半徑0.25D

(d)倒角半徑0.35D

從上圖可以看出，由于方柱具有固定的分離角，其旋渦脫落在高度方向一致，而帶倒角方柱具有明顯的胞體結(jié)構(gòu)，呈現(xiàn)明顯的三維特征。由于胞體結(jié)構(gòu)的存在說(shuō)明其旋渦脫落的相關(guān)性較弱，這也是導(dǎo)致其阻力系數(shù)小于方柱的原因之一。

2.4 邊界層分離

亞臨界區(qū)雷諾數(shù)下方柱繞流邊界層分離為層流分離，且邊界層內(nèi)液體運(yùn)動(dòng)時(shí)會(huì)伴有能量損失，本次模擬的方柱雷諾數(shù)Re=22000，正好處于亞臨界區(qū)。對(duì)于圓柱見(jiàn)圖6(a)，隨著雷諾數(shù)的增大，分離點(diǎn)會(huì)逐漸向后推移。在壓強(qiáng)最大的P點(diǎn)至柱上下A、B兩點(diǎn)的邊界層內(nèi)液流處于加速減壓狀態(tài)，在A、B兩點(diǎn)時(shí)壓強(qiáng)最小，流速最大，在這兩點(diǎn)之后液流處于減速增壓狀態(tài)，當(dāng)動(dòng)能減小至零無(wú)法提供壓能時(shí)，完成分離。對(duì)于方柱見(jiàn)圖6(b)分離點(diǎn)固定在柱前的兩個(gè)棱角位置。而對(duì)于帶倒角方柱見(jiàn)圖6(c)分離點(diǎn)則相對(duì)復(fù)雜，其分離點(diǎn)介于圓柱和方柱之間。

(a)圓柱

(b)方柱

(c)帶倒角方柱

3 結(jié)論

本文運(yùn)用OpenFOAM軟件，通過(guò)對(duì)亞臨界區(qū)Re=22000下單方柱和倒角半徑為0.15D、0.25D和0.35D的方柱進(jìn)行了數(shù)值模擬，得出以下結(jié)論：

(1)Re=22000時(shí)，單方柱的升力系數(shù)和阻力系數(shù)最大，隨著倒角半徑的增加，升阻力系數(shù)會(huì)減小，但當(dāng)?shù)菇前霃綖?.35D時(shí)，升力系數(shù)相對(duì)于倒角半徑半徑為0.15D和0.25D時(shí)會(huì)增加，但始終小于方柱的升力系數(shù)。

(2)帶倒角方柱的三維數(shù)值模擬會(huì)出現(xiàn)胞體結(jié)構(gòu)，呈現(xiàn)明顯的三維特性，柱體高度方向旋渦脫落的相關(guān)性較弱。

(3)在一定范圍內(nèi)，帶倒角方柱的斯托羅哈數(shù)會(huì)隨著倒角半徑的增大而逐漸增大。

(4)單方柱繞流存在固定的分離點(diǎn)，即柱前兩個(gè)棱角的位置，當(dāng)方柱出現(xiàn)倒角后，方柱的分離點(diǎn)發(fā)生改變，不再固定在之前的棱角位置，而是介于圓柱和方柱之間。隨著倒角半徑的增大，帶倒角方柱的流動(dòng)特性越來(lái)越近似于圓柱。