一类三阶变系数偏微分方程的格子Boltzmann模型

武芳芳, 王可心

(沈阳工业大学 理学院, 沈阳 110870)

0 引 言

格子Boltzmann方法(lattice Boltzmann method, LBM)在计算流体力学中应用广泛[1-2], 具有物理背景清晰、 边界易处理、 并行性良好和可扩展等优点[3]. 目前, LBM已成为一种求解各类非线性偏微分方程的有效数值方法[4-8].

考虑具有如下一般形式的一类三阶变系数偏微分方程：

(1)

其中a(t),b(t),c(t),d(t)均为t的任意函数,u是空间坐标x和时间坐标t的波幅函数, 参数p是非负常数,n为正整数.方程(1)在固体材料、 等离子体、 流体等物理学领域应用广泛[9-10].在等离子体和流体中, 方程(1)的特例如下：

1) 当a(t)=b(t)=0,n=1,F(u)=0时, 方程(1)为变系数KdV方程:

(2)

2) 当n=1,p=1,F(u)=0时, 方程(1)为修正的变系数KdV方程:

(3)

3) 当a(t)=0,b(t)=-3β,c(t)=2α,d(t)=1,n=1,p=1,F(u)=0时, 方程(1)为Gardner方程:

(4)

大多数变系数偏微分方程很难找到其解析解, 因此人们致力于寻找适当的数值方法求得这些复杂方程的数值解, 例如有限元法、 有限差分法和谱方法等[11-13].本文给出求解变系数偏微分方程(1)的高精度格子Boltzmann模型, 并通过数值算例验证该模型的有效性和精度.

1 格子Boltzmann模型

对方程(1)采用D1Q5速度模型, 其离散速度集合为{e0,e1,e2,e3,e4}={0,1,-1,2,-2}.格子Boltzmann模型[14]的局部粒子分布函数演化方程为

(5)

将式(5)等号左侧对空间x和时间t进行Taylor展开, 可得

再进行如下形式的多尺度Chapman-Enskog展开[14]：

(7)

其中ε为Knudsen数, 可表示为ε=l/L,l表示平均自由程,L表示特征长度.将式(7)代入式(6), 并比较ε的各阶系数, 可得:

定义宏观量u(x,t)为该点各方向的局部粒子分布函数之和, 即

综上所述，现代电气工程的建设和发展过程中，对于自动化技术及智能化已经取得了明显的进步，通过智能化技术已经提升了电气工程的发展效果，为了保障整体电气工程的技术应用，还应该注重对电气智能化控制技术进行分析，充分发挥电气工程智能化和自动化的优势，将电气智能化更好的应用到电气工程的各个领域，从PLC技术到故障诊断技术全面提升电气智能化的技术效果。

(12)

由质量守恒知, 局部平衡态分布函数应满足

(13)

则由式(7),(8),(13), 有

(14)

(15)

并且选取的补偿函数满足

(16)

由式(7),(16)可得

(17)

其中λ1,λ2,λ3,λ4为待定参数, 满足

a(t)=τΔtcλ1,b(t)=(2p+1)τΔtcλ2,c(t)=(n+1)τΔtcλ3,d(t)=Δt2c3(τ2-τ+1/6)λ4.

(18)

对式(9)关于i求和, 并结合式(14),(15), 有

(19)

对式(10)关于i求和, 并结合式(14),(15),(17),(19), 有

(20)

对式(11)关于i求和, 并结合式(14),(15),(17)～(20), 有

(21)

由式(19)×ε+式(20)×ε2+式(21)×ε3可恢复出具有三阶精度的宏观方程：

合并式(14),(15)可得一组平衡态分布函数:

根据式(16)可选取一组补偿函数为

2 数值算例

为验证LBM的有效性, 分别定义全局误差(GRE)和最大绝对误差(MAE)为

其中u(xi,t),uexact(xi,t)分别表示数值解和精确解.下面对变系数偏微分方程(1)进行数值模拟.宏观的初边值条件由算例方程的精确解确定, 初始的粒子分布函数直接等于平衡态分布函数, 边界处理采用非平衡态外推格式[16].

例1考虑方程(1), 当a(t)=b(t)=0,c(t)=6cos(2t),d(t)=cos(2t),n=1,F(u)=0时方程形式如下：

(22)

其单孤子情形的精确解[10]为

(23)

计算中, 选取参数r=0.5, 无量纲松弛时间τ=0.958, 时间步长Δt=10-5, 空间步长Δx=0.2, 计算域固定为[0,40].图1为例1中t=10和t=20两个时刻数值解与精确解的对比, 表1列出了例1不同时刻的GRE和MAE. 结果表明, LBM可以得到精度较高的数值解. 图2为例1在对数坐标系下不同时刻的GRE与格子数NX之间的关系.由图2可见, 直线斜率接近3.0, 表明本文提出的模型具有空间三阶精度, 与理论精度相符.

图1 不同时刻例1数值解和精确解的对比Fig.1 Comparison of numerical solutions and exact solutions of example 1 at different time

图2 例1不同时刻的精度分析Fig.2 Accuracy analysis of example 1 at different time

表1 例1不同时刻的误差

例2考虑方程(1), 当a(t)=b(t)=0,c(t)=6t2,d(t)=t2,n=1,F(u)=0时方程形式如下：

(24)

其双孤子情形的精确解为

(25)

计算中, 选取参数k1=0.5,k2=0.7, 无量纲松弛时间τ=0.958, 时间步长为Δt=10-6, 空间步长为Δx=0.25, 计算域固定为[-20,35].图3为例2不同时刻数值解与精确解的对比, 图4为例2数值解和精确解随时间的演化过程.由图3和图4可见, 数值解与精确解吻合较好.表2列出了例2不同时刻的GRE和MAE. 由表2可见, LBM可以得到精度较高的数值解.

图3 不同时刻例2数值解和精确解的对比Fig.3 Comparison of numerical solutions and exact solutions of example 2 at different time

图4 例2数值解和精确解随时间的演化过程Fig.4 Evolution process of numerical solutions and exact solutions with time of example 2

表2 例2不同时刻的误差

例3考虑方程(1), 当a(t)=0,b(t)=-3β,c(t)=2α,d(t)=1,n=1,p=1,F(u)=0时, 方程形式如下：

(26)

其单孤子解情形的精确解为

(27)

计算中, 选取参数β=0.5,α=5,γ=2, 无量纲松弛时间τ=0.957, 时间步长为Δt=10-5, 空间步长为Δx=0.2, 计算域固定为[-10,20].图5为例3不同时刻数值解与精确解的对比.由图5可见, 数值解与精确解吻合较好.表3列出了本文例3的LBM与文献[8]的LBM在不同时刻的误差对比. 由表3可见, 本文的LBM数值解更精确、 演化时间更长, 进一步验证了本文建立模型的有效性.

图5 例3不同时刻数值解和精确解的对比Fig.5 Comparison of numerical solutions and exact solutions of example 3 at different time

综上所述, 针对一类含有源项的变系数三阶微分方程, 本文通过适当地选取平衡态分布函数和补偿函数, 合理地构造了具有三阶精度的格子Boltzmann模型, 其中补偿函数可同时恢复出宏观方程的对流项和源项. 数值算例模拟结果验证了本文模型的有效性和精度. 该格子Boltzmann模型具有一般性, 还可用于求解三阶常系数偏微分方程.