REV尺度多孔介質(zhì)格子Boltzmann方法的數(shù)學(xué)模型及應(yīng)用的研究進(jìn)展.pdf

CHEMICAL INDUSTRY AND ENGINEERING PROGRESS 2016 年第 35 卷第 6 期 1698 化 工 進(jìn) 展 REV 尺度多孔介質(zhì)格子 Boltzmann 方法的數(shù)學(xué)模型及應(yīng)用的 研究進(jìn)展 張瀟丹 1，2，雍玉梅2，李文軍3，趙元生4，李媛媛2，楊巧文1，楊超2 （1中國礦業(yè)大學(xué)（北京）化學(xué)與環(huán)境工程學(xué)院，北京 100083；2中國科學(xué)院過程工程研究所綠色過程與工程重 點(diǎn)實(shí)驗(yàn)室，北京 100190；3華北科技學(xué)院環(huán)境工程學(xué)院，河北 廊坊 065201；4中國石油化工研究院渣油加氫實(shí) 驗(yàn)室，北京 102200） 摘要：綜述了多孔介質(zhì)表征體元尺度（REV）格子 Boltzmann 模型的研究進(jìn)展，根據(jù)對多孔介質(zhì)處理方式主要分 為部分反彈模型和阻力模型兩類，分析歸納了各類模型的優(yōu)缺點(diǎn)。由于阻力模型中滲流的廣義格子 Boltzmann 方程（GLBE）的作用力是基于 GUO 等的作用力模型，可以準(zhǔn)確得到宏觀方程，不存在離散誤差，且模型的平 衡分布函數(shù)和作用力項(xiàng)中都包含反應(yīng)介質(zhì)特性的孔隙率，因而應(yīng)用最為廣泛。本文還重點(diǎn)介紹了 REV 尺度多孔 介質(zhì) LBE 模型在流動、傳熱、傳質(zhì)、化學(xué)反應(yīng)及相變等過程中的具體應(yīng)用，認(rèn)為 REV 尺度多孔介質(zhì)內(nèi)的三傳一 反數(shù)學(xué)模型中需要加入孔隙尺度因素，在更大工程尺度上應(yīng)該考慮過程參數(shù)的各向異性，展望了 REV 尺度多孔 介質(zhì) LBE 模型的發(fā)展和應(yīng)用前景。 關(guān)鍵詞：多孔介質(zhì)；表征體元尺度；格子 Boltzmann 方法；流動；傳熱；傳質(zhì) 中圖分類號：TQ021.9 文獻(xiàn)標(biāo)志碼：A 文章編號：10006613（2016）06169815 DOI：10.16085/j.issn.1000-6613.2016.06.010 Models and application of lattice Boltzmann method at REV-scale in porous media ZHANG Xiaodan1 ，2， YONG Yumei2， LI Wenjun3， ZHAO Yuansheng4， LI Yuanyuan2， YANG Qiaowen1， YANG Chao2 （1School of Chemical & Environmental Engineering，China University of Mining & Technology (Beijing)，Beijing 100083， China； 2Key Laboratory of Green Process and Engineering， Institute of Process Engineering， Chinese Academy of Sciences，Beijing 100190，China；3School of Environmental Engineering，North China Institute of Science and Technology，Langfang 065201，Hebei，China；4Laboratory of Residue Hydrotreating，Research Institute of Petroleum Processing，PetroChina，Beijing 102200，China） Abstract： This paper discusses the lattice Boltzmann model at representative elementary volume (REV) scale for porous media. According to different treatments of porous media， the lattice Boltzmann model at REV-scale for porous media can be classified into two categories，the partially bouncing-back model and the resistance model. The advantages and disadvantages of various models are analyzed. The Generalized lattice Boltzmann equation (GLBE model) in the resistance model is most widely used. Firstly， the force item of the GLBM model is based on the method proposed by Guo et al， which can be 第一作者：張瀟丹（1989） ，女，碩士研究生，主要從事化學(xué)工程 數(shù)值模擬。聯(lián)系人：雍玉梅，研究員，主要從事化學(xué)工程數(shù)值模擬 E-mail ymyong 。 特約評述 收稿日期：2016-01-19；修改稿日期：2016-01-31。 基金項(xiàng)目：國家自然科學(xué)基金（21276256,21490584） 、國家重點(diǎn)基礎(chǔ)研究發(fā) 展計(jì)劃（2012CB224806）及中央高?；究蒲袠I(yè)務(wù)費(fèi)（3142013097）項(xiàng)目。 第 6 期 張瀟丹等：REV 尺度多孔介質(zhì)格子 Boltzmann 方法的數(shù)學(xué)模型及應(yīng)用的研究進(jìn)展 1699 accurately recovered to the macroscopic equation without discretization error. Secondly，the equilibrium distribution function and the force items involve porosity， which reflects the characteristics of porous media. This review introduces the application of the lattice Boltzmann model at REV-scale in porous media with flow， heat transfer， mass transfer， chemical reaction and phase transition. We should take the pore scale factors into consideration when studying momentum transport， heat transport， mass transport and reaction engineering in porous media at REV-scale. The anisotropic nature of process parameters should be considered at the larger engineering scale. It gives some predictions and perspectives of applications for the lattice Boltzmann model at REV-scale. Key words：porous media；representative elementary volume (REV)；lattice Boltzmann method (LBM)；flow；heat transfer；mass transfer 流體在多孔介質(zhì)內(nèi)的流動是許多領(lǐng)域都涉及的 一類復(fù)雜流動，在能源化工、生物醫(yī)藥、空氣凈化、 污水處理等領(lǐng)域都有廣泛應(yīng)用。多孔介質(zhì)內(nèi)的流動 是一種典型的多尺度問題，通常涉及 3 個(gè)尺度，即 孔隙尺度、表征體元尺度（REV）和宏觀尺度。格 子 Boltzmann 方法（LBM）是介于宏觀和微觀兩種 方法之間的介觀數(shù)值模擬方法，它既克服了宏觀方 法需要處理復(fù)雜邊界的劣勢，又克服了微觀方法計(jì) 算尺度受限的局限性，在模擬多孔介質(zhì)方面具有獨(dú) 特優(yōu)勢1。 由于 LBM 算法簡單， 計(jì)算網(wǎng)格容易生成， 可以比較方便處理流體與邊界之間、不同流體組分 之間等復(fù)雜的相互作用， 且程序代碼簡單容易并行， 應(yīng)用逐漸廣泛。格子 Boltzmann 方法模擬多孔介質(zhì) 滲流主要有孔隙尺度和 REV 尺度?？紫冻叨?LBM 的研究對象是單個(gè)或者若干個(gè)孔隙內(nèi)的流體，把介 質(zhì)的骨架作為流場的邊界，用來探索滲流機(jī)理和基 本規(guī)律。REV 尺度 LBM 是以一個(gè)控制體為研究對 象， 求解 REV 體積平均的宏觀流動量， 流體與介質(zhì) 的相互作用通過適當(dāng)?shù)奶幚矸椒▉砻枋觥?REV 尺度 LBM 方法相對于孔隙尺度 LBM 方法，優(yōu)點(diǎn)是前者 不需要介質(zhì)的細(xì)節(jié)結(jié)構(gòu)、僅僅依賴于介質(zhì)的統(tǒng)計(jì)參 數(shù)，因而計(jì)算效率較高，比較適合大區(qū)域的工程滲 流計(jì)算。已有文獻(xiàn)總結(jié)了 LBM 方法的基本要素2， 介紹了 LBM 的起源、連續(xù) Boltzmann 方程、邊界 條件和數(shù)值穩(wěn)定性，討論了 LBM 在流體流動、多 相流和多組分流、反應(yīng)擴(kuò)散中的應(yīng)用。CHEN 等3 則對 LBM 在湍流、多孔介質(zhì)流、燃燒、多孔介質(zhì) 內(nèi)燃燒的應(yīng)用進(jìn)行了綜述， 李楠等4分析了 LBM 在 模擬這些復(fù)雜流體問題時(shí)存在的優(yōu)缺點(diǎn)。SUKOP 等5通過比較兩種 REV 尺度 LBE 模型， 提出 Dardis 和McCloskey的LBE模型不能用于各向異性的多孔 介質(zhì)的觀點(diǎn)，但沒有詳細(xì)論述兩模型的具體差異。 目前還沒有 REV 尺度多孔介質(zhì) LBM 方法的綜述， 因而本文總結(jié)了 REV 尺度多孔介質(zhì) LBE 各類數(shù)學(xué) 模型，分析了各模型的優(yōu)缺點(diǎn)，綜述了該方法在 REV 尺度多孔介質(zhì)內(nèi)流動、傳熱、傳質(zhì)、化學(xué)反應(yīng) 和相變等領(lǐng)域的應(yīng)用，分析了未來應(yīng)用前景，并對 其發(fā)展趨勢提出建議。 1 REV 尺度 LBE 模型 REV 尺度 LBE 模型根據(jù)對多孔介質(zhì)的處理方 式不同主要分為兩類：部分反彈模型6-11，此種 模型對經(jīng)典的固壁回彈模型進(jìn)行修改來擬合達(dá)西公 式；阻力模型12-21，此模型通過施加外力公式， 得到 Stokes 或 Brinkman 方程。 1.1 部分反彈模型 部分反彈模型由 3 個(gè)步構(gòu)成，如式(1)。前兩個(gè) 步是碰撞步和流動步， 與標(biāo)準(zhǔn) LBE 模型相同。第三 個(gè)步是部分反彈步，該過程中流體分布函數(shù)根據(jù)該 點(diǎn)的固體密度重新定義。部分反彈模型中引入了多 孔介質(zhì)連續(xù)表示方法，并將反彈邊界的作用耦合到 演化方程中7。 (eq) 1 (,)( , )( , )( , )( , ) ii ttiiii fctftftftt += + xxxxx (1) 式中，i 是速度離散方向；t 是離散時(shí)間；f 為 分布函數(shù)； 為松弛時(shí)間；i (x,t)是流體粒子與固 體骨架之間的碰撞，表示多孔介質(zhì)對流體影響。 最早的部分反彈模型是 GAO 和 SHARMA6提 出的多孔介質(zhì)格子氣自動模型，可延伸到 LBM 中 （GS-LBE 模型)。該模型能夠模擬非均質(zhì)多孔介質(zhì) 內(nèi)的流動。該模型中 i (x,t)如式(2)。 s ( , )=( , )( , ) ii i tnftft xxx (2) 式中，ns是固體格點(diǎn)的平均散射密度，取值范 圍為 01，代表多孔介質(zhì)的固體結(jié)構(gòu)；i 是 i 的相 化 工 進(jìn) 展 2016 年第 35 卷 1700 反方向。 但是 GS-LBE 模型存在格子氣的統(tǒng)計(jì)噪聲。 1998 年， DARDIS 和 MCCLOSKEY7則提出了第一 個(gè)部分反彈 LBE 模型（DM-LBE 模型)，是對 Brinkman 方程的間接求解。與 GS-LBE 模型相比， DM-LBE模型用ns代替孔隙尺度多孔介質(zhì)的結(jié)構(gòu)參 數(shù)布爾變量，消除了格子氣的不穩(wěn)定性。滲透率依 賴于固體格點(diǎn)的平均散射密度 ns，并且對應(yīng)的有效 黏性系數(shù)與流體的黏性系數(shù)相同。 該模型中為 i (x,t) 如式(3)。 s ( , )(, )( , ) iiti i tnftft=+ xxcx (3) THORNE 和 SUKOP8根據(jù)部分反彈形式不同 提出了第二類部分反彈格式（TS-LBE 模型） ，該模 型流出分布函數(shù)反映臨近點(diǎn)流入分布函數(shù)碰撞過 程，該類反彈格式需計(jì)算每個(gè)格點(diǎn)的碰撞步來估算 流出分布函數(shù)密度。該模型中 i (x,t)如式(4)。 c*c* s ( , )(,)( ,) ii ti i tnf tft =+ xxcx (4) 式中 t*是碰撞前的時(shí)間，滿足 tt*t+1，fic 為碰撞后的分布函數(shù)。 DM-LBE模型和TS-LBE模型的部分反彈格式 需要根據(jù)臨近格點(diǎn)計(jì)算得到，計(jì)算復(fù)雜，不利于并 行計(jì)算。WALSH、BURWINKLE 和 SAAR9引入 了新的部分反彈方法（WBS-LBE 模型） ，該方法 在碰撞步前改變流入流體方向， 同時(shí)計(jì)算碰撞步和 反彈步，且計(jì)算中不包含臨近格點(diǎn)的流體分布函 數(shù)，使計(jì)算更容易并行，從而節(jié)省計(jì)算時(shí)間。與 DM-LBE模型和TS-LBE模型的部分反彈方法相比 較， 該方法在非均質(zhì)多孔介質(zhì)流動過程中滿足質(zhì)量 守恒，且提高了 LB 方法數(shù)值模擬浮力驅(qū)動及應(yīng)力 和濃度梯度引發(fā)的擴(kuò)散過程的精度。該模型中 i (x,t)如式(5)。 cc* s ( , )( , )( ,) ii i tnftft = xxx (5) WBS-LBE 模型雖然在非均質(zhì)多孔介質(zhì)內(nèi)滿足 質(zhì)量守恒，更加適用于浮力驅(qū)動的流動和擴(kuò)散過 程， 但是該模型在不同滲透特性的多孔介質(zhì)界面處 速度不連續(xù)，ns接近 1 時(shí)會導(dǎo)致數(shù)值不穩(wěn)定。 YEHYA、NAJI 和 SUKOP10修正了該模型，提出 新的模型（YNS-LBE 模型）。該模型也在碰撞前 改變流入流體方向，但是公式中使用 i 方向的分布 函數(shù)代替i 方向的分布函數(shù)， 使其在界面處速度連 續(xù)，如式(6)。 * s ( , )( , )( ,) c iii tnftft = xxx (6) WBS-LBE 模型在 i 方向上的流體分布函數(shù)是 碰撞后的分布函數(shù)，而i 方向上的分布函數(shù)則是在 碰撞前的分布函數(shù)，導(dǎo)致反彈回的流體不僅沿其反 方向流體，也會向其他方向流動，引起數(shù)值不穩(wěn)定 性。 ZHU 和 MA11重新分配同一時(shí)刻同一格點(diǎn)碰撞 后方向相反的一對流體質(zhì)點(diǎn)，每個(gè)方向上碰撞后的 流體質(zhì)點(diǎn)分為兩部分，一部分沿該方向流向下一格 點(diǎn)，另一部分由于固體壁面的作用沿相反方向反彈 回去，提出新的部分反彈模型（ZM-LBE 模型） ，該 模型在同一時(shí)間的格點(diǎn)僅改變方向相反的一對流體 質(zhì)點(diǎn)，因此該模型滿足質(zhì)量守恒，且數(shù)值穩(wěn)定性比 WBS-LBE 模型好。該模型的松弛時(shí)間接近 1 時(shí)， 有效黏度系數(shù)與流體黏度系數(shù)相等。 該模型中 i (x,t) 如式(7)。 cc* s ( , )( , )( ,) ii i tnftft = xxx (7) 部分反彈模型區(qū)別主要在于如何處理反彈邊 界的作用、有效黏度、平均散射密度及有效滲透 率，總結(jié)目前主要部分反彈模型見表 1。表中包 括 DM-LBE 模型、 TS-LBE 模型、 WBS-LBE 模型、 YNS-LBE 模型和 ZM-LBE 模型， 這 5 個(gè)模型的滲 透率根據(jù)每個(gè)格點(diǎn)的平均散射密度 ns計(jì)算得到。 前 4 個(gè)模型在 ns=0 時(shí)，均可簡化為標(biāo)準(zhǔn) LBE 模 型. 在 ns=0 時(shí)，ZM-LBE 模型不能自動滿足反彈 邊界，因?yàn)樵撃Ｐ驮谕粫r(shí)間格點(diǎn)僅改變方向相 反的一對流體質(zhì)點(diǎn)。WBS-LBE 模型、YNS-LBE 模型和 ZM-LBE 模型在 ns=1 時(shí)能夠滿足全反彈 LBE 模型，其余模型不滿足全反彈 LBE 模型，因 為當(dāng) ns=1 時(shí)，DM-LBE 部分反彈模型在不透水邊 界處有質(zhì)量損失，而 TS-LBE 部分反彈格式模型 僅反映一個(gè)時(shí)間步的流體函數(shù)而不是原則上的兩 個(gè)時(shí)間步。其次，GS-LBE 模型、WBS-LBE 模型、 YNS-LBE 模型和 ZM-LBE 模型能夠滿足質(zhì)量守 恒，DM-LBE 模型和 TS-LBE 模型在高滲透性和 低滲透性多孔介質(zhì)邊界之間不對稱而不滿足質(zhì)量 守恒，需要在格點(diǎn)的邊緣設(shè)置 ns來修正模型，只 有在密度恒定或者平均散射密度不變時(shí)才滿足質(zhì) 量守恒。CHEN 等 22比較了 GS-LBE 模型、 DM-LBE 模型和 TS-LBE 模型，分析了在松弛時(shí) 間 =1 時(shí)， 為了保證流體黏度為正值， 3 個(gè)模型中 ns的上限值。GS-LBE 模型中 ns (1/23)0.289 DM-LBE模型中 ns (571/28)0.234，TS-LBE 第 6 期 張瀟丹等：REV 尺度多孔介質(zhì)格子 Boltzmann 方法的數(shù)學(xué)模型及應(yīng)用的研究進(jìn)展 1701 表 1 REV 尺度多孔介質(zhì)內(nèi)部分反彈 LBE 模型（f為孔隙流體黏度）11 模型 ( , ) i tx 有效黏度 ns（=1） 有效滲透率 GS-LBM s ( , )( , ) i i nftft xx () ss s 121(1)(21) 6 12 nn n + 0.289 f s 2n DM-LBM s (, )( , ) i ti i nftft+ xcx () ssss s 121(1)(21)(12) 6 12 nnnn n + + 0.234 f s 2n TS-LBM c*c* s (,)( ,) ti i nftft i + xcx ()()() () () s s s 212611 1 6 121 n n n + + 1 f s 2n WBS-LBM c* s ( , )( ,) i i nftft xx ()() ()() ss ss 21121 6 1211 nn nn + + 1 () sf s 1 2 n n YNS-LBM c* s ( , )( ,) ii nftft xx 1 - ZM-LBM cc* s ( , )( ,) ii nftft xx ()()() ()() 2 ss ss 268112121 6 12112 nn nn + + 0.5 () sf s 12 2 n n 模型和WBS-LBE模型中 ns在01之間。ZM-LBE 模型的 ns上限為0.5， 說明ZM-LBE模型最多有一 半的流體質(zhì)點(diǎn)可以反彈回相反的方向。綜上， YNS-LBE模型在 ns=0時(shí)，均可簡化為標(biāo)準(zhǔn)LBE 模型，在 ns=1時(shí)能夠滿足全反彈LBE模型，能夠 滿足質(zhì)量守恒，數(shù)值穩(wěn)定性也較好，因此應(yīng)用比 較好。 1.2 阻力模型 REV尺度多孔介質(zhì)的另一種LBE模型是阻力 模型12-21，即多孔介質(zhì)的影響通過施加外力得到 Stokes或Brinkman方 程 ， 這 個(gè) 思 想 是 基 于 NITHIARASU等23提出的求解不同孔隙度滲流的 Navier-Stokes方程，用其他CFD方法（有限差分、 有限體積）也可以方便求解，采用LBM方法是眾 多CFD方法的一種。該類模型的實(shí)現(xiàn)途徑有兩類： 第一類是修正速度12-17，即將多孔介質(zhì)的阻礙作用 通過修正速度來體現(xiàn)；另一種是流體與多孔介質(zhì)構(gòu) 造體-固體的相互作用力直接加入分布函數(shù)的演化 方程中18-21。 1.2.1 修正速度的阻力模型 修正速度的REV尺度多孔介質(zhì)LBE阻力模型 是在碰撞步中修正瞬時(shí)速度來描述多孔介質(zhì)對流體 的影響，所以其演化方程與標(biāo)準(zhǔn)BGK演化方程相 同12，如式(8)。 (eq) 1 (,)( , )( , )( , ) iittii ftftftft += xcxxx (8) 1997年SPAID和PHELAN12建立了直接求解 BRINKMAN方程的LBE模型（SP-LBE模型） 。該 模型將流體-固體的相互作用力的影響通過修正平 衡分布函數(shù)的速度來實(shí)現(xiàn)，這種作用力處理方式對 應(yīng)的宏觀流動方程中存在一些誤差，忽略誤差后， SP-LBE對應(yīng)的宏觀方程在穩(wěn)態(tài)時(shí)就是Brinkman方 程。SP-LBE模型平衡態(tài)分布函數(shù)的計(jì)算采用修正 后的平衡態(tài)速度，修正的平衡態(tài)速度包含了多孔介 質(zhì)的阻力作用，而阻力大小又與當(dāng)?shù)厮俣认嚓P(guān)，如 式(9)式(11)。 () 2 (eq) (eq)(eq) 2 (eq) 242 sss () ( , )1 22 i i ii ft ccc =+ cu cuu x (9) (eq) F =+uu (10) = Fu (11) 式中， v/K，是滲透勢能參數(shù)，是黏度 v 與 滲透率 K 的比值。 與SP-LBE模型類同，F(xiàn)REED13把多孔介質(zhì)阻 力通過碰撞過程在平均速度中體現(xiàn)Freed-LBE模 型，式(13)，可以用于Mach數(shù)較大的流動過程。 與多孔介質(zhì)流動控制方程相比，F(xiàn)reed-LBE模型對 應(yīng)的宏觀方程存在人工多余項(xiàng)，需要通過具體問題 來分析。該模型比較適合用于具有裂隙及導(dǎo)電介質(zhì) 的巖石孔隙中流體及其他工程模擬。Freed-LBE模 型公式如式(12)式(14)。 ()2 2 (eq) 242 sss ( , )1 22 i i ii ft ccc =+ cucuu x (12) 11 = 1) + 22 F =+( uuuu (13) 化 工 進(jìn) 展 2016 年第 35 卷 1702 Gu =u (14) () =Fuu (15) 式中，u為碰撞后速度， u 碰撞前速度，u 為 平均速度，G是速度張量， 取決于流體的運(yùn)動黏度、 松弛因子和滲透系數(shù)。 Shan-Chen模型14-15最初是應(yīng)用于多相流和多 組分流動，把體積力 F 在平衡速度中體現(xiàn)，體積力 F 主要包括液液流體間相互作用力、液固相間相互 作用力和外力。LI等16-17基于Shan-Chen模型提出 單相單組分LBM方法（Li-LBE模型） ，其平衡分 布函數(shù)可以簡化成式(16)，平衡態(tài)速度為式(17)。該 模型與其他模型相比，不用計(jì)算體積力 F，節(jié)省計(jì) 算時(shí)間，提高計(jì)算效率。 () () (eq) eq 2 s ,1 i ii ft c =+ cu x (16) 8 (eq)0 2(12 ) ii i f = =+ c uu (17) 8 1 0.5 1 2 i it i t f K = + = + cG u (18) 式中， K 為滲透率， 為孔隙率， G 為外部體積力。 1.2.2 修正演化方程的阻力模型 修正演化方程的阻力模型是將多孔介質(zhì)效應(yīng)直 接加入到演化方程，其演化方程19為式(19)。 (eq) 1 (,)( , )( , )( , ) iittiitii ftftftftF += + xcxxx (19) 式中，F(xiàn)i是 F 的函數(shù)，根據(jù)其處理方式不同， 可以分成不同的阻力模型。 MARTYS18對SP-LBE模 型 進(jìn) 行 改 進(jìn) （Martys-LBE模型） ， 采用LUO24的作用力處理方 法描述介質(zhì)阻力，在演化方程中直接增加了一個(gè)作 用力項(xiàng)來消除部分離散誤差， 但是不能完全消除， 如式(20)。其平衡態(tài)分布函數(shù)與SP-LBE模型相同， 即式(9)，但是作用力不修正平衡態(tài)分布函數(shù)的速 度，而是直接修正演化方程。Martys-LBE模型比 SP-LBE模型的數(shù)值穩(wěn)定性好，因?yàn)镾P-LBE模型 修正平衡分布函數(shù)的速度，增加了非物理項(xiàng)。 Martys-LBE模型其作用力表達(dá)式如式(20)此時(shí) c =Fa， i i f=。 22 ss () iii ii FF cc =+ cucu c (20) GUO等19在2002年提出了求解包含線性阻 力、 黏性項(xiàng)和非線性阻力項(xiàng)的通用滲流模型的REV 尺度滲流的廣義LBE模型（GLBE模型） 。該模型 的平衡分布函數(shù)和作用力項(xiàng)中均中包含反映介質(zhì)特 性的孔隙率 ，如式(21)和式(22)，而且對滲流速度 大小沒有要求。此模型采用GUO等25提出的作用 力模型， 不存在離散誤差， 可得到準(zhǔn)確的宏觀方程。 其中，F(xiàn)是包含介質(zhì)阻力和外部體積力的總力，式 (23)右端第一和第二項(xiàng)是線性（Darcy）和非線性 （Forchheimer）介質(zhì)阻力。F是結(jié)構(gòu)參數(shù)，與孔隙 率有關(guān)，對由固體顆粒構(gòu)成的多孔介質(zhì)，如 公式(24)。 ()2 2 (eq) 242 sss ( ) ( , )1 22 i i ii ft ccc =+ cucu u x (21) 2 s 24 ss :()1 (1) 2 ii i ii c F cc =+ cFuFc cI (22) F KK +F =uu uG (23) 3 1.75 150 F = (24) 目前還有一種趨勢，將LBM與其他數(shù)值方法 結(jié)合起來模擬多孔介質(zhì)內(nèi)流動和傳遞過程。 ZARGHAMI等20把GLBE模型和有限體積法相結(jié) 合，用有限體積方法離散對流擴(kuò)散項(xiàng)，如式(25)， 并用校正因子來提高數(shù)值穩(wěn)定性。該模型通過 Poiseuille流、Couette流和平板驅(qū)動流的驗(yàn)證，用 于模擬均質(zhì)和非均質(zhì)多孔介質(zhì)管道流。 () eq neqneqneq ,1,1 , neqneqneq 1,11,1 neqneqneq 1,11,11,1 1 d 11 48 1 16 ii acbd iii x yxyx y x y iii xyx yxy iii xyxyxy ffA fff A fff fff + + + = + + + (25) 式中，Ax.y為 abcd 的面積，ac、ab、bd、bc 分 別為格子邊界。 ZHANG等21提出了Navier-Stokes方程的顯性格 式解的LBE模型（Zhang-LBE模型） ，該模型為通過 Chapman-Enskog和泰勒展開，省略高階項(xiàng)之后，得 到平均體積Navier-Stokes方程，在演化方程式(26) 中加入了附加項(xiàng) Pi，如式(28)。該模型只在平衡分布 第6期 張瀟丹等：REV尺度多孔介質(zhì)格子Boltzmann方法的數(shù)學(xué)模型及應(yīng)用的研究進(jìn)展 1703 函數(shù)中加入反應(yīng)介質(zhì)特性的孔隙率，如式(27)，但是 與GLBE模型中孔隙率加入方式不同， 作用力項(xiàng)處理 方法采用式(29)，F(xiàn)公式與GLBE模型一致， 見式(23)。 該模型能有效地解決多孔介質(zhì)流和兩相流。 (eq) 1 (,)( , )( , )( , ) ii ttiiit iti f xctfx tfx tfx tPF += + (26) ()2 2 (eq) 242 sss ( ) ( , )1 22 i i ii u ft ccc =+ cuc u x (27) 24 ss 1 (1) 2 ii iii Pp cc =+ cucu c (28) 24 ss 1 (1) 2 ii iii F cc =+ cucu c F (29) 兩種類型多孔介質(zhì)REV尺度LBE模型，應(yīng) 用較為廣泛的是Freed-LBE模型和GLBE模型。 如果考慮不可壓流時(shí)，GLBE模型更加精確；對 于穩(wěn)態(tài)不可壓低速多孔介質(zhì)流動問題，兩者的區(qū) 別在于GLBE模型加入了非線性阻力項(xiàng)，在低速 情況下，該項(xiàng)的作用較弱26。GLBE模型的作用 力是基于GUO等的作用力模型， 可以準(zhǔn)確得到宏 觀方程，不存在離散誤差；其次，該模型的平衡 分布函數(shù)中包含反應(yīng)介質(zhì)特性的孔隙率，而其余 REV尺度LBE模型的平衡分布函數(shù)與標(biāo)準(zhǔn)模型 類似，不能反映介質(zhì)特性；最后，GLBE模型可 以用于高速流動。 2 REV 尺度 LBE 模型應(yīng)用現(xiàn)狀 2.1 REV尺度LBE模型應(yīng)用于多孔介質(zhì)內(nèi)流動 各REV尺度LBE模型在多孔介質(zhì)內(nèi)流動中均 有運(yùn)用，但是運(yùn)用最廣泛的是GLBE模型。 WBS-LBE模型、YNS-LBE模型和ZM-LBE模 型均被用于模擬不同滲透特性多孔介質(zhì)內(nèi)的流 動10-11，如圖1。模擬結(jié)果顯示W(wǎng)BS-LBE模型在 界面處不連續(xù)，YNS-LBE模型和ZM-LBE模型在 界面處速度連續(xù)，而且在 ns較小時(shí)的連續(xù)性比較好 （ns1=ns3=0.2，ns2=0.1） ，且多孔介質(zhì)的孔隙度變化 不會影響中心處的速度。 REV尺度LBE模型不僅被用于分層多孔介質(zhì) 內(nèi)速度的連續(xù)性， 還被用于多孔介質(zhì)內(nèi)Klinkenberg 效應(yīng)等基礎(chǔ)研究。孔隙大小在納米到微米范圍內(nèi)， 克努森數(shù) （Kn 數(shù)， 分子平均自由程和流動的宏觀特 征長度之比）相對較大時(shí)，氣體不再滿足連續(xù)介質(zhì) 假設(shè)，導(dǎo)致氣體滑移流、瞬變流，甚至真空分子流。 圖1 不同滲透特性多孔介質(zhì)示意圖及其速度分布圖10-11 由于存在滑移效應(yīng)，測得表觀滲透率高于固有滲透 率，兩者之差隨著 Kn 數(shù)增大而增大。CHEN等27 基于四參數(shù)隨機(jī)生成方法（QSGS）并運(yùn)用基本元 素建模模型（EBB模型）重新構(gòu)造頁巖，采用 BESKOK和KARNIADAKIS28提出的表觀滲透率 公式，在 Kn 數(shù)為常數(shù)時(shí)根據(jù)固有滲透率計(jì)算表觀 滲 透 率 ， 運(yùn) 用GLBE模 型 研 究 了 二 維 含 有 Klinkenberg效應(yīng)的多孔介質(zhì)內(nèi)流動。高孔隙率 （=0.8）時(shí)，Klinkenberg效應(yīng)可以忽略，隨著孔隙 率減小，Klinkenberg效應(yīng)在低滲透多孔介質(zhì)比較明 顯， 并且隨著 Kn 數(shù)增大、Klinkenberg效應(yīng)越明顯， 化 工 進(jìn) 展 2016年第35卷 1704 隨著壓力的減小、Klinkenberg效應(yīng)更加明顯 （圖2） 。 隨后CHEN等29通過Darcy定律預(yù)測REV尺度的 整個(gè)計(jì)算區(qū)域內(nèi)的有效滲透率， 研究多孔介質(zhì)成分、 顆粒大小、 壓力、 滑移及吸附作用對REV尺度滲透 率的影響。 在致密儲氣巖層或頁巖層開采天然氣時(shí)， 如果不考慮Klinkenberg效應(yīng)， 用表觀滲透率代替固 有滲透率，不僅會導(dǎo)致定量錯(cuò)誤而且會導(dǎo)致定性 錯(cuò)誤。 圖2 含有Klinkenberg效應(yīng)和不含Klinkenberg效應(yīng)的 速度分布圖27（1psi=6894.76Pa） 除了用于致密巖石中滲流， 土壤滲流也是REV 尺度多孔介質(zhì)LBE模型的研究熱點(diǎn)。 樊火30通過繪 制滲流場幾何剖面圖或者讀入圖形文件，格點(diǎn)化流 場、生成LBM格式的邊界數(shù)據(jù)文件，基于GLBE 型開發(fā)了土壤滲流LBM_Seepage軟件平臺， 生成數(shù) 據(jù)文件，以直觀可視的等值線或者云圖實(shí)現(xiàn)滲流場 速度實(shí)時(shí)可視化，并對豎向具有不同滲透率的土體 進(jìn)行了模擬。結(jié)果表明，最上層土的相對的滲透能 力在滲透初期對滲流速度的傳播起控制作用，只有 經(jīng)過一段時(shí)間后土體才會表現(xiàn)出等效應(yīng)的滲透能力 （圖3） 。 申林芳等31進(jìn)一步運(yùn)用GLBE模型探討了 土體在壓力作用下孔隙率、滲透率及滲透壓力等影 響因素與滲流速度的相互關(guān)系。 Freed-LBE模型被推廣到非均勻網(wǎng)格中32，采 用HE等33提出的插入格式描述非均勻網(wǎng)格，該模 型在流動和碰撞過程之間采用插值法，不僅能保證 數(shù)值準(zhǔn)確性而且易于實(shí)施。該模型對多孔介質(zhì)、孔 隙和邊界的處理方式一樣，避免了不同尺度采用不 同數(shù)學(xué)方法而產(chǎn)生的誤差，與其他多尺度LB模型 相比，它消除了內(nèi)部邊界，而且采用阻力模型，因 此可以用于大尺度。KANG等32在裂隙周圍細(xì)化網(wǎng) 格研究裂隙對流動的影響，遠(yuǎn)離裂隙的地方采用粗 網(wǎng)格（圖4） ，模擬了復(fù)雜裂隙和簡化裂隙的多孔介 質(zhì)內(nèi)的流動，通過 u l K p = 預(yù)測REV尺度多孔介質(zhì) 內(nèi)的有效滲透率（表2） ，并提出可以把該模型用曲 線坐標(biāo)和復(fù)合網(wǎng)格的設(shè)想。 REV尺度LBE模型不僅可以模擬含裂隙的均 質(zhì)巖石中的流動，還能根據(jù)通過高分辨率的掃描圖 像確定空腔、黏土、長石、石英不同成分構(gòu)成的非 均質(zhì)砂巖結(jié)構(gòu)，如圖5(a)，基于GLBE模型獲得砂 巖樣品內(nèi)的流體流動速度，然后根據(jù)達(dá)西定律求得 整體的滲透率34-35。由圖5(b)可以看出，樣品大于 圖3 不同滲透率的的邊界示意圖及滲透率 K 由小到大和 K 由大到小中心線速度圖30 第6期 張瀟丹等：REV尺度多孔介質(zhì)格子Boltzmann方法的數(shù)學(xué)模型及應(yīng)用的研究進(jìn)展 1705 圖4 裂隙示意圖、非均勻網(wǎng)格示意圖、計(jì)算結(jié)果32 表2 滲透率值32（格子單位） 項(xiàng)目 上半部分 下半部分 裂隙 整體 傳統(tǒng) LBE 模型 0.64 0.64 35.31 1.94 Freed-LBE 模型 0.61 0.61 37.01 1.99 誤差 4.7% 4.7% 4.8% 2.6% 10mm時(shí)整體的滲透率趨于常數(shù)，對于砂巖樣品， 大于10mm時(shí)可以達(dá)到REV尺度的要求。 由圖5(c) 和圖5(d)可以看出，流場主要受無滲透性的石英和 空腔影響，由于各成分的孔隙率不同導(dǎo)致各成分邊 界處的流線不閉合。但是CT掃描大孔時(shí)容易忽略 狹小的孔喉結(jié)構(gòu)而造成誤差，因而掃描圖像的分辨 率必須要高36。 REV尺度多孔介質(zhì)模型除了可以通過滲流速 度預(yù)測多孔介質(zhì)的滲透率外，還可以分析多孔介質(zhì) 圖5 砂巖樣品 幾何結(jié)構(gòu)特性及外部條件對滲流速度的影響。 GLBE模型被用于在壓力驅(qū)動下的非牛頓流體的電 滲 流 過 程 37 ， 在 演 化 方 程 中 加 入 了 基 于 Herschel-Bulkley模型的非牛頓流體力項(xiàng)， 研究外部 電勢驅(qū)動流、離子濃度和電場強(qiáng)度、固體顆粒直徑 和孔隙度、冪律指數(shù)和屈服應(yīng)力、德拜長度對滲透 速度的影響。多孔介質(zhì)內(nèi)的流速隨著冪律指數(shù)和屈 化 工 進(jìn) 展 2016年第35卷 1706 服應(yīng)力的減小而增大，壁面處的電勢主要影響壁面 處流速，多孔介質(zhì)的電勢對通道中心處流速影響較 大。低冪律指數(shù)時(shí)，壁面處的電勢在對流速的影響 更加敏感，特別是在通道高度與德拜長度的比值較 小時(shí)。 HUANG等38對LBM兩相流模型及其應(yīng)用 做了詳細(xì)描述， 目前LBM兩相流模型已經(jīng)應(yīng)用于 REV尺度多孔介質(zhì)內(nèi)多相流的驅(qū)替過程中。 SPAID和PHELAN39對纖維多孔介質(zhì)的流動進(jìn)行 了計(jì)算，纖維內(nèi)部使用SP-LBE模型，外部使用 標(biāo)準(zhǔn)LBE模型，模擬圓形纖維內(nèi)樹脂驅(qū)替空氣多 組分流動過程。纖維間的流動比纖維內(nèi)的流動速 度快，導(dǎo)致纖維內(nèi)部含有空氣空穴，空穴最初是 長軸垂直于流動方向的橢圓形，隨后由于表面張 力的影響，從橢圓形變?yōu)閳A形（圖6） 。根據(jù)達(dá)西 定律預(yù)測的不飽和滲透率在定性上正確。但是需 要更多研究證明纖維內(nèi)的空洞隨時(shí)間消失模擬結(jié) 果是否準(zhǔn)確。SCHAAP等 40應(yīng)用阻力模型中 Shan-Chen格式多相流LBE模型，模擬了充滿玻 璃珠的多孔介質(zhì)內(nèi)水-空氣及水-石油精體系驅(qū)替 實(shí)驗(yàn)。該模型基于顯微層析構(gòu)造的多孔介質(zhì)幾何 結(jié)構(gòu)，忽略重力、慣性力和黏性力，只考慮毛細(xì) 作用力，設(shè)置合理的表面張力和接觸角，得到壓 力與飽和度的關(guān)系。水-空氣體系模擬值與測量值 相差不大，但是水-石油精體系模擬值與測量值差 異很大（圖7） ，原因是模擬條件與實(shí)驗(yàn)條件不一 致，如接觸角。該方法也可以用于微尺度界面現(xiàn) 象的模擬。PORTER等41也運(yùn)用該LBE模型估算 潤濕-非潤濕相界面的面積，模擬值與實(shí)驗(yàn)值 吻合。 圖6 不同時(shí)刻單行圓形纖維樹脂驅(qū)替空氣模擬結(jié)果圖39 （黑色：空氣；灰色：樹脂） 圖7 壓力與飽和度演化圖 40 阻力模型中的單松弛GLBE模型推廣到多松弛 模型，LIU等42運(yùn)用二維D2Q8 MRT-GLBE模型模 擬不可壓多孔介質(zhì)流，在低黏性時(shí)的數(shù)值穩(wěn)定性比 單松弛模型好。隨后該模型被用于解含有時(shí)間導(dǎo)數(shù) 和非線性對流項(xiàng)平均體積動量方程，模擬了在多孔 介質(zhì)圓柱靜止時(shí)上下壁面以恒定速度運(yùn)動和上下壁 面靜止時(shí)管道內(nèi)圓柱多孔介質(zhì)以恒定速度運(yùn)動的兩 種情況下的流動43。 結(jié)果表明R1和R2的本征相平 均速度吻合，二者的相平均速度在流體區(qū)吻合，但 是在多孔介質(zhì)區(qū)域內(nèi)有偏差圖8(b)，并且隨著孔 隙度的減小而增大。多孔介質(zhì)阻力對本征相平均速 度和相平均速度有類似規(guī)律。只有含有本征相平均 速度的宏觀方程滿足伽利略不變性。多松弛模型與 單松弛模型相比，克服了許多限制條件，比如運(yùn)動 黏度和體積黏度比不再是定值，但是多松弛模型增 加了計(jì)算量。 REV尺度LBM已經(jīng)應(yīng)用到工程滲流問題。何 瑩松44將REV尺度GLBE模型應(yīng)用于煤礦開采過 程中回采工作面瓦斯?jié)B流問題，該研究中認(rèn)為煤體 的孔隙度受應(yīng)力重新分布的影響很小，孔隙度基本 不發(fā)生變化，滲透率沿流動方向呈現(xiàn)分段函數(shù)的不 均勻分布，得到了不同時(shí)刻非均質(zhì)煤層中瓦斯壓力 分布和速度分布。這些數(shù)據(jù)能夠生動再現(xiàn)瓦斯流動 第6期 張瀟丹等：REV尺度多孔介質(zhì)格子Boltzmann方法的數(shù)學(xué)模型及應(yīng)用的研究進(jìn)展 1707 圖8 多孔介質(zhì)圓柱運(yùn)動示意圖及=0.7、Re=100通道中心 固有相平均速度、相平均速度圖42 過程，有助于指導(dǎo)煤礦安全生產(chǎn)。 2.2 REV尺度LBE模型應(yīng)用于多孔介質(zhì)內(nèi)流動與 傳熱耦合過程 流動、 傳熱和傳質(zhì)兩個(gè)及兩個(gè)以上過程耦合時(shí)， LBE模型可以分為3種；即多速模型、雙分布函數(shù) 模型（DDF）與差分方法相結(jié)合的混合模型。目前 文獻(xiàn)大多采用DDF模型來模擬多過程， 該模型使用 兩個(gè)（及以上）分布函數(shù)演化方程，分別用于速度 場和溫度場（濃度場） ，格子結(jié)構(gòu)簡單，且具有良好 的數(shù)值穩(wěn)定性。 等溫GLBE模型通過DDF模型推廣到多孔介 質(zhì)熱流動中（DDF-GLBE模型）45，速度場用等溫 GLBE模型模擬，而溫度場使用一個(gè)新的分布函數(shù) 描述。DDF-GLBE模型被用于多孔介質(zhì)方腔模擬孔 隙度恒定45-49、孔隙度變化50、含有熱源51的速度 場和溫度場，研究孔隙度、雷諾數(shù)（Re 數(shù)） 、達(dá)西 數(shù)（Da 數(shù)）等參數(shù)對速度場和溫度場的影響。 DDF-GLBE模型還被用于研究邊界充滿多孔介質(zhì) 和中心充滿多孔介質(zhì)兩種情況下管道內(nèi)的流動與傳 熱過程52-53。 REV尺度的GLBE模型可以直接將孔隙尺度模 擬得到的滲透率和有效熱導(dǎo)率帶入REV模型， 避免 了采用經(jīng)驗(yàn)公式可能得到的錯(cuò)誤結(jié)論54。模擬結(jié)果 表明，REV尺度得到的熱壁上的平均Nusselt數(shù)與 孔隙尺度下的趨勢一致， 證明REV尺度在定性上是 正確的。在瑞利數(shù)（Ra 數(shù)）低于107時(shí)，GLBE模 型能定量預(yù)測孔隙尺度下自然對流的宏觀規(guī)律，Ra 數(shù)大的時(shí)候只能定性預(yù)測孔隙尺度的宏觀規(guī)律（圖 9） 。REV尺度GLBE模型對滲透率如何影響多孔介 質(zhì)自然對流宏觀現(xiàn)象刻畫很好，與孔隙尺度模擬結(jié) 果一致，但是往往夸大了流固熱導(dǎo)率帶來的影響。 圖9 不同Ra數(shù)下熱壁上平均Nusselt數(shù)隨方柱個(gè)數(shù)N的變 化結(jié)果對比54 虛線REV；實(shí)線孔隙尺度 REV尺度的LBE模型為太陽墻工程應(yīng)用研究 提供參考， 通過含熱源項(xiàng)的LBE模型對太陽墻系統(tǒng) 中復(fù)雜流場及熱傳遞進(jìn)行模擬55-56，可以預(yù)測太陽 墻在不同工作條件下的狀況，并優(yōu)化太陽墻系統(tǒng)的 設(shè)計(jì)。 針對工程中常見的多孔介質(zhì)圓管內(nèi)傳熱問題， 榮 伏 梅 等 57-58 構(gòu) 造 了REV尺 度 軸 對 稱 熱 DDF-GLBE模型，包含弱可壓軸對稱LBE模型和 不可壓軸對稱LBE模型， 把刻畫多孔介質(zhì)的參數(shù)加 入到分布函數(shù)中并構(gòu)造合適的分布函數(shù)，使之能恢 復(fù)正確的宏觀方程，具有形式簡單、外力項(xiàng)不含復(fù) 雜梯度形式的優(yōu)點(diǎn)， 并應(yīng)用于非線性滲流數(shù)值模擬。 分別研究了覆蓋多孔介質(zhì)的方柱繞流、充滿多孔介 質(zhì)內(nèi)的變黏性Rayleigh-Benard對流、 以及圓管內(nèi)插 入多孔介質(zhì)強(qiáng)化對流傳熱非線性滲流問題。由于這 類問題的計(jì)算量都很大，單個(gè)CPU的計(jì)算效率較 低，因此采用GPU（圖形處理器）加速計(jì)算，可支 化 工 進(jìn) 展 2016年第35卷 1708 持復(fù)雜的計(jì)算任務(wù)59。 單松弛DDF-GLBE模型推廣到多松弛模型中， 建立REV尺度二維多松弛（MRT）DDF-GLBE方 程，用D2Q9 MRT-GLBE模型模擬流場，同時(shí)用 D2Q5 MRT-LBE模型模擬溫度場， 模擬結(jié)果和解 析解與文獻(xiàn)計(jì)算結(jié)果吻合60，而且在低黏度時(shí)， MRT-LBE模型的數(shù)值穩(wěn)定性比LBGK模型好，但 是計(jì)算量大。該直角坐標(biāo)系下的MRT-DDF-GLBE 模型延伸到軸對稱模型61。 2.3 REV尺度LBE模型應(yīng)用于多孔介質(zhì)內(nèi)流動、 傳質(zhì)和化學(xué)反應(yīng)耦合過程 與流動和傳熱過程類同， 仍使用DDF模型描述 流動和傳質(zhì)兩個(gè)過程，即運(yùn)用兩個(gè)分布函數(shù)分別模 擬流場和濃度場，并在濃度演化方程中加入反應(yīng)源 項(xiàng)， 用于模擬化學(xué)反應(yīng)。TIAN等62基于DDF-GLBE 模型模擬了CO2地下填埋地球化學(xué)反應(yīng)過程，研究 CO2注入過程中溶質(zhì)離子反應(yīng)及擴(kuò)散、巖石骨架 CaCO3溶解相互作用， 引起孔隙率和滲透率變化， 反作用于流場。如圖10，開始在一段時(shí)間內(nèi)出口處 的H+、Ca2+和HCO3離子濃度為0，隨著反應(yīng)的進(jìn) 行，3種離子濃度逐漸增加。由于反應(yīng)消耗和對流 擴(kuò)散的相互作用，在出口處反應(yīng)物H+離子濃度在 1105時(shí)出現(xiàn)變化。Ca2+離子濃度到達(dá)最大后保持平 衡，是因?yàn)镃a2+離子是可逆反應(yīng)。 該模型還用于研究有機(jī)廢水溶液中化學(xué)反應(yīng)、 燃料電池和質(zhì)子交換膜燃料電池中56 ，63-66。LIAO 等63和楊艷霞64研究了有機(jī)廢水中含擾流光合細(xì) 菌包埋顆粒的流動、傳質(zhì)及其內(nèi)部的光生化反應(yīng)過 程，分析光照強(qiáng)度、進(jìn)口流速、包埋顆粒的滲透率 及孔隙率對流場、濃度場及產(chǎn)氫性能的影響。XU 等65-66把該模型應(yīng)用于固體燃料電池多組分化學(xué)反 應(yīng)，計(jì)算濃差極化電勢，并研究了電池電極厚度、 電解質(zhì)濃度、孔隙率、燃料組成等結(jié)構(gòu)參數(shù)和工作 圖10 出口離子濃度圖62 溫度操作參數(shù)對濃差極化、超電勢、電化學(xué)效率、 電池電壓和功率密度的影響，LBM模擬結(jié)果與其他 模擬結(jié)果相比更接近實(shí)驗(yàn)值。該模型也被用于模擬 質(zhì)子交換膜