基于格子Boltzmann方法的船舶風力助推轉子繞流特性

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

(江蘇科技大學 能源與動力學院，江蘇 鎮江 212000)

0 引 言

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

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

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

1 數值方法

1.1 多松弛模型

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

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

速度矩m與分布函數f之間的轉換如下：

對應的平衡矩為：

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

松弛矩陣為：

根據粒子碰撞過程的守恒定律有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ν由流體粘度系數所決定。

1.2 復雜曲線處理

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

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

其中：

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

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

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

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

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

1.3 數值驗證

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

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

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

2 計算區域

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

圖2 計算域Fig.2 Computational domain

設置Re= 100，轉速比|α| = 0.5，1.0，1.5 和 2.0，令上圓柱順時針轉動，下圓柱逆時針轉動。選取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 結果及分析

3.1 靜止并列雙圓柱

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

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

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

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

3.2 旋轉并列雙圓柱

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

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

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

圖4 間距比L/D = 1.7時不同轉速比下升力系數、阻力系數的變化曲線與對應的渦量云圖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時不同轉速比下升力系數、阻力系數的變化曲線與對應的渦量云圖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 結 語

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

1）對于單圓柱，基于自編程序對靜止和旋轉狀態分別進行了模擬，將結果與其他文獻中的數據進行對比，驗證程序的準確性。

圖6 不同間距比的時均升力系數和時均阻力系數隨轉速比的變化曲線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）并列雙圓柱的對轉可以明顯減輕間隙流對尾流的影響，同時也可以有效地抑制圓柱后渦的生成和脫落，并且當|α|達到臨界轉速比后，流場開始變得穩定，渦街隨之消失，給出了不同間距比對應的臨界轉速比|αc|。

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