崔曉兵，季振林

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

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

目前FMBEM的應(yīng)用多局限于空氣介質(zhì)的單區(qū)域聲場計(jì)算，其在多區(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)的聲場計(jì)算(如內(nèi)插管消聲器，穿孔管消聲器等)成為可能。然而，對于FMBEM在吸聲材料介質(zhì)或多介質(zhì)復(fù)合的聲場中的應(yīng)用，由于復(fù)數(shù)形式波數(shù)的影響，其多極子展開式的適用性及求解精度還有待考察，復(fù)波數(shù)條件下展開式數(shù)值計(jì)算的參數(shù)選取仍需進(jìn)一步研究與分析。

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

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

1.1 邊界積分方程

圖1為邊界平滑的三維內(nèi)部聲場示意圖，包含剛性邊界、振動(dòng)邊界、吸收邊界三種邊界。假設(shè)穩(wěn)態(tài)聲場聲壓隨時(shí)間變化是簡諧的，取時(shí)間項(xiàng)為exp(－ jωt)，其中 j=，則該聲場滿足三維Helmholtz方程，邊界上P點(diǎn)處聲壓可由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é)點(diǎn)進(jìn)行分組，分級，按照每級中各組集的遠(yuǎn)近關(guān)系將聲場分為近場與遠(yuǎn)場，在近場仍用傳統(tǒng)BEM法求解，遠(yuǎn)場則采用快速多極子方法加速求解［2］，最后將兩部分結(jié)果相加即可得到整個(gè)聲場的系數(shù)矩陣和向量積。此種求解聲場問題的方法即為FMBEM(具體實(shí)現(xiàn)過程可參考文獻(xiàn)［10］)。

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

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

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

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

其中:

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

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

由于:

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

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

式中:

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

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

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

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

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

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

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

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

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

2.1 單級傳遞

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

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

2.2 多級傳遞

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

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

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

3 計(jì)算實(shí)例

圖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計(jì)算，則其根組邊長b至少為0.67 m，在2級交互組中，約為 0.87 m，可知當(dāng)ki約大于15時(shí)，即頻率大于820 Hz時(shí)，展開式的計(jì)算將進(jìn)入錯(cuò)誤值范圍，選用上節(jié)中提到的辦法，取α=0.8時(shí)的 Nc值，當(dāng) ki·rLM大于13時(shí)，將展開式中TLM(k)的值近似為0。為了實(shí)現(xiàn)FMBEM在多介質(zhì)復(fù)合聲場中的計(jì)算，應(yīng)用子結(jié)構(gòu)FMBEM將此模型按傳播介質(zhì)分為如圖Ω1與Ω2兩子結(jié)構(gòu)，計(jì)算其傳遞損失，并與子結(jié)構(gòu)BEM［16］比較其求解精度。圖12展示了應(yīng)用子結(jié)構(gòu)FMBEM與BEM計(jì)算其傳遞損失曲線，可見，當(dāng)頻率大于820Hz時(shí)，二者計(jì)算結(jié)果仍吻合良好，證實(shí)了該方法的有效性與可行性。

4 結(jié)論

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

鑒于此，對于FMBEM在吸聲材料聲場中的應(yīng)用，本文提出了兩種解決辦法:(1)當(dāng)ki·rLM大于13時(shí)，將展開式中TLM(k)項(xiàng)的貢獻(xiàn)視為0。(2)利用子結(jié)構(gòu)FMBEM，將大尺寸結(jié)構(gòu)劃分為若干小結(jié)構(gòu)分別計(jì)算，即可有效的避免產(chǎn)生過大的ki·rLM值，保證展開式有較高的計(jì)算精度。最后，通過對某膨脹腔阻性消聲器的傳遞損失計(jì)算，證實(shí)了本文所述方法與技術(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é)報(bào)，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ì)算及分析［J］.振動(dòng)工程學(xué)報(bào)，2005，18(4):453－457.