二維方柱繞流阻塞效應的大渦模擬

郜 陽， 全 涌， 顧 明

(同濟大學 土木工程防災國家重點實驗室， 上海 200092)

在對建筑結構進行風洞試驗研究時，阻塞效應可能導致風洞內的流場特性偏離真實流場，從而影響試驗結果的準確性.深入研究阻塞比(Blockage Ratio，以下簡稱R)對流場以及柱體氣動力的影響對保證風洞試驗結果的正確性是至關重要的.

研究人員[1-4]已采用風洞試驗的方法，就阻塞比對模型風壓分布、氣動力等氣動特性的影響展開了大量的研究.大多數的學者，如黃劍等[3-4]，均是通過對相同截面形式、不同縮尺比的模型進行風洞試驗，來研究阻塞效應.而少數學者[1]考慮到模型尺寸對風洞試驗結果的影響，通過改變在風洞內安裝的隔板的位置實現不同阻塞比，保證了模型的尺寸不變，從而對阻塞效應展開研究.盡管已有大量試驗研究，但對阻塞效應的定量認識仍然比較模糊.

利用計算流體力學(computational fluid dynamics，以下簡稱CFD)開展阻塞效應的研究最早開始于1980年代，Stafford[5]通過數值方法驗證了風洞洞壁形狀的改變對降低阻塞效應的可行性.之后，Davis[6]對兩種阻塞比(16.7%和25%)的矩形柱進行數值模擬，研究表明平均阻力系數以及Strouhal數隨著阻塞比的增大而增加.Okajima[7]、Sharma[8]和Patil[9]等通過數值模擬方法考察了方柱阻塞效應，研究了平均阻力系數和脈動升力系數隨阻塞比的變化規(guī)律.值得注意的是，這些關于阻塞效應的數值模擬均是在二維平面內建立計算域，并且局限于低雷諾數范圍內.相比于二維平面，三維空間的數值模擬更為復雜，但也更接近于實際.遺憾的是，到目前為止，僅有少數學者針對三維空間建模，系統(tǒng)地研究阻塞比的影響.基于Nakagawa[2]的方形截面柱風洞試驗，Kim[10]在三維空間建立流場，研究在約束條件下的流場特性，并與試驗結果吻合較好，但該文獻僅給出了阻塞比為20%一種阻塞比工況.

方柱作為研究鈍體繞流的常用模型，在工程中也是廣泛采用的截面形式.針對二維方柱繞流，無論是風洞試驗還是數值模擬都有很多已發(fā)表的文獻，尤其是在雷諾數Re=2.2×104條件下[11-16].本文選取同樣的截面形式及相近的雷諾數，采用三維大渦模擬(LES)，研究阻塞比(柱體迎風面面積與計算域截面面積之比)對方柱繞流的影響.在下文中，首先介紹了數值計算采用的湍流模型、邊界條件、求解器設置以及網格劃分等參數；然后通過與以往試驗及數值模擬的結果進行對比，驗證本文數值模擬計算結果的合理性；最后詳細分析了阻塞比對方柱繞流氣動力及流場的影響.

1 計算參數設置

1.1 湍流模型

本文采用大渦模擬進行數值模擬計算.通過對流體連續(xù)性方程及動量方程空間上的濾波，得到過濾后的Navier-Stocks控制方程，表達式如下：

(1)

(2)

(3)

式中：ui(i=1,2,3)是指3個方向的風速分量；p表示的壓力；v為運動粘性系數；帶有上劃線的量表示空間濾波后的流場變量；τij代表亞格子應力，在動態(tài)Smagorinsky亞格子應力模型中[17]，利用線性渦粘模型將亞格子應力的各項異性部分與最小求解尺度上的應變率聯系起來，如式(4)所示.

(4)

1.2 計算域及邊界條件

為了研究不同的阻塞比對柱體氣動力以及風場特性的影響，設置了5種阻塞比工況，不同工況之間只是計算域的寬度發(fā)生變化，計算域的寬度從24D減小到4D(D為模型迎風面寬度)，對應的阻塞比由4.2%增大25%，各工況的阻塞比以及對應的計算域尺寸設置如表 1所示.圖 1所示為計算域尺寸及邊界條件，上游入口到柱體迎風面距離為6D，柱體背風面到下游出口的距離為16D，計算域順流向長度總計為23D.以往的數值模擬，多數采用了上游4.5D～8D以及下游15D～16.5D的計算域長度[13-14,18-19]，本文的長度介于上述范圍.為了形成三維湍流結構，展向長度取4D[13,20].計算域寬度B的取值決定了阻塞比的大小，見表1.采用圖 1所示的網格區(qū)域劃分方法：在計算域橫向上劃分了多個不同的區(qū)域，其中I區(qū)為柱體附近區(qū)域，其寬度為2D，在此區(qū)域內網格比較密集，用于提高柱體附近區(qū)域渦旋捕捉的精度；II區(qū)為靠近約束邊界的區(qū)域，寬度為0.5D，在計算域的兩側設置了無滑移邊界條件，網格局部加密；通過改變III區(qū)的寬度改變計算域的寬度，得到不同的阻塞比.

表1 各阻塞比工況計算域尺寸

a 側視圖

b 俯視圖

對于5種阻塞比工況， I區(qū)和II區(qū)兩個區(qū)域內的網格劃分完全一致，以提高計算結果的可對比性.I區(qū)是柱體所在區(qū)域，在柱體表面周向均勻布置 200個網格.對于LES數值模擬，需要保證障礙物附近區(qū)域有足夠精度的網格分辨率，在柱體近壁面第1層網格高度約為D/2000，要小于0.1D/Re0.5=6.7×10-4D，向外延伸的增長率為1.13，計算得到的所有工況柱體近壁面Y+<2(Y+=uτy/v，uτ為壁面摩擦速度；y為近壁面第1層網格高度).柱體展向網格長度為0.1D.網格劃分全部采用結構化，圖 2為工況1的整體網格劃分以及柱體附近區(qū)域網格局部放大示意圖.由于精確控制每一個分區(qū)網格的尺寸，在分區(qū)分界線兩側的網格尺寸基本一致，避免相鄰網格尺寸突變過大而導致的計算誤差.

圖2 工況1計算域網格劃分示意圖

圖1也給出了計算域的各邊界條件.入口邊界條件設置為速度入口(velocity-inlet)，u=(U0,0,0)，湍流強度為0，雷諾數Re=2.2×104(根據入口平均風速以及柱體寬度尺寸計算得到)；柱體表面以及計算域的兩側邊界均采用無滑移邊界條件(Wall)；展向的邊界設置為周期邊界條件(Periodic)；出口采用自由出流邊界條件(Outflow).

數值計算采用有限體積法，采用SIMPLEC格式求解壓力速度耦合方程組，對流項采用數值耗散低的二階中心差分格式(bounded central differencing)，時間離散項都采用二階全隱式方法(bounded second-order implicit).為了保證計算最大柯朗數(Courant number)小于1，量綱為一計算時間步長Δt*為0.003 2(Δt*=Δt·U0/D，Δt為實際計算時間步，U0為來流風速).為加快計算的收斂速度，在進行LES計算之前，先采用雷諾平均定常計算(RANS)方法對流場進行初始化.在LES計算穩(wěn)定之后開始采樣，采樣時長均大于40個旋渦脫落周期，所有計算工況均采用以上求解設置.

2 結果分析

2.1 計算結果驗證

在進行不同阻塞比工況計算之前，針對工況1進行了兩組不同數量網格的計算，如表 2所示.同時在表 2中給出了各氣動力參數的計算結果.可以看出，對網格進行加密之后，各氣動力參數并沒有明顯的變化.因此本文以下工況采用網格2的計算結果，并且其他阻塞比工況在I區(qū)和III區(qū)域工況1保持一致的網格.

表2 工況1網格收斂性驗證

為了驗證本文數值模擬的準確性，將本文工況1工況計算結果與以往的試驗[11,21-22]結果以及數值模擬[13-16,23-24]結果進行對比.與本文不同的是，這些數值模擬中，兩個側面的邊界均不是無滑移邊界條件.表 3列舉了風洞試驗、LES模擬以及直接模擬(DNS)的結果，對比項包括：平均阻力系數(CD,mean)、脈動阻力系數(CD,rms)、脈動升力系數(CL,rms)、斯托羅哈數(St)、回流區(qū)長度(Lf)等.由表 3可見，本文的平均阻力系數、脈動阻力系數相比前人的研究結果略微偏大，而斯托羅哈數和回流區(qū)長度與各數值模擬以及試驗結果都比較接近.表 3中也給出了文獻中相應的模型長寬比(Ar)以及來流的湍流強度(Iu).

圖3給出了工況1沿柱體中心截面環(huán)向的風壓系數分布圖，并將DNS[15]和LES[16,24]模擬結果以及風洞試驗結果[21-22]一并示于圖中.圖 3a為平均風壓系數分布，可以看出在迎風面，本文結果與風洞試驗結果以及前人數值模擬結果吻合較好；在側面與背風面，本文結果與風洞試驗結果略有不同，但與Trias[15]和Cao[16]的數值模擬結果吻合較好.

圖3b為脈動風壓系數的分布，在迎風面以及背風面，本文計算結果與試驗結果、數值模擬結果吻合較好；在側面，本文計算結果要略大于試驗結果，而與其他數值模擬結果吻合較好.

由表3及圖3可見，不同的風洞試驗及數值模擬結果存在明顯的差異，這主要是由不同風洞實驗室的試驗條件、模型尺寸、來流湍流度等以及數值模擬的湍流模型、邊界條件、網格劃分等因素引起.而采樣時間不足也會引起結果的偏差[15].但整體上，本文的數值模擬結果與前人試驗結果及數值模擬結果均較為接近，是可靠有效的.

表3 整體氣動力特性參數對比

a 平均風壓系數

b 脈動風壓系數

2.2 阻塞效應分析

2.2.1阻塞比對柱體氣動力系數的影響

當阻塞比較大時，在方柱前緣邊角分離的氣流以及尾流會受到邊界的約束作用，從而對柱體的氣動力產生影響.柱體氣動力系數的統(tǒng)計值隨阻塞比的變化規(guī)律如圖 4所示.其中，圖 4a給出了平均阻力系數(CD,mean)、迎風面風力系數均值(Cf,mean)以及背風面風力系數均值(Cr,mean)隨阻塞比的變化規(guī)律.從中可以看出，隨著阻塞比R由4.2%增大到25%，平均阻力系數由2.35增加到3.0，增加了26.7%，同時迎風面平均風力系數略有減小，背風面風力系數(主要是吸力)均值絕對值則單調增加.這主要是由于阻塞比增大，氣流流通面積被壓縮，柱體尾流的流速增加，使得柱體背風面的吸力增加.圖 4a還給出了Kim[10]在阻塞比為20%時方柱繞流的LES計算結果，其平均阻力系數與本文工況4(R=20%)的結果完全一致.

圖 4b中給出了脈動阻力系數(CD,rms)、背風面風力系數脈動值(Cr,rms)和脈動升力系數(CL,rms)隨阻塞比的變化規(guī)律.可以看出，脈動阻力系數值與背風面風力系數脈動值基本重合，這表明，阻力的脈動特性主要是來自背風面吸力的脈動特性的貢獻.這種脈動特性會隨著阻塞比的增大略微有所增加.脈動升力系數在阻塞比由4.2%增大到25%的過程中，由1.56逐漸增加到2.24，增加了43.6%.在圖4b給出了Kim[10]在阻塞比為20%條件下方柱繞流的LES數值模擬結果，其脈動升力系數也與本文工況4(R=20%)的結果基本一致.

a 氣動力均值

b 氣動力脈動值

2.2.2阻塞比對Strouhal數的影響

柱體Strouhal數隨阻塞比的變化規(guī)律如圖 5所示，當阻塞比從4.2%增大到25%，Strouhal數從0.133增加到了0.175，增大了31%.這是由于洞壁對流場的約束作用加強，鈍體繞流的流速增加所導致的.對于阻塞比較小的兩組工況(R=4.2%和R=10%)，Strouhal數值相差很小，相差3%；同樣平均阻力系數以及脈動升力系數值也很接近，相差在3%以內，如圖 4所示.

圖5 阻塞比對Strouhal數的影響

圖6為各柱體阻力系數功率譜曲線.可以看出，對于阻塞比較小的前兩個工況1、工況2(R=4.2%和10%)，阻力系數的功率譜曲線的值偏小，沒有出現峰值，這時的阻力系數的脈動比較低.當阻塞比增大到一定值(R≥16.7%)，阻力系數的頻譜曲線會在大約2倍Strouhal數附近出現明顯的峰值，尤其工況5(R=25%)，其阻力系數頻譜曲線的峰值非常突出.黃劍等[25]對不同阻塞比的三維矩形柱進行風洞試驗時，也發(fā)現當阻塞比過大時，阻力系數頻譜曲線會在大約2倍Strouhal數處出現峰值.

圖6 阻塞比對阻力系數功率譜的影響

2.2.3阻塞比對風壓系數的影響

圖 7為各工況的柱體中心橫截面(z=0)的平均風壓系數以及脈動風壓系數分布.如圖 7所示，在柱體迎風面，隨著阻塞比的增大，除靠近邊緣區(qū)域的風壓系數均值略微減小而脈動略微增大外，內部區(qū)域的風壓系數幾乎沒有變化.在兩側面，隨著阻塞比的增大，負風壓系數的均值絕對值以及脈動值均顯著增加(見圖 7a、7b).同時可以看出，柱體兩側面，負風壓系數均值絕對值和脈動風壓系數都在靠近后緣附近出現極大值.極大值的位置對應柱體側面附近回流區(qū)渦心的位置.隨著阻塞比的增大，該極大值愈加明顯，回流區(qū)渦心附近區(qū)域的流速在增加，且渦心位置逐漸向上游移動.

a 平均風壓系數

b 脈動風壓系數

在緊鄰柱體背風面下游形成一對反向旋轉的回流渦，由此在背風面形成穩(wěn)定的負壓.隨著阻塞比的增大，負風壓系數均值絕對值和脈動值，均有較小的增加.而阻塞比R≤10%的兩個工況，洞壁對流場的約束不強烈，其平均風壓系數及脈動風壓系數的分布均比較接近.

圖 8為各工況沿計算域中軸線分布的時間平均風壓系數.隨著阻塞比的增大，尾流區(qū)的負風壓系數均值會顯著減小.在x/D=1附近風壓系數極小值位置處，當阻塞比R由4.2%增大到25%時，極小值由-2.23減小到-2.97，出口處的風壓系數均值由-0.33減小到-0.89.同時可以看出，對于阻塞比較小(R≤10%)的兩個工況，其沿軸線的平均風壓系數分布比較接近，阻塞比的影響比較小.

圖8 流場中軸線平均風壓系數分布隨阻塞比的變化規(guī)律

2.2.4阻塞比對柱體附近流場的影響

圖9a、9b所示為順流向平均速度分布，隨著橫向y的增加，順風向速度由零先減小，在y/D=0.6附近取得極小值之后迅速增加，在y/D=0.8附近增加到極大值，之后再逐漸減小，直到在邊界附近迅速減小到零(圖中未畫出).隨著阻塞比的增大，極大值的位置會逐漸向遠離柱體方向移動，且數值有明顯的增加.相比于上游(x/D=-0.125)，下游位置(x/D=0.125)的順流向平均速度的極大值出現的位置要更遠離柱體.通常定義邊界層厚度為壁面附近速度為遠場速度的99%的位置[26]，由此可見隨著阻塞比的增大，柱體壁面邊界層的厚度略有減小.

順流向雷諾正應力沿橫向分布如圖 9c、9d所示.對于上游位置處，隨著遠離柱體，雷諾正應力先增加后減小.而對于下游位置，雷諾正應力沿橫向分布的規(guī)律要更為復雜一些，要經過兩次先增加后減小的過程，最后逐漸減小到接近于零.隨著阻塞比的增大，雷諾正應力也在逐漸變大.對于上游位置處，在y/D<0.7范圍內，阻塞效應較為明顯，在y/D>0.7范圍內，不同阻塞比工況的雷諾正應力分布比較接近；對于下游位置處，這個分界點為y/D=0.8.

雷諾剪切應力沿橫向分布如圖 9e、9f所示.在上游的y/D<0.8以及下游的y/D<0.9區(qū)域內，隨著阻塞比的增大，雷諾剪切應力逐漸減??；在上游的y/D>0.8以及下游的y/D>0.9區(qū)域內，隨著阻塞比的增大，雷諾剪切應力有逐漸增加，隨著與柱體距離的增加，阻塞效應會逐漸減弱.另外當阻塞比R≤10%時，流場特性比較接近，阻塞比對流場特性的影響較小.

a x/D=-0.125b x/D=0.125

c x/D=-0.125d x/D=0.125

e x/D=-0.125f x/D=0.125

2.2.5阻塞比對流場瞬時渦量的影響

圖10給出了4種不同阻塞比(10%、16.7%、20%、25%)工況的瞬時渦量云圖.圖中所示平面為流場中軸面，瞬時渦量采用量綱為一的形式(ω*=ωD/U0).以逆時針旋轉為正，實線代表的是正向渦旋，虛線代表的是反向渦旋.圖中所有量綱為一渦量等值線最大值、最小值取值分別為±10，間隔為1.

對比4組不同阻塞比工況，可以看出隨著阻塞比的增大，尾流流場渦旋的復雜性在增強、脫落的旋渦間距減小.如圖 10a所示，當R≤10%時，洞壁對柱體繞流的約束作用較小，在柱體后緣脫落的旋渦向尾流耗散；由于洞壁給定的是無滑移邊界條件，因此在上、下邊界表面會形成一定厚度的邊界層.當阻塞比增大到R=16.7%時，洞壁對流場的約束效應已經比較明顯，柱體背風面附近的渦旋也更加復雜，尾流渦旋的耗散進一步減弱，在x/D=12處依舊可以看見比較明顯的大尺寸旋渦；值得注意的是，此時，上、下洞壁表面不僅有邊界層的出現，同時在此邊界層內部也會出現一些小的渦旋，并向下游脫落.當阻塞比進一步增大到20%和25%時，尾流渦旋的耗散進一步減弱，同時在洞壁邊界層內形成的渦旋的尺寸也在隨之增大，并與柱體脫落的旋渦相互摻混.尤其是R=25%時，在上側洞壁表面x/D=5以及x/D=13位置處有明顯的大尺寸渦旋.Kim[10]通過被動粒子追蹤模擬，也觀察到沿上、下表面釋放的粒子會形成二次旋渦，并向下游移動，這種現象會周期性地在上、下壁面交替出現.Okajima[22]在低雷諾數條件下，進行不同阻塞比的二維方柱的數值模擬也發(fā)現了相同的現象.

a 工況2(R=10%)

b 工況3(R=16.7%)

c 工況4(R=20%)

d 工況5(R=25%)

3 結論

對4.2%、10%、16.7%、20%和25%這5種不同阻塞比條件下的方柱繞流進行了CFD大渦模擬，將計算域橫向邊界設置為“wall”邊界條件，用于模擬洞壁對流場的約束作用.通過與以往風洞試驗及數值模擬結果的對比，驗證了本文計算結果的合理性.在此基礎上，分析了阻塞比對柱體氣動力、表面風壓分布、柱體附近流場特性以及尾流瞬時渦量的影響，得到了以下結論：

(1) 當阻塞比4.2%逐漸增大到25%時，方柱的平均阻力系數增大了27.6%，脈動升力系數增大了43.6%，Strouhal數增大了31%，且阻力系數頻譜曲線會在2倍Strouhal數處逐漸出現明顯的峰值.

(2) 隨著阻塞比的增大，方柱迎風面的風壓系數幾乎不變，而背風面和側面的負風壓系數均值絕對值以及脈動值均明顯增加；尾流中軸線的負風壓系數均值絕對值隨著阻塞比的增大明顯增大.

(3) 隨著阻塞比的增大，在分離剪切層內部、緊鄰柱體表面區(qū)域的氣流回流速度以及分離剪切層外部的順流向流速均會增加；阻塞比對雷諾應力的影響主要表現在分離剪切層內部區(qū)域，并且下游受到的影響明顯要強于上游位置處.

(4) 當阻塞比增大到一定程度(本文為R=16.7%)后，洞壁對流場約束作用的增強，在尾流區(qū)洞壁上、下表面附近會出現二次旋渦；隨著阻塞比的增大，該渦旋的尺寸增加，與方柱上脫落的渦旋摻混現象也會越加明顯.

本文重點考察當阻塞比在較大范圍內變化時，柱體氣動力及流場特性的變化.為了明確給出允許范圍內的有價值的阻塞度建議值以及具有普適性的修正方法，需要更多的計算模擬，同時必須進行大量的模型試驗，進行更深入地研究.