劉建鋒, 何 鑫, 薛福軍, 代晶晶, 楊建雄, 黃浩勇, 侯正猛

(1.四川大學(xué) 水利水電學(xué)院, 四川 成都 610065; 2.中國石油天然氣股份有限公司西南油氣田分公司 頁巖氣研究院, 四川 成都 610051; 3.四川大學(xué) 中德能源研究中心, 四川 成都 610065)

頁巖氣作為一種清潔的非常規(guī)油氣資源, 在全球非常規(guī)油氣的開發(fā)進(jìn)程中, 占據(jù)著重要地位。近年來, 北美頁巖氣的商業(yè)化開采使得世界能源格局發(fā)生變化[1–3], 我國的頁巖氣勘探開發(fā)也處于迅速發(fā)展階段。在當(dāng)前的頁巖氣開采實踐中, 常采用水力壓裂技術(shù)來刺激儲層裂縫[4–6], 從而形成水力壓裂區(qū)域與天然裂縫共同組成的改造儲層體積, 即Stimulate Reservoir Volume (SRV), 從而達(dá)到有效開采頁巖氣的目的。地質(zhì)勘探調(diào)查發(fā)現(xiàn), 頁巖儲層內(nèi)存在豐富的天然裂縫[7–9], 其對儲層滲透性有著不可忽視的影響[10]。通常情況下, 裂縫網(wǎng)絡(luò)的發(fā)育程度越高, 儲層的滲透性就越高, 儲層的開發(fā)效果也越佳[11–14], 而不同的天然裂縫特征, 會對儲層的滲透性產(chǎn)生較大影響, 從而影響頁巖儲層滲流行為及產(chǎn)氣量。因此, 綜合考慮多種天然裂縫特征對儲層滲透性的影響, 對于頁巖氣的開發(fā)利用具有重要意義。

近年來, 國內(nèi)外學(xué)者在天然裂縫特征對儲層滲流行為的影響方面進(jìn)行了諸多研究。李玉梅等[15]基于滲流–應(yīng)力耦合數(shù)值算法, 分析了壓裂裂縫長度、天然裂縫角度等因素對儲層滲流的影響; 趙猛等[16]基于天然裂縫網(wǎng)絡(luò)模型, 探究了天然裂縫間距、簇數(shù)、滲透率等裂縫特征對儲層滲流的影響,并對不同特征的影響大小進(jìn)行了排序; 杜旭林等[17]基于嵌入式離散裂縫模型, 并考慮了裂縫的閉合作用, 以此研究了天然裂縫寬度特征對儲層滲流的影響; 慎國強等[18]基于流固耦合作用下的雙重介質(zhì)滲流理論, 分析了天然裂縫長度、傾斜角以及連通性對儲層滲流的影響; 張志偉等[19]通過建立等效裂縫滲流模型, 探究了裂縫長度、寬度、滲透率對儲層滲流量的影響。雖然已有研究在天然裂縫特征影響儲層滲流方面已有較多具有指導(dǎo)意義的成果, 但以上成果尚未充分探討不同天然裂縫特征組合對儲層滲流行為的影響, 且鮮少涉及對于天然裂縫密集度特征的研究, 會影響頁巖氣藏水力壓裂施工方案設(shè)計、產(chǎn)量評估等。為解決以上問題, 筆者從天然裂縫的多種特征著手, 系統(tǒng)探究了天然裂縫特征組合變化對儲層滲流的綜合影響。

數(shù)值模擬在頁巖氣的開發(fā)研究中具有十分重要的意義, 通過數(shù)值模擬的方式, 能夠較為準(zhǔn)確地模擬頁巖儲層內(nèi)頁巖氣的賦存狀態(tài)、天然裂縫分布特征及縫網(wǎng)形態(tài)等, 從而對儲層產(chǎn)量進(jìn)行有效地分析評估, 這種方法使得研究者能夠更好地了解頁巖氣開采的相關(guān)條件, 為提高頁巖儲層產(chǎn)量提供有益參考[14,20–23]。COMSOL Multiphysics為一種基于有限元方法的數(shù)值模擬軟件, 具有很強的多物理場耦合計算能力, 包含廣泛的物理庫和多模塊功能,可自定義偏微分方程, 支持并行計算, 對于解決多場耦合下的基質(zhì)–縫網(wǎng)滲流問題有著很高的適用性。

筆者基于COMSOL軟件對特定頁巖儲層進(jìn)行了數(shù)值模擬, 旨在研究不同天然裂縫特征組合對儲層滲流的影響規(guī)律, 以期通過該研究對頁巖儲層壓裂施工方案設(shè)計、產(chǎn)量提升提供參考。

1 數(shù)學(xué)模型

在頁巖氣的開采過程中, 孔隙壓力的變化會導(dǎo)致儲層產(chǎn)生變形, 使得儲層的滲流參數(shù)發(fā)生改變,從而影響儲層的滲流量, 且儲層中的氣體流動又包括了基質(zhì)中的氣體擴(kuò)散以及縫網(wǎng)中的氣體滲流, 故而筆者基于此建立數(shù)學(xué)模型。

1.1 滲流場方程

頁巖氣在儲層中主要以吸附氣和游離氣的形式存在, 其滲流形式可分為在巖石基質(zhì)孔隙中的擴(kuò)散以及在裂縫中的流動。在基質(zhì)系統(tǒng)中, 頁巖氣的擴(kuò)散可以通過達(dá)西定律進(jìn)行描述, 頁巖氣的吸附特性可以通過朗格繆爾的吸附方程來考慮, 結(jié)合小變形原理及質(zhì)量守恒, 整理后可得到基質(zhì)系統(tǒng)內(nèi)的滲流方程:

式中,φm為基質(zhì)的孔隙度;km為基質(zhì)滲透率,10–3μm2;μg為氣體黏度, Pa·s;αm為比奧系數(shù);Vε為體積應(yīng)變;ρg為氣體密度, kg/m3;M為氣體摩爾質(zhì)量, kg/mol;R為普適氣體常數(shù), J/(mol·K);T為溫度, K;LV為朗格繆爾體積常數(shù), m3/kg;PL為朗格繆爾壓力常數(shù), Pa;ρs為基質(zhì)密度, kg/m3;ρgst為常態(tài)下的氣體密度, kg/m3;pmg為基質(zhì)內(nèi)的氣體壓力, Pa。

在裂縫系統(tǒng)中, 水力裂縫和天然裂縫是頁巖氣流動的主要通道, 可用達(dá)西定律結(jié)合質(zhì)量守恒定律描述其流動, 裂縫內(nèi)的滲流方程可以表示為

式中,kf為裂縫滲透率, 10–3μm2;pfg為裂縫內(nèi)的氣體壓力, Pa;ω為裂縫張開度, 指裂縫面在法向上的位移, m。

1.2 變形場方程

在頁巖儲層系統(tǒng)中, 變形場方程會受到孔隙壓力的影響[24–27], 忽略慣性力和黏性力的作用, 根據(jù)Biot理論, 其平衡方程為

式中,σ為總應(yīng)力, Pa;g為重力加速度, m/s2。

根據(jù)巖體小變形假設(shè), 儲層內(nèi)巖石基質(zhì)的應(yīng)力–應(yīng)變關(guān)系可以表示為

式中,σ′為基質(zhì)的有效應(yīng)力, Pa;C為彈性剛度張量;ε為應(yīng)變張量;I為單位張量。

將式(4), (5)代入式(3), 可得:

式中,u為位移分量。

為了表述儲層內(nèi)裂縫面的應(yīng)力–應(yīng)變關(guān)系, 引入式(7)[28–30]。

式中,σf′為裂縫面有效應(yīng)力, Pa;T為裂縫面剛度矩陣。

裂縫面不連續(xù)處加載的自然邊界條件可表示為

式中,n為裂縫表面的單位向量。

2 模型驗證

2.1 DFM方法驗證

為了驗證該數(shù)學(xué)模型的準(zhǔn)確性, 通過引入離散裂縫模型DFM(Discrete Fracture Model)進(jìn)行數(shù)值模擬, 將其結(jié)果與擴(kuò)展有限元法XFEM(Extended Finite Element Method)及擴(kuò)展有限體積法XFVM (Extended Finite Volume Method)的結(jié)果[31]進(jìn)行對比,以此來驗證模型的正確性。

離散裂縫模型的網(wǎng)格如圖1所示。

圖1 DFM網(wǎng)格及邊界條件Fig.1 DFM Mesh and boundary conditions

包含6 850個域單元和224個邊界單元, 模型高寬均為100 m, 采用自由三角形網(wǎng)格劃分, 模型正中央有一水平裂縫, 長度為40 m。模型左、右邊界上的壓力分別為1×107Pa與1×105Pa, 上、下兩側(cè)采用無流動邊界條件, 模型區(qū)域內(nèi)初始壓力為0。四周均采用固定約束邊界條件。模擬使用的數(shù)據(jù)見表1。

表1 數(shù)值模擬的主要參數(shù)Table 1 Main parameters for numerical simulation

使用COMSOL軟件引入上述數(shù)學(xué)模型, 并對DFM模型進(jìn)行數(shù)值模擬, 得到流體壓力分布, 如圖2所示。

圖2 流體壓力分布Fig.2 Fluid pressure distribution

由圖2可以明顯地看出, 由于邊界條件壓力差的緣故, 基質(zhì)部分表現(xiàn)出明顯的壓力梯度; 裂縫部分, 由于裂縫的滲透率相對于基質(zhì)更大, 因此裂縫內(nèi)部的流動速度更快, 表現(xiàn)出沿著裂縫方向的優(yōu)勢流動, 這說明數(shù)值模擬結(jié)果與實際情況相吻合。圖3為x及y方向上的位移場分布情況, 由圖3可知, 位移場分布大致沿裂縫對稱, 從y方向上的位移場分布可見, 裂縫上表面有著向下的位移趨勢, 裂縫下表面有著向上的位移趨勢, 因此得到由于裂縫內(nèi)外壓力差的緣故, 裂縫表現(xiàn)為趨向于閉合的趨勢。

圖3 x, y方向位移場云圖Fig.3 Displacement field maps in x, y directions

為了將不同裂縫位置的數(shù)值模擬結(jié)果與已有文獻(xiàn)結(jié)果進(jìn)行對比, 讀取裂縫上距離裂縫最左側(cè)點0, 5, 10, 15, 20, 25, 30, 35, 40 m共9個不同裂縫位置的上表面y方向位移大小, 由此得到數(shù)值模擬的結(jié)果, 并將其與已有文獻(xiàn)結(jié)果進(jìn)行對比, 如圖4所示。

圖4 模擬與已有文獻(xiàn)結(jié)果對比Fig.4 Comparison of simulation and literature results

由圖4可知, 數(shù)值模擬與文獻(xiàn)[31]的結(jié)果吻合較好, 由此可以驗證數(shù)學(xué)模型的正確性。

2.2 現(xiàn)場驗證

為了進(jìn)一步驗證該數(shù)學(xué)模型的正確性, 以美國密西西比Barnett頁巖氣藏現(xiàn)場生產(chǎn)資料[6]進(jìn)行驗證, 此處的頁巖氣資源豐富, 如圖5所示, 其開采經(jīng)歷了長達(dá)幾年的時間并且有系統(tǒng)的監(jiān)測數(shù)據(jù), 便于開展模型驗證。

圖5 Barnett頁巖氣藏研究區(qū)域與地質(zhì)概況Fig.5 Study area of Barnett shale and the geological characteristics

該頁巖氣藏主要通過水力壓裂、水平井開發(fā)的方式進(jìn)行開采, 為便于計算, 筆者將整個Barnett頁巖儲層中的單個水力壓裂區(qū)域建立模型進(jìn)行數(shù)值仿真, 由此得到滲流量, 并與現(xiàn)場實際生產(chǎn)數(shù)據(jù)進(jìn)行對比, 以驗證模型的正確性。

所使用的數(shù)值模擬驗證模型及邊界條件如圖6所示。Barnett頁巖儲層總尺寸為1 100 m×290 m,厚90 m, 筆者選擇其中一條水力壓裂裂縫及周邊區(qū)域作為模擬模型, 區(qū)域總數(shù)量56個, 故模型尺寸選擇為30.5 m×145 m的單個水力壓裂區(qū)域, 厚90 m,采用自由三角形網(wǎng)格劃分; 水力壓裂裂縫長度為47.2 m, 壓裂間距為30.5 m, 所使用的儲層參數(shù)來源于美國New Mexico東南部的Barnett頁巖[6,14], 見表2。

表2 數(shù)值模擬主要參數(shù)Table 2 Main parameters for numerical simulation of shale reservoirs

圖6 數(shù)值模擬驗證模型及邊界條件Fig.6 Numerical simulation validation model and boundary condition diagram

基于COMSOL軟件對頁巖氣儲層模型進(jìn)行數(shù)值模擬, 通過計算得到Barnett頁巖儲層內(nèi)不同時間裂縫內(nèi)的氣體壓力分布如圖7所示。由圖7可知, 由于水平井與頁巖儲層存在較大的壓力差, 初始階段氣體會通過水力裂縫快速流向水平井(圖7(a)), 待水力裂縫內(nèi)的氣壓與水平井基本一致后, 水力壓裂區(qū)內(nèi)的氣體將逐漸流入水力裂縫, 引起壓裂區(qū)域氣壓下降, 隨著時間的推移, 離壓裂區(qū)較遠(yuǎn)處基質(zhì)內(nèi)的氣體會逐漸流入水力壓裂區(qū), 進(jìn)而流入水平井,這一過程逐漸由水平井朝遠(yuǎn)離水平井的區(qū)域發(fā)生擴(kuò)散(圖7(b)), 當(dāng)基質(zhì)內(nèi)的氣壓接近開采壓力時, 氣體滲流過程會逐漸緩慢, 其氣壓分布如圖7(c)所示。

圖7 裂縫內(nèi)的氣體壓力分布Fig.7 Gas pressure distribution inside the crack

頁巖氣藏的數(shù)值模擬與實測結(jié)果對比驗證如圖8所示。由圖8可知, 模擬結(jié)果峰值與現(xiàn)場生產(chǎn)數(shù)據(jù)峰值幾乎一致。10~300 d內(nèi), 滲流量從2.18×105m3/d降低至5.39×104m3/d, 模擬結(jié)果與現(xiàn)場數(shù)據(jù)吻合較好; 300~1 600 d內(nèi), 滲流量從5.39×104m3/d降低至1.84×104m3/d, 模擬結(jié)果略高于現(xiàn)場數(shù)據(jù),大致吻合, 推測可能是由于基質(zhì)內(nèi)水分的影響。模擬結(jié)果與現(xiàn)場實際生產(chǎn)數(shù)據(jù)整體吻合較好, 由此可驗證模型的正確性。

圖8 模擬與實測結(jié)果對比Fig.8 Comparison of simulated and measured results

3 不同天然裂縫特征組合的影響

在頁巖氣實際生產(chǎn)過程中, 不同的儲層條件相差各異, 其天然裂縫特征情況也各不相同[32], 不同的天然裂縫滲透率、密集度、傾斜角等條件都會對儲層的滲流產(chǎn)生一定的影響。為研究不同天然裂縫特征變化對于儲層滲流的影響, 使用COMSOL軟件引入前文已經(jīng)驗證的數(shù)學(xué)模型, 通過生成隨機裂縫, 有梯度地改變特征因素進(jìn)行研究, 研究中將特征因素兩兩組合進(jìn)行數(shù)值模擬, 以探究不同的天然裂縫特征組合對儲層滲流產(chǎn)生的影響。模擬使用的儲層模型尺寸為60 m×140 m的單個水力壓裂區(qū)域, 厚20 m, 采用自由三角形網(wǎng)格劃分; 水力壓裂裂縫長度為80 m, 壓裂間距為60 m, 所使用的其余基本儲層參數(shù)與2.2節(jié)中參數(shù)相同, 見表2。

基本模型及邊界條件如圖9所示, 以長方形區(qū)域代表單個儲層區(qū)域, 長直線段代表水力壓裂裂縫, 其周圍區(qū)域為裂縫區(qū)域, 隨機的線段代表天然裂縫。

圖9 天然隨機裂縫模型及邊界條件Fig.9 Natural random fracture model and boundary condition diagram

3.1 滲透率及密集度的影響

為了研究天然裂縫滲透率與密集度組合變化對頁巖儲層滲流的影響, 以50×10–3μm2為梯度設(shè)置了(50~500)×10–3μm2共10組不同的天然裂縫滲透率, 并用頁巖儲層區(qū)域單位面積內(nèi)的天然裂縫總長度來表征密集度, 以0.08 m/m2為梯度設(shè)置了0~0.24 m/m2共4組不同的密集度, 通過數(shù)值模擬軟件進(jìn)行分別計算, 得到不同天然裂縫滲透率、密集度情況下的滲流量峰值, 并繪制3D曲面如圖10所示。由圖10可知, 該3D曲面沿著密集度的傾斜度較沿著滲透率的更大, 這說明天然裂縫密集度變化對頁巖儲層滲流的影響更為顯著。此外, 由圖10還可知, 滲流量峰值隨著滲透率及密集度的增加, 始終保持著增加趨勢, 但變化率卻在不同的滲透率及密集度下有所不同, 這也說明在不同滲透率下, 密集度的增加對儲層滲流的影響是不同的, 為了探究這種影響的大致規(guī)律, 選擇滲透率為50×10–3μm2,500×10–3μm2的情況進(jìn)行量化分析, 并且計算了比例, 結(jié)果見表3。由表3可知, 當(dāng)密集度增加時, 滲流量明顯提升, 滲透率為50×10–3μm2時的最大增加比例為3.5%, 500×10–3μm2時的最大增加比例為10.05%。由以上結(jié)果可知, 密集度的增加能夠增加頁巖儲層滲透性, 且隨著滲透率的增大, 密集度對頁巖儲層滲流的影響將會變得更大。

表3 滲透率及密集度組合時不同密集度的滲流量分析Table 3 Seepage volume analysis for different density with permeability and density combinations

圖10 不同裂縫密集度及滲透率情況下的滲流量峰值Fig.10 Surface plot of peak seepage volume for different fracture densities and permeabilities

同理, 在不同密集度情況下, 滲透率對儲層滲流的影響也不相同, 選擇裂縫密集度為0.08,0.24 m/m2時的情況進(jìn)行對比, 結(jié)果見表4。由表4可知, 當(dāng)滲透率提升時, 滲流量明顯增加, 裂縫密集度為0.08 m/m2時的最大增加比例為2.61%, 0.24 m/m2時的最大增加比例為6.32%。由以上結(jié)果可知,滲透率的增加能夠增加頁巖儲層滲透性, 且隨著密集度的增大, 滲透率對頁巖儲層滲流的影響將會增加。

表4 滲透率及密集度組合時不同滲透率的滲流量分析Table 4 Seepage volume analysis for different permeability with permeability and density combinations

3.2 滲透率及傾斜角的影響

生成傾斜角度一致的天然裂縫模型如圖11所示, 以30°為梯度設(shè)置了0°~90°共4組不同的天然裂縫傾斜角, 在不同的滲透率條件下, 通過數(shù)值模擬軟件進(jìn)行分別計算, 得到不同天然裂縫傾斜角、滲透率情況下的滲流量峰值, 并繪制3D曲面如圖12所示, 由圖12可知, 該3D曲面沿著傾斜角方向的變化率明顯較沿著滲透率方向的更大, 這說明裂縫傾斜角對頁巖儲層滲流的影響更大。

圖11 裂縫傾斜角模型Fig.11 Fracture tilt angle model

圖12 不同裂縫滲透率及傾斜角情況下的滲流量峰值Fig.12 Surface plot of peak seepage volume for different fracture permeabilities and tilt angles

與3.1節(jié)類似, 選擇對裂縫滲透率為50×10–3,500×10–3μm2的情況進(jìn)行對比, 裂縫傾斜角為0°,90°的情況進(jìn)行對比, 結(jié)果見表5, 6。由表5, 6數(shù)據(jù)可知, 滲透率為50×10–3μm2時的最小縮減比例為–3.99%, 500×10–3μm2時為–9.06%; 傾斜角為0°時的最大增加比例為5.92%, 90°時為0.33%。

表5 滲透率及傾斜角組合時不同傾斜角的滲流量分析Table 5 Seepage volume analysis for different tilt angle with permeability and tilt angle combinations

表6 滲透率及傾斜角組合時不同滲透率的滲流量分析Table 6 Seepage volume analysis for different permeability with permeability and tilt angle combinations

由以上結(jié)果可知, 傾斜角的增加會降低頁巖儲層滲透性, 并且隨著滲透率的增大, 傾斜角對頁巖儲層滲流的影響將會變得更大; 而隨著傾斜角的增大, 滲透率對頁巖儲層滲流的影響將會變得更小。

3.3 密集度及傾斜角的影響

在不同的密集度及傾斜角條件下, 通過數(shù)值模擬軟件進(jìn)行分別計算, 并繪制3D曲面如圖13所示。由圖13可知, 該3D曲面沿著密集度的變化率與沿著傾斜角的變化率大致近似, 需要通過具體計算進(jìn)一步確定影響的相對大小。與3.1節(jié)類似, 選擇密集度為0.08, 0.24 m/m2, 以及裂縫傾斜角為0°,90°的情況進(jìn)行對比, 結(jié)果見表7, 8。

表7 密集度及傾斜角組合時不同傾斜角的滲流量分析Table 7 Seepage volume analysis for different tilt angle with densities and tilt angle combinations

表8 密集度及傾斜角組合時不同密集度的滲流量分析Table 8 Seepage volume analysis for different densities with densities and tilt angle combinations

圖13 不同裂縫密集度及傾斜角情況下的滲流量峰值Fig.13 Surface plot of peak seepage volume for different fracture densities and tilt angles

由表7, 8可知, 密集度為0.08 m/m2時的最小縮減比例為–3.65%, 0.24 m/m2時為–4.76%; 傾斜角為0°時的最大增加比例為7.35%, 90°時為2.24%。

對比以上數(shù)據(jù)可知, 隨著密集度的增大, 傾斜角對頁巖儲層滲流的影響將會變得更大; 而隨著傾斜角的增大, 密集度對頁巖儲層滲流的影響將會變得更小。傾斜角與密集度的變化對頁巖儲層滲流的影響大致近似, 密集度的影響略大一些。

4 結(jié) 論

(1)根據(jù)頁巖儲層水力壓裂后的頁巖氣流動特征, 以質(zhì)量平衡方程為基礎(chǔ), 推導(dǎo)了儲層內(nèi)的滲流場及應(yīng)力場數(shù)學(xué)模型, 通過DFM模型模擬裂縫–基質(zhì)內(nèi)的流固耦合行為, 將數(shù)值模擬結(jié)果與已有文獻(xiàn)結(jié)果進(jìn)行對比, 驗證了模型的正確性。

(2)天然裂縫滲透率、密集度、傾斜角均對儲層滲透性有較明顯的影響, 增加天然裂縫滲透率會加大傾斜角、密集度對儲層滲流的影響, 對產(chǎn)氣量的影響約10%; 增加密集度會加大滲透率對儲層滲流的影響, 減弱傾斜角對儲層滲流的影響, 對產(chǎn)氣量的影響約6.3%; 增加傾斜角會同時減弱密集度、滲透率對儲層滲流的影響, 對產(chǎn)氣量的影響約7.4%。

(3)增加天然裂縫傾角導(dǎo)致的儲層滲透性降低現(xiàn)象可能與天然裂縫、水力裂縫交叉數(shù)量降低相關(guān), 所建立的數(shù)值模型能很好捕捉并說明這一宏觀現(xiàn)象, 對現(xiàn)場水力壓裂方案及參數(shù)的確定、提高頁巖氣的持續(xù)高產(chǎn)有重要指導(dǎo)意義。