格子Boltzmann方法在串列雙圓柱繞流數值模擬中的應用研究

周 凱，王 震，陳維山，龍曉軍

（1.山東農業大學 機械與電子工程學院，山東 泰安 271018；2.哈爾濱工業大學 機器人技術與系統國家重點實驗室，哈爾濱 150001）

0 引 言

格子Boltzmann方法是一種不同于傳統數值方法的流體建模和計算方法[1-3]。傳統的數值方法的基本思路是將宏觀控制方程進行離散，然后利用數值求解的方法去求解該方程，例如有限差分、有限體積、有限元法等。格子Boltzmann方法的是基于分子運動學和統計力學，在微觀粒子尺度上建立離散速度模型。在滿足質量、動量和能量守恒的基本條件下，通過粒子遷移和碰撞兩個微觀行為來反應宏觀特性，這使得算法簡化為對粒子遷移和碰撞的模擬。格子Boltzmann方法采用均勻的正交網格，使得網格算法的工作量大大簡化，特別是在處理復雜曲邊界時，結合更為準確的邊界處理算法和多塊網格耦合算法，這種優勢更加明顯。

多個圓柱的繞流問題在工程中是很常見的，例如海洋平臺支撐柱、橋墩和水底管路與流體之間的作用以及它們之間的相互作用。圓柱之間相互作用的存在使得此類問題比單圓柱繞流的情況要復雜得多。這類問題的早期研究多側重于對實驗現象的觀察和圓柱力學參數的測量。近代以來數值模擬方法得到了巨大發展，此類問題可以用數值模擬的方法進行研究，并取得了一定研究成果。

Zdravkovich[4-5]等對不同排列方式的雙圓柱繞流問題進行了深入的實驗研究。針對串列圓柱，研究發現了臨界間距的存在。當兩圓柱間距小于該臨界值時，上游圓柱沒有明顯的脫渦現象。這一臨界值約為3.5倍直徑。上下游圓柱尾流場速度分布在該臨界間距附近突然發生改變，上游圓柱出現渦脫落現象，兩柱之間速度和尾流速度都變大。隨著流體力學理論的完善和計算機技術的發展，計算流體力學成為解決流動問題的重要工具。對串列圓柱繞流的數值模擬發現了與實驗現象一致的結果，同時獲得了更多的流場細節，模擬結果也更接近實驗結果[6]。

文章將基于格子Boltzmann方法，對流體的流動和流固耦合作用進行數值模擬。對邊界的處理，采用經典的粒子反彈格式，即流體粒子在碰到固體邊界后速度反向。針對格子Boltzmann方法的均勻規則網格在處理曲邊界問題時的缺陷，采用更為精確的邊界處理算法，讓流體粒子與固體的碰撞反彈過程在固體格點處發生。為提高計算效率，滿足固體周圍相對較高的網格密度要求，開發分塊網格耦合算法。以單個靜止圓柱繞流的已有研究結果為依據，對本文的數值算法進行驗證。最后，對不同間距下靜止串列雙圓柱繞流問題進行數值模擬，研究不同情況下圓柱之間的相互作用以及尾流特征的變化。

1 數值方法

1.1 格子Boltzmann方法

在格子Boltzmann方法中，粒子的運動由波爾茲曼方程描述

其中：u 為流速速度，c=δx/δt，ωα是加權系數，

最終，可以得到粒子的動力學演化方程：

（1）碰撞過程：

（2）遷移過程：

粒子的運動分解為碰撞和遷移兩個過程。在遷移過程中，粒子沿eα方向運動，與相鄰網格點進行數據交換；在碰撞過程中，粒子與來自不同方向的粒子進行碰撞，在滿足質量守恒、動量守恒和能量守恒的條件下，重新構造粒子分布函數。通過（5）式可以看出，粒子的碰撞方程是顯式可解的，遷移過程的計算量也較小，算法的編制相對簡單。

在粒子的分布函數確定后，流場的宏觀量可以表示為

本文采用經典二維九速度模型（D2Q9），將粒子速度離散為9個方向，即α的值有9個，如圖1所示。不同方向的速度矢量表達式為

圖1 二維9速度（D2Q9）模型Fig.1 The 2D,nine-velocity lattice(D2Q9)model

1.2 多塊網格耦合算法

不同于傳統的貼體網格處理方法，格子Boltzmann通常采用均勻結構網格來劃分流場，網格算法相對更容易實現。然而，在處理流固耦合問題時，固體壁面周圍對網格的密度要求較高，而遠離壁面處這種要求相對較低，均勻網格顯然會造成計算量的浪費。因此，研究多塊網格的耦合算法是必要的。本文在保證質量連續和應力連續前提下，采用了多塊網格耦合的方法來處理流固耦合中的網格問題。

為簡單起見，文章以兩塊網格的耦合為例進行分析。粗細網格的尺寸比m，而且在網格塊的計算中滿足 δx=δy=δt，粗細網格交界處的網格結構如圖2所示。

前面已經提到，運動粘度和無量綱松弛時間應滿足（8）式。根據（8）式，為保證在交界處運動粘度的連續，不同網格塊的無量綱松弛時間應該滿足（9）式。

圖2 不同網格塊交界處網格結構Fig.2 Interface stucture between two blocks of different spacing

其中：下標f表示細網格，c表示粗網格。可以將速度分布函數分解為平衡部分和非平衡部分因為在網格交界面上必須滿足速度和密度的連續，因此

圖3 網格耦合算法計算流程圖Fig.3 Flow chart of the computational procedure in the multi-block method

文章采用的多塊網格尺寸比m=2，相鄰網格會有重合。粗網格和細網格有部分節點是重合的，重合節點處可以利用（18）式和（19）式直接進行數據的交換，其他未知點采用多點插值的方法得出，不同的插值方式對計算結果將產生不同的影響。計算流程如圖3所示，n表示當前的時間步，箭頭后表示的是該步驟后得到的計算參數。

1.3 曲邊界處理方法和作用力的提取

在格子Boltzmann方法中，邊界條件的處理方法將直接關系到計算的精度和穩定度。對于靜止的固體邊界，最常見的處理方法是切向速度為零的無滑移邊界條件。具體到格子Boltzmann方法中，入射粒子在碰到邊界后，速度方向變為相反方向，這種處理方法也被稱為反彈格式。反彈格式在處理靜止邊界問題上可以得到較好的效果。但是，由于格子Boltzmann方法多采用均勻結構網格，邊界處會出現鋸齒，只能采用加密網格的方法來提高對邊界的擬合精度，因此無法準確擬合曲邊界。郭照立等[8-9]結合平直邊界的非平衡態外推格式，得到了一種曲邊界處理格式，這種格式通過在固體格點位置處執行碰撞過程來確定碰撞后的分布函數，可以獲得較好的效果。

曲邊界處網格結構如圖4所示。以s點代表固體格點，f代表流體格點，w代表固體邊界點。本文所采用的方法假定碰撞在s處執行，分析以一個方向的處理為例，其他方向的處理方法類似。在固體點s處，將碰撞前的分布函數分解為平衡態feq和非平衡態fneq兩部分

圖4 曲邊界示意圖Fig.4 Schematic description of curved wall boundary

其中，平衡態部分采用粒子的平衡分布函數近似

對于（21）式中的非平衡部分，可表示為

至此，碰撞前固體點的分布函數就確定了。然后，根據LBGK碰撞模型進行標準碰撞，可以得到碰撞后的分布函數

在處理流固耦合問題中，固體的受力參數往往是問題的目標參數，因此流體與固體相互作用力的提取是必要的。本文采用動量交換法來計算流體與圓柱之間的相互作用力，即通過計算流體粒子與壁面發生碰撞前后的動量變化量，根據動量定理來計算粒子所受的流體作用力。在格子Boltzmann方法中，在與壁面碰撞前粒子的動量可以表示為

與壁面碰撞后粒子速度反向，其動量可以表示為

其中的負號表示與α方向的相反方向。

由此，可以計算出流體粒子所受的作用力。根據經典牛頓力學中相互作用力的性質，二者大小相

同而方向相反。因此，通過作用力沿壁面的積分就可以得到固體所受的作用力，具體形式為

進而可以得出圓柱的阻力系數Cd和升力系數Cl：

其中：FD和FL分別是圓柱所受的流體作用力在順流方向和橫流方向的分量，U為來流速度。

2 數值方法驗證

為驗證算法的可靠性，本文在耦合網格下對Re=200靜止單圓柱繞流問題進行了數值模擬，其中雷諾數Re=DU/v，D為圓柱直徑，U為來流速度，v是流體的運動粘度。計算區域的劃分和網格劃分情況如圖5所示。在網格密度要求較高的圓柱周圍采用尺寸較小的網格塊，其他區域采用尺寸相對較大的網格塊。 細網格區域取即一倍直徑D分為50個格子，粗網格區域網格尺寸是細網格區域的2倍，即網格尺寸比m=2。網格塊尺寸間距和時間步大小的確定兼顧了計算的精度和效率。靜止單個圓柱繞流問題已經有較多的理論研究、實驗測量和數值模擬研究成果，通過與已有結果對比可以驗證本文算法的可靠性。

本文數值模擬結果如圖6所示。圓柱表面發生了周期性的脫渦現象，正負渦從圓柱上下表面交替脫落。由圓柱升阻力系數曲線可以看出，其升阻力系數曲線呈現周期性震蕩，Cd的頻率是Cl頻率的兩倍，這是因為每個周期脫落了一個正渦和一個負渦。

圖5 計算區域的劃分Fig.5 The division of computational domain

圖6 單網格下圓柱繞流計算結果Fig.6 The results of single cylinder with single grid block

表1是本文結果與文獻研究成果的比較，可以看出本文結果和文獻結果符合得較好。其中的Strouhal數是表征圓柱脫渦情況的重要參數，定義為

其中：f為圓柱的脫渦頻率，U為自由來流速度。

表1 結果對比Tab.1 Comparison of results

3 串列雙圓柱相互影響研究

對于靜止單個圓柱的繞流問題，流場特征和圓柱受力情況主要依賴于雷諾數Re。對于串列雙圓柱繞流問題，除了雷諾數之外，兩圓柱中心間距比L/D的影響也是關鍵的（L為前后圓柱中心間距，D為圓柱直徑）。

3.1 計算區域和網格劃分

本文在耦合網格下對Re=200串列雙圓柱繞流問題進行了數值模擬，計算區域的劃分如圖7所示。計算區域分為兩部分，兩個圓柱周圍采用尺寸較小的細網格，離圓柱較遠的區域采用尺寸較大的粗網格。粗網格區域取粗細網格尺寸比 m=2。

圖7 計算區域的劃分Fig.7 The division of computational domain

3.2 計算結果分析

本文對圓柱中心間距比L/D=1.5、2.0、3.0和4.0的情況進行了數值模擬。不同間距比下的上下游圓柱升力系數曲線、阻力系數曲線、流場渦量等值線和流場流線如圖8～11所示。

圖8 L=1.5D情況下雙圓柱計算結果Fig.8 The results of two cylinders with L=1.5D

圖11 L=4.0D情況下雙圓柱計算結果Fig.11 The results of two cylinders with L=4.0D

本文將不同間距下上下游圓柱升力系數、阻力系數和Strouhal數進行了匯總，并和文獻結果進行了對比，如表2所示。這三個參數是表征圓柱受力情況的重要參數，可以作為判斷數值模擬結果可靠性的依據。可以看出，本文計算結果與文獻計算值基本相符。

表2 結果對比Tab.2 Comparison of results

綜合數值模擬結果可以看出臨界間距的存在（3.0＜L/D＜4.0）。當間距小于該臨界值時，上游圓柱的分離剪切層附著在下游圓柱上，只有下游圓柱有脫渦現象。此時，下游圓柱的阻力系數較小且小于零。當兩圓柱間距超過該臨界值時，上下游圓柱都產生明顯的渦脫落，且升力系數最大幅值都變大，平均阻力系數也突然變大，特別是對下游圓柱而言，由一個負值變到一個較大的正值。

具體來看，當L/D=1.5、2.0（小于臨界值）時，下游圓柱的阻力系數為負值，即下游圓柱并沒有受到與來流方向相同的阻力，而是受到了與來流方向相反的推力。L/D=1.5、2.0時，兩柱的升力系數幅值都比較小，尤其上游圓柱尾流受到限制，其升力系數幾乎為零。當間距比增大到4.0（大于臨界值）時，兩圓柱阻力系數和升力系數都突然增大。上游柱的升力系數振幅突然增大至0.8左右，下游圓柱的升力系數振幅增至更大的值。下游圓柱的阻力系數由負值變為較大的正值，但其依然比上游柱的阻力系數小，每個圓柱的阻力系數也都小于單柱繞流時的阻力系數。由升阻力系數變化曲線還可以看出，升力振動正負幅值相等，即圓柱所受升力合力為零。另外，上、下游圓柱的渦脫落頻率保持相等，但都較單圓柱繞流小。從Strouhal數的變化趨勢可以看出，當間距比跨越臨界間距值時Strouhal數會突然增大，而隨著兩圓柱間距的繼續增大，Strouhal數接近于單一圓柱繞流情況。

4 結 論

本文基于格子Boltzmann方法，結合多塊網格耦合算法和曲邊界處理算法，對雷諾數Re=200條件下的靜止串列雙圓柱繞流進行了數值模擬。數值模擬得到了和已有研究成果一致的結果，驗證了臨界間距的存在。當兩圓柱間距小于該臨界值時，上游圓柱沒有明顯脫渦現象，下游圓柱出現周期性脫渦；當間距大于該臨界值時，上游圓柱開始出現周期性脫渦，而下游圓柱所受的升、阻力系數明顯提高，并且下游圓柱阻力方向發生了突變。不同間距下兩圓柱受力參數和脫渦情況與已有研究成果符合得較好。

[1]Chen S,Doolen G D.Lattice Boltzmann method for fluid flows[J].Annual Review of Fluid Mechanics,1983,30:329-364.

[2]韓文驥,顏 開,褚學森,張 珂.LBM方法數值模擬飛濺流動的影響因素研究[J].船舶力學,2015,19(4):356-362.Han Wenji,Yan Kai,Chu Xueseng,Zhang Ke.Influence factors to numerical simulation of splash flow based on LBM[J].Journal of Ship Mechanics,2015,19(4):356-362.

[3]黃橋高,潘 光.疏水表面流體流動特性的格子Boltzmann方法模擬[J].船舶力學,2016,20(10):1211-1218.Huang Qiaogao,Pan Guang.Lattice Boltzmann simulation of liquid flow characteristics of hydrophobic surfaces[J].Journal of Ship Mechanics,2016,20(10):1211-1218.

[4]Zdravkovich M M.Review of flow interference between two circular cylinders in various arrangements[J].ASME Journal of Fluids Engineering,1977,99:618-633.

[5]Zdrvakoeieh M M.The effets of interference between circular cylinders in a cross flow[J].Journal of Fluids and Structures,1987,l:239-261.

[6]梁亮文.低雷諾數下圓柱橫向受迫振蕩和渦激運動的數值分析[D].上海:上海交通大學,2009.Liang Liangwen.Numerical analysis of forced oscillation and vortex-induced motion of circular cylinder in cross flow with low Reynolds number[D].Shanghai:Shanghai Jiao Tong University,2009.

[7]Yu D,Mei R,Shyy W.A multi-block lattice Boltzmann method for viscous fluid flows[J].International Journal for Numerical Methods in Fluids,2002,39:99-120.

[8]Guo Z L,Zheng C G,Shi B C.An extrapolation method for boundary conditions in lattice Boltzmann method[J].Physics of Fluids,2002,14:2007-2010.

[9]龔 帥,郭照立.流向振蕩圓柱繞流的格子Boltzmann方法模擬[J].力學學報,2011,43(1):11-17.Gong Shuai,Guo Zhaoli.Lattice Boltzmann simulation of the flow around a circular cylinder oscillating streamwisely[J].Chinese Journal of Theoretical and Applied Mechanics,2011,43(1):11-17.

[10]Chan C T,Anastasiou K.Solution of incompressible flows with or without a free surface using the finite volume method on unstructured triangular meshes[J].International Journal for Numerical Methods in Fluids,1995,29:35-57.

[11]Lecointe Y,Piquet J.On the use of several compact methods for the study of unsteady incompressible viscous flow round a circular cylinder[J].Computers&Fluids,1984,12(4):255-280.

[12]Rosenfeld M,Kwak D,Vinokur M.A solution method for the unsteady and in compressible Navier-Stokes equations in generalized coordinate systems[R].AIAA Paper 88-0718,1988.

[13]劉 松,符 松.串列雙圓柱繞流問題的數值模擬[J].計算力學學報,2000,17(3):260-266.Liu Song,Fu Song.Numerical simulation of flow past two cylinders in tandem arrangement[J].Chinese Journal of Computational Mechanics,2000,17(3):260-266.