1、 第第2 2章章 平面問題的有限元法平面問題的有限元法上次課內容回顧上次課內容回顧 平面矩形單元平面矩形單元平面平面6節點三角形單元節點三角形單元有限元法分析思路流程有限元法分析思路流程解綜合方程解綜合方程K= P求結構節點位移求結構節點位移計算結構內力和應力計算結構內力和應力系統分析系統分析(把單元剛度矩陣集合成結構剛度矩陣把單元剛度矩陣集合成結構剛度矩陣K形成等價節點荷載形成等價節點荷載P )離散（剖分）結構離散（剖分）結構為若干單元為若干單元單元分析單元分析(建立單元剛度矩陣建立單元剛度矩陣ke形成單元等價節點力形成單元等價節點力) 2.9 平面有限元分析實施步驟與注意事項(1) 將要計

2、算的彈性體劃分成三角形單元。對結點進行編號，列出結點坐標作為輸入信息。(2) 對單元進行編號，列出單元三個結點的號碼作為輸入信息。 (3) 計算載荷的等效結點力，把等效結點力作為輸入信息。(4) 計算各單元的常數， 再計算2。iijjmmbcbcbc、 、 、 、以三角形常應變單元解平面問題為例，具體歩驟可以三角形常應變單元解平面問題為例，具體歩驟可歸納如下歸納如下 一、實施步驟一、實施步驟 2.9 平面有限元分析實施步驟與注意事項(6) 計算各單元的剛度矩陣 。(7) 形成整體剛度矩陣。(8) 處理約束及消除剛體位移。(9) 解線性方程組，求結點位移 。(10)計算應力矩陣，再計算單元應力。

3、根據需要計算主應力和主方向。 2.9 平面有限元分析實施步驟與注意事項(1) 對稱性的利用。(2) 結點的選擇和單元劃分。 (3) 結點的編號。(4) 單元結點i、j、m的次序。(5) 邊界條件處理和整體剛度矩陣的修正(6) 應力計算結果的整理二、注意事項二、注意事項 2.9 平面有限元分析實施步驟與注意事項二、注意事項二、注意事項44434241343332312423222114131211KKKKKKKKKKKKKKKK2211vuvu1122xyxyFFFF邊界條件的處理：劃邊界條件的處理：劃0置置1法法44422422000100000001KKKK2211vuvu112112333

4、2411433yyFKKFKK11u23u邊界邊界條件條件 2.9 平面有限元分析實施步驟與注意事項二、注意事項二、注意事項邊界條件的處理：沖大數法邊界條件的處理：沖大數法151511111112131411212223241515231323334333241424344210101010yyKuKKKKFvKKKKuKKKKKvKKKKF 2.10 平面問題的計算機程序一、有限元分析程序基本框圖一、有限元分析程序基本框圖開始輸入基本數據計算單元剛度矩陣形成總體剛度矩陣形成結點荷載向量引入約束條件求解方程組，輸出結點位移計算單元應力，輸出結果結束1、輸入基本數據（結構描述）、輸入基本數據（結

5、構描述）：（1）控制數據：如結點總數、單元總數、約束條件總數等；（2）結點數據：如結點編號、結點坐標、約束條件等；（3）單元數據：如單元編號、單元結點序號、單元的材料特性、幾何特性等；（4）載荷數據：包括集中載荷、分布載荷等。 2.10 平面問題的計算機程序一、有限元分析程序基本框圖一、有限元分析程序基本框圖開始輸入基本數據計算單元剛度矩陣形成總體剛度矩陣形成結點荷載向量引入約束條件求解方程組，輸出結點位移計算單元應力，輸出結果結束（1）各單元的bi,ci(i,j,m) ， 面積A；（2）應變矩陣B，應力矩陣S；（3）單元剛度矩陣k；（4）單元等價載荷列向量F。2、單元分析 3、系統分析（1）

6、整體剛度矩陣K的組裝；（2）整體載荷列陣P的形成；關鍵問題：K的存儲； 2.10 平面問題的計算機程序一、有限元分析程序基本框圖一、有限元分析程序基本框圖開始輸入基本數據計算單元剛度矩陣形成總體剛度矩陣形成結點荷載向量引入約束條件求解方程組，輸出結點位移計算單元應力，輸出結果結束4、約束引入5、線性方程組求解 求解方法常用：GAUSS消元法，LDLT、QR分解法等。其程序在一些專著中列出（例如：徐士良編。FORTRAN常用算法程序集。清華大學出版社）。在此不作詳細介紹，其方法參閱有關書籍。 2.10 平面問題的計算機程序一、有限元分析程序基本框圖一、有限元分析程序基本框圖開始輸入基本數據計算單

7、元剛度矩陣形成總體剛度矩陣形成結點荷載向量引入約束條件求解方程組，輸出結點位移計算單元應力，輸出結果結束6 應力計算結果的整理應力計算結果的整理 計算結果包括位移和應力兩個方面。在位移方面，一般無須進行整理工作。應力結果則需要整理。通常認為計算出的應力是三角形單元形心處的應力。而相鄰單元之間的應力存在突變，甚至正、負符號都不相同。為了由計算結果推算出結構內某一點的接接實際的應力，必須通過某種平均計算。通常可采用兩單元平均法或繞結點平均兩單元平均法或繞結點平均法法。 2.10 平面問題的計算機程序一、有限元分析程序基本框圖一、有限元分析程序基本框圖開始輸入基本數據計算單元剛度矩陣形成總體剛度矩陣

8、形成結點荷載向量引入約束條件求解方程組，輸出結點位移計算單元應力，輸出結果結束6 應力計算結果的整理應力計算結果的整理 兩單元平均法兩單元平均法：把兩個相鄰單元中的常應力加以平均，用來表示公共邊界中點處的應力。 繞結點平均法繞結點平均法：把環繞某一結點的各單元常應力加以平均，用以表示該結點的應力。在內結點效果較好，而在邊界結點可能很差，一般改為應由內結點的應力外推計算出來。 2.10 平面問題的計算機程序二、總體剛度矩陣二、總體剛度矩陣K的存儲的存儲全矩陣存貯法：不利于節省計算機的存貯空間，很少采用。對稱三角存貯法：存貯上三角或下三角元素。半帶寬存貯法半帶寬存貯法 ：存貯上三角形（或下三角形）

9、半帶寬以內的元素 。一維壓縮存貯法一維壓縮存貯法 ：半帶寬存貯中仍包含了許多零元素。存貯每一行的第一個非零元素到主對角線元素。 000000000000000000等帶寬形式等帶寬形式UBWUBW行行 號號1 IR N1列列 號號JC行行 號號1 IR N1JC-(IR-1) 00000000方陣形式方陣形式1 1 半帶寬存貯法半帶寬存貯法方陣存貯和半帶寬存貯地址關方陣存貯和半帶寬存貯地址關系系存貯方式存貯方式行號行號列號列號方陣存貯方陣存貯IRJC等帶寬存貯等帶寬存貯IRJC-IR+1u 半帶寬計算：設結構單元網格中相鄰結點編號的半帶寬計算：設結構單元網格中相鄰結點編號的最大差值是最大差值是

10、d，則最大半帶寬為，則最大半帶寬為UBW：為一個結點的自由度數ndfndfdUBW,) 1(u結點編號：欲使最大半帶寬結點編號：欲使最大半帶寬UBW最小，必須注最小，必須注意結點編號方法，使直接聯系的相鄰節點的最大點意結點編號方法，使直接聯系的相鄰節點的最大點號差最小。號差最小。 例：計算下圖半帶寬。結點數N=91，總剛K中的元素總數為：（912）（91 2 ）=33124最大半帶寬UBW=(7+1) 2=16，半帶寬存儲矩陣元素總數為182 16=2912，約方陣元素的8.8%。2 變帶寬存貯變帶寬存貯（一維壓縮存貯）（一維壓縮存貯） 等帶寬存貯雖然已經節省了不少內存，但認真等帶寬存貯雖然已

11、經節省了不少內存，但認真研究半帶寬內的元素，還有相當數量的零元素。在研究半帶寬內的元素，還有相當數量的零元素。在平衡方程求解過程中，有些零元素只增加運算工作平衡方程求解過程中，有些零元素只增加運算工作量而對計算結果不產生影響。如果這些零元素不存、量而對計算結果不產生影響。如果這些零元素不存、不算，更能節省內存和運算時間，采用變帶寬存貯不算，更能節省內存和運算時間，采用變帶寬存貯可以實現（也稱一維數組存貯）可以實現（也稱一維數組存貯） 。變帶寬存貯編程。變帶寬存貯編程技巧要求較高，程序較長。技巧要求較高，程序較長。88787767665856554645443634332322141211000

12、00000000000000KKKKKKKKKKKKKKKKKKKK對對 稱稱u方陣形式的剛度矩陣方陣形式的剛度矩陣KUBW=4頂頂 線線頂線以上零元素無須存貯，僅頂線以下元素。頂線以上零元素無須存貯，僅頂線以下元素。 )()()()()()()()()()()()()()()()()()()()()(181916201712211310141161574852931AAAAAAAAAAAAAAAAAAAAAA8878776766585655464544363433232214121100000000000000000KKKKKKKKKKKKKKKKKKKK124610121618MAXA 2

13、2u一維數組一維數組A存貯剛度矩陣存貯剛度矩陣K(1)(2)(6)(3)(4)(7)(5)(8)(12)(9)(10)(13)(11)(14)(18)(15)(16)(19)(17)(20)(21)AAAAAAAAAAAAAAAAAAAAAA8878776766585655464544363433232214121100000000000000000KKKKKKKKKKKKKKKKKKKK135911151721MAXA 22u一維數組一維數組A存貯剛度矩陣存貯剛度矩陣K 變帶寬存貯：按列存貯方式。變帶寬存貯：按列存貯方式。從左到右，逐列存從左到右，逐列存放；對每一列，先存主對角線元素，然后由

14、下而上順放；對每一列，先存主對角線元素，然后由下而上順序存放，直到頂線下第一個元素為止。序存放，直到頂線下第一個元素為止。為避免混淆，為避免混淆，我們把存貯我們把存貯K的一維數組稱為的一維數組稱為A。 實現變帶寬存貯的關鍵問題是：總剛中元素實現變帶寬存貯的關鍵問題是：總剛中元素Kij在在一維數組一維數組A中的地址是什么？為此，需要知道主元中的地址是什么？為此，需要知道主元Kii在在A中的位置和相應列高中的位置和相應列高hi。u主元位置：采用一個一維數組主元位置：采用一個一維數組MAXA存主元在存主元在A中中位置。位置。 MAXA =1,2,4,6,10,12,16,18，22。u列高列高hj：

15、第：第j行的左帶寬。行的左帶寬。 從第從第j列的主對角線元素起到該列上方第一個非零列的主對角線元素起到該列上方第一個非零元素為止，所含元素的個數稱為第元素為止，所含元素的個數稱為第j列的列的列高列高，記為，記為hj ；如果把第如果把第j列上方列上方第第1個非零元素的行號個非零元素的行號記為記為mj，則第，則第j列的列高為列的列高為 hj = j - mj + 1其實，其實，hj就是第就是第j行的左帶寬，因而必有行的左帶寬，因而必有 UBW= max(hj) j=1,2, ,N 利用節點位移信息數組利用節點位移信息數組 （去約束后節點位移自由（去約束后節點位移自由度編碼），可容易地確定剛度矩陣度

16、編碼），可容易地確定剛度矩陣K任何一列的列任何一列的列高。高。 u主元主元在一維數組在一維數組A中的地址中的地址 數組數組MAXA的長度是的長度是K的行或列數加的行或列數加1（N+1）。）。 K的任何一個主對角元在一維數組的任何一個主對角元在一維數組A中的地址：中的地址：第第j列主對角線元素列主對角線元素Kjj在一維數組在一維數組A中的地址等于中的地址等于前前(j-1)列的列高之和加列的列高之和加1，即，即 確定確定第第j列列列高的辦法是：從列高的辦法是：從1號單元起，對所有號單元起，對所有單元逐個進行檢查。單元逐個進行檢查。 MAXA(j) = h1 + hj-2 - hj-1 + 1 =（

17、h1 + hj-2 + 1）+ hj-1 = MAXA（j - 1）+ hj-1因為永遠有因為永遠有 MAXA(1)= 1， MAXA(2)= 2故計算主元地址的公式可寫為故計算主元地址的公式可寫為 MAXA(j+1)= MAXA(j)+ hj式中，式中， j = 2，3，N； hj剛度矩陣剛度矩陣K第第j列的列高。列的列高。一維數組一維數組A的總長度的總長度(S) ，即剛度矩陣，即剛度矩陣K按變帶寬存貯按變帶寬存貯的總存貯量的總存貯量 S = MAXA(N+1)- MAXA(1) uKi,j在一維數組在一維數組A中的地址中的地址 記記Ki,j在一維數組在一維數組A中的地址為中的地址為AIJ。則由下圖可知，。則由下圖可知， A I J = MAXA J + J I其中，其中，I = mj，mj+1，J。j列列 0000000000000000K第第i 行行頂線下頂線下第第j列列MAXA(J) (j-i)個個元素元素AIJ第第j 行行A中地址中地址Kj,jKi,jKmj,jK中地址中地址作業: （二選一）1) 寫出計算方法中LDL