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

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

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

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

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

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

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

基金項(xiàng)目：國(guó)家自然科學(xué)

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

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

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

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

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

圖1　計(jì)算模型 Fig.1 Calculation model

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

1.1散射場(chǎng)構(gòu)造

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

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

(1a)

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

(1b)

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

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

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

(2)

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

(Δ+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ù)單層位勢(shì)理論，散射場(chǎng)可由面S上分布的虛擬波源產(chǎn)生，相應(yīng)位移和應(yīng)力可表達(dá)為：

x∈DE,x1∈S

(5a)

x∈DE, x1∈S

(5b)

(6a)

(6b)

由式(2)和式(4)，全空間中膨脹波、剪切波動(dòng)力格林函數(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個(gè)膨脹波和剪切波線源強(qiáng)度，T代表應(yīng)力。通過(guò)式(11)計(jì)算波源強(qiáng)度，代入式(5a)、(5b)即可求出無(wú)限域中任意一點(diǎn)的位移和應(yīng)力。

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

2彈性波散射FMM-MFS方法

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

2.1FMM-MFS基本原理

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

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

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

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

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

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

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

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

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

(12)

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

2.1.3多極展開(kāi)的傳遞(M2M)

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

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

(13)

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

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

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

(14)

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

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

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

2.1.5局部展開(kāi)的傳遞(L2L)

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

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

(15)

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

2.1.6FMM-MFS的實(shí)施

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

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

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

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

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

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

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

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ù)均已說(shuō)明。

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

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

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

(17)

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

圖5　符號(hào)定義 Fig.5 Symbol definition

(18)

(20)

Jm(kρ′)eimγ′=

(21)

(22)

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

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

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

3.1FMM-MFS計(jì)算精度及效率檢驗(yàn)

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

圖7為傳統(tǒng)基本解方法與快速多極子基本解方法的迭代收斂殘差與迭代步關(guān)系曲線，展開(kāi)截?cái)嚯A數(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時(shí)間結(jié)果比較 Fig.8 Comparisons between the total CPU time of FMM-MFS and that of MFS

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

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

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

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

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

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

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

為考察FMM-MFS方法對(duì)復(fù)雜的多體散射問(wèn)題求解的適應(yīng)性，圖11、12分別給出隨機(jī)分布的孔洞群對(duì)P、SV波的散射結(jié)果。選取與圖10所示相同范圍內(nèi)隨機(jī)分布的30個(gè)圓形孔洞為計(jì)算模型，水平左側(cè)入射，孔洞直徑取為100 m，入射頻率分別取1、5、10、20 Hz，介質(zhì)剪切波速1 000 m/s，泊松比0.25，換算無(wú)量綱頻率η=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è)水平入射下，受孔洞群之間散射波相互干涉影響，孔洞群對(duì)P、SV波的散射與單孔洞情況具有明顯差別，孔邊動(dòng)力放大效應(yīng)更為顯著，即便是在低頻η=0.1情況，位移放大系數(shù)接近2倍(單洞時(shí)則表現(xiàn)為弱散射)。而在孔洞群的右側(cè)，位移幅值很小，可以看出孔洞群對(duì)入射P、SV波的“屏障”效應(yīng)更為明顯。隨著入射波頻率增加，孔洞間散射波的相干效應(yīng)更為強(qiáng)烈。如η=2.0時(shí)，迎波面一側(cè)前兩列孔洞邊緣出現(xiàn)顯著的位移放大效應(yīng)，P、SV波入射下放大系數(shù)峰值分別達(dá)到4.5和3.2。值得指出的是，對(duì)于η=2.0情況，模型計(jì)算自由度約為6 000，計(jì)算時(shí)間僅為100 s，即實(shí)現(xiàn)了孔洞群對(duì)波散射高效精確求解。

4結(jié)論

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

參考文獻(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]崔曉兵, 季振林. 快速多極子邊界元法在吸聲材料聲場(chǎng)計(jì)算中的應(yīng)用[J]. 振動(dòng)與沖擊, 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), 蔣彥濤, 張志佳. 快速多極虛邊界元法對(duì)含圓孔薄板有效彈性模量的模擬分析[J]. 計(jì)算力學(xué)學(xué)報(bào), 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.