寬高比為0.62的矩形柱繞流特性研究

李 鰲，程文明，王書標，杜 潤

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

1 引言

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

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

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

2 湍流模型與計算方法

2.1 大渦模擬模型

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

(1)

(2)

(3)

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

(4)

(5)

2.2 計算流場及算法

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

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

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

圖1 計算流場域和邊界條件

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

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

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

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

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

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

4 結(jié)果與討論

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

4.1 氣動特性研究

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

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

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

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

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

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

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

4.2 流場特性研究

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

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

圖5 矩形柱中心線平均流向速度分布

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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