基于格子Boltzmann方法的固定方柱繞流研究

靳遵龍， 王元凱， 王永慶

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

0 引言

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

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

1 格子Boltzmann方法

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

格子Boltzmann方法中的基本變量是格點的分布函數fi，fi為沿i方向的粒子分布函數.本文采用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)，

由平衡態分布函數通過Chapman Enskog展開，可以導出模型的宏觀密度和宏觀速度分別為

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

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

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

2.1 非平衡態外推格式邊界

2.2 非平衡態反彈格式邊界

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

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

3 固定方柱繞流的數值模擬與分析

3.1 單方柱繞流模擬

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

圖3為不同雷諾數下流場渦量云圖.可以看出，當Re較小時，流場是定常的，方柱后面出現一對上下對稱的尾渦.隨著Re的增加，當Re=38時，流場逐漸轉變為非定常，尾渦也逐漸脫落，尾流區中形成周期性擺動和交錯的漩渦，即出現卡門渦街.隨著Re的進一步增加，當Re=238時，尾渦有紊亂的趨勢.進一步的模擬顯示，當Re大于238時，尾流區由層流狀態過渡到紊流狀態.

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

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

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

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

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

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

3.2 雙方柱繞流模擬

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

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

采用相同的邊界條件，X、Y方向格子數分別為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中可以看出，當分布間距較小(M/H=1)時，兩柱形成的漩渦彼此影響較小，還保留著部分單柱繞流的流場特征；隨著分布間距的增加，流場逐漸變得復雜，漩渦分布呈現不規則狀態，當M/H=2時，流體繞流方柱的渦流彼此影響最為明顯；分布間距繼續增加，當M/H大于3時，兩個渦的相互影響越來越小，方柱下游形成兩個反向同步脫落的渦街，各自渦形逐漸接近單方柱卡門渦街狀態.

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

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

4 小結

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