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

1、2021/2/111工程電磁場數(shù)值計算（4）（電磁場有限元法）2021/2/112第4章 電磁場有限元法（Finite Element Method, FEM） 有限元法可以基于變分原理導出,也可以基于加權余量法導出,本章以加權余量法作為有限元法的基礎,以靜電場問題的求解為例介紹有限元法的基本原理與實施步驟。并介紹有限元法中的一些特殊問題。2021/2/113第4章 電磁場有限元法（FEM）1.有限元基本原理與實施步驟:1D FEM2.有限元基本原理與實施步驟:2D FEM3.有限元方程組的求解4.二維有限元工程應用5.三維有限元原理與工程應用6.矢量有限元2021/2/114加權余量法回顧:

2、對算子方程用 作為該方程的近似解（試探解）:代入方程得余量:1. 有限元法基本原理與實施步驟:一維問題( )L uf1niiiu u( )RL uf 在有限元法中，基函數(shù)一般用 表示。采用Galerkin方案，取權函數(shù)與基函數(shù)相同。使與余量正交化：(, ) ( ) d0iiN RN L uf (1, 2, , )in, 1, 2, , iNin2021/2/115設L為線性算子，代入 ，得11 () d() d0nnijjijjjjN LNfNL Nf 1niiiuN或(1, 2, , )in1() ddnjijijN L NN f 記diibN f,() di jijKN L N得代數(shù)方程組

3、: Kb加權余量法回顧（續(xù)）(, ) ( ) d0iiN RN L uf 2021/2/116利用有限元法求解一維邊值問題:（1）單元剖分 如圖5個單元,6個節(jié)點（2）選取基函數(shù) 22d( ) 01d(0)(1)0uL uuxxxuu 111111 (, ) ( , )iiiiiiiiiiixxxxxxxNxxxxxxx1niiiuu N2021/2/117（3）方程離散（計算系數(shù)陣 K 和右端項 b） 基函數(shù) Ni 只是一階可導的,不能嚴格滿足微分方程,稱為“弱解”。diibN f,() di jijKN L N Kb,2222() dd(+) ddd d ddi jijjijjiijKN

4、L NNNNxNNN Nx 2021/2/1182222d ddddddddddddjijjiijixjixxxjjiixxNNxNNxxNNNNxxxx（3）方程離散第一項在 xj 處為0,在 xi 處的值被來自 (i-1) 單元的貢獻抵消,故只剩下第二項。由于基函數(shù) Ni 局域支撐,顯見只有 不為0。使用分步積分:(1)ji ,1,1, , i ii ii iKKKdiibN f,() di jijKN L N2021/2/119（3）方程離散故diibN f,() di jijKN L N2,2d d ddddddddjjiiji jiijxxjiijxxNKNN NxNNxN Nxxx

5、 (1)ji 類似,當 j = i 時1111,ddddddiiiixxiii iiixxNNKxN Nxxx 11ddiixiiixbN fN f x 右端項:2021/2/1110總體方程11121121222322323334334344454454555655656666KKubKKKubKKKubKKKubKKKubKKub強加邊界條件:u1 = 0, u6 = 01221222323323334344344454554555656100010uuKKKbuKKKbuKKKbuKKKbu2021/2/1111（4）求解方程思考:（1）有限元的解跟有限差分法的解有何根本不同?（2）有限

6、元的系數(shù)陣總是對稱的嗎?ux0000.20.03610.03600.40.06280.06250.60.07100.07080.80.05250.05231.000*u2021/2/111222d( ) 01d(0)(1)0uL uuxxxuu 與有限差分法（FDM）相比,有限差分法是對點的離散,得到一系列離散點上的解;而有限元（FEM）是對區(qū)域的離散（單元）,盡管所求的是節(jié)點上的自由度,但它的解在場域中每一個點上都有定義。 所以,即是有限元節(jié)點上的解是精確的,有限元的整個解仍然是近似的。好的數(shù)據(jù)處理技術可以從該近似解中提取更精確的分析結果。 線性單元中,如果所求的自由度是電位j,單元中的電場

7、 E是場量;節(jié)點上的 E 取鄰近單元的平均。一些補充說明:關于有限元的解2021/2/1113 計算系數(shù)陣是有限元分析的主要工作量。所涉及到的積分,如果不是解析可積的,通常要用到數(shù)值積分。其中最常用的數(shù)值積分方法是Gauss數(shù)值積分。一些補充說明: 高斯數(shù)值積分diibN f,() di jijKN L N111( )d( )niiiP xxw P x先將積分區(qū)間變換到-1,1上;按照固定的積分點計算若干函數(shù)值 P(xi), 以固定權值 wi 累加即可。具(2n+1)階精度。2021/2/1114n=4 x(1)= 0.8613d0 x(2)= 0.339981043584856d0 w(1)

8、= 0.347854845137454d0 w(2)= 0.6526d0n=5 x(1)= 0.9664d0 x(2)= 0.538469310105683d0 x(3)= 0.0d0 w(1)= 0.236926885056189d0 w(2)= 0.478628670499366d0 w(3)= 0.568888888888889d0n=6 x(1)= 0.932469514203152d0 x(2)= 0.6612d0 x(3)= 0.238619186083197d0 w(1)= 0.1770d0 w(2)= 0.3639d0 w(3)= 0.4679d0n=16 x(1) = 0.9

9、894003948d0 x(2) = 0.9445750231d0 x(3) = 0.8656312024d0 x(4) = 0.7554044084d0 x(5) = 0.6178762444d0 x(6) = 0.4580167777d0 x(7) = 0.2816035508d0 x(8) = 0.0950125098d0 w(1) = 0.0271524594d0 w(2) = 0.0622535239d0 w(3) = 0.0951585117d0 w(4) = 0.1246289713d0 w(5) = 0.1495959888d0 w(6) = 0.1691565194d0 w(

10、7) = 0.1826034150d0 w(8) = 0.1894506105d0一些Gauss積分點和權值:（關于x=0對稱,只給出一半）2021/2/111522d( ) 01d(0)(1)0uL uuxxxuu 為提高有限元分析精度,有兩種方法:其一:增加節(jié)點,細化網(wǎng)格稱為h方法。其二:增加有限元的階數(shù)稱為p方法。一些補充說明: 線性單元與高階單元2021/2/111622d( ) 01d(0)(1)0uL uuxxxuu 一些補充說明: 二階單元2021/2/111722d( ) 01d(0)(1)0uL uuxxxuu 一些補充說明: 三階單元2021/2/1118h方法和p方法的求

11、解精度By Jianming Jin. The Finite Element Method in Electromagnetics, 2nd Ed., 20022021/2/1119作業(yè):要獨立完成,凡雷同者沒分!22d( ) 01d(0)(1)0uL uuxxxuu 編寫有限元程序,計算一維邊值問題。改變剖分單元數(shù)目,觀察解的精度變化。（建議也同時做一個有限差分法的程序,比較二者的精度差別）2021/2/1120 以二維靜電場泊松方程的求解為例。2. 有限元法基本原理與實施步驟:二維問題目標:依據(jù)加權余量法,利用分域基,建立離散的代數(shù)方程組,即確定系數(shù)Kij 和bi。diibN f,() d

12、i jijKN L N122222( )uuL ufxyuguhn2021/2/1121 場域離散二維問題常使用三角形單元離散,便于處理復雜的場域形狀,容易實現(xiàn)。單元:互不重疊,覆蓋全部場域;每個單元內介質是 單一、均勻的。節(jié)點:網(wǎng)格的交點,待求變量的設置點。該步驟需要記錄的信息:節(jié)點編號、節(jié)點坐標節(jié)點屬性(激勵源、是否邊界等)單元編號單元節(jié)點編號單元介質2021/2/1122基函數(shù) 有限元采用分片逼近的思想,類似于一維情況下使用折線逼近一條任意曲線。 使用分域基Ni,基函數(shù)的個數(shù)等于節(jié)點的個數(shù);每個基函數(shù)Ni的作用區(qū)域是與該節(jié)點i相關聯(lián)的所有單元。2021/2/1123 三角形單元內的基函數(shù)

13、設三角形三個頂點處待求函數(shù)值分別為u1, u2, u3。如果單元足夠小,可以采用線性近似,將單元內任意p點的u(x,y)表示為 ( , )u x yabxcy312123( , )( , )( , )( , )x yx yx yu x yuuu代入三個頂點的坐標和函數(shù)值,可以解出a、b、c。得到2021/2/1124312123( , )( , )( , )( , )x yx yx yu x yuuu12312311112xxxyyy 1232311112xxxyyy 2131311112xxxyyy 3121211112xxxyyy 單元節(jié)點的編號按逆時針方向排列!其中,2021/2/112

14、5112233( , )u x yNNN記住我們的任務尋找基函數(shù)對比312123( , )( , )( , )( , )x yx yx yu x yuuu( 1, 2, 3 )i 可得( , )iix yN基函數(shù)Ni常被稱為插值函數(shù)或者形狀函數(shù)，具有以下性質：（1）是插值的；（2）（3）在相鄰單元的公共邊界上， Ni是連續(xù)的，從而通過Ni構造的逼近函數(shù)也是連續(xù)的。1 ()(,)0 ()ijjijN xyij2021/2/1126在積分 中,對于確定的 i,j的有效取值為i 本身以及與節(jié)點i相聯(lián)的周圍節(jié)點,積分的有效區(qū)域為以i、j 為公共節(jié)點的所有三角形單元 ,在這些單元中Ni、Nj才有交疊。(

15、)dijijKN L N計算系數(shù)陣 diibN f,() di jijKN L N2021/2/1127這些積分可以分單元進行。例如對右圖所示的局部編碼,K01、K00以及b0的計算公式為:1234560000()dKN L N 160101()dKN L N 12345600dbN f 計算系數(shù)陣 diibN f,() di jijKN L N2021/2/1128以下把單元e的貢獻記為( )( )( )()deeeeijijKNL N( )( )( )deeeeiibNf這樣，就有(1)(2)(3)(4)(5)(6)00000000000000KKKKKKK(1)(6)010101KKK(

16、1)(2)(3)(4)(5)(6)0000000bbbbbbb 每個 或 的計算都在具體的單元內單獨考慮（稱為單元分析）。( )eijK( )eib2021/2/1129單元分析:計算單元內積分對系數(shù)陣和右端項元素的貢獻。( )( )( )()deeeeijijKNL N系數(shù)陣元素: 當L為拉普拉斯算子時,由于Ni在單元內是(x, y)的線性函數(shù),經Laplace算子作用后值為0。 但是,在相鄰單元的邊界上, Ni是連續(xù)但是不光滑的,因此對積分的貢獻主要來自邊界。 為考慮單元邊界的影響,需要借助于格林公式:2021/2/1130故 ,( )2()d dd ddeeejeijijijiNKNNx

17、 yNNx yNn 2( )d dVSVS格林公式:1233223321()()()2x yx yyy xxxy 23321()()()2yyxxNij31132()()()2yyxxNij23313213122()()()()()()4yyyyxxxxNN( , )iix yN因：2021/2/1131寫成一般形式,若一個三角形三個頂點編號為i, j, m（逆時針順序）,)()()()()()4yyyyxxxxNN d d()()()()4eijimjmimjmNNx yyyyyxxxx從而2()()()()()()4ijimjmimjmNNyyyyxxxx 202

18、1/2/1132dejiNNn再看邊界部分:（1）在節(jié)點 i 的對邊jm上,Ni0,故積分貢獻為0;( )()()()()d d4eimjmimjmeijijyyyyxxxxKNNx y 結論:單元邊界對積分的貢獻為0。所以單元e對系數(shù)陣元素的貢獻為:（2）在節(jié)點 i 的鄰邊ij上，由于計算Kij時需要把具有公共鄰邊的單元的積分累加，此二單元的Ni是連續(xù)的；對于單一均勻媒質，要求相鄰單元滿足 ，故積分的貢獻相互抵消。12()()/eeijijunun 2021/2/1133由于單元很小,做單元分析時通?？梢匀?f (e) 為常數(shù)值（可以認為等于三個頂點上的平均值）。因此( )( )( )dee

19、eeiibNf右端項元素:123! ! !() () () d2(2)!elmnl m nNNNlmn 公式：( )( )( )( )d3eeeeeiibfNf 2021/2/1134 上述以節(jié)點為序的分析過程對于有限元原理的說明是易于理解的。而在實際編程中,更有效率的是以單元為序,逐個計算單元系數(shù)陣K(e),然后合成整體系數(shù)陣K。單元系數(shù)陣K(e)定義為設 i, j, m 是節(jié)點的整體編號,元素Kij在整體矩陣中的實際位置是第i行、j列;因此 必須合成到整體矩陣的第i行、j列元素上。( )( )( )( )( )( )( )( )( )( )eeeiiijimeeeejijjjmeeemim

20、jmmKKKKKKKKKK( )eijK單元矩陣:2021/2/1135整體矩陣合成:iiijimiijijjmjjjmimjmmmmKKKubKKKubKKKub( )( )( )( )( )( )( )( )( )( )eeeiiijimeeeejijjjmeeemimjmmKKKKKKKKKK( )( )( )( )d3eeeeeiibfNf 2021/2/1136通過上述過程,對于一個“正?！钡膬炔抗?jié)點就建立起了一個代數(shù)方程?！胺钦！钡墓?jié)點包括:媒質交界面銜接條件和場域邊界條件。2021/2/1137對于靜電場問題,媒質分界面銜接條件為媒質交界面銜接條件12uu1212uunndej

21、iNNn第一個條件是自動滿足的（Why?）,無須格外處理。20/ru 0/ru 1212rruunn對于第二個條件，前面計算單元邊界上積分 時，默認兩邊 u 的法向導數(shù)相等，使內邊界上的積分結果抵消。因此只要把泊松方程寫成 或滿足的條件將是 , 從而也無需另行處理。2021/2/1138由于有限元方法能夠自動滿足媒質交界面條件,因此有限元法特別適合于處理多層復雜媒質問題。這是其它方法無可比擬的。媒質交界面銜接條件2021/2/1139第一類邊界條件（強加邊界條件）第一類邊界節(jié)點是指邊界上函數(shù)值 已知。因此處理方法是，合成整體系數(shù)陣之后，將該節(jié)點所在行的主元素置1，其它元素均置零，同時將右端項中

22、對應元素設為已知函數(shù)值。iiug要保持對稱性;有更簡便的做法2021/2/1140第二類邊界條件（自然邊界條件）第二類邊界節(jié)點是指邊界上函數(shù)法向導數(shù) 已知。對于內部單元,相鄰單元邊界的積分 相互抵消。但是對于場域邊界,如果給定第二類邊界條件不為0,則積分結果要計入右端項中。但是若給定的是齊次第二類邊界條件,則積分結果為0,無需另行處理,非常方便。/jjunh dejiNNn這是ANSYS中自動滿足的邊界條件。2021/2/1141 有限元方法的推導過程雖然看起來有些復雜,但是最終結果是非常簡單而且優(yōu)美的。因為邊界條件的處理和媒質交界面條件的處理都非常方便,使得有限元方法在處理復雜媒質問題和復雜

23、場域問題時得心應手,獲得了廣泛的應用,成為最重要的數(shù)值分析手段,廣泛應用于各個領域。有人用“功蓋四方”來形容有限元,實不為過。 中國人在有限元的發(fā)明中有自己獨特的貢獻。2021/2/1142作業(yè):(2) 對于研究方向為數(shù)值計算的同學:編寫一個二維靜電場有限元程序,計算右圖所示問題,或其它自己找一個問題。（1）推導三角形單元的2次和3次插值函數(shù)。2021/2/11433. 有限元方程組的求解代數(shù)方程組求解方法概述所有的數(shù)值方法最終都歸結為求解一個代數(shù)方程組:Axb系數(shù)陣 A也稱系統(tǒng)矩陣或剛度矩陣。不同離散方法得到的系統(tǒng)矩陣具有不同的特點,方程組的解法也就不同?；谖⒎址匠蹋ㄈ鏔EM、FDM等）得

24、到的系統(tǒng)矩陣是稀疏的,有時還是對稱的;而基于積分方程得到的系統(tǒng)矩陣則是稠密的,如BEM、模擬電荷法等。2021/2/1144代數(shù)方程組的求解是數(shù)值計算（計算數(shù)學）研究的核心內容。求解代數(shù)方程組的方法歸納起來有兩類:直接法和迭代法。3. 有限元方程組的求解 直接法:直接法都是基于高斯消去法,經過確定次數(shù)的運算,理論上可以得到方程組的精確解。適用于小型、稠密方程組的計算。 迭代法:是一種間接方法,從某個預定的初值出發(fā),按照一定的迭代步驟,逐漸逼近方程組的真解。得到一個滿足給定精度要求的近似解。適用于大型、稀疏方程組的計算。2021/2/11453. 有限元方程組的求解直接法（LU分解算法）2021

25、/2/1146LU分解算法:回帶:消元:計算量:需要的乘除法次數(shù): O(n3)穩(wěn)定性:選主元2021/2/1147迭代法(1)( )kkxA xb迭代法的基本思想:AxbxA xb（等價方程組）從一組猜測的初值 開始迭代(0)x直至 不再變化為止,即為方程組的解（收斂）。( )kx好的迭代法應該:對初值不敏感;收斂速度快。2021/2/1148112233445566410100141010014100100410010141001014bbbbbb 例如:高斯賽德爾迭代（有限差分法常用）11122233344455566640101004101010401010041000104010101

26、4001010bbbbbb(1)( )kkxA xb2021/2/114911122233344455566601010010101001010011000104010101001010bbbbbb(1)( )kkxA xb高斯賽德爾迭代實際運算過程:24114b135224b 。這就是通常的高斯塞德爾迭代格式,矩陣中的零元素不參與運算,矩陣甚至可以不出現(xiàn)。大大減少了內存需求量和計算量。2021/2/1150共軛梯度法（Conjugate Gradient Method, CG法）共軛梯度法在原理上可以通過n步迭代得到方程的準確解,因而也稱為半直接法或半迭代法。把迭代法表示為更一般的形式:(1

27、)( )( )( )kkkkxxp 稱為步長，p 稱為搜索方向，用r 表示殘差。2021/2/1151預優(yōu)共軛梯度法（Preconditioned Conjugate Gradient Method, PCG法）當系數(shù)陣的特征值較均勻地分布在一個很長的區(qū)間上時,稱矩陣具有很大的條件數(shù);此時共軛梯度法的收斂速度可能很慢。解決的辦法是選取合適的非奇異矩陣C進行處理: 11CAxC bAxbA yb1TT1CACC xC bTC xyTMCC矩陣 稱為預處理矩陣。2021/2/11522021/2/1153預優(yōu)矩陣M應具有如下特性:稀疏性與A相近;矩陣 的特征值分布集中;形如 的方程組容易求解;易于

28、尋找。目前公認有效的方法是對系數(shù)陣A做不完全Cholesky分解,以M = LDLT 為預優(yōu)矩陣。這種方法稱為不完全分解預優(yōu)共軛不完全分解預優(yōu)共軛梯度法梯度法（Incomplete Cholesky decomposition preconditioned Conjugate Gradient Method,ICCG法）。AM1rMy ICCG法2021/2/1154稀疏矩陣技術沒有稀疏矩陣技術就沒有有限元的成功。 帶狀矩陣技術通過合適的節(jié)點編碼,可使系統(tǒng)矩陣的非零元素集中于主對角線附近的帶形區(qū)域內。在使用LU分解法求解方程組的過程中,帶形區(qū)域以外的0元素無需計算。缺點:節(jié)點編碼優(yōu)化;帶形區(qū)域

29、內仍然有大量零元素2021/2/1155 非零元素存儲技術只存儲矩陣的非零元素。一般用一個一維數(shù)組存儲矩陣的非零元素,用另一個索引數(shù)組存儲這些非零元素在原矩陣中的行和列值。例如:1.02.000003.004.00005.006.007.008.009.00010.011.0A1.02.03.04.05.06.07.08.09.010.011.0IA1 1223344555JA122435241452021/2/1156 非零元素存儲技術在迭代法中,系統(tǒng)矩陣參與的運算只有矩陣左乘以某個相量。如果采用上面的存儲方法,則實現(xiàn) q = A p 的算法如下:C Compute the matrixve

30、ctor product DO k=1, not i=IA(k) j=JA(k) q(i)=q(i)+A(k)*p(j) ENDDO2021/2/1157 非零元素存儲技術下列的二維存儲方案是我的發(fā)明:直接將矩陣擠扁。數(shù)組A有n行m列,m是A各行非零元素個數(shù)的最大值。對角元素放在該行的第一個位置上,且不管是否為0,都要存貯。用一個整型數(shù)組JA存放各元素的列號。JA中第一列可以用來記錄每行的非零元素個數(shù)。 2021/2/1158 本節(jié)更多的參考文獻:(1) 金建銘. 電磁場有限元方法,西安電子科技大學出版社,1998(2) 徐樹方,矩陣計算的理論與方法,北京大學出版社,1995 (3) 楊紹祺,

31、談根林,稀疏矩陣算法及其程序實 現(xiàn) ,高等教育出版社, 1985(4) 劉萬勛, 劉長學, 華伯浩等,大型稀疏線性方程組的解法. 國防工業(yè)出版社, 1981(5) R P 梯華森, 稀疏矩陣,科學出版社,19812021/2/1159 有限元分析精度的影響因素 靜電場問題 靜磁場問題 渦流問題 波的傳輸與散射4. 二維有限元的工程應用2021/2/1160有限元分析精度的影響因素（1）數(shù)學模型對工程問題的近似;（2）材料電磁參數(shù)的不確定性;（3）數(shù)學模型的有限元近似: 場域擬合精度單元大小、未知數(shù)個數(shù)與局部場的變化情況; 邊界擬合精度曲線邊界; 系數(shù)陣計算過程中數(shù)值積分精度; 方程求解精度、數(shù)

32、字誤差; 計算結果的數(shù)據(jù)處理;好的處理技術可以提高分析精度。2021/2/1161h方法、p方法;高次單元單元階數(shù)與計算精度的關系By Jin。2021/2/1162曲邊三角形單元:更好地擬合曲面邊界2021/2/1163數(shù)值積分三角元的高斯數(shù)值積分單元比較小的時候,單元內的函數(shù)可以近似認為是常數(shù),通?？梢垣@得滿意的精度。當單元內函數(shù)變化比較快,或者采用曲邊三角形單元、高次單元時時,都會用到數(shù)值積分。高斯-勒讓德數(shù)值積分公式:( )( )( )()deeeeijijKNL N( )( )( )deeeeiibNf1231231(,)d d(,)emeeeeeeeiiiiiF L L Lx yw

33、F LLL L1、L2、L3 是位置的面積坐標。 權值 wi 如下頁表格。2021/2/1164三角形單元高斯勒讓德數(shù)值積分點和權值。2021/2/1165三角形單元高斯勒讓德數(shù)值積分點和權值。2021/2/1166二維邊值問題的通用形式在媒質交界面12 上,DN()() () xyxxyynuuufxxyyupuuuqxyeee1211221122()()xxyynxxyynuuuuuuxyxyeeeeee2021/2/1167單元矩陣元素計算公式:( )d dejjeiiijxyijNNNNKN Nx yxxyy單元右端項計算公式:()()xyuuufxxyy二維邊值問題一般形式( )d

34、deeiibf Nx y在媒質交界面12 上的條件自動滿足;第一類邊界條件需要強加;第三類邊界條件的計算參看金建銘電磁場有限元方法。2021/2/1168二維靜電場:0D22()() 0 ()lnrruuxxyyupuurxyrrr無限遠邊界條件()()xyuuufxxyy靜電場分析lnln lnlnuAAruuArrrrrrr電力線平行與垂直邊界條件2021/2/1169軸對稱靜電場:0D22()() 0 ()rruurrrrrzzupuuRrzRR無限遠邊界條件()()xyuuufxxyy靜電場分析211 AuAAuuRRRR RR 2021/2/1170二維靜磁場:0D2211()() 0 ()lnzzzrrzz