， ， ， ，

(中國石油大學(xué)(華東) 石油工程學(xué)院，青島 266580)

1 引 言

多孔介質(zhì)中精確有效的流動模擬對于了解地下資源的動態(tài)，實(shí)施有效管理有著重要的意義，如地下水、地?zé)崮茉醇绊搸r氣開發(fā)項目等。這些大規(guī)模的地質(zhì)構(gòu)造，除了固有的非均質(zhì)性，往往還包含復(fù)雜的天然裂縫。據(jù)不完全統(tǒng)計，世界上裂縫性油藏約占已探明總儲量的一半；在我國，已探明的裂縫性油藏地質(zhì)儲量超過40億噸，占探明總儲量的28%以上[1]。而致密油氣和頁巖氣等非常規(guī)油氣資源，需壓裂后方可商業(yè)開采[2]，壓裂后也為裂縫性油氣藏。

裂縫對流動形態(tài)有重要影響，因此裂縫的描述直接影響最終模擬結(jié)果的正確性。目前，主要采用雙重介質(zhì)模型、等效連續(xù)介質(zhì)模型和離散裂縫模型作為裂縫性介質(zhì)的流動模型。雙重介質(zhì)模型[3,4]是目前廣泛應(yīng)用的流動模型，但當(dāng)存在數(shù)條控制流體流動方向和規(guī)模的大裂縫時，其計算結(jié)果誤差較大。離散裂縫模型顯示表征介質(zhì)中的裂縫?，F(xiàn)有的離散裂縫數(shù)值模型都是基于匹配型網(wǎng)格，由于裂縫分布復(fù)雜，往往基于非結(jié)構(gòu)網(wǎng)格劃分，其剖分過程繁瑣，計算量大。尤其當(dāng)裂縫相距較近時，網(wǎng)格剖分質(zhì)量較差，導(dǎo)致計算誤差。對此，Lee等[5-8]提出了嵌入式離散裂縫模型，對基巖直接進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分，然后將裂縫嵌入基巖網(wǎng)格系統(tǒng)中，并根據(jù)裂縫與基巖的相交情況形成裂縫網(wǎng)格，避免了上述復(fù)雜的非結(jié)構(gòu)化網(wǎng)格剖分過程，大幅降低了計算復(fù)雜度。

然而，現(xiàn)有的嵌入式離散裂縫模型是基于有限差分法，當(dāng)裂縫介質(zhì)中裂縫分布復(fù)雜時，往往需要進(jìn)行非結(jié)構(gòu)化網(wǎng)格剖分，此時有限差分法并不適用。而且在規(guī)模較大的裂縫性油藏中，裂縫具有強(qiáng)烈多尺度性，為了得到精細(xì)解，需要進(jìn)行較精細(xì)的網(wǎng)格劃分，即便采用嵌入式離散裂縫模型，仍然可達(dá)數(shù)百萬甚至數(shù)億個網(wǎng)格單元。如果基于有限元法或有限差分法等傳統(tǒng)數(shù)值方法求解此模型，巨大的計算量將超出當(dāng)今計算機(jī)的計算能力。當(dāng)前尺度升級法[9]廣泛用于降低計算量，但是此方法無法充分利用小尺度精細(xì)數(shù)據(jù)信息，模擬精度不高。為此，需要一種既有尺度升級法的計算效率，又有較高計算精度的方法，多尺度方法[10-13]應(yīng)運(yùn)而生。

多尺度方法在粗網(wǎng)格上求解控制方程，在細(xì)網(wǎng)格上求解局部流動方程獲得多尺度基函數(shù)，以捕捉小尺度流動特征。多尺度方法最初用于求解橢圓方程，其后應(yīng)用于油藏數(shù)值模擬領(lǐng)域[14-16]，并隨著該方法的推廣，其逐漸應(yīng)用于裂縫性油藏數(shù)值模擬[17,18]。且多尺度有限體積法[19]還用于求解嵌入式離散裂縫模型。但是多尺度有限體積法需要構(gòu)建雙網(wǎng)格結(jié)構(gòu)，小尺度映射過程復(fù)雜。本文在多尺度模擬有限差分法的框架下，使用模擬有限差分法求解嵌入式離散裂縫模型的多尺度基函數(shù)，避免了使用復(fù)雜的雙網(wǎng)格結(jié)構(gòu)。本文闡述了多尺度算法的基本原理，建立了嵌入式離散裂縫模型多尺度基函數(shù)的模擬有限差分計算格式，最后通過多尺度解和參考解的對比驗(yàn)證了方法的正確性和程序的魯棒性。

2 嵌入式離散裂縫模擬

2.1 數(shù)學(xué)模型

假設(shè)流體等溫滲流，基巖和流體不可壓縮，且不考慮重力和毛管力的影響。流體在裂縫和基巖中的流動滿足Darcy定律:

基巖系統(tǒng)

-·[(Km)·pm]=qm+(qm f/Vm)δm f

(1)

裂縫系統(tǒng)

(2)

裂縫與基巖間的竄流量可以寫為

qm f=-Tm f(pm-pf)

(3)

當(dāng)裂縫相交時,裂縫間的竄流量可以寫為

qf f=Tf f(pf i-pf j)

(4)

2.2 數(shù)值離散

基巖區(qū)域使用模擬有限差分法求解[20],該方法適用于任何復(fù)雜網(wǎng)格的計算,而且擁有良好的局部守恒性。由Darcy定律可知,對于任一個網(wǎng)格單元,其邊界上的法向速度可以寫為

vi=Ti·(eipi-πi)

(5)

式中Ti為傳導(dǎo)矩陣;vi=[v1,…,vm]T,m為網(wǎng)格單元界面數(shù);ei=[1,…,1]T;pi為單元壓力;πi為邊界面壓力。裂縫系統(tǒng)采用隱式差分求解,由方程(2)可得

Tξ i +1/2(pf i +1-pf i)-Tξ i -1/2(pf i-pf i -1)=

ffi+qm f i+qf f iδf f i

(6)

ff i=Vf iqf i

將基巖系統(tǒng)和裂縫系統(tǒng)耦合在一起可以得到嵌入式離散裂縫模型的模擬有限差分?jǐn)?shù)值計算格式：

(7)

式中Tm f i=[Tm f i]為第i條裂縫與基巖竄流系數(shù)矩陣,Tf f=[Tf f]為裂縫之間的竄流系數(shù)矩陣，Tf i和pfi分別為第i條裂縫的有限差分傳導(dǎo)系數(shù)矩陣和裂縫單元壓力列陣。

3 多尺度計算格式

多尺度計算方法包含兩套網(wǎng)格系統(tǒng),如圖1所示。其中,細(xì)網(wǎng)格系統(tǒng)包含巖石及流體性質(zhì)等參數(shù),粗網(wǎng)格單元由相互連接的細(xì)網(wǎng)格組成。

圖1 多尺度粗網(wǎng)格示意圖

Fig.1 Schematic of the coarse gridcells

(8)

(9)

(10)

式中σ(x)=trace[K(x)]/d。將速度基函數(shù)表示為局部區(qū)域速度值的列向量,將壓力基函數(shù)表示為局部區(qū)域壓力值的列向量。

當(dāng)粗網(wǎng)格單元內(nèi)包含裂縫時,采用嵌入式離散裂縫模型表征裂縫單元。通過求解方程(11，12)可以獲得多尺度基函數(shù),捕捉裂縫及基巖流動特征。

(11)

(12)

計算多尺度基函數(shù)后,為了得到粗網(wǎng)格方程,要將局部區(qū)域的多尺度基函數(shù)分為兩部分：

Ψi j=ΨHi j-ΨHj i

(13)

滿足

(14)

多尺度方法的一個顯著特點(diǎn)就是易于進(jìn)行小尺度映射,即細(xì)網(wǎng)格上的速度和壓力可由對應(yīng)的多尺度基函數(shù)經(jīng)過拓展得來。對于小尺度速度,vf≈Ψvc,對于小尺度壓力,

pf=Ip+ΦDλvc

(15)

粗網(wǎng)格表面壓力πc可表示為

(16)

由式(15,16)可得粗網(wǎng)格方程為

(17)

式中Bc=ΨTBfΨ,Cc=ΨTCfI,fc=ITqf,Dc=ΨTCfJ。求解粗網(wǎng)格方程之后,可通過式(15)得到小尺度解。

圖2 目標(biāo)粗網(wǎng)格單元的4個多尺度壓力基函數(shù)

Fig.2 Illustration of all four pressure basis functions of a coarse gridcell

4 數(shù)值算例

4.1 均質(zhì)嵌入式離散裂縫模型算例

考慮二維離散裂縫介質(zhì)包含兩條垂直相交的裂縫,如圖(3)所示。研究區(qū)域大小為1×1 m2。裂縫與基巖的滲透率比Kf/Km=105,基巖是均質(zhì)各向同性,孔隙度為0.2,裂縫開度為4 cm。區(qū)域上下邊界為不滲透邊界,左右邊界為定壓邊界。細(xì)網(wǎng)格劃分為40×40,粗網(wǎng)格劃分為10×10。

圖4(a)給出了40×40的網(wǎng)格劃分下使用模擬有限差分法得到的壓力分布,圖4(b)是10×10的粗網(wǎng)格下多尺度方法得到的壓力分布,圖4(c)是多尺度方法小尺度映射后的壓力分布。對比三個壓力分布可以看出,多尺度方法得到的壓力值與參考解呈現(xiàn)相同的分布形態(tài),小尺度映射后的壓力分布與參考解的壓力分布基本一致。同時多尺度方法的計算速度比傳統(tǒng)數(shù)值方法提高了3.4倍,大大節(jié)約了內(nèi)存,減少了計算時間。

圖3 裂縫模型及多尺度網(wǎng)格劃分示意圖

Fig.3 Illustration of the geometry and multiscale grids of the fracture model

4.2 非均質(zhì)嵌入式離散裂縫模型算例

算例2是一個非均質(zhì)介質(zhì),滲透率場如圖5(a)所示,包含兩條嵌入其中的裂縫,裂縫開度為 1 cm,網(wǎng)格劃分如圖5(b)所示。模型大小為10 m×5 m,小尺度細(xì)網(wǎng)格剖分包含40×20個網(wǎng)格,多尺度粗網(wǎng)格包含10×5個網(wǎng)格。區(qū)域上下邊界為不滲透邊界,左右邊界為定壓邊界,右端壓力為 1 MPa,左端壓力為0。

圖6(a)給出了40×20網(wǎng)格下的參考壓力分布,圖6(b)給出了20×10粗網(wǎng)格下多尺度解在小尺度映射后的壓力分布。可以看出,多尺度方法不僅可以反映裂縫的壓力,而且可以準(zhǔn)確反映介質(zhì)的非均質(zhì)性。圖7給出了介質(zhì)在y=2.25 m上的參考解和多尺度解的速度和壓力值對比。數(shù)值結(jié)果表明,多尺度方法在較粗的尺度上就可以獲得較高的精度,驗(yàn)證了多尺度方法在處理非均質(zhì)介質(zhì)時的正確性和魯棒性。

圖4 參考解和多尺度解的壓力分布

Fig.4 Pressure maps obtained by reference solution and multiscale method

圖5 非均質(zhì)油藏嵌入式離散裂縫模型的滲透率場及網(wǎng)格劃分

Fig.5 Permeability and grids in a 2D heterogeneous reservoir with an embedded fracture network

4.3 三維嵌入式離散裂縫模型算例

算例3是一個3D模型,如圖8所示,區(qū)域大小為10 m×10 m×4 m。均質(zhì)基巖滲透率為K=1 μm2,6條裂縫嵌入基巖中,裂縫開度均為1 cm。注入井位于左下角,生產(chǎn)井位于右上角。小尺度細(xì)網(wǎng)格剖分包含20×20×8個網(wǎng)格,多尺度粗網(wǎng)格包含5×5×4個網(wǎng)格。

圖9(a)給出了三維嵌入式離散裂縫模型的壓力分布參考解,圖9(b)是多尺度方法得到的壓力分布。采用相對L2范數(shù)表示壓力相對誤差：

ep=‖pf-pm s‖22/‖pf‖22

(18)

數(shù)值結(jié)果表明,多尺度方法與傳統(tǒng)數(shù)值方法得到的壓力相對誤差小于0.106,同時計算速度提高了 3.1倍。與傳統(tǒng)尺度升級方法不同,多尺度方法在大尺度上計算時,通過多尺度映射捕捉小尺度特征。圖9 證明了多尺度方法模擬三維裂縫性油藏的正確性和高效性,具有用于開發(fā)新一代數(shù)值模擬軟件的巨大潛力。

圖6 參考解和多尺度解的壓力分布

Fig.6 Pressure maps obtained by reference solution and multiscale method

圖7 參考解和多尺度解的壓力對比與速度對比

Fig.7 Velocity and pressure comparisons for fine -scale and MsMFEM solutions aty=2.25 m

圖8 三維裂縫模型及多尺度網(wǎng)格示意圖

Fig.8 Illustration of the geometry and multiscale grids of the fracture model

圖9 參考解和多尺度解的壓力分布

Fig.9 Pressure maps obtained by reference solution and multiscale method

5 結(jié) 論

(1) 離散裂縫模型顯示表征介質(zhì)中的每條裂縫,具有計算精度高、擬真性好的優(yōu)點(diǎn)。但是傳統(tǒng)的離散裂縫模型基于匹配性網(wǎng)格劃分,造成網(wǎng)格剖分復(fù)雜且計算誤差大。本文采用嵌入式離散裂縫模型,提高了計算效率。但是在大規(guī)模油藏數(shù)值模擬中,地質(zhì)模型可能包含數(shù)百萬甚至數(shù)億個網(wǎng)格,采用傳統(tǒng)的數(shù)值方法對其進(jìn)行求解,仍然計算量巨大,將超出當(dāng)今計算機(jī)的計算能力。本文提出嵌入式離散裂縫模型的多尺度模擬有限差分計算格式,在保證計算精度的同時大幅減小計算量。

(2) 在小尺度上使用模擬有限差分法求解嵌入式離散裂縫模型，準(zhǔn)確獲取了多尺度速度和壓力基函數(shù),并通過基函數(shù)映射得到小尺度解。得益于模擬有限差分的靈活性,可以處理任何復(fù)雜網(wǎng)格,為進(jìn)行復(fù)雜地層模擬奠定基礎(chǔ)。

(3) 多尺度模擬有限差分的基函數(shù)可以采用并行計算得到,進(jìn)一步減少了計算量。因此,本文方法對于裂縫性油藏數(shù)值模擬有很高的潛在價值。