一种流-固耦合的波场模拟方法、装置、设备及存储介质制造方法及图纸

本发明专利技术公开了一种流

【技术实现步骤摘要】
一种流
‑
固耦合的波场模拟方法、装置、设备及存储介质


[0001]本专利技术涉及地球物理勘探
，具体为一种流
‑
固耦合的波场模拟方法、装置、设备及存储介质。

技术介绍

[0002]地下油气储集层通常具有高度不均匀性，含有孔隙、裂缝等复杂微结构，并且孔隙、裂缝中可能充填油、气、水等多种流体。随着油气勘探越来越精细化，传统的基于单相声波/弹性波波动方程的波场模拟方法(比如有限差分法、有限元法等)，难以满足实际生产中面临的地震波在双相或者多相介质中传播机理的数值计算需求。
[0003]为了解决以上问题，一部分学者通过在经典有限差分法、有限元法等数值模拟方法的基础上，人为定义各种独特的流
‑
固耦合边界条件，并通过对流
‑
固界面附近的介质参数进行算术平均、几何平均等方式来提高流
‑
固耦合双相介质的波场计算精度(Moczo et al.,2002；van Vossen et al.,2002；Zhang et al,2004；Carcione et al,2018)；另一部分学者则提出了多种双相介质模型及与之匹配的双相介质波动方程，比如Biot模型(Biot，1956a，1956b；Ba等，2011)、喷射流模型(Mavko和Nur，1975，1979；Aki和Richards，2002；Liner，2012)、BISQ模型(Dvorkin和Nur，1993；杨宽德等，2002)等。
[0004]有限差分方法处理流r/>‑
固界面的两大缺陷：一是由于离散的原因，导致界面位置偏离实际情况，进而使得波的旅行时不对；二是由于不规则界面会被离散成阶梯状，从而导致突变点处出现人工假想或者噪声。除了有限差分方法，其他可选用的处理流
‑
固界面的方法，比如有限元法、有限体积法、边界元法等，计算量通常都比较大、且实现起来不容易。
[0005]此外，经典的Biot模型、BISQ模型等双相介质理论只是含流体孔隙介质声学的简化表达或者等效模型，故这类模型在实际生产应用中会受到双相介质模型本身以及控制波场传播机理的双相介质波动方程的一些连续性假设条件限制。基于这类模型对应的波动方程进行波场正演模拟，对实际地下复杂双相介质的适应性会受到限制、且数值计算的灵活度有待提高，这将导致其计算精度随着多孔介质复杂度的增加而降低。

技术实现思路

[0006]本专利技术的目的在于提供一种流
‑
固耦合的波场模拟方法、装置、设备及存储介质，以解决上述
技术介绍
中提出的问题。
[0007]为实现上述目的，本专利技术提供如下技术方案：
[0008]一种流
‑
固耦合的波场模拟方法，包括步骤如下：
[0009]S100、设定介质的物理参数及数值模拟参数，根据流
‑
固双相介质模型，自动判别当前计算节点介质属性；
[0010]S200、计算流体区域或固体区域的波场值，定义流
‑
固边界处的波场耦合条件，实现地震波在流体与固体区域之间双向传播；
[0011]S300、在规定时间内循环进行波场正演计算；
[0012]S400、输出所有流和固体区域中的波场数值。
[0013]进一步的，所述物理参数包括固体介质的纵横波速度、体积密度、Q值以及流体介质的纵波速度、体积密度、Q值物理参数，所述数值模拟参数包括时间采样间隔、空间采样间隔、震源参数、地震观测系统、地震记录时长正演模拟参数以及设定模拟的区域的外部边界条件。
[0014]进一步的，LBM离散模型和LSM模型采用“DdQq”形式表示，其中d表示空间维度，q表示LBM的离散速度个数；
[0015]在二维情况下，LBM可用的离散模型包括D2Q7模型、D2Q9模型、D2Q13模型，在三维情况下，LBM可用的离散模型包括D3Q14模型、D3Q19模型、D3Q27模型；
[0016]在二维情况下，LSM模型可用的模型包括D2Q8模型、D2Q12模型等，在三维情况下，LSM模型可用的模型包括D3Q14模型、D3Q26模型。
[0017]进一步的，自动判别当前计算节点介质属性包括根据流
‑
固双相介质模型，定义一个表征双相介质流体或者固体属性的相函数，通过相函数，当计算某个节点的波场值时，即可通过相函数自动判别计算中设计的相关节点的介质属性。
[0018]进一步的，根据相函数定义的介质属性，如当前节点为流体介质，则采用LBM进行粘滞声波数值模拟；如当前节点为固体介质，则采用LSM进行弹性波或者粘弹性波数值模拟。
[0019]进一步的，LBM进行粘滞声波数值模拟包括步骤：
[0020]S1、根据流体介质参数，更新各阶矩张量；
[0021]S2、计算平衡态分布函数；
[0022]S3、进行离散粒子的碰撞步与迁移步运算；
[0023]S4、采用插值处理获得粒子数分布函数；
[0024]S5、处理计算区域外边界与模型内部界面；
[0025]S6、判断循环是否终止，如是则输出计算结果，如否则返回到S1；
[0026]LSM进行弹性波或者粘弹性波数值模拟包括步骤：
[0027]A1、读入固体介质参数；
[0028]A2、根据参数建立离散弹簧模型；
[0029]A3、定义边界条件；
[0030]A4、施加震源或外力；
[0031]A5、根据Verlet算法，更新位移、速度以及作用力，并判断是否满足计算终止条件，如是则输出计算结果，如否则返回到A4。
[0032]进一步的，步骤S1采用：控制方程表示：f(x+c
i
dt,t+dt)
‑
f(x,t)＝
‑
M
‑1SM[f(x,t)
‑
f
(eq)
(x,t)]dt，其中，f(x,t)为粒子数分布函数向量，其中各个元素对应各个方向的粒子数分布函数为f
i
(x,t),i＝1,2,...,q，q为LBM离散模型的离散向量个数；f
(eq)
(x,t)为粒子数分布函数向量；c
i
为方向的离散速度；M为变换矩阵；S为松弛矩阵，其对角元素由松弛因子构成，dt是时间步长(或时间间隔)；当公式中松弛矩阵S的对角元素全都相等时，代表单松弛时间LBM模型；当松弛矩阵S的对角元素仅仅由两个数构成时，代表双松弛时间LBM模型；否则，当松弛矩阵S的对角元素由两个以上实数构成时，代表多松弛时间LBM模型；
[0033]步骤S2采用：方程(1)中沿着i方向运动的粒子的平衡态分布函数f
i(eq)
采用
Maxwell
‑
Boltzmann分布，其二阶精度形式为：
[0034][0035]其中，w
i
和c
s
分本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种流
‑
固耦合的波场模拟方法，其特征在于，包括步骤如下：设定介质的物理参数及数值模拟参数，根据流
‑
固双相介质模型，自动判别当前计算节点介质属性；计算流体区域或固体区域的波场值，定义流
‑
固边界处的波场耦合条件，实现地震波在流体与固体区域之间双向传播；在规定时间内循环进行波场正演计算；输出所有流体和固体区域中的波场数值。2.如权利要求1所述的波场模拟方法，其特征在于，所述物理参数包括固体介质的纵横波速度、体积密度、Q值以及流体介质的纵波速度、体积密度、Q值物理参数，所述数值模拟参数包括时间采样间隔、空间采样间隔、震源参数、地震观测系统、地震记录时长正演模拟参数以及设定模拟的区域的外部边界条件。3.如权利要求1所述的波场模拟方法，其特征在于，自动判别当前计算节点介质属性包括根据流
‑
固双相介质模型，定义一个表征双相介质流体或者固体属性的相函数，通过相函数，当计算某个节点的波场值时，即可通过相函数自动判别计算中设计的相关节点的介质属性。4.如权利要求3所述的波场模拟方法，其特征在于，根据相函数定义的介质属性，如当前节点为流体介质，则采用LBM进行粘滞声波数值模拟；如当前节点为固体介质，则采用LSM进行弹性波或者粘弹性波数值模拟。5.如权利要求4所述的波场模拟方法，其特征在于，LBM进行粘滞声波数值模拟包括步骤：S1、根据流体介质参数，更新各阶矩张量；S2、计算平衡态分布函数；S3、进行离散粒子的碰撞步与迁移步运算；S4、采用插值处理获得粒子数分布函数；S5、处理计算区域外边界与模型内部界面；S6、判断循环是否终止，如是则输出计算结果，如否则返回到S1；弹簧网络模型进行弹性波或者粘弹性波数值模拟包括步骤：A1、读入固体介质参数；A2、根据参数建立离散弹簧模型；A3、定义边界条件；A4、施加震源或外力；A5、根据Verlet算法，更新位移、速度以及作用力，并判断是否满足计算终止条件，如是则输出计算结果，如否则返回到A4。6.如权利要求5所述的波场模拟方法，其特征在于，步骤S1采用：控制方程表示为：f(x+c
i
dt,t+dt)
‑
f(x,t)＝
‑
M
‑1SM[f(x,t)
‑
f
(eq)
(x,t)]dt，其中，f(x,t)为粒子数分布函数向量，其中各个元素对应各个方向的粒子数分布函数为f
i
(x,t),i＝1,2,...,q，q为LBM的离散向量个数；f
(eq)
(x,t)为粒子数分布函数向量；c
i
为方向的离散速度；M为变换矩阵；S为松弛矩阵，其对角元素由松弛因子构成，dt是时间步长；当公式中松弛矩阵S的对角元素全都相等时，代表单松弛时间LBM模型；当松弛矩阵S的对角元素仅仅由两个数构成时，代表双松弛时
间LBM模型；否则，当松弛矩阵S的对角元素由两个以上实数构成时，代表多松弛时间LBM模型；步骤S2采用：方程(1)中沿着i方向运动的粒子的平衡态分布函数f
i(eq)
采用Maxwell
‑
Boltzmann分布，其二阶精度形式为：其中，w
i
和c
s
分别为各个方向的权系数和格子声速，其大小由LBM离散速度模型决定；宏观流体密度ρ和速度u分别对应粒子数分布函数对离散速度的零阶矩和一阶矩张量，表示如下：下：其中，n为粒子每次迁移的网格点数，即迁移速度；步骤S3运算采用：方程(1)描述的LBM粒子的运动状态包括碰撞步和迁移步，在某个时刻，同时运动到某个网格节点的粒子之间会发生碰撞，碰撞之后，粒子会根据碰撞后的离散速度迁移到相邻的网格节点，在迁移的过程中，若遇到波阻抗界面，则会根据以下公式定义的规则进行反射与透射：其中，R和T分别为反射与透射系数；n1和n2分别为界面两侧的粒子迁移速度；ρ1和ρ2分别为界面两侧的流体密度；步骤S4采用：当LBM中出现非整数的迁移速度n时，LBM离散粒子则会运动到非整数网格节点上，根据非整数网格节点上的粒子数分布函数通过插值算法获得整数网格节点上的粒子数分布函数；步骤S5确定采用：标准反弹、镜面反弹等对波场能量无吸收作用的处理方法，或应用完全匹配层、指数衰减吸收边界方法；步骤S6计算结果：通过在特定网格点施加扰动，然后循环进行离散粒子的碰撞步与迁移步运算，可得到流体区域在各个时刻的波场数值；步骤A2采用：对于各向同性介质，以质点i为中心的弹性体元中所蕴含的弹性能量E
i
可用下式表示：其中，K和c分别为线弹簧与角弹簧的弹性常数，θ
jik
为以质点i为顶点、ji和ki这两条线
共同构成的一个角，u
i
为质点i的位移，为由质点i指向相邻质点j的单位向量，而x
ij
则为由质点i指向相邻质点j的向量，u
j
为相邻质点j的位移，上式中，第一项为所有邻近质点j对中心位置的质点i的中心力作用或者线弹簧作用对弹性能量的总贡献，第二项为以中心质点i为顶点的所有角度上的键弯曲作用或者角弹簧作用对弹性能量的总贡献；通过将方程(6)进行线性化、进行泰勒展开保留至二阶精度，并与标准的弹性波波动方程进行对比系数，可得弹簧的弹性系数与拉梅常数之间的定量关系如下：与此同时，可得到纵横波速度的计算表达式：式中，ρ为固体介质的体密度，a为LSM离散网格边长；利用方程(6)得到能量密度，并对位移u
i
求偏导，可得LSM中心质点i所受到的相邻质点j的总的弹簧作用力F

【专利技术属性】
技术研发人员：夏木明，杨长春，王灿云，张文秀，田飞，李宗伟，
申请(专利权)人：中国科学院地质与地球物理研究所，
类型：发明
国别省市：