非线性偏微分方程的高阶格子BGK模型

1、中国科学 G 辑: 物理学 力学 天文学 2009 年 第 39 卷 第 7 期: 913 922 HYPERLINK / 中国科学杂志社SCIENCE IN CHINA PRESS非线性偏微分方程的高阶格子 BGK 模型赖惠林, 马昌凤*福建师范大学数学与计算机科学学院, 福州 350007* 联系人, E-mail: HYPERLINK mailto:macf macf收稿日期: 2008-09-25; 接受日期: 2009-03-30国家自然科学基金(批准号: 10661005)和福建省科技厅 K 类基金(编号: 2008F5019)资助项目摘要 考虑一维含源非线性偏微分方程: ut +

2、 uux + unux uxx + uxxx = F(u), 建立 D1Q5 带修正项的高阶格子 BGK 模型, 通过 Chapman-Enskog 多尺度展开 技术, 不同类型的非线性偏微分方程从连续的 Boltzmann 方程中得到了正确恢 复, 数值模拟结果表明该方法十分有效.关键词 非线性偏微分方程 格子 Boltzmann 模型 多尺度技术 泰勒级数展开Chapman-Enskog 技术展开非线性偏微分方程(nonlinear partial differential equations, NPDEs)在物理学和数学等不同领域里扮 演着非常重要的角色1,2. 在物理系统中, 许多很有

3、 趣且有用的特性隐含在它们的非线性行为里, 如果 它们非线性方程的解析解存在, 则通过这些解析解 我们可以更好地了解其复杂物理现象和动力学过程 的机理. 因此, 寻找和构造 NPDEs 的解析解显得越 来越重要36. 但是由于解析解一般只存在于某些比 较严格的条件下, 大部分 NPDEs 的研究是用近似的 数值方法来处理方程中的非线性项. 近十几年来, 许 多关于这些 NPDEs 的数值模拟方法发展起来, 包括 有限差分方法7,8、热平衡积分法9、有限元法10、谱 方法11、变分迭代法等12.作为一种新兴的数值方法, 格子 Boltzmann 方法 (lattice Boltzmann met

4、hod, LBM)不同于传统的数值 方法, 它是基于微观模型和细观运动论的介观方法. LBM 在求解非线性方程以及复杂系统的演化1316, 特别在流体力学的研究中取得了很大成果1719, 这是由于 LBM 具有物理图像清晰、边界处理容易、编 程实现简单等优点. 从计算的角度看, LBM 属于显示 时间推进方法，每个时间步的计算量为 O(MN) (M 为 离散速度数, N 为计算格点数), 其计算效率要高于一 般的数值方法. 由于模型所涉及的计算都是具有局 部性, 所需平衡态分布函数是同时进行计算的, 具有 天然的本质并行性, 非常适合在大规模并行计算机 上运行. 再加上“迁移”步在实际编程计算

5、中只是一个 赋值过程, 并不占计算时间, 所以其效率是比较高的, 这也是这个方法蓬勃发展的原因之一. LBM 提供了 联系宏观和微观的可能性和现实性, 除了在一般的 流体力学问题中得到了成功的验证之外, 在湍流、多 相流、多组分流、粒子悬浮流、量子力学以及磁流体 力学等相关领域也具有广阔的应用前景.近年来, LBM 成功应用于模拟某些复杂的非线 性演化方程, 如对流扩散方程20、反应扩散方程21、 Burgers 方程22、MKdV 方程23、KdV-Burgers 方程24. 但是, 在现有的许多格子 Boltzmann 模型中存在一些引用格式: 赖惠林, 马昌凤. 非线性偏微分方程的高阶格

6、子 BGK 模型. 中国科学 G 辑, 2009, 39(7): 913922Lai H L, Ma C F. A higher order lattice BGK model for simulating some nonlinear partial differential equations. Sci China Ser G, 2009, 52(7): 1053 1061, doi: 10.1007/s11433-009-0149-3赖惠林等: 非线性偏微分方程的高阶格子 BGK 模型问题 , 即如何构 造更高阶 精度模型和 如何导出 O(t5 ) 1 ff (0) t 2h t F (

7、u).(4)NPDEs 中更复杂的非线性项. 因此, 本文在文献25 27的启发下, 构造出一个带修正项和一个源项的高 阶格子 BGK(Bhatnager-Gross-Krook)模型, 提高了模iii5使用 Chapman-Enskog 多尺度展开技术,iiiiif f (0) f (1) 2 f (2) 3 f (3)型的精度, 使得模型具有五阶精度. 通过这个模型,i4 f (4)O(5 ),可以求解更一般的非线性偏微分方程, 包括广义 (5) 234KdV 方程, KdV-Burgers 方程, 组合 KdV-MKdV 方程,t t ttO(),Boussinesq 方程, 广义 Bu

8、rgers-Huxley 方程等. 数值 实验模拟结果表明, 该方法是十分有效的. 为以后更F (u) 12312 F (u),复杂和维数更高的非线性偏微分方程的数值模拟积其中是 Knudsen 数, 定义为 A / L,A 是平均自由累了经验.程, L 是特征长度, 可取为时间步长t, 同时非平衡i态分布函数 f (k ) (k 1, 2,) 满足以下守恒律1格子 BGK 模型考虑一维含源非线性偏微分方程如下:f (k ) 0 (k 1, 2, ).ii(6)u uuunuuuF (u),(1)令 t,并将方程(5)代入方程(4)中, 得txxxxxxx其中, , 和为实常数.我们采用 D1

9、Q5 模型, 使用的离散速度方向定义t tt1t 2 t2t3 t3cei x 为 (i 0, 1, 2, 3, 4)f (0) tf (1) t 2 f (2) t3 f (3) iiii e0 , e1, e2 , e3 , e4 0, 1, 1, 2, 2 ,21 22 带修正项和源项的格子 Boltzmann 方程为t2t(0) tf (1) iiit1tt2ceix fi (x cei t, t t) fi (x, t) 1 f (x, t) f (0) (x, t) t 2 h (x, t) t F (u),(2)ft 2 f (2) 3iii51 t3 t ceif (0) tf

10、 (1) ii其中 fi (x, t),(0)fi(x, t) 分别为分布函数和局域平衡态6 1 t1 4x 分布函数,hi (x, t) 为一修正项函数, c 为常数, 为弛t 4 cef (0) O(t5 )豫时间, 稳定性要求 0.5.241i x i (1) 2 (2) 3 (3) 4 (4) 宏观变量 u 满足如下守恒律:tfit fit fit fiu(x, t) fi (x, t) f (0) (x, t).(3)t 2h t3F (u).(7)iii那么, 通过选择适当的局域平衡态分布函数 f (0) (x, t),i51i比较方程(7)两端小参数t 的一阶项得 O(t):利用

11、 Chapman-Enskog 多尺度展开技术, 我们可以将c e f (0)1 f (1) ,i ii(8)方程(2)还原成宏观方程(1).事实上, 对方程(2)左边进行泰勒展开, 并保留由方程(8)我们有x至 O(t5 ) 项, 我们得到f (1)ce f (0) ,12 ixi ii i(9)2t t cei x fi 2 tt cei x fie f (1) ce2 f (0) ,i i13 3 1 4 4e2 f (1) c t cei6 txi 24 t t cei xxxe3 f (0) ,(10)(11)ffii ii i914中国科学 G 辑: 物理学 力学 天文学2009

12、年 第 39 卷 第 7 期e3 f (1)c e4 f (0) .e f (3) e Fe f (0)i ixi i(12)i i( i 1 ) 5t2i i比较方程(7)两端小参数t 的二阶项, 得 O(t2):c(21) e2 f (0) c2(e2h )i ii i f (0)c e f (1) 1 c22 2 (0) t1 xxit1xi ie f2x2i i c3 2 1 3 4 (0) 1 f (2) h .(13)3i ie f6.x(21)ii再由方程(14)得22把方程(10)代入(13)式得 f (0) c2 1 1 e2 f (0) 22 i2i i f (0) c2

13、1 1 e2 f (0) t12t1 xti2x2i i1 (2)1fh .(22) 11 f (2) h .(14)t1itiii比较(7)式两端小参数t 的四阶项得 O(t4):由方程(14)我们有 f (0) f (2) f (1) c e f (3) 2tititixi if (2) h f (0) c2 2 1 1 e2 f (0) ,(15)312iiti2x2i i1 2 1f (0) ce f (0) 12e f (2) (e h ) e f (0) 2 t 2it xi i1e fi ii ii iti ic e f (1) 1 c22 2 (2)2 c2 2 1 1 e3

14、f (0) ,(16)t1 x2x2i i13 3 3 (1) e2 f (2) (e2h ) e2 f (0) c6x3ei fii ii ii i24t1 1 c2 e2 f (0) 1 c4 e4 f (0) i ii i2 c2 2 1 1 e4 f (0) .(17)2t1 x224x42x2i i1 f (4) .(23) i比较方程(7)两端小参数 t 的三阶项得 O(t3):2把方程(12), (17), (21)和(22)代入(23)式, 得 f (0) f (1) c e f (2) 1 c2 e2 f (1) 1 1 titixi i2x2i if (0) 1f (2)

15、2 f (1)21ti2titi3312c e f (0) 1 c3 e3 f (0) c1 2t xi i6x3i i(e F ) c2 22 2e2 f (0) 11 f (3) F1 (u) .(18)5 x i 124 t1x2i ii5c21 (e2 h )2x2i i把方程(11)和(16)代入(18)式, 得 f (0) 2 1 f (1) c(e h )43c 3 227 1 1224 iii it2t1x34e4 f (0)1 h 1 f (4) .c3 2 1 e3 f (0) 1 f (3) 1 F (u).(19)x4i i2 tii(24)6 x3i ii5 11由方

16、程(19)我们有为恢复宏观方程(1), 我们选择局域平衡态分布 函数 f (0) (i 0, 1, 2, 3, 4) 满足:(3) (0) (1)2 ifiF1 fi5t(12)fitc(eihi )xe f (0) 0,e2 f (0) 0,21i ii i3 c3 2 1 e3 f (0) ,iie3 f (0) u,e4 f (0) 0,(25)6 x3i i(20)i ii iii915赖惠林等: 非线性偏微分方程的高阶格子 BGK 模型其中为待定参数.f (0) u,f (0) u,f (0) u,同时修正项函数 hi (x, t) (i 0,1, 2, 3, 4) 满足:01626

17、(33)2n1f (0) u,f (0) u,hi 0,eihi 1u2u,312412iie2 h u,(26)其中为唯一的自由参数, 用来调整局域平衡态分布i ii其中 1, 2 和为待定参数.对方程(14), (19)和(24)两边关于 i 分别求和, 并函数, 以提高模型的精度和稳定性. 由方程(26)可得修正函数 hi (x, t) 的表达式, 为了简单起见, 我们只给出其中一种情况:利用方程(25)和(26), 得h06u2 u,O(t2):1121n1h1 2u 2 1 4 u2u,2u 0,(27) hut111221 4 u2 1 un1,(34)O(t3):2222ut22

18、cuu (n 1)cunuh3 u ,2h4 u .1x2x6 c3 2 1 uF (u),(28)注: 当 , , , 和 F (u) 取不同的值时, 我们 O(t4):xxx1可以得到不同类型的非线性偏微分方程:() 当 F (u) 0 时, 方程(1)退化为 KdV-Bur-xxu c2 1 u0.(29)gers 方程.t32由(27)式t+(28)式t2+(29)式t3 得u 2ct 2uu (n 1)ct2unu t1x2x() 当 n 2, F (u) 0 时方程(1)退化为组合KdV-MKdV 方程, 若 0 , 则退化为广义 MKdV方程.2 132c t uxx c3 2

19、1 t 2uF (u).(30)( ) 当 F (u) 0 时方程(1) 退化为广义KdV 方程.6 xxx( ) 当 0, 1,并且 F (u) (1 un ) 为恢复宏观方程(1), 我们只需令22(un ) 时方程(1)退化为广义 Burgers-Huxley 方程,2ct 1,(n 1)ct 2 ,如果 F (u) u(1un ),则为广义 Burgers-Fisher 方程.c21 t3, c3 2 1 t 2,(31)26 由于此时 0 , 我们需要选择其他方法来求弛豫时这样我们就有间 , 我们在(31)式中令 0,则松弛时间由式子2 131 1 ,c t 确定, 这样计算所需参数

20、确定 2212c3t 2c2 1 t32 如下: (32)1 1 ,12ct 22.(n 1)ct2416c2t3(35)由(3)式和(25)式联立方程组, 可确定局域平衡态1 0,2(n 1)ct2分布函数的具体表达式如下:其中自由参数为.916中国科学 G 辑: 物理学 力学 天文学2009 年 第 39 卷 第 7 期2数值模拟为验证上述模型的有效性, 本节将给出几种不 同类型的非线性偏微分方程的数值例子, 并与相应 的解析解相比较, 其中包括 KdV-Burgers 方程, 组合 KdV-MKdV 方程, KdV 方程, Burgers-Huxley 方程.i我们通过设定局域平衡态分布

21、函数 f (0) (x, t) 在 t 0时的值来初始化分布函数 fi (x, t),宏观量 u(x, t) 可由初始条件来初始化. 对于边界处理, 我们统一采用 Guo 等人提出的非平衡态外推格式28. 除此之外, 为 测定模型误差精度, 我们定义总体相对误差(global relative error, GRE)为ii| u(x , t) u* (x , t) |图 1 算例 1, KdV-Burgers 方程不同时刻的模拟结果实线代表解析解GRE i,| u* (x , t) |(36)i表 1 算例 1 不同时刻 KdV-Burgers 方程数值解与解析解i*之间的整体相对误差 其中

22、u(xi , t),u (xi , t) 分别为数值解和解析解, 在所tGRE有格点进行求和.算例 1给出 KdV-Burgers 方程24:ut uux uxx uxxx 0,其解析解为u(x, t) 22,1e2( xt ) 2101.0213105501.83971051501.70181052501.45331053001.3539105其中,1062.25在模拟中, 我们取= 1, = 0, = 0.0009, = 0.00002, x = t = 0.01, = 1.473. 数值模拟区域为I 4, 4.模拟结果见图 1 和表 1.算例 2给出组合 KdV-MKdV 方程5:u u

23、uu2uu0,其解析解为txxxxxu(x, t) 6c2 tanh c2 , 0, c0,22 2图 2 算例 2, 组合 KdV-MKdV 方程不同时刻的模拟结果其中 c2 为常数,2x (2 4c )t/4.实线代表解析解在模拟中, 我们取 n = 2, = 10, = 60, = 0, = 1.0, c2 = 0.005, x = 0.1, t = 0.01, = 0.226. 数值模表 2 算例 2 不同时刻组合 KdV-MKdV 方程数值解与解析 解之间的整体相对误差拟区域为 I 200, 200.模拟结果见图 2 和表 2. t GRE 101.0215105算例 3给出 KdV

24、 方程9ut uux uxxx 0,503.8235105 906.761510 5917赖惠林等: 非线性偏微分方程的高阶格子 BGK 模型边界条件为初始条件为u(0, t) u(2, t) 0,t 0,u(x, 0) 3Csech2 ( Ax E), 0 x 2,其解析解为u(x, t) 3Csech2 ( Ax Bt E), 0 x 2,其中 A 12C ,B AC,C 和 E 为常数.在 模拟中 , 我们 取 = 10, = 0, = 0, = 0.000484, C = 0.3, E = 6.0. 数值模拟区域为 I = 0, 2.我们分别给出 3 个不同时刻的数值解与解析解 模拟结

25、果的比较:() 模拟 t = 0.00001 时刻, 我们取= 2.8, x = 0.001, t = 0.000001. 模拟结果见图 3.图 4 算例 3 (), t = 0.005 时刻 KdV 方程的模拟结果图 5 算例 3 (), t = 0.01 时刻 KdV 方程的模拟结果图 3 算例 3 (), t = 0.00001 时刻 KdV 方程的模拟结果1/ n() 模拟 t = 0.005 时刻, 我们取= 0.14, x =u(0, t) tanh A1 A2t22, t 0,1/ n0.001, t = 0.00001. 模拟结果见图 4.u(1, t) tanh A1 (1

26、A2t)22, t 0,() 模拟 t = 0.01 时刻, 我们取= 0.715, x =0.001, t = 0.0001. 模拟结果见图 5.初始条件为1/ n同时, 我们还给出表 3, 分别表示 3 种不同时刻 下数值解与解析解在各个节点处的比较.算例 4给出 Burgers-Huxley 方程29,30:解析解为u(x, 0) tanh A1x22,1/ nu unu uu(1un )(un ),u(x, t) 2 2 tanh A1 (x A2t),txxx0 x 1, t 0.边界条件为n 0, 0, (0,1), x 0,1,其中918中国科学 G 辑: 物理学 力学 天文学2

27、009 年 第 39 卷 第 7 期表 3 算例 3, 不同时刻 KdV 方程解析解与数值解对比t = 0.00001t = 0.005t = 0.01解析解数值解解析解数值解解析解数值解0.00.00.00.00.00.00.00.10.000266630.000266630.000256880.000256880.000247460.000247540.20.003209530.003209530.003092320.003082670.002979150.002978050.30.037944960.037944960.036585290.036554080.035270770.0355

28、18260.40.366156720.366156680.355746110.356638360.345516360.348507690.50.856279070.856279210.863056010.865335190.869319440.874636560.60.172022340.172022370.177874230.177329780.183912190.181397090.70.015681130.015681140.016271190.016301470.016884490.016796650.80.001311080.001311080.001360830.001365050

29、.001412570.001409880.90.000108810.000108810.000112940.000113310.000117240.000117061.00.000009030.000009030.000009370.000009400.000009720.000009711.10.000000750.000000750.000000780.000000780.000000810.000000811.20.000000060.000000060.000000060.000000060.000000070.000000071.30.000000010.000000010.0000

30、00010.000000010.000000010.000000011.40.00.00.00.00.00.02.00.00.00.00.00.00.0 xiA1 (1 n )( 2 4 (1 n)A 21n n n 2 4 (1 n),4(1n).2(1n)表 5 算例 4, 时间相对较大时 Burgers-Huxley 方程数值 解与解析解之间的整体相对误差 tGRE1001.804410510001.09721057在模拟中, 我们取 n = 2, = 0, = 0.1, = 1.0, = 0, = 0, = 0.001, = 0.0001, x = 0.001. 数值模拟区 域为 I

31、= 0,1.我们分别模拟了不同时间尺度的情况, 并与解 析解进行对比, 同时给出总体相对误差:() 当时间相对较小时, 我们取t = 0.0001, = 1.0105, 此时分别模拟了 t = 0.01, t = 0.2, t = 0.5, t =0.8 4 种情况, 模拟结果见图 6(a)(d)和表 4.表 4 算例 4, 时间相对较小时 Burgers-Huxley 方程数值 解与解析解之间的整体相对误差tGRE0.011.34651060.27.00191060.59.02331060.89.3523106() 当时间相对较大时, 我们取t = 0.01, = 2.2, 此时分别模拟了

32、t = 100, t = 1000, t = 10000, t =100000 四种情况, 模拟结果见图 6(e)(h)和表 5. 结100005.9951101000000果表明, 经过长时间的演化, 数值解与解析解仍然十 分吻合.3结论本文针对一类含源非线性偏微分方程, 构造了 具有五阶精度的高阶格子 BGK 模型. 通过 Chapman- Enskog 多尺度展开技术, 几类非线性偏微分方程从 连续的 Boltzmann 方程中得到还原. 通过几个与解析 解进行比较的数值实验, 我们验证了本文所提模型 的数值有效性. 在前面的 3 个数值实验中, 我们必须 确保 0.8, 以保证数值结果

33、的稳定性, 这由自由参 数来控制; 在第四个算例中, 由于= 0, 我们选择另外一种方法来确定弛豫时间, 即由 c2(1/ 2 )t3来确定, 此时的自由参数为, 基于实际的数 值模拟经验, 好的数值模拟结果需要尽量使方程(35)919赖惠林等: 非线性偏微分方程的高阶格子 BGK 模型图 6 算例 4 不同时间尺度 Burgers-Huxley 方程的模拟结果(a) t=0.01; (b) t=0.2; (c) t=0.5; (d) t=0.8; (e) t=100; (f ) t=1000; (g) t=10000; (h) t=100000, 其中实线为解析解920中国科学 G 辑: 物

34、理学 力学 天文学2009 年 第 39 卷 第 7 期中/c2t3 的值落在区域(0.5,1.5)内, 而这由来控制. 自由参数, 的选择使得模型更具有灵活性. 数值模 拟结果表明, 经过长时间的演化, 数值解与解析解仍然十分吻合. 本文的高阶格子 BGK 模型可推广到更 复杂和维数更高的模型, 我们将在今后继续展开相 关研究.参考文献 Dodd R K, Eilbeck J C, Gibbon J D, et al. Solitons and Nonlinear Wave Equations. London: Academic Press, 1982Ablowitz M J, Clarks

35、on P A. Solitons, Nonlinear Evolution Equations and Inverse Scattering. London: Cambridge University Press, 1991Fan E G. Extended tanh-function method and its applications to nonlinear equations. Phys Lett A, 2000, 277: 212219Helal M A. Solition solution of some nonlinear partial differential equati

36、ons and its applications in fluid mechanics. Chaos Solitons Fractals, 2002, 13: 19171929Fan E G. Uniformly constructing a series of explicit exact solutions to nonlinear equations in mathematical physics. Chaos Solitions Fractals, 2003, 16: 819839Yomba E. The extended Fans sub-equation method and it

37、s application to KdV-MKdV, BKK and variant Boussinesq equations. Phys Lett A, 2005, 336: 463476zer S, Kutluay S. An analytical-numerical method for solving the Koreweg-de Vries equation. Appl Math Comput, 2005, 164: 789 797Helal M A, Mehanna M S. A comparison between two different methods for solvin

38、g KdV-Burgers equation. Chaos Solitons Fractals, 2006, 28: 320326Kutluay S, Bahadir A R, zdes A. A small time solutions for the Korteweg-de Vries equation. Appl Math Comput, 2000, 107: 203 210Geyikli T, Kaya D. An application for a modified KdV equation by the decomposition method and finite element

39、 method. Appl Math Comput, 2005, 169: 971981Helal M A. A Chebyshev spectral method for solving Korteweg-de Vries equation with hydrodynamical application. Chaos Solitons Fractals, 2001, 12: 943950Soliman A A. A numerical simulation and explicit solutions of KdV-Burgers and Laxs seventh-order KdV equ

40、ations. Chaos Solitons Fractals, 2006, 29: 294302Doolen G D. Lattice gas and lattice Boltzmann for partial differential equation. Physica D, 1991, 47: l200Benzi R, Succi S, Vergasola M. The lattice Boltzmann equation: Theory and application. Phys Rep, 1992, 222(3): 145197Li Q, He Y L, Gao Y J. Coupl

41、ed double-distribution-function lattice Boltzmann model with an adjustable bulk viscosity. Int J Mod Phys C, 2008, 19: 19191938Lebowitz, Orsazag S A, Qlan Y H. Program of the lattice gas94 meeting at Princeton. J Stat Phys, 1995, 68: 1536Chen S, Doolen G D. Lattice Boltzmann method for fluid flows.

42、Annu Rev Fluid Mech, 1998, 30: 329364Guo Z L, Shi B C, Wang N C. Lattice BGK model for incompressible Navier-Stokes equation. J Comp Phys, 2000, 165: 288306Li Q, He Y L, Wang Y, et al. Coupled double-distribution-function lattice Boltzmann method for the compressible. Phys Rev E, 2007, 76: 056705Deng B, Shi B C, Wang G C. A new lattice Bhatnagar-Gross-Krook model for convection-diffusion equation with a source term. Chin