平面位勢(shì)問(wèn)題的簡(jiǎn)單解邊界元法

許文政，張耀明

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

邊界元法是分析科學(xué)和工程應(yīng)用問(wèn)題的強(qiáng)有力的數(shù)值模擬工具，廣泛應(yīng)用于位勢(shì)流，彈性靜、動(dòng)力學(xué)，聲傳播等線(xiàn)性問(wèn)題及各種非均勻材料與非線(xiàn)性問(wèn)題[1-6].但是，不同于域型方法，邊界元法需要計(jì)算奇異核積分，它的有效計(jì)算是邊界法成功實(shí)施的關(guān)鍵.時(shí)至今日，這方面的研究仍在進(jìn)行中[1-3].

現(xiàn)有許多有效計(jì)算奇異積分的方法[7-9].大致可分為直接計(jì)算方法和間接計(jì)算方法兩大類(lèi)[10-20].其中，Guiggiani等[10]提出的方法具有代表性，該方法在內(nèi)蘊(yùn)坐標(biāo)系的參數(shù)平面內(nèi)操作，能計(jì)算各種階奇異積分.該方法需要將被積函數(shù)中的每個(gè)量展成單元局部距離Ω的泰勒級(jí)數(shù)，因此需要許多復(fù)雜的推導(dǎo)工作和與單元有關(guān)的計(jì)算.高效偉[11]也提出了一種計(jì)算奇異積分的一般性算法，它將積分核中的非奇異部分表示成在內(nèi)蘊(yùn)坐標(biāo)系統(tǒng)下的局部距離參數(shù)的冪級(jí)數(shù)，然后解析計(jì)算.確定冪級(jí)數(shù)的系數(shù)時(shí)，需要求解線(xiàn)性方程組，而且這樣的操作每個(gè)單元都需要.間接計(jì)算方法主要是通過(guò)建立各種規(guī)則化邊界積分方程，來(lái)間接計(jì)算CPV和HFP積分.包括直接變量規(guī)則化邊界積分方程和間接變量規(guī)則化邊界積分方程[7-9,15-20].

不同于已有的工作，本文提出了一個(gè)新的計(jì)算強(qiáng)奇異積分的方法，稱(chēng)為簡(jiǎn)單解法.它利用簡(jiǎn)單解過(guò)程，獲得強(qiáng)奇異積分值，無(wú)需直接計(jì)算積分，也無(wú)需建立規(guī)則化邊界積分方程，因此方法具有理論簡(jiǎn)單、計(jì)算效率高、結(jié)果準(zhǔn)確等優(yōu)點(diǎn).此外，方法能夠計(jì)算任何邊界通量?u/?xi(i=1,2).該方法容易跟隨和拓展到其他如彈性力學(xué)問(wèn)題、Stokes問(wèn)題、Helmholtz問(wèn)題等.

1 基本概念

本文假定Ω是R2中的一個(gè)有界區(qū)域，Ωc是它的補(bǔ)域，Γ=?Ω是它們的邊界.n(x)是區(qū)域Ω的邊界Γ在x點(diǎn)處的單位外法向量.位勢(shì)問(wèn)題的基本解是

(1)

(2)

2 簡(jiǎn)單解法

Ω內(nèi)的位勢(shì)可表示為[17]

(3)

(4)

對(duì)方程(4)，令y→Γ，可獲得如下邊界積分方程

(5)

(6)

為了能夠更清楚地表述問(wèn)題的本質(zhì)，采用方程(6)的離散化形式.將邊界Γ分割成N段(Segment)，并假定y在第l段上，即y∈Γl，并記為yl，于是有

(7)

式(7)右端第一項(xiàng)中的積分均為正常積分，因此問(wèn)題的關(guān)鍵在于如何估計(jì)式(7)右端第二項(xiàng)積分.若記

方程(7)可寫(xiě)為

φ(yl)A,yl∈Γl

(8)

當(dāng)常元插值(常單元)被采用時(shí)，方程(8)可寫(xiě)為

yl∈Γl

(9)

這里φm表示密度函數(shù)φ在第m個(gè)節(jié)點(diǎn)上的值.為了確定(7)式中的A，提出一種簡(jiǎn)單、高效、準(zhǔn)確的簡(jiǎn)單解法：

(1)對(duì)于有限域Ω，引入如下簡(jiǎn)單解

u(x)=x1+x2+1

(10)

(2)對(duì)于無(wú)限域Ωc，引入如下簡(jiǎn)單解

(11)

這里(a,b)是Ωc外的任一點(diǎn).

注：原則上，簡(jiǎn)單解可以任意選取.對(duì)于內(nèi)邊值問(wèn)題，簡(jiǎn)單解只需滿(mǎn)足Laplace方程；對(duì)于外邊值問(wèn)題，簡(jiǎn)單解除滿(mǎn)足Laplace方程外，還需滿(mǎn)足無(wú)窮遠(yuǎn)邊界條件.

方程(4)的離散化形式

l=1,…,N

(12)

這里φm是邊界密度函數(shù)φ(x) 在第m個(gè)單元Γm(m=1,…,N)上的節(jié)點(diǎn)處的值.

對(duì)于有限(無(wú)限)域問(wèn)題，將簡(jiǎn)單解式(10)(或式(11))代入方程(12)中，可求得密度函數(shù)在N個(gè)節(jié)點(diǎn)處的值φm,m=1,…,N.然后，將所求得的密度函數(shù)解φm及簡(jiǎn)單解式(10)(或式(11))代入方程(9)中，即可求得A.

以上方程經(jīng)過(guò)適當(dāng)?shù)慕M合變形后，可適用于任何邊界條件的邊值問(wèn)題.

3 數(shù)值算例

利用三個(gè)數(shù)值算例，以驗(yàn)證本文算法的精確性、有效性及收斂性.所有算例中，邊界離散單元的邊界量采用常元插值，單元幾何邊界采用精確單元L2范數(shù)下的相對(duì)誤差定義為

Relative error =

(13)

算例1首先考察一個(gè)具有Dirichlet不連續(xù)邊界條件的的穩(wěn)定溫度場(chǎng)問(wèn)題.計(jì)算域是[0,1]×[0,1]，邊界條件是

(14)

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

u(x1,x2)=

(15)

其中

(16)

計(jì)算時(shí)，邊界單元數(shù)是80個(gè).圖1(a)和圖1(b)分別描述了解析法和本文數(shù)值法獲得的場(chǎng)位勢(shì)的等值圖，可觀察到，本文數(shù)值解和解析解比較吻合.為了進(jìn)一步考察本文方法的準(zhǔn)確性和有效性，選擇1 600個(gè)計(jì)算點(diǎn)均布在區(qū)域0.15≤x,y≤0.85內(nèi)，圖2給出了這些計(jì)算點(diǎn)上的場(chǎng)溫度數(shù)值解的絕對(duì)誤差曲面.從中可以看出，盡管具有很少的邊界單元數(shù)，本文獲得的數(shù)值結(jié)果仍然相當(dāng)準(zhǔn)確，表明本文給出的方法能夠有效地求解具有不連續(xù)Dirichlet邊界條件的邊值問(wèn)題.意味著，本文計(jì)算強(qiáng)奇異積分對(duì)角元的方法是非常有效的.

(a) 解析解 (b) 數(shù)值解圖1 場(chǎng)位勢(shì)解Fig.1 The field potential solution

圖2 邊界通量解的收斂曲線(xiàn)Fig.2 Convergence curves of boundary flux solutions

(17)

圖3 計(jì)算域Fig.3 Problem sketch

計(jì)算時(shí)，外邊界Γ1離散為36單元，內(nèi)邊界離散為24個(gè)單元.表1給出了沿x軸一些域內(nèi)點(diǎn)處溫度的數(shù)值結(jié)果；圖4(a)和(b)分部描述了邊界Γ1上的法向通量解和邊界Γ2的溫度解的相對(duì)誤差.從中可看出，即使使用很少的單元數(shù)，數(shù)值解已相當(dāng)?shù)販?zhǔn)確.

表1 內(nèi)點(diǎn)溫度解

Tab.1 Temperature solutions on internal points

內(nèi)點(diǎn)精確解數(shù)值解相對(duì)誤差(6.0, 0.0)0.1400000×1030.1399699×1032.149503×10-4(7.0, 0.0)0.1840000×1030.1839418×1033.163785×10-4(8.0, 0.0)0.2340000×1030.2339130×1033.719932×10-4(9.0, 0.0)0.2900000×1030.2898743×1034.333464×10-4

(a) 通量q

(b) 溫度u圖4 邊界量數(shù)值解的相對(duì)誤差Fig.4 Ralative error of boundary temperature uand flux q

算例3考慮不規(guī)則區(qū)域上，具有混合邊界條件的溫度場(chǎng)問(wèn)題，如圖5所示.邊界Γ由兩部分構(gòu)成，即Γ=Γ1∪Γ2，其中

圖5 計(jì)算域Fig.5 Problem sketch

問(wèn)題的解析解是如下振蕩調(diào)和函數(shù)

u(x1,x2)=ex1+3cosx2

(18)

圖6 給出了邊界Γ上的位勢(shì)函數(shù)u的圖像，從圖中可觀察到，u的最大值大于300，最小值小于-20，表明位勢(shì)邊界函數(shù)u具有很強(qiáng)的振蕩性.

圖6 邊界Γ上的位勢(shì)函數(shù)Fig.6 Potential boundary function on Γ

20個(gè)域內(nèi)計(jì)算點(diǎn)均布在中心在(0.5,0.5)，半徑是0.4的圓周上，如圖5所示.圖7描述了這些計(jì)算點(diǎn)上的溫度數(shù)值解u的收斂曲線(xiàn)，圖8給出了邊界Γ1上的通量和Γ2上的溫度數(shù)值解的收斂曲線(xiàn).從中可看出，不論是內(nèi)點(diǎn)還是邊界點(diǎn)，隨著邊界單元數(shù)的增加，相對(duì)誤差迅速減小，可見(jiàn)本文方法在求解具有振蕩解析的問(wèn)題時(shí)仍具有很好的收斂特征.表明本文計(jì)算強(qiáng)奇異對(duì)角元的方法是非常有效的.

圖7 內(nèi)點(diǎn)溫度數(shù)值解的收斂曲線(xiàn)Fig.7 Convergence curves of temperature solutions at internal points

圖8 邊界溫度u和通量q的收斂曲線(xiàn)Fig.8 Convergence curves of boundary temperature uand flux q

4 結(jié)束語(yǔ)

本文提出了計(jì)算強(qiáng)奇異積分的簡(jiǎn)單解法，適用于任何邊界通量的計(jì)算，具有數(shù)學(xué)理論簡(jiǎn)單，程序設(shè)計(jì)容易，計(jì)算準(zhǔn)確性高，容易跟隨和推廣等優(yōu)點(diǎn).數(shù)值算例驗(yàn)證了方法的有效性，表明本文方法能夠有效準(zhǔn)確地計(jì)算強(qiáng)奇異積分的對(duì)角元.