有限元和邊界元方法

計算物理3/lesson/ComputationalPhysics有限元和邊界元方法有限元和邊界元方法物理問題的變分原理泊松方程的有限元方法擴散方程的有限元方法波動方程的有限元方法邊界積分方程邊界元近似單一邊界下的邊界元法兩種介質的邊界元方法√物理問題的變分原理(1/3)有限元方法基于變分原理的離散化方法——部分逼近地離散化劃分整體區(qū)域為有限個基本塊(單元)在單元上插值逼近，得到結構簡單的函數集(有限元空間，是泛函

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)有限元方法的具體步驟劃分整體區(qū)域為有限個單元和編號劃分要點三角形的頂點相連避免鈍角(因引入較大誤差)每個三角形不跨越不同的介質每個三角形最多只有一條邊在G2

上(方便計算)三角形覆蓋盡量多的區(qū)域編號約定三角形單元的編號：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

的正方形區(qū)域G1G2G2G2xyO劃分整體區(qū)域為有限個單元和編號單元：①②③④頂點：123(23在G2上)結點：⑴⑵⑶⑷⑸(先內部和

G2

)③①②④121112223333⑴⑵⑶⑷⑸√泊松方程的有限元方法(8/11)構造線性插值函數G1G2G2G2xyO③①②④121112223333⑴⑵⑶⑷⑸√泊松方程的有限元方法(9/11)建立單元的矩陣建立頂點和結點的對應關系

V(e,i)G1G2G2G2xyO③①②④121112223333⑴⑵⑶⑷⑸√泊松方程的有限元方法(10/11)集成泛函和建立方程求解方程√泊松方程的有限元方法(11/11)例：如下圖的環(huán)形均勻帶電板，內徑

6，外徑

10，外圈

G1

的電勢為常數，內圈

G2

的電場為常數G1G2610考慮對稱性，取

1/4

環(huán)形區(qū)域以簡化計算G2G1①②③④⑤⑥⑦⑧⑨⑩123⑴⑵⑶⑷⑸⑹⑺⑻⑼⑽⑾⑿⒀⒁⒂√擴散方程的有限元方法(1/1)二維擴散方程離散化關于ai(t)

的常微分方程組初始條件求解方法：二級歐拉法√波動方程的有限元方法(1/1)二維波動方程離散化關于ai(t)

的常微分方程組初始條件求解方法：二級歐拉法√邊界積分方程(1/2)邊界元方法的特點基于邊界積分方程的近似方法結點僅分布在區(qū)域的邊界以邊界積分方程為控制方程，將邊界離散插值，轉化為代數方程組未知量的個數少，求解的計算量少面積分

邊界積分定理：如果

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

的計算：將區(qū)域

D劃分為有限個三角形單元例：三角區(qū)域的電勢和電量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)例：三角區(qū)域的電勢和電量主程序xOy11√單一邊界下的邊界元法(2/6)邊界離散化和初始化√單一邊界下的邊界元法(3/6)構造矩陣

H、G、A和

R√單一邊界下的邊界元法(4/6)非對角元Hij

和Gij√單一邊界下的邊界元法(5/6)求解線性方程√單一邊界下的邊界元法(6/6)計算內點作圖√兩種介質的邊界元方法(1/1)二維區(qū)域、無自由電荷的泊松方程G1G2GID1D2區(qū)域拆分和邊界元方程區(qū)域拆分