有限元程序設(shè)計(jì)方法PPT

1、第2章 有限元程序設(shè)計(jì)方法2.1 程序基本框圖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ù)：包括集中載荷、分布載荷等。開始輸入基本數(shù)據(jù)計(jì)算單元?jiǎng)偠染仃囆纬煽傮w剛度矩陣形成結(jié)點(diǎn)荷載向量引入約束條件求解方程組，輸出結(jié)點(diǎn)位移計(jì)算單元應(yīng)力，輸出結(jié)果結(jié)束2、單元分析 （1）各單元的bi,ci(i,j,m) ， 面積A；（2）應(yīng)變矩陣B，應(yīng)力矩陣S；（3）單元?jiǎng)偠染仃噆；（4）單元等價(jià)載荷列向量F。開始輸入基本數(shù)據(jù)計(jì)算單元?jiǎng)偠染仃囆纬?/p>

2、總體剛度矩陣形成結(jié)點(diǎn)荷載向量引入約束條件求解方程組，輸出結(jié)點(diǎn)位移計(jì)算單元應(yīng)力，輸出結(jié)果結(jié)束3、系統(tǒng)分析（1）整體剛度矩陣K的組裝；（2）整體載荷列陣P的形成；K的存儲(chǔ)；約束引入；求解總剛存貯全矩陣存貯法：不利于節(jié)省計(jì)算機(jī)的存貯空間，很少采用。Ki,j對(duì)稱三角存貯法：存貯上三角或下三角元素。半帶寬存貯法 ：存貯上三角形（或下三角形）半帶寬以內(nèi)的元素 。一維壓縮存貯法 ：半帶寬存貯中仍包含了許多零元素。存貯每一行的第一個(gè)非零元素到主對(duì)角線元素。 等帶寬形式UBWUBW行 號(hào)1 IR N1列 號(hào)JC行 號(hào)1 IR N1JC-(IR-1)方陣形式（1）半帶寬存貯法方陣存貯和半帶寬存貯地址關(guān)系存貯方式行

3、號(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ù)為： 82（912）（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)算工

4、作量而對(duì)計(jì)算結(jié)果不產(chǎn)生影響。如果這些零元素不存、不算，更能節(jié)省內(nèi)存和運(yùn)算時(shí)間，采用變帶寬存貯可以實(shí)現(xiàn)（也稱一維數(shù)組存貯） 。變帶寬存貯編程技巧要求較高，程序較長。對(duì) 稱方陣形式的剛度矩陣KUBW=4頂 線頂線以上零元素?zé)o須存貯，僅頂線以下元素。124610121618MAXA 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)j 。

5、 ( hj 列號(hào)j-第1個(gè)非零元素的行號(hào)mj 1)主元位置：采用一個(gè)一維數(shù)組MAXA存主元在A中位置。 MAXA =1,2,4,6,10,12,16,18，22。主元在一維數(shù)組A中的地址MAXA(J) 數(shù)組MAXA的長度是K的行或列數(shù)加1（N+1）。 K的任何一個(gè)主對(duì)角元在一維數(shù)組A中的地址：第j列主對(duì)角線元素Kjj在一維數(shù)組A中的地址等于前(j-1)列的列高之和加1，即主元地址： MAXA(j+1)= MAXA(j)+ hj （2-1）式中， j = 2，3，N； hj剛度矩陣K第j列的列高。 MAXA(1)= 1， MAXA(2)= 2Ki,j在一維數(shù)組A中的地址 AIJ記Ki,j在一維數(shù)

6、組A中的地址為AIJ。則由下圖可知， A I J = MAXA J + J I （2-2） 其中，I = mj，mj+1，J。圖5-12j列第i 行頂線下第j列MAXA(J)(j-i)個(gè)元素AIJ第j 行A中地址Kj,jKi,jKmj,jK中地址mj4、引入約束條件手算時(shí)采用去行列法 ，而計(jì)算機(jī)編程時(shí)采用乘大數(shù)法。即：指定結(jié)點(diǎn)位移對(duì)應(yīng)的主對(duì)角元素乘上一個(gè)大數(shù) ，同時(shí)將P中對(duì)應(yīng)元素?fù)Q為結(jié)點(diǎn)位移指定值與擴(kuò)大了的主對(duì)角線元素的乘積 。5、線性方程組求解求解方法常用：GAUSS消元法，QR分解法等。其程序在一些專著中列出（例如見：徐士良編。FORTRAN常用算法程序集。清華大學(xué)出版社）。在此不作詳細(xì)介

7、紹。6、單元應(yīng)力節(jié)點(diǎn)位移求單元應(yīng)力。首先整體節(jié)點(diǎn)位移變換成單元節(jié)點(diǎn)位移，然后再用物理方程求單元應(yīng)力。2.2 提高計(jì)算精度的方法 （1） 計(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)平均法。平均法整理單元應(yīng)力邊界上：采用兩單元平均法，即把兩個(gè)相鄰單元中的常應(yīng)力加以平均，用來表示公共邊界中點(diǎn)處的應(yīng)力。節(jié)點(diǎn)上：采用繞結(jié)點(diǎn)平均法，即把環(huán)繞某一結(jié)點(diǎn)的

8、各單元常應(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ì)算出來。（2）網(wǎng)格的細(xì)分通過網(wǎng)格的細(xì)分，使每個(gè)單元的面積縮小，那么盡管每個(gè)單元是應(yīng)變、常應(yīng)力單元，仍可較好地反映結(jié)構(gòu)中的應(yīng)力變化，使得到的解答收斂于問題的精確解。（3）網(wǎng)格合理布局 根據(jù)應(yīng)力梯度使網(wǎng)格的布局合理化。即在梯度大的區(qū)域網(wǎng)格密些，梯度小的區(qū)域應(yīng)稀些。密、稀網(wǎng)格之間應(yīng)逐步過渡。帶圓孔方板的網(wǎng)格劃分（4）改用高階單元 受集中力的懸臂梁，采用128個(gè)三結(jié)點(diǎn)三角形常應(yīng)變單元，以及3個(gè)八結(jié)點(diǎn)四邊形高階單元結(jié)果?？梢?，采用高階元的計(jì)算精度比常應(yīng)變?cè)叩枚唷?.3 有限分析軟件1、三角

9、形單元分析程序設(shè)計(jì) 教材P.26程序FEM1。2、ANSYS軟件 結(jié)構(gòu)、熱、流體、電磁學(xué)、聲學(xué)等通用分析軟件。3、SAP2000軟件 土木結(jié)構(gòu)分析軟件。2.4 算例分析算例1：對(duì)角受壓的正方形薄板，載荷沿厚度均勻分布，為2N/m。由于對(duì)稱性，取1/4部分作為計(jì)算對(duì)象，試用有限元程序進(jìn)行計(jì)算。 2N/m2N/m2m2mxy算例2：簡支梁，梁高3m,跨度18m，厚度1m，承受均布荷載10N/m2。已知 按平面應(yīng)力問題進(jìn)行計(jì)算。18m3mxy網(wǎng)格劃分考察點(diǎn)y(m)-1.25-0.75-0.250.250.751.25有限元結(jié)果19711436-36-114-197彈性力學(xué)結(jié)果22513444-44-134-225誤差28208-8-20-28考察點(diǎn)y(m)-1.25-0.75-0.250.250.751.25有限元結(jié)果16.231.237.233.720.73.6彈性力學(xué)結(jié)