運用格子Boltzmann方法求解空間分?jǐn)?shù)階方程

李 康, 張建影

(長春工業(yè)大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院, 吉林 長春 130012)

0 引 言

近年來,分?jǐn)?shù)階偏微分方程受到越來越廣泛的關(guān)注,學(xué)者對于不同方程的數(shù)值解進(jìn)行了大量研究,如分?jǐn)?shù)階擴散方程[1-6]、分?jǐn)?shù)階Navier-Stokes方程[7]、分?jǐn)?shù)階平流擴散方程[8-11]、分?jǐn)?shù)階反應(yīng)擴散方程[12]等。對于導(dǎo)數(shù)項是分?jǐn)?shù)階的偏微分方程,一般來說,求解方法主要為數(shù)值求解,因為只有極少數(shù)特殊的分?jǐn)?shù)階偏微分方程可以得到解析解[11],應(yīng)用范圍并不廣泛,因此主要通過數(shù)值方法得到方程的近似解。

對于計算分?jǐn)?shù)階偏微分方程,部分學(xué)者采用有限差分法[1,2,9],也有部分學(xué)者用有限元法[5]、有限體積法[4]。近年來,由于格子Boltzmann方法在求解其他的線性和非線性偏微分方程[13]取得了滿意的成果,所以部分學(xué)者開始對格子Boltzmann方法求解分?jǐn)?shù)階偏微分方程進(jìn)行研究[14-18]。從目前研究結(jié)果看,用格子Boltzmann方法構(gòu)建的模型可以準(zhǔn)確恢復(fù)出宏觀連續(xù)方程,并且通過數(shù)值模擬驗證了數(shù)值計算的準(zhǔn)確性和有效性,所以LBM在分?jǐn)?shù)階微分方程上的應(yīng)用還具有廣闊前景。

文中主要采用格子Boltzmann方法,對擴散項含有分?jǐn)?shù)階導(dǎo)數(shù)的情況進(jìn)行處理后,求解如下形式的一維Riesz空間分?jǐn)?shù)階偏微分方程:

(1)

式中:α----Riesz空間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù),1<α<2;

kα----擴散系數(shù);

s(y,t)----源項;

g(y)----初始條件;

h1(t),h2(t)----邊界條件。

1 Riesz空間分?jǐn)?shù)階Navier-Stokes方程的預(yù)處理

首先,需要將方程中Riesz空間分?jǐn)?shù)階導(dǎo)數(shù)項進(jìn)行處理。文中采取將Riesz空間分?jǐn)?shù)階導(dǎo)數(shù)進(jìn)行離散化處理的方法。

Riesz空間分?jǐn)?shù)階導(dǎo)數(shù)的定義式為

(2)

其中

(3)

(4)

分別是左側(cè)和右側(cè)的Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)。

方程(1)中Riesz空間分?jǐn)?shù)階導(dǎo)數(shù)項為

(5)

根據(jù)Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)定義,利用數(shù)值積分公式可將式(5)化簡為

(6)

式(6)中左積分又可以化簡為

(7)

令

則式(6)中左積分

同理,式(6)中右積分也可化簡為

(8)

令

則式(6)中右積分

通過以上處理,方程(1)中的擴散項可以變?yōu)?/p>

(9)

再令

則方程(1)可以變?yōu)?/p>

(10)

(11)

式中:H----源項;

β----擴散系數(shù)。

2 構(gòu)造格子Boltzmann模型

(12)

格子Boltzmann方程可以寫成

fα(y+εeα,t+ε)=fα(y,t)-

(13)

式中:τ----單松弛時間因子,是一個常數(shù);

L----特征長度。

我們假設(shè)在數(shù)值模擬中,克努森數(shù)ε等于時間步長Δt。

對式(13)進(jìn)行Taylor展開,可以得到

fα(y+εeα,t+ε)-fα(y,t)=

(14)

在小的克努森數(shù)ε的假設(shè)下,對上述式子中的fα(y,t)進(jìn)行Chapman-Enskog展開

(15)

為了討論不同時間尺度的變化,采用多尺度分析,引入t0,t1,t2,t3作為時間尺度。

tn=εnt,n=0,1,2,3,

(16)

并且

(17)

同時假設(shè)源項為

(18)

將上述展開式(12)、式(13)、式(15)、式(16)都代入格子Boltzmann方程中,可以得到關(guān)于ε階的方程

(19)

(20)

其中,系數(shù)Ci為

偏微分算子

將式(17)+式(18)×ε,并且令

可得

(21)

為了恢復(fù)出宏觀方程(1),我們假設(shè)平衡分布函數(shù)滿足以下條件

(22)

(23)

(24)

則可以推出平衡分布函數(shù)為

(25)

(26)

β=-ελC2,

(27)

式(21)可以變?yōu)?/p>

（二）從材料包含的知識深度上看，要做到“三個思考”，即是什么，為什么，怎么樣。深層次解讀信息，是什么（什么現(xiàn)象、問題的實質(zhì)）、為什么（原因、作用、危害等）、怎么樣（措施、建議、態(tài)度等）。當(dāng)然，不一定每題都同時回答三個層次，通常用在認(rèn)識、評價類試題，要根據(jù)問題并結(jié)合材料做到具體問題具體分析。

(28)

其中β=-λεC2是擴散系數(shù),

我們選擇

所以

3 數(shù)值模擬

為了說明該算法的有效性,我們用下面的數(shù)值算例去驗證。

算例1 兩平行板間的壓力驅(qū)動流動(Riesz空間分?jǐn)?shù)階Navier-Stoke方程)

(29)

這個方程的精確解是

u(y,t)=y2(1-y)2tα。

我們給出全局相對誤差(R)和最大絕對誤差(E)的定義

(30)

(31)

式中:u(yi,t)----在(yi,t)處u的數(shù)值解;

u*(yi,t)----在(yi,t)處u的精確解。

圖1 數(shù)值解與精確解對比

圖中空間步長h=0.02,時間步長Δt=0.025,空間分?jǐn)?shù)階α=1.6,Re=1 000,時間t=2時,數(shù)值解與精確解的對比。

從圖1可以看出,數(shù)值解與精確解比較吻合,該算法適用于Riesz空間分?jǐn)?shù)階方程。

該算例收斂階的對數(shù)圖像如圖2所示。

圖2 算法收斂階的對數(shù)圖像

圖2中橫坐標(biāo)是空間步長h的對數(shù),縱坐標(biāo)表示全局相對誤差R的對數(shù)。經(jīng)過數(shù)據(jù)擬合后得出的收斂階為1.286。改變參數(shù)時數(shù)值解的對比如圖3所示。

(a) 改變網(wǎng)格數(shù)m

從圖3可以看出,改變網(wǎng)格數(shù)m以及時間步長Δt數(shù)值解并沒有明顯變化。與傳統(tǒng)的數(shù)值方法(有限差分法、有限元法)求解Riesz空間分?jǐn)?shù)階方程相比,格子Boltzmann方法求解方程時物理意義更加清晰,邊界條件易于處理,可以模擬各種復(fù)雜的宏觀現(xiàn)象,同時編程比較容易,計算前后處理也比較簡單,避免了冗長復(fù)雜的網(wǎng)格劃分過程,可以提高計算效率。

圖3(a)中m取值分別為40,60,80,100,其余參數(shù)取值不變。圖(b)中Δt取值分別為0.01,0.02,0.05,其余參數(shù)取值不變,(b)中dt表示為Δt。

算例2 Riesz空間分?jǐn)?shù)階擴散方程

(32)

這個方程的解析解為

其中bn的表達(dá)式為

數(shù)值解與解析解對比如圖4所示。

圖4 數(shù)值解與解析解對比

圖4是網(wǎng)格數(shù)m為50,時間步長Δt為0.015,時間t為0.4,擴散系數(shù)kα為0.01,空間分?jǐn)?shù)階α為1.8時,通過該方法所得的數(shù)值解(RJ)與解析解(R)的對比圖,從圖中可以看出,該方法得到的數(shù)值解與解析解比較吻合,適用于求解此方程。

改變參數(shù)時,數(shù)值解的對比如圖5所示。

(a) 改變算法中的網(wǎng)格數(shù)m

圖5(a)中m取值分別為30,40,50,其余參數(shù)取值不變。圖5(b)中Δt取值分別為0.005,0.015,0.025,0.050,其余參數(shù)取值不變,(b)中dt表示為Δt。

從圖中可以看出,改變網(wǎng)格數(shù)m以及時間步長Δt數(shù)值解并沒有明顯變化。

4 結(jié) 語

通過格子Boltzmann方法求解一維Riesz空間分?jǐn)?shù)階偏微分方程,通過對Riesz空間分?jǐn)?shù)階進(jìn)行離散化處理,得到近似標(biāo)準(zhǔn)的擴散方程,并且通過LBM的D1Q3模型,恢復(fù)出宏觀方程,同時計算出平衡態(tài)分布函數(shù)。通過兩個帶有解析解的Riesz空間分?jǐn)?shù)階偏微分方程進(jìn)行數(shù)值驗證,對比發(fā)現(xiàn)數(shù)值解與解析解能夠很好吻合。格子Boltzmann方法求解Riesz空間分?jǐn)?shù)階方程時物理意義更加清晰,邊界條件易于處理,編程計算過程比較簡單,避免了冗長復(fù)雜的網(wǎng)格劃分過程,可以提高計算效率。