有限元數(shù)值模擬:兩種四邊形元的計(jì)算比較

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

中圖分類號：O17 文獻(xiàn)標(biāo)識碼：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 引言

有限元法，實(shí)質(zhì)上就是Ritz-Galerkin法。它和傳統(tǒng)的Ritz-Galerkin法的主要的區(qū)別在于，它應(yīng)用樣條函數(shù)方法提供了一種選取“局部基函數(shù)”或“分片多項(xiàng)式空間”的新技巧，從而在很大的程度上克服了Ritz-Galerkin法選取基函數(shù)的固有難度。有限元法成功地應(yīng)用于結(jié)構(gòu)力學(xué)、固體力學(xué)、流體力學(xué)、物理學(xué)和其它的工程學(xué)科。本文主要討論了雙線性形式下的有限元方法，介紹它在不同的兩種坐標(biāo)下的計(jì)算過程。

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

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

（1）等參元坐標(biāo)下的單元形狀函數(shù)的構(gòu)造、單元剛度矩陣和載荷矩陣的形成。

首先，在單元ei上：其次，在單元ei上來構(gòu)造單元形狀函數(shù)：

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

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

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

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

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

（不含的單元）

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

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

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

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

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

在用等參元坐標(biāo)進(jìn)行有限元的計(jì)算過程中要涉及到復(fù)雜的坐標(biāo)變換，涉及到雅可比行列式的計(jì)算，其計(jì)算積分采用的是數(shù)值積分方法，在這一系列的計(jì)算過程中，會因?yàn)樗惴ǖ脑蚨斐梢欢ǖ恼`差；而在面積元坐標(biāo)形式下的有限元的計(jì)算不需要進(jìn)行坐標(biāo)的變換，在計(jì)算單元的載荷矩陣和剛度矩陣時(shí)是采用的精確積分，這在算法上來看是不會有誤差產(chǎn)生。然而，因?yàn)樵谑褂糜?jì)算機(jī)進(jìn)行計(jì)算的過程中，因計(jì)算機(jī)字長的原因，在我們所取的算例中，等參元坐標(biāo)下的有限元計(jì)算的結(jié)果比在面積元坐標(biāo)下的計(jì)算結(jié)果要好一些，這也就說明，算法的好壞并不絕對的依賴于理論的分析，還要綜合考慮各種因素，如計(jì)算機(jī)的字長等因素。通過對計(jì)算實(shí)例的計(jì)算，我們可以發(fā)現(xiàn)，每加密一次后的范數(shù)意義下的誤差都是前一次的四分之一，這說明在范數(shù)意義下，兩種方法皆為平方收斂。

參考文獻(xiàn)

[1]李榮華，馮果忱編.微分方程數(shù)值解（第三版）.北京:高等教育出版社，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.艾金著.有限元法的應(yīng)用與實(shí)現(xiàn).張紀(jì)剛，郁衛(wèi)中，林翠虹，譯.北京:科學(xué)出版社，1992.8.

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