简单典型二阶系统_【文献选译】二阶弹性波动方程PML的简单实现

A simple implementation of PML for second-order elastic wave equations

Mingwei Zhuang（厦门大学）

Qiwei Zhan（杜克大学）

Jianyang Zhou（厦门大学）

Zichao Guo（厦门大学）

Na Liu（厦门大学）

Qing Huo Liu（杜克大学）

刊载：2019.08.13

地球物理局 地震波场模拟实验室 边界条件组 译

Zhuang, M., Zhan, Q., Zhou, J., Guo, Z., Liu, N., & Liu, Q. H. (2019). A simple implementation of PML for second-order elastic wave equations. Computer Physics Communications, 106867. doi:10.1016/j.cpc.2019.106867

---

摘 要

在模拟弹性波在无界空间中的传播时，标准完全匹配层（PML）对于一阶偏微分方程（PDEs）而言简单直接；相比之下，PML需要以二阶PDE的形式对控制方程进行大量的重构，但由于内存和时间消耗少得多，这种方法更可取。因此，探索一种二阶系统PML的简单实现势在必行。在本研究中，我们首先系统地将一阶近似PML(NPML)技术推广到二阶系统，用谱元法和时域有限差分算法实现。它具有以下优点：通过保持基于二阶PDE的控制方程基本相同，使得实现简单；通过引入一组辅助常微分方程（ODEs）来提高计算效率。在数学上，这种PML技术有效地混合了二阶PDEs和一阶ODEs，并在局部衰减了输出波，从而有效地避免了空间或时间上的全局卷积。数值实验表明，二阶PDE的NPML在吸收精度、实现复杂度、计算效率等方面对弹性、非弹性和各向异性介质都有很好的吸收性能。

---

1 引 言

---

2 控制方程

在非均匀弹性各向异性介质中，三维笛卡尔坐标下的线性波动方程为：

本构关系为：

应变的定义为：

其中

为位移向量，

表示应力张量，

为应变张量；

为弹性介质的四阶刚度张量，通常用Voigt表示法

表示，

为固体质量密度，而

表示外源力；

符号表示双张量收缩（double tensorial contraction）操作，上角标

表示矩阵转置。

然而，非弹性波的衰减通常由一个被称为品质因子的无量纲参数

来描述，以表征波的能量耗散。这将导致频率域中的弹性张量

成为复模量，其中相速度与弹性模量的实部相关，衰减与弹性模量的虚部相关。在时域中，可以用标准线性固体串联（SLS）来近似黏弹性效应，而弹性张量

可以重写为

^{[1] [2]}：

---

3 二阶波动方程的PML

3.1 标准完全匹配层

我们首先简要回顾一下笛卡尔坐标系中的标准PML。为了得到给定波动方程的PML公式，基于PML区域复坐标伸展的概念引入复坐标。将复坐标中的方程转化为原始的笛卡尔坐标^[3]得到：

其中

是复伸展函数，决定PML的特征。复伸展函数

的一般选择如下

^[3]：

其中

为衰减剖面，导致PML域内传播波场的衰减呈指数下降，

为频移因子，使得衰减取决于频率，因此提供了一个Butterworth型滤波器

^[4]，而

是比例因子，导致PML层内的材料各向异性，以及减少垂直于PML层的相速度

^[3]。

，

为角频率。

对（1）-（3）式进行傅里叶变换，然后利用复坐标拉伸变换空间坐标。由于外源力在PML区域应该是零，所以我们有：

注意指标

对应于

。通过将式（9）代入式（11），我们有

为了得到仅基于位移场的时域PML公式，我们需要对上述公式进行重新表述，将等式（12a）两边同时乘以

，并重新整理，将等式变换回时间域，得到：

其中

且

，

表示卷积积分，

为傅里叶反变换。根据Xie et al.

^[5]的方法，

和

的卷积项可以用递归卷积技术精确计算，但表达式仍然相对复杂（Xie et al.

^[5]式（18））。 因此，对于二阶位移波动方程，标准PML不容易实现。为了将PML技术引入二阶系统，需要对现有的数值方法进行大量的重构。

3.2 近似完全匹配层

根据Cummer^[6]，NPML的想法是将复伸展函数直接从空间偏导数外部偏移到物理场中。由于复伸展函数随空间变化，所以这种变化是不严格正确的，但这种变化并不会显著影响NPML的性能。因此，式（12）需要重写成NPML公式的形式，变为：

接下来，我们将证明NPML与标准PML基本等价。将式（14）两端乘以

，得到：

注意，我们总能找到两个与式(15)右端的空间偏导数无关的伸展坐标函数。然后我们可以重构NPML公式为：

将（16）与（12）比较，我们发现除了NPML区域的所有物理场都乘以三个方向的拉伸坐标函数外，它们的表达式是相似的。差别在于，标准PML是PML区中的

和

的方程，而NPML是

和

的方程。因为伸展坐标函数在物理区域与NPML之间的界面上，它确保了界面两边

和

的连续性。因此，NPML和标准PML公式在这种变化下等价。将方程（14）变换到时间域，并重新整理，我们得到：

其中

通过引入18个辅助变量，式（18）可以在NPML层中简单地求解。

---

【未完待续】

---

《刺客信条：英灵殿》

^[7]

封面图来源^[8]

参考

1. ^ Q. Zhan, M. Zhuang, Q. Sun, Q. Ren, Y. Ren, Y. Mao, Q.H. Liu, IEEE Geosci. Remote. S 55 (10) (2017) 5577–5584.
2. ^ Q. Liu, (Ph.D. thesis), California Institute of Technology, 2006.
3. ^^abcW. Zhang, Y. Shen, Geophysics 75 (4) (2010) T141–T154.
4. ^ R. Matzen, Internat. J. Numer. Methods Engrg. 88 (10) (2011) 951–973.
5. ^^ab Z. Xie, D. Komatitsch, R. Martin, R. Matzen, Geophys. J. Int. 198 (3) (2014) 1714–1747.
6. ^S.A. Cummer, IEEE Microw. Wirel. Co. 13 (3) (2003) 128–130.
7. ^https://www.google.com/url?sa=i&url=https%3A%2F%2Fgamingbolt.com%2Fassassins-creed-valhalla-gets-first-screenshots-details-season-pass-and-collectors-edition&psig=AOvVaw1AqSJeevhkLmyzujRojlnJ&ust=1588490770120000&source=images&cd=vfe&ved=0CAIQjRxqFwoTCOjQ_bTTlOkCFQAAAAAdAAAAABAQ
8. ^https://www.google.com/url?sa=i&url=http%3A%2F%2Fwww.pushsquare.com%2Fnews%2F2020%2F04%2Fof_course_assassins_creed_valhalla_brings_back_the_worst_part_of_the_series&psig=AOvVaw1AqSJeevhkLmyzujRojlnJ&ust=1588490770120000&source=images&cd=vfe&ved=0CAIQjRxqFwoTCOjQ_bTTlOkCFQAAAAAdAAAAABAX