【正文】 通過改變冪律指數(shù) n 來考察不同冪律流體的流動特性 (n1 為假塑性流體。格子 Boltzmann 算法 (LBM)流程圖如圖 所示，值得注意的是，非牛頓流體的松弛時間 ? 會隨著運動粘度的變化而變化，因此在每個節(jié)點都需要進行松弛時間等的計算。在 Dαβ的計算過程中，Dαβ估算基于之前所得的物理量，碰撞項中的 ? 是從表觀粘度 ),( tx? 計算得到。尤其，消除計算速度導(dǎo)數(shù)的需要成為一個顯著的優(yōu)點在計算復(fù)雜流動通過流道的固體障礙時。這樣，我們可以從粒子分布函數(shù)的非平衡部分直接計算局部剪切速率在二階精度上。下標 ? 和 ? 表示笛卡爾坐標。結(jié)合冪律流體，可以看出牛頓流體只是非牛頓流體的一種特殊形式。對于足夠大的 m 值，兩者近似，式 ()可以用來代替賓漢姆流體的本構(gòu)方程。其中 是 D1的不變量，通常定義 為 )(2 21Dtr??? () 實際上不能精確的確定 ，它的值可能依賴于用于確定它的儀器。當切應(yīng)力大于屈服應(yīng)力之后，流體開始流動，其剪應(yīng)力與剪切速率的關(guān)系與牛頓流體類似。上式稱為賓漢方程，符合賓漢方程的流體稱為賓漢流體，塑性流體也稱為賓漢流體。 賓漢流體的表觀粘度為： () 可以看出，賓漢流體的表觀粘度是隨流速梯度而變化的。當 n1時，流體表觀粘度隨剪切變形速率的增大而減小，稱這種流體為假塑性流體，也稱為剪切稀化流。當 n=1 時，流體即為牛頓流體。工程應(yīng)用中常將 ()式改寫成類似牛頓流體的本構(gòu)方程 ??? ?0? () 式中 ?? ? 10 ?? nK () 稱為冪律流體的表觀粘度。這類流體的本構(gòu)方程是 nK?? ?? () 式中， ? 為剪切應(yīng)力； K 為稠度系數(shù)，或稱冪律系數(shù)，是粘度的度量，其隨著粘度的增大而增大； n 代表流性指數(shù)，又稱冪律指數(shù)，反映偏離牛頓流體的程度。 冪律流體與賓漢姆流體本構(gòu)方程 許多種非牛頓流體的實驗表明，剪切應(yīng)力雖然與剪切變形速率不呈線性函數(shù)關(guān)系，但是可以表示為簡單的冪函數(shù)關(guān)系。 在 LBM 中，對于非牛頓流體，由于其運動粘度與剪切速率有關(guān)，所以每個格子節(jié)點的粘度都要根據(jù)與剪切速率關(guān)系進行計算，計算出粘度之后，才能得出各個節(jié)點的松弛時間。 在外力作用下，固體一般先產(chǎn)生彈性變形，當應(yīng)力超過屈服極限時就產(chǎn)生塑性變形，不是有任意微小剪切力作用，就會產(chǎn)生流動變形，它有一個屈服限，稱為屈服應(yīng)力。 非牛頓流體的彈塑性 非牛頓流體一般為液體，屬于不可壓縮流體。這里不同于牛頓流體，其中 a? 不是常數(shù)，而是剪切速率的函數(shù)。 非牛頓流體的物質(zhì)構(gòu)成要比牛頓流體的要復(fù)雜得多。 非牛頓流體的粘性 流體的分子結(jié)構(gòu)決定流體的粘性。非牛頓流體的種類很多，如血液、蛋白液、油漆、泥漿以及高分子聚合物的溶液等。物體受外力作用，當外力除去后能恢復(fù)原狀的變形，則稱為彈性變形，不能恢復(fù)原狀的變形稱為塑性變形。若流場屬于定常問題，則需要判斷相鄰時間步長的各物理量是 否達到收斂條件，若滿足收斂判據(jù)，則輸出結(jié)果，否則，繼續(xù)進行循環(huán)。至此主要計算完成。常用的邊界條件有周期邊界、反彈邊界、速度邊界、外推邊界等。和只進行遷移過程的邊界節(jié)點不同，內(nèi)部節(jié)點同時進行碰撞和遷移過程。在計算開始之前，需要選擇適合的網(wǎng)格類型。 格子 Boltzmann 方法的計算步驟 應(yīng)用 LB 方法模擬具體流場的算法程序流程圖如圖所示。uν i 2)(21 ti if ??? ???? ? () 在演化方程中增加作用力項 格子 Boltzmann 方法及其在非牛頓流體力學的應(yīng)用 17 這種作用力模型可用下式表示： ? ? ? ? ? ? ? ?? ? iteqiiittii Ftftftftf ???? ?????? ,1, )( xxxex () 其中， Fα表示作用在流體粒子上的總作用力 F 的離散形式，它不能直接對連續(xù)的作用力項進行離散處理。39。假設(shè)發(fā)生一次碰撞的平均時間是 τδt，則由作用力引起的液體動量變化為 Fτδt。因此，可通過適當修改平衡態(tài)粒子分布函數(shù)將狀態(tài)方程修改為 ??? )()( 2 ?? scp () 變換過程如下所示： ??????????????????0],2)()1(1[0],22 )()([22020224222)(icucpicucccpfssssssieqi????????? ueue ii () 這樣就得到了包含作用力 F 的 LBGK 模型。標準的 LBGK 方程對應(yīng)的宏觀流動方程如下： 格子 Boltzmann 方法及其在非牛頓流體力學的應(yīng)用 16 0)( ???? u????t () )2()()( Su ?????? ????????? puut () 其中 ,2?scp? ,)21(2 tsc ??? ?? T)( uuS ???? () 將作用力表示為勢函數(shù)的形式，即 ??? ???? aF 。 平衡態(tài)分布的壓力校正方法 應(yīng)用平衡態(tài)分布的壓力校正方法必須滿足的條件上，作用力可以表示成勢函數(shù)的形式，同時作用力導(dǎo)致的密度變化可以忽略不計。因此，如何在離散速度的框架下描述作用力是一個十分重要的問題。在 D2Q9 模型中詳細介紹該方法，如圖 ，假設(shè) 301 是邊界， 625 位于流場內(nèi)， 748 位于流場外， t 是當前時刻，已知 2 點粒子分布函數(shù)，求 0 點粒子分布函數(shù)，該方法將 0 點粒子分布函數(shù) ),0( tfi分為兩部分： ),0(),0(),0( tftftf n e qieqii ?? () 節(jié)點 2 粒子分布函數(shù) ),2( tfi ，宏觀速度 ),2( tu 及密度 ),2( t? 都已知，因此 2 點非平衡態(tài)粒子分布函數(shù)可以計算為 ),2(),2(),2( tftftf eqiin e qi ?? () 0 點平衡態(tài)部分近似為 )),0(),2((),0( ttftf eqieqi u?? () 通過非平衡外推可得 0 點粒子分布函數(shù) )],2(),2([),0(),0( tftftftf eqiieqii ??? () 非穩(wěn)態(tài)外推格式在時間和空間上都能夠達到二階精度，同時具有很好的數(shù)值穩(wěn)定性，并且兼顧了其他外推格式的諸多優(yōu)點，且其計算不復(fù)雜，較容易實現(xiàn)，因此適用范圍十分廣泛。 (4) Chen 格式 Chen 格式在時間和空間上都具有二階精度，很好 的保持了格子 Boltzmann 方法的整體精度，無需對邊界的流體分布函數(shù)作任何假設(shè)，因此適用非常廣泛，但也有其缺點，就是數(shù)值穩(wěn)定性較差。 Nobel 處理格式：通過構(gòu)造壓力約束來求解邊界上未知粒子的分布函數(shù)，其基本思想是：在邊界節(jié)點速度已知的情況下，根據(jù)邊界上分布函數(shù)與宏觀量之間的關(guān)系，確定邊界節(jié)點的分布函數(shù)， 不管是速度邊界還是壓力邊界，實際上就是對流體施加一個外力作用，引起動量變化， 非平衡反彈格式： 1997 年 Zou 和 He 對二維和三維的 LBGK 模型進行了研究，提出來這種新的邊界處理格式，在實驗研究過程中通過與其他處理格式比較，發(fā)現(xiàn)這種格式具有良好的數(shù)值穩(wěn)定性。反彈格式在處理復(fù)雜的不規(guī)則界面具有很大優(yōu)勢。該格式易于操作計算簡單，且能滿足系統(tǒng)的各物理量守恒，但是精度不高，只能達到一階精度，且對運動邊界等復(fù)雜邊界很難處理，為了改進這一缺點，人們提出了修正反彈格式和半步長反彈格式：將物理壁面移至格點中間，在邊界與格點之間應(yīng)用修正反彈格式和半步長反 彈格式。包括標準反彈格式、半步長反彈格式及修正反彈格式。 充分發(fā)展邊界處理格式：針對的是流體在管內(nèi)流動達到充分發(fā)展段后，此時流體的速度和密度在主流方向上都不再發(fā)生改變，因而可以采用這種數(shù)值穩(wěn)定性較好的邊界處理格式。周期性邊界處理格式相比于其他常用的邊界處理格式，在數(shù)值計算的 穩(wěn)定收斂性方面保持最好。 周期性邊界處理格式：數(shù)值模擬時，如果流場在空間內(nèi)呈現(xiàn)周期性變化或在某一方向無限大時，則可以將周期性單元取出作為模擬區(qū)域，并在相應(yīng)的那些邊界上采用周期性邊界。本文主要應(yīng)用的是平直邊界格式，下面介紹幾種常用的邊界格式。因此，需要從已知的宏觀上的邊界條件來試圖確定邊界節(jié)點上的粒子分布函數(shù)，在這個過程所應(yīng)用的方法稱為格子 Boltzmann方 法的邊界處理方法，所設(shè)計出的計算格式稱為邊界處理格式。邊界條件處理的方式對計算的效率、計算數(shù)值的穩(wěn)定性和計算結(jié)果的精度都有很大影響。以穩(wěn)態(tài)問題為例，當系統(tǒng)達到穩(wěn)定后，流場和溫度場與初始條件無關(guān)，而主要取決于邊界條件。另一種方法是迭代 計算方法 [14]，其主旨是借助格子 Boltzmann 方程求解初始壓力對應(yīng)的 Poisson 方程，并未相應(yīng)的分布函數(shù)賦初值，由此獲得與速度場相匹配的初始分布函數(shù)。 (1) D1Q3 模型： ? ?1,1,0 ??ie 3,31,32 210 cc s ??? ??? () (2) D1Q5 模型： ? ?2,2,1,1,0 ???ie cc s ???? ?? ,61,61,32 43210 ??? () (3) D2Q7 模型： ?????? ?????? 23 2123 210123 2123 210100ie 2,121,32 610 cc s ??? ??? () (4) D2Q9 模型： ?????? ??????? 111111111001100100ie 格子 Boltzmann 方法及其在非牛頓流體力學的應(yīng)用 11 3,361,91,94 85410 cc s ???? ?? ??? () (5) D3Q15 模型： ??????????????????????????111111111111111111111111100100010010001001000ie 3,721,91,92 147610 cc s ???? ?? ??? () (6) D3Q19 模型： ??????????????????????????110110110110101101101101011011011011100100010010001001000ie 3/,36/1,18/1,3/1 187610 cc s ???? ?? ??? () 初始條件及邊界條件處理 初始條件和邊界條件是流體力學研究中的重要課題，一般說來，初始條件和邊界條件總是按著宏觀量 (壓力、速度、溫度 )給出，但在 LBM 中基本變量是分布函數(shù)，雖然根據(jù)分布函數(shù)可以簡單的確定宏觀流動變量，但是根據(jù)宏觀量確定分布函數(shù)卻不簡單，因此如何根據(jù)給定的初始和邊界條件構(gòu)造離散分布函數(shù)的相應(yīng)條件，是 LBM 用于實際計算的前提，大量實驗研究表明，實現(xiàn)初始和邊界條件的不同方法對 LBM 的計算精度、數(shù)值穩(wěn)定性和計算效率都有很大的影響。換言之， LBGK 方法實際上是求解不可壓縮流體 NS 方程的一種人工壓縮法。依賴于所使用的格子類型，對于 D2Q9 模型相關(guān)參數(shù)設(shè)置如下： 4 9 01 9 1 41 3 6 5 8iiii????? ? ??? ??? () 3ccs ? () 此模型中的宏觀壓力與密度之間關(guān)系如下 ： 2),( Sctp x?? () 模型中宏觀密度 ),( tx? 和和速度 ),tu(x 可以分別由 ),( tfi x 定義： ),(),( 80 tft i i xx ???? () ),(),( 1), 8 0 tftt ii xexu (x i??? ? () 使用 ChapmanEnskog 方法展開便可推導(dǎo)出該模型所對