王博宇,閆廣武

(吉林大學(xué) 數(shù)學(xué)學(xué)院,長春 130012)

?

用于求解Poisson方程的格子Boltzmann模型

王博宇,閆廣武

(吉林大學(xué) 數(shù)學(xué)學(xué)院,長春 130012)

提出一個求解Poisson方程的格子Boltzmann模型.通過使用Chapman-Enskog展開和多尺度展開得到了在不同時間尺度下的系列偏微分方程及平衡態(tài)分布函數(shù)和具有三階截斷的誤差修正Poisson方程.用該模型計算Kolmogorov流和Green-Taylor渦流,并與解析解進(jìn)行比較,計算結(jié)果表明,數(shù)值結(jié)果與經(jīng)典解析結(jié)果基本相符.

格子Boltzmann模型；Poisson方程；Kolmogorov流

0 引 言

格子Boltzmann方法(LBM)目前已經(jīng)發(fā)展成為用于計算流體動力學(xué)(CFD)的可選擇性方法[1-2].LBM可用于模擬單組分水動力學(xué)問題,包括懸浮微粒的多組分流體問題、磁流體力學(xué)、反應(yīng)擴散系統(tǒng)、多孔介質(zhì)流動問題及其他復(fù)雜系統(tǒng).此外,格子Boltzmann模型廣泛應(yīng)用于模擬線性和非線性偏微分方程,如波動方程[3]、Burgers方程[4]、KDV方程[5]、非線性Sch?dinger方程[6]、Poisson方程[7]和復(fù)Ginzburg-Landau方程[8]等.

本文提出一種基于格子Boltzmann模型處理Poisson方程的新數(shù)值方法.其中Poisson方程為

(1)

在應(yīng)用LBM求解方程(1)時,有兩種選擇:1)使用一種與時間無關(guān)的空間多尺度技術(shù)恢復(fù)該方程；2)使用時間多尺度技術(shù)恢復(fù)該Poisson方程.本文選用第二種方法.引入方程：

(2)

其中u(x,t)滿足

(3)

則當(dāng)t→∞時,方程(2)的解恰好是方程(1)的解.

目前,模擬Poisson方程的方法主要有有限差分法、有限單元法以及最小二乘法等[9-11].新的求解技術(shù)有分辨算法和預(yù)調(diào)節(jié)技術(shù)[12-13]、Trefftz方法[14]、平面波方法[15]、非協(xié)調(diào)Galerkin方法[16]、最小二乘法[17]和無網(wǎng)格薄板樣條方法[18]等.本文提出一種一維和二維格子Boltzmann模型,并給出了數(shù)值算例.

1 格子Boltzmann模型

1.1不同時間尺度下的系列偏微分方程

定義fα(x,t)為在t時刻位于x處的單粒子分布函數(shù).考慮一維或二維格子Boltzmann模型,在該模型中,一維或二維格子上具有b+1個離散速度,除靜止粒子外,其他粒子沿路徑移動到其他相鄰格點上.如果考慮粒子到達(dá)某一個格點時發(fā)生了粒子間的碰撞,則格子Boltzmann方程為

其中:Ωα(fα(x,t))稱為碰撞算子；ωα(x,t)表示附加項.在標(biāo)準(zhǔn)格子Boltzmann模型中,ωα=0.本文選取

其中:ε為Knudsen數(shù);φα為獨立于α的函數(shù).一種簡單碰撞算子是選取一個線性形式的具有弛豫時間τ的算符：

(7)

Knudsen數(shù)ε定義為ε=l/L,其中:l是粒子自由程;L是特征長度,可以被選為時間步長Δt.因此,格子Boltzmann方程(4)可寫為

方程(8)中,Knudsen數(shù)ε被假設(shè)為小量,故可在小Knudsen數(shù)假設(shè)下進(jìn)行Chapman-Enskog展開：

(9)

t0,t1,t2表示不同的時間尺度,定義為

(10)

及

(11)

因此,可得一系列偏微分方程：

(12)

(13)

(14)

(15)

方程(12)～(15)又稱為不同時間尺度下的系列偏微分方程,它適用于一維、二維及三維情形,其中3個多項式系數(shù)與弛豫時間τ有關(guān):

(16)

(17)

(18)

式(16)～(18)為第二～第四個Chapman多項式,與文獻(xiàn)[19]結(jié)果完全一致.

1.2Poisson方程的恢復(fù)

(19)

(20)

(21)

(22)

(23)

特別地,對于一維情形,這些矩函數(shù)可記為

(24)

(25)

(26)

(27)

(28)

圖1 計算格子示意圖Fig.1 Diagrammatic sketches of lattice

對于一維模型,考慮5-bit格子.離散化的速度向量為e=(0,c,-c,2c,-2c),其中:α=0,1,2,3,4;c表示速度.格子的示意圖如圖1所示.5-bit格子模型的平衡態(tài)分布函數(shù)由式(24)～(28)解出,分別為

(29)

(30)

(31)

(32)

(33)

1.3一維Poisson方程的格子Boltzmann模型

一維Poisson方程具有如下形式：

(34)

選取

(35)

(36)

(37)

(38)

(39)

將式(12)+ε×式(13)+ε2×式(14)+ε3×式(15),并對α求和,可得：

如果選取(b+1)εφ=μ(ku-g),則修正后的一維Poisson方程為

方程(41)中,E2和E3由下式確定：

(42)

(43)

1.4二維Poisson方程

二維Poisson方程具有如下形式：

使用如圖1(B)所示的二維格子,則粒子速度向量為

eα={(0,0),(c,0),(0,c),(-c,0),(0,-c),(2c,0),(0,2c),(-2c,0),(0,-2c)},α=0,1,…,8.

假設(shè)各階矩為式(35)及

(45)

(46)

(47)

(48)

則平衡態(tài)分布函數(shù)記為

(49)

(51)

將式(12)+ε×式(13)+ε2×式(14)+ε3×式(15)并對α求和,得

如果選取(b+1)εφα=μ(ku-g)作為式(51)的源項,則修正后的二維Poisson方程為

其中:E2為式(42);

(54)

因此,可得具有三階截斷誤差精度的Poisson方程：

(55)

2 數(shù)值算例

下面應(yīng)用格子Boltzmann方法模擬一維和二維Poisson方程.

例1考慮Kolmogorov流:

(56)

邊界條件為

如果u獨立于y,則方程具有解析解：

(58)

其數(shù)值計算結(jié)果如圖2和圖3所示.其中:圖2(A)是LBM計算結(jié)果；圖2(B)是解析解；圖3(A)為LBM結(jié)果與解析解在y=0.2π處的比較；圖3(B)為y=0.2π處的絕對誤差曲線；圖3(C)為lg Errmax對lgε在y=0.2π處的無窮模曲線.

參數(shù)A0=1.0,c=168.0,μ=1.0,λ=c2;網(wǎng)格數(shù)M=101,Δ x=π/M=π/101,Δ t=ε=Δ x/c;迭代步數(shù)10 000.

參數(shù)c=168.0,μ=1.0,λ=0.43c2;網(wǎng)格數(shù)M=101×101,Δ x=Δ y=π/M=π/101,Δ t=ε=Δ x/c=Δ y/c.

例2考慮Green-Taylor渦流,其中流函數(shù)滿足Poisson方程：

2u(x,y)+8π2u(x,y)=0, 0≤x≤1, 0≤y≤1;

(59)

邊界條件為

(60)

方程(59)具有解析解：

(61)

其解析解和數(shù)值解如圖4所示.其中圖4(A)和圖4(B)分別為LBM數(shù)值解和解析解的表面圖.由圖4可見,LBM數(shù)值解與經(jīng)典解析結(jié)果基本吻合.圖5為絕對誤差最大值Err=|uN-uA|和Knudsen數(shù)ε的關(guān)系曲線.由圖5可見,LBM數(shù)值解的精度與網(wǎng)格密度成正比,這與模型的選取有關(guān).

參數(shù)A0=1.0,c=168.0,μ=1.0,λ=c2;格子規(guī)模M=101×101,Δ x=Δ y=π/M=π/101,Δ t=ε=Δ x/c=Δ y/c;迭代步數(shù)為10 000.

參數(shù)c=168.0,A0=1.0,μ=1.0,λ=0.43c2.

綜上,本文提出了求解Poisson方程具有高階精度的格子Boltzmann模型.通過使用Chapman-Enskog展開和多尺度展開技術(shù)得到了系列偏微分方程及平衡態(tài)分布函數(shù)的高階矩和修正后的具有三階截斷誤差的Poisson方程;給出了LBM數(shù)值解與Kolmogorov流和Green-Taylor渦流精確解的比較,結(jié)果表明,數(shù)值結(jié)果與經(jīng)典解析解基本吻合.

[1] CHEN Shiyi,Doolen G D.Lattice Boltzmann Method for Fluid Flows [J].Annu Fluid Mech,1998,30:329-364.

[2] Succi S.The Lattice Boltzmann Equation for Fluid Dynamics and Beyond [M].New York:Oxford University Press,2001.

[3] ZHANG Jianying,YAN Guangwu,SHI Xiubo.Lattice Boltzmann Model for Wave Propagation [J].Phys Rev E,2009,80:026706.

[4] Velivelli A C,Bryden K M.Parallel Performance and Accuracy of Lattice Boltzmann and Traditional Finite Difference Methods for Solving the Unsteady Two-Dimensional Burger’s Equation [J].Physica A:Statistical Mechanics and Its Applications,2006,362(1):139-145.

[5] YAN Guangwu,ZHANG Jianying.A Higher-Order Moment Method of the Lattice Boltzmann Model for the Korteweg-De Vries Equation [J].Math Comput Simulat,2009,79(5):1554-1565.

[6] ZHANG Jianying,YAN Guangwu.A Lattice Boltzmann Model for the Nonlinear Schr?dinger Equation [J].J Phys A:Math and Theore,2007,40(33):10393-10405.

[7] WANG Huimin,YAN Guangwu,YAN Bo.Lattice Boltzmann Model Based on the Rebuilding-Divergency Method for the Laplace Equation and the Poisson Equation [J].J Sci Comput,2011,46(3):470-484.

[8] ZHANG Jianying,YAN Guangwu.Lattice Boltzmann Model for the Complex Ginzburg-Landau Equation [J].Phys Rev E,2010,81:066705.

[9] Nabavi M,Kamran Siddiqui M H,Dargahi J.A New 9-Point Sixth-Order Accurate Compact Finite-Difference Method for the Helmholtz Equation [J].Journal of Sound and Vibration,2007,307(3):972-982.

[10] Babuska I,Osborn J.Handbook of Numerical Analysis [M].Vol.Ⅱ.Eigenvalue Problems.Amsterdam:North-Holland,1991:641-787.

[11] Chang C L.A Least Squares Finite Element Method for the Helmholtz Equation [J].Comput Methods Appl Mech Engrg,1990,83:1-7.

[12] Farhat C,Hetmaniuk U.A Fictious Domain Decomposition Method for the Solution of Partially Axisymmetric Acoustic Scattering Problems.Part Ⅰ:Dirichlet Boundary Conditions [J].Int J Numer Methods Engrg,2002,54(9):1309-1332.

[13] Farhat C,Macedo A,Lesoinne M.A Two-Level Domain Decomposition Method for the Iterative Solution of High Frequency Exterior Helmholtz Problems [J].Numer Math,2000,85(2):283-308.

[14] Li Z C.The Trefftz Method for the Helmholtz Equation with Degeneracy [J].Applied Numerical Mathematics,2008,58(2):131-159.

[15] JIN Bangti,Marin L.The Plane Wave Method for Inverse Problems Associated with Helmholtz Type Equations [J].Engineering Analysis with Boundary Elements,2008,32(3):223-240.

[16] Douglas J,Jr,Santos J E,Shenn D.Nonconforming Galerkin Methods for the Helmholtz Equation [J].Numer Methods Partial Differential Equations,2001,17(5):475-494.

[17] Monk P,WANG Daqing.A Least-Squares Method for the Helmholtz Equation [J].Comput Methods Appl Mech Engrg,1999,175(1/2):121-136.

[18] Bouhamidi A,Jbilou K.Meshless Thin Plate Spline Methods for the Modified Helmholtz Equation [J].Comput Methods Appl Mech Engrg,2008,197(45/46/47/48):3733-3741.

[19] Holdych D J,Noble D R,Georgiadis J G,et al.Truncation Error Analysis of Lattice Boltzmann Methods [J].J Comput Phys,2004,193(2):595-619.

(責(zé)任編輯：趙立芹)

LatticeBoltzmannModelforPoissonEquation

WANG Boyu,YAN Guangwu

(CollegeofMathematics,JilinUniversity,Changchun130012,China)

A lattice Boltzmann model for the Poisson equation was proposed.By means of the Chapman-Enskog expansion and the multi-scale time expansion,a series of partial differential equations in different time scales was obtained.Moreover,the equilibrium distribution functions and the modified partial differential equation of the Poisson equation with the third-order truncation error were obtained.In numerical examples,the Kolmogorov flow and Green-Taylor vortex flow were simulated,and the comparison between numerical results of the lattice Boltzmann models and exact solutions were given.The numerical results are acceptable.

lattice Boltzmann model;Poisson equation;Kolmogorov flow

10.13413/j.cnki.jdxblxb.2015.03.11

2014-09-29.

王博宇(1992—),男,漢族,碩士研究生,從事流體力學(xué)與復(fù)雜系統(tǒng)數(shù)值模擬的研究,E-mail:bywang14@mails.jlu.edu.cn.通信作者：閆廣武(1964—),男,漢族,博士,教授,博士生導(dǎo)師,從事復(fù)雜系統(tǒng)數(shù)值模擬的研究,E-mail:yangw_jlu@126.com.

國家自然科學(xué)基金(批準(zhǔn)號：11272133).

O351.2

：A

：1671-5489(2015)03-0407-07