靳遵龍， 王元凱， 王永慶

(鄭州大學(xué) 化工與能源學(xué)院 河南 鄭州 450001)

0 引言

實際生活中很多情形都可以歸結(jié)為鈍體繞流問題，如河流繞過橋墩的流動以及建筑物的風(fēng)載等.流體流經(jīng)鈍體時，由于過流界面減小，會產(chǎn)生一系列復(fù)雜的物理現(xiàn)象，特別是邊界層的分離進而在鈍體后產(chǎn)生周期脫落的渦.目前鈍體繞流研究主要以圓柱繞流為主，對于一般形狀柱體的研究還更多的處于實驗階段.

對于復(fù)雜流體繞流，傳統(tǒng)的CFD方法是從宏觀Navier-Stokes方程出發(fā)，利用各種數(shù)值方法進行模擬.文獻[1]分別使用有限差分法和離散渦方法計算了低雷諾數(shù)和高雷諾數(shù)下的方柱繞流；文獻[2]通過有限元法求解了2個圓柱左右并排時和前后排列時的繞流.格子Boltzmann方法作為一種介觀模擬方法，將宏觀運動和微觀運動的統(tǒng)計平均聯(lián)系起來，在流體力學(xué)問題的數(shù)值模擬方面具有廣闊的前景.格子Boltzmann方法不但可以方便地模擬復(fù)雜幾何形狀邊界的流體流動，還可以模擬系統(tǒng)的時間演化，并且只需要簡單的算法就可以在計算機上實現(xiàn)并行計算，具有更高的計算效率.本文利用格子Boltzmann方法對不同雷諾數(shù)下固定單方柱繞流流場進行分析，并對并列雙方柱不同分布間距的流場進行了模擬，驗證了格子Boltzmann方法邊界處理和數(shù)值模擬的正確性和便捷性.

1 格子Boltzmann方法

圖1 D2Q9模型離散速度Fig.1 Discrete velocities of D2Q9 model

格子Boltzmann方法中的基本變量是格點的分布函數(shù)fi，fi為沿i方向的粒子分布函數(shù).本文采用D2Q9模型，9個離散速度分布如圖1所示.其演化方程(不含外力項)為

式中：Δt為時間步長；τ為無量綱松弛時間，τ=τ0/Δt；離散速度ci表示粒子的運動方向，可以表示為

c0=(0,0)，c1,3,c2,4=(±c,0),(0,±c)，c5,6,c7,8=(±c,c),(c,±c)，

由平衡態(tài)分布函數(shù)通過Chapman Enskog展開，可以導(dǎo)出模型的宏觀密度和宏觀速度分別為

格子Boltzmann方法在物理空間上將系統(tǒng)粒子運動分為遷移和碰撞兩個相對獨立的過程，使其具備很好的并行特性以及較強的復(fù)雜邊界處理能力.

2 格子Boltzmann方法的初始條件和邊界格式處理

流體的運動總是在一定的初始條件和邊界條件下進行的，因此必須設(shè)置相應(yīng)的初始條件和邊界條件.格子Boltzmann方法中的基本變量是分布函數(shù).對于穩(wěn)態(tài)或者準穩(wěn)態(tài)流動，初始條件對最終計算結(jié)果影響不大，可以直接將初始分布函數(shù)設(shè)為其平衡態(tài)分布函數(shù).實際問題中，邊界條件往往基于宏觀物理量，如何根據(jù)宏觀量合理地確定分布函數(shù)是格子Boltzmann方法中的重要問題.常用的邊界主要分為啟發(fā)格式邊界、動力學(xué)格式邊界、外推格式邊界以及復(fù)雜格式邊界；根據(jù)實際問題又可分為速度邊界和壓力邊界等.

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

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

實際應(yīng)用中，邊界處的速度分量或者壓力分布往往是已知的，如管道流動等.文獻[4]在1997年提出的非平衡態(tài)反彈格式則可以很好地應(yīng)用于上述兩種邊界條件. 對于速度邊界條件，有

通過非平衡態(tài)反彈格式可以直接計算垂直于速度邊界上的3個未知平衡態(tài)分布函數(shù)，相似的方法也可以用于壓力邊界和其他邊界條件.

3 固定方柱繞流的數(shù)值模擬與分析

3.1 單方柱繞流模擬

圖2 單方柱繞流模型示意圖Fig.2 Schematic diagram of channel flow past a single square cylinder

圖3為不同雷諾數(shù)下流場渦量云圖.可以看出，當(dāng)Re較小時，流場是定常的，方柱后面出現(xiàn)一對上下對稱的尾渦.隨著Re的增加，當(dāng)Re=38時，流場逐漸轉(zhuǎn)變?yōu)榉嵌ǔ?，尾渦也逐漸脫落，尾流區(qū)中形成周期性擺動和交錯的漩渦，即出現(xiàn)卡門渦街.隨著Re的進一步增加，當(dāng)Re=238時，尾渦有紊亂的趨勢.進一步的模擬顯示，當(dāng)Re大于238時，尾流區(qū)由層流狀態(tài)過渡到紊流狀態(tài).

圖3 不同雷諾數(shù)下流場渦量云圖Fig.3 Vorticity contours at different Reynolds numbers

圖4顯示了Re=200時阻力系數(shù)和升力系數(shù)隨時間的變化曲線，當(dāng)前平均阻力系數(shù)和最大升力系數(shù)分別為1.51和0.31，與文獻[5]中實驗研究的數(shù)據(jù)相比具有一致性，可見利用格子Boltzmann方法模擬是可行的.

圖4 Re=200時阻力系數(shù)和升力系數(shù)隨時間的變化曲線Fig.4 Time history of drag coefficient and lift coefficient at Re=200

Re=35及Re=37時的流場流線分析如圖5所示，可以看出，上下對稱的尾渦逐漸擴大，上渦有分離的趨勢.因此，定常流失穩(wěn)的臨界Re在37左右，這與文獻[6]中的結(jié)論是相符合的.

圖5 不同雷諾數(shù)下流場流線圖Fig.5 Streamline patterns at different Reynolds numbers

與理論研究及文獻中圓柱繞流的結(jié)果對比后發(fā)現(xiàn)，方柱繞流整體上符合一般的鈍體繞流趨勢，并且由于方柱的特殊性，臨界Re又與圓柱繞流不同.由以上分析可知，方柱繞流中卡門渦街現(xiàn)象在Re為40～238時較為明顯.本文將采用Re=200進行下一步的模擬研究.

3.2 雙方柱繞流模擬

圖6 并列雙方柱繞流模型示意圖Fig.6 Schematic diagram of channel flow past parallel square cylinders

實際生活中更多的是流體繞過多個柱體流動的情形，柱體的個數(shù)、形狀、排列方式等均會影響流場的結(jié)構(gòu).本文研究雙方柱并列模型下不同分布間距對流場的影響，控制的變量為柱間距M與方柱特征尺寸H的比值(M/H)，并列雙方柱繞流模型示意圖如圖6所示.

采用相同的邊界條件，X、Y方向格子數(shù)分別為360和120，在Re=200的條件下分別對并列方柱M/H=1、M/H=2、M/H=3、M/H=4四種條件進行模擬，其渦量云圖如圖7所示.

圖8顯示了并列方柱在M/H=2和M/H=3 時的流場流線圖.從圖7和圖8中可以看出，當(dāng)分布間距較小(M/H=1)時，兩柱形成的漩渦彼此影響較小，還保留著部分單柱繞流的流場特征；隨著分布間距的增加，流場逐漸變得復(fù)雜，漩渦分布呈現(xiàn)不規(guī)則狀態(tài)，當(dāng)M/H=2時，流體繞流方柱的渦流彼此影響最為明顯；分布間距繼續(xù)增加，當(dāng)M/H大于3時，兩個渦的相互影響越來越小，方柱下游形成兩個反向同步脫落的渦街，各自渦形逐漸接近單方柱卡門渦街狀態(tài).

圖7 不同分布間距下流場渦量云圖Fig.7 Vorticity contours at different distribution distances

圖8 不同分布間距下流場流線圖Fig.8 Streamline patterns at different distribution distances

4 小結(jié)

通過對單方柱和并列雙方柱的模擬研究，驗證了格子Boltzmann方法邊界處理和數(shù)值模擬的正確性和便捷性.通過對不同Re下單方柱流場的模擬分析，探究了方柱產(chǎn)生卡門渦街的臨界雷諾數(shù)(Re=37)，并得到了產(chǎn)生明顯卡門渦街的雷諾數(shù)范圍(Re=40～238).此外，雷諾數(shù)Re=200下并列雙方柱不同分布間距(M/H=1、M/H=2、M/H=3、M/H=4)的流場模擬結(jié)果表明，當(dāng)柱間距為2倍方柱邊長時，流體繞流方柱的渦流彼此影響最為明顯.