空間分?jǐn)?shù)階對(duì)流方程的格子Boltzmann方法*

魏雪丹,李夢(mèng)軍,戴厚平

(吉首大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,湖南 吉首 416000)

對(duì)流是由流體粒子的整體運(yùn)動(dòng)和流體分子的隨機(jī)運(yùn)動(dòng)形成的熱量、質(zhì)量和動(dòng)量傳遞[1],經(jīng)典的對(duì)流擴(kuò)散問(wèn)題在不考慮擴(kuò)散過(guò)程時(shí)是純對(duì)流問(wèn)題.目前,對(duì)流方程廣泛應(yīng)用于泥沙運(yùn)移、污染物輸送等實(shí)際問(wèn)題中[2-5].但是在復(fù)雜不均勻介質(zhì)中,由于不同介質(zhì)之間的相互作用,傳統(tǒng)的整數(shù)階對(duì)流方程已不能準(zhǔn)確刻畫(huà)復(fù)雜流體的演化過(guò)程.隨著分?jǐn)?shù)階微分理論的發(fā)展,人們逐漸發(fā)現(xiàn)分?jǐn)?shù)階微積分算子的非局部性能夠很好地刻畫(huà)具有時(shí)間記憶和全局依賴特點(diǎn)的演化過(guò)程,于是分?jǐn)?shù)階對(duì)流方程應(yīng)運(yùn)而生.然而,也正是分?jǐn)?shù)階微分算子的全局相關(guān)性,使得大多數(shù)分?jǐn)?shù)階微分方程的精確解不能顯式表達(dá),從而數(shù)值解逼近精確解在分?jǐn)?shù)階微分方程求解方面越來(lái)越重要[6-11].對(duì)于空間Caputo型分?jǐn)?shù)階微分方程,陳雪娟等[12]運(yùn)用二次多項(xiàng)式樣條函數(shù)數(shù)值求解了一類空間分?jǐn)?shù)階Fisher方程.為了進(jìn)一步拓展空間Caputo型分?jǐn)?shù)階微分方程的數(shù)值解法,筆者擬考慮一類分?jǐn)?shù)階對(duì)流方程,利用積分中值定理和線性插值方法將分?jǐn)?shù)階對(duì)流方程轉(zhuǎn)化為標(biāo)準(zhǔn)對(duì)流方程的形式,運(yùn)用格子Boltzmann方法(Lattice Boltzmann Method,LBM)進(jìn)行數(shù)值求解.

1 分?jǐn)?shù)階對(duì)流方程的格子Boltzmann方法

1.1 分?jǐn)?shù)階對(duì)流方程

考慮如下形式的分?jǐn)?shù)階對(duì)流方程：

(1)

事實(shí)上,由于Caputo型分?jǐn)?shù)階導(dǎo)數(shù)與Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)二者可以相互轉(zhuǎn)化,因此方程(1)也可以表示為

1.2 積分離散化

1.3 格子Boltzmann方法

分布函數(shù)采用如下形式的演化方程：

(2)

(3)

根據(jù)Chapman-Enskog多尺度展開(kāi)分析,引入二階時(shí)間尺度和一階空間尺度,并將g(X,t)展開(kāi)到二階,可得

(4)

(5)

g(X,t)=ε2g(2)(X,t).

(6)

在小Knudsen數(shù)ε的假設(shè)下展開(kāi)fi(X,t),可得

(7)

將方程(3)～(7)代入方程(2),可得

(8)

比較方程(8)兩邊ε的各階系數(shù),ε的系數(shù)

(9)

ε2的系數(shù)

(10)

(11)

將方程(11)代入方程(10),可得

(12)

(13)

(14)

(15)

將方程(9)和方程(12)在i方向上求和,并恢復(fù)到原尺度,可得

對(duì)于一維分?jǐn)?shù)階對(duì)流方程,采用D1Q3 LB模型求解,由方程(13)～(15)可求得平衡態(tài)分布函數(shù)的表達(dá)式為

對(duì)于二維分?jǐn)?shù)階對(duì)流方程,采用D2Q5 LB模型求解,由方程(13)～(15)可求得平衡態(tài)分布函數(shù)的表達(dá)式為

2 數(shù)值算例

為了量化格子Boltzmann模型的精確性,引入誤差檢驗(yàn)公式,全局相對(duì)誤差(Global Relative Error，用EGRE表示)和最大誤差(Maximum Error，用EME表示)分別定義為

問(wèn)題1在有限區(qū)間[0,1]上考慮方程(1)的一維初邊值問(wèn)題，其初邊值條件為

源項(xiàng)為

驗(yàn)證該一維初邊值問(wèn)題的精確解u(x,t)=x2(1-x)2et.

取參數(shù)k=0.001,τ=1.95,Δt=0.001,N=64,此時(shí)全局相對(duì)誤差見(jiàn)表1.由表1可知：相同時(shí)間T,隨著分?jǐn)?shù)階α的增加,全局相對(duì)誤差逐漸增加；相同分?jǐn)?shù)階α,隨著時(shí)間T的增加,由于誤差的累計(jì),全局相對(duì)誤差逐漸增加,但誤差范圍基本保持在10-3以內(nèi).

表1 問(wèn)題1在N=64時(shí)的全局相對(duì)誤差

當(dāng)T=1,Δt=0.001時(shí)，不同分?jǐn)?shù)階α對(duì)應(yīng)的空間最大誤差及其收斂階見(jiàn)表2.由表2可知,LBM在空間方向收斂.

表2 問(wèn)題1在T=1,Δt=0.001時(shí)的最大誤差及其收斂階

圖1示出了α=0.4時(shí)不同時(shí)刻的數(shù)值解與精確解.由圖1可見(jiàn),數(shù)值解與精確解吻合較好,且可以長(zhǎng)時(shí)間保持穩(wěn)定.

圖1 α=0.4時(shí)的數(shù)值解與精確解

問(wèn)題2在區(qū)間[0,1]×[0,1]上考慮方程(1)的二維初邊值問(wèn)題,其初邊值條件為

源項(xiàng)為

驗(yàn)證該二維初邊值問(wèn)題精確解u(x,y,t)=x2y2sint.

取參數(shù)k1=k2=0.001,τ=1.25,Δt=10-4,T=0.3,N=100,當(dāng)α=0.5時(shí)數(shù)值解與精確解如圖2所示,數(shù)值解與精確解的絕對(duì)誤差如圖3所示.由圖2,3可見(jiàn),數(shù)值解與精確解的全局相對(duì)誤差為7.888 8×10-4,最大誤差為6.188 1×10-4,二者在任意時(shí)刻和位置的誤差達(dá)到10-4數(shù)量級(jí).

圖2 α=0.5時(shí)的數(shù)值解與精確解

圖3 數(shù)值解與精確解的絕對(duì)誤差

當(dāng)T=0.1時(shí),不同α和網(wǎng)格數(shù)下的全局相對(duì)誤差見(jiàn)表3.由表3可知：相同α?xí)r,隨著網(wǎng)格數(shù)的增加,全局相對(duì)誤差逐漸減小,說(shuō)明LBM在空間上收斂；相同網(wǎng)格數(shù)下,隨著α的增加,全局相對(duì)誤差逐漸增加,但控制在10-2以內(nèi).

表3 問(wèn)題2在T=0.1時(shí)的全局相對(duì)誤差

當(dāng)網(wǎng)格數(shù)N=32時(shí),不同時(shí)刻和α下的全局相對(duì)誤差見(jiàn)表4.由表4可知，數(shù)值解能較好逼近精確解,進(jìn)一步說(shuō)明了LBM的有效性.

表4 問(wèn)題2在不同時(shí)刻和α下的全局相對(duì)誤差

3 結(jié)語(yǔ)

運(yùn)用構(gòu)建的格子Boltzmann模型數(shù)值求解了一類空間Caputo型分?jǐn)?shù)階對(duì)流方程.運(yùn)用積分中值定理和線性插值方法將分?jǐn)?shù)階對(duì)流方程轉(zhuǎn)化為標(biāo)準(zhǔn)對(duì)流方程,使得分?jǐn)?shù)階算子產(chǎn)生的記憶部分保留在源項(xiàng)中,極大地簡(jiǎn)化了數(shù)值計(jì)算過(guò)程.數(shù)值解與精確解的實(shí)例比較結(jié)果表明,格子Boltzmann模型可以有效地求解Caputo型分?jǐn)?shù)階對(duì)流方程.本研究可以推廣到空間三維問(wèn)題及時(shí)間分?jǐn)?shù)階對(duì)流問(wèn)題,進(jìn)而拓展LBM在分?jǐn)?shù)階微分方程數(shù)值解求解方面的應(yīng)用.