基于边界元法的耳机声振耦合分析及声场优化

林永霖 黎亚军

摘  要：文章针对耳机主被动降噪问题，提出了一种声学边界元与有限元声振耦合算法，构建某型军用耳机的声场分析模型，完整建立耳罩内部各腔体以及振动单体，研究了耳机在不同耳罩结构时的内外声场分布，为耳机被动降噪提供了依据，为了进一步提高主动降噪性能，建立以耳罩声压频响为目标参数，指出了不同麦克风位置时的声压频响差异，该研究结果为耳机主动降噪工程研发提供了一种方法。

关键词：声学边界元;声振耦合;耳机;耳罩;主动降噪;被动降噪;有限元

中图分类号：TB535;O39   文献标识码：A   文章编号：2096-4706（2021）06-0042-06

Acoustic Vibration Coupling Analysis and Sound Field Optimization of

Earphone Based on Boundary Element Method

LIN Yonglin1，LI Yajun2

（1.No.701 Factory of PLA（N），Beijing  100015，China;2.University of Chinese Academy of Sciences，Beijing  100049，China）

Abstract：Aiming at the problem of active and passive noise reduction of earphones，this paper proposes an acoustic boundary element and finite element acoustic vibration coupling algorithm，constructs a sound field analysis model of a certain type of military earphone，and establishes the internal cavities and vibration monomer of the earmuffs completely，and studies the internal and external acoustic field distribution of earphones with different earmuff structures，which provides a basis for passive noise reduction of earphones. In order to further improve the active noise reduction performance，the target parameter of the sound pressure frequency response of the earmuffs is established，and the difference of sound pressure frequency response at different microphone positions is pointed out. The research results provide a method for the research and development of active noise reduction engineering of earphones.

Keywords：acoustic boundary element;acoustic vibration coupling;earphone;earmuff;active noise reduction;passive noise reduction;finite element

0  引  言

边界元法[1]（Boundary element method，BEM）对具有复杂几何构造的无限域问题，会大大降低处理工作量。从20世纪60年代开始，边界元法逐渐应用于结构声学、结构声辐射和声散射问题，它是将HelmholtZ[2]方程边值问题转化为边界积分方程，以格林函数为基本解，利用有限元法[3，4]计算声场域内的物理量[5]。同有限元法相比，边界元法有许多优点：首先，边界元法将流体域内的计算转化到边界上，使问题的维数降低了一维，从而减少了问题的自由度和原始信息量;其次，由于利用的微分方程的解析基本解作为边界积分方程的核函数，具有半解析半数值方法的特点，因而具有较高的精度;最后，对于无限域或半无限域问题，边界元法十分适合，无须在远场边界离散，所有计算都在结构表面进行，大大减小了计算域。由于边界元法的这些优点，边界元法在结构声学计算的数值方法中占据主导地位。结构声學边界元法还可分为直接边界元法[6]（Direct Boundary Element Method，DBEM）和间接边界元法[7]（Indirect Boundary Element Method，IBEM）。目前已有多种方法来克服基于边界元法外部解的非唯一性问题[8，9]，主要有两种方法，一是CHIEF方法[10-12]，但对于复杂模型和高频问题，CHIEF点很难选择。另一种方法是Burton-Miller方法[13]，它是由边界积分方程及其对源点法向导数方程线性组合而成，由于其便于构造并已在许多问题下获得验证，逐渐被广泛采用。基于Burton-Miller法的多重边界积分法（Dual BEM）[14]会导致解的奇异性。针对不同单元类型，目前已有很多方法，如Hadamard有限项积分[15]、高阶Galerkin方法[16]、核函数相减方法[17]及场相减方法等。

本文首先简述了边界元的基本理论，然后分析了边界元与有限元声振耦合算法[18，19]基本过程，此后基于此方法构建了某型军用耳机的声场分析模型，分析了耳机在不同耳罩结构时的内外声场分布，研究结果为类似工程问题提供了参考依据。

1  BEM基本理论

满足理想流体介质假设下声波波动方程为：

（1）

式中?2表示Laplace算子，p（r，t）表示t时刻位于域中的声压，Ωa表示流体域，对于稳态简谐声场，p（r，t）= p（r）ejωt，由此可得Helmholtz方程：

?2p（r）+k2p（r）=0                         （2）

式中k=ω/c为复波数，ω为角频率，c为声速。当存在外部源，波动方程满足：

（3）

为源强，j表示虚数单位，点声源作为声场中最典型的激励源，可用格林函数表示，G（x，y）为三维自由空间的格林函数，满足以下：

（4）

纯声场域内声学边界条件为三类即声压、法向速度和声阻抗边界;外场除了满足三类基本边界条件外，还要满足无限域的sommerfeld边界[20]。

将格林函数代入格林公式的第二类表达式，从而得到直接边界积分方程为：

（5）

为积分点，为源点，C（x）为积分点位置函数，n为单元法线方向，当x在域内、光滑边界上和域外分别取1、1/2和0。

对于开口边界或者同时求解内部和外部声场的情况，声学域在边界的内外同时存在，需要同时考虑内外边界，从而得到：

（6）

当  在Ωa上时为超奇异积分，采用Hadamard原理降低奇异性。

将边界离散为单元，在每个小单元内通过单元形函数建立单元内的变量与单元尺寸的关系，生成了单元代数方程组，然后由全局形函数，节点i属于哪一个单元，全局形函数就是那个节点的形函数，其他节点为0，建立单元求解，得到：

（7）

进一步得到全局形函数：

（8）

为1×na维，  和  为na×1未知节点自由度。对于未知自由度的求解，将全局形函数代入边界积分方程中，遍及边界上的所有节点，得到总体方程组：

（9）

[A]=[B]为na×na维，方程为na代数方程组，具有2na个未知自由度，由边界条件可得每个节点的预知声压或法向速度或者法向声阻抗率，因此可以求得na个未知数。对于场点的自由度求解，直接利用已求得的边界值，代入边界积分方程中可得：

（10）

间接边界元建模过程与直接边界元类似，不同的是间接边界元单元首先求解边界上的单层势和双层势，由于在某些位置上间接边界元的求解需要数值计算hadamard有限积分项，使得直接使用边界积分方程建模困难，为此，用一种变量公式，壁面超奇异积分，确保了边界单元方程矩阵的对称性。

将离散化后的边界模型代入基于双层势和单层势建立的积分方程中，可得到关于未知量的边界单元积分方程：

（11）

其中：

，

接下来，同直接边界元，对于域内不属于边界Ωa上任意场点的未知量，直接由边界上的已知量代入间接边界积分方程求得：

（12）

其中：

，

边界单元模型的计算精度主要依赖于单元个数和形函数，形函数的完整性和对称性。

内部声振耦合边界条件除了原来的三种，声学边界包含了弹性结构边界，还要保证表面的法向速度连续条件，即结构与流体耦合界面处的流体法向速度等于结构法向速度：

（13）

由声场压力引起结构的载荷变化，得到耦合结构动力学方程：

（14）

其中，[K]为刚度矩阵;[C]指阻尼矩阵;[M]指质量矩阵耦合矩阵为ns×na1：

（15）

相应地，结构运动位移引起耦合面上声域法向质点速度变化。从而得到耦合声场方程：

（16）

结合耦合动力学方程和声场边界积分方程可得耦合总体方程：

（17）

其中：

，

对于有限元与间接边界元的耦合，同理可将边界元表面看作两部分，一部分与结构有限元耦合nμ1，另一部分为已知法向速度nμ2，边界条件和直接边界元耦合边界条件已知，在耦合面上法向速度相等。结构有限元模型和声学间接边界元模型：

（18）

（19）

声学间接边界元方程：

（20）

其中：

，

，

合并結构动力学方程和间接边界元方程可得声振耦合方程：

（21）

其中：

，

特别的，当边界元表面全部由弹性结构所组成时，耦合方程变为：

（22）

其中：

综上，耦合有限元比边界元在计算内场时有更高的效率，比如汽车和航空，但是对于像扬声器设计和隔声墙之类的外场或开口声振耦合问题，边界元耦合比有限元更具优势，因为它自动满足外场sommerfeld远程辐射边界条件，在有限元中只能用一种近似的方式来模拟。

2  耳罩模型的数值验证

2.1  问题分析

主动降噪耳机作为近年来发展快速的声学器件，对其频响和三维声场的准确模拟一直以来都是研究的重点，Borwick等总结了耳机结构类比集总参数模型，给出了不同形状腔体和导管的仿真方法，Slotte等人在此基础上优化研究，但高频段由于结构尺寸不再远小于声波波长，集总参数模型仿真误差增大，另外由于耳机腔体内部不规则，耳内声场并不是单向行进的平面波，故传输线模型也有局限。

为此，本文提出了基于声学边界元的结构声振耦合分析，完整建立耳罩内部各腔体以及振动单体，振动部分采用有限元求解保证频响精度，然后基于声振耦合理论建立声场边界，进而可以获得宽频空间声场分布特性。

2.2  参数设置

建立耳罩式耳机模型如图1所示，计算平台为个人PC，配置为Intel Core（TM）i7-3840QM@2.8 GHz和RAM 16 GB。

将模型离散为三角形常单元，最大单元长度0.02 mm，可计算最高频率8 500 Hz，根据腔体结构组成，分析三个工况，分别为前腔+扬声器、前腔+内腔+扬声器、前腔+内腔+后腔+扬声器。主要求解参数如表1所示。

3  结果分析

3.1  结构模态

声振耦合分析首先需要求解结构有限元模型的振动参数，为此，本模型中以扬声器为声振耦合对象，分析了振膜的振动模态，表2为前10阶固有频率，图2为前6阶振型，根据理论解可以求得解析解，，看出与理论解非常一致。

3.2  声场分布

结构边界条件为振膜外沿固定约束，声学边界条件为在耳机正前方1 m处平面波，p（x）=pae-jkx幅值为1 Pa，初相位为0，计算声波传播以及与振膜耦合后总声场。分三种工况（如前述）进行计算，计算结束后可获得声场中任一点处的声压，为便于分析，提取耳机中心对称面上的内外声场云图。

另外，基于人耳对声波不同频率的响应灵敏度不同，依次选取为500 Hz、1 000 Hz、2 000 Hz、3 000 Hz的4个最典型频率。在只有前腔和内腔时，随着频率升高，耳罩壳体对声场的散射作用逐渐增强，前腔主要影响耳外特别是出声孔前方的声场，内腔主要影响耳内扬声器振膜背部的声场。在500～1 000 Hz两个较低频率处，声场比较均匀，声压最高值主要集中在扬声器振膜背部，但后两个频率处声压最高值主要集中在内腔，说明随着频率增高，声振耦合作用增强，振膜背部的振动使得后腔有限容积的声场共振，特别是2 000 Hz处，内腔共振非常明显，而3 000 Hz处声场指向性增强。

随着频率升高壳体对声场散射影响增强，前2个频率处外声场比较均匀，主要在扬声器振膜背部声压较大，后两个频率处声场指向性增强。不同的是相比于工况1，工况3下由于后壳为凸型结构，有利用声扩散从而使得声场最大声压由明显的降低，被动降噪效果显著，另外，由于后腔容积增大，使得耳腔内的共振频率向低频移动，相比于工况1，2 000 Hz处的声压明显降低23 dB，同时1 000 Hz处由2 dB增大，除此之外，还可以得出，由于后腔的不规则形状，使得耳内和耳外的声场都呈现出不对称分布。

在后腔和内腔共同作用下，更加复杂的声场分布，呈现兼具前两种工况的各自特点，基于工况2，由于增加了内腔，后腔被分解为两个腔体，相当与两个声容并联，使得腔内声场向高频转移，低频段声场分布更具工况2的特点，高频段分布更具工况1的特点，在500 Hz时由于增大了背腔，从而内腔声压级有一定降低，外场声压分布都近似为94 dB，在1 000 Hz时，由于外腔内部结构引起类似于工况2的内声场不均匀性，外声场依然均匀分布;在2 000 Hz时内腔的影响，使得内场声压有一定增高，同时出声孔处的声场降低;在3 000 Hz时由于使得耳内声压相比工况2增大5 dB，但相比工况1因后腔的影响而有2 dB的下降。

4  频响分析

为了分析人耳接收到的声压特性，提取三种工况下耳机出声孔中心点的声压频响如图3所示，结合上述声压分布云图，可以看出，在耳外声压频响谷值处恰为耳内声振耦合的共振点，对于只有内腔的耳外频响，在2 000 Hz为最小，耳内声压为最大，在后腔1 000 Hz附近耳外声压谷值处耳内声压同样最大，同时双腔时谷值进一步向低频移动。考虑到耳腔结构的不规则，将引起耳外出声孔处的声压不均匀，为此，将出声孔表面的声压平均，得到三种工况下的平均声压频响如圖4所示，可以看出，平均后的声压更能说明声场特性，也更接近人耳实际听到的声压，腔体对2 500 Hz以后的声压具有被动降噪作用，双腔的降噪效果大于单独内腔或后腔，在1 500～2 500 Hz频段降噪效果具有一定复杂性，在频段端两种腔体基本降噪效果都不明显，具体降噪频响如图5所示，可以看出，双腔体比单独内腔降噪效果明显，最大可达6 dB，在高频段平均降低1 dB。

5  结  论

本文研究了基于Burton-Miller法间接边界元（BEM）声振耦合分析在同时计算内外腔以及开口声学边界下的应用，对某型军用耳罩式耳机的腔体声场分布及人耳处频响进行了分析，分析说明了双腔体外壳相比于单壳能更好的起到被动降噪的作用，内腔的存在使得高频处的耳腔增加，外腔的存在有利于降低耳内共振，并使高频处的被动降噪增加，为了被动降噪效果，应该尽可能增加内腔容积或增加吸声材料，同时进行内腔结构设计优化，以降低人耳声压频响的不均匀性。本研究结果说明了间接边界元法在求解内外声场的有效性，为耳机声学研发提供了方法指导，关于腔体与人耳耦合以及更详细的声腔优化需要进一步的研究。

参考文献：

[1] MASUMOTO T，YASUDA Y，INOUE N，et al. Fast Calculation of Far-Field Sound Directivity Based on Fast Multipole Boundary Element Method [J].Journal of Theoretical and Computational Acoustics，2020，28（4）.

[2] GILVEY B，TREVELYAN J，HATTORI G. Singular enrichment functions for Helmholtz scattering at corner locations using the boundary element method [J]. International Journal for Numerical Methods in Engineering，2020，121（3）：519-533.

[3] 魯建锋，许孝卓，高岩，等.磁阻式悬浮平台侧向气隙控制研究 [J].电子科技，2020，33（10）：15-20+63.

[4] 许伯强，来锴，徐桂东，等.复合层合梁层裂损伤与超声导波的相互作用研究 [J].电子科技，2015，28（11）：55-60.

[5] 李俞霖，向宇，伍松.基于声全息技术的噪声源识别与DSP实现 [J].电子科技，2017，30（4）：27-31.

[6] BRITO W K F，MAIA C D C D，MENDONCA A V. Bending analysis of elastically connected Euler–Bernoulli double-beam system using the direct boundary element method [J].Applied Mathematical Modelling，2019，74：387-408.

[7] HAERI H. Experimental and numerical study on crack propagation in pre-cracked beam specimens under three-point bending [J].Journal of Central South University，2016，23（2）：430-439.

[8] LIU Y J，RIZZO F J. A weakly singular form of the hypersingular boundary integral equation applied to 3-D acoustic wave problems [J].North-Holland，1992，96（2）：271–287.

[9] KRESS R. On the condition number of boundary integral operators in acoustic and electromagnetic scattering [J].The Quarterly Journal of Mechanics and Applied Mathematics，1985，38（2）：323-341.

[10] LIU Y J. Analysis of shell-like structures by the Boundary Element Method based on 3-D elasticity：formulation and verification [J].International Journal for Numerical Methods in Engineering，1998，41（3）：541-558.

[11] 庞岳峰，马占顺，张佳宁，等.基于修正最小二乘法的测控设备同站引导算法 [J].电子科技，2019，32（3）：82-86.

[12] OTHMAN N Y，SALEH Z A，OMRAN Z A. Development of Stage–Distance–Discharge Relationship and Rating Curve using Least Square Method [J].Civil Engineering Journal，2019，5（9）：1959-1969.

[13] ZHANG J M，LIN W C，SHU X M，et al. A dual interpolation boundary face method for exterior acoustic problems based on the Burton–Miller formulation [J].Engineering Analysis with Boundary Elements，2020，113（C）：219-231.

[14] KAO J H，CHEN K H，CHEN J T，et al. Isogeometric Analysis of the Dual Boundary Element Method for the Laplace Problem With a Degenerate Boundary [J].Journal of Mechanics，2020，36（1）：35-46.

[15] FENG H，GAO Y，JYU L L，et al. A new collocation method for solving certain Hadamard finite-part integral equation [J].International Journal of Numerical Analysis and Modeling，2018，16（2）：240-254.

[16] FEDOROV V E. Error estimate of the nonstationary Galerkin method for a higher order equation with changing time direction [C]//AIP Conference Proceedings，2017，1907（1）：1-4.

[17] LI K M，WANG Y M. Rapid computation of waves diffracted by an absorbing plane with precise error bounds [J].Journal of the Acoustical Society of America，2016，139（4）：2010-2010.

[18] WANG Q，CHANG W D，LIU S H，et al. Study on noise-vibration coupling characteristics of premixed methane–air flame propagation in a tube with an acoustic absorption material [J].RSC Advances，2019，9（49）：28323-28329.

[19] 张冠军，朱翔，李天匀，等.双层加筋板水下声振耦合特性研究 [J].船舶力学，2019，23（1）：78-87.

[20] LUO W H，HUANG T Z，LI L，et al. Quadratic spline collocation method and efficient preconditioner for the Helmholtz equation with the Sommerfeld boundary conditions [J].Japan Journal of Industrial and Applied Mathematics，2016，33（3）：701-720.

作者简介：林永霖（1979—），男，汉族，福建漳州人，工程师，硕士，研究方向：信号处理、电子信息;黎亚军（1986—），男，汉族，甘肃灵台人，工程师，博士，研究方向：有源噪声控制。