史秀波,閆廣武

(1.桂林理工大學(xué) 理學(xué)院,廣西 桂林 541004； 2.吉林大學(xué) 數(shù)學(xué)學(xué)院,長春 130012)

格子Boltzmann方法(LBM)作為一種新的數(shù)值方法在計(jì)算流體力學(xué)、 非線性偏微分方程等領(lǐng)域應(yīng)用廣泛[1-6].閆廣武等[7-10]將這種方法應(yīng)用于波傳播問題,為研究其他波動(dòng)問題提供了可選擇的途徑.在科學(xué)工程領(lǐng)域,許多物理問題都可以借助變系數(shù)初邊值問題描述,這些線性和非線性模型及其解析解或數(shù)值解對應(yīng)用科學(xué)具有重要意義[11].本文用格子Boltzmann方法對如下變系數(shù)wave-like方程進(jìn)行模擬研究:

(1)

其中f(x),g(x)和h(x)是關(guān)于x,y,z的函數(shù).

方程(1)可以表示成如下形式：

(2)

這里:

Cs(x)=f(x)+g(x)+h(x);

u(x,t)的下一個(gè)時(shí)間步表達(dá)式[7]為

(3)

本文用格子Boltzmann方法為方程(2)建立了格子Boltzmann模型.通過使用Chapman-Enskog展開和多尺度技術(shù),獲得了系列格子Boltzmann偏微分方程[12]、 平衡態(tài)分布函數(shù)的高階矩以及宏觀變系數(shù)wave-like方程.數(shù)值實(shí)驗(yàn)表明,模擬結(jié)果與解析解吻合較好.

1 格子Boltzmann模型

選擇二維5-bit網(wǎng)格和三維7-bit網(wǎng)格,分布函數(shù)fα(x,t)定義為在某節(jié)點(diǎn)x上、t時(shí)刻、 具有速度eα(α=0,1,…,b)的粒子出現(xiàn)的概率,其中α=0表示靜止粒子.在二維空間中,b=4,粒子速度為

eα={(0,0),(c,0),(0,c),(-c,0),(0,-c)}；

三維空間中,b=6,粒子速度為

eα={(0,0,0),(c,0,0),(0,c,0),(0,0,c),(-c,0,0),(0,-c,0),(0,0,-c)},

其中c表示速率.定義宏觀量:

(4)

(5)

格子Boltzmann方程表示為

fα(x+eα,t+1)-fα(x,t)=Ωα+ωα,

(6)

選取Knudsen數(shù)ε作為數(shù)值模擬的時(shí)間步長和Chapman-Enskog展開的小參數(shù)[13],在該尺度上,方程(6)可寫為

fα(x+εeα,t+ε)-fα(x,t)=Ωα+ωα,

(7)

其中

ωα(x,t)=ε2θα(x,t).

(8)

通過Chapman-Enskog展開和時(shí)間多尺度可獲得不同時(shí)間尺度上的系列格子Boltzmann偏微分方程：

(9)

(10)

為了獲得平衡態(tài)分布函數(shù),選取高階矩為

(14)

當(dāng)選取

(15)

且假設(shè)附加分布函數(shù)θα和α無關(guān),并取

(16)

后,即得宏觀變系數(shù)wave-like方程(2).

結(jié)合式(5),(12),(13),可得平衡態(tài)分布函數(shù)的表達(dá)式為

其中D表示空間維數(shù).

進(jìn)行數(shù)值計(jì)算時(shí),迭代過程分為如下兩步：

1) 碰撞：

(19)

2) 流：

(20)

(21)

為了對模型進(jìn)行誤差分析,將式(9)+式(10)×ε+式(11)×ε2,并對兩端關(guān)于α求和有

(22)

由于附加分布函數(shù)θα不是一個(gè)正常尺度上的量,因此導(dǎo)致誤差出現(xiàn)反彈現(xiàn)象,使模型的精度最后降為一階.

2 數(shù)值模擬

為了驗(yàn)證模型的效果,本文分別對二維和三維變系數(shù)wave-like問題進(jìn)行數(shù)值模擬.對二維問題使用5-bit模型,三維問題使用7-bit模型.

2.1 二維初邊值問題

二維初邊值問題(initial boundary value problems,IBVP):

Neumann邊界條件為

初始條件為

u(x,y,0)=x4,ut(x,y,0)=y4;

(24)c

精確解[14]為

u(x,y,t)=x4cosht+y4sinht.

(24)d

圖1 t=0.6時(shí)初邊值問題(24)的LBM解(A)和精確解(B)Fig.1 LBM solution (A) and exact solution (B) for IBVP(24) at t=0.6

2.2 三維非齊次IBVP

三維非齊次IBVP:

Dirichlet邊界條件為

初始條件為

u(x,y,z,0)=0,ut(x,y,z,0)=x2+y2-z2;

(25)c

精確解[14]為

u(x,y,z,t)=(x2+y2)et+z2e-t-(x2+y2+z2).

(25)d

圖2 t=0.6時(shí)初邊值問題(24)在x=0.2處的相對誤差曲線(A)及 絕對誤差無窮模和Knudsen數(shù)ε的對數(shù)關(guān)系曲線(B)Fig.2 Curves of the relative error (A) and logarithmic relationship curves for IBVP(24) of the infinite norm of the absolute error Ea versus the Knudsen number ε at t=0.6,on x=0.2

參數(shù)設(shè)置為：格子尺寸50×50×50,c=5.0,τ=1.01,t=2.0.圖3(A)和(B)分別為三維非齊次IBVP(25)格子Boltzmann模型的數(shù)值解u和精確解u*.圖4(A)和(B)分別為三維非齊次IBVP(25)兩種結(jié)果在t=2時(shí)的相對誤差Er=|(u-u*)/u*|曲線及相對誤差的無窮模和Knudsen數(shù)ε的對數(shù)關(guān)系曲線.

圖3 t=2時(shí)初邊值問題(25)的LBM解(A)和精確解(B)Fig.3 LBM solution (A) and exact solution (B) for IBVP(25) at t=2

圖4 t=2時(shí)初邊值問題(25)在x=0.6,z=0.6處的相對誤差曲線(A)及 相對誤差無窮模和Knudsen數(shù)ε的對數(shù)關(guān)系曲線(B)Fig.4 Curves of the relative error (A) and logarithmic relationship curves (B) for IBVP(25) of the infinite norm of the relative error Er versus the Knudsen number ε at t=2,on x=0.6,z=0.6

由圖3和圖4可見,LBM解與精確解基本一致.由圖4(A)可見,誤差在(0,0.022 5)范圍內(nèi),表明數(shù)值解與精確解吻合較好.由圖4(B)可見,本文模型的誤差對網(wǎng)格數(shù)有依賴關(guān)系,網(wǎng)格越密,模型的誤差越小,表明LBM模型是收斂的.數(shù)值實(shí)驗(yàn)表明,LBM是用于模擬wave-like方程的一種有效方法.

[1] QIAN Yue-hong,D’Humieres D,Lallemand P.Lattice BGK Model for Navier-Stokes Equations [J].Europhysics Letters,1992,17(6): 479-484.

[2] Benzi R,Succi S,Vergassola M.The Lattice Boltzmann Equation: Theory and Applications [J].Physics Report,1992,222: 145-197.

[3] CHEN Shi-yi,Doolen G D.Lattice Boltzmann Method for Fluid Flow [J].Annual Review of Fluid Mechanics,1998,30: 329-364.

[4] Palpacelli S,Succi S.Quantum Lattice Boltzmann Simulation of Expanding Bose-Einstein Condensates in Random Potentials [J].Physical Review E,2008,77: 066708.

[5] SHI Yong,Peter L B,Ying W Y,et al.Accuracy of the Lattice Boltzmann Method for Low-Speed Noncontinuum Flows [J].Physical Review E,2011,83: 045701.

[6] Dellar P J,Lapitski D,Palpacelli S,et al.Isotropy of Three-Dimensional Quantum Lattice Boltzmann Schemes [J].Physical Review E,2011,83: 046706.

[7] YAN Guang-wu.A Lattice Boltzmann Equation for Waves [J].Journal of Computational Physics,2000,161(1): 61-69.

[8] ZHANG Jian-ying,YAN Guang-wu,SHI Xiu-bo.Lattice Boltzmann Model for Wave Propagation [J].Physical Review E,2009,80: 026706.

[9] SHI Xiu-bo,YAN Guang-wu,ZHANG Jian-ying.A Multi-energy-level Lattice Boltzmann Model for Two-Dimensional Wave Equation [J].International Journal for Numerical Methods in Fluids,2010,64(2): 148-162.

[10] ZHANG Jian-ying,YAN Guang-wu,YAN Bo,et al.A Lattice Boltzmann Model for Two-Dimensional Sound Wave [J].International Journal for Numerical Methods in Fluids,2011,67: 214-231.

[11] Allahviranloo T,Abbasbandy S,Rouhparvar H.The Exact Solutions of Fuzzy Wave-Like Equations with Variable Coefficients by a Variational Iteration Method [J].Applied Soft Computing,2011,11(2): 2186-2192.

[12] ZHANG Jian-ying,YAN Guang-wu.Lattice Boltzmann Method for One and Two- Dimensional Burgers Equation [J].Physica A: Statistical Mechanics and Its Applications,2008,387(19/20): 4771-4786.

[13] Chapman S,Cowling T G.The Mathematical Theory of Non-uniform Gas [M].Cambridge: Cambridge University Press,1939.

[14] Wazwaz A M,Gorguis A.Exact Solutions for Heat-Like and Wave-Like Equations with Variable Coefficients [J].Applied Mathematics and Computation,2004,149(1): 15-29.