樂(lè)勵(lì)華, 高 云, 劉唐偉

(東華理工大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院,江西撫州 344000)

偏微分方程求解的一種新穎方法
——格子Boltzmann模型

樂(lè)勵(lì)華, 高 云, 劉唐偉

(東華理工大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院,江西撫州 344000)

介紹了一種偏微分方程求解的一種新穎方法格子Boltzmann模型,詳細(xì)分析了它的基本理論和基本原理。并通過(guò)不可壓Navier-Stokes方程組和二維含源項(xiàng)擴(kuò)散方程的數(shù)值模擬計(jì)算實(shí)例,說(shuō)明格子Boltzmann方法的有效性,展示了廣闊的應(yīng)用前景,為今后更深入的研究和廣泛應(yīng)用提供參考.

格子Boltzmann方法;平衡態(tài)分布函數(shù);D2Q9模型;Navier-Stokes方程;對(duì)流—擴(kuò)散方程

1 引 言

偏微分方程廣泛應(yīng)用在物理、化學(xué)、生物、材料、金融和廣大工程技術(shù)領(lǐng)域,大學(xué)數(shù)學(xué)的數(shù)學(xué)物理方程課程也牽涉其中,對(duì)偏微分方程求解主要有兩種方法：解析方法和數(shù)值方法.一般情況下,嚴(yán)格求解偏微分方程是非常困難,僅僅一些具有簡(jiǎn)單邊界或者有比較嚴(yán)格物理限制的現(xiàn)象才能夠得到理論解析解.最廣泛使用的是數(shù)值解法,如有限元(FEM)和有限差分方法(FD)等.格子Boltzmann(LB)方法(也稱(chēng)作格子波爾茲曼方法)是近幾十年發(fā)展的偏微分方程模擬方法,首先由McNamara和Zanetti在1988年提出,它繼承了格子氣自動(dòng)機(jī)(Lattice Gas Automaton,L GA)的主要原理并對(duì)L GA作了改進(jìn),成功地進(jìn)行了流體力學(xué)模擬,對(duì)許多問(wèn)題建立的格子模型,取得了意想不到的效果;一些學(xué)者對(duì)Navier-Stokes方程,Euler方程,mKdV方程,對(duì)流-擴(kuò)散方程,反應(yīng)-擴(kuò)散方程,熱傳導(dǎo)方程等進(jìn)行了數(shù)值求解計(jì)算,也取得了較好的結(jié)果.Lattice Boltzmann(LB)方法是一種全新概念的數(shù)值求解方法,我們?cè)谶@里對(duì)它做一定的介紹,借以拓寬方程的求解思路,為以后的應(yīng)用提供良好的借鑒作用.

2 格子Boltzmann(LB)方法模型的建立

2.1 格子Boltzmann方法的主要思想.

這里以二維中常用的D2Q9模型為例,來(lái)介紹格子Boltzmann方法.

將流體存在的區(qū)域劃分為均勻網(wǎng)格,將流體想象成許多只有質(zhì)量沒(méi)有體積的微小粒子組成,在同一時(shí)刻同一網(wǎng)格節(jié)點(diǎn)上,粒子向九個(gè)方向運(yùn)動(dòng)(如圖1),移動(dòng)到最近的網(wǎng)格節(jié)點(diǎn).其中每個(gè)節(jié)點(diǎn)上允許一個(gè)靜止粒子存在,加上與其相鄰的有8個(gè)節(jié)點(diǎn),所以稱(chēng)為二維九點(diǎn)格子模型,記為D2Q9模型.到目前為止已建立的LBM模型有：D1Q3,D2Q9,D2Q7,D2Q13,D3Q15,D3Q18,D3Q27等(D指維數(shù),Q指粒子運(yùn)動(dòng)方向的總數(shù)).

其演化過(guò)程主要分兩個(gè)步驟：(a)遷移,粒子從一個(gè)節(jié)點(diǎn)在一個(gè)時(shí)間步長(zhǎng)內(nèi),以恒定的速度運(yùn)動(dòng)到相鄰節(jié)點(diǎn);(b)碰撞,在一個(gè)節(jié)點(diǎn)上從相鄰節(jié)點(diǎn)運(yùn)動(dòng)來(lái)的粒子發(fā)生碰撞,根據(jù)質(zhì)量、動(dòng)量和能量守恒規(guī)則改變粒子的速度,然后各個(gè)粒子又以改變后的速度遷移.這兩個(gè)步驟交替循環(huán),直到流場(chǎng)達(dá)到收斂.

圖1 二維九點(diǎn)格子模型速度方向以及網(wǎng)格劃分

設(shè)(X,t)代表在時(shí)刻t,位置X=(x,y)處的節(jié)點(diǎn),流體的密度為ρ=ρ(X,t),流體的速度為u=u(X,t),時(shí)間步長(zhǎng)為Δt,t=0,Δt,2Δt,…,mΔt,格子步長(zhǎng)為Δx,粒子遷移速率c=

(cosθi,sinθi)c,(θi=(i-5)π/2+π/4,i=5～8),粒子分布函數(shù)fi(X,t)表示節(jié)點(diǎn)(X,t)處運(yùn)動(dòng)速度為ci的粒子數(shù)量,i=0,1,…,8,則粒子分布函數(shù)演化過(guò)程如下: ,粒子離散速度c0=(0,0),ci= (cosθi,sinθi)c,(θi=(i-1)π/2,i=1～4),ci=

Ωi(f(X,t))稱(chēng)為碰撞算子,表示碰撞引起的變化.Higuera等在1989年作了一個(gè)非常重要的簡(jiǎn)化,假定流場(chǎng)的分布接近于局部平衡狀態(tài),并將碰撞操作過(guò)程Ωα(f(x,t))進(jìn)行線性化處理,用單一松弛時(shí)間使流場(chǎng)逐漸達(dá)到局部平衡狀態(tài)并且線性穩(wěn)定;后來(lái),Qian,Chen等學(xué)者提出,如果松弛形式采用Bhatnagar-Gross-Krook(BGK)操作就能再現(xiàn)宏觀N-S方程,這就是所謂的格子BGK(LBGK)模型,它所對(duì)應(yīng)的LB方程(LBE)是(也稱(chēng)為格子Boltzman模型的演化方程)：

其中τ是松弛時(shí)間尺度,控制達(dá)到平衡的速度(可根據(jù)需要進(jìn)行設(shè)置),由于穩(wěn)定性的原因,經(jīng)過(guò)實(shí)際測(cè)算τ必須大于1/2.

事實(shí)上不同的網(wǎng)格剖分有著不同的平衡分布函數(shù),LBM建立模型的核心問(wèn)題就是根據(jù)不同的網(wǎng)格確定對(duì)應(yīng)的平衡分布函數(shù)和格子Boltzman模型的演化方程,對(duì)D2Q9模型,我們就取上面的平衡分布函數(shù).

在節(jié)點(diǎn)(X,t)上根據(jù)質(zhì)量和動(dòng)量守恒規(guī)則,流體的宏觀密度、速度、壓強(qiáng)定義如下:

另外滿(mǎn)足動(dòng)量通量守恒法則

其中I為二階單位矩陣.

2.2 從Lattice Boltzmanm模型再現(xiàn)宏觀流場(chǎng)控制方程.

格子法建模的核心就是確定對(duì)應(yīng)網(wǎng)格的平衡分布函數(shù),而建模成功與否,看格子Boltzmann模型能否恢復(fù)宏觀流場(chǎng)控制方程.下面以Lattice BGK(LBGK)模型對(duì)Navier-Stokes方程(簡(jiǎn)稱(chēng)N-S方程)恢復(fù)為例加以說(shuō)明.

這里簡(jiǎn)單介紹多尺度分析Chpamna-Enskog展開(kāi),Chpamna-Enskog展開(kāi)法本質(zhì)上是一種多尺度方法,基本思想是將時(shí)空變量用多層次時(shí)空尺度來(lái)表示,而在各級(jí)尺度上,物理量的量級(jí)一致.我們進(jìn)行Chapman-enskog展開(kāi)：

假設(shè)系統(tǒng)雖然沒(méi)達(dá)到平衡態(tài),但在局部小體積中已離平衡態(tài)不遠(yuǎn),因此可假定分布函數(shù)的展開(kāi)

引入兩個(gè)宏觀時(shí)間尺度t1=εt,t2=ε2t,一個(gè)宏觀空間尺度x1=εx,則算子的展開(kāi)

(26)與(28)正是不可壓Navier-Stokes方程組(簡(jiǎn)稱(chēng)N-S方程).

2.3 格子Boltzmann方法的計(jì)算流程.

i)為流場(chǎng)內(nèi)的每一個(gè)格子設(shè)置初始條件(密度、速度、壓強(qiáng)等),選擇合適的時(shí)間松弛尺度τ;

ii)對(duì)于每個(gè)格子,用式(3)和式(4)計(jì)算它們的宏觀區(qū)域中的密度和速度的變化;

iii)用式(2)計(jì)算平衡態(tài)分布函數(shù);

iv)把格子分布函數(shù)和平衡態(tài)分布函數(shù)代入式(1);

v)根據(jù)步驟iv)的結(jié)果,對(duì)各個(gè)相鄰格子進(jìn)行相應(yīng)的調(diào)整;

vi)根據(jù)邊界條件來(lái)調(diào)整分布函數(shù);

vii)返回步驟ii)進(jìn)行循環(huán)執(zhí)行,直到流場(chǎng)達(dá)到收斂要求.

接下來(lái)我們以?xún)蓚€(gè)定解問(wèn)題為例,驗(yàn)證上面方法的正確性和可行性.

3 格子Boltzmann方法數(shù)值模擬不可壓Navier-Stokes方程組

設(shè)某流場(chǎng)流動(dòng)區(qū)域?yàn)?≤x≤2π,0≤y≤0.16π,流體的速度為u=(u(x,y,t),v(x,y,t)),流體的壓強(qiáng)p=p(x,y,t),滿(mǎn)足不可壓Navier-Stokes方程組,即·u=0,+u·u=-狆＋2u,其定解問(wèn)題如下:

其中A,B,k,p0為常數(shù),可以驗(yàn)證該問(wèn)題具有如下解析解:

用Lattice Boltzmann方法進(jìn)行數(shù)值模擬,模擬中不妨取A=B=k=1.0,p0=3.0,模擬三個(gè)不同時(shí)刻t=0.018277和t=0.18277,t=1.8277在y方向上的分量v(x,y,t)如圖2,圖上畫(huà)出v(x,y,t)在通過(guò)流場(chǎng)中心(y=0.08π)模擬結(jié)果與解析解,可以看到模擬結(jié)果與解析解在不同時(shí)刻都吻合得很好.

圖2 y方向速度分量v(x,y,t)在t=0.018277和t=0.18277,t=1.8277時(shí)刻的模擬解與解析解相比較

4 格子Boltzmann方法數(shù)值模擬含源項(xiàng)的擴(kuò)散方程

二維含源項(xiàng)對(duì)流-擴(kuò)散方程的一般形式如下:

其中u=(u1,u2)是一個(gè)2維的常向量,代表對(duì)流速度;α是擴(kuò)散系數(shù),是一常數(shù);X=(x,y)是空間的坐標(biāo),t代表時(shí)間;ρ(X,t)是物質(zhì)在空間點(diǎn)x和時(shí)刻t的密度;F(X,t)是在空間點(diǎn)x和時(shí)刻t的源項(xiàng);當(dāng)u=0時(shí),方程(29)就變?yōu)楹错?xiàng)的擴(kuò)散方程.

仍采用D2Q9模型,由于方程有改變,對(duì)應(yīng)的LB模型有所調(diào)整,這是關(guān)鍵,主要變化是格子Boltzman模型的演化方程(1).采用的演化方程為

用格子Boltzmann方法進(jìn)行數(shù)值模擬,模擬中不妨取l=1,α=0.001,模擬兩個(gè)不同時(shí)刻t=0.1和t=0.3在y=0.5處的密度ρ,如圖3,可以看到模擬結(jié)果與解析解在不同時(shí)刻都吻合得很好.

圖3 密度ρ在t=0.1和t=0.3時(shí)刻在y=0.5處的模擬解與解析解相比較

5 結(jié) 語(yǔ)

本文初步介紹了Lattice Boltzmann方法的理論,作為一種新的模擬和建模方法,其理論思想簡(jiǎn)單,入門(mén)門(mén)檻不高,具有微積分和一些場(chǎng)論基礎(chǔ)就行.格子波爾茲曼方法打破傳統(tǒng)的數(shù)學(xué)建模觀念,不需建立和求解復(fù)雜的偏微分方程,是一種高效建模數(shù)值方法,具有廣闊應(yīng)用前景,其研究已從純粹的理論試驗(yàn)領(lǐng)域邁向工程應(yīng)用,在生物流體、磁流體、交通流、燃燒、微尺度流動(dòng)、化工、圖像處理、量子力學(xué)、納米流體以及光學(xué)、聲學(xué)等眾多領(lǐng)域得到了應(yīng)用,形成了一個(gè)國(guó)際研究熱點(diǎn).目前格子Boltzmann方法還處于不斷發(fā)展之中,作為一種新生事物,人們對(duì)它還缺乏了解和比較陌生,需要大量的論證和實(shí)踐,挖掘其潛力,開(kāi)拓新的應(yīng)用領(lǐng)域,作為我們數(shù)學(xué)工作者是一個(gè)很好的切入點(diǎn),很值得我們?nèi)プ鲞M(jìn)一步的研究工作.

[1] McNamaraG R,Zanetti G.Use of the lattice Boltzmann equation to simulate lattice-gas automata[J].Phys.Rev. Lett.,1988,61：2332-2335.

[2] GuoZhaoli,Shi Baochang,Wang Nengchao.Lattice BGK model for incompressible Navier-Stokes equation[J].Journal of Computational Physics,2000,165：288-306.

[3] Deng Bin,Shi Bao-chang,Wang Guang-chao.A new lattice Bhatnagar-Gross-Krook model for the convection-diffusion equation with a source term[J].Chin.Phys.Lett,2005,22：267-270.

A Novel Method for the Partial Differential Equation to Solve—Lattice Boltzmann Model

L E L i-hua, GAO Yun, L IU Tang-wei
(School of Mathematics&Information Science,East China Institute of Technology Fuzhou,Jiangxi 344000,China)

We introduce a novel method for the partial differential equation to solve-Lattice Boltzmann model,The basic theory and basic principle of LBM is analyzed in detail.Two results illustrate that Lattice Boltzmann method is right,effective through numerical simulation on incompressible Navier-Stokes equation and the 2-dimensional diffusion equation with the source term,that shows a wide application prospect to us,which probably provides reference for intensive research and extensive applications later on.

lattice Boltzmann method;equilibrium distribution functions;D2Q9 model;Navier-Stokes equation; convection-diffusion equation

O241.82

A

1672-1454(2011)03-0075-08

2008-08-28

國(guó)家自然科學(xué)基金項(xiàng)目“復(fù)雜流動(dòng)的格子Boltzmann建模與計(jì)算機(jī)仿真”(60773195);中國(guó)科學(xué)院邊緣海地質(zhì)重點(diǎn)實(shí)驗(yàn)室開(kāi)放研究基金課題“孔隙裂隙介質(zhì)中多相流體逾滲模型及數(shù)值模擬”(2010.1—2011.12)