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

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

(武漢理工大學 機電工程學院，武漢 430000)

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

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

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

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

1 控制方程

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

對于不可壓縮粘性流體，其Navier-Stokes運動方程

(1)

(2)

亞格子應力τij計算采用Smagorinsky-Lilly的亞格子模型：

(3)

(4)

(5)

(6)

(7)

式中：Δ表示過濾的網格尺度，Cs為Smagorinsky常數，本文取0.2。

2 模型驗證

計算域大小為30D×10D×10D，D為圓柱直徑，D=0.01 m，圓柱總長為H0，長徑比H0/D=7，笛卡爾坐標系的原點位于圓柱中心，圓柱中心距入口為10D，距離出口為20D，x軸正方向為順流向，z軸正方向沿圓柱軸心向上，如圖1所示。

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

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

本文計算壓力速度耦合采用PISO算法，離散化采取二階迎風格式[4]。表1為計算模型網格獨立性及時間步長無關性驗證結果，本文計算結果與實驗以及數值仿真結果基本一致，可以認為case2網格密度、無量綱時間步長Δt*=0.001以及本次仿真所用的算法正確，可以運用于錯列雙波浪錐柱仿真實驗。

表1 仿真結果與已有文獻對比Tab.1 Simulation results compared with previous literature

3 計算模型

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

(8)

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

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

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

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

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

4 結果分析

4.1 升阻力系數

升力系數Cd和阻力系數Cl的計算公式如下:

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

(9)

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

(10)

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

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

(a)時均阻力系數

4.2 波浪錐柱表面壓力系數

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

時均壓力系數計算公式為

(11)

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

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

4.3 流動結構

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

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

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

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

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

4.4 渦量圖

4.4.1 三維渦量圖

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

式中：Ω表示渦量，S表示應變率張量，‖ ‖為取范數運算符。

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

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

4.4.2Z方向渦量圖

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

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

(a)截面示意圖

5 結 論

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

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

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

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

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