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

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

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

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

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

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

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

基金項目：國家自然科學

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

中圖分類號:Tp12；Tp13.3文獻標志碼：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

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

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

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

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

圖1　計算模型 Fig.1 Calculation model

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

1.1散射場構造

根據疊加原理，總位移場和應力場可分別表達為：

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

(1a)

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

(1b)

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

由Helmholtz矢量分解原理，二維平面應變下位移場可表達：

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

(2)

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

(Δ+h2)φ=0

(3a)

(Δ+k2)ψ=0

(3b)

式中，Δ為二維拉普拉斯算子，h和k分別為P波和SV波波數。由各向同性假設，應力可表達為：

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

(4)

根據單層位勢理論，散射場可由面S上分布的虛擬波源產生，相應位移和應力可表達為：

x∈DE,x1∈S

(5a)

x∈DE, x1∈S

(5b)

(6a)

(6b)

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

膨脹波線源:

(7a)

(7b)

γ=k/h

(7c)

γ=k/h

(7d)

(7e)

剪切波線源:

(8a)

(8b)

(8c)

(8d)

(8e)

1.2邊界條件及求解

根據孔洞邊界上零應力條件：

(9a)

(9b)

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

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

(10a)

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

(10b)

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

(11a)

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

(11b)

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

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

2彈性波散射FMM-MFS方法

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

2.1FMM-MFS基本原理

2.1.1模型邊界離散與樹結構生成

圖2　自適應樹結構 Fig.2 The adaptive tree structure

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

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

2.1.2核函數的多極展開

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

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

圖3　快速多極展開相關點 Fig.3 The related points for FMM

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

(12)

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

2.1.3多極展開的傳遞(M2M)

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

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

(13)

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

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

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

(14)

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

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

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

2.1.5局部展開的傳遞(L2L)

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

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

(15)

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

2.1.6FMM-MFS的實施

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

(1)將彈性體邊界離散化，這與常規MFS方法一致，進而自動生成自適應樹結構。

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

(3)根據需要，可將多級展開點平移至對應父結點(M2M)，計算式(13)構成傳遞關系。在這一步中可多次向上層父細胞傳遞，在滿足精度的條件下盡可能多的往上傳遞。

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

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

經過上述5個步驟，核函數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系數，其余函數均已說明。

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

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

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

(17)

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

圖5　符號定義 Fig.5 Symbol definition

(18)

(20)

Jm(kρ′)eimγ′=

(21)

(22)

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

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

3數值算例與分析

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

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

圖7為傳統基本解方法與快速多極子基本解方法的迭代收斂殘差與迭代步關系曲線，展開截斷階數取為20。從圖中可以看出，兩種方法迭代收斂速度相當； FMM-MFS方法下，自由度數不同其收斂速度也相當，從而證明了此方法的良好收斂性。圖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與常規MFS 的總CPU時間結果比較 Fig.8 Comparisons between the total CPU time of FMM-MFS and that of MFS

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

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

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

表1　位移幅值數值穩定性檢驗 (η=1.0)

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

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

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

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

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

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

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

4結論

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

參考文獻

[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]崔曉兵, 季振林. 快速多極子邊界元法在吸聲材料聲場計算中的應用[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]許強, 蔣彥濤, 張志佳. 快速多極虛邊界元法對含圓孔薄板有效彈性模量的模擬分析[J]. 計算力學學報, 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]王海濤. 快速多極邊界元法在二維彈性力學中的應用[D]. 北京：清華大學, 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.