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

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

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

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

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

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

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

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

四參數隨機生長法；格子Boltzmann方法；多孔介質；耦合傳熱

0 引言

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

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

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

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

1 四參數隨機生長法

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

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

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

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

由以上函數可以定義孔隙度ε為：

2 格子Boltzmann方法

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

2.1 速度場模型

速度場采用標準D2Q9模型，離散速度可表示為：

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

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

2.2 溫度場模型

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

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

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

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

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

3 問題描述

3.1 邊界條件

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

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

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

3.2 多孔介質模型簡化

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

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

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

4 結果與討論

4.1 模型驗證：

4.1.1 頂蓋驅動流

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

4.1.2 二維方腔自然對流傳熱

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

計算得到沿熱壁面的平均Nu數，并與Barakos[16]等的計算結果進行對比。從表1可以看出LBM的模擬結果和Barakos的計算結果吻合度較高。本文討論的多孔介質內的流體流動及傳熱問題，當流體在骨架間的孔隙內部流動時，孔隙間的骨架具有一定的溫度，其傳熱機理及處理方法與方腔自然對流類似，因此本算例說明了LBM處理多孔介質傳熱問題是可行的。

圖3 方腔內流函數示意圖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 結果與討論

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

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

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

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

表1 熱壁面上平均Nu數與基準解的比較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 多孔介質內流線圖(ε=0.75)Fig. 7 The contour map of velocity (ε=0.75)

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

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

圖9 平均速度與孔隙度(a)及壓力梯度(b)的關系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

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

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

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

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

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

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

5 結論

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

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

[2] 姚軍, 孫致學, 張凱, 等. 非常規(guī)油氣藏開采中的工程科學問題及其發(fā)展趨勢[J]. 石油科學通報, 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 方法的理論及應用[M]. 北京:科學出版社, 2009.[HE Y L, WANG Y, LI Q. Lattice Boltzmann Method: Theory and application[M]. Beijing: Science Press, 2009. ]

[6] 張磊, 姚軍, 孫海, 等. 基于數字巖心技術的氣體解析/擴散格子Boltzmann模擬[J]. 石油學報, 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]. 科學通報, 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方法的原理及應用[M].北京: 科學出版社, 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

吳子森, 董平川, 袁忠超, 周陰國, 曹耐, 許佳良. 多孔介質流動及傳熱的格子Boltzmann方法研究. 石油科學通報, 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