5有限元法.doc

第五章 有限元法有限元法是以變分原理為基礎(chǔ)，將要求解的微分方程型數(shù)學模型邊值問題，首先轉(zhuǎn)化為相應的變分問題，即泛函求極值問題；然后，利用剖分插值將變分問題離散化為普通多元函數(shù)的極值問題，最終歸結(jié)為一組多元的代數(shù)方程組，求解該方程組，從而獲得邊值問題的數(shù)值解。有限元法的核心是剖分插值，即將連續(xù)場分割為有限個單元，然后用較簡單的插值函數(shù)來表示每個單元的解。5.1 變分原理與尤拉方程 dxdsA(x1, y1)B(x2, y2)xy圖5.1 最速降線問題O在微積分學形成初期，以數(shù)學物理問題為背景，與多元函數(shù)的極值問題相對應，已在幾何、力學上提出了若干個求解泛函極值的問題。如圖5.1中的質(zhì)點最速降線問題所述，質(zhì)點A從定點自由下滑到定點B，試求使滑行時間最短的質(zhì)點下滑軌道。圖示滑行弧段所需時間為滑行總時間為（5-1）（5-1）式不僅取決于積分端點和，而且取決于的選取。取決于，所以是函數(shù)的函數(shù)，稱之為的泛函，記作。于是所述之最速降線問題，在數(shù)學上就歸結(jié)為研究泛函的極值問題，即 （5-2）泛函的極值（max或min）問題就稱為變分問題。對一般問題而言，可導出下列對應于一個自變量、單個函數(shù)及其導數(shù)的已知函數(shù) （5-3）式中為、和的已知函數(shù)。泛函的自變量不是一般的自變量，而是一個或幾個函數(shù)所屬的函數(shù)族。在端點和上分別等于給定值的無數(shù)個函數(shù)中，僅有一個能使定積分達到極小值，此函數(shù)稱為極值函數(shù)。因此，變分問題就在于尋求使泛函達到極值的該極值函數(shù)，即分析研究泛函的極值問題。泛函變分問題的經(jīng)典解法有兩種，一種稱之為直接解法，另一類是間接解法。直接解法是直接把泛函的極值問題近似地轉(zhuǎn)化為一般多元函數(shù)的極值問題，用有限維子空間中的函數(shù)去逼近無窮維空間中的極值函數(shù)，從而近似求得泛函的極值。間接解法是將變分問題轉(zhuǎn)化為尤拉方程（微分方程）的定解問題，即邊值問題來求解。我們以式（5-3）這種最簡形式來推導尤拉方程。設(shè)函數(shù)稍有變化，記作，稱之為的變分，它反映了整個函數(shù)的變化量。這樣泛函的值也應隨之變動，相應于變分的泛函增量為（5-4）將（5-4）式由多元函數(shù)的泰勒公式展開（5-5）式中作為泛函增量的線性主部為（5-6）稱為泛函的一次變分（簡稱變分）。而、分別是函數(shù)變分及其導數(shù)的二次、三次齊次式等的積分，依次稱為二次變分，三次變分令變分問題的解為，且設(shè)極值解稍有變動，且令（5-7）式中為任意給定的微量實參數(shù)，值就確定了函數(shù)族中的某一曲線，進而確定泛函之值；而是定義于區(qū)間且滿足齊次邊界條件的可微函數(shù)。于是泛函就成為變量的函數(shù)，且當時獲極值函數(shù)的解。在時取得極值的必要條件是（5-8）（5-9） 故 （5-10）簡寫為（5-11）（5-11）與（5-6）式比較，只差一個數(shù)值因子，故極值函數(shù)解必須滿足的必要條件（5-8）等同于（5-12）還可寫成（5-13）利用分部積分，并根據(jù)變分與微分順序可互換原理，（5-13）可寫為（5-14）在變分問題中，變分通常在端點保持為零，即于是（5-14）可寫為（5-15）由于（5-15）對任意均成立，故只有（5-16）方程（5-16）就稱為泛函（5-3）的極值問題的尤拉方程與上述過程類似，可繼續(xù)推導出各種復雜情況下的泛函極值存在的必要條件。例如在二維電磁場問題中，對應的泛函取決于一個二元函數(shù)，相應的泛函為（5-17）其極值存在的必要條件為偏微分方程（5-18）式中， （且，）。（5-19）在這里，就是對x的全偏導函數(shù)，此時應看作是固定的，而把、和看作是依賴于的，因此有 同理（5-20）5.2 有限元基本原理接5.1節(jié)有限元法的變分原理，通常有限元法的應用步驟為i 給出于待求邊值問題相應的泛函及其等效變分問題；ii 應用有限單元剖分場域，選取相應的插值函數(shù)iii 將變分問題離散化為一個多元函數(shù)的極值問題，導出一組聯(lián)立的代數(shù)方程；iv 選擇適當?shù)拇鷶?shù)解法，解有限元方程，即可得邊值問題的近似解（數(shù)值解）。1.電磁場邊值問題及其等價變分問題眾所周知，變分原理的應用實質(zhì)上是對物理學定律的一種重新描述，如電磁邊值問題中解電學的湯姆遜定律即是描述靜電現(xiàn)象的“最小作用原理”。湯姆遜定律指出：處于介質(zhì)中一個固定的帶電導體系統(tǒng)，其表面上電荷的分布，應使合成的靜電場具有最小的靜電能量。因此，任一由n個帶電導體構(gòu)成的二維靜電場問題的規(guī)律性可通過能量積分表示為（5-21）式中為媒質(zhì)介電系數(shù)，為靜電場能量，為靜電場場強，為靜電場電位，為研究域，且每一帶電導體的電位的邊值條件為（5-22）對照（5-21）和（5-17），靜電能量積分即是一類取決于二元位函數(shù)分布的泛函。因而，根據(jù)湯姆遜定理，二維靜電場的規(guī)律性就歸結(jié)為下述變分問題（5-23）（5-23）式的解答，即其極值函數(shù)的解答應滿足尤拉方程（5-18）。并在邊界上滿足相應的邊界條件（5-22）。將（5-23）式中函數(shù)以（5-18）尤拉方程中相關(guān)各項的運算，即得場與變分問題（5-23）對應的尤拉方程為（5-24）可見，由變分問題（5-23）給出的極值函數(shù)應滿足具有給定邊值（5-22）的拉普拉斯方程（5-24），顯然（5-24）和（5-22）一起構(gòu)成大家熟知的第一類邊值問題。與此相仿，通過尤拉方程，可知與下述變分問題（5-25）等價的邊值問題是（5-26）（5-26）式即泊松方程的第一類邊值問題。對第二、三類邊值問題(式中即為第二類邊值問題） （5-27）其等價變分問題可表示為（5-28）由（5-28）式可見，第二類或第三類邊界條件在變分問題中已被包含在泛函達到極值的要求之中，不必單獨寫出，是自動滿足的，不必另行處置，故稱此種邊界條件為自然邊界條件，其相應的變分問題稱為無條件變分問題。但對于第一類邊界條件，必須作為定解條件列出（式（5-25），故其變分問題求極值函數(shù)時必須在滿足這一類邊界條件的函數(shù)中去尋求。因此，稱這類邊界條件為強加邊界條件，其相應的變分問題（5-25）稱條件變分問題。在上述能量積分對應的泛函中，二次地依賴于函數(shù)及其偏導數(shù)，故又稱為函數(shù)的二次泛函。2. 波導場問題上述我們給出靜電場問題的變分問題。對時諧場，我們應從時變場著手去分析?，F(xiàn)以波導時諧場為例進行分析。分析波導中電磁波傳播問題時，為簡化分析，我們可設(shè)i 波導壁由純導體構(gòu)成；ii 波導中無自由電荷和傳導電流；iii 波導工作于匹配狀態(tài)，即只考慮向前傳播的入射波，無反射波。根據(jù)波導理論，由激勵源激發(fā)的形式不同，波導內(nèi)傳播的波可分為橫電波（TE波）和橫磁波（TM波）兩種，且一旦場量的向分量（縱向分量）和確定后，其相應的橫向分量和便可求出。因此，對波導場的分析，可歸結(jié)為定解縱向分量所對應的波動方程。對TE波應取為分析對象，由矢量波動方程在直角坐標系中的展開式（1-8）可知滿足如下波動方程（5-29）同理，對TM波則歸結(jié)為對相應的縱向分量的波動方程（5-30）由于、沿的行波特性，且令此向每單位長度中相位變化的相位系數(shù)為，則若以標記相應的和，則波導場分析可歸結(jié)為如下定義的波導橫截面平面內(nèi)的二位標量波動方程（亥姆霍茲方程）的解，即（5-31）式中（5-32）由此可知描述波導場的定解問題為（5-33）式中L為邊界波導壁。根據(jù)變分原理，與上述邊值問題對應的泛函為：（5-34）因此，與TE波波導場定解問題等價的無條件變分問題為：（5-35）而與TM波波導場定解問題等價的為條件變分。對泛函（5-35）取極值，即由泛函變分，再經(jīng)有限元離散化處理，便可求得相應的有限元方程。3.變分問題的離散化與有限元方程對平面域D進行離散化（剖分）處理時，可采用多種幾何剖分與相應的分片插值法。這里以常用的三角剖分及相應的三頂點線性插值為例，討論場域D的三角剖分問題，如圖5.2所示。將D域剖分為有限個互不重疊的三角形有限單元（簡稱三角元）。要求任意三角元的頂點必須也是其相鄰三角元的頂點，而不能是相鄰三角元邊上的內(nèi)點。當有不同媒質(zhì)的分界線時，不允許有跨越分界線的三角元。剖分一直延到邊界L，如邊界線為曲線，則應以三角元的一邊去逼近。三角元可大可小，應根據(jù)計算精度要求，確定剖分密度。對三角元頂點的編號，以壓縮存儲量，簡化計算程序及計算量為準。圖中給出一種三頂點編號。對任一三角元e（單元編號e=1,2,，e0）其三頂點的節(jié)點編號，以逆時針順序建立局部編碼序，標記為。i(xi, yi)LxyDL圖5.2 場域D的三角剖分示意圖m(xm, ym)j(xj, yj)基于上述剖分，在三角元e內(nèi)，分別給定對呈線性變化的插值函數(shù)（5-36）以此近似替代該三角元內(nèi)的待求變分問題的解，式（5-36）中待定系數(shù)、和，可由如下聯(lián)立方程求解：（5-37）式中，為三角元面積。而各系數(shù)可按下標順序置換而得。于是可得定義于三角元e上的線性插值函數(shù)為（5-38）式中稱三角元e上的線性插值基函數(shù)（或稱形狀函數(shù)），他取決于單元的形狀及其相應節(jié)點的配置，記作（5-39）由此，式（5-38）可簡潔的以矩陣形式寫為：（5-40）由于相關(guān)的三角元的公共邊及公共節(jié)點上的函數(shù)取值相同，故可以將每個三角元上構(gòu)造的線性插值函數(shù)進行拼合，使整個D域用拼合的分片線性插值函數(shù)描述。顯然，它取決于待求函數(shù)在各個節(jié)點上的值（為總結(jié)點數(shù)）。下面我們以二維拉普拉斯場的第一類邊值問題所對應的變分問題為例，說明變分問題的離散化過程。i 單元分析根據(jù)三角剖分，二次泛函可表示為遍及所有單元的能量積分的總和，即（5-41）式中表示三角元e所對應的能量積分，且由式（5-38）可得所以（5-42）式中 (5-43)同理可得（5-44）式中（5-45）因此（5-46）式中（5-47）此三階方陣式單元電場能的離散矩陣，稱為單元電場能系數(shù)矩陣，是對稱陣，其元素一般式為：（5-48）ii 總體合成由式（5-41）可知，為得到D域內(nèi)關(guān)于節(jié)點電位的離散表達時，首先應將式（5-46）中的擴充為，此矩陣式全部節(jié)點處電位值按節(jié)點編號順序排列成的一個階列陣；同時將擴充為，此系在式（5-47）所示的基礎(chǔ)上，按節(jié)點編號順序展成行與列，構(gòu)成階方陣，其中除行、列數(shù)分別為時有九個原的元素外，階方陣中，其余元素均為零。于是（5-46）式可改寫為：（5-49）總體能量積分二次泛函可離散為（5-50）式中稱為總電場能系數(shù)矩陣。因為，可見其元素應為（5-51）由以上可判定，是對稱陣（即）。由式（5-50）可將變分問題（5-23）離散化為：（5-52）根據(jù)函數(shù)極值條件，應有故由式（5-52）可得即或以矩陣表示為待解的多元線性代數(shù)方程組（5-53）此方程即為有限元方程。在獲得有限元方程后，可用關(guān)于此種代數(shù)方程組的各種計算機求解方法，如高斯消去法，列主元消去法，共軛梯度加速迭代法等進行求解。5.3 有限元應用實例作為典型示例，選取矩形波導BJ-100（）中TE波的截止波長的分布問題進行分析。式（5-33）中對TE波導場定解問題等價的無條件變分問題為（5-35）令而按5.2中和可得與式（5-46）相同的表達式其中與式（5-47）相同，其元素如（5-48）給出。經(jīng)過總體合成，及運用函數(shù)極值條件，可得如下有限元方程：（5-54）其中，對應于所選取的有限單元e，其各個單元矩陣的元素分別為（5-55）（5-56）和為形狀函數(shù)，由式（5-39）給定，式（5-54）是廣義代數(shù)特征值問題。式（5-54）中為對稱陣，為對稱正定矩陣，為求解式（5-54）應首先把廣義特征值問題變換為對稱陣的特征值問題，為此，將對稱正定陣利用平方根法（Cholesky）分解為下三角陣與其轉(zhuǎn)置的乘積，即（5-57）為便于書寫，將上式以黑體字寫為由式（5-54）可導得 (5-58)式中，令則（5-54）式的廣義特征值問題的求解可化為對稱陣的特征值問題，即的特征值問題，而原問題的特征向量現(xiàn)變換為，在求得特征向量后，再經(jīng)過的變換，方能求得原問題的特征向量。限于篇幅，這里不再贅述求解特征向量的全過程，只要將上述特征值問題利用豪斯豪爾法（HouseHolder）變換，把對稱陣化為對稱三角矩陣，再利用實對稱三角矩陣，經(jīng)過多次相似變換，再化做滿足指定精度的對角陣即可。因為對角陣的每個元素就是其特征值。從而求得特征向量，再經(jīng)過反變換求出原問題的特征向量。對應于由（5-54）求得的一系列特征值，其中非負的最小非零特征值就給出相應波導中最低型（主模）的截止波長的解答。數(shù)值解與理論值之間的對比結(jié)果如表（5-1）給處。表5.1 數(shù)值計算結(jié)果與理論值的比較波型截止波長(cm)理論值數(shù)值解相對誤差（）TE104.5724.5124.5694.5740.066TE202.2862.1932.2622.2751.0