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

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

?

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

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

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

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

縫洞型油藏；離散縫洞網(wǎng)絡模型；油水兩相流；重力分異；流相流動分數(shù)

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

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

1.1 數(shù)學模型

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

(1)

(2)

(3)

Φβ=pβ-ρβgD。

(4)

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

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

1.2 數(shù)值離散

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

(5)

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

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

2 溶洞單元的流相流動分數(shù)

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

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

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

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

(6)

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

(7)

(8)

圖1 流相流動分數(shù)取值的示意Fig.1 Schematic of flow fraction

圖2 單元i到單元j的流相流動分數(shù)Fig.2 Flow fraction of the connection from element i to j

3 算例分析

3.1 驗證算例

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

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

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

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

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

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

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

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

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

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

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

溶洞編號模擬結束時水相飽和度水相流動分數(shù)為1時水相飽和度10.9470.9520.9971.0030.8450.8540.9671.0050.7880.8060.9951.0070.9901.00

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

4 結 論

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

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

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

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

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

[5] 李超，李龍龍，汪洋，等.復雜裂縫性碳酸鹽巖油藏數(shù)值模擬方法研究[J].長江大學學報(自然科學版),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)絡流動數(shù)學模型[J].石油學報,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] 代亞非.縫洞型碳酸鹽巖油藏數(shù)值模擬[D].北京:北京大學,2014.

[17] 李小波，榮元帥，龍喜彬，等.縫洞型油藏強邊、底水竄進油井特征及機理研究[J].西南石油大學學報(自然科學版),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

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

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

TE344；O359

A

1671-8798(2017)03-0172-06