李春熠，胡安琦，张格格，白在桥，俞骁翀

(北京师范大学 物理学系，北京 100875)

随着计算技术的进步，数值仿真已经成为与理论和实验并列的第三种科学研究方法.大学物理教学实验的选题一般比较简单，理论分析相对完善，但有时也会遇到涉及复杂几何形状或边界条件的波动或流体方面的问题.此时简化的理论往往不能给出精确的结果或精细的物理图像.如果针对这些实验把数值计算引入教学，不仅可以加深学生对所涉及物理问题的理解，还可以为他们以后的学习和研究打下一定基础.本文以声波谐振管为例，介绍数值仿真方法在基础物理实验教学中的应用.

声波在管中传播的衰减很小，在低频时可以用一维声波近似，因此实验教学中常利用声波谐振管来研究波的反射、叠加和谐振等现象[1-3].但实际空气中的声波是由三维波动方程:

(1)

描述的.该方程主要的约束是在硬边界(比如管壁)处声压p的法向导数为零.虽然一维波模型可以给出谐振频率的计算公式，但公式中的管口修正项并不能从一维波模型得到.此外，一些非理想情况也需要借助三维波动方程理解.

有限元法(Finite Element Method，FEM)是一种非常有效的计算偏微分方程近似解的数值方法.有限元法起源于20世纪中期对弹性力学问题的数值求解[4,5].后来扩展到计算流体力学、热学、电磁学、声学等领域，在科学研究和工程中有十分广泛的应用[6].随着计算机辅助设计(CAD)、复杂几何形状的自动网格划分、稳健的求解器等方面的进步，性能优秀、使用方便的有限元计算软件包应运而生，COMSOL Multiphysics就是其中的一个重要代表[7].

本文将介绍如何利用COMSOL分析声波谐振管中的物理问题.

1 COMSOL使用方法

有限元方法的原理是将待求问题的几何区域网格化，每个单元是简单的几何形状.然后通过变分法最小化能量泛函或加权残差法离散化，将偏微分方程转换为包含所有单元的代数方程组，最后用数值方法求解得到微分方程的近似解[7].

本文选择COMSOL的“压力声学/频域(acpr)”物理场，通过求解频域的三维亥姆霍兹方程对声波谐振管实验进行数值仿真.下面以测量部分封闭管的谐振频率为例说明仿真计算的步骤.在这个实验中，设想闭管的端面上有一个孔，通过改变孔的面积，可以实现从闭管到开管的连续过渡.设计这个实验的初衷是想探究闭管密封不严对谐振频率的影响.为了验证仿真结果的可靠性，本文还完成了真实实验，实验装置参照图1.

谐振管

1.1 构建三维几何模型

本实验用到的几何模型包括谐振管及两端的附件、扬声器、有不同孔的有机玻璃挡片.除此之外，为了保证声场方程的可解性，还需要定义“无穷远”处的边界条件.构建方法分述如下：

1) 谐振管管体

在三维直角坐标系中，不妨设圆柱形谐振管的中轴线为z轴.以坐标原点(0, 0, 0)为底面中心点，沿z轴正方向构建高度900 mm、半径20 mm的圆柱体(参见图2).对圆柱体进行“分层”，并设置“层的厚度”(即管壁厚度)为3 mm.

2) 扬声器

在谐振管左端管口附近[坐标为(0, 0, -90)]构建圆柱体，半径与高度均为30 mm.在真实实验中，为防止扬声器反射的声波进入谐振管(从而在扬声器与管口之间发生共振)，扬声器不能正对管口[2].相应地，将圆柱体绕x轴旋转45°角(参见图2).

图2 谐振管管体和扬声器模型的三维图像

3) 挡片

在右端管口和端盖之间构建带有不同尺寸和形状的开孔有机玻璃挡片(参见图3).构建方法为以(0, 0, 900)点为底面中心沿z轴正方向构建圆柱体.没有孔的挡片选择不分层.圆形孔的挡片采用构建圆柱体再“分层”的方法构建，“层的厚度”为圆柱体的半径减去圆形孔的半径；正方形孔、长方形孔的挡片则在不分层圆柱体的基础上构建正方体和长方体；半圆形孔的挡片在不分层圆柱体的基础上构建圆柱体和长方体取“交集”，如图4所示.模型的具体尺寸参数为：“有机玻璃片”圆柱体高度2 mm，半径23.5 mm；“圆形孔”圆柱体高度2 mm，半径23.5 mm，层的厚度为7.5～21.3 mm，间隔2 mm；“正方形孔”长方体高度2 mm，宽度、深度相等，为3～21 mm，间隔3 mm；“长方形孔”长方体高度2 mm，宽度24 mm，深度为2～24 mm，间隔2 mm；“半圆形孔”圆柱体高度2 mm，半径为2～16 mm，间隔2 mm .

图3 右端管口的yOz平面框架图

闭管

4) 完美匹配层

为了使问题有确定的解，COMSOL仿真用充分大的“完美匹配层”(Perfectly Matched Layers, PMLs)包围求解区域. PMLs完全吸收声波，不会产生反射，作用相当于定义了“无穷远”处的边界条件.构建高度2 000 mm、半径500 mm、层厚度300 mm的圆柱体将谐振管、扬声器等几何体包含在内，并将该圆柱体定义为“完美匹配层” (参见图5).

图5 完美匹配层模型的三维图像

目前所建的几何模型关于yOz平面对称，建立一个大的长方体与前面的几何模型取“交集”，将中心截面(yOz平面)暴露出来，方便后续选择区域和边界，以及呈现内部的声压分布等情况(参见图6).

图6 对半“切割”后的几何模型的三维图像

1.2 定义物理场和驱动

选择管内区域、挡片开孔区域以及完美匹配层区域为求解的空气介质区域(参见图7),并将以上区域做“对称”操作.空气介质区域的边界(包括管体、扬声器、挡片等固体)会自动被定义为“硬声场边界”.硬声场边界可视为理想的刚体.

左端管口局部图 右端管口局部图

将材料库中的空气“Air”添加到物理场中，采用“线弹性”流体模型，不考虑空气的热传导、黏性等耗散因素. 空气温度设置为真实实验温度值296.95 K(23.8 ℃)，绝对压强设为1 atm.

定义驱动的方法是：“添加声激励”，给扬声器圆柱体面朝向左端管口的底面(参见图8)施加“法向位移”1 mm，使之做简谐运动，模拟扬声器的膜片振动.

图8 扬声器的振动面

最后根据物理场划分网格.网格划分需要设置“最小单元大小”“最大单元增长率”等多个参数或选项.对于本文的问题，从物理上考虑，管内单元的尺寸应该小于波长的1/6(参见图9).

左端管口局部图

1.3 求解

对于声波谐振管，最基本的测量对象是谐振频率.添加“特征频率”的“研究”，设置特征频率搜索基准值和所需特征频率数，点击“计算”，让求解器在基准值附近按照“最近的绝对值”的方式搜索本征频率.以一端开口、一端封闭的情形(闭管)为例，设800 Hz为基准值，计算结果见表1.需要注意的是，虽然模型忽略了空气的耗散，但完美匹配层的吸声作用使得系统仍存在耗散，所以谐振频率一般为复数，其虚部可理解为共振峰的宽度或耗散的大小(快慢).仿真结果给出了很多谐振频率，其中只有虚部很小的才有物理意义(表1中方括号标记).图10显示了谐振频率662.95+0.6137 9i对应谐振管内的声压分布.不同灰度代表不同的声压值，白色为声压最大值，黑色为声压最小值，黑白交替处的声压为零(波节).图中可见明显的3个波节，因此模式数n=4.

图10 闭管发生谐振时的声压分布图

实验中经常会研究谐振频率对某个变量的依赖关系.这时可以将该变量设置为参数.添加“参数化扫描”，设置起始值、终止值和变化步长，使求解器自动在参数的定义域内扫描求解.例如，在改变长方形孔的深度时，可将其深度设置为参数d，为了使其关于yOz平面对称，y的位置设为-d/2，使d从2 mm到22 mm变化，步长为2 mm.这样就可以得到不同宽度矩形孔对应的谐振频率.

表1 闭管时的谐振频率仿真结果

续表

还可以添加“频域”的“研究”，设置起始频率、终止频率和变化步长，点击“计算”，即可自动绘制出各频率(包括谐振频率和非谐振频率)下的声压分布图和声压级分布图. 可将任何绘图的数据以文本格式或者CSV格式导出，便于进一步分析.

下一节将结合具体实验对仿真结果进行更仔细的分析.

2 仿真实验及结果分析

这一节介绍几个用COMSOL仿真的实验.

2.1 管口附近的声压分布

在一维波近似下，两端都为开口的谐振管(称为开管)的谐振频率经验公式为[2]

(2)

式中c为声速，n=1,2,3, …为模式数，而等效管长为

L=L0+0.8d

(3)

其中L0为实际管长，d为管的内直径.0.8d称为管口修正(end correction).从表达式上看，管口修正可理解为：声波可以从管口延伸出来，在距离管口0.4d的位置形成声压的波节.

为了探究管口修正的真实物理意义，本文计算了开管情况下的声压分布.

以开管的模式数n=5，特征频率为936.47 Hz的情况为例.定义点(0, 0, -500)到点(0, 0, 1 500)的“三维截线”(即中心轴线).以该截线的z轴坐标为横坐标，声压为纵坐标绘制一维声压曲线，结果见图11.为便于比较，图中还画出来一维声波模型预言的声压曲线(正弦函数曲线).由图像可看出，在谐振管内部(0～900 mm)，仿真实验曲线与理论曲线基本重合.但在管口外，仿真曲线迅速衰减，与一维声波模型偏离较大.

图11 中心轴线上声压的仿真结果(实线)与一维模型所得的理论曲线(虚线)

为了进一步探究中心轴线声压的实验曲线与简化理论曲线的偏离细节，在右管口(z=900 mm)附近选取15个截面(xOy平面)，绘制截面上的声压分布图.截面的z坐标范围从893 mm 到907 mm，间隔为1 mm.为了便于通过颜色进行定性比较，将灰度图颜色范围固定为-0.1～1.5 Pa.15个截面上的声压分布如图12所示. 可以看出，在管内部远离开口处，声压均匀分布，符合一维波模型.在开口附近，不论是管内还是管外，声压分布不再均匀，呈现声压随着半径增大而减小的分布规律.

还可以绘制管口附近的声压等值面(图13).同样，在管内距离管口较远的位置，声压等值面基本垂直于中心轴，这意味着一维波模型比较符合实际.但越靠近管口，等压面越向外凸出为弧面(管口类似一个声源)，此时一维波模型不再符合实际情况.

图12 开管时管口附近不同位置的径向声压分布图：共15个截面，z坐标为893～907 mm，间隔1 mm

图13 开管情况管口附近的声压等值面图

2.2 麦克风对谐振频率的影响

在声波谐振管实验中，使用微型麦克风在管内移动以测量各处的声压.麦克风的存在必然会影响管内的声场分布，从而改变谐振频率.这里定量评估麦克风对谐振频率的影响，考虑的影响因素包括麦克风的半径和在管内的位置.为简单起见，只考虑开管情况.

构建圆柱体模拟麦克风，设置半径为3.5 mm，高度为10 mm(真实实验所用麦克风的尺寸).为使麦克风紧贴管内壁，圆柱体底面中心的坐标x=0 mm，y=-13.5 mm.z坐标设为参数，从0 mm到890 mm改变，步长为10 mm.图14画出了在振动模式数n=4、5时，谐振频率的移动量与麦克风的位置的关系(管内未放置麦克风时的谐振频率分别为

749.1 Hz和936.47 Hz).仿真结果显示，麦克风的确会影响谐振频率，但总体上谐振频率的移动量都不大(都不超过1 Hz). 有趣的是，谐振频率的移动也呈正弦周期性变化的规律：在声压波节附近，谐振频率向下移动；而在声压波腹附近，谐振频率向上移动.此外，麦克风在波节处(相比于波腹处)对谐振频率的影响要更大.

n=4

由于麦克风在波腹和波节处对谐振频率的影响最大，进一步探究当麦克风在波腹或波节处时，其半径对谐振频率的影响.把麦克风放置在管中央[底面中心坐标x=0 mm，z=445 mm，y=(-17+r)mm]，其中半径r设为参数，从1 mm到11 mm改变，步长为2 mm.由于开管情况下声压的分布具有对称性，当模式数n为偶数时，麦克风处于波节处；模式数n为奇数时，麦克风处于波腹处.计算不同半径的麦克风引起的谐振频率移动，结果见表2和图15.

表2 存在不同半径麦克风时的谐振频率仿真结果

图15 不同半径麦克风引起的谐振频率移动

仿真结果显示麦克风的半径r越大，谐振频率的移动量越大.当r较小时，频率移动量正比于r2(由于圆柱体高度固定，也正比于麦克风的体积).此外，同在波节(或波腹)处，模式数n越大，谐振频率移动越多.最后一点容易理解：模式数n越大，波长越短，麦克风的大小就越不可以忽略.以上现象在此前的教学中并没有注意.虽然实验中麦克风的半径较小，总体上对谐振频率的影响不大，但这仍是一个值得注意和进一步研究的问题.

2.3 一端部分封闭管的谐振频率

这个实验的意义前文已有说明，下面给出实验结果.表3列出了模式数n=5下的圆形孔、正方形孔、长方形孔和半圆形孔挡片对应谐振频率的实测值(fa)、仿真值(fs).可以看出，绝大部分仿真值与实测值的差在1 Hz以内，而实验测量的重复性误差也在1 Hz左右.这表明，对确定的物理模型，可以用虚拟仿真实验代替实际实验，得到准确的结果.

表3 部分封闭管的谐振频率实验结果与仿真结果

3 总结与讨论

本文介绍了COMSOL数值仿真在声波谐振管实验中一些应用，包括：对管口处声压分布的仿真揭示了管口声波的衰减形式，并给出一维波模型成立的大致范围；仿真麦克风对谐振频率的影响，发现谐振频率移动的方向和大小与麦克风的位置和振动模式数有关；部分封闭管的仿真结果与实测结果符合很好，验证了仿真结果的可靠性.以上结果表明：数值仿真作为真实实验的补充，可以揭示更多的不便实际测量的物理细节，为建立正确的物理图像提供了重要参考.同时，在一些情况下仿真实验可以设计成独立的实验，也可以代替真实的物理实验，得出准确的实验结果，从而大大降低实验成本并缩短实验时间.对于声波谐振管实验，已有利用有限元分析设计仿真实验的教学实践[8]. 此外，落球法测量液体黏性系数、声速测量[9]、动态法测量杨氏模量[10,11]、音叉共振[12]、弦振动[13]等基础物理实验也提供了可以进行数值仿真深入研究的问题.把数值仿真融入基础物理实验教学，可以让学生更深入地理解所涉及的物理问题，同时提高学生解决实际物理问题的能力.