耦合邊界元法和等效源法的穩(wěn)健CHIEF法

包英超, 向 宇, 陳 潔, 石梓玉

(1. 廣西科技大學(xué) 廣西汽車零部件與整車技術(shù)重點實驗室,廣西 柳州 545006;2. 廣西科技大學(xué) 機械與汽車工程學(xué)院,廣西 柳州 545006;3. 合肥工業(yè)大學(xué) 機械工程學(xué)院,合肥 230009)

在工程聲學(xué)問題計算中,邊界元法(boundary element method,BEM)[1-2]由于其計算精度高、系數(shù)矩陣主對角占優(yōu)、求解穩(wěn)定好、降維求解、適合無限域/半無限域問題的優(yōu)點,被廣泛應(yīng)用于飛機、列車、船舶、建筑以及環(huán)境的噪聲預(yù)測和控制[3-5]。

然而,常規(guī)邊界元法(conventional boundary element method, CBEM)存在兩個缺陷:(1) 解的非唯一性問題[6-7];(2) 奇異積分問題[8]。目前,有兩種經(jīng)典方法解決Helmholtz邊界積分方程的非唯一性問題。第一種是組合亥姆霍茲積分方程公式(combined helmholtz integral equation formulation,CHIEF)法,先通過在內(nèi)域中選取若干個CHIEF點形成補充方程,再與常規(guī)邊界元方程聯(lián)立形成超定方程組,利用最小二乘法求解該超定方程即可獲得唯一解。該方法原理簡單,易于實現(xiàn)。但是,若CHIEF點位于內(nèi)部本征問題的節(jié)點線上時,該方法失效。另一種是Burton-Miller方法,即將常規(guī)邊界積分方程與其法向?qū)?shù)方程線性耦合,當(dāng)耦合系數(shù)取為復(fù)數(shù)時就可得到全頻段唯一解。然而,法向求導(dǎo)增加了積分核函數(shù)的奇異性,在計算中需處理強奇異積分和超強奇異積分。雖然可以采用區(qū)間分割法[9]、解析與半解析法[10-11]和非線性變量轉(zhuǎn)換法(包括雙曲正弦變換[12]、距離變換[13]和指數(shù)變換[14-15])等方法消除或計算奇異性,但這些方法大多處理過程較繁瑣或者得到的邊界積分方程形式復(fù)雜(尤其對曲面或曲線單元),增加了應(yīng)用實施的難度。

對于Helmholtz外聲場問題的求解,Ochmann提出一種源強模擬技術(shù)[16-17],其基本思想源于疊加原理和Helmholtz外問題解的唯一性定理,即虛擬源強疊加產(chǎn)生的外聲場與振動體自身產(chǎn)生的輻射或散射聲場等效。因源強的作用位置、數(shù)量及類型的不同,又可分為單點-多極子法、多點-單極子法、多點-偶極子法、多點-多極子法等,在應(yīng)用中這類方法通稱等效源法[18-19](equivalent source method, ESM)。由于等效源點(虛擬源強)與真實邊界上的配點不重合,因此ESM完全避免了BEM中的奇異積分計算,且可以很容易采用單極子源和偶極子源線性組合的方法來獲得唯一解。此外,ESM的插值形函數(shù)在實邊界的分布精確滿足控制方程[20],其計算精度高于BEM。但ESM矩陣通常具有較大的條件數(shù),呈病態(tài)性,其數(shù)值求解穩(wěn)定性不如BEM。尤其是當(dāng)?shù)刃г次恢迷O(shè)置不當(dāng)時,嚴(yán)重病態(tài)的系數(shù)矩陣會導(dǎo)致求解過程對輸入誤差極度敏感,很難得到穩(wěn)定的數(shù)值解[21]。

綜上所述,CBEM存在解的非唯一性問題;CHIEF法有內(nèi)點選擇失效的風(fēng)險;Burton-Miller法存在超強奇異積分處理和計算精度下降的問題;ESM存在系數(shù)矩陣病態(tài),求解穩(wěn)定性和魯棒性較差的問題。為此,本文提出一種既能克服解的非唯一性問題,又無需奇異積分處理,且具有高計算精度和穩(wěn)定性的耦合CHIEF法。文中首先簡要闡述了CBEM和ESM的基本理論,然后詳細(xì)地推導(dǎo)了耦合CHIEF法的理論公式并分析了其作用機理,最后利用經(jīng)典的聲輻射和聲散射算例驗證了所提方法的有效性和優(yōu)越性。

1 基本理論

1.1 邊界元法

考慮一個三維空間聲輻射外問題,如圖1所示。其中,S為結(jié)構(gòu)體的實際表面邊界,D-為結(jié)構(gòu)體的內(nèi)部域,D+為結(jié)構(gòu)體的外部域(分析域)。在區(qū)域D+內(nèi),聲場控制方程為非齊次Helmholtz方程

圖1 邊界元法示意圖Fig.1 Schematic diagram of the boundary element method

(?2+k2)p(r)=-f(r) ,r∈D+

(1)

式中:p(r)為r處的聲壓;k=ω/c為波數(shù);ω為角頻率;c為聲速;f(r)為空間分布源。應(yīng)用格林第二恒等式可以將方程(1)表示為常規(guī)邊界積分方程(conventional boundary integral equation,CBIE)為

(2)

式中:r和rS分別為場點和源點位置矢量;pin(r)為場點處的入射聲壓;nS為結(jié)構(gòu)表面S上背離分析域的法向;C(r)為奇異性系數(shù),當(dāng)場點r∈S且S為光滑邊界時,C(r)=0.5;當(dāng)場點r∈D+時,C(r)=1;當(dāng)場點r∈D-時,C(r)=0。G(r,rS)為三維自由空間格林函數(shù)

將場點在實邊界面S上進(jìn)行配點[C(r)=0.5],離散式(2)可得到CBEM方程

ASP=BSQ+Pin

(3)

以常數(shù)單元為例,式(3)中,AS、BS的系數(shù)矩陣元素分別為

(4)

其中,

(5)

當(dāng)場點和源點位于同一個單元內(nèi)時,式(4)中的Green函數(shù)及其法向?qū)?shù)分別具有O(1/R)的弱奇異性和O(1/R2)的強奇異性,因此,CBEM在形成系數(shù)矩陣的過程中需要處理對角元素的奇異積分。

通常利用邊界積分方程(2)求解外部聲學(xué)問題時,對應(yīng)的內(nèi)部聲學(xué)問題的特征頻率處不能保證解的唯一性。Burton等提出一種邊界積分方程及其法向?qū)?shù)方程線性組合的方法,消除解的非唯一性問題。式(2)的法向?qū)?shù)邊界積分方程(normal derivative boundary integral equation,NDBIE)為

(6)

式中,n為邊界S上場點r處的內(nèi)法向,如圖1所示。將式(6)與式(2)進(jìn)行線性組合并離散,可得Burton-Miller法的求解方程

(7)

式中,ABM,BBM的系數(shù)矩陣元素分別為

式中,α是組合系數(shù)[22],當(dāng)α為復(fù)數(shù)時,式(7)可在全頻域內(nèi)克服解的非唯一性問題,但由于邊界積分方程的法向求導(dǎo)引入了超強奇異積分,增加了計算的難度。

1.2 等效源法

ESM是在聲源內(nèi)部布置一系列滿足Helmholtz方程和Sommerfeld輻射條件的虛擬源強來模擬聲源向外輻射的聲場。為了便于計算,可采用簡單點源,并將這些點源布置在與實邊界共形并向內(nèi)回縮距離為d的內(nèi)部虛擬面上,點源的數(shù)量和位置通常與實邊界配點一一對應(yīng),如圖2所示。其中,rE為等效源位置,nE為等效源的單位法向。根據(jù)等效源的思想,其邊界場點r處的聲壓可表示為

圖2 等效源法示意圖Fig.2 Schematic diagram of the equivalent source method

(8)

P=AEW+Pin

(9)

Q=BEW+Qin

(10)

式中,W為等效源源強列向量。AE和BE系數(shù)矩陣元素分別為

由于實邊界上的場點和等效源點不重合,因此,AE和BE中不存在奇異積分,其計算難度遠(yuǎn)低于BEM。此外,在ESM中,實邊界的聲壓是由內(nèi)部各個源點輻射的聲壓疊加形成,其分布精確滿足Helmholtz方程[23];而在BEM中,無論采用常數(shù)單元還是高階單元,其實邊界聲壓分布均為近似假設(shè),不滿足Helmholtz方程。因此,在相同的單元數(shù)和單元類型情況下,ESM的計算精度高于BEM。

2 耦合CHIEF法

2.1 公式推導(dǎo)

本文借鑒CHIEF法的思想,將ESM方程作為補充方程,并與CBEM矩陣方程(3)組合,建立一種耦合算法。

聯(lián)立式(9)和式(10),消去虛擬源強列向量W,經(jīng)整理后可得

(11)

在ESM中,矩陣AE和BE通常是病態(tài)的。為了獲得較高的計算穩(wěn)定性,一般需要避免直接計算單個矩陣的逆。分別利用AE和BE左乘式(11),使兩個矩陣成對出現(xiàn),得到

(12)

(13)

將式(12)和式(13)分別與式(3)聯(lián)立,可得以下兩種組合方式:

組合1

(14)

組合2

(15)

方程(14)和方程(15)即為本文提出的耦合CHIEF法的方程。本文的組合方程雖然在形式和原理上與CHIEF法相似,但ESM方程和CBEM方程是不相關(guān)的兩個獨立方程,其中,CBEM方程對應(yīng)實邊界面,ESM方程對應(yīng)虛擬等效源面,兩組方程有著不同的本征頻率,所以,它們各自的方程不會在同一波數(shù)下同時失效。因此,相對于傳統(tǒng)CHIEF法,本文方法具有更高的可靠性和魯棒性。

2.2 CBEM奇異矩陣的間接計算

Xiang等[24]利用系數(shù)矩陣的等價關(guān)系,提出了一種奇異矩陣的間接計算方法,用于間接求解CBEM的奇異矩陣BS。將CBEM的方程(3)改寫為

(16)

根據(jù)解的唯一性原理,通過對比式(12)和式(16)右側(cè)第一項的系數(shù)矩陣,可進(jìn)一步得到CBEM和ESM系數(shù)矩陣之間的等價耦合關(guān)系

ASAE=BSBE

(17)

矩陣BS可由式(17)間接計算

(18)

使用式(18)計算BS不需要直接計算奇異積分,且間接矩陣BS包含了ESM系數(shù)矩陣AE和BE的信息。因此,在一定程度上還可進(jìn)一步繼承ESM的高精度優(yōu)點。

3 數(shù)值算例

3.1 球面模型的聲輻射和聲散射

通過計算幾個典型球面模型的聲輻射和聲散射問題,驗證本文方法的有效性和優(yōu)越性。在算例中均采用三角形平面常數(shù)單元離散真實邊界,虛擬邊界縮進(jìn)距離d取為真實邊界平均單元大小的1.0倍,可獲得與真實邊界離散單元數(shù)量相同且共形的虛擬球面。球面模型的平面三角形離散如圖3所示(990個單元)。實邊界離散為青色網(wǎng)格,等效源布置在虛擬面上(取三角形單元質(zhì)心)。

(1) 聲輻射

情形1:半徑為a的脈動球距離球心r處的聲壓解析值為

(19)

情形2:半徑為a,徑向振速為vncosφ的振動球?qū)嚯x球心r處的聲壓解析值為

(20)

式中:vn為脈動球表面質(zhì)點法向振速幅值;ρ為介質(zhì)的平均密度;c為介質(zhì)中的聲速;φ為場點矢徑與z軸正方向的夾角。

(2) 聲散射

設(shè)有沿著z軸正向入射的平面波pin=eikz被一個中心位于坐標(biāo)原點、半徑為a的球體阻擋,考慮如下三種球面邊界:

情形1:當(dāng)球體表面為剛性邊界時(即表面總振速?(psc+pin)/?r|r=a=0),球面無量綱散射聲場的解析值為

(21)

情形2:當(dāng)球體表面為柔性邊界時(即表面總聲壓(psc+pin)|r=a=0),球面無量綱散射聲場的法向?qū)?shù)的解析值為

(22)

情形3:當(dāng)球體表面為阻抗邊界時(即表面總聲壓和表面總振速滿足如下關(guān)系: [psc+pin+ξ(?psc+?pin)/?r]|r=a=0),球面無量綱散射聲場的解析值為

(23)

3.2 本文方法的有效性

在下面的算例中,球面聲輻射與聲散射問題均僅計算真實邊界上的表面聲壓或振速,取r=a=1 m,聲阻抗取ξ=10,計算波數(shù)范圍設(shè)置為k=0～10。

對于平面常數(shù)單元,在形成CBEM和Burton-Miller的系數(shù)矩陣時,非對角元素采用25點Gauss-Legendre積分計算,AS的主對角元素為0.5,BS、ABM和BBM的主對角元素采用無奇異的線積分計算[25-26]

(24)

(25)

(26)

式中:R(θ)為單元邊界到配點ri的距離;θ為以配點ri為原點的局部極坐標(biāo)系中的角度變量。但需要注意的是,式(24)、式(25)和式(26)的奇異積分解析算法僅適用于平面單元,曲面單元一般采用奇異相消法,計算更為繁瑣。

3.2.1 解的非唯一性及計算誤差比較

首先對比CBEM、Burton-Miller法和本文方法的計算誤差。其中,本文方法分別考慮了替換BS矩陣前后的組合1和組合2。結(jié)果如圖4所示,聲壓相對誤差ε定義為

圖4 球面模型在不同邊界條件下的表面聲壓和聲壓法向?qū)?shù)的誤差對比Fig.4 Error comparison of surface sound pressure and normal derivative of sound pressure in spherical model under different boundary conditions

(27)

由圖4可知,CBEM在本征頻率處失效,Burton-Miller法和本文所提出的方法在整個計算波數(shù)內(nèi)均能有效克服解的非唯一性。從計算誤差來看,Burton-Miller法的誤差最大,其在非本征頻率處的誤差高于CBEM,說明超強奇異積分核函數(shù)影響了計算精度。而本文方法的誤差明顯低于Burton-Miller法,尤其是替換BS之后,其在整個波數(shù)范圍的誤差均低于CBEM。此外,雖然在輻射問題中,組合1和組合2的計算誤差相近,但在聲散射問題中,組合2的計算精度要顯著高于組合1。因此,在下文中將采用組合2+替換BS矩陣的方式進(jìn)行計算比較。

3.2.2 條件數(shù)比較

CBEM、ESM、Burton-Miller法和本文方法的系數(shù)矩陣條件數(shù),如圖5和圖6所示。由圖5和圖6可知,CBEM在本征頻率處出現(xiàn)奇異性的跳躍,其它三種方法均為平滑曲線。進(jìn)一步表明了ESM、Burton-Miller和本文所提方法均能克服解的非唯一性。

圖5 聲壓系數(shù)矩陣的條件數(shù)對比Fig.5 Comparison of the coefficient matrix condition numbers of sound pressure

圖6 聲壓法向?qū)?shù)系數(shù)矩陣的條件數(shù)對比Fig.6 Comparison of the coefficient matrix condition numbers of Normal derivative of sound pressure

在圖5的聲壓系數(shù)矩陣條件數(shù)中,本文方法的條件數(shù)雖然略大于CBEM和Burton-Miller中的條件數(shù),但數(shù)值在30以內(nèi),但遠(yuǎn)小于ESM。同樣,在圖6的聲壓法向?qū)?shù)系數(shù)矩陣條件數(shù)中,本文方法的條件數(shù)最小,其次是Burton-Miller和CBEM,ESM的條件數(shù)則大得多。由此可見,本文方法的求解穩(wěn)定性與BEM相當(dāng),顯著優(yōu)于ESM。

3.3 計算精度和計算效率對比

為了進(jìn)一步的對比本文方法和Burton-Miller法的計算精度,在不同單元數(shù)下分別采用這兩種方法計算振動球輻射和柔性球散射,其相對誤差如圖7和圖8所示。由圖7和圖8可知,本文方法的計算誤差要小于Burton-Miller方法。且隨單元數(shù)量的增加,本文方法的收斂速度略快于Burton-Miller方法。

圖7 不同單元數(shù)振動球表面聲壓誤差對比Fig.7 Comparison of dimensionless sound pressure errors on the surface of the oscillating sphere with different number of elements

圖8 不同單元數(shù)柔性球表面聲壓誤差對比Fig.8 Comparison of scattered sound pressure normal derivative errors on the flexible spherical surface with different number of elements

再來對比計算效率。以阻抗球散射為例,表1中記錄了Burton-Miller法、BEM、ESM和本文方法的計算時間,包括不同單元數(shù)下生成系數(shù)矩陣的時間,以及求解計算整個問題的總時間。其中,求解計算整個問題的總時間包括了輸入計算參數(shù)、生成系數(shù)矩陣以及求解矩陣方程的時間。

表1 不同方法對應(yīng)的計算時間Tab.1 The calculation time corresponding to the different methods

由表1可以看出,Burton-Miller法的矩陣計算時間略大于BEM的計算時間,而ESM的矩陣計算時間要遠(yuǎn)小于Burton-Miller法的計算時間。其次,對比BS矩陣替換前后的計算時間,可以看出,間接替換計算形成BS所需的時間遠(yuǎn)小于直接計算BS的時間。最后可以發(fā)現(xiàn),本文方法的系數(shù)矩陣的計算時間及問題求解的整體時間要顯著小于Burton-Miller法的計算時間。其主要原因是ESM中的AE和BE矩陣在形成的過程中不需要積分計算和無需處理對角元素的奇異積分,因此,其計算效率遠(yuǎn)高于Burton-Miller法或BEM中需要數(shù)值積分的矩陣。且替換計算形成BS矩陣是由AS,AE,BE的簡單運算得到,較于直接計算BS,無需通過數(shù)值積分逐個計算矩陣所有元素,故所需時間也非常短。綜上所述,本文方法在計算效率方面相比于Burton-Miller方法具有一定優(yōu)勢。

3.4 不同結(jié)構(gòu)的聲輻射和聲散射問題的計算結(jié)果

為了進(jìn)一步考察本文方法適應(yīng)任意形狀結(jié)構(gòu)體的情況,本節(jié)采用CBEM、Burton-Miller法以及組合2+替換BS計算光滑連續(xù)邊界的兩端帶球帽的圓柱模型、環(huán)面模型和“音叉”模型,模型及參數(shù)如圖9～圖11所示。在計算三種模型時,將等效源布置在由真實邊界向內(nèi)共形回縮的虛擬面上,縮進(jìn)距離d為1.0倍的真實邊界平均單元大小。考慮如下五種邊界條件:脈動輻射、振動輻射、剛性散射、柔性散射和阻抗散射,然后計算模型上P點的聲壓及其法向?qū)?shù)。由于兩端帶球帽的圓柱模型、環(huán)面模型和“音叉”模型難以獲得解析解,因此,將CBEM的計算結(jié)果作為參考和比較。結(jié)果如圖12～圖14所示。

圖9 兩端帶球帽的圓柱面的三維離散模型及三維視圖參數(shù)Fig.9 3D discrete model and 3D view parameters of cylinder with ball cap at both ends surface

圖10 環(huán)面的三維離散模型及三維視圖參數(shù)Fig.10 3D discrete model and 3D view parameters of torus surface

圖11 “音叉”面的三維離散及三維視圖參數(shù)Fig.11 3D discrete model and 3D view parameters of tuning fork surface

圖12 在不同邊界條件下兩端帶球帽的圓柱面模型表面P點的無量綱聲壓幅值及無量綱聲壓法向?qū)?shù)幅值Fig.12 Dimensionless sound pressure amplitude and normal derivative of sound pressure amplitude of P point on the surface of cylindrical model with ball cap at both ends under different boundary conditions

圖13 在不同邊界條件下環(huán)面模型表面P點的無量綱聲壓幅值及無量綱聲壓法向?qū)?shù)幅值Fig.13 Dimensionless sound pressure amplitude and normal derivative of sound pressure amplitude of P point on the surface of torus under different boundary conditions

圖14 在不同邊界條件下“音叉”面模型表面P點的無量綱聲壓幅值及無量綱聲壓法向?qū)?shù)幅值Fig.14 Dimensionless sound pressure amplitude and normal derivative of sound pressure amplitude of P point on the surface of tuning fork under different boundary conditions

由圖12～圖14可知,在特征頻率下,本文方法有效克服解的非唯一性。在非特征頻率下,本文所提方法的計算結(jié)果與CBEM方法的計算結(jié)果吻合較好,說明本文方法可以用于不同結(jié)構(gòu)模型的聲輻射及聲散射問題計算。

4 結(jié) 論

針對邊界元法的非唯一性問題和奇異積分問題,本文借鑒CHIEF法思想,將CBEM和ESM方程聯(lián)立形成超定方程組,并利用兩組方程系數(shù)矩陣的耦合等價關(guān)系進(jìn)行替換計算,提出了一種在全波數(shù)均具有唯一解的高精度、高穩(wěn)定性的聲學(xué)耦合CHIEF法。文中對所提方法進(jìn)行了詳細(xì)的推導(dǎo)和闡述,通過典型三維聲輻射和聲散射球源算例進(jìn)行了驗證,得到結(jié)論如下:

(1) 由于采用單極子和偶極子線性組合的ESM方程作為補充方程,所提方法在計算波數(shù)域內(nèi)均可獲得唯一解,有效避免了傳統(tǒng)CHIEF法中的內(nèi)點方程失效問題。

(2) 利用奇異矩陣的間接計算,不僅避免了直接對奇異積分的計算,而且還一定程度上繼承ESM高精度的特性。且所提方法僅需對一個系數(shù)矩陣進(jìn)行積分計算,與Burton-Miller法相比,減少了三個積分核函數(shù)的數(shù)值積分計算,可以大幅減少計算時間,提高計算效率。

(3) 本文方法的計算精度遠(yuǎn)高于Burton-Miller法,且具有更好的收斂性。特別是對于球面模型的聲輻射問題,在相同計算精度下,本文方法極大地減小了計算規(guī)模。

(4) 本文方法的系數(shù)矩陣條件數(shù)遠(yuǎn)低于ESM,其條件數(shù)值較好,具有良好的求解穩(wěn)定性。

此外,文中雖然只考慮了平面常數(shù)單元,但所提出的耦合算法及奇異矩陣間接替換的思想完全可以方便地應(yīng)用到曲面或曲線單元中。