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

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

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

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

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

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

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

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

變帶寬存貯：按列存貯方式。從左到右，逐列存放；對每一列，先存主對角線元素，然后由下而上順序存放，直到頂線下第一個元素為止。為避免混淆，我們把存貯[K]的一維數(shù)組稱為[A]。實現(xiàn)變帶寬存貯的關(guān)鍵問題是：總剛中元素Kij在一維數(shù)組A中的地址是什么？為此，需要知道主元Kii在A中的位置和相應(yīng)列高h(yuǎn)i。主元位置：采用一個一維數(shù)組MAXA存主元在A中位置。MAXA=[1,2,4,6,10,12,16,18，22]。列高h(yuǎn)j：第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因為永遠(yuǎn)有

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