兩套節(jié)點(diǎn)格林元嵌入式離散裂縫模型數(shù)值模擬方法1)

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

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

? (長(zhǎng)江大學(xué)石油工程學(xué)院,武漢 430100)

引言

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

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

本文從原始EDFM 的適用局限性出發(fā),利用格林元方法在處理非定常和非均質(zhì)問(wèn)題方面具有高精度的獨(dú)特優(yōu)勢(shì),建立一種能與EDFM 合理匹配的改進(jìn)格林元法,能精確求解基質(zhì)與裂縫間的非穩(wěn)態(tài)竄流問(wèn)題,拓展原始EDFM 技術(shù)應(yīng)用層面的局限性,使其既能求解各類復(fù)雜油藏滲流問(wèn)題,又能使求解結(jié)果具有明確的物理意義.

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

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

式中,p是計(jì)算域內(nèi)基質(zhì)的孔隙流體壓力,MPa;K和c為滲流介質(zhì)的屬性參數(shù)函數(shù),均為標(biāo)量;f是計(jì)算域內(nèi)的源匯項(xiàng)強(qiáng)度,若滲流介質(zhì)某基質(zhì)網(wǎng)格內(nèi)含有裂縫單元,則有,其中qomf表示該基質(zhì)網(wǎng)格與裂縫單元間的傳質(zhì)速度,m3/d;V是基質(zhì)網(wǎng)格體積,m3.

上式可以改寫為

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

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

式中,Mi為空間中的源點(diǎn)坐標(biāo),M為空間中任意一點(diǎn)的坐標(biāo),δ(M,Mi) 為Delta 函數(shù).

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

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

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

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

2 壓力流量?jī)商坠?jié)點(diǎn)格林元法

為此,本文提出了一種壓力流量?jī)商坠?jié)點(diǎn)的改進(jìn)格林元法,用一組節(jié)點(diǎn)包含表示壓力值的多邊形單元的所有頂點(diǎn)作為壓力節(jié)點(diǎn);用另一組包含表示法向流量值的網(wǎng)格邊的中心點(diǎn)作為流量節(jié)點(diǎn).其中心思想是通過(guò)將流量節(jié)點(diǎn)設(shè)置在網(wǎng)格邊的中點(diǎn),并作為常單元,以此滿足局部物質(zhì)守恒.該方法可適用于任意形狀的網(wǎng)格剖分,下述推導(dǎo)以三角形網(wǎng)格為例,如圖1 所示,改進(jìn)的格林元法具有壓力和流量?jī)商撞煌?jié)點(diǎn),其中壓力節(jié)點(diǎn)分布在三角形單元的三個(gè)頂點(diǎn),流量節(jié)點(diǎn)分布在三角形單元三條邊的中心點(diǎn).壓力節(jié)點(diǎn)分別標(biāo)記為1,2,3,流量節(jié)點(diǎn)分別標(biāo)記為a,b,c.

圖1 壓力和流量?jī)商坠?jié)點(diǎn)三角形單元示意圖Fig.1 Sketch of pressure and flux two sets of nodes in triangle cell

對(duì)三角形單元獲取如式(4)的邊界積分方程,其包含有與壓力有關(guān)的積分項(xiàng),需對(duì)壓力節(jié)點(diǎn)值進(jìn)行插值來(lái)估計(jì)壓力函數(shù)在整個(gè)單元上的取值,通常壓力節(jié)點(diǎn)仍然選擇在多邊形網(wǎng)格的頂點(diǎn),以保證插值函數(shù)之和在整個(gè)單元內(nèi)的取值均為1.每個(gè)壓力節(jié)點(diǎn)相應(yīng)的基函數(shù)表示為 φi,三角形單元內(nèi)的壓力值可由節(jié)點(diǎn)值和基函數(shù)的加權(quán)平均得到,即

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

式中

則單元內(nèi)的局部方程組的矩陣形式為

計(jì)算矩陣B的逆矩陣,將上式變形為

利用權(quán)重系數(shù)θ控制顯隱式程度,得到

簡(jiǎn)寫為

式中

3 改進(jìn)格林元法與EDFM 的耦合求解

基質(zhì)與裂縫之間的竄流量被認(rèn)為是基質(zhì)單元中的源項(xiàng)或匯項(xiàng),因此有必要明確哪個(gè)單元中存在源項(xiàng)或匯項(xiàng),而這和基質(zhì)與裂縫之間的幾何性質(zhì)有關(guān).本文所提出的改進(jìn)格林元法可適用于任意單元類型的網(wǎng)格剖分,如矩形單元和三角形單元,其中三角形單元的前處理,本文借助了開(kāi)源程序Distmesh 三角形非結(jié)構(gòu)化網(wǎng)格剖分的Matlab 代碼[30-31].相比于矩形結(jié)構(gòu)化網(wǎng)格,三角形非結(jié)構(gòu)化網(wǎng)格沒(méi)有規(guī)則的拓?fù)浣Y(jié)構(gòu),網(wǎng)格節(jié)點(diǎn)的分布更具靈活性,在諸多工程力學(xué)領(lǐng)域均有較好的應(yīng)用[32-33].如圖2 所示的二維EDFM,其中兩條離散裂縫分別以紅色和綠色的線段繪制,通過(guò)基質(zhì)單元的網(wǎng)格邊界將與其相交的裂縫離散為多個(gè)裂縫單元,并對(duì)每個(gè)裂縫單元進(jìn)行編號(hào).

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

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

此時(shí)Mfi是裂縫網(wǎng)格中心坐標(biāo).將式(12) 變形后,得到

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

對(duì)于相鄰的基質(zhì)網(wǎng)格,存在連續(xù)性條件: 共有的壓力節(jié)點(diǎn)具有相同的壓力,共有的流量節(jié)點(diǎn)處的流量相同,即兩個(gè)相鄰基質(zhì)網(wǎng)格在共有流量節(jié)點(diǎn)處的外法線流量的代數(shù)和等于0,這也是兩套節(jié)點(diǎn)格林元法的基本原理.如圖3 所示,以兩個(gè)相鄰三角形基質(zhì)單元為例,介紹耦合過(guò)程的細(xì)節(jié).

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

假設(shè)基質(zhì)單元e1內(nèi)的方程組為

單元e2內(nèi)的方程組為

則兩個(gè)相鄰基質(zhì)單元方程組耦合為

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

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

得到相應(yīng)的有限差分隱式格式,共nf個(gè)方程

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

整體上,全局方程組共有ne+nf個(gè)方程(其中第一部分方程有ne個(gè),第二部分方程有nf個(gè)) 和np+nf個(gè)未知量(包含np個(gè)基質(zhì)網(wǎng)格節(jié)點(diǎn)壓力和nf個(gè)來(lái)源于裂縫單元的未知量).因?yàn)閚e大于np,可得知由改進(jìn)格林元法獲得的該全局方程組為超定方程組,在封閉外邊界條件下,此時(shí)方程組可表示為

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

4 模型驗(yàn)證及討論

4.1 模型驗(yàn)證

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

表1 數(shù)值算例的輸入?yún)?shù)Table 1 Input parameters of numerical case

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

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

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

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

4.2 基質(zhì)-裂縫間瞬態(tài)流動(dòng)對(duì)比

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

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

圖9 井筒流量早期結(jié)果對(duì)比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 計(jì)算所得的井筒流量,m3/d;為由tNavigator?計(jì)算所得的井筒流量,m3/d;N為時(shí)間步數(shù).

4.3 復(fù)雜油藏邊界-縫網(wǎng)-SRV 分區(qū)模型

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

圖11 儲(chǔ)層域示意圖Fig.11 Sketch of reservoir domain

圖12 不同模型的網(wǎng)格剖分對(duì)比圖(續(xù))Fig.12 Comparison diagram of mesh division for different simulators (continued)

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

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

圖14 基礎(chǔ)模型壓力分布場(chǎng)圖Fig.14 Diagram of pressure distribution field for basic model

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

圖16 不同SRV 對(duì)應(yīng)的產(chǎn)量動(dòng)態(tài)特征對(duì)比Fig.16 Comparison of dynamic production characteristics corresponding to different SRV

5 結(jié)論

(1)本文提出了滿足局部物質(zhì)守恒且具有高精度、適用于任意網(wǎng)格剖分的壓力流量?jī)商坠?jié)點(diǎn)的格林元方法,能用于二階偏微分方程的穩(wěn)定求解;將兩套節(jié)點(diǎn)格林元法與EDFM 進(jìn)行耦合求解,能更精確地獲得瞬態(tài)壓力和流量分布,通過(guò)實(shí)例驗(yàn)證了本文模型的準(zhǔn)確性.

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