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

程海龍

（河南省機(jī)場集團(tuán)有限公司，河南 鄭州 450000）

0 引言

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

1 物理模型及相關(guān)假設(shè)

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

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

2 控制方程

2.1 D2Q9模型

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

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

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

為恢復(fù)宏觀方程，平衡態(tài)分布函數(shù)需滿足的矩方程見式（2）至式（4）。

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

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

3 計算流程

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

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

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

③執(zhí)行遷移，見式（11）。

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

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

4 數(shù)值方法

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

5 結(jié)果分析

5.1 速度驅(qū)動

按照上述設(shè)定進(jìn)行仿真模擬試驗，試驗結(jié)果見圖2、圖3。

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

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

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

5.2 壓力驅(qū)動

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

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

圖5 速度矢量圖

6 結(jié)語

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