非線性偏微分方程的高階格子BGK模型

1、中國科學 G 輯: 物理學 力學 天文學 2009 年 第 39 卷 第 7 期: 913 922 HYPERLINK / 中國科學雜志社SCIENCE IN CHINA PRESS非線性偏微分方程的高階格子 BGK 模型賴惠林, 馬昌鳳*福建師范大學數學與計算機科學學院, 福州 350007* 聯(lián)系人, 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 多尺度展開 技術, 不同類型的非線性偏微分方程從連續(xù)的 Boltzmann 方程中得到了正確恢 復, 數值模擬結果表明該方法十分有效.關鍵詞 非線性偏微分方程 格子 Boltzmann 模型 多尺度技術 泰勒級數展開Chapman-Enskog 技術展開非線性偏微分方程(nonlinear partial differential equations, NPDEs)在物理學和數學等不同領域里扮 演著非常重要的角色1,2. 在物理系統(tǒng)中, 許多很有

3、 趣且有用的特性隱含在它們的非線性行為里, 如果 它們非線性方程的解析解存在, 則通過這些解析解 我們可以更好地了解其復雜物理現象和動力學過程 的機理. 因此, 尋找和構造 NPDEs 的解析解顯得越 來越重要36. 但是由于解析解一般只存在于某些比 較嚴格的條件下, 大部分 NPDEs 的研究是用近似的 數值方法來處理方程中的非線性項. 近十幾年來, 許 多關于這些 NPDEs 的數值模擬方法發(fā)展起來, 包括 有限差分方法7,8、熱平衡積分法9、有限元法10、譜 方法11、變分迭代法等12.作為一種新興的數值方法, 格子 Boltzmann 方法 (lattice Boltzmann met

4、hod, LBM)不同于傳統(tǒng)的數值 方法, 它是基于微觀模型和細觀運動論的介觀方法. LBM 在求解非線性方程以及復雜系統(tǒng)的演化1316, 特別在流體力學的研究中取得了很大成果1719, 這是由于 LBM 具有物理圖像清晰、邊界處理容易、編 程實現簡單等優(yōu)點. 從計算的角度看, LBM 屬于顯示 時間推進方法，每個時間步的計算量為 O(MN) (M 為 離散速度數, N 為計算格點數), 其計算效率要高于一 般的數值方法. 由于模型所涉及的計算都是具有局 部性, 所需平衡態(tài)分布函數是同時進行計算的, 具有 天然的本質并行性, 非常適合在大規(guī)模并行計算機 上運行. 再加上“遷移”步在實際編程計算

5、中只是一個 賦值過程, 并不占計算時間, 所以其效率是比較高的, 這也是這個方法蓬勃發(fā)展的原因之一. LBM 提供了 聯(lián)系宏觀和微觀的可能性和現實性, 除了在一般的 流體力學問題中得到了成功的驗證之外, 在湍流、多 相流、多組分流、粒子懸浮流、量子力學以及磁流體 力學等相關領域也具有廣闊的應用前景.近年來, 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的啟發(fā)下, 構造出一個帶修正項和一個源項的高 階格子 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態(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) 分別為分布函數和局域平衡態(tài)6 1 t1 4x 分布函數,hi (x, t) 為一修正項函數, c 為常數, 為弛t 4 cef (0) O(t5 )豫時間, 穩(wěn)定性要求 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那么, 通過選擇適當的局域平衡態(tài)分布函數 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), 我們選擇局域平衡態(tài)分布 函數 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)其中為唯一的自由參數, 用來調整局域平衡態(tài)分布i ii其中 1, 2 和為待定參數.對方程(14), (19)和(24)兩邊關于 i 分別求和, 并函數, 以提高模型的精度和穩(wěn)定性. 由方程(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)式聯(lián)立方程組, 可確定局域平衡態(tài)1 0,2(n 1)ct2分布函數的具體表達式如下:其中自由參數為.916中國科學 G 輯: 物理學 力學 天文學2009 年 第 39 卷 第 7 期2數值模擬為驗證上述模型的有效性, 本節(jié)將給出幾種不 同類型的非線性偏微分方程的數值例子, 并與相應 的解析解相比較, 其中包括 KdV-Burgers 方程, 組合 KdV-MKdV 方程, KdV 方程, Burgers-Huxley 方程.i我們通過設定局域平衡態(tài)分布

21、函數 f (0) (x, t) 在 t 0時的值來初始化分布函數 fi (x, t),宏觀量 u(x, t) 可由初始條件來初始化. 對于邊界處理, 我們統(tǒng)一采用 Guo 等人提出的非平衡態(tài)外推格式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. 數值模擬區(qū)域為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 方程數值解與解析 解之間的整體相對誤差擬區(qū)域為 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. 數值模擬區(qū)域為 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 種不同時刻 下數值解與解析解在各個節(jié)點處的比較.算例 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. 數值模擬區(qū) 域為 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 多尺度展開技術, 幾類非線性偏微分方程從 連續(xù)的 Boltzmann 方程中得到還原. 通過幾個與解析 解進行比較的數值實驗, 我們驗證了本文所提模型 的數值有效性. 在前面的 3 個數值實驗中, 我們必須 確保 0.8, 以保證數值結果

33、的穩(wěn)定性, 這由自由參 數來控制; 在第四個算例中, 由于= 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 的值落在區(qū)域(0.5,1.5)內, 而這由來控制. 自由參數, 的選擇使得模型更具有靈活性. 數值模 擬結果表明, 經過長時間的演化, 數值解與解析解仍然十分吻合. 本文的高階格子 BGK 模型可推廣到更 復雜和維數更高的模型, 我們將在今后繼續(xù)展開相 關研究.參考文獻 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