CN113221409B - 一种有限元和边界元耦合的声波二维数值模拟方法和装置 - Google Patents

Abstract

本申请涉及一种有限元和边界元耦合的声波二维数值模拟方法、装置、计算机设备和存储介质。所述方法包括：通过获取待模拟区域信息，根据待研究区域信息和预设的频率域二维标量声波方程，对非均匀介质及复杂地形区域通过Galerkin有限元法建立有限元方程，对空气介质、地下均匀介质与有限元研究区域的边界区域通过格林公式建立边界元方程；根据有限元方程和边界元方程建立耦合方程组；求解耦合方程组，得到待研究区域声波波场的数值模拟值。本发明所提出的有限元和边界元耦合方结合了两种方法的优势，规避了两种方法的缺点，能够很好解决复杂波场的正演模拟问题。

Description

一种有限元和边界元耦合的声波二维数值模拟方法和装置

技术领域

本申请涉及计算机数值模拟领域，特别是涉及一种有限元和边界元耦合的声波二维数值模拟方法、装置、计算机设备和存储介质。

背景技术

地震勘探是一种应用非常广泛的探测方法，它主要是以地下介质的密度和速度差异为物质基础，以人工方式激发地震波，并采用检波器接收各个观测点的振动信号，对采集的资料进行处理、解释，推断地下地层岩性及含油气可能性。与其他物探方法相比，地震勘探具有分层详细准确、分辨率高、探测深度大等优点，并且已广泛应用于石油、天然气、固体资源勘探等领域。

目前，对于大多数地震问题，由于研究区域的复杂性，一般难以推导其解析表达式，因此，地震数值模拟方法显得尤为重要。在众多数值模拟方法中，有限单元法是将求解微分方程运用数学工具转变成求泛函的相对应极值问题，然后将研究区域进行网格剖分，每个单元采用插值函数，将各单元节点的函数扩展成整个研究区域的线性方程组，最后求解方程组。有限单元法适用性好，网格剖分灵活，能够解决非均匀介质和复杂地形数值模拟问题。然而，在实际波场正演模拟中，受硬件设备限制，只能模拟有限区域内波的传播。如果不设边界条件，边界反射会对模拟区域的波场产生强烈的干扰，从而影响数值模拟精度。为解决边界反射对研究区域波场模拟的影响，一般采用扩展边界、一次性边界和吸收边界条件。扩展边界增大研究范围才能取得好的效果，但随着研究区域的扩大，计算时间增加，影响求解的收敛速度，同时还会降低数值模拟精度，因此，该方法在地震中很少用到。Berenger提出了一种PML吸收边界条件，在电磁波中取得很好的效果，后来人们将该方法应用到地震边界处理中。以上边界处理方法改变了波场的动力学特性，极大地影响纵波速度和密度的反演成像效果。

边界元单元法根据偏微分方程所满足边值问题的通过数学工具转化为它所对应的边界型的积分方程，再对边界离散、插值、积分后相加，最后求解线性方程组，得到各节点函数值。与其他数值方法相比，边界单元法只对求解区域的边界进行剖分，降低了求解问题的维数；基本解能够自动满足无穷远处的边界条件，因此边界元适合处理无限域和半无限域问题；仅在边界剖分，能够很好处理不规则界面和复杂起伏地形问题。但由于边界单元法只在边界剖分，因此，对于非均匀介质求解比较繁琐；另外，边界单元法组成的系数矩阵是不对称满秩矩阵，当剖分网格节点较多时，计算效率低；对于复杂问题，基本解难以推导；编写程序比较复杂。

因此，现有技术存在计算效率低、效果不佳的问题。

发明内容

基于此，有必要针对上述技术问题，提供一种能够提高声波数值模拟效率的有限元和边界元耦合的声波二维数值模拟方法、装置、计算机设备和存储介质。

一种有限元和边界元耦合的声波二维数值模拟方法，所述方法包括：

获取待模拟区域信息，根据所述待研究区域信息和预设的频率域二维标量声波方程，对第一研究区域通过Galerkin有限元法建立有限元方程，对第二研究区域通过格林公式建立边界元方程；所述第一研究区域为非均匀介质及复杂地形区域，所述第二研究区域为空气介质与所述第一研究区域的边界区域，以及地下均匀介质与所述第一研究区域的边界区域；

根据所述有限元方程和所述边界元方程建立耦合方程组；

求解所述耦合方程组，得到所述待研究区域声波波场的数值模拟值。

在其中一个实施例中，还包括：根据所述待研究区域信息和预设的频率域二维标量声波方程，将频率域速度位满足的边值问题采用Galerkin法转换为有限元方程问题；

将所述第一研究区域剖分成多个四边形单元，再将所述四边形单元剖分成两个三角形，将所述三角形作为基单元；

采用线性插值算法求解所述基单元的第一单元积分系数；

将每个基单元的所述第一单元积分系数通过压缩存储方式合成稀疏矩阵，建立有限元方程。

在其中一个实施例中，还包括：根据所述待研究区域信息和预设的频率域二维标量声波方程，对第一研究区域通过Galerkin有限元法建立有限元方程为：

Ku＝f

其中，K为组合成的Na×Na维的刚度矩阵，Na为所述第一研究区域剖分的网格节点总数，Na＝N_x×N_z，N_x为所述第一研究区域水平方向剖分的节点个数，N_z为所述第一研究区域垂直方向剖分的节点个数；u为Na×1维的矩阵，即Na行的列矩阵，f为Na×1维的矩阵，为右端项。

在其中一个实施例中，还包括：根据所述待研究区域信息和预设的频率域二维标量声波方程，对第二研究区域通过格林公式建立边界元方程；所述第二研究区域包括空气介质与所述第一研究区域的第一边界区域，以及地下均匀介质与所述第一研究区域的第二边界区域。

在其中一个实施例中，还包括：根据所述待研究区域信息和预设的频率域二维标量声波方程，将频率域速度位满足的二维边值问题采用格林公式转换为边界上一维积分方程问题；

对所述第二研究区域进行网格剖分；

采用高斯积分方式计算第二单元积分系数；

将所述第二单元积分系数通过压缩存储方式合成稠密矩阵，建立边界元方程。

在其中一个实施例中，还包括：根据所述待研究区域信息和预设的频率域二维标量声波方程，对第二研究区域通过格林公式建立边界元方程，得到所述第一边界区域的边界元方程为：

其中，ω＝diag(ω_j)，u＝(u_i)^T，

C＝(C_i)^T，F＝(F_ji)^T，D＝(D_ji)^T，i,j＝1,…n，n即为所述第一边界区域上剖分的节点总数；C、F、D表示积分系数；

得到所述第二边界区域的边界元方程为：

其中，ω′＝diag(ω_j′)，u＝(u_i)^T，

Α＝(A_i)^T，G＝(G_ji)^T，E＝(E_ji)^T，i,j＝1,…m，m即为所述第二边界区域上剖分的节点总数；Α、G、H表示积分系数。

在其中一个实施例中，还包括：通过求解器PARDISO求解所述耦合方程组，得到所述待研究区域声波波场的数值模拟值。

一种有限元和边界元耦合的声波二维数值模拟装置，所述装置包括：

有限元方程和边界元方程建立模块，用于获取待模拟区域信息，根据所述待研究区域信息和预设的频率域二维标量声波方程，对第一研究区域通过Galerkin有限元法建立有限元方程，对第二研究区域通过格林公式建立边界元方程；所述第一研究区域为非均匀介质及复杂地形区域，所述第二研究区域为空气介质与所述第一研究区域的边界区域，以及地下均匀介质与所述第一研究区域的边界区域；

耦合方程组建立模块，用于根据所述有限元方程和所述边界元方程建立耦合方程组；

方程组求解模块，用于求解所述耦合方程组，得到所述待研究区域声波波场的数值模拟值。

一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现以下步骤：

获取待模拟区域信息，根据所述待研究区域信息和预设的频率域二维标量声波方程，对第一研究区域通过Galerkin有限元法建立有限元方程，对第二研究区域通过格林公式建立边界元方程；所述第一研究区域为非均匀介质及复杂地形区域，所述第二研究区域为空气介质与所述第一研究区域的边界区域，以及地下均匀介质与所述第一研究区域的边界区域；

根据所述有限元方程和所述边界元方程建立耦合方程组；

求解所述耦合方程组，得到所述待研究区域声波波场的数值模拟值。

一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：

获取待模拟区域信息，根据所述待研究区域信息和预设的频率域二维标量声波方程，对第一研究区域通过Galerkin有限元法建立有限元方程，对第二研究区域通过格林公式建立边界元方程；所述第一研究区域为非均匀介质及复杂地形区域，所述第二研究区域为空气介质与所述第一研究区域的边界区域，以及地下均匀介质与所述第一研究区域的边界区域；

根据所述有限元方程和所述边界元方程建立耦合方程组；

求解所述耦合方程组，得到所述待研究区域声波波场的数值模拟值。

上述有限元和边界元耦合的声波二维数值模拟方法、装置、计算机设备和存储介质，通过获取待模拟区域信息，根据待研究区域信息和预设的频率域二维标量声波方程，对非均匀介质及复杂地形区域通过Galerkin有限元法建立有限元方程，对空气介质、地下均匀介质与有限元研究区域的边界区域通过格林公式建立边界元方程；根据有限元方程和边界元方程建立耦合方程组；求解耦合方程组，得到待研究区域声波波场的数值模拟值。本发明所提出的有限元和边界元耦合方结合了有限单元法网格剖分灵活，适用于复杂物性问题和边界单元法基本解满足无穷远边界条件，适合处理边界问题的优势，规避两种方法的缺点，能够很好解决复杂波场的正演模拟问题。

附图说明

图1为一个实施例中有限元和边界元耦合的声波二维数值模拟方法的流程示意图；

图2为一个具体实施例中有限元和边界元耦合的声波二维数值模拟方法的流程示意图；

图3为一个具体实施例中有限元和边界元耦合模型示意图；

图4为一个具体实施例中有限元和边界元耦合边界示意图，其中4a为边界ABΓ₁CD示意图，4b为边界ABΓCD示意图；

图5为一个具体实施例中均匀全空间解析解、有限元、有限元与边界元耦合方法振幅等值线图；

图6为一个具体实施例中不同深度解析解、有限元、有限元与边界元耦合方法振幅及相对误差图，其中6a为深度为0m时解析解、有限元、有限元与边界元耦合方法振幅示意图，6b为深度为-100m时解析解、有限元、有限元与边界元耦合方法振幅示意图，6c为深度为0m时解析解、有限元、有限元与边界元耦合方法相对误差图，6d为深度为-100m时解析解、有限元、有限元与边界元耦合方法相对误差图；

图7为一个实施例中有限元和边界元耦合的声波二维数值模拟装置的结构框图；

图8为一个实施例中计算机设备的内部结构图。

具体实施方式

为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。

本申请提供的有限元和边界元耦合的声波二维数值模拟方法，可以应用于如下应用环境中。其中，终端执行一种有限元和边界元耦合的声波二维数值模拟方法。通过获取待模拟区域信息，根据待研究区域信息和预设的频率域二维标量声波方程，对非均匀介质及复杂地形区域通过Galerkin有限元法建立有限元方程，对空气介质、地下均匀介质与有限元研究区域的边界区域通过格林公式建立边界元方程；根据有限元方程和边界元方程建立耦合方程组；求解耦合方程组，得到待研究区域声波波场的数值模拟值。。其中，终端可以但不限于是各种个人计算机、笔记本电脑、平板电脑和便携式设备。

在一个实施例中，如图1所示，提供了一种有限元和边界元耦合的声波二维数值模拟方法，包括以下步骤：

步骤102，获取待模拟区域信息，根据待研究区域信息和预设的频率域二维标量声波方程，对第一研究区域通过Galerkin有限元法建立有限元方程，对第二研究区域通过格林公式建立边界元方程。

在数值模拟方法中，有限单元法是将求解微分方程运用数学工具转变成求泛函的相对应极值问题，然后将研究区域进行网格剖分，每个单元采用插值函数，将各单元节点的函数扩展成整个研究区域的线性方程组，最后求解方程组。有限单元法适用性好，网格剖分灵活，能够解决非均匀介质和复杂地形数值模拟问题，但边界条件影响正演效果。边界元单元法根据偏微分方程所满足边值问题的通过数学工具转化为它所对应的边界型的积分方程，再对边界离散、插值、积分后相加，最后求解线性方程组，得到各节点函数值。与其他数值方法相比，边界单元法只对求解区域的边界进行剖分，降低了求解问题的维数；基本解能够自动满足无穷远处的边界条件，适合处理无限域和半无限域问题；仅在边界剖分，能够很好处理不规则界面和复杂起伏地形问题。但由于边界单元法只在边界剖分，因此，对于非均匀介质求解比较繁琐。

本发明以非均匀介质及复杂地形区域为第一研究区域，采用有限单元法，以空气介质、地下均匀介质与非均匀介质及复杂地形区域的边界区域为第二研究区域，采用边界单元法。

步骤104，根据有限元方程和边界元方程建立耦合方程组。

将第二研究区域作为第一研究区域的一个特殊单元加入到有限元方程中，根据有限元方程和边界元方程建立耦合方程组。

步骤106，求解耦合方程组，得到待研究区域声波波场的数值模拟值。

耦合方程组是一个稀疏矩阵，它的带宽为耦合边界上剖分节点个数，并将节点上对应的系数加入到有限元方程中。采用求解器PARDISO求得各个节点的波场值。

上述有限元和边界元耦合的声波二维数值模拟方法中，通过获取待模拟区域信息，根据待研究区域信息和预设的频率域二维标量声波方程，对非均匀介质及复杂地形区域通过Galerkin有限元法建立有限元方程，对空气介质、地下均匀介质与有限元研究区域的边界区域通过格林公式建立边界元方程；根据有限元方程和边界元方程建立耦合方程组；求解耦合方程组，得到待研究区域声波波场的数值模拟值。本发明所提出的有限元和边界元耦合方结合了有限单元法网格剖分灵活，适用于复杂物性问题和边界单元法基本解满足无穷远边界条件，适合处理边界问题的优势，规避两种方法的缺点，能够很好解决复杂波场的正演模拟问题。

在其中一个实施例中，还包括：根据待研究区域信息和预设的频率域二维标量声波方程，将频率域速度位满足的边值问题采用Galerkin法转换为有限元方程问题；将第一研究区域剖分成多个四边形单元，再将四边形单元剖分成两个三角形，将三角形作为基单元；采用线性插值算法求解基单元的第一单元积分系数；将每个基单元的第一单元积分系数通过压缩存储方式合成稀疏矩阵，建立有限元方程为：

Ku＝f

其中，K为组合成的Na×Na维的刚度矩阵，Na为第一研究区域剖分的网格节点总数，Na＝N_x×N_z，N_x为第一研究区域水平方向剖分的节点个数，N_z为第一研究区域垂直方向剖分的节点个数；u为Na×1维的矩阵，即Na行的列矩阵，f为Na×1维的矩阵，为右端项。

在其中一个实施例中，还包括：根据待研究区域信息和预设的频率域二维标量声波方程，对第二研究区域通过格林公式建立边界元方程；第二研究区域包括空气介质与第一研究区域的第一边界区域，以及地下均匀介质与第一研究区域的第二边界区域。

在其中一个实施例中，还包括：根据待研究区域信息和预设的频率域二维标量声波方程，将频率域速度位满足的二维边值问题采用格林公式转换为边界上一维积分方程问题；对第二研究区域进行网格剖分；采用高斯积分方式计算第二单元积分系数；将第二单元积分系数通过压缩存储方式合成稠密矩阵，建立边界元方程。得到第一边界区域的边界元方程为：

其中，ω＝diag(ω_j)，u＝(u_i)^T，

C＝(C_i)^T，F＝(F_ji)^T，D＝(D_ji)^T，i,j＝1,…n，n即为第一边界区域上剖分的节点总数；C、F、D表示积分系数；

得到第二边界区域的边界元方程为：

其中，ω′＝diag(ω_j′)，u＝(u_i)^T，

Α＝(A_i)^T，G＝(G_ji)^T，H＝(H_ji)^T，i,j＝1,…m，m即为第二边界区域上剖分的节点总数；Α、G、H表示积分系数。

在其中一个实施例中，还包括：通过求解器PARDISO求解耦合方程组，得到待研究区域声波波场的数值模拟值。

应该理解的是，虽然图1的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。

在一个具体实施例中，如图2所示，提供了一种有限元和边界元耦合的声波二维数值模拟方法，包括以下步骤：

第一步.根据频率域二维标量声波方程，利用Galerkin有限元法建立相应的有限元方程；

其中，所述建立的有限元方程包括：

(1)将频率域速度位满足的边值问题采用Galerkin法转换为有限元方程问题；

(2)区域剖分，对于非均匀介质及复杂地形区域采用有限单元法见图3，二维模型首先剖分成多个四边形单元，然后将每个四边形剖分成两个三角形，三角形作为基单元；

(3)单元分析，采用线性插值的方式求解单元积分系数。

(4)总体合成，每个单元计算的系数矩阵采用压缩存储方式合成一个大的稀疏矩阵中；

则形成的方程组为：

Ku＝f                            (1)

式中，K为组合成的Na×Na维的刚度矩阵，Na为所述第一研究区域剖分的网格节点总数，Na＝N_x×N_z，N_x为所述第一研究区域水平方向剖分的节点个数，N_z为所述第一研究区域垂直方向剖分的节点个数；u为Na×1维的矩阵，即Na行的列矩阵，f为Na×1维的矩阵，为右端项。

第二步.根据频率域速度位满足的边界问题，利用格林公式建立相应地边界元方程；

其中，所述建立的边界元方程包括：

(1)将频率域速度位满足的二维边值问题采用格林公式转换为边界上一维积分方程问题；

(2)网格剖分，对于图3区域Ω₀空气介质和Ω₂均匀介质采用边界元法，只在研究区域边界上进行剖分网格，即图3区域Ω₀空气介质和Ω₂的边界上进行网格剖分；

(3)数值积分，采用高斯积分方式计算单元积分系数。

(4)合成线性方程组，每个单元计算的系数矩阵合成一个大的稠密矩阵中；

对于空气层与有限元区域耦合的边界ABΓ₁CD，如图4a虚线部分，可以得到边界方程可以表示为：

式中，

其中，ω_j是p点对区域Ω的张角。Ω₀中的基本解为

K₀是第二类零阶修正贝塞尔函数，r是源点到观测点之间的距离，k₀、k₁分别为空气和地下均匀介质中的波数，v₀、v₁分别为空气和地下介质的速度，ρ₀、ρ₁分别为空气和地下介质的密度，i为虚数单位，K₁是第二类一阶修正贝塞尔函数。

对每一个节点都可以得到如式(3)的方程，那么全部节点得到的线性方程组为，

式中，ω＝diag(ω_j)，u＝(u_i)^T，

C＝(C_i)^T，F＝(F_ji)^T，D＝(D_ji)^T，i,j＝1,…n，n即为边界ABΓ₁CD上剖分的节点总数；C、F、D表示积分系数；。

同理，对于有限元地下边界处理ABΓCD，如图4b虚线部分，同样采用高斯积分的方式，可得到类似的方程组：

其中，

全部节点得到的线性方程组为：

式中，ω′＝diag(ω_j′)，u＝(u_i)^T，

Α＝(A_i)^T，G＝(G_ji)^T，H＝(H_ji)^T，i,j＝1,…m，m即为边界ABΓCD上剖分的节点总数；Α、G、H表示积分系数。

第三步.耦合方程组

对分别构建的有限元和边界元方程，将边界元研究区域看作有限元域的一个特殊单元，加入到有限元方程组中，建立耦合方程组；

具体步骤如下：

对于有限元域建立的方程(1)，可以改写为，

[K_Na×Na][u_Na×1]＝[f_Na×1]                       (8)

对于Ω₀空气层建立的边界方程(4)，改写为，

将F和D分成三个部分AB为F_l、D_l，CD为F_r、D_r，

和

为耦合边界Γ₁。式(9)可以写为，

同理，对于Ω₂域建立的边界方程式(10)，改写为，

将G和H分成三个部分AB为G_d、H_d，CD为G_n、H_n，G_Γ和H_Γ为耦合边界Γ上的波场值及波场一阶导数。

则，式(11)可以表示为，

将式(8)、(10)、(12)方程组进行耦合，可直观地表示为，

式中：

b₁＝[C_l A_l]^T，b₂＝[C_r A_r]^T，b₂＝[C_r A_r]^T,BEM表示边界元方法，FEM表示有限元方法。

第四步.方程组求解

本发明有限元和边界元耦合矩阵，仍然是一个稀疏矩阵，它的带宽为耦合边界上剖分节点个数，并将节点上对应的系数加入到有限元方程中。最终组合的方程(9)采用直接求解器PARDISO求得各个节点的波场值。

在另一个具体实施例中，设计了一个有解析解的均匀全空间模型，研究区域大小为1000m×1000m，剖分网格大小为101×101；水平和垂直方向采样间隔均10m；背景速度为3000m/s；震源为主频30Hz的子波，震源坐标为(-500m,0m)，计算频率为10Hz。

图5为解析解、有限元和边界元耦合方法及有限元计算的均匀全空间振幅谱结果，从图5可以看出解析解与耦合方法振幅谱形态近似，并且无论是在研究区域中心还是边界，有限元和边界耦合方法与解析解均匀吻合较好，但是有限单元法与解析解在研究区域中心较好，随着向四周边界的逼近，解析解与有限元法吻合度越来越差，从而证明了有限元中采用的衰减边界等边界条件改变了波场的动力学特性，而采用有限元与边界元耦合方法处理边界，是一种严格的边界条件，在理论上是完备的。

为了更加清楚的研究在某一深度上波场信息，取不同深度测线振幅谱及相对误差，如图6所示，从图中可以看出不同深度解析解与有限元和边界元耦合算法振幅相对误差均远远小于1％，证明有限元和边界元耦合方法是正确的且精度较高。在有限元声波数值模拟中算例中已经验证了衰减边界会影响波场的传播规律，从图6a振幅谱中可以看出有限元法模拟结果边界处波场振幅与解析解出现偏差，且误差较大，取地面z＝0m和z＝-100m测线振幅谱相对误差可以看出，有限元与解析解振幅谱在边界处误差较大，在其他位置吻合较好，而有限元和边界元耦合方法计算的振幅谱与解析解吻合很好，误差较小。因此，有限元和边界元耦合方法对于边界元处理更有效，既能够抑制边界截断效应问题，并且不改变波的动力学特征。

在一个实施例中，如图7所示，提供了一种有限元和边界元耦合的声波二维数值模拟装置，包括：有限元方程和边界元方程建立模块702、耦合方程组建立模块704和方程组求解模块706，其中：

有限元方程和边界元方程建立模块702，用于获取待模拟区域信息，根据待研究区域信息和预设的频率域二维标量声波方程，对第一研究区域通过Galerkin有限元法建立有限元方程，对第二研究区域通过格林公式建立边界元方程；第一研究区域为非均匀介质及复杂地形区域，第二研究区域为空气介质与第一研究区域的边界区域，以及地下均匀介质与第一研究区域的边界区域；

耦合方程组建立模块704，用于根据有限元方程和边界元方程建立耦合方程组；

方程组求解模块706，用于求解耦合方程组，得到待研究区域声波波场的数值模拟值。

有限元方程和边界元方程建立模块702还用于将第一研究区域剖分成多个四边形单元，再将四边形单元剖分成两个三角形，将三角形作为基单元；采用线性插值算法求解基单元的第一单元积分系数；将每个基单元的第一单元积分系数通过压缩存储方式合成稀疏矩阵，建立有限元方程。

有限元方程和边界元方程建立模块702还用于根据待研究区域信息和预设的频率域二维标量声波方程，将频率域速度位满足的二维边值问题采用格林公式转换为边界上一维积分方程问题；对第二研究区域进行网格剖分；采用高斯积分方式计算第二单元积分系数；将第二单元积分系数通过压缩存储方式合成稠密矩阵，建立边界元方程。

方程组求解模块706还用于通过求解器PARDISO求解耦合方程组，得到待研究区域声波波场的数值模拟值。

关于有限元和边界元耦合的声波二维数值模拟装置的具体限定可以参见上文中对于有限元和边界元耦合的声波二维数值模拟方法的限定，在此不再赘述。上述有限元和边界元耦合的声波二维数值模拟装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。

在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图8所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种有限元和边界元耦合的声波二维数值模拟方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。

本领域技术人员可以理解，图8中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。

在一个实施例中，提供了一种计算机设备，包括存储器和处理器，该存储器存储有计算机程序，该处理器执行计算机程序时实现上述方法实施例中的步骤。

在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法实施例中的步骤。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(ROM)、可编程ROM(PROM)、电可编程ROM(EPROM)、电可擦除可编程ROM(EEPROM)或闪存。易失性存储器可包括随机存取存储器(RAM)或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM(SRAM)、动态RAM(DRAM)、同步DRAM(SDRAM)、双数据率SDRAM(DDRSDRAM)、增强型SDRAM(ESDRAM)、同步链路(Synchlink)DRAM(SLDRAM)、存储器总线(Rambus)直接RAM(RDRAM)、直接存储器总线动态RAM(DRDRAM)、以及存储器总线动态RAM(RDRAM)等。

以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。

Claims (8)

1.一种有限元和边界元耦合的声波二维数值模拟方法，其特征在于，所述方法包括：

获取待研究区域信息，根据所述待研究区域信息和预设的频率域二维标量声波方程，对第一研究区域通过Galerkin有限元法建立有限元方程，对第二研究区域通过格林公式建立边界元方程；所述第一研究区域为非均匀介质及复杂地形区域，所述第二研究区域包括空气介质与所述第一研究区域的第一边界区域，以及地下均匀介质与所述第一研究区域的第二边界区域；其中，所述边界元方程包括第一边界区域的边界元方程和第二边界区域的边界元方程；

所述第一边界区域的边界元方程为：

其中，ω＝diag(ω_j)，u＝(u_i)^T，u为Na×1维的矩阵，即Na行的列矩阵，Na为所述第一研究区域剖分的网格节点总数，Na＝N_x×N_z，N_x为所述第一研究区域水平方向剖分的节点个数，N_z为所述第一研究区域垂直方向剖分的节点个数，

C＝(C_i)^T，F＝(F_ji)^T，D＝(D_ji)^T，i,j＝1,…n1，n1即为所述第一边界区域上剖分的节点总数；C、F、D表示积分系数；

所述第二边界区域的边界元方程为：

其中，ω′＝diag(ω_j′)，u＝(u_i)^T，

Α＝(A_i)^T，G＝(G_ji)^T，H＝(H_ji)^T，i,j＝1,…m，m即为所述第二边界区域上剖分的节点总数；Α、G、H表示积分系数；

根据所述有限元方程和所述边界元方程建立耦合方程组；

求解所述耦合方程组，得到所述待研究区域声波波场的数值模拟值。

2.根据权利要求1所述的方法，其特征在于，根据所述待研究区域信息和预设的频率域二维标量声波方程，对第一研究区域通过Galerkin有限元法建立有限元方程，包括：

根据所述待研究区域信息和预设的频率域二维标量声波方程，将频率域速度位满足的边值问题采用Galerkin法转换为有限元方程问题；

将所述第一研究区域剖分成多个四边形单元，再将所述四边形单元剖分成两个三角形，将所述三角形作为基单元；

采用线性插值算法求解所述基单元的第一单元积分系数；

将每个基单元的所述第一单元积分系数通过压缩存储方式合成稀疏矩阵，建立有限元方程。

3.根据权利要求2所述的方法，其特征在于，根据所述待研究区域信息和预设的频率域二维标量声波方程，对第一研究区域通过Galerkin有限元法建立有限元方程，包括：

根据所述待研究区域信息和预设的频率域二维标量声波方程，对第一研究区域通过Galerkin有限元法建立有限元方程为：

Ku＝f

其中，K为组合成的Na×Na维的刚度矩阵，f为Na×1维的矩阵，为右端项。

4.根据权利要求3所述的方法，其特征在于，根据所述待研究区域信息和预设的频率域二维标量声波方程，对第二研究区域通过格林公式建立边界元方程，包括：

根据所述待研究区域信息和预设的频率域二维标量声波方程，将频率域速度位满足的二维边值问题采用格林公式转换为边界上一维积分方程问题；

对所述第二研究区域进行网格剖分；

采用高斯积分方式计算第二单元积分系数；

将所述第二单元积分系数通过压缩存储方式合成稠密矩阵，建立边界元方程。

5.根据权利要求1至4任意一项所述的方法，其特征在于，求解所述耦合方程组，得到所述待研究区域声波波场的数值模拟值，包括：

通过求解器PARDISO求解所述耦合方程组，得到所述待研究区域声波波场的数值模拟值。

6.一种有限元和边界元耦合的声波二维数值模拟装置，其特征在于，所述装置包括：

有限元方程和边界元方程建立模块，用于获取待研究区域信息，根据所述待研究区域信息和预设的频率域二维标量声波方程，对第一研究区域通过Galerkin有限元法建立有限元方程，对第二研究区域通过格林公式建立边界元方程；所述第一研究区域为非均匀介质及复杂地形区域，所述第二研究区域包括空气介质与所述第一研究区域的第一边界区域，以及地下均匀介质与所述第一研究区域的第二边界区域；其中，所述边界元方程包括第一边界区域的边界元方程和第二边界区域的边界元方程；

所述第一边界区域的边界元方程为：

其中，ω＝diag(ω_j)，u＝(u_i)^T，

C＝(C_i)^T，F＝(F_ji)^T，D＝(D_ji)^T，i,j＝1,…n1，n1即为所述第一边界区域上剖分的节点总数；C、F、D表示积分系数；

所述第二边界区域的边界元方程为：

其中，ω′＝diag(ω_j′)，u＝(u_i)^T，

Α＝(A_i)^T，G＝(G_ji)^T，H＝(H_ji)^T，i,j＝1,…m，m即为所述第二边界区域上剖分的节点总数；Α、G、H表示积分系数；

耦合方程组建立模块，用于根据所述有限元方程和所述边界元方程建立耦合方程组；

方程组求解模块，用于求解所述耦合方程组，得到所述待研究区域声波波场的数值模拟值。

7.根据权利要求6所述的装置，其特征在于，有限元方程和边界元方程建立模块还包括：

根据所述待研究区域信息和预设的频率域二维标量声波方程，将频率域速度位满足的边值问题采用Galerkin法转换为有限元方程问题；

将所述第一研究区域剖分成多个四边形单元，再将所述四边形单元剖分成两个三角形，将所述三角形作为基单元；

采用线性插值算法求解所述基单元的第一单元积分系数；

将每个基单元的所述第一单元积分系数通过压缩存储方式合成稀疏矩阵，建立有限元方程。

8.根据权利要求7所述的装置，其特征在于，有限元方程和边界元方程建立模块还包括：

根据所述待研究区域信息和预设的频率域二维标量声波方程，将频率域速度位满足的二维边值问题采用格林公式转换为边界上一维积分方程问题；

对所述第二研究区域进行网格剖分；

采用高斯积分方式计算第二单元积分系数；

将所述第二单元积分系数通过压缩存储方式合成稠密矩阵，建立边界元方程。

Images