張娜 王少椿 李立 孫乾 劉偉

摘要：擬有限差分方法因具有良好的局部守恒性且適用于任意復雜網(wǎng)格系統(tǒng)，在流體力學等領(lǐng)域得到應用。建立基于擬有限差分法的裂縫性介質(zhì)兩相流模型，詳細闡述擬有限差分方法的基本原理，并與離散裂縫技術(shù)相結(jié)合推導得到考慮毛管力的裂縫性油水兩相滲流問題擬有限差分計算格式；采用全隱式方法對其兩相流問題進行求解，克服傳統(tǒng)求解方法不能有效處理全張量滲透率的缺陷。為了驗證方法的正確性，對比不同濕潤狀態(tài)下的離散裂縫全隱式擬有限差分方法準確性和收斂性，并將該方法應用于油田現(xiàn)場實例。結(jié)果表明：當考慮毛管力時水驅(qū)替前緣緩慢穿過基質(zhì)和裂縫，然后到達生產(chǎn)井；忽略毛管力時水在基質(zhì)中的流動相對緩慢，流體在裂縫中流動迅速，快速到達生產(chǎn)井；該方法能夠正確地模擬不同潤濕狀態(tài)和不同裂縫位置下的毛管力和重力情況。

關(guān)鍵詞：自發(fā)滲吸； 擬有限差分； 毛管力； 全張量滲透率； 離散裂縫

中圖分類號：TE 122.14 文獻標志碼：A

引用格式：張娜，王少椿，李立，等.考慮毛管力的全張量裂縫性介質(zhì)兩相流擬有限差分模擬［J］.中國石油大學學報（自然科學版），2023，47（1）：98-105.

ZHANG Na， WANG Shaochun， LI Li， et al. A mimetic finite difference simulation method of multiphase flow in fractured reservoir considering capillarity and full tensor permeability［J］.Journal of China University of Petroleum（Edition of Natural Science），2023，47（1）：98-105.

A mimetic finite difference simulation method of multiphase flow in

fractured reservoir considering capillarity and full tensor permeability

ZHANG Na¹， WANG Shaochun¹， LI Li¹， SUN Qian²， LIU Wei¹

（1.College of Energy Resources ，Chengdu University of Technology， Chengdu 610059， China；

2.College of Petroleum Engineering， Shandong Institute of Petroleum and Chemical Technology， Dongying 257061， China）

Abstract：The mimetic finite difference method has been applied in the field of fluid mechanics， because it has good local conservation and is suitable for complex grids. In this paper， a two-phase flow model for fractured media based on the mimetic finite difference method was established， and the basic principles of the mimetic finite difference method were described in detail. In the new model， the capillarity was considered， and the corresponding scheme of fractured two-phase flow in fractured reservoirs was deduced. A full implicit scheme was used to solve the two-phase flow problem. The new method can overcome the shortcomings of the traditional method that cannot effectively deal with the full tensor permeability. In order to verify the correctness of the method and the robustness of the program， the accuracy and convergence of the mimetic finite difference method for fractured media under different wetting conditions were analyzed， and the method was also applied in a field case study. The results show that， when capillary pressure is considered， the water displacement front slowly passes through the matrix and fractures before reaching to the producing well， while the capillary pressure is ignored， the flow of water through the matrix is relatively slow， but the fluid flows rapidly through the fracture， reaching to the producing well quickly. The new method can accurately simulate the capillary pressure and gravity under different wetting conditions and at different fracture locations.

Keywords：spontaneous imbibition; mimetic finite difference; capillarity; full tensor permeability; discrete fractures

由于裂縫-基質(zhì)體系中多相流的復雜性，裂縫性儲層已成為多孔介質(zhì)多相流領(lǐng)域的一個重要研究課題^［1-2］。因此了解控制此類儲層可采烴量的潛在物理現(xiàn)象至關(guān)重要^［3-6］。模擬裂縫性儲層通常采用雙重孔隙^［7-9］和離散裂縫方法^［10-16］兩種模型。雙重孔隙模型主要假設(shè)采用統(tǒng)一的傳遞速率，即裂縫立即被潤濕相填充^［17-18］。離散裂縫方法依賴于通過數(shù)值逼近方法將裂縫和矩陣域離散化，從而在數(shù)值復雜性方面產(chǎn)生一個更簡單的問題，但它不能用于具有復雜的互聯(lián)裂縫網(wǎng)絡(luò)的系統(tǒng)。在基巖-裂縫框架中，通常采用有限元方法對該區(qū)域進行離散，并將其分為經(jīng)典有限元法和混合有限元法。兩者主要區(qū)別在于后者是局部保守的^［19］。然而混合有限元法生成的數(shù)學系統(tǒng)通常在數(shù)值上很難求解，因此引入了混合雜交有限元方法，但該方法基函數(shù)的計算僅適用于某些特定類型的單元。為了克服這些限制，在基函數(shù)計算的基礎(chǔ)上引入擬有限差分方法對高度非結(jié)構(gòu)化多邊形進行建模，拓寬了其適用范圍^［20］。筆者對不同巖石潤濕狀態(tài)下滲吸過程的流動動力學進行研究，同時將擬有限差分方法推廣至考慮滲吸過程的裂縫性介質(zhì)流動數(shù)值模擬中，詳細推導擬有限差分的非結(jié)構(gòu)網(wǎng)格計算形式，通過算例驗證方法的正確性。

1 滲流控制方程

根據(jù)質(zhì)量守恒、達西定律和狀態(tài)方程，給出多孔介質(zhì)中流體兩相流動的控制方程。假設(shè)擴散力很小，可以忽略不計。連續(xù)性方程為

式中，w為水相；o為油相;φ為基巖孔隙度，%；S_α為α相的飽和度，%；ρ_α為α相的密度，g/cm³；u_α為α相的達西速度，cm³/s；q為質(zhì)量流量，g/s。

每個階段的達西速度為

然后引入局部約束來定義系統(tǒng)，即

相流度表示為

式中，K_rα為α相的相對滲透率；μ_α為α相的黏度，mPa·s；p_α為α相的壓力，MPa；D為垂向深度，m；p_c為毛管壓力，MPa；λ_α為α相的流度；K為滲透率矩陣;g為重力加速度，m/s²。

1.1 擬有限差分數(shù)值離散

擬有限差分方法的出發(fā)點是方程中描述質(zhì)量平衡（式（1））和達西速度（式（2））的控制方程。利用向量基函數(shù)V_i來描述動量平衡?；瘮?shù)具有以下屬性：

（1）向量函數(shù)在界面i處的流通量為1，其他地方為0。

（2）基函數(shù)在單元上的散度為常數(shù)。

（3）速度場u可以用 表示，其中Q_i為界面流量，N_f為界面數(shù)。

基質(zhì)和裂縫（藍色部分）單元示意圖如圖1所示。

擬有限差分方法的關(guān)鍵元素是局部內(nèi)積W，可用于表征界面流量Q為

式中，Q=［Q₁，Q₂，…，Q_m］^T為界面流量向量；m為單元Ω_i的界面數(shù)；Λ為相對滲透率矩陣；e=（1，…，1）^T_1×m；X_i為單元Ω_i中點標量；σ為定義在單元Ω_i界面中心處的向量；W為正定矩陣，是擬有限差分法的核心部分。

線性壓力可以寫成p=X_ia+b的形式，其中a和b是常量向量和常量標量。這里考慮毛管力和重力，因此它們在中心和表面的數(shù)量滿足線性關(guān)系，即

式中，pi和pc，i分別為單元Ωi中心處的壓力和毛管壓力，Pa；σk，i為界面Ak的界面壓力，Pa；σc，k，i為界面Ak的界面毛管壓力，Pa；Di和Dz，k，i分別為單元Ωi和界面Ak中心處的深度，m；xi和xk分別為單元Ωi的中心和界面Ak的中心的坐標向量。

根據(jù)式（2），界面Ak處的流量為

式中，u_k為界面處達西流速，cm³/s；A_k為界面面積，cm²；_k為界面處單元法向量。

將式（7）～（9）和（10）代入式（6），則W滿足條件：

其中

文獻［21］中詳細討論了獲得對稱正定矩陣W的方法。在這里，W表達式為

其中

式中，d為空間維度；A為由界面面積構(gòu)成的對角矩陣；Y為單元中心到界面中心的距離組成的對角矩陣；tr（K）為K的跡。

計算得到矩陣W后，界面通量Q可以簡化為數(shù)值形式，即

其中

式中，po為油壓，Pa；pc，α為相α的毛管壓力，Pa；σo，j、σc，α，j、DD，j分別為相α在界面j上的油壓、毛管壓力和深度。

在擬有限差分方法中，需要通過增加一個由動量平衡方程表示的約束來保證界面上的流量連續(xù)性。如果某個單元的界面不與另一個單元相交共享（界面在邊界處），則假設(shè)該界面處的流量為零，即Qα，i=0。

但是，如果一個界面由兩個相鄰單元共享，則Qα，i=-Q′α，i。這意味著拉格朗日乘數(shù)在界面處也是相同的，則相α在交界面上的總通量Qt，α，i為

1.2 基于離散裂縫模型的裂縫離散化

對于離散裂縫處理，上述流量方程適用于較低空間維度的裂縫。因此裂縫體系的流量方程為

其中

式中，N_{f，F(xiàn)}為裂縫邊數(shù)；下角標F代表裂縫系統(tǒng)。

如圖1所示的單元，離散裂縫A_F被認為是單元Ω_i和Ω_j的交界面。因此裂縫單元壓力等于基質(zhì)界面壓力σ_i（σ_j）。這樣，在計算最終線性方程組時消除了裂縫界面壓力，只保留基巖單元壓力。

為了耦合基巖和裂縫兩個系統(tǒng)，定義兩個相鄰單元界面處的條件（界面A_F為單元Ω_i和Ω_j的交界面，如圖1所示）：①如果A_F是基巖單元，則流量滿足Q_i+Q_j=0；②如果A_F為裂縫，裂縫流量滿足質(zhì)量守恒方程Q_F=Q_i+Q_j+q_F，其中，q_F為裂縫的匯/源。

1.3 非線性求解

根據(jù)基質(zhì)和裂縫體系的無流動邊界條件和連續(xù)條件，可以得到殘差項為

式中，h_F為裂縫開度。

設(shè)主變量集合為y={p_o，σ_o，σ_c，σ_o，F(xiàn)，σ_c，_F，S_w，S_w，_F}，線性系統(tǒng)在每次非線性迭代的形式為

式中，R_y為質(zhì)量平衡和飽和方程的余量。因此可以通過全隱式法同時求解所有控制方程。

2 毛管力模型建立

根據(jù)Schmid等^［22］研究中提出的基質(zhì)巖石性質(zhì)的定義，基于冪律模型構(gòu)建相對滲透率曲線，表達式為

式中，K_romax和K_rwmax分別為油和水的最大相對滲透率；n和m為相對滲透率指數(shù)；S_wi為初始含水飽和度；S_or為殘余油飽和度；K_ro和K_rw分別為油和水的相對滲透率。

此外，通用混濕體系^［23］毛管壓力預測模型為

式中，p_c為毛管壓力；a和c為常量，表示排水或自吸過程。除了c_o，所有常數(shù)都是正的。毛管壓力曲線由兩條漸近分支組成。表1為本文中采用的不同濕潤情況下的測試參數(shù)。圖2為相應的相對滲透率和毛管曲線。

3 算例分析

3.1 模型驗證

考慮一個均勻基巖模型，20 m×20 m×100 m的三維域。最初該區(qū)域被等體積的油和水所飽和；水相位于油相之上。該區(qū)域完全密封，流體不與周圍相互作用，從而消除了任何外部流動影響。由于重力和油水相的毛管壓力決定了區(qū)域內(nèi)的流動，而且毛管壓力和重力最初不平衡，水趨向于垂直向下流動，直至平衡。水柱移動的高度為

式中，H為水柱移動高度，m。

為了驗證區(qū)域內(nèi)自吸流動行為，采用不同模型進行驗證（圖3）。三裂縫模型中兩條裂縫有交叉。基巖孔隙度為0.5，裂縫孔隙度為0.75，油和水的黏度分別為2.5和0.85 mPa·s，油和水的密度分別為875和1000 kg/m³，殘余油飽和度S_or為0.1。

水柱移動高度的數(shù)值解和解析解的對比如圖4所示，而圖5為通過擬有限差分方法在三維域內(nèi)模擬產(chǎn)生的平衡態(tài)含水飽和度剖面。由數(shù)值解計算得到的水柱高度與毛管壓力曲線的解析解吻合較好。由于毛管力和重力不平衡，水向下運動，當重力和毛管力達到平衡后，含水飽和度保持不變。結(jié)果表明，該方法能夠正確地模擬不同潤濕狀態(tài)和不同裂縫位置下的毛管力和重力情況。

3.2 實例應用

裂縫模型是根據(jù)W油田裂縫儲層構(gòu)建，由離散裂縫構(gòu)成，如圖6所示。模型尺寸為486 m×684 m×10 m，包含419個離散裂縫，共有16131個四面體單元和37654個界面，其中1260個三角形單元表示離散裂縫。以10000 m³/d的恒定速率生產(chǎn)和注入，相對滲透率和毛管壓力曲線如圖2中弱水濕對應的曲線。基巖和裂縫滲透率張量K_m和K_f分別為

圖7為注入0.55V_P（V_P為孔隙體積）時含水飽和度分布。當考慮毛管力時，水驅(qū)替前緣緩慢穿過基質(zhì)和裂縫，然后到達生產(chǎn)井。忽略毛管力時水在基質(zhì)中的流動相對緩慢，流體在裂縫中流動迅速，快速到達生產(chǎn)井。因此毛管壓力對現(xiàn)場原油采收率和產(chǎn)量的影響顯著。此外，引入毛管壓力并不影響生產(chǎn)井的最終含水率，而是控制水突破的速度，從圖8中的含水率曲線可以看到這一點，含水率的最終值受路徑和水突破點的影響。

此外，考慮毛管力作用和不考慮毛管力作用的平均計算時間分別為9.571×10⁵和3.346×10⁵s。當考慮毛管力作用時，方程和非線性迭代的規(guī)模增大；因此運行該算法所需的時間迅速增加約65%。毛管壓力作用對該問題造成了物理上的復雜性，尤其是存在裂縫的情況下，因此需要更多的計算來收斂到物理解。

4 結(jié) 論

（1）根據(jù)裂縫性油藏兩相流特征，建立了擬有限差分法和離散裂縫技術(shù)相結(jié)合的新型全隱式數(shù)值方法，并通過數(shù)值算例驗證方法的正確性。擬有限差分法自身良好的局部守恒性和對任意復雜網(wǎng)格的適用性，將其用于考慮毛管力的裂縫性油藏兩相流數(shù)值模擬，有效克服傳統(tǒng)方法不能處理全張量滲透率的弊端。

（2）當考慮毛管力時，水驅(qū)替前緣緩慢穿過基質(zhì)和裂縫，然后到達生產(chǎn)井；忽略毛管力時，水在基質(zhì)中的流動相對緩慢，流體在裂縫中流動迅速，快速到達生產(chǎn)井；該方法能夠模擬不同潤濕狀態(tài)和不同裂縫位置下的毛管力和重力情況。

參考文獻：

［1］王敬，劉慧卿，夏靜，等.裂縫性油藏滲吸采油機理數(shù)值模擬［J］.石油勘探與開發(fā)，2017，44（5）：761-770.

WANG Jing， LIU Huiqing， XIA Jing， et al. Mechanism simulation of oil displacement by imbibition in fractured reservoirs［J］. Petroleum Exploration and Development， 2017，44（5）：761-770.

［2］唐潮，陳小凡，杜志敏，等.基于有限體積法的裂縫性油藏兩相流動模型［J］.石油學報，2018，39（8）：924-936.

TANG Chao， CHEN Xiaofan， DU Zhimin， et al. Fractured reservoir two-phase flow model based on finite volume method［J］. Acta Petrolei Sinica，2018，39（8）：924-936.

［3］AUSTAD T， MILTER J. Spontaneous imbibition of water into low permeable chalk at different wettabilities using surfactants［J］. Oil Field， 1997，24（5）：257-266.

［4］何志剛.毛管力自吸采油［J］.科技導報，2011，29（4）：39-43.

HE Zhigang. Oil by capillary sopontanueous imbibitions［J］. Science & Technology Review，2011，29（4）：39-43.

［5］孟智強，凌浩川，朱志強，等.裂縫性油藏吸入毛管壓力曲線研究及應用［J］.西安石油大學學報（自然科學版），2019，34（2）：60-65.

MENG Zhiqiang， LING Haochuan， ZHU Zhiqiang， et al. Research and application of sucking capillary pressure curve of fracture reservoir［J］.Journal of Xi'an Shiyou University（Natural Science Edition），2019，34（2）：60-65.

［6］MORROWN R， MASON G. Recovery of oil by spontaneous imbibition［J］.Current Opinion in Colloid and Interface Science， 2001，4356（4）：321 - 337.

［7］BAARENBLATT G， ZHELTOV I P， KOCHINA I. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks［J］. Journal of Applied Mathematics and Mechanics， 1960，24（5）：1286-1303.

［8］WARREN J， ROOT P J. The behavior of naturally fractured reservoirs［J］. SPE Journal， 1963，3（3）：245-255.

［9］KAZEMI H， LS M， PORTERFIELD K， et al. Numerical simulation of water-oil flow in naturally fractured reservoirs［J］. SPE Journal， 1976，16（6）：317-326.

［10］MAHDI A， MOHAMMADM M， MOHAMMADS S， et al. Fluid flow in fractured reservoirs： exact analytical solution for transient dual porosity model with variable rock matrix block size［J］. Journal of Petroleum Science and Engineering，2018，164（3）：571-583.

［11］LOUIS C， WITTKE W. Experimental study of water flow in jointed rock massif［J］. Tachien Project Formosa Geotechnique， 1971，21（1）：29-35.

［12］WILSON C R， WITHERSPOON P A. Steady state flow in rigid networks of fractures［J］. Water Resources Research， 1974，10（2）：328-339.

［13］孫致學，姜傳胤，張凱，等.基于離散裂縫模型的CO₂增強型地熱系統(tǒng)THM耦合數(shù)值模擬［J］.中國石油大學學報（自然科學版），2020，44（6）：79-87.

SUN Zhixue， JIANG Chuanyu， ZHANG Kai， et al. Numerical simulation for heat extraction of CO₂-EGS with thermal-hydraulic-mechanical coupling method based on discrete fracture models［J］. Journal of China University of Petroleum （Edition of Natural Science）， 2020，44（6）：79-87.

［14］徐榮利，郭天魁，曲占慶，等.基于離散裂縫模型的頁巖油儲層壓裂滲吸數(shù)值模擬［J］.工程科學學報，2022，44（3）：451-463.

XU Rongli， GUO Tiankui， QU Zhanqing， et al. Numerical simulation of fractured imbibition in a shale oil reservoir based on the discrete fracture model［J］. Chinese Journal of Engineering， 2022，44（3）：451-463.

［15］WOOTEN S O， JONES T A. Graphical method to determine deviation from the cubic law for single-phase flow through natural fractures［R］. SPE 020163，1989.

［16］糜利棟，姜漢橋，李濤，等.基于離散裂縫模型的頁巖氣動態(tài)特征分析［J］.中國石油大學學報（自然科學版），2015，39（3）：126-131.

MI Lidong， JIANG Hanqiao， LI Tao， et al. Characterization and dynamic analysis of shale gas production based on discrete fracture model［J］. Journal of China University of Petroleum （Edition of Natural Science）， 2015，39（3）：126-131.

［17］張有天，劉中.降雨過程裂隙網(wǎng)絡(luò)飽和-非飽和、非恒定滲流分析［J］.巖石力學與工程學報，1997，16（2）：104-111．

ZHANG Youtian， LIU Zhong. Analysis of saturated/unsaturated and unsteady seepage flow in fracture network during rainfall［J］ Chinese Journal of Rock Mechanics and Engineering， 1997，16（2）：104-111．

［18］DONATOG D，BLUNT M J. Streamline-based dual-porosity simulation of reactive transport and flow in fractured reservoirs［J］. Water Resources Research， 2004，40（4）：285-296.

［19］蘇皓，雷征東，李俊超，等.儲集層多尺度裂縫高效數(shù)值模擬模型［J］.石油學報，2019，40（5）：587-593.

SU Hao， LEI Zhengdong， LI Junchao， et al. An effective numerical simulation model of multi-scale fractures in reservoir［J］. Acta Petrolei Sinica， 2019，40（5）：587-593.

［20］姚軍，張娜，黃朝琴，等.非均質(zhì)油藏多尺度混合有限元數(shù)值模擬方法［J］.石油學報，2012，33（3）：442-447.

YAO Jun， ZHANG Na， HUANG Zhaoqin， et al. Multiscale mixed finite element numerical simulation for heterogeneous reservoirs［J］. Acta Petrolei Sinica，2012，33（3）：442-447.

［21］BREZZI F， LIPNIKOV K， SHASHKOV M.Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes［J］. SIAM Journal on Numerical Analysis. 2005， 43（5）：1872-1896.

［22］VEIGA L B D， GYRYA V， LIPNIKOV K， et al. Mimetic finite difference method for the stokes problem on polygonal meshes［J］.Journal of Computational Physics， 2009，228（19）：7215-7232.

［23］SCHMID K S， ALYAFEI N， GEIGER S， et al. Analytical solutions for spontaneous imbibition： fractional-flow theory and experimental analysis［J］.SPE Journal， 2016，21（6）：2308-2316.

［24］SKJAEVELAND S， SIQVELAND L， KJOSAVIK A，? et al. Capillary pressure correlation for mixed-wet reservoirs［J］.SPE Reservoir Evaluation & Engineering， 2000，3（1）：60-67.

（編輯 李志芬）