压力声学·一

本节涉及部分压力声学的基本理论，comsol的简单介绍和压力升学接口的部分使用方式，以及使用comsol进行仿真的一个简单案例。本节主要内容来自于comsol视频课程中心的对应模块教程。

1. 压力声学的理论基础

1.1 压力声学的研究对象

在 COMSOL 的声学模块中，不同的物理假设对应不同声学接口。声-结构耦合适用于流体声场与固体结构振动相互作用的问题；气动声学适用于存在显著背景流动的声学问题；热粘性声学适用于微小尺度结构中粘性边界层与热边界层不可忽略的问题；几何声学则适用于几何尺度远大于波长、可以近似使用射线描述声传播的问题。压力声学则适用于流体中以压力波为主的传播问题。

具体来说，压力声学研究的是流体介质中的小幅声学扰动传播问题。流体包括空气、水、软组织、多孔材料中的等效流体区域，或者其他主要通过压缩与膨胀传递声能的连续介质。在这类问题中，声波被建模为纵波，即介质质点的振动方向与波传播方向一致。因此，压力声学的核心变量是声压扰动 $p$。

压力声学是根据其假设得到的简化模型。这种简化的优势在于计算效率高，主变量少，适合用于常规声场、共振腔、管道、扬声器辐射、声压分布和频域响应等问题的建模。然而，当声场附近存在显著边界层损耗、背景流动、柔性结构振动或固体弹性波传播时，单纯压力声学可能不再充分，则需借助上述其他接口。

1.2 连续流体的控制方程

声学方程来源于连续流体力学方程。对于一般流体，控制方程主要包括质量守恒方程、动量守恒方程、能量守恒方程以及状态方程：

{ \[\begin{aligned} & \frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf{u})=0,\\& \frac{\partial(\rho\mathbf{u})}{\partial t} +\nabla\cdot(\rho\mathbf{u}\mathbf{u}^{T}) = \nabla\cdot\boldsymbol{\sigma}+\mathbf{F},\\& \rho C_p\frac{\partial T}{\partial t} + \alpha_0T\frac{\partial p}{\partial t} = -\nabla\cdot(\mathbf{q}-\mathbf{u}\cdot\boldsymbol{\tau})+Q,\\& \rho=\rho(p,T). \end{aligned}\]

. $$

其中，第一式为质量守恒方程，描述流体密度随时间的变化与质量通量之间的关系。$\rho$ 表示流体密度，$\mathbf{u}$ 表示流体速度矢量。第二式为动量守恒方程，描述流体动量随时间的变化及其与应力和体力之间的关系。其中，$\boldsymbol{\sigma}$ 是应力张量，$\mathbf{F}$ 是体力项。第三式为能量守恒方程，用于描述温度、压力、热传导、粘性耗散和热源之间的关系。其中，$T$ 是温度，$C_p$ 是定压热容，$\alpha_0$ 是热膨胀系数，$\mathbf{q}$ 是热通量，$\boldsymbol{\tau}$ 是粘性应力张量，$Q$ 是热源项。第四式为状态方程，是密度与压力和温度的关系。

对于牛顿流体，应力张量可表示为

\[ \boldsymbol{\sigma} = -p\mathbf{I} + \boldsymbol{\tau}, \]

其中 $p$ 是压力，$\mathbf{I}$ 是单位张量，$\boldsymbol{\tau}$ 是粘性应力张量。粘性应力张量通常写为

\[ \boldsymbol{\tau} = \mu\left[\nabla\mathbf{u}+(\nabla\mathbf{u})^{T}\right] - \left(\frac{2}{3}\mu-\mu_{B}\right) (\nabla\cdot\mathbf{u})\mathbf{I}, \]

其中 $\mu$ 为动力粘度，$\mu_B$ 为体积粘度。热通量通常满足傅里叶热传导定律

\[ \mathbf{q}=-k\nabla T, \]

其中 $k$ 是热导率。

1.3 小扰动假设与线性声学

声波通常被看作叠加在背景状态上的小幅扰动。设背景状态为静态或稳态场，声学扰动远小于背景量，则总温度、总速度、总压力和总密度可分别表示为背景量与一阶扰动量之和：

$$ \[\begin{aligned} T_{\mathrm{tot}}(\mathbf{x},t) &= T_0(\mathbf{x})+T_1(\mathbf{x},t),\\ \mathbf{u}_{\mathrm{tot}}(\mathbf{x},t) &= \mathbf{u}_0(\mathbf{x})+\mathbf{u}_1(\mathbf{x},t),\\ p_{\mathrm{tot}}(\mathbf{x},t) &= p_0(\mathbf{x})+p_1(\mathbf{x},t),\\ \rho_{\mathrm{tot}}(\mathbf{x},t) &= \rho_0(\mathbf{x})+\rho_1(\mathbf{x},t). \end{aligned}\]

在频域分析中，扰动量进一步假设为简谐形式。例如，声压扰动可以写成

\[ p_1(\mathbf{x},t) = \mathrm{Re}\left\{p(\mathbf{x})e^{i\omega t}\right\}, \]

其中 $p(\mathbf{x})$ 是复声压幅值，$\omega$ 是角频率。类似地，温度、速度和密度扰动可写为

\[ \begin{aligned} T_1(\mathbf{x},t) &= \mathrm{Re}\left\{T(\mathbf{x})e^{i\omega t}\right\}, \\ \mathbf{u}_1(\mathbf{x},t) &= \mathrm{Re}\left\{\mathbf{u}(\mathbf{x})e^{i\omega t}\right\}, \\ \rho_1(\mathbf{x},t) &= \mathrm{Re}\left\{\rho(\mathbf{x})e^{i\omega t}\right\}. \end{aligned} \]

小扰动假设的核心条件是

\[ |p_1|\ll |p_0|,\qquad |\rho_1|\ll |\rho_0|,\qquad |T_1|\ll |T_0|. \]

1.4 线性 Navier-Stokes 声学

以上所述，当介质中存在背景流动时，声学扰动是在运动介质中传播的，此时总速度可写为

\[ \mathbf{u}_{\mathrm{tot}}(\mathbf{x},t) = \mathbf{u}_0(\mathbf{x}) + \mathbf{u}(\mathbf{x})e^{i\omega t}, \]

其中 $\mathbf{u}_0$ 是背景流速，$\mathbf{u}$ 是声学扰动速度。其他变量也可以写作类似像是。将这些形式代入 Navier-Stokes 方程并保留一阶扰动项，可得到线性 Navier-Stokes 声学模型。

$$ { \[\begin{aligned}& \frac{\partial \rho}{\partial t} +\nabla\cdot(\rho\mathbf{u}_0+\rho_0\mathbf{u}) =M, \\& \rho_0\left( \frac{\partial \mathbf{u}}{\partial t} +(\mathbf{u}\cdot\nabla)\mathbf{u}_0 +(\mathbf{u}_0\cdot\nabla)\mathbf{u} \right) +\rho(\mathbf{u}_0\cdot\nabla)\mathbf{u}_0 = \nabla\cdot\boldsymbol{\sigma} +\mathbf{F} -\mathbf{u}_0M,& \\& \rho_0C_p\left( \frac{\partial T}{\partial t} +\mathbf{u}\cdot\nabla T_0 +\mathbf{u}_0\cdot\nabla T \right) +\rho C_p(\mathbf{u}_0\cdot\nabla T_0) -\alpha_pT_0\left( \frac{\partial p}{\partial t} +\mathbf{u}\cdot\nabla p_0 +\mathbf{u}_0\cdot\nabla p \right) -\alpha_pT(\mathbf{u}_0\cdot\nabla p_0) = \nabla\cdot(k\nabla T)+\Phi+Q, \\& \rho = \rho_0(\beta_Tp-\alpha_pT). \end{aligned}\]

. $$

其中 $\beta_T$ 是等温压缩率，$\alpha_p$ 是热膨胀系数。线性 Navier-Stokes 声学比普通压力声学更一般，保留了背景流动、粘性和热效应对声波传播的影响。

在介质可以近似认为无平均流动时，即

\[ \mathbf{u}_0=0. \]

可进一步简化为无流动声学问题。

1.5 热粘性声学

如果背景流动可以忽略，但几何尺度很小，声波在壁面附近形成的粘性边界层和热边界层不可忽略，则需要使用热粘性声学。热粘性声学保留声压、速度和温度扰动，并考虑粘性耗散与热传导损耗。

在无背景流动条件下，热粘性声学的连续性方程可写为

\[ \frac{\partial \rho_t}{\partial t} + \nabla\cdot(\rho_0\mathbf{u}_t) = 0. \]

其动量方程可写为

\[ \rho_0\frac{\partial \mathbf{u}_t}{\partial t} = \nabla\cdot \left[ -p_t\mathbf{I} + \mu\left(\nabla\mathbf{u}_t+(\nabla\mathbf{u}_t)^T\right) - \left(\frac{2}{3}\mu-\mu_B\right) (\nabla\cdot\mathbf{u}_t)\mathbf{I} \right]. \]

能量方程可概括为

\[ \rho_0C_p \left( \frac{\partial T_t}{\partial t} + \mathbf{u}_t\cdot\nabla T_0 \right) - \alpha_0T_0 \left( \frac{\partial p_t}{\partial t} + \mathbf{u}_t\cdot\nabla p_0 \right) = \nabla\cdot(\kappa\nabla T_t)+Q. \]

状态关系可表示为

\[ \rho_t = \rho_0(\beta_Tp_t-\alpha_pT_t), \]

与压力声学相比，热粘性声学的变量更多，计算成本更高，但能够描述微型声学结构中的边界层耗散。例如微孔、MEMS 声学器件、窄缝、耳机泄漏缝隙、小型声腔等问题中，热粘性损耗可能决定系统响应。

1.6 从热粘性声学到压力声学

压力声学方程是由连续流体控制方程在特定近似条件下化简得到的。其基本假设包括：介质无平均流动、忽略粘性损耗、忽略热传导损耗，并认为声学扰动相对于背景量足够小。设背景介质均匀且静止，即

\[ \mathbf{u}_0=0,\qquad \rho_0=\text{const.},\qquad p_0=\text{const.} \]

在小扰动假设下，压力、密度和速度分别写为

\[ p_{\mathrm{tot}}=p_0+p, \]

\[ \rho_{\mathrm{tot}}=\rho_0+\rho, \]

\[ \mathbf{u}_{\mathrm{tot}}=\mathbf{u}. \]

其中 $p$、$\rho$ 和 $\mathbf{u}$ 分别表示声压扰动、密度扰动和声学质点速度。将这些变量代入质量守恒方程，并忽略二阶小量，可得线性化连续性方程：

\[ \frac{\partial \rho}{\partial t} + \rho_0\nabla\cdot\mathbf{u} = 0. \]

在无粘、无背景流动条件下，动量方程退化为线性欧拉方程：

\[ \rho_0 \frac{\partial \mathbf{u}}{\partial t} = -\nabla p. \]

同时，在绝热声学近似下，压力扰动与密度扰动满足

\[ p=c^2\rho, \]

即

\[ \rho=\frac{p}{c^2}, \]

其中 $c$ 是介质中的声速。

由此可以消去 $\rho$ 和 $\mathbf{u}$，得到只关于声压 $p$ 的方程。首先将

\[ \rho=\frac{p}{c^2} \]

代入线性化连续性方程，得到

\[ \frac{1}{c^2} \frac{\partial p}{\partial t} + \rho_0\nabla\cdot\mathbf{u} = 0. \]

对时间再求一次导数，有

\[ \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2} + \rho_0\nabla\cdot \frac{\partial \mathbf{u}}{\partial t} = 0. \]

由线性欧拉方程可得

\[ \frac{\partial \mathbf{u}}{\partial t} = -\frac{1}{\rho_0}\nabla p. \]

将其代入上式，得到

\[ \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2} - \nabla^2p = 0. \]

这就是无源、无损、无背景流动条件下的声压波动方程。若写成更适合有限元形式的散度形式，则为

\[ \nabla\cdot \left( -\frac{1}{\rho_0}\nabla p \right) + \frac{1}{\rho_0c^2} \frac{\partial^2p}{\partial t^2} = 0. \]

如果进一步考虑声源项，可以在方程右端加入单极源 $Q$，并在梯度项中加入偶极源 $\mathbf{q}_d$，于是得到带源项的时域压力声学方程：

\[ \nabla\cdot \left[ -\frac{1}{\rho_0} (\nabla p-\mathbf{q}_d) \right] + \frac{1}{\rho_0c^2} \frac{\partial^2p}{\partial t^2} = Q. \]

这里，$Q$ 表示单极声源，例如体积注入型声源；$\mathbf{q}_d$ 表示偶极声源，例如方向性力源或偶极型声源。若不考虑声源，则有

\[ Q=0,\qquad \mathbf{q}_d=\mathbf{0}. \]

所有变量都被消去，只剩下声压 $p$ 作为唯一主变量。

若考虑频域简谐声场，令

\[ p(\mathbf{x},t) = \mathrm{Re} \left\{ p(\mathbf{x})e^{i\omega t} \right\}, \]

则有

\[ \frac{\partial^2p}{\partial t^2} \rightarrow -\omega^2p. \]

将其代入无源时域波动方程，得到频域压力声学方程：

\[ \nabla\cdot \left( -\frac{1}{\rho_0}\nabla p \right) - \frac{\omega^2p}{\rho_0c^2} = 0. \]

这就是压力声学频域形式的 Helmholtz 方程。有源情况下，则可写为

\[ \nabla\cdot \left[ -\frac{1}{\rho_0} (\nabla p-\mathbf{q}) \right] - \frac{\omega^2p}{\rho_0c^2} = Q. \]

1.7 压力声学中的派生变量

虽然压力声学以声压 $p$ 为主变量，但其他声学量仍然可以通过声压推导得到。对于频域谐波问题，线性化动量方程给出声压与质点速度之间的关系。若采用 $e^{i\omega t}$ 的时间因子，则可写为

\[ i\omega\rho_0\mathbf{u} = \nabla p. \]

在绝热条件下，密度扰动与压力扰动之间满足

\[ \rho = \frac{p}{c^2}. \]

温度扰动也可通过热力学关系由声压推得。例如，在小扰动条件下可写为

\[ \rho_0C_pT = pT_0\alpha_0. \]

2. 特征频率与频域响应

2.1 无源分析与有源分析

压力声学分析通常可以分为无源分析和有源分析。无源分析关注系统自身的固有性质，例如特征频率和模态分布；有源分析关注系统在外部声源激励下的实际响应，例如某一频率下的声压分布、扫频响应、传递函数和声压级。

特征频率分析属于无源分析。它不需要定义外部声源，而是求解在给定几何、材料参数和边界条件下，声学系统自身允许存在的非零振动模式。频域响应分析则属于有源分析，需要给定声源、激励频率和边界条件，求解声场在外部强迫下的响应。

2.2 特征频率方程

在频域压力声学中，无源 Helmholtz 方程为

\[ \nabla\cdot \left( -\frac{1}{\rho_0}\nabla p \right) - \frac{\omega^2p}{\rho_0c^2} = 0. \]

在特征频率分析中，角频率 $\omega$ 不再是预先指定的参数，而是需要求解的特征值。COMSOL 中常将其写为与特征值 $\lambda$ 有关的形式：

\[ \nabla\cdot \left( -\frac{1}{\rho_0}\nabla p \right) + \frac{\lambda^2p}{\rho_0c^2} = 0. \]

若时间因子采用 $e^{i\omega t}$，则特征值与角频率之间存在相应关系。特征频率 $f$ 与角频率 $\omega$ 的关系为

\[ \omega=2\pi f. \]

特征频率方程中没有外部源项 $Q$，也没有偶极源项 $\mathbf{q}$。因此，特征频率不是由声源决定的，而是由几何形状、材料参数、边界条件和介质声速决定的。声源的作用是在频域响应分析中激发这些模态，而不是改变这些模态本身。

2.3 特征频率与声源

频域有源问题可以抽象写为

\[ \mathcal{L}(p,\omega)=Q, \]

其中，$\mathcal{L}$ 是由几何形状、材料参数和边界条件共同决定的线性算子，$p$ 是待求解的声压场，$\omega$ 是外部指定的激励角频率，$Q$ 是声源项。该方程表示：在给定声源和给定频率下，求解声学系统的受迫响应。

与之对应，特征频率问题是无源问题，可写为

\[ \mathcal{L}(p,\lambda)=0. \]

这里右端没有外部声源项，因此该方程是齐次方程。对于任意 $\lambda$，平凡解

\[ p=0 \]

总是存在，但该解没有物理意义。特征频率分析真正寻找的是使齐次方程存在非零解的特殊 $\lambda$，即

\[ p(\mathbf{x})\neq 0. \]

因此，特征频率问题的本质是寻找满足下列条件的特征值与特征模态：

\[ \mathcal{L}(p,\lambda)=0, \qquad p(\mathbf{x})\neq 0. \]

其中，$\lambda$ 是特征值，对应系统的固有频率；$p(\mathbf{x})$ 是对应的特征模态，表示该固有频率下声压的空间分布。换言之，特征频率并不是由外部声源激发出来的频率，而是在没有外部声源时，系统自身由几何、材料和边界条件所允许的自由振动频率。

有源频域响应与特征频率之间的关系可以从非齐次线性方程的解来理解。对于有源问题

\[ \mathcal{L}(p,\omega)=Q, \]

其解可以形式上表示为

\[ p=p_h+p_p, \]

其中 $p_h$ 是齐次解，满足

\[ \mathcal{L}(p_h,\omega)=0, \]

而 $p_p$ 是特解，满足

\[ \mathcal{L}(p_p,\omega)=Q. \]

在实际声学问题中，更常用的理解方式是将有源响应表示为各个特征模态的加权叠加：

\[ p(\mathbf{x},\omega) = \sum_n a_n(\omega)\phi_n(\mathbf{x}), \]

其中 $\phi_n(\mathbf{x})$ 是第 $n$ 阶特征模态，$a_n(\omega)$ 是该模态在当前声源和当前激励频率下的参与系数。该系数不仅与激励频率 $\omega$ 有关，也与声源的空间位置、相位和分布有关。一般来说，当激励频率接近某一阶特征频率时，该模态更容易被强烈激发；但如果声源位置与该模态的空间分布不匹配，例如声源位于该模态的节点附近，则该模态的响应也可能较弱。

因此，声源与特征频率之间的关系可以概括为：声源通常不决定系统的特征频率，但决定哪些特征模态会被激发，以及被激发的强弱。系统的特征频率由几何形状、材料参数和边界条件决定；有源频域响应则由系统特征、声源形式和激励频率共同决定。

例如，在房间声学模型中，房间的固有频率主要由房间尺寸、空气声速和墙面边界条件决定。改变扬声器的位置或声源幅值通常不会改变房间的固有频率，但会改变不同房间模态的激发程度，从而改变实际测得的声压响应曲线。

需要注意的是，上述结论成立的前提是声源被理想化为外部给定的压力源、速度源或体源项，而不改变系统本身的边界条件和动力学性质。如果声源本身具有明显机械阻抗，例如扬声器振膜、动铁单元或压电片，并且该结构与声腔发生强耦合，那么声源就不再只是外部激励，而会成为系统边界或结构的一部分。此时，系统的特征频率可能会发生改变，问题也应从纯压力声学问题扩展为声-结构耦合问题。

3. 典型理论问题

3.1 刚性圆柱声腔的特征频率

为了理解特征频率的物理含义，可以考虑一个刚性圆柱声腔。设圆柱半径为 $R$，长度为 $L$，腔内介质声速为 $c$，密度为 $\rho_0$。在无源、无损、无背景流动条件下，声压满足 Helmholtz 方程

\[ \nabla^2p+k^2p=0, \]

其中

\[ k=\frac{\omega}{c}. \]

如果圆柱壁面为声硬边界，则法向质点速度为零。对于压力声学，在均匀密度和无偶极源条件下，这等价于声压法向导数为零：

\[ \frac{\partial p}{\partial n}=0. \]

对于二维圆形截面，声压模态可在极坐标下写为

\[ p_{mn}(r,\theta) = J_m(k_{mn}r) \left[ A\cos(m\theta)+B\sin(m\theta) \right], \]

其中 $J_m$ 是 $m$ 阶第一类 Bessel 函数，$m$ 为周向模态阶数，$n$ 为径向模态阶数。刚性壁边界要求

\[ \left. \frac{\partial p}{\partial r} \right|_{r=R} = 0. \]

因此有

\[ J_m'(k_{mn}R)=0. \]

令 $\alpha'_{mn}$ 表示 $J_m'(x)$ 的第 $n$ 个零点，则

\[ k_{mn} = \frac{\alpha'_{mn}}{R}. \]

对应的特征频率为

\[ f_{mn} = \frac{c}{2\pi} \frac{\alpha'_{mn}}{R}. \]

对于三维有限长圆柱，还需要考虑轴向模态。若两端也为声硬边界，则轴向波数可写为

\[ k_z = \frac{l\pi}{L}, \]

其中 $l=0,1,2,\cdots$。此时总波数满足

\[ k_{mnl}^2 = \left(\frac{\alpha'_{mn}}{R}\right)^2 + \left(\frac{l\pi}{L}\right)^2. \]

因此三维圆柱刚性声腔的特征频率为

\[ f_{mnl} = \frac{c}{2\pi} \sqrt{ \left(\frac{\alpha'_{mn}}{R}\right)^2 + \left(\frac{l\pi}{L}\right)^2 }. \]

该解析结果可用于验证 COMSOL 的特征频率计算。若在 COMSOL 中建立相同半径、长度、材料参数和声硬边界条件的圆柱模型，则其前几个特征频率应与上述理论结果接近。偏差主要来自几何离散误差、网格误差、边界条件设置差异以及数值求解精度。

3.2 刚性圆柱声腔特征频率的理论计算

设三维刚性圆柱声腔的半径和长度分别为

\[ R=1\ \mathrm{m}, \qquad L=10\ \mathrm{m}. \]

若腔内介质为空气，并取声速为

\[ c=343\ \mathrm{m/s}, \]

刚性圆柱声腔的特征频率公式为

\[ f_{mnl} = \frac{c}{2\pi} \sqrt{ \left(\frac{\alpha'_{mn}}{R}\right)^2 + \left(\frac{l\pi}{L}\right)^2 }. \]

代入 $R=1\ \mathrm{m}$、$L=10\ \mathrm{m}$ 和 $c=343\ \mathrm{m/s}$，可得

\[ f_{mnl} = \frac{343}{2\pi} \sqrt{ (\alpha'_{mn})^2 + \left(\frac{l\pi}{10}\right)^2 }. \]

其中，

\[ \frac{343}{2\pi} \approx 54.59. \]

因此，

\[ f_{mnl} \approx 54.59 \sqrt{ (\alpha'_{mn})^2 + \left(\frac{l\pi}{10}\right)^2 }. \]

其中，$m$ 表示声压沿圆柱周向 $\theta$ 方向的变化阶数，$n$ 表示声压沿径向 $r$ 方向的变化阶数，$l$ 表示声压沿轴向 $z$ 方向的变化阶数。

对于刚性圆柱声腔，最低阶径向均匀模态满足

\[ \alpha'_{00}=0. \]

因此，当只考虑轴向模态，即 $m=0$、$n=0$ 时，特征频率公式化简为

\[ f_{00l} = \frac{343}{2\pi} \left( \frac{l\pi}{10} \right). \]

整理得

\[ f_{00l} = \frac{343l}{20}. \]

因此，前几个纯轴向模态为

\[ f_{001} = \frac{343}{20} = 17.15\ \mathrm{Hz}, \]

\[ f_{002} = \frac{343\times 2}{20} = 34.30\ \mathrm{Hz}, \]

\[ f_{003} = \frac{343\times 3}{20} = 51.45\ \mathrm{Hz}, \]

\[ f_{004} = \frac{343\times 4}{20} = 68.60\ \mathrm{Hz}. \]

需要注意的是，$(m,n,l)=(0,0,0)$ 对应

\[ f_{000}=0, \]

它表示整个腔体内压力均匀变化的零频模式，通常不作为实际声学共振频率讨论。

接下来考虑具有横向声压分布的模态。常用的 Bessel 函数导数零点包括

\[ \alpha'_{11}=1.8412, \]

\[ \alpha'_{21}=3.0542, \]

\[ \alpha'_{01}=3.8317. \]

对于模态 $(m,n,l)=(1,1,0)$，有

\[ f_{110} = \frac{343}{2\pi} \sqrt{ (1.8412)^2 + 0 }. \]

因此

\[ f_{110} = 54.59\times 1.8412 \approx 100.51\ \mathrm{Hz}. \]

对于模态 $(m,n,l)=(2,1,0)$，有

\[ f_{210} = \frac{343}{2\pi} \sqrt{ (3.0542)^2 + 0 }. \]

因此

\[ f_{210} = 54.59\times 3.0542 \approx 166.73\ \mathrm{Hz}. \]

对于模态 $(m,n,l)=(0,1,0)$，有

\[ f_{010} = \frac{343}{2\pi} \sqrt{ (3.8317)^2 + 0 }. \]

因此

\[ f_{010} = 54.59\times 3.8317 \approx 209.17\ \mathrm{Hz}. \]

由此可得，若 $R=1\ \mathrm{m}$、$L=10\ \mathrm{m}$，并取空气声速 $c=343\ \mathrm{m/s}$，则刚性圆柱声腔的部分特征频率为

模态 $(m,n,l)$	$\alpha'_{mn}$	特征频率 $f_{mnl}$
$(0,0,1)$	$0$	$17.15\ \mathrm{Hz}$
$(0,0,2)$	$0$	$34.30\ \mathrm{Hz}$
$(0,0,3)$	$0$	$51.45\ \mathrm{Hz}$
$(0,0,4)$	$0$	$68.60\ \mathrm{Hz}$
$(1,1,0)$	$1.8412$	$100.51\ \mathrm{Hz}$
$(2,1,0)$	$3.0542$	$166.73\ \mathrm{Hz}$
$(0,1,0)$	$3.8317$	$209.17\ \mathrm{Hz}$

从计算结果可以看出，由于圆柱长度 $L=10\ \mathrm{m}$ 明显大于半径 $R=1\ \mathrm{m}$，最低的非零特征频率首先表现为轴向模态。其最低非零特征频率为

\[ f_{001}=17.15\ \mathrm{Hz}. \]

而第一个具有明显横向声压变化的模态为

\[ f_{110}\approx 100.51\ \mathrm{Hz}. \]

这说明，对于细长圆柱声腔，低频模态主要由轴向长度决定；随着频率升高，径向和周向模态才逐渐出现。

3.2 左侧压力源与右侧开放边界的一维声传播

另一个基础问题是长度为 $L$ 的一维管道声传播。设左侧边界为压力源，右侧边界为开放或无反射边界，介质声速为 $c$，密度为 $\rho_0$。在频域中，波数为

\[ k=\frac{\omega}{c}. \]

如果左端给定压力源

\[ p(0)=p_s, \]

且右端为理想无反射边界，则管道内只存在向右传播的行波。此时声压可写为

\[ p(x)=p_s e^{-ikx}. \]

对应的质点速度为

\[ u(x)=\frac{p(x)}{\rho_0c}. \]

该情况下声波从左向右传播并离开计算区域，不形成明显驻波，因此不会出现由右端反射导致的强烈共振峰。理想无反射边界在数值仿真中通常通过辐射边界、端口边界或完美匹配层近似实现。

如果右端不是无反射边界，而是声软边界，即

\[ p(L)=0, \]

则解可写为

\[ p(x) = p_s \frac{\sin[k(L-x)]}{\sin(kL)}. \]

当分母趋近于零时，系统出现强烈响应，即

\[ \sin(kL)=0. \]

因此共振条件为

\[ kL=n\pi, \]

对应频率为

\[ f_n=\frac{nc}{2L}. \]

如果右端为声硬边界，即

\[ \left. \frac{\partial p}{\partial x} \right|_{x=L} = 0, \]

则声压可写为

\[ p(x) = p_s \frac{\cos[k(L-x)]}{\cos(kL)}. \]

共振条件变为

\[ \cos(kL)=0, \]

即

\[ kL=\frac{(2n+1)\pi}{2}. \]

对应频率为

\[ f_n = \frac{(2n+1)c}{4L}. \]

4. 使用COMSOL进行仿真

4.1 特征频率分析

首先使用 COMSOL 对刚性圆柱声腔进行仿真。

点击新建->模型向导->三维->声学->压力声学->压力声学，频域

点击选择物理场下侧添加，将该接口添加进 “添加到物理场接口”，点击研究。

由于本次研究特征频率，因此选择一般研究中的特征频率。点击完成，进入以下界面。

首先添加材料，在左侧组件1中右键材料，选择第一项，从库中添加材料：

选择用于仿真的材料，这里我们搜索并选择air。双击目标材料，材料将出现在左侧组件1的材料栏中。其中，材料明细中会显示此材料的各项属性值。例如，空气的平均摩尔质量被设置为0.2897Kg/mol等等。

然后增加几何。我们计划对刚性圆柱声腔进行仿真，因此不必使用其他建模软件建模。右键组件1中的几何1，选择圆柱体，指定半径和高度（这里我们选择R = 1，L = 10，与我们之前的假设对应），点击构建选定对象，圆柱体将显示在右侧。此时圆柱体的主轴为Z轴，如果希望完全仿照之前的理论计算，可以指定轴类型。由于不影响计算，我们保持。

计算前的最后一步是构建网格。构建网格有多种方式。这里简单起见，我们双击网格，选择最大网格单元大小控制参数，改为频率，并指定最大频率：

一般来说，最大频率和最大网格单元大小的关系是：

\[ 最大网格单元大小 = 声速/最大频率/5 \]

注意越大的网格会带来越大的计算量。根据之前的理论计算，主要关注250Hz之内的特征频率。

点击全部构建，完成网格构建。

点击研究1的步骤1：特征频率。可选栏中包括围绕偏移量搜索特征频率，这一项与下一选单共同决定哪些特征频率会被计算。所需特征频率数顾名思义。我们将后者改为10，点击计算。

在右侧结果中选择声压或声压级（声压的分贝表示），即可查看仿真结果。

可见部分结果与我们的理论计算相同，例如17.15Hz处的 (0,0,1) 模态。

4.2 频率分析

接下来我们尝试一些添加声源的频域分析效果。

首先右键模型开发器最上部的root栏，增加一个新研究，选择频域。或者也可以直接重新建一个文件并重复上述材料、几何、网络步骤。

在频域栏，指定频率：

注意不要超过之前网络中设置的最大频率，或重新设置网络。

我们首先尝试平面辐射波，其本质是一个开边界的声场输入。

右键压力声学，频域，

在辐射条件中选择平面辐射波，然后在边界中选择顶部和底部两个边界。

此时，这两个边界将会变成开边界，声学信号遇到此处边界将完全不会反射。

接下来，右键平面波辐射，点击入射压力场。设置入射压力长的压力幅值、声速和材料。

选择顶部边界或底部边界作为入射边界。

如此声源部分构建完成。在研究2中点击计算，即可得到结果。

最终的结果显而易见的，是一个声学纵波。

我们尝试其他边界，比如将原本出射的边界改为软声场边界。

或不使用入射压力场，改为使用一个球形声源：

模态 \((m,n,l)\)	\(\alpha'_{mn}\)	特征频率 \(f_{mnl}\)
\((0,0,1)\)	\(0\)	\(17.15\ \mathrm{Hz}\)
\((0,0,2)\)	\(0\)	\(34.30\ \mathrm{Hz}\)
\((0,0,3)\)	\(0\)	\(51.45\ \mathrm{Hz}\)
\((0,0,4)\)	\(0\)	\(68.60\ \mathrm{Hz}\)
\((1,1,0)\)	\(1.8412\)	\(100.51\ \mathrm{Hz}\)
\((2,1,0)\)	\(3.0542\)	\(166.73\ \mathrm{Hz}\)
\((0,1,0)\)	\(3.8317\)	\(209.17\ \mathrm{Hz}\)