西交大有限元原理及的應(yīng)用-大作業(yè)

1、有限元原理及工程應(yīng)用-大作業(yè)學(xué) 院：機(jī)械工程學(xué)院班 級(jí)：碩4002班小組成員：李 追3114001089陳 草31140010802015.5.19作業(yè)題目：禾I用有限元方法對(duì)簡(jiǎn)支梁?jiǎn)栴}進(jìn)行求解，梁的橫截面為矩形，其約束情況如圖已知梁的幾何尺寸和物理參數(shù)如下：(1)幾何尺寸：長(zhǎng)度L 40cm，截面尺寸b t 2cm 0.2cm；圖 1 .梁及其橫截面示意圖要求:(1)至少劃分五個(gè)節(jié)點(diǎn)(四個(gè)單元)；給出單元節(jié)點(diǎn)信息；給出單元?jiǎng)偠染仃嚭唾|(zhì)量矩陣；給出總剛度矩陣和總質(zhì)量矩陣；求出梁各界固有頻率及振型(五階)1所示。篥的橫戟面(2)物理參數(shù)：彈性模量E 70GPa，泊松比30.3，密度=2700kg/

2、m。(v)(h)(6)將所得結(jié)果與理論值進(jìn)行對(duì)比，驗(yàn)證方法的可行性。解：由有限元知識(shí)，根據(jù)Rayleigh-Ritz法，解有限元分為四步：建立離散化、單元分析、 成總體方程、解方程，具體步驟如下：（1）建立離散化這里我們將矩形截面簡(jiǎn)支梁等分四等分，即分為六節(jié)點(diǎn)的五個(gè)桿單元，如圖2所示:L 40每個(gè)單元尺寸丨-40cm=8cm,這里只考慮桿在豎直平面的彎曲，每個(gè)節(jié)點(diǎn)只有y方向55位移和繞z軸的旋轉(zhuǎn)自由度。(2)單元分析構(gòu)造一組Lagrange插值基函數(shù)，在本節(jié)點(diǎn)值為1，其他節(jié)點(diǎn)值為0。從Rayleigh-Ritz法可以看到，插值函數(shù)要p次可微，最高階導(dǎo)數(shù)出現(xiàn)在應(yīng)變能表達(dá)式中；同樣，我們使用有限元

3、逼近方法。梁的彎曲問(wèn)題，應(yīng)變能計(jì)算公式:U -LEI202 +2v x,t2dx可以這一原則適用于基函數(shù)的選擇以及形狀函數(shù),否則我們將無(wú)法正確計(jì)算應(yīng)變能當(dāng)我們(1-1)其中，E為彈性模量，Iz為截面慣性矩。從公式可知，位移函數(shù)必須連續(xù)，并且二階 導(dǎo)數(shù)平方可積。如圖3，是一維桿單元模型，每個(gè)節(jié)點(diǎn)兩個(gè)自由度，該單元含有四個(gè)自由度，即(Vi,zi,Vj,zj)。本題中我們采用三次多項(xiàng)式插值函數(shù)：u x12x3x24x3(1-2)-V因此，我們必須給出四個(gè)形函數(shù)(位移模式)。圖3一維桿單元模型1)構(gòu)造 Hermite 插值函數(shù)。選擇局部坐標(biāo)系x l (ni0,0 ,njl,0 )，其中l(wèi)是單元長(zhǎng)度，轉(zhuǎn)

4、角z是撓度值v的一階導(dǎo)數(shù)，定義邊界條件：V xivix(1-3)因此，我們給出變形的Hermite的多項(xiàng)式插值函數(shù):2VH(0)i 1節(jié)點(diǎn)位移值也可以得出同時(shí)，表達(dá)式（1-7）用矩陣表示為3）用能量表達(dá)式替代表達(dá)式中的動(dòng)能表達(dá)式:其中，H（0）和H分別滿足如下條件,對(duì)應(yīng)的圖形如圖4所示:dH(0)djijjij=12)NiN2Hi(0)3223(1-5)N3= Hi(0)=H20)=3N4二Hi=HTl(1-6)ViN2ziN3VjN4zj(1-7)Ve(1-8)T其中，VeVj,zi,Vj,zjNiN2N3N42ViHiVi 1(1-4)&eTMe&e(1-9)T IOAV將

5、（1-8）帶入（1-9），得到A0(1-10)從而獲得質(zhì)量矩陣：Me(1-11)帶入X l,15622l5413lAl22l4l213l3l2Me4205413l15622l2 213131222l 4l2(1-12)應(yīng)變能表達(dá)式:1Elo2v剛度矩陣表達(dá)式:帶入x l，可以得到（3） 形成總體方程(1-13)1 -Ve2KeElzl3Vel3(1-14)12 6l 12 6l2 2EIz6l 4l 6l 2l3l 12 6l 12 6l2 26l 2l 6l 4l(1-15)將每個(gè)桿單元的能量方程組裝。完整梁上的的總動(dòng)能和能量的和所做的總功梁上的外力作用，所有的自由度的位移矢量可以給出:&a

6、mp;eTMe&e(1-9)z6(1-16)V1z1V2z2V3z3V4z4V5z5V6將位移矢量轉(zhuǎn)換為全局坐標(biāo)系下的位移矢量，變換矩陣為:veaev（1-17）100000000000a1010000000000001000000000000100000000001000000000a20001000000000000100000000000010000000 000100000000 00001000000a3（1-18）0 000001000000 000000100000 0000010 00000 0000001 0000a40000000010000 0000000 010

7、0000000001000a50000000001000000000000100000000000 01分別以矩陣形式給出動(dòng)能和質(zhì)量矩陣:5TMaeMeaeTtotTaeae(1-19)(1-20)因此，總體質(zhì)量矩陣為總體剛度矩陣:126112610000000061412612120000000012612401261000000612120護(hù)6121200000000126124012610000KElz0061212081261212000013000012612401261000000612120812612120000000012612401261000000612120812612

8、120000000012611261000000006121261412求固有頻率。15622154131022141213131205413131205413131208121310054131312MAl0013131204200000540000131000000000000000000000000000000000013100000031200000005413100008121313120000131312054131003120812131312000541313120541310131312081213131200054131156221000131312221412總應(yīng)變能(1

9、-21)(4)解方程，有限元分析結(jié)果基于上述理論，我們獲得了采用MATLAB的有限元分析程序代碼。提交邊界條件、材料特性和幾何參數(shù)到上面的方程，使用MATLAB代碼我們得到以下結(jié) 果。(1)離散簡(jiǎn)支梁為5單元6節(jié)點(diǎn),那么我們得到的單元質(zhì)量矩陣和剛度矩陣，如下所示：2.18750.0875 -2.18750.087540.0875 0.0047 -0.08750.0023Ke 1.0 10-2.1875 -0.08752.1875 -0.08750.08750.0023 -0.08750.0047(2)總質(zhì)量矩陣和總剛度矩陣如下:0.00320.0000 0.0011-0.00000.00000

10、.0000 0.0000-0.00000.00110.00000.0032-0.0000-0.0000-0.0000 -0.00000.0000Me1Utot2aeTKeaeaeKeae(1-22)利用Lagrange方程，推導(dǎo)簡(jiǎn)支梁自由振動(dòng)方程:U&(1-23)這里,我們假設(shè)簡(jiǎn)支梁做簡(jiǎn)諧振動(dòng)，則(1-24)因此, 特征方程為:(1-25)其中,2， 為固有頻率。0.00320.00000.0011-0.000000 0000000.00000.00000.0000-0.000000 0000000.00110.00000.006400.0011-0.0000 000000-0.000

11、0 -0.0000 00.00000.000C)-0.0000 000000000.00110.0000 0.00640 0.0011-0.0000000000 -0.0000 -0.000000.0000 0.0000-0.00000000M0000 0.00110.00000.006400.0011-0.0000000000 -0.0000 -0.0000 00.00000.0000-0.000000000000 0.00110.00000.006400.0011-0.0000000000 -0.0000-0.000000.00000.0000-0.0000000000 000.00110

12、.00000.0032-0.0000000000 00-0.0000-0.0000-0.00000.00002.18750.0875-2.18750.0875000000000.0875 0.0047 -0.08750.002300000000-2.1875 -0.0875 4.37500-2.18750.08750000000.0875 0.002300.0093 -0.08750.002300000000-2.1875-0.08754.37500-2.18750.08750000000.08750.002300.0093-0.08750.00230000K 1.0 1040000-2.18

13、75-0.08754.37500-2.18750.08750000000.08750.002300.0093-0.08750.002300000000-2.1875 -0.08754.37500-2.18750.08750000000.0875 0.002300.0093-0.08750.002300000000-2.1875-0.08752.1875-0.0875000000000.08750.0023-0.08750.0047（3）簡(jiǎn)支梁的振動(dòng)分析表1列出了簡(jiǎn)支梁的五階固有振動(dòng)頻率。從表中我們可以看出，有限元模擬分析方法和 理論值在誤差允許范圍內(nèi)是比較吻合的。計(jì)算得離散為五單元下的簡(jiǎn)支梁固

14、有振動(dòng)頻率:計(jì)算值（104）理論值（104）誤差率（%）一階固有振動(dòng)頻率0.01810.01810.0107二階固有振動(dòng)頻率0.07250.07270.1657三階固有振動(dòng)頻率0.16320.16450.7942四階固有振動(dòng)頻率0.29010.29682.3037五階固有振動(dòng)頻率0.45330.503210.9918離散為30單元的簡(jiǎn)支梁的振動(dòng)模式如下圖所示五階振動(dòng)圖像六階振動(dòng)圖像3.采用Matlab編寫(xiě)的程序代碼%利用有限元方法求解簡(jiǎn)支梁的振動(dòng)問(wèn)題%clc ; clear all;syms x l rho b t EA = b*t;I = b*tA3/12;n = input (Please

15、 input the number of discrete elements n = ); %量n%定義形函數(shù)N1 = 1-3*(x/l)A2+2*(x/l)A3;N2 = (x/l-2*(x/l)A2+(x/l)A3)*l;N3 = 3*(x/l)A2-2*(x/l)A3;N4 = (x/l)A3-(x/l)A2)*l;%求解單元質(zhì)量矩陣、剛度矩陣以及總質(zhì)量矩陣和剛度矩陣N = N1,N2,N3,N4;Me0 =int(N*N,x,0,l);Me = rho*A*Me0Ke0 = int(diff(N.,x,2)*diff(N,x,2),x,0,l)Ke = (E*I)*Ke0;輸入離散化單

16、元的數(shù)M0 = zeros(2*(n+1),2*(n+1);K0 = zeros(2*(n+1),2*(n+1);for i=1:1:nae = zeros(4,2*(n+1);for j=1:1:4ae(j,2*i+j-2) = 1; %定義坐標(biāo)變換矩陣endM0=M0+ae.* Me * ae;K0 =K0+ae.* Ke *ae;enddisp(The element mass matrix Me = ); disp(Me);disp(The element stiffness matrix Ke =);disp(Ke);disp(The total mass matrix M = );

17、disp(M0);disp(The total stiffness matrix K = );disp(K0);%Me1 = matlabFunction(Me);Ke1 = matlabFunction(Ke);M1 = matlabFunction(M0); %將質(zhì)量符號(hào)矩陣轉(zhuǎn)換為代數(shù)矩陣K1 = matlabFunction(K0); %將剛度符號(hào)矩陣轉(zhuǎn)換為代數(shù)矩陣%輸入的幾何參數(shù)和物理常數(shù)L = 0.4; %梁的長(zhǎng)度, mb = 0.02; %梁的寬度, mt = 0.002; %梁的厚度, mE = 0.7*10X1; %梁的彈性模量,GParho = 2700; %梁的密度,kg/

18、mA3l = L/n;MeNumeric = Me1(b,l,t,rho);KeNumeric = Ke1(E,b,l,t);MNumeric =M1(b,l,t,rho);KNumeric =K1(E,b,l,t);disp(The numerical element mass matrix Me= );disp(MeNumeric); %陣disp(The numerical element stiffness matrix Ke = );disp(KeNumeric); %矩陣disp(The numerical total mass matrix M = );disp(MNumeric

19、); %disp(The numerical total stiffness matrix K = );disp(KNumeric); %簡(jiǎn)支梁的振動(dòng)求解%對(duì)于簡(jiǎn)支梁考慮約束條件：v 00,v L 0;%采用消元法：消去整體質(zhì)量矩陣的第1行，第1列；倒數(shù)第二行和倒數(shù)第二列輸出單元質(zhì)量矩輸出單元?jiǎng)偠容敵隹傎|(zhì)量矩陣輸出總剛度矩陣%消去整體剛度矩陣的第1行，第1列；倒數(shù)第二行和倒數(shù)第二列Mssb = MNumeric(2:2*n,2*(n+1),2:2*n,2*(n+1);Kssb = KNumeric(2:2*n,2*(n+1),2:2*n,2*(n+1);Xssb,lamdassb=eig(Kssb,Mssb);omegassb=sort(sqrt(diag(lamdassb);XXssb=zeros(n,n);XXssb(1,:)=0;XXssb(n,:)=0;for i=2:n-1for j=1:nXXssb(i,j)=Xssb(2*(i-1),j);endend%根據(jù)振動(dòng)力學(xué)方程得到簡(jiǎn)支梁固有頻率betalssb = ;omega_realssb = ;errorssb = ;for i=1:10betalssb(i) = pi*i;omega_realssb(i) = (be