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

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

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

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

引言

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

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

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

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

1 流固耦合數學模型

1.1 物理模型及基本假設

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

1.2 準靜態力學平衡方程

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

圖1 彈性靜力平衡狀態下的多孔介質Fig.1 Porous medium in elastostatic equilibrium state

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

1.3 油水滲流控制方程

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

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

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

1.4 基質-裂縫傳質方程

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

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

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

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

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

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

式中,qomF為從裂縫單元流入基質網格的油相流量,m3;下標nF代表裂縫單元個數;h為基質網格厚度,等矩陣形式見附錄.

1.5 井方程

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

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

2 多孔介質本構關系

2.1 修正的Biot 方程

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

式中,σm表示基質總應力張量,Pa; σf為裂縫總應力張量為基質有效應力表示裂縫有效應力,Pa; αm表示畢渥系數,小數; ε 為總應變張量,m; εm為基質應變張量,m; εf為裂縫應變張量,m;上述式中下標f 表示小尺度裂縫.

Ren 根據式(19)～ 式(22)及胡克定律,指出雙重介質有效應力的表達式為[22]

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

2.2 孔滲參數表征方法

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

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

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

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

3 耦合模型的數值離散

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

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

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

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

圖3 矩形網格離散示意圖Fig.3 Sketch of discretization for rectangle grids

式中,n表示與該基質網格相鄰的網格數量,ΔVi為該基質網格的體積,Δt為相鄰時間步.

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

同理,對于裂縫系統

3.2 擴展有限元離散力學方程

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

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

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

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

式中,δε 表示試應變,δu表示試位移,δu=δu+-δu-.

裂縫開度位移可以表示為

采用標準伽遼金有限元法對系統整體進行離散,可獲得一個線性平衡方程組如下

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

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

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

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

4 模型驗證分析

4.1 連續性介質流固耦合算例

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

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

4.2 離散裂縫單相流算例

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

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

圖7 對裂縫網格加密示意圖Fig.7 Sketch of fracture mesh refinement

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

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

圖9 誤差分析Fig.9 Error analysis

4.3 裂縫性介質應力強度因子算例

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

圖10 裂縫性介質示意圖Fig.10 Schematic of fractured medium

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

5 典型實例分析

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

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

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

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

圖16 兩種數值模擬器開發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 考慮與不考慮流固耦合作用的生產動態曲線對比Fig.18 Comparison of production dynamic curves with and without flow-geomechanics coupled effect

6 結論

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

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

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

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

附錄