高Re數(shù)圓柱繞流二維RANS模擬適用性分析

祝志文

(湖南大學(xué) 風(fēng)工程試驗研究中心，長沙 410082)

圓柱繞流即為經(jīng)典的流體力學(xué)問題，亦為流場配置最簡單、最基本的鈍體繞流。由于圓柱繞流隨流動Re數(shù)的變化呈現(xiàn)極復(fù)雜的變化特征，因此吸引諸多研究者以期從流動試驗或數(shù)值模擬上解釋其流動形態(tài)與機理［1］。以往因受計算資源限制，圓柱體繞流數(shù)值模擬通常在較低的Re數(shù)下進行［2］。雖CFD模擬可用于橋梁抗風(fēng)［3-4］，但斜拉橋拉索在高風(fēng)速下Re數(shù)一般在105以上，屬于高Re數(shù)范圍，相關(guān)數(shù)值模擬報導(dǎo)較少。眾所周知，圓柱繞流隨流動Re數(shù)、來流條件及圓柱表面粗糙度的變化呈現(xiàn)極復(fù)雜的變化特征，尤其阻力系數(shù)呈現(xiàn)較顯著的Re數(shù)效應(yīng)，如圖1所示［5］，其中Re＜2×105稱為亞臨界區(qū)，阻力系數(shù)的Re數(shù)效應(yīng)并不明顯。Re＞5×105稱為超臨界區(qū)，圓柱阻力系數(shù)隨Re數(shù)的增大而緩慢增長。2×105＜Re＜5×105通常被定義為圓柱的臨界Re數(shù)區(qū)。在斜拉索絕大部分工作狀態(tài)風(fēng)速下，其流動Re數(shù)可能均落在圓柱體的臨界Re數(shù)區(qū)，其阻力系數(shù)隨Re數(shù)的變化而快速變化，合理確定該阻力系數(shù)會有較大困難。

圖1 圓柱體阻力系數(shù)隨Re數(shù)變化曲線Fig.1 Drag coefficient of circular cylinder Vs.Re number

本文以直徑D=0.12 m圓柱為研究對象，通過數(shù)值模擬獲得圓柱氣動力參數(shù)與漩渦脫落St數(shù)?；贑FD模擬的質(zhì)量控制要求［6］，開展網(wǎng)格無關(guān)與計算時間步長無關(guān)檢查。通過獲得圓柱多個工作風(fēng)速下的阻力系數(shù)與St數(shù)，與相關(guān)文獻(xiàn)研究結(jié)果對比，評價二維RANS數(shù)值模型在獲得高Re數(shù)圓柱繞流氣動特性上的可行性，澄清CFD研究認(rèn)識誤區(qū)，為高Re數(shù)圓柱氣動特性數(shù)值研究提供建議。

1 數(shù)值方法

1.1 控制方程

繞圓柱斷面非定常二維不可壓流動的雷諾時均N-S方程為:

其中:ρ，μ分別為空氣密度與分子粘性;p，t分別為靜壓與時間;ui，u'i和uj，u'j(i=1，2;j=1，2)分別為氣流沿坐標(biāo)軸xi(i=1，2)的平均與脈動速度;為雷諾應(yīng)力。

上述雷諾應(yīng)力的引入使控制方程不封閉，需引入湍流模型以求解。若基于渦粘假設(shè)，可將雷諾應(yīng)力表示為:

其中:μt=ρCμk2/ε為湍流粘性;Cμ為經(jīng)驗常數(shù);k，ε分別為湍動能及耗散率，通過求解湍流模型方程確定。

標(biāo)準(zhǔn)k-ε模型會過高估計流動滯點區(qū)域的湍動能，不適用于風(fēng)工程問題的數(shù)值模擬［6］。本文采用綜合標(biāo)準(zhǔn)k-ε模型與k-ω湍流模型的SSTk-ε湍流模型。該模型結(jié)合上述兩湍流模型優(yōu)點，拋棄其缺點。即利用標(biāo)準(zhǔn)k-ε模型適合剪切層模擬而k-ω模型適合近壁區(qū)模擬優(yōu)點，從而通過設(shè)定一混合函數(shù)，使k-ω模型能在邊界層內(nèi)靠近壁面使用，而邊界層外使用kε模型求解。相關(guān)研究認(rèn)為，對分離點附近邊界層內(nèi)的非平衡流動，k-ωSST的模擬結(jié)果［7］優(yōu)于標(biāo)準(zhǔn)k-ε模型。

1.2 計算域與計算參數(shù)描述

采用如圖2所示形狀的計算域，以便能采用分區(qū)網(wǎng)格劃分方案并形成繞圓柱的高質(zhì)量正交四邊形單元。采用大的計算域，圖2中圓柱左側(cè)Z1區(qū)為半圓形計算域，該外邊界為計算域入口，入口至圓柱中心距離為26D。圓柱右側(cè)計算域由Z2，Z3兩個區(qū)組成，其中Z2采用與Z1類似的半圓形計算域。Z3區(qū)右側(cè)邊界為計算域出口，至圓柱中心距離為40D。因此上下邊界至圓柱中心距離同樣為26D，對應(yīng)的堵塞度小于2%。上述邊界至圓柱中心距離由試算獲得，可避免在外邊界施加的邊界條件對計算域內(nèi)數(shù)值計算影響。

Z1、Z2區(qū)均采用結(jié)構(gòu)化貼體四邊形網(wǎng)格，見圖3。Z3區(qū)為非結(jié)構(gòu)四邊形網(wǎng)格。圓柱周向等分為140個網(wǎng)格。為進行網(wǎng)格無關(guān)性檢查，最靠近圓柱表面第一個網(wǎng)格點至物面距離h0分別為0.000 02 m，0.000 01 m，0.0000 05 m，分別稱為粗網(wǎng)格 G1、細(xì)網(wǎng)格 G2、最細(xì)網(wǎng)格G3。沿物面法向，采用1.06的網(wǎng)格高度生長率，以保證在物面等流動變量變化梯度大的位置獲得高分辨率網(wǎng)格。Z2、Z3區(qū)公共邊兩側(cè)網(wǎng)格尺度保持接近。對三套網(wǎng)格系統(tǒng)，Z3區(qū)網(wǎng)格劃分完全相同，區(qū)別在于Z1、Z2區(qū)網(wǎng)格數(shù)量及網(wǎng)格分辨率。因此在圓柱周圍及尾跡區(qū)的大范圍內(nèi)能獲得較高質(zhì)量的正交四邊形網(wǎng)格。三套網(wǎng)格系統(tǒng)各自總單元數(shù)N見表1。

圖2 計算域分區(qū)Fig.2 Partition in computational domain

圖3 圓柱周圍網(wǎng)格Fig.3 Grids arrangement around cylinder

表1 網(wǎng)格劃分參數(shù)Tab.1 Grid key parameters

對計算域外邊界，在入口處采用速度入口條件，水平向速度等于來流速度，垂直水平向速度等于零，來流湍流度取0.25%。在計算域出口采用流動出口條件，即沿出口邊界的法線方向，速度梯度等于零。計算域的上、下邊界均采用對稱邊界條件，即垂直于該邊界方向，流動變量梯度為零。在靜止圓柱表面，采用無滑移邊界條件，速度與湍流度均為零。

圖4 G1網(wǎng)格圓柱表面的Y+值分布Fig.4Y+value on circular cylinder surface of G1

對圓柱繞流的非定常計算，控制方程時間離散采用二階隱式格式，空間離散采用二階迎風(fēng)格式。壓力方程與動量方程解耦采用SIMPLEC算法及欠松弛迭代。在網(wǎng)格無關(guān)及時間步長無關(guān)檢查中，來流風(fēng)速設(shè)為14.6 m/s，對應(yīng)的圓柱流動 Re數(shù)為1.2 ×105。計算0.000 05 s～0.000 5 s范圍內(nèi)多個時間步大小。對每一時間步數(shù)值計算，每一子步進行50次迭代，當(dāng)動量方程與湍流方程的殘差小于1×10-5時，認(rèn)為此時間步迭代計算已收斂。氣動力系數(shù)及其它參數(shù)的獲取，可用足夠多時間步進行數(shù)值計算，充分剔除初始計算影響，氣動力系數(shù)時程及趨勢基本穩(wěn)定后開始記錄。三套網(wǎng)格對應(yīng)的圓柱表面最大Yplus(Y+max)(表1)，因此均滿足SSTk-ω湍流模型對物面Yplus值要求。其中粗網(wǎng)格系統(tǒng)G1在圓柱表面的Yplus分布見圖4。本文數(shù)值模擬均在CFD專用軟件Fluent 3.2.26中實現(xiàn)。

1.3 網(wǎng)格與時間步長無關(guān)檢查

風(fēng)工程數(shù)值計算，尤其需重視計算模型的網(wǎng)格無關(guān)與計算時間步長無關(guān)檢查，目的為獲得與網(wǎng)格尺度、單元數(shù)量、及非定常計算時間步大小無關(guān)解，即確定數(shù)值模擬結(jié)果不再隨網(wǎng)格尺度與時間步長大小而變化的計算參數(shù)。此為風(fēng)工程數(shù)值模擬質(zhì)量保證的重要步驟。通常要求至少用三套網(wǎng)格進行繞流模擬，并比較不同網(wǎng)格與時間步長下數(shù)值模擬結(jié)果［6］。不同網(wǎng)格尺度主要集中在物面網(wǎng)格分辨率與離開物面的網(wǎng)格生長率(一般要求不大于1.2)。

定義二維圓柱阻力系數(shù)CD、升力系數(shù)CL、St數(shù)分別為:

其中:FD，F(xiàn)L分別對應(yīng)作用在圓柱截面上的阻力和升力，阻力順來流流向為正，升力垂直來流方向向上為正;U為計算域入口風(fēng)速;f為二維圓柱漩渦脫落頻率。

圖5、圖6分別為G2網(wǎng)格系統(tǒng)在時間步長dt=0.000 1 s下獲得的圓柱截面阻力與升力系數(shù)時程，均表現(xiàn)為非純諧波，說明除主頻率成分外，亦有高階諧波成分。但升力系數(shù)諧波特性較阻力系數(shù)好。升力系數(shù)主頻率為阻力系數(shù)的二分之一，也對應(yīng)圓柱的漩渦脫落頻率，由此可確定二維圓柱與漩渦脫落有關(guān)的St數(shù)。因阻力系數(shù)時程幅值波動較大，因此不同計算工況下的阻力系數(shù)平均值及根方差值可能會受所取時程位置及長度影響，反映渦脫頻率特征的St數(shù)更具參考價值。圖7、圖8為三套不同計算網(wǎng)格分別在不同時間步長下獲得的圓柱阻力系數(shù)平均值、根方差RMS及渦脫St數(shù)。

圖5 阻力系數(shù)時程Fig.5 Drag coefficient time history

圖6 升力系數(shù)時程Fig.6 Lift coefficient time history

圖7 不同工況阻力系數(shù)平均值和RMSFig.7 Mean and RMS values of drag coefficients in different simulation cases

圖8 不同工況St數(shù)Fig.8 Strouhal number in different simulation cases

由圖7、圖8知，在所有計算工況下，阻力系數(shù)RMS值均較小，而阻力系數(shù)平均值則呈現(xiàn)差別，隨時間步的增大，氣動參數(shù)差別亦增大。為有效分辨鈍體繞流非定常特征，非定常計算時間步長應(yīng)不小于渦脫周期的1/200～1/500。因此本文取時間步長不大于0.000 1 s的不同網(wǎng)格獲得的氣動參數(shù)結(jié)果。可見，三套網(wǎng)格計算所得阻力系數(shù)平均值及St數(shù)的差別很小，減小時間步長該值無明顯變化。因此，當(dāng)時間步長取不大于0.000 1 s時，三套網(wǎng)格的模擬結(jié)果無明顯變化，可認(rèn)為此時三套網(wǎng)格均能獲得與網(wǎng)格與時間步長無關(guān)的解。基于阻力系數(shù)平均值與St數(shù)在時間步長大于0.000 2 s后的變化及所用機時，本文確定細(xì)網(wǎng)格G2為后續(xù)圓柱不同風(fēng)速下的計算網(wǎng)格。

2 二維圓柱氣動特性模擬結(jié)果

為獲得二維圓柱在不同Re數(shù)下的氣動特性，據(jù)以上數(shù)值模型網(wǎng)格無關(guān)與時間步無關(guān)檢查結(jié)果，分別對該二維圓柱計算來流風(fēng)速為10 m/s、15 m/s、20 m/s、25 m/s、30 m/s、35m/s、40m/s 和 45m/s 下的繞流場。考慮時間步長影響，前四個風(fēng)速的時間步長為0.000 1 s，后四個風(fēng)速的時間步長為0.000 05 s，所有計算工況下圓柱表面的Y+值均滿足SSTk-ω湍流模型要求。圖9為計算所得圓柱截面在多個來流風(fēng)速下的阻力系數(shù)時程，可見經(jīng)初始計算波動及平穩(wěn)后阻力系數(shù)均逐漸呈規(guī)律振動。圖10為30 m/s來流風(fēng)速下圓柱周圍流場壓力等高線，從尾跡壓力等高線可觀察到圓柱尾跡的漩渦脫落情況。

圖9 不同風(fēng)速下的阻力系數(shù)時程Fig.9 Drag coefficients time histories under various wind speeds

圖10 圓柱周圍的壓力等高線云圖Fig.10 Pressure contours around cylinder

圖11為圓柱繞流在T/4、T/2、3T/4及T四個時刻的尾跡漩渦脫落圖畫，其中T/4，3T/4對應(yīng)升力系數(shù)時程分別達(dá)到負(fù)最大值與正最大值時刻，而T/2，T分別對應(yīng)升力系數(shù)時程回到零點及離開零點時刻?？梢姅?shù)值模擬能顯示圓柱尾跡漩渦產(chǎn)生、拉長、脫落及漂移過程。

圖11 一個周期內(nèi)圓柱尾跡漩渦脫落Fig.11 Vortex shedding in wake of circular cylinder during one cycle

數(shù)值模擬所得二維圓柱阻力系數(shù)平均值與渦脫St數(shù)隨Re數(shù)的變化分別見圖12、圖13。在計算的Re數(shù)范圍內(nèi)，圓柱阻力系數(shù)平均值及渦脫St數(shù)與相關(guān)文獻(xiàn)［9-10］結(jié)果變化趨勢相同，即本文數(shù)值模型能在一定程度上反映氣動參數(shù)的Re數(shù)效應(yīng)。數(shù)量上，無論阻力系數(shù)平均值或渦脫St數(shù)，均與相關(guān)文獻(xiàn)結(jié)果存在較大差異。Re＜3.55×105時，本文數(shù)值模擬所得阻力系數(shù)平均值小于風(fēng)洞試驗結(jié)果［10］，而 Re＞3.55 ×105時，本文數(shù)值模擬所得可能會大于風(fēng)洞試驗阻力系數(shù)平均值，表明本文CFD數(shù)值模型在預(yù)測圓柱定常氣動特性上的不足。另外，從圖13可見，本文CFD數(shù)值模型明顯高估圓柱渦脫St數(shù)［9］，即圓柱實際渦脫周期被明顯低估了，反映了基于RANS的CFD數(shù)值模型在預(yù)測圓柱非定常氣動特性上的不足，可能與RANS對非定常流動的時均處理，使部分非定常特性被平均化，導(dǎo)致難以準(zhǔn)確捕捉圓柱繞流的非定常特性。本文曾采用雷諾應(yīng)力模型(RSM)代替SSTk-ω湍流模型進行圓柱繞流模擬，研究結(jié)果基本相同。篇幅有限，不再贅述。

圖12 圓柱阻力系數(shù)平均值隨Re數(shù)的變化Fig.12 Cable drag coefficients vs.Re number

圖13 渦脫St數(shù)隨Re數(shù)的變化Fig.13 Strouhal number vs.Re number

3 結(jié)論

基于CFD專用軟件Fluent3.2.26與嚴(yán)格數(shù)值模擬方法及步驟，對二維RANS數(shù)值模型在預(yù)測圓柱定常與非定常氣動特性上展開研究，結(jié)論如下:

(1)本文給出的圓柱計算域分區(qū)網(wǎng)格劃分方法能在圓柱周圍及尾跡區(qū)大范圍內(nèi)形成高質(zhì)量正交四邊形網(wǎng)格。由此開展的網(wǎng)格無關(guān)與時間步長無關(guān)檢查是本文數(shù)值方法重要基礎(chǔ)。

(2)二維RANS數(shù)值模型能一定程度給出圓柱氣動參數(shù)的Re數(shù)效應(yīng)，也能給出圓柱尾跡漩渦產(chǎn)生、生長、脫落及漂移的周期性過程。

(3)二維RANS數(shù)值模型給出的圓柱阻力系數(shù)平均值與渦脫St數(shù)，均與相關(guān)文獻(xiàn)結(jié)果存在較大差異。尤其CFD模擬明顯高估了圓柱的渦脫St數(shù)，反映出二維RANS模型在預(yù)測圓柱定常與非定常氣動特性上的不足。

采用二維RANS模型難以合理預(yù)測圓柱的氣動力及渦脫特性。若研究圓柱定常與非定常特性，需嘗試采用三維模擬或其它相關(guān)方法。

［1］Nishimuraa H，Taniike Y.Aerodynamic characteristics of uctuating forces on a circular cylinder［J］.Journal of Wind Engineering and Industrial Aerodynamics，2001，89(7-8):713-723.

［2］Franke J，F(xiàn)rank W.Large eddy simulation of the flow past a circular cylinder at ReD=3 900［J］.Journal of Wind Engineering and Industrial Aerodynamics，2002，90(10):1191-1206.

［3］祝志文，陳政清，顧 明.橋梁氣彈響應(yīng)分析的時域氣動模型［J］.中國公路學(xué)報，2007，20(6):43-48.

ZHU Zhi-wen，CHEN Zheng-qing，GU Ming.Time-domain aerodynamic model in aeroelastic response analysis of bridges，China Journal of Highway and Transport，2007，20(6):43-48.

［4］祝志文，陳政清.數(shù)值模擬橋梁斷面的顫振導(dǎo)數(shù)和顫振臨界風(fēng)速［J］.中國公路學(xué)報，2004，15(4):41-50.

ZHU Zhi-wen，CHEN Zheng-qing.Numerical simulations for aerodynamic derivatives and critical flutter velocity of bridge deck［J］.China Journal of Highway and Transport，2004，15(4):41-50.

［5］Simiu E，Scanlan R H.Wind effects on structures-an introduction to wind engineering［M］.New York:Wiley，Third Edition，1996.

［6］Casey M，Wintergerste T.Best practice guidelines:ercoftac special interest group on“quality and trust in industrial CFD”［M］.Fluid Dynamics Laboratory，Sulzer Innotec，2000.

［7］ Franke J，Hirsch C，Jensen A G.，et al.Recommendations on the use of CFD in wind engineering［C］.Proceeding of International Conference on Urban Wind Engineering and Building Aerodynamics.Belgium，2004.

［8］Poulin S，Larsen A.Drag loading of circular cylinders inclined in the along-wind direction［J］.Journal of Wind Engineering and Industrial Aerodynamics，2007，95(9-11):1350-1363.

［9］Norberg C.Fluctuating lift on a circular cylinder:review and new measurements［J］.Journal of Fluids and Structures，2003，17(1):57-96.

［10］林志興，楊立波，李文勃.斜拉索順橋向風(fēng)阻系數(shù)的試驗研究［J］.鄭州大學(xué)學(xué)報(工學(xué)版)，2005，26(1):38-41.

LIN Zhi-xing，YANG Li-bo，LI Wen-bo.Experimental study on drag coefficientsofstay-cablescorresponding towind direction along the bridge centrallin［J］.Journalof Zhengzhou University(Engineering Science)，2005，26(1):38-41.