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