陳林沖

(重慶師范大學 數(shù)學科學學院， 重慶 401331)

?

二維Helmholtz方程的邊界點解法

陳林沖

(重慶師范大學 數(shù)學科學學院， 重慶 401331)

針對二維Helmholtz方程的混合邊值求解問題，采用邊界點方法(boundary node method)，在直接邊界積分方程的基礎(chǔ)上，建立了求解Helmholtz方程邊值問題的正則化形式，有效地避免了強奇異積分的計算，并且推導了弱奇異積分的計算公式。兩個數(shù)值算例表明本方法可取得較高的可行性和有效性。

二維Helmholtz方程；混合邊值問題；邊界點法；強奇異積分；弱奇異積分

Helmholtz方程在工程技術(shù)、電磁場理論、散射理論、力學等較多領(lǐng)域有著廣泛的應用，研究其數(shù)值解具有廣泛的實際意義和重要的理論價值，尤其是對混合邊值問題的處理[1]。傳統(tǒng)的數(shù)值計算方法有有限差分法、有限元法、邊界元法等。

有限差分法是將偏微分方程的定解問題按差分格式離散求解，但由于數(shù)值精度不穩(wěn)定，很少用于物理工程問題的處理。有限元法將問題所在的區(qū)域離散為有限單元，在單元上對求解變量建立基于節(jié)點的插值函數(shù)近似，以變分原理或加權(quán)殘數(shù)法建立控制方程，聯(lián)立數(shù)值積分建立相關(guān)的代數(shù)方程組，但是它需要在整個區(qū)域上進行離散，在處理大變形以及動態(tài)裂紋擴展等問題時，網(wǎng)格可能會發(fā)生畸變，此時必須進行網(wǎng)格重構(gòu)，從而損失了精度。邊界元法[2]是繼有限元法之后的一種基于邊界單元和節(jié)點的數(shù)值方法，是將偏微分方程的定解問題轉(zhuǎn)化為邊界積分方程，并對求解區(qū)域的邊界進行單元離散而形成的數(shù)值方法，必須要建立與插值節(jié)點拓撲相關(guān)的網(wǎng)格。對于復雜方程，不僅需要克服基本解的奇異性，還要進行必要的網(wǎng)格重構(gòu)，增加計算時間，降低精度。因此，為了克服基于網(wǎng)格(單元和節(jié)點)的數(shù)值方法對于單元或網(wǎng)格的依賴，保證計算精度,節(jié)約計算時間，產(chǎn)生了一種不需要與節(jié)點相關(guān)聯(lián)網(wǎng)格的新型數(shù)值計算方法——無網(wǎng)格方法，并且近年來得到了較大的發(fā)展[3-4]。

最早的無網(wǎng)格方法是1977年Lucy提出的光滑粒子法[5]，主要用于研究無邊界的天體問題和流體問題，但穩(wěn)定性不是很好。1981年，Lancaster等建立了移動最小二乘法[6]，采用完備多項式作為基函數(shù)，利用誤差的加權(quán)平方和構(gòu)造泛函，得到的逼近函數(shù)精度高、光滑性好，函數(shù)本身和其導函數(shù)都連續(xù)。20世紀90年代，Nayroles和Belytschko等將移動最小二乘法引入微分方程的邊值問題中[7]，后來Mukherjee等對移動最小二乘法進行了改進，提出了求解勢問題和彈性力學問題的邊界點法(Boundary Node Method)[8-11]。邊界點法是將移動最小二乘法導出的插值公式與邊界積分方程結(jié)合起來，通過邊界積分方程的邊值條件，建立邊界上的線性方程組,計算時間快，效率高。

本文致力于用邊界點方法數(shù)值求解二維Helmholtz方程。利用格林公式(Green)得到Helmholtz方程的直接邊界積分方程，相比間接的層位勢理論建立的間接邊界積分方程更適用于混合問題。由于二維Helmholtz方程基本解為特殊函數(shù)，將其展開成了一種級數(shù)形式，巧妙地利用二維Laplace方程的基本解的特性，建立了二維Helmholtz方程的正則化的邊界積分方程，有效地避免了強奇異積分的計算。針對其中含基本解的弱奇異積分項，詳細地推導了其數(shù)值計算公式[12]。最后給出了數(shù)值算例,驗證了該方法的可行性和有效性。

1 Helmholtz方程混合邊值問題

1.1 直接邊界積分方程及其正則化形式

考慮如下二維Helmholtz方程混合邊值問題：

利用格林公式得到Helmholtz方程的解的表達式和邊界積分方程[13]：

(1)

(2)

其中c≈0.577 215 664 9，為歐拉常數(shù)，則有

定理1 當kr很小時,二維Helmholtz方程的正則化形式為

(3)

證明 在式(2)中代入其展開級數(shù)得

所以有

因此

定理得證。

1.2 邊界點方法求解二維Helmholtz方程

其中，φl(x)為移動最小二乘法(MLS)的形函數(shù)[15-17]。

設(shè)在邊界Γu上N1有個邊界節(jié)點，則Γq上有N-N1個邊界節(jié)點，代入直接邊界條件形成線性方程組

在式(1)中，?y∈Ω有

(4)

1.3 積分的計算

將邊界Γ劃分成N個積分背景網(wǎng)格Γn，n=1,…,N。Γn選為直線段，可表示為

其中，

(xp,yp)∈Γn，所以有

則有

對于后一個積分I2，采用分部積分公式得

所以

2 數(shù)值算例

為了誤差估計和收斂階的研究，定義下面的評估方式[18]：

圖的誤差

從圖1中可以看出：當k取不同值時，函數(shù)本身和它的一階導函數(shù)的收斂階都在2左右；從圖2中可以看出：當k取不同值時，函數(shù)本身和它的一階導函數(shù)的收斂階都在3左右，都取得了較好的數(shù)值效果。

圖的誤差

3 結(jié)束語

在已有的直接邊界積分方程理論的基礎(chǔ)上，本文建立了二維Helmholtz方程混合邊值問題的正則化邊界積分方程，有效地避免了強奇異積分的計算，引入了弱奇異積分的處理方法。采用邊界點方法(Boundary Node Method)得到了較高的精度和收斂階，2個數(shù)值算例驗證了該方法的可行性和有效性。進一步地，可以發(fā)現(xiàn)由于正方形區(qū)域的特殊性，它所布置的節(jié)點剛好落在正方形區(qū)域的邊界上，沒有產(chǎn)生幾何誤差，因此其得到的收斂階要比圓域得到的收斂階更高一些。

[1] ERLANGGA Y A.A Robust and Efficient Iterative Method for the Numerical Solution of the Helmholtz Equation[D].Delft:Technische Universit Delft,2005.

[2] 嵇醒，臧躍龍，程玉民.邊界元進展及通用程序[M].上海：同濟大學出版社，1997.

[3] BELYTSCHKO T,KRONGAUZ Y,ORGAN D.Meshless methods:An overview and recent developments[J].Computer Methods in Applied Mechanics and Engineering,1996,139(1/4):3-47.

[4] YAGAWA G,FURUKAWA T.Recent developments of free mesh method[J].International Journal for Numerical Methods in Engineering,2000,47(8):1419-1443.

[5] LUCY L B.A numerical approach to the testing of the fission hypothesis[J].The Astron J,1977,8(2):1013-1024.

[6] LANCASTER P,SALKAUSKAS K.Surfaces generated by moving least square methods[J].Mathematics of Computation,1981,37:141-158.

[7] NAYROLES B.Generalizing the finite element method:Diffuse approximation and diffuse elements[J].Computational Mechanics,1992,10:307-318.

[8] MUKHERJEE Y X,MUKHERJEE S.The boundary node method for potential problems[J].International Journal for Numerical Methods in Engineering,1997,40(5):797-815.

[9] KOTHNOR VS,MUKHERJEE S,MUKHERJEE Y X.Two dimensional linear elasticity by the boundary node method[J].International Journal of Solids and Structures,1999,36(8):1129-1147.

[10]CHATI M K,MUKHERJEE S,MUKHERJEE Y X.The boundary node method for three-dimensional linear elasticity[J].International Journal for Numerical Methods in Engineering,1999,46(8):1163-1184.

[11]CHATI M K,MUKHERJEE S.The boundary node method for three-dimensional problems in potential theory[J].International Journal for Numerical Methods in Engineering,2000,47(9):1523-1547.

[12]SRINIVAS T,SUBRATA M.An extended boundary node method for modeling normal derivative discontinuities in potential theory across edges and corners[J].Engineering Analysis with Boundary Elements,2004,28:1099-1110.

[13]祝家麟，袁政強.邊界元分析[M].北京：科學出版社，2009.

[14]王竹溪，郭敦仁.函數(shù)論-特殊函數(shù)概論[M].北京：國防工業(yè)出版社，1983.

[15]張雄，劉巖.無網(wǎng)格法[M].北京：清華大學出版社，2004.

[16]孫新志，李小林.復變量移動最小二乘近似在Sobolev空間中的誤差估計[J].應用數(shù)學和力學,2016,37(4):416-425.

[17]程玉民.無網(wǎng)格方法[M].北京：科學出版社，2015.

[18]LI X L,ZHANG S G.Meshless analysis and applications of a symmetric improved Galerkin boundary node method using the improved moving least-square approximation [J].ELSEVIER.Applied Mathematical Modelling,2016(40):2875-2896.

(責任編輯 陳 艷)

Boundary Node Method for the 2-D Helmholtz Equation

CHEN Lin-chong

(School of Mathematical Science, Chongqing Normal University, Chongqing 401331, China)

Aiming at 2D Helmholtz equation mixed boundary problem, based on the boundary node method, the regularized form to 2D Helmholtz equation boundary value problem is given on the foundation of direct boundary integral equation. The method can effectively avoid strongly singular integral, and it also deduces a calculating formula to weakly singular integral in detail. Two numerical examples demonstrate that this method has the high feasibility and efficiency.

2D Helmholtz equation; mixed boundary problem; boundary node method; strongly singular integral; weakly singular integral

2017-04-18

國家自然科學基金資助項目(11471063)；重慶市基礎(chǔ)科學與前沿技術(shù)研究重點項目(cstc2015jcyjBX0083)

陳林沖(1988—)，男，碩士，主要從事計算數(shù)學研究，E-mail:794530653@qq.com。

陳林沖.二維Helmholtz方程的邊界點解法[J].重慶理工大學學報(自然科學)，2017(7):188-194.

format：CHEN Lin-chong.Boundary Node Method for the 2-D Helmholtz Equation[J].Journal of Chongqing University of Technology(Natural Science)，2017(7):188-194.

10.3969/j.issn.1674-8425(z).2017.07.030

O242.2

A

1674-8425(2017)07-0188-07