考虑毛管力的全张量裂缝性介质两相流拟有限差分模拟
2023-09-05张娜王少椿李立孙乾刘伟
张娜 王少椿 李立 孙乾 刘伟
摘要：擬有限差分方法因具有良好的局部守恒性且适用于任意复杂网格系统，在流体力学等领域得到应用。建立基于拟有限差分法的裂缝性介质两相流模型，详细阐述拟有限差分方法的基本原理，并与离散裂缝技术相结合推导得到考虑毛管力的裂缝性油水两相渗流问题拟有限差分计算格式；采用全隐式方法对其两相流问题进行求解，克服传统求解方法不能有效处理全张量渗透率的缺陷。为了验证方法的正确性，对比不同湿润状态下的离散裂缝全隐式拟有限差分方法准确性和收敛性，并将该方法应用于油田现场实例。结果表明：当考虑毛管力时水驱替前缘缓慢穿过基质和裂缝，然后到达生产井；忽略毛管力时水在基质中的流动相对缓慢，流体在裂缝中流动迅速，快速到达生产井；该方法能够正确地模拟不同润湿状态和不同裂缝位置下的毛管力和重力情况。
关键词：自发渗吸； 拟有限差分； 毛管力； 全张量渗透率； 离散裂缝
中图分类号：TE 122.14 文献标志码：A
引用格式：张娜，王少椿，李立，等.考虑毛管力的全张量裂缝性介质两相流拟有限差分模拟［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 Na1， WANG Shaochun1， LI Li1， SUN Qian2， LIU Wei1
（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
由于裂缝-基质体系中多相流的复杂性，裂缝性储层已成为多孔介质多相流领域的一个重要研究课题［1-2］。因此了解控制此类储层可采烃量的潜在物理现象至关重要［3-6］。模拟裂缝性储层通常采用双重孔隙［7-9］和离散裂缝方法［10-16］两种模型。双重孔隙模型主要假设采用统一的传递速率，即裂缝立即被润湿相填充［17-18］。离散裂缝方法依赖于通过数值逼近方法将裂缝和矩阵域离散化，从而在数值复杂性方面产生一个更简单的问题，但它不能用于具有复杂的互联裂缝网络的系统。在基岩-裂缝框架中，通常采用有限元方法对该区域进行离散，并将其分为经典有限元法和混合有限元法。两者主要区别在于后者是局部保守的［19］。然而混合有限元法生成的数学系统通常在数值上很难求解，因此引入了混合杂交有限元方法，但该方法基函数的计算仅适用于某些特定类型的单元。为了克服这些限制，在基函数计算的基础上引入拟有限差分方法对高度非结构化多边形进行建模，拓宽了其适用范围［20］。筆者对不同岩石润湿状态下渗吸过程的流动动力学进行研究，同时将拟有限差分方法推广至考虑渗吸过程的裂缝性介质流动数值模拟中，详细推导拟有限差分的非结构网格计算形式，通过算例验证方法的正确性。
1 渗流控制方程
根据质量守恒、达西定律和状态方程，给出多孔介质中流体两相流动的控制方程。假设扩散力很小，可以忽略不计。连续性方程为
式中，w为水相；o为油相;φ为基岩孔隙度，%；Sα为α相的饱和度，%；ρα为α相的密度，g/cm3；uα为α相的达西速度，cm3/s；q为质量流量，g/s。
每个阶段的达西速度为
然后引入局部约束来定义系统，即
相流度表示为
式中，Krα为α相的相对渗透率；μα为α相的黏度，mPa·s；pα为α相的压力，MPa；D为垂向深度，m；pc为毛管压力，MPa；λα为α相的流度；K为渗透率矩阵;g为重力加速度，m/s2。
1.1 拟有限差分数值离散
拟有限差分方法的出发点是方程中描述质量平衡（式（1））和达西速度（式（2））的控制方程。利用向量基函数Vi来描述动量平衡。基函数具有以下属性：
（1）向量函数在界面i处的流通量为1，其他地方为0。
（2）基函数在单元上的散度为常数。
（3）速度场u可以用 表示，其中Qi为界面流量，Nf为界面数。
基质和裂缝（蓝色部分）单元示意图如图1所示。
拟有限差分方法的关键元素是局部内积W，可用于表征界面流量Q为
式中，Q=［Q1，Q2，…，Qm］T为界面流量向量；m为单元Ωi的界面数；Λ为相对渗透率矩阵；e=（1，…，1）T1×m；Xi为单元Ωi中点标量；σ为定义在单元Ωi界面中心处的向量；W为正定矩阵，是拟有限差分法的核心部分。
线性压力可以写成p=Xia+b的形式，其中a和b是常量向量和常量标量。这里考虑毛管力和重力，因此它们在中心和表面的数量满足线性关系，即
式中，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的中心的坐标向量。
根据式（2），界面Ak处的流量为
式中，uk为界面处达西流速，cm3/s；Ak为界面面积，cm2；k为界面处单元法向量。
将式（7）～（9）和（10）代入式（6），则W满足条件：
其中
文献［21］中详细讨论了获得对称正定矩阵W的方法。在这里，W表达式为
其中
式中，d为空间维度；A为由界面面积构成的对角矩阵；Y为单元中心到界面中心的距离组成的对角矩阵；tr（K）为K的迹。
计算得到矩阵W后，界面通量Q可以简化为数值形式，即
其中
式中，po为油压，Pa；pc，α为相α的毛管压力，Pa；σo，j、σc，α，j、DD，j分别为相α在界面j上的油压、毛管压力和深度。
在拟有限差分方法中，需要通过增加一个由动量平衡方程表示的约束来保证界面上的流量连续性。如果某个单元的界面不与另一个单元相交共享（界面在边界处），则假设该界面处的流量为零，即Qα，i=0。
但是，如果一个界面由两个相邻单元共享，则Qα，i=-Q′α，i。这意味着拉格朗日乘数在界面处也是相同的，则相α在交界面上的总通量Qt，α，i为
1.2 基于离散裂缝模型的裂缝离散化
对于离散裂缝处理，上述流量方程适用于较低空间维度的裂缝。因此裂缝体系的流量方程为
其中
式中，Nf，F为裂缝边数；下角标F代表裂缝系统。
如图1所示的单元，离散裂缝AF被认为是单元Ωi和Ωj的交界面。因此裂缝单元压力等于基质界面压力σi（σj）。这样，在计算最终线性方程组时消除了裂缝界面压力，只保留基岩单元压力。
为了耦合基岩和裂缝两个系统，定义两个相邻单元界面处的条件（界面AF为单元Ωi和Ωj的交界面，如图1所示）：①如果AF是基岩单元，则流量满足Qi+Qj=0；②如果AF为裂缝，裂缝流量满足质量守恒方程QF=Qi+Qj+qF，其中，qF为裂缝的汇/源。
1.3 非线性求解
根据基质和裂缝体系的无流动边界条件和连续条件，可以得到残差项为
式中，hF为裂缝开度。
設主变量集合为y={po，σo，σc，σo，F，σc，F，Sw，Sw，F}，线性系统在每次非线性迭代的形式为
式中，Ry为质量平衡和饱和方程的余量。因此可以通过全隐式法同时求解所有控制方程。
2 毛管力模型建立
根据Schmid等［22］研究中提出的基质岩石性质的定义，基于幂律模型构建相对渗透率曲线，表达式为
式中，Kromax和Krwmax分别为油和水的最大相对渗透率；n和m为相对渗透率指数；Swi为初始含水饱和度；Sor为残余油饱和度；Kro和Krw分别为油和水的相对渗透率。
此外，通用混湿体系［23］毛管压力预测模型为
式中，pc为毛管压力；a和c为常量，表示排水或自吸过程。除了co，所有常数都是正的。毛管压力曲线由两条渐近分支组成。表1为本文中采用的不同湿润情况下的测试参数。图2为相应的相对渗透率和毛管曲线。
3 算例分析
3.1 模型验证
考虑一个均匀基岩模型，20 m×20 m×100 m的三维域。最初该区域被等体积的油和水所饱和；水相位于油相之上。该区域完全密封，流体不与周围相互作用，从而消除了任何外部流动影响。由于重力和油水相的毛管压力决定了区域内的流动，而且毛管压力和重力最初不平衡，水趋向于垂直向下流动，直至平衡。水柱移动的高度为
式中，H为水柱移动高度，m。
为了验证区域内自吸流动行为，采用不同模型进行验证（图3）。三裂缝模型中两条裂缝有交叉。基岩孔隙度为0.5，裂缝孔隙度为0.75，油和水的黏度分别为2.5和0.85 mPa·s，油和水的密度分别为875和1000 kg/m3，残余油饱和度Sor为0.1。
水柱移动高度的数值解和解析解的对比如图4所示，而图5为通过拟有限差分方法在三维域内模拟产生的平衡态含水饱和度剖面。由数值解计算得到的水柱高度与毛管压力曲线的解析解吻合较好。由于毛管力和重力不平衡，水向下运动，当重力和毛管力达到平衡后，含水饱和度保持不变。结果表明，该方法能够正确地模拟不同润湿状态和不同裂缝位置下的毛管力和重力情况。
3.2 实例应用
裂缝模型是根据W油田裂缝储层构建，由离散裂缝构成，如图6所示。模型尺寸为486 m×684 m×10 m，包含419个离散裂缝，共有16131个四面体单元和37654个界面，其中1260个三角形单元表示离散裂缝。以10000 m3/d的恒定速率生产和注入，相对渗透率和毛管压力曲线如图2中弱水湿对应的曲线。基岩和裂缝渗透率张量Km和Kf分别为
图7为注入0.55VP（VP为孔隙体积）时含水饱和度分布。当考虑毛管力时，水驱替前缘缓慢穿过基质和裂缝，然后到达生产井。忽略毛管力时水在基质中的流动相对缓慢，流体在裂缝中流动迅速，快速到达生产井。因此毛管压力对现场原油采收率和产量的影响显著。此外，引入毛管压力并不影响生产井的最终含水率，而是控制水突破的速度，从图8中的含水率曲线可以看到这一点，含水率的最终值受路径和水突破点的影响。
此外，考虑毛管力作用和不考虑毛管力作用的平均计算时间分别为9.571×105和3.346×105s。当考虑毛管力作用时，方程和非线性迭代的规模增大；因此运行该算法所需的时间迅速增加约65%。毛管压力作用对该问题造成了物理上的复杂性，尤其是存在裂缝的情况下，因此需要更多的计算来收敛到物理解。
4 結 论
（1）根据裂缝性油藏两相流特征，建立了拟有限差分法和离散裂缝技术相结合的新型全隐式数值方法，并通过数值算例验证方法的正确性。拟有限差分法自身良好的局部守恒性和对任意复杂网格的适用性，将其用于考虑毛管力的裂缝性油藏两相流数值模拟，有效克服传统方法不能处理全张量渗透率的弊端。
（2）当考虑毛管力时，水驱替前缘缓慢穿过基质和裂缝，然后到达生产井；忽略毛管力时，水在基质中的流动相对缓慢，流体在裂缝中流动迅速，快速到达生产井；该方法能够模拟不同润湿状态和不同裂缝位置下的毛管力和重力情况。
参考文献：
［1］王敬，刘慧卿，夏静，等.裂缝性油藏渗吸采油机理数值模拟［J］.石油勘探与开发，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］孙致学，姜传胤，张凯，等.基于离散裂缝模型的CO2增强型地热系统THM耦合数值模拟［J］.中国石油大学学报（自然科学版），2020，44（6）：79-87.
SUN Zhixue， JIANG Chuanyu， ZHANG Kai， et al. Numerical simulation for heat extraction of CO2-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］徐榮利，郭天魁，曲占庆，等.基于离散裂缝模型的页岩油储层压裂渗吸数值模拟［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］糜利栋，姜汉桥，李涛，等.基于离散裂缝模型的页岩气动态特征分析［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］张有天，刘中.降雨过程裂隙网络饱和-非饱和、非恒定渗流分析［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］苏皓，雷征东，李俊超，等.储集层多尺度裂缝高效数值模拟模型［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］姚军，张娜，黄朝琴，等.非均质油藏多尺度混合有限元数值模拟方法［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.
（编辑 李志芬）
