弓小影，閆濤紅，于春肖*

(1. 石家莊經(jīng)濟(jì)學(xué)院 數(shù)理學(xué)院，河北 石家莊 050031；2. 燕山大學(xué) 理學(xué)院，河北 秦皇島 066004)

0 引言

邊界元法是一種精確高效的工程數(shù)值計算方法.二維域上邊界元法的主要思想是將域內(nèi)的微分方程轉(zhuǎn)化為邊界上的積分方程，使得問題的維數(shù)降低一維.通常，轉(zhuǎn)化后的積分方程具有奇異性，為了解決這一問題，陳一鳴等用Galerkin-wavelet法去離散自然邊界積分方程[1-2]，從而使強(qiáng)奇異性減弱，還提出了一種強(qiáng)奇異積分方程的數(shù)值解法.

快速多極算法使邊界元法的大規(guī)模計算成為現(xiàn)實，并應(yīng)用于彈性力學(xué)、聲學(xué)等領(lǐng)域[3-6].姚振漢、王海濤引入指數(shù)展開的定義[7]，得到一種新的展開傳遞格式，使快速多極邊界元中的展開系數(shù)和傳遞關(guān)系得以簡化.Bapat等通過改進(jìn)“鄰居”和“相互作用列表”的概念來減少多極展開到局部展開的工作量[8].Greengard結(jié)合近遠(yuǎn)場劃分準(zhǔn)則對二維位勢問題的多極邊界元方法的多極展開和局部展開的截斷誤差進(jìn)行了分析[9].

本文以二維Stokes flow問題為研究背景，在基本解快速多極展開平移的基礎(chǔ)上，闡述了改進(jìn)“相互作用列表”以后的算法是如何在計算多極展開系數(shù)到局部展開系數(shù)的平移中加速計算的.最后給出了二維Stokes flow問題位移基本解多極展開的截斷誤差估計式.

1 核函數(shù)展開及平移過程

二維Stokes flow問題的邊界積分方程一般形式如下：

Tij(x,y)uj(y))dS(y),

其中：下標(biāo)i,j=1,2;S是求解域的邊界;x，y分別代表邊界上的源點和場點;uj和tj分別是邊界的位移與面力;系數(shù)Cij(x)是與邊界形狀相關(guān)的常數(shù);Uij(x,y),Tij(x,y)是二維Stokes flow問題的Kelvin基本解.

位移的復(fù)變函數(shù)形式的基本解可表達(dá)為：

(1)

平面Stokes flow問題位移基本解的多極展開式為：

(2)

其中：z0為源點,z為場點.利用泰勒級數(shù)在zc展開，zc靠近場點z，且滿足|z-zc|<|z0-zc|，多極展開系數(shù)為

(3)

(4)

(5)

其中Sc為積分域.

如果展開點從zc平移到一個新的點zc′時，得M2M平移公式

(6)

(7)

局部展開系數(shù)可以通過多極展開系數(shù)平移得到，假設(shè)zL是距離源點z0非常近的一點，且滿足|z0-zL|<|zc′-zL|，在傳統(tǒng)多極邊界元方法中，有

(8)

其中，Ll(zL)和Kl(zL)為局部展開系數(shù)，可以通過M2L平移得到

(9)

(10)

當(dāng)展開中心從zL移到zL′時，得L2L平移公式

(11)

(12)

2 主要計算步驟及計算效率簡析

2.1 主要計算步驟

二維Stokes flow多極邊界元方法的主要計算步驟可以分為以下五步：

第一步離散邊界，生成樹結(jié)構(gòu).

第三步下行遍歷計算多極展開系數(shù)的平移，由式(6)和式(7)對每個多極展開系數(shù)進(jìn)行平移，累加得到下一層結(jié)點的多極展開系數(shù)，依次往下對樹結(jié)構(gòu)的l≥2層結(jié)點進(jìn)行計算.

第四步利用式(8)～(10)，上行遍歷計算局部展開系數(shù).

在這一步，對“相互作用列表”進(jìn)行改進(jìn).在傳統(tǒng)算法中，結(jié)點x的“相互作用列表”是一組結(jié)點y的集合，見圖1(a)，結(jié)點y與x不相鄰，但它們的父結(jié)點與x的父結(jié)點相鄰.可以看出，一個結(jié)點的“相互作用列表”中最多有27個結(jié)點.改進(jìn)“相互作用列表”，見圖1(b)，用父結(jié)點z代替4個子結(jié)點y使結(jié)點總數(shù)減少到12，這樣可以降低所需的存儲量，提高計算效率.

第五步由式(11)和式(12)進(jìn)行局部展開系數(shù)的平移，得計算結(jié)果.

(a) (b)

2.2 改進(jìn)前后計算效率簡析

在二維Stokes flow快速多極算法中，第四步占有主要的計算量.在傳統(tǒng)的多極算法中第四步需要最多27次的M2L操作，根據(jù)式(7)，這27次操作的計算量級為O(27p2)，改進(jìn)“相互作用列表”以后，結(jié)點總數(shù)減少到12，總的計算量就可以減少到O(12p2)，提高了計算效率.

3 多極展開截斷誤差分析

3.1 近遠(yuǎn)場劃分準(zhǔn)則

考慮到計算的需要，在多極展開中，一般不會計算到展開式的無窮多項，而是對展開式截取到有限項，這樣就會產(chǎn)生誤差，在討論多極展開的截斷誤差之前首先給出二維問題近遠(yuǎn)場劃分準(zhǔn)則.

定義1 若存在兩組點x1,x2,…,xn∈Ox和y1,y2,…,ym∈Oy(見圖2)和實數(shù)R>0，且

|xi-x0|

|yj-y0|

|x0-y0|>3R,

則稱點集{xi}與{yj}是相互獨立的，或者稱點集{xi}與{yj}相距足夠遠(yuǎn).

圖2 近遠(yuǎn)場劃分示意

用Ox表示中心在x0，半徑為R的圓；Oy表示中心在y0，半徑為R的圓.根據(jù)三角不等式，對任意yj∈Oy,j=1,2,…，m，滿足

‖yj-x0‖≥2R≥2‖xi-x0‖,i=1,2,…,n.

(13)

在進(jìn)行二維Stokes flow多極邊界元方法核函數(shù)的展開平移計算時，對于源點z0，區(qū)域內(nèi)點zL滿足|zL-z0|

3.2 多極展開的截斷誤差分析

定理1 在二維Stokes flow問題核函數(shù)的多極展開式(1)中，當(dāng)截斷項數(shù)到p時，產(chǎn)生的誤差為

證明將式(1)中二維Stokes flow問題的位移基本解寫成

其中，

對于S3，由|z0-z|≤|z0-zc|+|zc-z|，得

z)Ik(z-zc)||t(z)|dS(z)+

(14)

(15)

根據(jù)近遠(yuǎn)場的劃分準(zhǔn)則可得，|z0-zc|≥2R，這樣就有ρ≥2，進(jìn)一步得到

(16)

同理，可得

(17)

于是由式(14)～(17)，可得

且

從而由上面的分析可得

證畢.

4 結(jié)論

通過對多極邊界元計算量級的分析，說明改進(jìn)“相互作用列表”以后的算法加速計算的原理.同時,通過對二維Stokes flow問題的位移基本解多極展開的截斷誤差進(jìn)行分析，驗證了截斷誤差與截斷項數(shù)p的選取有關(guān)，可由截斷項數(shù)p控制.