牛 駿, 蘇建政, 嚴(yán) 俠, 汪友平, 孫 海

(1.國(guó)家能源頁(yè)巖油研發(fā)中心， 中國(guó)石化石油勘探開(kāi)發(fā)研究院，北京 100083；2.中國(guó)石油大學(xué)(華東)石油工程學(xué)院，青島 266580)

頁(yè)巖油氣作為一種不可或缺的非常規(guī)油氣資源，是目前研究和開(kāi)發(fā)的熱點(diǎn)[1]。頁(yè)巖油氣藏儲(chǔ)層致密，其滲透率非常低，一般必須經(jīng)過(guò)大規(guī)模水力壓裂才可以進(jìn)行商業(yè)化開(kāi)采[2]。壓裂后的頁(yè)巖儲(chǔ)層，由于水力裂縫的存在使其在滲流機(jī)理和流動(dòng)特征上均不同于常規(guī)儲(chǔ)層，表現(xiàn)出明顯的應(yīng)力敏感特征，且具有較強(qiáng)的非均質(zhì)各向異性，常規(guī)基于修正絕對(duì)滲透率的應(yīng)力敏感模型已不適用，亟需建立滲流場(chǎng)和應(yīng)力場(chǎng)耦合的流固耦合模型[3]。然而，相比于常規(guī)連續(xù)性儲(chǔ)層，頁(yè)巖油藏由于水力裂縫的存在，導(dǎo)致流動(dòng)和力學(xué)性質(zhì)發(fā)生較大變化，流固耦合效應(yīng)更明顯，處理難度更大，常用的裂縫性儲(chǔ)層流固耦合數(shù)學(xué)模型有：等效連續(xù)介質(zhì)模型[4-5]、雙重介質(zhì)模型[6-7]和離散介質(zhì)模型[8-9]。

等效連續(xù)介質(zhì)模型將裂縫性介質(zhì)看作一個(gè)連續(xù)體，首先需要判斷表征單元體的存在性，并確定其大小，再計(jì)算等效滲流和力學(xué)參數(shù)表征裂縫對(duì)流固兩場(chǎng)的影響，在此基礎(chǔ)上根據(jù)成熟的連續(xù)介質(zhì)力學(xué)理論進(jìn)行耦合流動(dòng)分析，模型數(shù)據(jù)需求簡(jiǎn)單而且計(jì)算效率高，一直以來(lái)在裂隙介質(zhì)的水力學(xué)研究中占據(jù)重要地位。然而，由于等效連續(xù)介質(zhì)模型過(guò)于簡(jiǎn)化、等效參數(shù)計(jì)算不夠準(zhǔn)確(尤其對(duì)于多相流問(wèn)題)，無(wú)法準(zhǔn)確刻畫裂縫性介質(zhì)中的多種流動(dòng)模式共存特征，因此，在油藏實(shí)際應(yīng)用中其分析結(jié)果往往誤差較大；雙重介質(zhì)模型最早由Barenblatt等[6]建立，Warren等[10]對(duì)其進(jìn)行了簡(jiǎn)化和推廣，此后眾多學(xué)者將多孔彈性力學(xué)與雙重介質(zhì)模型相結(jié)合用于裂縫性介質(zhì)流固耦合研究，這類模型將裂縫系統(tǒng)和基巖系統(tǒng)之間的物質(zhì)交換考慮在內(nèi)，可以實(shí)際地模擬裂縫內(nèi)的優(yōu)先流現(xiàn)象。然而被裂縫分割成的基質(zhì)巖塊系統(tǒng)被假定具有規(guī)則的大小和形狀，過(guò)于簡(jiǎn)化卻不能充分表現(xiàn)出裂縫性介質(zhì)的不連續(xù)性、各向異性等特點(diǎn)，同時(shí)基巖和裂縫間的竄流系數(shù)也難以準(zhǔn)確獲取(尤其對(duì)于多相、非穩(wěn)態(tài)問(wèn)題)，此外該模型僅適用于高度發(fā)育、空間尺度跨越小的裂縫性介質(zhì)；離散介質(zhì)模型對(duì)裂縫顯示表征，能夠反應(yīng)裂縫對(duì)滲流場(chǎng)與應(yīng)力場(chǎng)的影響。由于裂縫開(kāi)度和油藏尺度相差極大，需要對(duì)裂縫進(jìn)行降維處理，以防裂縫周圍局部加密而增加網(wǎng)格量。然而裂縫兩側(cè)的位移場(chǎng)具有間斷性，裂縫降維后采用常規(guī)數(shù)值方法難以有效解決，對(duì)此學(xué)者提出了相應(yīng)的解決方法，主要分為以下三類：界面雙節(jié)點(diǎn)單元法[11]、位移不連續(xù)法[12]和擴(kuò)展有限元法[13]。界面雙節(jié)點(diǎn)單元法在常規(guī)有限單元法的基礎(chǔ)上，在裂縫單元節(jié)點(diǎn)處新增一個(gè)額外節(jié)點(diǎn)，反映裂縫兩側(cè)變量的不連續(xù)性。該方法求解較簡(jiǎn)單，但其網(wǎng)格劃分和編號(hào)存在一定難度，并且數(shù)值求解穩(wěn)定性和收斂性較差。位移不連續(xù)法作為一種邊界元方法，求解簡(jiǎn)單、網(wǎng)格量小、精度較高，但其僅適用于均質(zhì)儲(chǔ)層，對(duì)于實(shí)際強(qiáng)非均質(zhì)各向異性儲(chǔ)層難以推導(dǎo)解析解建立離散格式，并且該方法最終形成的系數(shù)矩陣為滿秩矩陣，對(duì)于大規(guī)模問(wèn)題的求解難度大。擴(kuò)展有限元法在標(biāo)準(zhǔn)有限元框架下，在連續(xù)區(qū)域內(nèi)采取標(biāo)準(zhǔn)有限元方法，在包含不連續(xù)邊界的區(qū)域內(nèi)，增加間斷形函數(shù)以修正有限元的位移近似函數(shù)，采取水平集方法來(lái)描述間斷界面，使間斷的描述獨(dú)立于有限元網(wǎng)格，避免了以裂縫為約束的復(fù)雜非結(jié)構(gòu)網(wǎng)格劃分，但擴(kuò)展有限元作為一種改進(jìn)的有限元方法，其方程離散為全局離散，不滿足單元的局部物質(zhì)守恒，因此對(duì)于滲流方程(尤其多相流)，其計(jì)算精度較差。

近些年來(lái)，Yan等[14]基于等效連續(xù)介質(zhì)模型和嵌套離散裂縫模型，提出了一種有效的數(shù)學(xué)模型來(lái)模擬含多尺度裂縫和溶洞的裂隙型多孔介質(zhì)中流體力學(xué)耦合問(wèn)題；張世明等[15]基于離散介質(zhì)模型，提出了一種求解裂縫性介質(zhì)等效滲透率的新方法；李友全等[16]基于離散介質(zhì)模型，對(duì)致密油藏壓裂直井進(jìn)行了動(dòng)態(tài)分析；孫若凡等[17]基于雙重介質(zhì)模型，對(duì)致密油藏壓裂水平井進(jìn)行了模擬分析；盛茂等[18]基于雙重介質(zhì)模型，建立了一種裂縫性頁(yè)巖氣藏的流固耦合模型并采用有限元法進(jìn)行數(shù)值求解。然而目前并沒(méi)有基于離散介質(zhì)模型并使用擴(kuò)展有限元法求解的研究。于是，將擴(kuò)展有限元法與嵌入式離散裂縫模型[19]相結(jié)合，建立一種適用于裂縫性頁(yè)巖油藏流固耦合數(shù)值模擬的新方法，該方法可以考慮裂縫開(kāi)度的變化以及裂縫內(nèi)流壓對(duì)開(kāi)度變化的影響，分別采用擴(kuò)展有限元和模擬有限差分[19-20]求解力學(xué)和滲流方程，模擬有限差方法適用于強(qiáng)非均質(zhì)各向異性儲(chǔ)層，且具有良好的局部守恒性，可以準(zhǔn)確求解壓力場(chǎng)和速度場(chǎng)。嵌入式離散裂縫模型網(wǎng)格劃分與擴(kuò)展有限元類似，無(wú)需考慮儲(chǔ)層內(nèi)的裂縫形態(tài)，大大降低了網(wǎng)格劃分難度，其計(jì)算精度與基于匹配性非結(jié)構(gòu)網(wǎng)格的常規(guī)離散裂縫模型基本一致。

1 流固耦合數(shù)學(xué)模型

基于二維單相流固耦合問(wèn)題來(lái)闡述本文數(shù)值模擬方法的原理與思想。假設(shè)在耦合流動(dòng)過(guò)程中溫度恒定，考慮基巖和流體的壓縮性，基巖與裂縫中的流體流動(dòng)均滿足達(dá)西定律，巖石變形為微小線彈性變形，不考慮慣性力。采用常規(guī)多孔介質(zhì)彈性力學(xué)符號(hào)約定，即拉伸為正，壓縮為負(fù)。

1.1 巖石變形數(shù)學(xué)模型

(1)

σ=σeff-αpmI

(2)

(3)

(4)

σ·nf=-pf·nfon Γf

(5)

圖1 裂縫性介質(zhì)及其內(nèi)、外邊界示意圖Fig.1 The sketch of a fractured porous medium and its boundaries

1.2 嵌入式離散裂縫模型

基巖系統(tǒng)流動(dòng)數(shù)學(xué)模型：

(6)

(7)

(8)

裂縫系統(tǒng)流動(dòng)數(shù)學(xué)模型：

(9)

(10)

式中:vm和vf分別為基巖和裂縫內(nèi)滲流速度，m/s；km和Kf分別為基巖和裂縫滲透率，m2；μ為流體黏度，mPa·s；ρl為流體密度，kg/m3；φ為基巖孔隙度；Ks和Kl分別為巖石顆粒和流體的體積模量，Pa；εv為基巖體積應(yīng)變；qm和qf分別為基巖和裂縫的源匯項(xiàng)，s-1；Vm和Vf分別為基巖和裂縫的單元體積，m3；bf0和bf分別為裂縫初始和目前的開(kāi)度，m；qmf為基巖與裂縫間的竄流量，m3/s；qff為相交裂縫間的竄流量，m3/s；δ為狄克拉函數(shù)；qmf、qff、δmf和δff的詳細(xì)計(jì)算格式和取值方法見(jiàn)參考文獻(xiàn)[19]；Γq為流量外邊界，Γp為定壓外邊界，如圖1所示。整個(gè)系統(tǒng)的初始?jí)毫閜0。

2 模型離散求解

圖2 網(wǎng)格單元數(shù)值離散及自變量分布示意圖Fig.2 The schematic of discretization and locations of the solution variables

2.1 擴(kuò)展有限元離散力學(xué)方程

擴(kuò)展有限元根據(jù)單位分解法[21](PUM)的思想，在常規(guī)有限元連續(xù)位移場(chǎng)的基礎(chǔ)上，增加不連續(xù)位移場(chǎng)(裂縫段采用Heaviside函數(shù)，裂縫尖端采用尖端漸進(jìn)函數(shù))來(lái)表征計(jì)算區(qū)域內(nèi)的位移間斷特性，因此，其單元位移的近似解可以表示為

H[φ(xi)]}ai+

(11)

(12)

H(ξ)=1, ?ξ>0 andH(ξ)=-1, ?ξ<0

(13)

F(r,θ)=

(14)

將式(2)代入式(1)，考慮邊界條件[式(4)和式(5)]，得到方程(1)的積分弱形式為

(15)

式(15)中：δu為位移試函數(shù)；δu+和δu-分別表示裂縫正負(fù)兩側(cè)的位移試函數(shù)，如圖1所示。

裂縫壁面發(fā)生相對(duì)位移后，裂縫開(kāi)度及滲透率(立方定律)表達(dá)式為

df=df0+(u+-u-)·nf

(16)

(17)

在擴(kuò)展有限元中，采用水平集方法[22]確定裂縫的幾何信息，同時(shí)采用高斯數(shù)值積分方法計(jì)算單元和邊界上的積分，而對(duì)于裂縫穿過(guò)的單元，常規(guī)的高斯積分方法不適用，需要根據(jù)裂縫形狀對(duì)這些單元進(jìn)行單元分解，再在每個(gè)分解子單元內(nèi)進(jìn)行數(shù)值積分。

2.2 模擬有限差分離散滲流方程

(18)

對(duì)于方程(7)，在基質(zhì)網(wǎng)格單元上積分且應(yīng)用散度定理得

(19)

考慮裂縫滲透率為標(biāo)量且為一維流動(dòng)，因此采用有限差分進(jìn)行離散。將方程(9)代入方程(10)，方程左右兩端同乘裂縫網(wǎng)格單元體積得

(20)

2.3 整體計(jì)算格式

以存在一條裂縫為例，考慮基巖網(wǎng)格單元邊界面上的速度連續(xù)，將力學(xué)方程和滲流方程的數(shù)值離散格式組裝到一起，得到嵌入式離散裂縫流固耦合模型數(shù)值計(jì)算格式為

(21)

(22)

(23)

式中:Ku為擴(kuò)展有限元位移剛度矩陣；Qum和Quf分別為基巖和裂縫的流固耦合系數(shù)矩陣；Km為基巖模擬有限差分系數(shù)矩陣；Tmf為基巖與裂縫竄流系數(shù)矩陣；Kf為嵌入式裂縫有限差分系數(shù)矩陣；Sm和Sf分別為基巖和裂縫的壓縮項(xiàng)系數(shù)矩陣；u、a和b分別為節(jié)點(diǎn)位移和額外未知量，整體對(duì)應(yīng)于系數(shù)矩陣A的第一行；vm、pm和πm分別為網(wǎng)格邊界流速、網(wǎng)格中心壓力和網(wǎng)格邊界壓力，整體對(duì)應(yīng)于系數(shù)矩陣A的第二行；pf為裂縫單元壓力；C為體積力載荷、邊界應(yīng)力載荷以及基巖和裂縫源匯項(xiàng)列陣。

對(duì)于方程(21)中的時(shí)間導(dǎo)數(shù)項(xiàng)，采用全隱式差分格式進(jìn)行離散。

(Adt+B)Xn+1=BXn+Cdt

(24)

式(24)中：n表示第n個(gè)時(shí)間步；dt為時(shí)間步長(zhǎng)。

3 算例分析

首先分別通過(guò)與一維多孔介質(zhì)和二維裂縫性介質(zhì)流固耦合問(wèn)題的解析解和精細(xì)網(wǎng)格數(shù)值解進(jìn)行對(duì)比，驗(yàn)證所提方法的正確性；然后分析不同條件下裂縫開(kāi)度的變化情況以及裂縫開(kāi)度變化對(duì)耦合流動(dòng)的影響；最后對(duì)于某一二維頁(yè)巖油藏彈性開(kāi)采進(jìn)行數(shù)值模擬，對(duì)比基巖滲透率為對(duì)角張量和全張量時(shí)的差異，說(shuō)明所提方法對(duì)于全張量滲透率的適用性。本文算例均采用平面應(yīng)變假設(shè)。

3.1 模型驗(yàn)證

考慮如圖3(a)所示的深度200 m的一維流固耦合模型，初始流壓和位移均為0，在上、下邊界同時(shí)施加100 kN/m2的恒定邊界載荷，左、右邊界封閉，上、下為邊界定壓(壓力為0)，模型參數(shù)如表1所示。采用所提方法對(duì)該一維模型進(jìn)行數(shù)值求解并與其解析解進(jìn)行對(duì)比，圖3(b)給出了不同時(shí)刻流體壓力的數(shù)值解與解析解對(duì)比結(jié)果(由于模型上、下對(duì)稱，圖中只給出了上半部分的壓力對(duì)比)，能夠看出數(shù)值解與解析解基本一致。

表1 模型基本參數(shù)Table 1 The model basic parameters

圖3 一維模型幾何示意圖及不同時(shí)間流體壓力數(shù)值解與解析解對(duì)比Fig.3 The sketch of a 1D model and the comparison of numerical solution and analytical solution of fluid pressure at different times

如圖4(a)所示，為一10 m×20 m 的裂縫性介質(zhì)模型，其中裂縫長(zhǎng)10 m、開(kāi)度0. 001 m，垂直位于模型區(qū)域中心。初始流壓和位移均為0，模型的力學(xué)邊界條件如圖4(a)所示，滲流邊界條件為：上邊界定壓(壓力為0)，其余邊界封閉。裂縫孔隙為1.0，初始滲透率為8.33×10-8m2，其他參數(shù)如表1所示。由于裂縫的存在，該模型難以推導(dǎo)解析解，因此，將裂縫顯示處理，通過(guò)在裂縫周圍進(jìn)行局部網(wǎng)格加密構(gòu)建精細(xì)網(wǎng)格，如圖4(b)所示，采用有限單元法求解作為參考解。圖4(c)為所提方法采用的計(jì)算網(wǎng)格，可以看出其網(wǎng)格劃分難度及網(wǎng)格數(shù)量都大大降低。

不考慮裂縫開(kāi)度變化，分別采用精細(xì)網(wǎng)格有限元方法和所提方法對(duì)該裂縫性介質(zhì)流固耦合模型進(jìn)行模擬，圖4(d)和圖4(e)分別為模擬10 000 s后兩種方法得到的壓力場(chǎng)圖，圖5(a)和圖5(b)分別給出了兩種方法計(jì)算得到的測(cè)試點(diǎn)A和B處壓力和y方向位移隨時(shí)間的變化曲線對(duì)比。通過(guò)對(duì)比可以看出，所提方法的計(jì)算結(jié)果與精細(xì)網(wǎng)格有限元方法參考解基本一致。因此，本節(jié)的兩個(gè)算例驗(yàn)證了所提方法的正確性。

圖4 裂縫性介質(zhì)模型示意圖、不同方法網(wǎng)格劃分結(jié)果及模擬10 000 s后壓力場(chǎng)圖Fig.4 The model and grids of a fractured medium, pressure distribution after 10 000 seconds for different methods

圖5 兩種方法計(jì)算得到的測(cè)試點(diǎn)壓力和位移隨時(shí)間變化曲線Fig.5 Pressure and displacement (y) the test point various time obtained from two methods

3.2 裂縫開(kāi)度隨時(shí)間變化

由于本文模型分別采用嵌入式離散裂縫和擴(kuò)展有限元中將裂縫的滲流和力學(xué)特征顯示表征，因此可以有效刻畫裂縫開(kāi)度隨時(shí)間的變化情況以及裂縫開(kāi)度變化對(duì)耦合流動(dòng)的影響?；趫D4(a)所示的裂縫多孔介質(zhì)模型，模型參數(shù)與上一個(gè)算例一致，不考慮裂縫內(nèi)流體壓力的影響，分別采用所提方法和精細(xì)網(wǎng)格有限元方法計(jì)算100、2 000、10 000 s時(shí)的裂縫開(kāi)度分布，結(jié)果如圖6(a)所示，能夠看出兩者結(jié)果基本一致，驗(yàn)證了所提方法用于計(jì)算裂縫開(kāi)度的正確性。圖6(b)給出了不考慮與考慮裂縫內(nèi)流體壓力影響兩種情況下所提方法計(jì)算得到的不同時(shí)間裂縫開(kāi)度分布，可以看出模擬100 s時(shí)，兩者裂縫開(kāi)度相差很大，隨著模擬時(shí)間增長(zhǎng)，兩者之間的差距逐漸減小，到模擬10 000 s時(shí)，兩者裂縫開(kāi)度基本一致。這是由于早期裂縫內(nèi)流體壓力較高，考慮其影響后對(duì)裂縫起到支撐作用，可以有效阻止裂縫閉合，隨著模擬時(shí)間增長(zhǎng)，裂縫內(nèi)的流體壓力逐漸減小，其支撐作用逐漸減弱。

圖6 不同時(shí)間裂縫開(kāi)度分布Fig.6 Fracture aperture distribution at different time

圖7 模擬10 000 s后的壓力場(chǎng)分布圖Fig.7 Pressure distribution after 10 000 seconds

在以上模型的基礎(chǔ)上，將上邊界載荷增加到200 kN/m2，其他模型參數(shù)和條件均不變。采用所提方法分別對(duì)裂縫開(kāi)度不變和變化兩種情況下的流固耦合進(jìn)行模擬，圖7(a)和7(b)分別給出了模擬10 000 s后兩種情況下壓力場(chǎng)分布，圖8給出了兩種情況下A點(diǎn)處壓力及兩者誤差隨時(shí)間的變化曲線，可以看出：考慮裂縫開(kāi)度變化后，隨著模擬時(shí)間增加，裂縫開(kāi)度和導(dǎo)流能力逐漸減小，流體流出速度減小，區(qū)域內(nèi)的流體壓力增加，兩者相對(duì)誤差逐漸增大。

圖8 兩種情況下A點(diǎn)處壓力及兩者相對(duì)誤差隨時(shí)間變化曲線Fig.8 Pressure and relative error at point A for various time in the two cases

3.3 裂縫性頁(yè)巖油藏彈性開(kāi)采

圖9 裂縫性頁(yè)巖油藏模型示意圖及網(wǎng)格劃分結(jié)果Fig.9 The sketch of fractured shale reservoir model and meshing results

表2 頁(yè)巖油藏基本參數(shù)Table 2 The reservoir parameters

圖10～圖12分別給出了基巖滲透率為對(duì)角張量和全張量?jī)煞N情況下，模擬10 d后的壓力場(chǎng)和x方向位移場(chǎng)圖、不同裂縫的開(kāi)度變化情況(裂縫1為最左側(cè)裂縫，裂縫3為從左側(cè)數(shù)第三條裂縫)以及累產(chǎn)隨時(shí)間變化情況?？梢钥闯?，考慮全張量滲透率后，壓力場(chǎng)、位移場(chǎng)、裂縫開(kāi)度和累積產(chǎn)量相比于對(duì)角張量滲透率條件下的結(jié)果有一定變化，說(shuō)明對(duì)于基巖滲透率為全張量的實(shí)際強(qiáng)各向異性儲(chǔ)層，不能采用有限差分或有限體積等忽略非主對(duì)角線值的常規(guī)數(shù)值方法。

圖10 對(duì)角張量滲透率下彈性開(kāi)采10 d時(shí)壓力場(chǎng)和x方向位移Fig.10 Pressure and displacement distribution after 10 days for diagonal-tensorpermeability

圖11 全張量滲透率下彈性開(kāi)采10 d時(shí)壓力場(chǎng)和x方向位移Fig.11 Pressure and displacement (x) distribution after 10 days for full-tensorpermeability

圖12 不同基巖滲透率條件下，模擬10 d裂縫開(kāi)度和累產(chǎn)對(duì)比Fig.12 Fracture apertures and cumulative production for different matrix permeability on different bedrock permeability conditions after 10 days

4 結(jié)論

通過(guò)耦合嵌入式離散裂縫模型和擴(kuò)展有限元提出一種適用于裂縫性頁(yè)巖油藏的流固耦合高效數(shù)值模擬方法，該方法的創(chuàng)新性以及所取得的結(jié)論如下。

(1)該方法對(duì)裂縫使用降維顯式處理，能夠準(zhǔn)確刻畫裂縫對(duì)滲流場(chǎng)和應(yīng)力場(chǎng)的互相影響，同時(shí)可以模擬裂縫開(kāi)度的動(dòng)態(tài)變化過(guò)程以及裂縫內(nèi)流體壓力對(duì)裂縫變形的影響。該模型劃分網(wǎng)格時(shí)不必考慮裂縫幾何形態(tài)，只需對(duì)基巖系統(tǒng)進(jìn)行正交網(wǎng)格剖分，使得網(wǎng)格劃分時(shí)的復(fù)雜度大幅降低。

(2)對(duì)于模型中流動(dòng)方程求解采取模擬有限差分方法，滿足局部質(zhì)量守恒并且可以有效處理全張量形式滲透率；力學(xué)方程求解采取擴(kuò)展有限元方法，可以準(zhǔn)確處理裂縫兩側(cè)的位移間斷性。

(3)分別通過(guò)與一維多孔介質(zhì)和二維裂縫性介質(zhì)流固耦合問(wèn)題的解析解和精細(xì)網(wǎng)格數(shù)值解進(jìn)行對(duì)比，驗(yàn)證所提方法的正確性，分析了不同條件下裂縫開(kāi)度的變化情況以及裂縫開(kāi)度變化對(duì)耦合流動(dòng)的影響。

(4)針對(duì)一、二維裂縫性頁(yè)巖油藏彈性開(kāi)采進(jìn)行了數(shù)值模擬，對(duì)比基巖滲透率為對(duì)角張量和全張量時(shí)的差異，說(shuō)明所提方法對(duì)于全張量滲透率的適用性。