KdV方程的格子Boltzmann模型求解

陳夢(mèng)涵,王希胤,李金

(華北理工大學(xué) 理學(xué)院,河北 唐山 063210)

波動(dòng)現(xiàn)象是自然界最普遍的現(xiàn)象之一,對(duì)于水波的研究也一直是科學(xué)和工程研究領(lǐng)域的重要課題。盡管波動(dòng)問題所涉及的領(lǐng)域不同,但是描述波動(dòng)現(xiàn)象的方程卻是相同的。由于水波千姿百態(tài),用肉眼就可以觀察到,因此很早就引起了人們的注意,可以說是人們最為熟悉的一種波。波動(dòng)是物質(zhì)運(yùn)動(dòng)的重要形式,廣泛存在于自然界。波動(dòng)中被傳遞的物理量的擾動(dòng)和振動(dòng)有多種形式,例如,弦線中的波、空氣或固體中的聲波、水波、電磁波,等等。為了更加具體地研究各種波動(dòng),就產(chǎn)生了各種形式的波動(dòng)方程,因此,淺水波方程也成為重要的研究對(duì)象之一。而Korteweg-de Vries(KdV)方程是1895年由荷蘭數(shù)學(xué)家科特韋格(Korteweg)和德弗里斯(de Vries)在研究淺水中小振幅長波運(yùn)動(dòng)時(shí)共同發(fā)現(xiàn)的一種單向運(yùn)動(dòng)淺水波偏微分方程。

在求解偏微分方程地過程中,我們經(jīng)常用到的數(shù)值計(jì)算方法有:有限元法(Finite Element Methods),有限差分法(Finite Diference Methods),有限體積法(Finite Volume Methods)和格子Bolzmann(LBM)方法等。其中,有限差分法雖然相對(duì)其他三種方法而言簡便易行,而且有豐富多樣的離散方法,但是它在求解問題時(shí)對(duì)求解區(qū)域的適應(yīng)性比較差。有限元法雖然采用的網(wǎng)格剖分更加靈活,從一定程度上講對(duì)求解區(qū)域具有更強(qiáng)的適應(yīng)性,但是它在求解間斷問題時(shí)會(huì)受到很大的限制,達(dá)不到有限差分法的效果。而有限體積法可以被視為是上述兩種方法的結(jié)合,雖然能夠充分利用有限元網(wǎng)格靈活性和克服差分法對(duì)網(wǎng)格適應(yīng)性差的缺陷,但是數(shù)值實(shí)驗(yàn)較難進(jìn)行[1]。作為一種新興的數(shù)值模擬方法,LBM基于Boltzmann方程的離散,是一種自下而上的求解方法。它描述了微觀粒子的碰撞和遷移,利用分布函數(shù)(一種概率密度分布函數(shù))來確定粒子的分布,即分布函數(shù)描述了流體的宏觀運(yùn)動(dòng)。近年來,由于LBM具有計(jì)算簡便、良好的并行性、處理不規(guī)則的復(fù)雜邊界容易且對(duì)于源項(xiàng)的考慮簡單等諸多優(yōu)勢(shì),已經(jīng)自然而然地發(fā)展成為了求解淺水波方程的一種新方法[2]。

在以往的研究中,英國利物浦大學(xué)的教授Zhou[3]較為全面地闡述了淺水波方程的LBM理論,包括外力的不同處理格式、湍流模型的構(gòu)造、多種邊界條件的處理方法以及對(duì)于許多經(jīng)典淺水波問題的驗(yàn)證。中山大學(xué)環(huán)境科學(xué)與工程學(xué)院的Li和Huang[4]進(jìn)行了對(duì)流-擴(kuò)散方程與淺水波方程耦合的研究,并采用LB的多松弛模型和雙松弛模型分別對(duì)流場(chǎng)和污染物場(chǎng)進(jìn)行了模擬。文獻(xiàn)[5]提出了一類粘性淺水方程的晶格Boltzmann (LB)模型,該模型采用源項(xiàng)的二階矩來恢復(fù)控制方程中的粘性,并消除Chapman-Enskog分析過程中產(chǎn)生的附加誤差。文獻(xiàn)[6]建立了一種適用于淺水方程的晶格玻爾茲曼模型(LABSWE),它用源項(xiàng)如床面坡度,床面摩擦力來求解方程。通過求解定常和非定常流動(dòng)問題,驗(yàn)證了該模型的有效性。

鑒于以上背景,文章首先給出了格子Boltzmann方法(LBM)的基本理論,然后利用經(jīng)典的一維五速度(D1Q5)的離散速度模型,給出KdV方程中含有修正項(xiàng)的格子Boltzmann(LB)模型推導(dǎo)公式,最后進(jìn)行數(shù)值模擬,將KdV方程的精確解和模擬解進(jìn)行比較,然后驗(yàn)證修正模型的精確性。實(shí)驗(yàn)結(jié)果表明,用格子Boltzmann方法對(duì)KdV方程進(jìn)行求解,其模擬解和精確解吻合度較高。

1 方法

格子玻爾茲曼方法(Lattice Boltzmann Method)是一種基于微觀介觀的流體力學(xué)計(jì)算方法,適用于二維或三維流體流動(dòng)問題的模擬[7]。圖1是格子玻爾茲曼方法的流程圖:其主要思想是離散化流體的分布函數(shù),通過對(duì)分布函數(shù)的演化來模擬流體的運(yùn)動(dòng)。近年來,LBM由于計(jì)算簡單、并行性好、易于處理復(fù)雜不規(guī)則的邊界及能簡單方便地考慮源項(xiàng)等優(yōu)勢(shì),已經(jīng)發(fā)展成為求解淺水波方程(SWEs)的一種新方法。下面是格子玻爾茲曼方法的求解步驟:

圖1 LBM方法流程圖

(1)確定格子和速度模型:首先需要確定流場(chǎng)的離散化格點(diǎn)和速度模型。通常情況下,將流體分成若干個(gè)小區(qū)域,每個(gè)小區(qū)域都對(duì)應(yīng)一個(gè)格點(diǎn),格點(diǎn)上有一組離散的速度向量。

(2)定義分布函數(shù):為了描述格點(diǎn)上流體的狀態(tài),引入一個(gè)分布函數(shù)g,用來表示在每個(gè)格點(diǎn)上,每個(gè)速度方向上的粒子數(shù)密度。它是時(shí)間和位置的函數(shù),通常用離散的速度和離散的時(shí)間步長表示。

(3)離散Boltzmann方程:基于Boltzmann方程,對(duì)分布函數(shù)進(jìn)行離散化,得到離散化的Boltzmann方程,它描述了分布函數(shù)的演化過程。在格子玻爾茲曼方法中,Boltzmann方程可以看成是一個(gè)簡單的微分方程,其左側(cè)是分布函數(shù)的時(shí)間導(dǎo)數(shù),右側(cè)是一個(gè)碰撞項(xiàng)和一個(gè)弛豫項(xiàng),用于描述粒子之間的相互作用和粒子與流體之間的相互作用。

(4)離散碰撞項(xiàng)和弛豫項(xiàng):將碰撞項(xiàng)和弛豫項(xiàng)進(jìn)行離散化,得到離散化的碰撞算子和弛豫算子。碰撞算子用于描述粒子之間的相互作用,而弛豫算子用于描述粒子與流體之間的相互作用。

(5)迭代求解:通過迭代求解離散化的Boltzmann方程,計(jì)算出每個(gè)格點(diǎn)上的分布函數(shù),從而得到流場(chǎng)的速度場(chǎng)和密度場(chǎng)。

(6)計(jì)算宏觀量:根據(jù)格點(diǎn)上的分布函數(shù),可以計(jì)算出宏觀量,如速度、密度、壓力等。

(7)處理邊界條件:對(duì)于邊界處的格點(diǎn),需要根據(jù)具體的物理問題設(shè)置邊界條件。

(8)模擬結(jié)束:當(dāng)達(dá)到預(yù)設(shè)的模擬時(shí)間或達(dá)到收斂條件時(shí),模擬結(jié)束。

2 模型簡介

2.1 離散速度模型

LBM中一維離散速度模型最常見的是D1Q3模型和D1Q5模型[8],具體如下:

(1)對(duì)于D1Q3模型(見圖2),模型參數(shù)如下:

圖2 D1Q3離散速度模型

(1)

其中,ωi為權(quán)重系數(shù),c=Δx/Δt為粒子遷移速度,cs是與當(dāng)?shù)芈曀傧嚓P(guān)的量。

(2)對(duì)于D1Q5模型(見圖3),模型參數(shù)如下:

圖3 D1Q5離散速度模型

(2)

(3)

其中,ωi為權(quán)重系數(shù),c=Δx/Δt為粒子遷移速度,cs是與當(dāng)?shù)芈曀傧嚓P(guān)的量。

2.2 KdV方程

將LB模型應(yīng)用于KdV方程中,需要將KdV方程離散化成網(wǎng)格上的方程組,然后通過LB模型求解這個(gè)方程組。具體來說,LB模型中的速度分布函數(shù)被定義為格點(diǎn)上的波高,通過計(jì)算速度分布函數(shù)在不同時(shí)間和空間的演化來模擬KdV方程的行為。與傳統(tǒng)的有限差分法和有限元法相比,LB模型具有計(jì)算效率高、適合并行計(jì)算等優(yōu)點(diǎn),因此在模擬非線性波等問題時(shí)得到了廣泛應(yīng)用。

考慮非線性偏微分方程一般形式[9]:

(4)

其中,u=u(x,t)是物質(zhì)在空間x處和時(shí)刻t時(shí)的密度,α,β,γ,δ為參數(shù)。

當(dāng)β=0,γ=0時(shí),方程(4)化為KdV方程:

(5)

3 模型推導(dǎo)

采用D1Q5模型給出KdV方程含有修正項(xiàng)的LB模型推導(dǎo),給出 的演化方程為:

(6)

(7)

(8)

其中,λ1,i,λ2,i和pi為調(diào)整參數(shù)。宏觀變量u(x,t)定義為[13-15]:

(9)

為了得到穩(wěn)定的宏觀變量u,假設(shè)分布函數(shù)gi也處于平衡狀態(tài),且有:

(10)

由(10)可得,

(11)

(12)

(13)

使用多尺度分析將方程恢復(fù)到宏觀方程,引入1個(gè)離散的時(shí)間尺度和3個(gè)連續(xù)的時(shí)間尺度,其具體表達(dá)形式為:

t0,t1=εt,t2=ε2t,t3=ε3t

(14)

對(duì)分布函數(shù)gi和時(shí)間導(dǎo)數(shù)進(jìn)行Chapman-Enskog展開,可得:

(15)

(16)

其中,ε表示任意小的參數(shù),在宏觀方程的推導(dǎo)過程中,不妨假設(shè)ε=Δt,將(6)式的左邊對(duì)時(shí)間和空間進(jìn)行泰勒展開,并保留Ο(ε4)項(xiàng),可得

(17)

將(15)和(16)代入(17)式中得,并對(duì)比左右兩邊可得ε的同階項(xiàng):可以得出,

O(ε0)系數(shù):

(18)

O(ε1)系數(shù):

(19)

O(ε2)系數(shù):

(20)

O(ε3)系數(shù):

(21)

(22)

(23)

同理,結(jié)合約束條件(12)和(13),將方程(20)式兩邊分別乘以ci后并對(duì)i求和,得出:

(24)

結(jié)合(7)(8)(9)和(19),將方程(16)兩邊對(duì)i求和,得出:

(25)

同理,結(jié)合(10)(11)(12)和(22)～(24),將方程(20)兩邊對(duì)i求和,得出:

(26)

將(3.19)×ε+(3.20)×ε2,可得:

(27)

將方程(27)和(4)對(duì)比可得:

α=2cτλ1ε,β=(n+1)cτλ2ε

(28)

(29)

可以得出,方程(27)就是一維KDV方程的LB模型。由方程(7)和(13)式,得出修正函數(shù)hi為:

(30)

(31)

4 數(shù)值模擬

將LB模型應(yīng)用于KdV方程中,需要將KdV方程離散化成網(wǎng)格上的方程組,然后通過LB模型求解這個(gè)方程組。具體來說,LB模型中的速度分布函數(shù)被定義為格點(diǎn)上的波高,通過計(jì)算速度分布函數(shù)在不同時(shí)間和空間的演化來模擬KdV方程的行為。與傳統(tǒng)的有限差分法和有限元法相比,LB模型具有計(jì)算效率高、適合并行計(jì)算等優(yōu)點(diǎn),因此在模擬非線性波等問題時(shí)得到了廣泛應(yīng)用。

(1)考慮如下的KdV方程:

取c1=2c,c2=c3=1,得出該方程的精確解為:

模擬結(jié)果如圖4、圖5所示。圖4和圖5分別給出了t=0.01和t=0.25時(shí)刻的LB模擬解和解析解的對(duì)比圖,從圖中可以看出:在t=0.25之前,模擬解和解析解吻合的程度較高,但是隨之時(shí)間的推移,模擬解與解析解存在一定的偏離,這主要是原因有擾動(dòng)項(xiàng)O(ε4),它在一定程度會(huì)對(duì)孤子高度、速度以及形狀有影響,且當(dāng)t>0.25時(shí),方程的模擬解和精確解差別較大。

圖4 t=0.01,LB模擬解和精確解對(duì)比 圖5 t=0.25,LB模擬解和精確解對(duì)比

(2)考慮如下的KdV方程:

其中,u=u(x,t),u為波動(dòng)地振幅,x為波橫向傳播的位移,t為時(shí)間。

初始條件為:

u(x,0)=3Asech2(Bx+C),x∈[0,2]

邊界條件為:

u(0,t)=u(2,t)=0,t>0

解析解為:

u(x,t)=3Asech2(Bx-Dt+C),x∈[0,2]

設(shè)置參數(shù)如下:Δx=0.001,Δt=5×10-4,α=1,δ=4.84×10-4,模擬結(jié)果如圖(6)和圖(7)所示.

圖6 t=0.000 5時(shí),模擬解與解析解對(duì)比 圖7 t=0.002 0時(shí),模擬解與解析解對(duì)比

圖6和圖7分別給出了t=0.0005和t=0.002時(shí)刻的LB模擬解和解析解的對(duì)比圖,通過對(duì)該方程的模擬結(jié)果進(jìn)行分析,發(fā)現(xiàn)在t=0.002之前,模擬解和解析解吻合的程度較高,但是隨著時(shí)間的推移,模擬解與解析解存在一定的偏離,主要的原因是含有擾動(dòng)項(xiàng)O(ε4),當(dāng)然也可能是該類波的傳播速度極快,長時(shí)間模擬就會(huì)產(chǎn)生偏差。表1給出了該方程在不同時(shí)刻的誤差。

表1 方程在不同模擬時(shí)刻的誤差比較

從表1可以看出,在LB模型下得到的模擬解和解析解非常逼近,無論是L∞誤差,還是均方根誤差L2和整體相對(duì)誤差GRE,其兩者之間的誤差數(shù)量級(jí)都達(dá)到了 ,說明該數(shù)值結(jié)果是比較理想的。

5 結(jié)論

現(xiàn)在,微分方程無處不在,各個(gè)科學(xué)領(lǐng)域的研究都伴隨著微分方程模型。由于實(shí)際生活中的微分方程模型形式日趨復(fù)雜,為了與實(shí)際問題相匹配,微分方程解的形式越來越多樣化。本文對(duì)兩個(gè)特殊的KDV方程,利用格子Boltzmann模型求解并與其精確解進(jìn)行比較,得出使用格子Boltzmann方法對(duì)非線性偏微分方程求解取得了較好的效果。在未來的工作中,將嘗試?yán)^續(xù)改進(jìn)格子Boltzmann模型,并對(duì)更加復(fù)雜的偏微分方程或者淺水波方程進(jìn)行模擬。希望本文可以為其他學(xué)者在求解偏微分方程方面的研究工作提供一定的參考價(jià)值。