逄澤慧，張耀明

(山東理工大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,淄博 255049)

1 引言

許多科學(xué)與工程問(wèn)題可歸結(jié)為求解位勢(shì)問(wèn)題或者與求解位勢(shì)問(wèn)題密切相關(guān)。如穩(wěn)定滲流、穩(wěn)定熱傳導(dǎo)、動(dòng)水壓力、薄膜平衡、彈性桿件的扭轉(zhuǎn)、波傳播、電磁場(chǎng)、非均質(zhì)及非線性問(wèn)題等[1,2]。求解上述問(wèn)題的主流數(shù)值方法是有限差分法、有限元法及邊界元法。這些方法需求助于網(wǎng)格框架，來(lái)構(gòu)造基本物理量的插值逼近函數(shù)和在其上進(jìn)行能量積分。然而，對(duì)于復(fù)雜計(jì)算域問(wèn)題，網(wǎng)格產(chǎn)生的前處理是相當(dāng)困難和耗時(shí)的。

為了消除或者緩解基于網(wǎng)格的數(shù)值方法的這些缺陷，無(wú)網(wǎng)格法應(yīng)運(yùn)而生。根據(jù)是否需要域內(nèi)配點(diǎn)，可分為區(qū)域型和邊界型兩大類(lèi)。邊界型無(wú)網(wǎng)格法相對(duì)更具競(jìng)爭(zhēng)力，因?yàn)槠淅^承了邊界型數(shù)值法,如邊界元法的降維及易于求解外邊值問(wèn)題的優(yōu)點(diǎn)。目前，已發(fā)展了多種邊界型無(wú)網(wǎng)格法，如基本解法MFS[3]、邊界配點(diǎn)法BCM[4]、邊界結(jié)點(diǎn)法BKM[5]、修改的基本解法MMFS[6]、規(guī)則化無(wú)網(wǎng)格法RMM[7]、邊界分布源法BDSM[8]、奇異邊界法SBM[9]、節(jié)點(diǎn)法BNM[10]及其各種衍生法[11]等。需指出的是，對(duì)于依賴(lài)邊界積分方程的邊界型無(wú)網(wǎng)格法，關(guān)鍵問(wèn)題是如何準(zhǔn)確估計(jì)影響矩陣的對(duì)角元，但是解決該問(wèn)題非常困難，因?yàn)樾枰幚砘窘饧捌鋵?dǎo)數(shù)引起的奇異性。

張耀明等[12-14]提出了一種新的邊界型無(wú)網(wǎng)格法，即平均源邊界節(jié)點(diǎn)法ASBNM。ASBNM法耦合平均源技術(shù)和完全的規(guī)則化邊界積分方程，不涉及單元和積分的概念，具有方法簡(jiǎn)單和精度高等優(yōu)點(diǎn)。本文改進(jìn)了文獻(xiàn)[13]的方法，提出了直接計(jì)算影響矩陣對(duì)角元的方法，將影響矩陣對(duì)角元的計(jì)算轉(zhuǎn)化為與求解區(qū)域的邊界有關(guān)的純幾何問(wèn)題，因而方法適用于任何二維邊值問(wèn)題。

2 邊界積分方程

假定Ω是R2的有界區(qū)域,Ωc是其補(bǔ)域，Γ=?Ω是兩者的共同邊界。

2.1 邊界值問(wèn)題

邊界條件

(2)

(3)

|u(x)|=O(1)

u*(x,y)=-ln|x-y|/2π

(5)

2.2 間接邊界積分方程

(6)

(y∈Γ)(7)

(y∈Γ)(8)

3 平均源邊界節(jié)點(diǎn)法

為了避免估計(jì)方程(6)的對(duì)角元時(shí)遇到的弱奇異性，將方程改寫(xiě)成規(guī)則化的形式為

(y∈Γ)(11)

基于平均源技術(shù)[12-14]，方程(6,8,11)可歸結(jié)為如下平均源邊界節(jié)點(diǎn)法公式

(12～14)

(15)

(16)

式中ri k是從第i個(gè)節(jié)點(diǎn)xi到第k個(gè)節(jié)點(diǎn)xk的向量，ri k是第i個(gè)節(jié)點(diǎn)xi與第k個(gè)節(jié)點(diǎn)xk之間的距離，nk是區(qū)域Ω的邊界Γ在第k個(gè)節(jié)點(diǎn)xk處的單位外法向量，J(xi)是第i個(gè)節(jié)點(diǎn)xi處的雅克比，κ(xi)是邊界Γ在節(jié)點(diǎn)xi處的曲率。式(16)的推導(dǎo)詳見(jiàn)附錄。

4 數(shù)值算例

通過(guò)不連續(xù)邊界條件及多聯(lián)通區(qū)域問(wèn)題的算例來(lái)驗(yàn)證方法的有效性和準(zhǔn)確性，并與邊界元法比較，邊界元法使用精確幾何單元[16]和常元插值計(jì)算。為了評(píng)估結(jié)果的準(zhǔn)確性，定義相對(duì)誤差為

(19)

算例1考慮方域[0,π]×[0,π]上的混合邊值問(wèn)題[7]，如圖1所示。邊界條件是

u(π,x2)=1,q(x1,0)=u(x1,π)=q(0,x2)=0

(20)

問(wèn)題的解析解是

式中

(21)

圖2分別給出了場(chǎng)位勢(shì)u及其通量?u/?x1的數(shù)值解的絕對(duì)誤差曲面，誤差曲面產(chǎn)生于均勻分布在方形區(qū)域[0.15,3.0]×[0.15,3.0]上的25×25個(gè)計(jì)算點(diǎn)，使用100個(gè)節(jié)點(diǎn)計(jì)算。圖3描述了場(chǎng)溫度u和場(chǎng)梯度?u/?x1的相對(duì)誤差隨邊界節(jié)點(diǎn)數(shù)的變化。

圖1 計(jì)算域

Fig.1 Problem sketch

圖2 場(chǎng)位勢(shì)u及其導(dǎo)數(shù)?u/?x1數(shù)值解的絕對(duì)誤差曲面(100常節(jié)點(diǎn))

Fig.2 Absolute error surfaces for field solutions and fluxes(100 nodes)

圖3 場(chǎng)位勢(shì)u和場(chǎng)梯度?u/?x1的收斂曲線

Fig.3 Relative errors for the field potentialuand its derivative ?u/?x1vs.the number of boundary nodes

算例2考察多聯(lián)通區(qū)域上的混合邊值問(wèn)題，如圖4所示。問(wèn)題的解析解是

u=r3cos(3θ)

(22)

求解時(shí)，各個(gè)邊界上的邊界節(jié)點(diǎn)數(shù)之比是Γ1∶Γ2∶Γ3∶Γ4∶Γ5∶Γ6=5∶1∶1∶1∶1∶1。邊界上均布 120個(gè) 節(jié)點(diǎn)，10個(gè)計(jì)算點(diǎn)均勻分布在以原點(diǎn)O為圓心，0.5 為半徑的圓周上，如圖5所示。表1列出了 10個(gè) 計(jì)算點(diǎn)處的場(chǎng)函數(shù)u的數(shù)值結(jié)果。圖6 描述了邊界Γ3～Γ6上通量?u/?n數(shù)值解的相對(duì)誤差隨邊界節(jié)點(diǎn)數(shù)的變化及邊界元解的比較。

圖4 計(jì)算域(R=2,r=0.25和a=1)

Fig.4 Problem sketch(R=2,r=0.25 anda=1)

圖5 邊界節(jié)點(diǎn)分布及計(jì)算點(diǎn)

Fig.5 Computational model,where symbols(○)and(●)represent calculation points and boundary nodes,respectively

表1 內(nèi)點(diǎn)位勢(shì)u數(shù)值結(jié)果Tab.1 Numerical results for field potential

圖6 邊界Γ3～Γ6上法向通量?u/?n的收斂曲線

Fig.6 Relative errors for the boundary quantities vs.the number of nodes: the normal flux onΓ3～Γ6

5 結(jié) 論

本文提出直接計(jì)算平均源邊界節(jié)點(diǎn)法影響矩陣對(duì)角元的新方法。該方法不使用與問(wèn)題有關(guān)的規(guī)則化邊界積分方程，而是將對(duì)角元計(jì)算的問(wèn)題轉(zhuǎn)化為一個(gè)解域邊界的純幾何問(wèn)題，因此該方法不僅簡(jiǎn)單易實(shí)施，而且可直接推廣到其他平面問(wèn)題，如彈性力學(xué)問(wèn)題和Stokes問(wèn)題等。

附 錄

因此

利用平均源技術(shù)，可得