用格子Boltzmann方法求解修正的时间分数阶方程

第61卷 第6期吉林大学学报(理学版)
V o l .61 N o .6
2023年11月
J o u r n a l o f J i l i nU n i v e r s i t y (
S c i e n c eE d i t i o n )N o v 2023
d o i :10.13413/j .c n k i .j
d x b l x b .2023074用格子B o l t z m a n n 方法求解
修正的时间分数阶方程
刘 鑫,张建影
(长春工业大学数学与统计学院,长春130012
)摘要:首先,根据T a y l o r 展开和C h a p m a n -E n s k o g 多尺度展开等技术,用格子B o l t z m a n n 方法准确地恢复所讨论的宏观方程,并推导D 1Q 3和D 2Q 9模型的平衡态分布函数的表达式.
其次,给出两个数值实例验证该方法的有效性.结果表明,用格子B o l t z m a n n 方法能求解
C a p
u t o 型修正的时间分数阶方程的数值解.关键词:C a p u t o 型分数阶方程;格子B o l t z m a n n 方法;反应扩散方程;分数阶方程离散化处理
中图分类号:O 242 文献标志码:A 文章编号:1671-5489(2023)06-1333-06
L a t t i c eB o l t z m a n n M e t h o d t o S o l v e
M o d i f i e dT i m eF r a c t i o n a l E q
u a t i o n L I U X i n ,Z H A N GJ i a n y i n g
(S c h o o l o f M a t h e m a t i c s a n dS t a t i s t i c s ,C h a n g c h u nU n i v e r s i t y o f T e c h n o l o g y ,C h a n g
c h u n 130012,C h i n a )A b s t r a c t :F i r s t l y ,b a s e do nt e c h n i q u e ss u c ha s T a y l o re x p a n s i o na n
d C h a p m a n -E n s k o g m u l t i -s c a l
e e x p a n s i o n ,t h e l a t t i c eB o l t z m a n n m e t h o d w a su s e dt oa c c u r a t e l y r e c o v e rt h ed i s c u s s e d m a c r o s c o p i c e q u a t i o n s ,a n dt h ee q u i l i b r i u m d i s t r i b u t i o n
f u n c t i o ne x p
r e s s i o n so fD 1Q 3a n d D 2Q 9m o d e l s w e r e d e r i v e d .S e c o n d l y ,t w on u m e r i c a le x a m p l e s w e r eu s e dt ov e r i f y t h ee f f e c t i v e n e s so ft h e p r o p o s e d m e t h o d .T h er e s u l t ss h o wt h a tt h el a t t i c eB o l t z m a n n m e t h o dc a nb eu s e dt os o l v et h en u m e r i c a l
s o l u t i o no f t h eC a p u t o t y p em o d i f i e d t i m e f r a c t i o n a l e q
u a t i o n .K e y w o r d s :C a p u t ot y p ef r a c t i o n a l e q u a t i o n ;l a t t i c eB o l t z m a n n m e t h o d ;r e a c t i o nd i f f u s i o ne q u a t i o n ;d i s c r e t i z a t i o no f f r a c t i o n a l e q
u a t i o n 收稿日期:2023-03-03.
第一作者简介:刘 鑫(1996 ),女,汉族,硕士研究生,从事计算数学的研究,E -m a i l :2459808210@q q
.c o m.通信作者简介:张建影(1982 ),女,汉族,博士,副教授,从事计算数学的研究,E -m a i l :z h a n g j y
_c c u t @126.c o m.基金项目:国家自然科学基金(批准号:11602033).
分数阶微积分方程是在传统微积分方程的基础上,用分数阶导数(积分)项代替传统微积分方程的时间或空间导数(积分)项.与整数阶方程相比,分数阶微积分方程在描述一些特殊问题时更有意义.
例如:多孔介质渗流问题[
1];非牛顿流体的非局域性问题[2];描述地下水在非饱和土中的渗流过程[3
];将分数阶方程与最优控制结合起来开发高效算法[4]
等.但由于分数阶算子的特殊性,使得寻找分数阶
微积分方程的精确解析解较困难.相比于求解精确解析解,数值解法更灵活,也适用于更多情况.目
前关于分数阶微积分方程数值解的研究已有一定进展,例如:M u k h t a r 等[5]用A d o m i a n 分解变换方法和q -
同伦分析变换方法,得到了分数阶方程多维模型的数值结果;A l i 等[6]用最优同伦渐近方法推导
出C a p u t o 型分数阶导数的半解析解;B u r q a n 等[7]利用L a p
l a c e 变换和残差幂级数方法求解了分数阶方程的级数解.
时间分数阶方程是分数阶方程的重要分支,如何开发其数值算法目前已得到广泛关注.相比于宏
观的有限差分方法,在求解时间分数阶微分方程时使用的复杂格式降低了算法的有效性.介于统计力学与流体力学之间的人工系统 格子B o l t z m a n n 模型,其更简单并适用于数值模拟.本文用格子
B o l t z m a n n 方法(L B M )
求解修正的时间序列分数阶方程.1 基本理论
1.1 修正的时间分数阶扩散方程
考虑一个修正的时间分数阶扩散方程[
8
]:∂u (x ,t )∂t =D ∂1-α∂t 1-α-∂1-β∂t 1-æèçöø÷β
∂2u (x ,t )∂x
2+g (x ,t ), 0<t ɤT , x ɪχ, αɪ(0,1), βɪ(
0,1),u (
x ,0)=0, x ɪχ,u (x ,t )=ψ(x ,t ), x ɪχìîíïïï
ï
ïï,
(1)其中χ是定义域,χ和χ分别是χ的内部和边界,ψ(
x ,t )是边界条件,g (x ,t )是源项,D 是扩散系数.式(1)是一个C a p u t o 型时间分数阶导数.根据C a p
u t o 型分数阶导数的定义[9]
,可将式(1)中关于时间的分数阶导数转化为
∂1-αu (x ,t )∂t 1-α+∂1-βu (x ,t )∂t 1-β
=1Γ(α)ʏt 0(t -ξ)α-1∂u (x ,ξ)∂ξd ξ+1Γ(β)ʏ
t 0(t -ξ)β-1∂u (x ,ξ)∂ξ
d ξ.(2)1.2 对修正的时间分数阶扩散方程进行离散
下面对式(2)中的时间进行离散.首先将时间段[0,t ]均匀划分为n 个区间,每个区间间隔为Δt =
t n
,然后将每个时间节点定义为t m (0<m ɤn ).于是t m =m
Δt ,引用线性插值对式(2)进行处理,可得1Γ(α)ʏ
t 0(t -ξ)α-1∂u (x ,ξ)∂ξd ξ+1Γ(β)ʏ
t 0(t -ξ)β-1∂u (x ,ξ)∂ξ
d ξ=Δt α-1(1+α)ðn
m =1
[u (x ,t m )-u (x ,t m -1)][(n -m +1)α-(n -m )α]+Δt β
-1Γ(1+β)ðn
m =1
[u (x ,t m )-u (x ,t m -1)][(n -m +1)β-(n -m )β],(3)整理可得
∂u (x ,t )∂t =D ∂2 u (x ,t )∂x
2
+g (x ,t ),(4
)其中
u (x ,t )=Δt α-1Γ(1+α)ðn
m =1
[u (x ,t m )-u (x ,t m -1)][(n -m +1)α-(n -m )α]+Δt β
-1Γ(1+β)ðn
m =1
[u (x ,t m )-u (x ,t m -1)][(n -m +1)β-(n -m )β].(5
)1.3 格子B o l t z m a n n 模型
下面用格子B o l t z m a n n 方法恢复宏观方程.定义
ðα
f α
(
x ,t )=u (x ,t ),(6
)这里f α(x ,t )表示在时间为t ㊁位置为x ㊁速度为e α的粒子的分布函数,
α为粒子速度方向.格子B o l t z m a n n 方程为
f
α(x +εe α,t +ε)-f α(x ,t )=-1τ
[f
α(x ,t )-f e q
α(x ,t )]+Ωα(x ,t ),(7
)4331 吉林大学学报(理学版) 第61卷
其中f e q
α(x ,t )是f α(x ,t )的平衡态分布函数,τ是松弛因子,Ωα(
x ,t )为源项.将方程(7)左端进行T a y
l o r 展开,可得f α(x +εe α,t +ε)-f α(x ,t )=ðɕ
n =1εn
n
!∂∂t +e α∂∂æèçöø÷x n
f α(
x ,t )+Ωα(x ,t ).(8)引入小K n u d s e n 数ε,定义ε=l /L =Δt ,其中l 为平均自由程,L 为宏观特征长度.在引入ε的情况
下,对分布函数f α(x ,t )进行C h a p m a n -E n s k o g 展开,有f α(x ,t )=f (0)α+εf (1)α+ε2f (2)α+ε3f (3)α+ε4f
(4)α+O (ε5
),(9)其中f e q
α恒为f 0α,于是有
ðα
f
n
α
=0, n ȡ1.(10)进行时间多尺度展开可得t =t 0+εt 1+ε2t 2+ε3t 3+ε4t 4+O (
ε5
),(11)从而有
∂∂t =∂∂t 0+ε∂∂t 1+ε2∂∂t 2+ε3∂∂t 3+ε4∂∂t 4
+O (
ε5
),(12)源项有如下表达形式:
Ωα=
εΩ1+ε2Ω2+ε3Ω3+ε4Ω4+O (ε5
).(13)将式(8),(9),(11)~(13)代入式(7
),并比较ε的各阶系数有ε:Δf
0α=-1τ
f 1α+Ω1
α,(14)ε2:∂∂t 1
f 0α+Δ2f 0αc 2+τΔΩ1α=-1τf 2α+Ω2
α,(15
)其中Δ=∂∂x e α+∂∂t 0
,c 2=
12-τ.根据质量守恒定律有ðα
f α
(
x ,t )=ðα
f e q
α
(x ,t ).
(16)设源项满足如下条件:
ðαΩ
1
α
=0,(17)ð
α
εΩ2
α
=g (x ,t ),(18)平衡态分布函数满足:
ðαf
e q
α
=u (
x ,t ),(19)M 0=ðα
f e q
αe
α=0,(20)π=ðα
f e q αe
αe α=λI u (x ,t ),(21
)其中λ是一个常数,I 是单位矩阵.根据式(5)可知 u (x ,t )是与u (x ,t )相关的量.对式(14)
求和,对式(15)两端乘ε后求和,再将两者相加,并根据式(19)~(21
)可得∂u ∂t +∂M 0∂x +εc 2∂2M 0∂t 0∂x +εc 2∂2
π∂x ∂x +O (ε2)=ðα
(Ω1α+εΩ2α-ετΔΩ1
α)
.(22)根据式(17)~(21)对式(22
)进行整理,即可恢复出宏观方程∂u (x ,t )∂t =D ∂2 u (x ,t )∂x
2
+g (x ,t )+O (ε2
),(23
)图1 D 1Q
3的格子模型F i g .1 L a t t i c em o d e l o fD 1Q
3这里D =-λεc 2.
下面推导平衡态分布函数.D 1Q 3的格子模型
如图1所示.
5
331 第6期 刘 鑫,等:用格子B o l t z m a n n 方法求解修正的时间分数阶方程
图2 D 2Q 9的格子模型F i g .2 L a t t i c em o d e l o fD 2Q
9根据式(19)~(21
)可得平衡态分布函数为f e q
α=λ u 2c
2, α=1,2,(24)f e q
0=u -λ u c
2,
(25
)其中c =e α是粒子沿各方向的移动速度.
D 2Q 9的格子模型如图2所示.
同理可得平衡态分布函数为
f e q
α=λ u 3c 2 (α=1,2,3,4), c =e α,(26)f e q
α=λ u 12c 2 (α=5,6,7,8), 2c =e α,(27)f e q
0=u -5λ u 3c
2.
(28
)2 数值计算
下面计算一个一维实例,因此所涉及的矢量x 这里均为标量x .定义u N (x k ,
t )为数值解,u E
(x k ,
t )为精确解.选择相对误差㊁最大绝对误差和全局相对误差进行验证,各误差表达式分别如下:G (x k ,t )=u N (x k ,t )-u E
(x k ,t )u E (x k ,t ),(29)M =m a x k
u N (x k ,t )-u E (x k ,
t ),(30)E =
ðk
u
N
(x k ,t )-u E
(x k ,
t )ðk
u E
(
x k
,t ).
(31
) 例1 考虑带有如下初边值条件的时间分数阶方程:
∂u (x ,t )∂
t =D ∂1-α∂t 1-α-∂1-β∂t 1-æèçöø÷β∂2u (x ,t )∂x 2
+g (x ,t ),
t ɪ[0,1], x ɪ[0,1], αɪ(0,1), βɪ(0,1),u (x ,0)=0, x ɪ[0,1],u (0,t )=t 2, u (1,t )=t 2e , x ɪ[0,1], t ɪ[0,1]ìî
íïïïïï
ï.(32
)方程(32)的源项为g
(x ,t )=2e x t -t 1+αΓ(2+α)+t 1+β
Γ(2+βæèç
öø
÷),(33)其精确解为
u (x ,t )=
t 2e x
.(34
)数值计算的参数选取:时间步长和空间步长分别为Δt =1650,Δx =165
,格子数为65ˑ65,松弛因子
τ=1.25,分数阶取值为β=0.01,α=0.03,扩散系数为D =0.001.
图3给出了例1的数值计算结果㊁误差分析和收敛阶.由图3(A )可见,数值解和精确解曲线吻合度较高,几乎完全重叠;由图3(B )可见,最大绝对误差为0.00910,最大相对误差为0.00452,两种误差都较小;由图3(C )可见,在一维情形下格子B o l t z m a n n 模型具有空间一阶收敛阶.
下面计算一个二维实例,故x ,y 是标量.同样选择相对误差㊁最大绝对误差和全局相对误差进行
验证,定义数值解为u N (x k ,y i ,t )㊁精确解为u E (x k ,y i ,t ).定义^u N (x k ,t ),^u E (x k ,t )为关于x k 和t 的函数,存在如下数量关系:
6331 吉林大学学报(理学版) 第61卷
图3 例1的数值计算结果(A )㊁误差分析(B )和收敛阶(C )
F i g .3 N u m e r i c a l c a l c u l a t i o n r e s u l t s (A ),e r r o r a n a l y s i s (B )a n d c o n v e r g e n c e o r d e r (C )o f e x a m p
l e 1ði
u
N
(x k ,y i ,t )=^u N (x k ,t ), ði
u E (x k ,y
i ,t )=^u E (x k ,t ).(35)相对误差㊁最大绝对误差和全局相对误差的表达式如下:
G (x k ,t )=^u N (x k ,t )-^u E (x k ,t )^u E (x k ,
t ),(36)M =m a x k ,
i
u N (x k ,y i ,t )-u E (x k ,y i ,t ),(37)E =
ðk ,i
u
N
(x k ,y i ,t )-u E (x k ,y i ,
t )ð
k ,i
u E (x k ,y
i ,t ).
(38
) 例2 考虑带有如下初边值条件的时间分数阶方程:
∂u (x ,y ,t )∂t =D ∂1-α∂t 1-α-∂1-β∂t 1-æèçöø÷β∂2u (x ,y ,t )∂x 2+∂2u (x ,y ,t )∂y
æèçöø÷2+g (x ,y ,t ), t ɪ[0,1], x ,y ɪ[0,1], αɪ(0,1), βɪ(0,1),u (x ,y ,0)=0, x ,y ɪ[0,1],u (0,y ,t )=t 2e y , u (1,y ,t )=t 2e 1
+y , u (x ,0,t )=t 2e x , u (x ,1,t )=
t 2e 1
+x , t ɪ[0,1]ìîíïïïïïïï
ï.(39
)方程(39)的源项为g (x ,y ,t )=2
e x +y t -t 1+αΓ(2+α)+t 1+β
Γ(2+βæèçöø
÷),(40)精确解为
u (x ,y ,
t )=t 2e x +y .(41
)数值计算的参数选取:时间步长和空间步长分别为
Δt =1650, Δx =165, Δy =1
65
,
格子数为65ˑ65,时间t =1.0,松弛因子τ=1.25,分数阶取值为β=0.001,α=0.001,扩散系数为
D =0.001.
图4给出了例2的二维数值计算结果㊁误差分析和收敛阶.由图4(A ),(B )可见,方程(39)的数值解和精确解基本一致;由图4(C )可见,最大相对误差是0.0122,最大绝对误差是0.00906;由图4(D )可见,在二维情形下格子B o l t z m a n n 模型具有空间一阶精度.
综上所述,本文基于格子B o l t z m a n n 方法,针对修正的时间分数阶方程进行了讨论.首先将方程
离散化处理,得到了标准的反应扩散方程.其次,根据T a y l o r 展开式和C h a p m a n -E n s k o g 多尺度展开等技术,用格子B o l t z m a n n 方法准确地恢复了宏观方程,并推导出各方向的平衡态分布函数表达式.
最后,用一个一维实例和一个二维实例进行数值计算,所得数值解与精确解吻合较好,误差在合理范围内.因此,格子B o l t z m a n n 模型在求解修正的时间分数阶方程的数值解上效果较好.
7
331 第6期 刘 鑫,等:用格子B o l t z m a n n 方法求解修正的时间分数阶方程
8331吉林大学学报(理学版)第61卷
图4例2的数值计算结果(A),(B)㊁误差分析(C)和收敛阶(D)
F i g.4N u m e r i c a l c a l c u l a t i o n r e s u l t s(A),(B),e r r o r a n a l y s i s(C)a n d c o n v e r g e n c e o r d e r(D)o f e x a m p l e2
参考文献
[1] P R A K A S HJ,S A T Y A N A R A Y A N AC.A x i s y mm e t r i cS l o w M o t i o no f aP o r o u sS p h e r i c a lP a r t i c l e i naV i s c o u s
F l u i dU s i n g T i m eF r a c t i o n a lN a v i e r-S t o k e sE q u a t i o n[J].C o l l o i d s a n d I n t e r f a c e s,2021,5(2):24-1-24-16.
[2] S U N H G,J I A N GLJ,X I A Y.L B MS i m u l a t i o n o fN o n-N e w t o n i a nF l u i dS e e p a g eB a s e d o nF r a c t i o n a l-D e r i v a t i v e
C o n s t i t u t i v eM o d e l[J].J o u r n a l o fP e t r o l e u mS c i e n c e&E n g i n e e r i n g,2022,213:110378-1-110378-9.
[3]朱帅润,李绍红,钟采伊,等.时间分数阶的非饱和渗流数值分析及其应用[J].应用数学和力学,2022,43(9):
966-975.(Z HUSR,L IS H,Z HO N G C Y,e t a l.N u m e r i c a lA n a l y s i so fT i m eF r a c t i o n a l-O r d e rU n s a t u r a t e d
F l o wa n d I t sA p p l i c a t i o n[J].A p p l i e d M a t h e m a t i c s a n d M e c h a n i c s,2022,43(9):966-975.)
[4]周兆杰,王方圆,王起明,等.分数阶方程约束最优控制问题数值算法研究进展[J].山东师范大学学报(自然科
学版),2022,37(1):1-9.(Z HO UZ J,WA N GFY,WA N GQ M,e t a l.A d v a n c e s i nN u m e r i c a lA l g o r i t h m s f o r O p t i m a lC o n t r o lP r o b l e m s G o v e r n e d b y F r a c t i o n a lD i f f e r e n t i a lE q u a t i o n s[J].J o u r n a lo fS h a n d o n g N o r m a l U n i v e r s i t y(N a t u r a l S c i e n c e),2022,37(1):1-9.)
[5] MU K H T A RS,S HA H R,N O O RS.T h eN u m e r i c a l I n v e s t i g a t i o n o f aF r a c t i o n a l-O r d e rM u l t i-d i m e n s i o n a lM o d e l
o fN a v i e r-S t o k e sE q u a t i o nv i aN o v e lT e c h n i q u e s[J].S y mm e t r y,2022,14(6):1102-1-1102-21.
[6] A L I Z,N I ASN,R A B I E IF,e t a l.AS e m i a n a l y t i c a lA p p r o a c h t o t h eS o l u t i o no fT i m e-F r a c t i o n a lN a v i e r-S t o k e s
E q u a t i o n[J].A d v a n c e s i n M a t h e m a t i c a l P h y s i c s,2021,2021:5547804-1-5547804-13.
[7] B U R Q A N A,E L-A J O U A,S A A D E H R,e ta l.A N e w E f f i c i e n tT e c h n i q u e U s i n g L a p l a c e T r a n s f o r m sa n d
S m o o t hE x p a n s i o n s t oC o n s t r u c t a S e r i e s S o l u t i o n t o t h eT i m e-F r a c t i o n a lN a v i e r-S t o k e sE q u a t i o n s[J].A l e x a n d r i a
E n g i n e e r i n g J o u r n a l,2022,61(2):1069-1077.
[8] C H E N Y,C H E NC M.N u m e r i c a l S c h e m ew i t hH i g hO r d e rA c c u r a c y f o r S o l v i n g aM o d i f i e dF r a c t i o n a l D i f f u s i o n
E q u a t i o n[J].A p p l i e d M a t h e m a t i c s&C o m p u t a t i o n,2014,244:772-782.
[9]吴强,黄建华.分数阶微积分[M].北京:清华大学出版社,2016:74-76.(WU Q,HU A N GJH.F r a c t i o n a l
C a l c u l u s[M].B e i j i n g:T s i n g h u aU n i v e r s i t y P r e s s,2016:74-76.)
(责任编辑:李琦)。