格子玻爾茲曼方法(LBM)及其在微通道繞流中的應用

馮俊杰，孫 冰，姜 杰，徐 偉，石 寧

(中國石化青島安全工程研究院化學品安全控制國家重點實驗室，山東青島 266071)

0 前言

微反應器在提高反應過程安全性、縮短反應時間、提高轉化率、靈活生產等方面具有獨特的優勢，實現微通道流動的精確測定和控制是微反應器發揮諸多優勢的保障和廣泛應用的基礎[1]。由于微通道內的流動具有尺度小、多尺度、相界面與邊界復雜的特點，傳統的計算流體力學(CFD)方法作為宏觀模擬方法存在著諸多不足，而格子玻爾茲曼方法(lattice Boltzmann method，LBM)突破了傳統計算方法的框架，直接從離散模型出發，通過粒子群的碰撞和遷移代替傳統的連續流體模型，更接近流動的微觀本質，在微流控領域具有明顯的優勢[2-3]。

格子玻爾茲曼方法的核心思想是將流體離散為在網格上運動的介觀粒子，通過計算粒子的碰撞和遷移規律得到粒子分布函數，進而統計計算得到宏觀變量如壓力、速度等分布規律，創造性地實現了模擬流體運動的連續介質模型向離散模型的轉變[4]。由于LBM方法基于非平衡統計物理學的Boltzmann方程，因而能成為聯系微觀分子尺度與宏觀尺度之間的紐帶[5-6]。傳統的CFD方法主要基于宏觀的連續介質假設，而難以計算那些不符合連續介質假設或者難以用宏觀方程描述的系統，對于這些體系往往需要借助微觀的分子動力學或氣體動理論來進行描述[7]。對于分子動力學來說必須同時跟蹤大量粒子的運動，實際求解的計算量非常大。在這種背景下，基于分子運動論和概率統計力學的LBM方法就成為一種有效方法，其具有更高的計算效率，并且容易實現并行計算[8-10]。

繞流是工程中以及自然界中常見的流體運動行為，化工反應器中的攪拌槳葉的轉動、液相流過催化劑層、擋板及其他內構件等現象都屬于繞流，這種廣泛發生的物理現象包含著深刻的傳遞機理[11-12]。研究者已針對繞流的數值模擬開展了很多卓有成效的工作，有限差分、有限體積以及有限元等離散方法，大渦模擬、雷諾平均納維-斯托克斯、離散渦等流場計算方法，被廣泛應用于繞流行為的模擬研究[13-14]。

由于微通道特征尺寸較小(一般在10～1 000 μm之間)，其內部流動機理與宏觀尺度存在顯著的差異，宏觀尺度下的許多規律不再適用于微通道中，目前對于微通道中的繞流行為理解仍不夠深入，實現微反應器的流動精確控制及科學設計仍舊存在困難。

本文闡述了格子玻爾茲曼方法的原理和計算方法，并應用該方法構建數值模型進行了宏觀尺度及微通道中的繞流計算模擬，結果可為相關多相流動過程提供理論指導與數據支撐。

1 LBM基本原理與方法

LBM的前身為“格子氣”(Lattice Gas Automata)或“元胞自動機”模型，Boltzmann方程是統計力學中用以描述非平衡態分布函數演化規律的方程，在Boltzmann方程的求解過程中，其難點在分布函數的非線性項——碰撞項，該非線性項還與分子間的具體作用力有關。格子氣自動機方法有存在隨機噪音、碰撞算子的指數復雜等缺點，20世紀80年代研究者通過引入平衡態分布函數將碰撞算子線性化，降低了碰撞算子的指數復雜性，直到90年代提出的單松弛時間法和LBGK模型，進一步簡化了碰撞算子，此后LBM逐漸發展成熟，近年來已成為流體力學領域研究熱點之一[15]。

BGK近似通過一個簡單的算子來替代碰撞項從而簡化Boltzmann方程，而格子玻爾茲曼方程是Boltzmann—BGK方程的進一步離散形式，這一離散形式包括了速度離散、時間離散、空間離散。至今LBM模型中，使用最為廣泛的單松弛模型或LBGK模型方程如下：

式中：τ——無量綱松弛時間；

t——時間；

fi——分布函數；

feq——局部平衡分布函數。

如果將數值粘度吸收到物理粘度中，則LBGK方程的空間和時間精度都是二階的。對上式根據特定時間步長Δt積分可得：

LBM的粒子分布的網格模型通常用DdQm模型表達，其中前者代表網格維度，后者代表粒子可能的運動方向，平衡態分布函數可統一表示為：

式中：wi——權重因子；

c——粒子速度，m/s；

cs——格子聲速，m/s。

目前，最常用的基本模型有D2Q7、D2Q9、D3Q15、D3Q19、D3Q27模型等，網格結構如圖1所示。

圖1 LBM方法中的主要網格模型

本文中使用的D2Q9模型的速度配置如下：

對于格子玻爾茲曼方程的初始化有兩種方法——非平衡態校正方法和迭代方法，前者基本思想是基于Chapman-Enskog展開，求出分布函數的高階項近似，進而得到分布函數的近似表達式；后者基本思路是求解初始壓力的Poisson方程，同時初始化分布函數，以求得到與速度場相一致的初始分布函數。此外，邊界處理是格子玻爾茲曼方法不同于傳統CFD數值模擬的重要內容，根據計算格式的不同主要可分為啟發式、動力學、外推以及其他復雜邊界處理格式。

綜上所述，完整的格子玻爾茲曼模型主要分為三個部分：格子，即離散速度模型；平衡態分布函數；分布函數的演化方程。其計算過程可以歸納為：根據研究對象特性，首先確定所有節點上的宏觀物理量，進而計算各節點在各方向上的平衡態分布函數，獲得初始場；求解離散后的控制方程，針對不同的邊界節點特征采用不同的邊界處理格式，基于格子玻爾茲曼模型的宏觀量的定義法則，計算各節點的物理量分布。

2 繞流過程的LBM模擬

2.1 宏觀尺度繞流模擬

本文首先對典型的宏觀尺度流體繞圓柱流動行為進行模擬，采用D2Q9模型與反彈格式邊界條 件(即粒子碰撞固體邊界后反彈回流動域，保證邊界的質量與動量守恒)。圓柱體半徑為0.25 m，高度為1 m，將其置于一個x方向風速為10 m/s的風洞中。對于模擬重點關注圓柱體周圍及尾流區域的格子分布進行了加密，如圖2所示。

圖2 圓柱體網格設置及其周圍網格密度變化示意

經過模擬得到圓柱周圍繞流流體速度分布(俯視圖)隨時間的變化如圖3所示。

圖3 宏觀尺度繞流流場隨時間變化規律

由模擬結果可知，繞流經過約2 s之后達到穩定狀態，障礙物兩側周期性地脫落出旋轉方向相反、并排列成有規則的雙列線渦，這兩列線渦分別保持自身的運動前進，進而互相干擾、吸引，形成了非線性的渦街。渦街的形成原理是由于繞流在任何時候左右都不會完全對稱，某一時刻必有一方壓力高，而另一方壓力低，所以必定產生相對運動，與此同時，流動使得低壓區也在發展，當其壯大到一定的程度后開始反向運動，至壓力平衡時移動停止，導致尾渦脫落。由于流體在不斷的流動，上述的過程往復進行，最終實現周期性的尾渦脫落現象。

流體力學中渦量是描寫旋渦運動強度和方向的最重要的物理量之一，其數學定義為速度場的旋度，渦量的大小是流體微團繞該點旋轉的平均角速度的兩倍，方向與微團的瞬時轉動軸線重合。本研究中繞流渦量隨時間變化規律如圖4所示，可以看到尾渦變化經歷了對稱——混沌——周期性穩定的過程，在本條件下，尾渦脫落周期約為0.1 s。

圖4 繞流渦量隨時間變化規律

本案例為圓柱繞流的經典體系，已有物理學家通過理論推導得出x方向上阻力系數理論值約為1，而y方向上以原點為中心呈周期性振蕩[16]。通過本文模擬結果計算障礙物的受力情況，可以得到繞流經過圓柱體時x、y兩個方向的阻力系數，分別如圖5顯示，可以發現模擬得到的x、y兩個方向的阻力系數均與理論值非常接近，驗證了本模擬模型的準確性。

圖5 繞流阻力系數隨時間變化

2.2 微通道中的繞流過程數值模擬

本節中使用2.1節建立的LBM模型，研究微通道尺度下的繞流行為，與2.1節算例保持繞流停留時間一致，圓柱體半徑為0.25 mm，高度為10 mm，將其一個水平風速為0.01 m/s的風洞中。模擬得到的速度云圖如圖6所示。

圖6 微通道繞流流場穩定后速度分布

本研究條件下微通道繞流流場與宏觀尺度相比，能夠迅速達到穩定狀態，而且不發生周期性變化。由圖可知，微通道繞流分布整體更加均勻對稱，障礙物對微通道繞流的阻礙影響范圍小于5倍直徑范圍，而且尾渦區域更接近層流的梯度分布，湍動效果不明顯。計算流場的渦量分布，如圖7所示。

圖7 微通道繞流流場穩定后渦量分布

由渦量云圖變化圖可以發現，并未出現渦街脫落現象。在2.1節的案例中，由于液相流速比較高，物體兩邊的渦開始互相影響，也就有了所謂的卡門渦街。而本條件下液相流速比較小，物體兩邊產生的渦效應較弱，互相影響較小，因此難以形成周期性的渦街，最終體現為障礙物與繞流之間的作用力更加穩定。

選取圓柱中心所在水平線(x方向)計算尾流流動信息，如圖8所示，其中原點位置為圓柱右側邊緣。可以發現障礙物后的速度分布呈遞增規律，且增幅逐漸放緩，在障礙物后3 mm附近(約5倍直徑距離)達到最大值，此后基本保持穩定。而通過渦量分布可以發現，在尾流中心線上渦量先增后減，在0.5 mm附近(約1倍直徑距離)達到最大值，此后逐漸下降，同樣是在障礙物后3 mm附近達到穩定值。

圖8 障礙物中心線上尾流流動變化

3 結論

本文介紹的LBM方法基本的計算變量是微觀上的粒子分布函數，通過求解離散玻爾茲曼方程，采用碰撞模型(如BGK模型)來模擬牛頓型流體流動，從而代替求解N-S方程，模擬得到的有限個粒子的流動和碰撞等相互作用成為宏觀粘性流體流動的縮影。對繞流過程中流場結構、固體阻力、脫落周期等變化規律的分析表明，格子玻爾茲曼方法計算穩定可靠，效率高，能夠用于科學地計算微反應器及微流控領域的數值模擬；同等停留時間條件下，微反應器中的圓柱繞流相對宏觀尺度湍動程度明顯降低，尾流未形成周期性的渦街，在實際生產過程中體現為流動更加均勻、可控，有助于實現化學反應的精確控制。