三維錯列雙波浪錐柱繞流流動特性數(shù)值仿真

鄒 琳，劉 健，吳偉男，閆豫龍，魏翼鷹

(武漢理工大學(xué) 機(jī)電工程學(xué)院，武漢 430000)

圓柱繞流的研究受到國內(nèi)外學(xué)者的廣泛關(guān)注，大部分學(xué)者研究目的是減阻減振，而對于某些工程領(lǐng)域，如渦激振動發(fā)電[1]，減阻增振研究就顯得尤為重要。

對于單圓柱研究，國內(nèi)外許多學(xué)者圍繞改變圓柱外形來控制圓柱周圍的流動結(jié)構(gòu)，來達(dá)到減阻或增振的目的。LAM等[2]和ZHANG Zhihao等[3]通過改變波浪形圓柱波長比和波幅比，發(fā)現(xiàn)波浪形圓柱的阻力系數(shù)相對于直圓柱能顯著降低。而存在自由端的圓柱，ZHANG Hui等[4]發(fā)現(xiàn)有限長圓柱時均阻力系數(shù)以及脈動升力系數(shù)遠(yuǎn)小于無限長圓柱。趙桂欣等[5]對12種不同波長的單自由端圓柱進(jìn)行了研究，發(fā)現(xiàn)對于有限長圓柱，波浪外形也能有很好的減阻效果，最大可達(dá)5.36%。楊耀宗[6]則通過在波浪形圓柱中引入錐度得到波浪錐形圓柱，發(fā)現(xiàn)波長比為1.75，波幅比為0.1，斜率為0.05的波浪錐形圓柱脈動升力系數(shù)相對于直圓柱增加約20%。

雙圓柱之間的流動模式十分復(fù)雜，ALAM等[7-8]通過實(shí)驗(yàn)觀察不同交錯角和間距比下的流動特性，將兩個圓柱之間的流動狀態(tài)進(jìn)行分類，同時發(fā)現(xiàn)在交錯角為10°時，間隙流作用于下游柱最強(qiáng)烈，下游柱脈動阻力達(dá)到最大。SUMNER等[9]通過實(shí)驗(yàn)重點(diǎn)研究了錯列圓柱近尾渦的流動模式，發(fā)現(xiàn)下游柱阻力系數(shù)較單圓柱低。LAM等[10]研究發(fā)現(xiàn)，間距比為3.5左右為串列波浪錐形圓柱臨界間距比，其他間距比減阻效果仍然存在。杜曉慶等[11]通過數(shù)值模擬研究錯列雙圓柱的脈動升力系數(shù)，發(fā)現(xiàn)在較大間距比及較小的交錯角下，下游柱的脈動升力系數(shù)較單圓柱能夠顯著提升。汪洋等[12]研究的錯列雙圓柱，發(fā)現(xiàn)在交錯角約為8°時，下游柱的脈動升力系數(shù)達(dá)到最大值。

綜上所述，在小交錯角以及中等間距比時，下游柱有很好的減阻增振效果。因此，本文采用商用軟件Fluent中的大渦模擬，在亞臨界雷諾數(shù)(Re=3 900)下，探究錯列波浪錐柱在間距比L/Dm=4、5下(L為兩個波浪錐柱的中心距離，Dm為波浪錐柱的平均直徑)，隨交錯角α(α=0～15°)變化的上下游柱之間流場機(jī)理，分析其升阻力變化原因，以達(dá)到減阻增振目的，為風(fēng)力俘能結(jié)構(gòu)列陣布局提供理論支持。

1 控制方程

本文采用大渦模擬進(jìn)行求解，大渦模擬通過濾波函數(shù)和截止尺度Δ將渦分為大尺度旋渦和小尺度旋渦兩種，大尺度旋渦通過Navier-Stokes方程直接求解，小尺度渦則通過SGS(Sub-gird-scale Stress,SGS)模型進(jìn)行求解[11]。

對于不可壓縮粘性流體，其Navier-Stokes運(yùn)動方程

(1)

(2)

亞格子應(yīng)力τij計(jì)算采用Smagorinsky-Lilly的亞格子模型：

(3)

(4)

(5)

(6)

(7)

式中：Δ表示過濾的網(wǎng)格尺度，Cs為Smagorinsky常數(shù)，本文取0.2。

2 模型驗(yàn)證

計(jì)算域大小為30D×10D×10D，D為圓柱直徑，D=0.01 m，圓柱總長為H0，長徑比H0/D=7，笛卡爾坐標(biāo)系的原點(diǎn)位于圓柱中心，圓柱中心距入口為10D，距離出口為20D，x軸正方向?yàn)轫樍飨?，z軸正方向沿圓柱軸心向上，如圖1所示。

(a)整體視角 (b)Z=0截面視角圖1 網(wǎng)格驗(yàn)證及計(jì)算域示意圖Fig.1 Grid model and calculation domain

用結(jié)構(gòu)化網(wǎng)格對整個流場進(jìn)行網(wǎng)格劃分，其中圓柱周圍有0.1D的O形塊區(qū)域來控制近壁面網(wǎng)格的高度，本次驗(yàn)證采用近壁面網(wǎng)格高度為0.001D，使圓柱壁面y+最大值小于1。流體為空氣，密度ρ=1.225 kg/m3，動力粘度μ=1.813 Pa·s，入口采用速度入口(Velocity-inlet)，大小為5.772 m/s，出口采用壓力出口(Pressure-outlet)，大小設(shè)置為0，側(cè)面以及底面采用對稱邊界(Symmetry)，圓柱表面采用無滑移壁面(No slip wall)。無量綱時間步長Δt*=U∞Δt/D，U∞為來流速度，Δt為時間步長,大渦模擬在計(jì)算中，要求時間步長Δt滿足庫朗數(shù)C(Courant number)在0.5～1.0之間，C=uΔt/Δx，Δx為最小網(wǎng)格尺寸。當(dāng)取Δt*=0.001時，相應(yīng)庫朗數(shù)約為0.981，滿足計(jì)算條件。

本文計(jì)算壓力速度耦合采用PISO算法，離散化采取二階迎風(fēng)格式[4]。表1為計(jì)算模型網(wǎng)格獨(dú)立性及時間步長無關(guān)性驗(yàn)證結(jié)果，本文計(jì)算結(jié)果與實(shí)驗(yàn)以及數(shù)值仿真結(jié)果基本一致，可以認(rèn)為case2網(wǎng)格密度、無量綱時間步長Δt*=0.001以及本次仿真所用的算法正確，可以運(yùn)用于錯列雙波浪錐柱仿真實(shí)驗(yàn)。

表1 仿真結(jié)果與已有文獻(xiàn)對比Tab.1 Simulation results compared with previous literature

3 計(jì)算模型

波浪錐柱如圖2所示，其幾何表達(dá)式為

(8)

式中：Dz表示波浪錐柱對應(yīng)高度Z處的直徑，Dm表示波浪錐柱的平均直徑，DH和DL分別表示波浪錐柱對應(yīng)斜率為k的直錐柱最大直徑和最小直徑，H為波浪錐柱總長，H=7Dm；a、λ分別表示波浪錐柱表面的波幅與波長。根據(jù)楊耀宗[6]的研究，選取振動性能最佳的一組參數(shù)，即k=0.05，a/Dm=0.1，λ/Dm=1.75。

本次仿真計(jì)算域?yàn)?30Dm+Lcosα)×(20Dm+Lsinα)×10Dm。原點(diǎn)位于上游波浪錐柱的中心，兩個波浪錐柱中心距為L，中心線與入流向夾角為α，即交錯角，波浪錐柱長徑比H/Dm=7。

圖2 波浪錐柱幾何模型示意圖Fig.2 Geometric model of wavy conical cylinder

入口采用速度入口(Velocity-inlet)，大小為5.772 m/s，出口采用壓力出口(Pressure-outlet)，背壓設(shè)置為零，波浪錐柱壁面采用無滑移邊界條件，計(jì)算域頂面、底面和側(cè)面均采用Symmetry壁面，如圖3所示。為了探究小交錯角下和中等間距下波浪錐柱的升、阻力系數(shù)特性，本次研究的間距比L/Dm為4和5，交錯角α為0°、5°、7.5°、10°、12.5°和15°。網(wǎng)格劃分方法和數(shù)值算法均參照驗(yàn)證中的Case2，網(wǎng)格模型如圖3所示。

(a)整體視角 (b)Z=0截面圖3 錯列波浪錐柱計(jì)算域及網(wǎng)格模型示意圖Fig.3 Schematic diagram of calculation domain and grid model of two staggered wavy conical cylinders

4 結(jié)果分析

4.1 升阻力系數(shù)

升力系數(shù)Cd和阻力系數(shù)Cl的計(jì)算公式如下:

Cd=2Fd/(ρU∞2Dm·H)

(9)

Cl=2Fl/(ρU∞2Dm·H)

(10)

式中Fd、Fl分別為順流向受到的阻力和橫流向受到的升力。

圖4 (a)為上游波浪錐柱和下游波浪錐柱(下文分別簡稱為上游柱和下游柱)的時均阻力系數(shù)隨交錯角變化曲線。當(dāng)交錯角為α=0°時，下游柱時均阻力系數(shù)遠(yuǎn)小于上游柱，隨交錯角增加，下游柱時均阻力系數(shù)急劇上升，逐漸逼近上游柱時均阻力系數(shù)；除間距比L/Dm=4，交錯角α=15°工況外，兩種間距比的上游柱阻力系數(shù)均略小于單直圓柱；對于下游柱，間距比L/Dm=5的時均阻力系數(shù)大于間距比L/Dm=4工況，最大可達(dá)24.7%。圖4(b)是升力系數(shù)的均方根值Clrms曲線。兩種工況下下游柱的Clrms均遠(yuǎn)大于上游柱，且隨交錯角的變化趨勢基本相同，均在交錯角α=10°時取得最大值，間距比L/Dm=4和L/Dm=5取得脈動升力系數(shù)最大值時相對單直圓柱分別提高約20.1倍和21.4倍。在α=10°時脈動升力系數(shù)取得最大值，這與ALAM等[8]和汪洋等[12]的研究基本一致。時均阻力系數(shù)的降低和脈動升力系數(shù)的顯著提高將有利于風(fēng)力俘能結(jié)構(gòu)的能量獲取[6]。

(a)時均阻力系數(shù)

4.2 波浪錐柱表面壓力系數(shù)

圓柱表面的時均壓力系數(shù)能夠清晰展示表面的壓力分布規(guī)律，為方便觀察其分布形式，將圓柱表面以角度展開成二維平面，如圖5(a)所示，展開后的壓力系數(shù)分布如圖5(b)所示。

時均壓力系數(shù)計(jì)算公式為

(11)

對于上游柱，其壓力系數(shù)分布隨間距比和交錯角變化，基本保持對稱，高壓區(qū)(>0.8)約在在0～25°和335～360°之間，由于波浪錐柱形狀的影響，沿z方向變化時區(qū)域有一定波動，低壓區(qū)中心(<0.6)約在75°和285°，受到波浪錐柱錐度和表面的影響，其分布形式類似三角形，同時也沿z方向區(qū)域存在波動，兩種間距比下高壓區(qū)和低壓區(qū)分布范圍基本相同，這可能是阻力系數(shù)隨交錯角變化時，兩種間距比相差不大的原因。對于下游柱，其時均壓力系數(shù)隨交錯角α的變化存在較大變化，當(dāng)α=0°時，兩種間距比下，圓柱表面正壓和負(fù)壓區(qū)域能保持對稱，但高壓范圍和低壓范圍都很小，這可能導(dǎo)致下游柱時均阻力系數(shù)較?。浑S著α的增加，高壓區(qū)和低壓區(qū)范圍逐漸擴(kuò)大，當(dāng)達(dá)到α=15°時，其分布形式已經(jīng)接近上游柱，這是時均阻力系數(shù)隨著交錯角增加而增加最終接近上游柱時均阻力系數(shù)的主要原因；同時可以觀察到當(dāng)間距比為L/Dm=5時，正壓區(qū)間和負(fù)壓區(qū)間面積略大于間距比L/Dm=4，這可能是導(dǎo)致時均阻力系數(shù)隨交錯角α變化時，L/Dm=5略大于L/Dm=4的原因。

圖5 波浪錐柱表面時均壓力系數(shù)分布Fig.5 Time-average pressure coefficient on the surface of the wavy conical cylinder

4.3 流動結(jié)構(gòu)

如圖6給出了間距比L/Dm=4和5，交錯角α=0°，5°，10°和15°的三維時均流線圖和對應(yīng)的Node1、Saddle1和Node2截面的時均流線圖。

兩種間距比下錯列雙波浪錐柱的流動模式，按照下洗(流經(jīng)上自由端部的流體在柱體后方向下運(yùn)動)和上洗(流經(jīng)下自由端部的流體在柱體后方向上運(yùn)動)在兩個圓柱中間形成的流動狀態(tài)可分為兩類：

第一類為上游柱上自由端的下洗作用占主導(dǎo)，一部分流線在上游柱的后上方形成X-Z方向的回流區(qū)，一部分則繼續(xù)向波浪錐柱下端流動在下游柱迎風(fēng)面的底部與上游柱下自由端的射流相互作用形成一個X-Z方向回流區(qū)，這兩個回流區(qū)中心向兩側(cè)發(fā)展出的流體與上游柱側(cè)面來流相互作用形成“拱形”渦。如圖6(a)，(g)所示，三維圖中的“拱形”渦對應(yīng)右側(cè)Node和Saddle截面波浪錐柱之間形成的回流區(qū)。

(a)L/Dm=4,α=0° (b)L/Dm=5,α=0°

第二類為上游柱下自由端的上洗作用占主導(dǎo)，一部分流線在上游柱的背風(fēng)面底部形成Y方向的回流區(qū)，一部分則繼續(xù)向上端發(fā)展與來自上游柱頂部自由端的射流相互作用在下游柱上端迎風(fēng)面形成一個X-Z方向回流區(qū)，同樣，兩個回流區(qū)中心也會向兩側(cè)發(fā)展出流體與來自上游柱側(cè)面的流體相互作用形成“拱形”渦，如圖6(b)，(c)，(d)，(e)，(f)和(h)所示。第一類流動中，上游柱上側(cè)形成的回流區(qū)和下游柱的迎風(fēng)面下側(cè)的回流區(qū)使這兩個區(qū)域壓力降低，這與圖5(b)中L/Dm=4，α=0°，15°壓力系數(shù)分布區(qū)域相符。同樣的，第二類流動中，上游柱的背風(fēng)面下端和下游柱迎風(fēng)面的上端能形成的兩個回流區(qū)使這兩個區(qū)域的壓力降低，這可能是圖5(b)中對應(yīng)區(qū)域壓力系數(shù)低的原因。

4.4 渦量圖

4.4.1 三維渦量圖

上游柱的尾流對下游柱表面及其尾渦有很大的影響，使得下游柱相對于上游柱有更大的脈動升力系數(shù)值，并且在不同間距比和交錯角下尾流作用形式各有不同，下面通過分析不同間距比下取得脈動升力系數(shù)峰值的兩組參數(shù)(L/Dm=4，α=10°，L/Dm=5，α=10°)的渦量圖，來分析下游柱取得脈動升力最大值的原因，Q定義[16]為

式中：Ω表示渦量，S表示應(yīng)變率張量，‖ ‖為取范數(shù)運(yùn)算符。

圖7中顯示的是Q=1×106時錯列波浪錐柱周圍的渦量圖，可以觀察到，錐柱之間的渦主要集中在波浪錐柱的中上部，由于波浪錐柱表面為波浪型，使得流過波浪錐柱表面的流體同時具有波浪型的特征，波浪錐柱后方渦出現(xiàn)按圓柱波長周期性分層現(xiàn)象。上游柱在后上方形成的大量渦團(tuán)，由于下游柱的存在，阻礙了其向正后方繼續(xù)發(fā)展，而是偏向一側(cè)，如圖7中View1，在下游柱迎風(fēng)面和側(cè)面均有大量渦聚集，這些渦部分形成于流體與下游柱表面的撞擊，部分來自上游柱向下發(fā)展的渦，因此在下游柱表面形成周期性的力，同時可以觀察到下游柱另一側(cè)面，如圖7中View2，也能觀察到周期性脫落的渦，同樣也能對下游柱表面產(chǎn)生周期性的力，兩種周期性的力相互耦合作用，極易產(chǎn)生大的脈動升力值。

圖7 α=10°時的Q準(zhǔn)則瞬時渦量圖(Q=1×106)Fig.7 Instantaneous vorticity of Q criterion at α=10°(Q=1×106)

4.4.2Z方向渦量圖

上游柱的尾流對下游柱作用非常復(fù)雜，為了更好觀察兩種間距比下的流動模式，選用圖8(a)中各個截面觀察流動情況。

由于波浪錐柱表面形狀的影響，流動結(jié)構(gòu)存在明顯分層，上游柱尾流對下游柱的作用展現(xiàn)為完全撞擊狀態(tài)，側(cè)面撞擊狀態(tài)和尾流干擾狀態(tài)。如圖8(b)所示(紅色線框表示完全撞擊狀態(tài)，橙色線框表示側(cè)面撞擊狀態(tài)，綠色線框表示尾流干擾狀態(tài))。完全撞擊狀態(tài)下，上游柱上側(cè)剪切層將下游柱迎風(fēng)面完全包裹，使下游柱表面阻力系數(shù)(前端壓力降低)和脈動升力系數(shù)(側(cè)面壓力降低)均有影響。而側(cè)面撞擊狀態(tài)下，上游柱尾流作用下游柱時，下游柱的側(cè)面受到上游柱的影響增強(qiáng)(側(cè)面壓力降低大)，而正面受到的影響則變?nèi)?正面壓力降低小)。側(cè)面壓力更大的增幅導(dǎo)致下游柱的脈動升力系數(shù)更大，如圖8(b)所示，間距比L/Dm=5相對于L/Dm=4擁有更多的側(cè)面撞擊狀態(tài)，這也表明該交錯角下的L/Dm=5的下游柱有更大的脈動升力系數(shù)。同時，間距比L/Dm=4相對于間距比L/Dm=5擁有更多的尾流包裹狀態(tài)，因此下游柱的時均阻力系數(shù)更低。

(a)截面示意圖

5 結(jié) 論

本文采用大渦模擬數(shù)值研究了亞臨界雷諾數(shù)(Re=3 900)下對間距比L/Dm=4和5，交錯角α=0°，5°，7.5°，10°，12.5°和15°錯列雙波浪錐柱繞流流動機(jī)理，探究了波浪錐柱平均阻力系數(shù)和脈動升力系數(shù)的變化規(guī)律，并對圓柱表面的時均壓力系數(shù)和流場結(jié)構(gòu)進(jìn)行了分析。主要結(jié)論如下：

1)由于上游柱的影響，下游柱時均阻力系數(shù)存在較大降幅，但隨著交錯角的增加，下游柱時均阻力系數(shù)也逐漸增大，趨近于上游柱。下游柱脈動升力系數(shù)較單圓柱均有大幅提升，且在交錯角α=10°時，兩種間距比下均能取得最大值，分別提升約20.1倍和21.4倍。

2)上游柱尾跡對下游柱的迎風(fēng)面和背風(fēng)面均有影響，當(dāng)交錯角較小時，正壓區(qū)和負(fù)壓區(qū)主要集中在柱體底部和受上游柱影響小的一側(cè)；當(dāng)交錯角逐漸增大時，正壓區(qū)和負(fù)壓區(qū)的范圍逐漸向圓柱上端擴(kuò)散，強(qiáng)度也有所增加。

3)波浪錐柱之間在Y截面存在兩個回流區(qū)，回流區(qū)中心向兩邊發(fā)展的渦與側(cè)面的渦相互作用形成“拱形”渦，存在回流區(qū)的波浪錐柱區(qū)域有較大的壓力系數(shù)變化。

4)在交錯角α=10°時，上游柱的下自由端上洗作用較上自由端下洗作用強(qiáng)，使得下游柱迎風(fēng)面上端有大量渦聚集；間距比L/Dm=5相對于L/Dm=4的下游柱有更多的側(cè)面撞擊狀態(tài)，脈動升力系數(shù)更大。本文的研究結(jié)果將對風(fēng)力俘能結(jié)構(gòu)布局提供有益理論支持。