離散縫洞網(wǎng)絡(luò)油水兩相流動(dòng)的數(shù)值模擬

宋世瀚，邸 元，劉中春，代亞非

(1.北京大學(xué) 工學(xué)院，北京 100871；2.中國石油化工股份有限公司 石油勘探開發(fā)研究院，北京 100083； 3. 遠(yuǎn)景能源(江蘇)有限公司，上海 200051)

?

離散縫洞網(wǎng)絡(luò)油水兩相流動(dòng)的數(shù)值模擬

宋世瀚1，邸 元1，劉中春2，代亞非3

(1.北京大學(xué) 工學(xué)院，北京 100871；2.中國石油化工股份有限公司 石油勘探開發(fā)研究院，北京 100083； 3. 遠(yuǎn)景能源(江蘇)有限公司，上海 200051)

基于溶洞內(nèi)油水兩相流體瞬時(shí)平衡和重力分異的假定，采用流相流動(dòng)分?jǐn)?shù)模擬離散縫洞網(wǎng)絡(luò)中油水兩相的流動(dòng)，提出一種離散縫洞網(wǎng)絡(luò)兩相流動(dòng)的數(shù)值模擬方法，并編制相應(yīng)的數(shù)值計(jì)算程序。分別采用計(jì)算流體力學(xué)的VOF(volume of fluid)模型和本文中提出的方法對(duì)縫洞網(wǎng)絡(luò)模型算例進(jìn)行計(jì)算，驗(yàn)證了本文方法的正確性。

縫洞型油藏；離散縫洞網(wǎng)絡(luò)模型；油水兩相流；重力分異；流相流動(dòng)分?jǐn)?shù)

碳酸鹽巖油氣藏在現(xiàn)今石油天然氣工業(yè)中具有重要的地位。碳酸鹽巖占全球沉積巖面積的20%，相應(yīng)油氣儲(chǔ)量約占總儲(chǔ)量的50%，而油氣產(chǎn)量則約占總產(chǎn)量的60%[1]。碳酸鹽巖油氣藏中的縫洞型油藏[2]，主要以大型洞穴、溶蝕孔洞和各種裂縫作為儲(chǔ)滲空間[3]，表現(xiàn)出極強(qiáng)的非均質(zhì)性，儲(chǔ)層中流體流動(dòng)狀態(tài)十分復(fù)雜[4-5]。離散縫洞網(wǎng)絡(luò)模型[6-7]是針對(duì)縫洞型油藏非均質(zhì)性極強(qiáng)的特點(diǎn)而提出的一種數(shù)值計(jì)算模型，它將儲(chǔ)層看作裂縫介質(zhì)和溶洞介質(zhì)構(gòu)成的網(wǎng)絡(luò)系統(tǒng)，采用Darcy方程模擬裂縫中流體的流動(dòng)和采用Stokes方程模擬溶洞中流體的流動(dòng)，并在界面上進(jìn)行耦合，即Darcy-Stokes耦合計(jì)算方法。當(dāng)表征單元體不存在時(shí)，這種模型能準(zhǔn)確地對(duì)縫洞網(wǎng)絡(luò)進(jìn)行模擬[8]。Popov等將網(wǎng)格粗化方法與離散縫洞網(wǎng)絡(luò)模型相結(jié)合，采用Stokes-Brinkmann方程作為控制方程，以避免滲流區(qū)和自由流動(dòng)區(qū)在界面上耦合帶來的數(shù)值困難[9-13]。采用Stokes方程模擬溶洞中的流體流動(dòng)，現(xiàn)有的離散縫洞網(wǎng)絡(luò)模型Darcy-Stokes耦合計(jì)算方法主要用于單相流體流動(dòng)問題的數(shù)值模擬。Chen等將Cahn-Hilliard公式引入Darcy-Stokes耦合系統(tǒng)，實(shí)現(xiàn)了對(duì)兩相流問題的數(shù)值模擬[14-15]，然而這僅限于處理單一界面上的耦合問題，其原因在于，采用該耦合系統(tǒng)模擬縫洞網(wǎng)絡(luò)各溶洞中的流體流動(dòng)并計(jì)算溶洞內(nèi)油水兩相流的流態(tài)時(shí)，計(jì)算量非常大。因此，Chen等提出的方法目前還難以用于縫洞網(wǎng)絡(luò)中兩相流體流動(dòng)問題的數(shù)值模擬。對(duì)溶洞中兩相流態(tài)的分析表明，在較大的溶蝕孔洞中，毛管力可以忽略，流速低時(shí)溶洞內(nèi)的油水兩相應(yīng)處于重力分異的狀態(tài)[16]。因此，筆者基于離散縫洞網(wǎng)絡(luò)模型和溶洞內(nèi)油水兩相流體重力分異的假定，采用流相流動(dòng)分?jǐn)?shù)描述離散縫洞中油水兩相的流動(dòng)，提出一種離散縫洞網(wǎng)絡(luò)油水兩相流動(dòng)的計(jì)算模型。

1 離散縫洞網(wǎng)絡(luò)油水兩相流動(dòng)模型

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

將縫洞型碳酸鹽巖儲(chǔ)層設(shè)定為等溫、處于熱力學(xué)瞬時(shí)平衡的狀態(tài)，并且包含油水兩相流體，水和油這兩種液體組分僅存在于其伴生相中。油水兩組分的物質(zhì)守恒方程由式(1)和式(2)來描述；β相流體(β為w時(shí)代表水，β為o時(shí)代表油)在壓力、重力和毛管力的作用下，按照式(3)的Darcy定律流動(dòng)。

(1)

(2)

(3)

Φβ=pβ-ρβgD。

(4)

式(1)～(4)中：φ為儲(chǔ)層的有效孔隙度；Sβ為β相的飽和度；ρβ為β相在油藏條件下的密度；vβ為β相的滲流速度；qβ為地層β組分每單位體積匯點(diǎn)/源點(diǎn)項(xiàng)；k為儲(chǔ)層的絕對(duì)滲透率；krβ為β相的相對(duì)滲透率；μβ為β相的黏度；Φβ為β相的流動(dòng)勢(shì)；pβ為β相的壓力；g為重力加速度；D為深度。

式(1)～(3)、飽和度約束方程和本構(gòu)方程構(gòu)成了基本的流動(dòng)控制方程。

1.2 數(shù)值離散

對(duì)式(1)和式(2)采用有限體積法進(jìn)行空間離散，采用向后一階差分進(jìn)行時(shí)間離散，可得：

(5)

式(5)中：M為β相的質(zhì)量；n表示前一時(shí)刻,n+1表示當(dāng)前時(shí)刻；Vi為單元i的體積；Δt為時(shí)間步長；ηi為由所有同單元i相連接的單元j組成的集合；Fβ,ij為單元i同單元j之間β相的質(zhì)量流動(dòng)項(xiàng)，計(jì)算質(zhì)量流動(dòng)項(xiàng)時(shí)相對(duì)滲透率kr β采用單點(diǎn)上游權(quán)處理法；Qβ i為單元i內(nèi)β相的源/匯項(xiàng)。

式(5)適用于單重或多重連續(xù)介質(zhì)的網(wǎng)格，也適用于離散縫洞網(wǎng)絡(luò)模型中的裂縫單元或溶洞單元，可通過Newton-Raphson方法進(jìn)行求解。

2 溶洞單元的流相流動(dòng)分?jǐn)?shù)

對(duì)于離散縫洞網(wǎng)絡(luò)模型油水兩相流動(dòng)問題，裂縫中的流動(dòng)可采用Darcy定律來描述；溶洞內(nèi)部兩相流體的流動(dòng)屬于自由流動(dòng)，是數(shù)值計(jì)算模型的關(guān)鍵問題。

設(shè)定碳酸鹽巖縫洞型儲(chǔ)層的較大尺度溶洞內(nèi)油水兩相流體瞬時(shí)平衡，并處于上油下水的重力分異狀態(tài)，兩相間界面清晰，不存在過渡帶且相間毛管力可以忽略?；诖酥亓Ψ之惖募俣?，可將溶洞單元處理成單一網(wǎng)格，并通過溶洞單元的流相飽和度追蹤溶洞中流相界面的變化。

基于重力分異假定計(jì)算溶洞內(nèi)多相流體運(yùn)動(dòng)時(shí)，式(1)、式(2)及式(5)仍然都成立，式(5)中流動(dòng)項(xiàng)Fβ,ij可表示為：

Fβ,ij=ρβλβ,ij+1/2γij[Φβ j-Φβ i]。

(6)

式(6)中：λβ,ij+1/2為β相的流度；γij為單元i和單元j交界面處的傳導(dǎo)系數(shù)。

(7)

(8)

圖1 流相流動(dòng)分?jǐn)?shù)取值的示意Fig.1 Schematic of flow fraction

圖2 單元i到單元j的流相流動(dòng)分?jǐn)?shù)Fig.2 Flow fraction of the connection from element i to j

3 算例分析

3.1 驗(yàn)證算例

假設(shè)一個(gè)雙洞串聯(lián)縫注縫采模型算例，兩個(gè)矩形溶洞的尺寸相同，由一裂縫相連，左側(cè)溶洞由裂縫連接到注入井，右側(cè)溶洞由裂縫連接到生產(chǎn)井。初始時(shí)刻模型油相飽和，注入井以恒定速率0.144 m3/d注水。

分別采用3種方法對(duì)此算例進(jìn)行模擬計(jì)算。第一種是傳統(tǒng)的基于連續(xù)介質(zhì)的數(shù)值模擬方法，即將溶洞劃分成細(xì)密的網(wǎng)格，并賦予超大的滲透率；第二種是采用計(jì)算流體力學(xué)的VOF(volume of fluid)兩相流動(dòng)模型進(jìn)行計(jì)算；第三種是本文基于重力分異假定、采用流相流動(dòng)分?jǐn)?shù)控制的計(jì)算方法。采用第一種和第二種計(jì)算方法模擬時(shí)，使用結(jié)構(gòu)化的正方形網(wǎng)格對(duì)研究區(qū)域進(jìn)行剖分，單元總數(shù)為19 440個(gè)。

圖3 雙洞串聯(lián)縫注縫采模型算例的含水率曲線Fig.3 Water-cut curve of dual-vug model

圖3給出了3種不同計(jì)算方法對(duì)應(yīng)的生產(chǎn)井含水率隨時(shí)間的變化曲線。圖4為由VOF油水兩相流模型和傳統(tǒng)基于連續(xù)介質(zhì)滲流的數(shù)值模擬方法計(jì)算得到的同一時(shí)刻的油相飽和度分布圖(t=0.048 h)，此時(shí)右側(cè)溶洞見水。計(jì)算結(jié)果表明，離散縫洞網(wǎng)絡(luò)兩相流動(dòng)模型的計(jì)算結(jié)果與VOF油水兩相流模型的計(jì)算結(jié)果基本一致，由此驗(yàn)證了該方法和模擬程序的正確性。傳統(tǒng)基于連續(xù)介質(zhì)的數(shù)值模擬方法不能反映油水兩相在溶洞內(nèi)流動(dòng)時(shí)重力引起的分異作用，無水采油期很短。VOF方法雖然能夠準(zhǔn)確地反映大溶洞中油水界面的變化，但是使用的網(wǎng)格數(shù)量較大，計(jì)算耗時(shí)較長，難以處理溶洞數(shù)量較多的算例。本文提出的計(jì)算方法只需用2個(gè)網(wǎng)格即可準(zhǔn)確模擬生產(chǎn)井含水率，以及采出程度的變化特征，這大大降低了計(jì)算量，因而計(jì)算效率較高。

圖4 油相飽和度分布Fig.4 Distribution of oil saturation

3.2 縫洞網(wǎng)絡(luò)模型算例

假設(shè)如圖5所示的縫洞網(wǎng)絡(luò)的油水兩相流動(dòng)問題，忽略基質(zhì)儲(chǔ)滲能力。模型初始為油相飽和，水從下部溶洞1定流量注入，每天的注水量為溶洞總體積的千分之一，從上部溶洞7定壓采出。

圖5 縫洞網(wǎng)絡(luò)模型算例示意Fig.5 Schematic of a fracture-vug network model

采用本文提出的方法對(duì)該縫洞網(wǎng)絡(luò)模型算例進(jìn)行數(shù)值模擬。圖6給出了計(jì)算而得的生產(chǎn)井含水率曲線和日采油率曲線。計(jì)算結(jié)果表明，與雙洞串聯(lián)縫注縫采模型算例一樣，縫洞網(wǎng)絡(luò)模型算例中同樣出現(xiàn)了無水采油期后生產(chǎn)井含水率短時(shí)間內(nèi)迅速上升到1.0的情況。這種暴性水淹現(xiàn)象在實(shí)際縫洞型油藏的開采過程中也較為常見[17]。

模擬結(jié)束時(shí)(約1 500 d)各溶洞內(nèi)水相飽和度如表1所示。表1中還給出了各溶洞流體流出方向的連接中水相流動(dòng)分?jǐn)?shù)為1.0時(shí)水相飽和度的最小值，同圖5縫洞網(wǎng)絡(luò)幾何構(gòu)形相對(duì)比，各溶洞內(nèi)水相飽和度完全滿足預(yù)先設(shè)定的溶洞內(nèi)油水兩相重力分異的假定。

圖6 縫洞網(wǎng)絡(luò)模型算例模擬結(jié)果Fig.6 Numerical result of the fracture-vug network model

溶洞編號(hào)模擬結(jié)束時(shí)水相飽和度水相流動(dòng)分?jǐn)?shù)為1時(shí)水相飽和度10.9470.9520.9971.0030.8450.8540.9671.0050.7880.8060.9951.0070.9901.00

計(jì)算結(jié)果也表明，在模擬結(jié)束時(shí)各裂縫內(nèi)僅有水相流動(dòng)，溶洞內(nèi)油相體積幾乎不再變化，剩余油均為無法被采出的閣樓油，這同文獻(xiàn)[18]由實(shí)驗(yàn)得到的結(jié)論相吻合。當(dāng)生產(chǎn)井見水后，由于裂縫尺度相比溶洞尺度小幾個(gè)數(shù)量級(jí)，生產(chǎn)井與溶洞的連接內(nèi)水相流動(dòng)分?jǐn)?shù)會(huì)很快隨著溶洞水相飽和度的微小增加而迅速接近或達(dá)到1.0，因此生產(chǎn)井出現(xiàn)了明顯的暴性水淹現(xiàn)象。

4 結(jié) 論

筆者基于溶洞內(nèi)油水兩相流體瞬時(shí)平衡和重力分異的假定，提出了一種離散縫洞網(wǎng)絡(luò)兩相流動(dòng)的數(shù)值模擬方法，采用流相流動(dòng)分?jǐn)?shù)來模擬離散縫洞網(wǎng)絡(luò)中油水兩相的流動(dòng)，并編制相應(yīng)的數(shù)值計(jì)算程序。分別采用了傳統(tǒng)基于連續(xù)介質(zhì)的數(shù)值模擬方法(對(duì)溶洞細(xì)分網(wǎng)格并設(shè)定超大滲透率)、計(jì)算流體力學(xué)的VOF兩相流動(dòng)模型和本文建議的算法，對(duì)一個(gè)雙洞串聯(lián)縫注縫采模型算例進(jìn)行了計(jì)算。結(jié)果表明，傳統(tǒng)基于連續(xù)介質(zhì)的數(shù)值模擬方法不能很好地模擬縫洞網(wǎng)絡(luò)的流體流動(dòng)特征，無水采油期很短。本文離散縫洞網(wǎng)絡(luò)兩相流動(dòng)模型的計(jì)算結(jié)果與VOF油水兩相流模型的計(jì)算結(jié)果基本一致，且只需用兩個(gè)網(wǎng)格即可準(zhǔn)確模擬采出井含水率及采出程度的變化特征，這大大降低了計(jì)算量。對(duì)縫洞網(wǎng)絡(luò)模型算例的計(jì)算結(jié)果表明，各溶洞中水飽和度的分布符合縫洞分布的幾何構(gòu)形，完全滿足預(yù)先設(shè)定的溶洞內(nèi)油水兩相重力分異的假定；無水采油期后，生產(chǎn)井含水率短時(shí)間內(nèi)迅速上升，出現(xiàn)了明顯的暴性水淹現(xiàn)象。

[1] 江懷友，宋新民，王元基，等.世界海相碳酸鹽巖油氣勘探開發(fā)現(xiàn)狀與展望[J].海洋石油,2008,28(4):7.

[2] 林忠民.塔河油田奧陶系碳酸鹽巖儲(chǔ)層特征及成藏條件[J].石油學(xué)報(bào),2002,23(3):23.

[3] 李江龍，陳志海，高樹生.縫洞型碳酸鹽巖油藏水驅(qū)油微觀實(shí)驗(yàn)?zāi)M研究：以塔河油田為例[J].石油實(shí)驗(yàn)地質(zhì),2009,31(6):637.

[4] 康志江，李江龍，張冬麗，等.塔河縫洞型碳酸鹽巖油藏滲流特征[J].石油與天然氣地質(zhì),2005,26(5):640.

[5] 李超，李龍龍，汪洋，等.復(fù)雜裂縫性碳酸鹽巖油藏?cái)?shù)值模擬方法研究[J].長江大學(xué)學(xué)報(bào)(自然科學(xué)版),2015,12(20):65.

[6] ARBOGAST T, BRUNSON D S. A computational method for approximating a Darcy-Stokes system governing a vuggy porous medium[J]. Computational Geosciences, 2007, 11(3):207.

[7] 姚軍，黃朝琴，王子勝，等.縫洞型油藏的離散縫洞網(wǎng)絡(luò)流動(dòng)數(shù)學(xué)模型[J].石油學(xué)報(bào),2010,31(5):815.

[8] BERKOWITZ B. Characterizing flow and transport in fractured geological media: a review[J]. Advances in Water Resources, 2002, 25(8):869.

[9] POPOV P, QIN G, BI L, et al. Multiphysics and multiscale methods for modeling fluid flow through naturally fractured carbonate karst reservoirs[J]. SPE Reservoir Evaluation and Engineering, 2009, 12(2): 218.

[10] POPOV P, EFENDIEV Y, QIN G. Multiscale modeling and simulations of flows in naturally fractured karst reservoirs[J]. Communications in Computational Physics, 2009, 6(1):162.

[11] GULBRANSEN A, HAUGE V L, LIE K A. A multiscale mixed finite element method for vuggy and naturally fractured reservoirs[J]. SPE Journal, 2009, 15(15):395.

[12] HE J, KILLOUGH J E, FADLElMULA F, et al. A unified finite difference model for the simulation of transient flow in naturally fractured carbonate karst reservoirs[C]//SPE Reservoir Simulation Symposium. Houston, Texas: Society of Petroleum Engineers, 2015.

[13] HE J, KILLOUGH J E, GAO S, et al. Confronting the simulation of fluid flow in naturally fractured carbonate karst reservoirs[C]//Abu Dhabi International Petroleum Exhibition and Conference. Abu Dhabi: Society of Petroleum Engineers, 2016.

[14] CHEN J, SUN S, WANG X P. A numerical method for a model of two-phase flow in a coupled free flow and porous media system[J]. Journal of Computational Physics, 2014, 268(2):1.

[15] XIE H, LI A, HUANG Z, et al. Coupling of two-phase flow in fractured-vuggy reservoir with filling medium[J]. Open Physics, 2017, 15(1):12.

[16] 代亞非.縫洞型碳酸鹽巖油藏?cái)?shù)值模擬[D].北京:北京大學(xué),2014.

[17] 李小波，榮元帥，龍喜彬，等.縫洞型油藏強(qiáng)邊、底水竄進(jìn)油井特征及機(jī)理研究[J].西南石油大學(xué)學(xué)報(bào)(自然科學(xué)版),2015,37(1):136.

[18] 劉中春.塔河縫洞型油藏剩余油分析與提高采收率途徑[J].大慶石油地質(zhì)與開發(fā),2015,34(2):62.

Numerical simulation of oil-water two-phase flow
in the discrete fracture-vug network

SONG Shihan1, DI Yuan1, LIU Zhongchun2, DAI Yafei3

(1. College of Engineering, Peking University, Beijing 100871, China; 2. Exploration and Production Research Institute, China Petroleum and Chemical Corporation, Beijing 100083, China; 3. Envision Energy Limited, Shanghai 200051, China)

Based on the assumption of transient equilibrium and gravity segregation of oil-water two-phase flow in major vugs, a numerical method is proposed to simulate the two-phase flow in the discrete fracture-vug network by using flow fraction, with the corresponding numerical calculation program developed. Several numerical examples of the fracture-vug network are calculated by applying the model of volume of fluid(VOF) and the proposed method of the article respectively, whose numerical results have proved validity of the proposed method.

fractured-vuggy reservoir; discrete fracture-vug network model; oil-water two-phase flow; gravity segregation; flow fraction

10.3969/j.issn.1671-8798.2017.03.003

2017-04-02

國家科技重大專項(xiàng)項(xiàng)目(2016ZX05014)；國家自然科學(xué)基金項(xiàng)目(51674010)

邸 元(1968— )，男，陜西省西安人，副教授，博士，主要從事油藏?cái)?shù)值模擬研究。E-mail：diyuan@mech.pku.edu.cn。

TE344；O359

A

1671-8798(2017)03-0172-06