崔曉兵

（中國(guó)艦船研究院 北京100192）

引 言

將快速多極子算法應(yīng)用到傳統(tǒng)邊界元法中，研究形成快速多極子邊界元法（FMBEM），憑借其計(jì)算效率高、存儲(chǔ)空間要求小和傳統(tǒng)邊界元法具有同等計(jì)算精度等特點(diǎn)，被廣泛應(yīng)用于大尺度聲場(chǎng)的三維計(jì)算中［1］，但對(duì)于復(fù)雜結(jié)構(gòu)聲場(chǎng)計(jì)算將存在不足，例如對(duì)于復(fù)雜結(jié)構(gòu)的排氣消聲器的內(nèi)部聲場(chǎng)計(jì)算，有以下問(wèn)題需要解決：

一是對(duì)于抗性消聲器，其往往具有進(jìn)出口內(nèi)插管、穿孔管、隔板等薄壁結(jié)構(gòu)，在離散求解時(shí)出現(xiàn)奇異性邊界，F(xiàn)MBEM無(wú)法直接應(yīng)用；二是對(duì)于阻抗復(fù)合型消聲器，由于結(jié)構(gòu)中存在吸聲阻尼材料，由于聲場(chǎng)波數(shù)的不同，聲學(xué)域需要按照不同介質(zhì)進(jìn)行劃分，各區(qū)域間的傳播與耦合問(wèn)題將使一般的FMBEM應(yīng)用受到限制。

為了解決上述問(wèn)題，將子結(jié)構(gòu)技術(shù)［2］應(yīng)用到快速多極子邊界元法中，發(fā)展形成子結(jié)構(gòu)快速多極子邊界元法，應(yīng)用于含有薄壁結(jié)構(gòu)和帶有阻性材料的聲場(chǎng)快速計(jì)算中，由于FMBEM直接計(jì)算各行的矩陣向量積，并不形成完整的系數(shù)矩陣，常用Krylov子空間迭代法來(lái)求解線性方程組。所以，子結(jié)構(gòu)快速多極子邊界元法的迭代收斂速度無(wú)疑是其在復(fù)雜結(jié)構(gòu)聲場(chǎng)中成功應(yīng)用的關(guān)鍵。本文將以消聲器內(nèi)部聲場(chǎng)計(jì)算為例，分析影響子結(jié)構(gòu)FMBEM的迭代次數(shù)的主要因素，對(duì)矩陣方程的建立進(jìn)行優(yōu)化，找到使其快速收斂的方法，通過(guò)與傳統(tǒng)邊界元法比較來(lái)驗(yàn)證其計(jì)算精度。為子結(jié)構(gòu)FMBEM在大尺度、高網(wǎng)格數(shù)聲學(xué)問(wèn)題中的應(yīng)用提供指導(dǎo)意見(jiàn)。

1 子結(jié)構(gòu)快速多極子邊界元法

圖1 聲學(xué)域的劃分

圖1（a）和（b）分別為矩形支管與膨脹腔消聲器示意圖。按照子結(jié)構(gòu)技術(shù)的解決思路，聲學(xué)域可劃分為如圖子結(jié)構(gòu)Ω1與子結(jié)構(gòu)Ω2。p1、v1為子結(jié)構(gòu)Ω1真實(shí)邊界上的聲壓與質(zhì)點(diǎn)振速，p2、v2分別為子結(jié)構(gòu)Ω2真實(shí)邊界上的聲壓與質(zhì)點(diǎn)振速。將兩個(gè)子結(jié)構(gòu)的共用面稱為交界面，子結(jié)構(gòu)Ω1的交界面處聲壓與質(zhì)點(diǎn)振速為p12與v12，對(duì)于子結(jié)構(gòu)Ω2則為p21與v21。按照這樣劃分，所構(gòu)成的兩個(gè)子結(jié)構(gòu)均具有完整的邊界，進(jìn)口、出口以及壁面的邊界條件分別為振速邊界、阻抗邊界和剛性邊界，交界面的邊界條件為連續(xù)性邊界，從而應(yīng)用FMBEM 到每個(gè)子結(jié)構(gòu)中［3］。

不同于傳統(tǒng)邊界元法，F(xiàn)MBEM迭代計(jì)算時(shí)直接求解矩陣行向量與列向量的積，不構(gòu)成完整系數(shù)矩陣，所以，子結(jié)構(gòu)FMBEM不能通過(guò)對(duì)子結(jié)構(gòu)矩陣的求逆運(yùn)算來(lái)構(gòu)建整體矩陣。本文通過(guò)疊加各子結(jié)構(gòu)矩陣各行向量積以計(jì)算結(jié)構(gòu)整體矩陣的向量積，如式（1）所示：

式中：Hmn與Gmn為離散的子結(jié)構(gòu)單元m中的邊界n所形成的系數(shù)矩陣。交界面滿足連續(xù)性條件如下：

聯(lián)立以上三式，即可求解所有邊界上離散節(jié)點(diǎn)處物理量。

2 迭代收斂的影響因素

2.1 矩陣構(gòu)建方式對(duì)迭代收斂的影響

鑒于FMBEM迭代求解時(shí)無(wú)法構(gòu)建完整的未知量系數(shù)矩陣，所以無(wú)法應(yīng)用高斯直接消元法。在迭代算法中，工程中應(yīng)用最為廣泛的是Krylov子空間方法，然而，其未知系數(shù)矩陣的條件數(shù)是影響其收斂速度的重要因素，因此引入恰當(dāng)?shù)木仃囶A(yù)條件處理技術(shù)對(duì)大型矩陣問(wèn)題求解就尤為需要，否則迭代次數(shù)將使FMBEM的計(jì)算速度優(yōu)勢(shì)大打折扣。對(duì)于子結(jié)構(gòu)FMBEM，其整體系數(shù)矩陣為帶狀，喪失了矩陣對(duì)角占優(yōu)的特點(diǎn)，同樣影響著迭代計(jì)算的收斂特性，而且，隨著模型的改變以及求解方程數(shù)目的變化，求解過(guò)程的迭代次數(shù)也極不穩(wěn)定，甚至?xí)霈F(xiàn)不能收斂的情況。

本文應(yīng)用BICGSTAB（l）子空間迭代法和ILUT預(yù)條件處理技術(shù)［4-5］計(jì)算子結(jié)構(gòu)FMBEM構(gòu)建的整體矩陣方程，下面就以圖1（a）所示矩形側(cè)支管道內(nèi)部聲場(chǎng)計(jì)算為例，分析其迭代收斂特性，并找到影響其收斂的幾個(gè)因素。

離散節(jié)點(diǎn)編號(hào)順序?yàn)椋海?）子結(jié)構(gòu)Ω1的實(shí)邊界面；（2）子結(jié)構(gòu)Ω1的交界面；（3）子結(jié)構(gòu)Ω2的實(shí)邊界面；（4）子結(jié)構(gòu)Ω2的交界面。取列向量構(gòu)建次序?yàn)椋篜1、P12、P2、V12。構(gòu)建矩陣方程如式（4）；取列向量構(gòu)建次序?yàn)椋篜1、P2、P12、V12。構(gòu)建矩陣方程如式（5）。

采用構(gòu)建列向量次序?yàn)椋篜1、P12、V12、P2。離散節(jié)點(diǎn)編號(hào)順序?yàn)椋海?）Ω1真實(shí)邊界面；（2）Ω1交界面；（3）Ω2交界面；（4）Ω2真實(shí)邊界面。構(gòu)建的矩陣方程如式（6）所示。另取離散節(jié)點(diǎn)編號(hào)順序?yàn)椋海?）Ω1真實(shí)邊界面；（2）Ω1交界面；（3）Ω2真實(shí)邊界面；（4）Ω2交界面。構(gòu)建形成的矩陣方程如式（7）所示。

選取計(jì)算頻率為500 Hz，以空氣為傳播介質(zhì)，模型結(jié)構(gòu)尺寸：l1= 0.35 m，l2= 0.08 m，c= 0.15 m，a1=a2=b=0.05 m，離散Ω1與Ω2得到總節(jié)點(diǎn)數(shù)Na為1 812。表1為求解各矩陣方程迭代次數(shù)的比較，可見(jiàn)，若要改善收斂速度，應(yīng)使盡量多的子塊矩陣的對(duì)角線元素處在總體矩陣的對(duì)角線上。此外發(fā)現(xiàn)，預(yù)處理技術(shù)中的數(shù)值丟棄閥值和填充的非零元素個(gè)數(shù)也一定程度影響著迭代收斂次數(shù)。

表1 求解各矩陣方程迭代次數(shù)

2.2 交界面離散單元數(shù)對(duì)迭代收斂的影響

以圖1（b）所示膨脹腔消聲器內(nèi)部聲場(chǎng)計(jì)算為例，分析頻率為500 Hz，傳播介質(zhì)為空氣，結(jié)構(gòu)尺寸為d= 0.049 m、D= 0.164 4 m、l= 0.257 2 m，按照上述要求構(gòu)建矩陣方程，選取ILUT（10-5，50）預(yù)處理方案，根據(jù)不同網(wǎng)格尺寸離散其邊界，應(yīng)用子結(jié)構(gòu)FMBEM分析迭代次數(shù)，并計(jì)算傳遞損失。

下頁(yè)圖2為迭代次數(shù)隨交界面離散單元數(shù)變化柱形圖。由圖2可見(jiàn)，交界面的離散單元個(gè)數(shù)與其迭代收斂次數(shù)息息相關(guān)，隨著單元個(gè)數(shù)的增加，迭代次數(shù)逐漸增多。為了提高計(jì)算效率，在劃分子結(jié)構(gòu)時(shí)，應(yīng)減少產(chǎn)生大尺寸的交界面；另外，在滿足足夠的計(jì)算精度的前提下，應(yīng)減少交界面的離散單元數(shù)。圖3為應(yīng)用子結(jié)構(gòu)FMBEM與傳統(tǒng)BEM在傳遞損失計(jì)算上的比較，結(jié)果展現(xiàn)良好的吻合性，也證實(shí)子結(jié)構(gòu)FMBEM的準(zhǔn)確性。

圖2 迭代次數(shù)隨交界面離散單元數(shù)變化柱形圖

圖3 膨脹腔消聲器傳遞損失曲線

2.3 吸聲材料對(duì)迭代收斂的影響

在消聲器的工程實(shí)際應(yīng)用中，吸聲材料的運(yùn)用十分廣泛，而對(duì)于含有多種傳播介質(zhì)的聲場(chǎng)計(jì)算，傳統(tǒng)的單一區(qū)域FMBEM將無(wú)法適用，因而子結(jié)構(gòu)FMBEM的發(fā)展使多區(qū)域多介質(zhì)的聲場(chǎng)計(jì)算成為可能。假設(shè)圖1（b）所示為膨脹腔阻性消聲器，子結(jié)構(gòu)Ω1中為空氣，Ω1中填充長(zhǎng)纖維玻璃絲綿，取分析頻率為500 Hz，結(jié)構(gòu)尺寸與3.2節(jié)相同，穿孔管只考慮其支撐吸聲材料的作用?？疾煳暡牧蠈?duì)迭代計(jì)算的影響，并計(jì)算其傳遞損失。

取空氣密度ρa(bǔ)=1.156 kg/m3，測(cè)量到長(zhǎng)纖維玻璃絲綿的特性阻抗和波數(shù)表達(dá)式為：

式中：za為空氣的特性阻抗，ka為空氣波數(shù)。當(dāng)材料填充密度為200 g/L時(shí)，實(shí)驗(yàn)測(cè)得此材料的流阻率為 17 378 Rayls/m［6］。

圖4針對(duì)消聲器聲場(chǎng)中有無(wú)吸聲材料比較其迭代次數(shù)。可見(jiàn)，當(dāng)聲場(chǎng)中填充吸聲材料后，由于其阻尼作用，其迭代次數(shù)普遍較低，可見(jiàn)將子結(jié)構(gòu)FMBEM應(yīng)用于阻性消聲器聲場(chǎng)的計(jì)算，更能體現(xiàn)其計(jì)算速度上的優(yōu)勢(shì)。圖5較好地展現(xiàn)了子結(jié)構(gòu)FMBEM與子結(jié)構(gòu)BEM兩者在頻域計(jì)算結(jié)果的一致性，證實(shí)子結(jié)構(gòu)FMBEM良好的計(jì)算精度。

圖4 迭代次數(shù)隨交界面離散單元數(shù)變化柱形圖

圖5 膨脹腔阻性消聲器傳遞損失曲線

3 結(jié) 論

本文通過(guò)應(yīng)用子結(jié)構(gòu)FMBEM對(duì)各種結(jié)構(gòu)消聲器內(nèi)部聲場(chǎng)的計(jì)算，在選取恰當(dāng)?shù)念A(yù)處理方法下，矩陣方程未知量矩陣向量的構(gòu)建次序，交界面的選取與離散以及吸聲材料的填充均對(duì)迭代次數(shù)具有一定影響。其中，模型求解過(guò)程構(gòu)建矩陣方程時(shí)要注意：一是安排對(duì)角占優(yōu)子矩陣盡量位于總體矩陣的長(zhǎng)對(duì)角線上；二是采取恰當(dāng)?shù)念A(yù)處理方法有效地減少迭代次數(shù)；三是過(guò)大的交界面會(huì)降低迭代收斂速度；四是當(dāng)求解結(jié)構(gòu)中帶有吸聲阻尼材料時(shí)，子結(jié)構(gòu)FMBEM的迭代次數(shù)將會(huì)明顯降低，且隨著頻率增大，其迭代特性更為穩(wěn)定。因而，子結(jié)構(gòu)FMBEM更適用于求解阻性和阻抗復(fù)合型消聲器三維聲場(chǎng)計(jì)算。此外，通過(guò)子結(jié)構(gòu)FMBEM與BEM在消聲器傳遞損失計(jì)算的比較，驗(yàn)證了該方法的正確性與計(jì)算精度。

［1］ Nishimura N.Fast multipole accelerated boundary integral equation methods［J］.Appl Mech Rev，2002（4）：299-324.

［2］ Ji Zhen-lin.Acoustic attenuation performance prediction and analysis of perforated tube dissipative silencers［J］.Journal of Vibration Engineering，2005（4）：453-457.

［3］ 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.

［4］ Saad Y ILUT：A dual threshold incomplete LU factorization［J］.Numer.Linear Algebra Appl.，1994（4）：387-402.

［5］ Saad Y.Iterative Methods for Sparse Linear Systems ［M］.Boston：PWS Publishing Company，1996.

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