等效源方法的適用性判定及源的優(yōu)化布置

魏應(yīng)三，王永生，沈 陽，靳栓寶，蔣國健

（1.海軍工程大學(xué)動(dòng)力工程學(xué)院，武漢 430033；2.海軍工程大學(xué)艦船綜合電子技術(shù)國防重點(diǎn)實(shí)驗(yàn)室，武漢 430033；3.中國科學(xué)院聲學(xué)所北海站，山東青島 266023）

等效源方法的適用性判定及源的優(yōu)化布置

魏應(yīng)三1，2，王永生1，沈 陽2，靳栓寶2，蔣國健3

（1.海軍工程大學(xué)動(dòng)力工程學(xué)院，武漢 430033；2.海軍工程大學(xué)艦船綜合電子技術(shù)國防重點(diǎn)實(shí)驗(yàn)室，武漢 430033；3.中國科學(xué)院聲學(xué)所北海站，山東青島 266023）

采用等效源方法建立了輻射體表面法向振速到場點(diǎn)聲壓的傳遞函數(shù)，實(shí)現(xiàn)了聲場的數(shù)值預(yù)報(bào)。以相對(duì)速度誤差為評(píng)價(jià)指標(biāo)，對(duì)等效源法的求解精度進(jìn)行了先驗(yàn)判定，提出了以1／2速度誤差限制線來確定等效源法適用的頻率范圍，進(jìn)一步以聲輻射效率為指標(biāo)對(duì)聲壓計(jì)算結(jié)果的合理性進(jìn)行后驗(yàn)判定。分別針對(duì)單個(gè)頻率點(diǎn)、多個(gè)頻率點(diǎn)對(duì)等效源的位置進(jìn)行優(yōu)化，結(jié)果表明存在最優(yōu)等效聲源位置以滿足輻射體表面法向振速的最佳匹配。將等效源的優(yōu)化方法應(yīng)用于具有聲反射面的脈動(dòng)球源及錐－柱組合殼體的輻射聲場預(yù)報(bào)，預(yù)報(bào)結(jié)果與解析解及邊界元法計(jì)算結(jié)果吻合較好，說明等效源法在工程上應(yīng)用是可行的。

等效源；優(yōu)化；聲輻射

邊界元法（BEM）是繼有限元法（FEM）之后發(fā)展起來的一種數(shù)值方法，可實(shí)現(xiàn)物理問題有效降維，但是相對(duì)于傳統(tǒng)的域方法如FEM、有限差分法，常規(guī)BEM求解方程的系數(shù)矩陣為滿陣，導(dǎo)致解題規(guī)模受到了很大限制［1］。等效源方法（ESM）的提出則促進(jìn)了BEM在工程上的應(yīng)用。Chen等［2］提出通過在輻射體表面布置等效源，并通過匹配輻射體邊界速度勢(shì)法向?qū)?shù)來預(yù)報(bào)聲場，但由于Chen將等效源布置在輻射體表面，存在積分奇異。Koopmann等［3－4］提出將等效源布置在輻射體內(nèi)部可有效克服積分奇異。于飛等［5－6］對(duì)ESM作了深入的研究，并成功應(yīng)用于聲場的全息重建。吳紹維等［7］采用無單元ESM研究了自輻射近似解析解問題，并提出等效源必須位于實(shí)際輻射體內(nèi)的一定范圍內(nèi)才能用于聲場計(jì)算。另外，ESM的求解精度與等效源數(shù)量、位置及分析頻率范圍也密切相關(guān)［8］。研究這些因素對(duì)ESM計(jì)算結(jié)果的影響的文獻(xiàn)較多，但絕大多數(shù)數(shù)值計(jì)算結(jié)果均采用相互檢驗(yàn)，較少有文獻(xiàn)對(duì)數(shù)值計(jì)算本身的可信性進(jìn)行判定，導(dǎo)致計(jì)算結(jié)果本身的可靠性又沒有判據(jù)校核。朱之墀等［9］在對(duì)管道聲場的數(shù)值預(yù)報(bào)中采用聲能守恒判據(jù)、動(dòng)量守恒判據(jù)和質(zhì)量守恒判據(jù)對(duì)計(jì)算結(jié)果進(jìn)行判定，這些判據(jù)對(duì)于指導(dǎo)數(shù)值計(jì)算具有一定借鑒作用，尤其在沒有解析解或?qū)嶒?yàn)值對(duì)比的情況下顯得更為重要。本文借鑒文獻(xiàn)［9］的思想，提出了“速度誤差小于1／2”的先驗(yàn)判據(jù)和“聲輻射效率小于1”的后驗(yàn)判據(jù)對(duì)ESM的數(shù)值計(jì)算結(jié)果的合理性進(jìn)行了判定。

1 等效源原理

等效源的主要思想是：任何輻射體的聲場可以由置于該輻射體內(nèi)若干個(gè)不同聲源強(qiáng)度的簡單源產(chǎn)生的聲場疊加［8］，這些源強(qiáng)可通過匹配輻射體表面的法向振速得到。當(dāng)?shù)刃г磸?qiáng)確定后，空間聲場任一點(diǎn)r的聲壓可由如下積分方程得到：

則當(dāng)輻射體表面的法向速度已知時(shí)，聲源強(qiáng)度表示為a＝H+Vns，此處H+為H的廣義逆。為保證聲源強(qiáng)度解的唯一性，須保證M≤N。由于等效源方法避免了聲源點(diǎn)與場點(diǎn)重合的現(xiàn)象，因此積分不會(huì)出現(xiàn)奇異。當(dāng)輻射體內(nèi)部的聲源分布確定后，根據(jù)式（1）可確定輻射體表面?D及外聲場空間D+任意一點(diǎn)rfj的聲壓，即：

對(duì)于半無限域聲輻射問題，將無限域中的格林函數(shù)G修正為半無限域中的格林函數(shù)即可。

2 等效源方法的適用性判定

2.1 先驗(yàn)判據(jù)

于飛等［8］提出了通過計(jì)算輻射體表面上的振速均方根誤差來評(píng)估求解精度，為保證ESM的求解精度提供了先驗(yàn)判據(jù)，但該方法無法確定ESM的適用頻段的截止頻率。本文在文獻(xiàn)［8］的基礎(chǔ)上，提出了一種相對(duì)振速誤差Verror來評(píng)估ESM求解精度，并確定ESM適用頻段的截止頻率，

式中，Vthero為給定的輻射體表面法向振速，Vpredict為根據(jù)等效聲源強(qiáng)度預(yù)報(bào)的輻射體表面法向振速。上式分母表示輻射體表面的等效輻射聲功率，分子表示“絕對(duì)誤差波”［8］所輻射的功率。通過求解不同頻率的振速誤差可得到誤差曲線隨頻率的變化，該曲線與“1／2誤差限制線”的交點(diǎn)所對(duì)應(yīng)的頻率即為截止頻率。本文通過對(duì)脈動(dòng)球源、一階球源、活塞球源聲輻射計(jì)算，結(jié)果表明當(dāng)Verror小于1／2時(shí)，數(shù)值解與理論解吻合較好。

球源基本參數(shù)為：球源半徑rs＝1 m，聲速1 500 m／s，密度1 000 kg／m3。等效源所在球面半徑為r0，且與脈動(dòng)球源同心，距離d＝rs－r0。分析脈動(dòng)球源時(shí)，球源以0.001 m／s的速度脈動(dòng)；分析一階球源時(shí)，球源沿x方向以0.001 m／s的速度振動(dòng)。在此算例中，球源內(nèi)部的簡單等效源的數(shù)量等于輻射體表面的節(jié)點(diǎn)數(shù)量。圖1（a）、（b）分別為脈動(dòng)球源、一階球源的速度誤差曲線，圖1（c）、（d）則為相應(yīng)的聲功率對(duì)比結(jié)果。以1／2速度誤差限制線來確定截止頻率，如圖1（a）、（b）所示。當(dāng)分析頻率遠(yuǎn)低于該截止頻率時(shí)，基于ESM計(jì)算的聲功率與理論值吻合較好，而當(dāng)分析頻率接近截止頻率時(shí)，受質(zhì)點(diǎn)速度與聲壓預(yù)報(bào)精度的雙重影響，輻射聲功率的預(yù)報(bào)值會(huì)出現(xiàn)3～6 dB的誤差，如圖1（c）中d＝0.92以及圖1（d）中d＝0.98的情況，說明采用ESM進(jìn)行聲場預(yù)報(bào)時(shí)，分析頻率上限不能太接近截止頻率。當(dāng)分析頻率超過截止頻率時(shí)，聲功率誤差顯著增加。

圖1 球源的速度誤差及輻射聲功率.（a）、（b）依次為脈動(dòng)球源、一階球源在等效源不同布置位置時(shí)的速度誤差；（c）、（d）為相應(yīng)的球源輻射的聲功率Fig.1 Velocity error（a），（b）and sound power（c），（d）of the pulsating sphere，first order sphere

Verror小于1／2的物理意義在于：當(dāng)ESM預(yù)報(bào)的振源表面等效輻射聲功率誤差小于3 dB的限制線時(shí)，認(rèn)為ESM求得的源強(qiáng)的精度能滿足下一步求解聲壓的要求，至于ESM求得的聲壓是否合理，則需根據(jù)后驗(yàn)判據(jù)來判定。

2.2 后驗(yàn)判據(jù)

對(duì)于規(guī)則輻射體聲輻射問題，其數(shù)值計(jì)算的準(zhǔn)確性可用解析解檢驗(yàn)。對(duì)于復(fù)雜結(jié)構(gòu)的聲輻射問題，在無解析解時(shí)除了用先驗(yàn)判據(jù)來判定EMS外，還需聲場的后驗(yàn)判據(jù)來檢驗(yàn)。本文提出的先驗(yàn)判據(jù)實(shí)際上是用均方振速來刻化，而當(dāng)聲壓與質(zhì)點(diǎn)速度都已知時(shí)則可用聲功率與等效輻射聲功率的比值，即輻射效率來表征。根據(jù)聲能守恒原則，輻射效率須小于1，故可用輻射效率作為后驗(yàn)判據(jù)。

以更復(fù)雜的活塞球源為例，球源基本參數(shù)同上，在球坐標(biāo)系（r，θ，φ）下，球源表面只有在張角φ≤60°范圍內(nèi)法向振速為0.001 m／s，其它范圍內(nèi)振速為0。圖2（a）所示為r0＝0.25時(shí)球源的速度誤差曲線，圖2（b）為活塞球源等效輻射聲功率、輻射聲功率理論值及預(yù)報(bào)值的對(duì)比結(jié)果。分析知，當(dāng)無量綱頻率k rs＞3.15時(shí)，ESM不滿足先驗(yàn)判據(jù)，如圖2（a），與此同時(shí)輻射聲功率預(yù)報(bào)值超過了等效輻射聲功率的理論值（123.7 dB），后驗(yàn)判據(jù)也不滿足，如圖2（b）；而當(dāng)k rs＜3.15時(shí)，ESM同時(shí)滿足先驗(yàn)判據(jù)和后驗(yàn)判據(jù)，預(yù)報(bào)值與理論值吻合較好。

圖2 活塞球源的速度誤差及輻射聲功率Fig.2 Velocity error（a）and sound power （b）of the piston source of the sphere

3 等效源的優(yōu)化布置

3.1 單目標(biāo)優(yōu)化

等效源的位置對(duì)ESM的計(jì)算結(jié)果有較大影響，不同等效源位置所對(duì)應(yīng)的截止頻率也不同。因此在工程應(yīng)用中就存在等效源位置的優(yōu)化，以使截止頻率盡可能的高。下面仍然以活塞球源為例，對(duì)給定某一頻率下的等效聲源位置進(jìn)行尋優(yōu)。假設(shè)簡單等效源所在球面與活塞球源同心，優(yōu)化變量為等效源所在的球面半徑r0，其半徑改變時(shí)，等效源的空間坐標(biāo)等比例地改變，其中等效源的數(shù)量保持不變。為避免等效聲源靠近輻射體邊界而出現(xiàn)奇異積分，此處令0＜r0／rs＜0.9。

考慮優(yōu)化問題：

令fc＝800 Hz，對(duì)應(yīng)無量綱頻率krs＝3.35，通過變r(jià)0，得到不同r0對(duì)應(yīng)的速度誤差如圖3所示。結(jié)果表明，當(dāng)0.33＜r0／rs＜0.77時(shí)速度誤差皆在1／2誤差限制線以內(nèi)，說明等效聲源布置在此范圍內(nèi)時(shí)截止頻率皆不低于3.35；當(dāng)r0／rs＝0.75時(shí)，速度誤差存在最小值，說明當(dāng)?shù)刃曉床贾迷?.75rs的球面上時(shí)可實(shí)現(xiàn)振速的最佳匹配，經(jīng)計(jì)算此時(shí)的無量綱截止頻率為3.5，說明分析頻率上限在截止頻率之前，滿足先驗(yàn)判據(jù)。為驗(yàn)證上述結(jié)論，分別取r0／rs為0.2、0.35、0.75、0.8，對(duì)活塞球源的輻射聲功率進(jìn)行預(yù)報(bào)如圖4所示。當(dāng)r0／rs取0.35（對(duì)應(yīng)的無量綱截止頻率為3.38）、0.75時(shí)，在分析頻率0＜krs＜3.35范圍內(nèi)預(yù)報(bào)值與理論值皆吻合較好，其中r0／rs取0.75時(shí)預(yù)報(bào)精度最高；當(dāng)r0／rs取0.2、0.8時(shí)，由于這兩個(gè)位置對(duì)應(yīng)的速度誤差皆大于1／2，相應(yīng)的無量綱截止頻率分別為3.02、2.94，達(dá)不到3.35，因此預(yù)報(bào)值在頻率范圍3.1＜krs＜3.35內(nèi)皆出現(xiàn)較大誤差。

圖3 krs＝3.35時(shí)單頻優(yōu)化結(jié)果Fig.3 Optimized result for single frequency krs＝3.35

圖4 不同等效源位置時(shí)的輻射聲功率Fig.4 Sound power for various locations of the source

3.2 多目標(biāo)優(yōu)化

仍然以3.1節(jié)問題為例，考慮如下優(yōu)化問題：

將該優(yōu)化問題轉(zhuǎn)化為單目標(biāo)優(yōu)化問題：

圖5 多個(gè)頻率點(diǎn)優(yōu)化結(jié)果Fig.5 Optimized result formulti frequencies

圖6 krf＝8.38時(shí)的聲指向性圖Fig.6 Sound directivity at krf＝8.38

4 應(yīng)用案例

4.1 具有聲反射面的脈動(dòng)球源的聲輻射

脈動(dòng)球源參數(shù)同2.1節(jié)中的脈動(dòng)球源，分析頻帶0～800 Hz，反射面距離球心5 m，聲反射系數(shù)為1，等效源位于與脈動(dòng)球源同心、半徑為0.01 m的球面上，輻射體采用三角形單元離散。經(jīng)計(jì)算在分析頻段內(nèi)速度誤差皆小于1／2，同時(shí)聲輻射效率小于1，滿足先驗(yàn)判據(jù)和后驗(yàn)判據(jù)。圖7為空間場點(diǎn)P（5，0，5）處的聲壓預(yù)報(bào)值與理論值，吻合較好；圖8為反射面在krs＝0.8時(shí)的聲壓分布，預(yù)報(bào)值與理論值也吻合較好，驗(yàn)證了方法的可信性。

圖7 P點(diǎn)處的聲壓曲線Fig.7 Sound pressure at point“P”

圖8 krs＝0.8時(shí)聲反射面上的聲壓分布Fig.8 Sound pressure at the reflecting plane at krs＝0.8

4.2 錐－柱組合結(jié)構(gòu)的聲輻射

圖9所示結(jié)構(gòu)由錐殼段與柱殼段組成，相關(guān)形狀參數(shù)、材料參數(shù)及激勵(lì)參數(shù)見文獻(xiàn)［10］。此例中仍然將等效聲源布置在球面上，其中球心限定在錐柱殼體的軸心線上，球面半徑以使球面不超過組合殼體為限，優(yōu)化變量為球面半徑、球心位置。經(jīng)計(jì)算，當(dāng)?shù)刃г此谇蛎姘霃饺?.15m，球心距離左端面0.3m時(shí)速度誤差達(dá)到最小值，且同時(shí)滿足先驗(yàn)判據(jù)和后驗(yàn)判據(jù)。分別采用ESM和BEM計(jì)算組合殼體的水下輻射聲功率如圖10所示，在整個(gè)分析頻段范圍內(nèi)二者吻合較好，說明本文提出的等效源優(yōu)化方法可用于水下輻射噪聲的工程預(yù)報(bào)。

圖9 組合結(jié)構(gòu)示意圖Fig.9 Scheme of the coupling structure

圖10 水下輻射聲功率Fig.10 Radiated underwater sound power

5 結(jié) 論

本文從基本的脈動(dòng)球源、一階球源、活塞球源聲輻射問題中經(jīng)驗(yàn)地提出了以1／2速度誤差限制線作為ESM適用性的先驗(yàn)判據(jù)，以聲輻射效率作為后驗(yàn)判據(jù)。最后通過兩個(gè)案例對(duì)ESM的適用性進(jìn)行了驗(yàn)證。得到如下結(jié)論：等效源位置對(duì)截止頻率影響較大，對(duì)于某一頻率（或頻段）的聲場計(jì)算，存在最優(yōu)聲源位置以滿足輻射體表面法向振速的最佳匹配，此時(shí)ESM聲場預(yù)報(bào)精度能夠滿足工程需要。

［1］Nishimura N.Fastmultipole accelerated boundary integral equationmethods［J］.Appl Mech Rev，2002，55（4）：299－324.

［2］Chen LH，Schweikert D G.Sound radiation from an arbitrary body［J］.The Journal of the Acoustical of Ameirca，1963，35（10）：1626－1632.

［3］Koopmann G H，Song L，F(xiàn)ahnline J.A Method for computing acoustic fields based on the principle of wave superposition［J］.The Journal of the Acoustical of Ameirca，1989，86 （5）：2433－2438.

［4］Song L，Koopmann G H，F(xiàn)ahnline J.Numerical errors associated with themethod of superposition for computing acoustic fields［J］.The Journal of the Acoustical of America，1991，89（6）：2625－2633.

［5］于飛，陳心昭，李衛(wèi)兵.空間聲場全息重建的波疊加方法研究［J］.物理學(xué)報(bào)，2004，53（8）：2607－2613.YU Fei，CHEN Xin-zhao，LIWei-bin.Investigatio on holographic reconstruction of sound field using wave superposition approach［J］.Acta Physica Sinica，2004，53（8）：2607－2613.

［6］李加慶，陳進(jìn)，楊超.基于波束形成和波疊加法的復(fù)合聲全息技術(shù)［J］.聲學(xué)學(xué)報(bào)，2008，33（2）：152－158.

LI Jia-qing，CHEN Jin，YANG Chao.A hybrid acoustic holography technique based on beamforming and wave superposition algorithm［J］.Acta Acustica，2008，33（2）：152－158.

［7］吳紹維，向陽，夏雪寶.基于無單元聲波疊加的自輻射近似解析表達(dá)研究［J］.振動(dòng)與沖擊，2014，33（7）：79－85.

WU Shao-wei，XIANG Yang，XIA Xue-bao.Approximate analytical expressions of self-radiation termsincluding acoustic pressure and velocity based on element free acoustic superposition［J］.Journal of Vibration and Shock，2014，33（7）：79－85.

［8］于飛.基于波疊加方法的聲全息技術(shù)與聲學(xué)靈敏度分析［D］.合肥：合肥工業(yè)大學(xué)，2005.

［9］朱之墀，李沛滋，王治國.管道聲傳播數(shù)值計(jì)算與實(shí)驗(yàn)結(jié)果比較［J］.聲學(xué)學(xué)報(bào)，1988，13（1）：1－8.

ZHU Zhi-xi，LI Pei-zi，WANG Zhi-guo.Numerical solution compared with experimental result for sound propagation in ducts［J］.Acta Acustica，1988，13（1）：1－8.

［10］魏應(yīng)三，王永生.錐－柱組合殼體聲輻射影響因素研究［J］.振動(dòng)與沖擊.2012，31（22）：174－178.

WEIYing-san，WANG Yong-sheng.Study of the effect factors on acoustical characteristics of coupled cylindrical-conical shell［J］.Journal of Vibration and Shock，2012，31（22）：174－178.

Applicability criterion of equivalent sourcesmethod and locating optim ization of sources

WEIYing-san1，2，WANG Yong-sheng1，SHEN Yang2，JIN Shuan-bao2，JIANGGuo-jian3

（1.College of Marine Power，Naval University of Engineering，Wuhan 430033，China；2.National Key Laboratery of Science and Technology on Vessel Integrated Power System，Naval University of Engineering，Wuhan 430033，China；3.Institute of Acoustics，The Chinese Academy of Sciences，Qingdao 266023，China）

Based on the equivalent sourcemethod（ESM），the transfer function from the normal vibration velocity of a radiation body's surface to the sound pressure of a field point was built and the prediction of a sound field was realized.To evaluate the accuracy of ESM in prior，the relative velocity error indexwas used and a 1／2 restrict-line of the velocity error was given to determine the frequency range for which ESM is applicable.The rationality of the sound pressure calculation results were evaluated via the sound radiation efficiency index afterward.Then the locations of the equivalent sourceswere optimized for single frequency points andmulti-frequency points，respectively.The results showed that there is an optimal sound source location to satisfy the bestmatch of the normal vibration velocity of the radiation surface.Finally，the ESM optimization algorithm was applied to predict the radiation sound fields of a pulsating sphere with a sound reflection surface and a coupled conical-cylindrical shell，the prediction results agreed wellwith those of the analyticalmethod and BEM.The results indicated that the ESM is applicable for engineering problems.

equivalent source；optimization；sound radiation

TB553；U663

A

10.13465／j.cnki.jvs.2015.21.035

國家自然科學(xué)基金項(xiàng)目（51409256，51307177，51309229）；海軍工程大學(xué)自然科學(xué)基金引導(dǎo)項(xiàng)目（HGDYDJJ13003）

2014－06－17 修改稿收到日期：2014－10－23

魏應(yīng)三男，博士，講師，1984年5月生

王永生男，博士，教授，1955年10月生郵箱：yongshengwang666＠126.com