運用格子Boltzmann方法求解空間分數階方程

李 康, 張建影

(長春工業大學 數學與統計學院, 吉林 長春 130012)

0 引 言

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

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

文中主要采用格子Boltzmann方法,對擴散項含有分數階導數的情況進行處理后,求解如下形式的一維Riesz空間分數階偏微分方程:

(1)

式中:α----Riesz空間分數階導數的階數,1<α<2;

kα----擴散系數;

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

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

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

1 Riesz空間分數階Navier-Stokes方程的預處理

首先,需要將方程中Riesz空間分數階導數項進行處理。文中采取將Riesz空間分數階導數進行離散化處理的方法。

Riesz空間分數階導數的定義式為

(2)

其中

(3)

(4)

分別是左側和右側的Riemann-Liouville分數階導數。

方程(1)中Riesz空間分數階導數項為

(5)

根據Riemann-Liouville分數階導數定義,利用數值積分公式可將式(5)化簡為

(6)

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

(7)

令

則式(6)中左積分

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

(8)

令

則式(6)中右積分

通過以上處理,方程(1)中的擴散項可以變為

(9)

再令

則方程(1)可以變為

(10)

(11)

式中:H----源項;

β----擴散系數。

2 構造格子Boltzmann模型

(12)

格子Boltzmann方程可以寫成

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

(13)

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

L----特征長度。

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

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

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

(14)

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

(15)

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

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

(16)

并且

(17)

同時假設源項為

(18)

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

(19)

(20)

其中,系數Ci為

偏微分算子

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

可得

(21)

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

(22)

(23)

(24)

則可以推出平衡分布函數為

(25)

(26)

β=-ελC2,

(27)

式(21)可以變為

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

(28)

其中β=-λεC2是擴散系數,

我們選擇

所以

3 數值模擬

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

算例1 兩平行板間的壓力驅動流動(Riesz空間分數階Navier-Stoke方程)

(29)

這個方程的精確解是

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

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

(30)

(31)

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

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

圖1 數值解與精確解對比

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

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

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

圖2 算法收斂階的對數圖像

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

(a) 改變網格數m

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

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

算例2 Riesz空間分數階擴散方程

(32)

這個方程的解析解為

其中bn的表達式為

數值解與解析解對比如圖4所示。

圖4 數值解與解析解對比

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

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

(a) 改變算法中的網格數m

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

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

4 結 語

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