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

李精偉,劉德民,張國(guó)梁,林玉婷

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

1 引言

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

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

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

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

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

2 單弛豫格子Boltzm ann模型

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

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

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

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

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

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

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

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

運(yùn)用Taylor公式將(3)式兩邊展開(kāi)至?t7項(xiàng)，則得到

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

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

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

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

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

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

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

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

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

2.3 廣義Ginzburg-Landau方程的截?cái)嗾`差

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

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

我們作如下計(jì)算

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

其中誤差項(xiàng)E4的表達(dá)式為

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

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

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

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

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

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

那么修正后的偏微分方程可以寫(xiě)成

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

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

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

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

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

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

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

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

3 數(shù)值算例

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

算例1 首先考慮如下線性問(wèn)題

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

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

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

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

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

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

算例2 考慮如下形式的非線性對(duì)流擴(kuò)散方程

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

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

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

這里選取g(u)=c(1?u2)2,c>0．該問(wèn)題的精確解為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時(shí)0.25T,0.5T,0.75T以及T時(shí)刻精確解與數(shù)值解比較圖，圖4顯示了這四個(gè)時(shí)刻的誤差曲線圖，其中高度方向?yàn)檎`差的對(duì)數(shù)值．

圖1:精確解

圖2:數(shù)值解

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

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

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

4 結(jié)論

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

致謝：作者真誠(chéng)感謝審稿人和編輯部提出寶貴的修改意見(jiàn)！

參考文獻(xià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]陳國(guó)旺，孫和生.廣義G inzburg-Landau型非線性高階拋物型方程的Cauchy問(wèn)題[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]王艷萍，陳國(guó)旺.人口問(wèn)題中一廣義Ginzbu rg-Landau模型方程的時(shí)間周期解[J].數(shù)學(xué)物理學(xué)報(bào)，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