一种预测多相孔隙介质弹性模量的高效数值模拟方法与流程

1.本发明涉及油气及煤层气勘探开发技术领域，具体涉及一种预测多相孔隙介质弹性模量的高效数值模拟方法。


背景技术：

2.在油气或煤层气勘探开发过程中，评估储层岩石中多相流体对岩石声学属性的影响是地震资料成像、解释及油气或者煤层气资源储量评估的重要约束条件。获得岩石中流体分布与岩石声学属性关系的方法主要包括实验室直接测量法和数值模拟预测法：实验室测量法优点在于承认了岩石的各向异性及非均匀性特征，测量结果直接、准确，是对数值模拟方法的重要约束和验证。然而，较高的实验成本使实验室直接测量方法无法在研究者中广泛推广；数值模拟法相对于高成本的实验法而言，具有广义性，计算成本低等特点。因此，在过去数十年，众多研究者针对其展开了数量庞大的研究。
3.white等首先提出了流体分布会造成岩石介观尺度的非均质性，从而导致岩石的声波速度出现频散和衰减。除此之外，white等认为气-水交互的层状模型和均匀球状气泡分布模型等价于大尺度非均质体对地震波属性产生的影响。dutta&ode等将white的球状模型从二维情况拓展为三维情况。对这类模型的求解表明：在物理上，低频biot纵波通过孔弹性介质的表现特征等价于通过更大尺度的粘弹性固体的表现特征。据此，一些研究者通过特定模型对孔弹性介质的地震响应进行了解析求解：steven r.pride和berryman的双孔模型；norris的层状模型；johnson的任意几何形状的斑块饱和模型；m
ü
ller的随机多孔介质模型；toms的随机离散和连续流体分布模型。
4.上述解析模型能很好的处理满足常规几何形状的非均质孔弹性模型。但是，对于较为复杂的非均质情况，人们通过数值模拟波传播过程进行研究。虽然数值模拟方法计算精度高，可以预测更加复杂的情况，但是，其在计算上昂贵，甚至有时无法完成计算。这是因为：需要非常精细的网格来表示非均质性；与地震波长相比，扩散长度非常小，因此在低频范围内，与流体压力平衡相关的扩散过程的分辨率是一个关键问题。为了克服限制，masson&pride提出了一种不同以往且有趣的数值计算方法。他们将随着时间变化的应力施加到非均质样品的边界上，并通过数值计算平均应力和应变场，确定其有效复模量。数值结果说明模拟过程的稳定性与岩石中流体粘度密切相关。然而，在时间域求解biot方程时无法避免小的时间积分步长问题，同时需要将方程按照硬孔和非硬孔分开，这使得整个计算过程处于冗余状态。此外，对于频变模型，除了动态渗透率，方程中所有弹性模量均与频率相关，这将会在时间域求解时造成额外的卷积计算。相反，在频率域求解biot方程则不需要考虑这类问题。在频率域可以避免卷积计算，使整个方程的求解过程变方便。rubino等人还提出了一种升尺度数值模拟方法，以获得多相流体饱和多孔弹性介质的频散和衰减特征。该方法基于空间频率域中经典biot方程的有限元解模拟振荡压缩性和剪切测试。虽然该数值方法大幅度提高了计算效率，但是在有限元计算过程中需要求取存在二阶导数的流体和固体位移方程，不利于有限元方法计算；除此之外，对于复杂的流体分布，尤其岩石模型中
存在细小气泡的三维模型，常规网格剖分导致资源浪费巨大，最终常常无法计算，获得有效结果。


技术实现要素：

5.为了克服现有技术方案的不足，本发明提供一种预测多相孔隙介质弹性模量的高效数值模拟方法，利用biot位移-压力方程进行有限元算法高效求解非均质模型频散衰减特征。为了实现上述技术目的，本发明采用如下技术方案：
6.一种预测多相孔隙介质弹性模量的高效数值模拟方法，包括如下步骤：
7.s1：通过流体分布图像获得第i相流体分布范围ωi，i＝1，2
…
n；n为混合流体的最大相数；
8.s2：定义流体密度ρf为：ρf(x，y，z)＝ρi，x，y，z∈ωi，i＝1，2
…
n，x，y，z为笛卡尔坐标系空间坐标；ρi是第i项流体密度；
9.s3：定义流体粘度为ηf：ηf＝ηi，x，y，z∈ωi，i＝1，2
…
n；ηi是第i相流体粘度；
10.s4：定义流体体积模量kf为：kf＝ki，x，y，z∈ωi，i＝1，2
…
n；
11.s5：构建频率域biot方程的位移-压力形式：
[0012][0013]
其中，s＝c：ε-αbpi；ρ
av
为饱和岩石密度；ρ
av
＝ρd φρf；ρd为岩石骨架密度，φ为孔隙度，ρc为流体复密度，τ为扭曲因子，κ
p
为岩石渗透率；s为总应力，c为干燥岩石骨架柔度，ε为应变张量，αb为biot-willis系数，p为流体压力，i为单位张量矩阵，f为周期震荡的外力；ω为角频率，u为固体位移张量；
[0014][0015]
其中，kd为岩石骨架体积模量；
[0016]
s6：约束边界条件：数值模拟的目标是实验室岩石样品，岩石样品轴向方向是z方向，垂直于z方向的平面上存在着正交的x、y方向，根据纵波模量的定义，横向应变ε
xx
＝ε
xy
＝ε
yy
＝ε
yx
＝0，及流体边界没有流体流出，因此边界条件为：
[0017][0018]
下边界为固定边界，需要同时满足两个条件：边界法向上流体加速度为0，固体位移是0，因此边界条件是：
[0019]
[0020]
上边界加载周期震荡应力，对应的边界条件是边界法向应力是σ
zz
，流体流动的加速度是0，因此边界条件是：
[0021][0022]
s7：全局粗化网格剖分，局部网格精细剖分求解；
[0023]
s8：使用有限元求解方程，获得应变固体位移场和流体压力场；
[0024]
s9：计算纵波模量m：
[0025][0026]
其中，σ
zz
是垂向应力，ε
zz
是垂向应变；
[0027]
则纵波速度v
p
和纵波衰减q
p
为：
[0028][0029]
其中，和是虚步和实步算子。
[0030]
优选地，步骤s1中流体分布图像通过ct扫描获得。
[0031]
优选地，步骤s7具体包括：
[0032]
s71：利用有限元算法使用粗化网格求解方程；
[0033]
s72：计算所有网格元素上的偏微分方程的残差；
[0034]
s73：估计所有网格元素的解中误差，其中第l方程的l2范数误差为h为局部网格尺寸，ρ
l
为第1个偏微分方程的残差，qi为偏微分方程稳定性估计导数阶数，默认值是2，由于参数c是任意常数，计算出的误差实际上是误差指示因子；
[0035]
s74：如果已进行了预设的优化次数或超过了最大有限元剖分数，则终止执行；
[0036]
s75：建立局部误差指示因子，其中a是网格的面积，根据局部误差指标的大小细化有限元；
[0037]
s76：重复步骤s71-s75。
[0038]
与现有技术相比，本发明的有益效果：
[0039]
本发明提供的预测多相孔隙介质弹性模量的高效数值模拟方法，相比于以前求解孔弹性biot方程，实现对部分饱和岩石样品频散和衰减预测的有限元数值模拟方法，具有同等计算精度下，占据内存少，计算时间少的特点，特别适用于模型中存在小非均质目标，可避免常规四面体网格剖分占用大量内存，计算占用大量时间的问题。
附图说明
[0040]
为了更清楚的说明本发明实施例或现有技术的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做简单的介绍，显而易见的，下面描述中的附图仅仅是本发明的一些实施例，对于本领域中的普通技术人员来说，在不付出创造性劳动的前提下，还可
根据这些附图获得其他附图。
[0041]
图1为气水两相流的双层模型；
[0042]
图2为本发明实施例频散和衰减图；
[0043]
图3为本发明实施例水/气的流体分布模式；
[0044]
图4为本发明实施例流体体积模量、密度和粘度；
[0045]
图5为本发明实施例横向边界条件加载位置；
[0046]
图6为本发明实施例震源加载位置；
[0047]
图7为本发明实施例固定边界条件加载位置；
[0048]
图8为本发明实施例网格初步剖分；
[0049]
图9为本发明实施例网格自适应剖分；
[0050]
图10为本发明实施例在频率1000hz时的流体压力分布图；
[0051]
图11为本发明实施例体积模量随频率变化图；
[0052]
图12为本发明实施例体积衰减随频率变化图。
具体实施方式
[0053]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0054]
一种预测多相孔隙介质弹性模量的高效数值模拟方法，包括如下步骤：
[0055]
s1：通过流体分布图像获得第i相流体分布范围ωi，流体分布图像通过ct扫描获得，i＝1，2
…
n；n为混合流体的最大相数；
[0056]
s2：定义流体密度ρf为：ρf(x，y，z)＝ρi，x，y，z∈ωi，i＝1，2
…
n，x，y，z为笛卡尔坐标系空间坐标；ρi是第i项流体密度；
[0057]
s3：定义流体粘度为ηf：ηf＝ηi，x，y，z∈ωi，i＝1，2
…
n；ηi是第i相流体粘度；
[0058]
s4：定义流体体积模量kf为：kf＝ki，x，y，z∈ωi，i＝1，2
…
n；
[0059]
s5：构建频率域biot方程的位移-压力形式：
[0060][0061]
其中，ρ
av
为饱和岩石密度；ρ
av
＝ρd φρf；ρd为岩石骨架密度，φ为孔隙度，ρc为流体复密度，τ为扭曲因子，κ
p
为岩石渗透率；s为总应力，c为干燥岩石骨架柔度，ε为应变张量，αb为biot-willis系数，p为流体压力，i为单位张量矩阵，f为周期震荡的外力；ω为角频率，u为固体位移张量；
[0062][0063]
其中，kd为岩石骨架体积模量；
[0064]
s6：约束边界条件：数值模拟的目标是实验室岩石样品，岩石样品轴向方向是z方向，垂直于z方向的平面上存在着正交的x、y方向，根据纵波模量的定义，横向应变ε
xx
＝ε
xy
＝ε
yy
＝ε
yx
＝0，及流体边界没有流体流出，因此边界条件为：
[0065][0066]
下边界为固定边界，需要同时满足两个条件：边界法向上流体加速度为0，固体位移是0，因此边界条件是：
[0067][0068]
上边界加载周期震荡应力，对应的边界条件是边界法向应力是σ
zz
，流体流动的加速度是0，因此边界条件是：
[0069][0070]
s7：全局粗化网格剖分，局部网格精细剖分求解，具体包括：
[0071]
s71：利用有限元算法使用粗化网格求解方程；
[0072]
s72：计算所有网格元素上的偏微分方程的残差；
[0073]
s73：估计所有网格元素的解中误差，其中第1方程的l2范数误差为h为局部网格尺寸，ρ
l
为第1个偏微分方程的残差，qi为偏微分方程稳定性估计导数阶数，默认值是2，由于参数c是任意常数，计算出的误差实际上是误差指示因子；
[0074]
s74：如果已进行了预设的优化次数或超过了最大有限元剖分数，则终止执行；
[0075]
s75：建立局部误差指示因子，其中a是网格的面积，根据局部误差指标的大小细化有限元；
[0076]
s76：重复步骤s71-s75；
[0077]
s8：使用有限元求解方程，获得应变固体位移场和流体压力场；
[0078]
s9：计算纵波模量m：
[0079][0080]
其中，σ
zz
为垂向应力，ε
zz
为垂向应变；
[0081]
则纵波速度v
p
和纵波衰减q
p
为：
[0082][0083]
其中，和是虚步和实步算子。
[0084]
本发明对实验中气/水部分饱和indiana石灰岩中流体流动引起的频散和衰减进行预测，表1为气和水的流体物性参数，表2为干燥indiana岩石物性参数，样品饱和度为88％，其中，步骤s1水/气的流体分布模式通过ct扫描获得，并显示在图3中。
[0085]
根据ct扫描结果，该样品对应的模型计算区域为圆柱形区域。为利用本发明方法研究目标样品中流体流动引发的频散和衰减现象，需要按照以下步骤进行计算：
[0086]
样品中为气/水两相，按照步骤s2-s4，获得图4(a)所示的流体体积模量分布函数、图4(b)所示的流体密度分布函数、图4(c)所示的流体黏度分布函数。
[0087]
按照步骤s6，约束边界条件为：
[0088]
如图5所示，横向边界条件为：
[0089][0090]
如图6所示，上边界加载周期震荡应力，对应的边界条件是边界法向应力及流体流动的加速度是0，边界条件为：
[0091][0092]
如图7所示，下边界为固定边界，需要同时满足两个条件：边界法向上流体加速度为0，固体位移是0，因此边界条件是：
[0093][0094]
为求解方程，首先利用图8所示的粗化网格进行计算，根据残差构建图9所示的自适应网格。求解方程，获得indiana灰岩在不同频率下的流体压力分布，图10展示了频率为1000hz时的流体压力分布。根据步骤s9，获得频变体积模量和衰减，如图11和图12所示。
[0095]
为说明本发明的效果，以附图1(a)层状模型为例，该模型上层为空气，下层为水，表1为对应的流体物性。表2为模型岩石物性信息。结合上述物性信息，利用常规有限元计算方法求解biot流-固耦合位移方程，为保证分层界面的精度，采用细化四面体网格剖分模型，如图1(b)所示，计算结果获得的流体压场如图1(c)所示，体积模量频散和衰减变化特征如图2中实线所示，计算过程总耗时为2分32秒。
[0096]
使用本发明方法计算时，首先将所需物性信息按照步骤s1-s4构建系数矩阵，结合步骤s5的biot位移-压力方程，粗化网格剖分，如图1(d)所示。在获得初步残差后，利用自适应算法细化小目标地质体网格，如图1(e)所示，可以看到，在边界层处网格充分细化，但其他位置网格仍然为粗化网格。据此计算的流体压力场为图1(f)。所得流体压力场图1(c)和图1(f)差距不大，其对应的频散和衰减曲线为图2虚线所示，对比常规算法结果(图2实线)，可以看到两者的计算结果基本一致，但本发明所使用的计算时间为32秒，比常规方法速度提高五倍之多。
[0097]
表1流体物性参数表
[0098]
流体属性数值水体积模量2.25e9pa气体积模量1e5pa水密度1000kg/m3气密度78kg/m3水黏度0.001pa
·
s气黏度1.5e-4pa
·s[0099]
表2岩石物性参数表
[0100]
固体属性数值岩石排水体积模量2.5e10pa骨架颗粒模量7.7e10pa孔隙度0.108干骨架密度2367.4kg/m3渗透率9.8692e-18m2岩石排水剪切模量1.52e10pa扭曲因子1模型长度0.4m模型高度0.4m