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

張烈輝 賈鳴郭晶晶

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 格子Boltzmann模型

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

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

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

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

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

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

3 算法驗(yàn)證

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

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

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

4 頁(yè)巖氣流動(dòng)模擬

4.1 頁(yè)巖氣滲流場(chǎng)分析

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

圖3 頁(yè)巖掃描電鏡圖

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

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

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

圖5 流場(chǎng)穩(wěn)定時(shí)速度分布場(chǎng)圖

圖6 流場(chǎng)穩(wěn)定時(shí)速度等值線圖

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

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

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

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

5 結(jié)論

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

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

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

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

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

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

4朱光普,姚軍,樊冬艷等.頁(yè)巖氣藏壓裂水平井試井分析.力學(xué)學(xué)報(bào),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é)報(bào)，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ì)觀滲流場(chǎng)的數(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方法的頁(yè)巖氣流動(dòng)數(shù)值模擬.力學(xué)與實(shí)踐,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