缝洞组合模型的构建方法及探测储层缝洞变形的方法

摘要

本发明属于油气田开发技术领域，具体涉及一种缝洞组合模型的构建方法及探测储层缝洞变形的方法。本发明提出的缝洞组合模型的构建方法不仅考虑了溶洞整体控制模型，而且考虑了裂缝整体控制模型，同时引入损伤参数对控制模型进行修正，使得缝洞组合模型更加接近于实际待分析的储层缝洞，使用本发明的缝洞组合模型对实际的储层缝洞进行分析，考虑了应力特征和流体渗流情况下的缝洞组合变形，可以根据油藏应力状况结合渗流特征，对油藏开发过程中引起的缝洞组合变形特征进行预测，为开发调整提供依据。

权利要求

1.一种缝洞组合模型的构建方法，其特征在于，包括：获得溶洞整体控制模型；

获得裂缝整体控制模型；

对所述溶洞整体控制模型和所述裂缝整体控制模型进行耦合叠加，获得缝洞组合模型；

其中，所述溶洞整体控制模型为：

其中，K表示渗透系数张量， 表示材料损伤时节点位移与节点力的整体关系矩阵，[Kup]表示孔隙压力与节点力的整体关系矩阵，[Kpp]表示孔隙压力与渗流节点载荷的整体关系矩阵，Δp表示孔隙压力增量，ΔU表示位移增量，{ΔRU}表示节点力增量，{ΔRp}表示渗流节点载荷增量；

获得所述裂缝整体控制模型的方法包括：获得裂缝流体流动数学模型；获得裂缝体积变形控制模型；获得孔隙介质的有效应力模型；获得裂缝渗透系数控制模型；根据所述裂缝流体流动数学模型、所述裂缝体积变形控制模型、所述孔隙介质的有效应力模型和所述裂缝渗透系数控制模型获得初始裂缝整体控制模型；还包括修正步骤，获得损伤参数，利用所述损伤参数对初始裂缝整体控制模型进行修正，得到修正后的裂缝整体控制模型；

所述修正后的裂缝整体控制模型为：

其中，K表示渗透系数张量， 表示孔隙压力与节点力的整体关系矩阵， 表示孔隙压力与渗流节点载荷的整体关系矩阵， 表示节点位移与节点力的整体关系矩阵，ΔPf表示孔隙压力增量，ΔUf表示位移增量， 表示节点力增量， 表示渗流节点载荷增量。

2.根据权利要求1所述的方法，其特征在于，进行所述耦合叠加的条件为：|p-pf|≤5％*pt和|v-vf|≤1％*v，或|p-pf|≤1％*p和|v-vf|≤1％*v，其中，p为所述溶洞整体控制模型在每个时间步得到的溶洞周边孔隙压力，pf为所述裂缝整体控制模型在与所述溶洞整体控制模型相同时间步得到的裂缝周边孔隙压力，pt为每个时间步的实际压力值，所述实际压力值通过查询开发过程中生产井的测压数据获得，v为所述溶洞整体控制模型在每个时间步得到的溶洞内流体流速，vf为所述裂缝整体控制模型在与所述溶洞整体控制模型相同时间步得到的裂缝内流体流速。

3.根据权利要求1所述的方法，其特征在于，所述裂缝流体流动数学模型为：sφevi＝kijJj，

其中，s为流体的饱和度，φe为岩石的有效孔隙率，vi为i向量方向的流体速度，Jj为j向量方向的测压水头的梯度，kij为渗透张量，记作[K]：

4.根据权利要求3所述的方法，其特征在于，获得所述裂缝体积变形控制模型包括以下步骤：获得外力压力作用下的裂缝体积第一应变函数；

获得裂缝整体变形的第二应变函数；

将所述第一应变函数和所述第二应变函数联立，获得所述裂缝体积变形控制模型。

5.根据权利要求4所述的方法，其特征在于，所述第一应变函数为：其中，σp为外部压力，Kv为多孔介质的体积模量，εv为体积应变量。

6.根据权利要求5所述的方法，其特征在于，所述第二应变函数为：其中，φ为岩石的孔隙度，Ks为固体颗粒的体积模量，Kf为孔隙流体的体积模量。

7.根据权利要求6所述的方法，其特征在于，所述裂缝体积变形控制模型为：

8.根据权利要求7所述的方法，其特征在于，所述孔隙介质的有效应力模型为：其中，σij为总应力，δij为Kronecker符号，σij'为有效应力， 为流体的平均压力，α为Biot系数，

9.根据权利要求8所述的方法，其特征在于，所述裂缝渗透系数控制模型为：k＝k0exp(-βσ')，其中k为渗透系数，σ'＝γh-pf，γh为上覆岩层的重力，pf为裂缝周边孔隙压力，β为常数，k0为上覆岩层的重力与裂缝周边孔隙压力相等时的渗透系数。

10.根据权利要求9所述的方法，其特征在于，所述损伤参数的表达式为：根据如下表达式计算存在损伤时单元的弹性模量：

Eξ＝E0(1-ξ)

其中，ξ表示损伤参数， 表示归一化的等效塑性应变，a表示材料常数，E0和Eξ分别表示不存在损伤时和存在损伤时单元的弹性模量。

11.一种探测储层缝洞变形的方法，其特征在于，所述方法包括：构建待分析的储层缝洞的模型，对所述储层缝洞的模型进行网格划分；

获取所述储层缝洞的边界条件；

获取所述储层缝洞的地层岩石材料的力学参数；

获取所述储层缝洞的应力状况；

根据所述边界条件、所述力学参数和所述应力状况，建立所述储层缝洞的组合模型；

建立溶洞垮塌模拟与裂缝闭合模拟的耦合条件，对所述储层缝洞的组合模型的变形状况进行模拟分析，其中，所述储层缝洞的组合模型是采用如权利要求1至10中任意一项所述的方法建立的。

12.根据权利要求11所述的方法，其特征在于，所述边界条件是根据研究区构造演化史、区域应力变化特征和生产动态状况获取的。

13.根据权利要求11所述的方法，其特征在于，所述力学参数是通过测井曲线计算或岩石物理实验获取的。

14.根据权利要求11所述的方法，其特征在于，所述应力状况是通过临近井测试获取或通过高精度地应力场获取。

15.根据权利要求14所述的方法，其特征在于，所述地应力场的模型的获得方法包括：根据研究区区域地质调查结果，建立地质模型；

结合岩石变形的力学特性，给定相应的岩石变形物理方程及对应的力学参数和边界条件，将所述地质模型转化为力学模型；

通过不断调整边界作用力方式和大小、方向，利用有限元方法对所述力学模型进行模拟运算，当模拟结果与已有地应力实测值、地震震源机制或井孔地应力达到最佳拟合时，所述力学模型为所述地应力场的模型。

说明书

缝洞组合模型的构建方法及探测储层缝洞变形的方法

技术领域

[0001] 本发明属于油气田开发技术领域，具体涉及一种缝洞组合模型的构建方法及探测储层缝洞变形的方法。

背景技术

[0002] 缝洞型碳酸盐岩油藏主要的储集空间为各种不同形态的大型溶洞，裂缝主要起导流作用，但孤立的溶洞、裂缝均较少，储集体主要以不同类型的缝洞组合形式存在。在油藏开采过程中随着地层压力的下降，溶洞及裂缝都会发生不同程度的变形，导致溶洞产生垮塌效应，裂缝会出现闭合情况，影响了油井的产量及最终采收率。因此，对缝洞组合在油藏条件下的变形规律分析具有重要意义。

[0003] 目前，相关研究主要包括：ZL201020259505.4提出了一种岩溶洞穴储层埋藏保存程度的预测装置，可以在无钻井区钻前预测岩溶洞穴储层埋藏保存程度；郑兴平等人(郑兴平，沈安江，寿建峰和潘文庆，埋藏岩溶洞穴垮塌深度定量图版及其在碳酸盐岩缝洞型储层地质评价预测中的意义，海相油气地质，14(4)：55-59,2009)通过岩溶洞穴垮塌的岩石力学原理及塔里木盆地奥陶系实钻资料的分析，明确了洞穴埋藏垮塌的控制因素，并得到了洞穴垮塌深度的定量图版；201510315477.3公开了一种用于确定缝洞型油藏中生产井钻遇溶洞发生坍塌事件的方法及系统；王超等人(王超，张强勇，刘中春，张岳和李兴军，缝洞型油藏裂缝宽度变化预测模型及其应用，中国石油大学学报(自然科学版)，40(1)：86-91,2016)，建立考虑内压变化的有限元力学模型，获得不同形态裂缝闭合过程及闭合规律；王超等人(王超，张强勇，刘中春，袁圣渤，张绪涛和王有法，缝洞型油藏溶洞临界垮塌深度预测模型及其应用，中国石油大学学报(自然科学版)，39(1)：103-110,2015)提出确定溶洞临界垮塌埋深的二分深度折减法，分析了溶洞垮塌对多因素的敏感性。

[0004] 现有技术主要集中于利用概念模型对单一溶洞的垮塌或裂缝的闭合进行预测，没有考虑流体渗流作用溶洞或裂缝变形的影响，且不能对缝洞组合的变形进行预测。

发明内容

[0005] 针对现有技术中的不足，本发明提出了缝洞组合模型的构建方法及探测储层缝洞变形的方法，本发明综合考虑应力特征与流体渗流情况，并可以根据油藏应力状况结合流体渗流特征，对油藏开发过程中引起的缝洞组合变形特征进行预测，为开发调整提供依据。

[0006] 本发明提出的缝洞组合模型的构建方法，包括：获得溶洞整体控制模型；获得裂缝整体控制模型；对所述溶洞整体控制模型和所述裂缝整体控制模型进行耦合叠加，获得缝洞组合模型。

[0007] 作为对该方法的进一步改进，进行所述耦合叠加的条件为：

[0008] |p-pf|≤5％*pt和|v-vf|≤1％*v，或

[0009] |p-pf|≤1％*p和|v-vf|≤1％*v，

[0010] 其中，p为所述溶洞整体控制模型在每个时间步得到的溶洞周边孔隙压力，pf为所述裂缝整体控制模型在与所述溶洞整体控制模型相同时间步得到的裂缝周边孔隙压力，pt为每个时间步的实际压力值，所述实际压力值通过查询开发过程中生产井的测压数据获得，v为所述溶洞整体控制模型在每个时间步得到的溶洞内流体流速，vf为所述裂缝整体控制模型在与所述溶洞整体控制模型相同时间步得到的裂缝内流体流速。

[0011] 在一个优选的实施例中，获得所述裂缝整体控制模型的方法包括：获得裂缝流体流动数学模型；获得裂缝体积变形控制模型；获得孔隙介质的有效应力模型；获得裂缝渗透系数控制模型；根据所述裂缝流体流动数学模型、所述裂缝体积变形控制模型、所述孔隙介质的有效应力模型和所述裂缝渗透系数控制模型获得初始裂缝整体控制模型。

[0012] 进一步，获得所述裂缝整体控制模型的方法还包括模型修正步骤，获得损伤参数，利用所述损伤参数对初始裂缝整体控制模型进行修正，得到修正后的裂缝整体控制模型。

[0013] 在一个优选的实施例中，所述裂缝流体流动数学模型为：

[0014] sφevi＝kijJj，

[0015] 其中，s为流体的饱和度，φe为岩石的有效孔隙率，vi为i向量方向的流体速度，Jj为j向量方向的测压水头的梯度，kij为渗透张量，记作[K]：

[0016]

[0017] 在一个优选的实施例中，获得所述裂缝体积变形控制模型包括以下步骤：获得外力压力作用下的裂缝体积第一应变函数；获得裂缝整体变形的第二应变函数；将所述第一应变函数和所述第二应变函数联立，获得所述裂缝体积变形控制模型。

[0018] 作为对裂缝流体流动数学模型的进一步改进，所述第一应变函数为：

[0019] 其中，σp为外部压力，Kv为多孔介质的体积模量，εv为体积应变量。

[0020] 作为对裂缝流体流动数学模型的进一步改进，所述第二应变函数为：

[0021] 其中，φ为岩石的孔隙度，Ks为固体颗粒的体积模量，Kf为孔隙流体的体积模量。

[0022] 作为对裂缝流体流动数学模型的进一步改进，所述裂缝体积变形控制模型为：

[0023]

[0024] 在一个优选的实施例中，所述孔隙介质的有效应力模型为：

[0025] 其中，σij为总应力，δij为Kronecker符号，σij'为有效应力， 为流体的平均压力，α为Biot系数，

[0026] 在一个优选的实施例中，所述裂缝渗透系数控制模型为：

[0027] k＝k0exp(-βσ')，其中k为渗透系数，σ'＝γh-pf，γh为上覆岩层的重力，pf为裂缝周边孔隙压力，β为常数，k0为上覆岩层重力与裂缝周边孔隙压力相等时的渗透系数。

[0028] 在一个优选的实施例中，所述损伤因子的表达式为：

[0029]

[0030] 根据如下表达式计算存在损伤时单元的弹性模量：

[0031] Eξ＝E0(1-ξ)

[0032] 其中，ξ表示损伤参数， 表示归一化的等效塑性应变，a表示材料常数，E0和Eξ分别表示不存在损伤时和存在损伤时单元的弹性模量。

[0033] 进一步，所述修正后的裂缝整体控制模型为：

[0034]

[0035] 其中，K表示渗透系数张量， 表示孔隙压力与节点力的整体关系矩阵，表示孔隙压力与渗流节点载荷的整体关系矩阵， 表示节点位移与节点力的整体关系矩阵，ΔPf表示孔隙压力增量，ΔUf表示位移增量， 表示节点力增量，表示渗流节点载荷增量。

[0036] 在一个优选的实施例中，所述溶洞整体控制模型为：

[0037]

[0038] 其中，K表示渗透系数张量， 表示材料损伤时节点位移与节点力的整体关系矩阵，[Kup]表示孔隙压力与节点力的整体关系矩阵，[Kpp]表示孔隙压力与渗流节点载荷的整体关系矩阵，Δp表示孔隙压力增量，ΔU表示位移增量，{ΔRU}表示节点力增量，{ΔRp}表示渗流节点载荷增量。

[0039] 本发明同时提出了一种探测储层缝洞变形的方法，其特征在于，所述方法包括：构建待分析的储层缝洞的模型，对所述储层缝洞的模型进行网格划分；获取所述储层缝洞的边界条件；获取所述储层缝洞的地层岩石材料的力学参数；获取所述储层缝洞的应力状况；根据所述边界条件、所述力学参数和所述应力状况，建立所述储层缝洞的组合模型；建立溶洞垮塌模拟与裂缝闭合模拟的耦合条件，对所述储层缝洞的组合模型的变形状况进行模拟分析；其中，所述储层缝洞的组合模型是采用上述所述的缝洞组合模型的构建方法建立的。

[0040] 在一个优选的实施例中，所述边界条件是根据研究区构造演化史、区域应力变化特征和生产动态状况获取的。

[0041] 在一个优选的实施例中，所述力学参数是通过测井曲线计算或岩石物理实验获取的。

[0042] 在一个优选的实施例中，所述应力状况是通过临近井测试获取或通过高精度地应力场获取。

[0043] 作为对地应力场获取方法的进一步改进，所述地应力场的模型的获取方法包括：根据研究区区域地质调查结果，建立地质模型；结合岩石变形的力学特性，给定相应的岩石变形物理方程及对应的力学参数和边界条件，将所述地质模型转化为力学模型；通过不断调整边界作用力方式和大小、方向，利用有限元方法对所述力学模型进行模拟运算，当模拟结果与已有地应力实测值、地震震源机制或井孔地应力达到最佳拟合时，所述力学模型为所述地应力场的模型。

[0044] 本发明提出的缝洞组合模型的构建方法不仅考虑了溶洞整体控制模型，而且考虑了裂缝整体控制模型，同时引入损伤参数对控制模型进行修正，使得缝洞组合模型更加接近于实际待分析的储层缝洞，使用本发明的缝洞组合模型对实际的储层缝洞进行分析，考虑了应力特征和流体渗流情况下的缝洞组合变形，可以根据油藏应力状况结合渗流特征，对油藏开发过程中引起的缝洞组合变形特征进行预测，为开发调整提供依据。

附图说明

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

[0046] 图1为本发明的构建缝洞组合模型的方法流程图；

[0047] 图2为获得溶洞整体控制模型的流程图；

[0048] 图3为获得裂缝整体控制模型的流程图；

[0049] 图4为本发明储层缝洞变形的分析方法的流程图；

[0050] 图5为实施例一中的垂直裂缝和圆形溶洞组合模型；

[0051] 图6a为L＝5m时，最大水平主应力分布示意图；

[0052] 图6b为L＝5m时，损伤区域的分布示意图；

[0053] 图6c为L＝5m时，裂缝最大宽度处a、c两点的闭合位移变化图；

[0054] 图7a为L＝10m时，最大水平主应力分布示意图；

[0055] 图7b为L＝10m时，损伤区域的分布示意图；

[0056] 图7c为L＝10m时，裂缝最大宽度处a、c两点的闭合位移变化图；

[0057] 图8a为L＝15m时，最大水平主应力分布示意图；

[0058] 图8b为L＝15m时，损伤区域的分布示意图；

[0059] 图8c为L＝15m时，裂缝最大宽度处a、c两点的闭合位移变化图；

[0060] 图9a为地层压力为59.99Mpa时，缝洞组合的状态示意图；

[0061] 图9b为地层压力为57.7Mpa时，缝洞组合的状态示意图；

[0062] 图9c为地层压力为53.9Mpa时，缝洞组合的状态示意图；

[0063] 图10为地应力场模拟方法示意图；

[0064] 图11为水平方向构造作用与水平面内均匀剪切构造作用叠加示意图。

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

具体实施方式

[0066] 以下将结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。

[0067] 同时，在以下说明中，出于解释的目的而阐述了许多具体细节，以提供对本发明实施例的彻底理解。然而，对本领域的技术人员来说显而易见的是，本发明可以不用这里的具体细节或者所描述的特定方式来实施。

[0068] 当前对碳酸盐岩缝洞变形的分析方法都仅仅是基于概念模型来对单一溶洞的垮塌或裂缝的闭合进行研究，以此分析静态参数对缝洞变形的影响。然而这些方法均没有考虑流体渗流作用溶洞或裂缝变形的影响，因此很难对缝洞变形进行准确地预判。

[0069] 针对上述问题，本发明提出了一种综合考虑应力特征与流体渗流情况下的缝洞变形分析方法，该方法是基于缝洞组合模型来对缝洞变形进行分析的，因此本发明首先提出了一种缝洞组合模型的构建方法，图1示出了本实施例中构建缝洞组合模型的流程图。

[0070] 如图1所示，该方法首先在步骤S11中获得溶洞整体控制模型。

[0071] 假设不考虑岩体温度的变化，岩体和流体均为微压缩介质，岩体变形属于小变形范畴，渗流服从达西定律。

[0072] 在本实施例中，获得溶洞整体控制模型的步骤如图2所示。

[0073] 如图2所示，首先在步骤S111中获得弹性平衡方程。

[0074] 本实施例中，根据有效应力原理，待分析储层的岩体内各点的应力函数可以利用如下表达式进行表示：

[0075] {σ}＝[D]{εe}-α{M}p  (1)

[0076] 其中，[D]表示刚度矩阵，{εe}表示岩体内任一点的弹性应变，α表示Biot系数，p表示孔隙压力，{M}表示Kronecker符号，{σ}表示岩体内任一点的应力。

[0077] 待分析储层的平衡模型用矩阵可以表示为：

[0078]

[0079] 其中， 表示算子矩阵，{σ}表示岩体内任一点的应力，{b}表示体力向量。

[0080] 这样，将表达式(1)代入表达式(2)中，得到：

[0081]

[0082] 对表达式(3)进行时间离散，可得：

[0083] [Kuu]{ΔU}+[Kup]{Δp}＝{ΔRu}  (4)

[0084] 其中，

[0085] [Kuu]＝∑[kuu]＝∑∫Ω[B]T[D][B]dΩ  (5)

[0086]

[0087] 其中，[D]、[B]、 分别表示刚度矩阵、几何矩阵和插值函数矩阵，α表示Biot系数，p表示孔隙压力，{M}表示Kronecker符号，[Kuu]表示节点位移与节点力的整体关系矩阵，[Kup]表示孔隙压力与节点力的整体关系矩阵，{ΔRu}表示节点力增量，[kuu]表示节点位移与节点力的单元关系矩阵，Ω表示积分区域，[kup]表示孔隙压力与节点力的单元关系矩阵，{ΔU}表示位移增量。

[0088] 由此便得到了待分析储层的弹性平衡方程(4)。

[0089] 其次，在步骤S112中获取待分析储层的流体质量守恒方程和岩石骨架质量守恒方程，并根据这两个方程确定出待分析储层的连续性方程。

[0090] 具体地，待分析储层的运动方程可以表示为：

[0091]

[0092] 其中，v表示流体流动速度，μ表示流体的动粘度，K表示渗透系数张量，p表示流体压力。

[0093] 流体压缩系数αf可以根据如下表达式计算得到：

[0094]

[0095] 其中，Vf表示流体体积。

[0096] 对表达式(8)进行积分，可以得到：

[0097]

[0098] 其中，ρf表示流体密度， 表示流体初始密度，p0表示初始孔隙压力。

[0099] 岩体骨架压缩系数αs可以根据如下表达式计算得到：

[0100]

[0101] 其中，Vs表示岩石骨架体积。

[0102] 对表达式(10)进行积分，可以得到：

[0103]

[0104] 其中，ρs表示岩石骨架密度， 表示岩石骨架初始密度。

[0105] 当Δp＝p-p0的取值不大时，表达式(9)和表达式(11)可以近似表示为：

[0106]

[0107]

[0108] 本实施例中，质量守恒方程包括：流体质量守恒方程和岩石骨架质量守恒方程。其中，流体质量守恒方程可以表示为：

[0109]

[0110] 其中，φ表示岩体的孔隙度，t表示时间，εv表示岩体体积应变，kx、ky和kz分别表示岩石在直角坐标系下x，y，z向的渗透系数，μ表示岩体的动力粘度。

[0111] 岩石骨架质量守恒方程可以表示为：

[0112]

[0113] 表达式(14)与表达式(15)相加可得：

[0114]

[0115] 其中，ki表示岩体对坐标的渗透系数，Kf和Ks分别表示流体体积模量和岩石骨架体积模量，p.ii表示孔隙压力P对坐标的二阶微分。

[0116] 对表达式(16)进行时间离散，可以得到：

[0117] [Kup]T{ΔU}+[Kpp]{Δp}＝{ΔRp}  (17)

[0118] 其中，存在：

[0119]

[0120]

[0121] 其中，K表示渗透系数张量，α表示Biot系数， 表示插值函数矩阵，[B]表示几何矩阵，{M}表示Kronecker符号，[Kup]表示孔隙压力与节点力的整体关系矩阵，[Kpp]表示孔隙压力与渗流节点载荷的整体关系矩阵，kp表示孔隙压力对应的单元刚度矩阵；kpp表示孔隙压力与渗流节点载荷的单元关系矩阵；αm表示Biot系数矩阵，kup表示孔隙压力与节点力的单元关系矩阵，θ表示计算系数，k表示单元刚度矩阵， 表示单元几何矩阵， 表示差值函数矩阵，ΔRp表示渗流节点载荷增量，Δp表示孔隙压力增量，ΔU表示位移增量，Ω表示积分区域。

[0122] 由此便构建得到了待分析储层的连续性方程式(17)。

[0123] 如图2所示，在步骤S113中获得初始溶洞整体控制模型。

[0124] 由于有限元平衡方程与连续方程具有耦合项{ΔU}和{Δp}，因此需要进行联立求解，因此应用本实施例所得到的弹性平衡方程(4)和连续性方程(17)构建待分析储层的溶洞整体控制方程。

[0125] 再次，获得的的总体控制方程可以表示为：

[0126]

[0127] 其中，K表示渗透系数张量，[Kup]表示孔隙压力与节点力的整体关系矩阵，[Kpp]表示孔隙压力与渗流节点载荷的整体关系矩阵，[Kuu]表示节点位移与节点力的整体关系矩阵，Δp表示孔隙压力增量，ΔU表示位移增量，{ΔRu}表示节点力增量，{ΔRp}表示渗流节点载荷增量。

[0128] 由此便构建得到了待分析储层的初始溶洞整体控制模型式(20)。

[0129] 如图2所示，在步骤S114中引入损伤参数，获得修正后的溶洞整体控制模型。

[0130] 塑性损伤耦合模型能够真实表达实际岩石的力学行为，即塑性应变和损伤两个变量在计算时会相互影响，当单元的塑性应变达到给定的损伤值时，单元开始发生损伤。

[0131] 材料在出现损伤后，弹性刚度发生退化，为了更加准确地反应这一特征，本实施例所提供的方法引入了损伤因子标量ξ来作为损伤标量，其中，损伤因子标量ξ满足：

[0132] 0≤ξ≤1  (21)

[0133] 因此，本实施例所提供的方法在构建得到待分析储层的溶洞整体控制模型后，通过引入损伤因子标量(即损伤参数)来对构建得到的溶洞整体控制模型进行修正。具体地，计算待分析储层的损伤参数。

[0134] 在损伤及孔隙水压的作用下，Mohr-Coulomb准则可以表示为：

[0135]

[0136] 其中，τ表示岩石某平面上的剪应力，c*表示粘聚力，σn表示该平面上的正应力，pw表示孔隙水压，φ*表示内摩擦角。

[0137] 当岩石单元的等效塑性应变超过极限塑性应变 时，则可以认为该岩石单元受到破坏，将等效塑性应变用一阶指数衰减函数进行归一化，可得：

[0138]

[0139] 其中， 为归一化的等效塑性应变，a为材料常数。

[0140] 在得到待分析储层的损伤参数后，利用该损伤参数对在步骤S113中得到的溶洞整体控制模型(即初始溶洞整体控制模型)进行修正，得到修正后的溶洞整体控制模型。

[0141] 具体地，本实施例中，首先根据损伤参数计算存在损伤时单元的弹性模量；随后利用存在损伤时单元的弹性模量来替换弹性平衡模型中刚度矩阵中的不存在损伤时单元的弹性模量，从而得到修正后的刚度矩阵；最后根据修正后的刚度矩阵来确定出修正后的溶洞整体控制模型。

[0142] 考虑材料损伤的情况下，材料弹性指标可以根据如下表达式计算得到：

[0143] Eξ＝E0(1-ξ)  (24)

[0144] 其中，E0表示未存在损伤的单元的弹性模量，Eξ表示存在损伤的单元的弹性模量。

[0145] 以存在损伤的单元的弹性模量Eξ替换刚度矩阵[D]中的E，则可以得到材料损伤时的刚度矩阵 这样便存在：

[0146]

[0147] 其中， 表示材料损伤时节点位移与节点力的整体关系矩阵， 表示材料损伤时节点位移与节点力的单元关系矩阵，[B]表示几何矩阵， 表示材料损伤时的刚度矩阵，Ω表示积分区域。

[0148] 这样也就可以得到存在损伤时渗流应力耦合有限元整体控制方程：

[0149]

[0150] 其中，K表示渗透系数张量。

[0151] 至此也就构建得到了基于岩石破裂损伤因子的溶洞整体控制模型(26)。

[0152] 需要指出的是，在本发明的其他实施例中，上述步骤S111与步骤S112还可以根据实际需要进行调整(例如首先确定连续性方程，再确定弹性平衡方程)，本发明不限于此。

[0153] 如图1所示，在步骤S12中获得裂缝整体控制模型。

[0154] 裂缝真实的流体流动过程中，空隙中流体承担的压力的变化，引起裂缝周围固体颗粒的有效应力变化，由此引起岩石渗透率、孔隙比等表征渗流的水力参数的变化，这些水力参数的变化又影响裂缝孔隙流体的流场分布。因此，该裂缝整体控制模型需要通过如图3所述的方法获得，图3示意出了本实施例中获得裂缝整体控制模型的方法。

[0155] 在本实施例中，如图3所示，首先在步骤S121中获得裂缝流体流动数学模型。

[0156] 裂缝中流体的流动遵循达西定律或福熙海默定律，低速流动一般使用线性关系的达西定律，较高流动速度时，使用非线性的福熙海默定律。

[0157] 达西定律用于描述粘滞力较大的层流通过多孔介质的流动速度与水力梯度满足的关系，在一维条件下满足以下关系式：

[0158]

[0159] 式中，v'为平均流动速度；A流体流过的面积，Q为流量；k为渗透系数；H为测压水头， z表示参考面之上的高度，ρw表示流体密度，pw表示孔隙流体压力

[0160] 式(27)表示的渗流速度是一种理想的状态，当考虑流体的饱和度s和岩石的有效孔隙率φe时，式(27)变形为：

[0161] sφev＝-kgradH＝kJ  (28)

[0162] 因此，各向异性岩土介质的达西定律可表示为：

[0163] sφevi＝kijJj  (29)

[0164] 式中vi为i向量方向的流体速度，Jj为j向量方向的测压水头的梯度，kij为渗透张量，记作[K]：

[0165]

[0166] 至此，也就得到裂缝流体流动数学模型(29)。

[0167] 如图3所示，在步骤S122中获得裂缝体积变形控制模型。

[0168] 假设裂缝中固体颗粒和孔隙中流体的变形是相互独立的，那么在外部压力σp的作用下，裂缝的体积应变满足第一应变函数：

[0169]

[0170] 其中Kv为这个多孔介质的体积模量，εv为体积应变量。

[0171] 由于裂缝的整体变形等于固体颗粒和孔隙流体的变形的总和，故有：

[0172]

[0173] 其中指标s和f分别代表固体颗粒和孔隙流体，Vtot代表总体积。

[0174] 当考虑岩石的孔隙度φ、固体颗粒的体积模量Ks和孔隙流体的体积模量Kf时，则式(31)变形为第二应变函数：

[0175]

[0176] 将第一应变函数(30)与第二应变函数(32)联立，得到三个体积模量之间的关系，如下：

[0177]

[0178] 至此，得到了裂缝体积变形控制模型(33)。

[0179] 如图3所示，在步骤S123中获得孔隙介质的有效应力模型。

[0180] 岩石是由固体骨架和相互连通的孔隙及充满于骨架孔隙中的流体组成的多孔介质。岩石中的流体承担或传递的压力称为孔隙压力，裂缝周边固体颗粒所承受的应力称为有效应力。

[0181] 采用Biot有效应力，存在：

[0182]

[0183] 式中，σij为总应力，δij为Kronecker符号，σij'为有效应力， 为流体的平均压力，α为Biot系数：

[0184]

[0185] 通常取α＝1。

[0186] 至此，得到了孔隙介质的有效应力模型(34)。

[0187] 如图3所示，在步骤S124中获得裂缝渗透系数控制模型。

[0188] 在外载荷作用下，孔隙介质的体积和几何形状发生了改变，从而孔隙比和渗透性发生了改变。因此渗透系数是应力、应变或者孔隙比的函数。

[0189] 岩土介质的变形分析中，孔隙比对岩石的渗透性影响是非常重要的。

[0190] 采用Casagrande提出的关系来表达渗透系数和孔隙比，则二者满足以下关系式：

[0191] k＝1.4k0.85e2  (35)

[0192] 式中，k为孔隙比为e时渗透系数，k0.85表示孔隙比为0.85时的渗透系数值，孔隙比e和孔隙度φ的转换关系表达式为

[0193] 通过室内实验或其它手段可以得出渗透系数与应力之间的关系，在本实施例中，二者满足以下关系式：

[0194] k＝k0exp(-βσ')  (36)

[0195] 式中，k为渗透系数，σ'＝γh-p，γh为上覆岩层的重力，p为裂缝周边孔隙压力，β为常数。

[0196] 至此，得到了裂缝渗透系数控制模型(36)。

[0197] 如图3所示，在步骤S125中获得初始裂缝整体控制模型，同时引入损伤损伤参数(23)，获得修正后的裂缝整体控制模型。

[0198] 具体地，联立式(29)(33)(34)(36)，同时引入损伤参数式(23)，建立考虑损伤的渗流应力耦合有限元整体控制方程：

[0199]

[0200] 其中，K表示渗透系数张量， 表示孔隙压力与节点力的整体关系矩阵，表示孔隙压力与渗流节点载荷的整体关系矩阵， 表示节点位移与节点力的整体关系矩阵，ΔPf表示孔隙压力增量，ΔUf表示位移增量， 表示节点力增量，表示渗流节点载荷增量。

[0201] 至此，得到了修正后的裂缝整体控制模型(37)。

[0202] 需要指出的是，在本发明的其他实施例中，上述步骤S121、步骤S122、步骤S123和步骤S124还可以根据实际需要进行调整，本发明不限于此。

[0203] 如图1所示，在获得溶洞整体控制模型(17)和裂缝整体控制模型(37)后，需要对二者进行耦合叠加，从而获得缝洞组合模型。

[0204] 在本实施例中，当溶洞整体控制模型和雷锋整体控制模型满足如下耦合条件时，将二者叠加获得缝洞组合模型，该耦合叠加的条件为：

[0205] |p-pf|≤5％*pt和|v-vf|≤1％*v，或

[0206] |p-pf|≤1％*p和|v-vf|≤1％*v，

[0207] 其中，p为所述溶洞整体控制模型在每个时间步得到的溶洞周边孔隙压力，pf为所述裂缝整体控制模型在与所述溶洞整体控制模型相同时间步得到的裂缝周边孔隙压力，pt为每个时间步的实际压力值，所述实际压力值通过查询开发过程中生产井的测压数据获得，v为所述溶洞整体控制模型在每个时间步得到的溶洞内流体流速，vf为所述裂缝整体控制模型在与所述溶洞整体控制模型相同时间步得到的裂缝内流体流速。

[0208] 当溶洞整体控制模型和裂缝整体控制模型满足以上条件时，计算机会继续模拟计算，直到达到裂缝闭合或溶洞垮塌状态，如果不满足上述条件，则需要返回修改边界参数，直至符合上述条件。

[0209] 需要指出的是，在本发明的其他实施例中，上述步骤S11与步骤S12还可以根据实际需要进行调整(例如首先确定裂缝整体控制模型，再确定溶洞整体控制模型)，本发明不限于此。

[0210] 在本实施例中，基于上述得到的缝洞组合模型，提供了待分析的储层缝洞的分析方法，如图4示出了该分析方法的流程图。

[0211] 首先，在步骤S201中构建待分析的储层缝洞的模型，本实施例中，构建的模型如图5所示，图5中示出了垂直裂缝和圆形溶洞的组合模型，模型尺寸为100m*100m，圆形溶洞直径为4m，裂缝宽度为2m，裂缝长度为L，在本实施例中，裂缝长度L设置了3个值，分别为5m、
10m和15m。模型建好后，对该模型进行网格划分，网格类型可以为三角形或四边形等。

[0212] 然后，在步骤S202中获取待分析储层缝洞的边界条件，边界条件可以根据研究区构造演化史、区域应力变化特征和生产动态状况获取。在本实施例中，设置的边界条件为，如图5所示：模型上边界31施加200Mpa垂直主应力，右边界41施加180Mpa水平最大主应力；下边界32约束垂直方向位移，左边界42约束水平方向位移。

[0213] 在步骤S203中，通过测井曲线计算或岩石物理实验获取该储层缝洞的地层岩石材料的力学参数，紧接着在步骤S204中，通过临近井测试获取或通过高精度地应力场获取该储层缝洞的应力状况。在本实施例中，给整个模型施加地层压力，地层压力F满足如下关系式F＝60-0.01*y(Mpa)。这样，裂缝和溶洞内表面施加初始地层压力为60Mpa，之后其按一定速率衰减至51Mpa。根据上述的构建缝洞组合模型的方法构建该储层缝洞的组合模型，根据构建的缝洞组合模型分析地层压力衰减过程中模型应力场分布及变形结果。具体地，总的模型参数及边界条件归纳如表1所示：

[0214] 表1模型边界条件

[0215]

[0216] 在考虑损伤的情况下，得到裂缝长度L分别为5m、10m、15m时溶洞与裂缝组合模型的应力与变形规律模拟图，由于边界条件设置合适，在模拟过程中，溶洞整体控制模型和裂缝整体控制模型始终满足耦合叠加条件。

[0217] 在本实施例中，如图6a为L＝5m时，最大水平主应力分布示意图，如图6b为L＝5m时，损伤区域的分布示意图，如图6c为L＝5m时，裂缝最大宽度处a、c两点的闭合位移；如图7a为L＝10m时，最大水平主应力分布示意图，如图7b为L＝10m时，损伤区域的分布示意图，如图7c为L＝10m时，裂缝最大宽度处a、c两点的闭合位移；如图8a为L＝15m时，最大水平主应力分布示意图，如图8b为L＝15m时，损伤区域的分布示意图，如图8c为L＝15m时，裂缝最大宽度处a、c两点的闭合位移。从图6a、图7a和图8a中可以看出，随着地层压力的降低，导致裂缝和溶洞附近的最大水平主应力出现拉应力，因此造成裂缝和溶洞的损伤，如图6b、图7b和图8b所示，在裂缝的上部尖端和溶洞左下和右下两侧出现损伤，其中L＝5m时，损伤面积是4.21m2；L＝10m时，损伤面积是6.71m2；L＝15m时，损伤面积是7.40m2，从而得出，在同样的边界条件下，裂缝越长，损伤面积越大，即损伤越严重。参照图6c、图7c和图8c，可以看出，随着压力降低，裂缝出现闭合，其中L＝5m时，a、c两点最大闭合位移是0.038m；L＝10m时，a、c两点最大闭合位移是0.024m；L＝15m时，a、c两点最大闭合位移是0.032m。

[0218] 图9a为地层压力为59.99Mpa时，缝洞组合的状态示意图；图9b为地层压力为57.7Mpa时，缝洞组合的状态示意图；图9c为地层压力为53.9Mpa时，缝洞组合的状态示意图。从图9a中可以看出，当地层压力从60Mpa开始降低时，溶洞和裂缝无变形，随着地层压力降低，当地层压力降低到57.7Mpa时，开始出现垮塌现象，当地层压力降低到53.9Mpa时，垮塌面积达到17.45m2。在本实施例中，定义17.45m2为最大垮塌面积，这里满足开始垮塌的极限地层压力为57.7Mpa。

[0219] 通过模型模拟，可以得知储层缝洞组合随压力变化的状态，得知开始垮塌的极限地层压力，从而对油藏开发过程中引起的缝洞组合变形进行预测，为油藏开发调整提供了依据。

[0220] 实施例二：

[0221] 与实施例一不同的是，在本实施例中，通过高精度地应力场获取储层缝洞的应力状况，以下将详细说明高精度地应力场的获得方法。

[0222] 区域应力场的总体规律可以在调查断层的活动特征、地应力实测的基础上得出初步认识。但要定量地反映区域应力场，找出应力集中部位，则需要在地应力实测基础上通过区域应力场的有限元数值模拟、反演来得到。根据现有已知的地应力实测点的应力资料和震源机制解资料来推求整个计算区域的地应力场。

[0223] 具体方法是：首先，根据区域地质调查结果，建立研究区的地质力学模型；然后，通过不断改变边界力作用方式和大小、方向来模拟计算区域应力场，使区域介质体内应力计算结果与已有地应力实测结果(最大主应力大小和方向)达到最佳拟合。最后，得出反映研究区现今应力——形变场的区域应力场的真实情况，如图10所示。

[0224] 下面结合图10详细说明地应力场模拟的具体步骤：

[0225] (1)根据研究区区域地质调查结果，特别是研究区的区域地球动力学背景、地层岩性、构造特征和构造演化过程等建立地质模型。

[0226] 地质模型的建立主要是通过三维可视化建模软件，按照建模行业标准实现。

[0227] (2)结合岩石变形的力学特性，给定相应的岩石变形物理方程及对应的力学参数和边界条件，从而将地质模型转化为力学模型。

[0228] 将建模软件构建的地质模型通过有限元软件转化为有限元计算模型，进而对目的层区块的三维地应力场分布规律进行计算和展布。

[0229] (3)利用有限单元方法对力学模型进行模拟运算，在不违背研究区所处区域地球动力学背景的前提下，通过不断适当调整边界作用力方式和大小，使模拟结果与已有地应力实测值、地震震源机制或井孔地应力等达到最佳拟合。

[0230] 在地质力学模型构建、单元离散划分的基础上，以单井实测地应力分析为基本约束，通过反演分析确定计算模型的边界作用荷载，进而计算分析研究工区地应力场的分布规律。边界作用荷载的作用方式及大小是地应力场反演分析的关键，数值模拟分析过程中，重力场可通过设置地层的容重来实现，而构造应力场则需分析具体的构造状态并通过合理设置计算模型的构造作用边界来实现。其中，远场边界构造作用可认为是如图11所示两种基本构造状态的叠加结果：第一种：沿X方向与Y方向的水平挤压或水平拉伸构造作用；第二种：水平面内的均匀剪切构造变形作用。

[0231] 边界面采用位移加载方式，则模型中各边界面上的位移作用可表示为式(38)。

[0232]

[0233] 式中， 分别为X方向、Y方向边界面水平加载位移矢量； 分别为X方向、Y方向水平挤压或拉伸作用位移矢量； 为水平剪切构造变形矢量。

[0234] 试算边界条件的过程为：施加一组挤压力和剪切力，对模型进行求解，提取已知井点位置的三个主应力大小以及最大水平主应力方向，将其与由水力压裂资料和成像测井资料得到的结果进行比对，如果比对结果在误差允许范围内，则认为该边界条件为真实地层边界条件，否则修改边界条件再次计算，直到达到与井点实测值误差小于5％。

[0235] (4)根据三维地应力场的真实情况，获取溶洞垮塌及裂缝闭合所需的应力值的大小。

[0236] 最后说明的是，以上实施例仅用于说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。