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

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

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

，面積A；（2）應(yīng)變矩陣[B]，應(yīng)力矩陣[S]；（3）單元?jiǎng)偠染仃嘯k]；（4）單元等價(jià)載荷列向量[F]。2、單元分析

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

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

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

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

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

；如果把第j列上方第1個(gè)非零元素的行號(hào)記為mj，則第j列的列高為

hj=j-mj

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

UBW=max(hj)

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

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

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

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

MAXA(j)=h1+…+hj-2-hj-1+1=（h1+…+hj-2+1）+hj-1=MAXA（j-1）+hj-1因?yàn)橛肋h(yuǎn)有

MAXA(1)=1，MAXA(2)=2故計(jì)算主元地址的公式可寫為

MAXA(j+1)=MAXA(j)+hj式中，j=2，3，…，N；

hj——?jiǎng)偠染仃嘯K]第j列的列高。一維數(shù)組A的總長(zhǎng)度(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