有限元法與程序-平面問題6

第2章平面問題的有限元法上次課內(nèi)容回顧

平面矩形單元平面6節(jié)點三角形單元有限元法分析思路流程解綜合方程[K]{⊿}={P}求結構節(jié)點位移{⊿}計算結構內(nèi)力和應力系統(tǒng)分析(把單元剛度矩陣集合成結構剛度矩陣[K]形成等價節(jié)點荷載{P})離散（剖分）結構為若干單元單元分析(建立單元剛度矩陣[k]e形成單元等價節(jié)點力)§2.9平面有限元分析實施步驟與注意事項將要計算的彈性體劃分成三角形單元。對結點進行編號，列出結點坐標作為輸入信息。對單元進行編號，列出單元三個結點的號碼作為輸入信息。計算載荷的等效結點力，把等效結點力作為輸入信息。計算各單元的常數(shù)，再計算2Δ。以三角形常應變單元解平面問題為例，具體歩驟可歸納如下一、實施步驟§2.9平面有限元分析實施步驟與注意事項計算各單元的剛度矩陣。形成整體剛度矩陣。處理約束及消除剛體位移。解線性方程組，求結點位移。計算應力矩陣，再計算單元應力。根據(jù)需要計算主應力和主方向?！?.9平面有限元分析實施步驟與注意事項對稱性的利用。結點的選擇和單元劃分。結點的編號。單元結點i、j、m的次序。邊界條件處理和整體剛度矩陣的修正應力計算結果的整理二、注意事項§2.9平面有限元分析實施步驟與注意事項二、注意事項邊界條件的處理：劃0置1法邊界條件§2.9平面有限元分析實施步驟與注意事項二、注意事項邊界條件的處理：沖大數(shù)法§2.10平面問題的計算機程序一、有限元分析程序基本框圖開始輸入基本數(shù)據(jù)計算單元剛度矩陣形成總體剛度矩陣形成結點荷載向量引入約束條件求解方程組，輸出結點位移計算單元應力，輸出結果結束1、輸入基本數(shù)據(jù)（結構描述）：（1）控制數(shù)據(jù)：如結點總數(shù)、單元總數(shù)、約束條件總數(shù)等；（2）結點數(shù)據(jù)：如結點編號、結點坐標、約束條件等；（3）單元數(shù)據(jù)：如單元編號、單元結點序號、單元的材料特性、幾何特性等；（4）載荷數(shù)據(jù)：包括集中載荷、分布載荷等?！?.10平面問題的計算機程序一、有限元分析程序基本框圖開始輸入基本數(shù)據(jù)計算單元剛度矩陣形成總體剛度矩陣形成結點荷載向量引入約束條件求解方程組，輸出結點位移計算單元應力，輸出結果結束（1）各單元的bi,ci(i,j,m)

，面積A；（2）應變矩陣[B]，應力矩陣[S]；（3）單元剛度矩陣[k]；（4）單元等價載荷列向量[F]。2、單元分析

3、系統(tǒng)分析（1）整體剛度矩陣[K]的組裝；（2）整體載荷列陣{P}的形成；關鍵問題：[K]的存儲；§2.10平面問題的計算機程序一、有限元分析程序基本框圖開始輸入基本數(shù)據(jù)計算單元剛度矩陣形成總體剛度矩陣形成結點荷載向量引入約束條件求解方程組，輸出結點位移計算單元應力，輸出結果結束4、約束引入5、線性方程組求解求解方法常用：GAUSS消元法，LDLT、QR分解法等。其程序在一些專著中列出（例如：徐士良編。FORTRAN常用算法程序集。清華大學出版社）。在此不作詳細介紹，其方法參閱有關書籍。§2.10平面問題的計算機程序一、有限元分析程序基本框圖開始輸入基本數(shù)據(jù)計算單元剛度矩陣形成總體剛度矩陣形成結點荷載向量引入約束條件求解方程組，輸出結點位移計算單元應力，輸出結果結束6應力計算結果的整理計算結果包括位移和應力兩個方面。在位移方面，一般無須進行整理工作。應力結果則需要整理。通常認為計算出的應力是三角形單元形心處的應力。而相鄰單元之間的應力存在突變，甚至正、負符號都不相同。為了由計算結果推算出結構內(nèi)某一點的接接實際的應力，必須通過某種平均計算。通?？刹捎脙蓡卧骄ɑ蚶@結點平均法?！?.10平面問題的計算機程序一、有限元分析程序基本框圖開始輸入基本數(shù)據(jù)計算單元剛度矩陣形成總體剛度矩陣形成結點荷載向量引入約束條件求解方程組，輸出結點位移計算單元應力，輸出結果結束6應力計算結果的整理兩單元平均法：把兩個相鄰單元中的常應力加以平均，用來表示公共邊界中點處的應力。繞結點平均法：把環(huán)繞某一結點的各單元常應力加以平均，用以表示該結點的應力。在內(nèi)結點效果較好，而在邊界結點可能很差，一般改為應由內(nèi)結點的應力外推計算出來。§2.10平面問題的計算機程序二、總體剛度矩陣[K]的存儲全矩陣存貯法：不利于節(jié)省計算機的存貯空間，很少采用。對稱三角存貯法：存貯上三角或下三角元素。半帶寬存貯法：存貯上三角形（或下三角形）半帶寬以內(nèi)的元素。一維壓縮存貯法：半帶寬存貯中仍包含了許多零元素。存貯每一行的第一個非零元素到主對角線元素。等帶寬形式UBWUBW行號1→IR

→N→1列號JC行號1→IR→N→1JC-(IR-1)方陣形式1半帶寬存貯法方陣存貯和半帶寬存貯地址關系存貯方式行號列號方陣存貯IRJC等帶寬存貯IRJC-IR+1半帶寬計算：設結構單元網(wǎng)格中相鄰結點編號的最大差值是d，則最大半帶寬為UBW：結點編號：欲使最大半帶寬UBW最小，必須注意結點編號方法，使直接聯(lián)系的相鄰節(jié)點的最大點號差最小。例：計算下圖半帶寬。結點數(shù)N=91，總剛[K]中的元素總數(shù)為：（91×2）×（91×2）=33124最大半帶寬UBW=(7+1)×2=16，半帶寬存儲矩陣元素總數(shù)為182×16=2912，約方陣元素的8.8%。2變帶寬存貯（一維壓縮存貯）等帶寬存貯雖然已經(jīng)節(jié)省了不少內(nèi)存，但認真研究半帶寬內(nèi)的元素，還有相當數(shù)量的零元素。在平衡方程求解過程中，有些零元素只增加運算工作量而對計算結果不產(chǎn)生影響。如果這些零元素不存、不算，更能節(jié)省內(nèi)存和運算時間，采用變帶寬存貯可以實現(xiàn)（也稱一維數(shù)組存貯）。變帶寬存貯編程技巧要求較高，程序較長。對稱方陣形式的剛度矩陣[K]UBW=4頂線頂線以上零元素無須存貯，僅頂線以下元素。124610121618MAXA

22一維數(shù)組[A]存貯剛度矩陣[K]135911151721MAXA

22一維數(shù)組[A]存貯剛度矩陣[K]

變帶寬存貯：按列存貯方式。從左到右，逐列存放；對每一列，先存主對角線元素，然后由下而上順序存放，直到頂線下第一個元素為止。為避免混淆，我們把存貯[K]的一維數(shù)組稱為[A]。實現(xiàn)變帶寬存貯的關鍵問題是：總剛中元素Kij在一維數(shù)組A中的地址是什么？為此，需要知道主元Kii在A中的位置和相應列高hi。主元位置：采用一個一維數(shù)組MAXA存主元在A中位置。MAXA=[1,2,4,6,10,12,16,18，22]。列高hj：第j行的左帶寬。從第j列的主對角線元素起到該列上方第一個非零元素為止，所含元素的個數(shù)稱為第j列的列高，記為hj

；如果把第j列上方第1個非零元素的行號記為mj，則第j列的列高為

hj=j-mj

+1其實，hj就是第j行的左帶寬，因而必有

UBW=max(hj)

j=1,2,…,N利用節(jié)點位移信息數(shù)組

（去約束后節(jié)點位移自由度編碼），可容易地確定剛度矩陣[K]任何一列的列高。

主元在一維數(shù)組[A]中的地址數(shù)組MAXA的長度是[K]的行或列數(shù)加1（N+1）。

[K]的任何一個主對角元在一維數(shù)組A中的地址：第j列主對角線元素Kjj在一維數(shù)組A中的地址等于前(j-1)列的列高之和加1，即確定第j列列高的辦法是：從1號單元起，對所有單元逐個進行檢查。

MAXA(j)=h1+…+hj-2-hj-1+1=（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列的列高。一維數(shù)組A的總長度(S)，即剛度矩陣K按變帶寬存貯的總存貯量

S=MAXA(N+1)-MAXA(1)Ki,j在一維數(shù)組[A]中的地址

記Ki,j在一維數(shù)組A中的地址為AIJ。則由下圖可知，

AIJ=MAXA[J]+J–I其中，I=mj，mj+1，…，J