矩形單元聲場波函數構造及其應用?

賀佐潦霜 向 宇 石梓玉 陳巖豪 陸 靜

(1 廣西科技大學 廣西汽車零部件與整車技術重點實驗室 柳州 545006)

(2 廣西科技大學機械與汽車工程學院 柳州 545006)

0 引言

在振動結構的輻射聲場計算中，邊界元法(Boundary element method,BEM)是一種被廣泛使用的計算方法[1-2]。通常使用的邊界元法包括直接邊界元法(Direct-BEM)[3]和間接邊界元法(Indirect-BEM)[4]，其在計算聲源外部輻射聲場時，由于需要數值積分的計算，因此計算效率較低，而且當積分面與振動體表面重合時，會產生奇異性問題[5]，需要使用額外的數值技術解決[6-7]。若將積分面設置在振動體內部(用于求解外問題)或者外部(用于求解內問題)的一個虛擬表面上，即可避免奇異積分的處理，這種方法被稱為波疊加法(Wave superposition method,WSM)[8-9]或修正邊界元法(Modified BEM,MBEM)[10-11]。該方法避免了奇異積分的處理，但在求解振動體輻射聲場時仍需對離散單元進行數值積分計算，當離散單元數量較多、求解規模較大時，單元的積分計算會耗費大量時間。

WSM 中的數值積分來自于對面源輻射聲場的計算，若將其用單極子點源代替，則可以得到一種無需積分的高效率計算方法，即等效源法(Equivalent source method,ESM)[12-13]。但ESM 中的單極子點源是對面源輻射聲場的過度簡化，其計算精度和穩定性較大程度依賴于等效源位置和數目的選擇。針對以上問題，相關學者提出了優化等效源配置[14-15]、采用具有指向性的射線波函數來代替單極子點源[16-17]等方法，并取得了一定的研究成果。但在面源簡化為點源的過程中，始終存在較大積分近似誤差[11]，一定程度上影響了計算精度。

針對上述缺陷，本文根據外問題中聲源產生的空間聲場，總可以展開成一系列不同階次的球Hankel 函數與球諧波函數乘積的加權和[18]，構造了一種波函數替代WSM中Green 函數在單元區域的積分，避免了在求解振動體外部輻射聲場中復雜的積分計算。以矩形常數單元為例，推導了替代矩形單元積分一般形式和內推形式的波函數，以及當離散單元為正方形時，可將波函數進一步簡化為圓形域內推波函數。文中對比了幾種波函數與直接積分的計算精度和計算效率，并通過簡支板聲源和立方箱體輻射聲源的數值仿真算例，驗證了圓形域內推波函數與ESM在聲場計算中的效果。

1 波疊加積分方程及其離散形式

根據WSM 理論，振動體外部r處的輻射聲壓可以利用式(1)計算：

其中，S′為振動體內部的一個封閉虛擬曲面，q(rE)為虛擬面S′上rE處的源強，G(r,rE)為自由空間Green函數：

其中，|r-rE|為場點r與虛擬面上rE處之間的距離，k為波數，i=

若將虛擬面離散成N個單元，式(1)可寫為

當離散的單元足夠小時，可以將每一個單元的源強qi(rE)均視為常數qEi，式(1)變為

其中，qEi為i個單元的源強。式(4)就是常用的WSM，該方法雖然計算精度較高，但其在計算振動體外部不同場點r處的聲壓時，需要對所有單元計算Green 函數的數值積分，導致整體的計算效率較低。

一種更為簡化的方法是用振動體內部若干個不同源強的等效源產生的聲場，代替振動體輻射的聲場，即將單元積分(面源)直接簡化為點源G(r,rEi)，得到在聲學計算中廣泛使用的ESM：

其中，rEi表示第i個等效源的位置，q(rEi)為第i個等效源的源強。

ESM 由于形式簡單且無需積分，因此計算效率相較于WSM 得到了很大提高。但該方法對物理模型過度簡化，只有當等效源點位置和數目合適時才可達到較高的計算精度，并且從面源簡化為點源的近似過程始終存在較大的積分誤差，影響了計算精度。

針對上述兩種方法缺陷，為了盡可能減小積分近似誤差，提高振動體外部輻射聲場計算效率，本文采用波函數來代替離散單元的數值積分。

2 矩形常數單元的波函數構造

2.1 矩形常數單元的波函數構造

在工程實際中，對于平板或圓柱等結構的外部輻射聲場求解，可采用WSM 設置規則的虛擬平面或圓柱面替代實際聲場，這種虛擬曲面通?？梢员粍澐譃榫匦螁卧?；在近場聲全息中，如果使用常規的平面矩形或正方形陣列，為了共形，也可將虛擬面劃分為矩形單元。因此本文主要考慮離散單元為矩形時的波函數構造。

在單元局部球坐標系下考慮一個矩形單元的積分，將該單元記為S，如圖1所示。其中，單元的局部坐標系的原點位于單元質心ξ，坐標系的x軸和y軸分別平行于矩形兩邊，將這兩邊的長度分別記為Lx和Ly。在單元區域內對Green 函數的積分可用一個函數K(r,θ,?)代替：

圖1 矩形單元S 及其外域球面示意圖Fig.1 Schematic diagram of rectangular element S and its outer sphere

其中，(r,θ,?)表示場點在坐標系的位置，r為矩形單元中心ξ到場點的距離，θ為對應俯仰角，?為對應方位角，r′為單元內部的場點。

由于式(6)是一個空間中的輻射聲場，因此，可以利用Helmholtz方程在單元局部球坐標系下的解將函數K(r,θ,?)近似為如下的波函數形式：

其中，Cnm為球面波展開系數(kr)為第二類n階球Hankel 函數，k為波數，(θ,?)為n階m次的球諧函數，它決定了聲壓在不同角度的輻射屬性，其表達式為

一旦確定了展開系數Cnm，矩形單元的輻射聲場也就隨之確定了。又由于單元積分和波函數K(r,θ,?)都滿足Helmholtz 方程和Sommerfeld 輻射條件，因此根據微分方程的定解理論，只要兩者在某一邊界上等價，那么它們在整個空間中的聲場分布都是等價的。

為計算展開系數Cnm，在單元的外域任選一個半徑為R的人工邊界球面。則該人工邊界上的聲場分布為

其中，R是人工邊界球面上的點。由球諧函數正交歸一性，將式(9)與式(7)聯立，即可求出展開系數：

式(11)可進一步表示為全場坐標變量形式：

其中，|r-r′i|為場點r與單元中心點之間的距離，zi為第i個單元在局域坐標系中z軸方向的單位向量。

由于式(11)中人工邊界球面是任意選取的，因此本文把該式構造的波函數K(r,θ,?)稱為矩形域一般形式波函數。但從式(10)可以看出，在計算展開系數Cnm時需要求解四重積分，因此，當虛擬面離散的單元不一致且數目較多時，將花費大量時間計算展開系數。為縮短計算展開系數Cnm的時間，本文考慮用矩形單元的遠場輻射聲壓解析解代替單元數值積分。

2.2 矩形域內推波函數的構造

若將人工邊界位于遠場，即人工邊界球面的半徑RF遠大于單元尺寸。此時，矩形單元的格林函數有近似解析表達式：

其中，(kx,ky)=(ksinθRFcos?RF,ksin?RFsinθRF)。將式(13)代入矩形單元積分可得

其中，sinc(x)=sin(x)/x。將式(14)代入式(10)可得

其中，d?RF=sin?RFdθRFd?RF，同理，將式(15)代入式(7)可得

式(16)可進一步表示為全場坐標變量形式：

由于式(16)選用的人工邊界位于遠場，因此本文把該式構造的波函數KF(r,θ,?)稱為矩形域內推波函數。由式(15)可見，當選用遠場人工邊界時，僅需二重積分即可求出展開系數Cnm，而一般形式波函數在求解展開系數Cnm需要計算四重積分。當離散的單元數目較多時，可大幅度縮短計算展開系數的時間，從而提高聲場計算效率。當離散的單元為正方形時，本文利用圓形單元替代正方形單元，進一步簡化波函數。

2.3 圓形域內推波函數的構造

當離散的單元為正方形單元S時，可用一個中心點相同、面積相等的圓形域近似代替，則有

如圖2(a)所示?？梢郧蟮?，邊長為L的正方形單元的等效圓域半徑為a=

圖2 正方形單元S 的等效圓域Fig.2 Equivalent circle regionof square element S

由于聲場關于z軸旋轉對稱，因此僅在xz平面上用Legendre 正交函數逼近實際指向性函數即可。如圖2(b)所示，半徑為R(R>a)的球形人工邊界上的實際指向性函數分布為(場點rp選在xz平面人工邊界上)：

其中，rp=(R,θ,0)為人工邊界xz平面上的場點，r′=(ρ,0,?)為單元內部場點。

由于對稱性，式(18)可進一步寫為

人工邊界xz平面上的聲場可用Legendre 正交多項式函數逼近展開：

其中，Pn(cosθ)為Legendre函數。利用Legendre函數的正交性可得展開系數：

與式(13)類似，當選用一個遠場人工邊界RF時，此時圓形單元的Green函數有近似解析表達式：

其中，J1(kasinθRF)為一階貝塞爾函數。與前述同理，可得到圓形域外部輻射聲場的等效波函數為

式(24)可進一步表示為全場坐標變量形式：

由于式(24)的選用的人工邊界位于遠場，因此本文把該式構造的波函數稱為圓形域內推波函數?？梢钥闯?r,θ)是在KF(r,θ,?)的基礎上做了進一步的簡化，構造了一個與m項無關的波函數。并且，從式(23)可以看出，展開系數Cnm由二重積分簡化為一重積分Cn，進一步縮短計算聲場的時間。

2.4 三種波函數的計算精度與計算效率對比

由于本文構造的波函數與WSM 類似，離散的單元位于真實邊界回縮的虛擬邊界上，其聲壓表達式不僅嚴格滿足Helmholtz 方程，且在真實邊界上關于離散結點的聲壓分布還是嚴格滿足形函數定義的高階形函數[19]，因而相較于BEM 所需的單元數更少，對于結構和聲場較為簡單的情況，有時每個波長僅需2 個單元即可。例如，當劃分的單元為正方形時，單元的邊長L必須小于或等于半波長，即L≤λ/2，也就是kL≤π。因此，本文比較3 種波函數在kL=π/4、π/2、π的情況下與單元直接積分的擬合效果。

為了對比構造的3 種波函數的計算精度和計算效率，算例選用位于xy平面上的正方形單元和圓形單元聲源。正方形單元邊長為L，圓形單元的半徑為a=近場人工球面半徑為L，遠場人工球面半徑為105L，單元中心點與人工球面中心點都位于坐標原點。計算半徑為2L的球面聲場，計算的場點數目為105個。算例中的積分均采用Gauss-Legendre 積分，其中，直接積分計算單元聲場時，每個積分變量使用4個高斯點，計算波函數展開系數時，每個積分變量使用20 個高斯點。值得指出的是，波函數的最高階數n選擇合適與否對計算結果的影響很大。圖3 為在kL=π/4、π/2、π 的情況下，不同階數的矩形域內推波函數對應的重建聲壓相對誤差，其中，相對誤差由式(26) 計算：

圖3 矩形域內推波函數不同階數對應的重建聲壓誤差Fig.3 Reconstructed sound pressure errors corresponding to different orders of pushwave function in rectangular domain

其中，pn為第n階波函數重建聲壓向量，p為直接積分聲壓向量。

由圖3 可以看出，在kL分別為π/4、π/2、π 的情況下，當波函數的階數n=1 時，計算聲壓相對誤差較大，但隨著波函數階數的增加，相對誤差呈下降趨勢；當n=4 時，在不同情況下相對誤差都小于0.5%，達到了足夠的計算精度。

進一步對比3 種波函數在不同情況下的計算精度，表1為?=π/2時，3種波函數在kL分別為π/4、π/2、π 的情況下與單元直接積分的擬合圖像與相對誤差。當kL分別為π/4、π/2 時，3 種波函數在計算聲場聲壓時，與直接積分的相對誤差在0.5%以下，達到了足夠的計算精度。當kL=π 時，相當于一個波長內只有兩個單元，此時單元內的聲場更為復雜。即便如此，矩形域一般形式波函數與內推波函數的計算精度達到了99.5%以上，即便是相對簡化的圓形域內推波函數計算精度也能達到98.9%。

表1 kL 分別為π/4、π/2、π 時，3 種波函數與直接積分計算精度對比Table 1 Comparison of calculation accuracy between three wave functions and direct integration at kL is π/4、π/2、π

表2 記錄了表1 中kL=π 時波函數與直接積分消耗的計算時間?？梢钥闯?，在計算場點聲壓時，3 種波函數的計算時間均遠低于直接積分。但是，矩形域一般形式波函數在求展開系數消耗時間較多，導致總時間還要超過直接積分。因此，這種波函數只有在離散的矩形單元完全相同時，才具有較高的聲場計算效率，否則其效率非常低。但采用遠場內推之后，矩形域內推波函數求展開系數的時間大幅縮短，其計算速度約為直接積分的4～5倍，尤其是圓形域內推波函數計算聲場的速度是直接積分的12～13 倍。因此，即便離散單元的形狀和大小不一致，內推波函數仍然具有較高的計算效率。

表2 當kL 為π 時3 種波函數與直接積分計算聲壓的CPU 耗時Table 2 Three kinds of wave functions and CPU time consumption of direct integration to calculate sound pressure at kL is π

為了進一步對比在不同數量場點下矩形域內推波函數、圓形域內推波函數與直接積分的計算效率，設置場點數量從106個增加至9×106個，對比直接積分和波函數所消耗的時間，結果如圖4 所示。由圖4 可見，直接積分所消耗的計算時間約為矩形域內推波函數的5～6 倍，為圓形域內推波函數的12～13 倍。說明波函數的計算效率高于直接積分，且計算的場點數目越多，波函數的優勢越明顯。

圖4 計算不同數量場點聲壓時直接積分和波函數的CPU 耗時Fig.4 CPU time consumption of direct integration and wave function when calculating sound pressure at different number of field points

3 圓形域內推波函數在聲場計算中的應用

3.1 基于圓形域內推波函數的聲場計算公式

對于傳統的ESM，空間任意r處的聲壓為

空間任意r處質點速度為

其中，ρ0為介質的密度，ω為角頻率。在本算例中，將ESM 中的G(r,rEi)采用3 種波函數中精度相對最低、計算效率最高的圓形域內推波函數替代，考察波函數方法在近場聲全息中的計算效果。

與ESM類似，空間任意r處的聲壓和質點速度分別為

其中，Wi為第i個單元等效波函數的源強。

3.2 圓形域內推波函數在近場聲全息中的應用

在實際工程中，板、殼等連續分布的結構振動聲源較為常見。因此，本文以四邊無限大障板的簡支正方形板為研究對象，對比分析圓形域內推波函數與ESM 求解的重建聲場。將簡支板左下角的頂點作為坐標原點建立坐標系，仿真參數如表3 所示。全息面位于簡支板正上方0.03 m 處，均勻分布20×20 個測點。重建面在簡支板正上方0.01 m 處，尺寸與簡支板大小相同，均勻分布40×40 個重建測點。虛源面尺寸與簡支板大小相同，均勻劃分為20×20 個正方形單元，單元長度為0.025 m，并用面積相同的圓形域代替正方形單元，在每一個單元中心布置一個等效源點，模型示意如圖5 所示。將虛源面分別置于重建面下方0.5倍、1倍、2倍單元長度位置，并在仿真中對全息面測量聲壓添加信噪比為30 dB 的高斯白噪聲，然后對比圓形域內推波函數與ESM重建聲場的相對誤差。

表3 仿真參數Table 3 Simulation parameters

圖5 仿真模型示意圖Fig.5 Schematic diagram of simulation model

圖6 為不同情況下ESM 與圓形域內推波函數重建聲壓和振速的相對誤差，其中，相對誤差的定義為

圖6 200 ～3000 Hz 頻率下ESM 與圓形域內推波函數重建誤差對比Fig.6 Comparison of reconstruction error between equivalent source method and circular domain extrapolation function at 200–3000 Hz frequency

其中，E為重建面的解析聲壓或解析振速向量，E′為算法重建結果向量。

由圖6 的計算結果可以看出，將虛源面置于重建面下方不同距離時，ESM與圓形域內推波函數的聲壓和振速重建誤差均隨著重建頻率的增加而增大。但在整個計算頻段，即便對全息面測量聲壓添加了信噪比為30 dB 的高斯白噪聲，波函數法的聲壓和振速重建誤差也低于ESM，且在虛源面與重建面距離較小時優勢更為明顯。

3.3 圓形域內推波函數在聲輻射計算中的應用

對于無解析解的立方箱體結構，在立方箱體內部放置若干單位源強的單極子源，以這些單極子源在箱體表面單元中心點產生的質點速度作為該箱體的速度邊界條件，這些單極子源在單元中心點處的輻射聲壓則作為標準聲壓，進一步驗證圓形域內推波函數的有效性。

以該立方箱體的中心點為坐標原點，邊長為1 m，每面均勻劃分為100 個四邊形單元，總共600個表面單元。計算時，選取每個表面單元中心點坐標乘以縮放系數進行同形縮放來獲得等效源點和圓形單元波函數中心點所在位置坐標(計算時，縮放系數取為0.7)，圓形單元面積為表面四邊形單元面積乘以縮放系數。立方箱體表面單元中心點與等效源點、圓形單元中心點如圖7所示。

圖7 立方箱體表面單元中心點與等效源點、圓形單元中心點示意圖Fig.7 Schematic diagram of unit center point,equivalent source point and circular unit center point on the surface of cubic box

算例選用3 個單位源強的單極子點源，分別布置在(0.3 m,0,0)、(0,0.3 m,0)、(0,0,0.3 m)處，聲速及介質密度等參數與3.2 節相同。本文對比了200～1200 Hz頻率下ESM與圓形域內推波函數計算立方箱體表面聲壓相對誤差。如圖8 所示，當分析頻率等于或接近該模型對應的特征頻率時，ESM與圓形域內推波函數的計算誤差都較大；除特征頻率外，兩種方法的計算誤差都很低，且隨著頻率的增加而增大。但是，在不同的頻率下，圓形域內推波函數的計算誤差始終低于ESM。

圖8 200 ～1200 Hz 頻率下ESM 與圓形域內推波函數計算聲壓相對誤差Fig.8 Relative error of sound pressure calculated by equivalent source method and push wave function in circular domain at 200–1200 Hz frequency

4 結論

為了避免WSM在求解振動體外部輻射聲場中復雜的積分計算，本文構造了一種波函數替代單元區域關于Green 函數的積分。并以矩形常數單元為例，推導了替代單元積分的一般形式和內推形式波函數，以及當單元為正方形時的簡化內推波函數，并對比了3種波函數與直接積分的計算精度和計算效率。最后，通過簡支板聲源和立方箱體輻射聲源的數值仿真算例，驗證了圓形域內推波函數與ESM在聲場計算中的效果，主要結論如下：

(1) 在計算單個單元外部聲場時，本文構造的3種波函數均與單元直接積分高度擬合，且矩形域一般形式和內推形式波函數的計算效率達到了直接積分的5～6 倍，圓形域內推波函數的計算效率達到了直接積分的12～13倍。

(2) 在虛源面劃分單元數目、回退距離相同的情況下，即便是3 種波函數中精度相對較低的圓形域內推波函數，在無噪聲和加入30 dB 高斯白噪聲時，重建簡支板外部聲壓與振速誤差仍低于ESM。

(3) 對于內部放置若干單極子點源的立方箱體，在求解箱體表面聲壓時，當分析頻率接近模型的特征頻率，ESM 與圓形域內推波函數的計算誤差較大。除特征頻率外，兩種方法的計算精度較高，且在不同頻率下，圓形域內推波函數的計算精度顯著高于ESM。