有限元法

　　把物體或結構整體所具有的域V劃分為有限多個被稱為單元的子域V_n，以求得近似解的一種數值計算方法。劃分整體的離散化思想，早在20世紀40年代就已經有人提出。50年代M.J.特納等人首先提出對連續體進行離散化的分析法。近年來，在電腦的配合下，它已不同程度地推廣應用到許多工程技術領域，成為處理工程計算問題的有效方法之一。

　　若單元是一根直的細長桿，可以為隻沿桿的軸線方向有一個度量尺寸時，稱為一維單元。圖1所示的等厚度平面板被劃分為28個小三角形平面板，它們很薄，其“中面”是二個度量尺寸的平面，稱為二維單元。除三角形外，二維單元還可以是矩形、梯形、四邊形或其他形狀。若單元是長方形體、四面體、六面體、或其他形狀的立體，具有三個度量尺寸，則稱為三維單元。

　　離散後的物體或結構中相鄰單元之間，通常隻在若幹點上相互連接，這些點稱為節點；相鄰單元在圖形上相吻合的邊不相接觸，而隻通過邊上的節點傳力。此外，由於計算工作的需要，單元內部也可設置內節點（圖2）。所有荷載包括體力和面力都要移置到節點上；支座也要設置在節點處。節點可視為鉸結或固結或其他狀況的連接。如上所述，有限元法所分析的是與原物體或原結構的材料相同、形狀接近而由許多個單元以一定的方式連接起來的離散物體，計算所得應力、變形等結果是近似的，但當劃分的單元合理且數量很多時，即可得到非常逼近於真實的解。

　　在有限元法的應用中，最常用的計算法是以節點位移為基本未知量的位移法，另外還有力法和其他一些特殊的方法。采用位移法計算時，首要的問題是采用什麼樣的插值函數來建立單元的位移模式。

　　位移模式和插值函數　物體或結構離散化後，在單元中各點位移的變化可采用逼近於其原函數的近似函數來描述，這種近似函數稱為位移模式，通常用插值多項式表示。如圖1所示最簡單的三節點三角形單元的位移模式，可表示為如下的線性插值多項式：

　(1)

式中 α_m( m＝1，2，…，6)是待定的參數，稱為廣義坐標。

　　取1、2、3三個節點為插值點，使在每個插值點上有：

　　(2)

式中 x _i、 y _i為節點的坐標； u _i、 v _i為節點線位移。由(2)式6個公式可定出6個廣義坐標，於是可將(1)式所示的插值多項式改寫為：

　　(3)

而位移模式列陣為：

式中δ ^θ＝[ u ₁ v ₁ u ₂ v ₂ u ₃ v ₃] ^T，稱為單元的節點位移矢量矩陣；/是二階的單位矩陣；而　　　　

　　　　(4)

式中 A為三角形單元的面積；

(5)

在單元圖中，三節點1、2、3的次序必須是逆時針轉向。

　　由(4)式可見N_i(i＝1，2，3)是坐標(x，y)的函數，在數學上稱為插值的基函數，在力學上由於它反映單元的位移形態，所以稱為位移的形態函數，簡稱形函數。

　　在多節點三角形單元中，最常用的是在三角形123中三個邊的中點增加三個節點4、5、6（圖3）。在三角形內一點P的位置，可用面積坐標L₁＝A₁/A、L₂＝A₂/A、L₃＝A₃/A表示。取6個節點為插值點，位移模式可表示為完全的二次插值多項式，並類似於(3)式將插值多項式寫為：

　　　　(6)

式中的形函數可用面積坐標表示：

　　　(7)

(2)、(3)和(6)式中所示的節點線位移 u _i、 v _i稱為節點參數。各節點的線位移的導數也可取作節點參數。單元中所有節點的節點參數的數目總和稱為該單元的自由度。位移模式中廣義坐標 α_m( m＝1，2，…)的數目須等於或大於單元的自由度的數目。如圖3所示6個節點的三角形單元，若1、2、3三個主節點各取如下6個節點參數：兩個線位移 u _i、 v _i及其導數 ；而4、5、6三個副節點各取如下4個節點參數：兩個線位移 u_j、 v_j和沿法線方向的兩個導數 ，則該單元自由度的總數為3×6+3×4=30。位移模式可取完全四次插值多項式，此時廣義坐標 α_m( m＝1，2，…，30)的數目恰等於30。

　　更如圖1所示三節點三角形單元，設若每個節點的節點參數都是u_i、v_i、

、 、 、 　( i=1，2，3)，則單元的自由度數等於3×6＝18，位移模式可取完全三次插值多項式：

(8)

式中廣義坐標 α_m( m＝1，2，…，20)的數目為20，大於單元自由度數18。為瞭得解，必須增加兩個包括廣義坐標 α_m的方程式稱為約束方程式。

　　設有一多節點的矩形單元(圖4)，稱為一般拉族元，四邊節點2(m+n)個，內節點(m-1)(n-1)個。取各節點為插值點，表示位移模式的插值多項式為：

　　　　　　(9)

(i=0，1，2，…，m；j=0，1，2，…，n)

式中插值的基函數（即形函數），一般用拉格朗日多項式：

N_ij＝L怢(ξ)L怹(η)　　　　　　　(10)

式中 L怢( ξ)為拉格朗日乘子函數：　

　　　　(11)

將上式中的 _i換成 η， i和 m分別換成 j和 n，即得(10)式中的 L怹( η)。

　　對於圖5所示的正方形單元，是最簡單的4個節點的方形單元，當用線性拉族元的公式，得與四個節點相應的形函數為：

　(12)

式中 ξ _i、 η _i是節點 i的局部坐標，而該單元的位移模式為：

　　　(13)

　　四節點正方形或矩形單元不能適應曲線邊界和非正交的直線邊界，也不便隨意改變單元的大小，為此將正方形或矩形單元改變為圖6a所示的任意四邊形單元則可避免上述的缺點。

　　通過把整體坐標(x，y)變換為局部坐標（_ξ，η），圖6b所示正方形單元的位移模式和形函數式分別如(13)和(12)式所示，同樣適用於圖6a所示的四邊形單元，前者是為瞭後者而引出的單元稱為基本單元或母單元，後者是計算時所實際采用的單元。

　　圖6中所示的局部坐標（ξ，η）可依下列變換式轉換為整體坐標：

　　　　(14)

式中 N _i( i=1，2，3，4)仍為如(12)式所列包含局部坐標 ξ、 η的形函數； x _i、 y _i( i＝1，2，3，4)為四個節點的整體坐標值。

　　由於圖6a、b所示的兩種單元的位移模式(13)和坐標變換式(14)中所用的形函數是等同的，所以實際采用的四邊形單元稱為等參數單元或同參數單元。

　　空間問題所采用的最簡單的單元為圖7所示的四面體單元，其位移模式類似於平面問題中三節點三角形單元的(3)至(5)式。

　　對於空間單元，可類似於圖6b，取圖8所示的八節點立方體作為基本單元，仿照式(13)和式(14)得位移模式和坐標變換式。

　　單元分析　如圖7所示的空間四面體單元中，各節點上所承受的力U_i、V_i、W_i(i＝1，2，3，4)稱為節點力，其矢量矩陣為

　　　　　　　　　　

^θ＝[ U ₁ V ₁ W ₁… U ₄ V ₄ W ₄] ^T

單元內的應力矢量矩陣為

單元內的應力和節點力與節點位移的關系式分別為

┿＝

δ ^θ　　　　　 　　　　　(15)

^θ＝ δ ^θ　　　　　 　　　　(16)

以上二式中的 和 分別稱為單元的應力矩陣和剛度矩陣。所謂單元分析就是計算這兩個矩陣以及將單元的體力和面力等效地移置到節點上。

　　根據彈性力學適用於小位移的幾何方程和物理方程得如下的兩個關系式：

ε＝

δ ^θ　　 　　　　　 　　　(17)

┿＝

ε 　 　　　　　　(18)

式中ε為單元內的應變矢量矩陣

； 和 分別為應變矩陣和彈性矩陣。單元的應力矩陣為

＝

　　　　　　　　　(19)

　　剛度矩陣和體力、面力移置到節點上的計量公式，可用如下最小勢能原理導出。

　　在滿足小位移幾何方程和位移邊界條件的所有允許的應變和位移中，實際的應變和位移必使彈性體的總勢能

為最小：

(20)

式中第一項為單元的應變能； =[ u v w] ^T， =[ X Y Z] ^T分別為單元內各點的位移和體力的矢量矩陣； _b和 分別為單元的某一邊界面上的位移和面力的矢量矩陣。運用變分原理，得計算單元的剛度矩陣 和體力 、面力 移置到節點上的節點荷載公式為

　　　　　　　　　(21)

　　　　　

　　　　　　　　　(22)

　　　　　

　　　　　　　　　(23)

　　對於平面問題，隻須將以上三式中體積分轉化為面積分，面積分轉化為線積分。

　　對於等參單元，剛度矩陣公式和體力的移置公式隻須將式(21)和式(22)中的dV變換為｜J｜dξdηdζ，式中

　　　　　　　　　(24)

此時積分符號應改為 。對邊界面上承受面力（如在 ζ＝±1的 ηζ面上承受勻佈壓力 q）的荷載移置公式為 　　

　　總體合成和邊界條件的處理　把結構劃分為許多單元，求出各個單元的剛度矩陣（簡稱單剛），再利用它們求出整個結構的剛度矩陣，稱為總體剛度矩陣（簡稱總剛），才能達到求出結構中各節點的位移的目的。這種由部分求總合的工作（包括對邊界條件的處理），稱為總體合成，其方法有自由度編碼法和變換矩陣法，前者的優點是在計算機上容易實現，缺點是當節點參數包括有線位移的導數時，則難適用。此時，當采用後者。

　　自由度編碼法　設有包括1、2、…、5等5個節點的平面結構被劃分為Ⅰ、Ⅱ、Ⅲ三個單元，根據結構的左豎邊和下橫邊分別沒有橫向和豎向位移的邊界條件，佈置橫向和豎向支承鏈桿(圖9a)。

每個單元都有 i、 j、 m3個頂點，6個節點位移 u _i、 v _i( i， j， m)，因而有6個自由度，編碼為1、2、…、6，寫作節點力 F和節點位移 δ的下角標（圖9b），於是各單元的節點位移矢量矩陣和節點力矢量矩陣分別為

δ^θ=[u_iv_iu_jv_ju_mv_m]^T=[δ₁δ₂δ₃δ₄δ₅δ₆]^T

^θ=[ U _i V _i U_j V_j U_m V_m] ^T=[ F ₁ F ₂ F ₃ F ₄ F ₅ F ₆] ^T

式中u、U為橫向的節點位移和節點力；v、V為豎向的節點位移和節點力；δ和F為不分橫、豎而是按自由度編碼的順序來排列的節點位移和節點力。

　　對於整個結構，除去沿支承鏈桿方向的線位移為零不計外，共有5個節點位移，因而有5個自由度，編碼為1、2、…、5，寫作節點位移Δ和節點荷載P的下角碼(圖9c)。

　　　　　　　 Δ＝[Δ₁Δ₂Δ₃Δ₄Δ₅]^T

＝[ P ₁ P ₂ P ₃ P ₄ P ₅] ^T

　　單元節點的節點力包括加在該節點上的荷載(即集中加在節點上的荷載、體力和面力移置到節點上的荷載)和其他匯交於該節點的單元加在該節點上的力，但對於整個結構來說，後者互相抵消，因而節點上隻承受各類節點荷載的總和。

式中

和 分別為單剛和總剛。單剛中第 r行第 s列，總剛中第 R行第 S列的元素分別寫作 k_rs和 K_RS。

　　總剛中的各元素都是由一個或若幹個單元的單剛

中的元素來貢獻。單剛元素的下角標 rs和總剛元素的下角標 RS所代表的自由度兩兩相當時，該單剛元素 k_rs便將貢獻給 K_RS。對照圖9b和c可得單剛和總剛中兩兩相當的自由度的數碼，如第Ⅰ、Ⅱ、Ⅲ單元中自由度的數碼6、4、2與總剛中的數碼2相當。總剛中的元素通過如上找出相當的自由度數碼的方法，從若幹個單元中相當的元素求和來一一算得，從而確定總剛。

　　變換矩陣法　設有一平面結構共有n個單元，把各單元的節點位移矢量矩陣δ

、δ 、…、δ 按順序排列(這些單元中相同的元素並未歸並，都分別依序列入)，得：

　　　　　　　

＝[δ ₁ ^eδ ₂ ^e…δ ₃ ^e] ^T　　　　　　　(27)

　　把各單元的單剛

₁、 ₂、…、 _n作為主對角線的子矩陣，其餘為零，得如下的矩陣：

　　　　　　　　(28)

　　若整個結構共有m個節點，δ₁、δ₂、…、δ_m為各節點的節點位移矢量矩陣，按順序排列起來得

δ＝[δ₁δ₂…δ_m]^T　　　　　　　　(29)

　　

和δ兩列陣之間的關系式為

＝ δ　　　　　 　　　　 　(30)

式中 稱為變換矩陣。

　　求出

_s和 後便可依下式列出總剛

　　　 　　

＝ ^T _s (31)

　　列出瞭總剛和節點荷載矢量矩陣

以後，並處理瞭邊界條件，就可通過解式(26)那樣的線性代數方程組把基本未知量──各節點的位移求出。於是各單元內的應變、位移和應力就可分別依前邊列出過的公式算出。

　　收斂準則　有限元法的解答要求在單元的尺寸逐步取小時能夠收斂於正確解答。

　　選取的位移模式必須滿足如下的收斂準則：①每個單元的位移一般包括本單元的形變引起的位移和本單元以外其他單元發生的形變而連帶引起的位移兩個部分。選用的位移模式必須能反映後一部分的位移，即位移與本單元的應變無關時，將使本單元產生剛體位移；②每個單元的應變一般包括與各點的位置坐標有關的所謂變量應變和與該坐標無關的所謂常量應變兩個部分。選用的位移模式必須能反映常量應變。

　　以上兩條是能夠收斂於正確解答的必要條件。

　　另一個條件是相鄰單元之間的位移是協調的，即不出現圖10所示的縫隙或重合。相鄰單元位移的協調性是能夠收斂於正確解答的充分條件。

　　特種方法　從有限元的基本方法派生出來的方法很多，如有限條法、邊界元法、雜交元法、非協調元法和擬協調元法等，用以解決特殊的問題。

　　

參考書目

　馮康、石鐘慈：《彈性結構的數學理論》，科學出版社，北京，1981。

　徐次達、華伯浩：《固體力學有限元理論、方法及程序》，水利電力出版社，北京，1983。

　徐芝倫：《彈性力學問題的有限單元法》，水利電力出版社，北京，1978。

　O.C.Zienkiewicz，The Finite Element Method，McGraw-Hill，London，1977.