基于XFEM-MBEM 的嵌入式離散裂縫模型流固耦合數(shù)值模擬方法1)

杜旭林 程林松 牛烺昱 方思冬 曹仁義 ,3)

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

? (中國石化石油勘探開發(fā)研究院,北京 100083)

引言

非常規(guī)裂縫性儲層中,巖石基質(zhì)提供油氣儲存空間,人工裂縫、誘導(dǎo)裂縫和天然裂縫提供其流動通道[1].目前油田實際歷史數(shù)據(jù)已表明,地質(zhì)力學(xué)特性對裂縫性油藏的產(chǎn)量起重要作用,會導(dǎo)致一些關(guān)鍵滲流參數(shù)發(fā)生改變,如基巖的孔隙度和滲透率;此外,致密油藏以“縫控儲量”為主[2],巖石形變會引起裂縫開度變化,這對于壓裂水平井產(chǎn)能評價的影響極大.傳統(tǒng)的應(yīng)力敏感模型無法揭示深層次的科學(xué)道理[3-4],為了準(zhǔn)確模擬此類裂縫性多孔介質(zhì)中的裂縫動態(tài)行為和流體流動機理,需要綜合考慮流體和巖石兩方面的力學(xué)性質(zhì).

巖石力學(xué)中有效應(yīng)力的概念最早由Terzaghi[5]和Biot[6]提出,并由Biot[7]發(fā)展為較完善的三維固結(jié)理論,這是研究應(yīng)力場與滲流場耦合的早期基礎(chǔ).為了將上述的理論基礎(chǔ)應(yīng)用于傳統(tǒng)的滲流模型,Barenblatt和Zheltov[8]首次提出了雙重孔隙度/雙重滲透率的概念,基巖與裂縫為兩種重疊的、不同的連續(xù)介質(zhì),其中基巖完全被裂縫所圍繞,用于表征裂縫性多孔介質(zhì)中的流體和巖石形變.后續(xù),Bai[9]提出了關(guān)于雙重介質(zhì)更為嚴(yán)格的物理解釋即雙孔有效應(yīng)力理論,用于分析不同尺度的裂縫和基質(zhì)的形變.針對多孔介質(zhì)流固耦合數(shù)學(xué)模型的建立,Wilson 和Aifantis[10]、Beskos 和Aifantis[11]和Khaled 等[12]建立了基于有限元法的雙重介質(zhì)流固耦合模型,該模型的主要問題在于其忽視了基質(zhì)與裂縫之間壓差導(dǎo)致的耦合流動作用;Lewis 和Schrefler[13]對于油藏中的流固耦合問題進(jìn)行了研究,分別考慮了單相流、多相流、線彈性變形以及彈塑性變形,并給出了相應(yīng)的有限元計算格式.在上述早期的研究中,采用連續(xù)介質(zhì)力學(xué)理論表征裂縫性油藏,將不連續(xù)的裂縫系統(tǒng)粗化為連續(xù)介質(zhì),能在一定程度上刻畫基巖和裂縫的流體流動及力學(xué)特征,但對于裂縫的幾何形狀及傳導(dǎo)率具有較強的約束,更適用于描述小尺度裂縫,對于大尺度裂縫處理精度較低,無法滿足以體積壓裂為開發(fā)手段的非常規(guī)油氣藏的礦場實際需求.

近年來,基于離散裂縫的流固耦合理論受到了廣泛關(guān)注,此類模型將裂縫離散化顯式表征,能準(zhǔn)確刻畫大尺度裂縫對滲流場和應(yīng)力場的影響,根據(jù)前處理算法的不同,可劃分為離散裂縫模型[14](DFM)與嵌入式離散裂縫模型[15-16](EDFM),其中嵌入式裂縫模型最早由Li 和Lee[17]開發(fā),旨在提高天然裂縫建模的精確程度.與使用復(fù)雜非結(jié)構(gòu)網(wǎng)格在空間上匹配裂縫幾何結(jié)構(gòu)的DFM 方法相比,EDFM 只需使用一組固定的結(jié)構(gòu)化矩形網(wǎng)格,裂縫在前處理中與網(wǎng)格分開并顯式精確刻畫,這是EDFM 相對于DFM的關(guān)鍵優(yōu)勢[18].流固耦合研究方面,由于裂縫開度與油藏尺度存在極大差異,需對裂縫進(jìn)行降維處理避免裂縫周圍局部加密以減小網(wǎng)格量,而降維后裂縫兩側(cè)的位移場存在間斷性[19].由于EDFM 采用非匹配性網(wǎng)格,使得其在表征裂縫力學(xué)間斷效應(yīng)方面極具挑戰(zhàn)性.目前主要處理方法有: 有限元法、位移不連續(xù)法和擴展有限元法[20](XFEM).其中,XFEM 與EDFM 具有相似的優(yōu)點,兩者均基于結(jié)構(gòu)化網(wǎng)格,分別相較于傳統(tǒng)的基于有限元法的應(yīng)力場、滲流場計算,避免了網(wǎng)格重構(gòu)和非結(jié)構(gòu)化網(wǎng)格剖分帶來的困難.基于此,Ren 等[21-22]、Yan 等[23]、Zeng 等[24]先后實現(xiàn)了將EDFM 與用于應(yīng)力場計算的擴展有限元法的耦合研究;Ren 等[25]進(jìn)一步開發(fā)了pEDFM-XFEM流固耦合數(shù)值模擬器,能更精確地模擬多相流情況下的巖石和流體的力學(xué)性質(zhì).

在上述前人的研究中,仍亟待解決的主要問題有兩點: 致密基質(zhì)的滲透率極低,大規(guī)模體積壓裂之后儲層中的大尺度水力壓裂縫與大量被激活的、小尺度的天然裂縫共存,對此類“人工裂縫性油藏”,單一數(shù)學(xué)模型無法進(jìn)行準(zhǔn)確表征,需要更一步發(fā)展綜合數(shù)學(xué)模型;目前現(xiàn)有數(shù)值模擬方法在計算包含裂縫單元的基質(zhì)網(wǎng)格內(nèi)的壓力分布時,采用了線性分布假設(shè),這導(dǎo)致了早期及非穩(wěn)態(tài)流動計算精度的不足.本文以上述突出的科學(xué)問題為導(dǎo)向,建立了致密油藏基質(zhì)-天然裂縫-水力壓裂縫綜合數(shù)學(xué)模型,基于單孔模型表征基質(zhì)滲流參數(shù)隨著壓力、應(yīng)變場的改變,基于雙重介質(zhì)模型表征天然裂縫滲流參數(shù)的變化,基于嵌入式離散裂縫模型表征水力壓裂縫的動態(tài)特征;并提出了XFEM-MBEM 的混合數(shù)值離散化方法,用于計算基巖和裂縫間的非穩(wěn)態(tài)竄流,提高了模擬的早期精度,可準(zhǔn)確表征致密油藏開采過程中的裂縫變形及流體流動機理,此研究為致密油氣開發(fā)領(lǐng)域提供了理論基礎(chǔ).

1 流固耦合數(shù)學(xué)模型

1.1 物理模型及基本假設(shè)

本文采用了黑油模型假設(shè),考慮基質(zhì)和流體的壓縮性,油水兩相流動均符合達(dá)西定律,等溫滲流,巖石變形考慮為微小線彈性變形,力學(xué)方程中采用常規(guī)彈性力學(xué)符號約定,即拉伸為正,壓縮為負(fù).

1.2 準(zhǔn)靜態(tài)力學(xué)平衡方程

圖1 為彈性靜力平衡狀態(tài)下的多孔介質(zhì),中間存在一條水力壓裂縫,pF為作用裂縫內(nèi)表面的縫內(nèi)孔隙壓力,pp為支撐劑顆粒的支撐力,表示裂縫表面兩側(cè) ΓF上的兩個單位法向量,nt表示基巖應(yīng)力外邊界 Γt上的單位外法向量,t為應(yīng)力外邊界上的牽引力,為位移外邊界 Γu上的平均固體位移.由此可建立準(zhǔn)靜態(tài)應(yīng)力平衡方程式(1)～ 式(4)

圖1 彈性靜力平衡狀態(tài)下的多孔介質(zhì)Fig.1 Porous medium in elastostatic equilibrium state

式中,σ 為總應(yīng)力張量,Pa;u為巖石固體位移,m;Pm為基巖孔隙壓力,Pa;ρb為多孔介質(zhì)基巖密度,kg/m3;g為重力加速度,N/kg.

1.3 油水滲流控制方程

基質(zhì)體積單元中的油水兩相質(zhì)量平衡方程

離散裂縫單元中的油水兩相質(zhì)量平衡方程

式中,k為滲透率,m2;kro為油相相對滲透率,小數(shù);krw為水相相對滲透率,小數(shù); μ 為流體黏度,mPa·s;s為流體飽和度,小數(shù);B為流體體積系數(shù),小數(shù);p為流體壓力,Pa;要補充的是,上述式中下標(biāo)o 和w 分別表示油相和水相,上標(biāo)m和F分別表示基巖和水力壓裂縫,上標(biāo)w表示井筒表示基巖與裂縫之間的竄流量表示基巖與井筒之間的流量;φ*表示為Lagrangen 孔隙度,后續(xù)進(jìn)行詳細(xì)闡述.

1.4 基質(zhì)-裂縫傳質(zhì)方程

EDFM 中存在4 類非相鄰連接關(guān)系 (NNRs): 相鄰基質(zhì)網(wǎng)格之間的連接;位于不同基質(zhì)網(wǎng)格中的相鄰裂縫單元之間的連接;位于同一基質(zhì)網(wǎng)格中的相鄰裂縫單元之間的連接;基質(zhì)網(wǎng)格與所包含裂縫單元之間的連接.其中,傳導(dǎo)率可表示為流度與其幾何因子的乘積,以油相為例,流度和幾何因子可表示為

式中,λij表示單元i和j之間的油相流度;kro(upstream)表示利用上游權(quán)格式取值的油相相對滲透率;Gij為單元i和j之間的幾何因子;TIij表示單元i和j之間的傳導(dǎo)系數(shù).

基質(zhì)與裂縫之間傳質(zhì)項的處理是基于EDFM 中的第4 類連接關(guān)系,原始EDFM 方法[26]在計算包含裂縫單元的基質(zhì)網(wǎng)格內(nèi)的壓力分布時,采用了Lee 等[27]及Praditia 等[28]提出的壓力與到裂縫面的垂向距離成正比的線性分布假設(shè)如下

式中ki表示單元i的滲透率,m2;A代表接觸面積,m2;〈d〉代表基質(zhì)網(wǎng)格中心點與裂縫面之間的距離,m;S為所在基質(zhì)網(wǎng)格的面積,m2.

但在實際物理情景中由于裂縫與基質(zhì)的滲透率級別相差極大,這種求取“平均法向距離”的近似方法在處理非穩(wěn)態(tài)傳質(zhì)時會導(dǎo)致一定的誤差.因此,本文基于Fang 等[29]與Cao 等[30]所提出的混合邊界元法,利用穩(wěn)態(tài)滲流控制方程的邊界積分方程推導(dǎo)得到雙孔基質(zhì)網(wǎng)格向裂縫網(wǎng)格傳質(zhì)量的新近似格式,以油相為例見式(16),EDFM 中MBEM 流量處理的連接關(guān)系示例見圖2 所示.

圖2 雙孔系統(tǒng)EDFM 中MBEM 流量項處理的連接關(guān)系示例Fig.2 Example of connection list for MBEM handling flux in EDFM with dual-porosity system

式中,qomF為從裂縫單元流入基質(zhì)網(wǎng)格的油相流量,m3;下標(biāo)nF代表裂縫單元個數(shù);h為基質(zhì)網(wǎng)格厚度,等矩陣形式見附錄.

1.5 井方程

經(jīng)典Peaceman 公式[31]可用于計算井筒與基質(zhì)或裂縫之間的流量.以油相為例,對于壓裂水平井井筒與裂縫單元相連接

式中,wF為水力壓裂縫的開度,m;re為等效補給井徑,m;rw為油井井徑,m;l1和l2分別為裂縫矩形單元的兩個邊長,m.

2 多孔介質(zhì)本構(gòu)關(guān)系

2.1 修正的Biot 方程

為準(zhǔn)確描述非常規(guī)儲層中基質(zhì)與小尺度裂縫間的流動和巖石形變,可采用由Bai[32]提出的雙孔有效應(yīng)力原理,對基質(zhì)和裂縫分別建立具有嚴(yán)格物理意義的有效應(yīng)力關(guān)系

式中,σm表示基質(zhì)總應(yīng)力張量,Pa; σf為裂縫總應(yīng)力張量為基質(zhì)有效應(yīng)力表示裂縫有效應(yīng)力,Pa; αm表示畢渥系數(shù),小數(shù); ε 為總應(yīng)變張量,m; εm為基質(zhì)應(yīng)變張量,m; εf為裂縫應(yīng)變張量,m;上述式中下標(biāo)f 表示小尺度裂縫.

Ren 根據(jù)式(19)～ 式(22)及胡克定律,指出雙重介質(zhì)有效應(yīng)力的表達(dá)式為[22]

式中,Dmf為彈性剛度張量,Pa;Cm為基質(zhì)柔度張量,Pa;Cf為裂縫柔度張量,Pa.

2.2 孔滲參數(shù)表征方法

孔隙度是關(guān)于體積應(yīng)變和流體壓力的函數(shù)[25].考慮線彈性微小形變下的基巖Lagrangen 孔隙度和大尺度裂縫Lagrangen 孔隙度可分別由下式表征

本文中將小尺度天然裂縫和基巖視為雙重連續(xù)介質(zhì),因此小尺度裂縫表征方法與基質(zhì)相同.對于大尺度水力壓裂縫,生產(chǎn)過程中由于縫內(nèi)流體壓力損失會導(dǎo)致裂縫面的閉合,然而當(dāng)存在支撐劑時,支撐劑發(fā)生形變會對裂縫面施加支撐力阻止其閉合.根據(jù)胡克定律,由支撐劑引起的支撐力可表示為

式中,Ep為支撐劑彈性模量,Pa;w為裂縫開度的位移,為水力壓裂縫的原始開度,m.

此過程中,水力壓裂縫的實際導(dǎo)流能力依賴于裂縫滲透率和開度的更新迭代[33]

3 耦合模型的數(shù)值離散

3.1 有限體積法離散流動方程

XFEM 與EDFM 具有相似的優(yōu)點,兩者均基于結(jié)構(gòu)化網(wǎng)格,網(wǎng)格劃分容易且物理意義清晰.因此本文采用正交矩形網(wǎng)格對基質(zhì)進(jìn)行幾何離散,針對EDFM 前三類非相鄰連接關(guān)系,采用滿足局部物質(zhì)守恒且簡單易行的塊中心有限體積方法來獲取滲流控制方程的離散格式,仍以油相為例,此時為公式簡潔,方程中去除可能存在的點源/匯項,對式(5)的兩側(cè)對時間和控制體積進(jìn)行積分

式中左側(cè)利用散度定理將體積分轉(zhuǎn)換為面積分,并用符合物理意義的相鄰網(wǎng)格之間的法向流量之和近似;右側(cè)對時間的積分精確求解,對體積分采用矩形法估計,由此式(29)可以寫為

對于相鄰基質(zhì)網(wǎng)格如圖3,式(30)取隱式格式可改寫為

圖3 矩形網(wǎng)格離散示意圖Fig.3 Sketch of discretization for rectangle grids

式中,n表示與該基質(zhì)網(wǎng)格相鄰的網(wǎng)格數(shù)量,ΔVi為該基質(zhì)網(wǎng)格的體積,Δt為相鄰時間步.

根據(jù)式(29)～ 式(31)的推導(dǎo)過程可以獲得油水相的全隱式離散方程,但對于裂縫與基質(zhì)連接的第4 類非相鄰連接關(guān)系即含傳質(zhì)項,需要進(jìn)行單獨考慮.這里將已推導(dǎo)的多相流情況下EDFM 第4 類連接中基質(zhì)網(wǎng)格與裂縫網(wǎng)格之間傳質(zhì)量的近似格式(16),耦合到基質(zhì)與裂縫之間多相流方程的有限體積離散格式中,這也是混合邊界元EDFM 的核心思想.仍以油相為例,對于基質(zhì)系統(tǒng)

同理,對于裂縫系統(tǒng)

3.2 擴展有限元離散力學(xué)方程

XFEM 是在常規(guī)的有限元位移模式中根據(jù)單位分解的思想加引進(jìn)一個跳躍函數(shù)和裂尖漸進(jìn)位移場以反映位移不連續(xù)性的數(shù)值方法.XFEM 中的不同節(jié)點和單元的類型見圖4,單元位移的近似解可以寫成

圖4 混合模型離散化方法示意圖Fig.4 Sketch of hybrid model discretization strategy

式中,I表示網(wǎng)格中所有節(jié)點的集合,L表示包含所有裂縫段所在單元的節(jié)點集合,M為裂縫與裂縫相交的節(jié)點集合,K表示裂縫尖端所在單元的節(jié)點集合;Ni為節(jié)點位移形函數(shù),ui為節(jié)點位移矢量,ai,j,ci,j,為單元增強節(jié)點的額外自由度,H表示Heaviside 函數(shù),J表示裂間連接函數(shù),F表示裂尖漸進(jìn)函數(shù),三者的表達(dá)式如下所示

將式(19)和式(20)代入式(1),考慮邊界條件式(2)～ 式(4),得到方程(1)的等效積分弱形式見式(38),弱形式方程包含不同尺度裂縫以及基巖流體壓力的相互作用

式中,δε 表示試應(yīng)變,δu表示試位移,δu=δu+-δu-.

裂縫開度位移可以表示為

采用標(biāo)準(zhǔn)伽遼金有限元法對系統(tǒng)整體進(jìn)行離散,可獲得一個線性平衡方程組如下

式中,K表示整體剛度矩陣,f表示整體力向量.

3.3 整體計算格式及求解策略

本文分別采用有限體積法和擴展有限元法求解滲流控制方程及力學(xué)平衡方程,主要變量包括基巖及裂縫中的流體壓力、含水飽和度以及分別設(shè)置在網(wǎng)格中心和角點的節(jié)點位移,對時間項的離散采用一階向后歐拉全隱式差分格式,對完全耦合的全局方程組采用Newton-Raphson 迭代求解,并利用自動微分算法計算迭代過程中的雅可比矩陣

式中系數(shù)矩陣是非對稱的,其中反對角項是主變量之間的耦合項.

4 模型驗證分析

4.1 連續(xù)性介質(zhì)流固耦合算例

此算例選取了經(jīng)典的一維彈性土固結(jié)的孔彈性問題,將本文模型與Terzaghi[5]所提出的解析解進(jìn)行對比.假設(shè)等溫多孔介質(zhì)由單相流體和固體組成,表現(xiàn)為線性孔彈性特征,模型的寬度為5 m,深度為10 m,頂部為排水邊界,其余邊界均不排水,底部為零位移邊界,左右邊界位移在水平方向上約束,于頂部施加20 MPa 的恒定載荷.巖土的孔隙度為0.4,滲透率為100 mD,彈性模量為20 GPa,泊松比為0.2,孔隙比為0.4,一維模型數(shù)值解與解析解對比結(jié)果如圖5 所示,從圖中可看出不同時間下本文提出的流固耦合數(shù)值解與Terzaghi 解析解的結(jié)果較為吻合,誤差在合理范圍之內(nèi),能證明該模型的準(zhǔn)確性.

圖5 一維模型數(shù)值解與解析解的對比Fig.5 Comparison between numerical solution and analytical solution of one dimensional model

4.2 離散裂縫單相流算例

此算例是通過多物理場耦合模擬的商業(yè)模擬器COMSOL 的標(biāo)準(zhǔn)有限元 (SFEM) 驗證本文建立的流動模塊.物理模型為矩形油藏,油藏中心有一條水力壓裂縫,設(shè)置裂縫半長為基礎(chǔ)長度將整個區(qū)域無因次化,裂縫無因次長度為4,油藏在X方向的長度為10 且Y方向長度為10,見圖6.基于拉氏空間下無因次滲流方程,將本文修正后的EDFM 和有限元法與Blasingame 精確解[34]進(jìn)行對比,其中裂縫無因次導(dǎo)流能力為200,有限元網(wǎng)格裂縫采用局部加密如圖7 所示,裂縫加密的個數(shù)為20,如圖8 所示.本文模型采用了非匹配性的矩形網(wǎng)格和混合邊界元法,不對裂縫進(jìn)行局部加密,分析離散裂縫壓力動態(tài)曲線中不同數(shù)值算法的適應(yīng)性.

圖6 單裂縫矩形油藏模型Fig.6 Single fracture in rectangular reservoir model

圖7 對裂縫網(wǎng)格加密示意圖Fig.7 Sketch of fracture mesh refinement

圖8 不同數(shù)值解與解析解的對比Fig.8 Comparison of different numerical solutions and analytical solution

模擬結(jié)果如圖7 所示,模擬無因次時間從10-9～10,覆蓋了有限導(dǎo)流能力的幾個階段,可以看出修正后的EDFM 方法與精確解無論從早期還是晚期都較為吻合,而有限元模型在早期與解析解的差距很大,中后期曲線重合,這在圖9 的誤差分析中也可以看出.原因在于有限元法處理網(wǎng)格累積項時,采用線性插值積分,在一個時間步內(nèi)是穩(wěn)態(tài)過程,而本文模型采用混合邊界元法處理時采用了非穩(wěn)態(tài)滲流的基本解,能精確描述初期的流動特征,根據(jù)對比結(jié)果可以說明本文模型的流動模塊具有較高的早期精度.

圖9 誤差分析Fig.9 Error analysis

4.3 裂縫性介質(zhì)應(yīng)力強度因子算例

此算例是通過經(jīng)典的Tada 解析解[35]驗證本文建立的力學(xué)模塊計算裂縫尖端應(yīng)力強度因子 (SIFs)的準(zhǔn)確性.考慮一個平面應(yīng)變模型 (10 m×10 m),上邊界為施加均勻的拉應(yīng)力,下邊界為固定位移的滾輪邊界,模型內(nèi)部中心存在一條傾斜裂縫,裂縫半長l為1 m,如圖10 所示.解析解中應(yīng)力強度因子是外部應(yīng)力和裂縫幾何尺寸和傾角的函數(shù),因此為方便對比,定義無因次應(yīng)力強度因子:在裂縫尺寸保持不變的條件下,計算不同裂縫傾角無因次應(yīng)力強度因子數(shù)值解與解析解的結(jié)果見圖11,對比顯示本文模型的數(shù)值解與經(jīng)典解析解的結(jié)果幾乎保持一致,具有較高的可信度.

圖10 裂縫性介質(zhì)示意圖Fig.10 Schematic of fractured medium

圖11 不同裂縫傾角無因次應(yīng)力強度因子數(shù)值解與解析解的對比Fig.11 Comparison of numerical and analytical solutions for dimensionless SIFs with different fracture dip angles

5 典型實例分析

如圖12 所示的致密油藏,考慮分為壓裂改造區(qū)(SRV) 與非壓裂改造區(qū) (USRV),SRV 內(nèi)部存在一口水平井和7 條水力壓裂縫,水力壓裂縫中考慮支撐劑的作用.采用平面應(yīng)變模型,油藏邊界條件如圖12所示.SRV 內(nèi)考慮微裂縫的影響,采用雙孔模型;USRV 不考慮微裂縫的影響,采用單孔模型.致密油藏的基本參數(shù)見表1,設(shè)置該水平井為定壓生產(chǎn)1000 d,井底流壓為10 MPa,并將模擬的結(jié)果與商業(yè)模擬器COMSOL 的標(biāo)準(zhǔn)有限元 (SFEM) 作對比,圖13 為兩種數(shù)值模擬器的網(wǎng)格剖分,其中COMSOL-SFEM 中對SRV 內(nèi)部及裂縫尖端進(jìn)行了局部加密處理,本文模型XFEM-MBEM 中采用非匹配性的矩形網(wǎng)格.圖14 為兩種數(shù)值模擬器在開發(fā)1000 d 時刻的壓力場對比,圖15 和圖16 分別為X-位移場與Y-位移場的對比,從上述場圖可看出兩種數(shù)值模擬器的結(jié)果基本一致.以圖12 中SRV 內(nèi)部最中間的那條水力壓裂縫為研究目標(biāo),分析其不同時刻的開度分布見圖17,圖中能看出XFEM-MBEM 與COMSOL-SFEM不同時刻裂縫開度的計算結(jié)果較為吻合,可以說明XFEM-MBEM 流固耦合模擬的準(zhǔn)確性.此外,從圖中可以隨著油藏的不斷開發(fā),裂縫呈不斷閉合的趨勢,定壓生產(chǎn)下靠近井筒附近的裂縫段閉合更加劇烈.圖18 為該模型在表1 參數(shù)條件下考慮流固耦合作用與不考慮流固耦合作用的生產(chǎn)動態(tài)曲線對比,可以看出應(yīng)力場所引起的滲流參數(shù)改變及裂縫開度降低對產(chǎn)能的影響很大,因此對致密油藏壓裂水平井進(jìn)行產(chǎn)能評價時,地質(zhì)力學(xué)效應(yīng)的影響不可忽略.

圖12 致密油藏示意圖Fig.12 Sketch of tight oil reservoir

表1 致密油藏基本參數(shù)Table 1 Basic parameters of tight oil reservoir

圖15 兩種數(shù)值模擬器開發(fā)1000 d 的X-位移場對比(單位: mm)Fig.15 Comparison of X-displacement distribution maps over 1000 days of two simulators (unit: mm)

圖16 兩種數(shù)值模擬器開發(fā)1000 d 的Y-位移場對比(單位: mm)Fig.16 Comparison of Y-displacement distribution maps over 1000 days of two simulators (unit: mm)

圖17 不同時刻裂縫開度分布Fig.17 Fracture aperture distribution at different time

圖18 考慮與不考慮流固耦合作用的生產(chǎn)動態(tài)曲線對比Fig.18 Comparison of production dynamic curves with and without flow-geomechanics coupled effect

6 結(jié)論

(1)本文提出了基于XFEM-MBEM 的嵌入式離散裂縫模型流固耦合數(shù)值模擬方法,采用正交的非匹配型網(wǎng)格進(jìn)行幾何離散,其優(yōu)勢在于裂縫獨立于計算網(wǎng)格,與常規(guī)離散裂縫模型基于裂縫形狀的匹配型網(wǎng)格相比,網(wǎng)格劃分難度大大降低,能避免當(dāng)裂縫間距或夾角非常小時,網(wǎng)格劃分質(zhì)量不佳等問題.本文通過幾個實例分別驗證了其力學(xué)模塊、流動模塊、耦合模塊的準(zhǔn)確性.

(2)傳統(tǒng)EDFM 方法對于第4 類非相鄰連接關(guān)系在計算包含裂縫單元的基質(zhì)網(wǎng)格內(nèi)的壓力分布時采用了線性分布假設(shè),導(dǎo)致了開發(fā)早期及非穩(wěn)態(tài)流動計算精度的不足.局部網(wǎng)格加密技術(shù)可以減小包含裂縫網(wǎng)格的基質(zhì)網(wǎng)格尺寸,從而有效解決這個問題,但會對前處理算法帶來巨大的復(fù)雜性.本文采用MBEM 方法獲取了雙孔基質(zhì)網(wǎng)格向裂縫網(wǎng)格傳質(zhì)量的新近似格式,能在不加密情況下提高早期計算精度.

(3) SRV 改造區(qū)內(nèi)的微裂縫對提高儲層流體的動用程度有重要影響.然而,用EDFM 和XFEM 進(jìn)行顯式表征是較為困難的,且會增加計算耗時.本文采用雙孔有效應(yīng)力原理耦合雙重介質(zhì)模型隱式表征SRV 區(qū)內(nèi)部的天然裂縫,能有效模擬具有復(fù)合分區(qū)特征的此類裂縫性油藏問題.

(4)致密油藏開發(fā)過程中隨著縫內(nèi)流體被采出,裂縫會出現(xiàn)不同程度的閉合,從而導(dǎo)致了整體導(dǎo)流能力降低,油藏采收率將降低.考慮到非常規(guī)儲層原油的產(chǎn)量對裂縫的依賴性如此之高.因此,考慮開發(fā)過程中流固耦合作用引起的參數(shù)變化對產(chǎn)量的評價具有重要意義.

附錄