KdV方程的格子Boltzmann模型求解

陳夢涵,王希胤,李金

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

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

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

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

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

1 方法

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

圖1 LBM方法流程圖

(1)確定格子和速度模型:首先需要確定流場的離散化格點和速度模型。通常情況下,將流體分成若干個小區域,每個小區域都對應一個格點,格點上有一組離散的速度向量。

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

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

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

(5)迭代求解:通過迭代求解離散化的Boltzmann方程,計算出每個格點上的分布函數,從而得到流場的速度場和密度場。

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

(7)處理邊界條件:對于邊界處的格點,需要根據具體的物理問題設置邊界條件。

(8)模擬結束:當達到預設的模擬時間或達到收斂條件時,模擬結束。

2 模型簡介

2.1 離散速度模型

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

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

圖2 D1Q3離散速度模型

(1)

其中,ωi為權重系數,c=Δx/Δt為粒子遷移速度,cs是與當地聲速相關的量。

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

圖3 D1Q5離散速度模型

(2)

(3)

其中,ωi為權重系數,c=Δx/Δt為粒子遷移速度,cs是與當地聲速相關的量。

2.2 KdV方程

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

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

(4)

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

當β=0,γ=0時,方程(4)化為KdV方程:

(5)

3 模型推導

采用D1Q5模型給出KdV方程含有修正項的LB模型推導,給出 的演化方程為:

(6)

(7)

(8)

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

(9)

為了得到穩定的宏觀變量u,假設分布函數gi也處于平衡狀態,且有:

(10)

由(10)可得,

(11)

(12)

(13)

使用多尺度分析將方程恢復到宏觀方程,引入1個離散的時間尺度和3個連續的時間尺度,其具體表達形式為:

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

(14)

對分布函數gi和時間導數進行Chapman-Enskog展開,可得:

(15)

(16)

其中,ε表示任意小的參數,在宏觀方程的推導過程中,不妨假設ε=Δt,將(6)式的左邊對時間和空間進行泰勒展開,并保留Ο(ε4)項,可得

(17)

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

O(ε0)系數:

(18)

O(ε1)系數:

(19)

O(ε2)系數:

(20)

O(ε3)系數:

(21)

(22)

(23)

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

(24)

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

(25)

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

(26)

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

(27)

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

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

(28)

(29)

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

(30)

(31)

4 數值模擬

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

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

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

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

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

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

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

初始條件為:

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]

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

圖6 t=0.000 5時,模擬解與解析解對比 圖7 t=0.002 0時,模擬解與解析解對比

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

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

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

5 結論

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