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

（1.廣州廣電計(jì)量檢測股份有限公司，廣州 510000；2.武漢理工大學(xué)能源與動(dòng)力工程學(xué)院，武漢 430063）

0 引 言

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

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

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

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

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

如圖1所示，波疊加法計(jì)算時(shí)，通過疊加N個(gè)輻射體內(nèi)部光滑曲面上聲源的輻射聲場，可得到表面S外部空間場點(diǎn)的聲壓為[13-15]

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

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

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

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

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

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

由式（1）可知，α、β的取值決定了波疊加法聲源類型。文獻(xiàn)[16]指出：當(dāng)?shù)刃曉礊閱螛O子源或偶極子源時(shí)，波疊加積分方程的分析波數(shù)為內(nèi)部Dirichlet條件或Neumann 條件所對應(yīng)的特征波數(shù)，聲場均會(huì)出現(xiàn)解的非唯一現(xiàn)象；而采用三極子波疊加法時(shí)，特征波數(shù)不會(huì)等于或接近內(nèi)部Robin 條件對應(yīng)的特征波數(shù)，該方法可以保證全波數(shù)域聲場解的唯一性。

2 附加源波疊加法

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

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

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

3 聲源位置優(yōu)化計(jì)算

波疊加法計(jì)算時(shí)，為保證速度邊界條件能精確地被虛擬單極子源等效，可通過給定的速度邊界條件求得聲源的源強(qiáng)后，對輻射體表面的體積速度進(jìn)行重構(gòu)，對比重構(gòu)前后體積速度的相對誤差。因此，本文以重構(gòu)前后體積速度的相對誤差為目標(biāo)函數(shù)，對結(jié)構(gòu)輻射體內(nèi)部單極子源所在曲面位置進(jìn)行優(yōu)化。表面離散為N個(gè)聲學(xué)網(wǎng)格輻射體的表面體積速度相對誤差定義為

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

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

本文采用附加源波疊加法計(jì)算聲場時(shí)，以分析波數(shù)范圍內(nèi)重構(gòu)前后體積速度相對誤差為目標(biāo)函數(shù)，以同形縮小系數(shù)為優(yōu)化變量對聲源所在曲面位置進(jìn)行優(yōu)化。優(yōu)化的數(shù)學(xué)模型為

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

4 數(shù)值算例

4.1 脈動(dòng)球源

4.1.1 單極子波疊加法球心位于坐標(biāo)原點(diǎn)的脈動(dòng)球源聲場解析解為[3]

式中，r為坐標(biāo)原點(diǎn)與空間場點(diǎn)間的距離，v為表面振速幅值，a為半徑?，F(xiàn)假設(shè)a=1 m，v=1 m/s，在球源內(nèi)部半徑為0.5 m（b=0.5a）的球面上均勻添加59 個(gè)單極子源。采用單極子波疊加法及解析解分別求解坐標(biāo)（0 m，0 m，1 m）處聲壓，聲壓實(shí)部與虛部計(jì)算結(jié)果如圖3 所示，分析波數(shù)范圍為0～20，步長為0.01。

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

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

4.1.2 附加源波疊加法

采用改進(jìn)的波疊加法（附加源波疊加法）進(jìn)行計(jì)算時(shí)，在4.1.1節(jié)基礎(chǔ)上，于坐標(biāo)（0 m，0 m，0 m）處即坐標(biāo)原點(diǎn)添加一個(gè)單極子源，計(jì)算結(jié)果如圖4所示。

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

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

4.1.3 計(jì)算耗時(shí)及計(jì)算精度對比

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

式中，p1,i為聲壓解析解，p2,i為聲壓數(shù)值解，N為輻射體表面節(jié)點(diǎn)個(gè)數(shù)。計(jì)算結(jié)果如圖5所示。

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

由求解的表面聲壓相對誤差結(jié)果可知，附加源波疊加法計(jì)算誤差幾乎為0，相應(yīng)的計(jì)算時(shí)間僅為三極子波疊加法的50%左右。

4.2 復(fù)雜幾何體

如圖6 所示，以總包絡(luò)尺寸1.5 m × 1m × 1m 復(fù)雜幾何結(jié)構(gòu)體為對象，建立該結(jié)構(gòu)的邊界元模型（表面共有600 個(gè)單元）。波疊加法計(jì)算時(shí)，該模型的單極子源位置由該結(jié)構(gòu)表面單元中心位置坐標(biāo)縮小后得到。表面速度邊界條件及表面標(biāo)準(zhǔn)聲壓則由位于該幾何體中心坐標(biāo)原點(diǎn)處的簡單源產(chǎn)生。

圖6 復(fù)雜幾何結(jié)構(gòu)體邊界元網(wǎng)格模型Fig.6 Boundary element model of complex volume

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

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

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

本文通過平均值法計(jì)算體積速度相對誤差俯視圖的灰度矩陣，再統(tǒng)計(jì)分析各cs值下灰度總值趨勢，進(jìn)而確定最優(yōu)cs值，結(jié)果如圖8 所示。由圖8 可得出，cs= 0.62 時(shí)，體積速度相對誤差灰度總值最小，即為計(jì)算精度的最優(yōu)解?；叶瓤傊第厔菖c體積速度相對誤差分布一致，呈現(xiàn)兩頭大中間小的趨勢。

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

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

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

5 結(jié) 語

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