熊濟時 吳崇健 徐志云 曾革委

中國艦船研究設(shè)計中心，湖北武漢 430064

三維數(shù)字空間法與波疊加法結(jié)合的近場聲輻射計算

熊濟時 吳崇健 徐志云 曾革委

中國艦船研究設(shè)計中心，湖北武漢 430064

傳統(tǒng)的波疊加法是邊界元法的有效替代方法，其具有計算速度快計算精度較高等優(yōu)點，但是也有復(fù)雜模型的建立較困難，以及對模型適應(yīng)性較差等缺點。為了克服傳統(tǒng)波疊加法的缺點，從而使其能夠計算較為復(fù)雜模型的聲輻射，將三維數(shù)字空間方法與傳統(tǒng)波疊加法結(jié)合，因為三維數(shù)字空間方法具有建模迅速，對模型適應(yīng)性很好等優(yōu)點。分別采用邊界表示法和體積表示法建模，且通過簡單算例計算表明，采用這兩種方法與傳統(tǒng)波疊加法結(jié)合計算外場聲輻射是完全可行的，并且結(jié)果的精確度相對較高，采用邊界表示法的效率相對較高。

聲輻射；波疊加法；邊界表示法；體積表示法

1 引言

對于振動體的外場聲輻射問題，目前最常用到的數(shù)值方法是有限元（FE）和邊界元（BE）方法［1－4］，這兩種方法極大地加速了聲輻射問題的工程發(fā)展。但它們有明顯的不足，即對不同階的奇異積分要做相應(yīng)數(shù)值處理，當(dāng)遇到具有解非唯一性的特征波數(shù)時，對該奇異積分的處理更加困難。為了尋找邊界元的有效替代方法，20世紀(jì)80年代末，Koopmann 等［5-6］和 Miller等［7］提出了基于簡單源替代的波疊加方法 （Wave Superposition Method，WSM）。其主要思想是：任何物體輻射的聲場可以由置于該輻射體內(nèi)部若干個簡單源產(chǎn)生的聲波場疊加代替，而這些簡單源的源強可以通過匹配輻射體表面上的法向振速得到。在聲輻射問題中，場中的聲壓和質(zhì)點振速必須同時滿足波動方程和輻射體表面上的邊界條件，而波疊加方法就是尋求近似解來滿足這樣的邊值問題。該方法通過在聲輻射體內(nèi)放置若干個滿足波動方程的聲源，來近似表面上的邊界條件，這種方法被證明等效于Helmholtz積分方程，也就是說和邊界元方法在理論上是一致的。該方法由于不需要求解邊界積分方程，從而避免了繁瑣的各階奇異積分處理，大大降低了數(shù)值實現(xiàn)的難度，且更易于理解和在工程上進(jìn)行推廣。

從理論上講，波疊加法適合于計算任意光滑表面物體的聲輻射，計算平臺通常選擇Matlab，因為Matlab是當(dāng)今世界上使用最為廣泛的數(shù)學(xué)軟件，它具有相當(dāng)強大的數(shù)值計算、數(shù)據(jù)處理、系統(tǒng)分析、圖形顯示，甚至符號運算功能，是一個完整的數(shù)學(xué)平臺，在這個平臺上，你只需寥寥數(shù)語就可以完成十分復(fù)雜的功能，大大提高了工程分析計算的效率［8］。對于形狀相對簡單的物體，可以很方便地在Matlab上進(jìn)行建模，但是如果物體的外形很復(fù)雜，在Matlab上建模則變得很困難。另外，如果物體外形變化，則相應(yīng)的計算程序要做較大的改動，即傳統(tǒng)的WSM方法的外形適應(yīng)性不好，為了解決這一問題，Zellers［9］和 Hwang［10］采用了一種簡單的三維數(shù)字空間方法來建模，然后運用基于實際邊界的波疊加方法來求解物體的聲輻射。這種方法的最大優(yōu)點是建模迅速，對模型外形的適應(yīng)性非常好，另外這種方法能夠較好的解決邊界奇異值的問題，但是計算精度沒有傳統(tǒng)的WSM方法高。為了充分發(fā)揮傳統(tǒng)WSM方法和三維數(shù)字空間方法的優(yōu)點，本文將傳統(tǒng)的WSM方法與三維數(shù)字空間方法結(jié)合來求解三維物體的聲輻射，即采用三維數(shù)字空間方法來建模，然后運用基于虛擬邊界的波疊加方法來求解物體的聲輻射。

2 傳統(tǒng)的波疊加方法

由線性歐拉方程和疊加積分，輻射體表面的速度可以用下式來表示：

式中，rs是S表面點的位置矢量；q（r0）是輻射體內(nèi)虛擬源分布在點 r0處的源強值；g（r，r0）為自由場的格林函數(shù)。對于復(fù)雜結(jié)構(gòu)而言，很難找到其解析解，因此有必要對其進(jìn)行離散。連續(xù)分布的虛擬源可以由N個離散的點源來取代，再在輻射體表面選取M個離散點，獲取其法向加速度值，則輻射體表面某點的法向速度可以表示如下：

式中，Qi是第i小段的源強，由N個簡單聲源可以構(gòu)建復(fù)雜振源的表面法向速度，而表面法向速度un（rs）是已知的，因此式（2）可用于計算每個簡單源的源強Qi。假設(shè)一共獲知N個點上的法向振速，那么將其構(gòu)成一個N維列向量U，并將虛球源上的N個點源的源強構(gòu)成一個N維列向量Q，將它們之間的關(guān)系寫成矩陣形式：

D是M×N的傳遞矩陣，其表達(dá)式如下：

當(dāng)M＞N時，即可獲得方程的唯一解，獲得虛擬點源源強之后，再通過式（5）可以計算得到任意點的聲壓：

其中，M為單極矩陣，其與自由格林函數(shù)成正比。

3 三維數(shù)字空間方法

數(shù)字表達(dá)意味著計算機圖形，所有的圖像都用離散的像素來表示，而無需真實形狀的數(shù)值或理論信息，這意味著為了在計算機域中表示一個物體，必須將物體轉(zhuǎn)化為離散的形式。三維數(shù)字空間方法即是以離散像素來表示具體的物體。具體實現(xiàn)的方法有兩種，一種是邊界表示法，即將表面表示成一系列的三角形或者四邊形的細(xì)小平面；另一種是體積表示法，即將物體的體積以細(xì)小的體積單元代替。如圖1所示，分別給出了一個圓球的兩種表示方法。

上圖給出了三維物體的數(shù)字表示信息，而在波疊加法中只需要給出物體表面的信息，因此首先需要描述物體的表面，與通常的有網(wǎng)格方法不一樣，邊界表示法只是提取物體表面的節(jié)點信息來描述物體的表面，而無需考慮表面單元的信息；而體積表示法通常有兩種描述方法，一種是直接從劃分網(wǎng)格的軟件導(dǎo)入，例如從ANSYS、PATRAN、I－DEAS等；另一種是采用理論方法來描述一個物體的表面，通常是一些比較規(guī)則的表面如球面、橢球面等。兩種方法都需要設(shè)置一個如圖2所示的閾值（閾值越大則選中的節(jié)點數(shù)越多，所描述的表面卻越不真實，閾值越小則選中的節(jié)點數(shù)越小，但是所描述的表面卻越真實，但是如果節(jié)點數(shù)太少，計算精度則可能受到影響），將符合閾值要求的節(jié)點認(rèn)定為表面節(jié)點。圖2給出了第二種方法確定物體表面的方法，圖中紅色帶*點組成的圖形用來近似描述球表面。Zellers［5］和Hwang采用了體積表示法來建模，然后運用基于實際邊界的波疊加方法來求解物體的聲輻射。

4 數(shù)值計算

由于脈動圓球的聲輻射有解析解，因此本文采用脈動圓球作為計算模型，所有計算均在MATLAB7.0平臺完成，圓球半徑取為1 m，表面脈動速度為1 m／s，計算了距離實際邊界表面1 m距離處的聲輻射 （本文定義的近場為距離實際邊界表面1 m處），分別采用邊界表示法和體積表示法來求解。

4.1 邊界表示法

由于在文獻(xiàn)［8］中已經(jīng)做過虛擬源半徑與實際邊界半徑之間關(guān)系的研究，研究表明虛擬源半徑取到實際邊界的0.2倍時，即可滿足計算要求。因此本文首先通過ANSYS軟件前處理器建立半徑分別為1 m和0.2 m的兩個物理模型，并將兩個模型分別劃分網(wǎng)格，分別在ANSYS中提取兩個模型的外表面節(jié)點坐標(biāo)，半徑為1圓球的節(jié)點坐標(biāo)定義為（X，Y，Z），它表示實際邊界的節(jié)點坐標(biāo)；半徑為 0.2 m 圓球的節(jié)點坐標(biāo)定義為 （x，y，z），而它則表示虛擬邊界的節(jié)點坐標(biāo)；然后分別將其保存為兩個文本文檔 1.txt和 2.txt；再將這些節(jié)點信息導(dǎo)入MATLAB；最后根據(jù)公式（4）計算得到虛擬源與振動源表面之間的振速傳遞系數(shù)矩陣D，根據(jù)式（3）計算得到虛擬源的源強值，根據(jù)式（6）得到單極系數(shù)矩陣，就可以根據(jù)式（5）算出外部輻射聲場。

4.2 體積表示法

本文與Zellers和Hwang一樣，也采用體積表示法建模，但是采用傳統(tǒng)的WSM方法來求解聲輻射。首先將三維空間劃分為細(xì)小的晶格，然后將在網(wǎng)格文件導(dǎo)入或者直接用理論方法描述的物體置于劃分好的三維空間當(dāng)中，設(shè)置閾值，找出符合閾值條件的節(jié)點；再將符合條件的所有節(jié)點按一定比例縮小，使新得到的節(jié)點位于實際物體邊界內(nèi)部，作為虛擬邊界節(jié)點；最后根據(jù)傳統(tǒng)的WSM方法生成矩陣并且求解得出外部輻射聲場。

4.3 結(jié)果及討論

首先我們給出脈動圓球外場聲輻射的解析解（無量綱聲壓）表達(dá)式如下：

式中，a為圓球的半徑；r為外場中某點與球心之間的距離（r＝2 m）；p（r）為距離 r處的聲壓；k 是波數(shù)。再分別按照上述步驟，采用邊界表示法和體積表示法也計算出脈動圓球的外場聲輻射。

假定數(shù)值計算得到的無量綱聲壓實部為per，無量綱聲壓虛部為pei，對應(yīng)的無量綱聲壓解析解分別為pr和pi，則實部誤差εr和虛部誤差εi分別為：

采用邊界表示法時，分別考慮了兩種情況，第一種情況是半徑為1 m球表面包含114個節(jié)點，半徑為0.2 m球表面也包含114個節(jié)點；第二種情況是半徑為1 m球表面包含114個節(jié)點，而半徑為0.2 m球表面包含78個節(jié)點，計算結(jié)果如圖3所示。圖3a為脈動圓球數(shù)值解與解析解之間的誤差隨波數(shù)變化實部，圖3b為脈動圓球數(shù)值解與解析解之間的誤差隨波數(shù)變化虛部，數(shù)值解的實部和虛部在波數(shù)為0～20范圍內(nèi)基本上都與解析解相吻合，誤差很小。虛擬邊界表面節(jié)點選取較多時，計算精度較高，但是在ka＝15.7≈5*pi時，最大實部誤差達(dá)到了1%左右，最大虛部的誤差更是達(dá)到了10%以上，這是因為脈動圓球的特征波數(shù)為pi的整數(shù)倍。

采用體積表示法時，將三維空間劃分為邊長為0.02 m的立方體格子，然后將半徑為1 m球置于三維空間當(dāng)中，設(shè)置閾值為10－6，找出符合閾值條件的表面節(jié)點，共150個，然后將這些表面節(jié)點按0.2的倍數(shù)縮減，重新得到150個小球的表面節(jié)點，最后將半徑為1 m球表面的節(jié)點作為實際邊界節(jié)點，而半徑為0.2 m球表面邊界節(jié)點作為虛擬邊界節(jié)點，分別代入傳統(tǒng)波疊加方法，求解出脈動圓球的聲輻射。圖4給出了采用體積表示法時，數(shù)值解和解析解之間誤差的實部和虛部。結(jié)果同樣示出了波數(shù)為0～20范圍內(nèi)的誤差值，從圖中可以看出，實部和虛部的誤差值都非常小，均在0.1%以下，結(jié)果的精度要顯著高于Hwang給出的結(jié)果。

從圖3和圖4來看，只要邊界節(jié)點選擇恰當(dāng)?shù)脑挘吔绫硎痉ê腕w積表示法均能得到較好的結(jié)果。邊界表示法只需要構(gòu)造輻射體的表面，在進(jìn)行數(shù)字空間構(gòu)造時網(wǎng)格處理量較少，而體積法則需要構(gòu)造出整個輻射體，在進(jìn)行空間構(gòu)造時處理的網(wǎng)格數(shù)通常非常大，構(gòu)造工作量相對較大（本文當(dāng)中的晶格數(shù)量為50×50×50＝125 000），因此邊界表示法比體積表示法效率更高。

5 結(jié)束語

傳統(tǒng)的WSM方法是一種聲輻射的無網(wǎng)格計算方法，其具有計算速度快，精度較高的優(yōu)點，但是其外形適應(yīng)性不是太好，對于復(fù)雜模型的建立難度較大，因此對于工程推廣應(yīng)用有一定的限制。三維數(shù)字空間方法則具有建模迅速，對模型外形的適應(yīng)性非常好等優(yōu)點。本文將兩者結(jié)合用于求解外場聲輻射，并以具有解析解的脈動圓球聲輻射算例來驗證了方法的可行性，以及此種計算方法的結(jié)果具有較高精度這一特性。

［1］譚林森，駱東平，吳崇健，等.潛水器動力艙振動與聲輻射［J］.華中理工大學(xué)學(xué)報，1999，27（11）：7－9.

［2］徐張明，沈榮瀛，華宏星.利用FEM／IBEM計算流體介質(zhì)中的殼體的結(jié)構(gòu)聲耦合問題 ［J］.振動工程學(xué)報，2002，15（3）：363－367.

［3］徐張明，汪玉，華宏星，等.船舶結(jié)構(gòu)的建模及水下振動和輻射噪聲的FEM ／BEM計算 ［J］.船舶力學(xué)，2002，6（4）：89－95.

［4］商德江，何祚鏞.加肋雙層圓柱殼振動聲輻射數(shù)值計算分析［J］.聲學(xué)學(xué)報，2001，26（3）：193－201.

［5］KOOPMANN G H，SONG L，F(xiàn)AHNLINE J B.A method for computing acoustic fields based on the principle of wave superposition［J］.Journal of Acoustic Society America，1989，86（6）：2433－2438.

［6］SONG L，KOOPMANN G H，F(xiàn)AHNLINE J B.Numeric errors associated with the method of superposition for computing acoustic fields ［J］.Journal of Acoustic Society America，1991，89（6）：2625－2633.

［7］MILLER R D.A comparison between the boundary element method and the wave superposition approach for the analysis of the scattered fields from rigid bodies and elastic shells ［J］.Journal of Acoustic Society America，1991， 89（2）：2185－2196.

［8］熊濟時，吳崇健，曾革委，等.基于波疊加法的近場聲輻射研究 ［C］//2009年第22屆全國振動與噪聲高技術(shù)及應(yīng)用學(xué)術(shù)會議，2009：135－140.

［9］ZELLERS B C.An acoustic superposition method for computing structural radiation in spatially digitized domain［D］.USA：The Pennsylvania State University， 2006.

［10］ HUANG Y S.A wave superposition method formulated in digital acoustic space ［D］.USA：The Pennsylvania State University，2009.

Near-Field Sound Radiation Numeration Based on Wave Superposition Method in the 3D Digital Space

Xiong Ji－shi Wu Chong－jian Xu Zhi－yun Zeng Ge－wei
China Ship Development and Design Center， Wuhan 430064， China

The traditional Wave Superposition Method （WSM） may be substituted for the Boundary Element Method （BEM） because of its capabilities of rapid computation and greater precision.However，with merely reliance on mesh geometry，it cannot accommodate fast shape changes in the design stage of a consumer＇s product or machinery， where a great amount of iterations of shape changes are required.In order to overcome these shortcomings of WSM，a new approach to representing geometry was introduced by constructing a uniform lattice in the 3D digital work space.With this method， geometry was represented with boundary nodes or small lattice that can easily adapt to shape changes and therefore is more suitable for shape optimization.The boundary representation and volume elements representation method respectively combined with WSM were validated by computing sound radiation of simple structures.The results show that both methods are entirely capable of handling this calculation，and have higher precision， in addition， the boundary representation method exhibits higher efficiency than volume element representation method.

sound radiation； wave superposition method； boundary representation method； volume elements representation method

O422.6

A

1673－3185（2011）01－41－05

10.3969/j.issn.1673-3185.2011.01.008

2010-03-12

船舶工業(yè)國防科技預(yù)研基金（08J1.4.3）

熊濟時（1977－），男，博士研究生。研究方向：噪聲與振動控制。E-mail：xiong26@163.com

吳崇?。?960－），男，研究員，博士生導(dǎo)師。研究方向：噪聲與振動控制