数值求解波动方程 [3]

2022-08-23 15:19 作者:nyasyamorina 0人读过 | 我要投稿

完美匹配层

之前提到的几种截断方式: Dirichlet 边界, Neumann 边界和循环边界, 这几种截断方式都有一个共同特征: 就是波会在边界处反射回来 (循环边界可以看作一边的波被"反射"到另一边). 但对于实际应用来说, 很多情况下要求边界不能产生反射, 以避免反射波对内部区域造成影响. 因为不存在反射波, 对于内部区域来说, 波好像消散在无限空间里.

有一种实现无反射"边界"的技术叫完美匹配层 (perfectly matched layer), 简称 PML. 准确来说这不是一个边界条件, 而是一个靠近边界的区域. 在对 PML 进行展开细说之前, 首先要对波动方程进行开刀.

引入一个新量 $%5Cvec%20v(x_1%2Cx_2%2C%5Ccdots%2Ct)$ 使得 $%5Cfrac%7B%5Cpartial%5Cvec%20v%7D%7B%5Cpartial%20t%7D%3A%3D%5Cvec%5Cnabla%20u$ , 代入波动方程方程 $%5Cfrac%7B%5Cpartial%5E2%20u%7D%7B%5Cpartial%20t%5E2%7D%3Dc%5E2%5Cvec%5Cnabla%5E2u$ 得 $%5Cfrac%7B%5Cpartial%5E2%20u%7D%7B%5Cpartial%20t%5E2%7D%3Dc%5E2%5Cvec%5Cnabla%5Ccdot%5Cfrac%7B%5Cpartial%5Cvec%20v%7D%7B%5Cpartial%20t%7D%3D%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20t%7D%5Cvec%5Cnabla%5Ccdot%5Cvec%20v$ , 即 $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D%3Dc%5E2%5Cvec%5Cnabla%5Ccdot%5Cvec%20v$ .

如此一元二阶微分方程 $%5Cfrac%7B%5Cpartial%5E2%20u%7D%7B%5Cpartial%20t%5E2%7D%3Dc%5E2%5CDelta%20u$ 分解为二元一阶微分方程组 $%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D%5Cfrac%7B%5Cpartial%5Cvec%20v%7D%7B%5Cpartial%20t%7D%3D%5Cvec%5Cnabla%20u%5C%5C%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D%3Dc%5E2%5Cvec%5Cnabla%5Ccdot%5Cvec%20v%5Cend%7Bmatrix%7D%5Cright.$ .

题外话: 可以记 $%5Cvec%20w%3D%5Cbegin%7Bbmatrix%7Du%5C%5C%5Cvec%20v%5Cend%7Bbmatrix%7D$ , 则有 $%5Cfrac%7B%5Cpartial%5Cvec%20w%7D%7B%5Cpartial%20t%7D%3D%5Chat%20D%5Cvec%20w$ , 其中 $%5Cvec%20D%3D%5Cbegin%7Bbmatrix%7D%26c%5E2%5Cvec%5Cnabla%5Ccdot%5C%5C%5Cvec%5Cnabla%26%5Cend%7Bbmatrix%7D$ 是反厄米算符. 一般地, 如果一个方程可以写为 $%5Cfrac%7B%5Cpartial%5Cvec%20w%7D%7B%5Cpartial%20t%7D%3D%5Chat%20D%5Cvec%20w%0A$ 形式, 并且 $%5Chat%20D$ 为反厄米算符, 那么 $%5Cvec%20w$ 就会产生 "类波" 解: 麦克斯韦方程组, 薛定谔方程, Lamé-Navier 方程都符合. 这也是为什么薛定谔方程又叫薛定谔波动方程, 有些人别整天揪着 "薛定谔波动方程" 不放了, 是薛定谔 "波动" 方程, 不是薛定谔 "波动方程".

PML 背后的思想就是变换坐标, 来看看变换坐标如何做到 "吸收" 波.

定义坐标变换 $x%5Cmapsto%5Ctilde%20x%3Dx%2Bif(x)$ , 其中 $%5Cforall%20x%3A%5C%3Af(x)%5Cgeq0$ . 当 $f(x)%3D0$ 时, 坐标变换为自身, 那么 $e%5E%7Bikx%7D%0A$ 的图像可以展示为

改变 $f(x)$ 后可以得到:

可以看到波在 $f(x)%3E0$ 的区域内呈指数下降, 如同波在这个区域内被 "吸收" 了一样. 实际上, $%20e%5E%7Bikx%7D$ 在施加坐标变换 $x%5Cmapsto%5Ctilde%20x%3Dx%2Bif(x)$ 得到 $e%5E%7Bikx%7D%5Cmapsto%20e%5E%7Bik(x%2Bif(x))%7D%3De%5E%7Bikx%7De%5E%7B-kf(x)%7D$ , 其中项 $e%5E%7B-kf(x)%7D$ 即表现为指数衰减.

当使用 PML 包围内部区域, 并且在 PML 外进行截断, 既然波在 PML 区域内呈指数衰减, 那么只要 PML 足够厚, 波在到达边界时可以变得足够小, 即使波在边界处被反射, 反射波也无法在内部区域产生足够大的影响.

PML 内使用的坐标变换称为 PML 变换. PML 变换与上面展示的坐标变换类似: 由 $%5Ctilde%20x%3Dx%2Bif(x)$ 得到 $d%5Ctilde%20x%3D%5Cleft(1%2Bi%5Cfrac%7Bdf%7D%7Bdx%7D%5Cright)dx$ . 在 PML 变换里记 $%5Cfrac%7Bdf%7D%7Bdx%7D%3D%5Csigma(x)$ , 但坐标变换为 $d%5Ctilde%20x%3D%5Cleft(1%2Bi%5Cfrac%7B%5Csigma(x)%7D%7B%5Comega%7D%5Cright)dx$ , 引入频率 $%5Comega%0A$ 的原因是上面求得的衰减项 $e%5E%7B-kf(x)%7D$ 与波数 $k$ 有关, 而频率和波数又有 $%5Comega%3Dc%7Ck%7C$ , 也就是说衰减速度于频率成正比, 这会造成低频分量无法有效地被吸收. 在引入频率后对应的衰减项为 $e%5E%7B-%5Cfrac%7B1%7D%7Bc%7D%5Cint%5Ex%5Csigma(x)dx%7D$ , 衰减速度与频率无关. 为了确保波在 PML 内总是衰减的, 应该有 $%5Cforall%20x%3A%5C%3A%5Csigma(x)%5Cgeq0$ .

下面以一维做例子, 二元一阶微分方程组版的一维波动方程为 $%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20t%7D%3D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D%5C%5C%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D%3Dc%5E2%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20x%7D%5Cend%7Bmatrix%7D%5Cright.$ . 初值条件 $%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Du(x%2C0)%5C%5C%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D(x%2C0)%5Cend%7Bmatrix%7D%5Cright.$ 可以根据关系变为 $%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Du(x%2C0)%5C%5Cv(x%2C0)%3D%5Cfrac%7B1%7D%7Bc%5E2%7D%5Cint%5Ex%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D(x'%2C0)dx'%5Cend%7Bmatrix%7D%5Cright.$ .

计算过程以 $%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20t%7D%3D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D$ 为例: 在 PML 内有 $v%3De%5E%7Bi(kx-%5Comega%20t)%7D$ 和 $%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20t%7D%3D-i%5Comega%20v$ , 即 $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D%3D-i%5Comega%20v$ . 施加 PML 变换 $%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D%5Cmapsto%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%5Ctilde%20x%7D%3D%5Cfrac%7B1%7D%7B1%2Bi%5Cfrac%7B%5Csigma%7D%7B%5Comega%7D%7D%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D$ , 得到 $%5Cfrac%7B1%7D%7B1%2Bi%5Cfrac%7B%5Csigma%7D%7B%5Comega%7D%7D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D%3D-i%5Comega%20v$ , 然后得到 $-i%5Comega%20v%3D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D-%5Csigma%20v$ , 即 $%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20t%7D%3D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D-%5Csigma%20v$ .

同理可得 $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D%3Dc%5E2%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20x%7D-%5Csigma%20u$ . 在内部区域里, 应该有 $%5Csigma%3D0$ , 这时 PML 版本波动方程退化为普通的波动方程, 而在 PML 区域内则有 $%5Csigma%3E0$ .

这里选取 $%5Cfrac%7Bu_%7Bx%2Ct%2B1%7D-u_%7Bx%2Ct%7D%7D%7B%5CDelta%20t%7D$ 作为 $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D$ 的离散化也不会造成太大问题 (节省内存), 那么 PML 版波动方程离散化为 $%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Dv_%7Bx%2Ct%2B1%7D%3D%5Cfrac%7B%5CDelta%20t%7D%7B2%5CDelta%20x%7D(u_%7Bx%2B1%2Ct%7D-u_%7Bx-1%2Ct%7D)%2B(1-%5Csigma_x%5CDelta%20t)v_%7Bx%2Ct%7D%5C%5Cu_%7Bx%2Ct%2B1%7D%3D%5Cfrac%7Bc_x%5E2%5CDelta%20t%7D%7B2%5CDelta%20x%7D(v_%7Bx%2B1%2Ct%7D-v_%7Bx-1%2Ct%7D)%2B(1-%5Csigma_x%5CDelta%20t)u_%7Bx%2Ct%7D%5Cend%7Bmatrix%7D%5Cright.$ . 并且这种实现不存在数值稳定的参数选择, 最多只能做到较小的时间数值不稳定.

需要注意的是, 在连续波动方程里, PML 是无论 $%5Csigma$ 如何选取都是没有反射波的, 但离散化后则不是这样, 为了确保离散化后的 PML 仍然没有反射, 需要满足 $%5Csigma$ 具有连续的二阶或三阶导数.

这里介绍一个常用具有连续二阶导数的 $%5Csigma$ : 以右边界为例, 假设 $x%5Cin(0%2C1)$ 为 PML, 那么可以选取 $%5Csigma(x)%3Ds%5Cleft(x-%5Cfrac%7B%5Csin(2%5Cpi%20x)%7D%7B2%5Cpi%7D%5Cright)$ , 其中 $s$ 为 $%5Csigma$ 的最大值, 对于不同位置不同厚度的 PML 可以使用坐标变换把这个 σ 映射过去. 刚刚说过, 尽管 PML 可以快速吸收波, 但截断的边界仍会产生微小的反射波返回内部区域, 当选取 Dirichlet 或 Neumann 边值条件时, 反射波的相对大小由 $R%3De%5E%7B-%5Cfrac%7BsL%7D%7Bc%7D%7D$ 给出, 其中 $L$ 为 PML 的厚度. 那么 PML 系统如下图所示: