張烈輝 賈鳴郭晶晶

(西南石油大學(xué)石油與天然氣工程學(xué)院，成都610500)

基于REV尺度格子Boltzmann方法的頁巖氣流動數(shù)值模擬

張烈輝 賈鳴1)E-mail:2585119904@qq.com郭晶晶

(西南石油大學(xué)石油與天然氣工程學(xué)院，成都610500)

結(jié)合頁巖掃描電鏡圖像，提出頁巖氣藏物理模型，采用表征單元體積(representative elementary volume,REV)尺度格子Boltzmann方法，考慮滑脫效應(yīng)，模擬頁巖氣在頁巖氣藏中的流動.模擬結(jié)果表明，頁巖氣主要沿著天然裂縫竄進(jìn)，但在有機(jī)質(zhì)和無機(jī)質(zhì)中也存在緩慢的流動，且有機(jī)質(zhì)中的流速要略大于無機(jī)質(zhì)中的流速.通過改變地層壓力，研究地層壓力對頁巖氣滲流特性的影響.研究結(jié)果表明，整個流場的速度和滲透率均隨著地層壓力的下降而增加.

格子Boltzmann方法,表征單元體積尺度,滑脫效應(yīng),頁巖氣,數(shù)值模擬

隨著常規(guī)油氣資源的減少和非常規(guī)資源開采技術(shù)的提高，頁巖氣在油氣田開發(fā)領(lǐng)域扮演的角色越來越重要[14].據(jù)前人研究表明，頁巖氣藏主要由基質(zhì)(有機(jī)質(zhì)、無機(jī)質(zhì))和天然裂縫組成，且絕大部分有機(jī)質(zhì)和無機(jī)質(zhì)中的孔隙均為納米級孔隙[5].頁巖氣在這樣微小孔道中流動，常常由于較大的克努森數(shù)而引起氣體分子在固體表面滑脫，從而引起滑脫效應(yīng).研究滑脫效應(yīng)對頁巖氣滲流影響對于頁巖氣藏的理論研究和工程實踐都具有重大的意義.

頁巖氣的流動數(shù)值模擬常常分為三個尺度：孔隙尺度，表征單元體積(representative elementary volume,REV)尺度和油田尺度[67].通常，在REV尺度上，基于控制偏微分方程的離散數(shù)值方法常被用來研究流動過程.孔隙尺度方法包括格子Boltzmann方法和直接蒙特卡洛方法.由于格子Boltzmann方法具有編程簡單，易于并行計算，且能夠處理復(fù)雜邊界條件的優(yōu)點，近幾年在頁巖氣數(shù)值模擬研究領(lǐng)域得到廣泛的應(yīng)用[8].但是，如果將孔隙尺度的格子Boltzmann方法直接運用到REV尺度上，計算量將會巨大[9].因此，研究REV尺度的格子Boltzmann方法具有很大的工程運用價值.

本文基于Guo等[7]提出的REV尺度格子Boltzmann方法，考慮滑脫效應(yīng)，研究頁巖氣在頁巖氣藏中的滲流特性.考慮實際生產(chǎn)過程中氣藏壓力不斷下降，本文還將研究壓力對頁巖氣藏滲流特性以及全局滲透率的影響.

1 考慮滑脫效應(yīng)的頁巖氣流動控制方程

Nithiarasu等[10]最初提出均質(zhì)多孔介質(zhì)瞬態(tài)滲流的Brinkman--Forchheimer--Darcy模型，后來Gao等[11]為了研究具有非均質(zhì)性的多孔介質(zhì)滲流問題，對上述模型進(jìn)行了修正，提出了非均質(zhì)模型.由于頁巖氣藏可以看作是由有機(jī)質(zhì)、無機(jī)質(zhì)和天然裂縫組成的非均質(zhì)場，故本文采用Gao等[11]的非均質(zhì)模型，其二維滲流的宏觀控制方程可以寫成

式中，▽2為拉普拉斯算子；u為流體的體積平均速度；t為時間;φi為第i個成分的孔隙度；ρ為流體密度；p為流體壓力；υe為有效黏度；F為流體在多孔介質(zhì)中所受的合力，包括線性項、非線性項和體積力項.為考慮頁巖氣藏中氣體流動的滑脫現(xiàn)象，在Gao的模型中引入視滲透率計算阻力

式中，υ為流體黏度，常與有效黏度υe取為相等；Fφi為多孔介質(zhì)的幾何形狀因子

Kai為第i個成分考慮滑脫效應(yīng)后的視滲透率，它與由多孔介質(zhì)結(jié)構(gòu)決定的絕對滲透率Ki和克努森數(shù)Kn有關(guān).其中，第i個成分的絕對滲透率Ki根據(jù)Kozeny--Carman(KC)方程[12]計算

其中，C為KC常數(shù)，由多孔介質(zhì)的固體顆粒直徑確定.根據(jù)Klinkenberg滑脫定律，視滲透率Ka為絕對滲透率K和滑脫因子fc的乘積.在本文中，滑脫因子fc由克努森數(shù)Kn確定，采用具有二階精度的Beskok--Karniadakis修正關(guān)系[13]計算

式中的b為滑移系數(shù)，稀疏因子α(Kn)采用Civan[14]提出的關(guān)系式計算

2 格子Boltzmann模型

本文采用Qian等[15]提出的D2Q9(二維九速)模型進(jìn)行計算.流體的密度分布函數(shù)fi(r,t)在格子空間的演化方程如下

式中，fi(r,t)為t時刻在格點r處沿i方向的流體密度分布函數(shù)；△t為格子空間的時間步長；τ是弛豫時間；ei為格子空間的離散速度，在D2Q9模型中，其表達(dá)式如下

式中c為格子速度，其值為格子步長和時間步長的比值，通常取為1.式(8)中為平衡態(tài)分布函數(shù)，在D2Q9模型中，當(dāng)考慮孔隙度的非均質(zhì)性時，計算如下

式中，wi為加權(quán)系數(shù)，并有w0=4/9，w1?4=1/9，w5?8=1/36；cs為格子聲速，且=c2/3；I為單位矩陣.Fi為外力在格子空間的離散項，本文依然采用Guo等[7]給出的表達(dá)形式

物理空間和格子空間之間的聯(lián)系由宏觀密度、宏觀速度和流體密度分布函數(shù)之間的關(guān)系體現(xiàn)，表示如下

根據(jù)式(3)可知，外力F的求解包含速度u，因此式(13)是一個非線性方程.Guo通過引入中間變量V，成功地給出了速度u的計算表達(dá)式，具體推導(dǎo)過程參考文獻(xiàn)[7].

3 算法驗證

圖1 經(jīng)典Poiseuille流模型示意圖

為了驗證本數(shù)值算法的可靠性，本文采用文獻(xiàn)[16]中的計算模型和計算參數(shù)，模擬壓力作用下的經(jīng)典Poiseuille流.文獻(xiàn)[16]中的模型如圖1所示，無量綱計算參數(shù)如下(計算量均為格子單位)：計算區(qū)域為L×H=200×50；孔隙度φ=1.0；滲透率K=∞；出入口壓力p由密度ρ確定，ρin=1.001，ρout=0.999；時間步長和格子步長均取一個格子單位；弛豫時間τ=1；左右邊界均為不滲透邊界(ux=uy=0)，在格子空間采用非平衡外推格式計算.當(dāng)流場達(dá)到穩(wěn)定的時候，繪制流場中間斷面流速uy的分布圖，并與經(jīng)典理論解和文獻(xiàn)[16]的結(jié)果對比.從圖2中可以看出，本文數(shù)值算法與理論解，以及文獻(xiàn)[8]的結(jié)果非常吻合，這驗證了本數(shù)值算法的正確性.

圖2 兩平板間y方向的速度剖面圖

4 頁巖氣流動模擬

4.1 頁巖氣滲流場分析

對某地區(qū)的頁巖巖心用掃描電鏡掃描，得到分辨率為5μm的圖像(圖3).從圖中可見，天然裂縫的寬度也不過在100nm左右.為了更好地研究頁巖氣在有機(jī)質(zhì)、無機(jī)質(zhì)和天然裂縫中的滲流情況，本文根據(jù)圖3提出了頁巖氣滲流多孔介質(zhì)模型(圖4，紅色表征天然裂縫，綠色表征有機(jī)質(zhì)，藍(lán)色表征無機(jī)質(zhì)).

圖3 頁巖掃描電鏡圖

圖4 頁巖多孔介質(zhì)模型圖

對圖4模型中的有機(jī)質(zhì)、無機(jī)質(zhì)和天然裂縫分別賦予不同的孔隙度，然后進(jìn)行格子Boltzmann數(shù)值模擬.其中，有機(jī)質(zhì)孔隙度φom=0.03，無機(jī)質(zhì)孔隙度φim=0.01，裂縫的孔隙度φf=0.1.其余的參數(shù)設(shè)定如下：計算區(qū)域為170μm×200μm；流體沿y正方向流動，在入口處的壓力pin=15.15MPa,出口的壓力pout=14.85MPa,其余區(qū)域壓力為pres= 15MPa，上下邊界為定壓邊界，左右邊界為封閉邊界；溫度T=323K；流體黏度υ=1.4×10?5Pa·s；氣體常數(shù)R=8.314J/(mol·K)；頁巖氣相對分子質(zhì)量Mr=16；克努森數(shù)Kn由分子平均自由程和有效孔隙半徑的比值確定；每個區(qū)域的視滲透率按照式(5)～式(7)確定.

當(dāng)計算達(dá)到收斂時，整個流場的速度場如圖5所示，可以看出，流體主要是沿著天然裂縫流動，在有機(jī)質(zhì)和無機(jī)質(zhì)中也存在緩慢流動.為了更好地區(qū)分流體在有機(jī)質(zhì)和無機(jī)質(zhì)中的差別，本文繪制了流速的等值線圖，如圖6.在等值線圖中，流速越大，等值線越密集，因此，流體在有機(jī)質(zhì)中的流速要稍微大于無機(jī)質(zhì)中的流速.圖7為y=146μm截面的速度剖面分布圖，其中，在4μm≤x≤28μm和141μm≤x≤150μm為有機(jī)質(zhì)，在60μm≤x≤65μm為天然裂縫，其余為無機(jī)質(zhì).該圖也印證了流體主要沿著天然裂縫竄進(jìn)，在有機(jī)質(zhì)中的速度略大于無機(jī)質(zhì)中流速這一結(jié)論.

圖5 流場穩(wěn)定時速度分布場圖

圖6 流場穩(wěn)定時速度等值線圖

4.2 地層壓力對流速和全局滲透率的影響

在頁巖氣藏的開采過程中，地層壓力會隨著開采過程而逐漸降低，因此很有必要研究地層壓力對滲流特性的影響.在4.1節(jié)給出的計算參數(shù)中，僅僅改變地層壓力和進(jìn)出口壓力，但保持壓差不變，其余參數(shù)均不改變，計算y=146μm處的速度剖面圖和全局滲透率.全局滲透率采用達(dá)西公式計算，即根據(jù)達(dá)到穩(wěn)定后出入口的流量和壓差計算.圖7為地層壓力取5MPa，10MPa，15MPa和20MPa時y=146μm處的速度剖面圖，從圖中可以得出，當(dāng)壓力逐漸減小的時候，各處的速度均在增大，這主要是因為壓力減小，導(dǎo)致克努森數(shù)Kn增大，滑脫現(xiàn)象變得嚴(yán)重，從而導(dǎo)致了速度增加.滑脫現(xiàn)象嚴(yán)重化也導(dǎo)致了全局滲透率的增加，如表1.這表明，隨著頁巖氣藏的開采，地層壓力降低將導(dǎo)致滑脫現(xiàn)象愈發(fā)明顯，從而導(dǎo)致全局滲透率增加.如果采用恒定的初始滲透率估算產(chǎn)量，將導(dǎo)致錯誤.

圖7 不同地層壓力下y=146μm處的速度剖面圖

表1 不同地層壓力下的全局滲透率

5 結(jié)論

(1)采用REV尺度的格子Boltzmann方法，考慮滑脫效應(yīng)和儲層非均質(zhì)性，并考慮有機(jī)質(zhì)、無機(jī)質(zhì)和天然裂縫差異，提出了頁巖氣流動數(shù)值模擬算法，并用經(jīng)典Poiseuille流動算例驗證了算法的正確性.

(2)頁巖氣在頁巖流動過程中主要沿著天然裂縫竄進(jìn)，在有機(jī)質(zhì)和無機(jī)質(zhì)孔隙中也存在緩慢流動.隨著地層壓力的降低，整體速度都存在上升的趨勢.

(3)隨著頁巖氣的開采，頁巖氣藏的地層壓力逐漸降低，這將導(dǎo)致地層中的克努森數(shù)增加，使滑脫效應(yīng)變得嚴(yán)重，從而使整個儲層的全局滲透率增加.如果采用初始的滲透率估算頁巖氣的產(chǎn)量，將導(dǎo)致錯誤.

1金衍，程萬，陳勉.頁巖氣儲層壓裂數(shù)值模擬技術(shù)研究進(jìn)展.力學(xué)與實踐，2016，38(1)：1-9

2朱維耀，鄧佳，楊寶華等.頁巖氣致密儲層滲流模型及壓裂直井產(chǎn)能分析.力學(xué)與實踐，2014,36(2)：157-160

3郭為,胡志明,左羅等.頁巖基質(zhì)解吸–擴(kuò)散–滲流耦合實驗及數(shù)學(xué)模型.力學(xué)學(xué)報,2015,47(6):916-922

4朱光普,姚軍,樊冬艷等.頁巖氣藏壓裂水平井試井分析.力學(xué)學(xué)報,2015,47(6):945-954

5 Wang JJ，Chen L，Kang QJ，et al.Apparent permeability prediction of organic shale with generalized lattice Boltzmann model considering surface di ff usion e ff ect.Fuel，2016，181：478-490

6張瀟丹，雍玉梅，李文軍等.REV尺度多孔介質(zhì)格子Boltzmann方法的數(shù)學(xué)模型及應(yīng)用的研究進(jìn)展.化工進(jìn)展，2016，35(6)：1698-1712

7 Guo ZL，Zhao TS.Lattice Boltzmann model for incompressible fl ows through porous media.Phys Rev E，2002，66:036304

8張磊，姚軍，孫海等.基于數(shù)字巖心技術(shù)的氣體解析/擴(kuò)散格子Boltzmann模擬.石油學(xué)報，2015，36(3)：161-365

9 Chen L，F(xiàn)ang WZ，Kang QJ，et al.Generalized lattice Boltzmann model for fl ow through tight porous media with Klinkenberg’s e ff ect.Phys Rev E,2015,91:033004

10 Nithiarasu P，Seetharamu KN，Sundararajan T.Natural convective heat transfer in a fl uid saturated variable porosity medium.International Journal of Heat and Mass Transfer，1997，40(16)：3955-3967

11 Gao JF，Xing HL，Tian ZW，et al.Lattice Boltzmann modeling and evaluation of fl uid fl ow in heterogeneous porous media involving multiple matrix constituents.Computers&Geosciences，2014，62：198-207

12 Carman PC.Fluid fl ow through granular beds.Trans Inst Chem Eng，1937，15:150-167

13 Beskok A，Karniadakis G.Report a model for fl ows in channels,pipes,and ducts at micro and nano scales.Microscale Therm Emg，1999，3(1)：43-77

14 Civan F.E ff ective correlation of apparent gas permeability in tight porous media.Transp Porous Med，2010，82：375-384

15 Qian YH，d’Humieres D，Lallemand P.Lattice BGK model for Navier--Stokes equation.Europhysics Letters，1992，17(6)：479-484

16申林方，王志良，李邵軍.基于格子博爾茲曼方法表征體元尺度土體細(xì)觀滲流場的數(shù)值模擬.巖土力學(xué)，2015，36(2)：689-694

(責(zé)任編輯:周冬冬)

NUMERICAL SIMULATION OF SHALE GAS FLOW BASED ON THE LATTICE BOLTZMANN METHOD AT REV SCALE

ZHANG LiehuiJIA Ming1)GUO Jingjing
(College of Petroleum Engineering,Southwest Petroleum University,Chengdu 610500,China)

During the seepage of the shale gas,its e ff ect is due to the large Knudson number,which makes the gas molecules slip from the solid surface.Based on the REV(representative elementary volume)scale lattice Boltzmann method,the fl ow of the shale gas is simulated with consideration of the slippage e ff ect.With the scanning electron microscope image,a physical model with organic rock,nonorganic rock and natural fractures is taken,and then the shale gas fl ow in the model is simulated.The simulated results indicate that the shale gas mainly fl ows in the natural fractures.However,the shale gas also fl ows slowly in the organic and nonorganic rocks,and the velocity in the organic rocks is larger than that in the nonorganic rocks.The in fl uence of the reservoir pressure on the shale gas seepage is also studied by varying the reservoir pressure.It is indicated that with the decrease of the reservoir pressure,the Knudson number increases and the slippage e ff ect intensi fi es, thus the velocity and the permeability of the whole fl uid fi eld increase.

lattice Boltzmann method,REV(representative elementary volume)scale,slippage e ff ect,shale gas,numerical simulation

TE371

A

10.6052/1000-0879-16-372

2016–11–15收到第1稿，2016–11–29收到修改稿.

張烈輝,賈鳴,郭晶晶.基于REV尺度格子Boltzmann方法的頁巖氣流動數(shù)值模擬.力學(xué)與實踐,2017,39(2):130-134 Zhang Liehui,Jia Ming,Guo Jingjing.Numerical simulation of shale gas flow based on the lattice Boltzmann method at REV scale.Mechanics in Engineering,2017,39(2):130-134