1、湖南科技大學土木工程學院湖南科技大學土木工程學院謝獻忠謝獻忠2014年年8月月 “用兵之法，十則圍之，倍則分之，少則能逃之。” “古之所謂善戰者，勝于易勝者也。” 懸臂梁力學模型懸臂梁有限元模型 早在四十年代初，Courant就提出了有限單元法的基本思想，但一直沒有引起人們的足夠重視。直到五十年代中期，Turner、Clough才利用這種思想對航空工程中的飛機結構進行分析。這種處理問題的思路，在1960年以后被廣泛用于求解彈性力學的平面問題，并首次提出了（Clough）“有限單元法”這一術語。之后，隨著電子計算機的飛速發展，有限單元法如虎添翼，經過50多年的發展，目前國內外已有許多大型通用的有

2、限元分析程序可供使用。 工程技術人員進行結構分析時的主要任務就是設法將復雜的工程實際問題加以簡化、建立合理的計算力學模型，然后再按所選用的程序的要求，準備好全部所需的數據和信息，運用計算機進行求解，最后再檢查計算結果的合理性，事實上，現在許多大型有限元分析軟件都已配備了功能很強的前后置處理程序，并已出現了將人工智能技術引入有限元分析軟件，形成了比較完善的專家系統，逐步實現了有限元分析的智能化。 幾十年來，有限元單元法已在各個工程領域得到了廣泛的應用，相應的大型軟件已成為現代工程設計中一個重要的不可缺少的計算工具。 平衡問題 穩定問題與動力問題 彈性問題 彈塑性與粘性問題、土力學與巖石力學 結構

3、計算 結構優化設計問題 固體力學 流體力學、流固耦合問題 、熱傳導 與熱應力、電磁場、建筑聲學與噪聲問題 平面應變問題平面應變問題Pxyxy平面應力問題平面應力問題一、離散化一、離散化 平面三角形單元平面三角形單元 ui (Ui ) um (Um ) uj (Uj ) vj (Vj ) vi (Vi ) um (Um ) j i m x y o xy1432u1u4u3u2v1v4v3v2 平面四邊形單元平面四邊形單元二、位移模式二、位移模式 由于在彈性體內，各點的位移變化情況非常復雜，很難在整個彈性體內選取一個恰當的位移函數來表示位移的復雜變化，但是如果將整個區域分割成許多小單元，那么在每個

4、單元的局部范圍內就可以采用比較簡單的函數來近似地表示單元的真實位移，將各單元的位移函數連起來，便可近似地表示整個區域的真實位移函數。這種化繁為簡、聯合局部逼近整體的思想，正是有限單元法的精妙之處。 uxyvxy123456 基于上述思想，我們可以選擇一種單元位移模式： ui (Ui ) um (Um ) uj (Uj ) vj (Vj ) vi (Vi ) um (Um ) j i m x y o 三、三、 單元剛度方程單元剛度方程 TmmjjiieVUVUVUR ui (Ui ) um (Um ) uj (Uj ) vj (Vj ) vi (Vi ) um (Um ) j i m x y o

5、 Tmmjjiievuvuvu eeekR四、四、 單元選用單元選用 彈性平面問題：二維實體單元，PLANE42、82等 彈性空間問題：三維實體單元，SOLID45、65等 平面桿件問題：二維桿單元、梁單元，LINK1、BEAM3等 空間桿件問題：三維桿單元、梁單元，LINK8、BEAM4等 板 殼 問 題 ：殼單元，SHELL63等 索結構問題 ：索單元，LINK等BEAM3PPPPP1234568791011121314iijxjuiujMjMiYiiYjjXjXivivjLINK1LINK8Qyj(vj) BEAM4ijMyi(yi)Qyi(vi)yMxi( xi)Ni(ui)Qxi (

6、Wi)Mzi( zi)zNj(uj)Mzj( zj)Qxj (Wj)Mxj( xj)xMyj(yj)BEAM188SHELL63xyxy1432u1u4u3u2v1v4v3v2 PLANE42BEAM188LINK10LINK10LINK8LINK8找形分析、應力剛化、預應力模態分析找形分析、應力剛化、預應力模態分析 eeekR ui (Ui ) um (Um ) uj (Uj ) vj (Vj ) vi (Vi ) um (Um ) j i m x y o 下面以具體實例闡述整體剛度矩陣的組集規律。2N2N2m2m2m2m由于對稱性，由于對稱性，只取四分之一只取四分之一作為研究對象作為研究對

7、象xy1N1234(1)(2)（一）整體編碼與局部編碼的對應關系（一）整體編碼與局部編碼的對應關系 劃分單元以后，建立單元和節點的整體和局部編碼，并建立兩種編碼的對應關系。局部局部編碼編碼整體編碼整體編碼單元單元(1)單元單元(2)ijm124234xy1N1234(1)(2)（二）單元剛度方程及其擴充形式（二）單元剛度方程及其擴充形式單元單元(1): i, j, m 1, 2, 4 xy1N1234(1)(2)(1)(1)432144424124222114121142142144424124222114121142100000000kkkkkkkkkRRRkkkkkkkkkRRRkkkkk

8、kkkkRRRmjimmmjmijmjjjiimijiimji單元單元(2)(2): i, j, m 2, 3, 4 432144434234333224232243243244434234333224232243200000000kkkkkkkkkRRRkkkkkkkkkRRRkkkkkkkkkRRRmjimmmjmijmjjjiimijiimji(2)(2)xy1N1234(1)(2)（三）總剛方程（三）總剛方程 將擴充后的單元剛度方程左右兩邊相加，得總剛方程：4321)21 (44)2(43)21 (42) 1 (41)2(34)2(33)2(32)21 (24)2(23)21 (22)

9、 1 (21) 1 (14) 1 (12) 1 (11)21 (4)2(3)21 (2) 1 (100kkkkkkkkkkkkkkRRRRxy1N1234(1)(2) kR tdxdypftdsqfGfRTTTceTe)( 由總剛方程可知，總的節點載荷列陣R是由彈性體的全部單元的等效節點力集合而成，而其中單元的等效節點力Re 則是由作用在單元上的集中力、表面力和體積力分別移置到節點上，再逐點加以合成求得。根據虛位移原理，等效節點力的大小，應按其所做的功與作用在單元上的三種力在任何虛位移上所做的功相等這一原則來確定。即 上式中等號的左邊表示單元的等效節點力Re 所做的虛功；等號右邊的第一項是集中

10、力G所做的虛功、第二項的積分是沿著單元的邊界進行，表示面力q所(a) ) ()()(tdxdypNtdsqNGNRTTTcTeeTe做的虛功、第三項的積分則是遍及整個單元，表示體積力p所做的虛功；t為單元的厚度，假定為常量。將上節中的 (c) 式代入上式，并注意到節點虛位移列陣* e 中的元素都是常量，這就可以把( *e ) T 提到積分號的外面，于是有注意到( * e ) T 的任意性，有 tdxdypNtdsqNGNRTTTce 下面我們將逐項進行討論。 A：按虛功相等原理移置單元載荷。：按虛功相等原理移置單元載荷。集中力的等效載荷列陣 y0 xijmiyRixRmxRmyRjyRjxRC

11、yGxGG集中力 （原結構上的載荷） yxcmmjjiimymxjyjxiyixeGGNNNNNNRRRRRRR000000其中： 為形函數在集中力作用點處的值。 cmcjciNNN， 表面力的等效載荷列陣如下圖所示的單元e，在ij邊上作用有表面力。假設ij邊的長度為l，其上任一點P距節點i的距離為s。根據形函數的性質，有lsNi1lsNj0mN 0100llmjiemejeietdsqlstdsqlstdsqNtdsqNtdsqNRRRRsldseemijPqtdsmijRjeRie 體積力的等效載荷列陣設在單元ijm上受有分布體力 。y0 xijmiyRixRmxRmyRjyRjxReYX

12、pp tdxdypNtdxdypNtdxdypNRRRRmjiemejeie dxdyYXNNNNNNtRRRRRRRmmjjiimymxjyjxiyixe000000或當單元體是均質、等厚、比重為 時，則X=0，Y= - 故有：10101030000000tdxdyNNNNNNtRmmjjiie單元面積B：按靜力等效原理移置單元載荷。：按靜力等效原理移置單元載荷。一、體力的移置： 303030WWWRey0 xijmcwuP 01010130303032cTuuuextPPPRcx密度gtW eR如果任意三角形單元ijk的重心c上受有自重 則按剛體靜力等效原理可把W直接移置到i，j，m三個節

13、點上而組成 ：2cuxtP如果單元的重心c受有慣性力Pu作用，且 ，則Pu移置到i，j，m節點上的等效結點力為：式中： 旋轉角速度 是單元重心處x坐標。二、面力的移置y0 xijmlq2ql2ql已知在ij邊受有面力q，則移置到i、j結點上的等效節點力為： TeqlR0010102y0 xijml0q20lq3l TelqR0320310020當某一邊上有三角形分布的面力時，可由剛體靜力等效直接寫出三、集中力的移置如集中力G做用于其一邊界上如圖：先將G分解為 ， 后，分別按線段的比例把 和 分別移置到i，j兩點上。即：xPyPxPyP 001122llPllPllPllPRyxyxe0 xij

14、mG2lxPyP1lly下面我們通過一個只有四個方程的簡單例子來說明。4321432144434241343332312423222114131211RRRRuuuuKKKKKKKKKKKKKKKK總剛方程在前面討論整體剛度矩陣時，已經提到，整體剛度矩陣的奇異性可以通過考慮邊界約束條件來排除彈性體的剛體位移，以達到求解的目的。3431414332312121432144422422000100000001KKRKKRuuuuKKKK假定該系統中節點位移u1 和u3分別已知為當引入這些節點的已知位移之后，總剛方程就變成為然后，就用這組維數不變的方程來求解所有的節點位移。顯然，其解答仍為原方程的解

15、答。 u1 = 1 ， u3 = 3方案(1)-置零法4321432144434241343332312423222114131211RRRRuuuuKKKKKKKKKKKKKKKK415333215111432144434241341533323124232221141312151110101010RKRKuuuuKKKKKKKKKKKKKKKK15111414313212115111010KuKuKuKuK事實上，該方程組的第一個方程為方案(2)-乘大數法4321432144434241343332312423222114131211RRRRuuuuKKKKKKKKKKKKKKKK假定該系

16、統中節點位移u1 和u3分別已知為當引入這些節點的已知位移之后，總剛方程就變成為u1 = 1 ， u3 = 3方案(3)-方程縮減法假定該系統中節點位移u1 和u3分別已知為當引入這些節點的已知位移之后，總剛方程就變成為u1 = 0 ， u3 = 04321432144434241343332312423222114131211RRRRuuuuKKKKKKKKKKKKKKKK424244422422RRuuKKKK314234321412RRuuKKKK求位移求反力可以證明，對于一個給定的位移模式，其剛度系數的數值比精確值要大。所以，在給定的載荷之下，有限元計算模型的變形將比實際結構的變形小。

17、因而，當單元網格分得越來越細時，位移的近似解將由下方收斂于精確解，即得到真實解的下界。 對于一個數值計算方法，一般總是希望隨著網格的逐步細分所得到的解答能夠收斂于問題的精確解。根據前面的分析，我們知道，在有限元分析中，一旦確定了單元的形狀之后，位移模式的選擇將是非常關鍵的。由于載荷的移置、應力矩陣和剛度矩陣的建立等等，都依賴于單元的位移模式，所以，如果所選擇的位移模式與真實的位移分布有很大的差別，那么就很難獲得良好的數值解。 位移模式必須能包含單元的常應變。每個單元的應變一般都是包含著兩個部分：一部分是與該單元中各點的坐標位置有關的應變（即所謂各點的變應變）；另一部分是與位置坐標無關的應變（即

18、所謂的常應變）。從物理意義上看，為了保證解答的收斂性，要求位移模式必須滿足以下三個條件，即 位移模式必須包含單元的剛體位移。也就是說，當節點位移是由某個剛體位移所引起時，彈性體內將不會產生應變。所以，位移模式不但要具有描述單元本身形變的能力，而且還要具有描述由于其它單元形變而通過節點位移引起單元剛體位移的能力。 xvv, yuu00剛體位移: 位移模式在單元內要連續、且在相鄰單元之間的位移必須協調。當選擇多項式來構成位移模式時，單元內的連續性要求總是得到滿足的，單元間的位移協調性，就是要求單元之間既不會出現開裂也不會出現重疊的現象。通常，當單元交界面上的位移取決于該交界面上節點的位移時，就可以

19、保證位移的協調性。 當單元尺寸無限縮小時，每個單元中的應變應該趨于常量。因此，在位移模式中必須包含有這些常應變，否則就不可能使數值解收斂于正確解。很顯然，在前面的位移模式(b)中，與2、3、5、6 有關的線性項就是提供單元中的常應變的。 在有限單元法中，把能夠滿足條件1和2的單元，稱為完備單元；滿足條件3的單元，叫做協調單元或保續單元。前面討論過的三角形單元和矩形單元，均能同時滿足上述三個條件，因此都屬于完備的協調單元。在某些梁、板及殼體分析中，要使單元滿足條件3比較困難，所以實踐中有時也出現一些只滿足條件1和2的單元，其收斂性往往也能夠令人滿意特別是放松條件3的單元，即完備而不協調的單元，已

20、獲得了很多成功的應用。不協調單元一般不象協調單元那樣剛硬（即比較柔軟），因此有可能會比協調單元收斂得快。在選擇多項式作為單元的位移模式時，其階次的確定，要考慮解答的收斂性，即單元的完備性和協調性要求。實踐證明，這兩項確實是所要考慮的重要因素，但并不是唯一的因素。選擇多項式位移模式階次時，需要考慮的另一個因素是，所選的模式應該與局部坐標系的方位無關，這一性質稱為幾何各向同性。對于線性多項式，各向同性的要求通常就等價于位移模式必須包含常應變狀態。對于高次位移模式，就是不應該有一個偏惠的坐標方向，也就是位移形式不應該隨局部坐標的更換而改變。經驗證明，實現幾何各向同性的一種有效方法是，根據巴斯卡三角形

21、來選擇二維多項式的各項。在二維多項式中，如果包含有對稱軸一邊的某一項，那么就必須同時包含有另一邊的對稱項。選擇多項式位移模式時，還應該要考慮的一個因素就是，多項式中的項數必須等于或稍大于單元邊界上的外節點的自由度數。通常是取項數與單元的外節點的自由度數相等，取過多的項數是不恰當的。122322343223454322345xyxxyyxx yxyyxx yx yxyyxx yx yx yxyy常數項線性項二次項三次項四次項五次項對稱軸巴斯卡三角形有限元法求解力學問題的具體步驟:根據工程實際情況確定力學模型，并按一定比例繪制結構圖、注明尺寸、載荷和約束情況等。將計算對象進行離散化，并對節點、單元

22、進行編號。確定全部節點的坐標值。 計算等效節點載荷。 計算單元剛度矩陣。 組集整體剛度矩陣。 處理約束，消除剛體位移。 求解線性方程組，得到節點位移。 計算單元應力。 整理計算結果。 如圖所示為一厚度t=1cm的均質正方形薄板，上下受均勻拉力q=106N/m，材料彈性模量為E，泊松比 ，不記自重，試用有限元法求其應力分量。3/1單元劃分單元劃分2myxq=106N/m正方形板正方形板123421xy2q2q.力學模型的確定力學模型的確定按平面應力問題處理，考慮到結構和載荷的對稱性，可取結構的1/4來研究。該1/4結構被離散為兩個三角形單元。解：解：011jimimjmjiyybyybyyb2.

23、編碼與坐標編碼與坐標.求單元的剛度矩陣求單元的剛度矩陣1) 計算單元的b、c常數及單元面積單元（單元（i、j、m1，2，3）節點節點坐標坐標xy123421xy2q2q局部局部編碼編碼整體編碼整體編碼單元單元1單元單元2ijm1233412) 計算單元的剛度矩陣 先計算用到的常數211101jimimjmjixxcxxcxxc16009)1 (431212EEt 31001160091131000131103110310131003111160091EEKii代入可得:3432323416009;03131016009;313131116009111EKEKEKjjimij1003116009;13131311600911EKEKmmjm單元單元1的剛度矩陣為的剛度矩陣為: 103113134313132