[发明专利]二阶波动方程的非分裂完全匹配层的构造方法

说明书

技术领域

本发明属于数值模拟仿真技术领域，具体为一种适用于二阶波动方程的非分裂完全吸收匹配层的构造方法。

背景技术

基于波动方程的数值仿真中，需采用数值方法离散化波动方程。由于电脑内存和计算时间的限制，需截断计算区域。由此将带来边界反射的问题。通常，在截断边界处可设置吸收边界条件，以消除或减弱由于截断计算区域所引起的边界反射波的影响。完全匹配层（PML）^[1]是目前仿真中运用最广的吸收边界之一。 PML是一种特殊的介质层，它具有这样的性质，当入射波从仿真区域进入其内部区域时，不会发生反射，而在其内部的传播过程中迅速衰减，进而达到消除反射波的目的。由于PML最先是为一阶波动方程而设计，以往对PML的应用多局限于一阶形式的波动方程。对二阶方程的PML研究则较少。一般情况下，一个二阶方程可分解为一阶方程组，然后便可使用PML。然而某些情况下必须使用二阶波动方程，如超声非线性成像模拟中运用广泛的全波模型和Westervelt 模型^[2],[3]，因此找到一种有效的适用于二阶波动方程的PML对此类仿真具有重要意义。另外，二阶形式的方程由于只涉及一个主要变量，在迭代和计算上也具有优势。目前很多商业有限元仿真软件是以二阶方程为核心，二阶形式的方程更便于采用这些软件直接进行仿真。若将二阶方程化为一阶方程组，则会破坏原来方程的结构，增加计算的复杂度。

现有的几种适用于二阶波动方程的PML主要有两种类型。第一种是基于场量分裂^[4]-[6]的方法。该方法将声压分为若干部分，然后分别处理各个部分。其主要缺点为所需的额外存储量较大，计算量大。第二种是一种混合形式（mixed form）的方法^[3],[7]，它形式上是二阶的，但执行以一阶的方式进行。该方法的主要缺点是破坏了原二阶方程的执行结构，处理不便，计算量较大。

本发明克服了已有方法的缺点，提供了一种执行效率高，存储量低的PML, 以解决二阶波动方程仿真中的边界问题。

发明内容

本发明的目的是提出一种执行效率高，存储量低的构造二阶波动方程的完全匹配层（PML）的方法，以解决基于二阶波动方程的仿真中的吸收边界问题。 

本发明提供的适用于二阶波动方程完全匹配层的构造方法，其步骤如下：

1、基于伸缩坐标变换，采用直接求导的方法获得二阶波动方程的PML的频域表达式；

2、以角频率为中心变量，进行因式分解，逐步引入辅助变量，构造辅助微分方程，得到时域表达式。

下面对各步骤进行具体描述。

步骤一、基于伸缩坐标变换，采用直接求导的方法获得二阶波动方程的PML的频域表达式；

二阶声场波动方程的频域表达式为： 

          （1）

以x方向为例，推导二阶波动方程的PML。根据伸缩变换法^[8]，PML的构造可通过如下变换实现：

                （2）

因为σ_x是关于x的函数。故有：

               （3）

通过式（1）和式（3）可直接得到二阶波动方程x方向上的频域表达式为：

  （4）

其中，为伸缩坐标变换算子， 为s_x关于x的导数，ω为角频率，j为虚数单位；为声压u的频率表示；c为声速。 为衰减系数，为其关于x的导数。

用相同的处理方式，即可得 y，z方向的PML频域表达式：

 （5）

（6）

其中分别为y, z方向上的衰减参数，分别为其关于y,z的导数。

步骤二、以角频率为中心变量，进行因式分解，逐步引入辅助变量，构造辅助微分方程，得到时域表达式。 

对（4）式进行因式分解可得：

 （7）

引入辅助变量u₁ ，将式（7）式写为如下两式： 

                    （8）

（9）

式（7）两边同时乘以jω+σ_x，可得：