工程電磁場(chǎng)數(shù)值計(jì)算4(有限元法1)剖析

工程電磁場(chǎng)數(shù)值計(jì)算4(有限元法1)剖析第一頁，共82頁。第4章電磁場(chǎng)有限元法

（FiniteElementMethod,FEM）

有限元法可以基于變分原理導(dǎo)出，也可以基于加權(quán)余量法導(dǎo)出，本章以加權(quán)余量法作為有限元法的基礎(chǔ)，以靜電場(chǎng)問題的求解為例介紹有限元法的基本原理與實(shí)施步驟。并介紹有限元法中的一些特殊問題。第二頁，共82頁。第4章電磁場(chǎng)有限元法（FEM）有限元基本原理與實(shí)施步驟：1DFEM有限元基本原理與實(shí)施步驟：2DFEM有限元方程組的求解二維有限元工程應(yīng)用三維有限元原理與工程應(yīng)用矢量有限元第三頁，共82頁。加權(quán)余量法回顧： 對(duì)算子方程用作為該方程的近似解（試探解）：代入方程得余量：1.有限元法基本原理與實(shí)施步驟：一維問題

在有限元法中，基函數(shù)一般用表示。采用Galerkin方案，取權(quán)函數(shù)與基函數(shù)相同。使與余量正交化：第四頁，共82頁。設(shè)L為線性算子，代入，得或記得代數(shù)方程組：加權(quán)余量法回顧（續(xù)）第五頁，共82頁。利用有限元法求解一維邊值問題：（1）單元剖分 如圖5個(gè)單元，6個(gè)節(jié)點(diǎn)（2）選取基函數(shù)

第六頁，共82頁。（3）方程離散（計(jì)算系數(shù)陣[K]和右端項(xiàng)[b]）

基函數(shù)Ni只是一階可導(dǎo) 的，不能嚴(yán)格滿足微分方 程，稱為“弱解”。第七頁，共82頁。（3）方程離散第一項(xiàng)在xj

處為0，在xi

處的值被來自(i-1)單元的貢獻(xiàn)抵消，故只剩下第二項(xiàng)。由于基函數(shù)Ni

局域支撐，顯見只有不為0。使用分步積分：第八頁，共82頁。（3）方程離散故類似，當(dāng)j=i

時(shí)右端項(xiàng)：第九頁，共82頁。總體方程強(qiáng)加邊界條件：u1=0,u6=0第十頁，共82頁。（4）求解方程思考：（1）有限元的解跟有限差分法的解有何根本不同？ （2）有限元的系數(shù)陣總是對(duì)稱的嗎？x0000.20.03610.03600.40.06280.06250.60.07100.07080.80.05250.05231.000第十一頁，共82頁。與有限差分法（FDM）相比，有限差分法是對(duì)點(diǎn)的離散，得到一系列離散點(diǎn)上的解；而有限元（FEM）是對(duì)區(qū)域的離散（單元），盡管所求的是節(jié)點(diǎn)上的自由度，但它的解在場(chǎng)域中每一個(gè)點(diǎn)上都有定義。所以，即是有限元節(jié)點(diǎn)上的解是精確的，有限元的整個(gè)解仍然是近似的。好的數(shù)據(jù)處理技術(shù)可以從該近似解中提取更精確的分析結(jié)果。線性單元中，如果所求的自由度是電位j，單元中的電場(chǎng)E是場(chǎng)量；節(jié)點(diǎn)上的E取鄰近單元的平均。一些補(bǔ)充說明：

關(guān)于有限元的解第十二頁，共82頁。計(jì)算系數(shù)陣是有限元分析的主要工作量。所涉及到的積分，如果不是解析可積的，通常要用到數(shù)值積分。其中最常用的數(shù)值積分方法是Gauss數(shù)值積分。一些補(bǔ)充說明：

高斯數(shù)值積分先將積分區(qū)間變換到[-1,1]上；按照固定的積分點(diǎn)計(jì)算若干函數(shù)值P(xi),以固定權(quán)值wi累加即可。具(2n+1)階精度。第十三頁，共82頁。n=4x(1)=0.8613d0x(2)=0.339981043584856d0w(1)=0.347854845137454d0w(2)=0.6526d0n=5x(1)=0.9664d0x(2)=0.538469310105683d0x(3)=0.0d0w(1)=0.236926885056189d0w(2)=0.478628670499366d0w(3)=0.568888888888889d0n=6x(1)=

0.932469514203152d0x(2)=

0.6612d0x(3)=

0.238619186083197d0w(1)=0.1770d0w(2)=0.3639d0w(3)=0.4679d0n=16x(1)=0.9894003948d0x(2)=0.9445750231d0x(3)=0.8656312024d0x(4)=0.7554044084d0x(5)=0.6178762444d0x(6)=0.4580167777d0x(7)=0.2816035508d0x(8)=0.0950125098d0w(1)=0.0271524594d0w(2)=0.0622535239d0w(3)=0.0951585117d0w(4)=0.1246289713d0w(5)=0.1495959888d0w(6)=0.1691565194d0w(7)=0.1826034150d0w(8)=0.1894506105d0一些Gauss積分點(diǎn)和權(quán)值：

（關(guān)于x=0對(duì)稱，只給出一半）第十四頁，共82頁。為提高有限元分析精度，有兩種方法： 其一：增加節(jié)點(diǎn)，細(xì)化網(wǎng)格——稱為h方法。 其二：增加有限元的階數(shù)——稱為p方法。一些補(bǔ)充說明：

線性單元與高階單元第十五頁，共82頁。一些補(bǔ)充說明：

二階單元第十六頁，共82頁。一些補(bǔ)充說明：

三階單元第十七頁，共82頁。h方法和p方法的求解精度ByJianmingJin.TheFiniteElementMethodinElectromagnetics,2ndEd.,2002第十八頁，共82頁。作業(yè)：要獨(dú)立完成，凡雷同者沒分！！

編寫有限元程序，計(jì)算一維邊值問題。改變剖分單元數(shù)目，觀察解的精度變化。（建議也同時(shí)做一個(gè)有限差分法的程序，比較二者的精度差別）第十九頁，共82頁。以二維靜電場(chǎng)泊松方程的求解為例。2.有限元法基本原理與實(shí)施步驟：

二維問題

目標(biāo)：依據(jù)加權(quán)余量法，利用分域基，建立離散的代數(shù)方程組，即確定系數(shù){Kij}

和{bi}。第二十頁，共82頁。

場(chǎng)域離散 二維問題常使用三角形單元離散，便于處理復(fù)雜的場(chǎng)域形狀，容易實(shí)現(xiàn)。

單元：互不重疊，覆蓋全部場(chǎng)域；每個(gè)單元內(nèi)介質(zhì)是 單一、均勻的。

節(jié)點(diǎn)：網(wǎng)格的交點(diǎn)，待求變量的設(shè)置點(diǎn)。

該步驟需要記錄的信息： 節(jié)點(diǎn)編號(hào)、節(jié)點(diǎn)坐標(biāo) 節(jié)點(diǎn)屬性(激勵(lì)源、是否邊界等)

單元編號(hào) 單元節(jié)點(diǎn)編號(hào) 單元介質(zhì) 第二十一頁，共82頁?；瘮?shù)

有限元采用分片逼近的思想，類似于一維情況下使用折線逼近一條任意曲線。使用分域基Ni，基函數(shù)的個(gè)數(shù)等于節(jié)點(diǎn)的個(gè)數(shù)；每個(gè)基函數(shù)Ni的作用區(qū)域是與該節(jié)點(diǎn)i相關(guān)聯(lián)的所有單元。 第二十二頁，共82頁。三角形單元內(nèi)的基函數(shù) 設(shè)三角形三個(gè)頂點(diǎn)處待求函數(shù)值分別為u1,u2,u3。如果單元足夠小，可以采用線性近似，將單元內(nèi)任意p點(diǎn)的u(x,y)表示為

代入三個(gè)頂點(diǎn)的坐標(biāo)和函數(shù)值，可以解出a、b、c。得到第二十三頁，共82頁。

單元節(jié)點(diǎn)的編號(hào)按逆時(shí)針方向排列！其中，第二十四頁，共82頁。記住我們的任務(wù)—尋找基函數(shù)對(duì)比可得基函數(shù)Ni常被稱為插值函數(shù)或者形狀函數(shù)，具有以下性質(zhì)：（1）是插值的；（2）（3）在相鄰單元的公共邊界上，

Ni是連續(xù)的，從而通過Ni構(gòu)造的逼近函數(shù)也是連續(xù)的。第二十五頁，共82頁。在積分中，對(duì)于確定的i，j的有效取值為i本身以及與節(jié)點(diǎn)i相聯(lián)的周圍節(jié)點(diǎn)，積分的有效區(qū)域?yàn)橐詉、j為公共節(jié)點(diǎn)的所有三角形單元，在這些單元中Ni、Nj才有交疊。計(jì)算系數(shù)陣

第二十六頁，共82頁。

這些積分可以分單元進(jìn)行。例如對(duì)右圖所示的局部編碼，K01、K00以及b0的計(jì)算公式為：

計(jì)算系數(shù)陣

第二十七頁，共82頁。以下把單元e的貢獻(xiàn)記為這樣，就有

每個(gè)或的計(jì)算都在具體的單元內(nèi)單獨(dú)考慮（稱為單元分析）。第二十八頁，共82頁。單元分析：計(jì)算單元內(nèi)積分對(duì)系數(shù)陣和右端項(xiàng)元素的貢獻(xiàn)。

系數(shù)陣元素：當(dāng)L為拉普拉斯算子時(shí)，由于Ni在單元內(nèi)是(x,y)的線性函數(shù)，經(jīng)Laplace算子作用后值為0。但是，在相鄰單元的邊界上，Ni是連續(xù)但是不光滑的，因此對(duì)積分的貢獻(xiàn)主要來自邊界。為考慮單元邊界的影響，需要借助于格林公式：第二十九頁，共82頁。故，格林公式：因：第三十頁，共82頁。

寫成一般形式，若一個(gè)三角形三個(gè)頂點(diǎn)編號(hào)為i,j,m（逆時(shí)針順序），則

從而第三十一頁，共82頁。再看邊界部分：

（1）在節(jié)點(diǎn)i的對(duì)邊Gjm上，Ni＝0，故積分貢獻(xiàn)為0；

結(jié)論：?jiǎn)卧吔鐚?duì)積分的貢獻(xiàn)為0。所以單元e對(duì)系數(shù)陣元素的貢獻(xiàn)為：

（2）在節(jié)點(diǎn)i

的鄰邊Gij上，由于計(jì)算Kij時(shí)需要把具有公共鄰邊的單元的積分累加，此二單元的Ni是連續(xù)的；對(duì)于單一均勻媒質(zhì)，要求相鄰單元滿足 ，故積分的貢獻(xiàn)相互抵消。第三十二頁，共82頁。

由于單元很小，做單元分析時(shí)通?？梢匀(e)為常數(shù)值（可以認(rèn)為等于三個(gè)頂點(diǎn)上的平均值）。因此

右端項(xiàng)元素：公式：第三十三頁，共82頁。上述以節(jié)點(diǎn)為序的分析過程對(duì)于有限元原理的說明是易于理解的。而在實(shí)際編程中，更有效率的是以單元為序，逐個(gè)計(jì)算單元系數(shù)陣[K(e)]，然后合成整體系數(shù)陣[K]。單元系數(shù)陣[K(e)]定義為

設(shè)i,j,m是節(jié)點(diǎn)的整體編號(hào)，元素Kij在整體矩陣中的實(shí)際位置是第i行、j列；因此必須合成到整體矩陣的第i行、j列元素上。

單元矩陣：第三十四頁，共82頁。

整體矩陣合成：第三十五頁，共82頁。

通過上述過程，對(duì)于一個(gè)“正?！钡膬?nèi)部節(jié)點(diǎn)就建立起了一個(gè)代數(shù)方程?！胺钦！钡墓?jié)點(diǎn)包括：媒質(zhì)交界面銜接條件和場(chǎng)域邊界條件。第三十六頁，共82頁。

對(duì)于靜電場(chǎng)問題，媒質(zhì)分界面銜接條件為媒質(zhì)交界面銜接條件

第一個(gè)條件是自動(dòng)滿足的（Why？），無須格外處理。

對(duì)于第二個(gè)條件，前面計(jì)算單元邊界上積分時(shí)，默認(rèn)兩邊u

的法向?qū)?shù)相等，使內(nèi)邊界上的積分結(jié)果抵消。因此只要把泊松方程寫成或 滿足的條件將是,從而也無需另行處理。第三十七頁，共82頁。

由于有限元方法能夠自動(dòng)滿足媒質(zhì)交界面條件，因此有限元法特別適合于處理多層復(fù)雜媒質(zhì)問題。這是其它方法無可比擬的。媒質(zhì)交界面銜接條件第三十八頁，共82頁。第一類邊界條件（強(qiáng)加邊界條件）第一類邊界節(jié)點(diǎn)是指邊界上函數(shù)值已知。因此處理方法是，合成整體系數(shù)陣之后，將該節(jié)點(diǎn)所在行的主元素置1，其它元素均置零，同時(shí)將右端項(xiàng)中對(duì)應(yīng)元素設(shè)為已知函數(shù)值。要保持對(duì)稱性；有更簡(jiǎn)便的做法第三十九頁，共82頁。第二類邊界條件（自然邊界條件）第二類邊界節(jié)點(diǎn)是指邊界上函數(shù)法向?qū)?shù)已知。對(duì)于內(nèi)部單元，相鄰單元邊界的積分相互抵消。但是對(duì)于場(chǎng)域邊界，如果給定第二類邊界條件不為0，則積分結(jié)果要計(jì)入右端項(xiàng)中。但是若給定的是齊次第二類邊界條件，則積分結(jié)果為0，無需另行處理，非常方便?！@是ANSYS中自動(dòng)滿足的邊界條件。第四十頁，共82頁。有限元方法的推導(dǎo)過程雖然看起來有些復(fù)雜，但是最終結(jié)果是非常簡(jiǎn)單而且優(yōu)美的。因?yàn)檫吔鐥l件的處理和媒質(zhì)交界面條件的處理都非常方便，使得有限元方法在處理復(fù)雜媒質(zhì)問題和復(fù)雜場(chǎng)域問題時(shí)得心應(yīng)手，獲得了廣泛的應(yīng)用，成為最重要的數(shù)值分析手段，廣泛應(yīng)用于各個(gè)領(lǐng)域。有人用“功蓋四方”來形容有限元，實(shí)不為過。中國人在有限元的發(fā)明中有自己獨(dú)特的貢獻(xiàn)。第四十一頁，共82頁。作業(yè)：(2)對(duì)于研究方向?yàn)閿?shù)值計(jì)算的同學(xué)： 編寫一個(gè)二維靜電場(chǎng)有限元程序，計(jì)算右圖所示問題，或其它自己找一個(gè)問題。（1）推導(dǎo)三角形單元的2次和3次插值函數(shù)。第四十二頁，共82頁。3.有限元方程組的求解代數(shù)方程組求解方法概述

所有的數(shù)值方法最終都?xì)w結(jié)為求解一個(gè)代數(shù)方程組：

系數(shù)陣A也稱系統(tǒng)矩陣或剛度矩陣。不同離散方法得到的系統(tǒng)矩陣具有不同的特點(diǎn)，方程組的解法也就不同。

基于微分方程（如FEM、FDM等）得到的系統(tǒng)矩陣是稀疏的，有時(shí)還是對(duì)稱的； 而基于積分方程得到的系統(tǒng)矩陣則是稠密的，如BEM、模擬電荷法等。第四十三頁，共82頁。

代數(shù)方程組的求解是數(shù)值計(jì)算（計(jì)算數(shù)學(xué)）研究的核心內(nèi)容。求解代數(shù)方程組的方法歸納起來有兩類：直接法和迭代法。 3.有限元方程組的求解直接法：直接法都是基于高斯消去法，經(jīng)過確定次數(shù)的運(yùn)算，理論上可以得到方程組的精確解。適用于小型、稠密方程組的計(jì)算。

迭代法：是一種間接方法，從某個(gè)預(yù)定的初值出發(fā)，按照一定的迭代步驟，逐漸逼近方程組的真解。得到一個(gè)滿足給定精度要求的近似解。適用于大型、稀疏方程組的計(jì)算。第四十四頁，共82頁。3.有限元方程組的求解直接法（LU分解算法）第四十五頁，共82頁。LU分解算法：回帶：消元：

計(jì)算量：需要的乘除法次數(shù)：

O(n3)

穩(wěn)定性：選主元第四十六頁，共82頁。迭代法迭代法的基本思想：（等價(jià)方程組）從一組猜測(cè)的初值開始迭代直至不再變化為止，即為方程組的解（收斂）。好的迭代法應(yīng)該：對(duì)初值不敏感；收斂速度快。第四十七頁，共82頁。例如：高斯—賽德爾迭代（有限差分法常用）第四十八頁，共82頁。高斯—賽德爾迭代實(shí)際運(yùn)算過程：

。。。。。這就是通常的高斯塞德爾迭代格式，矩陣中的零元素不參與運(yùn)算，矩陣甚至可以不出現(xiàn)。大大減少了內(nèi)存需求量和計(jì)算量。第四十九頁，共82頁。共軛梯度法（ConjugateGradientMethod,CG法） 共軛梯度法在原理上可以通過n步迭代得到方程的準(zhǔn)確解，因而也稱為半直接法或半迭代法。把迭代法表示為更一般的形式：

稱為步長，p稱為搜索方向，用r表示殘差。第五十頁，共82頁。預(yù)優(yōu)共軛梯度法（PreconditionedConjugateGradientMethod,PCG法）當(dāng)系數(shù)陣的特征值較均勻地分布在一個(gè)很長的區(qū)間上時(shí)，稱矩陣具有很大的條件數(shù)；此時(shí)共軛梯度法的收斂速度可能很慢。解決的辦法是選取合適的非奇異矩陣C進(jìn)行處理：矩陣稱為預(yù)處理矩陣。第五十一頁，共82頁。第五十二頁，共82頁。

預(yù)優(yōu)矩陣M應(yīng)具有如下特性：稀疏性與A相近；矩陣的特征值分布集中；形如的方程組容易求解；易于尋找。目前公認(rèn)有效的方法是對(duì)系數(shù)陣A做不完全Cholesky分解，以M=LDLT

為預(yù)優(yōu)矩陣。這種方法稱為不完全分解預(yù)優(yōu)共軛梯度法（IncompleteCholeskydecompositionpreconditionedConjugateGradientMethod，ICCG法）。ICCG法第五十三頁，共82頁。稀疏矩陣技術(shù)

沒有稀疏矩陣技術(shù)就沒有有限元的成功。

帶狀矩陣技術(shù)

通過合適的節(jié)點(diǎn)編碼，可使系統(tǒng)矩陣的非零元素集中于主對(duì)角線附近的帶形區(qū)域內(nèi)。在使用LU分解法求解方程組的過程中，帶形區(qū)域以外的0元素?zé)o需計(jì)算。

缺點(diǎn)：節(jié)點(diǎn)編碼優(yōu)化；帶形區(qū)域內(nèi)仍然有大量零元素第五十四頁，共82頁。

非零元素存儲(chǔ)技術(shù)

只存儲(chǔ)矩陣的非零元素。一般用一個(gè)一維數(shù)組存儲(chǔ)矩陣的非零元素，用另一個(gè)索引數(shù)組存儲(chǔ)這些非零元素在原矩陣中的行和列值。例如：第五十五頁，共82頁。

非零元素存儲(chǔ)技術(shù)

在迭代法中，系統(tǒng)矩陣參與的運(yùn)算只有矩陣左乘以某個(gè)相量。如果采用上面的存儲(chǔ)方法，則實(shí)現(xiàn)q=Ap的算法如下：

CComputethematrix–vectorproduct DOk=1,not i=IA(k) j=JA(k) q(i)=q(i)+A(k)*p(j) ENDDO第五十六頁，共82頁。

非零元素存儲(chǔ)技術(shù)

下列的二維存儲(chǔ)方案是我的發(fā)明：直接將矩陣擠扁。數(shù)組A有n行m列，m是A各行非零元素個(gè)數(shù)的最大值。對(duì)角元素放在該行的第一個(gè)位置上，且不管是否為0，都要存貯。用一個(gè)整型數(shù)組JA存放各元素的列號(hào)。JA中第一列可以用來記錄每行的非零元素個(gè)數(shù)。第五十七頁，共82頁。本節(jié)更多的參考文獻(xiàn)：金建銘.電磁場(chǎng)有限元方法，西安電子科技大學(xué)出版社，1998徐樹方，矩陣計(jì)算的理論與方法，北京大學(xué)出版社，1995楊紹祺，談根林，稀疏矩陣——算法及其程序?qū)崿F(xiàn),高等教育出版社,1985劉萬勛,劉長學(xué),華伯浩等，大型稀疏線性方程組的解法.國防工業(yè)出版社,1981RP梯華森,稀疏矩陣，科學(xué)出版社，1981第五十八頁，共82頁。

有限元分析精度的影響因素靜電場(chǎng)問題靜磁場(chǎng)問題渦流問題波的傳輸與散射4.二維有限元的工程應(yīng)用第五十九頁，共82頁。有限元分析精度的影響因素（1）數(shù)學(xué)模型對(duì)工程問題的近似；（2）材料電磁參數(shù)的不確定性；（3）數(shù)學(xué)模型的有限元近似：場(chǎng)域擬合精度——單元大小、未知數(shù)個(gè)數(shù)與局部 場(chǎng) 的變化情況；邊界擬合精度——曲線邊界；系數(shù)陣計(jì)算過程中數(shù)值積分精度；方程求解精度、數(shù)字誤差；計(jì)算結(jié)果的數(shù)據(jù)處理；好的處理技術(shù)可以提高分析 精度。第六十頁，共82頁。h方法、p方法；高次單元單元階數(shù)與計(jì)算精度的關(guān)系

ByJin。第六十一頁，共82頁。曲邊三角形單元：更好地?cái)M合曲面邊界第六十二頁，共82頁。數(shù)值積分——三角元的高斯數(shù)值積分 單元比較小的時(shí)候，單元內(nèi)的函數(shù)可以近似認(rèn)為是常數(shù)，通常可以獲得滿意的精度。當(dāng)單元內(nèi)函數(shù)變化比較快，或者采用曲邊三角形單元、高次單元時(shí)時(shí)，都會(huì)用到數(shù)值積分。

高斯-勒讓德數(shù)值積分公式：——L1、L2、L3

是位置的面積坐標(biāo)?！獧?quán)值wi如下頁表格。第六十三頁，共82頁。三角形單元高斯勒讓德數(shù)值積分點(diǎn)和權(quán)值。第六十四頁，共82頁。三角形單元高斯勒讓德數(shù)值積分點(diǎn)和權(quán)值。第六十五頁，共82頁。二維邊值問題的通用形式在媒質(zhì)交界面G12上，第六十六頁，共82頁。

單元矩陣元素計(jì)算