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

崔曉兵，季振林

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

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

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

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

1 復波數條件下的多極子展開式

1.1 邊界積分方程

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

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

其中，?/?n指該函數的外法向導數，G為Helmholtz方程的基本解，形式如下:

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

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

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

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

圖2 四點關系圖Fig.2 Geometry of the four points

其中:

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

格林函數外法向導數為:

由于:

其外法向導數可展開為:

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

式中:

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

1.3 多極子展開式的數值計算

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

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

2 復波數下多極子展開式的數值研究

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

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

圖3 四點三維關系圖Fig.3 Geometry of the four points in three dimensions

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

2.1 單級傳遞

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

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

2.2 多級傳遞

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

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

所以，對于FMBEM在吸聲材料中的應用，可選取α =0.8時的Nc值，當ki·rLM大于13時，不可再用原展開式計算，鑒于此時格林函數及其法向導數值為10－9或10－8的數量級，且會隨著ki·rLM的增大而越來越小，為了實現FMBEM的編程計算及滿足工程應用的精度要求，可將展開式中TLM(k)的值近似為0。除此方法外，應用子結構FMBEM，將大尺寸結構模型劃分為若干小尺寸結構分別計算，即可有效的避免產生過大的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，穿孔管只考慮其支撐吸聲材料的作用。若應用FMBEM計算，則其根組邊長b至少為0.67 m，在2級交互組中，約為 0.87 m，可知當ki約大于15時，即頻率大于820 Hz時，展開式的計算將進入錯誤值范圍，選用上節(jié)中提到的辦法，取α=0.8時的 Nc值，當 ki·rLM大于13時，將展開式中TLM(k)的值近似為0。為了實現FMBEM在多介質復合聲場中的計算，應用子結構FMBEM將此模型按傳播介質分為如圖Ω1與Ω2兩子結構，計算其傳遞損失，并與子結構BEM［16］比較其求解精度。圖12展示了應用子結構FMBEM與BEM計算其傳遞損失曲線，可見，當頻率大于820Hz時，二者計算結果仍吻合良好，證實了該方法的有效性與可行性。

4 結論

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

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

［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］王雪仁，季振林.快速多極子聲學邊界元法及其應用研究［J］.哈爾濱工程大學學報，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］.振動工程學報，2005，18(4):453－457.