葛欽欽，馬新華

(江蘇大學(xué) 流體機(jī)械工程技術(shù)研究中心，江蘇 鎮(zhèn)江 21203)

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

葛欽欽，馬新華

(江蘇大學(xué) 流體機(jī)械工程技術(shù)研究中心，江蘇 鎮(zhèn)江 21203)

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

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

此篇論文的研究工作是基于格子Boltzmann方法(LBM)開展的。這是一種近二十年快速發(fā)展起來(lái)的較為年輕的數(shù)值方法。它是基于粒子動(dòng)理論，由格子氣自動(dòng)機(jī)(LGA)理論演化而來(lái)，在系統(tǒng)一定量的離散格子進(jìn)行計(jì)算，通過(guò)對(duì)大量格子的統(tǒng)計(jì)平均值獲得宏觀結(jié)果。LBM屬于介觀模擬方法，即可以認(rèn)為該方法在宏觀上是離散，在微觀上是連續(xù)的。與傳統(tǒng)的流體力學(xué)計(jì)算方法相比，LBM具有顯著的優(yōu)勢(shì)，因?yàn)樗哂刑焐牟⑿刑匦?，邊界條件設(shè)置簡(jiǎn)單，程序易于實(shí)現(xiàn)等特點(diǎn)，尤其適用于微尺度流型，多相多組分流型和多孔介質(zhì)流型等許多傳統(tǒng)模擬方法難以勝任的領(lǐng)域。

1 數(shù)值算法

1.1 Boltzmann傳遞方程

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

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

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

(3)局部碰撞的動(dòng)力學(xué)行為不受外力的影響。

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

(1)

泰勒展開得：

(2)

根據(jù)氣體動(dòng)理論知識(shí)對(duì)于碰撞項(xiàng)做出數(shù)學(xué)表達(dá)[4]：

(3)

1.2 BGK模型

由方程(3)可知，由于碰撞項(xiàng)的存在，Boltzmann傳遞方程是一個(gè)復(fù)雜的微積分方程，給方程的求解帶來(lái)了很大的困難。所以需要對(duì)粒子間的碰撞找出近似合理的簡(jiǎn)化才能求解方程。

(4)

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

ν=1/τ

(5)

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

(6)

1.3 質(zhì)量統(tǒng)計(jì)與動(dòng)量統(tǒng)計(jì)

單個(gè)粒子的分布函數(shù)可以簡(jiǎn)單地認(rèn)為是代表各個(gè)粒子在某一狀態(tài)下(包括位置矢量和速度矢量)發(fā)生的頻率。系統(tǒng)的宏觀密度和速度可以由系統(tǒng)內(nèi)所有的粒子分布函數(shù)求得[9]。

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

(8)

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

(9)

系統(tǒng)動(dòng)量的表達(dá)式：

(10)

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

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

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

圖1 D2Q9離散速度模型

D2Q9速度模型下平衡態(tài)分布函數(shù)各參數(shù)：

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

圖2 D3Q19離散速度模型

1.5 邊界條件

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

圖3 管道進(jìn)口處的D2Q9邊界條件

由方程(8)可知

(11)

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

(12)

(13)

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

(14)

最終根據(jù)以上四個(gè)方程求解得出所需的未知邊界參數(shù)：

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

(16)

(17)

對(duì)速度矢量u正交分解得

由方程組

ux是進(jìn)口速度，為已知量。故可求：

在左邊界上有Zou -He假設(shè)：

(21)

同時(shí)滿足BGK模型的平衡態(tài)分布函數(shù)

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

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

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

設(shè)定管道流體的密度為常數(shù)，即不可壓縮流體。流體流動(dòng)的驅(qū)動(dòng)力恒定為F。

(25)

(26)

式中，L為管道的半徑，y是計(jì)算點(diǎn)處距離管道中心的距離，υ是流體的運(yùn)動(dòng)粘度。

2.2 D2Q9數(shù)值模擬結(jié)果

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

圖4 二維管道泊肅葉流的結(jié)構(gòu)示意圖

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

圖5 二維管道泊肅葉流平衡態(tài)呈現(xiàn)

與解析解結(jié)果對(duì)比如圖(6)所示：

圖6 泊肅葉流解析解與D2Q9模擬結(jié)果的對(duì)比

可以看出LBM理論的數(shù)值模擬結(jié)果基本吻合解析解的結(jié)果。我們通過(guò)計(jì)算相對(duì)誤差發(fā)現(xiàn)系統(tǒng)的最大誤差出現(xiàn)在靠近管壁處。

(27)

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

圖7 管道各位置解析解與模擬結(jié)果的相對(duì)誤差

這說(shuō)明對(duì)于管壁的邊界條件處理我們用的最為理想的剛體彈性方法并不能很精確的反應(yīng)真實(shí)的邊界情況。但是優(yōu)點(diǎn)是可以更加清晰簡(jiǎn)便的描述問(wèn)題，減小計(jì)算負(fù)擔(dān)。

2.3 D2Q9模型下的質(zhì)量守恒和能量守恒驗(yàn)算

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

取在不同的進(jìn)口速度，雷諾數(shù)下的計(jì)算結(jié)果匯總?cè)绫砀瘛?/p>

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

數(shù)值模擬結(jié)果表明LBM理論采用D2Q9速度模型，有著很好的穩(wěn)定性和精確度。

2.4 D3Q19數(shù)值模擬結(jié)果

結(jié)構(gòu)示意如圖所示：

圖8 三維管道泊肅葉流的結(jié)構(gòu)示意圖

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

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

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

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

2.5 D3Q19模型下的穩(wěn)定性驗(yàn)算

管道橫截面為正方形對(duì)稱結(jié)構(gòu)。根據(jù)對(duì)稱性原理，在算法穩(wěn)定的情況下，在對(duì)應(yīng)點(diǎn)上X-Y和X-Z切面的上對(duì)應(yīng)點(diǎn)的速度值也應(yīng)該相等[16]。我們截取管道在出口中心處的X-Y和X-Z切面上的速度分布進(jìn)行匹配。同時(shí)把該案例用D2Q9速度模型進(jìn)行模擬，與三維結(jié)果進(jìn)行比照。圖示如下：

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

圖12 泊肅葉流在二維和三維速度模型下的結(jié)果對(duì)比

[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方法的工程熱物理應(yīng)用[J]. 科學(xué)通報(bào),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方法的非常規(guī)顆粒兩相流的機(jī)理研究[D].武漢：華中科技大學(xué),2012.

[7] 胡安杰. 多相流動(dòng)格子Boltzmann方法研究[D].重慶大學(xué),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方法的二維多孔介質(zhì)滲流研究[D].大連理工大學(xué),2014.

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

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

[15] 李 元. 格子Boltzmann方法的應(yīng)用研究[D].合肥：中國(guó)科學(xué)技術(shù)大學(xué),2009.

[16] 王福軍.計(jì)算流體動(dòng)力學(xué)分析--CFD軟件的理論與應(yīng)用[M].北京：清華大學(xué)出版社，2004.

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

2016-04-15

TQ019

A

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