1、第六章 單元形函數的構造16.1形函數構造的一般原理16.2形函數的性質76.3用面積坐標表達的形函數86.4有限元的收斂準則106.5 等效結點載荷列陣116.5.1 單元載荷的移置116.5.2 結構整體載荷列陣的形成116.5.3載荷移置與靜力等效關系12習題14第六章 單元形函數的構造在有限單元法的基本理論中，形函數是一個十分重要的概念，它不僅可以用作單元的內插函數，把單元內任一點的位移用結點位移表示，而且可作為加權余量法中的加權函數，可以處理外載荷，將分布力等效為結點上的集中力和力矩，此外，它可用于后續的等參數單元的坐標變換等。根據形函數的思想，首先將單元的位移場函數表示為多項式的形

2、式，然后利用結點條件將多項式中的待定參數表示成場函數的結點值和單元幾何參數的函數，從而將場函數表示成結點值插值形式的表達式。在本節中，重點討論幾種典型單元的形函數插值函數的構造方式，它們具有一定的規律。然后以平面三角形單元為例，討論了形函數的性質，在此基礎上分析了有限元的收斂準則。6.1形函數構造的一般原理單元的類型和形狀決定于結構總體求解域的幾何特點、問題類型和求解精度。根據單元形狀,可分為一維、二維、三維單元。單元插值形函數主要取決于單元的形狀、結點類型和單元的結點數目。結點的類型可以是只包含場函數的結點值，也可能還包含場函數導數的結點值。是否需要場函數導數的結點值作為結點變量一般取決于單

3、元邊界上的連續性要求，如果邊界上只要求函數值保持連續，稱為C0型單元，若要求函數值及其一階導數值都保持連續，則是C1型單元。在有限元中，單元插值形函數均采用不同階次的冪函數多項式形式。對于C0型單元，單元內的未知場函數的線性變化僅用角（端）結點的參數來表示。結點參數只包含場函數的結點值。而對于C1型單元，結點參數中包含場函數及其一階導數的結點值。與此相對應，形函數可分為Lagrange型（不需要函數在結點上的斜率或曲率）和Hermite型（需要形函數在結點上的斜率或曲率）兩大類，而形函數的冪次則是指所采用的多項式的冪次，可能具有一次、二次、三次、或更高次等。另外，有限元形函數N是坐標x、y、z

4、的函數，而結點位移不是x、y、z的函數，因此靜力學中的位移對坐標微分時，只對形函數N作用，而在動力學中位移對時間t微分時，只對結點位移向量作用。 （1）一維一次兩結點單元 圖3-8 一維一次兩結點單元模型設位移函數u(x)沿x軸呈線性變化,即 (5.90)寫成向量形式為 (5.91)設兩個結點的坐標為；兩結點的位移分別為，可以代入上式并解出，得 (5.92)位移函數u(x)記作形函數與結點參數乘積的形式 (5.93)得到形函數為 (5.94)在自然坐標系內進行定義，則可得到形函數的標準化形式 (5.95)其中，自然坐標的變換公式為。圖3-9一維一次兩結點單元的局部坐標表達（2）二維一次三結點單

5、元（平面三角形單元）在總體坐標系統下，任一點的某一方向的位移是 (5.96)設三個結點的坐標是，為三個結點在某方向上的位移，具有如下關系 (5.97)得到形函數矩陣如下式 (5.98)上述推導可用如下MATLAB程序實現：clearv=sym('1, x,y')m=sym('1,x1,y1;1,x2,y2;1,x3,y3')mm=inv(m)N=v*mmsimplify(factor(N)（3）三維一次四結點單元（三維四面體單元）在總體坐標系統下，任一點的某一方向的位移是 (5.99)按相似的方法可以得到 (5.100)形函數矩陣如下式 (5.101)（4）一維

6、二次三結點單元（高次單元）圖3-10一維二次三結點單元模型設位移函數為 (5.102)用結點位移代入并求解, (5.103)得到 (5.104)上式等號右端第一項矩陣即為形函數。（5）一維三次四結點單元（Lagrange型） 圖3-11 一維三次四結點單元模型位移函數為三次方程 (5.105)需要四個結點參數才能唯一地確定其中的常系數。這四個結點可以分別取兩個端點和兩個三分點。類似地，可以得到如下形函數方程 (5.106)其中形函數中的各元素為，. (5.107)（6）一維三次二結點單元（Hermite型）（平面梁單元）圖3-12 一維三次二結點單元這類單元的位移函數為 (5.108)對應的轉

7、角方程為 (5.109)用結點參數代入求解，即 (5.110)得到 (5.111)其中形函數矩陣中各元素為， (5.112)上述結果可用MATLAB程序進行驗證：clearx=sym('x');j=0:3;v=x.j % v=1 x x2 x3;m=sym('1,x1,x12,x13;1,x2,x22,x23;0,1,2*x1,3*x12;0,1,2*x2,3*x22')mm=inv(m)N=v*mm;simplify(factor(N)（7）二維一次四結點單元（平面四邊形單元或矩形單元）用形函數表達的位移方程如下 (5.113)其中形函數矩陣的元素為，i=1,

8、2,3,4 (5.114)對于平面四邊形單元和矩形單元,可用局部坐標系統很好地加以解釋。局部坐標的范圍定義為-1+1,四個結點的值固定。局部坐標系下的形函數為 (5.115) 圖3-13 二維一次四結點單元（8）三維一次八結點單元在三維一次單元形函數中，函數值沿三坐標軸（x、y、z軸）呈線性變化。假設位移函數沿各坐標軸的線性變化可寫成 (5.116)假設在i結點的位移值為ui，并將數值代入上式，其他各結點(j,k,l,m,n,p,q)亦類推，共有8個式子，其中第1式如下 (5.117)可是以求得系數解 (5.118)則有 (5.119)最后得到形函數的表達式為 (5.120)（9）帕斯卡三角形

9、上述各種位移函數的構造有一定的規律，可以根據所謂的帕斯卡三角形加以確定，同時，這樣制定的位移模式，還能夠滿足有限元的收斂性要求。以下是幾種典型情況。一維兩結點單元的情況：圖3-14 一維兩結點單元的變量組成一維三結點單元的情況：圖3-15 一維三結點單元的變量組成二維高階單元的情況：常數項線性項二次項三次項四次項五次項圖3-16 二維高階單元的變量組成三維四結點單元的情況：圖3-17 三維四結點單元的變量組成6.2形函數的性質下面以平面三角形單元為例討論形函數的一些性質。平面三角形單元的形函數為， （i =1, 2 , 3） (a)其中，為三角形單元的面積，為與結點坐標有關的系數，它們分別等于

10、公式中的行列式的有關代數余子式，即a1 、b1 、c1 ，a2 、b2 、c2 和a3 、b3 、c3 分別是行列式中的第一行、第二行和第三行各元素的代數余子式。對于任意一個行列式, 其任一行（或列）的元素與其相應的代數余子式的乘積之和等于行列式的值，而任一行（或列）的元素與其他行（或列）對應元素的代數余子式乘積之和為零。因此有：第一，形函數在各單元結點上的值，具有“本點是1、它點為零”的性質，即在單元結點1上，滿足 (b)在結點2、3上，有 (c) (d)類似地有 (e)第二，在單元的任一結點上，三個形函數之和等于1，即 (f)簡記為 (g)這說明，三個形函數中只有二個是獨立的。第三，三角形

11、單元任意一條邊上的形函數，僅與該邊的兩端結點坐標有關、而與其它結點坐標無關。例如，在23 邊上有 (h)這一點利用單元坐標幾何關系很容易證明。根據形函數的這一性質可以證明，相鄰單元的位移分別進行線性插值之后，在其公共邊上將是連續的。例如，單元123和124具有公共邊12。由上式可知，在12邊上兩個單元的第三個形函數都等于0，即 (i)不論按哪個單元來計算，公共邊12上的位移均由下式表示 (j)可見，在公共邊上的位移u、v 將完全由公共邊上的兩個結點1、2的位移所確定，因而相鄰單元的位移是保持連續的。6.3用面積坐標表達的形函數為了能夠更好地理解形函數的概念，這里引入面積坐標。在如圖5.18所示

12、的三角形單元ijm中，任意一點P(x , y)的位置可以用以下三個比值來確定圖3-18 平面三角形單元的面積坐標 (5.121)式中，D三角形單元ijm的面積，Di 、Dj 、Dm 三角形Pjm、Pmi、Pij的面積。Li ，Lj ，Lm叫做P點的面積坐標。顯然，這三個面積坐標不是完全獨立的，這是由于Di +Dj +Dm =D (5.122)所以有Li +Lj +Lm =1 (3-123)對于三角形Pjm，其面積為 (5.124)故有 (5.125)類似地有 (5.126) (5.127)可見，前面講述的平面三角形單元的形函數Ni 、Nj 、Nm 等于面積坐標Li 、Lj 、Lm 。容易看出，

13、單元三個結點的面積坐標分別為結點 i： Li =1 Lj =0 Lm =0結點 j： Li =0 Lj =1 Lm =0 結點m： Li =0 Lj =0 Lm =1根據面積坐標的定義，平行于jm邊的某一直線上的所有各點都有相同的坐標Li，并且等于該直線至jm邊的距離與結點i至jm邊的距離之比，圖3-18中給出了Li的一些等值線。平行于其它邊的直線也有類似的情況。不難驗證，面積坐標與直角坐標之間還存在以下變換關系： (5.128)當面積坐標的函數對直角坐標求導時，有下列公式： (5.129)求面積坐標的冪函數在三角形單元上的積分時，有 (5.130)式中, a、b、g 為整常數。求面積坐標的冪

14、函數在三角形某一邊上的積分值時，有 (5.131)式中, l為該邊的長度。6.4有限元的收斂準則對于一個數值計算方法，一般總希望隨著網格的逐步細分所得到的解答能夠收斂于問題的精確解。根據前面的分析，在有限元中，一旦確定了單元的形狀，位移模式的選擇將是非常關鍵的。由于載荷的移置、應力矩陣和剛度矩陣的建立都依賴于單元的位移模式，所以，如果所選擇的位移模式與真實的位移分布有很大的差別，會將很難獲得良好的數值解。可以證明，對于一個給定的位移模式，其剛度系數的數值比精確值要大。所以，在給定的載荷之下，有限元計算模型的變形將比實際結構的變形小。因此細分單元網格，位移近似解將由下方收斂于精確解，即得到真實解

15、的下界。為了保證解答的收斂性，位移模式要滿足以下三個條件，即 位移模式必須包含單元的剛體位移。也就是，當結點位移由某個剛體位移引起時，彈性體內將不會產生應變。所以，位移模式不但要具有描述單元本身形變的能力，而且還要具有描述由于其它單元形變而通過結點位移引起單元剛體位移的能力。例如，平面三角形單元位移模式的常數項a1、a4 就是用于提供剛體位移的。 位移模式必須能包含單元的常應變。每個單元的應變一般包含兩個部分：一部分是與該單元中各點的坐標位置有關的應變，另一部分是與位置坐標無關的應變（即所謂的常應變）。從物理意義上看，當單元尺寸無限縮小時，每個單元中的應變應趨于常量。因此，在位移模式中必須包含

16、有這些常應變，否則就不可能使數值解收斂于正確解。很顯然，在平面三角形單元的位移模式中，與a2、a3、a5、a6 有關的線性項就是提供單元中的常應變的。 位移模式在單元內要連續、且在相鄰單元之間的位移必須協調。當選擇多項式來構成位移模式時，單元內的連續性要求總是得到滿足的，單元間的位移協調性，就是要求單元之間既不會出現開裂也不會出現重疊的現象。通常，當單元交界面上的位移取決于該交界面上結點的位移時，就可以保證位移的協調性。在有限單元法中，把能夠滿足條件1和2的單元，稱為完備單元；滿足條件3的單元，叫做協調單元或保續單元。前面討論過的三角形單元和矩形單元，均能同時滿足上述三個條件，因此都屬于完備的

17、協調單元。在某些梁、板及殼體分析中，要使單元滿足條件3會比較困難，實踐中有時也出現一些只滿足條件1和2的單元，其收斂性往往也能夠令人滿意。放松條件3的單元，即完備而不協調的單元，已獲得了很多成功的應用。不協調單元的缺點主要是不能事先確定其剛度與真實剛度之間的大小關系。但不協調單元一般不象協調單元那樣剛硬（即比較柔軟），因此有可能會比協調單元收斂得快。在選擇多項式作為單元的位移模式時，其階次的確定要考慮解答的收斂性，即單元的完備性和協調性要求。實踐證明，雖然這兩項確實是所要考慮的重要因素，但并不是唯一的因素。選擇多項式位移模式階次時，需要考慮的另一個因素是，所選的模式應該與局部坐標系的方位無關，

18、這一性質稱為幾何各向同性。對于線性多項式，各向同性的要求通常就等價于位移模式必須包含常應變狀態。對于高次位移模式，就是不應該有一個偏移的坐標方向，也就是位移形式不應該隨局部坐標的更換而改變。經驗證明，實現幾何各向同性的一種有效方法是，可以根據巴斯卡三角形來選擇二維多項式的各項。在二維多項式中，如果包含有對稱軸一邊的某一項，就必須同時包含有另一邊的對稱項。選擇多項式位移模式時，還應考慮多項式中的項數必須等于或稍大于單元邊界上的外結點的自由度數。通常取項數與單元的外結點的自由度數相等，取過多的項數是不恰當的。6.5 等效結點載荷列陣在結構有限元整體分析時，結構的載荷列陣R是由結構的全部單元的等效結

19、點力集合而成，而其中單元的等效結點力Re 則是由作用在單元上的集中力、表面力和體積力分別移置到結點上，再逐點加以合成求得。本節以平面三角形單元為例，討論集中力、表面力和體積力的等效移置方法以及如何形成結構等效載荷列陣，并與靜力等效進行了對比。6.5.1 單元載荷的移置根據虛位移原理，等效結點力所做的功與作用在單元上的集中力、表面力和體積力在任何虛位移上所做的功相等，由此確定等效結點力的大小。對于平面三角形單元，有 (5.132)式中，單元結點虛位移列陣，單元內任一點的虛位移列陣；等號左邊表示單元的等效結點力Re 所做的虛功；等號右邊第一項是集中力G所做的虛功，等號右邊第二項是面力q所做的虛功，

20、積分沿著單元的邊界進行；等號右邊第三項表示體積力p所做的虛功，積分遍及整個單元；t為單元的厚度，假定為常量。用形函數矩陣表示的單元位移模式方程為 (5.133)代入式(3-132)，注意到結點虛位移列陣d*e可以提到積分號的外面，于是有 (5.134)注意到(d * e ) T 的任意性，上式化簡為R e = F e +Q e +P e (5.135)其中 F e = N T G (5.136) (5.137) (5.138)式(5.134)右端括號中的第一項與結點虛位移相乘等于集中力所做的虛功，它是單元上的集中力移置到結點上所得到的等效結點力，它是一個6×1階的列陣，記為F e。同

21、理，式(5.134)右端括號中的第二項是單元上的表面力移置到結點上所得到的等效結點力，記為Qe；第三項是單元上的體積力移置到結點上所得到的等效結點力，記為P e。6.5.2 結構整體載荷列陣的形成結構載荷列陣由所有單元的等效結點載荷列陣疊加得到。注意到疊加過程中相互聯接的單元之間存在大小相等方向相反的作用力和反作用力，它們之間相互抵消，因此，結構載荷列陣中只有與外載荷有關的結點有值。下面逐項進行討論。（1）集中力的等效載荷列陣逐點合成各單元的等效結點力，并按結點號碼的順序進行排列，組成結構的集中力等效載荷列陣，即 (5.139)上式中，單元e的集中力的等效結點力為(記單元結點局部編號為i,j,

22、m) (5.140)式中 （i, j, m） (5.141)式中，(Ni )c 、(Nj )c 、(Nm )c 為形函數在集中力作用點處的值。（2）表面力的等效載荷列陣把作用在單元邊界上的表面力移置到結點上，得到各單元的表面力的等效結點力。按照結點號碼的順序進行排列，逐個結點疊加合成后，組成結構表面力的等效載荷列陣，即 (5.142)式中， (5.143)由于作用在單元邊界上的內力在合成過程中已相互抵消，上式中的結點力只由作用在結構邊界上的表面力所引起。（3）體積力的等效載荷列陣與表面力類似，體積力的等效載荷列陣也是由單元體積力的等效結點力按結點號碼順序排列，在各結點處合成得到 (5.144)式中，單元e的體積力的等效結點力為 (5.145)6.5.3載荷移置與靜力等效關系上述基于形函數的載荷等效所得到的結果與按照靜力學的平行力分解原理得到的結果完全一致。例如，如圖3-19所示的單元e，在ij邊上作用有表面力。假設ij邊的長度為l，其上任一點P距結點i的距離為s。根據面積坐標的概念，有， ， (a)代入式(5.137)，求得單元表面力的等效結點力 (b)可見，求得的結果與按照靜力等效原理將表面力q向結點i及j分解所得到的分力完全相同。圖3-19 表面力等效示意再如，從圖3-20所示的單元e的A點處取體積微