快速多極子邊界元法在吸聲材料聲場計算中的應(yīng)用

崔曉兵，季振林

(哈爾濱工程大學(xué) 動力與能源工程學(xué)院，哈爾濱 150001)

自 Rokhlin［1］于 1983年提出快速多極子法則(FMA)以來，F(xiàn)MA在分子動力學(xué)、電磁學(xué)、聲學(xué)、流體力學(xué)等領(lǐng)域得到了迅速發(fā)展［2－7］。在聲學(xué)方面，將其與邊界元法(BEM)相結(jié)合，形成了快速多極子邊界元法(FMBEM)，鑒于其在求解聲場問題中計算速度快，精度高，節(jié)約內(nèi)存等優(yōu)勢，使其成為處理大尺度聲場問題的有效工具，為邊界元法的發(fā)展帶來了新曙光。

目前FMBEM的應(yīng)用多局限于空氣介質(zhì)的單區(qū)域聲場計算，其在多區(qū)域多介質(zhì)復(fù)合聲場問題中的應(yīng)用成為聲學(xué)FMBEM發(fā)展的新方向，將子結(jié)構(gòu)技術(shù)應(yīng)用于FMBEM中，發(fā)展形成子結(jié)構(gòu)FMBEM［8］，能有效地解決多區(qū)域復(fù)合聲場問題，使帶有薄壁結(jié)構(gòu)的聲場計算(如內(nèi)插管消聲器，穿孔管消聲器等)成為可能。然而，對于FMBEM在吸聲材料介質(zhì)或多介質(zhì)復(fù)合的聲場中的應(yīng)用，由于復(fù)數(shù)形式波數(shù)的影響，其多極子展開式的適用性及求解精度還有待考察，復(fù)波數(shù)條件下展開式數(shù)值計算的參數(shù)選取仍需進一步研究與分析。

鑒于此，本文對復(fù)波數(shù)及實波數(shù)下的格林函數(shù)及其法向?qū)?shù)多極子展開式進行數(shù)值研究，建立四點單級和多極傳遞關(guān)系模型，通過與理論值比較，考察其計算精度，分析其誤差的影響因素及來源，并找到解決此問題的方法。最后以膨脹腔阻性消聲器傳遞損失計算為例，比較FMBEM與BEM的求解精度，驗證本文方法的有效性與可行性。

1 復(fù)波數(shù)條件下的多極子展開式

1.1 邊界積分方程

圖1為邊界平滑的三維內(nèi)部聲場示意圖，包含剛性邊界、振動邊界、吸收邊界三種邊界。假設(shè)穩(wěn)態(tài)聲場聲壓隨時間變化是簡諧的，取時間項為exp(－ jωt)，其中 j=，則該聲場滿足三維Helmholtz方程，邊界上P點處聲壓可由Kirchhoff-Helmholtz邊界積分方程［9］得到:

圖1 平滑邊界內(nèi)部聲場示意圖Fig.1 Illustration of an interior sound field with smooth boundary

其中，?/?n指該函數(shù)的外法向?qū)?shù)，G為Helmholtz方程的基本解，形式如下:

眾所周知，式(1)即為傳統(tǒng)BEM的核心控制方程，將該控制方程離散，通過某特定方法(1/8子組劃分法，1/2子組劃分法等)對邊界節(jié)點進行分組，分級，按照每級中各組集的遠近關(guān)系將聲場分為近場與遠場，在近場仍用傳統(tǒng)BEM法求解，遠場則采用快速多極子方法加速求解［2］，最后將兩部分結(jié)果相加即可得到整個聲場的系數(shù)矩陣和向量積。此種求解聲場問題的方法即為FMBEM(具體實現(xiàn)過程可參考文獻［10］)。

1.2 Helmholtz方程基本解的多極子展開式

為了在遠場應(yīng)用 FMA求解，需將邊界積分方程等效為多極子展開式進行計算，而Helmholtz方程基本解即格林函數(shù)的多極子展開是FMBEM順利實施的前提，其展開式的計算精度尤其成為FMBEM成功的關(guān)鍵。

如圖2所示四點關(guān)系，M為多極子展開點，L為本地展開點。根據(jù)Gegenbauer附加定理［11］及平面波展開公式［12］，格林函數(shù)可展開成:

圖2 四點關(guān)系圖Fig.2 Geometry of the four points

其中:

為保證式(3)成立，要求 rpL＜rLM且 rMQ＜rPM。，為波數(shù)矢量是單位球面積分向量為第一類球漢克函數(shù)，pl為勒讓德多項式。聲傳播介質(zhì)為空氣時，k=2πf/c0為實波數(shù)，當傳播介質(zhì)為吸聲材料時，假設(shè)吸聲材料介質(zhì)分布均勻，且聲波以簡諧波的形式傳播，則在該介質(zhì)中，密度與聲速均可等效為復(fù)數(shù)形式的聲參數(shù)，從而波數(shù)k=kr+j ki為復(fù)波數(shù)，代入式(3)可得:

格林函數(shù)外法向?qū)?shù)為:

由于:

其外法向?qū)?shù)可展開為:

在進行遠場多級影響系數(shù)計算時，根據(jù)快速多極子算法，格林函數(shù)需由下式計算:

式中:

其中，L為最低級級數(shù)(級數(shù)最高)，λmL為L級m組的中心點;λm'L為L級m組的某交互組的中心點(關(guān)于近場組、交互組等概念參見文獻［10］);級數(shù)I由p點、q點位置關(guān)系決定。由于2級之后才有交互組出現(xiàn)，因而上述兩式適用于L≥3的情況。同理格林函數(shù)法向?qū)?shù)的多極表達式為:

1.3 多極子展開式的數(shù)值計算

為了進行有效的數(shù)值計算，式(5)需用前Nc項和近似，考慮到復(fù)波數(shù)的影響，為使展開式在各種情況下得到足夠的精度，Nc可由下式計算:

可見，Nc及Nθ的取值直接影響著多極子展開式的求解精度。在空氣介質(zhì)條件下，即實波數(shù)條件下，Yasuda和Sakuma［14］對 Nc及 Nθ等參數(shù)的選取做了詳細的數(shù)值研究與論證，并提出了恰當?shù)慕?jīng)驗公式以保證展開式的計算誤差在允許的范圍。然而，當聲波在吸聲材料介質(zhì)中傳播時，復(fù)波數(shù)的介入對展開式計算精度的影響，以及經(jīng)驗公式在復(fù)波數(shù)環(huán)境下的適用性等許多問題仍需進一步研究與討論。

2 復(fù)波數(shù)下多極子展開式的數(shù)值研究

在復(fù)波數(shù)條件下，考慮到格林函數(shù)值若按式(6)計算，衰減項ki的介入會使其隨著多極子點間距離與ki的增大而迅速衰減，而若按式(3)計算，多項式的求和可能會使其衰減效應(yīng)大大減小，從而產(chǎn)生誤差。鑒于此，為了考察復(fù)波數(shù)對格林函數(shù)及其法向?qū)?shù)展開式計算精度的影響，取L點、M點、p點、q點各自之間距離最遠的情況，針對圖3所示單級和多極傳遞關(guān)系，數(shù)值計算各展開式值并與理論值比較，分析討論其計算誤差及產(chǎn)生原因。

以1/8子組劃分法為例，如圖3所示四點三維關(guān)系圖，(a)圖為單級傳遞關(guān)系圖，由于其最低級為2級，其格林函數(shù)及其法向?qū)?shù)的展開式值可直接由(3)式和(9)式計算。(b)圖為多級傳遞關(guān)系圖，由p、q點位置關(guān)系可知級數(shù)I為2，其最低級為4級，因而其格林函數(shù)及其導(dǎo)數(shù)的展開式值由式(10)和式(12)計算。圖中 q1即 λm4，q2即 λm3，M 即 λm2，L 即 λm2'，p2即 λm3'，p1即λm4'。吸聲材料為長纖維玻璃絲綿，實驗測得該吸聲材料的特性阻抗Zb和波數(shù)kb的表達式為:

圖3 四點三維關(guān)系圖Fig.3 Geometry of the four points in three dimensions

其中:za和ka分別為空氣的特性阻抗與波數(shù)，ρa為空氣密度，當材料的填充密度為200 g/L時，測得此材料的流阻率σ為17378 Rayls/m。吸聲材料流阻率的測量方法可參考文獻［15］。

2.1 單級傳遞

以圖3(a)所示四點單級傳遞關(guān)系為例，計算式(2)格林函數(shù)理論值與式(3)展開式值，得到圖4至圖6。觀察圖4和圖5發(fā)現(xiàn)，在空氣介質(zhì)中，無論式(13)中α取何值(0～1)，格林函數(shù)展開式值均能保證較高的計算精度，而對于吸聲材料介質(zhì)，當ki·rLM達到某值后，其展開式值會隨著ki·rLM的增大而與理論值愈發(fā)背離。圖5展示了當α取不同值時，展開式值與理論值的吻合情況，可見在兩種變化方案中(定頻率，變rLM;定rLM，變頻率)，α取0.8時吻合情況最好，失真范圍最小。而且，無論是改變rLM還是改變頻率，展開式值與理論值的分歧點均與ki·rLM值密切相關(guān)，約為13。圖6為在兩種變化方案中，α取不同值時Nc隨ki·rLM值的變化曲線，結(jié)合圖5可知，在ki·rLM小于13時，只有α=0.8的Nc值對復(fù)波數(shù)展開式求解最為合適，過大或過小均會對其計算精度產(chǎn)生影響。

取式(13)中α=0.8，同時變化頻率與rLM值，計算式(9)格林函數(shù)法向?qū)?shù)展開式，并與式(7)所示理論值比較，得到圖7和圖8。觀察圖7可見，在空氣介質(zhì)中，法向?qū)?shù)展開式值與理論值基本吻合，具有較高的計算精度，而對于吸聲材料介質(zhì)，則在ki·rLM約大于13后，開始與理論值相背離。圖8展示了兩種介質(zhì)下計算值實部和虛部的相對誤差，可見對于實波數(shù)，除了在k·rLM很小時相對誤差較大外，在其他范圍均有較高的精度。而對于復(fù)波數(shù)，只有在ki·rLM小于13時的誤差基本在允許范圍。

2.2 多級傳遞

以圖3(b)所示四點多級傳遞關(guān)系為例，取式(13)中α=0.8，同時變化頻率與rLM值，計算式(2)格林函數(shù)理論值與式(10)遠場展開式值，得到圖9和圖10。觀察兩圖發(fā)現(xiàn)，多級傳遞關(guān)系并未增加展開式在復(fù)波數(shù)條件下的失真范圍，均在ki·rLM約大于13時與理論值發(fā)生背離，其與單級傳遞具有相似的誤差精度。由此可斷定，式(3)中E(k)項的計算在實波數(shù)和復(fù)波數(shù)條件下均能保證較高的計算精度，誤差主要集中在TLM(k)項的計算中，分析其誤差原因為:由于TLM(k)疊加項中h(1)l(krLM)的幅值隨著l的增大而逐漸增大，且其相位交替變化，由于Nc值隨著ki·rLM的增大而不斷增高，過大的求和次數(shù)使函數(shù)幅值得以積累，導(dǎo)致最終解向極大值發(fā)散。

總之，實波數(shù)條件下，格林函數(shù)及其法向?qū)?shù)展開式的計算只在頻率或rLM值極低時有較大誤差，其余范圍具有較高的計算精度，而且其對求和項數(shù)Nc的選取適應(yīng)性強。復(fù)波數(shù)條件下，當ki·rLM超過某值時，展開式計算誤差會越來越大，展開式的計算精度對Nc值的選取比較敏感，過大或過小均對求解精度及可信值范圍有較大影響。

所以，對于FMBEM在吸聲材料中的應(yīng)用，可選取α =0.8時的Nc值，當ki·rLM大于13時，不可再用原展開式計算，鑒于此時格林函數(shù)及其法向?qū)?shù)值為10－9或10－8的數(shù)量級，且會隨著ki·rLM的增大而越來越小，為了實現(xiàn)FMBEM的編程計算及滿足工程應(yīng)用的精度要求，可將展開式中TLM(k)的值近似為0。除此方法外，應(yīng)用子結(jié)構(gòu)FMBEM，將大尺寸結(jié)構(gòu)模型劃分為若干小尺寸結(jié)構(gòu)分別計算，即可有效的避免產(chǎn)生過大的ki·rLM值，保證展開式有較高的計算精度。

3 計算實例

圖11所示為某膨脹腔阻性消聲器示意圖，其尺寸為 d=0.1 m，D=0.3 m，l=0.47 m，l1=l2=0.1 m，空氣中聲速c=344 m/s，吸聲材料為長纖維玻璃絲綿，材料填充密度為200 g/L，穿孔管只考慮其支撐吸聲材料的作用。若應(yīng)用FMBEM計算，則其根組邊長b至少為0.67 m，在2級交互組中，約為 0.87 m，可知當ki約大于15時，即頻率大于820 Hz時，展開式的計算將進入錯誤值范圍，選用上節(jié)中提到的辦法，取α=0.8時的 Nc值，當 ki·rLM大于13時，將展開式中TLM(k)的值近似為0。為了實現(xiàn)FMBEM在多介質(zhì)復(fù)合聲場中的計算，應(yīng)用子結(jié)構(gòu)FMBEM將此模型按傳播介質(zhì)分為如圖Ω1與Ω2兩子結(jié)構(gòu)，計算其傳遞損失，并與子結(jié)構(gòu)BEM［16］比較其求解精度。圖12展示了應(yīng)用子結(jié)構(gòu)FMBEM與BEM計算其傳遞損失曲線，可見，當頻率大于820Hz時，二者計算結(jié)果仍吻合良好，證實了該方法的有效性與可行性。

4 結(jié)論

本文通過對格林函數(shù)及其外法向?qū)?shù)展開式的數(shù)值計算與研究，發(fā)現(xiàn)吸聲材料介質(zhì)對FMBEM展開式的計算精度有重要影響，由四點單級傳遞和多極傳遞關(guān)系的計算，揭示了誤差的大小與復(fù)波數(shù)的虛部值和展開點間距離的乘積有重要關(guān)系。通過與空氣介質(zhì)中展開式的計算比較發(fā)現(xiàn)，復(fù)波數(shù)環(huán)境下展開式的計算對Nc值的選取十分敏感，取式(13)中α=0.8時的Nc值能保證較大范圍的求解精度，但當ki·rLM大于13時，展開式值開始與理論值相背離。

鑒于此，對于FMBEM在吸聲材料聲場中的應(yīng)用，本文提出了兩種解決辦法:(1)當ki·rLM大于13時，將展開式中TLM(k)項的貢獻視為0。(2)利用子結(jié)構(gòu)FMBEM，將大尺寸結(jié)構(gòu)劃分為若干小結(jié)構(gòu)分別計算，即可有效的避免產(chǎn)生過大的ki·rLM值，保證展開式有較高的計算精度。最后，通過對某膨脹腔阻性消聲器的傳遞損失計算，證實了本文所述方法與技術(shù)的有效性與可行性。

［1］Rokhlin V.Rapid solution of integral equations of classical potential theory［J］.Journal of Computional Physics，1983，60:187－207.

［2］Nishimura N.Fast multipole accelerated boundary integral equation methods［J］.American Society of Mechanical Engineers，2002，55(4):299 －324.

［3］王雪仁，季振林.快速多極子聲學(xué)邊界元法及其應(yīng)用研究［J］.哈爾濱工程大學(xué)學(xué)報，2007，28(7):752－757.

［4］Zhao J S，Chew W C.Three-dimensional multilevel fast multipole algorithm from static to electrodynamics［J］.Microwave Opt.Technol.Lett.，2000，26:43 －48.

［5］Liu Y J，Nishimur N，Otani Y，et al.A fast boundary element method for the analysis of fiber-reinforced composites based on a rigid-inclusion model［J］.ASME Journal of Applied Mechanics，2005，72(1):115－128.

［6］Liu Y J.A new fast multipole boundary element method for solving 2-D stokes flow problems based on a dual BIE formulation［J］. Engineering Analysis with Boundary Elements，2008，32(2):139 －151.

［7］Fukui T，Kozuka M.Analysis of sound reflection and diffraction in half space by fast multipole boundary element method［R］.Proc of 17th Japan Natl Symp on Boundary Element Methods，2000.

［8］Yasuda Y，Sakamoto S，Sakuma T.Application of the fast multipole boundary element method to sound field analysis using a domain decomposition approach［C］.Proceedings of INTER-NOISE 2006(Honolulu)，2006:624.

［9］Pierce A D.Acoustics:an introduction to its physical principles and applications［M］. McGraw-Hill:New York，1981.

［10］Sakuma T，Yasuda Y.Fast multipole boundary element method for large-scale steady-state sound field analysis.part I:setup and validation［J］.Acta Acustica United with Acustica，2002，88:513－525.

［11］Abramowitz M，Stegun I A.Handbook of mathematical functions［M］.Dover:New York，1965.

［12］Rokhlin V.Diagonal forms of translation operators for the Helmholtz equation in three dimensions［J］.Applied and Computational Harmonic Analysis，1993，1(1):82 －93.

［13］Koc S，Chew W C.Calculation of acoustical scattering from a cluster of scatterers［J］.The Journal of Acoustical Society of America，1998，103:721 －734.

［14］Yasuda Y，Sakuma T.Fast multipole boundary element method for large-scale steady-state sound field analysis.Part II:examination of numerical items［J］.Acta Acustica United with Acustica，2003，89:28－38.

［15］Standard test method for airflow resistance of acoustical materials［R］.American Society for Testing and Materials ASTM C 522－87，1993，Philadelphia PA.

［16］季振林.穿孔管阻性消聲器消聲性能計算及分析［J］.振動工程學(xué)報，2005，18(4):453－457.