《計(jì)算物理學(xué)》課件第11章

11.1邊界積分方程11.2邊界元近似11.3單一邊界下的邊界元法計(jì)算程序11.4兩種介質(zhì)情況下的邊界元方法習(xí)題十一11.1邊界積分方程11.1.1面積分化為邊界積分

定理11.1

若u(x,y)和u*(x,y)是定義在平面域D上的兩個(gè)任意函數(shù)，它們在邊界Γ外法線上的導(dǎo)數(shù)分別為

(如圖11.1所示)，則有

(11.1)

證明由格林公式可知

(11.2(a))

(11.2(b))

由式(11.2(b))~(11.2(a))即可得到式(11.1)。圖11.1邊界積分示意圖11.1.2關(guān)于δ函數(shù)的基本知識

1.δ函數(shù)的定義

滿足關(guān)系式

的函數(shù)稱為δ函數(shù)。δ函數(shù)的一個(gè)基本性質(zhì)是對于一個(gè)連續(xù)函數(shù)f(x)都有

2.定義在平面域上的函數(shù)

定義在平面域上的函數(shù)，有

其中,r為動(dòng)點(diǎn)到定點(diǎn)的距離(如圖11.2所示)。以下對上式作一簡單說明。圖11.2r定義示意圖在平面極坐標(biāo)系中有

(1)若r≠0，則;

(2)若r=0，則;

(3)計(jì)算，其中D是包含定點(diǎn)(i點(diǎn))的任意區(qū)域(如圖11.3所示)。圖11.3邊界的積分示意圖

圖11.4積分的局部示意圖

由圖11.4所示可知

(11.3)

由于

而

所以有

代入式(11.3)可得

(11.4)

根據(jù)式(11.1)~(11.3)的結(jié)果可知

(11.5(a))

(11.5(b))11.1.3邊界積分方程

設(shè)u*=lnr，則有

(11.6)

而根據(jù)泊松方程有

(11.7)將式(11.6)和式(11.7)代入式(11.1),有

(11.8)其中r=|rP-ri|，P為邊界上的動(dòng)點(diǎn)。令上式中等號右邊第二項(xiàng)為

(11.9)式(11.9)可以通過高斯求積法獲得，而對于式(11.8)中等號左邊可以表示為

(11.10)圖11.5

如圖11.5所示，上式中的系數(shù)Ci定義為

(11.11)

將式(11.9)和式(11.10)代入式(11.8)，可得

(11.12)

因?yàn)閡*和q*是已知的，所以上式給出了區(qū)域D和邊界Γ上任意一點(diǎn)與邊界Γ上的u、q

之間的關(guān)系。這就是與泊松方程對應(yīng)的邊界積分方程，它是邊界元法的基礎(chǔ)。11.2.1常數(shù)邊界元近似

為了計(jì)算(11.12)式中的邊界積分，將邊界Γ分成N段，假設(shè)在每一段上u、q近似為常數(shù)并等于該段中點(diǎn)處的值。對于給定的點(diǎn)i，(11.12)式可離散為

(11.13)11.2邊界元近似令

(11.14)

Hij和Gij均可以通過計(jì)算得到，式(11.13)可重寫為

(11.15)11.2.2邊界元方程

(1)若定點(diǎn)i在邊界Γ上(Ci=π)，則ui必定是{uj}中的一個(gè)。式(11.15)可轉(zhuǎn)化成

(11.16)可見上式為一個(gè)N×N的線性方程組。

以下針對式(10.9)~(10.11)所滿足的泊松方程定解問題進(jìn)行邊界元方法的介紹。這里Γ由Γ1和Γ2組成，，

?？梢詫⑹?11.16)表示為

(11.17)顯然,上式中當(dāng)j∈Γ1時(shí)，qj是未知的，而uj已知的；當(dāng)j∈Γ2時(shí)，qj是已知的，而uj是未知的。將上式中含未知量的項(xiàng)移到等號左端，而將已知項(xiàng)移到等號右端可得

(11.18)將式(11.18)可以寫成如下矩陣形式

AX=R

(11.19)

其中,

采用高斯消元法解方程(11.19)，可求得和。(2)若定點(diǎn)i在區(qū)域D內(nèi)(Ci=2π)，根據(jù)式(11.15)和已經(jīng)求出的u|Γ和q|Γ，可求得區(qū)域內(nèi)的ui。11.2.3矩陣H、G

和B的計(jì)算

1．Hii和Gii的計(jì)算

如圖11.6所示，在某一段邊界元Γi上取一段線元ds，其法向?yàn)閚。當(dāng)Γi較小時(shí)，可認(rèn)為n與該邊界元垂直，即圖中的r不隨n變化，故積分

圖11.6Γi示意圖因此有

Hii=-Ci=-π(i∈Γ光滑點(diǎn))(11.20)

(11.21)

其中,si為線元Γi的長度。

2．Hij和Gij的計(jì)算

如圖11.7所示，rd是由點(diǎn)i到線元Γj的距離，與n同向?yàn)檎聪驗(yàn)樨?fù)。規(guī)定外邊界Γ為逆時(shí)針走向，內(nèi)邊界Γ為順時(shí)針走向(n與Γ方向滿足右手螺旋法則)，如圖11.8所示。圖11.7法線方向與相關(guān)變量的定義

圖11.8邊界Γ的走向利用類似式(11.4)的推導(dǎo)可得

(11.22)

式中,θj是區(qū)域內(nèi)定點(diǎn)i到線元Γj的張角(即關(guān)于動(dòng)點(diǎn)j在線元始、末點(diǎn)所張開的角),d是rd的垂足到Γj起點(diǎn)的距離。

(11.23)由于，代入上式后可得

(11.24)式(11.23)和式(11.24)中的d、θj、sj、r1、r2、rd可由定點(diǎn)i坐標(biāo)(xi,yi)和Γj兩端點(diǎn)坐標(biāo)(x1,y1)、(x2,y2)計(jì)算求得，即

3．Bi的計(jì)算

把區(qū)域D剖分成有限個(gè)三角形單元，利用高斯求積法計(jì)算各單元的積分，即

(11.25)需要說明的是,當(dāng)泊松方程的邊界條件全部為第一類或第二類邊界條件時(shí)，可直接利用式(11.16)計(jì)算出邊界上未知的q或u，進(jìn)而根據(jù)式(11.15)計(jì)算出區(qū)域D內(nèi)的u。當(dāng)然，矩陣H、G的計(jì)算方法不變。

【例11.1】

設(shè)有泊松方程的混合邊值問題(如圖11.9所示)圖11.9例11.1示意圖

取每個(gè)邊為邊界元。圖中①、②、③三點(diǎn)為三角形各邊的中點(diǎn)。試寫出邊界積分方程并采用常數(shù)邊界法計(jì)算圖中點(diǎn)③的電勢和①、②兩點(diǎn)的電量。

解取每個(gè)邊即為邊界元。

(1)計(jì)算對角線元素Hii、Gii

。根據(jù)式(11.20)和(11.21),有

計(jì)算非對角線元素Hij、Gij(i≠j)。根據(jù)式(11.22)有

顯然Hij≠Hji。因此有

而根據(jù)式(11.24)有G31=(1-0)ln

-1+0.5×1.107=-0.335,Gij≠Gji，經(jīng)過計(jì)算可得

Gq(2)形成方程AX=R。根據(jù)式(11.19)有

A=其中A=X=R=即有

(3)用高斯消元法解矩陣方程，求出uq=(q1,q2,u3)=(-2.580,1.730,0.825)。

(4)對邊界數(shù)據(jù)進(jìn)行整理。

為下節(jié)編程需要，我們將最后求出的uq=(q1,q2,u3)中的u3和邊界上已知的uq0=(u1,u2,q3)中的q3互換，即將(q1,q2,q3)存入uq，(u1,u2,u3)存入uq0。仍然考慮式(11.9)~(11.11)所滿足的泊松方程定解問題，即

11.3單一邊界下的邊界元法計(jì)算程序?qū)τ趩我贿吔鐥l件下的泊松方程邊界元計(jì)算步驟如下：由輸入信息、計(jì)算矩陣、引入邊界條件，形成方程到解邊界元方程求邊界上的u、q，再到求區(qū)域D內(nèi)的u，最后輸出計(jì)算結(jié)果。

主程序流程如下:

(1)輸入信息。

dimensionxm(40),ym(40),xd,yd,ud(40),

*H(40,40),G(40,40),x(40),y(40),gama(40),uq0(40),uq(40)!(x,y)為邊界上各段端點(diǎn)坐標(biāo)，(xm,ym)為邊界上各段中點(diǎn)的坐標(biāo)，(xd,yd)為內(nèi)點(diǎn)坐標(biāo),!gama(j)=0(j∈Γ1),gama(j)=1(j∈Γ2);uq0=u0j(j∈Γ1),uq0=q0j(j∈Γ2)

(2)具體計(jì)算。

callmatrix(n,x,y,gama,uq0,G,uq)

do10i=1,n

callbi(xm(i),ym(i),B)

10

uq(i)=uq(i)+Bcallgauss(n,G,uq)

write(*,*)(uq0(i),i=1,n)

write(*,*)(uq(i),i=1,n)

callinter(n,xd,yd,x,y,uq0,uq,ud)

callbi(xd,yd,B)

ud=ud+B/6.2832

write(*,*)xd,yd,ud

(3)輸出結(jié)果。

輸出n個(gè)邊界中點(diǎn)處的坐標(biāo)及u，q值xm(j),ym(j),uq0(j),uq(j),j=1,2,…,n，輸出內(nèi)點(diǎn)坐標(biāo)和u值xd,yd,ud

子程序流程如下：

(1)計(jì)算矩陣，形成方程。

subroutinematrix(n,x,y,gama,uq0,G,uq)

dimensionxm(n),ym(n),H(n,n),G(n,n)

dimensionx(40),y(40),gama(40),uq0(40),uq(40)

x(n+1)=x(1)

y(n+1)=y(1)

do10j=1,n

xm(j)=0.5*(x(j)+x(j+1))

10ym(j)=0.5*(y(j)+y(j+1))

！計(jì)算H(i,i),G(i,i)

do30i=1,n

H(i,i)=-3.14159

s=((x(i+1)-x(i))2+(y(i+1)-y(i))2)1/2

G(i,i)=s*(lns-1.69315)

do30j=1,n

If((i-j).eq.0)goto30

callHgij(xm(i),ym(i),x(j),y(j),x(j+1),y(j+1),H(i,j),G(i,j))

30continue以下為形成方程AX=R，將A和R分別存入G和uq

當(dāng)j∈Γ2時(shí)，將-Hij存入Gij(形成A)，而將-Gij存入Hij(形成R

do50j=1,n

if(gama(j))50,50,40

40do55i=1,n

ch=G(i,j)

G(i,j)=-H(i,j)!對應(yīng)于形成的矩陣A

H(i,j)=-ch

55continue

50continue

do60i=1,n

uq(i)=0

do70j=1,n70uq(i)=uq(i)+H(i,j)*uq0(j)!對應(yīng)于形成的矩陣R

60continue

return

end

！計(jì)算H(i,j)，G(i,j)的子程序

subroutineHGij(xi,yi,x1,y1,x2,y2,H,G)

x21=x2-x1

y21=y2-y1

x1i=x1-xi

y1i=y1-yi

x2i=x2-xi

y2i=y2-yi

s=sqrt((x21*x21+y21*y21))

d=-(x1i*x21+y1i*y21)/s

r1=sqrt(x1i*x1i+y1i*y1i)

r2=sqrt(x2i*x2i+y2i*y2i)

rd=(x1i*y21-y1i*x21)/s

H=acos((x1i*x2i+y1i*y2i)/(r1*r2))

G=(s-d)*lnr2+d*lnr1-s+rd*H

return

end

(2)求區(qū)域D內(nèi)的u。

subroutineinter(n,xd,yd,x,y,uq0,uq,ud)

dimensionx(40),y(40),gama(40),uq0(40),uq(40)

do20j=1,n

if(gama(j))20,20,10

10ch=uq0(j)

uq0(j)=uq(j)

uq(j)=ch20continue

x(n+1)=x(1)

y(n+1)=y(1)

ud=0

do30j=1,n

callHGij(xd,yd,x(j),y(j),x(j+1),y(j+1),A,B)

30ud=ud+A*uq0(j)-B*uq(j)

ud=ud/6.2832

return

end如圖11.10(a)所示，區(qū)域1和2分別有介電常數(shù)為ε1和ε2的兩種均勻介質(zhì)。假設(shè)兩區(qū)域內(nèi)沒有自由電荷，則電勢u滿足拉普拉斯方程，2u=0，在邊界元方程中，對應(yīng)的邊界元方程為GQ=HU，如式(11.16)。11.4兩種介質(zhì)情況下的邊界元方法當(dāng)我們解邊界元方程時(shí)，不僅要利用外邊界Γ1和Γ2上的邊界條件，還要利用內(nèi)邊界ΓI上的邊界條件。為此，引入下列符號

圖11.10區(qū)域1與區(qū)域2的劃分

將兩個(gè)區(qū)域的邊界元方程寫成如下形式:

(11.26)

(11.27)各矩陣的上標(biāo)代表區(qū)域，下標(biāo)代表邊界。在兩接著的交界處，電勢u和電位移法向分量εr連續(xù)，所以

(11.28)

(11.29)利用式(11.28)、(11.29)和式(11.26)、(11.27)合成如下形式：

(11.30)這就是有兩種介質(zhì)的邊界元方程，它已經(jīng)包含了內(nèi)邊界(交界)條件。進(jìn)一步要利用外邊界條件，把已知量移到方程的右邊，未知量移到左邊,變AX=R的形式。這些過程與11.2節(jié)所討論的相同。但應(yīng)該注意編號規(guī)則。兩個(gè)區(qū)域的邊界元都按逆時(shí)針向編號，并且兩個(gè)區(qū)域的編號連起來，如圖11.10(b)所示。設(shè)邊界Γ1、Γ2、ΓI上的邊界元數(shù)分別為N1、N2、NI，則方程(11.30)中各個(gè)矩陣在總矩陣中的位置如表11.1所示，其中N=N1+N2+2NI。表11.1各個(gè)矩陣在總矩陣中的位置程序流程如下：

(1)輸入數(shù)據(jù):

subroutineinput(N,εr)

commonN1,N2,NI,KAMA(40),uq0(40),uq(40),x(40),y(40)

dataN1,N2,NI,εr/…/

N=N1+N2+2NI

兩區(qū)邊界元數(shù)，N<40邊界元端點(diǎn)坐標(biāo)和邊界條件：

x(j),y(j),KAMA(j),uq0(j)

j=1…,N1,…,N1+NI,…,N1+2NI,…,N

Γ1

Γ2

(在ΓI上，KAMA和uq0無值)

(2)計(jì)算矩陣，形成方程:

subroutinematrix(N,εr,xm,ym,G,H)

dimensionxm(N),ym(N),G(N,N),H(N,N),J0(2),J1(2),J2(2)

commonN1,N2,NI,KAMA(40),uq0(40),uq(40),x(40),y(40)

兩個(gè)區(qū)的始末編號

J1(1)=1,J2(1)=N1+NI,J1(2)=J2(1)+1,J2(2)=N

L0=2*(N1+NI)+1,x(N+1)=x(1),y(N+1)=y(1)

do10j=1,N

xm(j)=(x(j)+x(j+1))/210ym(j)=(y(j)+y(j+1))/2

do100K=1,2

!分兩個(gè)區(qū)算G和H

K1=J1(K),K2=J2(K)!兩個(gè)區(qū)的始末編號

do100i=K1,K2

s=((x(i+1)-x(i)))2+(y(i+1)-y(i))2)1/2,uq(i)=0

H(i,i)=-3.14159,G(I,i)=s(lns-1.69315)do100j=K1,K2

J0(1)=j-N1,J0(2)=N1+2NI+1-j

if(i-j)20,30,20

20callHgij(xm(i),ym(i),x(j),y(j),x(j+1),y(j+1),H(i,j),G(i,j))

利用邊界條件Γ1和Γ2，通過移項(xiàng)把方程(11.30)變成AX=R，但A和B分別存入G和uq。

30if(J0(K))40,40,70

若J0≤0，則j∈

,

40

if(KAMA(j))60,60,50

50

ch=G(i,j),G(i,j)=-H(i,j),H(i,j)=-ch

60

uq(i)=uq(i,j)+H(i,j)*uq0(j)

goto100

若J0>0，則j∈Γ1,Γ2

70

L=L0-j若j∈則L∈

if(K-1)80,80,90

把-存入的位置80

G(i,L)=-H(i,j)

goto100

把-

/εr存入的位置

90

G(i,L)=-G(i,j)/εr

把-存入的位置

G(i,j)=-H(i,j)

100continue

return

end

(3)輸出結(jié)果:

subroutineoutput(N,εr,xm,ym)

dimensionJ1(2),J2(2),xm(N),ym(N)

commonN1,N2,NI,KAMA(40),uq0(40),uq(40),x(40),y(40)

調(diào)用subroutinegauss(N,G,uq)后，uq存有求出的u和q，而uq0中存有已知的u和q。需要進(jìn)行整理，把u和q分別放在uq0和uq。整理Γ1和Γ2上的u和q

J1(1)=1,J2(1)=N1!Γ1的始末編號

J1(2)=N1+2NI+1,J2(2)=N

!Γ2的始末編號

do30K=1,2

K1=J1(K),K2=J2(K)

do30j=K1,K2

if(KAMA(j))30,30,20

20ch=uq0(j),uq0(j)=uq(j),uq(j)=ch

30continue整理Γ1上的u和q。由(11.30)式知，uq中求出的是和,需先給出和

I1=N1+1,I2=N1+NI,的始末編號