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

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

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

1 引言

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

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

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

2 湍流模型與計算方法

2.1 大渦模擬模型

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

(1)

(2)

(3)

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

(4)

(5)

2.2 計算流場及算法

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

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

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

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

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

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

圖2 計算流場域內XY平面內網格劃分示意圖

3 網格無關性驗證

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

表1 計算結果對比

4 結果與討論

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

4.1 氣動特性研究

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

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

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

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

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

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

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

4.2 流場特性研究

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

對比二維LES和三維LES結果發現，二維LES計算結果所得的回轉長度略小，說明二維LES計算阻力系數更大，與表1計算所得氣動參數相對應；最小流速 umin/U 相較于三維LES偏大，說明展向長度對umin/U有較大影響，對回轉長度影響較小，所得結果存在一定偏差。

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

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

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

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

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

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

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

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

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

5 結論

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

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

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

3)考慮到計算效率和計算成本，如果在只關心結構阻力系數以及渦脫頻率時，選用二維計算更合理。

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