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

知乎用户7QViRQ

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 控制方程

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

\rho \partial_{t t} \mathbf{u}=\nabla \cdot \boldsymbol{\tau}+\mathbf{f}，(1)

本构关系为：

\boldsymbol{\tau}=\mathbf{C}: \boldsymbol{\varepsilon}，(2)

应变的定义为：

\varepsilon=\frac{1}{2}\left[\nabla \mathbf{u}+(\nabla \mathbf{u})^{T}\right]，(3)

其中 \mathbf{u} 为位移向量， \boldsymbol{\tau} 表示应力张量， \boldsymbol{\varepsilon} 为应变张量； \mathbf{C} 为弹性介质的四阶刚度张量，通常用Voigt表示法 C_{i j}(i, j=1,2, \ldots, 6) 表示， \rho 为固体质量密度，而 \mathbf{f} 表示外源力； : 符号表示双张量收缩（double tensorial contraction）操作，上角标 T 表示矩阵转置。

然而，非弹性波的衰减通常由一个被称为品质因子的无量纲参数 Q_{i j} 来描述，以表征波的能量耗散。这将导致频率域中的弹性张量 C_{i j} 成为复模量，其中相速度与弹性模量的实部相关，衰减与弹性模量的虚部相关。在时域中，可以用标准线性固体串联（SLS）来近似黏弹性效应，而弹性张量 C_{i j} 可以重写为^[1]^[2]：

3 二阶波动方程的PML

3.1 标准完全匹配层

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

\frac{\partial}{\partial_{\tilde{k}}}=\frac{1}{s_{k}} \frac{\partial}{\partial_{k}} \quad(k=x, y, z)，(9)

其中 s_{k} 是复伸展函数，决定PML的特征。复伸展函数 s_{k} 的一般选择如下^[3]：

s_{k}=\beta_{k}(k)+\frac{d_{k}(k)}{\alpha_{k}(k)+i \omega}，(10)

其中 d_{k}(k) \geqslant 0 为衰减剖面，导致PML域内传播波场的衰减呈指数下降， \alpha_{k}(k) \geqslant 0 为频移因子，使得衰减取决于频率，因此提供了一个Butterworth型滤波器^[4]，而 \beta_{k}(k) \geqslant 1 是比例因子，导致PML层内的材料各向异性，以及减少垂直于PML层的相速度^[3]。 i^{2}=-1 ， \omega 为角频率。

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

\begin{aligned} &\rho \omega^{2} \tilde{u}_{i}=\sum_{j=1}^{3} \frac{\partial \tilde{\tau}_{i j}}{\partial \tilde{x}_{j}}\\ &\tilde{\tau}_{i j}=\sum_{k, l=1}^{3} C_{i j k l} \frac{\partial \tilde{u}_{k}}{\partial \tilde{x}_{l}} \end{aligned}，(11)

注意指标 \{i, j, k, l\} 对应于 \{x, y, z\} 。通过将式（9）代入式（11），我们有

\bbox[#EFF,5px,border:2px solid blue]{\rho \omega^{2} \tilde{u}_{i}=\sum_{j=1}^{3} \frac{1}{s_{j}} \frac{\partial \tilde{\tau}_{i j}}{\partial x_{j}}}，(12a)

\bbox[#EFF,5px,border:2px solid blue]{\tilde{\tau}_{i j}=\sum_{k, l=1}^{3} C_{i j k l} \frac{1}{s_{l}} \frac{\partial \tilde{u}_{k}}{\partial x_{l}}}，(12b)

为了得到仅基于位移场的时域PML公式，我们需要对上述公式进行重新表述，将等式（12a）两边同时乘以 s_{x} s_{y} s_{z} ，并重新整理，将等式变换回时间域，得到：

\begin{aligned} &\rho \omega^{2} s_{x} s_{y} s_{z} \tilde{u}_{i}=\sum_{j=1}^{3} \frac{\partial \tilde{\tau}_{i j}}{\partial x_{j}} \Longrightarrow \rho L(t) * u_{i}=\sum_{j=1}^{3} \frac{\partial \tau_{i j}}{\partial x_{j}}\\ &\tilde{\tau}_{i j}=\sum_{k, l=1}^{3} C_{i j k l} \frac{s_{x} s_{y} s_{z}}{s_{l} s_{k}} \frac{\partial \tilde{u}_{k}}{\partial x_{l}} \Longrightarrow \tau_{i j}=\sum_{k, l=1}^{3} \bar{c}_{i j k l} * \frac{\partial u_{k}}{\partial x_{l}} \end{aligned}，(13)

其中 L(t)=\mathcal{F}^{-1}\left[\omega^{2} s_{x} s_{y} s_{z}\right] 且 \bar{C}_{i j k l}=C_{i j k l} \mathcal{F}^{-1}\left[\frac{s_{x} s_{y} s_{z}}{s_{l} s_{k}}\right] ， \ast 表示卷积积分， \mathcal{F}^{-1} 为傅里叶反变换。根据Xie et al.^[5]的方法， L(t) 和 \bar{C}_{i j k l} 的卷积项可以用递归卷积技术精确计算，但表达式仍然相对复杂（Xie et al.^[5]式（18））。因此，对于二阶位移波动方程，标准PML不容易实现。为了将PML技术引入二阶系统，需要对现有的数值方法进行大量的重构。

3.2 近似完全匹配层

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

\bbox[#EFF,5px,border:2px solid blue]{\begin{aligned} &\rho \omega^{2} \tilde{u}_{i}=\sum_{j=1}^{3} \frac{\partial \frac{1}{s_{j}} \tilde{\tau}_{i j}}{\partial x_{j}}\\ &\tilde{\tau}_{i j}=\sum_{k, l=1}^{3} C_{i j k l} \frac{\partial \frac{1}{s_{l}} \tilde{u}_{k}}{\partial x_{l}} \end{aligned}}，(14)

接下来，我们将证明NPML与标准PML基本等价。将式（14）两端乘以 \frac{1}{s_{x} s_{y} s_{z}} ，得到：

\begin{aligned} &\frac{1}{s_{x} s_{y} s_{z}} \rho \omega^{2} \tilde{u}_{i}=\sum_{j=1}^{3} \frac{1}{s_{x} s_{y} s_{z}} \frac{\partial \frac{1}{s_{j}} \tilde{\tau}_{i j}}{\partial x_{j}}\\ &\frac{1}{s_{x} s_{y} s_{z}} \tilde{\tau}_{i j}=\sum_{k, l=1}^{3} C_{i j k l} \frac{1}{s_{x} s_{y} s_{z}} \frac{\partial \frac{1}{s_{l}} \tilde{u}_{k}}{\partial x_{l}} \end{aligned}，(15)

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

\begin{aligned} &\rho \omega^{2} \bar{\tilde{u}}_{i}=\sum_{j=1}^{3} \frac{1}{s_{j}} \frac{\partial \bar{\tilde{\tau}}_{i j}}{\partial x_{j}}\\ &\bar{\tilde{\tau}}_{i j}=\sum_{k, l=1}^{3} C_{i j k l} \frac{1}{s_{l}} \frac{\partial \bar{\tilde{u}}_{k}}{\partial x_{l}}\\ &\bar{\tilde{u}}_{i}=\frac{\tilde{u}_{i}}{s_{x} s_{y} s_{z}}\\ &\bar{\tilde{\tau}}_{i j}=\frac{\tilde{\tau}_{i j}}{s_{x} s_{y} s_{z}} \end{aligned}，(16)

将（16）与（12）比较，我们发现除了NPML区域的所有物理场都乘以三个方向的拉伸坐标函数外，它们的表达式是相似的。差别在于，标准PML是PML区中的 \tilde{u}_{i} 和 \tilde{\tau}_{i j} 的方程，而NPML是 \bar{\tilde{u}}_{i} 和 \bar{\tilde{\tau}}_{ij} 的方程。因为伸展坐标函数在物理区域与NPML之间的界面上，它确保了界面两边 \bar{\tilde{u}}_{i} 和 \bar{\tilde{\tau}}_{ij} 的连续性。因此，NPML和标准PML公式在这种变化下等价。将方程（14）变换到时间域，并重新整理，我们得到：

\begin{aligned} &\rho \frac{\partial^{2} u_{i}}{d t^{2}}=\sum_{j=1}^{3} \frac{\partial \bar{\tau}_{i j}}{\partial x_{j}}\\ &\tau_{i j}=\sum_{k, l=1}^{3} C_{i j k l} \frac{\partial \bar{u}_{k l}}{\partial x_{l}} \end{aligned}，(17)

其中

\begin{aligned} &\bar{\tau}_{i j}=\frac{1}{s_{j}} \tau_{i j}\\ &\bar{u}_{k l}=\frac{1}{s_{l}} u_{k} \end{aligned}，(18)

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

【未完待续】