5有限元法.doc

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