CN112001100B - 一种三维地震波场se-ibe耦合模拟方法 - Google Patents

一种三维地震波场se-ibe耦合模拟方法 Download PDF

Info

Publication number
CN112001100B
CN112001100B CN202010671252.2A CN202010671252A CN112001100B CN 112001100 B CN112001100 B CN 112001100B CN 202010671252 A CN202010671252 A CN 202010671252A CN 112001100 B CN112001100 B CN 112001100B
Authority
CN
China
Prior art keywords
ibe
subdomain
coupling
sub
unit
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
CN202010671252.2A
Other languages
English (en)
Other versions
CN112001100A (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.)
Tianjin University
Original Assignee
Tianjin University
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 Tianjin University filed Critical Tianjin University
Priority to CN202010671252.2A priority Critical patent/CN112001100B/zh
Publication of CN112001100A publication Critical patent/CN112001100A/zh
Application granted granted Critical
Publication of CN112001100B publication Critical patent/CN112001100B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Graphics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种三维地震波场SE‑IBE耦合模拟方法，其特征在于，构建包含有限断层震源、地壳层速度结构以及复杂场地的整体物理模型；在整体物理模型中建立包含沉积盆地的立方体模型，将立方体表面作为耦合边界面，位于立方体内的空间域设为SE子域，位于立方体外的空间域设为IBE子域；以SE子域模拟复杂场地区域，以IBE子域模拟复杂场地之外的层状半无限空间域；设SE子域内仅含散射波场，以谱单元法求解SE子域内波场；设IBE子域内同时包含输入波场和散射波场，以间接边界元法求解IBE子域内波场。本发明可以在保证精度的同时大大节省计算资源。

Description

一种三维地震波场SE-IBE耦合模拟方法
技术领域
本发明涉及一种地震波场模拟方法，特别涉及一种三维地震波场SE-IBE耦合模拟方法。
背景技术
目前，作为典型复杂场地，沉积盆地对地震动的显著放大作用(盆地效应)已在多次震害调查和强震观测中得以证实。因而，针对沉积盆地开展其地震效应研究，科学定量评估其地震地面运动，是确保工程抗震安全和推动韧性城市建设的重要课题。
开展大尺度沉积盆地地震效应研究对计算资源有着很高要求，精确高效的数值方法是研究顺利开展的保障。目前三维地震波场模拟方法主要包括域离散型方法(有限差分法、有限元法、谱单元法)和边界型方法(边界元法、边界积分方程法)等，其中域离散型方法存在人工边界问题，同时在处理大规模计算时效率较低，而边界型方法则难以处理非线性问题，同时在处理复杂模型方面也不够灵活。
发明内容
本发明为解决公知技术中存在的技术问题而提供一种高效的三维地震波场SE-IBE耦合模拟方法。
本发明为解决公知技术中存在的技术问题所采取的技术方案是：一种三维地震波场SE-IBE耦合模拟方法，构建包含有限断层震源、地壳层速度结构以及复杂场地的整体物理模型；在整体物理模型中建立包含沉积盆地的立方体模型，将立方体表面作为耦合边界面，位于立方体内的空间域设为SE子域，位于立方体外的空间域设为IBE子域；以SE子域模拟复杂场地区域，以IBE子域模拟复杂场地之外的层状半无限空间域；设SE子域内仅含散射波场，以谱单元法求解SE子域内波场；设IBE子域内同时包含输入波场和散射波场，以间接边界元法求解IBE子域内波场。
进一步地，立方体表面接近沉积盆地的轮廓面。
进一步地，该方法包括如下步骤：
步骤一，构建包含有限断层震源、地壳层速度结构以及复杂场地的整体物理模型；将模型划分成SE子域和IBE子域；
步骤二，对SE子域进行网格剖分，将其位于耦合边界面上的网格单元定义为SE耦合边界单元；对IBE子域进行网格剖分，将其位于耦合边界面上的网格单元定义为IBE耦合边界单元；
步骤三，基于频域内波动方程弱形式，获得谱单元频域离散矩阵，通过在SE耦合边界单元的耦合边界面施加一组密度不同的三个坐标方向的均布荷载，以荷载产生的波场之和建立SE子域内散射波场模型；
步骤四，基于傅里叶变换的波数积分方法得到IBE子域内输入波场；基于单层位势理论，对IBE耦合边界单元的耦合边界面施加另一组密度不同的三个坐标方向均布荷载，由荷载产生的波场之和建立IBE子域内散射波场模型；
步骤五，分别求得SE子域和IBE子域在各自耦合边界单元上的平均响应，基于SE耦合边界单元和IBE耦合边界单元共有耦合边界面，满足位移相容和牵引力平衡条件，形成最终求解矩阵。
进一步地，步骤一中，在SE子域完全包含复杂场地的情况下，IBE子域大于SE子域。
进一步地，步骤二中，对SE子域进行六面体网格剖分。
进一步地，步骤二中，对IBE子域的耦合边界面进行面单元离散。
进一步地，步骤三包括如下分步骤：
a.设耦合边界单元为标准参考单元，采用GLL积分获得标准参考单元的质量矩阵和刚度矩阵；将GLL积分节点取为插值点，在GLL积分节点上构造基函数，进而基于伽辽金法将基函数作为测试函数，得到标准参考单元控制方程；
b.集合所有标准参考单元控制方程，得到SE子域的离散线性系统矩阵；
c.将施加给每个耦合单元的均布荷载转换成节点荷载，设在耦合边界单元的耦合边界面上施加单位荷载产生的SE子域内动力响应定义为SE子域响应函数；由SE子域响应函数构造SE子域内散射波场模型。
进一步地，步骤四中，由均布荷载产生的动力格林函数来表征IBE子域内散射波场模型。
本发明具有的优点和积极效果是：本发明充分整合了谱单元(SEM)在高效处理复杂介质及复杂模型方面的优势和间接边界元(IBEM)在自动满足辐射条件及降维方面的优势，通过引入SE子域和IBE子域响应函数，实现了两子域的完全独立求解，较传统迭代耦合求解，具有两种方法互不污染的显著优势。该方法可以在保证精度的同时大大节省计算资源。本发明可用于科学定量评估沉积盆地的地震地面运动，为工程抗震安全提供参考。
附图说明
图1为大型沉积盆地整体模型示意图；
图2为SE子域波场剖分构建示意图；
图3为IBE子域波场剖分构建示意图；
图4为SE子域和IBE子域耦合示意图；
图5为斜面单元所在土层的上下层界面处于锁定状态的层状半空间中斜面均布荷载示意图；
图6为斜面单元所在土层的上下层界面处于解锁状态的层状半空间中斜面均布荷载示意图；
图7为斜面单元所在土层的上下层界面处于原始状态的层状半空间中斜面均布荷载示意图；
图8为t＝4s时刻的断层震源下沉积盆地地面运动的位移波场仿真图；
图9为t＝8s时刻的断层震源下沉积盆地地面运动的位移波场仿真图；
图10为t＝12s时刻的断层震源下沉积盆地地面运动的位移波场仿真图；
图11为t＝16s时刻的断层震源下沉积盆地地面运动的位移波场仿真图；
图12为t＝20s时刻的断层震源下沉积盆地地面运动的位移波场仿真图；
图13为t＝24s时刻的断层震源下沉积盆地地面运动的位移波场仿真图；
图14为t＝28s时刻的断层震源下沉积盆地地面运动的位移波场仿真图；
图15为t＝30s时刻的断层震源下沉积盆地地面运动的位移波场仿真图。
图2中，px、py、pz对应为SE子域侧对SE耦合边界单元的耦合边界面施加X、Y、Z坐标方向力。
图3中，qx、qy、qz对应为IBE子域侧对IBE耦合边界单元的耦合边界面施加X、Y、Z坐标方向力。
图4中，
Figure BDA0002582381300000031
对应为SE子域沿X、Y、Z坐标方向的位移平均响应；
Figure BDA0002582381300000032
对应为SE子域沿X、Y、Z坐标方向的牵引力平均响应；
Figure BDA0002582381300000033
对应为IBE子域沿X、Y、Z坐标方向的位移平均响应；
Figure BDA0002582381300000034
对应为IBE子域沿X、Y、Z坐标方向的牵引力平均响应；
Figure BDA0002582381300000035
对应为IBE子域沿X、Y、Z坐标方向在第i个耦合边界单元的位移平均响应；
Figure BDA0002582381300000036
对应为IBE子域沿X、Y、Z坐标方向在第i个耦合边界单元的牵引力平均响应。
图5至图7中：建立直角坐标系(x,y,z)，x、y、z对应为直角坐标系中的三个方向。px0、py0和pz0表示斜面均布荷载p0分解到x、y和z三个方向上后的单位均布荷载。px1、py1和pz1为使荷载所在层上顶面固定所需要施加的外力，px2、py2和pz2为使荷载所在层下底面固定所需要施加的外力。
具体实施方式
为能进一步了解本发明的发明内容、特点及功效，兹列举以下实施例，并配合附图详细说明如下：
本申请中的英文缩写的中文注释如下：
SE为Spectral Element的缩写，中文释义为谱单元。
IBE为Indirect Boundary Element的缩写，中文释义为间接边界元。
GLL为Gauss-Lobatto-Legendre缩写，中文释义为高斯-洛巴托-勒根德积分。
Lagrange中文释义为拉格朗日。
SEM为Spectral Element Method的缩写，中文释义为谱单元法。
Layer为层。
IBEM为Indirect Boundary Element Method的缩写，中文释义为间接边界元法。
请参见图1至图15，一种三维地震波场SE-IBE耦合模拟方法，其特征在于，构建包含有限断层震源、地壳层速度结构以及复杂场地的整体物理模型；在整体物理模型中建立包含沉积盆地的立方体模型，将立方体表面作为耦合边界面，位于立方体内的空间域设为SE子域，位于立方体外的空间域设为IBE子域；以SE子域模拟复杂场地区域，以IBE子域模拟复杂场地之外的层状半无限空间域；设SE子域内仅含散射波场，以谱单元法求解SE子域内波场；设IBE子域内同时包含输入波场和散射波场，以间接边界元法求解IBE子域内波场。
其中，输入波场指层状半空间有限断层震源地震反应产生的波场；散射波场指沉积盆地存在时震源地震反应引起的附加波场。
优选地，立方体表面接近沉积盆地的轮廓面。这样可使耦合边界划分尽量规则，SE子域在完全包含拟研究复杂场地的情况下尽量占有较少区域，从而在有限的资源下提高计算效率和降低存储量。
本发明分域模式的优势在于：IBE子域仅需离散规则耦合边界，无需离散自由地表和层交界面；相对于直接处理含有限断层震源的整体模型，合理分域在显著降低SE子域计算量和存储量的同时精确满足无穷远边界条件，实现有限计算资源下，精细刻画沉积盆地的同时提高模拟的频率。
谱单元法是利用快速傅立叶变换，在频域内求解连续质量分布体系的平衡微分方程，再利用傅立叶逆变换得到时域范围内结构的动力响应。
边界元法是基于控制微分方程的基本解来建立相应的边界积分方程，再结合边界的剖分而得到的离散算式。
进一步地，该方法可包括如下步骤：
步骤一，构建包含有限断层震源、地壳层速度结构以及复杂场地的整体物理模型；将模型划分成SE子域和IBE子域。
步骤二，对SE子域进行网格剖分，将其位于耦合边界面上的网格单元定义为SE耦合边界单元；对IBE子域进行网格剖分，将其位于耦合边界面上的网格单元定义为IBE耦合边界单元。
步骤三，基于频域内波动方程弱形式，获得谱单元频域离散矩阵，通过在SE耦合边界单元的耦合边界面施加一组密度不同的三个坐标方向的均布荷载，以荷载产生的波场之和建立SE子域内散射波场模型。
密度不同的三个坐标方向的均布荷载，即在X、Y、Z三个坐标方向上分别施加均布载荷，且在X、Y、Z三个坐标方向上施加的载荷密度不同。px、py、pz对应为SE子域侧对SE耦合边界单元的耦合边界面施加X、Y、Z坐标方向力；px、py、pz的载荷密度不同。
步骤四，基于傅里叶变换的波数积分方法得到IBE子域内输入波场；基于单层位势理论，对IBE耦合边界单元的耦合边界面施加另一组密度不同的三个坐标方向均布荷载，由荷载产生的波场之和建立IBE子域内散射波场模型。
qx、qy、qz对应为IBE子域侧对IBE耦合边界单元的耦合边界面施加X、Y、Z坐标方向力；qx、qy、qz的载荷密度不同。
步骤五，分别求得SE子域和IBE子域在各自耦合边界单元上的平均响应，基于SE耦合边界单元和IBE耦合边界单元共有耦合边界面，SE耦合边界单元和IBE耦合边界单元构成图2中所示的SE-IBE子域耦合边界单元，满足位移相容和牵引力平衡条件，形成最终求解矩阵。
进一步地，步骤一中，在SE子域完全包含复杂场地的情况下，IBE子域可大于SE子域。
进一步地，步骤二中，可对包含沉积盆地的立方体模型进行六面体网格剖分。六面体可为单位立方体。因为SE子域位于立方体模型内，所以相当于对SE子域进行六面体网格剖分，位于耦合边界面上的网格立方体为SE耦合边界单元。
对IBE子域通过边界元方法进行求解，而边界元方法的本质就是将虚拟单位荷载直接施加到感兴趣区域(如盆地和外部的交界面)的表面，所以可对IBE子域的耦合边界面进行面单元离散，比如可对IBE子域的耦合边界面进行正方形网格剖分，正方形的网格单元即设定为IBE耦合边界单元；正方形可为单位正方形。
用单位立方体对包含沉积盆地的立方体模型进行网格剖分，并对IBE子域的耦合边界面进行单位正方形网格剖分，这样IBE耦合边界单元的正方形的边长等于SE耦合边界单元的正立方体的棱长。
进一步地，步骤三可包括如下分步骤：
a.设耦合边界单元为标准参考单元，采用GLL积分获得标准参考单元的质量矩阵和刚度矩阵；将GLL积分节点取为插值点，在GLL积分节点上构造基函数，进而基于伽辽金法将基函数作为测试函数，得到标准参考单元控制方程。首先采用具有谱收敛性的Gauss-Lobatto-Legendre高阶基函数作为测试函数，然后采用有限元的数值离散方式，通过若干六面体网格剖分计算区域，最后在形成基函数对应的系数矩阵基础上，求解大型线性方程组，完成离散网格对应的场值计算。
b.集合所有标准参考单元控制方程，得到SE子域的离散线性系统矩阵。
c.将施加给每个耦合单元的均布荷载转换成节点荷载，设在在耦合边界单元的耦合边界面上施加单位荷载产生的SE子域内动力响应定义为SE子域响应函数；由SE子域响应函数构造SE子域内散射波场模型。
进一步地，步骤四中，可由均布荷载产生的动力格林函数来表征IBE子域内散射波场模型。将均布荷载产生的动力格林函数定义为IBE子域响应函数。
耦合边界单元的耦合边界面为水平方向或垂直方向，采用在水平方向或垂直方向的耦合边界面均布荷载格林函数，可借助格林函数间的递推性，从而显著降低IBE子域求解时间。
下面以本发明的一个优选实施例来进一步说明本发明的工作流程和工作原理：
我们将本发明方法应用于一个有限断层震源-地壳层波速结构-大型沉积盆地三维地震波场模拟。
一种三维地震波场SE-IBE耦合模拟方法，构建一个包含有限断层震源、地壳层速度结构以及杂场地整体物理模型；震源位于层状半无限空间中，断层面大小为15km×10km，其中心与自由表面的垂直距离为8km，与盆地中心的水平距离为14km。沉积盆地半径为6km，深度为2km，盆地及周围土层参数如表1所示。
表1.大型沉积盆地及周围土层参数
Figure BDA0002582381300000071
在整体物理模型中建立包含沉积盆地的立方体模型，将立方体表面作为耦合边界，位于立方体内的空间域设为SE子域，位于立方体外的空间域设为IBE子域；整体模型如图1所示。
构建有限断层震源-地壳层速度结构-复杂场地整体物理模型，开展三维地震波场模拟。具体工作流程如下：
1.采用较佳的分域方式对整体模型进行分域。
本发明对整体模型分域的特点是：以SE子域和IBE子域分别精细模拟复杂场地区域和复杂场地之外的层状半无限空间域，耦合边界划分尽量规则，SE子域在完全包含拟研究复杂场地的情况下尽量占有较少区域，从而在有限的资源下提高计算效率和降低存储量。
具体实施时，采用图1所示的分域模式(以沉积盆地为例)，在沉积盆地附近引入虚拟的立方体，立方体的表面即耦合边界。立方体之内为SE子域，而之外为IBE子域，考虑到有限断层长和宽可到几十甚至上百千米，同时深达数十千米，IBE子域将显著大于SE子域。本发明分域模式的优势在于：IBE子域仅需离散规则耦合边界，无需离散自由地表和层交界面，离散规则耦合边界面包括垂直于水平面的竖面和水平面，同时在竖面和水平面上均布荷载格林函数，可借助格林函数间的递推性，从而显著降低IBE子域求解时间；相对于直接处理含有限断层震源的整体模型，本发明的分域方法，在显著降低SE子域计算量和存储量的同时精确满足无穷远边界条件，实现在有限计算资源下，精细剖分沉积盆地的同时，还提高模拟的频率。
2.SE子域波场构建。
图2为SE子域波场构建示意图，对包含沉积盆地的立方体模型进行网格剖分，将其位于耦合边界面上的网格单元定义为耦合边界单元，SE子域内波场即可通过在所有耦合边界单元上施加密度不同的三个坐标方向均布荷载来模拟，耦合边界单元上施加单位荷载时SE子域内响应定义为响应函数。
首先基于频域内波动方程弱形式，推导谱单元频域求解的离散矩阵形式，进而通过引入响应函数构造SE子域内散射波场，具体方案如下：
2-1.谱单元法频域弱形式与离散：
由于设定SE子域内不含震源，因而不考虑体力，对频域内不含体力的波动方程的两侧乘以测试函数φm，并在物理单元域Ω上积分，物理单元域即具有有限个单元的有限空间域。可得下式(1)所示的控制方程的弱形式：
Figure BDA0002582381300000081
式中，ρ为质量密度，ω为角频率，φm为测试函数(m＝1,2,3,...,N)，τ为应力张量，
Figure BDA0002582381300000082
为位移张量，
Figure BDA0002582381300000083
为外法向矢量，Ω为物理单元域，S为域边界。
Figure BDA0002582381300000084
表示梯度算子，是一个向量微分算子。
然后将物理单元域Ω映射到标准参考单元上，并对标准参考单元采用Gauss-Lobatto-Legendre(GLL)法配制网格点，对于GLL网格点可由Lagrange多项式逼近并采用高斯积分求积。进而在伽辽金方法框架下将基函数取为测试函数，结合如式(2)的插值和GLL积分可得如式(3)所示的标准参考单元控制方程。
Figure BDA0002582381300000085
[ω2Me+Ke]Ue＝Pe (3)；
式中，Me为参考单元质量矩阵，Ke为参考单元刚度矩阵，Ue为参考单元位移向量，Pe为参考单元等效节点力向量；u1为使用GLL节点基函数在计算区域上展开的X方向的位移分量，u2为使用GLL节点基函数在计算区域上展开的Y方向的位移分量，u3为使用GLL节点基函数在计算区域上展开的Z方向的位移分量，φi(x,y,z)为参考单元坐标下的GLL节点中的第i个节点的基函数，
Figure BDA0002582381300000086
对应为第i个节点的基函数的X方向位移分量的展开系数，
Figure BDA0002582381300000087
对应为第i个节点的基函数的Y方向位移分量的展开系数，
Figure BDA0002582381300000088
对应为第i个节点的基函数的Z方向位移分量的展开系数，ω为角频率。
由于Lagrange插值和GLL积分采用相同的配置点，所以Me为对角矩阵，可大为降低计算难度，同时也易于实现并行计算。
集合所有单元，可得谱单元对于上述弱形式的离散线性矩阵系统如下式(4)：
[ω2M+K]U＝P (4)；
式中，M为全局质量矩阵，K为全局刚度矩阵，U为全局位移向量，P为全局等效节点力向量，ω为角频率。
2-2.SE子域波场构建：
SE子域内波场模型，可通过在所有耦合边界单元上施加密度不同的三个坐标方向均布荷载来模拟，定义在耦合边界单元上施加单位荷载时SE子域内响应为SE子域响应函数。
基于单层位势理论，如图2所示，通过在SE子域侧对耦合边界单元的耦合边界面施加X、Y、Z坐标方向力，以其产生波场之和模拟SE子域内散射波场。具体实施时需针对每个耦合单元将均布荷载转换成节点荷载，比如单元等效到节点，单元均布荷载转换成节点荷载时，设每个单元面有n个节点，则一个节点荷载为1/n单元面均布荷载；以三角形面单元为例，三角形面单元有3个节点则一个节点荷载为1/3三角形面单元均布荷载；当节点等效到单元时,一个节点位于n个单元上，即n个单元有一个共同的节点，则每个单元的荷载为1/n该节点的荷载。例如一个节点位于4个单元上则每个单元的荷载为1/4节点的荷载。结合式(4)，构成如式(5)所示的散射波场方程组。在耦合边界单元上施加单位荷载产生的SE子域内动力响应定义为SE子域响应函数。
Figure BDA0002582381300000091
其中，耦合单元在耦合边界面的节点为耦合节点，下标s表示SE子域的耦合节点，下标i表示SE子域的非耦合节点的其他内部节点，Kii为SE子域各内部节点形成的单元刚度矩阵，Kis为SE子域内部节点在耦合节点处形成的单元刚度矩阵，Ksi为SE子域耦合节点在内部节点处形成的单元刚度矩阵，Kss为SE子域各耦合节点形成的单元刚度矩阵，Ui为内部单元的位移向量，Us为耦合单元的位移向量。p为耦合边界单元上的荷载向量，其包括密度不同的三个坐标方向的均布荷载，R定义为均布荷载转换成节点荷载的等效转换矩阵，施加的荷载向量p将通过与IBE子域耦合求得。
基于以上方案，本发明推导的频域求解公式本质上是对惠更斯原理的实现；由于Lagrange插值和GLL积分采用相同配置点，使得全局质量矩阵为对角矩阵，由此可大大降低计算难度；由于推导中同样采用了GLL插值基函数，可在较少自由度情况下获得很高精度，使得计算效率大幅提高。
3.IBE子域波场构建。
图3为IBE子域波场构建示意图，IBE子域边界单元，即SE和IBE子域耦合边界面上单元，IBE子域内波场即可通过在所有边界单元上施加密度不同的三个坐标方向均布荷载来模拟。
层状半空间有限断层震源/点震源引起的地震反应作为IBE子域内输入波场，输入波场由基于Fourier变换的波数积分方法给出，散射波场则基于单层位势理论，如图3所示，通过IBE子域侧对耦合边界单元的耦合边界面施加X、Y、Z坐标方向力，以其产生波场之和模拟IBE子域内散射波场。基于设定的分域模式，耦合边界单元为包括垂直于水平面的竖面和平行于水平面的水平面，同时考虑地壳层状特性，拟选用层状半空间中水平面和竖面均布荷载动力格林函数作为基本解。水平面均布荷载动力格林函数可由动力刚度矩阵方法直接求得。请参见图5至图7，建立直角坐标系，直角坐标系为(x,y,z)，x、y、z对应为直角坐标系中的三个方向。px0、py0和pz0表示斜面均布荷载p0分解到x、y和z三个方向上后的单位均布荷载。px1、py1和pz1为使荷载所在层上顶面固定所需要施加的外力，px2、py2和pz2为使荷载所在层下底面固定所需要施加的外力。
竖面均布荷载的格林函数可由如下方法得到：
A.令斜面单元所在土层的上下层界面处于锁定状态，此时均布荷载只会对该固定层内部的位移和应力产生影响，对其他土层无影响。计算使土层顶面固定所需要施加的外力px1、py1和pz1，以及底面所需要施加的外力px2、py2和pz2。
B.在图6中，解除斜面单元所在土层上下层界面的锁定状态，将上步中求得的外力px1、py1、pz1、px2、py2和pz2以相反方向作用到层状半空间上，令其他层界面的外加荷载值为0，采用直接刚度法即可求得各土层交界面的位移。
C.将图5和图6的两种状态进行叠加，即可恢复到如图7所示的原始状态，这里的px0、py0和pz0表示均布荷载p0分解到x、y和z三个方向上后荷载幅值。由各土层交界面的位移可求得各土层内部的反射波幅值，进而可求得土层内任意一点的位移和应力。
层状半空间动力格林函数相对于全空间动力格林函数，具有无需离散自由表面和层交界面的优点，而均布荷载动力格林函数相对于集中荷载动力格林函数，则具有形成的IBEM不存在奇异性的优点。将均布荷载产生的动力格林函数定义为IBE子域响应函数。
4.SE子域和IBE子域耦合。
在完成SE子域和IBE子域波场构建后，分别求得SE子域和IBE子域在各自耦合边界单元上的平均响应，通过耦合边界单元满足位移相容和牵引力平衡条件，形成最终求解矩阵，该矩阵为稠密矩阵，实现SE子域和IBE子域的耦合。两子域耦合的总体原则为：通过引入SE和IBE子域响应函数，实现SE子域和IBE子域独立计算，避免两子域的迭代耦合、通过SE子域和IBE子域在耦合边界单元上满足位移相容和牵引力平衡实现两子域的耦合。
具体实施时，采用图4所示的SE-IBE耦合模式，右侧区域对应为IBE子域四边形的IBE耦合边界单元，左侧区域对应为SE子域六面体的SE耦合边界单元，SE耦合边界单元左侧为SE子域内部单元。px、py、pz对应为SE子域侧对耦合边界单元的耦合边界面施加X、Y、Z坐标方向力；qx、qy、qz对应为IBE子域侧对耦合边界单元的耦合边界面施加X、Y、Z坐标方向力；
Figure BDA0002582381300000111
对应为SE子域沿X、Y、Z坐标方向的位移平均响应；
Figure BDA0002582381300000112
Figure BDA0002582381300000113
对应为SE子域沿X、Y、Z坐标方向的牵引力平均响应；
Figure BDA0002582381300000114
对应为IBE子域沿X、Y、Z坐标方向的位移平均响应；
Figure BDA0002582381300000115
对应为IBE子域沿X、Y、Z坐标方向的牵引力平均响应；
Figure BDA0002582381300000116
对应为IBE子域沿X、Y、Z坐标方向在第i个耦合边界单元的位移平均响应；
Figure BDA0002582381300000117
对应为IBE子域沿X、Y、Z坐标方向在第i个耦合边界单元的牵引力平均响应。则有：
Figure BDA0002582381300000118
首先由式(5)，求解SE子域的位移响应函数矩阵和牵引力响应函数矩阵，设Ru为SE子域的位移响应函数矩阵，设Rt为SE子域的牵引力响应函数矩阵；在各耦合边界单元上分别施加三个坐标方向的均布荷载(以向量p表示)，得到SE耦合边界单元的平均响应，平均响应定义为SE耦合边界单元四个节点平均值，设u1为SE子域的位移平均响应，t1为SE子域侧的牵引力平均响应；然后求解IBE子域的位移和应力响应函数矩阵，设Gu为IBE子域的位移响应函数矩阵；设Gt为IBE子域的应力响应函数矩阵，在各IBE耦合边界单元上同样施加另一组三个坐标方向的均布荷载(以向量q表示)，得到IBE耦合边界单元的平均响应，设u2为IBE子域的位移平均响应，t2为IBE子域的应力平均响应；设ui为IBE子域对应第i个耦合边界单元位移平均响应，ti为IBE子域对应第i个IBE耦合边界单元应力平均响应。
最后结合IBE子域内的输入波场产生的各IBE耦合边界单元平均响应，通过满足两子域位移相容和牵引力平衡边界条件，可整理得式(6)所示耦合求解方程组：
Figure BDA0002582381300000119
基于以上分析，可看出本发明提出的SE-IBE耦合方法，通过引入SE子域响应函数，实现了SE子域和IBE子域独立计算，避免了两子域的迭代耦合，相比于传统耦合法的全域耦合方程，本文耦合方程仅涉及两子域在耦合边界面上的耦合单元，计算量显著减小。
图8至图15分别为不同时刻的断层震源下沉积盆地地面运动的位移波场仿真图，可以看出本发明较好地模拟了地震波传播过程。
以上所述的实施例仅用于说明本发明的技术思想及特点，其目的在于使本领域内的技术人员能够理解本发明的内容并据以实施，不能仅以本实施例来限定本发明的专利范围，即凡本发明所揭示的精神所作的同等变化或修饰，仍落在本发明的专利范围内。

Claims (7)

1.一种三维地震波场SE-IBE耦合模拟方法，其特征在于，构建包含有限断层震源、地壳层速度结构以及复杂场地的整体物理模型；在整体物理模型中建立包含沉积盆地的立方体模型，将立方体表面作为耦合边界面，位于立方体内的空间域设为SE子域，位于立方体外的空间域设为IBE子域；以SE子域模拟复杂场地区域，以IBE子域模拟复杂场地之外的层状半无限空间域；设SE子域内仅含散射波场，以谱单元法求解SE子域内波场；设IBE子域内同时包含输入波场和散射波场，以间接边界元法求解IBE子域内波场；
该方法包括如下步骤：
步骤一，构建包含有限断层震源、地壳层速度结构以及复杂场地的整体物理模型；将模型划分成SE子域和IBE子域；
步骤二，对SE子域进行网格剖分，将其位于耦合边界面上的网格单元定义为SE耦合边界单元；对IBE子域进行网格剖分，将其位于耦合边界面上的网格单元定义为IBE耦合边界单元；
步骤三，基于频域内波动方程弱形式，获得谱单元频域离散矩阵，通过在SE耦合边界单元的耦合边界面施加一组密度不同的三个坐标方向的均布荷载，以荷载产生的波场之和建立SE子域内散射波场模型；
步骤四，基于傅里叶变换的波数积分方法得到IBE子域内输入波场；基于单层位势理论，对IBE耦合边界单元的耦合边界面施加另一组密度不同的三个坐标方向均布荷载，由荷载产生的波场之和建立IBE子域内散射波场模型；
步骤五，分别求得SE子域和IBE子域在各自耦合边界单元上的平均响应，基于SE耦合边界单元和IBE耦合边界单元共有耦合边界面，满足位移相容和牵引力平衡条件，形成最终求解矩阵。
2.根据权利要求1所述的三维地震波场SE-IBE耦合模拟方法，其特征在于，立方体表面接近沉积盆地的轮廓面。
3.根据权利要求1所述的三维地震波场SE-IBE耦合模拟方法，其特征在于，步骤一中，在SE子域完全包含复杂场地的情况下，IBE子域大于SE子域。
4.根据权利要求1所述的三维地震波场SE-IBE耦合模拟方法，其特征在于，步骤二中，对SE子域进行六面体网格剖分。
5.根据权利要求1所述的三维地震波场SE-IBE耦合模拟方法，其特征在于，步骤二中，对IBE子域的耦合边界面进行面单元离散。
6.根据权利要求1所述的三维地震波场SE-IBE耦合模拟方法，其特征在于，步骤三包括如下分步骤：
a.设耦合边界单元为标准参考单元，采用GLL积分获得标准参考单元的质量矩阵和刚度矩阵；将GLL积分节点取为插值点，在GLL积分节点上构造基函数，进而基于伽辽金法将基函数作为测试函数，得到标准参考单元控制方程；
b.集合所有标准参考单元控制方程，得到SE子域的离散线性系统矩阵；
c.将施加给每个耦合单元的均布荷载转换成节点荷载，设在耦合边界单元的耦合边界面上施加单位荷载产生的SE子域内动力响应定义为SE子域响应函数；由SE子域响应函数构造SE子域内散射波场模型。
7.根据权利要求1所述的三维地震波场SE-IBE耦合模拟方法，其特征在于，步骤四中，由均布荷载产生的动力格林函数来表征IBE子域内散射波场模型。
CN202010671252.2A 2020-07-13 2020-07-13 一种三维地震波场se-ibe耦合模拟方法 Active CN112001100B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010671252.2A CN112001100B (zh) 2020-07-13 2020-07-13 一种三维地震波场se-ibe耦合模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010671252.2A CN112001100B (zh) 2020-07-13 2020-07-13 一种三维地震波场se-ibe耦合模拟方法

Publications (2)

Publication Number Publication Date
CN112001100A CN112001100A (zh) 2020-11-27
CN112001100B true CN112001100B (zh) 2022-07-29

Family

ID=73467596

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010671252.2A Active CN112001100B (zh) 2020-07-13 2020-07-13 一种三维地震波场se-ibe耦合模拟方法

Country Status (1)

Country Link
CN (1) CN112001100B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113221409B (zh) * 2021-05-07 2023-04-14 香港中文大学(深圳)城市地下空间及能源研究院 一种有限元和边界元耦合的声波二维数值模拟方法和装置

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105184010A (zh) * 2015-09-25 2015-12-23 天津城建大学 基于快速多极间接边界元法的高频地震波散射模拟方法
CN105869208A (zh) * 2016-03-28 2016-08-17 天津城建大学 三维饱和沉积盆地地震动模拟方法
WO2016187237A1 (en) * 2015-05-20 2016-11-24 Schlumberger Technology Corporation Inversion for tectonic stress
CN106991236A (zh) * 2017-04-05 2017-07-28 西南石油大学 一种基于四维地应力动态变化的重复压裂选井选层方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2504502A (en) * 2012-07-31 2014-02-05 Geco Technology Bv Processing wavefield data incorporating large timesteps and upscaled medium properties

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2016187237A1 (en) * 2015-05-20 2016-11-24 Schlumberger Technology Corporation Inversion for tectonic stress
CN105184010A (zh) * 2015-09-25 2015-12-23 天津城建大学 基于快速多极间接边界元法的高频地震波散射模拟方法
CN105869208A (zh) * 2016-03-28 2016-08-17 天津城建大学 三维饱和沉积盆地地震动模拟方法
CN106991236A (zh) * 2017-04-05 2017-07-28 西南石油大学 一种基于四维地应力动态变化的重复压裂选井选层方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Earthquake site effect modeling in the Granada basin using a 3-D indirect boundary element method;Jimin Lee;《Physics and Chemistry of the Earth》;20131231;正文第1-14页 *
复杂局部场地对地震波的散射间接边界积分方程-有限元耦合分析;黄磊;《中国优秀硕士学位论文全文数据库 (工程科技Ⅱ辑)》;20160915;正文第2、4章 *
点震源作用下三维沉积盆地地震动谱元模拟;刘中宪 等;《世界地震工程》;20200430;正文第200-207页 *

Also Published As

Publication number Publication date
CN112001100A (zh) 2020-11-27

Similar Documents

Publication Publication Date Title
Faccioli et al. 2D and 3D elastic wave propagation by a pseudo-spectral domain decomposition method
Taborda et al. Large-scale earthquake simulation: computational seismology and complex engineering systems
Edwards Unstructured, control-volume distributed, full-tensor finite-volume schemes with flow based grids
AU2007325904B2 (en) Electromagnetic imaging by four dimensional parallel computing
Genes Dynamic analysis of large-scale SSI systems for layered unbounded media via a parallelized coupled finite-element/boundary-element/scaled boundary finite-element model
CN106646645B (zh) 一种重力正演加速方法
CN108875195B (zh) 一种考虑接触的三维力学随机振动仿真模拟方法
CN107479092A (zh) 一种基于方向导数的频率域高阶声波方程正演模拟方法
CN112528546B (zh) 一种非结构化网格的重磁数据三维正反演方法
CN112001100B (zh) 一种三维地震波场se-ibe耦合模拟方法
Golubev et al. Compact grid-characteristic scheme for the acoustic system with the piece-wise constant coefficients
Cai et al. A mixed cover meshless method for elasticity and fracture problems
Li et al. Stability analysis of an explicit integration algorithm with 3D viscoelastic artificial boundary elements
Chen et al. A gradient stable node-based smoothed finite element method for solid mechanics problems
CN111638556A (zh) 基于地空分解策略的大地电磁正演方法及装置、存储介质
Tscherning Functional methods for gravity field approximation
Huang et al. A multi-block finite difference method for seismic wave equation in auxiliary coordinate system with irregular fluid–solid interface
CN105184010A (zh) 基于快速多极间接边界元法的高频地震波散射模拟方法
Giraldo et al. The spectral cell method in nonlinear earthquake modeling
Chen et al. Dynamic modelling of axisymmetric foundations
Yamazaki et al. Three‐dimensional cut‐cell modelling for high‐resolution atmospheric simulations
GB2554865A (en) Seismic modeling
CN102938017B (zh) 基于无网格模型计算周期结构板声学散射系数的方法
Luzón et al. Propagation of SH elastic waves in deep sedimentary basins with an oblique velocity gradient
CN112748471B (zh) 一种非结构化等效源的重磁数据延拓与转换方法

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