三維壓裂縫網不穩定壓力半解析求解方法1)

王志凱 程林松,2) 曹仁義,3) 王 進 賈 品 王選茹

* (中國石油大學(北京)石油工程學院,北京 102249)

? (中國石油天然氣股份有限公司長慶油田分公司,西安 710000)

引言

隨著常規油氣資源儲量不斷降低,全球非常規油氣資源勘探開發進入活躍期,致密油、頁巖油等非常規油氣資源探明儲量不斷上升[1-3].特別是近5年,中國非常規油氣進入了規模發現與開發上產階段,致密、頁巖油氣展現出了巨大的資源潛力[4-7].與常規儲層不同,致密儲層多為源儲一體,儲集層致密,未經壓裂儲層難以形成工業油流.開發中需采用體積壓裂技術將儲層“打碎”,形成“人工油氣藏”[8,9].通常認為壓裂縫沿最大主應力方向延伸[10-12],但受原始地應力及壓裂工藝影響,壓裂階段儲層應力狀態非常復雜,現場微地震數據及耦合應力場壓裂模擬結果表明,體積壓裂后儲層縫網展布復雜,存在不同的傾角、形態及穿透比[13-16].

對于上述儲層,其滲流規律較為復雜,單一線性流難以表征其流動狀態[17-19],樊冬艷等[20]、Al Rbeawi 等[21-22]基于點源函數疊加原理研究了相同傾斜方向下多段壓裂無限導流傾斜縫典型試井曲線特征.賈品等[23-25]、Teng 和Li[26]、萬義釗等[27]通過將裂縫離散為微元的方法實現不同裂縫形態的表征,同時Teng 和Li[28-29]、王蘇冉[30]通過引入空間傾斜面源實現單條傾斜縫壓力、產量瞬態分析.但上述研究所考慮的裂縫只能沿單一方向傾斜,在此條件下,可通過調整坐標系的方法降低模型復雜度實現模型有效求解.但唐慧瑩等[14]研究表明,即使在常規多段壓裂過程中,壓裂縫也會因縫間應力干擾發生裂縫面彎曲,出現裂縫傾斜方向不同的現象.因此為更準確表征水力壓裂縫空間展布,建立可表征任意方向傾斜、任意形態展布的多段壓裂井壓力瞬態分析模型[31-34]是非常必要的.

本文通過離散裂縫網絡表征不同裂縫形態,結合點源函數面積分與裂縫微元局部坐標系,構建了不受空間坐標系約束的封閉盒式儲層任意平面縫源函數,并驗證了其準確性.針對三重積分面源積分函數存在的計算效率較低的問題,提出了應用點源或特殊線源代替面源的方法實現了傾斜裂縫快速求解,探討了此類求解方法的準確性及可行性,基于此模型劃分了單段、多段傾斜壓裂縫流態,分析了裂縫導流能力、裂縫傾角、裂縫高度以及裂縫段間距等因素的影響.

1 多段壓裂井三維縫不穩定滲流模型

1.1 物理模型及模型假設

斜縫,裂縫半縫長為lf/2,m;沿裂縫面延伸高度為hf,m;裂縫寬度為w,m;各裂縫傾斜方向不同,每條裂縫內采用局部坐標 (eεI,eξI)進行定位.

如圖1 所示,均質封閉儲層 (xe,ye,h) 中存在一口體積壓裂大斜度井,沿井筒延伸方向存在M條傾儲層及裂縫孔隙度為 φm,φf;滲透率為km,kf,mD;綜合壓縮系數為Ctm,Ctf,MPa?1.流體黏度為μ,mPa·s;體積系數為B;油井定產qw生產,m3/d;井儲系數為C,m3/MPa;不考慮表皮影響.

圖1 多段壓裂井三維縫物理模型Fig.1 Physical model of multi-fractured well with 3D fractures

模型基本假設如下:

(1) 將滲流階段劃分為基質向裂縫滲流及裂縫向井筒流動兩部分,不考慮基質向井筒的直接滲流;

(2) 不考慮井筒沿程壓力降,忽略滲流過程中的重力影響;

(3) 裂縫有限導流,單條裂縫內孔滲數據保持一致;

(4) 儲層內流體為油單相,無毛管力影響.

1.2 裂縫內部流動模型

水力裂縫形態受儲層物性及壓裂工藝影響[35],常見裂縫形態有矩形縫、橢圓縫、菱形縫等.本文將裂縫離散為矩形微元,實現不同形態裂縫有效表征.如圖2 所示,井筒穿過橢圓縫中心,結合橢圓裂縫形態將其離散為17 個矩形微元,其中中心點微元(編號9)為井筒所在微元.

圖2 橢圓縫離散表征Fig.2 Discrete characterization of elliptical fracture

模型推導與計算過程中,裂縫微元相關參數定義如下:第I條裂縫微元總量為NfI,其中井筒所在微元編號NfwI;裂縫系統微元總量為Nf;第i個微元的長、寬為 Δ εi, Δ ξi,m;定義每個裂縫微元的滲透率為Ki,相應的可確定裂縫微元有限導流能力為Kiwi,通過給定不同的Ki值即可實現裂縫不同有限導流能力求解.

基于裂縫內二維不穩定滲流,建立裂縫內流動綜合控制方程.

考慮井筒與橢圓縫相交于微元9,以該微元為研究對象,采用有限差分方法將公式(1)離散處理

式中, β為單位轉換系數,取值0.0864;為第n時間步基質向第i個微元的供給,m3/d;為第n時間步第i條裂縫向井筒的 供給,m3/d;為第n時間步第i個微元的壓力,MPa;Δd(i,j)為微元i,j中心點連線的距離,m.

為簡化模型推導過程,引入無量綱參數

將無量綱參數代入式(2)

化簡可得

式中,TD(i,j)為i,j兩裂縫微元間的無因次傳導率,αDi為相鄰時間步間的影響.

將式(4)應用到離散裂縫系統所有微元上,可得包含Nf個裂縫微元的流動方程有限差分方程組

式中,矩陣T大小為Nf×Nf,T(i,j) =TD(i,j),表征微元間不穩定滲流影響;INf為Nf×Nf的單位矩陣,表征基質源匯項影響;矩陣0 大小為Nf×M的零矩陣;矩陣a大小為Nf×M,表征井筒源匯項影響;矩陣R1大小為Nf× 1,由第n?1 時間步參數決定

1.3 油藏流動模型

考慮上述離散方法所得第I個微元中心點坐標為 (xDI,yDI,zDI),其局部坐標系由矩形微元相鄰邊所在方向的單位向量 (eεI,eξI)構成,基于中心點坐標及局部坐標系可確定微元內任意點坐標

結合Newman 乘積與封閉邊界無限大平面源經典解[36]可得封閉儲層(x,y,z) 處單一點源以定產q生產時在點(x0,y0,z0)處產生壓力降

將無量綱參數代入式(14)

由疊加原理可知,空間面源可由空間點源沿微元面積分獲得,因此空間面源函數表達式

基于空間面源壓力分布,結合壓力疊加原理及杜哈美原理得第i個微元第n時間步在第j個微元處產生的壓力降

為簡化方程形式,定義式(16)中的積分項為參數

從而

基于式(19)建立基質向裂縫流動矩陣

式中,矩陣B大小為Nf× 2(Nf+M),表征基質向各微元流動的影響;矩陣R2大小為Nf× 1,由n? 1 時間步中壓力及流量參數構成參數求解過程中計算效率較低,因此考慮以下

由式(19)可知,隨著Nf及n的增大,式(18)的調用次數迅速升高,而式(18)中的三重積分會導致兩種求解方法.

(1) 線源代替面源

當多條裂縫法向量共面時,可建立坐標軸(x,y,z)使其中一條坐標軸與裂縫微元局部坐標系的一條向量平行,實現模型簡化,如圖3 所示.

圖3 裂縫微元局部坐標系示意圖Fig.3 Schematic diagram of local coordinate system of fracture elements

以局部坐標軸eε與坐標軸x平行為例,此時xξ,yε,zε均為0,式(12)中yD(ε,ξ),zD(ε,ξ)將轉化成ξ的函數

相應的參數可簡化為雙重積分形式

考慮參數由線源沿微元斜邊積分所得,通過調整ξ 方向上微元數量使單個微元內該方向上線源函數變化較小,引入積分中值定理,以微元中心點(xDI,yDI,zDI)處線源與Δξ 之積近似代替源函數積分實現式(24)的進一步簡化,獲得參數的一重積分形式

(2) 點源代替面源

參照方案一的思路,直接對三重積分形式下參數使用積分中值定理,以微元中心點(xDI,yDI,zDI)處點源與微元面積Δε × Δξ 之積近似代替源函數積分實現式(18)的簡化,該條件下需保證裂縫微元的尺寸在一定范圍內,后面將對其精度進行討論

將上述方案中的代入式(21)可實現簡化求解 條件下油藏流動模型求解矩陣的構建.

1.4 模型耦合及求解

由Peaceman 公式[37]可知井筒所在微元內邊界條件可由穩態徑向流公式處理

基于Hurst 等人[38]的研究,在模型中考慮井儲系數CwD的影響.

同時,井筒無限導流條件下,各裂縫處井筒壓力相等,即

本研究選擇武漢市洪山區光谷地區部分城鄉結合部2006年、2009年、2011年3年的遙感影像(3期遙感影像購買自中國遙感數據網，都為6月份的數據)作為初始數據，空間分辨率分別為2.4 m、2.4 m和2.0 m.為了滿足本研究需要，把研究區用地簡單分為城市用地和非城市用地，具體詳見表1，土地利用類型分類參考參考文獻[1]中的分類標準.

結合式(27)?(29)建立與井筒內流動相關滲流方程組

式中,矩陣C大小為2M× 2(Nf+M),表征了式(27)?(29)中第n時間步壓力及產量的系數項,矩陣R3大小為Nf× 1,由n?1 時間步相關參數構成.

該方程組中存在2(Nf+M)個方程,其中未知量包含Nf個裂縫微元的壓力pfD,Nf個裂縫微元的產量qfD,M條裂縫處的井底壓力pfwD以及M條裂縫向井筒的供給qfwD共計2(Nf+M)個,未知量數與方程數相等 ,方程組可封閉求解.

2 模型驗證分析

2.1 三重積分模型驗證

本模型中,考慮多段壓裂大斜度井壓裂縫為有限導流,且裂縫面沿任意方向展布,現有文獻相關研究較少.為驗證模型準確性,分別驗證有限導流單一傾斜縫與無限導流多段壓裂縫求解的準確性.針對圖4 所示部分貫穿傾斜縫(partially penetrating inclined fracture,PPIF),Teng 和Li[28]通過建立坐標軸1 (Coordinate 1)實現模型點源函數求解簡化.為驗證本文模型三重積分求解的準確性,將坐標軸1 沿z 軸旋轉45°建立新坐標軸2,增加模型復雜度.

圖4 坐標軸轉換示意圖Fig.4 Schematic diagram of coordinate axis transformation

模型數據與Teng 圖4 數據相同,具體無因次參數如下:xeD=yeD= 20,hD= 0.5,θ= 45°,hfD=CfD= 5,Cs= 10,CwD= 0,wD= 1.0 × 10?5,γ= 1.3 × 10?8.

坐標軸2 條件下微元內局部坐標系如下

圖5 本文模型與Teng 模型[28]結果對比Fig.5 Comparison of the result of our model with Teng model[28]

為研究多縫條件下模型求解準確性,選取Al Rbeawi 和Tiab[22]2013年提出的無限大地層多段壓裂傾斜縫圖版(A-1)進行對比驗證,選取裂縫傾角為15°,hD= 1,無因次段間距DD= 1.如圖6 所示.

圖6 本文模型與Al Rbeawi 模型結果對比Fig.6 Comparison of the result of our model with Al Rbeawi model

可以看出,不同時間下兩模型計算結果基本吻合,由于本模型采用封閉邊界,當tD>10 壓力響應到達封閉邊界后的壓力及壓力導數變化與Al Rbeawi模型不同由于該模型采用三重積分進行源函數計算,求解時間較長(數小時),適用性較差,考慮對點源函數求解進行進一步調整.

2.2 源函數調整準確性驗證

線源、點源簡化模型求解主要誤差來自于積分中值的選取,當裂縫微元沿積分方向的長度足夠小時,點源函數沿積分方向上單調變化,此時選取積分中點作為積分中值可取得較好的積分效果.圖7 為不同網格劃分下線源求解所得曲線變化對比.可以看出,當網格劃分較為粗糙時,線源法所得初期壓力及壓力導數曲線誤差較大.特別的壓力導數曲線存在較為明顯的回流特征,考慮是由于中心線所在線源與積分中值相差較大,各線源之間存在流向各微元中心線源的徑向流動所致(圖8).

圖7 網格劃分對線源求解效果影響Fig.7 Effect of meshing on line source solution

圖8 回流段流動示意圖Fig.8 Diagram of flow in flowback period

隨著劃分微元數目增多,簡化模型所得壓力及壓力導數曲線與面源積分所得壓力及壓力導數曲線趨于一致.大量模擬表明當所劃分微元無因次邊長在0.15 左右時,采用點源或線源代替面源的方法可獲得較高的精度.選取的單縫,劃分微元數目15 × 3,模擬獲得相應典型試井曲線如圖9所示.可以看出在tD大于10?7時間內采用點源或線源函數代替面源積分可取得較高計算精度,相同模型下,面源積分計算時長為435 443 s,而線源積分計算時長僅為17 s,點源積分計算時長僅為22 s,極大地提升了模型計算效率.

圖9 簡化模型求解效果示意圖Fig.9 Diagram of simplified model

3 典型曲線分析

3.1 流動階段劃分

通過分析典型試井曲線中的壓力及壓力導數曲線特征,可實現不同流動階段識別.選取無因次參數xeD=yeD= 20,hD= 0.5,θ= 45°,hfD=50,Cs= 10,CwD= 5,wD= 1.0 × 10?5,γ= 1.3 × 10?8.繪制考慮井儲的壓力及壓力導數曲線.從圖10 可以看出,單一傾斜縫典型試井曲線可劃分為8 個階段.

圖10 單一傾斜縫流動階段劃分Fig.10 Identification of the flow regimes of incline fracture

(1)井筒儲集及過渡階段:井儲階段壓力及壓力導數曲線重合,均為斜率等于1 的直線,該階段主要流體供給來自井筒;

(2)裂縫線性流階段:壓力導數曲線斜率為1/2,該階段基質暫未動用,與井相連的裂縫內部流體近線性流入井筒;

(3)雙線性流階段:壓力導數曲線斜率為1/4,該階段除裂縫內部流動外,近縫基質向裂縫的線性流動開始出現;

(4)基質線性流階段:壓力導數曲線斜率為1/2,該階段裂縫內壓力趨于穩定,近縫基質向裂縫的流動仍以線性流為主;

(5)早期徑向流階段:壓力導數曲線斜率為0,此時波及范圍擴大,基質向裂縫的供給由線性流轉變為徑向流;

(6)橢圓流階段:壓力導數曲線斜率近似為0.36,通常認為該階段為縫周基質徑向流向裂縫系統徑向流過渡的階段;

(7)晚期徑向流階段:壓力導數為斜率等于0 的直線,流體圍繞裂縫整體成徑向流;

(8)邊界控制流階段:流體流動到封閉邊界,壓力及壓力導數響應主要受邊界影響,壓力及壓力導數均為斜率等于1 的直線.

為研究裂縫條數對流動階段影響,取縫間距DD=1,對比裂縫條數分別為1,2,3 條件下壓力及壓力導數曲線(圖11)可以看出,隨著裂縫條數增加,主要流動階段仍為8 個,但與單條縫相比存在一定的差別.隨著裂縫段數增多,裂縫及近縫區域流動影響增強,井儲階段提前進入過渡階段,相應的雙線性流階段出現時間略有提前;由于縫間干擾的存在,裂縫條數越多,早期徑向流階段越不太明顯,過渡階段斜率更大,晚期徑向流及邊界控制流階段,流體相對于裂縫系統進行流動,不同裂縫條數下壓力及壓力導數曲線基本重合.

圖11 多段壓裂傾斜縫流動階段Fig.11 Flow regimes of multi-stage incline fractures

3.2 影響因素分析

本文所建模型在處理裂縫過程中主要進行了以下幾點考慮:多段壓裂縫、裂縫有限導流、裂縫以任意形態展布以及裂縫向任意方向傾斜,因此針對上述幾點考慮分別進行裂縫導流能力、裂縫傾角、井筒傾角、裂縫縱向高度以及裂縫段間距等參數敏感性分析.

(1)裂縫導流能力敏感性分析

圖12 給出了無因次裂縫導流系數分別為0.5,5,50 及500 時壓力及壓力導數曲線響應,可以看出隨著裂縫導流能力的升高,井儲階段及過渡階段更早結束,早期及晚期徑向流階段時間變短,特別是早期徑向流階段基本消失.

圖12 裂縫導流能力敏感性分析Fig.12 Sensitivity analysis of well storage coefficient

(2)裂縫傾角敏感性分析

考慮裂縫縱向高度保持一致,對比裂縫面傾角分別為15,30,45,60 以及75 時壓力及壓力導數曲線變化,由圖13 可以看出,裂縫面傾角主要影響雙線性流階段至橢圓流階段的壓力導數曲線,當裂縫傾角小于45°時,傾角影響較小,隨著傾角進一步增大,雙線性流作用時間增長,早期徑向流更為明顯,相應的橢圓流階段傾角增大.

圖13 裂縫傾角敏感性分析Fig.13 Sensitivity analysis of well storage coefficient

(3)裂縫高度敏感性分析

取裂縫傾角為45°,考慮裂縫無因次縱向高度在0.2–1.0 范圍內變化.由圖14 可以看出,當裂縫高度較小時,裂縫內線性流即雙線性流階段較不明顯,且基質線性流后期存在一定的回流現象,隨著裂縫高度的升高,隨著裂縫高度的增大,各階段流態逐漸出現.

圖14 裂縫高度敏感性分析Fig.14 Sensitivity analysis of fracture height

(4)裂縫段間距敏感性分析

取裂縫傾角為15°,考慮無因次裂縫段間距在1–4 范圍內變化,由圖15 可以看出,隨著裂縫段間距增大,早期徑向流至晚期徑向流階段的特征更為不明顯,表明隨著裂縫段間距的增大,縫間干擾越晚出現且系統徑向流出現時間越晚.

圖15 段間距敏感性分析Fig.15 Sensitivity analysis of fracturing interval

4 結論

(1)對于任意方向展布裂縫可通過裂縫離散表征,耦合裂縫內數值解與基質內面源函數解析解的方法進行求解,但存在計算效率低的問題.借助點源、特殊線源替代面源的方法可大大提升計算效率,且在裂縫微元劃分較為精細(微元無因次邊長小于0.15)時精度較高.

(2)傾斜縫典型試井曲線可劃分為井筒儲集及過渡流、裂縫線性流、雙線性流、基質線性流、早期徑向流、橢圓流、晚期徑向流及邊界控制流八個流動階段.隨著裂縫段數增多,縫間影響增強,井儲、過渡階段及徑向流階段特征減弱.

(3)裂縫傾角及裂縫段間距主要影響裂縫線性流到早期徑向流階段壓力及壓力導數曲線,特別的當裂縫高度較小時壓力導數存在回流現象.