謝晨 李銀山 霍樹浩

摘要 基于格子Boltzmann方法的原理和對邊界條件處理的便捷性，分別對單方柱和并列雙方柱繞流進行了研究。對于單方柱繞流，分別研究了不同雷諾數(shù)和阻塞比對它的影響，給出了不同雷諾數(shù)下單方柱繞流的流線圖和渦量圖，以及平均阻力系數(shù)和最大升力系數(shù)隨雷諾數(shù)的變化曲線，得到單方柱產(chǎn)生卡門渦街的臨界雷諾數(shù)，也給出了不同阻塞比下單方柱繞流的渦量圖，以及不同阻塞比對平均阻力系數(shù)和最大升力系數(shù)所產(chǎn)生的影響。對于并列雙方柱繞流，研究了不同間距比對它的影響，給出了不同間距比下的并列雙方柱繞流的流線圖和渦量圖，得到了流動從偏流型向?qū)ΨQ型轉(zhuǎn)換的臨界間距比以及使得兩個方柱后的渦流互相影響最為明顯的間距比。

關(guān) 鍵 詞 格子Boltzmann方法;方柱繞流;雷諾數(shù); 阻塞比;間距比

Abstract Based on the lattice Boltzmann method and the convenience of the boundary condition processing， the flow past a single square cylinder and the parallel square cylinder were studied respectively. For the single square cylinder flow， the effects of different Rayleigh numbers and blocking ratios on it are studied separately， the streamline diagram and vorticity diagram of the flow around a single square cylinder under different Rayleigh numbers， and the variation curve of the average drag coefficient and the maximum lift coefficient with the Rayleigh number are given. The critical Rayleigh number of the Karman vortex street generated by the single square cylinder is obtained. The vorticity diagrams of the flow around a single square cylinder with different blocking ratios and the effects of different blocking ratios on the average drag coefficient and the maximum lift coefficient are also given. For the parallel flow around the two square cylinder， the influence of different spacing ratios on the flow is studied. The flow diagram and vorticity diagram of the flow around the two square cylinder at different spacing ratios are given， and the flow from the bias flow to the symmetrical conversion is obtained. The eddy currents behind the two square cylinder affect each other's most obvious spacing ratio.

Key words lattice boltzmann method; flow past square cylinder; reynolds number;blocking ratio; spacing ratio

0 引言

鈍體繞流廣泛存在于實際生活中，例如流體經(jīng)過橋墩或者海洋平臺等。柱體繞流是鈍體繞流的經(jīng)典問題之一，在流體經(jīng)過柱體時，由于柱體的阻礙作用，在柱體后面會產(chǎn)生渦，而柱體后面的渦會對柱體產(chǎn)生很大的破壞力，當流體的速度達到某個臨界值時，在柱體后面就會出現(xiàn)卡門渦街。由于這類問題在工程上的重要性，所以越來越多的人開始研究柱體繞流。

許多人用數(shù)值模擬或者試驗的方法來研究不同的柱體繞流問題，其中數(shù)值模擬不失為一種有效的手段，如周強等[1]用LES方法對單方柱繞流進行了模擬;Okajima A[2]用有限差分法對單方柱繞流進行了模擬;程維賢等[3]用有限元模擬單方柱和串列雙方柱繞流;陳素琴等[4]用 MAC方法對并列雙方柱繞流進行了模擬。格子Boltzmann方法的思路和傳統(tǒng)的模擬流體的方法完全不同，它用粒子群來代替流體，粒子流沿著給定的方向遷移并在格子節(jié)點處發(fā)生碰撞重整，再對每個粒子新得到的分布函數(shù)進行統(tǒng)計，得到密度和速度等宏觀變量，求得的宏觀物理量能夠滿足宏觀偏微分方程[5-6]，對于復雜邊界條件[7-12]的問題處理起來相對容易，只需要簡單的算法就可以在計算機上實現(xiàn)并行[13]計算，具有更高的計算效率。本文對單方柱和并列雙方柱繞流進行了模擬，得到了不同雷諾數(shù)和阻塞比對單方柱繞流的影響以及兩方柱之間的間距比對并列雙方柱繞流造成的影響。

1 格子Boltzmann方法

格子Boltamann方法是由McNamara和Zaretti在1988年引入，用于克服格子氣自動機的缺點，從此，格子Boltzmann方法稱為一種求解流體動力學問題的功能強大的替代方法。格子Boltzmann方法可視為一種顯示方法，編程容易，并行性也很好。

2 格子Boltzmann方法的邊界格式處理

流體運動總會涉及到邊界條件，格子Boltzmann方法邊界處理格式有很多種，本文根據(jù)問題所給邊界用到非平衡態(tài)反彈格式和非平衡態(tài)外推格式這2種邊界處理方式。

2.1 非平衡態(tài)外推格式邊界

非平衡態(tài)外推格式邊界的核心是將邊界節(jié)點a上的分布函數(shù)分成平衡態(tài)和非平衡態(tài)兩部分，其中平衡態(tài)部分根據(jù)邊界條件的定義近似得到，非平衡態(tài)部分則用流體相鄰節(jié)點b處的非平衡態(tài)部分代替。表達式如下：

2.2 非平衡態(tài)反彈格式邊界

3 方柱繞流模擬

3.1 單方柱繞流模擬

3.1.1 物理模型及邊界條件

從圖2中可以看出計算區(qū)域為45D×15D，其中D為方柱邊長，入口速度為均勻來流U = 0. 1。其中方柱的特征尺寸D取8個格子單位，管道水平和豎直方向分別劃分為360和120個格子單位。入口和出口分別采用速度和壓力邊界，其他邊界用的都是非平衡態(tài)外推格式邊界。

3.1.2 利用格子Boltzmann方法估計柱體表面的力

利用格子Boltzmann方法估計柱體表面的力通常有2種辦法，應(yīng)力積分法[15]和動量交換法[16]。相比之下，動量交換法實現(xiàn)起來更加簡單，直接用分布函數(shù)進行計算，特別適合于用在格子Boltzmann方法中，這也是格子Boltzmann方法區(qū)別于其他方法的特色之一。所以，本文用動量交換法來計算柱體表面的升力系數(shù)和阻力系數(shù)。

3.1.3 雷諾數(shù)對方柱繞流的影響

本文對雷諾數(shù)Re為10～300范圍內(nèi)的方柱繞流進行模擬，雷諾數(shù)定義為[Re=UDv]，其中[ν]為運動粘性系數(shù)。圖3給出了不同雷諾數(shù)時方柱繞流的流線圖，從圖3可以看出，當Re = 10～89時，在方柱后緣發(fā)生流動分離現(xiàn)象，方柱后面產(chǎn)生一對對稱的回流漩渦，隨著雷諾數(shù)的增大漩渦也在逐漸擴大;當Re = 90時，尾流的對稱狀態(tài)破壞，開始出現(xiàn)漩渦脫落現(xiàn)象，因此，定常流失穩(wěn)的臨界Re在90左右。

圖4為不同雷諾數(shù)下流場渦量云圖。從圖4可以看出，當Re比較小的時候，在方柱后面形成上下對稱的渦，并且在一定范圍內(nèi)隨著Re的增加渦在變長，此時流場是定常的。隨著Re的持續(xù)增加，方柱后的渦開始脫落，當Re = 90時，方柱后面出現(xiàn)周期性擺動和交錯的卡門渦街，此時流場變?yōu)榉嵌ǔ?。當Re增加到300時，尾渦有了紊亂的趨勢。

Re對平均阻力系數(shù)[Cdavg]和最大升力系數(shù)[Clmax]的影響如圖5所示。當Re = 10～100時，平均阻力系數(shù)下降比較快，說明雷諾數(shù)在這個范圍內(nèi)的變化對平均阻力系數(shù)的影響比較大。此后，隨著Re的增大，平均阻力系數(shù)下降的變的相對平緩，Re對平均阻力系數(shù)的影響變小。當Re = 10～89時，最大升力系數(shù)為0，由于此時雷諾數(shù)比較小，對最大升力系數(shù)的影響還不足以顯現(xiàn)出來。此后，Re的變化對最大升力系數(shù)的影響慢慢凸顯出來，即隨著Re的增大，最大升力系數(shù)在逐漸變大。

3.1.4 阻塞比的影響

在研究阻塞比[B=DH]對方柱繞流的影響時，本文取Re = 200，除了改變阻塞比之外，其他條件不變。圖7給出了阻塞比B為0.083、0.125、0.167和0.25的方柱繞流的渦量圖。

由圖7可見，當阻塞比較小時，如當B = 0.083時，此時壁面效應(yīng)對渦脫落所產(chǎn)生的影響比較小，方柱后面的卡門渦街比較清楚、規(guī)則。隨著阻塞比的增大，上下壁面附近開始出現(xiàn)小渦，這些小渦與方柱后面產(chǎn)生的大渦相混合，一同向下游運動;這種現(xiàn)象隨阻塞比的增加越來越明顯（圖7c））;隨著阻塞比進一步增加（圖7d）），在方柱后形成了很長的回流區(qū)，渦脫落變得不再規(guī)則，剛剛脫落的大渦會立刻與壁面附近的小渦混合，尾跡里見不到明顯的卡門渦街。

圖8給出了方柱平均阻力系數(shù)和最大升力系數(shù)隨阻塞比的變化曲線。由圖可知，平均阻力系數(shù)隨阻塞比增大逐漸增加，最大升力系數(shù)則隨阻塞比增大在逐漸減小。這是因為隨著阻塞比的增大，上下邊界與方柱越來越靠近，流體的流動受到影響，這就導致了柱體表面的阻力系數(shù)和升力系數(shù)發(fā)生改變。

3.2 雙方柱繞流模擬

本文研究的是并列雙方柱之間的間距S與方柱特征尺寸D的比值（間距比）對流場的影響。計算示意圖如圖9所示。

計算區(qū)域為60D × 18D，其中方柱的特征尺寸D取10個格子單位，管道水平和豎直方向分別劃分為600和180個格子單位。與單方柱采用相同的邊界條件，在Re = 200的條件下分別對并列方柱間距比S/D=0.2、S/D = 1.0、S/D = 1.4、S/D = 1.8、S/D = 2.0、S/D = 4.0這6種條件進行模擬。

圖10給出了6種不同方柱間距比下的流場流線圖。從圖10中可以看出，剛開始方柱間距很小，間隙流不足以影響方柱后方的流動狀態(tài)，此時很難看到偏流現(xiàn)象，所以在方柱后方會形成共同的漩渦（圖10a））。隨著2個方柱之間的距離慢慢變大，2個方柱彼此之間的干擾在減弱，但是，間隙流對流場的影響開始顯現(xiàn)出來，每個方柱后面都產(chǎn)生了各自的漩渦，并且把它們共同產(chǎn)生的漩渦緩慢的推離方柱后方（圖10b）和圖10c）），偏流現(xiàn)象已經(jīng)出現(xiàn)。當間距比達到臨界值S/D =1.8時，偏流現(xiàn)象消失（圖10d）），間距比S/D >1.8時，2個方柱后面出現(xiàn)了對稱的尾流（圖10e）和圖10f）），此時流場轉(zhuǎn)變?yōu)閷ΨQ流。

圖11給出了6種不同方柱間距比下的流場渦量圖。從圖11中可以看出，當間距比較小S/D = 0.2時，2個方柱之間的影響較小，流場的特征和單方柱繞流相似;隨著間距比的增加，流場開始變復雜，2個方柱后的漩渦也變的不規(guī)則;當間距比S/D = 1.4時，2個方柱的渦流彼此影響最明顯;當間距比S/D > 1.8時，2個方柱各自形成的渦之間的影響越來越小，方柱后面形成2個反向同步脫落的渦街，每個方柱后面的渦形逐漸接近單方柱卡門渦街狀態(tài)。

4 結(jié)語

本文利用格子Boltzman方法對單方柱和并列雙方柱繞流進行了數(shù)值模擬。對于單方柱研究了雷諾數(shù)Re = 10～300范圍內(nèi)以及4種阻塞比即B = 0.083、0.125、0.167和0.25對單方柱繞流的影響，給出了不同雷諾數(shù)下單方柱繞流的流線圖和渦量圖，得到定常流失穩(wěn)的臨界Re和開始出現(xiàn)明顯卡門渦街的Re在90左右，也得到了當Re = 10～300時平均阻力系數(shù)隨Re的增大而減小，當Re = 10～89時，最大升力系數(shù)為0，此后，隨著Re的增大，最大升力系數(shù)在逐漸變大。當阻塞比較小時，壁面效應(yīng)對渦脫落所產(chǎn)生的影響比較小，方柱后面的卡門渦街比較清楚、規(guī)則，隨著阻塞比的增大，上下壁面附近開始出現(xiàn)小渦，這些小渦與方柱后面產(chǎn)生的大渦相混合，一同向下游運動，這種現(xiàn)象隨阻塞比的增加越來越明顯，隨著阻塞比的進一步增加，在方柱后形成了很長的回流區(qū)，渦脫落變得不再規(guī)則，剛剛脫落的大渦會立刻與壁面附近的小渦混合，尾跡里見不到明顯的卡門渦街。也得到了平均阻力系數(shù)和最大升力系數(shù)都隨阻塞比增大逐漸增加。對于并列雙方柱研究了2個方柱之間的間距比對并列雙方柱繞流的影響，給出了不同間距比下的流場流線圖和渦量圖，得到了間距比S/D = 1.8為流場由偏流轉(zhuǎn)化為對稱流的臨界值，也得到了使得雙方柱繞流的渦流彼此影響最為明顯的間距比為S/D = 1.4。

參考文獻：

[1]? ? 周強，廖海黎，曹曙陽. 高雷諾數(shù)下方柱繞流特性的數(shù)值模擬[J]. 西南交通大學學報，2018，53（3）：533-539.

[2]? ? OKAJIMA A. Numerical simulation of flow around rectangular cylinders[J]. Journal of Wind Engineering and Industrial Aerodynamics，1990，33（1/2）：171-180.

[3]? ? 穆維賢，賴國璋. 方柱和雙方柱繞流的數(shù)值模擬[J]. 水動力學研究與進展（A輯），1990，5（4）：76-80.

[4]? ? 陳素琴，顧明，黃自萍. 兩并列方柱繞流相互干擾的數(shù)值研究[J]. 應(yīng)用數(shù)學和力學，2000，21（2）：131-146

[5]? ? HE X Y，LUO L S. Lattice boltzmann model for the incompressible navier-stokes equation[J]. Journal of Statistical Physics，1997，88（3/4）：927-944.

[6]? ? MCNAMARA G R，ZANETTI G. Use of the Boltzmann equation to simulate lattice gas automata[J]. Physical Review Letters，1988，61（20）：2332-2335.

[7]? ? 王露，李天勻，朱翔，等. 基于改進的浸沒邊界-格子Boltzmann方法的圓柱繞流仿真計算[J]. 中國造船，2017，58（4）：150-159.

[8]? ? 顧娟，黃榮宗，劉振宇，等. 一種滑移區(qū)氣體流動的格子Boltzmann曲邊界處理新格式[J]. 物理學報，2017，66（11）：220-228.

[9]? ? 任曉飛，魏守水，劉飛飛，等. 改進的浸入邊界-晶格Boltzmann法的蠕動流分析[J]. 振動. 測試與診斷，2018，38（3）：454-460.

[10]? 顧娟，黃榮宗，劉振宇，等. 一種滑移區(qū)氣體流動的格子Boltzmann曲邊界處理新格式[J]. 物理學報，2017，66（11）：220-228.

[11]? 程志林，寧正福，曾彥，等. 格子Boltzmann方法模擬多孔介質(zhì)慣性流的邊界條件改進[J]. 力學學報，2019，51（1）：124-134.

[12]? DE ROSIS A. Ground-induced lift enhancement in a tandem of symmetric flapping wings：Lattice Boltzmann-immersed boundary simulations[J]. Computers & Structures，2015，153：230-238.

[13]? SCHMIESCHEK S，SHAMARDIN L，F(xiàn)RIJTERS S，et al. LB3D：a parallel implementation of the Lattice-Boltzmann method for simulation of interacting amphiphilic fluids[J]. Computer Physics Communications，2017，217：149-161.

[14]? CHEN H D，CHEN S Y，MATTHAEUS W H. Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method[J]. Physical Review A，1992，45（8）：r5339.

[15]? LADD A J C，VERBERG R. Lattice-boltzmann simulations of particle-fluid suspensions[J]. Journal of Statistical Physics，2001，104（5/6）：1191-1251.

[16]? HE X Y，DOOLEN G. Lattice boltzmann method on curvilinear coordinates system：flow around a circular cylinder[J]. Journal of Computational Physics，1997，134（2）：306-315.