基于邊界元法的耳機聲振耦合分析及聲場優(yōu)化

林永霖 黎亞軍

摘 ?要：文章針對耳機主被動降噪問題，提出了一種聲學邊界元與有限元聲振耦合算法，構(gòu)建某型軍用耳機的聲場分析模型，完整建立耳罩內(nèi)部各腔體以及振動單體，研究了耳機在不同耳罩結(jié)構(gòu)時的內(nèi)外聲場分布，為耳機被動降噪提供了依據(jù)，為了進一步提高主動降噪性能，建立以耳罩聲壓頻響為目標參數(shù)，指出了不同麥克風位置時的聲壓頻響差異，該研究結(jié)果為耳機主動降噪工程研發(fā)提供了一種方法。

關(guān)鍵詞：聲學邊界元;聲振耦合;耳機;耳罩;主動降噪;被動降噪;有限元

中圖分類號：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）對具有復雜幾何構(gòu)造的無限域問題，會大大降低處理工作量。從20世紀60年代開始，邊界元法逐漸應用于結(jié)構(gòu)聲學、結(jié)構(gòu)聲輻射和聲散射問題，它是將HelmholtZ[2]方程邊值問題轉(zhuǎn)化為邊界積分方程，以格林函數(shù)為基本解，利用有限元法[3，4]計算聲場域內(nèi)的物理量[5]。同有限元法相比，邊界元法有許多優(yōu)點：首先，邊界元法將流體域內(nèi)的計算轉(zhuǎn)化到邊界上，使問題的維數(shù)降低了一維，從而減少了問題的自由度和原始信息量;其次，由于利用的微分方程的解析基本解作為邊界積分方程的核函數(shù)，具有半解析半數(shù)值方法的特點，因而具有較高的精度;最后，對于無限域或半無限域問題，邊界元法十分適合，無須在遠場邊界離散，所有計算都在結(jié)構(gòu)表面進行，大大減小了計算域。由于邊界元法的這些優(yōu)點，邊界元法在結(jié)構(gòu)聲學計算的數(shù)值方法中占據(jù)主導地位。結(jié)構(gòu)聲學邊界元法還可分為直接邊界元法[6]（Direct Boundary Element Method，DBEM）和間接邊界元法[7]（Indirect Boundary Element Method，IBEM）。目前已有多種方法來克服基于邊界元法外部解的非唯一性問題[8，9]，主要有兩種方法，一是CHIEF方法[10-12]，但對于復雜模型和高頻問題，CHIEF點很難選擇。另一種方法是Burton-Miller方法[13]，它是由邊界積分方程及其對源點法向?qū)?shù)方程線性組合而成，由于其便于構(gòu)造并已在許多問題下獲得驗證，逐漸被廣泛采用?；贐urton-Miller法的多重邊界積分法（Dual BEM）[14]會導致解的奇異性。針對不同單元類型，目前已有很多方法，如Hadamard有限項積分[15]、高階Galerkin方法[16]、核函數(shù)相減方法[17]及場相減方法等。

本文首先簡述了邊界元的基本理論，然后分析了邊界元與有限元聲振耦合算法[18，19]基本過程，此后基于此方法構(gòu)建了某型軍用耳機的聲場分析模型，分析了耳機在不同耳罩結(jié)構(gòu)時的內(nèi)外聲場分布，研究結(jié)果為類似工程問題提供了參考依據(jù)。

1 ?BEM基本理論

滿足理想流體介質(zhì)假設(shè)下聲波波動方程為：

（1）

式中?2表示Laplace算子，p（r，t）表示t時刻位于域中的聲壓，Ωa表示流體域，對于穩(wěn)態(tài)簡諧聲場，p（r，t）= p（r）ejωt，由此可得Helmholtz方程：

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

式中k=ω/c為復波數(shù)，ω為角頻率，c為聲速。當存在外部源，波動方程滿足：

（3）

為源強，j表示虛數(shù)單位，點聲源作為聲場中最典型的激勵源，可用格林函數(shù)表示，G（x，y）為三維自由空間的格林函數(shù)，滿足以下：

（4）

純聲場域內(nèi)聲學邊界條件為三類即聲壓、法向速度和聲阻抗邊界;外場除了滿足三類基本邊界條件外，還要滿足無限域的sommerfeld邊界[20]。

將格林函數(shù)代入格林公式的第二類表達式，從而得到直接邊界積分方程為：

（5）

為積分點，為源點，C（x）為積分點位置函數(shù)，n為單元法線方向，當x在域內(nèi)、光滑邊界上和域外分別取1、1/2和0。

對于開口邊界或者同時求解內(nèi)部和外部聲場的情況，聲學域在邊界的內(nèi)外同時存在，需要同時考慮內(nèi)外邊界，從而得到：

（6）

當 ?在Ωa上時為超奇異積分，采用Hadamard原理降低奇異性。

將邊界離散為單元，在每個小單元內(nèi)通過單元形函數(shù)建立單元內(nèi)的變量與單元尺寸的關(guān)系，生成了單元代數(shù)方程組，然后由全局形函數(shù)，節(jié)點i屬于哪一個單元，全局形函數(shù)就是那個節(jié)點的形函數(shù)，其他節(jié)點為0，建立單元求解，得到：

（7）

進一步得到全局形函數(shù)：

（8）

為1×na維， ?和 ?為na×1未知節(jié)點自由度。對于未知自由度的求解，將全局形函數(shù)代入邊界積分方程中，遍及邊界上的所有節(jié)點，得到總體方程組：

（9）

[A]=[B]為na×na維，方程為na代數(shù)方程組，具有2na個未知自由度，由邊界條件可得每個節(jié)點的預知聲壓或法向速度或者法向聲阻抗率，因此可以求得na個未知數(shù)。對于場點的自由度求解，直接利用已求得的邊界值，代入邊界積分方程中可得：

（10）

間接邊界元建模過程與直接邊界元類似，不同的是間接邊界元單元首先求解邊界上的單層勢和雙層勢，由于在某些位置上間接邊界元的求解需要數(shù)值計算hadamard有限積分項，使得直接使用邊界積分方程建模困難，為此，用一種變量公式，壁面超奇異積分，確保了邊界單元方程矩陣的對稱性。

將離散化后的邊界模型代入基于雙層勢和單層勢建立的積分方程中，可得到關(guān)于未知量的邊界單元積分方程：

（11）

其中：

，

接下來，同直接邊界元，對于域內(nèi)不屬于邊界Ωa上任意場點的未知量，直接由邊界上的已知量代入間接邊界積分方程求得：

（12）

其中：

，

邊界單元模型的計算精度主要依賴于單元個數(shù)和形函數(shù)，形函數(shù)的完整性和對稱性。

內(nèi)部聲振耦合邊界條件除了原來的三種，聲學邊界包含了彈性結(jié)構(gòu)邊界，還要保證表面的法向速度連續(xù)條件，即結(jié)構(gòu)與流體耦合界面處的流體法向速度等于結(jié)構(gòu)法向速度：

（13）

由聲場壓力引起結(jié)構(gòu)的載荷變化，得到耦合結(jié)構(gòu)動力學方程：

（14）

其中，[K]為剛度矩陣;[C]指阻尼矩陣;[M]指質(zhì)量矩陣耦合矩陣為ns×na1：

（15）

相應地，結(jié)構(gòu)運動位移引起耦合面上聲域法向質(zhì)點速度變化。從而得到耦合聲場方程：

（16）

結(jié)合耦合動力學方程和聲場邊界積分方程可得耦合總體方程：

（17）

其中：

，

對于有限元與間接邊界元的耦合，同理可將邊界元表面看作兩部分，一部分與結(jié)構(gòu)有限元耦合nμ1，另一部分為已知法向速度nμ2，邊界條件和直接邊界元耦合邊界條件已知，在耦合面上法向速度相等。結(jié)構(gòu)有限元模型和聲學間接邊界元模型：

（18）

（19）

聲學間接邊界元方程：

（20）

其中：

，

，

合并結(jié)構(gòu)動力學方程和間接邊界元方程可得聲振耦合方程：

（21）

其中：

，

特別的，當邊界元表面全部由彈性結(jié)構(gòu)所組成時，耦合方程變?yōu)椋?/p>

（22）

其中：

綜上，耦合有限元比邊界元在計算內(nèi)場時有更高的效率，比如汽車和航空，但是對于像揚聲器設(shè)計和隔聲墻之類的外場或開口聲振耦合問題，邊界元耦合比有限元更具優(yōu)勢，因為它自動滿足外場sommerfeld遠程輻射邊界條件，在有限元中只能用一種近似的方式來模擬。

2 ?耳罩模型的數(shù)值驗證

2.1 ?問題分析

主動降噪耳機作為近年來發(fā)展快速的聲學器件，對其頻響和三維聲場的準確模擬一直以來都是研究的重點，Borwick等總結(jié)了耳機結(jié)構(gòu)類比集總參數(shù)模型，給出了不同形狀腔體和導管的仿真方法，Slotte等人在此基礎(chǔ)上優(yōu)化研究，但高頻段由于結(jié)構(gòu)尺寸不再遠小于聲波波長，集總參數(shù)模型仿真誤差增大，另外由于耳機腔體內(nèi)部不規(guī)則，耳內(nèi)聲場并不是單向行進的平面波，故傳輸線模型也有局限。

為此，本文提出了基于聲學邊界元的結(jié)構(gòu)聲振耦合分析，完整建立耳罩內(nèi)部各腔體以及振動單體，振動部分采用有限元求解保證頻響精度，然后基于聲振耦合理論建立聲場邊界，進而可以獲得寬頻空間聲場分布特性。

2.2 ?參數(shù)設(shè)置

建立耳罩式耳機模型如圖1所示，計算平臺為個人PC，配置為Intel Core（TM）i7-3840QM@2.8 GHz和RAM 16 GB。

將模型離散為三角形常單元，最大單元長度0.02 mm，可計算最高頻率8 500 Hz，根據(jù)腔體結(jié)構(gòu)組成，分析三個工況，分別為前腔+揚聲器、前腔+內(nèi)腔+揚聲器、前腔+內(nèi)腔+后腔+揚聲器。主要求解參數(shù)如表1所示。

3 ?結(jié)果分析

3.1 ?結(jié)構(gòu)模態(tài)

聲振耦合分析首先需要求解結(jié)構(gòu)有限元模型的振動參數(shù)，為此，本模型中以揚聲器為聲振耦合對象，分析了振膜的振動模態(tài)，表2為前10階固有頻率，圖2為前6階振型，根據(jù)理論解可以求得解析解，，看出與理論解非常一致。

3.2 ?聲場分布

結(jié)構(gòu)邊界條件為振膜外沿固定約束，聲學邊界條件為在耳機正前方1 m處平面波，p（x）=pae-jkx幅值為1 Pa，初相位為0，計算聲波傳播以及與振膜耦合后總聲場。分三種工況（如前述）進行計算，計算結(jié)束后可獲得聲場中任一點處的聲壓，為便于分析，提取耳機中心對稱面上的內(nèi)外聲場云圖。

另外，基于人耳對聲波不同頻率的響應靈敏度不同，依次選取為500 Hz、1 000 Hz、2 000 Hz、3 000 Hz的4個最典型頻率。在只有前腔和內(nèi)腔時，隨著頻率升高，耳罩殼體對聲場的散射作用逐漸增強，前腔主要影響耳外特別是出聲孔前方的聲場，內(nèi)腔主要影響耳內(nèi)揚聲器振膜背部的聲場。在500～1 000 Hz兩個較低頻率處，聲場比較均勻，聲壓最高值主要集中在揚聲器振膜背部，但后兩個頻率處聲壓最高值主要集中在內(nèi)腔，說明隨著頻率增高，聲振耦合作用增強，振膜背部的振動使得后腔有限容積的聲場共振，特別是2 000 Hz處，內(nèi)腔共振非常明顯，而3 000 Hz處聲場指向性增強。

隨著頻率升高殼體對聲場散射影響增強，前2個頻率處外聲場比較均勻，主要在揚聲器振膜背部聲壓較大，后兩個頻率處聲場指向性增強。不同的是相比于工況1，工況3下由于后殼為凸型結(jié)構(gòu)，有利用聲擴散從而使得聲場最大聲壓由明顯的降低，被動降噪效果顯著，另外，由于后腔容積增大，使得耳腔內(nèi)的共振頻率向低頻移動，相比于工況1，2 000 Hz處的聲壓明顯降低23 dB，同時1 000 Hz處由2 dB增大，除此之外，還可以得出，由于后腔的不規(guī)則形狀，使得耳內(nèi)和耳外的聲場都呈現(xiàn)出不對稱分布。

在后腔和內(nèi)腔共同作用下，更加復雜的聲場分布，呈現(xiàn)兼具前兩種工況的各自特點，基于工況2，由于增加了內(nèi)腔，后腔被分解為兩個腔體，相當與兩個聲容并聯(lián)，使得腔內(nèi)聲場向高頻轉(zhuǎn)移，低頻段聲場分布更具工況2的特點，高頻段分布更具工況1的特點，在500 Hz時由于增大了背腔，從而內(nèi)腔聲壓級有一定降低，外場聲壓分布都近似為94 dB，在1 000 Hz時，由于外腔內(nèi)部結(jié)構(gòu)引起類似于工況2的內(nèi)聲場不均勻性，外聲場依然均勻分布;在2 000 Hz時內(nèi)腔的影響，使得內(nèi)場聲壓有一定增高，同時出聲孔處的聲場降低;在3 000 Hz時由于使得耳內(nèi)聲壓相比工況2增大5 dB，但相比工況1因后腔的影響而有2 dB的下降。

4 ?頻響分析

為了分析人耳接收到的聲壓特性，提取三種工況下耳機出聲孔中心點的聲壓頻響如圖3所示，結(jié)合上述聲壓分布云圖，可以看出，在耳外聲壓頻響谷值處恰為耳內(nèi)聲振耦合的共振點，對于只有內(nèi)腔的耳外頻響，在2 000 Hz為最小，耳內(nèi)聲壓為最大，在后腔1 000 Hz附近耳外聲壓谷值處耳內(nèi)聲壓同樣最大，同時雙腔時谷值進一步向低頻移動?？紤]到耳腔結(jié)構(gòu)的不規(guī)則，將引起耳外出聲孔處的聲壓不均勻，為此，將出聲孔表面的聲壓平均，得到三種工況下的平均聲壓頻響如圖4所示，可以看出，平均后的聲壓更能說明聲場特性，也更接近人耳實際聽到的聲壓，腔體對2 500 Hz以后的聲壓具有被動降噪作用，雙腔的降噪效果大于單獨內(nèi)腔或后腔，在1 500～2 500 Hz頻段降噪效果具有一定復雜性，在頻段端兩種腔體基本降噪效果都不明顯，具體降噪頻響如圖5所示，可以看出，雙腔體比單獨內(nèi)腔降噪效果明顯，最大可達6 dB，在高頻段平均降低1 dB。

5 ?結(jié) ?論

本文研究了基于Burton-Miller法間接邊界元（BEM）聲振耦合分析在同時計算內(nèi)外腔以及開口聲學邊界下的應用，對某型軍用耳罩式耳機的腔體聲場分布及人耳處頻響進行了分析，分析說明了雙腔體外殼相比于單殼能更好的起到被動降噪的作用，內(nèi)腔的存在使得高頻處的耳腔增加，外腔的存在有利于降低耳內(nèi)共振，并使高頻處的被動降噪增加，為了被動降噪效果，應該盡可能增加內(nèi)腔容積或增加吸聲材料，同時進行內(nèi)腔結(jié)構(gòu)設(shè)計優(yōu)化，以降低人耳聲壓頻響的不均勻性。本研究結(jié)果說明了間接邊界元法在求解內(nèi)外聲場的有效性，為耳機聲學研發(fā)提供了方法指導，關(guān)于腔體與人耳耦合以及更詳細的聲腔優(yōu)化需要進一步的研究。

參考文獻：

[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] 魯建鋒，許孝卓，高巖，等.磁阻式懸浮平臺側(cè)向氣隙控制研究 [J].電子科技，2020，33（10）：15-20+63.

[4] 許伯強，來鍇，徐桂東，等.復合層合梁層裂損傷與超聲導波的相互作用研究 [J].電子科技，2015，28（11）：55-60.

[5] 李俞霖，向宇，伍松.基于聲全息技術(shù)的噪聲源識別與DSP實現(xiàn) [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] 龐岳峰，馬占順，張佳寧，等.基于修正最小二乘法的測控設(shè)備同站引導算法 [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—），男，漢族，甘肅靈臺人，工程師，博士，研究方向：有源噪聲控制。