吳國(guó)忠，袁兆成，齊晗兵，李 棟，劉 杰

(東北石油大學(xué) 土木建筑工程學(xué)院，黑龍江 大慶 163318)

近年來(lái)，格子Boltzmann 方法(LBM)在求解偏微分方程領(lǐng)域發(fā)展很迅速，特別是求解Navier －Stokes 方程獲得很大成功。由于Boltzmann 方法具有物理圖像清晰、邊界條件容易處理以及并行性能好等優(yōu)點(diǎn)，被廣泛應(yīng)用于微/納米尺度流［1］、多孔介質(zhì)流［2］、多相多質(zhì)流［3］和晶體生長(zhǎng)［4］等許多領(lǐng)域。

擴(kuò)散過(guò)程在物理、化學(xué)、生物等許多領(lǐng)域中起著很重要的作用。熱擴(kuò)散方程是描述傳熱過(guò)程的一個(gè)重要方程，但在復(fù)雜的邊界和初始條件下，解析求解是很困難的。許多學(xué)者利用格子Boltzmann 方法求解擴(kuò)散問(wèn)題，并取得了很多成果。劉慕仁等人給出了求解一維有源擴(kuò)散方程的格子Boltzmann 模型，確定了局部平衡函數(shù)Chapman－Enskog 展開(kāi)的待定系數(shù)［5］，他們還利用格子Boltzmann 方法求解了一維對(duì)流擴(kuò)散方程，確定了方法中的粘滯系數(shù)與對(duì)流系數(shù)的關(guān)系［6］。徐世英等人利用濃度分布的Chapman－Enskog 展開(kāi)及多尺度技術(shù)，得出了Tyson 反應(yīng)擴(kuò)散的反應(yīng)物和催化劑隨時(shí)間的濃度空間分布值［7］。本文在前人的研究基礎(chǔ)上，提出了格子Boltzmann 方法求解一維熱擴(kuò)散方程的有效數(shù)值方法。

1 熱擴(kuò)散方程的Boltzmann 模型

一維熱擴(kuò)散方程表達(dá)式

式中

T——溫度，要求T≥0;

u——速度;

λ/ρcp——熱擴(kuò)散系數(shù);

λ——導(dǎo)熱系數(shù);

ρ——密度;

cp——定壓比熱容。

1．1 平衡態(tài)分布函數(shù)

將一維空間離散成均勻的線性格子，每個(gè)節(jié)點(diǎn)與相鄰的節(jié)點(diǎn)相連。將速度ca=cea離散為三個(gè)方向，ea可取值－1，0，1，分別對(duì)應(yīng)粒子速度方向向左、不動(dòng)和向右三種情況。c=Δx/Δt 為粒子遷移速率，3 －速格子模型如圖1 所示。

圖1 3 －速格子模型示意圖

用fi(x，t)表示沿ea方向的粒子分布函數(shù)，表示在時(shí)間t 時(shí)，在位置x 處粒子出現(xiàn)的概率，應(yīng)有fi≥0，宏觀量T 定義為

式中 feq

i ——平衡態(tài)分布函數(shù)。

根據(jù)統(tǒng)計(jì)物理，分布函數(shù)遵守以下格子Boltzmann 方程

式中 τ——弛豫時(shí)間。

式(3)也稱LBGK(Lattice BGK，簡(jiǎn)稱LBGK)方程，為滿足穩(wěn)定性條件，要求τ≥0．5。

由式(3)可知，只要得到feq(x，t)就可以得出以后時(shí)間的分布函數(shù)f(x，t)。由LBGK 方程恢復(fù)宏觀方程的關(guān)鍵是選擇適當(dāng)?shù)钠胶鈶B(tài)分布函數(shù)。由統(tǒng)計(jì)力學(xué)可知，在局域平衡條件下，分布函數(shù)可表示為局域平衡量的函數(shù)，考慮到ea的對(duì)稱性，將平衡態(tài)分布函數(shù)表示為

式中 a0、a1、a2——待定系數(shù)。

1．2 格子Boltzmann 方程

對(duì)式(3)的左側(cè)項(xiàng)進(jìn)行空間和演化時(shí)間的Taylor 級(jí)數(shù)展開(kāi)，取其到二級(jí)項(xiàng)可得

應(yīng)用Chapman－Enskog 方法對(duì)時(shí)間系數(shù)、空間系數(shù)和分布函數(shù)進(jìn)行多尺度展開(kāi)，可得到

式中 ε——小參量。

將式(6)代入式(5)，在ε 一階小量下有

在ε 二階小量下有

將式(7)和式(8)分別對(duì)下標(biāo)求和得出

由此得到的式(9)、式(10)即為時(shí)間t1和t2尺度下的密度平衡方程。

1．3 宏觀熱擴(kuò)散方程

將式(4)代入式(2)，可以得到

為了在t1尺度上恢復(fù)方程的對(duì)流項(xiàng)，令

由式(9)、式(12)得到

由式(7)和式(14)得到將式(15)代入式(10)，得到

將式(14)乘以ε，式(16)乘以ε2，兩式相加便可回到宏觀尺度下的方程(1)。令，則

由式(11)、式(13)和式(17)可得出待定系數(shù)

2 數(shù)值實(shí)驗(yàn)

為驗(yàn)證得到的熱擴(kuò)散方程的Boltzmann 模型，采用文獻(xiàn)［6，8］中的算例，進(jìn)行了驗(yàn)證性模擬實(shí)驗(yàn)。

2．1 驗(yàn)證算例1

算例1 中熱擴(kuò)散方程的邊界條件為第一類邊界條件且為恒溫，具體的控制方程

算例1 的解析解在文獻(xiàn)［6］中已給出

算例1 中熱擴(kuò)散方程的格子Boltzmann 模型求解參數(shù):τ =1． 2，u =0． 2，v =0． 5，時(shí)間步長(zhǎng)Δt =0.01，空間步長(zhǎng)Δx=0．1，粒子的遷移速率為c=10，總長(zhǎng)度l=40，t 分別為280、290、300 和310。

圖2 給出了不同時(shí)刻t 下算例1 的模擬值與解析解。從圖中可以看出，t 為280、290、300 和310時(shí)，模擬值與解析解的最大絕對(duì)誤差分別為1．31E－04、8．70E－05、7．20E －05 和7．00E －06，從而說(shuō)明數(shù)值解與解析解吻合良好。

網(wǎng)格數(shù)目對(duì)模擬結(jié)果精度有一定的影響。本文分析網(wǎng)格數(shù)目為80、200 和400 時(shí)，對(duì)計(jì)算精度的影響情況。

圖2 不同時(shí)刻下算例1 的模擬值與解析解的比較

表1 為τ=1．2，u=0．2，v=0．5，時(shí)間步長(zhǎng)Δt=0．01，t =300，總長(zhǎng)度l =40，空間步長(zhǎng)Δx 分別為0.5，0．2 和0．1 時(shí)，各個(gè)節(jié)點(diǎn)在不同網(wǎng)格數(shù)目下的解析解與數(shù)值解及其絕對(duì)誤差。從表中可以看出，網(wǎng)格數(shù)目分別為80、200 和400 時(shí)的最大絕對(duì)誤差分別為3．02E－03、3．31E－04 和7．00E－05。由此可見(jiàn)，隨著網(wǎng)格的細(xì)化，模擬結(jié)果的絕對(duì)誤差越來(lái)越小。

2．2 驗(yàn)證算例2

算例2 中的熱方程的邊界條件是隨著時(shí)間t 變化的，具體方程如下該問(wèn)題的解析解在文獻(xiàn)［8］中已給出

表1 各節(jié)點(diǎn)在不同網(wǎng)格數(shù)目下的解析解與模擬值及其絕對(duì)誤差

圖3 給出了不同時(shí)刻t 下算例2 的模擬值與解析解。從圖中可以看出，t 分別為1、4、7 和10 時(shí)的最大絕對(duì)誤差分別為6．21E－04、4．62E－04、4．98E－04 和5．09E－04，從而說(shuō)明數(shù)值解與解析解吻合良好。

圖3 不同時(shí)刻下算例2 的模擬值與解析解的比較

3 結(jié)語(yǔ)

本文基于格子Boltzmann 方法，采用Chapman－Enskog 展開(kāi)和多尺度技術(shù)，建立了求解一維熱擴(kuò)散方程的D1Q3 模型，模型的數(shù)值解與文獻(xiàn)的解析解吻合良好，其兩者的誤差隨網(wǎng)格細(xì)化而大幅度減小，從而說(shuō)明了本文模型可用于求解一維熱擴(kuò)散方程。

［1］Raabe D． Overview of the lattice Boltzmann method for nano－ and microscale fluid dynamics in materials science and engineering［J］． Modeling and Simulation in Materials Science Engineering，2004，12(6):R13 －R46．

［2］Tang G H，Tao W Q，He Y L． Gas slippage effect on microscale porous flow using the lattice Boltzmann method［J］．Physical Review E，2005，72(5):056301．

［3］Grunau D，Chen S，Eggert K． A lattice Boltzmann model for multiphase fluid flows［J］． Physics of Fluids，1993，5(10):2557 －2562．

［4］Miller W，Succi S，Mansutti D． Lattice Boltzmann model for anisotropic liquid－solid phase transition［J］．Physical review letters，2001，86(16):3578 －3581．

［5］鄧敏藝，劉慕仁，李玨，等． 一維有源擴(kuò)散方程的格子Boltzmann 方法［J］． 廣西師范大學(xué)學(xué)報(bào):自然科學(xué)版，2000，18(1):9 －12．

［6］劉慕仁，何云，陳若航，等． 一維對(duì)流擴(kuò)散方程的格子Boltzmann 方法［J］．廣西科學(xué)，1999，6(3):168 －169．

［7］徐世英，衛(wèi)玉敏，吳春光，等．一維Tyson 反應(yīng)擴(kuò)散系統(tǒng)的格子Boltzmann 方法模擬［J］．計(jì)算機(jī)與應(yīng)用化學(xué)．2008，25(5):537 －540．

［8］李貴艷，羅東升． 對(duì)流擴(kuò)散方程的數(shù)值計(jì)算［J］． 數(shù)學(xué)雜志，2009，29(2):186 －190．