王露，李天勻，朱翔，郭文杰

1華中科技大學船舶與海洋工程學院，湖北武漢430074

2華中科技大學船舶和海洋水動力湖北省重點實驗室，湖北武漢430074

浸沒邊界格子Boltzmann方法的改進及轉(zhuǎn)動圓柱繞流模擬

王露1，2，李天勻1，2，朱翔1，2，郭文杰1，2

1華中科技大學船舶與海洋工程學院，湖北武漢430074

2華中科技大學船舶和海洋水動力湖北省重點實驗室，湖北武漢430074

針對浸沒邊界格子Boltzmann方法計算效率不足的問題，提出一種改進浸沒邊界格子Boltzmann方法。通過使用作用力格子Boltzmann模型，并簡化流固耦合作用力的計算方法，使計算流程得到簡化，同時，對比原方法轉(zhuǎn)動圓柱繞流計算結(jié)果，計算時間縮短一半以上，提高了數(shù)值計算效率。結(jié)合改進算法研究轉(zhuǎn)動圓柱繞流的流場升阻力系數(shù)、壓力系數(shù)、流場速度等流體參數(shù)場隨轉(zhuǎn)速比的關(guān)系，揭示了一定的轉(zhuǎn)動圓柱繞流動力學規(guī)律。研究結(jié)果表明，改進方法準確可靠。

流固耦合；格子Boltzmann方法；浸沒邊界法；轉(zhuǎn)動圓柱繞流

0 引 言

流體中的圓柱轉(zhuǎn)動是學術(shù)研究和實際工程中十分重要的研究對象。圓柱在流體中的運動過程包含在船舶艉軸承推進和軸系冷卻等實際工程中。在圓柱繞流的問題中，圓柱繞流一般會產(chǎn)生漩渦脫落并誘發(fā)流體的脈動載荷與結(jié)構(gòu)振動。渦激振動不僅對物體造成疲勞，還會對結(jié)構(gòu)造成破壞。針對這些情況，許多學者就如何抑制流體中的渦脫落提出了控制方法，其中圓柱轉(zhuǎn)動就是其中之一。本文將利用轉(zhuǎn)動圓柱繞流的渦脫落現(xiàn)象及升阻力現(xiàn)象的研究成果，開展水潤滑軸承的水膜動力學特性分析研究，為軸系振動控制提供基礎(chǔ)。

隨著數(shù)值仿真方法的興起與發(fā)展，對于轉(zhuǎn)動圓柱的流固耦合現(xiàn)象的模擬與研究越發(fā)深入。2004年，F(xiàn)eng等［8］首先提出了將格子Boltzmann方法（LBM）［9］與浸沒邊界法（IBM）［10］相結(jié)合的浸沒邊界格子Boltzmann方法（IB-LBM），并成功地進行了流固耦合仿真。接著Niu等［11］提出了基于動量交換的IB-LBM，Wu等［12］提出了一種隱式速度修正的IB-LBM，然而在一些傳統(tǒng)的IB-LBM中，有時在處理流體與固體間的作用力時對網(wǎng)格的尺寸限制較大，有時在處理流體與固體間節(jié)點相互作用力時十分復雜，這些局限使得IB-LBM的計算量會增大。Suzuki等［13］提出的相互迭代的IB-LBM方法在處理流固耦合作用力時具有簡單的特性，同時，在運動固體的流固耦合模擬中顯示出了較高的精度。但此方法在處理流固耦合作用力時需要迭代逼近來保證收斂，會增加計算成本。

針對轉(zhuǎn)動圓柱繞流問題，本文將通過改進流固耦合作用力的計算方式，提出一種簡單高效的浸沒邊界格子Boltzmann方法，并通過轉(zhuǎn)動圓柱繞流模型驗證改進方法的有效性及高效性。同時，為了研究轉(zhuǎn)動圓柱繞流的水動力學參數(shù)變化特性，本文分別研究Re為40和200時，轉(zhuǎn)動圓柱繞流隨轉(zhuǎn)速比變化的流場，及相關(guān)流體動力學參數(shù)變化的特性，可為工程實際應用提供參考。

1 浸沒邊界格子Boltzmann法

1.1 格子Boltzmann方法

格子Boltzmann方程描述的是具有離散速度的流體粒子分布函數(shù)在固定格點上的運動過程。對于經(jīng)典格子Boltzmann模型——LBGK模型［14］，在流體碰撞過程中，對應的格子Boltzmann方程表達式為

式中：ωi為離散速度的權(quán)系數(shù)；u為格點遷移速度；對應的宏觀量密度和動量為：

式中，n為離散速度方向的總個數(shù)。

式（2）～式（4）構(gòu)成了經(jīng)典的LBGK模型。

對于二維D2Q9模型，n為9，其模型離散的速度分布如圖1所示，對應的離散速度表達式為：

圖1 D2Q9離散速度分布圖Fig.1 Discrete speed distribution of D2Q9 model

1.2 浸沒邊界法

Peskin［10］提出的浸沒邊界法簡化了流固耦合作用力的處理（圖2）。在浸沒邊界法中，對于固體和流體分別采用拉格朗日坐標系和歐拉坐標系進行描述。首先建立描述固體的拉格朗日坐標系，對應的固體節(jié)點集為Gh，用X(q ，r，s，t)表示在時間t下的拉格朗日節(jié)點，F(xiàn)(q ，r，s，t)為相應時刻固體邊界受到的流體作用力，U(X ，t)為固體邊界的速度。相應地采用歐拉坐標系描述流場，對應的流體節(jié)點集為gh。

圖2 流體及固體邊界節(jié)點示意圖Fig.2 Schematic of fluid and solid boundary

1.3 改進方法

Suzuki等［13］的IB-LBM方法在處理流固耦合作用力時，需采用迭代過程，使得計算速度下降，增加計算時間。其主要原因為在流體演化過程中采用的是標準LBGK模型，并沒有考慮流體受到固體作用力的影響，從而使得演化后的流體速度分布函數(shù)與實際速度分布函數(shù)存在誤差。

為了提高計算時流場速度分布函數(shù)的計算精度，Guo等［15］在處理Boltzmann方程時，考慮到外力項對流體演化過程的影響，推導出具有二階精度的作用力格子Boltzmann模型，本文采用該作用力模型處理流場的演化過程，從而提高流場演化過程中的分布函數(shù)精度，其對應的方程為

對應的動量表達式為

對應的改進方法所采用的拉格朗日節(jié)點作用的計算方法如下：

式中：Ub為拉格朗日節(jié)點實際遷移速度；Sh為尺度因子。采用了作用力格子Boltzmann模型后，計算步驟如下：

1）給定流場的初始值 ρ，u，利用式（2）計算出對應的初始平衡分布函數(shù)計算拉格朗日節(jié)點運動變量

2）通過式（9），計算由t到t+δt時刻，流場發(fā)生碰撞和遷移過程后的分布函數(shù)運用式（3），（4）和（11），計算t+δt時刻對應的ρ，u*和u。

3）通過浸沒邊界法提供的式（12），計算插值出的拉格朗日邊界節(jié)點速度U*(X，t+δt)。

5）重復步驟2）～4），直至收斂。

式中：P(s,t)為瞬變電磁場分量；U(s,q)為轉(zhuǎn)換波場分量；s為波場分量對應的時間的場值；t為瞬變電磁場時間；q為波場類時間。

2 轉(zhuǎn)動圓柱繞流模擬

2.1 方法驗證

針對改進方法研究轉(zhuǎn)動圓柱繞流模擬，本文仿真過程采用的計算物理模型如圖3所示，其中，圓柱轉(zhuǎn)動方向為順時針，流體沿x正方向流入。流域尺寸為 30D×15D，圓柱圓心坐標為(6.25D，7.5D )，其中D為圓柱直徑?？紤]到計算精度，圓柱直徑D=24δx，仿真過程中選取的格子尺寸δx=δt=1。流域的邊界條件為：左邊入口u=U，v=0；右邊出口即沿x方向速度不變；上下邊界即沿 y方向速度不變，其中U=0.1，為流域入口速度。為了處理格子Boltzmann碰撞遷移過程的邊界條件，流場邊界采用的是非平衡態(tài)外推邊界條件［16］。對于圓柱邊界將其等分為120段。

圖3 仿真計算模型Fig.3 Simulation model

計算圓柱在流場中所受到的升阻力時，直接取為拉格朗日節(jié)點所受的流體作用力。相對應的阻力系數(shù)和升力系數(shù)可表示為：

式中：Cd，Cl分別為阻力系數(shù)和升力系數(shù)；Fx，F(xiàn)y分別為拉格朗日節(jié)點作用力F在 x方向和 y方向上的分量。

圖4為雷諾數(shù)Re=200，轉(zhuǎn)速比α分別為0.5，1，1.5時，圓柱升阻力系數(shù)關(guān)系曲線，對比左側(cè)的渦量云圖，可以發(fā)現(xiàn)圓柱尾部渦量隨著轉(zhuǎn)速比的增加圓柱尾渦逐漸下移。隨著轉(zhuǎn)速比的增加，圓柱的升阻力系數(shù)曲線逐漸變化為橢圓形，這一曲線變化現(xiàn)象與Mittal等［17］的有限元仿真結(jié)果一致。

圖4 Re=200，α分別為0.5，1，1.5時渦量云圖及圓柱升阻力曲線Fig.4 Vorticity contours and the drag coefficient-lift coefficient line of flow over rotating cylinder when Re is 200，αis 0.5,1.0,1.5

圖4右側(cè)的升阻力系數(shù)關(guān)系曲線，是由計算穩(wěn)定后，10 000時間步下的數(shù)值點構(gòu)成，由多周期曲線構(gòu)成。從圖中可以看出，不同周期下的節(jié)點在相同位置重合，表明改進算法計算結(jié)果穩(wěn)定。

表1 圓柱阻力系數(shù)和升力系數(shù)比較Tab.1 Comparison between drag coefficient and lift coefficient

表1為不同轉(zhuǎn)速比下，本文的模擬結(jié)果與鄔小軍等［18-19］采用格子Boltzmann方法仿真的升力系數(shù)和阻力系數(shù)對比。比較表1中給出的4組數(shù)據(jù)，改進算法采用的是浸沒邊界法計算固體邊界受力，而另外3種算法則采用的是3種不同的LB曲面邊界處理格式，分別為統(tǒng)一邊界法（YMS）、非平衡態(tài)外推格式法和反彈邊界條件法，由于這些邊界條件的處理方式均采用數(shù)值插值格式方式計算固體邊界流場的速度分布，在計算過程中會產(chǎn)生誤差，而固體表面受力計算過程是通過固體表面流場速度分布變化（浸沒邊界法）或固體表面流場動量分布變化（格子Boltzmann方法）計算，從而會進一步產(chǎn)生數(shù)值誤差，導致各文獻間的數(shù)據(jù)存在一定誤差。但從表中數(shù)據(jù)可發(fā)現(xiàn)，升力系數(shù)與阻力系數(shù)變化趨勢基本吻合，對應的幅值變化也基本一致，表明通過與退化為格子Boltzmann方法計算的結(jié)果對比，本文的改進方法在轉(zhuǎn)動圓柱繞流仿真的可行性得到驗證。

為進一步驗證本文的改進方法在轉(zhuǎn)動圓柱繞流中的有效性，將本方法計算流場結(jié)果與Yu等［20］基于格子Boltzmann方法提出的YMS的計算結(jié)果進行對比。圖5給出了轉(zhuǎn)速比α=0.5時，當流場穩(wěn)定后，半個周期內(nèi)不同時刻，在圓心高度上，從圓柱邊界到圓心2個直徑范圍內(nèi)的速度分布。從圖中可以看出，二者在x及y方向不同時間的速度分布結(jié)果十分接近，這也進一步驗證了本方法的在轉(zhuǎn)動圓柱繞流仿真中的有效性。

圖5 Re=200，α=0.5，圓柱后方速度分布曲線Fig.5 The velocity distribution of flow over rotating cylinder when Re is 200，α is 0.5

為比較改進計算方法與原方法的計算效率，本文模型仿真的電腦參數(shù)為CPU為E5-2650 V2 2.6 GHz，內(nèi)存為128 GB，在2012a版本的Matlab上仿真2種方法的轉(zhuǎn)動圓柱繞流。表2為Re為40和200，α=0.5這2種工況下，為保證計算收斂時仿真計算時間對比表。對比2組計算數(shù)據(jù)可以發(fā)現(xiàn)，相同循環(huán)次數(shù)下，改進計算方法所需計算時間縮短了一半以上，這充分體現(xiàn)了改進方法計算效率上的優(yōu)勢。

表2 Re為40和200計算模型對應的計算時間對比Tab.2 Comparison of the calculating time for flow over a cylinder when Re is 40 and 200

2.2 計算結(jié)果

圖6為不同轉(zhuǎn)速比下，Re=200時，圓柱阻力系數(shù)、升力系數(shù)隨時間變化的曲線。觀察圖6（a）和圖6（b）中曲線的變化過程，可以發(fā)現(xiàn)當轉(zhuǎn)速比α＜1.0時，阻力曲線后半周期中對應處的波谷會逐漸上移。當α＞1.0后圓柱阻力曲線周期會增大到約為原來的2倍，即原來的波谷出現(xiàn)消失。圓柱的阻力系數(shù)幅值明顯增大，平均阻力系數(shù)也明顯減小，隨著轉(zhuǎn)速比的增大圓柱受到的升力逐漸增大。隨著轉(zhuǎn)速比的增加，圓柱的升力系數(shù)曲線逐漸向上平移，幅值并沒有觀察到明顯變化。同時對比各條曲線的波谷波峰可以發(fā)現(xiàn)，各曲線的相位并不相同，這是由于圓柱轉(zhuǎn)動產(chǎn)生的相位滯后。

為研究圓柱表面正壓力在圓柱表面的分布與轉(zhuǎn)速比的關(guān)系，引入圓柱表面壓力系數(shù)Cp作為研究參數(shù)。其表達式為

式中：Pθ為圓柱表面θ處的壓力；P0為無窮遠處來流壓力。

圖7是Re為40和200時，本文和Suzuki［13］這2種計算方法計算的不同轉(zhuǎn)速比下的圓柱表面壓力系數(shù)分布。由圖中的峰值曲線可以看出，2種計算結(jié)果的壓力系數(shù)分布十分接近，即本計算方法有效。角度接近90°在圓柱正下方時圓柱表面壓力系數(shù)最小，且隨著轉(zhuǎn)速比的增大峰值出現(xiàn)明顯變化，這與實際的流場壓力分布變化一致。

圖7 Re為40和200的不同轉(zhuǎn)速比下，圓周表面壓力系數(shù)分布Fig.7 The pressure coefficient along the cylinder surface when Re is 40 and 200 with differentα

圖8為升阻力系數(shù)隨轉(zhuǎn)速比變化曲線。隨著α增大，當升力系數(shù)較小且α＜0.5時，阻力系數(shù)并沒有明顯改變，當1.0＜α＜2.4時，阻力系數(shù)呈近似線性減小，同時升力系數(shù)隨轉(zhuǎn)速比增加呈線性變化，但Re為200與40的計算結(jié)果相比較，阻力系數(shù)變化更劇烈，這一變化趨勢與Mittal等［17］的有限元仿真結(jié)果一致。

圖8 Re為40和200時圓柱升阻力系數(shù)隨轉(zhuǎn)速比變化Fig.8 The drag coefficient and lift coefficient of flow over a rotating cylinder when Re is 40 and 200 with differentα

3 結(jié) 語

本文通過結(jié)合作用力格子Boltzmann模型和Suzuki等［13］浸沒邊界法模型，提出了一種計算效率更高的浸沒邊界格子Boltzmann改進算法。通過轉(zhuǎn)動圓柱繞流的模型，對比不同文獻與本方法的仿真結(jié)果的阻力系數(shù)、升力系數(shù)及流場分布等流體動力學參數(shù)，證明了改進計算方法的有效

性。同時，對比不同雷諾數(shù)轉(zhuǎn)動圓柱繞流計算時間，發(fā)現(xiàn)改進計算方法節(jié)省一半以上的時間，充分說明改進算法高效。針對轉(zhuǎn)動圓柱繞流進行了仿真，并對隨轉(zhuǎn)速比變化的流體動力學參數(shù)變化進行了分析，其結(jié)果可為相關(guān)的工程應用提供參考。

［1］ PRANDTL L.The magnus effect and wind-powered ships［J］.Naturwissenschaften，1925（13）：93-108

［2］ TOKUMARU P T，DIMOTAKIS P E.The lift of a cylinder executing rotary motions in an uniform flow［J］. Journal of Fluid Mechanics，1993（255）：1-10.

［3］ CHEN Y M，OU Y R，PEARLSTEIN A J.Development of the wake behind a circular cylinder impulsively started into rotatory and rectilinear motion［J］.Journal of Fluid Mechanics，1993（253）：449-484.

［4］ CHEW Y T，CHENG M，LUO S C.A numerical study of flow past a rotating circular cylinder using a hybrid vortex scheme［J］.Journal of Fluid Mechanics，1995（299）：35-71.

［5］ KANG S，CHOI H，LEE S.Laminar flow past a rotating circular cylinder［J］.Physics of Fluids（1994 present），1999，11（11）：3312-3321.

［6］ 何穎，楊新民，陳志華，等.旋轉(zhuǎn)圓柱繞流的流場特性［J］.船舶力學，2015，19（5）：501-508. HE Ying，YANG Xinmin，CHEN Zhihua，et al.Flow field characteristics of flow past a rotating cylinder［J］. Journal of Ship Mechanics，2015，19（5）：501-508.

［7］ 赫鵬，李國棟，楊蘭，等.圓柱繞流流場結(jié)構(gòu)的大渦模擬研究［J］.應用力學學報，2012，29（4）：437-443. HAO Peng，LI Guodong，YANG Lan，et al.Large eddy simulation of the circular cylinder flow in different regimes［J］.Chinese Journal of Applied Mechanics，2012，29（4）：437-443.

［8］ FENG Z G，MICHAELIDES E E.The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems［J］.Journal of Computational Physics，2004，195（2）：602-628.

［9］ SUCCI S.The lattice Boltzmann equation：for fluid dynamics and beyond［M］.Oxford：Oxford University Press，2001.

［10］ PESKIN C S.The immersed boundary method［J］.Acta Numerica，2002，11：479-517.

［11］ NIU X D，SHU C，CHEW Y T，et al.A momentum exchange-basedimmersedboundary-latticeBoltzmann method for simulating incompressible viscous flows［J］.Physics Letters A，2006，354（3）：173-182.

［12］ WU J，SHU C.Implicit velocity correction-based immersed boundary-lattice Boltzmann method and its applications［J］.Journal of Computational Physics，2009，228（6）：1963-1979.

［13］ SUZUKI K，INAMURO T.A higher-order immersed boundary-lattice Boltzmann method using a smooth velocity field near boundaries［J］.Computers and Fluids，2013，76：105-115.

［14］ QIAN Y H，D'HUMIèRES D，LALLEMAND P.Lattice BGK models for Navier-Stokes equation［J］.Europhysics Letters，1992，17（6）：479-484.

［15］ GUO Z L，ZHENG C G，SHI B C.Discrete lattice effects on the forcing term in the lattice Boltzmann method［J］.Physical Review E，2002，65（4）：046308.

［16］ GUO Z L，ZHENG C G，SHI B C.Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method［J］. Chinese Physics，2002，11（4）：366-374.

［17］ MITTAL S，KUMAR B.Flow past a rotating cylinder［J］.Journal of Fluid Mechanics，2003，476：303-334.

［18］ 鄔小軍.靜止和轉(zhuǎn)動圓柱繞流的格子Boltzmann模擬［D］.武漢：華中科技大學，2009. WU Xiaojun.The lattice Boltzmann method simulation of the flow past static and rotating circular cylinder［D］.Wuhan：Huazhong University of Science and Technology，2009

［19］ YAN Y Y，ZU Y Q.Numerical simulation of heat transfer and fluid flow past a rotating isothermal cylinder-a LBM approach［J］.International Journal of Heat and Mass Transfer，2008，51（9/10）：2519-2536.

［20］ YU D Z，MEI R W，SHYY W.A unified boundary treatment in lattice Boltzmann method［C］//Proceedings of the 41st Aerospace Sciences Meeting and Exhibit.Reno，Nevada：AIAA，2003.

Improved immersed boundary lattice Boltzmann method and simulation of flow over rotating cylinder

WANG Lu1，2，LI Tianyun1，2，ZHU Xiang1，2，GUO Wenjie1，2

1 School of Naval Architecture and Ocean Engineering，Huazhong University of Science and Technology，Wuhan 430074，China

2 Hubei Key Laboratory of Naval Architecture and Ocean Engineering Hydrodynamics，Huazhong University of Technology and Science，Wuhan 430074，China

Based on Suzuki's immersed boundary lattice Boltzmann method,an improved approach is proposed.By adopting the force model of the lattice Boltzmann method and improving the technique of fluid structure interaction force calculation,the process of calculating the force of fluid structure interaction is simplified.Moreover，its calculation time is over 50%shorter than that of the original method,which means that efficiency is greatly improved.Based on the improved method,such features as lift coefficient,drag coefficient，pressure coefficient and flow field are discussed under different conditions of flow over a rotating cylinder.Compared with the related results,this new improved method is shown to be accurate and reliable.

fluid structure interaction；lattice Boltzmann method；immersed boundary method；flow around rotating cylinder

U661.1

A

10.3969/j.issn.1673-3185.2016.06.015

2015-12-17

時間：2016-11-18 15:19

國家自然科學基金資助項目（51379083，51479079，51579109）；高等學校博士學科點專項科研基金資助項目（20120142110051）

王露，男，1991年生，碩士生。研究方向：船舶與海洋結(jié)構(gòu)物設(shè)計制造。E-mail：1229398383@qq.com李天勻（通信作者），男，1969年生，教授，博士生導師。研究方向：結(jié)構(gòu)振動噪聲分析。E-mail：ltyz801@hust.edu.cn

http://www.cnki.net/kcms/detail/42.1755.tj.20161118.1519.030.html 期刊網(wǎng)址：www.ship-research.com

王露，李天勻，朱翔，等.浸沒邊界格子Boltzmann方法的改進及轉(zhuǎn)動圓柱繞流模擬［J］.中國艦船研究，2016，11（6）：97-103. WANG Lu，LI Tianyun，ZHU Xiang，et al.Improved immersed boundary lattice Boltzmann method and simulation of flow over rotating cylinde［rJ］.Chinese Journal of Ship Research，2016，11（6）：97-103.