熱波方程的格子Boltzmann模型

史秀波,閆廣武

(1. 桂林理工大學(xué) 理學(xué)院,廣西 桂林 541004；2. 吉林大學(xué) 數(shù)學(xué)學(xué)院,長春 130012)

格子Boltzmann方法(LBM)作為一種新的數(shù)值方法在計算流體力學(xué)、 非線性偏微分方程等領(lǐng)域受到廣泛關(guān)注[1-2]. 閆廣武等[3-4]將該方法應(yīng)用于波傳播問題,為研究其他波動問題提供了可選擇的途徑. 波動方程中的波速通常是一個常數(shù)值,通過引入動量變量ρuj(x,t)將波動方程轉(zhuǎn)換為小擾動Euler方程進行求解[5-6]. 本文用格子Boltzmann方法對熱波動方程進行模擬. 在該方程中,波速不再是一個常量,而是一個變量,其表達式為

(1)

其中Cs(x)表示波速,是關(guān)于x或y的函數(shù).u(x,t)的下一個時間步表達式為

(2)

本文提出熱波的格子Boltzmann模型,通過使用Chapman-Enskog展開和多尺度技術(shù),得到了系列格子Boltzmann偏微分方程、 平衡態(tài)分布函數(shù)的高階矩及二階精度宏觀熱波方程. 數(shù)值實驗將模型結(jié)果與變分迭代法獲得的解析解及經(jīng)典中心差分格式獲得的結(jié)果進行比較,結(jié)果表明,該方法所得結(jié)果與經(jīng)典方法所得結(jié)果相符.

1 格子Boltzmann模型

選擇一維3-bit網(wǎng)格和二維5-bit網(wǎng)格,分布函數(shù)fα(x,t)定義為在某節(jié)點x上、t時刻、 具有速度eα(α=0,1,…,b)的粒子出現(xiàn)的概率,其中α=0表示靜止粒子. 在一維空間中,b=2,粒子速度為eα={0,c,-c};二維空間中,b=4,粒子速度為

(3)

其中c表示速率. 定義宏觀量:

(4)

(5)

格子Boltzmann方程表示為

fα(x+eα,t+1)-fα(x,t)=Ωα+ωα,

(6)

選取Knudsen數(shù)ε作為數(shù)值模擬的時間步長和Chapman-Enskog展開的小參數(shù)[7],在該尺度上,方程(6)可寫為

fα(x+εeα,t+ε)-fα(x,t)=Ωα+ωα,

(7)

在方程(7)中,假設(shè)

ωα(x,t)=ε2φα(x,t).

(8)

運用Chapman-Enskog展開和多尺度技術(shù),并對方程(7)進行Taylor展開,保留余項到O(ε3)的精度,可得不同時間尺度上的系列格子Boltzmann方程[8]:

結(jié)合方程(5),(9),并假設(shè)

(12)

可得時間尺度t0上的守恒方程

(13)

式(9)+式(10)×ε并對α求和,同時假設(shè)

(14)

得到宏觀熱波方程為

(15)

結(jié)合方程(5),(12),(14)易得平衡態(tài)分布函數(shù)為

其中D表示空間維數(shù)(一維空間中D=1；二維空間中D=2).

2 數(shù)值模擬

為了驗證模型效果,分別對一維和二維熱波問題進行數(shù)值模擬. 一維熱波問題使用3-bit模型,二維熱波問題使用5-bit模型.

例1一維熱波方程

Dirichlet邊界條件：

u(0,t)=0,u(1,t)=1+sinht;

(18)b

初始條件：

u(x,0)=x,ut(x,0)=x2.

(18)c

變分迭代法獲得的精確解[9]為

u=x+x2sinht.

(18)d

選取參數(shù)：格子尺寸m=100,Δx=0.01,c=5.0,τ=1.2,t=1. 圖1(A)為t=1時LBM數(shù)值解和精確解的比較結(jié)果；圖1(B)為兩種結(jié)果的相對誤差Er=|(u-u*)/u*|,其中:u表示LBM數(shù)值解;u*表示精確解. 由圖1(B)可見,相對誤差在(1×10-3,9×10-3)內(nèi),數(shù)值解和精確解吻合較好.

圖1 一維熱波方程LBM數(shù)值解和精確解的比較(A)及相對誤差曲線(B)Fig.1 Comparison of LBM solution and the exact solution of one-dimensional thermal wave equation (A) and the curves of their relative error (B)

例2二維熱波方程

(19)a

Neumann邊界條件：

ux(0,y,t)=0,ux(1,y,t)=2sinht,

uy(x,0,t)=0,uy(x,1,t)=2cosht;

(19)b

初始條件：

u(x,y,0)=y2,ut(x,y,0)=0.

(19)c

選取參數(shù)：格子尺寸m×n=100×100,Δx=0.01,Δy=Δx,c=5,τ=1.01,ε=Δt=Δx/c,t=1. 圖2(A)為t=1時LBM的模擬結(jié)果；圖2(B)為t=1時經(jīng)典中心差分格式的數(shù)值解,將其作為精確解；圖2(C)為兩種結(jié)果在x=0.4處的相對誤差曲線. 由圖2(C)可見,誤差區(qū)域在(0.00,0.05)內(nèi),數(shù)值解與精確解吻合較好.

圖2 二維熱波方程LBM模擬結(jié)果(A)、 精確解(B)和兩種結(jié)果的相對誤差曲線(C)Fig.2 LBM result (A),the exact solution (B) and the curves of their relative error (C) for of two-dimensional thermal wave equation

綜上,本文提出了一個用于熱波方程的格子Boltzmann模型,可得如下結(jié)論：

1) 不同時間尺度的系列偏微分方程對構(gòu)建熱波格子Boltzmann模型非常重要,通過使用高階矩得到了平衡態(tài)分布函數(shù)的表達式；

[1] CHEN Shi-yi,Doolen G D. Lattice Boltzmann Method for Fluid Flow [J]. Annual Review of Fluid Mechanics,1998,30: 329-364.

[2] DUAN Ya-li,KONG Ling-hua,ZHANG Rui. A Lattice Boltzmann Model for the Generalized Burgers-Huxley Equation [J]. Physica A,2012,391(3): 625-632.

[3] YAN Guang-wu. A Lattice Boltzmann Equation for Waves [J]. Journal of Computational Physics,2000,161(1): 61-69.

[4] YAN Guang-wu,DONG Yin-feng. Application of the Lattice Bhatnagar-Gross-Krook Model to the Simulation of Seismic Pressure Wave [J]. Acta Mechanica Sinica,2005,37(2): 238-243. (閆廣武,董銀峰. 基于格子Bhatnagar-Gross-Krook模型的地震壓力波模擬 [J]. 力學(xué)學(xué)報,2005,37(2): 238-243.)

[5] ZHANG Jian-ying,YAN Guang-wu,SHI Xiu-bo. Lattice Boltzmann Model for Wave Propagation [J]. Physical Review E,2009,80(2): 026706.

[6] SHI Xiu-bo,YAN Guang-wu,ZHANG Jian-ying. A Multi-energy-level Lattice Boltzmann Model for Two-Dimensional Wave Equation [J]. International Journal for Numerical Methods in Fluids,2010,64(2): 148-162.

[7] Chapman S,Cowling T G. The Mathematical Theory of Non-uniform Gas [M]. Cambridge: Cambridge University Press,1939.

[8] ZHANG Jian-ying,YAN Guang-wu. Lattice Boltzmann Method for One and Two-Dimensional Burgers Equation [J]. Physica A,2008,387(19/20): 4771-4786.

[9] Wazwaz A M,Gorguis A. Exact Solutions for Heat-Like and Wave-Like Equations with Variable Coefficients [J]. Applied Mathematics and Computation,2004,149(1): 15-29.