《工程電磁場原理》課程講稿(有限元)

1、二維靜態(tài)電磁場的有限元方法(FEM)簡介靜電場： 特定條件才有解析解穩(wěn)恒磁場： 特定條件才有解析解，不適合工程應(yīng)用。第一節(jié) 加權(quán)余量法(Weighted Residual Method)-伽遼金法（Galerikin Mehtod）一 加權(quán)余量法(Weighted Residual Method) 有一邊值問題方程： 算符，對函數(shù)u的運(yùn)算，f是已知函數(shù)，求解u。為了求解u，有一系列線性無關(guān)函數(shù)u1,u2,ui，也叫基（序列）函數(shù)。 取前m項(xiàng)近似求u-即u的線性組合。 （當(dāng)m, ）則余量差（誤差）： 精確解：R=0;但是，如果在誤差允許范圍內(nèi)，滿足需要即可。滿足強(qiáng)制余量的加權(quán)（weight）積分

2、為零: Wi叫權(quán)函數(shù)序列，亦線性無關(guān)。二 伽遼金法（Galerikin Mehtod）-最常用方法若取權(quán)函數(shù)與基函數(shù)相等， 這種方法叫做伽遼金方程（Galerikin Mehtod）加權(quán)余量方法。第二節(jié) 有限元方法的基本思想1有限元方法的基本思想首先，將一個(gè)閉合場域進(jìn)行有限元剖分，也就是把一個(gè)閉合場域劃分為N個(gè)微小的有限單元（簡稱有限元或單元），即其次，在每個(gè)單元上構(gòu)造插值函數(shù)逼近真解，將待求函數(shù)用各單元上的表示為在單元上，進(jìn)一步地將用插值函數(shù)和節(jié)點(diǎn)待求函數(shù)值表示為其中，i為單元上節(jié)點(diǎn)序號(hào)，r為單元的總的節(jié)點(diǎn)數(shù)。第三，求各個(gè)單元上的加權(quán)余量方程，并將各個(gè)單元上的加權(quán)余量方程相加獲得代數(shù)方程組（

3、或?qū)⒚總€(gè)單元插值合成的總插值函數(shù)代泛定方程的等價(jià)泛函并求極值獲得代數(shù)方程組）。第四，求解代數(shù)方程組即得場域中的各節(jié)點(diǎn)函數(shù)值，從而完成函數(shù)的數(shù)值求解。進(jìn)一步求解其他相關(guān)問題。下面，對二維靜態(tài)電磁場的有限元方法進(jìn)行介紹。第三節(jié) 單元剖分與插值函數(shù)1單元剖分在單元剖分過程中，一般應(yīng)該遵守如下幾條規(guī)則：（1）場域是一個(gè)封閉區(qū)域（對于開區(qū)域問題，需特殊處理）；（2）單元不能跨越邊界或介質(zhì)交界；（3）單元上的點(diǎn)不能落在相鄰單元的邊或面上，只能與相鄰單元的點(diǎn)重合；（4）各個(gè)單元不能共交；（5）全部單元應(yīng)充滿整個(gè)場域。圖1電機(jī)的單元剖分圖。圖1 電機(jī)定子與轉(zhuǎn)子結(jié)構(gòu)2單純形單元（1）一維空間設(shè)是一維空間待求的解

4、函數(shù)，和分別為一維問題中場域的邊界點(diǎn)。首先，將場域剖分為如圖2所示的N個(gè)線段形的單元（e=1, 2, , N），即這些線段形單元的長短可以不同，在一維空間中，最簡單的單元形狀是線段，因此線段又被稱為一維空間中的單純形。在每個(gè)單元上構(gòu)造插值函數(shù)，并用單元兩個(gè)端點(diǎn)和上的節(jié)點(diǎn)函數(shù)插值和來逼近單元上的解，即 在上圖2 一維場域的單元剖分和線性逼近其次，將所有單元上的近似解合成，構(gòu)成整個(gè)場域的近似解，即 在上可以看出，在有限元方法中一個(gè)重要的問題就是插值函數(shù)的選取。通常將插值函數(shù)選為x的函數(shù)，即 在上其中稱為單元上節(jié)點(diǎn)i的插值函數(shù)。因此在單元上近似解可以表示為 在上顯然，在各單元之間的節(jié)點(diǎn)上應(yīng)該是連續(xù)的

5、。從上式可以推斷出，插值函數(shù)應(yīng)該具有下面的特性，當(dāng)時(shí)，時(shí)，（2）二維空間設(shè)是二維問題中場域上的解函數(shù)，場域的邊界為。首先，將場域剖分為如圖3所示的N個(gè)三角形的單元（e=1, 2, , N），即圖3 二維場域的單元剖分和線性逼近這些單元的形狀可以是多樣的，既可以是三角形也可以是四邊形或其他形狀等。在二維空間中，最簡單的單元形狀是三角形，因此三角形又被稱為二維空間中的單純形（三維空間的單純形為四面體）。本講義只討論單元為三角形的情況。將三角形單元上的三個(gè)節(jié)點(diǎn)函數(shù)值選為待求變量，則在單元上解函數(shù)可以被三個(gè)節(jié)點(diǎn)函數(shù)插值表示為 在上其中i, j, k分別單元上的三個(gè)節(jié)點(diǎn)編號(hào)，稱為三角形單元上節(jié)點(diǎn)i的插值

6、函數(shù)。其次，將各個(gè)單元上的近似解合成在一起，即得場域上的近似解函數(shù)為 在上顯然，和在各單元之間的節(jié)點(diǎn)上應(yīng)該是連續(xù)的。如以標(biāo)示三角形單元上節(jié)點(diǎn)的坐標(biāo)，由單元上的近似式，推得當(dāng)時(shí)， ，當(dāng)時(shí)，當(dāng)時(shí)，而在單元分析時(shí)，一個(gè)關(guān)鍵的問題是插值函數(shù)的選取及相應(yīng)的計(jì)算公式。下面，分節(jié)進(jìn)行詳細(xì)討論。3插值函數(shù)以上討論中引入的插值函數(shù)又常被稱為形狀函數(shù)，即稱為單元上節(jié)點(diǎn)的形狀函數(shù)。由以上討論知，它具有如下性質(zhì)不妨設(shè)，即，為單元上的節(jié)點(diǎn)總數(shù)，則有（歸一化特性）（1）一維空間在一維空間上，一個(gè)線段構(gòu)成一維空間的單純形。設(shè)單元落在區(qū)間，如圖4所示。圖4 一維形狀函數(shù)顯然，有寫成矩陣形式有由此，可以求出形狀函數(shù)和。易知顯然

7、，上述行列式就是線段單元的長度，即進(jìn)一步，得顯然，形狀函數(shù)就是與相應(yīng)節(jié)點(diǎn)線段的長度之比。當(dāng)時(shí)， ，當(dāng)時(shí)， ，在有限元方法中，常常用到如下形狀函數(shù)的積分運(yùn)算，即為計(jì)算該積分方便，考慮如下問題。由于，且，這就形成了與的線性變換。而可以被線性表示，即。顯然對于的積分不如對的積分具有普遍性，因?yàn)楹笳叩姆e分的上下限分別為0和1。其微元變換為當(dāng)時(shí)，；當(dāng)時(shí)，。因此，有為計(jì)算上式積分，討論如下典型積分所以利用上述公式，得（2）二維空間在二維空間上，三角形構(gòu)成單純形。設(shè)單元落在圖5所示的三角形上。則有圖5 二維形狀函數(shù)上式即確定了單元上的三個(gè)形狀函數(shù)，即、和。為求形狀函數(shù)方便，將上述方程組寫為矩陣形式，有為使行

8、列式取正值，本講義設(shè)三角形的三個(gè)節(jié)點(diǎn)編號(hào)1、2、3依次按逆時(shí)針順序編號(hào)。其行列式為可見，上式行列式為單純形單元面積的兩倍。下面，求形狀函數(shù)得從上式可見，單元上各節(jié)點(diǎn)的形狀函數(shù)為該節(jié)點(diǎn)對應(yīng)的三角形面積與單元三角形面積之比。故形狀函數(shù)又稱為面積坐標(biāo)，又由于它在0到1的區(qū)間內(nèi)取值，也稱為自然坐標(biāo)。今后，為編制程序方便，將節(jié)點(diǎn)1 置換為，將節(jié)點(diǎn)2置換為，將節(jié)點(diǎn)3置換為，則可以求出單元上節(jié)點(diǎn)的形狀函數(shù)為 最后一項(xiàng)三角形剖分式中其中，下標(biāo)按模3相加。在有限元方法中，常常用到以下積分計(jì)算，即為計(jì)算積分方便，將坐標(biāo)系變換到坐標(biāo)系，如圖6所示。因?yàn)?，所以，其微元變換為（雅可比變換-Jacobi Transfor

9、mation）圖6 二維坐標(biāo)變換積分計(jì)算如下所以，得第四節(jié) 單元分析與單元合成1二維靜態(tài)電磁場的泊松方程現(xiàn)在，對如下泊松方程邊值問題進(jìn)行有限元分析，即 對于靜電場問題，為標(biāo)量電位，為介電常數(shù)，為自由電荷體密度。其邊值問題為對于恒定磁場問題，為矢量磁位的z軸分量A，為磁阻率，為自由電流密度。其邊值問題為 2單元分析首先，選擇單元上的每個(gè)節(jié)點(diǎn)的形狀函數(shù)Nie為權(quán)函數(shù)ui（Galerkin Method），這里 , ,fj為待定常數(shù)。然后在單元上對泊松方程進(jìn)行加權(quán)積分，將泊松方程在場域內(nèi)任意一點(diǎn)滿足放松為在單元的加權(quán)積分下滿足。例如對于單元上的節(jié)點(diǎn)i而言，代入Galerikin 方程 又由格林定理（

10、Green Theorem）(p333,13) 做替換， 并代入 將上式整理為 在單元上，將待解函數(shù)用單元上的節(jié)點(diǎn)形狀函數(shù)和節(jié)點(diǎn)待解函數(shù)值近似為（加權(quán)將上式代入前面的方程，得 ， i=1，2，3（三角形三個(gè)頂點(diǎn)）上式被稱為單元上的有限元方程。3單元合成在上述單元有限元方程中，可以看出存在一個(gè)單元邊界積分項(xiàng)，即 顯然，單元的邊界無外乎下列三種情況。（1）為單元與其他單元的交界面。此時(shí)可以看出，上述恰為靜電場邊界條件和恒定磁場邊界條件中的相關(guān)項(xiàng)，又考慮到相鄰兩個(gè)單元交界面的法線方向相反，因此，如果將兩個(gè)相鄰單元的有限元方程相加，并利用邊界條件，可知相鄰單元的邊界積分將相互抵消，即不出現(xiàn)在相加后的兩

11、個(gè)有限元方程中。以此類推，將所有單元的有限元方程相加，則在有限元方程中僅僅剩下場域外邊界上的邊界積分。（2）單元邊界落在第二類邊界上。此時(shí)有由于在第二類邊界上給定，為已知量，直接進(jìn)行積分即可。為記述方便，引入如下符號(hào) i, j=1,2,3可以看出(互易reciprocal)。單元上的有限元方程就可以寫為如寫成矩陣方程為由于，所以上述有限元方程的矩陣為對稱矩陣?？梢钥闯?，上述討論主要是圍繞單元進(jìn)行的，所以我們又將這個(gè)過程稱為單元分析過程。將全部單元的有限元方程進(jìn)行合成，得總的有限元方程為將單元求和與節(jié)點(diǎn)求和進(jìn)行換序，并整理得 i, j=1,2,3上式給出的有限元方程可以簡寫為或?qū)憺榫仃囆问狡渲杏?/p>

12、于，所以，即S為對稱矩陣。（3）單元邊界落在第一類邊界上, 若第k個(gè)點(diǎn)是第一類邊界條件的結(jié)點(diǎn)，。此時(shí)為未知量，需做特殊處理. 總系數(shù)矩陣S和右端列向量G，F(xiàn)做如下處理：(i) 令;(ii) 余量做替換 (i=1,2,3n)。(iii) ;(iv)具體實(shí)施方法：在有限元方程中，對應(yīng)于第一類邊界的邊界積分為未知量，需做特殊處理。由于它是未知的，因此第一類邊界上節(jié)點(diǎn)對應(yīng)的有限元方程的右端項(xiàng)中也是未知的。此時(shí)，這個(gè)節(jié)點(diǎn)對應(yīng)的方程可以用該節(jié)點(diǎn)給定的值來替代，即。在有限元計(jì)算程序中，通常采用如下的強(qiáng)加方法處理。設(shè)原有限元方程為改為其中，m為一個(gè)確保（可以忽略其它項(xiàng)）遠(yuǎn)遠(yuǎn)大于矩陣元素的整數(shù)，以確保方程求解后

13、使得（第i個(gè)方程變?yōu)?。這種方法被稱為第一類邊界條件強(qiáng)加方法。對所有邊界點(diǎn)進(jìn)行強(qiáng)加處理后，再代入已完成的矩陣方程，即可進(jìn)行求解。可以看出，上述矩陣和列向量中的各個(gè)元素是在總矩陣和總列向量中通過各個(gè)單元上節(jié)點(diǎn)對應(yīng)的元素疊加獲得的。所以，上式疊加過程又被稱為單元合成過程。為了清楚計(jì)算過程，結(jié)合圖7對單元合成時(shí)的有限元方程的形成過程進(jìn)行討論。單元、和的單元方程矩陣方程和在單元合成后在總矩陣方程的位置分別如下。圖7 單元和節(jié)點(diǎn)編號(hào)（55 matrices）單元： 單元： 單元： 三個(gè)單元合成后總的矩陣方程為推廣到n個(gè)節(jié)點(diǎn)(n/3)個(gè)三角形網(wǎng)格5二維靜態(tài)電磁場的單元分析現(xiàn)在，以二維空間為例作進(jìn)一步的討論。選三角形單元為剖分（meshing）單元，則形狀函數(shù)/權(quán)函數(shù)為單元矩陣的元素計(jì)算公式為 單元矩陣方程的右端列向量為按(i-1,i,i+1)次序綜上，得（1）二維靜電場有限元方程。對于靜電場問題，。其有限元方程為 求解上式方程即得各個(gè)節(jié)點(diǎn)電位。各個(gè)單元上的電場強(qiáng)度為 （2）二維恒定磁場有限元方程。對于靜電場問，。其有限元方程為 求解上式方程即得各個(gè)節(jié)點(diǎn)的矢量磁位的縱向分