声波导管理论及仿真

本文首先介绍了声波导管的基础理论，并在此基础上利用COMSOL对导管进行模拟，以此验证理论的正确性。

---

前言

一般声源在无界空间中辐射的常常是波阵面逐步发散的球面波，现在将声的辐射约束在管子中，以此探讨自然管子的形状、尺寸以及管壁材料还有声源的状态对管中声波传播的影响。

在此展现部分仿真结果，后文将对其进行理论分析。

图1 导管的三维声压图和截面高次波

一、矩形声波导管理论分析

设有如图2所示的一矩形管，其宽度为 $l_y$，高度为 $l_x$，管长用 $z$坐标表示。已知三维波动方程为：
$$\frac{{\partial }^{2}p}{\partial x^2}+\frac{{\partial }^{2}p}{\partial y^2}+\frac{{\partial }^{2}p}{\partial z^2}=\frac{1}{c_0^2}\frac{{\partial }^{2}p}{\partial t^2}\text{\hspace{1em}(1)}$$

图2 矩阵导管

现令解为：
$$p=p_{\mathrm{a}}(x,y,z){\mathrm{e}}^{\mathrm{j}\omega t}\text{\hspace{1em}(2)}$$

代入方程(1)可得：
$$\frac{{\partial }^2{p}_{\mathrm{a}}}{\partial x^2}+\frac{{\partial }^2{p}_{\mathrm{a}}}{\partial y^2}+\frac{{\partial }^2{p}_{\mathrm{a}}}{\partial z^2}+k^2{p}_{\mathrm{a}}=0,\text{\hspace{1em}(3)}$$

式中，$k=\frac{\omega }{c_0}$，对方程(3)作分离变量，设：
$$p_{\mathrm{a}}(x,y,z)=X(x)Y(y)Z(z)\text{\hspace{1em}(4)}$$

于是可以得到三个独立坐标的常微分方程：
$$\begin{array}{c}\frac{{\partial }^{2}X(x)}{\partial x^2}+k_x^{2}X(x)=0,\\ \frac{{\partial }^{2}Y(y)}{\partial y^2}+k_y^{2}Y(y)=0,\\ \frac{{\partial }^{2}Z(z)}{\partial z^2}+k_z^{2}Z(z)=0,\end{array}\text{\hspace{1em}(5)}$$

其中 $k_x$，$k_y$，$k_z$为三个待定常数，它们之间满足以下关系：
$$k^2=k_x^2+k_y^2+k_z^2\text{\hspace{1em}(6)}$$

考虑到管子的 $x$，$y$ 方向有界，将存在驻波，因为方程(5)中第一和第二方程的解为以下形式：
$$\begin{array}{rlr}X(x)=& A_x\cos k_{x}x+B_x\sin k_{x}x& (7a)\\ Y(y)=& A_y\cos k_{y}y+B_y\sin k_{y}y& (7b)\end{array}$$

第三方程考虑到 $z$方向无限，没有反射波，因而取行波解为：
$$\begin{array}{rl}Z(z)=A_z{\mathrm{e}}^{-\mathrm{j}{k}_{z}z}& (7c)\end{array}$$

从(7a)与(7b)式可求得 $x$，$y$方向上的质点速度为：
$$\begin{array}{rl}{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)\left[\frac{\partial X(x)}{\partial x}\right]{\mathrm{e}}^{\mathrm{j}\omega t}\\ \text{=}& \frac{\mathrm{j}{k}_x}{{\rho }_0\omega }Y(y)Z(z)(-A_x\sin k_{x}x+B_x\cos k_{x}x){\mathrm{e}}^{\mathrm{j}\omega t}\\ v_y=& \frac{\mathrm{j}{k}_y}{{\rho }_0\omega }X(x)Z(z)(-A_y\sin k_{y}y+B_y\cos k_{y}y){\mathrm{e}}^{\mathrm{j}\omega t}\end{array}\text{\hspace{1em}(8)}$$

根据刚性管壁边界条件：
$$\begin{gathered} (v_x{)}_{(x=0,l_x)}=0 \\ (v_y{)}_{(y=0,l_y)}=0 \end{gathered}$$

可得：
$$\begin{gathered} B_x=0,k_x{l}_x=n_x\pi (n_x=0,1,2,\cdots ) \\ {B}_y=0,k_y{l}_y=n_y\pi (n_y=0,1,2,\cdots )\text{\hspace{1em}(9)} \end{gathered}$$

于是方程(2)可以转变为：
$$p_{n_x{n}_y}=A_{n_x{n}_y}\cos k_{x}x\cos k_{y}y{e}^{j(\omega t-k_{z}z)}\text{\hspace{1em}(10)}$$

这里 $p_{n_x{n}_y}$ 为每一组（$n_x$，$n_y$）数值对应的一个特解，它表示了在声波导管中可能存在的沿 $z$ 方向传播的一种声波。这种声波的圆频率为 $\omega$ ，传播速度为 $c_z=\frac{\omega }{k_z}$ ，振幅由 $A_{n_x{n}_y}\cos k_{x}x\cos k_{y}y$ 决定。将方程(6)和方程(9)联立，可以得到：
$$k_z={\left[\begin{array}{c}{k}^2-(k_x^2+k_y^2)\end{array}\right]}^{\frac{1}{2}}={\left(\frac{{\omega }^2}{c_0^2}-{\beta }_{n_x{n}_y}^2\right)}^{\frac{1}{2}}\text{\hspace{1em}(11)}$$

而：
$${\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\text{\hspace{1em}(12)}$$

我们可知仅当 $k_z$ 为实数时，在 $z$ 方向才表现有波的传播。而从方程(11)可以看到，这一 $k_z$并不在任何条件下都为实数，因此欲在 $z$ 方向传播声波就必须满足如下条件：
$$\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\text{\hspace{1em}(13)}$$

综上，当 $\frac{{\omega }^2}{c_0^2}<{\beta }_{n_{x^{n}y}}^2$，那么方程(11)应化为 $k_z=-\mathrm{j}{\alpha }_{n_{x^{n}y}}$ ，其中 ${\alpha }_{n_x{n}_y}=\sqrt{{\beta }_{n_x{n}_y}^2-\frac{{\omega }^2}{c_0^2}}$ 为正的实数，于是方程(10)就变为：
$$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}\text{\hspace{1em}(14)}$$

此式表明，沿 $z$ 轴传播的声波在做衰减得到整体振动。

由此，我们把管中产生沿 $z$ 方向传播声波的条件归结为：
$$f>f_{n_x{n}_y}\text{\hspace{1em}(15)}$$

这里
$$f_{n_x{n}_y}=\frac{c_0}{2}\sqrt{{\left(\frac{n_x}{l_x}\right)}^2+{\left(\frac{n_y}{l_y}\right)}^2}\text{\hspace{1em}(16)}$$

称为声波导管的简正频率。

分析公式(10)可知，对于不同的一组（$n_x$，$n_y$）数值将得到不同的波。我们称对应于（$n_x$，$n_y$）的波为（$n_x$，$n_y$）次的简正波。并称（0，0）次波为主波，除（0，0）次以外的波称高次波。从上面分析可知，只有当声源的激发频率 $f$ 比管中某个简正频率 $f_{n_x{n}_y}$ 高时，才能在管中激发出对应的（$n_x$，$n_y$）次波。若声源的频率低于管中除零以外的最低一个简正频率，那么管中所有的高次波都不能出现。因为（0，0）次简正频率为0，所以只要有声波存在，任何频率都总是大于零，因此这是管中只可能传播唯一的（0，0）次波。为之我们称除零以外的一个最低简正频率为声波导管的截止频率。

二、COMSOL仿真分析

2.1、几何模型构建

建模模型如图3所示，结合前文的理论分析我们可以得到对应的边长值 $l_x=0.08m$，$l_y=0.14m$。左侧平面的九个点为点源，它们将发射方向为（1，1，1），声强为 $1W$ 的声波。最右侧为完美匹配层，用于模拟无限延长的 $z$ ，主要作用是吸收入射波，避免反射波。

图3 矩阵导管建模

2.2、压力声学频域求解

将声源的频率分别设为 $1kHz$，$2kHz$，$3kHz$，得到的声压结果如下。

图4 1kHz声压场

图5 2kHz声压场

图6 3kHz声压场

从图4-6中可以看出，不同频率的声源在导管中传播的声压场不同，这是由于在低频时无法通过高次波，只能通过主波。随着频率的升高，通过的高次波不断增加，导致了声压场在管中呈现出不同的结果。

2.3、高次波分析

通过式(16)
$$f_{n_x{n}_y}=\frac{c_0}{2}\sqrt{{\left(\frac{n_x}{l_x}\right)}^2+{\left(\frac{n_y}{l_y}\right)}^2}\text{\hspace{1em}(16)}$$

我们可以得到不同波次的（$n_x$，$n_y$）的数值如下：
$$\{\begin{array}{l}{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{array}$$

从中可以得到，频率为 $1kHz$ 时只能通过主波，频率为 $2kHz$ 时能通过主波和（0，1）的高次波，频率为 $3kHz$ 时能通过主波和4种高次波。

在COMSOL中，我们截取横截面，通过分析横截面的波形态，来确定是否有高次波通过。

当声源频率为 $1kHz$ 时，通过如图7所示的波形。

图7 1kHz横截面声压

图7的声压场在 $x$ 和 $y$ 轴上无任何方向波动，说明该通过的波是主波。

当声源频率为 $2kHz$ 时，通过如图8所示的波形。

图8 2kHz横截面声压

从图8中可以看出，横截面通过两种波形，分别为主波和（0，1）的高次波。

当声源频率为 $3kHz$ 时，除主波外，还能通过如图9-10所示的四种高次波。

图9 3kHz横截面声压

图10 3kHz横截面声压

从左上角到右下角依次为（1，1），（0，2），（1，0），（0，1）的高次波。同理论计算结果一致。

---

总结

当声源的频率大于特定波次的简正频率时，该简正频率产生的波会沿轴向不衰减传播，并呈现出周期性变化。相反，若声源的频率小于特定波次的简正频率时，该简正频率产生波的声压将会呈现指数衰减的特点。当声源的频率大于多个简正频率时，多种波形将在传播的过程中叠加，呈现出不同的声压结果。