第一作者劉中憲男，副教授，博士后，碩士生導(dǎo)師，1982年生

基于快速多極子基本解方法(FMM-MFS)的彈性波二維散射模擬研究

劉中憲1,2，王冬1,2，梁建文3

(1 天津城建大學(xué)土木工程學(xué)院，天津300384； 2 天津市軟土特性與工程環(huán)境重點實驗室，天津300384;3.天津大學(xué)土木工程系，天津300384)

摘要：針對彈性波二維散射問題，發(fā)展一種新的快速多極子基本解方法(FMM-MFS)。方法基于單層位勢理論，通過在虛邊界上設(shè)置膨脹波線源和剪切波線源以構(gòu)造散射波場，從而避免了奇異性的處理和邊界單元離散；結(jié)合快速多極子展開技術(shù)(FMM)，大幅度降低了計算量和存儲量，突破了傳統(tǒng)方法難以處理大規(guī)模散射問題的瓶頸。以全空間孔洞對P、SV波的二維散射為例，給出了具體求解步驟，并在個人計算機(jī)上實現(xiàn)了上百萬自由度問題的快速精確計算。在方法效率和精度檢驗基礎(chǔ)上，分別以單孔洞和隨機(jī)孔洞群對平面波(P、SV波)的散射為例進(jìn)行計算模擬，揭示了孔洞(群)周圍彈性波散射的若干重要規(guī)律。

關(guān)鍵詞：基本解方法；快速多極子展開方法；快速多極子基本解方法(FMM-MFS)；彈性波散射

基金項目：國家自然科學(xué)

收稿日期：2013-12-30修改稿收到日期：2014-03-12

中圖分類號:Tp12；Tp13.3文獻(xiàn)標(biāo)志碼：A

FMM-MFS soluton to two-dimensinal scattering of elastic waves

LIUZhong-xian1,2,WANGDong1,2,LIANGJian-wen3(1. College of Civil Engineering, Tianjin Chengjian University, Tianjin 300384, China;2. Tianjin Municipal Key Laboratory of Soft Soils and Engineering Environment, Tianjin 300384, China; 3. Department of Civil Engineering, Tianjin University, Tianjin 300072, China)

Abstract:A new algorithm named the fast multipole fundamental solution method(FMM-MFS) was presented for calculating two-dimensional elastic wave scattering problems. The algorithm could avoid the singularity of matrix by placing the line sources of compressional wave and shear wave on a virtual boundary based on the single layer potential theory, and avoid elements discretization on the boundary. Combined with FMM, MFS can solve large-scale problems of wave scattering with greatly reducing computation and the memory requirement. Taking the two-dimensional scattering of P and SV waves around a cavity in elastic full-space as an example, the implement procedures were presented in detail, and up to millions of DOF’s scattering problems were solved successfully on a personal computer. Based on the tests of accuracy and efficiency of FMM-MFS, the scatterings of plane waves around a cavity and randomly distributed group cavities in full-space were solved. Several important conclusions about scattering of elastic waves around cavity(cavities) were obtained.

Key words:method of fundamental solutions (MFS); fast multipole expansion method (FMM); FMM-MFS; scattering of elastic wave

彈性波散射是機(jī)械、材料、地球物理、地震學(xué)等多個領(lǐng)域的一個研究熱點。在眾多彈性波動模擬方法中，基本解方法(MFS)，除了具有邊界元方法降維求解、自動滿足無限遠(yuǎn)輻射條件的優(yōu)勢，同時具有無奇異性、無網(wǎng)格特性以及極高的數(shù)值精度，近些年來在聲波、電磁波以及彈性波波場模擬中獲得了廣泛應(yīng)用[1-2]。該方法最初源于Kupradze等[3-4]的工作，其基本思路是在邊界附近假想一虛擬邊界，并在其上布置虛擬波源來構(gòu)造域內(nèi)散射波場，其本質(zhì)在于利用一系列格林函數(shù)解的線性組合給出偏微分方程的近似解，各基本解未知系數(shù)需首先根據(jù)邊界條件求出。因此MFS也可看做是一種特殊的間接邊界積分方程法，并和Treffez方法具有很大的相似性[5]，非常適合無限域中波動問題求解。

然而在實際工程應(yīng)用中，MFS的一個顯著弱點在于其求解方程矩陣的稠密特性嚴(yán)重影響了高自由度問題(如大尺度、高頻或多體散射)求解的計算效率。針對該問題，近些年來一些學(xué)者嘗試將快速多極子展開(FMM)技術(shù)運用到了MFS求解當(dāng)中。FMM通過將核函數(shù)表達(dá)式場點和源點分離，將原來的源點、場點間單粒相互作用轉(zhuǎn)換為粒組與粒組之間的作用，從而大幅度降低存儲量，并顯著提高計算效率[6-8]。最新發(fā)展的FMM能夠?qū)FS的計算量和存儲量均降低到O(N)量級，使大規(guī)模工程問題快速求解成為可能[9-11]。需注意的是，上述文獻(xiàn)分別針對位勢問題、靜力分析和聲場模擬，對于固體中彈性波散射問題FMM-MFS求解，據(jù)作者所知，目前國內(nèi)外相應(yīng)的研究較少。而需指出的是，相關(guān)的快速多極邊界元方法求解在國內(nèi)外已有一些報道[12-13]。

即針對無限域或半無限域中大規(guī)模彈性波動問題求解(二維平面問題)，基于單層位勢理論，發(fā)展一種新的快速多極子基本解方法。在格林函數(shù)矩陣建立和虛擬波源強(qiáng)度求解中，引入快速多極子展開。該方法整體求解思路比較直觀，且便于數(shù)值實現(xiàn)。以全空間孔洞周圍P、SV波散射為例給出了具體求解步驟，進(jìn)而通過典型算例驗證了方法的計算精度和求解效率。最后以全空間單孔洞對P、SV波的高頻散射和隨機(jī)分布孔洞群對P、SV波的散射為例，驗證方法對實際復(fù)雜問題的適應(yīng)性并得出了一些重要結(jié)論。

1彈性波散射常規(guī)基本解方法(MFS)

圖1　計算模型 Fig.1 Calculation model

為便于說明，首先以全空間圓形孔洞周圍P、SV波散射為例，闡述常規(guī)基本解方法(MFS)的具體實現(xiàn)步驟。計算模型如圖1所示，無限域中一足夠長圓柱形孔洞，假設(shè)平面P或SV波從左側(cè)水平入射(α=0°)，波面平行于孔洞縱軸，待求問題為平面應(yīng)變狀態(tài)下的彈性波散射?；趩螌游粍菰?，可在孔洞內(nèi)部一虛擬源面S上施加虛擬波源以構(gòu)建散射波場。根據(jù)孔洞邊界零應(yīng)力條件構(gòu)建方程以求解波源強(qiáng)度，進(jìn)而計算散射場位移，將其和自由場位移疊加即得到總場位移。根據(jù)試算經(jīng)驗，圖中虛擬波源面S的最優(yōu)半徑低頻時取0.4~0.6倍散射體半徑，高頻時則擴(kuò)大至0.7~0.9倍。該方法對于半無限空間彈性波散射問題也是同樣適用的，僅需在受散射波影響的半空間地表附近布置虛擬波源即可。

1.1散射場構(gòu)造

根據(jù)疊加原理，總位移場和應(yīng)力場可分別表達(dá)為：

u(x)=uf(x)+us(x)

(1a)

σ(x)=σf(x)+σs(x)

(1b)

式中:uf、σf分別表示平面波入射下自由場位移和應(yīng)力(無孔洞時)，us、σs表示散射場位移和應(yīng)力。

由Helmholtz矢量分解原理，二維平面應(yīng)變下位移場可表達(dá)：

u=φ+×(0,0,ψ)

(2)

式中:φ和ψ分別為P波、SV波勢函數(shù), 滿足下列運動方程：

(Δ+h2)φ=0

(3a)

(Δ+k2)ψ=0

(3b)

式中，Δ為二維拉普拉斯算子，h和k分別為P波和SV波波數(shù)。由各向同性假設(shè)，應(yīng)力可表達(dá)為：

σij=λuk,kδi,j+μ(ui,j+uj,i)

(4)

根據(jù)單層位勢理論，散射場可由面S上分布的虛擬波源產(chǎn)生，相應(yīng)位移和應(yīng)力可表達(dá)為：

x∈DE,x1∈S

(5a)

x∈DE, x1∈S

(5b)

(6a)

(6b)

由式(2)和式(4)，全空間中膨脹波、剪切波動力格林函數(shù)可表達(dá)如下：

膨脹波線源:

(7a)

(7b)

γ=k/h

(7c)

γ=k/h

(7d)

(7e)

剪切波線源:

(8a)

(8b)

(8c)

(8d)

(8e)

1.2邊界條件及求解

根據(jù)孔洞邊界上零應(yīng)力條件：

(9a)

(9b)

圓形孔洞情況，邊界法線方向余弦為(cosθ,sinθ),正向和切向應(yīng)力可表達(dá)如下：

σnn=σxxcos2θ+σyysin2θ+2σxysinθcosθ

(10a)

σnt=(-σxx+σyy)sinθcosθ+σxy(cos2θ-sin2θ)

(10b)

根據(jù)孔洞表面自由邊界條件，式(9a)、(9b)可表示為：

(11a)

xn∈B,xm∈S;n=1,…,N

(11b)

式中:cm、dm代表第m個膨脹波和剪切波線源強(qiáng)度，T代表應(yīng)力。通過式(11)計算波源強(qiáng)度，代入式(5a)、(5b)即可求出無限域中任意一點的位移和應(yīng)力。

在上述求解中，采用常規(guī)高斯消元法求解時，系數(shù)矩陣{Tij}為滿陣，解大規(guī)模問題時其效率低下且存儲量大。下面考慮結(jié)合廣義極小余量法(GMRES)[14]和快速多極子展開技術(shù)，突破常規(guī)MFS求解大規(guī)模問題的計算瓶頸。

2彈性波散射FMM-MFS方法

FMM-MFS使用樹結(jié)構(gòu)作為主要的存儲和運算對象，結(jié)合GMRES迭代算法，通過核函數(shù)的多極展開、多極展開傳遞(M2M)、多極展開向局部展開傳遞(M2L)、局部展開傳遞(L2L)以及局部展開5個步驟，避免了對系數(shù)矩陣的直接存儲和對線性方程組的直接求逆。與傳統(tǒng)MFS相比，對于高自由度問題，F(xiàn)MM-MFS大幅度提高了計算效率和經(jīng)濟(jì)性，可將存儲量和計算量降至O(N)量級。下面，首先闡述FMM-MFS的基本原理，進(jìn)而給出彈性波散射FMM-MFS解答過程。

2.1FMM-MFS基本原理

2.1.1模型邊界離散與樹結(jié)構(gòu)生成

圖2　自適應(yīng)樹結(jié)構(gòu) Fig.2 The adaptive tree structure

首先對計算模型邊界進(jìn)行數(shù)值離散(配點)和樹結(jié)構(gòu)生成。對于二維問題，這里使用四叉樹結(jié)構(gòu)。將圖1中彈性體邊界離散為N個邊界點，同時構(gòu)建虛擬荷載面S，同樣離散為N個點，稱之為源點。獲取各離散邊界點和源點的平面坐標(biāo)。注意整個過程并不需要引入單元網(wǎng)格。這與常規(guī)邊界元法一致。然后，將二維彈性體全部邊界放置于一足夠大的正方形中，此正方形代表樹結(jié)構(gòu)的根結(jié)點，稱作0層；將大正方形分成4個小正方形，稱作第1層；依次類推，將正方形所包含的邊界點數(shù)達(dá)到預(yù)設(shè)值為止，刪除不包含任何邊界點的正方形。每個正方形代表一個樹結(jié)點，葉子結(jié)點則不再包含任何子正方形。

圖2中展示了圖1計算模型的四叉樹結(jié)構(gòu)劃分(a)及其示意圖。其中方形四叉樹示意圖(b)為邊界點樹結(jié)構(gòu)，圓形四叉樹示意圖(c)為虛擬源點樹結(jié)構(gòu)。圖2中，源點樹結(jié)構(gòu)主要存儲多極展開系數(shù)、多極展開傳遞系數(shù)；邊界點樹結(jié)構(gòu)主要存儲局部展開傳遞、局部展開系數(shù)；多極展開向局部展開傳遞則由兩部分共同存儲。

2.1.2核函數(shù)的多極展開

下面對式(11)中形如Tij(xn,xm)的核函數(shù)進(jìn)行展開。從源點開始，將源點信息傳遞到對應(yīng)的葉子結(jié)點中心，這一過程稱為多極展開。

參考文獻(xiàn)[15]中展開方式，假設(shè)計算式(11)中核函數(shù)為K(x,y)。在進(jìn)行展開之前需要指出：本文中多極展開在源點Y處進(jìn)行，而局部展開則在邊界點X處進(jìn)行。

圖3　快速多極展開相關(guān)點 Fig.3 The related points for FMM

K(x,y)=N(x,y0)M(y0,y)

(12)

式中:M(y0,y)稱為展開點y0(源點葉子結(jié)點)的多極展開系數(shù)。這里，M(y0,y)只與y0變化有關(guān)，故只需計算一遍即可和滿足條件的任意點x進(jìn)行計算。

2.1.3多極展開的傳遞(M2M)

如圖3，當(dāng)展開點從y0移到較近點y1(源點父結(jié)點)時，可得到新的展開系數(shù)，即：

M(y1,y)=I(y1,y0)M(y0,y)

(13)

由式(11)可以看出，我們不需要重新計算展開系數(shù)，僅需將已計算的多極展開系數(shù)平移即可得到新的展開系數(shù)。其中，I(y1,y0)為多極展開傳遞系數(shù)。

2.1.4多極展開向局部展開的傳遞(M2L)

N(x,y0)=L(x,x0)H(x0,y0)

(14)

式中:L(x,x0)即為局部展開系數(shù)。與此同時，式(14)同樣滿足多極展開向局部展開的傳遞關(guān)系(M2L)，其中H(x0,y0)即為傳遞系數(shù)。

在進(jìn)行局部展開之前，需要了解“鄰居”和“相互作用列表”兩個概念。結(jié)點M的“鄰居”是指與M位于同一層且至少有一個角點共用，一個結(jié)點至多有8個“鄰居”。結(jié)點M的“相互作用列表”是指某結(jié)點的父結(jié)點與M的父結(jié)點互為“鄰居”且本身與結(jié)點M不為“鄰居”，一個結(jié)點至多有27個“相互作用列表”，如圖4所示。

圖4　結(jié)點“鄰居”與“相互作用列表” Fig.4 Adjacent cells and cells in the interaction list

2.1.5局部展開的傳遞(L2L)

如圖4，當(dāng)局部展開點從x0移到較近點x1(邊界父結(jié)點)時，可得到新的局部展開系數(shù)，即：

L(x,x1)=L(x,x0)J(x0,x1)

(15)

由式(15)可以看出，新的局部展開系數(shù)不需要重新計算，僅需將已計算的局部展開系數(shù)平移即可得到新的局部展開系數(shù)。其中，J(x0,x1)即為局部展開傳遞系數(shù)。

2.1.6FMM-MFS的實施

依據(jù)快速多極子基本解方法的基本原理，將其具體實施過程表述如下：

(1)將彈性體邊界離散化，這與常規(guī)MFS方法一致，進(jìn)而自動生成自適應(yīng)樹結(jié)構(gòu)。

(2)從源點y開始，將源點信息傳給葉子結(jié)點，通過式(12)在葉子中心展開形成多極展開系數(shù)。

(3)根據(jù)需要，可將多級展開點平移至對應(yīng)父結(jié)點(M2M)，計算式(13)構(gòu)成傳遞關(guān)系。在這一步中可多次向上層父細(xì)胞傳遞，在滿足精度的條件下盡可能多的往上傳遞。

(4)將源點樹結(jié)構(gòu)各層結(jié)點的多極展開系數(shù)傳遞給邊界點樹結(jié)構(gòu)同層相應(yīng)結(jié)點，形成局部展開系數(shù)(M2L)，式(14)構(gòu)建此關(guān)系。所謂“相應(yīng)”，即邊界父結(jié)點接收非“鄰居”的多極展開系數(shù)，邊界子結(jié)點(非最底層父結(jié)點)接收“相互作用列表”的多極展開系數(shù)。

(5)式(15)將邊界點樹結(jié)構(gòu)中各層結(jié)點的局部展開系 數(shù)層層傳遞直到葉子結(jié)點(L2L)。至此，將葉子結(jié)點處累加而來的所有局部展開系數(shù)通過局部展開傳給邊界點x，完成整個展開以及傳遞過程。剩余近場源點，即葉子結(jié)點的鄰居以及自身所包含的源點，直接采用常規(guī)MFS計算。

經(jīng)過上述5個步驟，核函數(shù)K(x,y)最終展開為：

K(x,y)=

L(x,x0)J(x0,x1)H(x1,y1)I(y1,y0)M(y0,y)

(16)

式中:H(x1,y1)為M2L系數(shù)，其余函數(shù)均已說明。

計算過程中，采用GMRES進(jìn)行迭代求解運算。為降低迭代次數(shù)，在GMRES求解前采用預(yù)處理技術(shù)[16]，這里可采用近場格林函數(shù)解構(gòu)造窄帶稀疏矩陣，對方程組進(jìn)行預(yù)處理。

2.2彈性波二維散射FMM-MFS解答

根據(jù)式(6a)和(6b)全空間勢函數(shù)形式，采用Graf加法定理[17]：

(17)

式中:Z可以是J、Y、H(1)、H(2)以及他們的線性組合，n為截斷數(shù)，其中幾何參數(shù)如圖5所示。

圖5　符號定義 Fig.5 Symbol definition

(18)

(20)

Jm(kρ′)eimγ′=

(21)

(22)

將上式分別代入式(7)即可求出位移和應(yīng)力。同理，可導(dǎo)出式(6a)展開式，僅需將式(22)中的剪切波波數(shù)k改為縱波波數(shù)h即可，此處不再贅述。

上述即為全空間中彈性波平面內(nèi)二維散射的快速多極基本解方法解答。對于半空間情況，可利用相應(yīng)的半空間格林函數(shù)進(jìn)行展開求解，或者仍采用全空間基本解，但需對受散射波影響的半空間表面進(jìn)行離散。

3數(shù)值算例與分析

3.1FMM-MFS計算精度及效率檢驗

首先以全空間孔洞為例，檢驗快速多極基本解方法的精度、迭代收斂速度及數(shù)值穩(wěn)定性。這里引入無量綱頻率η的概念。它定義為散射體等效直徑與入射波波長之比，即：η=2a/λ=ωa/πβ。分別利用快速多極子基本解方法和解析解即波函數(shù)展開法計算全空間圓形孔洞對平面P、SV波的散射。其孔洞直徑取10 m，入射波頻率取為100 Hz，介質(zhì)剪切波速取1 000 m/s，泊松比取為0.25，換算無量綱頻率η=1，圖6為全空間孔洞表面位移幅值結(jié)果對比情況。邊界配點數(shù)取N=1 000，即自由度數(shù)為2 000。容易看出，本文所發(fā)展快速多極子基本解方法同解析解結(jié)果吻合良好。

圖7為傳統(tǒng)基本解方法與快速多極子基本解方法的迭代收斂殘差與迭代步關(guān)系曲線，展開截斷階數(shù)取為20。從圖中可以看出，兩種方法迭代收斂速度相當(dāng)； FMM-MFS方法下，自由度數(shù)不同其收斂速度也相當(dāng)，從而證明了此方法的良好收斂性。圖8給出了P波入

圖6　全空間孔洞表面位移幅值FMM-MFS和解析解 Fig.6 Comparisons between the displacement amplitude around the cavity by FMM-MFS and that of analytic solution

圖7　迭代收斂曲線 Fig.7 The convergence of the iteration curve

圖8　FMM-MFS與常規(guī)MFS 的總CPU時間結(jié)果比較 Fig.8 Comparisons between the total CPU time of FMM-MFS and that of MFS

射時計算時間隨自由度之間的關(guān)系，入射波頻率取η=5.0，入射角度取α=0，即水平入射。其中父細(xì)胞劃分到7層，級數(shù)累加從-6到6；葉子細(xì)胞劃分到11層，級數(shù)累加從-1到1。從圖中可清晰看出，自由度超過5 000以后，常規(guī)算法計算時間急劇增加。10 000自由度時計算時間已經(jīng)達(dá)到約25 000 s，遠(yuǎn)遠(yuǎn)大于FMM-MFS計算160萬自由度計算時間。而且，常規(guī)方法的內(nèi)存需要量大，目前一般個人臺式電腦也只能算至大約兩萬自由度。然而，采用FMM-MFS方法，160萬自由度問題僅需約9 000 s，證明了FMM-MFS的計算效率遠(yuǎn)遠(yuǎn)高于傳統(tǒng)MFS，且突破了傳統(tǒng)MFS不能計算大規(guī)模問題的弱點。

為檢驗本文方法數(shù)值穩(wěn)定性，表1給出了不同邊界配點數(shù)下FMM-MFS位移幅值計算結(jié)果，配點數(shù)分別取2 000、20 000和100 000。由表可知，隨配點數(shù)增大，孔洞表面位移幅值相差僅在小數(shù)點后7位，位移結(jié)果收斂良好，反映了該方法具有良好的數(shù)值穩(wěn)定性。

以上的計算結(jié)果以及計算時間統(tǒng)計，均是在8G內(nèi)存、32位XP操作系統(tǒng)的個人電腦上(Dell Pretision T5500)，使用Matlab進(jìn)行編程運算得到。

表1　位移幅值數(shù)值穩(wěn)定性檢驗 (η=1.0)

3.2全空間單孔洞對平面P、SV波的高頻散射

圖9給出了高頻P、SV波入射下，全空間單孔洞面上的位移結(jié)果。圖10分別給出了高頻P、SV波入射下，全空間孔洞附近15倍孔洞半徑范圍內(nèi)的x、y方向位移幅值云圖?？锥粗睆饺?00 m，入射波頻率取為500 Hz，介質(zhì)剪切波速取1 000 m/s，泊松比取為0.25，換算為無量綱頻率η=50。入射角度取0°，即左側(cè)水平入射。

為充分提高計算精度和計算效率，孔洞表面離散5 000點(自由度10 000)。分別采用常規(guī)MFS和FMM-MFS計算，計算機(jī)耗時分別約為6 000 s和410 s，F(xiàn)MM-MFS計算時間不同于η=5.0是因為劃分層數(shù)不同所導(dǎo)致。結(jié)果表明本文方法能夠高效精確地求解高頻波散射問題。另外容易看出，高頻P、SV波水平入射，迎波面一側(cè)出現(xiàn)顯著的位移放大效應(yīng)(集中在一倍孔洞半徑范圍內(nèi))，位移幅值達(dá)到入射波幅值的2倍多，在孔洞右側(cè)則呈現(xiàn)很長的“陰影區(qū)”，表明孔洞對高頻P、SV波有很強(qiáng)的“屏障”效應(yīng)。因此，地下結(jié)構(gòu)抗爆設(shè)計尤需注意迎爆面的動應(yīng)力集中效應(yīng)。

圖9　P、SV波入射下孔洞表面位移幅值 (η=50) Fig.9 Displacement amplitude on the cavity surface for P 、SV waves incident(η=50)

圖10　P、SV波入射下孔洞附近位移幅值云圖 Fig.10 Displacement amplitude in the neighborhood of the cavity for P 、SV waves incident

3.3全空間孔洞群對平面P、SV波的散射

為考察FMM-MFS方法對復(fù)雜的多體散射問題求解的適應(yīng)性，圖11、12分別給出隨機(jī)分布的孔洞群對P、SV波的散射結(jié)果。選取與圖10所示相同范圍內(nèi)隨機(jī)分布的30個圓形孔洞為計算模型，水平左側(cè)入射，孔洞直徑取為100 m，入射頻率分別取1、5、10、20 Hz，介質(zhì)剪切波速1 000 m/s，泊松比0.25，換算無量綱頻率η=0.1、0.5、1.0、2.0。分別給出了不同頻率P、SV波水平入射下孔洞附近x、y方向位移幅值云圖。

圖11　P波入射下隨機(jī)多孔洞周圍位移云圖 Fig.11 Displacement amplitude in the neighborhood of random group cavities for P waves incident

圖12　SV波入射下隨機(jī)孔洞群周圍位移云圖 Fig.12 Displacement amplitude in the neighborhood of random group cavities for SV waves incident

結(jié)果表明：P、SV波左側(cè)水平入射下，受孔洞群之間散射波相互干涉影響，孔洞群對P、SV波的散射與單孔洞情況具有明顯差別，孔邊動力放大效應(yīng)更為顯著，即便是在低頻η=0.1情況，位移放大系數(shù)接近2倍(單洞時則表現(xiàn)為弱散射)。而在孔洞群的右側(cè)，位移幅值很小，可以看出孔洞群對入射P、SV波的“屏障”效應(yīng)更為明顯。隨著入射波頻率增加，孔洞間散射波的相干效應(yīng)更為強(qiáng)烈。如η=2.0時，迎波面一側(cè)前兩列孔洞邊緣出現(xiàn)顯著的位移放大效應(yīng)，P、SV波入射下放大系數(shù)峰值分別達(dá)到4.5和3.2。值得指出的是，對于η=2.0情況，模型計算自由度約為6 000，計算時間僅為100 s，即實現(xiàn)了孔洞群對波散射高效精確求解。

4結(jié)論

針對彈性波二維散射問題，基于單層位勢理論，結(jié)合快速多極子展開技術(shù)，發(fā)展了一種新的快速多極子基本解方法。數(shù)值分析表明：該方法具有高精度、良好的收斂性及數(shù)值穩(wěn)定性、無網(wǎng)格特性，在提高運算速度的同時能夠大幅度降低計算存儲量。通過算法優(yōu)化，可在目前主流計算機(jī)上實現(xiàn)數(shù)百萬自由度彈性波散射問題的快速精確求解。最后以全空間單孔洞和隨機(jī)分布孔洞群對P、SV波散射為例，研究了孔洞周圍高頻波散射和多體散射的若干規(guī)律，可為孔隙介質(zhì)波動分析、材料無損檢測等提供部分理論依據(jù)。另外，本文方法對于大規(guī)模聲波、電磁波散射分析等同樣具有參考價值。

參考文獻(xiàn)

[1]Golberg M A, Chen C S. The method of fundamental solutions for potential, Helmholtz and diffusion problems[J]. Boundary Integral Methods-Numerical and Mathematical Aspects, 1998: 103-176.

[2]Fairweather G, Karageorghis 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.

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

[4]Mathon R,Johnston R L. The approximate solution of elliptic boundary-value problems by fundamental solutions[J]. SIAM Journal on Numerical Analysis, 1977, 14(4): 638-650.

[5]Chen J T, Lee Y T, Yu S R, et al. Equivalence between the Trefftz method and the method of fundamental solution for the annular Green’s function using the addition theorem and image concept[J]. Engineering Analysis with Boundary Elements, 2009, 33(5): 678-688.

[6]Greengard L, Rokhlin V. A fast algorithm for particle simulations[J]. Journal of Computational Physics, 1997, 135(2): 280-292.

[7]姚振漢，王海濤. 邊界元法[M].北京：高等教育出版社，2009.

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

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

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

[10]許強(qiáng), 蔣彥濤, 張志佳. 快速多極虛邊界元法對含圓孔薄板有效彈性模量的模擬分析[J]. 計算力學(xué)學(xué)報, 2010, 27(3): 548-555.

XU Qiang,JIANG Yan-tao,ZHANG Zhi-jia. Fast multipole VBEM for analyzing the effective elastic moduli of a sheet containing circular holes[J].Chinese Journal of Computational Mechanics, 2010, 27(3): 548-555.

[11]Jiang X R, 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),1341008.

[12]Chen Y H, Chew W C, Zeroug S. Fast multipole method as an efficient solver for 2D elastic wave surface integral equations[J]. Computational Mechanics, 1997, 20(6): 495-506.

[13]Chaillat S, Bonnet M, Semblat J F. A new fast multi-domain BEM to model seismic wave propagation and amplification in 3-D geological structures[J]. Geophysical Journal International, 2009, 177(2): 509-531.

[14]Saad Y, Schultz M H. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[J]. SIAM Journal on Scientific and Statistical Computing, 1986, 7(3): 856-869.

[15]Yoshida K. Applications of fast multipole method to boundary integral equation method[D]. Dept. of Global EnvironmentEng., Kyoto Univ., Japan, 2001.

[16]王海濤. 快速多極邊界元法在二維彈性力學(xué)中的應(yīng)用[D]. 北京：清華大學(xué), 2002.

[17]Utsunomiya T, Watanabe E, Nishimura N. Fast multipole algorithm for wave diffraction/ radiation problems and its application to VLFS in variable water depth and topography[C]//Proc 20th Int Conf on Offshore Mech and Arcctic Engrg,2001.