有限元和邊界元方法

1、計(jì)算物理計(jì)算物理3/lesson/ComputationalPhysics3/lesson/ComputationalPhysics有限元和邊界元方法有限元和邊界元方法有限元和邊界元方法有限元和邊界元方法n物理問題的變分原理物理問題的變分原理n泊松方程的有限元方法泊松方程的有限元方法n擴(kuò)散方程的有限元方法擴(kuò)散方程的有限元方法n波動(dòng)方程的有限元方法波動(dòng)方程的有限元方法n邊界積分方程邊界積分方程n邊界元近似邊界元近似n單一邊界下的邊界元法單一邊界下的邊界元法n兩種介質(zhì)的邊界元方法兩種介質(zhì)的邊界元方法物理問題的變分原理物理問

2、題的變分原理(1/(1/3 3) )n有限元方法有限元方法n基于變分原理的離散化方法基于變分原理的離散化方法部分逼近地離散化部分逼近地離散化n劃分整體區(qū)域?yàn)橛邢迋€(gè)基本塊劃分整體區(qū)域?yàn)橛邢迋€(gè)基本塊( (單元單元) )n在單元上插值逼近，得到結(jié)構(gòu)簡單的函數(shù)集在單元上插值逼近，得到結(jié)構(gòu)簡單的函數(shù)集( (有限元空有限元空間，是泛函間，是泛函 J( (y) ) 的定義域的子集的定義域的子集) )n將邊值問題轉(zhuǎn)化為泛函的極值問題將邊值問題轉(zhuǎn)化為泛函的極值問題n在有限元空間中尋找泛函在有限元空間中尋找泛函 J( (y) ) 的極小值，作為近似解的極小值，作為近似解n物理中的變分物理中的變分n例：力學(xué)體系的最

3、小作用量原理例：力學(xué)體系的最小作用量原理n體系的特性可以用拉格朗日體系的特性可以用拉格朗日函數(shù)函數(shù) L( (q, , q, , t) ) 描寫描寫n在在時(shí)刻時(shí)刻 t1 1 和和 t2 2 之間之間體系按照以下積分取最小值的方式體系按照以下積分取最小值的方式運(yùn)動(dòng)運(yùn)動(dòng)( (即，運(yùn)動(dòng)軌道由泛函的的極小值決定即，運(yùn)動(dòng)軌道由泛函的的極小值決定) )得到運(yùn)動(dòng)方程，由0d),(21 SttqqLStt物理問題的變分原理物理問題的變分原理(2/(2/3 3) )n例：光學(xué)的費(fèi)馬原理例：光學(xué)的費(fèi)馬原理n光從點(diǎn)光從點(diǎn) A 到點(diǎn)到點(diǎn) B 的傳播路徑是使光程的傳播路徑是使光程 L 取極值取極值B A dsnLn由由

4、L 0 0 得到幾何光學(xué)的折射定律和反射定律得到幾何光學(xué)的折射定律和反射定律n例：電磁學(xué)的麥克斯韋方程組例：電磁學(xué)的麥克斯韋方程組n電磁場的電磁場的拉格朗日拉格朗日函數(shù)函數(shù) L 是空間積分是空間積分：體積元：標(biāo)勢和磁矢勢：電荷和電流密度，：電場和磁場強(qiáng)度導(dǎo)率，：介質(zhì)的介電常數(shù)和磁d , ,d )(21 22AujHE AjHELn電磁學(xué)的作用量是時(shí)間積分電磁學(xué)的作用量是時(shí)間積分21 dtttLS 運(yùn)動(dòng)方程由泛函的的極小值決定運(yùn)動(dòng)方程由泛函的的極小值決定( (即即 S 0 0 ) )物理問題的變分原理物理問題的變分原理(3/(3/3 3) )n例：靜電場的泊松方程例：靜電場的泊松方程n第一類邊界

5、條件第一類邊界條件),(),( ),(02222yxuyxuyxfyuxu 等價(jià)的變分問題為求解泛函的極值問題等價(jià)的變分問題為求解泛函的極值問題),(),( ,dd)()(21)(022yxuyxuyxf uyuxuuJ 泛函的求解必須在邊界條件下：條件變分問題泛函的求解必須在邊界條件下：條件變分問題n第二類和第三類邊界條件第二類和第三類邊界條件)( ,2222 unufyuxu 等價(jià)的變分問題為求解泛函的極值問題等價(jià)的變分問題為求解泛函的極值問題s u uyxf uyuxuuJd)21(dd)()(21)(222 邊界條件包含泛函中：自然邊界條件邊界條件包含泛函中：自然邊界條件泊松方程的有限

6、元方法泊松方程的有限元方法(1/(1/1111) )n靜電場中二維泊松方程的靜電場中二維泊松方程的有限元方法有限元方法2102222, ),(, ),(),(D, ),(21yxyxqnuyxyxuyxuyxyxfyuxu 1 1 2 2ABD 1 1 陽極陽極 2 2 1 1 陰極陰極 2 2n例：例：陰極射線管陰極射線管( (如右圖如右圖) ) ，在兩極上，在兩極上( (邊邊界界 1 1) )的電勢的電勢 u 是已知的，在左右兩側(cè)是已知的，在左右兩側(cè)( (邊界邊界 2 2) )的的 q 是已知的。如果管中的自由是已知的。如果管中的自由電荷密度分布電荷密度分布 ( (x, , y) ) 已知

7、，則已知，則qnuuuyxfyxu21 , ),(),(002 以上的泊松方程以上的泊松方程等價(jià)為求解以下泛函等價(jià)為求解以下泛函 J( (u) ) 的極值問題的極值問題0D212 ,ddd )(21)(uusq uyxf uuuJ泊松方程的有限元方法泊松方程的有限元方法(2/(2/1111) )n有限元方法的具體步驟有限元方法的具體步驟n劃分整體區(qū)域?yàn)橛邢迋€(gè)單元和編號(hào)劃分整體區(qū)域?yàn)橛邢迋€(gè)單元和編號(hào)n劃分要點(diǎn)劃分要點(diǎn)n三角形的頂點(diǎn)相連三角形的頂點(diǎn)相連n避免鈍角避免鈍角( (因引入較大誤差因引入較大誤差) )n每個(gè)三角形不跨越不同的介質(zhì)每個(gè)三角形不跨越不同的介質(zhì)n每個(gè)三角形最多只有一條邊每個(gè)三角形

8、最多只有一條邊在在 2 2 上上( (方便計(jì)算方便計(jì)算) )n三角形覆蓋盡量多的區(qū)域三角形覆蓋盡量多的區(qū)域n編號(hào)約定編號(hào)約定n三角形單元的編號(hào)：三角形單元的編號(hào)：e , , , , , ,n頂點(diǎn)頂點(diǎn)的的編號(hào)：按逆時(shí)針為編號(hào)：按逆時(shí)針為 1, 1, 2, 2, 3 3n頂點(diǎn)頂點(diǎn)的坐標(biāo)：的坐標(biāo)：( (x1 1, , y1 1), ), ( (x2 2, , y2 2), ), ( (x3 3, , y3 3) )n單元的單元的泛函：泛函：Je( (u) )e123en整體的整體的泛函：泛函：1)()(eeuJuJ泊松方程的有限元方法泊松方程的有限元方法(3/(3/1111) )3 (x3, y3)

9、e1 (x1, y1)2 (x2, y2)n構(gòu)造線性插值函數(shù)構(gòu)造線性插值函數(shù)n假設(shè)每個(gè)單元內(nèi)假設(shè)每個(gè)單元內(nèi) u( (x, , y) ) 是是 x 和和 y 的的線性線性函數(shù)函數(shù)n每個(gè)每個(gè)的的三個(gè)基函數(shù)三個(gè)基函數(shù)nu( (x, , y) ) 的的插值表達(dá)式中，插值表達(dá)式中，a, , b, , c, , d 可由三角形的頂點(diǎn)坐可由三角形的頂點(diǎn)坐標(biāo)確定，只標(biāo)確定，只剩余剩余 u1 1, , u2 2, , u3 3 未知未知泊松方程的有限元方法泊松方程的有限元方法(4/(4/1111) )n建立單元的矩陣建立單元的矩陣n單元的泛函單元的泛函123e2eeeesq uf uuuJ2ddd )(21)(

10、2n第一項(xiàng)積分與單元?jiǎng)偠染仃嚨谝豁?xiàng)積分與單元?jiǎng)偠染仃?( (z i j) )n第二項(xiàng)積分與單元矩陣第二項(xiàng)積分與單元矩陣 ( (r f j) )n第三項(xiàng)積分與單元矩陣第三項(xiàng)積分與單元矩陣 ( (r q j) )泊松方程的有限元方法泊松方程的有限元方法(5/(5/1111) )n建立頂點(diǎn)和結(jié)點(diǎn)的建立頂點(diǎn)和結(jié)點(diǎn)的( (V n) )對應(yīng)關(guān)系對應(yīng)關(guān)系n單元編號(hào)：有一條邊在單元編號(hào)：有一條邊在 2 2 上且上且 q 0 0 的單元編號(hào)為的單元編號(hào)為 1, 1, 2, 2, , , e1 1，其，其余的單元編號(hào)余的單元編號(hào)為為 e1 11,1, e1 12,2, , e0 0 2 2 1 1e1 1e1 1

11、11e0 0n頂點(diǎn)編號(hào)：用頂點(diǎn)編號(hào)：用 V( (e, , i) ) 表示，逆時(shí)針方向，表示，逆時(shí)針方向，2 2和和 3 3在在 2 2 上上n結(jié)點(diǎn)編號(hào)：內(nèi)部和結(jié)點(diǎn)編號(hào)：內(nèi)部和 2 2 上的結(jié)點(diǎn)編號(hào)為上的結(jié)點(diǎn)編號(hào)為 1, 1, 2, 2, , , n1 1， 1 1 上的結(jié)點(diǎn)編號(hào)為上的結(jié)點(diǎn)編號(hào)為 n1 11,1, n1 12,2, , , n0 0n建立頂點(diǎn)和結(jié)點(diǎn)的對應(yīng)關(guān)系：建立頂點(diǎn)和結(jié)點(diǎn)的對應(yīng)關(guān)系：V( (e, , i) ) nn集成泛函和建立方程集成泛函和建立方程n泛函的離散化泛函的離散化nK 為總體剛度矩陣，由單元?jiǎng)偠染仃嚍榭傮w剛度矩陣，由單元?jiǎng)偠染仃?( (z i j) ) 合成合成nR

12、f 由由單元矩陣單元矩陣 ( (r f j) ) 合成，合成，Rq 由由單元矩陣單元矩陣 ( (r q j) ) 合成合成nJ( (u) ) 被離散化為二次多元函數(shù)被離散化為二次多元函數(shù) J( (u1 1, , u2 2, , , , un0 0) )泊松方程的有限元方法泊松方程的有限元方法(6/(6/1111) )n有限元方程有限元方程( (關(guān)于關(guān)于 um 的線性方程組的線性方程組) )n求解方程求解方程泊松方程的有限元方法泊松方程的有限元方法(7/(7/1111) )n例：如右圖的邊長為例：如右圖的邊長為 1 1 的正方形區(qū)域的正方形區(qū)域 1 1 2 2 2 2 2 2xyOyxnuuux

13、yyuxu2 , 5 . 0) 1 , 0( , 4 . 0) 1 , 1 (2222n劃分整體區(qū)域?yàn)橛邢迋€(gè)單元和編號(hào)劃分整體區(qū)域?yàn)橛邢迋€(gè)單元和編號(hào)n單元：單元：n頂點(diǎn)：頂點(diǎn)：1 2 31 2 3( (2 32 3在在 2 2上上) )n結(jié)點(diǎn)：結(jié)點(diǎn)：( (先先內(nèi)部和內(nèi)部和 2 2 ) )1 12 21 11 11 12 22 22 23 33 33 33 3泊松方程的有限元方法泊松方程的有限元方法(8/(8/1111) )n構(gòu)造線性插值函數(shù)構(gòu)造線性插值函數(shù) 1 1 2 2 2 2 2 2xyO1 12 21 11 11 12 22 22 23 33 33 33 3泊松方程的有限元方法泊松方程的

14、有限元方法(9/(9/1111) )n建立單元的矩陣建立單元的矩陣n建立頂點(diǎn)和結(jié)點(diǎn)的對應(yīng)關(guān)系建立頂點(diǎn)和結(jié)點(diǎn)的對應(yīng)關(guān)系 V ( (e, , i) ) 1 1 2 2 2 2 2 2xyO1 12 21 11 11 12 22 22 23 33 33 33 3泊松方程的有限元方法泊松方程的有限元方法(10/(10/1111) )n集成泛函和建立方程集成泛函和建立方程n求解方程求解方程泊松方程的有限元方法泊松方程的有限元方法(11/(11/1111) )n例：如下圖的環(huán)形均勻帶電板，內(nèi)徑例：如下圖的環(huán)形均勻帶電板，內(nèi)徑 6 6，外徑，外徑 1010，外圈外圈 1 1 的電勢為常數(shù)，內(nèi)圈的電勢為常數(shù)，

15、內(nèi)圈 2 2 的電場為常數(shù)的電場為常數(shù)12 ,100 , 4212222nuuyuxu 1 1 2 26 61010n考慮對稱性，取考慮對稱性，取 1/41/4 環(huán)形環(huán)形區(qū)域以簡化計(jì)算區(qū)域以簡化計(jì)算 2 2 1 11 12 23 3擴(kuò)散方程的有限元方法擴(kuò)散方程的有限元方法(1/(1/1 1) )n二維二維擴(kuò)散方程擴(kuò)散方程0 ),()0 ,( ),(02uyxuyxuyxfuDtun離散化離散化ijjjiiiiiDiDiDiyxyxttyxufuDtu),( , , ),()(),(d dd為基函數(shù)n關(guān)于關(guān)于 i( (t) ) 的常微分方程組的常微分方程組DiiDjiijDjiijf fDKMf

16、KtMd ,d ,dddn初始條件初始條件DiuMd ,)0(0n求解方法：二級(jí)歐拉法求解方法：二級(jí)歐拉法波動(dòng)方程的有限元方法波動(dòng)方程的有限元方法(1/(1/1 1) )n二維二維波動(dòng)方程波動(dòng)方程0 ),( ),( ),(000222uyxRtuyxuuyxfuDtuttn離散化離散化ijjjiiiiiDiDiDiyxyxttyxufuDtu),( , , ),()(),(d dd22為基函數(shù)n關(guān)于關(guān)于 i( (t) ) 的常微分方程組的常微分方程組DiiDjiijDjiijf fDKMfKtMd ,d ,ddd22n初始條件初始條件DiDiRrurtMMd ,d ,d)0(d ,)0(0n求

17、解方法：二級(jí)歐拉法求解方法：二級(jí)歐拉法邊界積分方程邊界積分方程(1/(1/2 2) )n邊界元方法的特點(diǎn)邊界元方法的特點(diǎn)n基于基于邊界積分方程的近似方法邊界積分方程的近似方法n結(jié)點(diǎn)僅分布在區(qū)域的邊界結(jié)點(diǎn)僅分布在區(qū)域的邊界n以邊界積分方程為控制方程，將邊界離散插值，轉(zhuǎn)化以邊界積分方程為控制方程，將邊界離散插值，轉(zhuǎn)化為代數(shù)方程組為代數(shù)方程組n未知量的個(gè)數(shù)少，求解的計(jì)算量少未知量的個(gè)數(shù)少，求解的計(jì)算量少n面積分面積分邊界積分邊界積分n定理：如果定理：如果 u( (x, , y) ) 和和 v( (x, , y) ) 是定義在平面域是定義在平面域 D 上的兩個(gè)上的兩個(gè)任意函數(shù)，并且它們在邊界外法線上的

18、導(dǎo)數(shù)為任意函數(shù)，并且它們在邊界外法線上的導(dǎo)數(shù)為xyOddsnpnvqnu ,supvqvuuvd)(d)( D22那么n證明：利用格林公式證明：利用格林公式( (略略) )邊界積分方程邊界積分方程(2/(2/2 2) )n 函數(shù)函數(shù)n定義：定義：1d)( ,0, 10, 0)(xxxxxn性質(zhì)：性質(zhì)：)0(d)()(fxxxfn平面域：平面域：的距離到定點(diǎn)為動(dòng)點(diǎn) P ),(ln212irrrn邊界積分方程邊界積分方程nv 和和 p 是已知的，上式給出是已知的，上式給出 D 內(nèi)任意一點(diǎn)內(nèi)任意一點(diǎn)( (等等式左邊式左邊) )與與 上某點(diǎn)上某點(diǎn)( (等式右邊等式右邊) )的的 u 和和 q 之間的之

19、間的關(guān)系關(guān)系n上式稱為邊界積分方程，是邊界元法的基礎(chǔ)上式稱為邊界積分方程，是邊界元法的基礎(chǔ)邊界元近似邊界元近似(1/(1/4 4) )n邊界元方程邊界元方程n當(dāng)定點(diǎn)當(dāng)定點(diǎn) i 在在邊界邊界 上上n常數(shù)邊界離散化常數(shù)邊界離散化n邊界積分方程的右邊：將邊界邊界積分方程的右邊：將邊界 分成分成 N 段，以每段中點(diǎn)段，以每段中點(diǎn)的的 u 和和 q 近似該段的函數(shù)值近似該段的函數(shù)值n例：泊松方程例：泊松方程n當(dāng)定點(diǎn)當(dāng)定點(diǎn) i 不在不在邊界邊界 上：利用前面的結(jié)果，代入最上公式上：利用前面的結(jié)果，代入最上公式邊界元近似邊界元近似(2/(2/4 4) )n對角對角元元 Hii 和和 Gii 的的計(jì)算計(jì)算(

20、(定點(diǎn)定點(diǎn) i 在在 i 上上) )n非對角非對角元元 Hij 和和 Gij 的的計(jì)算計(jì)算( (定點(diǎn)定點(diǎn) i 不不在在 j 上上) )邊界元近似邊界元近似(3/(3/4 4) )nBi 的的計(jì)算：將區(qū)域計(jì)算：將區(qū)域 D 劃分為有限個(gè)三角形單元?jiǎng)澐譃橛邢迋€(gè)三角形單元eeirfrfBdlndlnDn例：三角區(qū)域的電勢和電量例：三角區(qū)域的電勢和電量0 ,2 , 02左邊斜邊、底邊nuyuuxOy11 泊松方程的混合邊界問題泊松方程的混合邊界問題n常數(shù)邊界常數(shù)邊界離散化離散化n對角對角元元 Hii 和和 GiiA 2 2 1 1 1 1n非對角非對角元元 Hij 和和 GijA 2 2 1 1 1 1A 2 2 1 1 1 1A 2 2 1 1 1 1A 2 2 1 1 1 1nBi ：因?yàn)椋阂驗(yàn)?f 0 0，所以，所以 Bi 0 0邊界元近似邊界元近似(4/(4/4 4) )n構(gòu)造構(gòu)造方程方程 AX Rn內(nèi)部內(nèi)部( (三角形中心三角形中心(1/3,(1/3, 0)0) )：u 0.6230.623單一邊界下的邊界元法單一邊界下的邊界元法(1/(1/6 6) )n例：三角區(qū)域的電勢和電量例：三角區(qū)域的電勢和電量0 ,2 , 02左邊斜邊、底邊nuyuun主程序主程序xOy11單一邊界下的邊界元法單一邊界下的邊界元法(2/(2/6 6) )n邊界邊界離