【文章內(nèi)容簡介】 、 Lallemand 和 Frisch 提出了四維面心立方模型，即 FCHC 模型。這兩種模型都成功克服了 HPP 模型對(duì)稱性不足的缺點(diǎn)，在一定條件下可以得到二維和三維的不可壓縮 NavierStokes 方程。 自格子氣元胞自動(dòng)機(jī)模型 (Lattice Gas Automata， LGA)出現(xiàn)以來，得到了一定的發(fā)展和進(jìn)步。 LGA 模型將流體視為駐留在一個(gè)規(guī)則格子上的大量假想粒子，這些例子按照一定的規(guī)格在格子上進(jìn)行碰撞和遷移，通過統(tǒng)計(jì)學(xué)的知識(shí)得到宏觀平均值。雖然 LAG方法物理清晰直觀，計(jì)算較為穩(wěn)定，但是仍然存在一些問題需要解決。例如結(jié)果中會(huì)含有統(tǒng)計(jì)噪聲，需要對(duì)時(shí)間和空間進(jìn)行數(shù)值平均，從而在很大程度上增加了計(jì)算量。同時(shí)推導(dǎo)出的動(dòng)量方程不滿足 Galiean 恒定性，并且算子具有指數(shù)復(fù)雜性，對(duì)存儲(chǔ)空間也有較為 嚴(yán)格的要求。為了解決這些問題， LAG 方法也在不斷被修正和改進(jìn)。 1988 年，McNamara 和 Zati 提出在 LGA中使用布爾變量的統(tǒng)計(jì)平均量粒子分布函數(shù)進(jìn)行演化以消除統(tǒng)計(jì)噪聲。 但是， LGA 方法的指數(shù)復(fù)雜性仍未得到解決。 隨后為了進(jìn)一步 對(duì) LGA方法進(jìn)行優(yōu)化 ， Higuere 和 Jinenez 提出了線性化碰撞算子模型， 近似代替計(jì)算量很大的指數(shù)型算子 。 在這種模型中，引入 分布函數(shù)距離其平衡態(tài)分布函數(shù) ()eqif ，即 格子 Boltzmann 方法及其在非牛頓流體力學(xué)的應(yīng)用 4 ( ) ( )eq neqi i if f f??。其中 ()eqif 是 FermiDirac 分布的展開形式， ()neqif 為非平衡部分。在低速條件下，可以進(jìn)一步將 ()eqif 按照速度項(xiàng)進(jìn)行展開： ? ?( ) ( ) , 0 ( ) ,1 ( ) , 2 3e q e q e q e qi i i if f f f O M a? ? ? ? () 其中 ( ),eqkif 為包含速度 uk 的項(xiàng)。 將碰撞算子在全距平衡態(tài) ( ),0eqif 處展開， ? ? ? ? ? ? ? ? ? ?0 ( ) , 0 ( ) , 0 ( ) , 012e q e q e qi i i j j i i j k j i k iΩ f Ω Ω ff Ω f f f f? ? ? ? ? ? () 略去高階小量可得簡便形式如下： ? ? ? ? ? ? ? ?( ) , 0 ( ) , 0 ( ) , 0 ( )12e q e q e q n e qi i j j i i j k j i k i i j iΩ f Ω ff Ω f f f f Ω f? ? ? ? ? ? () 取 ()eqiiff? ，同時(shí)略去高階小量，可得 ( ) ,1 ( ) ,1 ( ) ,11 02e q e q e qi j i i j k j kΩ f Ω ff?? () 于是得到了一個(gè)近似線性化的碰撞算子，其中 K 為線性化碰撞矩陣。由于采用了線性化手段， HJ 模型大大提高了計(jì)算效率，為 LBE 模型的發(fā)展奠定了重要基礎(chǔ)。 19911992 年，一些專家學(xué)者 又 相繼提出了單松弛模型 (SRT 模型 )和 BGK 模型。在格子 BGK 模型 (即 LBGK 模型 )中，碰撞過程用區(qū)域某一平衡態(tài)的松弛過程代替，而碰撞矩陣由松弛時(shí)間決定。這一模型完全克服了格子氣自動(dòng)機(jī)的一系列缺點(diǎn)，因此也是目前最基本且應(yīng)用最廣泛的計(jì)算模型 [11]。 格子 Boltzmann方法的基本結(jié)構(gòu) 格子 Boltzmann 方程模型包含三個(gè)要素： (1) 流體粒子的離散速度集合， (2) 格子結(jié)構(gòu)， (3) 演化方程。它描述了具有離散速度的流體粒子分布函數(shù)在一個(gè)固定格子上的演化運(yùn)動(dòng)過程。 ? ?( , ) ( , ) ,i i t t i if x c t f x t Ω xt??? ? ? ? () 其中 x 是格子上的一個(gè)格點(diǎn)， ci (i=1， 2， 3， 4， … b)是流體粒子的離散速度集合，δt 為離散時(shí)間步長， t 為當(dāng)前時(shí)間步； fi 是以速度 ci 運(yùn)動(dòng)的速度分布函數(shù)， Ωi 是碰撞算子，表示分子間的碰撞對(duì)速度分布函數(shù)的影響。流體的密度 ρ、速度 u 和內(nèi)能 e 可由離散分布函數(shù)的速度矩得到。 格子 Boltzmann 方法及其在非牛頓流體力學(xué)的應(yīng)用 5 ? ? 2122iiiiiiiifu c fD R Te c e f??????? ? ???? () 其中 D 為空間維數(shù)。 演化方程又稱為格子 Boltzmann方程，其中的碰撞算子反映了微觀流體的相互作用。因此碰撞算子對(duì)模型能否準(zhǔn)確刻畫出流體系統(tǒng)的物理規(guī)律起到了決定性的作用。目前LBE 方法一般采用線性化碰撞算子，即 ()eqi ij j ijΩ K f f??????? () 其中 K 是一個(gè) bb 的碰撞矩陣， ()eqif 是依賴于宏觀物理量的平衡態(tài)分布函數(shù)。因此，碰撞矩陣和平衡態(tài)分布函數(shù)完全決定了碰撞算子的特性，也就相應(yīng)決定了一個(gè) LBE 模型刻畫流動(dòng)問題的精度。 從上述分析中可以發(fā)現(xiàn)，格子 Boltzmann 方法是一種不同于傳統(tǒng)數(shù)值方法的流體計(jì)算和建模方法。從離散的網(wǎng)格來說， LBE 具有 Euler 方法的屬性，而從離散的粒子觀點(diǎn)來說， LBE 又具有 Lagrange 方法的特點(diǎn)。這種特性使得 LBE 方法具有以下幾個(gè)顯著特點(diǎn)： (1) 雖然 LBE模型主要用于模擬宏觀連續(xù)流動(dòng)，但它不是基于連續(xù)模型的 NS方程，而是基于介觀模型，本身沒有連續(xù)介質(zhì)條件的假設(shè)。因此，只要設(shè)計(jì)的參數(shù)合理，可以用于描述微觀尺度和稀薄流等非連續(xù)流動(dòng)問題。 (2) 格子 Boltzmann 方法可以較為直觀地處理流體內(nèi)部 以及流體和周圍環(huán)境的相互作用，從而給多組分、多相態(tài)系統(tǒng)、界面動(dòng)力學(xué)、滲流等復(fù)雜流動(dòng)現(xiàn)象的研究提供了方便； (3) 從計(jì)算的角度來看， LBE 方法計(jì)算簡單，易于用編程實(shí)現(xiàn)，具有較好的并行性和可擴(kuò)展性，對(duì)于大規(guī)模流動(dòng)問題的計(jì)算具有很大優(yōu)勢。 本文的主要工作 首先是通過對(duì) LBM 理論方法的學(xué)習(xí)，掌握 LBM 方法的整個(gè)計(jì)算流程，以及整個(gè)方法計(jì)算過程中各部分之間的聯(lián)系，包括粒子的碰撞遷移過程， 和 如何選取邊界條件，不同的邊界格式的適用條件對(duì)于程序計(jì)算 過程 穩(wěn)定性的影響 。 其次是非牛頓流體 LBM模型的構(gòu)建，包括冪律流體和 Bingham 流體 。 在 LBGK 模型的基礎(chǔ)上，結(jié)合非牛頓流體的本構(gòu)方程，建立了本文的非牛頓流體計(jì)算模型 。由于 Bingham 流體本構(gòu)方程本身的不連續(xù)性，對(duì)模擬造成了一定的困難，所以模擬中選取了 Papanastasiou 模型，通過模擬直通道中 冪律流體 的流動(dòng)，并且與解析解進(jìn)行了比較，數(shù)值解與解析解吻合良好。最后 格子 Boltzmann 方法及其在非牛頓流體力學(xué)的應(yīng)用 6 對(duì) Bingham 流體在 2:1 突擴(kuò)通道中的流動(dòng)過程進(jìn)行了模擬，固定 Bn 數(shù)，選取不同的 Re數(shù)來考察屈服區(qū)域和非屈服區(qū)域的分布情況以及流場情況，最后通過無量綱漩渦長度定量分析了模擬結(jié)果的正確性。 格子 Boltzmann 方法及其在非牛頓流體力學(xué)的應(yīng)用 7 2 LBGK模型及邊界處理 LBGK模型介紹 現(xiàn)在已經(jīng)發(fā)展完善的格子 Boltzmann模型中，一般包含了三個(gè)部分 [11]：離散速度型，平衡態(tài)分布函數(shù)，分布函數(shù)的演化方程。這三大部分是彼此相關(guān)相互聯(lián)系的，構(gòu)造格子Boltzmann 模型的關(guān)鍵是選擇合適的平衡態(tài)分布函數(shù)，而平衡態(tài)分布函數(shù)的具體形式又與構(gòu)造離散速度模型有關(guān)，離散速度的對(duì)稱性決定了相應(yīng)的格子 Boltzmann 模型能否還原到所要求解的宏觀方程，因此，選擇建立合理模型對(duì)準(zhǔn)確還原到宏觀方程至關(guān)重要。 單松弛模型或 LBGK 模型是至今應(yīng)用最為廣泛的 LBM 模型，由 Qian， Chen 等人于 1992 年分別獨(dú)立地提出的，是用一種簡單的松弛方法來簡化復(fù)雜的碰撞過程。其中最具代表性的是由他們提出 DnQb 模型 [12](n代表空間維數(shù)， b 代表離散速度數(shù) )，被視為格子 Boltzmann方法的基本模型，其他常用模型都是在此模型基礎(chǔ)上發(fā)展演化而形成的，下面詳細(xì)介紹一下格子 Boltzmann 方法常用的模型。 格子 BoltzmannLBGK方程 根據(jù)微觀氣體動(dòng)力學(xué)原理和統(tǒng)計(jì)學(xué)方法，可以描述微觀運(yùn)動(dòng)的格子方程，流體特性由粒子分布函數(shù) ),( tfi x 描述，該函數(shù)給出了虛擬流體質(zhì)點(diǎn)在離散速度 ie 上， t 時(shí)刻位于網(wǎng)格點(diǎn) x 處的粒子密度分布函數(shù)。格點(diǎn)上碰撞和遷移所引起的平衡態(tài)分布函數(shù)的演化過程是由離散的 Boltzmann 方程所決定的，如下所示 ),(),(),( ttftttf iii xxex ??????? i (i=0， 1， ...， 8) () 式中， t? 是時(shí)間步長，這里引入了一個(gè)碰撞算子 i? 是關(guān)于 ),( tfi x 的函數(shù)，表示的是碰撞造成粒子分布函數(shù)的變化所占比率，定量描述了分子間的微觀相互作用。在LBGK 模型中利用了一個(gè)簡單的碰撞算子 i? 代替 LBE 中原有碰撞項(xiàng)： )],(),([1),( )( tftft eqiii xxx ???? ? (i=0， 1， ...， 8) () 式中， ),( tfi x 為速度平衡態(tài)分布方程， ? 為無量綱的弛豫時(shí)間，只與流體的物性參數(shù) (如：粘性系數(shù)、熱傳導(dǎo)系數(shù)、流體的質(zhì)量擴(kuò)散系數(shù)等 )有關(guān)。 將 LBE 方程進(jìn)行離散，粒子的速度在相空間離散為有限維的速度空間，由此可以得到含有外力項(xiàng)的格子 Boltzmann 方程的有限差分格式： )],(),([1),(),( )( tftftftttf eqiiii xxxex ???????? ?i () 格子 Boltzmann 方法及其在非牛頓流體力學(xué)的應(yīng)用 8 此方程極大簡化了原來的 LBE 方程，而且適用于不同平衡態(tài)分布函數(shù)和多物理場問題，方程具有 Lagrange 特性，一階精度的 LBGK 方程實(shí)際上具有二階精度，令該模型的計(jì)算精度基本達(dá)到二階精度，使很多很復(fù)雜的邊界條件都能夠很容易得到解決，因此該模型具有較強(qiáng)的邊界問題處理能力。 DnQb離散模型 在 DnQb 系列模型中， n 代表空間維數(shù)， b 代表離散速度數(shù)，根據(jù) n 和 b 的不同，可以分為 D2Q D2Q D3Q1 D3Q19 模型等，這些模型中所采用的建模控制方程和平衡態(tài)分布函數(shù)的求解方程形式均相同，區(qū)別在于離散速度的具體配置情況。本文所采用的是 D2Q9 模型。 在此模型中， 我們將不可壓縮流體流動(dòng)和一個(gè)九速度模型表示在一個(gè)二維格子(D2Q9)上，網(wǎng)格采用標(biāo)準(zhǔn)的正方形格子，根據(jù)速度大小和方向，其離散速度配置如圖 所示。其中 D 代表空間維度，而 Q 則代表在一個(gè)節(jié)點(diǎn)的不同速度數(shù)量。每個(gè)格點(diǎn)只與相鄰的 8 個(gè)格點(diǎn)接觸， LBGK 方法的思想是，用粒子分布函數(shù)來代表各個(gè)格點(diǎn)的流體，在流體流動(dòng)過程中，在該區(qū)域內(nèi)的流體離散分向 9 個(gè) 方向，當(dāng) i=0 時(shí)， fi 停留在原地，而當(dāng) i=1~8 時(shí)，則粒子將沿如圖 所示的 8 個(gè)方向以速度 ei 運(yùn)動(dòng)，經(jīng)過時(shí)間步長 t? 到達(dá)相應(yīng)格點(diǎn)。經(jīng)歷了一個(gè) t? 后，在內(nèi)部相互的作用下，各格點(diǎn)所代表的粒子又匯聚成了新的流體群，上述流動(dòng)過程將會(huì)隨著時(shí)間繼續(xù)不斷發(fā)生。 圖 D2Q9 離散速度模型 速度配置如下： 格子 Boltzmann 方法及其在非牛頓流體力學(xué)的應(yīng)用 9 ( 0 , 0) 0c os ( 1 ) , sin ( 1 ) 1 4222 c os ( 2 1 ) , sin ( 2 1 ) 5 844iie c i i ic i i i????????? ??? ? ? ? ?? ? ? ? ?? ??? ? ? ?? ? ? ????? ??? ? ? ?? ? ? ? ???? ? ? ?? ? ? ?? ??? () 其中， txc ?? /? 為格子速度， x? 和 t? 分別為網(wǎng)格步長和時(shí)間步長，在通常情況下 x 方向和 y 方向的網(wǎng)格步長相同，即 x? = t? ，在模擬過程中兩者的取值均為 1。 在 DnQb 系列模型中，平衡態(tài)分布函數(shù) 有一個(gè)統(tǒng)一的表達(dá)形式： ]22 )(1[ 24 22)(sssieqi cccf 2ii uueue ?????? ?? () 其中， i? 是代表粒子速度函數(shù)權(quán)系數(shù)， sc 與聲速有關(guān) (或稱為格子聲速 )，這兩個(gè)參數(shù)是決定 LBGK 模型的關(guān)鍵參數(shù)。依賴于所使用的格子類型，對(duì)于 D2Q9 模型相關(guān)參數(shù)設(shè)置如下： 4 9 01 9 1 41 3 6 5 8iiii????? ? ??? ??? () 3ccs ? () 此模型中的宏觀壓力與密度之間關(guān)系如下 ：