壓裂油藏導(dǎo)流裂縫多相流動的高效有限差分模型

周方奇，施安峰，王曉宏

（中國科學技術(shù)大學熱科學和能源工程系）

0 引言

壓裂技術(shù)作為開發(fā)低滲透油藏的有效增產(chǎn)手段已被廣泛應(yīng)用[1-4]。壓裂中產(chǎn)生若干條尺度較大的導(dǎo)流裂縫，為準確預(yù)測壓裂后油藏中的流體流動狀況，需采用準確的導(dǎo)流裂縫數(shù)值模擬方法。

對于裂縫型油藏，雙重介質(zhì)模型[5-8]已被廣泛應(yīng)用，而含有大尺度導(dǎo)流裂縫的油藏并不滿足該模型的使用條件。早期學者們使用單重介質(zhì)模型模擬導(dǎo)流裂縫，但由于該模型需要極細密的網(wǎng)格描述裂縫和基質(zhì)，所以計算量很大、效率很低，難以應(yīng)用到實際的油藏數(shù)值模擬中。為了提高計算效率，一些學者提出離散裂縫模型[9-15]，假設(shè)N維的系統(tǒng)中裂縫內(nèi)流體的流動是N?1維的，在劃分網(wǎng)格時裂縫節(jié)點與基質(zhì)節(jié)點相互重合，并把裂縫及基質(zhì)方程疊加后進行計算，該模型可較準確地模擬單相流動，但無法應(yīng)用于多相流動問題。隨后一些學者對離散裂縫模型進行了改進[16-19]，使用基于達西定律的交換流量耦合裂縫及基質(zhì)方程，改進后的模型可以模擬裂縫油藏中的多相流動，但是仍然需要沿著裂縫的走勢劃分基質(zhì)網(wǎng)格。Lee等[20-22]提出了嵌入裂縫模型，把裂縫近似地看作基質(zhì)網(wǎng)格中的井源，采用類Peaceman公式[23-24]計算裂縫與基質(zhì)網(wǎng)格之間的物質(zhì)交換，雖然該模型可以提高網(wǎng)格劃分的效率，但在多相流動問題中，并不能準確描述裂縫和基質(zhì)之間的物質(zhì)交換。

本文考慮基質(zhì)相飽和度在裂縫兩側(cè)的不連續(xù)性，提出一種新的嵌入裂縫模型，以期提高計算效率和精度，準確模擬壓裂后油藏中的流動狀況。

1 數(shù)學模型

在人工壓裂低滲透油藏中，油藏介質(zhì)為流體的主要存儲空間，其滲透率較低，通常被稱為基質(zhì)。裂縫滲透能力強，成為流體的主要流通通道。本文考慮油水兩相流問題，以二維空間為例進行研究。由于裂縫的開度很小，其中流體的流動可以認為是一維的，設(shè)裂縫的開度為wf，單位時間內(nèi)從單位長度裂縫流入基質(zhì)的α（α=w, o）相流體流量為qα，則裂縫中α相流體的質(zhì)量守恒方程為：

把裂縫看作基質(zhì)中的源項，基質(zhì)方程可寫作：

其中α相流體的流速滿足達西定律：

由于裂縫開度很小，基質(zhì)中的壓力場在跨越裂縫時保持連續(xù)，但裂縫和基質(zhì)中多相流體流動特征不同，所以在兩者界面處相飽和度發(fā)生間斷，基質(zhì)的相飽和度場在越過裂縫時并不連續(xù)（見圖 1），因此裂縫和其兩側(cè)基質(zhì)的物質(zhì)交換并不相等，基質(zhì)與裂縫的交換流量為：

圖1 油藏基質(zhì)對角線（裂縫法向）上的油相飽和度分布曲線

裂縫流入兩側(cè)基質(zhì)的流量qα,1和qα,2滿足：

2 控制方程的數(shù)值求解

2.1 裂縫及基質(zhì)的網(wǎng)格劃分

對于含有復(fù)雜裂縫的油藏，使用不規(guī)則網(wǎng)格劃分基質(zhì)系統(tǒng)可以有效地貼合裂縫，但當裂縫分布較為密集時，這類網(wǎng)格的劃分過程十分復(fù)雜，且工作量巨大。本文使用正交結(jié)構(gòu)網(wǎng)格劃分油藏基質(zhì)，網(wǎng)格的空間步長為 Δx=Δy=1.0 m，將裂縫看作嵌入在基質(zhì)網(wǎng)格中，這樣可以對含有復(fù)雜裂縫的油藏進行高效網(wǎng)格劃分，以提高模擬效率。

2.2 裂縫網(wǎng)格的計算

嵌入基質(zhì)的裂縫被規(guī)則的基質(zhì)網(wǎng)格分割為Nf段一維裂縫網(wǎng)格，假設(shè)其開度都為wf，取第k段裂縫網(wǎng)格為研究對象，其長度為Lk。由于裂縫內(nèi)流體的流動被認為是一維的，所以可以使用常規(guī)的一維離散格式將（1）式全隱式離散為：

通過裂縫網(wǎng)格2個界面的流速可由下式得到：

對于裂縫與基質(zhì)之間的流量交換，無法使用常規(guī)格式進行離散。本文在裂縫兩側(cè)分別建立參考點A和B（見圖 2），其位于裂縫網(wǎng)格的垂直平分線上，與裂縫的距離分別為dA和dB。由于壓力場在裂縫附近連續(xù)分布，且其導(dǎo)數(shù)在單側(cè)基質(zhì)中也保持連續(xù)，所以參考點A和B上的壓力可以通過裂縫及基質(zhì)網(wǎng)格中的壓力單側(cè)插值予以確定，因此裂縫與兩側(cè)基質(zhì)的流量交換可由下式求出：

圖2 裂縫網(wǎng)格兩側(cè)參考點示意圖

2.3 基質(zhì)網(wǎng)格的處理

在油藏中，把空間坐標為（i，j）的網(wǎng)格作為研究對象，使用二維九點離散格式，（2）式可以被全隱式離散為：

其中λα=TKrα,mρα/μα，系數(shù)T滿足：

利用（12）式可以對常規(guī)基質(zhì)網(wǎng)格進行離散，由于裂縫穿過基質(zhì)網(wǎng)格，受其影響流體的某些物理量會在網(wǎng)格中發(fā)生間斷，導(dǎo)致該網(wǎng)格的參數(shù)無法準確定義，因此本文在對基質(zhì)方程進行離散時，這種網(wǎng)格將不參與實際計算，而是為其他網(wǎng)格的計算提供邊界條件。由于裂縫兩側(cè)基質(zhì)中壓力梯度分別保持連續(xù)，可以通過單側(cè)插值的方法由相鄰基質(zhì)網(wǎng)格及裂縫網(wǎng)格中的壓力確定被裂縫隔開的兩部分基質(zhì)中的壓力分布，然后得到由網(wǎng)格（i，j）進入網(wǎng)格（i－1，j）和（i，j+1）的流量qi－1/2,j和qi,j+1/2（見圖 3），從而為網(wǎng)格（i－1，j）和（i，j+1）的計算提供邊界條件。同理，該類網(wǎng)格與相鄰常規(guī)網(wǎng)格的流量交換都可以使用以上的方法求出，因此，圖 3中虛線界面可以被認為是整個油藏基質(zhì)中的內(nèi)邊界，且滿足第二類邊界條件。

圖3 內(nèi)邊界替代覆蓋裂縫的基質(zhì)網(wǎng)格

在求解裂縫控制方程時，若裂縫位于油藏內(nèi)部，則在兩端取絕流邊界條件；若裂縫延伸至油藏邊界，其兩端的邊界條件則由油藏的外邊界條件確定。

3 算例

本文用 2個算例驗證模型的可靠性，并將計算結(jié)果與改進后的離散裂縫模型（參考模型）及 Lee的模型進行比較。

首先考慮包含直裂縫（即裂縫與油藏邊界正交）的油藏，基質(zhì)和裂縫的滲透率分別為 1.0 μm2和10 000.0 μm2[22]，注入井和生產(chǎn)井井底壓力分別為10 MPa和 7 MPa，油相和水相的初始飽和度分別為 0.7和0.3。開井2 d后，油藏中基質(zhì)和裂縫的油相飽和度分布情況見圖4?？梢钥闯霰疚哪Ｐ团c參考模型的計算結(jié)果基本相同（見圖5），這是因為本文模型考慮了裂縫兩側(cè)相飽和度的間斷，分別計算裂縫與兩側(cè)基質(zhì)的流量交換，所以可以準確地描述裂縫對油藏的影響，得到較精確的結(jié)果；而 Lee模型并未考慮裂縫兩側(cè)基質(zhì)流場，所以對交換流量的描述存在誤差，從而影響了模擬結(jié)果的準確性。

圖4 包含直裂縫的油藏開井2 d后油相飽和度分布圖

對于含有斜裂縫的油藏，即裂縫與邊界并不正交的油藏，假設(shè)其中基質(zhì)和裂縫的滲透率分別為1.0×10?3μm2和20.0 μm2，注入井和生產(chǎn)井井底壓力分別為100 MPa和70 MPa[22]，通過數(shù)值模擬計算出開井100 d后基質(zhì)和裂縫中油相飽和度分布情況（見圖 6），與直裂縫油藏中的情況類似，無論是在基質(zhì)中還是在裂縫中，本文模型的結(jié)果都與參考結(jié)果保持一致，且優(yōu)于 Lee模型（見圖7）。

圖5 直裂縫中的油相飽和度分布曲線

圖6 包含斜裂縫的油藏開井100 d后油相飽和度分布圖

圖7 斜裂縫中的油相飽和度分布曲線

以上計算實例說明，對于油藏中不同形態(tài)的導(dǎo)流裂縫，本文模型都可以準確地進行模擬，而且計算精度與改進后的離散裂縫模型相同，高于Lee模型。

4 結(jié)論

本文采用規(guī)則網(wǎng)格劃分油藏區(qū)域，裂縫被覆蓋在基質(zhì)網(wǎng)格中，大幅度降低了網(wǎng)格劃分難度和計算量；同時考慮到流體相飽和度在裂縫兩側(cè)的不連續(xù)性，分別計算裂縫和其兩側(cè)基質(zhì)之間的物質(zhì)交換量；在離散基質(zhì)方程時，使用內(nèi)邊界替代含有裂縫的基質(zhì)網(wǎng)格參與計算，有效避免了因這種網(wǎng)格內(nèi)物理量被模糊定義而產(chǎn)生的誤差；算例表明，本文模型可以對含有導(dǎo)流裂縫的油藏進行高效準確的流動模擬。

本文只考慮了獨立的裂縫，應(yīng)對裂縫與裂縫之間的相互交叉與影響作進一步研究。

符號注釋：

φ——孔隙度，f；ρ——流體密度，kg/m3；S——流體飽和度，f；um——基質(zhì)中流體的流速，m/s；uf——裂縫中流體的流速，m/s；q——流量，m3/s；t——時間，s；r——原點到方程控制點的向量，m；τf——裂縫的切向向量，m；τf——裂縫切向向量的模，m；μ——流體的黏度，Pa·s；K——絕對滲透率，μm2；Kr——相對滲透率，f；p——流體的壓力，Pa；n——裂縫表面法向向量；Δt——時間步長，s；d——參考點與裂縫的距離，m。下標：α——相，α=w，o；β——介質(zhì)，β= m，f；m——基質(zhì)；f——裂縫；0——上時刻的變量；1，2——裂縫的2個表面。

[1]李根生, 盛茂, 田守 嶒, 等.水平井水力噴射分段酸壓技術(shù)[J].石油勘探與開發(fā), 2012, 39(1): 100-104.Li Gensheng, Sheng Mao, Tian Shouceng, et al.Multistage hydraulic jet acid fracturing technique for horizontal wells[J].Petroleum Exploration and Development, 2012, 39(1): 100-104.

[2]賈承造, 鄭民, 張永峰.中國非常規(guī)油氣資源與勘探開發(fā)前景[J].石油勘探與開發(fā), 2012, 39(2): 129-136.Jia Chengzao, Zheng Min, Zhang Yongfeng.Unconventional hydrocarbon resources in China and the prospect of exploration and development[J].Petroleum Exploration and Development, 2012,39(2): 129-136.

[3]黃昌武.頁巖油氣壓裂理論方法和技術(shù)體系初步形成[J].石油勘探與開發(fā), 2012, 39(4): 443.Huang Changwu.Theoretical methodology and technical system for shale reservoir fracturing[J].Petroleum Exploration and Development,2012, 39(4): 443.

[4]徐鵬, 程遠方, 劉新云, 等.低滲透油氣藏爆炸壓裂模擬試驗及裂縫分形特征[J].石油勘探與開發(fā), 2013, 40(5): 636-640.Xu Peng, Cheng Yuanfang, Liu Xinyun, et al.Explosive fracturing simulation experiment for low permeability reservoirs and fractal characteristics of cracks produced by explosive fracturing[J].Petroleum Exploration and Development, 2013, 40(5): 636-640.

[5]Barenblatt G I, Zheltov I P.Basic flow equations for homogeneous fluids in naturally fractured rocks[J].Doklady Akademii Nauk, 1960,13: 545-548.

[6]Warren J E, Root P J.The behavior of naturally fractured reservoirs[J].Society of Petroleum Engineers Journal, 1963, 3(3):245-255.

[7]Lim K T, Aziz K.Matrix-fracture transfer shape factors for dual-porosity simulators[J].Journal of Petroleum Science and Engineering, 1995, 13(3/4): 169-178.

[8]Kazemi H, Merrill L S, Porterfield K L, et al.Numerical-simulation of water-oil flow in naturally fractured reservoirs[J].Society of Petroleum Engineers Journal, 1976, 16(6): 317-326.

[9]Noorishad J, Mehran M.An upstream finite element method for solution of transient transport equation in fractured porous media[J].Water Resources Research, 1982, 18(3): 588-596.

[10]Baca R G, Arnett R C, Langford D W.Modelling fluid flow in fractured-porous rock masses by finite-element techniques[J].International Journal for Numerical Methods in Fluids, 1984, 4(4):337-348.

[11]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.

[12]Karimi-Fard M, Firoozabadi A.Numerical simulation of water injection in fractured media using discrete-fracture model and the Galerkin method[J].SPE Reservoir Evaluation & Engineering, 2003,6(2): 117-126.

[13]Kim J G, Deo M D.Comparison of the performance of a discrete fracture multiphase model with those using conventional methods[R].SPE 51928, 1999.

[14]Monteagudo J, 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):1-20.

[15]Monteagudo J, Firoozabadi A.Control-volume model for simulation of water injection in fractured media: Incorporating matrix heterogeneity and reservoir wettability effects[J].Society of Petroleum Engineers Journal, 2007, 12(3): 355-366.

[16]Slough K J, Sudicky E A, Forsyth P A.Numerical simulation of multiphase flow and phase partitioning in discretely fractured geologic media[J].Journal of Contaminant Hydrology, 1999, 40(2):107-136.

[17]Slough K J, Sudicky E A, Forsyth P.A.Grid refinement for modeling multiphase flow in discretely fractured porous media[J].Advances in Water Resources, 1999, 23(3): 261-269.

[18]Karimi-Fard M, Durlofsky L J, Aziz K.An efficient discrete-fracture model applicable for general-purpose reservoir simulators[J].Society of Petroleum Engineers Journal, 2004, 9(2): 227-236.

[19]Sandve T H, 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.

[20]Lee S H, Jensen C L, Lough M F.An efficient finite difference model for flow in a reservoir with multiple length-scale fractures[J].Society of Petroleum Engineers Journal, 1999, 5(3): 268-275.

[21]Li L, Lee S H.Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media[J].SPE Reservoir Evaluation & Engineering,2008, 11(4): 750-758.

[22]Moinfar A, Narr W, Hui M H, et al.Comparison of discrete-fracture and dual-permeability models for multiphase flow in naturally fractured reservoirs[R].SPE 142295, 2011.

[23]Peaceman D W.Interpretation of well-block pressures in numerical reservoir simulation[J].Society of Petroleum Engineers Journal,1978, 19(3): 183-194.

[24]Peaceman D W.Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability[J].Society of Petroleum Engineers Journal, 1983, 23(3):531-543.