黃朝琴，高 博，王月英，嚴(yán) 俠，姚 軍

(中國(guó)石油大學(xué)石油工程學(xué)院,山東青島266580)

裂縫作為最小的地質(zhì)構(gòu)造廣泛存在于地殼中,其尺度從微米級(jí)裂縫到千米級(jí)大斷層,可跨越多個(gè)數(shù)量級(jí)。裂縫既可作為高導(dǎo)流通道,也可成為流動(dòng)屏障。裂縫性介質(zhì)流動(dòng)傳輸問題一直是石油工業(yè)、巖石水力學(xué)以及核廢料處理工程等研究領(lǐng)域的熱點(diǎn)和難點(diǎn)[1]。目前,主要有雙重介質(zhì)、等效連續(xù)介質(zhì)和離散裂縫模型3種流動(dòng)數(shù)學(xué)模型。雙重介質(zhì)模型由Barrenblatt等[2]最早提出,由 Warren和 Root、Kazami等[3-4]進(jìn)一步發(fā)展。該模型認(rèn)為裂縫性介質(zhì)中存在裂縫和基巖兩個(gè)平行的滲流系統(tǒng),裂縫系統(tǒng)為流動(dòng)通道,基巖系統(tǒng)為儲(chǔ)集空間。該模型刻畫出了裂縫優(yōu)先滲流的特點(diǎn)。隨后,Pruess和Wu等[5-6]在雙重介質(zhì)模型基礎(chǔ)上,對(duì)基巖系統(tǒng)進(jìn)行網(wǎng)格細(xì)分,提出了MINC(Multiple Interaction Continua)模型,有效地提高了計(jì)算精度和適用性[5-6],然而雙重介質(zhì)模型中竄流函數(shù)的確定仍存在較大困難,尤其是對(duì)于兩相和多相流問題[7-8]。與雙重介質(zhì)模型不同,等效連續(xù)介質(zhì)模型將介質(zhì)視為一個(gè)連續(xù)性系統(tǒng),通過等效參數(shù)來表征其非均質(zhì)性；該模型計(jì)算效率高、參數(shù)需求簡(jiǎn)單,在巖石水力學(xué)中得到了長(zhǎng)足發(fā)展[9]。目前,對(duì)該模型的研究大都仍局限于單相流,對(duì)于兩相或多相流問題,如何獲取相應(yīng)的等效參數(shù),尚未有成熟的理論和方法[10-11]。離散裂縫模型對(duì)介質(zhì)中的每條裂縫予以顯示表征,具有計(jì)算精度高、擬真性好的優(yōu)點(diǎn),但計(jì)算量大。筆者將離散裂縫網(wǎng)絡(luò)模型視為離散裂縫模型的一種,即忽略基巖滲透性[12]。近年來,隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展,基于該模型的精細(xì)流動(dòng)模擬已成為可能；同時(shí),該模型可用于求解雙重介質(zhì)和等效連續(xù)介質(zhì)模型的相關(guān)參數(shù)[7,13-14]?，F(xiàn)有的離散裂縫流動(dòng)數(shù)值計(jì)算方法主要有有限體積法和有限元法兩類。前者在裂縫交叉處須進(jìn)行簡(jiǎn)化和等效處理,導(dǎo)致在大規(guī)模計(jì)算時(shí)計(jì)算精度降低[15-17]；后者則在守恒型計(jì)算格式構(gòu)造和計(jì)算穩(wěn)定性方面存在一定缺陷[18-22]。模擬有限差分法(Mimetic Finite Difference,MFD)作為一種新型數(shù)值計(jì)算方法,因其良好的局部守恒性和對(duì)復(fù)雜網(wǎng)格的適用性,在計(jì)算流體力學(xué)、電磁場(chǎng)和油藏?cái)?shù)值模擬等研究中得到了成功應(yīng)用[23-25]。筆者將該方法進(jìn)一步推廣至離散裂縫模型流動(dòng)數(shù)值模擬研究中,詳細(xì)闡述模擬有限差分方法的基本原理,建立相應(yīng)的離散裂縫數(shù)值計(jì)算格式,并采用IMPES方法對(duì)其兩相流問題進(jìn)行求解,最后通過算例驗(yàn)證方法的正確性和程序的魯棒性。

1 兩相滲流數(shù)學(xué)模型

為簡(jiǎn)單起見,本文中僅考慮不可壓縮油水兩相滲流問題。采用經(jīng)典的分流量數(shù)學(xué)模型,其中壓力方程為

式中,v為滲流總速度；K為滲透率張量；λ為總流度系數(shù),并定義分流函數(shù)fl=λl/λ；ρl為兩相流體的密度；g為重力加速度；z為油藏深度(向上為正)。

定義全局壓力p為

式中,pc為毛管力；Sw為水相飽和度。

相應(yīng)的水相飽和度方程為

式中,φ為孔隙度。假設(shè)基巖和裂縫中的流動(dòng)均滿足Darcy定律,因此,上述方程在整個(gè)裂縫性介質(zhì)區(qū)域上均適用。

本文中采用IMPES(Implicit Pressure and Explicit Saturation Scheme)方案順序求解方程(1)和(3)。其中,壓力方程(1)采用模擬有限差分法進(jìn)行求解,飽和度方程(3)則采用有限體積法進(jìn)行顯示求解。

為適應(yīng)離散裂縫模型的復(fù)雜幾何形狀,采用非結(jié)構(gòu)化網(wǎng)格剖分技術(shù)對(duì)研究區(qū)域進(jìn)行離散,如圖1所示。由于裂縫開度較小,對(duì)裂縫予以降維處理,即在二維問題中裂縫簡(jiǎn)化為裂縫線元,三維問題中則簡(jiǎn)化為裂縫面元。通過降維處理,能夠大大減少網(wǎng)格數(shù)量,提高計(jì)算效率,而裂縫開度僅需在具體的數(shù)值計(jì)算中考慮。

圖1 離散裂縫模型及其非結(jié)構(gòu)化網(wǎng)格剖分示意圖Fig.1 Schematic of discrete fractured porous medium and its unstructured meshing

2 壓力方程求解

2.1 模擬有限差分方法

假設(shè)研究區(qū)域Ω∈Rd被一組互不重疊的多邊形(d=2)或多面體(d=3)網(wǎng)格Ωh={Ωi}剖分。如圖2所示,取任意單元Ωi進(jìn)行分析,其中Ωj為其相鄰單元,交界面,為交界面Ak的面積加權(quán)法向量,k為其單位外法線向量。

首先,在單元中心點(diǎn)xi和邊界面中心點(diǎn)xk上,分別定義單元壓力pi和邊界面壓力πk,

若考慮重力項(xiàng),則上述壓力的定義應(yīng)視作流動(dòng)勢(shì)。由Darcy定律知邊界面上的法向滲流速度v可寫成：

式中,Ti為傳導(dǎo)矩陣；m為單元Ωi的邊界面數(shù)。

圖2 模擬有限差分網(wǎng)格單元分析示意圖Fig.2 Schematic of grid analysis of mimetic finite difference method

矩陣Ti的構(gòu)造是MFD方法的關(guān)鍵。假設(shè)壓力在單元上呈線性變化,即p=ax+b,則由Darcy定律可得

結(jié)合方程(6)和(7),考慮到pi-πk=a(xi-xk),可得

在此定義x(i)表示x的第i維笛卡爾坐標(biāo),則有

即,其中Ed為d階單位矩陣,故通過方程(8)求得傳導(dǎo)矩陣Ti為

其中T2X=0,為保證Ti矩陣逆的存在性,可應(yīng)用Brezzi-Lipnikov-Simoncini定理[23]構(gòu)造T2矩陣,表達(dá)式為

對(duì)于方程(1)中的連續(xù)性方程,直接在單元Ωi上積分,結(jié)合散度定理可得

考慮單元邊界面上的速度連續(xù)條件,結(jié)合方程(6)和(13),可得到MFD數(shù)值計(jì)算格式為

式中,v=[vk]為單元邊界面滲流速度列陣；p=[pi]為單元中心壓力列陣；π=[πk]為單元邊界面中心壓力列陣；g=[gk]為重力作用項(xiàng)；q=[qi]為單元Ωi的源匯項(xiàng)；f=[fi]為流量邊界條件,f=0表示不滲透邊界；Ne為網(wǎng)格單元總數(shù)。

方程(14)的第一行對(duì)應(yīng)于方程(1)中的Darcy定律,第二行對(duì)應(yīng)于方程(1)中連續(xù)性方程,第三行則是單元邊界面上的法向速度連續(xù)性條件。從上述推導(dǎo)過程可知：MFD方法僅基于單個(gè)網(wǎng)格單元構(gòu)造數(shù)值計(jì)算格式,具有良好的局部守恒性,這一點(diǎn)類似于混合有限元法[21]；對(duì)于復(fù)雜網(wǎng)格系統(tǒng),混合有限元數(shù)值計(jì)算格式的構(gòu)造存在困難,而MFD方法原則上適用于任意復(fù)雜網(wǎng)格,具有較大優(yōu)勢(shì)。

2.2 離散裂縫模型壓力方程求解

裂縫和基巖中的流動(dòng)均滿足Darcy定律,考慮封閉外邊界,則相應(yīng)的方程(14)變?yōu)?/p>

式中,下標(biāo)m和f分別表示基巖和裂縫。

本文中對(duì)裂縫進(jìn)行了降維處理,因此方程(16)比方程(15)的空間維數(shù)低一維。MFD離散裂縫數(shù)值計(jì)算格式的構(gòu)造,關(guān)鍵在于基巖和裂縫壓力方程的耦合。圖3所示為裂縫—基巖混合網(wǎng)格示意圖。

圖3 裂縫—基巖耦合流動(dòng)分析示意圖Fig.3 Schematic of flow analysis between fracture cell and matrix cells

由于裂縫網(wǎng)格單元可視為基巖網(wǎng)格單元的邊界面,裂縫單元壓力pf與相鄰基巖單元的邊界面壓力πm相等,因此在數(shù)值計(jì)算格式中僅需保留πm。對(duì)于方程(15)和(16)中的滲流速度項(xiàng),按照下述條件,在裂縫單元F=E∩E′上進(jìn)行耦合。

(1)若F為導(dǎo)流裂縫,則將裂縫單元與相鄰基巖單元間總的流量交換記為。對(duì)于裂縫單元,該流量可作為源匯項(xiàng)來處理,因此有

方程(17)的第二行對(duì)應(yīng)于裂縫單元的連續(xù)性方程。

(2)若F為流動(dòng)屏障,則按照不滲透邊界進(jìn)行處理。此時(shí),方程(15)和(16)可通過方程(17)耦合起來形成相應(yīng)的離散裂縫數(shù)值計(jì)算格式,

3 飽和度方程求解

3.1 有限體積法計(jì)算格式

本文中采用IMPES方法對(duì)方程(1)和(3)進(jìn)行順序求解。首先,應(yīng)用方程(18)求解離散裂縫模型的壓力方程(1),計(jì)算中,與飽和度相關(guān)的參數(shù)均取上一個(gè)時(shí)間步的數(shù)值；對(duì)于飽和度方程(3),則采用有限體積法進(jìn)行求解,在單元上直接對(duì)方程(3)進(jìn)行積分可得

為書寫方便,在此去掉了水相飽和度Sw的下標(biāo)w。對(duì)于時(shí)間維,應(yīng)用θ準(zhǔn)則可得到下述有限體積數(shù)值離散格式：

式中,上標(biāo)n表示時(shí)間步。

對(duì)飽和度方程(20)進(jìn)行顯示求解,即θ=0。為達(dá)到計(jì)算的穩(wěn)定性,時(shí)間步長(zhǎng)采用如下CFL條件,

式中,S*為歸一化后的水相飽和度；Swc為束縛水飽和度；Sor為殘余油飽和度。

3.2 裂縫交叉處飽和度計(jì)算

當(dāng)兩條或多條裂縫相交時(shí),交叉處飽和度的計(jì)算是離散裂縫流動(dòng)模擬的關(guān)鍵。目前主要有兩種處理方式：一種是基于Delta-Star傳導(dǎo)率計(jì)算的上游迎風(fēng)格式[15],該方法對(duì)交叉裂縫進(jìn)行簡(jiǎn)化和等效處理；另一種則是上游迎風(fēng)加權(quán)計(jì)算格式[20],該方法計(jì)算精度高,但需獲取交叉處各裂縫單元的真實(shí)滲流速度。本文中采用后者,如圖4所示,假設(shè)有NI個(gè)裂縫單元ei相交于I,每個(gè)裂縫單元相應(yīng)的分流函數(shù)為fw,ei,在交叉處的滲流速度為vf,ei；定義I處的流入和流出為

由質(zhì)量守恒定律可知

由上游迎風(fēng)計(jì)算格式定義,可得

因此,裂縫交叉處I的上游迎風(fēng)加權(quán)分流函數(shù)為

圖4 裂縫交叉處飽和度計(jì)算示意圖Fig.4 Schematics of saturation analysis for intersecting fractures

4 數(shù)值算例

首先給出兩個(gè)簡(jiǎn)單離散裂縫模型數(shù)值算例,并通過與實(shí)驗(yàn)結(jié)果的對(duì)比驗(yàn)證方法和程序的正確性；然后,通過復(fù)雜離散裂縫模型算例進(jìn)一步驗(yàn)證該方法的正確性和程序的魯棒性。

4.1 簡(jiǎn)單離散裂縫模型算例

考慮圖5、6所示的一注一采物理模型,其尺寸為1 m×1 m×0.025 m,可視為平面流動(dòng)問題。圖5為單條裂縫模型,圖6為兩條交叉裂縫模型,均采用玻璃砂(0.100～0.104 mm)結(jié)合環(huán)氧樹脂膠結(jié)壓實(shí)而成,然后由透明有機(jī)玻璃板封裝?；鶐r可視為均質(zhì)各向同性介質(zhì),其孔隙度φ≈0.4,滲透率km=10 μm2。模型制作時(shí),裂縫由超薄鋼片替代,待模型膠結(jié)后抽離,開度約為1 mm,滲透率kf=a2/12=8.33×104μm2。注水井的流量qin=0.01VP/min,采出井與大氣壓相連,其中VP表示孔隙體積。水的黏度μw=1 mPa·s,油的黏度μo=5 mPa·s,水的密度ρw=1000 kg/m3,油的密度ρo=800 kg/m3。

模型初始時(shí)刻飽和油,束縛水飽和度和殘余油飽和度均為零?；鶐r和裂縫的油相相對(duì)滲透率Kro=1-Sw,水相相對(duì)滲透率Krw=Sw。應(yīng)用本文中方法對(duì)上述兩個(gè)物理實(shí)驗(yàn)過程進(jìn)行了數(shù)值模擬,計(jì)算中忽略毛管力和重力的影響,相應(yīng)的Delaunay三角網(wǎng)格剖分和數(shù)值模擬結(jié)果如圖5和6所示。通過與試驗(yàn)流動(dòng)過程的對(duì)比可看到：數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果基本一致,從而驗(yàn)證了本文方法和程序的正確性。值得注意的是,圖5中實(shí)驗(yàn)結(jié)果沿左邊界的上半部分出現(xiàn)了快速流動(dòng)現(xiàn)象,這是由于在模型制作時(shí)側(cè)面的有機(jī)玻璃板和玻璃砂膠結(jié)模型并未達(dá)到完全密封造成的。

圖5 單裂縫模型及其含水飽和度結(jié)果的對(duì)比Fig.5 Single discrete fracture model and comparison of water saturation profiles between numerical results and experimental results

圖6 兩條交叉裂縫模型及其含水飽和度結(jié)果對(duì)比Fig.6 Intersecting discrete fractures model and comparison of water saturation profiles between numerical results and experimental results

4.2 復(fù)雜離散裂縫模型算例1

研究圖7所示的復(fù)雜裂縫性油藏模型,其中包含導(dǎo)流裂縫和淤泥質(zhì)充填斷層,油藏厚度為單位長(zhǎng)度。非均質(zhì)各向異性基巖的孔隙度φ=0.2,滲透率分布如圖7(b)所示,裂隙開度a=1 mm,滲透率kf=a2/12=8.33×107μm2。水的黏度μw=1 mPa·s,油的黏度μn=5 mPa·s,束縛水飽和度為0,殘余油飽和度為0.2。水相相對(duì)滲透率Krw=S2w,油相相對(duì)滲透率Krn=(1-Sw)2,初始?jí)毫?0 MPa,注水和采油端的速度均為qin=30 m3/d。

圖7 復(fù)雜裂縫性油藏模型Fig.7 Fractured reservoir within complex discrete fracture

假設(shè)模型為水濕性油藏儲(chǔ)層,考慮基巖和裂縫中毛管力的影響,假設(shè)兩者的毛管力均符合Brooks-Corey毛管力函數(shù)。對(duì)于基巖,閾壓值pd=10 kPa,λ取值2.0；對(duì)于裂縫,閾壓值pd=1 kPa,λ取值1.0,表達(dá)式為

如圖8所示,分別采用兩種不同計(jì)算網(wǎng)格系統(tǒng)進(jìn)行數(shù)值模擬,其中三角形網(wǎng)格包含881個(gè)節(jié)點(diǎn)和1615個(gè)單元,四邊形網(wǎng)格包含1859個(gè)節(jié)點(diǎn)1737個(gè)單元。圖8給出了不同時(shí)刻的含水飽和度分布,結(jié)果表明兩種不同網(wǎng)格的計(jì)算結(jié)果基本吻合,進(jìn)一步驗(yàn)證了方法的正確性和程序的魯棒性。同時(shí),基于三角網(wǎng)格系統(tǒng)應(yīng)用混合有限元方法進(jìn)行了模擬,MFD計(jì)算時(shí)間為356.13 s,混合有限元?jiǎng)t為351.21 s,因此兩者的計(jì)算量相當(dāng),但MFD對(duì)網(wǎng)格的要求更低,更適用于復(fù)雜離散裂縫模型的研究。

圖8 基于不同網(wǎng)格系統(tǒng)的含水飽和度計(jì)算結(jié)果Fig.8 Water saturation simulation results of different grid systems

4.3 復(fù)雜離散裂縫模型算例2

考慮圖9所示的復(fù)雜裂縫性介質(zhì)模型,研究區(qū)域?yàn)?00 m×50 m(x×y),圖中藍(lán)線代表裂縫,該模型基于某油田露頭資料生成。均質(zhì)各向同性基巖的孔隙度φ=0.2,滲透率km=10×10-3μm2,裂縫開度a=1 mm,滲透率kf=a2/12=8.33×104μm2,油水的物性參數(shù)與4.2算例一致。初始油藏壓力為10 MPa,初始含水飽和度為零,注水井和采油井的速度均為0.01VP/d。

圖9 復(fù)雜離散裂縫模型及其非結(jié)構(gòu)化網(wǎng)格剖分Fig.9 Complex discrete fracture model and corresponding Delaunay triangle meshing

基巖和裂縫的油相相對(duì)滲透率Kro=(1-Sw)2,水相相對(duì)滲透率Krw=Sw2,假設(shè)模型為水濕性儲(chǔ)層,考慮基巖和裂縫中毛管力的影響,假設(shè)兩者的毛管力均符合Brooks—Corey毛管力函數(shù),見式(25),相關(guān)參數(shù)與上述4.2算例相同。計(jì)算中采用如圖9所示的Delaunay三角網(wǎng)格。圖10給出了不同時(shí)刻的含水飽和度分布,其數(shù)值計(jì)算結(jié)果表明：注入水沿著大裂縫迅速竄進(jìn),裂縫的存在導(dǎo)致了介質(zhì)的強(qiáng)烈非均質(zhì)性；由于毛管力的存在,注水波及面積有所增加,但整體效果仍由宏觀大裂縫控制。

圖10 不同時(shí)刻含水飽和度分布Fig.10 Water saturation profiles at different time

5 結(jié)束語

基于模擬有限差分方法,建立了一套新的離散裂縫數(shù)值計(jì)算格式,并采用IMPES方法對(duì)其兩相滲流問題進(jìn)行了研究,通過算例驗(yàn)證了方法的正確性和程序的魯棒性。模擬有限差分法在構(gòu)造數(shù)值計(jì)算格式時(shí),僅基于單個(gè)網(wǎng)格的節(jié)點(diǎn)和面信息,具有良好的局部守恒性,理論上適用于任何復(fù)雜網(wǎng)格系統(tǒng),適用于離散裂縫流動(dòng)模擬研究。本文中僅對(duì)二維兩相流問題進(jìn)行了討論,對(duì)三維三相問題進(jìn)行研究時(shí)在模型中應(yīng)考慮非線性流動(dòng)問題,以適用于低滲透油藏和致密油氣藏開發(fā)的需求。

[1] National Research Council.Rock fractures and fluid flow： contemporary understanding and applications[M].Washington,DC：National Academies Press,1996：1-10.

[2] BARENBLATT G,ZHELTOV I P,KOCHINA I.Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks[J].Journal of Applied Mathematics and Mechanics,1960,24(5)：1286-1303.

[3] WARREN J,ROOT P J.The behavior of naturally fractured reservoirs[J].Old SPE Journal,1963,3(3)：245-255.

[4] KAZEMI H,LS M,PORTERFIELD K,et al.Numerical simulation of water-oil flow in naturally fractured reservoirs[J].Old SPE Journal,1976,16(6)：317-326.

[5] PRUESS K.A practical method for modeling fluid and heat flow in fractured porous media[J].Old SPE Journal,1985,25(1)：14-26.

[6] WU Y-S,PRUESS K.A multiple-porosity method for simulation of naturally fractured petroleum reservoirs[J].SPE Reservoir Engineering,1988,3(1)：327-336.

[7] KARIMI-FARD M,GONG B,DURLOFSKY L.Generation of coarse-scale continuum flow models from detailed fracture characterizations[J].Water Resources Research,2006,42(10)：W10423.

[8] BOURBIAUX B.Fractured reservoir simulation：a challenging and rewarding issue[J].Oil&Gas Science and Technology,2010,65(2)：227-238.

[9] 周志芳.裂隙介質(zhì)水動(dòng)力學(xué)原理[M].北京：高等教育出版社,2007.

[10] 王月英,姚軍,黃朝琴.裂隙巖體流動(dòng)模型綜述[J].大慶石油學(xué)院學(xué)報(bào),2011,35(5)：42-48.WANG Yueying,YAO Jun,HUANG Zhaoqin.A review of fluid flow models in fractured rock[J].Journal of Daqing Petroleum Institute,2011,35(5)：42-48.

[11] HUANG Z-Q,YAO J,WANG Y-Y.An efficient numerical model for immiscible two-phase flow in fractured Karst reservoirs[J].Communication in Computational Physics,2013,13(2)：540-558.

[12] HUANG Z,YAO J,WANG Y,et al.Numerical study on two-phase flow through fractured porous media[J].Science China Technological Sciences,2011,54(9)：2412-2420.

[13] GONG B,KARIMI-FARD M,DURLOFSKY L.Upscaling discrete fracture characterizations to dual-porosity,dual-permeability models for efficient simulation of flow with strong gravitational effects[J].SPE Journal,2008,13(1)：58-67.

[14] HUI M-H,MALLISON B,FYROZJAEE M,et al.The upscaling of discrete fracture models for faster,coarsescale simulations of IOR and EOR processes for fractured reservoirs[R].SPE 166075-MS,2013.

[15] KARIMI-FARD M,DURLOFSKY L,AZIZ K.An efficient discrete-fracture model applicable for general-purpose reservoir simulators[J].SPE Journal,2004,9(2)：227-236.

[16] SANDVE T,BERRE I,NORDBOTTEN J M.An efficient multi-point flux approximation method for discrete fracture—matrix simulations[J].Journal of Computational Physics,2012,231(9)：3784-3800.

[17] PANFILI P,COMINELLI A,SCOTTI A.Using embed-ded discrete fracture models(EDFMs)to simulate realistic fluid flow problems[C]//Second Workshop on Naturally Fractured Reservoirs,8-11 December,2013,Muscat Oman,European Association of Geoscientists&Engineering,c2013.

[18] KIM J G,DEO M D.Finite element,discrete-fracture model for multiphase flow in porous media[J].AIChE Journal,2000,46(6)：1120-1130.

[19] KARIMI-FARD M,FIROOZABADI A.Numerical simulation of water injection in fractured media using the discrete-fracture model and the Galerkin method[J].SPE Reservoir Evaluation&Engineering,2003,6(2)：117-126.

[20] MONTEAGUDO J E P,FIROOZABADI A.Control-volume method for numerical simulation of two-phase immiscible flow in two-and three-dimensional discrete-fractured media[J].Water Resources Research,2004,40(7)：W07405.

[21] HOTEIT H,FIROOZABADI A.An efficient numerical model for incompressible two-phase flow in fractured media[J].Advances in Water Resources,2008,31(6)：891-905.

[22] 黃朝琴,姚軍,王月英,等.基于離散裂縫模型的裂縫性油藏注水開發(fā)數(shù)值模擬[J].計(jì)算物理,2011,28(1)：41-49.HUANG Zhaoqin,YAO Jun,WANG Yueying,et al.Numerical simulation on water flooding development of fractured reservoirs in a discrete-fracture model[J].Chinese Journal of Computational Physics,2011,28(1)：41-49.

[23] BREZZI F,LIPNIKOV K,SIMONCINI V.A family of mimetic finite difference methods on polygonal and polyhedral meshes[J].Mathematical Models and Methods in Applied Sciences,2005,15(10)：1533-1551.

[24] LIE K A,KROGSTAD S,LIGAARDEN I S,et al.Open-source MATLAB implementation of consistent discretisations on complex grids[J].Computational Geosciences,2012,16(2)：297-322.

[25] LIPNIKOV K,MANZINI G,SHASHKOV M.Mimetic finite difference method[J].Journal of Computational Physics,2014,257：1163-1227.