劉 青，李世俊，武 偉,?

(1．西北工業(yè)大學(xué)航天學(xué)院，陜西西安 710072；2．寶雞市158信箱，陜西寶雞 721000)

邊界元法是近年來(lái)受到廣泛關(guān)注的一種數(shù)值方法[1]．它將控制方程轉(zhuǎn)化為邊界積分方程，通過(guò)求解邊界上的未知量獲得實(shí)際問(wèn)題的解．它將問(wèn)題降低一維，形成的方程自由度小，計(jì)算量降低，處理方便，求解精度高．因此，近年來(lái)邊界元法在許多領(lǐng)域得到了廣泛應(yīng)用[2]．然而，常規(guī)邊界元法得到的系數(shù)矩陣為滿陣，計(jì)算量隨自由度數(shù)目的增加而呈平方量級(jí)增長(zhǎng)，因此，難以應(yīng)用于大規(guī)模工作問(wèn)題中．為解決該問(wèn)題，近三十年來(lái)，人們提出了多種快速算法，如快速多極方法（Fast Multipole Method, FMM）[3-5]、H-矩陣法[5]、快速小波方法[6-7]、PFFT（Precorrected Fast Fourier Transform）[8]、ACA（Adaptive Cross Approximation）[9]等．這些算法可以應(yīng)用于加速邊界元法的求解，形成的快速邊界元法可以將邊界元法的計(jì)算復(fù)雜度由O(N2)降低到O(NlogαN)α≥0[10-14]．

快速多極方法因其優(yōu)點(diǎn)和影響力，已經(jīng)廣泛應(yīng)用于靜電學(xué)、聲學(xué)、電磁學(xué)、彈性動(dòng)力學(xué)等許多領(lǐng)域[10-13,15]．它的基本思路是，首先將邊界逐次剖分成越來(lái)越小的區(qū)域，然后對(duì)每個(gè)小區(qū)域，根據(jù)核函數(shù)的遠(yuǎn)場(chǎng)低階解析展開，依次采用S2M（Source to Multipole）、M2M（Multipole to Multipole）、M2L（Multipole to Local）、L2L（Local to Local）、L2T（Local to Target）變換，計(jì)算它與遠(yuǎn)場(chǎng)區(qū)域之間的相互作用，最后實(shí)現(xiàn)統(tǒng)一．其優(yōu)點(diǎn)是，可以達(dá)到任意指定的精度，而且計(jì)算效率高，可以將邊界元法的計(jì)算量降低到O(N)．但是，缺點(diǎn)是它的各步變換公式是由核函數(shù)的解析展開式推導(dǎo)而來(lái)，而不同的展開方式對(duì)計(jì)算效率有較大的影響，且不同物理問(wèn)題的核函數(shù)往往不同，解析展開式也不相同．因此，對(duì)于每一個(gè)新型物理問(wèn)題的核函數(shù)，都需要重新推導(dǎo)它的低階解析展開式以及各步變換的計(jì)算公式，以構(gòu)造FMM快速計(jì)算格式，這為FMM向不同領(lǐng)域的推廣帶來(lái)了很大的不便．

為解決上述基于核函數(shù)解析展開的FMM（下面簡(jiǎn)稱“解析FMM”）的這一缺陷，本文將以文獻(xiàn)[16]提出的核無(wú)關(guān)快速多極方法（Kernel Independent Fast Multipole Method, KIFMM）為例，研究探討KIFMM與解析FMM之間的關(guān)系，進(jìn)一步深入研究KIFMM中等效源近似與解析FMM中核函數(shù)解析展開的關(guān)系，從而將KIFMM統(tǒng)一到傳統(tǒng)解析FMM的理論框架下，并在此基礎(chǔ)上，開展應(yīng)力計(jì)算實(shí)例，驗(yàn)證本文理論的正確性．

1 解析FMM和KIFMM原理

基于等效源的KIFMM與解析FMM方法的結(jié)構(gòu)相同，它們都需要構(gòu)建八叉樹（對(duì)于二維問(wèn)題，構(gòu)建四叉樹）結(jié)構(gòu)，在八叉樹上進(jìn)行上行和下行運(yùn)算．首先，定義第0層的正方體格子C0，使所有源點(diǎn)和目標(biāo)點(diǎn)均在C0內(nèi)，然后，將C0均分成8個(gè)子格子，這樣就得到了八叉樹的下一層，即第1層的格子，C0即為這些格子的父格子．一直細(xì)分下去，直到八叉樹中每個(gè)葉子格子中源點(diǎn)個(gè)數(shù)和目標(biāo)點(diǎn)個(gè)數(shù)均不超過(guò)預(yù)先給定的值s．對(duì)于每一個(gè)格子C，定義它的近場(chǎng)NC為同層上與它相鄰格子以及它自身所在的區(qū)域，其余區(qū)域稱為它的遠(yuǎn)場(chǎng)FC．對(duì)格子與其遠(yuǎn)場(chǎng)之間的相互作用KIFMM和FMM都是采用上行、下行變換進(jìn)行快速近似計(jì)算，但KIFMM與FMM最顯著的不同之處在于，它采用上行等效源代替解析FMM中在源點(diǎn)附近的多極展開，采用下行檢測(cè)勢(shì)代替解析FMM中目標(biāo)點(diǎn)附近的局部展開．本節(jié)首先分析解析FMM和KIFMM的加速計(jì)算原理．

1.1 解析FMM

解析FMM的遠(yuǎn)場(chǎng)計(jì)算步驟如圖1所示．

圖1 解析 FMM 中的遠(yuǎn)場(chǎng)變換計(jì)算

假設(shè)需要計(jì)算某格子內(nèi)N個(gè)源點(diǎn){yj}點(diǎn)處的源強(qiáng)度{qj}對(duì)它某個(gè)遠(yuǎn)場(chǎng)格子內(nèi)M個(gè)目標(biāo)點(diǎn){xi}點(diǎn)處位勢(shì)的貢獻(xiàn){pi}，

其中G(x,y)為核函數(shù)．FMM首先將核函數(shù)分別在源點(diǎn)附近某點(diǎn)0y處進(jìn)行低階多極展開，在目標(biāo)點(diǎn)附近某點(diǎn)0x處進(jìn)行低階局部展開，可得：

其中，展開項(xiàng)數(shù)T()ε與選用的精度ε有關(guān)．這樣，式（1）可以通過(guò)如下3步計(jì)算：

1）根據(jù)源點(diǎn)的源強(qiáng)度qj，計(jì)算y0處的多極展開系數(shù)ms(y0)．

即FMM中的 S2M 變換．

2）根據(jù)y0處的多極展開系數(shù)ms(y0)計(jì)算x0處的局部展開系數(shù)lt(x0)．

即FMM中的M2L變換．

3）根據(jù)x0處的局部展開系數(shù)lt(x0)計(jì)算目標(biāo)點(diǎn)處的位勢(shì)pi．

比較兩種計(jì)算方法可知，采用式（1）直接計(jì)算{pi}，計(jì)算量為O(MN)，而若采用單層FMM計(jì)算，則計(jì)算量為O(T(ε)?(M+N+T(ε))．當(dāng)M、N遠(yuǎn)大于T(ε)時(shí)，采用單層FMM計(jì)算即可顯著減少計(jì)算量．對(duì)于邊界元法來(lái)說(shuō)，單層FMM最多可將O(N2)的計(jì)算量降低到O(N3/2)．對(duì)于大規(guī)模問(wèn)題，為獲得O(N)的計(jì)算復(fù)雜度，則需要采用多層FMM算法，即引入M2M變換和L2L變換，其中，M2M變換將低層格子上的多極展開系數(shù)變換為高層格子上的多極展開系數(shù)，L2L變換將高層格子上的局部展開系數(shù)變換為低層格子上的局部展開系數(shù)，它們可以采用插值的方法計(jì)算[17]．多層FMM的計(jì)算步驟如圖1所示．從圖1中也可以看出，上行運(yùn)算的S2M、M2M和下行運(yùn)算的L2T、L2L基本上是對(duì)稱的．

1.2 KIFMM

KIFMM的算法結(jié)構(gòu)與解析FMM類似，但遠(yuǎn)場(chǎng)作用采用等效源和檢測(cè)勢(shì)來(lái)近似計(jì)算．單層KIFMM計(jì)算過(guò)程如圖2所示．

圖2 單層KIFMM計(jì)算步驟（實(shí)心點(diǎn)表示等效點(diǎn)，空心圓環(huán)表示檢測(cè)點(diǎn)）

其中，T()ε為上行等效點(diǎn)的個(gè)數(shù)，T′()ε為上行檢測(cè)點(diǎn)的個(gè)數(shù)，它們均與選用的精度有關(guān)，寫為矩陣形式，得：

對(duì)應(yīng)于圖2中的第1步、第2步。

對(duì)應(yīng)于圖2中的第5步．

為保持KIFMM中上行運(yùn)算和下行運(yùn)算類似于解析FMM的對(duì)稱性，可以將第1步和第2步稱為S2M變換，它將格子C中的源變換為上行等效源，變換矩陣為：

矩陣維數(shù)為M×T(ε)．所以，采用這種方法計(jì)算量為O(T(ε)?(M+N+T(ε))，與解析FMM相同．

對(duì)于大規(guī)模問(wèn)題，KIFMM 同樣可以通過(guò)引入M2M變換和L2L變換，采用多層算法格式以獲得O(N)的計(jì)算復(fù)雜度．其中，M2M變換將低層格子上的上行等效源變換為高層格子上的上行等效源，計(jì)算公式與S2M類似，只需要將公式（6）中源點(diǎn)上的源替換為低層格子上的上行等效源即可．L2L變換將高層格子上的下行檢測(cè)勢(shì)變換為低層格子上的下行檢測(cè)勢(shì)，計(jì)算公式與L2T類似，只需要將公式（10）中目標(biāo)點(diǎn)上的位勢(shì)替換為低層格子上的下行檢測(cè)勢(shì)即可．多層KIFMM的遠(yuǎn)場(chǎng)計(jì)算變換步驟如圖3所示．

圖3 KIFMM的遠(yuǎn)場(chǎng)變換計(jì)算

2 KIFMM中等效源近似與核函數(shù)展開

本節(jié)研究KIFMM中等效源近似與解析FMM中核函數(shù)解析展開的關(guān)系，討論等效源近似對(duì)應(yīng)的核函數(shù)展開形式，將兩者統(tǒng)一到核函數(shù)展開理論的框架下．

首先觀察利用上行等效源計(jì)算C內(nèi)的源在其遠(yuǎn)場(chǎng)產(chǎn)生的位勢(shì)場(chǎng)，即圖2中第1步和第2步．C內(nèi)的任意一點(diǎn)y處單位源在遠(yuǎn)場(chǎng)任意點(diǎn)x處產(chǎn)生的位勢(shì)場(chǎng)p(x)為：

根據(jù)KIFMM，采用上行等效源計(jì)算遠(yuǎn)場(chǎng)處的位勢(shì)場(chǎng)p(x)，需要進(jìn)行如下3步計(jì)算．

所以，構(gòu)造上行等效源，實(shí)際上相當(dāng)于采用函數(shù)（20）對(duì)核函數(shù)在上行等效點(diǎn)處進(jìn)行T()ε項(xiàng)近似低階展開．

類似地，采用格子D的下行等效源，計(jì)算它遠(yuǎn)場(chǎng)DF中任一點(diǎn)y處單位源在D內(nèi)任一點(diǎn)x點(diǎn)處產(chǎn)生的位勢(shì)，可得：

所以，下行計(jì)算中，實(shí)際上相當(dāng)于采用函數(shù)（23）對(duì)核函數(shù)在下行等效點(diǎn)處進(jìn)行T()ε項(xiàng)近似低階展開．

綜合式（21）和（24），C中任一源點(diǎn)y上單位源在它遠(yuǎn)場(chǎng)格子D中任一目標(biāo)點(diǎn)x上產(chǎn)生的位勢(shì)可表示為：

對(duì)比式（25）和解析FMM中核函數(shù)的解析展開式（2）可知，KIFMM的等效源近似相當(dāng)于采用（25）式對(duì)核函數(shù)進(jìn)行近似低階展開．

3 算 例

本節(jié)以兩個(gè)彈性場(chǎng)的計(jì)算為例來(lái)驗(yàn)證本文理論和方法的正確性．

由于彈性動(dòng)力學(xué)核無(wú)關(guān)快速多極邊界元法（Kernel Independent Fast Multipole Boundary Element Method, KIFMBEM）的性能已經(jīng)在文獻(xiàn)[18]中進(jìn)行過(guò)詳細(xì)討論，因此本文對(duì)此不再重復(fù)，僅研究KIFMBEM計(jì)算應(yīng)力應(yīng)變場(chǎng)的精度．程序采用C++串行編程實(shí)現(xiàn)．Nystr?m方法離散得到的方程用GMRes求解器求解，采用對(duì)角塊預(yù)處理技術(shù)以加速收斂．方程組迭代收斂限與選用的算法精度ε相同．

首先采用一個(gè)有解析解的算例來(lái)驗(yàn)證本文理論和方法的正確性．假設(shè)在無(wú)限大彈性體中存在一個(gè)球形空腔，空腔半徑為1 m，內(nèi)壁面受幅值為100 MPa、頻率為2 672 Hz的簡(jiǎn)諧壓強(qiáng)作用，如圖4所示．彈性體彈性模量為200 GPa，密度為7 850 kg/m3，泊松比為0.31．計(jì)算彈性體內(nèi)的位移場(chǎng)、應(yīng)變場(chǎng)和應(yīng)力場(chǎng)．

圖4 受內(nèi)壓的球形空腔

首先在球面上劃分網(wǎng)格，求解位移邊界積分方程得到邊界上的位移和面力．創(chuàng)建外徑為3 m、內(nèi)徑為1 m的圓環(huán)面，在圓環(huán)面上按同樣單元大小劃分網(wǎng)格，計(jì)算網(wǎng)格點(diǎn)上的應(yīng)力應(yīng)變，并與解析解對(duì)比計(jì)算其誤差．采用不同大小的單元?jiǎng)澐忠唤M網(wǎng)格以觀察應(yīng)力應(yīng)變場(chǎng)計(jì)算結(jié)果的收斂性．球面網(wǎng)格單元數(shù)分別為48、140、542、1 138、1 920、2 916、11 864和48 192，圓環(huán)面上的節(jié)點(diǎn)數(shù)分別為32、92、356、732、1 032、1 664、7 296和29 186．由于在基于二次元的Nystr?m離散方法中，每個(gè)邊界單元上有6個(gè)Gauss點(diǎn)，每個(gè)Gauss點(diǎn)上有3個(gè)自由度，所以邊界上的自由度為單元數(shù)乘以18．由于應(yīng)力應(yīng)變張量的對(duì)稱性，在求解應(yīng)力邊界積分方程時(shí)，圓環(huán)面上每個(gè)點(diǎn)有6個(gè)自由度，所以圓環(huán)面上總自由度數(shù)為節(jié)點(diǎn)數(shù)乘以6．

本算例在一臺(tái)8核（Xeon 5355, 2.66 GHz）、32 GB內(nèi)存的工作站上計(jì)算．在KIFMM中采用不同的上行等效點(diǎn)數(shù)p（p= 4,5,6,7），計(jì)算了應(yīng)力應(yīng)變場(chǎng)，得到的誤差與邊界上自由度N的關(guān)系如圖5和圖6所示，其中，conventional是用未加速的傳統(tǒng)邊界元法計(jì)算的結(jié)果．可以看出，在選用合適的等效點(diǎn)數(shù)時(shí)，采用本文的KIFMM方法可以有效地保持傳統(tǒng)邊界元法的精度，表明了本文理論和方法的正確性．

圖5 圓環(huán)面上的應(yīng)力誤差

圖6 圓環(huán)面上的應(yīng)變誤差

4 結(jié) 論

本文推導(dǎo)證明了基于等效源展開的核無(wú)關(guān)快速多極方法（KIFMM）等同于采用了一種特殊的核函數(shù)展開，從而它與解析FMM可以統(tǒng)一到相同的理論框架下，證明了KIFMM中等效源近似與解析FMM中核函數(shù)解析展開相似．在此基礎(chǔ)上，將此理論應(yīng)用于彈性動(dòng)力學(xué)的應(yīng)力應(yīng)變場(chǎng)分析，驗(yàn)證了本文理論和方法的正確性．