本文主要是介绍声波传播方程-Helmholtz方程差分计算，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
1.声波传播方程程序
∇ 2 p + k 2 p = 0 (1) {\nabla ^2}p + {k^2}p = 0 \tag{1} ∇2p+k2p=0(1)
其中 ， k = ω / c , ω = 2 π f k = \omega /c,\omega = 2\pi f k=ω/c,ω=2πf，采用中心差分格式： ∇ 2 p = ∂ 2 p ∂ x 2 + ∂ 2 p ∂ y 2 {\nabla ^2}p = {{{\partial ^2}p} \over {\partial {x^2}}} + {{{\partial ^2}p} \over {\partial {y^2}}} ∇2p=∂x2∂2p+∂y2∂2p ，其中:
∂ 2 p ∂ x 2 = p i + 1 , j + p i − 1 , j − 2 p i , j ( Δ x ) 2 {{{\partial ^2}p} \over {\partial {x^2}}} = {{{p_{i + 1,j}} + {p_{i - 1,j}} - 2{p_{i,j}}} \over {{{(\Delta x)}^2}}} ∂x2∂2p=(Δx)2pi+1,j+pi−1,j−2pi,j ∂ 2 p ∂ y 2 = p i , j + 1 + p i , j − 1 − 2 p i , j ( Δ y ) 2 (2) {{{\partial ^2}p} \over {\partial {y^2}}} = {{{p_{i,j + 1}} + {p_{i,j - 1}} - 2{p_{i,j}}} \over {{{(\Delta y)}^2}}}\tag{2} ∂y2∂2p=(Δy)2pi,j+1+pi,j−1−2pi,j(2)
将(2)带入(1)式中得：
p i + 1 , j + p i − 1 , j − 2 p i , j ( Δ x ) 2 + p i , j + 1 + p i , j − 1 − 2 p i , j ( Δ y ) 2 + k 2 p i , j = 0 (3) {{{p_{i + 1,j}} + {p_{i - 1,j}} - 2{p_{i,j}}} \over {{{(\Delta x)}^2}}} + {{{p_{i,j + 1}} + {p_{i,j - 1}} - 2{p_{i,j}}} \over {{{(\Delta y)}^2}}} + {k^2}{p_{i,j}} = 0\tag{3} (Δx)2pi+1,j+pi−1,j−2pi,j+(Δy)2pi,j+1+pi,j−1−2pi,j+k2pi,j=0(3)
(3)式为Helmholtz方程差分计算的基本方程。因此将声压变量向量形式得到：
A P = 0 (4) AP = 0\tag{4} AP=0(4)
由于差分方程(3)中的变量数目为5，因此A是带宽为5的稀疏矩阵。 P = [ p 1 , 1 , p 1 , 2 , p 1 , 3 ⋯ p m , n ] T P = {[{p_{1,1}},{p_{1,2}},{p_{1,3}} \cdots {p_{m,n}}]^T} P=[p1,1,p1,2,p1,3⋯pm,n]T为声压列向量。由于矩阵A是满秩的，因此方程(4)仅有0解。在实际频域计算时往往存在一些激励，比如平面波入射，单极子、偶极子点源等。
假设位置 r 0 r0 r0处含有单极子电源，此时方程(1)可表示如下形式：
∇ 2 p + k 2 p = 4 π S δ ( r − r 0 ) (5) {\nabla ^2}p + {k^2}p = 4\pi S\delta (r - r0)\tag{5} ∇2p+k2p=4πSδ(r−r0)(5)
其中 r 0 r0 r0处的激励为 ， 此时(4)式右边存在与 r 0 r0 r0点对应非零项b，矩阵方程(4)转化为：
A P = b (6) AP = b\tag{6} AP=b(6)
其中 b b b为已知列向量，然后可以利用迭代求解得到未知量 。
2.边界条件
根据离散方程(3)的差分格式，我们注意到，角边界点只有三个变量，而非角边界只有4个变量。这是由于在形成矩阵的时候默认了最外层边界为0。实际上在单个点源激励时，无穷远处才会存在声压为0的现象，因此在计算区域限制的情况下，必然会出现反射问题。为了减小反射的影响，除了扩大计算域外，还可以添加吸收层，其中最常用的吸收层为完美匹配层(PML),其基本原理是将复数空间引入到Helmholtz方程之中，如同在边界上添加了一层阻尼，声压传递到外层时会呈指数下降，因此声压传递到最外层时已经大小趋近于0，反射对计算域产生的影响较小。
当然边界条件设定除了满足计算需求外，还可以简化计算，可以对称界面设定，比如以 p 1 , j {p_{1,j}} p1,j为对称点即 p 2 , j + p 0 , j − 2 p i , j ( Δ x ) 2 + p 1 , j + 1 + p 1 , j − 1 − 2 p i , j ( Δ y ) 2 + k 2 p 1 , j = 0 {{{p_{2,j}} + {p_{0,j}} - 2{p_{i,j}}} \over {{{(\Delta x)}^2}}} + {{{p_{1,j + 1}} + {p_{1,j - 1}} - 2{p_{i,j}}} \over {{{(\Delta y)}^2}}} + {k^2}{p_{1,j}} = 0 (Δx)2p2,j+p0,j−2pi,j+(Δy)2p1,j+1+p1,j−1−2pi,j+k2p1,j=0可转化为：
2 p 2 , j − 2 p i , j ( Δ x ) 2 + p 1 , j + 1 + p 1 , j − 1 − 2 p i , j ( Δ y ) 2 + k 2 p 1 , j = 0 (7) {{2{p_{2,j}} - 2{p_{i,j}}} \over {{{(\Delta x)}^2}}} + {{{p_{1,j + 1}} + {p_{1,j - 1}} - 2{p_{i,j}}} \over {{{(\Delta y)}^2}}} + {k^2}{p_{1,j}} = 0\tag{7} (Δx)22p2,j−2pi,j+(Δy)2p1,j+1+p1,j−1−2pi,j+k2p1,j=0(7)
p 2 , j {p_{2,j}} p2,j对应系数增加一倍。
以上便是频域差分方程求解的整体思路，要注意的是与时域的求解方式link不同，椭圆形的频域方程无法采用递进计算的方式求解。交流学习请联系qq2862274738或邮箱dynyanhao@163.com
这篇关于声波传播方程-Helmholtz方程差分计算的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
