劉 芳， 趙海龍

（長春工業(yè)大學(xué) 基礎(chǔ)科學(xué)學(xué)院，吉林 長春 130012）

0 引 言

格子Boltzmann方法（Lattice Boltzmann Method，LBM）是20世紀(jì)80年代誕生的一種新興的計算流體力學(xué)方法，它植根于物理學(xué)描述氣體運動的連續(xù)Boltzmann方程以及計算機理論的元胞自動機模型。與傳統(tǒng)方法不同，該方法是從微觀動力學(xué)角度出發(fā)，把宏觀物理量當(dāng)作微觀量的統(tǒng)計平均結(jié)果［1－2］。格子Boltzmann方法誕生以來，已在許多流體問題的數(shù)值計算上取得了很大成功。近年來，一些學(xué)者將格子Boltzmann方法用于求解其他的線性和非線性偏微分方程，也取得了滿意的成果。目前該方法已成為偏微分方程數(shù)值解領(lǐng)域除傳統(tǒng)方法外的又一個可供選擇的新興數(shù)值算法［3－9］。

文中考慮用LBM求解含有非線性源項和非線性擴散項的多孔介質(zhì)方程［10］：

式中，m和k是有理數(shù)，a和b是參數(shù)，u＝u（x，t），x和t分別表示空間和時間變量，D（u）＝um是非線性擴散項。這類方程經(jīng)常出現(xiàn)在熱和質(zhì)量傳遞、燃燒理論和多孔介質(zhì)流動等非線性問題中。例如，該方程可以描述在靜止的介質(zhì)中，溫度是以一個冪率函數(shù)變化的熱傳遞過程，一些學(xué)者已經(jīng)用其它方法研究過該方程，例如緊致有限差分方法［11］和Adomian分解法［12］等。但是到目前為止，還沒有應(yīng)用LBM對該方程的相關(guān)數(shù)值研究，文中將對方程（1）構(gòu)建一個簡單實用的LBM求解模型。

目前LBM已被應(yīng)用于解決大量的線性和非線性偏微分方程。然而，一般的格子Boltzmann模型在處理非線性擴散項上還有一定的困難。文中通過靈活的選取平衡分布函數(shù)，使這一問題得到了很好的解決。

1 多孔介質(zhì)方程的格子Boltzmann模型

文中考慮的多孔介質(zhì)方程（1）可以改寫為：

令局部粒子分布函數(shù)的演化方程形式如下：

Si（x，t）——源項分布函數(shù)；

i——離散速度集合｛c0，c1，c2，…，cN｝的指標(biāo)。

通過實施Taylor展開，Chapman－Enskog展開以及多尺度分析技術(shù)，可以從方程（3）恢復(fù)出具有二階精度的方程（2）。這一過程的實施需滿足以下幾個對平衡態(tài)分布函數(shù)和源項分布函數(shù)的限制條件：

由式（4）～式（6）解得平衡態(tài)分布函數(shù)分別為：

為了滿足式（7）和式（8），取

源項F（x，t）＝buk只針對文中研究的多孔介質(zhì)方程。方程（1）和方程（2）中的源項可以是u，x和t的其它函數(shù)形式。

2 數(shù)值模擬

為了檢驗?zāi)Ｐ偷木群陀行?，我們將對幾個多孔介質(zhì)方程進(jìn)行數(shù)值模擬。定義時間步t＝tj的全局相對誤差（GRE）和最大絕對誤差（MAE）為：

式中：unum（xi，tj），u（xi，tj）——分別為tj時刻點xi處的數(shù)值解和精確解；

宏觀的初邊值條件由精確解確定，對于微觀的初始條件，令局部粒子分布函數(shù)直接等于局部平衡分布函數(shù)，即。演化方程（3）中導(dǎo)數(shù)項的計算將采用顯示差分格式?tSi（x，t）＝［Si（x，t）－Si（x，t－Δt）］／Δt計算。邊界處理采用Z L Guo［13］的非平衡外推法來確定邊界處的粒子分布函數(shù)。

在非線性科學(xué)中，多孔介質(zhì)方程是一個非常重要的偏微分方程。下面的例子都是來自于實際物理問題。

例1［10］：取m＝1，k＝0，則方程（1）變成

其精確解為

現(xiàn)在取a＝1，b＝1，按照方程（2）改寫方程（9）為

例2［10］：取a＝1，m＝－1和b＝0，則方程（1）變成

其精確解為

其中c1和c2都為任意常數(shù)。簡單選取c1＝1，c2＝10，因此精確解變成

按照方程（2）改寫方程（10）為

例3［10］：取和b＝0，則方程（1）變成

其精確解為

其中，c1和c2為任意常數(shù)。簡單選取c1＝1，c2＝15，得到

按照方程（2）改寫方程（11），有

例4［10］：取和則方程（1）變成

其精確解為

按照方程（2）改寫方程（12），有

在數(shù)值模擬中，例1～例3選取計算域為［0，1］，例4選取計算域為［2，3］。4個算例的時間步長均為Δt＝0.000 1，空間步長均為Δx＝0.01，自由參數(shù)τ的選取依據(jù)使t＝1時刻GRE達(dá)到最優(yōu)的原則。例1～例4的τ分別為8.062，9.998，14.002，0.999。在時間t＝1，2，3，4時刻，4個算例的格子Boltzmann數(shù)值解與精確解之間的GRE和MAE分別見表1和表2。

表1 4個算例的不同時間的全局相對誤差比較

表2 4個算例的不同時間的最大絕對誤差比較

表中的數(shù)據(jù)表明數(shù)值解與精確解吻合的非常好，顯示了該模型的有效性。

為了檢驗該方案的實際空間精度，分別在網(wǎng)格數(shù)目為16，32，64和128（即Δx＝0.062 5，0.031 25，0.015 625，0.007 812 5）時對4個算例的方程進(jìn)行數(shù)值計算，且保持其它參數(shù)不變。t＝1時刻的GRE和MAE與空間步長Δx在對數(shù)坐標(biāo)系下的擬合圖如圖1所示。

圖1 t＝1時刻的空間精度檢測

圖中所有擬合直線的斜率都在2.0左右，這就驗證了該模型在實際計算中達(dá)到了空間二階精度。

3 結(jié) 語

對含源項的非線性多孔介質(zhì)方程構(gòu)造了一個簡單實用的格子Boltzmann模型，通過改造演化方程以及對局部平衡態(tài)分布函數(shù)的一些限制，恢復(fù)出具有二階精度的宏觀方程。除時間步長和空間步長外，計算格式中只有一個自由參數(shù)τ。采用D1Q3速度模型，編程簡單，計算速度快。對幾個典型的多孔介質(zhì)方程進(jìn)行了數(shù)值模擬，驗證了該模型的精度和有效性。該模型同樣適用于同類型的其他方程，且可以直接擴展到二維和三維問題。

［1］ R Benzi，S Succi，M Vergassola.The lattice Boltzmann equation：theory and applications［J］.Phys.，1992，222：145－197.

［2］ Y H Qian，S Succi，S A Orszag.Recent advances in lattice Boltzmann computing［J］.Annu.Rev.Comput.Phys.，1995（3）：195－242.

［3］ B C Shi，Z L Guo.Lattice Boltzmann model for nonlinear convection－diffusion equations［J］.Phys.Rev.E.，2009，79：16701.

［4］ I Ginzburg.Truncation errors，exact and heuristic stability analysis of tw－relaxation－times lattice Boltzmann schemes for anisotropic advection－diffusion equation［J］.Commun.Comput.Phys.，2012（11）：1439－1052.

［5］ F Liu，W P Shi.Numerical solutions of two－dimensional Burgers＇equations by lattice Boltzmann method［J］.Commun.Non.Sci.Numer.Simul.，2011，16：150－157.

［6］ H L Lai，C F Ma.A new lattice Boltzmann model for solving the coupled viscous Burgers＇equation［J］.Physica A，2014，395：445－457.

［7］ J Y Zhang，G W Yan.Numerical studies based on higher－order accuracy lattice Boltzmann model for the complex Ginzburg－Landau equation［J］.J.Sci.Comput.，2012，52：656－674.

［8］ B C Shi.Lattice Boltzmann simulation of some nonlinear complex equations［J］.Lect.Notes Comput.Sci.，2007，4487：818－825.

［9］ F F Wu，W P Shi，F(xiàn) Liu.A lattice Boltzmann model for the Fokker－Planck equation［J］.Commun.Non.Sci.Numer.Simul.，2012，17：2776－2790.

［10］ A D Polyanin，V F Zaitsev.Handbook of nonlinear partial differential equations［M］.Boca Raton：Chapman and Hall／CRC Press，2004.

［11］ N Pamuk.Series solution for porous medium equation with a sourceterm by adomian＇s decomposition method［J］.Appl.Math.Comput.，2006，178：480－485.

［12］ M Sari.Solution of the porous media equation by a compact finitedifference method［J］.Math.Prob.Eng.，2009，17：1－13.

［13］ Z L Guo，C G Zheng，B C Shi.Non－equilibrium extrapolation method forvelocity and pressure boundary conditions in the lattice Boltzmannmethod［J］.Chin.Phys.，2002，11：366－374.