謝晨 李銀山 霍樹(shù)浩

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

關(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 引言

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

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

1 格子Boltzmann方法

格子Boltamann方法是由McNamara和Zaretti在1988年引入，用于克服格子氣自動(dòng)機(jī)的缺點(diǎn)，從此，格子Boltzmann方法稱為一種求解流體動(dòng)力學(xué)問(wèn)題的功能強(qiáng)大的替代方法。格子Boltzmann方法可視為一種顯示方法，編程容易，并行性也很好。

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

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

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

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

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

3 方柱繞流模擬

3.1 單方柱繞流模擬

3.1.1 物理模型及邊界條件

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

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

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

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

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

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

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

3.1.4 阻塞比的影響

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

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

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

3.2 雙方柱繞流模擬

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

計(jì)算區(qū)域?yàn)?0D × 18D，其中方柱的特征尺寸D取10個(gè)格子單位，管道水平和豎直方向分別劃分為600和180個(gè)格子單位。與單方柱采用相同的邊界條件，在Re = 200的條件下分別對(duì)并列方柱間距比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種條件進(jìn)行模擬。

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

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

4 結(jié)語(yǔ)

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

參考文獻(xiàn)：

[1]? ? 周強(qiáng)，廖海黎，曹曙陽(yáng). 高雷諾數(shù)下方柱繞流特性的數(shù)值模擬[J]. 西南交通大學(xué)學(xué)報(bào)，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]? ? 穆維賢，賴國(guó)璋. 方柱和雙方柱繞流的數(shù)值模擬[J]. 水動(dòng)力學(xué)研究與進(jìn)展（A輯），1990，5（4）：76-80.

[4]? ? 陳素琴，顧明，黃自萍. 兩并列方柱繞流相互干擾的數(shù)值研究[J]. 應(yīng)用數(shù)學(xué)和力學(xué)，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]? ? 王露，李天勻，朱翔，等. 基于改進(jìn)的浸沒(méi)邊界-格子Boltzmann方法的圓柱繞流仿真計(jì)算[J]. 中國(guó)造船，2017，58（4）：150-159.

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

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

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

[11]? 程志林，寧正福，曾彥，等. 格子Boltzmann方法模擬多孔介質(zhì)慣性流的邊界條件改進(jìn)[J]. 力學(xué)學(xué)報(bào)，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.