李 鰲，程文明，王書(shū)標(biāo)，杜 潤(rùn)

(1. 西南交通大學(xué)機(jī)械工程學(xué)院，四川 成都 610031；2..軌道交通運(yùn)維技術(shù)與裝備四川省重點(diǎn)實(shí)驗(yàn)室，四川 成都 610031)

1 引言

鈍體繞流作為風(fēng)工程領(lǐng)域基礎(chǔ)問(wèn)題，一直深受?chē)?guó)內(nèi)外學(xué)者廣泛關(guān)注。在風(fēng)洞試驗(yàn)和數(shù)值仿真方面有大量的研究資料。Lyn，Einav等[1]應(yīng)用LDV(激光多普勒風(fēng)速測(cè)量技術(shù))對(duì)方柱繞流進(jìn)行風(fēng)洞試驗(yàn)。Norberg[2]在低速閉環(huán)風(fēng)洞對(duì)不同寬高比矩形柱繞流進(jìn)行了系統(tǒng)研究。秦浩，肖姚等[3]結(jié)合PIV(粒子圖像測(cè)速法)和POD(正交分解法)對(duì)方柱繞流尾流場(chǎng)進(jìn)行了解析和研究，結(jié)果與Lyn，Einav等[1]的吻合較好，同時(shí)應(yīng)用γ-Reθ和SST k-ω?cái)?shù)值模擬也得到了不錯(cuò)的結(jié)果。數(shù)值研究方面，國(guó)內(nèi)外學(xué)者就方柱繞流[4-8]運(yùn)用LES(大渦模擬)仿真，均獲得不錯(cuò)效果。Sohankar[9]采用LES對(duì)雷諾數(shù)為Re=105的寬高比為0.4～4的矩形進(jìn)行仿真分析，得到的結(jié)果和Norberg[2]實(shí)驗(yàn)結(jié)果相一致。Tian，Ong等[10]采用SST k-ω對(duì)寬高比0.05-1的矩形進(jìn)行分析，但對(duì)寬高比小于0.6時(shí)結(jié)果不理想。Islam，Zhou等[11]應(yīng)用ILBM(不可壓縮亞格子玻爾茲曼模型)對(duì)低雷諾數(shù)(Re≤250)矩形進(jìn)行分析，得到與高雷諾數(shù)不同的結(jié)論。Mannini，Soda等[12]對(duì)寬高比5.0的矩形討論了風(fēng)向角和雷諾數(shù)的影響，并與風(fēng)洞結(jié)果吻合較好。

實(shí)際生活中，矩形柱和類(lèi)矩形柱的結(jié)構(gòu)很多，例如起重機(jī)箱型梁結(jié)構(gòu)，針對(duì)其展開(kāi)研究意義明顯，而試驗(yàn)和仿真也均表明：矩形柱繞流在阻力系數(shù)上存在一個(gè)臨界值，即寬高比0.62附近時(shí)阻力系數(shù)最大。而對(duì)于其機(jī)理產(chǎn)生，也作了各種解釋和分析[13， 14]。但是單獨(dú)對(duì)其進(jìn)行分析和研究的討論略顯不夠。

因此，本文主要針對(duì)該臨界情況進(jìn)行深入分析。選擇Norberg[2]文中雷諾數(shù)Re=19000作為研究條件，采用大渦模擬方法作為主要手段，通過(guò)矩形柱氣動(dòng)特征，流場(chǎng)及尾渦特性等等進(jìn)行分析和討論，研究阻力和結(jié)構(gòu)關(guān)系和機(jī)理，并進(jìn)行詳細(xì)地分析與討論。

2 湍流模型與計(jì)算方法

2.1 大渦模擬模型

對(duì)不可壓縮流體控制方程N(yùn)-S方程而言，經(jīng)過(guò)濾波得到LES的連續(xù)性和動(dòng)量方程可以表示為

(1)

(2)

(3)

其中，νSGS為亞格子模型渦粘系數(shù)，表達(dá)式為

(4)

(5)

2.2 計(jì)算流場(chǎng)及算法

針對(duì)R=0.62矩形柱的數(shù)值計(jì)算仿真的流場(chǎng)域及相應(yīng)的邊界條件如圖1所示(針對(duì)二維模擬展向L為0)。

計(jì)算域在空間各方向無(wú)量綱長(zhǎng)度(即長(zhǎng)度與特征高度的比值)Lu，LD，h，H，L分別為12，30，1，20，4。本文展向長(zhǎng)度方向選取的是L/b=6.45(L/h=4)。針對(duì)展向長(zhǎng)度L/b的選取，國(guó)內(nèi)外學(xué)者對(duì)此也進(jìn)行了很多相關(guān)性實(shí)驗(yàn)和仿真研究。Vickery[15]，McLean & Gartshore[16]等通過(guò)實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證了方柱展向長(zhǎng)度相關(guān)性長(zhǎng)度大約在L/b=5.0，Tamura，Miyaji等[17]通過(guò)仿真計(jì)算也得到展向長(zhǎng)度在L/b=4.0～10.0時(shí)結(jié)果差異不大的結(jié)論，周強(qiáng)，廖海黎等[7]也針對(duì)方柱展向長(zhǎng)度進(jìn)行分析和討論。

計(jì)算邊界條件的選取如下：入口處采用inlet速度入口邊界，均勻來(lái)流速度條件：u=U， v=w=0；出口采用零壓力出口：擬壓為0；柱體壁面采用無(wú)滑移壁面：No-slip；上下端面采用對(duì)稱(chēng)邊界條件：Symmetry；展向邊界采用周期邊界條件：Periodic。

圖1 計(jì)算流場(chǎng)域和邊界條件

在空間離散上，采用有界的二階中心差分格式，時(shí)間離散上，采用有界的二階隱式法，對(duì)壓力速度耦合采用PISO求解算法；收斂精度設(shè)置為10-5，無(wú)量綱時(shí)間步長(zhǎng)為ΔtU / D=2.78×10-3，以確保庫(kù)朗數(shù)小于1，且迭代在20次內(nèi)完成。

流場(chǎng)域選用六面體網(wǎng)格進(jìn)行空間離散。三維流域在XY平面采用混合網(wǎng)格，即柱體表面 5h 范圍內(nèi)采用O型網(wǎng)格，外圍采用正交型網(wǎng)格，并沿著展向z方向進(jìn)行延伸。貼體O型網(wǎng)格采用較低網(wǎng)格增長(zhǎng)率(1.05)以考慮邊界層效應(yīng)；外圍則采用較大網(wǎng)格增長(zhǎng)率(1.1)在保證計(jì)算合理基礎(chǔ)上減少計(jì)算成本。網(wǎng)格在XY方向示意圖如圖2所示。二維流場(chǎng)及網(wǎng)格設(shè)置。

圖2 計(jì)算流場(chǎng)域內(nèi)XY平面內(nèi)網(wǎng)格劃分示意圖

3 網(wǎng)格無(wú)關(guān)性驗(yàn)證

在進(jìn)行繞流分析前一般需要進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證，這里選用3套疏密不同的網(wǎng)格，網(wǎng)格差異主要在于表面首層網(wǎng)格高度的不同，Case1、Case2 和Case3沿著展向方向解析度為Δz/h=0.1，首層網(wǎng)格高度設(shè)置有所差異，具體參數(shù)如表1所示。將三組不同網(wǎng)格計(jì)算得到的平均積分分量結(jié)果與其余相關(guān)實(shí)驗(yàn)和數(shù)值結(jié)果進(jìn)行了對(duì)比，結(jié)果如表1所示。對(duì)比三組網(wǎng)格Case1、Case2 和Case3，結(jié)果間差異不大，且與參考值間吻合較好。綜合考慮計(jì)算精度和成本，選用Case2結(jié)果作為參照。此外，表1中還列出了采用二維LES計(jì)算得到的結(jié)果，即Case4數(shù)據(jù)作為對(duì)比。

表1 計(jì)算結(jié)果對(duì)比

4 結(jié)果與討論

以高雷諾數(shù)Re=19000時(shí)，臨界寬高比R=0.62下的矩形柱為研究對(duì)象，結(jié)合通用的LES湍流模型對(duì)矩形柱結(jié)構(gòu)下的氣動(dòng)特征，繞流特性及其產(chǎn)生機(jī)理等進(jìn)行研究和分析，同時(shí)對(duì)二維和三維LES進(jìn)行對(duì)比研究。

4.1 氣動(dòng)特性研究

表1對(duì)比了計(jì)算所得三分力系數(shù)及頻率參數(shù)等與前人試驗(yàn)和仿真所得結(jié)果以及方柱繞流對(duì)比情況，所得三分力系數(shù)可直接通過(guò)計(jì)算所得時(shí)程曲線(xiàn)(如圖3所示LES計(jì)算時(shí)程曲線(xiàn))取平均值得到，頻率參數(shù)St 數(shù)需要對(duì)升力系數(shù)結(jié)果進(jìn)行FFT變換得到渦脫頻率 f，再根據(jù)計(jì)算公式St=fh/U 得到(h 為特征高度，U為來(lái)流速度)。

圖3 三維LES升力系數(shù) Cd 和阻力系數(shù) Cl 時(shí)程曲線(xiàn)

由表1結(jié)果可知，三維LES計(jì)算得到Case2的平均阻力系數(shù) Cd與參考仿真值[9， 10， 18]間的差異在1%，與實(shí)驗(yàn)值[2， 13]吻合也很好； St 數(shù)與實(shí)驗(yàn)值[2， 13]均比較吻合，相對(duì)誤差也在2%左右；而升力系數(shù)均方根 Clrms與LES仿真結(jié)果[9]吻合，遠(yuǎn)小于二維RANS仿真[10， 18]結(jié)果。而Case4為二維LES計(jì)算所得，計(jì)算所得Cd和St 數(shù)與三維LES結(jié)果相似，均和參考值間吻合很好；而升力系數(shù)均方根 Clrms高于Sohankar[9]采用三維LES所得結(jié)果(相對(duì)偏差28.4%)，同時(shí)又小于RANS計(jì)算[10， 18]所得結(jié)果(相對(duì)偏差分別為51.5%和39.1%)。

表1對(duì)比結(jié)果說(shuō)明氣動(dòng)力參數(shù)的預(yù)測(cè)在阻力系數(shù)和頻率上，采用各模型計(jì)算所得結(jié)果都比較好；而在升力系數(shù)均方根上，采用三維LES結(jié)果最好，其次是二維LES，采用二維RANS計(jì)算預(yù)測(cè)效果最差。

圖4為柱體截面表面平均壓力系數(shù) Cp的分布情況，同時(shí)與方柱表面壓力分布的實(shí)驗(yàn)[19]及仿真結(jié)果[7， 8]進(jìn)行了對(duì)比。由圖可知，柱體前端面AB面和上端面BC面，矩形柱和方柱壓力分布差異不大，差異主要發(fā)生在后端面CD面：矩形柱后端面壓力分布明顯小于方柱，且在D處差值最大。而后端面壓力分布反映的是繞流阻力的大小，也表明了該矩形柱阻力系數(shù)更大，與表1中Cd結(jié)果相呼應(yīng)。

同時(shí)，BC段上端面壓力系數(shù)分布可知，三維LES結(jié)果大于二維LES，即二維LES上端面負(fù)壓更大，而上下端面壓力分布反映的是繞流升力系數(shù)相關(guān)參數(shù)，表明二維仿真二維LES升力參數(shù)估值過(guò)大，與表1中Clrms結(jié)果相呼應(yīng)。也說(shuō)明了Clrms是一個(gè)與三維相關(guān)的參數(shù)，與周強(qiáng)等[7]方柱繞流結(jié)果相吻合。

圖4 矩形柱中截面表面壓力系數(shù)Cp分布

4.2 流場(chǎng)特性研究

圖5給出并對(duì)比了中心平面(y=0平面)上平均流向速度的分布情況(橫坐標(biāo)X-position=x/h，下同)。有圖可以看出，矩形柱和方柱一樣，流向速度在柱體后壁面速度為零，然后逐漸減小到最小速度Umin，再單調(diào)增加到與來(lái)流速度 U 相同為止。相對(duì)于方柱，矩形柱最小速度比方柱更小，最小值到結(jié)構(gòu)中心截面(x=0平面)距離更短(定義為繞流回轉(zhuǎn)長(zhǎng)度)，說(shuō)明矩形柱后的尾渦距離更短，即矩形柱背壓更大(如圖4)，因此阻力系數(shù)更大(如表1)。同時(shí)，尾流區(qū)內(nèi)，矩形柱流向速度在 x/h > 6時(shí)約為0.4U小于方柱的0.4U，說(shuō)明矩形柱尾流區(qū)內(nèi)尾渦區(qū)域更長(zhǎng)。

對(duì)比二維LES和三維LES結(jié)果發(fā)現(xiàn)，二維LES計(jì)算結(jié)果所得的回轉(zhuǎn)長(zhǎng)度略小，說(shuō)明二維LES計(jì)算阻力系數(shù)更大，與表1計(jì)算所得氣動(dòng)參數(shù)相對(duì)應(yīng)；最小流速 umin/U 相較于三維LES偏大，說(shuō)明展向長(zhǎng)度對(duì)umin/U有較大影響，對(duì)回轉(zhuǎn)長(zhǎng)度影響較小，所得結(jié)果存在一定偏差。

圖5 矩形柱中心線(xiàn)平均流向速度分布

圖6和圖7分別為尾流區(qū)內(nèi)距離柱體后壁面分別為x/h=1，2和3.5時(shí)流向速度，豎向速度的平均值及均方根，以及對(duì)應(yīng)的雷諾應(yīng)力。通過(guò)對(duì)比不同位置處速度及其脈動(dòng)量的差異可以反映尾流區(qū)內(nèi)尾渦變化情況，同時(shí)對(duì)比方柱和矩形柱在相同位置處速度及其脈動(dòng)值差異可以反映不同結(jié)構(gòu)對(duì)尾流區(qū)的影響，也可以對(duì)比二維和三維仿真在尾流區(qū)的差異情況。

圖6(a)為尾流區(qū)平均流向速度不同位置速度分布，速度分布以y=0處為中心呈現(xiàn)“U型”的對(duì)稱(chēng)分布，且隨著遠(yuǎn)離尾流區(qū)逐漸變緩；圖6(b)為平均豎向速度分布，分布以y=0處為中心呈現(xiàn)“N型”的反對(duì)稱(chēng)分布。與方柱試驗(yàn)及仿真結(jié)果對(duì)比發(fā)現(xiàn)，矩形柱在y=0附近速度相較于方柱更小，整體上“U型”更窄，“N型”更緩，這些說(shuō)明了尾渦區(qū)內(nèi)矩形柱的尾渦寬度更窄，間距更短，也即印證了矩形柱阻力系數(shù)更大，和表1結(jié)果相呼應(yīng)。同時(shí)，對(duì)比二維LES和三維LES結(jié)果二者平均速度差異不大，結(jié)果基本一致，說(shuō)明展向方向的長(zhǎng)度對(duì)平均速度流場(chǎng)影響不大。

圖6 尾流區(qū)不同位置平均速度分布(—，Present-Rectangular-3DLES；---，Present-Rectangular-3DLES；，Exp-Square(Lyn et，al，Re21400)；--×--，LES-Square(Zhou et，al，Re22000))

速度平均值反映的是尾渦中結(jié)構(gòu)的時(shí)均特性，而脈動(dòng)值和雷諾應(yīng)力則可以反映結(jié)構(gòu)的湍流特性，如圖7為尾流區(qū)內(nèi)的速度均方根及雷諾應(yīng)力分布。圖7(a)為來(lái)流速度均方根分布，由于卡門(mén)渦街的存在，存在兩個(gè)峰值，呈“M型”對(duì)稱(chēng)分布于中心點(diǎn)兩側(cè)，同時(shí)隨著遠(yuǎn)離后壁面而變緩并消失，說(shuō)明旋渦在下游逐漸減弱；圖7(b)為豎向速度均方根分布，在中心線(xiàn)處呈現(xiàn)“∩型”對(duì)稱(chēng)分布；圖7(c)為雷諾應(yīng)力分布，呈“N型”反對(duì)稱(chēng)分布。

對(duì)比圖7中方柱和矩形柱脈動(dòng)參數(shù)可以發(fā)現(xiàn)，矩形柱的流向速度均方根分布在x/h=3.5時(shí)“M型”雙峰依舊明顯，說(shuō)明旋渦在流向速度方向相較于方柱衰減更慢；豎向速度方向“∩型”相較于方柱在靠近后壁面時(shí)峰值稍大，而遠(yuǎn)離柱體后趨于一致；雷諾應(yīng)力上矩形柱和方柱在各個(gè)位置上表現(xiàn)差異不明顯。以上說(shuō)明方柱和矩形柱在湍流特性上的差異主要體現(xiàn)在流向方向上，矩形柱相較于方柱尾渦更大，且更接近于后壁面，同時(shí)尾流區(qū)更長(zhǎng)。

圖7 尾流區(qū)不同位置速度均方根分布

同時(shí)，圖7中二維LES和三維LES結(jié)果對(duì)比發(fā)現(xiàn)，二維LES的仿真結(jié)果變化幅值更大，趨勢(shì)更明顯，說(shuō)明二維LES仿真在尾渦模擬上會(huì)進(jìn)行一定程度上的放大，也說(shuō)明展向長(zhǎng)度對(duì)速度脈動(dòng)值模擬，即尾渦湍流特征的描述上有一定影響。這也解釋了表1中二維LES得到升力系數(shù)均方根偏大的原因。

通過(guò)對(duì)寬高比0.62的矩形柱的三分力系數(shù)，表面壓力系數(shù)分布，流向速度分布，尾流區(qū)不同位置時(shí)均速度，湍流特性等等結(jié)果的與方柱間的對(duì)比分析和討論，解釋說(shuō)明了該臨界寬高比下矩形柱的特殊性和相應(yīng)的產(chǎn)生機(jī)理。而二維和三維計(jì)算結(jié)果的對(duì)比也說(shuō)明進(jìn)行三維計(jì)算的必要性。

5 結(jié)論

針對(duì)高雷諾數(shù)Re=19000下，臨界寬高比為0.62的矩形柱，采用大渦模擬(LES)方法，就二維和三維模型進(jìn)行繞流數(shù)值模擬。在驗(yàn)證網(wǎng)格及數(shù)值結(jié)果準(zhǔn)確性的基礎(chǔ)上，詳細(xì)分析了矩形柱的氣動(dòng)力和流場(chǎng)特性，并與方柱結(jié)果進(jìn)行分析對(duì)比，得到以下結(jié)論：

1)采用三維大渦模擬算法，對(duì)寬高比為0.62的矩形柱進(jìn)行繞流分析，得到氣動(dòng)參數(shù)等和相關(guān)實(shí)驗(yàn)和仿真數(shù)據(jù)吻合很好。同時(shí)對(duì)比矩形柱表面壓力系數(shù)，中心線(xiàn)來(lái)流速度分布，以及不同位置處平均速度以及湍流特性的分布，可以很好地解釋臨界寬高比0.62氣動(dòng)參數(shù)的特殊性。

2)對(duì)比二維和三維大渦模擬計(jì)算結(jié)果，可以發(fā)現(xiàn)：二維大渦模擬計(jì)結(jié)果在升力系數(shù)均方根，柱體上下壁面壓力系數(shù)分布，尾流區(qū)最小速度，速度脈動(dòng)值分布等的預(yù)測(cè)上偏差較大，說(shuō)明這些參數(shù)對(duì)繞流結(jié)構(gòu)的展向長(zhǎng)度較敏感，三維效應(yīng)明顯。

3)考慮到計(jì)算效率和計(jì)算成本，如果在只關(guān)心結(jié)構(gòu)阻力系數(shù)以及渦脫頻率時(shí)，選用二維計(jì)算更合理。

至于采用三維計(jì)算時(shí)，展向長(zhǎng)度對(duì)矩形柱繞流特性的影響程度和規(guī)律，還需要后續(xù)深入研究探討。