波疊加法解的非唯一性問(wèn)題改進(jìn)及聲源位置優(yōu)化研究

夏雪寶，向 陽(yáng)，張 波，程 鵬

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

0 引 言

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

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

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

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

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

如圖1所示，波疊加法計(jì)算時(shí)，通過(guò)疊加N個(gè)輻射體內(nèi)部光滑曲面上聲源的輻射聲場(chǎng)，可得到表面S外部空間場(chǎng)點(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)及空間場(chǎng)點(diǎn)位置向量；R為聲源點(diǎn)rn與空間場(chǎng)點(diǎn)rm之間的距離；i為虛數(shù)單位；k為波數(shù)；?n為對(duì)聲源點(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為對(duì)邊界表面質(zhì)點(diǎn)求梯度。對(duì)式（2）兩邊進(jìn)行積分，可得輻射體邊界表面單元的體積速度為

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

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

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

2 附加源波疊加法

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

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

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

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

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

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

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

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

式中，k為分析波數(shù)；cs為同形縮小系數(shù)?？s小系數(shù)cs的取值范圍不能過(guò)大或過(guò)小，縮小系數(shù)cs過(guò)大或過(guò)小時(shí)，單極子源分別距結(jié)構(gòu)表面節(jié)點(diǎn)或單極子源相互之間距離過(guò)近，導(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)前后體積速度相對(duì)誤差為目標(biāo)函數(shù)，優(yōu)化得到最佳縮小系數(shù)cs，確定最優(yōu)聲源位置；第三步，添加適量附加單極子源來(lái)保證分析波數(shù)范圍內(nèi)聲場(chǎng)的唯一解。

4 數(shù)值算例

4.1 脈動(dòng)球源

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

式中，r為坐標(biāo)原點(diǎn)與空間場(chǎng)點(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，步長(zhǎng)為0.01。

圖3 采用單極子波疊加法及解析解求得不同坐標(biāo)處聲壓對(duì)比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π…）聲場(chǎng)會(huì)出現(xiàn)解非唯一性問(wèn)題，其它波數(shù)處聲場(chǎng)解精度較高。

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)處聲壓對(duì)比Fig.4 Comparison of sound pressures at different coordinates by additional sources wave superposition method and analytical solution

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

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

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

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

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

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

4.2 復(fù)雜幾何體

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

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

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

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

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

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

圖8 體積速度相對(duì)誤差灰度趨勢(shì)曲線(xiàn)Fig.8 Gray scale distribution curve of the volume velocity relative error

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

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

5 結(jié) 語(yǔ)

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