子結構快速多極子邊界元法聲學迭代計算收斂特性分析

崔曉兵

（中國艦船研究院 北京100192）

引 言

將快速多極子算法應用到傳統邊界元法中，研究形成快速多極子邊界元法（FMBEM），憑借其計算效率高、存儲空間要求小和傳統邊界元法具有同等計算精度等特點，被廣泛應用于大尺度聲場的三維計算中［1］，但對于復雜結構聲場計算將存在不足，例如對于復雜結構的排氣消聲器的內部聲場計算，有以下問題需要解決：

一是對于抗性消聲器，其往往具有進出口內插管、穿孔管、隔板等薄壁結構，在離散求解時出現奇異性邊界，FMBEM無法直接應用；二是對于阻抗復合型消聲器，由于結構中存在吸聲阻尼材料，由于聲場波數的不同，聲學域需要按照不同介質進行劃分，各區域間的傳播與耦合問題將使一般的FMBEM應用受到限制。

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

1 子結構快速多極子邊界元法

圖1 聲學域的劃分

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

不同于傳統邊界元法，FMBEM迭代計算時直接求解矩陣行向量與列向量的積，不構成完整系數矩陣，所以，子結構FMBEM不能通過對子結構矩陣的求逆運算來構建整體矩陣。本文通過疊加各子結構矩陣各行向量積以計算結構整體矩陣的向量積，如式（1）所示：

式中：Hmn與Gmn為離散的子結構單元m中的邊界n所形成的系數矩陣。交界面滿足連續性條件如下：

聯立以上三式，即可求解所有邊界上離散節點處物理量。

2 迭代收斂的影響因素

2.1 矩陣構建方式對迭代收斂的影響

鑒于FMBEM迭代求解時無法構建完整的未知量系數矩陣，所以無法應用高斯直接消元法。在迭代算法中，工程中應用最為廣泛的是Krylov子空間方法，然而，其未知系數矩陣的條件數是影響其收斂速度的重要因素，因此引入恰當的矩陣預條件處理技術對大型矩陣問題求解就尤為需要，否則迭代次數將使FMBEM的計算速度優勢大打折扣。對于子結構FMBEM，其整體系數矩陣為帶狀，喪失了矩陣對角占優的特點，同樣影響著迭代計算的收斂特性，而且，隨著模型的改變以及求解方程數目的變化，求解過程的迭代次數也極不穩定，甚至會出現不能收斂的情況。

本文應用BICGSTAB（l）子空間迭代法和ILUT預條件處理技術［4-5］計算子結構FMBEM構建的整體矩陣方程，下面就以圖1（a）所示矩形側支管道內部聲場計算為例，分析其迭代收斂特性，并找到影響其收斂的幾個因素。

離散節點編號順序為：（1）子結構Ω1的實邊界面；（2）子結構Ω1的交界面；（3）子結構Ω2的實邊界面；（4）子結構Ω2的交界面。取列向量構建次序為：P1、P12、P2、V12。構建矩陣方程如式（4）；取列向量構建次序為：P1、P2、P12、V12。構建矩陣方程如式（5）。

采用構建列向量次序為：P1、P12、V12、P2。離散節點編號順序為：（1）Ω1真實邊界面；（2）Ω1交界面；（3）Ω2交界面；（4）Ω2真實邊界面。構建的矩陣方程如式（6）所示。另取離散節點編號順序為：（1）Ω1真實邊界面；（2）Ω1交界面；（3）Ω2真實邊界面；（4）Ω2交界面。構建形成的矩陣方程如式（7）所示。

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

表1 求解各矩陣方程迭代次數

2.2 交界面離散單元數對迭代收斂的影響

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

下頁圖2為迭代次數隨交界面離散單元數變化柱形圖。由圖2可見，交界面的離散單元個數與其迭代收斂次數息息相關，隨著單元個數的增加，迭代次數逐漸增多。為了提高計算效率，在劃分子結構時，應減少產生大尺寸的交界面；另外，在滿足足夠的計算精度的前提下，應減少交界面的離散單元數。圖3為應用子結構FMBEM與傳統BEM在傳遞損失計算上的比較，結果展現良好的吻合性，也證實子結構FMBEM的準確性。

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

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

2.3 吸聲材料對迭代收斂的影響

在消聲器的工程實際應用中，吸聲材料的運用十分廣泛，而對于含有多種傳播介質的聲場計算，傳統的單一區域FMBEM將無法適用，因而子結構FMBEM的發展使多區域多介質的聲場計算成為可能。假設圖1（b）所示為膨脹腔阻性消聲器，子結構Ω1中為空氣，Ω1中填充長纖維玻璃絲綿，取分析頻率為500 Hz，結構尺寸與3.2節相同，穿孔管只考慮其支撐吸聲材料的作用。考察吸聲材料對迭代計算的影響，并計算其傳遞損失。

取空氣密度ρa=1.156 kg/m3，測量到長纖維玻璃絲綿的特性阻抗和波數表達式為：

式中：za為空氣的特性阻抗，ka為空氣波數。當材料填充密度為200 g/L時，實驗測得此材料的流阻率為 17 378 Rayls/m［6］。

圖4針對消聲器聲場中有無吸聲材料比較其迭代次數。可見，當聲場中填充吸聲材料后，由于其阻尼作用，其迭代次數普遍較低，可見將子結構FMBEM應用于阻性消聲器聲場的計算，更能體現其計算速度上的優勢。圖5較好地展現了子結構FMBEM與子結構BEM兩者在頻域計算結果的一致性，證實子結構FMBEM良好的計算精度。

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

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

3 結 論

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

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