用格子Boltzmann方法求解修正的時間分數階方程

劉 鑫,張建影

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

分數階微積分方程是在傳統微積分方程的基礎上,用分數階導數(積分)項代替傳統微積分方程的時間或空間導數(積分)項.與整數階方程相比,分數階微積分方程在描述一些特殊問題時更有意義.例如: 多孔介質滲流問題[1];非牛頓流體的非局域性問題[2];描述地下水在非飽和土中的滲流過程[3];將分數階方程與最優控制結合起來開發高效算法[4]等.但由于分數階算子的特殊性,使得尋找分數階微積分方程的精確解析解較困難.相比于求解精確解析解,數值解法更靈活,也適用于更多情況.目前關于分數階微積分方程數值解的研究已有一定進展,例如: Mukhtar等[5]用Adomian分解變換方法和q-同倫分析變換方法,得到了分數階方程多維模型的數值結果;Ali等[6]用最優同倫漸近方法推導出Caputo型分數階導數的半解析解;Burqan等[7]利用Laplace變換和殘差冪級數方法求解了分數階方程的級數解.

時間分數階方程是分數階方程的重要分支,如何開發其數值算法目前已得到廣泛關注.相比于宏觀的有限差分方法,在求解時間分數階微分方程時使用的復雜格式降低了算法的有效性.介于統計力學與流體力學之間的人工系統——格子Boltzmann模型,其更簡單并適用于數值模擬.本文用格子Boltzmann方法(LBM)求解修正的時間序列分數階方程.

1 基本理論

1.1 修正的時間分數階擴散方程

考慮一個修正的時間分數階擴散方程[8]:

(1)

(2)

1.2 對修正的時間分數階擴散方程進行離散

整理可得

(4)

其中

1.3 格子Boltzmann模型

下面用格子Boltzmann方法恢復宏觀方程.定義

(6)

這里fα(x,t)表示在時間為t、位置為x、速度為eα的粒子的分布函數,α為粒子速度方向.格子Boltzmann方程為

(7)

(8)

引入小Knudsen數ε,定義ε=l/L=Δt,其中l為平均自由程,L為宏觀特征長度.在引入ε的情況下,對分布函數fα(x,t)進行Chapman-Enskog展開,有

(9)

(10)

進行時間多尺度展開可得

t=t0+εt1+ε2t2+ε3t3+ε4t4+O(ε5),

(11)

從而有

(12)

源項有如下表達形式:

Ωα=εΩ1+ε2Ω2+ε3Ω3+ε4Ω4+O(ε5).

(13)

將式(8),(9),(11)～(13)代入式(7),并比較ε的各階系數有

(16)

設源項滿足如下條件:

平衡態分布函數滿足:

(19)

(20)

(21)

(22)

根據式(17)～(21)對式(22)進行整理,即可恢復出宏觀方程

(23)

圖1 D1Q3的格子模型Fig.1 Lattice model of D1Q3

這里＝－Dλεc2

下面推導平衡態分布函數.D1Q3的格子模型如圖1所示.

圖2 D2Q9的格子模型Fig.2 Lattice model of D2Q9

根據式(19)～(21)可得平衡態分布函數為

(24)

(25)

其中c=|eα|是粒子沿各方向的移動速度.D2Q9的格子模型如圖2所示.

同理可得平衡態分布函數為

2 數值計算

下面計算一個一維實例,因此所涉及的矢量x這里均為標量x.定義uN(xk,t)為數值解,uE(xk,t)為精確解.選擇相對誤差、最大絕對誤差和全局相對誤差進行驗證,各誤差表達式分別如下:

(29)

(30)

(31)

例1考慮帶有如下初邊值條件的時間分數階方程:

(32)

方程(32)的源項為

(33)

其精確解為

u(x,t)=t2ex.

(34)

圖3 例1的數值計算結果(A)、誤差分析(B)和收斂階(C)Fig.3 Numerical calculation results (A),error analysis (B) and convergence order (C) of example 1

圖3給出了例1的數值計算結果、誤差分析和收斂階.由圖3(A)可見,數值解和精確解曲線吻合度較高,幾乎完全重疊;由圖3(B)可見,最大絕對誤差為0.009 10,最大相對誤差為0.004 52,兩種誤差都較小;由圖3(C)可見,在一維情形下格子Boltzmann模型具有空間一階收斂階.

(35)

相對誤差、最大絕對誤差和全局相對誤差的表達式如下:

例2考慮帶有如下初邊值條件的時間分數階方程:

(39)

方程(39)的源項為

(40)

精確解為

u(x,y,t)=t2ex+y.

(41)

數值計算的參數選取: 時間步長和空間步長分別為

格子數為65×65,時間t=1.0,松弛因子τ=1.25,分數階取值為β=0.001,α=0.001,擴散系數為D=0.001.

圖4 例2的數值計算結果(A),(B)、誤差分析(C)和收斂階(D)Fig.4 Numerical calculation results (A),(B),error analysis (C) and convergence order (D) of example 2

圖4給出了例2的二維數值計算結果、誤差分析和收斂階.由圖4(A),(B)可見,方程(39)的數值解和精確解基本一致;由圖4(C)可見,最大相對誤差是0.012 2,最大絕對誤差是0.009 06;由圖4(D)可見,在二維情形下格子Boltzmann模型具有空間一階精度.

綜上所述,本文基于格子Boltzmann方法,針對修正的時間分數階方程進行了討論.首先將方程離散化處理,得到了標準的反應擴散方程.其次,根據Taylor展開式和Chapman-Enskog多尺度展開等技術,用格子Boltzmann方法準確地恢復了宏觀方程,并推導出各方向的平衡態分布函數表達式.最后,用一個一維實例和一個二維實例進行數值計算,所得數值解與精確解吻合較好,誤差在合理范圍內.因此,格子Boltzmann模型在求解修正的時間分數階方程的數值解上效果較好.