基于格子Boltzmann法的管道流模擬二維和三維離散模型的比較

葛欽欽，馬新華

(江蘇大學 流體機械工程技術研究中心，江蘇 鎮江 21203)

基于格子Boltzmann法的管道流模擬二維和三維離散模型的比較

葛欽欽，馬新華

(江蘇大學 流體機械工程技術研究中心，江蘇 鎮江 21203)

格子Boltzmann方法是近二十年迅速發展起來的一種介觀模擬方法。它是由格子氣自動機理論演化而來，有別于傳統模擬方法，具有宏觀上離散，微觀上連續的特點。本文著重介紹了格子Boltzmann法控制方程的理論基礎，基于LBGK理論的D2Q9和D3Q19離散速度模型邊界條件的定義以及算法演變情況。使用MATLAB軟件編寫程序進行二維和三維的泊肅葉流的數值模擬，對比分析D2Q9和D3Q19離散速度模型在LBGK理論下的邊界條件的影響。

格子Boltzmann法；LBGK理論模型，D2Q9離散模型；D3Q19離散模型；MATLAB軟件

此篇論文的研究工作是基于格子Boltzmann方法(LBM)開展的。這是一種近二十年快速發展起來的較為年輕的數值方法。它是基于粒子動理論，由格子氣自動機(LGA)理論演化而來，在系統一定量的離散格子進行計算，通過對大量格子的統計平均值獲得宏觀結果。LBM屬于介觀模擬方法，即可以認為該方法在宏觀上是離散，在微觀上是連續的。與傳統的流體力學計算方法相比，LBM具有顯著的優勢，因為它具有天生的并行特性，邊界條件設置簡單，程序易于實現等特點，尤其適用于微尺度流型，多相多組分流型和多孔介質流型等許多傳統模擬方法難以勝任的領域。

1 數值算法

1.1 Boltzmann傳遞方程

Boltzmann傳遞方程由Boltzmann，Mcnamara and Zanetti三位學者于1988年首先提出。其重要意義是引入統計概念，在宏觀現象和微觀基礎之間建立了橋梁[1]。Boltzmann方程的基本思想是不去研究每個粒子的具體運動狀態，而是研究大量粒子處在某一運動狀態下的概率，通過統計平均的方法得出系統的宏觀參數。波爾茲曼方程的推導是基于以下三個假設條件下的[2]。

(1)粒子相互碰撞時只考慮兩體碰撞，即三個或三個以上的粒子同時發生碰撞為小概率事件可忽略不計。

(2)粒子混沌假設。即各個粒子的速度分布是獨立的，不依賴于其他粒子。

(3)局部碰撞的動力學行為不受外力的影響。

為了刻畫粒子的分布情況，我們引入粒子分布函數f(s,u,t)，其中s(x,y,z)表示空間矢量，u(x,y,z)表示速度適量，t表示某一時刻。粒子的運動加速度表示為a。對于任意粒子在時間間隔dt內，運動變化是主體運動和分子間的相互碰撞綜合的結果。數學表達式為[3]：

(1)

泰勒展開得：

(2)

根據氣體動理論知識對于碰撞項做出數學表達[4]：

(3)

1.2 BGK模型

由方程(3)可知，由于碰撞項的存在，Boltzmann傳遞方程是一個復雜的微積分方程，給方程的求解帶來了很大的困難。所以需要對粒子間的碰撞找出近似合理的簡化才能求解方程。

(4)

式中feq(s,u)表示平衡態時的分布函數；比例系數ν和松弛時間τ成反比：

ν=1/τ

(5)

最后我們得到Boltzmann-BGK方程：

(6)

1.3 質量統計與動量統計

單個粒子的分布函數可以簡單地認為是代表各個粒子在某一狀態下(包括位置矢量和速度矢量)發生的頻率。系統的宏觀密度和速度可以由系統內所有的粒子分布函數求得[9]。

模型的宏觀密度表示如下：

(8)

系統的宏觀速度是微觀粒子速度矢量和方向密度乘積的平均值，表示如下：

(9)

系統動量的表達式：

(10)

1.4 LBM二維和三維離散速度模型

90年代初期，Qian等人在LGA理論的FHP離散模型上進行改進，提出DxQy形式(x維空間，y個離散速度)既適合LBM平衡態分布函數又具有完全對稱性的系列離散速度模型[10-12]。

其中二維模型中應用最為廣泛的是D2Q9模型，圖示如下。

圖1 D2Q9離散速度模型

D2Q9速度模型下平衡態分布函數各參數：

三維模型中應用最為廣泛的是D3Q15模型和D3Q19模型，由于D3Q19模型具有更加豐富的離散方向數，所得到模擬結果一般更加精準和可靠。同時更多的離散方向也會給計算提出更高的要求。一般在計算資源充足的條件下，建議采用D3Q19模型。模型給出如下：

圖2 D3Q19離散速度模型

1.5 邊界條件

管道流需要處理的邊界條件有進出口處和對管壁的邊界條件處理。以D2Q9速度模型在管道入口即左邊界為例分析LBM邊界條件，圖示如下：

圖3 管道進口處的D2Q9邊界條件

由方程(8)可知

(11)

由方程(11)，對速度矢量u正交分解成ux和uy得出：

(12)

(13)

此時對于左邊界情況，未知條件有f1，fs，f8和邊界條件處的 ，所以我們需要四個獨立方程來求解這四個未知數。現在已有方程(11)，(12)和(13)。Zou -He 邊界條件法給出了簡化，假設粒子處于非平衡態時鋼球反彈規則仍然適用，那么我們則得到了第四個需要的方程[13]：

(14)

最終根據以上四個方程求解得出所需的未知邊界參數：

對于三維模型D3Q19，我們依然以左邊界為例加以說明，未知量有以下六個：p,f1,f7,f9,fu,fb。同二維推到類似，根據方程

(16)

(17)

對速度矢量u正交分解得

由方程組

ux是進口速度，為已知量。故可求：

在左邊界上有Zou -He假設：

(21)

同時滿足BGK模型的平衡態分布函數

可由以上(21)(22)方程求得

2 兩種速度模型的泊肅葉流型的分析對比

2.1 基于Navier-Stokes方程的解析解

設定管道流體的密度為常數，即不可壓縮流體。流體流動的驅動力恒定為F。

(25)

(26)

式中，L為管道的半徑，y是計算點處距離管道中心的距離，υ是流體的運動粘度。

2.2 D2Q9數值模擬結果

泊肅葉流是粘性流體在管道中的典型流型，初始條件如圖(4)所示：

圖4 二維管道泊肅葉流的結構示意圖

數值模擬在Matlab軟件下進行。流域網格為150*50，雷諾數為100，松弛時間根據松弛系數公式omega = 1/ (3*nu+1/2))得出，其值算出得0.65。進出口邊界條件采用Zou-He邊界條件，管壁采用剛體反彈邊界條件處理。得出模擬結果在平衡態下結果如圖(5)所示：

圖5 二維管道泊肅葉流平衡態呈現

與解析解結果對比如圖(6)所示：

圖6 泊肅葉流解析解與D2Q9模擬結果的對比

可以看出LBM理論的數值模擬結果基本吻合解析解的結果。我們通過計算相對誤差發現系統的最大誤差出現在靠近管壁處。

(27)

管道各位置的相對誤差如圖(7)所示

圖7 管道各位置解析解與模擬結果的相對誤差

這說明對于管壁的邊界條件處理我們用的最為理想的剛體彈性方法并不能很精確的反應真實的邊界情況。但是優點是可以更加清晰簡便的描述問題，減小計算負擔。

2.3 D2Q9模型下的質量守恒和能量守恒驗算

為了檢驗LBM算法的穩定性，在管道的進出口處進行質量守恒和能量守恒的驗算。進出口的質量和能量計算如方程(28)和(29)

取在不同的進口速度，雷諾數下的計算結果匯總如表格。

150*50;u=0.1;Re=50MassEnergyinlet4.40990.0287outlet4.40990.0287

數值模擬結果表明LBM理論采用D2Q9速度模型，有著很好的穩定性和精確度。

2.4 D3Q19數值模擬結果

結構示意如圖所示：

圖8 三維管道泊肅葉流的結構示意圖

三維管道泊肅葉流采用D3Q19模型進行計算，結。初始條件，在管道進口處給一個水平速度( ， )，雷諾數為100，松弛時間設定為1。進出口邊界條件采用引入修正變量的Zou-He邊界條件，管壁采用剛體反彈邊界條件處理。得出模擬結果在平衡態下結果分別在Y-Z, X-Z 和 X-Y界面進行切片，圖示如下：

圖9 管道進口，中心和出口處在Y-Z切面上的速度分布

圖10 沿Y軸中心的X-Z切片上的速度分布

圖11 沿Z軸中心的X-Y切片上的速度分布

2.5 D3Q19模型下的穩定性驗算

管道橫截面為正方形對稱結構。根據對稱性原理，在算法穩定的情況下，在對應點上X-Y和X-Z切面的上對應點的速度值也應該相等[16]。我們截取管道在出口中心處的X-Y和X-Z切面上的速度分布進行匹配。同時把該案例用D2Q9速度模型進行模擬，與三維結果進行比照。圖示如下：

通過對比圖像，我們發現在D3Q9速度模型下，出口處的X-Y和X-Z切面的上的速度能夠很好的匹配，表明LBGK模型用D3Q19速度模型下的算法是穩定可靠的。同時與在二維的模擬中出口處所得到的速度分布進行對比，我們發現所得到的速度分布趨勢大致相同，特別是在流域中心處更為接近。但是D3Q19模型出口處各點的速度分布均略大于使用D2Q9模型下的速度。根據流體力學知識，靠近管壁處的層流內層速度梯度無限小接近于0，二維的模擬結果很好的符合這一特性，而三維的模擬結果靠近管壁處速度雖然很小但是沒有達到0。分析三維數值模擬出現偏差的原因可能是管壁邊界條件假設的過于理想化。無滑移剛體反彈的邊界條件是使用最為廣泛的，原因是把邊界條件做出理想化假設達到簡化的效果，從而易于實現和理解。在二維的模擬中可以做出定性的分析能保證一定的正確性。但是在三維較為復雜的情況下就會與真實情況出現偏差。

圖12 泊肅葉流在二維和三維速度模型下的結果對比

[1] A A. Lattice Boltzmann method fundamentals and fngineering applications with computer codes[M]. London: Springer London Dordrecht Heidelberg New York Springer-Verlag London Limited,2011.

[2] 何雅玲,李 慶,王 勇,等. 格子Boltzmann方法的工程熱物理應用[J]. 科學通報,2009,18:2638-2656.

[3] Sauro Succi. The lattice Boltzmann equation, for fluid dynamics and beyond[M]. Oxford ：Oxford Science Publications,2001.

[4] Chen S, Doolen G D. Lattice Boltzmann method for two phase flow modeling[J]. Annual Review of Fluid Mechanics, 1998, 30: 329-364.

[5] Wolf-Gladrow Dieter A. Lattice-gas cellular automata and lattice Boltzmann models:an introduction[M]. [S.I.]:Spring,2000.

[6] 王 亮. 基于格子Boltzmann方法的非常規顆粒兩相流的機理研究[D].武漢：華中科技大學,2012.

[7] 胡安杰. 多相流動格子Boltzmann方法研究[D].重慶大學,2015.

[8] He X，Chen S, Doolen G D. A novel thermal model for the Lattice Boltzmann method in incompressible limit[J]. Journal of Computational Physics, 1998, 146(1):282-300.

[9] Peng Y, Shu C, Chew Y T. Simplified thermal lattice Boltzmann model for incompressible thermal flows[J]. Physical Review E,2003, 68(2): 026701-026709.

[10] Aidun C K. Lattice-Boltzmann method for complex flows[J]. Annual Review of Fluid Mechanics, 2009, 42(1):439-472.

[11] He X, Luo L S. Lattice Boltzmann model for the incompressible navier-stokes equation[J].Journal of Statistical Physics, 1997,88:927-944.

[12] 鄧平平. 基于格子Boltzmann方法的二維多孔介質滲流研究[D].大連理工大學,2014.

[13] 陳 柔. 格子Boltzmann方法的若干應用[D].杭州：浙江師范大學,2013.

[14] ames Maxwell Buick.Lattice Boltzmann methods in interfacial wave modelling[D]. Edinburgh: University of Edinburgh,1997.

[15] 李 元. 格子Boltzmann方法的應用研究[D].合肥：中國科學技術大學,2009.

[16] 王福軍.計算流體動力學分析--CFD軟件的理論與應用[M].北京：清華大學出版社，2004.

(本文文獻格式：葛欽欽，馬新華.基于格子Boltzmann法的管道流模擬二維和三維離散模型的比較[J].山東化工，2016，45(12)：146-150，153.)

2016-04-15

TQ019

A

1008-021X(2016)12-0146-05