1、精選優(yōu)質文檔-傾情為你奉上有限元分析半開卷資料基本概念：位移a=u1v1u2v2unvnT；等效結點荷載R=R1xR1yR2xR2yRnxRnyTu=uv=Nae; Fe=kae; =xyxyT=Lu=LNae=Bae; =xyxyT=D=DBae=Sae第二章 平面彈性力學問題1.位移模式與收斂性條件三角形單元的位移模式：u=Niui+Njuj+Nmumv=Nivi+Njvj+Nmvm式中Ni=ai+bix+ciy2A i,j,m；ai=xjym-xmyj;bi=yj-ym; ci=-(xj-xm), A為單元面積A=121xiyi1xjyj1xmym (為了使面積A不成為負值，規(guī)定結點I,

2、 j, m的次序按逆時針轉向)矩形單元的位移模式：u=Niui+Njuj+Nmum+Npupv=Nivi+Njvj+Nmvm+Npvp式中Ni=141+ixa1+iyb; i=xixi; i=yiyi (i,j,m,p)位移模式需滿足的條件：a 位移模式必須能反映單元的剛體位移（與本單元的變形無關的位移） b. 位移模式必須能反映單元的常量應變（與位置坐標無關的應變） c. 位移模式應當盡可能反映位移的連續(xù)性（相鄰單元之間位移的連續(xù)性和單元內部位移也是連續(xù)的）以三角形單元為例說明，a1,a4,a5-a32反映了剛體移動和剛體轉動，a2,a6,a3+a5反映常量應變2.形函數(shù)及其性質形函數(shù)：即插

3、值基函數(shù)，反應單元的位移形態(tài)，因而也稱為位移的形態(tài)函數(shù)，簡稱形函數(shù)(NiNjNm)形函數(shù)的性質：（2） 在節(jié)點i上Ni=1，在其他節(jié)點上Ni=0，該性質可導出形函數(shù)在三角形單元上的積分和在某邊界上的積分為eNidxdy=13A,ijNids=12lije表示對單元積分，ij表示對單元的ij邊線積分（3） 在單元中，任意點形函數(shù)之和等于1（4） 形函數(shù)的值在0-1間變化3.面積坐標及三角形高次形函數(shù)的構造面積坐標：三角形單元中，任一點P(x,y)與其3個角點相連形成3個子三角形，其位置可以用三個比值來確定，即面積坐標。面積坐標只限于用在一個三角形單元內，因而是一種局部坐標。在平行于jm邊的一根直

4、線上的所有點，都具有相同的Li坐標，而且這個坐標就等于“該直線至jm邊的距離”與“結點i至jm邊的距離”的比值用直角坐標表示面積坐標的關系式：Li=(ai+bix+ciy)2A式中ai=xjym-xmyj,bi=yj-ym,ci=-xj+xm用面積坐標表示直角坐標的關系式：x=xiLi+xjLj+xmLmy=yiLi+yjLj+ymLm六結點三角形單元的位移模式：u=Niui+Njuj+Nmum+N1u1+N2u2+N3u3v=Nivi+Njvj+Nmvm+N1v1+N2v2+N3v3Ni=Li2Li-1 (i,j,m),N1=4LjLm (1,2,3)4.有限元支配方程的推導（結構力學法/變

5、分原理）根據(jù)平衡條件，各環(huán)繞單元對該結點作用的結點力之和應等于由各環(huán)繞單元移置而來的結點荷載之和，即eFi=eRi，再將結點力公式代入，變可得有限元的支配方程Ka=R,式中K為整體剛度矩陣，a為整體結點位移列陣，R為整體結點荷載列陣變分原理導出有限元支配方程證明：用最小勢能原理推導出有限元的求解方程=12TDtdxdy-uTftdxdy-SuTfds式中，t是平面彈性體的厚度，f是體積力，f是物體表面的面力離散成有限網(wǎng)格，其中=Bae, u=Nae代入上式可得=12e(ae)TeBTDBtdxdyae-e(ae)TeNTftdxdy-e(ae)TseNTftds利用ae=Cea, 則=12aT

6、Ka-aTR, 其中K=eCeTkCe;R=eCeTRe;k=eBTDBtdxdy; Re=eNTftdxdy+seNTftds根據(jù)最小勢能原理=0，即a=0, 這樣就得到有限元的求解方程Ka=R虛功原理建立有限元的支配方程假設單元發(fā)生了虛位移，其相應虛位移為u=uvT而該單元上各結點的相應虛位移為ae=uiviujvjumvmT按照靜力等效原理，即結點荷載與原荷載在上述虛位移上的虛功相等，有(ae)TRe=uTP將u=Nae 代入，得(ae)TRe=（ae）TNTP由于虛位移是任意的，得Re=NTP5.荷載列陣：單元到整體體力引起的等效結點荷載：Re=-bb-aaNTftdxdy面力引起的等

7、效結點荷載：Re=-bbNTftdy常見分布荷載產生的等效結點荷載1.單元自重 Re=-13gtA010 101T 2.在ij邊界上受x方向均布力q作用，邊界長度為l，Re=12qlt101 00 0T 3.在ij邊界上受三角形分布荷載作用，邊界長度為l，Re=12ql23013 00 0T 整體結點荷載列陣：確定每個單元的結點荷載列陣，然后根據(jù)各個單元的結點局部編碼與整體編碼的對應關系，將單元的結點荷載列陣中每個子矩陣疊加到R的相應位置上。6.剛度矩陣：單元到整體單元剛度矩陣（平面應變問題）k=-bb-aaBTDBtdxdy線性位移模式下，有k=BTDBt=kiikijkimkjikjjkj

8、mkmikmjkmmFe=UiViUjVjUmVm=kiixxkiixykiiyxkiiyyuiviujvjumvm kij表示j結點對i結點的剛度貢獻kiixx+kjixx+kmixx=0;kiiyx+kjiyx+kmiyx=0單元剛度矩陣的力學意義：單元剛度矩陣中任一個元素（如kijyx）表示當j結點x方向發(fā)生單位位移時，在i結點y方向產生的結點力單元剛度矩陣的性質：對稱性，奇異性，主元素恒正，單元均勻放大或縮小不會改變剛度矩陣的數(shù)值，單元水平或豎向移動不會改變剛度矩陣的數(shù)值平面應力問題 krs=Et41-v2Abrbs+1-v2CrCsvbrbs+1-v2CrCsvCrCs+1-v2br

9、bsCrCs+1+v2brbs (r=I,j,m s=I,j,m)平面應變問題krs=E1-vt41+v1-2vAbrbs+1-2v2(1-v)crcsv1-vbrcs+1-2v2(1-v)crbsv1-vcrbs+1-2v2(1-v)brcscrcs+1-2v2(1-v)brbs整體單元建立步驟：將K全部充零，逐個單元地建立單元 的剛度矩陣，然后根據(jù)單元結點的局部編碼與整體編碼的關系，將單元的剛度矩陣中每一個子矩陣疊加到K中的相應位置上。對所有的單元全部完成上述疊加步驟，就形成了整體剛度矩陣。7.簡單問題的有限元具體計算（利用已算好的剛度矩陣）8.計算結果的整理與分析計算成果主要包括兩個方面

10、：即位方面和應力方面。主要進行對應力方面的計算成果討論。把計算出來的常量應力作為單元形心處的應力，拉應力用箭頭表示，壓應力用平頭表示。由計算成果推出結構物內某一點的接近實際的應力，必須通過某種平均計算，通常可采用繞結點平均法或二單元平均法。繞結點平均法：把環(huán)繞某一結點的各單元中的常量應力加以平均，用來表征該結點處的應力。（在內結點處具有較好的表征性，但在邊界結點處表征性較差，因此邊界結點處的應力宜用內結點應力插值推算）二單元平均法：把兩個相鄰單元中的常量應力加以平均用來表征公共邊中點處的應力。PS：只容許對厚度及彈性常數(shù)都相同的單元進行平均計算9.網(wǎng)格劃分的注意事項a. 應根據(jù)所需精度，在合理

11、的計算時間內來決定單元的大小 b.為了便于整理和分析應力成果，往往采用直角三角形的單元 c.應合理利用結構的對稱性，減少計算量 d.為了使地基彈性對結構物中應力的影響能反映出來，必須把和結構物相連的那一部分地基也取為彈性體，和結構物一起作為計算對象e. 可進行二次計算在以下情況下應將單元尺寸減?。篴.計算對象的厚度或其彈性有突變之處，且應吧突變線作為單元的界線，b.計算對象受季度突變的分布荷載或受集中荷載 c.結構物具有凹槽或孔洞第三章 平面等參有限元1. 等參單元的概念將位移模式與坐標變換式具有相同的形式，即形函數(shù)相同，參數(shù)個數(shù)相同的單元稱為等參單元。2.單元形函數(shù)的構造方法整體坐標系為（x

12、,y），在每個單元上建立局部坐標系(,),只需雅克比行列式大于0即可保證坐標一一對應的實現(xiàn)。J=xyxy>0坐標轉換式：x=N1x1+N2x2+N3x3+N4x4 y=N1y+N2y2+N3y3+N4y4Ni=14(1+i)(1+i)3.位移模式及收斂性要求單元位移模式：u=N1u1+N2u2+N3u3+N4u4 v=N1v1+N2v2+N3v3+N4v4收斂性要求：完備性和連續(xù)性以4結點四邊形等參單元為例說明完備性，i=14Ni=1，i=14Nixi=x，i=14Niyi=yi=14Ni=14(1-)1-+141+1-+141+1+141-1+=121-+121+=1對于具有m各結點的

13、等參單元，為了位移模式滿足完備性要求，形函數(shù)必須滿足i=1mNi=1連續(xù)性：只需考察任意兩單元交界面上位移是否連續(xù)即可。為了保證整體坐標與局部坐標的一一對應關系，單元不能歪斜，單元的各邊長不能等于0。4.單元剛度矩陣和荷載列陣的計算Bi=Nix00NiyNiyNix; Si=E1-v2NixvNiyvNixNiy1-v2Niy1-v2Nix 需根據(jù)復合函數(shù)的求導法則來求出各形函數(shù)對整體坐標的導數(shù)Ni=Nixx+Niyy; Ni=Nixx+Niyy 即可解得Nix和Niy單元剛度矩陣：k=eBTDBtdA=-11-11BTDBtJdd體力引起的單元結點荷載Re=eNTftdA=-11-11NTf

14、tJdd面力引起的單元結點荷載，在=±1；Re=-11Ntft(x)2+(y)2d;在=±1；Re=-11Ntft(x)2+(y)2d以平行四邊形為例，常見荷載下的單元結點荷載a. 在自重f=0-g;Re=-14gtATb. 在=1;Re=12atqTc. 在=1上法向分布面力的集度為q, Re=-11NTy-xqtd5.高斯數(shù)值積分數(shù)值積分一般有兩類方法，一類是等間距數(shù)值積分，如辛普森方法等。另一類是不等間距數(shù)值積分，如高斯數(shù)值積分法。高斯積分法對積分點的位置進行了優(yōu)化處理，所以精度較高。例：求兩點高斯積分的積分點坐標和積分權系數(shù)二次多項式：P=(-1)(-2)積分點位置

15、由下式確定-11iPd=0 i=0,1； P=-1-2*(-n)當i=0時，-11-1-2d=23+212=0當i=1時，-11-1-2d=-231+2=0聯(lián)立方程后可解得1=-0.577；2=0.577積分權系數(shù)為：Hi=-11lin-1()d; lin-1=-1-2*-i-1-i+1*(-n)i-1i-2*i-i-1i-i+1*(i-n)可解得H1=H2=1第四章 空間彈性力學問題1.單元形函數(shù)的構造方法V為四面體的體積V=161xiyizi1xjyjzj11xmxpymypzmzp 式中 ai=xjyjzjxmymzmxpypzp; bi=-1yjzj1ymzm1ypzp; ci=-xj

16、1zjxm1zmxp1zp; ai=-xjyj1xmym1xpyp12.位移模式及收斂性要求位移模式：u=Niui+Njuj+Nmum+Npup ; v=Nivi+Njvj+Nmvm+Npvpw=Niwi+Njwj+Nmwm+Npwp; Ni=16Vai+bix+ciy+diz (i,j,m,p)3.單元剛度矩陣和荷載列陣的計算krs=E(1-v)361+v1-2vV*brbs+A2(crcs+drds)A1brcs+A2crbsA1brbs+A2drbsA1crbs+A2brcscrcs+A2(brbs+drds)A1crds+A2drcsA1drbs+A2brdsA1drcs+A2crds

17、drds+A2(crcs+brbs)若單元受自重作用，體積力f=00-gT;Re=-14gV1T,表面單元重量平均分到4個結點。若單元某邊界面如ijm面x方向受線性分力作用，設結點i的面力集度為q，結點就，p集度為0，則面力矢量f=Niq00T，Re=13qAijm0000T4.高斯數(shù)值積分第六章 有限元程序設計1.主元素序號指示矩陣MA由子程序CBAND來實現(xiàn)（1） 將MA（*）全部充零（2） 對全部單元循環(huán)，利用結構自由度序號矩陣JR（2，*），計算每個單元的結點自由度序號數(shù)組NN（8）單元i，NNx=單元ii，NNx=單元iii，NNx=單元iv，NNx=（3） 計算單元自由度序號中最小

18、的序號，賦給L（4） 計算與單元各自由度對應行的半帶寬。把各單元計算得到的半帶寬中的最大值賦給MA（*）數(shù)組, MA*=9（5） 對數(shù)組MA（*）逐個元素累加，便得到所需的主元素序號指示相量，此時MA*=2.半帶寬的計算K中每行從第一個非零元素到對角線元素的元素個數(shù)稱為該行的半帶寬。半帶寬是衡量結點編號好壞的一個指標，最優(yōu)的結點編號可以使半帶寬最小半帶寬等于該行行號（即自由度序號）減去與該自由度號關聯(lián)的最小自由度號再加1，即JPL=JP-L+1,JPL為半帶寬，JP為自由度序號，L為某一單元最小自由度序號。3.一維變帶寬的存儲將按一維變帶寬存儲的K所需的存儲量記為NH，存儲有效元素的一維數(shù)組記

19、為SK（NH），結構總自由度數(shù)，亦即需要建立的方程個數(shù)記為N。則對于本例來講，N=12，NH=66。通過指示向量MA（*）便可將一維數(shù)組SK（*）與二維矩陣K的元素一一對應起來一維數(shù)組與二維矩陣K有如下對應關系（1） K中第I行的主元素就是SKMA(I)（2） KI,J=SKMAI-(I-J)4.結點自由度序號矩陣JRJR2,12=5.體力引起的荷載列陣6.面力引起的荷載列陣7.荷載列陣的集成8.剛度矩陣的集成由子程序SKO來實現(xiàn)（1） 采用3個積分點的高斯數(shù)值積分來計算單元剛度矩陣，對高斯點局部坐標和權系數(shù)賦值（2） 在對單元循環(huán)的循環(huán)體內形成單元自由度數(shù)組NN（8）和單元結點坐標值XYZ（

20、2,4）。NN（8）是為了將單元剛度矩陣的元素累加到整體剛度矩陣的需要，XYZ（2,4）是為了計算單元剛度矩陣的需要。（3） 調用子程序STIF計算單元的剛度矩陣SKE（8,8）（4） 將每個單元的剛度矩陣累加到一維剛度矩陣數(shù)組SK（*）中去。9.讀懂源程序10在源程序基礎上修改某些功能三角形單元有關公式1. 位移模式（三角形單元）：u=Niui+Njuj+Nmumv=Nivi+Njvj+Nmvm式中Ni=ai+bix+ciy2A i,j,m；ai=xjym-xmyj;bi=yj-ym; ci=-(xj-xm), A為單元面積A=121xiyi1xjyj1xmym (為了使面積A不成為負值，規(guī)

21、定結點I, j, m的次序按逆時針轉向)2. 應變轉換矩陣B=12Abi0bj0bm00cicibi0cjcjbj0cmcmbm 3. 應力轉換矩陣平面應變問題Si=E2-v2Abivcivbici1-2v2(1-v)ci1-2v2(1-v)bi 對于平面應力問題，要把E換成E1-v2, v換成v1-v;Si=E(1-v)21+v1-2vAbiv1-vciv1-vbici1-2v2(1-v)ci1-2v2(1-v)bi4. 單元剛度矩陣平面應力問題 krs=Et41-v2Abrbs+1-v2CrCsvbrbs+1-v2CrCsvCrCs+1-v2brbsCrCs+1+v2brbs (r=I,j

22、,m s=I,j,m)平面應變問題krs=E1-vt41+v1-2vAbrbs+1-2v2(1-v)crcsv1-vbrcs+1-2v2(1-v)crbsv1-vcrbs+1-2v2(1-v)brcscrcs+1-2v2(1-v)brbs5. 等效結點荷載體力引起的等效結點荷載：Re=-bb-aaNTftdxdy面力引起的等效結點荷載：Re=-bbNTftdy常見分布荷載產生的等效結點荷載1.單元自重 Re=-13gtA010 101T 2.在ij邊界上受x方向均布力q作用，邊界長度為l，Re=12qlt101 00 0T 3.在ij邊界上受三角形分布荷載作用，邊界長度為l，Re=12ql23

23、013 00 0T 4結點四邊形等參單元的有關主要公式1. 位移模式u=uv=Nae；其中N=N10N20N30N400N10N20N30N4;ae=u1v1u2v2u3v3u4v4T;Ni=14(1+i)(1+i) 2. 坐標變換x=N1x1+N2x2+N3x3+N4x4 y=N1y+N2y2+N3y3+N4y43. 應變公式=xyxyT=Bae 其中B=B1B2B3B4; Bi=Nix00NiyNiyNix; NixNiy=J-1NiNi4. 應力公式=D; x=D1x+D2y;y=D2x+D1y; xy=D3xy 其中D=D1D20D2D1000D3; D1=1-vE(1+v)(1-2v

24、); D2=vE(1+v)(1-2v); D3=E2(1+v)5. 單元剛度矩陣k=-11-11BTDBJdd; k=k11k12k13k14k21k22k23k24k31k41k32k42k33k43k34k44kij=-11-11Nix0Niy0NiyNixD1D20D2D1000D3Nix00NiyNiyNixJdd 令每一子塊矩陣為kij=kiikijkjikjj 其中kii=-11-11(D1NixNjx+D3NiyNjy)Jdd ;kij=-11-11(D2NixNjy+D3NiyNjx)Jddkji=-11-11(D2NiyNjx+D3NixNjy)Jdd; kjj=-11-11

25、(D1NiyNjy+D3NixNjx)Jdd6. 單元等效結點荷載（1） 集中力P=PxPyT; Re=NTP; Rix=NiPx;Riy=NiPy （2） 自重體力f=0-gT; Re=-11-11NT0-gJdd ; Rix=0; Riy=-g-11-11NiJdd （3） 分布面力f=fxfyT 在=±1邊界，Re=-11Nf(x)2+(y)2d在=±1邊界，Re=-11Nf(x)2+(y)2d如果面力為法向壓力，則f=-qlqm(i) 在=-1邊界 l=-y(x)2+(y)2, m=x(x)2+(y)2 ; Re=-11NTy-xqd(ii) 在=1邊界 l=y(x)2+(y)2, m=-x(x)2+(y)2; Re=-11NTy-xqd(iii) 在=-1邊界l=y(x)2+(y)2, m=-x(x)2+(y)2；Re=-11NTy-xqd子程序的主要功能專心-專注-專業(yè)INPUT輸入原始數(shù)據(jù)CBAND形成主元素序號指示矩陣MA（600）SKO形成整體剛度矩陣KCONCR計算集中力引起的等效結點荷載ReBODYR計算自重體力引起的等效結點荷載ReFACER計算分布面力引起的等效結點荷載ReDECOP支配方程Crout直接解法中的分解和前代計算FOBAC