串列雙錐柱繞流流動特性的大渦模擬研究1)

劉 健 鄒 琳 陶 凡 左紅成 徐漢斌

(武漢理工大學機電工程學院,武漢 430070)

引言

風能作為一種重要的可再生清潔能源,對其開發利用受到國內外學者的廣泛關注[1-2].對于風能的獲取,如渦激振動發電[3]、風力壓電效應發電[4],增振就尤其重要.

圓柱作為繞流的經典研究對象,通常通過改變圓柱幾何形狀來改變周圍流動狀態以實現流動控制,Zhang 等[5]通過改變波浪外形圓柱的波長比以及波幅比,發現波幅比為0.2 時波浪圓柱受到的時均阻力系數顯著降低.Zhang 等[6]通過對比研究無限長和一端自由的圓柱,發現有限長圓柱的時均阻力系數和脈動升力系數遠小于無限長圓柱.趙桂欣[7]通過對比四種不同端面以及側面的兩端自由的有限長圓柱,發現改變外形對柱體阻力系數均有降低作用,在特定雷諾數下,脈動升力系數較單直圓柱有所提升.趙萌等[8]、喬永亮[9]的研究發現有限長圓柱長徑比會明顯改變圓柱的流場結構,存在一個臨界長徑比使得有限長圓柱的流動模式由自由端來流主導的無序渦脫模式轉向側面來流主導的卡門渦街渦脫模式.楊耀宗[10]的研究發現,有限長錐柱在斜率k=0.05 時錐柱振動性能有顯著提升.由此可見,對于有限長圓柱,通過改變其表面形狀對脈動升力系數有較大提升.

雙圓柱較單圓柱的流動模式更為復雜.Alam等[11-12]通過實驗,研究不同布置形式的雙圓柱,對兩個圓柱之間的尾流形態進行分類,研究發現在間距比為2.4～ 3.0,交錯角為10°時的錯列雙圓錐脈動阻力出現最大值.Zdravkovich[13]研究兩個串列圓柱間距比小于5 的情況,發現間距比小于3 時,上、下游圓柱之間不會產生明顯渦脫,在間距比大于3 時,上下游圓柱之間逐漸產生渦脫.杜曉慶等[14-16]對不同布置狀態下的雙圓柱進行了數值仿真,按照上、下游圓柱相互作用的不同形式將流態進行分類,解釋了下游圓柱在上游圓柱尾流的不同狀態下的升阻力特性,發現下游圓柱被上游包裹狀態下,下游圓柱的阻力系數顯著降低,脈動升力系數反而受到一定抑制,但是相對于單圓柱,下游圓柱脈動升力系數均有不同程度增幅.文獻[17-19]的研究發現,串列雙圓柱隨間距比增加,脈動升力系數和時均阻力系數先是緩慢降低,但是在間距比為3.0～ 4.0 附近發生突增跳躍到一個較高的水平,上游圓柱對下游圓柱的影響由剪切層包裹的抑制狀態逐漸變為渦脫撞擊狀態,使得下游圓柱在超過臨界間距比后脈動升力相對于單圓柱大幅提升.對于改變圓柱外形的串列雙圓柱,唐濤等[20]數值模擬了間距比為3 的不同長度比的梯形截面串列柱,發現兩圓柱之間的間隙流會隨著梯形長度比的改變而呈現不同的耦合模式,導致升力發生改變.因此上游圓柱的尾流結構對下游圓柱脈動升力的提升有著決定性的作用,但是對于有限長串列雙錐柱,其流動模式和升阻力特性缺乏研究.

綜上所述,下游柱體脈動升力系數有可能相對于單柱體有大幅提升,同時表面的形狀改變對于柱體脈動升力也有很大影響.根據課題組前期研究,發現錐形圓柱在斜率k=0.05 時相對于直圓柱增振最明顯[10],因此,本文將以提升脈動升力系數為目的,利用商用軟件Fluent 中的大渦模擬,在亞臨界雷諾數(Re=3900)下,對斜率k=0.05 的串列雙錐柱在間距比L/Dm=2～ 10 下進行數值模擬,通過分析上、下游錐柱之間的流動機理,探索下游錐柱時均阻力系數和脈動升力系數變化規律,為風力俘能結構群的列陣布局提供理論支持.

1 數值方法

1.1 控制方程

大渦模擬(large eddy simulation,LES)[21-22]采用空間濾波技術,通過截止尺度和濾波函數將渦識別為大渦和小渦,大渦直接解析,小渦則通過亞格子尺度應力模型(sub-grid-scale stress,SGS)封閉求解,過濾后的不可壓縮黏性流體Navier-Stokes 方程為

式中,xi,yi為笛卡兒坐標,i,j=x,y,z,ui,uj為濾波速度矢量,為流體所受壓力,ρ為流體密度,ν為流體運動黏度,τij為亞格力應力.

根據Smagorinsky[23]的基礎SGS 模型,設應力表達式

式中,Δi為網格在i方向的尺度,Cs為Smagorinsky尺度,本文取0.2[24].

1.2 控制方程離散

本文計算控制方程離散均采用有限體積法,壓力速度耦合方式為PISO,離散化均采用二階格式[6].

2 計算域及網格

2.1 模型簡介

基于多數人研究的直圓柱,本文加入斜率k,將直圓柱轉化為錐柱,錐柱如圖1 所示,幾何表達式如下

圖1 錐柱結構示意圖及表面網格劃分Fig.1 Schematic diagram of conical cylinder structure and meshing

式中,Dz為錐柱在高度為z處的橫截面直徑,Dm=0.01 m 為錐柱的平均直徑,k為錐柱的斜率,圓柱高度比H/Dm=7,Dmax和Dmin分別為錐柱最大直徑和最小直徑.

本文中出現的相關參數定義如表1.

表1 相關參數定義Table 1 Parameter definition

2.2 計算域及網格劃分

本文研究的串列雙錐柱計算域大小為(30Dm+L) × 20Dm× 10Dm,X軸正方向為順流向,Y軸為橫流向,Z軸為圓柱展向,笛卡爾坐標系圓心位于上游錐柱中心處.上游錐柱中心距離計算域入口10Dm,下游錐柱中心距離計算域出口20Dm,錐柱的上下端面距離計算域頂部和底部均為1.5Dm,如圖2 所示.入口采用速度入口(velocity-inlet),U∞=5.772 m/s,出口采用壓力出口(pressure-outlet),背壓設置為0 Pa,計算域兩個側面以及頂面采用對稱邊界(symmetry),計算域底面和圓柱體表面為無滑移壁面(noslip-wall).流體介質為空氣,密度ρ為1.225 kg/m3,運動黏性系數ν為1.48 × 10-5m2/s.

圖2 串列雙錐柱計算域及網格示意圖Fig.2 Schematic diagram of calculation domain and grid of tandem two conical cyliners

本文使用ICEM 對計算域進行結構化網格劃分,利用O-block 技術對圓柱周圍以及頂部網格進行加密,同時避免了片狀低質量網格的生成.大渦模擬對于邊界層的要求控制y+≤1,每個圓柱周圍有3Dm×3Dm指數形式增長的加密區域,同時也對兩個圓柱之間以及下游錐柱尾流區的網格進行了加密,如圖2 所示.

2.3 計算模型驗證

本文首先驗證使用單直圓柱計算模型,探究時間步長Δt、網格參數(圓周節點數,第一層邊界層高度)對圓柱的時均阻力系數(Cdmean),脈動升力系數(Clrms)和Strouhal 數(St)的影響,計算結果如表2.綜合考慮計算成本和計算精度的影響,本文選用Case2 的網格精度和時間步長進行下面仿真.

表2 時間步長、圓周節點數和第一層高度對時均阻力系數,脈動升力系數和Strouhal 數的影響Table 2 The influence of time step and grid parameters on the Cdmean,Clrms and St

由于自由端數目以及長徑比不同,本文單圓柱數值結果與仿真數據(num.)[25]以及實驗數據(exp.)[26]有一定誤差,為了進一步驗證本次仿真所用的算法,網格精度以及時間步長能夠適用于串列雙錐柱,對串列雙錐柱的流動結構進行煙線實驗,其實驗裝置如圖3 所示.風洞為直流式風洞(CSUWT2450-40),煙線控制器(SW-2) 通過充電放電加熱涂有石蠟油的銅絲,石蠟油在瞬間高溫下產生白色煙霧,風洞試驗段均勻來流經過銅絲時,白煙在風的作用下繞過固定在風洞低端串列布置的錐柱,通過相機(Canon EOS 850 D)捕捉串列雙波浪錐柱之間白煙形成的流場.對串列雙波浪錐柱四個截面Y=0,S1(z/H=3/4),S2(z=0),S3(z/H=-3/4)與同截面下仿真結果進行對比.

圖3 串列雙錐柱實驗臺架Fig.3 Experimental of two conical cylinders in tandem arrangement

實驗結果如圖4 所示,串列雙錐柱實現顯示,Y=0 截面可以看到,間距比L/Dm=5 的串列雙錐柱兩個自由端來流撞擊在下游波浪錐柱迎風面,這與仿真的結論一致,觀察S1,S2,S3 截面發現,實驗顯示錐柱后方回流區發展完成,剪切層已經收縮,尾流撞擊在下游錐柱表面,仿真結果與實驗結論基本一致.因此可以認為本文驗證中Case2 的算法,網格精度及時間步長可用于串列雙錐柱仿真研究.

圖4 間距比L/Dm=5 時串列雙錐形圓柱煙線實驗圖與仿真流線圖Fig.4 Smokeline diagrams and numerical streamline diagrams of two conical cylinders in tandem arrangement at spacing ratio L/Dm=5

3 結果分析

3.1 時均阻力系數和脈動升力系數

時均阻力系數和脈動升力系數是風力俘能結構獲取風能的兩個重要參數.本文與文獻[18]上游柱和下游柱及Case2 中單圓柱數據對比結果,如圖5所示,其中“UC”表示上游柱,“DC”表示下游柱.

圖5 串列雙錐柱時均阻力系數和脈動升力系數隨間距比變化Fig.5 Cdmean and Clrms of two conical cylinders in tandem arrangement change with the spacing ratios

本文上、下游錐柱與文獻中上、下游圓柱的Cdmean和Clrms隨間距比變化規律基本一致,但是數值大小存在很大差異,可能是因為本文錐柱自由端的影響,下文將通過流場結構分析其原因.本研究發現,Cdmean和Clrms在L/Dm=4～ 5 時存在突變,在這個區間可能發生流動狀態的改變;在臨界間距比之前上、下游錐柱Clrms呈下降趨勢,而Cdmean則隨間距比先增加后減小;在間距比超過臨界值之后,上游錐柱Clrms在臨界值之后基本保持不變,且與單直圓柱基本一致.下游錐柱Clrms在間距比L/Dm=5～ 6有一個略微降低然后突增的過程,在L/Dm=7 的時候出現峰值,此時較單直圓柱增長約20.7 倍,此階段可能存在流動結構改變.上游錐柱Cdmean在低于單直圓柱一定值處波動,下游錐柱Cdmean由一個較小值逐漸增大,趨近于上游錐柱,下游錐柱相對于單直圓柱Cdmean最大減小約90.7%,在L/Dm=7 時,減小約19.8%.

3.2 表面壓力系數

表面壓力系數Cp能夠直接反映出錐柱表面受力情況,為了方便觀察表面Cp分布情況,如圖6(a)所示,將錐柱的表面沿母線剪開,并以角度“θ”和高度“Z/H”作為橫、縱坐標鋪平,表面時均Cp分布如圖6(b)所示.

圖6 間距比L/Dm=2～ 9 下串列雙錐柱時均壓力系數分布圖Fig.6 Distribution diagram of time-average Cp of two conical cylinders in tandem arrangement with spacing ratio L/Dm=2～ 9

從整體來看,無論上游錐柱還是下游錐柱在不同間距比下的時均Cp分布基本關于180°對稱,說明來流對上、下游錐柱兩側的影響在時間平均上基本相同.根據云圖分布的不同,可將上游圓柱的Cp分為兩類,第一類為低壓區(壓力系數小于-0.6 的區域)集中在錐柱頂部(DIS1),如L/Dm=2,3,5,第二類分為高壓區(壓力系數大于0.6 的區域)主要集中在錐柱底部(DIS2),如L/Dm=4,6,7,8,9.可以看到上游錐柱除L/Dm=2,3,4,低壓區強度相對較大以外,其他間距比分布基本相近,這是Cdmean在臨界區域前端逐漸降低的主要原因;由于下游錐柱受到上游錐柱尾流的影響,在小間距比下壓力分布強度較弱,隨著間距比增加,高壓區和低壓區范圍和強度都有所增加,并且分布形式逐漸接近上游圓錐,這也是下游錐柱阻力系數隨間距比增加而增大的主要原因(臨界區域除外).錐柱的表面壓力分布主要是由其流場決定,下面從流場的本質對壓力分布進行分析.

3.3 時均流線圖

三維流場較為復雜,但串列雙錐柱的流場分布具有一定對稱性,下面取Y/Dm=0,S1,S2,S3 截面(見圖3(b))對串列雙錐柱的時均流線進行分析,如圖7 所示.

圖7 間距比L/Dm=2～ 9 下串列雙錐柱在Y=0,S1,S2,S3 截面時均流線圖Fig.7 Time-average streamline diagram of two conical cylinders in tandem arrangement with spacing ratio L/Dm=2～ 9 at Y=0,S1,S2,S3 sections

由于自由端的存在,錐柱上、下端部流體分別存在上洗(流經下自由端的空氣向上運動)和下洗(流經上自由端的空氣向下運動)現象,上洗、下洗作用會抑制側面來流在錐柱后方形成周期性渦脫落,從而導致脈動升力系數降低,這也是串列雙錐柱脈動升力遠小于串列無限長圓柱的主要原因[27-28].空氣流經上游錐柱在其后方形成一對Y截面的展向回流區(arch-vortex),同時空氣流過柱體側面,與arch-vortex 相互作用,在柱體后方形成Z截面的橫向回流區(side-vortex,以下簡稱回流區),除此之外,在柱體端面處能觀察到頂渦(tip-vortex)[29-30],如圖7(g).課題組前期研究發現[10],由于錐柱斜率的存在,導致自由端上洗作用強于下洗,下端形成的arch-vortex 擠壓上端形成的arch-vortex 使其變為橢圓形,對于串列雙錐柱,在L/Dm≥ 7 時,上游錐柱后方能夠明顯觀察到一對arch-vortex,但在間距比L/Dm=2～ 6 時,受到下游錐柱影響,上游錐柱自由端形成的渦結構發生明顯變化,根據上洗和下洗在兩個圓柱之間形成的arch-vortex 不同,將兩錐之間的流場形式分為兩類,第一類是由下洗作用占主導,流體經上自由端向下運動一部分錐柱上端形成一個渦,另一部分與下自由端的上洗流相作用在錐柱下端形成一個渦,如L/Dm=2,3,5,第二類則與第一類相反,如L/Dm=4,6,7,8,9.這兩種分類方式正好也解釋上一節中DIS1 和DIS2 的分布規律,占主導作用的一端,空氣流速更快,上游錐柱背風面對應位置壓力也就越低.而無限長串列圓柱上游柱背風面不存在自由端影響,背風面形成的回流區主要來自側面,壓力系數展向對稱分布[6].串列雙錐柱間流動結構隨間距比變化可分為三種狀態:剪切層包裹狀態,過渡狀態,尾流撞擊狀態.剪切層包裹狀態,上游錐柱回流區受下游錐柱影響并未完全發展,剪切層完全包裹住下游錐柱,如圖7 中L/Dm=2,3;尾流撞擊狀態,上游錐柱當前位置回流區已經完全發展并且大小隨間距比基本不發生改變,且下游錐柱已位于回流區之外,不再被上游剪切層包裹,如圖7 中間距比L/Dm=7,8,9;過渡狀態,即上游錐柱自由端來流作用到下游錐柱表面,而側面剪切層部分包裹下游錐柱,如圖7 中間距比L/Dm=4,5,6.隨著三種狀態變化,上游錐柱的上洗和下洗直接作用在下游錐柱表面的范圍逐漸減小,同時下游錐柱后方的回流區隨著間距比增加而逐漸擴大,這也是下游錐柱隨間距比增加,其高壓區和低壓區范圍逐漸擴大并增強的主要原因.

3.4 渦結構圖

渦的識別方法可以分為三代[31],其中第二代中的Q方法運用較為廣泛,表達式如下

式中,A為對稱張量,B為反對稱張量,為矩陣的Frobenius 范數.下面對串列雙錐柱的渦采用Q方法識別,如圖8 所示,每個間距比提供兩個視角view1,view2.

圖8 間距比L/Dm=2～ 9 下串列雙錐柱瞬時渦量圖(Q=1 × 104),view1 為整體視圖,view2 為俯視圖Fig.8 Instantaneous vorticity diagram of two conical cylinders in tandem arrangement with spacing ratio L/Dm=2～ 9 (Q=1 × 104).View1 is the general view;view2 is the top view

從渦量圖中可以看到,前面提到的兩種分類,DIS1 在這里顯示的是大量渦“下泄”集中在兩個錐柱之間的中下部,如圖8(a)(b)(d),DIS2 在這里反應的是大量渦“上瀉”集中在兩個錐柱之間的中上部,如圖8(c)～圖8(h).從渦的結構上來看,間距比較小時,上游錐柱脫落的渦碎且無序,與下游圓柱產生的渦相互作用,在下游圓柱后方形成“肋狀渦”(rib vortex),如圖8(a)中所示,當間距比逐漸增大時,上游錐柱脫落的渦逐漸變大并且有周期特性,下游錐柱后方“肋狀渦”也逐漸變大變長,如圖8(e)～圖8(h)所示.而對于無限長圓柱,不同間距比下均能觀察到卡門渦街現象[19],所以其脈動升力系數能明顯高于本文有限長雙串列錐柱.上游錐柱脫落的尾渦作用在下游錐柱表面,與來自下游錐柱端部以及側面來流相互作用,產生較大的升力,相較于低間距比串列雙錐柱,當間距比為L/Dm=7～ 9 時,下游錐柱處于上游錐柱回流區臨界位置或回流區下游,上游錐柱尾流能夠充分發展,形成類似卡門渦街特征,周期性的漩渦作用在下游錐柱表面,使下游圓錐表面更容易產生大的升力.間距比L/Dm=7 時,上游錐柱尾流發展最充分,其尾流對下游錐柱的作用更為強烈,這可能是下游錐柱Clrms相對于間距比L/Dm=8,9 更大的原因.下面通過分析幾個截面的渦量,來具體分析上下游錐柱之間的流動情況.

3.5 Z 截面渦量圖

為了方便觀察串列雙錐柱之間的流動形態,下面取S1,S2,S3(見圖3(b)) 三個截面的渦量圖,如圖9 所示.尾流之間的干擾可分為很多不同形式,根據錐柱截面直徑以及流場特性,參照文獻[12-13]對流場的分析,將Z截面渦量分為兩類,V1:上游錐柱尾流包裹下游錐柱,如圖9 中雙點劃紅線外框中渦量圖,上、下游錐柱共同形成的剪切層極不穩定,在下游錐柱后方形成的回流區相較于單錐柱[8]很小,在后方形成非周期性的渦脫;V2:下游錐柱位于上游錐柱回流區外端或恰好位于其回流區邊界處,如圖7 中虛線黑線外框中渦量圖,上游錐柱尾渦得到充分發展,尾流撞擊在下游錐柱表面,結合下游錐柱產生的渦在下游錐柱后方產生大量不規則碎渦,相較于V1 為下游錐柱提供更大的升力.結合圖8 中渦量圖顯示,上游錐柱的上洗和下洗作用強度不同,使得作用在下游錐柱前端的渦集中位置不同,如圖8(f),上游錐柱尾渦集中在錐柱中下部,導致圖9 中S1 截面兩個圓柱相互作用較弱,這也是圖6 中間距比L/Dm=4 的下游圓柱頂部時均壓力系數高的主要原因.圖9 中,當間距比L/Dm=2,3 時,三個截面流動模式均為V1,即4.3 節提到的剪切層包裹狀態,間距比L/Dm=4,5,6 中,三個截面中V1,V2 模式共同存在,即過渡狀態,而當間距比L/Dm=7,8,9 時,三個截面流動模式全部變為V2,即尾流撞擊狀態,這是間距比L/Dm=7,8,9 有更高的脈動升力值主要原因.從圖9 中可以觀察到,圖9(f) 相對于9(g) 和9(h)每個截面下游錐柱距離上游錐柱回流區更近,上游錐柱發展的周期性尾渦作用在下游錐柱表面更劇烈,這是下游錐柱在L/Dm=7 取得Clrms幅值的主要原因.

圖9 間距比L/Dm=2～ 9 下串列雙錐柱在S1,S2,S3 截面的渦量圖Fig.9 Z-vorticity of two conical cylinders in tandem arrangement with spacing ratio L/Dm=2～ 9 at S1,S2,S3 sections

4 結論

本文在Re=3900 下利用大渦模擬對間距比為L/Dm=2～ 10 的串列雙錐柱進行了固定繞流仿真.首先,驗證了Case2 時間步長以及網格達到收斂,通過煙線實驗證明本次網格精度,時間步長和算法可用于串列雙錐柱仿真;其次,分析了其升阻力系數隨間距比變化的特性;再次,分析了其時均流場信息分析其阻力變化的主要原因;最后,通過渦量圖分析了脈動升力系數變化的主要原因,得出以下結論.

(1)隨間距比增加,上游錐柱Clrms先減小逐漸趨近于單直圓柱,在L/Dm=4 之后基本保持不變;上游錐柱Cdmean先減小,在L/Dm=4 之后保持低于單直圓柱的波動狀態.

(2)下游錐柱Clrms隨間距比基本保持增加后減小趨勢(除臨界區間),在L/Dm=7 時,取得峰值,約為單直圓柱的21.7 倍;下游錐柱Cdmean由一個較小值隨間距比增加逐漸接近上游錐柱的Cdmean,下游錐柱相對于單直圓柱Cdmean最大減小約90.7%,在L/Dm=7 時,減小約19.8%.

(3)相較于無限長串列雙圓柱,串列雙錐柱的上游錐柱自由端的上洗和下洗作用能顯著影響側面來流在錐柱后方形成的渦脫的發展,抑制其產生規則卡門渦街現象;上洗和下洗作用強的一端會使上游錐柱背風面形成大且強的低壓區,同時使同端下游錐柱迎風面的高壓區減小.

(4)隨間距比增加,下游錐柱與上游錐柱回流區的關系有剪切層包裹狀態、過渡狀態、尾流撞擊狀態,其中尾流撞擊狀態時,上游錐柱尾流能夠充分發展,形成交替脫落的肋狀渦,撞擊在下游錐柱表面,使下游錐柱產生較大的脈動升力系數,間距繼續增加,這種相互作用的強度會隨之減弱,脈動升力系數隨之減小.本文的研究能為風力俘能結構的列陣布局提供理論支持.