兩套節點格林元嵌入式離散裂縫模型數值模擬方法1)

程林松 杜旭林 饒 翔 曹仁義 賈 品

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

? (長江大學石油工程學院,武漢 430100)

引言

對非常規裂縫性儲層的壓力和流體分布進行精確模擬,是許多能源工程技術人員面臨的一個熱點和挑戰性工作[1].目前為止,已有許多數值方法被用于解決這一物理問題,如等效連續介質模型[2]、離散裂縫模型[3]和嵌入式離散裂縫模型(embedded discrete fracture method,EDFM)[4-5].由于具有較高的關鍵地質特征可視化水平,其中EDFM 的應用最為廣泛.但由文獻[6-8]所提出的原始EDFM 在計算基質網格與裂縫單元間傳導系數的幾何因子時,采用了壓力與到裂縫面的垂向距離成正比的線性假設,由于非常規儲層裂縫與基質的滲透率極差較大,這種理想假設無法準確描述早期的非穩態竄流,會導致一定的誤差;此外,原始EDFM 多采用僅適用于矩形網格剖分的有限差分或有限體積法,難以考慮復雜油藏邊界,技術應用層面存在局限性.

針對EDFM 早期精度不足、網格適用性不佳的問題,目前國內外有兩個主要研究思路: 其一是優化EDFM 網格剖分的前處理算法[9-11];其二是尋求一種更加高效、適用性更強的數值解法.針對基質網格和裂縫單元之間傳質量的準確計算問題,由于邊界元法相較于有限元法具有半解析精度,可能會是一種較為有效的思路[12].然而,由于實際儲層的非均質性極強,原始邊界元法無法有效處理這一問題.格林元法(Green element method,GEM)本質上是邊界元法的一種變體,在剖分的每個網格中建立邊界積分方程,結合了有限元法的變分原理,更適用于解決非均勻介質的非線性問題.原始GEM 由Taigbenus等[13-17]提出,核心思想是將計算域由多邊形單元劃分,單元頂點被視為未知節點,通過用網格內部壓力近似表達式的法向導數估計邊界上的法向流量,但僅顯式考慮了節點的壓力值,因此總體精度不高,誤差會隨著多邊形單元尺寸的減小而增大.后續的學者對原始GEM 進行了改進,Archer 等[18-19]指出,采用較大網格范圍內構造的插值基函數能減少原始GEM 的誤差;Pecher等[20]和Lorinczi 等[21-24]提出了基于流量向量的格林元法;Taigbenus[25-26]后續發展了修正流量項的GEM,但始終未解決無法顯式求解的根本問題;方思冬等[27]結合了混合有限元法和邊界元法二者的優勢,將壓力節點和流量節點設置在網格邊的中心點,建立了物理意義明確的混合邊界元法;為提高GEM 的魯棒性,Rao 等[28]提出了模擬格林元法,該方法耦合了模擬有限差分法的核心思想,可滿足局部質量守恒,能穩定地求解二階偏微分方程;Rao 等[29]還提出了兩套節點格林元的基本思想,但沒有實現與EDFM 的耦合及相關應用.綜合上述研究來看,目前大部分格林元法在精度、收斂性及適應的網格類型上各有優劣性,但若要將格林元的優勢引入EDFM 的發展中,就必須建立一種同時保持高精度、強魯棒性且能適用于非結構網格的格林元法.

本文從原始EDFM 的適用局限性出發,利用格林元方法在處理非定常和非均質問題方面具有高精度的獨特優勢,建立一種能與EDFM 合理匹配的改進格林元法,能精確求解基質與裂縫間的非穩態竄流問題,拓展原始EDFM 技術應用層面的局限性,使其既能求解各類復雜油藏滲流問題,又能使求解結果具有明確的物理意義.

1 格林元法雛形及其局限性

本文研究的是二維孔隙空間中的單相等溫滲流問題,流動遵守Darcy 定律,并考慮封閉外邊界條件.原始GEM[13]的基本思想是將邊界積分方程應用至各個子單元中,離散格式與直接法推導的有限元法較為相似,被廣泛應用在如式(1)的二階偏微分方程中

式中,p是計算域內基質的孔隙流體壓力,MPa;K和c為滲流介質的屬性參數函數,均為標量;f是計算域內的源匯項強度,若滲流介質某基質網格內含有裂縫單元,則有,其中qomf表示該基質網格與裂縫單元間的傳質速度,m3/d;V是基質網格體積,m3.

上式可以改寫為

式中,ψ=lnK,ν=1/K,σ=cv.

拉氏方程在無限大計算域的基本解G(M,Mi) 為

式中,Mi為空間中的源點坐標,M為空間中任意一點的坐標,δ(M,Mi) 為Delta 函數.

根據格林第二公式和式(3),可以得到邊界積分方程如下

式中,λ 為邊界積分方程中所選取源點所處位置的特征角度,Ω和 ?Ω 分別表示計算域及計算域的邊界,n是計算域邊界上的外法向量,i是所選取源點的序號.

Taigbenu 等[13-17]利用了法向流量的概念,即令q=-K?p·n,則上式改寫為

從式(5)可以看出,式中存在與流體壓力和介質屬性參數有關的域積分,由此體現了GEM 與邊界元法的顯著差異性,其可以有效處理計算域非均質及參數非線性等更為復雜的問題.以前學者提出過的原始GEM[13-17]、基于流量向量的GEM[21-24]及流量修正的GEM[25],都是在上述GEM 雛形基礎上的進一步發展,此類經典方法存在的突出共性問題在于子單元頂點處的外法向流量 (qn=?p·n) 不具備整體連續性,使得局部物質不守恒.

2 壓力流量兩套節點格林元法

為此,本文提出了一種壓力流量兩套節點的改進格林元法,用一組節點包含表示壓力值的多邊形單元的所有頂點作為壓力節點;用另一組包含表示法向流量值的網格邊的中心點作為流量節點.其中心思想是通過將流量節點設置在網格邊的中點,并作為常單元,以此滿足局部物質守恒.該方法可適用于任意形狀的網格剖分,下述推導以三角形網格為例,如圖1 所示,改進的格林元法具有壓力和流量兩套不同節點,其中壓力節點分布在三角形單元的三個頂點,流量節點分布在三角形單元三條邊的中心點.壓力節點分別標記為1,2,3,流量節點分別標記為a,b,c.

圖1 壓力和流量兩套節點三角形單元示意圖Fig.1 Sketch of pressure and flux two sets of nodes in triangle cell

對三角形單元獲取如式(4)的邊界積分方程,其包含有與壓力有關的積分項,需對壓力節點值進行插值來估計壓力函數在整個單元上的取值,通常壓力節點仍然選擇在多邊形網格的頂點,以保證插值函數之和在整個單元內的取值均為1.每個壓力節點相應的基函數表示為 φi,三角形單元內的壓力值可由節點值和基函數的加權平均得到,即

但對于每個流量節點,可根據法向流量的分段連續性,將網格邊上的法向通量視為一個常量,這是與原始格林元法截然不同的地方.流量節點所在的邊記作Γz(z=a,b,c),按上述思想,式(4)的離散格式為

式中

則單元內的局部方程組的矩陣形式為

計算矩陣B的逆矩陣,將上式變形為

利用權重系數θ控制顯隱式程度,得到

簡寫為

式中

3 改進格林元法與EDFM 的耦合求解

基質與裂縫之間的竄流量被認為是基質單元中的源項或匯項,因此有必要明確哪個單元中存在源項或匯項,而這和基質與裂縫之間的幾何性質有關.本文所提出的改進格林元法可適用于任意單元類型的網格剖分,如矩形單元和三角形單元,其中三角形單元的前處理,本文借助了開源程序Distmesh 三角形非結構化網格剖分的Matlab 代碼[30-31].相比于矩形結構化網格,三角形非結構化網格沒有規則的拓撲結構,網格節點的分布更具靈活性,在諸多工程力學領域均有較好的應用[32-33].如圖2 所示的二維EDFM,其中兩條離散裂縫分別以紅色和綠色的線段繪制,通過基質單元的網格邊界將與其相交的裂縫離散為多個裂縫單元,并對每個裂縫單元進行編號.

圖2 基于改進格林元法的二維EDFMFig.2 Two-dimensional EDFM based on the modified GEM

針對滲流介質某基質網格內含有裂縫單元的情況,可將源點取在裂縫網格中心,由于該點是基質網格內點,其特征角度為2π,則得到相應的邊界積分方程離散格式如下

此時Mfi是裂縫網格中心坐標.將式(12) 變形后,得到

式(14)即為考慮壓力瞬態效應的基質網格與裂縫單元之間的傳質表達式,式中B-1,R,E,C僅僅與基質網格和裂縫單元的幾何參數有關,表征了裂縫單元分布在基質網格中的幾何特征.與傳統線性分布假設不同,該表達式中除了包含基質孔隙壓力和裂縫孔隙壓力,還包含了裂縫孔隙壓力對時間的導數項,即體現了裂縫孔隙壓力的瞬時變化率.因此,相較于傳統兩點線性流量估計格式,可提高針對低滲透/致密儲層中由于基質與裂縫兩者間滲透率差異極大所產生的早期非穩態竄流量的數值計算精度.

對于相鄰的基質網格,存在連續性條件: 共有的壓力節點具有相同的壓力,共有的流量節點處的流量相同,即兩個相鄰基質網格在共有流量節點處的外法線流量的代數和等于0,這也是兩套節點格林元法的基本原理.如圖3 所示,以兩個相鄰三角形基質單元為例,介紹耦合過程的細節.

圖3 兩套節點格林元方法中相鄰單元方程組耦合示意圖Fig.3 Sketch of coupling adjacent element equations in two sets of nodes-based GEM

假設基質單元e1內的方程組為

單元e2內的方程組為

則兩個相鄰基質單元方程組耦合為

此外,對于離散裂縫網絡中的流體流動問題,可采用有限差分法對流動方程進行離散.如圖4 所示,單條裂縫可離散成多個裂縫單元,離散裂縫交匯處的流量交換可采用星三角變換方法[34-35]進行處理.

圖4 裂縫交匯的星三角變換原理圖Fig.4 Schematics of star-delta transformation for fracture segments intersecting

得到相應的有限差分隱式格式,共nf個方程

假設基質單元節點數為np、基質單元邊的個數為ned、裂縫單元的個數為nf.由此可知,全局方程組由兩部分構成: 一部分是基于式(11),采用式(17)耦合方法獲取的表征基質網格間滲流的方程組;第二部分是由式(18)構成的表征裂縫網格間流動的方程組,且這兩部分方程組中的qomf用式(14)表示.

整體上,全局方程組共有ne+nf個方程(其中第一部分方程有ne個,第二部分方程有nf個) 和np+nf個未知量(包含np個基質網格節點壓力和nf個來源于裂縫單元的未知量).因為ne大于np,可得知由改進格林元法獲得的該全局方程組為超定方程組,在封閉外邊界條件下,此時方程組可表示為

基于泛函分析中的正交投影定理,該超定方程組的解可等價于計算式(20)的解.由此,所有的未知量均可獲得.

4 模型驗證及討論

4.1 模型驗證

將本文提出的兩套節點格林元法EDFM 與原始EDFM[6]及商業模擬軟件tNavigator?LGR 模塊進行對比,其中tNavigator?是由俄羅斯Rock Flow Dynamics (RFD)公司開發的高性能商業模擬軟件,其結果可以參考為精確解.該實例是采用單水平井和單壓裂縫的機理模型驗證本文模型對長期產能預測的準確性,低滲透儲層物性參數值見表1.圖5 展示了矩形儲層域的示意圖,地層中心有一口生產井和一條壓裂210 m 的主裂縫,水平井筒穿過裂縫中點.對儲層域進行離散化,tNavigator?LGR 模塊局部加密網格處理的裂縫必須沿網格線方向展布,而在EDFM 中裂縫走向不受限制.水平井采用10 MPa 恒定井底流壓制度進行生產,圖6 展示了三種不同模型計算1000 天的壓力分布場圖,圖7展示了三種不同模型所獲得的日產油量曲線.結果表明,三種不同模型的結果吻合較好,驗證了該方法對產能預測的整體精度.

表1 數值算例的輸入參數Table 1 Input parameters of numerical case

圖5 壓裂水平井示意圖Fig.5 Sketch of fractured horizontal well

圖6 三種不同模型壓力分布對比圖(1000 d)Fig.6 Comparison of pressure distribution maps for three various simulators (1000 d)

圖6 三種不同模型壓力分布對比圖(1000 d)(續)Fig.6 Comparison of pressure distribution maps for three various simulators (1000 d) (continued)

圖7 三種不同模型產油速度的對比圖(1000 d)Fig.7 Comparison of oil rate for three various simulators (1000 d)

4.2 基質-裂縫間瞬態流動對比

該算例旨在驗證兩套節點格林元法EDFM 能夠較好地反映基質-裂縫間的瞬態流動,算例的基本屬性與表1 中的參數值相同,討論中考慮了兩個重要因素: 網格大小和基質滲透率,將網格尺寸設置為2 m × 2 m 和10 m × 10 m,基質滲透率分別為0.1 mD,1 mD 和10 mD.計算時間設置為10 d.為了實現對早期瞬態流動的準確描述,與上一個算例相比,在模擬過程中使用了更小的時間步長.圖8 展示了三種模型(網格步長10 m × 10 m)模擬10 d 壓力分布的比較,從壓力場的結果來看,不同模型之間沒有顯著差異.圖9 比較了在不同網格尺寸(2 m ×2 m 和10 m × 10 m)和基質滲透率(0.1 mD,1 mD,10 mD) 條件下,由三種模型(原始EDFM、本文EDFM 和tNavigator?)計算所得的日產油量曲線.與tNavigator?結果相比,從生產早期可以明顯看出,本文EDFM比原始EDFM 具有更高的精度,這是因為本文提出的兩套節點格林元法能夠有效地反映局部基質網格和裂縫單元之間的瞬態流動,取代了原始EDFM 的線性流動假設,瞬態效應的程度與網格尺寸和基質滲透率密切相關.平均相對誤差可由式(21) 獲得,計算結果如圖10 所示.結果表明,本文EDFM 的平均相對誤差比原始EDFM 小得多.同時,隨著基質滲透率的降低和網格尺寸的增大,瞬態效應的影響增大,平均相對誤差增大.本文所提出的兩套節點格林元法實現了對基質-裂縫間瞬態流動更準確的表征,能提高模擬的早期精度.

圖8 三種不同模型壓力分布對比圖(10 d)Fig.8 Comparison of pressure distribution maps for three various simulators (10 d)

圖9 井筒流量早期結果對比Fig.9 Comparison of early-time results for wellbore flow rate

圖10 原始EDFM 和修正EDFM 之間的誤差比較Fig.10 Error comparison between the original EDFM and themodified EDFM

式中,qi表示由EDFM 計算所得的井筒流量,m3/d;為由tNavigator?計算所得的井筒流量,m3/d;N為時間步數.

4.3 復雜油藏邊界-縫網-SRV 分區模型

原始EDFM 僅支持矩形結構化網格剖分,在網格適用性方面存在局限性.由于格林元法的引入,本文改進后的EDFM 可支持非匹配性的非結構化網格剖分,能較完善地考慮復雜油藏邊界-縫網及SRV 分區的非常規油氣儲層體積壓裂井產量分析的復雜問題.本算例考慮基質滲透率為0.1 mD 的致密儲層,人工裂縫滲透率為5000 mD,其他參數值同表1,圖11(a)展示了儲層域的基本模型,圖11(b)和圖11(c) 分別考慮了圓形和矩形的SRV 改造區.COMSOL 商業模擬軟件以標準有限元(SFEM)為基礎,對DFM 的模塊化處理極為成熟.本文以SFEMCOMSOL 數值解為參考,以圖11(a)中的基本模型為例,對本文模型進行驗證,不同模型的網格剖分對比見圖12.水平井采用10 MPa 恒定井底流壓制度進行生產,圖13 對比了本文模型與SFEM-COMSOL兩種模型模擬1000 d 的日產油量曲線,結果對比本文模型的解與參考解較為吻合,驗證了本文模型在考慮復雜油藏邊界條件下產量預測的準確性.

圖11 儲層域示意圖Fig.11 Sketch of reservoir domain

圖12 不同模型的網格剖分對比圖(續)Fig.12 Comparison diagram of mesh division for different simulators (continued)

圖13 修正EDFM 與COMSOL 產油速度的對比圖(1000 d)Fig.13 Comparison of oil rate between modified EDFM and COMSOL (1000 d)

圖14 展示了本文模型對圖11(a)的基本模型不同開發階段得到的壓力分布情況.圖15 展示了本文模型對圖11(b)和圖11(c)兩種不同SRV 分區模擬100 天的壓力分布對比,從圖中能看出開發早期以縫網內部流體流動為主,生產井產量完全由裂縫供給;當考慮SRV 分區時,開發中期離散裂縫周圍的內區流體以垂直于裂縫面的方向流入裂縫,由于內外區基質滲透率級別存在差異,開發后期壓力波才慢慢向外區傳播.本文所提出的方法采用非穩態滲流控制方程的邊界積分形式推導了基質網格與裂縫網格之間傳質量的新格式,因此也可以將EDFM 的適用性拓展至油藏壓裂井早期生產動態特征分析.圖16 對比了不考慮分區、矩形SRV 和圓形SRV三種SRV 形態的生產動態特征,其中矩形SRV 對應的外區流動形態為復合線性流,圓形SRV 對應的為擬徑向流,不考慮分區的產量動態特征位于二者之間.

圖14 基礎模型壓力分布場圖Fig.14 Diagram of pressure distribution field for basic model

圖15 不同SRV 分區模型的壓力分布對比(100 d)Fig.15 Comparison of pressure distribution of different SRV zoning models (100 d)

圖16 不同SRV 對應的產量動態特征對比Fig.16 Comparison of dynamic production characteristics corresponding to different SRV

5 結論

(1)本文提出了滿足局部物質守恒且具有高精度、適用于任意網格剖分的壓力流量兩套節點的格林元方法,能用于二階偏微分方程的穩定求解;將兩套節點格林元法與EDFM 進行耦合求解,能更精確地獲得瞬態壓力和流量分布,通過實例驗證了本文模型的準確性.

(2)針對原始EDFM 僅適用于矩形網格剖分、難于考慮復雜油藏邊界的局限性,本文提出的修正EDFM 能適應任意多邊形單元的網格剖分,可適用于復雜油藏邊界、縫網及SRV 分區等復雜滲流模型的計算,是對原始EDFM 在網格適用性方面的拓展和升級.