第五章 線性三角形單元

1、三角形單元 3:212引 言 桿梁結(jié)構(gòu)：由于有自然的連接關(guān)系，可以憑一種直覺將其進(jìn)行自然的離散。 連續(xù)體：它的內(nèi)部沒有自然的連接節(jié)點，必須完全通過人工的方法進(jìn)行離散。三維問題平面問題平面應(yīng)力平面應(yīng)變平面問題平面應(yīng)力平面應(yīng)變離散3:213三節(jié)點平面三角形單元 x, u y, v 1 (x1, y1) (u1, v1) 2 (x2, y2) (u2, v2) 3 (x3, y3) (u3, v3) A fsx fsy 112233 euvuvuvd節(jié)點1的位移節(jié)點2的位移節(jié)點3的位移三節(jié)點三角形單元的位移函數(shù)可假設(shè)為：123456,u x yxyv x yxy“位移函數(shù)”也稱“位移模式”，是單元內(nèi)

2、部位移變化的數(shù)學(xué)表達(dá)式，是坐標(biāo)的函數(shù)。有限元分析必須事先給出（設(shè)定）位移函數(shù)。一般而論，位移函數(shù)選取會影響甚至嚴(yán)重影響計算結(jié)果的精度。彈性力學(xué)中，恰當(dāng)選取位移函數(shù)不是一件容易的事情。有限單元法中當(dāng)單元劃分得足夠小時，把位移函數(shù)設(shè)定為簡單的多項式也可得到相當(dāng)精確的結(jié)果。這正是有限單元法具有的重要優(yōu)勢之一。引入位移函數(shù)的概念：3:214平面三角形單元顯然，三角形三個節(jié)點的的位移可由下列方程給出，在各節(jié)點上的水平位移方程為： u1=1+2 x1 +3 y1 u2=1+2 x2 +3 y2 u3=1+2 x3 +3 y3解出11112223331,111xyuu x yxyxyuxyu11111222

3、23333111xyuxyuxyu3:215平面三角形單元11112322331111xyNNNxyxyxy假設(shè)求得1233223322311331133122112211()()21()()21()()2Nx yx yyy xxxyANx yx yyy xxxyANx yx yyyxxx yA11223311121xyAxyxy其中A是三角形的面積3:216平面三角形單元1 1223 3,v x yN vN vN v式中N1，N2和N3是坐標(biāo)的函數(shù)，反映了單元內(nèi)近似解的形態(tài)，稱為單元的形函數(shù)，數(shù)學(xué)上它反應(yīng)了由節(jié)點的場量對單元內(nèi)任意一點場量的插值，也叫做插值函數(shù)。三個函數(shù)其實描述的就是單元上近

4、似解的插值關(guān)系，它決定了近似解在單元上分布的形狀，所以稱它為形函數(shù)（shape function）。這里值得注意一下的是近似解，前面我們說過，假設(shè)位移模式是線性變化的，實際情況并不一定是線性變化的，所以我們通過所做假設(shè)得到的結(jié)果只能說是近似解，而不能說是精確解。 為什么叫形函數(shù)?1 12233,u x yN uN uN u同理( , )( , )ex yx yuNd3:217平面三角形單元12iiiiNab xc yA其中ijki = 1, 2, 3j = 2, 3, 1k = 3, 1, 2ijkmmjjjmmjiyxyxyxyxamjmiiyyyyb11mjjmixxxxc111211ii

5、jjmmxyAxyxy三角形的形函數(shù)可統(tǒng)一表示為：3:218形函數(shù)的性質(zhì)在單元任一點上三個形函數(shù)之和等于1(單位分解性)1. 三個形函數(shù)只有兩個是獨立的2. 當(dāng)三角形單元的三個結(jié)點的位移相等*ijmuuuu*( , )()iijjmmijmu x yN uN uN uNNNuu( , )( , )( , )1() ()()2ijmijmijmijmN x yN x yNx yaaabbb xcccyA第一列與它的代數(shù)余子式乘積之和第一列與第二列的代數(shù)余子式乘積之和第一列與第三列的代數(shù)余子式乘積之和2A0013:219形函數(shù)Ni 在節(jié)點i 上的值等于1，在其它節(jié)點上的值等于0。0),(0),(1

6、),(mmijjiiiiyxNyxNyxN0),(1),(0),(mmjjjjiijyxNyxNyxN1),(0),(0),(mmmjjmiimyxNyxNyxNNi =1ijmNj =1ijmNm =1ijm形函數(shù)的性質(zhì)3:2110在三角形單元的邊界ij上任一點（x，y），有: 形函數(shù)的性質(zhì)0),(),(1),(yxNxxxxyxNxxxxyxNmijijijiixxixjxyNi(xi，yi)j (xj，yj)m (xm，ym)Ni(x、y)1證jijixxxxyxN1),(ijijijiijijixxxxxxxxxxxxxxyxN1),( , )( , )( , )1ijmN x yN

7、x yNx y()jiiijimiimyyyyxxxxbyyxxcij方程3:2111形函數(shù)的性質(zhì)相鄰單元的位移在公共邊上是連續(xù)的ijpm0),(),(1),(yxNxxxxyxNxxxxyxNmijijijii形函數(shù)在單元上的面積分和邊界上的線積分公式為 ijijiAildlNAdxdyN213式中 為 邊的長度。 ijlijxxixjxyNi(xi，yi)j (xj，yj)m (xm，ym)Ni(x、y)1Ni =1ijm3:2112形函數(shù)的性質(zhì)完備性包含常應(yīng)變項和剛體位移項如果在勢能泛函中所出現(xiàn)的位移函數(shù)的最高階導(dǎo)數(shù)是m階，則選取的位移函數(shù)至少是m階完全多項式。協(xié)調(diào)性相鄰單元公共邊界保持

8、位移連續(xù)如果在勢能泛函中所出現(xiàn)的位移函數(shù)的最高階導(dǎo)數(shù)是m階，則位移函數(shù)在單元交界面上必須具有直至(m-1)階的連續(xù)導(dǎo)數(shù)，即Cm-1連續(xù)性。 如果在單元交界面上位移不連續(xù)，表現(xiàn)為當(dāng)結(jié)構(gòu)變形時將在相鄰單元間產(chǎn)生縫隙或重疊，這意味著將引起無限大的應(yīng)變，這時必然會發(fā)生交界面上的附加應(yīng)變能補充到系統(tǒng)的應(yīng)變能中去，有限元解就不可能收斂于真正解。收斂單元尺寸趨于零時，有限元解趨于真解3:2113形函數(shù)的性質(zhì)當(dāng)單元的位移函數(shù)滿足完備性要求時，稱單元是完備的（通常較容易滿足）。當(dāng)單元的位移函數(shù)滿足協(xié)調(diào)性要求時，稱單元是協(xié)調(diào)的。當(dāng)勢能泛函中位移函數(shù)的導(dǎo)數(shù)是2階時，要求位移函數(shù)在單元的交界面上具有C1或更高的連續(xù)性

9、，這時構(gòu)造單元的插值函數(shù)往往比較困難。在某些情況下，可以放松對協(xié)調(diào)性的要求，只要單元能夠通過分片試驗 (Patch test)，有限元分析的解答仍然可以收斂于正確的解。這種單元稱為非協(xié)調(diào)單元。分片試驗由B.M.Irons首先提出，已經(jīng)證明它給出了收斂性的充分條件。 3:2114單元應(yīng)變和應(yīng)力矩陣 0( , )( , )0 ( , )( , )( , )( , ) xyeexyxu x yx yx yx yx yv x yyyx uNdBd應(yīng)變矩陣123123 0 0 0 0( , ) 0 0 0 0 xNNNx yNNNyyx BN3:2115單元應(yīng)變和應(yīng)力矩陣123123123112233

10、0 0 01( , )0 0 0 2 bbbx ycccAcbcbcbBBBB312112233112233 0 0 00 0 0 bbbccccbcbcbBBB 由于 與x、y無關(guān)，都是常量，因此B矩陣也是常量。單元中任一點的應(yīng)變分量是B矩陣與單元節(jié)點位移的乘積，因而也都是常量。因此，這種單元被稱為常應(yīng)變單元。332211,cbcbcbA3:2116單元應(yīng)變和應(yīng)力矩陣21 0( , ) 1 0( , )11-0 0 2xxyyxyxyEx yx yD平面應(yīng)力：( , )( , )( , )( , )eex yx yx yx yDDBdSd 123123( , ) x y SDBD BBBSS

11、S2 2 (1)1-1- 22iiiiiiiibcEbcAcbSDB應(yīng)力矩陣平面應(yīng)變：用平面應(yīng)變彈性矩陣代入得到類似結(jié)果。3:2117單元應(yīng)變和應(yīng)力矩陣由于同一單元中的D、B矩陣都是常數(shù)矩陣，所以S矩陣也是常數(shù)矩陣。也就是說，三角形三節(jié)點單元內(nèi)的應(yīng)力分量也是常量。當(dāng)然，相鄰單元的E，A和bi、ci(i，j，m)一般不完全相同，因而具有不同的應(yīng)力，這就造成在相鄰單元的公共邊上存在著應(yīng)力突變現(xiàn)象。但是隨著網(wǎng)格的細(xì)分，這種突變將會迅速減小。3:2118單元分析 6 1ed 3 1 3 6B幾何關(guān)系位移函數(shù) 3 3D本構(gòu)關(guān)系 3 1平衡關(guān)系 6 1F 3 6SDB6 6?eK單元剛度矩陣3:2119單

12、元應(yīng)變能AAxyxyyyxxhdxdyhdxdyUT21)(21AAhdxdyhdxdyUDTT2121eAeAeehdxdyhdxdydDBBdDBdBdTTTT2121單元應(yīng)變能 U為：TTTT)(DDeBdijmxyh注意到彈性矩陣D的對稱性3:2120剛度矩陣引入剛度矩陣K：AehdxdyDBBKTeeeUdKdT21則：注意：hdxdy的實質(zhì)是任意的微體積dv，于是得 Ke的一般式：VeVdTDBBK3:2121單元外力功單元受到的外力一般包括體積力、表面力和集中力。自重屬于體積力范疇。表面力指作用在單元表面的分布載荷，如風(fēng)力、壓力，以及相鄰單元互相作用的內(nèi)力等。 ijmxyqVij

13、mxyqsijmxyfcVyVxVqqqsysxsqqqcycxcfff3:2122單元外力功（1）體積力所做的外力功ijmxyqVAVAVyVxVhdxdyhdxdyvquqWquTeuNdAVeVhdxdyWqNdTTVVeVVWdTTqNd3:2123單元外力功（2）面力所做的外力功lSlSySxShdlhdlvquqWquTeuNdlSeShdlWqNdTTtVeSWdTTqNdijmxyqsqs3:2124單元外力功（3）集中力所做的外力功cecWfdTijmxyfc當(dāng)結(jié)構(gòu)受到集中力時，通常在劃分單元網(wǎng)格時就把集中力的作用點設(shè)置為節(jié)點。于是單元集中力fc的勢能Vc為綜合以上諸式，單元

14、外力的總外力功V為 eecVVeceVeVeCSVttttWWWWfdfqNqNdfdqNdqNdTTTTTTTTTddddcVVettfqNqNfddTT3:2125系統(tǒng)勢能KddT121NeeUU擴(kuò)充疊加TT12UW d Kdd fT1NeeWWd fTeeeW d f擴(kuò)充疊加系統(tǒng)勢能eeeeUdKdT213:2126單元剛度矩陣的擴(kuò)充疊加T12310 0 0 0 0 00 0 0 0 0 00 0 01 0 0 020 0 eeeiiijimenneeejijjjmemimUKKKd KduuuuuKKKKK1231 00 0 0 0 0 0neenjmmuuuuuKmijmij單元編號

15、 ijmeeeeUdKdT213:2127單元等效節(jié)點載荷列陣的擴(kuò)充疊加T12310 0eieennjemWfd fuuuuuffmijTeeeW d f單元編號 ijm3:2128能量原理和系統(tǒng)平衡方程TT12UW d Kdd f系統(tǒng)勢能根據(jù)彈性力學(xué)能量原理：結(jié)構(gòu)處于穩(wěn)定平衡的必要和充分條件是總勢能有極小值。00fKdfKd 上式是從能量原理導(dǎo)出的系統(tǒng)平衡方程。這個方程表達(dá)了節(jié)點力與節(jié)點位移之間的關(guān)系。3:2129剛度矩陣單元剛度矩陣：VeVdTDBBK11122122 1 -1( )( ) -1 1 eeeeeeTeeeeelKKE Ax E Ax dxlKKKBBTT eeeiiijim

16、eeeejijjjmAeeemimjmmhdxdytAKKKKB DBB DBKKKKKK1D:2D:系統(tǒng)剛度矩陣：1eNeeKK彈性矩陣D的對稱性Ke對稱K 對稱3:2130剛度矩陣剛度矩陣K的詳細(xì)內(nèi)容為：ssNNjNiNNNjNjjjijjiNijiiiiNjiNjiNNKKKKKKKKKKKKKKKKKKKKKKKKKssssssssss21212122222211111211Ki、j是行列號，Ns為系統(tǒng)自由度數(shù)。3:2131剛度矩陣(1) 剛度矩陣中每個元素有明確的物理意義。例如，Kij表示當(dāng)節(jié)點位移中第j個元素為1(dj=1)其余元素為零時，引起的單元力中的第i個節(jié)點力fi。把平衡方

17、程寫開主對角線上元素Kii(i=1, Ns)恒為正值：位移和作用力同向ssssssssssssNjiNjiNNjNiNNNjNjjjijjiNijiiiiNjiNjifffffdddddKKKKKKKKKKKKKKKKKKKKKKKKK2121212121222222111112113:2132剛度矩陣(2) K的每一行或每一列元素之和為零以上式中第i行為例:iNiNjijiiiiifdKdKdKdKdKss2211f2i-1 =0 f2i =0f2j-1=0 f2j=0f2m-1 =0 f2m =0ijmxyijm11當(dāng)所有節(jié)點沿x向或y向都產(chǎn)生單位位移時，單元作平動運動，無應(yīng)變，也無應(yīng)力，

18、因而單元結(jié)點力為零（不含初應(yīng)力）。 所以有即，K的每一行元素之和為零。由于對稱性，每一列元素之和也為零。021siNijiiiiKKKKK3:2133剛度矩陣(3) 系統(tǒng)剛度矩陣是奇異矩陣（ 即K的行列式為零）(4) 系統(tǒng)剛度矩陣是常量矩陣系統(tǒng)的節(jié)點力和節(jié)點位移成線性關(guān)系是基于彈性理論的結(jié)果。剛度矩陣是在系統(tǒng)處于平衡狀態(tài)的前提下得出的。作用在它上面的外力必定是平衡力系。然而，研究系統(tǒng)平衡時沒有引入約束。承受平衡力系作用的無約束系統(tǒng)，其變形是確定的，但位移不是確定的。所以出現(xiàn)性質(zhì)（2）中的“平動問題”，即可以發(fā)生任意的剛體運動。從數(shù)學(xué)上講，系統(tǒng)平衡方程的解不是唯一的或不能確定的。由此，系統(tǒng)剛度矩

19、陣一定是奇異的。單元剛度矩陣也一定是奇異的。3:2134位移邊界條件的處理系統(tǒng)剛度矩陣是奇異矩陣，其物理原因是結(jié)構(gòu)缺少剛性位移的約束，實際的工程結(jié)構(gòu)都受有足夠的支承約束，排除了發(fā)生任何剛體位移的可能性，因此，必須引入位移約束。有限元中，位移約束都設(shè)置在節(jié)點處。這里，只討論剛性約束情況，即被約束的位移分量為零。設(shè)討論的結(jié)構(gòu)有Nn個節(jié)點，每個節(jié)點有ndf個自由度。則系統(tǒng)的總自由度為Ns，且dfnsnNN節(jié)點總位移列向量d中共包含Ns個分量T1sNidddd3:2135為了引入位移約束，把節(jié)點總位移列向量d分成兩部分。一部分是不受約束的位移分量，記為df。另一部分是受剛性約束的位移分量，記為dr。不

20、失一般性，設(shè)1N號位移分量是不受約束的；N+1N+Nr共Nr個分量是受剛性約束的。即：rNNNNNdddddd2121,rfdd位移邊界條件的處理3:2136位移邊界條件的處理 111rsNNNrfddd顯然不受約束的節(jié)點位移的總數(shù)N為 N=Ns-Nr 對方程中的剛度矩陣K和節(jié)點荷載向量列陣f也作相應(yīng)分割，則得到rfrfrrrffrffffddKKKK式中，ff是已知力邊界，fr是約束反力。rfddd3:2137位移邊界條件的處理按矩陣乘法規(guī)則得frfrffffdKdK每個受剛性支承約束的位移分量都等于零，即0dr從而得到rfrfrrrffrffffddKKKKrrrrfrffdKdKffff

21、fdKrfrffdK3:2138位移邊界條件的處理 Kff 引入約束后的約化的系統(tǒng)剛度矩陣。 這是一個非奇異 矩陣，它的逆矩陣Kff-1是存在的。fffffdK引入約束后的約化的系統(tǒng)平衡方程在分析計算時，從無約束的系統(tǒng)剛度矩陣K中刪去與受約束位移號對應(yīng)的行和列，再將矩陣壓縮排列成NN階方陣，即為約化后的結(jié)構(gòu)剛度矩陣。rfrfrrrffrffffddKKKK3:2139位移邊界條件的處理ssssssssssssNjiNjiNNjNiNNNjNjjjijjiNijiiiiNjiNjifffffdddddKKKKKKKKKKKKKKKKKKKKKKKKK21212121212222221111121

22、100010iiud iiuf 顯然iNiuddds0101iiud 3:2140位移邊界條件的處理111211111212222222121212sssssssssssijNijNiiiiijiNiijjjjijjjNNNNN iN jN NKKKKKdfKKKKKdfKKKKKdfdKKKKKdKKKKKsjNffiiud iiiiuKf顯然iiiNiNiiiiuKdKdKdKss11iiud iiiiiiuKdKiiK3:2141節(jié)點位移和約束反力f1ffffKd通過求解平衡方程即可解出全部未知的節(jié)點位移：把解出的d代入未經(jīng)修改的平衡方程，即可得到約束反力：關(guān)于上述方程的解算方法，一般不

23、采用求逆的方法求解，而是直接采用高斯消元法等求解線性方程組的方法求解求解。施加邊界條件后，得到修改后的平衡方程fdK(未約化的)fffffdK(約化的)或fKd1fKd 或rfrffdK3:2142單元應(yīng)變和應(yīng)力( , )( , )xyexyx yx yBd( , )( , )( , )xyexyx yx yx yDDBd根據(jù)三角形節(jié)點的位移，求出單元應(yīng)力應(yīng)變?yōu)閕jk如何求系統(tǒng)應(yīng)變能和節(jié)點應(yīng)力？3:2143有限元解的收斂性由于在有限元計算中引入了結(jié)構(gòu)離散和位移模式，導(dǎo)致有限元計算結(jié)果和真實解的偏差。單元劃分越小、位移模式取得越接近真實變形，解答越收斂于真實解。當(dāng)單元的位移模式采用解析的位移解時

24、，有限元的計算結(jié)果和解析結(jié)果是相同的。然而，許多情況無真實的位移模式可以借用，只能尋求其近似函數(shù)，不可避免帶來計算精度問題。實踐證明：只要位移模式滿足單元的完備性準(zhǔn)則和協(xié)調(diào)性條件，就保證了有限元的解答收斂于真實解。3:2144有限元計算過程框圖剖分結(jié)構(gòu)為有限個單元，對節(jié)點、單元編號構(gòu)建單元剛度矩陣和單元等價節(jié)點力向量組裝系統(tǒng)剛度矩陣并引入約束組裝整體等價節(jié)點荷載向量從節(jié)點平衡方程解未知節(jié)點位移計算結(jié)構(gòu)應(yīng)變、應(yīng)力3:2145解綜合方程Kd=f得結(jié)構(gòu)節(jié)點位移d從d中找單元位移de用公式用公式=Bde和和 =D，計算應(yīng)力計算應(yīng)力應(yīng)變應(yīng)變把單元剛度矩陣組裝成系統(tǒng)剛度矩陣K離散結(jié)構(gòu)為若干單元建立單元剛度

25、矩陣Ke形成等價節(jié)點荷載f形成單元等價節(jié)點力有限元計算流程圖3:2146關(guān)于三角形單元形函數(shù)的一點補充12223322332331111()()() 221xyAxyx yx yyy xxxyxy11eALA 2-3-P: 同樣, 3-1-PA2 1-2-PA322eALA33eALA i,1 j, 2 m, 3 x y P A1 面積坐標(biāo)的定義 在三角形內(nèi)任意一點P定義 123 1LLL3:2147關(guān)于三角形單元形函數(shù)的一點補充 i,1 j, 2 m, 3 x y P A1 面積坐標(biāo)與形函數(shù)的關(guān)系1111221121212ijjiiimmiiiiiijjjjjjmmmmmmxyxyab xc

26、 yxyLab xc yNLab xc yNLab xcyN1iijjmmiijjmmijmxx Lx Lx Lyy Ly Ly LLLL 面積坐標(biāo)與直角坐標(biāo)的關(guān)系3:21(0,0)A(0, 2)B(4,0)D(4,1)C1F q 考慮一個平面應(yīng)力問題如圖所示，假設(shè)厚度考慮一個平面應(yīng)力問題如圖所示，假設(shè)厚度h=1，材料為各項，材料為各項同性，楊氏模量為同性，楊氏模量為E=1，泊松比為，泊松比為=0，相關(guān)力和位移邊界條，相關(guān)力和位移邊界條件如圖中所示，問題左端為固定約束。試用兩個三角形單元件如圖中所示，問題左端為固定約束。試用兩個三角形單元分析此問題，三角形單元的網(wǎng)格劃分如圖所示。試求問題各分析

27、此問題，三角形單元的網(wǎng)格劃分如圖所示。試求問題各節(jié)點位移節(jié)點位移u、v和應(yīng)力和應(yīng)力x，y和和xy。1243 1 2例題3:2112323311212332132111223332231331211200000010000002ebbbyyyyyycccxxxxxxAcbcbcbxxyyxxyyxxyyB對于三角形單元，其對于三角形單元，其B B矩陣的表達(dá)式為：矩陣的表達(dá)式為：23322313211()()()2eAx yx yyy xxxy1243 1 2123TeehAkB DB2101001001011000.5002EvDFor e=1:1922472216.501215.5204020

28、141024116722492215.501216.5ek110201010400048410241eB1, (2)3, (1)2, (3)1, (2)2, (3)3, (1)例題3:2112323311212332132111223332231331211200000010000002ebbbyyyyyycccxxxxxxAcbcbcbxxyyxxyyxxyyB23322313211()()()2eAx yx yyy xxxy1243 1 2123TeehAkB DB2101001001011000.5002EvDFor e=2:216016404032032001601842414324330116002020404101ek200101010404004404101eB1, (3)3, (1)2, (4)1, (3)2, (4)3, (1)對于三角形單元，其對于三角形單元，其B B矩陣的表達(dá)式為：矩陣的表達(dá)式為：例題3:21組裝剛度矩陣 12112722424217.5215.5410172922400215.5216.50100124202001641641410340