声波传播方程-Helmholtz方程差分计算

1.声波传播方程程序

$${\nabla }^{2}p+k^{2}p=0\text{\hspace{1em}(1)}$$

其中 ， $k=\omega /c,\omega =2\pi f$，采用中心差分格式：${\nabla }^{2}p=\frac{{\partial }^{2}p}{\partial x^2}+\frac{{\partial }^{2}p}{\partial y^2}$ ，其中:
$$\frac{{\partial }^{2}p}{\partial x^2}=\frac{p_{i+1,j}+p_{i-1,j}-2p_{i,j}}{{(\Delta x)}^2}$$$$\frac{{\partial }^{2}p}{\partial y^2}=\frac{p_{i,j+1}+p_{i,j-1}-2p_{i,j}}{{(\Delta y)}^2}\text{\hspace{1em}(2)}$$
将(2)带入(1)式中得：
$$\frac{p_{i+1,j}+p_{i-1,j}-2p_{i,j}}{{(\Delta x)}^2}+\frac{p_{i,j+1}+p_{i,j-1}-2p_{i,j}}{{(\Delta y)}^2}+k^2{p}_{i,j}=0\text{\hspace{1em}(3)}$$

(3)式为Helmholtz方程差分计算的基本方程。因此将声压变量向量形式得到：
$$AP=0\text{\hspace{1em}(4)}$$
由于差分方程(3)中的变量数目为5，因此A是带宽为5的稀疏矩阵。$P=[p_{1,1},p_{1,2},p_{1,3}\cdots p_{m,n}{]}^T$为声压列向量。由于矩阵A是满秩的，因此方程(4)仅有0解。在实际频域计算时往往存在一些激励，比如平面波入射，单极子、偶极子点源等。
假设位置$r0$处含有单极子电源，此时方程(1)可表示如下形式：
$${\nabla }^{2}p+k^{2}p=4\pi S\delta (r-r0)\text{\hspace{1em}(5)}$$
其中$r0$处的激励为 ， 此时(4)式右边存在与$r0$点对应非零项b，矩阵方程(4)转化为：
$$AP=b\text{\hspace{1em}(6)}$$
其中$b$为已知列向量，然后可以利用迭代求解得到未知量 。

2.边界条件

根据离散方程(3)的差分格式，我们注意到，角边界点只有三个变量，而非角边界只有4个变量。这是由于在形成矩阵的时候默认了最外层边界为0。实际上在单个点源激励时，无穷远处才会存在声压为0的现象，因此在计算区域限制的情况下，必然会出现反射问题。为了减小反射的影响，除了扩大计算域外，还可以添加吸收层，其中最常用的吸收层为完美匹配层(PML),其基本原理是将复数空间引入到Helmholtz方程之中，如同在边界上添加了一层阻尼，声压传递到外层时会呈指数下降，因此声压传递到最外层时已经大小趋近于0，反射对计算域产生的影响较小。

当然边界条件设定除了满足计算需求外，还可以简化计算，可以对称界面设定，比如以$p_{1,j}$为对称点即$\frac{p_{2,j}+p_{0,j}-2p_{i,j}}{{(\Delta x)}^2}+\frac{p_{1,j+1}+p_{1,j-1}-2p_{i,j}}{{(\Delta y)}^2}+k^2{p}_{1,j}=0$可转化为：
$$\frac{2p_{2,j}-2p_{i,j}}{{(\Delta x)}^2}+\frac{p_{1,j+1}+p_{1,j-1}-2p_{i,j}}{{(\Delta y)}^2}+k^2{p}_{1,j}=0\text{\hspace{1em}(7)}$$
$p_{2,j}$对应系数增加一倍。
以上便是频域差分方程求解的整体思路，要注意的是与时域的求解方式link不同，椭圆形的频域方程无法采用递进计算的方式求解。交流学习请联系qq2862274738或邮箱dynyanhao@163.com