《计算物理学》课件第11章

11.1边界积分方程11.2边界元近似11.3单一边界下的边界元法计算程序11.4两种介质情况下的边界元方法习题十一11.1边界积分方程11.1.1面积分化为边界积分

定理11.1

若u(x,y)和u*(x,y)是定义在平面域D上的两个任意函数，它们在边界Γ外法线上的导数分别为

(如图11.1所示)，则有

(11.1)

证明由格林公式可知

(11.2(a))

(11.2(b))

由式(11.2(b))~(11.2(a))即可得到式(11.1)。图11.1边界积分示意图11.1.2关于δ函数的基本知识

1.δ函数的定义

满足关系式

的函数称为δ函数。δ函数的一个基本性质是对于一个连续函数f(x)都有

2.定义在平面域上的函数

定义在平面域上的函数，有

其中,r为动点到定点的距离(如图11.2所示)。以下对上式作一简单说明。图11.2r定义示意图在平面极坐标系中有

(1)若r≠0，则;

(2)若r=0，则;

(3)计算，其中D是包含定点(i点)的任意区域(如图11.3所示)。图11.3边界的积分示意图

图11.4积分的局部示意图

由图11.4所示可知

(11.3)

由于

而

所以有

代入式(11.3)可得

(11.4)

根据式(11.1)~(11.3)的结果可知

(11.5(a))

(11.5(b))11.1.3边界积分方程

设u*=lnr，则有

(11.6)

而根据泊松方程有

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

(11.8)其中r=|rP-ri|，P为边界上的动点。令上式中等号右边第二项为

(11.9)式(11.9)可以通过高斯求积法获得，而对于式(11.8)中等号左边可以表示为

(11.10)图11.5

如图11.5所示，上式中的系数Ci定义为

(11.11)

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

(11.12)

因为u*和q*是已知的，所以上式给出了区域D和边界Γ上任意一点与边界Γ上的u、q

之间的关系。这就是与泊松方程对应的边界积分方程，它是边界元法的基础。11.2.1常数边界元近似

为了计算(11.12)式中的边界积分，将边界Γ分成N段，假设在每一段上u、q近似为常数并等于该段中点处的值。对于给定的点i，(11.12)式可离散为

(11.13)11.2边界元近似令

(11.14)

Hij和Gij均可以通过计算得到，式(11.13)可重写为

(11.15)11.2.2边界元方程

(1)若定点i在边界Γ上(Ci=π)，则ui必定是{uj}中的一个。式(11.15)可转化成

(11.16)可见上式为一个N×N的线性方程组。

以下针对式(10.9)~(10.11)所满足的泊松方程定解问题进行边界元方法的介绍。这里Γ由Γ1和Γ2组成，，

。可以将式(11.16)表示为

(11.17)显然,上式中当j∈Γ1时，qj是未知的，而uj已知的；当j∈Γ2时，qj是已知的，而uj是未知的。将上式中含未知量的项移到等号左端，而将已知项移到等号右端可得

(11.18)将式(11.18)可以写成如下矩阵形式

AX=R

(11.19)

其中,

采用高斯消元法解方程(11.19)，可求得和。(2)若定点i在区域D内(Ci=2π)，根据式(11.15)和已经求出的u|Γ和q|Γ，可求得区域内的ui。11.2.3矩阵H、G

和B的计算

1．Hii和Gii的计算

如图11.6所示，在某一段边界元Γi上取一段线元ds，其法向为n。当Γi较小时，可认为n与该边界元垂直，即图中的r不随n变化，故积分

图11.6Γi示意图因此有

Hii=-Ci=-π(i∈Γ光滑点)(11.20)

(11.21)

其中,si为线元Γi的长度。

2．Hij和Gij的计算

如图11.7所示，rd是由点i到线元Γj的距离，与n同向为正，反向为负。规定外边界Γ为逆时针走向，内边界Γ为顺时针走向(n与Γ方向满足右手螺旋法则)，如图11.8所示。图11.7法线方向与相关变量的定义

图11.8边界Γ的走向利用类似式(11.4)的推导可得

(11.22)

式中,θj是区域内定点i到线元Γj的张角(即关于动点j在线元始、末点所张开的角),d是rd的垂足到Γj起点的距离。

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

(11.24)式(11.23)和式(11.24)中的d、θj、sj、r1、r2、rd可由定点i坐标(xi,yi)和Γj两端点坐标(x1,y1)、(x2,y2)计算求得，即

3．Bi的计算

把区域D剖分成有限个三角形单元，利用高斯求积法计算各单元的积分，即

(11.25)需要说明的是,当泊松方程的边界条件全部为第一类或第二类边界条件时，可直接利用式(11.16)计算出边界上未知的q或u，进而根据式(11.15)计算出区域D内的u。当然，矩阵H、G的计算方法不变。

【例11.1】

设有泊松方程的混合边值问题(如图11.9所示)图11.9例11.1示意图

取每个边为边界元。图中①、②、③三点为三角形各边的中点。试写出边界积分方程并采用常数边界法计算图中点③的电势和①、②两点的电量。

解取每个边即为边界元。

(1)计算对角线元素Hii、Gii

。根据式(11.20)和(11.21),有

计算非对角线元素Hij、Gij(i≠j)。根据式(11.22)有

显然Hij≠Hji。因此有

而根据式(11.24)有G31=(1-0)ln

-1+0.5×1.107=-0.335,Gij≠Gji，经过计算可得

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

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

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

(4)对边界数据进行整理。

为下节编程需要，我们将最后求出的uq=(q1,q2,u3)中的u3和边界上已知的uq0=(u1,u2,q3)中的q3互换，即将(q1,q2,q3)存入uq，(u1,u2,u3)存入uq0。仍然考虑式(11.9)~(11.11)所满足的泊松方程定解问题，即

11.3单一边界下的边界元法计算程序对于单一边界条件下的泊松方程边界元计算步骤如下：由输入信息、计算矩阵、引入边界条件，形成方程到解边界元方程求边界上的u、q，再到求区域D内的u，最后输出计算结果。

主程序流程如下:

(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)为边界上各段端点坐标，(xm,ym)为边界上各段中点的坐标，(xd,yd)为内点坐标,!gama(j)=0(j∈Γ1),gama(j)=1(j∈Γ2);uq0=u0j(j∈Γ1),uq0=q0j(j∈Γ2)

(2)具体计算。

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)输出结果。

输出n个边界中点处的坐标及u，q值xm(j),ym(j),uq0(j),uq(j),j=1,2,…,n，输出内点坐标和u值xd,yd,ud

子程序流程如下：

(1)计算矩阵，形成方程。

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))

！计算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

当j∈Γ2时，将-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)!对应于形成的矩阵A

H(i,j)=-ch

55continue

50continue

do60i=1,n

uq(i)=0

do70j=1,n70uq(i)=uq(i)+H(i,j)*uq0(j)!对应于形成的矩阵R

60continue

return

end

！计算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)求区域D内的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)所示，区域1和2分别有介电常数为ε1和ε2的两种均匀介质。假设两区域内没有自由电荷，则电势u满足拉普拉斯方程，2u=0，在边界元方程中，对应的边界元方程为GQ=HU，如式(11.16)。11.4两种介质情况下的边界元方法当我们解边界元方程时，不仅要利用外边界Γ1和Γ2上的边界条件，还要利用内边界ΓI上的边界条件。为此，引入下列符号

图11.10区域1与区域2的划分

将两个区域的边界元方程写成如下形式:

(11.26)

(11.27)各矩阵的上标代表区域，下标代表边界。在两接着的交界处，电势u和电位移法向分量εr连续，所以

(11.28)

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

(11.30)这就是有两种介质的边界元方程，它已经包含了内边界(交界)条件。进一步要利用外边界条件，把已知量移到方程的右边，未知量移到左边,变AX=R的形式。这些过程与11.2节所讨论的相同。但应该注意编号规则。两个区域的边界元都按逆时针向编号，并且两个区域的编号连起来，如图11.10(b)所示。设边界Γ1、Γ2、ΓI上的边界元数分别为N1、N2、NI，则方程(11.30)中各个矩阵在总矩阵中的位置如表11.1所示，其中N=N1+N2+2NI。表11.1各个矩阵在总矩阵中的位置程序流程如下：

(1)输入数据:

subroutineinput(N,εr)

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

dataN1,N2,NI,εr/…/

N=N1+N2+2NI

两区边界元数，N<40边界元端点坐标和边界条件：

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

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

Γ1

Γ2

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

(2)计算矩阵，形成方程:

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)

两个区的始末编号

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和H

K1=J1(K),K2=J2(K)!两个区的始末编号

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，通过移项把方程(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)输出结果:

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)

调用subroutinegauss(N,G,uq)后，uq存有求出的u和q，而uq0中存有已知的u和q。需要进行整理，把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,的始末编号