二维波动方程数值模拟(非均匀介质)

本文主要是介绍二维波动方程数值模拟(非均匀介质)，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

二维波动方程数值模拟(非均匀介质)

0. 前言

最近学习了波动方程在弹性介质中的传播，笔者对其中的面波、反射波、折射波的产生原理产生了兴趣，本文想要就波运动的本质方程 – 波动方程, 直接通过数值模拟的方法, 看看是否在介质边界是否能产生反射波、折射波，并观察其传播特性。

两层波速不同的介质实验结果如下:

在这里插入图片描述

1. 二维波动方程

$\frac{\partial^2 u}{\partial t^2} = v^2(\frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2}) \tag{1}$

2. 二维波动方程差分形式建立

2.1 一阶偏导差分形式

对于任意的函数 $f$ , 若其在 $x$ 邻域内可导，则
$\dfrac{\partial f}{\partial x} = \lim_{\Delta x \to 0} \frac{f(x+ \Delta x) - f(x)}{\Delta x}$
令 $x_k, \ \ x + \Delta x = x_{k + 1}$ , 那么有:
$\frac{\partial f}{\partial x} = \lim_{\Delta x \to0}\frac{f(x_{k+1}) - f(x_k)}{\Delta x}$
当 $\Delta x$ 充分小时, 可以近似的认为
$\frac{\partial f}{\partial x} = \frac{f(x_{k+1}) - f(x_k)}{\Delta x}$

2.2 二阶偏导的差分形式

同理, 对于任意函数 $f$ , 若其在 $x$ 邻域内二阶可导, 则可近似认为:
$\frac{\partial^2 f}{\partial x^2} = \frac{f'(x_{k+1}) - f'(x_k)}{\Delta x}= \frac{f(x_{k+1}) - 2f(x_{k}) + f(x_{k-1})}{\Delta x^2} \tag{2}$
式中 $f^{'}$ 为 $f$ 对 $x$ 的偏导。

2.3 二维均匀介质波动方程差分形式

对于均匀介质：

我们将二维空间 $x, y$ 分别分割为 $n, m, l$ 个等分段，且满足这些等分段张成的小正方体区域足够小, 每一个区域所在的坐标为 $x_i, y_j, t_k$ , 每小段长度分别为 $\Delta x, \Delta y, \Delta t$ , 那么可以将波场 $u$ 按上式作差分形式展开:

记号 $u_{i, j}^k$ 表达经数值差分得到的波场 $u$ 在网格点 $x_i, y_j, t_k)$ 处的值

$\frac{u_{i, j}^{k - 1} - 2{u_{i, j}^k} + u_{i, j}^{k + 1}}{\Delta t^2} = v^2\left[\frac{u_{i - 1, j}^{k } - 2{u_{i, j}^k} + u_{i + 1, j}^{k}}{\Delta x^2} + \frac{u_{i, j - 1}^{k} - 2{u_{i, j}^k} + u_{i, j + 1}^{k}}{\Delta y^2}\right]$

经化简得:
$u_{i, j}^{k + 1} = r_x{(u_{i -1, j}^k + u_{i+1, j}^k)} + r_y(u_{i, j-1}^k, u_{i, j +1}^k) + 2(1 - r_x - r_y) u_{i, j}^{k} - u_{i, j}^{k - 1} \tag{3}$
式中
$r_x = \frac{\Delta t^2}{\Delta x^2}v^2, \ r_y = \frac{\Delta t^2}{\Delta y^2}v^2$
根据 $(3)$ 式, 如果我们有第 $k$ 步和第 $k - 1$ 步的波, 就可以推算出 $k + 1$ 步的波场。

2.4 二维非均匀介质波动方程差分形式

对于非均匀介质, 其本质在于各点的波速不同, 观察 $(3)$ 式, 实际上对于不同的点只有 $r_x, r_y$ 不同。

我们令 $\rho(x, y)$

记记号 $v_{i, j}$ 表示 $\rho(x_i, y_j)$ 。

记号 $r_{x}^{i, j}$ 表示 $\dfrac{\Delta t^2}{\Delta x^2}v_{i, j}^2$ ,

记号 $r_{y}^{i, j}$ 表示 $\dfrac{\Delta t^2}{\Delta y^2}v_{i, j}^2$ .

对于 $(3)$ 式中波长的每一个点 $u_{i, j}$ , 将 $r_x, r_y$ 改写为 $r_x^{i, j}, r_y^{i, j}$ 。

得到:
$u_{i, j}^{k + 1} = r_x^{i-1, j}\cdot u_{i -1, j}^k + r_x^{i+1, j}\cdot u_{i+1, j}^k \\ \ + r_y^{i, j-1}\cdot u_{i, j-1}^k + r_y^{i, j+1}\cdot u_{i, j +1}^k\\ \ + 2(1 - r_x^{i, j} - r_y^{i, j}) u_{i, j}^{k} - u_{i, j}^{k - 1} \tag{4}$
得到非均匀介质的波动方程差分形式。

3. 初值条件及边界条件

3.1 初值条件

选用震源:
$ae^{-b\left[(x-x_0)^2 + (y-y_0)^2\right]}$
$x_0, y_0$ 为震源位置, $a, b$ 为系数。

特别的, 从 $u_{i, j}^0$ 递推到 $u_{i, j}^1$ 时, 需要用到 $u_{i, j}^{-1}$ 时的值, 但我们并没有当前时步的波场情况。这里我们使用外插法来获得 $- 1$ 时步的波长, 即
$u_{i, j}^{-1} = u_{i,j}^0 - 2\Delta t v_0(x_i, y_j)$
$v_0$ 为初始状态下每一点的波动速度。

另外, 三点差分格式必须满足收敛条件:
$\frac{4 v^2 \Delta t^2}{\Delta x^2 + \Delta y^2} \le 1$

3.2 边界条件

波在传播到实验区域边界后, 如果将边界的波场值设定恒定 $0$ , 那么会产生反射波，影响波长结果, 为此我们需要一种吸收边界来抵消边界产生的反射波。

其方程满足如下形式:
$\frac{\partial u}{\partial n} - \frac{1}{v}\frac{\partial u}{\partial t} = 0$
$n$ 为边界的法向量。

以顶部的边界条件为例, 其差分形式满足:
$\frac{u_{1, j}^k - u_{0, j}^k}{\Delta y} = \frac{1}{v}\frac{u^{k +1}_{0, j} - u^k_{0, j}}{\Delta t}$
即
$u_{0, j}^{k+1} = (1 - A_y) u_{0, j}^k + A_y \cdot u_{1, j}^k \tag{5}$

$A_y = \frac{\Delta t}{\Delta y}v$

我们用与 $(4)$ 同样的方法可以将 $A_y$ 代换为 $A_y^{i, j}$ 得到非均匀介质下的吸收条件。

其余三个边界的吸收条件与 $(5)$ 式类似，不再赘述。

4. 代码实现

4.1 约定

mask蒙层

在传递参数的过程中, 我们传递一个数组 $v$ 表示一组波速情况, 程序会构建一个蒙层 $ma s k$ 来记录每一个点的波速索引。其中 $ma s k$ 与波场具有同样的大小, $ma s k [i, j]$ 表示 $x_i, y_j)$ 处的波的波速所在 $v$ 数组的下标。

比如：

我们假设
$300]\\$
如果 $ma s k [1, 2] = 1$ , 表示 $x_1, y_2)$ 处的波速为 $v [1]$ , 即 $200$ ,

如果 $ma s k [5, 3] = 0$ , 表示 $x_5, y_3)$ 处的波速为 $v [0]$ , 即 $100$

sep函数

在传递参数的过程中需要传递一个 $se p$ 参数来指定生成蒙层的方法。 $se p$ 是一个函数, 该函数传入两个参数 $i, j$ 返回 $ma s k [i, j]$ 的值。

比如:

sep = lambda x, y: 0 if x > 2 else 1

表示当 $x_i > 2 时$ , $ma s k [i, j]$ 的值为 $0$ , 否则为 $1$ 。

4.2 代码

import osimport numpy as np
import matplotlib.pyplot as plt
import seaborn as snsclass WaveEquation2D:# initify horizion wave equation 2d simulator.def __init__(self, mt, mx, my, nt, nx, ny, px, py, v, sep):self.mt = mtself.mx = mxself.my = myself.nt = ntself.nx = nxself.ny = nyself.v = np.array(v)self.sep = sepself.dx = self.mx / self.nxself.dy = self.my / self.nyself.dt = self.mt / self.ntself.x = np.arange(0, self.mx + self.dx, self.dx)self.y = np.arange(0, self.my + self.dy, self.dy)self.t = np.arange(0, self.mt + self.dt, self.dt)self.u0 = lambda r, c: 0.2 * \np.exp(-((r - px)**2 + (c - py)**2) / 0.01)self.v0 = lambda r, c: 0r = 4 * max(self.v) ** 2 * self.dt ** 2 / (self.dx ** 2 + self.dy ** 2)assert r < 1, f"r should be less then 1, r is {r}."self.Ax = self.v * self.dt / self.dxself.Ay = self.v * self.dt / self.dyself.rx = self.Ax ** 2self.ry = self.Ay ** 2self.rxy = 1 - self.rx - self.ryself.u = np.zeros((self.nt + 1, self.nx + 1, self.ny + 1))for i in range(1, self.nx):for j in range(1, self.ny):self.u[0, i, j] = self.u0(self.x[i], self.y[j])self.mask = np.array([[sep(self.x[i], self.y[j]) for j in range(0, ny + 1)]for i in range(0, nx + 1)])def info(self):print("wave arguments:")print(f"dx: {self.dx}, dy: {self.dy}, dt:{self.dt}")print(f"ax: {self.Ax}, ay: {self.Ay}")print(f"rx: {self.rx}, ry: {self.ry}, rxy: {self.rxy}")def apply_boundary(self, t):# absort 1 dimfor i in range(self.nx + 1):self.u[t, i, 0] =\self.Ay[self.mask[i, 1]] * self.u[t-1, i, 1] + \(1 - self.Ay[self.mask[i, 0]]) * self.u[t-1, i, 0]self.u[t, i, self.ny] =\self.Ay[self.mask[i, self.ny - 1]] * \self.u[t-1, i, self.ny-1] + \(1-self.Ay[self.mask[i, self.ny]]) * self.u[t-1, i, self.ny]for j in range(self.ny + 1):self.u[t, 0, j] =\self.Ax[self.mask[1, j]] * self.u[t-1, 1, j] + \(1 - self.Ax[self.mask[0, j]]) * self.u[t-1, 0, j]self.u[t, self.nx, j] =\self.Ax[self.mask[self.nx-1, j]] * \self.u[t-1, self.nx-1, j] + \(1-self.Ax[self.mask[self.nx, j]]) * self.u[t-1, self.nx, j]def solve(self):# solve t == 1for i in range(1, self.nx):for j in range(1, self.ny):self.u[1, i, j] =\0.5 * self.rx[self.mask[i-1, j]] * self.u[0, i-1, j] + \0.5 * self.rx[self.mask[i+1, j]] * self.u[0, i+1, j] + \0.5 * self.ry[self.mask[i, j-1]] * self.u[0, i, j-1] + \0.5 * self.ry[self.mask[i, j+1]] * self.u[0, i, j+1] + \self.rxy[self.mask[i, j]] * self.u[0, i, j] + \self.dt * self.v0(i, j)for t in range(2, self.nt):print('\rsolving step {} / {}'.format(t + 1, self.nt), end="")# fills boundary conditiosself.apply_boundary(t)for i in range(1, self.nx):for j in range(1, self.ny):self.u[t, i, j] = \self.rx[self.mask[i-1, j]] * self.u[t-1, i-1, j] + \self.rx[self.mask[i+1, j]] * self.u[t-1, i+1, j] + \self.ry[self.mask[i, j-1]] * self.u[t-1, i, j-1] + \self.ry[self.mask[i, j+1]] * self.u[t-1, i, j+1] + \2 * self.rxy[self.mask[i, j]] * self.u[t-1, i, j] -\self.u[t-2, i, j]def plot(self, img_path="./image", show=True, save=True):if not os.path.exists(img_path):os.mkdir(img_path)for t in range(0, self.nt):sns.heatmap(self.u[t, :, :],vmin=-np.max(self.u)/10,vmax=np.max(self.u)/10,cmap='seismic',xticklabels=False,yticklabels=False)if save:plt.savefig(os.path.join(img_path, str(t)))if show:plt.pause(0.01)plt.cla()plt.clf()if __name__ == '__main__':we = WaveEquation2D(mt=2,mx=4,my=4,nt=300,nx=80,ny=80,px=0,py=2,v=[3.0, 5.0],sep=lambda r, c: 1 if r > 2 else 0)we.solve()we.plot(show=True, save=True)

附录

A. 收敛性条件

为了方便后面描述，定义波场精确值
$u_{i, j, k} = u(x_i, y_j, t_k)$
与波场的差分值 $u_{i, j}^{k}$ 区分。

和二阶差分算子

$\delta_{t}^{2} U = U_{i, j, k + 1} - 2U_{i, j, k} + U_{i, j, k-1}\\ \delta_{x}^{2} U = U_{i+1, j, k} - 2U_{i, j, k} + U_{i-1, j, k}\\ \delta_{y}^{2} U = U_{i, j+1, k} - 2U_{i, j, k} + U_{i, j-1, k}$

算子的右下部分表示所要二阶差分的坐标。

精确值的差分格式可以写为

$\dfrac{\delta^2_t u}{\Delta t^2}=v^2\left[ \dfrac{\delta^2_x u}{\Delta x^2} + \dfrac{\delta^2_y u}{\Delta y^2} \right]$

在差分格式下，我们取截断误差
$T_{i, j, k} = \dfrac{\delta^2_t u}{\Delta t^2}-v^2\left[ \dfrac{\delta^2_x u}{\Delta x^2} + \dfrac{\delta^2_y u}{\Delta y^2} \right]$

此时截断误差是当 $u$ 为精确值时，由于差分格式所引入的误差。

上式经过化简可得
$\begin{aligned} u_{i, j, k+1} &= r_x{(u_{i -1, j, k}+ u_{i+1, j, k})} \\ & + r_y(u_{i, j-1, k}, u_{i, j +1, k}) \\ &+ 2(1 - r_x - r_y) u_{i, j, k}\\ &- u_{i, j, k-1}\\ & - (\Delta t)^2 T_{i, j, k} \end{aligned} \tag{6}$

定义真实误差
$e_{i, j, k} = u_{i, j}^k - u_{i, j, k}$
也就是在网格点上函数真实值与差分估计值的差。

用 $(6)$ 式减去 $(3)$ 式，可得

$\begin{aligned} e_{i, j, k+1} &= r_x{(e_{i -1, j, k}+ e_{i+1, j, k})} \\ & + r_y(e_{i, j-1, k}, e_{i, j +1, k}) \\ &+ 2(1 - r_x - r_y) e_{i, j, k}\\ &- e_{i, j, k-1}\\ & - (\Delta t)^2 T_{i, j, k} \end{aligned} \tag{7}$

如果满足 $1-r_x-r_y \ge 0$ , 则上式所有与 $e$ 有关的系数均为整数，打绝对值以后可以提到绝对值号的外面。
由三角不等式有
$\begin{aligned} |e_{i, j, k+1}| & \le r_x{(|e_{i -1, j, k}|+ |e_{i+1, j, k}|)} \\ & + r_y(|e_{i, j-1, k}|+ |e_{i, j +1, k}|) \\ &+ 2(1 - r_x - r_y) |e_{i, j, k}|\\ & + |e_{i, j, k-1}|\\ & + |(\Delta t)^2 T_{i, j, k}| \end{aligned} \tag{8}$

将截断误差 $T$ 的定义式进行泰勒展开，可知 $O(\Delta t^2 + \Delta x^2 + \Delta y^2)$ , 那么必存在 $M > 0$ , 使得任意的 $T$ 都有
$\le M(\Delta t^2 + \Delta x^2 + \Delta y^2)$ .
再令
$E_k = \sup_{i, j}|e_{i, j, k}|$
代入 (8) 得
$\begin{aligned} |e_{i, j, k+1}| & \le r_x{(E_k+ E_k)} \\ & + r_y(E_k+ E_k) \\ &+ 2(1 - r_x - r_y) E_m\\ & + E_{k-1}\\ & + (\Delta t)^2 M(\Delta t^2 + \Delta x^2 + \Delta y^2) \end{aligned}$
于是有
$E_{k+1} \le E_k + E_{k-1} + (\Delta t)^2 M(\Delta t^2 + \Delta x^2 + \Delta y^2)$
递推可得
$E_{k} \le 2 E_0 + k (\Delta t)^2 M(\Delta t^2 + \Delta x^2 + \Delta y^2)$
由于在初值给定，必然有 $E_0 = 0$ ,
于是此时有
$E_{k} \le k (\Delta t)^2 M(\Delta t^2 + \Delta x^2 + \Delta y^2)$

当 $\Delta t, \Delta x, \Delta y \to 0$ 时， $E_k \to 0$ , 可以说明收敛性。

综上，在条件 $r_x - r_y \ge 0$ 时, 即
$r_x + r_y \le 1 \\ \dfrac{\Delta t^2 v^2}{\Delta x^2} + \dfrac{\Delta t^2 v^2}{\Delta y^2} \le 1\\ \dfrac{\Delta t^2 v^2 (\Delta x^2 + \Delta y^2)}{\Delta x^2 \Delta y^2} \le 1$