波疊加法解的非唯一性問題改進及聲源位置優化研究

夏雪寶，向 陽，張 波，程 鵬

（1.廣州廣電計量檢測股份有限公司，廣州 510000；2.武漢理工大學能源與動力工程學院，武漢 430063）

0 引 言

目前，邊界元（BEM）的方法廣泛用于任意形狀結構的聲輻射問題求解，但在求解計算時需對奇異積分進行數值處理，且當遇到會產生非唯一性問題的積分時，對奇異積分的數值處理會更加困難[1-3]。為避免聲學求解時遇到的奇異積分問題，Koopmann 等[4]提出一種可避免奇異積分的聲場求解方法——波疊加法，該方法是將若干虛擬簡單源（單極子源）置于輻射體內部光滑曲面上，通過給定的速度邊界條件求得聲源的源強后，線性疊加虛擬簡單源的輻射聲場以獲得輻射體的聲場解。

采用波疊加法求解時，單極子源均勻置于輻射體內部封閉光滑曲面上，其位置對計算的精度有很大影響[5]。由于單極子源位置為光滑封閉的曲面，當分析波數等于或接近內部Dirchlet 條件對應的特征波數時，就會出現解的非唯一性問題。針對波疊加法解的非唯一性問題，文獻[6]采用三極子源（單極子與偶極子組合）可有效解決非唯一性問題，但由于簡單源變為三極子，使得計算效率降低。文獻[7]將單極子源所在光滑曲面到輻射體中心點的矢徑模值取為復數，也能保證聲場解的唯一解，但由于矢徑模值的取值憑經驗選取，計算精度較難保證。

CHIEF 法（Combined Helmholtz Integral Equation Formulation）是一種改進的聲輻射邊界元法[8]。該方法通過在輻射體內部增加若干CHIEF點，聯立這些點處與邊界表面節點的Helmholtz積分方程[9]，可有效保證聲場的唯一解。文獻[10]提出一種用于求解聲場問題的基本解方法MFS（Method of Funda?mental Solutions），該方法通過疊加一組輻射體內部隨機布置的聲源聲場以獲得輻射體結構總的聲場解，數值算例表明該方法亦可保證全波數域內聲場的唯一解，但需對聲源位置進行優化以保證聲場解的計算精度[11]。

本文基于CHIEF 法和基本解法（MFS）的思想及前期研究成果[12]，通過在輻射體內部添加額外的單極子源（稱為附加源）來保證聲場的唯一解。波疊加法計算結果的總體精度與輻射體表面速度重構密切相關。Koopmann 對比了配置點法、最小二乘法及體積速度匹配法，指出采用體積速度匹配對重構表面速度邊界條件及計算輻射聲場具有最優的精度[13]。因此本文采用附加源波疊加法求解輻射聲場，以重構后體積速度相對誤差為目標函數，獲取單極子源的最優位置，算例表明，該方法的聲場計算精度滿足需求。

1 波疊加法及解的非唯一性

如圖1所示，波疊加法計算時，通過疊加N個輻射體內部光滑曲面上聲源的輻射聲場，可得到表面S外部空間場點的聲壓為[13-15]

圖1 波疊加法聲源位置示意圖Fig.1 Monopole source location of wave superposition method

式中，sn為等效聲源強度；G(rm|rn)= eikR/R為自由空間格林函數[12]；rn和rm分別為聲源點及空間場點位置向量；R為聲源點rn與空間場點rm之間的距離；i為虛數單位；k為波數；?n為對聲源點求梯度；nn為偶極子聲源單位法向量。α= 1，β= 0時表示聲源為單極子源；α= 0，β= i/k表示聲源為偶極子源。

理想流體介質（空氣）質點速度與聲壓關系為v(rm)=?p(rm)/ikρc，因此由式（1）可得

式中，ρc為空氣特性阻抗；nm為邊界表面位置單位法向量；?m為對邊界表面質點求梯度。對式（2）兩邊進行積分，可得輻射體邊界表面單元的體積速度為

式中，Sm為表面單元m的面積區域。式（3）可用矩陣的形式表示為

式中，u、s分別為體積速度向量及聲源強度向量。矩陣U中的元素為

由式（1）可知，α、β的取值決定了波疊加法聲源類型。文獻[16]指出：當等效聲源為單極子源或偶極子源時，波疊加積分方程的分析波數為內部Dirichlet條件或Neumann 條件所對應的特征波數，聲場均會出現解的非唯一現象；而采用三極子波疊加法時，特征波數不會等于或接近內部Robin 條件對應的特征波數，該方法可以保證全波數域聲場解的唯一性。

2 附加源波疊加法

波疊加法相比基本解法（MFS）的不同之處在于前者將聲源置于輻射體內部光滑曲面上，而非在輻射體內部任意放置。CHIEF 主要思想是通過在輻射體內部添加額外點與邊界表面聯立積分方程來保證聲場解的唯一性。因此，基于CHIEF 法和基本解法（MFS）法的思想，通過在輻射體內部添加若干個額外的單極子源來保證分析波數域內聲場的唯一解，添加的額外單極子源稱為附加源，如圖2 所示，相應的聲場解為

圖2 單極子源及附加源位置示意圖Fig.2 Monopole source and additional source location of wave superposition method

式中，A為附加源的個數。為發揮附加源的作用，附加源的位置要避免與單極子源所在光滑曲面重合。聲場整體求解精度仍取決于輻射體內部光滑曲面的位置及單極子源數量N，因此需對單極子源的數量及位置進行優化選取，以確保聲場計算的精度及效率最優。

3 聲源位置優化計算

波疊加法計算時，為保證速度邊界條件能精確地被虛擬單極子源等效，可通過給定的速度邊界條件求得聲源的源強后，對輻射體表面的體積速度進行重構，對比重構前后體積速度的相對誤差。因此，本文以重構前后體積速度的相對誤差為目標函數，對結構輻射體內部單極子源所在曲面位置進行優化。表面離散為N個聲學網格輻射體的表面體積速度相對誤差定義為

式中，u1,i為結構輻射體表面單元i重構后的體積速度，u0,i為給定或已知的表面單元i的體積速度。

實際采用波疊加法求解過程中，一定范圍內增加單極子源個數可提高聲場求解計算的精度，但單極子源個數過多時，單極子之間距離過近會導致計算精度降低，同時聲源數量過多會使得計算效率降低。通過建立邊界網格模型，確保表面邊界單元尺寸小于λ6（λ為聲波波長）以確定簡單聲源數量。而輻射體內部聲源所在曲面位置則可通過對結構表面單元網格節點縮小得到。

本文采用附加源波疊加法計算聲場時，以分析波數范圍內重構前后體積速度相對誤差為目標函數，以同形縮小系數為優化變量對聲源所在曲面位置進行優化。優化的數學模型為

式中，k為分析波數；cs為同形縮小系數。縮小系數cs的取值范圍不能過大或過小，縮小系數cs過大或過小時，單極子源分別距結構表面節點或單極子源相互之間距離過近，導致積分計算時計算精度低。由于cs范圍小，實際優化計算時，可以采用遍歷法得到最優cs。因此，本文采用改進的源波疊加法優化計算流程為：第一步，將輻射體表面劃分為聲學網格，確定單極子源個數；第二步，以重構前后體積速度相對誤差為目標函數，優化得到最佳縮小系數cs，確定最優聲源位置；第三步，添加適量附加單極子源來保證分析波數范圍內聲場的唯一解。

4 數值算例

4.1 脈動球源

4.1.1 單極子波疊加法球心位于坐標原點的脈動球源聲場解析解為[3]

式中，r為坐標原點與空間場點間的距離，v為表面振速幅值，a為半徑。現假設a=1 m，v=1 m/s，在球源內部半徑為0.5 m（b=0.5a）的球面上均勻添加59 個單極子源。采用單極子波疊加法及解析解分別求解坐標（0 m，0 m，1 m）處聲壓，聲壓實部與虛部計算結果如圖3 所示，分析波數范圍為0～20，步長為0.01。

圖3 采用單極子波疊加法及解析解求得不同坐標處聲壓對比Fig.3 Comparison of sound pressures at different coordinates by monopole sources wave superposition method and analytical solution

由圖3 可知，單極子聲源所在曲面為b=0.5a內部球面時，特征波數處（kb= π,2π,3π…）聲場會出現解非唯一性問題，其它波數處聲場解精度較高。

4.1.2 附加源波疊加法

采用改進的波疊加法（附加源波疊加法）進行計算時，在4.1.1節基礎上，于坐標（0 m，0 m，0 m）處即坐標原點添加一個單極子源，計算結果如圖4所示。

圖4 采用附加源波疊加法及解析解求得不同坐標處聲壓對比Fig.4 Comparison of sound pressures at different coordinates by additional sources wave superposition method and analytical solution

由圖4可知，脈動球源中心處添加一個單極子源求解的表面聲壓與解析解完全一致，有效克服了聲場解的非唯一性問題。

4.1.3 計算耗時及計算精度對比

對比附加源與三極子波疊加法的計算耗時和精度，以上述脈動球源為例，分別求解脈動球源表面聲壓相對誤差以及計算所需時間，其中表面聲壓相對誤差定義為

式中，p1,i為聲壓解析解，p2,i為聲壓數值解，N為輻射體表面節點個數。計算結果如圖5所示。

圖5 表面聲壓相對誤差及計算耗時Fig.5 Relative error of surface pressure and computing time

由求解的表面聲壓相對誤差結果可知，附加源波疊加法計算誤差幾乎為0，相應的計算時間僅為三極子波疊加法的50%左右。

4.2 復雜幾何體

如圖6 所示，以總包絡尺寸1.5 m × 1m × 1m 復雜幾何結構體為對象，建立該結構的邊界元模型（表面共有600 個單元）。波疊加法計算時，該模型的單極子源位置由該結構表面單元中心位置坐標縮小后得到。表面速度邊界條件及表面標準聲壓則由位于該幾何體中心坐標原點處的簡單源產生。

圖6 復雜幾何結構體邊界元網格模型Fig.6 Boundary element model of complex volume

采用單極子波疊加法求解不同縮小系數cs時體積速度相對誤差，分析不同cs時對應聲場解的非唯一性問題。cs的范圍為0.4～0.85，步長為0.01，分析的波數為0～20，步長為0.01。計算的結果如圖7所示。

由圖7 可以得出，不同的cs對應的特征波數點不一致，且無規律可循，圖中的尖峰或亮點為對應的特征波數。從圖中顏色分布可以看出，表面體積速度相對誤差隨cs呈兩頭大、中間小的趨勢。由于特征波數處存在解的非唯一性問題，因此不能采用分析波數范圍內體積速度相對誤差的總和或均值對計算精度進行評判，進而確定最優的cs值。

圖7 體積速度相對誤差Fig.7 Relative error of volume velocity

本文通過平均值法計算體積速度相對誤差俯視圖的灰度矩陣，再統計分析各cs值下灰度總值趨勢，進而確定最優cs值，結果如圖8 所示。由圖8 可得出，cs= 0.62 時，體積速度相對誤差灰度總值最小，即為計算精度的最優解。灰度總值趨勢與體積速度相對誤差分布一致，呈現兩頭大中間小的趨勢。

圖8 體積速度相對誤差灰度趨勢曲線Fig.8 Gray scale distribution curve of the volume velocity relative error

確定最佳cs值以后，采用附加源波疊加法對該模型進行計算，以克服聲場在k=6.21處解的非唯一性問題。確定縮小系數cs= 0.62 后，添加了5 個附加源的計算結果如圖9 所示，復雜幾何體內部添加的5個單極子源位置由軟件隨機確定，且經對比計算發現數量為5時效果最好。由圖9可知，cs= 0.62時，5 個附加源可有效解決復雜幾何體聲場特征波數（k=6.21）處解的非唯一性問題；同時整體計算誤差lg(δp%) < 0 即δp% < 1，計算精度非常高，說明以體積速度相對誤差為目標函數優化得到最佳cs值滿足計算精度要求。

圖9 表面聲壓相對誤差Fig.9 Relative error of surface pressure

5 結 語

對應復雜幾何體，波疊加法計算時以表面體積速度相對誤差為目標函數，優化得到最佳單極子源位置，再通過輻射體內部添加適量額外的單極子源可有效保證聲場的唯一解，且滿足計算精度需求。但對如何優化確定附加源的數量及位置，以保證全波數域聲場的唯一解有待后續進行深入的研究。