頁巖氣藏體積壓裂水平井產(chǎn)能有限元數(shù)值模擬

何易東，任嵐，趙金洲，李志強，鄧鵬

（1.西南石油大學油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，四川 成都 610500；2.中國石油塔里木油田分公司塔中油氣開發(fā)部，新疆 庫爾勒 841000）

頁巖氣藏體積壓裂水平井產(chǎn)能有限元數(shù)值模擬

何易東1，任嵐1，趙金洲1，李志強1，鄧鵬2

（1.西南石油大學油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，四川 成都 610500；2.中國石油塔里木油田分公司塔中油氣開發(fā)部，新疆 庫爾勒 841000）

考慮到壓裂過程中的多重復合作用，將壓后頁巖儲層分為支撐主裂縫、縫網(wǎng)波及區(qū)與未壓裂區(qū)?？紤]基巖納米孔隙中氣體吸附與解吸、Knudsen擴散、滑脫流、黏性流，以及水壓誘導裂縫應力敏感效應，建立了頁巖氣藏體積壓裂生產(chǎn)動態(tài)模擬的物理模型和滲流數(shù)學模型。結(jié)合Galerkin有限元方法，對基質(zhì)和裂縫滲流方程進行空間上的離散，推導了三角形單元有限元數(shù)值模型，給出了壓裂水平井二維滲流場內(nèi)、外邊界條件和水力裂縫處理方法，對時間域采用向后差分，最后順序求解裂縫和基質(zhì)壓力方程，模擬了頁巖氣藏體積壓裂水平井生產(chǎn)動態(tài)和壓力場分布。該研究為頁巖氣儲層體積壓裂產(chǎn)能評價提供了理論模型，對于有限元法模擬雙重介質(zhì)滲流場和產(chǎn)能預測具有現(xiàn)實意義。

頁巖氣；體積壓裂；雙重介質(zhì)；有限元；產(chǎn)能模擬

0 引言

近年來，頁巖氣藏的開發(fā)方興未艾，然而頁巖儲層物性較差，基巖滲透率一般在10-9μm2數(shù)量級，孔隙度一般在4%～6%。同時，儲層非均質(zhì)性也較為嚴重，天然裂縫發(fā)育。通過水力壓裂溝通儲層天然裂縫，形成較大范圍的網(wǎng)絡裂縫，是實現(xiàn)該類儲層有效開發(fā)的關鍵技術之一，油藏增產(chǎn)體積是衡量頁巖氣藏改造效果的關鍵參數(shù)［1］。頁巖氣的儲集、運移機理和常規(guī)天然氣存在較大差異，目前認為，頁巖氣主要以游離氣和吸附氣形式儲集于基質(zhì)納米孔隙中，一部分以游離氣形式儲集于天然裂縫系統(tǒng)中，還有一部分是以溶解氣形式儲集于有機質(zhì)干酪根中［2］。頁巖儲層中發(fā)育大量連通的納米孔隙、不連通或閉合的天然微裂縫以及水力壓裂形成的體積裂縫。頁巖氣在開采過程中將在納米孔隙—微裂縫—水壓體積裂縫多尺度空間進行流動傳輸。首先從納米孔隙介質(zhì)向微裂縫發(fā)生微觀分子級的擴散流動匯聚，然后從微裂縫向宏觀水力壓裂體積裂縫發(fā)生宏觀滲流［3］。頁巖氣開發(fā)過程中，流體傳輸機理和流動路徑復雜，流態(tài)形式多樣［4］，但目前對頁巖氣在儲層中滲流機理認識還不清，頁巖氣儲層體積壓裂產(chǎn)能模擬方法也存在嚴重滯后。

應用有限元法可以較好地處理內(nèi)外邊界條件，且其網(wǎng)格劃分的靈活性優(yōu)于有限差分法。目前有限元法已經(jīng)成功應用于模擬油氣藏單相和多相流體流動［5-10］，但是對于模擬頁巖這種雙重介質(zhì)氣藏的研究較少［11］，特別是如何采用Galerkin有限元法模擬頁巖氣在多尺度空間中的復雜流動，建立相應的有限元數(shù)值求解模型，以及頁巖儲層壓裂水平井考慮儲層中流體為二維流動的內(nèi)外邊界條件處理方法尤為少見。

1 物理模型

對于天然裂縫發(fā)育、低水平主應力差的頁巖儲層，由于低黏度壓裂液的攜砂能力差，在壓裂施工過程中形成一條高導流能力的支撐主裂縫，主裂縫之間出現(xiàn)相互連通的誘導裂縫網(wǎng)絡［12］，從而形成一個體積壓裂區(qū)。如圖1所示，該復合氣藏模型由支撐主裂縫、體積壓裂區(qū)和未壓裂區(qū)構(gòu)成［13］。為便于工程上的模擬應用，該復合模型可以用離散裂縫描述支撐主裂縫，而采用雙重介質(zhì)模型來描述體積壓裂區(qū)和未壓裂區(qū)［14］，也可以將未壓裂區(qū)考慮為單孔隙介質(zhì)系統(tǒng)。

2 頁巖氣滲流數(shù)學模型

2.1 連續(xù)性方程

隨著基質(zhì)納米孔隙半徑的減小，氣體分子平均自由程和流動空間的特征長度將越來越接近，此時，氣體分子和孔壁的碰撞將會更加頻繁，氣體分子更容易表現(xiàn)出微觀級的分子流態(tài)，這種效應使得儲層表觀滲透率進一步增加［15］?？紤]納米孔隙中氣體滑脫、Knudsen擴散和黏性流的表觀滲透率模型為

圖1 頁巖儲層體積壓裂產(chǎn)能模擬物理模型

頁巖儲層水力壓裂后，氣體將在由基質(zhì)納米孔隙、天然裂縫和人工裂縫組成的多尺度空間中運移，其流動機理復雜，主要包括基質(zhì)中的吸附、解吸、Knudsen擴散、黏性流，以及裂縫中的黏性流［16-17］。由此，頁巖基質(zhì)系統(tǒng)中氣體連續(xù)性方程為

頁巖氣在基質(zhì)納米孔隙中的吸附量可由Langmuir等溫吸附方程描述：

頁巖氣在天然裂縫系統(tǒng)中的連續(xù)性方程為

頁巖儲層水力壓裂后，主裂縫周圍的縫網(wǎng)一般處于未支撐或弱支撐狀態(tài)，隨著儲層孔隙壓力的下降，開啟的裂縫將會逐漸閉合，表現(xiàn)出較強的應力敏感性。表征體積壓裂區(qū)滲透率隨壓力變化的公式為

頁巖氣在主裂縫中的連續(xù)性方程為

2.2 初始條件

初始條件下，基質(zhì)系統(tǒng)和裂縫系統(tǒng)壓力相同，均為原始儲層壓力：

2.3 邊界條件

Γ1表示求解域外邊界，Γ2表示求解域內(nèi)邊界。假設對于給定的研究區(qū)域外邊界封閉，內(nèi)邊界定壓，則邊界條件為

3 數(shù)值離散模型

3.1 Galerkin有限元方法

將式（3）代入式（2）中，只考慮儲層中流體平面二維流動，將式（2）和式（4）展開得：

其中：Cft=Cg+Cf，Cmt=Cg+Cm。

對于頁巖氣滲流微分方程（9），由Galerkin有限元法在單元內(nèi)進行積分得［18-19］：

應用Green-Gauss公式得：

采用Galerkin加權(quán)殘量法來推導式（14）的有限元計算格式，采用三角形單元進行剖分，設插值函數(shù)為

將式（15），（16）代入式（14）并寫成矩陣形式，得：

則式（19）變?yōu)?/p>

類似于裂縫系統(tǒng)有限元方程的推導過程，可得基質(zhì)系統(tǒng)有限元方程：

上述單元矩陣方程Mf和Mm中，氣體密度、黏度、以及基質(zhì)和裂縫滲透率均是壓力的函數(shù)，考慮單元內(nèi)物性參數(shù)均勻。因此，可以用三角形單元形心處的值近似處理，即：

暫時不考慮方程（20）左邊最后一項（即壓力場邊界條件影響項，將在后面進行討論）。對時間域采用向后差分格式離散，并對與壓力有關的系數(shù)進行顯式處理［20］，最終得到天然裂縫系統(tǒng)任意單元矩陣平衡方程：

對單元有限元方程進行總體合成，獲得整體的有限元方程，并引入初始條件和邊界條件，在有限元整體平衡方程中，由于互含基質(zhì)和天然裂縫壓力項，因此，先求解裂縫系統(tǒng)壓力方程，再求解基質(zhì)系統(tǒng)壓力方程，從而可獲得地層中任一節(jié)點處壓力隨時間的變化規(guī)律，并求得任意時間步井底流量。

3.2 邊界條件處理

對于方程（20）中邊界影響項?！遥∟）TeρgvdΓ，若單元e為不含邊界的單元，則該項自然為0，單元特征分析后處理為［0，0，0］T向量。值得注意的是，對于封閉外邊界，v=0，該項依然處理為［0，0，0］T向量。因此，對于內(nèi)部單元和封閉外邊界，該項自然為0。對于內(nèi)邊界，由于井筒尺寸相比整個油藏小得多，因此，根據(jù)質(zhì)量守恒原理，可以將qgwell考慮為裂縫單元中的點匯項，這與采用有限差分法進行油藏模擬類似。將裂縫放寬處理并剖分為多個矩形單元，進一步將矩形單元剖分為三角形單元后，井點則位于三角形單元的某一節(jié)點上。在實際計算中，采用上一時間步的井點壓力計算的點源或點匯強度，代入總剛度矩陣合成的某一裂縫單元節(jié)點上，對于定產(chǎn)內(nèi)邊界，可以直接在單元節(jié)點上加上或減去一個源匯強度。對于定壓內(nèi)邊界：

4 實例模擬

4.1 基本參數(shù)

頁巖儲層厚度40 m，長度1 500 m，寬度400 m，井筒長度1 200 m，井底流壓10 MPa，原始地層壓力20 MPa，地層溫度70℃，基質(zhì)孔隙度5%，天然裂縫孔隙度0.1%，基質(zhì)孔隙半徑5 nm，孔隙介質(zhì)迂曲度為5，天然裂縫滲透率 3×10-5μm2，水力主裂縫導流系數(shù) 0.5 μm2·cm，次裂縫導流系數(shù) 0.03 μm2·cm，縫網(wǎng)長 120 m，縫網(wǎng)寬90 m，主裂縫數(shù)10條，天然裂縫間距3 m，天然裂縫應力敏感系數(shù)0.05 MPa-1，基質(zhì)壓縮系數(shù)2×10-4MPa-1，天然裂縫壓縮系數(shù) 0.05 MPa-1，Langmuir體積 2×10-3m3/kg，Langmuir壓力 5 MPa， 頁巖密度 2 600 kg/m3，甲烷分子質(zhì)量16 g/mol，標準狀況下頁巖氣摩爾體積0.022 4 m3/mol。

4.2 模擬結(jié)果分析

基于上述參數(shù)，模擬頁巖氣水平井體積壓裂后生產(chǎn)5，10 a的儲層壓力分布（見圖2），可以看出體積壓裂區(qū)內(nèi)的儲層壓力下降較多，而未壓裂區(qū)的儲層壓力幾乎沒有變化，表明該時間段主要采出了改造體積內(nèi)的游離氣和吸附氣。

圖2 頁巖氣生產(chǎn)過程中儲層壓力分布

圖3，4為頁巖儲層改造體積對產(chǎn)氣量的影響情況。由圖可以看出：在其他參數(shù)相同的情況下，改造體積SRV從36×104m3增大到180×104m3， 產(chǎn)氣量增加幅度較大；隨著改造體積的增加，日產(chǎn)氣量和累計產(chǎn)氣量都增加，但是增加幅度逐漸變小，對于給定儲層地質(zhì)條件，存在最佳增產(chǎn)改造體積。

圖3 頁巖儲層改造體積對日產(chǎn)氣量的影響

圖4 頁巖儲層改造體積對累計產(chǎn)氣量的影響

圖5是本文有限元方法與商業(yè)軟件模擬產(chǎn)氣量結(jié)果的對比（SRV=180×104m3）。由圖可以看出，兩者結(jié)果誤差為8%，這也說明了有限元方法模擬結(jié)果的準確性，同時商業(yè)軟件的模擬結(jié)果略小于有限元方法的模擬結(jié)果，這是因為本文有限元方法考慮了表觀滲透率的影響。

圖5 模擬結(jié)果對比

5 結(jié)論

1）基于三孔雙滲連續(xù)介質(zhì)模型，考慮基巖納米孔隙中氣體吸附解吸效應、Knudsen擴散和滑脫流、儲層應力敏感效應，建立了頁巖儲層考慮復雜流動行為的動態(tài)評價數(shù)學模型?；陧搸r儲層體積壓裂縫網(wǎng)擴展特征，劃分了支撐主裂縫、縫網(wǎng)波及區(qū)和未改造區(qū)，提出了采用雙重介質(zhì)和離散裂縫模型描述復雜縫網(wǎng)的物理模擬模型。

2）推導了基于Galerkin有限元法離散的頁巖氣雙重介質(zhì)二維滲流數(shù)值模型，給出了壓裂水平井內(nèi)外邊界處理方法，實現(xiàn)了應用有限元法求解頁巖氣壓裂井產(chǎn)能動態(tài)和壓力場分布，對于采用有限元法模擬雙重介質(zhì)氣藏滲流場具有重要的理論意義。

3）頁巖儲層改造體積是影響頁巖氣產(chǎn)量的重要因素，隨著儲層改造體積的增大，氣井產(chǎn)量逐漸增加，但是增加的幅度會下降。開展頁巖氣在多尺度儲層中的復雜流動行為的研究，建立相應的產(chǎn)能評價模型，分析影響產(chǎn)能的縫網(wǎng)參數(shù)和儲層參數(shù)對于認識完善頁巖氣滲流理論，弄清影響產(chǎn)能的主次因素，提高頁巖儲層體積壓裂設計的科學性具有重要的理論和現(xiàn)實意義。

6 符號注釋

Kapp為基巖表觀滲透率，10-3μm2；φ為孔隙度；τ為多孔介質(zhì)迂曲度；rn為頁巖平均基質(zhì)孔隙半徑，nm；Mg為氣體分子質(zhì)量，kg/mol；μ為黏度 mPa·s；R 為絕對氣體常數(shù)，8.314 J·mol/K；T 為絕對溫度，K；ρ為密度，kg/m3；p為壓力，MPa；α為切向動量適應系數(shù)，本文取值0.8；qmf為氣體從基質(zhì)向天然裂縫的流量，kg/（m3·s）；qa為單位體積的吸附量，kg/m3；ρs為頁巖密度，kg/m3；Vstd為氣體在標準狀況下的摩爾體積，m3/mol；VL為Langmuir體積，m3/kg；pL為 Langmuir壓力，MPa；qwell為氣體從裂縫流向井筒的竄流量，kg/（m3·s）；Lx，Ly分別為基質(zhì)巖塊在x，y方向的長度，m；σ為形狀因子；Kfi為裂縫初始條件下的滲透率，10-3μm2；df為天然裂縫應力敏感系數(shù)，MPa-1；pi為原始儲層壓力，MPa；qhf為氣體從主裂縫向天然裂縫的竄流量，kg/（m3·s）；n 為求解域外邊界曲面的外法線矢量；pwf為井底流壓，MPa；Cft，Cmt分別為裂縫、 基質(zhì)的總壓縮系數(shù)，MPa-1；Ni為單元形函數(shù)；Γ為區(qū)域Ω的邊界；QT為水平井總?cè)债a(chǎn)量，m3/d；W為裂縫寬度，m；req為等效井半徑，m；rw為井半徑，m；Bg為體積系數(shù)。下標 m，f，g，hf分別表示基質(zhì)、裂縫、氣體、主裂縫。

［1］ MAYERHOFER M J，LOLONE P，WARPINSKI N R，et al.What is stimulated reservoir volume （SRV）? ［J］.SPE Prod&Oper，2008，15（4）：473-485.

［2］ SWAMI V，SETTARI A T，TAVADPOUR F.A numerical model for multi-mechanism flow in shale gas reservoirs with application to laboratoryscaletesting［C］//TheEAGEAnnualConference&Exhibition incorporating SPE Europe，London，United Kingdom，2013.

［3］ CUO C H，WEI M Z，CHEN H W，et al.Improved numerical simulation for shale gas reservoirs［C］//The Offshore Technology Conference，Kuala Lumpur，Malaysia，2014.

［4］ 任嵐，舒亮，胡永全，等.納米尺度頁巖儲層的氣體流動行為分析［J］.西南石油大學學報（自然科學版），2014，36（5）：112-115.

［5］ HATTIMGH S K F，REDDY B D.A finite element approach for modelling single-phase compressible flow in dual porosity systems ［J］.Journal of Petroleum Science and Engineering，2009，69（1）：1-24.

［6］ PARK Y C，SUNG W M.Development of FEM reservoir modelequipped with effective permeability tensor and its application to naturally fractured reservoirs［C］//The International Oil and Gas Conference and Exhibition in China，Beijing，2000.

［7］ SHENG M，LI G，SHAH S N，et al.Extended finite element modeling of multi-scale flow in fractured shale gas reservoirs ［R］.SPE 159919，2012.

［8］ BHATIA K S，ADVANI S H，LEE J K.Finite element representation of two-phase fluid flow through a naturally fractured reservoir［R］.SPE 19069，1989.

［9］ 程林松，李春蘭，郎兆新.裂縫性底水油藏水平井三維油水兩相有限元數(shù)值模擬方法［J］.石油勘探與開發(fā)，1998，25（2）：41-45.

［10］盛茂，李根生，黃中偉，等.頁巖氣藏流固耦合滲流模型及有限元求解［J］.巖石力學與工程學報，2013，32（9）：1894-1900.

［11］徐兵祥，李相方，張磊，等.裂縫性頁巖氣藏水平井產(chǎn)能預測模型［J］.中國石油大學學報（自然科學版），2013，37（6）：92-99.

［12］ ZHAO Y L，ZHANG L H，LUO J X，et al.Performance of fractured horizontal well with stimulated reservoir volume in unconventional gas reservoir［J］.Journal of Hydrology，2014，512：447-456.

［13］ WU Y S，LI J F，DING D Y，et al.A generalized framework model for the simulation of gas production in unconventional gas reservoirs［J］.SPE Journal，2014，19（5）：845-857.

［14］ SWAMI V，CLARKSON C R，SETTARI A.Non darcy flow in shale Nan pores：do we have a final answer?［R］.SPE 162665，2012.

［15］ AZOM P N，JAVADPOUR F.Dual-continuum modeling of shale and tight gas reservoirs［R］.SPE 159584，2012.

［16］ SHABRO V，TORRES-ERDIN C，SEPEHRNOORI K.Forecasting gas production in organic shale with the combined numerical simulation of gas diffusion in Kerogen，Langmuir desorption from Kerogen surfaces，and advection in Nan pores［R］.SPE 159250，2012.

［17］董平川，牛彥良，李莉.各向異性油藏滲流的有限元數(shù)值模擬［J］.巖石力學與工程學報，2007，26（增刊1）：2634-2639.

［18］劉振宇.有限元法在油藏滲流中的理論和應用［D］.大慶：東北石油大學，2003.

［19］ 姚軍，孫海，樊東艷，等.頁巖氣藏運移機制及數(shù)值模擬［J］.中國石油大學學報（自然科學版），2013，37（1）：91-98.

［20］黃朝琴，姚軍，王月英，等.基于離散裂縫模型的裂縫性油藏注水開發(fā)數(shù)值模擬［J］.計算物理，2011，28（1）：41-49.

（編輯 孫薇）

Finite element numerical simulation of shale gas production of hydraulically fractured horizontal well with stimulated reservoir volume

HE Yidong1,REN Lan1,ZHAO Jinzhou1,LI Zhiqiang1,DENG Peng2
(1.State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation,Southwest Petroleum University,Chengdu 610500,China;2.Tazhong Oil and Gas Development Department,Tarim Oilfield Company,PetroChina,Korla 841000,China)

Affecting by multiple composite effects,shale gas reservoirs with complex fracture-network were divided into propped primary-fracture,stimulated area and unstimulated area.Considering adsorption and desorption of shale gas in matrix porosity,Knudsen diffusion,slippage effect,viscous flow,stress sensitivity of natural fractures,corresponding physical model for production simulation of volume fracturing and mathematical model of shale gas flow after hydraulic fracturing were established.Galerkin finite element has been used to discretize seepage equations for natural fracture and matrix with only taking two-dimensional flow into consideration.Finite element equation based on triangular element was derived,we also give the method for handling the boundary condition with considering two-dimensional seepage field for fractured horizontal well and main hydraulic fracture;backward difference was employed for time domain discretization,and pressure equation of matrix and fracture systems were solved sequentially to model pressure field and production of shale horizontal well with stimulated reservoir volume.This study provides a theoretical model for the evaluation of the volume fracturing of shale gas reservoirs,which has important theoretical and practical significance for the finite element method to simulate the flow field of the double porosity reservoir.

shale gas;volume fracturing;dual continuum;finite element;production simulation

TE357

A

國家自然科學基金青年科學基金項目“基于多尺度頁巖復雜流動行為的動態(tài)評價模型研究”（51404204）；國家科技重大專項“涪陵頁巖氣開發(fā)示范工程”（2016ZX05060）；中國石油化工股份有限公司項目“裂縫動態(tài)擴展影響因素研究”（31400026-16-ZC0607-0013）

10.6056/dkyqt201704025

2016-12-23；改回日期：2017-05-20。

何易東，男，1993年生，在讀碩士研究生，主要從事油氣田增產(chǎn)改造技術研究工作。E-mail：heyidong1993@163.com。

任嵐，男，1979年生，副教授，2003年畢業(yè)于西南石油學院石油工程專業(yè)，主要從事油氣田增產(chǎn)改造技術研究工作。E-mail：renlanswpu@163.com。

何易東，任嵐，趙金洲，等.頁巖氣藏體積壓裂水平井產(chǎn)能有限元數(shù)值模擬［J］.斷塊油氣田，2017，24（4）：550-556.

HE Yidong，REN Lan，ZHAO Jinzhou，et al.Finite element numerical simulation of shale gas production of hydraulically fractured horizontal well with stimulated reservoir volume［J］.Fault-Block Oil&Gas Field，2017，24（4）：550-556.