基于區域分解的快速多極基本解法預測發動機表面輻射聲場

第一作者陳長征男，教授，博士生導師，1964年生

通信作者張士偉男，博士生，1983年生

基于區域分解的快速多極基本解法預測發動機表面輻射聲場

陳長征1，張士偉1, 周勃2, 黃鶴艇3

(1. 沈陽工業大學機械工程學院，沈陽110870；2.沈陽工業大學建筑工程學院,沈陽110870；3.廈工機械股份有限公司, 福建 廈門361000)

摘要：針對發動機表面結構復雜，多狹長邊界，用通常算法求解源點分布無法保障基本解法計算精度問題，提出改進的計算源點分布算法，利用區域分解處理影響基本解計算精度的狹長邊界；結合快速多極算法，形成適于三維復雜表面輻射聲場預測基于區域分解的快速多極基本解法。利用區域分解求解發動機計算模型的源點分布，在源點分布已知基礎上利用快速多極基本解法預測發動機表面輻射聲場。結果表明，利用區域分解求得源點與配置點最優距離為34.5 mm，可提高發動機輻射聲場預報精度；截斷項數為15時可避免截斷項數過小引起的低頻不穩定，并保證計算效率；不同于快速多極邊界元用網格加密提高計算精度，基于區域分解的快速多極基本解法通過區域分解進行優化源點分布可提高精度，且不增加計算負擔。通過發動機表面輻射噪聲實驗，計算值與實驗值吻合較好，說明該解法可提高發動機表面輻射聲場的預測精度。

關鍵詞：快速多極基本解法；聲輻射；發動機；截斷誤差；聲場測量

基金項目：中國博士后科學基金(2014M560220)；國家科學支撐計劃(2013BAF07B04)

收稿日期：2014-12-18修改稿收到日期：2015-03-13

中圖分類號:TB532；TK4文獻標志碼：A

Prediction of acoustic radiation field on engine surface by domain decomposition based fastmultipole method of fundamental solution

CHENChang-zheng1,ZHANGShi-wei1,ZHOUBo2,HUANGHe-ting3(1. School of Mechanical Engineering, Shenyang University of Technology, Shenyang 110870, China；2.Architectural and Engineering Institute, Shenyang University of Technology, Shenyang 110870, China；3. Xiamen XGMA Machinery Co., Ltd, Xiamen 361000, China)

Abstract:Due to the complex structure of engine surface and a lot of long narrow borders, standard approaches to solve the distribution of source points cannot gain enough computational precision in using the method of fundamental solution(MFS). Thus, an improved algorithm for solving the distribution of source points was presented and used to deal with long narrow borders. Combining with the fast multipole algorithm, a domain decomposition based fast mutipole method of fundamental solution (FMMFS) was proposed. The improved algorithm was used to solve the distribution of source points, while the FMMFS was used to predict the acoustic radiation field on engine surface. The results show that the optimal distance between the collocation points and source points is 34.5 mm and the computational precision is improved. When the number of truncation terms is 15, the uncertainty in low frequency due to too few number of truncation terms can be avoided and the computational efficiency can be guaranteed. To improve the computational precision, the method uses domain decomposition to solve the distribution of source points, it differs from the fast multipole BEM which adopts mesh refinement. The experimental results are in good agreement with the computational results, which indicate that the domain decomposition based FMMFS can improve substanlially the computational precision of the prediction of acoustic radiation field on engine surface.

Key words:fast multipole method of fundamental solution; acoustic radiation; engine; truncation term; sound field measurement

發動機的空氣動力噪聲可由空氣濾清器、消聲器[1]有效控制，輻射噪聲成為發動機主要噪聲源。預測發動機表面輻射聲場特性對噪聲控制非常重要。結構輻射聲場預測的數值方法包括有限元、邊界元、基本解法等。預測輻射聲場通常用有限元及邊界元法[2]。有限元法可求解任意形狀結構的輻射聲場，對工程中無限域聲輻射問題通常利用有限區域截斷模擬聲場無限邊界，而因人為確定邊界條件，會導致計算誤差。邊界元法基本方程通解自動滿足Sommerfeld聲輻射條件，積分區域只含流固交界面的勢函數，能有效降低求解維數。對在特征波數處解不唯一性問題，需進行冗繁的插值求解，耗時長。當場點位于輻射體表面時需處理格林函數不同階次的奇異積分。Burton等[3]提出任意頻率下求得唯一解方法，然而Helmholtz 積分方程中法向導數引入了超奇異積分，導致計算困難；改進的 Burton-Miller方法能克服Burton-Miller方法不足，在任意頻率下可求得唯一解。

與邊界元相比，基本解法(Method of Fundamental Solutions,MFS)可通過源點及配置點合理配置避免核函數奇異積分。Aleksidze等[4]提出基本解法，將邊界離散點直接布置在自然邊界上,將源強點布置在求解域外，避開核函數的奇異積分，并用于諧函數理論的混合邊界問題。該方法用于聲場預測時形成的系統矩陣是滿秩、非對稱矩陣，存儲量大、求解速度慢，工程應用受到限制。Greengard等[5]提出快速多極算法(Fast Multipole Algorithm)，可解決傳統基本解法求解效率低的問題，并用于邊界元法中。Shen等[6]利用基于自適應樹狀結構的快速多極邊界元法計算發動機聲輻射。吳海軍等[7]提出基于Burton-Miller方程的三維多層聲學快速多極子邊界元算法，并用于車內聲場的分析計算。季振林等[8]將快速多極邊界元用于計算吸聲材料聲場。李善德[9]提出寬頻快速多極邊界元法求解大空間聲學問題。Liu等[10]將快速多極算法與基本解法結合，推導出針對二維勢問題快速算法的理論公式。Jiang等[11]采用快速多極算法中多極展開、多極展開轉移，實現用基本解法求解二維Helmholtz問題。陳劍等[12-13]將統計方法與波疊加法結合提出統計波疊加法，能對隨機聲場準確計算；引入快速多極子算法，提出快速多極基本解法(FMMFS，Fast Multipole Method of Fundamental Solutions)用于二維類車體輻射模型仿真分析。吳紹維等[14]采用奇異點挖去法、部分積分區域替換法及不變量嵌入法推導出自輻射聲壓項及速度項的非奇異表達，解決聲源點與接收點重合時格林函數的奇異性，實現中低頻聲場精確計算。

傳統基本解法預測發動機輻射聲場關鍵為如何布置源點，源點分布不同會直接影響計算精度[15]。傳統波疊加法中源點須置于輻射體內離表面一定距離處，否則會使算法中的格林函數產生奇異性，無法準確計算聲場。通常將源點布置在虛球面或與邊界類似的虛源面上，并通過求解兩個面之間最優距離獲得源點分布[16]。對形狀復雜邊界如發動機表面，存在較多狹長邊界，用此法布置源點會導致較大計算誤差。為此，本文提出改進算法用于處理狹長邊界計算源點分布，結合快速多極算法形成適于三維復雜表面輻射聲場預測基于區域分解的快速多極基本解法，并用于某發動機表面輻射聲場預測分析，提高發動機的表面輻射聲場預測精度，并實驗驗證。

1基于區域分解的快速多極基本解法

1.1快速多極基本解法

各向同性均勻介質中，穩態聲場Helmholtz方程為

(1)

圖1　求解域與源點布置示意圖 Fig.1 The diagram of the solution domain and the distribution of source points

對半空間聲輻射問題[17]基本解法中聲壓可近似為

(2)

對式(2)進行離散，xi(i=1,2,3,…,M)為邊界?Ω上配置點，從而有

(3)

式中：i=1, 2, 3, …,M。

得標準線性方程組為

(4)

即

Aq=b

式中：A為系數矩陣；q為未知向量；b為已知向量。

本文源點數目與配置點數目相等，即M=N。

快速多極子算法的主要為通過多極展開、局部展開以及相應展開轉移過程。考慮離配置點x較遠的源強點yi，將式(3)多極展開為

(5)

(6)

為將晶格Y1的多極展開系數轉移其父層晶格中心點Y0，多極展開系數轉移式(Moments to Moments,M2M)為

(7)

Wn,n′,m,m′,l=

式(3)沿X0局部展開，即

(8)

將晶格Y0的多極展開系數變換為晶格X0的局部展開系數，多極展開到局部的變換式(Moments to local translation, M2L)為

(9)

將晶格X0的局部展開系數轉移為父晶格X1局部展開系數，局部展開系數轉移式(Local to Local translation, L2L)為

(10)

1.2區域分解

源點與配置點最優距離計算步驟[19]為

步驟1: 設源點數目為N，配置點數目均為M，計算步長為δd。

步驟3:在集合ej中選最小值記為ej*，故dopt=jδd。

對求解域邊界簡單情形此方法較精確，但對復雜邊界(圖1)狹長區域，求解最優距離d時其取值會受狹長區域限制。由此產生的最優值會產生較大誤差。為提高計算精度分解整個區域，將影響計算精度的狹長區域作為一個區域，尋找最優距離時并非針對整個區域，而應針對離散后各狹長區域進行最優值求解獲得局部最優值，進而獲得全局最優值。

對每個含狹長邊界的立方體，利用求解最優距離算法進行求解獲得局部最優值，進而獲得全局最優值。據最優值dopt、各區域配置點及最優距離可求得源點，即在邊界?Ω上每個配置點做法線，沿法線行進距離dopt后即獲得該配置點源點。求得源點分布后利用快速多極基本解法求解輻射聲場，此即本文基于區域分解的快速多極基本解法，求解外場聲輻射步驟為：①據改進算法計算獲得源點分布，針對半空間模型重新對八叉樹結構分層。建立含全部源點及邊界配置點的立方體(晶格)，將該立方體作為樹結構第0層；八等分該立方體，對等分的每個子立方體進行相同操作，直到立方體所含單元數為給定值為止。②上行歷遍計算多極展開系數。先由式(6)計算所有葉子的多極展開系數，再利用式(7)逐層向上計算樹結構中所有父層晶格的多極展開系數。③下行遍歷計算局部展開系數。由式(9)、(10)求出所有晶格的局部展開系數。④迭代求解。對給定的初始向量q0，通過③～⑤采用廣義極小殘差法(GMRES)求解器迭代求出問題解q。⑤由式(2)求出場點聲壓。

2發動機聲輻射分析相關參數計算

基于區域分解的快速多極基本解法預測發動機輻射聲場，由配置點求出源點后據速度邊界條件即可求得發動機的輻射聲場。

2.1發動機模型速度邊界條件與源點布置

發動機模型為直列兩缸柴油機250，發動機表面法向振速分布云圖見圖2。其最大值分布為圖中右上部分。將該振速保存為PCH文件，以便為配置點提供初始邊界條件。

圖2　發動機表面法向振速分布云圖 Fig.2 Distribution of normal velocity on engine surface

鑒于發動機表面較復雜且存在狹長區域，據改進算法將發動機模型分割為兩個區域，見圖3(a)，區域1、2為發動機模型分割的兩區域。該分區理由為區域2的邊界表面振速較大(圖2法向振速分布云圖)，區域1法向振速較區域1小。計算時將區域1中表面狹長區域適當簡化，使區域1表面較平整。發動機表面配置點分布示意圖見圖3(b)，發動機模型劃分三角單元21 754個，配置點10 859個。用改進算法計算源點與配置點最優距離，據圖3中配置點及最優距離即可求得源點。此源點位于求解域外部即發動機內部，當源點與配置點重合時會出現奇異性問題，通過以上源點配置方法可有效克服用傳統基本解法出現的奇異性問題。

圖3　發動機表面配置點分布示意圖(左視圖) Fig.3 Distribution of collocation points on engine surface (left view)

2.2截斷項數選取

選擇含發動機模型所有邊界的最小立方體為樹結構第0層，將第0層等份成8個立方體為第1層，每個立方體按以上分法劃分直到每個立方體所含單元數小于25停止，本文樹結構層數向下分劃到第5層。

在發動機聲輻射計算中截斷項數取值取決于計算的頻率范圍，由經驗公式[20]確定，即

kDl

(11)

式中：k為計算波數；Dl為第l層劃分立方體直徑；rlm為多極擴展點到局部擴展點距離。

計算精度隨截斷項數增大而提高，但截斷項數并非越大越好，需據問題規模及精度要求選擇適當的截斷項數。發動機模型計算頻率為400 Hz時波數k為1.167，D5為5.74。由式(11)知，截斷項數T大于7即可。為避免因截斷項數取值過小，引起低頻不穩定及計算精度下降，本文截斷項數T的取值為14～16。工程計算分析中在滿足式(11)時截斷項數盡量取大些，粗略獲得截斷項數取值后據計算能力及精度取舍。

3結果分析與實驗驗證

3.1發動機聲輻射計算精度及結果分析

利用改進算法，步長δd=2 mm計算源點與配置點最優距離。設源點數目與配置點數目相等均為10 859，將距離d作為ej連續函數，最大誤差為

式中：j=1,2,…,N，N為源點數；M為配置點數。

據改進算法可求得發動機模型的源點分布。算法改進前后源點、配置點間距離d對應的最大誤差見圖4。由圖4可知，源點與配置點間距離對計算誤差影響較大，距離d變小時計算誤差變大；d趨于零時源點及配置點重合，此時因奇異性問題計算誤差較大；理論上源點與配置點距離d增大時，計算誤差呈指數趨于零。隨d增大式(4)中矩陣會呈病態導致計算誤差變大。改進前最優值d=55 mm；距離d繼續變大時計算誤差陡增，因發動機表面復雜且存在狹長區域，過大距離d直接導致計算誤差變大。為減小狹長區域對計算誤差影響，改進后算法采取分區域求解，所得最優距離d=34.5 mm，最優距離d較改進前最優距離小，且計算誤差有降低。

圖4 計算源點分布算法改進前后比較Fig.4Comparisonofstandardapproachandimprovedalgorithmforthedistributionofsourcepoints圖5 最優值d=34.5mm時最大誤差與頻率對應關系Fig.5Themaximumerrorandfrequencywiththeoptimizedvalued=34.5mm圖6 截斷項數對迭代誤差影響Fig.6Effectsoftruncationtermsoniterationerrors

最優值d=34.5 mm時最大誤差與頻率對應關系見圖5。由圖5看出，頻率400 Hz附近最大誤差較小。而頻率越大計算誤差亦越大。因此，計算針對400 Hz附近頻段進行。

用于計算發動機聲輻射預測的式(5)～式(10)在多級展開式中均舍棄若干項，必會影響計算的收斂性及精度，為便于研究比較采用已給固定截斷項數。頻率為400 Hz時分層樹結構不同項數對計算收斂性影響見圖6，該圖揭示出截斷項數對計算效率的影響規律。截斷項數T=14時迭代20次后收斂速度明顯放緩，收斂殘差保持在約7.497×10-4，收斂差、精度低；截斷項數T=15時35次迭代即達精度要求；截斷項數T=16時26次迭代達到精度要求。表明收斂速度隨截斷項數增加而提高，但并非截斷項數越大越好，此因截斷項數過大會增加計算復雜性。

為說明改進算法提高計算精度的有效性，取截斷項數T=15，用俯視圖比較采用改進算法前后發動機模型輻射聲場分布，見圖7。由圖7看出，聲壓級較大區域分別分布在發動機前后端對應的輻射面及發動機頂部輻射面；而聲壓級較低區域分布在發動機兩側對應的輻射面。與圖7(a)相比圖7(b)聲場分布更精細，計算精度更高，由于改進算法采用區域分解求解模型的源點分布，可減少影響模型中狹長區域對計算精度影響。

圖7　T=15時發動機模型輻射聲場(俯視圖) Fig.7 Acoustic radiation field of the model of engine with T=15(top view)

快速多極邊界元法與本文方法均采用快速多極子算法，均能極大提高計算效率。由于邊界元及基本解法具有不同的核函數與積分域，導致兩種方法用于輻射聲場預測精度、效率存在差異。快速多極邊界元法針對所求輻射體整個邊界面進行離散求解；而快速多極基本解法據求解域外虛源面的若干源點進行求解，能極大降低求解矩陣規模。對具有復雜邊界的輻射聲場求解，快速多極邊界元需更精細的網格劃分，加重計算負擔；而快速多極基本解法只需增加若干個源點即可，減少計算量。因此，對表面復雜的發動機輻射聲場預測，快速多極基本解法計算效率高。

快速多極邊界元通過網格加密及增加截斷項數提高求解精度。雖快速多極基本解法也可通過增加截斷項數提高求解精度，在相同截斷項數下與快速多極邊界元采用網格加密提高精度不同，快速多極基本解法通過區域分解進行優化源點分布提高精度。

基于快速多極邊界元法的聲學分析大多采用常量單元插值計算，可直接求解強奇異積分，但難以獲得高計算精度結果，尤其當求解區域邊界復雜多狹長邊界時計算精度受影響，而本文方法采用區域分解處理影響計算精度的狹長邊界，并未增加計算負擔，且能提高計算精度。

3.2實驗分析

測量選室外較空曠處、風小天氣下測量。平均大氣壓1 055.9 hPa，平均溫度25℃。測量發動機輻射噪聲儀器為B&K公司聲級計B&K2250及UA-1650 90 mm自動調節防風球。為降低背景噪聲影響，測量選夜間進行。發動機為直列兩缸柴油機250，為排除排氣噪聲影響，將排氣口引到較遠處半封閉空間。據國際測量標準ISO3744，測點布置見圖8。圖8(a)為正視圖，圖8(b)為俯視圖，發動機位于中心，所選19個測點均布于半球面。

截斷誤差T=15時實驗值與計算值比較見圖9。由圖9看出，實驗值與計算值吻合較好，輻射聲壓較大值出現在位于發動機前上方和左前上方的測點4、13、14、15、16，原因為發動機前上方振速較大(圖2)。最高處達71 mm/s，最大振速分布于該區域，輻射的噪聲與其它部位噪聲疊加到該測點所在區域，引起測點輻射聲壓較大。測點9輻射聲壓最低為77.1 dB，其所在區域輻射聲壓較低。圖9中實驗值與計算值總體上一致，從而證明本文基于區域分解的快速多極基本解法計算發動機聲輻射問題的正確性。

圖8　測點布置 Fig.8 Plot of measurement points

實驗值與計算值在測點14、15、16、18相差較大，甚至達到3 dB，因發動機前上方表面結構復雜，源點不同分布影響聲輻射問題計算精度。為說明截斷項數對發動機聲輻射影響，不同截斷誤差下計算值與實驗值比較見圖10。由圖10看出，隨截斷項數增加計算值逐漸向實驗值所在直線靠近，說明計算精度隨截斷項數增加而提高。

圖9 T=15時各測點計算值與實驗值比較Fig.9ComparisonofcomputationvaluesandexperimentalvaluesonmeasurementpointswithT=15圖10 不同截斷項數下計算值與實驗值比較Fig.10Comparisonofcomputationvaluesandexperimentalvaluesonmeasurementpointswithdifferenttruncationterms圖11 算法改進前后計算值與實驗值比較Fig.11Comparisonofcomputationvaluesofstandardapproachandimprovedalgorithmwithexperimentalvalues

T=16時計算值離實驗值所在中心線最近，表明截斷項數取該值時計算精度較高，但截斷項數過大會降低快速多極方法的計算效率；T=14截斷項數取值較小時計算效率更高，計算精度會變差，亦會引發低頻計算值不穩定問題。因此需據問題規模及精度要求選擇適當的截斷項數，本文取T=15。算法改進前后計算值與實驗值比較見圖11。由圖11看出，圖8中19個測點的實驗值分布在圖11橫、縱坐標的角平分線上，改進前、后測點處聲壓級值分布于兩側。改進后計算所得19個測點值分布在角平分線較近，而改進前計算所得測點值分布離角平分線較遠，說明改進后計算的各測點值與實驗值更接近，基于區域分解的快速多極基本解法較原算法計算發動機輻射聲場更精確。

4結論

針對求解域因存在狹長邊界損失計算精度問題提出改進算法計算基本解法的源點分布，結合快速多極基本解法形成基于區域分解的快速多極基本解法。用區域分解計算發動機源點分布，并將基于區域分解的快速多極基本解法用于預測發動機表面輻射聲場，結論如下：

(1)用改進算法即區域分解可求得發動機表面上配置點與源點的最優距離，不僅能避免該兩點距離過小引起積分奇異性、過大導致矩陣病態，且能提高發動機表面輻射聲場的預測精度。

(2)截斷項數過小會引起低頻不穩定，而過大，雖計算精度高，但會影響計算效率。

(3)相同截斷項數下，基于區域分解的快速多極基本解法通過區域分解進行優化源點分布，可提高發動機表面輻射聲場精度，且不增加計算負擔。

參考文獻

[1]馮仁華,張大鳴,鄧幫林,等. 一種高效的發動機輻射噪聲計算方法研究[J]. 振動與沖擊,2014,33(18):198-203.

FENG Ren-hua, ZHANG Da-ming, DENG Bang-lin, et al. A fast and high accuracy calculation method for engine radiated noise[J].Journal of Vibration and Shock, 2014,33(18):198-203.

[2]趙志高,黃其柏. 復雜結構的聲輻射解耦及其聲輻射效率分析[J]. 振動工程學報,2004,17(3):326-331.

ZHAO Zhi-gao, HUANG Qi-bai. Analysis of acoustical radiation decoupling and acoustical radiation efficiency of complex structure[J]. Journal of Vibration Engineering, 2004, 17(3): 326-331.

[3]Burton A J, Miller G F. The application of integral equation methods to the numerical solution of some exterior boundary value problems[J]. Proceedings of the Royal Society of London Series A, 1971, 323:201-210.

[4]Kupradze V D, Aleksidze M A. The method of functional equations for the approximate solution of certain boundary value problems[J].Computational Mathematics and Mathematical Physics,1964,4(4):82-126.

[5]Greengard L, Rokhlin V. A fast algorithm for particle simulations[J]. Journal of Computational Physics, 1987,73(2): 325-348.

[6]SHEN L, Liu Y J. An adaptive fast multiple boundary element method for three-dimensional acoustic wave problems based on the burton-miller formulation[J].Computational Mechanics, 2007, 40(3): 461-472.

[7]吳海軍,蔣偉康,魯文波. 三維聲學多層快速多極子邊界元及其應用[J]. 物理學報,2012,61(5):1-8.

WU Hai-jun, JIANG Wei-kang, LU Wen-bo. Multilevel fast multipole boundary element method for 3D acoustic problems and its applications[J]. Acta Phys. Sin.,2012, 61(5): 1-8.

[8]崔曉兵,季振林. 快速多極子邊界元法在吸聲材料聲場計算中的應用[J]. 振動與沖擊,2011,30(8):187-193.

CUI Xiao-bing, JI Zhen-lin. Application of FMBEM to calculation of sound fields in sound-absorbing material [J]. Journal of Vibration and Shock, 2011, 30(8):187-193.

[9]李善德.大規模聲學問題的快速多極邊界元方法研究[D].武漢：華中科技大學,2011.

[10]Liu Y J, Nishimura N, Yao Z H. A fast multiple accelerated method of fundamental solutions for potential problems [J].Engineering Analysis with Boundary elements, 2005,29(11):1016-1024.

[11]Jiang X, Chen W, Chen C S. A fast method of fundamental solutions for solving Helmholtz-type equations[J]. International Journal of Computational Methods,2013, 10(2): 89-108.

[12]陳劍,高煜,許濱,等. 計算結構隨機振動輻射聲場的統計波疊加方法[J]. 振動工程學報,2009,22(3):246-251.

CHEN Jian, GAO Yu, XU Bin, et al. Sound radiation in stochastic acoustic field using statistical wave superposition method[J]. Journal of Vibration Engineering, 2009,22(3):246-251.

[13]張炳榮,陳劍,陳立濤,等. 二維場預測的快速多極基本解法[J]. 聲學學報, 2014,39(3):347-352.

ZHANG Bing-rong, CHEN Jian, CHEN Li-tao, et al. A fast multipole method of fundamental solutions for two dimensional acoustic radiation problems[J]. Chinese Journal of Acoustics, 2014, 39(3):347-352.

[14]吳紹維,向陽,夏雪寶.基于無單元聲波疊加的自輻射近似解析表達研究[J].振動與沖擊,2014,33(7):79-85.

WU Shao-wei,XIANG Yang,XIA Xue-bao.Approximate analytical expressions of self-radiation terms inclading acoustic pressure and velocity based on element free acoustic wave superposition[J].Journal of Vibration and Shock,2014,33(7): 79-85.

[15]Fairweather G, Karageoghis A, Martin P A. The method of fundamental solutions for scattering and radiation problems [J].Engineering Analysis with Boundary Elements, 2003, 27(7):759-769.

[16]Karageorghis A. A prac tical algorithm for determining the optimal pseudo-boundary in the method of fundamental solutions[J]. Advances in Applied Mathematics and Mechanics, 2009,4(1):510-528.

[17]Seybert A F, Soenarko B. Radiation and scattering of acoustic waves from bodies of arbitrary shape in a three-dimensional half space[J]. Transactions of the ASME, Journal of Vibration, Acoustics, Stress, and Reliability in Design,1988,110:113-117.

[18]Gorzelańczyk P, Kolodziej J A. Some remarks concerning the shape of the source contour with application of the method of fundamental solutions to elastic torsion of prismatic rods[J]. Engineering Analysis with Boundary Elements,2008,32(1):64-75.

[19]Tankelevich R, Fairweather G, Karageorghis A, et al. Potential field based geometric modeling using the method of fundamental solutions[J]. International Journal for Numerical Methods in Engineering, 2006, 68(12):1257-1280.

[20]Tetsuya S,Yosuke Y. Fast multipole boundary element method for large-scale steady-state sound field analysis, part I: setup and validation[J]. Acta Acustica United with Acustica, 2002, 88(4):513-525.