有限元數值模擬:兩種四邊形元的計算比較

摘要文中綜述了單元剛度矩陣和載荷矩陣的形成，總體剛度矩陣和載荷矩陣的形成，以及對邊界條件的處理。比較了雙線性元在兩種不同坐標下的計算精度。

中圖分類號：O17 文獻標識碼：A

Finite Element Numerical Imitation：Tow Types of

Tetragonum Element Calculation Comparison

WU Tianzhi

(Pass College of Chingqing Technology and Business University， Chongqing 401520)

AbstractThis paper reviews the formation of element stiffness matrix and loading matrix， formation of total stiffness matrix and loading matrix， and the process of boundary conditions. Meanwhile， it compared the computational accuracy of bilinear elements under condition of different coordinate.

Key wordsfinite element; tetragonum; isoparmetric element; area coordinate; bilinear form; isoparametric coordinates

0 引言

有限元法，實質上就是Ritz-Galerkin法。它和傳統的Ritz-Galerkin法的主要的區別在于，它應用樣條函數方法提供了一種選取“局部基函數”或“分片多項式空間”的新技巧，從而在很大的程度上克服了Ritz-Galerkin法選取基函數的固有難度。有限元法成功地應用于結構力學、固體力學、流體力學、物理學和其它的工程學科。本文主要討論了雙線性形式下的有限元方法，介紹它在不同的兩種坐標下的計算過程。

1 單元剛度矩陣和單元載荷矩陣的形成

用有限元方法解二維問題，首先要對二維區域進行剖分。首先，通過四節點的等參元變換將凸四邊形區域變換到一個標準的正方形區域，然后對該正方形區域在橫向（分nelemi個單元）和縱向（分nelemj個單元）分別作剖分，對劃分出的每個小矩形，節點，小矩形的每個節點進行編號，最后，我們通過變換式將區域變換到原來的區域上。這樣就完成了對計算區域的剖分。

（1）等參元坐標下的單元形狀函數的構造、單元剛度矩陣和載荷矩陣的形成。

首先，在單元ei上：其次，在單元ei上來構造單元形狀函數：

在單元上ei，計算出v（看成u的形式）在單元ei上關于x，y的偏導數代入a (u，v)，這樣就得到雙線性形式下的單元剛度矩陣：

根據，在ei上進行積分得到：由此，得到單元載荷矩陣（不含的單元上）：在含的單元上只須再作一個線積分即可。在含的單元上的單元載荷矩陣為：到此，已經構造出雙線性形式下的單元剛度矩陣和單元載荷矩陣（參看[1]）。

（2）面積元坐標下的單元剛度矩陣和載荷矩陣的形成。

首先，定義g1，g2（如圖1），表示任意的四邊形的面積，用A1表示三角形123的面積，用A2表示三角形124的面積，g1 =A1/A，g2 =A2/A。（參見[2]、[3]）

定義p的面積元坐標為（L1，L2，L3，L4）=(a1/A，a2/A，a3/A，a4/A)（參見[2]），由此得到四邊形的每個頂點的坐標。已知在一個單元上的四個插值節點的坐標，類似地構造在該單元上的插值基函數，得到在面積元坐標下的單元形狀函數為： ，關于x，y求導數得到：在單元上進行面積元坐標下的積分（參看[3]）得到單元剛度矩陣和載荷矩陣：

（不含的單元）

（含的單元，如等參元所得）

2 由單元剛度矩陣和單元載荷矩陣組裝成總體剛度矩陣和總體載荷矩陣

在這一過程中，我們采用分單元計算的方法，從一號單元開始計算，按前面介紹的單元剛度矩陣和單元載荷矩陣的計算方法計算出該單元的剛度矩陣和載荷矩陣，然后就該單元上的剛度矩陣和載荷矩陣疊加到總體剛度矩陣和總體載荷矩陣中去，一直到將所有的單元循環完畢。就得到了我們所需要的總體剛度矩陣和總體載荷矩陣。

3 邊界條件的處理與精度分析

對第一類的邊值條件，為了保證總體剛度矩陣的對稱性，我們采用的處理方式是將有約束的節點在總體剛度矩陣中的對應的行和列的元素全置為零，將行和列的交叉處（即對角元）對應的元素置為1，將總體載荷矩陣中的對應行的元素置為約束的值，其余的元素減去約束值乘以該節點對應于剛度矩陣中列的對應元素，這樣將所有的約束節點處理完畢后就得到了有限元方程，這一方程的系數矩陣為對稱正定，保證了有限元的解的存在唯一性。由于對第二和第三邊值條件，我們可以在單元剛度矩陣和載荷矩陣的形成過程中來進行處理。在對有限元方程的求解過程中，我們將根據問題的規模和計算機的容量、速度，選取適當的求解方法。

在用等參元坐標進行有限元的計算過程中要涉及到復雜的坐標變換，涉及到雅可比行列式的計算，其計算積分采用的是數值積分方法，在這一系列的計算過程中，會因為算法的原因而造成一定的誤差；而在面積元坐標形式下的有限元的計算不需要進行坐標的變換，在計算單元的載荷矩陣和剛度矩陣時是采用的精確積分，這在算法上來看是不會有誤差產生。然而，因為在使用計算機進行計算的過程中，因計算機字長的原因，在我們所取的算例中，等參元坐標下的有限元計算的結果比在面積元坐標下的計算結果要好一些，這也就說明，算法的好壞并不絕對的依賴于理論的分析，還要綜合考慮各種因素，如計算機的字長等因素。通過對計算實例的計算，我們可以發現，每加密一次后的范數意義下的誤差都是前一次的四分之一，這說明在范數意義下，兩種方法皆為平方收斂。

參考文獻

[1]李榮華，馮果忱編.微分方程數值解（第三版）.北京:高等教育出版社，1996.

[2]Yuqiu Long，Yuxuan Li，Zhifei Long and Song Cen，‘Area coordinates used in quadrilateral elements’，Commun.Numer.Methods Eng.，15，533-545(1999).

[3]Zhifei Long，Juxuan Li，Song Cen and Yuqiu Long，‘Some basic formulae for area co-ordinatesin in quadrilateral elements’， Commun．Numer．Methods Eng.，15， 841-845(1999).

[4][美] J.E.艾金著.有限元法的應用與實現.張紀剛，郁衛中，林翠虹，譯.北京:科學出版社，1992.8.

[5]李人憲.有限元法基礎（第二版）.北京.國防工業出版社，2004.