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

熊濟(jì)時(shí) 吳崇健 徐志云 曾革委

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

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

熊濟(jì)時(shí) 吳崇健 徐志云 曾革委

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

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

4 數(shù)值計(jì)算

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

4.1 邊界表示法

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

4.2 體積表示法

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

4.3 結(jié)果及討論

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

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

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

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

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

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

5 結(jié)束語

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

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

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

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

［4］商德江，何祚鏞.加肋雙層圓柱殼振動(dòng)聲輻射數(shù)值計(jì)算分析［J］.聲學(xué)學(xué)報(bào)，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］熊濟(jì)時(shí)，吳崇健，曾革委，等.基于波疊加法的近場(chǎng)聲輻射研究 ［C］//2009年第22屆全國(guó)振動(dòng)與噪聲高技術(shù)及應(yīng)用學(xué)術(shù)會(huì)議，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è)國(guó)防科技預(yù)研基金（08J1.4.3）

熊濟(jì)時(shí)（1977－），男，博士研究生。研究方向：噪聲與振動(dòng)控制。E-mail：xiong26@163.com

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