多孔介質(zhì)流動及傳熱的格子Boltzmann方法研究

吳子森，董平川*，袁忠超，周陰國，曹耐，許佳良

1 中國石油大學(xué)石油工程教育部重點實驗室, 北京 102249

2 中海油研究總院開發(fā)研究院, 北京 100027

多孔介質(zhì)流動及傳熱的格子Boltzmann方法研究

吳子森1，董平川1*，袁忠超2，周陰國1，曹耐1，許佳良1

1 中國石油大學(xué)石油工程教育部重點實驗室, 北京 102249

2 中海油研究總院開發(fā)研究院, 北京 100027

基于隨機四參數(shù)生長法構(gòu)造不同孔隙度的多孔介質(zhì)，采用稀疏矩陣存儲方式存儲內(nèi)部流體節(jié)點、流固邊界節(jié)點及物理邊界節(jié)點。通過頂蓋驅(qū)動流和方腔自然對流驗證了格子Boltzmann方法處理多孔介質(zhì)的傳熱問題的可行性?；诟褡覤oltzmann耦合傳熱模型，計算得到復(fù)雜多孔介質(zhì)內(nèi)的速度及溫度分布云圖，對多孔介質(zhì)內(nèi)的速度、溫度與迭代時間及孔隙度的關(guān)系作了詳細(xì)的分析。結(jié)果表明：多孔介質(zhì)中速度和溫度分布受孔隙度影響，相同壓力梯度下，平均速度隨著孔隙度的增加而增大；孔隙度相同時，平均速度隨著壓力梯度的增加而增大。相同溫差條件下，隨著孔隙度的增大，多孔介質(zhì)中對流傳熱達到穩(wěn)態(tài)的時間逐漸減?。欢嗫捉橘|(zhì)內(nèi)最大溫度出現(xiàn)在高溫壁面處，靠近高溫壁面的孔隙度越大，高溫部分體積越大。相同條件下，孔隙度越小，多孔介質(zhì)的平均溫度越低。

四參數(shù)隨機生長法；格子Boltzmann方法；多孔介質(zhì)；耦合傳熱

0 引言

多孔介質(zhì)是由固體骨架和流體組成的復(fù)合介質(zhì)，多孔介質(zhì)內(nèi)的流動和傳熱過程在自然界、人類生產(chǎn)及生活中廣泛存在[1]。從學(xué)科發(fā)展的角度看，多孔介質(zhì)流動及傳熱問題已經(jīng)滲透到許多學(xué)科和領(lǐng)域。地下巖層中的石油、天然氣和水是自然界多孔介質(zhì)中一種復(fù)雜的多元體系，研究油、氣的開采、特別是石油的熱采技術(shù)，促使石油工程學(xué)對多孔介質(zhì)的傳熱傳質(zhì)進行系統(tǒng)的研究。地?zé)豳Y源的勘測評估和開發(fā)利用以及利用土壤巖層作為蓄熱蓄冷介質(zhì)，也需應(yīng)用類似的理論與技術(shù)。因此對多孔介質(zhì)流動及傳熱的研究，是一項具有重大學(xué)術(shù)價值、對學(xué)科發(fā)展和技術(shù)創(chuàng)新具有深遠意義的研究課題[2]。

自然界中實際多孔介質(zhì)形狀多變而復(fù)雜，國內(nèi)外許多學(xué)者對其進行了研究，主要可以歸納為“連續(xù)介質(zhì)模型”和“非連續(xù)介質(zhì)模型”[3]。上述模型在應(yīng)用中遇到的主要問題是，由于多孔介質(zhì)內(nèi)部孔隙和骨架之間的邊界極為復(fù)雜，有時過于簡化，因此即使對模型求解也得不到有價值的信息。隨著計算機技術(shù)的發(fā)展，編制程序生成滿足目標(biāo)孔隙度的方法也越來越多，其中Wang M.R[4]等提出的四參數(shù)隨機生長法(QSGS)重構(gòu)多孔介質(zhì)，通過相關(guān)性比較其與實際多孔介質(zhì)符合很好，QSGS方法也成重構(gòu)多孔介質(zhì)的主要方法。

格子Boltzmann方法是上個世紀(jì)末期發(fā)展起來的一種流體系統(tǒng)建模及介觀尺度模擬的數(shù)值方法，與傳統(tǒng)的計算方法相比：其主要特點是物理圖像清晰、邊界條件容易處理及并行性能好等[5]。因而迅速在單相流[6-7]、湍流[8]、化學(xué)反應(yīng)流[9]、燃燒[10]、磁流體[11]、晶體生長[12]等領(lǐng)域得到了應(yīng)用。

本文采用QSGS重構(gòu)出不同孔隙度的多孔介質(zhì)，應(yīng)用稀疏矩陣存儲方式將孔隙相節(jié)點、流固邊界節(jié)點及物理邊界節(jié)點進行存儲。通過對頂蓋驅(qū)動流和二維方腔自然對流的模擬分析了格子Boltzmann方法在處理流體和固相耦合傳熱問題可行性。基于格子Boltzmann方法，分析了復(fù)雜多孔介質(zhì)內(nèi)部的速度、溫度與孔隙度和迭代時間的關(guān)系。

1 四參數(shù)隨機生長法

Wang M R[4]等人提出的QSGS生長步驟為：①根據(jù)給定的生長核分布概率，隨機播種生長核，②基于給定的生長概率Di(i表示方向)，生長核按照指定方向生長，如果相鄰節(jié)點上的隨機數(shù)不大于Di，則將該節(jié)點生長為骨架，③重復(fù)步驟②，直到多孔介質(zhì)孔隙度等于設(shè)定的孔隙度，④輸出結(jié)構(gòu)數(shù)據(jù)。本文根據(jù)上述方法在計算空間中生成滿足目標(biāo)孔隙度的二維隨機多孔介質(zhì)，如圖1所示：

圖1 多孔介質(zhì)重構(gòu)圖：ε=0.45(左)、ε=0.8(右)Fig. 1 Porous medium constructed by QSGS:ε=0.45(left)、ε=0.8(right)

圖1中紅色區(qū)域表示骨架部分，用1表示，藍色區(qū)域表示孔隙部分，用0表示。根據(jù)以上結(jié)構(gòu)圖，引入函數(shù)：

式中：表示空間位置，s指紅色骨架區(qū)域部分。

由以上函數(shù)可以定義孔隙度ε為：

2 格子Boltzmann方法

針對多孔介質(zhì)內(nèi)部流動及傳熱問題，本文采用雙分布函數(shù)模型：即對速度和溫度分別采用不同的分布函數(shù)進行描述。

2.1 速度場模型

速度場采用標(biāo)準(zhǔn)D2Q9模型，離散速度可表示為：

式中c表示格子速度，，Δx和Δt分別表示格子步長和時間步長，α表示9個速度方向。

無外力項的Boltzmann方程離散為如下格式：

2.2 溫度場模型

溫度場本文仍采用D2Q9模型，溫度場的分布函數(shù)演化方程描述如下：

其中g(shù)α表示時刻點處溫度分布函數(shù)，表示平衡態(tài)溫度分布函數(shù)，τg表示弛豫時間。

宏觀溫度T可根據(jù)分布函數(shù)對離度速度矢量的零階矩求得：

通過Chapman-Enskog多尺度展開可以得到溫度擴散方程：

式中χ表示有效熱擴散系數(shù)，方程(8)、(9)、(14)構(gòu)成了完整的流動及傳熱數(shù)學(xué)描述。

3 問題描述

3.1 邊界條件

文中模擬區(qū)域網(wǎng)格大小為200×200，內(nèi)部是孔隙度從0.45～0.80的多孔介質(zhì)，多孔介質(zhì)上下邊界為骨架，采用反彈格式處理，左、右邊界分別為入、出口，采用壓力邊界條件，入口溫度設(shè)為293 K，出口溫度設(shè)為0 K。邊界條件描述如下：

對于上述邊界均采用郭照立等[13]提出的非平衡外推格式：

式中：為邊界節(jié)點，為與邊界相鄰層節(jié)點。

3.2 多孔介質(zhì)模型簡化

一般地，流體在孔隙內(nèi)的流動仍然滿足N-S方程，將多孔介質(zhì)的固相看作是內(nèi)部流場的邊界，孔隙中的流體與固相邊界的相互作用通過邊界條件來實現(xiàn)。對于多孔介質(zhì)問題而言，固相部分占有很大的比重，同時固相內(nèi)部沒有流體流動。如果采用完全數(shù)組的數(shù)據(jù)結(jié)構(gòu)，將造成內(nèi)存的巨大損耗，并且也會因為對固體格點的循環(huán)操作和邏輯運算而降低計算效率。為此我們采用如圖2所述的處理方法，應(yīng)用一定的算法將多孔介質(zhì)內(nèi)部的流固邊界的第一層進行識別，從而在計算的時候只需對圖2中的藍色和綠色部分進行計算，紅色的部分直接初始化為0即可。這樣節(jié)省了計算機的內(nèi)存損耗，同時也降低了計算時間。

圖2 多孔介質(zhì)內(nèi)部節(jié)點存儲算法Fig. 2 Internal node storage algorithm in porous media

對于流固邊界，圖2中的綠色節(jié)點，采用碰撞反彈邊界處理格式：f1,2,5,6=f3,4,7,8，研究表明[14]，對于低馬赫數(shù)問題，其能嚴(yán)格的保證系統(tǒng)的質(zhì)量和動量守恒。

4 結(jié)果與討論

4.1 模型驗證：

4.1.1 頂蓋驅(qū)動流

采用頂蓋驅(qū)動流的算例，驗證格子Boltzmann 方法的可行性。設(shè)定Reynolds數(shù)Re=100，網(wǎng)格劃分為65×65，松弛時間τf=0.54，初始條件設(shè)定為最頂部節(jié)點的無量綱速度u=0.1，圖3為方腔內(nèi)流線圖，圖4為截面上y方向無量綱速度分布。從圖4看出LBM模擬結(jié)果與與Ghia等人[15]吻合的很好。

4.1.2 二維方腔自然對流傳熱

通過格子Botzmann方法模擬二維方腔自然對流，參數(shù)設(shè)定為：左側(cè)壁面無量綱溫度為20，右側(cè)壁面無量綱溫度為1，Prandtl數(shù)Pr=1，Rayleigh數(shù)Ra=103～105，網(wǎng)格劃分為100× 100，松弛時間τf=τg=0.54。圖5為流體達到穩(wěn)態(tài)時流線和等溫線圖。

計算得到沿?zé)岜诿娴钠骄鵑u數(shù)，并與Barakos[16]等的計算結(jié)果進行對比。從表1可以看出LBM的模擬結(jié)果和Barakos的計算結(jié)果吻合度較高。本文討論的多孔介質(zhì)內(nèi)的流體流動及傳熱問題，當(dāng)流體在骨架間的孔隙內(nèi)部流動時，孔隙間的骨架具有一定的溫度，其傳熱機理及處理方法與方腔自然對流類似，因此本算例說明了LBM處理多孔介質(zhì)傳熱問題是可行的。

圖3 方腔內(nèi)流函數(shù)示意圖Fig. 3 Schematic diagram of fl ow function

圖4x/L=0.5截面上沿Y方向無量綱速度分布Fig. 4 Non-dimensional velocity distribution along theYdirection (x/L=0.5)

4.2 結(jié)果與討論

從圖6和圖7可以看出，應(yīng)用格子Boltzmann方法可以得出多孔介質(zhì)內(nèi)部的詳細(xì)流動信息：包括孔隙間的流線圖、不同時間步長的溫度分布云圖，這是傳統(tǒng)計算方法在處理多孔介質(zhì)流動及傳熱問題時無法解決的難題。從圖6可以看出，多孔介質(zhì)內(nèi)最大溫度出現(xiàn)在高溫壁面處，靠近高溫壁面的孔隙率越大，高溫部分體積越大。

從圖8可以看出孔隙度不同，收斂時間步長不同，總體表現(xiàn)為孔隙度越小，收斂得越慢，原因是隨著孔隙度的減小，流體的流速減慢，同時骨架份額的增多導(dǎo)致孔隙區(qū)域變得曲折，許多區(qū)域流體需要一定的時間才能到達。另一方面是隨著孔隙度的減小，流固邊界節(jié)點越來越多，需要處理的邊界也越來越多，這無疑增大了計算量因此收斂的時間步長隨著孔隙度的減小增大。

平均速度是整個流場速度的統(tǒng)計平均，在統(tǒng)計上消除了局部最大速度帶來的偶然性。圖9(a)為平均速度與孔隙度的關(guān)系曲線，從圖上可以看出相同壓力梯度下，平均速度隨著孔隙度的增加而增大，這是由于孔隙度越大，流體的份額越多，骨架份額減少，骨架對流體流動的阻礙作用越低，使得流體的流速增加。圖9(b)為平均速度與壓力梯度的關(guān)系曲線，可以看出孔隙度不變時，平均速度隨著壓力梯度的增加而增大，這是因為壓力梯度增大，驅(qū)動力變大，促進了流體的流動。

圖5 不同Ra數(shù)下的流場(左)和等溫線(右)Fig. 5 Streamlines (left) and isotherms (right) for differentRavalues, from top to botton,Ra= 103, 104, 105

表1 熱壁面上平均Nu數(shù)與基準(zhǔn)解的比較Table 1 Comparison of average Nusselt numbers on the hot wall between the LBM results and other published data

圖6 不同時間步長的溫度分布云圖(ε=0.65)Fig. 6 The distribution of temperature fi eld at different time steps (ε=0.65)

圖7 多孔介質(zhì)內(nèi)流線圖(ε=0.75)Fig. 7 The contour map of velocity (ε=0.75)

從圖10、11可以看出，多孔介質(zhì)進出口的速度總體呈現(xiàn)進口速度大于出口速度的趨勢，但當(dāng)ε=0.8時，呈現(xiàn)出口速度大于進口速度的現(xiàn)象，這是由于出口端固體骨架的比重小于入口端固體骨架比重。從圖12可以看出，當(dāng)孔隙度較小時，壓降曲線在多孔介質(zhì)前半部分波動較為明顯，因為此時固體骨架所占比重比較大，骨架之間流道狹窄，流體流動阻力增大，導(dǎo)致壓降變化劇烈。

圖8 不同孔隙度下最大相對誤差(左)及平均速度(右)與迭代時間的關(guān)系Fig. 8 The relationship between the maximum relative error (left), average velocity (right) and iteration time with different porosity

圖9 平均速度與孔隙度(a)及壓力梯度(b)的關(guān)系Fig. 9 The relationship between average velocity、porosity (a) and pressure gradient (b)

圖10 不同孔隙度進出口速度：進口速度(左)、出口速度(右)Fig. 10 Velocity of inlet (left) and outlet (right) with different porosity

模擬孔隙度不同時，多孔介質(zhì)平均溫度隨時間的變化規(guī)律。從圖13可以看出孔隙度不同時，平均溫度隨迭代時間的變化曲線不同。在相同的迭代時間內(nèi)，隨著孔隙度增大，曲線斜率增大，整體的平均溫度也隨之升高，主要是因為隨著孔隙度的增大，骨架份額較少，相應(yīng)的孔隙相的比重增大，孔隙中的傳熱由于流體流動的存在屬于強化傳熱，同時孔隙度越大，平均流速越大，強化傳熱效果得到加強，因此使得曲線斜率增大及整體的平均溫度升高。

為了分析多孔介質(zhì)孔隙度對溫度的影響，進一步模擬了不同孔隙度下，多孔介質(zhì)的溫度變化。從圖14中可以看出孔隙度越小，多孔介質(zhì)的平均溫度越低，這是因為多孔介質(zhì)孔隙度越小，固相的份額越多，相應(yīng)的孔隙相的份額越低，而溫度由高溫壁面向低溫壁面?zhèn)鬟f的過程中，固相內(nèi)部只有熱傳導(dǎo)，而孔隙相則既有熱傳導(dǎo)，又有熱對流，使得傳熱程度相對較為劇烈，故而孔隙度越小的平均溫度比孔隙度較大的多孔介質(zhì)溫度低。

圖11ε=0.7時進出口速度Fig. 11 Velocity of inlet and outlet whenε=0.7

圖12 不同孔隙度下的壓力降落曲線Fig. 12 Pressure drop curves with different porosity

圖13 孔隙度不同時平均溫度與迭代時間步長的關(guān)系Fig. 13 The relationship between average temperature and iteration time step with different porosity

圖14 平均溫度與孔隙度的關(guān)系Fig. 14 Relationship between average temperature and porosity

5 結(jié)論

基于格子Boltzmann方法耦合傳熱模型，計算得到復(fù)雜多孔介質(zhì)內(nèi)的速度及溫度分布云圖，對多孔介質(zhì)內(nèi)的速度、溫度與迭代時間及孔隙度的關(guān)系作了詳細(xì)的分析。結(jié)果表明：流體的速度場和溫度場受多孔介質(zhì)孔隙度影響。相同壓力梯度下，平均速度隨著孔隙度的增加而增大；孔隙度相同時，平均速度隨著壓力梯度的增加而增大。相同溫差條件下，隨著孔隙度的增大，多孔介質(zhì)中對流傳熱達到穩(wěn)態(tài)的時間逐漸減小；多孔介質(zhì)內(nèi)最大溫度出現(xiàn)在高溫壁面處，靠近高溫壁面的孔隙度越大，高溫部分體積越大。相同條件下，孔隙度越小，多孔介質(zhì)的平均溫度越低。格子Boltzmann方法能處理多孔介質(zhì)的流動及傳熱問題，并能得出相關(guān)物理量隨孔隙度的變化關(guān)系，可以為處理更為復(fù)雜的多孔介質(zhì)傳熱及流動問題提供理論指導(dǎo)。

[1] DULLIEN F. Porous media fl uid transport and pore structure[M]. Utah: Academic Press, 1992.

[2] 姚軍, 孫致學(xué), 張凱, 等. 非常規(guī)油氣藏開采中的工程科學(xué)問題及其發(fā)展趨勢[J]. 石油科學(xué)通報, 2016, 1(1): 128-142. [YAO J, SUN Z X, ZHANG K, et al. Scienti fi c engineering problems and development trends in unconventional oil and gas reservoirs[J]. Petroleum Science Bulletin, 2016, 1(1): 128-142.]

[3] BEAR J. Dynamics of fl uids in porous media[M]. New York: Dover Publications, 1989.

[4] WANG M, WANG J, PAN N, et al. Mesoscopic predictions of the effective thermal conductivity for microscale random porous media[J]. Physical Review E, 2007, 75(2): 036 702-036 711.

[5] 何雅玲, 王勇, 李慶. 格子 Boltzmann 方法的理論及應(yīng)用[M]. 北京:科學(xué)出版社, 2009.[HE Y L, WANG Y, LI Q. Lattice Boltzmann Method: Theory and application[M]. Beijing: Science Press, 2009. ]

[6] 張磊, 姚軍, 孫海, 等. 基于數(shù)字巖心技術(shù)的氣體解析/擴散格子Boltzmann模擬[J]. 石油學(xué)報, 2015, 36(3): 361-365.[ZHANG L, YAO J, SUN H, et al. Lattice Boltzmann simulation of gas desorption and diffusion based on digital core technology[J]. Acta Petrolei Sinica, 2015, 36(3): 361-365.]

[7] 張磊, 康欽軍, 姚軍, 等. 頁巖壓裂中壓裂液返排率低的孔隙尺度模擬與解釋[J]. 科學(xué)通報, 2014, 59(32): 3197-3203.[ZHANG L, KANG Q J, YAO J, et al. The explanation of low recovery of fracturing fl uid in shale hydraulic fracturing by pore-scale simulation[J]. Chinese Science Bulletin, 2014, 59(32): 3197-3203.]

[8] CHEN H, SATHEESH K, STEVEN O, et al. Extended Boltzmann kinetic equation for turbulent fl ows[J]. Science, 2003, 301(5633): 633-636.

[9] CHEN S Y, DAWSON S P, DOOLEN G D, et al. Lattice Method and their applications to reaction systems[J]. Computers & Chemical Engineering, 1995, 19(4):395-408.

[10] YU H D, LUO L S, GIRIMAJI S. Scalar mixing and chemical reaction simulations using Lattice Boltzmann Method[J]. International Journal of Computational Engineering Science, 2012, 3(1):73-87.

[11] CHEN S, CHEN H, MARTNEZ D, et al. Lattice Boltzmann model for simulation of magnetohydro dynamics[J]. Physical Review Letters, 1991, 67(27): 3776-3779.

[12] MILLER W, SUCCI S, MANSUTTI D. Lattice Boltzmann model for anisotropic liquid-solid phase transition[J]. Physical Revive Letters, 2001, 86(16): 3578-3581.

[13] 郭照立,鄭楚光. 格子Boltzmann方法的原理及應(yīng)用[M].北京: 科學(xué)出版社, 2009.[GUO Z L, ZHENG C G. Theory and applications of Lattice Boltzmann Method[M]. Beijing: Science Press, 2009.]

[14] ZIEGLER D P. Boundary conditions for Lattice Boltzmann simulations[J]. Journal of Statistical Physics, 1993, 71(5): 1171-1177.

[15] GHIA U, GHIA K N, SHIN C T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[J]. Journal of Computational Physics, 1982, 48(3): 387-411.

[16] BARAKOS G, MITSOULIS E, ASSIMACOPOULOS D. Natural convection fl ow in a square cavity revisited: Laminar and turbulent models with wall functions[J]. International Journal for Numerical Methods in Fluids, 1994, 18(7): 695-719.

Study of flow and heat transfer in porous media based on the Lattice Boltzmann Method

WU Zisen1, DONG Pingchuan1, YUAN Zhongchao2, ZHOU Yinguo1, CAO Nai1, XU Jialiang1
1 MOE Key Laboratory of Petroleum Engineering, China University of Petroleum-Beijing, Beijing 102249, China
2 CNOOC Research Institute, Beijing 100027, China

A porous medium with different porosity is modelled based on the random four parameters growth method and the Sparse Matrix Storage Formats is used to store the internal fl uid node, the fl uid solid boundary node and the physical boundary node. The simulation of lid driven fl ow and natural convection in a square cavity are studied and the correctness of the Lattice Boltzmann model is verified. Based on the lattice Boltzmann method, the velocity and temperature field distribution in the complex porous medium are calculated, and the relationship between velocity or temperature and each of iteration time and porosity are analyzed in detail respectively. The results show that the distributions of fl uid velocity and temperature are affected by the porosity of the porous medium. With an increase of porosity, the average velocity increases under the same pressure gradient; with an increase of pressure gradient, the average velocity increases under the same porosity. At the same temperature, with the increase of porosity, the time convection heat transfer reaches a steady-state in the porous medium is gradually reduced. The maximum temperature appears at the high temperature surface, and the bigger the porosity near the high temperature wall, the greater the volume of the high temperature part. With a decrease of porosity, the average temperature of the porous medium decreases under the same condition.

quartet structure generation method; Lattice Boltzmann Method; porous media; coupled heat transfer

10.3969/j.issn.2096-1693.2017.01.008

(編輯 馬桂霞)

*通信作者, dpcofm@163.com

2016-07-09

吳子森, 董平川, 袁忠超, 周陰國, 曹耐, 許佳良. 多孔介質(zhì)流動及傳熱的格子Boltzmann方法研究. 石油科學(xué)通報, 2017, 01: 76-85

WU Zisen, DONG Pingchuan, YUAN Zhongchao, ZHOU Yinguo, CAO Nai, XU Jialiang. Study of fl ow and heat transfer in porous media based on the Lattice Boltzmann Method. Petroleum Science Bulletin, 2017, 01: 76-85. doi: 10.3969/j.issn.2096-1693.2017.01.008