基于格子Boltzmann方法的泊肅葉流數值模擬研究

程海龍

（河南省機場集團有限公司，河南 鄭州 450000）

0 引言

格子Boltzmann方法（Lattice Boltzmann Method，LBM）被學者提出來已有近30 年的時間［1］。經過幾十年的發展，LBM 理論及工程應用研究等方面都有新的進展，并取得了卓越的成果。近年來，LBM 理論及應用研究成為研究熱點之一，受到國內外眾多學者的重視。其基本思想是將系統中流體離散成大量的流體粒子，在簡單劃定的網格上，有規律地進行反復碰撞和遷移，給出一個簡化的動力模型。然后通過統計來求取平均值，得到宏觀物理量（如密度、速度）的流體動力學方程［2］。在實際應用中，有很多以壓力進行驅動的不可壓縮流體有周期性流動的現象［3］，特別是在工業生產及醫學研究中，泊肅葉流的應用效果顯著。最初泊肅葉流動是為了對血管中的血液流動進行描述的，在生物力學中有著重要的地位，且許多微尺度換熱器中的流體流動也是用泊肅葉定律進行描述的。

1 物理模型及相關假設

如圖1 所示，兩個平行的平板間充滿流體，以沿x方向的壓強梯度為動力，推動整個流體流動，這就叫泊肅葉流［4］。該流動是流體力學中十分經典的流動問題，且可用Navier-Stokes 直接算出其解析解［5］，因而對基于格子Boltzmann方法所開發的程序，經常使用該算例進行驗證［6］。假定所述研究對象滿足以下假設：平行板之間的流動為不可壓流動；流動為二維流動，且只考慮XY方向；滿足剛性固體壁的邊界條件。

圖1 二維平板間的流動模型圖

2 控制方程

2.1 D2Q9模型

二維定常流動一般采用二維九速的D2Q9模型進行描述，該模型中的每個節點上都有9 個不同的速度方向。其中，每個節點上允許有一個靜止粒子存在，加上與其相鄰的8 個節點，記為D2Q9 模型。模型采用的是D2Q9模型平衡態分布函數，見式（1）。

式中：f eqα為平衡態分布函數；eα為速度向量；cs為格子聲速；ωα為權系數；ρ為宏觀密度；u為宏觀速度。在D2Q9模型中，cs=c=δx/δt，δx為網格步長，δt為時間步長，假定x和y方向的網格步長均相等，且等于時間步長，單位均為1，即δx=δy=δt，可得c= 1。

2.2 宏觀方程及格子Boltzmann-BGK 方程

為恢復宏觀方程，平衡態分布函數需滿足的矩方程見式（2）至式（4）。

式中：δij=ei·ej，且p=ρc2s；u為宏觀速度；δij為一個矢量；p為宏觀壓力。因此，模型的宏觀密度、速度、壓力定義見式（5）至式（8）。

式（8）是考慮外力項影響時的格子Boltzmann-BGK 方程。式中：Fα(x,t) 是外力分布函數；τ=τ0/δt為無量綱松弛時間，它直接影響流體的黏性和宏觀屬性。對上述方程進行求解就是格子Boltzmann方法的核心。

3 計算流程

本研究以visual studio 平臺為整個程序的運行環境，編程語言選用C++。用LBM 模型對泊肅葉流進行數值模擬計算。該方法在商業軟件數值模擬中具有優勢，可有針對性地進行計算，通過簡單程序可循環進行碰撞和遷移兩個過程，直到符合設置的精度，因而程序的編寫和實施容易實現，且對計算機系統的要求不高。本研究的程序計算流程如下所示。

①初始化整個流場，確定各個節點的宏觀量，并同時初始化分布函數，見式（9）。

②在時刻t執行碰撞，見式（10）。

③執行遷移，見式（11）。

④計算宏觀量，見式（12）。

⑤重復步驟②③④，直到滿足終止條件。

4 數值方法

本研究編寫了基于LBM 模型的程序，通過運行程序，將模擬值、相應的解析值與已有的相關文獻給出的結果進行分析對比，以驗證所選模型和開發出的程序的正確性。用格子Boltzmann 方法對模擬二維速度驅動下的泊肅葉流動和用商業軟件有限容積法模擬得出的結果進行對比，二者的結果非常吻合。本研究采用D2Q9 模型，分別對速度驅動和壓力驅動進行模擬，并對模擬結果進行對比分析［7］。與現有的研究相比［8］，模擬中均采用格子單位。在速度驅動下，一是進口處給定恒定速度uin= 0.01，而壓力驅動則沒有初速度；二是給定一個沿x方向的力，大小為1.0×10-6，y軸方向的力為0，兩種情況均采用均勻網格，網格數為Nx×Ny=500×50，pr數為0.71，粒子速度分布函數對應的松弛時間、進口邊界上的宏觀參數速度已知，并假定進口處的分布函數滿足相應的平衡態分布，右端出口為充分發展邊界，粒子速度分布函數和溫度等于臨近的流體節點對應的值，上下壁面為無滑移壁面，并采用非平衡外推格式對壁面進行處理。

5 結果分析

5.1 速度驅動

按照上述設定進行仿真模擬試驗，試驗結果見圖2、圖3。

圖2 X=40截面的速度分布

圖2是X=40截面處的速度分布圖，經過比較可得，在X=40 截面處速度取得最大值，并且可以明顯看出，在流動達到穩定后，沿著平板方向的截面速度呈拋物線分布。由流體力學知識可知，在此條件下的最大速度是平均速度的1.5倍［9］。由圖2可知，模擬結果為0.014 698，與解析解最大誤差為2.01%，可以看出該方法比文獻解更接近解析解，說明格子Boltzmann 方法在處理二維平板間的不可壓縮流動問題具有一定的精確度和可行性［10］。中心線上的速度大小在x方向是逐漸減小的，這是因為壓強是線性變化的，壓強梯度是流動的動力，當達到穩定后，速度大小保持不變。在以X/Lx為橫坐標、（p-pout）/（pin-pout）為縱坐標的曲線上，壓力變化如圖3所示，入口和出口處稍有偏差，這是因為采用的邊界處理格式不同，即假定入口處分布函數完全滿足給定條件下的平衡態，而出口處采用的是充分發展格式，和實際還是存在偏差的。即在x方向上存在壓力梯度。壓力梯度也是導致平板間流體流動的原因，壓力梯度是平板間流體流動的驅動力。

圖3 截面Y/Ly=0.5的壓力曲線

5.2 壓力驅動

壓力驅動和速度驅動的最大區別在于，壓力驅動沒有給定初速度，只給了一個沿x方向的力，當流動達到穩定后，其流動情況類似于速度驅動。在y方向上幾乎沒有速度，而x方向上的速度呈對稱分布，并且從邊緣處到中間是逐漸增大的，這是因為驅動力較小時，要考慮流體的黏度和邊界層因素的影響。在邊界層的影響下，同一垂直界面上的流體的不同位置的速度不同，中間位置離邊界層最遠，受邊界層的影響最小，因此速度較大（見圖4）。圖5 給出了流動穩定后所獲得的速度場矢量圖，從圖5 可以看出，當流動達到穩定后，速度矢量均指向x方向，速度大小隨著y值的變化而變化，且呈拋物線分布，但這并不是在入口處就形成的。這是因為當流體以某一較小速度均勻進入平板間時，受黏性的影響，緊挨著壁面的流體會完全粘在壁面上，流體的速度逐漸降為0。根據流量連續原理，流過各個截面上的流量是相等的，因而越靠近平板的中心軸線，流速就越大，當流動達到穩定后，各截面上的流速分布規律就不再變化，呈拋物線形狀。

圖4 截面X/Lx=0.5的速度分布

圖5 速度矢量圖

6 結語

本研究采用格子Boltzmann 方法，選取基本的D2Q9 模型、標準的碰撞遷移，對速度驅動和壓力驅動的泊肅葉流進行數值模擬，根據泊肅葉流的特征，邊界處理采用非平衡態外推格式和充分發展格式兩種方式。試驗結果表明，基于LBM 模型得出的模擬解與泊肅葉流速度的解析解吻合緊密。同時，不同縱截面上的速度分布以及整體變化趨勢同現有的研究成果也保持一致，這說明格子Boltzmann 方法完全適合處理壓力驅動與速度驅動類的流動問題，并且具有很高的數值精度和數值穩定性，與其他方法相比，格子Boltzmann 方法具有很多優勢。通過分析速度云圖和矢量圖可知，壓力梯度是流體流動的驅動力，并且通過對比速度驅動和壓力驅動，對泊肅葉流動機理進行分析，為進一步研究流體的流動與傳熱奠定理論基礎。