顆粒流體系統(tǒng)的格子 Boltzmann 數(shù)值方法研究進展

王利民， 付少童

(1.中國科學院過程工程研究所 多相復雜系統(tǒng)國家重點實驗室，北京 100190；2.中國科學院大學 化學工程學院，北京 100049)

1 引 言

顆粒流體系統(tǒng)廣泛存在于自然界和工業(yè)過程，由于其涉及復雜的多相流動過程和時空多尺度結(jié)構(gòu)，通過實驗或理論分析存在較大困難，因此基于數(shù)值模擬展開顆粒流體系統(tǒng)的研究已成為不可或缺的重要手段。根據(jù)對流體相和顆粒相描述精度的不同，顆粒流體系統(tǒng)的模擬方法通常分為顆粒解析直接數(shù)值模擬PR-DNS(particle-resolved direct numerical simulation)、離散顆粒模擬DPS(discrete particle simulation)和雙流體模型TFM(two -fluid model)三類[1]。

顆粒流體系統(tǒng)模擬中,TFM的顆粒相處理為連續(xù)介質(zhì)的擬流體，PR-DNS和DPS的顆粒相通常是離散處理，根據(jù)牛頓第二定律對顆粒進行 Lagrange 跟蹤。目前在顆粒流體系統(tǒng)的PR-DNS中，常忽略或簡單處理顆粒之間的作用。然而，大多數(shù)體系中顆粒間的相互作用至關(guān)重要，需要準確和高效地對其進行描述。顆粒碰撞處理主要有軟球模型、硬球模型和直接Monte Carlo方法DSMC(direct simulation Monte Carlo)。軟球模型的典型代表是離散單元法DEM(discrete element method)[2]，該方法基于彈性阻尼和滑移等機制計算顆粒間的碰撞力，再根據(jù)牛頓第二定律采用時間驅(qū)動算法更新顆粒的速度和位置；硬球模型將顆粒碰撞處理成二體碰撞，其不求解顆粒受力，基于動量守恒原理直接獲得顆粒碰撞后的速度信息；DSMC[3]同樣采用硬球模型更新碰撞后顆粒速度與角速度，但運用概率抽樣確定顆粒是否碰撞，因此該方法也可視為硬球模型的一種變化形式。其中，軟球模型由于能夠獲得顆粒詳盡的受力信息而得到廣泛應(yīng)用。

顆粒流體系統(tǒng)中的流體相主要基于連續(xù)介質(zhì)假設(shè)框架下的N-S(Navier-Stokes)方程或其簡化形式，采用常規(guī)的有限差分、有限體積和有限元等方法進行求解。由于該類方法流程復雜且包含隱式算法，在并行化和大規(guī)模計算上還存在困難[4]。格子Boltzmann方法LBM(lattice Boltzmann method)[5]是20世紀80年代末逐漸發(fā)展起來的一種流體系統(tǒng)建模和模擬的新方法。從離散的網(wǎng)格來說，其具有歐拉方法的屬性；從離散的粒子觀點來看，又具備拉格朗日方法的特點。這種特性使得該方法的研究思路與傳統(tǒng)的流體模擬方法完全不同，其屬于介于宏觀連續(xù)介質(zhì)模型和微觀分子動力學模型之間的介觀模型，且物理背景清晰，受到國內(nèi)外學者的廣泛關(guān)注。隨著LBM基礎(chǔ)理論不斷完善，目前其研究已經(jīng)從理論向工程應(yīng)用發(fā)展，并在許多復雜物理現(xiàn)象和過程進行了大量的應(yīng)用[6]，如湍流DNS、大渦模擬和耦合各種湍流模型的模擬等；針對多組分系統(tǒng)和不混溶的多流體系統(tǒng)，分別建立了顏色模型、偽勢模型、自由能模型和動力學模型等；針對化學反應(yīng)系統(tǒng)，建立擴散和燃燒模型。此外，LBM已拓展到模擬非牛頓流體、磁流體、聚合物流動、血液流、光子運動、聲場和微流動等復雜體系[6]。

工業(yè)流化床反應(yīng)器的最大流動結(jié)構(gòu)可達到米量級，但這些結(jié)構(gòu)卻受發(fā)生在毫米量級的顆粒-顆粒間碰撞和顆粒-流體間相互作用的直接影響。目前很難通過某一具體的模擬方法涵蓋所有的時間和空間尺度，根據(jù)關(guān)注重點問題的不同，應(yīng)用格子Boltzmann框架可在不同時空尺度和精度上開展顆粒流體系統(tǒng)模擬(圖1)，如在微觀層次，LB -based PR-DNS擬獲得的詳細信息深入探索顆粒流體界面的流動、傳遞和反應(yīng)以及兩相相互作用的本構(gòu)關(guān)系；在介觀層次，利用 LB -based DPS探索實驗室規(guī)模的顆粒流體流動，建立更為準確的宏觀多相流動模型；在宏觀層次，發(fā)展格子Boltzmann雙流體模型(LB -based TFM)有望實現(xiàn)工業(yè)級別的運算。

圖1 基于格子Boltzmann的顆粒流體系統(tǒng)計算框架

本文將重點介紹LBM在顆粒流體系統(tǒng)模擬方面的研究進展，包括 LB -based PR-DNS方法、LB -based DPS方法及其最新應(yīng)用，以及關(guān)于開發(fā) LB -based TFM研究的一些進展。

2 顆粒解析直接數(shù)值模擬

顆粒流體系統(tǒng)PR-DNS與單相湍流DNS相類似，其中流體相的控制方程直接采用數(shù)值計算求解，無任何湍流模型，流場中顆粒周圍計算網(wǎng)格縮小到顆粒尺度以下，能夠分辨包含Kolmogorov耗散尺度以內(nèi)的所有空間尺度，而顆粒的受力則通過對其表面的黏性力與壓力進行積分獲得，不引入經(jīng)驗模型[7]。根據(jù)流體相的描述方式不同，顆粒流體系統(tǒng)的PR-DNS又可以分成三類(圖2)，一是基于N-S方程的PR-DNS方法；二是基于格子的PR-DNS方法；三是基于粒子的PR-DNS方法。LB -based PR-DNS方法屬于基于格子的PR-DNS方法。根據(jù)流固作用處理方式不同，基于LBM的顆粒流體系統(tǒng)PR-DNS分為邊界鏈法和格點法。

圖2 顆粒流體系統(tǒng)的顆粒解析直接數(shù)值模擬

2.1 邊界鏈法

Ladd[8]率先提出了基于有限體積顆粒的LBM模擬液固懸浮，通過修正的回彈邊界保證運動邊界的無滑移，通過動量交換計算流體和固體顆粒之間的作用力。Ladd方法中固體顆粒實際是由一個階梯狀的殼組成，固體顆粒邊界像許多流體-固體格子鏈組成的，因此稱為邊界鏈法。

Ladd方法[8]的流固作用通過流體粒子和固體顆粒之間的動量交換實現(xiàn)，其優(yōu)點是計算直觀簡潔，而且保證了流體的質(zhì)量和動量局部守恒。但這種階梯狀的殼近似表示真實固體顆粒，往往會引入誤差，而且在大多數(shù)情況下，這種誤差會導致顆粒的動力學半徑略大于其真實半徑。只有在固體顆粒半徑遠大于LBM網(wǎng)格時，固體顆粒的動力學半徑才和其真實半徑相近。此外，Ladd方法中，固體顆粒內(nèi)部由流體粒子占據(jù)，要求固體與流體密度比大于1。針對這一情況，Aidun等[9]提出不使用顆粒內(nèi)部流體的方法，能夠處理固體密度小于流體密度的問題。在Ladd方法的基礎(chǔ)上，有學者進一步提出了滿足局部質(zhì)量守恒的方法[10]，并模擬了多顆粒沉降問題。然而，以Ladd為基礎(chǔ)的方法均存在著邊界形狀與真實邊界并不相符的問題，為了解決階梯狀近似，不少學者提出了改進方法，大部分都是采用處理曲面邊界的方法,如曲面邊界插值[11]。而這些處理方法雖然提高了邊界處理的精確性，但在計算效率上卻付出了巨大的代價，通常處理顆粒數(shù)目較少。但基于邊界鏈法已可實現(xiàn)模擬含17845個顆粒的湍管流[12]。

2.2 格點法

Feng等[13]拋棄了固定顆粒邊界的想法，把顆粒的邊界當作可變形，結(jié)合浸入邊界法，提供了一種新的處理運動流固邊界的方法，這種方法避開了顆粒邊界無法精確描述的困難，但對顆粒表面剛性系數(shù)的選擇具有隨意性。Noble等[14]提出了浸入運動邊界法來處理流固耦合，通過在LBM演化方程中加入附加碰撞項，實現(xiàn)了對固體顆粒邊界相對準確的描述，同時沿用了LBM的計算體系，保留了原算法的優(yōu)點。Noble等在浸入運動邊界法中引入了格子控制體和格子固含率的概念。格子固含率為格子控制體內(nèi)顆粒體積占格子控制體體積的比重，固含率越高表示格子內(nèi)包含的顆粒組分越多，相應(yīng)的流體格子與顆粒的相互作用越強。將格子固含率作為重要參數(shù)進一步得出權(quán)重系數(shù),實現(xiàn)對顆粒邊界的光滑描述。權(quán)重系數(shù)存在多種構(gòu)建形式[15-17]，均在一定程度上改進了模擬的精度。

Cook等[18]耦合浸入運動邊界法和DEM模擬了顆粒相互作用占主導的顆粒流。馮云田等[19]在此基礎(chǔ)上，引入大渦模擬LES(large eddy simulation)，成功地對高雷諾數(shù)顆粒流體系統(tǒng)進行了模擬。王利民等[20]發(fā)展了一種基于粒子-格子耦合方法的PR-DNS方法，該方法中顆粒模型采用結(jié)合硬球模型和軟球模型各自優(yōu)點的時間驅(qū)動硬球模型TDHS(time -driven hard-sphere)[21]進行求解，該算法將顆粒/流體密度比從2.5提高到了1500，成功復現(xiàn)了散式流態(tài)化和聚式流態(tài)化的典型現(xiàn)象[20]。在此基礎(chǔ)上，周國峰等[22]成功模擬了二維和三維DKT(Drafting-Kissing-Tumbling)過程，驗證了該算法的有效性，此外，還實現(xiàn)了Lee-Edwards邊界條件與基于格點 LB -based PR-DNS的耦合[23]，克服了邊界鏈接法無法滿足伽利略不變性的問題，并成功預測了低雷諾數(shù)下顆粒懸浮體系的表觀粘度。熊勤鋼等[7]借鑒基于粒子-格子耦合方法的顆粒流體系統(tǒng)PR-DNS的思想，使用簡化DEM代替硬球算法，利用中國科學院過程工程研究所研制的千萬億次超級計算系統(tǒng) Mole -8.5，使用672個GPU進行了大規(guī)模并行計算，實現(xiàn)了全球最大規(guī)模顆粒流體系統(tǒng)PR-DNS，其中二維模擬包含1166400個顆粒，三維模擬包含129024個顆粒(圖3)。

圖3 顆粒懸浮系統(tǒng)大規(guī)模直接模擬粒子和流場分布[7]

周國峰等[24]分析了非均勻氣固系統(tǒng)中顆粒受到的曳力情況，對氣固流態(tài)化系統(tǒng)進行了三維大規(guī)模PR-DNS，成功復現(xiàn)了流動的非均勻結(jié)構(gòu)，揭示了沒有考慮非均勻結(jié)構(gòu)的Wen &Yu曳力公式預測偏高，且受力方向存在偏差。

劉曉雯等[25,26]在中國科學院過程工程研究所研制的千萬億次超級計算系統(tǒng) Mole -8.5E 上，使用192個Tesla K80 GPU運行17天，對4.83×109個流體網(wǎng)格和115200個固體顆粒的周期懸浮進行了直接數(shù)值模擬(圖4)，探索了介尺度結(jié)構(gòu)對顆粒統(tǒng)計屬性、氣固相間作用力以及傳遞特性的影響。

圖4 大規(guī)模氣固懸浮系統(tǒng)直接模擬瞬態(tài)結(jié)果[26]

劉曉雯等[25,26]發(fā)現(xiàn)顆粒脈動速度和局部平均無量綱曳力具備很強的尺度依賴性，提出了一種雙峰分布形式以修正顆粒脈動速度分布，闡明了顆粒脈動速度偏離Maxwell分布的內(nèi)在機理，并從微觀層次上揭示了介尺度曳力的尺度依賴性，提出了以Fr數(shù)反映結(jié)構(gòu)效應(yīng)的曳力關(guān)聯(lián)式[26]。

此外，耦合直接力浸入邊界法的LBM應(yīng)用于三維顆粒懸浮PR-DNS，如1024個熱顆粒在流化床中的運動[27]；Rubinstein等[28]模擬了2000個顆粒在流化床中的運動；Eshginejadfard等[29]對305個顆粒槽道流進行了模擬。LB -based PR-DNS也應(yīng)用于分析多顆粒沉降[30,31]或顆粒自旋及振動對流動的影響[32-34]、圓柱與柔性體耦合[35]、管道泄漏時土壤的流化現(xiàn)象[36]、泥石流[37]以及耦合傳熱的顆粒流問題[38]。通過對LBM中的弛豫時間進行修正，LB -based PR-DNS擴展到研究非牛頓流體顆粒流體系統(tǒng)中的顆粒流動規(guī)律[39-44]。通過重構(gòu)顆粒的形式和顆粒間作用力，基于浸入邊界法的LBM可以非常容易地處理諸如多邊形顆粒[45]、可變形顆粒[46]和粘性顆粒[47]的遷移等，甚至可以分析巖土工程中的水力壓裂問題[48]。以上均顯示了基于 LB -based PR-DNS的計算策略在模擬顆粒流體系統(tǒng)中的巨大優(yōu)勢和潛力。

3 離散顆粒模擬

離散顆粒模擬(或顆粒軌道模型)將流體相處理為連續(xù)介質(zhì)，固相處理為離散顆粒。在歐拉坐標系考察流體運動，在拉格朗日坐標系以及顆粒層次分析固相運動[49]。由流固耦合的實現(xiàn)方式不同，離散顆粒模擬一般可分三類，(1)單向耦合(one-way coupling)，只考慮流體對顆粒的作用，不考慮顆粒對流體的影響；(2)雙向耦合(two-way coupling)，考慮流體和顆粒之間的相互作用；(3)四向耦合(four-way coupling)，除了考慮流體和顆粒之間的相互作用外，還考慮顆粒間碰撞?，F(xiàn)在通常指的離散顆粒模擬，都是四向耦合的確定性顆粒軌道模型，在考慮顆粒間相互作用的情況下，采用曳力公式對顆粒與流體的相互作用進行描述。同時，基于牛頓第三定律可實現(xiàn)流場與顆粒的相互關(guān)聯(lián)，即在各個控制體中，顆粒受到流體的作用力等于顆粒對流體施加的作用力。

在傳統(tǒng)的離散顆粒模擬中，流體相的控制方程是基于體積平均Navier-Stokes (volume-averaged Navier-Stokes,VANS)方程或者Navier-Stokes方程的簡化形式，通常采用SIMPLE(semi-implicit method for pressure-linked equations)算法等求解，這些處理方式都不易實現(xiàn)并行化,難以應(yīng)用于大規(guī)模的計算中[4]?；贚BM的算法具有良好的并行性，王利民等[1,50]將LBM應(yīng)用于離散顆粒模擬中流體相的流動(圖5)，提出了基于格子 Boltzmann 的離散顆粒模擬( LB -based DPS)，其使用一種既考慮孔隙率又考慮氣固相對滑移速度的修正LBM求解流體流動，取代了基于VANS方程的求解傳統(tǒng)，固體運動采用TDHS求解，即固體顆粒的碰撞為硬球碰撞，在固定的時間步長內(nèi)檢索顆粒間是否發(fā)生碰撞以及更新顆粒運動。由于在顆粒濃度較高或顆粒溫度較低的場景下應(yīng)用硬球模型還有一定的應(yīng)用局限性，王利民等[51]進一步發(fā)展了一種介尺度LBM-DEM耦合模型。熊勤鋼等[52]也提出了一種快速的LBM-DEM耦合模型，但對力源項的添加方式有所不同。

圖5 D2Q9計算流體格點和固體顆粒[1]

(1)

式中v為格點處的流體速度，vs為顆粒在當前格點的平均速度，-i表示與i方向相反。修正后的格子Boltzmann方程為[1]

(2)

式中εs和εg分別為格子控制體內(nèi)部固體顆粒和氣體的體積分率，且εs+εg=1。

LB -based DPS成功模擬了二維單孔射流鼓泡床[1,50,52]和三維鼓泡床(圖6)的流態(tài)化過程，得到了與經(jīng)驗關(guān)聯(lián)式和實驗數(shù)據(jù)相吻合的模擬結(jié)果。李斌等[53,54]將 LB -based DPS進一步應(yīng)用于二維噴動床和鼓泡床的流化過程，獲得了與實驗數(shù)據(jù)相吻合的結(jié)果。以上均驗證了 LB -based DPS能夠得到流固系統(tǒng)更多的細節(jié)信息和內(nèi)在特性，此外，LB -based DPS在計算精度、耗費時間和效率之間能達到很好的平衡,可以直接擴展到與其他固相求解器耦合計算。

圖6 三維流化床離散顆粒模擬流動結(jié)構(gòu)的演化

4 雙流體模型

目前 LB -based TFM研究文獻較少，傳統(tǒng)LBM的顏色模型、偽勢模型、自由能模型和動力學模型只適用于不混溶的二元流體模擬，不能應(yīng)用于流體相互滲透的多相流雙流體層次模擬。王鐵峰等[55]率先開展利用LBM構(gòu)建雙流體模型的研究，其利用無量綱方法變換格子Boltzmann方程并添加壓力修正項來實現(xiàn)，但其恢復的宏觀方程中，密度為兩相的合密度而非分相密度；Sankaranarayanan等[56]從連續(xù)Boltzmann方程出發(fā)，考慮顆粒溫度方程，提出了一種隱式格子Boltzmann模型，此模型也沒有真正解決問題。實現(xiàn) LB -based TFM的關(guān)鍵在于構(gòu)造能恢復為VANS方程的格子Boltzmann模型，郭照立等[57]構(gòu)造了一個可用于求解VANS方程的格子Boltzmann模型，該模型的平衡分布函數(shù)和作用力項都包含孔隙率，作用力的計算需要給定組成多孔介質(zhì)的顆粒粒徑，適合模擬均勻孔隙率的多孔介質(zhì)內(nèi)部流動，對于具有孔隙率梯度的復雜多相流問題，顯式添加的離散力項通常導致模型不易收斂。Sungkon等[58]將VANS方程中帶有體積分數(shù)的項均移到方程右端作為力源項，變換得到帶有源項的N-S方程，并構(gòu)造能恢復帶有源項的N-S方程的LB方程，但是所得LB方程只有一階精度。宋飛飛等[59]類似地將VANS方程轉(zhuǎn)為帶有源項的N-S方程，利用帶有源項的作用力LBM模型發(fā)展了具有二階精度的隱式格子Boltzmann方程，處理較復雜，不易并行。王利民等[1]考慮了體積分數(shù)和相間滑移速度對流體相流動的影響，基于浸入邊界的思想，提出了修正的格子Boltzmann方程，用于描述氣相流動，嘗試拓展到雙流體模擬。張金鳳等[4]提出基于VANS方程的格子Boltzmann模型，通過Chapman-Enskog分析證明提出模型確實能夠恢復到VANS方程。Blais等[60]指出體系存在孔隙率梯度的情況下，直接添加力源項修正壓力梯度的LBM形式是不穩(wěn)定的，并通過重構(gòu)LBM中的碰撞算子以還原成VANS方程，其在大孔隙率梯度的情況下具備更好的穩(wěn)定性。H?cker等[61]則是通過修正遷移步驟以提高穩(wěn)定性，結(jié)合對流擴散的LBM實現(xiàn)了Rayleigh-Taylor不穩(wěn)定性的模擬，但如何解決TFM下顆粒的緊密堆積還需進一步研究。雖然存在多種思路構(gòu)建VANS方程的格子Boltzmann模型，但在使用上均存在一定局限性，將其應(yīng)用到兩相流的研究仍需進一步探索。

5 格子多相流體力學模擬軟件

LB -based PR-DNS和DPS基于的浸入運動邊界法以固含率將流體、顆粒和邊界統(tǒng)一描述，保留了LB求解的并行性和遷移的簡單性；同時格子Boltzmann方程基于笛卡爾網(wǎng)格求解能夠處理復雜邊界，避免了傳統(tǒng)CFD繪制復雜構(gòu)體網(wǎng)格的困難，在大規(guī)模并行條件下仿真規(guī)模已具備由實驗室尺度向工業(yè)尺度的可行性。基于此，中國科學院過程工程研究所EMMS團隊集成相關(guān)算法開發(fā)了以顆粒流仿真和多GPU并行為特色的高性能格子多相流體力學LMFD(lattice-based multi-fluids dynamics)模擬軟件，如圖7所示。LMFD基于C/C++開發(fā)，目前穩(wěn)定版本支持二維 LB -based PR-DNS和DPS，以及三維簡單構(gòu)型的 LB -based PR-DNS計算。

圖7 LMFD軟件界面

顆粒流體仿真的主流商業(yè)軟件有Barracuda，F(xiàn)luent和開源的軟件如MFiX和OpenFOAM。以上產(chǎn)品均基于傳統(tǒng)的有限體積法開發(fā)，在并行性上受到很多制約。而LB基的主流求解器，如PowerFlow，XFlow和Palabos均以求解單相流為主，而以多相流仿真為主的LBM求解器還比較少。Seil等[62]耦合了Palabos與LIGGGHTS兩個開源軟件，實現(xiàn)了 LB -based PR-DNS求解器LBDEMcoupling。G?tz等[63]在具備自適應(yīng)網(wǎng)格加密的單相流程序WaLBerla基礎(chǔ)上實現(xiàn)了 LB -based PR-DNS，用于研究河床中的沙丘形成。但以上程序的應(yīng)用場景較為有限，且均無法實現(xiàn) LB -based DPS。LMFD軟件通過構(gòu)建通用的LB基高性能氣固兩相流離散模擬平臺，再逐步形成完備的多相流求解模塊，將有望彌補這一空缺，在化工、冶金、能源、水利、環(huán)保和航空航天等相關(guān)領(lǐng)域具有廣泛的應(yīng)用前景。

6 結(jié) 論

LBM作為一種近幾十年發(fā)展的流體系統(tǒng)建模及模擬新方法，為顆粒流體系統(tǒng)的建模及計算提供了一種新思路。本文對LBM應(yīng)用于顆粒流體系統(tǒng)的進展作了較為系統(tǒng)總結(jié)，針對不同層次的顆粒流體計算框架，未來重點的研究方向如下。

(1) 對于格子Boltzmann顆粒解析直接數(shù)值模擬，隨著計算機硬件的發(fā)展，高雷諾數(shù)顆粒懸浮、顆粒湍流作用、復雜顆粒、耦合傳熱、傳質(zhì)及化學反應(yīng)等復雜過程等基礎(chǔ)研究將會受到廣泛關(guān)注。

(2) 對于格子Boltzmann離散顆粒模擬，擴大計算規(guī)模應(yīng)用于工業(yè)反應(yīng)器是未來發(fā)展趨勢，如顆粒相采用多相流網(wǎng)格質(zhì)點群(MP-PIC)方法中parcel的處理方式，或采用GPU等硬件加速計算等。

(3) 對于格子Boltzmann雙流體模擬，隨著LBM基礎(chǔ)理論研究深入開展，LBM代替?zhèn)鹘y(tǒng)有限元或有限體積方法求解雙流體模型方程在不久的將來肯定會得以實現(xiàn)。

致謝:感謝李靜海院士對顆粒多相流模擬工作的長期支持與指導；感謝葛蔚研究員和王小偉研究員在格子Boltzmann模擬方面長期友好合作；感謝曾參與本文涉及工作的熊勤鋼、周國峰、劉曉雯、魏民、張博和張金鳳等畢業(yè)同學；感謝陳建華博士對本文工作的積極建議和幫助。