不同格子Boltzmann模型在毛細管流動模擬中的應用

馬瑞程 辛顯康 李宜強 管錯 郭虎

1. 中國石油大學(北京)提高采收率研究院；2. 中國石油大學(北京)石油工程學院；3. 中國石油大學(北京)繼續教育學院

流體在多孔介質中的流動模擬是油氣田開發的重要研究方向。目前，多數流動模擬方法依據連續性方程與達西定律中平均流速和壓力梯度的線性關系，在宏觀尺度上對求解空間和時間進行離散求解。其中，連續性方程的本質仍然是Navier-Stokes方程的特殊形式之一。實際多孔介質的孔隙空間微結構形狀復雜、尺寸不一，流體在多孔介質中的實際流動情況也更加復雜。近年來，從微觀層面更精細地研究滲流規律成為一個熱門話題。格子Boltzmann方法(LBM)邊界處理簡單，在邊界幾何形狀復雜的多孔介質中具有明顯的優勢［1-5］。利用LBM在孔隙尺度層面上研究流體流動過程，如測定多孔介質絕對滲透率等，能夠進一步較精確地解決真實復雜的多孔介質流動問題。

1 不同格子Boltzmann模型以及邊界條件

1.1 標準不可壓縮格子Boltzmann模型

對不可壓縮流體的流動模擬，單松弛模型(LBGK)是目前應用最廣泛的格子模型。其中，文獻［6］提出的DnQb模型(n代表維數，b代表速度離散方向的個數)使用最普遍，常用的模型包括D2Q9、D3Q15和D3Q19等。在這些模型中，平衡態分布函數可以表示為

對于D2Q9模型有

式中，fieq為i方向平衡態分布函數，ωi為權系數，ρ為流體密度，ci為離散速度的方向矢量，cs為格子聲速，u為流體流速矢量，c為格子速度。

格子速度數值上等于離散后的空間步長和時間步長的比值。在確定了模型參數和輸入參數后，即可進行流動模擬。重復以下3個步驟，即可進行不可壓縮Navier-Stokes方程的模擬。

(1)碰撞步驟。

(2)遷移步驟。

(3)計算宏觀量。

式中，fin為第n時間步i方向的分布函數；fieq為i方向平衡態分布函數；τ為無量綱松弛時間，其大小和流體的運動黏度有關，通常取值范圍為0.5～1.5［7］；t為時間變量；δt為時間步長；x為節點坐標矢量；p為流體壓力；υ為流體運動黏度。

1.2 改進的不可壓縮格子Boltzmann模型

盡管LBGK模型應用廣泛，然而其在模擬不可壓縮流體流動時，需要在流動馬赫數以及壓力變化等參數合適時才能夠保證足夠的精度，模型的實用性大打折扣。為了克服這一問題，人們基于LBGK模型提出了一系列的改進模型，如Zou-Hou模型、He-Luo模型以及D2G9模型等［7］。重點研究二維模擬和三維模擬均適用的Zou-Hou模型以及He-Luo模型。

1.2.1 Zou-Hou模型

Zou-Hou模型的基本理論［8］：如果將D2Q9模型的平衡態分布函數表達式更改為

同時，壓力和運動黏度的計算公式保持不變，而宏觀密度和宏觀速度的計算表達式定義為

通過Chapman-Enskog展開可以得到宏觀微分方程為

可以發現，定常情況下，上述方程組與定常不可壓縮N-S方程組一致。

1.2.2 He-Luo模型

該模型的基本思想是利用數學方法消除平衡態分布函數中由密度變化引起的高階馬赫數相關項［9-10］。在此模型中，新的平衡態分布函數為

引入壓力函數pi

構建出以壓力為平衡態的分布函數為

式中，ρ0為不可壓縮流體的真實密度，pi為i方向上的壓力分量，pieq為i方向平衡態壓力，p0為原始狀態下的基準壓力。

壓力表達式為

因此得到基于壓力的演化方程為

基于壓力分布函數模型的宏觀變量計算公式為

1.3 LBM邊界條件

在使用LBM進行流動模擬時，除了設定初始狀況的參數，還需要設定邊界條件。將多孔介質壁面設定為反彈格式邊界，同時將出入口端面設定為恒速邊界或定壓邊界。

(1)反彈格式邊界。標準反彈格式是處理靜止無滑移壁面的一類常用格式(常用于固液邊界)。如圖1所示，該格式假設粒子遷移到固體邊界后，將會在下一時間步沿遷移的反方向彈回［11］。碰撞后，壁面上的分布函數為

式中，

j

為邊界反彈方向，

為邊界反彈后的分布函數，

為邊界反彈前的分布函數，

x

b

為壁面格點坐標，

x

f

為流體所在格點坐標

為壁面反彈后的速度矢量，

cj

為壁面反彈前的速度矢量。

圖1 反彈邊界格式Fig. 1 Sketch of rebound boundary format

(2)非平衡外推格式。除了用于壁面的反彈格式，流動入口和出口的邊界條件仍需要設定。針對管流流動和滲流流動的入口和出口，定流速邊界和定壓力邊界是常用邊界。采用非平衡外推格式得到邊界表達式［12］。

恒速邊界為式中，ub為流體在邊界的速度，ρ(xf,t)為t時鄰近邊界流體域格點密度，ρb為邊界節點的密度。

2 LBM模型程序實現及模擬結果

2.1 程序編寫及單位轉換

研究中LBM模擬器采用MATLAB進行編寫，其主要流程如圖2所示。

圖2 LBM程序執行流程Fig. 2 Execution flow of LBM program

在程序進行初始化時，要將輸入的國際單位制物理參數轉化為格子單位。lu代表格子長度，ts代表時間步長，mu代表格子質量，依據量綱導出格子系統中其他物理量的基本單位。

不可壓縮流體在毛細管軸向剖面上的流動是一個二維問題，因此采用D2Q9模型進行求解。如圖3所示，假設毛細管長度為L，半徑為R，流動流體黏度為μ，在壓差Δp下為黏滯性滲流，滲流入口速度為Uin。流體潤濕毛細管壁的情況下，管壁出的液體流速為0，在管中軸線處流速最大。

圖3 流動模擬情境Fig. 3 Flow simulation scenario

考慮到黏滯力和驅替壓力的平衡關系，得出Poseuille流動的中的流速表達式為

在沿徑向的剖面上，流動平均速度可以表達為

式中，r為毛細管中某一點距毛細管中軸線的徑向距離，m；為距離毛細管中軸線的徑向距離為r處的流速，m/s；pin為入口壓力，Pa；pout為出口壓力，Pa；R為毛細管半徑，m；μ為流體黏度，Pa· s ；L為毛細管長度，m;為毛細管中流體平均流速，m/s。

通過上式可以看出，入口速度剖面應當為一條拋物線，而平均流速和驅替壓差呈線性關系。

將輸入參數進行轉化后即可進行Poseuille流動模擬。其中，模擬需要的松弛時間和驅替壓差折算等效密度為

式中，μD為流體黏度，ρD為無量綱流體密度，pDin為入口壓力，ρDin為入口等效密度，pDout為出口壓力，ρDout為出口等效密度。

2.2 程序模擬結果及誤差分析

編寫模擬器并運行后，可得到二維Poseuille流動時速度分布圖像，從圖4可以看出，在求解定常問題時，隨著時間步數的增加，速度場的分布逐漸均勻，最終趨于收斂。此外，徑向上的速度剖面為一拋物線，且毛細管壁處的速度為0，這與Poseuille解析解的基本假設一致。為了確認模擬結果是否正確，設毛管半徑為3 μm，毛細管長度為7.5 μm，流體密度為1 000 kg/m3，流體黏度為1 mPa · s，兩端的驅替壓差為700 Pa。經過計算，LBM的平均模擬結果為0.034 4 m/s，而解析解計算結果為0.035 m/s，二者相對誤差為1.71%。求解收斂時的速度場分布圖和毛細管中各個位置的速度剖面情況如圖5和圖6所示。計算結果表明，LBM在模擬液體在毛細管中的低雷諾數流動上具有足夠的準確性，加之LBM在編程上較為簡潔且易于并行計算，因此，其在處理微觀層面上的多孔介質流動問題具有明顯優勢。

圖4 不同時間步數無量綱速度場分布Fig. 4 Distribution of dimensionless velocity field of different time steps

圖5 真實算例無量綱速度場分布Fig. 5 Distribution of dimensionless velocity field of a real case

圖6 真實不同位置出口速度剖面Fig. 6 Exit velocity section of different real positions

為了進一步對3種LBM模型的性能進行評價，選取相同的Poseuille流動算例對3種不同模型的計算結果與解析解結果進行比對。在計算過程中，為便于評價和分析，算例中的輸入參數全部為已經預先處理好的無量綱參數，模擬輸入參數及模擬結果如圖7所示，可以看出，隨著驅替壓差折算后的當量密度不斷增加，3種模型計算出的無量綱平均速度盡管也隨之增加，但都與解析解結果產生了誤差。其中，Zou-Hou模型與解析解契合最好，He-Luo模型次之，LBGK模型最差。

為了對其精確性和穩定性進行分析，作出不同模型在不同馬赫數下計算結果的相對誤差曲線，結果如圖8所示，可以看出，3種模型的相對誤差均隨著馬赫數的增加而增加，誤差增幅從小到大排序依次為：Zou-Hou模型＜He-Luo模型＜LBGK模型。

圖7 不同模型及解析解計算結果對比Fig. 7 Comparison between the calculation results of different models and the analytical solutions

圖8 不同模型相對誤差隨馬赫數變化曲線Fig. 8 Variation of the relative error of different models with the Mach number

在標準LBM模型中，動量方程與不可壓縮NS模型不完全一致［13］。如果流動馬赫數極小時，則LBGK模型對應的宏觀方程與標準的不可壓縮Navier-Stokes方程近似。在流體力學中，常用馬赫數作為衡量流體壓縮程度的無量綱準數，其定義為流體流速和聲速的比值。流體力學原理表明，當馬赫數大于0.3時，流體的壓縮性將不可忽略。而LBGK方法是求解不可壓縮N-S方程的一種人工壓縮解法，其相對誤差會隨著馬赫數的增加而逐漸增加，因此LBGK模型適用的馬赫數范圍較小。Zou-Hou模型和He-Luo模型的誤差盡管也會隨著馬赫數的增加而增加，但二者均可以通過數學方法消除與馬赫數有關的誤差項，相對誤差增幅相對較低。

在LBM中，時間步長往往和格子長度滿足一定的相關關系［14］。在多孔介質流動模擬中，如果直接對物理問題進行折算，時間步長通常在10-12s到10-10s數量級之間，模擬效率很低。如圖9所示，為了提高運行效率，需要保證雷諾數不變，提高馬赫數對物理問題進行放縮從而增大時間步長［15］。此時，可以放縮的最大倍數受到模型精確性的制約。3種模型的運算效率從高到低依次為：Zou-Hou模型＞He-Luo模型＞LBGK模型。因此，在編寫二維和三維通用模擬器時，如果求解問題為定常問題，推薦使用Zou-Hou模型進行求解；而如果求解問題為非定常問題，則推薦使用He-Luo模型進行求解。

圖9 系統放縮前后對比Fig. 9 Comparison before and after the system zooming

3 結論

（1）格子Boltzmann方法可以有效解決低雷諾數、低流速條件下，流體在毛細管中的流動模擬問題且具有足夠精度。LBM編程簡單，結果可靠，在模擬微觀層面的流體流動上具有明顯優勢。

（2）在相同馬赫數條件下，模型的相對誤差從大到小排序為：LBGK模型＞He-Luo模型＞Zou-Hou模型。在相對誤差相同條件下對同一問題進行求解，程序循環步數由大到小依次為：LBGK模型＞He-Luo模型＞Zou-Hou模型。

（3）在二維及三維問題的模擬中，如果求解問題為定常問題，則Zou-Hou模型效果更佳；如果求解非定常問題，則推薦使用He-Luo模型。