基于格子Boltzmann方法的船舶風(fēng)力助推轉(zhuǎn)子繞流特性

穆 鑫，王 蛟，蘇石川，耿珊珊，Verendiaikin German，李 毅

(江蘇科技大學(xué) 能源與動力學(xué)院，江蘇 鎮(zhèn)江 212000)

0 引 言

風(fēng)帆曾是船舶最主要的推進方式，隨著全球能源短缺、環(huán)境污染等問題的日益凸顯，船舶航運利用風(fēng)能進行輔助推進再次獲得廣泛關(guān)注。船用風(fēng)力助推轉(zhuǎn)子，利用馬格努斯效應(yīng)，根據(jù)風(fēng)速和風(fēng)向的變化調(diào)整轉(zhuǎn)子的速度和旋轉(zhuǎn)方式，形成沿行駛方向的推力以達到降低推力成本目的[1]。其作為國際海事組織指定的節(jié)能技術(shù)，具有節(jié)能效果好，適用船型廣的特點，成為相關(guān)行業(yè)的研究熱點。國外已經(jīng)有數(shù)艘實船應(yīng)用的案例，國內(nèi)尚處于理論研究階段[2]。

船用風(fēng)力助推轉(zhuǎn)子通常是多個轉(zhuǎn)子組合使用，對船用風(fēng)力助推轉(zhuǎn)子的研究簡化為流體力學(xué)經(jīng)典問題—— 圓柱繞流。目前，對單個或多個圓柱繞流的流動特性已進行大量研究。雙圓柱作為多圓柱的基礎(chǔ)形式，其研究結(jié)果對多圓柱繞流有重要的理論意義。近年來，何穎等[3]對在亞臨界區(qū)均勻來流作用下的旋轉(zhuǎn)單圓柱繞流進行了數(shù)值模擬，結(jié)果證明圓柱旋轉(zhuǎn)可以大幅提高其升阻比。周凱等[4]對二維靜止串聯(lián)雙圓柱進行了數(shù)值模擬，重點研究圓柱間距對其升阻力及尾流特征的影響。巴悅等[5]基于Fluent對串聯(lián)旋轉(zhuǎn)雙圓柱的氣動性能進行了模擬，給出了較優(yōu)的排列旋轉(zhuǎn)方式。聶德明等[6]對靜止并列雙圓柱繞流進行了數(shù)值模擬，給出了圓柱間距對尾流及其升力、阻力的影響。Yoon等[7]對不同間距下的旋轉(zhuǎn)并列雙圓柱流動特性進行了模擬，探究了圓柱間距以及轉(zhuǎn)速比對繞流的影響，并且給出了不同間距下的臨界轉(zhuǎn)速。

格子Boltzmann方法（LBM）是20世紀80年代中期出現(xiàn)的一種新流體計算方法，具備流體相互作用描述簡單、易于編程、天然并行性等優(yōu)點。由于采用均勻規(guī)則的網(wǎng)格，LBM多被用于方形鈍體繞流的研究，對圓柱等復(fù)雜曲面的研究依然欠缺。近年來，隨著LBM理論的發(fā)展，能夠處理復(fù)雜曲面使得LBM得以快速被推廣。采用格子Boltzmann方法對并列旋轉(zhuǎn)雙圓柱進行模擬，對平直邊界采用經(jīng)典的反彈格式，而對于曲線邊界采用的是具有2階精度的反彈插值格式。研究雙圓柱在不同間距比、不同轉(zhuǎn)速比、不同轉(zhuǎn)向下流體的相互干涉效應(yīng)，以及這種干涉對雙圓柱升力、阻力系數(shù)的影響，為風(fēng)力轉(zhuǎn)子助推的控制與設(shè)計提供相關(guān)參考。

1 數(shù)值方法

1.1 多松弛模型

為提高LBM的模擬穩(wěn)定性，采用多松弛格子Boltzmann模型，速度離散模型采用文獻[8]提出的D2Q9模型。分布函數(shù)f通過線性變換到矩空間m，從而建立速度空間與矩空間的聯(lián)系。每一個速度矩代表一個實際的物理量，不同的物理量采用不用的松弛參數(shù)。MRT-LBM的演化方程為：

式中：f為分布函數(shù)；e為離散速度；M為變換矩陣；具體形式參見文獻[9]。S為松弛參數(shù)對角矩陣；下標i，j表示離散方向（下同）；m為分布函數(shù)對應(yīng)的速度矩；meq為平衡矩。

速度矩m與分布函數(shù)f之間的轉(zhuǎn)換如下：

對應(yīng)的平衡矩為：

式中：ρ為密度；ux和uy分別為速度u在x和y坐標方向上的分量。

松弛矩陣為：

根據(jù)粒子碰撞過程的守恒定律有sρ（ρ-ρ） = 0，sj（jx-jx） = 0 和sj（jy-jy） = 0，因此sρ和sj可以選取任意值，選取sρ=sj= 0，se= 1.2，sε= 1.1，sq= 1.1，sν由流體粘度系數(shù)所決定。

1.2 復(fù)雜曲線處理

如圖1所示的任意曲線邊界，實心圓點xw代表的是固體邊界與格子連線的交點，空心圓點xf代表流體節(jié)點，空心方形xb代表固體內(nèi)的邊界節(jié)點（固體內(nèi)臨近邊界的節(jié)點），空心菱形xs代表固體節(jié)點。定義q為：

可知 0 ≤q≤1，其中xf，xb，xw分別表示網(wǎng)格節(jié)點f，w，b的坐標。定義由流體節(jié)點流向邊界節(jié)點的粒子速度方向為ei，則反方向為eˉi。為了確定xf點處的分布函數(shù)，必出先對xb點進行插值處理。采用在碰撞之后，遷移之前對xb點進行線性插值[10]：

其中：

式中：wi為權(quán)重因子；f‘為碰撞后的分布函數(shù)；ubf為虛擬速度；uf為粒子位于xf點處的速度；uw為運動固體邊界xw處速度；χ為由q確定的加權(quán)。定義轉(zhuǎn)速比|α| =ΩD/（2U），其中Ω表示圓柱的角速度，U為來流速度。

采用動量轉(zhuǎn)換法計算圓柱受到的作用力。其具體形式為：

圖1 空間格子和曲線壁面的分布Fig.1 Distribution of spatial lattices and curved wall

圓柱的阻力系數(shù)CD以及升力系數(shù)CL表達式如下:

式中：FD和FL分別為圓柱所受流體作用力F在順流方向和法線方向上的分量。

1.3 數(shù)值驗證

為了驗證程序的正確性以及計算結(jié)果的精度，首先模擬了靜止及旋轉(zhuǎn)狀態(tài)下的單圓柱繞流。定義Strouhal數(shù)St=Dfv/U，其中fv為渦脫落頻率。表1給出了在Re=100的情況下計算結(jié)果，并與參考文獻中的數(shù)據(jù)進行對比。

表1 模擬結(jié)果與參考文獻結(jié)果對比Tab.1 Comparison of simulation results with others

由表1可知，靜止狀態(tài)時本文模擬所得平均阻力系數(shù)和St數(shù)與實驗結(jié)果擬合較好，誤差控制在2%以內(nèi)。St數(shù)與同樣采用模擬手段文獻[13]的結(jié)果基本相同，平均阻力系數(shù)誤差也小于3%。而當圓柱旋轉(zhuǎn)狀態(tài)時，本文模擬得到的平均阻力系數(shù)與文獻數(shù)據(jù)依然吻合較好，與Stojkovi?等模擬結(jié)果最大偏差為6.3%。同時，St數(shù)與參考數(shù)據(jù)基本相同，最大誤差僅為3.5%。從而可知該自編程序?qū)τ谀M圓柱繞流具有良好的可靠性。

2 計算區(qū)域

圓柱直徑D為流動特征長度，令兩圓柱圓心距為L，如圖2所示，計算區(qū)域為20D×50D，兩圓柱距入口長度為10D，距出口長度為40D。

圖2 計算域Fig.2 Computational domain

設(shè)置Re= 100，轉(zhuǎn)速比|α| = 0.5，1.0，1.5 和 2.0，令上圓柱順時針轉(zhuǎn)動，下圓柱逆時針轉(zhuǎn)動。選取L/D=1.2， 1.7， 2.5和4.0四個典型間距比。圖中U為來流速度，流向速度和法向速度分別定義為u和v。上下壁面采用無滑移邊界條件，進出口邊界條件定義如下：

進口u=U=0.1，v=0;

出口?u/?x=0，?v/?x=0。

3 結(jié)果及分析

3.1 靜止并列雙圓柱

圖3給出了不同間距比并列雙圓柱在靜止狀態(tài)（|α|=0）下的升阻力系數(shù)變化曲線及對應(yīng)的渦量云圖?？梢钥吹介g距比對并列雙圓柱繞流的尾流有明顯的影響。從圖3（a）中渦量云圖可以看到，當間距比很小時，間隙流對圓柱后尾流幾乎沒有影響。在雙圓柱的上下兩側(cè)交替生成漩渦脫落，并在其尾流只形成一個渦街。相應(yīng)的流動穩(wěn)定后圓柱升阻力系數(shù)隨無量綱時間（t=tU/D，下同）呈現(xiàn)周期性變化。

隨著間距比的增大（1.2＜L/D＜2.2），從圖3（b）中可知間隙流對圓柱后尾流的影響增強，導(dǎo)致單渦街結(jié)構(gòu)被破壞，尾流不再具有周期性，間隙流交替地偏向一側(cè)，呈現(xiàn)出偏流模式，同時可見在兩圓柱后形成寬的和窄的尾流結(jié)構(gòu)。此模式下兩圓柱的升阻力系數(shù)的變化不再具有周期性，可見兩圓柱的阻力系數(shù)交替變化，間隙流偏向的圓柱其阻力系數(shù)較大。

當間距比持續(xù)增大（L/D＞2.2），由圖3（c）和圖3（d）可知尾流再次呈現(xiàn)出周期性，并且2列旋渦同步脫落，稱為對稱流模式。對稱流模式又分為同步同相和同步反相。如圖3（c）所示，對于同步同相模式，兩圓柱表面形成相位相同的旋渦，在距離圓柱一段距離后，符號相同的渦會融合成尺度更大的渦，說明此時尾流中的渦結(jié)構(gòu)是不穩(wěn)定的。相應(yīng)圓柱的升力系數(shù)隨時間的變化曲線相位是相同的。而從圖3（d）可見圓柱的升力系數(shù)變化曲線相位差180°，表明此時旋渦同步反相脫落。從渦量云圖中可以看到兩圓柱后緣生成同步反相的旋渦，在圓柱后同時存在2列關(guān)于流場中線對稱的平行反相渦街。不同于同步同相模式的是此時2列渦街在很長一段距離內(nèi)一直保持穩(wěn)定，并且與單圓柱繞流的尾流結(jié)構(gòu)相似。上述模擬結(jié)果與Sumner D.等[18]的結(jié)論相符。

圖3 靜止狀態(tài)下不同模式的升力系數(shù)和阻力系數(shù)變化曲線及對應(yīng)的渦量云圖Fig.3 The curves of lift and drag coefficient of different modes under stationary state and the corresponding vorticity cloud diagram

3.2 旋轉(zhuǎn)并列雙圓柱

圓柱旋轉(zhuǎn)將帶動其周圍流體一起運動，對于均勻來流的流場，將在圓柱兩側(cè)對流體分別產(chǎn)生正的和負的作用。根據(jù)馬格努斯效應(yīng)，由流速差引起的壓力差將會在圓柱上產(chǎn)生與來流方向垂直的橫向力。本文的并列雙圓柱采用上順下逆的對轉(zhuǎn)方式，可知上圓柱的升力方向為垂直向上，下圓柱的升力方向則為垂直向下。從圖3可見，對于不同轉(zhuǎn)速比下上圓柱的升力系數(shù)大于零，而下圓柱的升力系數(shù)小于零，即上下圓柱所受升力方向相反，與上述分析結(jié)論相符，同時也表明此種對轉(zhuǎn)方式兩圓柱之間存在排斥力。

圖4為間距比L/D= 1.7時，升力系數(shù)與阻力系數(shù)的變化曲線及對應(yīng)的渦量云圖。對比圖3（b）和圖4（a）可知，對于偏流模式，轉(zhuǎn)速比增加到1.0時，圓柱后的尾流被拉長，并且在圓柱后出現(xiàn)單渦街流型，尾流由偏流模式轉(zhuǎn)變?yōu)閱螠u街模式。同時從升阻力系數(shù)曲線圖可以看到，在流動穩(wěn)定后，圓柱的升阻力系數(shù)曲線呈現(xiàn)出周期性。當轉(zhuǎn)速比增大到|α| = 1.2時，從圖4（b）可以看到2個圓柱后尾流渦變得狹長并且不再有漩渦脫落，同時可以看到兩圓柱的阻力系數(shù)幾乎重疊以及升阻力系數(shù)曲線的脈動值近乎為零，流動穩(wěn)定后的升阻力系數(shù)曲線接近一條直線，可知間距比為1.7時的臨界轉(zhuǎn)速比|αc| = 1.2。由圖4（c）的渦量云圖可以看得，當轉(zhuǎn)速比進一步增大（超過臨界轉(zhuǎn)速比）時，兩圓柱之間的間隙流被完全削弱，對尾流結(jié)構(gòu)不再有影響作用，流場變得穩(wěn)定不再有漩渦產(chǎn)生和脫落。升、阻力系數(shù)變化曲線的趨勢與圖4（b）一致。

圖5給出了間距比L/D= 2.5時，升力系數(shù)與阻力系數(shù)的變化曲線以及對應(yīng)的渦量云圖。由前文可知圓柱靜止時的流場為對稱流中的同步同相模式，從圖5（a）渦量云圖可以明顯看到當轉(zhuǎn)速比為1.0時圓柱后的尾流呈現(xiàn)同步反相模式，從對應(yīng)的升阻力系數(shù)曲線也可以得到證實。之后升阻力系數(shù)曲線以及渦量云圖隨轉(zhuǎn)速比增大的變化趨勢與間距比為1.7時相似。結(jié)合圖4和圖5可以知道，并列圓柱的旋轉(zhuǎn)可以有效地抑制漩渦的產(chǎn)生和脫落，甚至徹底消除渦街的生成。間距比L/D為1.2和4.0隨轉(zhuǎn)速比變化而變化的趨勢與上述兩者類似。根據(jù)本次的模擬結(jié)果證實了不同間距下臨界轉(zhuǎn)速比的存在，對于選取的間距比L/D= 1.2，1.7，2.5和4.0，對應(yīng)的臨界轉(zhuǎn)速比分別為|αc| = 1.2，1.2，1.3和 1.7。

圖4 間距比L/D = 1.7時不同轉(zhuǎn)速比下升力系數(shù)、阻力系數(shù)的變化曲線與對應(yīng)的渦量云圖Fig.4 The curves of lift and drag coefficient and the corresponding vorticity cloud diagram at different speed ratios when L/D=1.7

圖5 間距比L/D = 2.5時不同轉(zhuǎn)速比下升力系數(shù)、阻力系數(shù)的變化曲線與對應(yīng)的渦量云圖Fig.5 The curves of lift and drag coefficient and the corresponding vorticity cloud diagram at different speed ratios when L/D=2.5

4 結(jié) 語

風(fēng)力助推轉(zhuǎn)子能夠合理有效地利用風(fēng)能，進一步完善船舶推進技術(shù)。本文基于格子Boltzmann方法，對不同間距比和轉(zhuǎn)速比下的并列轉(zhuǎn)子繞流進行了研究。文中主要對升力系數(shù)、阻力系數(shù)以及尾流流型等數(shù)據(jù)進行處理和分析，得出如下結(jié)論：

1）對于單圓柱，基于自編程序?qū)o止和旋轉(zhuǎn)狀態(tài)分別進行了模擬，將結(jié)果與其他文獻中的數(shù)據(jù)進行對比，驗證程序的準確性。

圖6 不同間距比的時均升力系數(shù)和時均阻力系數(shù)隨轉(zhuǎn)速比的變化曲線Fig.6 The curves of the mean lift and mean drag coefficient with different rotational speeds at different spacing ratios

2）對于靜止并列雙圓柱，隨著間距比的不同存在3種典型的尾流流型： 1.0＜L/D≤1.2時，單渦街流型；1.2＜L/D≤2.2時，偏流流型；L/D＞2.2時，同步對稱流型，其包括同步同相和同步反相2種流型。

3）并列雙圓柱的對轉(zhuǎn)可以明顯減輕間隙流對尾流的影響，同時也可以有效地抑制圓柱后渦的生成和脫落，并且當|α|達到臨界轉(zhuǎn)速比后，流場開始變得穩(wěn)定，渦街隨之消失，給出了不同間距比對應(yīng)的臨界轉(zhuǎn)速比|αc|。

4）兩圓柱時均升力系數(shù)的絕對值隨轉(zhuǎn)速比的增大而增大，而時均阻力系數(shù)則隨著轉(zhuǎn)速比的增大而減小。時均阻力系數(shù)的脈動值隨著轉(zhuǎn)速比達到|αc|后出現(xiàn)變小的趨勢。