CN111814364A - 一种岩溶储层演化数值模拟方法 - Google Patents

Abstract

本发明公开了一种岩溶演化数值模拟方法，包括以下步骤：步骤1、建立裂缝与基质耦合地质模型；步骤2、建立多类型储层流‑固‑化耦合数学模型；步骤3、形成快速、稳定的全耦数值求解算法。利用嵌入式离散裂缝方法并结合裂缝力学本构模型，处理裂缝或断层导流和形变过程，不仅能提高计算效率，还能准确反映裂缝溶蚀和力学形变对岩溶过程的控制作用。本发明能模拟复杂岩溶环境下储层孔隙度和渗透率的演化规律，为揭示碳酸盐岩缝洞型储层的形成机制提供理论基础。

Description

一种岩溶储层演化数值模拟方法

技术领域

本发明涉及油气田开发地质领域，具体涉及一种岩溶储层演化数值模拟方法。

背景技术

碳酸盐岩缝洞型油藏的开发方案的制定，很大程度上依赖于储层内裂缝和溶洞的空间尺寸、空间形态、连通程度和组合关系。因此，弄清碳酸盐岩储集体的缝洞组合形态、认识岩溶演化规律，是实现碳酸盐岩缝洞型油藏高效开发的关键技术之一。通过岩溶储层形成过程的演化模拟，反演碳酸盐岩储集体缝洞组合特征与规律，是描述储层的重要手段与途径。

目前，与碳酸盐岩岩溶储层演化相关的数值模拟模型主要有两类，一是针对溶洞形成过程的管流模型，二是针对酸化压裂过程的双尺度蚓孔模型。

1999年Kaufmann G.，Braun J.在Water Resources Research期刊上发表了题为Karst Aquifer Evolution in Fractured Rocks的论文，建立了溶洞形成过程的管流模型。2009年Kaufmann G.又发表了题为Modelling Karst Geomorphology on DifferentTime Scales的论文，对模型进行了完善。该模型将不同形状和大小的溶蚀洞等效为具有不同半径、相互连接的圆柱形导管，采用哈根-泊肃叶定律计算导管中流体流动，通过求解带化学反应的溶质运移-扩散方程，得到流体流动速度、溶质浓度和化学反应速率。然后，建立导管半径与流体流速和反应速率的关系，确定当前流速和反应速率情况下导管形状和尺寸变化及延展规律，进而模拟储层在岩溶作用下溶洞的演化过程。该方法未考虑基质中的流动，也未考虑裂缝在整个岩溶过程中的控制作用。此外还忽略了地应力对基质和裂缝，以及溶洞中流体流动的影响。

2005年Panga M.等在AIChE Journal上发表题为Two-scale continuum modelfor simulation of wormholes in carbonate acidization的论文，建立了碳酸盐岩酸化的连续介质模型。2007年Kalia N.和Balakotaiah V.在Chemical Engineering Science上发表题为Modeling and Analysis of Wormhole Formation in Reactive Dissolutionof Carbonate Rocks的论文，建立了碳酸盐岩酸化虫洞形成的模型。该类模型通过模拟酸化虫洞，研究其对近井地带渗透性的改善，以达到提高单井产量的目的。首先，将基质、溶孔、裂缝和溶洞统一处理为具有一定孔隙度和渗透率的单一基质的一类方法，通过求解带化学反应的溶质运移-扩散方程，得到流体流动速度、溶质浓度和化学反应速率。然后，建立孔隙度和渗透率对流速、溶质浓度和化学反应速率的数学依赖关系，模拟储层在高压、强腐蚀的压裂液作用下，孔隙度和渗透率的演化过程。目前，此方法模拟酸化过程，没有考虑应力对基质、裂缝孔隙度和渗透率的影响，也未能考虑裂缝对化学反应前缘运移速度的控制作用。酸化过程的时间尺度和空间尺度较小，只考虑百米范围内的近井地带，在以小时为计量单位的储层孔隙度和渗透率变化过程，不能用于数公里甚至更大范围内，在百年甚至更长时间内储层中溶孔、溶缝和溶洞的演化规律。

碳酸盐岩储层中溶缝、溶孔、溶洞及其复合储集体形成过程极其复杂，目前国内外尚未形成考虑缝洞型储集体岩溶演化过程的数学模型。

发明内容

为了模拟复杂地质环境下溶蚀孔洞储层的演化形成过程，实现缝洞型储集体形成过程的动态表征，揭示碳酸盐岩缝洞型储层的形成机制。本发明提出了一种基于嵌入式离散裂缝的流-固-化全耦合的岩溶演化数值模拟方法，综合考虑了流体流动、溶质运移与扩散、化学反应以及岩体变形过程，实现了对缝洞型储集体形成过程的动态表征。

本发明公开一种岩溶储层演化数值模拟方法，包括以下步骤：

步骤1、建立裂缝与基质耦合地质模型；

步骤2、建立多类型储层流-固-化耦合数学模型；

步骤3、形成快速、稳定的全耦数值求解算法。

在一个实施方式中，步骤1包括以下子步骤：

步骤1.1、根据模拟区域几何边界数据和各天然裂缝的顶点数据，对基质和裂缝进行网格剖分；

步骤1.2、建立基质与裂缝、裂缝与裂缝之间的几何连接关系；

步骤1.3、计算裂缝与裂缝的交线长度、裂缝与基质的相交面积以及相应的等效距离，并结合裂缝和基质的渗透率计算传导率系数；

步骤1.4、建立联通表，包括基质与基质、基质与裂缝及裂缝与裂缝三种连接关系，得到裂缝与基质耦合地质模型。

在一个实施方式中，所述步骤1.1中采用笛卡尔正交网格或角点网格进行网格剖分。

在一个实施方式中，步骤2包括以下子步骤：

步骤2.1、根据岩溶演化的主要物理、化学过程及其耦合关系，推导得到流-固-化耦合系统的控制方程；

步骤2.2、根据步骤2.1中的所述控制方程，结合不同类型储层的化学反应动力学方程，建立多类型储层流-固-化耦合数学模型。

在一个实施方式中，所述流-固-化耦合系统的控制方程包括：式(1)所示的流体质量守恒方程，式(2)所示的带有化学反应的溶质对流-扩散方程，和式(3)所示的力学平衡方程，

其中，式(2)中的孔隙度由式(4)演化得到，

式中：p为流体压力，ε为应变张量，D为四阶弹性矩阵张量，δ为克罗内克张量，ρ_f和ρ_b分别为流体密度和岩石骨架与孔隙内流体体积平均密度，

为达西流速，C_f和C_s分别为孔隙内的钙离子浓度和固-液交界面上的钙离子浓度，D_e为弥散系数张量，k_c为局部质量交换系数，a_v为单位体积中用于化学反应的表面积，ε_V为体应变，g重力加速度，b为比奥系数，K_s为基质颗粒体积模量，μ为粘度，k为渗透率，B为流体体积系数，φ和φ₀分为当前孔隙度和初始孔隙度，t为模拟时间，q_f为质量源/汇项。

在一个实施方式中，步骤3包括以下子步骤：

步骤3.1、对步骤2中所述的多类型储层流-固-化耦合数学模型进行数值离散；

步骤3.2、基于步骤1.4中所述的联通表构建雅克比矩阵结构，形成非线性/线性矩阵求解格式；

步骤3.3、采用有限体积方法和有限元的方法对多类型储层流-固-化耦合数学模型进行离散求解。

在一个实施方式中，所述有限体积方法用于求解流体压力和溶质浓度，所述有限元方法用于求解应力、应变和位移。

在一个实施方式中，其中步骤3.1中所述数值离散得到式(19)：

式中：R^u、R^p和R^C分别为所述力学平衡方程、所述流体质量守恒方程和所述带有化学反应的溶质对流-扩散方程的残差，K为刚度矩阵，Δt为时间步长，L和M为多孔弹性耦合矩阵，

与流体压力有关的压缩矩阵，

为与流体应力有关的压缩矩阵，F为与流体压力有关的传导系数矩阵，Q为与流体浓度有关的传导系数矩阵，上标u、p、C、n和k分别表示位移、流体压力、浓度、时间步和迭代步，下标i/j表示单元，下标a/b表示节点，Δu、Δp、ΔC分别表示牛顿迭代步内位移、流体压力、钙离子浓度的变化，Cⁿ表示时间步n初始时刻的钙离子浓度。

在一个实施方式中，其中步骤3.3中对所述数值离散得到的式(19)采用Newton-Raphson方法进行求解，得到式(29)：

式中：X为待求解的变量系统，δ为变量增量符号，R为某一牛顿迭代步内的残差，P和S分别表示主变量和次要变量，n和k分别表示时间步和迭代步。

在一个实施方式中，所述采用Newton-Raphson方法进行求解包括以下步骤：

首先，采用Newton-Raphson方法统一求解式(1)、式(2)和式(3)，得到流体压力、位移、应力、应变和溶质浓度；

然后，根据式(4)计算孔隙度，并代入Carman-Kozeny方程求解孔喉半径、渗透率和单位体积中用于化学反应的表面积；

最后，通过孔隙结构与物理化学参数之间的关系，求得弥散系数、流速和密度。

上述技术特征可以各种技术上可行的方式组合以产生新的实施方案，只要能够实现本发明的目的。

附图说明

在下文中将基于实施例并参考附图来对本发明进行更详细的描述。其中：

图1是本发明实施例中的基于嵌入式离散裂缝的流-固-化耦合岩溶演化数值模拟方法流程示意图；

图2是本发明实施例中的流-固-化耦合系统原理图；

图3是本发明实施例中的有限体积有限元方法的插值函数；

图4是本发明实施例中岩溶演化模型实现子模块工程；

图5是本发明实施例中的岩溶演化构造裂缝模型；

图6是本发明实施例中的岩溶演化构造裂缝模型模拟得到了100年后储层孔隙度的分布及其变化情况；

图7是本发明实施例中的岩溶演化构造裂缝模型模拟得到了200年后储层孔隙度的分布及其变化情况；

图8是本发明实施例中的岩溶演化构造裂缝模型模拟得到了300年后储层孔隙度的分布及其变化情况。

在附图中，相同的部件使用相同的附图标记。附图并未按照实际的比例。

具体实施方式

下面将结合附图对本发明作进一步说明。

如图1所式为本发明是实例中的岩溶演化数值模拟方法流程示意图，下面结合对本发明中的岩溶演化数值模拟方法做一个具体的介绍。

本发明提出的是一种基于嵌入式离散裂缝的流-固-化全耦合的岩溶演化数值模拟方法，主要包括三步：步骤1、建立裂缝与基质耦合地质模型；步骤2、建立多类型储层流-固-化耦合数学模型；步骤3、形成快速、稳定的全耦数值求解算法。

下面分别对三个步骤进行具体的说明。

步骤1、建立裂缝与基质耦合地质模型。

裂缝系统与基质系统作为两个独立的系统耦合处理，裂缝系统采用显式离散方式处理，而不是传统双重介质方法所采用的等效方式处理。该方法的优势在于不需要对裂缝进行等效隠式处理，而是保持天然裂缝原有几何形态显示处理，并且具有与经过局部加密处理的离散裂缝方法(DFN)十分接近的计算精度。无需对裂缝所在区域进行加密处理，可以在简单的结构化网格中模拟裂缝高导流通道的效果，具有很高的计算效率。因此，基于嵌入式离散裂缝方法(EDFM)建立裂缝与基质耦合地质模型的优势可概括为：1)裂缝建模方便，研究区域可采用简单的结构化网格剖分，给出每条裂缝的四个顶点坐标，所实现的前处理程序能自动生成联通表用于后续计算；2)无需在裂缝附近区域进行网格加密处理，计算效率高；3)计算精度高，能准确模拟化学反应前缘。

该离散裂缝建模方法的具体实施步骤，即步骤1包括的子步骤如下：

步骤1.1、根据模拟区域几何边界数据和各天然裂缝的顶点数据，对基质和裂缝进行网格剖分(其中顶点可不在模拟区域范围内，非平面天然裂缝可以通过多个平面裂缝拼接)；

步骤1.2、建立基质与裂缝、裂缝与裂缝之间的几何连接关系(即建立非相邻联通表)；

步骤1.3、计算裂缝与裂缝的交线长度、裂缝与基质的相交面积以及相应的等效距离，并结合裂缝和基质的渗透率计算传导率系数；

步骤1.4、建立联通表，包括基质与基质、基质与裂缝及裂缝与裂缝三种连接关系，得到裂缝与基质耦合地质模型。

在一个优选的实施例中，步骤1.1中采用笛卡尔正交网格或角点网格进行网格剖分。

步骤2、建立多类型储层流-固-化耦合数学模型。

裂缝和基质耦合地质模型的建立是流-固-化耦合模拟方法的基础，接下来需要推导流-固-化耦合系统的控制方程(即地质模型背后的数学模型)。首先，根据岩溶演化的主要物理、化学过程及其耦合关系，推导得到流-固-化耦合系统的控制方程；然后，结合不同类型储层的化学反应动力学方程，建立多类型储层流-固-化耦合数学模型。

岩溶演化涉及化学反应、岩石变形、溶质运移、流体对流和扩散等过程，每个过程都可以基于一定的假设和简化，概化出相应的物理数学模型。模型概化的基本原则是：假设合理，简化得当，能准确反映真实的物理过程。图2所示为流-固-化耦合系统原理图，盐岩体的化学反应改变岩石的孔隙结构和比表面积，从而使孔隙度、渗透率、反应速率发生相应变化，并最终对流体流动、溶质运移产生显著影响；流体流动伴随溶质运移；溶质运移改变局部反应物浓度，影响反应速率。图2中展示的是岩溶演化包含的各个过程及各过程间的相互耦合关系：(1)溶有CO₂的淡水及地下水与碳酸盐岩发生化学反应，形成不同的溶蚀结构类型，该过程取决于盐岩的矿物组成、化学反应速率、流动类型、储层非均质性、流体滤失进入基质的速度等；(2)流体在压力势的驱动下，在基质、裂缝及溶蚀洞内发生渗流或自由流动；(3)溶质(Ca²⁺，Mg²⁺，H⁺，CO₃ ^2-等)随流体流动产生对流，在浓度势驱动下发生扩散；(4)在构造应力、自重应力等条件下，岩体发生变形，包括孔隙体积的变化和裂缝开度的变化。岩体在应力作用下产生变形，会改变基质的孔隙度和渗透率，同时改变裂缝的开度，进而影响流体的流动和溶质运移。裂缝开度变化主要与应力大小、裂缝面上的矿物溶解、运输、沉淀等过程有关。其中，化学溶解改变岩体的孔隙结构，影响孔隙的比表面积、渗透率和孔隙度，从而改变流体的流动速度及溶质分布情况；溶质运移改变溶质浓度，进而影响化学溶解过程；流体流动引起浓度的改变，浓度改变反过来又会影响化学溶解过程。因此，整个系统通过耦合作用共同影响岩溶演化过程。

下面对系统的控制方程进行具体的介绍，所述流-固-化耦合系统的控制方程包括：式(1)所示的流体质量守恒方程，式(2)所示的带有化学反应的溶质对流-扩散方程，和式(3)所示的力学平衡方程：

其中，式(2)中的孔隙度由式(4)演化得到，式(4)表征的是溶蚀及压力和应力变化引起的孔隙度变化；

以上各式中：p为流体压力，ε为应变张量，D为四阶弹性矩阵张量，δ为克罗内克张量，ρ_f和ρ_b分别为流体密度和岩石骨架与孔隙内流体体积平均密度，

为达西流速，C_f和C_s分别为孔隙内的钙离子浓度和固-液交界面上的钙离子浓度，D_e为弥散系数张量，k_c为局部质量交换系数，a_v为单位体积中用于化学反应的表面积，ε_V为体应变，g重力加速度，b为比奥系数，K_s为基质颗粒体积模量，μ为粘度，k为渗透率，B为流体体积系数，φ和φ₀分为当前孔隙度和初始孔隙度，t为模拟时间，q_f为质量源/汇项。

除了以上的控制方程外，还包括一些其他的辅助方程，例如应变-位移方程、应力-应变方程、裂缝本构方程、达西渗流方程等。

溶质在孔隙网络中的质量传输由以下方程确定：

式中，S_h为Sherwood数或无量纲质量传输系数，k_c为局部质量交换系数，r_P为平均孔隙半径，D_m为酸在任意剪切速度下有效分子扩展系数，

为渐进Sherwood数，m_p为孔隙长度和半径的比值，

为孔隙尺度的雷诺数，为惯性力与粘性力的比值，υ为动力粘度，

为Schmidt数，α_os为依赖于孔隙连通性的常系数，D_eX为纵向(X方向)弥散系数，D_eT为横向(Y或Z方向)弥散系数，λ_x和λ_T为依赖于孔隙结构的常系数，φ为当前孔隙度。

达西模型中岩石的物理化学参数(局部孔隙尺寸、孔隙表面积和渗透率)取决于孔隙结构，可通过建立孔隙结构与物理化学参数之间的关系(Carman-Kozeny方程)来考虑物理化学参数的变化，具体如式(8)和式(9)所示：

式中：I为二阶单位张量，k₀为初始平均值渗透率，φ₀为初始孔隙度，r₀为初始平均孔隙半径，a₀为初始比表面积，β为扩孔系数，k、φ、r_P、a_v分别为当前的渗透率、当前孔隙度、平均孔隙半径和单位体积中用于化学反应的表面积。

步骤3、形成快速、稳定的全耦数值求解算法。

在步骤2中得到多类型储层流-固-化耦合数学模型以后，需要对该数学模型进行求解。步骤3可以包括以下子步骤：

步骤3.1、对步骤2中所述的多类型储层流-固-化耦合数学模型进行数值离散；

步骤3.2、基于步骤1.4中所述的联通表构建雅克比矩阵结构，形成非线性/线性矩阵求解格式；

步骤3.3、采用有限体积和有限元的方法对多类型储层流-固-化耦合数学模型进行离散求解。

优选的，其中有限体积方法用于求解流体压力和溶质浓度，有限元方法用于求解应力、应变和位移。

如图3所示，为有限体积有限元方法的插值函数。包括孔隙压力、浓度和位移的插值函数，其中，位移变量位于块节点，孔隙压力和浓度位于块中心。在图3中，其中(a)为单元内位移插值函数数值分布，(b)为单元内孔隙压力和浓度插值函数数值分布；(c)为主变量在单元内对应的位置示意图。

所述插值函数的具体表达式如式(10)所示：

其中，位移插值函数ψ_i满足在节点i处值为1，在其他所有节点处为0，在单元内由式(11)所示的六面体单元有限元插值函数插值得到：

孔隙压力和浓度的插值函数

和η_l满足在单元a和l处值为1，其他单元为0，该插值函数在单元内处处为1或0，与位移插值函数存在区别，具体如图3所示。将位移插值函数代入式(3)所示的控制方程，并整理成残差迭代格式，则式(3)所示的力学平衡方程可整理为如式(14)所示，

其中，

将孔隙压力和浓度的插值函数分别代入式(1)和式(2)所示的控制方程，并整理成残差迭代格式，则式(1)所示的流体质量守恒方程可整理为如式(17)所示：

式(2)所示的带有化学反应的溶质对流-扩散方程可整理为如式(18)所示：

进一步将式(14)、式(17)和式(18)整理得到如式(19)所示的整个系统的牛顿迭代矩阵形式：

式中：R^u、R^p和R^C分别为所述力学平衡方程、所述流体质量守恒方程和所述带有化学反应的溶质对流-扩散方程的残差，K为刚度矩阵，Δt为时间步长，L和M为多孔弹性耦合矩阵，

与流体压力有关的压缩矩阵，

为与流体应力有关的压缩矩阵，F为与流体压力有关的传导系数矩阵，Q为与流体浓度有关的传导系数矩阵，上标u、p、C、n和k分别表示位移、流体压力、浓度、时间步和迭代步，下标i/j表示单元，下标a/b表示节点，Δu、Δp、ΔC分别表示牛顿迭代步内位移、流体压力、钙离子浓度的变化，Cⁿ表示时间步n初始时刻的钙离子浓度。

式(19)中的各个矩阵块如式(20)至式(28)所示：

对多场耦合的岩溶过程，本发明拟采用全耦合的求解策略实施求解(如图1所示)，即在同一个时间步内，同时求解主变量压力、位移和浓度，其他变量(包括孔隙度、渗透率、弥散数、密度和流速)均可表示为主变量的函数，通过函数关系可以求得次要变量的值，并做更新，进入下一个时间步的计算。整个系统采用Newton-Raphson方法对偏微分方程进行线性化求解：

式中，X为待求解的变量系统，δ为变量增量符号，R为某一牛顿迭代步内的残差，P和S分别表示主变量和次要变量，n和k分别表示时间步和迭代步。

该求解系统以流体压力、位移和溶质浓度为主要变量，孔隙度、渗透率、比表面积、弥散系数、流速、密度等为次要变量。

优选的，采用Newton-Raphson方法包括以下步骤：

首先，采用Newton-Raphson方法统一求解式(1)、式(2)和式(3)得到流体压力、位移、应力、应变和溶质浓度；

然后，根据式(4)计算孔隙度，并代入Carman-Kozeny方程求解孔喉半径、渗透率和单位体积中用于化学反应的表面积；

最后，通过孔隙结构与物理化学参数之间的关系，求得弥散系数、流速和密度。

根据岩溶演化数学模型特点，设计相应的C++程序设计框架，按照模块化设计理念，各个功能模块既相互独立，又遵循代码共享的特点，达到代码重复利用率高、模块可扩展性强的设计目标。如图4所示，整个岩溶演化模拟系统KarstEvolSys按功能分为三个部分：数据处理部分、模拟器驱动部分和模型耦合模块部分。其中，数据处理部分由三个子模块组成，包括数据文件输入输出管理的FileIO模块、负责数据存储与查询功能的Database模块以及负责输出一维曲线、二维曲面、三维等值面的3D_Visualization模块；模拟器驱动模块包括线性求解器模块LinearSolver(本模型采用Intel数学核心函数库MKL)、处理弹性变形和弹塑性变形的本构模型模块ConstituModel、根据网格自动生成两点或多点流动格式联通表的模块TP-MPFA模块以及计算化学性质的模块ChemProps；第三部分为模型耦合模块，包括流体流动模块FluidFlow、岩石变形模块GeoMech以及化学反应模块ChemReac，三个模块通过耦合矩阵联立为一个整体求解系统，从而能在同一时间步将主变量流体压力p、位移u和溶质浓度C_f，同时在每个迭代步将次要变量渗透率k、孔隙度φ、单位体积中用于化学反应的表面积a_V、弥散系数S_h等完成更新。

本发明提出了一种流-固-化全耦合的岩溶储层演化数值模拟方法，该方法通过嵌入式离散裂缝模型显示模拟裂缝对岩溶演化的控制作用，同时考虑应力对岩体渗透率和孔隙度变化的影响。采用混合有限体积有限元方法(mixed FV/FE)离散整个系统，能够模拟复杂环境下裂缝和溶洞的演化过程，能够为碳酸盐缝洞型储层储集体建模提供理论依据。其方法可以通过软件编程实现。软件系统流程可描述如下：

(1)工程设置：设置模拟岩溶演化过程需要的文件、参数；

(2)计算实例设置：设置流量、流体压力、浓度和应力等初始和边界条件；

(3)根据构造裂缝几何参数建立裂缝与基质的耦合地质模型；

(4)建立多类型储层岩溶演化数学模型；

(5)根据流-固-化耦合方法初始化流体流动、溶质浓度和原位地应力；

(6)采用全耦合求解算法求解流-固-化耦合的岩溶演化系统；

(7)根据需要输出需要的参数，如孔隙度、渗透率等；

(8)采用Tecplot绘图软件对模拟结果进行处理和结果展示。

使用本发明提出的方法模拟含天然构造裂缝的储层演化过程，主要的模拟流程是：①建立数模网格模型，对三维地质模型的网格数据进行粗化；②根据天然构造裂缝采用前面介绍的嵌入式离散裂缝方法实施预处理，生成非相邻联通表，构建雅克比矩阵结构；③准备输入文本文件，并对流-固-化耦合模型进行初始化；④运行模拟器得到模拟结果，并对结果进行处理。如图5所示，整个模拟区域范围为500m×500m×10m(沿垂直纸面的水平方向为10m)，基质孔隙度为5％，渗透率为5mD，模拟区域包含一组共轭天然构造缝，裂缝渗透率均为500mD，年降雨量为20mm，地表温度为25℃，地表大气压为1个标准大气压(即0.1MPa)，流体压力和温度分布分别按照静水压力和地温梯度(按照25℃/km)分别计算得到，二氧化碳分压为0.05个标准大气压(即0.005MPa)，二氧化碳平衡浓度为2.1mol/m³，然后按照二氧化碳浓度平衡理论计算二氧化碳浓度分布，最终模拟得到了100年(图6)、200年(图7)和300年(图8)后。模拟结果显示了不同时间储层孔隙度的分布及其变化情况。从结果可知本发明可模拟天然裂缝性储层随时间的岩溶演化过程，得到了很好的展示效果，能够为碳酸盐岩缝洞型油藏高效开发方案的制定提供指导。

虽然已经参考优选实施例对本发明进行了描述，但在不脱离本发明的范围的情况下，可以对其进行各种改进并且可以用等效物替换其中的部件。尤其是，只要不存在结构冲突，各个实施例中所提到的各项技术特征均可以任意方式组合起来。本发明并不局限于文中公开的特定实施例，而是包括落入权利要求的范围内的所有技术方案。

Claims (10)

1.一种岩溶储层演化数值模拟方法，包括以下步骤：

步骤1、建立裂缝与基质耦合地质模型；

步骤2、建立多类型储层流-固-化耦合数学模型；

步骤3、形成快速、稳定的全耦数值求解算法。

2.根据权利要求1所述的岩溶储层演化数值模拟方法，其特征在于，步骤1包括以下子步骤：

步骤1.1、根据模拟区域几何边界数据和各天然裂缝的顶点数据，对基质和裂缝进行网格剖分；

步骤1.2、建立基质与裂缝、裂缝与裂缝之间的几何连接关系；

步骤1.3、计算裂缝与裂缝的交线长度、裂缝与基质的相交面积以及相应的等效距离，并结合裂缝和基质的渗透率计算传导率系数；

步骤1.4、建立联通表，包括基质与基质、基质与裂缝及裂缝与裂缝三种连接关系，得到裂缝与基质耦合地质模型。

3.根据权利要求2所述的岩溶储层演化数值模拟方法，其特征在于，所述步骤1.1中采用笛卡尔正交网格或角点网格进行网格剖分。

4.根据权利要求1所述的岩溶储层演化数值模拟方法，其特征在于，步骤2包括以下子步骤：

步骤2.1、根据岩溶演化的主要物理、化学过程及其耦合关系，推导得到流-固-化耦合系统的控制方程；

步骤2.2、根据步骤2.1中的所述控制方程，结合不同类型储层的化学反应动力学方程，建立多类型储层流-固-化耦合数学模型。

5.根据权利要求4所述的岩溶储层演化数值模拟方法，其特征在于，所述流-固-化耦合系统的控制方程包括：式(1)所示的流体质量守恒方程，式(2)所示的带有化学反应的溶质对流-扩散方程，和式(3)所示的力学平衡方程，

其中，式(2)中的孔隙度由式(4)演化得到，

式中：p为流体压力，ε为应变张量，D为四阶弹性矩阵张量，δ为克罗内克张量，ρ_f和ρ_b分别为流体密度和岩石骨架与孔隙内流体体积平均密度，

为达西流速，C_f和C_s分别为孔隙内的钙离子浓度和固-液交界面上的钙离子浓度，D_e为弥散系数张量，k_c为局部质量交换系数，a_v为单位体积中用于化学反应的表面积，ε_V为体应变，g重力加速度，b为比奥系数，K_s为基质颗粒体积模量，μ为粘度，k为渗透率，B为流体体积系数，φ和φ₀分为当前孔隙度和初始孔隙度，t为模拟时间，q_f为质量源/汇项。

6.根据权利要求5所述的岩溶储层演化数值模拟方法，其特征在于，步骤3包括以下子步骤：

步骤3.1、对步骤2中所述的多类型储层流-固-化耦合数学模型进行数值离散；

步骤3.2、基于步骤1.4中所述的联通表构建雅克比矩阵结构，形成非线性/线性矩阵求解格式；

步骤3.3、采用有限体积和有限元的方法对多类型储层流-固-化耦合数学模型进行离散求解。

7.根据权利要求6所述的岩溶储层演化数值模拟方法，其特征在于，所述有限体积方法用于求解流体压力和溶质浓度，所述有限元方法用于求解应力、应变和位移。

8.根据权利要求6所述的岩溶储层演化数值模拟方法，其特征在于，其中步骤3.1中所述数值离散得到式(19)：

式中：R^u、R^p和R^C分别为所述力学平衡方程、所述流体质量守恒方程和所述带有化学反应的溶质对流-扩散方程的残差，K为刚度矩阵，Δt为时间步长，L和M为多孔弹性耦合矩阵，

与流体压力有关的压缩矩阵，

为与流体应力有关的压缩矩阵，F为与流体压力有关的传导系数矩阵，Q为与流体浓度有关的传导系数矩阵，上标u、p、C、n和k分别表示位移、流体压力、浓度、时间步和迭代步，下标i/j表示单元，下标a/b表示节点，Δu、Δp、ΔC分别表示牛顿迭代步内位移、流体压力、钙离子浓度的变化，Cⁿ表示时间步n初始时刻的钙离子浓度。

9.根据权利要求8所述的岩溶储层演化数值模拟方法，其特征在于，其中步骤3.3中对所述数值离散得到的式(19)采用Newton-Raphson方法进行求解，得到式(29)：

式中：X为待求解的变量系统，δ为变量增量符号，R为某一牛顿迭代步内的残差，P和S分别表示主变量和次要变量，n和k分别表示时间步和迭代步。

10.根据权利要求9所述的岩溶储层演化数值模拟方法，其特征在于，所述采用Newton-Raphson方法进行求解包括以下步骤：

首先，采用Newton-Raphson方法统一求解式(1)、式(2)和式(3)得到流体压力、位移、应力、应变和溶质浓度；

然后，根据式(4)计算孔隙度，并代入Carman-Kozeny方程求解孔喉半径、渗透率和单位体积中用于化学反应的表面积；

最后，通过孔隙结构与物理化学参数之间的关系，求得弥散系数、流速和密度。

Images