矩形柱體氣動性能的流場機理研究

杜曉慶，林偉群，代 欽

(1. 上海大學(xué) 土木工程系，上海 200444；2. 上海大學(xué) 風(fēng)工程和氣動控制研究中心，上海 200444；3. 上海大學(xué) 上海市應(yīng)用數(shù)學(xué)和力學(xué)研究所，上海 200444)

超高層建筑的抗風(fēng)設(shè)計是結(jié)構(gòu)設(shè)計中需要重點考慮的因素.由于矩形截面在超高層建筑中經(jīng)常被采用，二維矩形柱的氣動性能和繞流場特性受到廣泛關(guān)注[1-16].研究表明，在一定的雷諾數(shù)下，矩形柱的氣動性能和流場結(jié)構(gòu)受截面寬厚比B/D(其中B為順風(fēng)向?qū)挾?、D為橫風(fēng)向厚度)的控制[1-4].矩形柱的平均阻力系數(shù)會在B/D=0.6附近達到峰值[1]，Strouhal數(shù)(St數(shù))則會在B/D=2.8和B/D=6附近會發(fā)生不連續(xù)的跳躍現(xiàn)象[1-3].

研究者主要通過風(fēng)洞試驗對矩形柱的氣動性能進行研究，對矩形柱繞流場的研究較少.隨著計算機性能的提高和計算方法的改進，CFD方法逐漸被用于研究矩形柱繞流問題.YU和KAREEM較早開展了對矩形柱的研究[5-6]，現(xiàn)有研究中對B/D=1的方柱[7-8]和B/D=5的矩形柱[9-10]的數(shù)值模擬研究相對較多.最近SOHANKAR[11]和YU等[12]采用大渦模擬方法分別對寬厚比為B/D=0.4～4.0和B/D=0.3～7.0的矩形柱進行了較為系統(tǒng)的研究.與雷諾平均法相比，大渦模擬雖然計算精度較高，但其計算量大.SHIMADA和ISHIHARA[13]采用雷諾平均法和修正的k-ε湍流模型研究了寬厚比B/D=0.6～8的矩形柱，其數(shù)值模擬結(jié)果與試驗結(jié)果吻合較好.不少研究發(fā)現(xiàn)[11-13]，矩形柱的St數(shù)在B/D=2附近達到最小值，并且B/D=2矩形柱的流場特性也與其他寬厚比矩形柱有明顯差異，其確切的機理尚未得到澄清.

本文采用綜合了標準k-ε模型和k-ω湍流模型優(yōu)點的k-ωSST湍流模型，在雷諾數(shù)Re=2.2×104時，研究了寬厚比B/D=1～4矩形柱的氣動性能和流場特性，探究了流場結(jié)構(gòu)、表面風(fēng)壓和氣動力之間的內(nèi)在關(guān)系，基于表面摩擦力系數(shù)估算了B/D=3和4的分離泡長度，并重點分析了寬厚比B/D=2矩形柱的瞬態(tài)流場特性，澄清其St數(shù)較其他寬厚比矩形柱體小的流場機理.

1 數(shù)值方法和計算模型

1.1 控制方程和湍流模擬

本文采用基于k-ωSST湍流模型的雷諾平均法N-S方程進行模擬.

連續(xù)性方程：

(1)

動量方程：

(2)

(3)

式中：ut=ρCμk2/εμt=ρCμk2/ε為湍流黏性；Cμ為經(jīng)驗常數(shù)；k和ε分別為湍動能及耗散率，需要通過求解湍流模型方程來確定.

本文采用綜合了標準k-ε模型和k-ω湍流模型的k-ωSST湍流模型.因為標準k-ε模型適合剪切層模擬而k-ω模型適合近壁區(qū)模擬，所以k-ωSST湍流模型通過設(shè)定一個混合函數(shù)，使得k-ω模型能在邊界層內(nèi)靠近壁面使用，而邊界層外則利用k-ε模型求解.

1.2 計算模型

本文采用四種寬厚比的矩形柱體，分別為B/D=1, 2, 3和4.來流為均勻來流，風(fēng)速為Uo，風(fēng)攻角為0°，雷諾數(shù)為Re=2.2×104(根據(jù)矩形柱厚度D和來流風(fēng)速Uo計算得到).

圖1為本文采用的計算域及邊界條件.本文計算采用結(jié)構(gòu)化網(wǎng)格，以矩形柱的中心做為半圓域的圓心，半圓域的直徑為40D，順風(fēng)向計算域的總長為45D，展向長度為B，阻塞率為2.5%.網(wǎng)格在靠近矩形柱體處加密，第一層網(wǎng)格高度為0.001D，近壁面網(wǎng)格增長率為1.08.計算模型局部網(wǎng)格如圖2所示.

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

圖2 計算模型的平面網(wǎng)格Fig.2 Plane computation grid scheme

計算域采用速度入口邊界條件，自由出口邊界條件，展向采用周期性邊界條件，矩形柱表面采用無滑移壁面邊界條件.本文采用SIMPLEC算法來求解壓力-速度場耦合方程；速度和壓強的離散采用二階迎風(fēng)格式；時間離散方法采用二階隱式格式.

定義矩形柱體的升力系數(shù)、阻力系數(shù)為

(4)

式中：FL和FD分別為單位長矩形柱的氣動升力和氣動阻力，ρ為空氣密度.

矩形柱表面的風(fēng)壓系數(shù)和摩擦力系數(shù)分別為

(5)

式中：p-po為當(dāng)?shù)仫L(fēng)壓和遠前方上游壓力之差，τw為表面粘性切應(yīng)力.

2 計算結(jié)果及分析

2.1 表面風(fēng)壓系數(shù)

圖3為各種寬厚比矩形柱表面的平均風(fēng)壓系數(shù)分布圖，圖中也列出了文獻值[6][13]進行對比，由圖可見，本文結(jié)果與文獻結(jié)果的總體變化趨勢是相同的，結(jié)果吻合較好.

各寬厚比矩形柱的迎風(fēng)面(a-b)風(fēng)壓相近，但其側(cè)面(b-c)和背風(fēng)面(c-d)風(fēng)壓有明顯差異.寬厚比B/D=1和2時，矩形柱側(cè)面的風(fēng)壓分布比較平穩(wěn)，側(cè)面負壓絕對值隨著寬厚比的增大有減小趨勢；而寬厚比B/D=3和4時，側(cè)面的風(fēng)壓分布有較大的變化，靠近前角(b點)的側(cè)面負壓較強，而靠近后角(c點)的側(cè)面負壓有逐漸減弱的趨勢，出現(xiàn)壓力恢復(fù)區(qū).而矩形柱背風(fēng)面的負壓絕對值則隨著寬厚比的增大逐漸減小，風(fēng)壓系數(shù)從B/D=1的-1.3逐漸變?yōu)锽/D=4的-0.5左右.

2.2 氣動力系數(shù)

圖4和圖5分別為不同寬厚比矩形柱體的平均阻力系數(shù)和脈動升力系數(shù)，圖中也列出了文獻中的風(fēng)洞試驗值[1,3,17-18]和數(shù)值模擬值[5-6,13,19]進行比較.

圖3 平均風(fēng)壓系數(shù)Fig.3 Mean pressure coefficient distribution

圖4 矩形柱體的平均阻力系數(shù)Fig.4 Mean drag coefficient of the rectangular cylinders

由圖4可見，本文計算得到的平均阻力系數(shù)與文獻結(jié)果的總體變化趨勢是相同的.除了B/D=1的方柱的平均阻力系數(shù)稍小于文獻結(jié)果，這可能是本文采用的湍流模型造成的；其他寬厚比矩形柱體的平均阻力系數(shù)都處在之前文獻值的區(qū)間范圍內(nèi).從總體上看，隨著寬厚比B/D的增大，矩形柱平均阻力系數(shù)有逐步減小的趨勢.由于矩形柱的阻力主要是其迎風(fēng)面和背風(fēng)面的壓差造成的，因而阻力系數(shù)的這種變化趨勢也可以從圖3的表面風(fēng)壓分布看出：隨著寬厚比的增大，矩形柱迎風(fēng)面的風(fēng)壓基本保持不變，而背風(fēng)面的負壓絕對值卻逐漸減小.

由圖5可見，本文計算得到的脈動升力系數(shù)與文獻總體變化的趨勢是相同的.隨著寬厚比的增大，脈動升力系數(shù)會迅速下降，其均方根值會從寬厚比B/D=0.6的2.8減小至B/D=4的0.2左右.但從圖中也可以看出：對于寬厚比為B/D=2～4的矩形柱，不同研究者得到的脈動升力系數(shù)偏差較大，這可能是因為這一寬厚比區(qū)間正好是兩種不同流態(tài)發(fā)生突變的范圍，因而對于試驗和數(shù)值模擬的條件比較敏感，特別是來流湍流度和雷諾數(shù).

圖5 矩形柱體的均方根升力系數(shù)Fig.5 RMS lift coefficient of the rectangular cylinders

2.3 Strouhal數(shù)

圖6為不同寬厚比矩形柱體的St數(shù)(St=fD/Uo，其中f為渦脫頻率)，圖中也列出了文獻中的風(fēng)洞試驗值[2-3]和數(shù)值模擬值[5-6,13,19].

由圖可見，本文計算值隨寬厚比的變化趨勢與文獻結(jié)果是一致的和接近的，但數(shù)值上總體有偏小的趨勢，這可能是k-ωSST湍流模型的局限性造成的.隨著寬厚比的增大，St數(shù)的變化比較復(fù)雜.總體來說，St數(shù)在B/D=2.5左右出現(xiàn)最小值，而在B/D=3附近達到最大值，即St數(shù)在B/D=2～3之間會發(fā)生突變.

圖6 不同寬厚比矩形柱體的Strouhal數(shù)均風(fēng)壓系數(shù)Fig.6 Strouhal numbers of the rectangular cylinders with different aspect ratios

圖7 時間平均流線圖Fig.7 Time-averaged streamline

2.4 平均流線圖流場

圖7為不同寬厚比矩形柱的時間平均流線圖，并基于計算結(jié)果繪制了流場結(jié)構(gòu)圖.

由圖7可見，不同寬厚比矩形柱的側(cè)面和背風(fēng)面均存在兩個尺度較大的回流區(qū)，但這兩個回流區(qū)隨著寬厚比的增大會發(fā)生顯著變化.當(dāng)B/D=1和2時，分離的剪切層并沒有再附到矩形柱的側(cè)面上，側(cè)面和背風(fēng)面的兩個回流區(qū)在矩形柱的后角處相互連接，并在側(cè)面后角有一尺度很小、方向相反的回流；B/D=2的背風(fēng)面回流區(qū)的長度和寬度均大于其他工況，這是其St數(shù)遠低于其他寬厚比矩形柱的原因.當(dāng)B/D=3和4時，側(cè)面和背風(fēng)面的兩個回流區(qū)被后角分隔開，形成兩個相互獨立的回流區(qū)，并且側(cè)面回流區(qū)的尺度很大，而背風(fēng)面的回流區(qū)很?。浑S著寬厚比從B/D=3增大到4，側(cè)面回流中心會向上游移動，而尾流回流中心則有向下游移動的趨勢.

為了進一步定量研究流場結(jié)構(gòu)，圖8給出了不同寬厚比矩形柱側(cè)面的平均摩擦力系數(shù)分布.從圖中可見，矩形柱側(cè)面的摩擦力系數(shù)大部分為負值，在靠近后角處有一大小不一的正值區(qū)域，但摩擦力系數(shù)出現(xiàn)正負交替的原因不同.對于B/D=1和2的矩形柱，摩擦力系數(shù)正負交替是由于側(cè)面后角處的小回流區(qū)造成的，回流區(qū)的長度分別為0.04D和0.05D.對于B/D=3和4的矩形柱，摩擦力系數(shù)的正負交替則是因為存在剪切層再附造成的，在側(cè)面回流區(qū)的壁面摩擦力為負值，在再附點后側(cè)的摩擦力則為正值；根據(jù)TAFTI和VANKA的定義[20]，平均摩擦力系數(shù)由負到正的位置即是時間平均的分離流再附點，從平均流線圖中也可明顯地看到存在分離泡和再附流；當(dāng)B/D=3時，再附點到前角分離點的長度(即為分離泡的長度)約為2.90D，而B/D=4時約為3.82D.

正是由于B/D=3和4的矩形柱的側(cè)面存在分離泡，分離泡內(nèi)的側(cè)面風(fēng)壓強度較強，在后角分離流再附區(qū)出現(xiàn)風(fēng)壓的恢復(fù)區(qū)；分離泡的出現(xiàn)還會引起其旋渦脫落強度降低和尾流變窄，這也是其平均阻力系數(shù)和脈動升力系數(shù)較小、而St數(shù)較大的原因.

圖8 矩形柱體側(cè)面摩擦力系數(shù)分布圖Fig.8 Friction coefficient distribution on the lateral side of rectangular cylinders

2.5 瞬態(tài)流場

圖9為不同寬厚比矩形柱繞流場的典型時刻瞬時渦量圖(升力系數(shù)達到最大的時刻).從圖中可以看出，所有寬厚比的矩形柱都會有旋渦從兩側(cè)交替脫落，在下游形成一條交錯排列的渦街，但渦街的特性會隨著寬厚比發(fā)生劇烈變化.B/D=1的尾流寬度最大，這也說明其渦脫強度最強，B/D=3和4的由于存在分離再附現(xiàn)象因而其尾流最窄；B/D=2時尾流中相鄰旋渦的間距最長，而B/D=3時的渦間距最短，這也反映出前者的St數(shù)最小而后者最大；此外，B/D=2矩形柱尾流中的上下兩側(cè)剪切層的相互作用范圍明顯長于其他工況，呈現(xiàn)出獨特的流場特性，目前尚未有明確的解釋.

圖9 瞬時渦量圖(升力最大時刻)Fig.9 Instantaneous vorticity field at the maximum lift

為了進一步理解B/D=2時獨特流場特性的發(fā)生機理，圖11給出了4個典型時刻B/D=2矩形柱局部流場的瞬時渦量圖、瞬時風(fēng)壓場和瞬時流線圖.這4個采樣時刻在升力系數(shù)時程曲線上的位置見圖10.從圖11可見，在t1和t4時刻，B/D=2的矩形柱上側(cè)面后角有一個逆時針旋轉(zhuǎn)的旋渦；在t7和 t8時刻，B/D=2的矩形柱下側(cè)面后角有一個順時針旋轉(zhuǎn)的旋渦.這是因為矩形柱的一側(cè)剪切層在角點分離后，會在背風(fēng)面回流到矩形柱上，并在另一側(cè)后角的角點處再次分離，隨后又再附到另一側(cè)后角，在另一側(cè)后角形成一個小旋渦，旋渦旋轉(zhuǎn)的方向如局部流線圖所示，這也解釋了圖7中寬厚比B/D=2矩形柱的平均流線圖側(cè)面后角出現(xiàn)了一個小回流區(qū)的原因.從風(fēng)壓場中可以看出，模型的迎風(fēng)面除了靠近角點處均為正壓，模型的側(cè)面和背風(fēng)面為負壓.在后角旋渦的位置，局部的負壓會比其它地方明顯增大，這是后角旋渦導(dǎo)致的.

從圖11b的流線圖(第三列)中可見，在t4時刻，矩形柱上側(cè)尾流會形成一個二次渦(見箭頭所指處)；同樣，在t7時刻(見圖11(d)第三列)，下側(cè)尾流中也會形成一二次渦.這些二次渦并非是在矩形柱上脫落的渦形成的.

為了研究二次渦的形成過程，圖12進一步在升力時程曲線上選取更多的瞬時時刻(t2～t6)流線來進一步分析，這些瞬時時刻的具體位置見圖10.在t2時刻，下側(cè)尾流中已形成一個即將脫離矩形柱的旋渦，而矩形柱下側(cè)剪切層會沿著背風(fēng)面回流到上側(cè)面，下側(cè)剪切層回流的方向與上側(cè)的剪切層方向相反(見圖12(a)).在t3時刻，在下側(cè)尾流中的旋渦和上側(cè)剪切層的相互作用下，在上側(cè)尾流中會形成二次渦(見圖12(b)箭頭所在處).在t4和t5時刻，二次渦有所發(fā)展，并有朝著來流方向運動的趨勢.順時針旋轉(zhuǎn)的二次渦與下側(cè)尾流中逆時針旋轉(zhuǎn)的旋渦會發(fā)生相互作用，由于兩者的運動方向相反，二次渦的出現(xiàn)會阻止矩形柱下側(cè)尾流旋渦向下游的移動及其脫落過程(見圖12(c)(d)).二次渦在t6時刻消失(見圖12(e)).

從上述分析可知，二次渦的出現(xiàn)會延緩矩形柱尾流渦脫的速度，同時也會使矩形柱尾流中的旋渦被拉長，從而大大增大了旋渦脫落的周期，從而使得St數(shù)明顯低于其他的矩形柱.

圖10 升力系數(shù)時程曲線上的采樣時刻Fig.10 Samples times on the time history of lift coefficient

圖11 B/D=2時的瞬時渦量場、風(fēng)壓場和流線圖Fig.11 Instantaneous vorticity field, pressure field, streamline

圖12 二次渦形成過程的瞬態(tài)流線圖Fig.12 Instantaneous streamline of the formation of the secondary vortex

3 結(jié)論

本文采用雷諾平均法和k-ωSST湍流模型，對四種寬厚比的矩形柱體繞流場進行了CFD計算，并對氣動力、Strouhal數(shù)、壓力分布以及流場特性進行了分析，得到下述結(jié)論：

(1)計算得到的不同寬厚比矩形柱的表面風(fēng)壓系數(shù)、氣動力系數(shù)和Strouhal數(shù)和文獻結(jié)果吻合較好，說明采用k-ωSST湍流模型可較準確地模擬矩形柱體的繞流問題.

(2)寬厚比B/D=1和2的矩形柱沒有分離剪切層再附現(xiàn)象；B/D=3和4的矩形柱會發(fā)生分離剪切層的再附現(xiàn)象，并在側(cè)面形成分離泡，分離泡的出現(xiàn)會引起尾流寬度變窄和渦脫強度減弱，并導(dǎo)致平均阻力系數(shù)和脈動升力系數(shù)的減小.

(3)根據(jù)矩形柱側(cè)面的平均摩擦力系數(shù)，可以估算B/D=1和2矩形柱側(cè)面后角回流區(qū)的寬度分別為0.04D和0.05D；而B/D=3和4矩形柱的分離泡長度分別為2.90D和3.82D；隨著寬厚比從B/D=3增大到4，分離泡的中心位置會向上游移動.

(4)對于寬厚比B/D=2的矩形柱，其在一側(cè)分離的剪切層會在背風(fēng)面回流，并在另一側(cè)的后緣角點再次分離，矩形柱一側(cè)的尾流旋渦會與另一側(cè)的剪切層相互作用，在矩形柱后側(cè)形成二次渦，二次渦的出現(xiàn)會抑制脫落旋渦的速度，增大旋渦脫落的周期，從而使得B/D=2的矩形柱的St數(shù)明顯低于其他寬厚比矩形柱.