有限元和边界元方法

计算物理3/lesson/ComputationalPhysics有限元和边界元方法有限元和边界元方法物理问题的变分原理泊松方程的有限元方法扩散方程的有限元方法波动方程的有限元方法边界积分方程边界元近似单一边界下的边界元法两种介质的边界元方法√物理问题的变分原理(1/3)有限元方法基于变分原理的离散化方法——部分逼近地离散化划分整体区域为有限个基本块(单元)在单元上插值逼近，得到结构简单的函数集(有限元空间，是泛函

J(y)

的定义域的子集)将边值问题转化为泛函的极值问题在有限元空间中寻找泛函

J(y)

的极小值，作为近似解物理中的变分例：力学体系的最小作用量原理体系的特性可以用拉格朗日函数L(q,

q',

t)

描写在时刻t1

和t2

之间体系按照以下积分取最小值的方式运动(即，运动轨道由泛函的的极小值决定)√物理问题的变分原理(2/3)例：光学的费马原理光从点A到点B的传播路径是使光程L

取极值由d

L=

0

得到几何光学的折射定律和反射定律例：电磁学的麦克斯韦方程组电磁场的拉格朗日函数L

是空间积分电磁学的作用量是时间积分运动方程由泛函的的极小值决定(即

d

S

=

0

)√物理问题的变分原理(3/3)例：静电场的泊松方程第一类边界条件等价的变分问题为求解泛函的极值问题泛函的求解必须在边界条件下：条件变分问题第二类和第三类边界条件等价的变分问题为求解泛函的极值问题边界条件包含泛函中：自然边界条件√泊松方程的有限元方法(1/11)静电场中二维泊松方程的有限元方法G1G2ABDG1阳极G2G1阴极G2例：阴极射线管(如右图)

，在两极上(边界G1)的电势u

是已知的，在左右两侧(边界G2)的q

是已知的。如果管中的自由电荷密度分布r(x,y)

已知，则以上的泊松方程等价为求解以下泛函J(u)

的极值问题√泊松方程的有限元方法(2/11)有限元方法的具体步骤划分整体区域为有限个单元和编号划分要点三角形的顶点相连避免钝角(因引入较大误差)每个三角形不跨越不同的介质每个三角形最多只有一条边在G2

上(方便计算)三角形覆盖尽量多的区域编号约定三角形单元的编号：e

=

①,②,③,…顶点的编号：按逆时针为

1,

2,

3顶点的坐标：(x1,

y1),

(x2,

y2),

(x3,

y3)单元的泛函：Je(u)③①②e123e整体的泛函：√泊松方程的有限元方法(3/11)3(x3,y3)e1(x1,y1)2(x2,y2)构造线性插值函数假设每个单元内

u(x,y)

是

x和

y的线性函数每个的三个基函数u(x,y)

的插值表达式中，a,

b,

c,

d

可由三角形的顶点坐标确定，只剩余u1,

u2,

u3

未知√泊松方程的有限元方法(4/11)建立单元的矩阵单元的泛函123eG2第一项积分与单元刚度矩阵

(zij)第二项积分与单元矩阵

(rfj)第三项积分与单元矩阵

(rqj)√泊松方程的有限元方法(5/11)建立顶点和结点的(V-n)对应关系单元编号：有一条边在

G2

上且

q

0

的单元编号为

1,

2,…,

e1，其余的单元编号为e1+1,

e1+2,…,e0G2G1③①②e1e1+1e0顶点编号：用

V(e,i)

表示，逆时针方向，2和

3在

G2

上结点编号：内部和

G2

上的结点编号为

1,

2,…,

n1，G1

上的结点编号为

n1+1,

n1+2,…,

n0建立顶点和结点的对应关系：V(e,i)

=n集成泛函和建立方程泛函的离散化K为总体刚度矩阵，由单元刚度矩阵

(zij)

合成Rf

由单元矩阵

(rfj)

合成，Rq

由单元矩阵

(rqj)

合成J(u)

被离散化为二次多元函数J(u1,

u2,…,

un0)√泊松方程的有限元方法(6/11)有限元方程(关于

um

的线性方程组)求解方程√泊松方程的有限元方法(7/11)例：如右图的边长为

1

的正方形区域G1G2G2G2xyO划分整体区域为有限个单元和编号单元：①②③④顶点：123(23在G2上)结点：⑴⑵⑶⑷⑸(先内部和

G2

)③①②④121112223333⑴⑵⑶⑷⑸√泊松方程的有限元方法(8/11)构造线性插值函数G1G2G2G2xyO③①②④121112223333⑴⑵⑶⑷⑸√泊松方程的有限元方法(9/11)建立单元的矩阵建立顶点和结点的对应关系

V(e,i)G1G2G2G2xyO③①②④121112223333⑴⑵⑶⑷⑸√泊松方程的有限元方法(10/11)集成泛函和建立方程求解方程√泊松方程的有限元方法(11/11)例：如下图的环形均匀带电板，内径

6，外径

10，外圈

G1

的电势为常数，内圈

G2

的电场为常数G1G2610考虑对称性，取

1/4

环形区域以简化计算G2G1①②③④⑤⑥⑦⑧⑨⑩123⑴⑵⑶⑷⑸⑹⑺⑻⑼⑽⑾⑿⒀⒁⒂√扩散方程的有限元方法(1/1)二维扩散方程离散化关于ai(t)

的常微分方程组初始条件求解方法：二级欧拉法√波动方程的有限元方法(1/1)二维波动方程离散化关于ai(t)

的常微分方程组初始条件求解方法：二级欧拉法√边界积分方程(1/2)边界元方法的特点基于边界积分方程的近似方法结点仅分布在区域的边界以边界积分方程为控制方程，将边界离散插值，转化为代数方程组未知量的个数少，求解的计算量少面积分

边界积分定理：如果

u(x,

y)

和

v(x,

y)

是定义在平面域

D上的两个任意函数，并且它们在边界外法线上的导数为xyOdsds√证明：利用格林公式(略)边界积分方程(2/2)d

函数定义：性质：平面域：边界积分方程v和p

是已知的，上式给出D+

G

内任意一点(等式左边)与G

上某点(等式右边)的u

和q

之间的关系上式称为边界积分方程，是边界元法的基础√边界元近似(1/4)边界元方程当定点

i在边界G

上常数边界离散化边界积分方程的右边：将边界G

分成N

段，以每段中点的u

和q

近似该段的函数值例：泊松方程当定点

i不在边界G

上：利用前面的结果，代入最上公式√边界元近似(2/4)对角元Hii

和Gii

的计算(定点

i在

Gi

上)非对角元Hij

和Gij

的计算(定点

i不在

Gj

上)√边界元近似(3/4)Bi

的计算：将区域

D划分为有限个三角形单元例：三角区域的电势和电量xOy11泊松方程的混合边界问题常数边界离散化对角元Hii

和GiiA①G2G1G1②③非对角元Hij

和Gij①②③①②③A①G2G1G1②③A①G2G1G1②③A①G2G1G1②③A①G2G1G1②③√Bi

：因为f

=

0，所以Bi

=

0边界元近似(4/4)构造方程AX

=

R内部(三角形中心(1/3,

0))：u

=

0.623√单一边界下的边界元法(1/6)例：三角区域的电势和电量主程序xOy11√单一边界下的边界元法(2/6)边界离散化和初始化√单一边界下的边界元法(3/6)构造矩阵

H、G、A和

R√单一边界下的边界元法(4/6)非对角元Hij

和Gij√单一边界下的边界元法(5/6)求解线性方程√单一边界下的边界元法(6/6)计算内点作图√两种介质的边界元方法(1/1)二维区域、无自由电荷的泊松方程G1G2GID1D2区域拆分和边界元方程区域拆分