高雷諾數下方柱繞流特性的數值模擬

周 強, 廖海黎, 曹曙陽

(1. 西南交通大學風工程四川省重點實驗室, 四川 成都 610031; 2. 同濟大學土木工程防災國家重點實驗室, 上海 200092)

長期以來,方柱繞流一直是風工程、鈍體空氣動力學、海洋工程等研究領域重要的、基礎性的熱點問題.在實際工程應用中,如橋梁的橋塔和橋墩、高層建筑、海洋平臺等具有方柱形狀的構件,在風或者水流作用下,會產生有規律的大尺度旋渦脫落現象,致使結構發生風致振動,進而導致其受損或發生疲勞破壞.因此,研究在方柱繞流問題及其尾流結構對于實際工程設計和建造具有十分重要的意義.

不同雷諾數下,流體繞過方柱呈現出明顯不同的流動特性.其中在高雷諾數下,流體經方柱前端邊角后發生明顯的流動分離,并伴隨著回流或轉捩,形成大尺度和不穩定的旋渦脫落,即與低雷諾數下流場相比,高雷諾數下的方柱繞流特性更加復雜,且表現出更加顯著的三維特性.鑒于此,若采用雷諾時均方法,很難撲捉到準確的高雷諾數下方柱繞流特性;由于計算量太大,目前直接數值模擬還只限于雷諾數較低的流動.因此,對于高雷諾數方柱繞流問題,主要的數值研究方法就是大渦模擬(large eddy simulation, LES)方法.

到目前為止,國內外學者針對方柱繞流問題開展了大量而豐富的實驗和數值模擬研究.在實驗研究方面,Lyn等[1]通過水洞實驗,采用LDV(laser doppler velocimetry)技術研究了Re=21 400下的方柱繞流特性,詳細分析了其尾流結構,同時其研究結果也廣泛被作為數值模擬的參考標準;此后國外學者Saha等[2]在風洞中對高雷諾數下方柱的繞流特性進行了深入研究;近年來,Nishimura和Taniike[3]通過試驗詳細研究了不同攻角下方柱繞流特性;Yen等[4]利用粒子圖像測速技術(particle image velocimetry, PIV)研究了方柱在4 000

綜上可見,對于高雷諾數下方柱繞流的數值模擬,因其流動的復雜性和明顯的三維特性,此時,采用LES是比較合適的選擇.然而如上文所述,大多采用LES方法模擬中展向長度的選取等細節問題還是存在不夠明確的地方,導致一些流場速度、氣動力等關鍵參數的脈動項與實驗結果存在一定差別.為此,本文基于開源CFD軟件OpenFoam平臺,采用基于動態亞格子模型的大渦模擬方法,在Re=22 000情況下,對均勻來流條件下的不同網格系統、展向長度的方柱繞流問題進行了三維數值模擬,首先對數值結果的準確性進行驗證,然后重點分析斷面周圍網格以及長徑比對尾流區流場結構的影響,并給出不同長徑比下氣動力的展向空間相關性.

1 控制方程與建模

1.1 控制方程

經過濾波后,大渦模擬的連續性和N-S方程為

(1)

(2)

(3)

式中：τkk為應力，k=1,2,3;δij為 Kroneker符號;Sij為變形速率張量;vSGS為亞格子模型渦粘系數,

(4)

(5)

其中:Lkk為亞格子Leonard 應力;Lij和Mij的表達式分別為

其中，Δ為格子過濾尺寸標準.

1.2 建模與數值算法

(1) 計算區域

本文共計算3個算例,其中Case 1和Case 2的計算區域為30D×18D×4D(D為方柱斷面寬度),Case 3則為30D×18D×8D.方柱模型上游來流區域長度為10D,下游尾流區長度為20D,即x方向上的長度為30D;模型離上下邊界各為9D,即y方向上的寬度為18D;為了研究長徑比的影響,展向長度分別取4D和8D兩種長度,如圖1所示.

圖1 計算區域示意圖及邊界條件Fig.1 Computational domain and boundary conditions

(2) 邊界條件(如圖1所示)

方柱表面邊界條件:采用無滑移光滑固壁邊界;

入口邊界條件:采用均勻來流速度條件,U=10 m/s;擬壓(pseudo-pressure)φ采用Neumann條件;

出口邊界條件:速度滿足格式為

?/?t+c?/?x=0

的對流邊界條件,擬壓取為0,這里c為任意常數;

上下邊界條件:采用自由流動邊界條件;

展向邊界條件:采用周期邊界條件.

(3) 網格劃分

(4) 數值算法

在空間離散上,本文采用NVD(normalized variable diagram)格式,從而避免中心差分引起的非物理性振蕩.對于時間的離散格式,采用二階完全隱式法,以獲得穩定并準確的結果.為避免計算中的數值振蕩問題,本文在數值算法上采用了MIM(momentum interpolation method)方法[14];對N-S方程組采用Piso求解算法,并將連續性方程的收斂值設置為1×10-6;無量綱時間步長設置為Δt*=ΔtU/D=2.5×10-4,以確保柯朗數小于1,且保證每步計算可在20次迭代內完成.

圖2 整個計算區域內網格劃分情況Fig.2 Grid system

(a) Case 1(b) Case 2和Case 3圖3 方柱表面周圍網格劃分情況Fig.3 Grid around square cylinder

2 數值驗證

本節針對基于時間積分的平均積分分量,與其他實驗和數值結果進行了比較,以確保本文數值計算的準確性.

如表1所示,本文3個算例得到的St均與Lyn等[1]和Norberg[15]的實驗值比較吻合;平均阻力系數CD與Sohankar等[5]的LES結果間誤差小于1%,兩者結果十分接近.

需要指出的是,表1中計算氣動力系數所采用的特征長度均為方柱斷面寬度D.由此可見,本文進行的數值模擬是準確和可靠的.

表1 平均積分分量對比表Tab.1 Comparison of integral parameters

3 結果與討論

3.1 氣動力特性

同時由表1可以發現,與Case 1中各項氣動力參數的結果相比,Case 2和Case 3結果更加準確,即在方柱表面附近先布置厚度均勻的正交網格,可有效改善數值結果.因此以下將主要針對后兩者的結果進行說明.

(a) 平均值

(b) 脈動值圖4 方柱表面壓力分布Fig.4 Surface pressure distribution on square cylinder

由圖4可見,本文結果與其他實驗和數值結果基本吻合;且與Oka等[17]的數值結果相比,本文更加接近Nishmura等[3]的實驗值.如圖4(b)所示,與Case 2相比,Case 3中的脈動值與實驗值更加吻合,特別是在方柱背面區域和靠近后方邊角附近區域,也就是旋渦主要作用區域.這再次表明氣動力系數的脈動值對展向長度更敏感.

3.2 流場特性

圖5給出并比較了中心平面(y=0的平面)上平均流向速度的分布情況.

圖5 中心線上的平均流向速度分布Fig.5 Stream-wise velocity at centreline

在尾流區,本文的結果較Wang等[8]的結果更加與實驗值一致.但可以發現各個結果間的Lr偏差較小,而Umin/U比值偏差更大.其原因在于Umin/U較Lr對長徑比更加敏感,因此,Case 2和Case 3的Umin/U比值也因為展向長度的不同而出現不同的結果.與此類似,Norberg[18]在研究圓柱繞流時,對比不同的長徑比和阻塞率的試驗結果后,也發現Umin對長徑比更加敏感.

流動速度的脈動值和雷諾應力分布都是流場結構中重要的湍流特征,因此可通過對速度脈動值和雷諾應力的分析,研究方柱尾流結構中的湍流特性.

(a) 平均流向速度(b) 平均豎向速度圖6 尾流區不同位置處平均速度分布Fig.6 Mean velocity profile in the wake

由于旋渦的交替出現使得尾流區流向速度脈動值的峰值不在中心處,而是在兩側出現雙峰值,如圖7(a)所示.這說明在尾流區存在明顯卡門渦街現象.然而隨著流動向下游發展,峰值對應的尖角消失,變得更圓滑,甚至消失.這說明旋渦在下游發展中逐漸減弱.如圖7(b)所示,尾流區不同位置處的脈動豎向速度均在在中心線上達到峰值,即存在單個峰值.因此表明交替出現的旋渦對流向和豎向上速度脈動值的影響是不同的.如圖7(c)所示,由于交替旋渦的出現,雷諾應力(u′v′)呈現關于y=0軸的反對稱性.此外以上三者均在位置x/D=1.5處達到峰值,然后隨著流向而逐漸衰減,這表明旋渦強度最大值對應的位置在x/D=1.5附近,即回轉長度大約為1.5D,再次證明了表1中得到的回轉長度是合理準確的.

同時可以發現,本文Case 2和Case 3結果均與其他結果基本吻合,但Case 3的結果更加接近Lyn等[1]實驗值.由此可見,速度的脈動值對展向長度,即長徑比更加敏感,因此為保證得到準確的湍流流場特性,需要選擇合理長徑比.

為了進一步研究展向長度(長徑比)對流場的影響,本文在方柱表面不同位置上,沿展向布置了一系列監測點,并通過分析得到其壓力的相關性系數,如圖8所示.

由于展向采用周期邊界條件,因此壓力的相關性系數在兩側均等于1.同時可以發現隨著長徑比的增加,背面和側面的壓力相關性均減弱.背面壓力相關性的減弱導致方柱的阻力系數的平均值和脈動值均有所減小,如表1所示;而側面壓力的減弱使得升力的脈動值明顯減小,從而使得其值與實驗結果更加吻合.因此,在對方柱繞流進行數值模擬時,為保證獲取準確合理的平均和脈動氣動力,以及準確的流場湍流特性,取8D的展向長度是十分必要的.

(a) 脈動流向速度(b) 脈動豎向速度(c) 雷諾應力(u′v′)圖7 尾流區不同位置處的脈動速度分布Fig.7 Fluctuating wake velocity

(a) 方柱背面不同位置處(b) 方柱側面不同位置處圖8 方柱表面壓力相關性系數Fig.8 Pressure correlation on cylinder surface

4 結 論

在典型高雷諾數Re=22 000下,本文基于開源軟件Openfoam平臺,采用大渦模擬方法,對不同網格劃分體系、不同展向長度下的方柱繞流進行三維數值模擬研究,在驗證了數值結果的準確性基礎上,深入分析了其尾流的流場結構和氣動力特性,并給出了展向長度(長徑比)對流場的影響,得到以下結論:

(1) 在方柱表面附近先布置厚度均勻的正交網格,可獲得更加準確的氣動力結果.

(2) 在方柱繞流中,回轉長度是尾流結構的重要參考指標,其值大小直接影響著方柱氣動力的表現,特別是斷面的阻力系數.

(3) 高雷諾數下,展向長度對回轉長度影響較小,但對尾流流向速度分布中的Umin/U比值影響較大,即Umin/U比值對長徑比更為敏感.

(4) 長徑比的大小直接影響方柱氣動力的表現,特別是展向相關性,這使得氣動力的脈動值對長徑比也更加敏感,故為保證計算結果的準確性,特別為了獲取更準確的湍流特性,展向長度取為8D是更為合理和有效.

[1] LYN D, EINAV S, RODI W, et al. A laser-Doppler velocimetry study of ensemble-averaged characteristics of the turbulent near wake of a square cylinder[J]. Journal of Fluid Mechanics, 1995, 304: 285-319.

[2] SAHA A, MURALIDHAR K, BISWAS G. Experimental study of flow past a square cylinder at high Reynolds numbers[J]. Experiments in Fluids, 2000, 29(6): 553-563.

[3] NISHIMURA A, TANIIKE Y. Fluctuating wind forces of a stationary two dim. square prism[C]∥Proceedings of 16th National Symposium on Wind Engineering. Tokyo: Springer, 2000: 255-260.

[4] YEN S C, YANG C W. Flow patterns and vortex shedding behavior behind a square cylinder[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2011, 99(8): 868-878.

[5] SOHANKAR A, NORBERG C, DAVIDSON L. Numerical simulation of flow past a square cylinder[C]∥The 3rd ASME/JSME Joint Fluids Engineering Conference. San Francisco: American Society of Mechanical Engineers, 1999: 18-23.

[6] OKAJIMA A. Strouhal numbers of rectangular cylinders[J]. Journal of Fluid Mechanics, 1982, 123: 379-398.

[7] SRINIVAS Y, BISWAS G, PARIHAR A, et al. Large-eddy simulation of high Reynolds number turbulent flow past a square cylinder[J]. Journal of Engineering Mechanics, 2006, 132(3): 327-335.

[8] WANG G, VANKA S. Large-eddy simulations of high Reynolds number turbulent flow over a square cylinder[J]. Dept. of Mechanical and Industrial Engineering Rep. No. CFD, 1996: 96-02.

[9] RODI W, FERZIGER J, BREUER M, et al. Status of large eddy simulation: results of a workshop[J]. Transactions-American Society of Mechanical Engineers Journal of Fluids Engineering, 1997, 119: 248-262.

[10] VOKE P R. Flow past a square cylinder: test case LES2[J]. Direct and Large Eddy Simulation Ⅱ, Spring, 1997, 5: 335-373

[11] GRIGORIADIS D, BARTZIS J, GOULAS A. LES of the flow past a rectangular cylinder using the immersed boundary concept[J]. International Journal for Numerical Methods in Fluids, 2003, 41(6): 615-632.

[12] 張偉,葛耀君. 方柱繞流粒子圖像測速試驗與數值模擬[J]. 同濟大學學報:自然科學版,2009,37(7): 857-861,892.

ZHANG Wei, GE Yaojun. Particle image velocimetry study and numerical simulation of turbulent near wake of square cylinder[J]. Journal of Tongji University: Natural Science, 2009, 37(7): 857-861, 892.

[13] LILLY D K. A proposed modification of the Germano subgrid-scale closure method[J]. Physics of Fluids A: Fluid Dynamics, 1992, 4(3): 633-635.

[14] CHOI S K. Note on the use of momentum interpolation method for unsteady flows[J]. Numerical Heat Transfer: Part A: Applications, 1999, 36(5): 545-550.

[15] NORBERG C. Flow around rectangular cylinders: pressure forces and wake frequencies[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1993, 49(1/2/3): 187-196.

[16] 周強,曹曙陽,周志勇. 亞臨界雷諾數下圓柱體尾流結構的數值模擬[J]. 同濟大學學報:自然科學版,2013,41(1): 33-38.

ZHOU Qiang, CAO Shuyang, ZHOU Zhiyong. Numerical studies of wake characteristics on a circular cylinder at sub-critical reynolds number[J]. Journal of Tongji University: Natural Science, 2013, 41(1): 33-38.

[17] OKA S, ISHIHARA T. Numerical study of aerodynamic characteristics of a square prism in a uniform flow[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2009, 97(11): 548-559.

[18] NORBERG C. An experimental investigation of the flow around a circular cylinder: influence of aspect ratio[J]. Journal of Fluid Mechanics, 1994, 258: 287-316.