本文首先介绍了声波导管的基础理论,并在此基础上利用COMSOL对导管进行模拟,以此验证理论的正确性。
前言
一般声源在无界空间中辐射的常常是波阵面逐步发散的球面波,现在将声的辐射约束在管子中,以此探讨自然管子的形状、尺寸以及管壁材料还有声源的状态对管中声波传播的影响。
在此展现部分仿真结果,后文将对其进行理论分析。
![](https://i-blog.csdnimg.cn/direct/1cbab0737c664fa98db2378feeffdaea.png)
![](https://i-blog.csdnimg.cn/direct/d27e29441d894f7c9eeb52f43ea23c55.png)
一、矩形声波导管理论分析
设有如图2所示的一矩形管,其宽度为
l
y
l_{y}
ly,高度为
l
x
l_{x}
lx,管长用
z
z
z坐标表示。已知三维波动方程为:
∂
2
p
∂
x
2
+
∂
2
p
∂
y
2
+
∂
2
p
∂
z
2
=
1
c
0
2
∂
2
p
∂
t
2
(1)
\frac{\partial^2p}{\partial x^2}+\frac{\partial^2p}{\partial y^2}+\frac{\partial^2p}{\partial z^2}=\frac1{c_0^2}\frac{\partial^2p}{\partial t^2} \tag{1}
∂x2∂2p+∂y2∂2p+∂z2∂2p=c021∂t2∂2p(1)
现令解为:
p
=
p
a
(
x
,
y
,
z
)
e
j
ω
t
(2)
p=p_\mathrm{a}(x,y,z)\mathrm{e}^{\mathrm{j}\omega t} \tag{2}
p=pa(x,y,z)ejωt(2)
代入方程(1)可得:
∂
2
p
a
∂
x
2
+
∂
2
p
a
∂
y
2
+
∂
2
p
a
∂
z
2
+
k
2
p
a
=
0
,
(3)
\frac{\partial^2p_\mathrm{a}}{\partial x^2}+\frac{\partial^2p_\mathrm{a}}{\partial y^2}+\frac{\partial^2p_\mathrm{a}}{\partial z^2}+k^2p_\mathrm{a}=0, \tag{3}
∂x2∂2pa+∂y2∂2pa+∂z2∂2pa+k2pa=0,(3)
式中,
k
=
ω
c
0
k=\frac\omega{c_0}
k=c0ω,对方程(3)作分离变量,设:
p
a
(
x
,
y
,
z
)
=
X
(
x
)
Y
(
y
)
Z
(
z
)
(4)
p_\mathrm{a}(x,y,z)=X(x)Y(y)Z(z) \tag{4}
pa(x,y,z)=X(x)Y(y)Z(z)(4)
于是可以得到三个独立坐标的常微分方程:
∂
2
X
(
x
)
∂
x
2
+
k
x
2
X
(
x
)
=
0
,
∂
2
Y
(
y
)
∂
y
2
+
k
y
2
Y
(
y
)
=
0
,
∂
2
Z
(
z
)
∂
z
2
+
k
z
2
Z
(
z
)
=
0
,
(5)
\begin{gathered} \frac{\partial^2X(x)}{\partial x^2}+k_x^2X(x)=0, \\ \frac{\partial^2Y(y)}{\partial y^2}+k_y^2Y(y)=0, \\ \frac{\partial^2Z(z)}{\partial z^2}+k_z^2Z(z)=0, \end{gathered} \tag{5}
∂x2∂2X(x)+kx2X(x)=0,∂y2∂2Y(y)+ky2Y(y)=0,∂z2∂2Z(z)+kz2Z(z)=0,(5)
其中
k
x
k_{x}
kx,
k
y
k_{y}
ky,
k
z
k_{z}
kz为三个待定常数,它们之间满足以下关系:
k
2
=
k
x
2
+
k
y
2
+
k
z
2
(6)
k^2=k_x^2+k_y^2+k_z^2 \tag{6}
k2=kx2+ky2+kz2(6)
考虑到管子的
x
x
x,
y
y
y 方向有界,将存在驻波,因为方程(5)中第一和第二方程的解为以下形式:
X
(
x
)
=
A
x
cos
k
x
x
+
B
x
sin
k
x
x
(
7
a
)
Y
(
y
)
=
A
y
cos
k
y
y
+
B
y
sin
k
y
y
(
7
b
)
\begin{aligned}X(x)=&A_x\cos k_xx+B_x\sin k_xx&(7a)\\Y(y)=&A_y\cos k_yy+B_y\sin k_yy&(7b)\end{aligned}
X(x)=Y(y)=Axcoskxx+BxsinkxxAycoskyy+Bysinkyy(7a)(7b)
第三方程考虑到
z
z
z方向无限,没有反射波,因而取行波解为:
Z
(
z
)
=
A
z
e
−
j
k
z
z
(
7
c
)
\begin{aligned}Z(z)=A_z\mathrm{e}^{-\mathrm{j}k_zz}&(7c)\end{aligned}
Z(z)=Aze−jkzz(7c)
从(7a)与(7b)式可求得
x
x
x,
y
y
y方向上的质点速度为:
v
x
=
−
1
j
ρ
0
ω
∂
p
∂
x
=
−
1
j
ρ
0
ω
Y
(
y
)
Z
(
z
)
[
∂
X
(
x
)
∂
x
]
e
j
ω
t
=
j
k
x
ρ
0
ω
Y
(
y
)
Z
(
z
)
(
−
A
x
sin
k
x
x
+
B
x
cos
k
x
x
)
e
j
ω
t
v
y
=
j
k
y
ρ
0
ω
X
(
x
)
Z
(
z
)
(
−
A
y
sin
k
y
y
+
B
y
cos
k
y
y
)
e
j
ω
t
(8)
\begin{aligned} v_{x}=& \frac{-1}{\mathrm{j}\rho_0\omega}\frac{\partial p}{\partial x}=\frac{-1}{\mathrm{j}\rho_0\omega}Y(y)Z(z)\biggl[\frac{\partial X(x)}{\partial x}\biggr]\mathrm{e}^{\mathrm{j}\omega t} \\ \text{=}& \frac{\mathrm{j}k_x}{\rho_0\omega}Y(y)Z(z)(-A_x\sin k_xx+B_x\cos k_xx)\mathrm{e}^{\mathrm{j}\omega t} \\ v_{y}=& \frac{\mathrm{j}k_y}{\rho_0\omega}X(x)Z(z)(-A_y\sin k_yy+B_y\cos k_yy)\mathrm{e}^{\mathrm{j}\omega t} \\ \end{aligned} \tag{8}
vx==vy=jρ0ω−1∂x∂p=jρ0ω−1Y(y)Z(z)[∂x∂X(x)]ejωtρ0ωjkxY(y)Z(z)(−Axsinkxx+Bxcoskxx)ejωtρ0ωjkyX(x)Z(z)(−Aysinkyy+Bycoskyy)ejωt(8)
根据刚性管壁边界条件:
(
v
x
)
(
x
=
0
,
l
x
)
=
0
(
v
y
)
(
y
=
0
,
l
y
)
=
0
(v_{x})_{(x=0,l_{x})}=0\\(v_{y})_{(y=0,l_{y})}=0
(vx)(x=0,lx)=0(vy)(y=0,ly)=0
可得:
B
x
=
0
,
k
x
l
x
=
n
x
π
(
n
x
=
0
,
1
,
2
,
⋯
)
B
y
=
0
,
k
y
l
y
=
n
y
π
(
n
y
=
0
,
1
,
2
,
⋯
)
(9)
B_{x}=0,k_{x}l_{x}=n_{x}\pi\quad(n_{x}=0,1,2,\cdots)\\B_{y}=0,k_{y}l_{y}=n_{y}\pi\quad(n_{y}=0,1,2,\cdots) \tag{9}
Bx=0,kxlx=nxπ(nx=0,1,2,⋯)By=0,kyly=nyπ(ny=0,1,2,⋯)(9)
于是方程(2)可以转变为:
p
n
x
n
y
=
A
n
x
n
y
cos
k
x
x
cos
k
y
y
e
j
(
ω
t
−
k
z
z
)
(10)
p_{n_xn_y}=A_{n_xn_y}\cos k_xx\cos k_yye^{j(\omega t-k_zz)} \tag{10}
pnxny=Anxnycoskxxcoskyyej(ωt−kzz)(10)
这里
p
n
x
n
y
p_{n_xn_y}
pnxny 为每一组(
n
x
n_{x}
nx,
n
y
n_{y}
ny)数值对应的一个特解,它表示了在声波导管中可能存在的沿
z
z
z 方向传播的一种声波。这种声波的圆频率为
ω
\omega
ω ,传播速度为
c
z
=
ω
k
z
c_z=\frac\omega{k_z}
cz=kzω ,振幅由
A
n
x
n
y
cos
k
x
x
cos
k
y
y
A_{n_xn_y}\cos k_xx\cos k_yy
Anxnycoskxxcoskyy 决定。将方程(6)和方程(9)联立,可以得到:
k
z
=
[
k
2
−
(
k
x
2
+
k
y
2
)
]
1
2
=
(
ω
2
c
0
2
−
β
n
x
n
y
2
)
1
2
(11)
k_z=\begin{bmatrix}k^2-(k_x^2+k_y^2)\end{bmatrix}^{\frac12}=\left(\frac{\omega^2}{c_0^2}-\beta_{n_xn_y}^2\right)^{\frac12} \tag{11}
kz=[k2−(kx2+ky2)]21=(c02ω2−βnxny2)21(11)
而:
β
n
x
n
y
2
=
[
(
n
x
l
x
)
2
+
(
n
y
l
y
)
2
]
π
2
(12)
\beta_{n_xn_y}^2=\left[\left(\frac{n_x}{l_x}\right)^2+\left(\frac{n_y}{l_y}\right)^2\right]\pi^2 \tag{12}
βnxny2=[(lxnx)2+(lyny)2]π2(12)
我们可知仅当
k
z
k_{z}
kz 为实数时,在
z
z
z 方向才表现有波的传播。而从方程(11)可以看到,这一
k
z
k_{z}
kz并不在任何条件下都为实数,因此欲在
z
z
z 方向传播声波就必须满足如下条件:
ω
2
c
0
2
>
β
n
x
n
y
2
=
[
(
n
x
l
x
)
2
+
(
n
y
l
y
)
2
]
π
2
(13)
\frac{\omega^{2}}{c_{0}^{2}}>\beta_{{n_{x}n_{y}}}^{2}=\left[\left(\frac{n_{x}}{l_{x}}\right)^{2}+\left(\frac{n_{y}}{l_{y}}\right)^{2}\right]\pi^{2} \tag{13}
c02ω2>βnxny2=[(lxnx)2+(lyny)2]π2(13)
综上,当
ω
2
c
0
2
<
β
n
x
n
y
2
\frac{\omega^2}{c_0^2}<\beta_{n_{x^ny}}^2
c02ω2<βnxny2,那么方程(11)应化为
k
z
=
−
j
α
n
x
n
y
k_{z}=-\mathrm{j}\alpha_{n_{x^{n}y}}
kz=−jαnxny ,其中
α
n
x
n
y
=
β
n
x
n
y
2
−
ω
2
c
0
2
\alpha_{n_xn_y}=\sqrt{\beta_{n_xn_y}^2-\frac{\omega^2}{c_0^2}}
αnxny=βnxny2−c02ω2 为正的实数,于是方程(10)就变为:
p
n
x
n
y
=
A
n
x
n
y
cos
n
x
π
l
x
x
cos
n
y
π
l
y
y
e
−
a
n
x
n
y
z
e
j
ω
t
(14)
p_{n_{x}n_{y}}=A_{n_{x}n_{y}}\cos \frac{n_{x}\pi}{l_{x}}x\cos \frac{n_{y}\pi}{l_{y}}y \mathrm{e}^{-a_{n_{x}n_{y}}z} \mathrm{e}^{j\omega t} \tag{14}
pnxny=Anxnycoslxnxπxcoslynyπye−anxnyzejωt(14)
此式表明,沿 z z z 轴传播的声波在做衰减得到整体振动。
由此,我们把管中产生沿
z
z
z 方向传播声波的条件归结为:
f
>
f
n
x
n
y
(15)
f>f_{n_xn_y} \tag{15}
f>fnxny(15)
这里
f
n
x
n
y
=
c
0
2
(
n
x
l
x
)
2
+
(
n
y
l
y
)
2
(16)
f_{n_xn_y}=\frac{c_0}2\sqrt{\left(\frac{n_x}{l_x}\right)^2+\left(\frac{n_y}{l_y}\right)^2} \tag{16}
fnxny=2c0(lxnx)2+(lyny)2(16)
称为声波导管的简正频率。
分析公式(10)可知,对于不同的一组( n x n_{x} nx, n y n_{y} ny)数值将得到不同的波。我们称对应于( n x n_{x} nx, n y n_{y} ny)的波为( n x n_{x} nx, n y n_{y} ny)次的简正波。并称(0,0)次波为主波,除(0,0)次以外的波称高次波。从上面分析可知,只有当声源的激发频率 f f f 比管中某个简正频率 f n x n y f_{n_xn_y} fnxny 高时,才能在管中激发出对应的( n x n_{x} nx, n y n_{y} ny)次波。若声源的频率低于管中除零以外的最低一个简正频率,那么管中所有的高次波都不能出现。因为(0,0)次简正频率为0,所以只要有声波存在,任何频率都总是大于零,因此这是管中只可能传播唯一的(0,0)次波。为之我们称除零以外的一个最低简正频率为声波导管的截止频率。
二、COMSOL仿真分析
2.1、几何模型构建
建模模型如图3所示,结合前文的理论分析我们可以得到对应的边长值
l
x
=
0.08
m
l_{x}=0.08m
lx=0.08m,
l
y
=
0.14
m
l_{y}=0.14m
ly=0.14m。左侧平面的九个点为点源,它们将发射方向为(1,1,1),声强为
1
W
1W
1W 的声波。最右侧为完美匹配层,用于模拟无限延长的
z
z
z ,主要作用是吸收入射波,避免反射波。
2.2、压力声学频域求解
将声源的频率分别设为 1 k H z 1kHz 1kHz, 2 k H z 2kHz 2kHz, 3 k H z 3kHz 3kHz,得到的声压结果如下。
从图4-6中可以看出,不同频率的声源在导管中传播的声压场不同,这是由于在低频时无法通过高次波,只能通过主波。随着频率的升高,通过的高次波不断增加,导致了声压场在管中呈现出不同的结果。
2.3、高次波分析
通过式(16)
f
n
x
n
y
=
c
0
2
(
n
x
l
x
)
2
+
(
n
y
l
y
)
2
(16)
f_{n_xn_y}=\frac{c_0}2\sqrt{\left(\frac{n_x}{l_x}\right)^2+\left(\frac{n_y}{l_y}\right)^2} \tag{16}
fnxny=2c0(lxnx)2+(lyny)2(16)
我们可以得到不同波次的(
n
x
n_{x}
nx,
n
y
n_{y}
ny)的数值如下:
{
n
x
=
0
,
n
y
=
0
;
f
=
0
n
x
=
1
,
n
y
=
0
;
f
=
2143.75
n
x
=
0
,
n
y
=
1
;
f
=
1225
n
x
=
1
,
n
y
=
1
;
f
=
2469.06
n
x
=
2
,
n
y
=
0
;
f
=
4287.5
n
x
=
0
,
n
y
=
2
;
f
=
2450
\begin{cases} n_x=0, n_y=0; f=0\\ n_x=1,n_y=0; f=2143.75\\ n_x=0,n_y=1; f=1225\\ n_x=1,n_y=1; f=2469.06\\ n_x=2,n_y=0; f=4287.5\\ n_x=0,n_y=2; f=2450\\ \end{cases}
⎩
⎨
⎧nx=0,ny=0;f=0nx=1,ny=0;f=2143.75nx=0,ny=1;f=1225nx=1,ny=1;f=2469.06nx=2,ny=0;f=4287.5nx=0,ny=2;f=2450
从中可以得到,频率为 1 k H z 1kHz 1kHz 时只能通过主波,频率为 2 k H z 2kHz 2kHz 时能通过主波和(0,1)的高次波,频率为 3 k H z 3kHz 3kHz 时能通过主波和4种高次波。
在COMSOL中,我们截取横截面,通过分析横截面的波形态,来确定是否有高次波通过。
当声源频率为
1
k
H
z
1kHz
1kHz 时,通过如图7所示的波形。
图7的声压场在 x x x 和 y y y 轴上无任何方向波动,说明该通过的波是主波。
当声源频率为 2 k H z 2kHz 2kHz 时,通过如图8所示的波形。
![](https://i-blog.csdnimg.cn/direct/45dae686d0a648e3a154893f5c8aa7cb.png)
![](https://i-blog.csdnimg.cn/direct/9eaabb5274884662afb6001f234d9fab.png)
从图8中可以看出,横截面通过两种波形,分别为主波和(0,1)的高次波。
当声源频率为 3 k H z 3kHz 3kHz 时,除主波外,还能通过如图9-10所示的四种高次波。
![](https://i-blog.csdnimg.cn/direct/764c4ebf1f374943a131576d1c82eb70.png)
![](https://i-blog.csdnimg.cn/direct/4ae67231c25e4c9d908aa351edbf5d18.png)
![](https://i-blog.csdnimg.cn/direct/03b552161009457d8f1bcfaf36a5a924.png)
![](https://i-blog.csdnimg.cn/direct/6991c77e2dbe4cf99fca7d56f8d2906a.png)
从左上角到右下角依次为(1,1),(0,2),(1,0),(0,1)的高次波。同理论计算结果一致。
总结
当声源的频率大于特定波次的简正频率时,该简正频率产生的波会沿轴向不衰减传播,并呈现出周期性变化。相反,若声源的频率小于特定波次的简正频率时,该简正频率产生波的声压将会呈现指数衰减的特点。当声源的频率大于多个简正频率时,多种波形将在传播的过程中叠加,呈现出不同的声压结果。