廣義Ginzburg-Landau方程的格子Boltzmann模型?

李精偉,劉德民,張國梁,林玉婷

(新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，烏魯木齊 830046)

1 引言

格子Boltzmann模型是近二十年來發(fā)展起來的一類流體系統(tǒng)建模和仿真的新方法，在理論和應(yīng)用方面已取得了迅速的發(fā)展，并逐漸成為相關(guān)領(lǐng)域的研究熱點之一．格子Boltzmann模型是介于流體的微觀分子動力學(xué)模型和宏觀連續(xù)模型之間的介觀模型，它兼具二者的優(yōu)點，目前在計算流體動力學(xué)等問題中得到廣泛的應(yīng)用[1]．

本文研究廣義Ginzburg-Landau方程，其數(shù)學(xué)形式為

這里u=u(x,t)表示人口密度，是關(guān)于空間變量x與時間變量t≥0的函數(shù)；g(u)是給定的非線性函數(shù)，稱為動力項或反應(yīng)項，典型的取法有：g(u)=c(1?u2)2,c>0,g(u)=0，或g(u)=|u|p?1u；α,β,γ是正實數(shù)；f=f(x,t)表示人口擾動函數(shù)．該方程首先在研究人口增長與彌散過程中提出[2]，并進一步在非線性方程初邊值理論研究中得到了廣泛的應(yīng)用[3-7]；文獻[3–5]證明了方程(1)的Cauchy問題和初邊值問題整體廣義解和整體古典解的存在唯一性；文獻[6]證明了方程(1)時間周期問題廣義解和古典解的存在唯一性；文獻[7]建立了方程(1)最優(yōu)邊界控制問題解的存在性；文獻[8]研究了方程(1)在g(u)=c(1?u2)2,c>0情形下的行波解在H2攝動意義下的非線性不穩(wěn)定性；另外，文獻[9]利用廣義Herm ite函數(shù)作為基函數(shù)，建立了方程(1)的譜逼近格式及誤差估計理論．

注意到方程(1)是具有高階耗散性的非線性對流-擴散-傳導(dǎo)問題，如果采用有限元方法離散，需要選取高階形函數(shù)；如果采用有限差分法求解，會導(dǎo)致所得代數(shù)方程帶寬較大，兩者都需要較大的求解自由度，另外方程(1)是非線性方程，需要采用較好的非線性迭代方法才能滿足計算要求．理論上，格子Boltzmann模型可以恢復(fù)到任意微分方程，而模型本身卻是顯式推進的，因此可以大大簡化理論分析和程序?qū)崿F(xiàn)[1]．基于如上考慮，本文擬構(gòu)造方程(1)的高階形式的格子Boltzmann方法[10]，這里的高階即選擇多個速度方向，通過采用Chapman-Enskog多尺度展開可以恢復(fù)出方程的高階導(dǎo)數(shù)項和非線性項．

本文安排如下：在第2節(jié)建立了帶有附加項的單弛豫格子Boltzmann模型，并通過采用Chapman-Enskog多尺度展開方法及Taylor公式得到了平衡態(tài)分布函數(shù)與附加項所滿足的方程，通過選取不同的參數(shù)，給出平衡態(tài)分布函數(shù)與附加項的具體形式，得到截斷誤差和穩(wěn)定性條件；第3節(jié)給出了部分?jǐn)?shù)值結(jié)果．

2 單弛豫格子Boltzm ann模型

考慮一維情形，即x∈[a,b]，通過引入一系列均勻分布的格點xj,j=0,1,···,N，將區(qū)間[a,b]均分成一系列格子[xj,xj+1]，網(wǎng)格節(jié)點間的連線構(gòu)成向量集(e0,e1,e2,e3,e4)=(0,1,?1,2,?2)．進一步引進粒子分布函數(shù)fi(x,t)，代表t時刻在格點x處ei方向上的粒子分布函數(shù)．規(guī)定宏觀物理量u滿足如下限定關(guān)系

利用公式(2)，對物理量u求解轉(zhuǎn)化為尋求fi(x,t)隨時間的演化．本文引入如下形式的單弛豫格子Boltzmann模型

其中(x,t)是局部平衡態(tài)分布函數(shù)；?t是時間步長，滿足0≤ ?t? 1；τ是馳豫時間，表示分布函數(shù)fi(x,t)趨近平衡態(tài)分布函數(shù)(x,t)的速度，穩(wěn)定性要求τ≥0.5；c=?x/?t為常數(shù)；hi(x,t)與gi(x,t)為附加函數(shù)，分別用來恢復(fù)非線性項與人口擾動函數(shù)．

2.1 廣義Ginzburg-Landau方程的恢復(fù)

為了有效的恢復(fù)出廣義Ginzburg-Landau方程，首先引入一個空間尺度變量x=x，以及五個時間尺度變量t1= εt,t2= ε2t,t3= ε3t,t4= ε4t,t5= ε5t，其中ε為克努森數(shù)，并假定克努森數(shù)等于時間步長，即ε=?t(參見文獻[8])．對分布函數(shù)fi(x,t)引入Chapman-Enskog多尺度展開，即

這里是粒子分布函數(shù)相對于局域平衡態(tài)分布函數(shù)的偏差，均為變量(x,t)的函數(shù)，在長波低頻情形下滿足條件

利用(2)式及(5)式可得

由鏈導(dǎo)法則可知

運用Taylor公式將(3)式兩邊展開至?t7項，則得到

將(4),(7)帶入(8)式，然后分別令左右兩端的?t的各階冪次的系數(shù)相等，可以得到關(guān)于時間ti的系列方程．

利用(6)式，并限定平衡態(tài)分布函數(shù)f(0)i及附加函數(shù)hi,gi滿足條件

上式就是格子Boltzmann模型對應(yīng)的宏觀方程．

為了由(17)式恢復(fù)出廣義Ginzburg-Landau方程，只需

成立，即可得到宏觀的廣義Ginzburg-Landau方程(1)．方程組(18)是參數(shù)τ,λ,η,θ,λ1的超定方程組．

當(dāng)β?=0時，選取λ為自由未知量，求解可得

當(dāng)β=0時，選取η是自由未知量，求解可得

2.2 平衡態(tài)分布函數(shù)及附加函數(shù)的確定

注意到(13)式中包含三個方程組，分別是平衡態(tài)分布函數(shù)及附加函數(shù)hi,gi所滿足的方程，其中只有平衡態(tài)分布函數(shù)的方程組是適定方程組，其它兩個為欠定方程組．通過求解可得

2.3 廣義Ginzburg-Landau方程的截斷誤差

為了獲得廣義Ginzburg-Landau方程的截斷誤差，我們需要從系列方程出發(fā)，利用高階矩，再恢復(fù)到廣義Ginzburg-Landau方程，便可得到誤差項．我們進行如下的計算：對(8)式繼續(xù)展開，可得到?t5的系列方程．

根據(jù)(21)式，得到相應(yīng)f(0)i,hi,gi的高階矩方程

我們作如下計算

并利用(13),(18)和(23)式，可得到修正的廣義Ginzburg-Landau方程

其中誤差項E4的表達式為

由Weierstrass定理可知，非齊次項可以用多項式逼近，多項式的首項為aub，即

a是一個常數(shù)，b是整數(shù)，則(26)式化簡可得

這樣，我們構(gòu)造了四階精度的廣義Ginzburg-Landau方程

2.4 模型的穩(wěn)定條件

式(29)可以看做是修正后的廣義Ginzburg-Landau方程，其中E4為余項．

為了不失一般性，使用Warm ing和Hyett[11]描述的方法來消除高階導(dǎo)數(shù)項，可以得到等同于數(shù)值格式的修正后的偏微分方程．廣義Ginzburg-Landau方程有如下格式

那么修正后的偏微分方程可以寫成

在上式中，Rs(u)和Rp(u)分別為數(shù)值格式Boltzmann法的數(shù)值耗散項和數(shù)值色散項，由于方程(1)中右端的每一項都屬于耗散項．故我們只考慮數(shù)值耗散項，其有如下的形式

其中m是正整數(shù)，ν2m是2m階的數(shù)值耗散項系數(shù)．選擇Rs(u)中的主要項，則化為

式中的λ為波數(shù)[12]．

根據(jù)Hirt啟示性穩(wěn)定性理論[13]，我們認(rèn)為，其主要項的形式為

該數(shù)值方法具有Hirt啟示性穩(wěn)定性的條件為

本文中的LBM中取m0=1時，(29)式中的二階耗散余項的主項為

則Hirt啟示性穩(wěn)定條件為

因此，格子Boltzmann方程(3)的控制條件為(37)式．

3 數(shù)值算例

本節(jié)給出三個算例說明算法的有效性．由(19)和(20)式知道參數(shù)λ與η是自由未知量，程序?qū)崿F(xiàn)時這兩個參數(shù)由自適應(yīng)方式給出，以保證參數(shù)τ∈(0.5,1)，并滿足(37)式的條件，對于Dirch let邊界條件的處理采用非平衡態(tài)外推方法處理[1]，取初始分布函數(shù)分別為問題的精確解與格子Boltzmann解．定義離散的lq誤差為另外當(dāng)時間步長充分小時，可以認(rèn)為誤差主部僅依賴于空間步長，從而定義lq誤差的的收斂階Rq為Rq=中為第i次實驗的lq誤差．

算例1 首先考慮如下線性問題

該算例是具有四階耗散項的拋物問題，問題的精確解為u=sin(4x)cos(4t)，右端項f可由精確解給出．在計算中選取α=0.000001,β=0.02．

表1是取固定時間步長?t=1/2000時，對于不同的空間步長?x=1/M所得數(shù)值結(jié)果．我們知道，當(dāng)時間步長充分小時，誤差的主要來源是空間離散的誤差．在表1中，隨著空間步長逐次加密，即?x=1/100→ 1/800，各類誤差均不斷減小，并且收斂階Rq不斷上升，至少達到4階；但隨著空間步長繼續(xù)加密，即?x=1/900，l1誤差可以達到8階，同時所對應(yīng)的誤差可以達到7.8階；當(dāng)?x=1/1000時，l1誤差仍然可以保證8階收斂，l2誤差卻出現(xiàn)掉階；l∞誤差在?x=1/900,1/1000的收斂階Rq出現(xiàn)負(fù)值．通過檢查程序結(jié)果發(fā)現(xiàn)，這兩種情形l∞誤差均出現(xiàn)在緊鄰右側(cè)邊界的一列位置上，進一步研究發(fā)現(xiàn)主要是因為我們所選取的附加函數(shù)hi中包含u3項，該項的作用是是用來恢復(fù)原方程中的非線性對流項，在使用非平衡態(tài)外推邊界時需要用到兩側(cè)邊界以外位置即虛擬邊界的值，程序?qū)崿F(xiàn)時用到了外推的方式，所以精度損失較多．

表1:T=1時，不同空間步長?x=1/M所對應(yīng)的誤差

表2是取固定空間步長?x=1/800時，對于不同的時間步長?t=1/N所得數(shù)值結(jié)果．我們知道，當(dāng)空間步長充分小時，誤差的主要來源是時間離散誤差．但通過表2發(fā)現(xiàn)，隨著時間步長逐次加密，即?t=1/600→1/1200，l1和l2誤差均不斷減小，對應(yīng)的收斂階Rq均可增加到5階；同時對應(yīng)的l∞也在減小，變化不大，但隨著時間步長?t繼續(xù)減小，l1和l2誤差均明顯增加，l∞誤差卻變化不大，對應(yīng)的收斂階均掉階并出現(xiàn)負(fù)階．結(jié)合表1可知，此時空間離散誤差起主導(dǎo)作用．

表2:T=1時，不同時間步長?t=1/N所對應(yīng)的誤差

結(jié)合表1和表2可知，當(dāng)空間步長和時間步長充分小時，對于算例1這樣的線性問題，格子Boltzmann方法的空間與時間l1和l2誤差均很好，對應(yīng)的空間收斂階可以達到8階，時間收斂階可以達到5階；同時l∞誤差保持較好，對應(yīng)的空間收斂階可以到到3階，時間收斂階可以達到2階．

算例2 考慮如下形式的非線性對流擴散方程

此問題的精確解為u=t ex，關(guān)于變量x具有指數(shù)增長性，右端項f可由方程直接求出．計算中選取β=0.02,γ=0.002,?t=1/2000，從表3可知各類誤差的收斂階均達到2階，具有很好的收斂性．

表3:T=0.5時，不同空間步長?x=1/M所對應(yīng)的誤差

算例3 考慮如下帶有高階耗散性的非線性傳導(dǎo)擴散問題

這里選取g(u)=c(1?u2)2,c>0．該問題的精確解為u=sin(πx)cos(t)．算例中α=0.0002,β=0.005,c=0.005,?x=1/60,?t=1/80．圖1與圖2給出了精確解與數(shù)值解的曲面圖，圖3給出了當(dāng)T=1時0.25T,0.5T,0.75T以及T時刻精確解與數(shù)值解比較圖，圖4顯示了這四個時刻的誤差曲線圖，其中高度方向為誤差的對數(shù)值．

圖1:精確解

圖2:數(shù)值解

圖3:精確解與數(shù)值解的比較

圖4:精確解與數(shù)值解的誤差比較

從圖1至圖4可知數(shù)值解具有很好的逼近性質(zhì)，從圖4可知邊界附近的誤差較大，這其中的主要原因是g(u)具有較高的非線性，邊界處理時采用的非平衡態(tài)外推方法理論上只能達到二階精度，因此造成邊界附近誤差較大．

4 結(jié)論

本文提出了求解廣義Ginzburg-Landau方程單弛豫格子Boltzmann模型，通過理論分析確定了模型的正確性．?dāng)?shù)值算例1說明了對于線性問題，所提出的格子Boltzmann模型具有高效性；數(shù)值算例2表明所提出的模型同樣適用于非線性對流擴散方程；算例3表明所提出模型在高階耗散傳導(dǎo)方程數(shù)值模擬中的有效性；同時發(fā)現(xiàn)當(dāng)方程的非線性較強時，即使采用具有二階精度的非平衡態(tài)外推邊界處理方法，在靠近邊界時，誤差仍然較大，因此格子Boltzmann方法邊界處理的精度對數(shù)值解的穩(wěn)定性與收斂性有較大影響．

致謝：作者真誠感謝審稿人和編輯部提出寶貴的修改意見！

參考文獻：

[1]郭照立，鄭楚光.格子Boltzm ann方法的原理及應(yīng)用[M].北京：科學(xué)出版社，2009 Guo Z L,Zheng C G.Theory and App lications of Lattice Boltzm ann method[M].Beijing:Science Press,2009

[2]Cohen D S,murray JD.A generalized diff usion model for grow th and dispersal in a popu lation[J].Journal of mathem atical Biology,1981,12(2):237-249

[3]Liu B P,Pao C V.Integral rep resentation of generalized d iff usion model in popu lation prob lem s[J].Journal of Integral Equation,1984,6(2):175-185

[4]Chen G W.C lassical global solu tion of the initial boundary value prob lem s for a class of non linear parabolic equations[J].Comm entationes mathem aticae Universitatis Carolinae,1994,35(3):431-443

[5]陳國旺，孫和生.廣義G inzburg-Landau型非線性高階拋物型方程的Cauchy問題[J].數(shù)學(xué)年刊，1994,15A(4):485-492 Chen G W,Sun H S.Cauchy prob lem of a generalized Ginzbu rg-Landau higher order non linear parabolic equation[J].Chinese Annals of mathem atics,1994,15A(4):485-492

[6]王艷萍，陳國旺.人口問題中一廣義Ginzbu rg-Landau模型方程的時間周期解[J].數(shù)學(xué)物理學(xué)報，2006,26(2):258-266 Wang Y P,Chen G W.Tim e-periodic solution to a generalized G inzburg-Landau model equation in popu lation prob lem s[J].Acta mathem atica Scientia,2006,26(2):258-266

[7]Zhao X P,Duan N,Liu B.Op tim al control prob lem of a generalized Ginzbu rg-Landau model equation in popu lation prob lem s[J].mathem atical methods in the App lied Sciences,2014,37(3):435-446

[8]Liu C C.Instability of traveling waves for a generalized diff usion model in popu lation prob lem s[J].E lectronic Journal of Qualitative Theory of differential Equations,2004,18:1-10

[9]X iang X M,W ang ZQ.Generalized Herm ite spectralm ethod and itsapp lications to prob lem s in unbounded dom ains[J].SIAM Jou rnal on Num erical Analysis,2010,48(4):1231-1253

[10]Yan G W.A higher-order mom ent method of the lattice Boltzm ann model for the conservation law equation[J].App lied mathem atical modelling,2010,34(2):481-494

[11]W arm ing R F,Hyett B J.The mod ified equation app roach to the stability and accuracy analysis of finite diff erencem ethod[J].Journal of Com putational Physics,1974,14(2):159-179

[12]W ang J,Liu R X.A new app roach to design high-order schem es[J].Jou rnal ofComputational and App lied mathem atics,2001,134(1):59-67

[13]Hirt C W.Heu ristic stability theory for finite-d iff erence equations[J].Jou rnal of Com putational Physics,1968,2(4):339-355