CN110688785B - 基于平面波震源的Krauklis波数值模拟方法及设备 - Google Patents

基于平面波震源的Krauklis波数值模拟方法及设备 Download PDF

Info

Publication number
CN110688785B
CN110688785B CN201910768689.5A CN201910768689A CN110688785B CN 110688785 B CN110688785 B CN 110688785B CN 201910768689 A CN201910768689 A CN 201910768689A CN 110688785 B CN110688785 B CN 110688785B
Authority
CN
China
Prior art keywords
seismic source
particle vibration
vibration
displacement
acceleration
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910768689.5A
Other languages
English (en)
Other versions
CN110688785A (zh
Inventor
丁拼搏
刘海浩
狄帮让
魏建新
李向阳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum Beijing
Original Assignee
China University of Petroleum Beijing
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Petroleum Beijing filed Critical China University of Petroleum Beijing
Priority to CN201910768689.5A priority Critical patent/CN110688785B/zh
Publication of CN110688785A publication Critical patent/CN110688785A/zh
Application granted granted Critical
Publication of CN110688785B publication Critical patent/CN110688785B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明实施例提供了一种基于平面波震源的Krauklis波数值模拟方法及设备。其中,所述方法包括:构建粘性流体饱和单裂缝数字模型,对所述粘性流体饱和单裂缝数字模型进行单元网格剖分,得到离散后的有限元方程,确定数字模拟所需边界条件参数,根据所述边界条件参数,简化所述离散后的有限元方程,得到最终有限元离散方程;确定震源辅助矩阵和震源辅助向量,结合隐式迭代法,求解所述最终有限元离散方程,实现平面波震源的Krauklis波数值模拟。本发明实施例提供的基于平面波震源的Krauklis波数值模拟方法及设备,可以实现平面波震源的Krauklis波数值模拟。

Description

基于平面波震源的Krauklis波数值模拟方法及设备
技术领域
本发明实施例涉及Krauklis波研究技术领域,尤其涉及一种基于平面波震源的Krauklis波数值模拟方法及设备。
背景技术
地下裂缝广泛分布于地壳和上地幔空间,裂缝是地下流体的储集空间和运移通道,在地震学、地热学、地震勘探、油气开发、煤炭开采、垃圾填埋、核废料处理以及二氧化碳存储等都有重大的影响。地震波是研究裂缝的重要媒介,但是裂缝的存在会使得地震波场极为复杂,目前利用地震波检测裂缝的常用方法中,无论是叠后属性类方法,还是各向异性或者横波分裂类方法,只能定性识别裂缝,都难以提取裂缝的准确信息,特别是裂缝的几何信息,包括裂缝的开度和长度。
Krauklis波是一种裂缝空间中传播的导波。但是有限长裂缝中Krauklis波数值模拟始终是一个难点,这其中的困难之一就是裂缝发育的深部地下,地震波的波前面更接近于平面波,而目前弹性波数字模拟主要基于球面波进行;此外,裂缝介质中Krauklis波的数字模拟是一个典型的多尺度问题,模型尺度、子波尺度、裂缝尺度量级不同。必须选择合理的参数,在保证精度的前提下,尽可能的提高求解效率。因此,如何获取一种能够妥善处理模型尺度、子波尺度和裂缝尺度量级的各项异性,并且是基于平面波震源的Krauklis波数值模拟方法,就成为业界亟待解决的技术问题。
发明内容
针对现有技术存在的上述问题,本发明实施例提供了一种基于平面波震源的Krauklis波数值模拟方法及设备。
第一方面,本发明的实施例提供了一种基于平面波震源的Krauklis波数值模拟方法,包括:构建粘性流体饱和单裂缝数字模型,对所述粘性流体饱和单裂缝数字模型进行单元网格剖分,得到离散后的有限元方程,确定数字模拟所需边界条件参数,根据所述边界条件参数,简化所述离散后的有限元方程,得到最终有限元离散方程;确定震源辅助矩阵和震源辅助向量,结合隐式迭代法,求解所述最终有限元离散方程,实现平面波震源的Krauklis波数值模拟。
进一步地,在上述方法实施例内容的基础上,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟方法,所述构建粘性流体饱和单裂缝数字模型,包括:
Figure BDA0002172846440000021
Figure BDA0002172846440000022
Figure BDA0002172846440000023
其中,K’为体积模量;μ为剪切模量;ρ为密度;η为粘滞系数;ux为质点振动在x方向上的位移分量;uy为质点振动在y方向上的位移分量;fx为力源载荷的x分量;fy为力源载荷的y分量。
进一步地,在上述方法实施例内容的基础上,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟方法,所述离散后的有限元方程,包括:
Figure BDA0002172846440000024
其中,
Figure BDA0002172846440000025
为离散后的质点振动位移;M为质量矩阵;C为阻尼矩阵;K为刚度矩阵;F为力源载荷;R为与边界条件相关的项。
进一步地,在上述方法实施例内容的基础上,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟方法,所述确定数字模拟所需边界条件参数,相应地,所述边界条件参数,包括:地震子波、子波主频、子波延续时长、采样时长间隔和采样步长。
进一步地,在上述方法实施例内容的基础上,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟方法,所述最终有限元离散方程,包括:
Figure BDA0002172846440000031
其中,
Figure BDA0002172846440000032
为离散后的质点振动位移;M为质量矩阵;C为阻尼矩阵;K为刚度矩阵。
进一步地,在上述方法实施例内容的基础上,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟方法,所述确定震源辅助矩阵和震源辅助向量,结合隐式迭代法,求解所述最终有限元离散方程,包括:在子波延续时长范围内,采用一时刻质点振动的位移、速度和加速度,计算第一中间变量;根据所述第一中间变量,在子波延续时长结束后,采用另一时刻质点振动的位移、速度和加速度,计算第二中间变量。
进一步地,在上述方法实施例内容的基础上,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟方法,根据所述第二中间变量,求解所述另一时刻之后时刻的质点振动位移、速度和加速度,包括:
Figure BDA0002172846440000033
Figure BDA0002172846440000034
Figure BDA0002172846440000035
其中,
Figure BDA0002172846440000036
为k+1时刻质点振动位移;
Figure BDA0002172846440000037
为k+1时刻质点振动速度;
Figure BDA0002172846440000038
为k+1时刻质点振动加速度;β和γ为常数;Δt为采样时长间隔;
Figure BDA0002172846440000039
Figure BDA00021728464400000310
为所述第二中间变量;M为质量矩阵;C为阻尼矩阵;K为刚度矩阵。
第二方面,本发明的实施例提供了一种基于平面波震源的Krauklis波数值模拟装置,包括:
有限元方程获取模块,用于构建粘性流体饱和单裂缝数字模型,对所述粘性流体饱和单裂缝数字模型进行单元网格剖分,得到离散后的有限元方程,确定数字模拟所需边界条件参数,根据所述边界条件参数,简化所述离散后的有限元方程,得到最终有限元离散方程;
解析模块,用于确定震源辅助矩阵和震源辅助向量,结合隐式迭代法,求解所述最终有限元离散方程,实现对平面波震源的Krauklis波数值模拟。
第三方面,本发明的实施例提供了一种电子设备,包括:
至少一个处理器;以及
与处理器通信连接的至少一个存储器,其中:
存储器存储有可被处理器执行的程序指令,处理器调用程序指令能够执行第一方面的各种可能的实现方式中任一种可能的实现方式所提供的基于平面波震源的Krauklis波数值模拟方法。
第四方面,本发明的实施例提供了一种非暂态计算机可读存储介质,非暂态计算机可读存储介质存储计算机指令,计算机指令使计算机执行第一方面的各种可能的实现方式中任一种可能的实现方式所提供的基于平面波震源的Krauklis波数值模拟方法。
本发明实施例提供的基于平面波震源的Krauklis波数值模拟方法及设备,通过对粘性流体饱和单裂缝数字模型进行离散,并加载边界条件,得到最终有限元离散方程,在确定震源辅助矩阵和震源辅助向量后,采用隐式迭代法,求解最终有限元离散方程,就可以实现平面波震源的Krauklis波数值模拟。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图做一简单的介绍,显而易见地,下面描述中的附图是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的基于平面波震源的Krauklis波数值模拟方法流程图;
图2为本发明实施例提供的外部平面波震源激发的裂缝内部x分量情况示意图;
图3为本发明实施例提供的外部平面波震源激发的裂缝内部y分量情况示意图;
图4为本发明实施例提供的基于平面波震源的Krauklis波数值模拟装置结构示意图;
图5为本发明实施例提供的电子设备的实体结构示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。另外,本发明提供的各个实施例或单个实施例中的技术特征可以相互任意结合,以形成可行的技术方案,但是必须是以本领域普通技术人员能够实现为基础,当技术方案的结合出现相互矛盾或无法实现时,应当认为这种技术方案的结合不存在,也不在本发明要求的保护范围之内。
本发明实施例提供了一种基于平面波震源的Krauklis波数值模拟方法,参见图1,该方法包括:
101、构建粘性流体饱和单裂缝数字模型,对所述粘性流体饱和单裂缝数字模型进行单元网格剖分,得到离散后的有限元方程,确定数字模拟所需边界条件参数,根据所述边界条件参数,简化所述离散后的有限元方程,得到最终有限元离散方程;具体地,所述粘性流体饱和单裂缝数字模型的尺寸会先行确定,并且会确定流体类型、声速、粘性,裂缝开度、长度,背景岩石的纵横波速度、密度等岩石物理参数。对所述粘性流体饱和单裂缝数字模型进行单元网格剖分,得到离散后的有限元方程,需要选取合适的网格类型、及相关参数,进行有限元建模,剖分有限单元网格。在单元网格剖分过程中,网格参数的选取原则是:在保证裂缝尖端得到精细离散的前提下,总的单元个数尽可能少,以提高计算效率并节省内存。
102、确定震源辅助矩阵和震源辅助向量,结合隐式迭代法,求解所述最终有限元离散方程,实现平面波震源的Krauklis波数值模拟。隐式迭代的方法,具体可以参见(4)、(5)和(6)式。
基于上述方法实施例的内容,作为一种可选的实施例,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟方法,所述构建粘性流体饱和单裂缝数字模型,包括:
Figure BDA0002172846440000061
Figure BDA0002172846440000062
Figure BDA0002172846440000063
其中,K’为体积模量;μ为剪切模量;ρ为密度;η为粘滞系数;ux为质点振动在x方向上的位移分量;uy为质点振动在y方向上的位移分量;fx为力源载荷的x分量;fy为力源载荷的y分量;
Figure BDA0002172846440000064
为质点振动的x分量速度;
Figure BDA0002172846440000065
为质点振动的x分量加速度;
Figure BDA0002172846440000066
为质点振动的y分量速度;
Figure BDA0002172846440000067
为质点振动的y分量加速度。
基于上述方法实施例的内容,作为一种可选的实施例,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟方法,所述离散后的有限元方程,包括:
Figure BDA0002172846440000068
其中,
Figure BDA0002172846440000069
为离散后的质点振动位移;M为质量矩阵;C为阻尼矩阵;K为刚度矩阵;F为力源载荷;R为与边界条件相关的项。具体地,
Figure BDA00021728464400000610
为离散后的质点振动速度;
Figure BDA00021728464400000611
为离散后的质点振动加速度;F表征外力的影响。其中,(2)式是在(1)式的基础上,进行离散化后得到的有限元离散方程。
基于上述方法实施例的内容,作为一种可选的实施例,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟方法,所述确定数字模拟所需边界条件参数,相应地,所述边界条件参数,包括:地震子波、子波主频、子波延续时长、采样时长间隔和采样步长。具体地,确定数字模拟时所用的地震子波w、子波主频fw,以及子波延续时长Nw,再结合模型参数、网格参数确定合适的采样时长间隔Δt,以及采样步长N。因为方程组求解过程中采用的是隐式求解格式,迭代格式是无条件稳定的,因此可以采用较大的时间采样间隔,以提高计算效率。
基于上述方法实施例的内容,作为一种可选的实施例,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟方法,所述最终有限元离散方程,包括:
Figure BDA0002172846440000071
其中,
Figure BDA0002172846440000072
为离散后的质点振动位移;M为质量矩阵;C为阻尼矩阵;K为刚度矩阵。具体地,(3)式成立的前提是需要假定没有外力输入,则略去(2)式中的载荷项F。因为载荷项表征的是系统中外力的大小、方向等信息,当没有外力时F=0;此外,在引入边界条件以后,(2)式就被化简为(3)式的形式,(3)式表征的是没有外力条件下,地震波在模型系统内部的传播状态。
基于上述方法实施例的内容,作为一种可选的实施例,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟方法,所述确定震源辅助矩阵和震源辅助向量,结合隐式迭代法,求解所述最终有限元离散方程,包括:在子波延续时长范围内,采用一时刻质点振动的位移、速度和加速度,计算第一中间变量;根据所述第一中间变量,在子波延续时长结束后,采用另一时刻质点振动的位移、速度和加速度,计算第二中间变量。具体地,质点振动的位移、速度和加速度向量
Figure BDA0002172846440000073
初始值均为0。采用隐式Newmark算法推导了(3)式的迭代求解格式,采样时长间隔为Δt,在子波延续时长范围以内,即0≤k≤Nw时,采用k时刻质点振动的位移、速度和加速度向量
Figure BDA0002172846440000074
计算第一中间变量
Figure BDA0002172846440000075
具体计算公式如下:
Figure BDA0002172846440000081
Figure BDA0002172846440000082
在子波延续时长结束以后,即Nw<k≤N时,采用k时刻质点振动的位移、速度和加速度向量
Figure BDA0002172846440000083
计算第二中间变量
Figure BDA0002172846440000084
具体计算公式如下:
Figure BDA0002172846440000085
其中,(4)式和(5)式中,矩阵BC和向量b0中节点的排列顺序与质点振动的位移、速度、加速度向量相匹配。
基于上述方法实施例的内容,作为一种可选的实施例,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟方法,根据所述第二中间变量,求解所述另一时刻之后时刻的质点振动位移、速度和加速度,包括:
Figure BDA0002172846440000086
其中,
Figure BDA0002172846440000087
为k+1时刻质点振动位移;
Figure BDA0002172846440000088
为k+1时刻质点振动速度;
Figure BDA0002172846440000089
为k+1时刻质点振动加速度;β和γ为常数;Δt为采样时长间隔;
Figure BDA00021728464400000810
Figure BDA00021728464400000811
为所述第二中间变量;M为质量矩阵;C为阻尼矩阵;K为刚度矩阵。需要说明的是,(6)中的中间变量
Figure BDA00021728464400000812
采用的是(5)中的第二中间变量
Figure BDA00021728464400000813
作为常数,β=0.25,γ=0.5。
外部平面波震源激发的裂缝内部x分量情况可以参见图2,y分量情况可以参见图3。沿着裂缝,在裂缝内部布设检波器,记录内部X分量(如图2所示)和Y分量(如图3所示)的波场。从X分量和Y分量地震记录中,可以清晰的识别出Krauklis波同相轴,如图2和图3中分别所示。
本发明实施例提供的基于平面波震源的Krauklis波数值模拟方法,通过对粘性流体饱和单裂缝数字模型进行离散,并加载边界条件,得到最终有限元离散方程,在确定震源辅助矩阵和震源辅助向量后,采用隐式迭代法,求解最终有限元离散方程,就可以实现平面波震源的Krauklis波数值模拟。
本发明各个实施例的实现基础是通过具有处理器功能的设备进行程序化的处理实现的。因此在工程实际中,可以将本发明各个实施例的技术方案及其功能封装成各种模块。基于这种现实情况,在上述各实施例的基础上,本发明的实施例提供了一种基于平面波震源的Krauklis波数值模拟装置,该装置用于执行上述方法实施例中的基于平面波震源的Krauklis波数值模拟方法。参见图4,该装置包括:
有限元方程获取模块401,用于构建粘性流体饱和单裂缝数字模型,对所述粘性流体饱和单裂缝数字模型进行单元网格剖分,得到离散后的有限元方程,确定数字模拟所需边界条件参数,根据所述边界条件参数,简化所述离散后的有限元方程,得到最终有限元离散方程;
解析模块402,用于确定震源辅助矩阵和震源辅助向量,结合隐式迭代法,求解所述最终有限元离散方程,实现平面波震源的Krauklis波数值模拟。
本发明实施例提供的基于平面波震源的Krauklis波数值模拟方法,采用有限元方程获取模块和解析模块,通过对粘性流体饱和单裂缝数字模型进行离散,并加载边界条件,得到最终有限元离散方程,在确定震源辅助矩阵和震源辅助向量后,采用隐式迭代法,求解最终有限元离散方程,就可以实现平面波震源的Krauklis波数值模拟。
需要说明的是,本发明提供的装置实施例中的装置,除了可以用于实现上述方法实施例中的方法外,还可以用于实现本发明提供的其他方法实施例中的方法,区别仅仅在于设置相应的功能模块,其原理与本发明提供的上述装置实施例的原理基本相同,只要本领域技术人员在上述装置实施例的基础上,参考其他方法实施例中的具体技术方案,通过组合技术特征获得相应的技术手段,以及由这些技术手段构成的技术方案,在保证技术方案具备实用性的前提下,就可以对上述装置实施例中的装置进行改进,从而得到相应的装置类实施例,用于实现其他方法类实施例中的方法。例如:
基于上述装置实施例的内容,作为一种可选的实施例,本发明实施例中提供的基于平面波震源的Krauklis波数值模拟装置,还包括:中间变量计算模块,用于在子波延续时长范围内,采用一时刻质点振动的位移、速度和加速度,计算第一中间变量;根据所述第一中间变量,在子波延续时长结束后,采用另一时刻质点振动的位移、速度和加速度,计算第二中间变量。
本发明实施例的方法是依托电子设备实现的,因此对相关的电子设备有必要做一下介绍。基于此目的,本发明的实施例提供了一种电子设备,如图5所示,该电子设备包括:至少一个处理器(processor)501、通信接口(Communications Interface)504、至少一个存储器(memory)502和通信总线503,其中,至少一个处理器501,通信接口504,至少一个存储器502通过通信总线503完成相互间的通信。至少一个处理器501可以调用至少一个存储器502中的逻辑指令,以执行如下方法:构建粘性流体饱和单裂缝数字模型,对所述粘性流体饱和单裂缝数字模型进行单元网格剖分,得到离散后的有限元方程,确定数字模拟所需边界条件参数,根据所述边界条件参数,简化所述离散后的有限元方程,得到最终有限元离散方程;确定震源辅助矩阵和震源辅助向量,结合隐式迭代法,求解所述最终有限元离散方程,实现平面波震源的Krauklis波数值模拟。
此外,上述的至少一个存储器502中的逻辑指令可以通过软件功能单元的形式实现并作为独立的产品销售或使用时,可以存储在一个计算机可读取存储介质中。基于这样的理解,本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来,该计算机软件产品存储在一个存储介质中,包括若干指令用以使得一台计算机设备(可以是个人计算机,服务器,或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。例如包括:构建粘性流体饱和单裂缝数字模型,对所述粘性流体饱和单裂缝数字模型进行单元网格剖分,得到离散后的有限元方程,确定数字模拟所需边界条件参数,根据所述边界条件参数,简化所述离散后的有限元方程,得到最终有限元离散方程;确定震源辅助矩阵和震源辅助向量,结合隐式迭代法,求解所述最终有限元离散方程,实现平面波震源的Krauklis波数值模拟。而前述的存储介质包括:U盘、移动硬盘、只读存储器(ROM,Read-Only Memory)、随机存取存储器(RAM,Random AccessMemory)、磁碟或者光盘等各种可以存储程序代码的介质。
以上所描述的装置实施例仅仅是示意性的,其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的,作为单元显示的部件可以是或者也可以不是物理单元,即可以位于一个地方,或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性的劳动的情况下,即可以理解并实施。
通过以上的实施方式的描述,本领域的技术人员可以清楚地了解到各实施方式可借助软件加必需的通用硬件平台的方式来实现,当然也可以通过硬件。基于这样的理解,上述技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来,该计算机软件产品可以存储在计算机可读存储介质中,如ROM/RAM、磁碟、光盘等,包括若干指令用以使得一台计算机设备(可以是个人计算机,服务器,或者网络设备等)执行各个实施例或者实施例的某些部分所述的方法。
附图中的流程图和框图显示了根据本发明的多个实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。基于这种认识,流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分,所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意,在有些作为替换的实现方式中,方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如,两个连续的方框实际上可以基本并行地执行,有时也可以按相反的顺序执行,这依所涉及的功能而定。也要注意的是,框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合,可以用执行规定的功能或动作的专用的基于硬件的系统来实现,或者可以用专用硬件与计算机指令的组合来实现。
在本专利中,术语"包括"、"包含"或者其任何其它变体意在涵盖非排它性的包含,从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素,而且还包括没有明确列出的其它要素,或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下,由语句"包括……"限定的要素,并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
最后应说明的是:以上实施例仅用以说明本发明的技术方案,而非对其限制;尽管参照前述实施例对本发明进行了详细的说明,本领域的普通技术人员应当理解:其依然可以对前述各实施例所记载的技术方案进行修改,或者对其中部分技术特征进行等同替换;而这些修改或者替换,并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (7)

1.一种基于平面波震源的Krauklis波数值模拟方法,其特征在于,包括:
构建粘性流体饱和单裂缝数字模型,对所述粘性流体饱和单裂缝数字模型进行单元网格剖分,得到离散后的有限元方程,确定数字模拟所需边界条件参数,根据所述边界条件参数,简化所述离散后的有限元方程,得到最终有限元离散方程;
确定震源辅助矩阵和震源辅助向量,结合隐式Newmark算法,求解所述最终有限元离散方程,实现平面波震源的Krauklis波数值模拟;
所述构建粘性流体饱和单裂缝数字模型,包括:
Figure FDA0003013925090000011
其中,
Figure FDA0003013925090000012
Figure FDA0003013925090000013
其中,K’为体积模量;μ为剪切模量;ρ为密度;η为粘滞系数;ux为质点振动在x方向上的位移分量;uy为质点振动在y方向上的位移分量;fx为力源载荷的x分量;fy为力源载荷的y分量;
Figure FDA0003013925090000014
为质点振动的x分量速度;
Figure FDA0003013925090000015
为质点振动的x分量加速度;
Figure FDA0003013925090000016
为质点振动的y分量速度;
Figure FDA0003013925090000017
为质点振动的y分量加速度;
所述确定震源辅助矩阵和震源辅助向量,结合隐式Newmark算法,求解所述最终有限元离散方程,实现平面波震源的Krauklis波数值模拟,包括:
根据一时刻质点振动的位移、速度和加速度,利用震源辅助矩阵和震源辅助向量,求得中间变量;根据所述中间变量,求解另一时刻的质点振动位移、速度和加速度,以此迭代,实现平面波震源的Krauklis波数值模拟;
其中,在子波延续时长范围内,采用一时刻质点振动的位移、速度和加速度,利用公式(1)计算中间变量;
Figure FDA0003013925090000021
其中,
Figure FDA0003013925090000022
为中间变量;BC为震源辅助矩阵;b0为震源辅助向量;β和γ为常数;Δt为采样时长间隔;w为地震子波;
Figure FDA0003013925090000023
为k时刻质点振动的位移;
Figure FDA0003013925090000024
为k时刻质点振动的速度;
Figure FDA0003013925090000025
为k时刻质点振动的加速度向量;Nw为波延续时长;
在子波延续时长结束后,采用一时刻质点振动的位移、速度和加速度,利用公式(2)计算中间变量;
Figure FDA0003013925090000026
Figure FDA0003013925090000027
其中,
Figure FDA0003013925090000028
为中间变量;BC为震源辅助矩阵;β和γ为常数;Δt为采样时长间隔;
Figure FDA0003013925090000029
为k时刻质点振动的位移;
Figure FDA00030139250900000210
为k时刻质点振动的速度;
Figure FDA00030139250900000211
为k时刻质点振动的加速度向量;Nw为波延续时长。
2.根据权利要求1所述的基于平面波震源的Krauklis波数值模拟方法,其特征在于,所述离散后的有限元方程,包括:
Figure FDA00030139250900000212
其中,
Figure FDA00030139250900000213
为离散后的质点振动位移;M为质量矩阵;C为阻尼矩阵;K为刚度矩阵;F为力源载荷;R为与边界条件相关的项,包括地震子波、子波主频、子波延续时长、采样时长间隔和采样步长;
Figure FDA00030139250900000214
为离散后的质点振动速度;
Figure FDA00030139250900000215
为离散后的质点振动加速度。
3.根据权利要求1所述的基于平面波震源的Krauklis波数值模拟方法,其特征在于,所述最终有限元离散方程,包括:
Figure FDA00030139250900000216
其中,
Figure FDA0003013925090000031
为离散后的质点振动位移;M为质量矩阵;C为阻尼矩阵;K为刚度矩阵;
Figure FDA0003013925090000032
为离散后的质点振动速度;
Figure FDA0003013925090000033
为离散后的质点振动加速度。
4.根据权利要求1所述的基于平面波震源的Krauklis波数值模拟方法,其特征在于,所述根据所述中间变量,求解另一时刻的质点振动位移、速度和加速度,包括:
Figure FDA0003013925090000034
Figure FDA0003013925090000035
Figure FDA0003013925090000036
其中,
Figure FDA0003013925090000037
为k+1时刻质点振动位移;
Figure FDA0003013925090000038
为k+1时刻质点振动速度;
Figure FDA0003013925090000039
为k+1时刻质点振动加速度;β和γ为常数;Δt为采样时长间隔;
Figure FDA00030139250900000310
Figure FDA00030139250900000311
为中间变量;M为质量矩阵;C为阻尼矩阵;K为刚度矩阵。
5.一种基于平面波震源的Krauklis波数值模拟装置,其特征在于,包括:
有限元方程获取模块,用于构建粘性流体饱和单裂缝数字模型,对所述粘性流体饱和单裂缝数字模型进行单元网格剖分,得到离散后的有限元方程,确定数字模拟所需边界条件参数,根据所述边界条件参数,简化所述离散后的有限元方程,得到最终有限元离散方程;
解析模块,用于确定震源辅助矩阵和震源辅助向量,结合隐式Newmark算法,求解所述最终有限元离散方程,实现平面波震源的Krauklis波数值模拟;
所述构建粘性流体饱和单裂缝数字模型,包括:
Figure FDA00030139250900000312
其中,
Figure FDA00030139250900000313
Figure FDA0003013925090000041
其中,K’为体积模量;μ为剪切模量;ρ为密度;η为粘滞系数;ux为质点振动在x方向上的位移分量;uy为质点振动在y方向上的位移分量;fx为力源载荷的x分量;fy为力源载荷的y分量;
Figure FDA0003013925090000042
为质点振动的x分量速度;
Figure FDA0003013925090000043
为质点振动的x分量加速度;
Figure FDA0003013925090000044
为质点振动的y分量速度;
Figure FDA0003013925090000045
为质点振动的y分量加速度;
所述确定震源辅助矩阵和震源辅助向量,结合隐式Newmark算法,求解所述最终有限元离散方程,实现平面波震源的Krauklis波数值模拟,包括:
根据一时刻质点振动的位移、速度和加速度,利用震源辅助矩阵和震源辅助向量,求得中间变量;根据所述中间变量,求解另一时刻的质点振动位移、速度和加速度,以此迭代,实现平面波震源的Krauklis波数值模拟;
其中,在子波延续时长范围内,采用一时刻质点振动的位移、速度和加速度,利用公式(1)计算中间变量;
Figure FDA0003013925090000046
其中,
Figure FDA0003013925090000047
为中间变量;BC为震源辅助矩阵;b0为震源辅助向量;β和γ为常数;Δt为采样时长间隔;w为地震子波;
Figure FDA0003013925090000048
为k时刻质点振动的位移;
Figure FDA0003013925090000049
为k时刻质点振动的速度;
Figure FDA00030139250900000410
为k时刻质点振动的加速度向量;Nw为波延续时长;
在子波延续时长结束后,采用一时刻质点振动的位移、速度和加速度,利用公式(2)计算中间变量;
Figure FDA00030139250900000411
Figure FDA00030139250900000412
其中,
Figure FDA00030139250900000413
为中间变量;BC为震源辅助矩阵;β和γ为常数;Δt为采样时长间隔;
Figure FDA0003013925090000051
为k时刻质点振动的位移;
Figure FDA0003013925090000052
为k时刻质点振动的速度;
Figure FDA0003013925090000053
为k时刻质点振动的加速度向量;Nw为波延续时长。
6.一种电子设备,其特征在于,包括:
至少一个处理器、至少一个存储器、通信接口和总线;其中,
所述处理器、存储器、通信接口通过所述总线完成相互间的通信;
所述存储器存储有可被所述处理器执行的程序指令,所述处理器调用所述程序指令,以执行如权利要求1至4任一项权利要求所述的方法。
7.一种非暂态计算机可读存储介质,其特征在于,所述非暂态计算机可读存储介质存储计算机指令,所述计算机指令使所述计算机执行如权利要求1至4中任一项权利要求所述的方法。
CN201910768689.5A 2019-08-20 2019-08-20 基于平面波震源的Krauklis波数值模拟方法及设备 Active CN110688785B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910768689.5A CN110688785B (zh) 2019-08-20 2019-08-20 基于平面波震源的Krauklis波数值模拟方法及设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910768689.5A CN110688785B (zh) 2019-08-20 2019-08-20 基于平面波震源的Krauklis波数值模拟方法及设备

Publications (2)

Publication Number Publication Date
CN110688785A CN110688785A (zh) 2020-01-14
CN110688785B true CN110688785B (zh) 2021-05-28

Family

ID=69108572

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910768689.5A Active CN110688785B (zh) 2019-08-20 2019-08-20 基于平面波震源的Krauklis波数值模拟方法及设备

Country Status (1)

Country Link
CN (1) CN110688785B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113408161B (zh) * 2021-03-26 2023-05-16 中国石油大学(北京) 水力压裂激发Kraukls波模拟方法及系统
CN113341455B (zh) * 2021-06-24 2024-02-09 中国石油大学(北京) 一种粘滞各向异性介质地震波数值模拟方法、装置及设备

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102116869A (zh) * 2011-02-12 2011-07-06 中国石油大学(华东) 高精度叠前域最小二乘偏移地震成像技术
CN103760603A (zh) * 2014-01-28 2014-04-30 中国石油大学(北京) 转换波地震数据的叠前时间偏移方法及装置
DE102014213972A1 (de) * 2014-06-06 2015-12-17 Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. Vorrichtung und Verfahren zur Bestimmung von Rissparametern
CN108416182A (zh) * 2018-05-31 2018-08-17 长安大学 一种基于定量分析的明沟隔振的设计方法
CN109444954A (zh) * 2018-12-26 2019-03-08 中国矿业大学(北京) 裂缝数值的模拟方法、装置、电子设备及存储介质

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102116869A (zh) * 2011-02-12 2011-07-06 中国石油大学(华东) 高精度叠前域最小二乘偏移地震成像技术
CN103760603A (zh) * 2014-01-28 2014-04-30 中国石油大学(北京) 转换波地震数据的叠前时间偏移方法及装置
DE102014213972A1 (de) * 2014-06-06 2015-12-17 Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. Vorrichtung und Verfahren zur Bestimmung von Rissparametern
CN108416182A (zh) * 2018-05-31 2018-08-17 长安大学 一种基于定量分析的明沟隔振的设计方法
CN109444954A (zh) * 2018-12-26 2019-03-08 中国矿业大学(北京) 裂缝数值的模拟方法、装置、电子设备及存储介质

Also Published As

Publication number Publication date
CN110688785A (zh) 2020-01-14

Similar Documents

Publication Publication Date Title
Fukuda et al. Development of a GPGPU‐parallelized hybrid finite‐discrete element method for modeling rock fracture
Jin et al. An edge‐based strain smoothing particle finite element method for large deformation problems in geotechnical engineering
Franceschini et al. Block preconditioning for fault/fracture mechanics saddle-point problems
WO2019238451A1 (en) A method and a system for modelling and simulating a fractured geological structure
CA2947410A1 (en) Fast viscoacoustic and viscoelastic full-wavefield inversion
Mohammadnejad et al. An extended finite element method for fluid flow in partially saturated porous media with weak discontinuities; the convergence analysis of local enrichment strategies
CN110688785B (zh) 基于平面波震源的Krauklis波数值模拟方法及设备
EP3096252A2 (en) Adaptive multiscale multi-fidelity reservoir simulation
AU2010363352A1 (en) Systems and methods for generating updates of geological models
US20210096276A1 (en) Model for Coupled Porous Flow and Geomechanics for Subsurface Simulation
van Wees et al. 3-D mechanical analysis of complex reservoirs: a novel mesh-free approach
Lloberas-Valls et al. Strain injection techniques in dynamic fracture modeling
EP3400546B1 (en) Effective permeability upscaling for a discrete fracture network
CN110954950B (zh) 地下横波速度反演方法、装置、计算设备及存储介质
Li et al. Calibration of the discrete element method and modeling of shortening experiments
Jolfaei et al. Sensitivity analysis of effective parameters in borehole failure, using neural network
Ramos et al. Working with dynamic earthquake rupture models: A practical guide
CN110687589B (zh) 裂缝介质中横波激发Krauklis波的数值模拟方法及设备
Doonechaly et al. 3D hybrid tectono-stochastic modeling of naturally fractured reservoir: Application of finite element method and stochastic simulation technique
Sahin et al. Computational modeling of quasi static fracture using the nonlocal operator method and explicit phase field model
Chen et al. Numerical simulation of unstable suction transients in unsaturated soils: the role of wetting collapse
Wang et al. Implementation of efficient low-storage techniques for 3-D seismic simulation using the curved grid finite-difference method
Ashworth et al. An open source numerical framework for dual-continuum geomechanical simulation
Payne et al. Insights into pulverized rock formation from dynamic rupture models of earthquakes
Zakian et al. Spectral finite element simulation of seismic wave propagation and fault dislocation in elastic media

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant