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

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

0 前言

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

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

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

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

本文闡述了格子玻爾茲曼方法的原理和計算方法，并應(yīng)用該方法構(gòu)建數(shù)值模型進(jìn)行了宏觀尺度及微通道中的繞流計算模擬，結(jié)果可為相關(guān)多相流動過程提供理論指導(dǎo)與數(shù)據(jù)支撐。

1 LBM基本原理與方法

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

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

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

t——時間；

fi——分布函數(shù)；

feq——局部平衡分布函數(shù)。

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

LBM的粒子分布的網(wǎng)格模型通常用DdQm模型表達(dá)，其中前者代表網(wǎng)格維度，后者代表粒子可能的運動方向，平衡態(tài)分布函數(shù)可統(tǒng)一表示為：

式中：wi——權(quán)重因子；

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

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

目前，最常用的基本模型有D2Q7、D2Q9、D3Q15、D3Q19、D3Q27模型等，網(wǎng)格結(jié)構(gòu)如圖1所示。

圖1 LBM方法中的主要網(wǎng)格模型

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

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

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

2 繞流過程的LBM模擬

2.1 宏觀尺度繞流模擬

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

圖2 圓柱體網(wǎng)格設(shè)置及其周圍網(wǎng)格密度變化示意

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

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

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

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

圖4 繞流渦量隨時間變化規(guī)律

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

圖5 繞流阻力系數(shù)隨時間變化

2.2 微通道中的繞流過程數(shù)值模擬

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

圖6 微通道繞流流場穩(wěn)定后速度分布

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

圖7 微通道繞流流場穩(wěn)定后渦量分布

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

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

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

3 結(jié)論

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